Harmonic Caching for Walk on Spheres

We present a variance reduction technique for Walk on Spheres (WoS) that solves elliptic partial differential equations (PDEs) by combining overlapping harmonic expansions of the solution, each estimated using unbiased Monte Carlo random walks. Our method supports Laplace and screened Poisson equations with Dirichlet, Neumann, and Robin boundary conditions in both 2D and 3D. By adaptively covering the domain with local expansion regions and reconstructing the solution inside each region using an infinite Fourier series of the harmonic function, our method achieves over an order of magnitude lower error than traditional pointwise WoS in equal time. While low-order truncations of the series typically introduce limited bias, we also introduce a stochastic truncation scheme that eliminates this bias in the reconstructed solution. Compared to recently developed caching algorithms for WoS, such as Boundary and Mean Value Caching, our approach yields solutions with lower error and fewer correlation artifacts.