Meta-Learning Runge-Kutta
Authors: Nadine Behrmann
Patrick Schramowski
Kristian Kersting
Venue:  (2020)

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. In this paper, 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. The source code can be found at https://www.dropbox.com/sh/rkctdfhkosywnnx/AABKadysCR8-aHW_0kb6vCtSa?dl=0

Under review as a conference paper at ICLR 2020
META- RUNGE-KUTTA
Anonymous authors Paper under double-blind review
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-Kutte Kutta methods. They are, however, sensitive to the strategy the step size is adapted during integration, which has to be chosen by the exper- imenter. Here, we show how the design of 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 an automatic way. The key ingredients for the resulting Meta-Learnin
Runge-Kutta Kutta (MLRK) are the development of 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 mitigate sudden spikes in the local error function by faster adaptation of the step size. More importantly, the additional information in the form of partial derivatives and function values leads a substantial improvement perfor- mance. The source code can be found at https: www. wwwp:wwwwwooox.ooooooo com/ sh/ Besttd/hhkkkyymmx//AABaaea
aHW 1tW__00kb6oCttaadll1
INTRODUCTION
Differential equations in their general form cover 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, (973). 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 (llea et al. 2013). diabetes and the glucose metabolism (Esna-Ashari et al.) 2017) as well as for pharmaceutical drug design (Deufthard, 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, 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- 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 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 learning problem, c.f. Fig. [] The key ingredients of the resulting Meta-Leerring Runge-Kutte (MIRK) are the identification of good performance measure and appropriate inputs.
Related Work. Various approaches to control the step sizes in Runge-Kutte Kutta methods have been proposed. While some rely on signal processing principles where the goal is to produce 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. 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 single step close tolerance parameter.

---

Under review as a conference paper at ICLR 2020
Eu M
logha
loght
Runge-Kutte update
Runge-Kutar update
Optimizee
Model Input
Model Input
log Fn
log
Lin
Lin
Optimizer
LSTM
LSTM
Z
2x++1
Figure Meta-Learning Kutte- The model represented by the blue block determines the step size adjustment log ho-1 using LSTM and linear layer, c.f. Eq. (1O). The Runge-Kutte update then determines the new step size log h = log hn + log fn-1- The new step size used perform the next Runge- Kutte-Kuttaaaaa step which computes the approximation yn++1 time step 4u+1 and an error estimate err,+1 according to (12). Then the next input Un+1 for the model is computed.
Specifically, the design of the step size control algorithm is learned within MLRK. Indeed, casting algorithm design as learning problem not new, Andrychowicz et al. (2016) learned gradient- based 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 projection- free convex optimization algorithm, which relies on conditional gradients. Finally, Chen et al. (2017)
cast the design of gradient-free black-bock optimization
learning problem.
Furthermore, other approaches to solve differential equations using neural networks have been proposed, too. Lagaris al. (1998) suggested to use 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 deep neural network directly approximate the solution of a high-dimensional partial differential equation. [E & Yu (2018) proposed the Deep Ritz method
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 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 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, 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-KUTTH METHODS
We give a brief introduction to initial value problems and Kunge-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' = =gt y), y(tu) = yo with a function g [a, b) I be very hard to find, numerical methods such as Runge- Kutta methods have been developed. Their mathematical underpinning is quite evolved (Butcher) 2008
RT m N. As closed-fom solutions
(Hiirer et al. 2000
Hairer & Wanner 2002). A brief recap can be found in the appendix.
2

---

Under review as a conference paper at ICLR 2020
For iterative solvers, the local truncation error the error committed in 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 + is defined as e (h) lly(t,, + h) Mn++il C(t,)) + O
(1)
where Un+1 (t Vu, h)) is single Runge-Kutta step of p-th order method using step size h. We call (t.,) the principal error term.
We will also denote the local truncation enror local error and refer to the appendix for simple ways to estimate We denote the error estimates as emr,
The Objective of Step Size Control. As numerical integration is used approximate solution, high accuracy is a desired property. However, in practice, computational resources might be limited and efficient use of these resourees 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 ((H) and the work cost W (H)) produced by a sequence of step sizes H = (hn N-1 E(H) + À W(H). Butcher (008) uses two integrals describing the error E(H) and the work cost W WH). Under some assumptions one can show that optimal step sizes w.r.t. this objective causes local errors equal to the Lagrange multiplier A. which is often referred to as the tolerance parameter. Namely, we need to assume that the global error equal the sum of local errors and all step sizes are small enough so that the loss can be approximated by the integrals suggested by BBtcher 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 equal to the sum of the local errors is not true in general.
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 eиror for step size h at time is approximately emr C(t)APf)
and, accordingly. the desired optimal step size hope is given by hopt sJ h (tol/err) p+1, This formula then used to adapt the step size
tol
p+1
hn+1 = Fi
Tn
max a, min 3,7
err
where safety factor 0 < < and minimal and maximal factor a, 3 R I 8 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 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) dis- (PID) controller. The classical step size control mechanism in Eq. (2) can be regarded as an I-controller. They then demon- strated 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. (I988) suggested a Pl-controller, which in addition to the integral term used in standard step size control, also includes a proportional term and can be expressed as
cussed step size control in the context of a
tol err,+t1
tol err,
))
h = /
F'
max
a, min
8,7
(3)
Predictive Controller. Both the classical step size controller, see Eq. ( and the Pl-controller in Eq. (3) rely on the assumption that the local truncation error can be described as function of h that remains independent
This is certainly not true for all cases. Therefore, Gustafsson Therefore,GGustafsonn (1994). e.g.. discussed a controller based on 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
hn
tol
CIT
hn-1
h
(4)
Tn--1
eIT,+1
eIr,
3

---

Under review as a conference paper at ICLR 2020
3 META-LEARNING RUNGE- KUTTA: NEURAL STEP-SIZE SIEE-SIZE CONTROLLER
Step size controllers for Runge-Kutta are typically still designed by hand. Formulas such as Eq. (2) (3) and 44) require parameter fine-tuning, In this work, we take 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 controller learning problem is to determine good performance measure and appropriate inputs.
Performance measure. As the first step we need to identify 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, to into account as well, which leads to the following performance measure:
(t.) yn + tol tx to In contrast, if we follow the concept of most step size controllers, optimal step size is achieved if the local error is equal to the tolerance, which suggests the performance measure
[Lagramge (t Yn: err,))
(5)
plocal Vne eIr,)) = ||eTr, - toll|-
(6)
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. 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 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 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 update the step size. Specifi- cally, tol. 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,
tol
v (errm+ = log
eer
Here, vo () 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 eτrors can be expressed with elementary differentials of the ODE, off. (Hairer et al. 2000, Section II.3). These elementary differentials are formed of function values and partial derivatives of 9 and allow a Taylor approximation of the local error function e (h) h)) f+1 + (h*+2), where p(tn, h) is a polynomial in h with coefficients formed from partial derivatives of g up to order k > p. The roots of the polynomial ((n,, h+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. it clear that the optimal step size, real root of p(h), depends continuously on the partial derivatives and function values of 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 ( (tt, Vn)) al= ..... Ög(tns Wn)), g( yn), tol). This suggests that the partial derivatives of y can also be an appropriate input for our controller,
vp (err,+1;
log
tol
(0 gtn Fn) (alpp ..... Ig(tn, Ba), g(tns (lnneBB))))))))))) Un))
(8)
err
4

