Planning under Noise with Neural Operators

Michael Chen, Hong Chul Nam, Xi Deng, Anima Anandkumar

arXiv:2610.28623v1 · This workspace stays on this version. Check arXiv for newer versions

Affiliations and author notes

Computing and Mathematical Sciences, California Institute of Technology
Center for Data Science, New York UniversityEqual contribution

Some source details need review 3 to check

Original formulas, table values and figures are preserved. Comparison can miss differences; use the original PDF to check details.

Reading view checked: 802 math expressions, 0 rendering failures; 18 of 18 images available.

Compared 251 of 288 source passages against 28 of 28 PDF pages.

Some graphics are laid out differently in the HTML source than in the paper (text or formulas outside their frame); check the original PDF. (1)
S1_F2: the arXiv HTML draws text or formulas of this graphic outside its frame and it could not be corrected safely (too-small); the image shows the HTML layout, which differs from the PDF.
This long paper was only partially compared; the original sources remain available.

18 of 18 figure, table and graphic images captured. Open original PDF

Abstract

Planning under uncertain control execution noise requires reliable feedback policies across diverse environments. We introduce Planning Under Noise (PUN), a neural operator that learns noise-aware value fields for planning in unseen geometries. Building on an established screened-Poisson formulation of stochastic optimal control, we address two challenges: specifying boundary conditions that support navigation across varying geometries and efficiently computing fields with reliable gradients. Zero-Dirichlet obstacle boundaries terminate sampled paths at walls, weakening the goal signal along long routes and making it harder to estimate reliably. We instead impose reflecting Neumann conditions on obstacles and a Dirichlet condition at the goal, allowing sampled paths to continue exploring free space until they reach the target. To solve the resulting PDE without a configuration-space mesh, we extend a Monte Carlo (MC) solver, Walk on Stars (WoSt), to higher dimensions. However, local MC fluctuations can reverse estimated gradient directions and cause planning failures, making direct planning require costly sampling. PUN amortizes this computation by learning geometry-to-field mappings from noisy, partially converged Monte Carlo labels, avoiding expensive convergence for every training example. At inference, test-time scaling uses bounded local search guided by the predicted global field. We evaluate PUN on planning tasks across dimensions. Across four noise levels, PUN matches or outperforms all compared planners, achieving 95–100% success on the 4D replanning benchmark while maintaining the lowest runtime. At the highest tested noise level in 4D ( $\lambda=4$ ), PUN achieves 10.7–60.6% relative gains in success rate over other methods while running 3.29–32.15 $\times$ faster.

Planning (Latombe, 1991) is one of the most fundamental problems in robotics. Given a robot’s configuration space (LaValle, 2006), planning computes a path from the current configuration to a goal that is optimal under a given cost function. Real robots do not execute a commanded motion exactly (Thrun, 2002). These errors accumulate along a motion, so a route planned with a few centimetres of clearance can be executed into a wall. Sampling-based planners such as PRM (Kavraki et al., 1996) and RRT (LaValle and Kuffner, 2001), and graph search such as A⋆ (Hart et al., 1968), plan for a nominal system that executes commands exactly and leave the noise to a tracking controller. Sampling-based methods such as MPPI (Williams et al., 2017a) and PI $^{2}$ -CMA (Stulp and Sigaud, 2012) optimize control sequences or policy parameters using sampled trajectory costs. Optimal feedback under execution noise can be derived from a value function, which gives the minimum expected cost of reaching the goal from each configuration and satisfies the nonlinear Hamilton–Jacobi–Bellman (HJB) equation (Fleming and Soner, 2006). Solving this equation on a grid scales exponentially with the dimension of the configuration space (Bellman, 1957; Han et al., 2018).

While the problem can be converted into a linear screened Poisson equation (Theodorou et al., 2010) under control–noise matching (Assumption 1), a central difficulty in Poisson-based planning is specifying boundary conditions that produce a reliable value field across varying geometries.

Fig. 1: **(a) Fast replanning.** PUN has the lowest mean runtime on unseen 4D maps: 10.8 ms/start, compared with 32.8–387.1 ms/start for the baselines. Runtimes are averaged equally over the four execution-noise levels in Table 5 and plotted on a logarithmic axis. **(b) Resolution transfer.** FNO maintains over 90% local-search success across all five resolutions without retraining, reaching 92.50% at $4096\times 4096$ after training at $256\times 256$ . Missing points indicate out-of-memory failures. **(c, d) Continued goal seeking.** In the illustrated geometry, the solution induced by the Dirichlet boundary cannot propagate goal information further since most trajectories are absorbed by obstacles. Consequently, it is more difficult for planners to find the target than with the mixed boundary condition.

Conventional mesh-based numerical methods can solve the screened Poisson equation efficiently in two and three dimensions. However, a robot with $d$ degrees of freedom (DoF) induces a $d$ -dimensional configuration space, making volumetric discretizations and numerical PDE solution increasingly expensive. Planning must also accommodate complex and previously unseen obstacle geometries, for which constructing a new mesh and solving the PDE from scratch can be costly. Mesh-free Monte Carlo methods (Muller, 1956; Sawhney and Crane, 2020; Sawhney et al., 2023; Jambon et al., 2026) support PDE solution on complex geometries, motivating their use for planning in higher-dimensional configuration spaces.

Among these, Walk-on-Spheres (WoS) (Muller, 1956) estimates the solution at a point by jumping repeatedly to a uniformly sampled point on the largest ball fitting in free space, until it reaches a boundary of known value, which requires every reachable boundary to be Dirichlet. Walk-on-Stars (WoSt) (Sawhney et al., 2023) accommodates reflecting obstacles by replacing that ball with a star-shaped region that may contain part of the obstacle boundary, keeping every boundary point visible from the center so that one ray continues the walk instead of branching at each contact. The screened kernels for Walk-on-Stars have been derived only in two and three dimensions (Sawhney et al., 2023). Moreover, directly planning with Monte Carlo estimates is unreliable: trajectories depend on local value gradients, and even modest estimator noise can corrupt or reverse the estimated direction, leading to planning failure.

Neural operators (Kovachki et al., 2023) can amortize expensive PDE solves by learning mappings from input geometries directly to solution fields. Moreover, prior work has shown that fully converged labels are unnecessary when training from Monte Carlo simulations: neural operators can recover accurate predictions from noisy, partially converged estimates (Viswanath et al., 2026). This enables training of neural operators with inexpensive Monte Carlo supervision while providing fast, smooth value-field predictions at inference time.

We propose Planning Under Noise (PUN), shown in Fig. 2, an operator-learning framework that maps environment geometry to value fields that account for control noise. Boundary conditions can be understood through a heat analogy. With unit desirability at the goal and zero desirability on obstacle boundaries, the goal acts like a hot source, while the obstacle walls act like cold sinks. The zero-Dirichlet condition fixes the field to zero along each obstacle surface, making the entire obstacle boundary a heat sink. In a long or narrow passage, heat spreading from the goal therefore leaks into the walls along the route. The farther the goal’s influence must travel, the more it is attenuated, leaving a weak signal in distant parts of free space that discretization, sampling, or learning errors can obscure. Narrow passages make this particularly pronounced because the absorbing walls remain close on both sides. This attenuation is visible in Fig. 1(c), where the weak goal signal provides insufficient directional guidance for starting points far from the goal, making it difficult to navigate toward the target. Moreover, the screened Poisson equation already includes decay throughout free space through its screening term; absorbing obstacle boundaries introduce an additional source of attenuation. We therefore keep a Dirichlet condition at the goal but impose homogeneous Neumann (zero-flux) conditions on obstacle boundaries. The walls now act like insulation: they prevent heat from escaping through the boundary without forcing the desirability itself to vanish there. Heat can spread along a corridor without continually draining into its walls, allowing the goal’s influence to reach farther into free space. Screening still causes decay, but the additional loss through obstacles is removed. This is precisely the effect illustrated in Fig. 1(d): compared with the absorbing Dirichlet field in Fig. 1(c), the mixed-boundary field extends the goal’s influence much farther through the environment, preserving a stronger signal for planning along long routes. The Neumann condition also makes the value gradient parallel to smooth obstacle boundaries; therefore, gradient-based planning directs motion along walls rather than into them. This permits guidance along walls in narrow passages without the strong repulsion induced by zero-Dirichlet conditions. Under our regularity assumptions, the value field and control remain bounded at obstacles, whereas zero-Dirichlet desirability makes both diverge there (Appendix C.4).

We first derive dimension-general Poisson kernels and Monte Carlo weights, extending mesh-free supervision beyond three dimensions for scalability. This enables supervision generation without a volumetric configuration-space mesh. We then train PUN using noisy, partially converged Monte Carlo estimates in the weak supervision paradigm. By amortizing the PDE solution across geometries, PUN produces spatially coherent value fields that guide gradient following and local search without requiring a fully converged Monte Carlo solve for every environment. We evaluate PUN across 2D–4D environments and compare its planning success and runtime with nominal planners and stochastic controllers under execution noise. In the experiments varying the execution noise level $\lambda$ , we plan a path from the current state at each control step, follow it briefly, and repeat from the updated state; we refer to this as replanning. Our contributions can be summarized as follows:

•
We formulate Poisson-based path planning as a mixed boundary-value problem, imposing a Dirichlet condition at the goal and homogeneous Neumann conditions on obstacle boundaries.
•
We introduce an operator-learning framework that maps environment geometry to noise-aware value fields, enabling zero-shot planning in unseen complex geometries without per-instance retraining or numerical PDE solves.
•
We generalize the mesh-free Monte Carlo solver WoSt beyond three dimensions by deriving the dimension-dependent Green’s function and step weights for the screened Poisson equation.
•
We introduce a test-time scaling strategy that evaluates reachable candidate states on a local lattice using a fixed global learned field, without retraining, improving planning success over direct gradient following in the reported 2D–4D evaluations (Tables 1 and 11).

Fig. 2: Planning Under Noise: **(A) Planning:** a neural operator maps geometry, goal, cost, and noise inputs to a desirability field. Grid-based local search selects feasible routes with high endpoint desirability, executes a short prefix, and replans. **(B) Training:** Walk-on-Stars (WoSt) estimates provide pointwise supervision; neural caching supplies learned continuation at the walk horizon.

On 2D application benchmarks, zero-shot WoSt-supervised FNO reaches 93.62% local-search success, up from 83.87% with gradient following on the same fields (Table 1). With the same local-search, FNO-based PUN achieves approximately $941\times$ and $58\times$ speedups over direct WoSt on the 2D zero-shot application benchmarks and 4D manipulator test set, respectively (Table 12). In 4D replanning, PUN achieves 95.00–100% success at 10.51–11.66 ms/start across the tested noise levels (Table 5; Fig. 1(a)). At $\lambda=4$ , PUN is $3.29\times$ faster than RRT and $5.09$ – $32.15\times$ faster than A⋆, Dijkstra, MPPI, and PI $^{2}$ -CMA. FNO also maintains over 90% local-search success across resolutions up to $4096\times 4096$ ( $256\times$ as many grid points as used in training) without retraining (Fig. 1(b)).

In this section we derive the elliptic boundary value problem that our operator learns, and fix the notation used throughout.

Let $\Omega\subset\mathbb{R}^{d}$ and let $\xi\in\overline{\Omega}$ be the starting point. The controlled process $X^{\xi,\mathbf{u}}=(X_{t}^{\xi,\mathbf{u}})_{t\geq 0}$ has interior dynamics

\mathrm{d}X^{\xi,\mathbf{u}}_{t}=\big(\mathbf{b}(X^{\xi,\mathbf{u}}_{t})+G(X^{\xi,\mathbf{u}}_{t})\,\mathbf{u}_{t}\big)\,\mathrm{d}t+\Sigma^{1/2}(X^{\xi,\mathbf{u}}_{t})\,\mathrm{d}\mathbf{W}_{t},\qquad X^{\xi,\mathbf{u}}_{0}=\xi,

with control $\mathbf{u}_{t}\in\mathbb{R}^{k}$ and standard Brownian motion $\mathbf{W}$ ; the process reflects normally at obstacles and workspace boundaries ( $\Gamma_{N}$ ) and stops at the goal ( $\Gamma_{D}$ ). We minimize

J(x,\mathbf{u})=\mathbb{E}\left[\int_{0}^{\tau_{D}^{\mathbf{u}}}\Big(q(X^{\xi,\mathbf{u}}_{t})+\tfrac{1}{2}\mathbf{u}_{t}^{\top}R\,\mathbf{u}_{t}\Big)\mathrm{d}t\;+\;\phi(X^{\xi,\mathbf{u}}_{\tau_{D}^{\mathbf{u}}})\middle|\xi=x\right],

where $\tau_{D}^{\mathbf{u}}:=\inf\{t\geq 0:X_{t}^{\xi,\mathbf{u}}\in\Gamma_{D}\}$ is the first goal hit, $q\geq 0$ , $R\succ 0$ , and $\phi\geq 0$ is the terminal cost. Write $g=e^{-\phi/\lambda}$ for the Dirichlet desirability; operator learning fixes $\phi\equiv 0$ , hence $g\equiv 1$ , while the PDE identities retain general positive $g$ . Where the value $V(x)=\inf_{\mathbf{u}}J(x,\mathbf{u})$ is $C^{2}$ it satisfies the stationary HJB equation, which is nonlinear: eliminating the optimal control leaves a term quadratic in $\nabla V$ (Appendix C.2, Eq. 29).

Assumption 1 (Control–noise matching).

There exists $\lambda>0$ with $\Sigma(x)=\lambda\,G(x)R^{-1}G(x)^{\top}$ for all $x\in\Omega$ .

This is the standard condition of linearly solvable stochastic optimal control (Kappen, 2005; Todorov, 2006; Todorov, 2009; Theodorou et al., 2010): a direction that is cheap to actuate is also one along which the system is uncertain, and $\lambda$ sets the exchange rate between the two. This matching condition enables the exact linearization of the HJB equation used here. The Cole–Hopf change of variable $\Psi(x)\triangleq\exp(-V(x)/\lambda)$ (Hopf, 1950; Cole, 1951; Fleming, 1977) cancels the quadratic term and turns the HJB equation into a linear elliptic equation for $\Psi$ , as shown in Appendix C.2.

For a fully actuated, passive system we take $\mathbf{b}=0$ , $G=I$ , and $R=I$ , so Assumption 1 forces $\Sigma=\lambda I$ . Writing the state cost as $q(x)=\tfrac{1}{2}c(x)^{2}$ with $c>0$ , the linearized equation becomes the screened Poisson equation:

\Delta\Psi-\kappa^{2}(x)\,\Psi=0,\qquad\kappa(x)\triangleq\frac{c(x)}{\lambda}.

Throughout, we take $c$ constant within each environment, so $\kappa$ is constant.

Let $D$ be the workspace, $\mathcal{S}$ the obstacle interiors, and $\mathcal{G}$ an open goal region. The free domain is $\Omega=D\setminus(\overline{\mathcal{S}}\cup\overline{\mathcal{G}})$ , assumed connected, with boundary partitioned into $\Gamma_{D}=\partial\mathcal{G}$ and $\Gamma_{N}=\partial\mathcal{S}\cup\partial D$ . We impose

\Psi=g\ \text{ on }\Gamma_{D}\qquad\text{(absorbing goal)},\qquad\partial_{n}\Psi=0\ \text{ on }\Gamma_{N}\qquad\text{(reflecting)}.

The goal is the only Dirichlet set; obstacles and the outer domain boundary are reflecting in the PDE model. Under Assumption 2 the mixed problem has a unique solution with $\Psi\in(0,1]$ and $\inf_{\overline{\Omega}}\Psi>0$ , so $\log\Psi$ is bounded and the log parameterization is well defined (Proposition 3). See Appendix C.4 for a detailed analysis.

A neural operator learns a map from problem data to a solution field (Li et al., 2021). Here we learn the desirability operator

\mathcal{T}:\big(\ \underbrace{\mathbf{1}_{\mathcal{S}}}_{\text{obstacles}},\ \underbrace{\mathbf{1}_{\mathcal{G}}}_{\text{goal}},\ \underbrace{c}_{\text{cost}},\ \underbrace{\lambda}_{\text{noise}}\big)\;\longmapsto\;\Psi,

Both geometric inputs are masks on the workspace: $\mathbf{1}_{\mathcal{S}}$ is the occupancy map of the obstacles and $\mathbf{1}_{\mathcal{G}}$ marks the goal region; goal data $g\equiv 1$ is fixed. A single trained model $\mathcal{T}_{\theta}$ predicts the field for new environments. Its spatial gradients provide control through Eq. 11, and its values support local search (Section 4.4).

