Differential Walk on Spheres

• Differential Walk on Spheres is a mesh-free, Monte Carlo method for estimating PDE solution sensitivities with respect to both geometric and boundary parameters.
• It employs a pathwise, unbiased estimator that efficiently couples primal and derivative computations without requiring grid discretization.
• The method scales efficiently for high-dimensional inverse problems and shape optimization, handling complex, mixed boundary conditions via star-shaped domains.

Differential Walk on Spheres (Differential WoS) is a mesh-free, Monte Carlo method for efficiently estimating the derivatives of solutions to partial differential equations (PDEs) with respect to problem parameters, including both geometric domain deformation and boundary condition variation. Like the classical Walk on Spheres (WoS) algorithm, Differential WoS avoids volumetric solvers and grid discretization, enabling scalable, pointwise evaluation of high-dimensional gradients in complex and topologically dynamic geometries. Its major innovation lies in a pathwise, expectation-based estimator for PDE parameter sensitivity that is unbiased and supports joint estimation with respect to all parameters at cost nearly independent of parameter count, making it well suited for PDE-constrained inverse problems, shape optimization, and training scenarios in computational science and graphics (Miller et al., 2024).

1. Mathematical Formulation and Target Problems

Differential WoS is primarily applied to screened Poisson boundary value problems on parameterized domains $\Omega(\pi)\subset\mathbb{R}^3$ with Dirichlet or mixed boundary conditions: $$\Delta u(x) - \lambda u(x) = f(x), \quad x\in\Omega(\pi)$$

$u(x) = g(x; \pi), \quad x\in\partial\Omega(\pi)$

Here, $\lambda\geq 0$ is the screening coefficient, $f$ is a source term, and $g$ may depend on a general parameter vector $\pi\in\mathbb{R}^N$ encoding geometry (such as mesh vertices, implicit surface controls) or boundary data (e.g., color, temperature, etc.). The formulation naturally accommodates parameters controlling deformations (pose, scale, translation) and boundary values, as well as problems with Robin or Neumann boundaries. For the Dirichlet–Robin case,

$$\frac{\partial u}{\partial n} - r u = h \text{ on } \partial\Omega_N$$

with Dirichlet data on $\partial\Omega_D$ (Miller et al., 2024).

2. Probabilistic Representation and Walk-on-Spheres Fundamentals

WoS exploits a probabilistic representation of the PDE solution via Brownian motion: $$u(x) = \mathbb{E}_x\left[ g(X_\tau) \right] + \mathbb{E}_x\left[ \int_0^{\tau} f(X_t) dt \right]$$ where $X_t$ is killed Brownian motion with killing rate $\lambda$ and exit time $\tau$ from $\Omega(\pi)$. The pathwise observables—boundary data and accumulated source terms—can be sampled by discrete random walks, hopping between spheres entirely contained within $\Omega(\pi)$. At each step, the radius is determined by distance to the boundary, and the next location is sampled uniformly on the corresponding sphere. Russian roulette strategies handle exponential kill for $\lambda > 0$. Implementation is agnostic to how the boundary is represented, requiring only distance and projection queries (Miller et al., 2024).

The "Walk on Stars" (WoSt) generalization further supports arbitrary mixed boundary operators by using star-shaped walk regions and importance sampling of Neumann/Robin boundary contributions, reflecting first-order (differential) geometry in the sampling process (Sawhney et al., 2023).

3. Pathwise Differentiation and the Differential PDE

Shape derivative theory shows that the derivative $u'(x) = \frac{\partial u}{\partial\pi}(x;\pi)$ itself solves a PDE: $$\Delta u'(x) - \lambda u'(x) = 0,\quad x\in\Omega(\pi)$$

$$u'(x) = g'(x;\pi) - v_n(x;\pi)\frac{\partial u}{\partial n}(x), \quad x\in\partial\Omega(\pi)$$

where $g'$ is the parameter derivative of the boundary data, and $v_n(x;\pi)$ is the boundary normal velocity with respect to $\pi$ (Miller et al., 2024). The boundary condition, involving the normal derivative of $u$, indicates that the parameter sensitivity can be represented as

$$u'(x) = \mathbb{E}_x\left[ g'(X_\tau) - v_n(X_\tau)\frac{\partial u}{\partial n}(X_\tau) \right]$$

which requires pathwise coupling of the primal variable ($u$) and its derivative, and accurate estimation of normal derivatives at the boundary.

The WoSt algorithm simulates reflecting Brownian motion and handles Neumann/Robin contributions via star-shaped domains, using silhouette distances and boundary normals to capture the necessary first-order information (Sawhney et al., 2023).

4. Differential Walk on Spheres Estimator

