1 Introduction↩︎

The walk on spheres (WoS) algorithm is a grid free way to solve some differential equations at points \(\boldsymbol{z}\) in a closed region \(\Omega\subset\mathbb{R}^d\). It dates back to [1], with many early results included in the monograph of [2]. It has historically been used less often than the finite element method (FEM) but there has been a resurgence of interest in it, with [3] being a notable recent publication. They list several advantages of WoS over FEM, including: WoS can provide a solution at just the most important points in \(\Omega\) without requiring one to solve the equations everywhere, WoS saves the cost of generating an FEM grid which can be quite high when \(\Omega\in\mathbb{R}^3\) is bounded by millions of triangles, and as a Monte Carlo method, WoS is easily parallelized.

In this paper we will use randomized quasi-Monte Carlo (RQMC) in WoS computations. There was earlier work using quasi-Monte Carlo (QMC) to solve WoS problems by [4] and [5]. The WoS integrands typically have infinite variation in the sense of Hardy and Krause. We say that they are ‘not in BVHK’. This infinite variation complicates the theoretical treatment of their QMC error. In RQMC, by contrast, we only need the integrand to be in \(L_2\) in order to ensure a mean squared error of \(o(1/n)\) on \(n\) function evaluations.

What we typically see in our examples is a mean squared error (MSE) for RQMC that, for \(n\) sample points, decreases in proportion to a power of \(n\) slightly below \(-1\), commonly comparable to \(n^{-1.1}\) or \(n^{-1.2}\). A rate of \(n^{-1.2}\) can yield a variance reduction of five to ten fold without using extremely large sample sizes. The WoS integrands are discontinuous. It is common for them to be indicator functions of sets that do not have axis parallel boundaries, which is why they are not in BVHK. To explain the variance reductions we see, we draw on and extend some work of [6] and [7]. These papers show improved convergence rates when the integrand is the indicator function of a set \(\Theta\subset[0,1)^d\) whose boundary \(\partial\Theta\) has finite \(d-1\) dimensional Minkowski content, or when it is a sufficiently regular function multiplied by such an indicator. For a \(d\) dimensional integrand they get a variance of \(O(n^{-1-1/d})\). A major part of this paper is devoted to giving sufficient conditions for the sets in WoS to have finite Minkowski content. In our setting, the integrand also has unbounded dimension, though it is commonly truncated to finite dimension.

The path to finite Minkowski content is long because even if \(\partial\Omega\) is made up of infinitely differentiable curves there can be severe pathologies as shown by [8] whose work we describe below. They rule out some of those pathologies using real analytic curves. In our analysis of the WoS setting, the distance to subsets of \(\partial\Omega\) is critically important. Under our conditions, that function will be piecewise \(C^\infty\) but it will only be Lipschitz globally and this complicates some of the manifolds that we study.

This paper is organized as follows. Section 2 gives a basic account of WoS to orient the reader and a brief note on the RQMC methods we use. Section 3 illustrates the use of RQMC on the WoS problem for a gasket example taken from [9]. There we see that RQMC methods have an MSE that decays faster than \(n^{-1}\) for \(n\leqslant 2^{17}\). Section 4 provides background on Minkowski content, rectifiability and the medial axes of domains in \(\mathbb{R}^2\). Section 5 lists some theorems from differential geometry that we need. Section 6 has our main theorem. Letting \(\Theta_k\subset[0,1)^k\) represent the set of RQMC points that cause WoS to reach a given portion of \(\partial\Omega\subset\mathbb{R}^2\) in exactly \(k\) steps, we give sufficient conditions for \(\partial \Theta_k\) to have finite \(k-1\) dimensional Minkowski content. Section 7 discusses how the main theorem applies to RQMC problems. Section 8 briefly describes more numerical examples including three dimensional domains that are not covered by our theorems. For four RQMC methods and five examples we see modest variance reduction factors ranging from \(1.8\) to \(10.7\) at \(n=2^{17}\). The code used for our numerical examples is publicly available at https://github.com/hoval58/RQMC-WoS. Section 9 has some final comments.

2 Background on WoS and RQMC↩︎

Here we describe WoS and RQMC. For an (R)QMC readership we present WoS in a way that shows that a random walk algorithm is very natural and then we use RQMC in that random walk. For this readership we only emphasize a few RQMC details that are specific to this project.

Here is some notation that we use throughout the paper. The positive integers are denoted by \(\mathbb{N}\) and \(\mathbb{N}_0=\mathbb{N}\cup\{0\}\). For a set \(\Theta\subset\mathbb{R}^d\) we use \(\Theta^\circ\) for its interior, \(\overline{\Theta}\) for its closure and \(\partial\Theta\) for its boundary. If \(\Theta\) is Lebesgue measurable, then \(\mathrm{vol}(\Theta)\) is its measure. For nonempty \(\Theta\subset\mathbb{R}^d\) and \(\boldsymbol{z}\in\mathbb{R}^d\) we define the distance function \[\mathrm{dist}(\boldsymbol{z},\Theta) = \inf_{\tilde{\boldsymbol{z}}\in\Theta}\Vert\boldsymbol{z}-\tilde{\boldsymbol{z}}\Vert.\] For fixed \(\Theta\), \(\mathrm{dist}(\boldsymbol{z},\Theta)\) is a Lipschitz function of \(\boldsymbol{z}\). If \(\Theta\) is closed, then the infimum is attained, i.e., \(\mathrm{dist}(\boldsymbol{z},\Theta)=\Vert\boldsymbol{z}-\bar\boldsymbol{z}\Vert\) where \(\bar\boldsymbol{z}\in\Theta\) is a (not necessarily unique) projection of \(\boldsymbol{z}\) onto \(\Theta\). We use \(\boldsymbol{1}_\Theta\) for the function taking the value \(1\) on \(\Theta\) and \(0\) on \(\Theta^c\).

For \(\varepsilon\geqslant 0\), we will use balls and spheres \(B(\boldsymbol{z},\varepsilon) =\{\tilde{\boldsymbol{z}}\in\mathbb{R}^d\mid \Vert\tilde{\boldsymbol{z}}-\boldsymbol{z}\Vert\leqslant\varepsilon\}\) and \(S(\boldsymbol{z},\varepsilon)=\{\tilde{\boldsymbol{z}}\in\mathbb{R}^d\mid\Vert\tilde{\boldsymbol{z}}-\boldsymbol{z}\Vert =\varepsilon\}=\partial B(\boldsymbol{z},\varepsilon)\), respectively. We follow [8] in allowing \(\varepsilon=0\), and then \(B(\boldsymbol{z},0)=S(\boldsymbol{z},0)=\{\boldsymbol{z}\}\). When we need to indicate the dimension, then we use \(B_d\) and \(S_d\).

2.1 Basic walk on spheres↩︎

The problem in WoS is to solve a differential equation at some point \(\boldsymbol{z}_0\) inside a bounded domain \(\Omega\subset\mathbb{R}^d\). In the simplest setting, we want to compute \(u(\boldsymbol{z}_0)\) for a harmonic function \(u:\Omega\to\mathbb{R}\) at a point \(\boldsymbol{z}_0\in\Omega^\circ\), and we are given a boundary function \(h\), so that \(u(\boldsymbol{z})=h(\boldsymbol{z})\) for \(\boldsymbol{z}\in\partial\Omega\). This is the Dirichlet boundary condition. The Neumann boundary condition specifies the outward normal derivative on \(\partial \Omega\), and hybrid combinations of these two types also exist. However, in this paper we only consider boundary value problems (BVPs) with a Dirichlet condition.

A harmonic function \(u\) has Laplacian \(\Delta u(\boldsymbol{z}) =\sum_{j=1}^d\partial^2u(\boldsymbol{z})/\partial z_j^2=0\). For standard coordinate vectors \(\boldsymbol{e}_j\), writing \[0=\Delta u(\boldsymbol{z}_0)\approx \frac{1}{\varepsilon^2}\sum_{j=1}^d \bigl(u(\boldsymbol{z}_0+\varepsilon \boldsymbol{e}_j)-2u(\boldsymbol{z}_0)+u(\boldsymbol{z}_0-\varepsilon\boldsymbol{e}_j)\bigr)\] shows that \(u(\boldsymbol{z}_0)\) is approximately the average of \(u\) at \(2d\) neighbors. We could pick one of those neighbors (call it \(\boldsymbol{z}_1\)) at random and then \(u(\boldsymbol{z}_1)\) is a nearly unbiased estimate of \(u(\boldsymbol{z}_0)\). The value \(u(\boldsymbol{z}_1)\) is nearly the average of its \(2d\) neighbors and we could keep randomly selecting neighbors until we get \(\boldsymbol{z}_k\) in (or very close to) \(\partial\Omega\) where \(u\) is known. Now \(u(\boldsymbol{z}_k)\) is a nearly unbiased estimate of \(u(\boldsymbol{z}_0)\).

As \(\varepsilon\to0\) the random walk converges to a Brownian motion in \(\Omega\) that first hits \(\partial\Omega\) at some time \(\tau>0\). Then \(u(\boldsymbol{z}_\tau)\) is an unbiased estimate of \(u(\boldsymbol{z}_0)\). That is, \(u(\boldsymbol{z}_0)=\mathbb{E}(h(\boldsymbol{z}_\tau))\) where \(\tau = \inf\{t\in[0,\infty)\mid \boldsymbol{z}_t\in\partial\Omega\}\).

