# A Framework for Machine Learning of
Model Error in Dynamical Systems

Matthew E. Levine  Dept. of Computing and Mathematical Sciences
California Institute of Technology
Pasadena, CA 91125, USA<EMAIL_ADDRESS>and  Andrew M. Stuart  Dept. of
Computing and Mathematical Sciences
California Institute of Technology
Pasadena, CA 91125, USA<EMAIL_ADDRESS>
(Date: July 13, 2021.)

###### Abstract.

The development of data-informed predictive models for dynamical systems is of
widespread interest in many disciplines. We present a unifying framework for
blending mechanistic and machine-learning approaches to identify dynamical
systems from noisily and partially observed data. We compare pure data-driven
learning with hybrid models which incorporate imperfect domain knowledge,
referring to the discrepancy between an assumed truth model and the imperfect
mechanistic model as model error. Our formulation is agnostic to the chosen
machine learning model, is presented in both continuous- and discrete-time
settings, and is compatible both with model errors that exhibit substantial
memory and errors that are memoryless.

First, we study memoryless linear (w.r.t. parametric-dependence) model error
from a learning theory perspective, defining excess risk and generalization
error. For ergodic continuous-time systems, we prove that both excess risk and
generalization error are bounded above by terms that diminish with the square-
root of $T$, the time-interval over which training data is specified.

Secondly, we study scenarios that benefit from modeling with memory, proving
universal approximation theorems for two classes of continuous-time recurrent
neural networks (RNNs): both can learn memory-dependent model error, assuming
that it is governed by a finite-dimensional hidden variable and that,
together, the observed and hidden variables form a continuous-time Markovian
system. In addition, we connect one class of RNNs to reservoir computing,
thereby relating learning of memory-dependent error to recent work on
supervised learning between Banach spaces using random features.

Numerical results are presented (Lorenz ’63, Lorenz ’96 Multiscale systems) to
compare purely data-driven and hybrid approaches, finding hybrid methods less
data-hungry and more parametrically efficient. We also find that, while a
continuous-time framing allows for robustness to irregular sampling and
desirable domain-interpretability, a discrete-time framing can provide similar
or better predictive performance, especially when data are undersampled and
the vector field defining the true dynamics cannot be identified. Finally, we
demonstrate numerically how data assimilation can be leveraged to learn hidden
dynamics from noisy, partially-observed data, and illustrate challenges in
representing memory by this approach, and in the training of such models.

###### Key words and phrases:

Dynamical Systems, Model Error, Statistical Learning, Random Features,
Recurrent Neural Networks, Reservoir Computing.

###### 2020 Mathematics Subject Classification:

Primary 68T30, 37A30, 37M10; Secondary 37M25, 41A30

The authors are grateful to David Albers, Oliver Dunbar, Ian Melbourne, Fei
Lu, Ivan D. Jimenez Rodriguez, and Yisong Yue for helpful discussions. The
work of MEL and AMS was supported by NIH RO1 LM012734 “Mechanistic Machine
Learning”. MEL is also supported by the National Science Foundation Graduate
Research Fellowship under Grant No. DGE‐1745301. AMS is also supported by NSF
(award AGS-1835860), NSF (award DMS-1818977), the Office of Naval Research
(award N00014-17-1-2079), and the AFOSR under MURI award number
FA9550-20-1-0358 (Machine Learning and Physics-Based Modeling and Simulation).
The computations presented here were conducted in the Resnick High Performance
Center, a facility supported by Resnick Sustainability Institute at the
California Institute of Technology.

††copyright: ©:

##  1\. Introduction

###  1.1. Background and Literature Review

The modeling and prediction of dynamical systems and time-series is an
important goal in numerous domains, including biomedicine, climatology,
robotics, and the social sciences. Traditional approaches to modeling these
systems appeal to careful study of their mechanisms, and the design of
targeted equations to represent them. These carefully built mechanistic models
have impacted humankind in numerous arenas, including our ability to land
spacecraft on celestial bodies, provide high-fidelity numerical weather
prediction, and artificially regulate physiologic processes, through the use
of pacemakers and artificial pancreases, for example. This paper focuses on
the learning of model error: we assume that an imperfect mechanistic model is
known, and that data are used to improve it. We introduce a framework for this
problem, focusing on distinctions between Markovian and non-Markovian model
error, providing a unifying review of relevant literature, developing some
underpinning theory related to both the Markovian and non-Markovian settings,
and presenting numerical experiments which illustrate our key findings.

To set our work in context, we first review the use of data-driven methods for
time-dependent problems, organizing the literature review around four themes
comprising Sections 1.1.1, 1.1.2, 1.1.3 and 1.1.5; these are devoted,
respectively, to pure data-driven methods, hybrid methods that build on
mechanistic models, non-Markovian models that describe memory, and
applications of the various approaches. Having set the work in context, in
Section 1.2 we detail the contributions we make, and describe the organization
of the paper.

####  1.1.1. Data-Driven Modeling of Dynamical Systems

A recent wave of machine learning successes in data-driven modeling,
especially in imaging sciences, has shown that we can demand even more from
existing models, or that we can design models of more complex phenomena than
heretofore. Traditional models built from, for example, low order polynomials
and/or linearized model reductions, may appear limited when compared to the
flexible function approximation frameworks provided by neural networks and
kernel methods. Neural networks, for example, have a long history of success
in modeling dynamical systems [Narendra and Parthasarathy, 1992, González-
García et al., 1998, Krischer et al., 1993, Rico-Martinez et al., 1994, Rico-
Martines et al., 1993, Rico-Martínez et al., 1992, Lagaris et al., 1998] and
recent developments in deep learning for operators continue to propel this
trend [Lu et al., 2020, Bhattacharya et al., 2020, Li et al., 2021b, a].

The success of neural networks arguably relies on balanced expressivity and
generalizability, but other methods also excel in learning parsimonious and
generalizable representations of dynamics. A particularly popular methodology
is to perform sparse regression over a dictionary of vector fields, including
the use of thresholding approaches (SINDy) [Brunton et al., 2016] and
$L_{1}$-regularized polynomial regression [Tran and Ward, 2017, Schaeffer et
al., 2017, 2018, 2020]. Non-parametric methods, like Gaussian process models
[Rasmussen and Williams, 2006], have also been used widely for modeling
nonlinear dynamics [Wang et al., 2005, Frigola et al., 2014, Kocijan, 2016,
Chkrebtii et al., 2016]. While a good choice of kernel is often essential for
the success of these methods, recent progress has been made towards automatic
hyperparameter tuning via parametric kernel flows [Hamzi and Owhadi, 2021].
Successes with Gaussian process models were also extended to high dimensional
problems by using random feature map approximations [Rahimi and Recht, 2008a]
within the context of data-driven learning of parametric partial differential
equations (PDEs) and solution operators [Nelsen and Stuart, 2020].
Advancements to data-driven methods based on Koopman operator theory and
dynamic mode decomposition also offer exciting new possibilities for
predicting nonlinear dynamics from data [Tu et al., 2014, Korda et al., 2020,
Alexander and Giannakis, 2020].

It is important to consider whether to model in discrete- or continuous-time,
as both have potential advantages. The primary positive for continuous-time
modeling lies in its flexibility and interpretability. In particular,
continuous-time approaches are more readily and naturally applied to
irregularly sampled timeseries data, e.g. electronic health record data
[Rubanova et al., 2019], than discrete-time methods. Furthermore, this
flexibility with respect to timestep enables simple transferability of a model
learnt from discrete-time data at one timestep, to new settings with a
different timestep and indeed to variable timestep settings; the learned
right-hand-side can be used to generate numerical solutions at any timestep.
On the other hand, applying a discrete-time model to a new timestep either
requires exact alignment of subsampled data or some post-processing
interpolation step. Continuous-time models may also provide greater
interpretability than discrete-time methods when the right-hand-side of the
ordinary differential equation (ODE) is a more physically interpretable object
than the $\Delta t$-solution operator (e.g. for equation discovery, [Kaheman
et al., 2019]).

Traditional implementations of continuous-time learning require accurate
estimation of time-derivatives of the state, but this may be circumvented
using approaches that leverage autodifferentiation software [Ouala et al.,
2020, Rubanova et al., 2019, Kaheman et al., 2022] or methods which learn from
statistics derived from time-series, such as moments or correlation functions
[Schneider et al., 2021]. Keller and Du [2021], Du et al. [2021] provide
rigorous analysis demonstrating how inference of a continuous-time model from
discrete-time data must be conducted with great care; they prove how stable
and consistent linear multistep methods for continuous-time integration may
not possess the same guarantees when used for the inverse problem, i.e.
discovery of dynamics. Queiruga et al. [2020] provide pathological
illustrations of this phenomenon in the context of Runge-Kutta methods.

Discrete-time approaches, on the other hand, are easily deployed when train
and test data sample rates are the same. For applications in which data
collection is easily configured (e.g. simulated settings, available automatic
sensors, etc.), discrete-time methods are typically much easier to implement
and test than continuous-time methods. Moreover, they allow for “non-
intrusive” model correction, as additions are applied outside of the numerical
integrator; this may be relevant for practical integration with complex
simulation software. In addition, discrete-time approaches can be preferable
when there is unavoidably large error in continuous-time inference Chorin and
Lu [2015], Lu et al. [2016].

Both non-parametric and parametric model classes are used in the learning of
dynamical systems, with the latter connecting to the former via the
representer theorem, when Gaussian process regression [Rasmussen and Williams,
2006] is used [Burov et al., 2020, Gilani et al., 2021, Harlim et al., 2021].

####  1.1.2. Hybrid Mechanistic and Data-Driven Modeling

Attempts to transform domains that have relied on traditional mechanistic
models, by using purely data-driven (i.e. _de novo_ or “learn from scratch”)
approaches, often fall short. Now, there is a growing recognition by machine
learning researchers that these mechanistic models are very valuable [Miller
et al., 2021], as they represent the distillation of centuries of data
collected in countless studies and interpreted by domain experts. Recent
studies have consistently found advantages of hybrid methods that blend
mechanistic knowledge and data-driven techniques; Willard et al. [2021]
provide a thorough review of this shift amongst scientists and engineers. Not
only do these hybrid methods improve predictive performance [Pathak et al.,
2018b], but they also reduce data demands [Rackauckas et al., 2020] and
improve interpretability and trustworthiness, which is essential for many
applications. This is exemplified by work in autonomous drone landing [Shi et
al., 2018] and helicopter flying [Rackauckas et al., 2021], as well as
predictions for COVID-19 mortality risk [Sottile et al., 2021] and COVID-19
treatment response [Qian et al., 2021].

The question of how best to use the power of data and machine learning to
leverage and build upon our existing mechanistic knowledge is thus of
widespread current interest. This question and research direction has been
anticipated over the last thirty years of research activity at the interface
of dynamical systems with machine learning [Rico-Martinez et al., 1994, Wilson
and Zorzetto, 1997, Lovelett et al., 2020], and now a critical mass of effort
is developing. A variety of studies have been tackling these questions in
weather and climate modeling [Kashinath et al., 2021, Farchi et al., 2021];
even in the imaging sciences, where pure machine learning has been
spectacularly successful, emerging work shows that incorporating knowledge of
underlying physical mechanisms improves performance in a variety of image
recognition tasks [Ba et al., 2019].

As noted and studied by Ba et al. [2019], Freno and Carlberg [2019] and
others, there are a few common high-level approaches for blending machine
learning with mechanistic knowledge: (1) use machine learning to learn
additive residual corrections for the mechanistic model [Saveriano et al.,
2017, Shi et al., 2018, Kaheman et al., 2019, Harlim et al., 2021, Farchi et
al., 2021, Lu et al., 2017, Lu, 2020, Yin et al., 2021]; (2) use the
mechanistic model as an input or feature for a machine learning model [Pathak
et al., 2018b, Lei et al., 2020, Borra et al., 2020]; (3) use mechanistic
knowledge in the form of a differential equation as a final layer in a neural
network representation of the solution, or equivalently define the loss
function to include approximate satisfaction of the differential equation
[Raissi et al., 2019, 2018, Chen et al., 2021a, Smith et al., 2020]; and (4)
use mechanistic intuition to constrain or inform the machine learning
architecture [Haber and Ruthotto, 2017, Maulik et al., 2021]. Many other
successful studies have developed specific designs that further hybridize
these and other perspectives [Hamilton et al., 2017, Freno and Carlberg, 2019,
Yi Wan et al., 2020, Jia et al., 2021]. In addition, parameter estimation for
mechanistic models is a well-studied topic in data assimilation, inverse
problems, and other mechanistic modeling communities, but recent approaches
that leverage machine learning for this task may create new opportunities for
accounting for temporal parameter variations [Miller et al., 2020] and unknown
observation functions [Linial et al., 2021].

An important distinction should be made between physics-informed _surrogate_
modeling and what we refer to as _hybrid_ modeling. Surrogate modeling
primarily focuses on replacing high-cost, high-fidelity mechanistic model
simulations with similarly accurate models that are cheap to evaluate. These
efforts have shown great promise by training machine learning models on
expensive high-fidelity simulation data, and have been especially successful
when the underlying physical (or other domain-specific) mechanistic knowledge
and equations are incorporated into the model training [Raissi et al., 2019]
and architecture [Maulik et al., 2021]. We use the term hybrid modeling, on
the other hand, to indicate when the final learned system involves interaction
(and possibly feedback) between mechanism-based and data-driven models [Pathak
et al., 2018b].

In this work, we focus primarily on hybrid methods that learn residuals to an
imperfect mechanistic model. We closely follow the discrete-time hybrid
modeling framework developed by Harlim et al. [2021], while providing new
insights from the continuous-time modeling perspective. The benefits of this
form of hybrid modeling, which we and many others have observed, are not yet
fully understood in a theoretical sense. Intuitively, nominal mechanistic
models are most useful when they encode key nonlinearities that are not
readily inferred using general model classes and modest amounts of data.
Indeed, classical approximation theorems for fitting polynomials, fourier
modes, and other common function bases directly reflect this relationship by
bounding the error with respect to a measure of complexity of the target
function (e.g. Lipschitz constants, moduli of continuity, Sobolev norms, etc.)
[DeVore and Lorentz, 1993][Chapter 7]. Recent work by E et al. [2019] provides
a priori error bounds for two-layer neural networks and kernel-based
regressions, with constants that depend explicitly on the norm of the target
function in the model-hypothesis space (a Barron space and a reproducing
kernel Hilbert space, resp.). At the same time, problems for which mechanistic
models only capture low-complexity trends (e.g. linear) may still be good
candidates for hybrid learning (over purely data-driven), as an accurate
linear model reduces the parametric burden for the machine-learning task; this
effect is likely accentuated in data-starved regimes. Furthermore, even in
cases where data-driven models perform satisfactorily, a hybrid approach may
improve interpretability, trustworthiness, and controllability without
sacrificing performance.

Hybrid models are often cast in Markovian, memory-free settings where the
learned dynamical system (or its learned residuals) are solely dependent on
the observed states. This approach can be highly effective when measurements
of all relevant states are available or when the influence of the unobserved
states is adequately described by a function of the observables. This is the
perspective employed by Shi et al. [2018], where they learn corrections to
physical equations of motion for an autonomous vehicle in regions of state
space where the physics perform poorly— these residual errors are _driven_ by
un-modeled turbulence during landing, but can be _predicted_ using the
observable states of the vehicle (i.e. position, velocity, and acceleration).
This is also the perspective taken in applications of high-dimensional
multiscale dynamical systems, wherein sub-grid closure models parameterize the
effects of expensive fine-scale interactions (e.g. cloud simulations) as
functions of the coarse variables [Grabowski, 2001, Khairoutdinov and Randall,
2001, Tan et al., 2018, Brenowitz and Bretherton, 2018, O’Gorman and Dwyer,
2018, Rasp et al., 2018, Schneider et al., 2021, Beucler et al., 2021]. The
result is a hybrid dynamical system with a physics-based equation defined on
the coarse variables with a Markovian correction term that accounts for the
effects of the expensive fine scale dynamics.

####  1.1.3. Non-Markovian Data-Driven Modeling