---

Under review as a conference paper at ICLR 2020
However, providing partial derivatives requires additional computation whereas the 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:
tol
+ (eIIx+1,))
log
(1)
(p)
eff
where gh
9ri denotes function values of the initial value problem.
In contrast the traditional input 77, our novel inputs and ( both provide local information about the structure of the problem at hand, which allows one to keep the local emror 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 or as shown next, does so automatically and designs novel Kunge-Kutta methods.
Meta-Learner. With 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. | Let the input for the optimizer be 2n + (eIr,, -). We parameterize the optimizer using Φ and update the step size as follows, log ha+1 = log ht c, (e'(ern,+1:)) :), 6) An approach in log space has the advantage that that we do not need to constrain the output of to be positive. Since C2 corresponds to log r,, 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 optimizer e, 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 hn It can be expressed by
log (ogTn-1)
(Inz hn Ø)
(10)
hx+1
The subsequent Sunge-KKuter update in the yellow block of Fig. [first updates the step size
log hn log hn-1 + logrn-1. :
(11)
and then uses execute the next Runge-Kutta step,
In 1 bn + hn
yn+1 eer
(bon Moy hm)
(12)
Here, t denotes the Runge- Kutta step from Yn to yu+11 as well as the error estimate emf 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 described by function g, therefore we can represent class of initial value problems with a distribution over the functions Thus, an optimizer can be considered optimal for class of problems, if minimizes the expected loss:
nー Σ>,"/
N
[(o) R.g
1(Ins Yns (rr,),re)))))
(13)
where in. yn, err, were computed according Eqs. (TO).
12p. Here, I denotes the loss function defined for an approximation Yn. at time point fn and the local error estimate err,- 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 generalize to new examples of the same class and new, unseen classes. Specifically, since the provided performance measure I is differentiable a.c.,
optimize the learning objective Eq. ( 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).
de
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,
5

---

