Walk on Spheres: Monte Carlo PDE Solver

• Walk on Spheres is a grid-free Monte Carlo method defined by random jumps on ball boundaries to solve Dirichlet and Poisson problems without volumetric discretization.
• It leverages stochastic Brownian motion exit distributions to enable adaptive step sizes and variance reduction techniques like caching and RQMC for complex geometries.
• Modern extensions address mixed boundary conditions, fractional operators, and surface PDEs by integrating neural and quasi-Monte Carlo methods for enhanced performance.

Walk on Spheres (WoS) is a grid-free Monte Carlo method for elliptic boundary value problems. In its classical form, it solves Dirichlet problems for the Laplace or Poisson family by replacing volumetric discretization and global linear solves with a recursive sequence of random jumps: from the current point, one moves to a random point on the boundary of the largest ball centered at that point and still contained in the domain. The method is grounded in stochastic representations of PDE solutions by Brownian motion, and modern variants extend the same geometric-probabilistic idea to screened Poisson equations, mixed Neumann–Dirichlet problems, fractional operators driven by $\alpha$-stable Lévy processes, surface and manifold PDEs, and learning-based solvers (Kyprianou et al., 2016, Sawhney et al., 2023, Nam et al., 2024).

1. Probabilistic foundation

The classical setting is the Dirichlet problem

$$\Delta u(x)=0,\quad x\in D,\qquad u(x)=g(x),\quad x\in \partial D.$$

Its stochastic representation is

$$u(x)=\mathbb{E}_x[g(W_{\tau_D})],\qquad \tau_D=\inf\{t>0: W_t\notin D\},$$

with $W$ Brownian motion. This is the basic reason WoS works: instead of discretizing Brownian paths in time, it exploits the fact that Brownian motion started at the center of a ball exits that ball uniformly on the sphere boundary (Kyprianou et al., 2016).

For Poisson-type problems, the same idea is combined with a local Green representation on balls. Recent formulations write the solution on a ball as a boundary term plus a volume term involving the Green’s function, so a WoS trajectory carries both a boundary contribution and an accumulated source contribution. In neural formulations, this appears as a recursive decomposition

$$u(\xi)=\mathbb{E}\!\left[u(X_{\tau_0}^\xi)-F_{\tau_0}^\xi\mid \xi\right],$$

which is then iterated over a sequence of nested balls (Nam et al., 2024). Surface and embedded-manifold generalizations retain the same ball-based Green representation, but interpret off-surface samples through closest-point projection (Sugimoto et al., 2024, Hui et al., 20 Jun 2026).

The central conceptual distinction is therefore between pathwise diffusion simulation and sphere-exit simulation. WoS belongs to the latter class: it preserves the exact exit distribution from each local ball while avoiding fine time stepping.

2. Classical algorithm and numerical behavior

The basic WoS update is geometrically simple. Starting from $x=\rho_0\in D$, one repeatedly finds the largest sphere centered at the current point and contained in the domain, samples the next point uniformly on that sphere, and repeats. In the Brownian Dirichlet case, if $r_n$ is the radius of the largest inscribed sphere centered at $\rho_{n-1}$ and

$S_n=\{y\in\mathbb{R}^d: |y-\rho_{n-1}|=r_n\},$

then $\rho_n$ is sampled uniformly on $$\Delta u(x)=0,\quad x\in D,\qquad u(x)=g(x),\quad x\in \partial D.$$0 (Kyprianou et al., 2016).

Because Brownian motion approaches the boundary continuously, practical implementations stop not at the exact boundary but in an $$\Delta u(x)=0,\quad x\in D,\qquad u(x)=g(x),\quad x\in \partial D.$$1-skin,

$$\Delta u(x)=0,\quad x\in D,\qquad u(x)=g(x),\quad x\in \partial D.$$2

and use the boundary value at a nearest boundary point. Monte Carlo consistency is then obtained by averaging independent trajectories, with the standard strong-law estimator

$$\Delta u(x)=0,\quad x\in D,\qquad u(x)=g(x),\quad x\in \partial D.$$3

(Kyprianou et al., 2016).

A recurring practical advantage is that the steps are adaptive and typically large. Modern treatments emphasize progressive refinement, trivial parallelization, pointwise evaluation only where needed, and sublinear scaling with geometric detail (Sawhney et al., 2023). The same literature also records the standard tradeoff: the $$\Delta u(x)=0,\quad x\in D,\qquad u(x)=g(x),\quad x\in \partial D.$$4-shell prevents degenerate late-stage walks but introduces stopping bias. Reported asymptotics include an $$\Delta u(x)=0,\quad x\in D,\qquad u(x)=g(x),\quad x\in \partial D.$$5 bias for basic Dirichlet estimators and an $$\Delta u(x)=0,\quad x\in D,\qquad u(x)=g(x),\quad x\in \partial D.$$6 bias for certain pointwise WoS/WoSt estimators used in cache-based evaluation (Czekanski et al., 2024, Miller et al., 2023).

