arxiv: 2604.26933 · v1 · submitted 2026-04-29 · 🧮 math.DS

Recognition: unknown

Data-driven discovery of polynomial ODEs with provably bounded solutions

Albert Alcalde, Giovanni Fantuzzi

Pith reviewed 2026-05-07 09:21 UTC · model grok-4.3

classification 🧮 math.DS

keywords data-driven discoverypolynomial ODEsLyapunov functionsbounded trajectoriessystem identificationnonlinear dynamical systemsoptimization methods

0 comments

The pith

A data-driven method recovers polynomial ODEs with provably bounded trajectories by jointly finding the vector field and a certifying Lyapunov function.

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

The paper aims to show that it is possible to learn polynomial ordinary differential equations from data in such a way that the learned models are guaranteed to have all trajectories bounded. This is achieved by simultaneously determining the polynomial right-hand side and a polynomial function that serves as a Lyapunov function, thereby defining a compact absorbing set. A sympathetic reader would care because data-driven models without such guarantees can lead to simulations that diverge or behave unrealistically, limiting their reliability in scientific and engineering contexts. The approach solves the resulting nonconvex optimization by alternating between convex subproblems for the vector field coefficients and the Lyapunov coefficients, starting from a data-based initialization. This extends earlier techniques limited to quadratic forms and ellipsoidal bounds to general polynomials.

Core claim

The paper establishes that polynomial ODEs with provably bounded trajectories can be discovered from data through a well-posed nonconvex optimization problem that is solved via an alternating block-coordinate scheme with convex subproblems, enabled by a novel initialization that identifies a candidate Lyapunov function directly from the data.

What carries the argument

The alternating block-coordinate optimization scheme that alternately solves for the coefficients of the polynomial vector field and the polynomial Lyapunov function, with feasibility guaranteed by the initialization procedure.

Load-bearing premise

The load-bearing premise is that the data come from a system whose behavior is well approximated by a low-degree polynomial vector field for which a polynomial Lyapunov function of modest degree can be found through the alternating optimization.

What would settle it

Finding a dynamical system where the true trajectories are bounded but the method either fails to recover a model or returns one without a certifying Lyapunov function, or conversely recovers a model that the method certifies as bounded but which actually has unbounded trajectories.

Figures

Figures reproduced from arXiv: 2604.26933 by Albert Alcalde, Giovanni Fantuzzi.

**Figure 2.** Figure 2: Learned vector fields and corresponding absorbing sets (shaded in green) for Lyapunov function degree dv = 2 and model degrees df = 2, 3, 4, 5. Panels in the bottom row show a zoomed-in view of the box [−3, 3]2 . limit cycle and a certified absorbing ellipsoid. Thus, while the true system does not admit a quadratic Lyapunov function, it can be approximated by one that does, which SILAS recovers. Increasing… view at source ↗

**Figure 3.** Figure 3: Learned vector fields and corresponding absorbing sets (shaded in green) for Lyapunov function degree dv = 4 and model degrees df = 2, . . . , 9. The bottom row shows zoomed-out views of the panels in the middle row. dimension up to n = 6 and no reported unbounded variables. Of these 109 ODEs, 82 have a polynomial vector field. Some of them, such as the Lorenz system, are known to admit absorbing sets. Oth… view at source ↗

**Figure 4.** Figure 4: Pointwise 2-norm errors between the exact vector field of the cubic system in Section 3.1 and the vector fields learned with Lyapunov function degree dv = 4 and model degree df = 2, . . . , 9. fixing Λ = 0.1Λ0 and µ = 0.1µ0, where Λ0 and µ0 are selected such that the function z 7→ Λ −1 0 (z − µ0) maps the hypercube [−1, 1]n to a box bounding the training data with a 25% buffer. Precisely, if x max and x mi… view at source ↗

**Figure 5.** Figure 5: Error distributions (curves) and mean errors (vertical bars) for 82 polynomial ODEs in the dysts library. Colors identify models of degree 2 (■), 3 (■), 4 (■), 5 (■), 6 (■), and 7 (■). Green shading indicates errors smaller than 1%. error is smaller than 1% for 71 of the 82 polynomial systems and for 10 of the 27 nonpolynomial ones. Thus, SILAS works robustly for systems governed by polynomial equations, … view at source ↗