Table 2: 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- Under review as a conference paper at ICLR 2020
Table I: On the test set of class (low-Ireq.. MLRK was considerably faster (less many steps) than the baseline controller on average, while causing only 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--reep.
low-freq
med-freg
interval
steps
error
steps
emfor
Baseline
MLRK
Baseline
MLRK
Baseline
MLRK
Baseline
MLRK
3.42
3.15
0.000004
0.000009
26.24
7.95
0.001588
0,008253
3
7.59
6.05
0.000017
0.000119
87.16
29.63
0.001686
0.01 1236
11.76
8.22
0.000032
0.000366
148.05
53.44
0.001739
0.012904
15.80
10.23
0.000048
0.000668
208.95
77.54
0.001735
0.014222
10
21.92
13.15
0.000073
0.001
300.32
113.82
0.001816
0.016546
high-froq
interval
steps
error
Bascline
MLRK
Baseline
MLRK
47, 15
12.08
0.026415
0.085082
157.58
53.42
0.023223
0.081219
268.03
96.48
0,025230
0.091109
378.42
139.69
0.026177
0.094129
crement in the global error
10
544.05
204.57
0.024858
0.094562
instances of class (high-freq).
investigated our suggested loss functions and our theory that suggests that providing higher order partial derivatives of g can be used as means 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 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 Lagrannnga loss Eq, 55) as well as the number of steps needed for an integration interval and the corresponding error separately. Since the Lagrangian represents 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 fow-freq). The solid lines show the mean loss in step n. the colored areas mark the standard deviation. MLRK achieves smaller mean loss than the baseline, which confirms its ability to generalize to new problem instances 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, []] the necessity for a trade-off between the two objectives becomes clear. While the baseline needs on average more steps complete the integration, it achieves smaller error, the contrary is 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 MIRK maintains 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 these problems. Tab. reveals that MLRK uses 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 MIRK generalizes to problem

---

