# Discretize-Optimize vs. Optimize-Discretize for Time-Series Regression and
Continuous Normalizing Flows

Derek Onken<EMAIL_ADDRESS>
Lars Ruthotto<EMAIL_ADDRESS>
Department of Computer Science
Department of Mathematics
Emory University
Atlanta, GA 30322, USA

###### Abstract

We compare the discretize-optimize (Disc-Opt) and optimize-discretize (Opt-
Disc) approaches for time-series regression and continuous normalizing flows
(CNFs) using neural ODEs. Neural ODEs are ordinary differential equations
(ODEs) with neural network components. Training a neural ODE is an optimal
control problem where the weights are the controls and the hidden features are
the states. Every training iteration involves solving an ODE forward and
another backward in time, which can require large amounts of computation,
time, and memory. Comparing the Opt-Disc and Disc-Opt approaches in image
classification tasks, Gholami et al. (2019) suggest that Disc-Opt is
preferable due to the guaranteed accuracy of gradients. In this paper, we
extend the comparison to neural ODEs for time-series regression and CNFs.
Unlike in classification, meaningful models in these tasks must also satisfy
additional requirements beyond accurate final-time output, e.g., the
invertibility of the CNF. Through our numerical experiments, we demonstrate
that with careful numerical treatment, Disc-Opt methods can achieve similar
performance as Opt-Disc at inference with drastically reduced training costs.
Disc-Opt reduced costs in six out of seven separate problems with training
time reduction ranging from 39% to 97%, and in one case, Disc-Opt reduced
training from nine days to less than one day.

Keywords: neural ODEs, optimal control, Optimize-Discretize, Discretize-
Optimize, normalizing flows

##  1 Introduction

Coined by Chen et al. (2018), neural ODEs are ordinary differential equations
(ODEs) with neural network components (Section 2). Neural ODEs were developed
as continuous limits of discrete Residual Networks (ResNets) (He et al.,
2016), which can be seen as forward Euler discretizations of a continuous
time-dependent ODE (E, 2017; Haber and Ruthotto, 2017). Since ResNets
demonstrate impressive performance in applications ranging across image
classification, segmentation (Zhang et al., 2017), and deblurring (Nah et al.,
2017), this result motivated several follow-up works. Examples include
continuous models for PolyNet, FractalNet, and RevNet (Lu et al., 2018), the
use of higher-order ODE solvers (Zhu et al., 2018), architectures for
convolutional ResNets motivated from partial differential equations (PDEs)
(Ruthotto and Haber, 2019), and continuous versions of normalizing flows
(Rezende and Mohamed, 2015; Chen et al., 2018). The last work shows neural
ODEs’ potential for scientific machine learning applications that involve
dynamical systems, which differ from classification and other data science
problems traditionally solved with ResNets. This interplay between physics,
differential equations, and neural networks is apparent in Raissi (2018),
Brunton et al. (2019), Köhler et al. (2019), and Rackauckas et al. (2020). We
note earlier versions of continuous neural networks reminiscent of neural ODEs
can be found in Rico-Martínez et al. (1992) and González-García et al. (1998),
which also show the benefits of continuous models for time sequence data and
solutions of partial differential equations.

Training neural ODEs consists of minimizing a regularized loss over the
network weights subject to the nonlinear ODE constraint. Therefore, training
can be seen as an optimal control problem. Applying concepts from optimal
control theory has emerged as an active research area and more insight has
been gained in recent years. For example, Pontryagin’s maximum principle has
been used to efficiently train networks with discrete weights (Li et al.,
2017), multigrid methods have been proposed to parallelize forward propagation
during training (e.g., Günther et al., 2019), and analyzing the convergence on
the continuous and discrete level has led to novel architectures (Benning et
al., 2019).

In the language of optimal control, the adaptive adjoint-based approaches in
Chen et al. (2018), Rackauckas et al. (2019), and Grathwohl et al. (2019) can
be viewed as Optimize-Discretize (Opt-Disc) approaches since they optimize the
continuous ODE and discretize the optimal dynamics after training. A well-
known alternative is the Discretize-Optimize (Disc-Opt) approach, where one
first discretizes the continuous problem and then solves a finite-dimensional
optimization problem (Gunzburger, 2003). Gholami et al. (2019) give a thorough
discussion of the trade-off between Opt-Disc and Disc-Opt as well as numerical
examples for neural ODEs in image classification. Our goal in this paper is to
extend this discussion and perform similar experiments for time-series
regression and continuous normalizing flows using neural ODEs.

Common challenges in both Opt-Disc and Disc-Opt are the computational costs
associated with solving the ODE during forward propagation and, when using
gradient-based optimization, the calculation of the gradient via
backpropagation. Accurately solving the forward propagation requires vast
quantities of memory and floating point operations, which can make training
prohibitively expensive. As data in scientific machine learning applications
is often noisy, using a low-accuracy solver to speed up computations appears
tempting. However, when the forward propagation and adjoint are not solved
accurately, the quality of the gradient can deteriorate, an important drawback
of the Opt-Disc approach (Gholami et al., 2019). In Disc-Opt, the accuracy of
gradients does not depend on the accuracy of the forward propagation. This
crucial advantage of Disc-Opt motivates its use to improve the efficiency of
training; Gholami et al. (2019) provide detailed analysis and examples on
image classification tasks. Disc-Opt approaches include the original
formulations of ResNets (He et al., 2016) as well as the first continuous
ResNets (E, 2017; Haber and Ruthotto, 2017). Finally, Disc-Opt approaches are
easy to implement, especially when using automatic differentiation (Nocedal
and Wright, 2006), while Opt-Disc methods require the numerical solution of
the adjoint equation (Bliss, 1919).

Both approaches require storing or recomputing intermediate hidden features
from the forward propagation to differentiate the loss function wit respect to
the neural network’s weights during gradient-based optimization. While Chen et
al. (2018) propose to avoid these memory costs by recomputing the neural ODE
backward in time during the adjoint solve, time reversal of ODEs is known to
be prone to numerical instabilities (see example in Gholami et al., 2019).
Avoiding the storage altogether is possible when using forward-backward stable
networks, e.g., the Hamiltonian network (Chang et al., 2018; Ruthotto and
Haber, 2019), but not necessarily for any architecture. For general network
structures, checkpointing schemes can be used to reduce the memory footprint
without loss of stability (Gholami et al., 2019).

In this paper, we compare the Opt-Disc and Disc-Opt neural ODE applications in
time-series regression and continuous normalizing flows (CNFs). Here,
meaningful models must also satisfy additional requirements, e.g., the
invertibility of the CNF. As the continuous model satisfies these requirements
by design, Opt-Disc approaches may appear advantageous. However, due to the
prevalence, importance, and interest in neural ODEs for scientific
applications, we believe that a comparison of Opt-Disc and Disc-Opt comparison
is warranted. Through our numerical experiments, we demonstrate that with
careful numerical treatment, Disc-Opt methods can achieve competitive
performance to Opt-Disc at inference with drastically reduced training
costs.111Open-source code is provided at
https://github.com/EmoryMLIP/DOvsOD_NeuralODEs. To this end, we decouple the
weights and layers, which allows us to train the neural ODE with a reduced
number of layers and satisfy the additional requirements through re-
discretization. For this approach to be effective, we highlight the importance
of using sufficiently small step sizes in the Disc-Opt approach. We show that
a too coarse discretization can compromise the invertibility. Therefore, our
results align with and extend the image classification results achieved in
Gholami et al. (2019) to these learning problems using dynamical systems.

Figure 1: Illustration of the Disc-Opt approach for the Gaussian mixture
problem in Section 4.2.1. As common in optimal control, we use different
discretization points for the weights and the states. Following the
terminology in optimal control, we refer to the time points as control and
state layers, respectively. We discretize the state equation using a Runge-
Kutta 4 scheme with constant step size $h=1$. Image is not to scale.

##  2 Background

We define neural ODEs, place them in the context of discrete neural networks,
and provide a brief groundwork of optimal control theory for the Disc-Opt and
Opt-Disc comparison.

###  2.1 Neural ODEs

Neural ODEs, coined by Chen et al. (2018), are differential equations with
neural network components. The neural ODE performs a nonlinear transformation
of the feature space. We denote the transformation of a feature vector
$\mathbf{y}_{0}\in\mathbb{R}^{n_{f}}$ by model
$f(\boldsymbol{\theta},\mathbf{y}_{0})=\mathbf{y}(T)$ where $T$ is some
arbitrary final time and $\mathbf{y}\colon[0,T]\to\mathbb{R}^{n_{f}}$
satisfies the initial value problem

$\displaystyle\partial_{t}\mathbf{y}(t)$
$\displaystyle=\ell(\boldsymbol{\theta}(t),\mathbf{y}(t),t)\,,\quad\text{for}\quad
t\in(0,T]$ (1) $\displaystyle\mathbf{y}(0)$ $\displaystyle=\mathbf{y}_{0}.$

Here, $\boldsymbol{\theta}\colon[0,T]\rightarrow\mathbb{R}^{n_{p}}$ are the
weights of the neural ODE, and
$\ell\colon\mathbb{R}^{n_{p}}\times\mathbb{R}^{n_{f}}\times[0,T]\to\mathbb{R}^{n_{f}}$
is a neural network. We formulate layer $\ell$ to accept time $t$ as a
parameter though only in rare cases (Section 4.2) does a neural network layer
operate on time itself.

A prominent class of neural ODEs are continuous limits of ResNets (He et al.,
2016). An $N$-layer ResNet can be written as

$\mathbf{y}_{j+1}=\mathbf{y}_{j}+h\,\ell\left(\boldsymbol{\theta}_{j},\mathbf{y}_{j},t_{j}\right),\quad\text{where}\quad
j=0,1,\ldots,N{-}1,$ (2)

