Learning advanced mathematical computations from examples
Year: 2021

Abstract:
Using transformers over large generated datasets, we train models to learn mathematical properties of differential systems, such as local stability, behavior at infinity and controllability. We achieve near perfect prediction of qualitative characteristics, and good approximations of numerical features of the system. This demonstrates that neural networks can learn to perform complex computations, grounded in advanced theory, from examples, without built-in mathematical knowledge.

heading: 1 Introduction
text: Scientists solve problems of mathematics by applying rules and computational methods to the data at hand. These rules are derived from theory, they are taught in schools or implemented in software libraries, and guarantee that a correct solution will be found. Over time, mathematicians have developed a rich set of computational tools that can be applied to many problems, and have been said to be “unreasonably effective” (Wigner, 1960).
Deep learning, on the other hand, learns from examples and solves problems by improving a random initial solution, without relying on domain-related theory and computational rules. Deep networks have proven to be extremely efficient for a large number of tasks, but struggle on relatively simple, rule-driven arithmetic problems (Saxton et al., 2019; Trask et al., 2018; Zaremba and Sutskever, 2014).
Yet, recent studies show that deep learning models can learn complex rules from examples. In natural language processing, models learn to output grammatically correct sentences without prior knowledge of grammar and syntax (Radford et al., 2019), or to automatically map one language into another (Bahdanau et al., 2014; Sutskever et al., 2014). In mathematics, deep learning models have been trained to perform logical inference (Evans et al., 2018), SAT solving (Selsam et al., 2018) or basic arithmetic (Kaiser and Sutskever, 2015). Lample and Charton (2020) showed that transformers can be trained from generated data to perform symbol manipulation tasks, such as function integration and finding formal solutions of ordinary differential equations.
In this paper, we investigate the use of deep learning models for complex mathematical tasks involving both symbolic and numerical computations. We show that models can predict the qualitative and quantitative properties of mathematical objects, without built-in mathematical knowledge. We consider three advanced problems of mathematics: the local stability and controllability of differential systems, and the existence and behavior at infinity of solutions of partial differential equations. All three problems have been widely researched and have many applications outside of pure mathematics. They have known solutions that rely on advanced symbolic and computational techniques, from formal differentiation, Fourier transform, algebraic full-rank conditions, to function evaluation, matrix inversion, and computation of complex eigenvalues. We find that neural networks can solve these problems with a very high accuracy, by simply looking at instances of problems and their solutions, while being totally unaware of the underlying theory. In one of the quantitative problems
∗ Equal contribution, names in alphabetic order.
where several solutions are possible (predicting control feedback matrix), neural networks are even able to predict different solutions that those generated with the mathematical algorithms we used for training.
After reviewing prior applications of deep learning to related areas we introduce the three problems we consider, describe how we generate datasets, and detail how we train our models. Finally, we present our experiments and discuss their results.

heading: 2 Related work
text: Applications of neural networks to differential equations have mainly focused on two themes: numerical approximation and formal resolution. Whereas most differential systems and partial differential equations cannot be solved explicitly, their solutions can be approximated numerically, and neural networks have been used for this purpose (Lagaris et al., 1998; 2000; Lee and Kang, 1990; Rudd, 2013; Sirignano and Spiliopoulos, 2018). This approach relies on the universal approximation theorem, that states that any continuous function can be approximated by a neural network with one hidden layer over a wide range of activation functions (Cybenko, 1989; Hornik et al., 1990; Hornik, 1991; Petersen and Voigtlaender, 2018; Pinkus, 1999). This has proven to be especially efficient for high dimensional problems.
For formal resolution, Lample and Charton (2020) proposed several approaches to generate arbitrarily large datasets of functions with their integrals, and ordinary differential equations with their solutions. They found that a transformer model (Vaswani et al., 2017) trained on millions of examples could outperform state-of-the-art symbolic frameworks such as Mathematica or MATLAB (Wolfram-Research, 2019; MathWorks, 2019) on a particular subset of equations. Their model was used to guess solutions, while verification (arguably a simpler task) was left to a symbolic framework (Meurer et al., 2017). Arabshahi et al. (2018a;b) proposed to use neural networks to verify the solutions of differential equations, and found that Tree-LSTMs (Tai et al., 2015) were better than sequential LSTMs (Hochreiter and Schmidhuber, 1997) at generalizing beyond the training distribution.
Other approaches investigated the capacity of neural networks to perform arithmetic operations (Kaiser and Sutskever, 2015; Saxton et al., 2019; Trask et al., 2018) or to run short computer programs (Zaremba and Sutskever, 2014). More recently, Saxton et al. (2019) found that neural networks were good at solving arithmetic problems or at performing operations such as differentiation or polynomial expansion, but struggled on tasks like prime number decomposition or on primality tests that require a significant number of steps to compute. Unlike the questions considered here, most of those problems can be solved by simple algorithmic computations.

heading: 3 Differential systems and their stability
text: A differential system of degree n is a system of n equations of n variables x1(t), ..., xn(t),
dxi(t)
dt = fi
( x1(t), x2(t), ..., xn(t) ) , for i = 1...n
or, in vector form, with x ∈ Rn and f : Rn → Rn, dx(t)
dt = f
( x(t) ) Many problems can be set as differential systems. Special cases include n-th order ordinary differential equations (letting x1 = y, x2 = y′, ... xn = y(n−1)), systems of coupled differential equations, and some particular partial differential equations (separable equations or equations with characteristics). Differential systems are one of the most studied areas of mathematical sciences. They are found in physics, mechanics, chemistry, biology, and economics as well as in pure mathematics. Most differential systems have no explicit solution. Therefore, mathematicians have studied the properties of their solutions, and first and foremost their stability, a notion of paramount importance in many engineering applications.

heading: 3.1 Local stability
text: Let xe ∈ Rn be an equilibrium point, that is, f(xe) = 0. If all solutions x(t) converge to xe when their initial positions x(0) at t = 0 are close enough, the equilibrium is said to be locally stable (see Appendix B for a proper mathematical definition). This problem is well known, if f is differentiable in xe, an answer is provided by the Spectral Mapping Theorem (SMT) (Coron, 2007, Theorem 10.10): Theorem 3.1. Let J(f)(xe) be the Jacobian matrix of f in xe (the matrix of its partial derivatives relative to its variables). Let λ be the largest real part of its eigenvalues. If λ is positive, xe is an unstable equilibrium. If λ is negative, then xe is a locally stable equilibrium.
Predicting the stability of a given system at a point xe is our first problem. We will also predict λ, which represents the speed of convergence when negative, in a second experiment. Therefore, to apply the SMT, we need to:
1. differentiate each function with respect to each variable, obtain the formal Jacobian J(x)
f(x) =
( cos(x2)− 1− sin(x1)
x21 − √ 1 + x2
) , J(x) = ( − cos(x1) − sin(x2)
2x1 −(2 √ 1 + x2) −1 ) 2. evaluate J(xe), the Jacobian in xe (a real or complex matrix)
xe = (0.1, ...0.1) ∈ Rn, J(xe) = ( − cos(0.1) − sin(0.1)
0.2 −(2 √ 1 + 0.1)−1
) ,
3. calculate the eigenvalues λi, i = 1...n of J(xe)
λ1 = −1.031, λ2 = −0.441 4. compute λ = −max(Real(λi)) and return the stability (resp. λ the speed of
convergence)
λ = 0.441 > 0→ locally stable with decay rate 0.441

