Pith. sign in

REVIEW 2 major objections 3 minor 34 references

Hessian-augmented Supervised Learning for Hamilton-Jacobi-Bellman PDEs

T0 review · 2 major / 3 minor · reviewed 2026-06-26 · grok-4.3

Challenge this review Re-run · record.json Download PDF Read on arXiv ↗ Review JSON

Pith's one-line read Hessian data from Riccati equations along PMP trajectories cuts the samples needed for accurate value-function approximation by up to an order of magnitude.

desk verdict Hessian data from the Riccati equation along PMP trajectories cuts the samples needed for accurate hyperbolic-cross polynomial approximation of HJB value functions by up to 10x on the control-affine examples. read the letter →

arxiv 2606.23827 v1 pith:FQ5N3FGZ submitted 2026-06-22 math.OC cs.LGcs.NAmath.NA

Matías Gómez-Aedo , Behzad Azmi , Yuyang Huang , Dante Kalise , Karl Kunisch This is my paper · ORCID
classification math.OCcs.LGcs.NAmath.NA
keywords Hamilton-Jacobi-BellmanequationsoptimalcontrolPontryaginMaximumPrincipleRiccatiequationsupervisedlearningpolynomialregressionvaluefunctionapproximationfeedbacksynthesis
verification ladder T0 review T1 audit T2 compute T3 formal

The pith

A machine-rendered reading of the paper's core claim, the machinery that carries it, and where it could break.

The reading

The paper presents a supervised learning method that generates training data for Hamilton-Jacobi-Bellman value functions by solving the Pontryagin Maximum Principle system from many initial states. Each trajectory supplies not only the value and gradient but also the Hessian, computed via an auxiliary matrix Riccati equation. These second-order quantities are inserted as additional linear equations into a weighted least-squares fit over sparse polynomial bases on hyperbolic cross index sets. The resulting polynomial yields an analytic feedback law. Numerical tests on control-affine systems show that the extra derivative information produces lower approximation error and better closed-loop performance than value-only or gradient-only regression at the same sample count.

What carries the argument

Augmentation of the weighted least-squares system with gradient and Hessian rows derived from PMP trajectories and the associated Riccati equation, performed over sparse polynomial bases on hyperbolic cross index sets.

What would settle it

On a fixed high-dimensional test problem, run the regression with and without the Hessian equations using identical numbers of trajectories and measure whether the closed-loop cost and L2 error of the value function differ by more than the observed variance across random initial-condition sets.

Watch

Extended reading notes

Core claim

Incorporating Hessian information obtained from the matrix Riccati equation along optimal trajectories into a polynomial regression reduces the number of required training samples by up to an order of magnitude while improving both pointwise accuracy of the learned value function and the performance of the recovered feedback law.

Load-bearing premise

The matrix Riccati equation solved along PMP trajectories supplies accurate second derivatives of the true value function.

Editorial extensions

If this is right

  • Feedback laws are recovered by direct differentiation of the learned polynomial without additional optimization.
  • A partial Hessian strategy keeps data-generation cost manageable when the state dimension is large.
  • Second-order augmentation yields measurable gains in both approximation accuracy and closed-loop performance on problems of increasing dimension.
  • The method applies to deterministic control-affine dynamics and produces globally defined polynomial approximations.

Reading between the lines

Editorial extensions of the paper, not claims the author makes directly.

  • The same derivative-augmentation idea could be tested on other basis families that admit cheap evaluation of gradients and Hessians.
  • If the Riccati-derived Hessians remain reliable under moderate noise, the approach might extend to nearby stochastic problems where a similar linearization exists.
  • The reduction in sample count suggests that second-order information could be leveraged in other regression-based HJB solvers that currently rely only on value or gradient data.
Share X Bluesky LinkedIn Reddit HN

Signed reviews

No signed human review yet.

Editorial analysis

A structured set of objections, weighed in public.

Desk editor's note, referee report, simulated authors' rebuttal, and a circularity audit.

Referee Report

2 major / 3 minor

Summary. The manuscript develops a data-driven method for approximating value functions of deterministic optimal control problems with nonlinear control-affine dynamics. It solves the Pontryagin Maximum Principle optimality system from multiple initial conditions to generate training data consisting of values, gradients, and Hessians of the value function (Hessians obtained via the matrix Riccati equation along optimal trajectories). These quantities augment a weighted least-squares regression over sparse polynomial bases on hyperbolic cross index sets, with gradients and Hessians supplying additional linear equations per sample. Feedback laws are recovered analytically from the learned value function. In high dimensions a partial Hessian strategy is used to control data-generation cost. The approach is validated on control-affine problems of increasing state dimension, where second-order augmentation is reported to improve accuracy and closed-loop performance with up to an order-of-magnitude reduction in required training samples relative to value-only regression.

