Maxime Bouchereau

Philippe Chartier

Mohammed Lemou

Florian Méhats

Machine Learning Methods for Autonomous Ordinary Differential Equations

Keywords: modified equation, ordinary differential equation, neural network, numerical method, convergence analysis

Ordinary Differential Equations are generally too complex to be solved analytically. Approximations thereof can be obtained by general purpose numerical methods. However, even though accurate schemes have been developed, they remain computationally expensive: In this paper, we resort to the theory of modified equations in order to obtain "on the fly" cheap numerical approximations. The recipe consists in approximating, prior to that, the modified field associated to the modified equation by neural networks. Elementary convergence results are then established and the efficiency of the technique is demonstrated on experiments.

Introduction

Ordinary Differential Equations are ubiquitous in the modelling of systems in domains of science as diverse as biology, dynamics, fluid mechanics, quantum mechanics, thermodynamics or weather forecasting. In most situations, these differential equations can not be solved analytically and necessitate the use of numerical methods so as to compute accurate approximations [START_REF] Hairer | Geometric Numerical integration: structure-preserving algorithms for ordinary differential equations[END_REF][START_REF] Hairer | Solving ordinary differential equations[END_REF].

Generally speaking, the aforementioned methods are the result of a trade-off between accuracy and computational cost. In simple words, a small approximation error requires long computations. Of course, there are various ways in which one may soften the computational constraints, e.g. by 1 raising the order, taking into account the structure of the problem or optimising the coefficients of the method. Excellent numerical methods are abundant in the literature and we refer to reference books [START_REF] Beyn | Current Challenges in Stability Issues for Numerical Differential Equations[END_REF][START_REF] Butcher | Numerical methods for ordinary differential equations[END_REF][START_REF]Advances in Differential Equations and Applications[END_REF][START_REF] Hairer | Geometric Numerical integration: structure-preserving algorithms for ordinary differential equations[END_REF][START_REF] Hairer | Solving ordinary differential equations[END_REF][START_REF] Leimkuhler | Simulating hamiltonian dynamics[END_REF] for them.

In this paper, we choose to derive as much as possible prior information from the knowledge of the vector field of the equation in order to accelerate the solving of the equation. This is thus a two-step process: (i) first, generate data that will be used during the effective computation of the approximate solution; (ii) second, compute, from an initial value, an approximation of the solution "on the fly" as accurately as possible by using available pre-computed values. As such, the process remains too vague to become practical: this is where the crucial ingredient of the technique comes into play, that is to say modified equations.

The theory of modified equations has emerged in the context of ordinary differential equations as a powerful tool (referred to as backward error analysis) to explain the excellent behaviour of structure-preserving numerical methods [START_REF] Hairer | Geometric Numerical integration: structure-preserving algorithms for ordinary differential equations[END_REF]. More recently, modified fields have also been used in a dual manner, as a technique to raise the order of any existing numerical method by twisting appropriately the original vector field [START_REF] Beyn | Current Challenges in Stability Issues for Numerical Differential Equations[END_REF]. The idea is that, by adding ad-hoc perturbation terms to the vector field and then solving the associated differential equation, one may compute higher-order approximations. It can be indeed proved that, at least at a formal level, for any numerical method, there exists a modified equation whose solution by the aforementioned method coincides with the exact solution. The perturbation terms to be added may be obtained in analytical form as elementary differentials involving various derivatives of the original vector field. Complete expansions are for instance available with the help of representations by trees (known as B-series). However, computing the corresponding expressions analytically would be an extremely tedious process and this is the reason why the technique has remained confined to specific situations [START_REF] Chartier | Numerical integrators based on modified differential equations[END_REF][START_REF] Hairer | Geometric Numerical integration: structure-preserving algorithms for ordinary differential equations[END_REF][START_REF] Hairer | Preprocessed discrete Moser-Veselov algorithm for the full dynamics of a rigid body[END_REF] of limited practical interest.

The main idea of this work is thus to combine the theory of modified equations with machine learning techniques and more precisely neural networks. Given a differential equation and a numerical method, the ad-hoc perturbations of the vector field are first learnt by extensive simulations and then approximated by inference from a neural network. Once this representation of the modified vector field has been obtained, it is used to solve the original equation with the same numerical method for any initial value prescribed in a learnt domain from the phase-space. The combined computational work ids by far greater than for any usual reasonable numerical method. Nevertheless, if one omits the time spend to learn the perturbation, an accurate solution can be obtained very cheaply as compared to well-established schemes such as those of Dormand & Prince [START_REF] Prince | High order embedded Runge-Kutta formulae[END_REF].

Scope of the paper

The article is divided into two main sections: Section 2 is devoted to the exposition of the technique and its convergence analysis, while Section 3 presents numerical experiments illustrating the performances and properties of the schemes we analysed. Subsection 2.1 exposes the general strategy which is adopted. A specific structure of the neural network is selected, in agreement with the structure of the modified field. Moreover, the details of the machine leaning method for the learning of the modified field are given here, concerning in particular the choice of the training data set and the calibration of the parameters of the neural network (approximating the modified field via loss-minimization).

Subsection 2.2 establishes a convergence result for the numerical integration resulting from the combined use of neural networks and classical schemes. The essential difference with standard convergence results is reflected in the multiplicative constant of the local error which turns out to be smaller than for a direct application to the original vector field. A special focus is put on explicit Runge-Kutta methods in the same subsection, where more precise estimates are given. Moreover, two specific schemes are studied in this subsection: the forward Euler method and a Runge-Kutta method of order 2, which are two simple occurrences of explicit methods. A short convergence analysis is undertaken for each of them, leading to improved bounds. The same methods are then used for numerical experiments which illustrate two main contributions of this work: (i) The modified field can be learnt efficiently through a neural network, as is illustrated in 3.1; (ii) The resulting numerical methods have far greater efficiency, as is illustrated on convergence curves in Subsection 3.2. A comparison with the well-known DOPRI5 method in Subsection 3.3 is particularly enlightening with this respect.

Related work

The link between differential equations and machine learning has been already explored in several publications. Two approaches are prominent. On the one hand, the ODE vector field can be learnt by the technique of MSE Loss, which can be applied to both ODEs or PDEs [START_REF] Raissi | Multistep neural networks for data-driven discovery of nonlinear dynamical systems[END_REF]. Let us notice also a paper of Burton et al. [START_REF] Brunton | Discovering governing equations from data by sparse identification of nonlinear dynamical systems[END_REF] which proposes symbolic regression for the learning of complex dynamical systems. On the other hand, the ODE vector field can be learnt by statistical methods. for instance, Expectation-Maximization is used in a paper of Nguyen et al. [START_REF] Nguyen | Em-like learning chaotic dynamics from noisy and partial observations[END_REF], while a paper of Raissi et al. [START_REF] Raissi | Machine learning of linear differential equations using Gaussian processes[END_REF] proposes Gaussian processes for linear differential equations.

Links with modified equations. Links between machine learning and modified equations have been established more recently, e.g. in the paper [START_REF] Zhu | On Numerical Integration in Neural Ordinary Differential Equations[END_REF] where the theory of modified equation is used fora rigorous analysis. Offen et al. [START_REF] Offen | Symplectic integration of learned Hamiltonian systems[END_REF] consider numerical methods for Hamiltonian ODEs. The methods used therein are statistical methods, namely Gaussian processes.

Structure of Neural Networks. In order to preserve geometric properties, specific neural networks, adapted to the structure of the differential equations under consideration have been developed. A first example are Hamiltonian equations for which Hamiltonian neural networks have been developed (see [START_REF] David | Symplectic learning for Hamiltonian neural networks[END_REF] or [START_REF] Jin | SympNets: Intrinsic structurepreserving symplectic networks for identifying Hamiltonian systems[END_REF] where irregular time observed data can be used). A second example are Poisson systems for which Jin et al. [START_REF] Jin | Learning Poisson systems and trajectories of autonomous systems via Poisson neural networks[END_REF] developed Poisson neural networks. Another example of specific ODEs is studied in another recent paper [START_REF] Zhu | VPNets: Volume-preserving neural networks for learning source-free dynamics[END_REF], where VPNets are used to learn the volume-preserving flows.

ODEs methods for Neural Networks. In the same way as neural networks are used to solve ODEs, the reciprocal strategy can be pursued: neural networks can indeed be modelled by ODEs and their properties deduced from the corresponding ODE properties. For instance, Lu et al. [START_REF] Lu | Beyond finite layer neural networks: Bridging deep architectures and numerical differential equations[END_REF] have developed new neural networks which are discretizations of ODEs by various numerical methods and Haber et al. [START_REF] Haber | Stable architectures for deep neural networks[END_REF] have developed new structures of neural networks depending on the stability properties of the ODEs. Finally, Chen et al. [START_REF] Chen | Neural ordinary differential equations[END_REF] have derived new optimization methods of the loss function from the properties of the corresponding ODEs (the neural network is here again obtained through the discretization of the ODE).

Approximation by neural networks. In order to properly approximate functions by neural networks, error estimates have been established. Anastassiou [START_REF] Anastassiou | Quantitative approximations[END_REF] stated rates of convergence for approximations of functions by networks, according to the number of parameters and the dimension. By considering neural network spaces as functional spaces, Gribonval et al. [START_REF] Gribonval | Approximation spaces of deep neural networks[END_REF] have obtained inclusions of theses spaces in Besov or Lebesgue spaces, according to the number of parameters and a given rate of convergence. In a separate work, Bach [START_REF] Bach | Learning Theory from First Principles[END_REF] has given bounds on the approximation error in a Hilbertian setting. Finally, a problem of approximation [START_REF] Champagne | Les réseaux de neurones multi-couches, le comment et le pourquoi[END_REF] which is underlined is the curse of dimensionality, where high-dimensional vector fields are approximated with a slower rate of convergence than low-dimensional vector fields, i.e. for high dimensions, more parameters and more data will be required in order to get a satisfying learning.

Improving the accuracy of numerical methods with machine learning

Consider an autonomous ordinary differential equation of the form

𝑦(𝑡) = 𝑓 (𝑦(𝑡)) ∈ R 𝑑 , 𝑡 ∈ [0, 𝑇] 𝑦(0) = 𝑦 0 ,
where 𝑓 : R 𝑑 → R 𝑑 is assumed to be smooth enough. By Cauchy-Lipschitz theorem, we have existence and uniqueness of a solution for any given initial value 𝑦 0 ∈ R 𝑑 . We wish to approximate the solution over [0, 𝑇] at times 𝑡 𝑛 = 𝑛ℎ, 0 ⩽ 𝑛 ⩽ 𝑁, where ℎ = 𝑇 𝑁 is the time-step and 𝑁 is the number of discretization points.