heading: 3.2 Control theory
text: One of the lessons of the spectral mapping theorem is that instability is very common. In fact, unstable systems are plenty in nature (Lagrange points, epidemics, satellite orbits, etc.), and the idea of trying to control them through external variables comes naturally. This is the controllability problem. It has a lot of practical applications, including space launch and the landing on the moon, the US Navy automated pilot, or recently autonomous cars (Bernhard et al., 2017; Minorsky, 1930; Funke et al., 2016). Formally, we are given a system
dx dt = f
( x(t), u(t) ) , (1)
where x ∈ Rn is the state of the system. We want to find a function u(t) ∈ Rp, the control action, such that, beginning from a position x0 at t = 0, we can reach a position x1 at t = T (see Appendix B). The first rigorous mathematical analysis of this problem was given by Maxwell (1868), but a turning point was reached in 1963, when Kalman gave a precise condition for a linear system (Kalman et al., 1963), later adapted to nonlinear system: Theorem 3.2 (Kalman condition). Let A = ∂xf(xe, ue) and B = ∂uf(xe, ue), if
Span{AiBu : u ∈ Rm, i ∈ {0, ..., n− 1}} = Rn, (2)
then the system is locally controllable around x = xe, u = ue.
When this condition holds, a solution to the control problem that makes the system locally stable in xe is u(t) = ue +K(x(t)− xe) (c.f. Coron (2007); Kleinman (1970); Lukes (1968)
and appendix B.4 for key steps of the proof), where K is the m×n control feedback matrix:
K = −Btr ( e−AT [∫ T 0 e−AtBBtre−A trtdt ] e−A trT )−1 . (3)
In the non-autonomous case, where f = f(x, u, t) (and A and B) depends on t, (2) can be replaced by:
Span{Diu : u ∈ Rm, i ∈ {0, ..., 2n− 1} = Rn}, (4)
where D0(t) = B(t) and Di+1(t) = D′i(t) − A(t)Di(t). All these theorems make use of advanced mathematical results, such as the Cayley-Hamilton theorem, or LaSalle invariance principle. Learning them by predicting controllability and computing the control feedback matrix K is our second problem. To measure whether the system is controllable at a point xe, we need to:
1. differentiate the system with respect to its internal variables, obtain A(x, u) 2. differentiate the system with respect to its control variables, obtain B(x, u) 3. evaluate A and B in (xe, ue) 4. calculate the controllability matrix C with (2) (resp. (4) if non-autonomous) 5. calculate the rank d of C, if d = n, the system is controllable 6. (optionally) if d = n, compute the control feedback matrix K with (3)
In: f(x, u) = ( sin(x21) + log(1 + x2) + atan(ux1) 1+x2
x2 − ex1x2
) , xe = [0.1] ue = 1 , Out:  n− d = 0 System is controllable K = (−22.8 44.0)
A step by step derivation of this example is given in Section A of the appendix.

heading: 3.3 Stability of partial differential equations using Fourier Transform
text: Partial Differential Equations (PDEs) naturally appear when studying continuous phenomena (e.g. sound, electromagnetism, gravitation). Over such problems, ordinary differential systems are not sufficient. Like differential systems, PDEs seldom have explicit solutions, and studying their stability has many practical applications. It is also a much more difficult subject, where few general theorems exist. We consider linear PDEs of the form
∂tu(t, x) + ∑ |α|≤k aα∂ α x u(t, x) = 0, (5)
where t, x ∈ Rn, and u(t, x) are time, position, and state. α = (α1, ..., αn) ∈ Rn is a multi-index and aα are constants. Famous examples of such problems include the heat equation, transport equations or Schrodinger equation (Evans, 2010). We want to determine whether a solution u(t, x) of (5) exists for a given an initial condition u(0, x) = u0, and if it tends to zero as t→ +∞. This is mathematically answered (see appendix B.4 and Evans (2010); Bahouri et al. (2011) for similar arguments) by: Proposition 3.1. Given u0 ∈ S ′(Rn), the space of tempered distribution, there exists a solution u ∈ S ′(Rn) if there exists a constant C such that
∀ξ ∈ Rn , ũ0(ξ) = 0 or Real(f(ξ)) > C, (6)
where ũ0 is the Fourier transform of u0 and f(ξ) is the Fourier polynomial associated with the differential operator Dx = ∑ |α|≤k aα∂ α x . In addition, if C > 0, this solution u(t, x) goes to zero when t→ +∞.
Learning this proposition and predicting, given an input Dx and u0, whether a solution u exists, if so, whether it vanishes at infinite time, will be our third and last problem.
To predict whether our PDE has a solution under given initial conditions, and determine its behavior at infinity, we need to: find the Fourier polynomial f(ξ) associated to Dx; find the Fourier transform ũ0(ξ) of u0; minimize f(ξ) on F ; output (0,0) if this minimum is infinite, (1,0) is finite and negative, (1,1) if finite and positive. Optionally, output F . A step by step example is given in Appendix A.
In: Dx = 2∂2x0 + 0.5∂ 2 x1 + ∂ 4 x2 − 7∂ 2 x0,x1 − 1.5∂x1∂ 2 x2 ,
Out:(1, 0) → there exists a solution u ; it does not vanish at t→ +∞