Significance. If the reported empirical gains hold under broader testing, the method provides a concrete way to reduce sample complexity when learning HJB value functions by systematically incorporating second-order information generated from standard PMP+Riccati theory. This is potentially useful for moderate-to-high-dimensional deterministic control problems where trajectory data is expensive to obtain. The construction is parameter-free in the sense that no new constants are fitted to define the result, and the use of hyperbolic-cross polynomial bases together with the analytic feedback recovery are clear strengths.

major comments (2)
  1. [§5] §5 (numerical experiments): the reported order-of-magnitude sample-complexity reduction is shown only relative to value-only regression on the same polynomial basis; no comparison is supplied against other established HJB approximation techniques (e.g., deep neural network solvers or other sparse-grid methods), which limits assessment of whether the second-order augmentation yields a genuine advance beyond the immediate baseline.
  2. [§4.2] §4.2 (regression formulation): the weighting factors that balance the value, gradient, and Hessian equations in the least-squares objective are introduced without an explicit rule or sensitivity study; because the central claim rests on the benefit of the additional linear equations, the lack of justification for these weights is load-bearing for reproducibility of the reported gains.
minor comments (3)
  1. [Abstract] The abstract states quantitative gains but the main text should cross-reference the specific tables or figures that contain the error metrics and sample counts supporting those claims.
  2. [§3] Notation for the hyperbolic-cross index sets and the partial-Hessian selection mask could be introduced earlier (e.g., in §3) rather than only in the experimental section to improve readability.
  3. [Figures in §5] Figure captions should explicitly state the number of Monte-Carlo repetitions or random seeds used for each regression run so that variability of the reported accuracy improvements can be assessed.

Simulated Author's Rebuttal

2 responses · 0 unresolved

We thank the referee for the careful reading, positive assessment of the contribution, and constructive comments. We address the two major comments point by point below.

read point-by-point responses

  1. Referee: [§5] §5 (numerical experiments): the reported order-of-magnitude sample-complexity reduction is shown only relative to value-only regression on the same polynomial basis; no comparison is supplied against other established HJB approximation techniques (e.g., deep neural network solvers or other sparse-grid methods), which limits assessment of whether the second-order augmentation yields a genuine advance beyond the immediate baseline.

    Authors: We agree that the experiments isolate the benefit of Hessian augmentation relative to value-only regression on the identical polynomial basis and hyperbolic-cross index set. This design choice was deliberate to highlight the effect of the additional linear equations supplied by gradients and Hessians. Broader benchmarking against deep neural network HJB solvers or other sparse-grid techniques would indeed provide a more complete picture of the method's standing, but would have required substantial extra implementation and compute that lay outside the scope of the present work. In the revision we will add an explicit limitations paragraph in §5 acknowledging this and identifying such comparisons as a natural direction for follow-up research. revision: partial

  2. Referee: [§4.2] §4.2 (regression formulation): the weighting factors that balance the value, gradient, and Hessian equations in the least-squares objective are introduced without an explicit rule or sensitivity study; because the central claim rests on the benefit of the additional linear equations, the lack of justification for these weights is load-bearing for reproducibility of the reported gains.

    Authors: The referee correctly notes that an explicit, reproducible rule for the weighting coefficients is needed. In the current manuscript the weights are chosen to normalize each residual type (value, gradient, Hessian) by its typical magnitude across the training set so that the three blocks contribute on comparable scales; however, this procedure is only described informally. We will revise §4.2 to state the normalization rule in full detail and will add a short sensitivity study (varying the relative weights by factors of 0.1–10) demonstrating that the reported accuracy gains remain stable within a reasonable range around the chosen values. revision: yes

Circularity Check

0 steps flagged · score 0.0 of 10

No significant circularity in derivation chain

full rationale

The paper generates training data for value, gradient, and Hessian of the value function by solving the PMP optimality system and the associated matrix Riccati equation along trajectories; these are standard results from optimal control theory and are used only to supply linear equations for a subsequent weighted least-squares regression over polynomial bases. The regression itself is a conventional fit with no parameters whose values are defined by the target result. The reported improvement in sample complexity is an empirical observation on benchmark problems rather than a mathematical identity obtained by construction. No self-definitional steps, fitted-input predictions, load-bearing self-citations, or ansatz smuggling appear in the derivation.

