# Non-intrusive Nonlinear Model Reduction via Machine Learning Approximations
to Low-dimensional Operators

Zhe Bai
Computational Research
Lawrence Berkeley National Lab
Berkeley, CA 94720
& Liqian Peng
Facebook AI Applied Research
Facebook
Menlo Park, CA 94025

###### Abstract

Although projection-based reduced-order models (ROMs) for parameterized
nonlinear dynamical systems have demonstrated exciting results across a range
of applications, their broad adoption has been limited by their intrusivity:
implementing such a reduced-order model typically requires significant
modifications to the underlying simulation code. To address this, we propose a
method that enables traditionally intrusive reduced-order models to be
accurately approximated in a non-intrusive manner. Specifically, the approach
approximates the low-dimensional operators associated with projection-based
reduced-order models (ROMs) using modern machine-learning regression
techniques. The only requirement of the simulation code is the ability to
export the velocity given the state and parameters as this functionality is
used to train the approximated low-dimensional operators. In addition to
enabling nonintrusivity, we demonstrate that the approach also leads to very
low computational complexity, achieving up to a three order of magnitude
reduction in run time. We demonstrate the effectiveness of the proposed
technique on two types of PDEs.

_Keywords_ model reduction, machine learning, low-dimensional operators,
dynamical systems

##  1 Introduction

Modern computational architectures have enabled the detailed numerical
simulation of incredibly complex physical and engineering systems at a vast
range of scales in both space and time (Lee et al., 2013). Even with improved
high-performance computing, the iterative numerical solutions required for
design and optimization may quickly become intractable; the computational
demands for real-time feedback control are even more challenging (Brunton and
Noack, 2015). Fortunately, many high-dimensional dynamical systems, such as a
discretized simulation of a fluid flow, are characterized by energetic
coherent structures that evolve on an attractor or manifold with a lower
dimensional intrinsic dimension (Aubry et al., 1988; Berkooz et al., 1993;
Holmes et al., 1996). This observed low-dimensionality has formed the basis of
reduced-order modeling, which now plays a central role in modern design,
optimization, and control of complex systems. Despite the growing importance
of model reduction, the challenges of nonlinearity, identifying an effective
low-dimensional basis, and multiscale dynamics have limited its widespread use
across the engineering and natural sciences (Brunton and Kutz, 2018). A
leading method for reduced-order modeling involves the Galerkin projection of
known governing equations onto a mathematical or empirical basis, although
this approach is intrusive and challenging to implement. The leading
alternative involves black-box system identification, which is purely data-
driven, but generally yields uninterpretable models. In this work, we
investigate and compare several emerging techniques from machine learning,
i.e. applied data-driven optimization, for non-intrusive reduced-order
modeling.

Intrusive model reduction methods, based on a working and decomposable
numerical simulation of the governing equations, provide the most general and
widely used set of techniques. Foremost in this arsenal is the Galerkin
projection of the governing equations onto a low-dimensional linear subspace,
usually spanned by orthogonal modes, such as Fourier, or data-driven modes
from proper orthogonal decomposition (POD) (Aubry et al., 1988; Berkooz et
al., 1993; Holmes et al., 1996; Noack et al., 2003). These approaches benefit
from a strong connection to the underlying physics, the ability to include
constraints to enforce preservation of the underlying dynamic structure, and
adaptivity (Amsallem et al., 2009; Carlberg et al., 2013, 2015, 2017). In
addition, there are several extensions around hyper-reduction to improve
efficiency for systems with parametric dependence or higher order
nonlinearities, based on the empirical interpolation method (EIM) (Barrault et
al., 2004) and the discrete empirical interpolation method (DEIM)
(Chaturantabut and Sorensen, 2010). However, these models are limited to the
dynamic terms of the governing equations and the linear projection basis may
not be optimal for capturing the dynamics. The main drawback of intrusive
methods is that they are challenging to implement, requiring considerable
human effort and knowledge to manipulate the governing equations and
simulation code.

In contrast to intrusive model reduction, data-driven techniques are becoming
increasingly prevalent, catalyzed by the increasing availability of
measurement data, the lack of governing equations for many modern systems, and
emerging methods in machine learning and data-driven optimization.
Collectively, these data-driven techniques form the field of system
identification (Ljung, 1999; Juang, 1994; Nelles, 2013). Many techniques in
system identification use regression to identify linear models, such as the
eigensystem realization algorithm (ERA) (Juang and Pappa, 1985) and dynamic
mode decomposition (DMD) (Schmid, 2010; Tu et al., 2014; Kutz et al., 2016;
Bai et al., 2020); recently, both techniques have been connected to nonlinear
systems via the Koopman operator (Mezić, 2005; Rowley et al., 2009; Brunton et
al., 2017). Another classic technique in nonlinear system identification is
the NARMAX algorithm (Kukreja and Brenner, 2007; Billings, 2013; Semeraro et
al., 2016). Recently, powerful techniques in machine learning are re-defining
what is possible in system identification, leveraging increasingly large
datasets and more powerful optimization and regression techniques. Deep
learning, based on multi-layer neural networks, is increasingly used to model
fluid dynamics and obtain closure models (Milano and Koumoutsakos, 2002; Ling
et al., 2016; Zhang and Duraisamy, 2015; Singh et al., 2017; Duraisamy et al.,
2018; Wang et al., 2017; Kutz, 2017). More generally, machine learning is
having a considerable impact in the modeling of dynamical systems and physics
(Bongard and Lipson, 2007; Schmidt and Lipson, 2009; Brunton et al., 2016;
Raissi and Karniadakis, 2018), for example relevant work in cluster-based
reduced order models (Kaiser et al., 2014), long-short time memory networks
(LSTMs) (Wan et al., 2018; Vlachas et al., 2018), and Gaussian process
regression (Raissi and Karniadakis, 2017; Wan and Sapsis, 2017). Of particular
interest are techniques based on the principle of _parsimony_ that seek the
simplest models, with the fewest terms necessary to describe the observed
dynamics (Bongard and Lipson, 2007; Schmidt and Lipson, 2009; Brunton et al.,
2016; Rudy et al., 2017; Loiseau and Brunton, 2018; Loiseau et al., 2018;
Schaeffer, 2017). For example, the sparse identification of nonlinear dynamics
(SINDy) algorithm (Brunton et al., 2016) uses sparse regression to identify
the few active terms in a differential equation. The resulting models balance
model complexity with descriptive power, avoiding overfitting and remaining
both interpretable and generalizable.

Data-driven techniques in system identification and machine learning are
currently being employed for advanced, non-intrusive, and reduced-order
modeling, removing the need for full-order model equations or a modifiable
numerical code. Loiseau and Brunton (Loiseau and Brunton, 2018) extended the
SINDy modeling framework to enforce known constraints, such as energy
conserving quadratic nonlinearities in incompressible flows. This so-called
_Galerkin regression_ also enables higher order nonlinearities than are
possible with standard Galerkin projection, providing effective closure for
the dynamics of truncated modes. Swischuk et al. (Swischuk et al., 2018)
develop parametric models between input parameters and POD coefficients for
physics and engineering problems, comparing many leading methods in machine
learning. Peherstorfer and Willcox (Peherstorfer and Willcox, 2016) develop
time-dependent reduced order models for partial differential equations, by
non-intrusively inferring the reduced-order model operator output as a
function of initial conditions, inputs, trajectories of states, and full-order
model outputs, without requiring access to the full-order model. They prove
that for dynamics that are linear in the state or are characterized by low-
order polynomial nonlinearities, the ROM converges to the projection of the
FOM onto a low-order subspace. Carlberg et al. (Carlberg et al., 2019) propose
dimensionality reduction on autoencoders to learn dynamics for recovering
missing CFD data. Other non-intrusive ROMs have been developed based on neural
networks (Wang et al., 2016; Kani and Elsheikh, 2017; Xie et al., 2019;
Hesthaven and Ubbiali, 2018; Wang et al., 2018), radial basis functions (Xiao
et al., 2015), and kernel methods (Vaerenbergh, 2010; Wirtz et al., 2015).
Hybrid methods are also promising, for example by modeling the error of a
Galerkin model to produce closures (Xie et al., 2018).

Despite the considerable recent progress in leveraging machine learning for
non-intrusive reduced-order modeling, current approaches still have a number
of shortcomings. First, many methods are limited to low-order polynomial
nonlinearities. All of the above approaches learn the low-dimensional
operators from full-system trajectories, which is expensive and limits the
predictive accuracy, as the ROM trajectory is likely to quickly deviate from
the full-system trajectory. Current methods also have limited consideration
for stability and the interplay between regression methods and time
integration has not been explored. Finally, current approaches do not provide
a framework for model selection, as each method will likely be best suited to
a different problem class.

In this work, we continue to develop and explore machine learning techniques
for non-intrusive model reduction, addressing many of the shortcomings
described above. In particular, we develop and compare numerous leading
techniques in machine learning to produce accurate and efficient reduced order
models. We focus primarily on the class of Galerkin projection models,
although these approaches may be extended to more general models. To study the
dynamical system in the latent space, we map from the low-dimensional state at
the current time instance to the next, in explicit or implicit time
integration schemes. We investigate a variety of regression models for the
mapping, including k-nearest neighbors (Altman, 1992), support vector
regression (Smola and Schölkopf, 2004), random forest (Breiman, 2001), boosted
trees (Hastie et al., 2009), the vectorial kernel orthogonal greedy algorithm
(VKOGA) (Wirtz and Haasdonk, 2013; Wirtz et al., 2015) and SINDy (Brunton et
al., 2016). In the following, we explicitly consider stability and explore the
connection between various regression methods and their suitability with
different time integration schemes; for example, non-differentiable machine
learning models are only amenable to explicit time integration. We also
explore the modeling procedures of varying complexities on two examples, the
1D Burgers’ equation and the 2D convection-diffusion equation using both
explicit and implicit time integrators.

##  2 Problem formulation

###  2.1 General nonlinear dynamical systems

This work considers parameterized nonlinear dynamical systems characterized by
the system of nonlinear ODEs

$\displaystyle\dot{\boldsymbol{x}}$
$\displaystyle=\boldsymbol{f}(\boldsymbol{x};t,\boldsymbol{\mu}),\qquad\boldsymbol{x}(0;\boldsymbol{\mu})$
$\displaystyle=\boldsymbol{x}^{0}(\boldsymbol{\mu}),$ (1)

where $t\in[0,T]$ denotes time with $T\in\mathbb{R}_{+}{}$ representing the
final time, $\boldsymbol{x}\in\mathbb{R}^{N}$ denotes the state,
$\boldsymbol{\mu}\in\mathcal{D}\subseteq\mathbb{R}^{p}$ denotes the system
parameters, $\boldsymbol{x}^{0}:\mathcal{D}\rightarrow\mathbb{R}^{N}$ is the
parameterized initial condition, and
$\boldsymbol{f}:\mathbb{R}^{N}\times\mathbb{R}_{+}{}\times\mathcal{D}\rightarrow\mathbb{R}^{N}$
denotes the velocity.

###  2.2 Galerkin projection

We assume that we we have low-dimensional trial-basis matrix
$\boldsymbol{V}\in\mathbb{R}_{\star}^{N\times n}$ (with $n\ll N$) computed,
e.g., via proper orthogonal decomposition (POD), such that the solution can be
approximated as
$\boldsymbol{x}(t,\boldsymbol{\mu})\approx\tilde{\boldsymbol{x}}(t,\boldsymbol{\mu})=\bar{\boldsymbol{x}}(\boldsymbol{\mu})+\boldsymbol{V}\hat{\boldsymbol{x}}(t,\boldsymbol{\mu})\approx\boldsymbol{x}$
with $\hat{\boldsymbol{x}}\in\mathbb{R}^{n}$ denoting the reduced state. Then,
the Galerkin ROM can be expressed as

