# Mori-Zwanzig latent space Koopman closure for nonlinear autoencoder

Priyam Gupta
Department of Aeronautics
Imperial College London
London SW7 2AZ, United Kingdom
<EMAIL_ADDRESS>
&Peter J. Schmid
Department of Mechanical Engineering
KAUST
23955 Thuwal, Saudi Arabia
&Denis Sipp
DAAA
Onera
92190 Meudon, France
&Taraneh Sayadi
Institut Jean le Rond d’Alembert, Sorbonne University
75005, Paris, France
Institute for Combustion Technology, Aachen University
52062 Aachen, Germany &Georgios Rigas
Department of Aeronautics
Imperial College London
London SW7 2AZ, United Kingdom

###### Abstract

The Koopman operator presents an attractive approach to achieve global
linearization of nonlinear systems, making it a valuable method for
simplifying the understanding of complex dynamics. While data-driven
methodologies have exhibited promise in approximating finite Koopman
operators, they grapple with various challenges, such as the judicious
selection of observables, dimensionality reduction, and the ability to predict
complex system behaviours accurately. This study presents a novel approach
termed Mori-Zwanzig autoencoder (MZ-AE) to robustly approximate the Koopman
operator in low-dimensional spaces. The proposed method leverages a nonlinear
autoencoder to extract key observables for approximating a finite invariant
Koopman subspace and integrates a non-Markovian correction mechanism using the
Mori-Zwanzig formalism. Consequently, this approach yields a closed
representation of dynamics within the latent manifold of the nonlinear
autoencoder, thereby enhancing the precision and stability of the Koopman
operator approximation. Demonstrations showcase the technique’s ability to
capture regime transitions in the flow around a circular cylinder. It also
provided a low dimensional approximation for chaotic Kuramoto-Sivashinsky with
promising short-term predictability and robust long-term statistical
performance. By bridging the gap between data-driven techniques and the
mathematical foundations of Koopman theory, MZ-AE offers a promising avenue
for improved understanding and prediction of complex nonlinear dynamics.

_K_ eywords Nonlinear Model Reduction $\cdot$ Autoencoders $\cdot$ Koopman
operator $\cdot$ Mori-Zwanzig formalism

##  1 Introduction

Nonlinear systems are ubiquitous, spanning from the unsteady fluid flows and
the evolution of epidemics to intricate neural interactions in the brain.
These systems frequently exhibit a level of dimensionality that is exceedingly
high for any practical computational analysis. As a pragmatic solution,
reduced order models (ROM) are sought which provide tractable and accurate
dynamics for a small set of quantities of interest (resolved observables).
Although the intricate dynamics exhibited by these systems arise from the
nonlinear spatio-temporal interactions between multiple scales, it often makes
it difficult to find their analytical solutions. The Koopman operator [1, 2]
offers an alternative view through the lens of “dynamics of observables” which
facilitates a global linearisation for these inherently nonlinear systems.
This linear characteristic holds significant appeal across a spectrum of
applications, including nonlinear system identification [3] and nonlinear
control [4]. However, the Koopman operator’s definition within an infinite-
dimensional Hilbert space contradicts the fundamental goal of constructing a
ROM. Therefore, we seek a finite low-dimensional representation which entails
identifying a specific set of observables that span a finite invariant Koopman
subspace. This task is often challenging and requires approximation [5]. When
approximating, one must anticipate the accumulation of errors in the generated
state-space trajectories over time. The Mori-Zwanzig formalism [6, 7] enables
a closed-form representation of the infinite-dimensional Koopman operator by
providing a non-Markovian correction to its finite approximation. This
effectively closes the gap between the practical implementation of the Koopman
operator and its idealized form.

Dynamic mode decomposition (DMD) has shown a remarkable performance for
approximating the Koopman operator for fluid flows [8, 9, 10], or many other
dynamical systems, for example, epidemic evolution [11]. Nonetheless, it
struggles with the issue of choosing a proper set of observables. This is
mainly because deriving the closed-form solution in the Koopman theory
involves carefully resolving the key observables enriched in the nonlinear
information of the system and embed them into a linear dynamics. One approach
to identify these observables is to search for them within a pre-defined
dictionary (referred to as a feature map) of linear and nonlinear functions of
the state variables, as proposed by Williams  _et al._ [12]. This Extended DMD
(EDMD) method is equivalent to using a high-order Taylor series expansion
around equilibrium points as compared to merely a linear expansion by standard
DMD [13]. However, this approach relies on a priori knowledge of the behaviour
of the dynamical system. Another caveat of the same approach is the limited
representational capacity of the dictionary which can lead to overfitting due
to insufficient data.

An alternative technique for approximating the Koopman operator is to exploit
the universal approximation capability of the neural networks to learn the
observables from data. This method is generally referred to as the Deep
Koopman approach [14, 15, 16, 17] where, the state variables are passed
through a nonlinear autoencoder to produce a small set of observables enriched
with the nonlinearities of the dynamical system. To ensure the observables lie
in the Koopman invariant subspace, an approximate linear operator is learned
on the nonlinear manifold (latent space) of the autoencoder. The motivating
idea behind this approach is a two-step identification where (i) the
autoencoder learns energetically dominant modes, and (ii) the Koopman matrix
learns dynamically important features. Otto _et al._ [14] used a Linear
Recurrent Neural Network framework where the error of the learned Koopman
operator is minimized over multiple timesteps. Lusch _et al._ [15] extended
this work to dynamical systems with continuous frequency spectra. They
obtained the parametric dependence of the Koopman operator on the continuously
varying frequency using an auxiliary network. Pan _et al._ [18] proposed a
probabilistic Koopman learning framework based on Bayesian neural networks for
continuous dynamical systems while offering a stability constraint on their
Koopman parameterization.

Data-driven Koopman learning methods are founded on the assumption that a non-
trivial finite-dimensional Koopman invariant subspace exists [19]. Even if
this assumption holds true, it has proven to be exceedingly challenging to
resolve this finite set of observables that completely closes the dynamics
[5]. In order to obtain a closed dynamics, we need to account for the effects
of the unresolved observables that complete the invariant Koopman subspace.
Mori and Zwanzig introduced a general framework for the closed equations of
the resolved observables. They demonstrated that the interactions between
resolved and unresolved observables manifest themselves as non-Markovian non-
local effects onto the resolved observables. To accommodate these
interactions, it decomposes the dynamics into three parts – a Markovian term,
a non-Markovian or memory term, and a noise term – which together form a so-
called Generalised Langevin Equation (GLE). In this decomposition, the memory
and the noise terms are responsible for the effects of the unresolved
observables. While the evolution equations obtained for resolved observables
are methodically exact, it does not provide reduced computational complexity
without approximations. This is primarily because deriving the analytical form
of the memory kernel which accounts for the non-Markovian effect is an arduous
task. However, the GLE provides an excellent starting point to model closure
terms in a non-Markovian form.

It has been shown that a higher-order correction to the approximate Koopman
operator can be obtained using the Mori-Zwanzig formalism by accounting for
the residual dynamics through the non-Markovian term. Lin _et al._ [20]
proposed a data-driven method for this purpose that recursively learns the
memory kernels. Curtis _et al._ [21] used the popular optimal prediction
framework[22] to provide higher-order correction terms for DMD. This was
further extended in [23] where $t$-model [24] was utilized for memory
approximation. The primary challenge for these methods is the judicious choice
of their observables. They are selected from a predefined dictionary of
functions that has the same shortcomings as seen before, such as overfitting
and the need of a priori knowledge of the system as in Extended DMD.

This work proposes a data-driven framework termed Mori-Zwanzig autoencoder
(MZ-AE), that exploits the Mori-Zwanzig formalism and approximates the
invariant Koopman subspace in the latent manifold of a nonlinear autoencoder.
A higher-order non-Markovian correction is provided to the approximate Koopman
operator which guides it back to the true trajectory upon deviation. Through
this approach, we tackle the following challenges:

  * •

