Walk on Decomposed Subdomains

Blog post by Clément Jambon. 1Elliptic PDEs and Boundary Value Problems Many real-world phenomena are governed by elliptic partial differential equations (PDEs): heat conduction, electrostatics, path planning, steady-state potential flow, and more. These are often cast as boundary value problems (BVPs), where values are prescribed on the boundary of a domain and we seek the solution to the PDE in the interior. Consider for example the Laplace equation with Dirichlet boundary conditions: $$ \begin{cases} \Delta u = 0 & \text{in } \Omega \\ u = g & \text{on } \partial\Omega_D \end{cases} $$ Feel free to play with the interactive figure below to get an intuition for what this does: Brush value +1.00 −1 (cold)+1 (hot) Presets Click and paint inside the outer boundary band to paint Dirichlet values $g$. The interior satisfies $\Delta u=0$. Dirichlet problem. Solving the Laplace equation with Dirichlet boundary conditions on a square. We can make things more interesting by considering more complex geometries and boundary conditions. For example, people are often interested in solving mixed boundary value problems with Neumann boundary conditionsNote that we will restrict ourselves to zero-Neumann conditions.: $$ \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} $$ The interactive figure below illustrates this. Notice how the isolines bend to meet the zero-Neumann obstacle at right angles to satisfy $\frac{\partial u}{\partial n} = 0$. Brush value +1.00 −1 (cold)+1 (hot) Scene Boundary presets Click and paint the Dirichlet band as before. The interior obstacle (dashed, grey) enforces $\partial u/\partial n = 0$ — isolines bend to meet it at right angles. Mixed problem. Laplace equation with Dirichlet values painted on the outer boundary and zero-Neumann geometry inside. If you look closely, you'll see that the solution is actually "pixelated". This is because it is computed with finite differences on a grid. Finite differences are very easy to understand and to implement but they don't deal very well with complex geometries, often requiring extreme grid refinement. Another common alternative is to use finite elements. The problem is that finite elements require careful mesh generation, which can be particularly challenging and time-consuming for complex geometries Imagine designing a car and having to regenerate the mesh every time you tweak the design. That would be a nightmare — and it is! , such as the city shown below. Wind around a city. This scene contains hundreds of buildings with intricate geometry. Our method characterizes their influence on wind patterns, visualized here as steady-state potential streamlines. 2Grid-Free Monte Carlo Methods Luckily, the computer graphics community has recently revived an old idea: grid-free Monte Carlo methods. The canonical algorithm underpinning this approach is the Walk on Spheres (WoS) algorithm. The intuition, from stochastic calculus, is that if you were to simulate a Brownian motion (a continuous random walk) starting from an interior point $x$, it would eventually hit the boundary at some random location $Z_\tau$. The expected value $\mathbb{E}[g(Z_\tau)]$ of the boundary condition at that random location is exactly the solution to the Dirichlet problem given by . Brownian motion and Laplace equation. A particle starts at the interior point $x$ and diffuses until it is absorbed at a random boundary location $Z_\tau$. The boundary is colored by $g$; averaging $g(Z_\tau)$ over many such walks recovers $u(x)$. Doing this independently at every interior point gives us the whole solution. With only a few walks per point the estimate is noisy, but as we average more and more, variance progressively vanishes and the smooth harmonic solution emerges as shown in the figure below. Refresh rate: 1 walks/frame slowfast low $u$ high $u$ Every interior pixel runs independent Walk-on-Spheres estimates of $u(x)=\mathbb{E}[g(Z_\tau)]$. Averaging more walks per pixel denoises the field into the harmonic solution. It runs continuously; the slider sets how fast. Number of samples: 1 Progressive Monte Carlo estimate. The same square Dirichlet problem, now solved at every pixel by Monte Carlo. Each pixel averages many random walks: with few samples the interior is noisy, and the noise gradually fades as the number of samples grows. However, simulating Brownian motion is computationally expensive and often biasedYou typically have to discretize it, for example with fixed time steps with Euler–Maruyama.. The Walk on Spheres algorithm overcomes this by observing that the exit distribution of Brownian motion on a sphere is exactly uniform. This lets us instead hop from one sphere boundary to the next, where each sphere is inscribed in the domain and chosen to be as large as possible to converge fasterSince walks cannot reach the boundary exactly, we introduce a thin $\epsilon$-shell around it where walks are absorbed when they reach it.. The entire process is illustrated in the interactive figure below. For Neumann boundary conditions, it turns out that there is a special variant of WoS called Walk on Stars (WoSt)As its name suggests, this method "walks on stars" to better handle Neumann (reflective) boundary conditions. Note that from the interactive figure below, it isn't so obvious. The reason is that with rectangles, due to corners, star-shaped domains look basically like slices of a sphere.. If you're interested in learning more about Monte Carlo methods for PDEs, there's a great course and resources available here. SolverWoSt WoS Neumann obstacles: 10 Walk speed: 1.00× ε-shell: 1e−2 Avg steps: — Dirichlet (absorbing) Neumann (reflecting) Show jump positions Clean mode (animated segments) Walk-length histogram (log) Click anywhere inside to start a walk; many more walks are run in the background to fill the histogram. Walk on Spheres/Stars. Monte Carlo methods proceed recursively by sampling spheres (or stars) until they hit the boundary. When more Neumann obstacles are present, walks get "trapped" and take a long time to hit the boundary, leading to high variance and slow convergence. The histogram shows the distribution of walk lengths; as Neumann obstacles are added, it shifts to the right. Monte Carlo methods are truly magic! However, as you can see by playing with the interactive figure above, random walks take a lot of steps before they hit the boundary. This is particularly true for problems with complex geometries and Neumann-dominated boundariesOur manuscript showcases a lot of these examples: the warehouse in Figure 1, the city in Figure 4 or the maze in Figure 15.. The main consequence is variance and slow convergence, which render them quite impractical for scenarios where reliable solutions are requiredAnd I believe this may be why people have been hesitant to adopt these methods in practice.. In our work, we address this issue in two complementary waysI should add that we're not the first to tackle this problem. There have been many ways of solving it. The key contribution and novelty of our approach is that rather than simply improving the efficiency of Monte Carlo estimators or caching solutions, we propose a way to connect grid-free Monte Carlo methods with deterministic grid-based solvers, and leverage the nice properties of the latter. . First, as shown in , we can make walks shorter by decomposing the domain into smaller subdomains. Second, as introduced later in , we can "kill" variance, at the cost of a (controllable) discretization bias, by coupling all subdomains together with a deterministic solver, recovering the nice "variance-free" properties of deterministic grid-based methods. 3Divide to Conquer: Shorter Walks When a problem is complicated, the natural solution is to break it down into smaller, more manageable pieces. This is a common strategy in numerical methods: domain decomposition, multigrid, hierarchical matrices, and so on. And this is precisely the path we also followed in our paper. First, observe that a beautiful property of $\Delta u = 0$ in is that this holds everywhere, including on subdomains of the entire domain $\Omega$. As such, a key intuition is that walks can naturally be made shorter by decomposing the domain into smaller subdomains. More formally, we propose to decompose the domain into a partition of smaller non-overlapping subdomains $\mathcal{D}=\{\Omega_i\}$, for example regular tiles. On each tile, the Dirichlet boundary is the union of (a) physical Dirichlet pieces inherited from $\partial\Omega_D$ and (b) artificial Dirichlet pieces — the interfaces with neighboring tiles. Note that this decomposition does not require any meshing. Subdomains can be totally arbitrary, are free to intersect Neumann boundaries, and can cover parts outside of the domain $\Omega$. Walks within a tile now stop at the tile's boundary as shown in the interactive figure below. Feel free to adjust the tiling resolution and see how it affects walk lengths in the histogram. SolverWoSt WoS Tiling resolution: 1×1 Walk speed: 1.00× Avg steps: — "Physical" Dirichlet (absorbing) "Artificial" interfaces (absorbing) Neumann (reflecting) Show artificial interfaces Show walks Walk-length histogram (log) Walks terminate at tile interfaces; shrink the tiles and the histogram moves to the left. Walks in Decomposed Subdomains. Rather than executing random walks across the entire domain, we decompose it into smaller subdomains. This leads to much shorter walks with lower variance. This strategy gives us Walks in Decomposed Subdomains. So why does our paper title say Walks on Decomposed Subdomains? There's still some way to go. On their own, these walks don't tell us much yet: we still don't know the values at the interfaces between subdomains. This is where the connection to grid-based solvers will come into play. But before that, we need to understand