The key insight in WoS is that when a Brownian motion exits the ball \(B(\boldsymbol{z}_k,r)\) for \(r>0\), it exits with a uniform distribution over \(\partial B(\boldsymbol{z}_k,r)\). Instead of generating a Brownian motion starting at \(\boldsymbol{z}_0\), we can simply take \(\boldsymbol{z}_{k}\sim\mathbf{U}(S(\boldsymbol{z}_{k-1},r_k))\) for integers \(k\geqslant 1\) where \(r_k=\mathrm{dist}(\boldsymbol{z}_{k-1},\partial\Omega)\).

For most geometries, the WoS algorithm will never produce \(\boldsymbol{z}_k\in\partial\Omega\). Instead we stop the WoS when \(\boldsymbol{z}_k\in\partial\Omega_{\varepsilon}:=\{\boldsymbol{z}\in\Omega\mid \mathrm{dist}(\boldsymbol{z},\partial\Omega)<\varepsilon\}\). This happens at step \(\tau = \min\{k\in\mathbb{N}_0\mid \boldsymbol{z}_k\in\partial\Omega_{\varepsilon}\}\). We then return the value \(h(\bar\boldsymbol{z}_{\tau})\) where for \(\boldsymbol{z}\in\Omega\), \(\bar\boldsymbol{z}\) is the projection of \(\boldsymbol{z}\) onto \(\partial\Omega\) with an arbitrary tie-breaker rule when that projection is not unique. The \(i\)-th run in an MC approach to WoS generates a trajectory \(\boldsymbol{z}_{i,k}\) for \(1\leqslant k\leqslant\tau_i\) starting at \(\boldsymbol{z}_{i,0}=\boldsymbol{z}_0\). Then a Monte Carlo WoS estimate of \(u(\boldsymbol{z}_0)\) is \[\begin{align} \label{eq:mcwos} \hat{u}(\boldsymbol{z}_0)=\frac{1}{n}\sum_{i=1}^n h(\bar \boldsymbol{z}_{i,\tau_i}). \end{align}\tag{1}\]

A recent innovation of [10] pools information over a set of nearby starting points. When \(\tilde{\boldsymbol{z}}_0\in B(\boldsymbol{z}_0,r_1)\), the authors apply an importance sampling weight to the distribution of \(\boldsymbol{z}_1\) and then reuse the entire walk for \(\boldsymbol{z}_0\) in their estimate of \(u(\tilde{\boldsymbol{z}}_0)\). They do not need to apply weights for \(\boldsymbol{z}_k\) with \(k\geqslant 2\).

2.2 Bias and running time↩︎

A WoS algorithm stops when \(\mathrm{dist}(\boldsymbol{z}_k,\partial\Omega)<\varepsilon\). This raises questions about how the bias from stopping short of the boundary and how the number of steps to reach the boundary both depend on \(\varepsilon\). [11] give a very complete analysis of the number of steps to reach \(\partial\Omega\subset\mathbb{R}^d\). They have a notion of \(\alpha\)-thickness for \(\alpha\in[0,d]\) where \(\alpha>0\) allows for fractal behavior of \(\partial\Omega\). For any \(\alpha<2\), the expected number of steps to get close to \(\partial\Omega\) is \(O(\log(1/\varepsilon))\) as \(\varepsilon\to0\). For \(\alpha=2\) the rate is \(O( \log^2(1/\varepsilon))\). For \(\alpha>2\) it is \(O((1/\varepsilon)^{2-4/\alpha})\). The domains in our examples are all \(0\)-thick. That can be verified by applying Definition 1 of [11]. As a result, the number of steps taken are all \(O( \log(1/\varepsilon))\). In addition to a bound on the expected number of steps, they also consider the tail probability. Their Remark \(3\) gives an exponential decay for the probability that more than \(k\) steps are required.

The present understanding of the bias in WoS is less well developed than the convergence time. There is extensive empirical work in [4] where the bias was seen to be \(O(\varepsilon)\) in every one of a range of examples. They reason that the problem has distorted \(\partial\Omega\) by \(O(\varepsilon)\) and in that case a bias of \(O(\varepsilon)\) is natural. [3] report on the small bias in WoS in some large scale computations.

2.3 Nonzero source term↩︎

A more general boundary value problem has \(u(\boldsymbol{z})\) known for \(\boldsymbol{z}\in\partial\Omega\) and subject to \(\Delta u(\boldsymbol{z})=g(\boldsymbol{z})\) for \(\boldsymbol{z}\in\Omega\). When the source term \(g\) is nonzero, then the solution at \(\boldsymbol{z}\in\Omega^\circ\) involves some Green’s functions \(G_d^B\) that we define below. We follow the derivation in [12]. They use \(\Delta\) to represent the negative semi-definite Laplacian, so their \(\Delta u\) is our \(-\Delta u\). In order to use their presentation we replace our source function \(g\) by \(-g\). The equation (11) of [12] gives us \[\begin{align} u(\boldsymbol{z}_0) = \frac{1}{\mathrm{vol}(S(\boldsymbol{z}_0,r_1))}\int_{S(\boldsymbol{z}_0,r_1)}u(\boldsymbol{z})\,\mathrm{d}\boldsymbol{z}+ \int_{B(\boldsymbol{z}_0,r_1)}G_d^{B(\boldsymbol{z}_0,r_1)}(\boldsymbol{z}_0,\boldsymbol{z})(-g)(\boldsymbol{z})\,\mathrm{d}\boldsymbol{z}. \end{align}\] That justifies the estimate \[\begin{align} \label{eq:woswithsource2} \hat{u}(\boldsymbol{z}_0) = h(\bar{\boldsymbol{z}}_{\tau}) -\sum_{k=1}^{\tau} \mathrm{vol}(B(\boldsymbol{z}_{k-1},r_k))\, G_d^{B(\boldsymbol{z}_{k-1},r_k)}(\boldsymbol{z}_{k-1},\boldsymbol{w}_k)\, g(\boldsymbol{w}_k) \end{align}\tag{2}\] that we average over \(n\) independent replicates using independent vectors \(\boldsymbol{w}_k\sim \mathbf{U}(B(\boldsymbol{z}_{k-1},r_k))\) and \(\boldsymbol{z}_{k}\sim \mathbf{U}(S(\boldsymbol{z}_{k-1},r_k))\) each time.

Our examples have \(d\in\{2,3\}\) and for \(\boldsymbol{w}\in B(\boldsymbol{z},r)\), Appendix A.1 of [12] gives \[G_2^{B(\boldsymbol{z},r)}(\boldsymbol{z},\boldsymbol{w})= \frac{\log(r\Vert\boldsymbol{w}-\boldsymbol{z}\Vert^{-1})}{2\pi} \quad\text{and}\quad G_3^{B(\boldsymbol{z},r)}(\boldsymbol{z},\boldsymbol{w})= \frac{1}{4\pi}\Bigl(\frac{1}{\Vert\boldsymbol{w}-\boldsymbol{z}\Vert}-\frac{1}{r}\Bigr).\] Our \(r\) is their \(R\) and their \(r\) is our \(\Vert \boldsymbol{w}-\boldsymbol{z}\Vert\).

An interesting special case arises when the source term \(g\) takes a constant value \(\nu\). Then instead of sampling \(\boldsymbol{w}_k\), we may use \[\begin{align} \label{eq:wosconstsource} \hat{u}(\boldsymbol{z}_0) &=h(\bar{\boldsymbol{z}}_{\tau}) -\nu\sum_{k=1}^{\tau} \int_{B(\boldsymbol{z}_{k-1},r_k)} G_d^{B(\boldsymbol{z}_{k-1},r_k)}(\boldsymbol{z}_{k-1},\boldsymbol{w})\,\mathrm{d}\boldsymbol{w}\notag\\ &=h(\bar{\boldsymbol{z}}_{\tau}) -\frac{\nu}{2d}\sum_{k=1}^{\tau} r_k^2 \end{align}\tag{3}\] for \(d=2,3\) by equations (27) and (28) of [12], with \(\boldsymbol{z}_{k}\sim \mathbf{U}(S(\boldsymbol{z}_{k-1},r_k))\). This simplifies further when \(u\) vanishes on \(\partial\Omega\) which we see for a dumbbell-shaped domain in Section 8.4.

2.4 Randomized Quasi-Monte Carlo↩︎

We assume familiarity with the basic notions of QMC and RQMC. This background can be found in [13], [14] or [15].

We use RQMC points \(\boldsymbol{x}_1,\dots,\boldsymbol{x}_n\in[0,1]^d\) to approximate \(\mu = \int_{[0,1]^d}f(\boldsymbol{x})\,\mathrm{d}\boldsymbol{x}\) by \((1/n)\sum_{i=1}^nf(\boldsymbol{x}_i)\). The function \(f\) includes the transformation of uniform random variables to the WoS trajectories and sample points \(\boldsymbol{w}\) as well as the evaluations of \(h\) and \(g\). The RQMC points individually satisfy \(\boldsymbol{x}_i\sim\mathbf{U}[0,1]^d\) and that gives \(\mathbb{E}(\hat{\mu})=\mu\). They are more evenly distributed through \([0,1]^d\) than IID Monte Carlo (MC) points would be as quantified by discrepancy measures, which is a reason to believe that they will provide more accuracy than MC points do.

We include several RQMC algorithms from QMCPy [16]. These are scrambled Sobol’ points, scrambled Niederreiter points, scrambled Halton points and randomly shifted lattice rules. The Sobol’ points use the direction numbers from [17] and are given the matrix scramble and digital shift from [18]. The Niederreiter points and the Halton points are similarly given those scrambles. The default lattice rule in QMCPy at the time of our computation is lattice-33002-1024-1048576.9125 from Frances Kuo’s site: https://web.maths.unsw.edu.au/~fkuo/lattice/. That lattice was designed for sample sizes \(n=2^m\) for \(10\leqslant m\leqslant 20\). We got good results for smaller values of \(n\) but anybody getting poor results for \(n<2^{10}\) should not view that as a weakness of the lattice.

