∎ .tiffpng.pngconvert #1 .tiff

11institutetext: J. Baker, E. Cherkaev 22institutetext: Department of
Mathematics
University of Utah
22email<EMAIL_ADDRESS>22email<EMAIL_ADDRESS>33institutetext: A.
Narayan, B. Wang 44institutetext: Department of Mathematics, and Scientific
Computing and Imaging (SCI) Institute
University of Utah
44email<EMAIL_ADDRESS>44email<EMAIL_ADDRESS>
# Learning Proper Orthogonal Decomposition of Complex Dynamics Using Heavy-
ball Neural ODEs

Justin Baker  Elena Cherkaev  Akil Narayan  Bao Wang

(Received: date / Accepted: date)

###### Abstract

Proper orthogonal decomposition (POD) allows reduced-order modeling of complex
dynamical systems at a substantial level, while maintaining a high degree of
accuracy in modeling the underlying dynamical systems. Advances in machine
learning algorithms enable learning POD-based dynamics from data and making
accurate and fast predictions of dynamical systems. This paper extends the
recently proposed heavy-ball neural ODEs (HBNODEs) [Xia et al. NeurIPS, 2021]
for learning data-driven reduced-order models (ROMs) in the POD context, in
particular, for learning dynamics of time-varying coefficients generated by
the POD analysis on training snapshots constructed by solving full-order
models. HBNODE enjoys several practical advantages for learning POD-based ROMs
with theoretical guarantees, including 1) HBNODE can learn long-range
dependencies effectively from sequential observations, which is crucial for
learning intrinsic patterns from sequential data, and 2) HBNODE is
computationally efficient in both training and testing. We compare HBNODE with
other popular ROMs on several complex dynamical systems, including the von
Kármán Street flow, the Kurganov-Petrova-Popov equation, and the one-
dimensional Euler equations for fluids modeling.

###### Keywords:

Neural ODE Momentum Reduced-order modeling Deep learning

###### MSC:

65P99 68T07

##  1 Introduction

Numerical long-time simulation of full-order models (FOMs) of complex
dynamical systems is computationally costly. This is particularly true for
physical systems that contain a wide range of spatial and temporal scales,
including direct numerical simulation (DNS) moin1998direct  or large eddy
simulation (LES) in fluid mechanics germano1991dynamic ; you2007dynamic ;
craster2009dynamics  and chaotic systems 10.1143/PTPS.64.346 ; 10.2307/2100687
; SIVASHINSKY19771177 . Reduced-order models (ROMs) have been utilized as
alternative scientific simulation tools, which are computationally much more
efficient than FOMs and retain comparable accuracy for simulating complex
dynamical systems. ROMs have played crucial roles in designing, optimizing,
and controlling dynamical systems gugercin_survey_2004 ;
antoulas_approximation_2005 ; antoulas_interpolatory_2010 ; benner_survey_2015
.

Several data-driven numerical algorithms have been proposed for reduced-order
modeling, including dynamic mode decomposition (DMD) schmid_2010  and proper
orthogonal decomposition (POD) benner2015survey . These models leverage some
FOM simulation data to construct low-dimensional simplified models that
describe the underlying dynamics, with the goal of using these simplified
models in generalization regimes to predict the unseen dynamics. Classical
projection-based reduced-order modeling techniques (of which DMD and POD are
examples) are among the most popular approaches for constructing ROMs of
dynamical systems. This approach transforms the simulation results of FOM into
a suitable low-dimensional subspace that preserves the largest variance of the
training data. In, e.g., POD, classical numerical algorithms (such as Galerkin
methods), are subsequently used to rewrite the state variable in the governing
equation of the underlying dynamics into a system of ODEs, resulting in a
substantially reduced degree of freedom for describing the complex dynamics.
Both DMD and POD have been widely used in scientific simulations, particularly
for fluid simulations.

ROMs generated from projection-based approaches can preserve crucial physical
structures of the dynamics system. However, inappropriate truncation of the
POD modes in governing equations can severely degrade modeling accuracy and
result in unexpected, unphysical predictive results. Moreover, the precise
strategy for mode truncation is task-dependent and is typically limited to
explicit and closed definitions of the mathematical models SAN2019271 .
Another drawback of direct projection-based approaches is that they require
knowledge of governing equations that model the dynamical system, and this
information is often absent for real-world problems. As such, data-driven
reduced-order modeling has drawn significant recent attention. For instance,
the learning of closure models to compensate for information loss due to mode
truncation SAN2018681 ; san2017neural ; mou2020data ;
https://doi.org/10.1002/fld.4684 , and data-driven reduced basis
representations have been learned from simulation data that provides
significantly improved predictive performance of the dynamics compared to
classical models murata_fukami_fukagata_2020 ; lui2019construction ; 7572934 .
More recently, “vanilla” versions of machine learning approaches such as
neural ODEs (NODEs) and recurrent neural networks (RNNs) have been used to
learn temporal coefficients of the POD of a given complex dynamical system
rom_node ; kani2017dr ; kani2019reduced . A well-known issue of the vanilla
NODEs and RNNs is that they lack the capability of capturing long-range
dependencies from data, making these machine learning models fail to learn the
intrinsic physics of complex systems bengio1994learning ; HBNODE:2021 .

###  1.1 Our contribution

We employ the recently developed _heavy-ball neural ODE_ (HBNODE) HBNODE:2021
, an extension of NODE chen2018neural , to learn the temporal coefficients of
the POD of complex physical systems with a focus on time-dependent simulations
in scientific computing. In particular, our examples include the von Kármán
Street (VKS) flow, the Kurganov-Petrova-Popov (KPP) equation, and the one-
dimensional Euler equations for fluids modeling. We provide numerical
validation on the above three benchmark tasks and a detailed empirical and
analysis of why HBNODEs are beneficial for learning the dynamics of POD modes.
Our numerical results show the adjoint state of HBNODEs does not vanish,
confirming that HBNODEs do learn long-range dependencies, which results in
remarkable performance gain over the baseline NODEs. Moreover, our
experimental results show significant computational advantages in training and
testing HBNODEs over the baseline ROM models.

###  1.2 Related work

There is a healthy amount of recent work on learning POD mode dynamics using
deep neural networks, particularly RNNs and vanilla NODEs. Perhaps the most
related papers to this article are rom_node ; dutta2021neural ; dutta2021data
, which study the NODE framework for learning ROMs. In rom_node , the authors
developed a POD-NODE ROM framework for learning POD coefficients, which starts
from FOM snapshots and then uses an autoencoder to encode the POD
representations of FOM snapshots, followed by NODE training and forecasting.
The POD-NODE ROM framework achieves appealing results for learning reduced
dynamics of the VKS model, and it significantly outperforms the direct
application of a long short-term memory (LSTM) network for sequential
learning. In dutta2021neural ; dutta2021data , the authors study the
effectiveness of NODEs for reduced-order modeling and predicting environment
hydrodynamics. On the one hand, they find that NODEs provide an elegant
framework for the stable and accurate evolution of latent-space dynamics with
promising generalizability. On the other hand, they noticed that in order to
facilitate the widespread adoption of NODEs for large-scale systems,
significant effort needs to be directed at accelerating training time. This
limitation motivates this article’s study and utilization of HBNODEs
HBNODE:2021 , which is the machine-learning backbone of the reduced-order
modeling pipeline proposed in this work. There are three major advantages of
learning PODs using HBNODEs over the existing deep learning approaches:

  * •

HBNODEs are a class of continuous-depth neural networks, and they are suitable
for learning irregularly-sampled simulation data or physical observations.
Hence, observation protocols that entail missing or sparse data are easily
tackled in this framework.

  * •

Certain spectral properties of HBNODEs enable them to capture long-range
dependencies from sequential data, which is crucial for learning PODs of
complex dynamics.

  * •

Both HBNODEs and their adjoint ODEs are computationally much more efficient
than baseline NODEs.

In addition to the NODE paradigm of continuous-depth neural networks for
reduced-order modeling, the RNN — a natural sequential deep learning model —
has also been successfully used for learning-assisted model reduction. Many
advanced RNN algorithms can also be leveraged to enhance learning ROMs, e.g.,
LSTM networks LSTM . RNN-based ROMs have achieved remarkable success in many
applied domains, including multiphase flow simulation kani2017dr ;
kani2019reduced , learning advection-dominated systems maulik2021reduced ,
learning chaotic dynamics ma2018model , and learning nonlinear aeroelastic
models mannarino2014nonlinear . Compared to NODEs for learning ROMs, RNNs
cannot learn irregularly-sampled time series effectively and can even depart
from the underpinning physics due to their discrete nature.

###  1.3 Organization

We organize the paper as follows: In Sections 2 and 3, we briefly review the
POD-based reduced-order modeling and HBNODE for continuous-depth deep
learning, respectively. We present the benchmark physical models of the
complex dynamical systems and full-order modeling for data generation in
Section 4. Section 5 shows the detailed deep learning model and pipeline for
learning POD-based ROMs. We verify the efficacy of our proposed machine
learning models and contrast them with several baseline models in Section 6,
followed by concluding remarks.

##  2 POD-based Reduced-order Modeling

In this section, we briefly review key ideas and procedures of POD-based
reduced-order modeling.

###  2.1 Notation

We denote vectors and matrices by lower- and upper-case boldface letters,
respectively. For a vector
${\bm{x}}=(x_{1},\cdots,x_{d})^{\top}\in\mathbb{R}^{d}$, where
$(x_{1},\cdots,x_{d})^{\top}$ denotes the transpose of the row vector
$(x_{1},\cdots,x_{d})$, we use
$\|{\bm{x}}\|={(\sum_{i=1}^{d}x_{i}^{2})^{1/2}}$ to denote its $\ell_{2}$
norm, and use $\mathbf{0}$ to denote the zero vector. In cases when $d=2$ and
$\bm{x}$ is a spatial vector, we will write the components instead as
$\bm{x}=(x,y)^{\top}$. For a matrix ${\bm{A}}$, we use ${\bm{A}}^{\top}$,
${\bm{A}}^{-1}$, and $\|{\bm{A}}\|$ to denote its transpose, inverse, and
spectral norm, respectively. We use ${\bm{I}}$ to denote the identity matrix,
whose size will be clear based on context.

