Low Residuals, Wrong Value: Learning the Viscosity Solution of Hamilton-Jacobi-Isaacs Equations with PINNs
Authors: Raul Jimenez, Pablo Tejerina-Pérez, Leonid Sarieddine, Pedro Tarancón-Álvarez, Marisol Traforetti, Pau Solé-Vilaró, Oleksii Ignatenko, Nazar Tkach, Pavlos Protopapas
Organizations: Institut de Ciències del Cosmos (ICCUB), Universitat de Barcelona, Barcelona, Spain · Institució Catalana de Recerca i Estudis Avançats (ICREA), Barcelona, Spain · Departament de Física Quàntica i Astrofísica, Universitat de Barcelona, Barcelona, Spain · Faculty of Applied Sciences, Ukrainian Catholic University, Lviv, Ukraine · Department of Computer Science, National University of Kyiv-Mohyla Academy, Kyiv, Ukraine · Institute for Applied Computational Science, Harvard University, Cambridge, MA, USA
Physics-informed neural networks (PINNs) and other residual-trained solvers judge a candidate solution by the size of its sampled PDE residual. For first-order Hamilton-Jacobi equations this is not sufficient: the same equation and boundary data admit Lipschitz almost-everywhere solutions that differ from the viscosity solution by O(1) while violating the equation only on sets of measure zero. We study whether residual minimization selects such solutions in a fixed-horizon pursuit-evasion game whose Hamilton-Jacobi-Isaacs value function is known in closed form on the circle S1 and on the flat tori T2 and T4. We train four neural parameterizations with ten seeds on each geometry, 120 runs in total. A standard network with periodic Fourier features converges to folded almost-everywhere solutions in all 30 runs: small residuals and the expected capture-layer gradient, yet values that are wrong over most of the domain. A network with periodicity imposed through the loss recovers the viscosity solution in one and two dimensions but becomes unreliable in four. Two parameterizations that receive the geodesic distance do substantially better: a single network fed the wrapped displacement and its distance recovers the viscosity solution in all 30 runs, while a minimum over periodic image charts succeeds in 27 of 30 and introduces a wrong-image failure mode. Finally, for the convex Hamiltonian considered here, we use the semiconcavity of the viscosity solution to build a mesh-free finite-difference test that needs only forward network evaluations; it detects the folded solutions and the harmful min-of-charts failures without access to the exact value function. A small sampled residual thus does not by itself validate a learned Hamilton-Jacobi solution; the parameterization and independent structural diagnostics can be essential for selecting the viscosity solution.
Figures & tables
Figure 1: The value function V∗ of Eq. ( 7 ) for ε=0.08 , Δv=0.2 . Left: V∗ on T2 in the (d,τ) plane. We show the capture set (black band), the cone dc=ε+Δvτ (green line), the capture crease d=ε (a convex kink on the boundary), the horizon at which the cone reaches the point of T2 farthest from the pursuer (dotted line at τ=6.67 ), the nearest points of the cut locus (dashed d=1 ) and the distances spanned by the cut locus (hatched band). Right: A time-slice V∗(Δ,τ=4) on the fundamental domain of T2 . The dashed, red circle marks the cone dc=ε+Δvτ ; inside it the value function grows linearly with d , while outside it is V∗=τ .
Figure 2: The chart selector m(Δ,τ) on T2 , drawn on [0,2)2 with the pursuer at the origin (0,0) and its periodic-images at the other corners. The faces defining the cut locus of the pursuer on the torus, which are Δ~i=±1 on the fundamental domain, appear as the midlines Δi=1 . Left. Analytic nearest-image index m∗(Δ,τ) . The four colors index each of the charts. Arrows show the optimal direction of motion of the pursuer, from the nearest periodic image, towards the evader. The charts meet on the ridge (cut locus). In the fundamental domain [−1,1)2 (coordinates of the text), the nearest image is always in the chart with bm=0 . In the representative coordinates of the figure Δ∈[0,2)2 all signs are positive, so zm=Δ−Lbm , the images sit at Lbm , and the four colors are bm={(0,0),(1,0),(0,1),(1,1)} . Middle. Selector m^(Δ,τ=4) learned by the min-of-charts NN (seed 0). Horizon-limited states, where all charts tie, are shown in grey. On the distance-limited states, the learned m^(Δ) agrees with the analytic value (colors in left plot) everywhere. Right. The nearest-image distance d=minm∥zm∥ , equal to ∥Δ∥ in the fundamental domain. It is continuous, with a concave kink at the ridge, where the index m∗ jumps, and ∥∇d∥=1 on both sides.
Configuration
Best validation loss
RMS error against V∗
∥∇V∥ at capture
Within 5% of V∗
Circle S1
Fourier-feature
1.5[0.8,4.6]×10−4
1.36[0.93,1.86]
5.00[4.99,5.01]
44[19,60]%
soft-periodic
1.5[1.2,1.9]×10−5
8.7[7.4,9.5]×10−3
5.00[5.00,5.01]
100%
min-of-charts
2.2[1.2,9.0]×10−6
4.4[3.6,6.8]×10−3
5.00[5.00,5.00]
100%
geometric
2.1[1.6,2.6]×10−6
4.3[3.9,4.8]×10−3
5.00[5.00,5.00]
100%
Torus T2
Table 1: Metrics for the four configurations on the three geometries Tn ( n=1,2,4 ), ten seeds each, with the same training budget within a geometry. ( 5×105 iterations on S1 , 106 on T2 and T4 ). Entries are the median over seeds with the seed range in brackets (in the units of the entry). Best validation loss is Eq. ( 11 ) on a fixed validation set; for the soft-periodic configuration it includes the extra term of Eq. ( 15 ). The RMS is computed against the analytic V∗ on 2×105 uniform states ( τ∼U(0.05,8) , d>ε ); ∥∇V∥ at capture is the mean gradient norm in the band close to the capture set, d∈(ε,1.5ε] , τ∼U(0.5,7.5) , whose analytic value is 5 ; the last column is the share of the same 2×105 states within 5% relative error of V∗ (with the caveat in the item 5 at the beginning of the section). Bold marks the best median within each geometry (lowest loss and error, highest share within 5% ), with ties all marked.
Figure 3: Validation loss on S1 , T2 and T4 (left to right) for ten seeds per configuration, computed from Eq. ( 11 ) evaluated on a fixed random sampling. The Fourier-feature NN (red) ends 70 – 180 times above the min-of-charts (orange) and geometric (blue) NNs, which descend together. The soft-periodic NN (green), whose loss includes the additional term ( 15 ), ends between the two groups on S1 and T2 and above the Fourier-feature NN on T4 . All losses are still decreasing at the end of the budget.
Figure 4: Value function (left) and absolute error (right) of every configuration, against the analytic V∗ (gray), evaluated at τ=3 of an example seed (seed 0). We show S1 (top), a cut along the diagonal of T2 (middle) and a cut along the main diagonal of T4 (bottom). The dotted line is the cone dc=ε+Δvτ , beyond which V∗=τ .
Figure 5: Value function V averaged over directions: at each geodesic distance d , V is evaluated on all directions of the sphere ∥Δ∥=d inside the fundamental domain ( 2 on S1 , 720 on T2 , 512 on T4 ), and the mean and standard deviation are computed. The solid vertical line marks d=1 , beyond which only part of the sphere lies inside the domain. The dotted vertical line marks the cone, dc . Same run and panel disposition as Fig. 4 . Left: mean over directions with a ±1σ band. The analytic V∗ depends on d alone (not on direction), so any spread in the 1σ band is error due to angular variation of the learned solutions. Right : mean absolute error with its ±1σ band.
Figure 6: Structure maps on T2 at τ=4.7 (seed 0): from left to right columns, learned value VNN , gradient norm ∥∇ΔV∥ , pointwise residual ∣r∣ from Eq. ( 10 ) (log scale) and absolute error ∣VNN−V∗∣ , as a function of the wrapped coordinate Δ . In each row, we show the four configurations. The dashed red circle is the cone dc . The Fourier-feature NN has a gradient ∥∇ΔV∥ of correct magnitude 1/Δv almost everywhere, but a lattice of folds instead of a cone in the value V , with the residual concentrated on the fold lines and an error that grows from the capture set outward reaching O(1) values. The soft-periodic NN has the correct shape of the value V , with maximum absolute error of 10−1 concentrated around the cone. The min-of-charts and geometric NNs are within 3⋅10−2 of V∗ everywhere.
Figure 7: Equivalent structure maps as in Fig. 6 , but for T4 , on the slice Δ3=Δ4=0 at τ=4.7 (seed 0). The Fourier-feature NN shows the same fold lattice as Fig. 6 , with again O(1) absolute errors over most of the slice, increasing as one moves away from the capture set. The soft-periodic NN has the right overall shape, but errors between 0.5 and 1, largest around the cone, especially along the line Δ1=0 and, more weakly, along part of the line Δ2=0 . The min-of-charts and geometric NNs are within about 5⋅10−2 of V∗ on the whole slice, with the largest errors along the cone, and much lower values inside. In particular, notice that the first two configurations lack to recover the angular symmetry.
Figure 8: Evolution of Fourier-feature NN into almost-everywhere, folded solutions. The value function V(Δ,τ=3) of the Fourier-feature NN on S1 (seed 0) is shown as a red line at different stages of a 105 –iteration run, from the initialization. The analytic V∗ is shown as a gray line. The gray band is the capture set ∥Δ∥≤ε .
Appendix figures & tables2 assets
Supplementary material from the paper’s appendix.
Appendix
candidate ( τ=4 )
value RMS
flagged %
max ratio λ^max/C∗
circle S1
V∗ (analytic)
—
0.0
0.0
planted sawtooth
—
8.2
100
Fourier-feature
1.36[0.93,1.86]
11.6[5.2,18.8]
212[88,568]
soft-periodic
8.7[7.4,9.5]×10−3
0.0[0.0,0.0]
0.8[0.4,1.1]
min-of-charts
4.4[3.6,6.8]×10−3
0.0[0.0,3.6]
1.8[0.4,5.9]
Appendix
Table 2: Monte Carlo semiconcavity test at τ=4 ( 4×104 probes; eight random directions each and one on S1 ; h=0.02 ; the region d≤3ε and concave kinks are masked; a probe is flagged when λ^max>2C∗ ), applied to V∗ , to a planted sawtooth (a τ -independent triangle wave in the first coordinate Δx , of period 0.5 and amplitude 1.25 , so with slope of magnitude 1/Δv on its flanks and convex kinks at its valleys, which shows what a clearly folded function scores) and to all 120 trained networks; median and range over the ten seeds. The value RMS column repeats Table 1 ( τ∈(0.05,8) , not only τ=4 ). On S1 the curvature of V∗ vanishes away from its kinks, so its ratio is 0 rather than 1 . Bold marks the lowest median among the trained networks within each geometry, with ties all marked.
Figure 9: Seed census of the 120 runs; one marker per run, colour by configuration, marker shape by topology. Left. Value RMS against V∗ (Table 1 ) grouped by topology. Right. The semiconcavity certificate maxλ^max/C∗ (Table 2 ) against the same RMS; the dashed line is the ratio 10 . Folded and relabelled runs lie above it and correct runs below it, while the T4 soft-periodic failures, smooth creases rather than folds, straddle it.
Many core problems in nonlinear systems analysis and control can be recast as solving partial differential equations (PDEs) such as Lyapunov and Hamilton-Jacobi-Bellman (HJB) equations. Physics-informed neural networks (PINNs) have emerged as a promising mesh-free approach for approximating their solutions, but in most existing works there is no rigorous guarantee that a small PDE residual implies a small solution error. This paper develops verifiable error bounds for approximate solutions of Lyapunov and HJB equations, with particular emphasis on PINN-based approximations. For both the Lyapunov and HJB PDEs, we show that a verifiable residual bound yields relative error bounds with respect to the true solutions as well as computable a posteriori estimates in terms of the approximate solutions. For the HJB equation, this also yields certified upper and lower bounds on the optimal value function on compact sublevel sets and quantifies the optimality gap of the induced feedback policy. We further show that one-sided residual bounds already imply that the approximation itself defines a valid Lyapunov or control Lyapunov function. We illustrate the results with numerical examples.
Jun Liu
Department of Applied Mathematics, Faculty of Mathematics, University of Waterloo, Waterloo, Ontario N2L 3G1, Canada
Physics-Informed Neural Networks (PINNs) are a common class of machine learning-based partial differential equation (PDE) solvers which train a network to represent a solution by minimizing a residual loss that encodes the PDE. Despite their successes, they are known to fail on certain simple equations, converging to an incorrect solution despite low loss. These failure modes have garnered significant attention in the literature over the past several years, motivating both architectural and optimization based solutions. By directly visualizing the residual, we show that failure modes are the result of overfitting: the loss is minimized on the collocation points, but not elsewhere. Applying regularization causes the failure modes to vanish. Finally, we extend double backpropagation over the full set of residuals, and use it to achieve state-of-the-art performance on four standard failure mode equations with up to 23× fewer collocation points and a vanilla architecture.
Nigel T. Andersen, Takashi Matsubara
Graduate School of Information Science and Technology Hokkaido University Sapporo, Japan · RIKEN Center for Advanced Intelligence Project (AIP) Tokyo, Japan
Classical finite-difference solvers remain reliable tools for partial differential equations, but their efficiency depends on where mesh resolution is placed. Uniform refinement can waste degrees of freedom when solution difficulty is localised near sharp gradients, fronts, oscillations, or constraint-sensitive regions. This paper studies a hybrid strategy in which a physics-informed neural network (PINN) is used not as the final solver, but as an off-grid residual probe for adaptive mesh refinement. The PINN residual is sampled over the domain, converted into cellwise indicators, and used to guide refinement before the final approximation is computed by a finite-difference solver. The method is evaluated on three benchmarks. The main full-solver validation uses the one-dimensional viscous Burgers equation with a nonuniform finite-difference solve on the adapted meshes. PINN-threshold refinement attains final relative L2 error 0.021067 with 60 degrees of freedom, compared with 0.022617 for uniform refinement with 192 degrees of freedom. At matched mesh size, PINN-threshold reduces the error by about 67.5%. PINN-Dörfler refinement gives similar performance, with error 0.021264 using 58 degrees of freedom. A gradient indicator remains slightly more accurate, so the result supports usefulness rather than universal superiority. Manufactured 2D and 3D proxy tests, based on a nonlinear Schrödinger equation and an incompressible Navier--Stokes system, show that PINN residuals can organise structured refinement and improve over random refinement, although they do not consistently outperform gradient or uniform baselines. The results support PINN-guided AMR as a residual-indicator strategy for transferring physics-informed diagnostic information into finite-difference mesh adaptation while preserving the classical solver as the final approximation engine.
Henry Kasumba, Ronald Katende
Department of Mathematics · Makerere University · 7062, Kampala, Uganda +2