# Region of Interest Detection in Melanocytic Skin Tumor Whole Slide Images

Yi Cui
Department of Economics
University of North Carolina at Chapel Hill
Chapel Hill, NC 27514
<EMAIL_ADDRESS>
Yao Li
Department of Statistics & Operations Research
University of North Carolina at Chapel Hill
Chapel Hill, NC 27514
<EMAIL_ADDRESS>
Jayson R. Miedema
School of Medicine
University of North Carolina at Chapel Hill
Chapel Hill, NC 27514
<EMAIL_ADDRESS>
Sherif Farag
School of Medicine
University of North Carolina at Chapel Hill
Chapel Hill, NC 27514
<EMAIL_ADDRESS>
J.S. Marron
Department of Statistics & Operations Research
University of North Carolina at Chapel Hill
Chapel Hill, NC 27514
<EMAIL_ADDRESS>
Nancy E. Thomas
School of Medicine
University of North Carolina at Chapel Hill
Chapel Hill, NC 27514
<EMAIL_ADDRESS>
Corresponding author.

###### Abstract

Automated region of interest detection in histopathological image analysis is
a challenging and important topic with tremendous potential impact on clinical
practice. The deep-learning methods used in computational pathology help us to
reduce costs and increase the speed and accuracy of regions of interest
detection and cancer diagnosis. In this work, we propose a patch-based region
of interest detection method for melanocytic skin tumor whole-slide images. We
work with a dataset that contains 165 primary melanomas and nevi Hematoxylin
and Eosin whole-slide images and build a deep-learning method. The proposed
method performs well on a hold-out test data set including five TCGA-SKCM
slides (accuracy of 93.94% in slide classification task and intersection over
union rate of 41.27% in the region of interest detection task), showing the
outstanding performance of our model on melanocytic skin tumor. Even though we
test the experiments on the skin tumor dataset, our work could also be
extended to other medical image detection problems, such as various tumors’
classification and prediction, to help and benefit the clinical evaluation and
diagnosis of different tumors.

##  1 Introduction

From the study of Pollack et al. (2011), melanocytic skin tumor that has
metastasized (distant stage) or is thicker than 4.00 mm has a poor prognosis
(5-year survival: 15.7% and 56.6%). Also, Gerger et al. (2005); Argenziano et
al. (2010); Seidenari et al. (2012) show that established methods have been
proved insufficient to diagnose melanoma at an early stage. One solution may
be histopathological images, which have long been utilized in treatment
decisions and prognostics for cancer. Deep learning-based predictors trained
on annotated data and non-annotated data could be a potential efficient
technology to improve the early detection and diagnosis of such tumors and
overall survival rates of treatment.

The goal of this work is to develop a deep neural network-based Region of
Interest (ROI) detection method that can precisely detect the ROIs in
melanocytic skin tumors through Whole Slide Images (WSIs) and at the same
time, diagnose the type of skin tumor (melanoma vs. nevus) accurately.
However, most deep learning architectures cannot work on large WSIs directly
due to their large sizes. Recent studies (Lerousseau et al., 2020; Lu et al.,
2021) have shown some progress in dealing with this issue by multiple instance
learning. They propose to break down the large images into small patches and
extract features from these patches, then slide-level features are generated
by pooling the patch features together. Following that direction, we also
break the large WSIs into patches and train our model with them. Besides, we
leverage some partial information from annotations that are not made for ROI
detection to help improve the performance of the model. Also, we show that the
proposed method can work with various training set sizes by decreasing the
size of the training set to different levels and testing the model performance
on the same test set. Figure 1 illustrates the overview of our method.

Figure 1: Overview of the proposed detection framework. 1) Extract melanoma,
nevus and other patches from training data. 2) Train a three-class patch
classifier based on extracted patches. 3) For each slide, generate predicted
scores for all patches and calculate patch as well as slide classification
accuracy. 4) Rank all patches from a slide based on the corresponding
predicted scores in the context of melanoma or nevus, depending on the slide
classification result. 5) Generate visualization results based on predicted
scores.

##  2 Data and Methods

Data pre-processing. WSIs (165 in total) from private melanocytic tumor
dataset cohort and TCGA-SKCM cohort are randomly split and assigned to 80% for
training and 20% for testing. Given a WSI, patches are extracted based on the
slide-level label and annotations (see subfigure 1 of figure 1). If the slide-
level label is nevus, all patches inside the annotated regions will be labeled
as nevus. If the slide-level label is melanoma, all patches inside the
annotated regions will be labeled as melanoma. Besides patches from annotated
regions, some patches outside those regions are also extracted and labeled as
other. The reason why we add the class of other is that if we just train a two
class classifier model, it will wrongly classify some other patches as
melanoma or nevus and those wrongly classified patches will also affect the
slide classification result. However, since not all ROIs are annotated by
pathologists, there could be melanoma and nevus patches outside annotated
regions. To avoid labeling those patches as other, we manually extract patches
of other classes from regions. To show that our method can work with small
amount of training samples, we conduct the experiments with 80% (107 WSIs),
60% (80 WSIs), 40% (54 WSIs), and 20% (33 WSIs) of training samples, and test
the models on the same test set (31 WSIs).

Training data preparation and model. A three-class patch classification model
(PCLA-3C) is trained on the labeled patches with VGG16, see Simonyan and
Zisserman (2015), as base architecture. In the testing stage, all patches from
a WSI will first be fed into the trained patch classifier. Ignoring patches
predicted as other, slide-level prediction is done by majority vote based on
patches predicted as melanoma and nevus. For a WSI classified as melanoma, all
the patches from this slide will be ranked by melanoma predicted scores.
Otherwise, all the patches will be ranked by nevus predicted scores. Three
visualization maps are generated based on the predicted scores (See subfigure
5 of figure 1).

Model assessment. To evaluate the performance of ROI detection, the annotated
ratio is measured to calculate Intersection over Union (IoU) for each slide.
Given a slide, annotated ratio $\beta$ is calculated by the number of patches
in the annotated regions divided by the number of patches extracted from the
slide:

$\displaystyle\beta=\frac{A_{p}}{C_{p}},$

where $A_{p}$ is the number of patches in A (annotated region) and $C_{p}$ is
the number of patches in C (WSI). Then, the top $n\beta$ patches based on
predicted scores are classified as ROI, where $n$ is the total number of
patches from a slide. The performance will be measured by the IoU, which
compares the annotated regions and predicted ROI regions. Since the framework
is patch-based, IoU is calculated by the number of patches in the intersection
region (the region in both annotated and predicted regions) divided by the
number of patches in the union of the annotated and predicted ROI regions:

$\displaystyle\text{IoU}=\frac{\underline{AB}_{p}}{\overline{AB}_{p}},$

where $\underline{AB}_{p}$ shows the number of patches in the region of
$(A\cap B)$ and $\overline{AB}_{p}$ shows the number of patches in the region
of $(A\cup B)$. A is annotated region and B is the predicted/highlighted
region.

##  3 Results and discussion

Table 1: Robustness performance of patch classification accuracy, slide
classification and IoU by PCLA-3C and CLAM using different splits of the
original training set. Since CLAM does not do patch classification, it does
not have patch classification accuracy.

| 20% split | 40% split
---|---|---
| PCLA-3C | CLAM | PCLA-3C | CLAM
| Mean | 95% CI | Mean | 95% CI | Mean | 95% CI | Mean | 95% CI
Patch classification accuracy | 0.6397 | [0.5193, 0.7601] | - | - | 0.7887 | [0.7536, 0.8238] | - | -
Slide classification accuracy | 0.7406 | [0.6627, 0.8185] | 0.6710 | [0.6386, 0.7033] | 0.8430 | [0.8043, 0.8817] | 0.6976 | [0.6619, 0.7333]
Intersection over Union | 0.3026 | [0.2394, 0.3327] | 0.0427 | [0.0342, 0.0512] | 0.3402 | [0.3057, 0.3784] | 0.0524 | [0.0297, 0.0751]
| 60% split | 80% split
| PCLA-3C | CLAM | PCLA-3C | CLAM
| Mean | 95% CI | Mean | 95% CI | Mean | 95% CI | Mean | 95% CI
Patch classification accuracy | 0.8191 | [0.7766, 0.8616] | - | - | 0.8210 | [0.7949, 0.8471] | - | -
Slide classification accuracy | 0.8721 | [0.8458, 0.8985] | 0.7097 | [0.6830, 0.7364] | 0.8885 | [0.8607, 0.9163] | 0.7258 | [0.7117, 0.7399]
Intersection over Union | 0.3652 | [0.3369, 0.3934] | 0.0621 | [0.0428, 0.0814] | 0.3710 | [0.3335, 0.4084] | 0.1103 | [0.0529, 0.1677]

Two methods are tested on the melanocytic skin tumor dataset to finish ROI
detection and slide classification: 1) CLAM (clustering-constrained attention
multiple instance learning) from Lu et al. (2021), 2) PCLA-3C (the proposed
patch-based classification
model)111https://github.com/cyMichael/ROI_Detection. Both methods are trained
on the training set, and the performances on both training and testing sets
are evaluated. The relatively low IoUs come from the fact that ground-truth
labels are not available for evaluation. We can see more visualization results
in Appendix.