with some function $\ell$ and step size $h=T/N$. We refer to the
$\boldsymbol{\theta}_{j}$ as control layers and the $\mathbf{y}_{j}$ as state
layers (Figure 1). Interpreting the weights $\boldsymbol{\theta}_{j}$ as
evaluations of $\boldsymbol{\theta}$ at the time points $t_{j}=j\cdot h$ and
taking the limit as $N\to\infty$, we see that the $N$-layer ResNet (2)
converges to the neural ODE (1). Infinitely deep networks such as this have
outperformed shallower networks in recent years (Sonoda and Murata, 2019).
Analogously, one can view the $N$-layer ResNet as a forward Euler
discretization of the neural ODE (E, 2017; Haber and Ruthotto, 2017). This
observation led to follow-up works on analyzing and improving ResNets via
insights from the continuous model. Since ResNets are discrete processes and
ultimately a discrete network is desired, some caution must be used when
studying their continuous limits (Ascher, 2020).

In this paper, we consider a different class of neural ODEs arising in time-
series regression and CNFs. Unlike ResNets, this class of learning problems is
phrased in the continuous setting, and training aims at obtaining the
continuous neural ODE model for prediction or inference. These classes also
impose additional requirements on meaningful models, e.g., the ability to
reverse the neural ODE in time. While the continuous neural ODE satisfies
these requirements by design, discrete versions may not.

###  2.2 Learning and Optimal Control

Training the neural ODE (1) can be phrased as an infinite-dimensional optimal
control problem

$\displaystyle\min_{\boldsymbol{\theta},\mathbf{y}}$
$\displaystyle\left\\{\mathcal{J}(\boldsymbol{\theta},\mathbf{y}):=\frac{1}{S}\sum_{i=1}^{S}\mathcal{L}(\boldsymbol{\theta},\mathbf{y}^{(i)})+\mathcal{R}(\boldsymbol{\theta})\right\\}$
(3) subject to $\displaystyle\mathbf{y}^{(i)}\text{ solves
}\text{\eqref{eq:cont_ode}}\text{ with initial value }\mathbf{y}_{0}^{(i)}$

where $\mathbf{y}^{(1)},\ldots,\mathbf{y}^{(S)}$ satisfy the neural ODEs for
initial values given by the training data; loss functional $\mathcal{L}$ and
regularization functional $\mathcal{R}$ are chosen to model a given learning
task. For regularization, we consider Tikhonov regularization (often called
weight decay, Krogh and Hertz, 1992)

$\mathcal{R}(\boldsymbol{\theta})=\frac{\alpha}{2}\int_{0}^{T}\|\boldsymbol{\theta}(t)\|^{2}\,\mathrm{d}t,$
(4)

where we assume that the regularization parameter, $\alpha>0$, is chosen and
kept fixed.

We will consider two types of loss functionals; regression loss and
likelihood. The former is used when training a neural ODE to approximate a
given function $\mathbf{u}\colon[0,T]\to\mathbb{R}^{n_{f}}$ and reads

$\mathcal{L}(\boldsymbol{\theta},\mathbf{y})=\int_{0}^{T}L\big{(}\mathbf{y}(t),\mathbf{u}(t)\big{)}\,\mathrm{d}t,\quad
L(\mathbf{y},\mathbf{u})=\frac{1}{2}\|\mathbf{y}-\mathbf{u}\|^{2}.$ (5)

When the $\mathbf{u}$ is known only at some time points, a discretized version
of this functional $\mathcal{L}$ is used (e.g., the time-series example in
Section 4.1). Describing the likelihood function used in CNFs is more
involved, and we postpone its discussion to Section 4.2.

Existing training approaches for neural ODEs almost exclusively consider the
reduced version of (3), which is obtained by solving the ODE (1) for fixed
weights. Denoting the unique solution of the state equation (1) for weights
$\boldsymbol{\theta}$ by $\mathbf{y}(\boldsymbol{\theta})$, we obtain the
reduced optimization problem

$\min_{\boldsymbol{\theta}}\quad\mathcal{J}(\boldsymbol{\theta},\mathbf{y}(\boldsymbol{\theta})),$
(6)

which is an infinite-dimensional unconstrained problem in the controls
$\boldsymbol{\theta}\colon[0,T]\to\mathbb{R}^{n_{p}}$.

To transform the variational problem (6) into a finite-dimensional problem, it
is most common to discretize the control on a control grid (e.g.,
$\boldsymbol{\theta}$ at $t=0,1$ in Figure 1) and numerically approximate the
integrals in $\mathcal{J}$ using quadrature rules. We use this strategy and in
Section 3 focus on handling the state equation (1). To this end, two main
approaches exist; the optimize-discretize (Opt-Disc) approach optimizes the
differential equation and discretizes the ODE only for the optimal weights,
while the discretize-optimize (Disc-Opt) approach uses a numerical time
integrator to obtain a fully discrete version of (6) that is then solved to
determine the optimal weights.

##  3 Discretize-Optimize and Optimize-Discretize

The Disc-Opt approach to performing time integration discretizes the ODE (1)
in time and then optimizes over the discretization. This approach is common in
neural networks and easy to implement, especially when automatic
differentiation can be used. For example, a discrete ResNet with input
features (or initial conditions) $\mathbf{y}_{0}$ follows forward propagation
(2) as the explicit Euler method on a uniform time discretization with step
size $h$. The choice of fixed step size $h$ trades off the accuracy of the
solution with the amount of computation and overall training time.

In applications that require the neural ODE for prediction or inference, such
as time-series regression and CNFs, the step size must be chosen judiciously
to ensure the trained discrete model captures the properties of the continuous
model. Choosing a too large value of $h$, may, for example, yield a discrete
flow model with an inaccurate or suboptimal inverse. Choosing a too small
value of $h$ leads to greater computational costs. In contrast to selection of
other training hyperparameters, the step size choice provides the advantage of
monitoring the accuracy of the ODE solver. While we use a fixed step size $h$
for simplicity, we can obtain even more efficient approaches by combining
adaptive discretization schemes with the backpropagation used in Disc-Opt.

Following similar steps used in Gholami et al. (2019), we expose a crucial
difference in the gradient computation between the Disc-Opt and Opt-Disc
approaches. For ease of presentation, we use and assume that the forward
propagation is the forward Euler scheme (2) and that the objective function
consists of the time-series regression loss defined by (5) with no
regularization. Accurate gradients are critical in ensuring the efficiency of
gradient-based optimization algorithms, including Stochastic Gradient Descent
(SGD) (Robbins and Monro, 1951) and ADAM (Kingma and Ba, 2014). The updates
for $\boldsymbol{\theta}$ in such methods depend on the gradient of the
objective function in (6).

In Disc-Opt, the backpropagation of the discrete ResNet (2) computes the
gradients. Automatic differentiation, which traverses the computational graph
backward in time, is used commonly in machine learning frameworks. The
discretization of the forward propagation completely determines this process.
For discrete objective function

$J(\boldsymbol{\theta})=h\sum_{i=1}^{N}L(\mathbf{y}_{i},\mathbf{u}_{i})$ (7)

and using auxiliary variable $\mathbf{z}$, the backpropagation through the
forward Euler discretization in (2) is

$\begin{split}\nabla_{\boldsymbol{\theta}_{j}}J(\boldsymbol{\theta})&=h\nabla_{\boldsymbol{\theta}}\,\ell(\boldsymbol{\theta}_{j},\mathbf{y}_{j},t_{j})\,\,\mathbf{z}_{j},\quad\text{where}\\\
\mathbf{z}_{N}&=h\,\nabla_{\mathbf{y}}L(\mathbf{y}_{N},\mathbf{u}_{N})\quad\text{and}\\\
\mathbf{z}_{j}&=\mathbf{z}_{j+1}+h\left(\nabla_{\mathbf{y}}\ell(\boldsymbol{\theta}_{j},\mathbf{y}_{j},t_{j})\mathbf{z}_{j+1}+\nabla_{\mathbf{y}}L(\mathbf{y}_{j},\mathbf{u}_{j})\right),\\\
\end{split}$ (8)

and the computations are backward through the layers, i.e.,
$j=N{-}1,N{-}2,\ldots,1$.

In Opt-Disc the gradients are computed by numerically solving the adjoint
equation

$\begin{split}\nabla_{\boldsymbol{\theta}(t)}\mathcal{J}(\boldsymbol{\theta})&=\mathbf{z}(t)\nabla_{\boldsymbol{\theta}}\ell\big{(}\boldsymbol{\theta}(t),\mathbf{y}(t),t\big{)},\quad\text{where}\\\
\mathbf{z}(T)&=\nabla_{\mathbf{y}}L\big{(}\mathbf{y}(T),\mathbf{u}(T)\big{)}\quad\text{and}\\\
-\partial_{t}\mathbf{z}&=\nabla_{\mathbf{y}}\ell\big{(}\boldsymbol{\theta}(t),\mathbf{y}(t),t\big{)}\mathbf{z}(t)+\nabla_{\mathbf{y}}L\big{(}\mathbf{y}(t),\mathbf{u}(t)\big{)}.\\\
\end{split}$ (9)

As indicated by the notation $-\partial_{t}$, this final value problem is
solved backward in time. This equation can be derived from the Karush-Kuhn-
Tucker (KKT) optimality conditions of the continuous learning problem
(Appendix A).

Comparison of the backpropagation (8) and the adjoint computation (9) shows
that both depend on the intermediate states, which need to be stored or
recomputed. However, the two differ because flexibility exists for choosing
the numerical scheme to discretize the adjoint computation in (9), whereas the
computation in (8) is determined by the discrete forward propagation. In fact,
the backpropagation shown in (8) can be seen as a discretization of the
adjoint equation; however, the standard backward Euler scheme reads

