Category: Research

Project Mentor: Bailey Miller

1 Introduction

A classic problem in computer graphics is the interpolation of boundary values into a region’s interior. More formally, given a domain \( \Omega \) and values prescribed on its boundary \( \partial \Omega \), we seek a smooth function \( u: \mathbb{R}^n \to \mathbb{R} \) that satisfies the boundary value problem

$$\begin{align}
\Delta u(x) = 0 \text{ in } \Omega, \
\quad u(x) = g(x) \text{ on } \partial \Omega,
\end{align}$$

where \(g(x)\) is a given boundary function. This is the well-known Laplace’s equation with Dirichlet boundary conditions. A standard approach for solving such PDEs is the finite element method (FEM), which discretizes the domain into a mesh and computes a weak solution. However, the accuracy of FEM depends critically on the quality of the mesh and how well it approximates the domain geometry. In contrast, Monte Carlo methods offer an appealing alternative — particularly the walk on spheres (WoS) algorithm [Muller, 1956], which avoids meshing altogether. WoS exploits the Mean Value Property of harmonic functions (i.e., functions satisfying \( \Delta u = 0) \), which states:

$$\begin{align}
u(x) = \frac{1}{|\partial B(x, r)|} \int_{\partial B(x, r)} u(y) \ \mathrm{d}y,
\end{align} \quad \text{(*)}$$

for any ball \( B(x, r) \subseteq \Omega \). This property implies that \( u(x) \) can be estimated by recursively sampling points on the boundary of spheres centered at \( x \), continuing until a point lands on \( \partial \Omega \), where the value is known from \( g(x) \). Compared to FEM, WoS scales much better with geometric complexity and does not require certain properties of the mesh like watertightness, free of self-intersections, manifoldness, etc. While WoS does not generalize to all classes of PDEs, it offers many advantages even beyond geometric robustness (See Sawhney and Crane 2020, Section 1).

While WoS is an unbiased estimator and is largely agnostic to geometric quality, it can suffer from high variance, often requiring many samples to converge and thus incurring high computational cost. A common way to reduce variance is to use a control variate, and prior work has explored applying them to WoS. The main challenge in using control variates is finding a suitable reference function (one with a known integral) that also closely matches the target function. One research direction to address this challenge is neural control variates—constructing the control variate function using neural networks (See Li et al. 2024, who applied this approach to WoS). While their method leverages a clever technique known as automatic integration [Li et al. 2024] to perform closed-form integration, it still suffers from substantial overhead during both training and inference due to the use of MLPs. In contrast, we propose harmonic control variates, where harmonicity is built into the control variate representation by design, enabling closed-form integration directly as a function evaluation. Moreover, our representation is fitted using a simple least squares procedure, which is significantly more efficient than training a neural network via gradient descent. It also allows for much faster evaluation at inference time.

Specifically, our approach consists of two main steps. (1) We construct a harmonic representation as a weighted sum of harmonic kernels (e.g., Green’s functions), and fit the weights/coefficients using a sparse set of points in the domain along with their WoS-based estimates. This fitting reduces to solving a linear system of the form \( \Phi \sigma = u \) for the weight vector \( \sigma \). The resulting function is harmonic by construction, and although it provides a useful denoised estimate, it remains biased with respect to the true solution \( u(x) \) almost everywhere in the domain. (2) To correct for this bias, we use the harmonic representation as a control variate, thus termed harmonic control variate. Because the representation is harmonic, its mean over any sphere can be computed in closed form via simple function evaluation. We then define a residual Laplace problem, where the boundary condition is the difference between the true boundary data and our representation. Solving this residual problem with WoS and adding the result back to the harmonic representation yields an unbiased estimator for \( u(x) \). This control variate estimator reduces variance whenever the representation closely approximates \( u(x) \) on the boundary.

Our harmonic control variate achieves significant variance reduction compared to naive WoS across a range of geometries, yielding an average 40% reduction in RMSE compared to the naive baseline at 128 walks. In addition to its effectiveness, our method is simple to implement: the entire pipeline can be reproduced in Python with minimal code, leveraging just a few function calls to the underlying integrated C++ library Zombie [Sawhney and Miller 2023]. To illustrate its flexibility, we also extend our approach to settings with Neumann boundary conditions and with an alternative Monte Carlo scheme Walk on Stars [Sawhney, Miller, et al. 2023]. We believe that this classically inspired yet practical control variate representation (a sum of harmonic kernels) offers useful insights for the broader computational science and graphics communities.

2 Background

2.1 Walk on Spheres

The Walk on Spheres (WoS) algorithm [Muller 1956] can be motivated in two complementary ways. One is via Kakutani’s principle, which states that the solution value \( u(x) \) of a harmonic function at any point \( x \in \Omega \) is equal to the expected boundary value \( u(x) = \mathbb{E}[g(y)] \), where \( y \in \partial \Omega \) is the first boundary point reached by a Brownian motion starting at \( x \). We focus on a second, more constructive motivation coming from the mean value property of harmonic functions (see eq. (*)). This property immediately suggests a single-sample Monte Carlo estimator:

$$\begin{align}
\langle u(x) \rangle = u(y),
\quad y \sim \mathscr{U}(\partial B(x,r)),
\end{align}$$

where \( \mathscr{U}(\partial B(x,r)) \) denotes the uniform distribution on the sphere of radius \( r \) centered at \( x \). Specifically, the mean value property suggests that this estimator is unbiased:

$$ \begin{align}
\mathbb{E}[\langle u(x) \rangle]
\,=\, \mathbb{E}_{y \sim \mathscr{U}(\partial B(x,r))}[u(y)] \,=\, \int_{\partial B(x,r)} p(y) u(y) \,\mathrm{d}y
\,=\, \frac{1}{|\partial B(x,r)|}\int_{\partial B(x,r)} u(y) \,\mathrm{d}y
\,=\, u(x),
\end{align} $$

where \( p(y) = 1/|\partial B(x,r)| \) denotes the uniform probability density on \( \partial B(x, r) \).

