# Enhancing Computational Efficiency in Multiscale Systems using Deep Learning
of Coordinates and Flow Maps

Asif Hamid  Department of Electrical Engineering, Islamic University of
Science and Technology, Awantipora,192122, Jammu & Kashmir, India  Danish
Rafiq  Shahkar Ahmad Nahvi  Mohammad Abid Bazaz  Department of Electrical
Engineering, National Institute of Technology, Srinagar, 190006, Jammu &
Kashmir, India

###### Abstract

Complex systems often show macroscopic coherent behavior due to the
interactions of microscopic agents like molecules, cells, or individuals in a
population with their environment. However, simulating such systems poses
several computational challenges during simulation as the underlying dynamics
vary and span wide spatiotemporal scales of interest. To capture the fast-
evolving features, finer time steps are required while ensuring that the
simulation time is long enough to capture the slow-scale behavior, making the
analyses computationally unmanageable. This paper showcases how deep learning
techniques can be used to develop a precise time-stepping approach for
multiscale systems using the joint discovery of coordinates and flow maps.
While the former allows us to represent the multiscale dynamics on a
representative basis, the latter enables the iterative time-stepping
estimation of the reduced variables. The resulting framework achieves state-
of-the-art predictive accuracy while incurring lesser computational costs. We
demonstrate this ability of the proposed scheme on the large-scale Fitzhugh
Nagumo neuron model and the 1D Kuramoto–Sivashinsky equation in the chaotic
regime.

###### keywords:

Multiscale modeling , Deep learning , Data-driven modeling , Large scale
systems , Scientific computing

##  1 Introduction

In the field of computational sciences and scientific computing, the numerical
solution of partial differential equations (PDEs) is essential for
understanding complex physical phenomena across a diverse range of
disciplines, from fluid dynamics and electromagnetism to structural mechanics
and quantum physics [1, 2]. However, complex PDEs often exhibit strong
nonlinear behavior that makes it difficult or impossible to obtain analytical
closed-form solutions. Therefore, numerical techniques are commonly employed
to obtain approximate solutions of PDEs. For instance, the finite element
analysis (FEA) [3] or the finite difference methods (FDM) [4] are used to
approximate the derivative term in the PDEs by discretizing the spatial
domains of interest into a grid or mesh. However, a fine discretization on
complicated boundaries results in a large system of ordinary differential
equations (ODEs) that can be computationally expensive to simulate.
Additionally, PDEs governing physical laws often exhibit a multiscale nature,
i.e., their behavior evolves on multiple distinct timescales that further
increase the computational costs during numerical solutions [5]. A common
situation is when two FEA models of different scales exist. In this case,
evaluating a single element or time step of the large-scale FEA model requires
at least one full simulation of the small-scale FEA model. This incurs
extensive runtimes and significant computational costs.

While the computational costs for such systems remain high, their online
evaluations are a major concern in many scientific and engineering
applications requiring many system evaluations within a many-query framework.
Some prominent applications include parameter and design optimization [6],
multi-objective optimization [7], uncertainty propagation [8], high-throughput
computing [9], and reinforcement learning (RL) [10, 11]. Even though expensive
computational simulations are essential to resolve all the relevant scales,
the important quantities of interest can be frequently characterized by a
coarse-grained or average behavior [12]. However, determining the appropriate
degree of freedom for such coarse-grained representations has always been
challenging. As a result, various research efforts have been put into
developing appropriate combinations of coarse-grained and fine-scale
simulations that can accelerate the intense simulations at a controlled
accuracy. Prominent hybrid methods include the equation-free framework (EFF)
[13, 14], the heterogeneous multiscale method (HMM) [15], and the flow
averaged integrator (FLAVOR) method [16]. These methods distinguish between a
fine (micro) scale, which is expensive to simulate, and a coarse (macro)
scale. While EFF and HMM use coarse time steppers to simulate the evolution of
macroscopic variables through a microscopic simulation, the FLAVOR method uses
averaging flows to study the behavior. However, the accuracy of these methods
is highly dependent on the separation of time scales, the time integrator
used, and how well the information is captured at each scale. In addition to
these methods, the classical Fourier and wavelet-based techniques can easily
separate the coarse-grain and fine-grain features in both space and time [17];
however, these methods focus on identifying exclusively temporal or
exclusively spatial coherence, which restricts their use for discovering
models and predicting system states directly from data [18]. Nevertheless,
these techniques are highly used and serve as the mathematical foundation of
multiresolution analysis (MRA) [19].

Recently, there has been growing interest in predicting the multiscale
behavior of systems directly from observational data using data-driven methods
based on machine learning (ML) [20]. In this direction, various neural
network-based schemes have been employed to forecast the solution of a
multiscale system. For instance, ML methods have been used for various
applications such as modeling dynamical systems [21, 22, 23, 24], in
mesoscopic material modeling [25], in biological systems [26, 27], and in
molecular kinematics, [28, 29]. The EFF framework has been extended with ML
algorithms that learn the time interrogators and the transfer operators in an
equation-free manner [30, 31]. In addition, spare regression algorithms have
been proposed to discover Poincare maps for multiscale problems [32] to
discover slow timescale dynamics from signal exhibiting multiscale timescales
[12]. Besides, the idea of learning a flow map of dynamical systems using
neural networks [33, 34, 35] has been extended to multiscale problems in [36],
wherein a novel multiscale hierarchical time-stepping (HiTS) method has been
proposed. In particular, the authors have suggested using a hierarchy of
neural network time stepper (NNTS) models that are trained at different time
steps and then coupled for multiscale purposes. However, as we will show in
this paper, the multiscale HiTS framework potentially suffers from huge
computational expenses for systems governed by PDEs.

Thus, the aim of the paper is two-fold: (i) to discover a coordinate system in
which the underlying dynamics of multiscale signals can be represented
accurately on a few representative bases, and (ii) to provide a time-stepping
scheme that marches the reduced coordinates forward in time. While we apply
deep autoencoders (AEs) to accomplish the former goal, we make use of the
multiscale HiTS estimation process to achieve the latter goal. In this way, we
extend the applicability of the multiscale HiTS method to PDEs. Thus, the
proposed method of latent hierarchical time-stepping (L-HiTS) learns an end-
to-end mapping that accurately yet efficiently predicts the multiscale
behavior of the PDEs from data. Although the multiscale HiTS scheme has been
demonstrated for PDEs [36], the need for real-time applications and tasks for
non-stationary dynamics calls for more efficient computations that the
proposed scheme intends to address. The result is that the L-HiTS method
maintains the same accuracy level as the multiscale HiTS method while
providing significant computational speedups. To summarize, the key
contributions of this paper are as follows:

  1. 1.

We present a novel data-driven time-stepping scheme for systems governed by
time-dependent PDEs involving multiscale physics. The proposed L-HiTS method
is based on the joint discovery of the coarse-grained variables and flow maps
inside the latent space, allowing for parallelization.

  2. 2.

The L-HiTS scheme incurs lesser computational costs both at the training stage
and testing stage than the multiscale HiTS method while delivering state-of-
the-art prediction accuracy.

  3. 3.

We demonstrate the applicability of the L-HiTS scheme for systems exhibiting
multiscale and chaotic behaviors, representing real-world phenomena.

  4. 4.

To make our work more accessible, we have provided open-source Python scripts
that can be used to reproduce the results.