$\begin{split}\nabla_{\boldsymbol{\theta}_{j}}J(\boldsymbol{\theta})&=h\nabla_{\boldsymbol{\theta}}\ell\big{(}\boldsymbol{\theta}(t_{j}),\mathbf{y}(t_{j}),t_{j}\big{)}\mathbf{z}_{j},\quad\text{where}\\\
\mathbf{z}_{N}&=h\nabla_{\mathbf{y}}L\big{(}\mathbf{y}(t_{N}),\mathbf{u}(t_{N})\big{)}\quad\text{and}\\\
\mathbf{z}_{j}&=\mathbf{z}_{j+1}+h\big{(}\nabla_{\mathbf{y}}\ell\big{(}\boldsymbol{\theta}(t_{j+1}),\mathbf{y}(t_{j+1}),t_{j+1}\big{)}\mathbf{z}_{j+1}+\nabla_{\mathbf{y}}L\big{(}\mathbf{y}(t_{j+1}),\mathbf{u}(t_{j+1})\big{)}\big{)}.\\\
\end{split}$ (10)

The only difference between (8) and (10) is a shift of the indices of the
intermediate weights and features. Thus, the gradients obtained using both
methods differ unless both equations are solved accurately, such as when $h$
converges to zero. Hence, the Opt-Disc approach in Chen et al. (2018), which
uses adaptive time integrators for the forward and adjoint equations, may
provide inaccurate gradients when the time steps differ between both solvers
and the tolerance of those solvers is not sufficiently small. This comparison
is fairly standard and is also performed in Gholami et al. (2019).

The differences of the gradients computed using Disc-Opt and Opt-Disc affect
the convergence of neural ODEs for image classification (Gholami et al.,
2019). Since the gradient in the Disc-Opt method is an accurate estimate of
the gradient of the discrete objective function for any accuracy, the
differences are most striking when solving the forward propagation
inaccurately. Data science and scientific applications may desire an
inaccurate forward propagation since the data may be far less accurate than
the tolerance setting necessary for rendering a useful gradient in the Opt-
Disc approach. In the context of time-series regression and continuous
normalizing flows, the potential savings are limited since the discretization
trained in Disc-Opt must capture the relevant properties of the continuous
model.

##  4 Numerical Experiments

We compare the Disc-Opt approach with the Opt-Disc approach for two use cases:
time-series regression and density estimation through continuous normalizing
flows. For the former, we use the test data and experimental setup of the Opt-
Disc in Rackauckas et al. (2019) to obtain a direct comparison. For the
latter, we use a synthetic data set and five sizable public data sets using
CNFs as in Grathwohl et al. (2019).

###  4.1 Time-Series Regression

Figure 2: Time-series regression training iterations 100 and 300 comparing the
Disc-Opt and Opt-Disc approaches with the ground truth (11). Comparative
convergence video of the two methods is available at
https://imgur.com/nWxwVoe. Discrepancies between convergence behavior of the
approaches vary with initial parameterization (Section 4.1.1).

We consider the time-series regression problem in Chen et al. (2018) and
Rackauckas et al. (2019) to perform a direct comparison between Opt-Disc and
Disc-Opt. Given time-series data
$\mathbf{u}_{1}=\mathbf{u}(t_{1}),\ldots,\mathbf{u}_{n}=\mathbf{u}(t_{n})$
obtained from some unknown function $\mathbf{u}$, the goal is to tune the
neural ODE weights $\boldsymbol{\theta}$ in (1) such that
$\mathbf{y}(t_{k})\approx\mathbf{u}_{k}$ for $k=0,\ldots,n$. The initial value
of the neural ODE is $\mathbf{y}(0)=\mathbf{y}_{0}=\mathbf{u}_{0}$.

We use the data set generated in Rackauckas et al. (2019), where the data is
obtained from the ODE model

$\begin{cases}&\partial_{t}\mathbf{u}=\mathbf{A}\mathbf{u}^{\circ 3}\\\
&\mathbf{u}_{0}=\left[\begin{array}[]{c}2\\\ 0\end{array}\right]\\\
\end{cases}\quad,\qquad\text{where
\,}\mathbf{A}=\left[\begin{array}[]{cc}-0.1&2\\\ -2&-0.1\end{array}\right]$
(11)

where $\mathbf{u}^{\circ 3}$ denotes the element-wise cubic and $t\in[0,1.5]$.
Using the adaptive Tsitouras 5/4 Runge-Kutta method (dopri5), the ODE is
solved at 30 equidistantly-spaced time points (Figure 2). As a loss function,
we use the regression loss in (5).

We use the same neural ODE as in Rackauckas et al. (2019) that is based on the
layer

$\ell(\boldsymbol{\theta},\mathbf{y}_{0},t)\,\mathrel{\hbox
to0.0pt{\raisebox{1.29167pt}{$\cdot$}\hss}\raisebox{-1.29167pt}{$\cdot$}}=\,\mathbf{D}_{\mathbf{\theta}_{2}}\circ\tanh\circ\,\mathbf{D}_{\mathbf{\theta}_{1}}(\mathbf{y}_{0}^{\circ
3}).$ (12)

Here, $\mathbf{D}_{\mathbf{\theta}_{i}}$ denotes a single linear layer that
performs an affine transformation of the inputs and is parameterized by
$\boldsymbol{\theta}_{i}$ for $i=1,2$. The neural ODE is trained using 300
steps of the ADAM optimizer with step size $0.1$, starting with a Glorot
initialization (Glorot and Bengio, 2010). The Opt-Disc and Disc-Opt approaches
use the same initializations (i.e., same random seed).

The Opt-Disc approach in Rackauckas et al. (2019) uses the default dopri5
solver. We compare this with a Disc-Opt approach that uses Runge-Kutta 4 with
a fixed step size $h$ for training the neural ODE. For the step size $h$, we
use the spacing of the data points, i.e., $h=\frac{1.5}{29}$. With this
choice, the evaluation of the loss function does not require any
interpolation.

The Disc-Opt approach substantially reduces the cost of training. We attribute
these savings to the following two reasons. First, the runtime per iteration
is 97% lower than the Opt-Disc approach’s runtime per iteration (Figure 3).
This speedup mainly comes from the fewer function evaluations, pre-determined
by the step size choice, used in the Disc-Opt approach. Since the Opt-Disc
method uses adaptive time-stepping, the number of function evaluations, and
thus, cost per iteration varies during the optimization. Second, the Disc-Opt
approach also converges in roughly one third of the total iterations (Figure
3) with predictions at iteration 100 drastically closer to the ground truth
than the Opt-Disc approach (Figure 2).

Figure 3: For time-series regression, the Disc-Opt method converges in fewer
iterations, and each of its iterations requires less time. The mean iteration
clocktimes are 2.0 ms for Disc-Opt and 80.5 ms for Opt-Disc. Around iteration
24, the ODE solvers struggle to improve the training (cf. Figure 4), but the
Opt-Disc approach appears to suffer more.

To explain the fewer iterations needed by Disc-Opt, we numerically test the
quality of the gradients provided by both approaches through use of Taylor’s
theorem. Let $\mathbf{g}\in\mathbb{R}^{n}$ denote the gradient of the
objective function $\mathbf{F}\colon\mathbb{R}^{n}\rightarrow\mathbb{R}$ at a
point $\boldsymbol{\theta}$ and let $\mathbf{v}\in\mathbb{R}^{n}$ be a
randomly chosen direction. Then, by Taylor’s theorem, we have that

$\displaystyle E_{0}(h)$ $\displaystyle\mathrel{\hbox
to0.0pt{\raisebox{1.29167pt}{$\cdot$}\hss}\raisebox{-1.29167pt}{$\cdot$}}=\lVert\mathbf{F}(\boldsymbol{\theta}+h\mathbf{v})-\mathbf{F}(\boldsymbol{\theta})\rVert=\mathcal{O}(h\|\mathbf{v}\|)\quad\text{and}$
(13) $\displaystyle E_{1}(h)$ $\displaystyle\mathrel{\hbox
to0.0pt{\raisebox{1.29167pt}{$\cdot$}\hss}\raisebox{-1.29167pt}{$\cdot$}}=\lVert\mathbf{F}(\boldsymbol{\theta}+h\mathbf{v})-\mathbf{F}(\boldsymbol{\theta})-h\mathbf{g}^{\top}\mathbf{v}\rVert=\mathcal{O}(h^{2}\|\mathbf{v}\|).$

We observe the decay of $E_{0}(h)$ and $E_{1}(h)$ as $h\to 0$ for a fixed
(randomly chosen) direction $\mathbf{v}$ around the current network weights
for two different iterations of the training (Figure 4). We scale both axes
logarithmically. As expected, $E_{0}$ decays almost perfectly linearly, and
the decay of $E_{1}$ using the Disc-Opt method is approximately twice as steep
for large $h$, ultimately leveling off due to rounding errors and conditioning
(Figure 4 bottom row). While the derivative obtained in the Opt-Disc approach
is correct at iteration 14 (Figure 4 top left), it fails to provide an
accurate gradient at iteration 24 (Figure 4 top right), where the error
$E_{1}$ is greater than $E_{0}$ for all $h$ that we tested.

Figure 4: The derivative check (13) for time-series regression iterations 14
and 24 shown on log-log plot. As expected, the gradients of the Disc-Opt
approach (bottom row) are correct; note the faster decay of $E_{1}(h)$ (red
line) compared to $E_{0}(h)$ (blue line) as $h$ decays. The gradient in the
Opt-Disc approach (top row) is correct for iteration 14, but not for iteration
24. In this case (top right), the function $E_{1}$ (red line) is greater than
$E_{0}$ (blue line).

####  4.1.1 Extrapolation and Different Initial Conditions