However, the difficulty is that \( u(y) \) is unknown unless \( y \) lies on the boundary \( \partial \Omega \). To address this, the mean value property is applied recursively: each sampled point becomes the new center for the next sphere, and the process continues until the boundary is reached. Since naïve recursion would never terminate, the algorithm introduces a tolerance \( \epsilon \). Specifically, once a point \( x_k \) lies within distance \( \epsilon \) of the boundary, we project it to the nearest boundary point \( \bar{x}_k \in \partial \Omega \) and evaluate the boundary condition \( g(\bar{x}_k) \). This serves as the base case of the recursion.

Formally, letting \( \Omega_\epsilon \) denote the \( \epsilon\)-shell adjacent to the boundary (illustrated in Figure 3), the recursive rule is

$$\begin{align}
\langle u(x_k) \rangle =
\begin{cases}
g(\bar{x}_k), & x_k \in \Omega\epsilon, \\
u(x_{k+1}), & \text{otherwise},
\end{cases}
\end{align}$$

with the next sample \( x_{k+1} \sim \mathscr{U}(\partial B(x_k, r_k)) \), where \( r_k \) is the largest radius such that the ball remains fully contained in \( \Omega \). As the recursion telescopes, the final estimator is simply

$$\begin{align}
\langle u(x_0) \rangle = g (\bar{x}_{n_i}^{(i)}), \quad x_n \in \Omega_{\epsilon},
\end{align}$$

where \( \bar{x}_n \) is the projected boundary point at termination.

By performing \( N \) independent walks and averaging their outcomes, we obtain a multi-sample Monte Carlo estimator of the solution:

$$\begin{align}
\langle u(x_0) \rangle_N = \frac{1}{N} \sum_{i=1}^N g\left( \bar{x}_{n_i}^{(i)}\right),
\end{align}$$

where each \( \bar{x}_{n_i}^{(i)} \in \partial \Omega \) is the endpoint of the \( i \)-th walk, reached on its \( n_i \)-th step. This estimator remains unbiased as \( \epsilon \to 0 \), but its efficiency is governed by the variance of the boundary condition \( g(x) \). In fact, the variance of the estimator is directly proportional to that of \( g(x) \):

$$\begin{align}
\mathrm{Var}[\langle u(x_0) \rangle_N]
= \frac{1}{N}\,\mathrm{Var}\left[g\left(\bar{x}_{n_i}^{(i)}\right)\right].
\end{align}$$

Intuitively, if \( g(x) \) varies significantly across the boundary region where walks terminate, the estimates \( \langle u(x_0) \rangle \) will fluctuate sharply between different exit points, leading to a noisy approximation of the solution. Conversely, if \( g(x) \) is nearly constant around the relevant boundary region, most walks yield similar values, and the estimator stabilizes. This observation is central to our two-pass variance-reduction strategy (see Section 3.2).

2.2 Control Variates

To build intuition, let us first consider the simple task of integrating a scalar function in one dimension. Suppose we want to estimate

$$\begin{align}
F = \int_{a}^{b} f(x)\, \mathrm{d}x,
\end{align}$$

using Monte Carlo (MC) sampling. A naive single-sample Monte Carlo estimator is

$$\begin{align}
\langle F \rangle = |b-a|\, f(x), \quad x \sim \mathscr{U}[a,b].
\end{align}$$

The idea of a control variate is to introduce a function (c(x)) that (i) closely approximates (f(x)) and (ii) has a known closed-form integral, denoted (C). We can rewrite the target integral as

$$\begin{align}
F = \int_a^b c(x)\, \mathrm{d}x + \int_a^b f(x) – c(x) \, \mathrm{d}x
= C + \int_a^b f(x) – c(x) \, \mathrm{d}x.
\end{align}$$

This yields a single sample Monte-Carlo estimator with control variate:

$$\begin{align}
\langle F \rangle^{(\text{cv})} = |b-a| \big(f(x) – c(x)\big) + C, \quad x \sim \mathscr{U}[a,b].
\end{align}$$

It is clear that it remains unbiased. Moreover, it achieves lower variance whenever the residual \( r(x) := f(x) – c(x) \) is less variable than the original function \( f(x) \):

$$\begin{align}
\mathrm{Var} [\langle F \rangle^{(\text{cv})} ]
<
\mathrm{Var}[\langle F \rangle]
\quad \iff \quad
\mathrm{Var}[f(x) – c(x)] < \mathrm{Var}[f(x)].
\end{align}$$

Control Variates for Walk-on-Spheres. Motivated by this principle, prior work [Li et al. 2024] incorporated control variates directly into the mean value formulation of WoS. At each step of the walk, they replace the boundary condition with a residual correction, and by eq. (*),

$$\begin{align} u(x) = \frac{1}{\lvert \partial B(x,r)\rvert}\int_{\partial B(x,r)} c(y)\,\mathrm{d}y + \frac{1}{\lvert \partial B(x,r)\rvert}\int_{\partial B(x,r)} \big(u(y)-c(y)\big)\,\mathrm{d}y \end{align}$$

where \( C(x) = \frac{1}{\lvert \partial B(x,r)\rvert}\int_{\partial B(x,r)} c(y)\,\mathrm{d}y \) and \( c : \mathbb{R}^d \to \mathbb{R} \) is the control variate. This introduces the challenge of designing a control variate \(c\) that is both expressive and admits a tractable integral \(C(x)\). Li et al. [2024] addressed this by parameterizing \( c \) with a neural network and using automatic integration to evaluate \(C(x)\). However, this approach requires costly neural network evaluations at every step of the walk and additional training via gradient descent. In contrast, we take the insight that if \( c \) is also harmonic, then \( C(x) = c(x) \) by eq. (*). Our method only evaluates the control variate twice per walk and fits it directly via least squares, making it significantly more efficient while still achieving variance reduction.

3 Method

Our method introduces a harmonic representation that serves as a control variate for variance reduction. We first describe the formulation of this harmonic representation and how it can be efficiently solved in Section 3.1. Next, we explain how this representation functions as a control variate to reduce variance in Monte Carlo estimation in Section 3.2. Finally, we present implementation details in Section 3.3. An overview of our method is illustrated in Figure 5.

3.1 Mixture of Harmonic Kernels Representation