The Differential WoS estimator operates recursively, combining primal WoS estimation with finite-difference-like estimation of boundary-normal derivatives: $$\widehat{u'}(x_k) = \begin{cases} g'(x_k) - \frac{v_n(x_k)}{\Delta}\left[g(x_k) - \hat{u}(x_k - \Delta n)\right] & d_k < \epsilon, \ \alpha_k\,\widehat{u'}(x_{k+1}) & \text{otherwise} \end{cases}$$ where $n$ is the unit outward normal, $\Delta$ is a small, fixed offset for backward-difference approximation, and $d_k$ is the distance to the boundary at step $k$ (Miller et al., 2024). Primal estimation $\hat u$ involves its own, potentially nested, WoS walk. Both $u'(x)$ and the full Jacobian with respect to $\pi$ are produced in a single Monte Carlo walk, providing efficient, simultaneous parameter sensitivity evaluation.

Bias in estimation comes from the finite $\epsilon$-shell termination (first exit within $\epsilon$ of the boundary) and finite-difference step $\Delta$; both can be decreased to tradeoff bias and computational cost.

In WoSt, the estimator incorporates star-shaped domains and samples Neumann or Robin boundary contributions explicitly. Star regions ensure the unbiasedness of the estimator in domains with complex (including non-convex) or mixed boundary conditions (Sawhney et al., 2023).

5. Theoretical Properties: Error, Complexity, and Scaling

The estimator is unbiased as $\epsilon, \Delta \to 0$, with statistical variance decaying as $O(1/M)$ for $M$ independent samples, leading to root-mean-square error $O(M^{-1/2})$. Bias introduced by $\epsilon$ and $\Delta$ is linear for sufficiently regular solutions. The overall mean-square error is the sum of squared bias and variance, and can be tuned by adjusting these parameters (Miller et al., 2024).

Key computational properties are summarized as follows:

Property	Primal WoS/Differential WoS	Consequence
Mesh/grid-free	Yes	Boundary-agnostic, no remeshing
Parameter scalability	Cost $\approx$ constant in $N$	Efficient for high-dimensional $\pi$
Parallelism	All walks independent	Trivial data parallelism (CPU/GPU)
Output sensitivity	Evaluates at selected points only	No global volumetric solve

Each sample step requires $O(\log M)$ time (with $M$ boundary elements) if using spatial acceleration structures for queries. For Dirichlet-dominated problems, the expected walk length is $O(\log(1/\epsilon))$; for Neumann-dominated with regularization, $O(1/\sqrt{\sigma})$, where $\sigma$ is the Neumann Tikhonov parameter (Sawhney et al., 2023).

6. Practical Implementation Aspects

Differential WoS and WoSt accept any boundary representation as long as required geometric queries are supported. Efficient implementations employ acceleration (e.g., BVH, ray tracing) for distance and normal queries. Walkers may use meshes, splines, or implicit surfaces, facilitating model interchange in geometry-rich domains.

Storing derivatives per sample requires $O(N)$ memory per walk, which is manageable even for large $N$ compared to storing a Jacobian in classical numerical PDE solvers. Unlike finite-difference-based PDE sensitivity, where each parameter requires a distinct solve, a single Differential WoS walk provides the entire derivative vector (Miller et al., 2024).

WoSt requires silhouette computation and normal cones stored within spatial hierarchies to determine star-shaped subdomains and facilitate Neumann and Robin boundary sampling (Sawhney et al., 2023).

7. Application Domains and Illustrative Examples

Differential WoS enables inverse optimization and sensitivity analysis in scenarios where traditional solvers are computationally prohibitive or ill-suited for complex geometries. Notable application classes include:

• Pose estimation: $\pi$ encoding SE(3) transformations, optimizing object pose from diffusive image data.
• Shape-from-diffusion: reconstructing geometry from external or boundary measurements, with $\pi$ controlling mesh vertices; employs stochastic preconditioners (e.g., Laplacian smoothing).
• Thermal design with mixed boundaries: optimal heating element geometry, solved by differential WoSt for mixed Dirichlet–Neumann BVPs, leveraging deferred shading for partial image loss.
• Image-space curve inflation: inverse Poisson inpainting, $\pi$ controlling boundaries via RBFs; boundary integrals addressed by $\epsilon$-shell fattening.
• Inverse diffusion curves: joint optimization of shape and boundary data in diffusion-curve images, with regularizers estimated via closed-form Monte Carlo integrals (Miller et al., 2024).

The parameter independence of sampling cost and output sensitivity (ability to evaluate only at points of interest) are central advantages in high-resolution or geometrically complex domains.

8. Differential WoS as a General Template

The walk on stars methodology generalizes the classical WoS approach by embedding first-order boundary differential operators (e.g., normal derivatives in Neumann or Robin BVPs) directly into the sampling process. Augmenting spatial hierarchies (e.g., BVHs) with normal cones discretizes boundary operators across scales. This strategy provides a grid-free Monte Carlo template capable of incorporating arbitrary linear boundary operators, via star-shaped domains, directional sampling, and structured kernel weighting for unbiased estimation (Sawhney et al., 2023).

In summary, Differential Walk on Spheres constitutes a robust, mesh-agnostic tool for computing PDE solution sensitivities with respect to arbitrary parameters. Its combination of parameter scalability, boundary-representation agnosticism, perfect parallelism, and unbiasedness positions it as a central building block for high-dimensional inverse problems and PDE-constrained optimization in scientific computing, computational geometry, and computer graphics (Miller et al., 2024, Sawhney et al., 2023).