General strategy

As explained in the Introduction section, we shall approximate the modified vector field with the help of a neural network.

Modified field

Let us consider Φ 𝑓 ℎ (•) the numerical flow associated to a given numerical method (ℎ is the time-step of the method) and to the vector field 𝑓 , and assume that it is of order 𝑝, in the sense that1

𝑀𝑎𝑥 0⩽𝑛⩽𝑁 Φ 𝑓 ℎ 𝑛 (𝑦 0 ) -𝜑 𝑓 𝑛ℎ (𝑦 0 ) ⩽ 𝐶ℎ 𝑝 (1)
for some constants 𝐶 > 0. If we modify the field 𝑓 used in Φ ℎ , i.e. if we apply the numerical flow Φ ℎ with the modified field fℎ instead of 𝑓 , we may obtain a higher-order approximation. In fact, the theory of modified equations states that it is possible to construct fℎ as a series of powers of ℎ multiplied by appropriately chosen functions (at least as a formal series), in such a way that (Φ fℎ ℎ ) 𝑛 (𝑦 0 ) coincides exactly with 𝑦(𝑡 𝑛 ). The structure of this modified field writes (see [START_REF] Hairer | Geometric Numerical integration: structure-preserving algorithms for ordinary differential equations[END_REF])

fℎ (𝑦) = 𝑓 (𝑦) + ℎ 𝑝 +∞ ∑ 𝑗=1 ℎ 𝑗-1 𝑓 [ 𝑗] (𝑦) = 𝑓 (𝑦) + ℎ 𝑝 𝑘-1 ∑ 𝑗=1 ℎ 𝑗-1 𝑓 [ 𝑗] (𝑦) + ℎ 𝑘+𝑝-1 𝑅(𝑦, ℎ) (2)
1Here and in the sequel, | • | denotes a norm on R 𝑑 .

where the coefficient-functions 𝑓 [ 𝑗] are built upon derivatives of 𝑓 . It can be shown rigorously that the truncation of this formal series (2) obtained by neglecting the O ℎ 𝑘+𝑝-1 -terms, leads to

(Φ fℎ ℎ ) 𝑛 (𝑦 0 ) = 𝜑 𝑓 𝑛ℎ (𝑦 0 ) + O ℎ 𝑘+𝑝-1
where 𝑘 denotes the number of terms kept in fℎ .

Machine learning methods

The main idea of this paper consists in approximating the modified field fℎ by a neural network approximation 𝑓 𝑎 𝑝 𝑝 (•, ℎ) whose structure mimics the structure of the theoretical modified field [START_REF] Bach | Learning Theory from First Principles[END_REF]. More precisely, we shall approximate separately each function 𝑓 [ 𝑗] in (2) with a neural network. As could be anticipated, the truncation of (2) will be echoed by a similar truncation

𝑓 𝑎 𝑝 𝑝 (𝑦, ℎ) = 𝑓 (𝑦) + ℎ 𝑝 𝑁 𝑡 -1 ∑ 𝑗=1 ℎ 𝑗-1 𝑓 𝑗 (𝑦) + ℎ 𝑁 𝑡 +𝑝-1 𝑅 𝑎 (𝑦, ℎ), (3)
where 𝑓 𝑖 , for 1 ⩽ 𝑖 ⩽ 𝑁 𝑡 -1, and 𝑅 𝑎 , are Multi Layer Perceptrons. An obvious advantage of this choice is that the corresponding numerical method remains consistent. Let us further notice that learning the perturbation in this way is also well-adapted to the situations where the original vector field possesses a specific structure [START_REF] David | Symplectic learning for Hamiltonian neural networks[END_REF][START_REF] Jin | Learning Poisson systems and trajectories of autonomous systems via Poisson neural networks[END_REF][START_REF] Jin | SympNets: Intrinsic structurepreserving symplectic networks for identifying Hamiltonian systems[END_REF][START_REF] Raissi | Multistep neural networks for data-driven discovery of nonlinear dynamical systems[END_REF][START_REF] Zhu | VPNets: Volume-preserving neural networks for learning source-free dynamics[END_REF].

The complete numerical procedure can be decomposed into three main steps : Firstly, data are collected by simulating very accurately the exact flow at various points of the domain. A high number of simulations and a high accuracy are prerequisite for a good approximation of the modified field. Secondly, the different neural networks are trained separately by minimising a prescribed Loss-function. Eventually, given an initial data 𝑦 0 , an approximation of the exact solution is obtained by applying the same numerical scheme as the one used for the training to the field 𝑓 𝑎 𝑝 𝑝 . In more details, we follow the three stages:

1. Construction of the data set: 𝐾 initial data 𝑦 (𝑘) 0 are randomly selected into a compact set Ω ⊂ R 𝑑 (where we wish to simulate the solution) with uniform distribution. Then, for all 0 ⩽ 𝑘 ⩽ 𝐾 -1, we compute a very accurate approximation of the the exact flow at times ℎ (𝑘) with initial condition 𝑦 (𝑘) 0 , denoted 𝑦 (𝑘) 1 . Time steps ℎ (𝑘) are selected in the domain [ℎ -, ℎ + ] (we actually pick up the value log ℎ (𝑘) randomly in the domain [log ℎ -, log ℎ + ] with uniform distribution).

Training of the neural networks:

We minimize the Mean Squared Error (MSE), denoted 𝐿𝑜𝑠𝑠 𝑇𝑟𝑎𝑖𝑛 , which measures the difference between predicted data ŷ(𝑘,ℓ) 1 and "exact data" 𝑦 (𝑘,ℓ) 1 by computing the optimal NN's parameters over 𝐾 0 data (where 1 ⩽ 𝐾 0 ⩽ 𝐾 -1) by resorting to a gradient method:

𝐿𝑜𝑠𝑠 𝑇𝑟𝑎𝑖𝑛 = 1 𝐾 0 𝐾 0 -1 ∑ 𝑘=0 1 ℎ (𝑘) 2 𝑝+2 Φ 𝑓 𝑎 𝑝 𝑝 (•,ℎ (𝑘 ) ) ℎ (𝑘 ) 𝑦 (𝑘) 0 = ŷ1 (𝑘 ) -𝜑 𝑓 ℎ (𝑘 ) 𝑦 (𝑘) 0 =𝑦 (𝑘 ) 1 2 (4)
At the same time, we compute the value of another MSE, denoted 𝐿𝑜𝑠𝑠 𝑇 𝑒𝑠𝑡 , which measures the difference between predicted data ŷ(𝑘,ℓ) 1 and "exact data" 𝑦 (𝑘,ℓ) 1 for a subset of the initial values which have not been used to train the NNs. The objective of this step is to estimate the performance of the training for "unknown" initial values :

𝐿𝑜𝑠𝑠 𝑇 𝑒𝑠𝑡 = 1 𝐾 -𝐾 0 𝐾-1 ∑ 𝑘=𝐾 0 1 ℎ (𝑘) 2 𝑝+2 Φ 𝑓 𝑎 𝑝 𝑝 (•,ℎ (𝑘 ) ) ℎ (𝑘 ) 𝑦 (𝑘) 0 -𝜑 𝑓 ℎ (𝑘 ) 𝑦 (𝑘) 0 2 (5)
If 𝐿𝑜𝑠𝑠 𝑇𝑟𝑎𝑖𝑛 and 𝐿𝑜𝑠𝑠 𝑇 𝑒𝑠𝑡 exhibit the same decay pattern, one considers that there is no overfitting, that is to say that the neural network model does not fit exactly against its training data and remains able to perform accurately against unseen data, which is its main purpose.

Numerical approximation:

At the end of the training, an accurate approximation 𝑓 𝑎 𝑝 𝑝 (•, ℎ) of fℎ is available. It is then used to compute the successive values of (Φ

𝑓 𝑎 𝑝 𝑝 (•,ℎ) ℎ
) 𝑛 (𝑦 0 ) for 𝑛 = 0, . . . , 𝑁.

Error analysis

In this subsection, we analyse the error resulting from the procedure described in previous Subsection. More specifically, we state estimates of the global error for any standard numerical method.

Theorem 1. Let us denote

Φ 𝑓 𝑎 𝑝 𝑝 (•,ℎ) ℎ
, the flow of a given numerical scheme Φ ℎ of order 𝑝, applied to the modified field 𝑓 𝑎 𝑝 𝑝 (•, ℎ) and let us consider the global error

𝑒 𝑛 := Φ 𝑓 𝑎 𝑝 𝑝 (•,ℎ) ℎ 𝑛 (𝑦 0 ) -𝜑 𝑓 ℎ 𝑛 (𝑦 0 ), (6)
at times 𝑡 𝑛 = 𝑛ℎ for 0 ⩽ 𝑛 ⩽ 𝑁. Denoting the learning error by

𝛿 := Max (𝑦,ℎ)∈Ω×[ℎ -,ℎ + ] fℎ (𝑦, ℎ) -𝑓 𝑎 𝑝 𝑝 (𝑦, ℎ) ℎ 𝑝 (7)
and assuming that (i) For any pair smooth vector fields 𝑓 1 and 𝑓 2 , we have

∀0 ≤ ℎ ≤ ℎ + , Φ 𝑓 1 ℎ -Φ 𝑓 2 ℎ 𝐿 ∞ (Ω) ⩽ 𝐶ℎ || 𝑓 1 -𝑓 2 || 𝐿 ∞ (Ω) (8)
for some positive constant 𝐶, independent of 𝑓 1 and 𝑓 2 ;

(ii) For any smooth vector field 𝑓 , there exists a constant 𝐿 > 0 such that

∀0 ≤ ℎ ≤ ℎ + , ∀(𝑦 1 , 𝑦 2 ) ∈ Ω 2 , Φ 𝑓 ℎ (𝑦 1 ) -Φ 𝑓 ℎ (𝑦 2 ) ⩽ (1 + 𝐿ℎ) |𝑦 1 -𝑦 2 | . ( 9
)
Then there exist two constants C, L > 0 such that:

𝑀𝑎𝑥 0⩽𝑛⩽𝑁 |𝑒 𝑛 | ⩽ 𝐶𝛿ℎ 𝑝 L 𝑒 L𝑇 -1 (10)
Proof. The arguments of the proof are completely standard ans thus omitted.

■

Remarks. (i)

The vector fields 𝑓 and fℎ are smooth, respectively by assumption and by construction. As for 𝑓 𝑎 𝑝 𝑝 (•, ℎ), it is smooth as well given that it is obtained through the composition of affine functions 𝐴 1 , • • • , 𝐴 𝐿+1 and nonlinear functions Σ 1 , • • • , Σ 𝐿 (the so-called activation functions). The output of the NN thus appears to be of the form

