HJ-Gauss: A Monte-Carlo HJ Reachability Scheme

Backward reachable tubes (BRTs), computed via grid-based levelset methods for viscous Hamilton-Jacobi (HJ) PDEs, provide principled safety certificates for learned controllers and planning algorithms in control and learning-enabled systems. However, classical grid-based HJ solvers require $O(M^n)$ memory footprint for $M$ grid points per $n$ state dimension. This renders them impractical for high-dimensional systems. We address this bottleneck with a local PDE linearization that enables a frozen-coefficient sampling scheme for the viscous HJ PDE: a generalized Cole-Hopf-type transformation reduces the nonlinear HJ equation to a sequence of linear heat equations, which admits Gaussian heat-kernel representations via the Feynman-Kac formula. The value function and its spatial gradient are then recovered via roll-outs of Monte Carlo expectations on Gaussian densities, yielding a storage-free and grid-free algorithm that scales as $N\cdot n$ for $N$ samples. This decoupling of memory from dimensionality enables reachability analysis on large-scale problems: safety analysis on European starlings’ ($\mathit{\text{sturnus}}\mathit{\text{vulgaris}}$) emergent behavior validated on $\mathbf{100,000}$ simulated starlings motion – modeled as 4D aerial Dubins vehicles. We prove a finite-sample concentration bound $O(N^{-1\mathrm{/}2})$ error, conditional linear convergence rates, and establish robustness properties for our introduced scheme. Numerical validation on pursuit-evasion games against the grid-based levelset method demonstrates relative $L_{\text{rel}}^2$ errors of $0.03-0.20$, with $14-26$ second wall-clock times per 2D slice on a CPU; and with validation on $n=45$-dimensional multi-agent 2D rocket games. Our numerical results demonstrate real scalability of HJ reachability safety verification on large scale multi-agent systems.