$\displaystyle\dot{\hat{\boldsymbol{x}}}$
$\displaystyle=\boldsymbol{V}^{T}\boldsymbol{f}(\boldsymbol{V}\hat{\boldsymbol{x}},t;\boldsymbol{\mu})$
(2) $\displaystyle\hat{\boldsymbol{x}}(0)$
$\displaystyle=\boldsymbol{V}^{T}\boldsymbol{x}^{0}\left(\boldsymbol{\mu}\right).$
(3)

Critically, note that the ROM is defined by the low-dimensional mapping

$\displaystyle\boldsymbol{f}_{r}:$
$\displaystyle(\hat{\boldsymbol{y}},\tau;\boldsymbol{\nu})\mapsto\boldsymbol{V}^{T}\boldsymbol{f}(\boldsymbol{V}\hat{\boldsymbol{y}},\tau;\boldsymbol{\nu})$
(4) $\displaystyle:$
$\displaystyle\mathbb{R}^{n}\times\mathbb{R}_{+}{}\times\mathbb{R}^{p}\rightarrow\mathbb{R}^{n}.$
(5)

The reduced velocity $\boldsymbol{f}_{r}$ is thus simply a function that maps
the reduced state and inputs to a low-dimensional vector. Thus, we can rewrite
Equation (2) as

$\dot{\hat{\boldsymbol{x}}}=\hat{\boldsymbol{f}}_{r}(\hat{\boldsymbol{x}},t;\boldsymbol{\mu})$
(6)

However, this approach is intrusive to implement in computational-mechanics
codes, as it requires querying the full-order model code to compute
$\boldsymbol{f}(\boldsymbol{V}\hat{\boldsymbol{x}},t;\boldsymbol{\mu})$ for
every instance of $\hat{\boldsymbol{x}}$ and $\boldsymbol{\mu}$ encountered
during the ROM simulation; without hyper-reduction (e.g., DEIM or gappy POD),
it is also computationally expensive, as computing
$\boldsymbol{V}^{T}\boldsymbol{f}$ given $\boldsymbol{f}$ incurs
$\mathcal{O}(Nn)$ flops. Even if hyper-reduction is present, however, the
approach remains intrusive, as sampled entries of the velocity must be
computed for every instance of $\hat{\boldsymbol{x}}$ and $\boldsymbol{\mu}$
encountered during the ROM simulation.

##  3 Regression-based reduced operator approximation

This section outlines learning the low dimensional operator through regression
methods. We commence describing the general formulation of the framework in
Section 3.1, and then Section 3.2 discusses the proposed regression models and
their computational complexity. In Section 3.3, we derive the boundedness and
stability of of the approximated discrete-time dynamical system.

###  3.1 Mathematical formulation

The objective is to develop a data-driven, non-intrusive approximation to the
reduced velocity with reduced-order models of different dynamical systems.
$\boldsymbol{f}_{r}$. In particular, we aim to devise an approximate mapping

$\displaystyle\hat{\boldsymbol{f}}_{r}:$
$\displaystyle\mathbb{R}^{n}\times\mathbb{R}_{+}{}\times\mathbb{R}^{p}\rightarrow\mathbb{R}^{n}.$
(7)

such that

$\displaystyle\hat{\boldsymbol{f}}_{r}(\hat{\boldsymbol{x}},t;\boldsymbol{\mu})\approx\boldsymbol{f}_{r}(\hat{\boldsymbol{x}},t;\boldsymbol{\mu}),\quad\forall(\hat{\boldsymbol{x}},t;\boldsymbol{\mu})\in\mathbb{R}^{n}\times\mathbb{R}_{+}{}\times\mathbb{R}^{p}.$
(8)

As the domain and codomain of $\hat{\boldsymbol{f}}_{r}$ exhibit no
specialized structure, we can consider general high-dimensional regression
techniques (e.g., kernel regression, random-forest regression).

Similar to Equation (6), the reduced order model is given by

$\dot{\hat{\boldsymbol{x}}}=\hat{\boldsymbol{f}}_{r}(\hat{\boldsymbol{x}},t;\boldsymbol{\mu}),$
(9)

with the initial condition

$\hat{\boldsymbol{x}}(0)=\boldsymbol{V}^{T}\boldsymbol{x}^{0}\left(\boldsymbol{\mu}\right).$
(10)

We assume that the sequences of reduced velocity $\boldsymbol{f}_{r}$ can be
studied from the Markovian discrete-time dynamical system

$\hat{\boldsymbol{x}}^{j+1}=\boldsymbol{g}(\hat{\boldsymbol{x}}^{j},\boldsymbol{\mu}),\quad
j=0,\cdots,N_{t},$ (11)

for discrete-time velocity
$\boldsymbol{g}:\mathbb{R}^{n}\times\mathbb{R}^{p}\rightarrow\mathbb{R}^{n}$.
If such a reduced velocity exists and we can compute the operator, the whole
time sequence of the reduced states can then be estimated from Equation (11)
by specifying the initial (reduced) state. By considering the mapping between
the current reduced state and the next, we cast the problem into a set of
regression models.

###  3.2 Surrogate ROM

In this work, we collect the reduced states from the Galerkin projection as in
Section 2.2 and approximate the operator $\boldsymbol{g}$ by regressing each
component of the reduced velocity. To be specific, we construct individual
regression $\hat{g}_{i}\approx g_{i}$ for $i=1,\cdots,n$, the $i$th equation
of $\boldsymbol{g}$ is:

$\hat{x}_{i}^{j+1}=\hat{g}_{i}(\hat{\boldsymbol{x}}^{j},\boldsymbol{\mu}),\quad
j=0,\cdots,N_{t},\quad i=1,\cdots,n,$ (12)

where $\hat{\boldsymbol{x}}=[x_{1}\cdots x_{n}]$ and
$\boldsymbol{g}=[g_{1}\cdots g_{n}]$. Equation (12) shows that we can map the
reduced state at time step $\hat{\boldsymbol{x}}^{n}$ to each component of the
state at the next time step $\hat{x}_{i}^{n+1}$ for $j=0,\cdots,N_{t}-1$.

In the offline stage, we obtain training data from the existing Galerkin
projection model. By querying both features ($\hat{\boldsymbol{x}}^{j}$) and
responses ($\hat{x}_{i}^{j+1}$) at every time instance, we generate a data set
$\mathcal{T}_{train,i}={(\hat{\boldsymbol{x}}^{j},x_{i}^{j+1})}$ for model
training and cross validation.

We consider a variety of conventional statistical models, including support
vector regression (SVR) with kernel functions (Gaussian, polynomials) (Smola
and Schölkopf, 2004), tree-based methods (boosted decision trees (Hastie et
al., 2009), random forest (Breiman, 2001)) and $k$-nearest neighbors (Altman,
1992), for regressing the feature-response pairs. Three types of kernel
functions are explored in SVR. Specifically, we refer to the SVR model using
the 2nd, 3rd order polynomial kernel as SVR2, SVR3 respectively. When the
Gaussian radial basis function (RBF) kernel is used, we refer to the model as
SVRrbf. In addition, we investigate the vectorial kernel orthogonal greedy
algorithm (VKOGA) and Sparse identification of nonlinear dynamics (SINDy) as
advanced regression models. More details of VKOGA and SINDy can be found in
Appendix A.1. We compare the computational complexity of all the proposed
regression models in Appendix A.1. For all the candidate models, we employ
cross validation to optimize the hyperparameters for model selection, aiming
for a balanced bias and variance. We also explore the training and validation
error as varying the number of samples in Appendix A.2, and we select a fixed
sample size for performance comparison between models.

For time integration along the trajectory of the dynamical system, we
investigate the appropriate time step for each regression model. We implement
the $4$th-order Runge-Kutta as the explicit method, as well as Newton-Raphson
and fixed-point iteration both in backward Euler as the implicit methods. We
report the numerical results, with respect to the Galerkin reduced-order model
(ROM) and the full-order model (FOM) in a sequence of time in Section 4.

###  3.3 Error analysis

Assuming the considered regression model generates bounded in reduced space,
we examine the boundedness of the surrogate FOM on the time evolution of the
states along the trajectory. Let $J_{t}=[0,T]$ denote the time domain,
$J_{\mu}\subset\mathbb{R}^{p}$ be a compact space of $\boldsymbol{\mu}$,
$\boldsymbol{x}:J_{t}\to\mathbb{R}^{n}$ denote the state variable, and
$\boldsymbol{f}(\boldsymbol{x},t;\boldsymbol{\mu})$ be the vector field. Let
$\boldsymbol{f}_{r}(\hat{\boldsymbol{y}},t;\boldsymbol{\mu})=\boldsymbol{f}(\boldsymbol{V}\hat{\boldsymbol{y}},t;\boldsymbol{\mu})$
represent the Galerkin reduced vector field. Let $\hat{\boldsymbol{x}}(t)$
denote an approximate reduced vector field constructed by a machine learning
method and $\hat{\boldsymbol{f}}_{r}(\hat{\boldsymbol{x}},t;\boldsymbol{\mu})$
be the corresponding vector field.

The error of the reduced model can be defined as
$\boldsymbol{e}(t):=\hat{\boldsymbol{x}}(t)-\boldsymbol{x}(t)$. Let
$P=\boldsymbol{V}\boldsymbol{V}^{*}$. Let
$\boldsymbol{e}_{o}(t):=(I_{n}-P)\boldsymbol{e}(t)$, which denotes the error
component orthogonal to range($\boldsymbol{V}$), and
$\boldsymbol{e}_{i}(t):=P\boldsymbol{e}(t)$, which denotes the component of
error parallel to range($\boldsymbol{V}$). Let
$\boldsymbol{\tilde{x}}:=P\boldsymbol{x}$ denote the direct projection of
$\boldsymbol{x}$. Thus, we have

$\boldsymbol{e}_{o}(t)=\boldsymbol{\tilde{x}}(t)-\boldsymbol{x}(t).$ (13)

However, since the system is evolutionary with time, further approximations of
the projection-based reduced model result in an additional error
$\boldsymbol{e}_{i}(t)$, and we have

$\boldsymbol{e}_{i}(t)=\boldsymbol{V}\hat{\boldsymbol{x}}(t)-\boldsymbol{\tilde{x}}(t).$
(14)

Although $\boldsymbol{e}_{i}(t)$ and $\boldsymbol{e}_{o}(t)$ are orthogonal to
each other, they are not independent.

lemma Consider the dynamical system as Equation (1) over the interval $t\in J$
and $\boldsymbol{\mu}\in J_{\mu}$. Suppose
$\boldsymbol{f}(\boldsymbol{x},t;\boldsymbol{\mu})$ is a uniformly Lipschitz
function of $\boldsymbol{x}$ with constant $K$ and a continuous function of
$t$ for all $(\boldsymbol{x},t,\boldsymbol{\mu})\in
B_{b}(\boldsymbol{x}_{0})\times J_{t}\times J_{\mu}$. Suppose
$\hat{\boldsymbol{f}}_{r}(\hat{\boldsymbol{x}},t;\boldsymbol{\mu})$ is a
uniformly Lipschitz function of $\hat{\boldsymbol{x}}$ and a continuous
function of $t$ for all $(\boldsymbol{x},t;\boldsymbol{\mu})\in
B_{b}(x_{0})\times J_{t}\times J_{\mu}$. Suppose the data-driven approximation
satisfies
$\|\hat{\boldsymbol{f}}_{r}(\hat{\boldsymbol{x}},t;\boldsymbol{\mu})-\boldsymbol{V}^{*}\boldsymbol{f}(\boldsymbol{V}\hat{\boldsymbol{x}},t;\boldsymbol{\mu})\|\leq
C$ for all $(\hat{\boldsymbol{x}},t,\boldsymbol{\mu})\in
B_{b}(\hat{\boldsymbol{x}}(\times_{0}))\times J_{t}\times J_{\mu}$. Then the
error $\boldsymbol{e}(t)=\hat{\boldsymbol{x}}(t)-\boldsymbol{x}(t)$ in the
infinity norm for the interval $J_{t}$ is bounded by