Choice of observables: A nonlinear autoencoder provides a data-driven approach
to finding the key observables for best approximation of the Koopman operator.
This relaxes the requirement of a priori knowledge of the system while
providing a suitable coordinate transformation into an approximate Koopman
invariant subspace.

  * •

Low-dimensional approximation: There is a need to have a tractable number of
observables that evolve accurately with time on a low-dimensional subspace.
Excessive number of observables can lead to spurious eigenvalues which are
non-physical and impede interpretability[19]. A nonlinear autoencoder with the
help of the nonlinear correction term allows an aggressively low-dimensional
model.

  * •

Predictability: It is highly challenging to obtain finite linear
approximations for a chaotic system [25]. The latter fail to achieve good
long-term predictions and require large number of observables to unfold the
dynamics further [14] defeating the purpose of a ROM. Through the Mori-Zwanzig
formalism we provide a memory correction term to ensure maximum predictability
and stability.

Organization. §2 presents the relevant background theory, briefly discussing
model-order reduction, Koopman formalism, and Mori-Zwanzig decomposition. §3
provides the details of the proposed MZ-AE algorithm. §4 sets up the numerical
experiments and discusses the results obtained. Finally §5, provides
conclusions and a discussion on future work.

##  2 Background theory

This section lays down the mathematical background building upto the framework
proposed in this work. We start with Poincare’s state space view of the
dynamical system in order to set up the operator theoretic method of Koopman.
Subsequently, we motivate our approach from the reduced order model
perspective since we seek a low-dimensional representation of the infinite
Koopman operator. Finally, we discuss the Mori-Zwanzig decomposition of the
observable dynamics and provide a rationale for the proposed MZ-AE framework
to close the finite Koopman approximation for nonlinear autoencoders.

###  2.1 Data-driven model order reduction

We consider an autonomous dynamical system evolving an $N$ dimensional vector
$\boldsymbol{\Phi}$ on the state space $\mathcal{M}\subset\mathbb{R}^{N}$,

$\frac{d\boldsymbol{\Phi}(t)}{dt}=S(\boldsymbol{\Phi}(t)),$ (1)

where $\boldsymbol{\Phi}$ is a set of physical space variables or state
variables which characterises the state of the system and S is an evolution
operator. Specifically, we are interested in a discrete autonomous dynamical
system sampled at fixed time step $\Delta t$,

$\boldsymbol{\Phi}_{n+1}=\mathbf{S}_{t}(\boldsymbol{\Phi}_{n}),$ (2)

where $\mathbf{\Phi}_{n}=\mathbf{\Phi}(n\Delta t)$. We refer the evolution
operator S in equation 2 as full order model (FOM) which is typically
characterized by nonlinearity, while the state variables represented by
$\mathbf{\Phi}$ can exhibit high dimensionality. These properties make it
prohibitively expensive to obtain the solution for equation 2 and only limited
insights can be gained into the inherent system dynamics. To address this, we
seek an accurate and amenable ROM for a tractable number of relevant
quantities $\mathbf{u}(\mathbf{\Phi})\in\mathbb{R}^{r}$ such that $r\ll N$.

A data-driven ROM extracts this low-dimensional representation of the full-
order model from the data snapshots. This is done by projecting onto a data-
driven expansion basis as in Proper Orthogonal Decomposition (POD) [26]. In
this method the basis functions (or modes) are ranked according to their
energy content and the dimension reduction entails decomposing the state
variable into first $r$ dominant energy modes,

${\Phi}(\mathbf{x},t)=\sum_{j=0}^{r}a_{j}(t){u}_{j}(\mathbf{x})+\mathbf{\epsilon},$
(3)

where $\mathbf{x}$ represents the spatial coordinates, $u_{j}$ is the spatial
mode, $a_{j}$ represents the time-dependent amplitude of the spatial modes and
$\mathbf{\epsilon}$ represent the residual error from unresolved modes. The
dynamical system can then be projected on to these modes using Galerkin
projection to obtain a reduced order model [27]. Alternatively, a data- driven
model, such as a Recurrent Neural Network (RNN), can be used to learn the
projected non-linear dynamics [28]. There is also a need to model effects of
unresolved modes on the dynamics of resolved modes and provide a closure[29,
30]. The key limitation of these methods is that the POD modes tend to
intertwine spatial and temporal frequencies which complicates their physical
relevance [31]. As a solution to this limitation, Dynamical mode decomposition
based methods provide a data-driven approach to extract dynamically relevant
modes representing spatio-temporal coherent structures by regressing a linear
evolution operator motivated by Koopman formulation.

###  2.2 Koopman formalism

The Koopman operator is an infinite-dimensional mapping that accomplishes
global linearization of a nonlinear dynamical system by transforming it into a
linear system of observables. Let us define a Hilbert space $\mathcal{H}$ of
the complex-valued functions on the smooth manifold $\mathcal{M}$. An
observable is any linear or nonlinear measurement function of the state
variables $\mathbf{g}(\mathbf{\Phi})$ that lie in $\mathcal{H}$. Starting from
the initial condition $\mathbf{\Phi_{0}},$ the observables at any time $t$
attain the value $\mathbf{g}(\mathbf{\Phi}(\mathbf{\Phi}_{0},t))$ which we
abbreviate by $\mathbf{g}(\mathbf{\Phi_{0}},t)$. The nonlinear dynamical
system $(\ref{eq:Discrete_Dynamical_System})$ can be represented linearly
through the infinite-dimensional Koopman operator $\mathcal{K}$ such that

$\mathbf{g}(\mathbf{\Phi}_{0},t)=\mathcal{K}^{t}\mathbf{g}(\mathbf{\Phi}_{0}).$
(4)

For the discrete system (2), the Koopman operator takes the form

$\mathbf{g}_{n+1}=\mathcal{K}_{\Delta t}\mathbf{g}_{n},$ (5)

where $\mathbf{g}_{n}=\mathbf{g}(\mathbf{\Phi}_{n})$. Such a linear operation
is possible only if the observables span an invariant Koopman subspace.
However, an infinite-dimensional subspace is computationally intractable. We
seek a finite-dimensional subspace which is spanned by a set of $r$
observables
$\mathbf{\hat{g}}=\\{\hat{g}_{i}\\}_{i=1}^{r}\in\mathcal{\hat{H}}\subset\mathcal{H}$
such that their evolution always lie within this subspace
$(\mathcal{K}\mathbf{\hat{g}}\in\mathcal{\hat{H}})$. In practice, it is
challenging to find such a finite closed invariant subspace, and we have to
work with an approximation of the invariant subspace. This results in the
accumulation of errors in the trajectory over time due to loss of dynamical
information leading to a closure problem [5]. Mori-Zwanzig provides a
framework to obtain a closed-form equation for finite approximation of the
invariant Koopman subspace by decomposing the infinite Hilbert space into a
finite resolved subspace $(\hat{\mathbf{g}}\in\mathcal{\hat{H}})$ and an
infinite unresolved subspace
$(\tilde{\mathbf{g}}\in\mathcal{\tilde{H}}\subset\mathcal{H})$.

###  2.3 Mori-Zwanzig decomposition

We begin by obtaining a linear representation for the nonlinear system (1)
through the linear partial differential equation of the observables
$\mathbf{g}(\mathbf{\Phi_{0}},t)$. This is achieved by using the Liouville
operator ($\mathcal{L}$),

$\frac{\partial}{\partial
t}\mathbf{g}(\mathbf{\Phi_{0}},t)=\mathcal{L}\mathbf{g}(\mathbf{\Phi_{0}},t),$
(6)

where
$\mathcal{L}=\sum_{i=1}^{N}S_{i}(\mathbf{\Phi_{0}})\frac{\partial}{\partial\mathbf{\Phi_{0}}_{i}}$.
This expression for the Liouville operator can be obtained by taking the chain
derivative of the observables through time. The solution of this Liouville
equation generates a Koopman operator following