**Figure 6.** Figure 6: Error distributions (curves) and mean errors (vertical bars) for 27 non-polynomial ODEs in the dysts library. Colors identify models of degree 2 (■), 3 (■), 4 (■), 5 (■), 6 (■), and 7 (■). Green shading indicates errors smaller than 1% view at source ↗

**Figure 7.** Figure 7: Attractors for selected systems from the dysts library. Top row: The true attractor. Middle row: Discovered attractors. Bottom row: Discovered attractors and absorbing sets certifying the boundedness of the discovered ODEs. the discovered models are also shown. In all cases, the discovered attractors are topologically similar to the exact ones, showing that SILAS can work well even when the true system d… view at source ↗

**Figure 8.** Figure 8: Attractors for selected systems from the dysts library, discovered with a refined dataset. Top row: True attractors. Middle row: Discovered attractors. Bottom row: Discovered attractors and absorbing sets certifying the boundedness of the discovered ODEs. there are also cases in which our bounded ODE models are unable to accurately reproduce the true system dynamics, such as the non-polynomial systems Sp… view at source ↗

**Figure 9.** Figure 9: Modeled vs. true POD-coefficient dynamics for the PDE example in Section 3.3, projected onto the (x1, x2) plane. The reduced-order models are for n = 3 (left), 4 (middle), and 5 (right) POD modes. Panels in each column show the PDE dynamics (top), the reducedorder ODE dynamics (middle), and the reduced-order ODE dynamics with their absorbing ellipsoid (bottom). with df = 3 for n = 5 view at source ↗

**Figure 10.** Figure 10: Reconstructed vs. true r-component of the PDE solution in Section 3.3, obtained from reduced-order models with n = 3 (top), 4 (middle), and 5 (bottom) POD modes. In each row, panels show the reference PDE simulation (left), the reconstruction from the learned ODE model (middle), and the absolute error between the two (right). Lyapunov functions characterizing compact absorbing sets. Concretely, SILAS iden… view at source ↗

read the original abstract

We introduce SILAS, a data-driven framework for discovering polynomial ordinary differential equations (ODEs) with provably bounded trajectories. Boundedness is certified by compact absorbing sets defined via polynomial Lyapunov functions. We jointly identify the ODE vector field and the Lyapunov function using a well-posed nonconvex optimization problem built using polynomial optimization tools. We solve this problem using an alternating block-coordinate optimization scheme with convex subproblems, whose feasibility is ensured by a novel model-agnostic initialization that identifies a candidate Lyapunov function from data. Our methods extend prior approaches for quadratic ODEs with absorbing ellipsoids to a significantly broader class of ODEs and absorbing sets. A suite of over 100 examples demonstrates that SILAS can recover accurate and provably bounded ODE models for a broad range of nonlinear dynamical systems.

Editorial analysis

A structured set of objections, weighed in public.

Desk editor's note, referee report, simulated authors' rebuttal, and a circularity audit. Tearing a paper down is the easy half of reading it; the pith above is the substance, this is the friction.

Desk Editor's Note private letter to a colleague

SILAS gives a practical alternating-optimization route to polynomial ODEs plus Lyapunov certificates, with decent empirical coverage but no convergence theory.

read the letter

The core contribution is a named procedure, SILAS, that simultaneously fits a polynomial vector field and a polynomial Lyapunov function to time-series data. It uses block-coordinate descent on a nonconvex polynomial-optimization problem, with a model-agnostic initialization that starts the Lyapunov block from a feasible point derived directly from the data. This extends earlier quadratic-ODE work to higher-degree polynomials and more general compact absorbing sets. The paper reports success on more than 100 synthetic and real examples, recovering models whose trajectories stay inside the certified sets. That empirical breadth is the strongest part of the manuscript; the examples span several standard nonlinear benchmarks and the recovered degrees are modest enough to be useful in control applications. The initialization trick is also a clean idea that removes one common source of infeasibility in these joint problems. The soft spot is the lack of any convergence or basin analysis for the alternation. Each subproblem is convex, but the joint problem is not, so stationary points can exist where the data fit is poor or the Lyapunov certificate is only marginally valid before rounding. The manuscript does not quantify how often the procedure lands in such points, nor does it supply error bars on the recovered coefficients or a systematic study of degree-selection sensitivity. The 100-example suite therefore shows feasibility rather than reliability across the broader class claimed in the abstract. Readers working on data-driven system identification or certified nonlinear control will find the method worth trying; the code and examples appear reproducible enough to test quickly. The work is coherent on its own terms and engages the relevant polynomial-optimization literature, so it deserves a serious referee. I would send it out, with the expectation that reviewers will press on the optimization guarantees and on how the method behaves when the underlying system is only approximately polynomial.

