BLUF (Bottom Line Up Front) - (1/2)

Grid-resolved Hamilton-Jacobi (HJ) BRTs yield certified safety sets at a $O(M^n)$ memory cost.

• Background: LevelSetPy resolves backward reachable sets/tubes by discretizing the HJ PDE (upwinding + Lax-Friedrichs + TVD-RK) on a grid — certified, but exponential in state dimension $n$.

• Idea: a Cole-Hopf transformation turns the viscous HJ PDE into a linear heat equation; Feynman-Kac evaluates it as a Gaussian expectation. No grid.

HJ-Gauss replaces the grid with Gaussian Monte-Carlo, dropping the footprint to $O(N\cdot n)$.

BLUF (Bottom Line Up Front) - (2/2)

HJ-Gauss replaces the grid with Gaussian Monte-Carlo, dropping the footprint to $O(N\cdot n)$.

• Generality: a frozen-coefficient Picard iteration handles the nonconvex, state-dependent Hamiltonians of reachability.

• Rigor: $O(N^{-1/2})$ concentration, contraction convergence, a linearization residual bound, and a conservative safety certificate.

• For this room: the windowed reachable tube is a dynamics-aware conflict predicate for MAPF; the value gradient is a score for diffusion planners.

This Gaussian kernel is the diffusion kernel, the bridge to learning-based planners.

Roadmap — Main Talk (Parts 1-5)

Continue: Parts 6-9 on the next slide.

Roadmap — Main Talk (Parts 6-9) & Appendix

Part	Theme
6	Application: multi-agent path finding (AMFS)
7	The diffusion connection
8	Dirty Laundry: limits & boundaries
9	Conclusions, future work, references
A–C	Appendix: HJ/viscosity theory (A) · reachability foundations (B) · LevelSetPy grid pipeline (C)

Background — notation plus HJ PDE & viscosity theory (Part A), reachability foundations (Part B), and the LevelSetPy grid pipeline (Part C) — is collected in the Appendix after the closing slide.

The Curse of Dimensionality, Quantified

Grid memory is $O(M^n)$ with $M$ points per dimension:

State dim $n$	$M=100$ cells	Feasibility (double array)
2	$10^{4}$	trivial
3	$10^{6}$	easy
4	$10^{8}$	~0.8 GB, heavy
5	$10^{10}$	~80 GB
6	$10^{12}$	~8 TB
45 (our multi-agent game)	$10^{90}$	exceeds atoms in the observable universe

• GPU acceleration lowers wall-clock, not the exponent.

Reachability's most valuable regime — many interacting agents, rich dynamics — is where grids are impossible.

Attempt 1: PINN Reachability (DeepReach)

• DeepReach trains a neural network to minimize the HJ PDE residual directly (physics-informed loss), no reference grid.

• **Pro:**Trades grid storage for network weights; sidesteps $O(M^n)$ memory.

• Con: Scales only to moderate dimensiond (reported ~9-10D);

  • And accuracy degrades as $n$ grows; training is a nonconvex optimization with no viscosity-solution guarantee.

Catch-22: Learning-based solvers help, but do not deliver a certified, dimension-robust representation.
HJ-Gauss keeps the viscosity-solution semantics.

Attempt 2: Convex-Duality Formulas (Hopf / Lax-Oleinik)

HJ-Gauss trades their exactness for generality, at the quantified price of a linearization residual bound (Part 4).

• Hopf and Lax-Oleinik formulas¹ evaluate the value pointwise via convex optimization — grid-free at a point.

• Limitation: they require a convex (or state-independent) Hamiltonian.

• Reachability's $H(x,p)=\max_u\min_w\langle p,f\rangle$ is state-dependent and nonconvex → outside this class.

It handles the nonconvex, state-dependent Hamiltonian head-on.

¹ Darbon-Osher 2016; Chow-Darbon-Osher-Yin 2017; Kirchner 2018

Attempts 3-4: Decomposition, Path Integral, Sampling

The open gap: a grid-free, storage-light scheme for nonconvex, state-dependent, adversarial Hamiltonians. That is HJ-Gauss.

¹ Chen-Herbert · ² Kappen; Theodorou · ³ Summers-Lygeros; Lesser-Oishi

Where HJ-Gauss Sits

• Keeps the viscosity-solution semantics (unlike PINN heuristics).

• Handles nonconvex, state-dependent $H$ (unlike Hopf/Lax-Oleinik).

• Indifferent to coupling structure (complements decomposition).

• Certifies adversarial worst-case safety (unlike probabilistic sampling).

• Memory $O(N\cdot n)$, grid- and discretization-free.

The rest is how: a Cole-Hopf transformation, a Feynman-Kac expectation, and a frozen-coefficient Picard iteration.

The Theoretical Roadmap of the Method

1. Start from the viscous HJ PDE (Part A).
2. Apply a generalized Cole-Hopf transform $\omega=e^{-cv}$ → a Linear heat equation.
3. Solve the heat equation by its Gaussian heat kernel → a Feynman-Kac expectation.
4. Recover value (log-sum-exp) and gradient from that expectation.
5. Because $c$ depends on the unknown $Dv$, iterate: freeze, solve, update (Picard) = Algorithm 1.
6. Analyze exactness (quadratic case) vs the quasi-linearization residual (general case).

Each arrow is a theorem or proposition in the paper; we walk them one by one.

The Transformation, and the Closed-Form/Approximate Split

For a general $H$, define the spatially-varying coefficient

$$c(t;x)=\frac{2}{\delta}\cdot\frac{H^\delta}{|Dv^\delta|^2},\qquad H^\delta:=H(t;x,Dv^\delta),$$

and set $\omega^\delta:=\exp(-c\,v^\delta)$.

• Quadratic case $H=\tfrac12|p|^2$: $c=1/\delta$ constant, $\omega^\delta$ solves the homogeneous heat equation with zero residual — a genuine Cole-Hopf identity.

• General $H$: the transform induces a residual $R=R_{\mathrm{alg}}+R_{\mathrm{der}}$;

  • $c$'s choice eliminates the algebraic part $R_{\mathrm{alg}}$, leaving the derivative part $R_{\mathrm{der}}$.

