分解子域上的随机游走

# 在分解子域上行走 来源：https://clementjambon.github.io/wods/index.html 博客文章作者：Clément Jambon (https://clementjambon.github.io/) ## 1 椭圆型偏微分方程与边值问题 许多现实世界中的现象都由椭圆型偏微分方程（PDE）所支配：热传导、静电学、路径规划、稳态势流等等。它们通常被表述为边值问题（BVP），即在区域边界上给定值，然后我们求解区域内部的 PDE 解。 以[拉普拉斯方程](https://en.wikipedia.org/wiki/Laplace%27s_equation)及 Dirichlet 边界条件为例： $$ \begin{cases} \Delta u = 0 & \text{in } \Omega \\ u = g & \text{on } \partial\Omega_D \end{cases} $$ 你可以随意操作下面的交互式图形，以直观感受它的作用： 笔刷值+1.00 −1（冷）+1（热） 预设 点击并在外部边界带上绘制以涂画 Dirichlet 值 $g$。内部满足 $\Delta u=0$。 **Dirichlet 问题。**在方形区域上求解带 Dirichlet 边界条件的拉普拉斯方程。 我们可以考虑更复杂的几何形状和边界条件，让问题更有趣。例如，人们常常对混合边值问题感兴趣，其中包含 Neumann 边界条件请注意，我们这里仅限制在零 Neumann 条件。： $$ \begin{cases} \Delta u = 0 & \text{in } \Omega \\ u = g & \text{on } \partial\Omega_D \\ \frac{\partial u}{\partial n} = 0 & \text{on } \partial\Omega_N \end{cases} $$ 下面的交互式图形对此进行了说明。请注意等值线是如何弯曲，以便与零 Neumann 障碍物成直角相交，从而满足 $\frac{\partial u}{\partial n} = 0$。 笔刷值+1.00 −1（冷）+1（热） 场景 边界预设 像之前一样点击并绘制 Dirichlet 带。内部障碍物（虚线、灰色）强制执行 $\partial u/\partial n = 0$ —— 等值线弯曲并以直角与其相交。 **混合问题。**拉普拉斯方程，在外边界上涂画 Dirichlet 值，内部为零 Neumann 几何形状。 如果仔细观察，你会发现解实际上是“像素化的”。这是因为它是在网格上用[有限差分](https://en.wikipedia.org/wiki/Finite_difference_method)计算的。有限差分非常容易理解和实现，但它不太擅长处理复杂几何形状，通常需要极度精细的网格划分。另一种常见的选择是使用[有限元](https://en.wikipedia.org/wiki/Finite_element_method)。问题在于有限元需要仔细的网格生成，对于复杂几何形状尤其困难且耗时试想设计一辆汽车时，每次调整设计都必须重新生成网格。那将是一场噩梦——事实也的确如此！，例如下图所示的城市。 **城市周围的风。**这个场景包含数百座具有复杂几何形状的建筑物。我们的方法描述了它们对风模式的影响，这里可视化为稳态势流线。 ## 2 无网格蒙特卡洛方法 幸运的是，计算机图形学界最近复兴了一个古老的想法：无网格蒙特卡洛方法。支撑此方法的核心算法是*Walk on Spheres*（WoS）算法。 其直觉来自随机微积分：如果你模拟一个从内部点 $x$ 出发的[布朗运动](https://en.wikipedia.org/wiki/Brownian_motion)（一种连续随机游走），它最终会在某个随机位置 $Z_\tau$ 击中边界。在该随机位置处边界条件的期望值 $\mathbb{E}[g(Z_\tau)]$ 恰好就是给出的 Dirichlet 问题的解。 **布朗运动与拉普拉斯方程。**一个粒子从内部点 $x$ 出发，扩散直到它被吸收在某个随机边界位置 $Z_\tau$。边界按 $g$ 着色；对许多次这样的游走求 $g(Z_\tau)$ 的平均值即可恢复 $u(x)$。 在每个内部点独立地执行此操作，我们就能得到整个解。每个点只需少量游走时估计是有噪声的，但随着我们平均越来越多的游走，方差逐渐消失，平滑的调和解就会出现，如下图所示。 刷新率：1每个/帧 慢快 低 $u$高 $u$ 每个内部像素运行独立的 Walk-on-Spheres 估计 $u(x)=\mathbb{E}[g(Z_\tau)]$。对每个像素的更多游走进行平均，会将场降噪为调和解。它会持续运行；滑块设置运行速度。 样本数：1 **渐进蒙特卡洛估计。**同样的方形 Dirichlet 问题，现在通过蒙特卡洛在每个像素处求解。每个像素平均许多随机游走：样本少时内部有噪声，随着样本数量增加噪声逐渐消失。 然而，模拟布朗运动计算量很大，而且通常有偏差你通常需要离散化它，例如使用固定时间步长的 Euler–Maruyama。Walk on Spheres 算法通过观察到布朗运动在球面上的退出分布恰好是均匀分布这一事实，克服了这个问题。这让我们可以在各个球面边界之间跳跃，其中每个球都内切于区域，并尽可能大地选取以加快收敛由于游走无法精确到达边界，我们在其周围引入一个薄的 $\epsilon$ 外壳，当游走到达该外壳时即被吸收。. 整个交互式过程如下图所示。对于 Neumann 边界条件，事实证明有一种特殊的 WoS 变体，称为 *Walk on Stars*（WoSt）正如其名，这种方法“在星形上行走”，以更好地处理 Neumann（反射）边界条件。请注意，从下面的交互式图中这一点并不那么明显。原因是对于矩形，由于角点的存在，星形域看起来基本上就是球面的切片。。如果你有兴趣了解更多关于 PDE 的蒙特卡洛方法，这里（https://rohan-sawhney.github.io/mcgp-resources/）有一门很棒的课程和资源。 Neumann 障碍物：10 行走速度：1.00×ε-外壳：1e−2 平均步数： — Dirichlet（吸收）Neumann（反射） 显示跳跃位置清洁模式（动画段） 游走长度直方图（对数） 点击任意位置开始游走；后台会运行更多次游走以填充直方图。 **Walk on Spheres/Stars。**蒙特卡洛方法通过递归地采样球面（或星形），直到它们击中边界。当存在更多 Neumann 障碍物时，游走会“被困住”，需要很长时间才能到达边界，从而导致高方差和慢收敛。直方图显示了游走长度的分布；随着 Neumann 障碍物的增加，它会向右移动。 蒙特卡洛方法真是神奇！然而，正如你通过操作上面的交互式图形所看到的，随机游走在到达边界之前需要很多步。对于具有复杂几何形状和以 Neumann 为主的边界的问题尤其如此我们的手稿展示了许多这样的例子：图 1 中的仓库、图 4 中的城市或图 15 中的迷宫。。主要后果是方差和慢收敛，这使得它们在实际需要可靠解的场景中相当不实用我相信这可能就是人们一直对在实践中采用这些方法犹豫不决的原因。. 在我们的工作中，我们通过两种互补的方式来解决这个问题我应该补充一点，我们并不是第一个解决这个问题的人。已经有很多解决方法。我们方法的关键贡献和新颖之处在于，我们不是简单地提高蒙特卡洛估计器的效率或缓存解，而是提出一种将无网格蒙特卡洛方法与确定性网格求解器连接起来的方法，并利用后者的优良性质。。首先，如所示，我们可以通过将区域分解为更小的子域来缩短游走。其次，如稍后在中所介绍的，我们可以“消除”方差，但代价是（可控的）离散偏差，方法是将所有子域与确定性求解器耦合在一起，从而恢复确定性网格方法的优良“无方差”性质。 ## 3 分而治之：更短的游走 当问题复杂时，自然的解决方案是将其分解为更小、更易于管理的部分。这是数值方法中的常见策略：区域分解、多重网格、层级矩阵等等。这正是我们在论文中所采取的路径。 首先，观察到中 $\Delta u = 0$ 的一个优美性质是它在任何地方都成立，包括在整个区域 $\Omega$ 的子域上。因此，一个关键的直觉是，通过将区域分解为更小的子域，可以自然地使游走变短。 更正式地说，我们建议将区域划分为一组较小的*不重叠*子域 $\mathcal{D}=\{\Omega_i\}$，例如规则瓦片。在每个瓦片上，Dirichlet 边界是 (a) 继承自 $\partial\Omega_D$ 的物理 Dirichlet 部分和 (b) 人工 Dirichlet 部分——即与相邻瓦片的*界面*的并集。请注意，这种分解不需要任何网格划分。子域可以完全任意，可以自由地与 Neumann 边界相交，并且可以覆盖区域 $\Omega$ 外部的部分。 瓦片内的游走现在会在瓦片的边界处停止，如下图所示。你可以随意调整瓦片划分的分辨率，看看它如何影响直方图中的游走长度。 瓦片分辨率：1×1 行走速度：1.00× 平均步数： — “物理”Dirichlet（吸收）“人工”界面（吸收）Neumann（反射） 显示人工界面显示游走 游走长度直方图（对数） 游走在瓦片界面处终止；缩小瓦片，直方图向左移动。 **在分解子域中*行走*。**我们不是在整个区域上执行随机游走，而是将其分解为更小的子域。这导致更短的游走和更低的方差。 这一策略给我们带来了*Walks**in**Decomposed Subdomains*。那么为什么我们的论文标题写的是*Walks**on**Decomposed Subdomains*？还有一段路要走。单独这些游走还不会告诉我们太多：我们仍然不知道子域之间界面处的值。 这就是与网格求解器建立联系的地方。但在此之前，我们需要理解泊松核和解算子。 ## 4 泊松核与解算子 在单个瓦片 $\Omega_i$ 内部，也满足，因此知道边界值就可以完全确定内部解。具体来说，我们可以将其编码为一个线性算子 $\mathcal{H}_i$，它将边界值映射为内部值： $$\mathcal{H}_i : \partial\Omega_i \to \Omega_i $$ 换句话说，对于所有 $x \in \Omega_i$，有 $u(x) = \mathcal{H}_i[u](x)$。但我们如何得到 $\mathcal{H}_i$？毕竟它局部地就是 PDE 的解…… 思路是观察到 $\mathcal{H}_i$ 可以写成积分形式 $$ u(x) = \int_{\partial\Omega_i} P_{\Omega_i}(x, z)\, u(z)\, dz $$ 其中 $P_{\Omega_i}(x, z)$ 称为*泊松核*。关键的观察是，$P_{\Omega_i}(x, z)$ 恰好是从 $x$ 出发的布朗运动在 $z$ 处撞击边界的首达概率密度。等等——这不正是 Walk on Spheres 所计算的吗？ 完全正确，这给了我们一个简单的近似方法：对于点 $x \in \Omega_i$，我们从 $x$ 运行多次随机游走，并将它们在 $\partial\Omega_i$ 上的出口位置进行分箱。下面的交互式图形可视化了几种不同区域的泊松核，这些核就是使用这一策略制表的。 场景 正在估计 P(x, ·)… 拖动 $x$ 移动源点。拖动虚线 Neumann 障碍物以移动它；或使用滑块调整形状参数。边界上的直方图显示了泊松核 $P(x, z)$，即沿 $\partial\Omega_D$ 的首达概率密度。 **泊松核。**对于选定的内部源点 $x$，沿 $\partial\Omega_D$ 的直方图显示 $P(x, z)$，即从 $x$ 出发的布朗游走的首达密度。相应的统计量使用蒙特卡洛方法实时估计。 泊松核有一个对偶解释：不是固定 $x$ 并询问*游走在何处退出*，而是固定边界点 $z$ 并询问*哪些内部点最有可能将游走送到那里*。换句话说，Dirichlet 边界上 $z$ 处的单位点源如何影响内部——这正是解算子！ 场景 源宽：0.10 正在预计算... 将 $z$ 沿着 $\partial\Omega_D$ 拖动到任意位置。内部就是泊松核 $P(\cdot, z)$ —— 对 $z$ 处点源的响应。“源宽”平滑了源点以强调效果。 **解算子。**同一个核 $P(x, z)$，只是从 $z$ 的角度而非 $x$ 来观察。对于每个场景，核矩阵通过蒙特卡洛样本即时预计算；随后拖动 $z$ 即可立即查询它的某一列。 求解器 WoStWoS 瓦片分辨率：6×6 正在估计 P(x, ·)… 显示源点 $x$ 和直方图 点击左侧的瓦片选择一个子域；在右侧拖动源点 $x$。右侧面板显示该子域的局部泊松核。 **子域泊松核。**在分解中选择一个瓦片（左），并估计其局部首达解算子（右）。每个瓦片边——无论是物理 Dirichlet 边界还是人工界面——都作为吸收出口，而内部障碍物保持为 Neumann。 求解器 WoStWoS 瓦片分辨率：4×4分箱分辨率：4×4样本/桶：1000 选择一个子域，然后点击“预计算”。 点击左侧瓦片选择子域，然后预计算其每个内部桶的边界首达直方图（Neumann 几何内部的桶将被跳过）。拖动 $x$ —— 它会吸附到最近的桶——以检查每个预计算的核。在你点击“预计算”之前不会运行任何东西。 **子域的预计算分箱解算子。**选择分解中的一个瓦片（左）；其内部被划分为一个 $R\times R$ 的桶网格（右）。对于每个桶，蒙特卡洛短游走估计每个边上 $R$ 个边界分箱上的首达分布。拖动 $x$ 会吸附到一个桶并显示其预计算的核——即子域的离散解算子，一次显示一行。 如上文交互式图形所暗示的，有一种自然且简单的方式可以将解算子预计算并离散化为一个简单的矩阵 $\mathbf{H}$请注意，我们并不声称这是最有效的方式。这里有巨大的改进空间！。如下图中的逐行所示，$H_{ij}$ 代表从 $x_i$ 到边界的首达概率；逐列看，$H_{ij}$ 给出内部点 $x_i$ 对位于 $z_j$ 的单位源点的响应。 离散解算子 H 将包含 z_j 的边界面板映射到内部点 x_i **离散解算子。**用配置点 $\{x_i\} \subset \Omega$ 离散内部，并将 Dirichlet 边界 $\partial\Omega_D$ 划分成面板 $\{\Gamma_j\}$（配置点为 $\{z_j\}$），得到一个近似解算子的矩阵 $\mathbf{H}$。逐行看，$H_{ij}$ 是在 $x_i$ 处发起的游走通过面板 $\Gamma_j$ 退出的首达概率；逐列看，$H_{ij}$ 给出内部点 $x_i$ 对位于 $z_j$ 的局部单位源点的响应。 通过为每个子域制作一个离散解算子 $\mathbf{H}_i$