$\left\|\boldsymbol{e}\right\|_{\infty}\leq
e^{KT}\left\|\boldsymbol{e}_{o}\right\|_{\infty}+e^{KT}\left\|{\boldsymbol{e}_{i}(0)}\right\|+\frac{C}{K}(e^{KT}-1)$
(15)

proof: For notation simplification, we fix $\boldsymbol{\mu}$ and do not
explicitly denote it in vector fields.

Substituting Equation (1) and (9) into the differentiation of
$\boldsymbol{e}_{o}(t)+\boldsymbol{e}_{i}(t)=\boldsymbol{V}\hat{\boldsymbol{x}}(t)-\boldsymbol{x}(t)$
yields

$\dot{\boldsymbol{e}}_{o}+\dot{\boldsymbol{e}}_{i}=\boldsymbol{V}\hat{\boldsymbol{f}}_{r}(\hat{\boldsymbol{x}},t)-\boldsymbol{f}(\boldsymbol{x},t).$
(16)

Left multiplying (16) by $P$ and recognizing that
$P\boldsymbol{V}=\boldsymbol{V}$ gives

$\displaystyle\dot{\boldsymbol{e}}_{i}(t)$
$\displaystyle=\boldsymbol{V}\hat{\boldsymbol{f}}_{r}(\hat{\boldsymbol{x}},t)-P\boldsymbol{f}(\boldsymbol{x},t)$
$\displaystyle=\boldsymbol{V}\hat{\boldsymbol{f}}_{r}(\hat{\boldsymbol{x}},t)-P\boldsymbol{f}(\boldsymbol{V}\hat{\boldsymbol{x}},t)+P\boldsymbol{f}(\boldsymbol{V}\hat{\boldsymbol{x}},t)-P\boldsymbol{f}(\boldsymbol{x}+\boldsymbol{e}_{o},t)+P\boldsymbol{f}(\boldsymbol{x}+\boldsymbol{e}_{o},t)-P\boldsymbol{f}(\boldsymbol{x},t)$
$\displaystyle=\boldsymbol{V}(\hat{\boldsymbol{f}}_{r}(\hat{\boldsymbol{x}},t)-\boldsymbol{V}^{*}\boldsymbol{f}(\boldsymbol{V}\hat{\boldsymbol{x}},t))+P(\boldsymbol{f}(\boldsymbol{x}+\boldsymbol{e}_{o}+\boldsymbol{e}_{i},t)-\boldsymbol{f}(\boldsymbol{x}+\boldsymbol{e}_{o},t))$
$\displaystyle\quad+P(\boldsymbol{f}(\boldsymbol{x}+\boldsymbol{e}_{o},t)-\boldsymbol{f}(\boldsymbol{x},t))$

Using this equation by expanding $\|\boldsymbol{e}_{i}(\boldsymbol{x}+h)\|$
and applying triangular inequality yields

$\displaystyle\|\boldsymbol{e}_{i}(t+h)\|$
$\displaystyle=\|\boldsymbol{e}_{i}(t)+h\dot{\boldsymbol{e}}_{i}(t)\|+\mathcal{O}({h^{2}})$
$\displaystyle\leq\|\boldsymbol{e}_{i}(t)\|+\|h\boldsymbol{V}(\hat{\boldsymbol{f}}_{r}(\hat{\boldsymbol{x}},t)-\boldsymbol{V}^{*}\boldsymbol{f}(\boldsymbol{V}\hat{\boldsymbol{x}},t))\|+\|hP(\boldsymbol{f}(\boldsymbol{x}+\boldsymbol{e}_{o}+\boldsymbol{e}_{i},t)$
$\displaystyle\quad-\boldsymbol{f}(\boldsymbol{x}+\boldsymbol{e}_{o},t))\|+\|hP(\boldsymbol{f}(\boldsymbol{x}+\boldsymbol{e}_{o},t)-\boldsymbol{f}(\boldsymbol{x},t))\|+\mathcal{O}(h^{2})$
$\displaystyle\leq\|\boldsymbol{e}_{i}(t)\|+h\|\boldsymbol{V}\|\cdot\|\hat{\boldsymbol{f}}_{r}(t,\hat{\boldsymbol{x}})-\boldsymbol{V}^{*}\boldsymbol{f}(\boldsymbol{V}\hat{\boldsymbol{x}},t)\|+h\|P\|\cdot\|\boldsymbol{f}(\boldsymbol{x}+\boldsymbol{e}_{o}+\boldsymbol{e}_{i},t)$
$\displaystyle\quad-\boldsymbol{f}(\boldsymbol{x}+\boldsymbol{e}_{o},t)\|+h\|P\|\cdot\|\boldsymbol{f}(\boldsymbol{x}+\boldsymbol{e}_{o},t)-\boldsymbol{f}(\boldsymbol{x},t)\|+\mathcal{O}(h^{2}).$

Using $\|\boldsymbol{V}\|=\|P\|=1$ and
$\|\hat{\boldsymbol{f}}_{r}(t,\hat{\boldsymbol{x}})-\boldsymbol{V}^{*}\boldsymbol{f}(t,\boldsymbol{V}\hat{\boldsymbol{x}})\|\leq
C$, the last inequality gives

$\displaystyle\|\boldsymbol{e}_{i}(t+h)\|$
$\displaystyle\leq\|\boldsymbol{e}_{i}(t)\|+hC+h\|\boldsymbol{f}(\boldsymbol{x}+\boldsymbol{e}_{o}+\boldsymbol{e}_{i},t)-\boldsymbol{f}(\boldsymbol{x},+\boldsymbol{e}_{o},t)\|$
$\displaystyle\quad+h\|\boldsymbol{f}(\boldsymbol{x}+\boldsymbol{e}_{o},t)-\boldsymbol{f}(\boldsymbol{x},t)\|+\mathcal{O}(h^{2}).$

Rearranging this inequality and applying the Lipschitz conditions gives

$\frac{\|\boldsymbol{e}_{i}(t+h)\|-\|\boldsymbol{e}_{i}(t)\|}{h}\leq
K\|\boldsymbol{e}_{i}(t)\|+K\|\boldsymbol{e}_{o}(t)\|+C+\mathcal{O}(h).$

Since $\mathcal{O}(h)$ can be uniformly bounded independent of
$\boldsymbol{e}_{i}(t)$, using the mean value theorem and letting $h\to 0$
give

$\frac{d}{dt}\|\boldsymbol{e}_{i}(t)\|\leq
K\|\boldsymbol{e}_{i}(t)\|+K\|\boldsymbol{e}_{o}(t)\|+C.$

Rewriting the above inequality into integral form,
$\left\|\boldsymbol{e}_{i}(t)\right\|\leq\alpha(t)+K\int_{0}^{t}\|\boldsymbol{e}_{i}(\tau)\|d\tau$,
where
$\alpha(t):=\left\|\boldsymbol{e}_{i}(0)\right\|+K\int_{0}^{t}{\|\boldsymbol{e}_{o}(\tau)\|d\tau}+Ct$,
and using Gronwall’s lemma, we obtain

$\left\|{{\boldsymbol{e}_{i}}(t)}\right\|\leq\alpha(t)+\int_{0}^{t}{\alpha(s)K\exp\left({\int_{s}^{t}{Kd\tau}}\right)ds}.$

By definition,
${\left\|\boldsymbol{e}_{o}\right\|_{\infty}}\geq\left\|\boldsymbol{e}_{o}(t)\right\|$
for any $t\in J_{t}$. It follows that
$\alpha(t)\leq\left\|{{\boldsymbol{e}_{i}}(0)}\right\|+Kt\left\|\boldsymbol{e}_{o}\right\|_{\infty}+Ct$.
Simplifying the integral of the right-hand side of the above inequality gives

$\|\boldsymbol{e}_{i}(t)\|\leq(e^{Kt}-1)\left(\|\boldsymbol{e}_{o}\|_{\infty}+\frac{C}{K}\right)+e^{Kt}\|\boldsymbol{e}_{i}(0)\|,$

for any $t\in J_{t}$. If follows that

$\|\boldsymbol{e}_{i}\|_{\infty}\leq(e^{KT}-1)\left(\|\boldsymbol{e}_{o}\|_{\infty}+\frac{C}{K}\right)+e^{KT}\|\boldsymbol{e}_{i}(0)\|,$

Combining the above inequality with
$\left\|\boldsymbol{e}\right\|_{\infty}\leq\left\|{\boldsymbol{e}_{i}}\right\|_{\infty}+\left\|{\boldsymbol{e}_{o}}\right\|_{\infty}$,
one can obtain Equation (15).

Remark: The above lemma provides a bound for $\|\boldsymbol{e}_{i}(t)\|$ in
terms of $\|\boldsymbol{e}_{o}\|_{\infty}$ and $\|\boldsymbol{e}_{i}(0)\|$. We
have $\left\|\boldsymbol{e}_{i}(0)\right\|=0$ when the initial condition of
the reduced model is given by
$\hat{\boldsymbol{x}}(0)=\boldsymbol{V}^{*}\boldsymbol{x}_{0}$ for Equation
(9). In this situation, Equation (15) becomes
${\left\|\boldsymbol{e}\right\|_{\infty}}\leq
e^{KT}{\left\|{{\boldsymbol{e}_{o}}}\right\|_{\infty}}+\frac{C}{K}(e^{KT}-1)$.

##  4 Numerical experiments

To assess the proposed non-intrusive ROM, we consider two parameterized PDEs:
(i) 1D inviscid Burgers’ equation, and (ii) 2D convection–diffusion equation.
We implement explicit integrator, including 4th-order Runge-Kutta solvers, and
Newton-Raphson, fixed-point iteration in backward Euler as the implicit
methods.

###  4.1 1D Burgers’ equation

The experiment first employs a 1D parameterized inviscid Burgers’ equation.
The input parameters $\boldsymbol{\mu}=(a,b)$ in the space
$[1.5,2]\times[0.02,0.025]$. In the current setup, the parameters for online
test are fixed to be $\boldsymbol{\mu}=(1.8,0.0232)$. In the FOM, the problem
is often solved using a conservative finite-volume formulation and Backward
Euler in time. The 1D domain is discretized using a grid with 501 nodes,
corresponding to $x=i\times(100/500),i=0,\cdots,500$. The solution $u(x,t)$ is
computed in the time interval $t\in[0,25]$ using different time step sizes
considering the convergence in each time integrator.

$\displaystyle\frac{\partial u(x,t)}{\partial
t}+\frac{1}{2}\frac{\partial(u^{2}(x,t))}{\partial x}$
$\displaystyle=0.02e^{bx},$ (17) $\displaystyle u(0,t)$
$\displaystyle=a,\forall t>0,$ (17a) $\displaystyle u(x,0)$
$\displaystyle=1,\forall x\in[0,100]$ (17b)

The solution $u(\boldsymbol{x},t)$ is computed in the time interval
$t\in[0,25]$ using a uniform computational time-step size $\Delta t$.

####  4.1.1 Data collection

\begin{overpic}[width=411.93767pt]{figures_1st/Fig0_RK} \put(3.0,32.0){(a)}
\end{overpic}

\begin{overpic}[width=411.93767pt]{figures_1st/Fig0_BE} \put(3.0,32.0){(b)}
\put(55.0,31.0){\vector(1,-1){4.5}} \put(57.0,24.0){$\Delta t$} \end{overpic}