No freezing of $c$'s values removes its derivatives. Hence for nonlinear $H$ this is a quasi-linearization, not an identity — and the leftover $R_{\mathrm{der}}$ is the residual we will bound.

Picard Quasi-Linearization: The Surrogate View

• At iteration $k$, freeze $c^{(k)}$ at the current iterate and solve the linear heat equation in closed form.

• Equivalently: replace the true Hamiltonian by the locally-matched quadratic surrogate

$$\tilde H^{(k)}(p)=\tfrac{\delta}{2}\,c^{(k)}\,|p|^2,$$

solve the resulting viscous HJ equation in closed form, obtain $v^{(k+1)}$, update $Dv^{(k+1)}$ and $c^{(k+1)}$.

• The iteration's limit is the fixed point of the surrogate solve map — equal to the viscous solution only when $R_{\mathrm{der}}\equiv 0$ (the quadratic case).

"Successive quadratic matching": each step is a locally closed-form quadratic-Hamiltonian solve, refined by re-estimating the co-state.

Reduction to the Heat Equation (Proposition)

With $c^{(k)}$ frozen, $\omega^\delta$ solves the heat initial-value problem

$$\omega_t^\delta-\tfrac{\delta}{2}\Delta\omega^\delta=0\ \text{ in }\mathbb{R}^n\times(0,T],\qquad \omega^\delta(0,x)=\exp(-c\,g(x)).$$

• When $H=\tfrac12|p|^2,\ c=1/\delta$: closed-form, no residual.

• For general $H$: discards the derivative residual (Lemma, appendix) — the quasi-linearization.

A nonlinear, nonconvex, state-dependent HJ PDE has become the heat equation — the most classical linear parabolic PDE, with an explicit Gaussian solution.

The Gaussian Heat Kernel & Feynman-Kac

The heat IVP has the unique bounded Green's-convolution solution

$$\omega^\delta(t,x)=\frac{1}{(\sqrt{2\pi\delta t})^{n}}\int_{\mathbb{R}^n} e^{-\frac{|x-y|^2}{2\delta t}}\,e^{-c\,g(y)}\,dy,$$

which, by Feynman-Kac, is a Gaussian expectation:

$$\omega^\delta(t,x)=\mathbb{E}_{y\sim\mathcal N(x,\,\delta t I_n)}\big[e^{-c\,g(y)}\big].$$

• The free-space kernel integrates to one over $\mathbb{R}^n$ (and no proper subset) — every estimator is posed on all of $\mathbb{R}^n$.

• Sampling $y\sim\mathcal N(x,\delta t I_n)$ and averaging $e^{-cg(y)}$ is an unbiased estimator of $\omega^\delta$.

The solution is a Gaussian roll-out average of the terminal cost. This is the entire computational payload — no grid, no marching.

Recovering the Value (Lemma)

Invert the transform:

$$v^\delta(t,x)=-\frac{1}{c^{(k)}}\log\,\mathbb{E}_{y\sim\mathcal N(x,\delta t I_n)}\big[e^{-c^{(k)} g(y)}\big].$$

• A log-sum-exp of Gaussian-sampled terminal costs.

• Numerically evaluated in a log-domain-stable form (log-mean-exp) to avoid overflow/underflow.

• $M$ arbitrary evaluation states $x_1,\dots,x_M$ — laid on a grid or scattered — each with its own fresh $N$ samples.

The value at a query point is a scalar reduction over $N$ samples. Embarrassingly parallel across both $M$ states and $N$ samples.

Recovering the Gradient (Corollary)

The spatial gradient (co-state) has a closed Gaussian-expectation form:

$$Dv^\delta=\frac{1}{t\,\delta\,c^{(k)}}\left(x-\frac{\mathbb{E}_{y\sim\mathcal N(x,\delta t I_n)}\!\big[y\,e^{-c^{(k)} g(y)}\big]}{\mathbb{E}_{y\sim\mathcal N(x,\delta t I_n)}\!\big[e^{-c^{(k)} g(y)}\big]}\right).$$

• A weighted-mean shift: the gradient points from $x$ toward the exponentially-weighted centroid of the samples.

• No finite differencing, no stored field — the co-state is a byproduct of the same samples.

This gradient is (up to constants) a score function $\nabla\log\omega$. Hold that: it is the entire bridge to diffusion in Part 7.

Admissible Data & Unbiasedness (Recap as a Guarantee)

• Datum $\omega^\delta(0,\cdot)=e^{-cg}$ is continuous, strictly positive, bounded above by $e^{-c g_{\min}}$, decaying at infinity.

• ⇒ the heat solution is the unique bounded one; the integral converges absolutely; the Gaussian expectation is genuine and unbiased.

• Only $g\ge g_{\min}$ is required (bounded target); no upper bound on $g$.

The estimator drawing $y\sim\mathcal N(x,\delta t I_n)$ over all of $\mathbb{R}^n$ is unbiased for the quantity inside the logarithm — the foundation for the concentration bound of Part 4.

Sign and the Removable Zero of the Coefficient

• $c\sim C_H/(\delta|Dv^\delta|)$: it grows without bound in flat regions and vanishes and flips sign where optimal dynamics run tangent to the level set (the barrier).

• The zero of $c$ is a removable singularity of the estimator: as $c\to 0$,

$$-\tfrac{1}{c}\log\mathbb{E}[e^{-c g}]\ \to\ \mathbb{E}[g],$$

i.e. the log-sum-exp gracefully degenerates to the Gaussian mean of the terminal data — the correct pure-diffusion limit ($H\to0$ reduces the PDE to the heat equation).

• In implementation: evaluate in a log1p-stable form and clip $|c|$ to $[c_{\min},c_{\max}]$.

The coefficient's pathologies are analytically understood and numerically tamed — clipping is itself a bounded, priced perturbation (next slide).

Coefficient Regularization (Lemma)

Run Algorithm 1 with the regularized coefficient

$$c_\eta(t;x)=\frac{2}{\delta}\cdot\frac{H^\delta}{|Dv^\delta|^2+\eta},$$