We propose a simple, kernel-based representation for approximating a target harmonic function \( u(x) \), subject to Dirichlet boundary conditions \( u(x) = g(x) \) on \( \partial \Omega \). Specifically, we approximate \( u(x) \) using a mixture of harmonic kernels \( \phi_j(\,\cdot\,) \):

$$\begin{align}
c(x) = \sum_{j=1}^{M} \sigma_j \,\phi_j(x),
\end{align}$$

where each \( \phi_j(x) \) is itself harmonic. By linearity, the mixture \( c(x) \) is also harmonic. To ensure that \( c(x) \) remains free of singularities inside the domain \( \Omega \), the kernels are chosen carefully. The coefficients \( \{\sigma_j\}_{j=1}^M \) are free parameters that must be fitted to approximate \( u(x) \).

There are many possible choices of harmonic kernels. A canonical example in two dimensions is the Green’s function,

$$\begin{align}
\phi_j(x) = G^{\mathbb{R}^2}(x, y_j) = \frac{1}{2\pi} \log\bigl(\|x – y_j\|_2\bigr),
\end{align}$$

where the centers \( \{y_j\} \) are placed outside of \( \Omega \) to avoid singularities within the domain. We will introduce another kernel choice in Section 4.

Fitting the coefficients: To determine the coefficients \( \sigma_j \), we first sample a sparse set of interior points \( \{x_i\}_{i=1}^N \subset \Omega \). At each point, we estimate \( \langle u(x_i) \rangle \) using Walk-on-Spheres (WoS), which leverages the known boundary condition. This produces training tuples \( \{(x_i, \langle u(x_i) \rangle) \}_{i=1}^N \), where \( c(x_i) \) should approximate \( \langle u(x_i) \rangle \). Each tuple gives a linear constraint on \( \sigma_j \), leading to the linear system:

Solving this system yields the optimal coefficients {\( \sigma_j \)}. This can also be written as an optimization problem (the least square problem):

$$\begin{align}
\underset{\sigma \in \mathbb{R}^{M}}{\mathrm{argmin}} \,
\|\Phi \sigma – \langle u \rangle \|_{2}^{2},
\label{eq:optimize_lin_sys} \quad \text{(***)}
\end{align}$$

such that we can introduce some regularization terms (see Section 3.3). An example of this fitting procedure using the Green’s kernel is illustrated in Figure 6.

Intuitively, this representation projects the set of sparse and noisy estimates of \( u(x_i) \), i.e., the vector \( \langle u \rangle \), onto the subspace spanned by the harmonic kernels \( \{\phi_j(x)\} \). Since this fit is based only on a limited set of points and kernels and does not enforce correctness elsewhere, this representation is a biased estimator. Therefore, while this representation \( c(x) \) serves as a useful denoiser due to the smoothness of the kernels, it is not a sufficiently accurate approximation of the true solution \( u(x) \), see Figure 5. In the next section, we show how it can be leveraged as a control variate to reduce variance and obtain an unbiased estimate of \( u(x) \) throughout \( \Omega \).

3.2 Two-Pass with Control Variates

The main advantage of our representation \( c(x) \) being harmonic is that it satisfies \( \Delta c(x) = 0 \) in \( \Omega \). This allows us to define a residual PDE:

$$\begin{align}
\Delta\left[ u(x) – c(x) \right] = 0 \text{ in } \Omega, \
\quad u(x) – c(x) = g(x) – c(x) \text{ on } \partial \Omega,
\end{align}$$

which corresponds to solving for the harmonic residual function \( r(x) := u(x) – c(x)\). Since the Dirichlet boundary condition for \( r(x) \) is known, namely \( r(x) = g(x) – c(x) \) on \( \partial \Omega \), we can apply the Walk-on-Spheres (WoS) algorithm to estimate \( \langle r(x) \rangle \). Our control variate estimator for \( u(x) \) is then given by

$$\begin{align}
\langle u(x) \rangle^{(\text{cv})}
:= c(x) + \langle r(x) \rangle.
\end{align}$$

The argument for why this estimator is still unbiased and has smaller variance is similar to the 1D case outlined in Section 2.2. This estimator is unbiased since

$$\begin{align}
\mathbb{E}[\langle u(x) \rangle^{(\text{cv})}]
= \mathbb{E}[c(x) + \langle r(x) \rangle]
= c(x) + \mathbb{E}[\langle r(x) \rangle]
= c(x) + r(x)
= u(x),
\end{align}$$

where we used the fact that \( \mathbb{E}[\langle r(x) \rangle] = r(x) \) due to the unbiasedness of WoS. Therefore, as the number of random walks increases, the control variate estimator \( \langle u(x) \rangle^{(\text{cv})} \) converges to the true solution \( u(x) \). For variance reduction, observe that

$$\begin{align}
\mathrm{Var}[\langle u(x) \rangle^{(\text{cv})}]
= \mathrm{Var}[c(x) + \langle r(x) \rangle]
= \mathrm{Var}[\langle r(x) \rangle],
\end{align}$$

since \( c(x) \) is deterministic. Thus, we have:

$$\begin{align}
\mathrm{Var}[\langle u(x) \rangle^{(\text{cv})}] < \mathrm{Var}[\langle u(x) \rangle]
\quad \text{iff} \quad
\mathrm{Var}[\langle r(x) \rangle] < \mathrm{Var}[\langle u(x) \rangle].
\end{align}$$

As discussed earlier in the end of Section 2.1, variance is reduced when \( c(x) \) provides a good approximation to \( u(x) \) on the boundary. In this case, the residual boundary condition \( g(x) – c(x) \) has a smaller range, resulting in reduced variance in the WoS estimates for \( r(x) \) and, by extension, for \( u(x) \).

3.3 Implementation Details

We next describe several implementation choices in setting up and solving the linear system, as well as in implementing control variates.

Linear System Setup. Our formulation yields an \( N \times M \) linear system, where \( N \) is the number of constraints (sampled interior points) and \( M \) is the number of kernels. We deliberately set \( N > M \), making the system over-constrained. This avoids ill-posedness and eliminates the possibility of infinitely many solutions. To further enhance stability, we adopt Kd-tree based adaptive sampling, which distributes interior points evenly across the domain rather than allowing clusters in small regions. This reduces the risk of nearly colinear rows in the system matrix, which would otherwise harm numerical conditioning. Finally, we add regularization to eq. (***) to improve robustness near the boundary. Without it, the fitted representation tends to produce large coefficient magnitudes, leading to spikes and non-smooth artifacts. We experimented with both \( \ell_2 \) (ridge) regularization, which discourages large weights, and \( \ell_1 \) (lasso) regularization, which can additionally drive coefficients to zero. Both were effective, but we ultimately used \( \ell_2 \) in our experiments.