Figure 1: Timestep verification for 1D inviscid Burgers’ equation: (a) Runge-
Kutta: we select a time step of $1.25e{\text{-}}1$, as it leads to an
approximated ($<1\%$ error) convergence rate of $4$ and an averaged
approximated error of $6e{\text{-}}5$ for the selected state. (b) backward-
Euler integrator: we select a time step of $3.12e{\text{-}}2$, as it leads to
an approximated convergence rate of 1 and an approximated error of
$4e{\text{-}}4$ for the selected state.

We investigate time step verification on choosing an appropriate $\Delta t$ of
the time integrator for the problem. We collect the solutions under an
increasing number of time steps $N_{t}=[25,50,100,200,400,800,1600,3200,6400]$
using both Runge-Kutta as well as backward Euler integrator. Throughout the
paper, we select the time step at $99\%$ of the asymptotic rate of
convergence. The verification results in Figure 1 show that $N_{t}=200$ is a
reasonable number of time steps to use for the 4th-order Runge-Kutta and
$nT=800$ for backward Euler method. During the offline stage, we run four full
simulations corresponding to the corner parameters of the space
$[1.5,2]\times[0.02,0.025]$. Then, we sample the training data from Latin-
hypercube, . In the sampling, $N_{\text{training}}$ and
$N_{\text{validation}}$ instances of the state, time and parameters are
generated following the criterion that the minimum distances between the data
points are maximized. For this study, the default size of the training set is
$N_{\text{training}}=1000$ and the default size of the validation set is
$N_{\text{validation}}=500$. The reduced vector field $\boldsymbol{f}_{r}$ is
computed for each input pairs $(\hat{\boldsymbol{x}},t;\boldsymbol{\mu})$.
Note that both the training and validation stage only involves pure machine
learning. Then in the test stage, we evaluate the proposed ROM. The parameters
are fixed to be $\mu=(1.8,0.0232)$ for testing purpose.

####  4.1.2 Model validation

We use SVR with kernel functions (2nd, 3rd order polynomials and radial basis
function), kNN, Random Forest, Boosting, VKOGA, and SINDy as regression models
to approximate reduced velocity. In particular, for each regression method, we
change the model hyperparameters and plot the relative training error and
validation error. The relative error is defined by

$err=\frac{\|\hat{\boldsymbol{f}}_{r}(\hat{\boldsymbol{x}},t;\boldsymbol{\mu})-\boldsymbol{f}_{r}(\hat{\boldsymbol{x}},t;\boldsymbol{\mu})\|}{\|\boldsymbol{f}_{r}(\hat{\boldsymbol{x}},t;\boldsymbol{\mu})\|}.$
(18)

We then plot the learning curve of each regression method and compare the
performance of each model on training and validation data over a varying
number of training instances in Appendix A.2. By properly choosing the
hyperparameters and the number of training instances, our regression models
can effectively balance bias and variance.

####  4.1.3 Simulation of the surrogate ROM

We can now solve the problem using the surrogate model along the trajectory of
the dynamical system. After applying time integration to the regression-based
ROM, we compute the relative error of the proposed models as a function of
time. We investigate both Newton-Raphson, fixed-point iteration in backward
Euler and 4th-order Runge-Kutta in explicit methods. Let $\boldsymbol{x}(t)$,
$\hat{\boldsymbol{x}_{b}}(t)$ and $\hat{\boldsymbol{x}}(t)$ denote the
solution of the FOM, the Galerkin ROM, and non-intrusive ROMs respectively. We
define the relative error with respect to FOM
$\boldsymbol{e}_{\textit{FOM}}(t)$ and Galerkin ROM
$\boldsymbol{e}_{\textit{ROM}}(t)$ as

$\boldsymbol{e}_{\textit{FOM}}(t)=\frac{\|(\boldsymbol{V}\hat{\boldsymbol{x}}(t)-\boldsymbol{x}(t))\|}{\|\boldsymbol{x}(t)\|},$
(19)
$\boldsymbol{e}_{\textit{ROM}}(t)=\frac{\|(\hat{\boldsymbol{x}}(t)-\hat{\boldsymbol{x}_{b}}(t))\|}{\|\hat{\boldsymbol{x}_{b}}(t)\|},$
(20)

The corresponding averaged relative error over the entire time domain
$\boldsymbol{e}_{\textit{FOM}}$ and $\boldsymbol{e}_{\textit{ROM}}$ can be
computed as

$\boldsymbol{e}_{\textit{FOM}}=\frac{1}{T}\int_{t=0}^{T}\boldsymbol{e}_{\textit{FOM}}(t)dt,$
(21)
$\boldsymbol{e}_{\textit{ROM}}=\frac{1}{T}\int_{t=0}^{T}\boldsymbol{e}_{\textit{ROM}}(t)dt.$
(22)

Let $t_{FOM}$, $t_{ROM}$, and $\tau$ denote the running time of FOM, Galerkin
ROM, and non-intrusive ROM respectively. Define the relative running time with
respect to the FOM and the Galerkin ROM:

$\tau_{\textit{FOM}}=\frac{\tau}{t_{\textit{FOM}}},$ (23)

and

$\tau_{\textit{ROM}}=\frac{\tau}{t_{\textit{ROM}}}.$ (24)

The following are the simulation results from backward Euler method with
$N_{t}=800$. Figure 2 plots the state-space error with respect to the FOM and
ROM using the backward Euler integrator. As validation results predict, SVR3
and SINDy behave better than the other models, achieving a relative ROM error
below $1e\text{-}4$ over the entire time domain, and the relative error in
terms of FOM is well bounded.

\begin{overpic}[width=195.12767pt]{figures_1st/FomErrVsTime_BE}
\put(-1.0,79.0){(a)} \put(49.0,-0.1){\hbox{\pagecolor{white}{Time}}}
\put(-2.0,32.0){\hbox{\pagecolor{white}\rotatebox{90.0}{{Relative error}}}}
\end{overpic}\begin{overpic}[width=195.12767pt]{figures_1st/RomErrVsTime_BE}
\put(-1.0,79.0){(b)} \put(49.0,-0.1){\hbox{\pagecolor{white}{Time}}}
\put(-2.0,32.0){\hbox{\pagecolor{white}\rotatebox{90.0}{{Relative error}}}}
\end{overpic}

Figure 2: Backward Euler for 1D inviscid Burgers’ equation: time evolution of
relative error: (a) $\boldsymbol{e}_{\textit{FOM}}(t)$ in FOM; and (b)
$\boldsymbol{e}_{\textit{ROM}}(t)$ in ROM.

Figure 3 plots the Pareto frontier error as a function of the relative running
time using the backward Euler integrator. For differential models, the
relative time is calculated using the less expensive approach, i.e. Newton’s
method. By comparison, SINDy requires much less relative time than SVR3, at a
comparable level of relative error in both FOM and ROM.

\begin{overpic}[width=195.12767pt]{figures_1st/Pareto_FOM_BE}
\put(-1.0,79.0){(a)}
\put(-2.0,32.0){\hbox{\pagecolor{white}\rotatebox{90.0}{{Relative error}}}}
\put(32.0,-0.1){\hbox{\pagecolor{white}{Relative time w.r.t. FOM}}}
\end{overpic}\begin{overpic}[width=203.80193pt]{figures_1st/Pareto_ROM_BE}
\put(-1.0,76.0){(b)}
\put(-2.0,32.0){\hbox{\pagecolor{white}\rotatebox{90.0}{{Relative error}}}}
\put(32.0,-0.1){\hbox{\pagecolor{white}{Relative time w.r.t. ROM}}}
\end{overpic}

Figure 3: Pareto frontier of relative error with respect to the relative
running time using backward Euler for 1D inviscid Burgers’ equation: (a)
$\boldsymbol{e}_{\textit{FOM}}$ vs. $\tau_{\textit{FOM}}$ in FOM; (b)
$\boldsymbol{e}_{\textit{ROM}}$ vs. $\tau_{\textit{ROM}}$ in ROM.

Table 2 summarizes (i) online running time of all methods, (ii) mean time-
integrated error with respect to the FOM, and (iii) mean time-integrated error
with respect to ROMs using the backward Euler integrator. For differentiable
models, we show the computational time of both Newton’s method and fixed point
iteration in backward Euler. For those models that are non-differentiable,
i.e. random forest, boosting and kNN, only the running time of the fixed-point
iteration method is reported. For a fair comparison, all the ROM and FOM
solutions are computed at the verified backward Euler time step
$\Delta{t}=3.12e{\text{-}}2$. Note that the non-intrusive ROM, e.g. SINDy with
Newton’s method can accelerate the solver by $10.4\times$ relative to the FOM
and $2.5\times$ compared to the Galerkin ROM at a relative error of $0.0346$
and $4.36e{\text{-}}5$ respectively.

Table 1: Comparison of different machine learning methods using backward Euler methods for 1D inviscid Burgers’ equation: Newton’s method (N); fixed-point iteration (FP). The running time of FOM and Galerkin is $3.886$s and $0.942$s respectively. Method | Online running time (s) | R Err (w.r.t. FOM) | R Err (w.r.t. Galerkin)
---|---|---|---
SVR 2 FP | 93.365 | 0.0346 | 1.62e-4
SVR 2 N | 53.772 | 0.0346 | 1.62e-4
SVR 3 FP | 102.873 | 0.0346 | 1.02e-5
SVR 3 N | 78.935 | 0.0346 | 1.02e-5
SVR rbf FP | 104.349 | 0.0353 | 8.12e-3
SVR rbf N | 104.548 | 0.0353 | 8.12e-3
Random Forest | 374.201 | 0.0756 | 0.0638
Boosting | 587.473 | 0.0980 | 0.0873
kNN | 1.108 | 0.0558 | 0.0427
VKOGA FP | 2.505 | 0.0365 | 0.0122
VKOGA N | 1.217 | 0.0365 | 0.0122
SINDy FP | 0.534 | 0.0346 | 4.36e-5
SINDy N | 0.373 | 0.0346 | 4.36e-5
Galerkin | 0.942 | 0.0346 | 0

We examine simulation results from 4th-order Runge-Kutta method with
$N_{t}=200$. Figure 4 shows the state-space error with respect to the FOM and
ROM using the Runge-Kutta integrator. SVR2, SVR3 and SINDy have a comparable
performance, and result in a bigger ROM error relative to the backward Euler
solver. We notice that the random forest model begins to diverge after $t>10$
in the explicit scheme. This can be explained by the corresponding performance
in model evaluation in Appendix A.1.

\begin{overpic}[width=195.12767pt]{figures_1st/FomErrVsTime_RK4}
\put(-1.0,78.0){(a)} \put(49.0,-0.1){\hbox{\pagecolor{white}{Time}}}
\put(-2.0,32.0){\hbox{\pagecolor{white}\rotatebox{90.0}{{Relative error}}}}
\end{overpic}\begin{overpic}[width=195.12767pt]{figures_1st/RomErrVsTime_RK4}
\put(-1.0,78.0){(b)} \put(49.0,-0.1){\hbox{\pagecolor{white}{Time}}}
\put(-2.0,32.0){\hbox{\pagecolor{white}\rotatebox{90.0}{{Relative error}}}}
\end{overpic}

Figure 4: Runge-Kutta for 1D inviscid Burgers’ equation: time evolution of
relative error: (a) $\boldsymbol{e}_{\textit{FOM}}(t)$ in FOM; and (b)
$\boldsymbol{e}_{\textit{ROM}}(t)$ in ROM.

Figure 5 plots the Pareto frontier error with respect to the relative running
time using the backward Euler integrator. VKOGA has the smallest relative time
in both the FOM and ROM comparison. SINDy requires slightly more running time
while the accuracy outperforms VKOGA.