We will consider the approximation of a space-time function
${\bm{u}}={\bm{u}}({\bm{x}},t)$ where ${\bm{x}}$ is a spatial vector
(typically of 1 or 2 dimensions) and $t$ is a scalar on $[0,T]$ for some fixed
and finite terminal time $T$. The function ${\bm{u}}$ may be vector-valued. In
much of our discussion, we will take the concrete example of ${\bm{u}}$ being
a solution to a discretized VKS problem, whose details are given in Section
4.1. For the VKS problem, ${\bm{u}}\in{\mathbb{R}}^{2}$ contains the
horizontal ($x)$ and vertical ($y)$ components of a fluid velocity field. We
will write ${\bm{u}}=(u_{x},u_{y})$ to denote these two components.

###  2.2 POD snapshots

Figure 1: POD pipeline: We first pre-process the data, from experimental
observation or FOM simulation, by subtracting the mean. Then we apply spectral
decomposition of the covariance matrix and only keep the first $r$ eigenmodes,
resulting in the reduced representation —
$\sum_{i=1}^{r}\alpha_{i}(t)\psi_{i}({\bm{x}})$ — of
${\bm{u}}^{\prime}({\bm{x}},t)$, where $\alpha(t)$s are the temporal
coefficients and $\psi(x)$s are the eigenmodes.

POD shares a similar spirit and implementation as the celebrated principal
component analysis (PCA), the latter of which has been a very popular tool for
data analysis pearson1901liii ; LIANG2002527 . The key idea of PCA is to
project high-dimensional data into a lower-dimensional space that is spanned
by the eigenvectors corresponding to the leading eigenvalues of the covariance
matrix of the data. PCA preserves the largest variance of training data and
thus contains the most important information contained in the originally high-
dimensional data. POD has been introduced in accelerating fluids simulation
and reducing the complexity of fluid models since the pioneering work of
Berkooz et al. berkooz1993proper . Once the leading eigenmodes are obtained
via analysis of training data, it is possible to reduce the order (the
computational complexity and degrees of freedom) of the complex FOMs. The POD-
based dimension reduction approach starts with some training samples of
physically observed or numerically simulated snapshots of dynamics. These
sample snapshots are aggregated into an ensemble matrix ${\bm{Y}}$, where each
row contains the state of a dynamical system at a fixed time step. Next, we
compute the covariance matrix of the rows of the matrix ${\bm{Y}}$, and the
eigenvectors — sorted according to the corresponding decreasing-ordered
eigenvalues — are used as the new orthogonal basis for representing the ROM.
Below we summarize the crucial steps of identifying the low-dimensional
representations via the POD approach, which has also been visualized in Fig.
1.

  * •

_Data generation._ We simulate the FOM, which is computationally expensive,
for a short time to obtain the training data at time steps
$t_{1},t_{2},\cdots,t_{N_{t}}$. For the VKS problem, when our solution $u$
contains two components ${\bm{u}}=(u_{x},u_{y})$ that depend on the two-
dimensional spatial variable $\bm{x}=(x,y)$, we assume that FOM snapshots
$\bm{u}_{x}(t_{j})$, $\bm{u}_{y}(t_{j})$ are vectorized representations of $N$
spatial degrees of freedom. Then we have

$\displaystyle\bm{u}_{x}(t_{j}),\bm{u}_{y}(t_{j})$
$\displaystyle\in\mathbb{R}^{N},$ $\displaystyle j$
$\displaystyle=1,2,\ldots,N_{t}.$

For real-world dynamical systems for which we do not know the exact governing
equation, we sample the true dynamics via experimental measurements as
training data. In either case, we assume that training data is available to us
(as the snapshots above), and our goal is to efficiently leverage this data to
learn reduced-order dynamics without recourse to the FOM, which we assume is
unknown.

  * •

_Linearly center the data dynamics._ With our VKS-centric notation above,
according to the Reynolds decomposition of the flow, we have for fixed
$t_{j}$,

$\bm{u}_{x}=\bar{\bm{u}}_{x}+\bm{u}_{x}^{\prime};\quad\bm{u}_{y}=\bar{\bm{u}}_{y}+\bm{u}_{y}^{\prime},$
(1)

where $\bar{\bm{u}}_{x}$ and $\bar{\bm{u}}_{y}$ are the temporal mean of the
solutions, computed over our $N_{t}$ snapshots. The components
$\bm{u}_{x}^{\prime}$ and $\bm{u}_{y}^{\prime}$ are the fluctuating components
of the data.

  * •

_Data assembling._ Concatenate the simulated and centered FOM snapshots into
the following matrix ${\bm{Y}}$,

${\bm{Y}}=\left(\begin{array}[]{ccccccc}\rule[2.15277pt]{10.76385pt}{0.5pt}&(\bm{u}^{\prime}_{x}(t_{1}))^{\top}&\rule[2.15277pt]{10.76385pt}{0.5pt}&&\rule[2.15277pt]{10.76385pt}{0.5pt}&(\bm{u}^{\prime}_{y}(t_{1}))^{\top}&\rule[2.15277pt]{10.76385pt}{0.5pt}\\\
\rule[2.15277pt]{10.76385pt}{0.5pt}&(\bm{u}^{\prime}_{x}(t_{2}))^{\top}&\rule[2.15277pt]{10.76385pt}{0.5pt}&&\rule[2.15277pt]{10.76385pt}{0.5pt}&(\bm{u}^{\prime}_{y}(t_{2}))^{\top}&\rule[2.15277pt]{10.76385pt}{0.5pt}\\\
&\vdots&&&&\vdots&\\\
\rule[2.15277pt]{10.76385pt}{0.5pt}&(\bm{u}^{\prime}_{x}(t_{N_{t}}))^{\top}&\rule[2.15277pt]{10.76385pt}{0.5pt}&&\rule[2.15277pt]{10.76385pt}{0.5pt}&(\bm{u}^{\prime}_{y}(t_{N_{t}}))^{\top}&\rule[2.15277pt]{10.76385pt}{0.5pt}\end{array}\right),$
(2)

so that row $j$ contains the concatenated snapshot
$\bm{u}^{\prime}_{x}(t_{j}),\bm{u}^{\prime}_{y}(t_{j})$, i.e., the two
flattened velocity components at time step $j$. The size of the matrix
${\bm{Y}}$ is $N_{t}\times 2N$.

  * •

_Perform a spectral decomposition of the covariance matrix._ We construct the
covariance matrix ${\bm{K}}$ of the rows of ${\bm{Y}}$ and compute its
eigendecomposition:

$\displaystyle{\bm{K}}$ $\displaystyle={\bm{Y}}{\bm{Y}}^{\top},$
$\displaystyle{\bm{K}}$ $\displaystyle=\bm{A}\bm{\Lambda}\bm{A}^{T},$
$\displaystyle\bm{A}$
$\displaystyle=\left(\bm{\alpha}_{1},\ldots,\bm{\alpha}_{N_{t}}\right),$ (3)

where $\bm{\alpha}_{j}$ is the $j$th eigenvector, and the matrix
$\bm{\Lambda}$ is diagonal containing entries $\lambda_{j}$, the associated
non-negative eigenvalues of $\bm{K}$. We assume the eigenvalues are listed in
non-increasing order, $\lambda_{j}\geq\lambda_{j+1}$.

  * •

_Identify reduced-order modes and truncate._ Larger eigenvalues of ${\bm{K}}$
are directly related to the dominant characteristics of the dynamical system,
while small eigenvalues correspond to small perturbations of the dynamical
behavior. The matrix ${\bm{K}}$ has $N_{t}$ eigenvalues, and we choose the
order of the reduced-order model to be $r\ll N_{t}$ by inspecting a relative
information content $I(r)$, defined as follows

$I(r)=\frac{\sum_{i=1}^{r}\lambda_{i}}{\sum_{i=1}^{N_{t}}\lambda_{i}},$ (4)

so that $1-I(r)$ is a relative Frobenius norm error between ${\bm{K}}$ and its
rank-$r$ spectral approximation. As we will see in Section 6, $I(r)$ is close
to one for practical problems, even for very small $r$. As output of the
procedure, we can construct the following (discretized) ROM of the fluctuating
component of the dynamics

$\bm{u}^{\prime}_{*}(t_{j})\approx\sum_{i=1}^{r}(\alpha_{i})_{j}\bm{\psi}_{*,i},$
(5)

where $*\in\\{x,y\\}$, and
$(\bm{\psi}_{x,i}^{\top},\bm{\psi}_{y,i}^{\top})^{\top}\in\mathbb{R}^{2N}$ is
a vector denoting a discretized spatial function; the entries of
$\bm{\psi}_{*,i}$ correspond to the $N$ degrees of freedom in the snapshots
$\bm{u}_{\ast}$, and is a subvector of the $i$-th right-singular vector of
${\bm{Y}}$. Equivalently, it is defined as,

$\bm{\psi}_{*,i}=\frac{1}{\lambda_{i}}\bm{Y}^{\top}\bm{\alpha}_{i}=\frac{1}{\lambda_{i}}\sum_{j=1}^{N_{t}}(\alpha_{i})_{j}\bm{u}^{\prime}_{\ast}(t_{j}).$

A “standard” POD approach would next project the (assumed known) dynamical
model onto $\mathrm{span}\\{\bm{\psi}_{*,i}\\}_{i=1}^{r}$. We will proceed to
assume that such a dynamical model is unknown to us, and will instead use
machine learning models to predict dynamics.

###### Remark 1

With the training data $\alpha_{i}$ available at time steps
$t_{1},t_{2},\cdots,t_{N_{t}}$ through the above procedure, extrapolation of
the FOM dynamics or experimental measurements amounts to predicting the POD
coefficients $\alpha_{i}(t)$ for future time $t$ accurately, in our case using
machine learning models. Notice that $\alpha_{i}(t)$ is observed sequentially
and has a continuous profile, indicating the potential advantages of using
NODE for learning $\alpha_{i}(t)$, as we describe next.

##  3 Heavy-ball Neural ODEs

In this section, we briefly review NODE and HBNODE and algorithms for their
training and testing. Moreover, we provide some simple analysis of why HBNODE
is better for learning POD coefficients for reduced-order modeling leveraging
insights from the acceleration theory of the classical momentum methods.

###  3.1 Neural ODEs

NODEs chen2018neural  are a class of continuous-depth (-time) neural networks
rosenblatt1961principles ; cohen1983absolute . The continuous-time nature of
NODEs makes them particularly suitable for learning complex dynamics from
irregularly-sampled sequential data, see, e.g., chen2018neural ; latentODE ;
NEURIPS2019_21be9a4b ; massaroli2020dissecting ; norcliffe2020_sonode .
Mathematically, a NODE is formulated as the following first-order ODE:

$\frac{d{{\bm{h}}}(t)}{dt}={\bm{f}}({{\bm{h}}}(t),t,\theta),$ (6)

where ${\bm{f}}({{\bm{h}}}(t),t,\theta)\in{\mathbb{R}}^{d}$ is specified by a
neural network parameterized by $\theta$, e.g., a two-layer feed-forward
neural network. Starting from the input ${{\bm{h}}}(0)$, NODEs learn the
representation and perform prediction by solving (6) from $t=0$ to $T$ using a
numerical integrator with a given error tolerance, often with adaptive step
size solver or adaptive solver for short DORMAND198019 . Solving (6) from
$t=0$ to $T$ in a single pass with an adaptive solver requires evaluating
${\bm{f}}({{\bm{h}}}(t),t,\theta)$ at various timestamps, with computational
complexity measured by the number of function evaluations in a time-forward
sweep (“forward NFEs”) chen2018neural .

The _adjoint sensitivity method_ , or the adjoint method adjoint , is a
memory-efficient method for training NODEs through optimization of $\theta$.
We regard the output ${{\bm{h}}}(T)$ as the prediction and denote the loss
between the prediction ${{\bm{h}}}(T)$ and the ground truth as $\mathcal{L}$.
Let ${{\bm{a}}}(t):={\partial\mathcal{L}}/{\partial{{\bm{h}}}(t)}$ be the
adjoint state, then we have (see chen2018neural ; adjoint  for details)

$\frac{d\mathcal{L}}{d\theta}=\int_{0}^{T}{{\bm{a}}}(t)^{\top}\frac{\partial{\bm{f}}({{\bm{h}}}(t),t,\theta)}{\partial\theta}dt,$
(7)

with ${{\bm{a}}}(t)$ satisfying the following adjoint ODE

$\frac{d{{\bm{a}}}(t)}{dt}=-{{\bm{a}}}(t)^{\top}\frac{\partial}{\partial{{\bm{h}}}}{\bm{f}}({{\bm{h}}}(t),t,\theta),$
(8)

which is solved numerically from $t=T$ to $0$ and also requires the evaluation
of the right-hand side of (8) at various timestamps, with the number of NFEs
during this time-backward sweep (“backward NFEs”) measuring the computational
complexity.

There are several critical problems with NODEs, including (i) Given an error
tolerance, the NFEs required in a single forward pass can be excessive.
Moreover, solving the adjoint ODE (8) often requires more NFEs than solving
the forward ODE (6). (ii) In training NODEs, the adjoint state ${{\bm{a}}}(t)$
often vanishes, i.e., the norm of ${{\bm{a}}}(t)$ tends to $0$, impeding NODEs
from learning long-range dependencies lechner2020learning , resulting in poor
predictive performance.

###  3.2 Heavy-ball neural ODEs

The authors of HBNODE:2021 ; wang2021does  proposed HBNODEs and their
generalized version, named generalized HBNODEs (GHBNODEs). HBNODEs are
motivated by ideas from momentum-accelerated gradient descent polyak1964some
and they can be regarded as the continuous limit of the MomentumRNN model
MomentumRNN . Mathematically, the HBNODE is a special second-order neural ODE
of the following form

$\frac{d^{2}{\bm{h}}(t)}{dt^{2}}+\gamma\frac{d{\bm{h}}(t)}{dt}={\bm{f}}\big{(}{\bm{h}}(t),t,\theta\big{)},$
(9)

where $\gamma\geq 0$ is the damping parameter, which can be set as a tunable
or a learnable hyperparameter with positivity constraint. In the trainable
case, we adopt the one used in HBNODE:2021 , that is
$\gamma=\epsilon\cdot\text{sigmoid}(\omega)$ for a trainable
$\omega\in\mathbb{R}$ and a fixed tunable upper bound $\epsilon$, e.g.,
$\epsilon=1$. The HBNODE (9) can be rewritten as the following system of
first-order NODEs

$\frac{d{\bm{h}}(t)}{dt}={\bm{m}}(t);\quad\frac{d{\bm{m}}(t)}{dt}=-\gamma{\bm{m}}(t)+{\bm{f}}({\bm{h}}(t),t,\theta).$
(10)

####  3.2.1 Computational advantages of HBNODE vs. NODE

To show why HBNODE enjoys computational efficiency in training and testing,
let us first consider the adjoint equation of (9), which will again be solved
using adaptive numerical ODE solvers. First, the following theoretical result
HBNODE:2021  shows that the adjoint of an HBNODE is also an HBNODE.

######  Proposition 1 (Adjoint equation for HBNODE HBNODE:2021 )

The adjoint state ${\bm{a}}(t):=\partial\mathcal{L}/\partial{\bm{h}}(t)$ for
the HBNODE (9) satisfies the following HBNODE with the same damping parameter
$\gamma$ as that in (9),

$\frac{d^{2}{\bm{a}}(t)}{dt^{2}}-\gamma\frac{d{\bm{a}}(t)}{dt}={\bm{a}}(t)\frac{\partial{\bm{f}}}{\partial{\bm{h}}}({\bm{h}}(t),t,\theta).$
(11)

Notice that we solve the adjoint equation (11) from $t=T$ to $0$ via backward
propagation. By letting $\tau=T-t$ and ${\bm{b}}(\tau)={\bm{a}}(T-\tau)$, we
can rewrite (11) as follows,

$\frac{d^{2}{\bm{b}}(\tau)}{d\tau^{2}}+\gamma\frac{d{\bm{b}}(\tau)}{d\tau}={\bm{b}}(\tau)\frac{\partial{\bm{f}}}{\partial{\bm{h}}}({\bm{h}}(T-\tau),T-\tau,\theta).$
(12)

Therefore, the adjoint of the HBNODE is also an HBNODE and they have the same
damping parameter.

The above result indicates that the adjoint problem for HBNODE is of the same
type as the forward problem, accelerating backward propagation provided the
forward propagation is accelerated.

Next, we provide theoretical insights into the computational efficiency of
HBNODEs. Leveraging the acceleration theory of the heavy-ball momentum method
in taming the oscillation of the optimization trajectory by reducing the
condition number of the underlying problem. Furthermore, it is worth
mentioning that the heavy-ball method only provably accelerates strongly
convex quadratic optimization problems, while in practice, it accelerates
general optimization problems. In analogy to the acceleration theory of the
heavy-ball method, we consider the following two linearized high-dimensional
ODE systems

$\frac{d{\bm{h}}(t)}{dt}={\bm{A}}{\bm{h}}(t),$ (13)

and

$\begin{aligned} \frac{d{\bm{h}}(t)}{dt}&={\bm{m}}(t)\\\
\frac{d{\bm{m}}(t)}{dt}&=-\gamma{\bm{m}}(t)+{\bm{A}}{\bm{h}}(t)\end{aligned}\Leftrightarrow\frac{d}{dt}\begin{pmatrix}{\bm{h}}(t)\\\
{\bm{m}}(t)\end{pmatrix}=\underbrace{\begin{pmatrix}{\bf 0}&\ \ {\bm{I}}\\\
{\bm{A}}&\ \
(-\gamma{\bm{I}})\end{pmatrix}}_{:={\bm{B}}}\begin{pmatrix}{\bm{h}}(t)\\\
{\bm{m}}(t)\end{pmatrix},$ (14)

where we assume ${\bm{A}}$ is positive definite to simplify our analysis and
reveal intuition of the advantages of HBNODE over NODE. Let the eigenvalues
and eigenvectors of ${\bm{A}}$ be given by $\lambda_{i},{\bm{v}}_{i}$
respectively. Following the proof of the acceleration of heavy-ball
momentum111see http://www.math.utah.edu/~bwang/mathds/Lecture8.pdf for
details., we can show that
${|\tilde{\lambda}_{\max{}}|}/{|\tilde{\lambda}_{\min{}}|}\leq\sqrt{{|\lambda_{\max{}}|}/{|\lambda_{\min{}}|}}$,
where $\lambda_{\max{}}$ and $\lambda_{\min{}}$ are the largest and smallest
eigenvalues. Similarly, $\tilde{\lambda}_{\max}$ and $\tilde{\lambda}_{\min}$
are the largest and smallest eigenvalues, in magnitude, of ${\bm{B}}$.

Notice that the ratio ${|\lambda_{\max{}}|}/{|\lambda_{\min{}}|}$ and
${|\tilde{\lambda}_{\max{}}|}/{|\tilde{\lambda}_{\min{}}|}$ are the stiffness
of the linear ODE model (13) and the corresponding linear HBNODE counterpart
(14). Thus the heavy-ball NODE can be much less stiff than the original NODE.
If the stiffness of the ODE model is $\kappa$, using the heavy-ball model
results in stiffness of at most $\sqrt{\kappa}$, which is a substantial
reduction.

Recall that we use the adaptive step size explicit solver to solve both
forward and backward ODEs, from $t=0$ to $T$, in training NODEs and HBNODEs. A
less stiff model allows the adaptive solver to use a much large step size and
thus can significantly reduce NFEs. Moreover, Proposition 1 indicates that the
adjoint equation of an HBNODE is also an HBNODE, and therefore we reap the
computational advantages of relaxed stiffness in both forward and backward
propagation phases. Our previous analysis only considers very simple linear
ODE models. How to extend the analysis to the neural network is a very
interesting future direction. One particular idea is analyzing the NODE and
HBNODE when they are overparameterized, in which case one could leverage
neural tangent kernel theory jacot2018neural .

####  3.2.2 Generalized HBNODEs (GHBNODEs)

Compared to vanilla NODEs, high-order NODEs include HBNODEs usually suffer
from the uncontrolled aggregation of the hidden state, deteriorating model
performance at best, and blowing up training at worst. To alleviate this
issue, in HBNODE:2021  the authors propose the following generalized HBNODE

$\displaystyle\frac{d{\bm{h}}(t)}{dt}=\sigma({\bm{m}}(t));\quad\frac{d{\bm{m}}(t)}{dt}=-\gamma{\bm{m}}(t)+{\bm{f}}({\bm{h}}(t),t,\theta)-\xi{\bm{h}}(t),$
(15)