heading: 4 Datasets and models
text: To generate datasets, we randomly sample problems and compute their solutions with mathematical software (Virtanen et al., 2020; Meurer et al., 2017) using the techniques described in Section 3. For stability and controllability, we generate differential systems with n equations and n+ q variables (i.e. n random functions, q > 0 for controllability).
Following Lample and Charton (2020), we generate random functions by sampling unarybinary trees, and randomly selecting operators, variables and integers for their internal nodes and leaves. We use +,−,×, /, exp, log, sqrt, sin, cos, tan, sin-1, cos-1, tan-1 as operators, and integers between −10 and 10 as leaves. When generating functions with n+ q variables, we build trees with up to 2(n+ q + 1) operators.
Generated trees are enumerated in prefix order (normal Polish notation) and converted into sequences of tokens compatible with our models. Integers and floating point reals are also represented as sequences: 142 as [INT+, 1, 4, 2], and 0.314 as [FLOAT+, 3, DOT, 1, 4, E, INT-, 1]. A derivation of the size of the problem space is provided in appendix D.4.
Local stability Datasets for local stability include systems with 2 to 6 equations (in equal proportion). Functions that are not differentiable at the equilibrium xe and degenerate systems are discarded. Since many of the operators we use are undefined at zero, setting xe = 0 would result in biasing the dataset by reducing the frequency of operators like division, square root, or logarithms. Instead, we select xe with all coordinates equal to 0.01 (denoted as xe = [0.01]). This is, of course, strictly equivalent mathematically to sampling systems with equilibrium at the origin or at any other point.
When predicting overall stability, since stable systems become exponentially rare as dimension increases, we use rejection sampling to build a balanced dataset with 50% stable systems. When predicting convergence speed, we work from a uniform (i.e. unbalanced) sample. The value of λ at xe is expressed as a floating point decimal rounded to 4 significant digits. For this problem, we generate two datasets with over 50 million systems each.
Control theory Datasets for automonous control include systems with 3 to 6 equations, and 4 to 9 variables (1 to 3 control variables). In the non-autonomous case, we generate systems with 2 or 3 equations. As above, we discard undefined or degenerate systems. We also skip functions with complex Jacobians in xe (since the Jacobian represents local acceleration, one expects its coordinates to be real). We have xe = [0.5] or [0.9].
In the autonomous case, more than 95% of the systems are controllable. When predicting controllability, we use rejection sampling to create a balanced dataset. In the non-autonomous case, we use a uniform sample with 83% controllable cases. Finally, to predict feedback matrices, we restrict generation to controllable systems and express the matrix as a sequence of floating point decimals. All 3 datasets have more than 50 million examples each.
Stability of partial differential equations using Fourier Transform We generate a differential operator (a polynomial in ∂xi) and an initial condition u0. u0 is the product of n functions f(ajxj) with known Fourier transforms, and d operators exp(ibkxk), with 0 ≤ d ≤ 2n and aj , bk ∈ {−100, . . . , 100}. We calculate the existence of solutions, their behavior when t → +∞, and the set of frequencies, and express these three values as a sequence of 2 Booleans and floating point decimals. Our dataset is over 50 million examples.
Models and evaluation In all experiments, we use a transformer architecture with 8 attention heads. We vary the dimension from 64 to 1024, and the number of layers from 1 to 8. We train our models with the Adam optimizer (Kingma and Ba, 2014), a learning rate of 10−4 and the learning rate scheduler in Vaswani et al. (2017), over mini-batches of 1024 examples. Additional information can be found in appendix D.1. Training is performed on 8 V100 GPUs with float16 operations. Our qualitative models (predicting stability, controllability and existence of solutions) were trained for about 12 hours, but accuracies close to the optimal values were reached after about 6 hours. Learning curves for this problem can be found in appendix D.3. On quantitative models, more training time and examples were needed: 76 hours for convergence speed, 73 hours for control matrices.
Evaluation is performed on held-out validation and test sets of 10000 examples. We ensure that validation and test examples are never seen during training (given the size of the problem space, this never happens in practice). Model output is evaluated either by comparing it with the reference solution or using a problem-specific metric.

heading: 5 Experiments

heading: 5.1 Predicting qualitative properties of differential systems
text: In these experiments, the model is given n functions f : Rn+p → R (n ∈ {2, . . . , 6}, p = 0 for stability, p > 0 for controllability) and is trained to predict whether the corresponding system is stable, resp. controllable, at a given point xe. This is a classification problem.
To provide a baseline for our results, we use fastText (Joulin et al., 2016), a state-of-the-art text classification tool, which estimates, using a bag of words model, the probability of a qualitative feature (stability) conditional to the distribution of tokens and of small fixed sequences (N-grams of up to five tokens). Such a model can detect simple correlations between inputs and outputs, such as the impact on stability of the presence of a given operator, or the number of equations in the system. It would also find out obvious solutions, due to the specifics of one problem or glitches in the data generator. FastText was trained over 2 million examples from our dataset (training over larger sets does not improve accuracy).
A 6-layer transformer with 512 dimensions correctly predicts the system stability in 96.4% of the cases. Since the dataset is balanced, random guessing would achieve 50%. FastText achieves 60.6%, demonstrating that whereas some easy cases can be learnt by simple text classifiers, no trivial general solution exists for this dataset. Prediction accuracy decreases with the degree, but remains high even for large systems (Table 1).
For autonomous controllability over a balanced dataset, a 6-layer transformer with 512 dimensions correctly predicts 97.4% of the cases. The FastText baseline is 70.5%, above the 50% chance level. Whereas accuracy increases with model size (dimension and number of layers), even very small models (dimension 64 and only 1 or 2 layers) achieve performance over 80%, above the FastText baseline (Table 2).
For non-autonomous systems, our dataset features systems of degree 2 and 3, 83% controllable. FastText achieves 85.3%, barely above the chance level of 83%. This shows that text classifiers have difficulty handling difficult problems like this one, even in low dimensions. Our model achieves 99.7% accuracy. Again, small models, that would be unsuitable for natural language processing, achieve near perfect accuracy (Table 3).

heading: 5.2 Predicting numerical properties of differential systems
text: Speed of convergence In these experiments, the model is trained to predict λ, the convergence speed to the equilibrium, up to a certain precision. Here, we consider predictions to be correct when they fall within 10% of the ground truth. Further experiments with different levels of precision (2, 3 or 4 decimal digits) are provided in Appendix C.
A model with 8 layers and a dimension of 1024 predicts convergence speed with an accuracy of 86.6% overall. While reasonably good results can be achieved with smaller models, the accuracy decrease quickly when model size falls under a certain value, unlike when qualitative properties were predicted. Table 4 summarizes the results.
Control feedback matrices In these experiments, we train the model (6 layers, 512 dimensions) to predict a feedback matrix ensuring stability of an autonomous system. We use two metrics to evaluate accuracy:
1) prediction within 10% of all coefficients in the target matrix K given by (3) and provided in the training set,
2) verifying that the model outputs a correct feedback matrix K1, i.e. that all eigenvalues in A+BK1 have negative real parts. This makes more mathematical sense, as it verifies that the model provides an actual solution to the control problem (like a differential equation, a feedback control problem can have many different solutions).
Using the first metric, 15.8% of target matrices K are predicted with less than 10% error. Accuracy is 50.0% for systems with 3 equations, but drops fast as systems becomes larger. These results are very low, although well above chance level (<0.0001%). With the second metric (i.e. the one that actually matters mathematically), we achieve 66.5% accuracy, a much better result. Accuracy decreases with system size, but even degree 6 systems, with 1× 6 to 3× 6 feedback matrices, are correctly predicted 41.5% of the time. Therefore, while the model fails to approximate K to a satisfactory level, it does learn to predict correct solutions to the control problem in 66.5% of the cases. This result is very surprising, as it suggests that a mathematical property characterizing feedback matrices might have been learned.