For time-series regression on $t\in[0,1.5]$, the neural networks appear to
have modeled the ground truth ODE quite well (Section 4.1). We ask if the
models extrapolate correctly to time $t$ outside the training period
$[0,1.5]$. We extend the prediction of the trained models to the time period
$t\in[0,6]$ (Figure 5) in which we see that both the Disc-Opt and Opt-Disc
approaches learn behaviors that represent smooth ODEs that are different from
the ground truth used to create the data. This presents an example where
neural ODEs extrapolate poorly.

We repeat the time-series regression for several initial conditions to
demonstrate that the results span a broad set of convergent behaviors. We plot
three other random seeds (Figure 6). Sometimes the Disc-Opt and Opt-Disc
approaches converge at a similar rate per iteration and can have similar
extrapolations (Figure 6a). Some cases exist where the loss skyrockets, and
the models fail to fully converge in the 300 iterations (Figure 6b). Still
other cases appear where Opt-Disc converges to a lower loss than Disc-Opt
after 300 iterations, but extrapolation shows that the two models learn a
similar ODE (Figure 6c). In such cases, we observe that the Disc-Opt model
still trains more efficiently in terms of time. In fact, out of ten random
initial conditions selected, the Disc-Opt approach used fewer iterations to
converge than the Opt-Disc approach in nine of the comparisons. In the one
model (Figure 6c) where the Opt-Disc converges in fewer iterations than Disc-
Opt, the Disc-Opt model converges to the same loss in 466 iterations and still
reduced training time by 94%. Training costs reductions due to the Disc-Opt
model across these ten initial conditions ranged from 94% to 99% with an
average reduction of 97% (a 20x speedup).

Figure 5: Extrapolation of time-series regression models. After training each
model on the time interval $t\in[0,1.5]$, we visualize how the trained neural
networks extrapolate up to $t=6$. Each neural ODE models behaves as a smooth
ODE, by construction, though the learned models oscillate much more quickly
than the ground truth ODE.

(a) Different seed 1

(b) Different seed 2

(c) Different seed 3

Figure 6: Time-series regression convergence for different initial
parameterizations (cf. Figure 3, 5)

###  4.2 Density Estimation Using Continuous Normalizing Flows

Density estimation aims at constructing the probability density function from
a finite number of randomly drawn samples. Many approaches exist for this
problem including probabilistic models, variational autoencoders (Kingma and
Welling, 2014), and autoregressive flows (Gregor et al., 2014). Here, we
consider the normalizing flows (Rezende and Mohamed, 2015) approach, which can
be formulated using neural ODEs.

Normalizing flows parameterize a path for the samples from the unknown density
to the target density. The latter is usually chosen to be a standard normal
density, hence the name normalizing flow. In continuous normalizing flows
(CNFs), the model is a dynamical system formulated by a neural ODE in
artificial time (Chen et al., 2018). A key motivation for using a continuous
flow model is to obtain a generative model when the network is invertible.
Computing the flow backward in time, one can compute the density of and
compute samples from the unknown density. Furthermore, the continuous view
allows for free-form network architectures (Grathwohl et al., 2019).

The CNF is a smooth bijective mapping between unknown density $\rho_{0}$ and
known density $\rho_{1}$, is modeled by a neural ODE (2), and is trained using
samples. In the discrete setting, normalizing flows map samples
$\mathbf{y}_{0}\sim\rho_{0}(\mathbf{y}_{0})$ to points $\mathbf{y}_{N}$ that
follow $\rho_{1}$. In the continuous framework (1), the CNF maps
$\mathbf{y}(0)$ to $\mathbf{y}(T)=f(\boldsymbol{\theta},\mathbf{y}(0))$. For
implementation ease, we forward propagate $f$ from unknown density $\rho_{0}$
to the standard normal density $\rho_{1}$ (Figure 1).

Consider the likelihood $p(\mathbf{y}(t))$ where
$p\colon\mathbb{R}^{n_{f}}\rightarrow\mathbb{R}$. The likelihood
$p(\mathbf{y}(t))$ varies continuously through time with features $\mathbf{y}$
following the neural ODE (1). The initial likelihood $p(\mathbf{y}(0))$ is the
unknown density $\rho_{0}(\mathbf{y}(0))$, and we want
$p(\mathbf{y}(T))=p(f({\boldsymbol{\theta}},\mathbf{y}(0)))$ to approximate
$\rho_{1}(\mathbf{y}(T))$. From this initial condition and the instantaneous
change of variables (Chen et al., 2018), we obtain the initial value problem

$\displaystyle\partial_{t}[\,\log
p(\mathbf{y}(t))-\log\rho_{0}(\mathbf{y}(0))\,]$
$\displaystyle=-\textrm{Tr}\,\,\nabla_{\mathbf{y}(t)}\ell\big{(}\boldsymbol{\theta}(t),\mathbf{y}(t),t\big{)}$
(14) $\displaystyle p(\mathbf{y}(0))-\rho_{0}(\mathbf{y}(0))$
$\displaystyle=0.$

To calculate how the likelihood varies over the time domain $[0,T]$, we
integrate the log-likelihood ODE (14) through time and rearrange to

$\log
p(f(\boldsymbol{\theta},\mathbf{y}(0)))=\log\rho_{0}(\mathbf{y}(0))-\int_{0}^{T}\textrm{Tr}\,\,\nabla_{\mathbf{y}(t)}\ell\big{(}\boldsymbol{\theta}(t),\mathbf{y}(t),t\big{)}\,\mathrm{d}t.$

We want $\log p(f(\boldsymbol{\theta},\mathbf{y}(0)))$ to match $\log\rho_{1}$
where $\rho_{1}$ is the standard normal density. Hence, the loss function,
motivated by the Kullback–Leibler divergence (Rezende and Mohamed, 2015), is
given by the negative log-likelihood of the final state as

$\mathcal{L}(\boldsymbol{\theta},\mathbf{y})=\frac{n_{f}}{2}\log(2\pi)+\frac{1}{2}\|\mathbf{y}(T)\|^{2}-\int_{0}^{T}\textrm{Tr}\,\,\nabla_{\mathbf{y}(t)}\ell\big{(}\boldsymbol{\theta}(t),\mathbf{y}(t),t\big{)}\,\mathrm{d}t.$
(15)

We augment our initial ODE formulation (1) with the time integration (14).
Specifically, we vertically concatenate the states $\mathbf{y}$ with the log
likelihood $p(\mathbf{y}(0))$ and solve both as a single initial value problem
with the ODE solver (Grathwohl et al., 2019). As before, we minimize the loss
$\mathcal{L}(\boldsymbol{\theta},\mathbf{y})$ with respect to
$\boldsymbol{\theta}$. A meaningful model must be invertible, i.e., the
reverse mode must also push-forward $\rho_{1}$ to $\rho_{0}$, which is not
enforced in training. We stress the importance of checking model invertibility
(Tables 1, 2).

Data Set | $T$ | Model | Testing | $\frac{\text{Time}}{\text{Epoch}}$ (s) | Average | Inverse
---|---|---|---|---|---|---
Loss | NFE-F | Error
Gaussian | 0.5 | Opt-Disc (dopri5) | 2.83 | 0.79 | 65 | 8.88e-7
Mixture | Disc-Opt (rk4 $h=0.05$) | 2.79 | 0.48 | 40 | 1.82e-8
Disc-Opt (rk4 $h=0.25$) | 2.83 | 0.10 | 8 | 1.47e-2
Table 1: CNF results for the Gaussian Mixture problem. The Disc-Opt approach
with $h=0.25$ has the lowest training time and achieves a spuriously low loss
value; however, because the step size is too large, the discrete model loses
its invertibility (Section 4.2.2). Decreasing the time step to $h=0.05$ still
leads to a substantial reduction in the number of function evaluations,
comparable loss, and comparable inverse error than the Opt-Disc approach.

####  4.2.1 Normalizing Flow for Gaussian Mixture

To help visualize CNFs, consider the synthetic test problem where $\rho_{0}$
is the Gaussian mixture obtained by averaging eight bivariate Gaussians
situated in a circular pattern about the origin (Figure 1). For the Disc-Opt
example, we use a Runge-Kutta 4 model with constant step size of $h{=}1$ and
$T{=}2$. We discretize the control layers at $\boldsymbol{\theta}(0)$ and
$\boldsymbol{\theta}(1)$, while the state layers are discretized at one-
quarter intervals as determined by the Runge-Kutta 4 scheme. Overall, the
forward pass requires eight evaluations of $\ell$ and the weights at
intermediate time points are obtained by interpolating the two control layers
(illustration in Figure 1).

We use the neural ODE proposed in Grathwohl et al. (2019), where $\ell$ is
given by

$\ell(\boldsymbol{\theta}(t),\mathbf{y}(t),t)=c_{64,2}(\boldsymbol{\theta}_{4}(t),\cdot,t)\circ
c_{64,64}(\boldsymbol{\theta}_{3}(t),\cdot,t)\circ
c_{64,64}(\boldsymbol{\theta}_{2}(t),\cdot,t)\circ
c_{2,64}(\boldsymbol{\theta}_{1}(t),\mathbf{y}(t),t).$

The so-called concatsquash layer
$c_{i,j}(\boldsymbol{\theta}(t),\mathbf{y}(t),t)$ maps features
$\mathbf{y}(t)\in\mathbb{R}^{i}$ to outputs in $\mathbb{R}^{j}$ and is defined
as

$c_{i,j}(\boldsymbol{\theta},\mathbf{y},t)=(\mathbf{D}_{\mathbf{\theta}_{2}}\mathbf{y})(\sigma\circ\mathbf{D}_{\mathbf{\theta}_{1}}t)+(\mathbf{D}_{\mathbf{\theta}_{0}}t).$
(16)

Here, the nonlinear activation function $\sigma$ is the hyperbolic tangent and
$\mathbf{D}_{\mathbf{\theta}_{0}}$ uses no bias. This layer accepts and
operates on the space features and the time.