The scrambled Sobol’, Niederreiter and Halton points all satisfy \(\mathrm{var}(\hat{\mu}) = o(1/n)\) as \(n\to\infty\) while even an adversarially chosen integrand in \(L_2\) yields \(\mathrm{var}(\hat{\mu})\leqslant\Gamma_d\sigma^2/n\) for finite \(n\) and some gain coefficient \(\Gamma_d<\infty\) where \(\sigma^2/n\) is the variance we would have had using Monte Carlo (MC) instead of QMC. The value of \(\Gamma_d\) grows only logarithmically in \(d\) for the scrambled Halton points [19] while it may grow exponentially in \(d\) for scrambled nets. There is also a strong law of large numbers for scrambled nets such as those of Sobol’ or Niederreiter, when \(f\in L_{1+\varepsilon}\) [20].

We will show numerical results for all these methods. The RQMC rules typically all outperform MC and have very nearly the same performance as each other. There was one example where the Niederreiter points performed notably worse than the other methods.

While lattice rules do not have a bounded gain coefficient \(\Gamma\), we took a special interest in them because WoS integrands are periodic and lattice rules are especially well suited to periodic integrands [21]. The best results for lattice rules also depend on smoothness, but the WoS integrands are not smooth.

We will need some results from [6] on the accuracy of RQMC estimation for integrands \(f(\boldsymbol{x}) =\boldsymbol{1}_\Theta(\boldsymbol{x})\) and the extensions in [7] to \(\boldsymbol{1}_\Theta(\boldsymbol{x})g(\boldsymbol{x})\) where \(\Theta\subset[0,1]^k\) is measurable and \(g\) satisfies a boundary condition. Note that this \(g\) is distinct from our source function. The function \(\boldsymbol{1}_{\Theta}\) is typically not in BVHK [22]. It is of course in \(L_2\) and so RQMC has variance \(o(1/n)\) for it but we will quote some sharper rates. Those results use the notion of Minkowski content which we present in Section 4 along with other geometric definitions.

Theorem 3.5 of [6] shows that if \(g\) is in BVHK and \(\partial\Theta\) has finite \(k-1\) dimensional Minkowski content, then scrambled net integration of \(\boldsymbol{1}_\Theta g\) has a mean squared error of \(O(n^{-1-1/(2k-1)}\log(n)^{2(k-1)/(2k-1)})\). Theorem 4.4 of that same paper shows that the scrambled net integration error of \(\boldsymbol{1}_\Theta\) in that case is \(O(n^{-1-1/k})\). Their upper bounds include leading constants. [7] shows that the rate is \(O(n^{-1-1/k})\) when \[\biggl|\frac{ \partial^{|u|}g(\boldsymbol{z})}{\prod_{j\in u}\partial z_j}\biggr| \leqslant\prod_{j=1}^k\min(z_j,1-z_j)^{-A_j-1_{\{j\in u\}}}\] holds for all non-empty \(u\subseteq\{1,2,\dots,k\}\) and \(\max_jA_j<1/2\). In particular this holds when all the above partial derivatives of \(g\) are bounded.

2.5 The WoS integrand in RQMC↩︎

Each step of the WoS algorithm consumes some number \(s\geqslant 1\) of uniform random variables. Let \(\psi_0:[0,1)^{s_0}\to S(\boldsymbol{0},1)\) with \(\psi_0(\boldsymbol{z})\sim\mathbf{U}(S(\boldsymbol{0},1))\) when \(\boldsymbol{z}\sim\mathbf{U}([0,1)^{s_0})\) and \(\psi_1:[0,1)^{s_1}\to B(\boldsymbol{0},1)\) with \(\psi_1(\boldsymbol{z})\sim\mathbf{U}(B(\boldsymbol{0},1))\) when \(\boldsymbol{z}\sim\mathbf{U}([0,1)^{s_1})\). For harmonic \(u\), we can then use \(s=s_0\) uniform variables to take one step of WoS. The first step is \(\boldsymbol{z}_1 \gets \boldsymbol{z}_0 + r_1\psi_0(\boldsymbol{x})\) for \(\boldsymbol{x}\in[0,1)^{s_0}\) and to allow for \(K\) steps requires \(Ks_0\) uniform variables. We then use an RQMC point set with \(\boldsymbol{x}_i\in[0,1)^{Ks_0}\) for \(i=1,\dots,n\). When the source \(g\) is nonzero, we use \(s=s_0+s_1\) uniform variables to take the step and also sample \(\boldsymbol{w}_k\). It then requires RQMC points \(\boldsymbol{x}_i\in[0,1)^{K(s_0+s_1)}\) to generate the WoS estimate. The RQMC point sets we use are all high dimensional and support very large \(K\) and the WoS convergence is fast enough that most samples have small \(\tau\).

There are standard mappings for \(\psi_0\) and \(\psi_1\). The ones in [23] have \(s_0=d-1\) and \(s_1=d\). Our examples in this paper have \(d=2\) or \(3\).

For \(d=2\), we always take \(\psi_0(x) = \theta(x) := (\cos(2\pi x),\sin(2\pi x))\sim\mathbf{U}(S_2(\boldsymbol{0},1))\) for \(x\sim\mathbf{U}[0,1)\). We use \(\psi_1(\boldsymbol{x})=(\sqrt{x_1}\cos(2\pi x_2),\sqrt{x_1}\sin(2\pi x_2))\) for \(\boldsymbol{x}\sim\mathbf{U}([0,1)^2)\) to sample \(\mathbf{U}(B_2(\boldsymbol{0},1))\).

For \(d=3\), we sample \(\boldsymbol{x}\sim\mathbf{U}([0,1)^2)\), define a latitude \(\lambda(\boldsymbol{x})=2x_1-1\), a longitude \(\theta(\boldsymbol{x})= 2\pi x_2\) and use the hat box transformation \[\begin{align} \label{eq:hatbox} \psi_0(\boldsymbol{x})=\bigl(\sqrt{1-\lambda^2(\boldsymbol{x})}\cos(\theta(\boldsymbol{x})),\,\sqrt{1-\lambda^2(\boldsymbol{x})}\sin(\theta(\boldsymbol{x})),\,\lambda(\boldsymbol{x})\bigr) \end{align}\tag{4}\] to get \(\psi_0(\boldsymbol{x})\sim\mathbf{U}(S_3(\boldsymbol{0},1))\). None of our examples with \(d=3\) had a nonzero source function, but for that case we would use \(\boldsymbol{w}=\psi_1(\boldsymbol{x}) = x_1^{1/3}\psi_0((x_2,x_3))\) for \(x\sim \mathbf{U}([0,1)^3)\).

The case of a harmonic function on \(\Omega\subset\mathbb{R}^2\) is especially interesting as \(\boldsymbol{z}_k\) is now a periodic function of the point \(\boldsymbol{x}\in[0,1)^k\) from which the first \(k\) WoS steps are generated. As mentioned earlier, periodic functions are especially favorable for lattice sampling.

Algorithm [alg:wos] has pseudo-code to compute the WoS estimate of \(u(\boldsymbol{z}_0)\) by RQMC when the source term \(g=0\). For a nonzero source term, the algorithm should be extended to include a second mapping \(\psi_1\) to generate the point \(\boldsymbol{w}_k\) at step \(k\) from \(s_1\) uniform inputs, and the RQMC points should have dimension \(K(s_0+s_1)\), with \(K\) being the maximum number of steps the WoS can take per trajectory.

3 The gasket example↩︎

We illustrate RQMC-WoS with an example from [9]. The authors consider the thermal conduction in a cylinder head gasket illustrated in Figure 1. In this problem \(u(\boldsymbol{z})\) is the temperature of the gasket at \(\boldsymbol{z}\in\Omega\). The purpose of the head gasket is to avoid any leakage of oil or coolant from the combustion engine into the cylinders where combustion occurs. Leaks are more likely to occur where the gasket is hotter.

They are interested in approximating the solution to the BVP with \(\Delta u(\boldsymbol{z})=0\) for \(\boldsymbol{z}\in\Omega^\circ\) and boundary values \[u(\boldsymbol{z}) = \sum_{r\in\{\mathrm{coolant},\,\mathrm{outer},\,\mathrm{oil},\,\mathrm{oil\;return},\,\mathrm{bore}\}} T_r\,\boldsymbol{1}_r(\boldsymbol{z}) \quad \textrm{ for } \boldsymbol{z}\in\partial\Omega\] where \(T_r\) is a constant temperature (in degrees Celsius) specific to boundary component \(r\) and \(\boldsymbol{1}_r(\boldsymbol{z})\) is the indicator of component \(r\). We computed the WoS estimator of \(u\) at point \(\boldsymbol{z}_0 =(0.240999,0.3)\) located above the center of the third borehole and roughly half way to the edge of the gasket as indicated in Figure 1. This point is sufficiently interior to avoid immediate absorption, yet close enough to several boundary components that typical WoS trajectories require multiple steps and the exit location is nontrivial. For each sample size, we consider Monte Carlo sampling and the four RQMC algorithms described in Section 2.4, and run \(100\) independent replicates of each method. The variance curves are shown in Figure 2. The reference curve was fit by a regression of log variance on \(\log(n)\), pooling data from all four RQMC methods for \(n\geqslant 2^7\). We use that definition for all the reference curves in this paper. The regression coefficients for individual RQMC methods are given in Table 1 of Section 8.6 for five numerical examples.

