PageSourceSearch

https://clementjambon.github.io/wods/index.html

html clementjambon.github.io collected 2026-10-03 08:48:16 UTC 88,400 bytes, 1,466 lines download raw bytes

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 &ldquo;faked&rdquo; 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$ &amp; 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 &amp; 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) &amp;\;=\; \sum_{j \in \mathcal{S}} P_{ij}\, u(j) &amp;&amp; \text{on } \mathcal{T}, \\
817      u(a) &amp;\;=\; g(a) &amp;&amp; \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} &amp; \mathbf{R} \\ \mathbf{0} &amp; \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              &nbsp;·&nbsp;
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> &nbsp;·&nbsp;
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 &mdash; 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> &mdash; add bias! (This is
1301        blasphemy
1302        in MC circles :^) ). I do want to keep the good parts of MC methods &mdash; they are damned robust and are
1303        possible to implement correctly &mdash; my MC code does not bomb on weird untweaked inputs &mdash; 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.