heading: 5.3 Predicting qualitative properties of PDEs
text: In this setting, the model is given a differential operator Dx and an initial condition u0. It is trained to predict if a solution to ∂tu+Dxu = 0 exists and, if so, whether it converges to 0 when t→ +∞. The space dimension (i.e. dimension of x) is between 2 and 6. In a first series of experiments the model is only trained to predict the existence and convergence of solutions. Overall accuracy is 98.4%. In a second series, we introduce an auxiliary task by adding to the output the frequency bounds F of u0. We observe it significantly contributes to the stability of the model with respect to hyper-parameters. In particular, without the auxiliary task, the model is very sensitive to the learning rate scheduling and often fails to converge to something better than random guessing. However, in case of convergence, the model reaches the same overall accuracy, with and without auxiliary task. Table 6 details the results.

heading: 6 Discussion
text: We studied five problems of advanced mathematics from widely researched areas of mathematical analysis. In three of them, we predict qualitative and theoretical features of differential systems. In two, we perform numerical computations. According to mathematical theory, solving these problems requires a combination of advanced techniques, symbolic and numerical, that seem unlikely to be learnable from examples. Yet, our model achieves more than 95% accuracy on all qualitative tasks, and between 65 and 85% on numerical computations.
When working from synthetic data, a question naturally arises about the impact of data generation on the results of experiments. In particular, one might wonder whether the model is exploiting a defect in the generator or a trivial property of the problems that allows for easier solutions. We believe this is very unlikely. First, because our results are consistent over different problems, using datasets generated with different techniques. Second, because a trivial solution would be found by the bag of words model we use as a baseline. And finally, because we build our datasets by direcly sampling problems from a distribution that includes all possible functions (up to the basis operators and the random number generator). This eliminates the biases that can result from sampling special instances or solutions (Yehuda et al., 2020). It also means that the training set is an extremely tiny sample of the whole problem space (over the 50 million examples generated, we did not get a single duplicate).
Learning from very large samples often raises questions about overfitting and generalisation out of the training distribution. Due to the size of the problem space, it is very unlikely that the model could memorize a large number of cases and interpolate between them. Note that because the space of functions from Rn to Rn has infinite dimension, the universal approximation theorem does not apply here. Note also that for some of our problems (e.g. local stability), mathematical theory states that solutions cannot be obtained by simple interpolation. To investigate out-of-distribution generalization, we modified our data generator to produce 10 new test sets for stability prediction. We changed the distribution of operators and variables, and experimented with systems with longer expressions and more equations. Table 7 (see Appendix C.2 for a detailed analysis) summarizes our key results.
Changes in the distribution of operators and variables have very little impact on accuracy, demonstrating that the model can generalize out of the training distribution. Our trained model also performs well on systems with longer expressions than the training data. This is interesting because generalizing to longer sequences is a known limitation of many sequence to sequence architectures. Finally, a model trained on systems with 2 to 5 equations predicts the stability of systems of 6 equations to high accuracy (78%). Being able to generalize to a larger problem space, with one additional variable, is a very surprising result, that tends to confirm that some mathematical properties of differential systems have been learned.
It seems unlikely that the model follows the same mathematical procedure as human solvers. For instance, problems involving more computational steps, such as non-autonomous controllability, do not result in lower accuracy. Also, providing at train time intermediate results that would help a human calculator (frequencies for PDE, or Jacobians for stability) does not improve performance. Understanding how the model finds solutions would be very interesting, as no simpler solutions than the classical mathematical steps are known.
To this effect, we tried to analyze model behavior by looking at the attention heads and the tokens the models focus on when it predicts a specific sequence (following Clark et al. (2019)). Unfortunately, we were not able to extract specific patterns, and found that each head in the model, from the first layer onwards, attends many more tokens than in usual natural language tasks (i.e. attention weights tend to be uniformly distributed). This makes interpretation very difficult.
These results open many perspectives for transformers in fields that need both symbolic and numerical computations. There is even hope that our models could help solve mathematical problems that are still open. On a more practical level, they sometimes provide fast alternatives to classical solvers. The algorithmic complexity of transformer inference and classical algorithms for the problems we consider here is discussed in appendix E.1. However, for the problems in our dataset, the simpler and parallelizable computations used by transformers allow for 10 to 100 times shorter evaluation times (see Appendix E.2).

heading: 7 Conclusion
text: In this paper, we show that by training transformers over generated datasets of mathematical problems, advanced and complex computations can be learned, and qualitative and numerical tasks performed with high accuracy. Our models have no built-in mathematical knowledge, and learn from examples only. However, solving problems with high accuracy does not mean that our models have learned the techniques we use to compute their solutions. Problems such as non-autonomous control involve long and complex chains of computations, which some of the smaller models we used could probably not handle.
Most probably, our models learn shortcuts that allow them to solve specific problems, without having to learn or understand their theoretical background. Such a situation is common in everyday life. Most of us learn and use language without understanding its rules. On many practical subjects, we have tacit knowledge and know more than we can tell (Polanyi and Sen (2009)). This may be the way neural networks learn advanced mathematics. Understanding what these shortcuts are, how neural networks discover them, and how they can impact mathematical practice, is a subject for future research.

heading: A Examples of computations
text: A.1 Step by step example : autonomous control
To measure whether the system
dx1(t)
dt = sin(x21) + log(1 + x2) + atan(ux1) 1 + x2
dx2(t)
dt = x2 − ex1x2 ,
is controllable at a point xe, with asymptotic control ue, using Kalman condition we need to
1. differentiate the system with respect to its internal variables, obtain the Jacobian A(x, u)
A(x, u) =
( 2x1 cos(x 2 1) + u(1+x2)−1
1+u2x21 (1 + x2) −1 − atan(ux1) (1+x2)2
−x2ex1x2 1− x1ex1x2 ) 2. differentiate the system with respect to its control variables, obtain a matrix B(x, u)
B(x, u) =
( x1((1 + u 2x21)(1 + x2)) −1
0 ) 3. evaluate A and B in xe = [0.5], ue = 1
A(xe, ue) = ( 1.50 0.46 −0.64 0.36 ) , B(xe, ue) = ( 0.27 0 ) 4. calculate the controllability matrix given by (2).
C = [B,AB]((xe, ue)) = [( 0.27 0 ) , ( 1.50 0.46 −0.64 0.36 )( 0.27 0 )] = ( 0.27 0.40 0 −0.17 ) 5. output n−d, with d the rank of the controllability matrix, the system is controllable
if n− d = 0
n− rank(C) = 2− 2 = 0 : System is controllable in (xe = [0.5], ue = 1)
6. (optionally) if n− d = 0, compute the control feedback matrix K as in (3)
K = (−22.8 44.0) .
A.2 Step by step example: stability of linear PDE
To find the existence and behavior at infinite time of a solution, given a differential operator Dx and an initial condition u0 we proceed as follows
1. find the Fourier polynomial f(ξ) associated to Dx Dx = 2∂ 2 x0 + 0.5∂ 2 x1 + ∂ 4 x2 − 7∂ 2 x0,x1 − 1.5∂x1∂ 2 x2 ,
f(ξ) = −4πξ20 − πξ21 + 2πξ42 + 14πξ0ξ1 + 3iπξ1ξ22
2. find the Fourier transform ũ0(ξ) of u0 u0(x) = e −3ix2x−10 sin(x0)e 2.5ix1e−x 2 2 ,
ũ0(ξ) = π 3/21[−(2π)−1,(2π)−1](ξ0)δ0(ξ1 − 2.5(2π)−1)e−π
2(ξ2+3(2π) −1)2
3. find the set F of frequency ξ for which ũ0(ξ) 6= 0 F = [−(2π)−1, (2π)−1]× {2.5(2π)−1} × (−∞,+∞)
4. minimize f(ξ) on F minF (f(ξ)) = −22.6
5. output (0,0) if this minimum is infinite, (1,0) is finite and negative, (1,1) if finite and positive. (optionally) output F Out = (1, 0) : there exists a solution u ; it does not vanish at t→ +∞