\begin{overpic}[width=195.12767pt]{figures_1st/Pareto_FOM_RK4}
\put(-1.0,79.0){(a)}
\put(-2.0,32.0){\hbox{\pagecolor{white}\rotatebox{90.0}{{Relative error}}}}
\put(32.0,-0.1){\hbox{\pagecolor{white}{Relative time w.r.t. FOM}}}
\end{overpic}\begin{overpic}[width=195.12767pt]{figures_1st/Pareto_ROM_RK4}
\put(-1.0,79.0){(b)}
\put(-2.0,32.0){\hbox{\pagecolor{white}\rotatebox{90.0}{{Relative error}}}}
\put(32.0,-0.1){\hbox{\pagecolor{white}{Relative time w.r.t. ROM}}}
\end{overpic}

Figure 5:  Pareto frontier of relative error with respect to the relative
running time using 4th-order Runge-Kutta for 1D inviscid Burgers’ equation:
(a) $\boldsymbol{e}_{\textit{FOM}}$ vs. $\tau_{\textit{FOM}}$ in FOM; (b)
$\boldsymbol{e}_{\textit{ROM}}$ vs. $\tau_{\textit{ROM}}$ in ROM.

Table 2 summarizes the online running time of all methods, the mean time-
integrated error with respect to the FOM, and the mean time-integrated error
with respect to the Galerkin ROM using the Runge-Kutta solver. For a fair
comparison, all the ROM and FOM solutions are computed at the verified Runge-
kutta time step. As shown in the table, SVR based models, e.g. SVR2 and SVR3,
yield the smallest relative errors, however the computational cost is more
expensive than the FOM. Note that the non-intrusive ROM (VKOGA) can speed up
the solver by $6.9\times$ relative to the FOM and $3.3\times$ over the
Galerkin ROM at a relative error of $0.0353$ and $0.0164$ respectively.

Table 2: Comparison of different machine learning methods using $4$th-order Runge-Kutta for 1D inviscid Burgers’ equation. The running time of FOM and Galerkin is $0.401$s and $0.194$s respectively. Method | Online running time (s) | R Err (w.r.t. FOM) | R Err (w.r.t. Galerkin)
---|---|---|---
SVR 2 | 1.854 | 0.0309 | 2.42e-3
SVR 3 | 1.744 | 0.0310 | 2.42e-3
SVR rbf | 1.773 | 0.0321 | 5.03e-3
Random Forest | 6.513 | 0.0817 | 0.0726
Boosting | 10.076 | 0.0904 | 0.0858
kNN | 0.244 | 0.0467 | 0.0336
VKOGA | 0.058 | 0.0353 | 0.0164
SINDy | 0.063 | 0.0310 | 2.42e-3
Galerkin | 0.194 | 0.0304 | 0

###  4.2 2D convection–diffusion

\begin{overpic}[width=195.12767pt]{figures_2nd/FOM1} \put(-1.0,79.0){(a)}
\end{overpic}\begin{overpic}[width=195.12767pt]{figures_2nd/FOM2}
\put(-1.0,79.0){(b)} \end{overpic}

Figure 6: Solution profile of 2D convection–diffusion equation, $u(x,y,t=2)$
with input parameter $\mu_{1}=\mu_{2}=9.5$.

\begin{overpic}[width=411.93767pt]{figures_2nd/Fig0_RK_2D} \put(3.0,32.0){(a)}
\end{overpic}

\begin{overpic}[width=411.93767pt]{figures_2nd/Fig0_BE_2D} \put(3.0,32.0){(b)}
\end{overpic}

Figure 7: Timestep verification for 2D convection–diffusion equation: (a)
Runge-Kutta: we select a time step of $1.25e{\text{-}}1$, as it leads to an
approximated convergence rate of $4$ and an averaged approximated error of
$2e{\text{-}}5$ for the selected state. (b) backward Euler integrator: we
select a time step of $3.12e{\text{-}}2$, as it leads to an approximated
convergence rate of 1 and an approximated error of $5e{\text{-}}4$ for the
selected state.

We consider a 2D parametric nonlinear hear equation. Given a state variable
$u=u(x,y,t)$, the governing equation is described as

$\displaystyle\frac{\partial u(x,y,t)}{\partial t}$
$\displaystyle=-\mu_{0}\bigtriangledown^{2}u-\frac{\mu_{0}\mu_{1}}{\mu_{2}}\left(e^{\mu_{2}u}-1\right)+\cos(2\pi
x)\cos(2\pi y),$ (25) $\displaystyle u(x,y,0)$ $\displaystyle=0.$ (25a)

The parameters are given by $\mu_{0}=0.01$, and
$(\mu_{1},\mu_{2})\in[9,10]^{2}$. The spatial domain is $[0,1]^{2}$ and
Dirichlet boundary conditions are applied. The FOM uses a finite difference
discretization with $51\times 51$ grid points. The full time domain is $[0,2]$
and we evaluate both the backward Euler and Runge-Kutta methods for time
integration with uniform time steps. Figure 6 shows the solution profile at
$t=2$ with input parameter $(\mu_{1},\mu_{2})=(9.5,9.5)$.

####  4.2.1 Data collection

Similar to the 1D case, first we investigate the appropriate time step $\Delta
t$ for solving the ODE. We collect the solutions of a sequence number of time
steps $N_{t}=[25,50,100,200,400,800,1600,3200,6400]$ for (i) explicit Runge-
Kutta and (ii) implicit backward-Euler integrator. The verification results in
Figure 7 shows that $N_{t}=200$ is a reasonable number of time steps to use
for Runge-Kutta and $N_{t}=800$ for backward Euler method. During the offline
stage, we run four full simulations corresponding to the corner parameters of
the space $[9,10]^{2}$. Then, the training data are sampled from a Latin-
hypercube for better covering the parameter space. In the sampling,
$N_{\text{training}}$ and $N_{\text{validation}}$ instances of the state, time
and parameters are generated following the criterion that the minimum
distances between the data points are maximized. We use the default size of
the training set $N_{\text{training}}=1000$ and of the validation set
$N_{\text{validation}}=500$. Then the reduced vector field
$\boldsymbol{f}_{r}$ is computed for each input pairs
$(\hat{\boldsymbol{x}},t;\boldsymbol{\mu})$. In the training and validation
stage, we regress the reduced vector field $\boldsymbol{f}_{r}$ by the input
$(\hat{\boldsymbol{x}},t;\boldsymbol{\mu})$; in the test stage, we evaluate
the ROM. The parameters are fixed to be $(\mu_{1},\mu_{2})=(9.5,9.5)$.

####  4.2.2 Model validation

We report the performance of SVR(2nd, 3rd poly and rbf), kNN, Random Forest,
Boosting, VKOGA, and SINDy as regression models to approximate reduced
velocity. In particular, as in Section 4.1.2, for each regression method, we
change the model hyperparameters and plot the relative training and validation
error. The relative error is defined by Equation (18). Similarly, we plot the
learning curve of each regression method and compare the performance of each
model on training and validation data over a varying number of training
instances in Appendix A.2. We aim to balance bias and variance in each
regression model, by properly choosing hyperparameters and the number of
training instances.

####  4.2.3 Simulation of the surrogate ROM

We can now solve this 2D problem using the surrogate model along the
trajectory of the dynamical system. After applying time integration to the
regression-based ROM, we compute the relative error of the proposed models as
a function of time. As in Section 4.1.3, we investigate both the backward
Euler and Runge-Kutta integrators. The following are the simulation results
from backward Euler method with $N_{t}=800$. Figure 8 plots the state-space
error with respect to FOM and ROM using the backward Euler integrator. VKOGA
outperforms the other models, achieving a relative ROM error below
$6e\text{-}2$ over the entire time domain, and the accuracy is closest to
Galerkin ROM.

\begin{overpic}[width=195.12767pt]{figures_2nd/FomErrVsTime_BE}
\put(-1.0,78.0){(a)} \put(49.0,-0.1){\hbox{\pagecolor{white}{Time}}}
\put(-2.0,32.0){\hbox{\pagecolor{white}\rotatebox{90.0}{{Relative error}}}}
\end{overpic}\begin{overpic}[width=195.12767pt]{figures_2nd/RomErrVsTime_BE}
\put(-1.0,78.0){(b)} \put(49.0,-0.1){\hbox{\pagecolor{white}{Time}}}
\put(-2.0,32.0){\hbox{\pagecolor{white}\rotatebox{90.0}{{Relative error}}}}
\end{overpic}

Figure 8: Backward Euler for 2D convection–diffusion equation: time evolution
of relative error: (a) $\boldsymbol{e}_{\textit{FOM}}(t)$ in FOM; and (b)
$\boldsymbol{e}_{\textit{ROM}}(t)$ in ROM.

Figure 9 plots the Pareto frontier error as a function of the relative running
time using the backward Euler integrator. VKOGA performs best in terms of both
accuracy and time efficiency, when comparing with the Galerkin ROM.

\begin{overpic}[width=195.12767pt]{figures_2nd/Pareto_FOM_BE}
\put(-1.0,78.0){(a)}
\put(-2.0,32.0){\hbox{\pagecolor{white}\rotatebox{90.0}{{Relative error}}}}
\put(32.0,-0.1){\hbox{\pagecolor{white}{Relative time w.r.t. ROM}}}
\end{overpic}\begin{overpic}[width=195.12767pt]{figures_2nd/Pareto_ROM_BE}
\put(-1.0,78.0){(b)}
\put(-2.0,32.0){\hbox{\pagecolor{white}\rotatebox{90.0}{{Relative error}}}}
\put(32.0,-0.1){\hbox{\pagecolor{white}{Relative time w.r.t. ROM}}}
\end{overpic}

Figure 9: Pareto frontier of relative error with respect to the relative
running time using backward Euler for 2D convection–diffusion equation: (a)
$\boldsymbol{e}_{\textit{FOM}}$ vs. $\tau_{\textit{FOM}}$ in FOM; (b)
$\boldsymbol{e}_{\textit{ROM}}$ vs. $\tau_{\textit{ROM}}$ in ROM.

Table 3 represents online running time of all methods, the mean time-
integrated error with respect to FOM, and the mean time-integrated error with
respect to the Galerkin ROM using the backward Euler integrator for the 2D
convection–diffusion equation. For a fair comparison, all the ROM and FOM
solutions are computed at the verified backward Euler time step. Note that the
non-intrusive ROM, e.g. VKOGA with Newton’s method, can improve the solve time
by three orders of magnitude over the FOM and $111.2\times$ compared to the
Galerkin ROM at a relative error of $0.0059$ and $0.0041$ respectively.

Table 3: Comparison of different machine learning methods using the Backward Euler integrator for 2D convection–diffusion equation: Newton’s method (N); FP-fixed-point iteration (FP). The running time of FOM and Galerkin is $209.393$s and $16.837$s respectively. Method | Online running time (s) | R Err (w.r.t. FOM) | R Err (w.r.t. Galerkin)
---|---|---|---
SVR 2 FP | 13.510 | 0.0074 | 0.0061
SVR 2 N | 10.510 | 0.0074 | 0.0061
SVR 3 FP | 13.873 | 0.0045 | 0.0028
SVR 3 N | 10.135 | 0.0045 | 0.0028
SVR rbf FP | 13.349 | 0.0181 | 0.0169
SVR rbf N | 10.405 | 0.0181 | 0.0169
Random Forest | 120.674 | 0.0332 | 0.0326
Boosting | 190.632 | 0.0147 | 0.0135
kNN | 0.435 | 0.0149 | 0.0142
VKOGA FP | 2.505 | 0.0059 | 0.0041
VKOGA N | 0.151 | 0.0059 | 0.0041
SINDy FP | 0.534 | 0.0066 | 0.0054
SINDy N | 0.236 | 0.0066 | 0.0054
Galerkin | 0.942 | 0.0029 | 0

We examine the simulation results from the 4th-order Runge-Kutta method with
$N_{t}=200$. Figure 10 shows the state-space error with respect to FOM and ROM
using the Runge-Kutta integrator. SVR2, SVR3, VKOGA, and SINDy have a
comparable performance, and result in a smaller ROM error relative to in the
backward Euler solver. We notice that the random forest, boosting models and
kNN begin to diverge quickly in the second half of the time domain.