To solve the problem, the CNF needs to be trained to push the Gaussian mixture
$\rho_{0}$ to the center to form $\rho_{1}$ (Figure 1). To discretize the loss
function, we use samples drawn from the Gaussian mixture $\rho_{0}$. The
reverse mode of the model will separate the large Gaussian back into the eight
terms of the Gaussian mixture. What makes the problem challenging is that the
true inverse is discontinuous at the origin.

To gauge the accuracy of the trained flow model, we compute the inverse error

$\Delta_{\rm inv}(\mathbf{y}_{0})=\lVert
f^{-1}\left(\boldsymbol{\theta},f(\boldsymbol{\theta},\mathbf{y}_{0})\right)-\mathbf{y}_{0}\rVert,$
(17)

where $\mathbf{y}_{0}$ is sampled from $\rho_{0}$. This calculates the
Euclidean distance between an initial point $\mathbf{y}_{0}$ and the result
from mapping $\mathbf{y}_{0}$ forward then back. We report the average
Euclidean distance for the samples in one batch (Tables 1, 2); the Disc-Opt
and Opt-Disc approaches use the same batch.

A relatively small network model consisting of one control layer accurately
solves the problem (Table 1). For a fixed 10,000 iterations, the Disc-Opt
approach with $h=0.05$ demonstrates a slight speedup per iteration and
achieves a similar testing loss as the Opt-Disc approach; also, the Disc-Opt
solution has a comparable inverse error.

####  4.2.2 Too Coarse Time Stepping

Using a too large step size, $h=0.25$, the Disc-Opt model may achieve a low
testing loss while losing the invertibility. In fact, this flow fails to
accurately align the densities.

A sufficiently small step size is needed to ensure the discrete forward
propagation adequately represents the neural ODE (Ascher and Greif, 2011);
however, a too small step size leads to excessive and unnecessary computation
and training time. Hence, there is a trade-off between accuracy and speed. For
the Gaussian mixture problem, we trained two Disc-Opt networks with step sizes
$h=0.05$ and $h=0.25$ (Table 1). The $h=0.25$ step size showed an inverse
error much higher than the other two because the larger step size defines a
coarser state grid, which leads to a less-accurate approximation for the ODE
(1).

The coarse model ($h=0.25$) achieves a loss on par with the fine grid models
and an inverse error that is accurate for a couple decimal points (Table 1).
However, visually, the forward propagation leaves artifacts and does not
create a smooth Gaussian (Figure 7b). Also, the inverse flow does not fully
match the initial distribution. Therefore, we observe issues with a coarse
discrete model that are not indicated by the loss. Re-discretizing the coarse
model after training, we can switch to a finer solver for testing (Figure 7c),
which smooths some of the artifacts in the central Gaussian and curves the
inverse plot. The coarse model without changes (Figure 7b) may be sufficient
for some application because the inverse error could be sufficiently small
(Figure 7d).

(a) Disc-Opt RK4 with $h=0.05$

(b) Disc-Opt RK4 with $h=0.25$

(c) Testing for Disc-Opt RK4 with $h=0.05$ (but the model was trained on RK4
with $h=0.25$)

(d) Invertibility of Disc-Opt RK4 with $h=0.25$ for one point
$\hat{\mathbf{y}}(0)$

Figure 7: Using a step size not sufficiently small enough in the Disc-Opt
method can result in lack of invertibility. (a) The fine Disc-Opt approach
demonstrates a good model. (b) The coarse model fails to make a smooth
Gaussian and has some inaccuracies with the inverse. (c) We can switch the
testing solver for that coarse model to mitigate some of the issues. (d) A
representation of the point $\mathbf{y}(0)=(-2,-2)$ through the propagation
shows the scale of the inverse error. The inverse error of this point is less
than the average for the model (Table 1) so we expect that other points
demonstrate a greater inverse error.

####  4.2.3 High-Dimensional CNF Examples

We now consider some common higher-dimensional instances of the CNF problem
arising in density estimation; the public data sets we use are listed in Table
2. A detailed description of the data sets is given in Papamakarios et al.
(2017). In short, Miniboone, Power, Hepmass, and Gas are data sets compiled
for various tasks and housed by the University of California, Irvine (UCI)
Machine Learning Repository. The contexts include neutrinos, electric power
consumption of houses, particle-producing collisions, and gas sensors,
respectively. The Berkeley Segmentation Data Set (Bsds300) contains image
patches for image segmentation and boundary detection. Papamakarios et al.
(2017) preprocessed these data sets for use in density estimation. These data
sets have also been used to validate the Opt-Disc method in Grathwohl et al.
(2019).

Given the samples in the data sets, the goal in CNFs is to estimate their
underlying density $\rho_{0}$ by constructing a mapping to the standard normal
density with density $\rho_{1}$. We quantify this using the loss (15) computed
on a hold-out testing set. Since the goal is to invert the CNF to characterize
$\rho_{0}$, we numerically evaluate the inverse error (Table 2).

For each data set, we define an $\ell$ via concatenations of concatsquash
layer $c_{i,j}(\boldsymbol{\theta}(t),\mathbf{y}(t),t)$ (16) that are more
complicated than for the Gaussian mixture problem. For hyperparameters, we can
select an activation function $\sigma$, number of hidden layers $n_{h}$,
number of flow steps $s$, and the number of hidden dimensions, which we choose
based on some multiplier $m$ of the number of input features (Table 3). Our
general layer $\ell$ then is

$\ell(\boldsymbol{\theta}(t),\mathbf{y}(t),t)=\left[c_{mn_{f},n_{f}}(\boldsymbol{\theta}(t),\cdot,t)\circ
c_{mn_{f},mn_{f}}(\boldsymbol{\theta}(t),\cdot,t)^{n_{h}{-}1}\circ
c_{n_{f},mn_{f}}(\boldsymbol{\theta}(t),\mathbf{y}(t),t)\right]^{s}.$

Each of the $s$ blocks starts with one concatsquash layer that maps from
$\mathbb{R}^{n_{f}}$ to $\mathbb{R}^{mn_{f}}$, then $n_{h}{-}1$ more
concatsquash layers that maintain the dimensionality, then a final
concatsquash layer that maps back to $\mathbb{R}^{n_{f}}$. Depending on the
difficulty of the data set, different hyperparameters are selected; we mostly
followed the hyperparameters (Table 3) tuned by FFJORD (Grathwohl et al.,
2019).

The large data sets have higher dimensionality. As such, computing the full
trace (14) during training becomes onerous. Using Hutchinson’s trace estimator
(Hutchinson, 1989)

$\textrm{Tr}(A)=\mathbb{E}_{\phi(\epsilon)}\left[\epsilon^{\top}A\epsilon\right]$

for square matrix $A$ and noise vector $\epsilon$ with density
$\phi(\epsilon)$ where $\mathbb{E}[\epsilon]=0$ and
$\textrm{Cov}(\epsilon)=I$, we can make training less computationally
expensive and less time-consuming (Grathwohl et al., 2019). During testing,
the model computes the full trace without estimation. The Opt-Disc approach
uses trace estimation and a noise vector $\epsilon$ resampled for each batch.
For fair comparison, we adhere to these design decisions for Disc-Opt.

The performance on the testing sets of all data sets demonstrates that the
Opt-Disc and Disc-Opt approaches result in similar negative log-likelihood
loss (Table 2). However, Disc-Opt does so more quickly on most of the data
sets but with greater inverse error. Although the automatic differentiation
used in Disc-Opt is faster than the adjoint-based recalculation used in Opt-
Disc for batches of the same size, the memory requirements of storing
intermediates in automatic differentiation restrict the Disc-Opt batch sizes
to be much smaller than the Opt-Disc batch sizes. Checkpointing can also
handle the memory constraints. Nonetheless, the Opt-Disc method for solving
the Hepmass density estimation trains faster than Disc-Opt because of the
larger batches. Since we use the same models, code, and hyperparameters (Table
3) as FFJORD (Grathwohl et al., 2019), we witness similar testing losses
(Table 4).

We see Disc-Opt converge in less time for most data sets, but no timing payoff
for Hepmass (Table 2). Even so, on average for the five high-dimensional data
sets, Disc-Opt offers a 6x speedup in training over the Opt-Disc method. All
Opt-Disc and Disc-Opt models are competitive in performance with other density
estimation models (Table 4).