Referee Report

2 major / 2 minor

Summary. The manuscript introduces SILAS, a data-driven framework that jointly identifies a polynomial vector field and a polynomial Lyapunov function from trajectory data via a nonconvex optimization problem solved by alternating block-coordinate descent on convex subproblems. Boundedness is certified by the existence of a compact absorbing set defined by the Lyapunov function. The method is initialized with a model-agnostic procedure that finds an initial feasible Lyapunov function, and the approach is validated empirically on a suite of over 100 examples spanning nonlinear dynamical systems.

Significance. If the alternating procedure reliably recovers accurate polynomial models together with valid Lyapunov certificates, the work would meaningfully extend data-driven discovery beyond quadratic ODEs with ellipsoidal bounds, providing explicit safety guarantees useful for control and verification. The scale of the numerical demonstration and the convex-subproblem structure are positive features; however, the absence of convergence analysis or failure-mode quantification limits the strength of the central claim that the method works for the broader class invoked in the abstract.

major comments (2)

[optimization procedure (alternating scheme)] The alternating block-coordinate scheme for the joint nonconvex problem (described in the optimization section following the problem formulation) lacks any convergence analysis, basin characterization, or quantification of how often it reaches points where the data-fit residual remains large or the Lyapunov certificate becomes invalid after rounding. This is load-bearing for the claim that SILAS recovers accurate and provably bounded models on the reported suite, because each block subproblem being convex does not guarantee that alternation avoids stationary points that are feasible yet inaccurate.
[numerical results / examples] No error bars, sensitivity analysis with respect to polynomial degree selection, or statistics on how often the post-hoc degree choice affects certificate validity are provided in the numerical results section. This weakens the assertion of success across more than 100 examples, as the central claim relies on the procedure producing both accurate dynamics and valid certificates without quantifying variability or failure rates.

minor comments (2)

[problem formulation] The notation for the polynomial bases and the precise definition of the absorbing-set radius should be stated more explicitly when first introduced, to avoid ambiguity when comparing to prior quadratic-ODE work.
[figures] A few figure captions in the examples section could more clearly indicate which trajectories are training data versus validation trajectories.

Simulated Author's Rebuttal

2 responses · 2 unresolved

We thank the referee for the constructive review and positive remarks on the framework's scope and numerical demonstration. We respond to each major comment below, indicating where revisions are planned.

read point-by-point responses

Referee: The alternating block-coordinate scheme for the joint nonconvex problem (described in the optimization section following the problem formulation) lacks any convergence analysis, basin characterization, or quantification of how often it reaches points where the data-fit residual remains large or the Lyapunov certificate becomes invalid after rounding. This is load-bearing for the claim that SILAS recovers accurate and provably bounded models on the reported suite, because each block subproblem being convex does not guarantee that alternation avoids stationary points that are feasible yet inaccurate.

Authors: We acknowledge the absence of a theoretical convergence analysis for the alternating scheme. Each block subproblem is convex and solved to global optimality, while the model-agnostic initialization explicitly constructs a feasible starting Lyapunov function from data. The joint problem remains nonconvex, so global convergence or basin characterization is not available. In revision we will expand the optimization section with a dedicated paragraph discussing these limitations, reporting observed iteration counts and residual behavior across the example suite, and clarifying that the method relies on empirical reliability rather than provable convergence to accurate models. revision: partial
Referee: No error bars, sensitivity analysis with respect to polynomial degree selection, or statistics on how often the post-hoc degree choice affects certificate validity are provided in the numerical results section. This weakens the assertion of success across more than 100 examples, as the central claim relies on the procedure producing both accurate dynamics and valid certificates without quantifying variability or failure rates.

Authors: The numerical results section reports success on more than 100 examples without error bars, degree-sensitivity studies, or aggregate statistics on post-rounding certificate validity. We agree that such quantification would strengthen the presentation. In the revision we will add a short subsection (or appendix) containing sensitivity results for a representative subset of systems, including the effect of degree choice on certificate validity, and report the fraction of examples in which the rounded Lyapunov function remains valid. Full error bars and exhaustive failure-mode statistics for every example would require substantial additional computation beyond the current scope. revision: partial