heading: A.3 Examples of inputs and outputs

heading: A.3.1 Local stability
text: System Speed of convergenceat xe = [0.01] d dtx0 = − x1 atan (8x0x2) + 0.01atan (0.0008) d dtx1 = − cos (9x0) + cos (0.09)
d dtx2 = x0 − √ x1 + x2 − 0.01 + 0.1 √ 2
−1250
 d dtx0 = − 2x2 x0−2x2(x1−5) + 0.182 d dtx1 = (x1 + (x2 − e x1) (tan (x0) + 3)) (log (3) + iπ) +3.0 log (3) + 3.0iπ
d dtx2 = asin ( x0 log ( − 4x1 )) − asin (0.06 + 0.01iπ)
−0.445
 d dtx0 = e x1+e − sin (x0−e2) − 1.01ee − sin (0.01−e2) d dtx1 = 0.06− 6x1
d dtx2 = −201 + x0+2 x20x2
6.0 (locally stable)
 d dtx0 = x2e −x1 sin (x1)− 9.9 · 10−5 d dtx1 = 7.75.10 −4 − e x2 atan (atan (x1)) 4ex2+9 d dtx2 = (x1 − asin (9)) e − x0 log (3)+iπ
− (0.01− asin (9)) e− 0.01 log (3)+iπ
−0.0384
 d dtx0 = − x0(7− 4 √ 7 √ i) 9 − x1 + 0.0178− 0.00111 4 √ 7 √ i d dtx1 = −0.000379 + e − 63 cos ((x2−9) atan (x1))+7 d dtx2 = −x0 − x1 + asin ( cos (x0) + x2 x0 ) −1.55 + 1.32i
3.52.10−11 (locally stable)
A.3.2 Controllability: autonomous systems
Autonomous system Dimension of
uncontrollable space at xe = [0.5], ue = [0.5] dx0 dt = − asin ( x1 9 − 4 tan (cos (10)) 9 ) − asin ( 4 tan (cos (10)) 9 − 0.0556 ) dx1 dt = u− x2 + log ( 10 + tan (x1)u+x0 ) − 2.36
dx2 dt = 2x1 + x2 − 1.5
0 (controllable)
 dx0 dt = u− asin (x0)− 0.5 + π 6 dx1 dt = x0 − x1 + 2x2 + atan (x0)− 1.46
dx2 dt = 5x2 cos (x2) − 2.85
1
 dx0 dt = 6u+ 6x0 − 6x1 x0 dx1 dt = 0.75 + x 2 1 − cos (u− x2)
dx2 dt = −x 2 0 + x0 + log (e x2)− 0.75
2
 dx0 dt = +x0 ( cos ( u x0+2x2 ) + asin (u)x1 ) −0.5 cos ( 1 3 ) − π6 dx1 dt = πx1 4(x2+4) − π36
dx2 dt = 2.5− 108e 0.5 − 12x0x2 + x1 + 108eu
0 (controllable)
 dx0 dt = −10 sin ( 3x0 log (8) − 22 ) − 6.54 dx1 dt = sin ( 9 + −x1−48x2 ) − 1
dx2 dt = 4 tan ( 4x0 u ) − 4 tan (4)
1
A.3.3 Controllability: non-autonomous systems
Non-autonomous system Local controllabilityat xe = [0.5], ue = [0.5] dx0 dt = (x2 − 0.5) e − asin (8) dx1 dt = e t+0.5 − et+x1 + −x1+e x0 u x2 + 1− 2e
dx2 dt = t(x2 − 0.5)
( asin (6) + √ tan (8) ) False  dx0 dt = atan ( √ x2) x0−1 − 2 atan (√ 2 2 ) dx1 dt = − u −√x0x1+3 + x2 + log (x0) + log(2)− 0.5 + (1/(6− √ 2))
dx2 dt = −70t(x0 − 0.5)
False
 dx0 dt = x0+7 sin (x0eu)+3 dx1 dt = − 9x2e − sin ( √ log (x1))
x0 dx2 dt = t+ asin (tx2 + 4)
False
 dx0 dt = 0.5− x2 + tan (x0)− tan (0.5) dx1 dt = t x1(t+cos (x1(t+u)))
− t0.5(t+cos (0.5t+0.25)) dx2 dt = 2.75− x0 (u+ 4)− x0
True
 dx0 dt = u (u− x0 − tan (8)) + 0.5(tan (8)) dx1 dt = − 6t(−2+π2 ) x0x1 − 12t (4− π) dx2 dt = −7(u− 0.5)− 7 tan (log (x2))
+7 tan (log (0.5))
True
A.3.4 Stability of partial differential equations using Fourier transform
PDE ∂tu+Dxu = 0 and initial condition Existence of a solution, u→ 0 at t→ +∞ Dx = 2∂x0 ( 2∂4x0∂ 4 x2 + 3∂ 3 x1 + 3∂ 2 x1 ) u0 = δ0(−18x0)δ0(−62x2)e89ix0−8649x 2 1+89ix1−59ix2 False , False
 Dx = −4∂ 4 x0 − 5∂ 3 x0 − 6∂ 2 x0∂ 2 x1∂ 2 x2 + 3∂ 2 x0∂x1 − 4∂ 6 x1
