Computation of Riesz $\alpha$-capacity $\mathrm{C}_{\alpha}$ of general sets in $\mathbb{R}^{d}$ using stable random walks^††thanks: Submitted to the editors August 19, 2025.

Mansfield et al. [22] introduced an algorithm for computing the 2-capacity in the special case where $\alpha=2$ and $d=3$ based on hitting the set $K$ with simple random walk trajectories. As noted above, these paths start uniformly on a launch sphere that contains $K$. We next describe a natural extension of this algorithm to general $\alpha$ and dimensionality $d$. Throughout we assume $\alpha<d$ to guarantee that the Lévy flights are transient (see Appendix for definition). The method requires two basic elements: a way to simulate isotropic stable random vectors and to be able to compute the probability of hitting a ball given that we start at an exterior point. The Appendix describes how to do both of these for $\alpha$-stable motion in $\mathbb{R}^{d}$, $0<\alpha\leq 2$.

Let $\gamma>0$ be a scale variable for the step size, i.e. each step $\bm{Y}$ of the random walk will have characteristic function

This isotropic stable process is described in Chapter 3 of Sato [36]. Let $\epsilon\geq 0$ be a hitting tolerance, i.e. we will consider that the walk has hit $K$ if the distance between a point and $K$ is less than or equal $\epsilon$. Typically $\epsilon\ll\gamma\ll R_{\text{launch}}$; the first inequality implies that the $\epsilon$ “thickening” of $K$ is small compared to the typical step size and the second inequality implies that the typical step size is small compared to the size of the object. We generate a path according to the following algorithm, see Figure 1(b).

1. 1.

   Pick a starting point $\bm{x}$. If $\alpha=2$, this means pick $\bm{x}$ uniformly on the launch sphere; if $\alpha<2$, this means pick $\bm{x}$ according to the equilibrium distribution on the launch ball using equation (9) below.

2. 2.

   Generate a step $\bm{y}$ with distribution given by equation (2) and compute a new location $\bm{x}\leftarrow\bm{x}+\bm{y}$.

3. 3.

   If $|\bm{x}|>R_{\text{escape}}$, compute the probability $p$ of a walk hitting the launch ball starting at $\bm{x}$. (Recall that the process is transient.) With probability $1-p$, declare that this walk will escape to infinity without ever hitting the launch ball (and therefore not hit $K$) and we stop this path. Otherwise, compute a new location $\bm{x}$ by simulating a hitting location in the launch ball by a walk starting at the current location as described in the Appendix.

4. 4.

   If $\bm{x}$ is within $\epsilon$ of $K$, declare the current location $\bm{x}$ to be a hit. Otherwise, go to step 2.

We note that since a stable process generally jumps into the interior of a set, we can set $\epsilon=0$ for the kinds of sets considered here. For thin sets, say fractals, we may want to use an $\epsilon>0$ to effectively “thicken” $K$, but the influence of this on the final results should be carefully examined.

This process is repeated $M$ times and a count of number of hits of $K$ and their locations are saved.

When $\alpha=2$, Muller [26] introduced a faster method to simulate from the hitting distribution of a general set $K$. The continuity of Brownian paths guarantees that a path can only leave a ball by crossing its boundary, a sphere. Furthermore, the rotational symmetry means that the distribution of exit points is uniform on the sphere. So rather than simulate many steps, this method chooses random exit points on a sequence of (random) spheres, see Figure 1(c). This is called the Walk-On-Spheres method, abbreviated WOS. Here are the steps in the algorithm to generate a single path.

1. 1.

   Pick a point $\bm{x}$ uniformly from the launch sphere.

2. 2.

   Compute $r\leftarrow|\bm{x}-K|=\inf_{\bm{k}\in K}|\bm{x}-\bm{k}|=$ distance between $\bm{x}$ and $K$.

3. 3.

   Pick a new point $\bm{y}$ uniformly on the sphere of radius $r$ centered on $\bm{x}$

4. 4.

   If $|\bm{y}-K|\leq\epsilon$, return $\bm{y}$ as the hitting location.

5. 5.

   Otherwise, set $\bm{x}\leftarrow\bm{y}$.