\begin{overpic}[width=195.12767pt]{figures_2nd/FomErrVsTime_RK4}
\put(-1.0,78.0){(a)}
\put(-2.0,32.0){\hbox{\pagecolor{white}\rotatebox{90.0}{{Relative error}}}}
\put(49.0,-0.1){\hbox{\pagecolor{white}{Time}}}
\end{overpic}\begin{overpic}[width=195.12767pt]{figures_2nd/RomErrVsTime_RK4}
\put(-1.0,78.0){(b)}
\put(-2.0,32.0){\hbox{\pagecolor{white}\rotatebox{90.0}{{Relative error}}}}
\put(49.0,-0.1){\hbox{\pagecolor{white}{Time}}} \end{overpic}

Figure 10: Runge-Kutta for 2D convection–diffusion equation: time evolution of
relative error: (a) $\boldsymbol{e}_{\textit{FOM}}(t)$ in FOM; and (b)
$\boldsymbol{e}_{\textit{ROM}}(t)$ in ROM.

Figure 11 plots the Pareto frontier error as a function of the relative
running time using the backward Euler integrator. VKOGA and SINDy outperform
the other models in terms of both computation accuracy and time cost. The
relative error compared to Galerkin ROM and FOM is below $1e{\text{-}}2$.

\begin{overpic}[width=195.12767pt]{figures_2nd/Pareto_FOM_RK4}
\put(-1.0,78.0){(a)}
\put(-2.0,32.0){\hbox{\pagecolor{white}\rotatebox{90.0}{{Relative error}}}}
\put(32.0,-0.1){\hbox{\pagecolor{white}{Relative time w.r.t. FOM}}}
\end{overpic}\begin{overpic}[width=195.12767pt]{figures_2nd/Pareto_ROM_RK4}
\put(-1.0,78.0){(b)}
\put(-2.0,32.0){\hbox{\pagecolor{white}\rotatebox{90.0}{{Relative error}}}}
\put(32.0,-0.1){\hbox{\pagecolor{white}{Relative time w.r.t. ROM}}}
\end{overpic}

Figure 11: Pareto frontier of relative error with respect to the relative
running time using 4th-order Runge-Kutta for 2D convection–diffusion equation:
(a) $\boldsymbol{e}_{\textit{FOM}}$ vs. $\tau_{\textit{FOM}}$ in FOM; (b)
$\boldsymbol{e}_{\textit{ROM}}$ vs. $\tau_{\textit{ROM}}$ in ROM.

Table 4 presents online running time of all methods, the mean time-integrated
error with respect to the FOM, and the mean time-integrated error with respect
to the Galerkin ROM using the Runge-Kutta integrator for the 2D
convection–diffusion equation. For a fair comparison, all the ROM and FOM
solutions are computed at the verified backward Euler time step. Note that the
computational efficiency of the non-intrusive ROM performs significantly
better than Galerkin ROM. VKOGA using Runge-Kutta can speed up the solver
$3406.9\times$ that of the FOM and $468.3\times$ that of the Galerkin ROM at a
relative error of $0.0083$ and $0.0069$. SINDy using Runge-Kutta can
accelerate the solver $2524.1\times$ over the FOM and $347.0\times$ compared
to the Galerkin ROM at a relative error of $0.0077$ and $0.0066$ respectively.

Table 4: Comparison of different machine learning methods using $4$th-order Runge-Kutta for 2D convection–diffusion equation. The running time of FOM and Galerkin is $138.319$s and $19.012$s respectively. Method | Online running time (s) | R Err (w.r.t. FOM) | R Err (w.r.t. Galerkin)
---|---|---|---
SVR 2 | 1.850 | 0.0087 | 0.0076
SVR 3 | 1.740 | 0.0065 | 0.0052
SVR rbf | 1.777 | 0.0185 | 0.0174
Random Forest | 20.313 | 0.0426 | 0.0419
Boosting | 32.668 | 0.0177 | 0.0166
kNN | 0.175 | 0.0153 | 0.0146
VKOGA | 0.041 | 0.0083 | 0.0069
SINDy | 0.055 | 0.0077 | 0.0066
Galerkin | 19.012 | 0.0029 | 0

##  5 Conclusions

In this study, we demonstrate the effectiveness of a data-driven model
reduction method for solving two types of parametric PDEs non-intrusively, in
both explicit and implicit time integration schemes. The approach successfully
avoids the cost and intrusiveness of querying the full-order model and
reduced-order model in the simulation, by approximating low-dimensional
operators using regression methods. In the offline stage, we train the
regression models using state-of-the art techniques from machine learning and
specific dynamical learning methods in a reduced order architecture; in the
online stage, the surrogate ROM then solves problems more efficiently (orders
of magnitude in speedup) using the learned reduced operators. The proposed
approach speeds up the full-order model simulation by one order of magnitude
in the 1D Burgers’ equation and three order of magnitude in the 2D convection-
diffusion equation. Among the machine learning models we evaluate, VKOGA and
SINDy distinguish themselves from the other models, delivering superior cost
vs. error trade-off.

Further work involves a number of important extensions and directions that
arise out of this work. First, it will be interesting to investigate the
effectiveness on solving a different types of PDE or a more complex nonlinear
dynamical systems, e.g. Euler’s equation. Morerover, for nonlinear Lagrangian
dynamical systems, we need to develop structure-preserving approximations for
all reduced Lagrangian ingredients in the model reduction. Rather than apply a
Galerkin projection to obtain the ROM, one can alternatively employ a least-
square Petrov-Galerkin (LSPG) (Carlberg et al., 2011, 2013, 2017), which
requires a regression method predicting non-negative values. The growing
intersection of dynamical systems and data science are driving innovations in
estimating the prediction of complex systems (Brunton and Kutz, 2019). With a
deeper understanding of neural networks, it may be possible to generalize the
non-intrusive approach to another level of accuracy. Furthermore, physics-
informed models (Raissi et al., 2019) may accelerate learning from data and
add interpretability to existing approaches. This framework can also be
combined with distributed systems to enable parallelism for extreme-scale high
performance computing.

## Acknowledgments

We would like to thank Kevin Carlberg, Steve Brunton and Youngsoo Choi for
valuable discussions. ZB acknowledges support by the U.S. Department of
Energy, Office of Science, under Award Number DE-AC02-05CH11231. ZB and LP
acknowledge support for preliminary work from Sandia National Laboratories.
Any subjective views or opinions that might be expressed in the paper do not
necessarily represent the views of the U.S. Department of Energy or the United
States Government. Sandia National Laboratories is a multimission laboratory
managed and operated by National Technology & Engineering Solutions of Sandia,
LLC, a wholly owned subsidiary of Honeywell International Inc., for the U.S.
Department of Energy’s National Nuclear Security Administration under contract
DE-NA0003525.