Control Variates with Zombie [Sawhney and Miller 2023]. We build on the Zombie library, which provides efficient C++ implementations of Monte Carlo geometry processing methods along with Python bindings. In the baseline setting, Zombie takes as input the geometry together with boundary conditions sampled on a grid. During Walk-on-Spheres (WoS), Zombie interpolates these grid values whenever a trajectory hits the boundary, and then outputs estimated interior values on a user-specified grid resolution. Our variance reduction technique using control variates requires only a minimal change to this setup. Instead of supplying the original boundary condition, we provide the boundary condition of the residual function. Zombie then performs WoS to estimate the residual inside the domain. The final solution is obtained by simply adding back the evaluation of \( c(x) \) on the grid to these residual estimates. This simple modification makes the method straightforward to implement while still delivering effective variance reduction.

4 Experiments

We explore two types of harmonic kernels for our harmonic representation. The first is the 2D fundamental solution (Green’s function): for \( x,y\in\mathbb{R}^2 \),

$$
\begin{align}
G^{\mathbb{R}^2}(x,y)=\frac{1}{2\pi}\log\|x-y\|_2,
\end{align}
$$

which satisfies \( \Delta_x G^{\mathbb{R}^2}(x,y)=\delta(x-y)\) and thus harmonic in \(x\) away from the singularity at \(y\); it underlies methods such as the Mixture of Fundamental Solutions and boundary-element approaches.

The second is the angle kernel obtained by integrating an angular density along a directed curve. For 2D polylines, the kernel admits a closed form: writing the (i)-th directed segment as \( \ell_i=[a_i,b_i] \), the associated angle kernel is defined as

$$
\begin{align}
\theta_i(x) = \operatorname{atan2}
\left(
\frac{\det(b_i-x,a_i-x)}{(b_i-x)\!\cdot\!(a_i-x)}
\right),
\end{align}
$$

with \( \det(u,v)=u_1v_2-u_2v_1 \). The \( \operatorname{atan2} \) formulation is numerically robust and returns the signed angle in \((-\pi,\pi]\). One can verify that the kernel is harmonic for \(x\notin\ell_i\); across the segment the angle exhibits \(\pm 1\) jump discontinuities. This angle representation is inspired by generalized winding-number and boundary-integral formulations (see Jacobson, Kavan, and Sorkine [2013])

To verify our 2-pass method, we compared it against a reference solution obtained by running WoS with 32,768 walks, which we take as ground truth since this number is sufficient for convergence. For evaluation, we tested using between 4 and 1024 walks, and compared three methods—naive WoS, control variates with the Green’s function, and control variates with the angle kernel. We measured the median RMSE across 20 unique domains relative to the reference solution.

For all experiments, we used the following settings. The grid resolution was \( 256 \times 256 \) and the epsilon shell for WoS was \( 0.001 \). The number of basis functions for the denoiser was \( 80 \) and the number of walks for the samples for the denoisers was \( 128 \). We used \( \ell_2 \) regularization with \( \alpha = 0.1 \). The Dirichlet boundary condition we used is given by:

$$
g(x, y) = \sin\big(5\pi y\big) + \cos\big(2\pi x\big)
$$

Based on RMSE, we observed that our control variate method for both kernels outperformed the naive WoS on average. The control variate with Green’s theorem showed significant linear improvement as the number of walks increased. The angle kernel did not yield as good results as Green’s; however, in the vast majority of experiments, the angle kernel control variate still achieved better RMSE than naive WoS. The only exception we observed was on more complex domains with a high number of walks. In contrast, for the Green’s kernel, we did not encounter any domain or walk count in which naive WoS outperformed the control variate method.

For computation time, we observed that the denoising process has a fixed amount of overhead. Therefore, as the number of walks increases, the denoising time becomes a negligible factor in the overall computation time. The bottleneck is then essentially the same as the naive method: the random walks.

5 Extensions

5. 1 N-Pass Gradient Descent Approach

Building upon the 2-pass harmonic control variate method, we experimented with an n-pass gradient descent approach to iteratively refine the control variate using Monte Carlo estimates as supervision. This method was motivated by the Neural Walk-on-Spheres framework [Nam, Berner, and Anandkumar 2024], which uses a general-purpose neural network \( u_\theta(x) \) as the learned control variate. Our method instead optimizes the coefficient weights \( \sigma ∈ ℝ^{M} \) of a harmonic control variate \( c(x) \), adhering to the structure of the underlying PDE.

The weights \( \sigma \) are initialized from a small set of sample points via the linear system (**), and stochastic gradient descent aims to then minimize the discrepancy between the control variate \( c(x) \) and Monte Carlo estimates \( \langle u(x) \rangle \).

Our optimization loop is implemented as follows:

1. Sample a batch of interior points \( x_i \) and compute Monte Carlo estimates \( \langle u(x_i) \rangle \) for each batch point using WoS
2. Evaluate the current control variate at each \( x_i \), and compute the mean square error (MSE) loss $$\mathscr{L} = \frac{1}{B} \sum_{i=1}^{B} \left( c(x_i) – \langle u(x_i) \rangle \right)^2$$
3. Perform a gradient step on \( \mathscr{L} \) with respect to \( \sigma \) via Adam optimizer

For the same total number of walks, Figure 11 illustrates a significant improvement in RMSE compared to naive WoS, with the contours indicating a much less noisy interpolation of the boundary for our n-pass solution.

Further improvements include accelerating the repeated computation of \( c(x) \), adaptively adjusting the learning rate, and experimenting with different weight initializations. For example, when artifacts appeared in the denoiser due to a large number of monopoles, these artifacts remained present even after the optimization, since weights were initialized based on the denoiser.

5.2 Neumann Boundary Conditions

