Global Lyapunov functions: a long-standing open problem in mathematics, with symbolic transformers
Authors: Alberto Alfarano
Francois Charton
Amaury Hayat
Venue: NeurIPS 2024 poster (2024)
Keywords: mathematics
Lyapunov
transformers
control
AI for science
AI for maths
reasoning
Primary area: machine_learning_for_other_sciences_and_fields
TL;DR: Transformers can be trained from synthetic data to find Lyapunov functions, a long-standing open problem in mathematics

Abstract: Despite their spectacular progress, language models still struggle on complex reasoning tasks, such as advanced mathematics.
We consider a long-standing open problem in mathematics: discovering a Lyapunov function that ensures the global stability of a dynamical system. This problem has no known general solution, and algorithmic solvers only exist for some small polynomial systems.
We propose a new method for generating synthetic training samples from random solutions, and show that sequence-to-sequence transformers trained on such datasets perform better than algorithmic solvers and humans on polynomial systems, and can discover new Lyapunov functions for non-polynomial systems.

Global Lyapunov functions: a long-standing open problem in mathematics, with symbolic transformers
Alberto Alfarano"
François Charton"
FAIR, Meta albealfa@meta. com
FAIR, Meta CERMICS, Ecole des Ponts
Echarton@meta..com
Amaury Hayat
CERMICS, Ecole des Ponts Institut Polytechnique de Paris amaury.. hayatSenpc.frr
Abstract
Despite their spectacular progress, language models still struggle on complex rea- soning tasks, such advanced mathematics. We consider long-standing open problem in mathematics: discovering Lyapunov function that ensures the global stability dynamical system. This problem has no known general solution, and algorithmic solvers only exist for some small polynomial systems. We propose a new method for generating synthetic training samples from random solutions, and show that sequence-tor sequence-to-seqeeme to-sequence transformers trained on such datasets perform bet- ter than algorithmic solvers and humans on polynomial systems, and can discover new Lyapunov functions for non-polynomial systems.
Introduction
As large language models achieve human-level performance over broad set of tasks 4 35 45
their capability to reason becomes focus discussion and research. There no single definition of reasoning, and work in this area encompasses factuality, real world alignment, compositionality, the discovery and following of rules, &c. Still, mathematics are considered as one the purest, and most demanding, forms of reasoning [17. As such, solving research-level mathematical problems is a major milestone demonstrating the reasoning capabilities of language models. Such an advance in AI would also transform mathematical practice.
There is little research on applying language models to open problems of mathematics. Except few papers combinatorial optimization and graph theory [34, 39], most prior works focus on problems with known solutions (37 23 g We believe this lack of results is due to two main reasons. First, research problems may require specialized work by mathematicians (6) before they can be handed to language models. Second, most math transformers are trained on sets of problems and solutions which are hard to generate in the case of open problems, when generic method for finding a solution is known.
In this paper, we focus on a long-standing yet easy to formalize, open problem in mathematics: discovering the Lyapunov functions that control the global stability of dynamical systems the boundedness of their solutions when time goes to infinity with respect to an equilibrium or an orbit. A famous instance of this problem the three-body problem: the long-term stability of a system of three celestial bodies subjected to gravitation. The stability problem was studied by Newton, Lagrange and Poincaré. Lyapunov discovered that stability is guaranteed if an entropy-lite function for the system
"Equal contribution
38th Conference on Neural Information Processing Systems (NcurIPS 2024).

---

-the Lyapunov function- can be found. Unfortunately, no method is known for deriving Lyapunov functions in the general case, and Lyapunov functions are only known for a small number of systems.
We propose a new technique for generating training data from randomly sampled Lyapunov functions. Sequence-to-sequene transformers trained these datasets achieve near perfect accuracy (99%) on held-out test sets, and very high performance (73%) on out-of-distribution test sets, We show that higher accuracies (84%) can be achieved by enriching the training set with a small number (300) of easier examples that can be solved with existing algorithmic methods. These enriched models greatly outperform state-of of-the-art techniques and human performance on a variety of benchmarks.
Finally, we test the capability of our models to discover yet unknown Lyapunov functions on randomly generated systems. On polynomial systems, the only ones current methods solve, our models find Lyapunov function for 10.1% or systems, vs 2.1% for mate-of-theear techniques. On non-polynomia systems, where algorithm is known, our best models discover new Lyapunov functions for 12.7% of systems. Our research demonstrates that generative models can be used to solve research-level problems in mathematics, by providing mathematicians with guesses of possible solutions. The solutions proposed by the black-box model are explicit and their mathematical correctness can be verified. We believe this research is an Al-driven blueprint for solving open problems in mathematics.
Related works
Most classical methods for finding Lyapunov rely on parameterized families of candidate solutions, and attempt to derive conditions on the parameters [1] [4]. Additional techniques such as backstep- ping or forwarding II Chap. 12], were introduced to leverage the specifics of particular systems (see also physics-basee methods 444). These techniques are limited to specific, or simple, systems. The global Lyapunov functions of polynomial systems that are sums of squares of polynomials of given degree can be found by computational-intensii algorithmic tools, such as SOSTOOLS 32 33) which leverage the fact that the Lyapunov function belongs to a dimee-dimeniional space.
Methods involving neural networks have been proposed in recent years 7 15 24 They train feed-forward networks to approximate Lyapunov functions of given system, and use a Satisfiability Modulo Theories (SMT) solver as a verifier which proposes potential counter-examples. This approach, very different from ours, was shown be successful for several well-studied high dimensional systems. However, it only finds local or semi-global Lyapunov functions (see Definition A.3)). Since the Lyapunov functions that are found are implicit, it would be hard for mathematicians to check whether they are global Lyapunov functions not. Semi-global Lyapunov functions are useful in many engineering fields such as robotics, where one wants system to be robust to small perturbations. In other fields, like epidemics, being resilient to large perturbations is central, and global Lyapunov functions are required.
Transformers trained on synthetic datasets have been proposed for many problems of mathematics, including arithmetic (28) linear algebra [0)) symbolic integration 22 symbolic regression S
Shortest Vector Problem 40, Gröbner basis computation [19] and theorem proving (30).
investigate problem related to ours: the local stability of dynamical systems. Different architectures
were used to solve hard problems in combinatorial optimisation 34 and graph theory [39].
2 System stability and Lyapunov functions
The stability of dynamical systems is a hard mathematical question, which intrigued many generations of mathematicians, from Newton and Lagrange in the 18th century, to Poincaré in the 20th in the context of the bhree-body problem. The main mathematical tool for assessing stability was proposed by Lyapunov, who showed in 1892 that system is stable if a decreasing entropy- like function the Lyapunov function- can be found 20 26). Later, the existence of a Lyapunov function was shown to be necessary condition for the stability of large classes of systems [2] 18). Unfortunately. these very strong results provide no clue on how to find Lyapunov functions, or just proving their existence for particular system. In fact, 130 years later, systematic derivations of global Lyapunov functions are only known in a few special cases, and their derivation in the general case remains a well-known open problem.
In mathematical terms, we consider the dynamical system
i == f(z),
(I)

---