Fig. 3: Replanning trajectories on the same 4D manipulator, start, and goal at $\lambda=0.25$ (top) and $\lambda=2$ (bottom). Columns show successive execution states; PUN’s trajectory is orange. PUN reaches the goal in both illustrated rollouts. Aggregate success and runtime are reported in Table 5.

Classical planners divide into two families. Sampling-based methods such as PRM (Kavraki et al., 1996), RRT (LaValle and Kuffner, 2001), and RRT⋆ (Karaman and Frazzoli, 2011) build a graph of collision-free motions, while search-based methods such as A⋆ (Hart et al., 1968), Dijkstra (Dijkstra, 1959), and fast marching (Sethian, 1996) compute a cost-to-go over a discretization; both plan for a deterministic system and defer execution noise to tracking and replanning. Sampling-based methods such as MPPI (Williams et al., 2017a; Williams et al., 2017b) and PI $^{2}$ -CMA (Stulp and Sigaud, 2012) optimize control sequences or policy parameters using sampled trajectory costs. Our implementations and their treatment of execution noise are described in Appendix D.1. We compare against these planners in Table 5. Planning Neural Operator (PNO) (Matada et al., 2025) learns the solution operator of the Eikonal equation from labels produced by a deterministic solver, but it is limited to deterministic systems. NTFields (Ni and Qureshi, 2023) and EikoNet (Smith et al., 2021) train a physics-informed network on the Eikonal residual, which fits a single geometry per model and must be retrained for each new environment. In contrast, PUN is trained over a family of geometry classes with approximate Monte Carlo supervision across diverse noise ranges to support planning with different noise spectra while achieving reliable stability.

Monte Carlo estimators for PDEs can be categorized depending on the type of PDE, include Walk-on-Spheres (Sawhney and Crane, 2020), Walk-on-Boundaries (Sugimoto et al., 2023), and Backward Stochastic Differential Equations (BSDE) (Zhang, 2017). Recent works combine Monte Carlo methods with neural networks to further reduce variance (Nam et al., 2024; Li et al., 2023a), which train a network by regressing on noisy Monte Carlo targets and thereby amortize repeated pointwise Monte Carlo estimation into a forward pass. Two of their ingredients we adopt directly: a deterministic cap on the number of walk steps, and evaluation of the partially trained model when a walk exceeds that cap. The diffusion loss (Nüsken and Richter, 2021) truncates an SDE simulation and supplements it with a boundary physics loss.

Neural operators parameterize maps between function spaces and can be evaluated on different discretizations, with numerical accuracy depending on the representation (Kovachki et al., 2023). The architectures we compare differ in how they propagate information: FNO (Li et al., 2021), UNO (Rahman et al., 2023), LocalNO (Liu-Schiaffini et al., 2024), DISCO (Ocampo et al., 2023), GINO (Li et al., 2023b), UPT (Alkin et al., 2024) and Transolver (Wu et al., 2024). We include U-Net (Ronneberger et al., 2015) as a fixed-resolution reference. Physics-informed operators such as PINO (Li et al., 2024) add a PDE residual to the objective; we use walk labels as the primary signal with regularization on gradient terms.

PUN learns a desirability field from Walk-on-Stars supervision and uses the predicted field to guide local search. Labels are constructed directly in $\Psi$ -space using finite-horizon walks, with learned continuation for unfinished trajectories. The neural operator is trained online on these labels. Algorithms 1 and 2 describe a single walk and the training loop, respectively; walks and terminal evaluations are batched.

Let $\sigma=\inf\{k:\xi_{k}\in\Gamma_{D}\}$ be the first exact goal-boundary hit, possibly infinite. The ideal chain stops at an exact hit, where the terminal value is $g(\xi_{\sigma})$ ; practical walks stop in a goal shell or at a finite horizon. For a spatial location $x$ , the star radius is

\begin{gathered}r(x)=\min\{d_{D}(x),\,d_{\mathrm{sil}}(x)\},\\ \mathrm{St}(x)=\text{the adjacent component of }B(x,r(x))\cap\Omega.\end{gathered}

Here $d_{\mathrm{sil}}$ is the distance to the closest visibility silhouette of $\Gamma_{N}$ . An admissible star step additionally requires that this component be star-shaped about $x$ and that its boundary consist only of the sphere and $\Gamma_{N}$ . Appendix D.1 describes the shell stopping rule and radius floor used in place of exact absorption.

Write $\xi_{k}$ for the WoSt position after $k$ steps, with $\xi_{0}=\xi$ the starting point and $\xi_{k}$ random for $k\geq 1$ . Given $\xi_{k}$ , construct $\mathrm{St}(\xi_{k})$ with radius $r_{k}=r(\xi_{k})$ . Draw a direction uniformly on the sphere in the interior or on the inward hemisphere at a smooth reflecting boundary point, then move to the first positive hit of the sphere or an obstacle/domain boundary:

\mathbf{v}_{k}\sim\mathrm{Unif}(\mathbb{S}^{d-1}),\qquad\rho_{k}=\text{first-hit distance},\qquad\xi_{k+1}=\xi_{k}+\rho_{k}\mathbf{v}_{k}.

The hemisphere restriction replaces the first distribution at $\xi_{k}\in\Gamma_{N}$ . We assume measurable geometry queries and positive hit distances almost surely.

On an admissible star, the direction law has surface density

p(x,\mathrm{d}\mathbf{z})=\frac{n_{\mathbf{z}}\cdot(\mathbf{z}-x)}{\alpha(x)\omega_{d}\|\mathbf{z}-x\|^{d}}\,\mathrm{d}S(\mathbf{z}),\qquad\int_{\partial\mathrm{St}(x)}p(x,\mathrm{d}\mathbf{z})=1,

where $\omega_{d}=2\pi^{d/2}/\Gamma(d/2)$ and $\alpha=1$ in the interior or $1/2$ at a smooth boundary pole. Tangencies have zero density. Where a step is not admissible, Eq. 7 still defines a law on first hits, but Eq. 8 no longer describes it.

Screening discounts each step by $W_{k}=Q_{d}^{\kappa}(\rho_{k},r_{k})\in(0,1)$ , given in closed form for every $d$ in Appendix C.5. For constant $\kappa$ , Green’s identity on each admissible active star gives

\mathbb{E}[W_{k}\Psi(\xi_{k+1})\mid\mathcal{F}_{k}]=\Psi(\xi_{k}),\qquad k<\sigma.

Multiplying by the preceding weights makes $M_{k}=(\prod_{j<k\wedge\sigma}W_{j})\Psi(\xi_{k\wedge\sigma})$ a bounded martingale. This gives the finite-horizon identity below without requiring an exact goal hit in finitely many steps. The identity concerns the weighted ray sampler; its steps need not follow Brownian exit distributions.

For a finite walk cap $K$ , write $\tau=\sigma\wedge K$ . Repeated conditional expectation in Eq. 9, keeping the chain fixed after a goal hit, gives

\mathbb{E}_{x}\!\left[\Big(\prod_{j<\tau}W_{j}\Big)\Psi(\xi_{\tau})\middle|\xi_{0}=x\right]=\Psi(x).

The terminal factor is $g(\xi_{\sigma})$ on a goal hit and the true continuation $\Psi(\xi_{K})$ otherwise. With true continuation and admissible steps, averaging $N$ independent replicas gives an unbiased label. Training substitutes learned continuation as described below. See Appendix C.6 for the log-space alternative and its bias correction.

Practical walks stop within a goal shell of width $\varepsilon$ and use a radius floor. Under the smooth-domain and admissibility assumptions, shell stopping with true continuation changes the mean by at most $C_{\Psi}\varepsilon$ , where $C_{\Psi}=\sup_{\overline{\Omega}}\|\nabla\Psi\|$ . The radius floor can violate admissibility and introduces an additional error of unknown sign (Sawhney et al., 2023). Appendix D.1 gives the implementation and the finite-horizon shell bound.

At the horizon $\Psi$ is not available, so training substitutes a detached EMA copy of the model. The resulting bias is that copy’s own error at the endpoint, weighted by the product of the step weights over the $K$ steps, so it is discounted wherever steps are long compared with the screening length.

For $\lambda>0$ and a positive differentiable prediction, the desirability field defines the feedback control through its logarithm: $^{1}$

\Psi_{\theta}=\mathcal{T}_{\theta}[a],\qquad V_{\theta}=-\lambda\log\Psi_{\theta},\qquad\mathbf{u}_{\theta}=\lambda R^{-1}G^{\top}\nabla\log\Psi_{\theta}.

For a coordinate-differentiable decoder, the gradient is obtained by automatic differentiation. A grid-output backbone requires a specified differentiable interpolant or a spatial differentiation rule. The input $a$ is defined in Eq. 5, with $g\equiv 1$ fixed.

We combine value and goal supervision with gradient loss to train fields for downstream planning:

\mathcal{L}_{\mathrm{val}}=\mathbb{E}\big[(\Psi_{\theta}-\widehat{\Psi})^{2}\big],\quad\mathcal{L}_{\mathcal{G}}=\mathbb{E}\big[(\Psi_{\theta}-g)^{2}\big],\quad\mathcal{L}_{\mathrm{grad}}=\frac{\mathbb{E}\big[\|\nabla\Psi_{\theta}-\nabla\widehat{\Psi}\|^{2}\big]}{\mathrm{sg}\big[\mathbb{E}\big[\|\nabla\widehat{\Psi}\|^{2}\big]\big]}.

The target gradient $\nabla\widehat{\Psi}$ is computed by finite differences of the label field. The displayed gradient loss normalization assumes positive target gradient energy.

The objective is their weighted sum,

\mathfrak{L}(\theta)=\mathcal{L}_{\mathrm{val}}+\beta_{\mathcal{G}}\,\mathcal{L}_{\mathcal{G}}+\beta_{\mathrm{grad}}\,\mathcal{L}_{\mathrm{grad}}.

The goal loss carries the exact boundary value $g$ rather than a walk label. The gradient term is included because a small value error need not give a small gradient error, and the control in Eq. 11 is computed from the gradient: $\mathcal{L}_{\mathrm{grad}}$ normalizes gradient error by target gradient energy, treating the denominator as fixed during differentiation. Direct- $\Psi$ supervision avoids log-transform bias and weights absolute field error, giving small absolute errors in low-desirability regions less weight than log-space supervision. Appendix C.6 reports the log-space alternative, in which errors are measured on a relative scale.

PUN supports test-time scaling by increasing the local search budget without modifying the learned global field. More depth lets the search evaluate more reachable endpoints, while using the field for guidance. The same endpoint search can be run on other fields; in PUN, it is run on a global noise-aware field that is reused across replans.

Paths are extracted from $\Psi_{\theta}$ on a lattice. Let $\mathcal{X}_{h}$ be the free cells at spacing $h$ . The feasible neighbors of $x$ are $\mathcal{N}(x)=\{x^{\prime}\in\mathcal{X}_{h}:\|x^{\prime}-x\|_{2}=h\}$ , the axis-adjacent free cells. Write $\mathcal{R}_{k}(x)$ for the states reachable from $x$ in at most $k$ such moves, including $x$ itself. A plan is a feasible lattice path $\pi=(x_{0},\ldots,x_{\ell})$ with $\ell\leq L$ , scored by the predicted desirability $\Psi_{\theta}(x_{\ell})$ at its endpoint. For $\lambda>0$ and $\Psi_{\theta}>0$ , maximizing this score is equivalent to minimizing $V_{\theta}=-\lambda\log\Psi_{\theta}$ at the endpoint. The score itself is well defined for any finite prediction.

A controller is specified by a depth $L$ and commit length $e\leq L$ . With one edge committed per replan, the next move is

x_{1}\in\operatorname*{arg\,max}_{y\in\mathcal{N}(x_{0})}\ \max_{z\in\mathcal{R}_{L-1}(y)}\Psi_{\theta}(z),

a step towards the largest desirability reachable through a candidate neighbor. For $e>1$ , a path maximizing endpoint desirability is constructed, its first $\min(e,\ell)$ edges are executed, and the search is repeated from the new state. Execution stops early upon reaching the goal. See Appendix D.1 for further details.

A key idea is that evaluating reachable endpoints permits an intermediate decrease in desirability when it leads to a better state, so local fluctuations need not determine the next move. This might happen when the field is not smooth and there are local valleys, so local search allows jumps over imperfect predictions. Appendix C.8 derives the graph-search rule and its large-depth limit, and Section 5 evaluates planning performance.

Optionally, local search can include a soft wall penalty,

C_{r}(\pi)=\sum_{j=1}^{\ell}\left[\frac{\big(a_{r}(x_{j})-d_{\mathrm{wall}}(x_{j})/h\big)_{+}}{r}\right]^{2},

where $(u)_{+}=\max(u,0)$ , $r$ is the preferred clearance in cells, and $d_{\mathrm{wall}}$ measures distance to obstacles or the outer boundary. Normally $a_{r}=r$ ; near a goal reachable within the local depth, it decreases with the remaining free-edge distance to that goal. Among paths retaining at least half the best positive local field gain, we prefer the smallest penalty. Narrow passages and starts or goals near walls are still possible.

We evaluate how backbone, supervision, training objective, and controller affect planning, then test transfer across grid resolutions and replanning under execution noise. Table 1 compares controllers on the same fields, while Table 5 compares complete planners under common execution rules. Models are pretrained on generated environments and evaluated on unseen maps without fine-tuning. The main 2D application suite contains 592 maps and 37,888 starts across robotics, terrain, marine, urban, and game environments (Appendix E). Success measures goal-reaching rate; relative $L_{2}$ measures field error against FEM references. Unless specified, local search (LS) denotes local search (10). Bold marks the best value in each comparison, including ties. Experiments use NVIDIA B300 and B200 GPUs; training and evaluation protocols are in Appendix D.1.

Original Table 1 — Table 1: 2D application-benchmark evaluation of neural backbones trained with WoSt supervision and direct WoSt inference: planning success and relative $L_{2}$ error against the FEM ground-truth field.

Table 1
Model	Gradient (%)	Local search (4) (%)	Local search (10) (%)	Relative $L_{2}$
UPT	2.48	2.78	3.57	0.5533
Transolver	3.26	5.34	15.30	0.5158
WoSt	2.71	6.30	36.45	0.1852
GINO	46.56	51.24	53.11	0.2628
U-Net	74.21	84.22	88.98	0.1605
DISCO	82.49	87.60	91.08	0.1041
LocalNO	80.34	86.63	91.35	0.0915
U-NO	78.08	86.11	91.68	0.0792
FNO	83.87	91.25	93.62	0.0615

Table 1 compares FNO (Li et al., 2021), U-NO (Rahman et al., 2023), U-Net $^{2}$ (Ronneberger et al., 2015), GINO (Li et al., 2023b), LocalNO (Liu-Schiaffini et al., 2024), DISCO (Ocampo et al., 2023), UPT (Alkin et al., 2024), and Transolver (Wu et al., 2024), all trained with WoSt supervision. FNO has the highest success under all three controllers and the lowest relative $L_{2}$ error (0.0615); its local-search (10) success is 93.62%, compared with 91.68% for the next-best U-NO. These rankings apply to the tested implementations under the training-compute policy in Appendix D.1. Direct WoSt inference reaches only 2.71% gradient and 36.45% local-search (10) success, versus 83.87% and 93.62% for FNO. Thus, at the reported solver settings, learning from approximate WoSt targets yields more useful planning fields than direct inference. We use FNO in the remaining ablations.

With FNO fixed, Table 4 compares WoSt targets, Deep Ritz training (E and Yu, 2018), and PDE-residual training with PINO (Li et al., 2024). WoSt achieves 93.62% local-search success, versus 91.46% for Deep Ritz and 59.69% for PINO; gradient success follows the same ordering (83.87%, 81.41%, and 32.50%). Deep Ritz has slightly lower field error than WoSt (0.0608 versus 0.0615), but lower planning success. Field error is therefore not a reliable indicator of how useful these learned fields are for control.

Original table — Table 2 FNO supervision on 2D application benchmarks: success and relative error.
Table 3 Walk-horizon ablation on the generated 2D test set: planning success (%).
Table 4 FNO loss ablations on 2D application benchmarks: planning success (%).

Paper table
Source	Grad. (%)	LS (%)	Rel. err.
WoSt	83.87	93.62	0.0615
Deep Ritz	81.41	91.46	0.0608
PINO	32.50	59.69	0.5644

Paper table
Steps	Grad. (%)	LS (%)
10	32.18	51.89
100	40.55	66.80
400	38.06	69.69
800	37.84	69.66

Paper table
Model	Grad. (%)	LS (%)
Full model	83.87	93.62
$-$ value loss	67.71	83.26
$-$ goal loss	82.93	93.11
$-$ gradient loss	80.17	92.04