Many real-world PDE problems involve boundary behavior that is best described by flux, rather than a set of fixed values. Neumann conditions prescribe the normal derivative of the solution along the boundary. In the case of mixed boundary conditions, we seek a function \( u \) satisfying:

$$
\Delta u = 0 \text{ in } \Omega, \quad u = g \text{ on } \partial\Omega_D, \quad \frac{\partial u}{\partial n} = h \text{ on } \partial\Omega_N
$$

where the boundary \( \partial \Omega = \partial \Omega_D \cup \partial\Omega_N \) is partitioned into Dirichlet and Neumann components. A 2-pass approach in this setting would first construct a control variate \( c(x) \approx u(x) \), and estimate the residual \( u(x) – c(x) \) which satisfies:

$$\Delta(u – c) = 0 \text{ in } \Omega, \quad (u – c) = (g – c) \text{ on } \partial\Omega_D, \quad \frac{\partial}{\partial n}(u – c) = h – \frac{\partial c}{\partial n} \text{ on } \partial\Omega_N$$

In particular, when we choose Green’s function for our control variate, the Neumann boundary condition is corrected via the Poisson kernel,

$$\frac{\partial}{\partial n} (u – c)(x) = h(x) – \sum_i \frac{\partial}{\partial n} P(x, y_i) \, \sigma_i,$$

as the Poisson kernel is the normal derivative of the Green’s function with respect to the boundary point:

$$P(x, y) = -\frac{\partial G}{\partial n_y}(x, y), \quad y \in \partial\Omega$$

By using this residual-based formulation, we apply the control variate to both value and flux constraints, enabling variance reduction even in domains with mixed boundaries.

To handle Neumann boundary conditions, the original WoS algorithm needs to be modified to address the behavior of walks near Neumann boundaries, the radius of closed balls around each point, and conditions for walk termination. The Walk on Stars algorithm (Sawhney et al., 2023) does exactly this. Walk on Stars generalizes WoS by replacing ball-shaped regions with stars — a star being any shape where every boundary point can be “seen” from the centre. At Neumann boundaries, Walk on Stars simulates reflecting Brownian motion — when a walk reaches a Neumann boundary, it is reflected according to the boundary’s normal vector and picks up a contribution of the Neumann condition, \( h \).

Figure 12 illustrates our results with a zero-valued Neumann boundary condition.The 2-pass method still reduces the RMSE compared to Naive WoS with the same number of walks. However, the contours indicate somewhat increased noise near the reflecting boundary compared to the absorbing boundary on our solution.

5.3 Grouping Techniques to Reduce Overfitting and Improve Stability

As detailed in Section 3.3, a major consideration when estimating control variates is the conditioning and placement monopoles. To further address this, we experimented with grouping techniques that compress nearby or similar samples and assign them a shared influence, namely:

• Greedy grouping: Choose a point and assign all other points within a fixed radius to one group. Each groups shares either a weight or a single basis function, effectively reducing the number of active degrees of freedom in the system. Our experiments show a reduction in RMSE via grouping for the angle kernel (see Figure 13), but contours indicate a slight increase in noise in some regions.
• DBSCAN: This density-based clustering technique automatically identifies dense clusters of points. Each cluster is treated as a group, and only a representative or averaged contribution is retained, while the rest are collapsed into it.

6 Conclusions and Future Directions

We introduced harmonic control variates: a simple, closed-form control variate for Walk-on-Spheres built from mixtures of harmonic kernels. By fitting a lightweight least-squares model once and efficiently evaluating it in closed form during sampling, we turn WoS into a two-pass estimator that remains unbiased while substantially lowering variance. Empirically, our method delivers consistent accuracy gains across diverse geometries—e.g., an average 40% RMSE reduction at 128 walks—while adding only minimal overhead compared to previous neural control variate approaches; as walks dominate runtime, the extra cost is effectively negligible. The result is a method that is easy to implement (a small change to a standard WoS pipeline), robust to geometric complexity, and markedly more sample-efficient.

Future work. We plan to release a plug-in that drops into existing WoS codebases and to continue improving our method for even stronger variance reduction at essentially the same runtime.

References

Jacobson, Alec, Ladislav Kavan, and Olga Sorkine [2013]. “Robust Inside-Outside Segmentation using Generalized Winding Numbers”. In: ACM Trans. Graph. 32.4.

Li, Zilu et al. [2024]. “Neural Control Variates with Automatic Integration”. In: ACM SIGGRAPH 2024 Conference Papers. SIGGRAPH ’24. Denver, CO, USA: Association for Computing Machinery. isbn: 9798400705250. doi: 10.1145/3641519. 3657395. url: https://doi.org/10.1145/3641519.3657395.

Madan, Abhishek et al. [July 2025]. “Stochastic Barnes-Hut Approximation for Fast Summation on the GPU”. In: Proceedings of the Special Interest Group on Computer Graphics and Interactive Techniques Conference Conference Papers. SIGGRAPH Conference Papers ’25. ACM, pp. 1–10. doi: 10.1145/3721238.3730725. url: http://dx.doi.org/10.1145/3721238.3730725.

Muller, Mervin E. [1956]. “Some Continuous Monte Carlo Methods for the Dirichlet Problem”. In: The Annals of Mathematical Statistics 27.3, pp. 569–589. doi: 10.1214/aoms/1177728169. Url: https://doi.org/10.1214/aoms/1177728169.

Nam, Hong Chul, Julius Berner, and Anima Anandkumar [2024]. “Solving Poisson Equations Using Neural Walk-on-Spheres”. In: Forty-first International Conference on Machine Learning

Sawhney, Rohan and Keenan Crane [Aug. 2020]. “Monte Carlo geometry processing: a grid-free approach to PDE-based methods on volumetric domains”. In: ACM Trans. Graph. 39.4. issn: 0730-0301. Doi: 10.1145/3386569.3392374. url: https://doi.org/10.1145/3386569.3392374.

Sawhney, Rohan and Bailey Miller [2023]. Zombie: Grid-Free Monte Carlo Solvers for Partial Differential Equations. Version 1.0.

Sawhney, Rohan, Bailey Miller, et al. [July 2023]. “Walk on Stars: A Grid-Free Monte Carlo Method for PDEs with Neumann Boundary Conditions”. In: ACM Trans. Graph. 42.4. issn: 0730-0301. doi: 10.1145/3592398. url: https://doi.org/10.1145/3592398.

