Title: META-LEARNING RUNGE-KUTTA

Abstract: Initial value problems, i.e. differential equations with specific, initial conditions, represent a classic problem within the field of ordinary differential equations (ODEs). While the simplest types of ODEs may have closed-form solutions, most interesting cases typically rely on iterative schemes for numerical integration such as the family of Runge-Kutta methods. They are, however, sensitive to the strategy the step size is adapted during integration, which has to be chosen by the experimenter. Here, we show how the design of a step size controller can be cast as a learning problem, allowing deep networks to learn to exploit structure in the initial value problem at hand in an automatic way. The key ingredients for the resulting Meta-Learning Runge-Kutta (MLRK) are the development of a good performance measure and the identification of suitable input features. Traditional approaches suggest the local error estimates as input to the controller. However, by studying the characteristics of the local error function we show that including the partial derivatives of the initial value problem is favorable. Our experiments demonstrate considerable benefits over traditional approaches. In particular, MLRK is able to mitigate sudden spikes in the local error function by a faster adaptation of the step size. More importantly, the additional information in the form of partial derivatives and function values leads to a substantial improvement in performance.

1. INTRODUCTION:

Differential equations in their general form cover an extremely wide variety of disciplines: While many applications are rather intuitive, as for instance simple Newtonian physics and engineering, other more exotic use cases include the governing of price evolution in economics (Black & Scholes, 1973) , the study of rates in chemical reactions (Scholz & Scholz, 2014) , and the modeling of population growths in biology (Lotka, 1925; Volterra, 1926) . In medicine, differential equations may be used to model cancer growth (Ilea et al., 2013) , diabetes and the glucose metabolism (Esna-Ashari et al., 2017) as well as for pharmaceutical drug design (Deuflhard, 2000) . Recently, differential equations have also been used as a way to design neural networks (Chen et al., 2018) . Unfortunately, finding an analytical solution in closed form is in many cases very difficult, if not impossible. Therefore, a variety of numerical integration methods have been developed to obtain accurate, but approximate solutions. Arguably, the most prominent ones are Runge-Kutta methods, a family of integration methods for initial value problems. However, setting up Runge-Kutta involves several design choices, one of which is the step size controller. Using an adaptive step size strategy instead of a constant step size can often increase efficiency by several orders of magnitude, c.f . (Söderlind, 2006) . Their performance is hampered by the fact that they only make use of hand-designed features.

Contribution. We show how to cast the design of a step size controller for Runge-Kutta as a learning problem, c.f . Fig. 1 . The key ingredients of the resulting Meta-Learning Runge-Kutta (MLRK) are the identification of a good performance measure and appropriate inputs.

Related Work. Various approaches to control the step sizes in Runge-Kutta methods have been proposed. While some rely on signal processing principles where the goal is to produce a smooth step size sequence in conjunction with acceptable local errors, others are based on the assumption that step sizes should be adapted to a prescribed function of the solution, e.g. to preserve structure in geometric integration (Söderlind, 2006) . In this work we will focus on control theoretic approaches which aim to keep the local error associated with a single step close to a tolerance parameter. x n x n+1 log r n 1 log r n

Optimizee:



Optimizer:

Figure 1 : Meta-Learning Runge-Kutta: The model represented by the blue block determines the step size adjustment log r n-1 using a LSTM and a linear layer, c.f . Eq. ( 10). The Runge-Kutta update then determines the new step size log h n = log h n-1 + log r n-1 . The new step size is used to perform the next Runge-Kutta step which computes the approximation y n+1 at time step t n+1 and an error estimate err n+1 according to (12) . Then the next input x n+1 for the model is computed.

Specifically, the design of the step size control algorithm is learned within MLRK. Indeed, casting algorithm design as a learning problem is not new. Andrychowicz et al. (2016) learned a gradientbased optimizer for nonlinear unconstrained optimization problems. Wichrowska et al. (2017) extended this work by introducing a hierarchical RNN architecture which improves the generalization ability of the optimizer. Schramowski et al. (2018) demonstrated the benefit of a learned projectionfree convex optimization algorithm, which relies on conditional gradients. Finally, Chen et al. (2017) cast the design of gradient-free black-bock optimization as a learning problem.

Furthermore, other approaches to solve differential equations using neural networks have been proposed, too. Lagaris et al. (1998) suggested to use a feed-forward neural network to approximate the solution of an ordinary or partial differential equation with initial or boundary conditions. Inspired by the Galerkin method, Sirignano & Spiliopoulos (2018) used a deep neural network to directly approximate the solution of a high-dimensional partial differential equation. E & Yu (2018) proposed the Deep Ritz method as a means to solve variational problems that arise from partial differential equations. Han et al. (2018) cast the problem of solving semilinear parabolic partial differential equations as a learning problem by using the reformulation of these differential equations as backward stochastic differential equations. In all of the above approaches, a neural network is used to help approximate the solution of the differential equation. The learning process corresponds to the numerical computation of a solution of the differential equation. In contrast, MLRK learns to improve the numerical integration process itself.

We proceed as follows. As the first step, i.e. identifying a good performance measure, we consider the general objective of step size control as well as the simplified and more practical objective that many controllers are based on. Then, we identify useful inputs by analyzing existing step size control algorithms. Furthermore, we show how local information about the ODE can be used as additional favorable inputs. Before concluding, we demonstrate empirically how our proposed controller can be used to improve step size control and investigate the benefit of different inputs.

2. INITIAL VALUE PROBLEMS AND RUNGE-KUTTA METHODS:

We give a brief introduction to initial value problems and Runge-Kutta methods, a numerical method for solving these problems. Furthermore, we will recap standard step size controllers and point out their underlying assumptions.

Basics of Initial Value Problems and Runge-Kutta Methods. An ordinary differential equation describes a system that depends on one variable, often referred to as time. An initial value problem additionally provides initial values: y = g(t, y) , y(t 0 ) = y 0 , with a function g :

[a, b] × R m → R m , m ∈ N.