equivalent to running it unregularized for the perturbed Hamiltonian $H_\eta=H\,|Dv^\delta|^2/(|Dv^\delta|^2+\eta)$, which satisfies a uniform bound

$$\|H_\eta-H\|_\infty\le \frac{C_H\sqrt\eta}{2},$$

so the induced value-function perturbation is at most $T C_H\sqrt\eta/2$.

Regularization cures the flat-interior blow-up at a quantified, uniform price — no hidden instability.

Algorithm 1: Quasi-Linearization (Cole-Hopf)

$$\text{Fix } \varepsilon>0,\qquad v^{(0)}(t,x)=g(x),\qquad c^{(0)}=\dfrac{2\,H(t,x,Dg)}{\delta\,\lvert Dg\rvert^{2}}.$$

$$\begin{aligned} &\textbf{for }\ k=0,1,2,\dots\text{:}\\[1pt] &\quad\text{(1) freeze }c^{(k)}\text{ at the current iterate;}\\[1pt] &\quad\text{(2) solve }\ \omega_t=\tfrac{\delta}{2}\Delta\omega,\ \ \omega^{(k)}(0,x)=e^{-c^{(k)}g(x)};\\[1pt] &\qquad\quad\Rightarrow\ \omega^{(k+1)}(t,x)=\mathbb{E}_{y\sim\mathcal N(x,\delta t I)}\!\big[e^{-c^{(k)}g(y)}\big]\ \ (\text{Monte-Carlo});\\[1pt] &\quad\text{(3) recover }\ v^{(k+1)}=-\tfrac{1}{c^{(k)}}\log\omega^{(k+1)};\\[1pt] &\quad\text{(4) update }Dv^{(k+1)},\ \ c^{(k+1)}=\dfrac{2\,H(t,x,Dv^{(k+1)})}{\delta\,\lvert Dv^{(k+1)}\rvert^{2}}\ \ \text{clip to }[c_{\min},c_{\max}];\\[1pt] &\quad\text{(5) stop when }\ \lVert v^{(k+1)}-v^{(k)}\rVert/\lVert v^{(k)}\rVert<\varepsilon. \end{aligned}$$

• Steps (2)-(4) are all Gaussian-expectation reductions over $N$ samples.

• Convergence typically in ≤ 20 iterations (empirically 12-15).

Every line is grid-free. The only state carried between iterations is the coefficient field $c^{(k)}$ at the $M$ evaluation states.

Evaluation States vs Monte-Carlo Samples (Remark)

• $N$ = i.i.d. Gaussian draws $y_1,\dots,y_N\sim\mathcal N(x,\delta t I_n)$ used independently at each evaluation state $x$ to form the estimators.

• Resampled fresh at every evaluation state and every Picard iteration.

• $M$ = the evaluation states $x_1,\dots,x_M$ where $v$ is queried — carry no randomness, may be a uniform grid or scattered arbitrarily.

• Total randomness per iteration: $M\times N$ draws.

• The zero level set is recovered by evaluating $v^\delta$ at the $M$ states and isocontouring at level zero (marching cubes) — a deterministic post-process consuming no further samples.

The solve is grid-free; the grids that appear in results are evaluation windows for visualization and error measurement, not sampling constraints.

Memory: $O(N\cdot n)$ vs $O(M^n)$

• Example: $N=14{,}000,\ n=3$ → ≈ 0.6 MB/iteration; $n=45$ → ≈ 7.2 MB/iteration (vs $10^{90}$ grid cells).

Memory now scales with sample budget and dimension, not with a grid raised to the dimension. This is the whole contribution in one line.

The Weight-Degeneracy Problem

• The estimators are a ratio of exponential-weight expectations, e.g. the gradient's numerator/denominator both carry $e^{-cg}$.

• In high dimensions (or small $\delta$), $c\,g$ is large over the kernel support, so the denominator $\mathbb{E}[e^{-cg}]$ is carried by a rare event.

• The weights collapse onto a few samples → the ratio's variance explodes; the effective sample size crashes.

Naive Gaussian sampling is unbiased but can be catastrophically high-variance right where we need it (high $n$, small $\delta$). We fix it with importance sampling.

Tilted Proposal & Density-Ratio Reweighting

• Replace $y\sim\mathcal N(x,\sigma^2 I_n)$ ($\sigma^2=\delta t$) by a shifted proposal $q_\theta=\mathcal N(x+\theta,\sigma^2 I_n)$.

• Reweight by the density ratio

$$w_\theta(y)=\exp\!\Big(\tfrac{|\theta|^2-2\langle y-x,\theta\rangle}{2\sigma^2}\Big),$$

which leaves every expectation unbiased.

We move the samples to where the integrand has mass, then correct via the ratio. Bias is untouched; variance can collapse by orders of magnitude.

The Zero-Variance Shift Is One Gradient Step

• The zero-variance proposal is proportional to the integrand $\varphi_x(y)\,e^{-cg(y)}$ itself.

• Its first-order Laplace / Gaussian approximation gives the shift

$$\theta^\star=-\sigma^2\,c^{(k)}\,Dv^{(k)}(t,x).$$

• That is, samples are pushed one preconditioned gradient step along the descent direction of the running value iterate.

The previous Picard iterate hands the sampler its drift for free. Reachability and importance sampling become the same computation — the value gradient is the optimal sampling direction.

Effective Sample Size as a Self-Diagnostic

• Monitor the effective sample size of the self-normalized weights $\tilde w_i=w_\theta(y_i)e^{-cg(y_i)}$:

$$\mathrm{ESS}=\frac{(\sum_i \tilde w_i)^2}{\sum_i \tilde w_i^2}\in[1,N].$$

• $\mathrm{ESS}\approx N$: weights well-spread → estimate certified.

• $\mathrm{ESS}\to 1$: weight degeneracy → flag, increase $N$, raise $\delta$, or re-tilt.

ESS is a free, per-query health check: it tells you when to trust the estimate and why it failed when it does.

Four Error Sources, One Budget

The distance from the computed field $\hat v^{K,N,\delta}$ to the true inviscid viscosity solution $v$ decomposes into four independently-controlled sources:

Three of four shrink with computational effort; the residual is a property of the surrogate and sets the floor. We now bound each.

Finite-Sample Concentration (Theorem)

Fix phase $(t,x)$ and frozen $c>0$; assume $g_{\min}\le g\le g_{\max}$ on the sampling support. With $Z=e^{-cg(\zeta)}$, $\mu=\mathbb E[Z]$, $v_c=-\tfrac1c\log\mu$, $\alpha=e^{-c g_{\max}}$, $\beta=e^{-c g_{\min}}$, and $\hat v_{c,N}=-\tfrac1c\log\bar Z_N$:

$$\mathbb P\!\left(|\hat v_{c,N}-v_c|\ge\varepsilon\right)\ \le\ 2\exp\!\left(-\frac{2N\mu^2(1-e^{-c\varepsilon})^2}{(\beta-\alpha)^2}\right).$$

• Bounded weights $\alpha\le Z_i\le\beta$ ⇒ Hoeffding.

• Rate is the standard $O(N^{-1/2})$, independent of state dimension $n$ — the crux of the scalability claim.

Dimension enters the cost per sample, not the number of samples for a target accuracy. That is how we beat $O(M^n)$.

Explicit Sample Size, and Its Worst Case (Corollary + Remark)

• Since $\mu\ge\alpha$, it suffices to take

$$N\ \ge\ \frac{(\beta-\alpha)^2}{2\alpha^2(1-e^{-c\varepsilon})^2}\log\frac{2}{\alpha}$$

to guarantee error $\ge\varepsilon$ with probability $\le\alpha$.

• Worst-case caveat: with $c=1/\delta$, $\beta/\alpha=e^{(g_{\max}-g_{\min})/\delta}$, so this bound grows like $e^{2(g_{\max}-g_{\min})/\delta}$ — exponential in $1/\delta$, while viscosity error is only $O(\sqrt\delta)$.

The accuracy gain from shrinking $\delta$ and the sample cost pull in opposite directions. This Hoeffding bound is loose; the next slide tightens it, and experiments sit far below it.

Tightening: Bernstein & Jensen

• Bernstein: when $\mathrm{Var}(Z)\ll(\beta-\alpha)^2/4$, a variance-aware tail bound gives strictly tighter concentration than Hoeffding.

• Jensen: the loose lower bound $\mu\ge\alpha=e^{-c g_{\max}}$ can be replaced by $\mu\ge e^{-c\,\mathbb E[g(\zeta)]}$, substantially reducing the required $N$.

• Importance sampling (Part 3) attacks the same variance directly.

The pessimistic $e^{1/\delta}$ is a worst-case artifact of Hoeffding + crude $\mu$ bound. In practice $N\in[14\text{k},20\text{k}]$ suffices at $\delta\in[0.08,0.1]$.

The Iteration Operator (Contraction Setup)

• Fix evaluation states $x_1,\dots,x_M$ and time $t$; let $\mathcal G$ be a stable gradient-reconstruction operator, $p_m(v)=(\mathcal G v)_m$.

• Coefficient-update map: $(\Gamma(v))_m=2H(t,x_m,p_m(v))/(\delta|p_m(v)|^2)$.

• Frozen-coefficient heat-kernel map:

$$(\Phi(c))_m=-\tfrac{1}{c_m}\log\mathbb E_{\zeta\sim\mathcal N(x_m,\delta t I_n)}\big[e^{-c_m g(\zeta)}\big].$$

• Algorithm 1 is the iteration $v^{(k+1)}=\Lambda(v^{(k)})$, $\Lambda=\Phi\circ\Gamma$.

The whole method is a fixed-point iteration of a composed map on $\mathbb R^M$ with the sup-norm. Contraction ⇒ convergence.

Assumptions for Contraction (Assumption A)

On a closed admissible set $\mathcal A$ (sup-norm), assume:

1. $|g|\le G$;
2. non-degenerate gradient & bounded coefficient: $m_0\le|p_m|\le P_*$, $c_{\min}\le(\Gamma v)_m\le c_{\max}$;
3. Hamiltonian Lipschitz in co-state: $|H|\le H_*$, $|H(\cdot,p)-H(\cdot,q)|\le L_H|p-q|$;
4. reconstruction Lipschitz: $\|\mathcal G v-\mathcal G w\|_\infty\le L_D\|v-w\|_\infty$;
5. invariance + contraction constant

$$q=\frac{2G}{c_{\min}}\cdot\frac{2L_D}{\delta}\Big(\frac{L_H}{m_0^2}+\frac{2H_*P_*}{m_0^4}\Big)<1.$$

Non-degeneracy (item 2) holds on the tube band $\{|v^\delta|\le\eta\}$ where the gradient is bounded away from zero — where the barrier lives. Outside, use the regularized coefficient.

Contraction Convergence (Theorem)

$\Lambda$ is a contraction on $\mathcal A$ ⇒ unique fixed point $v^\star$, and for any $v^{(0)}\in\mathcal A$:

$$\|v^{(k+1)}-v^\star\|_\infty\le q\,\|v^{(k)}-v^\star\|_\infty\le q^{\,k+1}\|v^{(0)}-v^\star\|_\infty,$$

with the a posteriori estimate

$$\|v^{(k)}-v^\star\|_\infty\le \frac{q}{1-q}\,\|v^{(k)}-v^{(k-1)}\|_\infty.$$

• Iteration count for tolerance $\varepsilon$: $K=O(\log(1/\varepsilon))$ — geometric convergence.

The residual you watch (successive-iterate change) bounds the error you cannot see (distance to $v^\star$). That is the a posteriori estimate — a runtime stopping certificate.

Scope of the Contraction (Remark)

• The theorem is about the frozen-coefficient numerical map, not the original nonlinear PDE.

• The fixed point $v^\star$ solves the viscous HJ equation with $H$ replaced by its quadratic surrogate at $v^\star$.

• $v^\star=v^\delta$ iff $H=\tfrac12|p|^2$ (then $c\equiv1/\delta$, $R_{\mathrm{der}}\equiv0$).

