Abstract
The Walk on Spheres (WoS) algorithm is a mesh-free and highly flexible Monte Carlo method for solving partial differential equations, but its practical applicability is limited by slow $O(N^{-1/2})$ convergence. While prior variance reduction techniques exploit spatial correlations through integral properties of the PDE, they do not fully utilize the intrinsic Markov structure of the WoS process. We introduce a new variance reduction framework based on reusing intermediate states along each random walk. Leveraging the Markov property, we show that every point visited by a WoS trajectory provides a valid unbiased estimator, but its direct use is hindered by the complex distribution induced by dynamically generated spheres. To resolve this, we propose the Walk on Probes (WoP) algorithm, which replaces dynamic spheres with a set of fixed, pre-distributed spherical probes inside the domain. This converts the intractable distribution of path points into samples on fixed boundaries, enabling efficient evaluation through the Poisson integral formula. We further develop a specialized method that combines control variates with self-normalization to further reduce variance. Together, these components substantially improve sample efficiency while preserving the flexibility of WoS.