As closed-form solutions can be very hard to find, numerical methods such as Runge-Kutta methods have been developed. Their mathematical underpinning is quite evolved (Butcher, 2008; Hairer et al., 2000; Hairer & Wanner, 2002) . A brief recap can be found in the appendix.

For iterative solvers, the local truncation error-the error committed in a single step using step size h-constitutes an important tool to control the step size as it allows to access the fitness of a step size.

Definition 2.1. The local truncation error in step n + 1 is defined as

EQUATION

where y n+1 = Φ(t n , y n , h) is a single Runge-Kutta step of a p-th order method using step size h. We call C(t n ) the principal error term.

We will also denote the local truncation error as local error and refer to the appendix for simple ways to estimate it. We denote the error estimates as err n .

The Objective of Step Size Control. As numerical integration is used to approximate a solution, high accuracy is a desired property. However, in practice, computational resources might be limited and efficient use of these resources can be crucial in some applications. Hence, an efficient step size controller should maximize the accuracy of the approximate solution while minimizing the computational cost. These are competing goals and require a trade-off, which is achieved by considering the Lagrangian of the error E(H) and the work cost W (H) produced by a sequence of step sizes Butcher (2008) uses two integrals describing the error E(H) and the work cost W (H). Under some assumptions one can show that optimal step sizes w.r.t. this objective causes local errors equal to the Lagrange multiplier λ, which is often referred to as the tolerance parameter. Namely, we need to assume that the global error is equal to the sum of local errors and all step sizes are small enough so that the loss can be approximated by the integrals suggested by Butcher. However, these are quite strong assumptions. The optimal step size depends on the tolerance parameter, and hence, might be too large. Furthermore, the assumption that the global error is equal to the sum of the local errors is not true in general.

H = (h n ) N -1 n=0 , E(H) + λ • W (H) .

Classical:

Step Size Control. Most common step size control mechanisms aim to keep the local errors equal to a tolerance parameter so that a standard controller is based on the following idea: For a Runge-Kutta method of order p, the local error for step size h at time t is approximately err ≈ C(t)h p+1 , and, accordingly, the desired optimal step size h opt is given by h opt ≈ h (tol/err) 1 p+1 . This formula is then used to adapt the step size as

EQUATION

where a safety factor 0 < γ < 1 and a minimal and maximal factor α, β ∈ R ≥0 , α ≤ β are included to avoid excessive step rejection and ensure a smooth step size control. Söderlind (2002) pointed out that the underlying assumptions of Eq. ( 2) are rather strong. It assumes slow variations in C and requires h to be sufficiently small so that it exhibits its theoretical asymptotic behavior. Both are not true in general.

Control Theory on:

Step Size Control in Runge-Kutta Methods. Gustafsson et al. (1988) discussed step size control in the context of a proportional-integral-derivative (PID) controller. The classical step size control mechanism in Eq. ( 2) can be regarded as an I-controller. They then demonstrated the oscillatory behavior of this controller when applied to certain problems. They report the poor stabilizing capability of an I-controller as the origin of these oscillations, which are further accentuated by a large integration gain. To overcome this, they Gustafsson et al. (1988) suggested a PI-controller, which in addition to the integral term used in a standard step size control, also includes a proportional term and can be expressed as

h n+1 = r n h n , r n = max α, min β, γ tol err n+1 n1 tol err n n2 .

(3)

A Predictive Controller. Both the classical step size controller, see Eq. ( 2), and the PI-controller in Eq. (3) rely on the assumption that the local truncation error can be described as a function of h that remains independent of t. This is certainly not true for all cases. Therefore, Gustafsson (1994) , e.g., discussed a controller based on a prediction of the principal error term. A simple model that assumes a constant linear trend in log C(t) was suggested. Implementing the predicted principal error leads to

h n+1 = h n h n-1 tol err n+1 n1 err n err n+1 n2 h n .

(4)

3 META-LEARNING RUNGE-KUTTA: NEURAL STEP-SIZE CONTROLLER

Step size controllers for Runge-Kutta are typically still designed by hand. Formulas such as Eq. ( 2), (3) and (4) require parameter fine-tuning. In this work, we take a different tack and instead suggest to replace these hand-designed update rules with a learned controller, hereinafter referred to as optimizer. The key steps to cast the design of a controller as a learning problem is to determine a good performance measure and appropriate inputs.

Performance measure. As the first step we need to identify an appropriate performance measure.

As discussed, a trade-off between numerical accuracy and computational cost of the solution must be made, e.g. by considering the Lagrangian of error and work cost. The work cost is generally correlated with the number of integration steps. It is crucial to take the length of the integration interval t nt 0 into account as well, which leads to the following performance measure:

EQUATION

In contrast, if we follow the concept of most step size controllers, an optimal step size is achieved if the local error is equal to the tolerance, which suggests the performance measure

EQUATION

)

Both types of loss functions offer various strengths and weaknesses over one another.

The loss Eq. ( 5) incorporates the essential quantities that an efficient algorithm should minimize. However, the true global error of numerical integration is generally unknown and can only be approximated. In general this can be computationally demanding and as a result make the training process very expensive. Furthermore, when the global error is employed, optimal step sizes may depend on the integration interval as well, e.g. the choice of step size in the first few steps may have different effects on the global error at different points in the integration interval.

Classic step size controllers are typically based on local errors alone without regards to global strategies. The objective of these controllers is reflected by the loss Eq. ( 6). Although this objective minimizes the Lagrangian of the global error and the number of steps only under strong assumptions it does qualify as an appropriate practical objective for step size controllers.The advantage here is that estimates of the local error associated with a step size are available and no additional computation to determine the loss of a step size is required. This allows efficient training even when the analytical solution of an initial value problem is unknown.

Input Features. Existing step size controllers rely on error estimates to update the step size. Specifically, tol errn is important, as evident in Eq. ( 2), ( 3) and (4). We intuitively expect a data driven approach to be able to utilize these features. This gives rise to the first set of input features,

EQUATION