• For general $H$, the two differ by the quasi-linearization residual $\mathcal E_{\mathrm{ql}}=\|v^\star-v^\delta\|_\infty$.

We do not claim unconditional global convergence to the true HJ solution. We claim geometric convergence to a surrogate fixed point, plus a bound on the gap — next.

The Quasi-Linearization Residual (Theorem, via Duhamel)

With converged coefficient $c^\star$, coefficient-variation bound $L_c$, value bound $G_v$:

$$\mathcal E_{\mathrm{ql}}=\|v^\star-v^\delta\|_\infty\ \le\ \frac{T}{c_{\min}}\,e^{2c_{\max}G_v}\big(\bar R_{\mathrm{alg}}+\bar R_{\mathrm{der}}\big),$$

• $\bar R_{\mathrm{alg}}\le c_{\max}(L_H+2H_*P_*/m_0^2)\,\|Dv^\delta-Dv^\star\|_\infty$ — contracts as the surrogate gradient tracks the true one.

• $\bar R_{\mathrm{der}}\le L_c G_v(1+\tfrac\delta2)+\tfrac\delta2 L_c^2 G_v^2+\delta(1+c_{\max}G_v)L_c P_*$ — governed by coefficient variation $L_c$.

Proof idea: $\tilde\omega=e^{-c^\star v^\delta}$ solves the heat equation with a source $\omega B$ (Lemma); Duhamel + heat-semigroup $L^\infty$-contraction bounds $\|\tilde\omega-\omega^\star\|$; mean-value theorem for $\log$ converts back to $v$.

The bound vanishes in the quadratic case ($L_c=0$, $Dv^\star=Dv^\delta$) and is largest where $H/|Dv|^2$ turns over sharply — the usable-part boundary.

Reading the Residual

• Vanishes when $H=\tfrac12|p|^2$: recovers the closed-form Cole-Hopf.

• Governed by $L_c$, the variation rate of the converged coefficient — largest near the barrier where the Hamiltonian-to-gradient ratio shocks.

• Algebraic part $\propto\|Dv^\delta-Dv^\star\|$: the residual contracts as the Picard iteration refines the co-state.

• Certifies smallness where $c^\star$ is slowly varying; offers no reprieve where it is not.

This is why experimental errors concentrate at the zero-level-set boundary — the theory predicts the location of the worst error.

The Total Error Bound (Theorem) & Bias-Variance

Chaining all four sources by the triangle inequality:

$$\|\hat v^{K,N,\delta}-v\|_\infty\ \le\ C_1\rho^K+\frac{C_2}{\sqrt N}\sqrt{\log(1/\delta_p)}+\mathcal E_{\mathrm{ql}}+C_3\sqrt\delta.$$

• Bias-variance tradeoff in $\delta$: viscosity bias $\downarrow$ as $\delta\downarrow$, but MC variance $\uparrow$.

• Balancing the two nontrivial $\delta$-dependent terms yields the optimal

$$\delta\sim N^{-1/3}\ \Rightarrow\ \text{overall rate}\ O(N^{-1/6}),$$

slower than plain MC $O(N^{-1/2})$ but dimension-robust and scalable.

A single decomposition tells you how to spend budget: iterate to kill $\rho^K$, sample to kill $N^{-1/2}$, choose $\delta\sim N^{-1/3}$, and accept the residual floor.

The Conservative Safety Certificate (Corollary)

Let $E$ be the total error budget. Declare state $x$:

• SAFE only if $\hat v^{K,N,\delta}(x)>E$;

• UNSAFE only if $\hat v^{K,N,\delta}(x)<-E$;

• UNDETERMINED otherwise ($|\hat v|\le E$).

With probability $\ge 1-\delta_p$: no unsafe state is certified safe and no safe state is certified unsafe; all classification error is confined to a declared band of width $2E$ about the boundary.

• The certificate errs on refusal, never on admission.

• This is the accuracy quotation: worst-case sign-correctness with an explicit abstention band, not an $L^2$ average that flatters the interior.

Robustness to Model Uncertainty

• The scheme is stable under small perturbations of the Hamiltonian and terminal cost (robustness theorems, appendix): perturb $H$ by $\Delta H$ in sup-norm ⇒ value perturbed by $\le T\|\Delta H\|_\infty$ (Duhamel/comparison).

• Consequence: perfectly known dynamics and terminal costs are not required — real-world models with bounded error still yield certified-with-margin sets.

• The coefficient-clipping and regularization perturbations are of this priced form.

The method degrades gracefully and quantifiably under model error — a prerequisite for deployment on real robots with imperfect models.

Experimental Setup & Statistics BLUF - (1/2)

• Hardware: single CPU (Intel i7-14700K, 20 cores, 31 GiB RAM, Ubuntu 22.04); JAX on CPU backend

  — Consistent with the memory-frugality claim.

• Replication: Each experiment re-solved with 30 independent Monte-Carlo seeds;

  • Evaluation points and the LevelSetPy reference held fixed - only sampler randomness varies. Tables report mean ± 1 s.d.

Figures show one representative seed for clarity; tables report full 30-seed statistics. Multiplicity-controlled reporting throughout.

Experimental Setup & Statistics BLUF - (2/2)

• Significance: Three families of Holm-Bonferroni-corrected tests at $\alpha=0.05$:

  • (a) paired Wilcoxon of 30-seed-averaged MC field vs grid reference;

  • (b) one-sided $t$-test that per-seed $L^2_{\mathrm{rel}}$ lies below the Crandall-Lions bound $\sqrt\delta$;

  • (c) cross-condition Mann-Whitney (Rockets vs Dubins; speed regimes).

Figures show one representative seed for clarity; tables report full 30-seed statistics. Multiplicity-controlled reporting throughout.

Benchmark 1: Two-Vehicle Pursuit-Evasion Games

• Rockets launch game¹: two identical rockets — pursuer $P$, evader $E$ — on the $x$-$z$ plane, thrusts $(u_p,u_e)$, relative orientation $\theta=u_p-u_e$. Capture when $\|PE\|<r$.