We use CLAM and PCLA-3C by randomly selecting different proportions of
original training data. We find that by utilizing 60% of original training
data, our accuracy is 0.8191 (95% CI, 0.7766-0.8616) at the patch level, and
accuracy is 0.8721 (95% CI, 0.8458-0.8985) at the slide level, see table 1 for
detailed values of each split. To calculate the 95% CI and mean value for 60%
training data, for instance, we randomly split the training data by 10 times
and then train the model on these independent subsamples. Our true testing
data are kept unchanged since these data include true annotations. Table 1’s
results are based on the testing set, and we also test whether ROIs can be
accurately detected in the testing set. As in the private melanocytic tumor
dataset cohort, our pathology team annotated the TCGA-SKCM cohorts to mark
ROIs. One figure in the appendix shows one representative sample of
annotations (black dot regions) as well as without annotations. There is a
high agreement between the predictions of the ROI and the PCLA-3C,
demonstrating the generality of our automatic ROI detection. We perform CI
estimation to enhance model performance. The unannotated and annotated models
achieve an accuracy of 0.7406 (95% CI: 0.6627-0.8185) at the slide level and
an accuracy of 0.6397 (95% CI: 0.5193-0.7601) at the patch level by only
utilizing 20% of original training samples. For the IoU value, we can see that
we have a mean of 37.10% for 80% of training data, which is much higher than
the results of CLAM. As in the PCLA-3C, the improvements in patch
classification accuracy, slide classification accuracy and IoU show the
importance of annotations in the training of deep learning classifiers for
prediction. We also test whether melanoma, nevus or other status classifiers
can directly predict certain tumor types. This is important as accurate tumor
type is the clinical biomarker for future treatment. In summary, our deep-
learning-based framework outperforms the state-of-the-art ROI detection
methods like Lu et al. (2021), which is quite crucial in medical imaging
fields and related treatment recommendations.

## References

  * Ankerst et al. (1999) Ankerst, M., Breunig, M.M., Kriegel, H.P., Sander, J., 1999\.  Optics: Ordering points to identify the clustering structure.  ACM Sigmod record 28, 49–60.
  * Argenziano et al. (2010) Argenziano, G., Kittler, H., Ferrara, G., Rubegni, P., Malvehy, J., Puig, S., Cowell, L., Stanganelli, I., De Giorgi, V., Thomas, L., Bahadoran, P., Menzies, S.W., Piccolo, D., Marghoob, A.A., Zalaudek, I., 2010.  Slow-growing melanoma: A dermoscopy follow-up study.  British Journal of Dermatology 162, 267–273.  doi:10.1111/j.1365-2133.2009.09416.x.
  * Gerger et al. (2005) Gerger, A., Koller, S., Kern, T., Massone, C., Steiger, K., Richtig, E., Kerl, H., Smolle, J., 2005\.  Diagnostic applicability of in vivo confocal laser scanning microscopy in melanocytic skin tumors.  Journal of Investigative Dermatology 124, 493–498.  URL: http://dx.doi.org/10.1111/j.0022-202X.2004.23569.x, doi:10.1111/j.0022-202X.2004.23569.x.
  * Lerousseau et al. (2020) Lerousseau, M., Vakalopoulou, M., Classe, M., Adam, J., Battistella, E., Carré, A., Estienne, T., Henry, T., Deutsch, E., Paragios, N., 2020\.  Weakly Supervised Multiple Instance Learning Histopathological Tumor Segmentation.  Lecture Notes in Computer Science (including subseries Lecture Notes in Artificial Intelligence and Lecture Notes in Bioinformatics) 12265 LNCS, 470–479.  doi:10.1007/978-3-030-59722-1{\\_}45.
  * Lu et al. (2021) Lu, M.Y., Williamson, D.F., Chen, T.Y., Chen, R.J., Barbieri, M., Mahmood, F., 2021\.  Data-efficient and weakly supervised computational pathology on whole-slide images.  Nature Biomedical Engineering 5, 555–570.  doi:10.1038/s41551-020-00682-w.
  * Pollack et al. (2011) Pollack, L.A., Li, J., Berkowitz, Z., Weir, H.K., Wu, X.C., Ajani, U.A., Ekwueme, D.U., Li, C., Pollack, B.P., 2011.  Melanoma survival in the United States, 1992 to 2005\.  Journal of the American Academy of Dermatology 65, S78.e1–S78.e10.  URL: http://dx.doi.org/10.1016/j.jaad.2011.05.030, doi:10.1016/j.jaad.2011.05.030.
  * Seidenari et al. (2012) Seidenari, S., Bassoli, S., Borsari, S., Ferrari, C., Giusti, F., Ponti, G., Tomasini, C., Magnoni, C., 2012\.  Variegated dermoscopy of in situ melanoma.  Dermatology 224, 262–270.  doi:10.1159/000338696.
  * Simonyan and Zisserman (2015) Simonyan, K., Zisserman, A., 2015\.  Very deep convolutional networks for large-scale image recognition.  3rd International Conference on Learning Representations, ICLR 2015 - Conference Track Proceedings , 1–14.

##  Appendix A Appendix

###  A.1 Potential Negative Societal Impact

To the best of our knowledge, we do not foresee many potential negative
societal impacts to our work. Furthermore, the deep learning-based methods
does not conduct any direct diagnosis but rather aids doctors with those ROI
regions.

###  A.2 Supplemental materials

Computational configuration. All analyses were used by Python. Images were
analyzed and processed using OpenSlide. All the computational tasks were
finished on Cluster with Linux (Tested on Ubuntu 18.04), NVIDIA GPU (Tested on
Nvidia GeForce RTX 3090 on local workstations). NVIDIA GPUs supports were
followed to set up and configure CUDA (Tested on CUDA 11.3) and the torch
version should be greater than or equal to 1.7.1.

Explanations for the figures. Figure 2 shows one representative sample of
annotations (black dot regions) as well as without annotations. Visualization
maps of two samples from the melanocytic skin tumor dataset are shown in
figure 3 and figure 4. Figure 3 shows the three types of figures of the
largest ROI region from melanoma found by the proposed method. Figure 4 is the
same for the nevus sample. The overlap map highlights top-ranked patches in a
WSI and masks other area with a transparent blue color. The percentage of
highlighted patches equals $\beta$ (the annotated ratio). The boundary map
shows the boundary of the largest ROI cluster based on the highlighted
patches, where the highlighted patches are clustered by OPTICS algorithm from
Ankerst et al. (1999). The last one is a heatmap where red covers regions that
have high predicted scores and blue covers regions that have low predicted
scores.

Figure 2: One representative sample of annotations (black dot regions) as well
as without annotations. Some ROIs on the right side are not annotated.

(a) Boundary of the largest predicted ROI region

(b) Overlay

(c) Heatmap

Figure 3: Visualization results for a melanoma sample.

(a) Boundary of the largest predicted ROI region

(b) Overlay

(c) Heatmap

Figure 4: Visualization results for a nevus sample.


# Energy-Preserving Reduced Operator Inference
for Efficient Design and Control

Tomoki Koike111PhD Candidate, School of Aerospace Engineering, AIAA Student
Member<EMAIL_ADDRESS>and Elizabeth Qian222Assistant Professor, Schools
of Aerospace Engineering & Computational Science and Engineering, AIAA Member,
<EMAIL_ADDRESS>Georgia Institute of Technology, Atlanta, Georgia, 30332

###### Abstract

Many-query computations, in which a computational model for an engineering
system must be evaluated many times, are crucial in design and control. For
systems governed by partial differential equations (PDEs), typical high-
fidelity numerical models are high-dimensional and too computationally
expensive for the many-query setting. Thus, efficient surrogate models are
required to enable low-cost computations in design and control. This work
presents a physics-preserving reduced model learning approach that targets
PDEs whose quadratic operators preserve energy, such as those arising in
governing equations in many fluids problems. The approach is based on the
Operator Inference method, which fits reduced model operators to state
snapshot and time derivative data in a least-squares sense. However, Operator
Inference does not generally learn a reduced quadratic operator with the
energy-preserving property of the original PDE. Thus, we propose a new energy-
preserving Operator Inference (EP-OpInf) approach, which imposes this
structure on the learned reduced model via constrained optimization. Numerical
results using the viscous Burgers’ and Kuramoto-Sivashinksy equation (KSE)
demonstrate that EP-OpInf learns efficient and accurate reduced models that
retain this energy-preserving structure.

##  1 Nomenclature

$n$ = |  full model state dimension
---|---
$r$ = |  reduced model state dimension
$\mathbf{A},\hat{\mathbf{A}}$ = |  full and reduced linear operators in $\mathbb{R}^{n\times n}$, $\mathbb{R}^{r\times r}$
$\mathbf{H},\hat{\mathbf{H}}$ = |  full and reduced quadratic operators in $\mathbb{R}^{n\times n^{2}}$, $\mathbb{R}^{r\times r^{2}}$
$\mathbf{F},\hat{\mathbf{F}}$ = |  full and reduced quadratic operators in $\mathbb{R}^{n\times n(n+1)/2}$, $\mathbb{R}^{r\times r(r+1)/2}$, containing unique columns of $\mathbf{H},\hat{\mathbf{H}}$
$K$ = |  total number of data points
$t$ = |  time variable
$\mathbb{R}^{m\times n}$ = |  Euclidean space with dimension $m$-by-$n$
$\mathbf{V}_{r}$ = |  POD basis in $\mathbb{R}^{n\times r}$
$\mathbf{x},\hat{\mathbf{x}}$ = |  full and reduced system states for linear-quadratic ODE in $\mathbb{R}^{n},~{}\mathbb{R}^{r}$
$\mathbf{X},\hat{\mathbf{X}}$ = |  full and reduced state snapshot data matrices in $\mathbb{R}^{n\times K},~{}\mathbb{R}^{r\times K}$
$\hat{\mathbf{X}}_{\otimes}$ = |  reduced state quadratic data matrix in $\mathbb{R}^{r^{2}\times K}$
$\dot{\mathbf{X}},\dot{\hat{\mathbf{X}}}$ = |  full and reduced state time derivative data matrix in $\mathbb{R}^{n\times K}$, $\mathbb{R}^{r\times K}$
$\omega$ = |  spatial variable
$\otimes$ = |  Kronecker product
$\langle\cdot,\cdot\rangle$ = |  Euclidean inner product
$\|\cdot\|,\|\cdot\|_{F}$ = |  Euclidean and Frobenius norm
$\text{vech}(\cdot)$ = |  half-vectorization
$\text{Cov}[\cdot,\cdot]$ = |  covariance
$\text{Var}[\cdot]$ = |  variance