where $\sigma(\cdot)$ is a nonlinear activation, which is set as $\tanh$ by
default. The positive hyperparameters $\gamma,\xi>0$ are two tunable or
learnable hyperparameters. In the trainable case, we let
$\gamma=\epsilon\cdot\text{sigmoid}(\omega)$ as in HBNODE, and
$\xi=\text{softplus}(\chi)$ to ensure that $\gamma,\xi\geq 0$. Compared to
HBNODEs, GHBNODEs integrate two ideas to improve the neural network
architecture design: (i) Incorporating the gating mechanism $\sigma$ used in
LSTM hochreiter1997long  and GRU cho2014learning , which can suppress the
aggregation of ${\bm{m}}(t)$; (ii) Following the idea of skip connections
he2016identity , HBNODEs add the term $\xi{\bm{h}}(t)$ into the governing
equation of ${\bm{m}}(t)$, which benefits training and generalization of
GHBNODEs. It has been extensively verified that GHBNODE can indeed control the
growth of ${\bm{h}}(t)$ effectively, which significantly improve the
performance of machine learning models on various sequential learning tasks.

Another interesting result is that though the adjoint state of the GHBNODE
does not satisfy the exact heavy-ball ODE, it also significantly reduces the
backward NFEs in practice. We observe that sometimes GHBNODEs are
computationally more efficient than HBNODEs.

####  3.2.3 (G)HBNODEs learn long-range dependencies effectively

Learning long-range dependencies is crucial for the success of deep learning
for sequential data, and vanishing and exploding gradients are two bottlenecks
for training RNNs to learn long-range dependencies bengio1994learning ;
pascanu2013difficulty . The exploding gradients issue can be effectively
resolved via gradient clipping, training loss regularization, etc
pascanu2013difficulty . The vanishing gradient phenomenon in training RNNs
materializes in continuous-depth neural networks as vanishing of the adjoint
state HBNODE:2021 . In particular, we consider
${\bm{a}}(t):=\partial\mathcal{L}/\partial{\bm{h}}(t)$, and when the vanishing
gradient phenomenon occurs, ${\bm{a}}(t)$ goes to ${\bm{0}}$ quickly as $T-t$
increases, so that $d\mathcal{L}/d\theta$ in (7) will be essentially
independent of ${\bm{a}}(t)$ for larger $T-t$. We have the following
expressions for the adjoint states of the NODE and HBNODE (see HBNODE:2021
for details):

  * •

For NODE, we have

$\frac{\partial\mathcal{L}}{\partial{\bm{h}}_{t}}=\frac{\partial\mathcal{L}}{\partial{\bm{h}}_{T}}\frac{\partial{\bm{h}}_{T}}{\partial{\bm{h}}_{t}}=\frac{\partial\mathcal{L}}{\partial{\bm{h}}_{T}}\exp\Bigg{\\{}-\int_{T}^{t}\frac{\partial{\bm{f}}}{\partial{\bm{h}}}({\bm{h}}(s),s,\theta)ds\Bigg{\\}}.$
(16)

  * •

For GHBNODE222HBNODE can be seen as a special GHBNODE with $\xi=0$ and
$\sigma$ be the identity map., we have

$\displaystyle\begin{bmatrix}\frac{\partial\mathcal{L}}{\partial{\bm{h}}_{t}}&\frac{\partial\mathcal{L}}{\partial{\bm{m}}_{t}}\end{bmatrix}$
$\displaystyle=\begin{bmatrix}\frac{\partial\mathcal{L}}{\partial{\bm{h}}_{T}}&\frac{\partial\mathcal{L}}{\partial{\bm{m}}_{T}}\end{bmatrix}\begin{bmatrix}\frac{\partial{\bm{h}}_{T}}{\partial{\bm{h}}_{t}}&\frac{\partial{\bm{h}}_{T}}{\partial{\bm{m}}_{t}}\\\
\frac{\partial{\bm{m}}_{T}}{\partial{\bm{h}}_{t}}&\frac{\partial{\bm{m}}_{T}}{\partial{\bm{m}}_{t}}\\\
\end{bmatrix}$ (17)
$\displaystyle=\begin{bmatrix}\frac{\partial\mathcal{L}}{\partial{\bm{h}}_{T}}\
\frac{\partial\mathcal{L}}{\partial{\bm{m}}_{T}}\end{bmatrix}\exp\Bigg{\\{}-\underbrace{\int_{T}^{t}\begin{bmatrix}{\bf
0}&\frac{\partial\sigma}{\partial{\bm{m}}}\\\
\big{(}\frac{\partial{\bm{f}}}{\partial{\bm{h}}}-\xi{\bm{I}}\big{)}&-\gamma{\bm{I}}\end{bmatrix}ds}_{:={\bm{M}}}\Bigg{\\}}.$

For the matrix ${\bm{M}}$, we have the following useful property about its
spectrum.

######  Proposition 2 (HBNODE:2021 )

The eigenvalues of $-{\bm{M}}$ can be paired so that the sum of each pair
equals $(t-T)\gamma$.

Following the argument in HBNODE:2021 , Proposition 2 can be used to show that
the adjoint state of NODE in (16) may vanish when $T-t$ is large, but the
adjoint state of (G)HBNODEs in (17) will not vanish. This property supports
the claim that HBNODEs benefit in learning long-range dependencies, which in
turn further boosts the accuracy in learning POD of complex dynamical systems.

In Section 6, we will validate the above theoretical merits of HBNODEs over
NODEs using the benchmark problems listed in Section 4 below.

##  4 Benchmarks and Data Preparation

In this section, we will present some details of the three benchmark physical
models — VKS, KPP, and Euler equations — used for validating the efficacy of
learning POD with HBNODEs. Our training data for the KPP and Euler equations
are generated by solving FOMs; the training data for the VKS dataset is
adopted from a publicly available dataset.

###  4.1 VKS model

The von Kármán vortex street (VKS) is a fluid dynamics phenomenon where
vortices appear in a periodic fashion in the wake of flow past a blunt object,
frequently a cylinder. A very small Reynolds number results in a laminar
smooth flow past the cylinder, and very large Reynolds numbers result in a
turbulent flow. In an appropriate middle regime, the VKS phenomenon appears
and can be simulated. The associated dynamical model is the two-dimensional
Navier-Stokes equations; with ${{\bm{u}}}=(u_{x},u_{y})$ the fluid velocity,
these equations read,

$\displaystyle\frac{\partial\bm{u}}{\partial
t}=\nu\nabla^{2}\bm{u}-(\bm{u}\cdot\nabla)\bm{u}-\frac{1}{\rho_{0}}\nabla p,$

where $\rho_{0}$ is the spatially uniform pressure, and $\nu$ is the kinematic
viscosity, which is inversely related to the Reynolds number. Our experimental
setup concerns flow past a cylinder in two spatial dimensions with conditions
that result in steady-state VKS flow after an initial transient period. We
follow the experimental setting used in rom_node  to acquire simulation data.

###  4.2 KPP model

The Kurganov-Petrova-Popov (KPP) model is a scalar, two-dimensional
conservation law, first proposed in kurganov_adaptive_2007 . This system is
difficult to simulate since it features a non-convex flux, and is given by,

$\displaystyle\frac{\partial u}{\partial t}+\nabla\cdot\bm{f}(u)=0,t>0,\qquad
x\in[-2,2],\hskip 5.0pty\in\left[-\frac{5}{2},\frac{3}{2}\right],$

where $\bm{f}(u)=\left(\sin u,\cos u\right)^{T}$ and
$\nabla:=\left(\frac{\partial}{\partial x},\frac{\partial}{\partial
y}\right)$. Our setup mirrors that in kurganov_adaptive_2007 , so that we use
the following initial data

