Global Lyapunov functions: a long-standing open problem in mathematics, with symbolic transformers

Thank you for your appreciation of the paper.

there are no equilibrium points in the three-body problem

We 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.

Stable and globally stable are different, which should be precise as an academic paper.

“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.

finding a Lyapunov function for a nonlinear system is not a very hard problem

We 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]).

As shown in Neural Lyapunov Control, one can easily use a one-hidden feedforward neural network [...] to learn a Lyapunov function.

Neural 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.

the Lyapunov function candidates they found using transformers still lack correctness

All 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).

the authors called SOSTOOLS which was proposed more than two decades ago

SOSTOOLS 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.

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.

We 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.

We 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.

this transformer method needs a very long training time and an extremely large amount of data.

Our 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.

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.

Our 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.

people barely solve Lyapunov functions by hand. The comparison with so-called 'human mathematicians' is not reasonable

When 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.

Since 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.

In line 189 [...] when is it 'diverse' enough to span the space?

We 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.

It's not surprising to see that after adding some data from forward datasets to the backward ones, the performances are much better

What 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.

It's very likely that there are some issues with the backward datasets

The 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).

the comparisons and examples are mainly for polynomial systems

Our 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.

Thank you for your interest and your in-depth reading of the paper.

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

We 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.

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.

For 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.

Concerning 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.

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.

Thank 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:

Differences with Lample and Charton, 2019 - Charton, Hayat, Lample, 2020

• 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.
• We develop a bespoke technique for generating the backward dataset, that mitigates the risks discussed by Yehudah et al
• 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.

Differences with Prajna et al. and other Sum-Of-Squares (SOS) approaches

• 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.
• Our methods allow to handle non-polynomial systems and polynomial systems that have no sum-of-square Lyapunov functions
• On average on the different datasets, our models have better performances and shorter evaluation/inference time.

Differences with recent neural methods (“Neural Lyapunov Function”, FOSSIL 2.0, LyZNet):

• 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.
• 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).
• We provide symbolic expressions for the Lyapunov functions, whereas existing neural methods learn an implicit “black box” Lyapunov function
• 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)

Thank you for your review and comments.

Backward generation. it would be better to elaborate on the method a little more.

We 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.

I encourage the authors to discuss the range of the Lyapunov function and 𝑓

Thank 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.

The 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).

Still, 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.

Another 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}$.

Following your suggestion, we’ll add a discussion on the class of functions that can be covered and the limitations in our revised version.

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.

Thank 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).

The 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.

There are four cases:

• The verifier is too restrictive. Polynomial systems may have non-polynomial Lyapunov functions, but SOS checkers always assume sum-of-square Lyapunov functions.
• The verifier aborts because of timeout or memory overflow, while checking a correct solution.
• The verifier aborts because of timeout or memory overflow, while checking a wrong solution, for a system that is globally stable
• The verifier aborts because the system is not globally stable (so no Lyapunov function exists)

The 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).

Discriminating 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.

Discriminating 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.

What do the "pre-defined set of increasing functions" and "pre-defined set of bounded-functions" refer to at [lines 479, 485]?

Currently 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.

Minor comments

Thank you for spotting these. We corrected them, and added the references you suggested.

Thank you for your review and comments. Here are a few answers and clarifications.

The algorithmic approach is not new

We agree. See section A of the authors’ response for a clarification of our novel contributions. Summarizing,

• 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,
• we show that mixing a tiny number of forward examples into the backward training set helps mitigate the OOD problem.

I find all the different dataset variants included to be a bit confusing

We will change this in the revised version, and limit the main paper to four datasets:

• Two forward sets: Fbarrier and FLyap, for two different problems: barrier functions and Lyapunov functions,
• Two backward sets: BPoly and BNonPoly, for polynomial and non-polynomial systems.

Results on other datasets will be moved to the Appendix.

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"

We 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.

In particular, we propose the following rewrite of the end of the discussion:

Our 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.

Questions

1- findylap solves the same number of problems as SOSTOOLs, right?

Yes. 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.

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?

Yes. 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).

The 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).

The 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.

3- are the polynomial systems for FSOSTOOLS homogenous, nonhomogenous, etc?

This 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.

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

Thank 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.

We 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!

5- Is there a real-world (or "real math") application of finding Lyapunov functions?

Yes, practical applications of stability analysis exist in fields as diverse as supply chain, space industry, car flow regulation on highways, chemical processes etc.

A 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.

In 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. However, 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]).