$\mathbf{g}(\mathbf{\Phi_{0}},t)=e^{t\mathcal{L}}\mathbf{g}(\mathbf{\Phi}_{0})=\mathcal{K}^{t}\mathbf{g}(\mathbf{\Phi_{0}}),$
(7)

where $\mathcal{K}^{t}=e^{t\mathcal{L}}$.

For a reduced-order model, we seek the evolution of only a few quantities of
interest (limited to $r$ with $r\ll N$); these are the resolved observables
$\hat{\mathbf{g}}$ which lie in the image of a projection operator
$\mathcal{P}$ such that $\mathcal{P}\mathbf{g}=\hat{\mathbf{g}}$. Naturally,
there exists a complementary projection $\mathcal{Q}=\mathcal{I}-\mathcal{P}$
with the set of unresolved observables ($\tilde{\mathbf{g}}$) as its image.
Here, $\mathcal{I}$ represents the identity operator. An appropriate
projection operator needs to be chosen such as Mori’s projection operator [7]
which comes with a specific commutative property and will be used in obtaining
the projected dynamics of resolved observables. The mathematical details of
these projection operators are discussed in [32]. The Liouville dynamics
(described by equation (6)) is then decomposed into a resolved and an
unresolved subspace such that

$\frac{\partial}{\partial
t}e^{t\mathcal{L}}\mathbf{g}(\mathbf{\Phi}_{0})=e^{t\mathcal{L}}(\mathcal{P}+\mathcal{Q})\mathcal{L}\mathbf{g}(\mathbf{\Phi}_{0})=e^{t\mathcal{L}}\mathcal{P}\mathcal{L}\mathbf{g}(\mathbf{\Phi}_{0})+e^{t\mathcal{L}}\mathcal{Q}\mathcal{L}\mathbf{g}(\mathbf{\Phi}_{0}).$
(8)

The second term on the right-hand side of equation (8) represents the
unresolved dynamics that need to be processed to obtain the memory and noise
terms. This is achieved by applying Dyson’s identity (see §C),

$e^{t\mathcal{L}}=\int_{0}^{t}e^{(t-s)\mathcal{L}}\mathcal{P}\mathcal{L}e^{s\mathcal{Q}\mathcal{L}}ds+e^{t\mathcal{Q}\mathcal{L}},$
(9)

on the unresolved dynamics
$\mathcal{Q}\mathcal{L}\mathbf{g}(\mathbf{\Phi}_{0})$ within equation (8).
Then, using equation 7, we obtain the governing equation which represents the
exact evolution of the observable $\mathbf{g}(\mathbf{\Phi}_{0},t)$,

$\frac{\partial}{\partial
t}\mathbf{g}(\mathbf{\Phi}_{0},t)=\underbrace{\mathcal{P}\mathcal{L}\mathbf{g}(\mathbf{\Phi}_{0},t)}_{\text{Markov}}+\underbrace{\int_{0}^{t}\mathcal{P}\mathcal{L}e^{s\mathcal{Q}\mathcal{L}}\mathcal{Q}\mathcal{L}\mathbf{g}(\mathbf{\Phi}_{0},t-s)ds}_{\text{Memory}}+\underbrace{e^{t\mathcal{Q}\mathcal{L}}\mathcal{Q}\mathcal{L}\mathbf{g}(\mathbf{\Phi}_{0})}_{\text{Noise}}.$
(10)

This expression is a formal rewriting of equation (6). It decomposes the high-
dimensional dynamics into a resolved part (Markov term), an unresolved part
(noise term), and the interaction between them (memory term). The first term
which accounts for the resolved dynamics relies solely on their instantaneous
values, hence it is referred to as the Markov term
($\mathbf{M}=\mathcal{PL}$). The noise term, which is
$\mathbf{F}(t)=\mathcal{Q}\mathcal{L}e^{t\mathcal{Q}\mathcal{L}}\mathbf{g}(\mathbf{\Phi}_{0})$,
incorporates the dynamics of the observables in the unresolved subspace. It is
apparent that the memory term corresponds to a convolution in time of
observables in the past $\mathbf{g}(t-s)$ through a Kernel
$\mathbf{\Omega}(s)=\mathcal{P}\mathcal{L}e^{s\mathcal{Q}\mathcal{L}}\mathcal{Q}\mathcal{L}$,
in which we recognize part of the noise term $\mathbf{F}(s)$. Consequently, it
accounts for the interaction of the resolved and unresolved observables. Eq.
(10) also shows that, given a linear projection operator, the convolution term
is in principle linear with respect to the past measurements
$\mathbf{g}(\mathbf{\Phi}_{0},t-s)$ and that the noise term is orthogonal to
$\mathcal{P}$ (since $\mathcal{P}\mathcal{Q}=0$). Using the above notations,
we rewrite equation 10 in a simplified form as

$\frac{\partial}{\partial
t}\mathbf{g}(\mathbf{\Phi}_{0},t)=\mathbf{M}\mathbf{g}(\mathbf{\Phi}_{0},t)+\int_{0}^{t}\mathbf{\Omega}(s)\mathbf{g}(t-s)ds+\mathbf{F}(t).$
(11)

Since we are interested in a closed subspace of only resolved observables, we
obtain the dynamics projected onto the resolved observables. To this end we
apply the projection operator on equation (11) and use their property to
commute through the evolution operator [32] yielding

$\frac{\partial}{\partial
t}\mathbf{\hat{g}}(\mathbf{\Phi}_{0},t)=\mathbf{M}\mathbf{\hat{g}}(\mathbf{\Phi}_{0},t)+\int_{0}^{t}\mathbf{\Omega}(s)\mathbf{\hat{g}}(t-s)ds.$
(12)

This equation represents the true evolution of the resolved observables. The
discrete-in-time evolution equation of the resolved observables reads

$\mathbf{\hat{g}}_{n+1}=\mathcal{\tilde{K}}_{\Delta
t}\mathbf{\hat{g}}_{n}+\sum_{k=1}^{n-1}\mathbf{\Omega}_{k}\mathbf{\hat{g}}_{n-k},$
(13)

where $\tilde{\mathcal{K}}_{\Delta t}$ is a discrete finite Koopman operator.
When $\mathbf{\hat{g}}$ resides precisely within a finite invariant Koopman
subspace, there is no unresolved subspace, and consequently, the memory term
becomes null. This notion of finite invariant Koopman subspace is often untrue
and we have to work with a finite approximation of the infinite-dimensional
Koopman operator. The Markov term, $\tilde{\mathcal{K}}_{\Delta t}$, in
equation (13) is analogous to the finite approximation of the Koopman dynamics
while the memory term, $\mathbf{\Omega}_{k}$, provides closure by accounting
for the effects of unresolved observables.

##  3 Mori-Zwanzig formulation for nonlinear autoencoder (MZ-AE)

The proposed MZ-AE framework is based on a nonlinear autoencoder for
identifying the observables $\mathbf{\hat{g}}$ for an optimal approximation of
a linearly invariant finite Koopman subspace. Simultaneously, the dynamics in
this approximate invariant subspace, described by $\mathcal{\tilde{K}}_{\Delta
t}$, are closed by accounting for the effect of unresolved observables through
the memory term $\mathbf{\Omega}_{k}$ (see equation 13), which is
parameterised using a long short-term memory (LSTM) network. In this section,
we lay out the architecture of the different components used in MZ-AE and
formulate the objective loss function.

###  3.1 Nonlinear Autoencoder

A nonlinear autoencoder (AE) provides a nonlinear coordinate transformation to
a low dimensional latent space. It primarily consists of an encoder and a
decoder. The encoder encodes a given input space
$\boldsymbol{\Phi}\in\mathbb{R}^{N\times N_{t}}$ with spatial and temporal
dimension of $N$ and $N_{t}$ respectively, into the latent space
$\mathbf{g}\in\mathbb{R}^{r\times N_{t}}$, where $r$ is the latent space size
and $r\ll N$. The decoder remaps this latent space back to the input space by
minimizing an appropriate norm (such as $L_{2}$), giving a reconstruction
$\boldsymbol{{\Phi^{\prime}}}$. This process of encoding and reconstructing
using the encoder $\mathcal{E}$ and decoder $\mathcal{D}$ can be formulated as