𝐴 𝐿+1 • Σ 𝐿 • 𝐴 𝐿 • • • • Σ 1 • 𝐴 1 (for 𝐿 layers).
Hence, if the activation functions are smooth, then so is 𝑓 𝑎 𝑝 𝑝 (•, ℎ). This is the case for instance if the Σ 𝑖 's are the hyperbolic tangent functions.

(ii) Assumption ( 8) is straightforwardly satisfied for all known consistent methods.

(iii) A similar error estimate holds for a variable step-size implementation of the numerical method Φ: if we indeed use the sequence of steps 0

≤ ℎ 𝑗 ≤ ℎ + , then 𝑇 = ℎ 0 + • • • , +ℎ 𝑁-1 and 𝑦 * 𝑛+1 = Φ 𝑓 𝑎 𝑝 𝑝 (•,ℎ 𝑛 ) ℎ 𝑛
𝑦 * 𝑛 , then there exist C, L > 0 such that:

𝑀𝑎𝑥 0⩽𝑛⩽𝑁 Φ 𝑓 𝑎 𝑝 𝑝 (•,ℎ 𝑛-1 ) ℎ 𝑛-1 • . . . • Φ 𝑓 𝑎 𝑝 𝑝 (•,ℎ 0 ) ℎ 0 (𝑦 0 ) -𝜑 𝑓 ℎ 𝑛-1 • . . . • 𝜑 𝑓 ℎ 0 (𝑦 0 ) ⩽ C𝛿ℎ 𝑝 L 𝑒 L𝑇 -1 , where ℎ = Max 0⩽ 𝑗 ⩽𝑁-1 ℎ 𝑗 .
We now focus on numerical schemes belonging to the class of explicit Runge-Kutta methods, as this allows to specify some of the constants of previous Theorem. Following Remark (i), we shall assume that 𝑓 𝑎 𝑝 𝑝 (•, ℎ) is well-defined and smooth on the compact set Ω × [𝑡 -, 𝑡 + ], so that it is Lipschitz with Lipschitz constant

𝜆 := 𝑀𝑎𝑥 ℎ∈𝐻 𝑑𝑓 𝑎 𝑝 𝑝 (•, ℎ) 𝐿 ∞ (Ω) . (11)
Corollary 1. Suppose that the numerical scheme Φ ℎ from Theorem 1 is the Runge-Kutta method with Butcher tableau ( 𝐴, 𝑏) where 𝐴 = (𝑎 𝑖, 𝑗 ) 1⩽ 𝑗 ⩽𝑖⩽𝑠 ∈ M 𝑠 (R) and 𝑏 = (𝑏 𝑗 ) 1⩽ 𝑗 ⩽𝑠 ∈ R 𝑠 . Assume further that 𝐴 is strictly lower triangular, so that the scheme is explicit, and that it is of order 𝑝.

Then inequality [START_REF] De Ryck | On the approximation of functions by tanh neural networks[END_REF] of Theorem 1 holds with C = 𝛼 and L = 𝛼𝜆 where

𝛼 = ||𝑏|| 1 1 + 𝜆ℎ + || 𝐴|| ∞ 𝑒 𝜆ℎ + || 𝐴|| ∞ .
Remarks.

1.

As the approximation 𝑓 𝑎 𝑝 𝑝 (•, ℎ) of fℎ contains an O (ℎ 𝑝 )-error term, the order of the new numerical procedure coincides with the order of the underlying scheme Φ ℎ . However, as soon as the NN becomes large, 𝛿 is small enough for the combined procedure to be significantly more accurate then the simple application of Φ ℎ .

2.

For 𝑁 𝑡 = 1, with 1 hidden layer, we have density of MLP's in C 1 (Ω) for the Sobolev norm

𝑊 1,∞ [18]. If we have Max ℎ∈𝐻 | fℎ (•)-𝑓 𝑎 𝑝 𝑝 (•,ℎ) | ℎ 𝑝 𝑊 1,∞ (Ω)
⩽ 𝛿, then we get