Unobserved and unmodeled processes are often responsible for model errors that
cannot be represented in a Markovian fashion within the observed variables
alone. This need has driven substantial advances in memory-based modeling. One
approach to this is the use of delay embeddings [Takens, 1981]. These methods
are inherently tied to discrete time representations of the data and, although
very successful in many applied contexts, are of less value when the goal of
data-driven learning is to fit continuous-time models; this is a desirable
modeling goal in many settings.

An alternative to understanding memory is via the Mori-Zwanzig formalism,
which is a fundamental building block in the presentation of memory and hidden
variables and may be employed for both discrete-time and continuous-time
models. Although initially developed primarily in the context of statistical
mechanics, it provides the basis for understanding hidden variables in
dynamical systems, and thus underpins many generic computational tools applied
in this setting [Chorin et al., 2000, Zhu et al., 2018, Gouasmi et al., 2017].
It has been successfully applied to problems in fluid turbulence [Duraisamy et
al., 2019, Parish and Duraisamy, 2017] and molecular dynamics [Li et al.,
2017, Hijón et al., 2010]. Lin and Lu [2021] demonstrate connections between
Mori-Zwanzig and delay embedding theory in the context of non-linear
autoregressive models using Koopman operator theory. Indeed, Gilani et al.
[2021] shows a correspondence between the Mori-Zwanzig representation of the
Koopman operator and Taken’s delay-embedding flow map. Studies by Ma et al.
[2019], Wang et al. [2020] demonstrate how the Mori-Zwanzig formalism
motivates the use of recurrent neural networks (RNNs) [Rumelhart et al., 1986,
Goodfellow et al., 2016] as a deep learning approach to non-Markovian closure
modeling. Harlim et al. [2021] also use the Mori-Zwanzig formalism to deduce a
non-Markovian closure model, and evaluate RNN-based approximations of the
closure dynamics. Closure modeling using RNNs has recently emerged as a new
way to learn memory-based closures [Kani and Elsheikh, 2017, Chattopadhyay et
al., 2020b, Harlim et al., 2021].

Although the original formulation of Mori-Zwanzig as a general purpose
approach to modeling partially observed systems was in continuous-time [Chorin
et al., 2000], many practical implementations adopt a discrete-time picture
[Darve et al., 2009, Chorin and Lu, 2015, Lin and Lu, 2021]. This causes the
learned memory terms to depend on sampling rates, which, in turn, can inhibit
flexibility and interpretability.

Recent advances in continuous-time memory-based modeling, however, may be
applicable to these non-Markovian hybrid model settings. The theory of
continuous-time RNNs (i.e. formulated as differential equations, rather than a
recurrence relation) was studied in the 1990s [Funahashi and Nakamura, 1993,
Beer, 1995], albeit for equations with a specific additive structure. This
structure was exploited in a continuous-time reservoir computing (RC) approach
by Lu et al. [2018] for reconstructing chaotic attractors from data.
Comparisons between RNNs and RC (a subclass of RNNs with random parameters
fixed in the recurrent state) in discrete-time have yielded mixed conclusions
in terms of their relative efficiencies and ability to retain memory [Pyle et
al., 2021, Gauthier et al., 2021, Vlachas et al., 2020, Chattopadhyay et al.,
2020a]. Recent formulations of continuous-time RNNs have departed slightly
from the additive structure, and have focused on constraints and architectures
that ensure stability and accuracy of the resulting dynamical system [Chang et
al., 2019, Erichson et al., 2020, Niu et al., 2019, Rubanova et al., 2019,
Sherstinsky, 2020, Ouala et al., 2020]. In addition, significant theoretical
work has been performed for linear RNNs in continuous-time [Li et al., 2020].
Nevertheless, these various methods have not yet been formulated within a
hybrid modeling framework, nor has their approximation power been carefully
evaluated in that context. A recent step in this direction, however, is the
work by Gupta and Lermusiaux [2021], which tackles non-Markovian hybrid
modeling in continuous-time with neural network-based delay differential
equations (DDEs).

####  1.1.4. Noisy Observations and Data Assimilation

For this work we consider settings in which the observations may be both noisy
and partial; the observations may be partial either because the system is
undersampled in time, or because certain variables are not observed at all. We
emphasize that ideas from statistics can be used to smooth and/or interpolate
data to remove noise and deal with undersampling Craven and Wahba [1978] and
to deal with missing data Meng and Van Dyk [1997]; and ideas from data
assimilation Asch et al. [2016], Law et al. [2015], Reich and Cotter [2015]
can be used to remove noise and to learn about unobserved variables Chen et
al. [2021b], Gottwald and Reich [2021b], Gottwald and Reich [2021a]. In some
of our experiments we will use noise-free data in continuous-time, to clearly
expose issues separate from noise/interpolation; but in other experiments we
will use methodologies from data assimilation to enhance our learning Chen et
al. [2021b].

####  1.1.5. Applications of Data-Driven Modeling

In order to deploy hybrid methods in real-world scenarios, we must also be
able to cope with noisy, partial observations. Accommodating the learning of
model error in this setting, as well as state estimation, is an active area of
research in the data assimilation (DA) community [Pulido et al., 2018, Farchi
et al., 2021, Bocquet et al., 2020]. Learning dynamics from noisy data is
generally non-trivial for nonlinear systems—there is a chicken-and-egg problem
in which accurate state estimation typically relies on the availability of
correct models, and correct models are most readily identified using accurate
state estimates. Recent studies have addressed this challenge by attempting to
jointly learn the noise and the dynamics. Gottwald and Reich [2021b] approach
this problem from a data assimilation perspective, and employ an Ensemble
Kalman Filter (EnKF) to iteratively update the parameters for their dynamics
model, then filter the current state using the updated dynamics. A recent
follow-up to this work applies the DA-approach to partially-observed systems,
and learns a model on a space of discrete-time delay-embeddings [Gottwald and
Reich, 2021a]. Similar studies were performed by Brajard et al. [2021], and
applied specifically in model error scenarios [Brajard et al., 2020, Farchi et
al., 2021, Wikner et al., 2020]. Ayed et al. [2019] focus on learning a
continuous-time neural network representation of an ODE from partial
observations, and learn a separate encoder neural network to map a historical
warmup sequence to likely initial conditions in the un-observed space. Kaheman
et al. [2019] approach this problem from a variational perspective, performing
a single optimization over all noise sequences and dynamical
parameterizations. Nguyen et al. [2019] use an Expectation-Maximization (EM)
perspective to compare these variational and ensemble-based approaches, and
further study is needed to understand the trade-offs between these styles of
optimization. Chen et al. [2021b] study an EnKF-based optimization scheme that
performs joint, rather than EM-based learning, by running gradient descent on
an architecture that backpropagates through the data assimilator.

We note that data assimilators are themselves dynamical systems, which can be
tuned (using optimization and machine learning) to provide more accurate state
updates and more efficient state identification. However, while learning
improved DA schemes is sometimes viewed as a strategy for coping with model
error [Zhu and Kamachi, 2000], we see the optimization of DA and the
correction of model errors as two separate problems that should be addressed
individually.

When connecting models of dynamical systems to real-world data, it is also
essential to recognize that available observables may live in a very different
space than the underlying dynamics. Recent studies have shown ways to navigate
this using autoencoders and dynamical systems models to jointly learn a latent
embedding and dynamics in that latent space [Champion et al., 2019]. Proof of
concepts for similar approaches primarily focus on image-based inputs, but
have potential for applications in medicine [Linial et al., 2021] and
reduction of nonlinear PDEs [Maulik et al., 2021].

###  1.2. Our Contributions

Despite this large and recent body of work in data-driven learning methods and
hybrid modeling strategies, many challenges remain for understanding how to
best combine mechanistic and machine-learned models; indeed, the answer is
highly dependent on the application. Here, we construct a mathematical
framework that unifies many of the common approaches for blending mechanistic
and machine learning models; having done so we provide strong evidence for the
value of hybrid approaches. Our contributions are listed as follows:

  1. (1)

We provide an overarching framework for learning model error from (possibly
noisy) data in dynamical systems settings, studying both discrete- and
continuous-time models, together with both memoryless (Markovian) and memory-
dependent representations of the model error. This formulation is agnostic to
choice of mechanistic model and class of machine learning functions.

  2. (2)

We study the Markovian learning problem in the context of ergodic continuous-
time dynamics, proving bounds on excess risk and generalization error.

  3. (3)

We present a simple approximation theoretic approach to learning memory-
dependent (non-Markovian) model error in continuous-time, proving a form of
universal approximation for two families of memory-dependent model error
defined using recurrent neural networks.

  4. (4)

We describe numerical experiments which: a) demonstrate the utility of
learning model error in comparison both with pure data-driven learning and
with pure (but slightly imperfect) mechanistic modeling; b) compare the
benefits of learning discrete- versus continuous-time models; c) demonstrate
the utility of auto-differentiable data assimilation to learn dynamics from
partially observed, noisy data; d) explain issues arising in memory-dependent
model error learning in the (typical) situation where the dimension of the
memory variable is unknown.

In Section 2, we address contribution 1. by defining the general settings of
interest for dynamical systems in both continuous- and discrete-time. We then
link these underlying systems to a machine learning framework in Sections 3
and 4; in the former we formulate the problem in the setting of statistical
learning, and in the latter we define concrete optimization problems found
from finite parameterizations of the hypothesis class in which the model error
is sought. Section 5 is focused on specific choices of architectures, and
underpinning theory for machine learning methods with these choices: we
analyze linear methods from the perspective of learning theory in the context
of ergodic dynamical systems (contribution 2.); and we describe an
approximation theorem for continuous-time hybrid recurrent neural networks
(contribution 3.). Finally, Section 6 presents our detailed numerical
experiments; we apply the methods in Section 5 to exemplar dynamical systems
of the forms outlined in Section 2, and highlight our findings (contribution
4.).

##  2\. Dynamical Systems Setting

In the following, we use the phrase _Markovian model error_ to describe model
error expressible entirely in terms of the observed variable at the current
time, the memoryless situation; _non-Markovian model error_ refers to the need
to express the model error in terms of the past history of the observed
variable.

We present a general framework for modeling a dynamical system with Markovian
model error, first in continuous-time (Section 2.1) and then in discrete-time
(Section 2.2). We then extend the framework to the setting of non-Markovian
model error (Section 2.3), including a parameter $\varepsilon$ which enables
us to smoothly transition from scale-separated problems (where Markovian
closure is likely to be accurate) to problems where the unobserved variables
are not scale-separated from those observed (where Markovian closure is likely
to fail and memory needs to be accounted for).

It is important to note that the continuous-time formulation necessarily
assumes an underlying data-generating process that is continuous in nature.
The discrete-time formulation can be viewed as a discretization of an
underlying continuous system, but can also represent systems that are truly
discrete.

The settings that we present are all intended to represent and classify common
situations that arise in modeling and predicting dynamical systems. In
particular, we stress two key features. First, we point out that mechanistic
models (later referred to as a vector field $f_{0}$ or flow map $\Psi_{0}$)
are often available and may provide predictions with reasonable fidelity.
However, these models are often simplifications of the true system, and thus
can be improved with data-driven approaches. Nevertheless, they provide a
useful starting point that can reduce the complexity and data-hunger of the
learning problems. In this context, we study trade-offs between discrete- and
continuous-time framings. While we begin with fully-observed contexts in which
the dynamics are Markovian with respect to the observed state $x$, we later
note that we may only have access to partial observations $x$ of a larger
system $(x,y)$. By restricting our interest to prediction of these
observables, we show how a latent dynamical process (e.g. a RNN) has the power
to reconstruct the correct dynamics for our observables.

###  2.1. Continuous-Time

Consider the following dynamical system

(2.1) $\displaystyle\dot{x}=f^{\dagger}(x),\quad x(0)=x_{0},$

and define $\mathsf{X}_{s}:=C([0,s];\mathbbm{R}^{d_{x}})$. If $f^{\dagger}\in
C^{1}(\mathbbm{R}^{d_{x}};\mathbbm{R}^{d_{x}})$ then (2.1) has solution
$x(\cdot)\in\mathsf{X}_{T}$ for any
$T<T_{\rm{max}}=T_{\rm{max}}(x_{0})\in\mathbbm{R}^{+}$, the maximal interval
of existence.

The primary model error scenario we envisage in this section is one in which
the vector field $f^{\dagger}$ can only be partially known or accessed: we
assume that

$f^{\dagger}=f_{0}+m^{\dagger}$

where $f_{0}$ is known to us and $m^{\dagger}$ is not known. For any $f_{0}\in
C^{1}(\mathbbm{R}^{d_{x}};\mathbbm{R}^{d_{x}})$ (regardless of its fidelity),
there exists a function $m^{\dagger}(x)\in
C^{1}(\mathbbm{R}^{d_{x}};\mathbbm{R}^{d_{x}})$ such that (2.1) can be
rewritten as

(2.2) $\dot{x}=f_{0}(x)+m^{\dagger}(x).$

However, for this paper, it is useful to think of $m^{\dagger}$ as being small
relative to $f_{0}$; the function $m^{\dagger}$ accounts for _model error._
While the approach in (2.2) is targeted at learning residuals of $f_{0}$,
$f^{\dagger}$ can alternatively be reconstructed from $f_{0}$ through a
different function $m^{\dagger}(x)\in
C^{1}(\mathbbm{R}^{2d_{x}};\mathbbm{R}^{d_{x}})$ using the form

(2.3) $\dot{x}=m^{\dagger}(x,f_{0}(x)).$

Both approaches are defined on spaces that allow perfect reconstruction of
$f^{\dagger}$. However, the first formulation hypothesizes that the missing
information is additive; the second formulation provides no such indication.
Because the first approach ensures substantial usage of $f_{0}$, it has
advantages in settings where $f_{0}$ is trusted by practitioners and model
explainability is important. The second approach will likely see advantages in
settings where there is a simple non-additive form of model error, including
coordinate transformations and other (possibly state-dependent) nonlinear
warping functions of the nominal physics $f_{0}$. Note that the use of $f_{0}$
in representing the model error in the augmented-input setting of (2.3)
includes the case of not leveraging $f_{0}$ at all. It is, hence, potentially
more useful than simply adopting an $x-$dependent model error; but it requires
learning a more complex function.

The augmented-input method also has connections to model stacking [Wolpert,
1992] or bagging [Breiman, 1996]; this perspective can be useful when there
are $N$ model hypotheses:

$\dot{x}=m^{\dagger}\big{(}x,f^{(1)}_{0}(x),\
\ldots\,f^{(N)}_{0}(x);\theta\big{)}.$

The residual-based design in (2.2) relates more to model boosting [Schapire,
1990].

Our goal is to use machine learning to approximate these corrector functions
$m^{\dagger}$ using our nominal knowledge $f_{0}$ and observations of a
trajectory $\\{x(t)\\}_{t=0}^{T}\in\mathsf{X}_{T}$, for some
$T<T_{\rm{max}}(x_{0})$, from the true system (2.1). In this work, we consider
only the case of learning $m^{\dagger}(x)$ in equation (2.2). For now the
reader may consider $\\{x(t)\\}_{t=0}^{T}$ given without noise so that, in
principle, $\\{\dot{x}(t)\\}_{t=0}^{T}$ is known and may be leveraged. In
practice this will not be the case, for example if the data are high-frequency
but discrete in time; we address this issue in what follows.

###  2.2. Discrete-Time

Consider the following dynamical system

(2.4) $x_{k+1}=\Psi^{\dagger}(x_{k})$

and define
$\mathsf{X}_{K}:=\ell^{\infty}\bigl{(}\\{0,\dots,K\\};\mathbbm{R}^{d_{x}}\bigl{)}$.
If $\Psi^{\dagger}\in C(\mathbbm{R}^{d_{x}};\mathbbm{R}^{d_{x}})$, the map
yields solution
$\\{x_{k}\\}_{k\in\mathbb{Z}^{+}}\in\mathsf{X}_{\infty}:=\ell^{\infty}\bigl{(}\mathbb{Z}^{+};\mathbbm{R}^{d_{x}}\bigl{)}.$111Here
we define $\mathbb{Z}^{+}=\\{0,\dots,\\}$, the non-negative integers,
including zero. As in the continuous-time setting, we assume we only have
access to an approximate mechanistic model $\Psi_{0}\in
C(\mathbbm{R}^{d_{x}};\mathbbm{R}^{d_{x}})$, which can be corrected using an
additive residual term $m^{\dagger}\in
C(\mathbbm{R}^{d_{x}};\mathbbm{R}^{d_{x}})$:

(2.5) $x_{k+1}=\Psi_{0}(x_{k})+m^{\dagger}(x_{k}),$

or by feeding $\Psi_{0}$ as an input to a corrective warping function
$m^{\dagger}\in C(\mathbbm{R}^{2d_{x}};\mathbbm{R}^{d_{x}})$

$x_{k+1}=m^{\dagger}(x_{k},\Psi_{0}(x_{k}));$

we focus our experiments on the additive residual framing in (2.5).

Note that the discrete-time formulation can be made compatible with
continuous-time data sampled uniformly at rate $\Delta t$ (i.e. $x(k\Delta
t)=x_{k}$ for $k\in\mathbbm{N}$). To see this, let
$\Phi^{\dagger}(x_{0},t):=x(t)$ be the solution operator for (2.1) (and
$\Phi_{0}$ defined analogously for $f_{0}$). We then have

(2.6a) $\displaystyle\Psi^{\dagger}(v)$
$\displaystyle:=\Phi^{\dagger}(v,\Delta t)$ (2.6b) $\displaystyle\Psi_{0}(v)$
$\displaystyle:=\Phi_{0}(v,\Delta t),$

which can be obtained via numerical integration of $f^{\dagger},f_{0}$,
respectively.

###  2.3. Partially Observed Systems (Continuous-Time)

The framework in Sections 2.1 and 2.2 assumes that the system dynamics are
Markovian with respect to observable $x$. Most of our experiments are
performed in the fully-observed Markovian case. However, this assumption
rarely holds in real-world systems. Consider a block-on-a-spring experiment
conducted in an introductory physics laboratory. In principle, the system is
strictly governed by the position and momentum of the block (i.e. $f_{0}$),
along with a few scalar parameters. However (as most students’ error analysis
reports will note), the dynamics are also driven by a variety of external
factors, like a wobbly table or a poorly greased track. The magnitude,
timescale, and structure of the influence of these different factors are
rarely known; and yet, they are somehow encoded in the discrepancy between the
nominal equations of motion and the (noisy) observations of this multiscale
system.

Thus we also consider the setting in which the dynamics of $x$ is not
Markovian. If we consider $x$ to be the observable states of a Markovian
system in dimension higher than $d_{x}$, then we can write the full system as

(2.7a) $\displaystyle\dot{x}$ $\displaystyle=f^{\dagger}(x,y),\quad
x(0)=x_{0}$ (2.7b) $\displaystyle\dot{y}$
$\displaystyle=\frac{1}{\varepsilon}g^{\dagger}(x,y),\quad y(0)=y_{0}.$

Here $f^{\dagger}\in
C^{1}(\mathbbm{R}^{d_{x}}\times\mathbbm{R}^{d_{y}};\mathbbm{R}^{d_{x}})$,
$g^{\dagger}\in
C^{1}(\mathbbm{R}^{d_{x}}\times\mathbbm{R}^{d_{y}};\mathbbm{R}^{d_{y}})$, and
$\varepsilon>0$ is a constant measuring the degree of scale-separation (which
is large when $\varepsilon$ is small). The system yields solution 222With
$\mathsf{Y}_{T}$ defined analogously to $\mathsf{X}_{T}.$
$x(\cdot)\in\mathsf{X}_{T},y(\cdot)\in\mathsf{Y}_{T}$ for any
$T<T_{\rm{max}}(x(0),y(0))\in\mathbb{R}^{+},$ the maximal interval of
existence. We view $y$ as the complicated, unresolved, or unobserved aspects
of the true underlying system.

For any $f_{0}\in C^{1}(\mathbbm{R}^{d_{x}};\mathbbm{R}^{d_{x}})$ (regardless
of its fidelity), there exists a function $m^{\dagger}(x,y)\in
C^{1}(\mathbbm{R}^{d_{x}}\times\mathbbm{R}^{d_{y}};\mathbbm{R}^{d_{x}})$ such
that (2.7) can be rewritten as

(2.8a) $\displaystyle\dot{x}$ $\displaystyle=f_{0}(x)+m^{\dagger}(x,y)$ (2.8b)
$\displaystyle\dot{y}$ $\displaystyle=\frac{1}{\varepsilon}g^{\dagger}(x,y).$

Now observe that, by considering the solution of equation (2.8b) as a function
of the history of $x$, the influence of $y(\cdot)\in\mathsf{Y}_{t}$ on the
solution $x(\cdot)\in\mathsf{X}_{t}$ can be captured by a parameterized
(w.r.t. $t$) family of operators
$m^{\dagger}_{t}\colon\mathsf{X}_{t}\times\mathbbm{R}^{d_{y}}\times\mathbbm{R}^{+}\mapsto\mathbbm{R}^{d_{x}}$
on the historical trajectory $\\{x(s)\\}_{s=0}^{t}$, unobserved initial
condition $y(0)$, and scale-separation parameter $\varepsilon$ such that

(2.9) $\displaystyle\dot{x}(t)$
$\displaystyle=f_{0}\bigl{(}x(t)\bigr{)}+m^{\dagger}_{t}\bigl{(}\\{x(s)\\}_{s=0}^{t};\
y(0),\ \varepsilon\bigr{)}.$

Our goal is to use machine learning to find a Markovian model, in which $x$ is
part of the state variable, using our nominal knowledge $f_{0}$ and
observations of a trajectory $\\{x(t)\\}_{t=0}^{T}\in\mathsf{X}_{T}$, for some
$T<T_{\rm{max}}(x_{0},y_{0})$, from the true system (2.7); note that
$y(\cdot)$ is not observed and nothing is assumed known about the vector field
$g^{\dagger}$ or the parameter $\varepsilon.$

Note that equations (2.7), (2.8) and (2.9) are all equivalent formulations of
the same problem and have identical solutions. The third formulation points
towards two intrinsic difficulties: the unknown “function” to be learned is in
fact defined by a family of operators $m^{\dagger}_{t}$ mapping the Banach
space of path history into $\mathbbm{R}^{d_{x}}$; secondly the operator is
parameterized by $y(0)$ which is unobserved. We will address the first issue
by showing that the operators $m^{\dagger}_{t}$ can be arbitrarily well-
approximated from within a family of differential equations in dimension
$\mathbbm{R}^{2d_{x}+d_{y}}$; the second issue may be addressed by techniques
from data assimilation [Asch et al., 2016, Law et al., 2015, Reich and Cotter,
2015] once this approximating family is learned. We emphasize, however, that
we do not investigate the practicality of this approach to learning non-
Markovian systems and much remains to be done in this area.

It is also important to note that these non-Markovian operators
$m^{\dagger}_{t}$ can sometimes be adequately approximated by invoking a
Markovian model for $x$ and simply learning function $m^{\dagger}(\cdot)$ as
in Section 2.1. For example, when $\varepsilon\to 0$ and the $y$ dynamics,
with $x$ fixed, are sufficiently mixing, the averaging principle [Bensoussan
et al., 2011, Vanden-Eijnden and others, 2003, Pavliotis and Stuart, 2008] may
be invoked to deduce that

$\lim_{\varepsilon\to 0}m^{\dagger}_{t}\bigl{(}\\{x(s)\\}_{s=0}^{t};\
y(0),\varepsilon\bigr{)}=m^{\dagger}(x(t))$

for some $m^{\dagger}$ as in Section 2.1. This fact is used in section 3 of
[Jiang and Harlim, 2020] to study the learning of closure models for linear
Gaussian stochastic differential equations (SDEs).

It is highly advantageous to identify settings where Markovian modeling is
sufficient, as it is a simpler learning problem. We find that learning
$m^{\dagger}_{t}$ is necessary when there is significant memory required to
explain the dynamics of $x$; learning $m^{\dagger}$ is sufficient when memory
effects are minimal. In Section 6, we show that Markovian closures can perform
well for certain tasks even when the scale-separation factor $\varepsilon$ is
not small. In Section 3 we demonstrate how the family of operators
$m^{\dagger}_{t}$ may be represented through ODEs, appealing to ideas which
blend continuous-time RNNs with an assumed known vector field $f_{0}.$

###  2.4. Partially Observed Systems (Discrete-Time)

The discrete-time analog of the previous setting considers a mapping

(2.10a) $\displaystyle x_{k+1}$
$\displaystyle=\Psi_{1}^{\dagger}(x_{k},y_{k})$ (2.10b) $\displaystyle
y_{k+1}$ $\displaystyle=\Psi_{2}^{\dagger}(x_{k},y_{k})$

with $\Psi_{1}^{\dagger}\in
C(\mathbbm{R}^{d_{x}}\times\mathbbm{R}^{d_{y}};\mathbbm{R}^{d_{x}})$,
$\Psi_{2}^{\dagger}\in
C(\mathbbm{R}^{d_{x}}\times\mathbbm{R}^{d_{y}};\mathbbm{R}^{d_{y}})$, yielding
solutions $\\{x_{k}\\}_{k\in\mathbb{Z}^{+}}\in\mathsf{X}_{\infty}$ and
$\\{y_{k}\\}_{k\in\mathbb{Z}^{+}}\in\mathsf{Y}_{\infty}.$ We assume unknown
$\Psi_{1}^{\dagger},\Psi_{2}^{\dagger}$, but known approximate model
$\Psi_{0}$ to rewrite (2.10) as

(2.11a) $\displaystyle x_{k+1}$
$\displaystyle=\Psi_{0}(x_{k})+m^{\dagger}(x_{k},y_{k})$ (2.11b)
$\displaystyle y_{k+1}$ $\displaystyle=\Psi_{2}^{\dagger}(x_{k},y_{k}).$

We can, analogously to (2.9), write a solution in the space of observables as

(2.12) $\displaystyle x_{k+1}$
$\displaystyle=\Psi_{0}\bigl{(}x_{k}\bigr{)}+m^{\dagger}_{k}\bigl{(}\\{x_{s}\\}_{s=0}^{k},\
y_{0}\bigr{)}$

with
$m^{\dagger}_{k}\colon\mathsf{X}_{k}\times\mathbbm{R}^{d_{y}}\to\mathbbm{R}^{d_{x}}$,
a function of the historical trajectory $\\{x_{s}\\}_{s=0}^{k}$ and the
unobserved initial condition $y_{0}$. If this discrete-time system is computed
from the time $\Delta t$ map for (2.1) then, for $\varepsilon\ll 1$ and when
averaging scenarios apply as discussed in Section 2.3, the memoryless model in
(2.5) may be used.

##  3\. Statistical Learning for Ergodic Dynamical Systems

Here, we present a learning theory framework within which to consider methods
for discovering model error from data. We outline the learning theory in a
continuous-time Markovian setting (using possibly discretely sampled data),
and point to its analogs in discrete-time and non-Markovian settings.

In the discrete-time settings, we assume access to discretely sampled training
data $\\{x_{k}=x(k\Delta t)\\}_{k=0}^{K}$, where $\Delta t$ is a uniform
sampling rate and we assume that $K\Delta t=T.$ In the continuous-time
settings, we assume access to continuous-time training data
$\\{\dot{x}(t),x(t)\\}_{t=0}^{T}$; Section 6.2.1 discusses the important
practical question of estimating $\dot{x}(t),x(t)$ from discrete (but high
frequency) data. In either case, consider the problem of identifying
$m\in\mathcal{M}$ (where $\mathcal{M}$ represents the model, or hypothesis,
class) that minimizes a loss function quantifying closeness of $m$ to
$m^{\dagger}.$ In the Markovian setting we choose a measure $\mu$ on
$\mathbbm{R}^{d_{x}}$ and define the loss

$\mathcal{L}_{\mu}(m,m^{\dagger}):=\int_{\mathbbm{R}^{d_{x}}}\|m(x)-m^{\dagger}(x)\|^{2}_{2}d\mu(x).$

If we assume that, at the true $m^{\dagger}$, $x(\cdot)$ is ergodic with
invariant density $\mu$, then we can exchange time and space averages to see,
for infinitely long trajectory $\\{x(t)\\}_{t\geq 0}$,

$\displaystyle\mathcal{I}_{\infty}(m)$
$\displaystyle:=\lim_{T\to\infty}\frac{1}{T}\int_{0}^{T}\|m(x(t))-m^{\dagger}(x(t))\|^{2}_{2}dt$
$\displaystyle=\int_{\mathbbm{R}^{d_{x}}}\|m(x)-m^{\dagger}(x)\|^{2}_{2}d\mu(x)$
$\displaystyle=\mathcal{L}_{\mu}(m,m^{\dagger}).$

Since we may only have access to a trajectory dataset of finite length $T$, it
is natural to define

$\mathcal{I}_{T}(m):=\frac{1}{T}\int_{0}^{T}\|m(x(t))-m^{\dagger}(x(t))\|^{2}_{2}dt$

and note that, by ergodicity,

$\lim_{T\to\infty}\mathcal{I}_{T}(m)=\mathcal{L}_{\mu}(m,m^{\dagger}).$

Finally, we can use (2.2) to get

(3.1)
$\mathcal{I}_{T}(m)=\frac{1}{T}\int_{0}^{T}\|\dot{x}(t)-f_{0}(x(t))-m(x(t))\|^{2}_{2}dt.$

This, possibly regularized, is a natural loss function to employ when
continuous-time data is available, and should be viewed as approximating
$\mathcal{L}_{\mu}(m,m^{\dagger}).$ We can use these definitions to frame the
problem of learning model error in the language of statistical learning
[Vapnik, 2013].

If we let $\mathcal{M}$ denote the hypothesis class over which we seek to
minimize $\mathcal{I}_{T}(m)$ then we may define

$m^{\ast}_{\infty}=\operatorname*{argmin}_{m\in\mathcal{M}}\mathcal{L}_{\mu}(m,m^{\dagger})=\operatorname*{argmin}_{m\in\mathcal{M}}\mathcal{I}_{\infty}(m),\quad
m^{\ast}_{T}=\operatorname*{argmin}_{m\in\mathcal{M}}\mathcal{I}_{T}(m).$

The risk associated with seeking to approximate $m^{\dagger}$ from the class
$\mathcal{M}$ is defined by
$\mathcal{L}_{\mu}(m^{\ast}_{\infty},m^{\dagger})$, noting that this is $0$ if
$m^{\dagger}\in\mathcal{M}.$ The risk measures the intrinsic error incurred by
seeking to learn $m^{\dagger}$ from the restricted class $\mathcal{M}$, which
typically does not include $m^{\dagger};$ it is an approximation theoretic
concept which encodes the richness of the hypothesis class $\mathcal{M}$. The
risk may be decreased by increasing the expressiveness of $\mathcal{M}$. Thus
risk is independent of the data employed. _Empirical risk minimization_ refers
to minimizing $\mathcal{I}_{T}$ (or a regularized version) rather than
$\mathcal{I}_{\infty},$ and this involves the specific instance of data that
is available. To quantify the effect of data volume on learning $m^{\dagger}$
through empirical risk minimization, it is helpful to introduce the following
two concepts. The excess risk is defined by

(3.2)
$R_{T}:=\mathcal{I}_{\infty}(m^{\ast}_{T})-\mathcal{I}_{\infty}(m^{\ast}_{\infty})$

and represents the additional approximation error incurred by using data
defined over a finite time horizon $T$ in the estimate of $m^{\dagger}.$ The
generalization error is

(3.3)
$G_{T}:=\mathcal{I}_{T}(m^{\ast}_{T})-\mathcal{I}_{\infty}(m^{\ast}_{T})$

and represents the discrepancy between training error, which is defined using
a finite trajectory, and idealized test error, which is defined using an
infinite length trajectory (or, equivalently, the invariant measure $\mu$),
when evaluated at the estimate of the function $m^{\dagger}$ obtained from
finite data. We return to study excess risk and generalization error in the
context of linear (in terms of parametric-dependence) models for
$m^{{\dagger}}$, and under ergodicity assumptions on the data generating
process, in Section 5.2.

We have introduced a machine learning framework in the continuous-time
Markovian setting, but it may be adopted in discrete-time and in non-Markovian
settings. In Section 4, we define appropriate objective functions for each of
these cases.