Here, ψ (err n+1 , •) denotes the input for the optimizer with optional variables. Next, we want to investigate additional input features that may be beneficial. The principal error term as well as higher order terms of the local errors can be expressed with elementary differentials of the ODE, c.f . (Hairer et al., 2000, Section II.3 ). These elementary differentials are formed of function values and partial derivatives of g and allow a Taylor approximation of the local error function err

n+1 (h) = φ(t n , h)h k+1 + O(h k+2 )

, where φ(t n , h) is a polynomial in h with coefficients formed from partial derivatives of g up to order k > p. The roots of the polynomial φ(t n , h)h k+1 -tol approximate the optimal step size. Moreover, it is reasonable to assume that the local error takes values greater tol for some h > 0, implying the existence of real roots. Since the complex roots of polynomials are known to depend continuously on the coefficients of the polynomial, c.f . (Rahman & Schmeisser, 2002, Theorem 1.3.1), it is clear that the optimal step size, a real root of p(h), depends continuously on the partial derivatives and function values of g as well as the tolerance tol. From the Universal Approximation Theorem (Hornik, 1991) , it then follows that the optimal step size function can be approximated by a neural network with input (∂ α g(t n , y n )) |α|=k , . . . , ∂g(t n , y n ), g(t n , y n ), tol . This suggests that the partial derivatives of g can also be an appropriate input for our controller,

EQUATION

However, providing partial derivatives requires additional computation whereas the p function values that were computed during the Runge-Kutta step are conveniently available. As a trade-off between "perfect" information in the form of higher order partial derivatives and the computational effort to compute them, we also propose to use the following input:

EQUATION

where g

(1)

n , . . . , g

n denotes p function values of the initial value problem. In contrast to the traditional input (7), our novel inputs ( 8) and ( 9) both provide local information about the structure of the problem at hand, which allows one to keep the local error close to the tolerance parameter. Moreover, hand-designed controller and existing Runge-Kutta methods do not make use of these information, and it is difficult-if not impossible-to do so. Training an optimizer with input (8) or ( 9), as shown next, does so automatically and designs novel Runge-Kutta methods.

Meta-Learner. With a good performance measure and appropriate inputs at hand, we can now solve the step-size control problem as a meta-learning problem as sketched already in Fig. 1 . Let the input for the optimizer be x n = ψ (err n , •). We parameterize the optimizer c using φ and update the step size as follows, log h n+1 = log h n + c n (ψ(err n+1 , •), φ) . An approach in log space has the advantage that that we do not need to constrain the output of c to be positive. Since c n corresponds to log r n , we will adjust the notation in a similar fashion. Due to their natural ability to handle sequential tasks, an LSTM was chosen in our experiments as optimizer c, and its prediction determines the step size which will be used in the next Runge-Kutta iteration. We denote the model of the optimizer, represented by the blue block in Fig. 1 , by m, its parameters by φ and the hidden state by ĥn . It can be expressed by

EQUATION

The subsequent Runge-Kutta update in the yellow block of Fig. 1 first updates the step size

EQUATION

and then uses it to execute the next Runge-Kutta step,

EQUATION

Here, Φ denotes the Runge-Kutta step from y n to y n+1 as well as the error estimate err n+1 that is computed during that step. Afterwards, the next input for the model will be computed.

Learning Objective. Different initial value problems can require very different step sizes. Using an optimizer that is specialized for a certain class of problems allows to exploit the structure of these problems. The behavior of an initial value problem is described by a function g, therefore we can represent a class of initial value problems with a distribution over the functions g. Thus, an optimizer can be considered optimal for a class of problems, if it minimizes the expected loss:

EQUATION

)

where t n , y n , err n were computed according to Eqs. ( 10), ( 11), ( 12). Here, l denotes the loss function defined for an approximation y n at time point t n and the local error estimate err n . Note that our training loss directly corresponds to the performance measure we are interested in. The model m can then learn the behavior of the given class of problems and use this knowledge to generalize to new examples of the same class and new, unseen classes. Specifically, since the provided performance measure l is differentiable a.e., we can optimize the learning objective Eq. ( 13) using gradient descent on the parameters φ. An estimate of the gradient ∂L(φ) ∂φ can be computed by sampling a function g from the distribution of the class of initial value problems and applying backpropagation, c.f . (Rumelhart et al., 1986) .

4. EXPERIMENTAL EVIDENCE:

Our intention here is to investigate the benefits of Meta-Learning Runge-Kutta (MLRK). To this end we conducted experiments with different classes of initial value problems. Specifically, we Table 1: On the test set of class (low-freq), MLRK was considerably faster (less many steps) than the baseline controller on average, while causing only a mildly larger mean global error at the end of the integration interval. Both the baseline controller and MLRK showed a gradually increasing mean global error at the end of the integration interval over the test set consisting of 1500 harmonic oscillators of class (med-freq). investigated our suggested loss functions and our theory that suggests that providing higher order partial derivatives of g can be used as a means to anticipate the evolution of the local error and, hence, may lead to an improved performance of the optimizer. For a class of initial value problems we assume a parametric form of g with a distribution over the parameters; details can be found in the appendix. Our datasets are obtained by sampling from these distributions. Finally, we compare our method to the baseline controller given by Eq. ( 2).

Harmonic Oscillator. First, we started with a simple class of ODEs. We trained MLRK on low varied harmonic oscillators. The corresponding results for training on high frequencies are similar and can be found in the appendix. As inputs we only used the error estimates, Eq. ( 7). We started with the dataset of harmonic oscillators with low frequencies (low-freq): the training set contained 30,000 instances, the validation set and test set each 1,500 instances. To measure the performance we considered the L 1 -Lagrangian loss Eq. ( 5) as well as the number of steps needed for an integration interval and the corresponding error separately. Since the Lagrangian represents a trade-off between the two it is a particularly good measure.

Fig. 2 (a) shows the mean loss during the integration of the test set, which consists of 1500 harmonic oscillators of class (low-freq). The solid lines show the mean loss in step n, the colored areas mark the standard deviation. MLRK achieves a smaller mean loss than the baseline, which confirms its ability to generalize to new problem instances of the same class of problem. When we compare the mean number of steps and the mean global error of MLRK and the baseline controller in Tab. 1, the necessity for a trade-off between the two objectives becomes clear. While the baseline needs on average more steps to complete the integration, it achieves a smaller error, the contrary is true for our optimizer. This makes it hard to compare the two methods based on just these values.