standing simulated objections not resolved

Theoretical convergence analysis or basin characterization for the nonconvex alternating block-coordinate scheme
Comprehensive error bars, sensitivity analysis, and failure-rate statistics across the entire suite of more than 100 examples

Circularity Check

0 steps flagged

No significant circularity in the derivation chain

full rationale

The paper presents SILAS as a new algorithmic framework that formulates a joint nonconvex optimization problem over polynomial vector-field coefficients and a polynomial Lyapunov function, then solves it via alternating convex block-coordinate steps initialized by a model-agnostic procedure. No step in the abstract or described method reduces a claimed prediction, uniqueness result, or boundedness certificate to a fitted parameter or self-citation by construction; the extension of prior quadratic-ODE work is stated as a generalization rather than a load-bearing premise that defines the new result. The central contribution is therefore an empirically validated procedure whose validity rests on the optimization formulation and numerical examples, not on any tautological re-expression of its inputs.

Axiom & Free-Parameter Ledger

1 free parameters · 1 axioms · 0 invented entities

The central claim rests on the existence of a polynomial vector field and a polynomial Lyapunov function of modest degree that can be recovered by the proposed nonconvex program; the abstract does not enumerate explicit free parameters or invented entities.

free parameters (1)

polynomial degree for vector field and Lyapunov function
Chosen to balance expressivity and tractability; not specified numerically in abstract.

axioms (1)

domain assumption The underlying system admits a polynomial approximation whose trajectories remain inside a compact absorbing set certified by a polynomial Lyapunov function.
Invoked implicitly when claiming provable boundedness for the discovered models.

pith-pipeline@v0.9.0 · 5425 in / 1254 out tokens · 39489 ms · 2026-05-07T09:21:18.118074+00:00 · methodology

discussion (0)

Reference graph

Works this paper leans on

64 extracted references · 46 canonical work pages · 2 internal anchors

[1]

Abyaneh and H.-C

A. Abyaneh and H.-C. Lin. Learning Lyapunov-stable polynomial dynamical systems through imitation. In J. Tan, M. Toussaint, and K. Darvish, editors,Proc. 7th Conf. Robot Learn., volume 229 ofProc. Mach. Learn. Res., pages 2642–2662. PMLR, 2023. DATA-DRIVEN DISCOVERY OF BOUNDED POLYNOMIAL ODES 21

2023
[2]

Jay Shah, Ganesh Bikshandi, Ying Zhang, Vijay Thakkar, Pradeep Ramani, and Tri Dao

A. Abyaneh, M. S. Guzmán, and H.-C. Lin. Globally stable neural imitation policies. In2024 IEEE Int. Conf. Robot. Autom. (ICRA), pages 15061–15067. IEEE, 2024. doi: 10.1109/ICRA57147.2024.10610791

work page doi:10.1109/icra57147.2024.10610791 2024
[3]

A. A. Ahmadi and B. E. Khadir. Learning dynamical systems with side information. SIAM Review, 65(1):183–223, 2023. doi: 10.1137/20M1388644

work page doi:10.1137/20m1388644 2023
[4]

ApS.MOSEK Optimization Toolbox for MATLAB 11.1.11, 2026

M. ApS.MOSEK Optimization Toolbox for MATLAB 11.1.11, 2026

2026
[5]

Brenner , doi =

Y.Bar-Sinai, S.Hoyer, J.Hickey, andM.P.Brenner. Learningdata-drivendiscretizations for partial differential equations.Proc. Natl. Acad. Sci. USA, 116(31):15344–15349, 2019. doi: 10.1073/pnas.1814058116

work page doi:10.1073/pnas.1814058116 2019
[6]

Beck and L

A. Beck and L. Tetruashvili. On the convergence of block coordinate descent type meth- ods.SIAM J. Optim., 23(4):2037–2060, 2013. doi: 10.1137/120887679

work page doi:10.1137/120887679 2037
[7]

D. P. Bertsekas.Nonlinear programming. Athena Scientific Optimization and Computa- tion Series. Athena Scientific, Belmont, MA, second edition, 1999. ISBN 1-886529-00-0

1999
[8]