## References

  * Lee et al. (2013) Myoungkyu Lee, Nicholas Malaya, and Robert D Moser.  Petascale direct numerical simulation of turbulent channel flow on up to 786k cores.  In _Proceedings of SC13: International Conference for High Performance Computing, Networking, Storage and Analysis_ , page 61. ACM, 2013.
  * Brunton and Noack (2015) S. L. Brunton and B. R. Noack.  Closed-loop turbulence control: Progress and challenges.  _Applied Mechanics Reviews_ , 67:050801–1–050801–48, 2015\.
  * Aubry et al. (1988) Nadine Aubry, Philip Holmes, John L Lumley, and Emily Stone.  The dynamics of coherent structures in the wall region of a turbulent boundary layer.  _Journal of Fluid Mechanics_ , 192:115–173, 1988.
  * Berkooz et al. (1993) 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.
  * Holmes et al. (1996) P. J. Holmes, J. L. Lumley, and G. Berkooz.  _Turbulence, coherent structures, dynamical systems and symmetry_.  Cambridge Monographs in Mechanics. Cambridge University Press, Cambridge, England, 1996.
  * Brunton and Kutz (2018) S. L. Brunton and J. N. Kutz.  _Data-Driven Science and Engineering: Machine Learning, Dynamical Systems, and Control_.  Cambridge University Press, 2018.
  * Noack et al. (2003) B. R. Noack, K. Afanasiev, M. Morzynski, G. Tadmor, and F. Thiele.  A hierarchy of low-dimensional models for the transient and post-transient cylinder wake.  _Journal of Fluid Mechanics_ , 497:335–363, 2003.
  * Amsallem et al. (2009) David Amsallem, Julien Cortial, Kevin Carlberg, and Charbel Farhat.  A method for interpolating on manifolds structural dynamics reduced-order models.  _International Journal for Numerical Methods in Engineering_ , 80(9):1241–1258, Nov 2009.  ISSN 1097-0207.  doi: 10.1002/nme.2681.  URL http://dx.doi.org/10.1002/nme.2681.
  * Carlberg et al. (2013) Kevin Carlberg, Charbel Farhat, Julien Cortial, and David Amsallem.  The gnat method for nonlinear model reduction: effective implementation and application to computational fluid dynamics and turbulent flows.  _Journal of Computational Physics_ , 242:623–647, 2013\.
  * Carlberg et al. (2015) Kevin Carlberg, Ray Tuminaro, and Paul Boggs.  Preserving lagrangian structure in nonlinear model reduction with application to structural dynamics.  _SIAM Journal on Scientific Computing_ , 37(2):B153–B184, 2015.
  * Carlberg et al. (2017) Kevin Carlberg, Matthew Barone, and Harbir Antil.  Galerkin v. least-squares petrov–galerkin projection in nonlinear model reduction.  _Journal of Computational Physics_ , 330:693–734, 2017\.
  * Barrault et al. (2004) Maxime Barrault, Yvon Maday, Ngoc Cuong Nguyen, and Anthony T Patera.  An ‘empirical interpolation’method: application to efficient reduced-basis discretization of partial differential equations.  _Comptes Rendus Mathematique_ , 339(9):667–672, 2004.
  * Chaturantabut and Sorensen (2010) Saifon Chaturantabut and Danny C Sorensen.  Nonlinear model reduction via discrete empirical interpolation.  _SIAM Journal on Scientific Computing_ , 32(5):2737–2764, 2010.
  * Ljung (1999) L. Ljung.  _System Identification: Theory for the User_.  Prentice Hall, 1999.
  * Juang (1994) J. N. Juang.  _Applied System Identification_.  Prentice Hall PTR, Upper Saddle River, New Jersey, 1994.
  * Nelles (2013) Oliver Nelles.  _Nonlinear system identification: from classical approaches to neural networks and fuzzy models_.  Springer Science & Business Media, 2013.
  * Juang and Pappa (1985) J. N. Juang and R. S. Pappa.  An eigensystem realization algorithm for modal parameter identification and model reduction.  _Journal of Guidance, Control, and Dynamics_ , 8(5):620–627, 1985.
  * Schmid (2010) P. J. Schmid.  Dynamic mode decomposition of numerical and experimental data.  _Journal of Fluid Mechanics_ , 656:5–28, August 2010.  ISSN 0022-1120.
  * Tu et al. (2014) J. H. Tu, C. W. Rowley, D. M. Luchtenburg, S. L. Brunton, and J. N. Kutz.  On dynamic mode decomposition: theory and applications.  _Journal of Computational Dynamics_ , 1(2):391–421, 2014.
  * Kutz et al. (2016) J. N. Kutz, S. L. Brunton, B. W. Brunton, and J. L. Proctor.  _Dynamic Mode Decomposition: Data-Driven Modeling of Complex Systems_.  SIAM, 2016.
  * Bai et al. (2020) Zhe Bai, Eurika Kaiser, Joshua L Proctor, J Nathan Kutz, and Steven L Brunton.  Dynamic mode decomposition for compressive system identification.  _AIAA Journal_ , 58(2):561–574, 2020.
  * Mezić (2005) Igor Mezić.  Spectral properties of dynamical systems, model reduction and decompositions.  _Nonlinear Dynamics_ , 41(1-3):309–325, 2005\.
  * Rowley et al. (2009) C. W. Rowley, I. Mezić, S. Bagheri, P. Schlatter, and D.S. Henningson.  Spectral analysis of nonlinear flows.  _J. Fluid Mech._ , 645:115–127, 2009.
  * Brunton et al. (2017) S. L. Brunton, B. W. Brunton, J. L. Proctor, E. Kaiser, and J. N. Kutz.  Chaos as an intermittently forced linear system.  _Nature Communications_ , 8(19):1–9, 2017.
  * Kukreja and Brenner (2007) Sunil L Kukreja and Martin J Brenner.  Nonlinear system identification of aeroelastic systems: A structure-detection approach.  In _Identification and Control_ , pages 117–145. Springer, 2007.
  * Billings (2013) Stephen A Billings.  _Nonlinear system identification: NARMAX methods in the time, frequency, and spatio-temporal domains_.  John Wiley & Sons, 2013.
  * Semeraro et al. (2016) Onofrio Semeraro, Francois Lusseyran, Luc Pastur, and Peter Jordan.  Qualitative dynamics of wavepackets in turbulent jets.  _arXiv preprint arXiv:1608.06750_ , 2016.
  * Milano and Koumoutsakos (2002) Michele Milano and Petros Koumoutsakos.  Neural network modeling for near wall turbulent flow.  _Journal of Computational Physics_ , 182(1):1–26, 2002.
  * Ling et al. (2016) Julia Ling, Andrew Kurzawski, and Jeremy Templeton.  Reynolds averaged turbulence modelling using deep neural networks with embedded invariance.  _Journal of Fluid Mechanics_ , 807:155–166, 2016.
  * Zhang and Duraisamy (2015) Ze Jia Zhang and Karthikeyan Duraisamy.  Machine learning methods for data-driven turbulence modeling.  In _22nd AIAA Computational Fluid Dynamics Conference_ , page 2460\. AIAA, 2015.
  * Singh et al. (2017) Anand Pratap Singh, Shivaji Medida, and Karthik Duraisamy.  Machine-learning-augmented predictive modeling of turbulent separated flows over airfoils.  _AIAA Journal_ , 2017.
  * Duraisamy et al. (2018) Karthik Duraisamy, Gianluca Iaccarino, and Heng Xiao.  Turbulence modeling in the age of data.  _To appear in Annual Reviews of Fluid Mechanics (arXiv preprint arXiv:1804.00183)_ , 2018.
  * Wang et al. (2017) Jian-Xun Wang, Jin-Long Wu, and Heng Xiao.  Physics-informed machine learning approach for reconstructing reynolds stress modeling discrepancies based on dns data.  _Physical Review Fluids_ , 2(3):034603, 2017.
  * Kutz (2017) J Nathan Kutz.  Deep learning in fluid dynamics.  _Journal of Fluid Mechanics_ , 814:1–4, 2017.
  * Bongard and Lipson (2007) Josh Bongard and Hod Lipson.  Automated reverse engineering of nonlinear dynamical systems.  _Proceedings of the National Academy of Sciences_ , 104(24):9943–9948, 2007.
  * Schmidt and Lipson (2009) Michael Schmidt and Hod Lipson.  Distilling free-form natural laws from experimental data.  _Science_ , 324(5923):81–85, 2009.
  * Brunton et al. (2016) S. L. Brunton, J. L. Proctor, and J. N. Kutz.  Discovering governing equations from data by sparse identification of nonlinear dynamical systems.  _Proceedings of the National Academy of Sciences_ , 113(15):3932–3937, 2016.
  * Raissi and Karniadakis (2018) Maziar Raissi and George E. Karniadakis.  Hidden physics models: Machine learning of nonlinear partial differential equations.  _Journal of Computational Physics_ , 357:125–141, 2018\.
  * Kaiser et al. (2014) E. Kaiser, B. R. Noack, L. Cordier, A. Spohn, M. Segond, M. Abel, G. Daviller, J. Osth, S. Krajnovic, and R. K. Niven.  Cluster-based reduced-order modelling of a mixing layer.  _J. Fluid Mech._ , 754:365–414, 2014.
  * Wan et al. (2018) Zhong Yi Wan, Pantelis Vlachas, Petros Koumoutsakos, and Themistoklis Sapsis.  Data-assisted reduced-order modeling of extreme events in complex dynamical systems.  _PloS one_ , 13(5):e0197704, 2018.
  * Vlachas et al. (2018) Pantelis R Vlachas, Wonmin Byeon, Zhong Y Wan, Themistoklis P Sapsis, and Petros Koumoutsakos.  Data-driven forecasting of high-dimensional chaotic systems with long short-term memory networks.  _Proc. R. Soc. A_ , 474(2213):20170844, 2018.
  * Raissi and Karniadakis (2017) Maziar Raissi and George Em Karniadakis.  Machine learning of linear differential equations using gaussian processes.  _arXiv preprint arXiv:1701.02440_ , 2017.
  * Wan and Sapsis (2017) Zhong Yi Wan and Themistoklis P Sapsis.  Reduced-space Gaussian process regression for data-driven probabilistic forecast of chaotic dynamical systems.  _Physica D: Nonlinear Phenomena_ , 345:40–55, 2017.
  * Rudy et al. (2017) S. H. Rudy, S. L. Brunton, J. L. Proctor, and J. N. Kutz.  Data-driven discovery of partial differential equations.  _Science Advances_ , 3(e1602614), 2017.
  * Loiseau and Brunton (2018) J.-C. Loiseau and S. L. Brunton.  Constrained sparse Galerkin regression.  _Journal of Fluid Mechanics_ , 838:42–67, 2018.
  * Loiseau et al. (2018) J.-C. Loiseau, B. R. Noack, and S. L. Brunton.  Sparse reduced-order modeling: sensor-based dynamics to full-state estimation.  _Journal of Fluid Mechanics_ , 844:459–490, 2018.
  * Schaeffer (2017) Hayden Schaeffer.  Learning partial differential equations via data discovery and sparse optimization.  _Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences_ , 473(2197):20160446, 2017.
  * Swischuk et al. (2018) Renee Swischuk, Laura Mainini, Benjamin Peherstorfer, and Karen Willcox.  Projection-based model reduction: Formulations for physics-based machine learning.  _Computers & Fluids_, 2018.  ISSN 0045-7930.  URL http://www.sciencedirect.com/science/article/pii/S0045793018304250.
  * Peherstorfer and Willcox (2016) Benjamin Peherstorfer and Karen Willcox.  Data-driven operator inference for nonintrusive projection-based model reduction.  _Computer Methods in Applied Mechanics and Engineering_ , 306:196 – 215, 2016.  ISSN 0045-7825.  URL http://www.sciencedirect.com/science/article/pii/S0045782516301104.
  * Carlberg et al. (2019) Kevin T Carlberg, Antony Jameson, Mykel J Kochenderfer, Jeremy Morton, Liqian Peng, and Freddie D Witherden.  Recovering missing cfd data for high-order discretizations using deep neural networks and dynamics learning.  _Journal of Computational Physics_ , 395:105–124, 2019\.
  * Wang et al. (2016) M. Wang, H. Li, X. Chen, and Y. Chen.  Deep learning-based model reduction for distributed parameter systems.  _IEEE Transactions on Systems, Man, and Cybernetics: Systems_ , 46(12):1664–1674, Dec 2016.  ISSN 2168-2216.
  * Kani and Elsheikh (2017) J.Nagoor Kani and Ahmed H. Elsheikh.  Dr-rnn: A deep residual recurrent neural network for model reduction.  _Arxiv_ , 09 2017.  URL https://arxiv.org/pdf/1709.00939.
  * Xie et al. (2019) Xuping Xie, Guannan Zhang, and Clayton G. Webster.  Non-intrusive inference reduced order model for fluids using deep multistep neural network.  _Mathematics_ , 7(8), 2019.  ISSN 2227-7390.  doi: 10.3390/math7080757.  URL https://www.mdpi.com/2227-7390/7/8/757.
  * Hesthaven and Ubbiali (2018) J.S. Hesthaven and S. Ubbiali.  Non-intrusive reduced order modeling of nonlinear problems using neural networks.  _Journal of Computational Physics_ , 363:55 – 78, 2018\.  ISSN 0021-9991.  doi: https://doi.org/10.1016/j.jcp.2018.02.037.  URL http://www.sciencedirect.com/science/article/pii/S0021999118301190.
  * Wang et al. (2018) Qian Wang, Jan S. Hesthaven, and Deep Ray.  Non-intrusive reduced order modeling of unsteady flows using artificial neural networks with application to a combustion problem.  _Journal of Computational Physics_ , 2018.  URL http://infoscience.epfl.ch/record/255708.
  * Xiao et al. (2015) D. Xiao, F. Fang, C. Pain, and G. Hu.  Non-intrusive reduced-order modelling of the navier-stokes equations based on rbf interpolation.  _International Journal for Numerical Methods in Fluids_ , 79(11):580–595, Jul 2015.  ISSN 0271-2091.
  * Vaerenbergh (2010) Steven Van Vaerenbergh.  Kernel methods for nonlinear identification, equalization and separation of signals kernel methods for nonlinear identification, equalization and separation of signals.  _Kernel methods for nonlinear identification, equalization and separation of signals_ , 2010.
  * Wirtz et al. (2015) D Wirtz, N Karajan, and B Haasdonk.  Surrogate modeling of multiscale models using kernel methods.  _International Journal for Numerical Methods in Engineering_ , 101(1):1–28, 2015.
  * Xie et al. (2018) X. Xie, M. Mohebujjaman, L. Rebholz, and T. Iliescu.  Data-driven filtered reduced order modeling of fluid flows.  _SIAM Journal on Scientific Computing_ , 40(3):B834–B857, 2018.  doi: 10.1137/17M1145136.  URL https://doi.org/10.1137/17M1145136.
  * Altman (1992) Naomi S Altman.  An introduction to kernel and nearest-neighbor nonparametric regression.  _The American Statistician_ , 46(3):175–185, 1992\.
  * Smola and Schölkopf (2004) Alex J Smola and Bernhard Schölkopf.  A tutorial on support vector regression.  _Statistics and computing_ , 14(3):199–222, 2004\.
  * Breiman (2001) Leo Breiman.  Random forests.  _Machine learning_ , 45(1):5–32, 2001.
  * Hastie et al. (2009) Trevor Hastie, Robert Tibshirani, and Jerome Friedman.  The elements of statistical learnin.  _Cited on_ , page 33, 2009.
  * Wirtz and Haasdonk (2013) Daniel Wirtz and Bernard Haasdonk.  A vectorial kernel orthogonal greedy algorithm.  _Dolomites Research Notes on Approximation_ , 6(Special_Issue), 2013.
  * Carlberg et al. (2011) K. Carlberg, C. Farhat, and C. Bou-Mosleh.  Efficient nonlinear model reduction via a least-squares Petrov-Galerkin projection and compressive tensor approximations.  _International Journal for Numerical Methods in Engineering,_in press__ , 86(2):155–181, October 2011.
  * Brunton and Kutz (2019) Steven L Brunton and J Nathan Kutz.  _Data-driven science and engineering: Machine learning, dynamical systems, and control_.  Cambridge University Press, 2019.
  * Raissi et al. (2019) Maziar Raissi, Paris Perdikaris, and George E Karniadakis.  Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations.  _Journal of Computational Physics_ , 378:686–707, 2019\.
  * Tibshirani (1996) Robert Tibshirani.  Regression shrinkage and selection via the lasso.  _Journal of the Royal Statistical Society: Series B (Methodological)_ , 58(1):267–288, 1996.

##  Appendix A Regression models and tuning