• Dubins two-car game (Merz): relative $(x_1,x_2,\theta)$, symmetric turn-rate bound.

• Both are 3D relative-state games; the BRT is a 3D tube, visualized as 2D $(x,z)$ slices at fixed $\theta$.

• Target: $\ell_2$-ball of capture radius $r=1.5$; horizon $(0,1]$; $\delta=0.08$; $N=14{,}000$ (rockets) / $20{,}000$ (Dubins) samples per iteration.

These are the canonical HJ-Isaacs benchmarks — and the Dubins game is the pairwise MAPF conflict of Part 6.

¹ Dreyfus 1966

The Two-Rockets Problem (Schematic)

Relative geometry: evader $E$ fixed at the origin (accel. $a_e$, control $u_e$, gravity $g$); pursuer $P$ at relative orientation $\theta=u_p-u_e$ on the $(x,z)$ plane. — HJ-Gauss (Molu et al., 2026).

The Dubins Air3D Problem (Schematic)

Air3D relative coordinates with the evader fixed at the origin: relative position $(x_1,x_2)$ and relative heading $\psi$; turn-rate-bounded pursuer vs evader. — LevelSetPy (Molu, ACM TOMS 2025 / IEEE CDC 2024).

Rockets Pursuit-Evasion BRT

$(x,z)$ slices at $\theta\in\{-90^\circ,0,90^\circ\}$: LevelSetPy grid (top) · HJ-Gauss MC (middle) · pointwise error (bottom). — HJ-Gauss (Molu et al., 2026).

Dubins Two-Car BRT

LevelSetPy grid (top) · HJ-Gauss MC (middle) · error (bottom). — HJ-Gauss (Molu et al., 2026); LevelSetPy (Molu, TOMS 2025 / CDC 2024).

Benchmark 1: Quantitative Results

• Reference: LevelSetPy $45^3$ grid interpolated to a $40\times40$ evaluation grid.

• Holm-Bonferroni one-sided $t$-test rejects $L^2_{\mathrm{rel}}\ge\sqrt\delta=0.283$ at $p_{\mathrm{holm}}<10^{-50}$ in every condition.

Sampling-plus-iteration error is statistically well inside the viscosity budget — not merely numerically.

Interpreting the Errors: Different Currencies

• The Crandall-Lions $O(\sqrt\delta)\approx0.283$ bounds the inviscid-vs-viscous sup-distance; the table measures MC-vs-grid in $L^2/L^\infty$ — different currencies, compared with care.

• Averaging does not equal agreement: a paired Wilcoxon of the 30-seed-averaged MC field vs the grid rejects equality at $p_{\mathrm{holm}}<10^{-8}$ everywhere — a systematic quasi-linearization residual remains (Theorem, Part 4).

• $\theta=0$ Dubins reaches $L^2_{\mathrm{rel}}=0.024$ (smooth interior); $\theta=\pm\pi/2$ and rockets sit ~$0.09$-$0.13$ (share of grid near the coefficient-turnover boundary).

• $L^\infty\sim0.7$-$1.4$ driven by the zero-level-set boundary where $|Dv^\delta|$ is maximal — where the residual $L_c$ term predicts.

The report is Corollary (conservative certificate): sign-correctness outside an abstention band; $L^2$ flatters the interior, $L^\infty$ is dominated by the declared-undetermined band.

A Physical Check: The Gravity Asymmetry

• Rockets are asymmetric between $\theta=-90^\circ$ and $\theta=+90^\circ$: the gravity term $(g-a-a\sin\theta)$ gives a $g-2a$ drift at $+90^\circ$ but $g$ at $-90^\circ$.

• Holm-Bonferroni test over 30 seeds detects this asymmetry for Rockets ($p_{\mathrm{holm}}<10^{-9}$) but not for the gravity-free Dubins vehicle ($p_{\mathrm{holm}}=0.53$).

The method resolves a genuine physical effect, and the statistics correctly find it present where physics says it should be and absent where it should not. Validation, not just fitting.

Benchmark 2: A 45-Dimensional Pursuit-Evasion Game

• 15 rockets (14 pursuers, 1 evader), state $x\in\mathbb R^{45}$: each agent $(x_i,y_i,\theta_i)$, control $u_i\in[-1,1]$.

• Dynamics $\dot x_i=a_i\cos\theta_i,\ \dot y_i=a_i\sin\theta_i,\ \dot\theta_i=u_i$; target $\phi(x)=\min_{i\le14}\|x_i^{\mathrm{pos}}-x_{15}^{\mathrm{pos}}\|-r_{\mathrm{capture}}$.

• Three regimes: evader faster ($2.0$ vs $1.0$), balanced ($1.0$/$1.0$), pursuers faster ($1.0$ vs $2.0$).

• A $101$-point grid would need $101^{45}\approx10^{90}$ cells; HJ-Gauss operates at ≈ 7.2 MB per iteration.

This is the flagship: a reachability computation grid solvers cannot even represent, done on a single CPU in seconds.

Benchmark 2: 45D Results

• Memory constant at 7.2 MB/iteration; residuals below $10^{-2}$ (indeed ~$10^{-4}$) → near-zero Picard residual floor.

• Min-max Hamiltonian evaluates in closed form (box-constrained control ⇒ sign structure), so each iteration is $O(N\cdot n)$ arithmetic and memory.

• Pairwise Holm-Bonferroni Mann-Whitney: no significant residual difference between speed regimes ($p_{\mathrm{holm}}\ge0.93$) → floor set by $(N,\delta)$, not game parameters.

At $n=45$ there is no grid reference; the report is measurable quantities — iteration stability, memory, wall-clock — i.e. scalability, not certified accuracy.

Why the 45D Residual Floor Is Lower Than 3D

• 3D benchmarks floor at $0.02$-$0.06$; the 45D game floors at $0.0002$-$0.0006$ — two orders of magnitude lower.

• Reason: the residual floor tracks proximity to the coefficient-turnover region (the shock-prone barrier), not dimension per se.

• 45D evaluation states are drawn uniformly over $[-100,100]$ per coordinate — with small probability of landing near the boundary — so the average residual is small.