The remainder of the manuscript is as follows. In Section 2, we briefly
introduce the multiscale analysis for predicting the solution of PDEs in a
data-driven fashion. Next, in Section 3, we present how deep AEs can be used
to learn a coordinate transformation that represents the actual dynamics onto
a latent space and then how the existing multiscale HiTS scheme can be
implemented in this latent space to achieve the task of prediction. We
describe the idea and present a numerical algorithm. Thereafter, in Section 4,
we show the application of the proposed scheme on two canonical PDEs. Finally,
in Section 5, we discuss the sensitivity analysis for different choices of
network parameters and present the conclusion of the study along with some
drawbacks of the method that can be taken up in a future study.

##  2 Problem setup

This section provides an overview of the general methodology for predicting
multiscale PDEs. To proceed with the idea, consider a nonlinear, continuous-
time, multiscale dynamical system of the form:

$\mathbf{u}_{t}=\mathcal{N}(\mathbf{u},\mathbf{u}_{x},\mathbf{u}_{xx},\dots,\mathbf{x}),$
(1)

where $t$ is the time $\in[0,T]$, $x\in[0,n]$ is the spatial variable and
$\mathcal{N}(.)$ is the nonlinear function of $\mathbf{u}(x,t)$ and its
derivatives. Using spatial discretization schemes such as the FEA, FDM, or the
spectral method [3], the PDE system (1) can be approximated by a system of
ODEs given by:

$\mathbf{\dot{x}}(t)=\mathcal{N}(\mathbf{x}(t;\boldsymbol{\mu})),,\quad\mathbf{x}(t_{0})=\mathbf{x}_{0},$
(2)

where $\mathbf{x}\in\mathbb{R}^{n}$ of dimension $n$ represents the state
vector, and $\boldsymbol{\mu}\in\mathbb{R}^{p}$ denotes the system parameters,
such as varying initial or boundary conditions. For a fine discretization, $n$
is typically high, resulting in a large-scale initial value problem to be
solved. Traditionally, the time solution of the ODE system (2) is obtained by
employing different time-stepping integration schemes, including the
explicit/implicit, fixed/adaptive, and one-step/multi-step methods with
varying degrees of stability [37]. These numerical methods take successive
time steps $\Delta t$ into the future in an iterative manner to construct the
discrete solutions for $x(t)$. Given the initial condition
$\mathbf{x}(t_{0})$, these methods approximate the discrete-time flow map

$\mathbf{x}(t+\Delta t)=\mathcal{N}(\mathbf{x}(t),\Delta
t)\triangleq\int_{\tau=t}^{t+\Delta
t}\mathcal{N}(\mathbf{x}(\tau),\tau)d\tau,$ (3)

usually via the Taylor-series expansion [38, 39, 40]. However, the accuracy of
these methods depends upon the local step size $\Delta t$ employed during the
integration, making such time discretization local in nature. For instance,
the classical ode45 method uses a fixed-step, fourth-order Runga-Kutta scheme
that incurs a local truncation error of $\mathcal{O}(\Delta t^{5})$ at every
step of the solution and a global truncation error of $\mathcal{O}(\Delta
t^{4})$. In contrast, neural networks-based time-stepper models like those in
[36, 33] learn a discrete flow map $\mathcal{N}()$ directly that progresses
the state variables forward in time without being restricted by local step
size limitations. Of the various neural network architectures available,
residual neural networks (ResNets) are commonly used because they have lesser
training times and work similar to that of a fixed-step, first-order Euler’s
scheme (see, e.g., Refs. [41, 33, 42, 43]).

ResNets are a special class of deep neural networks constructed by stacking
multiple residual blocks [43]. Each residual block comprises simple
feedforward layers that include identity skipped connections, allowing
information to bypass the residual layers (cf. Figure 1). ResNets have been
shown to outperform traditional feedforward neural networks (FNNs) in some
applications, particularly in computer vision and image classification tasks
[43]. The key advantage of ResNets is that they address the issue of vanishing
gradients, a common problem in deep FNNs [33].

Layer 1Layer 2Layer $\displaystyle
L$$\mathbf{x}(t_{0})$$\mathcal{N}(\mathbf{x}(t;\textbf{W}))$$\mathbf{x}(t_{0}+\Delta
t)=\mathbf{x}(t_{0})+\mathcal{N}(\mathbf{x}(t);W)$ResNet
Block$\mathbf{x}(t_{0})$$\mathbf{x}(t_{0}+\Delta
t)$$\mathcal{N}(\mathbf{x}(t;\mathbf{W}))$ Figure 1: Residual neural network
architecture

ResNets incorporate the identity operator into the network to have the network
effectively estimate the residual of the input-output mapping. The ResNet
block comprises $L$ hidden layers, along with an identity operator that adds
the input $\mathbf{x}(t)$ back into the output. The architecture describes the
following mapping;

$\mathbf{x}(t+\Delta
t)=\mathbf{x}(t_{0})+\mathcal{N}(\mathbf{x}(t);\mathbf{W}),\quad\mathbf{x}(t_{0})=\mathbf{x}_{0},$
(4)

where $\mathbf{W}$ represents the weights of the neural networks, and
$\mathcal{N}(:;\mathbf{W}):\mathbb{R}^{n}\rightarrow\mathbb{R}^{n}$ is the
nonlinear flow map learned by the network given as:

$\mathcal{N}(\mathbf{x}(t);\mathbf{W})=\mathbf{a}_{L}(\mathbf{W}_{L}(\dots\mathbf{a}_{1}(\mathbf{W}_{1}))\dots)),$
(5)

where $a_{k}~{}(k=1,..,L)$ is the activation function used at layer $k$, and
$W_{k}$ corresponds to weights for the $k^{th}$ layer. The weights are
obtained by minimizing the prediction loss function $\mathbf{W}^{*}_{pred}$
given as:

$\mathbf{W}^{*}_{pred}=\arg\min_{\mathbf{W}}\frac{1}{p\times
T}\sum_{i=1}^{p}\sum_{j=1}^{T}\left[\mathcal{L}_{pred}(\mathbf{\hat{x}}^{i}_{t+j\Delta
t},\mathbf{x}^{i}_{t+j\Delta t})\right],$ (6)

where $\mathcal{L}_{pred}$ denotes the discrepancy error between the true
state $\mathbf{x}(t)$ and the predicted state $\mathbf{\hat{x}}(t)$. It has
been well observed that ResNets work much similarly to a fixed step explicit
Euler scheme [41, 33]. This connection has been further exploited to develop
continuous normalizing flows [42].

However, when dealing with systems exhibiting multiscale phenomena, a single
discrete flow map operator $\mathcal{N}()$ may not provide an accurate
estimation across various time scales. As shown in [36], an NNTS model trained
on coarse time steps fails to capture fine-scale information, while the
prediction accuracy of the one trained on finer steps reduces rapidly due to
the autoregression-based estimation process. This issue is addressed by using
a combination of neural network models trained at different timescales, which
can effectively capture both coarse and fine-scale behaviors. However, this
approach requires extensive offline computational costs for training and
cross-validating the model combinations, especially when dealing with large
datasets. To overcome this drawback in the multiscale HiTS scheme, we propose
the L-HiTS framework. This data-driven scheme remains computationally
inexpensive while maintaining prediction accuracy at par with the multiscale
HiTS scheme. This can have practical applications for systems that involve
large datasets, such as PDEs. In the following, we describe the general
methodology and the algorithm for the L-HiTS method.