In this example we see that using RQMC points brings a better convergence rate for the variance, with a slope near \(-1.11\) in a log-log plot. As mentioned earlier, the different RQMC samplers have very nearly equal performance. Here, for \(n=2^{17}\), the RQMC points bring a roughly \(5\)-fold variance reduction, compared to standard MC. The performance of RQMC varies spatially in this domain. We describe some other results in Section 8.

4 Geometric quantities↩︎

It seems intuitively clear that the sets of RQMC points that lead to the outcomes of interest in WoS problems should be quite regular and not, for instance, of a fractal nature. However, finding sufficient conditions to justify this belief runs into some complications. Here, we present some background results on Hausdorff measure, Minkowski content, rectifiability of curves, domains in \(\mathbb{R}^2\) and analyticity of distances to plane curves.

4.1 Hausdorff measure and Minkowski content↩︎

The way we will control Minkowski content is via Theorem 3.2.39 of [24] quoted in Section 4.2 which gives conditions that make it equal to Hausdorff measure. Then we will apply theorems that control Hausdorff measure.

We use the definition of Hausdorff measure from [25] because it is more easily stated than some others. An arbitrary set \(C\subset\mathbb{R}^k\) has diameter \(\mathrm{diam}(C)=\sup_{\boldsymbol{x},\tilde{\boldsymbol{x}}\in C}\Vert\boldsymbol{x}-\tilde{\boldsymbol{x}}\Vert\), with \(\mathrm{diam}(\varnothing)=0\) by convention. Then for \(\Theta\subset\mathbb{R}^k\), \(0\leqslant m<\infty\) and \(0<\delta\leqslant\infty\) we let \[\mathcal{H}_{m,\delta}(\Theta) = \frac{\pi^{m/2}}{\Gamma(m/2+1)}\times \inf\Biggl\{\sum_{\ell=1}^\infty \Bigl(\frac{\mathrm{diam}(C_\ell)}{2}\Bigr)^m \biggm| \Theta\subset\bigcup_{\ell=1}^\infty C_\ell, \mathrm{diam}(C_\ell)\leqslant\delta \Biggr\}\] where \(\Gamma(\cdot)\) is the Gamma function. When \(m\) is an integer the leading factor above is \(\mathrm{vol}(B_m(\boldsymbol{0},1))\). The \(m\)-dimensional Hausdorff measure of \(\Theta\) is \[\begin{align} \label{eq:defhaus} \mathcal{H}_m(\Theta) = \lim_{\delta\to0}\mathcal{H}_{m,\delta}(\Theta) =\sup_{\delta>0}\mathcal{H}_{m,\delta}(\Theta). \end{align}\tag{5}\] We only need \(\mathcal{H}_m\) for integer \(m\). For \(m=0,1,2,3\), we get a measure that matches the usual notions of cardinality, length, area and volume. The thickness \(\alpha\) in [11] is defined using a closely related quantity called Hausdorff content.

We follow Chapter 3.3 of [26] to define Minkowski content. For \(m\in\{0,1,\dots,k\}\), the upper \(m\)-dimensional Minkowski content of non-empty \(\Theta\subset\mathbb{R}^k\) is \[\overline{\mathbb{M}}_{m}(\Theta) = \limsup_{\varepsilon\to0^+} \frac{\mathrm{vol}(\{\boldsymbol{x}\mid \mathrm{dist}(\boldsymbol{x},\Theta)<\varepsilon\})}{\mathrm{vol}(B_{k-m}(\boldsymbol{0},\varepsilon))}.\] The corresponding lower \(m\)-dimensional Minkowski content is \[\underline{\mathbb{M}}_{m}(\Theta) = \liminf_{\varepsilon\to0^+} \frac{\mathrm{vol}(\{\boldsymbol{x}\mid \mathrm{dist}(\boldsymbol{x},\Theta)<\varepsilon\})}{\mathrm{vol}(B_{k-m}(\boldsymbol{0},\varepsilon))}.\] When \(\underline{\mathbb{M}}_m(\Theta)=\overline{\mathbb{M}}_m(\Theta)\), we let \(\mathbb{M}_m(\Theta)\) denote their common value. If \(\mathbb{M}_m(\Theta)<\infty\), then we say that \(\Theta\) has \(m\) dimensional Minkowski content.

The way Minkowski content is used in RQMC is to partition the set \([0,1)^k\) into \(N\) sets \(T_1,\dots,T_N\). We write our integrand as \(f(\boldsymbol{x}) = 1\{\boldsymbol{x}\in\Theta\} = \sum_{j=1}^Nf_j(\boldsymbol{x})\) where \(f_j(\boldsymbol{x}) = 1\{\boldsymbol{x}\in \Theta\cap T_j\}\). Then our estimate of \(\mathrm{vol}(\Theta)\) is \[\widehat\mathrm{vol}(\Theta) = \frac{1}{n} \sum_{i=1}^n\sum_{j=1}^Nf_j(\boldsymbol{x}_i)\] for \(n\) sample points \(\boldsymbol{x}_i\in[0,1)^k\).

For scrambled \((t,m,k)\)-nets in base \(b\) it is convenient to take \(T_j\) to be elementary intervals of volume \(b^{t-m}\) with side lengths as nearly equal as possible. For the commonly used Sobol’ points with \(b=2\), these are dyadic hyperrectangles. Then \((1/n)\sum_{i=1}^nf_j(\boldsymbol{x}_i) = \int_{[0,1)^k}f_j(\boldsymbol{x})\) if either \(T_j\subset\Theta\) or \(T_j\subset\Theta^c\). Let \[\mathcal{B}=\bigl\{ j\in\{1,2,\dots,N\}\mid T_j\cap \Theta\ne\varnothing\;\&\;T_j\cap \Theta^c\ne\varnothing\bigr\}\] be the set of ‘boundary’ elementary intervals. Then \(\widehat\mathrm{vol}(\Theta) -\mathrm{vol}(\Theta)\) is the error in estimating the integral of \[f_\mathrm{bdy}(\boldsymbol{x})=\sum_{j\in\mathcal{B}}f_j(\boldsymbol{x}).\] The RQMC variance is no more than some gain factor \(\Gamma<\infty\) times the Monte Carlo variance. For \(n\geqslant 2\), that Monte Carlo variance is no more than \(|\mathcal{B}|/n^2\) because \(f_j(\boldsymbol{x})\in\{0,1\}\) takes the value \(1\) with probability at most \(1/n\) giving it variance at most \(1/n\). Then the RQMC variance is at most \(\Gamma|\mathcal{B}|/n^2\). If \(\overline{\mathbb{M}}_{k-1}(\partial\Theta)<\infty\), then \(|\mathcal{B}| = O(N^{(k-1)/k})\) and for digital nets \(N = n/b^t\). Then the RQMC variance is \(O(n^{-2+(k-1)/k})=O(n^{-1-1/k})\).

The RQMC argument only needs a finite upper Minkowski content for \(\partial\Theta\). Because the RQMC points are uniformly distributed we could reduce the boundary set to \[\begin{align} \tilde{\mathcal{B}}=\bigl\{ j\in\{1,2,\dots,N\}\mid 0<\mathrm{vol}(T_j\cap\Theta)<\mathrm{vol}(T_j)\bigr\}. \end{align}\] We are not aware of any proofs in the literature that use \(\tilde{\mathcal{B}}\) instead of \(\mathcal{B}\).

4.2 Rectifiability from Federer (1969)↩︎

Here are the definitions from Section 3.2.14 of [24]. He references a measure \(\phi\) that is commonly taken to be \(m\) dimensional Hausdorff measure. In each definition \(m\) is a positive integer and \(E\) is a subset of a metric space that we take to be \(\mathbb{R}^d\) for some \(d\geqslant 1\). For \(\phi\) to measure \(\mathbb{R}^d\) in Federer’s sense, it must be a countably subadditive function on all subsets of \(\mathbb{R}^d\), what we now call an outer measure.

We will need to make a somewhat subtle use of Definition 1 part (1). The bounded subset of \(\mathbb{R}^m\) that we use might have strictly lower dimension than \(m\). We will use the fact that any subset \(S\) of an \(m\) rectifiable set \(E\) is also \(m\) rectifiable, even for a set \(S\) with a fractal appearance. This holds because we can write \(S = f(f^{-1}(S)\cap B)\) where \(B\) is bounded and \(E=f(B)\). We need the following elementary result on finite unions of rectifiable sets. We suspect this is well known, but we could not find a reference.

The conclusion to Lemma 1 does not hold for countable unions. The next theorem is one of the few ways to get a conclusion about finite Minkowski content.

Federer’s definition of \(m\) rectifiability is stricter than we see in some more recent works. For example, \(m\)-rectifiability of Definition 15.3 of [27] is \((\mathcal{H}_m,m)\) rectifiability of Definition 1. That definition of rectifiability is easier for us to establish in our WoS problem but it is not strong enough to use in Theorem 1.

4.3 The domain, its boundary and distance functions↩︎

Our study of the WoS algorithm requires smoothness of \(z\mapsto \mathrm{dist}(\boldsymbol{z},\partial\Omega)\) almost everywhere in \(\Omega^\circ\). The required smoothness fails to hold at points \(\boldsymbol{z}\) whose projection onto \(\partial\Omega\) is not unique. We discuss such medial points below as well as some points that project onto places where two smooth subcurves of \(\partial\Omega\) meet. Under our conditions \(\mathrm{dist}(\boldsymbol{z},\partial\Omega)\) is analytic everywhere except on a finite union of analytic curves. The results we need to study the smoothness of the function \(\mathrm{dist}(\cdot,\partial\Omega)\) are scattered over different publications with multiple definitions of the necessary concepts.