Bou-Sakr-El-Tayar, J

M. Bou-Sakr-El-Tayar, J. J. Bramburger, and M. J. Colbrook. Weighted Birkhoff aver- ages accelerate data-driven methods. arXiv:2511.17772 [math.DS], 2025

work page arXiv 2025
[9]

J. J. Bramburger and G. Fantuzzi. Auxiliary functions as Koopman observables: data- driven analysis of dynamical systems via polynomial optimization.J. Nonlinear Sci., 34 (1):Paper No. 8, 38, 2024. doi: 10.1007/s00332-023-09990-2

work page doi:10.1007/s00332-023-09990-2 2024
[10]

J. J. Bramburger, S. L. Brunton, and J. Nathan Kutz. Deep learning of conjugate mappings.Phys. D, 427:Paper No. 133008, 16, 2021. doi: 10.1016/j.physd.2021.133008

work page doi:10.1016/j.physd.2021.133008 2021
[11]

M. M. Bronstein, J. Bruna, T. Cohen, and P. Veličković. Geometric deep learning: Grids, groups, graphs, geodesics, and gauges. arXiv:2104.13478 [cs.LG], 2021

work page internal anchor Pith review arXiv 2021
[12]

K. J. Burns, G. M. Vasil, J. S. Oishi, D. Lecoanet, and B. P. Brown. Dedalus: A flexible framework for numerical simulations with spectral methods.Phys. Rev. Res., 2:023068, Apr 2020. doi: 10.1103/PhysRevResearch.2.023068

work page doi:10.1103/physrevresearch.2.023068 2020
[13]

Cunis and J

T. Cunis and J. Olucak. CaΣos: A nonlinear sum-of-squares optimization suite. In2025 American Control Conf. (ACC), pages 1659–1666, 2025. doi: 10.23919/ACC63710.2025. 11107794

work page doi:10.23919/acc63710.2025 2025
[14]

C. R. Doering and J. D. Gibbon.Applied analysis of the Navier-Stokes equations. Cam- bridge Texts in Applied Mathematics. Cambridge University Press, Cambridge, 1995. ISBN 0-521-44557-4; 0-521-44568-X. doi: 10.1017/CBO9780511608803

work page doi:10.1017/cbo9780511608803 1995
[15]

T. A. Driscoll, N. Hale, and L. N. Trefethen.Chebfun Guide. Pafnuty Publications, 2014

2014
[16]

Elezgaray and A

J. Elezgaray and A. Arneodo. Crisis-induced intermittent bursting in reaction-diffusion chemical systems.Phys. Rev. Lett., 68:714–717, 1992. doi: 10.1103/PhysRevLett.68.714

work page doi:10.1103/physrevlett.68.714 1992
[17]

Fantuzzi, D

G. Fantuzzi, D. Goluskin, D. Huang, and S. I. Chernyshenko. Bounds for deterministic and stochastic dynamical systems using sum-of-squares optimization.SIAM J. Appl. Dyn. Syst., 15(4):1962–1988, 2016. doi: 10.1137/15M1053347

work page doi:10.1137/15m1053347 1962
[18]

W. Gilpin. Chaos as an interpretable benchmark for forecasting and data-driven mod- elling. In35th Conf. Neural Inform. Process. Syst. Datasets and Benchmarks Track (Round 2), 2021

2021
[19]

Goertzen, S

A. Goertzen, S. Tang, and N. Azizan. Learning chaotic PDEs with boundedness guar- antees. InNeurIPS 2025 AI for Science Workshop, 2025

2025
[20]

Goluskin

D. Goluskin. Bounding extrema over global attractors using polynomial optimisation. Nonlinearity, 33(9):4878–4899, 2020. doi: 10.1088/1361-6544/ab8f7b

work page doi:10.1088/1361-6544/ab8f7b 2020
[21]

Goyal and P

P. Goyal and P. Benner. Discovery of nonlinear dynamical systems using a Runge- Kutta inspired dictionary-based sparse regression approach.Proc. A, 478(2262):Paper No. 20210883, 24, 2022. doi: 10.1098/rspa.2021.0883. 22 DATA-DRIVEN DISCOVERY OF BOUNDED POLYNOMIAL ODES

work page doi:10.1098/rspa.2021.0883 2022
[22]

Goyal, I