6. 6.

   If $|\bm{x}|>R_{\text{escape}}$, compute the probability $p$ of a walk starting at $\bm{x}$ hitting the launch ball. With probability $1-p$, declare that this walk will escape to infinity without ever hitting the launch ball (and therefore not hit $K$). Otherwise, go to step 2.

When $0<\alpha<2$, the Walk-On-Spheres method cannot be used because the stable process has jumps and will exit a ball not on the boundary, but on some point exterior to the ball. Blumenthal et al. [4] give expressions for the distribution of these exit point for $\alpha$-stable processes with $0<\alpha<2$. They also give the probability of hitting a ball when starting from some exterior point and what the hitting location is on the interior of the ball. These expressions are given in the Appendix.

The approach below extends the work of Kyprianou et al. [19]. That work considers the interior problem: a stable process starts at a point $\bm{x}\in K$, and the paths are followed until they exit $K$. In contrast, we are concerned with the exterior problem: in which we start at a point $\bm{x}\not\in K$, and propogate the path until it hits $K$. There are several key differences between [19] and the method proposed here; we describe them below.

Our approach generates a sequence of points that are outside a sequence of balls, with a check to avoid drifting too far from the object that sometimes returns us to the launch ball $R_{\text{launch}}\mathbb{B}$. We call this method the Walk-In-and-Out-of-Balls method, abbreviated WIOB. The steps of the method are given below, see Figure 1(d) to visualize the steps.

1. 1.

   Pick a point $\bm{x}$ on the launch ball according to the equilibrium distribution in equation (9).

2. 2.

   Compute the distance $r=|\bm{x}-K|$ = distance between $\bm{x}$ and $K$.

3. 3.

   Simulate an exit location $\bm{y}$ from the ball of radius $r$ centered on $\bm{x}$.

4. 4.

   If $|\bm{y}-K|\leq\epsilon$, return $\bm{y}$ as the hitting location.

5. 5.

   Otherwise, set $\bm{x}\leftarrow\bm{y}$.

6. 6.

   If $|\bm{x}|>R_{\text{escape}}$, compute the probability $p$ of a walk starting at $\bm{x}$ hits the launch ball. With probability $1-p$, declare that this walk will escape to infinity without ever hitting the launch ball (and therefore not hit $K$) and stop this iteration. Otherwise, simulate reentering the launch ball at a new point $\bm{x}$ and go to step 2.

There are two advantages of this method. First, like the WOS method for Brownian motion, it is generally more efficient. Instead of many small steps, the focus on leaving a ball means each WIOB step is larger and faster than many small steps. Second, WIOB allows a miss of the object $K$ that could occur with a jump that the WOS method does not allow. The chance of this happening is small, because small jumps are more likely, but the probability of a large step increases as $\alpha$ decreases.

We now describe differences between [19] and WIOB. Both methods rely on the key results in [4], but have different objectives. [19] is interested in simulating paths to compute a solution of a fractional diffusion equation and do not consider capacity, which is the main focus of this paper. They fix a starting location $\bm{x}$, whereas we start at a random location, which we specified above, chosen to exactly simulate hitting from infinity in step 1. They do not have to deal with transience and the possibility that the process never hits $K^{C}$, whereas we deal with this possibility in step 6. Simulating the return to the launch ball uses the density $f_{BGR}$ from the Appendix. The straightforward way of doing this with simple rejection is increasingly less efficient as the point becomes closer and closer to the launch ball. Devroye and Nolan [8] have developed algorithms for exactly simulating locations for entering or leaving a ball that are uniformly bounded for all starting locations $\bm{x}$. We note that these methods work for all $0<\alpha\leq 2$ and any dimension $d\geq 2$, enabling calculating standard capacity and Riesz capacity. This enables new applications in dimension not equal to 3.

For electrostatic capacity, e.g. $\alpha=2$ and $d=3$, then for any $|\bm{x}|>1$,

Based on this result, it can be shown that if $K$ is any compact set,

The sample analog is straightforward: run a large number of random walks and count the number of hits of $K$ and divide by the number of paths simulated. This is the essential idea underlying the ZENO program described in Mansfield et al. [22] and Mansfield and Douglas [21].

When $0<\alpha<2$, the situation is more complicated. For a ball, there is an expression for the hitting probability that is the analog of equation (3), see equation (8) in the Appendix below. However, there is no known relationship between hitting probability and $\alpha$-capacity for a general set $K$.