Model | Training | Inference | Total Training | Training | Inference | Testing | Inverse
---|---|---|---|---|---|---|---
Solver | Solver | Time (s) | Avg NFE-F | NFE-F | Loss | Error
POWER, $T{=}5$
Opt-Disc | dopri5 | dopri5 | 272K | 511 | 649 | -0.42 | 1.98e-7
Disc-Opt | rk4 $h{=}0.10$ | dopri5 | 56.0K | 200 | 2066 | -0.25 | 1.58e-5
Disc-Opt | rk4 $h{=}0.10$ | rk4 $h{=}0.10$ | 56.0K | 200 | 200 | -0.33 | 4.23e-3
GAS, $T{=}5$
Opt-Disc | dopri5 | dopri5 | 205K | 454 | 527 | -10.53 | 2.71e-7
Disc-Opt | rk4 $h{=}0.10$ | dopri5 | 121K | 200 | 437 | -10.27 | 4.71e-7
Disc-Opt | rk4 $h{=}0.10$ | rk4 $h{=}0.10$ | 121K | 200 | 200 | -10.27 | 6.63e-3
HEPMASS, $T{=}2$
Opt-Disc | dopri5 | dopri5 | 186K | 739 | 866 | 16.61 | 9.26e-7
Disc-Opt | rk4 $h{=}0.10$ | dopri5 | 281K | 400 | 765 | 16.60 | 1.09e-6
Disc-Opt | rk4 $h{=}0.10$ | rk4 $h{=}0.10$ | 281K | 400 | 400 | 16.60 | 2.25e-4
MINIBOONE, $T{=}1$
Opt-Disc | dopri5 | dopri5 | 36.1K | 113 | 132 | 10.64 | 1.85e-7
Disc-Opt | rk4 $h{=}0.25$ | dopri5 | 3.77K | 16 | 118 | 10.50 | 2.17e-7
Disc-Opt | rk4 $h{=}0.25$ | rk4 $h{=}0.25$ | 3.77K | 16 | 16 | 10.45 | 1.45e-3
BSDS300, $T{=}2$
Opt-Disc | dopri5 | dopri5 | 767K | 377 | 544 | -142.91 | 1.06e-7
Disc-Opt | rk4 $h{=}0.05$ | dopri5 | 49.7K | 160 | 511 | -146.23 | 1.82e-7
Disc-Opt | rk4 $h{=}0.05$ | rk4 $h{=}0.05$ | 49.7K | 160 | 160 | -146.14 | 3.75e-4
Table 2: CNF density estimation on large public data sets. Negative log-likelihood loss (lower is better) in nats (natural unit of information). The number of functions evaluations for the forward propagation (NFE-F) are averaged over all the batches and epochs. Performance is competitive with other methods (Table 4). The Opt-Disc approach for Bsds300 was terminated prematurely after nearly nine days of training. While the discrete model in Disc-Opt can have a large inverse error with the step size used in training, re-discretizing in the inference phase can reduce this error substantially. In summary, Disc-Opt reduces the training times and yields invertible models with comparable performance to Opt-Disc.  Data Set | Model | step size | batchsize | nonlinearity | # layers | hidden dim | # flow
---|---|---|---|---|---|---|---
multiplier | steps
POWER | Opt-Disc | - | 30000 | tanh | 3 | 10 | 5
Disc-Opt | 0.10 | 10000
GAS | Opt-Disc | - | 5000 | tanh | 3 | 20 | 5
Disc-Opt | 0.10 | 5000
HEPMASS | Opt-Disc | - | 10000 | softplus | 2 | 10 | 10
Disc-Opt | 0.10 | 5000
MINIBOONE | Opt-Disc | - | 5000 | softplus | 2 | 20 | 1
Disc-Opt | 0.25 | 5000
BSDS300 | Opt-Disc | - | 10000 | softplus | 3 | 20 | 2
Disc-Opt | 0.05 | 500
Table 3:  CNF hyperparameters. Mostly a replication of Grathwohl et al. (2019)’s table. All used initial learning rate=1e-3. Disc-Opt often used smaller batch sizes so that the batch would fit into the Titan X GPU (memory limit: 12GB). | POWER | GAS | HEPMASS | MINIBOONE | BSDS300
---|---|---|---|---|---
Real NVP | -0.17 | -8.33 | 18.71 | 13.55 | -153.28
Glow | -0.17 | -8.15 | 18.92 | 11.35 | -155.07
FFJORD | -0.46 | -8.59 | 14.92 | 10.43 | -157.40
our FFJORD (Opt-Disc) | -0.42 | -10.53 | 16.61 | 10.64 | -142.91
our FFJORD (Disc-Opt) | -0.25 | -10.27 | 16.60 | 10.50 | -146.23
MADE | 3.08 | -3.56 | 20.98 | 15.59 | -148.85
MAF | -0.24 | -10.08 | 17.79 | 11.75 | -155.69
TAN | -0.48 | -11.19 | 15.12 | 11.01 | -157.03
MAF-DDSF | -0.62 | -11.96 | 15.09 | 8.86 | -157.73
Table 4: CNF comparison with other methods. Replicated from Grathwohl et al.
(2019) with the addition of the two models we trained. Our Opt-Disc models
should be the exact same as FFJORD; minor differences attributed to the
randomness of batching and stochastic subgradient optimization methods.

####  4.2.4 Re-Discretizing the Trained Flows

Re-discretization in the inference phase presents a straightforward idea for
reducing the inversion error of Disc-Opt trained models.

Using the Miniboone data set, we consider training a CNF with the Disc-Opt
approach with Runge-Kutta 4 and $h{=}0.25$. Training takes less than 3,800
seconds, roughly one-tenth of the time required by the Opt-Disc (Table 2).
After training (e.g., during inference or model deployment) we re-discretize
the model using step size $h{=}0.05$ (Table 5). This maintains a similar
testing loss of 10.50 while substantially reducing the inverse error.
Demonstrating the flexible choice of solver, we also consider the adaptive
dopri5 solver in the evaluation (which has 118 function evaluations).

Since the Disc-Opt model was trained using a fixed step size, the
generalization with respect to re-discretization is remarkable; note the low
loss observed for all time integrators. As to be expected, the inverse error
is reduced for more accurate time integrators. Less accurate time integrators
can use their error to their advantage and result in lower loss values while
sacrificing invertibility. Therefore, both loss and invertibility should be
used to evaluate CNF performance.

| Trained via Disc-Opt | Trained via Opt-Disc
---|---|---
| Testing | NFE-F | Inverse | Testing | NFE-F | Inverse
| Loss | Error | Loss | Error
MINIBOONE, $T=1$ |  |  |  |  |  |
dopri5 rtol 1e-1 atol 1e-3 | 10.075 | 28 | 4.25e-2 | 10.108 | 29 | 3.16e-2
dopri5 rtol 1e-2 atol 1e-4 | 10.529 | 38 | 3.55e-3 | 10.669 | 38 | 1.61e-2
dopri5 rtol 1e-6 atol 1e-8 | 10.502 | 118 | 2.17e-7 | 10.636 | 132 | 1.85e-7
rk4 $h$=0.05 | 10.502 | 80 | 3.40e-7 | 10.636 | 80 | 5.22e-7
rk4 $h$=0.25 | 10.454 | 16 | 1.45e-3 | 10.590 | 16 | 2.57e-3
rk4 $h$=0.50 | 10.122 | 8 | 2.97e-2 | 10.239 | 8 | 4.79e-2
rk4 $h$=1.0 | 8.916 | 4 | 2.77e-1 | 9.082 | 4 | 3.70e-1
Table 5: Re-discretization of the trained CNF in the inference phase. Model
trained using the solver associated with bold values. All other results in the
column result from changing the solver settings (and thus the state
discretization) for the forward propagation when running on the test data. The
inverse error $\Delta_{\rm inv}$ averaged over one batch via (17). Since the
Disc-Opt problem is solved with a fixed step size and not using adaptive
integration used in Opt-Disc, the generalization to other discretization is
remarkable.

####  4.2.5 Multilevel Training

Models trained with finer grids tend to have better convergence than models
trained with coarser grids; however, these fine grid models consume much more
training time from the large NFE. Taking a multilevel approach, we train the
initial epochs with a coarse grid (few state layers), then add state layers to
make the grid finer for the final training epochs.

We test the multilevel strategy for the Disc-Opt and Opt-Disc approaches on
the Miniboone data set (Figure 8). For both approaches, we train the first 500
epochs with a coarse grid and initial learning rate 1e-3 then the next 1500
epochs with a finer grid (using the ODE solver in Table 2) while lowering the
learning rate to 5e-4. For Disc-Opt, we start with a Runge-Kutta 4 scheme
using $h{=}0.50$ then switch the step size to $h{=}0.25$. For Opt-Disc, we
start with dopri5 with relative tolerance 1e-3 and absolute tolerance 1e-5,
then switch to the default dopri5 (dividing each tolerance by $10^{3}$). For
both approaches, when we switch grids at epoch 500, we witness a large uptick
in the loss, which quickly decays. Overall, the performance of the multilevel
schemes is comparable in both cases, and a similar reduction of the training
time is observed in Figure 8.

Model | Total Training Time (s) | Testing Loss | Inverse Accuracy
---|---|---|---
Opt-Disc | 23.8K | 10.64 | 1.85e-7
Opt-Disc Multilevel | 20.2K | 10.43 | 1.22e-7
Disc-Opt | 2.73K | 10.45 | 1.45e-3
Disc-Opt Multilevel | 2.38K | 10.38 | 1.32e-3

Figure 8: Training CNFs on Miniboone with and without multilevel. For
multilevel, train first 500 epochs with a quicker coarse grid, then switch to
a fine grid for 1500 epochs. The switch to the finer grid results in immediate
uptick in loss, which quickly recovers. The convergence by epoch remains
similar to non-multilevel approaches, but the overall time is reduced in the
multilevel approach.

##  5 Discussion

We compare the discretize-optimize (Disc-Opt) and optimize-discretize (Opt-
Disc) approaches for time-series regression and continuous normalizing flows
(CNFs) using neural ODEs. The Disc-Opt approach achieves similar or superior
values of the validation loss at reduced computational costs, similar to the
classification results in Gholami et al. (2019). In our time-series regression
example, Disc-Opt achieves an average 20x speedup over Opt-Disc. For the
synthetic CNF problem of a two-dimensional Gaussian mixture model, Disc-Opt
requires about half the number of function evaluations. Similar results are
observed on more realistic large-scale applications. Consider, e.g., the data
set Bsds300. The Disc-Opt approach finished training in less than 14 hours,
while the Opt-Disc trained for more than 8 days. In other data sets, runtime
difference was less pronounced. However, the Disc-Opt method led to speedups
on all except Hepmass, with an average speedup of 6x in total. Since the
inverse error in this case is comparably small, we could perhaps use a larger
time step or adaptive time stepping, directions for future work. We could also
combine the numerical benefits of Disc-Opt with regularized neural ODEs, shown
to speedup training 2.8x in the Opt-Disc approach (Finlay et al., 2020).