######  Remark 3.1.

The developments we describe here for learning in ODEs can be extended to the
case of learning SDEs; see [Bento et al., 2011, Kutoyants, 2004]. In that
setting, consistency in the large $T$ limit is well-understood. It would be
interesting to build on the learning theory perspective described here to
study statistical consistency for ODEs; the approaches developed in the work
by McGoff et al. [2015], Su and Mukherjee [2021] are potentially useful in
this regard. $\square$

##  4\. Parameterization of the Loss Function

In this section, we define explicit optimizations for learning (approximate)
model error functions $m^{\dagger}$ for the Markovian settings, and model
error operators $m^{\dagger}_{t}$ for the non-Markovian settings; both
continuous- and discrete-time formulations are given. We defer discussion of
specific approximation architectures to the next section. Here we make a
notational transition from optimization over (possibly non-parametric)
functions $m\in\mathcal{M}$ to functions parameterized by $\theta$ that
characterize the class $\mathcal{M}$.

In all the numerical experiments in this paper, we study the use of
continuous- and discrete-time approaches to model data generated by a
continuous-time process. The set-up in this section reflects this setting, in
which two key parameters appear: $T$, the continuous-time horizon for the
data; and $\Delta t$, the frequency of the data. The latter parameter will
always appear in the discrete-time models; but it may also be implicit in
continuous-time models which need to infer continuous-time quantities from
discretely sampled data. We relate $T$ and $\Delta t$ by $K\Delta t=T.$ We
present the general forms of $\mathcal{J}_{T}(\theta)$ (with optional
regularization terms $R(\theta)$). Optimization via derivative-based
methodology requires either analytic differentiation of the dynamical system
model with respect to parameters, or the use of autodifferentiable ODE solvers
[Rubanova et al., 2019].

###  4.1. Continuous-Time Markovian Learning

Here, we approximate the Markovian closure term in (2.2) with a parameterized
function $m(x;\ \theta)$. Assuming full knowledge of $\dot{x}(t),x(t)$, we
learn the correction term for the flow field by minimizing the following
objective function of $\theta$:

(4.1) $\displaystyle\mathcal{J}_{T}(\theta)$
$\displaystyle=\frac{1}{T}\int_{0}^{T}\big{\|}\dot{x}(t)-f_{0}(x(t))-m(x(t);\
\theta)\big{\|}^{2}dt+R(\theta)$

Note that $\mathcal{J}_{T}(\theta)=\mathcal{I}_{T}\bigl{(}m(\ \cdot\ ;\
\theta)\bigr{)}+R(\theta);$ thus the proposed methodology is a regularization
of the empirical risk minimization described in the preceding section.

Notable examples that leverage this framing include: the paper [Kaheman et
al., 2019], where $\theta$ are coefficients for a library of low-order
polynomials and $R(\theta)$ is a sparsity-promoting regularization defined by
the SINDy framework; the paper [Yin et al., 2021], where $\theta$ are
parameters of a deep neural network (DNN) and $L_{2}$ regularization is
applied to the weights; the paper [Shi et al., 2018], where $\theta$ are DNN
parameters and $R(\theta)$ encodes constraints on the Lipschitz constant for
$m$ provided by spectral normalization; and the paper [Watson, 2019] which
applies this approach to the Lorenz ’96 Multiscale system using neural
networks with an $L_{2}$ regularization on the weights.

###  4.2. Discrete-Time Markovian Learning

Here, we learn the Markovian correction term in (2.5) by minimizing:

(4.2) $\displaystyle\mathcal{J}_{T}(\theta)$
$\displaystyle=\frac{1}{K}\sum_{k=0}^{K-1}\big{\|}x_{k+1}-\Psi_{0}(x_{k})-m(x_{k};\
\theta)\big{\|}^{2}+R(\theta)$

This is the natural discrete-time analog of (4.1) and may be derived
analogously, starting from a discrete analog of the loss
$\mathcal{L}_{\mu}(m,m^{\dagger})$ where now $\mu$ is assumed to be an ergodic
measure for (2.4). If a discrete analog of (3.1) is defined, then
parameterization of $m$, and regularization, leads to (4.2). This is the
underlying model assumption in the work by Farchi et al. [2021].

###  4.3. Continuous-Time Non-Markovian Learning

We can attempt to recreate the dynamics in $x$ for (2.9) by modeling the non-
Markovian residual term. A common approach is to augment the dynamics space
with a variable $r\in\mathbbm{R}^{d_{r}}$ leading to a model of the form

(4.3a) $\displaystyle\dot{x}$ $\displaystyle=f_{0}(x)+f_{1}(r,x;\ \theta)$
(4.3b) $\displaystyle\dot{r}$ $\displaystyle=f_{2}(r,x;\ \theta).$

We then seek a $d_{r}$ large enough, and then parametric models
$\\{f_{j}(r,x;\cdot)\\}_{j=1}^{2}$ expressive enough, to ensure that the
dynamics in $x$ are reproduced by (4.3). Note that, although the model error
in $x$ is non-Markovian, as it depends on the history of $x$, we are seeking
to explain observed $x$ data by an enlarged model, including hidden variables
$r$, in which the dynamics of $[x,r]$ is Markovian.

When learning hidden dynamics from partial observations, we must jointly infer
the missing states $r(t)$ and the, typically parameterized, governing dynamics
$f_{1},f_{2}$. Furthermore, when the family of parametric models is not closed
with respect to translation of $r$ it will also be desirable to learn $r_{0}$;
when $x$ is observed noisily, it is similarly important to learn $x_{0}$.

To clarify discussions of (4.3) and its training from data, let $u=[x,r]$ and
$f$ be the concatenation of the vector fields given by $f_{0},f_{1},f_{2}$
such that

(4.4) $\dot{u}=f(u;\ \theta),$

with solution $u(t;\ v,\theta)$ solving (4.4) (and, equivalently, (4.3)) with
initial condition $v$ (i.e. $u(0;\ v,\theta)=v$). Consider observation
operators $H_{x},H_{r}$, such that $x=H_{x}u$, and $r=H_{r}u$, and further
define noisy observations of $x$ as

$z(t)=x(t)+\eta(t),$

where $\eta$ is i.i.d. observational noise. We now outline three optimization
approaches to learning from noisily, partially observed data $z$.

####  4.3.1. Optimization; Hard Constraint For Missing Dynamics

Since (4.3) is deterministic, it may suffice to jointly learn parameters and
initial condition $u(0)=u_{0}$ by minimizing [Rubanova et al., 2019]:

(4.5) $\displaystyle\mathcal{J}_{T}(\theta,u_{0})$
$\displaystyle=\frac{1}{T}\int_{0}^{T}\big{\|}z(t)-H_{x}u(t;\
u_{0},\theta)\big{\|}^{2}dt+R(\theta)$

A similar approach was applied in [Ayed et al., 2019], but where initial
conditions were learnt as outputs of an additional DNN encoder network that
maps observation sequences (of fixed length and temporal discretization) to
initial conditions.

####  4.3.2. Optimization; Weak Constraint For Missing Dynamics

The hard constraint minimization is very sensitive for large $T$ in settings
where the dynamics is chaotic. This can be ameliorated, to some extent, by
considering the objective function

(4.6) $\displaystyle\mathcal{J}_{T}\big{(}\theta,u(t)\big{)}$
$\displaystyle=\frac{1}{T}\int_{0}^{T}\big{\|}z(t)-H_{x}u(t)\big{\|}^{2}dt$
$\displaystyle\quad+\lambda\frac{1}{T}\int_{0}^{T}\big{\|}\dot{u}(t)-f(u(t);\
\theta)\big{\|}^{2}dt.$

This objective function is employed in [Ouala et al., 2020], where it is
motivated by the weak-constraint variational formulation (4DVAR) arising in
data assimilation Law et al. [2015].

####  4.3.3. Optimization; Data Assimilation For Missing Dynamics

The weak constraint approach may still scale poorly with $T$ large, and still
relies on gradient-based optimization to infer hidden states. To avoid these
potential issues, we follow the recent work of [Chen et al., 2021b], using
filtering-based methods to estimate the hidden state. This implicitly learns
initializations and it removes noise from data. It allows computation of
gradients of the resulting loss function back through the filtering algorithm
to learn model parameters. We define a filtered state

$\hat{u}_{t,\tau}:=\hat{u}_{t}\Big{(}\tau;\
\hat{v},\theta_{\textrm{DYN}},\theta_{\textrm{DA}},\big{\\{}z(t+s)\big{\\}}_{s=0}^{\tau}\Big{)}$

as an estimate of $u(t+\tau)|\\{z(t+s)\\}_{s=0}^{\tau}$ when initialized at
$\hat{u}_{t,0}=\hat{v}$. 333In practice we have found that setting $\hat{v}=0$
works well. In this formulation, we distinguish $\theta_{\textrm{DYN}}$ as
parameters for modeling dynamics via (4.3), and $\theta_{\textrm{DA}}$ as
hyper-parameters governing the specifics of a data assimilation scheme.
Examples of $\theta_{\textrm{DA}}$ are the constant gain matrix $K$ that must
be chosen for 3DVAR, or parameters of the inflation and localization methods
deployed within Ensemble Kalman Filtering. By parameterizing these choices as
$\theta_{\textrm{DA}}$, we can optimize them jointly with model parameters
$\theta_{\textrm{DYN}}$. To do this, let
$\theta=[\theta_{\textrm{DYN}},\theta_{\textrm{DA}}]$ and minimize

(4.7) $\displaystyle\mathcal{J}_{T}(\theta)$
$\displaystyle=\frac{1}{(T-\tau_{1}-\tau_{2})\tau_{2}}\int_{t=0}^{T-\tau_{1}-\tau_{2}}\int_{s=0}^{\tau_{2}}\big{\|}z(t+\tau_{1}+s)-H_{x}u(s;\
\hat{u}_{t,\tau_{1}},\theta)\big{\|}^{2}ds\ dt.$

Here, $\tau_{1}$ denotes the length of assimilation time used to estimate the
state which initializes a parameter-fitting over window of duration
$\tau_{2};$ this parameter-fitting leads to the inner-integration over $s$.
This entire picture is then translated through $t$ time units and the
objective function is found by integrating over $t$. Optimizing (4.7) can be
understood as a minimization over short-term forecast errors generated from
all assimilation windows. The inner integral takes a fixed start time $t$,
applies data assimilation over a window $[t,t+\tau_{1}]$ to estimate an
initial condition $\hat{u}_{t,\tau_{1}}$, then computes a short-term
($\tau_{2}$) prediction error resulting from this DA-based initialization. The
outer integral sums these errors over all available windows in long trajectory
of data of length $T$.

In our work, we perform filtering using a simple 3DVAR method, whose constant
gain can either be chosen as constant, or can be learnt from data. When
constant, a natural choice is $K\propto H_{x}^{T}$, and this approach has a
direct equivalence to standard warmup strategies employed in RNN and RC
training [Vlachas et al., 2020, Pathak et al., 2018b]. The paper [Chen et al.,
2021b] suggests minimization of a similar objective, but considers more
general observation operators $h$, restricts the outer integral to non-
overlapping windows, and solves the filtering problem with an EnKF with known
state-covariance structure.

######  Remark 4.1.

To motivate learning parameters of the data assimilation we make the following
observation: for problems in which the model is known (i.e.
$\theta_{\textrm{DYN}}$ is fixed) we observe successes with the approach of
identifying 3DVAR gains that empirically outperform the theoretically derived
gains in [Law et al., 2013]. Similar is to be expected for parameters defining
inflation and localization in the EnKF. $\square$

######  Remark 4.2.

Specific functional forms of $f_{1},f_{2}$ (and their corresponding parameter
inference strategies) reduce (4.3) to various approaches. For the continuous-
time RNN analysis that we discuss in Section 5 we will start by considering
settings in which $f_{1}$ and $f_{2}$ are approximated from expressive
function classes, such as neural networks. We will then specify to models in
which $f_{1}$ is linear in $r$ and independent of $x$, whilst $f_{2}$ is a
single layer neural network. It is intuitive that the former may be more
expressive and allow a smaller $d_{r}$ than the latter; but the latter
connects directly to reservoir computing, a connection we make explicitly in
what follows. Our numerical experiments in Section 6 will be performed in both
settings: we will train models from the more general setting; and by carefully
designed experiments we will shed light on issues arising from over-
parameterization, in the sense of choosing to learn a model in dimension
higher than that of the true observed-hidden model, working in the setting of
linear coupling term $f_{1},$ depending only on $r$. $\square$

######  Remark 4.3.

The recent paper [Gupta and Lermusiaux, 2021] proposes an interesting, and
more computationally tractable, approach to learning model error in the
presence of memory. They propose to learn a closure operator $m_{\tau}(\cdot\
;\ \theta)\colon\mathsf{X}_{\tau}\to\mathbbm{R}^{d_{x}}$ for a DDE with finite
memory $\tau$:

(4.8)
$\dot{x}(t)=f_{0}\big{(}x(t)\big{)}+m_{\tau}\big{(}\\{x(t-s)\\}_{s=0}^{\tau};\
\theta\big{)};$

neural networks are used to learn the operator $m_{\tau}.$ Alternatively,
Gaussian processes are used to fit a specific class of stochastic delay
differential equation (SDDE) (4.8) in [Schneider et al., 2021]. However,
although delay-based approaches have seen some practical success, in many
applications they present issues for domain interpretability and Markovian ODE
or PDE closures are more desirable. $\square$

###  4.4. Discrete-Time Non-Markovian Learning

In similar spirit to Section 4.3, we can aim to recreate discrete-time
dynamics in $x$ for (2.12) with model

(4.9a) $\displaystyle x_{k+1}$
$\displaystyle=\Psi_{0}(x_{k})+\Psi_{1}(r_{k},x_{k};\ \theta)$ (4.9b)
$\displaystyle r_{k+1}$ $\displaystyle=\Psi_{2}(r_{k},x_{k};\ \theta)$

and objective function

(4.10) $\displaystyle\mathcal{J}_{T}(\theta,r_{0})$
$\displaystyle=\frac{1}{K}\sum_{k=0}^{K-1}\big{\|}x_{k+1}-\Psi_{0}(x_{k})-\Psi_{1}(r_{k},x_{k};\
\theta)\big{\|}^{2}+R(\theta)$
$\displaystyle\textrm{s.t.}\quad\\{r_{k}\\}_{k=1}^{K-1}\quad\textrm{solves
\eqref{eq:ml_disc_mem2}.}$

Observe that estimation of initial condition $r_{0}$ is again crucial, and the
data assimilation methods discussed in Section 4.3 can be adapted to this
discrete-time setting. The functional form of $\Psi_{1},\Psi_{2}$ (and their
corresponding parameter inference strategies) reduce (4.9) to various
approaches, including recurrent neural networks, latent ODEs, and delay-
embedding maps (e.g. to get a delay embedding map, $\Psi_{2}$ is a shift
operator). Pathak et al. [2018b] use reservoir computing (a random features
analog to RNN, described in the next section) with $L_{2}$ regularization to
study an approach similar to (4.9), but included $\Psi_{0}(x_{k})$ as a
feature in $\Psi_{1}$ and $\Psi_{2}$ instead of using it as the central model
upon which to learn residuals. The data-driven super-parameterization approach
in [Chattopadhyay et al., 2020b] also appears to follow the underlying
assumption of (4.9). Harlim et al. [2021] evaluate hybrid models of form (4.9)
both in settings where delay embedding closures are employed and where RNN-
based approximations via LSTMs are employed.

##  5\. Underpinning Theory

In this section we identify specific hypothesis classes $\mathcal{M}.$ We do
this using random feature maps [Rahimi and Recht, 2008a] in the Markovian
settings (Section 5.1), and using recurrent neural networks in the memory-
dependent setting. We then discuss these problems from a theoretical
standpoint. In Section 5.2 we study excess risk and generalization error in
the context of linear models (a setting which includes the random features
model as a special case). And we conclude by discussing the use of RNNs
[Goodfellow et al., 2016][Chapter 10] for the non-Markovian settings
(discrete- and continuous-time) in Section 5.3; we present an approximation
theorem for continuous-time hybrid RNN models. Throughout this section, the
specific use of random feature maps and of recurrent neural networks is for
illustration only; other models could, of course, be used.