A subtle but important message: error concentrates geometrically (at the barrier), not dimensionally. High dimension is not inherently high error.

Cost Accounting: Time vs Memory

• 3D isosurface over $15{,}625$ points: 129.4 s (rockets, single CPU) vs 1.3 s for the dense $45^3=91{,}125$-point LevelSetPy grid.

• Per-iteration memory: ≈ 0.6 MB (MC, $N=14$k, $n=3$) vs ≈ 1.5 MB (grid value+gradient).

• In low dimensions ($n\le4$) grids win — smaller constant, no sampling variance.

• The Monte-Carlo overhead becomes favorable where grid storage becomes prohibitive ($n\ge5$).

The two paradigms are complementary, not competing. Use grids where they fit in memory; use HJ-Gauss where they cannot.

Benchmark 3: Safety at Population Scale ($10^5$ Starlings)

• Certify 100,000 European starlings (Sturnus vulgaris) as 4D aerial Dubins vehicles under predator attack.

• Structure (explicit): the murmuration is partitioned into flocks; each flock has its own value function; the population safe set is the aggregation of per-flock zero sublevel sets — many coupled low-D games, not one joint high-D PDE.

• Because the solve is grid-free, the value-function solve is independent of bird count; the population enters only through parallel per-bird certification — the operational content of $O(N\cdot n)$ at $10^5$ agents.

Two curses dodged at once: dimensionality (grid-free solve) and agent cardinality (per-flock decoupling + parallel evaluation).

The Certified Population Snapshot

2,000 of 100,000 birds (colored by heading) with 7 flock centers (stars) and their capture discs (dashed). Black contour is the certified safe-set boundary at $\tau=0$: an annular cordon around a protected core. — HJ-Gauss (Molu et al., 2026).

Watch the Safe Set Change Shape

• Backward time $\tau$ runs outward from the capture set: the certified safe set of the attacked flock grows and changes topology.

• Green crosses and dashed circles are the seven attacking predators with their capture discs.

• Blue interior = certified safe; the black curve is the $v=0$ boundary.

A safety certificate for $10^5$ agents, as one evolving surface.

Cordon, Then Collapse

Same solve, six instants. A hole persists in the safe set ($\chi=0,\ \beta_1=1$, cordon) through $\tau=1.30$; by $\tau=1.39$ it has closed ($\chi=1,\ \beta_1=0$, cohesion). — HJ-Gauss (Molu et al., 2026).

Topology as a Safety Instrument

The certificates recover the field-documented collective repertoire as topological events of the reachable set, compressed to three integers per time step $(\chi,\beta_1,n_{\mathrm{comp}})$:

The triple $(\chi,\beta_1,n_{\mathrm{comp}})$ tells an operator not merely that safety is being lost but how — nucleation vs cordon vs split. Safety posture as a 3-integer summary.

The Three Integers, Over Time

• $\chi$ steps $0\to1$ while $\beta_1$ drops $1\to0$ at the same $\tau=1.39$: one event, two independent witnesses.

• $n_c$ holds at $1$ throughout, so this flock never fragments.

• Vacuole nucleation is the reverse event: a predator penetration attaches a 1-handle and drops $\chi$ by one.

Three integers per tick say whether safety is being lost, and how.

The Repertoire Is Field-Documented

• Real starling responses to predation: cordon, tube, funnel, and split.

• Each one appears in our certificates as a topological signature of the safe set, and none of them is a hand-coded rule.

• The certificate reproduces a documented repertoire from dynamics and the game alone — which is the reason to trust it on a floor, where the repertoire is congestion, yielding, and deadlock instead.

Fragmentation event, $n_c:1\to2$. Field murmuration imagery as reproduced in HJ-Gauss (Molu et al., 2026).

The Floor, as MAPF Sees It

Storage blocks separated by single-cell travel corridors, task endpoints on the border columns — the RHCR fulfillment-warehouse map family (Li et al., AAAI 2021). CBS/ECBS resolve the two classical conflict types marked above.

What Counts as a Conflict

• Vertex conflict $\langle a_i,a_j,x,t\rangle$: two agents occupy the same cell at the same timestep.

• Edge (swap) conflict $\langle a_i,a_j,x\!\to\!y,t\rangle$: two agents traverse the same edge in opposition — the head-on case in a one-cell corridor.

• CBS branches on a conflict, adds a constraint to one agent, and re-plans that agent's shortest path; ECBS bounds suboptimality; RHCR applies this inside a rolling window.

• Every one of these predicates is a statement about cells and timesteps.

Nothing in this vocabulary mentions turn radius, actuation lag, or localization error.

CCBS Fixes Time, Not Dynamics

CCBS (Andreychuk et al., IJCAI 2019) replaces integer timesteps with continuous time: agents are discs moving along straight edges with arbitrary durations, and conflicts become unsafe time intervals. Turn radius, actuation lag, and disturbance remain outside the model.

The Gap in Classical MAPF

• Multi-Agent Path Finding (MAPF): move a team of agents to goals without collision. Solvers: CBS/ECBS (conflict-based search), PIBT, LaCAM/LaCAM*, and lifelong RHCR.

• All operate on a discrete grid with simplified kinematics (unit moves / rotate-in-place) and no disturbance model.

• A "conflict" is a shared vertex/edge (or swept-disc overlap) — dynamics-blind.

• Consequence: realized collisions on a real floor (turn-radius limits, localization noise, actuation lag) are out of scope.

The discrete plan is provably conflict-free on the graph — yet the continuous robot can still collide. That is the gap HJ reachability fills.

The Scaling Anchor: RHCR¹

• Rolling-Horizon Collision Resolution: decompose lifelong MAPF into a sequence of Windowed MAPF instances; resolve collisions only within a bounded time window $w$, ignore beyond.

• Scales to 1,000 agents on warehouse maps; co-authored by J. Durham (Amazon Robotics).

• Current frontier: LaCAM* (10,000+ agents, near-optimal, anytime), PIBT (thousands in ms), learning-based lifelong (imitation, guidance-graph, MAPF-GPT).

• All still discrete-grid, deterministic.