Under review as a conference paper at ICLR 2020
MIRK
Baseline
MLRK
Baseline
MLRK
Baseline
0,0040
0.040
0.6
0.0035
0.035
0.5
0.0030
0.030
0.4
0.0025
0.3
0.0020
0.020
0.2
0.0015
0.015
0.1
0.0010
0.010
0.0
0.0005
0.005
0.1
0.0000
20
40
60
100
20
40
60
80
100
20
40
80
100
(a) (low-treq)
(b) (mod-froq
(c) (high-froo
Figure The mean L1 -Lgarannninnn loss ( of the approximation (tn, yn) 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) 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-freg MLRK shows mean loss similar to that of the baseline controller.
Baseline
partial
grad
0.003500
local error
0.002167
0.000833
0.000500
10
40
70
100
10
40
70
100
10
40
70
100
70
100
Figure 3: The mean local error over the test set of van der Pol ODEs 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). Providing local information such as those in input (8) and results in increased responsiveness of the optimizer.
instances of this class as well. In Tab. 2 we observe 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 MIRK 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 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 Li 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 following dataset of van der Pol oscillators of class vd 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 being damped. As the behavior of 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 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-Kutte step, Eq. (9).
der Pol oscillators. We considered the
The results are summarized in Fig. Using err shows reoccurring high spikes similar to the baseline, although the spikes of the optimizer are lower 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 smooth step size sequence
7

---

Under review as a conference paper at ICLR 2020
Baseline
partial
grad
0.002000
local ernrr
0.0013333
0.000667
0.0000001
70
100
70
100 40 100 10
100
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.
similar to the optimizer partial. This shows that the optimizer with additional input-informant from are actually outperforming the baseline. The reoecurring spikes in the local error seem to decrease in amplitude with time. This is largely due to 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, cff. 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 MIRK avoid these spikes by adjusting the step size appropriately, we consider the L1 local loss as in the previous experiment and a dataset of double pendulums of class (pendulum). The rraining 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. and reveal the baseline controllers inability to produce steady local errors. Here, the benefit of learned update rule over static one becomes evident. The optimizer err that is only provided with the local error estimates 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.
CONCLUSION
We have shown how to cast the design Runge- Kutta methods for solving initial value problems as learning problem. We established appropriate performance measures and useful inputs for the controller, which constitute the key ingredients of the resulting Meta-learning Runge- (MLRK), which learns step size controllers that are specialized Our experimental results demonstrate that MLRK can indeed learn to design novel Runge-Kuree Kutta methods that perform better than hand-designed Runge- approach. Furthermore, 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.
particular class of initial value problems.
observed
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 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 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 single large error can in general not be compensated for, even subsequent step sizes are chosen very small. this case, step is usually rejected and repeated using smaller step size. Here the challenge is 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.
8

---

Under review as a conference paper at ICLR 2020
REFERENCES
Marcin Andrychowicz, Masha Denil, Sergio Gómez Colmenarejo, Matthew W. Hoffman, David Pfau, Tom Schaul, Brendan Shillingford, and Nando de Freitas. Learning to learn by gradient descent by gradient descent. In D. D. Lee, M. Sugiyama, U. V. Luxburg, I. Guyon, and R. Garnett (eds.), Advances in Neural Information Processing Systems, volume 29, pp. 3981-3989, Barcelona, Spain, December 2016.
Fischer Black and Myron Scholes. The pricing of options and corporate liabilities. Journal of Political Economy, 81(3):637-65 May-Jum 1973.
John C. Butcher. Numerical Methods for Ordinary Differential Equations. Wiley, Auckland, New Zealand, 2nd edition, 2008
Ricky T. Q. Chen, Yulia Rubanova, Jesse Bettencourt, and David K. Duvenaud. Neural ordinary differential equations. In S. Bengio, H. Wallach, H. Larochelle, K. Grauman, N. Cesa- Desa-Bianchi, and R. Garnett (eds.), Advances Neural Information Processing Systems, volume 31, pp. Montreal, Canada, December 2018.
Yutian Chen, Matthew W. Hoffman, Sergio Gómez Colmenarejo, Misha Denil, Timothy P. Lillicrap, Matt Botvinick, and Nando de Freitas. Learning learn without gradient descent by gradient descent. Proceedings of the 34th International Conference on Machine Learning, pp. 748-756, Sydney, Australia, August 2017.
Peter Deufthard. Differential equations in technology and medicine. Computational concepts, adaptive algorithms, and virtual labs. In V. Capasso, H.W. Engl, and J. Periaux (eds.), Computational Mathematics Driven by Industrial Applications, volume 1739 of Lecture Notes in Mathematics, pp. 69-125. Springer, Berlin, Heidelberg, Germany, 2000.
Weinan E and Bing Yu. The deep Ritz method: A deep learning-base numerical algorithm for solving variational problems. Communications in Mathematics and Statistics, 6(1):1-12, March 2018.
Mojgan Esna-Ashari, Maryam Zekri, Masood Askari, and Noushin Khalili. Predictive control of the blood glucose level in type diabetic patient using delay differential equation Wang model. Journal of Medical Signals and Sensors, 7(1):8 (-0,;820,.000 January-
2017.
Kjell Gustafsson.
techniques for stepsize selection in implicit Runge- Kutta methods. ACM Transactions on Mathematical Software, 20(4):496-5) 517, December 1994.
Kjell Gustafsson, Michael Lundh, and Gustaf Söderlind. API stepsize control for the numerical solution of ordinary differential equations. BIT Numerical Mathematics, 28(2):270-25) 287, June 1988.
Ernst Hairer and Gerhard Wanner. Solving Ordinary Differential Equations Stiff and Differential- Algebraic Problems, volume 14 of Springer Series in Computational Mathematics. Springer, Berlin, Heidelberg, Germany, 2nd edition, 2002.
Erst Hairer, Syvert P. Nørsett, and Gerhard Wanner. Solving Ordinary Differential Equations I Nonstiff Problems, volume Springer Series Computational Mathematics. Springer, Berlin, Heidelberg, Germany, 2nd edition. 2000.
Jiequn Han, Arnulf Jentzen, and Weinan E, Solving high-dimensional partial differential equations using deep learning. Proceedings of the National Academy of Sciences, (15(34):8505-8510 August 2018
Kurt Homik. Approximation capabilities of multilayer feedforward networks. Neural Networks, 4(2): 251-257, 1991.
Mihai Ilea, Marius Turnea, and Mariana D. Rotariu. Differential equations with applications in cancer diseases. Medical-Surgical Journal of the Society of Physicians and Naturalists lasi, 117 (2):572-577, April-June 2013.

---

Under review as a conference paper at ICLR 2020
Isaac E. Lagaris, Aristidis Likas, and Dimitrios I. Fotiadis. Artificial neural networks for solving ordinary and partial differential equations. IEEE Transactions on Neural Networks, 95)0987-1000, September 1998.
Alfred J. Lotka. Elements of physical biology. Nature, 116(2917), September 1925.
Qazi Ibadur Rahman and Gerhard Schmeisser. Analytic Theory of Polynomials. Number 26 in London Mathematical Society Monographs. Oxford University Press, Oxford, UK, 2002.
David E. Rumelhart, Geoffrey E. Hinton, and Ronald J. Williams. Learning representations by back-propagating errors. Nature, 323(6088):333-536 1986.
Gudrun Scholz and Fritz Scholz. First-order differential equations in chemistry. ChemTexts, 1(1), November 2014.
Patrick Schramowski, Christian Bauckhage, and Kristian Kersting. Neural conditional gradients. 1803.0400300000 March 2018.
Justin Sirignano and Konstantinos Spiliopoulos. DGM: A deep learning algorithm for solving partial differential equations. Journal of Computational Physics, 375: 1339-1444, December 201
Gustaf Söderlind. Automatic control and adaptive time-stepping. Numerical Algorithms, 31(1): 281-310, December 2002.
Gustaf Söderlind. Time-step selection algorithms: Adaptivity, control, and signal processing, Applied Numerical Mathematics,
March-April 2006.
Vito Volterra. Variazioni fluttuazioni del numero d' d'individui in specie animali conviventi. Memoria della Reale Accademia Nazionale dei Lincei, 2:31-113, 1926.
Olga Wichrowska, Niru Maheswaranathan, Matthew W. Hoffman, Sergio Gómez Colmenarejo, Misha Denil, Nando de Freitas, and Jascha Sohl-Dickstein. Learned optimizers that scale and generalize. In Proceedings of the 34th International Conference on Machine Learning, pp. 3751-3760, Sydney, Australia, August 2017.
A APPENDIX
RUNGE- KUTTA METHODS
The mathematical theory Runge- Kutta-Kutta methods is quite evolved. We will only recap few basics here, for further details see for example (Butcher 2008 Hairer et al. 2000 Hairer & Wanner 2002).
Definition A.1. Let E N, aij
E IR for i,j (1,..., s). E R for i E (1,...,8) and b; R for Let g IR Rm > T function that describes an initial value problem. The method defined by
i E (1.....8). $).
Un+1 Un h
beke.
(14)
k; = gtn + hc,, yo + h 2aighi)
(15)
for step size h is called an s-stage Runge-Kuuue Kutta method. If asj = 0 for i S j. the method called explicir, otherwise it is called implicit.
Definition A.2. A Runge-Kutte Kutta method is of order if for sufficiently smooth initial value problems
(tv + h)) vill < KhP+1
holds.
10