###  5.1. Markovian Modeling with Random Feature Maps

In principle, any hypothesis class can be used to learn $m^{\dagger}$.
However, we focus on models that are easily trained on large-scale complex
systems and yet have proven approximation power for functions between finite-
dimensional Euclidean spaces. For the Markovian modeling case, we use random
feature maps; like traditional neural networks, they possess arbitrary
approximation power [Rahimi and Recht, 2008c, b], but further benefit from a
quadratic minimization problem in the training phase, as do kernel or Gaussian
process methods. In our case studies, we found random feature models
sufficiently expressive, we found optimization easily implementable, and we
found the learned models generalized well. Moreover, their linearity with
respect to unknown parameters enables a straightforward analysis of excess
risk and generalization error in Section 5.2. Details on the derivation and
specific design choices for our random feature modeling approach can be found
in Section 8.4, where we explain how we sample $D$ random feature functions
$\varphi:\mathbbm{R}^{d_{x}}\to\mathbbm{R}$ and stack them to form a vector-
valued feature map $\phi\colon\mathbbm{R}^{d_{x}}\to\mathbbm{R}^{D}$. Given
this random function $\phi$, we define the hypothesis class

(5.1) $\mathcal{M}=\\{m\colon\mathbbm{R}^{d_{x}}\to\mathbbm{R}^{d_{x}}\ |\
\exists\ C\in\mathbbm{R}^{d_{x}\times D}:m(x)=C\phi(x)\\}.$

####  5.1.1. Continuous-Time

In the continuous-time framing, our Markovian closure model uses hypothesis
class (5.1) and thus takes the form

$\displaystyle\dot{x}$ $\displaystyle=f_{0}(x)+C\phi(x(t)).$

We rewrite (4.1) for this particular case with an $L_{2}$ regularization
parameter $\lambda\in\mathbbm{R}^{+}$:

(5.2)
$\mathcal{J}_{T}(C)=\frac{1}{2T}\int_{0}^{T}\left\|\dot{x}(t)-f_{0}(x(t))-C\phi\big{(}x(t)\big{)}\right\|^{2}dt+\frac{\lambda}{2}\|C\|^{2}.$

We employ the notation $A\otimes B:=AB^{T}$ for the outer-product between
matrices $A\in\mathbbm{R}^{m\times n},B\in\mathbbm{R}^{l\times n}$, and the
following notation for time-average

$\overline{A}_{T}:=\frac{1}{T}\int_{0}^{T}A(t)dt$

of $A\in L^{1}([0,T],\mathbbm{R}^{m\times n}).$ The objective function in
(5.2) is quadratic and convex in $C$ and thus is globally minimized for the
unique $C^{\ast}$ which makes the derivative of $\mathcal{J}_{T}$ zero.
Consequently, the minimizer $C^{\ast}$ satisfies the following linear equation
(derived in Section 8.5):

(5.3) $(Z+\lambda I)(C^{\ast})^{T}=Y.$

Here, $I\in\mathbbm{R}^{D\times D}$ is the identity and

(5.4) $\displaystyle Z$
$\displaystyle=\overline{[\phi\otimes\phi]}_{T}\in\mathbbm{R}^{D\times D},$
$\displaystyle Y$ $\displaystyle=\overline{[\phi\otimes
m^{\dagger}]}_{T}\in\mathbbm{R}^{D\times d_{x}}.$

Of course $m^{\dagger}$ is not known, but
$m^{\dagger}(t)=\dot{x}(t)-f_{0}(x(t))$ can be computed from data.

To summarize, the algorithm proceeds as follows: 1) create a realization of
random feature vector $\phi$; 2) compute the integrals in (5.4) to obtain
$Z,Y$; and 3) solve the linear matrix equation (5.3) for $C^{\ast}$. Together
this leads to our approximation $m^{\dagger}(x)\approx m^{\ast}_{T}(x;\
\theta):=C^{\ast}\phi(x)$.

####  5.1.2. Discrete-Time

In discrete-time, our Markovian closure model is

$\displaystyle x_{k+1}$ $\displaystyle=\Psi_{0}(x_{k})+C\phi(x_{k}),$

and is learnt by minimizing

(5.5)
$\mathcal{J}_{T}(\theta)=\frac{1}{K}\sum_{k=0}^{K-1}\big{\|}x_{k+1}-\Psi_{0}(x_{k})-C\phi\big{(}x(t)\big{)}\big{\|}^{2}+\frac{\lambda}{2}\|C\|^{2}.$

The objective function in (5.5) is quadratic in $C$ and thus globally
minimized at $C^{\ast}$. As in Section 5.1.1, we can compute $Z,Y$ and solve a
linear system for $C^{\ast}$ to approximate $m^{\dagger}(x)\approx
m^{\ast}_{T}(x;\ \theta):=C^{\ast}\phi(x)$. This formulation closely mirrors
the fully data-driven linear regression approach in [Gottwald and Reich,
2021b].

###  5.2. Learning Theory for Markovian Models with Linear Hypothesis Class

In this subsection we provide estimates of the excess risk and generalization
error in the context of learning $m^{\dagger}$ in (2.2) from a trajectory over
time horizon $T$. We study ergodic continuous-time models in the setting of
Section 4.1. To this end we consider the very general linear hypothesis class
given by

(5.6) $\mathcal{M}=\\{m\colon\mathbbm{R}^{d_{x}}\to\mathbbm{R}^{d_{x}}\ |\
\exists\
\theta\in\mathbbm{R}^{p}:m(x)=\sum_{\ell=1}^{p}\theta_{\ell}f_{\ell}(x)\\};$

we note that if the $\\{f_{\ell}\\}$ are i.i.d. draws of function $\phi$ in
the case $D=d_{x}$ then this too reduces to a random features model, but that
our analysis in the context of statistical learning does not rely on the
random features structure. In fact our analysis can be used to provide
learning theory for other linear settings, where $\\{f_{\ell}\\}$ represents a
dictionary of hypothesized features whose coefficients are to be learnt from
data. Nonetheless, universal approximation for random features [Rahimi and
Recht, 2008a] provides an important example of an approximation class for
which the loss function $\mathcal{I}_{\infty}$ may be made arbitrarily small
by choice of $p$ large enough and appropriate choice of parameters, and the
reader may find it useful to focus on this case. We also note that the theory
we present in this subsection is readily generalized to working with
hypothesis class (5.1).

We make the following ergodicity assumption about the data generation process:

######  Assumption A1.

Equation (2.2) possesses a compact attractor $\mathcal{A}$ supporting
invariant measure $\mu.$ Furthermore the dynamical system on $\mathcal{A}$ is
ergodic with respect to $\mu$ and satisfies a central limit theorem of the
following form: for all Hölder continuous
$\varphi:\mathbb{R}^{d_{x}}\mapsto\mathbb{R}$, there is
$\sigma^{2}=\sigma^{2}(\varphi)$ such that

(5.7)
$\sqrt{T}\Biggl{(}\frac{1}{T}\int_{0}^{T}\varphi\bigl{(}x(t)\bigr{)}dt-\int_{\mathbb{R}^{d_{x}}}\varphi\bigl{(}x\bigr{)}\mu(dx)\Biggr{)}\Rightarrow
N(0,\sigma^{2})$

where $\Rightarrow$ denotes convergence in distribution with respect to
$x(0)\sim\mu$. Furthermore a law of the iterated logarithm holds: almost
surely with respect to $x(0)\sim\mu$,

(5.8) ${\rm limsup}_{T\to\infty}\Bigl{(}\frac{T}{\log\log
T}\Bigr{)}^{\frac{1}{2}}\Biggl{(}\frac{1}{T}\int_{0}^{T}\varphi\bigl{(}x(t)\bigr{)}dt-\int_{\mathbb{R}^{d_{x}}}\varphi\bigl{(}x\bigr{)}\mu(dx)\Biggr{)}=\sigma.$

######  Remark 5.1.

Note that in both (5.7) and (5.8) $\varphi(\cdot)$ is only evaluated on
(compact) $\mathcal{A}$ obviating the need for any boundedness assumptions on
$\varphi(\cdot)$. In the work of Melbourne and co-workers, Assumption A1 is
proven to hold for a class of differential equations, including the Lorenz ’63
model at, and in a neighbourhood of, the classical parameter values: in
[Holland and Melbourne, 2007] the central limit theorem is established; and in
[Bahsoun et al., 2020] the continuity of $\sigma$ in $\varphi$ is proven.
Whilst it is in general very difficult to prove such results for any given
chaotic dynamical system, there is strong empirical evidence for such results
in many chaotic dynamical systems that arise in practice. This combination of
theory and empirical evidence justify studying the learning of model error
under Assumption A1. Tran and Ward [2017] were the first to make use of the
theory of Melbourne and coworkers to study learning of chaotic differential
equations from time-series. $\square$

Given $m$ from hypothesis class $\mathcal{M}$ defined by (5.6) we define

(5.9) $\theta^{\ast}_{\infty}={\rm
argmin}_{\theta\in\mathbb{R}^{p}}\mathcal{I}_{\infty}\bigl{(}m(\cdot\
;\theta)\bigr{)}\ ={\rm
argmin}_{\theta\in\mathbb{R}^{p}}\mathcal{L}_{\mu}\bigl{(}m(\cdot\
;\theta)\bigr{)}$

and

(5.10) $\theta^{\ast}_{T}={\rm
argmin}_{\theta\in\mathbb{R}^{p}}\mathcal{I}_{T}\bigl{(}m(\cdot\
;\theta)\bigr{)}.$

(Regularization is not needed in this setting because the data is plentiful—a
continuous-time trajectory—and the number of parameters is finite.) Then
$\theta^{\ast}_{\infty},\theta^{\ast}_{T}$ solve the linear systems

$A_{\infty}\theta^{\ast}_{\infty}=b_{\infty},\quad
A_{\scriptsize{T}}\theta^{\ast}_{T}=b_{\scriptsize{T}}$

where

$\displaystyle(A_{\infty})_{ij}$
$\displaystyle=\int_{\mathbbm{R}^{d_{x}}}\bigl{\langle}f_{i}(x),f_{j}(x)\bigr{\rangle}\,\mu(dx),\quad$
$\displaystyle(b_{\infty})_{j}$
$\displaystyle=\int_{\mathbbm{R}^{d_{x}}}\bigl{\langle}m^{\dagger}(x),f_{j}(x)\bigr{\rangle}\,\mu(dx),$
$\displaystyle(A_{\scriptsize{T}})_{ij}$
$\displaystyle=\frac{1}{T}\int_{0}^{T}\bigl{\langle}f_{i}\bigl{(}x(t)\bigr{)},f_{j}\bigl{(}x(t)\bigr{)}\bigr{\rangle}\,dt,\quad$
$\displaystyle(b_{\scriptsize{T}})_{j}$
$\displaystyle=\frac{1}{T}\int_{0}^{T}\bigl{\langle}m^{\dagger}\bigl{(}x(t)\bigr{)},f_{j}\bigl{(}x(t)\bigr{)}\bigr{\rangle}\,dt.$

These facts can be derived analogously to the derivation in Section 8.5. Given
$\theta^{\ast}_{\infty}$ and $\theta^{\ast}_{T}$ we also define

$m^{\ast}_{\infty}=m(\cdot\ ;\theta^{\ast}_{\infty}),\,\,m^{\ast}_{T}=m(\cdot\
;\theta^{\ast}_{T}).$

Recall that it is assumed that $f^{\dagger},f_{0}$, and $m^{\dagger}$ are
$C^{1}.$ We make the following assumption regarding the vector fields defining
hypothesis class $\mathcal{M}.$

######  Assumption A2.

The functions $\\{f_{\ell}\\}_{\ell=0}^{p}$ appearing in definition (5.6) of
the hypothesis class $\mathcal{M}$ are Hölder continuous on
$\mathbbm{R}^{d_{x}}$. In addition, the matrix $A_{\infty}$ is invertible.

######  Theorem 5.2.

Let Assumptions A1 and A2 hold. Then the scaled excess risk $\sqrt{T}R_{T}$ in
(3.2) (resp. scaled generalization error $\sqrt{T}|G_{T}|$ in (3.3)) is
bounded above by $\|\mathcal{E}_{R}\|$ (resp. $\|\mathcal{E}_{G}\|$), where
random variable $\mathcal{E}_{R}\in\mathbb{R}^{p}$ (resp.
$\mathcal{E}_{G}\in\mathbb{R}^{p+1}$) converges in distribution to
$N(0,\Sigma_{R})$ (resp. $N(0,\Sigma_{G}$)) w.r.t. $x(0)\sim\mu$ as
$T\to\infty.$ Furthermore, there is constant $C>0$ such that, almost surely
w.r.t. $x(0)\sim\mu$,

${\rm limsup}_{T\to\infty}\Bigl{(}\frac{T}{\log\log
T}\Bigr{)}^{\frac{1}{2}}\bigl{(}R_{T}+|G_{T}|\bigr{)}\leq C.$

The proof is provided in Section 8.1.

######  Remark 5.3.

The convergence in distribution shows that, with high probability with respect
to initial data, the excess risk and the generalization error are bounded
above by terms of size $1/\sqrt{T}.$ This can be improved to give an almost
sure result, at the cost of the factor of $\sqrt{\log\log T}$. The theorem
shows that (ignoring log factors and acknowledging the probabilistic nature of
any such statements) trajectories of length ${\mathcal{O}(\epsilon^{-2})}$ are
required to produce bounds on the excess risk and generalization error of size
${\mathcal{O}(\epsilon)}$.

The bounds on excess risk and generalization error also show that empirical
risk minimization (of $\mathcal{I}_{T}$) approaches the theoretically
analyzable concept of risk minimization (of $\mathcal{I}_{\infty}$) over
hypothesis class (5.6). The sum of the excess risk $R_{T}$ and the
generalization error $G_{T}$ gives

$E_{T}:=\mathcal{I}_{T}(m^{\ast}_{T})-\mathcal{I}_{\infty}(m^{\ast}_{\infty}).$

We note that $\mathcal{I}_{T}(m^{\ast}_{T})$ is computable, once the
approximate solution $m^{\ast}_{T}$ has been identified; thus, when combined
with an estimate for $E_{T}$, this leads to an estimate for the risk
associated with the hypothesis class used.

If the approximating space $\mathcal{M}$ is rich enough, then approximation
theory may be combined with Theorem 5.2 to estimate the trajectory error
resulting from the learned dynamical system. Such an approach is pursued in
Proposition 3 of [Zhang et al., 2020] for SDEs. Furthermore, in that setting,
knowledge of rate of mixing/decay of correlations for SDEs may be used to
quantify constants appearing in the error bounds. It would be interesting to
pursue such an analysis for chaotic ODEs with known mixing rates/decay of
correlations. Such results on mixing are less well-developed, however, for
chaotic ODEs; see discussion of this point in [Holland and Melbourne, 2007],
and the recent work [Bahsoun et al., 2020].

Work by Zhang et al. [2021] demonstrates that error bounds on learned model
error terms can be extended to bound error on reproduction of invariant
statistics for ergodic SDEs. Moreover, E et al. [2019] provide a direction for
proving similar bounds on model error learning using nonlinear function
classes (e.g. two-layer neural networks).

Finally we remark on the dependence of the risk and generalization error
bounds on the size of the model error. It is intuitive that the amount of data
required to learn model error should decrease as the size of the model error
decreases. This is demonstrated numerically in Section 6.2.3 (c.f. Figures
2(a) and 2(b)). Here we comment that Theorem 5.2 also exhibits this feature:
examination of the proof in Appendix 8.1 shows that all upper bounds on terms
appearing in the excess and generalization error are proportional to
$m^{\dagger}$ itself or to $m_{\infty}^{*}$, its approximation given an
infinite amount of data; note that $m_{\infty}^{*}=m^{\dagger}$ if the
hypothesis class contains the truth. $\square$

###  5.3. Non-Markovian Modeling with Recurrent Neural Networks

Recurrent Neural Networks (RNNs) are one of the _de facto_ tools for modeling
systems with memory. Here, we show straightforward residual implementations of
RNNs for continuous- and discrete-time, with the goal of modeling non-
Markovian model error.