u0 = (162x0x2) −1 (ei(−25x0+96x2) sin (54x0) sin (3x2)) True , False Dx = ∂x1 ( 4∂5x0∂x1 + 4∂ 2 x0 − 9∂x0∂ 6 x2 +2∂3x1∂ 5 x2 − 4∂ 3 x1∂ 4 x2 − 2∂x2 ) u0 = (33x0) −1 ( e86ix0−56ix1−16x 2 2+87ix2 sin (33x0)
) True , False  Dx = −6∂7x0∂ 2 x2 + ∂ 5 x0∂ 6 x2 − 9∂ 4 x0∂ 2 x1 − 9∂ 4 x0∂ 4 x2 +7∂2x0∂ 6 x2 + 4∂ 2 x0∂ 5 x2 − 6∂ 6 x1
u0 = δ0(88x1)e −2x0(2312x0+15i)
True , True

heading: B Mathematical definitions and theorems

heading: B.1 Notions of stability
text: Let us consider a system dx(t)
dt = f(x(t)). (7)
xe is an attractor, if there exists ρ > 0 such that
|x(0)− xe| < ρ =⇒ lim t→+∞ x(t) = xe. (8)
But, counter intuitive as it may seem, this is not enough for asymptotic stability to take place.
Definition B.1. We say that xe is a locally (asymptotically) stable equilibrium if the two following conditions are satisfied:
(i) xe is a stable point, i.e. for every ε > 0, there exists η > 0 such that
|x(0)− xe| < η =⇒ |x(t)− xe| < ε, ∀ t ≥ 0. (9)
(ii) xe is an attractor, i.e. there exists ρ > 0 such that
|x(0)− xe| < ρ =⇒ lim t→+∞ x(t) = xe. (10)
In fact, the SMT of Subsection 3.1 deals with an even stronger notion of stability, namely the exponential stability defined as follows:
Definition B.2. We say that xe is an exponentially stable equilibrium if xe is locally stable equilibrium and, in addition, there exist ρ > 0, λ > 0, and M > 0 such that
|x(0)− xe| < ρ =⇒ |x(t)| ≤Me−λt|x(0)|.
In this definition, λ is called the exponential convergence rate, which is the quantity predicted in our first task. Of course, if xe is locally exponentially stable it is in addition locally asymptotically stable.

heading: B.2 Controllability
text: We give here a proper mathematical definition of controllability. Let us consider a nonautonomous system
dx(t)
dt = f(x(t), u(t), t), (11)
such that f(xe, ue) = 0.
Definition B.3. Let τ > 0, we say that the nonlinear system (11) is locally controllable at the equilibrium xe in time τ with asymptotic control ue if, for every ε > 0, there exists η > 0 such that, for every (x0, x1) ∈ Rn × Rn with |x0 − xe| ≤ η and |x1 − xe| ≤ η there exists a trajectory (x, u) such that
x(0) = x0, x(τ) = x1
|u(t)− ue| ≤ ε, ∀ t ∈ [0, τ ]. (12)
An interesting remark is that if the system is autonomous, the local controllability does not depend on the time τ considered, which explains that it is not precised in Theorem 3.2.

heading: B.3 Tempered distribution
text: We start by recalling the multi-index notation: let α = (α1, ..., αn) ∈ Nn, x ∈ Rn, and f ∈ C∞(Rn), we denote
xα = xα11 × · · · × xαnn ∂αx f = ∂ α1 x1 . . . ∂ αn xn f.
(13)
α is said to be a multi-index and |α| = ∑n i=1 |αi|. Then we give the definition of the Schwartz functions: Definition B.4. A function φ ∈ C∞ belongs to the Schwartz space S(Rn) if, for any multi-index α and β,
sup x∈Rn
|xα∂βxφ| < +∞. (14)
Finally, we define the space of tempered distributions: Definition B.5. A tempered distribution φ ∈ S ′(Rn) is a linear form u on S(Rn) such that there exists p > 0 and C > 0 such that
|〈u, φ〉| ≤ C ∑
|α|,|β|<p
sup x∈Rn
|xα∂βxφ|, ∀ φ ∈ S(Rn). (15)

heading: B.4 Proofs of theorems

heading: B.4.1 Analysis of Problem 2
text: The proofs of Theorem 3.2, of validity of the feedback matrix given by the expression (3), and of the extension of Theorem 3.2 to the non-autonomous system given by condition (4) can be found in Coron (2007). We give here the key steps of the proof for showing that the matrix K given by (3) is a valid feedback matrix to illustrate the underlying mechanisms:
• Setting V (x(t)) = x(t)trC−1T x(t), where x is solution to x′(t) = f(x, ue+K.(x−xe)), and
CT =
( e−AT [∫ T 0 e−AtBBtre−A trtdt ] e−A trT ) . (16)
• Showing, using the form of CT , that d
dt (V (x(t))) = −|BtrC−1T x(t)| 2 − |Btre−TA tr C−1T x(t)| 2
• Showing that, if for any t ∈ [0, T ], |BtrC−1T x(t)|2 = 0, then for any i ∈ {0, ..., n− 1},
xtrC−1T A iB = 0, ∀ t ∈ [0, T ].
• Deducing from the controllability condition (2), that
x(t)trC−1T = 0, ∀ t ∈ [0, T ].
and therefore from the invertibility of C−1T ,
x(t) = 0, ∀ t ∈ [0, T ].
• Concluding from the previous and LaSalle invariance principle that the system is locally exponentially stable.