##  3 Latent Hierarchical Time-Stepping (L-HiTS) scheme for multiscale PDEs

In order to accurately and efficiently predict dynamics, we need to complete
two machine learning tasks. Firstly, we need to find a coordinate
transformation that can represent the original multiscale dynamics onto a
smaller subspace, using only a few latent variables. To achieve this, we use
deep AE networks. Secondly, we need to discover the discrete-time flow maps
for predictions within the latent space. This is done by using the multiscale
HiTS scheme on the latent variables. Both tasks are described in more detail
below.

###  3.1 Discovery of coordinates using deep autoencoders

Autoencoders are a special class of neural networks that learn to compress the
input data into a lower-dimensional representation and then reconstruct the
original data from the learned representation. In this study, we incorporate
deep AE models to learn a mapping
$\mathbf{f}:\mathbf{x}(t)\rightarrow\mathbf{\hat{x}}(t)$ from the input state
vector $\mathbf{x}(t)$ to its predicted output $\mathbf{\hat{x}}(t)$ such that
$\mathbf{\hat{x}}(t)\approx\mathbf{x}(t)$, where
$\mathbf{f}:\mathbb{R}^{n}\rightarrow\mathbb{R}^{n}$. We use the encoder
function of the autoencoder as a restricting operator that takes the input
state measurements $\mathbf{x}(t)$ from the dynamical system of dimension $n$
and map it to a latent state vector $\mathbf{z}(t)\in\mathbb{R}^{z}~{}(z<<n)$
given as:

$\mathbf{z}(t)=\mathbf{f}_{E}(\mathbf{x}(t);\mathbf{W}_{E}),$ (7)

where $\mathbf{f}_{E}()$ is the nonlinear function learned by the encoder, and
$\mathbf{W}_{E}$ are the weights of the encoder network. The latent vector
$\mathbf{z}(t)$ is then passed through a decoder function that acts as a
lifting operator, i.e., the decoder maps the latent vector $\mathbf{z}(t)$
back to the high-dimensional reconstructed state vector $\mathbf{\hat{x}}(t)$
given as:

$\mathbf{\hat{x}}(t)=\mathbf{f}_{D}(\mathbf{z}(t);\mathbf{W}_{D}),$ (8)

where $\mathbf{f}_{D}(\mathbf{z}(t);\mathbf{W}_{D})$ is the nonlinear function
learned by the decoder network, and $\mathbf{W}_{D}$ represents its weights.
Thus, an AE network combines the restricting and lifting mechanisms to produce
the reconstructed state $\mathbf{\hat{x}}$ given as:

$\mathbf{\hat{x}}(t)=\mathbf{f}_{D}(\mathbf{f}_{E}(\mathbf{x}(t);\mathbf{W}_{E}));\mathbf{W}_{D}).$
(9)

The reconstructed state is obtained by finding the optimal weight matrix
$\mathbf{W}=\mathbf{W}_{E},\mathbf{W}_{D}$ that minimizes the reconstruction
loss function $\mathbf{W}^{*}_{recon}$, given as:

$\mathbf{W}^{*}_{recon}=\arg\min_{\mathbf{W}}\frac{1}{N}\sum_{i=1}^{N}\left[\mathcal{L}_{recon}(\mathbf{\hat{x}}(t_{i}),\mathbf{x}(t_{i}))\right],$
(10)

where $\mathcal{L}_{recon}$ denotes the discrepancy error between
$\mathbf{x}(t)$ and $\mathbf{\hat{x}}(t)$. The condition $z<<n$ serves as a
regularization term that prevents the AE model from learning an identity
function. The schematic diagram of an AE is shown in Figure 2.

$\mathbf{x}(t)\in\mathbb{R}^{n}$$\mathbf{f}_{E}(\mathbf{x}(t);\mathbf{W}_{E})$$\hat{\mathbf{x}}(t)\in\mathbb{R}^{n}$$\mathbf{z}(t)\in\mathbb{R}^{z}$$\mathbf{f}_{D}(\mathbf{z}(t);\mathbf{W}_{D})$Encoder
functionDecoder functionRestricting OperatorLifting Operator Figure 2:
Schematic diagram of an autoencoder. The high dimensional input state
$\mathbf{x}(t)$ is restricted using an encoder function $\mathbf{f}_{E}$ to
obtain latent vector $\mathbf{z}(t)$, which is then passed through the lifting
mechanism of the decoder function $\mathbf{f}_{D}$ to obtain reconstructed
state vector $\hat{\mathbf{x}}(t)$.

The use of an autoencoder model to identify the coarse-grained dynamics has
several practical advantages. Firstly, the neural network-based design of the
model enables it to be scaled up to a large size, which can be difficult with
linear embedding-based techniques [44]. Secondly, the AE model describes both
the encoding and decoding transformations employed during the process. This
can be used within the EFF framework for multiscale applications, as
demonstrated in a recent study [30].

###  3.2 Discovery of multiscale flow maps using deep residual neural networks

Once the latent coordinates are obtained using a deep AE model, we proceed to
discover the flow-map representations of the latent variables. We accomplish
this by using the multiscale HiTS scheme (described in [36]) within the latent
space. The multiscale HiTS scheme involves learning different flow-map
representations across various temporal scales. This is done using a hierarchy
of deep neural network times steppers. The scheme uses deep ResNet models as
building blocks to learn different flow-maps that advance the system dynamics
forward in time. Given the data collected from the system at a unit time-step
$\Delta t$ sub-sampled uniformly at different multiples of $\Delta t$, i.e.,
$\Delta t_{m}=2^{m-d}\Delta t~{}(d=1,2,\dots,m)$, the multiscale HiTS scheme
learns multiple flow maps $\mathcal{N}_{m}(,;\Delta t_{m})$ that captures the
system dynamics across timescales. The NNTS model with a large temporal gap
captures the slow scale behavior, while the NNTS model with the smallest
temporal gap captures the fast scale features of the system.

After learning the different flow map representations, the scheme couples the
output of different NNTS models RN1, RN2,$\dots$, RNm using an iterative
vectorized fashion (see Algorithm 1 of [36]). The process involves using the
NNTS model with the largest temporal gap first to make predictions. These
predictions are then passed on to the next level NNTS model, and its outputs
are stacked together with the previous one. This process continues until the
smallest NNTS model is used. After rearranging the predictions, an
interpolation scheme is used to determine the predictions between time steps.
By doing this, the error propagation from individual NNTS models can be
eliminated. Before coupling the individual NNTS models for multiscale
applications, the HiTS scheme uses a cross-validation step to filter out the
best possible combination of NNTS models. This step is essential as the
accuracy of individual NNTS models may vary, which could lead to inaccurate
predictions after coupling. However, this process can be computationally
expensive, especially for large-scale systems where the sampled data is
extensive. To address this shortcoming, an adaptive hierarchical time stepping
(AHiTS) method has been proposed [45] which adapts the time-step, i.e., the
choice of NNTS model during the simulation, based on system dynamics.

Figure 3: Block diagram of the proposed L-HiTS method

###  3.3 Combining deep AEs and multiscale HiTS scheme for multiscale
predictions

The proposed L-HiTS scheme combines the use of coordinates and flow maps for
multiscale predictions within the latent space. We begin by collecting $p$
trajectories of the dynamical system state $\mathbf{x}(t)$ at $t_{f}$
instances of the unit time step $\Delta t$ and stack inside a matrix
$\mathbf{X}$ given as:

$\mathbf{X}^{i}=[\mathbf{x}^{i}(t_{0})~{}\quad\mathbf{x}^{i}(t_{0}+\Delta
t)~{}\dots~{}\mathbf{x}^{i}({t_{0}}+t_{f}\Delta t)~{}],$ (11)

where $i\in 1,\dots,p$. The data is divided into training, validation, and
testing datasets. The training data is compressed using the restricting action
of the encoder function within the AE model to obtain the latent variables
$\mathbf{z}(t)$. These latent variables are then fed to individual NNTS models
RN1 to RNm as per the respective step sizes. For instance, the RN1 model uses
the first $t_{f}-1$ entries of $\mathbf{z}(t)$, i.e.,
$[\mathbf{z}^{i}(t_{0}),...,\mathbf{z}^{i}({t_{0}}+(t_{f}-1)\Delta t)]$ as the
input and one $\Delta t$ ahead entries, i.e., $[\mathbf{z}^{i}(t_{0}+\Delta
t),\dots,\mathbf{z}^{i}({t_{0}}+t_{f}\Delta t)]$ as the output. By doing so,
each NNTS model learns its respective flow map
$\mathcal{N}_{m}({\mathbf{z}(t)};{\Delta t_{m}})$ to obtain the future
predictions of the latent variables from a given initial condition $t_{0}$,
which is given as:

$\mathbf{z}({t_{0}+{\Delta
t_{m}}})=\mathbf{z}(t_{0})+\mathcal{N}_{m}({\mathbf{z}(t)};{\Delta t_{m}}).$
(12)

After the individual responses are obtained, we use the validation dataset to
cross-validate the NNTS models for coupling, as described above. Using
iterative vectorized computations from [36], we obtain the future predictions
in the latent space. The solution thus obtained is finally lifted using the
decoder function within the AE model. Figure 3 provides a summary of the
L-HiTS scheme. We can use the L-HiTS scheme to obtain future state predictions
for any testing initial condition $\mathbf{x}_{0}$. For an unknown system
parameter $\mathbf{\mu}_{test}$, we pass the testing initial condition through
the encoder block to obtain the initial condition in the latent space
$\mathbf{z}_{0}$. We then utilize the multiscale HiTS framework to initialize
the iterative one-step ahead predictions within the latent space for a given
time span $[t_{0},T]$. The predictions within the latent space are then
obtained directly using the multiscale HiTS scheme. Finally, we pass the
solutions obtained through the decoder function to reconstruct the original
state variable $\hat{\mathbf{x}(t)}$. The pseudocode for the L-HiTS procedure
can be found in Algorithm 1. In the next section, we will demonstrate the
application of the L-HiTS scheme using numerical simulations.

Algorithm 1 : L-HiTS Algorithm

Input: ResNet models: RN1,.., RNm, initial condition $\mathbf{x}_{0}$
Output: Reconstructed states $\hat{\mathbf{x}}.$

1: $\mathit{sort}$(RN models) $\circ$ in descending order as per trained step
sizes.

2: $\mathbf{Z}=list\verb|{}|$

3: $\mathbf{z}=f_{E}(\mathbf{x}_{0})$ $\circ$ Encode the initial condition

4: $\mathit{append}(\mathbf{Z},\mathbf{z})$ $\circ$ Append the encoded initial
condition to $\mathbf{z}$

5: for k=1,2,…,m do

6: $Q_{k}:=$ number of forward steps.

7: $\mathbf{z}_{c}=\mathit{stack}(\mathbf{z})$ $\circ$ stack states together

8: $\mathbf{z}_{n}=RN_{k}.\mathit{forward}(\mathbf{z}_{c},\mathbf{z})$ $\circ$
use Algorithm 1 of [36] to obtain $\mathbf{z}_{n}$

9: $\mathit{append}(\mathbf{z}_{n},\mathbf{z})$ $\circ$ update the predicted
states

10: end for

11: $\mathit{rearrange}(\mathbf{z})$

12: $\mathit{Interpolate}(\mathbf{z})$

13: $\hat{\mathbf{X}}=f_{D}(\mathbf{z})$ $\circ$ Decode the reduced states

14: return($\hat{\mathbf{x}}$) $\circ$ Return the reconstructed states

##  4 Simulation results

To demonstrate the application of the L-HiTS scheme, we consider two benchmark
PDEs. The FitzHugh-Nagumo (FHN) model [46, 47] and the 1D Kuramoto-Sivashinsky
(KS) equation [48, 49]. To facilitate comparison with the multiscale HiTS
method, all the model parameters and data structures are kept the same. To
test the accuracy of the methods, we use the mean squared error (MSE) as a
metric for comparison. All the simulations are performed using the open source
Python API for PyTorch framework running on hp workstation Z1 with $11^{th}$
Gen Intel(R) Core (TM) $i7$-$11700$<EMAIL_ADDRESS>CPU. However, the code that we
accompany with this manuscript is written to take advantage of faster
computations via GPU if available.

###  4.1 The FitzHugh-Nagumo model

Figure 4: Results for the FHN model: (Plot A:) plot of the MSEs for
autoencoder and PCA for testing dataset against varying latent dimensions
(Plot B-C:) evolution of latent variables $\mathbf{z}_{1}(t)$ and
$\mathbf{z}_{2}(t)$ obtained from AE and L-HiTS scheme, (Plot D:) true state
$\mathbf{u}(\mathbf{x},t)$, (Plot E:) predicted state
$\mathbf{\hat{u}}(\mathbf{x},t)$ from the L-HiTS scheme, (Plot F:) MSE between
$\mathbf{\hat{u}}(\mathbf{x},t)$ and $\mathbf{\hat{u}}(\mathbf{x},t)$, (Plot
G:) true state $\mathbf{v}(\mathbf{x},t)$, (Plot H:) predicted state
$\mathbf{\hat{v}}(\mathbf{x},t)$ from the L-HiTS scheme, (Plot I:) MSE between
$\mathbf{\hat{v}}(\mathbf{x},t)$ and $\mathbf{\hat{v}}(\mathbf{x},t)$.

The FHN model is a simplified 2D version of the Hodgkin-Huxley model that
describes the activation and deactivation dynamics of a spiking neuron. The
membrane potential is described by a single variable,
$\mathbf{u}(\mathbf{x},t)$, known as the inhibitor, and the
activation/deactivation dynamics of the voltage-gated ion channels are
represented by another variable, $\mathbf{v}(\mathbf{x},t)$, known as the
activator. The governing equation consists of two coupled nonlinear
differential equations that describe the time evolution of
$\mathbf{u}(\mathbf{x},t)$ and $\mathbf{v}(\mathbf{x},t)$ given as;