####  5.3.1. General Case

Equation (4.3b), and its coupling to (4.3a), constitute a very general way to
account for memory-dependent model error in the dynamics of $x$. In fact, for
$f_{1},f_{2}$ sufficiently expressive (e.g. random feature functions, neural
networks, polynomials), and $d_{r}\geq d_{y}$, solutions to (4.3) can
approximate solutions to (2.8) arbitrarily well. We make this type of
universal approximation theorem concrete in Theorems 5.4 and 5.6. We start by
proving Theorem 5.4, which rests on the following assumptions:

######  Assumption A3.

Functions $f^{\dagger},g^{\dagger},f_{0},f_{1},f_{2}$ are all globally
Lipschitz.

Note that this implies that $m^{\dagger}$ is also globally Lipschitz.

######  Assumption A4.

Fix $T>0$. There exist $\rho_{0}\in\mathbbm{R},\rho_{T}\in\mathbbm{R}$ such
that, for equation (2.8), $\bigl{(}x(0),y(0)\bigr{)}\in B(0,\rho_{0})$ implies
that $\bigl{(}x(t),y(t)\bigr{)}\in B(0,\rho_{T})$ $\forall\ t\in[0,T]$.

######  Assumption A5.

The hidden state in (4.3), $r\in\mathbbm{R}^{d_{r}}$, has the same dimension
as the true hidden state $y$ in (2.8); that is $d_{r}=d_{y}$.

######  Assumption A6.

Let functions $f_{1}(\cdot\ ;\ \theta)\in
C^{1}(\mathbbm{R}^{d_{x}}\times\mathbbm{R}^{d_{y}};\ \mathbbm{R}^{d_{x}})$ and
$f_{2}(\cdot\ ;\ \theta)\in
C^{1}(\mathbbm{R}^{d_{x}}\times\mathbbm{R}^{d_{y}};\ \mathbbm{R}^{d_{y}})$ be
parameterized 444Here we define $\mathbb{N}=\\{1,2,\dots,\\}$, the strictly
positive integers. by $n\in\mathbb{N}$ and $\theta\in\mathbbm{R}^{n}$. Then,
for any $\delta>0$, there exists $n>0$ and $\theta\in\mathbbm{R}^{n}$ such
that

$\sup_{x,y\in B(0,\rho_{T})}\|f^{\dagger}(x,y)-f_{1}(x,y;\
\theta)\|\leq\delta$

and

$\sup_{x,y\in B(0,\rho_{T})}\|g^{\dagger}(x,y)-f_{2}(x,y;\
\theta)\|\leq\delta$

Note that A6 can be satisfied by any parametric function class possessing a
universal approximation property for maps between finite-dimensional Euclidean
spaces, such as neural networks, polynomials and random feature methods. The
next theorem transfers this universal approximation property for maps between
Euclidean spaces to a universal approximation property for representation of
model error with memory; this is a form of infinite dimensional approximation
since, via its own dynamics, the memory variable $r$ maps the past history of
$x$ into the model error correction term in the dynamics for $x$.

######  Theorem 5.4.

Let Assumptions A3-A6 hold. Fix any $T>0$ and $\rho_{0}>0,$ let
$x(\cdot),y(\cdot)$ denote the solution of (2.8) with $\varepsilon=1$ and let
$x_{\delta}(\cdot),r_{\delta}(\cdot)$ denote the solution of (4.3) with
parameters $\theta\in\mathbb{R}^{n}$. Then, for any $\delta>0$ and any $T>0$,
there is a parameter dimension $n=n_{\delta}\in\mathbbm{N}$ and
parameterization $\theta=\theta_{\delta}\in\mathbbm{R}^{n_{\delta}}$ with the
property that, for any initial condition $\bigl{(}x(0),y(0)\bigr{)}\in
B(0,\rho_{0})$ for (2.8), there is initial condition
$(x_{\delta}(0),r_{\delta}(0))\in\mathbbm{R}^{d_{x}+d_{y}}$ for (4.3), such
that

$\sup_{t\in[0,T]}\|x-x_{\delta}\|\leq\delta.$

The proof is provided in Section 8.2; it is a direct consequence of the
approximation power of $f_{1},f_{2}$ and the Gronwall Lemma.

######  Remark 5.5.

Note that this existence theorem also holds for $d_{r}>d_{y}$ by freezing the
dynamics in the excess dimensions and initializing it at, for example, $0$.
However it is possible for augmentations with $d_{r}>d_{y}$ to introduce
numerical instability when imperfectly initialized in the excess dimensions,
despite their provable correctness when perfectly initialized (see Section
6.4). Nevertheless, we did not encounter such issues when training the general
model class on the examples considered in this paper – see Section 6.3).
$\square$

####  5.3.2. Linear Coupling

We now study a particular form RNN in which the coupling term $f_{1}$
appearing in (4.3) is linear and depends only on the hidden variable:

(5.11a) $\displaystyle\dot{x}$ $\displaystyle=f_{0}(x)+Cr$ (5.11b)
$\displaystyle\dot{r}$ $\displaystyle=\sigma(Ar+Bx+c).$

Here $\sigma$ is an activation function. The specific linear coupling form is
of particular interest because of the connection we make (see Remark 5.9
below) to reservoir computing. The goal is to choose $A,B,C,c$ so that output
$\\{x(t)\\}_{t\geq 0}$ matches output of (2.8), without observation of
$\\{y(t)\\}_{t\geq 0}$ or knowledge of $m^{\dagger}$ and $g^{\dagger}.$ As in
the general case from the preceding subsection, inherent in choosing these
matrices $A,B,C$ and vector $c$ is a choice of embedding dimension for
variable $r$ which will typically be larger than dimension of $y$ itself. The
idea is to create a recurrent state $r$ of sufficiently large dimension
$d_{r}$ whose evolution equation takes $x$ as input and, after a final linear
transformation, approximates the missing dynamics $m^{\dagger}(x,y).$

There is existing approximation theory for discrete-time RNNs [Schäfer and
Zimmermann, 2007] showing that a discrete-time analog of our linear coupling
set-up can be used to approximate discrete-time systems arbitrarily well; see
also Theorem 3 of [Harlim et al., 2021]. There is also a general approximation
theorem using continuous-time RNNs proved in [Funahashi and Nakamura, 1993],
but it does not apply to the linear-coupling setting. We thus extend the work
in these three papers to the context of residual-based learning as in (5.11).
We state the theorem after making three assumptions upon which it rests:

######  Assumption A7.

Functions $f^{\dagger},g^{\dagger},f_{0}$ are all globally Lipschitz.

Note that this implies that $m^{\dagger}$ is also globally Lipschitz.

######  Assumption A8.

Let $\sigma_{0}\in C^{1}(\mathbbm{R};\mathbbm{R})$ be bounded and monotonic,
with bounded first derivative. Then $\sigma(u)$ defined by
$\sigma(u)_{i}=\sigma_{0}(u_{i})$ satisfies $\sigma\in
C^{1}(\mathbbm{R}^{p};\mathbbm{R}^{p}).$

######  Assumption A9.

Fix $T>0$. There exist $\rho_{0}\in\mathbbm{R},\rho_{T}\in\mathbbm{R}$ such
that, for equation (2.8), $\bigl{(}x(0),y(0)\bigr{)}\in B(0,\rho_{0})$ implies
that $\bigl{(}x(t),y(t)\bigr{)}\in B(0,\rho_{T})$ $\forall\ t\in[0,T]$.

######  Theorem 5.6.

Let Assumptions A7-A9 hold. Fix any $T>0$ and $\rho_{0}>0,$ let
$x(\cdot),y(\cdot)$ denote the solution of (2.8) with $\varepsilon=1$ and let
$x_{\delta}(\cdot),r_{\delta}(\cdot)$ denote the solution of (5.11) with
parameters $\theta\in\mathbb{R}^{n}$. Then, for any $\delta>0$ and any $T>0$,
there is embedding dimension $d_{r}\in\mathbbm{N}$, parameter dimension
$n=n_{\delta}\in\mathbbm{N}$ and parameterization
$\theta=\theta_{\delta}=\\{A_{\delta},\ B_{\delta},\ C_{\delta},c_{\delta}\\}$
with the property that, for any initial condition
$\bigl{(}x(0),y(0)\bigr{)}\in B(0,\rho_{0})$ for (2.8), there is initial
condition $(x_{\delta}(0),r_{\delta}(0))\in\mathbbm{R}^{d_{x}+d_{r}}$ for
(5.11), such that

$\sup_{t\in[0,T]}\|x-x_{\delta}\|\leq\delta.$

The complete proof is provided in Section 8.3; here we describe its basic
structure. Define $m(t):=m^{\dagger}\bigl{(}x(t),y(t)\bigr{)}$ and, with the
aim of finding a differential equation for $m(t)$, recall (2.8) with
$\varepsilon=1$ and define the vector field

(5.12)
$h^{\dagger}(x,y):=\nabla_{x}m^{\dagger}(x,y)[f_{0}(x)+m^{\dagger}(x,y)]+\nabla_{y}m^{\dagger}(x,y)g^{\dagger}(x,y).$

Since $\dot{m}(t)$ is the time derivative of
$m^{\dagger}\big{(}x(t),y(t)\big{)}$, when $(x,y)$ solve (2.8) we have

$\dot{m}=h^{\dagger}(x,y).$

Motivated by these observations, we now introduce a new system of autonomous
ODEs for the variables
$(x,y,m)\in\mathbbm{R}^{d_{x}}\times\mathbbm{R}^{d_{y}}\times\mathbbm{R}^{d_{x}}$:

(5.13a) $\displaystyle\dot{x}$ $\displaystyle=f_{0}(x)+m$ (5.13b)
$\displaystyle\dot{y}$ $\displaystyle=g^{\dagger}(x,y)$ (5.13c)
$\displaystyle\dot{m}$ $\displaystyle=h^{\dagger}(x,y).$

To avoid a proliferation of symbols we use the same letters for $(x,y)$
solving equation (5.13) as for $(x,y)$ solving equation (2.8). We now show
$m=m^{\dagger}(x,y)$ is an invariant manifold for (5.13); clearly, on this
manifold, the dynamics of $(x,y)$ governed by (5.13) reduces to the dynamics
of $(x,y)$ governed by (2.8). Thus $m(t)$ must be initialized at
$m^{\dagger}\bigl{(}x(0),y(0)\bigr{)}$ to ensure equivalence between the
solution of (5.13) and (2.8).

The desired invariance of manifold $m=m^{\dagger}(x,y)$ under the dynamics
(5.13) follows from the identity

(5.14)
$\frac{d}{dt}\Bigl{(}m-m^{\dagger}(x,y)\Bigr{)}=-\nabla_{x}m^{\dagger}(x,y)\bigl{(}m-m^{\dagger}(x,y)\bigr{)}.$

The identity is derived by noting that, recalling (5.12) for the definition of
$h^{\dagger}$, and then using (5.13),

$\displaystyle\frac{d}{dt}m$ $\displaystyle=h^{\dagger}(x,y)$
$\displaystyle=\nabla_{x}m^{\dagger}(x,y)[f_{0}(x)+m^{\dagger}(x,y)]+\nabla_{y}m^{\dagger}(x,y)g^{\dagger}(x,y)$
$\displaystyle=\nabla_{x}m^{\dagger}(x,y)[f_{0}(x)+m)]+\nabla_{y}m^{\dagger}(x,y)g^{\dagger}(x,y)$
$\displaystyle\quad\quad\quad-\nabla_{x}m^{\dagger}(x,y)\bigl{(}m-m^{\dagger}(x,y)\bigr{)}$
$\displaystyle=\frac{d}{dt}m^{\dagger}(x,y)-\nabla_{x}m^{\dagger}(x,y)\bigl{(}m-m^{\dagger}(x,y)\bigr{)}.$

We emphasize this calculation is performed under the dynamics defined by
(5.13).

The proof of the RNN approximation property proceeds by approximating vector
fields $g^{\dagger}(x,y),h^{\dagger}(x,y)$ by neural networks and introducing
linear transformations of $y$ and $m$ to rewrite the approximate version of
system (5.13) in the form (5.11). The effect of the approximation of the
vector fields on the true solution is then propagated through the system and
its effect controlled via a straightforward Gronwall argument.

######  Remark 5.7.

The details of training continuous-time RNNs to ensure accuracy and long-time
stability are a subject of current research [Chang et al., 2019, Erichson et
al., 2020, Ouala et al., 2020, Chen et al., 2021b] and in this paper we
confine the training of RNNs to an example in the general setting, and not the
case of linear coupling. Discrete-time RNN training, on the other hand, is
much more mature, and has produced satisfactory accuracy and stability for
settings with uniform sample rates that are consistent across train and
testing scenarios [Harlim et al., 2021]. The form with linear coupling is
widely studied in discrete time models. Furthermore, sophisticated variants on
RNNs, such as Long-Short Term Memory (LSTM) RNNs [Hochreiter and Schmidhuber,
1997] and Gated Recurrent Units (GRU) [Cho et al., 2014a], are often more
effective, although similar in nature RNNs. However, the potential
formulation, implementation and advantages of these variants in the
continuous-time setting [Niu et al., 2019] is not yet understood. We refer
readers to [Goodfellow et al., 2016] for background on discrete RNN
implementations and backpropagation through time (BPTT). For implementations
of continuous-time RNNs, it is common to leverage the success of the automatic
BPTT code written in PyTorch and Tensorflow by discretizing (5.11) with an ODE
solver that is compatible with these autodifferentiation tools (e.g.
torchdiffeq by [Rubanova et al., 2019], NbedDyn by [Ouala et al., 2020], and
AD-ENKF by [Chen et al., 2021b]). This compatibility can also be achieved by
use of explicit Runge-Kutta schemes [Queiruga et al., 2020]. Note that the
discretization of (5.11) can (and perhaps should) be much finer than the data
sampling rate $\Delta t$, but that this requires reliable estimation of
$x(t),\dot{x}(t)$ from discrete data. $\square$

######  Remark 5.8.

The need for data assimilation [Asch et al., 2016, Law et al., 2015, Reich and
Cotter, 2015] to learn the initialization of recurrent neural networks may be
understood as follows. Since $m^{\dagger}$ is not known and $y$ is not
observed (and in particular $y(0)$ is not known) the desired initialization
for (5.13), and thus also for approximations of this equation in which
$g^{\dagger}$ and $h^{\dagger}$ are replaced by neural networks, is not known.
Hence, if an RNN is trained on a particular trajectory, the initial condition
that is required for accurate approximation of (2.8) from an unseen initial
condition is not known. Furthermore the invariant manifold
$m=m^{\dagger}(x,y)$ may be unstable under numerical approximation. However if
some observations of the trajectory starting at the new initial condition are
used, then data assimilation techniques can potentially learn the
initialization for the RNN and also stabilize the invariant manifold. Ad hoc
initialization methods are common practice [Haykin et al., 2007, Cho et al.,
2014b, Bahdanau et al., 2016, Pathak et al., 2018a], and rely on forcing the
learned RNN with a short sequence of observed data to synchronize the hidden
state. The success of these approaches likely rely on RNNs’ abilities to
emulate data assimilators [Härter and Velho, 2012]; however, a more careful
treatment of the initialization problem may enable substantial advances.
$\square$

######  Remark 5.9.

Reservoir computing (RC) is a variant on RNNs which has the advantage of
leading to a quadratic optimization problem [Jaeger, 2001, Lukoševičius and
Jaeger, 2009, Grigoryeva and Ortega, 2018]. Within the context of the
continuous-time RNN (5.11) they correspond to randomizing $(A,B,c)$ in (5.11b)
and then choosing only parameter $C$ to fit the data. To be concrete, this
leads to

$r(t)=\mathcal{G}_{t}\bigl{(}\\{x(s)\\}_{s=0}^{t};\ r(0),A,B,c\bigr{)};$

here $\mathcal{G}_{t}$ may be viewed as a random function of the path-history
of $x$ upto time $t$ and of the initial condition for $r.$ Then $C$ is
determined by minimizing the quadratic function

$\mathcal{J}_{T}(C)=\frac{1}{2T}\int_{0}^{T}\left\|\dot{x}(t)-f_{0}(x(t))-Cr(t)\right\|^{2}dt+\frac{\lambda}{2}\|C\|^{2}.$