Our work was carried out under the mentorship of Yi Zhuo from Adobe Research, whose guidance was instrumental and very much appreciated. We also thank Sanghyun Son , the first author of the original DMesh and DMesh++ , for generously sharing his insights and providing invaluable input.

Article Composed By :
Gabriel Isaac Alonso Serrato
in collaboration with
Amruth Srivathsan

Click Below to Try :

INTRODUCTION

This article introduces Differentiable Mesh++ (DMesh++), exploring what a mesh is, why discrete mesh connectivity blocks gradient flow, how probabilistic faces restore differentiability, and how a local, GPU-friendly tessellation rule replaces a global weighted Delaunay step. We close with a practical optimization workflow and representative applications.

WHAT IS A MESH?

A triangle mesh combines geometry and topology. Geometry is given by vertex coordinates (named points in space such as A, B, C, and D) while topology specifies which vertices are connected by edges and which triplets form triangular faces (e.g., ABC, ACD). Together, these define shape and structure in a compact representation suitable for rendering and computation on complex surfaces.

Why Traditional Meshes Fight Gradients?

Classical connectivity is discrete: a triangle either exists or it does not. Gradient-based learning, however, requires smooth changes. When connectivity flips on or off, derivatives vanish or become undefined through those choices, preventing gradients from flowing across topology changes and limiting end-to-end optimization

Differentiable Meshes

Gradients are essential for learning-based pipelines, where losses may be defined in image space, perceptual space, or 3D metric space.

Differentiable formulations make it possible to propagate error signals from rendered observations or physical constraints back to vertex positions, enabling tasks such as inverse rendering, shape optimization, and end-to-end 3D reconstruction. This bridges the conceptual gap between discrete geometry and continuous optimization, allowing mesh-based representations to function as both modeling primitives and trainable structures.

How DMesh++ Works

DMesh++ makes meshing a differentiable, GPU-native step. Each point carries learnable features (position, a “realness” score ψ, and optional channels).

A soft tessellation function turns any triplet into a triangle weight rather than a hard on/off face.

During training, losses back-propagate to both positions and features, so connectivity can evolve smoothly.

The “One-Line” Idea in DMesh++

In DMesh++, we replace WDT ( Λ_wdt )

with a local Minimum-Ball test ( Λ_min ) when scoring each candidate face. Instead of rebuilding a global structure, we only check the smallest ball touching the face’s vertices and see if any other point lies inside it. This can be done quickly with k-nearest neighbor search on the GPU. The face probability is now:

Λ(F) = Λ_min(F) × Λ_real(F)

• Λ_min(F) says, “Is this face even eligible as part of a good tessellation?” (now via the Minimum-Ball test)

• Λ_real(F) says, “Is this face on the real surface?” (a differentiable min over the per-point realness ψ)

This swap reduces the effective evaluation cost from O(N) (WDT) to O(log N) with k-NN and parallelization, in practice, achieving speeds up to 16× faster in 2D and 32× faster in 3D, with significant memory savings.

Minimum-Ball, Explained

Pick a candidate face F (segment in 2D, triangle in 3D). Compute its minimum bounding ball B_F ​: the smallest ball that touches all vertices of F.

Now check: is any other point strictly inside that ball? If no, F passes. If yes, reject. That’s the Minimum-Ball condition. (It’s a subset of Delaunay, so you inherit good triangle quality and avoid self-intersections.)

DMesh++ makes this soft and differentiable. It computes a signed distance between B_F​ and the nearest other Point and runs it through a sigmoid to get Λ_min :

Minimum-Ball Playground:

Click here to try:

How to use it?

Click on the canvas to drop points. Every pair of points becomes a candidate edge F. We draw the minimum ball B_F​ (the smallest circle touching the two endpoints) and compute a soft eligibility via a sigmoid of its clearance to the nearest third point:

Because faces are probabilistic, the topology can change cleanly during optimization. You don’t hand-edit connectivity; it emerges from the features.

Recap: A Practical Optimization Workflow

Start: Begin with points (positions + realness ψ)

Candidates: Form nearby faces (via k-NN)

Minimum-Ball: For each face, find the smallest ball touching its vertices (if no point is inside, it’s valid) → gives Λₘᵢₙ

Weight by Realness: Multiply by vertex ψ to get final probability Λ = Λₘᵢₙ × Λᵣₑₐₗ

Optimize: Compute the loss (e.g., Chamfer), and backpropagate to update the points and ψ

Output: Keep high-Λ faces → final differentiable mesh

Some Applications of DMesh++

DMesh++ demonstrates how differentiable meshes can power real applications. Here are some ideas:

Generative Design: Train models that evolve shapes and topology for creative or engineering tasks without manual remeshing.

Adaptive Level of Detail: Automatically refine or simplify geometry based on learning signals or rendering needs.

Feature-Aware Reconstruction: Learn geometric structure, merge fragments, or fill missing regions directly through optimization.

Neural Simulation: Enable physics or visual simulations where mesh shape and motion are optimized end-to-end with gradients.

Future Directions

Looking ahead, future work on DMesh++ could focus on improving usability and scalability, for example, by developing more accessible APIs, real-time visualization tools, and integration with common 3D software. On the algorithmic side, enhancing the tessellation functions for higher precision and stability, expanding support for non-manifold and hybrid volumetric surfaces, and further optimizing GPU kernels for large-scale datasets could make differentiable meshing faster, more robust, and easier to adopt across geometry processing, simulation, and machine learning workflows.

Introduction and Motivation

One of the biggest bottlenecks in the modern model training pipelines for neuroimaging is diverse, well-labeled data. Real clinical MRI scans that could have been used, are under constraints of privacy, scarce annotations, and chaotic protocols. Because open-source MRI data comes from many scanners with different setups, the data varies a lot. This variation causes domain shifts and makes generalization much harder. All of these problems reinforce us to look for other solutions. Synthetic data directly handles these bottlenecks, once we have a robust pipeline that can generate unlimited, labelled data, we can expose models to the full range of scans (contrasts, noise levels, artifacts, orientations). With synthetic data, we are taking a step forward of solving the biggest bottlenecks in model training.