##  2 Introduction

Many-query computations, which require repeated evaluations of a predictive
model, arise in many different engineering contexts, including inverse
problems, state estimation, uncertainty propagation, design, optimization, and
control. In these many-query settings, high-dimensional models, such as those
emerging from the spatial discretization of a PDE, can be prohibitively
expensive to evaluate. Projection-based model reduction addresses this
challenge by deriving an approximate model with a smaller dimension by
projecting the original model governing equations onto a low-dimensional
subspace [1], thereby reducing computational cost. Model reduction methods
have enabled fast and cheap predictive simulations in many practical
applications such as turbulent aerodynamic flow, weather forecasting, and
large-scale structural analysis [2, 3, 4].

Model reduction methods are divided into intrusive and non-intrusive
approaches. Intrusive methods are applicable when there is direct access to
both the model and its underlying computational code; this access is used to
intrusively compute reduced operators of the projected governing equations.
Intrusive strategies include system theoretic methods like balanced truncation
[5] and rational interpolation methods [6, 7, 8], as well as sampling-based
methods such as Proper Orthogonal Decomposition (POD) [9, 10, 11], Balanced
Proper Orthogonal Decomposition [12], and reduced basis methods [13, 14, 15].
Because intrusive reduced models are constructed by projecting the original
governing equations onto a reduced subspace, intrusive methods can preserve
some types of mathematical structure in the original model (stability,
symmetries, invariances, etc.) that reflect physical principles governing the
original dynamics.

However, in many real-world systems, we encounter cases where we can compute
state trajectories from given initial conditions but lack access to the full
model operators or the code itself. This limitation necessitates the
development of non-intrusive model reduction methods which derive the reduced
model operators from data rather than via intrusive projection. Non-intrusive
methods which infer reduced matrix operators include Dynamic Mode
Decomposition [16], Operator Inference (OpInf) [17, 18, 19], and Lift & Learn
[20]. Several works have also developed approaches that use deep neural
networks to learn low-dimensional dynamics [21, 22, 23, 24]. An alternative
family of methods uses sparse regression to identify from data the terms of a
system’s governing equations [25, 26]. These approaches, leveraging available
data and some understanding of system dynamics and structures, make surrogate
modeling feasible for problems where the intrusive access to the model code is
unavailable. However, because non-intrusive reduced model operators are not
obtained from direct projection of the full model operators, they often fail
to retain mathematical structure that the intrusive operators would inherit.
This can yield reduced models that may lack accuracy, robustness, and
generalizability [27]. There is therefore a need to develop non-intrusive
model reduction methods that preserve these types of structure.

Many structure-preserving reduced model learning methods have been developed,
each incorporating different physical structures into their model learning
process in different ways. Among these is a physics-informed dynamic mode
decomposition method introduced in [27], which learns linear operators via a
manifold optimization to preserve certain structures including circulant and
symmetric structures. However, this method is applicable only to linear
operators. For Hamiltonian systems, the works [28, 29] constrain the Operator
Inference approach to retain Hamiltonian structure. Further contributing to
this field, [30] developed a structure-preserving OpInf method that ensures
the inferred system operators are Lagrangian and can preserve symmetry or
symmetric positive definite structure. These structure-preserving non-
intrusive methods have been demonstrated on applications ranging from the
linear wave equation to Schrödinger’s equation to soft robots [27, 28, 29,
30].

In this paper, we consider quadratic systems in which the quadratic term of
the governing equations is energy-preserving: this means that the contribution
of the quadratic term to the time derivative of the system energy is zero
[31]. Several works have already incorporated this energy-preserving structure
into different non-intrusive model learning methods. For instance, [32]
modified the sparse regression approach [25, 26] to promote learning quadratic
terms with this structure [33]. The work [34] proposes an OpInf method for
enforcing this energy-preserving structure on the quadratic term in port-
Hamiltonian systems by formulating a constrained OpInf optimization where the
constraint exploits the port-Hamiltonian structure. In contrast, in this
paper, we propose an OpInf method that enforces this quadratic energy-
preserving structure that is applicable to non-port-Hamiltonian systems.

We introduce Energy-Preserving OpInf (EP-OpInf), an OpInf method that enforces
energy-preserving structure in the learned quadratic operator via a
constrained optimization. Numerical experiments are conducted for the viscous
Burgers’ equation and the Kuramoto-Sivashinsky equation (KSE). Numerical
results demonstrate that EP-OpInf learned models achieve accuracy similar to
that of the unconstrained data-fit model learned via standard OpInf while non-
intrusively retaining the energy-preserving structure that intrusive models
inherit. This lays the foundation for applying stability analyses developed
for intrusive reduced models to the EP-OpInf learned models.

The remainder of the paper is organized as follows. Section 3 introduces the
quadratic systems that we consider and provides an overview of both intrusive
model reduction via POD and the standard non-intrusive OpInf method. Section 4
then presents our proposed EP-OpInf method and discusses some implementation
considerations. Section 5 demonstrates the efficacy of our proposed method
through a series of numerical experiments. Conclusions are presented in
Section 6.

##  3 Background

In this section, we introduce the quadratic system with energy-preserving
quadratic term that we consider (Section 3.1), the POD method for model
reduction (Section 3.2), and the standard OpInf method for learning reduced
models non-intrusively (Section 3.3).

###  3.1 Energy Preserving Quadratic Nonlinearities in PDEs

We consider an $n$-dimensional ordinary differential equation (ODE) which is
linear and quadratic in the state $\mathbf{x}$. Such an ODE often arises from
spatially discretizing a PDE and is given by

$\dot{\mathbf{x}}(t)=\mathbf{A}\mathbf{x}(t)+\mathbf{H}\left(\mathbf{x}(t)\otimes\mathbf{x}(t)\right)$
(1)

where $\mathbf{x}(t)\in\mathbb{R}^{n}$ is the system state vector over
$t\in[0,T_{\text{final}}]$, and $\otimes$ denotes the Kronecker product. The
operators $\mathbf{A}\in\mathbb{R}^{n\times n}$ and
$\mathbf{H}\in\mathbb{R}^{n\times n^{2}}$ are the linear and quadratic
operators, respectively. In our setting, $n$ is large, so simulating the
system (1) is computationally expensive.

The quadratic operator $\mathbf{H}$ is called ‘energy-preserving’ if

$\langle\mathbf{x},\mathbf{H}(\mathbf{x}\otimes\mathbf{x})\rangle=\mathbf{x}^{\top}\mathbf{H}(\mathbf{x}\otimes\mathbf{x})=0,\qquad\text{
for all }\mathbf{x}\in\mathbb{R}^{n}.$ (2)

This condition is derived by setting the quadratic term in the time derivative
of the energy, $\frac{1}{2}\|\mathbf{x}\|^{2}$, to zero. Note that this
condition can be generalized to energies defined by weighted norms or inner
products; our EP-OpInf approach also would generalize to such settings
following the weighted OpInf approach of [18]. Such energy-preserving
quadratic terms commonly arise from nonlinear convection terms in the
governing equations for fluid problems, including the Navier-Stokes equations.
Some types of stability guarantees are possible for models with energy-
preserving quadratic nonlinearities [31], making it desirable to enforce this
property in learning reduced models that will be used for design and control.

###  3.2 Proper Orthogonal Decomposition

Proper orthogonal decomposition [9, 10, 11] originated from the analysis of
turbulent flows in aerodynamics, and it has become one of the most widespread
projection-based model reduction methods. POD reduces the model by projecting
it onto a reduced subspace defined to be the span of basis vectors that
optimally represent a set of simulation or experimental data.

In POD, we begin by collecting snapshots of state trajectory time series data
by simulating the original full model ODE (1) with $K$ timesteps. We define
the state snapshot data matrix as follows:

$\boldsymbol{\mathbf{X}}=\begin{bmatrix}|&|&&|\\\
\boldsymbol{\mathbf{x}}(t_{1})&\boldsymbol{\mathbf{x}}(t_{2})&\cdots&\boldsymbol{\mathbf{x}}(t_{K})\\\
|&|&&|\\\ \end{bmatrix}\in\mathbb{R}^{n\times K}~{}.$ (3)

More generally, the state snapshot matrix can contain state data from multiple
simulations, e.g., from different initial conditions or using different
parameters. Let $\mathbf{X}=\mathbf{V\Sigma W}^{\top}$ denote the singular
value decomposition of the state snapshot matrix (3). To reduce the dimension
of the large-scale model, we denote by $\mathbf{V}_{r}\in\mathbb{R}^{n\times
r}$ the first $r\ll n$ columns of $\mathbf{V}$; this is called the POD basis.
Then, we approximate the state $\mathbf{x}$ in the subspace spanned by the POD
basis, $\mathbf{x}\approx\mathbf{V}_{r}\hat{\mathbf{x}}$ where
$\hat{\mathbf{x}}\in\mathbb{R}^{r}$ is called the reduced state. If we
substitute this approximation into (1) and enforce the Galerkin orthogonality
condition that the approximation residual be orthogonal to the span of
$\mathbf{V}_{r}$, we arrive at a POD-Galerkin reduced model of the form

$\dot{\hat{\mathbf{x}}}(t)=\hat{\mathbf{A}}\hat{\mathbf{x}}(t)+\hat{\mathbf{H}}(\hat{\mathbf{x}}(t)\otimes\hat{\mathbf{x}}(t)),$
(4)

where the reduced operators are
$\hat{\mathbf{A}}=\mathbf{V}^{\top}_{r}\mathbf{AV}_{r}\in\mathbb{R}^{r\times
r}$ and
$\hat{\mathbf{H}}=\mathbf{V}^{\top}_{r}\mathbf{H}(\mathbf{V}_{r}\otimes\mathbf{V}_{r})\in\mathbb{R}^{r\times
r^{2}}$. Since the size of the reduced operators scales with the reduced
dimension $r$ and not the original dimension $n$, we are able to evaluate the
reduced model efficiently in many-query computations arising in design and
control. The computation of the reduced operators in traditional POD-Galerkin
model reduction is intrusive, that is, we require access to the full operators
$\mathbf{A}$ and $\mathbf{H}$ in order to compute their projections. The
energy-preserving property of the original matrix $\mathbf{H}$ (2) is
inherited in the reduced matrix $\hat{\mathbf{H}}$ through this intrusive
projection process.