so = - -7x8 4x0a² - 5.0
in 7xd 321 22²
8 8r9 9
V(x) = 2n0 + 20 + +30g + 2x + +
-10
do =-10 + IDJ1
[1 = =-11
V(x) In(1 + 5 + 2
-30 -10
Figure 2: Two stable systems and associated Lyapunov functions discovered by model. The second, poly- nomial system with a non-polynomial Lyapunov func- tion, was studied in 1.
Figure 1: Dynamic of a stable system: trajectories may be complicated but as long as they start in the red ball they remain the blue ball.
where x E RR" f E C (IR") and dr. We want to know the system has a stable equilibrium around a point x such that f(x*) O. We assume, without loss of generality, that x= 0.
de
Definition 2.1. The system (T) is stable when, for any > O, there exists > such that, [m<(0]]| < 17, the system ① with initial condition (0) has unique solution æ E C²O(O, (20((0+000)) +o0)) and
[=((t] s, AE +00
In other words, a system is stable if solution that begins close to the origin (]){())] < 11) stays close the origin at all time (|]((t)] S =). Lyapunov proved that the stability is related the existence of what is now called a Lyapunov function.
Definition 2.2. The function V E CR""", R.) ) is said to be a (global) Srapanov,functionn for the system 1 if the following condition are satisfied
V(0) = 0, lim V(x) = +00,
||=00->000
(xx) > O, VV(x)- f(z)= S for + O.
Theorem 2.3 (Lyapunov 1892). If the system has Lyapunov function, then it is stable.
In fact, the existence of a Lyapunov function is more powerful and provides additional information. Theorem 2.4 (LaSalle, 1961). If the system (1) has a Lyapunov functiion V, then all the solutions of converge to the largest invariant set of (fT) VV(x) = 03.
In many cases this largest invariant set is reduced to {x* 0} and the system is said globally asymptotically stable (all solutions converge to the equilibrium, see Appendix AA).
Most dynamical systems are unstable. For instance, the solutions of the simple system i(t) = x()) grow exponentially with time, and the solutions of a(t) 1+ x(t)² (x R) always blow up before t = n. No Lyapunov functions can be found for these systems.
On the other hand, stable systems can have an infinite number of Lyapunov functions. The system
(t) (t) (b (1) -x (1)
has ((x) = 2g2g + ardه Lyapunov function for any choice of ag 0 and a1 0.
In the general case, there is no systematic way of discovering a Lyapunov function, or even showing that one exist. Tools exist for small polynomial systems with special "sum of squares" (SOS) Lyapunov functions, but they need a lot of resources, do not always find a solution, and fail once the systems involve more than few variables.
We also consider related, but easier, problem: finding nontrivial V which are semi-definite positive, ie. V verifying V(x) 2 0 instead of Vix) > 0 in Equation B These functions, called barrier

---

functions, form "barriers" that divide R" into two subspaces. A solution starting inside the barrier must remains in the same subspace, which an invariant set of the system [3 4T. For polynomial systems, barrier functions are slightly easier to find using SOS solvers.
3 Experimental settings
In this work, we train eqqencce-o--eeuee given system, when it exists. We frame the problem as a translation task: problems and solutions are represented as sequences of symbolic tokens, and the model is trained from generated pairs of systems and Lyapunov functions to minimize the cross-entropy between the predicted sequence and the correct solution. We train transformers with 8 layers, 10 attention heads and an embedding dimension of 640 (ablation studies on different model sizes can be found in Appendix (C. on batches of 16 examples, using the Adam optimizer [2]] with learning rate of 10-4, an initial linear warm-up phase of 10,000 optimization steps, and inverse square root scheduling. All experiments run on V100 GPU with 32 GB of memory, for 3 or 4 epochs of 2.4 million examples per epoch. Training time is between 12 to 15 hours per GPU.
transformers 38 to predict Lyapunov function for
Tokenization. Model inputs are systems of the form (#; = fi(z1, In) lie (1....)). represented by the n functions f.- Model outputs are single functions (( I.). As in [22]. functions are represented trees, with operators in their internal nodes, and variables or constants their leaves. Trees are then enumerated in Polish (pre-order) notation to produce sequences of tokens that can be processed by the transformer.
All operators and variables are tokenized as single symbols (c.g. 'cos' or 'x1'). "). Integer constants
tokenized sequences of "digits" in base 1000 (e.g. 1024 as the sequence [+, 1, 24]), and real constants, in scientific notation, as pairs of two integers (mantissa and exponent, c.g. -3.14 as [-, (1,314, 10", -; 2)). For instance:
*
sin
2g = cos(2.12)) + 22) ④ sin (3.11 + 2)
cos
is represented as
*
2
2.1 To
enumerated the sequences: [*, cos, *, 2,1, mo; +, 2] and [sin, +, 3, 01 2). and finally tok- enized as |*, cos, *, 21, 10", -, 1, Jp. +, 31, 2, SEP, sin, +, *, 3, æ1, 2) (using SEP as a separator).
Evaluation. Trained models are evaluated on sets of stable systems. Since systems have an infinite number of Lyapunov functions, we cannot check the model predictions by comparing them to the solutions from the test set, and need to there exists a small positive polynomial P such that VV f and V P are sum of squares (SOS) of polynomials (with P 0 for barrier functions), using Python solver based on SumOfSquares [43). For non-poli polynoolll systems, we also use a verifier based on shgo that checks ( numerically. To further ensure correctness we also verify the symbolic solutions using Satisfiability Modulo Theories (SMT) solvers, relying on dReal [3 for verification through interval analysis. This guarantees that equations (3) hold, at least in chosen ball around the origin. The performances of the two verifiers (numerical solver and SMT) are similar, a comparison provided in Table Both the SOS and SMT verifiers sometimes fail return answer. In that case, classify the solution as wrong, even though it might have been correct. As result, model accuracies may be underestimated.
external verifier. For polynomial systems, we verify that
Model predictions use beam search with early stopping, normalizing log-likelihood scores by their sequence length. We report results with beam size (greedy decoding) and beam size 50. With beam size 50, we consider the model to be correct if one Lyapunov function is found among the 50 guesses.
4 Data generation
Our models are trained and tested on large datasets of pairs of stable systems and associated Lyapunov functions. Sampling such stable systems raises two difficulties. First, most dynamical systems are unstable, and general method exists for deciding whether a system is stable. Second, once stable

---

system is sampled, there is no general technique for finding a Lyapunov function, except in particular cases. In this paper, we rely on Backward generation [22]. sampling solutions and generating associated problems, for the general case, and forward generation, sampling systems and calculating their solutions with a solver, for the tractable polynomial systems of small degree.
4.1 Backward generation
Backward generation methods, sampling problems from their solutions, are only useful if the model can be prevented from learning reverse the generation procedure, or from "reading" the solutions in the generated problems. For instance, when training a model to solve the hard problem of finding the roots of an integer polynomial [9] one can easily generate polynomial from its roots, i.e. from the roots 3, 5 and 7, generate the polynomial:
P(X) 2(x + 1)(X - ( - (x
However, the model is trained from factorized form of P(X), ), will learn to read the roots in the problem, instead of computing them. On the other hand, the developed and simplified form
P(X)= 2x 30X* + 144X - 240X + 142X - 210
offers no clues. A second difficulty of backward generation is that sampling solutions instead of problems biases the training distribution. A model trained on backward-generated data may not perform well on a forward-generated test set. Finally, prior work [42] observed that, for hard problems, backward generation methods sometimes focus on easier sub-problems (see, for instance, our comment below about choosing f = VV in step 2).
We propose a procedure for generating stable system S from random Lyapunov function V. The rationale is the following. Since V must be positive with a strict minimum in 0, and tend to infinity at infinity ), we first generate V Vpruper Venus where Vimoper belongs class of functions with a guaranteed strict minimum in zero and Venos to a larger class of non-negative functions, valued at the origin, but with no guarantee of strict minimum (step and Appendix [B. From V. need to generate f so that the third condition of (3) is met. A naive solution would be f = VV since f. VV S 0 would hold. But this would severely limit the systems we create, and turn the Lyapunov function discovery problem (find V from f) into an easier integration problem (find V from - VV). Instead, starting from fo = VV.. we apply the following transformations:
multiply each coordinate of fo by random non-negative functions h? (step 4) and call it fo- generate random function & = E 9:(z)) (i)) (steps and 3), where et are orthogonal
to (zx), and set y fo- We have db VV = 0 and (Φ + fo) VV S.
These transformations guarantee that all conditions in 3 met. On the other hand, they allow f to span a very large set of systems, since any I satisfying (x=) f(x) < 0 can be written as the sum of function collinear to (z) and function orthogonal Vz
Specifically, the procedure can be summarized as follows (see Appendix for more details).
Step 1 Generate random function V, satisfying V(x) > V(0), Vr E R" \ {0). and
Vzz) +oc when ||>||| - +00.
Step 2 Compute the gradient VV(x) and denote Hy = {2 E R" 2 VV(x) = 0} the hyperpland eyyerplannfoorhhooonaal to (V(x), for any x E R.
orthhggona
Step 3 Select S S at random and sample vectors (e(())) from hyperplane
H.g. Generate real- functions (g))ie(1....)))
Step Select < < n at random, generate bn random real-valued functions (ha))e(t......)) set h, = 0 for ka +1 +lsisn n.
Step 5 Build the n functions
f(x) = -
+ EZ
with a random permutation of [1....,
Step 6 Simplify the functions fi, obscuring patterns from the generative process.
ir (z) = 0 this is the whole space instead, but this does not change the method.
5

---

SMT Solver
SOS Solver
Timeout
10 minutes
60 minutes
minutes
60 minutes
Correct Lyap function
82.6
94.1
89.6
95.3
Solver Timcouts
17.4
5.9
10.4
4.7
Incorrect Lyap function
0
0
0
Table 1: SMT and SOS timeout and error rates, benchmarked on correct Lyapunov functions.
This method produces a stable system S à f()), with V as its Lyapunov function. The difficulty of inferring V from S hinges on a careful choice of the vectors e'. For instance, if we naively select et as an orthonormal basis of Hz. computed from (z(x) by Craam-chmiitt orthogonalization, prefactors like 1/IlVV((x) [xr)] appear at step 3, and are unlikely to simplify away step 6. This provides the model with a shorteut: reading [>)]l in S, and using it to recover VV and then V, not a trivial task, but easier one than discovering Lyapunov functions. To counter this, we relax the orthonormality condition on e² (æ), so that 1/ (>))ll never appears, yet keep the et e'(x) simple enough for VV- -sppciificc patterns in E; (x))e (z) simplify away at step 6. We also want to ensure that the e span all of Ha. or the systems generated will not be diverse enough.
In our experiments, we slightly modify this procedure, by running steps Lyapunov function V created at step 1. As result, 5 systems are generated that share the same Lyapunov function discussion of this choice can be found in Appendix C.J). From a mathematical point of view, Lyapunov function describes a hidden quantity in a system, and we believe that providing the model with several systems that share this hidden quantity should help it learn the parts of the system that contribute to this hidden quantity, and therefore learn a I Lyapunov function.
five times for each
This procedure
be tuned to generate specific classes of systems. By choosing V. Di and h, in particular classes, we can constrain the system functions fi to be polynomials, polynomials of functions (e.g. trigonometrie polynomials), or more general functions (see Appendix B.4 for more).
The Lyapunov functions obtained here are correct by design. Nevertheless, we still performed an evaluation of the solutions both safeguard and to benchmark the failure and timeout rates of the SMT and SOS solvers on correct solutions, which we report in Table
4.2 Forward generation
Whereas the stability problem is unsolved in the general case, methods exist to calculate Lyapunoy functions of polynomial systems, when they exist and can be written als (see Section D]. These algorithms, of polynomial complexity, are very efficient for small systems, but their CPU and memory requirements explode as the size of the systems grows. We leverage them to generate forward datasets, as follows.
of squares of polynomi-
Step 1 Generate polynomial system at random
Step 2 Use a routine to find a polynomial sum-of-squares (SOS) Lyapunov function. Step 3 Keep the system if such function exists, restart from step I otherwise.
This approach has several limitations. First, since most polynomial systems are not stable, and the computation of SOS Lyapunov function involves complicated search K3 it is slow and limited to small systems of polynomials with small degree. Second, because not all stable polynomial systems have polynomial SOS Lyapunov functions [1, it can only generate subset of stable polynomial
systems.
Finally, SOS routines process the constraints in Equation 3 by solving semi-definite programming (SDP) problems. This guarantees that V is sum-of-squares, hence we have V(x) > 0. but not necessarily Vxx) > 0, for * 0. As a result, these methods can only discover barrier functions. State-of-the the art methods circumvent this by introducing the stronger constraint V(x) 2 >i_1 e with E; small [32]. V then has a unique minimum in x 0, which makes it Lyapunov function, but this further restricts the class of polynomial systems that the method can solve.
4.3 Datasets
We generate backward and forward datasets for training and evaluation purpose, and one smaller forward dataset for evaluation purposes (see Table in Appendix B.6 for list).

---

Backward datasets Our main backward set, BPoly, features 1 million non -oneeneeeeetee polynomial systems S with integer coefficients, and to 5 equations (in equal proportions). We also create BNonPoly, dataset of 1 million non-degenerate non-polynomial systems with 2 to 5 equations. In this dataset, the coordinates of f are polynomials of general functions, e.g. trigonometrie polynomials, or functions such as 3 cos( + e= For such general systems, no method for discovering Lyapunov function is known.
Forward datasets All 2 forward datasets are generated using solver derived from the SumOf- Squares package in Python, and implementing techniques similar those used in SOSTOOLS (see Appendix B.5) All systems in these datasets are non con-zero integer polynomials with 2 to 3 equations, and integer polynomial Lyapunov functions the only systems these methods can solve. We create FLyap, dataset of 100,000 systems having a non-homogeneous polynomial as a Lyapunov function. We also have dataset focusing on barrier functions (see the end of section [.2]p FBarr features 300,000 systems having a non-homogeneous polynomial as a barrier function. The small size of these datasets due the computational cost of SOS methods, and the difficulty of discovering Lyapunov barrier functions.
To allow for comparison with SOSTOOL, the state-of-the-art method for discovering Lyapunov functions of polynomial systems, we also generated a test set of 1,500 polynomial systems with integer coefficients that SOSTOOLS can solve (FSOSTOOLS).
5 Results
Our models trained on different datasets achieve near perfect accuracy on held-out test sets, and very high performances on dfiss--dittiibutiin test sets, especially when enriching the training set with small number of forward examples. They greatly outperform ttate-of-the-art techniques and also allow to discover Lyapunov functions for new systems. These results are detailed below.
5.1 In and out-of-distribution accuracy
In this section, we present the performance of models trained on the 4 datasets. All models achieve high in-domain accuracy when tested on held-out test sets from the datasets they were trained on (Table On the forward datasets, barrier functions are predicted with more than 90% accuracy, and Lyapunov functions with more than 80%. On backward datasets, models trained on BPoly achieve close to 100% accuracy. We note that beam search, i.e. allowing several guesses at the solution, brings significant increase in performance (7 to 10% with beam size 50, for the low-performing models). We use beam size 50 in all further experiments.
Accuracy
Accuracy
Backward datasets
bs=1
bs=50
Forward datasets
bs=1
bs=50
BPoly (polynomial)
100
FBarr (barrier)
93
98
BNonPoly (non-poly)
77
87
FLyap (Lyapunov)
81
88
Table 2: In-domain accuracy of models. Beam size (bs) I and 50.
The litmus test for models trained generated data is their ability to generalize nut-of-distribution (OOD). 3resentss evaluations of backward models on forward-generated sets (and the other way around). All backward models achieve high accuracy (73 to 75%) when tested on forward- generated random polynomial systems with sum-of-squares Lyapunov functions (FLyap). The best performances are achieved by non-polynon pol-oolyloiial systems (BNonPoly), the most diverse training set. The lower accuracy of backward models on forward- generatde sets of systems with barrier functions (FBarr) may be due to the fact that many barrier functions are not necessarily Lyapunov functions. On those test sets, backward models must cope with a different distribution and (slightly) different task. Forward models, the other hand, achieve low performance on backward test sets. This is possibly due to the small size these training set.
Overall, these results seem to confirm that backward-traine models are not learning to invert their generative procedure. If were the case, their performance on the forward test sets would be close to zero. They also display good OOD accuracy.
7

---

Backward datasets
FLyap
FBarr
Forward datasets
BPoly
BPoly (polynomial)
73
35
FBarr (barrier)
15
Bonnoll (non-poly)
75
24
FLyap (Lyapunov)
10
Table 3: Out-of-domain accuracy of models. Beam size 50. Columns are the test sets.
5.2
Enriching training distributions for improved performance
To improve the OOD performance of backward models, we add to their training set a tiny number of forward-generated examples, as [66]. Interestingly, this brings a significant increase performance 4). Adding 300 examples from FBarr to BPoly brings accuracy on FBarr from 35 to 89% (even though the proportion of forward examples in the training set only 0.03%) and increases OOD accuracy on FLyap by more than 10 points. Adding examples from FLyap brings less improvement.
These results indicate that the OOD performance of models trained backward-generated data can be greatly improved by adding to the training set a small number of examples (tens or hundreds) that we know how solve. Here, the additional examples solve weaker but related problem: discovering barrier functions, The small number of examples needed to boost performance makes this technique especially cost-effective.
Forward datasets
Examples added (IM training set)
Flyap
FBarr
No mixture
0
73
35
FBar
30
75
61
300
83
89
3,000
85
93
30,000
89
95
FLyap
10
75
25
100
82
29
1,000
83
37
10,000
86
Table 4: Mixing backward data (BPoly) with a small number of forward examples. Beam size 50.
5.3 Comparing with state-of-the-art baselines
To provide a baseline for our models, we developed f findlyap, Python counterpart to the MATLAB Lyapunov function finder from SOSTOOLS (see Appendix B.5). We also introduce FSOSTOOLS.
test set of 1,500 polynomial systems with integer coefficients that SOSTOOLS can solve. We also tested Al-based tools (see Appendix E for the full list parameters sweeps we used for each of these methods), such as Fossil 2 [12). ANLC v2 [15) and LyzNet (25). These methods achieve low accuracies on our test sets. This might be due to the fact that these tools are designed solve different problem: discovering local or semi-globar Lyapunov function (and potentially finding control function), while we target global Lyapunov functions.
SOSTOOL findlyap
Existing Al methods
Models
Test sets
Fossil
ANLC
LyzNet
PolyMixture
FBar
FLyap
BPoly
FSOSTOOLS
32
30
46
84
80
53
54
FBam
12
89
FLyap
42
32
66
83
93
73
BPoly
15
10
24
99
15
10
Table 5: Performance comparison different test sets. Beam size 50. PolyMixture BPoly + 300 FBam.
Table 5 compares findlyap and Al-based tools to our models on available test sets. A model trained on BPoly complemented with 500 systems from FBarr (PolyMixture) achieves 84% on FSOS- TOOLS, confirming the high OOD accuracy of mixture models. On generated test sets, PolyMix- ture achieves accuracies over 84% whereas f findlyap achieves 15% the backward generated test

---

set. This demonstrates that, on polynomial systems, transformers trained from backward-generated data achieve very strong results compared the previous state of the art.
On average Transformer-hase models are also much faster than SOS methods. When trying to solve random polynomial system with to 5 equations (as used in Section findlyap takes an average of 935.2s (with timeout of 2400s). For our models, inference and verification of one system takes 2.6s on average with greedy decoding, and 13.9s with beam size 50.
5.4 Into the wild discovering new mathematics
Our ultimate goal to discover new Lyapunov functions. To test our models' ability to do so, we generate three datasets of random systems: polynomials systems with 2 or 3 equations (Poly3), polynomial systems with to 5 equations (Poly5), and non-polynomial systems with 2 or 3 equations (NonPoly). For each dataset, we generate 100,000 random systems and eliminate those that are trivially locally exponentially unstable in x* 0, because the Jacobian the system has an eigenvalue with strictly positive real part (20]. We compare f findlyap and AI based methods with two models trained on polynomial systems, FBarr, and PolyM(xture) a mixture of BPoly and 300 examples from FBarr- and one model trained on a mixture of BPoly, BNonPoly and 300 examples from FBarr (NonPolyM).
Table 6 presents the percentage of correct solutions found by our models. On the polynomial datasets, our best model (PolyM) discover Lyapunov functions for 11.8 .8 and 10.1% of the (degree 3 and degree 5) systems, ten times more than findlyap. For non-polynomial systems, Lyapunov functions are found for 12. 7.7% of examples. These results demonstrate that language model trained from generated datasets of systems and Lyapunov function can indeed discover yet unknown Lyapunov functions and perform at a much higher level that state-of-the-an SOS solvers.
Sample
SOSTOOL findlyap
Existing Al methods ANLC
Forward FBarr
Backward models
Test set size
Fossil 2
LyzNet
PolyM
NonPolyM
Poly3 65,215
1.1
0.9
0.6
4.3
11.7
11.8
11.2
Poly5 60,412
0.7
0.3
2.1
10.1
NonPoly 19,746
1.0
0.6
3.5
12.7
Table 6: Discovering Lyapunov comparison for random systems. Beam size 50. PolyM is BPoly + 300
FBarr NonPolyM is BNonPoly BPoly + 300 FBamr
5.5
Expert iteration
Given the performance on our model in Table we can use the newly solved problems to further fine-tune the model. Specifically, we create sample of verified model predictions for polynomial systems, WindoTheWild, we add to the original training sample and we continue training the model.
We test different strategy to finetune the model and we report performance on forward benchmarks and "into the wild" in Table
n/: Add 20,600 samples from BPoly (20,000), FBarr (50), FLyap (50) and FIntoTThe Will (500) to keep similar proportion used during pretraining
n2: Add 2,000 samples from FLyap (1,000) and WintoTheWild (1,000) to improve on both forward benchmark and in the wild
n3: Add 50 samples from FInto Will to show that this indeed helps
n4: Add 1,000 samples from Whelol
n5: Add 2,000 samples from FintoTheWWld
nб: Add 5,000 samples from FintoThe Wild to see if there are benefits to add more samples
We also retrain a model (n7) from scratch using a mixture of BPoly (IM), FBarr (500), FLyap (500) and FhetoTheWild (2,000).
We notice that the addition of 1,000 verified predictions our training set of million improves performance on the "into to wild" test sets by about 15%, while not affecting the other test sets (n4). Adding more examples seems to be detrimental, as it decreases the performance on other benchmarks (n5 and nō). We also notice that finetuning with mixed data from other distributions is not efficient (n/ and n2) and a small contribution already help to get some improvements (result n3). Finally, it's not efficient to pretrain the model from scratch using data from Finto (n?).

---

Forward benchmark
Regenerated Thel Poly3
Strategy
FBarr
FLyap
Poly5
Baseline
93
84
11.7
9.6
nl
94
84
10.3
9.6
n2
90
12.2
11.3
n3
92
84
12.4
10.1
n4
92
84
13.5
11.9
n5
89
79
13.5
11.9
n6
85
72
13.5
11.9
n7
93
12.1
10.0
Table 7: Expert iteration using Into The Wild correct guesses. The Poly3 and Poly5 test sets are regenerated, prevent data contamination.
6
Discussion
We have shown that models can be trained from generated datasets to solve long-standing open problem in mathematics: discovering the Lyapunov functions of stable dynamical systems. For random polynomial systems, our best models can discover Lyapunov functions in five times more cases than state-of-the-art methods. They can also discover Lyapunov functions of non-polynomial systems, for which no algorithm is yet known, and were able to re-discover non non-polynomial Lyapunov function of a polynomial systems discovered by [[] (Appendix (F).
The backward generation method introduced in section is the key innovation in this paper. The main problem with such approaches is their tendency generate training sets with very specific distributions, which prevent models from generalizing to general instances of the problem. Our models can generalize out of their training distributions (Table and can improve their performance by adding to their training set a tiny number of systems that we know how to solve (Table
While our models exceed the algorithmic state of the art, one might wonder how they compare to human mathematicians. To this effect, we proposed 75 problems from the FSOSTOOLS dataset (polynomial systems with 2 or 3 equations) as an examination for 25 first year Masters students in mathematics, following a course on the subject. Each student was given systems chosen at random and had a total of 30 min. Their performance was 9.33%, significantly lower than our models (84%).
Our work has a number of limitations. Because there is no known way to tell whether a random system is stable, we lack good benchmark on non- non-polynomial systems. Also, all the systems studied in this paper are relatively small, at most 5 equations for polynomial systems and 3 for non-polynomial. We believe that scaling to larger models should help tackle larger, and more complex, systems. Finally, this work could extended to take into account the domain of definition of non-polynomial systems.
The broader implications of our work extend into two directions: the capability of transformers to reason, and the potential role of AI in scientific discovery. While large language models perform at human level on a broad set of tasks, they are still embarrassingly clumsy on many simple problems of logic and reasoning, may be an inherent limitation of auto-regressive transformer architectures. Our results suggest that transformers can indeed be trained to discover solutions to hard problem of symbolic mathematics that humans solve through reasoning, and that this is enabled by careful selection of training examples, instead of a change of architecture. We do not claim that the Transformer reasoning but it may instead solve the problem by kind of "super-intuition" that stems from a deep understanding of a mathematical problem.
the point that was suggested that planning and high level reasoning
From a mathematical point of view, we propose new, Al-based, procedure for finding Lyapunov functions, for broader class of systems than were previously solvable using current mathematical theories. While this systematic procedure remains a black box, and the "thought process" of the transformer cannot be elucidated, the solutions are explicit and their mathematical correctness can be verified. This suggests that generative models can already be used to solve research-level problems in mathematics, by providing mathematicians with guesses possible solutions. While small minority of mathematicians is currently using deep-learning tools, we believe generative models have the potential to foster tremendous progress on number of research subjects, and may eventually become a central component in the future landscape of mathematical practice.
10

---

Acknowledgements
This work was performed in part using HPC resources from GENCI-IDRIS (Grant 2023- AD01 101452).... The authors also acknowledge the Office of Advanced Research Computing (OARC) Rutgers. The State University of New Jersey. The authors would also like to thank the Master students of the Mathematics and Computer Science Department of the Fcole des Ponts IP Paris from the year 2023-2024 who attended the course Control of Dynamical Systems. The authors also wish to thank Guillaume Lample for fruitful discussion.
References
[1] Amir Ali Ahmadi, Miroslav Kistic, and Pablo A. Parrilo. A globally asymptotically stable polynomial vector field with no polynomial Iyapunov function. In 2017 50th IEFE Conference on Decision and Control and European Control Conference, pages 7579-7580, 2011.
[2] Amir Ali Ahmadi, Miroslav Krstic, and Pablo A. Parrilo. globally asymptotically stable polynomial vector field with no polynomial Iyарunov function. In 2011 50th IEEE Conference on Decision and Control and European Control Conference, pages 7579-7580, 2011.
[3] Martin S. Andersen, Joachim Dahl, and Lieven Vandenberghe. Cvxopt. https:// //cxppt org/ 2023. Version 1.3.2; Free software package for convex optimization based Python. Last updated August 9, 2023.
[4] Anton Bakhtin, Noam Brown, Emily Dinan, Gabriele Farina, Colin Flaherty, Daniel Fried, Andrew Goff, Jonathan Gray, Hengyuan Hu, Athul Paul Jacob, Mojtaba Komeili, Karthik Konath, Minae Kwon, Adam Lerer, Mike Lewis, Alexander H. Miller, Sasha Mitts, Adithya Renduchintala, Stephen Roller, Dirk Rowe, Weiyan Shi, Joe Spisak, Alexander Wei, David Wu, Hugh Zhang, and Markus Zijlstra. Human-level play in the game of <i>diplomacy</i> by combining language models with strategic reasoning. Science, 37866624)10077-17777 2022.
[5] Luca Biggio, Tommaso Bendinelli, Alexander Neitz, Aurelien Lucchi, and Giambattisa Paras- candolo. Neural symbolic regression that scales, 2021.
[6] Kevin Buzzard, Johan Commelin, and Patrick Massot. Formalising perfectoid spaces. In Proceedings of the 9th ACM SIGPLAN International Conference on Certified Programs and Proofs. ACM, jan 2020.
[7] Ya-Chien Chang, Nima Roohi, and Sicun Gao. Neural Lyapunov Control. In Advances in Neural Information Processing Systems, volume 32. Curran Associates, Inc., 2019.
[8] François Charton, Amaury Hayat, and Guillaume Lample. Learning advanced mathematical computations from examples. arXiv preprint arXiv:2006. 204::2000000662,, 2020.
[9] François Charton. Computing the roots of polynomials, 2022. https: chharton.githl htppooyyomiall-rooe
[10] François Charton. Linear algebra with transformers, 2022.
[11] Jean-Michel Coron. Control and nonlinearity. American Mathematical Soc., 2007.
[12] Alec Edwards, Andrea Peruffo, and Alessandro Abate. Fossil 2.0: Formal certificate synthesis for the verification and control of dynamical models. In Proceedings of the 27th ACM Interna- tional Conference on Hybrid Systems: Computation and Control, HSCC '24, New York, NY, USA, 2024. Association for Computing Machinery.
[13] Sicun Gao, Soonho Kong, and Edmund M. Clarke. dreal: An smt solver for nonlinear theories over the reals. In Maria Paola Bonacina, editor, Automated Deduction - CADE-24, pages Berlin, Heidelberg, 2013. Springer Berlin Heidelberg.
[14] Peter Giesl. Construction of global Lyapunov functions using radial basis functions, volume 1904. Springer, 2007.
[15] Davide Grande, Andrea Peruffo, Enrico Anderlini, and Georgios Salavasidis. Augmented Neural Lyapunov Control. IEEE Access. 11:67979-67986, 2023.

---

[16] Samy Jelassi, Stéphane d' Ascoli, Carles Domingo-Emnnn Enrich, Yuhuai Wu, Yuanzhi Li, and François Charton. Length generalization in arithmetic transformers, arXiv preprint (Dt::3066114000 2023.
[17] Immanuel Kant. Critique of pure reason. 1787.
[18] Christopher M Kellett, Classical converse theorems in lyapunov's second method. Discrete & Continuous Dynamical Systems-Series B, 20(8), 2015.
[19] Hiroshi Kera, Yuki Ishihara, Yuta Kambe, Tristan Vaccon, and Kazuhiro Yokoyama. Learning to compute grl"" obner bases. arXiv preprint arXiv:2311.2290, 2023.
[20] Hassan K. Khalil. Nonlinear systems. Macmillan Publishing Company, New York, 1992.
[21] Diederik P Kingma and Jimmy Ba. Adam: A method for stochastic optimization. arXiv preprint arXiv:1412.6980, 2014.
[22] Guillaume Lample and François Charton. Deep learning for symbolic mathematics. arXiv preprint 192.0411 2019.
[23] Guillaume Lample, Marie- Anne Lachaux, Thibaut Lavril, Xavier Martinet, Amaury Hayat, Gabriel Ebner, Aurélien Rodriguez, and Timothée Lacroix. Hyper'Tree Proof Search for Neural Theorem Proving. Advances in neural information processing systems, 2022.
[24] Jun Liu, Maxwell Fitzsimmons, Ruikun Zhou, and Yiming Meng. Formally verified physics- informed neural control Yyapunov functions. arXiv preprint 2xXi::2409.20528, 2024.
[25] Jun Liu, Yiming Meng, Maxwell Fitzsimmons, and Ruikun Zhou. Lyznet: A lightweight python tool for learning and verifying neural іyapunov functions and regions of attraction, 2024.
[26] Aleksandr Mikhailovich Lyapunov. The general problem of the stability of motion. PhD thesis, Kharkov University, 1892.
[27] Jose Luis Massera. On liapounoff's conditions of stability. Annals of Mathematics, pages 705-721, 1949.
[28] Rodrigo Frassetto Nogueira, Zhiying Jiang, and Jimmy Lin. Investigating the limitations of the transformers with simple arithmetic tasks. CoRR, 0os/2102.13011 2021.
[29] KP Persidskii. On a theorem of liapunov. In CR (Dokl.) Acad. Sci. URSS, volume 14, pages 541-543, 1937.
[30] Stanislas Polu and Ilya Sutskever. Generative language modeling for automated theorem proving, 2020.
[31] Stephen Prajna, Ali Jadbabaie, and George J Pappas. A framework for worst-case and stochas- tic safety verification using barrier certificates. IEEE Transactions on Automatic Control, 52(8): 2007.
[32] Stephen Prajna, Antonis Papachristodoullou and Pablo A Parrilo. Introducing sostools: A gen- eral purpose sum squares programming solver. In Proceedings of the 41st IEEE Conference on Decision and Control, 2002., volume 1, pages 741-746. IEEE, 2002.
[33] Stephen Prajna, Antonis Papachristodoulou, Peter Seiler, and Pablo A Parrilo. Sostools and its control applications. Positive polynomials in control, pages 273-292, 2005.
[34] Bernardino Romera-Paredes Mohammadamin Barekatain, Alexander Novikov, Matej Balog, M. Pawan Kumar, Emilien Dupont, Francisco J. R. Ruiz, Jordan S. Ellenberg, Pengming Wang, Omar Fawzi, Pushmeet Kohli, and Alhussein Fawzi. Mathematical discoveries from program search with large language models. Nature, 625(7995):468-444 2024.
12

---

[35] Baptiste Rozière, Jonas Gehring, Fabian Gloeckle, Sten Sootla, Itai Gat, Xiaoqing Ellen Tan, Yossi Adi, Jingyu Liu, Romain Sauvestre. Tal Remez, Jérémy Rapin, Artyom Kozhevnikov, Ivan Evtimov, Joanna Bitton, Manish Bhatt, Cristian Canton Ferrer, Aaron Grattafiori, Wenhan Xiong, Alexandre Défosse7, Jade Copet, Faisal Azhar, Hugo Touvron, Louis Martin, Nicolas Usunier, Thomas Scialom, and Gabriel Synnaeve. Code Ilama: Open foundation models for code, 2024.
[36] Guillaume Sagnol and Maximilian Stahlberg. PICOS: A Python interface to conic optimization solvers. Journal of Open Source Software, 7(70):3915, February 2022.
[37] Trieu H. Trinh, Yuhuai Wu, Quoc V. Le, He He, and Thang Luong. Solving olympiad geometry without human demonstrations. Nature, 62777995):47--88, 2024.
[38] Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones, Aidan N Gomez, Lukasz Kaiser, and Illia Polosukhin. Attention is all you need. In Advances neural information processing systems, pages
2017.
[39] Adam Zsolt Wagner. Constructions in combinatorics via neural networks. arXiv preprint arXiv:2104. 14516, 2021.
[40] Emily Wenger, Mingjic Chen, Francois Charton, and Kristin E Lauter. Salsa: Attacking lattice cryptography with transformers. Advances in Neural Information Processing Systems, 35:34981-34999 2022.
[41] Xiangru Xu, Paulo Tabuada, Jessy W Grizzle, and Aaron D Ames. Robustness of control barrier functions for safety critical control. PaAC-PapersOnLine, 48(27))54-61, 2015.
[42] Gal Yehuda, Moshe Gabel, and Assaf Schuster. It's what machines can learn, it's what cannot teach. arXiv preprint arXiv:2002..99398 2020.
[43] Chenyang Yuan. Sumofsquares.py. 2024. Software version 1.2.0; Available at http://giiilll.ooomyymmeeennnnns
[44] Ruo-Shi Yuan, Yi-An Ma, Bo Yuan, and Ping Ao. Lyapunov function as potential function: A dynamical equivalence. Chinese Physics B, 23(1):010505, 2013.
[45] Yongehao Zhou, Andrei Ioan Muresanu, Ziwen Han, Keiran Paster, Silviu Pitis, Harris Chan, and Jimmy Ba. Large language models are human-level prompt engineers, 2023.
13

---

Appendix
A Mathematical definitions
In this Appendix, we recall several mathematical definitions and theorems related to the Lyapunov function problem. We first introduce the notion of global asymptotic stability (GAS).
Definition A.I. We say that the (equilibrium x² = 0 of the) system ① is globally asymptotically stable if it is stable and for any zo E IR" there exists unique solution z E +o0); RT) to which satisfies in addition
lim xx)) 0. t-+
This notion translates the fact that the equilibrium x* is robust even to large perturbations. This notion related to the existence of a Lyapunov function thanks, for instance, to LaSalle Invariance Principle:
Theorem A.2 (LaSalle Invariance Principle (global)), Assume there exists a Lyapunov function for the system ① and let S be the largest subset of { (VV= }} that is invariant by the dynamics of ①1 IfS = {0), then the system T globally asymptotically stable.
Note that if (x) f(z) Lyapunov function is challenging mathematical problem, and still an open problem in general, weaker notions exists.
for any I # 0 then necessarily S {0) Because finding a (global)
Definition A.3. The function V E C² (IR", RR,. is said to be a semi-global Lyapunov function for the system ① if there exists such that the following condition are satisfied
((0) 0, V(x) > O, (v(a) f(x) S 0 for | S
(5)
Finding a semi-global Lyapunov function is usually easier than finding global Lyapunov function. A semi-global Lyapunov function is enough to show that the equilibrium x" 0 is robust to small perturbations which, for several engineering applications, is enough. More specifically.
Definition A.4. We say that the (equilibrium x² 0 of the) system ( locally asymptotically
stable if it is stable and if there exists r > 0 such that for any ||z0||| S r there exists unique solution x E ([0, +oo); IR*) to (1) which satisfies addition
lim x(t)= 0.
(6)
Similarly to global Lyapunov function, the existence a semi-global Lyapunov function useful to ensure local asymptotic stability
Theorem A.5 (LaSalle Invariance Principle (local)). Assume there exists semi-global Lyapunov function V, and let S be the largest subset of (vV(x) f(x) = 0} invariant by the dynamics of (1.
If {0} then the system T is locally asymptotically stable.
B Generation procedure
B.1 Function generation
To generate random functions sample random trees with unary and binary internal nodes, and then randomly select operators for these nodes, and variables and integers for leaves (as in [2 g). Our binary operators are the four operations and the power function. Unary operators exp, log, sqrt, sin, cos, tan.
To generate polynomials, we randomly sample given number of monomials, with integer real coefficients. The number of monomials, range of the coefficients, and the powers and number of terms of each monomial, are randomly selected between bounds, provided as hyperparameters.
14

---

B.2 Backward generation
We build globally stable systems by first generating Lyapunov function V at random, and then building a dynamic system which has V as a Lyapunov function. The procedure is:
Step la: We generate V V Venoss + Vpoper where Vmnger belongs to given class of positive defi- nite function and Veess belongs to a larger class, but of non-negative functions only with Veress (0) = 0. More specifically, generate
Veross (I)
Σ (x)),
(7)
with m a random integer, and p; random functions verifying (o( = The nature of functions
P4 depends the systems we want generate (polynomial not). Clearly Veross (ar) 2 0 and Veross (0) = O. We similarly generate
n
Veroper (z)
Σ on a 1
(8)
with n a random integer, 8; random positive integers and A (oos.) definite matrix, with given probability of being diagonal. As consequence, Vproper is strictly minimal in x = 0. When generating barrier functions, we can optionally set Ypooeer = O.
random positive
Step 1b: In this step. increase the class of functions that be sampled for Veross and Viroper by several transformations:
1. Composition of Vrroper with probability Pi,c, replace
Veroper (x) + I(V (Vproper (x))
(9)
with selected at random from pre-defined set of increasing (Appendix B.4).
2. Product 1 propper with probability Pi,m. replace
Tnoper(1) + (Vprope((x)) (Onope(0)))((()))))
(10)
with g selected random from a pre-defined set of positive-functions (Appendix B.4).
Botts Boiss for q n and h sub-expression of Vpopeer i.e. h(x)) Σ (0) Fo(s)
permutation of nI.....n).
3. Composition of V Venas: for every E {I,..., m). with probability P2, replace (x + by + pt())),
(11)
with real function that is bounded from below with minimum (not necessarily unique) in & and chosen at random from pre-defined set of bounded-functiins (Appendix B.4). Recall that p; are the functions appearing in Veross.
Step 1e: Gathering the functions Vproper and Venas together, we define the Lyapunov function (candidate) V(x) V1noss + V Proppe(I).. Overall, we have
V(x) /
Σ a
S uz T(0) g Σ 100
Bos Bott) o()) a(j)
m
be + Pp(x)))
where 7 is the identity with probability 1 P1,e. is the constant function I with probability 1 P1.m and by (x) = x² with probability 1 - P2. Such Lyapunov function satisfies
V(x) V(0), Ve T \ (0).
(12)
Indeed,
m
V(0)= Σ (E())
and V(x)= > be(EE) for any x E TR \ (0), since g positive function, increasing and (JE...) positive definite.
15

---

Step 2: In this step we create the random vectors orthogonal to VV that will be useful in the generation of the system f (see Section 4.1 Taking advantage of the form the condition (3 for any E R, denote
Hl = E Rd |
(v) = =0
the hyperplane orthogonal to (z) Then, for a random p E (1....,] n). generate p random real-valued functions (I...e) and p vectors (e* (E(I.....)) from this hyperplane follows:
Ax if ri) (,) j= T2(()
(13)
otherwise,
where A (V()) and A, refers to the j th component of the vector and Ti and 72 random functions from (0) into (1,...,) such that (1) r2().. This implies that for any E (1,...,)
VV(x) - (VV(z)) (VV())))
(VV(z)) (VV = 0. (14)
Note that, so long (v o O, one can use this process to construct generative family of Hz, and the e¹ span the whole Hz. If VV(x) then Hy =K".
I."."
Step 3: Generate at random k real-valued functions (h+)sE(1....d)) where 1 S k1 S n chosen at random. Set h; 0 for hn <i
Step 4: Build the system
p
f(x) = (1
kE (e......))) +
Σ
(()))e((r),
(15)
with æ a random permutation of [1,...,n..
Overall, the function f satisfies
VV(x)-((0) - (()))-())--33 a - 50,
(16)
hence V is a Lyapunov function of the system a(t) = f(x(t))
Step 5: Expand and simplify the equations of (using Sympy), in order to eliminate obvious patterns due to the generation steps (that the model could recognize and leverage), eliminate duplicated systems in the training set, and limit the length of training sequences. All polynomial systems are expanded into normal form.
B.3 Backward generation modes
generate polynomial systems with sum-of-squarr Lyapunov functions to allow for easy comparison with existing methods such as SOSTOOLS polynomials with zero-order term and Pı,c Pi.m P2 0. Also, fo and zi are polynomials (Appendix I. We generate fi with a degree lower or equal to half the maximal degree of and maximal value of coefficients of the order of the square root of the maximal value of gi- Since the ft are squared in the final system, this allows f and to have the same order, and prevents the transformer from inferring unwanted additional information by looking at the higher degree monomial.
Polynomial generation:
this case, all P,
Generic generation: P. is generated as P.(x) Q,(x) (0), where Q((x) is random function generated as per Appendix and f; and 9s are also generated as per Appendix B.T Optionally the functions can be generated polynomials of non-polynomial operators taken from a dre-defined set of operators.
Other generation modes: we have other generation modes corresponding to interesting particular cases: gradient flow systems, systems where the 2-norm (resp. a weighted 2-norm) is Lyapunov function, etc.
16

---

B.4 Generation design parameters
Our generator allows us to generate generic stable systems and yet to have a large control on the distribution. For polynomials, for instance, we have control on the maximal and average degree, number of monomials, power and number variables the monomials, coefficients, etc. We can also specify whether the coefficients are integers, floats, with which precision. Overall we have total of 36 generation hyper- seeer--aameters that influence the distribution of the synthetic data created. The main generation design parameters are:
base: encoding base for integers
max int: Maximum integer value
precision: Float numbers precision
prob int: Probability of sampling integers variables (for non con-polynomial expressions)
min dim: minimal number of equations the system
max ___diimmoaa maximal number of equations
max max_degrre: maximal degree polynomial terms in a Lyapunov function
n t_terms: maximal number of terms in polynomials for the Lyapunov function
nb ab_ops_proper: maximal number of operators
Vproper (non polynomial generation)
nb me__op__lapp maximal number of operators in V Pproper (non polynomial generation)
operators_ Jyyp: list of operators to be considered (non polynomial generation)
polynomial V: if true generated expressions are polynomials of (potentially oon-polynomia) operators
pure_ pateeepoooooalll... generate polynomial systems only
cross term: Veross if False.
max_ni nbb _e___nb_roos____________________ cllss_roosttreeee term: bound in Veross
mooba_diagonal with this probability, the positive definite form of Yproper is imposed to be diagonal
only 1oy_22_ncc 2_norm: if True, the Lyapunov function is the 2-norm.
strict: True, generates strict Lyapunov function (i.e. VV << 0)
proper: if set to false, Vproper 0 and V is only a barrier function.
float resolution cicolutionppoly::: poly: float resolution of the polynomials generated by gener-
ate bounded polynomial.
generate gradient flow: When set True, the backward generation only generates gradient flows systems.
gen gen_weight: exponential weight which bias the sampling of and p. the number of components of non-zero hę and fi-
max_order Border_purr____ pale::ppllylooel maximal polynomial order of h;
max ____ee_r_fwwdd maximal number of terms in each equations in the fwd generation
SOS SoS_checker: True, uses a SOS verifier to evaluate the candidate Lyapunov function (if
False uses the verifier based shgo)
SMT checker: True, uses an SMT verifier to evaluate the candidate Lyapunov function (if
False uses the verifier based shgo)
multigen: number different system generated per Lyapunov function.
increasing fune: the of increasing functions used in the generation (see Step 1b). Default is (exp, In(1 + x²),, ) /1 +x)- x).
positive_ fositive_func: the set of positive functions used in the generation (see Step 1b). Default is
fexp, 1 + (x) + sin(z))-
bounded_ fune: the of bounded functions used in the generation (see Step lb). Default is {cos, sin}.
B.5 Forward SOS solver
SOSTOOLS is one of the most famous toolbox for sum-of-square optimization, in particular for finding SOS Lyapunov functions [32 33 It is natively available in MATLAB and relies an underlying SDP solver that can be chosen. In Python an analogous toolbox is the package SumOf. Squares 4 which relies on the same principle, however does not have specific functionalities for Lyapunov functions. a consequence we implemented these functionalities our codebase based on the MATLABB implementations in SOSTOOLS We implemented a function c0S_checkerr which takes in input system of equations in sympy and a candidate Lyapunov function and checks SOS conditions on (zx and - VV(x)) f(),, and function f findlyap. analogous the findlyap
17

---

function in SOSTOOLS, which takes system of equations in sympy and either returns a function satisfying SOS conditions on (a) and (v f(r) returns false no such function exists, or returns none if it fails to provide an answer. SumOffquares relies itself on picos 33 and we use the
default solver evxopt [
B.6 List of datasets
Dataset
Description
Size (000)
Resources (CPU.hourr)
BPoly
Backward polynomial systems,
1,000
210
BNonPoly
Backward non-pon pol-poooyoilll systems, non- non-zero
1,000
220
FBam
Forward, non-homogeneous polynomial barrier functions
300
9,670
FLyap FSOSTOOLS
Forward, homogeneous polynomial Lyapunov functions
4,620
Forward, SOSTOOLS solved systems
1.5
Table 8: Datasets generated. Backward systems are degree 2 to 5, forward systems degree to 3. All forward
systems polynomial.
C Additional results
C.1 Impact of multigeneration
In the backward generation procedure, after sampling one random V, it is possible to generate any number different systems fs such that V the Lyapunov function for each of the systems ft. We call the maximal number of system generated per Lyapunov function the multigen parameter. The actual number of systems generated per Lyapunov function is chosen random for each Lyapunov function between and multig mugtigen. In Section [eee reported results using mult mignigen equal to 5. Here we report the in-domain and out-of-domain performance of the models trained on backward BPoly datasets of size million varying the parameter mul tigen.
In-domain BPoly
OOD
Multigen
FLyap
95
58
5
100
73
10
100
25
100
76
50
100
70
100
100
68
Table 9: In-domain and out-of-domain accuracy of models. Beam size
Table | shows that generating a moderate amount different systems with the same Lyapunov function actually improves the model capability generalize out-of-domain. This suggests that the model learning, at least partially, to separate the parts of the system which contribute to the Lyapunov function. Above certain multigen threshold, model performances start decline. This may be due to the low diversity present in the dataset, i.c. the limited number of different Lyapunov functions the model trained on (the total number of systems the training remains constant so the total number of Lyapunov function decreases with the value of the parameter multigen).
C.2 Performance of smaller transformer models
In Section we report results using transformer with encoder and decoder layers, attention heads and an embedding dimension of 640. We also trained smaller models with encoder and decoder layers, attention heads and an embedding dimension 512. Tables TO report the main results. Results are in line with what we showed in section
18

---

In-domain
OOD
In-domain
OOD
Backward datasets
FLyap
8.4/9.9 Forward datasets
BPoly
BPoly (polynomial)
100
70
FBarr (barrier)
97
13
BNonPoly (non-poly)
71
FLyap (Lyapunov)
86
Table 10: In-domain and out-of-domain accuracy of models. Beam size 50.
Forward
Mixing proportion
Into the wild
datasets
FBarr
FLyap
Poly3
Poly5
No mixture
0%
31
70
2.3
1.6
FBarr
0.01 %
60
72
2.5
1.7
0.1% 11%
93
72
9.2
6.4
FLyap
0.01 %
29
73
2.8
1.6
0.1% 1
19
76
2.9
1.7
Table 11: Performance mixing backward data (BPoly) with small number of forward examples on forward benchmark and "into the wild". Beam size 50.
D Comparison of SOS, SMT and shgo
We compare model performance when we employ them discover Lyapunov function. We report performances with dReal SMT and SOS verifiers for Poly and dReal SMT and shgo for NonPoly distributions, respectively. Table 12 shows that SMT results are slightly lower, because of timeouts (which we report in Table [3] but comparable. Note that the performances on polynomial systems were already theoretically guaranteed thanks to the former SOS verifier.
Model (by training distribution)
FBarr
Bixtur
NonPolyMixture
Poly3
10.5/11.7 6.5/8.0
11.1/11.8 8.7/10.1
10.6/11.2
Poly5
NonPoly
8.3/12.7
Table 12: Results of SMT with SOS and shgo verifiers for Poly and NonPoly systems, respectively.
Test
Into the wild
Timcout
10 minutes
minutes
Correct Lyap function
87.3
92.2
SMT Timeouts
10.8
5.8
Incorrect Lyap function
1.9
2.0
Table 13: SMT timeout and rates. Most SMT failures due timeout.
E AI method sweep
To report the Al-based tools results on the seven benchmarks (BPoly, BNonPoly, FLyap, FBar,, Poly3, Poly5, NonPoly) we did a hyperparameter sweep. To get the best hyperparameter setting, sweep on FLyap and then fix these hyperparameters for the different datasets. In bold we show the chosen parameters, selected to maximize the correctness on FLyap, subject to the 20 minutes timeout.
Lyznet 25
Ir [3 10 5 10 4, 4. 10=4]
points 1100,000. 300,000, 1,000,000]
layer width = [(2,20), (3,6), (6,2)]
epoch [1, 25]
net type = [None, Poly]
19

---

Fossil 2.0 [12
iters [10, 250]
activations = [(x²), (,?_,²), (sigmoid), (sigmoid, sigmoid), () (poly,, poly,))
hidden neurons [6, 10, 20)
data = [500, 1000, 2000]
Ir [0.01, 0.03, 0.1]
ANLC v2 [15]
iters [10, 50, 250]
activations [(x x. x), (z², x² 2). (x (x2.=², x, 2), (a3 21))
hidden neurons [6. 10, 20]
max data [500, 1000, 2000]
Ir [0.01, 0.03, 0.1]
F
Some examples
To understand the model performance and compare against the SOSTOOL performance, perfonnnc,,,,,ae,, we manually inspect some systems with 2 or 3 equations where the following conditions hold: (1) the Jacobian of the system has the maximum eigenvalue with real part equal to 0 (i.e. tools like the spectral mapping theorem cannot decide on the stability), (2) weighted 2-norm functions can be a Lyapunov function, (3) f findlyap times out after hours. We show some examples below.
F.1
A polynomial system with polynomial solution
System
Tig + TD0
-1
It's known that there is no polynomial Lyapunov function for this system [2]. Our poly models and findlyap failed, as expected. Nonetheless, one of our non-poly models with beam search of beam size 100 proposed (xx) = Im(1 + 5xz)) + xf similar to the one that was recently found in [2].
It's clear that ((0) 0 and V(x) > 0 for all x * 0. Also
V(x) f(x) = 100x 10x 2z (1 + 5x/)
+ 5x
(17)
-52 - 5
5(20 - 20 )² 2x
< 0
1 1550
as desired.
F22 A system that has no diagonal Lyapunov function
System
i 2x
[ = -10x1
Model inference: Our model recovers V V(x) = 10x3 + 2x002 + 3xt + 6z7.
Clearly V(0) 0 and (xx = ((.)) + (No + x²)² + 2(x7)² + 6z² 0 for all x # 0. Also
VV(x) f(x) (1162 + 120) <
Non existence of Diagonal Lyapunov function: Suppose for the sake of contradiction that there exists function Vi which satisfies and can be expressed as
Vi = Σ
a
+
b b
20

---

Clearly (1( 0. Given that (ro, 0) 0 for 20 * 0, it follows that n is even and on > 0. Also "
m
we know that (v) f(x) = Σ
00
10
jbyay S 0 for all choices of (1). If Σ
j=1
m
we 11 we obtain (vi f(x)
Σ
10
Σ
jby. This expression can be seen
i=1
j=1
as polynomial g(0o) with real coefficients and odd degree - 1. The leading coefficient, 2na,, is
positive because on and This means that lim g(zo) +0, meaning that there exists 20 #00+00
an Ty such that gxx)) > 0. This contradicts
F.3 A system with equations and a higher degree
System
ž - 7 4xjaa - 5
d1 7x
3z1
20²
i2
82
9x2
Model inference: Our model recovers different solutions. Here we show two of them
(xx 4x4 + 10m + 20 + 10m 4m + 20x + 10x + 4x 2x1 - 88 + 4m2.
(2(x) 2x/ + +2rzz 32 2zf +22
We checked with SumOfSquares that V > 0, > 0, VV f < 0 and VV f S O.
F.4
Other examples
System
Lyapunov function
40 500 2xorr
(xx) = 6zØ 7 xu 10ng + 8n
41 9r1
30
1x1
4i =1 1x
9x0at + 3xga
V(x) = zt + 9x. + 3x7
Ø1 = Bagar
10x0x1 + 3xp2s 7
4u 3xd
3roD2 920
fi d 5xz + 5x
V(x) =t + Zugu3 + 353 + 4xpw 3af 2ng + 10z
422 -x
4i 8 2 10x1
F1
8r + 3x9 821
V(x) = 4x 2ngor + 6zz + 4xf + x3
12
Table 14: Some additional generate from our models.
21

---

NeurIPS Paper Checklist
1. Claims
Question: Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?
Answer: [Yes]
Justification: The claims in the abstract and introduction summarize the results of section "Main results".
Guidelines:
The answer NA means that the abstract and introduction do not include the claims made the paper.
The abstract and/or introduction should clearly state the claims made, including the contributions made in the paper and important assumptions and limitations. A No or
NA answer to this question will not be perceived well by the reviewers.
The claims made should match theoretical and experimental results, and reflect how much the results can be expected to generalize to other settings.
is fine to include aspirational goals as motivation as long as is clear that these goals not attained by the paper.
2. Limitations
Question: Does the paper discuss the limitations the work performed by the authors?
Answer: [Yes]
Justification: Yes, in the discussion.
Guidelines:
The answer NA means that the paper has no limitation while the answer No means that the paper has limitations, but those are not discussed in the paper.
The authors are encouraged to create a separate "Limitations" section their paper.
The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (e.g.. independence assumptions, noiseless settings, model well-specification, asymptotic approximations only holding locally). The authors should reflect how these assumptions might be violated in practice and what the implications would be.
The authors should reflect on the scope of the claims made, e.g., if the approach was only tested on few datasets or with few runs. In general, empirical results often depend on implicit assumptions, which should be articulated.
The authors should reflect on the factors that influence the performance of the approach. For example, facial recognition algorithm may perform poorly when image resolution is low images are taken in low lighting. Or speech-to-text system might not be used reliably to provide closed captions for online lectures because it fails to handle technical jargon.
The authors should discuss the computational efficiency of the proposed algorithms and how they scale with dataset size.
If applicable, the authors should discuss possible limitations of their approach to address problems of privacy and fairness.
While the authors might fear that complete honesty about limitations might used by reviewers as grounds for rejection, worse outcome might be that reviewers discover limitations that aren't 't acknowledged judgment and recognize that individual actions in favor of transparency play impor- tant role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations.
the paper. The authors should use their best
3. Theory Assumptions and Proofs
Question: For each theoretical result, does the paper provide the full set of assumptions and complete (and correct) proof?
Answer: [Yes]
22

---

Justification:
This paper is mostly experimental, the only theoretical results is the stability specific systems and we do verify the assumptions and state the theorem we use about Lyapunov functions.
Guidelines:
The answer NA means that the paper does not include theoretical results. All the theorems, formulas, and proofs in the paper should be numbered and cross- referenced.
All assumptions should be clearly stated or referenced the statement of any theorems.
The proofs can either appear in the main paper or the supplemental material, but they appear in the supplemental material, the authors are encouraged provide short proof sketch to provide intuition.
Inversely, any informal proof provided in the core of the paper should be complemented by formal proofs provided in appendix or supplemental material.
Theorems and Lemmas that the proof relies upon should properly referenced.
4. Experimental Result Reproducibility
Does the paper fully disclose all the information needed to reproduce the main ex- perimental results of the paper to the extent that affects the main claims and/or conclusions of the paper (regardless of whether the code and data are provided or not)?
Question:
Answer: [Yes]
Justification: Yes, we detail the experiments, the parameters and the data generator. The code will eventually be provided.
Guidelines:
The answer NA means that the paper does not include experiments. If the paper includes experiments, a No answer this question will not perceived well by the reviewers: Making the paper reproducible is important, regardless of whether the code and data are provided or not.
the contribution is a dataset and/or model, the authors should describe the steps taken to make their results reproducible or verifiable.
Depending on the contribution, reproducibility can be accomplished in various ways. For example, if the contribution is a novel architecture, describing the architecture fully might suffice, or if the contribution a specific model and empirical evaluation, may be necessary to either make it possible for others to replicate the model with the same dataset, or provide access to the model. In general. releasing code and data often one good way to accomplish this, but reproducibility can also be provided via detailed instructions for how to replicate the results, access to a hosted model (e.g., in the case of large language model), releasing of a model checkpoint, or other means that are appropriate to the research performed.
While NeurIPS does not require releasing code, the conference does require submis- sions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. For example
(a) If the contribution is primarily a new algorithm, the paper should make it clear how reproduce that algorithm.
(b) If the contribution is primarily new model architecture, the paper should describe the architecture clearly and fully.
(c) the contribution is a new model (e.g., a large language model), then there should either be way access this model for reproducing the results way to reproduce the model (e.g., with an open-source dataset or instructions for how to construct the dataset).
(d) We recognize that reprodueibility may be tricky in some cases, in which authors are welcome to describe the particular way they provide for reproducibility. In the case of closed-source models, it may be that access to the model is limited in some way (e.g., to registered users), but it should be possible for other researchers have some path reproducing verifying the results.
5. Open access to data and code
23

---

Question: Does the paper provide open access to the data and code, with sufficient instruc- tions faithfully reproduce the main experimental results, as described in supplemental material?
Answer: [No]
Justification: The code and datasets currently pending approval from some of the authors' institution for being open-sourced. They should be by the conference dates.
Guidelines:
The answer NA means that paper does not include experiments requiring code.
Please see the NeurlPS code and data submission guidelines (https://niis.cc.
cublicc/qudees/Codeuubbss catpssionoolyyyy for more details.
While encourage the release of code and data, understand that this might not be possible, so "No" is an acceptable answer. Papers cannot be rejected simply for not including code, unless this is central to the contribution (e.g., for new open-source benchmark).
The instructions should contain the exact command and environment needed to run to reproduce the results. See the NeurlPS code and data submission guidelines (https wtt.ee/eeeoSSSeeese..eeeee..
/ //nips. wwtps/cccppplgguuuiaaeeeae
www.sssoiicgcyo for more details.
The authors should provide instructions on data access and preparation, including how to access the raw data, preprocessed data, intermediate data, and generated data, etc.
The authors should provide scripts to reproduce all experimental results for the proposed method and baselines. only subset of experiments are reproducible, they should state which ones are omitted from the script and why.
At submission time, to preserve anonymity, the authors should release anonymized versions (if applicable).
Providing as much information as possible in supplemental material (appended the paper) is recommended, but including URLs to data and code is permitted.
6. Experimental Setting/Details
Question: Does the paper specify all the training and test details (e.g., data splits, hyper- parameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?
Answer: [Yes]
Justification:
We specify the main training and test details.
Guidelines:
The answer NA means that the paper does not include experiments.
The experimental setting should be presented in the core of the paper to a level of detail that is necessary to appreciate the results and make sense of them.
The full details can be provided either with the code, in appendix, as supplemental material.
7. Experiment Statistical Significance
Question: Does the report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?
Answer: [No]
Justification: We provide the number of elements in each sets, well the accuracies in percentage. Having error bars is not adapted when dealing with the infinite dimensional spaces of functions involved.
Guidelines:
The answer NA means that the paper does not include experiments.
The authors should answer "Yes" if the results are accompanied by emor bars, confi- dence intervals, or statistical significance tests, least for the experiments that support the main claims of the paper.
24

---

The factors of variability that the error bars are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall run with given experimental conditions).
The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.)
The assumptions made should be given (e.g., Normally distributed errors).
It should be clear whether the error bar is the standard deviation the standard error of the mean.
is OK to report 1-sigma error bars, but one should state it. The authors should preferably report a 2-sigma error bar than state that they have 96% CI, if the hypothesis of Normality of errors is not verified.
For asymmetric distributions, the authors should be careful not to show in tables figures symmetric error bars that would yield results that are out range (c.g. negative error rates).
If error bars reported in tables or plots, The authors should explain in the text how they were calculated and reference the corresponding figures or tables in the text.
8. Experiments Compute Resources
For each experiment, does the paper provide sufficient information on the com- puter resources (type of compute workers, memory, time of execution) needed to reproduce the experiments?
Question:
Answer: [Yes]
Justification: We provide the type of compute workers, the type GPU, the memory, and the GPU and CPU time needed for training the models and generating the datasets.
Guidelines:
The answer NA means that the paper does not include experiments.
The paper should indicate the type of compute workers CPU or GPU, internal cluster, or cloud provider, including relevant memory and storage.
The paper should provide the amount compute required for each of the individual experimental rums well estimate the total compute.
The paper should disclose whether the full research project required more compute than the experiments reported in the paper (e.g., preliminary or failed experiments that didn't make into the paper).
9. Code Of Ethics
Question: Does the research conducted in the paper conform, in every respect, with the NeurlPS Code of Ethics
Answer: [Yes]
Justification: This paper is about training neural networks find solutions to mathematical problems. As such, it's societal aspects (safety, security, discrimination, surveillance, deception and harassment, environment, human rights, Bias and fairness) are either irrelevant respect the NeurlPS Code of Ethics. The only data coming from human sources are answers to mathematics exercises from 25 students as part one of their course. All students were all offered the choice of allowing or not the use of these data for the purpose research, the answers of students who agreed were then anonymized and aggregated.
Guidelines:
The answer NA means that the authors have reviewed the NeurIPS Code of Ethics.
the authors answer No, they should explain the special circumstances that require deviation from the Code of Ethics.
The authors should make to preserve anonymity (e.g., if there is special consid- eration due to laws or regulations their jurisdiction).
10. Broader Impacts
Question: Does the paper discuss both potential positive societal impacts and negative socictal impacts of the work performed?
25

---

Answer:
[NA]
Justification: While there no direct societal impact, there is an impact on the research in mathematics and science in general and we mention it.
Guidelines:
The answer NA means that there is no societal impact of the work performed.
If the authors answer NA or No, they should explain why their work has no societal impact or why the paper does not address societal impact.
Examples of negative societal impacts include potential malicious or unintended uses (e.g.. disinformation, generating fake profiles, surveillance), faimess considerations (e.g., deptoyment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations.
The conference expects that many papers will be foundational research and tied to particular applications, let alone deployments. However, if there is direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster.
The authors should consider possible harms that could arise when the technology is being used intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional unintentional) misuse of the technology.
If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms to monitor how a system learns from feedback over time, improving the efficiency and accessibility of ML).
11. Safeguards
Does the paper describe safeguards that have been put in place for responsible release of data models that have high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?
Question:
Answer: [NA]
Justification: The models small language models trained solve specific mathematical question.
Guidelines:
The answer NA means that the paper poses no such risks.
Released models that have high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere to usage guidelines or restrictions to access model or implementing safety filters.
Datasets that have been scraped from the Internet could pose safety risks. The authors
should describe how they avoided releasing unsafe images.
We recognize that providing effective safeguards is challenging, and many papers do not require this, but we encourage authors to take this into account and make best faith effort.
12. Licenses for existing assets
Are the creators or original owners of assets (c.g., code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected?
Question:
Answer: [Yes]
Justification: With the exception of common Python's packages for which only give the name, we cite the specific packages we are relying on, when relevant. No data or model from another source is used.
Guidelines:
26

---

The answer NA means that the paper does not use existing assets.
The authors should cite the original paper that produced the code package dataset.
The authors should state which version of the asset used and, if possible, include URL.
The name of the license (e.g., CC-BY 4.0) should be included for each asset.
For scraped data from a particular source (e.g.. website). the copyright and terms of service of that source should be provided.
If assets are released, the license, copyright information, and terms of use in the package should be provided. For popular datasets, paaassttssppersssaaes paperswi thcode. con/datasess has curated licenses for some datasets. Their licensing guide can help determine the license of dataset.
For existing datasets that are re-packaged, both the original license and the license of the derived asset (if it has changed) should be provided.
If this information is not available online, the authors are encouraged to reach out to the asset's 's creators.
13. New Assets
Question: Are new assets introduced in the paper well documented and the documentation provided alongside the assets?
Answer: [NA]
Justification:
No new assets are released yet.
Guidelines:
The answer NA means that the paper does release new assets.
Researchers should communicate the details the datasct/code/model as part of their submissions via structured templates. This includes details about training, license, limitations, etc.
The paper should discuss whether and how consent was obtained from people whose asset used.
At submission time, remember to anonymize your assets (if applicable). You can either create anonymized URL or include an anonymized zip file.
14. Crowdsourcing and Research with Human Subjects
Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given participants and screenshots, if applicable, as well as details about compensation (if any)?
Answer: [No]
Justification: The data coming from human sources are only an element of the discussion and not the main point of the paper. The framework and instructions are given as well in the discussion.
Guidelines:
The answer NA means that the paper does involve crowdsourcing nor research with human subjects.
Including this information in the supplemental material is fine, but if the main contribu- tion of the paper involves human subjects, then as much detail as possible should be included in the main paper.
According to the NeurlPS Code of Ethies, workers involved in data collection, curation, or other labor should be paid at least the minimum wage in the country of the data collector.
15. Institutional Review Board (IRB) Approvals or Equivalent for Research with Human Subjects
Question: Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or an equivalent approval/review based on the requirements of your country or institution) were obtained?
27

---

Answer: [Yes]
Justification: The only data coming from human sources are answers to mathematics exercises from 25 students part of one of their course. All students were all offered the choice of allowing or not the use of these data for the purpose of research, the answers of student who agreed were then anonymized and aggregated. This done in coordination with the administration of the mathematics department.
Guidelines:
The answer NA means that the paper does not involve crowdsoureing nor research with human subjects.
Depending the country in which research is conducted, IRB approval (or equivalent)
may be required for any human subjects research. If you obtained IRB approval, you should clearly state this in the paper.
We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors to adhere to the NeurIPS Code of Ethics and the guidelines for their institution.
For initial submissions, do not include any information that would break anonymity (if applicable), such the institution conducting the review.
28


Reviews:
[{"review_id": "3LUk6iiuDd", "reviewer": "Reviewer_uYdz", "summary": "This study addresses the problem of designing the Lyapunov function of dynamical systems via learning. The existence of the Lyapunov function is a sufficient condition of the dynamical system's stability. However, its design is not established except for systems with sum-of-square polynomials. The Authors propose a novel dataset construction framework on which Transformer models can be trained. Extensive experiments demonstrate that Transformer models successfully learn to design Lyapunov functions on several types of datasets with moderate generalization ability across datasets. It is also shown that the injection of out-distribution samples in the training set significantly boosts the success rate. Yet restricted to small systems, this study suggests that the learning approach is a promising direction for addressing hard problems of symbolic mathematics.", "questions": "I'd like the Authors to answer the weaknesses raised above. Plus, \n- What do the \"pre-defined set of increasing functions\" and \"pre-defined set of bounded-functions\" refer to at [lines 479, 485]?", "limitations": "Limitations are adequately presented. For the potential improvements, see Weakness and Questions.", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 2, "contribution": 3, "strengths": "- This study clearly presents the problem and the approach. Overall, the writing is clean and accessible to diverse readers. \n- A novel backward dataset generation algorithm is proposed to train Transformer models for Lyapunov function design. \n- The experiments validate the proposed approach from various aspects. Particularly, the generalization across datasets is tested well.", "weaknesses": "I appreciate the overall contributions of this study. I'd like to raise two weaknesses, one for the methodology (backward generation) and one for the experiments. \n\n**Major comments**\n\nBackward generation. While the Author claims that the backward generation is the key innovation, the explanation of the generation process is very limited in the main text. I understand the restriction of pages, but it would be better to elaborate on the method a little more to give intuition to the readers. While the details are presented in Appendix A.2, it mainly provides the procedures, and the rationale is not much given. I was not able to fully understand the procedures either because of the unclarity of several notations (see Minor comments below). Plus, many mapping and transformations appear, but their intention and necessity are unclear to me. I encourage the authors to discuss the range of the Lyapunov function and $f$ covered by the proposed generation method and whether this is large enough or not, including, if any, the important class of functions that cannot be covered for now.\n\nAs for the experiments, assuming a practical use, it is important to know whether the Lyapunov function given by Transformers is really the Lyapunov function of the targeted system, and if not, whether this is because of a simple failure or the existence of such function. Discussion and experiments seem missing in this part. \n\n**Minor comments**\n\n*(Markdown math rendering somehow fails in the first item below, so I wrote it in a bit tricky way.)*\n- Some notations are unclear and inconsistent. For example, vectors $\\mathbf{z}$ and $\\mathbf{z}_{j}^{i}$ are in bold font, while other vectors such as $x$ are not. The subscript notations $A$ [line 177] \nand $B$ in Eq. (10) are not defined. Assuming that it refers to the $i$-th entry of vector and $C$ is $D$, the first term of the equation in [line 177] appears to be a scalar while the second term is a vector. The Authors may represent a vector by $E$ in the first term, but still, this is a row vector while the second term is a column vector. \n\n$A = (\\nabla V)_i(x)$\n\n$B = (\\nabla V(x))_{\\tau_2(i)}$\n\n$C = (\\nabla V)_i(x)$\n\n$D = (\\nabla V(x))_i$\n\n$E = (h_{\\pi(i)})_i$\n\n- Forward/backward methods and Forward/backward generations are both used. I recommend that the Authors keep terminology consistent.\n- [line 183] and the is --> and there is \n- [line 224] .In -> . In\n- [line 489] What does $Id$ mean (the identity map?)?\n- [line 492] \"$g$ is positive\" souds a bit wierd as $g$ is a function.\n- [line 489] Written $g= 1$ but $g$ is supposed to be a function.\n- [Eq. 10] Does $\\mathbf{e}^i_j$ refers to the $j$-th entry of $\\mathbf{e}^i$? Then, it should be unbold as it is a scalar.\n\nI also encourage the Authors to include more related works that exploit Transformers to address hard math problems. \n\n**Shortest Vector Problem (a series of works by the same group)**\n- \"SALSA: Attacking Lattice Cryptography with Transformers,\" Emily Wenger, Mingjie Chen, François Charton, Kristin Lauter\n\n**Gröbner basis computation**\n- \"Learning to Compute Gröbner Bases,\" Hiroshi Kera, Yuki Ishihara, Yuta Kambe, Tristan Vaccon, Kazuhiro Yokoyama\n\nIf not restricted to Transformers, \n\n**Detection of terminal singularity**\n- \"Machine learning detects terminal singularities,\" Tom Coates, Alexander M. Kasprzyk, Sara Veneziale\n\n**Integer programming**\n- \"Machine Learning for Cutting Planes in Integer Programming: A Survey,\" Arnaud Deza, Elias B. Khalil", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your review and comments.\n\n**Backward generation. it would be better to elaborate on the method a little more.**\n\nWe agree. See section B of the author rebuttal for a better presentation of our methods, explaining the motivation of the different steps. We will update the paper in this direction, and add further explanations in the appendix.\n\n**I encourage the authors to discuss the range of the Lyapunov function and 𝑓**\n\nThank you for bringing this up. One of the main advantages of the generation method we propose is that it allows us to parametrize the functions $V$ and $f$, by selecting the functions $V_\\text{proper}$, $V_\\text{cross}$, $g_i$ and $h_i$ in different classes. This enables us to generate polynomial or non polynomial Lyapunov functions, for polynomial or non polynomial systems. At present, we generate functions defined by the four operations, power, and usual functions (exp, ln, cos sin…), but this could be extended to any other functions, so long we can specify their gradient. As a consequence, our framework can generate any function $f$ that has a closed formula.\n\nThe main limitation is the restriction on $V$. According to definition 3.2, any non negative function with a strict minimum of 0 in the origin, and infinite at infinity would do, but sampling such a function in the general case is likely a NP-hard problem. Instead, we sample $V$ in the smaller class of functions that can be written $V = V_\\text{proper}+V_\\text{cross}$, with $V_\\text{cross}$ a product of non-negative functions, such as sum-of-squares, not necessarily polynomial, composed with bounded functions, and $V_\\text{proper}$ either a positive definite homogeneous polynomial (in the polynomial case) or a function of the form $V_\\text{proper} = \\sum_{i=1}^{n} f_{i}(x_{i})$ with $f_i$  a composition of an even polynomial with a strictly increasing function (currently limited to $\\exp$, $\\ln(1+x^2)$, $\\sqrt{1+x}$ but this could be extended). \n\nStill, this is a larger class than the functions considered in SOS approaches where $V_\\text{proper}$ must be the second form with $f_i$ even polynomials and $V_\\text{cross}$ a sum of squares. It also covers all the examples we have seen in textbook examples.\n\nAnother limitation is that, for now, we reject functions $f$ that are defined on a smaller domain than $\\mathbb{R}^{n}$ because we are interested in the **global** stability problem. Future improvement would include operators with a restricted domain of definition and focus on the stability problem on a restriction of $\\mathbb{R}^{n}$.\n\nFollowing your suggestion, we’ll add a discussion on the class of functions that can be covered and the limitations in our revised version.\n\n**it is important to know whether the Lyapunov function given by Transformers is really the Lyapunov function of the targeted system, and if not, whether this is because of a simple failure or the existence of such function.**\n\nThank you for the comment. All model predictions are passed to a verifier that checks whether it is a correct Lyapunov function for the input system. For polynomial systems, the verifier is a theoretical checker (based on sum-of-squares), for non-polynomial systems it was a numerical checker. In the new version we added a theoretical checker for the non-polynomial case (based on SMT, see our author's rebuttal, section D).\n\nThe use of an external checker guarantees that all verified model predictions are indeed correct Lyapunov functions. It may happen, however, that correct model predictions are verified as incorrect, because of a failure in the checker. \n\nThere are four cases:\n* The verifier is too restrictive. Polynomial systems may have non-polynomial Lyapunov functions, but SOS checkers always assume sum-of-square Lyapunov functions. \n* The verifier aborts because of timeout or memory overflow, while checking a correct solution.\n* The verifier aborts because of timeout or memory overflow, while checking a wrong solution, for a system that is globally stable\n* The verifier aborts because the system is not globally stable (so no Lyapunov function exists)\n\nThe first case can be identified by changing the verifier, we believe it is rare (polynomial systems with non-polynomial functions are a very special case).\n\nDiscriminating between the second and two last cases amounts to telling whether the verifier aborted for “good or bad” reasons. This cannot be done in a systematic way, but we can estimate the frequency of these failures by counting the number of verifier failures on the generated datasets, for which we know the Lyapunov functions. Following your suggestion, we estimated these failure rates and will report them in the updated paper. For non-polynomial systems with known Lyapunov functions, we observe that the SMT solver fails in about 18% of the cases.\n\nDiscriminating between the two last cases amounts to determining whether a system is globally stable. This is an open question, in fact the only way to verify this is to find a Lyapunov function. The only estimate we can have is for special classes of functions that we know are globally stable (polynomials SOS methods can solve, gradient flow systems deriving from a potential). We will report the accuracy on gradient flow systems in the revised version.\n\n**What do the \"pre-defined set of increasing functions\" and \"pre-defined set of bounded-functions\" refer to at [lines 479, 485]?**\n\nCurrently the pre-defined set of increasing functions is ($\\exp$, $\\sqrt{1+x)}$, $\\ln(1+x)$) and the pre-defined set of bounded-functions is (cos, sin). We didn’t initially specify them in the procedure because these sets are generation parameters that can be user-defined (one only has to specify the function and its gradient), but we have included them now in the revised version, thanks for the comment.\n\n**Minor comments** \n\nThank you for spotting these. We corrected them, and added the references you suggested.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "mGglADdZmQ", "reviewer": "Reviewer_LBjn", "summary": "This paper trains a transformer to predict Lapunov functions for both polynomial and nonpolynompial dynamical systems. They generate a dataset combining \"forward-generated\" problems — solutions discovered from existing solvers on randomly generated problems — with \"backward-generated\" problems — problems generated by generating a random solution, and then building a function solved by the solution. The trained transformer achieves great performance. On a test dataset of randomly generated problems, it solves 12% of problems, compared to 1% solved by the state of the art solver. In addition, the model solves 84% of 75 problems in the FSOTOOLS dataset where human experts only solved 9.33%.", "questions": "1. findylap solves the same number of problems as SOSTOOLs, right? \n2. aren't the other forward datasets also datasets that SOSTOOLS can solve, since its used to generate them? if so, what decision went into the characteristics of the test set over those used for training? \n3. are the polynomial systems for FSOSTOOLS homogenous, nonhomogenous, etc? I'm not an expert on the math at all, but those adjectives were used for the train set. I'm trying to understand, because the performance on the test set is the main way of knowing how good a ML model is. but if the test set is too similar to the train set, then the results aren't as meaningful. of course, you have the results on the \"random systems\" dataset.\n4. One idea could be to try to do some sort of expert iteration with your model: use it as the base solver when generating a forward dataset, and use newly solved problems (since it's better at solving than `findylap` on the random systems dataset) to further fine-tune the model. this probably isn't that useful for increasing performance, but it would be a really cool proof of concept to see the performance keep going up on the random systems dataset with some iterations!\n5. Is there a real-world (or \"real math\") application of finding Lyapunov functions? for example, maybe there are some applications. What kind of things do the existing tools get used for, and how might this model be used as a superior version?", "limitations": "The limitations are addressed well in the discussion section.", "rating": 8, "confidence": 3, "soundness": 4, "presentation": 4, "contribution": 4, "strengths": "- The paper is well-written, and clearly described, so that I could understand it even without having previously encountered Lyapunov functions. \n- The related work section is good\n- The dataset generation is very carefully described, and lots of important design decisions are included. There is even more detail in the Appendix. Despite not being an expert, it seems like the authors really know what they're doing, and have rigorously tested different interesting choices of what to include in the dataset. \n- The results are really cool! The authors test out of domain generalization well, by showing performance on a randomly generated dataset. Indeed, some of the models evaluated were only trained using \"ground truth\" solver solutions from a different randomly generated dataset. But  the model gets 10x higher performance than the ground truth solver on this random dataset.  \n- The big question for papers like this is how well the models generalize out side the train set. Due to the difficult nature of the problem, it's hard to have really good test sets, but the authors are aware of the importance of this, and do lots of OOD generalization tests of their systems, and report positive results.\n- Again, the authors are quite thorough with their documentation and evaluation of their system. The paper feels very polished and professional and the results are impressive.\n- Progress on this problem could lead to other work using deep learning for mathematics.", "weaknesses": "1. The algorithmic approach is not new — training a transformer to do math problems that a large dataset can be generated for — but this is okay given the careful effort needed and given to the dataset creation process, and the importance and performance of the problem achieved. \n\n2. I find all the different dataset variants included to be a bit confusing. The paper could be easier to read with a reduced the number of variants presented in the main section of the paper. \n\n3. I think the discussion section goes a little far in making suggestions about LLM reasoning ability given success on this dataset. The sentences are true, but it's not necessarily true that performance here is equivalent to reasoning. the transformers are probably doing sometihng more like \"intuition\" given the amount of data they are seeing and the quickness at which they generate an answer.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your review and comments. Here are a few answers and clarifications.\n\n**The algorithmic approach is not new**\n\nWe agree. See section A of the authors’ response for a clarification of our novel contributions. Summarizing,\n* we demonstrate the applicability of backward generation to an open problem (instead of one that we know how to solve, like integration). [Yehudah et al, ICML, 2020] explained why this is not straightforward. The adjustments we had to implement (section B off the author response) prove their point,\n* we show that mixing a tiny number of forward examples into the backward training set helps mitigate the OOD problem.\n  \n**I find all the different dataset variants included to be a bit confusing**\n\nWe will change this in the revised version, and limit the main paper to four datasets: \n* Two forward sets: Fbarrier and FLyap, for two different problems: barrier functions and Lyapunov functions,\n* Two backward sets: BPoly and BNonPoly, for polynomial and non-polynomial systems.\n\nResults on other datasets will be moved to the Appendix.\n\n**the discussion section goes a little far in making suggestions about LLM reasoning ability given success on this dataset. (...) The transformers are probably doing something more like \"intuition\"**\n\nWe agree, and we will clarify this paragraph. We used the word “reasoning” because it is common in current LLM literature, but we agree that its overuse makes the claim confusing. We changed it to intuition as you suggest.\n\nIn particular, we propose the following rewrite of the end of the discussion: \n\nOur results suggest that transformers can indeed be trained to discover solutions to a hard problem of symbolic mathematics that humans solve through reasoning, and that this is enabled by a careful selection of training examples, instead of a change of architecture. We do not claim that the Transformer is reasoning but it may instead solve the problem by a kind of “super-intuition” that stems from a deep understanding of a mathematical problem.\n\n**Questions**\n\n**1- findylap solves the same number of problems as SOSTOOLs, right?**\n\nYes. findlyap is a python analog of SOSTOOLS, that we needed because our computing cluster cannot run Matlab. It solves the same problems as SOSTOOLS, up to differences in speed and memory usage, which may cause different timeouts and failures.\n\n**2- aren't the other forward datasets also datasets that SOSTOOLS can solve (...)? if so, what decision went into the characteristics of the test set over those used for training?**\n\nYes. Since SOS methods are the main techniques known to find the symbolic global Lyapunov functions we need for our forward training sets, in this paper, the word “forward” is essentially equivalent to “solvable with  SOSTOOLS” (for this reason, findlyap accuracies on forward datasets, 100% by design, are not reported in table 4).\n\nThe three forward datasets include systems of increasing difficulty for SOSTOOLS. FHom and FNHom both focus on the easier “barrier Lyapunov function discovery” problem. FNHomProper, focuses on the harder problem of discovering a Lyapunov function for a general polynomial system. Number of equations, degree and range of coefficients were selected so that findlyap could find a Lyapunov function in a reasonable amount of time (a few minutes per problem).\n\nThe forward test sets, being held-out samples of the forward training sets, have the same distribution. We ensure that no overlap exists between train and test sets. \n\n**3- are the polynomial systems for FSOSTOOLS homogenous, nonhomogenous, etc?**\n\nThis was indeed missing, thank you for pointing it out. FSOSTOOLS contains non-homogeneous proper examples that correspond to the Lyapunov function problem with the most generic distribution of polynomials. \n\n**4- try to do some sort of expert iteration with your model: use it as the base solver when generating a forward dataset, and use newly solved problems to further fine-tune the model**\n\nThank you very much for this suggestion! We experimented with it during the rebuttal period. Specifically, we created a sample of verified model predictions, for polynomial systems, added it to the original training sample, and continued training the model.\n\nWe note that adding about 1,000 correct predictions to the 1 million original training sample improves our performance on the “into the wild” test sets (section 6.4), from 11.7 to 13.5 for Poly3, and 9.6 to 11.9% for Poly5 - a 15% increase in accuracy! The performance on other test sets is unaffected. We are still experimenting with this, and will add our results in the updated paper. Thanks again for suggesting this!\n\n**5- Is there a real-world (or \"real math\") application of finding Lyapunov functions?**\n\nYes, practical applications of stability analysis exist in fields as diverse as supply chain, space industry, car flow regulation on highways, chemical processes etc.\n\nA typical engineering application is the stabilization problem. Given a dynamical system that can be acted upon, e.g. hydraulic gates controlling flow on a river, nozzles controlling the trajectory of satellites, we want to select actions that keep the system stable (i.e. water levels remain constant, satellite stays on orbit). Mathematically, this amounts to setting a parameter of the dynamical equations, so that the system has a Lyapunov function.\n\nIn many cases, a global Lyapunov function cannot be found, and one must settle for a local or semi-global solution, which guarantees that the system remains stable so long it is subjected to small perturbations. This is typical in robotics. \nHowever, when a global Lyapunov function can be found, more efficient controls can be designed. This is particularly important in biological systems (see for instance [Agbo Bidi, Almeida, Coron, 2024] for the regulation of mosquitoes population), epidemiology (see  [Adda &Bichara, IJPAM 2012]) or traffic flow regulation (see [Hayat, Piccoli, Truong, SIAP, 2023]).", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "8NSqSXmMbx", "reviewer": "Reviewer_f8Nq", "summary": "In this paper, the authors propose a new method of generating synthetic training samples from random solutions, i.e. Lyapunov functions for dynamical systems in the forms of ordinary differential equations. They demonstrate that the sequence-to-sequence transformer training on such data can produce Lyapunov functions better than both existing algorithmic solvers and humans. More speficially, for polynomial systems, the method proposed in the paper can find Lyapunov functions for 10.1% of the randomly generated polynomial systems where the state-of-art algorithm can only find 0.7%, furthermore, the method can produce a Lyapunov function for 12,7% of randomly generated non-polynomial systems.", "questions": "It is not clear how the forward and backward generation methods are different from what proposed in the past, for example those in  Lample & Charton (2019), Prajna et al., 2002, and Prajna et al., 2005.", "limitations": "As I have pointed out, the paper does clearly describe the difference between what it proposes to some of existing methods, thus make it quite difficult to see the main reasons for the rather remarkable experiments results.", "rating": 5, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "The problem of Lyapunov function problem is an important problem. The performance of the method proposed in the paper is impressive and promising.", "weaknesses": "The novelty of the data generating methods, those described in Sec. 5, is not clearly stated and explained. Currently, it reads as if the results are obtained accidentally.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your interest and your in-depth reading of the paper.\n\n**The novelty of the data generating methods, those described in Sec. 5, is not clearly stated and explained. Currently, it reads as if the results are obtained accidentally**\n\nWe agree, thank you for raising this. We provide a rationale for the generation method we use in the author rebuttal (section B). We will update the paper accordingly.\n\n**It is not clear how the forward and backward generation methods are different from what proposed in the past, for example those in Lample & Charton (2019), Prajna et al., 2002, and Prajna et al., 2005.**\n\nFor relation with Lample and Charton, and the novelty of our approach, this is detailed in  section A of the authors’ rebuttal and will be made precise in the revised version. \n\nConcerning Prajna et al., 2002, and Prajna et al., 2005, our approach is very different: because Prajna et al. design a deterministic mathematical algorithm, there is no notion of data generation: their algorithm is thought of as a tool for mathematicians and engineers given an input system. We use the algorithm they designed (and subsequent more recent versions using the same approach) as a rejector for rejection sampling and as verifier of the Lyapunov functions provided by the model in polynomial cases.\n\n**As I have pointed out, the paper does not clearly describe the difference between what it proposes to some of existing methods, thus make it quite difficult to see the main reasons for the rather remarkable experiments results.**\n\nThank you for pointing this out, we will clarify this in the revised paper. We believe that the main differences with existing methods are the following:\n\nDifferences with Lample and Charton, 2019 - Charton, Hayat, Lample, 2020\n* These papers address problems with known solutions, for which forward generation of large training sets is possible, we address an open problem still unsolved in mathematics. This raises new difficulties : backward generation and OOD generalization. To our knowledge this is the first application of transformers, trained on synthetic data to solve an open problem end to end.\n* We develop a bespoke technique for generating the backward dataset, that mitigates the risks discussed by Yehudah et al\n* We introduce a new approach to prevent performance from dropping when the model is tested on distributions different from the training distribution. This relies on injecting a tiny proposition (0.05%) of forward examples in the backward dataset.\n\nDifferences with Prajna et al. and other Sum-Of-Squares (SOS) approaches\n* On the method side, SOS approaches are deterministic solvers relying on Semidefinite Programming, very different from our approach, we use these or similar algorithms for testing the output of our model.\n* Our methods allow to handle non-polynomial systems and polynomial systems that have no sum-of-square Lyapunov functions\n* On average on the different datasets, our models have better performances and shorter evaluation/inference time.\n\nDifferences with recent neural methods (“Neural Lyapunov Function”, FOSSIL 2.0, LyZNet):\n* Our approach is very different: we use a Transformer and treat the mathematical problem end-to-end symbolically. Neural methods rely on numerical approximation/simulations of the dynamics, provide an implicit solution, and need to be re-trained on every system.\n* Our methods have much higher performances both on polynomial and non-polynomial datasets (see section 6.3 and our new results in the author rebuttal). \n* We provide symbolic expressions for the Lyapunov functions, whereas existing neural methods learn an implicit “black box” Lyapunov function\n* Our method learn and provide global Lyapunov functions on the whole $\\mathbb{R}^{n}$ (rather than semi-global Lyapunov function on a compact set of $\\mathbb{R}^{n}$ which is an easier problem, although still challenging)", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "K68T9C7Kiz", "reviewer": "Reviewer_S4kD", "summary": "This paper presents a method using sequence-to-sequence transformers to find Lyapunov functions for dynamical systems, mainly polynomial systems. Even though it illustrates some new advances in LLMs in solving mathematical problems, the reviewer doubts the mechanism and the necessity of using large models such as transformers to learn Lyapunov functions for simple polynomial systems, as well as the fairness of the comparisons made in the paper. This work could be a very interesting and good workshop paper to show the advances of LLMs, but the reviewer doesn't think there are significant academic contributions to the math community.", "questions": "There are quite a couple of interesting questions that can be asked:\n1. It's pretty shocking and interesting to see that the authors called SOSTOOLS which was proposed and developed more than two decades ago. Similar to the basic Lyapunov analysis explained above, it seems that the authors were unclear about the latest development of the toolbox or numerical methods for finding Lyapunov functions for nonlinear systems. Even in this year's HSCC conference, there are two tool papers on finding Lyapunov functions for general nonlinear systems: [FOSSIL 2.0](https://dl.acm.org/doi/10.1145/3641513.3651398) and [LyZNet](https://dl.acm.org/doi/10.1145/3641513.3650134). Both of them can handle the non-polynomial cases. The authors did a comparison with FOSSIL, but they claimed that they just achieved <5% accuracy. I don't think they give the toolbox a fair enough try.\n\n2. With bullet point #1 being said, this transformer method needs a very long training time and an extremely large amount of data. However, with the aforementioned two toolboxes, for some systems, including the simple polynomial systems shown in Appendix C, valid Lyapunov functions can be found within a few seconds, and they are formally verified with SMT solvers. The same for SOSTOOLS, the obtained Lyapunov functions are guaranteed to be valid ones. On the contrary, the performances of the proposed method on some datasets are not satisfactory, such as Table 2 and Table 5. \n\n3. Either in the control and system community or in the applied mathematics community, people barely solve Lyapunov functions by hand, the same as adding two very large numbers. The comparison with so-called 'human mathematicians' is not reasonable.\n\n4. In line 189, '$e^i$ are diverse enough to span all of $\\mathcal{H}_x$ '. When is it 'diverse' enough to span the space?\n\n5. It's not surprising to see that after adding some data from forward datasets to the backward ones, the performances are much better. It's very likely that there are some issues with the backward datasets, while the forward datasets are guaranteed to be correct with SOS.", "limitations": "The fairness of the comparisons and the correctness of the learned Lyapuonv functions, as discussed above. Also, even though the authors claimed that their methods worked well with non-polynomial system, the comparisons and examples are mainly for polynomial systems.", "rating": 3, "confidence": 5, "soundness": 2, "presentation": 2, "contribution": 1, "strengths": "This paper finds a way to explore or show the capacity of transformers, or LLMs in general, to solve existing mathematical problems, that is, finding the Lyapunov functions for nonlinear dynamical systems in this paper. Also, the backward data generation method is interesting.", "weaknesses": "Even though the authors tried to make the story interesting and appealing, it seems that a deeper understanding on Lyapunov stability analysis is needed when improving this paper in the future. For example, Lyapunov functions are useful for providing stability guarantees for equilibrium points for dynamical systems, but there are no such equilibrium points in the three-body problem (line 29). Please see the introduction of chapter 3, Lyapunov stability in 'Nonlinear systems' by Khalil, H. K., one of the references of this work. \n\nMoreover, it seems that the authors didn't understand the stability concept fully when writing this paper. In line 27, they stated that the goal is \"discovering the Lyapunov functions that control the global stability of dynamical systems\", while in Def. 3.2, they had the definition of 'stable'. Stable and globally stable are different, which should be precise as an academic paper. Also, that's different from asymptotically stable. Meanwhile, 'control' is definitely not the correct wording here.\n\nLast but not least, finding a Lyapunov function (candidate) for a nonlinear system is not a very hard problem now, even for high-dimensional (polynomial) systems, but how to verify it is indeed a valid Lyapunov function satisfying the Lyapunov conditions is the bottleneck, which is difficult to address for non-polynomial or high-dimensional systems. As shown in the paper ['Neural Lyapunov Control'](https://papers.nips.cc/paper_files/paper/2019/hash/2647c1dba23bc0e0f9cdf75339e120d2-Abstract.html), one can easily use a one-hidden feedforward neural network with 6 hidden neurons to learn a Lyapunov function for many 2- or higher-dimensional nonlinear systems which are not necessary to be polynomials. That said, the results shown in this work are not impressive at all, and the Lyapunov function candidates they found using transformers still lack correctness, which needs to be verified against the Lyapunov conditions, given a new nonlinear system.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your appreciation of the paper.\n\n**there are no equilibrium points in the three-body problem**\n\nWe agree, we mention this problem to emphasize the importance of stability in general. This said, there are periodic orbits which, up to a change of variables, are equilibrium points of a non-autonomous system. \n\n**Stable and globally stable are different, which should be precise as an academic paper.**\n\n“Globally stable” meant “has a global Lyapunov function”. This is indeed different from asymptotic stability (LAS or GAS) or exponential stability for which we make no claim, but it coincides with GAS when the largest invariant set of $x\\in\\mathbb{R}^{n} | f(x)\\cdot \\nabla V(x) = 0$ is the equilibrium (LaSalle). The Lyapunov functions we generate are proper, we will add it in our definition and update the paper to achieve a level of rigor comparable to what we submit to math journals.\n\n**finding a Lyapunov function for a nonlinear system is not a very hard problem**\n\nWe disagree, at least for non-polynomial systems: finding a global Lyapunov function for some particular examples may be easy, but in general it is a very hard problem and still an active area of research in mathematics (see for instance [Adda & Bichara, IJPAM, 2012] or [Agbo Bidi, Almeida, Coron, 2024]).\n\n**As shown in Neural Lyapunov Control, one can easily use a one-hidden feedforward neural network [...] to learn a Lyapunov function.**\n\nNeural Lyapunov Function works well on their showcase examples, but often does not converge for the systems we consider. It achieves 2-22% acc. on polynomial test sets, and <1% on non polynomials. This may be because it was designed for stabilization problems, and performs worse when there is no control feedback.\n\n**the Lyapunov function candidates they found using transformers still lack correctness**\n\nAll the Lyapunov functions we find for polynomial systems are guaranteed true using a SOS checker. For non-polynomial systems, following your suggestion, we now use an SMT solver to theoretically ensure correctness (with no change in performance).\n\n**the authors called SOSTOOLS which was proposed more than two decades ago**\n\nSOSTOOLS is still one of the most widely recognized toolboxes to find Lyapunov functions. While it is two decades old, the underlying SDP solvers -at the core of the computation- have evolved and we are using up-to-date SDP solvers.\n\n**The authors were unclear about the latest development [...] for finding Lyapunov functions [...] FOSSIL 2.0 and LyZNet. [...] I don't think they give the toolbox a fair enough try.**\n\nWe did provide a comparison to FOSSIL. Following your comments, we experimented extensively with the very recent FOSSIL 2.0 and LyZNet (see rebuttal PDF). We were able to reproduce their results and we confirm they have very low accuracy on our polynomial and non-polynomial test sets (see the rebuttal PDF). We will share our code and datasets. \n\nWe believe the low accuracy of these models on our test sets is due to the fact that they address a different problem:  discovering **local or semi-global** Lyapunov functions, while we target global Lyapunov functions. Overall, we believe our work is related, but complementary. We will update the discussion accordingly.\n\n**this transformer method needs a very long training time and an extremely large amount of data.**\n\nOur model needs to be trained once from *synthetic* data, before it can be used at very little cost for any system, a common paradigm in machine learning. In contrast, the smaller models in FOSSIL and LyZNet must be retrained for every system.\n\n**with the two toolboxes [...] Lyapunov functions [...] are formally verified with SMT solvers. The same for SOSTOOLS [...] the performances of the proposed method on some datasets are not satisfactory, such as Table 2 and Table 5.**\n\nOur Lyapunov functions are now also formally verified with an SMT solver (this makes little difference in the results, see authors rebuttal).  On Table 2 and 5 (out-of-distribution and discovery) our methods perform better than SOSTOOLS and the toolboxes. \n\n**people barely solve Lyapunov functions by hand. The comparison with so-called 'human mathematicians' is not reasonable**\n\nWhen classical Lyapunov functions and classical techniques (LMI, BMI SOS, etc.) don’t work, mathematicians (in our community as least) do need to find the form of the Lyapunov functions by hand (see for instance [Adda & Bichara, IJPAM, 2012], [Friedrich et al, SICON, 2023] for ODE systems or [Hayat, Shang, JMPA, 2021] for an example although on PDE). Our models may help “guessing” candidate Lyapunov functions in these situations.\n\nSince finding Lyapunov function is often a textbook exercise in control theory courses, comparing our models with master students in such courses seemed to provide a fair evaluation of the hardness of the problems.\n\n**In line 189 [...] when is it 'diverse' enough to span the space?**\n\nWe mean that $e_{i}(x)$, $1 \\leq i \\leq n$ should be a generative family of the $\\mathcal{H}_{x}$. Thanks for the comment, we’ll clarify. \n\n**It's not surprising to see that after adding some data from forward datasets to the backward ones, the performances are much better**\n\nWhat is surprising is the proportion involved: performance on both distributions improve, after adding 50 forward samples to 1 million backward (0.05%). This is an important result, since OOD is a difficult problem in machine learning.\n\n**It's very likely that there are some issues with the backward datasets**\n\nThe backward datasets are guaranteed correct by design: the steps of the process ensure the conditions from Def. 3.2 hold (see section 5.1 and Appendix A.2). \n\n**the comparisons and examples are mainly for polynomial systems**\n\nOur framework can handle both polynomial and non-polynomial systems. We provide more comparisons on polynomial systems because it is a strong baseline with a long history. We will add more comparisons for non-polynomial systems.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:
We thank the reviewers for their comments an suggestions, which helped improve the paper. Here, we cover questions asked by several reviewers, and present the new experiments we ran for this rebuttal. The results are in the appended PDF file.

**A- On novelty, and comparison with Lample and Charton (2019)**

Lample and Charton introduced the backward approach as a fast data generation technique for symbolic integration. They noticed that it creates out-of-distribution (OOD) generalization issues. Models trained on backward data perform badly on forward-generated test sets.
[Yehudah et al, ICML, 2020] suggested that backward-generation, when applied to hard problems, tends to target easier subproblems. This limitation comes on top of OOD issues. 

Ours is the first attempt to actually use a backward approach to train a model to solve an open problem, for which only a very small number of forward examples are available. We make two novel contributions.
* 1-  We introduce complex generation techniques, like the orthogonal decomposition (step 3), to prevent backward generation from collapsing into solving easier subproblems. While these tricks are problem specific, they show that the limitations described by Yehudah et al. can be mitigated.
* 2- We show that OOD generalization issues can be greatly reduced if a tiny number of forward generated examples is added to the backward-generated training set (50 examples in a million, 0.05%). This approach, related to the priming technique (Jelassi et al 2023) for length generalization, is, to our knowledge, novel.

**B- Rationale of backward generation**

Our models are trained on pairs of problems and solutions. When a solver exists, we can use to generate a training set, by computing the solutions of random problems. This is the forward method. This approach fails for open problems, for which solvers do not exist, or are costly and limited to simple cases. The backward method creates a training set by sampling random solutions (Lyapunov function), and deriving associated problems (dynamical systems).

Lyapunov functions must satisfy three conditions (def. 3.2). One depends on the system and two are intrinsic: being strictly minimum in 0, tending to infinity at infinity.

In Step 1, we sample a function V that satisfies the two intrinsic conditions. To do so, we write $V=V_\text{proper}+V_\text{cross}$, where $V_\text{proper}$ belongs to a class with a guaranteed strict minimum in zero – e.g. sums of one-variable functions, positive definite polynomials – and $V_\text{cross}$ belongs to a large class of non-negative functions, valued 0 at the origin, but with no guarantee of a strict minimum. This class of functions is larger than those considered in usual SOS methods.

In Steps 2 to 4, we address the third condition of Def. 3.2:  $\nabla V(x)\cdot f(x)\leq 0$ for any $x\in\mathbb{R}^{n}$. 

A naive solution would be $f(x) = - \nabla V(x)$, but this would greatly reduce the class of systems we generate, and turn the Lyapunov function discovery problem (find $V$ from $f$) into an easier integration problem (find $V$ from $-\nabla V$). This is the point made by Yehadah et al.

To prevent this, we transform the naive solution $f = - \nabla V(x)$ in two ways:
* Step 4: we multiply it, along each coordinate axis by random non-negative functions $h_i^2$. Because, this does not change the sign of $f \cdot \nabla V$, Def. 3.2 still holds.
* Step 2 and 3: we add a random function, orthogonal to $\nabla V$, $\sum^p_{i=1} e_i(x) g_i(x)$, with $e_i$ random vectors orthogonal to $\nabla V$, and $g_i$ random functions. Since $e_i \cdot \nabla V=0$, Def. 3.2 will still hold.

This guarantees that the resulting system $f$ spans a very large set of systems, because any $f$ satisfying $\nabla V(x)\cdot f(x)\leq 0$ can be written as the sum of a function collinear to $\nabla V(x)$ and a function orthogonal to $\nabla V(x)$. This mitigates Yehuda's limitations, by preventing the model from guessing the solution by reversing the generative process. 

**C- Comparison with neural Lyapunov methods**

Reviewer S4kD raised concern about recent neural methods already solving this problem. We consider our approach to be complementary to theirs.

First, we focus on a different problem: proving that a system has a global Lyapunov function (valid in the full space) in symbolic form. Neural Methods, on the other hand, find implicit (black box)  semi-global Lyapunov functions. While finding semi-global functions, and the corresponding region of attraction, is an important problem in many engineering applications, it is a different problem: we aim at proving stability over the full space. 

Second, in our experiments, we note that our methods perform much better than neural Lyapunov models (including the most recent Fossil 2.0 and LyzNet, which we added to our baselines), see the rebuttal PDF. We believe this is due to the fact that neural methods solve a different problem, as underlined above.

Third, the methods are fairly different. We train a model only once, on generated data, and then use it to predict explicit Lyapunov functions. Neural methods need to be retrained on every system considered, and provide implicit “black box” solutions.

**D- New experiments**

* We added an SMT, verifier to theoretically guarantee the output of the models in the non-polynomial case. It confirms the performances of our models.
* We provide an extensive comparison to neural methods such as FOSSIL 2.0 and LyZNet (see pdf)
* We tested a suggestion by reviewer  LBjn (thanks again!): fine-tuning a model on its verified predictions. We notice that the addition of 1000 verified predictions to our training set of 1 million, improves performance on the “into to wild” test sets by about 15%, while not affecting the other test sets. Adding more examples seems to be detrimental, as it decreases the performance on other benchmarks. 

Experiments are ongoing, we will update the paper with the new results.

Meta-review: The paper proposes to apply transformers for finding Lyapunov functions. 

I find that the proposed technique is interesting and intriguing and may find practical applications 
in finding the global Lyapunov function for the given (polynomial) dynamical system, which
is often performed by manual experimentation. 
The authors designed the whole pipeline, including synthetic dataset creation and 
I therefore accept the paper for the NeurIPS conference.

Several crucial points must be incorporated into the final version of the paper:
* address Reviewer S4kD's concerns about clearly stating the problem setup, definitions, 
Moreover, discuss differences from the mentioned existing methods FOSSIL and LyZNet.
* include the additional experiments reported in the rebuttal pdf,
* share the created dataset, benchmark, and training/evaluation pipeline in an open-source environment.
Decision: Accept (poster)