Nevertheless, another important thing is resolution. Most MRI scans are relatively low-resolution. At this scale, thin brain structures blur and different tissues can get mixed in a single voxel, so edges look soft and details disappear. This especially hurts areas like the cortical ribbon, ventricle borders, and small nuclei.

When models are trained only on these blurry data, they learn what they see: oversmoothed segmentations and weak boundary localization. Even with strong architectures and losses, you can’t recover details that simply don’t exist in the initial data that was fed to the model.

We are tackling this problem by generating higher-resolution synthetic data from clean geometric surfaces. Training with these high-resolution, label-accurate volumes gives models a much clearer understanding of what boundaries look like.

From MRI images to Meshes

The generation of synthetic data begins by using an already available brain MRI scan and replicating it into a collection of MRI data such that all the files in this collection represent the same MRI scan. However, each file uses a different set of contrasts to represent different regions of the brain. This ensures that our data collection is representative of the different models of MRI scanning devices available as each device can use a specific set of contrasts different from the ones used by other devices. Furthermore, we can make the orientation of the brain random to allow the data to mimic the randomness of how a person may position their head inside the MRI device. This process can be done easily using the already implemented python tool SynthSeg.

Using SynthSeg

By running the SynthSeg Brain Generator using an available MRI scan, we can get random synthesized MRI data. For illustration, we can generate 3 samples that have different contrasts but the same orientation by setting the number of channels to 3. Here is a cross section of the output of these samples.

The output of the Brain Generator for each instance are two files: an image file and a label file. The image file is an MRI data file that contains all of the scan’s details. On the other hand, the label file is a file that only indicates the boundaries of different structures indicated in the image file.

Mesh Generation

The label files can be useful as they only contain data describing the boundaries of each different region in the MRI data. This makes them more suitable for brain surface mesh generation. Using the SimNIBS package, different meshes can be generated via a label file. Each mesh represents a different region as dictated by the label file. After that, these meshes can be combined together to get a full mesh of the surface of the brain. Here are two meshes generated from the same label file.

From Meshes to MRI Data

Once the MRI scans are transformed into meshes, the next step is to convert these geometric representations back into volumetric data. This process, called voxelization, allows us to represent the mesh on a regular 3D grid (voxels), which closely resembles the structure of MRI images.

Voxelization Process

Voxelization works by sampling the 3D mesh into a uniform grid, where each voxel stores whether it lies inside or outside the mesh surface. This results in a binary (or occupancy) grid that can later be enriched with additional attributes such as intensity or material properties.

Input: A mesh (constructed from MRI segmentation or surface extraction).

Output: A voxel grid representing the same anatomical structure.

Resolution and Grid Density

The quality of the voxelized data is heavily dependent on the resolution of the grid:

• Low-resolution grids are computationally cheaper but may lose fine structural details.
• High-resolution grids capture detailed geometry, making them particularly valuable for training machine learning models that require precise anatomical features.

Bridging Mesh and MRI Data

By voxelizing the meshes at higher resolutions than the original MRI scans, we can synthesize volumetric data that goes beyond the limitations of the original acquisition device. This synthetic MRI-like data can then be fed into training pipelines, augmenting the available dataset with enhanced spatial fidelity.

Conclusion

In conclusion, this pipeline can produce high resolution synthetic brain MRI data that can be used for various purposes. The results obtained from this process are representative as the contrasts of the produced MRI data are randomized to mimic the different types of MRI machines. Due to this, the pipeline has the potential to sustain large amounts of data suitable for the training of various models.

Written By

• Youssef Ayman
• Ualibyek Nurgulan
• Daryna Hrybchuk

Special thanks to Dr. Karthik Gopinath for guiding us through this project!

SGI Fellows: Melisa Oshaini Arukgoda, David Uzor, Lydia Madamopouloua, Joana Owusu-Appiah
Mentors: Oliver Gross, Quentin Becker
SGI Volunteer: Denisse Garnica

In Week 4 of SGI, we worked on a simulation project that built on the previous work of our Mentors Oliver Gross and Quentin Becker: motion from shape change and inverse geometric locomotion – to put it simply, identifying the sequence of shapes (gait) that a body will need to morph into in order to follow a specific path, without the influence of external forces. In nature, we see this behaviour in many different organisms: Stingrays, snails, Cats (when adjusting their position to allow them to land on their feet after falling), and snakes. Unlike standard physics-based simulations that rely on heavy Newtonian calculations of forces and accelerations, the approach we worked with reformulated the problem as an optimization over shape changes alone. By leveraging conservation laws through geometric mechanics, we could avoid full force-based simulations and instead compute efficient gaits directly — a shift that makes the process both faster and more generalizable.

Motion from Shape Change

Snake-like robots are fascinating because they don’t rely on wheels or legs for movement. Instead, they move by changing their shape over time. This type of motion is called geometric locomotion. By coordinating joint angles in a specific sequence, the robot can “slither” forward, turn, or even back up. Mathematically, we describe this by mapping changes in shape (joint angles) to displacements in space. The beauty of this approach is that the same principles apply whether the robot is operating on flat ground or navigating through more complex environments.

To simulate this, we used Python code that defines the joint angles, applies forward integration, and produces the resulting trajectory. By running sequences of positive and negative angles, we could generate gaits (movement patterns) and measure their effectiveness. Parameters like amplitude, phase shift, and anisotropy ratio all influence how efficient or smooth the motion becomes.

Inverse Geometric Locomotion

While forward motion planning tells us what happens when we execute a gait, inverse geometric locomotion helps us answer the opposite question: what shape sequence do we need to reach a target position or avoid an obstacle? This is formulated as an optimization problem. We define objectives such as “reach this point,” “avoid that obstacle,” or “turn left”, and then optimize the gait until the resulting motion matches the objective.

Here, the gradient descent optimization algorithm is applied to adjust the sequence of joint angles. Our experiments demonstrated that even when constraints are applied (such as broken joints or limited power), the robot can often still achieve meaningful locomotion by adapting its gait.

Breaking and Fixing Robots

