Walking on half-spaces

A Monte Carlo approach to solve fractional PDEs

Tags: Fractional analysis, Dirichlet problem, Feynman-Kac, Martingales, \alpha-stable-process, Laplace equation, fractional Laplacian.

Vanilla Laplacian and Brownian motion

Every mathematician is familiar with the Laplace equation.

$$\Delta f=0.$$

where $\Delta f = \nabla \cdot \nabla f$  is called the Laplacian operator and $\nabla f = \left ( \frac{\partial f }{\partial x_1} , \ldots , \frac{\partial }{\partial f x_d} \right )$ (if you’re not a mathematician or if you need some help remembering, see references 1 and 2).

This second order Partial Differential equation equation is a cornerstone of classical mathematics and physics. Among its applications, it is used to describe the behavior of scalar fields such as gravitational, electric, and fluid potentials in a state of equilibrium. In probability, it is closely related with the concept of continuous martingales, and hence it is an important concept in mathematical finance.

Adentrémonos un poco más en el lado probabilistico,

The Dirichlet problem for Laplace's equation, an ubiquitous problem in mathematical modelling, consists of finding a solution $f$ to the Laplace equation on some domain $D$, such that $f$ on the boundary of $D$ is equal to some given function $g$. i.e.,

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

Explicit solutions to this problem are difficult to find, and frequently, numerical methods are employed to approximate its solutions. The Finite Element method is one of the most used techniques in this matter. However, there is a Monte Carlo approach, based on a probabilistic interpretation of the problem that allows us to solve numerically without big complications.

If $X=(X_t:t\geq 0)$ with law $\mathbb{P}$, is a Brownian Motion in $\mathbb{R}^d$, and $\tau_D=\inf\{t>0:X_t\notin D\}$, and $g$ is a bounded continuous function defined on $D\cup\partial D$, then for $x \in D\cup\partial D$

$$u(x)=\mathbb{E}\left[ g(X_{\tau_D}) |X_{0}=x\right]$$

is the only continuous solution to the Dirichlet problem.

By the ‘law of the unconscious statistician, the function can be written as

$$u(x)=\int_{\partial D} \mathbb{P}(X_{\tau_D}\in {\rm d}y | X_0=0) g(y),$$

and, this representation this gives us a way of writing every possible solution of the Dirichlet problem. The distribution $\mathbb{P}(X_{\tau_D}\in {\rm d}y | X_0=0)$, is a key object in the theory of harmonic functions (those who solve the Laplace equation on $D$), and it is known as the Poisson kernel of $D.$

Fractional Laplacian and the $\alpha$-stable process

Many real-world systems exhibit dynamics that cannot be adequately captured by integer-order derivatives. Once again, in probability, there are discontinuous processes for which a ‘martingale’ behaviour leads us to the consideration of an oddbod known as the fractional Laplacian, and an analogous fractional Laplace equation. (Yes, fractional calculus and fractional analysis are a thing)

The fractional Laplacian extends the traditional Laplacian operator to non-integer orders of differentiation. It is often denoted as $(-\Delta)^{\alpha/2},\, 0<\alpha<2,$ and, by capturing the notion of derivatives of non-integer order, it provides a more nuanced and flexible approach to modeling various phenomena (See reference 3).

Its applications are far-reaching, impacting fields such as anomalous diffusion processes, image processing, finance, and quantum mechanics.

In this blog entry I introduce the fractional Laplace equation and a method to solve it numerically in:

1. The exterior of a bounded convex domain, this with the walk-on-spheres algorithm (see reference 4)

2. The interior of a bounded convex domain using the walk-on-half-spaces

The Generator and the Isotropic $\alpha$-stable process

The crucial link between the probabilistic interpretation and the analytic problem is that the fractional Laplacian is the infinitesimal generator of a well-known stochastic process: the isotropic $\alpha$-stable Lévy process. Just as the standard Brownian motion has the classical Laplacian operator as its generator, leading to connections between Brownian motion and harmonic functions, the $\alpha$-stable process (with its discontinuous, jump-like paths) governs the solutions to the fractional Laplace equation. This deep connection, explored thoroughly in potential analysis (see, for example, Bogdan, Burdzy et al., "Potential Analysis of Stable Processes and its Extensions"), allows us to use simulated paths of the $\alpha$-stable process to evaluate boundary value problems numerically.

1. Walk-on-spheres for the exterior domain

The Walk-on-spheres (WoS) algorithm, originally devised for the classical Dirichlet problem using Brownian motion, can be extended to the fractional setting for computing solutions outside a bounded domain. Instead of moving continuously and generating a spherical surface jump when hitting the local sphere, the $\alpha$-stable process "jumps" across space. For the exterior problem, we sample the exit position of the $\alpha$-stable process from successive spheres that do not intersect the boundary. By chaining these jumps, the random walk eventually escapes to the target boundary region, allowing us to approximate the expectation via Monte Carlo integration (Kyprianou et al., 2018).

2. Walk-on-half-spaces for the interior domain

While Walk-on-spheres relies on isotropic jumps from the center of successive maximal spheres, the Walk-on-half-spaces (WoHS) approach constructs the random walk by utilizing half-spaces to bound the steps. This is particularly advantageous for the interior of bounded convex domains. At each iteration, the algorithm samples the first hitting point of an appropriately chosen bounding half-space by an $\alpha$-stable process. Because the domain is convex, we can bound it with tangent hyperplanes, and the analytic knowledge of the jump distribution of the stable process crossing a hyperplane makes this sampling efficient.

References

1. Gilbert Strang. 2015. Laplace Equation [Video]. YouTube. URL

2. Khan Academy. 2016. Laplacian Intuition [Video]. YouTube. URL

3. Kwaśnicki, M. (2017). Ten equivalent definitions of the fractional Laplace operator. Fractional Calculus and Applied Analysis, 20(1), 7-51.

4. Kyprianou, A. E., Osojnik, A., & Shardlow, T. (2018). Unbiased ‘walk-on-spheres’ Monte Carlo methods for the fractional Laplacian. IMA Journal of Numerical Analysis, 38(3), 1550-1578

5. Qi, Y., Seyb, D., Bitterli, B., & Jarosz, W. (2022, July). A bidirectional formulation for Walk on Spheres. In Computer Graphics Forum (Vol. 41, No. 4, pp. 51-62).