𝜆 ⩽ 𝑀𝑎𝑥 ℎ∈𝐻 𝑑 fℎ 𝐿 ∞ (Ω) + 𝛿ℎ 𝑝 + . (12

An alternative method for parallel training

In this subsection, we show how to learn the modified field in an alternative way. The main idea consists in training separately each term (say for instance of the modified field for the forward Euler method), by creating different data sets for different time steps

ℎ 1 < • • • < ℎ 𝑁 ℎ : 𝑦 𝑗 = 𝑦 0 + ℎ 𝑗 𝑓 (𝑦 0 ) + ℎ 2 𝑗 𝑓 [1] (𝑦 0 ) + • • • + ℎ 𝑁 𝑡 𝑗 𝑓 [𝑁 𝑡 -1] (𝑦 0 ) + ℎ 𝑁 𝑡 +1 𝑗 𝑅(𝑦 0 , ℎ 𝑗 ). (13)
Note that in contrast with previous method, the step-size is not chosen at random for each initial value. We then obtain a linear system which can be solved by using the generalized inverse of a matrix. The solution of this linear system encompasses the values of 𝑓 [1] (𝑦 0 ), ..., 𝑓 [𝑁 𝑡 -1] (𝑦 0 ) and 𝑅(𝑦 0 , ℎ 1 ), ..., 𝑅(𝑦 0 , ℎ 𝑁 ℎ ) which correspond to data usable for learning each term of the modified field separately.

The main advantages of this method are a shorter training-time (at least on a parallel machine), thus allowing for a larger number of data, and a smaller computational time (again on a parallel machine) when it comes to obtaining the numerical solution from trained values.

Numerical experiments

In order to illustrate our theoretical results, we have tested the method given in Section 2.1 for two simple dynamical systems used from simple physics:

1. The Non-linear Pendulum: This system describes the movement of a pendulum under the influence of gravity. It is governed by the equations

𝑦 1 = -sin(𝑦 2 ) 𝑦 2 = 𝑦 1 ,
where 𝑦 2 denotes the angle of the pendulum with respect to the vertical and 𝑦 1 its angular velocity. Note that the system is Hamiltonian, see [START_REF] David | Symplectic learning for Hamiltonian neural networks[END_REF][START_REF] Hairer | Geometric Numerical integration: structure-preserving algorithms for ordinary differential equations[END_REF][START_REF] Jin | SympNets: Intrinsic structurepreserving symplectic networks for identifying Hamiltonian systems[END_REF]. Parameters are given in Appendix B.2.1 for the forward Euler method, Appendix B.2.3 for the Runge-Kutta 2 method and Appendix B.2.4 for the midpoint rule.

2. The Rigid Body system: This is a three-dimensional system which describes the angular rotation of a solid in the physical space

             𝑦 1 = 1 𝐼 3 -1 𝐼 2 𝑦 2 𝑦 3 𝑦 2 = 1 𝐼 1 -1 𝐼 3 𝑦 1 𝑦 3 𝑦 3 = 1 𝐼 2 -1 𝐼 1 𝑦 1 𝑦 2 ,
where 𝑦 1 , 𝑦 2 and 𝑦 3 denote the angular momenta, and 𝐼 1 , 𝐼 2 and 𝐼 3 the momenta of inertia [START_REF] Hairer | Geometric Numerical integration: structure-preserving algorithms for ordinary differential equations[END_REF] (we take here 𝐼 1 = 1, 𝐼 2 = 2, 𝐼 3 = 3). It possesses two invariants, the so-called Casimir 𝐶 (𝑦) = For the Forward Euler and RK2 methods, the modified field (2) can be computed by recursive formulae based on various derivatives of 𝑓 , see for instance [START_REF] Chartier | Numerical integrators based on modified differential equations[END_REF][START_REF] Hairer | Geometric Numerical integration: structure-preserving algorithms for ordinary differential equations[END_REF]. It is represented by a series whose general term 𝑓 [ 𝑗] has an explicit -though complicated-expression, which can be compared with its numerical counterpart, obtained by learning it from the data set. For all 𝑦 ∈ R 𝑑 , 1 ⩽ 𝑗 ⩽ 𝑘 -1, we have on the one hand 𝑓 [1] (𝑦) =

1 2 𝑑𝑓 (𝑦) 𝑓 (𝑦) (14
)

𝑓 [ 𝑗] (𝑦) = 1 𝑗 + 1 𝑑𝑓 [ 𝑗-1] (𝑦) 𝑓 (𝑦) (15)
for the Forward Euler method, and on the other hand 𝑓 [1] (𝑦) =

1 24 𝑑 (𝑑𝑓 • 𝑓 )(𝑦) 𝑓 (𝑦) + 1 8 𝑑𝑓 (𝑦) 2 𝑓 (𝑦) (16)
𝑓 [2] (𝑦) =

1 24 𝑑 (𝑑 (𝑑𝑓 • 𝑓 ) • 𝑓 ) 𝑓 (𝑦) - 1 2
𝑑𝑓 (𝑦) 𝑓 [1] (𝑦) -

1 2
𝑑𝑓 [1] (𝑦) 𝑓 (𝑦).

for the Runge-Kutta 2 method. Formulas associated to the midpoint method are given in [START_REF] Chartier | Numerical integrators based on modified differential equations[END_REF].

Truncating the formal power series (2) then gives an approximation of the theoretical modified field, which serves as a reference

fℎ (𝑦) = 𝑓 (𝑦) + ℎ 𝑝 𝑘-1 ∑ 𝑗=1 ℎ 𝑗-1 𝑓 [ 𝑗] (𝑦) + O ℎ 𝑘+𝑝-1 . (18)
Here, 𝑝 is the order of the numerical method under consideration.

Approximation of the modified field

In this subsection, we study the approximation error between the learned modified field (3) and the theoretical modified field (2) for the nonlinear Pendulum. We observe the learning error w.r.t. both space and time step variables. More precisely, we plot the function

𝑔 𝑘 ℎ : 𝑥 ↦ → 1 ℎ 𝑝 fℎ 𝑘 (𝑥) -𝑓 𝑎 𝑝 𝑝 (𝑥, ℎ) (19)
for 𝑘 = 4 over the domain Ω = [-2, 2] 2 in order to study the learning error in space, where the O-term in fℎ (𝑦) is simply neglected. We furthermore represent Max Ω 𝑔 𝑘 ℎ for several values of time steps ℎ, in order to study the learning error in function of the the time step. Note that we clearly get the expected order of convergence of 𝑓 𝑎 𝑝 𝑝 (•, ℎ) towards the modified field fℎ , with the exception of a plateau for small values of ℎ . ℎ is globally constant at the center of the domain and grows near its boundaries (see [START_REF] David | Symplectic learning for Hamiltonian neural networks[END_REF] where a similar behaviour is observed).

Altogether, these experiments confirm that the modified field can be appropriately learned with our neural network. Error between learned field and modified field with RK2 𝑘 and 𝑓 𝑎 𝑝 𝑝 (•, ℎ) for 1 ⩽ 𝑘 ⩽ 4. Note that owing to the structure of the modified field for midpoint method (see [START_REF] Chartier | Numerical integrators based on modified differential equations[END_REF][START_REF] Hairer | Geometric Numerical integration: structure-preserving algorithms for ordinary differential equations[END_REF]), terms for odd powers of ℎ vanish.

Error between learned field and modified field with Forward Euler h h h h 2 h h

Loss decay and Integration of ODE's

Now, in order to compare, for a given method, the integration of a dynamical system with the original field and with the learned modified field, we will solve the nonlinear Pendulum with the Forward Euler, Runge-Kutta 2 and midpoint methods and the Rigid Body system with the Forward Euler method. However, prior to that, we study the decays of the 𝐿𝑜𝑠𝑠-functions for the training and testing data sets (𝐿𝑜𝑠𝑠 𝑇𝑟𝑎𝑖𝑛 and 𝐿𝑜𝑠𝑠 𝑇 𝑒𝑠𝑡 ). Their similarity is a good indication that there is no overfitting (the size of the training data set is thus appropriately estimated). As the MSE 𝐿𝑜𝑠𝑠 is used, it gives an idea of the value of the square of the learning error. Figures 4, 5 and 6 show a more accurate numerical integration by using the corresponding learned modified field 𝑓 𝑎 𝑝 𝑝 (•, ℎ) than using 𝑓 . Moreover, exact flow and numerical flow with 𝑓 𝑎 𝑝 𝑝 (•, ℎ) seem identical due to the small numerical error.

Local error

Comparison of local errors

|( f h ) n (y0) f nh (y0)| |( fapp h ) n (y0) f nh (y0)|

Local error

Comparison of local errors As the midpoint method is a symmetric and symplectic method, it is known to preserve accurately the geometric properties of the model. In order to evaluate the extent to which this feature persists in our context, we simply plot the value of the Hamiltonian along the numerical solution obtained from learned data. We test this method for the pendulum system, which is hamiltonian. Figure 7 shows a smaller error for integration with 𝑓 𝑎 𝑝 𝑝 (•, ℎ) by using the midpoint method than integration with 𝑓 . Moreover, the hamiltonian function of the pendulum system, given by

|( f h ) n (y0) f nh (y0)| |( fapp h ) n (y0) f nh (y0)|
|( f h ) n (y0) f nh (y0)| |( fapp h ) n (y0) f nh (y0)|
𝐻 : 𝑦 ↦ → 1 2 𝑦 2 1 + (1 -cos(𝑦 2 )) (20)
is preserved by the midpoint method with 𝑓 𝑎 𝑝 𝑝 (•, ℎ) with smaller oscillations than midpoint with 𝑓 . Preservation is better than DOPRI5 too, which is a non-symplectic method, as shown in Figure 8.

Local error

Comparison of local errors Eventually, we study the global error between the exact flow and the approximation obtained from the original field, as well as the error between the exact flow and the numerical flow obtained form the learned modified field. The errors are plot as functions of the step-size and the curves are in perfect agreement with the estimates of previous theorems (for the Forward Euler, Runge-Kutta 2 and midpoint methods).

|( f h ) n (y0) f nh (y0)| |( fapp h ) n (y0) f nh (y0)|

Global error

Error between trajectories with Forward Euler Forward Euler -f Forward Euler -fapp Figure 10: Integration errors (green: integration with 𝑓 , red: integration with 𝑓 𝑎 𝑝 𝑝 (•, ℎ)). Left: midpoint method for the nonlinear Pendulum. Right: Forward Euler method for the Rigid Body System.

Note that the estimate of Theorem 1 is confirmed by Figures 9 and10 , with a smaller multiplicative constant for the RK2 method.

Computational times for explicit methods

In this subsection, we plot efficiency curves (global error w.r.t. computational time). As the main goal of this paper is to design cheaper and/or more accurate solvers, we shall compare our results with the state-of-the-art Dormand & Prince methods [START_REF] Hairer | Solving ordinary differential equations[END_REF]. Figures 11 and 12 show numerical errors for explicit methods (Forward Euler and Runge-Kutta 2) with 𝑓 𝑎 𝑝 𝑝 (•, ℎ) can be smaller than numerical errors for DOPRI5 with 𝑓 , especially for large time steps. Comparison between computational time and integration error (red: numerical method with 𝑓 , green: integration with 𝑓 𝑎 𝑝 𝑝 (•, ℎ), yellow: integration with DOPRI5) for Rigid Body system with Forward Euler method.

Besides, we compare the integration error with the learned modified field and the integration error using truncated modified field at the order 𝑘 fℎ 𝑘 . For the forward Euler method, use this numerical method with fℎ 𝑘 will give a numerical method of order 𝑘, called modified Euler [START_REF] David | Symplectic learning for Hamiltonian neural networks[END_REF][START_REF] Hairer | Geometric Numerical integration: structure-preserving algorithms for ordinary differential equations[END_REF]. Figure 13 shows for both systems that integration with fℎ 𝑘 is more expensive in time or less accurate that integration with 𝑓 𝑎 𝑝 𝑝 (•, ℎ). Moreover, we always conserve an interesting integration error whereas the other methods are not accurate at all for small computational time.

Evaluation of alternative method

In this subsection, a comparison between the two learning techniques is illustrated. For the comparison to be fair, we adapt the volume of data so that the training times are nearly identical.

Figure 14 shows that both methods are able to generalize for smaller time steps than those used for training. In the two cases, traditional and alternative method with parallel training give approximately identical results

Conclusions

The numerical experiments presented in this paper demonstrate that learning the modified vector field is very beneficial in terms of efficiency. Clearly, the training of the network is an overload which should be considered separately for real-time applications. Another interesting outcome of our study is the fact that even when explicit formulae of the various terms of the modified field are available, it is advantageous to use its learned counterpart. Eventually, even though the technique used here is not fitted to situations where invariants should be conserved (learning directly the modified Hamiltonian would definitely be a better option), the learned vector field approximately retains the properties of its exact counterpart. It remains to be emphasized that the problems that we solved here are of small dimension and further studies with larges systems (for instance originating from PDEs) are needed.

A Proof of Theorem 1

Let us consider, for time step ℎ and for all 𝑛 ∈ [[0, 𝑁]] (where ℎ = 𝑇/𝑁), the numerical scheme

𝑦 * 0 = 𝑦 0 , 𝑦 * 𝑛+1 = Φ 𝑓 𝑎 𝑝 𝑝 (•,ℎ) ℎ 𝑦 *

𝑛

The consistency error is of the form 𝐿 𝑓 𝑎 𝑝 𝑝 (•,ℎ) . A discrete Grönwall lemma then leads to

𝑒 * 𝑛 ⩽ 𝐶𝛿ℎ 𝑝+1 𝑛-1 ∑ 𝑗=0
𝑒 𝐿(𝑛-𝑗-1)ℎ ⩽ 𝐶𝛿ℎ 𝑝+1 𝑒 𝐿𝑛ℎ -1 𝑒 𝐿ℎ -1

⩽ 𝐶𝛿ℎ 𝑝 𝐿 𝑒 𝐿𝑇 -1 .

B Choice of the parameters B.1 Link between learning error and parameters

The influence of the number of parameters over the learning error has been studied. In particular, we have studied the effect of the number of neurons and hidden layers. The dynamical system chosen for the test is the non-linear pendulum while the numerical method is simply the Forward Euler method. In order to approximate the learning error, we compute the value

𝛿 ≈ 𝑀𝑎𝑥 ℎ∈𝐻 *
where (𝑥 𝑖, 𝑗 ) 0⩽𝑖, 𝑗 ⩽40 is a uniform grid on the square Ω = [-2, 2] 2 , 𝐻 * = (𝑒 ℎ * 𝑗 ) 0⩽ 𝑗 ⩽14 where (ℎ * 𝑗 ) 0⩽ 𝑗 ⩽14 is a uniform discretization of [log(ℎ -), log(ℎ + )] and fℎ 4 corresponds to the field (2) for 𝑘 = 4 with 𝑅 neglected, computed via the formulas given in [START_REF] Chartier | Numerical integrators based on modified differential equations[END_REF][START_REF] Hairer | Geometric Numerical integration: structure-preserving algorithms for ordinary differential equations[END_REF].

As observed in Figure 15, a plateau appears when the number of parameters is large. This plateau has a lower value for a larger number of data. Moreover, we observe that deep networks are more efficient than shallow networks in order to learn the good vector field.

B.2 Parameters selected in our numerical experiments

For the training, the optimizer Adam of Pytorch is used. Besides, the mini-batching option is activated as it appears to be more efficient. Hyperbolic tangent tanh was selected as activation function.

B.2.1 Nonlinear Pendulum -Forward Euler

Figures 1 ,

 1 Figures 1,2 and 3 show that the error 𝑔 4ℎ is globally constant at the center of the domain and grows near its boundaries (see[START_REF] David | Symplectic learning for Hamiltonian neural networks[END_REF] where a similar behaviour is observed).

Figure 1 :

 1 Figure 1: Forward Euler method. Left: Difference between fℎ 4 and 𝑓 𝑎 𝑝 𝑝 (•, ℎ) for ℎ = 0.1. Right:Error between fℎ 𝑘 and 𝑓 𝑎 𝑝 𝑝 (•, ℎ) for 1 ⩽ 𝑘 ⩽ 4.

Figure 2 :

 2 Figure 2: Runge-Kutta 2 method. Left: Difference between fℎ 4 and 𝑓 𝑎 𝑝 𝑝 (•, ℎ) for ℎ = 0.1. Right: Error between fℎ 𝑘 and 𝑓 𝑎 𝑝 𝑝 (•, ℎ) for 1 ⩽ 𝑘 ⩽ 4.

Figure 3 :

 3 Figure 3: Midpoint method. Left: Difference between fℎ 4 and 𝑓 𝑎 𝑝 𝑝 (•, ℎ) for ℎ = 0.05. Right: Error between fℎ𝑘 and 𝑓 𝑎 𝑝 𝑝 (•, ℎ) for 1 ⩽ 𝑘 ⩽ 4. Note that owing to the structure of the modified field for midpoint method (see[START_REF] Chartier | Numerical integrators based on modified differential equations[END_REF][START_REF] Hairer | Geometric Numerical integration: structure-preserving algorithms for ordinary differential equations[END_REF]), terms for odd powers of ℎ vanish.

Figure 4 :

 4 Figure 4: Comparison between 𝐿𝑜𝑠𝑠 decays (green: 𝐿𝑜𝑠𝑠 𝑇𝑟𝑎𝑖𝑛 , red: 𝐿𝑜𝑠𝑠 𝑇 𝑒𝑠𝑡 ), trajectories (dashed dark: exact flow, red: numerical flow with 𝑓 , green: numerical flow with 𝑓 𝑎 𝑝 𝑝 (•, ℎ)) and local error (blue: exact flow and numerical flow with 𝑓 , yellow: exact and numerical flow with 𝑓 𝑎 𝑝 𝑝 (•, ℎ) ) for the nonlinear pendulum with Forward Euler method.

Figure 5 :

 5 Figure 5: Comparison between 𝐿𝑜𝑠𝑠 decays (green: 𝐿𝑜𝑠𝑠 𝑇𝑟𝑎𝑖𝑛 , red: 𝐿𝑜𝑠𝑠 𝑇 𝑒𝑠𝑡 ), trajectories (dashed dark: exact flow, red: numerical flow with 𝑓 , green: numerical flow with 𝑓 𝑎 𝑝 𝑝 (•, ℎ)) and local error (blue: exact flow and numerical flow with 𝑓 , yellow: exact and numerical flow with 𝑓 𝑎 𝑝 𝑝 (•, ℎ) ) for the nonlinear pendulum with Runge-Kutta 2 method.

Figure 6 :

 6 Figure 6: Comparison between 𝐿𝑜𝑠𝑠 decays (green: 𝐿𝑜𝑠𝑠 𝑇𝑟𝑎𝑖𝑛 , red: 𝐿𝑜𝑠𝑠 𝑇 𝑒𝑠𝑡 ), trajectories (dashed dark: exact flow, red: numerical flow with 𝑓 , green: numerical flow with 𝑓 𝑎 𝑝 𝑝 (•, ℎ)) and local error (blue: exact flow and numerical flow with 𝑓 , yellow: exact and numerical flow with 𝑓 𝑎 𝑝 𝑝 (•, ℎ) ) for the Rigid Body system with Forward Euler method.

Figure 7 :Figure 8 :

 78 Figure 7: Comparison between 𝐿𝑜𝑠𝑠 decays (green: 𝐿𝑜𝑠𝑠 𝑇𝑟𝑎𝑖𝑛 , red: 𝐿𝑜𝑠𝑠 𝑇 𝑒𝑠𝑡 ), trajectories (dashed dark: exact flow, red: numerical flow with 𝑓 , green: numerical flow with 𝑓 𝑎 𝑝 𝑝 (•, ℎ)) and local error (blue: exact flow and numerical flow with 𝑓 , yellow: exact and numerical flow with 𝑓 𝑎 𝑝 𝑝 (•, ℎ) ) for the nonlinear pendulum with midpoint method.

Figure 9 :

 9 Figure 9: Integration errors (green: integration with 𝑓 , red: integration with 𝑓 𝑎 𝑝 𝑝 (•, ℎ)). Left: Nonlinear Pendulum with Forward Euler. Right: Nonlinear Pendulum with Runge-Kutta 2.

Figure 11 :

 11 Figure 11: Comparison between computational time and integration error (red: numerical method with 𝑓 , green: integration with 𝑓 𝑎 𝑝 𝑝 (•, ℎ), yellow: integration with DOPRI5). Left: Nonlinear Pendulum with Forward Euler. Right: Nonlinear Pendulum with Runge-Kutta 2.

Figure 12 :

 12 Figure 12: Comparison between computational time and integration error (red: numerical method with 𝑓 , green: integration with 𝑓 𝑎 𝑝 𝑝 (•, ℎ), yellow: integration with DOPRI5) for Rigid Body system with Forward Euler method.

Figure 14 :

 14 Figure 14: Comparison between forward Euler, traditional and alternative method for simple pendulum. Left: Comparison of local errors. Right: Comparison of integration errors.

  𝜀 *𝑛 := 𝑦(𝑡 𝑛+1 ) -Φ𝑓 𝑎 𝑝 𝑝 (•,ℎ) ℎ (𝑦(𝑡 𝑛 )) = 𝑦(𝑡 𝑛+1 ) -Φ fℎ ℎ (𝑦(𝑡 𝑛 )) =0 (Modified field) + Φ fℎ ℎ (𝑦(𝑡 𝑛 )) -Φ 𝑓 𝑎 𝑝 𝑝 (•,ℎ) ℎ (𝑦(𝑡 𝑛 )) = Φ fℎ ℎ (𝑦(𝑡 𝑛 )) -Φ 𝑓 𝑎 𝑝 𝑝 (•,ℎ) ℎ (𝑦(𝑡 𝑛 ))so that, taking[START_REF] Champagne | Les réseaux de neurones multi-couches, le comment et le pourquoi[END_REF] and (8) into account, we have𝜀 * 𝑛 ⩽ 𝐶𝛿ℎ 𝑝+1 .Upon using (6), (21) and (21), the local truncation error can then be written as𝑒 * 𝑛+1 = Φ 𝑓 𝑎 𝑝 𝑝 (•,ℎ) ℎ (𝑦 * 𝑛 ) -Φ 𝑓 𝑎 𝑝 𝑝 (•,ℎ) ℎ (𝑦(𝑡 𝑛 )) -𝜀 * 𝑛 ,and from (9), we get𝑒 * 𝑛+1 ⩽ 1 + 𝐿ℎ 𝑒 * 𝑛 + 𝜀 * 𝑛 ⩽ 1 + 𝐿ℎ 𝑒 * 𝑛 + 𝐶𝛿ℎ 𝑝+1where 𝐿 := 𝑀𝑎𝑥 ℎ∈[ℎ -,ℎ + ]

1 ℎ fℎ 4 (

 14 𝑥 𝑖, 𝑗 ) -𝑓 𝑎 𝑝 𝑝 (𝑥 𝑖, 𝑗 , ℎ) 𝑙 ∞

  weights (100000 data) w w 1/2 learning error (HL = 1) learning error (1 HL 10), = 10 learning error (1 HL 10), = 20 learning error (1 HL 10), = 50 learning error (1 HL 10), = 80 learning error (1 HL 10), = 100 learning error (1 HL 10), = 150 learning error (1 HL 10)weights (500000 data) w w 1/2 learning error (HL = 1) learning error (1 HL 10), = 10 learning error (1 HL 10), = 20 learning error (1 HL 10), = 50 learning error (1 HL 10), = 80 learning error (1 HL 10), = 100 learning error (1 HL 10), = 150 learning error (1 HL 10), = 200

Figure 15 :

 15 Figure 15: Learning error 𝛿 versus number of parameters 𝑤 in the neural network (number of neurons 𝜁 and hidden layers 𝐻 𝐿). Shallow network is plot with red points whereas deep network is plot with green points (light green for 10 neurons to dark green for 200 neurons). Bias are neglected. Learning error is plot for 100 000 data (left) and 500 000 data (right). 200 epochs are used. [ℎ -, ℎ + ] = [10 -2 , 10 -1 ]. The curve of 𝑤 ↦ → 𝑤 1 2 has been added for comparison purposes.

  The domain Ω used for training is thus chosen accordingly. Parameters are given in Appendix B.2.2.

	intersection of the sphere |𝑦| 2 = |𝑦(0)| 2 and of the ellipsoïd	𝑦 2 1 𝐼 1 +	𝑦 2 2 𝐼 2 +	𝑦 2 3 𝐼 3 = 2𝐻 (𝑦(0)).
	2	𝑦 2 1 𝐼 1 +	𝑦 2 2 𝐼 2 +	𝑦 2 3 𝐼 3 . Hence, the solution lies at the

1

2 |𝑦| 2 and the energy 𝐻 (𝑦) = 1

  Figure 13: Computational time versus integration error for the Forward Euler method (red: with

		10 1 10 0 Computation time vs Global error -AI method vs modified field Forward Euler -fapp Forward Euler 0rder 2 Forward Euler 0rder 3 Forward Euler 0rder 4	10 1 Computation time vs Global error -AI method vs modified field 10 0 Forward Euler -fapp Forward Euler 0rder 1 Forward Euler 0rder 2 Forward Euler 0rder 3 Forward Euler 0rder 4
	Global error	10 3 10 2				Global error	10 3 10 2
		10 4
		10 5					10 4
		10 2	10 1 Computation time (s)	10 0	10 3	10 2 Computation time (s) 10 1	10 0
	𝑓 , yellow: with fℎ	2 , magenta: with fℎ	3 , cyan: with fℎ

4

, green: with 𝑓 𝑎 𝑝 𝑝 (•, ℎ). Left: Nonlinear Pendulum. Right: Rigid Body system.

𝑇 -1 , for respectively Forward Euler and RK2 methods.

Acknoledgements

The authors of this paper would like thank Pierre Navaro for his advises in the implementation aspects.

B.3 Nonlinear Pendulum -Comparison between traditional and alternative method

For the standard method, data are created by simulating several solutions with the same initial condition for different time steps. For the alternative method with parallel training, the total computational time for training correspond to the maximum of all training times for each term of the modified field.


Sylvain Kahane
email: sylvain@kahane.fr

Ziqian Peng
email: ziqian.peng@universite-paris-saclay.fr

Kim Gerdes
email: kim.gerdes@universite-paris-saclay.fr

Word order flexibility: a typometric study

This paper introduces a typometric measure of flexibility, which quantifies the variability of head-dependent word order on the whole set of treebanks of a language or on specific constructions. The measure is based on the notion of head-initiality and we show that it can be computed for all of languages of the Universal Dependency treebank set, that it does not require ad-hoc thresholds to categorize languages or constructions, and that it can be applied with any granularity of constructions and languages. We compare our results with Bakker's (1998) categorical flexibility index. Typometric flexibility is shown to be a good measure for characterizing the language distribution with respect to word order for a given construction, and for estimating whether a construction predicts the global word order behavior of a language.

Introduction

For half a century, research in typology centers on the discussion of word order parameters, pioneered by [START_REF] Greenberg | Some universals of grammar with particular reference to the order of meaningful elements[END_REF][START_REF] Greenberg | Language Universals[END_REF], and elaborated by such authors as [START_REF] Hawkins | Word Order Universals[END_REF], [START_REF] Dryer | The Greenbergian word order correlations[END_REF] and [START_REF] Nichols | Linguistic Diversity in Space and Time[END_REF]. [START_REF] Bakker | Flexibility and consistency in word order patterns in the languages of Europe[END_REF] proposes a seminal study on word order flexibility, which we pursue here. First off, it is clear that languages differ in the flexibility of word order: Greek or Russian are more flexible than English or Chinese. Secondly, constructions differ in their flexibility across the diversity of languages: The relation between an adposition and its complement is less flexible than the relation between a verb and its direct object.

To give a first idea how this can be seen on a typometric scatter plot, consider Fig. 1 where each point corresponds to a language, with its x-value indicating the percentage of nominal dependents of adpositions on the right of the adposition and its y-value indicating the percentage of nominal object dependents of verbs to the right of the verb. As a first observation, note that a large majority of languages have x=0 or x=100, leaving a rather empty space in the middle of the scatter plot. This shows that most languages have almost only prepositions or almost only postpositions; languages that mix prepositions and postpositions equally are rare. 1 However, many more languages accept both pre-verbal and post-verbal direct objects (Guzmán [START_REF] Naranjo | Quantitative word order typology with UD[END_REF][START_REF] Gerdes | Typometrics: From implicational to quantitative universals in word order typology[END_REF]. Nonetheless, the postpositional languages on the left appear to be less strictly 1 The scatter plot view understates the fact that many languages cluster around the bottom left (0,0) and the top right (100,100) of the plot: Many strictly postpositional languages also have their nominal objects on the left, and, inversely, many strictly prepositional languages also have their nominal objects to the right of their verbal governor. This is a well-established observation in typology since [START_REF] Greenberg | Some universals of grammar with particular reference to the order of meaningful elements[END_REF], and is at the base of our choice of SUD treebanks (see Section 2 for details).

postpositional than the prepositional languages on the right are prepositional. This motivates the definition of flexibility in Section 4. See Section 3 for more details on how to compute and understand these plots.

How can we measure the flexibility of languages and constructions? What are the properties of flexibility across languages and constructions? We will try to give answers to these questions in Section 6.

Most classical approaches to typology, including [START_REF] Bakker | Flexibility and consistency in word order patterns in the languages of Europe[END_REF] and previous works, are categorical in the sense that languages are grouped into categories based on their order constraints, and often only one basic word order is assumed per language from which other word orders are derived by movement, dislocation, or similar operations.

We propose a typometric approach (also called token-based typology by Levshina 2019): With the availability of a wide range of uniformly annotated treebanks in the Universal Dependencies (UD) project, it has become possible to approach these questions empirically. Syntactic typology outgrows the need for ad hoc categories and measures of distribution of languages across empirical observations become the center of interest (Futtrell et al. 2020[START_REF] Levshina | Corpus-based typology: Applications, challenges and some solutions[END_REF]. In [START_REF] Gerdes | Typometrics: From implicational to quantitative universals in word order typology[END_REF], quantitative universals describe empty or sparsely populated spaces in unidimensional or multidimensional spaces instead of qualitative universals that are claiming rare or impossible combinations of language features based on categories. [START_REF] Tesnière | Eléments de syntaxe structurale[END_REF] proposed a classification of languages based on the dependency direction referring to [START_REF] Steinthal | Die Classification der Sprachen dargestellt als die Entwicklung der Sprachidee[END_REF] and [START_REF] Schmidt | Die Sprachfamilien und Sprachenkreise der Erde: Atlas von 14 Karten[END_REF]. He opposes strict word order, when head-daughter relations mostly go in one direction, to mitigated when the head is amidst its dependents going out in both directions. Among languages with mitigated word order, there are languages with free order, as well as languages with mixed word order, where word order is quite strict in most constructions but inconsistent between constructions. This is what flexibility measures.

In this paper we propose measures of flexibility that can be applied to dependency treebanks and discuss the distributions of these measures compared to other observations on dependency treebanks. Similar measures have been first introduced by [START_REF] Futrell | Quantifying word order freedom in dependency corpora[END_REF] under the name of word-order entropy and have been studied by [START_REF] Levshina | Token-based typology and word order entropy: A study based on Universal Dependencies[END_REF].

In this paper, we try to characterize the distribution of all languages of our sample in terms of word order direction for each construction C: We compute for each language L, the number of head-initial realization of the construction C in L, what we call the headinitiality of language L under C (Section 3). We deduce from head-initiality a second measure we call flexibility and study the relation between head-initiality and flexibility for all languages in our sample, distinguishing flexible languages from mixed word order languages (Section 4). The typometric measure of flexibility we introduce is compared with Bakker's (1998) categorical measure of flexibility, as well as a more typometric measure à la Bakker (Section 5). We show that the distribution of head-initiality for every construction C can be characterized by the average head-initiality of C and the flexibility of C (Section 6). In Section 7, we explore the question of the predictability of word order distribution from one construction to another.

2

Dependency syntax and word order Dependency syntax encodes constructions by relations between words representing combinations between larger units [START_REF] Tesnière | Eléments de syntaxe structurale[END_REF][START_REF] Hudson | Word Grammar[END_REF]. A dependency relation goes from one word to another, from governor to dependent. There is no a priori assumption on locality of a relation, and a long distance dependency, for example, does not need any special encoding in a dependency tree, which makes dependency treebanks the obvious choice when attempting to measure tendencies in word order across languages [START_REF] Liu | Dependency distance as a metric of language comprehension difficulty[END_REF].

A syntactic relation is a class of combinations of the same type, having similar properties. Dependency syntax makes the assumption that most constructions are asymmetric, with a head element controlling the distribution of the combination. In some languages, constructions are very rigid and combinations of a certain type tend to always have the same word order between the governor and the dependent. Examples of such rigid relations are the subject and the object relation in English. 2 Subject and Object are different constructions and therefore are annotated as different relations.

Although the main criteria of distinguishing one relation from another is valency and morphology (for example an accusative case can be the main criterion for delimiting the direct object relation), in some cases the definition of a dependency relation relies on the word order itself and thus the relations have by definition a strict word order. 3 Just as in the initial typometrics paper Gerdes et al. ( 2021), we rely for our measures on the Surface Syntactic Dependency (SUD) version of UD [START_REF] Gerdes | SUD or Surface-Syntactic Universal Dependencies: An annotation scheme near-isomorphic to UD[END_REF], in order to make our work comparable with previous work on word order typology, thus preserving "cross-category harmony" (Hawkings 1983) and avoiding complications in particular concerning adpositions and auxiliary verbs that are analyzed in an unusual manner in the original UD annotation scheme. 4Choosing SUD rather than UD has very little impact on the computation of the flexibility measures introduced in this paper. 5SUD's comp relation corresponds to UD's aux, ccomp, iobj, obj, obl:arg, xcomp, cop, mark, and case; mod corresponds to UD's first case. Clearly, the measures we end up with will always depend on the annotation choices of each treebank.

3

As an example, consider the annotation choices for Cantonese and Mandarin reported in [START_REF] Wong | Quantitative Comparative Syntax on the Cantonese-Mandarin Parallel Dependency Treebank[END_REF]: Any element to the left of the verb is considered as "dislocated" even if it fills the verb's object slot.

As SUD is obtained by a conversion of UD without any addition of information, the granularity remains similar, see Section 4. It only impacts locally some relations such as the subject, which, in SUD, is attached to the auxiliary rather than the content verb and whose direction can change in some cases (for instance, in German, where the subject can be between the auxiliary and the verb).

advcl, acl, advmod, amod, nmod, nummod, and obl:mod. SUD's subj combines UD's csubj and nsubj. The relations dislocated, det, and clf remain unchanged between SUD and UD.

We base our work on the latest SUD version 2.11 which includes 243 treebanks in 138 languages in total. For our study, treebanks of the same language are combined and taken as one data point. 65 UD languages cover Indo-European languages. Afro-Asiatic, Uralic, and Tupian languages have 11 languages each. Turkic covers another 6 languages. Of the remaining languages only Basque, Chinese, Classical Chinese, Indonesian, Japanese, Korean, and Naija have more than 100k tokens. 21 of the UD languages are very small (less than 1000 tokens), which falls beneath our threshold for most of our measures.

Typometrics and scatter plots

A typometric analysis does not assume a basic word order or any threshold for categorizing languages or construction. Our basic observation is the measure of head-initiality defined for a language L and a construction C involving a unique dependency as follows:6 head_initiality(L, C) = % of occurrences of C in L that are head-initial (governor < dependent)

In most cases the construction C limited to a dependency is defined as a gov-rel-dep triple (governor's POS, dependency relation, dependent's POS). In some cases the construction is defined as the sum of a series of gov-rel-dep triples. Note that any variable of the triplet (gov, rel, or dep) can be equal to all, denoting no restriction on this variable. A head-initiality of 0 for a given language and construction shows a strictly head-final construction, a head-initiality of 100 indicates a strictly head-initial construction. Measuring head-initiality across the UD languages for the combination of core dependency relations (to be defined in the next section) gives the unidimensional scatter plot of Figure 2, where, unsurprisingly, Japanese is the most head-final language, and Arabic the most head-initial language of our language sample set. 7 A second two-dimensional scatter plot (Fig. 3) opposes the head-initiality of nominal subjects in the x-axis (VERB-subj-NOUN| PROPN) and nominal direct objects in the yaxis (VERB-comp:obj-NOUN|PROPN). We allow both the UD POS noun and proper noun as arguments. We observe a typical triangular shape of the resulting distribution indicating that nearly all languages have the tendency to have direct objects more to the right than subjects. Put differently, hardly any language has a higher head-initiality on subjects than on direct objects. See [START_REF] Gerdes | Typometrics: From implicational to quantitative universals in word order typology[END_REF] for a discussion of how this observation generalizes to the well-known absence of OVS languages.

7

Colors and shapes of the language points follow the original typometrics paper with colored triangles for the different subgroups of Indo-European languages, plus signs for agglutinating languages, orange x signs for Afroasiatic and Semitic languages, and circles and stars for other groups. Data, scatterplots, and detailed captions are on https://typometrics.elizia.net/. Note that only some languages are labeled. This has no semantics and is done automatically to increase readability. The value of flexibility(L,C) ranges from 0 to 100 and measures the distance of C from a strictly head-initial or head-final construction. A very similar measure, word order entropy, has been proposed by [START_REF] Levshina | Token-based typology and word order entropy: A study based on Universal Dependencies[END_REF], inspired by [START_REF] Futrell | Quantifying word order freedom in dependency corpora[END_REF].8 She also considers the entropy for couples of dependencies, such as the relative position of subjects and objects.

Flexibility of languages

For a given language L, we can compute the weighted average of flexibility(L,C) for a relevant set S of constructions C, which will be discussed below. A measure very similar to flexibility(L) has been introduced by [START_REF] Futrell | Quantifying word order freedom in dependency corpora[END_REF], using conditional entropy. In information theory, the conditional entropy H(Y|X) quantifies the amount of information needed to describe the outcome of the random variable Y given that the value of the random variable X is known. The more H(Y|X) is close to 1, the more Y is independent from X, H(Y|X) being equal to 0. In [START_REF] Futrell | Quantifying word order freedom in dependency corpora[END_REF], X is used to select a set S of constructions, while Y describes the word order on S. In other words, entropy, like flexibility, measures the extent to which word order choices depend on syntactic constructions. 9 Let us discuss our choices of S for headinitiality and flexibility. The computation of head_initiality(L) and flexibility(L) is sensitive to the range D of data considered. Unlike headinitiality(L), the computation of flexibility(L) is sensitive to the granularity of the partition of D into a set S of constructions: the finer the partitioning S, the higher the yield of flexibility(L). In our case, we have adopted a rather fine granularity, as we consider any govrel-dep triplet as a different construction, where gov is the POS of the governor, dep is the POS of the dependent, and rel is the relationship between them. We could have used an even finer granularity, by taking into account certain features, for example by distinguishing relative pronouns (PronType= Rel) from personal pronouns (PronType=Prs) or by isolating demonstratives (PronType= Dem). 10 Moreover, when we have a direct complement of the verb, we will distinguish nominal complement (dep=NOUN) and pronominal complement (dep=PRON), but not when it is a prepositional complement (dep= ADP). Sometimes the granularity can be excessive, as when UD/SUD distinguishes proper nouns (PROPN) and common nouns (NOUN).11 It must also be remarked that Levshina (2018) restricts her computation for verbal constructions to verbs that are roots, arguing that word order can be quite different between main and subordinate clauses in some languages (German and Wolof for instance).

9

The entropy view of flexibility is very elegant, but, as mentioned by [START_REF] Levshina | Token-based typology and word order entropy: A study based on Universal Dependencies[END_REF], [START_REF] Futrell | Quantifying word order freedom in dependency corpora[END_REF] gives "one aggregate score" for each language, rather than considering individual constructions before aggregating them. 10 Levshina (2019) also considers constructions restricted to one dependent word form. This is only possible if the corpus contains enough occurrences of the word, which commonly implies for many languages to parse raw corpora that are bigger than the manually annotated corpora of UD.

Our preference is to keep all occurrences, but it is certainly interesting to do a partition between main clauses and subordinate clauses. 12 Unlike flexibility(L), the computation of head_initiality(L) is obviously very dependent on the choice made for the head of each construction. It is this question that motivated us not to work with UD, but to choose the SUD variant where adpositions, subordinating conjunctions, and auxiliaries are chosen as heads. 13 For auxiliaries, the question is delicate, because while for Indo-European languages, they are clearly heads, this is less obvious in languages where they are particles. On the other hand, the wh-words of Indo-European languages are treated as pronouns in both UD and SUD, even though they also have a head role, which explains in part their peculiar placement.

For head-initiality(L), we decide to consider the relations of type comp, mod, udep, subj, dislocated, that we call the core relations. We have included the dislocated relation, because the boundary between governed and dislocated elements is not always well defined. 14 We have eliminated the det relation because the direction of the determiner-noun relation is controversial (see the discussion around the DP-hypothesis since [START_REF] Hudson | Word Grammar[END_REF][START_REF] Abney | The English noun phrase in its sentential aspect[END_REF], as well as clf (for classifiers), which is used inconsistently. For flexibility(L), we could keep det and clf, because the choice of the governor of a given relation does not play any role:

Flexibility only measures the proportion of dependencies going in the same direction and remains the same when all dependencies are inverted. However the det and clf have only a small influence on the final result (cf. Table D) and, to be precise, we decided to use the core relations for the computation of flexibility(L) as well. Other SUD/UD relations are of little interest for our study as their direction is fixed in the UD annotation guidelines. This includes conj, fixed, goeswith, etc. It should also be noted that we have not considered the relations 12 About the relations between constructions in main and subordinate clauses, see [START_REF] Schachter | Focus and relativization[END_REF]. 13 For instance, [START_REF] Futrell | Quantifying word order freedom in dependency corpora[END_REF] based on computation on UD indicates that French and Italian are mostly head-final, while, based on SUD, they are head-initial at 76% and 77% respectively! 14 For instance, in a pro-drop language such as Chinese, it is difficult to decide if preverbal objects are dislocated or not. See Note 3.

between co-dependents at all. Yet, some languages with a very strict head-final order, such as Japanese or Korean, can have much greater freedom in the placement of codependents, which is not taken into account in the present study.

Lastly, we chose to give each construction a weight equivalent to its frequency in the corpora, unlike [START_REF] Bakker | Flexibility and consistency in word order patterns in the languages of Europe[END_REF], who gives the same weight to each of the 10 constructions he considers (as well as [START_REF] Levshina | Token-based typology and word order entropy: A study based on Universal Dependencies[END_REF], who considers quantitative values for each construction but does an average with equal weight).

Figure 4 shows the head-initiality of SUD 2.11 languages in the x-axis and the flexibility in the y-axis. For treebanks with at least 1000 core relations, we observe that Ancient Greek, Tupinambá, Emerillon(Teko), Turkish-German (code switching corpus) and Old East Slavic are the most flexible languages (with flexibilities 59.3, 56.7, 55.2, 52.9, and 51.7 respectively), while Japanese, Hindi, Xibe, Kazakh and Telugu (flexibilities of around 0.5, 1.6, 2.1, 2.4, and 2.4 respectively) are the least flexible languages.

Languages with head-initiality equal to 0 or 100 have flexibility 0 and the closer they are to 50 the more likely they are to be flexible. But there are languages L with head_initiality(L) close to 50 and flexibility(L) close to 0, such as Bambara: these are mixed order languages. Languages with high flexibility(L), such as Ancient Greek, are free order languages. 15 5 Comparison with Bakker's flexibility [START_REF] Bakker | Flexibility and consistency in word order patterns in the languages of Europe[END_REF] proposes a computation of flexibility based on the same principles but does not take into account the greater or lesser flexibility of each construction for each language. In Bakker's computation, a language is either flexible or completely rigid. flexibility [Bakker](C,L) = 0 if the construction C is completely rigid in L, 1 if it is flexible flexibility [Bakker](L)= (equal-weighted) average over 10 constructions C of flexibility [Bakker](C,L) 15 Our measures are also dependent on the corpus chosen for the calculus and its genre. The flexibility measure of Ancient Greek is certainly increased by the fact that the corpus contains poetry and theater. Bakker gives the flexibility values for a sample of 86 European languages 47 of which are Indo-European. We propose to compare Bakker's values with a typometric index of the same type. Some of Bakker's constructions can be directly translated into SUD corpus queries, others can be approximated. For example, his "Adj/N" translates directly into the typometric measure NOUN-mod-ADJ. The Verb-Recipient relation (V/R) can only be approximated by VERB-comp:obl-ADP| NOUN (cf. Table A2 in Annex). The complete list of Bakker's constructions and their translation into typometric measures are provided in the construction flexibility Table A3 of the Annex. 1616 Bakker (1998: 393ff) introduces another measure which he calls consistency and which is very dependent on the set of considered constructions, which are still the 10 same constructions (see Section 5):

consistency[Bakker](L) = | #{ C / C is head-initial for L } -#{ C / C is head-final for L } |.
It seems to us that the consistency of a language L is well captured by our head_initiality(L), which is not dependent on the partitioning P into constructions undertaken by the linguist. Moreover, Bakker (1998: 401-2) notes that flexibility [Bakker] and consistency [Bakker] are correlated, but this is obvious as soon as there are languages whose head-initiality is close to 0 or 100. Likewise, our flexibility and head-initiality are related, since flexibility(L) ≤ 2 min ( ✕ head_initiality(L), 100 -head_initiality(L) ) , which we have visualized as a triangle in Figure 4. As Bakker's flexibility measure is categorical per construction, we have to arbitrarily set a threshold at 5%, considering that languages with less than 5% variation for a given construction C are inflexible for C.

We can compare those 3 measures across the languages that we also find in UD: 1. Bakker's flexibility, 2. our recomputation of the flexibility à la Bakker, as a non-weighted average over Bakker's 10 constructions, with the 5% threshold as indicated above, and 3. our typometric flexibility. 17 The scatterplot of Figure 5 shows the strong correlation between the Bakker-like measure of flexibility and the typometric flexibility. The Bakker-like flexibility is also significantly correlated with Bakker's flexibility (Fig. C2), while there is only a weak correlation between the typometric flexibility and Bakker's flexibility (Fig. C1). See the complete results in Tables B1 andB2 in the Annex.

Flexibility of constructions

Having compared the overall flexibility of the languages, we can now see how the languages in a given sample S are distributed for each construction C and compare the constructions. Specifically, we are interested in how the headinitiality of the languages in our sample (the 138 languages in UD 2.11) is distributed for the different constructions C. Our hypothesis is 17 Some of Bakker's constructions, such as Dem/N or Pro/N, involve features that are not present in all UD treebanks (PronType=Dem and Poss=Yes in these two cases). Our computation is restrained to languages with all the required features. For a uniform distribution, flexibility = 50, head-initiality = 50, head-initial-weight = 1, and head-final-weight = 1. When head-finalweight > 1, the distribution is drawn towards 0 and when head-final-weight < 1, it is pushed away. The reverse holds for head-initialweight. Our postulate is that the distribution of head-initiality is similar to a uniform distribution that has been distorted by stretching it from both sides. 18 Our head-initial and head-final weights give us an estimate of the strength of the force at each end.

We observe that for all the most frequent C constructions, both head_initial_weightS(C) > 1 and head_final_weightS(C) > 1 (see Table A3 in the Annex where all but one of the weights for the 10 Bakker constructions are greater than 1), which means that languages are attracted on both ends.

To give an idea of the different distributions we encounter, the three scatter plots below show three head-initiality distributions on the treebanks: 1. Num/N (NOUN-any-NUM), 19 2. aux-v (AUX-comp:aux-VERB), and 3. Adj/N (NOUN-mod-ADJ). The first (Num/N) distribution tends towards 0 and is pushed away from 100 (weights 4 and 0.6), while the two other distributions are attracted both by 0 and 100, with a bigger attraction to 0 for Adj/N (weights 4.2 and 2) and to 100 for Aux/V (weights 1.3 and 2.6). To be precise, SUD uses a special feature ExtPos, indicating the external POS of a word. Numerals, all categorized NUM in UD, are nummod or det when they are used as a quantifier (my 7 cats). In other uses, they work as a proper noun (line 7, page 7, year 2023) and receive the feature ExtPos=PROPN. It is this feature, rather than the POS, that is used in all our computations. It remains that the use of nummod is not consistent across all treebanks.

Bakker restricts his study to a sample S of 86 European languages, 16 of them having enough data in UD to be compared. 20We find that V/O is the most flexible construction, followed by V/R and Adj/N. Two constructions do not behave at all as in Bakker's sample: Aux/V appears as the most flexible construction after V/O, while Rel/N appears as extremely inflexible.

Bakker also compares the flexibility of languages with head-initial and head-final basic order. Again we can introduce typometric Bakker-like measures. We consider that a language has head-initial basic order if more than 50% of core dependencies are headinitial. Bakker (1998: 392) "observed that headmodifier orders are, on the whole, more flexible than modifier-head orders." We have completely different results with our measures (see Table A3 in the Annex): Only for adpositions, languages with head-initial core order are more flexible than languages with head-final core order.

Predictivity

With these notions in place, it is now possible to measure which construction predicts best the overall core flexibility of a given language. For this, we measure the Spearman correlation between the distribution of flexibility(L,C) for various couples of construction C (Figure 8). We are particularly interested in the correlation with the core construction. Among the 10 Bakker constructions, the best predictors of the core flexibility is the V/O construction, unlike Adj/N. Note also some notable correlations: Aux/V and V/O that have a correlation of 0.59 as well as subj and comp that have a correlation of 0.53. Bakker (1998: 392) however states that "the best predictors of overall flexibility are the flexibility of the recipient, genitive, numeral and relative clause. On the other hand, no prediction whatsoever can be drawn from the behavior of adpositions and articles." This is not backed up by our data with our weighted flexibility measure: The typometric genitive flexibility has a correlation of only 0.18 with the core flexibility, numerals have a correlation of 0.01, and relatives of 0.07. Nonetheless, our data agrees with a common typologist view, most notably [START_REF] Dryer | The Greenbergian word order correlations[END_REF], that sees V/O as a good predictor for word order constraints on other constructions.

Note also that the constructions with the biggest flexibility are the best predictors. This was expected because constructions with low flexibility tend to gather all languages around 0 and 100.

Conclusion

We have introduced a first measure, headinitiality, which measures the variable headdependent word order on a language's treebank or on specific constructions. Based on this, we develop an operational notion of flexibility that renders the intuition that the average headinitiality can be far from 0 or 100 while the languages are strict per given construction.

We then show that our empirical notion of flexibility can be compared to previous definitions of flexibility of word order, notably to Bakker's work. Our notion of flexibility has the advantage that it can directly be computed from treebanks, that it does not require ad-hoc thresholds to categorize languages or constructions, and that it can be applied with any granularity of constructions.

Finally, we show which construction predicts overall word order flexibility of a language. For this, we rely on Spearman's rank correlation coefficient, which allows us to calculate a correlation between two distributions. We show that over UD's language sample, the highest correlation is obtained for nominal objects (V/O construction).

Since the Spearman correlation is a symmetric measure, we would like to continue our study by proposing an asymmetric measure that allows us to decide if one distribution can predict another. Our hypothesis, to be confirmed, is that constructions with the most uniform distribution, thus being flexible and well-balanced, provide better predictions of the behavior of other constructions. The V/O construction, which many authors take as a basic construction (see in particular the study of [START_REF] Dryer | The Greenbergian word order correlations[END_REF]

Figure 1 :

 1 Figure 1: Two-dimensional scatter plots with ADPcomp:obj-NOUN/PROPN in the x-axis and VERBcomp:obj-NOUN/PROPN in the y-axis.

Figure 3 :

 3 Figure 3: Two-dimensional scatter plots with verbal subjects and objects.

  For a language L, flexibility measures the distance of a construction C from a rigid construction. In this paper, we only consider constructions involving a governor G and a dependent D by a particular relation. The construction has a wider or narrower range depending on whether the relation between G and D or the categories of G and D are more or less constrained. flexibility(L,C) = 2 min( head_initiality(L,C), ✕ 100-head_initiality(L,C) ) = twice the smallest distance of head_ initiality(L,C) to 0 or to 100

  head_initiality(L) = weighted average of head_initiality(L, C) on constructions C flexibility(L) = weighted average of flexibility(L, C) on constructions C.

Figure 4 :

 4 Figure 4: Head-initiality vs Flexibility (core)

Figure 5 :

 5 Figure 5: Typometric vs Bakker-like flexibility

  that this distribution is reasonably well described by the following two values: average_head_initiality S (C) = average of head_initiality(L,C) over the Ls in S. average_flexibility S (C) = average of flexibility(L,C) over the Ls in S. The less flexible C is on average, the more languages are attracted to 0 and 100. The average-head-initiality indicates whether 0 or 100 attracts more to one than the other. We propose two other values that will help us to better understand this attraction towards 0 and 100. head_initial_weightS(C) = average_head_initiality S (C) / average_flexibility S (C) head_final_weightS(C) = (100 -average_head_initiality S (C)) / average_flexibility S (C))

Figure 6 :

 6 Figure 6: Head-initiality language distribution for three constructions Num/N, Adj/N, Aux/V. Again, we can compare our flexibility results with two measures: the flexibility measure proposed by Bakker (1998) and a Bakker-like measure that we calculate from our sample. flexibility[Bakker](C)= % of languages in Bakker's sample that are flexible for C. flexibility[Bakker_like] S (C) = % of Ls in S with flexibility(L,C) > 5.

  19

  head_initial_flexibility[Bakker_like]S(C) = % of Ls in S with head_initiality(L,C)>50 that have flexibility(L,C)>5. head_final_flexibility[Bakker_like]S(C) = % of Ls in S with head_initiality(L,C)<50 that have flexibility(L,C)>5.

Figure 7 :

 7 Figure 7: Heatmap of Spearman correlation between the distributions of flexibility(L,C) for various couples of construction C. See Table D of the Annex for detailed values.

Table A1 . Various flexibility measures for languages where a treebank-based verification of Bakker's measures is available as described in footnote 20.

 A1 is thus an excellent candidate.

	Annex
	Languages	Bakker-flexibility Bakker_like_flexibility typometric_flexibility
	Armenian	40	60	33
	Belarusian	-	50	26
	Bulgarian	60	60	19
	Czech	-	50	30
	Danish	30	20	26
	English	40	40	13
	French	10	40	14
	German	40	30	33
	Greek	60	40	16
	Hindi	-	20	2
	Icelandic	40	90	35
	Italian	30	50	18
	Latin	90	80	48
	Naija	_	10	6
	Norwegian	40	20	21
	OldEastSlavic	-	80	52
	OldFrench	-	60	31
	Polish	60	50	35
	Russian	70	70	29
	Slovak	50	50	37
	Spanish	30	40	17
	Swedish	40	20	21
	Ukrainian	-	70	26
	WesternArmenian	-	50	37
		Bakker's 10 relations	Corresponding construction
		V/O	VERB-comp:obj-NOUN/PROPN
		Adj/N	NOUN-mod-ADJ
		Pro/N	NOUN-any-[Poss=Yes]
		V/R	VERB-comp:obl-ADP/NOUN
		Gen/N	NOUN-mod[gen]-ADP/NOUN
		Rel/N	NOUN-mod@relcl-VERB
		Adpos	ADP-comp-NOUN
		Num/N	NOUN-any-NUM
		Dem/N	NOUN-any[PronType=Dem]
		Aux/V	AUX-comp:aux-VERB

Table A2 : the 10 Bakker's relation and their corresponding constructions

 A2

	Measures	V/O Adj/N Pro/N V/R Gen/N Rel/N Adpos Num/N Dem/N Aux/V
	freqSample	3.3	3.6	0.7	0.6	1.5	0.4	5	0.9	0.7	2.5
	concerned_languages	113	96	51	59	68	42	103	85	64	87
	typometric_flexibility%	26.2 16.2 15.9	31	20.2	2.5	5.5	22	10.3	25.7
	Bakker-like_flexibility%	62.5 62.5 37.5 66.7	62.5	8.3	8.3	66.7	50	54.2
	Bakker-like-flexibility(S)% 48.2 35.4 29.4 62.7	45.6	7.1	12.6	51.8	29.7	46
	Bakker-like-head_initial(S)% 46.3 18.8 14.6 56.2	37	0	23.3	51.8	26.3	43.5
	Bakker-like-head_final(S)% 49.3 68.8	90	65.1	51.2	7.9	8.2	50	57.1	46.9
	head_initiality%	61	32	19	68	57	89	71	13	15	68
	head_initial_weight	2.3	2	1.2	2.2	2.8	35.8 12.9	0.6	1.4	2.6
	head_final_weight	1.5	4.2	5.1	1	2.1	4.3	5.3	4	8.3	1.3

Table A3 . Measures for the 10 constructions considered by Bakker (1998). Among them Bakker-like flexibility is normalized over the 24 languages in Table A1, Others are normalized with the amount of languages in the row 'concerned languages'Table B .

 A3B Spearman correlation (left) and Pearson correlation (right) between Bakker's flexibility, Bakker-like flexibility and typometrics flexibility

Widely discussed exceptions to the rigidity of subject and object in English include the relative pronoun (a person who I never met) and marginal cases of dislocation such as Chocolate I adore! As in other typological studies, we restrict our object measures to nominal objects, thus excluding the

Guzmán[START_REF] Naranjo | Quantitative word order typology with UD[END_REF], for example, find that UD's case relation stands out in their directional correlation measures.

Head-initiality is introduced in Gerdes et al.(2019, 2021), where it is called direction.

Precisely, entropy(L,C) = -p*log2(p) -(1-p)*log2(1-p), with p = head_initiality(L,C). This value also ranges from 0 to 100%, with value 0 for p=0 or 100% and 100% for p=50%. The only difference with our calculation is that entropy smoothes values for p around the 50% mark.

On the other hand, UD does not usually distinguish prepositional dependents of the verb, which are all rel=obl, whether they are arguments or modifiers. This distinction is made only in a few treebanks, notably the native SUD treebanks (with the comp and mod labels). SUD uses the udep relation, for underspecified obl dependencies when the distinction in argument and modifier is not encoded.

[START_REF] Levshina | Token-based typology and word order entropy: A study based on Universal Dependencies[END_REF], like us, uses the mean head initiality and the standard deviation to characterize the distribution of a head initiality for a given construction. The standard deviation is relevant for characterizing Gaussian distributions, but not for "stretched" distributions as here, particularly when elements tend to move away from the center and when these movements are asymmetric, with one end more attractive than the other.

When looking at a particular construction C, we only consider a language L if the treebanks of L have at least 50 occurrences of C. For Bakker-like measures to be calculated for L, the threshold of 20 must be reached for each of the 10 constructions considered.

Acknowledgements

We would like to thank the three anonymous reviewers of the Gurt-Syntaxfest 2023 for their careful and patient examination as well as for the many valuable comments they made.

This research was supported by the French National Research Project (ANR) Autogramm.