In Table 4, we vary the walk cap using a smaller training budget of 16 maps, eight goals per map, four trajectories per query, and 40 epochs, evaluated on generated 2D maps. Local-search success increases from 51.89% at 10 steps to 69.69% at 400 steps, but is 69.66% at 800 steps. Gradient success instead peaks at 40.55% with 100 steps. Longer walks therefore do not uniformly improve planning: the horizon changes both sampled propagation and reliance on the learned continuation (Section 4.2). Further details of this study are in Appendix D.2. These success rates should not be compared directly with the full-budget application results.

In Table 4, we remove one loss term at a time from Eq. 13. Removing value supervision causes the largest decrease: gradient success falls from 83.87% to 67.71%, and local-search success from 93.62% to 83.26%. Removing goal or gradient supervision gives smaller local-search decreases, to 93.11% and 92.04%, respectively. All three terms help in this setting, with value supervision making the largest measured contribution. These ablations use direct desirability supervision.

On the same learned FNO fields, Table 1 shows success rising from 83.87% with gradient following to 91.25% with local search (4) and 93.62% with local search (10), a total gain of 9.75 percentage points. This test-time refinement requires no retraining. In 2D, the two local-search variants use $(L,e)=(5,1)$ and $(10,3)$ , so both depth and committed prefix length change. Local search (10) also exceeds gradient success for every other backbone, although it reaches only 36.45% on direct WoSt fields. Looking beyond immediate neighbors can therefore recover useful routes, but does not eliminate dependence on field quality. Figure 7 illustrates the two search neighborhoods on one fixed field; controller details are in Appendix D.1.

In Fig. 1(b), we test models trained at $256\times 256$ on grids up to $4096\times 4096$ , without retraining. FNO maintains 90.62–96.25% local-search (10) success across the five resolutions and reaches 92.50% at $4096\times 4096$ , compared with 71.25% for U-NO and 57.81% for U-Net. U-NO leads at $512\times 512$ (95.31% versus 94.06% for FNO); FNO leads at the other four resolutions. GINO runs out of memory above $256\times 256$ . The quantitative comparison supports FNO’s transfer across these discretizations; Fig. 4 separately illustrates predicted fields and paths at three resolutions.

Table 5 compares PUN with A⋆ (Hart et al., 1968), Dijkstra (Dijkstra, 1959), RRT (LaValle and Kuffner, 2001), MPPI (Williams et al., 2017a), and PI $^{2}$ -CMA (Stulp and Sigaud, 2012) on 2D and 3D application benchmarks and 4D manipulators. Obstacle hits count as failures (Appendix D.1). PUN has the highest reported success in all 12 dimension–noise settings, tying A⋆ and Dijkstra at $\lambda=0.25$ in 4D. At $\lambda=4$ , PUN exceeds the strongest baseline, RRT, by 9.00 percentage points in 2D (45.67% versus 36.67%), 4.00 in 3D (16.00% versus 12.00%), and 9.17 in 4D (95.00% versus 85.83%). PUN also has the lowest reported runtime throughout: 9.55–14.51 ms/start in 2D, 14.30–15.64 in 3D, and 10.51–11.66 in 4D. At $\lambda=4$ in 4D, its 11.66 ms/start is 3.29 times faster than the fastest baseline, RRT (38.31 ms/start); Fig. 1(a) summarizes mean 4D runtimes. The advantage is relative: PUN’s 3D success still falls from 92.00% to 16.00% as noise increases. These comparisons concern the reported implementations and settings. Figure 3 shows PUN reaching the goal in two illustrative 4D rollouts at $\lambda=0.25$ and $2$ .

Original Table 5 — Table 5: Replanning under execution noise on 2D and 3D application benchmarks and 4D manipulators. Each entry reports success (%) / runtime (ms/start). Bold marks the highest success and lowest runtime within each dimension and noise level, including ties.

Table 5
Dim.	Method	$\lambda=0.25$	$\lambda=1$	$\lambda=2$	$\lambda=4$
2D	A*	59.00 / 43.63	42.33 / 47.09	35.00 / 51.18	23.67 / 53.76
Dijkstra	59.00 / 154.94	42.33 / 174.95	35.00 / 203.41	23.67 / 239.97
RRT	82.00 / 98.45	68.67 / 103.87	52.33 / 112.27	36.67 / 124.71
MPPI	84.00 / 1734.11	61.00 / 1773.31	45.67 / 1838.73	31.67 / 2126.12
PI $^{2}$ -CMA	78.33 / 1330.76	60.00 / 1484.91	48.00 / 1524.70	35.00 / 1606.82
PUN (zero-shot)	91.67 / 9.55	74.00 / 11.14	62.33 / 14.07	45.67 / 14.51
3D	A*	36.67 / 73.40	23.33 / 60.47	12.00 / 73.85	6.00 / 65.25
Dijkstra	37.33 / 73.46	22.67 / 63.04	12.00 / 76.99	6.00 / 67.73
RRT	53.33 / 238.31	35.33 / 169.90	25.33 / 122.85	12.00 / 126.23
MPPI	62.67 / 532.31	34.67 / 577.84	18.67 / 554.11	7.33 / 704.31
PI $^{2}$ -CMA	47.33 / 452.25	22.00 / 435.01	21.33 / 475.98	8.67 / 606.55
PUN (zero-shot)	92.00 / 15.64	68.67 / 15.34	40.00 / 14.30	16.00 / 14.43
4D	A*	100.00 / 16.81	96.67 / 17.08	72.50 / 78.02	59.17 / 62.44
Dijkstra	100.00 / 17.17	96.67 / 17.34	72.50 / 66.75	59.17 / 59.39
RRT	98.33 / 24.20	95.83 / 25.80	90.00 / 43.00	85.83 / 38.31
MPPI	95.00 / 411.24	92.50 / 398.18	81.67 / 364.22	71.67 / 374.83
PI $^{2}$ -CMA	96.67 / 379.79	92.50 / 357.69	86.67 / 355.44	78.33 / 356.72
PUN	100.00 / 10.63	100.00 / 10.51	100.00 / 10.55	95.00 / 11.66

Several extensions are possible. A spatially varying cost $c(x)$ is the most useful: with constant $c$ , the state-cost term penalizes elapsed time and the objective also penalizes control effort. Reflection carries no additional contact cost; choosing $c(x)$ as a function of obstacle clearance would add a clearance-dependent running cost. Transient problems replace the stationary equation by its time-dependent counterpart, giving a finite-horizon value function rather than a first-exit one, which suits tasks with deadlines or moving goals. Partial information is the most direct use of our amortization: since the operator maps geometry to a field in one forward pass, replanning as a robot discovers obstacles requires only a forward pass on the updated geometry instead of a new solve. The open question is how to condition the operator on a partially observed occupancy map.

The value field is derived for reflecting boundaries. Local search allows its use in collision-failure tasks, but nominal graph feasibility does not guarantee collision-free execution under noise. Deriving the exact value function for that execution model would require a corresponding boundary condition. Assumption 1 is a further restriction. The control–noise matching assumption enables exact HJB linearization and is standard in classical path-integral control (Kappen, 2005; Theodorou et al., 2010). Our experiments are also limited to settings where this assumption is reasonable: we test geometric planning problems in up to four dimensions and do not include full robot models with contact. Evaluating on humanoid robots with contact would test both the assumption and the capacity of the operator. Finally, adversarial geometries such as mazes remain challenging, since long corridors can lead to weak or poorly resolved learned gradients.

We introduced Planning Under Noise (PUN), a neural operator framework that learns noise-aware planning fields from approximate Walk-on-Stars supervision. Dimension-general screening kernels support mesh-free supervision beyond three dimensions, and local search converts predicted fields into paths on unseen geometries. Experiments across 2D–4D environments demonstrate zero-shot planning and effective learning from stochastic labels. In replanning comparisons under execution noise, PUN achieves the lowest reported runtime among the tested methods, while matching or exceeding the success of every compared planner at each tested noise level across 2D–4D.

In this work, we used generative AI tools to assist with parts of method implementation. Additionally, we used generative AI tools to polish writing and proofs and identify relevant literature. We did not use generative AI tools for other applicable tasks requiring disclosure. All AI-assisted work was reviewed and verified by the authors. We take responsibility for the final content of this work, including text, claims, and artifacts produced with the aid of generative AI.

Anima Anandkumar and Xi Deng are both supported in part by the Bren endowed chair and the AI2050 senior fellow program at Schmidt Sciences. Hong Chul Nam is supported by a PhD fellowship for the Center for Data Science at NYU.

References

Alkin et al. (2024) B. Alkin, A. Fürst, S. Schmid, L. Gruber, M. Holzleitner, and J. Brandstetter Universal physics transformers: a framework for efficiently scaling neural operators. In Advances in Neural Information Processing Systems (NeurIPS), Cited by: §3, §5.
Bellman (1957) R. Bellman Dynamic programming. Princeton University Press, Princeton, NJ. Cited by: §1.
Cole (1951) J. D. Cole On a quasi-linear parabolic equation occurring in aerodynamics. Quarterly of Applied Mathematics 9 (3), pp. 225–236. Cited by: §2.1.
Dijkstra (1959) E. W. Dijkstra A note on two problems in connexion with graphs. Numerische Mathematik 1 (1), pp. 269–271. Cited by: §3, §5.
E and Yu (2018) W. E and B. Yu The Deep Ritz method: a deep learning-based numerical algorithm for solving variational problems. Communications in Mathematics and Statistics 6 (1), pp. 1–12. External Links: Document Cited by: §5.Source
Fleming and Soner (2006) W. H. Fleming and H. M. Soner Controlled markov processes and viscosity solutions. 2 edition, Stochastic Modelling and Applied Probability, Vol. 25, Springer, New York. Cited by: §1.
Fleming (1977) W. H. Fleming Exit probabilities and optimal stochastic control. Applied Mathematics and Optimization 4 (1), pp. 329–346. Cited by: §2.1.
Gilbarg and Trudinger (2001) D. Gilbarg and N. S. Trudinger Elliptic partial differential equations of second order. 2nd edition, Springer. Cited by: §C.1, §C.3.
Han et al. (2018) J. Han, A. Jentzen, and W. E Solving high-dimensional partial differential equations using deep learning. Proceedings of the National Academy of Sciences 115 (34), pp. 8505–8510. Cited by: §1.
Hart et al. (1968) P. E. Hart, N. J. Nilsson, and B. Raphael A formal basis for the heuristic determination of minimum cost paths. IEEE Transactions on Systems Science and Cybernetics 4 (2), pp. 100–107. Cited by: §1, §3, §5.
Hopf (1950) E. Hopf The partial differential equation $u_{t}+uu_{x}=\mu u_{xx}$. Communications on Pure and Applied Mathematics 3 (3), pp. 201–230. External Links: Document Cited by: §2.1.Source
Jambon et al. (2026) C. Jambon, M. S. Nabizadeh, and M. Konaković Luković Walk on decomposed subdomains: a hybrid monte carlo–deterministic solver for elliptic pdes. ACM Transactions on Graphics (TOG) 45 (4), pp. 132:1–132:22. Cited by: §1.
Kappen (2005) H. J. Kappen Path integrals and symmetry breaking for optimal control theory. Journal of Statistical Mechanics: Theory and Experiment 2005 (11), pp. P11011. Cited by: §2.1, §6.
Karaman and Frazzoli (2011) S. Karaman and E. Frazzoli Sampling-based algorithms for optimal motion planning. The International Journal of Robotics Research 30 (7), pp. 846–894. Cited by: §3.
Kavraki et al. (1996) L. E. Kavraki, P. Svestka, J. Latombe, and M. H. Overmars Probabilistic roadmaps for path planning in high-dimensional configuration spaces. IEEE Transactions on Robotics and Automation 12 (4), pp. 566–580. Cited by: §1, §3.
Kovachki et al. (2023) N. Kovachki, Z. Li, B. Liu, K. Azizzadenesheli, K. Bhattacharya, A. Stuart, and A. Anandkumar Neural operator: learning maps between function spaces with applications to pdes. Journal of Machine Learning Research 24 (89), pp. 1–97. Cited by: §1, §3.
Latombe (1991) J. Latombe Robot motion planning. The Kluwer International Series in Engineering and Computer Science, Vol. 124, Kluwer Academic Publishers. Cited by: §1.
LaValle and Kuffner (2001) S. M. LaValle and J. J. Kuffner Randomized kinodynamic planning. The International Journal of Robotics Research 20 (5), pp. 378–400. Cited by: §1, §3, §5.
LaValle (2006) S. M. LaValle Planning algorithms. Cambridge University Press, Cambridge, UK. Cited by: §1.
Li et al. (2023a) Z. Li, G. Yang, X. Deng, C. De Sa, B. Hariharan, and S. Marschner Neural caches for monte carlo partial differential equation solvers. In SIGGRAPH Asia 2023 Conference Papers, pp. 1–10. Cited by: §3.
Li et al. (2021) Z. Li, N. Kovachki, K. Azizzadenesheli, B. Liu, K. Bhattacharya, A. Stuart, and A. Anandkumar Fourier neural operator for parametric partial differential equations. In International Conference on Learning Representations (ICLR), Cited by: §2.2, §3, §5.
Li et al. (2023b) Z. Li, N. Kovachki, C. Choy, B. Li, J. Kossaifi, S. Otta, M. A. Nabian, M. Stadler, C. Hundt, K. Azizzadenesheli, and A. Anandkumar Geometry-informed neural operator for large-scale 3d pdes. In Advances in Neural Information Processing Systems (NeurIPS), Cited by: §3, §5.
Li et al. (2024) Z. Li, H. Zheng, N. Kovachki, D. Jin, H. Chen, B. Liu, K. Azizzadenesheli, and A. Anandkumar Physics-informed neural operator for learning partial differential equations. ACM / IMS Journal of Data Science 1 (3), pp. 1–27. Cited by: §3, §5.
Liu-Schiaffini et al. (2024) M. Liu-Schiaffini, J. Berner, B. Bonev, T. Kurth, K. Azizzadenesheli, and A. Anandkumar Neural operators with localized integral and differential kernels. In International Conference on Machine Learning (ICML), Cited by: §3, §5.
Matada et al. (2025) S. Matada, L. Bhan, Y. Shi, and N. Atanasov Generalizable motion planning via operator learning. In International Conference on Learning Representations (ICLR), Cited by: §3.
Muller (1956) M. E. Muller Some continuous monte carlo methods for the dirichlet problem. The Annals of Mathematical Statistics 27 (3), pp. 569–589. Cited by: §1, §1.
Nam et al. (2024) H. C. Nam, J. Berner, and A. Anandkumar Solving poisson equations using neural walk-on-spheres. In Proceedings of the 41st International Conference on Machine Learning, Proceedings of Machine Learning Research, Vol. 235. Cited by: §3.
Ni and Qureshi (2023) R. Ni and A. H. Qureshi NTFields: neural time fields for physics-informed robot motion planning. In International Conference on Learning Representations (ICLR), Cited by: §3.
Nüsken and Richter (2021) N. Nüsken and L. Richter Interpolating between BSDEs and PINNs: deep learning for elliptic and parabolic boundary value problems. Note: arXiv preprint arXiv:2112.03749 Cited by: §3.
Ocampo et al. (2023) J. Ocampo, M. A. Price, and J. D. McEwen Scalable and equivariant spherical cnns by discrete-continuous (disco) convolutions. In International Conference on Learning Representations (ICLR), Cited by: §3, §5.
Rahman et al. (2023) M. A. Rahman, Z. E. Ross, and K. Azizzadenesheli U-no: u-shaped neural operators. Transactions on Machine Learning Research. Cited by: §3, §5.
Ronneberger et al. (2015) O. Ronneberger, P. Fischer, and T. Brox U-net: convolutional networks for biomedical image segmentation. In Medical Image Computing and Computer-Assisted Intervention (MICCAI), pp. 234–241. Cited by: §3, §5.
Sawhney and Crane (2020) R. Sawhney and K. Crane Monte carlo geometry processing: a grid-free approach to pde-based methods on volumetric domains. ACM Transactions on Graphics 39 (4). Cited by: §1, §3.
Sawhney et al. (2023) R. Sawhney, B. Miller, I. Gkioulekas, and K. Crane Walk on stars: a grid-free monte carlo method for pdes with neumann boundary conditions. ACM Transactions on Graphics 42 (4). Cited by: §B.4, §B.4, §B.5, §C.4, item 1, §D.1, §1, §1, §4.2.
Sethian (1996) J. A. Sethian A fast marching level set method for monotonically advancing fronts. Proceedings of the National Academy of Sciences 93 (4), pp. 1591–1595. Cited by: §3.
Smith et al. (2021) J. D. Smith, K. Azizzadenesheli, and Z. E. Ross EikoNet: solving the eikonal equation with deep neural networks. IEEE Transactions on Geoscience and Remote Sensing 59 (12), pp. 10685–10696. Cited by: §3.
Stulp and Sigaud (2012) F. Stulp and O. Sigaud Path integral policy improvement with covariance matrix adaptation. In International Conference on Machine Learning, Cited by: §D.1, §1, §3, §5.
Sugimoto et al. (2023) R. Sugimoto, T. Chen, Y. Jiang, C. Batty, and T. Hachisuka A practical walk-on-boundary method for boundary value problems. ACM Transactions on Graphics 42 (4). Cited by: §3.
Theodorou et al. (2010) E. Theodorou, J. Buchli, and S. Schaal A generalized path integral control approach to reinforcement learning. The Journal of Machine Learning Research 11, pp. 3137–3181. Cited by: §1, §2.1, §6.
Thrun (2002) S. Thrun Probabilistic robotics. Communications of the ACM 45 (3), pp. 52–57. Cited by: §1.
Todorov (2006) E. Todorov Linearly-solvable Markov decision problems. In Advances in Neural Information Processing Systems (NIPS), Vol. 19, pp. 1369–1376. Cited by: §2.1.
Todorov (2009) E. Todorov Efficient computation of optimal actions. Proceedings of the National Academy of Sciences 106 (28), pp. 11478–11483. Cited by: §2.1.
Viswanath et al. (2026) H. Viswanath, H. C. Nam, X. Deng, J. Berner, A. Anandkumar, and A. Bera Operator learning using weak supervision from walk-on-spheres. arXiv preprint arXiv:2603.01193. Cited by: §1.
Williams et al. (2017a) G. Williams, A. Aldrich, and E. A. Theodorou Model predictive path integral control: from theory to parallel computation. Journal of Guidance, Control, and Dynamics 40 (2), pp. 344–357. Cited by: §D.1, §1, §3, §5.
Williams et al. (2017b) G. Williams, N. Wagener, B. Goldfain, P. Drews, J. M. Rehg, B. Boots, and E. A. Theodorou Information theoretic MPC for model-based reinforcement learning. In IEEE International Conference on Robotics and Automation (ICRA), pp. 1714–1721. Cited by: §3.
Wu et al. (2024) H. Wu, H. Luo, H. Wang, J. Wang, and M. Long Transolver: a fast transformer solver for pdes on general geometries. In International Conference on Machine Learning (ICML), Cited by: §3, §5.
Zhang (2017) J. Zhang Backward stochastic differential equations. In Backward Stochastic Differential Equations: From Linear to Fully Nonlinear Theory, pp. 79–99. Cited by: §3.