Generalization Capability of Optimizer based on Low Frequency Harmonic Oscillators. By evaluating the controller designed by MLRK on different test sets, we investigated its ability to transfer knowledge to problem instances of other classes. Fig. 2 (b) shows that MLRK maintains a smaller mean loss on 1500 harmonic oscillators of class (med-freq). These oscillators are both higher in frequency and amplitude than the problem instances contained in the training set. This indicates that MLRK is able to generalize to these problems. Tab. 1 reveals that MLRK uses on average about a third of the integration steps the baseline needs for the tested integration intervals. The resulting mean global errors are indeed larger than those of the baseline controller, but still very reasonable. Moreover, evaluating MLRK on 1500 harmonic oscillators of class (high-freq) showed a mean loss similar to that of the baseline controller in Fig. 2(c ). This suggests that MLRK generalizes to problem The optimizer err was trained with input (7), optimizer partial was trained with input (8) and grad was trained with input (9). Providing local information such as those in input ( 8) and ( 9) results in increased responsiveness of the optimizer.

instances of this class as well. In Tab. 2 we observe a similar situation as in the previous experiment.

The global error of the baseline controller stays approximately the same for different interval lengths, while the global error of MLRK increases very gradually, while using only a fraction of iterations.

van der Pol. "I have a theory that whenever you want to get in trouble with a method, look for the van der Pol equation" -P. Zadunaisky, 1982, c.f . (Hairer et al., 2000) .

Next we conducted an experiment regarding the local L 1 loss (6). Experiments with harmonic oscillators revealed the optimizers ability to keep the local errors close to the tolerance parameter (c.f . appendix). Keeping the local error equal to a constant tolerance is not a hard problem for simple harmonic oscillators, therefore we turned towards van der Pol oscillators. We considered the following dataset of van der Pol oscillators of class vdP(0, 1): The training set contained 50,000 instances, the validation set and the test set each 1,500 instances. The local errors of van der Pol oscillators vary drastically (c.f . Fig. 5 , appendix). High spikes occur when the solution changes from being driven to being damped. As the behavior of van der Pol equations are more complex compared to harmonic oscillators, we expect that providing additional information about the problem such as partial derivatives or function values of g can be beneficial. To investigate this, we considered MLRK with different inputs. One optimizer is only provided with the error estimate Eq. ( 7), another is provided with the error estimate and all partial derivatives Eq. ( 8) and the last is provided with the error estimate and the function evaluations that were computed during the Runge-Kutta step, Eq. ( 9).

The results are summarized in Fig. 3 . Using err shows reoccurring high spikes similar to the baseline, although the spikes of the optimizer are lower on average. In contrast, MLRK with additional information on the partial derivatives (partial) shows a much smoother sequence of local errors. This demonstrates that the additional information in the partial derivatives of g improves the ability of MLRK to anticipate the variations of the local errors and respond with adequate step sizes. MLRK with error estimates as well as function values of g (grad) produces a smooth step size sequence Under review as a conference paper at ICLR 2020 The optimizer err was trained with input (7), optimizer partial was trained with input (8) and grad was trained with input (9). In all three cases MLRK is able to reduce the spikes in the local errors drastically, the use of local information reduces the variance of the local errors even further.

similar to the optimizer partial. This shows that the optimizer with additional input-information from g are actually outperforming the baseline. The reoccurring spikes in the local error seem to decrease in amplitude with time. This is largely due to an averaging effect -for different values of σ these spikes occur at different integration steps. All three methods reveal a reduction in both the mean number of steps and the mean local errors, c.f . Tab. 3, appendix.

Double Pendulum. To demonstrate the benefit of MLRK in an application we conducted experiments with double pendulum ODEs (c.f . appendix). The chaotic behaviour of the double pendulum leads to sudden spikes in the local error, a particularly challenging problem for step size controllers. To investigate the ability of MLRK to avoid these spikes by adjusting the step size appropriately, we consider the L 1 local loss as in the previous experiment and a dataset of double pendulums of class (pendulum). The training set contained 50,000 instances, the validation set and the test set each 1,500 instances. Again, the different inputs Eqs. ( 7), ( 8) and ( 9) were examined.

The results are shown in Fig. 4 and reveal the baseline controllers inability to produce steady local errors. Here, the benefit of a learned update rule over a static one becomes evident. The optimizer err that is only provided with the local error estimates is able to reduce the spikes in the local errors to a remarkable degree. Moreover, local information about the problem allows the optimizers partial and grad to reduce the variance in the local errors even further.

5. CONCLUSION:

We have shown how to cast the design Runge-Kutta methods for solving initial value problems as a learning problem. We established appropriate performance measures and useful inputs for the controller, which constitute the key ingredients of the resulting Meta-learning Runge-Kutta (MLRK), which learns step size controllers that are specialized to a particular class of initial value problems.

Our experimental results demonstrate that MLRK can indeed learn to design novel Runge-Kutta methods that perform better than a hand-designed Runge-Kutta approach. Furthermore, we observed a remarkable degree of generalization to other classes of initial value problems. More importantly, examining the effect of different inputs demonstrated that the additional information contained in the function values and partial derivatives of g leads to a substantial improvement in performance of the automatically designed solvers.

There are several interesting avenues for future work. While MLRK generalizes well to problem instances of the same class and even to problems of similar classes, when confronted with a very different type of problem, MLRK does not generalize well yet. Wichrowska et al. (2017) showed how a carefully chosen network architecture and a diverse training set improves generalization of their optimizer to many different classes of optimization problems. A similar approach may lead to improved generalization of MLRK. Another common aspect of Runge-Kutta is the need for step rejection in case the local error exceeds the tolerance. The problem here is that a single large error can in general not be compensated for, even if subsequent step sizes are chosen very small. In this case, a step is usually rejected and repeated using a smaller step size. Here the challenge is to choose a step size small enough to meet the accuracy requirements but at the same time not too small, since the choice of step size influences the entire step size sequence to come.