---

Under review as a conference paper at ICLR 2020
Note that higher order method yields more accurate results, consequently high order is desired characteristic of a Runge-Kutte Kutta method. However, higher orders can only be achieved by the use of more stages. For an s-stage Rumge-Kuue Kutta method the order is bounded by the number of stages, S 8. Up to order 4 there exist methods with 8, for order and higher is strictly smaller than c.f- (Butcher 2008 Theorem 324B). The order of s-stage Runge- Kutta method depends on the coefficients Gij. bi, C.. Order conditions for these coefficients have been developed, which ensure 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. En 3n etr, and possibly additional characteristics g. we must choose a step size h, which is then used to evaluate g at different points. These points are determined by h22 see Definition A.1 The resulting stages in Eq. (15) are then used to compute the approximation yn+1 of y( + hm)) c.f' Eq. (14). The problem is to choose h2 in way such that the approximation Un+1 fulfills some desired properties. For example we may require that the local error of Vn+1 is close some tolerance parameter. In fact, this is the standard objective of most common step size controllers.
ERROR ESTIMATION
The local emrors 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 (f,,.) of order á and ğ respectively, that is
Ün s(.,) + O(hBt1),
In = y(nn) + O(h8tl),
where v()) is the solution of
y' v((1) s((,v((()))))))))))))))(())))s v(hn-1) yn-1,
and the step size h that was used to obtain both approximations is given by
h t fn-1-
If >
In Dn = y((n) - In + O(h+22)
can be used as an approximation for the local truncation error of the p-th order Runge-Kutta method:
enr, (h) ll... Gmll
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-K Kutta methods compute ÿ and ÿ. however this is quite costly, especially for implicit methods. A more efficient approach to obtain error estimates are the embedded Kunge-Kutta methods. Here both approximation use the same stages k; in Eq. (15) but two different pairs of coefficients b; by in Eq. (14). By doing so, the stages can be reused and an error estimate can be obtained with cheap extra linear combination of the k; 's. We denote the error estimates by
e 1l In ll.
(16)
CLASSES OF INITIAL VALUE PROBLEMS
For loss functions such as the one in Eq, (5) 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 of interest. For these loss functions it is interesting to consider mitial value problems for which the local emrors vary drastically during the integration such van der Pol oscillators.
11

---

Under review as a conference paper at ICLR 2020
SIMPLE HARMONIC OSCILLATOR
simple harmonic oscillator is a harmonic oscillator that is neither driven nor damped and can be characterized by
mx"(t) (t) =k
where m is the mass of the oscillator, x the position and describes the restoring force that is applied, when displaced from its equilibrium position. It can be transformed to
3 (t) y2(t):
(1)
m
k: v1())
The initial values provide the initial position x=fo) of x and the initial velocity z'(t)),
8/1 (to) x(ta);
y2((o) (to).
The solution of the simple harmonic oscillator takes the form
m
<(t) A cos + (p),
K
where and o are uniquely determined by the initial values. The frequency of the oscillations is determined by w, A is the magnitude of the oscillations and yo is phase-shift.
class of this kind of initial value problem can described by distribution over m and the initial values alternatively distribution over A, w and sp The following classes of initial value problems were used in the experiments
A (0, 5),
(0, 10), i (0,((0 1), ~
(high-freq)
A ~ (0, 5),
P~~ 1),
(higher-freq)
A ~ (O, 1),
(0,1), 8-~ U0,11) 1),
(low-freq) (med-freq)
~~ U115), 8~ (0,, 1), where (a, b) denotes the uniform distribution over the interval (a, b).
A ~ U11 5),
LINEAR CONSTANT COBFFICIENT
A linear differential equation with constant coefficients can be described by
y (1) v(), y(4n) = Vo
where A E mm The solution of this differential equation is known and can be expressed in terms of the eigenvalues .g...., .a.....am Am with corresponding eigenvectors U1 ......, Em of the matrix A:
y(t)
Σ j=1
Cj Axt Uy:
with coefficients Cj E R that are uniquely determined by the initial values y(to) 9o- The initial value problem is stable, if Re(A,) 0 for all j E (1....,m) and unstable if (e()) > O for some j {1....,m)) m). The oscillatory behavior of the problem is determined by Im(A,). Note that
the simple harmonic oscillators are linear differential equations with constant coefficients, where 0 for j {1,2) and Im(A1) Im())) / i.
A ? Re(A,)
For the experiments used differential equations with A 3x3 The matrix A either has three real eigenvalues λ1, >2, Ag E R or one real eigenvalue X E IR and two complex eigenvalues bi, with a, b E R. In this case we choose 02 = C3 to obtain real solution.
d + bi, A3 We can specify a class of linear constant coefficient differential equations by distribution on the eigenvalues of A and the coefficients cj. The classes that were used for the experiments were described by the distributions
A1 ~ 2-5.5, -0.1),
a~ U(-20(15 -0.1), b~ U(0,5),
(stable)
A1 1),1),
~~ (25,, -0.1,, a (0,
b b~ (0,5),,
(unstable)
Åι ~ U(-1,0),
b~ (1,2)),
(osc-increasing)
where U(a b) denotes the uniform distribution on the interval [a. b).
12

---

Under review as a conference paper at ICLR 2020
v(t)
-20
15
25
30
0.0045
errit
0.0030
0.0015
0.0000
10
15
20
25
30
Figure The upper plot shows the numerically approximated solution yı (t) of 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.
VAN DER POL OSCILLATOR
Van der Pol oscillators arise in the study of limit cycles. A van der Pol oscillator can be described by
1/1 (t) y2
1/ (0) 2,
16(() = o11 -f)0z V1-
3/2 (0) 0,
with σ 0. For 2 an approximate solution is displayed in Figure Small oscillations are
amplifies and large oscillations damped, c.f. (Hairer et al., 20000, 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
o ~ U0,
(vdP 1)) (vdP(1,
U1122.
DOUBLE PENDULUM
double pendulum is pendulum attached to another pendulum, c.f. Fig. 6 The length the first and second pendulum given by L, and L2, respectively, their masses by Mj and M2 and the angle by 0 and 0g. The equations of motion for the double pendulum can be written in different forms, chose the following formulation in terms of the angular acceleration:
0f sin(0:
02) cos (2 01) + sin cos(62 01)
01
L, (M, + M2 M2 cos? 11))
M2L203 sin (02 01) - (M1 + singgsin(())
(M M2 M2 cos" (02 01))
(Mi M2) (s))n(0)) cos(02 01) L, 0f sin(02 - 01) - ginn((2))
82
L2(M] 1 + M2 M2 2 01))
M2L sin 02 02) cos (22 01)
L2 + M2 M2 cos (02 01))
The initial angular position and velocity 0,, 02 0, 02 determine the trajectory of the double pendulum. For the experiments we used double pendulums with the following distribution
M₁ U0.55 1), M2 U(0.5, 2), Ln U0..5 1), L2 U0.5, 1).
(pendulum)
ADDITIONAL EXPERIMENTS
ON VAN DER PoL.
Table shows the mean number of steps, local error and wall clock time over 1500 van der Pol equations of class (v
13

---

Under review as a conference paper at ICLR 2020
01
L1
M1
L,2
02
M2
Figure 6: A double pendulum is one pendulum attached to another pendulum. The masses M1 and M2 of these pendulums may vary as do the lengths L, and Lz. The initial configuration 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. ^ Thus we chose the following datasets:
Training set: 30000 Harmonic oscillators of class (high-freq.
Validation set:
1000 Harmonic oscillators of class (high-
Test set: 1000 Harmonic oscillators of class (high-freqq.
As loss we chose the La-Lagrangann Eq. (5). As inputs we used only the error estimates, Eq. 77) The loss of the learned controller and the baseline controller, averaged over the test dataset, is shown in Fig. ((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 (In, yv). 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 the baseline controller, which demonstrates the controllers generalization capability to new problems of the same class.
Tab. 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 to 10. All values are averaged over the test dataset. Here, the baseline controller needs much more steps than the learned controller, but maintains smaller error. In fact, the етror 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 longer integration interval.
While the general assumption that the global error increases with the length of the integration interval often fulfilled, the specific assumption of a linear accumulation of the global error, in our loss, usually not fulfilled. This is a good example for why the assumptions leading to the widely used strong.
objective to keep the local errors constant
14

---

Under review as a conference paper at ICLR 2020
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 Figure 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 enror.
int
steps
baseline
CI
partial
grad
21.59
16.42
12.40
12.09
3
33.74
29.12
25.16
24.74
5
45.43
41.07
36.42
36.02
7
56.84
52.57
48.34
47.97
10
73.46
69.41
65.40
65.04
int
mean local error
baseline
сT
partial
grad
7.17e 4
6.58e
4.0le
3.74e
3
6.40e 4
5.45e
2.49e
2.28e
5
5.e
4
4.47e 4
1.95e
1.75e 4
7
4.97e
4.16c
1.57e
1.39e 4
10
4.59c
4
3.82c
1.32e
1.18c
int
time
bascline
сT
partial
grad
0.0255
0.0263
0.0254
0.0221
3
0.0405
0.0375
0.0460
0.0403
5
0.0591
0.0517
0.0681
0.0596
7
0.0858
0.0742
0.1036
0.0897
10
0.0971
0.0825
0.1201
0.1065
interval
steps
error
Baseline
Optimizer
Baseline
Optimizer
1
45.50
13.73
0.023092
0.030499
3
152.11
43.12
0.022266
0.040429
5
258.73
74.40
0,022731
0.049529
365.30
106.31
0.022011
0.055821
10
525.18
154.25
0.023353
0.070451
Table 4: The mean values over the test set consisting of 1000 harmonic oscillators class (high-freq) shown. While the baseline needs much more steps than trained model, it achieves smaller global error. The global error of the baseline stays approximately the same for all integration intervals.
also interesting 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 test set consisting of 1500 harmonic oscillators of class (higher-freq). This class contains oscillators of class (high-freg) but also includes higher frequency ones. The learned controller has lower mean loss than the bascline controller which indicates that able to transfer knowledge from problems of class (nigh-freq) to problems of class (higher-freq).
When we compare the number of steps and the global error at the end the integration interval, see Tab. observe similar situation in Tab. 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.
15

---

Under review as a conference paper at ICLR 2020
LSTM
Baseline
LSTM
Baseline
0.5
12
10
0.4
0.3
0.2
loss
loss
0.1
0.0
0.1
0.2
20
40
60
80
100
20
40
60
80
100
(a) high-freq
(b) higher-freg
Figure 7: The red and blue line show the mean loss Eq. ( of the approximation (fn, yn) in step over the test set. The red and blue area indicate the standard deviation. [a] The learned model LSTM attains smaller mean loss on new instances of the same class of problems (high-freg) it was trained on, this especially pronounced in the steps n Z (b) The test 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 lower than the mean loss of the baseline controller which indicates that method generalizes to these different problem instances.
interval
steps
error
Baseline
Optimizer
Baseline
Optimizer
1
116.76
22.26
0.453698
0.511174
3
386.40
75.64
0.442521
0.620779
5
656.06
141.29
0.426773
0.770457
925.80
212.19
0.418467
0.925463
10
1330.31
319.85
0.413902
1.427778
Table 5: The mean number of steps and the mean global error at the end of the integration interval computed over 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 lower global error, which stays approximately the same for different interval lengths.
HARMONIC OSCILLATOR
We learn controller optimized for simple harmonic oscillators. We choose the following data sets:
Training set: 30000 Harmonic oscillators class flow-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 test set consisting of the same class 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, 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.
find that our method generalizes to these problem instances well. In fact,
16

---

Under review as a conference paper at ICLR 2020
LSTM
Baseline
LSTM
Baseline
tol
0.0025
0.0014
0.0013
0.0020
0.0015
error
local errar
0.0011
0.0010
0.0009
0.0005
0.0008
0.0000
0.0007
-0. 0.0000
20
30
40
50
60
70
80
100
0.0000
100
(a) (low-freq)
(b) (high-freq)
Figure 8: The mean local error of our learned controller and the baseline controller is shown. Testing on 1500 harmonic oscillators of class (low-freq) shows that both the baseline and the learned controller able 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 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 similar results.
LSTM
Baseline
tol
Baseline
0.0015
0.0035
0.0011
0.0013
0.0025
local error
0.0020
0.0012
0.0015
local
0.0011
0.0010
0.0005
0.0009
0.0000
0.000g
20
30
40
50
60
70
60
90
100
-0.0005
20
30
40
60
100
(a) (osc-increasing)
(b) (v,P 1)
Figure 9: The mean local error different test sets is shown. (a) The test set consists of 1500 linear differential equations with constant coefficients of class (osc-increasing). 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 of 1500 van der Pol oscillators of type (vdP(0, 1) reveals 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
17

---

Under review as a conference paper at ICLR 2020
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-increasinn). the performance our controller is similar to the baseline controller. The learned controller obtains local erross which are closer the tolerance in some steps. Figure [() shows 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 van der Pol oscillator changes from being driven to being damped as demonstrated in Figure Our controller does not react quickly to the sudden changes the local error of van der Pol equations. This behavior is very different from harmonic oscillators, thus generalization can not be expected.
Figure (a) Here,
18


Reviews:
[{"review_id": "Syl61bA3KB", "reviewer": "AnonReviewer1", "summary": "", "questions": "", "limitations": "", "rating": 3, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "Updated review: Thanks to the authors for their response to my comments. I believe the strong point of this paper is the novel idea, however, I find the justification for that idea incomplete as author's seems to suggest that the proposed method is probably computationally more expensive (which is opposite to the original motivation of the paper).\n\n---------------------------------------------------------------------------------------------------------------------------------------------------------------------\n\nSummary: This paper casts the problem of step-size tuning in the Runge-Kutta method as a meta learning problem. The paper gives a review of the existing approaches to step size control in RK method. Deriving knowledge from these approaches the paper reasons about appropriate features and loss functions to use in the meta learning update. The paper shows that the proposed approach is able to generalize sufficiently enough to obtain better performance than a baseline. \n\nI think this paper, in general, is clear and well-written. I believe the idea of the paper is interesting too. \n\nThe paper argues that the main challenge of solving the step size control problem for the RK method is balancing the computation vs accuracy trade-off. Existing methods tackle this problem in different ways and this paper proposes to solve it via meta-learning. However, the paper does not mention how and why meta-learning is expected to tackle this challenge?\nSo a couple of comments on what set of problems do we expect meta-learning to better tackle this trade-off than the existing methods would have been useful. I am wondering if it is even possible to say something about this in principle? \n\nThe paper argues that the idea behind using meta-learning is to learn behaviour from a given class of problems and then generalize to new unseen problems (from the same or different classes). \nHow do we know that these problems are even from same distributions? \nWon't the proposed approach fail spectacularly when the problems are not from the same distribution? It would have been nice if the paper made this distinction even if empirically. \n\nIn the experiments section, I could not find/understand what exactly is the baseline the paper is comparing to. \n\nI was more interested in a study that compared the performance of MLRK as the number of instances of the training problems are varied. \nThis again makes me come back to the original point of computational cost vs accuracy. What is the computational cost of collecting data on 30000 instances of problems? Should we not worry about this cost?\nAlso, what is the computational cost of the proposed approach and why are we not comparing it to existing approaches/baseline?\n\nminor comments: \nwhat is tol? it tol the same tolerance as lambda.\n \n", "extra_scores": {"experience_assessment": "I do not know much about this area.", "review_assessment:_checking_correctness_of_experiments": "I assessed the sensibility of the experiments.", "review_assessment:_thoroughness_in_paper_reading": "I read the paper at least twice and used my best judgement in assessing the paper.", "review_assessment:_checking_correctness_of_derivations_and_theory": "I assessed the sensibility of the derivations and theory."}, "extra_text": {}}, {"review_id": "Hke_qEKjtB", "reviewer": "AnonReviewer2", "summary": "", "questions": "", "limitations": "", "rating": 3, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "## Summary ##\n\nThe authors present a method for learning a step-size adaptation strategy for Runge-Kutta methods. The principal contributions of the paper are:\n\n1. They define a loss function that better captures global performance of the controller, rather than just local behavior.\n2. 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 $g$.\n3. They describe a recurrent architecture for adapting step size over time.\n\nRunge-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.\n\nWhile 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.\n\n## Specific Comments ##\n\n1. 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,\n\n   * Inputs: t, y(t), g  # Also possibly intermediate values from previous time steps.\n   * Select a step size h(t)\n   * Evaluate g at p different points based on h(t).\n   * Use these evaluations to compute a value of y_(t + h_t)\n\n   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.\n\n2. 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.\n\n3.  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?\n\n4. There are a number of issues with the experiments that I think could use clarification or improvement:\n\n  (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?\n\n  (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?\n\n  (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.\n\n  (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.\n\n  (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.", "extra_scores": {"experience_assessment": "I have read many papers in this area.", "review_assessment:_checking_correctness_of_derivations_and_theory": "N/A", "review_assessment:_checking_correctness_of_experiments": "I assessed the sensibility of the experiments.", "review_assessment:_thoroughness_in_paper_reading": "I read the paper at least twice and used my best judgement in assessing the paper."}, "extra_text": {}}, {"review_id": "BJeZoWc19r", "reviewer": "AnonReviewer3", "summary": "", "questions": "", "limitations": "", "rating": 3, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "The paper proposes to learn the step size for a Runge-Kutta numerical integrator for solving ordinary differential equations initial value problems. The authors frame the stepsize control problem as a learning problem, based on different performance measures, on ODE dependent inputs and on a LSTM for predicting the next step coefficient. Experiments are performed on 3 ODEs by training and testing in different contexts.\nThe problem of designing adaptive controllers for ODE numerical schemes is interesting and is probably a new application for ML. The paper makes an effort to introduce the necessary background and for reviewing some classical adaptive controller techniques. The description of the method is relatively clear, but could however be largely improved in order to make it more accessible to the audience.  Many of the arguments and proposed ideas come without justification, some definitions should be made more precise. The construction of the training and test sets should be better explained. The experiments show that the proposed approach leads to fewer evaluations but larger mean errors.  The graphics also show that the local error is smaller for the proposed method than for the baselines which is in contradiction with the global error behavior. This should be clarified – the relations between the two error types should be made clear.  The baseline is not defined in the text so that it is difficult to judge the performance. Why not comparing to several adaptive baselines?\n\nIn conclusion, this is an interesting topic, the paper proposes new ideas. A more careful writing and especially a better comparison with sota baselines would greatly improve the paper. \n\n\n------ post rebuttal -----\nThanks for the answers. I still think that the ideas are interesting but that the experiments  do not demonstrate enough of the proposed method. I will keep my score.\n", "extra_scores": {"experience_assessment": "I have read many papers in this area.", "review_assessment:_checking_correctness_of_experiments": "I assessed the sensibility of the experiments.", "review_assessment:_thoroughness_in_paper_reading": "I read the paper at least twice and used my best judgement in assessing the paper.", "review_assessment:_checking_correctness_of_derivations_and_theory": "I assessed the sensibility of the derivations and theory."}, "extra_text": {}}]

Author rebuttal:

Meta-review: Summary: This paper casts the problem of step-size tuning in the Runge-Kutta method as a meta learning problem. The paper gives a review of the existing approaches to step size control in RK method. Deriving knowledge from these approaches the paper reasons about appropriate features and loss functions to use in the meta learning update. The paper shows that the proposed approach is able to generalize sufficiently enough to obtain better performance than a baseline. 


The paper was lacking in advocates for its merits, and needs better comparisons with other baselines before it is ready to be published.
Decision: Reject