$\mathbf{g}=\mathcal{E}\left(\boldsymbol{\Phi};\mathbf{\theta}_{\mathcal{E}}\right),\quad\boldsymbol{\Phi}\approx\boldsymbol{{\Phi}^{\prime}}=\mathcal{D}\left(\mathbf{g};\mathbf{\theta}_{\mathcal{D}}\right).$
(14)

In this study, the encoder and decoder architectures are multilayer perceptron
(MLP, see §A) networks, parameterised with weights and biases
$\theta_{\mathcal{E}}$ and $\theta_{\mathcal{D}}$, respectively.

###  3.2 MZ-AE framework

Figure 1: Schematic for the MZ-AE framework.

The architecture of MZ-AE is shown in Figure 1. We consider a set of $m$
discrete snapshots of the resolved observables
$\hat{\mathbf{g}}:\mathbb{R}^{N}\rightarrow\mathbb{R}^{r}$ which are nonlinear
functions of $N$ dimensional state variables $\mathbf{\Phi}$. We find these
resolved observables in the latent manifold of a nonlinear autoencoder,
$\hat{\mathbf{g}}_{n}=\mathcal{E}(\mathbf{\Phi}_{n})$, sampled at time $\Delta
t$. We obtain the dynamics of these resolved observables using the discrete
projected Mori-Zwanzig decomposition (equation 13). Here, we learn the
approximate Koopman operator ($\tilde{\mathcal{K}}$) and the memory term using
data-driven methods. The approximate Koopman matrix is obtained as the linear
operator in the least squares sense over the resolved observables
$\mathbf{\hat{g}}$. For robustness, we express the memory convolution as a
nonlinear function of the past resolved observables. An LSTM (see §B) is used
to regress the residual of the linear operator over time as an estimate of the
memory term along with any numerical evolution error. The latent space
dynamics are obtained as

$\hat{\mathbf{g}}_{n+1}=\tilde{\mathcal{K}}_{\Delta
t}\hat{\mathbf{g}}_{n}+\Omega(\hat{\mathbf{g}}_{n-1},...,\hat{\mathbf{g}}_{n-q}).$
(15)

The LSTM model recurs over past $q$ number of observables to produce the
residual ($\mathbf{\xi}_{n}$) for the next time step,

$\mathbf{\xi}_{n}=\Omega(\hat{\mathbf{g}}_{n-1},...,\hat{\mathbf{g}}_{n-q};\theta_{\Omega}).$
(16)

The evolved resolved observables can be reconstructed back to original state
space using the non-linear decoder $(\mathcal{D})$. We primarily use three
loss functions to enforce these dynamics into the neural network architecture:

  * •

Autoencoder reconstruction error: Minimizing this error allows the autoencoder
to learn a low-dimensional manifold over which an approximate Koopman subspace
can be enforced. It is the squared $L_{2}$-norm of difference between the
reconstructed and the actual state variables, stated as

$\mathcal{J}_{rec}=\frac{1}{m}\sum_{n=0}^{m}\left\|\boldsymbol{\Phi}_{n}-\mathcal{E}{\mathcal{D}(\boldsymbol{\Phi}_{n})}\right\|_{2}^{2}.$
(17)

  * •

Linear evolution error: This error enforces the Koopman invariant subspace in
the low-dimensional manifold of the autoencoder. Starting from a snapshot
$\hat{\mathbf{g}}_{n}$, a trajectory over $\eta$ timesteps is produced using
just the linear operator $\\{\tilde{\mathcal{K}}_{\Delta
t}\hat{\mathbf{g}}_{n},...,\tilde{\mathcal{K}}^{\eta}_{\Delta
t}\hat{\mathbf{g}}_{n}\\}$. The residual w.r.t. the true trajectory
$\\{\hat{\mathbf{g}}_{n},\hat{\mathbf{g}}_{n+1},...,\hat{\mathbf{g}}_{n+\eta}\\}$
in the latent manifold is then collected using

$\mathbf{R}_{n}=\hat{\mathbf{g}}_{n+1}-\tilde{\mathcal{K}}_{\Delta
t}\hat{\mathbf{g}}_{n}.$ (18)

The linear evolution error is obtained as the mean squared error of the
residual $\mathbf{R}$ over time,

$\mathcal{J}_{R}=\frac{1}{\eta(m-\eta)}\sum_{n=0}^{m-\eta}\sum_{i=0}^{\eta-1}|\mathbf{R}_{n+i}|^{2}.$
(19)

  * •

Residual error: This error regresses the memory model output,
$\mathbf{\xi}_{n}$ (from equation 16), with the linear evolution residual
($\mathbf{R}_{n}$),

$\mathcal{J}_{\mathbf{\Omega}}=\frac{1}{\eta(m-\eta)}\sum_{n=0}^{m-\eta}\sum_{i=0}^{\eta-1}|\mathbf{R}_{n+i}-\mathbf{\xi}_{n+i}|^{2},$
(20)

where
$\mathbf{\xi}_{n+i}=\Omega(\hat{\mathbf{g}}_{n+i-1},...,\hat{\mathbf{g}}_{n+i-q})$.

The overall objective function is a weighted sum of these loss functions,

$\mathcal{J}=\alpha_{rec}\mathcal{J}_{rec}+\alpha_{R}\mathcal{J}_{R}+\alpha_{\Omega}\mathcal{J}_{\Omega}.$
(21)

where the weights $\alpha_{rec}$, $\alpha_{R}$ and $\alpha_{\Omega}$ are tuned
for optimal prediction. The learnable weights corresponding to this parametric
model include the encoder ($\theta_{\mathcal{E}}$), decoder
($\theta_{\mathcal{D}}$), memory model ($\theta_{\Omega}$) and the MZ-AE
linear operator ($\theta_{\mathcal{\tilde{K}}}$). These weights are learnt by
backpropagating through the objective loss function $\mathcal{J}$,

$\mathbf{\theta}^{*}=\arg\min_{\mathbf{\theta}}\mathcal{J}(\mathbf{\Phi};\mathbf{\theta}_{\mathcal{E}},\mathbf{\theta}_{\mathcal{D}},\mathbf{\theta}_{\Omega},\mathbf{\theta}_{\mathcal{\tilde{K}}}).$
(22)

The parameters are then updated in the optimal direction using the ADAM[33]
optimization algorithm. The training strategy proposed in this work allows the
MZ-AE linear operator to learn the dynamics to its maximum potential while
constraining the memory model to the residual dynamics. Further, minimizing
the residual of the linear operator over multiple timesteps allows it to
identify low energy modes which may be insignificant for short-term dynamics.
However, they may amplify several steps ahead in the future effecting the
long-term dynamics [14].

###  3.3 Model Selection

The proposed model has several hyper-parameters which need to be carefully
chosen for optimal performance. Primarily grid search method was utilized for
determining the neural network architecture of the encoder, decoder, and LSTM.
The first key parameter is the number of resolved observables $r$. In
practice, it is decided empirically by evaluating the prediction error of the
model for increasing the number of observables until it converges. In this
study, we searched for the number of observables close to the inherent
dimensions of the model. The prediction horizon ($\eta$) is the number of
timesteps over which the approximate Koopman operator is allowed to predict
before the memory model learns the obtained residual. It was also chosen using
the grid search method.

For the LSTM, the two important hyper-parameters are the window length ($q$)
and the number of hidden units ($N_{hu}$). The window length ($q$) of the LSTM
is defined as the number of past observables that the LSTM model recurs over
to predict the residual of the approximate Koopman operator. Different
dynamical systems have different memory lengths, therefore window length of
the LSTM needs to be carefully chosen. A small value of $q$ could lead to loss
of information while a large value could increase the computational complexity
leading to difficulty in convergence while training. The number of hidden
units is the size of the hidden state and cell state vectors in LSTM that are
responsible for short-term and long-term information, respectively. A grid
search is done over different values to converge onto a suitable pair of $q$
and $N_{hu}$.