However, cases exist where access to the original operators is limited or
impossible, but we know the model structure, such as linear, quadratic, etc.,
based on the application. In such situations, Operator Inference [17],
described in Section 3.3, can non-intrusively learn reduced operators solely
from state trajectory and input data.

###  3.3 Operator Inference

The goal of Operator Inference is to non-intrusively obtain a reduced model of
the form (4). To do so, we will fit reduced operators $\hat{\mathbf{A}}$ and
$\hat{\mathbf{H}}$ to the reduced data in a least-squares sense. In addition
to the state trajectory data (3), we also require paired state time derivative
data:

$\displaystyle\dot{\mathbf{X}}=\begin{bmatrix}|&|&&|\\\
\dot{\mathbf{x}}(t_{1})&\dot{\mathbf{x}}(t_{2})&\cdots&\dot{\mathbf{x}}(t_{K})\\\
|&|&&|\end{bmatrix}\in\mathbb{R}^{n\times K}\quad\text{ where
~{}~{}}\dot{\mathbf{x}}(t_{i})=\mathbf{A}\mathbf{x}(t_{i})+\mathbf{H}(\mathbf{x}(t_{i})\otimes\mathbf{x}(t_{i}))~{}.$

This time derivative data can come directly from the simulation of the high-
dimensional ODE or can be approximated numerically from the snapshot data. We
use the POD basis $\mathbf{V}_{r}$ to compute reduced state and time
derivative data as follows: let
$\hat{\mathbf{x}}_{i}=\mathbf{V}_{r}^{\top}\mathbf{x}(t_{i})$, and
$\dot{\hat{\mathbf{x}}}_{i}=\mathbf{V}_{r}^{\top}\dot{\mathbf{x}}(t_{i})$ for
$i=1,\ldots,K$. Then, define

$\hat{\textbf{X}}=\begin{bmatrix}|&|&&|\\\
\boldsymbol{\hat{\mathbf{x}}}_{1}&\boldsymbol{\hat{\mathbf{x}}}_{2}&\cdots&\boldsymbol{\hat{\mathbf{x}}}_{K}\\\
|&|&&|\\\ \end{bmatrix}\in\mathbb{R}^{r\times
K},\qquad\text{and}\qquad\dot{\hat{\textbf{X}}}=\begin{bmatrix}|&|&&|\\\
\dot{\hat{\mathbf{x}}}_{1}&\dot{\hat{\mathbf{x}}}_{2}&\cdots&\dot{\hat{\mathbf{x}}}_{K}\\\
|&|&&|\\\ \end{bmatrix}\in\mathbb{R}^{r\times K}.$ (5)

Additionally, we define the matrix formed by the quadratic terms of the state
data

$\hat{\textbf{X}}_{\otimes}=\begin{bmatrix}|&|&&|\\\
(\hat{\mathbf{x}}_{1}\otimes\hat{\mathbf{x}}_{1})&(\hat{\mathbf{x}}_{2}\otimes\hat{\mathbf{x}}_{2})&\cdots&(\hat{\mathbf{x}}_{K}\otimes\hat{\mathbf{x}}_{K})\\\
|&|&&|\\\ \end{bmatrix}\in\mathbb{R}^{r^{2}\times K}.$ (6)

This allows us to formulate the following minimization for finding the reduced
operators $\hat{\mathbf{A}}$ and $\hat{\mathbf{H}}$:

$\textbf{Standard OpInf}:\qquad\min_{\hat{\mathbf{A}}\in\mathbb{R}^{r\times
r},~{}\hat{\mathbf{H}}\in\mathbb{R}^{r\times
r^{2}}}\sum_{i=1}^{K}\left\|\dot{\hat{\mathbf{x}}}_{i}-\hat{\mathbf{A}}\hat{\mathbf{x}}_{i}-\hat{\mathbf{H}}(\hat{\mathbf{x}}_{i}\otimes\hat{\mathbf{x}}_{i})\right\|^{2}_{2}=\min_{\mathbf{O}\in\mathbb{R}^{r\times(r+r^{2})}}\|\mathbf{D}\mathbf{O}^{\top}-\dot{\hat{\mathbf{X}}}\|_{F}^{2},$
(7)

where
$\mathbf{D}=[\hat{\mathbf{X}}^{\top},~{}\hat{\mathbf{X}}_{\otimes}^{\top}]\in\mathbb{R}^{K\times(r+r^{2})}$
and
$\mathbf{O}=[\hat{\mathbf{A}},~{}\hat{\mathbf{H}}]\in\mathbb{R}^{r\times(r+r^{2})}$.
We note that the redundancy in the Kronecker product means that the minimizer
$\hat{\mathbf{H}}$ of (7) is non-unique; this non-uniqueness can be eliminated
by requiring its entries to satisfy a symmetry condition (see Section 4.2 for
a discussion on implementation). The optimization (7) is a linear least-
squares problem and can be solved efficiently using standard numerical linear
algebra packages. In fact, the problem (7) decouples into $r$ independent
least-squares problems, one for each row of the reduced operators with
$r+r^{2}$ degrees of freedom each.

##  4 Energy-Preserving Operator Inference

Section 4.1 introduces our Energy-Preserving Operator Inference method.
Section 4.2 discusses efficient implementation of the method.

###  4.1 Energy-Preserving Constraint

The Operator Inference approach in the preceding section finds operators that
best fit the reduced data. However, it does not necessarily find a quadratic
operator $\mathbf{H}$ with the energy-preserving property (2). Thus, we now
introduce Energy-Preserving Operator Inference (EP-OpInf), which learns a
quadratic reduced model with an energy-preserving quadratic operator.

We first re-write the energy-preservation property (2) as a constraint on the
entries of $\hat{\mathbf{H}}$. Let $\hat{\mathbf{H}}[l,m]$ correspond to an
element of $\hat{\mathbf{H}}$ at the row index $l$ and column index $m$, and
define

$\hat{h}_{ijk}=\hat{\mathbf{H}}[i,r(k-1)+j].$ (8)

We assume without loss of generality that the entries of $\hat{\mathbf{H}}$
are symmetric in the last two indices, i.e., $\hat{h}_{ijk}=\hat{h}_{ikj}$ (in
Section 4.2 we discuss how we enforce this in implementation). Then, the
energy-preserving condition can be rewritten in terms of the entries of
$\hat{\mathbf{H}}$ as [32, 31]:

$\hat{h}_{ijk}+\hat{h}_{jik}+\hat{h}_{kij}=0,\quad 1\leq i,j,k\leq r~{}.$ (9)

Due to the symmetry of the last two indices of $\hat{h}_{ijk}$, this defines a
total of $\frac{r(r+1)(r+2)}{6}$ unique constraints. We emphasize that
submatrices of size $r^{\prime}<r$ extracted from the $\hat{\mathbf{H}}$
corresponding to $i,j,k=1,\ldots,r^{\prime}$ also satisfy the energy-
preserving condition, since the constraint 9 is true for all index
combinations.

To impose this energy-preserving structure on the operator, we propose EP-
OpInf. For this method, we incorporate the constraint (9) into optimization
(7) and formulate a constrained minimization as follows:

$\textbf{EP-
OpInf}:\qquad\min_{\mathbf{O}\in\mathbb{R}^{r\times(r+r^{2})}}\|\mathbf{D}\mathbf{O}^{\top}-\dot{\hat{\mathbf{X}}}\|_{F}^{2}\quad\text{
subject to }\quad\hat{h}_{ijk}+\hat{h}_{jik}+\hat{h}_{kij}=0,\quad 1\leq
i,j,k\leq r~{}.$ (10)

We note that the original Operator Inference approach (7) was an unconstrained
linear least-squares problem that could be efficiently decoupled into
individual rows of the operators and solved with standard numerical linear
algebra packages. With the addition of our constraints, the minimization can
no longer be decoupled into rows, and the minimization must be solved using
constrained optimization solvers. There is thus a trade-off between the
simplicity of implementing the Operator Inference approach and preserving the
energy-preserving structure of the quadratic term.

###  4.2 Implementation Considerations

For the efficient computation of the optimization (10), we reformulate the
quadratic operator and associated data matrix constructed by the Kronecker
product. The code implementations for this reformulation are available
online333GitHub repository: https://github.com/smallpondtom/LiftAndLearn.jl.
We note that the matrix $\hat{\mathbf{X}}_{\otimes}$ (6) contains redundant
rows due to the symmetry of the Kronecker product. Thus, the solution of the
minimization (10) is not unique unless we enforce symmetry in the last two
indices of the learned operator $\hat{\mathbf{H}}$. To do so, we will re-
express the optimization in terms of an operator
$\hat{\mathbf{F}}\in\mathbb{R}^{r\times r(r+1)/2}$, which removes the
redundant columns of $\hat{\mathbf{H}}$. This has the added advantage of
reducing the dimension of the optimization. This operator $\hat{\mathbf{F}}$
is defined as follows. Let $\hat{f}_{ijk}$ denote the entry of the $i$-th row
in column $(r-k/2)(k-1)+j$, that is:

$\hat{f}_{ijk}=\hat{\mathbf{F}}[i,(r-k/2)(k-1)+j]~{},\quad\text{for
}~{}~{}i=1,\ldots,r;~{}j=1,\ldots,r;~{}k=1,\ldots,j.$ (11)

Then, the entries of $\hat{\mathbf{H}}$ and $\hat{\mathbf{F}}$ can be related
as follows:

$\hat{f}_{ijk}=\begin{cases}\hat{h}_{ijk}&\textrm{ if }j=k\\\
\hat{h}_{ijk}+\hat{h}_{ikj}&\textrm{ otherwise }\end{cases},\quad
i=1,\ldots,r;~{}j=1,\ldots,r;~{}k=1,\ldots,j~{}.$ (12)

With this relation and considering the permutations of indices, we rewrite the
quadratic energy-preserving constraint as

$\delta_{jk}\hat{f}_{ijk}+\delta_{ik}\hat{f}_{jik}+\delta_{ij}\hat{f}_{kij}=0~{},\quad\textrm{where}\quad\delta_{jk}=\begin{cases}1&\textrm{
if }~{}j=k\\\ 1/2&\textrm{ otherwise }\end{cases},\quad
i=1,\ldots,r;~{}j=1,\ldots,i;~{}k=1,\ldots,j~{}.$ (13)

Denote $\mathbf{x}^{[2]}=\text{vech}(\mathbf{xx}^{\top})$ (‘half-
vectorization’ of $\mathbf{xx}^{\top}$), to be obtained by eliminating the
supradiagonal elements of $\mathbf{xx}^{\top}$ and vectorizing it, i.e.,
$x^{[2]}_{(r-k/2)(k-1)+j}=x_{j}x_{k}$ [35]. Using this operation, let

$\hat{\textbf{X}}^{[2]}=\begin{bmatrix}|&|&&|\\\
\hat{\mathbf{x}}_{1}^{[2]}&\hat{\mathbf{x}}_{2}^{[2]}&\cdots&\hat{\mathbf{x}}_{K}^{[2]}\\\
|&|&&|\\\ \end{bmatrix}\in\mathbb{R}^{(r(r+1)/2)\times K}.$ (14)

Then, we can reformulate the optimization (10) with matrices
$\mathbf{\widetilde{D}}=[\hat{\mathbf{X}}^{\top},\hat{\mathbf{X}}^{[2]\top}]\in\mathbb{R}^{r\times(r+r(r+1)/2)}$
and
$\mathbf{\widetilde{O}}=[\hat{\mathbf{A}},\hat{\mathbf{F}}]\in\mathbb{R}^{r\times(r+r(r+1)/2)}$
as

$\min_{{\mathbf{\widetilde{O}}}\in\mathbb{R}^{r\times(r+r(r+1)/2)}}\|\mathbf{\widetilde{D}}\mathbf{\widetilde{O}^{\top}}-\mathbf{R}\|^{2}_{F}\quad\text{
subject to
}\quad\delta_{jk}\hat{f}_{ijk}+\delta_{ik}\hat{f}_{jik}+\delta_{ij}\hat{f}_{kij}=0~{},\quad
i=1,\ldots,r;~{}j=1,\ldots,i;~{}k=1,\ldots,j~{}.$ (15)

This reformulation decreases the degrees of freedom in the least squares
problem associated with the quadratic operator from $r^{2}$ to $r(r+1)/2$, and
also reduces the number of constraints in the minimization problem, making it
more efficient than a naïve implementation of (10).

##  5 Numerical Experiments

This section demonstrates the proposed EP-OpInf method through its application
to the viscous Burgers’ equation (Section 5.1) and the Kuramoto-Sivashinksy
equation (Section 5.2).

###  5.1 Viscous Burgers’ Equation

In our first numerical experiment, we consider the viscous Burgers’ equation,
a simplified fluid dynamics model that describes the motion of a one-
dimensional viscous fluid [36]. It is given as:

$\frac{\partial x}{\partial
t}(\omega,t)=\mu\frac{\partial^{2}x}{\partial\omega^{2}}(\omega,t)-x(\omega,t)\frac{\partial
x}{\partial\omega}(\omega,t)~{},$ (16)

where $x(\omega,t)$ is the velocity of the fluid, $t$ is time, $\omega$ is the
spatial coordinate, and $\mu$ is the kinematic viscosity coefficient. The
viscous Burgers’ equation can be derived from the Navier-Stokes equations
under certain assumptions, including the assumption of one-dimensional flow
and the neglect of pressure gradient driving the flow. The Burgers’ equation
is suitable for testing our EP-OpInf method since it is expressed in the form
of (1) where the energy-preserving quadratic nonlinearity holds [37].

Figure 1: Visualization of Burgers’ equation state evolution and flow field
with $\mu=0.1$ and reduced dimension of 10 with initial condition
$x(\omega,0)=0.8\sin(2\pi\omega-0.125)$.

####  5.1.1 Setup

We solve Burgers’ equation on the unit domain with periodic boundaries,
$\omega\in[0,1)$, for $t\in[0,T]$, and for $\mu=0.1$. We use a sinusoidal
initial condition $x(\omega,0)=A\sin(2\pi f\omega+\phi)$ with three varied
parameters: the amplitude $A$, frequency $f$, and phase $\phi$. We collect
training data at $A=[0.8,0.9,1.0,1.1,1.2]$, $f=[1,2,3]$, and
$\phi=[-0.25,-0.125,0,0.125,0.25]$. In total, $5\times 3\times 5=75$ different
parameterized initial conditions define the training data set.

This system is discretized on a uniform spatial grid with a grid spacing
$\Delta\omega=2^{-7}$, leading to an ODE of the form (1). With this grid size,
the state dimension is $n=128$. The ODE is solved using a semi-implicit Euler
scheme with a time step size of $\Delta t=10^{-4}$ and a final time of $T=1$.
Figure 1 shows the behavior of the state for sample parameter values.

In the model learning phase, data at every 100th timestep are collected from
the simulation described above, and these data are used to construct the data
snapshot matrices and train an inferred reduced model. To assess the trained
operators, we conduct two tests where we generate additional data with
different initial conditions:

  * •

Test 1 (Interpolation): 50 different initial conditions generated uniformly
within the training region:
$A\in[0.8,1.2]$, $f\in\\{1,2,3\\}$, and $\phi\in[-0.25,0.25]$

  * •

Test 2 (Extrapolation): 50 different initial conditions generated uniformly
outside the training region:
$A\in[0.5,0.8)\cup(1.2,1.5]$, $f\in\\{4,5,6\\}$, and
$\phi\in[-0.5,-0.25)\cup(0.25,0.5]$.

To determine the number of POD modes to retain, we compute the relative energy
lost by a POD basis of size $r$, defined to be

$1-\frac{\sum_{i=r+1}^{n}\sigma_{i}^{2}}{\sum_{i=1}^{n}\sigma_{i}^{2}}~{}.$
(17)

In this equation, $\sigma_{i}$ represents the singular values of the snapshot
data matrix. The energy spectrum for the Burgers’ equation in Figure 2 shows
that when $r_{\text{max}}=10$ we capture 99.99999% of the energy in the data.
We therefore set $r_{\text{max}}=10$ in our numerical experiments. For both
standard OpInf and EP-OpInf, the reduced operators for dimensions
$r<r_{\text{max}}$ are retrieved by extracting the appropriate submatrices
from $r_{\text{max}}$-dimensional model operators. As mentioned in Section
4.1, the energy-preserving structure holds even for extracted operators.

Figure 2: Relative energy lost from truncated modes for viscous Burgers’
equation training data.

Once we train the inferred models, we use the learned models to reconstruct
the training data by simulating them from the same initial condition that
generated the training data. We report the relative state error of the reduced
model reconstruction, defined as

$\text{average relative state
error}=\frac{1}{J}\sum_{j=1}^{J}\frac{\left\|\mathbf{X}_{j}-\overline{\mathbf{X}}_{j}\mathbf{V}_{r}^{\top}\right\|^{2}_{F}}{\left\|\mathbf{X}_{j}\right\|_{F}^{2}}~{},$
(18)

where $J$ is the total number of initial conditions and
$\overline{\mathbf{X}}$ contains the trajectory data obtained by integrating
the surrogate model defined by either the intrusive (4) or inferred operators
with either standard OpInf (7) or EP-OpInf (10).

Additionally, because a key goal of this work is to learn models whose
quadratic terms have the energy-preserving property, we also report the
following measure of the violation of the energy-preserving property:

$\text{energy-preserving constraint violation}=\sum_{1\leq i,j,k\leq
r}|\hat{h}_{ijk}+\hat{h}_{jik}+\hat{h}_{kji}|~{},$ (19)

which is defined to be the sum of constraint residuals.

####  5.1.2 Results

The left plot in Figure 3 shows the relative state error for the training data
computed from the intrusive, standard OpInf, and EP-OpInf models. Figure 4
plots the energy-preserving constraint violations vs. the size of the reduced
model. From Figure 3, we see the accuracy of all methods in reconstructing the
training trajectories is similar. From the constraint violations in Figure 4,
it is evident that the intrusive and EP-OpInf models satisfy the energy-
preserving condition up to machine precision, indicating that the quadratic
operators of the two models possess the energy-preserving property. We
emphasize that EP-OpInf does this non-intrusively, whereas the intrusive
method requires access to the original model code. Additionally, despite the
learning process for EP-OpInf being constrained by the energy-preserving
quadratic structure, EP-OpInf achieves a level of model accuracy that is very
close to both that of the unconstrained standard OpInf method and the
intrusive method.

The right plot in Figure 3 plots the test errors for the intrusive, standard
OpInf, and EP-OpInf models, for both the test data sets considered. For
different initial conditions generated within the training region, the
accuracy is comparable to that of the training data results shown in the plot
on the left, with all three methods yielding similar levels of accuracy. Even
if we extrapolate from the training region as in Test 2, the standard OpInf
and EP-OpInf errors are on par with the intrusive method. We emphasize that
both the standard and EP-OpInf models non-intrusively achieve accuracy
comparable to that of the intrusive model, and that EP-OpInf does this while
enforcing the energy-preserving structure.

Figure 3: Comparison of intrusive, standard OpInf, and EP-OpInf methods of
relative state error for training data 18 (left) and test cases (right).
Figure 4: Constraint violation of all three methods for Burgers’ equation over
reduced dimensions $r$.