We use regularity conditions on \(\Omega\) taken from [8]. First, \(\Omega\) is the closure of a connected and bounded open subset of \(\mathbb{R}^2\). The boundary \(\partial\Omega\) is the union of finitely many disjoint simple closed curves. Those are described as embeddings of the unit circle into \(\mathbb{R}^2\). Such a curve is closed, continuous and does not intersect itself. It does not have to be diffeomorphic to the unit circle like the embeddings from differential topology [28]. For example, it can have corners. When there are \(G+1\) of these curves, then the set \(\Omega\) has one outer boundary curve and \(G\) ‘holes’ cut out of it by inner boundary curves. It is then said to have genus \(G\). The gasket example has \(G=50\). Many problems have genus \(0\).

Each of the \(G+1\) boundary curves is the union of a finite number of pieces we call arcs. Each of those arcs is a real analytic curve defined on a closed interval \([a,b]\). Those in turn are the restrictions of a real analytic curve over an open interval such as \((a-\varepsilon,b+\varepsilon)\) to \([a,b]\).

[8] have to exclude the possibility that \(\partial\Omega\) is actually a circle. That case is easy for our WoS theory and can also be sampled without WoS by one dimensional integration of \(h\) times a Poisson kernel over \(\partial\Omega\), so this exclusion does not cause us difficulty.

Analyticity may seem like an extremely strong assumption, but even \(C^\infty\) curves can show pathologies. See Figures 1 and 2 of [8] where sets with \(C^\infty\) boundaries have medial axes consisting of a countably infinite number of curves with infinite total length. They also remark that many applications define boundaries via non-uniform rational B-splines (NURBS) which satisfy their conditions, along with more commonly considered line segments and circular arcs.

For Dirichlet boundary conditions and a harmonic function \(u\), we will examine the case where \(u(\boldsymbol{z})\) takes the value \(0\) on all of \(\partial\Omega\) except for one arc \(\mathcal{A}\). Such cases form a basis for more general problems. We define the distances \[\begin{align} r_{\partial\Omega}(\boldsymbol{z}) &= \mathrm{dist}(\boldsymbol{z},\partial\Omega) = \min_{\tilde{\boldsymbol{z}}\in\partial\Omega}\Vert\boldsymbol{z}-\tilde{\boldsymbol{z}}\Vert,\quad\text{and}\\ r_\mathcal{A}(\boldsymbol{z}) &= \mathrm{dist}(\boldsymbol{z},\mathcal{A}) = \min_{\tilde{\boldsymbol{z}}\in\mathcal{A}}\Vert\boldsymbol{z}-\tilde{\boldsymbol{z}}\Vert. \end{align}\] Both of these are minima not just infima because they are distances to closed sets. These functions are both \(1\)-Lipschitz from the triangle inequality.

Here is the definition of a geometric graph from [8].

Next, we consider smoothness of the distance functions. For \(k\geqslant 2\), the distance to a \(C^k\) curve \(\mathcal{C}\) is \(C^k\) over \(U\setminus\mathcal{C}\) where \(U\) is a neighborhood of \(\mathcal{C}\) [29]. Our regularity condition requires greater smoothness and requires it over a larger set than just locally near our curves. We use the following definitions taken from Definition 1 of [30] to give a result on global smoothness.

There can be points in \(\mathrm{unpp}(M)\setminus\mathcal{O}(M)\). These are non-medial points to which some sequence of medial points converges. Example 1.2 of [31] has \(M = (\mathbb{R}\times \{-1,1\}) \setminus \{(0,1)\}\). Then as \(n\to\infty\), the medial points \((1/n,0)\) converge to the non-medial point \((0,0)\).

[8] define the medial axis of a CCM boundary \(\Omega\) in a way that is similar to \((\mathrm{unpp}(\partial\Omega))^c\) but with a few differences.

The medial axis \(\mathrm{MA}(\Omega)\) in Definition 6 differs from the set of points in Definition 5 with non-unique projections onto \(\Omega\). Let \(\widetilde{\mathrm{MA}}(\Omega)=(\mathrm{unpp}(\Omega))^c\) be that definition of medial points from [30]. Then \(\widetilde{\mathrm{MA}}(\Omega)\setminus\mathrm{MA}(\Omega)\) is made up of points in the exterior of \(\Omega\) with non-unique projections. These do not affect the WoS algorithm running inside \(\Omega^\circ\). The set \(\mathrm{MA}(\Omega)\setminus \widetilde{\mathrm{MA}}(\Omega)\) contains corner points of \(\partial\Omega\) where a disk of radius \(0\) is maximal but only touches \(\partial\Omega\) at one point. For example, if \(\Omega\) is a rectangle then \(\mathrm{MA}(\Omega)\) has the diagonals including vertices while \(\widetilde{\mathrm{MA}}(\Omega)\) excludes the vertices. Theorem 6.2 of [8] shows that a CCM domain has only finitely many of these points.

The most important result from [8] for our present purposes is the following theorem.

4.4 Distance to a piece of \(\partial\Omega\)↩︎

The distance to an arc \(\mathcal{A}\) that is just one analytic piece of one of the closed curves bounding \(\Omega\) is much more complicated than the combination of [31] and [8]. Such a set \(\mathcal{A}\) has endpoints that are not allowed by Theorem 2 and it is not the boundary of a compact set, which Theorem 3 requires. We did not see any treatment of this case in the literature that is as comprehensive as the study of CCM domains in [8].

Here we study \(r_\mathcal{A}\) for some special cases and then sketch why we think its analyticity away from a geometric graph should hold more generally.

In the special cases we consider next, \(\mathcal{A}\) is either a line segment or a portion of a circular arc. We also assume that any straight line in \(\mathbb{R}^2\) that intersects \(\Omega\) does so in a finite set of disjoint closed intervals. Then analyticity holds for the distance from these curves, apart from points in a finite union of geometric graphs. Many useful domains \(\Omega\) can be constructed using line segments and circular arcs as pieces of \(\partial\Omega\). This includes all the two dimensional examples we consider in this paper.

4.5 \(\mathcal{A}\) is a line segment↩︎

Suppose that \(\mathcal{A}\) is a finite line segment in \(\mathbb{R}^2\). Without loss of generality, \(\mathcal{A}= [0,1]\times\{0\}\). Then \[r_\mathcal{A}(\boldsymbol{z}) = \bigl(\mathrm{dist}(z_1,[0,1])^2+ z_2^2\bigr)^{1/2}\] which is analytic everywhere except \(\mathcal{A}\cup (\{0,1\}\times\mathbb{R})\). This exceptional set is the union of a line segment and two lines. Its intersection with the convex hull of \(\Omega\) is a geometric graph. The holes in \(\Omega\) may split that geometric graph into a finite set of smaller geometric graphs but under our assumption they cannot yield an infinite set of geometric graphs.

4.6 \(\mathcal{A}\) is part of circular arc↩︎

Suppose instead that \(\mathcal{A}\) is a portion of a circular arc. Without loss of generality \(\mathcal{A}= \{ \phi(t) \mid |t|\leqslant c\}\) for \(\phi(t) = (\cos(t),\sin(t))\) and \(0<c<\pi\). Write \(\boldsymbol{z}\) in polar coordinates with radius \(r(\boldsymbol{z})\geqslant 0\) and angle \(\theta(\boldsymbol{z})\in[-\pi,\pi)\). Define regions \[\begin{align} S_{\mathrm{In}} & =\{\boldsymbol{z}\in\mathbb{R}^2\mid |\theta(\boldsymbol{z})|<c,\;0<r(\boldsymbol{z})<1\},\\ S_{\mathrm{Out}} & =\{\boldsymbol{z}\in\mathbb{R}^2\mid |\theta(\boldsymbol{z})|<c,\;r(\boldsymbol{z})>1\},\\ S_{\mathrm{Up}} &=\{\boldsymbol{z}\in\mathbb{R}^2\mid c<\theta(\boldsymbol{z})<\pi\}, \quad\text{and}\\ S_{\mathrm{Down}}&=\{\boldsymbol{z}\in\mathbb{R}^2\mid -\pi <\theta(\boldsymbol{z})<-c\}, \end{align}\] depicted in Figure 3. A point \(\boldsymbol{z}\) in \(S_{\mathrm{In}}\) or \(S_{\mathrm{Out}}\) projects onto \(\phi(\theta(\boldsymbol{z}))\) and has \(r_\mathcal{A}(\boldsymbol{z})=|r(\boldsymbol{z})-1|\) which is analytic on those two sets. If \(\boldsymbol{z}\in S_{\mathrm{Up}}\) then \(\boldsymbol{z}\) projects to \(\phi(c)\) and \(r_{\mathcal{A}}(\boldsymbol{z})=\Vert\boldsymbol{z}-\phi(c)\Vert\) is analytic on \(S_{\mathrm{Up}}\). Similarly, \(r_\mathcal{A}\) is analytic on \(S_{\mathrm{Down}}\). The exceptional sets where \(r_\mathcal{A}\) is not analytic are \(\mathcal{A}\), the point \((0,0)\) which projects non-uniquely onto \(\mathcal{A}\), the ray \(R_0=\{(-t,0)\mid 0< t<\infty\}\) of points with \(r_{\phi(c)}(\boldsymbol{z})=r_{\phi(-c)}(\boldsymbol{z})=\mathrm{dist}(\boldsymbol{z},\mathcal{A})\), and the two rays \(R_a=\{r(\cos(a),\sin(a))\mid 0< r<\infty\}\) for \(a=\pm c\) which each separate a region projecting to \(\mathcal{A}^\circ = \phi( (-c,c))\) from a region projecting to an endpoint of \(\mathcal{A}\).