A APPENDIX RUNGE-KUTTA METHODS:

The mathematical theory of Runge-Kutta methods is quite evolved. We will only recap a few basics here, for further details see for example (Butcher, 2008; Hairer et al., 2000; Hairer & Wanner, 2002) .

Definition A.1. Let s ∈ N, a ij ∈ R for i, j ∈ {1, . . . , s}, c i ∈ R for i ∈ {1, .

. . , s} and b i ∈ R for i ∈ {1, . . . , s}. Let g : R × R m → R m be a function that describes an initial value problem. The method defined by

EQUATION

)

EQUATION

for step size h is called an s-stage Runge-Kutta method. If a ij = 0 for i ≤ j, the method is called explicit, otherwise it is called implicit.

Definition A.2. A Runge-Kutta method is of order p if for sufficiently smooth initial value problems

y(t 0 + h) -y 1 ≤ Kh p+1

holds.

Note that a higher order method yields more accurate results, consequently a high order is a desired characteristic of a Runge-Kutta method. However, higher orders can only be achieved by the use of more stages. For an s-stage Runge-Kutta method the order p is bounded by the number of stages, p ≤ s. Up to order 4 there exist methods with p = s, for order 5 and higher p is strictly smaller than s, c.f . (Butcher, 2008, Theorem 324B) . The order of an s-stage Runge-Kutta method depends on the coefficients a ij , b i , c i . Order conditions for these coefficients have been developed, which ensure a certain order of the method.

The problem of step size control can be described in the following way: Based on some input values, e.g. t n , y n , err n and possibly additional characteristics of g, we must choose a step size h n ,which is then used to evaluate g at p different points. These points are determined by h n , see Definition A.1. The resulting stages in Eq. ( 15) are then used to compute the approximation y n+1 of y(t n + h n ), c.f . Eq. ( 14). The problem is to choose h n in a way such that the approximation y n+1 fulfills some desired properties. For example we may require that the local error of y n+1 is close to some tolerance parameter. In fact, this is the standard objective of most common step size controllers.

ERROR ESTIMATION:

The local errors are of particular importance for step size control and can fortunately be estimated very efficiently. To that end, we take a look at a simple idea as described in (Butcher, 2008, p. 198 ): Suppose we have two approximations for y(t n ) of order p and p respectively, that is

ŷn = y(t n ) + O(h p+1 ) , ỹn = y(t n ) + O(h p+1 ) ,

where y(t) is the solution of

y (t) = g(t, y(t)) , y(t n-1 ) = y n-1 ,

and the step size h that was used to obtain both approximations is given by

h = t n -t n-1 . If p > p, ŷn -ỹn = y(t n ) -ỹn + O(h p+2 )

can be used as an approximation for the local truncation error of the p-th order Runge-Kutta method:

err n (h) ≈ ŷnỹn .

To obtain error estimates each step has to be carried out with two methods of different orders. Here, the easiest approach is to use two separate Runge-Kutta methods to compute ŷ and ỹ, however this is quite costly, especially for implicit methods. A more efficient approach to obtain error estimates are the embedded Runge-Kutta methods. Here both approximation use the same stages k i in Eq. ( 15) but two different pairs of coefficients b i , bi in Eq. ( 14). By doing so, the stages can be reused and an error estimate can be obtained with a cheap extra linear combination of the k i 's. We denote the error estimates by

EQUATION

CLASSES OF INITIAL VALUE PROBLEMS:

For loss functions such as the one in Eq. ( 5) it is useful to consider initial value problems with known solutions so that we can compare the numerical approximation to the true function values of the solution. For this reason we considere simple harmonic oscillators and linear differential equations with constant coefficients. If we employ loss functions as Eq. ( 6), local rather than global errors are of interest. For these loss functions it is interesting to consider initial value problems for which the local errors vary drastically during the integration such as van der Pol oscillators.

SIMPLE HARMONIC OSCILLATOR:

A simple harmonic oscillator is a harmonic oscillator that is neither driven nor damped and can be characterized by m x (t) = -k x(t), where m is the mass of the oscillator, x the position and k describes the restoring force that is applied, when displaced from its equilibrium position. It can be transformed to

y 1 (t) = y 2 (t) , y 2 (t) = - m k y 1 (t) .

The initial values provide the initial position x(t 0 ) of x and the initial velocity x (t 0 ),

y 1 (t 0 ) = x(t 0 ), y 2 (t 0 ) = x (t 0

). The solution of the simple harmonic oscillator takes the form

x(t) = A cos(ωt + ϕ), ω = m k ,

where A and ϕ are uniquely determined by the initial values. The frequency of the oscillations is determined by ω, A is the magnitude of the oscillations and ϕ is a phase-shift.

A class of this kind of initial value problem can be described by a distribution over m k and the initial values or alternatively a distribution over A, ω and ϕ. The following classes of initial value problems were used in the experiments

EQUATION

A ∼ U (0, 5), ω ∼ U (0, 20), ϕ ∼ U (0, 1), (higher-freq)

EQUATION

A ∼ U (1, 5), ω ∼ U (1, 5), ϕ ∼ U (0, 1), (med-freq) where U (a, b) denotes the uniform distribution over the interval [a, b] .

LINEAR CONSTANT COEFFICIENT:

A linear differential equation with constant coefficients can be described by y (t) = A y(t) , y(t 0 ) = y 0 , where A ∈ R m×m . The solution of this differential equation is known and can be expressed in terms of the eigenvalues λ 1 , . . . , λ m with corresponding eigenvectors v 1 , . . . , v m of the matrix A:

y(t) = m j=1

c j e λj t v j , with coefficients c j ∈ R that are uniquely determined by the initial values y(t 0 ) = y 0 . The initial value problem is stable, if Re(λ j ) < 0 for all j ∈ {1, . . . , m} and unstable if Re(λ j ) > 0 for some j ∈ {1, . . . , m}. The oscillatory behavior of the problem is determined by Im(λ j ). Note that the simple harmonic oscillators are linear differential equations with constant coefficients, where

