Walk-on-spheres method

In mathematics, the walk-on-spheres method (WoS) is a numerical probabilistic algorithm, or Monte-Carlo method, used mainly in order to approximate the solutions of some specific boundary value problem for partial differential equations (PDEs).^[1][2] The WoS method was first introduced by Mervin E. Muller in 1956 to solve Laplace's equation,^[1] and was since then generalized to other problems.

It relies on probabilistic interpretations of PDEs, and simulates paths of Brownian motion (or for some more general variants, diffusion processes), by sampling only the exit-points out of successive spheres, rather than simulating in detail the path of the process. This often makes it less costly than "grid-based" algorithms, and it is today one of the most widely used "grid-free" algorithms for generating Brownian paths.

Let ${\displaystyle \Omega }$ be a bounded domain in ${\displaystyle \mathbb {R} ^{d}}$ with a sufficiently regular boundary ${\displaystyle \Gamma }$, let h be a function on ${\displaystyle \Gamma }$, and let ${\displaystyle x}$ be a point inside ${\displaystyle \Omega }$.

Consider the Dirichlet problem:

${\displaystyle {\begin{cases}\Delta u(x)=0&{\mbox{if }}x\in \Omega \\u(x)=h(x)&{\mbox{if }}x\in \Gamma .\end{cases}}}$

It can be easily shown^[a] that when the solution ${\displaystyle u}$ exists, for ${\displaystyle x\in \Omega }$:

${\displaystyle u(x)=\mathbb {E} _{x}[h(W_{\tau })]}$

where W is a d-dimensional Wiener process, the expected value is taken conditionally on {W₀ = x}, and τ is the first-exit time out of Ω.

To compute a solution using this formula, we only have to simulate the first exit point of independent Brownian paths since with the law of large numbers:

${\displaystyle \mathbb {E} _{x}[h(W_{\tau })]\sim {\frac {1}{n}}\sum _{i=1}^{n}h(W_{\tau }^{i})}$

The WoS method provides an efficient way of sampling the first exit point of a Brownian motion from the domain, by remarking that for any (d − 1)-sphere centred on x, the first point of exit of W out of the sphere has a uniform distribution over its surface. Thus, it starts with x⁽⁰⁾ equal to x, and draws the largest sphere ${\displaystyle {\mathcal {S}}_{0}}$ centered on x⁽⁰⁾ and contained inside the domain. The first point of exit x⁽¹⁾ from ${\displaystyle {\mathcal {S}}_{0}}$ is uniformly distributed on its surface. By repeating this step inductively, the WoS provides a sequence (x⁽ⁿ⁾) of positions of the Brownian motion.

According to intuition, the process will converge to the first exit point of the domain. However, this algorithm takes almost surely an infinite number of steps to end. For computational implementation, the process is usually stopped when it gets sufficiently close to the border, and returns the projection of the process on the border. This procedure is sometimes called introducing an ${\displaystyle \varepsilon }$-shell, or ${\displaystyle \varepsilon }$-layer.^[4]

Choose ${\displaystyle \varepsilon >0}$. Using the same notations as above, the Walk-on-spheres algorithm is described as follows:

1. Initialize : ${\displaystyle x^{(0)}=x}$
2. While ${\displaystyle d(x^{(n)},\Gamma )>\varepsilon }$:
   1. Set ${\displaystyle r_{n}=d(x^{(n)},\Gamma )}$.
   2. Sample ${\displaystyle \gamma _{n}}$ a vector uniformly distributed over the unit sphere, independently from the preceding ones.
   3. Set ${\displaystyle x^{(n+1)}:=x^{(n)}+r_{n}\gamma _{n}}$
3. When ${\displaystyle d(x^{(n)},\Gamma )\leq \varepsilon }$:
4. ${\displaystyle x_{f}:=p_{\Gamma }(x^{(n)})}$, the orthogonal projection of ${\displaystyle x^{(n)}}$ on ${\displaystyle \Gamma }$
5. Return ${\displaystyle x_{f}}$

The result ${\displaystyle x_{f}}$ is an estimator of the first exit point from ${\displaystyle \Omega }$ of a Wiener process starting from ${\displaystyle x}$, in the sense that they have close probability distributions (see below for comments on the error).

In some cases the distance to the border might be difficult to compute, and it is then preferable to replace the radius of the sphere by a lower bound of this distance. One must then ensure that the radius of the spheres stays large enough so that the process reaches the border.^[1]

As it is a Monte-Carlo method, the error of the estimator can be decomposed into the sum of a bias, and a statistical error. The statistical error is reduced by increasing the number of paths sampled, or by using variance reduction methods.