This may be viewed as a random feature approach on the Banach space
$\mathsf{X}_{T}$; the use of random features for learning of mappings between
Banach spaces is studied by Nelsen and Stuart [2020], and connections between
random features and reservoir computing were introduced by Dong et al. [2020].
In the specific setting described here, care will be needed in choosing
probability measure on $(A,B,c)$ to ensure a well-behaved map
$\mathcal{G}_{t}$; furthermore data assimilation ideas [Asch et al., 2016, Law
et al., 2015, Reich and Cotter, 2015] will be needed to learn an appropriate
$r(0)$ in the prediction phase, as discussed in Remark 5.8 for RNNs. $\square$

##  6\. Numerical Experiments

In this section, we present numerical experiments intended to test different
hypotheses about the utility of hybrid mechanistic and data-driven modeling.
We summarize our findings in Section 6.1. We define the overarching
experimental setup in Section 6.2.1, then introduce our criteria for
evaluating model performance in Section 6.2.2. In the Lorenz ’63 (L63)
experiments (Section 6.2.3), we investigate how a simple Markovian random
features model error term can be recovered using discrete and continuous-time
methods, and how those methods scale with the magnitude of error, data
sampling rate, availability of training data, and number of learned
parameters. In the Lorenz ’96 Multiscale (L96MS) experiments (Section 6.2.4),
we take this a step further by learning a Markovian random features closure
term for a scale-separated system, as well as systems with less scale-
separation. As expected, we find that the Markovian closure approach is highly
accurate for a scale-separated regime. We also see that the Markovian closure
has merit even in cases with reduced scale-separation. However, this situation
would clearly benefit from learning a closure term with memory, a topic we
turn to in Section 6.3, where we demonstrate that non-Markovian closure models
can be learnt from noisy, partially observed data; for low-dimensional cases
(e.g. L63), our method of training converges to return a good model with high
short-term accuracy and long-term statistical validity. For higher-dimensional
cases (e.g. L96MS), we find the method to hold promise, but further research
is required in this general area. In Section 6.4, we demonstrate why non-
Markovian closures must be carefully initialized and/or controlled (e.g. via
data assimilation) in order to ensure their long-term stability and short-term
accuracy.

###  6.1. Summary of Findings from Numerical Experiments

  1. (1)

We find that hybrid modeling has better predictive performance than purely
data-driven methods in a wide range of settings (see Figures 2(a) and 2(b) of
Section 6.2.3): this includes scenarios where $f_{0}$ is highly accurate (but
imperfect) and scenarios where $f_{0}$ is highly inaccurate (but nevertheless
faithfully encodes much of the true structure for $f^{\dagger}$).

  2. (2)

We find that hybrid modeling is more data-efficient than purely data-driven
approaches (Figure 3 of Section 6.2.3).

  3. (3)

We find that hybrid modeling is more parameter-efficient than purely data-
driven approaches (Figure 4 of Section 6.2.3).

  4. (4)

Purely data-driven discrete-time modeling can suffer from instabilities in the
small timestep limit $\Delta t\ll 1$; hybrid discrete-time approaches can
alleviate this issue when they are built from an integrator $\Psi_{0}$, as
this will necessarily encode the correct parametric dependence on $\Delta t\ll
1$ (Figure 5 of Section 6.2.3).

  5. (5)

In order to leverage standard supervised regression techniques, continuous-
time methods require good estimates of derivatives $\dot{x}(t)$ from the data.
Figure 5 of Section 6.2.3 quantifies this estimation as a function of data
sample rate.

  6. (6)

Non-Markovian model error can be captured by Markovian terms in scale-
separated cases. Section 6.2.4 demonstrates this quantitatively in Figure 6,
and qualitatively in Figure 7. Beyond the scale-separation limit, Markovian
terms will fail for trajectory forecasting. However, Markovian terms may still
reproduce invariant statistics in dissipative systems (for example, in cases
with short memory-length). Section 6.2.4 demonstrates this quantitatively in
Figure 6; Figure 7 offers intuition for these findings.

  7. (7)

Non-Markovian description of model error is needed to accurately represent
problems where the hidden dynamics is not scale-separated from the observed
dynamics. Section 6.3 shows how partial and noisy observations can be
exploited by augmented ODE models of form (4.3) when the noise and hidden
dynamics are learnt implicitly by auto-differentiable data assimilation. We
observe high-quality reconstruction of the L63 system along its first-
component when choosing a correct (Figure 8) or overly enlarged (Figure 9)
hidden dimension. We also observe promising reconstruction of the L96MS system
in its slow components (Figure 10); however, long-time solutions to the learnt
model exhibited instabilities inconsistent with the true system.

  8. (8)

Non-Markovian models must be carefully initialized, and indeed data
assimilation is needed, in order to ensure accuracy (Section 6.4) of invariant
statistics (Figure 12), long-term stability (Figure 13), and accurate short-
term predictions (Figure 14). We explain observed phenomena in terms of the
properties of the desired lower-dimensional invariant manifold which is
embedded within the higher dimensional system used as the RNN’s basis of
approximation.

###  6.2. Learning Markovian Model Errors from Noise-Free Data

####  6.2.1. Experimental Set-Up

In the Markovian error modeling experiments described in Sections 6.2.3 and
6.2.4, whether using continuous- or discrete-time models, we train a random
features model on noise-free trajectories from the true system (an ODE). The
problems we study provably have a compact global attractor and are provably
(L63) or empirically (L96MS) ergodic; the invariant distribution is supported
on the global attractor and captures the statistics of long-time trajectories
which, by ergodicity, are independent of initial condition. The data
trajectories are generated using scipy’s implementation of Runge-Kutta 5(4)
(via solve_ivp) with absolute and relative tolerances both $10^{-9}$ and
maximum step size $10^{-4}$ [Dormand and Prince, 1980, Virtanen et al., 2020].
In order to obtain statistical results, we create 5 training trajectories from
the true system of interest with initial conditions sampled independently from
its attractor. Note that each training trajectory is long enough to explore
the attractor, and each is used to train a separate model; the purpose is to
observe the variance in learnt models with respect to randomly sampled paths
through the attractor. We use the same sampling procedure to generate short
independent validation and testing trajectories—we use 7 validation
trajectories and 10 testing trajectories (these are short because we only use
them to evaluate a model’s short term forecast performance; when assessing
long-term statistics of a learnt model, we compare to very long simulations
from the true system). All plots use error bars to represent empirical
estimates of the mean and standard deviation of the presented performance
metric, as computed by ensembling the performance of the 5 models (one per
training trajectory) over the 10 testing trajectories for a total of 70 random
performance evaluations.

Each training procedure also involves an independent draw of the random
feature functions as defined in (8.9). A validation step is subsequently
performed to optimize the hyperparameters $\omega,\beta$, as well as the
regularization parameter $\lambda$. We automate this validation using Bayesian
Optimization [Mockus, 1989, Nogueira, 2014], and find that it typically
identifies good hyperparameters within 30 iterations. The entire process of
entraining a model to a single, long training trajectory (including
hyperparameter validation) typically takes approximately $30$ minutes on a
single core of a 2.1GHz Skylake CPU with an allocated 1GB RAM. Given a
realization of random features and an optimal $\lambda$, we obtain the
minimizer $C^{\ast}$ using the Moore-Penrose Pseudoinverse implemented in
scipy (pinv2). This learned $C^{\ast}$, paired with its random feature
realization, is then used to predict 10 unseen testing trajectories (it is
given the true initial condition for each of these testing trajectories).

When implementing in continuous-time, given high frequency but discrete-time
data, two computational issues must be addressed: (i) extrapolation of the
data to continuous-time; (ii) discretization of the resulting integrals. The
approach we adopt avoids “inverse crimes” in which favourable behaviour is
observed because of agreement between the data generation mechanism (with a
specific integrator) and the approximation of the objective functions [Colton
and Kress, 2013, Kaipio and Somersalo, 2005, Wirgin, 2004]; see Queiruga et
al. [2020] for further illustration of this issue and Keller and Du [2021], Du
et al. [2021] for a rigorous analysis of this inversion process in the context
of linear multistep integration methods for deep learning. We interpolate the
data with a spline, to obtain continuous-time trajectories, and then
discretize the integrals using a simple Riemann sum; this strikes a desirable
balance between robustness and efficiency and avoids inverse crimes. The
discrete-time approaches, however, are able to learn not only model-
discrepancy, but also integrator-based discrepancies; hence, the discrete-time
methods may artificially appear to outperform continuous-time approaches,
when, in fact, their performances might simply be considered to be comparable.

####  6.2.2. Evaluation Criteria

Models are evaluated against the test set for their ability to predict
individual trajectories, as well as invariant statistics (the invariant
measure and the auto-correlation function).

Trajectory Validity Time: Given threshold $\gamma>0$, we find the first time
$t_{\gamma}$ at which the norm of discrepancy between true and approximate
solutions reaches $\gamma$:

$t_{\gamma}={\rm argmin}_{t\in[0,T]}\bigg{\\{}t\colon\
\|x(t)-x_{m}(t)\|\geq\gamma\overline{\|x(t)\|}\bigg{\\}},$

where $x(t)$ is the true solution to (2.2), $x_{m}(t)$ is the learned
approximation, and the normed time average $\overline{\|x(t)\|}$ is
approximated from training data. If the threshold is not violated on $[0,T]$,
we define $t_{\gamma}:=T$; this is rare in practice. We take $\gamma=0.05$
(i.e. $5\%$ relative divergence).

Invariant Distribution: To quantify errors in our reconstruction of the
invariant measure, we consider the Kullback-Leibler (KL) divergence [Kullback
and Leibler, 1951] between the true invariant measure $\mu$ and the invariant
measure produced by our learned model $\mu_{m}$. We approximate the divergence

$d_{\text{KL}}(\mu,\mu_{m}):=\int_{\mathbbm{R}}\log\bigg{(}\frac{d\mu}{d\mu_{m}}\bigg{)}d\mu$

by integrating kernel density estimates with respect to the Lebesgue measure.

Autocorrelation: We compare the autocorrelation function (ACF) with respect to
the invariant distribution of the true and learned models. We approximate the
ACF using a fast-fourier-transform for convolutions Seabold and Perktold
[2010], and compare them via a normalized $L_{2}$ norm of their difference.

####  6.2.3. Lorenz ’63 (L63)

##### Setting

The L63 system [Lorenz, 1963] is described by the following ODE

(6.1) $\displaystyle\dot{u}_{x}$ $\displaystyle=a(u_{y}-u_{x})$
$\displaystyle\dot{u}_{y}$ $\displaystyle=bu_{x}-u_{y}-u_{x}u_{z}$
$\displaystyle\dot{u}_{z}$ $\displaystyle=-cu_{z}+u_{x}u_{y}$

whose solutions are known to exhibit chaotic behavior for parameters $a=10,\
b=28,\ c=\frac{8}{3}$. We align these equations with our framework, starting
from equation (2.1), by letting $x=(u_{x},u_{y},u_{z})^{T}$ and defining
$f^{\dagger}(x)$ to be the vector field appearing on the right-hand-side in
(6.1). We define a discrete solution operator $\Psi^{\dagger}$ by numerical
integration of $f^{\dagger}$ over a fixed time window $\Delta t$ corresponding
to a uniform data sampling rate, so that the true system is given by (2.1) in
continuous-time and (2.6a) in discrete-time.

To simulate scenarios in which our available physics are good, but imperfect,
we assume there exists additive unknown model error of form

(6.2) $m^{\dagger}(x)=\epsilon\ m_{1}(x)$

with function $m_{1}$ determining the structure of model error, and scalar
coefficient $\epsilon$ determining its magnitude. Recall that
$f^{\dagger}=f_{0}+m^{\dagger}$ and we assume $f_{0}$ is known to us. Our task
is then to learn $f^{\dagger}$ by learning $m^{\dagger}$ and adding it to
$f_{0}$. The discrete solution operator $\Psi_{0}$ is obtained as in (2.6b) by
numerical integration of $f_{0}$ over a fixed time window $\Delta t$.

To simplify exposition, we explicitly define $m^{\dagger}$, then let
$f_{0}:=f^{\dagger}-m^{\dagger}$. We first consider the setting where

(6.3) $m_{1}(x):=\ \begin{bmatrix}0\\\ bu_{x}\\\ 0\end{bmatrix}$

(as in [Pathak et al., 2018b]) and modulate $\epsilon$ in (6.2) to control the
magnitude of the error term. In this case, $f_{0}$ can be viewed as the L63
equations with perturbed parameter $\tilde{b}=b(1-\epsilon)$, where $b$ is
artificially decreased by $100\epsilon\%$.

Then, we consider a more general case of heterogeneous, multi-dimensional
residual error by drawing $m_{1}$ from a zero-mean Gaussian Process (GP) with
a radial basis kernel (lengthscale 10). We form a map from $\mathbbm{R}^{3}$
into itself by constructing three independent draws from a scalar-valued GP on
$\mathbbm{R}^{3}.$ The resulting function is visualized in two-dimensional
projections in Figure 1.

Figure 1.  Here we visualize an example of the function $m_{1}$ in (6.2),
which is obtained as a single random draw from a zero-mean Gaussian Process
mapping $\mathbbm{R}^{3}\to\mathbbm{R}^{3}$. We have plotted its output
surface as three scalar functions (left to right) of the first two inputs (the
plot axes) with the third input component fixed at 0.

Observe that in the continuous-time framing, changes to $\epsilon$ do not
impact the complexity of the learned error term; however, it does grow the
magnitude of the error term. In the discrete-time framing, larger values of
$\epsilon$ can magnify the complexity of the discrepancy
$\Psi_{0}(x)-\Psi^{\dagger}(x)$.

##### Results

We perform a series of experiments with the L63 system in order to illustrate
key points about using data to learn model errors in dynamical systems. First,
we demonstrate that hybrid modeling tends to outperform data-only and physics-
only methods in terms of prediction quality. We first consider model error as
in (6.3); see Figure 2(a) in which we study performance (validity time) versus
model error amplitude ($\epsilon$), using random feature maps with $D=200$,
and a single trajectory of length $T=100$ sampled at timestep $\Delta
t=0.001$. Unless otherwise specified, this is also the configuration used in
subsequent experiments.

We see identical trends in Figure 2(b) for a more general case with the non-
parametric model error term constructed from Gaussian processes.
Interestingly, we see that for small and moderate amounts of model error
$\epsilon$, the hybrid methods substantially outperform data-only and physics-
only methods. Eventually, for large enough model discrepancy, the hybrid-
methods and data-only methods have similar performance; indeed the hybrid-
method may be outperformed by the data-only method at large discrepancies. For
the simple model error this appears to occur when the discrepancy term is
larger in magnitude than $f_{0}$ (e.g. for $b=28$ and $\epsilon=2$, the model
error term $\epsilon bu_{x}$ can take on values larger than $f^{\dagger}$
itself).

Figure 2(b) also shows that a continuous-time approach is favored over
discrete-time when using data-only methods, but suggests the converse in the
hybrid modeling context. We suspect this is an artifact of the different
integration schemes used in data generation, training, and testing phases; the
data are generated with a higher-fidelity integrator than the one available in
training and testing. For the continuous-time method, this presents a
fundamental limitation to the forecast quality (we chose this to avoid having
artificially high forecast validity times). However, the discrete-time method
can overcome this by not only learning the mechanistic model discrepancy, but
also the discrepancy term associated with a mis-matched integrator. This
typically happens when a closure is perfectly learnable and deterministic
(i.e. our Lorenz ’63 example); in this case, the combination of physics-based
and integrator-sourced closures can be learned nearly perfectly. In later
experiments with a multiscale system, the closures are considered approximate
(they model the mean of a noisy underlying process) and the discrete- and
continuous-time methods perform more similarly, because the inevitable
imperfections of the learned closure term dominate the error rather than the
mis-specified integrator. Note that approximate closures driven by scale-
separation are much more realistic; thus we should not expect the hybrid
discrete-time method to dramatically outperform hybrid continuous-time methods
unless other limitations are present (e.g. slow sampling rate).