heading: B.4.2 Analysis of Problem 3
text: In this section we prove Proposition 3.1. We study the problem ∂tu+ ∑ |α|≤k aα∂ α x u = 0 on R+ × Rn, (17)
with initial condition u(0, ·) = u0 ∈ S ′(Rn), (18)
and we want to find a solution u ∈ C0([0, T ],S ′(Rn)). Denoting ũ the Fourier transform of u with respect to x, the problem is equivalent to
∂tũ(t, ξ) + ∑ |α|≤k aα(iξ) αũ(t, ξ) = 0, (19)
with initial condition ũ0 ∈ S(Rn). As the only derivative now is with respect to time, we can check that
ũ(t, ξ) = ũ0(ξ)e −f(ξ)t, (20) where f(ξ) = ∑ |α|≤k aα(iξ)
α, is a weak solution to (19) belonging to the space C0([0,+∞),D′(Rn)). Indeed, first of all we can check that for any t ∈ [0,+∞), ξ → exp (−f(ξ)t) is a continuous function and ũ0 belongs to S ′(Rn) ⊂ D′(Rn), thus ũ(t, ·) belongs to D′(Rn). Besides, t → e−f(ξ)t is a C∞ function whose derivative in time are of the form P (ξ)e−f(ξ)t where P (ξ) is a polynomial function. ũ is continuous in time and ũ ∈ C0([0,+∞),D′(Rn)). Now we check that it is a weak solution to (19) with initial condition ũ0. Let φ ∈ C∞c ([0,+∞) × Rn) the space of smooth functions with compact support, we have
− 〈ũ, ∂tφ〉+ ∑ |α|≤k aα(iξ) α〈ũ, φ〉+ 〈ũ0, φ〉
=− 〈ũ0, ∂t(e−f(ξ)tφ)〉 − 〈ũ0, f(ξ)e−f(ξ)tφ〉+ 〈ũ0, e−f(ξ)tf(ξ)φ〉+ 〈ũ0, φ〉 =0.
(21)
Hence, u defined by (20) is indeed a weak solution of (19) in C0([0,+∞),D′(Rn)). Now, this does not answer our question as this only tells us that at time t > 0, u(t, ·) ∈ D′(Rn) which is a less regular space than the space of tempered distribution S ′(Rn). In other words, at t = 0, ũ = ũ0 has a higher regularity by being in S ′(Rn) and we would like to know if equation (19) preserves this regularity. This is more than a regularity issue as, if not, one cannot define a solution u as the inverse Fourier Transform of ũ because such function might not exist. Assume now that there exists a constant C such that
∀ξ ∈ Rn , ũ0(ξ) = 0 or Re(f(ξ)) > C. (22)
∀ ξ ∈ Rn, 1supp(ũ0)e −f(ξ)t ≤ e−Ct. (23)
This implies that, for any t > 0, ũ ∈ S ′(Rn). Besides, defining for any p ∈ N, Np(φ) = ∑
|α|,|β|<p
sup ξ∈Rn
|ξα∂βξ φ(ξ)|, (24)
then for t1, t2 ∈ [0, T ], Np((e−f(ξ)t1 − e−f(ξ)t2)φ) = ∑
|α|,|β|<p
sup ξ∈Rn
|ξαPβ(ξ, φ)|, (25)
where Pβ(ξ, φ) is polynomial with f(ξ), φ(ξ), and their derivatives of order strictly smaller than p. Besides, each term of this polynomial tend to 0 when t1 tends to t2 on supp(ũ0), the set of frequency of u0. Indeed, let β1 be a multi-index, k ∈ N, and Qi(ξ) be polynomials in ξ, where i ∈ {0, ..., k}.∣∣∣∣∣1supp(u0)∂β1ξ φ(ξ) ( k∑ i=0 Qi(ξ)t i 1e −f(ξ)t1 −Qi(ξ)ti2e−f(ξ)t2
)∣∣∣∣∣ ≤
k∑ i=0 max supp(ũ0) ∣∣∣ti1e−f(ξ)t1 − ti2e−f(ξ)t2∣∣∣max ξ∈Rn ∣∣∣∂β1ξ φ(ξ)Qi(ξ, t)∣∣∣ . (26)
From (22), the time-dependant terms in the right-hand sides converge to 0 when t1 tends to t2. This implies that u ∈ C0([0, T ],S ′(Rn)). Finally let us show the property of the behavior at infinity. Assume that C > 0, one has, for any φ ∈ S(Rn)
〈ũ(t, ·), φ〉 = 〈ũ0,1supp(ũ0)e −f(ξ)tφ〉. (27)
Let us set g(ξ) = e−f(ξ)tφ(ξ), one has for two multi-index α and β
|ξα∂βξ g(ξ)| ≤ |ξ αQ(ξ)e−f(ξ)t|, (28)
where Q is a sum of polynomials, each multiplied by φ(ξ) or one of its derivatives. Thus ξαQ(ξ) belongs to S(Rn) and therefore, from assumption (22),
|ξα∂βξ g(ξ)|1supp(u0) ≤ max ξ∈Rn |ξαQ(ξ)|e−Ct, (29)
which goes to 0 when t→ +∞. This imply that ũ(t, ·)→ 0 in S ′(Rn) when t→ +∞, and hence u(t, ·)→ 0. This ends the proof of Proposition 3.1.
Let us note that one could try to find solutions with lower regularity, where u is a distribution of D′(R+ × Rn), and satisfies the equation
∂tu+ ∑ |α|≤k aα∂ α x u = δt=0u0 on R+ × Rn. (30)
This could be done using for instance Malgrange-Erhenpreis theorem, however, studying the behavior at t→ +∞ may be harder mathematically, hence this approach was not considered in this paper.

heading: C Additional experiments
text: C.1 Prediction of speed of convergence with higher precision
In Section 5.1, λ is predicted with a 10% margin error. Prediction of λ to better accuracy can be achieved by training models on data rounded to 2, 3 or 4 significant digits, and measuring the number of exact predictions on the test sample. Overall, we predict λ with two significant digits in 59.2% of test cases. Table 8 summarizes the results for different precisions (for transformers with 6 layers and a dimensionality of 512).

heading: C.2 Out-of-distribution generalization
text: In all our experiments, trained models are tested on held-out samples generated with the same procedure as the training data, and our results prove that the model can generalize out of the training data. However, training and test data come from the same statistical distribution (iid). This would not happen in practical cases: problems would come from some unknown distribution over problem space. Therefore, it is interesting to investigate how the model performs when the test set follows a different statistical distribution. This provides insight about how learned properties generalize, and may indicate specific cases over which the model struggles.
To this purpose, we modified the data generator to produce new test datasets for end to end stability prediction (section 5.1). Four modifications were considered:
1. Unary operators: varying the distribution of operators in the system. In the training data, unary operators are selected at random from a set of nine, three trigonometric functions, three inverse trigonometric functions, logarithm and exponential, and square root (the four basic operations are always present). In this set of experiments, we generated four test sets, without trigonometric functions, without logs and exponentials, only with square roots, and with a different balance of operators (mostly square roots).
2. Variables and integers: varying the distribution of variables in the system. In the training data, 30% of the leaves are numbers, the rest variables. We changed this probability to 0.1, 0.5 and 0.7. This has no impact on expression length, but higher probabilities make the Jacobians more sparse.
3. Expression lengths: making expressions longer than in the train set. In the training data, for a system of n equations, we generate functions with 3 to 2n+ 3 operators. In this experiments, we tried functions between n+ 3 and 3n+ 3 and 2n+ 3 and 4n+ 3. This means that the test sequences are, on average, much longer that those seen at training, a known weakness of sequence to sequence models.
4. Larger degree: our models were trained on systems with 2 to 5 equations, we tried to test it on systems with 6 equations. Again, this usually proves difficult for transformers.
Note that the two first sets of experiments feature out-of-distribution tests, exploring different distributions over the same problem space as the training data. The two last sets, on the other hand, explore a different problem space, featuring longer sequences.
Table 9 presents the results of these experiments. Changing the distribution of operators, variables and integers has little impact on accuracy, up to two limiting cases. First, over systems of degree five (the largest in our set, and more difficult for the transformers) change in operator distribution has a small adverse impact on performance (but not change in variable distribution). Second, which the proportion of integers become very large, and therefore Jacobians become very sparse, the degree of the systems has less impact on performance. But overall results remain over 95%, and the model proves to be very resistant to changes in distribution over the same problem space.
Over systems with longer expressions, overall accuracy tends to decreases. Yet, systems of two or three equations are not affected by a doubling of the number of operators (and sequence length), compared to the training data. Most of the loss in performance concentrates on larger degrees, which suggests that it results from the fact that the transformer is presented at test time with much longer sequences that what it saw at training. In any case, all results but one are well above the fastText baseline (60.5%).
When tested on systems with six equations, the trained model predicts stability in 78.7% of cases. This is a very interesting result, where the model is extrapolating out of the problem space (i.e. no system of six equations have been seen during training) with an accuracy well above chance level, and the fastText baseline.