###  A.1 Special regression models

In this section, we detail two non-transitional machine learning models, VKOGA
and SINDy, in addition to the SVR, kNN, random forest and boosting.

####  A.1.1 Vectorial kernel orthogonal greedy algorithm (VKOGA)

Developed from kernel-based methods, VKOGA makes the assumption that a common
subspace can be used for every component of vector-valued response. With such,
it significantly reduces the training and evaluation time when the number of
kernels in the resulting regression model is large. VKOGA constructs the
regression model

$\boldsymbol{h}(\boldsymbol{x})=\sum_{i=1}^{n_{\text{VKOGA}}}{\boldsymbol{\alpha}_{i}}K(\boldsymbol{x}_{i},\boldsymbol{x}),$
(A.1)

where $K:\mathbb{R}^{N_{x}}\times\mathbb{R}^{N_{x}}\rightarrow\mathbb{R}$
denotes a kernel function, for each individual regression model $h_{i}$,
$i=1,\cdots,N_{x}$, and
$\alpha_{i}\in\mathbb{R}^{\boldsymbol{x}_{x}},i=1,\cdots,N_{x}$ are vector-
valued basis functions. VKOGA uses a greedy algorithm to compute kernel
functions $K(\bar{\boldsymbol{x}},\cdot)$. The greedy algorithm determines
kernel centers from
$\Omega={\bar{\boldsymbol{x}}^{1},\cdots,\bar{\boldsymbol{x}}^{n_{train}}}$.
Let the initial state $\Omega_{0}=\emptyset$, and then at stage $m$, choose

$\boldsymbol{x}_{m}=\operatorname*{argmax}_{\boldsymbol{x}\in\Omega\backslash\Omega_{m-1}}|\langle\boldsymbol{h},\phi_{x}^{m-1}\rangle|,$
(A.2)

where $\phi_{x}^{m-1}$ denotes the orthogonal remainder of
$K(\boldsymbol{x},\cdot)$ with respect to the reproducing kernel Hilbert space
spanned by
${K(\boldsymbol{x}_{1},\cdot)},\cdots,K(\boldsymbol{x}_{m-1},\cdot)$. Then the
kernel centers are updated $\Omega_{m}=\Omega_{m-1}\cup{\boldsymbol{x}_{m}}$.
At the next step, the basis function $\boldsymbol{\alpha}_{i}$,
$i=1,\cdots,n_{\text{VKOGA}}$ are determined by a least-squares approximation
to fit the training data. In the numerical experiments, we apply the Gaussian
RBF kernel
$K(\bar{\boldsymbol{x}}^{i},\bar{\boldsymbol{x}}^{j})=exp(-\gamma\|\bar{\boldsymbol{x}}^{i},\bar{\boldsymbol{x}}^{j}\|)$,
where the hyperparameter is the number of kernel functions used.

####  A.1.2 Sparse identification of nonlinear dynamics (SINDy)

The sparse identification of nonlinear dynamics (SINDy) method [Brunton et
al., 2016] identifies nonlinear dynamical system from measurement data,
seeking to approximate the vector field $\boldsymbol{f}$ by a generalized
linear model. This architecture avoids overfitting by identifying a
parsimonious model, which is more explainable or predictable than black-box
machine learning. Although SINDy was originally devised for continuous-time
dynamics, it can be extended to discrete-time systems. In particular, SINDy
yields a model

$\boldsymbol{h}(\boldsymbol{x},\theta)=\sum_{i=1}^{n_{\text{{SINDy}}}}p_{i}(\boldsymbol{x})\theta_{i},$
(A.3)

where $p_{i}:\mathbb{R}^{N_{x}}\rightarrow\mathbb{R}$,
$i=1,\cdots,n_{\text{SINDy}}$ denotes a "library" of candidate basis
functions. The relevant terms that are active in the dynamics are solved using
sparse regression that penalizes the number of functions in the library
$\boldsymbol{P}(\boldsymbol{x})$:

$\theta_{i}=\operatorname*{argmin}_{\theta^{\prime}_{i}}\|\boldsymbol{h}-\boldsymbol{P}(\boldsymbol{x})\theta^{\prime}_{i}\|_{2}+\alpha\|\theta^{\prime}_{i}\|_{1},$
(A.4)

where $\|\cdot\|_{1}$ promotes sparsity in the coefficient vector $\theta_{i}$
and the parameter $\alpha$ balances low model complexity with accuracy. The
least absolute shrinkage and selection operator (LASSO) [Tibshirani, 1996] or
the sequential threshold least-squares (STLS) [Brunton et al., 2016] can be
used to determine a sparse set of coefficients
($\theta_{1},\cdots,\theta_{n_{\text{SINDy}}}$). In the numerical experiments,
we consider linear and quadratic functions in the library, i.e.,
$p_{i}\in\\{x_{1},\cdots,x_{N_{x}},x_{1}x_{1},x_{1}x_{2},\cdots,x_{N_{x}}x_{N_{x}}\\}$.
Hyperparameters in this approach consist of the prescribed basis functions.

A.1 summarizes the computational complexities of all the considered regression
models. Note that the hyperparameter varies in different models, i.e.
$N_{\text{trees}}$ denotes the number of decision tress in random forest,
$N_{\text{learners}}$ denotes the number of weak learners in the boosting
method, $N_{\text{functions}}$ is the numner of kernel functions in VKOGA, and
$N_{\text{bases}}$ represents the number of bases in SINDy with 2nd order
polynomials.

Table A.1: Theoretical complexity per time step, where $n$ represents the subspace dimension, $p$ represents the parameter dimension, and $N_{\text{training}}$ represents the number of training points. Method | Complexity | Comments
---|---|---
SVR 2 | O($(n+p)nN_{\text{training}}$) |
SVR 3 | O($(n+p)nN_{\text{training}}$) |
SVR rbf | O($(n+p)nN_{\text{training}}$) |
Random Forest | O($nN_{\text{trees}}N_{\text{training}}\log(N_{\text{training}})$) |  $N_{\text{trees}}$: number of decision trees
Boosting | O($nN_{\text{learners}}$) |  $N_{\text{learners}}$: number of weak learners
kNN | O($(n+p)nN_{\text{training}}+KnN_{\text{training}}$) |  $K$: number of nearest neighbors
VKOGA | O($(n+p)N_{\text{functions}}$) |  $N_{\text{functions}}$: number of kernel functions
SINDy2 | O($nN_{\text{bases}}$) |  $N_{\text{bases}}$: number of bases, $N_{\text{bases}}<n^{2}$

###  A.2 Cross-validation and hyperparameter tuning

For each regression model proposed in Section 3.2, we use cross-validation to
prevent overfitting and select hyperparameters.

\begin{overpic}[width=429.28616pt]{AppendixFig3.pdf} \end{overpic} Figure A.1:
Hyperparameter selection for all considered regression models. Relative mean
squared errors are reported: (a) SVR2, (b) SVR3, (c) SVR rbf, (d) random
forest, (e) boosting, (f) kNN, (g) VKOGA. The training and validation set
correspond to the default size: $N_{\text{training}}=1000$ and
$N_{\text{validation}}=500$. Blue curves represent training errors and red
curves represent validation errors.
\begin{overpic}[width=429.28616pt]{AppendixFig4.pdf} \end{overpic} Figure A.2:
Relative mean squared error vs. training size: (a) SVR2, $\epsilon=10^{-4}$,
(b) SVR3, $\epsilon=10^{-4}$, (c) SVR rbf, $\epsilon=10^{-4}$, (d) random
forest, (e) boosting, (f) kNN, (g) VKOGA. Note that the validation set is
fixed, $N_{\text{validation}}=500$. Blue curves represent training errors and
red curves represent validation errors.

We report the hyperparameter selection of each model for the 2D convection-
diffusion equation in Figure A.1. Figure A.1(a) shows that with the
sensitivity margin $\epsilon\leq 0.1$, both the training error and validation
errors are very small, and thus we select the sensitivity margin
$\epsilon=0.1$ for SVR2. Similarly, $\epsilon=10^{-3}$ is selected for SVRrbf.
Figure A.1(d) shows that within 40 trees, the validation error of the random
forest model is still relatively large. Moreover, there is a big gap between
the training error and test error that implies the random forest model may
suffer from overfitting for the given data. We select $N_{\text{trees}}=15$ in
this study. As shown in Figure A.1(e), the validation error in Boosting
approaches the minimum at $N_{\text{learners}}=40$, and thus we select $40$
weak learners. Figure A.1(f) shows that by choosing the number of neighbors
$K\approx 4$, the validation error is minimized. Thus, $K=4$ is a good
hyperparameter to balance bias and variance. In Figure A.1(g) we observe that
the more kernel functions we use, the smaller value of training and validation
error we obtain. In this study, we select 500 kernel functions in VKOGA.

We examine the performance of the models as an increasing size of the training
sample in Figure A.2. Figure A.2(a) shows that with $200$ training instances,
the SVR2 model has very small validation error, while SVR3 model requires a
bigger training set. For SVRrbf, the more training instances we have, the
smaller validation error we obtain. This is because the resolution of the
radial basis function is determined by the density of input parameter space.
Figure A.2(d) shows that the more data are input, the less the relative error.
Intuitively, more data give higher density in the input parameters space; and
thus each leaf of a decision tree provides finer coverage for test data.
Similarly, for the kNN model as in Figure A.2(f), close neighbors give better
representation of test data. As a result, more data can effectively reduce the
variance and solve the overfitting issue. Figure A.2(e) shows that with
approximately $1000$ training samples, the training error and validation error
begin to approach $2e\text{-}3$. This implies that a relative large number of
data are required in order to constrain the variance of the boosting method.
Figure A.2(g) shows that with $1000$ training instances, the training and
validation errors begin to converge for VKOGA, while SINDy (2nd order
polynomials) only requires a small number of training instances to get a very
accurate prediction, as seen in A.2(h). Similarly, we select the following
hyperparameters: $\epsilon=10^{-4}$ for SVR2, $\epsilon=10^{-5}$ for SVR3,
$\epsilon=10^{-3}$ for SVRrbf, $N_{\text{trees}}=15$ for random forest,
$N_{\text{learners}}=40$ for boosting, $K=6$ for $k$-nearest neighbours, and
$500$ kernel functions for VKOGA in the 1D Burgers’ implementation.

\begin{overpic}[width=197.29788pt]{figures_1st/barML} \put(-0.3,77.0){(a)}
\end{overpic}\begin{overpic}[width=195.12767pt]{figures_2nd/plotBar_ML}
\put(-0.4,78.0){(b)} \end{overpic}

Figure A.3: Relative errors with $N_{\text{training}}=1000$ and
$N_{\text{validation}}=500$: (a) 1D Burgers’ equation; and (b) 2D convection-
diffusion equation. Blue bars represent training errors and red bars represent
validation errors.

Figure A.3 summarizes the overall performance of the above models for the 1D
Burgers’ and 2D convection-diffusion equation. All the proposed models are
validated using the selected hyperparameters. We notice that SVR with the
kernel of 3nd order polynomials outperforms the other models for the 1D
Burgers’ equation. This can be related to the discretization structure of
quadratic nonlinearties in the PDE solver. VKOGA has the a better
reperformance than RF, Boosting, kNN, with a relative error less than
$1e{\text{-}}4$. SINDy reaches a relative error of $3e{\text{-}}8$. For the 2D
convection-diffusion problem, SVR models again perform better than RF,
Boosting and kNN, and SVR3 outperforms the other two kernel functions. VKOGA
has the best accuracy with relative error $\approx 1e{\text{-}}5$, while SINDy
reaches a relative error of $1e{\text{-}}3$. Note that to be consistent, we
select fixed training and validation sample sizes for all the considered
regression models. However, as illustrated in Figure A.2(h), SINDy can
approach the same level of accuracy with a small number of training and
validation samples.