P. Goyal, I. Pontes Duff, and P. Benner. Guaranteed stable quadratic models and their applications in SINDy and operator inference.Phys. D, 483:Paper No. 134893, 16, 2025. doi: 10.1016/j.physd.2025.134893

work page doi:10.1016/j.physd.2025.134893 2025
[23]

Greydanus, M

S. Greydanus, M. Dzamba, and J. Yosinski. Hamiltonian neural networks. InProc. 33rd Int. Conf. Neural Inform. Proc. Syst., Red Hook, NY, USA, 2019. Curran Associates Inc

2019
[24]

Grippo and M

L. Grippo and M. Sciandrone. Globally convergent block-coordinate techniques for un- constrained optimization.Optim. Methods Softw., 10(4):587–637, 1999. doi: 10.1080/ 10556789908805730

1999
[25]

A. L. Heide, S.-C. Liao, S. Castiblanco-Ballesteros, G. B. Jacobs, P. Seiler, and M. S. Hemati. Data-driven nonlinear aerodynamics models with certifiably optimal bounded- ness properties. arXiv:2508.16800 [physics.flu-dyn], 2025

work page arXiv 2025
[26]

Henrion and M

D. Henrion and M. Korda. Convex computation of the region of attraction of polynomial controlsystems.IEEE Trans. Automat. Control, 59(2):297–312, 2014. doi: 10.1109/TAC. 2013.2283095

work page doi:10.1109/tac 2014
[27]

D. Hilbert. Ueber die Darstellung definiter Formen als Summe von Formenquadraten. Math. Ann., 32(3):342–350, 1888. doi: 10.1007/BF01443605

work page doi:10.1007/bf01443605
[28]

Holmes, J

P. Holmes, J. L. Lumley, G. Berkooz, and C. W. Rowley.Turbulence, coherent structures, dynamical systems and symmetry. Cambridge Monographs on Mechanics. Cambridge University Press, Cambridge, second edition, 2012. ISBN 978-1-107-00825-0. doi: 10. 1017/CBO9780511919701

2012
[29]

A. A. Kaptanoglu, J. L. Callaham, A. Aravkin, C. J. Hansen, and S. L. Brunton. Pro- moting global stability in data-driven models of quadratic nonlinear dynamics.Phys. Rev. Fluids, 6(9):094401, 2021. doi: 10.1103/PhysRevFluids.6.094401

work page doi:10.1103/physrevfluids.6.094401 2021
[30]

S. M. Khansari-Zadeh and A. Billard. Learning stable nonlinear dynamical systems with gaussian mixture models.IEEE Trans. Robotics, 27(5):943–957, 2011. doi: 10.1109/ TRO.2011.2159412

work page arXiv 2011
[31]

S. Klus, P. Koltai, and C. Schütte. On the numerical approximation of the Perron- Frobenius and Koopman operator.J. Comput. Dyn., 3(1):51–79, 2016. doi: 10.3934/ jcd.2016003

2016
[32]

S. Klus, F. Nüske, S. Peitz, J.-H. Niemann, C. Clementi, and C. Schütte. Data-driven approximation of the Koopman generator: Model reduction, system identification, and control.Phys. D, 406:132416, 15, 2020. doi: 10.1016/j.physd.2020.132416

work page doi:10.1016/j.physd.2020.132416 2020
[33]

M. Korda. Computing controlled invariant sets from data using convex optimization. SIAM J. Control Optim., 58(5):2871–2899, 2020. doi: 10.1137/19M1305835

work page doi:10.1137/19m1305835 2020
[34]

Korda and I

M. Korda and I. Mezić. On convergence of extended dynamic mode decomposition to the Koopman operator.J. Nonlinear Sci., 28(2):687–710, 2018. doi: 10.1007/ s00332-017-9423-0

2018
[35]

Korda, D

M. Korda, D. Henrion, and C. N. Jones. Convex computation of the maximum controlled invariant set for polynomial control systems.SIAM J. Control Optim., 52(5):2944–2969,
[36]

doi: 10.1137/130914565

work page doi:10.1137/130914565
[37]

J. B. Lasserre. Global optimization with polynomials and the problem of moments. SIAM J. Optim., 11(3):796–817, 2000/01. doi: 10.1137/S1052623400366802

work page doi:10.1137/s1052623400366802 2000
[38]

Imperial College Press, London, 2010