Table 6 collects the shared symbols used across the paper.

Original Table 6 — Table 6: Notation, with the equation, assumption, or section where each symbol is defined.

Table 6
Symbol	Definition	Meaning
Control problem
$\xi,\ x$	§2.1	starting point of the process or walk; conditioned as $\mid\xi=x$ where the start is not otherwise fixed
$X^{\xi,\mathbf{u}}_{t}$	Eq. 1	controlled process at time $t$ , from the prescribed start $\xi$ under control $\mathbf{u}$
$\mathbf{b}$	Eq. 1	drift of the controlled system ( $\mathbf{b}=0$ in the specialization)
$G,\ R$	Eq. 1, Eq. 2	control matrix and control penalty ( $G=R=I$ in the specialization)
$\Sigma$	Eq. 1	diffusion tensor; $\Sigma=\lambda I$ when $G=R=I$ under matching
$q,\ c$	Eq. 2, Eq. 3	state cost $q=\tfrac{1}{2}c^{2}$ , with $c>0$ the local slowness
$\phi,\ g$	Eq. 2, §2.1	terminal cost and its transform $g=e^{-\phi/\lambda}\in(0,1]$ ; $g\equiv 1$ by default
$\lambda$	Ass. 1	temperature; sets the noise–control exchange rate
$V,\ \Psi$	Eq. 2, §2.1	value function and desirability $\Psi=e^{-V/\lambda}$
Geometry and the elliptic problem
$\mathcal{G},\ \Gamma_{D}$	Eq. 4	open goal region and its absorbing boundary $\Gamma_{D}=\partial\mathcal{G}$
$\Gamma_{N}$	Eq. 4	obstacles $\partial\mathcal{S}$ and outer boundary, reflecting
$\kappa$	Eq. 3	screening coefficient $\kappa=c/\lambda$
$\omega_{d}$	Eq. 15	total solid angle, $2\pi^{d/2}/\Gamma(d/2)$
$\alpha(x)$	Eq. 8	fraction of solid angle subtended by $\Omega$ at $x$ ; $1$ inside, $\tfrac{1}{2}$ on $\Gamma_{N}$
The walk
$\xi_{k}$	Eq. 7	WoSt point after $k$ steps, from $\xi_{0}=\xi$ ; random for $k\geq 1$
$r(x),\ \mathrm{St}(x)$	Eq. 6	radius and star centered at the current spatial location $x$ ; $r_{k}=r(\xi_{k})$
$d_{D},\ d_{\mathrm{sil}}$	Eq. 6	distance to $\Gamma_{D}$ ; distance to the nearest visibility silhouette
$\rho_{k}$	Eq. 7	distance travelled by step $k$ ; $\rho_{k}\leq r_{k}$
$\sigma$	§4.1	first exact goal-hit step $\inf\{k:\xi_{k}\in\Gamma_{D}\}$ , possibly infinite; practical walks use shell stopping and a horizon
$W_{k}$	Eq. 25	step weight $Q^{\kappa}_{d}(\rho_{k},r_{k})$
$\gamma_{d}$	Eq. 22	full-sphere weight $Q^{\kappa}_{d}(r,r)$
Labels and early stopping
$\widehat{\Psi}$	Eq. 10	Monte Carlo estimate
$K,\ N,\ m$	§4.2, Alg. 2	step budget; walks per query point; query points per gradient step
Planning lattice
$h,\ \mathcal{X}_{h}$	§4.4	lattice spacing and free grid cells
$L,\ e$	§4.4	depth and number of edges committed per replan
$\mathcal{N},\ \mathcal{R}_{k}$	§4.4	feasible neighbors of a state; states reachable in at most $k$ moves
$B_{k}$	Eq. 37	largest desirability reachable within $k$ moves
Operator and training
$\mathcal{T}$	Eq. 5	exact solution operator $a\mapsto\Psi$ , which $\mathcal{T}_{\theta}$ approximates
$\mathcal{T}_{\theta},\ \theta^{-}$	Eq. 11, Alg. 2	the learned operator $a\mapsto\Psi_{\theta}$ ; its EMA target parameters
$\mathfrak{L}$	Eq. 13	training objective
Kernels
$\nu$	Eq. 25	Bessel order $\tfrac{d}{2}-1$
$G^{\kappa,B}$	Eq. 18	screened Green’s function of the ball
$P^{\mathbb{R}^{d}}$	Eq. 24	unscreened free-space Poisson kernel, the sampling density

This appendix derives the kernels used by the estimator of Section 4.1 in general dimension. Throughout, $d\geq 2$ , $\kappa>0$ , and

\nu:=\frac{d}{2}-1,\qquad r:=\|x-\mathbf{y}\|,\qquad\omega_{d}=\frac{2\pi^{d/2}}{\Gamma(d/2)}.

In Appendix B.1 the radial variable is the source separation $r=\|x-\mathbf{y}\|$ ; from Appendix B.2 onward $r$ is the radius of the ball or star and $\rho:=\|x-\mathbf{y}\|$ denotes the source separation. We seek Green’s functions for the screened Laplacian

\mathcal{L}_{\kappa}:=-\Delta+\kappa^{2},\qquad\mathcal{L}_{\kappa}\,G(x,\mathbf{y})=\delta(x-\mathbf{y}).

Away from the pole every case below solves the same radial ODE,

u^{\prime\prime}(r)+\frac{d-1}{r}\,u^{\prime}(r)-\kappa^{2}u(r)=0,\qquad\text{basis }\ \big\{\,r^{-\nu}I_{\nu}(\kappa r),\ \ r^{-\nu}K_{\nu}(\kappa r)\,\big\},

with $I_{\nu}$ regular at the origin and $K_{\nu}$ singular. Cases differ only in which combination the outer boundary condition selects; the overall constant is fixed once and for all by the flux jump at $r=0$ .

In free space the radial solution of Eq. 16 that decays at infinity is the $K_{\nu}$ branch, and the flux jump at the pole fixes its constant:

G_{\mathbb{R}^{d}}(r)=\frac{\kappa^{\nu}}{(2\pi)^{d/2}}\,\frac{K_{\nu}(\kappa r)}{r^{\nu}}

As a check, for $d=2$ ( $\nu=0$ ) this gives $G=K_{0}(\kappa r)/2\pi$ , and for $d=3$ ( $\nu=\tfrac{1}{2}$ ), using $K_{1/2}(z)=\sqrt{\pi/2z}\,e^{-z}$ gives the Yukawa potential $G=e^{-\kappa r}/4\pi r$ .

The Neumann part of $\partial\mathrm{St}$ requires the normal derivative of Eq. 17. Using $\tfrac{d}{dz}\big[z^{-\nu}K_{\nu}(z)\big]=-z^{-\nu}K_{\nu+1}(z)$ ,

\partial_{r}G_{\mathbb{R}^{d}}(r)=-\frac{\kappa^{\nu+1}}{(2\pi)^{d/2}}\,\frac{K_{\nu+1}(\kappa r)}{r^{\nu}}\;<\;0

which is strictly negative for all $r>0$ , so the kernel is monotonically decreasing in every dimension.

Let $B_{R}=\{\|x\|<R\}$ with the pole at the center, as in every walk step. The solution is radial and takes the form

G_{B_{R}}(r)=\frac{\kappa^{\nu}}{(2\pi)^{d/2}}\,\frac{1}{r^{\nu}}\Big[K_{\nu}(\kappa r)-A\,I_{\nu}(\kappa r)\Big],\qquad 0<r<R,

with the same prefactor as Eq. 17, since the flux jump at the origin is unchanged by the boundary. Only $A$ depends on the boundary condition. Writing $z=\kappa r$ and using the scaled-derivative identities

\frac{\mathrm{d}}{\mathrm{d}z}\big[z^{-\nu}I_{\nu}(z)\big]=z^{-\nu}I_{\nu+1}(z),\qquad\frac{\mathrm{d}}{\mathrm{d}z}\big[z^{-\nu}K_{\nu}(z)\big]=-z^{-\nu}K_{\nu+1}(z),

the radial derivative of Eq. 18 is proportional to $-K_{\nu+1}(z)-A\,I_{\nu+1}(z)$ . The two cases we use follow:

\begin{aligned}\displaystyle\text{Dirichlet }G(R)=0: & & \displaystyle A_{D} & \displaystyle=\frac{K_{\nu}(\kappa R)}{I_{\nu}(\kappa R)},\end{aligned}

\begin{aligned}\displaystyle\text{Neumann }\partial_{r}G(R)=0: & & \displaystyle A_{N} & \displaystyle=-\frac{K_{\nu+1}(\kappa R)}{I_{\nu+1}(\kappa R)}.\end{aligned}

Note the order shift in $A_{N}$ : since $\tfrac{\mathrm{d}}{\mathrm{d}z}[z^{-\nu}K_{\nu}]=-z^{-\nu}K_{\nu+1}$ and likewise for $I_{\nu}$ , imposing $\partial_{r}G(R)=0$ on $r^{-\nu}[K_{\nu}-AI_{\nu}]$ produces the order- $(\nu+1)$ ratio rather than the naive $K^{\prime}_{\nu}/I^{\prime}_{\nu}$ , for which the derivative vanishes only when $\nu=0$ .

The walk uses the Poisson kernel, $P^{\kappa}_{B_{R}}(\mathbf{0},\mathbf{z})=-\partial_{n_{\mathbf{z}}}G^{D}_{B_{R}}$ for $\mathbf{z}\in\partial B_{R}$ , which by radial symmetry is constant on the sphere. Its total mass has a closed form, and it is exactly the factor $\gamma_{d}$ appearing in Eq. 22.

Proposition 1 (Total mass of the screened Poisson kernel).

Let $u$ solve $\Delta u-\kappa^{2}u=0$ in $B_{R}$ with $u\equiv 1$ on $\partial B_{R}$ . Then

u(r)=\Big(\frac{R}{r}\Big)^{\nu}\,\frac{I_{\nu}(\kappa r)}{I_{\nu}(\kappa R)},\qquad\int_{\partial B_{R}}P^{\kappa}_{B_{R}}(\mathbf{0},\mathbf{z})\,\mathrm{d}S(\mathbf{z})\;=\;u(\mathbf{0})\;=\;\gamma_{d}(\kappa R),

where the full-sphere weight is

\gamma_{d}(z):=\frac{(z/2)^{\nu}}{\Gamma(\nu+1)\,I_{\nu}(z)}=Q^{\kappa}_{d}(r,r)\big|_{z=\kappa r}.

Consequently

P^{\kappa}_{B_{R}}(\mathbf{0},\mathbf{z})=\frac{\gamma_{d}(\kappa R)}{\omega_{d}\,R^{d-1}}=\frac{\Gamma(d/2)\,\gamma_{d}(\kappa R)}{2\pi^{d/2}R^{d-1}}.

Proof.

Regularity at the origin selects the $I_{\nu}$ branch of Eq. 16, and the stated normalization gives $u(R)=1$ . Taking $r\to 0$ with $I_{\nu}(z)\sim(z/2)^{\nu}/\Gamma(\nu+1)$ ,

u(\mathbf{0})=\frac{R^{\nu}(\kappa/2)^{\nu}}{\Gamma(\nu+1)\,I_{\nu}(\kappa R)}=\frac{(\kappa R/2)^{\nu}}{\Gamma(\nu+1)\,I_{\nu}(\kappa R)}=\gamma_{d}(\kappa R).

The mass identity follows since $u(\mathbf{0})=\int P^{\kappa}\,u\,\mathrm{d}S$ with $u\equiv 1$ on the boundary, and Eq. 23 follows from constancy on the sphere. ∎

Proposition 1 is the $\rho=r$ case of the star-step weight Eq. 25, and Appendix C.5 shows the general weight reduces to it. It also gives the two asymptotic regimes. Expanding $I_{\nu}$ and using $4(\nu+1)=2d$ ,

\gamma_{d}(z)=\Big[1+\frac{z^{2}}{2d}+O(z^{4})\Big]^{-1}=1-\frac{z^{2}}{2d}+O(z^{4}),\qquad\gamma_{d}(z)\sim\frac{(z/2)^{\nu}\sqrt{2\pi z}}{\Gamma(\nu+1)}\,e^{-z}\ \ (z\to\infty),

so a short step is nearly weight-preserving with deficit $\kappa^{2}R^{2}/2d$ , while a step much longer than the screening length is damped exponentially; hence steps beyond a few screening lengths contribute negligibly.

Appendix B.3 treated the Poisson kernel at the center of a ball, which corresponds to a full-length sphere step. A star step may instead terminate against $\Gamma_{N}$ at a distance $\rho<r$ , so we need the kernel at interior points of the ball.

Following Sawhney et al. (2023), the screened ball Green’s functions are

\begin{aligned}\displaystyle G^{\kappa,B}_{2}(x,\mathbf{y}) & \displaystyle=\frac{1}{2\pi}\Big[K_{0}(\kappa\rho)-I_{0}(\kappa\rho)\tfrac{K_{0}(\kappa r)}{I_{0}(\kappa r)}\Big], & \displaystyle G^{\kappa,B}_{3}(x,\mathbf{y}) & \displaystyle=\frac{1}{4\pi}\,\frac{\sinh\big((r-\rho)\kappa\big)}{\rho\,\sinh(\kappa r)},\end{aligned}

with $\rho=\|\mathbf{y}-x\|$ , which agree with the general- $d$ expression Eq. 18 at $d=2,3$ . Both vanish at $\rho=r$ , which is the cancellation that removes the unknown-flux term on the spherical part of the boundary.

The corresponding Poisson kernels factor as a scalar times the free-space kernel,

P^{\kappa,B}_{d}(x,\mathbf{z})=Q^{\kappa}_{d}(\rho,r)\,P^{\mathbb{R}^{d}}(x,\mathbf{z}),\qquad P^{\mathbb{R}^{d}}(x,\mathbf{z})=\frac{n_{\mathbf{z}}\cdot(\mathbf{z}-x)}{\omega_{d}\,\rho^{d}},