A ∈ R 2×2 , Re(λ j ) = 0 for j ∈ {1, 2} and Im(λ 1 ) = -Im(λ 2 ) = m k i.

For the experiments we used differential equations with A ∈ R 3×3 . The matrix A either has three real eigenvalues λ 1 , λ 2 , λ 3 ∈ R or one real eigenvalue λ 1 ∈ R and two complex eigenvalues λ 2 = a + bi, λ 3 = abi, with a, b ∈ R. In this case we choose c 2 = c 3 to obtain a real solution. We can specify a class of linear constant coefficient differential equations by a distribution on the eigenvalues of A and the coefficients c j . The classes that were used for the experiments were described by the distributions 

EQUATION

EQUATION

EQUATION

VAN DER POL OSCILLATOR:

Van der Pol oscillators arise in the study of limit cycles. A van der Pol oscillator can be described by y 1 (t) = y 2 , y 1 (0) = 2, y 2 (t) = σ(1y 2 1 )y 2y 1 , y 2 (0) = 0, with σ > 0. For σ = 2 an approximate solution is displayed in Figure 5 . Small oscillations are amplifies and large oscillations damped, c.f . (Hairer et al., 2000, pp. 111) . The local errors of a van der Pol oscillator show high variations. For the experiments we used van der Pol oscillators with the following distributions σ ∼ U (0, 1), (vdP(0, 1))

σ ∼ U (1, 2). (vdP(1, 2)) DOUBLE PENDULUM

A double pendulum is a pendulum attached to another pendulum, c.f . Fig. 6 . The length of the first and second pendulum are given by L 1 and L 2 , respectively, their masses by M 1 and M 2 and the angle by θ 1 and θ 2 . The equations of motion for the double pendulum can be written in different forms, we chose the following formulation in terms of the angular acceleration:

θ1 = M 2 L 1 θ2 1 sin(θ 2 -θ 1 ) cos(θ 2 -θ 1 ) + M 2 g sin(θ 2 ) cos(θ 2 -θ 1 ) L 1 (M 1 + M 2 -M 2 cos 2 (θ 2 -θ 1 )) + M 2 L 2 θ2 2 sin(θ 2 -θ 1 ) -(M 1 + M 2 )g sin(θ 1 ) L 1 (M 1 + M 2 -M 2 cos 2 (θ 2 -θ 1 )) θ2 = (M 1 + M 2 ) g sin(θ 1 ) cos(θ 2 -θ 1 ) -L 1 θ2 1 sin(θ 2 -θ 1 ) -g sin(θ 2 ) L 2 (M 1 + M 2 -M 2 cos 2 (θ 2 -θ 1 )) - M 2L 2 θ2 2 sin(θ 2 -θ 1 ) cos(θ 2 -θ 1 ) L 2 (M 1 + M 2 -M 2 cos 2 (θ 2 -θ 1 ))

The initial angular position and velocity θ 1 , θ 2 , θ1 , θ2 determine the trajectory of the double pendulum.

For the experiments we used double pendulums with the following distribution M 1 ∼ U (0.5, 1), M 2 ∼ U (0.5, 2), L 1 ∼ U (0.5, 1), L 2 ∼ U (0.5, 1).

(pendulum)

ADDITIONAL EXPERIMENTS:

ON VAN DER POL Table 3 shows the mean number of steps, local error and wall clock time over 1500 van der Pol equations of class (vdP(0, 1)).

L 1 L 2 M 1 M 2 θ 1 θ 2

Figure 6 : A double pendulum is one pendulum attached to another pendulum. The masses M 1 and M 2 of these pendulums may vary as do the lengths L 1 and L 2 . The initial configuration of angular position and velocity determines the trajectory of the double pendulum.

HIGH FREQUENCY HARMONIC OSCILLATORS:

In this experiment, we trained a controller specialized for simple harmonic oscillators of class (high-freq) which were described in Sec. A. Thus we chose the following datasets:

Training set: 30000 Harmonic oscillators of class (high-freq). Validation set: 1000 Harmonic oscillators of class (high-freq).

Test set: 1000 Harmonic oscillators of class (high-freq).

As loss we chose the L 1 -Lagrangian Eq. ( 5). As inputs we used only the error estimates, Eq. ( 7). The loss of the learned controller and the baseline controller, averaged over the test dataset, is shown in Fig. 7(a) . Both controllers are used to execute 100 integration steps for each instance of the test dataset. Subsequently, the loss Eq. ( 5) is computed for each approximation (t n , y n ). The solid lines show the mean loss in step n, the colored areas mark the standard deviation. Here, the loss of our learned controller is on average smaller than the loss of the baseline controller, which demonstrates the controllers generalization capability to new problems of the same class.

Tab. 4 summaizres the global error at the end of the integration interval and the number of steps needed during the integration for different interval lengths in range 1 to 10. All values are averaged over the test dataset. Here, the baseline controller needs much more steps than the learned controller, but maintains a smaller error. In fact, the error caused by the baseline is approximately the same for all shown interval lengths, while the error caused by the trained model increases with interval length. This is due to the assumption that the global error is the sum of the local errors. In particular, the loss Eq. ( 5) rates step size sequences with a higher global error better, if they can balance the higher error by achieving a longer integration interval.

While the general assumption that the global error increases with the length of the integration interval is often fulfilled, the specific assumption of a linear accumulation of the global error, as in our loss, is usually not fulfilled. This is a good example for why the assumptions leading to the widely used objective to keep the local errors constant are too strong. It is also interesting to evaluate the learned controller on other classes of initial value problems. Fig. 7 (b) shows the mean loss of the learned controller and the baseline controller over a test set consisting of 1500 harmonic oscillators of class (higher-freq). This class contains oscillators of class (high-freq) but also includes higher frequency ones. The learned controller has a lower mean loss than the baseline controller which indicates that it is able to transfer knowledge from problems of class (high-freq) to problems of class (higher-freq).