One reason for these differences in training performance is the potential
inaccuracies of the gradients in Opt-Disc approaches, as also addressed in
Gholami et al. (2019). We review this known issue by comparing backpropagation
(used in Disc-Opt) to the adjoint method (used in Opt-Disc). In the latter,
gradients can become inaccurate when state and adjoint equations are not
solved well enough (Gholami et al., 2019). In Disc-Opt, the gradients are
accurate for any accuracy, allowing the user to guide solver accuracy by the
accuracy of the data; e.g., the noise level in scientific applications. We
also show this property numerically.

Our results are in line with Gholami et al. (2019), who showed the improved
convergence of Disc-Opt over the Opt-Disc approach on image classification
tasks. Here, we show similar properties for time-series regression and CNFs
using several numerical examples. These applications differ from
classification in that the inference in these applications requires the
continuous model. For example, in CNFs, the inverse of the model is required
for inference. The continuous model is invertible by design, which may
motivate one to prefer Opt-Disc over Disc-Opt. Remarkably, in our examples,
Disc-Opt approaches can achieve competitive loss values and (upon re-
discretization) low inversion errors. Notably, the Disc-Opt models in our
examples were faster to train. We also show that both approaches can be
accelerated in similar ways using multilevel training.

In conclusion, we advocate for Disc-Opt approaches for training neural ODEs in
time-series regression and continuous normalizing flows. In our view, the
inconvenience of choosing a sufficiently accurate time integration technique
is compensated by the computational savings; since the Disc-Opt settings
already exist in many codebases, implementation costs are minimal. Also, in
our experiments, the models trained with Disc-Opt generalize with respect to
re-discretization.

Acknowledgments

Supported by the National Science Foundation award DMS 1751636, AFORS grant
20RT0237, Binational Science Foundation Grant 2018209, and a GPU donation by
NVIDIA Corporation.

##  A  Derivation of Adjoint Equations

We provide derivations behind the continuous backpropagation (9) and discrete
backpropagations (8) and (10).

###  A.1 Continuous Adjoint

Consider the continuous ResNet optimization problem (5) subject to (1),

$\min_{\boldsymbol{\theta}}\,\int_{0}^{T}L\left(\mathbf{y}(t),\mathbf{u}(t)\right)\mathrm{d}t\quad\text{s.t.}\quad\partial_{t}\mathbf{y}(t)=\ell\big{(}\boldsymbol{\theta}(t),\mathbf{y}(t),t\big{)},\quad\mathbf{y}(0)=\mathbf{y}_{0}$
(18)

where $\mathbf{u}$ is the ground truth and $\mathbf{y}$ depends on
$\boldsymbol{\theta}$.

From here, we calculate the Lagrangian $\mathcal{G}$, with adjoint variable
(and Lagrangian multiplier) $\mathbf{z}$

$\displaystyle\mathcal{G}(\boldsymbol{\theta},\mathbf{y},\mathbf{z})$
$\displaystyle=\int_{0}^{T}L(\mathbf{y}(t),\mathbf{u}(t))\mathrm{d}t\,+\,\int_{0}^{T}\mathbf{z}(t)^{\top}\,\big{(}\ell(\boldsymbol{\theta}(t),\mathbf{y}(t),t)-\partial_{t}\mathbf{y}(t)\big{)}\mathrm{d}t.$
(19)
$\displaystyle=\int_{0}^{T}L(\mathbf{y}(t),\mathbf{u}(t))\,+\,\mathbf{z}(t)^{\top}\ell(\boldsymbol{\theta}(t),\mathbf{y}(t),t)-\mathbf{z}(t)^{\top}\partial_{t}\mathbf{y}(t)\,\,\mathrm{d}t.$

From optimization theory, we know that all the variations of $\mathcal{G}$
have to vanish at an optimal point. Here, we derive the strong form of the
optimality system. Before doing so, we simplify the third term using
integration by parts and simplify using $\mathbf{z}(0)=0$ (due to the initial
condition of the neural ODE):

$\int_{0}^{T}\mathbf{z}(t)^{\top}\partial_{t}\mathbf{y}(t)\,\mathrm{d}t=\mathbf{z}(t)^{\top}\mathbf{y}(t)\biggr{\rvert}_{0}^{T}-\int_{0}^{T}\partial_{t}\mathbf{z}(t)^{\top}\mathbf{y}(t)\,\mathrm{d}t=\mathbf{z}(T)^{\top}\mathbf{y}(T)-\int_{0}^{T}\partial_{t}\mathbf{z}(t)^{\top}\mathbf{y}(t)\,\mathrm{d}t.$
(20)

After substituting (20) into (19), the Lagrangian now becomes

$\displaystyle\mathcal{G}(\boldsymbol{\theta},\mathbf{y},\mathbf{z})$
$\displaystyle=\int_{0}^{T}L(\mathbf{y}(t),\mathbf{u}(t))\,+\,\mathbf{z}(t)^{\top}\ell(\boldsymbol{\theta}(t),\mathbf{y}(t),t)+\partial_{t}\mathbf{z}(t)^{\top}\mathbf{y}(t)\,\mathrm{d}t\,-\,\mathbf{z}(T)^{\top}\mathbf{y}(T).$
(21)

For some time $t\in[0,1)$, the variational derivative of $\mathcal{G}$ with
respect to $\mathbf{y}(t)$ is

$\displaystyle\partial_{\mathbf{y}(t)}\mathcal{G}(\theta,\mathbf{y},\mathbf{z})$
$\displaystyle=\nabla_{\mathbf{y}}L(\mathbf{y}(t),\mathbf{u}(t))\,+\,\nabla_{\mathbf{y}}\ell(\boldsymbol{\theta}(t),\mathbf{y}(t),t)\,\mathbf{z}(t)+\,\partial_{t}\mathbf{z}(t)$
(22)

Setting this equal to zero gives the backward in time ODE

$-\partial\mathbf{z}(t)=L(\mathbf{y}(t),\mathbf{u}(t))\,+\,\nabla_{\mathbf{y}}\ell(\boldsymbol{\theta}(t),\mathbf{y}(t),t)\mathbf{z}(t).$
(23)

Using that this equation holds for all $t\in[0,T)$, we see that the
variational derivative of $\mathcal{G}$ with respect to $\mathbf{y}(T)$ is

$\nabla_{\mathbf{y}(T)}\mathcal{G}(\theta,\mathbf{y},\mathbf{z})=\nabla_{\mathbf{y}}L(\mathbf{y}(T),\mathbf{u}(T))-\mathbf{z}(T).$
(24)

Setting this equal to zero gives the final time condition

$\mathbf{z}(T)=\nabla_{\mathbf{y}}L(\mathbf{y}(T),\mathbf{u}(T))$ (25)

Finally, we note that, via chain rule, the variation of $\mathcal{J}$ with
respect to $\boldsymbol{\theta}$ is

$\nabla_{\boldsymbol{\theta}(t)}\mathcal{J}(\boldsymbol{\theta})=\nabla_{\boldsymbol{\theta}}\ell(\boldsymbol{\theta}(t),\mathbf{y}(t),t)\mathbf{z}(t).$
(26)

###  A.2 Discrete Adjoints

The backpropagations presented in (8) and (10) come from discretizing the
continuous adjoint (9).

We consider a deep discrete ResNet (forward Euler scheme) and calculate the
Disc-Opt discretization for the backpropagation (8). The forward mode (2)
follows

$\displaystyle\mathbf{y}_{1}$
$\displaystyle=\mathbf{y}_{0}+h\,\ell(\boldsymbol{\theta}_{0},\mathbf{y}_{0},t_{0})$
$\displaystyle\mathbf{y}_{2}$
$\displaystyle=\mathbf{y}_{1}+h\,\ell(\boldsymbol{\theta}_{1},\mathbf{y}_{1},t_{1})$
$\displaystyle\vdots$ $\displaystyle\quad\qquad\vdots$
$\displaystyle\mathbf{y}_{N}$
$\displaystyle=\mathbf{y}_{N{-}1}+h\,\ell(\boldsymbol{\theta}_{N{-}1},\mathbf{y}_{N{-}1},t_{N{-}1}).$

From outputs $\mathbf{y}_{1},\dots,\mathbf{y}_{N}$, we calculate a discrete
sum as the loss

$J(\boldsymbol{\theta})=\sum_{i=1}^{N}h\,L(\mathbf{y}_{i},\mathbf{u}_{i}).$
(27)

We then backpropagate to calculate the gradients to update model parameters.
Using auxiliary term $\mathbf{z}$ to accumulate the gradient of $J$ with
respect to the states $\mathbf{y}$ as we step backwards in time,

$\displaystyle\mathbf{z}_{N}$
$\displaystyle=h\,\nabla_{\mathbf{y}_{N}}L(\mathbf{y}_{N},\mathbf{u}_{N})$
$\displaystyle\mathbf{z}_{j}$
$\displaystyle=h\,\nabla_{\mathbf{y}_{j}}L(\mathbf{y}_{j},\mathbf{u}_{j})+h\,\nabla_{\mathbf{y}_{j}}L(\mathbf{y}_{j+1},\mathbf{u}_{j+1})$
$\displaystyle=h\,\nabla_{\mathbf{y}_{j}}L(\mathbf{y}_{j},\mathbf{u}_{j})+\nabla_{\mathbf{y}_{j}}\mathbf{y}_{j+1}\,\,\,h\,\nabla_{\mathbf{y}_{j+1}}L(\mathbf{y}_{j+1},\mathbf{u}_{j+1})$
$\displaystyle=h\,\nabla_{\mathbf{y}_{j}}L(\mathbf{y}_{j},\mathbf{u}_{j})+(I+h\nabla_{\mathbf{y}_{j}}\ell(\boldsymbol{\theta}_{j},\mathbf{y}_{j},t_{j}))\,\mathbf{z}_{j+1}$
$\displaystyle=\mathbf{z}_{j+1}+h\big{(}\nabla_{\mathbf{y}}\ell(\boldsymbol{\theta}_{j},\mathbf{y}_{j},t_{j})\mathbf{z}_{j+1}+\nabla_{\mathbf{y}}L(\mathbf{y}_{j},\mathbf{u}_{j})\big{)}$