with $Q^{\kappa}_{d}$ given for all $d$ by Eq. 25.

Lemma 1 (General-dimension step weight).

For every $d\geq 2$ and $0<\rho\leq r$ ,

Q^{\kappa}_{d}(\rho,r):=\frac{\partial_{\rho}G^{\kappa,B}}{\partial_{\rho}G^{B}}=\frac{2}{\Gamma(\nu+1)}\Big(\frac{\kappa\rho}{2}\Big)^{\nu+1}\left[K_{\nu+1}(\kappa\rho)+\frac{K_{\nu}(\kappa r)}{I_{\nu}(\kappa r)}I_{\nu+1}(\kappa\rho)\right],\qquad\nu=\tfrac{d}{2}-1,

and $Q^{\kappa}_{d}(r,r)=\gamma_{d}(\kappa r)$ .

Proof.

Write $z=\kappa\rho$ and $\alpha=K_{\nu}(\kappa r)/I_{\nu}(\kappa r)$ , so that by Eq. 18 and Eq. 20,

G^{\kappa,B}=\kappa^{2\nu}(2\pi)^{-d/2}\big[z^{-\nu}K_{\nu}(z)-\alpha\,z^{-\nu}I_{\nu}(z)\big].

The scaled-derivative identities Eq. 19 give

\partial_{\rho}G^{\kappa,B}=-\frac{\kappa^{\nu+1}}{(2\pi)^{d/2}}\,\rho^{-\nu}\big[K_{\nu+1}(\kappa\rho)+\alpha\,I_{\nu+1}(\kappa\rho)\big],

using $2\nu+1-\nu=\nu+1$ for the power of $\kappa$ . For $d>2$ , the denominator has $G^{B}=c_{d}(\rho^{2-d}-r^{2-d})$ with $c_{d}=\Gamma(\nu)/(4\pi^{\nu+1})$ , hence $\partial_{\rho}G^{B}=-2\nu\,c_{d}\,\rho^{-2\nu-1}=-\Gamma(\nu+1)(2\pi^{\nu+1})^{-1}\rho^{-2\nu-1}$ , where $2\nu\Gamma(\nu)=2\Gamma(\nu+1)$ . Dividing the preceding display by this and collecting the powers of $2$ and $\pi$ gives Eq. 25. Setting $\rho=r$ makes the bracket $\big[I_{\nu}K_{\nu+1}+I_{\nu+1}K_{\nu}\big]/I_{\nu}$ , and the Bessel identity $I_{\nu}(z)K_{\nu+1}(z)+I_{\nu+1}(z)K_{\nu}(z)=1/z$ reduces it to $1/(zI_{\nu}(z))$ , leaving $(z/2)^{\nu}/(\Gamma(\nu+1)I_{\nu}(z))=\gamma_{d}(z)$ . For $d=2$ , use $G^{B}(\rho)=(2\pi)^{-1}\log(r/\rho)$ and $\partial_{\rho}G^{B}=-1/(2\pi\rho)$ directly; the same ratio follows. The sphere reduction follows from the standard Wronskian $I_{0}K_{1}+I_{1}K_{0}=1/z$ . We verified Eq. 25 against the tabulated $d=2,3$ forms and the reduction against $\gamma_{d}$ at $d\in\{2,3,4,5,7,10\}$ to fourteen digits. ∎

Specializing Eq. 25 to the two dimensions tabulated by Sawhney et al. (2023),

\begin{aligned}\displaystyle Q^{\kappa}_{2}(\rho,r) & \displaystyle=\Big[K_{1}(\kappa\rho)+I_{1}(\kappa\rho)\,\tfrac{K_{0}(\kappa r)}{I_{0}(\kappa r)}\Big]\,\kappa\rho,\end{aligned}

\begin{aligned}\displaystyle Q^{\kappa}_{3}(\rho,r) & \displaystyle=e^{-\kappa\rho}\big(\kappa\rho+1\big)+\big(\kappa\rho\cosh(\kappa\rho)-\sinh(\kappa\rho)\big)\,\tfrac{e^{-\kappa r}}{\sinh(\kappa r)}.\end{aligned}

Table 7 gives representative values. Because of this factorization, the estimator is inexpensive: the geometric part $P^{\mathbb{R}^{d}}$ is exactly the signed solid angle subtended by $\partial\mathrm{St}$ at $\mathbf{z}$ , so sampling directions uniformly gives unit importance weights, and screening enters only through the scalar $Q^{\kappa}_{d}$ .

Original Table 7 — Table 7: Star-step weight $Q^{\kappa}_{3}(\rho,r)$ at $r=1$ . The weight is $1$ for a zero-length step and decreases with the realized hit distance $\rho$ , reaching $\gamma_{3}(\kappa r)=\kappa r/\sinh(\kappa r)$ at $\rho=r$ . A step cut short by an obstacle loses less weight than a full sphere step.

Table 7
$\rho/r$	$0$	$1/6$	$1/3$	$1/2$	$2/3$	$5/6$	$1$
$\kappa r=1$	1.0000	0.9880	0.9593	0.9232	0.8880	0.8614	0.8509
$\kappa r=4$	1.0000	0.8558	0.6157	0.4086	0.2629	0.1766	0.1466

The kernels of Appendices B.1–B.4 are constructed for free space and for a ball, yet Eq. 28 integrates them over $\partial\mathrm{St}$ , which is not a sphere. The following identity justifies this.

Proposition 2 (Admissibility).

Let $A\subseteq\Omega$ be a bounded Lipschitz domain with $x\in A$ , let $\Psi$ solve $\mathcal{L}_{\kappa}\Psi=0$ on $A$ , and let $G$ be any function with $(-\Delta+\kappa^{2})G(x,\cdot)=\delta_{x}$ on a neighborhood of $\overline{A}$ , so that $x$ is its only singularity in $\overline{A}$ . Then with $P:=-\partial_{n}G$ and $n$ the outward normal,

\Psi(x)=\int_{\partial A}\Big[P(x,\mathbf{z})\,\Psi(\mathbf{z})+G(x,\mathbf{z})\,\partial_{n}\Psi(\mathbf{z})\Big]\mathrm{d}S(\mathbf{z}).

Proof.

Green’s second identity on $A$ for the pair $(\Psi,G)$ with the operator $\mathcal{L}_{\kappa}$ gives

\int_{A}\big(\Psi\,\mathcal{L}_{\kappa}G-G\,\mathcal{L}_{\kappa}\Psi\big)=\int_{\partial A}\big(G\,\partial_{n}\Psi-\Psi\,\partial_{n}G\big)

the zeroth-order terms $\kappa^{2}\Psi G$ canceling identically between the two products. The second integrand on the left vanishes by hypothesis, and the first collapses to $\Psi(x)$ by the defining property of $G$ . Rearranging gives Eq. 28. ∎

The chain’s one-step identity is the case $G=G^{\kappa,B}$ , the screened Dirichlet Green’s function of the ball $B=B(x,r)$ , applied on an admissible star $A=\mathrm{St}(x)$ . Since $G^{\kappa,B}$ vanishes on the spherical part of $\partial A$ and $\partial_{n}\Psi=0$ on the reflecting part, the flux term drops and $\Psi(x)=\int_{\partial A}P^{\kappa,B}(x,\mathbf{z})\Psi(\mathbf{z})\,\mathrm{d}S(\mathbf{z})$ with $P^{\kappa,B}=-\partial_{n}G^{\kappa,B}$ . At a smooth boundary pole $x\in\Gamma_{N}$ , the same steps on $A$ minus a half-ball of radius $\eta$ about $x$ , with $\eta\to 0$ , give the identity with $\Psi(x)$ replaced by $\tfrac{1}{2}\Psi(x)$ , because the small hemisphere carries half of the singular flux; this is the factor $\alpha$ in Eq. 8. The direction law is the solid-angle density $P^{B}\,\mathrm{d}S/\alpha$ with $P^{B}=-\partial_{n}G^{B}$ , and the weight is the flux ratio $Q_{d}^{\kappa}=P^{\kappa,B}/P^{B}$ of Eq. 34, so $\mathbb{E}[W_{k}\Psi(\xi_{k+1})\mid\xi_{k}=x]=\alpha(x)^{-1}\int_{\partial A}P^{\kappa,B}\Psi\,\mathrm{d}S=\Psi(x)$ , the one-step property of Section 4.1. Multiplying by the earlier weights gives the martingale property before $\sigma$ ; the stopped sequence is constant afterwards.

The hypotheses require neither that $G$ be the Green’s function of $A$ nor that it satisfy any boundary condition on $\partial A$ . All that is needed is a fundamental solution whose singularity inside $A$ is the single pole at $x$ . Since $\mathrm{St}(x,r)\subseteq B(x,r)$ , the ball kernel $G^{\kappa,B}$ has that property, and so do the free-space kernel and the ball kernel of any radius $r^{\prime}\geq r$ . Each of these kernels gives a valid representation of $\Psi(x)$ , so we must choose among them.

The choice is governed by the unknown flux $\partial_{n}\Psi$ in Eq. 28. On $\partial\mathrm{St}$ that flux is zero on $\partial\mathrm{St}_{N}\subseteq\Gamma_{N}$ but unknown on $\partial\mathrm{St}_{B}\subseteq\partial B(x,r)$ , which lies in the interior of $\Omega$ where no boundary condition applies. A usable identity therefore needs a kernel that vanishes on $\partial\mathrm{St}_{B}$ while remaining finite on $\partial\mathrm{St}_{N}$ , and that is exactly the ball’s Dirichlet kernel: $G^{\kappa,B}(x,\cdot)\equiv 0$ on $\partial B(x,r)\supseteq\partial\mathrm{St}_{B}$ , and it is finite and nonzero at the interior points constituting $\partial\mathrm{St}_{N}$ . The free-space kernel is admissible but unsuitable here: $G^{\kappa,\mathbb{R}^{d}}>0$ on $\partial\mathrm{St}_{B}$ , so the unknown flux term does not vanish.

The walk nonetheless samples directions by the free-space solid angle. For the Laplace operator the two Poisson kernels agree on all of $\partial\mathrm{St}$ , because the boundary corrector $G^{B}-G^{\mathbb{R}^{d}}$ for a centered pole is the constant radial homogeneous solution and so has zero normal derivative; directional sampling then importance-samples the kernel exactly and the step weight is $1$ . Screening replaces that constant by the $I_{\nu}$ branch of Eq. 16, and the resulting $\rho$ -dependent corrector is precisely the second term of Eq. 25. This explains why $Q^{\kappa}_{d}$ depends on the star radius $r$ and not only on the hit distance $\rho$ : the $\rho$ -dependence comes from the free-space fundamental solution, and the $r$ -dependence comes from the ball corrector, whose normal derivative is nonzero only when $\kappa>0$ .

Proposition 2 is exact, so the ball-on-star transfer contributes no error of its own. One further term is specific to the identity: for nonzero Neumann data the second integrand of Eq. 28 requires $G^{\kappa,B}$ at interior points of $B$ , which is available in closed form but is estimated with a separate sample, since reusing the walk direction gives an estimator that is zero on flat boundaries and high-variance elsewhere (Sawhney et al., 2023). With $h\equiv 0$ that term drops entirely, which is the case throughout this paper.

Original Table 8 — Table 8: Screened Poisson kernels by dimension. All entries follow from Eq. 17, Eq. 18, and Eqs. 20–21 with $\nu=d/2-1$ .

Table 8
$d$	$G_{\mathbb{R}^{d}}(r)$	$G^{D}_{B_{R}}(r)$	$\gamma_{d}(z)$
$2$	$\dfrac{K_{0}(\kappa r)}{2\pi}$	$\dfrac{1}{2\pi}\!\left[K_{0}(\kappa r)-\dfrac{K_{0}(\kappa R)}{I_{0}(\kappa R)}I_{0}(\kappa r)\right]$	$\dfrac{1}{I_{0}(z)}$
$3$	$\dfrac{e^{-\kappa r}}{4\pi r}$	$\dfrac{1}{4\pi r}\!\left[e^{-\kappa r}-e^{-\kappa R}\dfrac{\sinh(\kappa r)}{\sinh(\kappa R)}\right]$	$\dfrac{z}{\sinh z}$
general	$\dfrac{\kappa^{\nu}}{(2\pi)^{d/2}}\dfrac{K_{\nu}(\kappa r)}{r^{\nu}}$	$\dfrac{\kappa^{\nu}}{(2\pi)^{d/2}r^{\nu}}\!\left[K_{\nu}(\kappa r)-\dfrac{K_{\nu}(\kappa R)}{I_{\nu}(\kappa R)}I_{\nu}(\kappa r)\right]$	$\dfrac{(z/2)^{\nu}}{\Gamma(\nu+1)I_{\nu}(z)}$

Table 8 collects the dimension-specific forms. For $d>2$ , $K_{\nu}(z)\sim 2^{\nu-1}\Gamma(\nu)z^{-\nu}$ gives the unscreened limit

G_{\mathbb{R}^{d}}(r)\;\longrightarrow\;\frac{\Gamma(d/2-1)}{4\pi^{d/2}}\,r^{2-d},\qquad\gamma_{d}(\kappa R)\;\longrightarrow\;1,

which is $1/(4\pi r)$ in three dimensions. In two dimensions the free-space Green’s function diverges by an additive constant as $\kappa\to 0$ ; subtracting that constant gives $-(2\pi)^{-1}\log r$ . In the Dirichlet ball kernel the constants cancel, giving $(2\pi)^{-1}\log(R/r)$ . In all dimensions $Q_{d}^{\kappa}(\rho,r)\to 1$ and $\gamma_{d}\to 1$ , recovering the unscreened weights. A purely Neumann unscreened problem needs a separate compatibility and normalization condition.

Assumption 2 (Regularity).

The free domain $\Omega$ is bounded and connected with $C^{2,\alpha}$ boundary, for $0<\alpha<1$ . Its Dirichlet and Neumann sets are disjoint unions of boundary components, each with positive surface measure. Hence they are separated by a positive distance. This excludes Dirichlet–Neumann junctions, near which the solution can have singular gradients. The goal data satisfies $g\in C^{2,\alpha}(\Gamma_{D})$ and $0<g_{\min}\leq g\leq 1$ . The cost satisfies $c\in C^{\alpha}(\overline{\Omega})$ and $0<c_{\min}\leq c\leq c_{\max}$ . The noise parameter $\lambda>0$ is fixed for each problem.

Under Assumption 2, the mixed elliptic problem has a unique $\Psi\in C^{2,\alpha}(\overline{\Omega})$ (Gilbarg and Trudinger, 2001). These regularity statements concern the ideal smooth domain. Voxel edges and corners are not covered by this assumption; the smooth-boundary factor $\alpha=1/2$ and global gradient bounds do not automatically extend to them. The exact PDE identities for the chain additionally require constant screening and admissible stars; regularity does not guarantee that the radius-clamped sampler is exact. Write $\mathcal{L}_{\kappa}=-\Delta+\kappa^{2}$ .

Where $V$ is $C^{2}$ it satisfies the interior stationary HJB equation

0=\min_{\mathbf{u}}\Big\{\nabla V^{\top}(\mathbf{b}+G\mathbf{u})+q+\tfrac{1}{2}\mathbf{u}^{\top}R\,\mathbf{u}\Big\}+\tfrac{1}{2}\operatorname{tr}\!\big(\Sigma\nabla^{2}V\big),

whose inner minimization is attained at $\mathbf{u}^{\star}=-R^{-1}G^{\top}\nabla V$ . Substituting that minimizer leaves a quadratic nonlinearity in $\nabla V$ ,

0=\nabla V^{\top}\mathbf{b}-\tfrac{1}{2}\nabla V^{\top}GR^{-1}G^{\top}\nabla V+q+\tfrac{1}{2}\operatorname{tr}\!\big(\Sigma\nabla^{2}V\big),

and this term cancels under the desirability transform.

Lemma 2 (Cole–Hopf linearization is exact precisely under matching).

Let $V\in C^{2}$ satisfy Eq. 29 and set $\Psi=e^{-V/\lambda}$ . Then

\tfrac{1}{2}\operatorname{tr}(\Sigma\nabla^{2}\Psi)+\mathbf{b}^{\top}\nabla\Psi-\tfrac{q}{\lambda}\Psi=\frac{1}{2\Psi}\,\nabla\Psi^{\top}\Big[\Sigma-\lambda\,GR^{-1}G^{\top}\Big]\nabla\Psi.

In particular the right side vanishes for every $\Psi$ with $\nabla\Psi\neq 0$ if and only if $\Sigma=\lambda GR^{-1}G^{\top}$ as quadratic forms, which is Assumption 1.

Proof.