The expected step count is logarithmic in the boundary tolerance in several classical settings. For convex domains, one bound is

$$\Delta u(x)=0,\quad x\in D,\qquad u(x)=g(x),\quad x\in \partial D.$$7

and related logarithmic complexity bounds reappear in moving-sphere variants for Brownian and Ornstein–Uhlenbeck exit problems (Kyprianou et al., 2016, Deaconu et al., 2014, Herrmann et al., 2019).

3. Boundary conditions, operators, and process generalizations

The most direct classical use of WoS is for pure Dirichlet problems. Modern work repeatedly describes this case as the regime where WoS is simplest and most natural (Sawhney et al., 2023). Extending the method beyond that setting has produced several distinct lines of development.

For mixed Neumann–Dirichlet problems, "Walk on Stars" replaces the spherical subdomains of WoS with star-shaped regions that may contain visible Neumann boundary pieces. The PDE is written as

$$\Delta u(x)=0,\quad x\in D,\qquad u(x)=g(x),\quad x\in \partial D.$$8

with $$\Delta u(x)=0,\quad x\in D,\qquad u(x)=g(x),\quad x\in \partial D.$$9. The key insight is that reflecting Brownian motion can be simulated by allowing the local domain to include visible Neumann boundary while preserving an exact Monte Carlo estimator for the Neumann flux term. The radius is controlled by both the distance to Dirichlet boundary and the distance to the closest point on the visibility silhouette of the Neumann boundary, and the implementation uses a spatialized normal cone hierarchy for the new silhouette query (Sawhney et al., 2023).

A different extension concerns exit times rather than only exit locations. "Walk on Moving Spheres" and related WoMS/WOMS constructions replace the fixed sphere by a moving sphere whose radius is a carefully chosen function of time. In Brownian and Bessel-process settings, the moving boundary is engineered so that the hitting-time law is explicit, allowing efficient joint simulation of exit location and exit time without time discretization; Ornstein–Uhlenbeck variants use generalized spheroids obtained from the Brownian time-change representation (Deaconu et al., 2011, Deaconu et al., 2014, Herrmann et al., 2019).

For constant-potential Schrödinger equations, including the Yukawa and Helmholtz cases, the Duffin correspondence reformulates

$$u(x)=\mathbb{E}_x[g(W_{\tau_D})],\qquad \tau_D=\inf\{t>0: W_t\notin D\},$$0

as a Laplace equation in one higher dimension, so that classical WoS can be applied to the transformed problem rather than to a modified weighted or killing scheme (Yang et al., 2015).

For nonlocal operators, the classical Brownian walk is replaced by an isotropic $$u(x)=\mathbb{E}_x[g(W_{\tau_D})],\qquad \tau_D=\inf\{t>0: W_t\notin D\},$$1-stable Lévy process. In the homogeneous fractional Dirichlet problem,

$$u(x)=\mathbb{E}_x[g(W_{\tau_D})],\qquad \tau_D=\inf\{t>0: W_t\notin D\},$$2

the process exits balls by jumps rather than by continuous boundary hits. That change has several structural consequences recorded in the literature: the algorithm terminates almost surely without $$u(x)=\mathbb{E}_x[g(W_{\tau_D})],\qquad \tau_D=\inf\{t>0: W_t\notin D\},$$3-truncation in the homogeneous convex case, disconnected domains can be handled, and the estimator is unbiased (Kyprianou et al., 2016). High-dimensional fractional Poisson variants add an interior source sample from a Green-function-based density, modified sampling schemes in hyperspherical coordinates, and neural amortization for dimensions up to $$u(x)=\mathbb{E}_x[g(W_{\tau_D})],\qquad \tau_D=\inf\{t>0: W_t\notin D\},$$4 (Jiao et al., 2022, Guo et al., 30 Jan 2026).

Partially reflecting random walk on spheres methods introduce finite-difference-based replacement rules near Robin boundaries and transmission interfaces. In electrical impedance tomography this yields a solver that uses standard WoS in constant-coefficient subdomains, but near the boundary replaces pure absorption by a stochastic mixture of reflection, absorption, and transmission (Maire et al., 2015).

4. Sample reuse, variance reduction, and learned surrogates