heading: D Model and problem space

heading: D.1 Model architecture
text: The training loss is the cross entropy between the model predicted output and actual result from the dataset. During training, we use the Adam optimizer, with a learning rate of 0.0001 and scheduling (as in Vaswani et al. (2017)). Mini-batch size varies from one problem to the other, typically between 32 and 128 examples.
During training, we use 8 GPU. The model is distributed across GPUs, so that all of them have access to the same shared copy of the model. At each iteration, every GPU processes an independently generated batch, and the optimizer updated the model weights using the gradients accumulated by all GPU. Overall, this is equivalent to training on a single GPU, but with 8 times larger batches.
D.2 Model behavior and attentions heads
We tried to analyze model behavior by looking at the attention heads and the tokens the models focus on when it predicts a specific sequence. As each head attends many more tokens than in usual natural language tasks, and to improve visualization, we tried to reduce the number of hidden states a head can attend by using a top-k on the attention weights, but this deteriorated the performance, and we did not investigate more in this direction. We also ran a sequence-to-sequence model without attention, so that each input equation is mapped to a fixed-sized representation. We then fed a set of input equations into the model, and used a t-SNE visualization to see whether we can observe clusters of equations. What we observed is mainly that equations with nearby representations have similar length / tokens. However, even embeddings in similar locations can lead to different decoded sequences. The relevance of the representations built in the encoder depends on how the computation is split between the encoder and the decoder. If the decoder does the majority of the work, encoder representations become less meaningful.

heading: D.3 Learning curves
text: Although all generated datasets included more than 50 million examples, most models were trained on less. Figure 1 shows how performance increases with the number of training examples, for the end to end stability problem (i.e. predicting whether systems of degree 2 to 5 are stable). There are twelve curves corresponding to as many experiments over shuffled versions of the dataset (i.e. different experiments used different parts of the dataset).
Overall, less than 10 million examples are needed to achieve close to optimal accuracy. Learning curves from different experiments are close, which proves the stability of the learning process.

heading: D.4 Size of the problem space
text: Lample and Charton (2020) provide the following formula to calculate the number of functions with m operators:
E0 = L
E1 = (q1 + q2L)L
(m+ 1)Em = (q1 + 2q2L)(2m− 1)Em−1 − q1(m− 2)Em−2 Where L is the number of possible leaves (integers or variables), and q1 and q2 the number of unary and binary operators. In the stability and controllability problems, we have q1 = 9, q2 = 4 and L = 20 + q, with q the number of variables.
Replacing, we have, for a function with q variables and m operators E0(q) = 20 + q
E1(q) = (89 + 4q)(20 + q)
(m+ 1)Em(q) = (169 + 8q)(2m− 1)Em−1 − 4(m− 2)Em−2
In the stability problem, we sampled systems of n functions, with n variables, n from 2 to 6. Functions have between 3 and 2n+ 2 operators. The number of possible systems is
PSst = 6∑ n=2 ( 2n+2∑ m=3 Em(n) )n > E14(6) 6 ≈ 3.10212
(since Em(n) increases exponentially with m and n, the dominant factor in the sum is the term with largest m and n)
In the autonomous controllability problem, we generated systems with n functions (n between 3 and 6), and n + p variables (p between 1 and n/2). Functions had between n + p and 2n+ 2p+ 2 operators. The number of systems is
PSaut = 6∑ n=3 n/2∑ p=1 2(n+p+1)∑ m=n+p Em(n+ p) n > E20(9)6 ≈ 4.10310 For the non-autonomous case, the number of variables in n+ p+ 1, n is between 2 and 3 and p = 1, therefore
PSnaut = 3∑ n=2  2(n+2)∑ m=n+1 Em(n+ 2) n > E10(5)3 ≈ 5.1074 Because expressions with undefinite or degenerate jacobians are skipped, the actual problem space size will be smaller by several orders of magnitude. Yet, problem space remains large enough for overfitting by memorizing problems and solutions to be impossible.

heading: E Computation efficiency

heading: E.1 Algorithmic complexity
text: Let n be the system degree, p the number of variables and q the average length (in tokens) of functions in the system. In all problems considered here, we have p = O(n). Differentiating or evaluating an expression with q tokens is O(q), and calculating the Jacobian of our system is O(npq), i.e. O(n2q).
In the stability experiment, calculating the eigenvalues of the Jacobian will be O(n3) in most practical situations. In the autonomous controllability experiments, construction of the n×np Kalman matrix is O(n3p), and computing its rank, via singular value decomposition or any equivalent algorithm, will be O(n3p) as well. The same complexity arise for feedback matrix computations (multiplication, exponentiation and inversion are all O(n3) for a square n matrix). As a result, for controllability, complexity is O(n4). Overall, the classical algorithms have a complexity of O(n2q) for Jacobian calculation, and O(n3) (stability) and O(n4) (controllability) for the problem specific computations.
Current transformer architectures are quadratic in the length of the sequence, in our case nq, so a transformer will be O(n2q2) (in speed and memory usage). Therefore, the final comparison will depend on how q, the average length of equations, varies with n, the number of parameters. If q = O(1) or O(log(n)), transformers have a large advantage over classical methods. This means sparse Jacobians, a condition often met in practice. For controllability, the advantage remains if q = O(n1/2), and the two methods are asymptotically equivalent if q = O(n).
However, current research is working on improving transformer complexity to log-linear or linear. If this happened (and there seem to be no theoretical reason preventing it), transformers would have lower asymptotic complexity in all cases.
E.2 Computation time versus evaluation time
Table 10 compares the average time needed to solve one problem, for a trained transformer running on a GPU, and a Python implementation of the algorithms, running on a MacBook Pro.

Review summaries:
This paper shows that transformer models can be used to accurately learn advanced mathematical computations from millions of examples.  The problems are drawn from the fields of differential equations and control theory.  The selected problems are ones that are solvable using known algorithms; however, these algorithms involve a sequence of advanced mathematical operations (e.g., differentiation, calculating the rank of a matrix, calculating the eigenvalues of a matrix, etc), for which no known simple shortcuts exist. For the experiments in this paper, for each problem a large number (50 million) of training examples are randomly generated, and are then used to train a transformer model.  Across these problems, the paper shows that the neural network is able to solve these problems at high accuracy (96-99.7% accuracy).