From $\Psi=e^{-V/\lambda}$ we get $\nabla V=-\lambda\nabla\Psi/\Psi$ and $\nabla^{2}V=-\lambda\frac{\nabla^{2}\Psi}{\Psi}+\lambda\frac{\nabla\Psi\otimes\nabla\Psi}{\Psi^{2}}$ Substitute into Eq. 29 term by term. The diffusion term becomes

\tfrac{1}{2}\operatorname{tr}(\Sigma\nabla^{2}V)=-\frac{\lambda}{2}\frac{\operatorname{tr}(\Sigma\nabla^{2}\Psi)}{\Psi}+\frac{\lambda}{2}\frac{\nabla\Psi^{\top}\Sigma\nabla\Psi}{\Psi^{2}}

and the quadratic control term becomes

-\tfrac{1}{2}\nabla V^{\top}GR^{-1}G^{\top}\nabla V=-\frac{\lambda^{2}}{2}\frac{\nabla\Psi^{\top}GR^{-1}G^{\top}\nabla\Psi}{\Psi^{2}}.

Adding the two quadratic contributions gives

\frac{\lambda}{2\Psi^{2}}\,\nabla\Psi^{\top}\big[\Sigma-\lambda GR^{-1}G^{\top}\big]\nabla\Psi,

which is Eq. 30 up to the stated factor. The remaining terms are $-\lambda\mathbf{b}^{\top}\nabla\Psi/\Psi+q$ and $-\tfrac{\lambda}{2}\operatorname{tr}(\Sigma\nabla^{2}\Psi)/\Psi$ ; multiplying the whole equation by $-\Psi/\lambda$ yields the stationary backward Kolmogorov equation

\tfrac{1}{2}\operatorname{tr}\!\big(\Sigma\nabla^{2}\Psi\big)+\mathbf{b}^{\top}\nabla\Psi-\frac{q}{\lambda}\Psi=0

exactly when Eq. 31 vanishes.

For universal cancellation of the transformed differential expression, put $E=\Sigma-\lambda GR^{-1}G^{\top}$ . The symmetric quadratic form $v^{\top}Ev$ vanishes for every $v$ if and only if $E=0$ . Smooth positive test fields can have an arbitrary gradient at a point, so matching is necessary for this operator identity. For a particular solution $\Psi$ , the residual can vanish even when matching fails. ∎

For the fully actuated reflected model under Assumption 2, let $\Psi$ solve the mixed problem and set $V=-\lambda\log\Psi$ and $\mathbf{u}^{*}=-\nabla V$ . Assume admissible controls, including $\mathbf{u}^{*}$ , give well-defined dynamics with finite expected goal-hitting time and control energy. Completing the HJB square and applying Itô’s formula up to $\tau_{D}\wedge T$ , then taking $T\to\infty$ , yields

J(x,\mathbf{u})=V(x)+\frac{1}{2}\mathbb{E}_{x}\!\int_{0}^{\tau_{D}}\|\mathbf{u}_{t}-\mathbf{u}^{*}(X_{t})\|^{2}\,\mathrm{d}t.

The Neumann condition removes the reflection term; bounded $V$ and $V=\phi$ at the goal justify the stopped limit. Thus $J(x,\mathbf{u})\geq V(x)$ , with equality for $\mathbf{u}^{*}$ , proving optimality.

Part (ii) below fixes the scale on which $\log\Psi$ varies. For constant $c$ the screening length is $\lambda/c$ , but additional length scales arise from the domain size, goal size, and obstacle geometry, so exponential decay with distance holds only asymptotically.

Proposition 3 (Positivity and range).

Under Assumption 2:

(i)
$0<\Psi\leq 1$ on $\overline{\Omega}$ , with $\Psi_{\min}=\min\Psi>0$ ;
(ii)
for the fixed annulus $a<|x|<D$ , $0<a<D$ , constant $\kappa$ , $\Psi(a)=1$ , and $\Psi^{\prime}(D)=0$ , the log range is $-\log\Psi(D)=\kappa(D-a)+O(1)$ as $\kappa\to\infty$ .

Proof.

For (i), the maximum principle excludes an interior maximum above one or a nonpositive interior minimum. At a reflecting boundary, the Hopf boundary lemma (Gilbarg and Trudinger, 2001) excludes either extremum unless the solution is constant. A constant violating the bounds contradicts the positive Dirichlet data. Continuity and compactness then give a strictly positive minimum.

For (ii), put $\nu=d/2-1$ . The radial solution is

\Psi(r)=\frac{r^{-\nu}[K_{\nu}(\kappa r)+AI_{\nu}(\kappa r)]}{a^{-\nu}[K_{\nu}(\kappa a)+AI_{\nu}(\kappa a)]},\qquad A=\frac{K_{\nu+1}(\kappa D)}{I_{\nu+1}(\kappa D)}.

Integrating $(r^{d-1}\Psi^{\prime})^{\prime}=\kappa^{2}r^{d-1}\Psi$ backward from $\Psi^{\prime}(D)=0$ shows $\Psi^{\prime}<0$ for $a<r<D$ . Thus the range is $-\log\Psi(D)$ . At fixed $a,D$ , the large-argument Bessel asymptotics give $\Psi(D)=2(a/D)^{(d-1)/2}e^{-\kappa(D-a)}(1+O(\kappa^{-1}))$ , which proves (ii). ∎

Reflection lets the process remain in the free domain and continue toward the goal after reaching an obstacle boundary. Its corresponding condition is $\partial_{n}\Psi=0$ , as in Eq. 4 and the reflecting walks of Sawhney et al. (2023). This choice also keeps $\log\Psi$ , and so the control, finite up to the boundary under Assumption 2.

Suppose instead that obstacles are absorbing, with zero terminal desirability. Denote the resulting solution by $\Psi_{\mathrm{abs}}$ , retaining positive goal data and reflection at the outer workspace boundary. Then

\Psi_{\mathrm{abs}}=0\text{ on }\partial\mathcal{S},\qquad\log\Psi_{\mathrm{abs}}(x)\to-\infty,\qquad V_{\mathrm{abs}}(x)=-\lambda\log\Psi_{\mathrm{abs}}(x)\to+\infty

as $x$ approaches an obstacle from the free domain. Thus zero desirability represents an infinite collision penalty. The singularity is in the logarithmic representation and value function; the linear PDE solution $\Psi_{\mathrm{abs}}$ can remain smooth.

The control also becomes singular at a smooth obstacle surface. With the smoothness and separation assumptions of Assumption 2, now replacing the obstacle condition by zero Dirichlet data, fix $z\in\partial\mathcal{S}$ and write $x=z-sn(z)$ , where $s>0$ measures distance into the free domain. The Hopf boundary lemma gives $a_{z}=-\partial_{n}\Psi_{\mathrm{abs}}(z)>0$ . Taylor expansion and $\mathbf{u}^{\star}=-\nabla V$ , for $G=R=I$ , yield

\begin{split}\Psi_{\mathrm{abs}}(z-sn(z))&=a_{z}s+O(s^{2}),\\ V_{\mathrm{abs}}(z-sn(z))&=-\lambda\log(a_{z}s)+O(s),\\ \mathbf{u}_{\mathrm{abs}}^{\star}(z-sn(z))&=-\frac{\lambda}{s}n(z)+O(1),\qquad s\downarrow 0.\end{split}

The leading control points away from the obstacle and grows as inverse distance. A bounded log-output model cannot approximate this target uniformly all the way to the obstacle surface. Clipping the desirability to a positive floor removes the divergence numerically but introduces a boundary approximation.

For the reflected problem, Proposition 3 gives $\Psi_{\min}>0$ on the closed free domain. Consequently $V$ and $\log\Psi$ are bounded, $\partial_{n}V=0$ , and $\|\nabla V\|_{\infty}\leq\lambda\|\nabla\Psi\|_{\infty}/\Psi_{\min}<\infty$ under the stated regularity assumptions. The reflected dynamics enforce the boundary constraint; the Neumann condition does not impose a positive clearance from obstacles.

We therefore use reflection because it permits continued goal seeking with a bounded log target; other formulations need not use it. Strictly positive Dirichlet obstacle data $\Psi=e^{-\phi_{\mathrm{obs}}/\lambda}$ with finite $\phi_{\mathrm{obs}}$ avoids this logarithmic singularity, but models termination at an obstacle with a finite penalty. The singularity claim applies specifically to zero Dirichlet data for $\Psi$ , not to Dirichlet conditions in general.

Throughout, $G^{B}$ and $G^{\kappa,B}$ denote the unscreened and screened Dirichlet Green’s functions of $B=B(x,r)$ with pole at the center; both are radial in $\rho=\|\mathbf{y}-x\|$ . By Lemma 1, the step weight is the flux ratio

Q_{d}^{\kappa}(\rho,r)=\frac{\partial_{\rho}G^{\kappa,B}}{\partial_{\rho}G^{B}}.

Both Green’s functions are radial about the center, so the outward normal derivative at $\mathbf{z}$ factors as $\partial_{n}G=(\partial_{\rho}G)(n_{\mathbf{z}}\cdot\hat{e}_{\rho})$ and the geometric factor cancels in the ratio, which is why Eq. 34 is also the ratio of the two Poisson kernels in Eq. 24.

For $d\geq 2$ , $\kappa>0$ and $0<\rho\leq r$ the step weight is bounded on both sides,

\gamma_{d}(\kappa r)\;\leq\;Q_{d}^{\kappa}(\rho,r)\;<\;1,

with $\gamma_{d}$ the full-sphere weight of Eq. 22. By the upper bound, the weighted chain of Section 4.1 remains bounded; the lower bound gives $Q_{d}^{\kappa}\geq\gamma_{d}(\kappa r_{\max})$ when the star radius is at most $r_{\max}$ , which yields the replica bounds used in Appendix C.7.

Bounds on $Q_{d}^{\kappa}$ .

Positivity is immediate from Eq. 25, since $K_{\nu+1}$ , $I_{\nu+1}$ and $K_{\nu}/I_{\nu}$ are all positive. For the upper bound, set $D:=G^{B}-G^{\kappa,B}$ , which is radial and vanishes at $\rho=r$ . Subtracting the two defining equations cancels the point sources and leaves $-\Delta D=\kappa^{2}G^{\kappa,B}$ in $B$ , so in flux form, the constant vanishing because the right side is locally integrable at the origin,

-\rho^{d-1}\partial_{\rho}D=\kappa^{2}\!\!\int_{0}^{\rho}\!t^{d-1}G^{\kappa,B}(t)\,\mathrm{d}t,\qquad-\rho^{d-1}\partial_{\rho}G^{B}=\frac{1}{\omega_{d}}.

Dividing the first by the second gives the weight-deficit identity

1-Q^{\kappa}_{d}(\rho,r)=\frac{\partial_{\rho}D}{\partial_{\rho}G^{B}}=\kappa^{2}\!\!\int_{B(x,\rho)}\!G^{\kappa,B}(x,\mathbf{y})\,\mathrm{d}\mathbf{y},

The integrand is strictly positive on the interior of $B$ , so the right side is strictly positive for $\rho>0$ and $Q^{\kappa}_{d}<1$ . For the lower bound, $B(x,\rho)\subseteq B(x,r)$ and the integrand is nonnegative, so the deficit is at most $\kappa^{2}\!\int_{B}G^{\kappa,B}$ ; integrating $(-\Delta+\kappa^{2})G^{\kappa,B}=\delta_{x}$ over $B$ gives $\kappa^{2}\!\int_{B}G^{\kappa,B}=1-\int_{\partial B}P^{\kappa}_{B}\,\mathrm{d}S=1-\gamma_{d}(\kappa r)$ by Proposition 1. Hence $\gamma_{d}(\kappa r)\leq Q^{\kappa}_{d}(\rho,r)<1$ , with equality on the left at $\rho=r$ . This follows directly from $B(x,\rho)\subseteq B(x,r)$ ; monotonicity in $\rho$ also follows from the same identity but is not needed. ∎

Alternatively, $\log\Psi_{\theta}$ can be fitted to log labels. There the weights are accumulated as $\Lambda^{(n)}=\sum_{j<\sigma^{(n)}\wedge K}\log W_{j}^{(n)}$ plus the log of the terminal value, which avoids underflow for long products, but $\log\widehat{\Psi}$ is a biased estimate of $\log\bar{\Psi}$ , since the logarithm of a sample mean is not the mean of logarithms. The correction below removes that bias to $O(N^{-2})$ .

For i.i.d. replicas with fixed bounds $0<z_{\min}\leq e^{\Lambda^{(n)}}\leq z_{\max}<\infty$ and variance $s^{2}$ ,

\mathbb{E}[\log\widehat{\Psi}]=\log\bar{\Psi}-\frac{s^{2}}{2N\bar{\Psi}^{2}}+O(N^{-2}).

The constants can depend on the replica bounds. For a fixed environment, such bounds follow from a finite horizon and a bounded target log field (Appendix C.7).

For $N\geq 2$ , compute all ratios from shifted exponentials. Set $b=\max_{n}\Lambda^{(n)}$ , $z_{n}=e^{\Lambda^{(n)}-b}$ , $\bar{z}=N^{-1}\sum_{n}z_{n}$ , and $s_{z}^{2}=(N-1)^{-1}\sum_{n}(z_{n}-\bar{z})^{2}$ . Then

\widehat{v}=\frac{s_{z}^{2}}{N\bar{z}^{2}},\qquad\widehat{\Lambda}=\log\widehat{\Psi}+\frac{\widehat{v}}{2}=b+\log\bar{z}+\frac{\widehat{v}}{2}.

Here $s_{z}^{2}$ is the sample variance of the shifted replicas $z_{n}$ . The random common shift makes those replicas dependent; unbiasedness is invoked only for the unshifted replicas below. Writing $\widehat{s}^{2}$ for the unbiased sample variance of the unshifted replicas $e^{\Lambda^{(n)}}$ , the ratio is unchanged by the common shift, $\widehat{v}=s_{z}^{2}/(N\bar{z}^{2})=\widehat{s}^{2}/(N\widehat{\Psi}^{2})$ , so Eq. 36 adds the bias correction $\widehat{s}^{2}/(2N\widehat{\Psi}^{2})$ . The correction removes the leading log bias, leaving $O(N^{-2})$ under the stated assumptions.

Proposition 4 (Bias of a log sample mean).

Let $Z_{1},\ldots,Z_{N}$ be i.i.d. with $Z_{i}\in[z_{\min},z_{\max}]$ , $0<z_{\min}\leq z_{\max}<\infty$ , mean $\mu$ , and variance $s^{2}$ . For $N\geq 2$ , write $\widehat{\Psi}=N^{-1}\sum_{i}Z_{i}$ and let $\widehat{s}^{2}$ be the unbiased sample variance. Then

\mathbb{E}[\log\widehat{\Psi}]=\log\mu-\frac{s^{2}}{2N\mu^{2}}+O(N^{-2}),\qquad\mathbb{E}\!\left[\log\widehat{\Psi}+\frac{\widehat{s}^{2}}{2N\widehat{\Psi}^{2}}\right]=\log\mu+O(N^{-2}).

The constants depend on the replica bounds.

Proof.

Set $U=\widehat{\Psi}-\mu$ . Taylor’s theorem through degree three gives $\log(\mu+U)=\log\mu+U/\mu-U^{2}/(2\mu^{2})+U^{3}/(3\mu^{3})+R_{4}$ with $|R_{4}|\leq CU^{4}$ because $\widehat{\Psi}\geq z_{\min}$ . Independence gives $\mathbb{E}U=0$ , $\mathbb{E}U^{2}=s^{2}/N$ , $\mathbb{E}U^{3}=O(N^{-2})$ , and $\mathbb{E}U^{4}=O(N^{-2})$ , proving the first expansion.

For the correction, expand $(\mu+U)^{-2}=\mu^{-2}-2\mu^{-3}U+R_{2}$ , where $|R_{2}|\leq CU^{2}$ . Since $\widehat{s}^{2}$ is bounded, $\mathbb{E}[\widehat{s}^{2}R_{2}]=O(N^{-1})$ . Also $\mathbb{E}\widehat{s}^{2}=s^{2}$ and $\mathbb{E}[\widehat{s}^{2}U]=m_{3}/N$ , with $m_{3}$ the third central moment of $Z_{i}$ . Hence $\mathbb{E}[\widehat{s}^{2}/\widehat{\Psi}^{2}]=s^{2}/\mu^{2}+O(N^{-1})$ , which proves the corrected expansion. Jensen’s inequality gives the nonpositive bias of the uncorrected logarithm for every $N$ . ∎