Assumptions & free parameters 0 free parameters · 1 assumptions · 0 invented entities

The central claim depends on the standard existence and smoothness assumptions that guarantee the PMP optimality system and the associated Riccati equation correctly recover the value function and its derivatives; no new entities are postulated and no free parameters are introduced in the abstract description.

assumptions (1)
  • domain assumption Solutions to the PMP optimality system exist and are sufficiently smooth for the matrix Riccati equation to yield the exact Hessian of the value function along optimal trajectories.
    Data generation step relies on this property to produce accurate second-order information.

how reviews work

0 comments
Cite this review

Pith. "Pith review of Hessian-augmented Supervised Learning for Hamilton-Jacobi-Bellman PDEs." pith.science (2026). https://pith.science/paper/FQ5N3FGZ

@misc{pith2026260623827,
  author       = {Pith},
  title        = {Pith review of: Hessian-augmented Supervised Learning for Hamilton-Jacobi-Bellman PDEs},
  year         = {2026},
  howpublished = {\url{https://pith.science/paper/FQ5N3FGZ}},
  note         = {Machine review of arXiv:2606.23827}
}
read the original abstract

A data-driven method is developed for approximating value functions in deterministic optimal control problems with nonlinear control-affine dynamics. The Pontryagin Maximum Principle optimality system is solved from multiple initial conditions to generate training data consisting of values, gradients, and Hessians of the value function, where Hessian information is obtained from a matrix Riccati equation along optimal trajectories. These quantities augment a weighted least-squares regression over sparse polynomial bases on hyperbolic cross index sets, with gradients and Hessians contributing additional linear equations per sample and substantially reducing sample complexity compared to value-only regression. Feedback laws are recovered analytically from the learned value function. In high dimensions, a partial Hessian strategy controls the cost of data generation. The approach is validated on problems of increasing state dimension, where second-order data augmentation is shown to improve approximation accuracy and closed-loop performance, with up to an order-of-magnitude reduction in the number of training samples required relative to lower-order methods.

Figures

Figures reproduced from arXiv: 2606.23827 by the authors.

Figure 1
Figure 1. Approximation of the Styblinski–Tang function [PITH_FULL_IMAGE:figures/full_fig_p003_1.png] view at source ↗
Figure 2
Figure 2. Relative L 2 , H1 , and H2 errors for the approximation of the 12-dimensional test (4.1) as a function of N, using IHC(4) with q = 3,482. Zeroth-, first-, and second￾order regressions are compared, with partial Hessian fractions ρ ∈ {0.3, 0.5, 0.8, 1.0}. Curves are averaged over 10 independent training sets. Note that each panel reports a different error metric, identified by the subscript of Err, and that the ordin… view at source ↗
Figure 3
Figure 3. Polynomial approximation of the value function for the 2D Van der Pol prob [PITH_FULL_IMAGE:figures/full_fig_p012_3.png] view at source ↗
Figures from the paper (6 more)
Figure 4
Figure 4. Figure 4: Relative validation errors as functions of [PITH_FULL_IMAGE:figures/full_fig_p012_4.png]
Figure 5
Figure 5. Figure 5: Closed-loop trajectories for N = 15 training samples. Each row corresponds to an initial condition x0 ∈ {(2, −1),(2, 1)}. Left: phase portrait. Centre: state norm ∥y(t)∥2. Right: control signal u(t), clipped to [−25, 25]. We next study how approximation quality depends…
Figure 6
Figure 6. Figure 6: Relative approximation errors as a function of basis cardinality [PITH_FULL_IMAGE:figures/full_fig_p014_6.png]
Figure 7
Figure 7. Figure 7: Relative L 2 , H1 , H2 errors and HJB residual for the 6D satellite problem as a function of N, for bases IHC(4) ∪ ITD(5) (q = 663, top row) and IHC(4) ∪ ITD(7) (q = 1,800, bottom row). Hessian fractions ρ ∈ {0.3, 0.5, 0.8, 1.0} included for second￾order regression. Cu…
Figure 8
Figure 8. Figure 8: Closed-loop trajectories for the 6D satellite using partial Hessian regression [PITH_FULL_IMAGE:figures/full_fig_p018_8.png]
Figure 9
Figure 9. Figure 9: Relative validation errors for the 19-dimensional Allen–Cahn system on Ω = [PITH_FULL_IMAGE:figures/full_fig_p019_9.png]

Discussion (0). Continue with ORCID to comment.

Reference graph

Works this paper leans on

34 extracted references · 3 canonical work pages

  1. [1]

    Adcock, S

    B. Adcock, S. Brugiapaglia, and C. G. Webster,Sparse Polynomial Approximation of High- Dimensional Functions, Comput. Sci. Eng., SIAM, Philadelphia, PA, 2022

    2022

  2. [2]

    Adcock and Y

    B. Adcock and Y. Sui,Compressive Hermite interpolation: Sparse, high-dimensional approx- imation from gradient-augmented measurements, Constr. Approx., 50(1):167–207, 2019

    2019

  3. [3]

    G. Albi, S. Bicego, and D. Kalise,Gradient-augmented supervised learning of optimal feedback laws using state-dependent Riccati equations, IEEE Control Syst. Lett., 6:836–841, 2022

    2022

  4. [4]

    B. Azmi, D. Kalise, and K. Kunisch,Optimal feedback law recovery by gradient-augmented sparse polynomial regression, J. Mach. Learn. Res., 22(48):1–32, 2021

    2021

  5. [5]

    Bardi and I

    M. Bardi and I. Capuzzo-Dolcetta,Optimal Control and Viscosity Solutions of Hamilton– Jacobi–Bellman Equations, Systems Control Found. Appl., Birkh¨ auser, Boston, MA, 1997

    1997

  6. [6]

    Breiten, K

    T. Breiten, K. Kunisch, and J. Schr¨ oder,Numerical realization of the Mortensen observer via a Hessian-augmented polynomial approximation of the value function, SIAM J. Sci. Comput., 47(1), 2025

    2025

  7. [7]

    Cannarsa and H

    P. Cannarsa and H. Frankowska,From pointwise to local regularity for solutions of Hamilton– Jacobi equations, Calc. Var. Partial Differ. Equ., 49(3):1061–1074, 2014

    2014

  8. [8]

    Cannarsa, H

    P. Cannarsa, H. Frankowska, and T. Scarinci,Second-order sensitivity relations and regularity of the value function for Mayer’s problem in optimal control, SIAM J. Control Optim., 53(6):3642–3672, 2015

    2015

Show all 34 references
  1. [9]

    Caroff and H

    N. Caroff and H. Frankowska,Conjugate points and shocks in nonlinear optimal control, Trans. Amer. Math. Soc., 348:3133–3153, 1996

    1996

  2. [10]

    Caroff and H

    N. Caroff and H. Frankowska,Optimality and characteristics of Hamilton–Jacobi–Bellman equations, in Optim., Optim. Control Partial Differ. Equ., Internat. Ser. Numer. Math. 107, pp. 169–180, Birkh¨ auser, Basel, 1992

    1992

  3. [11]

    F. H. Clarke and R. B. Vinter,Local optimality conditions and Lipschitzian solutions to the Hamilton–Jacobi equation, SIAM J. Control Optim., 21:856–870, 1983

    1983

  4. [12]

    F. H. Clarke and R. B. Vinter,The relationship between the maximum principle and dynamic programming, SIAM J. Control Optim., 25(5):1291–1311, 1987. HESSIAN-AUGMENTED SUPERVISED LEARNING FOR HJB PDES21

    1987

  5. [13]

    N. A. Corbin and B. Kramer,Scalable computation ofH ∞ energy functions for polynomial control-affine systems, IEEE Trans. Autom. Control, 70(5):3088–3100, 2025

    2025

  6. [14]

    N. A. Corbin and B. Kramer,Computing solutions to the polynomial-polynomial regulator problem, inProc. 63rd IEEE Conf. Decis. Control, pp. 3073–3078, 2024

    2024

  7. [15]

    W. M. Czarnecki, S. Osindero, M. Jaderberg, G. Swirszcz, and R. Pascanu,Sobolev training for neural networks, Adv. Neural Inf. Process. Syst., vol. 30, 2017

    2017

  8. [16]

    Darbon and S

    J. Darbon and S. Osher,Algorithms for overcoming the curse of dimensionality for certain Hamilton–Jacobi equations arising in control theory and elsewhere, Res. Math. Sci., 3:19, 2016

    2016

  9. [17]

    Dolgov, D

    S. Dolgov, D. Kalise, and K. Kunisch,Tensor decomposition methods for high-dimensional Hamilton–Jacobi–Bellman equations, SIAM J. Sci. Comput., 43(3):A1625–A1650, 2021

    2021

  10. [18]

    Dolgov, D

    S. Dolgov, D. Kalise, and L. Saluzzi,Data-driven tensor train gradient cross approximation for Hamilton–Jacobi–Bellman equations, SIAM J. Sci. Comput., 45(5):A2153–A2184, 2023

    2023

  11. [19]

    Ehring, B

    T. Ehring, B. Azmi, and B. Haasdonk,Recovery of the optimal control value function in reproducing kernel Hilbert spaces from verification conditions, arXiv:2512.07477, 2025

    2025

  12. [20]

    Ehring and B

    T. Ehring and B. Haasdonk,Hermite kernel surrogates for the value function of high- dimensional nonlinear optimal control problems, Adv. Comput. Math., 50(3):36, 2024

    2024

  13. [21]

    Esteve-Yag¨ ue, R

    C. Esteve-Yag¨ ue, R. Tsai, and A. Massucco,Finite-difference least square methods for solving Hamilton–Jacobi equations using neural networks, arXiv:2406.10758, 2024

    2024

  14. [22]

    Gr¨ une and A

    L. Gr¨ une and A. Rantzer,On the infinite horizon performance of receding horizon controllers, IEEE Trans. Autom. Control, 53(9):2100–2111, 2008

    2008

  15. [23]

    J. Han, A. Jentzen, and W. E,Solving high-dimensional partial differential equations using deep learning, Proc. Natl. Acad. Sci. USA, 115(34):8505–8510, 2018

    2018

  16. [24]

    Kalise and K

    D. Kalise and K. Kunisch,Polynomial approximation of high-dimensional Hamilton–Jacobi– Bellman equations and applications to feedback control of semilinear parabolic PDEs, SIAM J. Sci. Comput., 40(2):A629–A652, 2018

    2018

  17. [25]

    Kang and L

    W. Kang and L. C. Wilcox,A causality free computational method for HJB equations with application to rigid body satellites, inAIAA Guid. Navig. Control Conf., pp. 1–10, 2015

    2015

  18. [26]

    Kang and L

    W. Kang and L. C. Wilcox,Mitigating the curse of dimensionality: Sparse grid characteristics method for optimal feedback control and HJB equations, Comput. Optim. Appl., 68:289– 315, 2017

    2017

  19. [27]

    Kunisch and D

    K. Kunisch and D. Walter,Semiglobal optimal feedback stabilization of autonomous systems via deep neural network approximation, ESAIM Control Optim. Calc. Var., 27:16, 2021

    2021

  20. [28]

    T. Meng, Z. Zhang, J. Darbon, and G. Karniadakis,SympOCnet: Solving optimal control problems with applications to high-dimensional multiagent path planning problems, SIAM J. Sci. Comput., 44(6):B1341–B1368, 2022

    2022

  21. [29]

    Nakamura-Zimmerer, Q

    T. Nakamura-Zimmerer, Q. Gong, and W. Kang,Adaptive deep learning for high-dimensional Hamilton–Jacobi–Bellman equations, SIAM J. Sci. Comput., 43(2):A1221–A1247, 2021

    2021

  22. [30]

    L. S. Pontryagin, V. G. Boltyanskii, R. V. Gamkrelidze, and E. F. Mishchenko,The Mathe- matical Theory of Optimal Processes, Interscience, Wiley, New York, 1962

    1962

  23. [31]

    Sirignano and K

    J. Sirignano and K. Spiliopoulos,DGM: A deep learning algorithm for solving partial differ- ential equations, J. Comput. Phys., 375:1339–1364, 2018

    2018

  24. [32]

    Sperl, L

    M. Sperl, L. Saluzzi, D. Kalise, and L. Gr¨ une,Separable approximations of optimal value functions and their representation by neural networks, arXiv:2502.08559, 2025

    2025

  25. [33]

    N. N. Subbotina,The method of characteristics for Hamilton–Jacobi equations and applica- tions to dynamical optimization, J. Math. Sci., 135:2955–3091, 2006

    2006

  26. [34]

    Yegorov and P

    I. Yegorov and P. M. Dower,Perspectives on characteristics-based curse-of-dimensionality- free numerical approaches for solving Hamilton–Jacobi equations, Appl. Math. Optim., 81:1167–1218, 2021

    2021

Pith tools

Reviewed June 26, 2026 · model on record in the stance chip above.