A standard limitation of classical WoS is that each query point is estimated independently. Several recent methods address this by reusing information across space.

Boundary Value Caching samples boundary points and source points, estimates unknown boundary values or normal derivatives there using WoS or WoSt, and then evaluates the boundary integral representation at interior points through Monte Carlo sums. The cache itself introduces no additional statistical bias: by the formulation in that work, unbiased cached boundary estimates imply an unbiased interior estimator via linearity of expectation. The same cache can also be reused for gradient evaluation (Miller et al., 2023).

A different reuse strategy exploits Brownian path continuity. For Laplace’s equation with Dirichlet data, a walk started at one point can be reweighted by a Poisson-kernel factor to produce an unbiased estimator for nearby points inside the same first sphere. This leads to “talking to neighbors” estimators, variance bounds controlled by the relative distance to the sphere radius, and a deterministic fixed-size cache construction of size

$$u(x)=\mathbb{E}_x[g(W_{\tau_D})],\qquad \tau_D=\inf\{t>0: W_t\notin D\},$$5

for the interior region $$u(x)=\mathbb{E}_x[g(W_{\tau_D})],\qquad \tau_D=\inf\{t>0: W_t\notin D\},$$6 (Czekanski et al., 2024).

Randomized quasi-Monte Carlo has also been coupled to WoS. One RQMC formulation reports variance reduction factors ranging from $$u(x)=\mathbb{E}_x[g(W_{\tau_D})],\qquad \tau_D=\inf\{t>0: W_t\notin D\},$$7 to $$u(x)=\mathbb{E}_x[g(W_{\tau_D})],\qquad \tau_D=\inf\{t>0: W_t\notin D\},$$8 at $$u(x)=\mathbb{E}_x[g(W_{\tau_D})],\qquad \tau_D=\inf\{t>0: W_t\notin D\},$$9 and median variance rates slightly better than $W$0 over several examples (Ho et al., 8 May 2026). An Array-RQMC formulation, based on stepwise low-dimensional randomized point sets together with Hilbert sorting of the current walker cloud, reports variance reduction factors ranging from $W$1-fold to $W$2-fold at $W$3 trajectories and empirical rates between $W$4 and $W$5 (Ho et al., 13 May 2026).

Differential WoS differentiates the PDE solution map with respect to problem parameters rather than differentiating a particular Monte Carlo trace. For screened Poisson equations, the sensitivity $W$6 satisfies the same interior operator with a boundary condition involving the normal derivative of the primal solution, so derivatives with respect to many parameters can be estimated jointly at arbitrary points without a global solve (Miller et al., 2024).

Learning-based methods use WoS as supervision rather than only as a direct solver. Neural Walk-on-Spheres trains a neural network with the loss

$W$7

thereby amortizing pointwise Monte Carlo estimates across the whole domain (Nam et al., 2024). "Walk-on-Spheres Neural Operator" uses noisy WoS estimates as weak supervision for operator learning, avoiding higher-order derivatives in the loss and removing the need for expensive precomputed datasets; reported gains include up to $W$8 improvement in $W$9-error, up to $$u(\xi)=\mathbb{E}\!\left[u(X_{\tau_0}^\xi)-F_{\tau_0}^\xi\mid \xi\right],$$0 improvement in training speed, and up to $$u(\xi)=\mathbb{E}\!\left[u(X_{\tau_0}^\xi)-F_{\tau_0}^\xi\mid \xi\right],$$1 reduction in GPU memory consumption relative to standard physics-informed training schemes (Viswanath et al., 1 Mar 2026). Fractional neural variants similarly combine simplified fractional WoS with neural surrogates and buffered supervision, with tests on irregular domains and dimensions up to $$u(\xi)=\mathbb{E}\!\left[u(X_{\tau_0}^\xi)-F_{\tau_0}^\xi\mid \xi\right],$$2 (Guo et al., 30 Jan 2026).

5. Surface and manifold formulations

Classical WoS is volumetric. Projected Walk on Spheres extends the method to PDEs posed on surfaces by combining WoS with closest-point extension. The closest-point map is

$$u(\xi)=\mathbb{E}\!\left[u(X_{\tau_0}^\xi)-F_{\tau_0}^\xi\mid \xi\right],$$3

and the closest-point extension is

$$u(\xi)=\mathbb{E}\!\left[u(X_{\tau_0}^\xi)-F_{\tau_0}^\xi\mid \xi\right],$$4