##  4 Numerical Experiments

The proposed reduced-order modelling framework presented in §3 is demonstrated
on the two-dimensional flow over a cylinder, which is a popular benchmark for
reduced order models. Subsequently, we evaluate the performance of the model
on a chaotic Kuramoto-Sivashinsky system. The model parameters for each test
case are tabulated in §D. To assess the impact of the memory model, we also
examine the performance of MZ-AE in comparison to an approximate Koopman
operator learned without the memory model.

###  4.1 Two Dimensional Cylinder Flow

(a)

(b)

Figure 2: Results obtained with two observables for the transition dynamics of
cylinder flow, (a) Comparison of Kinetic Energy prediction by the MZ-AE and
the linear operator without memory model. (b) $L_{2}$ Norm (Scaled by max
value) of the observable dynamics contribution by MZ-AE linear operator
($\mathcal{\tilde{K}}_{\Delta t}\hat{\mathbf{g}}$) and memory model
($\mathbf{\xi}$).

In order to probe the workings of the proposed model-reduction architecture,
we first consider the Karman vortex street in the wake of the circular
cylinder. It is a suitable first test case since it represents an
intrinsically low-dimensional system [34] expressed in high-dimensional
snapshots. The trajectory transitions from an unstable equilibrium to a limit
cycle, enabling us to analyze the behavior of the system between these two
states. Through our proposed model, we aim to compress the dynamics to a
minimum number of Koopman observables while ensuring satisfactory predictive
accuracy over multiple timesteps. The training data is obtained using a high-
fidelity solver – Nektar++ [35] based on a spectral/hp element framework for
Reynolds number $Re=100$. The baseflow for the simulations was obtained using
the adaptive selective frequency damping (ASFD) method [36] in Nektar++. The
flow was simulated for a cylinder of unit diameter ($D$) in a domain of length
$-5D$ to $15D$ in the streamwise direction and $-5D$ to $5D$ in the normal
direction. The velocity snapshots were created by uniformly sampling the
streamwise and normal components of the velocity from a smaller subdomain with
a streamwise and normal length of $-2.5D$ to $10D$ and $-2.5D$ to $2.5D,$
respectively. The velocity components were then concatenated into a state
vector with $7738$ components. A total of $1600$ such velocity snapshots were
collected at time intervals of $0.125D/U_{\infty}$, where $U_{\infty}$ is the
free stream velocity of the flow. For training, $800$ snapshots are sampled at
a time interval of $0.25D/U_{\infty}$, and the remaining snapshots were used
for testing the model. This allowed us to test the prediction capability of
the model over the complete trajectory from the unstable equilibrium to the
attractive limit cycle.

In the supercritical regime, vortex shedding occurs at the limit cycle which
is inherently two-dimensional. A third dimension is often required to capture
the transient dynamics [34]. In this study, we used the nonlinear autoencoder
to extract only two modes (i.e. $r=2$) from the velocity snapshots and expect
the memory term to capture the transient dynamics in the residual. The model
was allowed to learn the trajectory from the unstable equilibrium to the
stable limit cycle. The prediction horizon $\eta$ of 10 timesteps and window
size of $q=8$ timesteps was found to be optimal for this case. The MZ-AE model
is able to forecast the dynamics through the transition regime till the
attracting limit cycle and way beyond the training data horizon, as shown in
Figure 2(a). It is observed that the memory correction term allows the MZ-AE
linear operator to transition. As expected, the linear operator without
correction is unable to transition since the two modes are insufficient to
learn the transition dynamics. This suggests that the MZ-AE provides an
aggressive reduction in the dimension compared to other data-driven Koopman
approaches [15, 14] that require a higher number of observables to capture
transition.

(a) MZ-AE linear operator phase portrait

(b) MZ-AE memory model phase portrait

(c) Time evolution of MZ-AE linear operator contribution

(d) Time evolution of memory model contribution ($\mathbf{\xi}$)

Figure 3: Time Evolution of the contribution by MZ-AE linear operator
($\mathcal{\tilde{K}}_{\Delta t}\hat{\mathbf{g}}$) and memory model
($\mathbf{\xi}$) to the observable dynamics.

The $L_{2}$ norm of the contributions of the MZ-AE linear operator and the
memory model to the observable dynamics are reported in Figure 2(b). It shows
that the memory model is primarily active during the transition regime while
it provides minor corrections in the limit cycle. This fact is further
consolidated in Figure 3 which shows the phase portrait of the MZ-AE linear
operator and the memory correction. Moreover, the magnitude of the
contribution to observables by the MZ-AE linear operator is much higher than
that by the memory model as can be seen in Figure 3. It is evident that the
linear operator captures accurately the limit-cycle dynamics while the memory
model enables the accurate transition from the unstable equlibrium.

Figure 4: Eigenvalues of the MZ-AE linear operator for cylinder flow.

Finally, we compare the eigenspectrum of the MZ-AE linear operator and the
linear operator obtained without memory correction in Figure (4). The
continuous time eigenvalues ($\mu$) are obtained from discrete-time
eigenvalues ($\lambda$) using the relation $\mu=log(\lambda)/\Delta t$. Both
the models are successfully able to capture the dominant vortex shedding
frequency of the limit cycle at $Re=100$ where the Strouhal number $St=\mu
D/2\pi U_{\infty}=0.178$. We observe that the MZ-AE linear operator has a
negative growth rate which again suggests that the instability required for
the transition is supplied by the memory term. However, the Koopman operator
without memory correction has a positive growth rate, which can be attributed
to the linear model’s attempt to learn the transition dynamics.

###  4.2 Kuramoto-Sivashinsky

The second case we study is the chaotic Kuramoto-Sivashinsky (KS) system which
has been developed as a model for flame front instabilities [37] and has also
been used to model weak fluid turbulence [38]. Again, it is a popular test
case for reduced-order models since obtaining DNS data is relatively
inexpensive, while the equation exhibits a rich and diverse dynamics at
relatively low values of the bifurcation parameter. We consider the one-
dimensional form of the system where the length of the domain acts as the
bifurcation parameter. The KS-system is governed by a fourth-order partial
differential equation,

$u_{t}+\lambda u_{xx}+u_{xxxx}+\frac{1}{2}(u_{x})^{2}=0,\quad 0\leq x\leq L,$
(23)

with periodic boundary conditions over the domain length $L$,

$u(0,t)=u(L,t),u_{x}(0,t)=u_{x}(L,t),$ (24)

and a user-specified initial condition,

$u(x,t=0)=u_{0}.$ (25)

The subscripts indicate temporal or spatial derivatives. In equation (23),
$u_{xx}$ and $u_{xxxx}$ correspond to a one-dimensional production and
dissipation of energy, while $u_{x}^{2}$ introduces nonlinearity to the
system. We choose a domain length of $L=22,$ which induces chaotic behavior
with a leading Lyapunov exponent of $\lambda_{1}=0.043$ and a Kaplan-Yorke
dimension $D_{KY}$ of $5.198$ [39]. The domain was spatially discretized using
$N=256$ Fourier modes, and the resulting system was integrated in time using a
fourth-order Runge-Kutta implicit-explicit time stepper [40] with a timestep
of $\Delta t=0.025$. We sample the generated solutions at $\Delta t=0.25$ and
eliminate initial $5000$ snapshots to exclude any transient data. A total of
$8\times 10^{4}$ snapshots were created which were divided into training,
validation, and testing datasets of size $5\times 10^{4}$, $1\times 10^{4}$,
and $2\times 10^{4},$ respectively.