For a fixed environment, horizon $K$ , and target bounded by $0<\Psi^{-}_{\min}\leq\Psi_{\theta^{-}}\leq\Psi^{-}_{\max}<\infty$ , let $W_{\min}>0$ be a lower bound on the step weight, which exists by Eq. 35 when the star radius is bounded. The walk replicas then lie in $[z_{\min},z_{\max}]$ with $z_{\min}=W_{\min}^{K}\min(g_{\min},\Psi^{-}_{\min})$ and $z_{\max}=\max(1,\Psi^{-}_{\max})$ . The target bounds hold because the log field is continuous on the compact domain. They are not assumed uniform over all trained parameters or environments. The shifted calculation in Eq. 36 preserves the variance ratio exactly in real arithmetic and avoids underflow from a common small scale.

Recall the lattice of Section 4.4, with vertices $\mathcal{X}_{h}$ and edges between axis-adjacent free cells. We take this graph to be finite and $\Psi_{\theta}$ finite at every vertex.

The largest desirability reachable in $k$ moves is computed by

B_{0}(x)=\Psi_{\theta}(x),\qquad B_{k+1}(x)=\max_{y\in\{x\}\cup\mathcal{N}(x)}B_{k}(y).

Indeed, $\mathcal{R}_{0}(x)=\{x\}$ and

\mathcal{R}_{k+1}(x)=\mathcal{R}_{k}(x)\cup\bigcup_{y\in\mathcal{N}(x)}\mathcal{R}_{k}(y),

so induction gives $B_{k}(x)=\max_{z\in\mathcal{R}_{k}(x)}\Psi_{\theta}(z)$ . Equation 14 therefore compares $B_{L-1}$ at candidate neighbors. Moving to that neighbor uses the first of the $L$ moves.

In each propagation layer, the previous score is kept unless a neighbor has a strictly higher score, in which case the corresponding edge is stored. Starting at horizon $L$ , follow the stored edges while decreasing the remaining horizon. A layer with no edge is skipped because its maximum was already reachable at a shorter horizon. Induction on the remaining horizon gives a feasible route of at most $L$ edges ending at a cell with score $B_{L}(x_{0})$ . Executing a prefix preserves graph feasibility, although the next replan may select a different endpoint. Continuous-motion feasibility additionally requires valid edges; collision-free sampled endpoints are insufficient. The implemented lattice uses cells marked free in the saved occupancy mask, connected through face-adjacent free cells. Continuous segments from off-center states are checked across all traversed voxels. All edges are checked for collisions. For $e=1$ , Eq. 14 instead selects a neighbor anew from the fixed $B_{L-1}$ map at each visited state.

The reachable sets are nested, hence $B_{k+1}(x)\geq B_{k}(x)$ . On a connected component $C$ of graph diameter $D$ , every $k\geq D$ gives $B_{k}(x)=\max_{z\in C}\Psi_{\theta}(z)$ for all $x\in C$ . Once $L-1\geq D$ , the candidate scores in Eq. 14 coincide and direction depends on tie-breaking. For a reconstructed route, maximizing the terminal score likewise does not distinguish the lengths of paths to the same endpoint. Thus increasing depth enlarges the set of candidate endpoints, but does not guarantee a shortest path or a successful rollout.

The pretrained 2D application benchmark contains 592 maps (37,888 starts); the 2D generated and 4D manipulator test sets each contain 512 maps (32,768 starts). Every map uses eight goals and every goal uses eight starts, chosen at the centers of eight equal goal-distance quantile bins. Success is evaluated on the matched saved queries, including failed rollouts in the denominator. Field prediction or solve, controller preprocessing, and rollout/search all contribute to reported test-time cost. These costs are amortized over the corresponding starts; they include more than the neural forward time. FEM solutions provide the numerical ground-truth fields for reconstruction-error evaluation.

Table 5 compares nominal planners, sampling-based controllers, and PUN at execution-noise levels $\lambda\in\{0.25,1,2,4\}$ . Before boundary handling, the proposed transition is $\widetilde{x}_{k+1}=x_{k}+a_{k}+\sqrt{\lambda\Delta t}\,z_{k}$ , with independent $z_{k}\sim\mathcal{N}(0,I_{d})$ , $\Delta t=1/65{,}536$ , and $\|a_{k}\|_{2}\leq 1/(8n)$ , where $n$ is the largest grid dimension. Each control interval contains eight microsteps. Coordinates occupy a unit square in 2D, a $1\times 1\times\tfrac{1}{2}$ domain in 3D, and normalized joint coordinates in 4D. The noise-increment covariance is $\lambda\Delta tI_{d}$ , consistent with the PDE specialization $\Sigma=\lambda I_{d}$ . In these replanning evaluations, all compared methods use the same execution and collision-checking rules. Collision checks cover the entire executed motion segment, and an obstacle hit is counted as a rollout failure. Thus, success requires reaching the goal without an obstacle hit or an outer-domain boundary violation. The reflected PDE provides the planning field, whose effectiveness is evaluated here under collision-failure execution. Joint limits are not modeled in these experiments. Moves outside the outer workspace or configuration-domain boundary are also counted as rollout failures. For each environment–goal query at a given noise level, PUN predicts the entire field in a single forward pass and reuses it throughout the rollout. The reported replanning comparisons use positive noise levels, consistent with the $\lambda>0$ scope of the screened-Poisson derivation. Success requires entry into the saved goal-cell mask before failure. Methods replan from the current noisy state every eight microsteps. Execution also stops on a collision failure, an unavailable route, or the controller’s stopping condition.

The MPPI and PI $^{2}$ -CMA implementations share a horizon of $H=32$ steps, $K=128$ sampled trajectories per optimization iteration, and one iteration per planning call. Both use double precision, componentwise action bounds $[-0.05,0.05]$ , initial exploration standard deviation $0.025$ , and optimization temperature $\eta=0.5$ . This temperature is distinct from the execution-noise parameter $\lambda$ . A fresh call initializes every action to $(x_{\mathcal{G}}-x)/H$ , where $x_{\mathcal{G}}$ is the goal center. Thus each call evaluates $KH=4{,}096$ sampled transitions, followed by a rollout of the selected nominal sequence. The wrapper supports resuming interrupted optimization from saved parameters and random-generator state; ordinary calls initialize a new optimizer rather than carrying its solution between replans. The configuration uses $(H,K,I)=(32,128,1)$ , where $I$ is the iteration count. These factors describe nominal sampled work $IKH$ . No validation-based parameter sweep was performed. We use standard values for temperature, exploration scale, action bounds, the objective weights.

A⋆ and Dijkstra use standard C++ implementations with face-neighbor connectivity and unit edge costs. Grids are the native 2D maps, $64\times 64\times 32$ in 3D, and $17^{4}$ in 4D. A⋆ uses Manhattan distance to the goal-mask bounding box; Dijkstra uses zero heuristic. CPU RRT uses OMPL 2.0.1, goal bias 0.05, uniform domain sampling, and maximum extension 0.4 in normalized coordinates. Search stops on a goal solution; there is no fixed sample or tree-size cap.

Both sampling-based optimizers use the following task cost, with stage costs evaluated after each predicted transition:

\begin{aligned}\displaystyle\ell_{t} & \displaystyle=\mathbf{1}_{\{x_{t+1}\notin\mathcal{G}\}}\left(1+10\|x_{t+1}-x_{\mathcal{G}}\|_{2}^{2}+\|u_{t}\|_{2}^{2}\right)+100I_{t},\end{aligned}

\begin{aligned}\displaystyle\phi_{H} & \displaystyle=100\,\mathbf{1}_{\{x_{H}\notin\mathcal{G}\}}\|x_{H}-x_{\mathcal{G}}\|_{2}^{2}.\end{aligned}

Here $I_{t}$ flags a proposed segment rejected by the geometry-validity routine when its starting state is not already in the goal. Before calling the geometry step routine, the wrapper replaces rejected actions and actions starting in the goal with zero displacement. This is a penalized planning surrogate; executed collisions are failures under the evaluation protocol above. Goal absorption is delegated to the geometry step routine. After optimization, the wrapper returns the predicted state path, truncated at its first recorded goal state. Execution noise is applied outside the optimization wrapper; its action-exploration distribution should not be identified with the execution-noise distribution. Geometry routines traverse all crossed voxels, including incident cells at edge/corner crossings, to detect goal entry, obstacle contact, and domain exit. Collision takes precedence over simultaneous goal entry in first-contact scoring. Actions are displacements in the normalized coordinates specified above. Predicted trajectories use deterministic geometry; the outer rollout adds execution noise. Returned paths are subdivided into bounded microsteps.

We use the MPPI class from pytorch_mppi $^{3}$ (Williams et al., 2017a), with $\texttt{noise\_sigma}=0.025^{2}I_{d}$ , $\texttt{lambda\_}=0.5$ , the task costs above, and the stated action bounds. The one call to command refines the nominal sequence at the same state with shift_nominal_trajectory=False. We retain the library’s MPPI update and control-perturbation cost in addition to the supplied task cost, and generate the returned path from the final nominal action sequence. The implementation uses the original MPPI update and control-perturbation cost.

Our adaptation of PI $^{2}$ -CMA (Stulp and Sigaud, 2012) represents the action sequence by 16 uniformly spaced knots per coordinate, linearly interpolated to the 32-step horizon with aligned endpoints and clipped to the action bounds. A single parameter vector is sampled for each trajectory from $\mathcal{N}(\theta,\Sigma)$ , with initial covariance $\Sigma=0.025^{2}I_{16d}$ . For sampled parameters $\theta_{k}$ , let $S_{kt}$ be the remaining task cost from step $t$ , including the terminal cost once. We compute

P_{kt}=\frac{\exp[-(S_{kt}-\min_{j}S_{jt})/\eta]}{\sum_{j}\exp[-(S_{jt}-\min_{l}S_{lt})/\eta]},\qquad w_{k}=\sum_{t=0}^{H-1}\frac{H-t}{H(H+1)/2}\,P_{kt},

then update $\theta^{+}=\sum_{k}w_{k}\theta_{k}$ and $\Sigma^{+}=\sum_{k}w_{k}(\theta_{k}-\theta)(\theta_{k}-\theta)^{\top}$ , centering the covariance at the old mean. This implements temporal averaging of the per-step updates for exploration that is constant over a trajectory. The implementation symmetrizes the covariance and samples using an eigendecomposition, checking negative eigenvalues against a numerical tolerance before clipping them to zero; it imposes no positive exploration floor. We use a fixed temperature on unnormalized costs, rather than the cost-range normalization used in the original paper’s experiments. The custom wrapper implements knot-based PI $^{2}$ -CMA, without CMA-ES evolution paths, rank-based selection, or step-size adaptation.

Both CPU and CUDA implementations were evaluated; the CPU version was faster in the implementation comparison. The optimizer’s wall-clock timer includes initialization, the one optimization iteration, logging, the final nominal rollout, and conversion of the returned path to host memory. The CUDA implementation also synchronizes after each iteration and before returning. These are per-call measurements; the tables report per-start costs and require the outer execution loop’s accounting across planning calls. Table 5 uses CPU baseline timings in 2D/3D/4D since they are faster than GPU, and GPU PUN timings throughout. Hardware is an Intel Xeon Platinum 8480+ or an NVIDIA B200. CPU software is Python 3.11.11, PyTorch 2.7.1+cpu, and NumPy 1.26.4; GPU software is Python 3.9.21, PyTorch 2.7.1+cu128, and NumPy 2.0.2.

Timing excludes loading and warmup and uses three observations per selected pair, method, and noise level for the benchmark sets. Variability is minimal. Per-start measurements include initialization, all replanning calls, geometry, and execution, including failed rollouts; the number of calls depends on the realized trajectory. GPU work is synchronized at timing boundaries.

PUN timings include fresh field prediction and preparation for each start, with the field reused during replanning.

Table 5 uses 500 maps in 2D, 250 in 3D, and 480 in 4D, with three goals per map and one start per goal. Each query uses one seeded execution-noise combination at each noise level; methods share the saved queries and Gaussian noise sequence. Rollouts receive equal weight, which also gives maps equal weight because each contributes three queries. One trained model per dimension is evaluated. We compute 95% intervals using 50,000 paired map-level bootstrap resamples, retaining all queries within each sampled map. Nondegenerate intervals are the estimate $\pm 1.96$ bootstrap standard errors. Runtime intervals use separate samples of maps, respectively.

Original Table 9 — Table 9: 95% confidence intervals for replanning success rates (%) corresponding to Table 5.

Table 9
Dim.	Method	$\lambda=0.25$	$\lambda=1$	$\lambda=2$	$\lambda=4$
2D	A*	[55.92, 62.08]	[39.16, 45.50]	[31.90, 38.10]	[20.89, 26.45]
Dijkstra	[55.81, 62.19]	[39.15, 45.51]	[31.99, 38.01]	[21.08, 26.26]
RRT	[79.37, 84.63]	[65.68, 71.66]	[49.14, 55.52]	[33.46, 39.88]
MPPI	[81.70, 86.30]	[57.77, 64.23]	[42.45, 48.89]	[28.82, 34.52]
PI $^{2}$ -CMA	[75.76, 80.90]	[56.81, 63.19]	[44.76, 51.24]	[31.95, 38.05]
PUN (zero-shot)	[89.90, 93.44]	[71.27, 76.73]	[59.15, 65.51]	[42.55, 48.79]
3D	A*	[32.37, 40.97]	[19.38, 27.28]	[9.10, 14.90]	[4.04, 7.96]
Dijkstra	[33.00, 41.66]	[18.98, 26.36]	[9.05, 14.95]	[3.73, 8.27]
RRT	[48.85, 57.81]	[30.86, 39.80]	[21.28, 29.38]	[8.69, 15.31]
MPPI	[58.44, 66.90]	[30.43, 38.91]	[14.98, 22.36]	[4.93, 9.73]
PI $^{2}$ -CMA	[42.92, 51.74]	[18.36, 25.64]	[17.51, 25.15]	[6.07, 11.27]
PUN (zero-shot)	[89.60, 94.40]	[64.44, 72.90]	[35.62, 44.38]	[12.71, 19.29]
4D	A*	[100.00, 100.00]	[95.39, 97.95]	[69.54, 75.46]	[55.94, 62.40]
Dijkstra	[100.00, 100.00]	[95.68, 97.66]	[69.66, 75.34]	[56.06, 62.28]
RRT	[97.53, 99.13]	[94.65, 97.01]	[87.97, 92.03]	[83.82, 87.84]
MPPI	[93.68, 96.32]	[90.59, 94.41]	[79.15, 84.19]	[68.91, 74.43]
PI $^{2}$ -CMA	[95.43, 97.91]	[90.70, 94.30]	[84.40, 88.94]	[75.77, 80.89]
PUN	[100.00, 100.00]	[100.00, 100.00]	[100.00, 100.00]	[93.49, 96.51]

Original Table 10 — Table 10: 95% confidence intervals for replanning runtime (ms/start) corresponding to Table 5.

Table 10
Dim.	Method	$\lambda=0.25$	$\lambda=1$	$\lambda=2$	$\lambda=4$
2D	A*	[24.99, 62.27]	[26.25, 67.93]	[28.10, 74.26]	[29.44, 78.08]
Dijkstra	[79.37, 230.51]	[88.65, 261.25]	[102.34, 304.48]	[123.40, 356.54]
RRT	[69.41, 127.49]	[72.10, 135.64]	[78.18, 146.36]	[82.46, 166.96]
MPPI	[1227.92, 2240.30]	[1257.88, 2288.74]	[1279.80, 2397.66]	[1400.54, 2851.70]
PI $^{2}$ -CMA	[964.06, 1697.46]	[1053.67, 1916.15]	[1067.48, 1981.92]	[1079.98, 2133.66]
PUN (zero-shot)	[7.79, 11.31]	[8.80, 13.48]	[10.71, 17.43]	[10.19, 18.83]
3D	A*	[50.74, 96.06]	[34.34, 86.60]	[42.19, 105.51]	[40.32, 90.18]
Dijkstra	[50.08, 96.84]	[35.29, 90.79]	[43.41, 110.57]	[41.51, 93.95]
RRT	[155.27, 321.35]	[93.18, 246.62]	[76.54, 169.16]	[94.56, 157.90]
MPPI	[382.10, 682.52]	[408.91, 746.77]	[433.13, 675.09]	[526.69, 881.93]
PI $^{2}$ -CMA	[332.42, 572.08]	[305.61, 564.41]	[371.44, 580.52]	[482.06, 731.04]
PUN (zero-shot)	[14.95, 16.33]	[14.61, 16.07]	[13.26, 15.34]	[13.69, 15.17]
4D	A*	[12.32, 21.30]	[12.25, 21.91]	[35.26, 120.78]	[28.18, 96.70]
Dijkstra	[12.51, 21.83]	[12.37, 22.31]	[31.21, 102.29]	[27.73, 91.05]
RRT	[16.26, 32.14]	[16.20, 35.40]	[23.09, 62.91]	[28.09, 48.53]
MPPI	[332.13, 490.35]	[329.20, 467.16]	[297.34, 431.10]	[300.44, 449.22]
PI $^{2}$ -CMA	[300.61, 458.97]	[291.21, 424.17]	[284.19, 426.69]	[288.82, 424.62]
PUN	[10.50, 10.76]	[10.32, 10.70]	[10.37, 10.73]	[11.43, 11.89]