The WoS theoretically provides exact (or unbiased) simulations of the paths of the Brownian motion. However, as it is formulated here, the ${\displaystyle \varepsilon }$-shell introduced to ensure that the algorithm terminates also adds an error, usually of order ${\displaystyle {\mathcal {O}}(\varepsilon )}$.^[4] This error has been studied, and can be avoided in some geometries by using Green's Functions First Passage method:^[5] one can change the geometry of the "spheres" when close enough to the border, so that the probability of reaching the border in one step becomes positive. This requires the knowledge of Green's functions for the specific domains. (see also Harmonic measure)

When it is possible to use it, the Green's function first-passage (GFFP) method is usually preferred, as it is both faster and more accurate than the classical WoS.^[4]

It can be shown that the number of steps taken for the WoS process to reach the ${\displaystyle \varepsilon }$-shell is of order ${\displaystyle {\mathcal {O}}(|\log(\varepsilon )|)}$.^[2] Therefore, one can increase the precision to a certain extent without making the number of steps grow notably.

As it is commonly the case for Monte-Carlo methods, this algorithm performs particularly well when the dimension is higher than ${\displaystyle 3}$, and one only needs a small set of values. Indeed, the computational cost increases linearly with the dimension, whereas the cost of grid dependant methods increase exponentially with the dimension.^[2]

This method has been largely studied, generalized and improved. For example, it is now extensively used for the computation of physical properties of materials (such as capacitance, electrostatic internal energy of molecules, etc.). Some notable extensions include:

The WoS method can be modified to solve more general problems. In particular, the method has been generalized to solve Dirichlet problems for equations of the form ${\displaystyle \Delta u=cu+f}$ ^[6] (which include the Poisson and linearized Poisson−Boltzmann equations) or for any elliptic partial differential equation with constant coefficients.^[7]

More efficient ways of solving the linearized Poisson–Boltzmann equation have also been developed, relying on Feynman−Kac representations of the solutions.^[8]

Again, within a regular enough border, it possible to use the WoS method to solve the following problem :

${\displaystyle {\begin{cases}\partial _{t}u(x,t)+{\frac {1}{2}}\Delta _{x}u(x,t)=0&{\mbox{if }}x\in \Omega {\mbox{ and }}t<T\\u(x,T)=h(x,T)&{\mbox{if }}x\in {\bar {\Omega }}\\u(x,t)=h(x,t)&{\mbox{if }}x\in \Gamma .\end{cases}}}$

of which the solution can be represented as:^[9]

${\displaystyle u(x,t)=\mathbb {E} _{t,x}(h(X_{T\wedge \tau },T\wedge \tau ))}$,

where the expectation is taken conditionally on ${\displaystyle X_{t}=x}$

To use the WoS through this formula, one needs to sample the exit-time from each sphere drawn, which is an independent variable ${\displaystyle \tau _{0}}$ with Laplace transform (for a sphere of radius ${\displaystyle R}$):^[10]

${\displaystyle \mathbb {E} (\exp(-s\tau _{0}))={\frac {R{\sqrt {2s}}}{\sinh(R{\sqrt {2s}})}}}$

The total time of exit from the domain ${\displaystyle \tau }$ can be computed as the sum of the exit-times from the spheres. The process also has to be stopped when it has not exited the domain at time ${\displaystyle T}$.

The WoS method has been generalized to estimate the solution to elliptic partial differential equations everywhere in a domain, rather than at a single point.^[11]

The WoS method has also been generalized in order to compute hitting times for processes other than Brownian motions. For example, hitting times of Bessel processes can be computed via an algorithm called "Walk on moving spheres".^[12] This problem has applications in mathematical finance.

The WoS can be adapted to solve the Poisson and Poisson–Boltzmann equation with flux conditions on the boundary.^[13]

Finally, WoS can be used to solve problems where coefficients vary continuously in space, via connections with the volume rendering equation.^[14]

• Feynman–Kac formula
• Stochastic processes and boundary value problems
• Euler–Maruyama method to sample the paths of general diffusion processes

• Sabelfeld, Karl K. (1991). Monte Carlo methods in boundary value problems. Berlin [etc.]: Springer-Verlag. ISBN 9783540530015.

• Some continuous Monte-Carlo methods for the Dirichlet problem The paper in which Marvin Edgar Muller introduced the method.
• Brownian Motion by Peter Mörters and Yuval Peres. See Chapter 3.3 on harmonic measure, Green's functions and exit-points.