###  5.2 Kuramoto-Sivashinsky Equation

Originally introduced in the context of flame propagation, the Kuramoto-
Sivashinsky equation has since found applications in various fields, including
fluid dynamics, pattern formation, and combustion theory [38, 39]. The
equation is given by:

$\frac{\partial x}{\partial
t}(\omega,t)+\frac{\partial^{2}x}{\partial\omega^{2}}(\omega,t)+\mu\frac{\partial^{4}x}{\partial\omega^{4}}(\omega,t)+x(\omega,t)\frac{\partial
x}{\partial\omega}(\omega,t)=0~{}.$ (20)

Here, $x(\omega,t)$ represents the fluid velocity, $t$ denotes time, $\omega$
is the spatial coordinate, and $\mu$ is the kinematic viscosity coefficient.
Like the viscous Burgers’ equation, the KSE satisfies the energy-preserving
quadratic nonlinearity [37]. However, the KSE poses challenges for model
reduction. This is because the KSE models localized turbulence in which
coherent spatial structures and temporal chaos coexist [40]. The temporal
chaos indicates the exponential separations of nearby state trajectories of
the KSE solutions; thereby, even small errors in the initial condition can be
amplified, making the long-time dynamics unpredictable [41]. This
unpredictability increases the error when reconstructing the state
trajectories from the reduced models. In contrast, the coherent spatial
structures indicate inertial manifolds, i.e., invariant manifolds containing
global attractors for dissipative dynamical systems [42, 43, 44]. The
solutions of KSE are attracted asymptotically to this manifold, and it is
important for OpInf to capture this finite-dimensional manifold accurately.

Several studies in the literature have explored model reduction approaches for
the KSE. For instance, one study employed an intrusive Galerkin projection
method with regularization to derive a reduced model and compared the state
trajectories to the direct numerical simulations [45]. Other studies
implemented non-intrusive methods such as OpInf and NARMAX (Nonlinear
AutoRegressive Moving Average with eXogenous input) [46] to explore
forecasting capabilities with reduced models [47, 43]. In contrast, our work
implements a physical structure-preserving non-intrusive model learning
method.

####  5.2.1 Setup

We solve the KSE on the spatial domain $\omega\in[0,L)$ and the time domain
$t\in[0,T]$, with periodic boundaries. We set $\mu=1$, and consider the
initial condition $x(\omega,0)=a\cos(2\pi\omega/L)+b\cos(4\pi\omega/L)$, where
$a$ and $b$ are parameters. For our experiment, we chose the spatial domain
size of $L=22$, which, according to [44], is a sufficiently large domain size
for the KSE’s dynamics to transition to an inertial manifold characterized as
a chaotic attractor [40]. Considering the increase in error for the intrusive
model beyond $t>150$, as illustrated in Figure 6, we chose $T=300$ for the
time domain to adequately capture the progression of dynamics from the
transient phase into the chaotic regime.

Furthermore, we discretize this PDE on a uniform spatial grid of $n=512$ grid
points, yielding an ODE in the form of (1). To integrate the ODE, we implement
a semi-implicit scheme where we use the Crank-Nicolson scheme for the linear
term and the Adams-Bashforth scheme for the quadratic term, maintaining a
constant time step size of $\Delta t=10^{-3}$. We collect training data at
$a=[0.8,1.0,1.2]$ and $b=[0.2,0.4,0.6]$ for a total of 9 initial conditions,
and retain the state snapshot at every 100th timestep in the training data
set.

Figure 5 plots the energy spectrum (17) of the KSE data set. We note that
$r=9$ is required to retain 99% of the snapshot energy, while $r=24$ retains
99.999999%. Therefore, for both standard OpInf and EP-OpInf, the models are
learned using $r_{\text{max}}=24$, and we extract submatrices of sizes $r>9$
and $r<r_{max}$. We note that the POD basis sizes required to retain most of
the data energy are notably higher than those for the Burgers’ equation. This
is expected due to the chaotic/turbulent nature of the KSE.

After learning the models, we test them on two different test data sets not
seen in training:

  * •

Test 1 (Interpolation): 50 uniformly sampled initial conditions within the
training parameter region: $a\in[0.8,1.2]$ and $b\in[0.2,0.6]$

  * •

Test 2 (Extrapolation): 50 uniformly sampled initial conditions outside the
training parameter region: $a\in[0,0.8)\cup(1.2,2.0]$ and
$b\in[0,0.2)\cup(0.6,0.8]$.

Figure 5: Relative energy lost from truncated modes for Kuramoto-Sivashinksy
equation training data. Figure 6: Training data flow field plots of the full
and reduced models for an initial condition with $a=0.8,~{}b=0.2$ and a
reduced dimension of $r=24$. The first column shows the state flow fields of
the full model and reconstructed state flow fields for each method. The second
column illustrates the errors of flow fields compared to the full model
states.

####  5.2.2 Measures of reduced model performance

The chaotic nature of the KSE dynamics means that small perturbations can lead
to large differences in the state at a later time. Therefore, instead of
assessing the reduced model performance using pointwise error metrics, we
instead consider the reduced model predictions of the system statistics. We
follow the work [43] which uses the state autocorrelation function as a
statistical measure for assessing the performance of surrogate models. The
normalized autocorrelation function is defined as follows:

$\rho(k)=\frac{\mathrm{Cov}[x(t),x(t+k)]}{\sqrt{\mathrm{Var}[x(t)]\mathrm{Var}[x(t+k)]}}~{},$
(21)

where $k$ is the time lag and the covariance and variance are taken over
infinite time [48]. The autocorrelation quickly decays to zero as the time lag
increases, a phenomenon known as correlation splitting which indicates the
statistical independence of the system state at a given time $t$ from its past
and future trajectories and is a key indicator of chaos [49, 50, 51, 52, 53].

In the case of finite time series data of a stationary process, we use the
sample normalized autocorrelation function

$\hat{\rho}(k)=\frac{c_{k}}{c_{0}}\quad\text{ where
~{}~{}}c_{k}=\frac{1}{T}\sum_{t=1}^{T-k}\left(x(t)-\bar{x}\right)\left(x(t+k)-\bar{x}\right)~{}.$
(22)

Here, $\bar{x}$ is the sample mean or time-averaged state values [48]. To
quantitatively evaluate our reduced models’ performance, we introduce the
normalized autocorrelation error, averaged over $J$ initial conditions:

$\text{average normalized autocorrelation
error}=\frac{1}{J}\sum_{j=1}^{J}\frac{\|\hat{\rho}_{j}(k)-\bar{\hat{\rho}}_{j}(k)\|^{2}_{2}}{\|{\hat{\rho}}_{j}(k)\|^{2}_{2}},$
(23)

where $\hat{\rho}_{j}(k)$ and $\bar{\hat{\rho}}_{j}(k)$ represent the sample
normalized autocorrelation of the full and reduced models, respectively.

We thus report the autocorrelation (21) and autocorrelation error (23) for
both the learned and intrusive models in our results. Additionally, to
qualitatively assess the model predictions, we present flow field plots
alongside their corresponding errors from the full model.

Figure 7: Normalized autocorrelation averaged over all training initial
conditions. The states of each normalized autocorrelation function are
reconstructed with reduced dimensions of $r=9,12,20,24$.

####  5.2.3 Results

In the first column of Figure 6, the flow fields predicted by the full model
are presented alongside those reconstructed by the intrusive, standard OpInf,
and EP-OpInf models with a reduced dimension of $r=24$. The second column
plots the errors between the full and reduced model flow fields. We observe
that all four models yield qualitatively similar predictions. During the
transient time period $t\in[0,100]$, the errors are small relative to the rest
of the domain, demonstrating that all three methods can accurately predict
transient trajectories in the pre-chaotic regime. However, the errors begin to
grow as the system transitions into the chaotic regime beyond this interval.
During the interval $t\in[100,200]$, the EP-OpInf method exhibits smaller
errors than the standard OpInf, and its performance is on par with the
intrusive methods. Figure 7 illustrates the average normalized autocorrelation
as a function of time lag, each at different reduced model dimensions
$r=9,12,20,24$, together with the autocorrelation of the full model. Each
column represents the results for the training, Test 1 (interpolation), and
Test 2 (extrapolation) data. For all three results, we observe that as the
model dimension $r$ increases, the reduced models converge toward the full
model’s autocorrelation pattern and better capture the correlation splitting.
The OpInf method generally aligns with the full model’s trend but exhibits
some discrepancies, particularly at lower dimensions, while the EP-OpInf is
closer to the intrusive method. For Tests 1 and 2, we see that all three
methods have smaller discrepancies among one another compared to the training
data. This observation could be attributed to the increased number (50) of
trajectories for which the error is averaged in the test cases in contrast to
the 9 trajectories available for training.

The energy-preserving constraint violation plot in Figure 8 validates that EP-
OpInf allows the quadratic operator to satisfy the energy-preserving
structure. Furthermore, we observe a decreasing trend in the autocorrelation
error plots of Figure 9, indicating the autocorrelation functions of the
reduced models converging to that of the full model. Despite the three methods
showing comparable error results, the error margins for Test 2 are slightly
larger than those for Test 1 since the initial conditions are sampled outside
the training parameter regimes.

Finally, Figure 10 provides a visual comparison of flow fields and their
corresponding errors for our two distinct tests. In Test 1, where the initial
conditions are uniformly drawn from the same parameter space as the training
set, the state reconstructions by EP-OpInf are on par with those by the
standard OpInf method. The intrusive method, on the other hand, demonstrates
very small errors. Notably, the errors in this test for standard OpInf and EP-
OpInf do not show a large discrepancy from the training results demonstrated
in Figure 6, indicating consistent model performance within the training
parameter space. In contrast, for Test 2, which involves parameters outside
the training regime, OpInf and EP-OpInf exhibit larger errors across an
expanded domain, but the errors of the two non-intrusive approaches are
comparable to each other. Overall, EP-OpInf results demonstrate our proposed
method’s capability of capturing the statistical properties and model accuracy
comparable to the standard OpInf and intrusive approaches while satisfying the
energy-preserving constraint and being a non-intrusive method.