The experiments use the following estimator and numerical choices.

1.
Absorption and radius floor. A walk reaches $\Gamma_{D}$ only in the limit, so absorption is triggered on the shell $\{d_{D}<\varepsilon\}$ and the terminal value is $g$ at the nearest point of $\Gamma_{D}$ ; $\varepsilon$ is chosen below the tubular-neighborhood radius of $\Gamma_{D}$ and its distance to $\Gamma_{N}$ . The star radius also carries a floor, $r(x)=\min\{d_{D}(x),\max(r_{\min},d_{\mathrm{sil}}(x))\}$ , so that steps do not become arbitrarily short near a silhouette. The shell costs at most $C_{\Psi}\varepsilon$ (Eq. 38); the floor can break star-shapedness, and the error it introduces has unknown sign (Sawhney et al., 2023).
2.
Weighted survival. The solid-angle sampling law in Eq. 8 leaves the scalar screening factor $Q_{d}^{\kappa}$ as the per-step weight. We multiply these factors rather than sample killing events or apply Russian roulette. Conditioning out the killing randomness removes that source of variance and avoids zero replicas caused by stochastic killing.
3.
Fixed-horizon continuation. Walks stop at the prescribed step cap or the goal. Capped paths contribute their accumulated weight times the detached EMA prediction at the endpoint, as in Section 4.2. Goal-reaching paths use the known boundary value. Every requested trajectory remains in the sample mean, including capped paths.
4.
Log accumulation and bias correction. The trace records the sum of log step weights independently of their floating-point product, so product underflow does not erase the log evidence. The log-label variants use log-sum-exp and the sample-variance correction in Eq. 36. The direct-label variants use the weighted arithmetic mean without this correction. Direct and log targets are constructed from the same sampled paths.
5.
Reuse of model-independent walks. Cached endpoints, weights, and completion states are independent of the neural parameters. The EMA continuation is reevaluated during training rather than frozen into the cache.
6.
Stable screening kernels. The voxel backend for dimensions three through seven evaluates the full-sphere factor using the positive series
$\gamma_{d}(z)=\left[\sum_{j=0}^{\infty}a_{j}\right]^{-1},\qquad a_{0}=1,\qquad a_{j}=a_{j-1}\frac{z^{2}}{4j(\nu+j)},\qquad\nu=\frac{d}{2}-1.$
The CUDA implementation sums in double precision and uses a logarithmic large-argument asymptotic branch. For interior reflecting hits, it evaluates Eq. 25 with $e^{-z}I_{\nu}(z)$ and $e^{z}K_{\nu}(z)$ , combining nonnegative terms with an exponential factor no greater than one. With $t=\kappa\rho$ , $T=\kappa r$ , and $\delta=T-t$ , the near-sphere branch uses
$Q_{d}^{\kappa}(\rho,r)\approx\gamma_{d}(T)\left[1+\frac{\delta^{2}}{2}-\frac{d-1}{6T}\delta^{3}\right].$
This reproduces the sphere endpoint without two separately rounded Bessel evaluations. The CUDA branches use this expansion for $\delta\leq 2\cdot 10^{-6}(1+T)$ and the limiting value one for $t\leq 10^{-7}$ or $T=0$ .
7.
Voxel boundary queries. The same backend uses bounding-volume hierarchies for goal faces and view-dependent reflecting silhouettes, and voxel traversal for first boundary hits. This avoids triangulating the voxel boundary. Grid-plane and corner handling preserves the appropriate incident free cell and resolves boundary hits before accepting a sphere endpoint. These geometric operations implement the first-hit chain in Eq. 7.

Let $\sigma_{\varepsilon}=\inf\{k:d_{D}(\xi_{k})<\varepsilon\}$ and $\tau=\sigma_{\varepsilon}\wedge K$ for finite $K$ . Assume admissible active steps, the smooth-domain regularity assumptions, and an $\varepsilon$ below the goal’s tubular-neighborhood radius and its separation from reflecting boundaries. Define $\Psi_{\varepsilon,K}(x)$ as the expected product of weights through $\tau$ times $g(\pi_{D}\xi_{\tau})$ if $\sigma_{\varepsilon}\leq K$ , or the true continuation $\Psi(\xi_{K})$ otherwise. Then

\|\Psi_{\varepsilon,K}-\Psi\|_{\infty}\leq C_{\Psi}\varepsilon,\qquad C_{\Psi}=\sup_{\overline{\Omega}}\|\nabla\Psi\|.

To prove this, bounded stopping gives $\Psi(x)=\mathbb{E}_{x}[(\prod_{j<\tau}W_{j})\Psi(\xi_{\tau})]$ . Only the shell event changes the terminal value. Since $g(\pi_{D}\xi_{\tau})=\Psi(\pi_{D}\xi_{\tau})$ , the projection segment lies in the closed free domain, and its length is below $\varepsilon$ , the terminal difference is bounded by $C_{\Psi}\varepsilon$ . The product of weights is at most one. This proves the bound without an infinite-walk termination assumption. If $\sigma_{\varepsilon}<\infty$ almost surely for each start, bounded convergence also yields the same bound for the uncapped shell estimator. Replacing true continuation by a fixed bounded target adds its discounted endpoint error on $\{\sigma_{\varepsilon}>K\}$ , as described in Section 4.2. The radius floor is not covered: a clamped step can violate admissibility and introduce an error of unknown sign (Sawhney et al., 2023).

Antithetic directions, first-step gradient control variates, weight-floor stopping, expected-clock stopping, and randomized horizons are disabled in the reported pretrained models. In particular, the full-sphere gradient control variate is not transferred to arbitrary WoSt stars.

FNO uses 200 epochs. Each alternative backbone has a maximum of 200 epochs and a registered training-compute allowance of 66,856.484 seconds; an epoch is not started if it is not expected to fit within the remaining allowance. Each completed epoch uses 64 training maps and eight goals per map. The 200-epoch cosine schedule is retained when a run stops at the compute cap. The allowance was derived from 1.5 times an estimated individual FNO training cost; shared work in the reference run means it is not an exact equal-GPU-time comparison. Completed training epochs were 200 for FNO, 200 for U-NO, 185 for U-Net, 13 for GINO, 200 for LocalNO, 200 for DISCO, 146 for UPT, and 134 for Transolver. Our overall method is architecture-agnostic; we chose FNO because it gives the highest planning success among the tested implementations under this compute policy, but our overall method is compatible with other neural-operator architectures. Planning benefits from balancing sensitivity to local obstacle geometry with a global view of the environment.

The primary WoSt rows fit direct desirability with value, goal, and gradient loss terms. The loss ablations use 200 epochs and remove one term at a time. All predictions are evaluated through the stated controller; a small field error does not guarantee path success.

Unless a depth is stated, Local search denotes Local search (10). The finite-lattice controllers use legal axis-adjacent moves, with cycle detection and a step cap. Greedy uses $(L,e)=(1,1)$ and requires strict improvement. Local search (4) uses $(L,e)=(5,1)$ in 2D, 3D, and 4D. Local search (10) uses $(L,e)=(10,3)$ with route reconstruction in 2D, and $(L,e)=(11,1)$ in 3D and 4D. In the single-move variants, four or ten propagation steps are counted after the candidate first move, explaining the extra edge in $L$ . Direct models use their raw finite desirability predictions, without clipping or enforcing a goal score during search. For the convention $\ell_{\theta}=\log\Psi_{\theta}$ in Appendix C.6, the maximizing search score is $\ell_{\theta}$ . If a model instead predicts $-\log\Psi_{\theta}$ or $V_{\theta}$ , its output must be negated to preserve this ranking (for $\lambda>0$ in the latter case). The implemented log-label model predicts $\ell_{\theta}=-\operatorname{softplus}(z_{\theta})$ , reconstructs desirability as $\exp(\ell_{\theta})$ , and uses $-z_{\theta}$ as its search score. Maximizing $-z_{\theta}$ gives the same ranking as maximizing $\ell_{\theta}$ , consistent with this convention. Success requires entering the saved goal mask. Ties in the single-move rule favor the candidate’s own field score, then fewer previous visits, then a fixed lexicographic coordinate order (ascending in 2D, descending in 3D/4D). Route reconstruction retains an existing choice on equal scores, using a fixed axis traversal order for competing improvements. The 2D ten-step controller stops when the best endpoint does not improve the current score by more than $10^{-12}$ . Missing legal moves, a fifth visit to a cell under local search, or exhaustion of the rollout budget are recorded as failures. The local diagrams below use one stored field to display these rules, with no new benchmark outcomes inferred from the illustration.

Details of the walk-horizon study and qualitative resolution-transfer examples are provided below.

The controlled horizon study in Table 4 uses 16 training maps, eight goals per map, four trajectories per query, and 40 epochs, with evaluation on the generated 2D test set. With the intermediate 100- and 400-step settings, gradient / local-search (10) success is 40.55% / 66.80% and 38.06% / 69.69%, respectively. Local-search success is highest at 400 steps (69.69%) and 0.03 percentage points lower at 800 steps, whereas gradient success is highest at 100 steps. A separate study of budget allocation without pretraining is reported in Appendix F.1.

Original Table 11 — Table 11: FNO (WoSt) planning success (%) with gradient following and local search.

Table 11
Test set	Gradient	Local search (10)
3D application benchmarks	65.95	89.46
4D manipulators	97.50	99.48

Original Table 12 — Table 12: Full test-time cost (ms/start) on the 2D zero-shot application-inference subset and 4D manipulator test set. Both methods use local search in 2D and in 4D.

Table 12
Method	2D zero-shot	4D manipulators
WoSt + FNO	1.20	2.79
Direct WoSt	1128.60	160.59

Figure 4 shows FNO predictions on held-out maps at $128\times 128$ , $200\times 200$ , and $256\times 256$ . Each field supports paths from multiple starts to one goal. These examples illustrate field reuse at different discretizations; the aggregate comparison across architectures and larger grids is Fig. 1(b).

Fig. 4: Zero-shot FNO field predictions and extracted paths at $128\times 128$ , $200\times 200$ , and $256\times 256$ (left to right). Color shows predicted desirability $\Psi_{\theta}$ , gray denotes obstacles, circles mark starts, and the star marks the shared goal in each map.

AI4Arctic Sea Ice (Arctic sea-ice maps from satellite observations); Fields2Benchmark (agricultural field geometries); HiRISE DTM (terrain models of the surface of Mars); HydroLAKES ASV (lake-shoreline geometries for autonomous surface vessels); LOLA LDEM (lunar elevation maps from laser altimetry); MovingAI BG512 (Baldur’s Gate maps); MovingAI Street (urban street-map grids); NOAA ENC Harbor (harbor geometries from navigational charts); Open Forest Observatory (forest structure maps from drone imagery); PNO2D (synthetic obstacle environments); TartanDrive 2 (terrain from off-road driving data).

Dryad SIC Foam (porous foam geometries); Falkor MAR AUV (Mid-Atlantic Ridge seafloor for autonomous underwater vehicles); Flybo UAV (scenes for aerial-robot navigation); HVaCT Liver Tip (liver and vessel anatomy for navigation); Industrial Pipe Routing (industrial layouts for pipe routing); Lung Needle (lung anatomy for needle-path planning); MOSAiC Under-Ice AUV (Arctic under-ice geometry for autonomous underwater vehicles); PNO3D iGibson/HouseExpo (indoor scene and floor-plan geometries); ReplicaCAD Indoor UAV (indoor CAD scenes for aerial navigation); Trabecular Bone (porous internal bone structure); Underwater AT40-02 MAR (seafloor surveyed for the Mid-Atlantic Ridge); Axial North Rift (seafloor terrain at Axial Seamount’s north rift); Escanaba (seafloor terrain in Escanaba Trough); Loihi Pele’s Pit (a submarine volcanic crater); Squid5 EDR (Eastern Dry Rocks coral reef topography from underwater imagery); USGS Tahoe (Lake Tahoe lakebed topography).

eCobra600 (Omron industrial robot arm); Kinova MICO (robot arm); OpenManipulator (ROBOTIS open-source robot arm); qArm (Quanser robot arm). Each 4D state represents a manipulator configuration in joint coordinates.

Shifting from 64 walks per query and 8 goals per map to 4 walks and 128 goals reduces 2D test success from 62.6% to 51.5% for DISCO and from 61.6% to 56.3% for FNO. In 3D, the same shift raises FNO success from 25.0% to 28.8%; Transolver success changes from 6.3% to 7.0%, while label-generation cost rises from 244.35 to 621.46 GPU hours. Equal nominal label budgets therefore need not imply equal computation or a common preferred allocation.

The following examples illustrate the geometry-to-field mapping, stochastic labels, and path extraction. The sampling budget for each field and walk figure is indicated in its panel.

Figure 5 illustrates the reuse of a goal-conditioned field across multiple initial states: the field represents a map and goal, while path extraction starts from each selected location.

Fig. 5: Goal-conditioned field prediction and reuse across initial states.

Figure 6 brings together walk geometry, label variability, and spatial sampling cost.

Fig. 6: WoSt walk and label statistics: (a) walk-length variability and boundary treatment; (b) label and gradient variability on a fixed map, with differences measured against a finite Monte Carlo reference; (c) spatial variation in mean walk length.

Fig. 7: Local search on one fixed urban-map field: (A) full map and routes; (B,C) enlarged neighborhoods for local search (4) and (10). Greater depth exposes endpoints farther along the passage. Gray cells are obstacles; overlays mark candidate regions and routes. Table 1 reports aggregate controller success.

Planning under Noise with Neural Operators

Abstract

References

Meow summary

The problem

The boundary-condition trick

Learn a field, then search

Results and limits

Make notes of your own

Abstract

1 Introduction

2 Preliminaries

2.1 Stochastic Formulation of Path Planning

Control problem.

Path-planning specialization.

Boundary conditions.

2.2 Neural Operator

3 Related Works

Path Planning.

Monte Carlo Methods.

Neural Operator.

4 Planning under Noise (PUN) with Neural Operator

4.1 The Walk-on-Stars Chain

Geometry and stopping.

The estimator.

4.2 Label Construction

Bootstrapping.

4.3 Operator Parameterization and Objective

4.4 Grid-Based Local Search

5 Experiments

Backbone and amortization.

Supervision.

Walk horizon and learned continuation.

Training objective.

Path extraction.

Transfer across discretizations.

Replanning under execution noise.

6 Discussion

Extensions.

Limitations.

7 Conclusion

AI use statement

Acknowledgements

References

Appendix A Notation Table

Appendix B Screened Poisson Kernels in ddd Dimensions

B.1 Free Space

Radial derivative.

B.2 The Ball: Dirichlet, Neumann, and Robin

B.3 The Screened Poisson Kernel and the Step Weight

B.4 The Screened Ball Poisson Kernel and the Star-Step Weight

B.5 Why a Ball Kernel Is Admissible on a Star Domain

B.6 Dimension-Specific Forms and the Laplace Limit

Appendix C Proofs

C.1 Standing Regularity Assumptions

C.2 Exact Linearization

C.3 Positivity and Bounded Log-Desirability

C.4 Why We Use Reflection at Obstacles

Zero Dirichlet data produces a logarithmic singularity.

Reflection preserves a bounded target.

C.5 The Chain Kernel

C.6 Log-Space Variant

Bias correction.

C.7 Delta-Method Bias of the Log Label

C.8 Bounded Local Search

Reachable endpoint scores.

Recovering a route.

Saturation at large depth.

Appendix D Supplementary Experiments

D.1 Evaluation and Training Protocol

Evaluation scope.

Execution model for replanning comparisons.

Baseline configuration.

Planning objective and geometry.

MPPI implementation.

PI2^{2}2-CMA implementation.

Hardware and timing.

Statistical reporting for replanning.

Implementation of the WoSt estimator.

Shell error.

Backbone compute policy.

Supervision parameterization.

Local-search implementation.

D.2 Walk Horizon and Label Budget

Label budget and planning success.

D.3 Resolution Generalization

Appendix E Benchmarks

2D (zero-shot; 592 maps).

3D (zero-shot; 278 maps).