When we compare the number of steps and the global error at the end of the integration interval, see Tab. 5, we can observe a similar situation as in Tab. 4. The learned controller needs less steps than the baseline controller. The global error however stays approximately the same for the baseline controller while it increases with interval length for the learned controller. Figure 7 : The red and blue line show the mean loss Eq. ( 5 The mean number of steps and the mean global error at the end of the integration interval was computed over a test set of 1500 harmonic oscillators of class (higher-freq). While the baseline controller needs much more steps than our learned controller, it has on average a lower global error, which stays approximately the same for different interval lengths.

HARMONIC OSCILLATOR:

We learn a controller optimized for simple harmonic oscillators. We choose the following data sets:

Training set: 30000 Harmonic oscillators of class (low-freq). Validation set: 1500 Harmonic oscillators of class (low-freq).

Test set: 1500 Harmonic oscillators of class (low-freq).

First, we evaluate our method on a test set consisting of the same class of initial value problems the controller was trained on. Figure 8 (a) shows the mean local error of the learned controller and the baseline on this test set. Both controllers are capable to keep the local errors close to the tolerance.

Our method is slightly closer to the tolerance. Next, we evaluate our learned controller on different test sets. When we test the learned controller on harmonic oscillators with higher frequencies and higher amplitude, we find that our method generalizes to these problem instances as well. In fact, both the learned and the baseline controller are able to keep the local errors close to the tolerance as displayed in Figure 8 (b). Similar results were found for harmonic oscillators of the class (med-freq). Hence, our learned controller is able to transfer knowledge from oscillators with low frequencies to oscillators with higher frequencies. For these problem instances the learned controller demonstrates its ability to keep the mean local error close to the tolerance. In fact, the local errors of our controller appear to be closer to the tolerance than the local errors of the baseline for some part of the integration. (b) Testing our controller on a set of 1500 van der Pol oscillators of type (vdP(0, 1)) reveals a similar performance to that of the baseline controller.

Both controllers show high spikes in the local errors. These occur whenever the solution changes from being driven to being damped, see Figure 5 .

We also want to evaluate the generalization capability of our method to problem instances of different classes of initial value problems. Our method shows an acceptable performance on linear differential equations with constant coefficients of the class (osc-increasing), see Figure 9 (a). Here, the performance of our controller is similar to the baseline controller. The learned controller obtains local errors which are closer to the tolerance in some steps. Figure 9 (b) shows a similar performance of our controller to the baseline controller on van der Pol oscillators. Both methods show high spikes in the local errors. These occur whenever the solution of a van der Pol oscillator changes from being driven to being damped as demonstrated in Figure 5 . Our controller does not react quickly to the sudden changes in the local error of van der Pol equations. This behavior is very different from harmonic oscillators, thus generalization can not be expected.

:



Figures and Tables:

Figure2: The mean L 1 -Lagrangian loss (5) of the approximation (t n , y n ) in step n over three different test sets. (a) MLRK achieved a smaller mean loss, thus, it generalizes well to new problem instances of the same class of problems (low-freq) it was trained on. (b) MLRK achieves a lower loss on instances with both higher amplitude and frequency than the problem instances used during training. Thus, MLRK is able to generalize to different problems. (c) Even on higher frequencies than the ones in (med-freq), MLRK shows a mean loss similar to that of the baseline controller.

Figure 4: The mean local error over the test set of double pendulums for different methods is shown.The optimizer err was trained with input (7), optimizer partial was trained with input (8) and grad was trained with input (9). In all three cases MLRK is able to reduce the spikes in the local errors drastically, the use of local information reduces the variance of the local errors even further.

) where U (a, b) denotes the uniform distribution on the interval [a, b].

Figure5: The upper plot shows the numerically approximated solution y 1 (t) of a van der Pol oscillator with σ = 2. The lower plot displays the local error estimates for h = 0.1. The error estimates show high spikes whenever the solution switches from being driven to being damped.

Figure 8: The mean local error of our learned controller and the baseline controller is shown. (a) Testing on 1500 harmonic oscillators of class (low-freq) shows that both the baseline and the learned controller are able to keep the local errors close to the tolerance. The mean local error of the learned controller is slightly closer to the tolerance. (b) For a test set consisting of 1500 harmonic oscillators of class (high-freq) the two controllers show a similar performance, both are able to keep the local error close to the tolerance. Testing the controllers on problem instances of class (med-freq) yields a similar results.



Table content: <table><tr><td>t n , y n</td><td>t n ,y n</td><td>t n+1 ,y n+1</td><td>t n+1 ,y n+1</td></tr><tr><td/><td>logh n 1</td><td/><td>logh n</td></tr></table>



Table content: <table><tr><td>: The mean global er-rors of the baseline stay ap-proximately the same. The optimizer trained on problem instances of class (low-freq) shows only a very gradual in-crement in the global error on instances of class (high-freq).</td></tr></table>

The mean number of steps, local error and wall clock time over 1500 van der Pol equations are shown for different lengths of the integration interval. Our method err is slightly faster and uses less steps than the baseline while producing smaller local errors during the integration, see also Figure3. While partial and grad reduce the number of steps even further, wall time is increased. However one can clearly see that partial and grad outperform both the baseline and err regarding the local error.