Figure 8: Energy-preserving constraint violation (19) of the KSE. Figure 9:
Normalized autocorrelation error results of the training and two tests.
Training and tests are averaged over 9 and 50 initial conditions,
respectively.

Figure 10: Flow field plots generated using one of the uniformly generated
initial condition parameters for each Test 1 and Test 2. The parameters are
(Test 1) $a=1.1141,~{}b=0.5135$ and (Test 2) $a=0.5055,~{}b=0.0308$, and the
reduced dimension is $r=24$. The first column shows the flow fields of the
full model and reconstructed state flow fields for each reduced-order method.
The second column illustrates the errors of the flow fields compared to the
full model states.

##  6 Conclusion

In this paper, we introduced Energy-Preserving Operator Inference, a new
reduced model learning method for quadratic dynamical systems whose quadratic
term preserves energy. Such systems arise naturally in fluid dynamics and this
structure enables certain types of stability guarantees. Learning models from
data which preserve this structure is therefore desirable when the learned
models will be used for design and control. The EP-OpInf method incorporates
an equality constraint (9) into the model learning problem to learn a reduced
quadratic operator which has energy-preserving structure.

We conducted numerical experiments on model problems drawn from fluid
dynamics, namely the viscous Burgers’ and Kuramoto-Sivashinksy (KSE)
equations. In these experiments, EP-OpInf achieved similar accuracy to the
standard OpInf approach as well as to the intrusive reduced model. We
emphasize that the EP-OpInf and standard OpInf approaches do this non-
intrusively, which is an advantage in settings where intrusive access to the
model code is not possible. Additionally, while the standard OpInf method
fails to learn energy-preserving quadratic operators, our new EP-OpInf method
succeeds at non-intrusively enforcing this structure. Directions for future
work include enforcing other types of structure within the OpInf framework as
well as building on the current EP-OpInf method to develop stability
guarantees for learned models with energy-preserving quadratic terms.