Importantly, the parameter regime for which hybrid methods sustain advantage
over the imperfect physics-only method is substantial; the latter has
trajectory predictive performance which drops off rapidly for very small
$\epsilon$. This suggests that an apparently uninformative model can be
efficiently modified, by machine learning techniques, to obtain a useful model
that outperforms a _de novo_ learning approach.

(a)

(b)

Figure 2. These plots shows the temporal length of the forecast validity of
our learnt models of L63 (6.1), each as a function of model error, as
parameterized by $\epsilon$ (6.2) (with $D=200$, $T=100$, and $\Delta
t=0.001$). Continuous-time methods are shown in blue, discrete-time approaches
in orange. Dotted lines indicate purely data-driven methods to learn the
entire vector field defining the dynamics; solid lines indicate methods that
learn perturbations to the imperfect mechanistic models $f_{0}$ or $\Psi_{0}$.
Integration using the imperfect mechanistic model, without recourse to data,
is shown in green. In Figure 2(a), we employ the linear form of model error
$m_{1}$ defined in (6.3). In Figure 2(b), we let $m_{1}$ be a single draw from
a Gaussian Process, whose structure is shown in Figure 1. Here, we plot means,
with error bars as 1 standard deviation.

Next, we show that hybrid methods simplify the machine learning task in terms
of complexity of the learned function and, consequently, the amount of data
needed for the learning. Figure 3 examines prediction performance (validity
time) as a function of training data quantity using random feature maps with
$D=2000$ and a fixed parametric model error ($\epsilon=0.2$ in (6.2)) and
sampling rate $\Delta t=0.01$. We see that the hybrid methods substantially
outperform the data-only approaches in regimes with limited training data. For
the continuous-time example, we see an expected trend, where the data-only
methods are able to catch up to the hybrid methods with the acquisition of
more data. The discrete-time models do not exhibit this behavior, but we
expect the data-only discrete-time model to eventually catch up, albeit with
additional training data and number of parameters. Note that greater
expressivity is also required from data-only methods—our choice of a large
$D=2000$ aims to give all methods ample expressivity, and thus test
convergence with respect to training data quantity alone. These results
demonstrate that the advantage of hybrid modeling is magnified when training
data are limited and cannot fully inform _de novo_ learning. Figure 4 further
studies the impact of expressivity by again fixing a parametric model error
($\epsilon=0.05$ in (6.2)), training length $T=100$, and sampling rate $\Delta
t=0.001$. We see that all methods improve with a larger number of random
features, but that relative superiority of hybrid methods is maintained even
for $D=10000$.

Figure 3. Here we examine the performance of the proposed methods as a
function of the length of the interval over which the training data is
provided, where $\Delta t=0.01$, ($\epsilon=0.2$ in (6.2)), and $D=2000$ are
held constant for the L63 example (6.1). See description of Figure 2 for
explanation of legend. We observe that all methods improve with increasing
training lengths. We see that, in continuous-time, the primary benefit in
hybrid modeling is when the training data are limited. Figure 4. Here we
examine the performance of the proposed methods as a function of model
complexity, where $\Delta t=0.001$, $\epsilon=0.05$, and $T=100$ are held
constant for the L63 example (6.1). See description of Figure 2 for
explanation of legend. We observe that all methods improve with increasing
number of parameters, and that hybrid methods are especially beneficial when
available complexity is limited.

Finally, we study trade-offs between learning in discrete- versus continuous-
time for the L63 example (6.1). Figure 5 examines prediction performance
(validity time) as a function of data sampling rate $\Delta t$ using random
feature maps with $D=200$ with a fixed parametric model error ($\epsilon=0.05$
in (6.2)) and an abundance of training data $T=1000$. We observe that for fast
sampling rates ($\Delta t<0.01$), the continuous-time and discrete-time hybrid
methods have similar performance. For $\Delta t>0.01$, derivatives become
difficult to estimate from the data and the performance of the continuous-time
methods rapidly decline. However, the discrete-time methods sustain their
predictive performance for slower sampling rates ($\Delta t\in(0.01,0.1)$). At
some point, the discrete-time methods deteriorate as well, as the discrete map
becomes complex to learn at longer terms because of the sensitivity to initial
conditions that is a hallmark of chaotic systems. Here, the discrete-time
methods begin to fail around $\Delta t=0.2$; note that they can be extended to
longer time intervals by increasing $D$ and amount of training data, but
returns diminish quickly.

Figure 5. This shows temporal forecast validity as a function of the step size
of training data for the tested methods in the L63 example (6.1). We hold
fixed $D=200$, $\epsilon=0.05$, and $T=1000$. See description of Figure 2 for
explanation of legend. We see that while purely data-driven discrete-time
methods struggle at short time steps, the hybrid version thrives in this
scenario. All approaches, of course, eventually decay as large time steps
create more complex forward maps, due to sensitivity to initial conditions. We
also see continuous-time methods work well for small time steps, then
deteriorate in tandem with quality of estimated derivatives.

####  6.2.4. Lorenz ’96 Multiscale (L96MS) System

##### Setting

Here, we consider the multiscale system [Lorenz, 1996] of form (2.7), where
each variable $X_{k}\in\mathbbm{R}$ is coupled to a subgroup of fast variables
$Y_{k}\in\mathbbm{R}^{J}$. We have $X\in\mathbbm{R}^{K}$ and
$Y\in\mathbbm{R}^{K\times J}$. For $k=1\dots K$ and $j=1\dots J$, we write

(6.4a) $\displaystyle\dot{X}_{k}$
$\displaystyle=f_{k}(X)+h_{x}\overline{Y}_{k}$ (6.4b)
$\displaystyle\dot{Y}_{k,j}$
$\displaystyle=\frac{1}{\varepsilon}r_{j}(X_{k},Y_{k})$ (6.4c)
$\displaystyle\overline{Y}_{k}$
$\displaystyle=\frac{1}{J}\sum_{j=1}^{J}Y_{k,j}$ (6.4d) $\displaystyle
f_{k}(X)$ $\displaystyle=-X_{k-1}(X_{k-2}-X_{k+1})-X_{k}+F$ (6.4e)
$\displaystyle r_{j}(X_{k},Y_{k})$
$\displaystyle=-Y_{k,j+1}(Y_{k,j+2}-Y_{k,j-1})-Y_{k,j}+h_{y}X_{k}$ (6.4f)
$\displaystyle X_{k+K}$ $\displaystyle=X_{k},\quad Y_{k+K,j}=Y_{k,j},\quad
Y_{k,j+J}=Y_{k+1,j}$

where $\varepsilon>0$ is a scale-separation parameter,
$h_{x},h_{y}\in\mathbbm{R}$ govern the couplings between the fast and slow
systems, and $F>0$ provides a constant forcing. We set $K=9,\ J=8,\
h_{x}=-0.8,\ h_{y}=1,\ F=10$; this leads to chaotic dynamics for $\varepsilon$
small. When studying scale-separation, we consider $\varepsilon\in\\{2^{-7},\
2^{-5},\ 2^{-3},\ 2^{-1}\\}$.

We consider the setting in which we learn Markovian random features models in
variable $X$ alone, from $X$ data generated by the coupled $(X,Y)$ system.
Large scale-separation between the observed ($X$) and unobserved ($Y$) spaces
can simplify the problem of accounting for the unobserved components; in
particular, for sufficient scale-separation, we expect a Markovian term to
recover a large majority of the residual errors. In fact, we further simplify
this problem by learning a scalar-valued model error $M$ that is applied to
each $X_{k}$ identically in the slow system:

$\dot{X}_{k}\approx f_{k}(X)+M(X_{k}).$

This choice stems from observations about statistical interchangeability
amongst the slow variables of the system; these properties of the L96MS model
in the scale-separated regime are discussed in [Fatkullin and Vanden-Eijnden,
2004]. We can directly align our reduction of (6.4) with the Markovian hybrid
learning framework in (2.2) as follows:

$\displaystyle\dot{X}$ $\displaystyle\approx f_{0}(X)+m(X)$ $\displaystyle
f_{0}(X)$ $\displaystyle:=[f_{1}(X),\ \cdots\ ,f_{K}(X)]^{T}$ $\displaystyle
m(X)$ $\displaystyle:=[M(X_{1}),\ \cdots\ ,M(X_{K})]^{T}.$

##### Results

We plot the performance gains of our hybrid learning approaches in Figure 6 by
considering validity times of trajectory forecasts, estimation of the
invariant measure, and ACF estimation. In all three metrics (and for all
scale-separations $\varepsilon$), _de novo_ learning in discrete
($\Psi^{\dagger}\approx m$) and continuous-time ($f^{\dagger}\approx m$) is
inferior to using the nominal mechanistic model $f_{0}$. We found that the
amount of data used in these experiments is insufficient to learn the full
system from scratch. On the other hand, hybrid models in discrete
($\Psi^{\dagger}\approx\Psi_{0}+m$) and continuous-time ($f^{\dagger}\approx
f_{0}+m$) noticeably outperformed the nominal physics.

Surprisingly, Figure 6 shows that the Markovian closure methods still
qualitatively reproduce the invariant statistics even for large $\varepsilon$
settings where we would expect substantial memory effects. Figure 6 also
demonstrates this quantitatively using KL-divergence between invariant
measures and mean-squared-error between ACFs. It seems that for this
dissipative system, memory effects tend to average out in the invariant
statistics. However, the improvements in validity time for trajectory-based
forecasting deteriorate for $\varepsilon=2^{-1}$.

To visualize this non-Markovian structure, and how it might be exploited, we
examine the residuals from $f_{0}$ in Figure 7 and observe that there are
discernible trajectories walking around the Markovian closure term. For small
$\varepsilon$, these trajectories oscillate rapidly around the closure term.
For large $\varepsilon$ (e.g. $2^{-1}$), however, we observe a slow-moving
residual trajectory around the Markovian closure term. This indicates the
presence of a stronger memory component, and thus would benefit from a non-
Markovian approach to closure modeling.

Jiang and Harlim [2020] show that the memory component in this setting with
$\varepsilon=2^{-1}$ can be described using a closure term with a simple delay
embedding of the previous state at lag $0.01$. They learn the closure using a
kernel method cast in an RKHS framework, for which random feature methods
provide an approximation.

Figure 6. This figure shows the performance of different approaches to
modeling the L96MS slow subsystem (6.4). In $f^{\dagger}\approx f_{0}$, we
only use the nominal physics $f_{0}$. In $\Psi^{\dagger}\approx m$ and
$f^{\dagger}\approx m$, we try to learn the entire right-hand-side using only
data (in discrete- and continuous-time settings, respectively). In
$\Psi^{\dagger}\approx\Psi_{0}+m$ and $f^{\dagger}\approx f_{0}+m$, we focus
on learning Markovian residuals for the known physics (in discrete- and
continuous-time settings, respectively). The residual-based correctors
substantially outperform the nominal physics and purely data-driven methods
according to all presented metrics: invariant measure (shown qualitatively in
the first row and quantitatively in the third row), ACF (shown qualitatively
in the second row and quantitatively in the fourth row), and trajectory
forecasts (shown in the final row). The boxplots show the distributions of
quantitative metrics (e.g. KL-divergence, squared errors, validity time),
which come from different models, each trained on a different trajectory, and
generated using an independent random feature set. Notably, the Markovian
residual-based methods’ performance deteriorates for small scale-separation
($\varepsilon=2^{-1}$), where the Markovian assumption breaks down. Figure 7.
This figure shows the observed and estimated residuals of the nominal physics
model $f_{0}$ for the L96MS slow subsystem (6.4) at different scale-separation
factors. The first row shows the density of these residuals (yellow is high
density, blue is low), as well as the fit of our closure terms in continuous-
(blue) and discrete- (orange) time (the discrete model was normalized by
dividing by $\Delta t$). The second row shows temporal structure in the errors
of our residual fit by superimposing a short (T=1) one-dimensional trajectory
(this represents $\sim 0.1\%$ of training data).

###  6.3. Learning From Partial, Noisy Observations

In this section, we focus on the non-Markovian setting outlined in Section
2.3, and attempt to model the dynamics of the observable using (4.3), with
$f_{1}$, $f_{2}$ given by two-layer, fully connected neural networks with GeLU
activations [Hendrycks and Gimpel, 2016], and perform the learning by
minimizing (4.7) from Section 4.3.3, using 3DVAR for the data assimilation Law
et al. [2015, 2013], with the ADAM optimizer [Kingma and Ba, 2015]. The
learning rate was initialized at $0.01$ and tuned automatically using a
scheduler that halved the learning rate if the training error had not
decreased over 10 (mini-batched) epochs. Data were sampled at $\Delta t=0.01$
in all cases, and normalized to have mean zero and unit variance. Numerical
integration was performed with the torchdiffeq implementation of the Dormand-
Prince adaptive fifth-order Runge-Kutta method: for the L63 example, simple
backpropagated autodifferentiation was performed through this solver; for the
L96MS example, we used the adjoint method provided by [Rubanova et al., 2019].

####  6.3.1. Lorenz ’63

We first consider modeling the dynamics of the first-component of the L63
system in (6.1), where we noisily observe the first-component – that is, we
observe a noisy trajectory of $u_{x}$ (i.i.d. additive zero-mean, variance-one
Gaussian), but do not observe the remaining components $u_{y},u_{z}$. We
jointly trained on 100 trajectories, each of length $T=10$ and randomly
initialized from a box around the attractor; we chose this approach to ensure
that we had data coverage both on and off the attractor although we note that
similar success is obtained with a single trajectory of length $T=1000$.). The
neural network had width $50.$ We chose an assimilation time of $\tau_{1}=3$
and a forecast time of $\tau_{2}=0.1$. The optimization ran for approximately
200 epochs, and took roughly 24hrs on a single GPU. Adequate results were
obtained using a fixed 3DVAR gain matrix $K=[0.5,0,0]^{T}$. However, we
present results using the algorithm in which $K=\theta_{\textrm{DA}}$ is
jointly learned along with parameters $\theta_{\textrm{DYN}}$, as described in
Section 4.3.3; this demonstrates that the gain need not be known _a priori_.

First, we present results using knowledge that the correct hidden dimension
$d_{r}=2$: in Figure 8(a), we show an example of the trained model being
assimilated (using 3DVAR with learnt $K$) during the first 3 time units, then
predicting for another 7 time units; recall that training was performed using
only a $\tau_{2}=0.1$ forecasting horizon, but we evaluate on a longer time
horizon to provide a robust test metric. Observe that the learnt hidden
dynamics in gray are synchronized with the data, then used to perform the
forecast. In Figures 8(b) and 8(c), we show that by solving the system for
long time (here, $T=10^{4}$), we are able to accurately reproduce invariant
statistics (invariant measure and autocorrelation, resp.) for the true system.
In Figure 8(d), we show the evolution of the learnt $K$.

Next, we let $d_{r}=10$, exceeding the true dimension of the hidden states;
thus we are able to explore issues caused by learning an overly expressive (in
terms of dimension of hidden states) dynamical model. Figure 9 shows dynamics
for a learnt model in this setting; we found its reproduction of invariant
statistics to be similar to the cases in Figures 8(b) and 8(c), but omit the
plots for brevity. This success aligns with the approximation theory, as
discussed in Remark 5.5, and provides empirical reassurance that the
methodology can behave well in situations where the dimension of the hidden
variable is unknown and dimension $d_{r}$ used in learning exceeds its true
dimension. Nevertheless, we construct an example in Section 6.4 in which a
specific embedding of the true dynamics in a system of higher dimension can
lead to poor approximation; this is caused by an instability in the model
which allows departure from the invariant manifold on which the true dynamics
is accurately captured. However, we emphasize that this phenomenon is not
observed empirically in the experiment reported here with $d_{r}=10.$
Nonetheless we also note expected decreases in efficiency caused by over-
estimating the dimension of the hidden variable, during both model training
and testing; thus determining the smallest choice for $d_{r}$, compatible with
good approximation, is important. Recent research has addressed this challenge
in the discrete-time setting by applying manifold learning to a delay-
embedding space, then using the learnt manifold to inform initialization and
dimensionality of LSTM hidden states [Kemeth et al., 2021].

Note that our early attempts at achieving these numerical results, using the
optimization ideas in Sections 4.3.1 and 4.3.2, yielded unstable models that