1<!doctype html> 2<html lang="en"> 3 4<head> 5 <meta charset="utf-8" /> 6 <meta name="viewport" content="width=device-width, initial-scale=1" /> 7 <title>Walk on Decomposed Subdomains</title> 8 <link rel="preconnect" href="https://fonts.googleapis.com"> 9 <link rel="preconnect" href="https://fonts.gstatic.com" crossorigin> 10 <link 11 href="https://fonts.googleapis.com/css2?family=Inter:wght@400;500;600;700&family=JetBrains+Mono:wght@400;500&display=swap" 12 rel="stylesheet"> 13 14 <link rel="icon" type="image/svg+xml" href="favicon.svg" /> 15 16 <link rel="stylesheet" href="css/theme.css" /> 17 <link rel="stylesheet" href="css/style.css" /> 18 19 <link rel="stylesheet" href="https://cdn.jsdelivr.net/npm/[email protected]/dist/katex.min.css" /> 20
20<script defer src="https://cdn.jsdelivr.net/npm/[email protected]/dist/katex.min.js"></script>
20 21
21<script defer src="https://cdn.jsdelivr.net/npm/[email protected]/dist/contrib/auto-render.min.js" 22 onload="renderMathInElement(document.body, {delimiters:[{left:'$$',right:'$$',display:true},{left:'$',right:'$',display:false}]});"></script>
22 23</head> 24 25<body> 26 27 <header class="paper"> 28 <h1>Walk on Decomposed Subdomains: A Hybrid Monte CarloâDeterministic Solver for Elliptic PDEs</h1> 29 <div class="venue"><i>ACM Transactions on Graphics (SIGGRAPH) 2026</i></div> 30 <div class="award"><a class="award-badge" 31 href="https://blog.siggraph.org/2026/05/siggraph-2026-technical-papers-awards-best-papers-honorable-mentions-and-test-of-time.html/" 32 target="_blank" rel="noopener noreferrer">Best Paper Award</a></div> 33 <div class="authors"><a href="https://clementjambon.github.io/" target="_blank" rel="noopener noreferrer">Clément 34 Jambon</a> · <a href="https://sinabiz.github.io/" target="_blank" rel="noopener noreferrer">Mohammad Sina 35 Nabizadeh</a> · <a href="https://people.csail.mit.edu/mina/" target="_blank" rel="noopener noreferrer">Mina 36 KonakoviÄ LukoviÄ</a></div> 37 <div class="affil">Massachusetts Institute of Technology</div> 38 <div class="links"> 39 <a href="refs/wods.pdf" target="_blank" rel="noopener noreferrer">Paper</a> 40 <a href="#bibtex">BibTeX</a> 41 <a href="#blogpost">Blog post</a> 42 </div> 43 <picture> 44 <source srcset="assets/teaser.webp" type="image/webp" /> 45 <img class="teaser" src="assets/teaser.png" alt="Warehouse scene with streamlines (Fig. 1 of the WoDS paper)" /> 46 </picture> 47 <div class="abstract"> 48 <p><span class="abstract-label">Abstract.</span> Elliptic partial differential equations are ubiquitous in 49 graphics and engineering, but remain challenging to solve on complex or evolving geometries. Traditional 50 discretization schemes (e.g., FEM/FDM) provide stable, globally coupled solutions but require heavy meshing or 51 extreme refinement to accurately resolve geometric detail. In contrast, grid-free Monte Carlo methods (e.g., 52 Walk on Spheres/Stars) adapt naturally to arbitrary geometry and offer massive parallelism, but rely on long 53 random walks whose variance grows rapidly, particularly in the presence of Neumann boundaries, leading to slow 54 convergence. We introduce a hybrid approach that combines the geometric flexibility of Monte Carlo estimation 55 with deterministic global solves that do not introduce additional stochastic error. Our method decomposes the 56 domain into simple, regular subdomains and uses Monte Carlo to estimate local first-passage solution operators 57 (Poisson kernels), where walk lengths and variance are inherently controlled by the reduced spatial scale. These 58 local operators are assembled into a sparse global system whose solution is obtained via a deterministic linear 59 solve that exactly replaces simulating discrete random walks through the domain. This global solve trades 60 stochastic variance for a fixed, resolution-dependent discretization bias, yielding stable and reusable solution 61 operators. As a result, our method attains accurate, geometry-aware solutions even on coarse discretizations, 62 and enables efficient solves and re-solves by computing and updating only the local operators affected by the 63 geometry and its changes. We evaluate the approach on complex two-dimensional domains, benchmarking accuracy and 64 convergence against standard grid-free and grid-based baselines, and demonstrate applications to microstructure 65 simulation and flow-based path planning and streamline visualization.</p> 66 </div> 67 68 <h2 class="paper-section-heading">Citation</h2> 69 70 <pre class="bibtex" id="bibtex"><button class="copy-btn" data-role="copy-bib">copy</button><code>@article{wods, 71 author = {Jambon, Cl\'{e}ment and Nabizadeh, Mohammad Sina and Konakovi\'{c} Lukovi\'{c}, Mina}, 72 title = {Walk on Decomposed Subdomains: A Hybrid Monte CarloâDeterministic Solver for Elliptic PDEs}, 73 year = {2026}, 74 issue_date = {July 2026}, 75 publisher = {Association for Computing Machinery}, 76 address = {New York, NY, USA}, 77 volume = {45}, 78 number = {4}, 79 url = {https://doi.org/10.1145/3811340}, 80 doi = {10.1145/3811340}, 81 journal = {ACM Trans. Graph.}, 82 month = jul, 83 articleno = {132}, 84 numpages = {22} 85}
85</code></pre> 86 87 <p class="acknowledgments">Acknowledgments: We thank the 88 reviewers for their insightful feedback. We are grateful to Pavle KonakoviÄ for his help with the design and 89 rendering of 3D scenes. We thank Rohan Sawhney, Bailey Miller and Hamid Kamkari for valuable discussions. This 90 work was supported 91 by the Wistron Corporation, the Siebel Scholars program and the MIT Generative AI Impact Consortium.</p> 92 </header> 93 94 <article id="blogpost"> 95 96 <p class="byline">Blog post by <a href="https://clementjambon.github.io/" target="_blank" 97 rel="noopener noreferrer">Clément Jambon</a>.</p> 98 99 <aside class="callout"> 100 <p><span class="callout-label">Disclaimer.</span> The goal of this blog post is to walk you through the main 101 intuitions and ideas behind our work. To this end, it takes a completely different approach from the exposition 102 in the paper. It also takes a few technical and non-rigorous shortcuts. If you want a more formal treatment, 103 please refer directly to <a href="refs/wods.pdf" target="_blank" rel="noopener noreferrer">the paper</a>. Note 104 also that the code behind this blogpost shouldn't be treated as a reference implementation: the 105 visualizations are for illustrative purposes only<fn>Some things are “faked” to keep the webpage 106 lightweight! 107 </fn>.</p> 108 </aside> 109 110 <!-- ===== SECTION 1 ===== --> 111 <h2><span class="secnum">1</span>Elliptic PDEs and Boundary Value Problems</h2> 112 113 114 <p>Many real-world phenomena are governed by elliptic partial differential equations (PDEs): heat conduction, 115 electrostatics, path planning, steady-state potential flow, and more<cite data-key="evans2022partial"></cite>. 116 These are often cast as boundary value 117 problems (BVPs), where values are prescribed on the boundary of a domain and we seek the solution to the PDE in 118 the interior.</p> 119 120 <p>Consider for example the <a href="https://en.wikipedia.org/wiki/Laplace%27s_equation" target="_blank" 121 rel="noopener noreferrer">Laplace equation</a> with Dirichlet boundary conditions:</p> 122 123 <div class="equation" data-label="laplace"> 124 $$ \begin{cases} \Delta u = 0 & \text{in } \Omega \\ u = g & \text{on } \partial\Omega_D \end{cases} $$ 125 </div> 126 127 <p>Feel free to play with the interactive figure below to get an intuition for what this does:</p> 128 129 <figure class="figure" id="i0"> 130 <div class="figure-body"> 131 <div class="canvas-wrap"> 132 <div class="canvas-row"> 133 <div class="canvas-col"><canvas width="320" height="320" class="diagram"></canvas></div> 134 <div class="controls"> 135 <label>Brush value <span class="brush-val" data-role="brush-val">+1.00</span></label> 136 <input type="range" class="colormap-slider" data-role="brush" min="-1" max="1" step="0.01" value="1" /> 137 <div class="slider-ends"><span>â1 (cold)</span><span>+1 (hot)</span></div> 138 <button data-role="reset">Reset boundary</button> 139 <label style="margin-top:14px;">Presets</label> 140 <div class="presets" data-role="presets"></div> 141 <p style="font-size:0.78rem; color:var(--color-text-muted); margin-top:14px;">Click and paint inside the 142 outer boundary band to paint Dirichlet values $g$. The interior satisfies $\Delta u=0$.</p> 143 </div> 144 </div> 145 </div> 146 <div class="figure-caption"><b>Dirichlet problem.</b> Solving the Laplace equation with Dirichlet boundary 147 conditions on a square.</div> 148 </div> 149 </figure> 150 151 <p>We can make things more interesting by considering more complex geometries and boundary conditions. For example, 152 people are often interested in solving mixed boundary value problems with Neumann boundary conditions<fn>Note that 153 we will restrict ourselves to zero-Neumann conditions.</fn>:</p> 154 155 <div class="equation" data-label="mixed"> 156 $$ \begin{cases} \Delta u = 0 & \text{in } \Omega \\ u = g & \text{on } \partial\Omega_D \\ \frac{\partial 157 u}{\partial n} = 0 & \text{on } \partial\Omega_N \end{cases} $$ 158 </div> 159 160 <p>The interactive figure below illustrates this. Notice how the isolines bend to meet the zero-Neumann obstacle at 161 right angles to satisfy $\frac{\partial u}{\partial n} = 0$.</p> 162 163 <figure class="figure" id="i0b"> 164 <div class="figure-body"> 165 <div class="canvas-wrap"> 166 <div class="canvas-row"> 167 <div class="canvas-col"><canvas width="320" height="320" class="diagram"></canvas></div> 168 <div class="controls"> 169 <label>Brush value <span class="brush-val" data-role="brush-val">+1.00</span></label> 170 <input type="range" class="colormap-slider" data-role="brush" min="-1" max="1" step="0.01" value="1" /> 171 <div class="slider-ends"><span>â1 (cold)</span><span>+1 (hot)</span></div> 172 <div style="display:flex; gap:6px; flex-wrap:wrap;"> 173 <button data-role="reset">Reset boundary</button> 174 <button data-role="toggle-iso">Hide isolines</button> 175 </div> 176 <label style="margin-top:14px;">Scene</label>
177 <div class="presets" data-role="scene-presets"></div> 178 <label style="margin-top:14px;">Boundary presets</label> 179 <div class="presets" data-role="presets"></div> 180 <p style="font-size:0.78rem; color:var(--color-text-muted); margin-top:14px;">Click and paint the 181 Dirichlet band as before. The interior obstacle (dashed, grey) enforces $\partial u/\partial n = 0$ â 182 isolines bend to meet it at right angles.</p> 183 </div> 184 </div> 185 </div> 186 <div class="figure-caption"><b>Mixed problem.</b> Laplace equation with Dirichlet values painted on the outer 187 boundary and zero-Neumann geometry inside.</div> 188 </div> 189 </figure> 190 191 <p>If you look closely, you'll see that the solution is actually "pixelated". This is because it is computed with 192 <a href="https://en.wikipedia.org/wiki/Finite_difference_method" target="_blank" rel="noopener noreferrer">finite 193 differences</a> on a grid. Finite differences are very easy to 194 understand and to implement but they don't deal very well with complex geometries, often requiring extreme grid 195 refinement. Another common alternative is to use <a href="https://en.wikipedia.org/wiki/Finite_element_method" 196 target="_blank" rel="noopener noreferrer">finite elements</a>. The 197 problem is that finite elements require 198 careful mesh generation, which can be particularly challenging and time-consuming for complex geometries<fn> 199 Imagine designing a car and having to regenerate the mesh every time you tweak the design. That would be a 200 nightmare â and it is! 201 </fn>, such as 202 the city shown below. 203 </p> 204 205 <figure class="figure"> 206 <div class="figure-body"> 207 <!-- Definite width (not width:100%): the figure-body is 208 shrink-to-fit, and a percentage-width lazy image contributes 209 zero intrinsic width until it loads â the body collapsed to 210 ~8px and the caption wrapped into a 1000px-tall column. --> 211 <img src="assets/city.jpg" alt="" width="1440" height="810" loading="lazy" decoding="async" style="width:672px; max-width:100%; height:auto;" /> 212 <div class="figure-caption"><b>Wind around a city.</b> This scene contains hundreds of buildings with intricate 213 geometry. 214 Our method characterizes their influence on wind patterns, visualized here as steady-state potential 215 streamlines.</div> 216 </div> 217 </figure> 218 219 <!-- ===== SECTION 2 ===== --> 220 <h2><span class="secnum">2</span>Grid-Free Monte Carlo Methods</h2> 221 222 <p>Luckily, the computer graphics community has recently revived an old idea: grid-free Monte Carlo methods. The 223 canonical algorithm underpinning 224 this approach is the <i>Walk on Spheres</i> (WoS) algorithm<cite 225 data-key="muller1956continuous,sawhney2020mcgp"></cite>.</p> 226 227 <p>The intuition, from stochastic calculus, is that if 228 you were to simulate a <a href="https://en.wikipedia.org/wiki/Brownian_motion" target="_blank" 229 rel="noopener noreferrer">Brownian motion</a> (a continuous 230 random walk) starting from an interior point $x$, it would 231 eventually hit the boundary at some random location $Z_\tau$. The expected value $\mathbb{E}[g(Z_\tau)]$ of the 232 boundary condition at that 233 random 234 location is exactly the solution to the Dirichlet problem given by <a class="ref" data-ref="laplace"></a>.</p> 235 236 <figure class="figure" id="ibrown"> 237 <div class="figure-body"> 238 <div class="canvas-wrap"> 239 <div class="canvas-row solo"> 240 <div class="canvas-col" style="position:relative;"> 241 <canvas width="380" height="420" class="diagram"></canvas> 242 <div data-role="x-label" 243 style="position:absolute; pointer-events:none; font-size:0.95rem; transform:translate(-50%,-50%); white-space:nowrap; filter:drop-shadow(0 0 2px var(--color-surface)) drop-shadow(0 0 2px var(--color-surface)) drop-shadow(0 0 2px var(--color-surface));"> 244 </div> 245 <div data-role="hit-label" 246 style="position:absolute; pointer-events:none; font-size:0.95rem; opacity:0; transition:opacity 120ms; transform:translate(-50%,-50%); white-space:nowrap; filter:drop-shadow(0 0 2px var(--color-surface)) drop-shadow(0 0 2px var(--color-surface)) drop-shadow(0 0 2px var(--color-surface));"> 247 </div> 248 <div data-role="eq-label" 249 style="position:absolute; left:50%; bottom:6px; transform:translateX(-50%); font-size:0.95rem; white-space:nowrap;"> 250 </div> 251 </div> 252 </div> 253 <div class="studio-only" style="text-align:center; margin-top:12px;"> 254 <button data-role="drop" class="panel-btn">â¶ Drop particle</button> 255 <button data-role="clear" class="panel-btn" style="margin-left:8px;">â Clear</button> 256 </div> 257 </div> 258 <div class="figure-caption"><b>Brownian motion and Laplace equation.</b> A particle starts at the interior 259 point $x$ and diffuses until it is absorbed at a random boundary location $Z_\tau$. The boundary is colored 260 by $g$; averaging $g(Z_\tau)$ over many such walks recovers $u(x)$.</div> 261 </div> 262 </figure> 263 264 <p>Doing this independently at every interior point gives us the whole solution. With only a few walks per 265 point the estimate is noisy, but as we average more and more, variance progressively vanishes and the smooth 266 harmonic solution emerges as shown in the figure below.</p> 267 268 <figure class="figure" id="idenoise"> 269 <div class="figure-body"> 270 <div class="canvas-wrap"> 271 <div class="canvas-row" style="justify-content:center;"> 272 <div class="canvas-col"><canvas width="380" height="380" class="diagram"></canvas></div> 273 <!-- Studio: full authoring panel (play/pause, restart, rate). --> 274 <div class="controls studio-only"> 275 <label>Refresh rate: <span data-role="rate">1</span> walks/frame</label> 276 <input type="range" data-role="rate-slider" min="0.25" max="4" step="0.25" value="1" /> 277 <div class="slider-ends"><span>slow</span><span>fast</span></div> 278 <button data-role="play">⸠Pause</button> 279 <button data-role="restart">
279â» Restart</button> 280 <div class="legend"> 281 <span><i style="background:var(--hm-low)"></i> low $u$</span> 282 <span><i style="background:var(--hm-high)"></i> high $u$</span> 283 </div> 284 <p style="font-size:0.78rem; color:var(--color-text-muted); margin-top:14px;">Every interior pixel runs 285 independent Walk-on-Spheres estimates of $u(x)=\mathbb{E}[g(Z_\tau)]$. Averaging more walks per pixel 286 denoises the field into the harmonic solution. It runs continuously; the slider sets how fast.</p> 287 </div> 288 </div> 289 <!-- Public page: a read-only slider that self-advances to show the 290 number of Monte-Carlo samples accumulated so far in the current 291 (looping) denoise. The pace is fixed; there is no manual control. --> 292 <div class="controls blog-only" style="max-width:300px; margin:14px auto 0;"> 293 <label>Number of samples: <span data-role="sample-count">1</span></label> 294 <input type="range" class="progress-slider" data-role="progress" min="0" max="120" step="1" value="1" 295 style="pointer-events:none;" tabindex="-1" aria-hidden="true" /> 296 </div> 297 </div> 298 <div class="figure-caption"><b>Progressive Monte Carlo estimate.</b> The same square Dirichlet problem, now 299 solved at every pixel by Monte Carlo. Each pixel averages many random walks: with few samples the interior 300 is 301 noisy, and the noise gradually fades as the number of samples grows.</div> 302 </div> 303 </figure> 304 305 <p>However, simulating Brownian 306 motion is computationally expensive and often biased<fn>You typically have to discretize it, for example with 307 fixed time 308 steps with <a href="https://en.wikipedia.org/wiki/EulerâMaruyama" target="_blank" 309 rel="noopener noreferrer">EulerâMaruyama</a>.</fn>. The Walk on Spheres algorithm overcomes this by observing 310 that the 311 exit 312 distribution of Brownian motion on a sphere is exactly uniform. 313 This lets us instead hop from one sphere boundary to the next, where each sphere is inscribed in 314 the domain and chosen to be as large as possible to converge faster<fn>Since walks cannot reach the boundary 315 exactly, we introduce a thin $\epsilon$-shell around it where walks are absorbed when they reach it.</fn>.</p> 316 317 <p>The entire process is illustrated in the interactive figure below. For Neumann boundary conditions, it turns out 318 that there is a special variant of WoS called <i>Walk on Stars</i><cite data-key="sawhney2023wost"></cite> 319 (WoSt)<fn>As its name suggests, this method 320 "walks on stars" to better handle Neumann (reflective) boundary conditions. Note that from the interactive 321 figure 322 below, it isn't so obvious. The reason is that with rectangles, due to corners, star-shaped domains look 323 basically like slices of a sphere.</fn>. If you're interested in learning more about Monte Carlo methods for 324 PDEs, there's a great course and resources 325 available <a href="https://rohan-sawhney.github.io/mcgp-resources/" target="_blank" 326 rel="noopener noreferrer">here</a><cite data-key="sawhney2025star"></cite>. 327 </p> 328 329 <figure class="figure" id="i1"> 330 <div class="figure-body"> 331 <div class="canvas-wrap"> 332 <div class="canvas-row"> 333 <div class="canvas-col"><canvas width="380" height="380" class="diagram"></canvas></div> 334 <div class="controls"> 335 <div class="studio-only"> 336 <label>Solver</label> 337 <div class="solver-toggle" data-role="solver-toggle"> 338 <label><input type="radio" name="i1-solver" value="wost" checked><span>WoSt</span></label> 339 <label><input type="radio" name="i1-solver" value="wos"><span>WoS</span></label> 340 </div> 341 </div> 342 <label style="margin-top:14px;">Neumann obstacles: <span data-role="obstacles-label">10</span></label> 343 <input type="range" data-role="obstacles" min="0" max="12" step="1" value="10" /> 344 <div class="studio-only"> 345 <label style="margin-top:14px;">Walk speed: <span data-role="speed-label">1.00Ã</span></label> 346 <input type="range" data-role="speed" min="-2" max="0" step="0.02" value="
3460" /> 347 <label style="margin-top:14px;">ε-shell: <span data-role="eps-label">1eâ2</span></label> 348 <input type="range" data-role="epsilon" min="-3" max="-1" step="1" value="-2" /> 349 </div> 350 <div class="stat" data-role="avg">Avg steps: â</div> 351 <div class="legend"> 352 <span><i style="background:var(--color-dirichlet)"></i> Dirichlet (absorbing)</span> 353 <span><i style="height:0; background:none; border-top:2px dashed var(--color-neumann);"></i> Neumann 354 (reflecting)</span> 355 </div> 356 <div class="studio-only"> 357 <div class="checks"> 358 <label><input type="checkbox" data-role="show-jumps"> Show jump positions</label> 359 <label><input type="checkbox" data-role="clean-mode"> Clean mode (animated segments)</label> 360 </div> 361 </div> 362 <label style="margin-top:14px;">Walk-length histogram (log)</label> 363 <canvas width="240" height="160" class="hist"></canvas> 364 <button data-role="reset">Reset histogram</button> 365 <p style="font-size:0.78rem; color:var(--color-text-muted); margin-top:14px;">Click anywhere inside to 366 start a walk; many more walks are run in the background to fill the histogram.</p> 367 </div> 368 </div> 369 </div> 370 <div class="figure-caption"><b>Walk on Spheres/Stars.</b> Monte Carlo methods proceed recursively by sampling 371 spheres (or stars) until they hit the boundary. When more Neumann obstacles are present, walks get 372 "trapped" 373 and take a long time to hit the boundary, leading to high variance and slow convergence. The histogram shows 374 the distribution of walk lengths; as Neumann obstacles are added, it shifts to the right. 375 </div> 376 </div> 377 </figure> 378 379 <p>Monte Carlo methods are truly magic! However, as you can see by playing with the interactive figure above, random 380 walks take a lot of steps before they hit the boundary. This is particularly true for problems with complex 381 geometries and Neumann-dominated boundaries<fn>Our manuscript showcases a lot of these examples: the warehouse in 382 Figure 1, the city in Figure 4 or the maze in Figure 15.</fn>. The main consequence is variance and slow 383 convergence, which render them quite impractical for scenarios where reliable solutions are required<fn>And I 384 believe this may be why people have been hesitant to adopt these methods in practice.</fn>.</p> 385 386 <p>
386In our work, we address this issue in two complementary ways<fn>I should add that we're not the first to tackle 387 this problem. There have been many ways of solving it. The key contribution and novelty of our approach is that 388 rather than simply improving the efficiency of Monte Carlo estimators or caching solutions, we propose a way 389 to connect grid-free Monte Carlo methods with deterministic grid-based solvers, and leverage the nice 390 properties of the latter. 391 </fn> 392 . First, as shown in <a class="ref" data-ref="decompose"></a>, we can make 393 walks shorter by decomposing the domain into smaller subdomains. Second, as introduced later in <a class="ref" 394 data-ref="coupling"></a>, we 395 can "kill" variance, at the cost of a (controllable) discretization bias, by coupling all subdomains together with 396 a 397 deterministic solver, 398 recovering the nice "variance-free" properties of deterministic grid-based methods.</p> 399 400 <!-- ===== SECTION 3 ===== --> 401 <h2 data-label="decompose"><span class="secnum">3</span>Divide to Conquer: Shorter Walks</h2> 402 403 <p> 404 When a problem is complicated, the natural solution is to break it down into smaller, more manageable pieces. This 405 is a common strategy in numerical methods: domain decomposition, multigrid, hierarchical matrices, and so on. 406 And this is precisely the path we also followed in our paper. 407 </p> 408 409 <p>First, observe that a beautiful property of $\Delta u = 0$ in <a class="ref" data-ref="laplace"></a> is that this 410 holds 411 everywhere, 412 including on subdomains of 413 the entire domain $\Omega$. 414 As such, a key intuition is that walks can naturally be made shorter by decomposing the domain into smaller 415 subdomains. 416 </p> 417 418 <p>More formally, we propose to decompose the domain into a partition of smaller <i>non-overlapping</i> subdomains 419 $\mathcal{D}=\{\Omega_i\}$, for example regular tiles. On each tile, the Dirichlet boundary is the union of (a) 420 physical Dirichlet pieces 421 inherited from $\partial\Omega_D$ and (b) artificial Dirichlet pieces â the <i>interfaces</i> with neighboring 422 tiles. Note that this decomposition does not require any meshing. Subdomains can be totally arbitrary, are free 423 to intersect Neumann 424 boundaries, and can cover parts outside of the domain $\Omega$.</p> 425 426 <p>Walks within a tile now stop at the tile's 427 boundary as shown in the interactive figure below. Feel free to adjust the tiling resolution and see how it 428 affects walk lengths in the histogram. 429 </p> 430 431 <figure class="figure" id="i2"> 432 <div class="figure-body"> 433 <div class="canvas-wrap"> 434 <div class="canvas-row"> 435 <div class="canvas-col"><canvas width="380" height="380" class="diagram"></canvas></div> 436 <div class="controls"> 437 <div class="studio-only"> 438 <label>Solver</label> 439 <div class="solver-toggle" data-role="solver-toggle"> 440 <label><input type="radio" name="i2-solver" value="wost" checked><span>WoSt</span></label> 441 <label><input type="radio" name="i2-solver" value="wos"><span>WoS</span></label> 442 </div> 443 </div> 444 <label style="margin-top:14px;">Tiling resolution: <span data-role="tiles-label">1Ã1</span></label> 445 <input type="range" data-role="tiles" min="1" max="12" step="1" value="6" /> 446 <div class="studio-only"> 447 <label style="margin-top:14px;">Walk speed: <span data-role="speed-label">1.00Ã</span></label> 448 <input type="range" data-role="speed" min="-2" max="0" step="0.02" value="
4480" /> 449 </div> 450 <div class="stat" data-role="avg">Avg steps: â</div> 451 <div class="legend"> 452 <span><i style="background:var(--color-dirichlet)"></i> "Physical" Dirichlet (absorbing)</span> 453 <span><i style="background:var(--color-interface)"></i> "Artificial" interfaces (absorbing)</span> 454 <span><i style="height:0; background:none; border-top:2px dashed var(--color-neumann);"></i> Neumann 455 (reflecting)</span> 456 </div> 457 <div class="studio-only"> 458 <div class="checks"> 459 <label><input type="checkbox" data-role="show-interfaces" checked> Show artificial interfaces</label> 460 <label><input type="checkbox" data-role="show-walks" checked> Show walks</label> 461 </div> 462 </div> 463 <label style="margin-top:14px;">Walk-length histogram (log)</label> 464 <canvas width="240" height="160" class="hist"></canvas> 465 <button data-role="reset">Reset histogram</button> 466 <p style="font-size:0.78rem; color:var(--color-text-muted); margin-top:14px;">Walks terminate at tile 467 interfaces; shrink the tiles and the histogram moves to the left.</p> 468 </div> 469 </div> 470 </div> 471 <div class="figure-caption"><b>Walks <i>in</i> Decomposed Subdomains.</b> Rather than executing random walks 472 across the 473 entire domain, we decompose it into smaller subdomains. This leads to much shorter walks with lower variance. 474 </div> 475 </div> 476 </figure> 477 478 <p>This strategy gives us <i>Walks <b>in</b> Decomposed Subdomains</i>. So why does our paper title say <i>Walks 479 <b>on</b> Decomposed Subdomains</i>? There's still some way to go. On their own, these 480 walks don't tell us much 481 yet: we still don't 482 know the values at the interfaces between subdomains.</p> 483 484 <p>This is where the connection to grid-based solvers will come into play. But before that, we need to understand 485 Poisson kernels and solution operators.</p> 486 487 <!-- ===== SECTION 4 ===== --> 488 <h2 data-label="poisson-kernel"><span class="secnum">4</span>Poisson Kernels and Solution Operators</h2> 489 490 <p>Within a single tile $\Omega_i$, <a class="ref" data-ref="laplace"></a> is also satisfied, so knowing the 491 boundary 492 values fully determines 493 the interior solution. 494 Concretely, we can encode this as a linear operator $\mathcal{H}_i$ that maps boundary values to interior values: 495 </p> 496 497 <div class="equation" data-label="H-operator"> 498 $$\mathcal{H}_i : \partial\Omega_i \to \Omega_i $$ 499 </div> 500 501 <p>In other words, $u(x) = \mathcal{H}_i[u](x)$ for all $x \in \Omega_i$. But how do we get $\mathcal{H}_i$? It's 502 locally the solution of the PDE after all...</p> 503 504 <p>The idea is to observe that $\mathcal{H}_i$ can be written in integral form</p> 505 506 <div class="equation" data-label="poisson-integral"> 507 $$ u(x) = \int_{\partial\Omega_i} P_{\Omega_i}(x, z)\, u(z)\, dz $$ 508 </div> 509 510 <p>where $P_{\Omega_i}(x, z)$ is called the <i>Poisson kernel</i>. The key observation is that $P_{\Omega_i}(x, z)$ 511 is 512 exactly the 513 first-passage probability density of a Brownian motion starting at $x$ and hitting the boundary at $z$. Wait â 514 isn't that precisely what Walk on Spheres computes?</p> 515 516 <p>Exactly, and that gives us an easy recipe to approximate it: for a point $x \in \Omega_i$, we run many random 517 walks from 518 $x$ and bin where they exit on $\partial\Omega_i$. The interactive figure below visualizes the Poisson kernel for 519 various domains tabulated using this strategy.</p> 520 521 <figure class="figure" id="i3"> 522 <div class="figure-body"> 523 <div class="canvas-wrap"> 524 <div class="canvas-row"> 525 <div class="canvas-col"><canvas width="480" height="480" class="diagram"></canvas></div> 526 <div class="controls"> 527 <div class="studio-only"> 528 <label>Solver</label> 529 <div class="solver-toggle" data-role="solver-toggle"> 530 <label><input type="radio" name="i3-solver" value="wost" checked><span>WoSt</span></label> 531 <label><input type="radio" name="i3-solver" value="wos"><span>WoS</span></label> 532 </div> 533 </div> 534 <label style="margin-top:14px;">Scene</label>
535 <div class="presets" data-role="shape-presets"></div> 536 <div data-role="shape-params" style="margin-top:10px;"></div> 537 <div class="stat" data-role="status" style="margin-top:10px;">Estimating P(x, ·)â¦</div> 538 <button data-role="recompute">Add more MC samples</button> 539 <p style="font-size:0.78rem; color:var(--color-text-muted); margin-top:14px;">Drag $x$ to move the source. 540 Drag the dashed Neumann obstacle to move it; or adjust shape parameters with the sliders. 541 Histograms on the boundary show the Poisson kernel $P(x, z)$, i.e., the first-passage probability 542 density 543 along 544 $\partial\Omega_D$. 545 </p> 546 </div> 547 </div> 548 </div> 549 <div class="figure-caption"><b>Poisson kernel.</b> For a chosen interior source $x$, the 550 histograms along $\partial\Omega_D$ show $P(x, z)$, the first-passage density of a Brownian walk from $x$. 551 The corresponding statistics are estimated in real-time using Monte Carlo.</div> 552 </div> 553 </figure> 554 555 <p>The Poisson kernel has a dual interpretation: instead of fixing $x$ and asking <i>where the walk exits</i>, fix a 556 boundary point $z$ and ask <i>which interior points are most likely to send walks there</i>. In other words, how 557 a unit point source at $z$ on the Dirichlet boundary affects the interior â this is precisely the solution 558 operator!</p> 559 560 <figure class="figure" id="i3b"> 561 <div class="figure-body"> 562 <div class="canvas-wrap"> 563 <div class="canvas-row"> 564 <div class="canvas-col"><canvas width="480" height="480" class="diagram"></canvas></div> 565 <div class="controls"> 566 <div class="studio-only"> 567 <label>Solver</label> 568 <div class="solver-toggle" data-role="solver-toggle"> 569 <label><input type="radio" name="i3b-solver" value="wost" checked><span>WoSt</span></label> 570 <label><input type="radio" name="i3b-solver" value="wos"><span>WoS</span></label> 571 </div> 572 </div> 573 <label style="margin-top:14px;">Scene</label> 574 <div class="presets" data-role="scene-presets"></div> 575 <label style="margin-top:14px;">Source width: <span data-role="sigma-label">0.10</span></label> 576 <input type="range" data-role="sigma" min="0.03" max="0.30" step="0.005" value="0.10" /> 577 <div class="stat" data-role="status" style="margin-top:10px;">Precomputingâ¦</div> 578 <button data-role="recompute">Add more MC samples</button> 579 <p style="font-size:0.78rem; color:var(--color-text-muted); margin-top:14px;">Drag $z$ anywhere along 580 $\partial\Omega_D$. The interior is the Poisson kernel $P(\cdot, z)$ â the response to a point source at 581 $z$. The "source width" smooths the source to emphasize the effect.</p> 582 </div> 583 </div> 584 </div> 585 <div class="figure-caption"><b>Solution operator.</b> The same kernel 586 $P(x, z)$, viewed in $z$ instead of $x$. For each scene, the kernel matrix is precomputed on-the-fly with 587 Monte 588 Carlo samples; 589 dragging $z$ then queries a column of it instantly.</div> 590 </div> 591 </figure> 592 593 <figure class="figure studio-only" id="subkernel"> 594 <div class="figure-body"> 595 <div class="canvas-wrap"> 596 <div class="canvas-row" style="justify-content:center;"> 597 <div class="canvas-col"><canvas width="360" height="360" class="diagram" data-role="decomp"></canvas></div> 598 <div class="canvas-col"><canvas width="360" height="360" class="diagram" data-role="kernel"></canvas></div> 599 <div class="controls"> 600 <label>Solver</label> 601 <div class="solver-toggle" data-role="solver-toggle"> 602 <label><input type="radio" name="subk-solver" value="wost" checked><span>WoSt</span></label> 603 <label><input type="radio" name="subk-solver" value="wos"><span>WoS</span></label> 604 </div> 605 <label style="margin-top:14px;">Tiling resolution: <span data-role="tiles-label">6Ã6</span></label> 606 <input type="range" data-role="tiles" min="2" max="10" step="1" value="6" /> 607 <div class="stat" data-role="status">Estimating P(x, ·)â¦</div> 608 <button data-role="recompute">Add more MC samples</button> 609 <div class="checks"> 610 <label><input type="checkbox" data-role="show-kernel" checked> Show source $x$ & histogram</label> 611 </div> 612 <p style="font-size:0.78rem; color:var(--color-text-muted); margin-top:14px;">Click a tile on the 613 left to pick a subdomain; drag the source $x$ on the right. The right panel shows that subdomain's 614 local Poisson kernel.</p> 615 </div> 616 </div> 617 </div> 618 <div class="figure-caption"><b>Subdomain Poisson kernel.</b> Picking one tile of the decomposition (left) 619 and estimating its local first-passage solution operator (right). Every tile edge â whether physical 620 Dirichlet boundary or an artificial interface â acts as an absorbing exit, while interior obstacles stay 621 Neumann.</div> 622 </div> 623 </figure> 624 625 <figure class="figure studio-only" id="binop"> 626 <div class="figure-body"> 627 <div class="canvas-wrap"> 628 <div class="canvas-row" style="justify-content:center;"> 629 <div class="canvas-col"><canvas width="360" height="360" class="diagram" data-role="decomp"></canvas></div> 630 <div class="canvas-col"><canvas width="360" height="360" class="diagram" data-role="operator"></canvas></div> 631 <div class="controls" 632 style="flex:1 1 480px; max-width:none; display:grid; grid-template-columns:1fr 1fr; gap:0 22px; align-items:start;"> 633 <div> 634 <label>Solver</label> 635 <div class="solver-toggle" data-role="solver-toggle"> 636 <label><input type="radio" name="binop-solver" value="wost" checked><span>WoSt</span></label> 637 <label><input type="radio" name="binop-solver" value="wos"><span>WoS</span></label> 638 </div> 639 <label style="margin-top:14px;">Tiling resolution: <span data-role="tiles-label">4Ã4</span></label> 640 <input type="range" data-role="tiles" min="2" max="10" step="1" value="4" /> 641 <label style="margin-top:14px;">Binning resolution: <span data-role="res-label">4Ã4</span></label> 642 <input type="range" data-role="res" min="2" max="10" step="1" value="4" /> 643 <label style="margin-top:14px;">Samples / bucket: <span data-role="samples-label">1000</span></label> 644 <input type="range" data-role="samples" min="100" max="10000" step="100" value="1000" /> 645 </div> 646 <div> 647 <button data-role="precompute">Precompute operator</button> 648 <div data-role="progress" style="display:none; margin-top:10px;"> 649 <div style="height:8px; background:var(--color-rule); border-radius:4px; overflow:hidden;"> 650 <div data-role="progress-fill" style="height:100%; width:0%; background:var(--color-accent);"></div> 651 </div> 652 <div data-role="progress-text" 653 style="font-size:0.72rem; color:var(--color-text-muted); margin-top:4px;"></div> 654 </div> 655 <div class="stat" data-role="status" style="margin-top:10px;">Pick a subdomain, then click âPrecomputeâ. 656 </div> 657 <div style="display:flex; gap:6px; flex-wrap:wrap; margin-top:10px;"> 658 <button data-role="toggle-grid">Hide subgrid</button> 659 <button data-role="toggle-kernel">Hide source & kernel</button> 660 </div> 661 <p style="font-size:0.78rem; color:var(--color-text-muted); margin-top:14px;">Click a tile on the left 662 to 663 choose a subdomain, then Precompute its boundary first-passage histogram for every interior bucket 664 (buckets inside the Neumann geometry are skipped). Drag $x$ â it snaps to the nearest bucket â to 665 inspect each precomputed kernel. Nothing runs until you click Precompute.</p> 666 </div> 667 </div> 668 </div> 669 </div> 670 <div class="figure-caption"><b>Precomputed binned solution operator of a subdomain.</b> Pick a tile of the 671 decomposition (left); its interior is split into an $R\times R$ grid of buckets (right). For each bucket, 672 short Monte-Carlo walks estimate the first-passage distribution over the $R$ boundary bins per edge. 673 Dragging $x$ snaps to a bucket and shows its precomputed kernel â the subdomain's discrete solution 674 operator, one row at a time.</div> 675 </div> 676 </figure> 677 678 <p>As implied by the interactive figures above, there's a natural and simple way of precomputing and discretizing 679 solution operators as a simple matrix $\mathbf{H}$<fn>Note that we do not claim that this is the most efficient 680 way. There's huge room for improvement here!</fn>. As shown in the figure below, row-wise, $H_{ij}$ represents 681 first- 682 passage probabilities from $x_i$ to the boundary; column-wise, $H_{ij}$ gives the 683 interior response to a unit source at $z_j$.</p> 684 685 <figure class="figure" id="hmatrix-fig"> 686 <div class="figure-body"> 687 <div class="canvas-wrap" style="display:flex; justify-content:center;"> 688 <!-- Definite width for the same shrink-to-fit reason as the 689 city image above. --> 690 <img src="assets/discrete_solution_operator.svg" 691 alt="Discrete solution operator H mapping a boundary panel c
691ontaining z_j to interior point x_i" 692 width="398" height="392" loading="lazy" decoding="async" 693 style="width:420px; max-width:100%; height:auto;" /> 694 </div> 695 <div class="figure-caption"><b>Discrete solution operator.</b> 696 Discretizing the interior with collocation points $\{x_i\} \subset \Omega$ and the Dirichlet boundary 697 $\partial\Omega_D$ into panels $\{\Gamma_j\}$ with collocation points $\{z_j\}$ yields a matrix 698 $\mathbf{H}$ that approximates the solution operator. Row-wise, $H_{ij}$ is the first-passage 699 probability that a walk launched at $x_i$ exits through panel $\Gamma_j$; column-wise, $H_{ij}$ gives 700 the interior response at $x_i$ to a unit source localized at $z_j$.</div> 701 </div> 702 </figure> 703 704 <p>By tabulating one discrete solution operator $\mathbf{H}_i$ for each tile $\Omega_i$, we can solve the discrete 705 Dirichlet problem within each tile by a simple matrix-vector 706 multiplication. However, 707 this still doesn't tell us how to find the values at the interfaces between tiles. Enter absorbing Markov chains! 708 </p> 709 710 <!-- ===== SECTION 5 ===== --> 711 <h2 data-label="markov"><span class="secnum">5</span>Absorbing Markov Chains</h2> 712 713 <p>An <a href="https://en.wikipedia.org/wiki/Absorbing_Markov_chain" target="_blank" 714 rel="noopener noreferrer">absorbing Markov chain</a> is also a random 715 walk, but this time on a discrete (and finite) set of states 716 $\mathcal{S}$. 717 Some states $\mathcal{T}$ are called <b>transient</b> (the walker may pass through them, possibly many times) and 718 the rest $\mathcal{A}$ are 719 <b>absorbing</b> (once entered, the walker never leaves), such that $\mathcal{S} = \mathcal{T} \sqcup 720 \mathcal{A}$. The figure below provides a simple example. 721 </p> 722 723 <figure class="figure" id="amc-anim"> 724 <div class="figure-body"> 725 <div class="canvas-wrap"> 726 <div class="amc-stage" data-role="stage"></div> 727 <!-- Filled by auto_markov.js at (lazy) init; present in the 728 markup so the row's height is reserved before then. --> 729 <div class="amc-legend" data-role="legend"></div> 730 </div> 731 <div class="figure-caption"><b>Absorbing Markov chain.</b> Three transient states 732 $\mathcal{T}=\{t_1,t_2,t_3\}$ sit between two absorbing states $\mathcal{A}=\{a_L,a_R\}$. At each step a 733 walker hops to a random neighbor; the absorbing states carry a self-loop, so once a walker lands there it 734 stays forever.</div> 735 </div> 736 </figure> 737 738 <p>A key observation is that we can also define a boundary value problem on the absorbing Markov chain analogous to 739 <a class="ref" data-ref="laplace"></a>. Concretely, as shown in the interactive example below, we can prescribe 740 values on absorbing states, run walks from transient states until they are absorbed, and average the absorbing 741 values to get a solution defined on the transient states, i.e., 742 $$u(t) = \mathbb{E}[g(A_\tau) \mid X_0 = t]$$ 743 for any transient state $t \in \mathcal{T}$, where $A_\tau \in \mathcal{A}$ is the absorbing state where the 744 walk ends and $g$ holds the prescribed values at absorbing states. This is very reminiscent of the Walk on Spheres 745 algorithm, isn't it? 746 </p> 747 748 <p>The interactive figure below provides an example for a simple Markov chain linking two 749 absorbing endpoints.</p> 750 751 <figure class="figure" id="imc"> 752 <div class="figure-body"> 753 <div class="canvas-wrap"> 754 <div style="display:flex; flex-direction:column; align-items:center; gap:18px;"> 755 <div style="display:flex; align-items:center; gap:10px;"> 756 <div 757 style="display:flex; flex-direction:column; align-items:center; gap:4px; font-family:var(--font-mono); font-size:0.72rem; color:var(--color-text-muted);"> 758 <span>+1</span> 759 <span class="colormap-slider-vwrap"> 760 <input type="range" class="colormap-slider" data-role="gL" min="-1" max="1" step="0.01" value="-1" /> 761 </span>
762 <span>â1</span> 763 </div> 764 <div class="canvas-col"><canvas width="580" height="170" class="diagram"></canvas></div> 765 <div 766 style="display:flex; flex-direction:column; align-items:center; gap:4px; font-family:var(--font-mono); font-size:0.72rem; color:var(--color-text-muted);"> 767 <span>+1</span> 768 <span class="colormap-slider-vwrap"> 769 <input type="range" class="colormap-slider" data-role="gR" min="-1" max="1" step="0.01" value="1" /> 770 </span> 771 <span>â1</span> 772 </div> 773 </div> 774 <div class="controls" 775 style="display:grid; grid-template-columns: repeat(2, 1fr); gap:24px; align-items:start; max-width:580px; width:100%;"> 776 <div style="display:flex; flex-wrap:wrap; gap:6px; align-content:flex-start;"> 777 <div class="stat" data-role="status" style="flex-basis:100%;">walks: 0</div> 778 <button data-role="play" style="width:auto;">Run walks</button> 779 <button data-role="step" style="width:auto;">Launch one walk</button> 780 <button data-role="show-exact" style="width:auto;">Show exact solve</button> 781 <button data-role="reset" style="width:auto;">Reset MC</button> 782 <button data-role="reset-bias" style="width:auto;">Reset probabilities</button> 783 </div> 784 <p style="font-size:0.78rem; color:var(--color-text-muted); margin:0;">Drag the vertical sliders next 785 to each absorbing endpoint to set its boundary value. Click on any transient node to start a walk. 786 Drag the slider beneath each node to change its transition probability. The <b>inner color</b> of 787 each interior circle is the running Monte Carlo estimate. You can launch many walks with "Run 788 walks" or show the exact solution with "Show exact solve".</p> 789 </div> 790 </div> 791 </div> 792 <div class="figure-caption"><b>Discrete boundary value problem.</b> By launching walks from transient states and 793 accumulating the values obtained at absorbing states, we can approximate the solution to a discrete 794 boundary value problem.</div> 795 </div> 796 </figure> 797 798 <p>Note how we're absolutely free to choose arbitrary transition probabilities and how they influence the 799 solution<fn>Spoiler: in the next section, we'll choose the geometrically-informed values given by Poisson kernels 800 as transition probabilities! 801 </fn>.</p> 802 803 <p>You may have also noticed that, even with a handful of states, you need to launch quite a few walks before the 804 Monte Carlo estimate stops wiggling around. So here is the natural question: do we really have to <i>simulate</i> 805 random walks?</p> 806 807 <p>Look at any transient state $i$. By the <a href="https://en.wikipedia.org/wiki/Memorylessness" target="_blank" 808 rel="noopener noreferrer">memoryless 809 property</a> of Markov chains, a walker sitting at $i$ takes a single step to a 810 neighbor $j$ with probability $P_{ij}$, and from there the rest of the walk is statistically identical to a 811 fresh walk launched at $j$. In other words, the expected absorbed value at $i$ is just the weighted average of 812 the expected absorbed values at its neighbors:</p> 813 814 <div class="equation"> 815 $$ \begin{aligned} 816 u(i) &\;=\; \sum_{j \in \mathcal{S}} P_{ij}\, u(j) && \text{on } \mathcal{T}, \\ 817 u(a) &\;=\; g(a) && \text{on } \mathcal{A}. 818 \end{aligned} $$ 819 </div> 820 821 <p>This is exactly a discrete analog of the Laplace equation <a class="ref" data-ref="laplace"></a>: <i>each 822 interior value is the average of its neighbors</i>, with prescribed values on the boundary<fn>In the continuous 823 case, this mean-value property is locally captured by the Laplacian operator being zero. But note that it also
824 holds in an integral sense as $u(x) = \frac{1}{|\partial B(x, R)|}\int_{\partial B(x, R)} u(z)\,dz$ for a 825 ball centered at $x$ of radius $R$.</fn>. 826 Random walks 827 have quietly disappeared: we're left with a system of linear equations where the unknowns are the transient 828 values.</p> 829 830 <p>To make this concrete, split the transition matrix into transient-to-transient and transient-to-absorbing 831 blocks,</p> 832 833 <div class="equation"> 834 $$ \mathbf{P} \;=\; \begin{bmatrix} \mathbf{Q} & \mathbf{R} \\ \mathbf{0} & \mathbf{I} \end{bmatrix}, $$ 835 </div> 836 837 <p>and collect the unknown transient values into a vector $\mathbf{u}_{\mathcal{T}}$ and the prescribed absorbing 838 values into $\mathbf{g}$. The averaging identity above becomes 839 $\mathbf{u}_{\mathcal{T}} = \mathbf{Q}\,\mathbf{u}_{\mathcal{T}} + \mathbf{R}\,\mathbf{g}$, i.e.,</p> 840 841 <div class="equation"> 842 $$ (\mathbf{I} - \mathbf{Q})\, \mathbf{u}_{\mathcal{T}} \;=\; \mathbf{R}\, \mathbf{g}. $$ 843 </div> 844 845 <p>That's it. As long as every walk is eventually absorbed (which it is, with probability one), $\mathbf{I} - 846 \mathbf{Q}$ is invertible, and a single linear solve hands us the <i>exact</i> expected absorbed value at every 847 transient state at once: no sampling, no variance, no waiting for the estimator to settle. Try it in the figure 848 above with "Show exact solve": the colors should snap to their final values immediately!</p> 849 850 851 <p>One last fun thing for the road! It turns out that you can see $\mathbf{I}-\mathbf{P}$ as a random-walk 852 Laplacian<cite data-key="chung1997spectral"></cite>. In 1D, if we choose a symmetric random walk, the 853 corresponding random-walk Laplacian should be very 854 familiar to you: it's exactly the <a 855 href="https://en.wikipedia.org/wiki/Compact_stencil#Three_Point_Stencil_Example" target="_blank" 856 rel="noopener noreferrer">three-point finite difference Laplacian</a> in 1D (up to a scaling factor).</p> 857 858 <figure class="figure" id="rwlap-fig"> 859 <div class="figure-body"> 860 <div class="canvas-wrap" style="display:flex; justify-content:center;"> 861 <div class="svg-label-wrap" style="max-width:460px;"> 862 <svg viewBox="0 0 460 150" style="display:block; width:100%; height:auto;" 863 xmlns="http://www.w3.org/2000/svg"> 864 <!-- baseline --> 865 <line x1="40" y1="95" x2="420" y2="95" stroke="#888" stroke-width="1" /> 866 <!-- nodes --> 867 <g> 868 <circle cx="80" cy="95" r="6" fill="#bbb" stroke="#444" /> 869 <circle cx="160" cy="95" r="6" fill="#bbb" stroke="#444" /> 870 <circle cx="230" cy="95" r="7" fill="#e07a5f" stroke="#444" /> 871 <circle cx="300" cy="95" r="6" fill="#bbb" stroke="#444" /> 872 <circle cx="380" cy="95" r="6" fill="#bbb" stroke="#444" /> 873 </g> 874 <!-- arrows: 1/2 left, 1/2 right from node i --> 875 <defs> 876 <marker id="arr" markerWidth="8" markerHeight="8" refX="7" refY="4" orient="auto"> 877 <path d="M0,0 L8,4 L0,8 z" fill="#2a7" /> 878 </marker> 879 </defs> 880 <path d="M222,80 Q195,45 168,80" fill="none" stroke="#2a7" stroke-width="1.6" marker-end="url(#arr)" /> 881 <path d="M238,80 Q265,45 292,80" fill="none" stroke="#2a7" stroke-width="1.6" marker-end="url(#arr)" /> 882 </svg> 883 <!-- KaTeX labels overlaid as positioned HTML (Safari-safe; avoids foreignObject in scaled SVG) --> 884 <div class="svg-label" style="left:17.39%; top:78.67%;">$i-2$</div> 885 <div class="svg-label" style="left:34.78%; top:78.67%;">$i-1$</div> 886 <div class="svg-label" style="left:50%; top:78.67%;">$i$</div> 887 <div class="svg-label" style="left:65.22%; top:78.67%;">$i+1$</div> 888 <div class="svg-label" style="left:82.61%; top:78.67%;">$i+2$</div> 889 <div class="svg-label" style="left:42.39%; top:92%; font-size:12px; color:var(--color-text-muted);">$h$ 890 </div> 891 <div class="svg-label" style="left:57.61%; top:92%; font-size:12px; color:var(--color-text-muted);">$h$ 892 </div> 893 <div class="svg-label" style="left:42.39%; top:21.33%; color:#2a7;">$\tfrac{1}{2}$</div> 894 <div class="svg-label" style="left:57.61%; top:21.33%; color:#2a7;">$\tfrac{1}{2}$</div> 895 </div> 896 </div> 897 <div class="equation"> 898 $$ \begin{aligned} 899 \big[(\mathbf{I}-\mathbf{P})\,\mathbf{u}\big]_i 900 &\;=\; u_i - \tfrac{1}{2}u_{i-1} - \tfrac{1}{2}u_{i+1} \\ 901 &\;=\; -\tfrac{h^2}{2}\,\underbrace{\frac{u_{i-1} - 2u_i + u_{i+1}}{h^2}}_{\Delta_h u_i}. 902 \end{aligned} $$ 903 </div> 904 <div class="figure-caption"><b>1D random-walk Laplacian.</b> For a symmetric random walk on a uniform grid of 905 spacing $h$, the operator $\mathbf{I}-\mathbf{P}$ is exactly the standard 906 three-point finite-difference Laplacian $\Delta_h$, up to the rescaling $-h^2/2$.</div> 907 </div> 908 </figure> 909 910 <p>And in 2D, the random-walk Laplacian for a symmetric walk on a regular grid recovers the standard <a 911 href="https://en.wikipedia.org/wiki/Five-point_stencil" target="_blank" rel="noopener noreferrer">five-point 912 finite-difference stencil</a>.</p> 913 914 <figure class="figure" id="rwlap-fig-2d"> 915 <div class="figure-body"> 916 <div class="canvas-wrap" style="display:flex; justify-content:center;"> 917 <div class="svg-label-wrap" style="max-width:460px;"> 918 <svg viewBox="0 0 460 320" style="display:block; width:100%; height:auto;" 919 xmlns="http://www.w3.org/2000/svg"> 920 <!-- grid lines --> 921 <g stroke="#ddd" stroke-width="1"> 922 <line x1="60" y1="60" x2="400" y2="60" />
923 <line x1="60" y1="160" x2="400" y2="160" /> 924 <line x1="60" y1="260" x2="400" y2="260" /> 925 <line x1="90" y1="40" x2="90" y2="280" /> 926 <line x1="230" y1="40" x2="230" y2="280" /> 927 <line x1="370" y1="40" x2="370" y2="280" /> 928 </g> 929 <!-- arrows: 1/4 to each of 4 neighbors --> 930 <defs> 931 <marker id="arr2d" markerWidth="8" markerHeight="8" refX="7" refY="4" orient="auto"> 932 <path d="M0,0 L8,4 L0,8 z" fill="#2a7" /> 933 </marker> 934 </defs> 935 <path d="M222,160 Q175,135 100,160" fill="none" stroke="#2a7" stroke-width="1.6" 936 marker-end="url(#arr2d)" /> 937 <path d="M238,160 Q285,135 360,160" fill="none" stroke="#2a7" stroke-width="1.6" 938 marker-end="url(#arr2d)" /> 939 <path d="M230,152 Q205,105 230,70" fill="none" stroke="#2a7" stroke-width="1.6" 940 marker-end="url(#arr2d)" /> 941 <path d="M230,168 Q255,215 230,250" fill="none" stroke="#2a7" stroke-width="1.6" 942 marker-end="url(#arr2d)" /> 943 <!-- nodes --> 944 <g> 945 <circle cx="90" cy="60" r="5" fill="#bbb" stroke="#444" /> 946 <circle cx="230" cy="60" r="6" fill="#bbb" stroke="#444" /> 947 <circle cx="370" cy="60" r="5" fill="#bbb" stroke="#444" /> 948 <circle cx="90" cy="160" r="6" fill="#bbb" stroke="#444" /> 949 <circle cx="230" cy="160" r="7" fill="#e07a5f" stroke="#444" /> 950 <circle cx="370" cy="160" r="6" fill="#bbb" stroke="#444" /> 951 <circle cx="90" cy="260" r="5" fill="#bbb" stroke="#444" /> 952 <circle cx="230" cy="260" r="6" fill="#bbb" stroke="#444" /> 953 <circle cx="370" cy="260" r="5" fill="#bbb" stroke="#444" /> 954 </g> 955 </svg> 956 <!-- KaTeX labels overlaid as positioned HTML (Safari-safe; avoids foreignObject in scaled SVG) --> 957 <div class="svg-label" style="left:50%; top:14.06%;">$i,\,j+1$</div> 958 <div class="svg-label" style="left:16.52%; top:50.94%; transform:translate(-100%,-50%);">$i-1,\,j$</div> 959 <div class="svg-label" style="left:53.04%; top:55.31%; transform:translate(0,-50%);">$i,\,j$</div> 960 <div class="svg-label" style="left:83.48%; top:50.94%; transform:translate(0,-50%);">$i+1,\,j$</div> 961 <div class="svg-label" style="left:50%; top:86.56%;">$i,\,j-1$</div> 962 <div class="svg-label" style="left:34.78%; top:94.06%; font-size:12px; color:var(--color-text-muted);">$h$ 963 </div> 964 <div class="svg-label" style="left:65.22%; top:94.06%; font-size:12px; color:var(--color-text-muted);">$h$ 965 </div> 966 <div class="svg-label" style="left:34.78%; top:39.69%; color:#2a7;">$\tfrac{1}{4}$</div> 967 <div class="svg-label" style="left:65.22%; top:39.69%; color:#2a7;">$\tfrac{1}{4}$</div> 968 <div class="svg-label" style="left:42.39%; top:30.31%; color:#2a7;">$\tfrac{1}{4}$</div> 969 <div class="svg-label" style="left:57.61%; top:67.81%; color:#2a7;">$\tfrac{1}{4}$</div> 970 </div> 971 </div> 972 <div class="equation"> 973 $$ \begin{aligned} 974 \big[(\mathbf{I}-\mathbf{P})\,\mathbf{u}\big]_{i,j} 975 &\;=\; u_{i,j} - \tfrac{1}{4}\big(u_{i-1,j} + u_{i+1,j} + u_{i,j-1} + u_{i,j+1}\big) \\ 976 &\;=\; -\tfrac{h^2}{4}\,\underbrace{\frac{u_{i-1,j} + u_{i+1,j} + u_{i,j-1} + u_{i,j+1} - 977 4u_{i,j}}{h^2}}_{\Delta_h u_{i,j}}. 978 \end{aligned} $$ 979 </div> 980 <div class="figure-caption"><b>2D random-walk Laplacian.</b> For a symmetric random walk on a uniform square 981 grid of spacing $h$, the operator $\mathbf{I}-\mathbf{P}$ is exactly the standard five-point 982 finite-difference Laplacian $\Delta_h$, up to the rescaling $-h^2/4$.</div> 983 </div> 984 </figure> 985 986 <p>From there, you probably see the pattern. What if instead of these canonical probabilities, we considered more 987 general transition probabilities based on the geometry of the subdomains?</p> 988 989 990 <!-- ===== SECTION 6 ===== --> 991 <h2 data-label="coupling"><span class="secnum">6</span>Coupling Tiles via an Absorbing Markov Chain</h2> 992 993 <p>
993Now we have almost all the pieces to recover values at interfaces between subdomains! The trick is to <i>Walk 994 <b>on</b> 995 Decomposed Subdomains</i>, or more precisely on their interfaces. 996 </p> 997 998 <p> 999 The last problem is that, so far, we've only seen 1000 in <a class="ref" data-ref="poisson-kernel"></a> how to jump from the 1001 inside of a subdomain to one of its interfaces and not from one interface to another. 1002 The fix is to view each interface through <i>another decomposition</i><fn>This time it is an overlapping 1003 decomposition.</fn> of the domain, one where the 1004 interface sits in the interior of a subdomain. On a grid, that subdomain is just the two cells sharing the 1005 edge, and its boundary is made of the neighboring interfaces (or pieces of the global boundary). We call 1006 these <i>âco-edgeâ subdomains</i>. 1007 </p> 1008 1009 <figure class="figure" id="coedge-fig"> 1010 <div class="figure-body"> 1011 <div class="canvas-wrap"> 1012 <div class="canvas-row" style="justify-content:center;"> 1013 <div class="canvas-col"><canvas width="400" height="400" class="diagram"></canvas></div> 1014 <!-- Studio: full controls (sliders, buttons, full legend, hint). --> 1015 <div class="controls studio-only"> 1016 <label>Tiling resolution: <span data-role="tiles-label">6Ã6</span></label> 1017 <input type="range" data-role="tiles" min="2" max="10" step="1" value="6" /> 1018 <label style="margin-top:14px;">Walk speed: <span data-role="speed-label">1Ã</span></label> 1019 <input type="range" data-role="speed" min="0" max="5" step="1" value="2" /> 1020 <div style="display:flex; gap:6px; flex-wrap:wrap; margin-top:14px;"> 1021 <button data-role="new">â» New walk</button> 1022 <button data-role="play">⸠Pause</button> 1023 <button data-role="release" style="display:none;">â Release</button> 1024 </div> 1025 <div class="legend" style="margin-top:10px;"> 1026 <span><i style="background:var(--color-interface)"></i> Tile interfaces</span> 1027 <span><i 1028 style="width:12px; height:12px; background:rgba(42,95,184,0.14); border:2px solid var(--color-accent); border-radius:1px;"></i> 1029 Current co-edge subdomain</span> 1030 <span><i style="height:0; background:none; border-top:2px dashed var(--color-neumann);"></i> Neumann 1031 obstacles</span> 1032 </div> 1033 <p style="font-size:0.78rem; color:var(--color-text-muted); margin-top:14px;">The walker hops from 1034 interface to interface. At each hop it sits on a tile interface; the two tiles sharing that edge form 1035 a <b>co-edge subdomain</b> (blue). It then jumps to a uniformly random exit on that subdomain's 1036 boundary â rejection-sampling any point inside the obstacles â which lands on a new interface and 1037 re-centers the next co-edge subdomain, until it reaches the global boundary $Z_\tau$. (A naive, 1038 educational cartoon: uniform over exits, not the true harmonic measure.) Click an interface to launch 1039 walks from that spot; <b>Release</b> resumes random starts.</p> 1040 </div> 1041 </div> 1042 <!-- Public page: just the co-edge subdomain legend, centered below. --> 1043 <div class="controls blog-only" style="margin-top:14px;"> 1044 <div class="legend" style="align-items:center;"> 1045 <span><i 1046 style="width:12px; height:12px; background:rgba(42,95,184,0.14); border:2px solid var(--color-accent); border-radius:1px;"></i> 1047 Current co-edge subdomain</span> 1048 </div> 1049 </div> 1050 </div> 1051 <div class="figure-caption"><b>Walk on Co-edge Subdomains.</b> From a point $x$ on an interface, we take a 1052 step by walking on its <i>co-edge subdomain</i> â the union of the two tiles that share that interface. 1053 Each exit lands on a new interface, with a new subdomain for the next step, until the walk is 1054 absorbed on the global Dirichlet boundary at $Z_\tau$.</div> 1055 </div> 1056 </figure> 1057 1058 <p>With this, we can define a discrete Markov chain where the interfaces are the transient states, the global 1059 Dirichlet boundary gives the absorbing states, and the first-passage probabilities between interfaces define the 1060 transitions and are estimated using the strategy defined in <a class="ref" data-ref="poisson-kernel"></a> 1061 applied to the co-edge subdomains. 1062 The Markov chain formalism from <a class="ref" data-ref="markov"></a> lets us solve for interface values directly 1063 through a deterministic (sparse) linear solve â sidestepping the relatively long, high-variance random walks that 1064 would 1065 otherwise be required. 1066 </p> 1067 1068 <p>Once interface values are known, the boundary of every tile is fully determined, and we recover each tile's 1069 interior with a single matrixâvector product using precomputed solution operators for the interior, i.e., 1070 $\mathbf{H}_i$ as defined in <a class="ref" data-ref="poisson-kernel"></a>.
1071 These per-tile 1072 reconstructions are independent and thus trivially parallelizable.</p> 1073 1074 <p>The interactive figure below summarizes all steps of the pipeline. Feel free to slide through the various 1075 stages and play with the different parameters.</p> 1076 1077 <figure class="figure" id="pipeline-fig"> 1078 <div class="figure-body"> 1079 <div class="canvas-wrap"> 1080 <div class="canvas-row"> 1081 <div class="canvas-col"> 1082 <canvas width="380" height="380" class="diagram"></canvas> 1083 <div class="pipeline-stage-bar"> 1084 <label>Pipeline stage: <span class="brush-val" data-role="stage-label">0 â Domain</span></label> 1085 <input type="range" class="colormap-slider stage-slider" data-role="stage" min="0" max="4" step="1" 1086 value="0" list="pipeline-stage-ticks" /> 1087 <datalist id="pipeline-stage-ticks"> 1088 <option value="0" label="0"></option> 1089 <option value="1" label="1"></option> 1090 <option value="2" label="2"></option> 1091 <option value="3" label="3"></option> 1092 <option value="4" label="4"></option> 1093 </datalist> 1094 <div class="stage-tick-labels"> 1095 <span>0</span><span>1</span><span>2</span><span>3</span><span>4</span> 1096 </div> 1097 <p class="stage-desc" data-role="stage-desc"></p> 1098 </div> 1099 </div> 1100 <div class="controls"> 1101 <label>Tiles per axis T: <span data-role="tiles-label">4</span></label> 1102 <input type="range" data-role="tiles" min="1" max="8" step="1" value="4" /> 1103 1104 <label style="margin-top:10px;">Subtile resolution B: <span data-role="sub-label">4</span></label> 1105 <input type="range" data-role="sub" min="1" max="8" step="1" value="4" /> 1106 <div class="stat" data-role="res-label">N = T Ã B = 16</div> 1107 1108 <label style="margin-top:14px;">Neumann obstacle</label> 1109 <div class="presets" data-role="scene-presets"></div> 1110 1111 <label style="margin-top:14px;">Boundary presets</label> 1112 <div class="presets" data-role="presets"></div> 1113 1114 <div class="legend" style="margin-top:10px;"> 1115 <span><i style="background:var(--color-dirichlet)"></i> Dirichlet boundary</span> 1116 <span><i style="height:0; background:none; border-top:2px dashed var(--color-neumann);"></i> Neumann 1117 obstacle</span> 1118 <span><i style="background:var(--color-interface)"></i> Tile interfaces</span> 1119 </div> 1120 <p style="font-size:0.78rem; color:var(--color-text-muted); margin-top:14px;">Slide through the 1121 pipeline. 1122 </p> 1123 </div> 1124 </div> 1125 </div> 1126 <div class="figure-caption"><b>The Walk on Decomposed Subdomains (WoDS) pipeline.</b> 1127 (1) Partition $\Omega$ into tiles $\{\Omega_i\}$ separated by interfaces. 1128 (2) For each tile, estimate first-passage transition probabilities between its interfaces with short 1129 local Walk-on-Stars walks; assemble these into per-tile operators $\mathbf{H}_i$. 1130 (3) Stitch the corresponding probabilities into one global absorbing Markov chain over all interfaces and 1131 recover interface values via a single sparse solve 1132 $(\mathbf{I}-\mathbf{Q})\,\mathbf{u}_{\mathcal{T}} = \mathbf{R}\,\mathbf{g}$. 1133 (4) With every tile's boundary now known, reconstruct the interior in parallel by applying local interior 1134 solution operators, i.e., 1135 $\mathbf{H}_i$. 1136 </div> 1137 </div> 1138 </figure> 1139 1140 <!-- ===== SECTION 7 ===== --> 1141 <h2 data-label="benefits"><span class="secnum">7</span>Additional Benefits</h2> 1142 1143 1144 <p>One thing that I haven't mentioned is that in practice, we needn't compute solution operators for every tile of 1145 the domain. If a tile does not intersect geometry, its solution operator is the same everywhere and we can 1146 precompute it only once across all scenes<fn>It actually has a closed-form expression.</fn>. In other words, 1147 the number of solution operators that must be estimated grows with the geometric complexity and not domain size 1148 (i.e., 1149 area in 2D). This has another significant advantage: if only a subset of the scene changes â as is common in 1150 many 1151 design or simulation workflows â only the affected tiles need to be re-estimated.</p> 1152 1153 <figure class="figure" id="locality-fig"> 1154 <div class="figure-body"> 1155 <div class="canvas-wrap"> 1156 <div class="canvas-row"> 1157 <div class="canvas-col"> 1158 <canvas width="380" height="380" class="diagram"></canvas> 1159 </div> 1160 <div class="controls"> 1161 <label>Tiles per axis $T$: <span data-role="tiles-label">10Ã10</span></label> 1162 <input type="range" data-role="tiles" min="2" max="64" step="1" value="10" /> 1163 1164 <label style="margin-top:14px;">Scene</label>
1165 <div class="presets" data-role="scene-presets"></div> 1166 1167 <div class="studio-only" style="margin-top:14px;"> 1168 <button data-role="toggle-interfaces">Hide interfaces</button> 1169 </div> 1170 1171 <div class="stat" style="margin-top:14px;"> 1172 MC-estimated tiles: <span data-role="mc-label">â</span><br /> 1173 Precomputed (shared): <span data-role="pre-label">â</span> 1174 </div> 1175 1176 <div class="legend" style="margin-top:10px;"> 1177 <span><i style="background:rgba(255,196,0,0.7)"></i> Intersects geometry â needs Monte Carlo</span> 1178 <span><i style="background:rgba(94,176,122,0.7)"></i> Empty tile â shared operator</span> 1179 </div> 1180 <p style="font-size:0.78rem; color:var(--color-text-muted); margin-top:14px;">Pick a scene and adjust 1181 $T$. Only the highlighted tiles must have their solution operator estimated; all empty tiles share a 1182 single closed-form operator computed once.</p> 1183 </div> 1184 </div> 1185 </div> 1186 <div class="figure-caption"><b>Locality of solution-operator estimation.</b> Our method inherits the locality 1187 of 1188 grid-based approaches by requiring the estimation of solution operators only in regions near the geometry. 1189 As 1190 a result, the actual Monte Carlo cost scales with geometric complexity rather than domain size. 1191 </div> 1192 </div> 1193 </figure> 1194 1195 <p>Another beautiful thing is that our approach exposes an adjustable trade-off between the cost of stochastic 1196 Monte 1197 Carlo 1198 estimation of solution operators (which is highly parallelizable within each individual tile) and the cost of 1199 the 1200 deterministic global solve on interfaces. As can be seen in the interactive figure below, at fixed output 1201 resolution $N$, increasing the size $B$ of each tile requires more Monte Carlo effort per tile, but the 1202 deterministic global solve becomes cheaper because the interfaces it couples are fewer and farther apart. In 1203 other 1204 words, we can directly amortize the $O(N^2)$ degrees of freedom of the global solve to $O(N^2/B)$ degrees of 1205 freedom if we can afford more Monte Carlo<fn>You may argue that what matters most is not the number of degrees 1206 of 1207 freedom but the sparsity or conditioning of the matrix. It turns out that our matrices are sparse because they 1208 only couple neighboring interfaces. In other words, the 1209 global connectivity of the system closely resembles 1210 finite-difference systems, except that each interface-to-interface connection may carry multiple "transition 1211 edges" proportional to the per-tile discretization parameter $B$.</fn> 1212 <fn>There's another subtle point, which is that our matrices aren't symmetric as we describe and discuss in 1213 Section 8 of the 1214 paper.</fn>. This is a significant result, which I 1215 hope 1216 we can build upon in the 1217 future because Monte Carlo is embarrassingly parallel and can run on accelerated hardware (e.g., GPUs)<fn>What 1218 I'm 1219 not saying here is that there is not only a cost in samples but also in memory! Larger tiles mean more memory 1220 usage, 1221 and actually quite a lot! But as we discuss in the paper, I believe our community has developed many 1222 adjacent tools to mitigate this issue, and I hope I or others will explore them in the future. 1223 </fn>. 1224 </p> 1225 1226 <figure class="figure" id="tradeoff-fig"> 1227 <div class="figure-body"> 1228 <div class="canvas-wrap"> 1229 <div class="canvas-row" style="flex-wrap:nowrap; justify-content:center;"> 1230 <div class="canvas-col" style="width:320px;"> 1231 <canvas width="320" height="320" class="diagram" data-role="domain"></canvas> 1232 <p style="font-size:0.78rem; color:var(--color-text-muted); margin-top:8px; width:320px;"> 1233 Interface values at fixed $N$. Tile interfaces are the only unknowns of the global 1234 solve. 1235 </p> 1236 </div> 1237 <div class="canvas-col" style="width:320px; position:relative;"> 1238 <canvas width="320" height="320" class="diagram" data-role="plot"></canvas>
1239 <span data-role="x-axis-label" 1240 style="position:absolute; left:181px; top:303px; transform:translate(-50%,0); font-size:12px; color:var(--color-text); white-space:nowrap;">Subtile 1241 resolution $B$</span> 1242 <span data-role="y-axis-label" 1243 style="position:absolute; left:8px; top:149px; transform:translate(-50%,-50%) rotate(-90deg); font-size:12px; color:var(--color-text); white-space:nowrap;">System 1244 size $N^2 / B$</span> 1245 <p style="font-size:0.78rem; color:var(--color-text-muted); margin-top:8px; width:320px;"> 1246 System size of the deterministic global solve as a function of subtile resolution $B$. 1247 </p> 1248 </div> 1249 </div> 1250 <div class="controls" style="margin-top:14px; width:660px; max-width:100%;"> 1251 <label style="font-variant-numeric: tabular-nums;">Subtile resolution 1252 $B$ = <span data-role="b-label" style="display:inline-block; min-width:2ch; text-align:right;">4</span> 1253 · 1254 Tiles per axis $T = N/B$ = <span data-role="t-label" 1255 style="display:inline-block; min-width:2ch; text-align:right;">16</span> · 1256 Global-solve DoFs $\approx N^2/B$ = <span data-role="dofs-label" 1257 style="display:inline-block; min-width:4ch; text-align:right;">1024</span></label> 1258 <input type="range" class="colormap-slider plain-slider" data-role="b-slider" min="0" max="6" step="1" 1259 value="2" /> 1260 <div class="slider-ticks" aria-hidden="true"> 1261 <span></span><span></span><span></span><span></span><span></span><span></span><span></span> 1262 </div> 1263 <div style="display:flex; justify-content:space-between; font-size:0.72rem;"> 1264 <span style="color:#6b7d54;">Grid-like ($T{=}N$, $B{=}1$)</span> 1265 <span style="color:#b07a3a;">Pure operator ($T{=}1$, $B{=}N$)</span> 1266 </div> 1267 </div> 1268 </div> 1269 <div class="figure-caption"><b>Trade-off between local Monte Carlo and global solve.</b> 1270 At fixed output resolution $N$, sliding $B$ from 1 to $N$ continuously interpolates between a grid-like 1271 regime â where the global solve carries all $O(N^2)$ degrees of freedom â and a pure solution-operator 1272 regime â where the global solve is trivial ($O(1)$) but every tile requires a fully tabulated interior 1273 operator. 1274 </div> 1275 </div> 1276 </figure> 1277 1278 1279 <!-- ===== SECTION 8 ===== --> 1280 <h2 data-label="reflect"><span class="secnum">8</span>Reflections and Future Work</h2> 1281 1282 <p>I am very excited by our work and the doors it opens for future research. One thing people often bring up is 1283 the 1284 analogy to <a href="https://en.wikipedia.org/wiki/Radiosity_(computer_graphics)" target="_blank" 1285 rel="noopener noreferrer">radiosity</a> and Monte Carlo path tracing in rendering<fn>It turns 1286 out our approach has even more connections to previous rendering works and in particular volume rendering. See 1287 for example these works<cite data-key="zhao2013modular,blumer2016reduced"></cite>.</fn>. 1288 There was once a time when people used finite elements for rendering and were actually very reluctant to adopt 1289 Monte Carlo methods because of their noise. 1290 </p> 1291 1292 <p>
1292In this context, Peter Shirley famously said in an <a 1293 href="https://www.realtimerendering.com/resources/RTNews/html/rtnv10n2.html#art6" target="_blank" 1294 rel="noopener noreferrer">email thread</a> in 1997<cite data-key="shirley1997whatswrong"></cite>:</p> 1295 1296 <blockquote class="email-quote"> 1297 <p>In summary, pure MCPT has only two advantages — it is so dumb that it doesn't get hit by big scenes, 1298 and 1299 it is easy to implement. [â¦]</p> 1300 <p>I think the solution is <span class="quote-highlight">hybrid methods</span> — add bias! (This is 1301 blasphemy 1302 in MC circles :^) ). I do want to keep the good parts of MC methods — they are damned robust and are 1303 possible to implement correctly — my MC code does not bomb on weird untweaked inputs — tell me 1304 with 1305 a straight face that is true of most non-MC implementations. However, you are right that the results are too 1306 noisy!!! We can keep these benefits and reduce noise <span class="quote-highlight">if we add bias the right 1307 way 1308 (not that I know what that right way is)</span>.</p> 1309 </blockquote> 1310 1311 <p>One thing I highlighted here is that Peter Shirley was right: to be adopted, Monte Carlo methods needed "hybrid 1312 methods". And the solution was <a href="https://blogs.nvidia.com/blog/what-is-denoising/" target="_blank" 1313 rel="noopener noreferrer">denoising</a>! Without denoising, existing visual effects and animated 1314 films would probably not be path traced! People even <a href="https://www.oscars.org/sci-tech/ceremonies/2025" 1315 target="_blank" rel="noopener noreferrer">won an Academy Award</a> for that.</p> 1316 1317 <p>So, to be competitive with classical numerical methods, do Monte Carlo PDE solvers also need 1318 hybrid approaches? Probably yes! But is denoising the right answer? Maybe not.</p> 1319 1320 <p>
1320In our paper, we take a different route and instead try to reconcile grid-free Monte Carlo methods with 1321 grid-based solvers. One thing that I learned with this project is that grid-based methods just work so well: if 1322 you want 1323 convergence, you simply cannot beat them! And it's not a surprise, ask anyone doing numerical analysis and 1324 they 1325 will tell you something like:</p> 1326 1327 <blockquote class="email-quote"> 1328 <p>Monte Carlo is an extremely bad method; it should be used only when all alternative methods are 1329 worse.<cite data-key="sokal1997monte"></cite></p> 1330 </blockquote> 1331 1332 <p>My take is: lean on the well-known benefits of grid-based methods as much as possible (e.g., strong and fast 1333 convergence guarantees), and bring in grid-free Monte Carlo where grid-based methods suffer most (e.g., complex 1334 geometry).</p> 1335 1336 <p>That said, there is still a long way to go: generalizing to other PDEs or boundary conditions, better 1337 discretization schemes, dedicated global solvers, improved Monte Carlo estimators, etc. Rest assured that we are 1338 working 1339 on 1340 that...</p> 1341 1342 <!-- ===== References ===== --> 1343 <h2>References</h2> 1344 <div id="bibliography"></div> 1345 1346 </article> 1347 1348 <footer> 1349 <a href="https://adg.csail.mit.edu/" class="footer-logo" target="_blank" rel="noopener noreferrer" 1350 aria-label="Algorithmic Design Group"><img src="assets/logo_adg.png" alt="Algorithmic Design Group" width="1262" height="728" loading="lazy" decoding="async" /></a> 1351 <span>Webpage by <a href="https://clementjambon.github.io/" target="_blank" rel="noopener noreferrer">Clément 1352 Jambon</a> and <a href="https://claude.ai/" target="_blank" rel="noopener noreferrer">Claude</a>.</span> 1353 </footer> 1354 1355 <!-- ===== Scripts ===== --> 1356
1356<script src="js/config.js"></script>
1356 1357
1357<script src="js/util.js"></script>
1357 1358
1358<script src="js/laplace.js"></script>
1358 1359
1359<script src="js/solver.js"></script>
1359 1360
1360<script src="js/scenes.js"></script>
1360 1361
1361<script src="js/interactive_dirichlet.js"></script>
1361 1362
1362<script src="js/interactive_neumann.js"></script>
1362 1363
1363<script src="js/interactive_brownian.js"></script>
1363 1364
1364<script src="js/interactive_denoise.js"></script>
1364 1365
1365<script src="js/interactive_walks.js"></script>
1365 1366
1366<script src="js/interactive_tiles.js"></script>
1366 1367
1367<script src="js/interactive_poisson_kernel.js"></script>
1367 1368
1368<script src="js/interactive_poisson_solution.js"></script>
1368 1369
1369<script src="js/interactive_subdomain_kernel.js"></script>
1369 1370
1370<script src="js/interactive_binned_operator.js"></script>
1370 1371
1371<script src="js/interactive_coedge.js"></script>
1371 1372
1372<script src="js/interactive_markov.js"></script>
1372 1373
1373<script src="js/interactive_pipeline.js"></script>
1373 1374
1374<script src="js/interactive_tradeoff.js"></script>
1374 1375
1375<script src="js/interactive_locality.js"></script>
1375 1376
1376<script src="js/auto_markov.js"></script>
1376 1377
1377<script src="js/references.js"></script>
1377 1378
1378<script src="js/labels.js"></script>
1378 1379
1379<script src="js/capture.js"></script>
1379 1380
1380<script> 1381 document.addEventListener('DOMContentLoaded', () => { 1382 WoDS.interactiveDirichlet(document.getElementById('i0')); 1383 WoDS.interactiveNeumann(document.getElementById('i0b')); 1384 WoDS.interactiveBrownian(document.getElementById('ibrown')); 1385 WoDS.interactiveDenoise(document.getElementById('idenoise')); 1386 WoDS.interactiveWalks(document.getElementById('i1')); 1387 WoDS.interactiveTiles(document.getElementById('i2')); 1388 WoDS.interactivePoissonKernel(document.getElementById('i3')); 1389 WoDS.interactivePoissonSolution(document.getElementById('i3b')); 1390 // subkernel + binop are studio-only (presentation) figures: run them 1391 // on studio.html but not on the public page. build-studio.js still 1392 // matches these calls and re-emits them (unguarded) into studio.html. 1393 if (WoDS.inStudio) WoDS.interactiveSubdomainKernel(document.getElementById('subkernel')); 1394 if (WoDS.inStudio) WoDS.interactiveBinnedOperator(document.getElementById('binop')); 1395 WoDS.interactiveCoedge(document.getElementById('coedge-fig')); 1396 WoDS.autoMarkov(document.getElementById('amc-anim')); 1397 WoDS.interactiveMarkov(document.getElementById('imc')); 1398 WoDS.interactivePipeline(document.getElementById('pipeline-fig')); 1399 WoDS.interactiveTradeoff(document.getElementById('tradeoff-fig')); 1400 WoDS.interactiveLocality(document.getElementById('locality-fig')); 1401 1402 // Scale display equations to fit their column. KaTeX renders 1403 // equations at their natural width; on narrow viewports the 1404 // wider ones (multi-line aligned blocks, long stencil formulas) 1405 // overflow. We measure each .katex-display and apply a 1406 // proportional transform so it fits without horizontal scroll. 1407 function fitEquations() { 1408 document.querySelectorAll('.katex-display').forEach(el => { 1409 const inner = el.querySelector('.katex'); 1410 if (!inner) return; 1411 // Reset before measuring. 1412 inner.style.fontSize = ''; 1413 const containerW = el.clientWidth; 1414 if (!containerW) return; 1415 // Measure the widest of scrollWidth and the rendered 1416 // .katex-html bounding rect â for aligned/mtable layouts 1417 // scrollWidth alone can under-report (it only catches 1418 // right-edge overflow, but centered content overflows 1419 // both sides). 1420 let innerW = inner.scrollWidth; 1421 const html = inner.querySelector('.katex-html'); 1422 if (html) innerW = Math.max(innerW, html.scrollWidth, html.getBoundingClientRect().width); 1423 if (innerW > containerW) { 1424 // KaTeX uses em units throughout, so shrinking font-size 1425 // proportionally shrinks the layout (not just visually, 1426 // unlike `transform: scale`) and the equation re-centers 1427 // naturally via the parent's text-align: center. We 1428 // scale relative to the *current* computed font-size 1429 // (KaTeX's own stylesheet sets .katex to 1.21em) â using 1430 // a percentage would clobber KaTeX's baseline and make 1431 // the equation much smaller than needed. 1432 const s = containerW / innerW; 1433 const cur = parseFloat(getComputedStyle(inner).fontSize); 1434 inner.style.fontSize = (cur * s) + 'px'; 1435 } 1436 }); 1437 } 1438 // Run after KaTeX has rendered (auto-render runs on its script 1439 // load, which may finish after DOMContentLoaded), and again 1440 // once KaTeX's web fonts finish loading â math glyphs swap 1441 // from fallback to the real font and widths change. 1442 window.addEventListener('load', fitEquations); 1443 if (document.fonts && document.fonts.ready) { 1444 document.fonts.ready.then(fitEquations); 1445 } 1446 let _eqResizeT; 1447 window.addEventListener('resize', () => { 1448 clearTimeout(_eqResizeT); 1449 _eqResizeT = setTimeout(fitEquations, 100); 1450 }); 1451 1452 // Bibtex copy 1453 const btn = document.querySelector('[data-role="copy-bib"]'); 1454 if (btn) btn.addEventListener('click', () => { 1455 const code = document.querySelector('#bibtex code').innerText; 1456 navigator.clipboard.writeText(code).then(() => { 1457 btn.textContent = 'copied'; 1458 setTimeout(() => btn.textContent = 'copy', 1500); 1459 }); 1460 }); 1461 }); 1462 </script>
1462 1463 1464</body> 1465 1466</html>
Line numbers count LF bytes from the start of the resource, as the search results do. Vendor segments are library code the classifier recognised; they are stored but not indexed. Bytes are shown as Latin1 characters, one per byte.