$\displaystyle u(x,y,0)$
$\displaystyle=\left\\{\begin{array}[]{rl}\frac{14\pi}{4},&x^{2}+y^{2}<1,\\\
\frac{\pi}{4},&\mathrm{else}\end{array}\right.$

We employ a finite volume scheme utilizing a Lax-Friedrichs flux with a 5th-
order WENO reconstruction over the two-dimensional rectangular domain with a
Cartesian mesh up to time $T=10$. The simulation uses a tensorial grid with
$N_{x}=50$, $N_{y}=50$ (corresponding to $N=N_{x}N_{y}=2500$ total spatial
degrees of freedom), and $N_{t}=1250$.

###  4.3 Euler equations for fluids modeling

The one-dimensional Euler equations of gas dynamics are a system of
conservation laws. We consider the simulation of a parameterized shock-entropy
problem from this differential equation, whose setup is given by,

$\displaystyle\frac{\partial\bm{u}}{\partial
t}+\frac{\partial\bm{f}(\bm{u})}{\partial x}$ $\displaystyle=0,\quad
t>0,x\in[-5,5],$ $\displaystyle\mbox{with}\ \bm{f}(u)$
$\displaystyle=\left(\begin{array}[]{c}\rho u\\\ \rho u^{2}+p\\\
(E+p)u\end{array}\right)^{T},$

where $\bm{u}:=(\rho\ \rho u\ E)^{\top}\in{\mathbb{R}}^{3}$ is the unknown
with $(\rho,u,p,E)$ denoting the gas density, velocity, pressure, and energy,
respectively. The system is closed via the following relationship between $E$
and $p$:

$\displaystyle p=(\gamma-1)\left(E-\frac{1}{2}\rho u^{2}\right),$

where $\gamma$ is the heat capacity ratio, a gas-dependent constant.333This
$\gamma$ is distinct from the $\gamma$ discussed in Section 3.2. We take
boundary conditions at $x=\pm 5$ as those given by the initial data. The
shock-entropy problem features smoothly oscillating as well as discontinuous
features.

We again employ a finite volume scheme to solve the Euler equations, using a
Harten-Lax-van Leer (HLL) flux, which is an approximated Riemann solver
harten_upstream_1983 . Our simulations integrate up to terminal time $T=1.8$,
with a uniform grid having $N=1000$ degrees of freedom in the scalar spatial
variable $x$.

This last example differs from the previous two in that we consider this a
parametric equation, where
$\bm{\eta}=(\eta_{u},\eta_{\rho})\in{\mathbb{R}}^{2}$ is a parameter for the
initial conditions. We initialize the dynamics using the parameter $\bm{\eta}$
as follows, where $\eta_{u}$ varies on the interval $[2,3]$ and $\eta_{\rho}$
varies on the interval from $[3,4]$. The parametric initial data
$(u(x,0),\rho(x,0),p(x,0))=(u_{0},\rho_{0},p_{0})$ are given by

$u_{0}=\begin{cases}\eta_{u}&x<-4\\\
0&\text{else}\end{cases},\quad\rho_{0}=\begin{cases}\eta_{\rho}&x<-4\\\
1+0.2\sin(\pi x)&\text{else}\end{cases},\quad
p_{0}=\begin{cases}\frac{31}{3}&x<-4\\\ 1&\text{else}\end{cases}$

We generate training data by gathering an ensemble of trajectories for the
above problem over a grid of $\bm{\eta}$ values and attempt to learn dynamics
on unseen values of $\bm{\eta}$. Thus, in this example we not only seek to
predict to future times, but also trajectories on parameter values not in the
training set.

##  5 Learning Pipeline

In this section, we describe the detailed pipeline of using deep learning for
reduced-order modeling accompanied by the baseline ROMs.

###  5.1 Learning-based reduced-order modeling

Our machine learning-based reduced-order modeling framework is flexible for
machine learning model selection, e.g., using either HBNODE or NODE as shown
in Fig. 2 and Fig. 3, respectively. In our learning-based reduced-order
modeling framework, we first apply POD outlined in Section 2 on the training
data to extract (discretized) temporal coefficients $\bm{\alpha}(t)$’s and the
eigenmodes $\bm{\psi}({\bm{x}})$’s following (5). Next, we will use machine
learning models to predict future dynamics $u({\bm{x}},t)$ leveraging these
coefficients and modes. In particular, the main task is an extrapolation of
the temporal coefficients $\bm{\alpha}(t)$’s using NODEs or HBNODEs.

To predict future values of the POD data, we consider two different machine
learning architectures, shown in Fig. 2 and Fig. 3, respectively. The first
architecture is a one-to-one architecture that predicts the value at $t_{k+1}$
based on the data at $t_{k}$. The second architecture is a sequence-to-
sequence architecture that uses sequence data points to predict the following
sequence of data points. The overlap in the sequence prediction can be
adjusted so that the predicted sequence is entirely new or that only one new
data point is predicted.

The first architecture under our study is adapted from rom_node , which was
originally used to compare the performance of NODE and LSTM in model
reduction. We replace the vanilla NODE used in rom_node  with the HBNODE, and
we depict the modified architecture in Fig. 2. Compared to the pipeline used
in rom_node , after the RNN encoding of the temporal coefficients we have to
sample both ${\bm{h}}$ and ${\bm{m}}$ to accommodate learning using HBNODE. In
contrast, the vanilla NODE used in rom_node  only needs to sample the state
${\bm{h}}$. The above encoding and sampling procedure is accomplished via a
variational autoencoder kingma2013auto1 .

Figure 2: The pipeline of predicting the temporal coefficients for a single
step forward using HBNODE leveraging a variational autoencoder. We first use a
RNN encoder to encode the input data and then sample the states ${\bm{h}}$ and
${\bm{m}}$ and evolve them using an HBNODE. Finally, we apply an RNN decoder
to the final representation to get the prediction.

We plot the second architecture in Fig. 3, where the vanilla NODE can be
replaced with (generalized) HBNODE. For the second architecture, i.e., the
sequence-to-sequence architecture, takes a sequence of length $n$ inputs and
predicts a sequence of outputs, we encode the input sequence
$\\{\alpha_{i}(t_{j})\\}_{j=0}^{n-1}$ into the latent sequence
$\\{z_{i}(t_{j})\\}_{j=0}^{n-1}$ using an RNN encoder, then we use NODE or
HBNODE to evolve the latent sequence to get the desired representation,
followed by an RNN decoder to get the final long-term prediction
$\\{\alpha_{i}(t_{j})\\}_{j=n}^{N}$.

Figure 3: The pipeline of predicting the temporal coefficients for multi-steps
ahead. First, we encode the input sequence
$\\{\alpha_{i}(t_{j})\\}_{j=0}^{n-1}$ using a RNN encoder to obtain the latent
sequence $\\{z_{i}(t_{j})\\}_{j=0}^{n-1}$. Second, we use a NODE to learn a
“good” representation of the input sequence by evolving the latent sequence
$\\{z_{i}(t_{j})\\}_{j=0}^{n-1}$. Third, we apply a RNN decoder to the “good”
representation to get the final prediction. Notice that NODE can be replaced
with (generalized) HBNODE, in which case we need to obtain another sequence of
momentum states from the RNN encoder.

###  5.2 A baseline comparison: Dynamic Mode Decomposition (DMD)

We employ DMD as another baseline model reduction method to demonstrate the
effectiveness of learning-based model reduction using HBNODEs. In this part,
we briefly review the idea of DMD for reduced-order modeling. To compare DMD
to the learning-based reduced-order modeling using HBNODE, we consider only
modeling the fluctuating components ${\bm{u}}^{\prime}$ of the snapshots, see
(1). The predictions of DMD are generated by a linear operator $\bm{A}$
corresponding to a linear difference equation
${\bm{u}}^{\prime}_{k+1}={\bm{A}}{\bm{u}}^{\prime}_{k}$, where $\bm{A}$ must
be learned. In DMD, dominant eigenvalues and eigenvectors of ${\bm{A}}$ are
computed via the singular value decomposition (SVD). Although the true
underlying dynamics may be nonlinear, the Koopman operator formalism concludes
that a lifted version of the dynamics is indeed linear. For nonlinear
problems, DMD attempts to learn these lifted linear dynamics.

Let ${\bm{U}}^{(k)}$ be the snapshot matrix for the time interval
$t_{0},\ldots,t_{\text{train}-1}$ and ${\bm{U}}^{(k+1)}$ be the snapshot
matrix for the time interval $t_{1},\ldots,t_{\text{train}}$, i.e., column $j$
of $\bm{U}^{(k)}$ corresponding to time snapshot $t_{j-1}$. In particular, let
${\bm{U}}^{(k+1)}\approx{\bm{A}}{\bm{U}}^{(k)}$ where ${\bm{U}}^{(k)}$ is
given by the SVD ${\bm{U}}^{(k)}={\bm{X}}\Sigma{\bm{V}}^{*}$. We further
denote $\tilde{\bm{X}}$, $\tilde{\bm{V}}$, and $\tilde{\bm{\Sigma}}$ as the
rank-$r$ truncation of $\bm{X}$, $\bm{V}$, and $\bm{\Sigma}$, respectively.
Then we may compute an approximation $\tilde{{\bm{A}}}$ directly from
${\bm{U}}^{(k+1)}$ by the following,

$\tilde{{\bm{A}}}=\tilde{{\bm{X}}}^{*}{\bm{U}}^{(k+1)}\tilde{{\bm{V}}}\tilde{\bm{\Sigma}}^{-1}$
(18)

The reduced matrix $\tilde{{\bm{A}}}$ is composed of the dominant $r$
eigenvalues $\lambda_{1},\ldots,\lambda_{r}$ and eigenvectors
$\bm{\phi}_{1},\ldots,\bm{\phi}_{r}$. These eigenvectors are also known as the
DMD modes. Given training data on the training interval
$t_{0},\ldots,t_{\text{train}}$, the matrix $\tilde{{\bm{A}}}$ is formulated
by partitioning the snapshot matrix ${\bm{u}}^{\prime}$ into two time
intervals. Validation data on the interval
$t_{\text{train}+1},\ldots,t_{\text{valid}}$ is generated by solving
$u^{\prime}(t_{\text{train}+k})=\tilde{{\bm{A}}}^{k}u^{\prime}(t_{\text{train}})$.
We depict DMD-based reduced-order modeling in Fig. 4. More details of DMD can
be found at e.g., schmid_2010 .

Figure 4: DMD pipeline: The data is pre-processed by subtracting the mean to
capture the fluctuating components of the data. In addition, the following
lifts $\\{\cos({\bm{x}}),\sin({\bm{x}}),{\bm{x}}^{2},{\bm{x}}^{3}\\}$ were
applied to the data, and then vectorized along the snap-shot axis. We then
generate the full-order model according to the spectral decomposition of the
linear transformation between the two snapshot matrices at subsequent time
intervals. To reduce the order only the dominant $r$ eigenmodes are selected,
resulting in a representation
$u^{\prime}({\bm{x}},t)=\sum^{r}_{i=1}\phi_{i}\lambda_{i}^{k-1}b_{i}$ for the
lifted data $u^{\prime}({\bm{x}},t)$.

##  6 Experimental Results

In each experiment below, we contrast the performance of HBNODE-based ROM to
two baseline ROMs, namely, NODE-based and DMD-based ROMs. We observe
consistently improved predictive performance of HBNODE over baseline ROMs. We
interpret the improved performance using HBNODEs by inspecting the stiffness
and adjoint state of HBNODEs, confirming the theoretical results. Animated
comparisons of the data reconstructions can be found at github-animation .

|
---|---
(a) Transient state $t=50$ | (b) Steady-state $t=150$
Figure 5: Comparison of transient and steady-state phases of the VKS dataset.
The steady-state phase contains quasi-periodic solutions conducive to machine
learning. The transient phase does not contain such well-behaved dynamics.

###  6.1 Transient and steady-state VKS

The VKS dataset is obtained by simulating the FOM in Section 4.1 on the time
interval $[0,400]$, containing two different regimes. When $t<100$, the
dynamics lie in the transient state and approach the steady state as $t$
increases; while the dynamics maintain a steady state when $t\geq 100$, as
shown in Fig. 5.

---
Figure 6: The steady-state POD modes are highly oscillatory with quasi-periodic patterns. |
---|---
(a) POD decay | (b) DMD decay
Figure 7: Comparison of the relative information decay for POD and DMD over
the steady-state VKS. The POD modes decay far more rapidly than those of DMD.
Therefore, one can expect better results with a smaller order of ROM using POD
than DMD.

##### ROMs for steady-state dynamics.

We contrast different ROMs for simulating VKS flow in the steady-state regime.
In particular, both the DMD and POD training is taken over the time interval
from $t=100$ to $400$. The POD modes for the steady-state flow oscillate
quasi-periodically, see Fig. 6, and the relative information content $I(r)$ in
(4) decays rapidly in $r$. The POD relative information content for $8$
leading modes is $\sim 99\%$, as illustrated in Fig. 7 (a). In contrast, the
lifted DMD model, using the lifts
$\\{\cos({\bm{x}}),\sin({\bm{x}}),{\bm{x}}^{2},{\bm{x}}^{3}\\}$ with
${\bm{x}}=(x,y)$, requires 24 modes to achieve $\sim 99\%$ relative
information content, shown in Fig. 7 (b). The quasi-periodic nature of the POD
modes indicates that a model with high training accuracy will continue to
perform well on the validation data.

##### ROMs for transient to steady state dynamics.

ROMs behave very differently over the entire time interval from $t=0$ to
$400$. The POD modes do not oscillate over the entire interval but only over
the steady-state phase. For both POD and DMD, the relative information decays
much slower. The POD relative information decreases to $\sim 96\%$ for the
dominant $8$ modes. While for the dominant $24$ lifted DMD modes, the relative
information is reduced to $\sim 94\%$ over the full dynamics. This suggests
that a machine learning-based ROM which is able to train on the transient
phase to predict the steady-state phase accurately captures the intrinsic
patterns of the underlying dynamics.

Hyperparameter |  Value
---|---
Latent dimension | 6
Layers encoder | 4
Units encoder | 10
Layers ODE | 12
Units decoder | 41
Layers decoder | 4
Learning rate | .00153
Epochs | 2000
Table 1: The hyperparameters for the VAE architecture — shown in Fig. 2 — for
NODE and HBNODE-based ROMs. The parameters are tuned to the best NODE
specification.

##### Learning steady-state dynamics.

In this task, we train the pipeline shown in Fig. 2 for single-input-single-
output dynamics prediction. Following the baseline in rom_node , we train over
the steady-state dynamics starting from $t=100$ using the dominant $8$ POD
modes. The training data consists of the POD modes from $t=100$ to $174$, and
the training labels consist of the POD modes from $t=101$ to $175$. The
validation data consists of the POD modes from $t=175$ to $199$, with the
objective to predict the POD modes at time steps from $t=176$ to $200$. We use
the mean squared error to measure the loss between the labeled data and the
predictions. We utilized an AdamW optimizer to train the network based on this
loss criteria. For the black-box integration method, we selected DOPRI-5
DORMAND198019  with a relative tolerance of $1\mathrm{e}{-8}$. The model’s
hyperparameters are tuned to best the NODE as outlined in rom_node  and
restated in Table 1.

Hyper-parameter |  Value
---|---
Layers | 12
Hidden layers | 64
Sequence length | 9
Learning rate | .001
Epochs | 500
Table 2: Hyperparameters of NODE and HBNODE for learning ROMs from transient
to steady-state VKS dynamics.

##### Learning transient to steady-state dynamics.

In this task, we train the pipeline outlined in Fig. 3 for multi-input-single-
output dynamics prediction. The objective of this task is to capture the phase
transition at $t=100$. The data consisted of the dominant $8$ POD modes for
the time interval from $t=0$ to $400$. The data was sequenced in a multi-
input-single-output structure so that $9$ preceding time steps were used to
predict the $10$-th time step. The training data consists of the POD modes for
the transient time interval from $t=0$ to $79$. The training labels consisted
of the POD modes from $t=10$ to $80$. The validation data utilizes the POD
modes from steady-state time interval $t=80$ to $119$, and the validation
labels consist of data from $t=90$ to $120$. The other experimental settings
follow the above single-input-single-output scenario. The model’s
hyperparameters are the same for NODE and HBNODE components and are given in
Table 2.

|
---|---
(a) Training loss | (b) Validation loss
Figure 8: Contrasting NODE and HBNODE-based ROMs for learning steady-state VKS
dynamics. HBNODE outperforms NODE in both training and validation loss.

####  6.1.1 Results and comparison to existing ROMs

##### Results of learning steady-state dynamics.

We contrast HBNODE and NODE-based ROMs in Fig. 8 and Fig. 9. Figure 8 shows
that HBNODE-based ROM not only achieves remarkably smaller training loss but
also significantly smaller validation loss than NODE-based ROM. In terms of
the predictive performance, we see that HBNODE performs better at capturing
several of the peaks of the oscillatory modes as shown in Fig. 9.

|
---|---
(a) VAE-NODE dominant mode | (b) VAE-HBNODE dominant mode
Figure 9: Comparison of modes reconstruction for NODE and HBNODE-based ROMs
for learning steady-state VKS dynamics. HBNODE captures the peaks of the
dominant POD modes better than NODE. Before and after the vertical blue line
stands for training and validation, respectively.

##### Results of learning transient to steady-state dynamics.

|
---|---
(a) Training loss | (b) Validation loss
Figure 10: Contrasting training and validation loss of NODE and HBNODE for
learning ROMs of transient to steady-state VKS dynamics. The progress made by
the NODE is significantly reduced by the vanishing gradient. HBNODE has a much
slower decaying gradient and is able to continue to make progress in both
training and validation sets.

Compared to learning steady-state VKS dynamics, HBNODE achieves more
significant performance gain over NODE for learning transient to steady-state
dynamics in terms of training and validation loss, as shown in Fig. 10. Since
we are doing sequential learning, one interpretation of the improvement in
learning dynamics is the effective learning of long-range dependencies.
Indeed, the criterion of learning long-range dependencies has been widely used
in measuring the efficacy of sequential learning models pascanu2013difficulty
; HBNODE:2021 . In NODE and HBNODE, the effectiveness of learning long-range
dependencies can be measured by whether the adjoint state vanishes quickly or
not. We visualize the evolution of the magnitude of the adjoint states of NODE
and HBNODE in Fig. 11, which support the theoretical result in Section 3.2.3.
In particular, we see that the adjoint state of NODE vanishes much more
rapidly than that of HBNODE as $T-t$ increases. A more detailed connection
between the adjoint state and learning long-range dependencies is provided in
HBNODE:2021 .

|
---|---
(a) NODE adjoint state | (b) HBNODE adjoint state
Figure 11: Comparison of the adjoint states for the NODE and HBNODE in learning multi-input-single-output. The NODE adjoint state vanishes substantially faster than HBNODE. |
---|---
(a) VAE-NODE modes | (b) VAE-HBNODE modes
Figure 12: Comparison of modes reconstruction for NODE and HBNODE-based ROMs
for learning multi-input-single-output task. NODE is able to reliably learn
the first two dominant modes but is unable to capture the steady-state
dynamics, having a much larger frequency and introducing substantial lag into
the oscillation frequency. HBNODE is able to more reliably capture the steady-
state dynamics with slightly larger frequency and developing lag much later
than the NODE component. Before and after the vertical blue line stands for
training and validation, respectively.

In terms of the predictive performance, as shown in Fig. 12, the HBNODE
predictor captures the peaks of the oscillatory dynamics better than NODEs,
especially in the first two modes $\alpha_{1}$ and $\alpha_{2}$. Moreover, the
prediction error using NODE is much larger than that of HBNODE, and the
prediction error amplifies as the prediction time goes, in particular, for
modes $\alpha_{3}$ and $\alpha_{4}$.

Another primary advantage of HBNODE over NODE-based ROMs lies in computational
efficiency, which is theoretically supported by the discussion in Section
3.2.1. As shown in Fig. 13 (a), the forward NFE required in each forward pass
by HBNODE is consistent smaller than that of NODE. We also monitor the
stiffness of both NODE and HBNODE during the learning process, and Fig. 13 (b)
shows that the stiffness of NODE oscillates and maintains much larger than
HBNODE.

|
---|---
(a) Forward NFE comparison | (b) Stiffness comparison
Figure 13: Comparison of NFEs and stiffness for the NODE and HBNODE in learning transient to steady-state VKS dynamics. NODE requires more NFEs in each forward pass than HBNODE as the NODE is much stiffer than HBNODE. The stiffness of NODE varies sharply as training goes on, while the stiffness of HBNODE decays during the training.  |
---|---
(a) POD mode decay | (b) DMD mode decay
Figure 14: The decay of relative information content for POD and DMD of the
KPP dataset. The rapid decay in the modes indicates the problem is suitable
for POD and DMD.

###  6.2 KPP model

We obtain the KPP dataset by simulating the FOM presented in Section 4.2 for
$1000$ timesteps ($N_{t}=1000$). The KPP model is well-suited for reduced-
order modeling due to the rapidly decaying eigenvalues in both POD and DMD,
seeing Fig. 14. However, we found in our experiments that it is particularly
difficult to capture the dynamics using machine learning architectures due to
the slow decaying ROM dynamics depicted in Fig. 15.

---
KPP POD modes
Figure 15: The dominant POD modes for the KPP data. The oscillations of the
modes decay rapidly and have a very low frequency, making it challenging for
learning compared to the quasi-periodic oscillations of the VKS dataset.

In our experiments, we note that the non-lifted DMD continuously deforms the
center of mass in a way that defies the physical constraints of the system. A
comparison of lifted and non-lifted DMD predictions and POD predictions are
shown in github-animation . To lift DMD, we utilized the lifting functions
$\\{\cos({\bm{x}}),\sin({\bm{x}}),{\bm{x}}^{2},{\bm{x}}^{3}\\}$ with
${\bm{x}}=(x,y)$. Figure 14 shows that the POD modes decay faster than the
lifted DMD modes. The dominant $24$ DMD modes correspond to $97\%$ of the
relative information content; in contrast, the dominant $8$ POD modes
correspond to $99\%$ of the relative information value.

We train the pipeline depicted in Fig. 3 for learning multi-input-single-
output dynamics. The data is constructed from the $8$ dominant POD modes on
the time interval from $t=0$ to $1000$. The data is sequenced so that every
$4$ preceding time step is used to predict the $5$-th time step. The training
data consists of the POD modes from $t=0$ to $799$ and the training labels
consist of POD modes from $t=5$ to $800$. The validation data utilizes data
from $t=800$ to $999$ and the validation labels consist of data from $t=805$
to $1000$. The model’s hyperparameters are the same for NODE and HBNODE
components and are given in Table 3.

Hyper-parameter |  Value
---|---
Layers | 2
Hidden layers | 64
Sequence length | 4
Learning rate | .01
Epochs | 500
Table 3: The hyperparameters of NODE and NODE for the learning ROMs of KPP
model.

####  6.2.1 Results and comparison to existing ROMs

|
---|---
(a) NODE mode recapture | (b) HBNODE mode recapture
Figure 16: Comparison of the KPP modes prediction using NODE and HBNODE. HBNODE is better at predicting the dynamics for the validation set than NODE. Before and after the vertical blue line stands for training and validation, respectively. |
---|---
(a) Training loss | (b) Validation loss
Figure 17: Comparison of the training and validation loss of NODE and HBNODE
for learning ROMs for the KPP model. NODE is unable to make progress due to a
rapidly vanishing gradient, impeding learning long-range dependencies.

We compare the prediction of NODE and HBNODE against ground truth in Fig. 16,
and we see that HBNODE performs remarkably better than NODE in predicting the
dynamics. In particular, HBNODE is able to properly capture the oscillation
dynamics of the modes, unlike NODE. Figure 17 shows that HBNODE has a much
smaller training and validation loss than that of NODE.

###  6.3 Euler equations for fluids modeling

Figure 18: The Euler Equations data at time step $t=50$ with two different
initial parameters, $\eta_{0}$ and $\eta_{50}$. Variations in the parameter
$\eta$ produce widely varying dynamics and as a result varying POD modes. The
objective of this task is to predict the dynamics for unseen parameters $\eta$
using a set of training parameters $\eta_{\mathrm{train}}$.

We further consider learning reduced-order models for simulating the Euler
equations, where the dataset is obtained by simulating the full-order model
presented in Section 4.3 with a discrete ensemble of parameters
$\bm{\eta}_{1},\ldots,\bm{\eta}_{M}$ with $M=100$, over 180 timesteps. Two
different parameters $\bm{\eta}$ values can produce rather different dynamics,
as evidenced in Fig. 18. The Euler equations data is unique in the sense that
it may be segmented based on these initial conditions. The ROM is generated by
taking the dominant $8$ POD modes for each parameter $\bm{\eta}_{i}$ on the
time interval from $t=0$ to $180$. This data is shuffled randomly among the
initial parameter $\bm{\eta}_{i}$ to no longer increase sequentially. The
average relative information content across all $100$ values of $\bm{\eta}$ is
$\sim 95\%$.

In this task, we train the machine learning pipeline shown in Fig. 3 for
learning multi-input-multi-output dynamics. The training dataset comprises the
dominant $8$ POD modes for each of the training parameters among
$\eta_{1},\ldots,\eta_{\mathrm{train}}$. We use $90$ of the $100$ parameters,
$\eta_{1},\ldots,\eta_{90}$, for training and the rest for validation. The
training input consists of the dominant $8$ POD modes for each of the $90$
training parameters on the time interval from $t=0$ to $150$. The training
labels consist of the dominant $8$ POD modes for each of the $90$ training
parameters time steps from $t=151$ to $180$. The validation dataset is
composed of the validation parameters $\eta_{91},\ldots,\eta_{100}$. The
validation input and labels are segmented using the same intervals as the
training data.

The model uses a GHBNODE component with a hyperbolic tangent activation
function. All other experimental settings are the same as in the KPP dataset,
and the tuned hyper-parameters are listed in table 4. The NODE and HBNODE
models are trained and validated over the same data shuffling.

Hyper-parameter |  Value
---|---
Layers | 6
Hidden layers | 16
Learning rate | .01
Epochs | 100
Table 4: The hyperparameters of NODE and NODE for the learning ROMs of Euler
equations.

####  6.3.1 Results and comparison to existing models

We compare the prediction of NODE and GHBNODE for a randomly selected
parameter $\eta_{\mathrm{train}}$ from the training set and a randomly
selected parameter $\eta_{\mathrm{valid}}$ from the validation set. The modes
for the training parameter $\eta_{\mathrm{train}}$ are shown in Fig. 19, and
the modes for the validation parameter $\eta_{\mathrm{valid}}$ are shown in
Fig. 20. We observe that the POD modes for the parameter
$\eta_{\mathrm{train}}$ in Fig. 19 differ from those for
$\eta_{\mathrm{valid}}$ primarily in amplitude rather than shape. As a result,
a poor prediction model will have a sudden discontinuity between the input and
prediction values. The transition point between the input and the prediction
is indicated in Fig. 19 and Fig. 20 by the vertical blue line.

Figure 19 shows the dominant $4$ POD modes for $\eta_{\mathrm{train}}$ from
the training set. The predictive capabilities of GHBNODE significantly
outperform the NODE model. In particular, for $\alpha_{2}$ and $\alpha_{3}$,
we observe that the NODE has a large jump discontinuity at the transition
between the input and prediction. The GHBNODE modes are smoother in the
transition region, which indicates the ability of the GHBNODE model to
distinguish between separate parameters.

|
---|---
(a) NODE Training Mode Recapture | (b) GHBNODE Training Mode Recapture
Figure 19: For a randomly selected initial parameter $\eta_{\mathrm{train}}$,
we plot the dominant four modes and their predicted data using NODE and
GHBNODE. The blue vertical line separates the input and output data. The
ground truth is in black, while the prediction data is in dashed red. GHBNODE
can learn the dynamics of $\alpha_{2}$ and $\alpha_{1}$ for the parameterized
data significantly better than NODE. This is evidenced by the fact that NODE
predicts the inflection point of $\alpha_{2}$ too early.

In Figure 20, we observe the same characteristics for NODE and GHBNODE. NODE
is unable to accurately predict the output for the parameter
$\eta_{\mathrm{valid}}$ from the validation set. In particular, for
$\alpha_{3}$ of the validation parameter, the NODE prediction is even less
smooth than that of the training parameter shown in Fig. 19. In this
experiment, we observe that NODE is unable to distinguish data with varying
parameters as accurately as GHBNODE.

|
---|---
(a) NODE Validation Mode Recapture | (b) GHBNODE Validation Mode Recapture
Figure 20: For the randomly selected initial parameter
$\eta_{\mathrm{valid}}$, we plot the dominant four modes and their predicted
data using NODE and GHBNODE. The blue vertical line separates the input and
output data. The ground truth is in black, while the prediction data is in
dashed red. HBNODE can predict the data for unseen parameterizations much
better than NODE.

##  7 Concluding Remarks

This paper employs the recently developed HBNODEs and their generalization for
learning POD coefficients for model reduction. We analyze through simple
linearized models and empirically verify the advantages of HBNODEs over
existing NODEs. In particular, HBNODEs enjoy the following advantages that
imply practical benefits for learning POD-based ROMs, including 1) The deep
learning model is continuous-depth, providing flexibility in learning
irregularly-sampled time series and faithful to the continuous profiling of
the underlying physical models. 2) Both the forward and adjoint ODEs of
HBNODEs are of the heavy-ball style, accelerating both training and testing of
the machine learning procedure. And 3) HBNODEs can learn long-range
dependencies effectively, capturing intrinsic patterns from data. There are
numerous avenues for future works, and two particular interesting directions
in our mind are 1) Improving HBNODEs, particularly replacing the fine-tuned or
learned damping parameter with an adaptive one that are motivated by certain
optimization algorithms with adaptive momentum wang2020scheduled ;
sun2021training ; wang2020stochastic , and 2) Applying HBNODE-based ROMs to
model reduction arising from scientific challenges, especially when we do not
have the ground truth governing equation of the dynamical systems.