The proposed solution is to run $n$ simulations of the WIOB method to obtain hitting locations $\bm{x}_{1},\ldots,\bm{x}_{n}$ in the set $K$. Using these hitting locations gives an empirical distribution $\widehat{\mu}_{K}$ that is an approximation to the equilibrium measure of the set. The straightforward approach is to evaluate the energy integral (1) with a sum, i.e.

where $\binom{n}{2}$ is the number of distinct non-diagonal pairs in $\{\bm{x}_{i}-\bm{x}_{j}\}$.

Unfortunately, there are two problems with this seemingly straightforward approach. The first problem is that this sum is a numerically unreliable way to approximate the energy integral because the integrand in (1) has a singularity on the diagonal $\bm{x}=\bm{y}$. With a sample, if $\bm{x}_{i}\approx\bm{x}_{j}$, one or more terms in the sum above can dominate the sum, leading to a poor estimate of $I_{\alpha}(\mu)$. The second problem is that the terms in this sum are dependent, so it is not straightforward to get a confidence interval for such an estimate. Below, we describe a way to split the sample into two groups and use different estimators form the groups in a way that gives less volatile estimates of capacity and valid confidence intervals.

An equivalent way of viewing this is to define the random variable $W=\|\bm{X}-\bm{Y}\|^{\alpha-d}$, where $\bm{X}$ and $\bm{Y}$ are independent copies with distribution $\widehat{\mu}_{K}$. Then $I_{\alpha}(\widehat{\mu}_{K})=EW$, the expectation of $W$. Points that are close to each other on $K$ correspond to large values of $W$. In probabilistic terms, $W$ is a heavy-tailed variable and therefore estimating the mean can be unreliable. And contrary to most statistical estimates, the larger the sample, the more likely it is that we will see very large values of $W$ and get a less reliable estimate of $EW$. Furthermore, in cases where the $\mathrm{C}_{\alpha}(K)=0$, $EW=\infty$, but the sample mean will never lead to a value of $\infty$.

We propose a method to approximate $EW$ that detects the $EW=\infty$ case, gives reliable estimates of $I_{\alpha}(\widehat{\mu}_{K})$ when it is finite, and in this later case, gives a large sample confidence interval for $I_{\alpha}(\mu)$, and thus a confidence interval for $\mathrm{C}_{\alpha}(K)$. Start with a sample of $n$ hitting locations $\bm{x}_{1},\ldots,\bm{x}_{n}$ and define $w_{ij}=\|\bm{x}_{i}-\bm{x}_{j}\|^{\alpha-d}$, $i\neq j$. Pick a threshold $\tau>0$ for determining which of these values is large. A convenient choice is to define $\tau$ to be an upper quantile of $W$, e.g. $P(W\leq\tau)=p_{\tau}$ for some value $p_{\tau}\approx 1$. Split the integral for $EW$ into two pieces:

We will estimate $I_{1}$ using $n_{1}$ of the $n$ samples and $I_{2}$ by a tail estimator based on the remaining $n_{2}=n-n_{1}$ values.

For $I_{1}$, define $W^{*}=\min(W,\tau)$. This is a bounded random variable, so necessarily has a finite mean and variance. Then, define the sample analogue $w^{*}_{ij}=\min(w_{i,j},\tau)$ for $1\leq i\neq j\leq n_{1}$ and the following terms:

Then $\widehat{I}_{1}$ is a U-estimator of $I_{1}$. The variance term $v_{1}$ is used below to get a confidence interval for $I_{1}$.

For $I_{2}$, we assume that the upper tail of $W$ has a power decay:

To estimate $\nu$, select values $z_{1},\ldots,z_{n_{3}}$ from $w_{ij}$, $n_{1}<i,j\leq n$ where (a) $w_{ij}>\tau$, (b) there is at most one value selected from each row of the $w_{ij}$ matrix, and (c) there is at most one value selected from each column of the $w_{ij}$ matrix. Since we are only choosing the $z_{i}$ from $i,j>n_{1}$, the $z_{i}$ are independent of the $w_{ij}$ used in estimating $I_{1}$, the values $z_{i}$ are independent of $I_{1}$ and because of conditions (b) and (c) above, the $z_{i}$ values are independent of each other. Hence, the $z_{i}$ values are an i.i.d. sample from the upper tail of $W$. We use the Hill estimator for the tail decay $\nu$:

In the simulations we have run, we’ve observed $\widehat{\nu}\leq 2$ when $d\geq 3$ or ($d=2$ and $\alpha\leq 1$), in which case Var$(W)=\infty$. Assume that the tail approximation (5) is exact for $W>\tau$. Then $W$ has an upper tail density $c\nu w^{-\nu-1}$ and necessarily $c=(1-p_{\tau})\tau^{\nu}$. Then the value of $I_{2}$ is $\int_{\tau}^{\infty}(w-\tau)c\nu w^{-\nu-1}=(1-p_{\tau})\tau/(\nu-1)$ and the estimate of $I_{2}$ is

The energy integral estimate is $\widehat{I}=\widehat{I}_{1}+\widehat{I}_{2}$, and the $\alpha$-capacity estimate is $\widehat{\mathrm{C}_{\alpha}}(K)=1/\widehat{I}$. Note that if $\nu\leq 1$, the sum $I_{1}+I_{2}=\infty$ and $\widehat{\mathrm{C}_{\alpha}}(K)=0$.

The ZENO program has a method of deriving large sample confidence intervals when $\alpha=2$ and $d=3$. Large sample confidence intervals based on the normal distribution are described next for arbitrary $d\geq 2$ and $0<\alpha\leq 2$.

We assume $\nu>1$, so $I<\infty$. First, the variance of $\widehat{I}$ will be estimated. By the way the sample is split above, the two estimators $\widehat{I}_{1}$ and $\widehat{I}_{2}$ are independent, so Var$(\widehat{I})$ = Var$(\widehat{I}_{1})$ + Var$(\widehat{I}_{2})$. The standard theory of U-statistics, e.g. Chapter 12 of van der Vaart [43], gives an estimator of the variance of $I_{1}$ that is sometimes negative. (This possibility is mentioned on page 419 of Schucany and Bankson [38]; surprisingly, this can happen with large samples, e.g. $n_{1}=10000$.) To avoid that problem, we use the modified variance estimator of Sen [39]: namely $\sigma_{1}^{2}=\text{Var}(\widehat{I}_{1})=4v_{1}/n_{1}.$

With appropriate regularity conditions, the Hill estimator $\widehat{\nu}$ based on the $n_{3}$ largest values is approximately N$(\nu,\nu^{2}/n_{3})$, see page 304 of Resnick [32]. Using the delta method on $g(\nu)=(1-p_{\tau})\tau/(\nu-1)$ shows that $\widehat{I}_{2}$ is approximately N$(I_{2},\sigma_{2}^{2})$, with $\sigma_{2}^{2}=g^{\prime}(\nu)^{2}/n_{3}$. Thus $\widehat{I}=\widehat{I}_{1}+\widehat{I}_{2}$ is approximately N$(I,\sigma_{I}^{2})$, where $\sigma_{I}^{2}=\sigma_{1}^{2}+\sigma_{2}^{2}$. Finally, using the delta method again, the estimate of $\alpha-$capacity $\widehat{\mathrm{C}_{\alpha}}=1/\widehat{I}$ is approximately N$(\mathrm{C}_{\alpha},\sigma_{\mathrm{C}_{\alpha}}^{2})$, where $\sigma_{\mathrm{C}_{\alpha}}^{2}=\sigma^{2}_{I}/I^{4}$. Hence, a large sample $(1-\delta)$ confidence interval is given by $\widehat{I}\pm z_{\delta/2}\sigma_{I}/\widehat{I}^{2}$.

The only case where an exact formula is known for the $\alpha$-capacity of a $d$-dimensional set is for a ball; it is given in equation (7) below. Here we focus on the unit ball in dimension $d=2,3,4,5$. This first simulation used the WIOB method to estimate $\alpha$-capacity and then compares those estimates to the theoretical values. Figure 2 compares these numerical estimates to the exact theoretical values in equation (7) in the Appendix. Since the exact $\alpha$-capacity of a sphere is known, we may gauge the accuracy of the simulation by the deviation between estimates based on simulation and analytic, exact values. For example, when $d=3$, the confidence interval width varies from 0.0948 when $\alpha=0.25$ to 0.00318 when $\alpha=2$. Increasing the sample size will improve the accuracy.