For chaotic systems, the Kaplan-Yorke dimension ($D_{KY}$) gives a good
estimate of the dimension of the system’s attractor. We searched for the
number of resolved observables near $D_{KY}=5.198$ and found 8 observables to
be optimal. The prediction horizon ($\eta$) of 20 timesteps and an optimal
window length $q$ of $15$ timesteps (corresponding to $\approx 0.13$ LTUs) was
taken. We present the prediction of the model over $2.5$ LTUs of unseen test
data in Figure 5 which indicates that the model is able to capture dominant
characteristics and length scales of the system for more than 2 LTUs. In
particular, it is able to predict the interactions of the waves as can be seen
from the contour plots. Naturally, the chaotic nature of the flow leads to
higher absolute errors as it evolves beyond 1.25 LTU.

Figure 5: Comparison of (a) DNS solution and (b) MZ-AE solution with respect
to (c) the absolute error of the DNS and MZ-AE solution.

(a)

(b)

(c)

Figure 6: Comparison of statistics for Kuramoto-Sivashinsky system for MZ-AE
prediction, linear operator prediction without memory and DNS solution of (a)
kinetic energy spectrum, (b) kinetic energy probability distribution, (c)
cross-correlation function of solution at locations $x=0$ and $x=1$.

In order to investigate the long-term behavior of the model we examine the
statistics of the system over 70 LTUs of unseen test data which are shown in
Figure 6. The time-averaged kinetic energy spectrum of the predictions is
presented in Figure 6(a) which provides the distribution of energy in
different scales of the flow. We can observe that the MZ-AE is able to
accurately capture energy on large scales while the model without memory
corrections deviates from the actual spectrum. Further, the MZ-AE model is
able to capture the long-term statistics of the flow as can be seen from the
kinetic energy distribution in Figure 6(b) and the cross-correlation function
(CCF) of the solution at coordinates $x=0$ and $x=1$ in Figure 6(c). The model
without memory corrections significantly fails to capture the kinetic energy
distribution showing its inability to predict the long-term dynamics for
chaotic systems.

The $L_{2}$ norm of the contribution to the observables is presented in Figure
7(b) which shows that the MZ-AE linear operator is primarily responsible for
the dynamics while the memory model provides corrections to keep it on the
desired trajectory. We analyze the eigenvalues of the MZ-AE linear operator
and the linear operator without memory correction in Figure 8. It is observed
that the majority of the eigenvalues of the MZ-AE linear operator have a
negative growth rate and are marginally stable. At $L=22$, the system has
three pairs of unstable eigenvalues [41]. The MZ-AE linear operator appears to
learn the effects of the unstable modes through the marginally stable modes.
We also notice that the eigenvalues of the linear operator without memory are
more damped (larger negative $Re(\mu)$) than the MZ-AE linear operator.
Furthermore, the former has only two pairs of complex conjugate eigenvalues,
whereas the latter has three pairs. This observation suggests that the former
may require additional modes to adequately capture all the dynamics.
Additionally, the disparity in the spectrum of the two linear operators
indicates that memory correction has an impact on the nature of the learned
Koopman modes. Further analysis of this effect is necessary and will be a part
of our future work.

(a)

(b)

Figure 7: (a) Comparison of Kinetic Energy of trajectories of DNS solution and
MZ-AE prediction for Kuramoto-Sivashinsky System at L = 22 with 8 observables,
(b) corresponding $L_{2}$ Norm of the contribution to the observable dynamics
by MZ-AE linear operator ($\mathcal{\tilde{K}}_{\Delta t}\hat{\mathbf{g}}$)
and memory model ($\mathbf{\xi}$). Figure 8: Eigenvalues of the MZ-AE linear
operator for Kuramoto-Sivashinsky system.

##  5 Conclusion

This work proposed a theoretical and computational framework, motivated by the
Mori-Zwanzig formalism, to furnish a closed-form approximation of the Koopman
operator for a nonlinear autoencoder. The framework provides an interpretable
reduced-order model while ensuring accurate reconstruction and evolution of
the dynamics. We demonstrated that, by introducing the memory correction term,
the proposed technique is able to learn regime transitions in flow past a
circular cylinder with only two observables which capture the limit cycle
dynamics, while the memory correction term enables the accurate prediction of
the unstable equilibrium. Furthermore, the model is able to extract the
fundamental frequencies of the flow, thus demonstrating its capability to
learn the essential dynamical features while remaining interpretable. In
another application, we tested the model on the Kuramoto-Sivashinsky system in
the chaotic parameter regime where it showed promising short-term
predictability and good long-term statistical performance for considerably low
model orders. Further, the proposed training strategy constrains the memory
model to the residual dynamics. This ensures that the approximate Koopman
operator has maximum contribution to the evolution, giving us dynamically
relevant eigenmodes. The future work involves analyzing the complete spectrum
of the dynamics in the latent manifold of the autoencoder and the application
of the proposed framework to fully turbulent regimes.

## Acknowledgments

This work was supported by the Imperial College London - CNRS PhD Joint
Program.