$\begin{split}\epsilon\frac{\partial\mathbf{u}}{\partial
t}(\mathbf{x},t)&=\epsilon^{2}\frac{\partial\mathbf{u}^{2}}{\partial\mathbf{x}^{2}}(\mathbf{x},t)+\mathbf{f}(\mathbf{u}(\mathbf{x},t)-\mathbf{v}(\mathbf{x},t)+0.05,\\\
\frac{\partial\mathbf{v}}{\partial
t}(\mathbf{x},t)&=0.5\mathbf{u}(\mathbf{x},t)-2\mathbf{v}(\mathbf{x},t)+0.05,\end{split}$
(13)

where $\mathbf{f}(\mathbf{u})=\mathbf{u}(\mathbf{u}-0.1)(1-\mathbf{u})$. The
model exhibits multiscale behavior with activator $\mathbf{u}(\mathbf{x},t)$
showing fast dynamics, whereas inhibitor $\mathbf{v}(\mathbf{x},t)$ has slower
dynamics. The difference is controlled by the parameter $\epsilon=0.015$. The
initial and boundary conditions are taken from [50]. Using Lattice Boltzmann
(LB) method [51], (13) is discretized using $n=202$, grid points in space and
$T=5120$ grid points in time. We utilize the LB method to collect data from
six different initial conditions in order to derive the true solution of the
system similar to [50, 30]. We consider four initial conditions for training
purposes, one for validation, and a different initial condition is used for
testing the accuracy of the L-HiTS method. To identify the latent space, we
use the principle component analysis (PCA) method and AE by comparing the MSE
of reconstruction on the test data. Figure 4, plot A shows the comparison in
terms of MSE against the number of latent dimensions $z$. As can be seen, AE
outperforms PCA by identifying a nonlinear subspace that effectively reduces
the dimensionality to $z=2$, and then the latent dimension size remains
constant. This indicates that the underlying solution space has an intrinsic
dimensionality of order two. The various network parameters for this case are
enlisted in Appendix A, Table 5. The plots B and C of the same figure show the
evolution of the two latent variables obtained from the AE model and the
corresponding predictions from the L-HiTS scheme. As can be seen, the true and
predicted dynamics are graphically indistinguishable. Numerically, the
accuracy of the L-HiTS scheme for varying latent dimensions for both latent
and reconstructed variables can be seen in Table 1. The results show that the
MSE plateaus after the latent dimension of $z=2$, which is also reflected in
the MSEs of the output predictions and we use this for initialization of the
L-HiTS scheme. The plots D-I of Figure 4 show the corresponding space-time
plots for both the variables $\mathbf{u}(\mathbf{x},t)$ and
$\mathbf{v}(\mathbf{x},t)$ which shows a good agreement between the true and
the predicted solution. The computational costs of the L-HiTS method are shown
in Table 2, where the efficiency of the L-HiTS scheme against the multiscale
HiTS scheme can be observed both during the training and prediction stages.

###  4.2 The 1D Kuramoto-Sivashinsky equation

Figure 5: Results for the 1D KS equation: (Plot A:) plot of MSEs for
autoencoder and PCA for testing dataset against varying latent dimensions,
(Plot B-I:) evolution of the latent variables
$\mathbf{z}_{1}(t),\dots,\mathbf{z}_{8}(t)$ obtained from the AE model and
L-HiTS scheme, (Plot J:) true solution of the state
$\mathbf{u}(\mathbf{x},t)$, (Plot K:) predicted state
$\mathbf{\hat{u}}(\mathbf{x},t)$ from the L-HiTS scheme, (Plot L:) MSE between
$\mathbf{u}(\mathbf{x},t)$ and $\mathbf{\hat{u}}(\mathbf{x},t)$

To test the proposed L-HiTS scheme on a more challenging system, we take the
1D KS equation that exhibits spatiotemporal chaos and is used to describe
various physical phenomena such as chemical reactions, turbulence, and flame
fonts. The system dynamics are described by a fourth-order nonlinear PDE
defined as:

$\mathbf{u}_{t}=-\mathbf{u}\mathbf{u}_{x}-\frac{1}{2}\mathbf{u}_{xx}-\mathbf{u}_{xxxx}$
(14)

for $x\in[-L/2,L/2]$ with $n=120$ grids points and a periodic boundary
condition $\mathbf{u}(0;t)=\mathbf{u}(L;t)$ and some given initial conditions.
We simulate the KS equation using a fourth-order stiff time-stepping (ETDRK4)
method [52] with $t\in[0,5121]$ with a step size of $\Delta t=0.05$. We
consider the case of $L=22$ as considered in [53, 30] for which the system
exhibits a structurally stable chaotic attractor,i.e., an inertial manifold in
which the long-term behavior of the system lies. For the training phase of the
L-HiTS method, we take $\mu_{train}=10$ random, periodic initial conditions,
$\mu_{val}=5$ initial conditions for the validation, and $\mu_{test}=5$ random
initial conditions for testing purposes. For the identification of the latent
space, we compare the PCA and AE models in terms of reconstruction MSE in the
test data as a function of latent dimension $z$, graphically shown in plot A
of Figure 5 and numerically enlisted in Table 1. We observe that the MSE
plateaus after $z=8$, in agreement with [53, 30]. Thus, we chose a latent
dimension of $z=8$ for the L-HiTS framework. We also compared the different
latent variables obtained by the L-HiTS scheme against the testing data
encoded within the latent space. These are shown in plots B-I of Figure 5, and
the corresponding hyper-parameters, including the network architectures of the
AE and different ResNet models, are enlisted in Appendix A, Table 6. As
observed, the L-HiTS method correctly estimates the evolution of the latent
variables in time.

To predict the chaotic dynamics within the latent space, we use $m=11$ NNTS
models trained on the outputs from the encoder. After coupling the individual
NNTS models, the output is reconstructed using the decoder function within the
AE model. The original and reconstructed space-time plots are shown in plots J
and K of Figure 5 respectively, with the error in plot L. The results show
that the L-HiTS scheme satisfactorily reconstructs the true chaotic dynamics
of the KS equation. Furthermore, the computational costs of the L-HiTS scheme
against the multiscale HiTS scheme are presented in Table 2. As expected, the
L-HiTS method yields almost one order of magnitude saving than the multiscale
HiTS scheme in predicting dynamics and a lesser training time while
maintaining the same accuracy. This substantiates our observation that
learning multiple flow maps of the system within the latent space yields an
efficient state prediction in the original solution space.

Table 1: MSEs for the FHN model and KS equation obtained using the L-HiTS method for varying dimensions of latent space Latent dimension | FHN model | KS equation
---|---|---
Latent space | Reconstruction | Latent space | Reconstruction
| (MSE) | (MSE) | (MSE) | (MSE)
1 | $2.34\text{\times}{10}^{0}$ | $2.44\text{\times}{10}^{-1}$ | $6.79\text{\times}{10}^{1}$ | $2.26\text{\times}{10}^{0}$
2 | $1.47\text{\times}{10}^{-3}$ | $1.37\text{\times}{10}^{-4}$ | $5.01\text{\times}{10}^{1}$ | $2.41\text{\times}{10}^{0}$
4 | $2.42\text{\times}{10}^{-2}$ | $2.20\text{\times}{10}^{-3}$ | $4.73\text{\times}{10}^{1}$ | $1.30\text{\times}{10}^{0}$
6 | $2.79\text{\times}{10}^{-2}$ | $3.60\text{\times}{10}^{-3}$ | $1.90\text{\times}{10}^{-2}$ | $4.05\text{\times}{10}^{-3}$
8 | $2.13\text{\times}{10}^{-3}$ | $2.88\text{\times}{10}^{-4}$ | $1.35\text{\times}{10}^{-3}$ | $1.50\text{\times}{10}^{-3}$
10 | $5.96\text{\times}{10}^{-3}$ | $9.97\text{\times}{10}^{-4}$ | $3.16\text{\times}{10}^{-3}$ | $1.42\text{\times}{10}^{-3}$
Table 2: Comparison of accuracy and computational cost Technique | FHN model | KS equation
---|---|---
Training | Prediction |  $\ell_{2}$-norm | Training | Prediction |  $\ell_{2}$-norm
| time | time | error | time | time | error
Multiscale-HiTS [36] | $1767s$ | $8.08s$ |  $2\times 10^{-4}$ | $1982s$ | $20.88s$ | $4\times 10^{-3}$
L-HiTS (proposed) | $1248s$ | $3.52s$ |  $2\times 10^{-4}$ | $1624s$ | $2.86s$ | $4\times 10^{-3}$

##  5 Discussions and conclusion

In this section, we provide some additional comments regarding the sensitivity
of the L-HiTS scheme for different network choices, followed by some
concluding remarks toward the end.

###  5.1 Effect of different latent dimensions

To evaluate the sensitivity of the L-HiTS method for different latent sizes,
we calculate the mean square errors for both the latent and reconstruction
states, as shown in Table 1. Additionally, we plot the MSEs for different time
steps during the forecasting stage, as presented in Figure 6. Our results
indicate that for both test cases, increasing the latent dimension of the
autoencoder model reduces the mean square error. This means that a larger
dimension of the latent space implies that there are more basis functions to
express the encoded information about the system, which is in agreement with
[23]. However, we found that the MSE saturates beyond $z=2$ for the FHN system
and $z=8$ for the case of the KS equation, indicating that no further
improvement in forecasting accuracy can be achieved.

In the context of machine learning, the relationship between accuracy in
latent space and reconstruction accuracy is complex and can vary depending on
the model architecture and the specific application. In general, there isn’t a
direct correlation where higher accuracy in latent space straightforwardly
results in higher reconstruction accuracy. However, our findings are aligned
with [23], wherein authors have established a similar relationship between the
latent dimensions and reconstruction accuracy. The authors therein report a
reduced $\mathcal{L}^{2}$ reconstruction error as the latent dimension of the
autoencoder model is increased. Besides, there are contexts and specific
methods where improvements in latent representation can enhance reconstruction
quality, such as the use of sparse and denoising autoencoders [54, 53],
variational autoencoders [55], and deep embedding clustering[56]. While our
framework currently incorporates a basic version of the autoencoder, it can
easily incorporate these variations of autoencoders with only minor
adjustments. This will ensure that the relationship between latent space and
reconstruction will remain consistent.

010002000300040005000$10^{-5}$$10^{-3}$$10^{-1}$Time step $(t)$MSE
($\mathbf{x}-\hat{\mathbf{x}}$)$\mathbf{z}=1$$\mathbf{z}=2$$\mathbf{z}=4$$\mathbf{z}=6$$\mathbf{z}=8$$\mathbf{z}=10$
(a) FHN model

010002000300040005000$10^{-3}$$10^{-2}$$10^{-1}$$10^{0}$Time step $(t)$MSE
($\mathbf{x}-\hat{\mathbf{x}}$) (b) KS equation

Figure 6: Sensitivity analysis: mean square error of state predictions at different time-steps for varying latent space dimension $z$. Table 3: FHN model: comparison of MSE and CPU times across various NNTS models at a latent dimension of $z=2$ Model | MSE ($x-\hat{x}$) | CPU times (seconds)
---|---|---
RN1 | $3.73\times 10^{-1}$ | 1.22
RN2 | $3.55\times 10^{-1}$ | 0.66
RN4 | $2.36\times 10^{-1}$ | 0.29
RN8 | $3.52\times 10^{-3}$ | 0.14
RN16 | $2.26\times 10^{-3}$ | 0.07
RN32 | $8.07\times 10^{-3}$ | 0.04
RN64 | $6.04\times 10^{-2}$ | 0.01
RN128 | $2.58\times 10^{-1}$ | 0.01
RN256 | $2.92\times 10^{-1}$ | 0.006
RN512 | $3.14\times 10^{-1}$ | 0.003
RN1024 | $2.85\times 10^{-1}$ | 0.002
L-HiTS | $1.37\times 10^{-4}$ | 3.52
Table 4: KS equation: comparison of MSE and CPU times across various NNTS models at a latent dimension of $z=8$. Model | MSE ($x-\hat{x}$) | CPU times (seconds)
---|---|---
RN1 | $1.97$ | 0.79
RN2 | $2.77$ | 0.37
RN4 | $2.68$ | 0.19
RN8 | $2.47$ | 0.09
RN16 | $1.67$ | 0.04
RN32 | $1.10$ | 0.02
RN64 | $1.09\times 10^{-2}$ | 0.01
RN128 | $8.74\times 10^{-2}$ | 0.009
RN256 | $4.04\times 10^{-1}$ | 0.0055
RN512 | $9.69\times 10^{-1}$ | 0.0037
RN1024 | $1.29$ | 0.0029
L-HiTS | $1.50\times 10^{-3}$ | 2.86

###  5.2 Effect of individual and multiple models within the latest space

To check the effect of individual NNTS models on the quality of output
prediction and computational costs, we use each individual NNTS model within
the latent space and compare that with the coupling framework used in the
L-HiTS scheme. The results are shown in Tables 3 and 4. As observed from the
results, models with small time steps yield accurate short-term predictions;
however, the error quickly accumulates due to autoregressive predictions,
making the overall reconstruction quite poor. In contrast, networks with
larger temporal gaps yield better long-term predictions but fail to capture
information between steps. Thus, models with intermediate time steps, such as
RN16 in the case of the FHN model and RN64 in the case of the KS equation,
balance out the two factors, resulting in a better MSE that the proposed
method effectively captures, thus outperforming any individual NNTS model. The
computation time, on the other hand, accelerates as the step size grows, and
the proposed method provides a good trade-off between accuracy and efficiency.

###  5.3 Concluding remarks

In this study, we propose a novel framework that enhances the computational
efficiency of multiscale PDEs. Our strategy is based on the idea of learning
coordinates and flow maps using deep neural networks that enable rapid and
accurate predictions of future states. The proposed L-HiTS method uses deep
autoencoder networks to learn suitable latent coordinates that represent the
dynamics on a low-dimensional basis. Then, we make use of the existing
multiscale HiTS method within the discovered latent space to obtain multiple
flow maps that march the latent variables in time. The predictions are then
lifted to the original space using the decoder function of the AE model. The
method outperforms any individual neural network time stepper model in terms
of accuracy and the original multiscale HiTS method in terms of computational
costs. We demonstrate the application of the proposed scheme on two large-
scale PDEs with multiscale physics.

Although the proposed method works well, there are some potential drawbacks
that can be taken up in a future study. Firstly, the proposed method assumes
full-state observations, meaning that all the time-dependent states are
available across a diverse range of operating conditions. However, in many
practical scenarios, only partial measurements are available. To address this
issue, a time delay embedded of the partial measurements can be used ( see,
e.g., [57, 58]) within the L-HiTS framework to make it more practical.
Secondly, the L-HiTS framework uses fixed step size-based predictions of the
latent states within the individual ResNet models. This can be computationally
expensive if some of the latent variables evolve faster than the rest. To
address this issue, an adaptive time-step-based strategy can be used within
the latent space, where the time step can be adjusted based on the evolving
dynamics.

## Credit authorship contribution statement

Asif Hamid: Software, Methodology, Data-curation, Validation, Writing-
original draft,
Danish Rafiq: Conceptualization, Methodology, Funding acquisition,
Investigation, Writing- Reviewing and Editing,
Shahkar A. Nahvi: Supervision, Reviewing,
Mohammad A. Bazaz: Supervision, Reviewing

## Declaration of competing interest

The authors declare no potential financial or non-financial competing
interests

## Acknowledgments

Author I would like to acknowledge the Doctoral fellowship from the Islamic
University of Science and Technology, Kashmir, India, via Grant No.
IUST0119013135, and Author II, acknowledges the financial assistance from the
Science and Engineering Research Board (SERB), a statutory body of the
Department of Science and Technology (DST), Government of India, via Grant No.
PDF/2022/002081.

## Code availability

The source code is available at https://github.com/erasifhamid/L-HiTS

## References

  * [1] J. H. Mathews, Numerical methods for mathematics, science and engineering, Vol. 10, Prentice-Hall International, 1992.
  * [2] J. Bezanson, A. Edelman, S. Karpinski, V. B. Shah, Julia: A fresh approach to numerical computing, SIAM review 59 (1) (2017) 65–98.
  * [3] S. D. Conte, C. De Boor, Elementary numerical analysis: an algorithmic approach, SIAM, 2017.
  * [4] J. W. Thomas, Numerical partial differential equations: finite difference methods, Vol. 22, Springer Science & Business Media, 2013.
  * [5] E. Weinan, Principles of multiscale modeling, Vol. 1, Cambridge University Press, 2011.
  * [6] O. Ghattas, K. Willcox, Learning physics-based models from data: perspectives from inverse problems and model reduction, Acta Numerica 30 (2021) 445–554.
  * [7] W. Gong, Q. Duan, J. Li, C. Wang, Z. Di, Y. Dai, A. Ye, C. Miao, Multi-objective parameter optimization of common land model using adaptive surrogate modeling, Hydrology and Earth System Sciences 19 (5) (2015) 2409–2425.
  * [8] E. A. Groen, R. Heijungs, E. A. Bokkers, I. J. De Boer, Methods for uncertainty propagation in life cycle assessment, Environmental Modelling & Software 62 (2014) 316–325.
  * [9] M. Taufer, E. Deelman, R. F. da Silva, T. Estrada, M. Hall, M. Livny, A roadmap to robust science for high-throughput applications: The developers’ perspective, in: 2021 IEEE International Conference on Cluster Computing (CLUSTER), IEEE, 2021, pp. 807–808.
  * [10] G. Novati, H. L. de Laroussilhe, P. Koumoutsakos, Automating turbulence modelling by multi-agent reinforcement learning, Nature Machine Intelligence 3 (1) (2021) 87–96.
  * [11] S. Verma, G. Novati, P. Koumoutsakos, Efficient collective swimming by harnessing vortices through deep reinforcement learning, Proceedings of the National Academy of Sciences 115 (23) (2018) 5849–5854.
  * [12] J. J. Bramburger, D. Dylewsky, J. N. Kutz, Sparse identification of slow timescale dynamics, Physical Review E 102 (2) (2020) 022204.
  * [13] Kevrekidis, et al., Equation-free, coarse-grained multiscale computation: enabling microscopic simulators to perform system-level analysis, Commun. Math. Sci 1 (4) (2003) 715–762.
  * [14] I. G. Kevrekidis, C. W. Gear, G. Hummer, Equation-free: The computer-aided analysis of complex multiscale systems, AIChE Journal 50 (7) (2004) 1346–1355.
  * [15] E. Weinan, B. Engquist, The heterognous multiscale methods, Communications in Mathematical Sciences 1 (1) (2003) 87–132.
  * [16] M. Tao, H. Owhadi, J. E. Marsden, Nonintrusive and structure preserving multiscale integration of stiff odes, sdes, and hamiltonian systems with hidden slow dynamics via flow averaging, Multiscale Modeling & Simulation 8 (4) (2010) 1269–1324.
  * [17] R. Ganesan, T. K. Das, V. Venkataraman, Wavelet-based multiscale statistical process monitoring: A literature review, IIE transactions 36 (9) (2004) 787–806.
  * [18] J. N. Kutz, X. Fu, S. L. Brunton, Multiresolution dynamic mode decomposition, SIAM Journal on Applied Dynamical Systems 15 (2) (2016) 713–735.
  * [19] S. G. Mallat, A theory for multiresolution signal decomposition: the wavelet representation, IEEE transactions on pattern analysis and machine intelligence 11 (7) (1989) 674–693.
  * [20] R. Vinuesa, S. L. Brunton, Enhancing computational fluid dynamics with machine learning, Nature Computational Science 2 (6) (2022) 358–366.
  * [21] Z. Zhang, R. Rai, S. Chowdhury, D. Doermann, Midphynet: Memorized infusion of decomposed physics in neural networks to model dynamic systems, Neurocomputing 428 (2021) 116–129.
  * [22] N. A. K. Doan, W. Polifke, L. Magri, Physics-informed echo state networks, Journal of Computational Science 47 (2020) 101237.
  * [23] V. Oommen, K. Shukla, S. Goswami, R. Dingreville, G. E. Karniadakis, Learning two-phase microstructure evolution using neural operators and autoencoder architectures, npj Computational Materials 8 (1) (2022) 190.
  * [24] T. Gracious, S. Gupta, A. Kanthali, R. M. Castro, A. Dukkipati, Neural latent space model for dynamic networks and temporal knowledge graphs, in: Proceedings of the AAAI Conference on Artificial Intelligence, Vol. 35, 2021, pp. 4054–4062.
  * [25] J. Fish, G. J. Wagner, S. Keten, Mesoscopic and multiscale modelling in materials, Nature materials 20 (6) (2021) 774–786.
  * [26] M. Meier-Schellersheim, I. D. Fraser, F. Klauschen, Multiscale modeling for biologists, Wiley Interdisciplinary Reviews: Systems Biology and Medicine 1 (1) (2009) 4–14.
  * [27] M. Garbey, M. Rahman, S. Berceli, A multiscale computational framework to understand vascular adaptation, Journal of computational science 8 (2015) 32–47.
  * [28] C. Wehmeyer, F. Noé, Time-lagged autoencoders: Deep learning of slow collective variables for molecular kinetics, The Journal of chemical physics 148 (24) (2018) 241703.
  * [29] P. R. Vlachas, J. Zavadlav, M. Praprotnik, P. Koumoutsakos, Accelerated simulations of molecular systems through learning of effective dynamics, Journal of Chemical Theory and Computation 18 (1) (2021) 538–549.
  * [30] P. R. Vlachas, G. Arampatzis, C. Uhler, P. Koumoutsakos, Multiscale simulations of complex systems by learning their effective dynamics, Nature Machine Intelligence 4 (4) (2022) 359–366.
  * [31] I. Kičić, P. R. Vlachas, G. Arampatzis, M. Chatzimanolakis, L. Guibas, P. Koumoutsakos, Adaptive learning of effective dynamics for online modeling of complex systems, Computer Methods in Applied Mechanics and Engineering 415 (2023) 116204.
  * [32] J. J. Bramburger, J. N. Kutz, Poincaré maps for multiscale physics discovery and nonlinear floquet theory, Physica D: Nonlinear Phenomena 408 (2020) 132479.
  * [33] T. Qin, K. Wu, D. Xiu, Data driven governing equations approximation using deep neural networks, Journal of Computational Physics 395 (2019) 620–635.
  * [34] L. Ying, E. J. Candes, The phase flow method, Journal of Computational Physics 220 (1) (2006) 184–215.
  * [35] S. Recanatesi, M. Farrell, G. Lajoie, S. Deneve, M. Rigotti, E. Shea-Brown, Predictive learning as a network mechanism for extracting low-dimensional latent space representations, Nature communications 12 (1) (2021) 1417.
  * [36] Y. Liu, J. N. Kutz, S. L. Brunton, Hierarchical deep learning of multiscale differential equation time-steppers, Philosophical Transactions of the Royal Society A 380 (2229) (2022) 20210200.
  * [37] K. Atkinson, W. Han, D. E. Stewart, Numerical solution of ordinary differential equations, John Wiley & Sons, 2011.
  * [38] J. C. Butcher, The numerical analysis of ordinary differential equations: Runge-Kutta and general linear methods, Wiley-Interscience, 1987.
  * [39] J. Guckenheimer, P. Holmes, Nonlinear oscillations, dynamical systems, and bifurcations of vector fields, Vol. 42, Springer Science & Business Media, 2013\.
  * [40] S. Wiggins, S. Wiggins, M. Golubitsky, Introduction to applied nonlinear dynamical systems and chaos, Vol. 2, Springer, 2003.
  * [41] B. Chang, L. Meng, E. Haber, F. Tung, D. Begert, Multi-level residual networks from dynamical systems view, arXiv preprint arXiv:1710.10348 (2017).
  * [42] R. T. Chen, Y. Rubanova, J. Bettencourt, D. K. Duvenaud, Neural ordinary differential equations, Advances in neural information processing systems 31 (2018).
  * [43] K. He, X. Zhang, S. Ren, J. Sun, Deep residual learning for image recognition, in: Proceedings of the IEEE conference on computer vision and pattern recognition, 2016, pp. 770–778.
  * [44] T. Simpson, N. Dervilis, E. Chatzi, Machine learning approach to model order reduction of nonlinear systems via autoencoder and lstm networks, Journal of Engineering Mechanics 147 (10) (2021) 04021061.
  * [45] A. Hamid, D. Rafiq, S. A. Nahvi, M. A. Bazaz, Hierarchical deep learning-based adaptive time-stepping scheme for multiscale simulations, arXiv preprint arXiv:2311.05961 (2023).
  * [46] R. FitzHugh, Impulses and physiological states in theoretical models of nerve membrane, Biophysical journal 1 (6) (1961) 445–466.
  * [47] T. Nakamura, K. Fukami, K. Hasegawa, Y. Nabae, K. Fukagata, Convolutional neural network and long short-term memory based reduced order surrogate for minimal turbulent channel flow, Physics of Fluids 33 (2) (2021) 025116.
  * [48] Y. Kuramoto, Diffusion-induced chaos in reaction systems, Progress of Theoretical Physics Supplement 64 (1978) 346–367.
  * [49] G. I. Sivashinsky, Nonlinear analysis of hydrodynamic instability in laminar flames—i. derivation of basic equations, Acta astronautica 4 (11) (1977) 1177–1206.
  * [50] S. Lee, M. Kooshkbaghi, K. Spiliotis, C. I. Siettos, I. G. Kevrekidis, Coarse-scale pdes from fine-scale observations via machine learning, Chaos: An Interdisciplinary Journal of Nonlinear Science 30 (1) (2020).
  * [51] I. V. Karlin, S. Ansumali, C. E. Frouzakis, S. S. Chikatamarla, Elements of the lattice boltzmann method i: Linear advection equation, Commun. Comput. Phys 1 (4) (2006) 616–655.
  * [52] H. Montanelli, N. Bootland, Solving periodic semilinear stiff pdes in 1d, 2d and 3d with exponential integrators, Mathematics and Computers in Simulation 178 (2020) 307–327.
  * [53] P. Cvitanović, R. L. Davidchack, E. Siminos, On the state space geometry of the kuramoto–sivashinsky flow in a periodic domain, SIAM Journal on Applied Dynamical Systems 9 (1) (2010) 1–33.
  * [54] P. Vincent, H. Larochelle, Y. Bengio, P.-A. Manzagol, Extracting and composing robust features with denoising autoencoders, in: Proceedings of the 25th international conference on Machine learning, 2008, pp. 1096–1103.
  * [55] D. P. Kingma, J. Ba, Adam: A method for stochastic optimization, arXiv preprint arXiv:1412.6980 (2014).
  * [56] J. Xie, R. Girshick, A. Farhadi, Unsupervised deep embedding for clustering analysis, in: International conference on machine learning, PMLR, 2016, pp. 478–487.
  * [57] M. Kamb, E. Kaiser, S. L. Brunton, J. N. Kutz, Time-delay observables for koopman: Theory and applications, SIAM Journal on Applied Dynamical Systems 19 (2) (2020) 886–917.
  * [58] S. Pan, K. Duraisamy, On the structure of time-delay embedding in linear models of non-linear dynamical systems, Chaos: An Interdisciplinary Journal of Nonlinear Science 30 (7) (2020).

##  Appendix A Hyperparameters for L-HiTS method

Table 5: Various hyper-parameters for the FHN model Parameters | Values
---|---
Spatial dimension $n$ | 202
Temporal dimension $T$ | 5120
Time step $\Delta t$ | 0.01
Latent dimensions $z$ | 1, 2, 4, 6, 8, 10
AE architecture | 100: 100: 100: $z$: 100: 100: 100
No. of epochs for training AE models | 5000
Activation function | ReLU
Batch size | 32
Learning rate | $1e-03$
No. of epochs for training ResNet models | 20000
Step sizes for ResNet models |  $\Delta$t, 2$\Delta$t, 4$\Delta$t, 8$\Delta$t, 16$\Delta$t, 32$\Delta$t, 64$\Delta$t, 128$\Delta$t, 256$\Delta$t, 512$\Delta$t, 1024$\Delta$t
ResNet architecture | 128: 128: 128: 128: 128: 128
Training data size | 4 $\times$ 5120$\times$ 202
Validation data size | 1 $\times$ 5120$\times$ 202
Testing data size | 1 $\times$ 5120$\times$ 202
Table 6: Various hyper-parameters used for the 1D KS equation Parameters | Values
---|---
Spatial dimension $n$ | 120
Temporal dimension $T$ | 5121
Time step $\Delta t$ | 0.05
Latent dimensions $z$ | 1, 2, 4, 6, 8, 10
AE architecture | 120: 120: 100: $z$: 100: 120: 120
No. of epochs for training AE models | 5000
Activation function | ReLU
Batch size | 32
Learning rate | $1e-03$
No. of epochs for training ResNet models | 20000
Step sizes for ResNet models |  $\Delta$t, 2$\Delta$t, 4$\Delta$t, 8$\Delta$t, 16$\Delta$t, 32$\Delta$t, 64$\Delta$t, 128$\Delta$t, 256$\Delta$t, 512$\Delta$t, 1024$\Delta$t
ResNet architecture | 1024: 1024: 1024
Training data size | 10 $\times$ 5121$\times$ 120
Validation data size | 5 $\times$ 5121$\times$ 120
Testing data size | 5 $\times$ 5121$\times$ 120