First let $K$ be the solid cube $[0,1]^{d}$. Figure 3 shows estimated $\mathrm{C}_{\alpha}$ of a these cubes in dimensions 2, 3, 4, and 5 as a function of $\alpha$. In particular, the estimated value of this capacity when $d=3$ and $\alpha=2$ is $0.66134\pm 9.7\times 10^{-4}$ using WIOB with $n=10000$. This compares to previous estimates of 0.663 in Mansfield et al. [22], $0.660675\pm 10^{-5}$ using 4.7 billion trajectories in Given et al. [12] and $0.6606782\pm 10^{-7}$ in Hwang et al. [14]. The particular virtue of the path integration method for estimating capacity is that it allows for highly precise estimates in comparison with other methods such as finite element methods when a large number of trajectories is employed.

Separate simulations, not shown here, demonstrate that the method introduced here satisfies the scaling property given in equation (6) in the Appendix: $\mathrm{C}_{\alpha}(rK)=r^{d-\alpha}\mathrm{C}_{\alpha}(K)$ for $r>0$. So these capacity values for unit cubes give the capacity of any cube.

Since non-Gaussian stable processes exhibit discontinuous path displacements, their hitting locations are not restricted to the surface of an object. This class of erratic path processes can thus penetrate through the boundary before they hit $K$. So unlike 2-capacity, the $\alpha$-capacity of a hollow object can differ from the $\alpha$-capacity of a solid object with the same bounding surface. To examine this effect in 2-dimensions, we look at a solid and hollow squares. Figure 4 shows these three regions: the first (upper left) is a single, solid square $[-1/2,1/2]^{2}$ (top left), the next is a hollow square with a thick wall of width 1/3 (top right), and one with a thin wall of width 1/100 (bottom left). The bottom right plot shows the estimated $\alpha$-capacity of the three objects for a range of $\alpha$ in (0,2]. The black curve for the solid cube is mostly covered by the thick walled curve in red, showing that it is difficult to get through a thick barrier when most of the jumps of a stable process are small. For the thin-walled cube, the green curve does diverge from the other two curves, and this difference is larger when $\alpha$ is small. In theory, this allows one to distinguish between solid and hollow objects with the same surface, something not possible with Brownian motion.

Here we consider an anisotropic example in $\mathbb{R}^{3}$: $K_{R}=[-R,R]\times[-1,1]^{2}$. This is a solid bar with length $2R$ and cross section a square of area 4, so it has volume $8R$. As $R$ varies away from 1 in either direction, the set $K_{R}$ gets more and more anisotropic. Since the volume of $K_{R}$ increases with $R$, we normalize each capacity estimate by dividing it by the $\alpha$-capacity of a ball with the same volume as $K_{R}$. This type of dimensionless ratio, the relative capacity, provides valuable information about the shape of a region, e.g., Polya and Szego [29]. Figure 5 shows the results of simulating hitting locations and computing $\alpha-$capacity for a range of $\alpha$ and $R$. The increase in relative capacity near $R=0$ is caused by the fact that $\mathrm{C}_{\alpha}(K_{R})$ decreases to 0 more slowly than $\mathrm{C}_{\alpha}$(ball with the same volume as $K_{R})$. The subplot in the figure shows that the WIOB method agrees closely with the ZENO program when $\alpha=2$.

Recently, there has been great interest in developing a ”fingerprint” of the shape of objects through the calculation of the eigenvalues of the region based on solving the wave equation for objects of general shape, e.g. Reuter et al. [33] and Peincke et al. [28]. The eigenvalue spectrum is clearly advantageous to other shape functionals as it involves the calculation of an infinite (large number) of numbers to specify object shape. The $\alpha-$capacity likewise involves a “spectrum” of values that might better quantify object shape. We note that the relative $\alpha-$capacity obeys the same isoperimetric property (see the Appendix), i.e., the $\alpha$-capacity quantity is minimized by a sphere of all regions K of a fixed volume. It is this isoperimetric property that makes the relatively capacity of interest in shape classification and discrimination problems, as discussed for the 2-capacity and other shape functionals having this fundamental isoperimetric property.