##  8 Data Availability Statement

All data and code related to this paper are available at
https://github.com/JustinBakerMath/pod_hbnode/.

##  9 Acknowledgement

This material is based on research sponsored by NSF grants DMS-1848508,
DMS-1924935, DMS-1952339, DMS-2110145, DMS-2111117, DMS-2152762, and
DMS-2208361, DOE grant DE-SC0021142 and DE-SC0023490, and AFOSR
FA9550-20-1-0338. We also acknowledge support from a seed grant from the
College of Science at the University of Utah.

## References

  * (1) A. Antoulas.  Approximation of Large-Scale Dynamical Systems.  Advances in Design and Control. Society for Industrial and Applied Mathematics, January 2005.
  * (2) Athanasios C. Antoulas, Christopher A. Beattie, and Serkan Gugercin.  Interpolatory Model Reduction of Large-Scale Dynamical Systems.  In Javad Mohammadpour and Karolos M. Grigoriadis, editors, Efficient Modeling and Control of Large-Scale Systems, pages 3–58. Springer US, Boston, MA, 2010.
  * (3) Justin Baker, Elena Cherkaev, Akil Narayan, and Bao Wang.  Learning pod of complex dynamics usingheavy-ball neural odes: Animations.  https://www.github.com/JustinBakerMath/pod_hbnode/blob/master/README.md#animations.
  * (4) Yoshua Bengio, Patrice Simard, and Paolo Frasconi.  Learning long-term dependencies with gradient descent is difficult.  IEEE Transactions on Neural Networks, 5(2):157–166, 1994.
  * (5) P. Benner, S. Gugercin, and K. Willcox.  A Survey of Projection-Based Model Reduction Methods for Parametric Dynamical Systems.  SIAM Review, 57(4):483–531, January 2015.
  * (6) Peter Benner, Serkan Gugercin, and Karen Willcox.  A survey of projection-based model reduction methods for parametric dynamical systems.  SIAM review, 57(4):483–531, 2015.
  * (7) Gal Berkooz, Philip Holmes, and John L Lumley.  The proper orthogonal decomposition in the analysis of turbulent flows.  Annual review of fluid mechanics, 25(1):539–575, 1993.
  * (8) L. Bittner.  L. S. Pontryagin, V. G. Boltyanskii, R. V. Gamkrelidze, E. F. Mishechenko, the mathematical theory of optimal processes. VIII + 360 S. New York/London 1962. John Wiley & Sons. Preis 90/–.  ZAMM - Journal of Applied Mathematics and Mechanics / Zeitschrift für Angewandte Mathematik und Mechanik, 43(10-11):514–515, 1963\.
  * (9) 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-Bianchi, and R. Garnett, editors, Advances in Neural Information Processing Systems, volume 31. Curran Associates, Inc., 2018.
  * (10) Kyunghyun Cho, Bart Van Merriënboer, Caglar Gulcehre, Dzmitry Bahdanau, Fethi Bougares, Holger Schwenk, and Yoshua Bengio.  Learning phrase representations using rnn encoder-decoder for statistical machine translation.  arXiv preprint arXiv:1406.1078, 2014.
  * (11) Michael A. Cohen and Stephen Grossberg.  Absolute stability of global pattern formation and parallel memory storage by competitive neural networks.  IEEE Transactions on Systems, Man, and Cybernetics, SMC-13(5):815–826, 1983.
  * (12) Richard V Craster and Omar K Matar.  Dynamics and stability of thin liquid films.  Reviews of modern physics, 81(3):1131, 2009.
  * (13) J.R. Dormand and P.J. Prince.  A family of embedded runge-kutta formulae.  Journal of Computational and Applied Mathematics, 6(1):19–26, 1980\.
  * (14) Emilien Dupont, Arnaud Doucet, and Yee Whye Teh.  Augmented neural odes.  In Advances in Neural Information Processing Systems, volume 32. Curran Associates, Inc., 2019.
  * (15) Sourav Dutta, Peter Rivera-Casillas, Orie M Cecil, Matthew W Farthing, Emma Perracchione, and Mario Putti.  Data-driven reduced order modeling of environmental hydrodynamics using deep autoencoders and neural odes.  arXiv preprint arXiv:2107.02784, 2021.
  * (16) Sourav Dutta, Peter Rivera-Casillas, and Matthew W Farthing.  Neural ordinary differential equations for data-driven reduced order modeling of environmental hydrodynamics.  arXiv preprint arXiv:2104.13962, 2021.
  * (17) Massimo Germano, Ugo Piomelli, Parviz Moin, and William H Cabot.  A dynamic subgrid-scale eddy viscosity model.  Physics of Fluids A: Fluid Dynamics, 3(7):1760–1765, 1991.
  * (18) Serkan Gugercin and Athanasios C. Antoulas.  A Survey of Model Reduction by Balanced Truncation and Some New Results.  International Journal of Control, 77(8):748–766, May 2004.
  * (19) Amiram Harten, Peter D. Lax, and Bram van Leer.  On Upstream Differencing and Godunov-Type Schemes for Hyperbolic Conservation Laws.  SIAM Review, 25(1):35–61, 1983.
  * (20) Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun.  Identity mappings in deep residual networks.  In European Conference on Computer Vision, pages 630–645, 2016\.
  * (21) S. Hochreiter and J. Schmidhuber.  Long short-term memory.  Neural Computation, 9(8):1735–1780, 1997.
  * (22) Sepp Hochreiter and Jürgen Schmidhuber.  Long short-term memory.  Neural Computation, 9(8):1735–1780, 1997.
  * (23) Arthur Jacot, Franck Gabriel, and Clément Hongler.  Neural tangent kernel: convergence and generalization in neural networks.  In Proceedings of the 32nd International Conference on Neural Information Processing Systems, pages 8580–8589, 2018.
  * (24) J Nagoor Kani and Ahmed H Elsheikh.  Dr-rnn: A deep residual recurrent neural network for model reduction.  arXiv preprint arXiv:1709.00939, 2017.
  * (25) J Nagoor Kani and Ahmed H Elsheikh.  Reduced-order modeling of subsurface multi-phase flow models using deep residual recurrent neural networks.  Transport in Porous Media, 126(3):713–741, 2019.
  * (26) Diederik P Kingma and Max Welling.  Auto-encoding variational bayes.  arXiv preprint arXiv:1312.6114, 2013.
  * (27) Yoshiki Kuramoto.  Diffusion-Induced Chaos in Reaction Systems.  Progress of Theoretical Physics Supplement, 64:346–367, 02 1978\.
  * (28) Alexander Kurganov, Guergana Petrova, and Bojan Popov.  Adaptive Semidiscrete Central-Upwind Schemes for Nonconvex Hyperbolic Conservation Laws.  SIAM Journal on Scientific Computing, 29(6):2381–2401, January 2007\.
  * (29) Mathias Lechner and Ramin Hasani.  Learning long-term dependencies in irregularly-sampled time series.  arXiv preprint arXiv:2006.04418, 2020.
  * (30) Y.C. Liang, H.P. Lee, S.P. Lim, W.Z. Lin, K.H. Lee, and C.G. Wu.  Proper orthogonal decomposition and its applications—part i: Theory.  Journal of Sound and Vibration, 252(3):527–544, 2002.
  * (31) Hugo FS Lui and William R Wolf.  Construction of reduced-order models for fluid flows using deep feedforward neural networks.  Journal of Fluid Mechanics, 872:963–994, 2019.
  * (32) Chao Ma, Jianchun Wang, et al.  Model reduction with memory and the machine learning of dynamical systems.  arXiv preprint arXiv:1808.04258, 2018.
  * (33) Andrea Mannarino and Paolo Mantegazza.  Nonlinear aeroelastic reduced order modeling by recurrent neural networks.  Journal of Fluids and Structures, 48:103–121, 2014.
  * (34) Stefano Massaroli, Michael Poli, Jinkyoo Park, Atsushi Yamashita, and Hajime Asama.  Dissecting neural odes.  In H. Larochelle, M. Ranzato, R. Hadsell, M. F. Balcan, and H. Lin, editors, Advances in Neural Information Processing Systems, volume 33, pages 3952–3963. Curran Associates, Inc., 2020.
  * (35) Romit Maulik, Bethany Lusch, and Prasanna Balaprakash.  Reduced-order modeling of advection-dominated systems with recurrent neural networks and convolutional autoencoders.  Physics of Fluids, 33(3):037106, 2021.
  * (36) M. Mohebujjaman, L.G. Rebholz, and T. Iliescu.  Physically constrained data-driven correction for reduced-order modeling of fluid flows.  International Journal for Numerical Methods in Fluids, 89(3):103–122, 2019.
  * (37) Parviz Moin and Krishnan Mahesh.  Direct numerical simulation: a tool in turbulence research.  Annual review of fluid mechanics, 30(1):539–578, 1998.
  * (38) Changhong Mou, Honghu Liu, David R Wells, and Traian Iliescu.  Data-driven correction reduced order models for the quasi-geostrophic equations: A numerical investigation.  International Journal of Computational Fluid Dynamics, 34(2):147–159, 2020.
  * (39) Takaaki Murata, Kai Fukami, and Koji Fukagata.  Nonlinear mode decomposition with convolutional neural networks for fluid dynamics.  Journal of Fluid Mechanics, 882:A13, 2020.
  * (40) Tan Nguyen, Richard Baraniuk, Andrea Bertozzi, Stanley Osher, and Bao Wang.  MomentumRNN: Integrating momentum into recurrent neural networks.  In H. Larochelle, M. Ranzato, R. Hadsell, M. F. Balcan, and H. Lin, editors, Advances in Neural Information Processing Systems, volume 33, pages 1924–1936. Curran Associates, Inc., 2020.
  * (41) Alexander Norcliffe, Cristian Bodnar, Ben Day, Nikola Simidjievski, and Pietro Lió.  On second order behaviour in augmented neural odes.  In H. Larochelle, M. Ranzato, R. Hadsell, M. F. Balcan, and H. Lin, editors, Advances in Neural Information Processing Systems, volume 33, pages 5911–5921. Curran Associates, Inc., 2020.
  * (42) Razvan Pascanu, Tomas Mikolov, and Yoshua Bengio.  On the difficulty of training recurrent neural networks.  In International Conference on Machine Learning, pages 1310–1318, 2013.
  * (43) Karl Pearson.  Liii. on lines and planes of closest fit to systems of points in space.  The London, Edinburgh, and Dublin philosophical magazine and journal of science, 2(11):559–572, 1901.
  * (44) Boris T Polyak.  Some methods of speeding up the convergence of iteration methods.  USSR Computational Mathematics and Mathematical Physics, 4(5):1–17, 1964.
  * (45) Carlos J. G. Rojas, Andreas Dengel, and Mateus Dias Ribeiro.  Reduced-order Model for Fluid Flows via Neural Ordinary Differential Equations.  arXiv:2102.02248 [physics], February 2021.  arXiv: 2102.02248.
  * (46) Frank Rosenblatt.  Principles of neurodynamics. perceptrons and the theory of brain mechanisms.  Technical report, Cornell Aeronautical Lab Inc Buffalo NY, 1961.
  * (47) Yulia Rubanova, Ricky T. Q. Chen, and David K Duvenaud.  Latent ordinary differential equations for irregularly-sampled time series.  In Advances in Neural Information Processing Systems, volume 32. Curran Associates, Inc., 2019.
  * (48) Omer San and Romit Maulik.  Neural network closures for nonlinear model order reduction.  arXiv preprint arXiv:1705.08532, 2017.
  * (49) Omer San and Romit Maulik.  Machine learning closures for model order reduction of thermal fluids.  Applied Mathematical Modelling, 60:681–710, 2018.
  * (50) Omer San, Romit Maulik, and Mansoor Ahmed.  An artificial neural network framework for reduced order modeling of transient flows.  Communications in Nonlinear Science and Numerical Simulation, 77:271–287, 2019.
  * (51) Peter J. Schmid.  Dynamic mode decomposition of numerical and experimental data.  Journal of Fluid Mechanics, 656:5–28, 2010.
  * (52) G. I. Sivashinsky.  On flame propagation under conditions of stoichiometry.  SIAM Journal on Applied Mathematics, 39(1):67–82, 1980.
  * (53) G.I. Sivashinsky.  Nonlinear analysis of hydrodynamic instability in laminar flames—i. derivation of basic equations.  Acta Astronautica, 4(11):1177–1206, 1977.
  * (54) Tao Sun, Huaming Ling, Zuoqiang Shi, Dongsheng Li, and Bao Wang.  Training deep neural networks with adaptive momentum inspired by the quadratic optimization.  arXiv preprint arXiv:2110.09057, 2021.
  * (55) Bao Wang, Tan M Nguyen, Andrea L Bertozzi, Richard G Baraniuk, and Stanley J Osher.  Scheduled restart momentum for accelerated stochastic gradient descent.  arXiv preprint arXiv:2002.10583, 2020.
  * (56) Bao Wang, Hedi Xia, Tan Nguyen, and Stanley Osher.  How does momentum benefit deep neural networks architecture design? a few case studies.  arXiv preprint arXiv:2110.07034, 2021.
  * (57) Bao Wang and Qiang Ye.  Stochastic gradient descent with nonlinear conjugate gradient-style adaptive momentum.  arXiv preprint arXiv:2012.02188, 2020.
  * (58) Mingliang Wang, Han-Xiong Li, Xin Chen, and Yun Chen.  Deep learning-based model reduction for distributed parameter systems.  IEEE Transactions on Systems, Man, and Cybernetics: Systems, 46(12):1664–1674, 2016.
  * (59) Hedi Xia, Vai Suliafu, Hangjie Ji, Tan Nguyen, Andrea Bertozzi, Stanley Osher, and Bao Wang.  Heavy ball neural ordinary differential equation.  In Advances in Neural Information Processing Systems, volume 34. Curran Associates, Inc., 2021.
  * (60) Donghyun You and Parviz Moin.  A dynamic global-coefficient subgrid-scale eddy-viscosity model for large-eddy simulation in complex geometries.  Physics of Fluids, 19(6):065110, 2007.