J.B.Lasserre.Moments, positive polynomials and their applications, volume1ofImperial College Press Optimization Series. Imperial College Press, London, 2010. ISBN 978-1- 84816-445-1; 1-84816-445-9

2010
[39]

J. B. Lasserre, D. Henrion, C. Prieur, and E. Trélat. Nonlinear optimal control via occupation measures and LMI-relaxations.SIAM J. Control Optim., 47(4):1643–1666, DATA-DRIVEN DISCOVERY OF BOUNDED POLYNOMIAL ODES 23
[40]

doi: 10.1137/070685051

work page doi:10.1137/070685051
[41]

M. Laurent. Sums of squares, moment matrices and optimization over polynomi- als. InEmerging applications of algebraic geometry, volume 149 ofIMA Vol. Math. Appl., pages 157–270. Springer, New York, 2009. ISBN 978-0-387-09685-8. doi: 10.1007/978-0-387-09686-5_7

work page doi:10.1007/978-0-387-09686-5_7 2009
[42]

S.-C. Liao, A. Leonid Heide, M. S. Hemati, and P. J. Seiler. A convex optimization approach to compute trapping regions for lossless quadratic systems.Internat. J. Robust Nonlinear Control, 35(6):2425–2436, 2025. doi: 10.1002/rnc.7807

work page doi:10.1002/rnc.7807 2025
[43]

J. Lofberg. YALMIP: A toolbox for modeling and optimization in MATLAB. In2004 IEEE Int. Conf. Robotics & Automation, pages 284–289, 2004. doi: 10.1109/CACSD. 2004.1393890

work page doi:10.1109/cacsd 2004
[44]

Mauroy and A

A. Mauroy and A. Sootla. Estimation of regions of attraction with formal certificates in a purely data-driven setting. In62nd IEEE Conf. Decis. Control (CDC), pages 4682–4687,
[45]

doi: 10.1109/CDC49753.2023.10383454

work page doi:10.1109/cdc49753.2023.10383454 2023
[46]

Mauroy, A

A. Mauroy, A. Sootla, and I. Mezić. Koopman framework for global stability analysis. InThe Koopman operator in systems and control—concepts, methodologies and applica- tions, volume 484 ofLect. Notes Control Inf. Sci., pages 35–58. Springer, Cham, 2020. ISBN 978-3-030-35712-2. doi: 10.1007/978-3-030-35713-9_2

work page doi:10.1007/978-3-030-35713-9_2 2020
[47]

K. G. Murty and S. N. Kabadi. Some NP-complete problems in quadratic and nonlinear programming.Math. Programming, 39(2):117–129, 1987. doi: 10.1007/BF02592948

work page doi:10.1007/bf02592948 1987
[48]

Nemirovski

A. Nemirovski. Advances in convex optimization: conic programming. InInternational Congress of Mathematicians. Vol. I, pages 413–444. Eur. Math. Soc., Zürich, 2007. ISBN 978-3-03719-022-7. doi: 10.4171/022-1/17

work page doi:10.4171/022-1/17 2007
[49]

Nesterov

Y. Nesterov. Squared functional systems and optimization problems. InHigh perfor- mance optimization, volume 33 ofAppl. Optim., pages 405–440. Kluwer Acad. Publ., Dordrecht, 2000. ISBN 0-7923-6013-3. doi: 10.1007/978-1-4757-3216-0_17

work page doi:10.1007/978-1-4757-3216-0_17 2000
[50]

Nesterov and A

Y. Nesterov and A. Nemirovskii.Interior-point polynomial algorithms in convex pro- gramming, volume 13 ofSIAM Studies in Applied Mathematics. Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, 1994. ISBN 0-89871-319-6. doi: 10.1137/1.9781611970791

work page doi:10.1137/1.9781611970791 1994
[51]

On the Practical Implementation of a Sequential Quadratic Programming Algorithm for Nonconvex Sum-of-squares Problems

J. Olucak and T. Cunis. Sequential quadratic sum-of-squares programming for nonlinear control systems. arXiv:2602.02394 [math.OC], 2025

work page internal anchor Pith review Pith/arXiv arXiv 2025
[52]

S. E. Otto, N. Zolman, J. N. Kutz, and S. L. Brunton. A unified framework to enforce, discover, and promote symmetry in machine learning.J. Mach. Learn. Res., 26:Paper No. [248], 83, 2025