HJ-Gauss does not compete with LaCAM's combinatorial search; it certifies the plan under continuous dynamics — a different, complementary axis.

¹ Li et al., AAAI 2021

The Mapping: HJ-Gauss Constructs → MAPF Roles

The Dubins pursuit-evasion BRT from Part 5 is the canonical pairwise MAPF conflict, computed rigorously.

The Dynamics-Aware Conflict Set (Pairwise Dubins BRT)

Black = certified inevitable-collision boundary; dashed = naive disc. The BRT bulges beyond the disc in a heading-dependent way — what a dynamics-blind check cannot see.

What the Certified Set Looks Like

• Sweeping the relative heading $x_3$ makes the certified set rotate and deform.

• A near head-on geometry grows a long lobe; a near-parallel one stays compact.

• The dashed red circle is the dynamics-blind disc, fixed for every heading.

This is the object a swept-disc test approximates with a circle.

Reading the Heading Sweep

• Only about 2% of the relative-state box is certified inevitable-collision at $w=0.6$ s.

• So the gate is selective, not blanket-conservative: it fires on geometry that actually traps the pair.

• Where the tube exceeds the disc, a geometric planner is blind; where it falls inside, the disc is needlessly cautious.

Falsifiable Hypotheses (Stated to Be Proven Wrong)

Naming: HJ = Hamilton-Jacobi (the method); HB-n = hypothesis $n$.

The AMFS Pipeline: World and Planner

• World: an RHCR-style structured warehouse grid — storage blocks, single-cell corridors, border endpoints — with lifelong task reassignment.

• Planner: windowed prioritized planning (space-time A*), the RHCR rolling horizon, shared by both policies.

• Executor: a continuous unicycle/Dubins rollout with bounded turn rate and bounded disturbance, where realized collisions are detected in continuous space.

The planner is held fixed on purpose. Only the safety layer varies.

The AMFS Pipeline: The Two Policies

• Geometric — the classical dynamics-blind vertex/edge conflict test, unshielded.

• HJ — the same planner plus the windowed-BRT runtime shield.

• Both run under Common Random Numbers: identical spawns, tasks, priorities, and disturbance realizations for a given seed.

• BRT precompute: the JAX HJ sampler (Algorithm B) on the Dubins-relative Hamiltonian, cached once; the online loop is numpy-only $O(1)$ lookups.

One clean contrast: the only difference between the two arms is whether the certificate is actuated.

The Simulated Floor

RHCR-style structured floor: storage blocks (obstacles), border endpoints (pickup/dropoff), 14 lifelong agents (circles) heading to goals (stars).

Watch It Run: Same Seed, Same Noise, One Difference

Left: geometric conflict check. Right: the HJ-Gauss windowed-BRT shield. Amber ring = shield firing (agent fleeing), red = realized contact, stars = current goals. Seed 1028 of the 30-seed harness: 98 → 66 collisions at identical throughput (0.825). Animated in the HTML deck; the PDF shows the opening frame.

The Same Episode, as a Filmstrip

Four instants from the run above: top row geometric, bottom row HJ-shield. Counters under each panel accumulate realized collisions and shield activations. Identical spawn, task, and disturbance streams throughout (CRN).

The Key Scientific Finding (Why the Naive Version Fails)

• Planning-only BRT predicate made things worse (HJ 84 vs geometric 53 collisions in the first trial). Two principled reasons:
  1. the pursuit-evasion BRT is a worst-case adversarial set, but agents are cooperative — they never apply the BRT's safe control, so plan-time membership does not predict cooperative-tracking collisions;
  2. plan-time checks see only discrete nodes at integer times; collisions happen in continuous time under disturbance.
• The certificate only bites when actuated. Fix: use the windowed BRT as a runtime shield — the lower-priority agent in an in-BRT pair flees (a braking "spin-in-place" bug had it get hit).

This is the arc: a negative result diagnosed to a principle, then the correct mechanism. Actuation, not inspection.

The Multiseed Statistical Harness

Mirrors the WIP-forecast harness; numpy-only:

• ≥30 seeds (configurable);

• MSER-5 warm-up truncation of the per-tick collision series (+ Welch running-mean cross-check);

• Bootstrap 95% CIs on every headline metric;

• Common Random Numbers (CRN) across policies (spawn/task/disturbance streams keyed to the seed) for a paired comparison;

• Paired bootstrap CI for $\mathrm{mean}(\text{hj}-\text{geom})$ + Holm-Bonferroni step-down FWER control across the HB family.

Single-seed results were noisy and even inverted; the harness is what turns a direction into a defensible claim.

Results: 30 CRN-Paired Seeds (14 Agents)

• ~32% fewer realized collisions, significant under Holm-Bonferroni, at no throughput/flowtime cost in this model.

• Aggregated sweeps: 18-37% reduction across densities and disturbance levels.

HB-1/HB-2 rejected in the safe direction; HB-3/HB-4 show no throughput penalty. The certificate, actuated, works.

Headline: ~32% Fewer Realized Collisions

30 CRN-paired seeds, bootstrap 95% CI. Geometric 97.7 → HJ-shield 66.1; paired diff −31.6 [−38.0, −25.4], Holm p≈0.

Per-Tick Collisions & MSER-5 Warm-Up

HJ shield (teal) sits below geometric (red) at every steady-state tick; the dashed line is the MSER-5 warm-up cut discarded before averaging.

The Effect Holds as the Fleet Grows

Collision reduction persists across agent densities (8-20 agents) — the certificate's value does not wash out with congestion.

Verdict Table

ID	Prediction	Result	Verdict
HB-1	Shield is sound	HJ never certifies a colliding plan safe; collisions ↓
HB-2	Dynamics-aware advantage	31.6 fewer collisions per run, paired, Holm p≈0
HB-3	$O(N\cdot n)$ memory	7.2 MB at $n=45$; flat in resolution	(paper)
HB-4	Little throughput cost	Throughput/flowtime diff = 0 in this model	(with caveat)
HB-5	Fast enough	BRT precompute ~10 s; $O(1)$ online lookup
HB-6	Diffusion-guidance helps	Deferred to Option C