The basic modification is to walk inside a volumetric neighborhood around the surface, sample a point on the boundary of a valid ball, and project that point back to the surface after each step. The admissible neighborhood is controlled by local feature size and by the distance to Dirichlet boundaries extended along surface normals; the method also includes a mean value filtering strategy for efficiency on meshes and point clouds (Sugimoto et al., 2024).

For screened Poisson equations on embedded manifolds, projected WoS has been developed in a tubular neighborhood $$u(\xi)=\mathbb{E}\!\left[u(X_{\tau_0}^\xi)-F_{\tau_0}^\xi\mid \xi\right],$$5 with closest-point projection

$$u(\xi)=\mathbb{E}\!\left[u(X_{\tau_0}^\xi)-F_{\tau_0}^\xi\mid \xi\right],$$6

The ambient Euclidean equation is corrected by a compensation term

$$u(\xi)=\mathbb{E}\!\left[u(X_{\tau_0}^\xi)-F_{\tau_0}^\xi\mid \xi\right],$$7

which measures the discrepancy between the ambient Laplacian of the closest-point extension and the Laplace–Beltrami operator away from the manifold. The resulting stochastic recursion uses projected boundary samples, local Green representations in Euclidean balls, and an adaptive radius

$$u(\xi)=\mathbb{E}\!\left[u(X_{\tau_0}^\xi)-F_{\tau_0}^\xi\mid \xi\right],$$8

Under approximate projection and compensation errors, mean-square error estimates take the form

$$u(\xi)=\mathbb{E}\!\left[u(X_{\tau_0}^\xi)-F_{\tau_0}^\xi\mid \xi\right],$$9

in the boundary case, and analogous estimates hold on closed manifolds (Hui et al., 20 Jun 2026).

These projected formulations preserve the meshfree, pointwise, and highly parallelizable character of WoS while replacing direct Laplace–Beltrami discretization by ambient-space sampling plus geometric projection.

6. Applications, competing methods, and limitations

The application range recorded in recent work is broad. Mixed-boundary WoSt has been demonstrated on a toaster heating bread, a lizard with spatially varying heat absorption, and a lung model for gas diffusion, all on large, detailed geometry without preprocessing into a volume mesh (Sawhney et al., 2023). Surface PWoS has been applied to diffusion curves, geodesic distance computation, and wave propagation animation, including surfaces of mixed codimension (Sugimoto et al., 2024). In robotics, WoS has been used for PDE-based path planning in configuration spaces; the reported experiments indicate trivial parallelization, empirical $x=\rho_0\in D$0 convergence in the number of walks, and validation on the RR platform (Muchacho et al., 2024). In electrical impedance tomography, partially reflecting random walk on spheres yields an embarrassingly parallel estimator for the voltage-to-current map (Maire et al., 2015).

The method is also frequently compared with mesh-based and boundary-based alternatives. Relative to FEM and BEM, the recurring advantages claimed for WoS-family solvers are the absence of volumetric meshing, pointwise evaluation only where needed, robustness to complex or imperfect geometry, and natural support for progressive Monte Carlo refinement (Sawhney et al., 2023, Miller et al., 2023). At the same time, the literature is explicit about classical WoS limitations: the $x=\rho_0\in D$1-shell causes intrinsic near-boundary bias, exact evaluation on the boundary is not natural, and Neumann- or Robin-dominated problems are substantially more delicate. "Walk on Boundary" is presented as an alternative that samples actual boundary intersections, naturally supports Dirichlet, Neumann, Robin, and mixed conditions for both interior and exterior domains, and does not suffer from the intrinsic bias of WoS near the boundary (Sugimoto et al., 2023).

Within the WoS lineage itself, the pattern is similar. Pure Dirichlet problems remain the setting in which classical WoS is simplest and often strongest, whereas Neumann-dominated problems can induce long walks, motivating star-shaped local domains, minimum-radius heuristics, or cached evaluation schemes (Sawhney et al., 2023). Fractional variants replace continuous sphere hitting by jump exits and thereby change the stopping theory fundamentally (Kyprianou et al., 2016). Learning-augmented versions trade unbiased pointwise sampling for amortized inference and reduced memory or runtime at evaluation time (Nam et al., 2024, Viswanath et al., 1 Mar 2026).

Taken together, these developments define Walk on Spheres not as a single algorithm but as a geometric Monte Carlo paradigm: local boundary-value problems are solved on adaptively chosen balls, and global PDE information is accumulated through recursive stochastic transport. The paradigm remains anchored in Brownian or Lévy exit representations, but its modern forms encompass reflecting and partially reflecting walks, moving spheres, star-shaped domains, surface projection, quasi-Monte Carlo sampling, differentiable estimators, and weakly supervised neural operators.