Here we consider a thin “coin” $K=\{\bm{x}\in\mathbb{R}^{3}:x_{2}^{2}+x_{3}^{2}\leq 1,-\epsilon\leq x_{1}\leq\epsilon\}$. Since $\epsilon$ will be small, $K$ is approximately a two dimensional set sitting inside of $\mathbb{R}^{3}$. This example illustrates two topics: using capacity to find Hausdorff dimension and subordination.

For the first topic, we calculate the $\alpha$-capacity of the coin in three dimensions for a range of $\alpha$ values and see where the capacity reaches 0. Figure 6 shows the result, where $\alpha=1$ defines a condition at which the $\alpha$-capacity reaches 0. This behavior arises because the Hausdorff dimension of the coin is $d-\alpha=3-1=2$. (This conclusion follows from work by Taylor [41], using a lemma of Frostman [11] for the calculation of Hausdorff dimension of Brownian motion paths and other “fractal” sets. That work shows that the value of $\alpha$ where $\mathrm{C}_{\alpha}(K)$ reaches 0 rigorously defines the Hausdorff dimension of $K$.) Correspondingly, the 2-capacity of the thin coin vanishes in four dimensions in the case of Brownian motion, since the thin ”coin” and Brownian motion have the same Hausdorff dimension, see Taylor [41], Dvoretzky et al. [10], and McKean [23].

For the other topic, consideration of the hitting locations of Brownian motion ($\alpha=2$) in $\mathbb{R}^{3}$ on an infinitely thin “coin” provides an example of subordination, the reduction of Brownian motion to a Lévy type process by “decimating” the Brownian motion at random time points. In particular, the times at which the random walk in three dimensions hits a disc defines a random Cantor set of time points of fractal dimension 1/2 (the subordinating set) and the path process connecting the original random walk path positions at these time points defines a Cauchy process ($\alpha=1$) where the paths are then limited by construction to the unit disk in $\mathbb{R}^{2}$. We use the WOS method for Brownian motion to generate a sample hitting distribution $\bm{x}_{1},\ldots,\bm{x}_{n}\in\mathbb{R}^{3}$ of the unit disk. We then perform a separate WIOB simulation for the unit disk in $\mathbb{R}^{2}$ with an $\alpha=1$ path process and determine hitting locations $\bm{y}_{1},\ldots,\bm{y}_{n}\in\mathbb{R}^{2}$. To compare these to each other, the radial symmetry says it suffices to compare the distributions of univariate radii $r_{i}=\sqrt{x_{i2}^{2}+x_{i3}^{2}}$ (ignoring the small first coordinate) to univariate $s_{i}=\sqrt{y_{i1}^{2}+y_{i2}^{2}}$. Figure 7 shows the cumulative distributions of these radii. Also on the plot is a third curve, which is the exact cumulative distribution of the radii of the Cauchy equilibrium distribution derived from (9) with $d=2$ and $\alpha=1$. The view that the two dimensional Cauchy process is just a subordinated version of the three dimensional Brownian path process, as described above, implies that the hitting distributions of the disk should then be the same. We find that this might arise naturally when viewing the dynamics of a higher dimensional Brownian process in some space of reduced dimension. We do not claim this is a universal explanation, but this may be a common mechanism.

To assess the performance of the confidence interval method of Section 3.3, we simulated $M=500$ runs with $K$ the unit ball in dimension $d=3$ and three $\alpha$ values. The exact value of the $\alpha$-capacity is known in this case, so we counted how many of the 500 simulated 95% confidence intervals contained the exact value. For $\alpha=0.5$, the observed coverage probability was 0.992; for $\alpha=1.2$, the observed coverage probability was 0.966; and for $\alpha=1.9$, the observed coverage probability was 0.958. The other parameters are $p_{\tau}=0.99$, $n=10000$, $n_{1}=n/2=5000$, $R_{\text{launch}}=2$, $R_{\text{escape}}=4$. The observed coverage probabilities are high, apparently due to the large values of $v_{1}$. It appears that even with the truncation at values of $w_{ij}$ above $\tau$, the large values of $w_{ij}^{*}$ inflate the estimator $v_{1}$. This will be investigated in future work. Presently, observe that the confidence intervals are conservative, and that a larger simulation can be used to get smaller confidence intervals.