Using chain rule, we can now calculate the gradient with respect to the
parameters

$\displaystyle\nabla_{\boldsymbol{\theta}_{j}}J(\boldsymbol{\theta})$
$\displaystyle=\nabla_{\boldsymbol{\theta}}\,\mathbf{y}_{j}\,\,\nabla_{\mathbf{y}_{j}}J(\boldsymbol{\theta})$
(28)
$\displaystyle=h\nabla_{\boldsymbol{\theta}}\,\ell(\boldsymbol{\theta}_{j},\mathbf{y}_{j},t_{j})\,\,\mathbf{z}_{j}.$

The Opt-Disc discretization (10) follows the backward Euler scheme. It
calculates the gradient at time $t_{j{+}1}$ instead of at $t_{j}$ as Disc-Opt
does.

## References

  * Ascher (2020) U. M. Ascher.  Discrete processes and their continuous limits.  _Journal of Dynamics & Games_, 7:123, 2020.  doi: 10.3934/jdg.2020008.
  * Ascher and Greif (2011) U. M. Ascher and C. Greif.  _A First Course in Numerical Methods_.  SIAM, 2011.
  * Benning et al. (2019) M. Benning, E. Celledoni, M. J. Ehrhardt, B. Owren, and C.-B. Schönlieb.  Deep learning as optimal control problems: Models and numerical methods.  _Journal of Computational Dynamics_ , 6(2):171, 2019.
  * Bliss (1919) G. A. Bliss.  The use of adjoint systems in the problem of differential corrections for trajectories.  _JUS Artillery_ , 51:296–311, 1919.
  * Brunton et al. (2019) S. L. Brunton, B. R. Noack, and P. Koumoutsakos.  Machine learning for fluid mechanics.  _Annual Review of Fluid Mechanics_ , 52, 2019.
  * Chang et al. (2018) B. Chang, L. Meng, E. Haber, L. Ruthotto, D. Begert, and E. Holtham.  Reversible architectures for arbitrarily deep residual neural networks.  In _AAAI Conference on Artificial Intelligence_ , pages 1–8, 2018\.
  * Chen et al. (2018) T. Q. Chen, Y. Rubanova, J. Bettencourt, and D. K. Duvenaud.  Neural ordinary differential equations.  In _Advances in Neural Information Processing Systems (NeurIPS)_ , pages 6571–6583, 2018.
  * E (2017) W. E.  A proposal on machine learning via dynamical systems.  _Comm. Math. Statist._ , 5(1):1–11, 2017.
  * Finlay et al. (2020) C. Finlay, J.-H. Jacobsen, L. Nurbekyan, and A. M. Oberman.  How to train your neural ode.  _arXiv:2002.02798_ , 2020.
  * Gholami et al. (2019) A. Gholami, K. Keutzer, and G. Biros.  ANODE: Unconditionally accurate memory-efficient gradients for neural odes.  _International Joint Conferences on Artificial Intelligence (IJCAI)_ , 2019.
  * Glorot and Bengio (2010) X. Glorot and Y. Bengio.  Understanding the difficulty of training deep feedforward neural networks.  In _International Conference on Artificial Intelligence and Statistics (AISTATS)_ , pages 249–256, 2010.
  * González-García et al. (1998) R. González-García, R. Rico-Martínez, and I. G. Kevrekidis.  Identification of distributed parameter systems: A neural net based approach.  _Computers Chem Engn._ , 22:S965–S968, Mar. 1998.
  * Grathwohl et al. (2019) W. Grathwohl, R. T. Chen, J. Betterncourt, I. Sutskever, and D. Duvenaud.  FFJORD: Free-form continuous dynamics for scalable reversible generative models.  _International Conference on Learning Representations (ICLR)_ , 2019\.
  * Gregor et al. (2014) K. Gregor, I. Danihelka, A. Mnih, C. Blundell, and D. Wierstra.  Deep autoregressive networks.  In _International Conference on Machine Learning (ICML)_ , pages 1242–1250, 2014.
  * Günther et al. (2019) S. Günther, L. Ruthotto, J. B. Schroder, E. C. Cyr, and N. R. Gauger.  Layer-Parallel Training of Deep Residual Neural Networks.  _SIAM Journal on the Mathematics of Data Science_ , pages 1–23, 2019\.
  * Gunzburger (2003) M. D. Gunzburger.  _Perspectives in flow control and optimization_ , volume 5.  SIAM, 2003.
  * Haber and Ruthotto (2017) E. Haber and L. Ruthotto.  Stable architectures for deep neural networks.  _Inverse Problems_ , 34:014004, 2017.
  * He et al. (2016) K. He, X. Zhang, S. Ren, and J. Sun.  Deep residual learning for image recognition.  In _IEEE Conference on Computer Vision and Pattern Recognition (CVPR)_ , pages 770–778, 2016.
  * Hutchinson (1989) M. F. Hutchinson.  A stochastic estimator of the trace of the influence matrix for laplacian smoothing splines.  _Communications in Statistics-Simulation and Computation_ , 18(3):1059–1076, 1989.
  * Kingma and Ba (2014) D. P. Kingma and J. Ba.  Adam: A method for stochastic optimization.  _arXiv:1412.6980_ , 2014.
  * Kingma and Welling (2014) D. P. Kingma and M. Welling.  Auto-encoding variational bayes.  _International Conference on Learning Representations (ICLR)_ , 2014\.
  * Köhler et al. (2019) J. Köhler, L. Klein, and F. Noé.  Equivariant flows: sampling configurations for multi-body systems with symmetric energies.  _arXiv:1910.00753_ , 2019.
  * Krogh and Hertz (1992) A. Krogh and J. A. Hertz.  A simple weight decay can improve generalization.  In _Advances in Neural Information Processing Systems (NeurIPS)_ , pages 950–957, 1992.
  * Li et al. (2017) Q. Li, L. Chen, C. Tai, and E. Weinan.  Maximum principle based algorithms for deep learning.  _The Journal of Machine Learning Research (JMLR)_ , 18(1):5998–6026, 2017.
  * Lu et al. (2018) Y. Lu, A. Zhong, Q. Li, and B. Dong.  Beyond finite layer neural networks: Bridging deep architectures and numerical differential equations.  In _International Conference on Machine Learning (ICML)_ , Proceedings of Machine Learning Research, pages 3276–3285. PMLR, 10–15 Jul 2018\.
  * Nah et al. (2017) S. Nah, T. H. Kim, and K. M. Lee.  Deep multi-scale convolutional neural network for dynamic scene deblurring.  In _IEEE Conference on Computer Vision and Pattern Recognition (CVPR)_ , pages 257–265, July 2017.  doi: 10.1109/CVPR.2017.35.
  * Nocedal and Wright (2006) J. Nocedal and S. Wright.  _Numerical optimization_.  Springer Science & Business Media, 2006.
  * Papamakarios et al. (2017) G. Papamakarios, T. Pavlakou, and I. Murray.  Masked autoregressive flow for density estimation.  In _Advances in Neural Information Processing Systems (NeurIPS)_ , pages 2338–2347, 2017.
  * Rackauckas et al. (2019) C. Rackauckas, M. Innes, Y. Ma, J. Bettencourt, L. White, and V. Dixit.  DiffEqFlux.jl - A julia library for neural differential equations.  _arXiv:1902.02376_ , 2019.
  * Rackauckas et al. (2020) C. Rackauckas, Y. Ma, J. Martensen, C. Warner, K. Zubov, R. Supekar, D. Skinner, A. Ramadhan, and A. Edelman.  Universal differential equations for scientific machine learning.  _arXiv:2001.04385_ , 2020.
  * Raissi (2018) M. Raissi.  Deep hidden physics models: Deep learning of nonlinear partial differential equations.  _The Journal of Machine Learning Research (JMLR)_ , 19(1):932–955, 2018.
  * Rezende and Mohamed (2015) D. J. Rezende and S. Mohamed.  Variational inference with normalizing flows.  In _International Conference on Machine Learning (ICML)_ , pages 1530–1538, 2015.
  * Rico-Martínez et al. (1992) R. Rico-Martínez, K. Krischer, and I. G. Kevrekidis.  Discrete- vs. Continuous-time Nonlinear Signal Processing of Cu Electrodissolution Data.  _Chemical Engineering Communications_ , 118(1):25–48, 1992.
  * Robbins and Monro (1951) H. Robbins and S. Monro.  A stochastic approximation method.  _Ann. Math. Statist._ , 22(3):400–407, 09 1951\.  doi: 10.1214/aoms/1177729586.
  * Ruthotto and Haber (2019) L. Ruthotto and E. Haber.  Deep neural networks motivated by partial differential equations.  _Journal of Mathematical Imaging and Vision_ , Sep 2019.  doi: 10.1007/s10851-019-00903-1.
  * Sonoda and Murata (2019) S. Sonoda and N. Murata.  Transport analysis of infinitely deep neural network.  _The Journal of Machine Learning Research (JMLR)_ , 20(1):31–82, 2019.
  * Zhang et al. (2017) Q. Zhang, Z. Cui, X. Niu, S. Geng, and Y. Qiao.  Image segmentation with pyramid dilated convolution based on resnet and u-net.  In _International Conference on Neural Information Processing_ , pages 364–372. Springer, 2017.
  * Zhu et al. (2018) M. Zhu, B. Chang, and C. Fu.  Convolutional neural networks combined with Runge-Kutta methods.  _arXiv:1802.08831_ , Feb 2018.