## References

  * Benner et al. [2015] Benner, P., Gugercin, S., and Willcox, K., “A Survey of Projection-Based Model Reduction Methods for Parametric Dynamical Systems,” _SIAM Review_ , Vol. 57, No. 4, 2015, p. 483–531.  10.1137/130932715, URL http://dx.doi.org/10.1137/130932715.
  * Rowley and Dawson [2017] Rowley, C. W., and Dawson, S. T., “Model Reduction for Flow Analysis and Control,” _Annual Review of Fluid Mechanics_ , Vol. 49, No. 1, 2017, p. 387–417.  10.1146/annurev-fluid-010816-060042, URL http://dx.doi.org/10.1146/annurev-fluid-010816-060042.
  * Horenko et al. [2008] Horenko, I., Klein, R., Dolaptchiev, S., and Schütte, C., “Automated Generation of Reduced Stochastic Weather Models I: Simultaneous Dimension and Model Reduction for Time Series Analysis,” _Multiscale Modeling & Simulation_, Vol. 6, No. 4, 2008, p. 1125–1145.  10.1137/060670535, URL http://dx.doi.org/10.1137/060670535.
  * Besselink et al. [2013] Besselink, B., Tabak, U., Lutowska, A., van de Wouw, N., Nijmeijer, H., Rixen, D., Hochstenbach, M., and Schilders, W., “A comparison of model reduction techniques from structural dynamics, numerical mathematics and systems and control,” _Journal of Sound and Vibration_ , Vol. 332, No. 19, 2013, p. 4403–4422.  10.1016/j.jsv.2013.03.025, URL http://dx.doi.org/10.1016/j.jsv.2013.03.025.
  * Moore [1981] Moore, B., “Principal component analysis in linear systems: Controllability, observability, and model reduction,” _IEEE Transactions on Automatic Control_ , Vol. 26, No. 1, 1981, p. 17–32.  10.1109/tac.1981.1102568, URL http://dx.doi.org/10.1109/TAC.1981.1102568.
  * Grimme [1997] Grimme, E. J., _Krylov projection methods for model reduction_ , University of Illinois at Urbana-Champaign, 1997.
  * Antoulas [2009] Antoulas, A. C., _Approximation of large-Scale Dynamical Systems_ , SIAM, 2009.
  * Antoulas et al. [2010] Antoulas, A. C., Beattie, C. A., and Gugercin, S., _Interpolatory Model Reduction of Large-Scale Dynamical Systems_ , Springer US, 2010, p. 3–58.  10.1007/978-1-4419-5757-3_1, URL http://dx.doi.org/10.1007/978-1-4419-5757-3_1.
  * Lumley [1967] Lumley, J. L., “The structure of inhomogeneous turbulent flows,” _Atmospheric turbulence and radio wave propagation_ , 1967, pp. 166–178.
  * Sirovich [1987] Sirovich, L., “Turbulence and the dynamics of coherent structures. I. Coherent structures,” _Quarterly of applied mathematics_ , Vol. 45, No. 3, 1987, pp. 561–571.
  * Berkooz et al. [1993] Berkooz, G., Holmes, P., and Lumley, J. L., “The proper orthogonal decomposition in the analysis of turbulent flows,” _Annual review of fluid mechanics_ , Vol. 25, No. 1, 1993, pp. 539–575.
  * Willcox and Peraire [2002] Willcox, K., and Peraire, J., “Balanced Model Reduction via the Proper Orthogonal Decomposition,” _AIAA Journal_ , Vol. 40, No. 11, 2002, p. 2323–2330.  10.2514/2.1570, URL http://dx.doi.org/10.2514/2.1570.
  * Peterson [1989] Peterson, J. S., “The Reduced Basis Method for Incompressible Viscous Flow Calculations,” _SIAM Journal on Scientific and Statistical Computing_ , Vol. 10, No. 4, 1989, p. 777–786.  10.1137/0910047, URL http://dx.doi.org/10.1137/0910047.
  * Veroy and Patera [2005] Veroy, K., and Patera, A. T., “Certified real-time solution of the parametrized steady incompressible Navier–Stokes equations: rigorous reduced-basis a posteriori error bounds,” _International Journal for Numerical Methods in Fluids_ , Vol. 47, No. 8–9, 2005, p. 773–788.  10.1002/fld.867, URL http://dx.doi.org/10.1002/fld.867.
  * Hesthaven et al. [2016] Hesthaven, J. S., Rozza, G., Stamm, B., et al., _Certified reduced basis methods for parametrized partial differential equations_ , Vol. 590, Springer, 2016.
  * Schmid [2010] Schmid, P. J., “Dynamic mode decomposition of numerical and experimental data,” _Journal of Fluid Mechanics_ , Vol. 656, 2010, p. 5–28.  10.1017/s0022112010001217, URL http://dx.doi.org/10.1017/S0022112010001217.
  * Peherstorfer and Willcox [2016] Peherstorfer, B., and Willcox, K., “Data-driven operator inference for nonintrusive projection-based model reduction,” _Computer Methods in Applied Mechanics and Engineering_ , Vol. 306, 2016, pp. 196–215.
  * Qian et al. [2022] Qian, E., Farcaş, I.-G., and Willcox, K., “Reduced Operator Inference for Nonlinear Partial Differential Equations,” _SIAM Journal on Scientific Computing_ , Vol. 44, No. 4, 2022, p. A1934–A1959.  10.1137/21m1393972, URL http://dx.doi.org/10.1137/21M1393972.
  * McQuarrie et al. [2021] McQuarrie, S. A., Huang, C., and Willcox, K. E., “Data-driven reduced-order models via regularised Operator Inference for a single-injector combustion process,” _Journal of the Royal Society of New Zealand_ , Vol. 51, No. 2, 2021, pp. 194–211.  10.1080/03036758.2020.1863237, URL https://doi.org/10.1080/03036758.2020.1863237, publisher: Taylor & Francis _eprint: https://doi.org/10.1080/03036758.2020.1863237.
  * Qian et al. [2020] Qian, E., Kramer, B., Peherstorfer, B., and Willcox, K., “Lift & Learn: Physics-informed machine learning for large-scale nonlinear dynamical systems,” _Physica D: Nonlinear Phenomena_ , Vol. 406, 2020, p. 132401.  10.1016/j.physd.2020.132401, URL https://linkinghub.elsevier.com/retrieve/pii/S0167278919307651.
  * Swischuk et al. [2019] Swischuk, R., Mainini, L., Peherstorfer, B., and Willcox, K., “Projection-based model reduction: Formulations for physics-based machine learning,” _Computers & Fluids_, Vol. 179, 2019, p. 704–717.  10.1016/j.compfluid.2018.07.021, URL http://dx.doi.org/10.1016/j.compfluid.2018.07.021.
  * Bhattacharya et al. [2021] Bhattacharya, K., Hosseini, B., Kovachki, N. B., and Stuart, A. M., “Model Reduction And Neural Networks For Parametric PDEs,” _The SMAI journal of computational mathematics_ , Vol. 7, 2021, p. 121–157.  10.5802/smai-jcm.74, URL http://dx.doi.org/10.5802/smai-jcm.74.
  * Hesthaven and Ubbiali [2018] Hesthaven, J., and Ubbiali, S., “Non-intrusive reduced order modeling of nonlinear problems using neural networks,” _Journal of Computational Physics_ , Vol. 363, 2018, p. 55–78.  10.1016/j.jcp.2018.02.037, URL http://dx.doi.org/10.1016/j.jcp.2018.02.037.
  * Gao et al. [2020] Gao, H., Wang, J.-X., and Zahr, M. J., “Non-intrusive model reduction of large-scale, nonlinear dynamical systems using deep learning,” _Physica D: Nonlinear Phenomena_ , Vol. 412, 2020, p. 132614.  10.1016/j.physd.2020.132614, URL http://dx.doi.org/10.1016/j.physd.2020.132614.
  * Brunton et al. [2016] Brunton, S. L., Proctor, J. L., and Kutz, J. N., “Discovering governing equations from data by sparse identification of nonlinear dynamical systems,” _Proceedings of the National Academy of Sciences_ , Vol. 113, No. 15, 2016, p. 3932–3937.  10.1073/pnas.1517384113, URL http://dx.doi.org/10.1073/pnas.1517384113.
  * Rudy et al. [2017] Rudy, S. H., Brunton, S. L., Proctor, J. L., and Kutz, J. N., “Data-driven discovery of partial differential equations,” _Science Advances_ , Vol. 3, No. 4, 2017.  10.1126/sciadv.1602614, URL http://dx.doi.org/10.1126/sciadv.1602614.
  * Baddoo et al. [2023] Baddoo, P. J., Herrmann, B., McKeon, B. J., Nathan Kutz, J., and Brunton, S. L., “Physics-informed dynamic mode decomposition,” _Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences_ , Vol. 479, No. 2271, 2023, p. 20220576.  10.1098/rspa.2022.0576, URL https://royalsocietypublishing.org/doi/full/10.1098/rspa.2022.0576, publisher: Royal Society.
  * Sharma et al. [2022] Sharma, H., Wang, Z., and Kramer, B., “Hamiltonian operator inference: Physics-preserving learning of reduced-order models for canonical Hamiltonian systems,” _Physica D: Nonlinear Phenomena_ , Vol. 431, 2022, p. 133122.  10.1016/j.physd.2021.133122, URL https://www.sciencedirect.com/science/article/pii/S0167278921002682.
  * Gruber and Tezaur [2023] Gruber, A., and Tezaur, I., “Canonical and noncanonical Hamiltonian operator inference,” _Computer Methods in Applied Mechanics and Engineering_ , Vol. 416, 2023, p. 116334.  10.1016/j.cma.2023.116334, URL http://dx.doi.org/10.1016/j.cma.2023.116334.
  * Sharma and Kramer [2022] Sharma, H., and Kramer, B., “Preserving Lagrangian structure in data-driven reduced-order modeling of large-scale dynamical systems,” , Nov. 2022.  10.48550/arXiv.2203.06361, URL http://arxiv.org/abs/2203.06361, arXiv:2203.06361 [cs, math].
  * Schlegel and Noack [2015] Schlegel, M., and Noack, B. R., “On long-term boundedness of Galerkin models,” _Journal of Fluid Mechanics_ , Vol. 765, 2015, p. 325–352.  10.1017/jfm.2014.736.
  * Kaptanoglu et al. [2021a] Kaptanoglu, A. A., Callaham, J. L., Aravkin, A., Hansen, C. J., and Brunton, S. L., “Promoting global stability in data-driven models of quadratic nonlinear dynamics,” _Phys. Rev. Fluids_ , Vol. 6, 2021a, p. 094401.  10.1103/PhysRevFluids.6.094401, URL https://link.aps.org/doi/10.1103/PhysRevFluids.6.094401.
  * Kaptanoglu et al. [2021b] Kaptanoglu, A. A., Morgan, K. D., Hansen, C. J., and Brunton, S. L., “Physics-constrained, low-dimensional models for magnetohydrodynamics: First-principles and data-driven approaches,” _Physical Review E_ , Vol. 104, No. 1, 2021b.  10.1103/physreve.104.015206, URL http://dx.doi.org/10.1103/PhysRevE.104.015206.
  * Goyal et al. [2023] Goyal, P., Duff, I. P., and Benner, P., “Guaranteed Stable Quadratic Models and their applications in SINDy and Operator Inference,” , 2023.  10.48550/ARXIV.2308.13819, URL https://arxiv.org/abs/2308.13819.
  * Magnus and Neudecker [1980] Magnus, J. R., and Neudecker, H., “The Elimination Matrix: Some Lemmas and Applications,” _SIAM Journal on Algebraic Discrete Methods_ , Vol. 1, No. 4, 1980, pp. 422–449.  10.1137/0601049, URL https://doi.org/10.1137/0601049.
  * Burgers [1948] Burgers, J. M., “A mathematical model illustrating the theory of turbulence,” _Advances in applied mechanics_ , Vol. 1, 1948, pp. 171–199.
  * Aref and Daripa [1984] Aref, H., and Daripa, P. K., “Note on Finite Difference Approximations to Burgers’ Equation,” _SIAM Journal on Scientific and Statistical Computing_ , Vol. 5, No. 4, 1984, pp. 856–864.  10.1137/0905060, URL http://epubs.siam.org/doi/10.1137/0905060.
  * Kuramoto [1978] Kuramoto, Y., “Diffusion-induced chaos in reaction systems,” _Progress of Theoretical Physics Supplement_ , Vol. 64, 1978, pp. 346–367.
  * Sivashinsky [1977] Sivashinsky, G. I., “Nonlinear analysis of hydrodynamic instability in laminar flames—I. Derivation of basic equations,” _Acta astronautica_ , Vol. 4, No. 11-12, 1977, pp. 1207–1221.
  * Hyman et al. [1986] Hyman, J. M., Nicolaenko, B., and Zaleski, S., “Order and complexity in the Kuramoto-Sivashinsky model of weakly turbulent interfaces,” _Physica D: Nonlinear Phenomena_ , Vol. 23, No. 1, 1986, pp. 265–292.  10.1016/0167-2789(86)90136-3, URL https://www.sciencedirect.com/science/article/pii/0167278986901363.
  * Strogatz [2018] Strogatz, S. H., _Nonlinear Dynamics and Chaos_ , CRC Press, 2018.  10.1201/9780429492563, URL http://dx.doi.org/10.1201/9780429492563.
  * Constantin et al. [1989] Constantin, P., Foias, C., Nicolaenko, B., and Teman, R., _Integral Manifolds and Inertial Manifolds for Dissipative Partial Differential Equations_ , Springer New York, 1989.  10.1007/978-1-4612-3506-4, URL http://dx.doi.org/10.1007/978-1-4612-3506-4.
  * Lu et al. [2017] Lu, F., Lin, K. K., and Chorin, A. J., “Data-based stochastic model reduction for the Kuramoto–Sivashinsky equation,” _Physica D: Nonlinear Phenomena_ , Vol. 340, 2017, pp. 46–57.  10.1016/j.physd.2016.09.007, URL https://doi.org/10.1016/j.physd.2016.09.007.
  * Cvitanović et al. [2010] Cvitanović, P., Davidchack, R. L., and Siminos, E., “On the State Space Geometry of the Kuramoto–Sivashinsky Flow in a Periodic Domain,” _SIAM Journal on Applied Dynamical Systems_ , Vol. 9, No. 1, 2010, pp. 1–33.  10.1137/070705623, URL https://epubs.siam.org/doi/10.1137/070705623, publisher: Society for Industrial and Applied Mathematics.
  * Sabetghadam and Jafarpour [2012] Sabetghadam, F., and Jafarpour, A., “$\alpha$ Regularization of the POD-Galerkin dynamical systems of the Kuramoto–Sivashinsky equation,” _Applied Mathematics and Computation_ , Vol. 218, No. 10, 2012, pp. 6012–6026.  10.1016/j.amc.2011.11.083, URL https://www.sciencedirect.com/science/article/pii/S0096300311014299.
  * Chen and Billings [1989] Chen, S., and Billings, S. A., “Representations of non-linear systems: the NARMAX model,” _International Journal of Control_ , Vol. 49, No. 3, 1989, p. 1013–1032.  10.1080/00207178908559683, URL http://dx.doi.org/10.1080/00207178908559683.
  * Almeida et al. [2022] Almeida, J. L. S., Pires, A. C., Cid, K. F. V., and Nogueira, A. C., “Non-intrusive operator inference for chaotic systems,” _IEEE Transactions on Artificial Intelligence_ , 2022, pp. 1–14.  10.1109/TAI.2022.3207449, URL https://ieeexplore.ieee.org/abstract/document/9893912, conference Name: IEEE Transactions on Artificial Intelligence.
  * Box et al. [1994] Box, G. E. P., Jenkins, G. M., and Reinsel, G. C., _Time Series Analysis: forcasting and control_ , 3rd ed., Pearson, Upper Saddle River, NJ, 1994.
  * Claus [1968] Claus, S., “Billingsley, P.: Ergodic Theory and Information. John Wiley & Sons, Inc., New York 1965. XIII + 193 S. Abb., Tab. Preis 64 s.” _Biometrische Zeitschrift_ , Vol. 10, No. 1, 1968, pp. 84–85.  https://doi.org/10.1002/bimj.19680100113, URL https://onlinelibrary.wiley.com/doi/abs/10.1002/bimj.19680100113.
  * Cornfeld et al. [1982] Cornfeld, I. P., Fomin, S. V., and Sinai, Y. G., _Ergodic Theory_ , Springer New York, 1982.  10.1007/978-1-4615-6927-5, URL https://doi.org/10.1007/978-1-4615-6927-5.
  * Zaslavskii [1985] Zaslavskii, G. M., _Chaos in Dynamic Systems_ , Harwood Academic, Amsterdam, Netherlands, 1985.
  * Sinai [1981] Sinai, Y., “The Stochasticity of Dynamical Systems,” Vol. 1, No. 1, 1981, pp. 100–119.
  * Anishchenko et al. [2003] Anishchenko, V. S., Vadivasova, T. E., Okrokvertskhov, G. A., and Strelkova, G. I., “Correlation analysis of dynamical chaos,” _Physica A: Statistical Mechanics and its Applications_ , Vol. 325, No. 1, 2003, pp. 199–212.  10.1016/S0378-4371(03)00199-7, URL https://www.sciencedirect.com/science/article/pii/S0378437103001997.