Table content: <table><tr><td>int</td><td/><td/><td>steps</td><td/></tr><tr><td/><td>baseline</td><td/><td>err</td><td>partial</td><td>grad</td></tr><tr><td>1 3 5 7 10</td><td>21.59 33.74 45.43 56.84 73.46</td><td/><td>16.42 29.12 41.07 52.57 69.41</td><td>12.40 25.16 36.42 48.34 65.40</td><td>12.09 24.74 36.02 47.97 65.04</td></tr><tr><td>int</td><td/><td/><td>mean local error</td><td/></tr><tr><td/><td>baseline</td><td/><td>err</td><td>partial</td><td>grad</td></tr><tr><td>1 3 5 7 10</td><td>7.17e -4 6.40e -4 5.18e -4 4.97e -4 4.59e -4</td><td/><td>6.58e -4 5.45e -4 4.47e -4 4.16e -4 3.82e -4</td><td>4.01e -4 2.49e -4 1.95e -4 1.57e -4 1.32e -4</td><td>3.74e -4 2.28e -4 1.75e -4 1.39e -4 1.18e -4</td></tr><tr><td>int</td><td/><td/><td>time</td><td/></tr><tr><td/><td>baseline</td><td/><td>err</td><td>partial</td><td>grad</td></tr><tr><td>1 3 5 7 10</td><td>0.0255 0.0405 0.0591 0.0858 0.0971</td><td/><td>0.0263 0.0375 0.0517 0.0742 0.0825</td><td>0.0254 0.0460 0.0681 0.1036 0.1201</td><td>0.0221 0.0403 0.0596 0.0897 0.1065</td></tr><tr><td/><td>interval</td><td>steps</td><td/><td>error</td></tr><tr><td/><td/><td>Baseline</td><td>Optimizer</td><td>Baseline</td><td>Optimizer</td></tr><tr><td/><td>1 3 5 7 10</td><td>45.50 152.11 258.73 365.30 525.18</td><td>13.73 43.12 74.40 106.31 154.25</td><td>0.023092 0.022266 0.022731 0.022011 0.023353</td><td>0.030499 0.040429 0.049529 0.055821 0.070451</td></tr></table>

The mean values over the test set consisting of 1000 harmonic oscillators of class (high-freq) are shown. While the baseline needs much more steps than our trained model, it achieves a smaller global error. The global error of the baseline stays approximately the same for all integration intervals.

Table content: <table/>

) of the approximation (t n , y n ) in step n over the test set. The red and blue area indicate the standard deviation. (a) The learned model LSTM attains a smaller mean loss on new instances of the same class of problems (high-freq) it was trained on, this is especially pronounced in the steps n ≥ 5. (b) The test set consists of 1500 harmonic oscillators of class (higher-freq), which means that it also contains higher frequency oscillators. The mean loss of our learned controller is lower than the mean loss of the baseline controller which indicates that our method generalizes to these different problem instances.

Table content: <table><tr><td>interval</td><td>steps</td><td/><td>error</td><td/></tr><tr><td/><td>Baseline</td><td>Optimizer</td><td>Baseline</td><td>Optimizer</td></tr><tr><td>1 3 5 7 10</td><td>116.76 386.40 656.06 925.80 1330.31</td><td>22.26 75.64 141.29 212.19 319.85</td><td>0.453698 0.442521 0.426773 0.418467 0.413902</td><td>0.511174 0.620779 0.770457 0.925463 1.427778</td></tr></table>

Peer review (iclr2020):
## Summary ## The authors present a method for learning a step-size adaptation strategy for Runge-Kutta methods. The principal contributions of the paper are: 1. They define a loss function that better captures global performance of the controller, rather than just local behavior. 2. They propose a set of input features to a step size controller. It includes the intermediate Runge-Kutta evaluation values, which allow the controller to approximate the derivatives of the Jacobian function . 3. They describe a recurrent architecture for adapting step size over time. Runge-Kutta methods are a workhorse of ordinary differential equations and choosing step size is one of the central challenges involved in their application. Better methods for step size selection would definitely be of broad interest. As the authors point out, existing methods often consist of hand-tuned heuristics---a feature that often suggests machine learning could provide significant improvements. While the premise of the paper is very promising, I don*t think it is ready to be accepted to ICLR at this time. Most significantly, the experimental results are not particularly compelling. I believe the authors should refine their method, aim for better experimental results and resubmit. I have included more specific comments below. ## Specific Comments ## 1. I think the paper would benefit from a clearer description of the RK step-size selection problem. For instance, for a p^th order RK solver, at each time step, * Inputs: t, y(t), g # Also possibly intermediate values from previous time steps. * Select a step size h(t) * Evaluate g at p different points based on h(t). * Use these evaluations to compute a value of y_(t + h_t) For those that aren*t familiar with these methods (at ICLR there will be many!) I think this would help explain where the authors* method (and the other methods your describe) fits into the larger algorithm. 2. I think a short explanation of error estimation would be helpful in addition to the reference to Butcher. This estimation is critical to step size adaptation. In particular I think this would be clearer if the authors expanded the paragraph at the bottom of page 4, where they describe error in a polynomial in (t_n, h) whose coefficients are derivatives of g. 3. The authors* proposed loss function (Eq. 5) includes the true value of y at t_n, y(t_n). They acknowledge that this may make it computationally prohibitive, but I think this point warrants further discussion. Does this mean that their loss function can only be used on problems for which we have a closed form solution (such as harmonic oscillators)? I noticed that in their van der Pol experiments, the authors switch to a more standard loss (Eq. 6). Is this because Eq. 5 is intractable in this example? Is there a reasonable approximation to Eq. 5 that could be used when a closed form solution is not known? 4. There are a number of issues with the experiments that I think could use clarification or improvement: (a) The authors compare to a *baseline* but I don*t believe this baseline is defined anywhere. Is it one of the adaptation methods described in Section 2? (b) In Table 1, the baseline method achieves lower error, while MLRK uses fewer steps. It is difficult to assess if this is an improvement since this this cost-accuracy tradeoff is at the heart of the Lagrangian formulation. Ideally, shouldn*t we be able to adjust *tol* to trace out some Pareto frontier for the cost-accuracy tradeoff? In this case, wouldn*t we hope for MLRK to be able to achieve better accuracy given the same computational budget? (c) In the van der Pol experiments, the authors switch to local L1 loss (Eq. 6). Was the intention to experiment with local L1 loss and van der Pol provided an interesting class of examples? Or was the reasoning that Eq. 5 is intractable for van der Pol so they had to use local L1 loss? If Eq. 5 can still be evaluated on these experiments, this might be a more convincing comparison. (d) The number of steps required was provided for the harmonic oscillator experiments, but not the van der Pol ones. This would be helpful for comparing the methods. (e) What is the computational overhead of running an RNN alongside your solver? Although it doesn*t tell the whole story, it would be informative to report wall clock time along with the number of steps required for each method.