Scalable Fixed-Point Framework for High-Dimensional Hamilton-Jacobi Equations

Technical Summary: Scalable Fixed-Point Framework for High-Dimensional Hamilton-Jacobi Equations

Problem Statement
Hamilton-Jacobi (HJ) partial differential equations are fundamental to optimal control, mechanics, and dynamic systems. However, solving these equations in high-dimensional state spaces presents a severe challenge known as the "curse of dimensionality." Traditional grid-based methods (e.g., WENO, semi-Lagrangian) require discretizing the state space, leading to exponential growth in computational cost and memory as dimensions increase, rendering them infeasible for problems exceeding ten dimensions. Conversely, recent deep learning approaches, while scalable, often lack rigorous convergence guarantees, require extensive training, and incur high computational costs due to repeated derivative calculations via automatic differentiation.

Methodology
The authors propose a novel, mesh-free, and gradient-free fixed-point framework for computing viscosity solutions of high-dimensional HJ equations. The core of the method relies on the Hopf-Lax formula, which characterizes the solution $u(x, t)$ as the infimum of a variational problem:
$$u(x, t) = \inf_{y \in \mathbb{R}^n} \left\{ t H^*\left(\frac{x - y}{t}\right) + g(y) \right\}$$
where $H^*$ is the Legendre transform of the Hamiltonian $H$, and $g$ is the initial condition.

Instead of solving the variational problem via global optimization or grid discretization, the authors reformulate the condition for the minimizer $y^*$ as a fixed-point equation:
$$y = x - t \nabla H(\nabla g(y))$$
This relation motivates a Picard iteration scheme:
$$y^{(k+1)} = x - t \nabla H(\nabla g(y^{(k)}))$$
Key features of the methodology include:

• Mesh-Free and Causality-Free: The method evaluates the solution at any point $(x, t)$ independently, without relying on neighboring grid points or time-stepping from $t=0$.
• Derivative-Free (in the context of the solver): While the formula involves gradients of $H$ and $g$, the iteration itself does not require differentiating the solution $u$ or using backpropagation through a neural network.
• Multiple Initialization Strategy: To handle non-smooth solutions where characteristics intersect (creating "kinks") and the fixed point is no longer unique, the authors propose sampling multiple random initializations. The algorithm converges to all local fixed points, and the one yielding the minimum Hopf-Lax energy is selected as the viscosity solution. Theoretical analysis (Theorem 2.4) establishes that this approach is probabilistically complete, ensuring the global minimizer is found with probability approaching one as the number of initializations increases.

Key Contributions and Theoretical Guarantees

1. Convergence and Error Bounds: Under assumptions of strict convexity for $H$ and Lipschitz continuity for $\nabla H$ and $\nabla g$, the authors prove (Theorem 2.1) that the fixed-point operator is a contraction mapping, guaranteeing linear convergence to a unique global minimizer.
2. Quantitative Error Estimates: Theorem 2.2 provides explicit bounds relating the residual of the iteration to the errors in the minimizer, the solution value, and the solution gradient. This ensures that controlling the iteration residual controls both the value function and the associated optimal control.
3. Iteration Complexity: Theorem 2.3 demonstrates that the number of iterations required to reach a prescribed accuracy scales as $O(\log(1/\epsilon))$, largely independent of the problem dimension $d$.
4. Control Recovery: The method directly yields the optimal control $q^* = (x - y^*)/t$ without additional post-processing or characteristic tracing.

Experimental Results
Numerical experiments were conducted in dimensions ranging from 1 to 100, comparing the proposed method against WENO (grid-based) and ImplicitHJ (deep learning-based).

• Accuracy: The fixed-point method consistently achieved errors near machine precision ($10^{-12}$ to $10^{-17}$) across all dimensions. In contrast, WENO errors grew significantly with dimension, and WENO became computationally infeasible beyond 3D. ImplicitHJ showed increasing errors in high dimensions (e.g., $L_2$ error of $0.137$ in 100D) and required massive training times (hundreds of thousands of seconds).
• Efficiency: The proposed method exhibited near dimension-independent computational time (e.g., 0.05s for 1D vs. ~0.05s for 100D in the Burgers' test) and minimal memory usage (0.001 MB).
• Non-Smooth Solutions: The method successfully captured sharp features (kinks) in solutions arising from non-smooth initial data (e.g., $g(x) = x|x|$) and moving kinks, provided the multiple initialization strategy was employed.
• Control Problems: The framework was validated on Linear Quadratic Regulator (LQR) problems and problems with cubic Hamiltonians, maintaining high accuracy and low computational cost up to 100 dimensions.

Significance and Claims
The paper claims to offer a scalable, provably accurate, and computationally efficient alternative for high-dimensional HJ equations. By avoiding grids, characteristic tracing, and the heavy training overhead of neural networks, the method achieves a computational cost that is largely independent of dimensionality. The authors emphasize that their approach complements existing classical and deep learning schemes by providing a theoretically sound framework with guaranteed convergence and error bounds, making it particularly suitable for real-world applications in optimal control, robotics, and finance where both accuracy and tractability in high dimensions are critical. The method is presented as a direct local evaluation tool that eliminates the approximation errors associated with spatial discretization and neural network training.