As in the case of a line segment, \(r_\mathcal{A}\) is analytic on the convex hull of \(\Omega\) apart from \(\mathcal{A}\) itself and three line segments through the origin. Then analyticity holds apart from a finite union of geometric graphs and the problem is WoS regular.

4.7 More general \(\mathcal{A}\)↩︎

The relevant literature for this case considers medial axes, point/curve bisectors and self bisectors of curves. Lacking results like those in [8], we sketch our reasons for believing that the distance to an analytic arc \(\mathcal{A}\) will be analytic apart from a geometric graph quite generally and not just for special cases like line segments and circular arcs. We draw on a result in [32] which uses parameterized arcs. We write our real analytic arc as \(\mathcal{A}= \{\phi(t) \mid 0\leqslant t\leqslant 1\}\) and we assume that \(\phi'(t)\ne\boldsymbol{0}\) at any \(t\in[0,1]\).

With this parameterization, we define an internal curve \(\mathcal{A}^\circ = \phi( (0,1))\) and endpoints \(\boldsymbol{a}=\phi(0)\) and \(\boldsymbol{b}=\phi(1)\). Then \[\begin{align} \label{eq:componentdistances} r_\mathcal{A}(\boldsymbol{z}) = \min\bigl( r_{\boldsymbol{a}}(\boldsymbol{z}), r_{\mathcal{A}^\circ}(\boldsymbol{z}),r_{\boldsymbol{b}}(\boldsymbol{z})\bigr) = \min\bigl( \Vert \boldsymbol{z}-\boldsymbol{a}\Vert, r_{\mathcal{A}^\circ}(\boldsymbol{z}),\Vert \boldsymbol{z}-\boldsymbol{b}\Vert\bigr). \end{align}\tag{6}\] The endpoint distances \(\Vert\boldsymbol{z}-\boldsymbol{a}\Vert\) and \(\Vert\boldsymbol{z}-\boldsymbol{b}\Vert\) are analytic at all but \(\boldsymbol{a}\) and \(\boldsymbol{b}\) respectively. The function \(r_\mathcal{A}\) can fail to be analytic at points where \(r_{\mathcal{A}^\circ}\) is not analytic or at points where two or more of the distances in 6 are equal. We consider those equalities first and then describe how smoothness of \(r_{\mathcal{A}^\circ}\) is the gap that remains.

The set of points \(\boldsymbol{z}\) where \(\Vert\boldsymbol{z}-\boldsymbol{a}\Vert=\Vert\boldsymbol{z}-\boldsymbol{b}\Vert\) is a line \(\ell_{\boldsymbol{a}\boldsymbol{b}}\) and so its intersection with \(\Omega\) is contained within a finite union of line segments. Hence it is a finite union of geometric graphs.

We can never have \(\mathrm{dist}(\boldsymbol{z},\boldsymbol{a})<\mathrm{dist}(\boldsymbol{z},\mathcal{A}^\circ)\). Those distances are equal on the set \[E_{\boldsymbol{a}} = \bigcap_{0<t<1}\bigl\{\boldsymbol{z}\mid \Vert\boldsymbol{z}-\boldsymbol{a}\Vert <\Vert \boldsymbol{z}-\phi(t)\Vert\bigr\}\] and there is an analogous set \(E_{\boldsymbol{b}}\) at the other endpoint of \(\mathcal{A}\). The set \(E_{\boldsymbol{a}}\) is an intersection of open half-spaces so it is convex. For any \(\boldsymbol{z}\in E_{\boldsymbol{a}}^{\circ}\setminus\ell_{\boldsymbol{a}\boldsymbol{b}}\), there is an open set \(U\) with \(\boldsymbol{z}\in U\subset E_{\boldsymbol{a}}\setminus\ell_{\boldsymbol{a}\boldsymbol{b}}\) on which \(\mathrm{dist}(\boldsymbol{z},\mathcal{A})=\mathrm{dist}(\boldsymbol{z},\boldsymbol{a})=\Vert\boldsymbol{z}-\boldsymbol{a}\Vert\) which is analytic. That argument does not work on \(\partial E_{\boldsymbol{a}}\) which [32] discuss.

The set of points where \(\Vert\boldsymbol{z}-\boldsymbol{a}\Vert =\mathrm{dist}(\boldsymbol{z},\mathcal{A})\) is part of the point/curve bisector between \(\mathcal{A}\) and \(\boldsymbol{a}\). Point/curve bisectors where the point is part of the curve are a degenerate case and Section 3.2 of [32] considers the bisector of a curve and one of its endpoints. They say that the bisector is a subset of three curves, defined by their equations (a), (b) and (c). Curve (c) is the normal line to \(\mathcal{A}\) at \(\boldsymbol{a}\). Curve (b) is the linear end-tangent extension to \(\mathcal{A}\) at the endpoints. Curve (a) is defined using the unit normal curve \(\nu(t)\) where \(\nu(t)\) is \(\phi'(t)/\Vert\phi'(t)\Vert\) rotated through 90 degrees. Then curve (a) has the expression \[a(t)= \phi(t) + \nu(t)\frac{|\boldsymbol{a}-\phi(t)|^2}{2(\boldsymbol{a}-\phi(t))^\mathsf{T}\nu(t)}\] and we interpret \(a(0)\) as \(\boldsymbol{a}\). Now if \(\phi\) and \(\nu\) are both analytic on \(\mathcal{A}^\circ\), then so is \(a\) on intervals where the denominator does not vanish. As \(t\) approaches an interior point where the denominator vanishes, \(\Vert a(t)\Vert\) diverges to infinity, thereby leaving the domain \(\Omega\). As a result we expect the curve \(a(t)\) to intersect \(\Omega\) only in a finite union of geometric graphs.

The gap is that we did not find conditions in the literature to ensure that \(r_{\mathcal{A}^\circ}\) is analytic apart from a finite union of geometric graphs. We would need the complement of \((\mathrm{unpp}(\mathcal{A}^\circ))^\circ\) to be a finite union of geometric graphs. Theorem 3.9 of [33] is similar to what we need but we could not determine whether it applies to \(\mathcal{A}^\circ\) nor whether the resulting one dimensional simplicial complex it describes would qualify as a finite union of geometric graphs.

5 Differential geometry theorems↩︎

We need some results from differential geometry. These are not commonly used in the RQMC literature so we state them here. We will not need the most general versions of these theorems, just the ones with statements that match our uses.

We need the notions of regular and critical values of a function as well as regular and critical points. For a smooth function \(f:M\to N\) between manifolds the point \(\boldsymbol{x}\in M\) is a critical point of \(f\) if the Jacobian of \(f\) at \(\boldsymbol{x}\) has less than full rank. Then \(\boldsymbol{y}= f(\boldsymbol{x})\) is a critical value of \(f\). If \(\boldsymbol{x}\) is not a critical point of \(f\) then it is a regular point of \(f\). If every \(\boldsymbol{x}\) with \(f(\boldsymbol{x})=\boldsymbol{y}\) is a regular point of \(f\) then \(\boldsymbol{y}\) is a regular value of \(f\). A regular level set is the level set of a regular value.

Smooth means \(C^\infty\) throughout [34] and this is enough for our purposes. The level set of \(f:M\to N\) has dimension equal to the dimension of \(M\) minus the dimension of \(N\), which in our case will be the dimension of \(M\) minus \(1\). [28] call their version of this result, the pre-image theorem. Because the level set is an embedded manifold, it is smooth [34]. We will not reference the properness of this embedding.

6 Minkowski content for WoS↩︎

For \(d=2\) our RQMC algorithm constructs \(\boldsymbol{z}_k\) from \(\boldsymbol{z}_0\) and a point \(\boldsymbol{x}\in[0,1)^k\). Here we consider the set of input points in \([0,1)^k\) that cause \(\mathrm{dist}(\boldsymbol{z}_k,\partial\Omega)<\varepsilon\) along with \(u(\bar\boldsymbol{z}_k)=1\). We give conditions under which that set has \(\overline{\mathbb{M}}_{k-1}<\infty\) for almost all \(\varepsilon\). A separate analysis in Section 7 accounts for the possibility that the walk may have terminated at \(k'<k\) steps.

We assume that \(\Omega\subset \mathbb{R}^2\) is a CCM domain. Then \(\partial\Omega\) is a finite union of closed analytic arcs one of which is \(\mathcal{A}\). The target function \(u:\partial\Omega\to\{0,1\}\) takes the value \(1\) on \(\mathcal{A}\) and is \(0\) in \(\partial\Omega\setminus\mathcal{A}\). We abbreviate \(r_{\mathcal{A}}\) to \(r_1\), \(r_{\partial\Omega\setminus\mathcal{A}}\) to \(r_0\) and \(r_{\partial\Omega}\) to \(r\).

Now let \[\begin{align} \label{eq:defcp} \mathcal{P}= \bigl\{\boldsymbol{z}\in\Omega^\circ\mid \text{r_1 is not analytic}\bigr\} \bigcup\, \bigl\{\boldsymbol{z}\in\Omega^\circ\mid \text{r is not analytic}\bigr\} \end{align}\tag{7}\] be the collection of points in \(\Omega^\circ\) where at least one of our distance functions fails to be analytic. For a WoS regular problem, \(\mathcal{P}\) is a finite union of geometric graphs.

The WoS problems we consider start at \(\boldsymbol{z}_0\in\Omega\) with \(r(\boldsymbol{z}_0)>\varepsilon\). A point \(\boldsymbol{x}\in[0,1)^k\) generates steps \(\boldsymbol{z}_j(\boldsymbol{x})\) for \(j=1,\dots,k\) via \[\boldsymbol{z}_{j}(\boldsymbol{x}) = \boldsymbol{z}_{j-1}(\boldsymbol{x}) + r(\boldsymbol{z}_{j-1}(\boldsymbol{x})) \theta(x_j),\quad\text{for \theta(x)= (\cos(2\pi x),\sin(2\pi x))^\mathsf{T}}.\] By convention, \(\boldsymbol{z}_0(\boldsymbol{x})=\boldsymbol{z}_0\) a constant function on \([0,1)^k\). The function \(\theta\) above is Lipschitz continuous with a constant of \(2\pi\) equal to the norm of its derivative.

We will use the set \[\mathcal{E}_k=\bigcup_{j=1}^k \boldsymbol{z}_j^{-1}(\mathcal{P})\subset\mathbb{T}^k\] to represent all walks that ever hit a problematic point in their first \(k\) steps. Under RQMC sampling, \(\Pr(\boldsymbol{x}\in\mathcal{E}_k)=0\).

We will use the following ‘chamber regularity’ assumption on \(\mathcal{E}_k\). This assumption is illustrated by Figure 4 with a context described in the proof of Theorem 10.

This definition does not allow for two of the manifolds in \(\mathcal{E}_k\) to intersect infinitely often within any neighborhood of \(\boldsymbol{p}\). That conclusion might possibly follow from arguments using analyticity similar to those used in [8] and we think that WoS regularity is almost enough to give chamber regularity, but exploring these points is outside the scope of this article.

It is very convenient that \(\boldsymbol{z}_k\) is a periodic function on \(\mathbb{R}^k\). We can therefore choose the domain to be the flat torus \(\mathbb{T}^k=\mathbb{R}^k/\mathbb{Z}^k\). This is a smooth manifold without boundary. We will need to apply the implicit function theorem to a function defined on a flat torus. That theorem is stated for functions on Euclidean space. It concerns local properties of the function. A small neighborhood \(U\) of \(\boldsymbol{x}\in\mathbb{T}^k\) can be identified with a small neighborhood \(\tilde{U}\) around a representative of \(\boldsymbol{x}\) in Euclidean space. We can then view a function on \(U\subset\mathbb{T}^k\) as a function on \(\tilde{U}\subset [-1,2)^k\subset\mathbb{R}^k\).

The \(k\)-step distance to \(\mathcal{A}\) is \[r_{k,1}(\boldsymbol{x}) = \mathrm{dist}(\boldsymbol{z}_k(\boldsymbol{x}),\mathcal{A}) =r_1(\boldsymbol{z}_k(\boldsymbol{x})).\] For RQMC we want \(\overline{\mathbb{M}}_{k-1}(\{\boldsymbol{x}\in[0,1)^k\mid r_{k,1}(\boldsymbol{x})=\varepsilon\})<\infty\).

7 Randomized quasi-Monte Carlo for WoS↩︎

Here we use the prior results on Minkowski content to understand how RQMC works on WoS problems. For integers \(k\geqslant 1\) we can use RQMC points to estimate \[\begin{align} \eta_{k,\ell} &= \int_{[0,1)^k} 1\{ r_{k,\ell}(\boldsymbol{x}) < \varepsilon\}\,\mathrm{d}\boldsymbol{x}\quad\text{by}\quad \hat{\eta}_{k,\ell} = \frac{1}{n}\sum_{i=1}^n 1\{ r_{k,\ell}(\boldsymbol{x}_i) < \varepsilon\} \end{align}\] for \(\ell=0,1\). By the finite \((k-1)\)-dimensional Minkowski content of the level sets established in the previous subsection, we can apply Theorem 4.4 of [6] to show for almost all \(\varepsilon\) that \(\mathrm{var}( \hat{\eta}_{k,\ell} ) = O( n^{-1-1/k})\). We may not know in practice whether our \(\varepsilon\) value is a favorable one. By the same token, it would be an odd coincidence for the geometry of \(\mathcal{A}\) and \(\partial\Omega\) to make \(\varepsilon\) an unfavorable value. Similarly, the results from [7] support this rate for \(h(\bar{\boldsymbol{z}}_k(\boldsymbol{x}))1\{r_{k,\ell}(\boldsymbol{x})<\varepsilon\}\) when \(h(\bar{\boldsymbol{z}}_k(\boldsymbol{x}))\) is sufficiently smooth.

Next, we estimate \[\mu_{k,1} =\int_{[0,1)^k} 1\{ r_{k,1}(\boldsymbol{x}) < \varepsilon\} \times\prod_{j=1}^{k-1} 1\{ r_{j,1}(\boldsymbol{x}) \geqslant\varepsilon\} \times\prod_{j=1}^{k-1} 1\{ r_{j,0}(\boldsymbol{x}) \geqslant\varepsilon\} \,\mathrm{d}\boldsymbol{x}\] by \[\hat{\mu}_{k,1} = \frac{1}{n}\sum_{i=1}^n 1\{ r_{k,1}(\boldsymbol{x}_i) < \varepsilon\} \times\prod_{j=1}^{k-1} 1\{ r_{j,1}(\boldsymbol{x}_i) \geqslant\varepsilon\} \times\prod_{j=1}^{k-1} 1\{ r_{j,0}(\boldsymbol{x}_i) \geqslant\varepsilon\}.\] The integrand in \(\mu_{k,1}\) is the indicator of the intersection of \(2k-1\) sets. The boundary of the intersection is contained in the union of those \(2k-1\) boundaries each of which has \(k-1\) dimensional Minkowski content for almost all \(\varepsilon\). As a result, this set either has \(k-1\) dimensional Minkowski content or has upper \(k-1\) dimensional Minkowski content of zero and so we get the same rate.

Our estimate of \(\mu_1 = \mathbb{E}( u(\boldsymbol{z}_0))\) is then \(\hat{\mu}_1=\sum_{k=1}^K \hat{\mu}_{k,1}\) for large \(K\). This large value of \(K\) introduces a small bias in addition to the \(O(\varepsilon)\) bias that comes from using \(\varepsilon>0\). The \(k\)-th term in \(\hat{\mu}_1\) has variance \(O(n^{-1-1/k})\). An empirical rate like \(n^{-1.1}\) can be interpreted as the integrand being effectively \(10\) dimensional. This doesn’t mean that \(\hat{\mu}_{10,1}\) has the most variance. The variance of the sum includes faster and slower terms. We expect that the variance includes a contribution larger than \(cn^{-1-1/K}\) for some \(c>0\) whenever \(n\geqslant N\) for a possibly very large \(N\). We saw no sign of that slowing in our examples with \(n\leqslant 2^{17}\).

8 Further examples and algorithms↩︎

Here we present some further numerical examples to explore different aspects of the RQMC for WoS problem. We consider examples with a known solution that lets us investigate bias, and we vary the gasket example to make it less favorable to RQMC. We also explore examples not covered by our theory, having nonzero source terms or having \(\Omega\subset\mathbb{R}^3\). What we see in those examples is that our sufficient conditions do not appear to be necessary. We also considered some alternative samplers designed to mitigate the non-smoothness of the WoS problem.

8.1 Unit disk example↩︎

[4] consider the Dirichlet BVP where \(\Omega=B_2(\boldsymbol{0},1)\) with \(\Delta u(\boldsymbol{z})=0\) for \(\boldsymbol{z}\in\Omega\) and \[\begin{align} \label{eq:mhdisk}u(\boldsymbol{z})=\frac{1}{2}\ln[(z_1-2)^2+z_2^2],\textrm{ for } z\in \partial\Omega. \end{align}\tag{9}\] Because \(\Delta u(\boldsymbol{z})=0\) for \(\boldsymbol{z}\ne(2,0)\) the formula in 9 also works for \(\boldsymbol{z}\in\Omega\). We can use that to estimate bias and mean squared error (MSE).

Starting from \(\boldsymbol{z}_0=(0,0.5)\), we run the RQMC-WoS algorithm, terminating once the walk hits the \(\varepsilon\)-shell of \(\Omega\), for \(\varepsilon=10^{-4}\). We run \(100\) replicates for each sample size. The results are shown in Figure 5. The bias in this example is negligible as expected and so the variance is essentially the MSE.

Figure 2 of [5] shows results from a QMC solution on \(n=20\) trajectories attaining an error of roughly \(10^{-8}\). They do not specify the QMC points in use or the starting point \(\boldsymbol{z}_0\), and their Figure 3 shows essentially no difference between MC and QMC solutions for \(n\) ranging to about \(10^5\) in another test problem where they do give \(\boldsymbol{z}_0\) and name the sequence. As a result, we do not think that the reported accuracy in their Figure 2 is directly comparable to the errors shown in our Figure 5.

8.2 Gasket revisited↩︎

We tried another starting point, midway between the two central bore holes, at \(\boldsymbol{z}_0=(0.0030105,0.002839)\). This point is of interest because it is very close to two of the hottest parts of \(\partial\Omega\). What we found was that the RQMC algorithms gave essentially the same accuracy as MC did. Most walks from this starting point hit one of the two nearby bore holes. Since these two bore holes have the same boundary temperature, deviations from that common value mainly come from walks that reach more distant boundary components, which typically require larger values of \(k\). Since RQMC is less favorable for larger \(k\), this can weaken the observed variance rate. To test this explanation, we lowered the temperature of the left bore hole to \(140^\circ\)C while leaving the right bore hole temperature at \(160^\circ\)C. This increased the importance of the walks hitting the nearby bore holes, and then we saw an RQMC variance rate of about \(n^{-1.1}\).

We also investigated the bias in our first gasket example. We averaged estimates from \(100\) MC replicates with \(n=2^{18}\) trajectories each and used that as a reference value. We found very small bias and the variance and MSE were quite close.

8.3 Major circular sector example↩︎

We consider another two-dimensional example from [4]. In polar coordinates, \(\Omega=\{(r,\theta)\mid 0\leqslant r\leqslant 1, -3\pi/2\leqslant\theta\leqslant 0\}\), resembling the 1980s video game character Pac-Man rotated through \(\pi/4\) radians.

This example has a nonzero source term, \(\Delta u(r,\theta)=-(2-r^2)e^{-r^2/2}\) and so we use the WoS update from Section 2.3. That update takes \(3\) uniform variables per step, one for the angle to sample on the circle and two to sample in the disk. Our theory does not cover this case because we have not studied the RQMC error for an integrand that evaluates \(\Delta u\) uniformly over the WoS disk at step \(k\), so it is of interest to see how RQMC performs.

The boundary conditions are \(u(r,0)=e^{-r^2/2}\), \(u(r,-3\pi/2)=-r^{1/3}+e^{-r^2/2}\) and \(u(1,\theta)= \sin(\theta/3)+e^{-1/2}\). The analytic solution is \(u(r,\theta)=r^{1/3}\sin(\theta/3)+e^{-r^2/2}.\)

Using RQMC points in this WoS problem gives the results shown in Figure 6. We again see a slightly better rate of convergence for the MSE of the RQMC-WoS estimator. This example has an anomaly that we have not seen in any other: RQMC with the Niederreiter points performs worse than plain MC sampling does until about \(n=2^{14}\) where it has nearly equal performance. It shows a rate comparable to the other RQMC methods but has a larger constant factor.

8.4 Dumbbell example↩︎

This example is inspired by [35] and [36]. They respectively study the number of critical points of \(u\) and their distance to the boundary of their domain \(\Omega\). Their \(u\) is the solution to the BVP with \[\label{eq:dumbbell} \text{\Delta u(\boldsymbol{z}) = -2,\,\;for \boldsymbol{z}\in \Omega}, \quad\text{and}\quad u(\boldsymbol{z}) = 0, \,\;\text{for \boldsymbol{z}\in \partial\Omega}.\tag{10}\] Their \(\Omega=B_2((-L,0),R)\cup([-L,L]\times [-w,w])\cup B_2((L,0),R)\) is a dumbbell-shaped region, the union of a \(2L\times 2w\) bar with tiny \(w>0\) and two circles of radius \(R<L\). See Figure 7. We will use RQMC-WoS to evaluate \(u\) at a point \(\boldsymbol{z}_0=(L-R,0)\) marked in the figure. Because this problem includes a source term it is not covered by our theorems, and instead provides insight into whether the conditions there are necessary. The interpretation in [36] is that \(u(\boldsymbol{z})\) is the flow velocity of a fluid through a pipe whose cross-sectional shape is \(\Omega\).

We can use equation 3 on this problem. Then an MC or RQMC based estimate of \(u(\boldsymbol{z}_0)\) takes the form \[\hat{u}(\boldsymbol{z}_0)=\frac{1}{2n}\sum_{i=1}^n\sum_{k=1}^{\tau_{i}} \mathrm{dist}(\boldsymbol{z}_{i,k-1},\partial\Omega)^2\] where \(\tau_i = \tau_{i,\varepsilon}=\min\{k\in\mathbb{N}_0\mid \mathrm{dist}(\boldsymbol{z}_{i,k},\partial\Omega)<\varepsilon\}\).

Using MC and RQMC samples, we get the results presented in Figure 8. Again, we observe that using RQMC points brings a better rate of convergence for the variance of the WoS estimator.

8.5 Unit ball example↩︎

Next, we use a three dimensional example from [4], with \(\Omega=B_3(\boldsymbol{0},1)\). The BVP has \(\Delta u(\boldsymbol{z})=0\) for \(z\in \Omega\) with boundary condition \(u(\boldsymbol{z})=[(z_1-2)^2+z_2^2+z_3^2]^{-1/2}\) for \(\boldsymbol{z}\in\partial\Omega\). The analytic solution has the same formula as the boundary condition. To generate a uniform sample on the sphere, we use the hatbox transformation 4 described in Section 2.5

The hatbox transformation \(\psi_0\) is not Lipschitz, so the analysis in Section 6 does not apply to this update. Figure 9 shows an apparent RQMC variance of \(O(n^{-1.14})\) which we take as evidence that our sufficient conditions are not necessary conditions.

8.6 Summary of examples↩︎

Table 1 shows regression fits of log variance (or log MSE when available) versus \(\log(n)\) for our sampling methods and examples. Most of the curves looked to be linear with a little noise, but some of them showed nonlinearities at small \(n\). We used least squares on the data with \(7\leqslant\log_2(n)\leqslant 17\). The Niederreiter points had quite bad performance for the Pac-Man example. While they had the best slope they had a very large constant term.

Table 2 computes variance reduction factors corresponding to the slopes and intercepts in Table 1. The variance reduction factors are modest though it is clear that they improve somewhat for larger \(n\).

8.7 Alternative samplers↩︎

For \(\Omega\subset\mathbb{R}^2\), the position \(\boldsymbol{z}_k\) is a periodic function on \([0,1)^k\), which is an advantage for the lattice sampling methods. The decision to stop or not to stop at step \(k\) introduces a discontinuity which violates the smoothness that lattice samplers need for high accuracy. Additionally, the function \(u\) is not necessarily a smooth function on \(\partial\Omega\).

We ran the algorithm with fixed large values of \(k\) which then give small \(\mathrm{dist}(\boldsymbol{z}_k,\partial\Omega)\) with high probability. This strategy has been used by [37]. After the \(k\) steps, we project \(\boldsymbol{z}_k\) to \(\partial\Omega\) and evaluate \(h\) at that projected point. This removes the discontinuity due to the stopping rule, leaving a fixed-step integrand based on the Lipschitz map \(x\mapsto \boldsymbol{z}_k(x)\). We did not see the lattice sampler outperform other RQMC methods for this algorithm except for very small values, such as \(k=3\), where the bias is too large for the method to be useful. We considered the unit disk case taking \(k=20\) steps toward the boundary, where \(h(\boldsymbol{z})\) is a smooth function. This is the smoothest example we considered and the lattice rule did not show an advantage over the other RQMC methods.

A harmonic function has \(u(\boldsymbol{z}_0)=\mathbb{E}(u(\boldsymbol{z}))\) for \(\boldsymbol{z}\sim\mathbf{U}(S(\boldsymbol{z}_0,r))\). By averaging over radii, the same is true for \(\boldsymbol{z}\sim\mathbf{U}(B(\boldsymbol{z}_0,r))\). [37] include an algorithm where the first step is uniform over \(B(\boldsymbol{z}_0,r_0)\) to make the estimate a smoother function of \(\boldsymbol{z}_0\). We investigated such steps but they did not produce an algorithm that made lattice methods dominant.

9 Discussion↩︎

We have replaced MC sampling by RQMC sampling in the WoS algorithm for a collection of BVP problems in \(\mathbb{R}^2\) and \(\mathbb{R}^3\). We see rates that are close to \(n^{-1-1/\tilde{k}}\) for some ‘effective’ dimension \(\tilde{k}\) over values of \(n\) up to \(2^{17}\). The actual error has contributions from many values of \(k\), and there is likely to be a nonzero contribution from \(k=K\), leading to a rate of only \(O(n^{-1-1/K})\), but there is no evidence that such a pessimistic rate will be relevant for sample sizes near the ones we have used. The unit disk example has a very smooth boundary function and illustrates the theorem from [7].

The Pac-Man example uses \(3k\) uniform variables to take \(k\) steps because it also samples disks at each step. The inferior convergence rate we saw for it could be due to this increased dimension. The integrand is only discontinuous in \(k\) of those variables used to define whether the walk terminates in \(k\) steps.

Materials and methods↩︎

We used ChatGPT, Claude and Gemini in parts of this work. The versions we used were from the early months of 2026. Those AIs helped us find relevant literature much more quickly than we could have done with a search engine. Many suggested references were not actually useful and had to be discarded but others were critical. The AIs found weaknesses in some early versions of our proofs and made suggestions. Quite often the AI suggestions were naive about subtleties such as the varying definitions of rectifiability and we rejected many AI suggestions. On the other hand, it was an AI that noticed that the second claim in Theorem 4 follows from the same argument as the first claim. The AIs made some suggestions about writing style and we adopted a few of those for clarity but we wrote everything ourselves. AI found some typos and incomplete sentences that we corrected. They also generated, following our instructions, the Python code used to draw the gasket and dumbbell domains shown in Figure 1 and 7 respectively, and the R code used to produce Figure 4, which illustrates the chambers. An AI suggested the term ‘chamber’. The AIs also pointed us toward useful theorems on manifolds, after which we identified versions that fit the paper without requiring substantial additional background. When we were investigating an apparent sign error in our computations for the dumbbell example, an AI recognized that [12] uses the negative semidefinite Laplacian. The AIs also assisted with LaTeX syntax in resolving formatting issues and with improving several BibTeX entries.

Acknowledgments↩︎

We thank Rohan Sawhney and Yang Liu for some helpful comments. The second author visited Bob Carpenter’s group at the Flatiron Institute and had productive discussions there with Misha Padidar, Michael Czekanski, Dan Fortunato and Charles Epstein.

References↩︎