## References

  * [1] Bernard O. Koopman and J. V. Neumann.  Dynamical systems of continuous spectra.  Proceedings of the National Academy of Sciences, 18(3):255–263, 1932.
  * [2] Igor Mezić.  Spectral properties of dynamical systems, model reduction and decompositions.  Nonlinear Dynamics, 41:309–325, 2005.
  * [3] Stefan Klus, Feliks Nüske, Sebastian Peitz, Jan-Hendrik Niemann, Cecilia Clementi, and Christof Schütte.  Data-driven approximation of the koopman generator: Model reduction, system identification, and control.  Physica D: Nonlinear Phenomena, 406:132416, 2020.
  * [4] Eurika Kaiser, J. Nathan Kutz, and Steven L. Brunton.  Data-driven discovery of Koopman eigenfunctions for control.  Machine Learning: Science and Technology, 2(3):035023, 2021.
  * [5] Steven L. Brunton, Bingni W. Brunton, Joshua L. Proctor, and J. Nathan Kutz.  Koopman invariant subspaces and finite linear representations of nonlinear dynamical systems for control.  PloS one, 11(2):e0150171, 2016.
  * [6] Hazime Mori.  Transport, collective motion, and brownian motion.  Progress of theoretical physics, 33(3):423–455, 1965.
  * [7] Robert Zwanzig.  Nonlinear generalized langevin equations.  Journal of Statistical Physics, 9(3):215–220, 1973.
  * [8] Jacob Page and Rich R. Kerswell.  Koopman analysis of burgers equation.  Physical Review Fluids, 3(7):071901, 2018.
  * [9] Taraneh Sayadi, Peter J. Schmid, Joseph W. Nichols, and Parviz Moin.  Reduced-order representation of near-wall structures in the late transitional boundary layer.  Journal of fluid mechanics, 748:278–301, 2014.
  * [10] Milan Korda, Mihai Putinar, and Igor Mezić.  Data-driven spectral analysis of the Koopman operator.  Applied and Computational Harmonic Analysis, 48(2):599–629, 2020.
  * [11] Joshua L. Proctor and Philip A. Eckhoff.  Discovering dynamic patterns from infectious disease data using dynamic mode decomposition.  International health, 7(2):139–145, 2015.
  * [12] Matthew O. Williams, Ioannis G. Kevrekidis, and Clarence W. Rowley.  A data–driven approximation of the koopman operator: Extending dynamic mode decomposition.  Journal of Nonlinear Science, 25:1307–1346, 2015.
  * [13] Peter J. Schmid.  Dynamic mode decomposition and its variants.  Annual Review of Fluid Mechanics, 54:225–254, 2022.
  * [14] Samuel E. Otto and Clarence W. Rowley.  Linearly recurrent autoencoder networks for learning dynamics.  SIAM Journal on Applied Dynamical Systems, 18(1):558–593, 2019.
  * [15] Bethany Lusch, J. Nathan Kutz, and Steven L. Brunton.  Deep learning for universal linear embeddings of nonlinear dynamics.  Nature communications, 9(1):1–10, 2018.
  * [16] Naoya Takeishi, Yoshinobu Kawahara, and Takehisa Yairi.  Learning koopman invariant subspaces for dynamic mode decomposition.  Advances in Neural Information Processing Systems, 30, 2017.
  * [17] Enoch Yeung, Soumya Kundu, and Nathan Hodas.  Learning deep neural network representations for Koopman operators of nonlinear dynamical systems.  In 2019 American Control Conference (ACC), pages 4832–4839. IEEE, 2019.
  * [18] Shaowu Pan and Karthik Duraisamy.  Physics-informed probabilistic learning of linear embeddings of nonlinear dynamics with guaranteed stability.  SIAM Journal on Applied Dynamical Systems, 19(1):480–509, 2020.
  * [19] Matthew J. Colbrook, Lorna J. Ayton, and Máté Szőke.  Residual dynamic mode decomposition: robust and verified koopmanism.  Journal of Fluid Mechanics, 955:A21, 2023.
  * [20] Yen Ting Lin, Yifeng Tian, Daniel Livescu, and Marian Anghel.  Data-driven learning for the Mori–Zwanzig formalism: A generalization of the koopman learning framework.  SIAM Journal on Applied Dynamical Systems, 20(4):2558–2601, 2021.
  * [21] Christopher W. Curtis and Daniel Jay Alford-Lago.  Dynamic-mode decomposition and optimal prediction.  Physical Review E, 103(1):012201, 2021.
  * [22] Alexandre J. Chorin, Ole H. Hald, and Raz Kupferman.  Optimal prediction and the Mori–Zwanzig representation of irreversible processes.  Proceedings of the National Academy of Sciences, 97(7):2968–2973, 2000.
  * [23] Aleksandr Katrutsa, Sergey Utyuzhnikov, and Ivan Oseledets.  Extension of dynamic mode decomposition for dynamic systems with incomplete information based on t-model of optimal prediction.  Journal of Computational Physics, 476:111913, 2023.
  * [24] Alexandre J. Chorin, Ole H. Hald, and Raz Kupferman.  Optimal prediction with memory.  Physica D: Nonlinear Phenomena, 166(3-4):239–257, 2002.
  * [25] Marko Budišić, Ryan Mohr, and Igor Mezić.  Applied Koopmanism.  Chaos: An Interdisciplinary Journal of Nonlinear Science, 22(4), 2012.
  * [26] JL Lumley, AM Yaglom, and VI Tatarski.  Atmospheric turbulence and radio wave propagation.  Journal of computational Chemistry, 23(13):1236–1243, 1967.
  * [27] Clarence W. Rowley and Scott TM Dawson.  Model reduction for flow analysis and control.  Annual Review of Fluid Mechanics, 49:387–417, 2017.
  * [28] Romit Maulik, Arvind Mohan, Bethany Lusch, Sandeep Madireddy, Prasanna Balaprakash, and Daniel Livescu.  Time-series learning of latent-space dynamics for reduced-order model closure.  Physica D: Nonlinear Phenomena, 405:132368, 2020.
  * [29] Qian Wang, Nicolò Ripamonti, and Jan S. Hesthaven.  Recurrent neural network closure of parametric POD-Galerkin reduced-order models based on the Mori-Zwanzig formalism.  Journal of Computational Physics, 410:109402, 2020.
  * [30] Emmanuel Menier, Michele Alessandro Bucci, Mouadh Yagoubi, Lionel Mathelin, and Marc Schoenauer.  CD-ROM: Complemented Deep-Reduced order model.  Computer Methods in Applied Mechanics and Engineering, 410:115985, 2023.
  * [31] Kilian Oberleithner, Moritz Sieber, Christian Navid Nayeri, Christian Oliver Paschereit, Christoph Petz, H-C Hege, Bernd R. Noack, and I Wygnanski.  Three-dimensional coherent structures in a swirling jet undergoing vortex breakdown: stability analysis and empirical mode construction.  Journal of fluid mechanics, 679:383–414, 2011.
  * [32] Jason M. Dominy and Daniele Venturi.  Duality and conditional expectations in the Nakajima-Mori-Zwanzig formulation.  Journal of Mathematical Physics, 58(8), 2017.
  * [33] Diederik P. Kingma and Jimmy Ba.  Adam: A method for stochastic optimization.  arXiv preprint arXiv:1412.6980, 2014.
  * [34] Bernd R. Noack, Konstantin Afanasiev, Marek Morzyński, Gilead Tadmor, and Frank Thiele.  A hierarchy of low-dimensional models for the transient and post-transient cylinder wake.  Journal of Fluid Mechanics, 497:335–363, 2003.
  * [35] Chris D. Cantwell, David Moxey, Andrew Comerford, Alessandro Bolis, Gabriele Rocco, Gianmarco Mengaldo, Daniele De Grazia, Sergey Yakovlev, J-E Lombard, Dirk Ekelschot, et al.  Nektar++: An open-source spectral/hp element framework.  Computer physics communications, 192:205–219, 2015.
  * [36] Bastien E. Jordi, Colin J. Cotter, and Spencer J. Sherwin.  An adaptive selective frequency damping method.  Physics of Fluids, 27(9), 2015.
  * [37] G. Sivashinsky.  Nonlinear analysis of hydrodynamic instability in laminar flames—i. derivation of basic equations.  In Dynamics of Curved Fronts, pages 459–488. Elsevier, 1988.
  * [38] James M Hyman, Basil Nicolaenko, and Stéphane Zaleski.  Order and complexity in the Kuramoto-Sivashinsky model of weakly turbulent interfaces.  Physica D: Nonlinear Phenomena, 23(1-3):265–292, 1986.
  * [39] Russell A. Edson, Judith E. Bunder, Trent W. Mattner, and Anthony J. Roberts.  Lyapunov exponents of the Kuramoto–Sivashinsky pde.  The ANZIAM Journal, 61(3):270–285, 2019.
  * [40] Benjamin Pachev, Jared P. Whitehead, and Shane A. McQuarrie.  Concurrent multiparameter learning demonstrated on the Kuramoto–Sivashinsky equation.  SIAM Journal on Scientific Computing, 44(5):A2974–A2990, 2022.
  * [41] Predrag Cvitanović, Ruslan L. Davidchack, and Evangelos Siminos.  On the state space geometry of the Kuramoto–Sivashinsky flow in a periodic domain.  SIAM Journal on Applied Dynamical Systems, 9(1):1–33, 2010.
  * [42] Sepp Hochreiter and Jürgen Schmidhuber.  Long short-term memory.  Neural computation, 9(8):1735–1780, 1997.
  * [43] Louis-S Bouchard.  Mori-Zwanzig Equations With Time-Dependent Liouvillian.  arXiv preprint arXiv:0709.1358, 2007.

##  Appendix A Multi Layer Perceptron

An Artificial Neural Network is a complex web of interconnected artificial
neurons that develops a nonlinear mapping between two spaces using a linear
combination of simple nonlinear functions. A multi-layer perceptron is the
most fundamental form of an artificial neural network. It consists of multiple
neurons intertwined together in a complex network. Neurons, being the basic
building blocks perform nonlinear input-output regression mapping. A single
neuron takes $n$ inputs $x_{1},x_{2},x_{3}…,x_{n}$, and its corresponding $n$
weights
$w\textsubscript{1},w\textsubscript{2},w\textsubscript{3}…,w\textsubscript{n}$.
The weighted output of a single neuron is expressed as

$z=\sum_{i=1}^{n}w\textsubscript{i}x\textsubscript{i}+b,$ (26)

where b is the bias term. This acts as the signal to the next neuron and to
introduce nonlinearity in the system, the output is obtained by multiplying
the weighted sum with a nonlinear activation function. The combined output
from a neuron is given by

$a(x)=\sigma\left(\sum_{i=1}^{n}w\textsubscript{i}x\textsubscript{i}+b)\right),$
(27)

where $a(x)$ is the activated output and $z$ is the signal to the neuron. The
neural network has at least three layers comprising input layers, hidden
layers and output layers. Each layer receives input from the last layer and
after activation, passes the output to the next layer. This process of
obtaining the activation and feeding them forward into the network is known as
a feed-forward network. MLP utilizes a supervised learning algorithm known as
backpropagation. Learning occurs after each piece of data is processed, and
weights along with bias are updated to minimize the losses.