2025
[53]

Papachristodoulou and S

A. Papachristodoulou and S. Prajna. A tutorial on sum of squares techniques for systems analysis. InProc. 2005 American Control Conf., pages 2686–2700 vol. 4, 2005. doi: 10.1109/ACC.2005.1470374

work page doi:10.1109/acc.2005.1470374 2005
[54]

P. A. Parrilo. Semidefinite programming relaxations for semialgebraic problems.Math. Program., 96(2, Ser. B):293–320, 2003. doi: 10.1007/s10107-003-0387-5

work page doi:10.1007/s10107-003-0387-5 2003
[55]

Polynomialoptimization, sumsofsquares, andapplications

P.A.Parrilo. Polynomialoptimization, sumsofsquares, andapplications. InSemidefinite optimization and convex algebraic geometry, volume 13 ofMOS-SIAM Ser. Optim., pages 47–157. SIAM, Philadelphia, PA, 2013. ISBN 978-1-611972-28-3

2013
[56]

M. Peng, A. A. Kaptanoglu, C. J. Hansen, J. Stevens-Haas, K. Manohar, and S. L. Brun- ton. Extending the trapping theorem to provide local stability guarantees for quadrati- cally nonlinear models.Phys. Fluids, 37(10):107115, 10 2025. doi: 10.1063/5.0287432

work page doi:10.1063/5.0287432 2025
[57]

Rindler.Calculus of variations

F. Rindler.Calculus of variations. Universitext. Springer, Cham, 2018. ISBN 978-3-319- 77636-1; 978-3-319-77637-8. doi: 10.1007/978-3-319-77637-8. 24 DATA-DRIVEN DISCOVERY OF BOUNDED POLYNOMIAL ODES

work page doi:10.1007/978-3-319-77637-8 2018
[58]

Schlegel and B

M. Schlegel and B. R. Noack. On long-term boundedness of Galerkin models.J. Fluid Mech., 765:325–352, 2015. doi: 10.1017/jfm.2014.736

work page doi:10.1017/jfm.2014.736 2015
[59]

Tacchi, Y

M. Tacchi, Y. Lian, and C. N. Jones. Robustly learning regions of attraction from fixed data.IEEE Trans. Automat. Control, 70(3):1576–1591, 2025. doi: 10.1109/tac.2024. 3462528

work page doi:10.1109/tac.2024 2025
[60]

S. Tang, T. Sapsis, and N. Azizan. Learning dissipative chaotic dynamics with bound- edness guarantees. arXiv:2410.00976 [eess.SY], 2024

work page arXiv 2024
[61]

Herzog, T

U. Topcu, A. Packard, and P. Seiler. Local stability analysis using simulations and sum- of-squares programming.Automatica J. IFAC, 44(10):2669–2675, 2008. doi: 10.1016/j. automatica.2008.03.010

work page doi:10.1016/j 2008
[62]

M. O. Williams, I. G. Kevrekidis, and C. W. Rowley. A data-driven approximation of the Koopman operator: extending dynamic mode decomposition.J. Nonlinear Sci., 25 (6):1307–1346, 2015. doi: 10.1007/s00332-015-9258-5

work page doi:10.1007/s00332-015-9258-5 2015
[63]

Ye.Interior Point Algorithms: Theory and Analysis

Y. Ye.Interior Point Algorithms: Theory and Analysis. John Wiley & Sons, Ltd, 1997. ISBN 9780471174202. doi: 10.1002/9781118032701. AppendixA.Absorbing sets for the ODE in Section 3.1 We show that the two-dimensional ODE˙x=f ∗(x)withf ∗ given by (3.1) has a compact absorbing set described by a quartic polynomial Lyapunov function, but no absorbing ellipso...

work page doi:10.1002/9781118032701 1997
[64]

We then conclude thatv(x)≥x 2 1 +x 2 2, sovis coercive. To prove thatvsatisfies the Lyapunov inequality (2.4) forα= 1and a suitably largeb, it suffices to show that the leading-order homogeneous part of the polynomial b−v− ⟨f,∇v⟩is coercive. A straightforward computation gives b−v− ⟨f,∇v⟩= 1 3 x6 1 +x 5 1x2 +x 4 1x2 2 +x 2 1x4 2 −x 1x5 2 +x 6 2 +lower-deg...