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.^[1][2] The WoS method was first introduced by M. 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, by simulating paths of Brownian motion (or for some more general variants, diffusion processes), and it is today one of the most widely used "grid-free" algorithms for generating Brownian paths.

Informal description

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

Let us consider the Dirichlet problem:

$\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^{[lower-alpha 1]} that when the solution $u$ exists, for $x \in \Omega$:

$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:

$\mathbb{E}_x[ h(W_\tau) ] \sim \frac{1}{n} \sum_{i=1}^n h(W^i_\tau)$

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. Rather than simulating in detail the path of the process, it samples only the exit-points out of successive spheres, often making it less costly than "grid-based" algorithms.

The WoS algorithm consists in drawing the largest sphere $\mathcal{S}_0$ centered on x and contained inside the domain. The first point of exit x₁ from $\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 $\varepsilon$-shell, or $\varepsilon$-layer.^[4]

Formulation of the method

Illustration of a run of the Walk-on-spheres algorithm on an arbitrary domain $\Omega$ with an $\varepsilon$-shell

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

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

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

Comments and practical considerations

Radius of the spheres

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]

Bias in the method and GFFP

The Walk-on-spheres method is used until the process gets $\delta$-close to the border. Then the sphere is replaced by its "intersection" with the boundary of the domain.

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 $\varepsilon$-shell introduced to ensure that the algorithm terminates also adds an error, usually of order $\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 GFFP method is usually preferred, as it is both faster and more accurate than the classical WoS.^[4]

Complexity

It can be shown that the number of steps taken for the WoS process to reach the $\varepsilon$-shell is of order $\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 $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]

Variants, extensions

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:

Elliptic equations

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 $\Delta u = cu + f$ ^[6] (which include the Poisson and linearized Poisson−Boltzmann equations) or for any elliptic 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]

Time dependency

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

$\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]

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

where the expectation is taken conditionally on $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 $\tau_0$ with Laplace transform (for a sphere of radius $R$):^[10]

$\mathbb{E}(\exp(- s \tau_0))) = \frac{R \sqrt{2s}}{\sinh(R \sqrt{2s})}$

The total time of exit from the domain $\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 $T$.

Other extensions

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.

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

See also

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

Notes

1. The link was first established by Kakutani for the 2-dimensional Brownian motion,^[3] it can now be seen as a trivial case of the Feynman−Kac formula.

References

1. Muller, Mervin E (Sep 1956). "Some continuous Monte-Carlo Methods for the Dirichlet Problem". The Annals of Mathematical Statistics 27 (3): 569–589. doi:10.1214/aoms/1177728169.
2. Sabelfeld, Karl K. (1991). Monte Carlo methods in boundary value problems. Berlin [etc.]: Springer-Verlag. ISBN 3540530010.
3. Kakutani, Shizuo (1944). "Two-dimensional Brownian motion and harmonic functions". Proceedings of the Imperial Academy 20 (10): 706–714. doi:10.3792/pia/1195572706.
4. Mascagni, Michael; Hwang, Chi-Ok (June 2003). "ϵ-Shell error analysis for "Walk On Spheres" algorithms". Mathematics and Computers in Simulation 63 (2): 93–104. doi:10.1016/S0378-4754(03)00038-7.
5. Given, James A.; Hubbard, Joseph B.; Douglas, Jack F. (1997). "A first-passage algorithm for the hydrodynamic friction and diffusion-limited reaction rate of macromolecules". The Journal of Chemical Physics 106 (9): 3761. doi:10.1063/1.473428.
6. Elepov, B.S.; Mikhailov, G.A. (January 1969). "Solution of the dirichlet problem for the equation Δu − cu = −q by a model of "walks on spheres"". USSR Computational Mathematics and Mathematical Physics 9 (3): 194–204. doi:10.1016/0041-5553(69)90070-6.
7. Booth, Thomas E (February 1981). "Exact Monte Carlo solution of elliptic partial differential equations". Journal of Computational Physics 39 (2): 396–404. doi:10.1016/0021-9991(81)90159-5.
8. Hwang, Chi-Ok; Mascagni, Michael; Given, James A. (March 2003). "A Feynman–Kac path-integral implementation for Poisson’s equation using an h-conditioned Green’s function". Mathematics and Computers in Simulation 62 (3-6): 347–355. doi:10.1016/S0378-4754(02)00224-0.
9. Gobet, Emmanuel (2013). Méthodes de Monte-Carlo et processus stochastiques du linéaire au non-linéaire. Palaiseau: Editions de l'Ecole polytechnique. ISBN 978-2-7302-1616-6.
10. Salminen, Andrei N. Borodin; Paavo (2002). Handbook of Brownian motion : facts and formulae (2. ed.). Basel [u.a.]: Birkhäuser. ISBN 978-3-7643-6705-3.
11. Booth, Thomas (August 1982). "Regional Monte Carlo solution of elliptic partial differential equations". Journal of Computational Physics 47 (2): 281–290. doi:10.1016/0021-9991(82)90079-1.
12. Deaconu, Madalina; Herrmann, Samuel (December 2013). "Hitting time for Bessel processes—walk on moving spheres algorithm (WoMS)". The Annals of Applied Probability 23 (6): 2259–2289. doi:10.1214/12-AAP900.
13. Simonov, Nikolai A. (2007). "Random Walks for Solving Boundary-Value Problems with Flux Conditions". Numerical Methods and Applications 4310: 181–188. doi:10.1007/978-3-540-70942-8_21.

Further reading

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

External links

• 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.