##  Appendix B Long short-term memory network

Figure 9: Unfolded computational graph for a LSTM.

Recurrent Neural Networks (RNN) are a class of neural networks that are
recursively applied over a time series to learn the temporal patterns in the
dynamical system. An RNN model ($\Omega$) with learnable parameters
$\mathbf{\theta}$ takes the input $\mathbf{g}\in{\mathbb{R}^{r}}$ and
processes it into an internal or hidden state
$\mathbf{h}\in\mathbb{R}^{N_{hu}}$. These hidden states are then recurrently
evolved through time from $\mathbf{h}_{0}$ to $\mathbf{h}_{n}$ using equation

$\mathbf{h}_{n}=\mathbf{F}(\mathbf{h}_{n-1},\mathbf{g}_{n};\mathbf{\theta}_{\Omega_{1}}).$
(28)

The final hidden state $\mathbf{h}_{n}$ is then passed through another dense
layer to obtain the desired output $\mathbf{\xi}_{n}$

$\mathbf{\xi}_{n}=\mathbf{G}(\mathbf{h}_{n};\mathbf{\theta}_{\Omega_{2}}).$
(29)

Here, $\mathbf{\theta}_{\Omega_{1}}$ and $\mathbf{\theta}_{\Omega_{2}}$ are
weights of the entire RNN model $\mathbf{\theta}_{\Omega}$. In principle, a
RNN is supposed to have infinite memory retention but in practice, they often
suffer from long-term memory loss due to vanishing gradients. LSTM networks
introduced by Hochreiter _et al._ [42] are a type of gated RNN that are
designed to deal with these problems. They use a set of specialized gates and
memory cells to store and manipulate information over time. The memory cells
employed are; (i) a cell state, $\mathbf{C}_{n}$, that carries the long-term
correlations and (ii) a hidden state, $\mathbf{h}_{n}$, that is responsible
for short-term trends. Firstly, the input vector $\mathbf{g}$ is mapped to the
size of the hidden state using a dense layer. Subsequently, the gates filter
information through the following set of equations,

$\displaystyle\mathbf{i}_{n}$
$\displaystyle=\sigma(\mathbf{W}_{i}\cdot[\mathbf{h}_{n-1},\mathbf{g}_{n}]+\mathbf{b}_{i}),$
(30) $\displaystyle\mathbf{f}_{n}$
$\displaystyle=\sigma(\mathbf{W}_{f}\cdot[\mathbf{h}_{n-1},\mathbf{g}_{n}]+\mathbf{b}_{f}),$
$\displaystyle\mathbf{o}_{n}$
$\displaystyle=\sigma(\mathbf{W}_{o}\cdot[\mathbf{h}_{n-1},\mathbf{g}_{n}]+\mathbf{b}_{o}),$
$\displaystyle\mathbf{\tilde{C}}_{n}$
$\displaystyle=\tanh(\mathbf{W}_{c}\cdot[\mathbf{h}_{n-1},\mathbf{g}_{n}]+\mathbf{b}_{c}),$
$\displaystyle\mathbf{C}_{n}$
$\displaystyle=(1-\mathbf{f})_{n}\cdot\mathbf{C}_{n-1}+\mathbf{i}_{n}\cdot\mathbf{\tilde{C}}_{n},$
$\displaystyle\mathbf{h}_{n}$
$\displaystyle=\mathbf{o}_{n}\cdot\tanh(\mathbf{C}_{n}),$

where sigmoid ($\sigma$) and tanh are nonlinear functions. The three gates
that are used in an LSTM to control the data flow are: (i) forget gate
($\mathbf{f}_{n}$) which filters out information to be retained from the cell
state, (ii) input gate ($\mathbf{i}_{n}$) which decides what information
should be added to the cell state and (iii) output gate ($\mathbf{o}_{n}$)
decides the quantity of the current cell state that should be a part of the
current hidden state. The flow information through these gates is
schematically shown in Figure (9). The set of equations (30) are recurred over
until the desired time step ($n\Delta t$) is reached which produces the hidden
state $\mathbf{h}_{n-1}$.

##  Appendix C Dyson identity

For convenience, we provide the derivation of Dyson’s identity [43]. We begin
by considering the operator

$W(t)=e^{-\mathcal{L}t}e^{\mathcal{Q}\mathcal{L}t}.$ (31)

Differentiating in time gives

$\displaystyle\dot{W}(t)$
$\displaystyle=-e^{-\mathcal{L}t}\mathcal{L}e^{\mathcal{Q}\mathcal{L}t}+e^{-\mathcal{L}t}\mathcal{QL}e^{\mathcal{Q}\mathcal{L}t}$
(32)
$\displaystyle=-e^{-\mathcal{L}t}\mathcal{(}{I}-\mathcal{Q})\mathcal{L}e^{\mathcal{QL}t}$
$\displaystyle=-e^{-\mathcal{L}t}\mathcal{P}e^{\mathcal{QL}t}.$

Integrating equation 32 with respect to time, from $0$ to $t$, gives

$W(t)=W(0)-\int_{0}^{t}e^{-\mathcal{L}s}\mathcal{P}e^{\mathcal{QL}s}ds.$ (33)

We note from equation 31 that $W(t=0)=1$. Utilising this fact and substituting
equation 33 in equation 31, we obtain

$e^{-\mathcal{L}t}e^{\mathcal{Q}\mathcal{L}t}=1-\int_{0}^{t}e^{-\mathcal{L}s}\mathcal{P}e^{\mathcal{QL}s}ds.$
(34)

This gives us the Dyson-identity,

$e^{\mathcal{Q}\mathcal{L}t}=e^{\mathcal{L}t}-\int_{0}^{t}e^{(t-s)\mathcal{L}}\mathcal{P}e^{\mathcal{QL}s}ds.$
(35)

##  Appendix D Hyper-parameters for MZ-AE

###  D.1 Cylinder Flow

Table 1: Training Parameters for Cylinder Wake Hidden and Cell state size | 40
---|---
Batch Size | 16
Objective Function Weights ($\alpha_{rec}$, $\alpha_{R}$, $\alpha_{\Omega}$) | 10e+2, 1, 1
Epochs | 3000
Learning Rate (Adam Optimiser) | 10e-5
Prediction horizon ($\eta$) | 10
Table 2: Encoder and Decoder parameters for Cylinder Flow. Encoder | Decoder
---|---
Input Layer | 7738 x 512 | Bottleneck Layer | 8 x 64
Hidden Layer 1 | 512 x 256 | Hidden Layer 1 | 64 x 128
Hidden Layer 2 | 256 x 128 | Hidden Layer 2 | 128 x 256
Hidden Layer 3 | 128 x 64 | Hidden Layer 3 | 256 x 512
Bottleneck Layer | 64 x 8 | Output Layer | 512 x 7738
Activation Function $\sigma$ | Relu

###  D.2 Kuramoto-Sivashinsky

Table 3: Training Parameters for Kuramoto-Sivashinsky Hidden and Cell state size | 100
---|---
Batch Size | 128
Objective Function Weights ($\alpha_{rec}$, $\alpha_{R}$, $\alpha_{\Omega}$) | 10e+2, 1, 1
Epochs | 10000
Learning Rate (Adam Optimiser) | 10e-5
Prediction horizon ($\eta$) | 20
Table 4: Encoder and Decoder parameters for Kuramoto-Sivashinsky. Encoder | Decoder
---|---
Input Layer | 256 x 512 | Bottleneck Layer | 8 x 64
Hidden Layer 1 | 512 x 256 | Hidden Layer 1 | 64 x 128
Hidden Layer 2 | 256 x 128 | Hidden Layer 2 | 128 x 256
Hidden Layer 3 | 128 x 64 | Hidden Layer 3 | 256 x 512
Bottleneck Layer | 64 x 8 | Output Layer | 512 x 256
Activation Function $\sigma$ | Relu