One of the challenges with snake robots is robustness. What happens if a joint breaks? We simulated this by disabling certain joints and observing the effect. Unsurprisingly, the robot’s gait efficiency drops, but not all joints are equally important. Some failures severely reduce mobility, while others can be compensated for.

To “fix” the robot, we rerun the optimization process under the new constraints. By adapting the gait, the robot can still move, albeit with reduced performance. This kind of resilience is critical for real-world applications, where hardware failures are inevitable.

Avoiding Obstacles

Another core part of our experiments was obstacle avoidance. Instead of using explicit geometry, we represented obstacles as implicit functions (mathematical descriptions where the function is positive outside, negative inside, and zero at the boundary). This allows smooth and flexible definitions of obstacles of any shape.

The robot’s motion planner then ensures that the trajectory stays in the safe (positive) region. We tested this with different obstacle shapes and verified that the snake robot can successfully reconfigure its gait to slither around barriers.

Applications

This computational framework for inverse geometric locomotion opens doors for several interesting applications, particularly in the fields of animation and robotics. The video below shows the snake robot in action (courtesy of the MiniRo Lab at the University of Notre Dame), loaded with the shape sequences from the simulation.

Thanks to Ioannis Mirtsopoulos, Elshadai Tegegn, Diana Avila, Jorge Rico Esteban, and Tinsae Tadesse for invaluable assistance with this article.

This article is a continuation of A Derivation of the Force Density Method, and will be continued by another article outlining results and next steps. Recall that we have constructed a function \[F:\mathscr Q \to \mathscr C,\] where \(\mathscr Q\) is the space of all possible force-length ratios and \(\mathscr C\) the space of all possible lists of coordinate pairs, such that the image of this function is the set of points in static equilibrium. Recall also that this function is a rational function (a ratio of two polynomial functions) with denominator given by a determinant. The goal of this article is to compute a simple description of the image of this function; in particular, we would like to compute polynomial equations which define the image.

Elimination Ideals and Groebner Bases

The technique we’ll use to compute the polynomial equations which define the image of this function is called an elimination ideal. Recall that the graph of a (polynomial) function \(f:\mathbb R^n\to \mathbb R^m\) is the set \(\{(x, f(x))\subset \mathbb R^n\times \mathbb R^m\}\). Note that there is an inclusion function \(\iota:x\mapsto (x, f(x))\) from \(\mathbb R^n\) to the graph of \(f\), and a projection function \(\pi\) from the graph of \(f\) to \(\mathbb R^m\) given by \((x, y)\mapsto y\), such that \(\pi\circ\iota = f\). Elimination theory first computes the polynomial equations for the graph of \(f\), and then computes the image of this graph under the projection function, thus computing the polynomial equations defining the image of \(f\).

In fact, this explanation has missed an important nuance: not every set can be defined by polynomial equations. If, however, the set cannot be defined by polynomial equations, this method will compute the smallest set which can be defined by polynomial equations and contains the image of \(f\).

The theory of Groebner bases and elimination ideals is rich and far exceeds the scope of this article; Eisenbud’s Commutative Algebra with a View Towards Algebraic Geometry is an excellent resource, in particular Chapter 15 on Groebner Bases, Section 15.10.4 being the relevant resource for this discussion. I will endeavor to sketch the techniques we use for this computation in layman’s terms, but it will be difficult to obtain a technical understanding of this theory without background in commutative algebra.

Localization

Recall first that the function \(F\) is a rational function, that is, it has a polynomial denominator and is thus only defined away from where that denominator attains the value zero. We need a way to ‘fix’ this, and the appropriate tool from commutative algebra is called localization. At a high level, localization is the process of formally inverting a polynomial, and considering the set of polynomials with that formally inverted polynomial. This not only allows us to work with this rational function, but also speeds up our computations down the line. There is a corresponding formal process which ensures that we only keep track of the solutions for this entire system which lie away from the set where that inverted polynomial is zero. For more information about this, look into the theory of affine schemes; the idea is that all of the information about the geometry we want to study is contained in the set of polynomials (called a ‘ring’) which are defined on that geometry. By inverting one of these polynomials, we remove the set where that polynomial is zero from the geometry of our space.

The Graph of \(F\) and projection

Computing the graph of the map \(F\) is not difficult; we simply take the set defined by polynomials of the form \(x_i – F_j\), where \(F_j\) is the polynomial defining the \(j\)th component of \(F\). Computing the image of the projection is more difficult; we must find a way of eliminating the first variables \(x_i\). To do this, we fix an ordering on the monomials of our polynomials, and use that ordering to perform a type of polynomial division – similar to the standard kind in one variable, but a little more complicated. For more information, see (for example) Chapter 15 of Eisenbud or Chapter 9.6 in Dummit and Foote’s Abstract Algebra. This leaves us with a set of polynomials, some of which will be in only the variables we want. It is a theorem (Proposition 15.29 in Eisenbud) that, so long as we pick our ordering correctly (it must be the case that whenever the first-ordered term in a polynomial contains only the variables we want to project onto, the entire polynomial contains only the variables we want to project onto), picking only the polynomials in the variables we want will give us polynomials who’s zero set contains every zero set of polynomials containing the image of this function.

Additional Thoughts

It’s important to realize that the result of this computation is not actually the image of the FDM map. This is a good thing, because the image of the FDM map cannot contain any values where nodes lying on the same line overlap, because this would result in a division by zero when computing force-length ratios. However, these are perfectly physical configurations, represented by a smaller network, and we want to be able to model them easily. The way to formally fix this is to introduce, instead of a single real variable for each force-length ratio, a point along the projective line. This allows the length to take the value zero, representing a force-length ratio of ‘infinity’, and yields a parameter space which consists of a copy of the projective line for every edge in the graph. This variety would be what is called a projective variety, for technical reasons; this would also imply that the variety was proper, in particular universally closed, another technical condition which implies that the projection map would take the zero sets of polynomials to the zero sets of polynomials (in the terms of algebraic geometry, the base-change of the map from the force-density space to the ground field along the map from the node-position space to the ground field would then be closed). It would be a little more complicated to apply our elimination theory method to this setup, but it would yield the assurance that the image of our map is exactly the set defined by the polynomials we obtain.

In a later blog post I will share the implementation of this algorithm in sage, and the results I obtained.