# Tensor-based flow reconstruction from optimally located sensor measurements

Mohammad Farazmand and Arvind K. Saibaba  Department of Mathematics, North
Carolina State University, 2311 Stinson Dr., Raleigh, NC 27695

###### Abstract

Reconstructing high-resolution flow fields from sparse measurements is a major
challenge in fluid dynamics. Existing methods often vectorize the flow by
stacking different spatial directions on top of each other, hence confounding
the information encoded in different dimensions. Here, we introduce a tensor-
based sensor placement and flow reconstruction method which retains and
exploits the inherent multidimensionality of the flow. We derive estimates for
the flow reconstruction error, storage requirements and computational cost of
our method. We show, with examples, that our tensor-based method is
significantly more accurate than similar vectorized methods. Furthermore, the
variance of the error is smaller when using our tensor-based method. While the
computational cost of our method is comparable to similar vectorized methods,
it reduces the storage cost by several orders of magnitude. The reduced
storage cost becomes even more pronounced as the dimension of the flow
increases. We demonstrate the efficacy of our method on three examples: a
chaotic Kolmogorov flow, in-situ and satellite measurements of the global sea
surface temperature, and 3D unsteady simulated flow around a marine research
vessel.

##  1 Introduction

Numerical simulations of fluid flows are carried out with ever growing spatial
resolution. In contrast, observational data is limited to relatively coarse
sensor measurements. This dichotomy inhibits efficient integration of
experimental data with existing high-fidelity computational methods to enable
detailed and accurate flow analysis or prediction.

As we review in §1.1, several methods have been developed to address this
disconnect. In particular, flow reconstruction methods seek to leverage
offline high-resolution simulations to estimate the entire flow field from
coarse online observations. All existing methods vectorize the simulation data
by stacking different spatial dimensions on top of each other. Vectorization
is convenient since it enables one to use familiar linear algebra techniques.
However, this approach inevitably leads to loss of information encoded in the
inherent multidimensional structure of the flow.

Here, we propose a tensor-based sensor placement and flow reconstruction
method which retains and exploits the multidimensional structure of the flow.
Our method significantly increases the accuracy of the reconstruction compared
to similar vectorized methods. We quantify this accuracy by deriving an upper
bound for the reconstruction error and show that the resulting approximation
exactly interpolates the flow at the sensor locations (assuming the
measurements are noise-free). Additionally, the proposed method has a smaller
memory footprint compared to the similar vectorized methods, and is scalable
to large datasets using randomized techniques. We demonstrate the efficacy of
our method on three examples: direct numerical simulation of a turbulent
Kolmogorov flow, in situ and satellite sea surface temperature (SST) data, and
3D unsteady simulated flow around a marine research vessel. We emphasize that
our focus here is only on flow reconstruction and not temporal prediction;
nonetheless, the reconstructed flow can be subsequently used as input to high-
fidelity or reduced-order predictive models.

###  1.1 Related work

Due to the broad applications of flow estimation from sparse measurements,
there is an expansive body of work on this subject (see Callaham et al. (2019)
for a thorough review). Here, we focus on the so-called library-based methods.
These methods seek to reconstruct the flow field by leveraging the sparse
observational measurements to interpolate a pre-computed data library
comprising high-fidelity numerical simulations. More specifically, consider a
scalar quantity $g(\boldsymbol{x},t)$ which we would like to reconstruct from
its spatially sparse measurements. For instance, this quantity may be a
velocity component, a vorticity component, pressure, or temperature. The data
library is a matrix $\Phi\in\mathbb{R}^{N\times T}$ whose columns are formed
from vectorized high-resolution simulations. Here, $N$ denotes the number of
collocation points used in the simulations. The columns of $\Phi$ may coincide
with the quantity of interest $g$ or be derived from this quantity, e.g.,
through proper orthogonal decomposition (POD) or dynamic mode decomposition
(DMD). The observational data $\boldsymbol{y}\in\mathbb{R}^{r}$ is a vector
containing $r$ measurements of the quantity $g$ at a particular time. Library-
based methods seek to find a map $F:\mathbb{R}^{N\times
T}\times\mathbb{R}^{r}\to\mathbb{R}^{N}$ such that $\boldsymbol{g}\simeq
F(\Phi,\boldsymbol{y})$. Here, $\boldsymbol{g}\in\mathbb{R}^{N}$ is a vector
obtained by stacking the quantity of interest $g(\boldsymbol{x},t)$ at the
collocation points.

Library-based methods differ in their choice of the data matrix $\Phi$ and the
methodology for finding the map $F$. A common choice for the columns of the
data matrix is the POD modes (Bui-Thanh et al., 2004; Willcox, 2006), although
DMD modes (Kramer et al., 2017; Dang et al., 2021) and flow snapshots (Clark
et al., 2021) have also been used. Bui-Thanh et al. (2004) use the gappy POD
algorithm of Everson & Sirovich (1995) to reconstruct the flow. They obtain
the map $F$ by solving a least squares problem which seeks to minimize the
discrepancy between the observations $\boldsymbol{y}$ and the reconstructed
flow $F(\Phi,\boldsymbol{y})$ at the sensor locations (also see Willcox
(2006)).

Discrete Empirical Interpolation Method (DEIM) takes a similar approach, but
the map $F$ is a suitable oblique projection on the linear subspace spanned by
the columns of $\Phi$. DEIM was first developed by Chaturantabut & Sorensen
(2010) for efficient reduced-order modeling of nonlinear systems and was later
used for flow reconstruction (Drmač & Gugercin, 2016; Wang et al., 2021).
Several subsequent modifications to DEIM have been proposed, e.g., to lower
its computational cost (Peherstorfer et al., 2014) and to generalize it for
use with weighted inner products (Drmač & Saibaba, 2018).

It is well-known that both gappy POD and DEIM suffer from overfitting
(Peherstorfer et al., 2020). Consequently, if the sensor data $\boldsymbol{y}$
is corrupted by significant observational noise, the reconstruction error will
be large. Callaham et al. (2019) use sparsity promoting techniques from image
recognition to overcome this problem (also see Chu & Farazmand (2021)). Their
reconstruction map $F$ is obtained by solving a sparsity-promoting
optimization problem with the constraint that the reconstruction error is
below a prescribed threshold. The resulting method is robust to observational
noise. However, unlike gappy POD and DEIM, the reconstruction map $F$ cannot
be expressed explicitly in terms of the training data $\Phi$.

Yet another flow reconstruction method is to represent the map $F$ with a
neural network which is trained using the observations $\boldsymbol{y}$ and
the data matrix $\Phi$. For instance, Yu & Hesthaven (2019) use an autoencoder
to represent the reconstruction map $F$ (also see Fukami et al. (2019);
Carlberg et al. (2019); Erichson et al. (2020)). Unlike DEIM, where the
reconstruction map $F$ is a linear combination of modes, neural networks can
construct nonlinear maps from the observations $\boldsymbol{y}$ and the
library $\Phi$. These machine learning methods have shown great promise;
however, the resulting reconstructions are not explicit, or even
interpretable, since they are only available as a complex neural network.

With the notable exception of convolutional neural networks (Fukami et al.,
2019; Carlberg et al., 2019), almost all existing methods treat the data as a
vector by stacking different spatial dimensions on top of each other. This
inevitably leads to loss of information encoded in the inherent
multidimensionality of the flow. Here, we propose a tensor-based method which
retains and exploits this multidimensional structure. Our method is similar to
the tensor-based DEIM which was recently proposed by Kirsten (2022) for model
reduction; but we use it for flow reconstructions which is the focus of this
paper. Numerical experiments show that the resulting reconstructions are more
accurate compared to similar vectorized methods, because of our method’s
ability to capture and exploit the inherent multidimensional nature of the
data. The computational cost of our tensor-based method is comparable to the
vectorized methods and can be further accelerated using randomized methods.
Furthermore, the tensor-based method requires much less storage compared to
the vectorized methods. This is especially important in large-scale 3D flows
where the storage costs can be substantial (Gelss et al., 2019). Although our
method is a tensorized version of DEIM, a similar tensor-based approach can be
applied to other flow reconstruction methods such as gappy POD, sparsity-
promoting methods, and autoencoders.

##  2 Tensor-based flow reconstruction

###  2.1 Set-up and preliminaries

Let $g(\boldsymbol{x},t)$ denote the quantity of interest at time $t$ that we
would like to reconstruct. This quantity may for instance be a velocity
component, a vorticity component, pressure, or temperature. The spatial
variable is denoted by $\boldsymbol{x}\in\Omega\subset\mathbb{R}^{d}$, where
$\Omega$ is the flow domain with $d=2$ or $d=3$ for two- and three-dimensional
flows, respectively. In numerical simulations, the quantity of interest $g$ is
discretized on a spatial grid of size $N_{1}\times N_{2}\times\cdots\times
N_{d}$ and saved as a tensor
$\boldsymbol{\mathcal{G}}\in\mathbb{R}^{N_{1}\times\cdots\times N_{d}}$.

We denote entries of $\boldsymbol{\mathcal{G}}$ by $g_{i_{1},\dots,i_{d}}$
where $1\leq i_{n}\leq N_{n}$ and $1\leq n\leq d$. There are $d$ different
matrix unfoldings of $\boldsymbol{\mathcal{G}}$, also called matricizations,
which we denote by
$\boldsymbol{G}_{(n)}\in\mathbb{R}^{N_{n}\times(\prod_{j\neq n}N_{j})}$. The
mode-$n$ product of a tensor $\boldsymbol{\mathcal{G}}$ with a matrix
$\boldsymbol{M}\in\mathbb{R}^{R\times N_{n}}$ is denoted as
$\boldsymbol{\mathcal{Y}}=\boldsymbol{\mathcal{G}}\times_{n}\boldsymbol{M}$
with entries

$y_{i_{i},\dots,j,\dots,i_{N}}=\sum_{k=1}^{R}g_{i_{1},\dots,k,\dots,i_{N}}m_{jk}\qquad
1\leq j\leq R.$ (1)

In terms of matrix unfoldings, it can be expressed as
$\boldsymbol{Y}_{(n)}=\boldsymbol{M}\boldsymbol{G}_{(n)}$. For matrices
$\boldsymbol{A}$ and $\boldsymbol{B}$ of compatible dimensions
$\boldsymbol{\mathcal{G}}\times_{m}\boldsymbol{A}\times_{n}\boldsymbol{B}=\boldsymbol{\mathcal{G}}\times_{n}\boldsymbol{B}\times_{m}\boldsymbol{A}$
if $m\neq n$ and
$\boldsymbol{\mathcal{G}}\times_{n}\boldsymbol{A}\times_{n}\boldsymbol{B}=\boldsymbol{\mathcal{G}}\times_{n}\boldsymbol{BA}$.
Associated with every tensor $\boldsymbol{\mathcal{G}}$ is a multirank
$(R_{1},\dots,R_{d})$ where $R_{n}=\text{rank}(\boldsymbol{G}_{(n)})$. The
Frobenius norm of a tensor is
$\|\boldsymbol{\mathcal{G}}\|_{F}^{2}=\sum_{i_{1},\dots,i_{d}}g_{i_{1},\dots,i_{d}}^{2}$.
We refer to Kolda & Bader (2009) for a detailed review of tensor operations.

###  2.2 Vectorized POD-DEIM

We first review the vectorized form of POD-DEIM from which our tensor-based
method is derived. We refer to this method as vector-DEIM, for short. In
vector-DEIM approach for sensor placement (Manohar et al., 2018a; Clark et
al., 2018), the training data is constructed as the snapshot matrix,
$\boldsymbol{G}=\begin{bmatrix}\boldsymbol{g}_{1}&\dots&\boldsymbol{g}_{T}\end{bmatrix}\in\mathbb{R}^{N\times
T}$, where each column $\boldsymbol{g}_{j}\in\mathbb{R}^{N}$ represents a
vectorized snapshot of the quantity of interest $g(\boldsymbol{x},t_{j})$. We
assume that the vectors are centered, which means that the mean is subtracted
from each column. We would like to pick $r\leq\min\\{N,T\\}$ number of sensor
locations at which to collect data (see Fig. 1).

We first compute the truncated singular value decomposition (SVD)
$\boldsymbol{G}\approx\boldsymbol{U}_{r}\boldsymbol{\Sigma}_{r}\boldsymbol{V}_{r}^{\top}$,
where the columns of $\boldsymbol{U}_{r}\in\mathbb{R}^{N\times r}$ coincide
with the POD modes. In DEIM, we compute the column pivoted QR factorization
(Golub & Van Loan, 2013, Section 5.4.2) of $\boldsymbol{U}_{r}$; that is, we
compute
$\boldsymbol{U}_{r}^{\top}\begin{bmatrix}\boldsymbol{P}_{1}&\boldsymbol{P}_{2}\end{bmatrix}=\boldsymbol{Q}_{1}\begin{bmatrix}\boldsymbol{R}_{11}&\boldsymbol{R}_{12}\end{bmatrix}$,
where $\boldsymbol{Q}_{1}\in\mathbb{R}^{r\times r}$ is orthogonal,
$\boldsymbol{R}_{11}\in\mathbb{R}^{r\times r}$ is upper triangular, and
$\boldsymbol{R}_{12}\in\mathbb{R}^{r\times N}$. The matrix
$\boldsymbol{P}=\begin{bmatrix}\boldsymbol{P}_{1}&\boldsymbol{P}_{2}\end{bmatrix}$
is a permutation matrix. Suppose we have
$\boldsymbol{P}_{1}=\begin{bmatrix}\boldsymbol{I}(:,{i_{1}})&\dots&\boldsymbol{I}(:,i_{r})\end{bmatrix}$,
where $\boldsymbol{I}$ is the $N\times N$ identity matrix. Then, the matrix
$\boldsymbol{P}_{1}\in\mathbb{R}^{N\times r}$ contains columns from the
identity matrix indexed by the set $\mathcal{I}=\\{i_{1},\dots,i_{r}\\}$. The
spatial locations corresponding to the index set $\mathcal{I}$ are used as the
optimal sensor locations, which may not be uniquely determined.

Suppose we want to reconstruct the flow field corresponding to the vector
$\boldsymbol{f}\in\mathbb{R}^{N}$. We collect measurements at the indices
corresponding to $\mathcal{I}$, given by the vector
$\boldsymbol{f}(\mathcal{I})$. In other words, $\boldsymbol{f}(\mathcal{I})$
is the available sensor measurements of the vector $\boldsymbol{f}$. To
reconstruct the full flow field $\boldsymbol{f}$, we use the approximation

$\boldsymbol{f}\approx\boldsymbol{U}_{r}\boldsymbol{U}_{r}(\mathcal{I},:)^{-1}\boldsymbol{f}(\mathcal{I})=\boldsymbol{U}_{r}(\boldsymbol{S}^{\top}\boldsymbol{U}_{r})^{-1}\boldsymbol{S}^{\top}\boldsymbol{f},$
(2)

where $\boldsymbol{S}=\boldsymbol{P}_{1}$. This approximation provides insight
into how the indices $\mathcal{I}$ should be selected. Since
$\boldsymbol{U}_{r}$ has rank $r$, it is guaranteed to have $r$ linearly
independent rows and $\boldsymbol{S}^{\top}\boldsymbol{P}_{1}$ is invertible.
The index set $\mathcal{I}$ is chosen in such a way that the corresponding
rows of $\boldsymbol{U}_{r}$ are well-conditioned. The error in the training
set takes the form

$\displaystyle\|\boldsymbol{G}-\boldsymbol{U}_{r}\boldsymbol{U}_{r}(\mathcal{I},:)^{-1}\boldsymbol{G}(\mathcal{I},:)\|_{F}\leq$
$\displaystyle\>\|(\boldsymbol{S}^{\top}\boldsymbol{U}_{r})^{-1}\|_{2}\left(\sum_{j=r+1}^{\min\\{N,T\\}}\sigma_{j}^{2}(\boldsymbol{G})\right)^{1/2},$
(3)

where $\sigma_{j}(\boldsymbol{G})$ represents the singular values of
$\boldsymbol{G}$. In the above expression for the error,
$\left(\sum_{j=r+1}^{\min\\{N,T\\}}\sigma_{j}^{2}(\boldsymbol{G})\right)^{1/2}=\|\boldsymbol{G}-\boldsymbol{U}_{r}\boldsymbol{U}_{r}^{\top}\boldsymbol{G}\|_{F}$
represents the error in the POD approximation due to the truncated singular
values, which is amplified by the factor
$\|(\boldsymbol{S}^{\top}\boldsymbol{U}_{r})^{-1}\|_{2}$ which arises due to
the DEIM approximation. The error in the test dataset can be obtained using
Lemma 3.2 of Chaturantabut & Sorensen (2010).

###  2.3 Tensor-based POD-DEIM

In vector-DEIM, the snapshots are treated as vectors meaning that the inherent
multidimensional structure of the flow is lost. In order to fully exploit this
multidimensional structure, we use tensor-based methods. We refer to the
resulting tensorized version of POD-DEIM as tensor-DEIM, for short. We
consider the collection of snapshots in the form of a tensor
$\boldsymbol{\mathcal{G}}\in\mathbb{R}^{N_{1}\times\dots\times N_{d}\times T}$
of order $d+1$, where $d$ represents the number of spatial dimensions and
$\prod_{j=1}^{d}N_{j}=N$ is the total number of grid points (see Fig. 1). With
this notation, the snapshot matrix $\boldsymbol{G}$ in vector-DEIM can be
expressed as $\boldsymbol{G}_{(d+1)}^{\top}$; that is the transpose of the
mode-$(d+1)$ unfolding.

Figure 1: Schematic representation of the workflow in vector-DEIM and tensor-
DEIM. The horizontal and vertical black bars mark the entries that are used
for sensor placement (black circles). In tensor-DEIM, some sensors may fall on
land and will be discarded (e.g., the one marked by the red circle).

Suppose we wanted to collect data at $r=\prod_{n=1}^{d}r_{n}$ sensor
locations. In tensor-DEIM, we first compute the truncated SVD of the first $d$
mode unfoldings. That is, we compute
$\boldsymbol{G}_{(n)}\approx\boldsymbol{U}_{r}^{(n)}\boldsymbol{\Sigma}_{r}^{(n)}(\boldsymbol{V}_{r}^{(n)})^{\top}$
where $\boldsymbol{U}_{r}^{(n)}\in\mathbb{R}^{N_{n}\times r_{n}}$. For ease of
notation, we define $\boldsymbol{\Phi}_{n}:=\boldsymbol{U}_{r}^{(n)}$. Next,
we compute the column-pivoted QR factorization of $\boldsymbol{\Phi}_{n}$ as

$\boldsymbol{\Phi}_{n}^{\top}\begin{bmatrix}\boldsymbol{P}_{1}^{(n)}&\boldsymbol{P}_{2}^{(n)}\end{bmatrix}=\boldsymbol{Q}_{1}^{(n)}\begin{bmatrix}\boldsymbol{R}_{11}^{(n)}&\boldsymbol{R}_{12}^{(n)}\end{bmatrix}\qquad
1\leq n\leq d.$ (4)

Here
$\begin{bmatrix}\boldsymbol{P}_{1}^{(n)}&\boldsymbol{P}_{2}^{(n)}\end{bmatrix}$
is a permutation matrix. Once again, for ease of notation, we set
$\boldsymbol{S}_{n}:=\boldsymbol{P}_{1}^{(n)}$; this matrix contains the
columns from the $N_{n}\times N_{n}$ identity matrix corresponding to the
indices $\mathcal{I}_{n}=\\{i_{1}^{(n)},\dots,i_{r_{n}}^{(n)}\\}$. Then the
sensors can be placed at the spatial locations corresponding to the index set
$\mathcal{I}_{1}\times\dots\times\mathcal{I}_{d}$. Note that when $d=1$
tensor-DEIM reduces to vector-DEIM.

Given a new data
$\boldsymbol{\mathcal{F}}\in\mathbb{R}^{N_{1}\times\dots\times N_{d}}$, such
that $\boldsymbol{f}=\text{vec}(\boldsymbol{\mathcal{F}})$, we only need to
collect data at the indices corresponding to
$\mathcal{I}_{1}\times\dots\times\mathcal{I}_{d}$, that is we measure
$\boldsymbol{\mathcal{F}}(\mathcal{I}_{1},\dots,\mathcal{I}_{d})$. To
reconstruct the flow field from these measurements, we compute

$\boldsymbol{\mathcal{F}}\approx\boldsymbol{\mathcal{F}}(\mathcal{I}_{1},\dots,\mathcal{I}_{d})\times_{n=1}^{d}\boldsymbol{\Phi}_{n}(\boldsymbol{S}_{n}^{\top}\boldsymbol{\Phi}_{n})^{-1}=\boldsymbol{\mathcal{F}}\times_{n=1}^{d}\boldsymbol{\Phi}_{n}(\boldsymbol{S}_{n}^{\top}\boldsymbol{\Phi}_{n})^{-1}\boldsymbol{S}_{n}^{\top}.$
(5)

The second expression, while equivalent to the first, is more convenient for
the forthcoming error analysis. The approximation just derived satisfies the
interpolation property; that is, the approximation exactly matches the
function $\boldsymbol{\mathcal{F}}$ at the sensor locations assuming the
measurements are noise-free. To see this, denote
$\boldsymbol{\mathcal{F}}_{\rm
TDEIM}:=\boldsymbol{\mathcal{F}}\times_{n=1}^{d}\boldsymbol{\Phi}_{n}(\boldsymbol{S}_{n}^{\top}\boldsymbol{\Phi}_{n})^{-1}\boldsymbol{S}_{n}^{\top}.$
Then

$\displaystyle\boldsymbol{\mathcal{F}}_{\rm
TDEIM}(\mathcal{I}_{1},\dots,\mathcal{I}_{d})=$
$\displaystyle\>\boldsymbol{\mathcal{F}}_{\rm
TDEIM}\times_{n=1}^{d}\boldsymbol{S}_{n}^{\top}=\boldsymbol{\mathcal{F}}\times_{n=1}^{d}\boldsymbol{S}_{n}^{\top}\boldsymbol{\Phi}_{n}(\boldsymbol{S}_{n}^{\top}\boldsymbol{\Phi}_{n})^{-1}\boldsymbol{S}_{n}^{\top}$
(6) $\displaystyle=$
$\displaystyle\>\boldsymbol{\mathcal{F}}\times_{n=1}^{d}\boldsymbol{S}_{n}^{\top}=\boldsymbol{\mathcal{F}}(\mathcal{I}_{1},\dots,\mathcal{I}_{d}).$

The following theorem provides an expression for the error in the
approximation as applied to the training dataset.

###### Theorem 1

Suppose $\boldsymbol{\Phi}_{n}\in\mathbb{R}^{N_{n}\times r_{n}}$ is computed
from the truncated rank-$r_{n}$ SVD of the mode unfolding
$\boldsymbol{G}_{(n)}$ and $\boldsymbol{S}_{n}$ is obtained by computing
column pivoted QR of $\boldsymbol{\Phi}_{n}^{\top}$ such that
$\boldsymbol{S}_{n}^{\top}\boldsymbol{\Phi}_{n}$ is invertible for $1\leq
n\leq d$. Define
$\boldsymbol{\Pi}_{n}:=\boldsymbol{\Phi}_{n}(\boldsymbol{S}_{n}^{\top}\boldsymbol{\Phi}_{n})^{-1}\boldsymbol{S}_{n}^{\top}$
and assume $1\leq r_{n}<N_{n}$ for $1\leq n\leq d$. Then

$\|\boldsymbol{\mathcal{G}}-\boldsymbol{\mathcal{G}}\times_{n=1}^{d}\boldsymbol{\Pi}_{n}\|_{F}\leq\left(\prod_{n=1}^{d}\|(\boldsymbol{S}_{n}^{\top}\boldsymbol{\Phi}_{n})^{-1}\|_{2}\right)\left(\sum_{n=1}^{d}\sum_{k>r_{n}}\sigma_{k}^{2}(\boldsymbol{G}_{(n)})\right)^{1/2}.$
(7)

The proof of this theorem is given in Appendix A. The interpretation of this
theorem is as follows: the term
$\left(\sum_{n=1}^{d}\sum_{k>r_{n}}\sigma_{k}^{2}(\boldsymbol{G}_{(n)})\right)^{1/2}$
represents the error due to the truncated SVD in each mode, and
$\left(\prod_{n=1}^{d}\|(\boldsymbol{S}_{n}^{\top}\boldsymbol{\Phi}_{n})^{-1}\|_{2}\right)$
represents the amplification due to the selection operator across each mode.
The upper bound (7) is similar to the upper bound (3) for vector-DEIM.
However, it is difficult to establish which bound is tighter a priori. As will
be shown in §3, numerical evidence strongly suggests that the error due to
tensor-DEIM is much lower than vector-DEIM.

The error in the test sample can be determined using Proposition 1 from
Kirsten (2022), which gives

$\|\boldsymbol{\mathcal{F}}-\boldsymbol{\mathcal{F}}\times_{n=1}^{d}\boldsymbol{\Pi}_{n}\|_{F}\leq\left(\prod_{n=1}^{d}\|(\boldsymbol{S}_{n}^{\top}\boldsymbol{\Phi}_{n})^{-1}\|_{2}\right)\|\boldsymbol{\mathcal{F}}-\boldsymbol{\mathcal{F}}\times_{n=1}^{d}\boldsymbol{\Phi}_{n}\boldsymbol{\Phi}_{n}^{\top}\|_{F}.$
(8)

If strong rank-revealing QR (sRRQR) algorithm (Gu & Eisenstat, 1996, Algorithm
4) with parameter $f=2$ is used to compute the selection operators
$\boldsymbol{S}_{n}$, then the bound in Theorem 1 simplifies to

$\|\boldsymbol{\mathcal{G}}-\boldsymbol{\mathcal{G}}\times_{n=1}^{d}\boldsymbol{\Pi}_{n}\|_{F}\leq\left(\prod_{n=1}^{d}\sqrt{1+4r_{n}(N_{n}-r_{n})}\right)\left(\sum_{n=1}^{d}\sum_{k>r_{n}}\sigma_{k}^{2}(\boldsymbol{G}_{(n)})\right)^{1/2}.$
(9)

See (Drmač & Saibaba, 2018, Lemma 2.1) for details.

####  2.3.1 Storage cost

We only need to store the bases
$\boldsymbol{\Phi}_{n}\in\mathbb{R}^{N_{n}\times r_{d}}$, which costs
$\sum_{n=1}^{d}N_{n}r_{n}$ entries. Compare this with vector-DEIM which
requires $r\prod_{n=1}^{d}N_{n}=rN$ entries. Assuming $N_{1}=\dots=N_{d}$ and
$r_{1}=\dots=r_{d}$, the ratio of storage cost of tensor-DEIM to that of
vector-DEIM is

$\text{ratio}_{\text{stor}}:=\frac{\sum_{n=1}^{d}r_{n}N_{n}}{rN}=\frac{dr_{1}N_{1}}{r_{1}^{d}N_{1}^{d}}=\frac{d}{r_{1}^{d-1}N_{1}^{d-1}}.$
(10)

Therefore, the compression available using tensors can be substantial when the
dimension $d$, the grid size $N_{1},$ and/or the number of sensors $r_{1}$ are
large. As an illustration in three spatial dimensions, $d=3$, let
$N_{1}=N_{2}=N_{3}=1024$ grid points, and the target rank
$r_{1}=r_{2}=r_{3}=25$; the fraction of the storage cost of tensor-DEIM bases,
compared to vector-DEIM, is $3/(25^{2}\cdot 1024^{2})\times 100\%\approx
4.6\times 10^{-7}\%$. Similar savings in terms of storage costs were also
reported in Gelss et al. (2019) who used tensor-based methods for data-driven
discovery of governing equations.

####  2.3.2 Computational Cost

The cost of vector-DEIM is essentially the cost of computing an SVD on a
$N\times T$ matrix, which is $\mathcal{O}(NT^{2})$ floating point operations
(flops) assuming $T\leq N$. The cost of computing the tensor-DEIM bases is
$\mathcal{O}(T\prod_{j=1}^{d}N_{j}\sum_{n=1}^{d}N_{n})$ flops. Tensor-DEIM is
slightly more expensive since it has to compute $d$ different SVDs compared to
vector-DEIM. This computational cost can be amortized by using the
sequentially truncated higher-order SVD of Vannieuwenhoven et al. (2012). The
cost of computing the indices that determine the sensor locations is
$\mathcal{O}(\sum_{n=1}^{d}N_{n}r_{n}^{2})$ flops for tensor-DEIM which is
much cheaper than vector-DEIM which costs $\mathcal{O}(Nr^{2})$ flops.

When the training dataset is large, computing the truncated SVD approximation
can be expensive. One way to accelerate this computation is to use randomized
methods as in (Minster et al., 2020). Suppose the randomized higher order SVD
is used, then the computational cost is
$\mathcal{O}(T\prod_{j=1}^{d}N_{j}\sum_{n=1}^{d}r_{n})$ flops. This cost is
substantially less than the cost of both vector-DEIM and tensor-DEIM.

##  3 Results and discussion

In the numerical experiments, we use QR with column pivoting as implemented in
MATLAB for both vector-DEIM and tensor-DEIM.

###  3.1 Kolmogorov flow

Kolmogorov flow refers to a turbulent flow with periodic boundary conditions
and a sinusoidal forcing. Here we consider the two-dimensional Kolmogorov
flow,

$\partial_{t}\omega+\mathbf{u}\cdot\nabla\omega=\nu\Delta\omega-n\cos(ny),$
(11)

where $\mathbf{u}=(\partial_{y}\psi,-\partial_{x}\psi)$ is the fluid velocity
field with the stream function $\psi(x,y,t)$ and $\omega=-\Delta\psi$ is the
vorticity field. We consider a two-dimensional domain
$\boldsymbol{x}=(x,y)\in[0,2\pi]\times[0,2\pi]$ with periodic boundary
conditions. The forcing wavenumber is $n=4$ and $\nu=\mbox{Re}^{-1}$ is the
inverse of the Reynolds number Re.

Figure 2: Comparing vector-DEIM and tensor-DEIM on the Kolmogorov flow. RE
denotes the relative error of the reconstruction.

We numerically solve equation (11) using a standard pseudo-spectral method
with $128\times 128$ modes and $2/3$ dealiasing. The temporal integration is
carried out with the embedded Runge–Kutta scheme of Dormand & Prince (1980).
The initial condition is random and is evolved long enough to ensure that the
initial transients have died out before any data collection is performed. Then
$10^{3}$ vorticity snapshots are saved, each $\Delta t=5$ time units apart.
The time increment $\Delta t$ is approximately 10 times the eddy turnover time
$\tau_{e}\simeq 0.5$ of the flow, ensuring that the snapshots are not strongly
correlated. First $75\%$ of the data are used for training. The remaining 25%
are used for testing. The training data forms the data tensor
$\boldsymbol{\mathcal{G}}\in\mathbb{R}^{128\times 128\times 750}$.

Figure 3: Optimal sensor locations obtained by vector-DEIM (blue circles) and
tensor-DEIM (red circles) for the Kolmogorov flow. The number of sensors are
(a) 25, (b) 100, and (c) 400 as in Fig. 2.

We have verified that $10^{3}$ snapshots are adequate for the results to have
converged. For instance, changing the number of snapshots to $800$ did not
significantly alter the results reported below. Furthermore, choosing the
training snapshots at random, instead of the first $75\%$, did not affect the
reported results.

Figure 2 compares reconstruction results using the conventional vector-DEIM
and our tensor-based DEIM. These reconstructions are performed for a vorticity
field in the testing data set. As the number of sensors increases, both
reconstructions improve. For the same number of sensors, our tensor-based
method always returns a more accurate reconstruction compared to vector-DEIM.
Furthermore, the storage cost for tensor-DEIM is much lower. For instance, for
reconstruction from 400 sensors, storing the tensor bases only requires
$0.07\%$ of the memory required by vector-DEIM. Figure 3 compares the location
of sensors used in vector-DEIM and tensor-DEIM. For relatively small number of
25 sensors, the optimal sensor locations corresponding to vector-DEIM and
tensor-DEIM are significantly different. However, as the number of sensors
increases the difference diminishes.

Figure 4: Relative flow reconstruction error for the Kolmogorov flow. The
number of sensors is denoted by $n_{s}$.

Figure 4 shows the relative reconstruction error as a function of the number
of sensors. For each snapshot in the testing data set, we compute this error
separately. The symbols in Fig. 4 mark the mean relative error taken over the
250 snapshots in this data set. The error bars show one standard deviation of
the error. In every case, tensor-DEIM outperforms its vectorized counterpart
as assessed by the mean relative error. In addition, the standard deviation of
the relative error is smaller when using tensor-DEIM as compared to vector-
DEIM.

###  3.2 Sea surface temperature

As the second test case, we consider the reconstruction of global ocean
surface temperature. The dataset is publicly available at NOAA Optimum
Interpolation SST V2, noaa.oisst.v2 (Reynolds et al., 2002). The temperature
distribution is affected by the complex ocean flow dynamics resulting in
seasonal variations. This dataset is in the form of a time series in which a
snapshot is recorded every week in the span of 1990-2016 and data is available
at a resolution of $1^{\degree}\times 1^{\degree}$. In total there are $1688$
snapshots, which we split into a training set of $1200$ (roughly $71\%$) and a
testing set of $488$ snapshots.

Figure 5: Comparing vector-DEIM and tensor-DEIM on the sea surface temperature
dataset on March 3, 2013. The white dots indicate sensor locations and RE
represents relative error. The color bar represents temperature in degrees
centigrade.

Figure 5 shows the reconstruction using tensor-DEIM and vector-DEIM compared
to the ground truth. The tensor-DEIM places sensors in a rectangular array and
some sensors may fall within the land surface. These sensors are discarded,
and only the ones on the ocean surface are retained. In this problem instance,
there are $764$ sensors. To measure the accuracy of the reconstruction, we use
the relative error of the fields centered around the mean sea surface
temperature. As is seen from the figure, the reconstruction error from tensor-
DEIM is superior to that of vector-DEIM. Storing the tensor bases only
requires $0.04\%$ of the memory footprint required by vector-DEIM.

Figure 6 shows the relative error as a function of the number of sensors. As
in the previous experiment, we compute the error over each snapshot in the
test dataset and display the mean over the $488$ snapshots with the error bars
indicating one standard deviation of the error. As can be seen, once again
tensor-DEIM outperforms its vectorized counterpart both in terms of having a
lower mean and standard deviation. The superiority of tensor-DEIM becomes more
and more pronounced as the number of sensors increases. Note that the error in
both methods does not decrease monotonically with more sensors. This is
because the error has two contributions: the approximation of the snapshot by
the basis and the amplification factor due to the interpolation; see Eq. (7).
While the first contribution is non-increasing, the second contribution may
increase with an increasing number of sensors; hence the overall error may
increase.

Figure 6: Comparing the contributions to the error in vector-DEIM and tensor-
DEIM on a training set with first $71\%$ of the data. (left) Projection error
normalized by the norm of the snapshot (center) Amplification factor, (right)
Overall relative error.

To explain why the interpolation error for vector-DEIM gets worse with an
increasing number of sensors, we plot the two different contributions to the
error. By Lemma 3.2 of Chaturantabut & Sorensen (2010), the testing error of
vector-DEIM is bounded as

$\|\boldsymbol{f}-\boldsymbol{U}_{r}(\boldsymbol{S}^{\top}\boldsymbol{U}_{r})^{-1}\boldsymbol{S}^{\top}\boldsymbol{f}\|_{2}\leq\|(\boldsymbol{S}^{\top}\boldsymbol{U}_{r})^{-1}\|_{2}\|\boldsymbol{f}-\boldsymbol{U}_{r}\boldsymbol{U}_{r}^{\top}\boldsymbol{f}\|_{2}.$
(12)

The error bound has two contributions:
$\|\boldsymbol{f}-\boldsymbol{U}_{r}\boldsymbol{U}_{r}^{\top}\boldsymbol{f}\|_{2}$
which we call the projection error, and
$\|(\boldsymbol{S}^{\top}\boldsymbol{U}_{r})^{-1}\|_{2}$ which we call as the
amplification factor. A similar classification can be done for the error due
to tensor-DEIM, using (8). The left panel of Fig. 6 shows the projection error
normalized by $\|\boldsymbol{f}\|_{2}$ and the center panel shows the
amplification factor. We see that while the projection error decreases with an
increasing number of sensors, the amplification factor increases, resulting in
an overall increased error on average. On the other hand, for tensor-DEIM the
projection error is higher compared to vector-DEIM. However, the amplification
factor is nearly constant, and the overall error for tensor-DEIM decreases
with an increasing number of sensors.

To explore this further, we repeat the experiment but with a randomly
generated training and testing split; more precisely, we randomly choose $\sim
71\%$ of the data to be the training set and the remaining to be the test set.
The results are shown in Fig. 7. Similar to Fig. 6, we plot the contributions
to the error in the left and center panels and the overall error in the right
panel. Qualitatively, we see a similar trend as in the previous experiment.
However, one major difference is that the overall mean error of vector-DEIM is
now much closer to that of tensor-DEIM. But tensor-DEIM still has a lower mean
error and a much smaller standard deviation.

Note that this is in contrast with the Kolmogorov flow data, where
randomization had no significant effect on the results. We attribute this to
the fact that the Kolmogorov data contains a large and well-separated set of
snapshots. Hence, the splitting of the data into training and test subsets
does not play a major role on the sampling of the attractor.

Figure 7: Comparing the contributions to the error in vector-DEIM and tensor-
DEIM on a randomly selected training set. (left) Projection error normalized
by the norm of the snapshot (center) Amplification factor, (right) Overall
relative error.

Finally, we turn to the problem of El Niño-Souther Oscillation (ENSO), i.e.,
cycles of warm (El Niño) and cold (La Niña) water temperature in the Pacific
ocean near the equator. It is known that reconstructing ENSO features from
sparse measurements is challenging (Manohar et al., 2018b; Maulik et al.,
2020). Here, we focus on the El Niño event of the winter of 1997-1998, which
is well-known for its intensity off the coast of Peru. In particular, Fig. 8
(left panel) shows the SST in December 1997. This snapshot lies in the test
dataset when the sets are chosen at random as opposed to sequentially. We see
that the overall reconstruction error using vector-DEIM is larger compared to
tensor-DEIM. More importantly, the spatial error near the El Niño oscillation
is significantly larger when using vector-DEIM. Therefore, tensor-DEIM more
successfully reconstructs the El Niño patterns. This is quite counter-
intuitive since vector-DEIM tends to place more sensors in the El Niño region
(see the right panel of figure 5 off the coast of Peru).

Figure 8: Comparing the spatial distributions of the error. Left: true SST
from December 1997. Center: spatial distribution of error (reconstructed SST
minus true SST) for vector-DEIM. Right: spatial distribution of error for
tensor-DEIM.

###  3.3 Three-dimensional unsteady flow

We consider a dataset obtained by the simulation of an incompressible 3D flow
around a CAD model of the research vessel Tangaroa (Popinet et al. (2004)).
This dataset is obtained by large-eddy simulation using the Gerris flow solver
(Popinet, 2004). The flow variables are resampled onto a regular grid in the
spatial region of interest, $[-0.35,0.65]\times[-0.3,0.3]\times[-0.5,-0.3]$.
The reported spatial variables are dimensionless, normalized with the
characteristic length scale $L=276\,$m, four times the ship length. Flow
velocity is normalized by the inflow velocity $U=1\,$ m/s. We consider only
the $u$ component of the velocity $(u,v,w)$ and subsampled the data to
consider a grid size $150\times 90\times 60$. We split the available $201$
snapshots into a training dataset with $150$ snapshots ($\sim 75\%$) and a
test dataset with $61$ snapshots.

Figure 9: (left) The true u-velocity from the test dataset. (right)
Reconstruction using tensor-DEIM with $125$ sensors. Black dots indicate the
position of the sensors. RE denotes the relative error.

Due to the enormous size of the resulting tensor, we used randomized SVD to
compute the factor matrices $\boldsymbol{\Phi}_{n}$ for $1\leq n\leq 3$ via
the MATLAB command `svdsketch`. We choose $5\times 5\times 5=125$ sensors to
reconstruct the flow field. The true $u$-velocity for a snapshot in the test
dataset is plotted in the left panel of Fig. 9. The reconstruction using
tensor-DEIM is displayed in the right panel of the same figure. The relative
error in the reconstruction is around $9\%$ suggesting that tensor-DEIM is
adequately reconstructing the snapshot. The corresponding error for the same
number of sensors using vector-DEIM is $20.08\%$. As in the previous examples,
it is seen that tensor-DEIM is far more accurate than vector-DEIM for the same
number of sensors. The average relative error over the entire test data is
$9.07\%$ for tensor-DEIM and $41.79\%$ for vector-DEIM. For some snapshots,
the error in vector-DEIM is as high as $75\%$ and the error appears to
increase for the later snapshots. Finally, the cost of storing the tensor-DEIM
factor matrices is $0.0015\%$ of the cost of storing the vector-DEIM basis.
Once again, tensor-DEIM proves to be more accurate and far more storage
efficient compared to vector-DEIM.

##  4 Conclusions

Our results make a strong case for using tensor-based methods for sensor
placement and flow reconstruction. In particular, our tensor-based method
significantly increases the reconstruction accuracy and reduces its storage
cost. Our numerical examples show that the relative reconstruction errors are
comparable to vectorized methods when the number of sensors is small. However,
as the number of sensors increases, our tensor-based method is 2-3 times more
accurate than its vectorized counterpart. The improvements in terms of the
storage cost are even more striking: tensor-DEIM requires only $0.0015-0.07\%$
of the memory required by vector-DEIM.

Future work could include the application of the method to reacting flows. In
such flows, the multidimensional nature of our tensor-based method allows for
separate optimal sensor placement for the flow field and each chemical
species. Although our method is a tensorized version of DEIM, a similar
tensor-based approach can be applied to other flow reconstruction methods such
as gappy POD, sparsity-promoting methods, and autoencoders.

## Acknowledgements

The authors would like to acknowledge the NOAA Optimum Interpolation (OI) SST
V2 data provided by the NOAA PSL, Boulder, Colorado, USA, from their website
at https://psl.noaa.gov. We would also like to thank the Computer Graphics Lab
(ETH Zurich) for making the 3D unsteady flow data used in section 3.3 publicly
available.

## Funding

This work was supported, in part, by the National Science Foundation through
the awards DMS-1821149 and DMS-1745654, and Department of Energy through the
award DE-SC0023188.

##  Appendix A Proof of Theorem 1

We will need the following notation for the proof. We express
$\boldsymbol{\mathcal{Y}}=\boldsymbol{\mathcal{G}}\times_{1}\boldsymbol{A}_{1}\dots\times\boldsymbol{A}_{d}$
in terms of unfoldings as
$\boldsymbol{\mathcal{Y}}_{d}=\boldsymbol{A}_{d}\boldsymbol{Y}_{(d)}(\boldsymbol{A}_{d-1}\otimes\dots\otimes\boldsymbol{A}_{1})^{\top}$.
Here we use $\otimes$ to denote the Kronecker product of two matrices.

The proof is similar to (Kirsten, 2022, Proposition 1). Using the properties
of matrix-unfoldings and Kronecker products, we can write an equivalent
expression for the error

$\|\boldsymbol{\mathcal{G}}-\boldsymbol{\mathcal{G}}\times_{n=1}^{d}\boldsymbol{\Pi}_{n}\|_{F}=\|(\boldsymbol{I}-\otimes_{n=d}^{1}\boldsymbol{\Pi}_{n})\boldsymbol{G}_{(d+1)}^{\top}\|_{F}.$

From
$\boldsymbol{\Pi}_{n}\boldsymbol{\Phi}_{n}\boldsymbol{\Phi}_{n}^{\top}=\boldsymbol{\Phi}_{n}\boldsymbol{\Phi}_{n}^{\top}$,
we get
$(\boldsymbol{I}-\otimes_{n=d}^{1}\boldsymbol{\Pi}_{n})(\boldsymbol{I}-\otimes_{n=d}^{1}\boldsymbol{\Phi}_{n}\boldsymbol{\Phi}_{n}^{\top})=\boldsymbol{I}-\otimes_{n=d}^{1}\boldsymbol{\Pi}_{n}.$.
Using this result and submultiplicativity

$\displaystyle\|(\boldsymbol{I}-\otimes_{n=d}^{1}\boldsymbol{\Pi}_{n})\boldsymbol{G}_{(d+1)}^{\top}\|_{F}=$
$\displaystyle\>\|(\boldsymbol{I}-\otimes_{n=d}^{1}\boldsymbol{\Pi}_{n})(\boldsymbol{I}-\otimes_{n=d}^{1}\boldsymbol{\Phi}_{n}\boldsymbol{\Phi}_{n}^{\top})\boldsymbol{G}_{(d+1)}^{\top}\|_{F}$
$\displaystyle\leq$
$\displaystyle\>\|(\boldsymbol{I}-\otimes_{n=d}^{1}\boldsymbol{\Pi}_{n})\|_{2}\|(\boldsymbol{I}-\otimes_{n=d}^{1}\boldsymbol{\Phi}_{n}\boldsymbol{\Phi}_{n}^{\top})\boldsymbol{G}_{(d+1)}^{\top}\|_{F}.$

Since $\boldsymbol{\Pi}_{n}$ is an oblique projector, so are
$\otimes_{j=1}^{d}\boldsymbol{\Pi}_{n}$ and
$\boldsymbol{I}-\otimes_{j=1}^{d}\boldsymbol{\Pi}_{n}$. By (Szyld, 2006,
Theorem 2.1), which applies since $\otimes_{j=1}^{d}\boldsymbol{\Pi}_{n}$ is
neither zero nor the identity,
$\|\boldsymbol{I}-\otimes_{j=1}^{d}\boldsymbol{\Pi}_{n}\|_{2}=\|\otimes_{j=1}^{d}\boldsymbol{\Pi}_{n}\|_{2}=\prod_{n=1}^{d}\|\boldsymbol{\Pi}_{n}\|_{2}$.
In the last step, we have used the fact that the largest singular value of a
Kronecker product is the product of the largest singular values. Therefore,

$\begin{aligned}
\|(\boldsymbol{I}-\otimes_{n=d}^{1}\boldsymbol{\Pi}_{n})\boldsymbol{G}_{(d+1)}^{\top}\|_{F}\leq&\>\left(\prod_{n=1}^{d}\|\boldsymbol{\Pi}_{n}\|_{2}\right)\|(\boldsymbol{I}-\otimes_{n=d}^{1}\boldsymbol{\Phi}_{n}\boldsymbol{\Phi}_{n}^{\top})\boldsymbol{G}_{(d+1)}^{\top}\|_{F}\\\
=&\>\left(\prod_{n=1}^{d}\|(\boldsymbol{S}_{n}^{\top}\boldsymbol{\Phi}_{n})^{-1}\|_{2}\right)\|\boldsymbol{\mathcal{G}}-\boldsymbol{\mathcal{G}}\times_{n=1}^{d}\boldsymbol{\Phi}_{n}\boldsymbol{\Phi}_{n}^{\top}\|_{F}.\end{aligned},$

since $\boldsymbol{\Phi}_{n}$ and $\boldsymbol{S}_{n}$ have orthonormal
columns, so
$\|\boldsymbol{\Pi}_{n}\|_{2}=\|(\boldsymbol{S}_{n}^{\top}\boldsymbol{\Phi}_{n})^{-1}\|_{2}$.
The result follows from (Vannieuwenhoven et al., 2012, Corollary 5.2).

## References

  * Bui-Thanh et al. (2004) Bui-Thanh, T., Damodaran, M. & Willcox, K. 2004 Aerodynamic data reconstruction and inverse design using proper orthogonal decomposition. AIAA Journal 42 (8), 1505–1516.
  * Callaham et al. (2019) Callaham, J. L., Maeda, K. & Brunton, S. L. 2019 Robust flow reconstruction from limited measurements via sparse representation. Phys. Rev. Fluids 4, 103907.
  * Carlberg et al. (2019) Carlberg, K. T., Jameson, A., Kochenderfer, M. J., Morton, J., Peng, L. & Witherden, F. D. 2019 Recovering missing CFD data for high-order discretizations using deep neural networks and dynamics learning. Journal of Computational Physics 395, 105–124.
  * Chaturantabut & Sorensen (2010) Chaturantabut, S. & Sorensen, D. C. 2010 Nonlinear model reduction via discrete empirical interpolation. SIAM Journal on Scientific Computing 32 (5), 2737–2764.
  * Chu & Farazmand (2021) Chu, B. & Farazmand, M. 2021 Data-driven prediction of multistable systems from sparse measurements. Chaos 31 (6), 063118\.
  * Clark et al. (2018) Clark, E., Askham, T., Brunton, S. L. & Kutz, J. N. 2018 Greedy sensor placement with cost constraints. IEEE Sensors Journal 19 (7), 2642–2656.
  * Clark et al. (2021) Clark, E., Brunton, S. L. & Kutz, J. N. 2021 Multi-fidelity sensor selection: Greedy algorithms to place cheap and expensive sensors with cost constraints. IEEE Sensors Journal 21 (1), 600–611.
  * Dang et al. (2021) Dang, F., Nasreen, S. & Zhang, F. 2021 DMD-based background flow sensing for AUVs in flow pattern changing environments. IEEE Robotics and Automation Letters 6 (3), 5207–5214.
  * Dormand & Prince (1980) Dormand, J. R. & Prince, P. J. 1980 A family of embedded Runge–Kutta formulae. J. Comp. App. Math. 6 (1), 19 – 26.
  * Drmač & Gugercin (2016) Drmač, Z. & Gugercin, S. 2016 A new selection operator for the discrete empirical interpolation method—Improved a priori error bound and extensions. SIAM Journal on Scientific Computing 38 (2), A631–A648.
  * Drmač & Saibaba (2018) Drmač, Z. & Saibaba, A. K. 2018 The discrete empirical interpolation method: Canonical structure and formulation in weighted inner product spaces. SIAM Journal on Matrix Analysis and Applications 39 (3), 1152–1180.
  * Erichson et al. (2020) Erichson, N. B., Mathelin, L., Yao, Z., Brunton, S. L., Mahoney, M. W. & Kutz, J. N. 2020 Shallow neural networks for fluid flow reconstruction with limited sensors. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 476 (2238), 20200097.
  * Everson & Sirovich (1995) Everson, R. & Sirovich, L. 1995 Karhunen–Loève procedure for gappy data. J. Opt. Soc. Am. A 12 (8), 1657–1664.
  * Fukami et al. (2019) Fukami, K., Fukagata, K. & Taira, K. 2019 Super-resolution reconstruction of turbulent flows with machine learning. Journal of Fluid Mechanics 870, 106–120.
  * Gelss et al. (2019) Gelss, P., Klus, S., Eisert, J. & Schütte, C. 2019 Multidimensional Approximation of Nonlinear Dynamical Systems. Journal of Computational and Nonlinear Dynamics 14 (6), 061006.
  * Golub & Van Loan (2013) Golub, G. H. & Van Loan, Charles F. 2013 Matrix Computations, 4th edn. Baltimore: The Johns Hopkins University Press.
  * Gu & Eisenstat (1996) Gu, M. & Eisenstat, S. C. 1996 Efficient algorithms for computing a strong rank-revealing QR factorization. SIAM Journal on Scientific Computing 17 (4), 848–869.
  * Kirsten (2022) Kirsten, G. 2022 Multilinear POD-DEIM model reduction for 2D and 3D semilinear systems of differential equations. Journal of Computational Dynamics 9 (2), 159.
  * Kolda & Bader (2009) Kolda, T. G. & Bader, B. W. 2009 Tensor decompositions and applications. SIAM Rev. 51 (3), 455–500.
  * Kramer et al. (2017) Kramer, B., Grover, P., Boufounos, P., Nabi, S. & Benosman, M. 2017 Sparse sensing and DMD-based identification of flow regimes and bifurcations in complex flows. SIAM Journal on Applied Dynamical Systems 16 (2), 1164–1196.
  * Manohar et al. (2018a) Manohar, K., Brunton, B. W, Kutz, J. N. & Brunton, S. L. 2018a Data-driven sparse sensor placement for reconstruction: Demonstrating the benefits of exploiting known patterns. IEEE Control Systems Magazine 38 (3), 63–86.
  * Manohar et al. (2018b) Manohar, K., Brunton, B. W., Kutz, J. N. & Brunton, S. L. 2018b Data-driven sparse sensor placement for reconstruction: Demonstrating the benefits of exploiting known patterns. IEEE Control Systems Magazine 38 (3), 63–86.
  * Maulik et al. (2020) Maulik, R., Fukami, K., Ramachandra, N., Fukagata, K. & Taira, K. 2020 Probabilistic neural networks for fluid flow surrogate modeling and data recovery. Phys. Rev. Fluids 5, 104401.
  * Minster et al. (2020) Minster, R., Saibaba, A. K. & Kilmer, M. E. 2020 Randomized algorithms for low-rank tensor decompositions in the Tucker format. SIAM Journal on Mathematics of Data Science 2 (1), 189–215.
  * Peherstorfer et al. (2014) Peherstorfer, B., Butnaru, D., Willcox, K. & Bungartz, H.-J. 2014 Localized discrete empirical interpolation method. SIAM Journal on Scientific Computing 36 (1), A168–A192.
  * Peherstorfer et al. (2020) Peherstorfer, B., Drmač, Z. & Gugercin, S. 2020 Stability of discrete empirical interpolation and gappy proper orthogonal decomposition with randomized and deterministic sampling points. SIAM Journal on Scientific Computing 42 (5), A2837–A2864.
  * Popinet (2004) Popinet, S. 2004 Free computational fluid dynamics. ClusterWorld 2 (6).
  * Popinet et al. (2004) Popinet, Stéphane, Smith, Murray & Stevens, Craig 2004 Experimental and numerical study of the turbulence characteristics of airflow around a research vessel. Journal of Atmospheric and Oceanic Technology 21 (10), 1575–1589.
  * Reynolds et al. (2002) Reynolds, R. W., Rayner, N. A., Smith, T. M., Stokes, D. C. & Wang, W. 2002 An improved in situ and satellite SST analysis for climate. Journal of Climate 15 (13), 1609 – 1625.
  * Szyld (2006) Szyld, Daniel B 2006 The many proofs of an identity on the norm of oblique projections. Numerical Algorithms 42 (3), 309–323.
  * Vannieuwenhoven et al. (2012) Vannieuwenhoven, N., Vandebril, R. & Meerbergen, K. 2012 A new truncation strategy for the higher-order singular value decomposition. SIAM Journal on Scientific Computing 34 (2), A1027–A1052.
  * Wang et al. (2021) Wang, Y., Ding, X., Hu, K., Fang, F., Navon, I.M. & Lin, G. 2021 Feasibility of DEIM for retrieving the initial field via dimensionality reduction. Journal of Computational Physics 429, 110005.
  * Willcox (2006) Willcox, K. 2006 Unsteady flow sensing and estimation via the gappy proper orthogonal decomposition. Computers & Fluids 35 (2), 208–226.
  * Yu & Hesthaven (2019) Yu, J. & Hesthaven, J. S. 2019 Flowfield reconstruction method using artificial neural network. AIAA Journal 57 (2), 482–498.


# Repetitive Infection Spreading and Directed Evolution in the Susceptible-
Infected-Recovered-Susceptible Model

Hidetsugu Sakaguchi and Keito Yamasaki  Interdisciplinary Graduate School of
Engineering Sciences, Kyushu University, Kasuga, Fukuoka 816-8580, Japan

###### Abstract

We study two simple mathematical models of the epidemic. At first, we study
the repetitive infection spreading in a simplified SIRS model including the
effect of the decay of the acquired immune. The model is an intermediate model
of the SIRS model including the recruitment and death terms and the SIR model
in which the recovered population is assumed to be never infected again. When
the decay rate $\delta$ of the immune is sufficiently small, the multiple
infection spreading occurs in spikes. The model equation can be reduced to be
a map when the decay rate $\delta$ is sufficiently small, and the spike-like
multiple infection spreading is reproduced in the mapping. The period-doubling
bifurcation and chaos are found in the simplified SIRS model with seasonal
variation. The nonlinear phenomena are reproduced by the map. Next, we study
coupled SIRS equations for the directed evolution where the mutation is
expressed with a diffusion-type term. A kind of reaction-diffusion equation is
derived by the continuum approximation for the infected population $I$. The
reaction-diffusion equation with the linear dependence of infection rate on
the type space has an exact Gaussian solution with a time-dependent average
and variance. The propagation of the Gaussian pulse corresponds to the
successive transitions of the dominant variant.

##  I Introduction

The nonlinear dynamics of epidemic models has been studied by many authors,
since Kermack and McKendrick proposed a mathematical model for epidemics in
1927 KM . The model describes the population dynamics of $S$ (susceptible),
$I$ (infected), and $R$ (recovered), and therefore, it is called the SIR
model. If the population $E$ of the exposed population before the appearance
of symptoms is included, the SIR model is generalized to the SEIR model
Anderson . In the SIR model, the total population is conserved. In a lengthy
period, the total population can grow or decay. The endemic SIR model includes
the effect of the growth and decay Hethcote . If the transition term from the
recovered population to the susceptible population is added to the endemic SIR
model, the SIRS model (susceptible-infected-recovered-susceptible model) is
derived Hethcote . Furthermore, various modified SIR models such as the SIS
model (susceptible-infected-susceptible model), SIRD model (susceptible-
infectious-recovered-deceased model), SEIS (susceptible-exposed-infected-
susceptible model), and so on have been studied for the analysis of the spread
of disease Murray ; Diekmann ; Capasso . Nobel performed a numerical
simulation of the geographic and temporal development of plagues using the
mathematical model with diffusion terms Nobel . The infection spreading on
various lattices and networks has been a topic of statistical mechanics and
complex network science Katori ; Pastor ; Odor . We studied the slow decay of
infection in the SIR model with random infection rate on lattices Sakaguchi .
The epidemic models were also used to study vaccination strategy Shulgin ;
Laguzet . The recent COVID-19 pandemic has excited the study of mathematical
models of epidemic Science . The behaviors in the post-pandemic period are
predicted using various epidemic models Kissler ; Are .

The SIR model is expressed as KM ; Murray .

$\displaystyle\frac{dS}{dt}$ $\displaystyle=$ $\displaystyle-\beta SI,$
$\displaystyle\frac{dI}{dt}$ $\displaystyle=$ $\displaystyle\beta SI-\gamma
I,$ $\displaystyle\frac{dR}{dt}$ $\displaystyle=$ $\displaystyle\gamma I,$ (1)

where $S$, $I$, and $R$ are respectively susceptible, infected, and recovered
populations. The parameters $\beta$ and $\gamma$ denote infection and
recovering rates. The total population $S+I+R$ is conserved. In this model,
the infection spreading occurs only once, and the infected population decays
monotonously after the infection spreads. Multiple spreading often occurs in
the pandemic such as COVID-19. There are several reasons for the multiple
infection spreading. Regulations such as the lockdown and the relaxation might
be related to the multiple spreading. The generation of new variants is
another reason for the waves Dutta . Here, we consider the effect of the decay
of the acquired immune.

The SIRS model including the effect of the day of the acquire immune is
expressed as

$\displaystyle\frac{dS}{dt}$ $\displaystyle=$ $\displaystyle b-\mu S-\beta
SI+\delta R,$ $\displaystyle\frac{dI}{dt}$ $\displaystyle=$
$\displaystyle\beta SI-\gamma I,$ $\displaystyle\frac{dR}{dt}$
$\displaystyle=$ $\displaystyle\gamma I-\delta R,$ (2)

where $b$ is the recruitment rate of the susceptible population, $\mu$ is the
death rate, and $\delta$ is the rate at which the recovered population becomes
susceptible because the effect of the acquired immune wears off. In this
model, $S+I+R$ is not conserved.

In this paper, we study a simplified SIRS model with $b=\mu=0$. The simplified
SIRS model is an intermediate model of the SIR model and the SIRS model. The
simplified SIRS model exhibits nontrivial nonlinear dynamics such as the
spike-like infection waves. In Sec. 2, we study the repetitive waves of
infection spreading. The simplified SIRS model is reduced to the SIR model
when the transition rate $\delta$ from the recovery population to the
susceptible population is zero. Although the infection spreading occurs only
once in the SIR model, the repetitive infection spreading occurs, and the peak
amplitude of $I(t)$ decays in time in the simplified SIRS model even if the
parameter $\delta$ is sufficiently small. The interval between the successive
spikes increases and the decay rate of peak amplitude of $I(t)$ decreases when
$\delta$ approaches 0. We will show that the nonlinear dynamics of the
differential equation can be approximated by a discrete map. In Sec. 3, we
consider the simplified SIRS model in case the infection rate changes
periodically. The repetitive infection waves is maintained owing to the
periodic change of the infection rate, and the periodic and chaotic dynamics
appear. The nonlinear dynamics can also be approximated by the discrete map.
In Sec. 4, we study a one-dimensional lattice model that takes the form of
coupled equations of the SIRS model. We find the directed evolution to higher
infectivity and lower virulence occurs in this simple model. We show that the
partial differential equation obtained by the continuum approximation has a
Gaussian solution with a time-dependent average and variance. In Sec. 5, we
summarize the results.

##  II Repetitive infection spreading in a simplified susceptible-infected-
recovered-susceptible model

Our simplified SIRS model without the recruitment and death terms is expressed
as

$\displaystyle\frac{dS}{dt}$ $\displaystyle=$ $\displaystyle-\beta SI+\delta
R,$ $\displaystyle\frac{dI}{dt}$ $\displaystyle=$ $\displaystyle\beta
SI-\gamma I,$ $\displaystyle\frac{dR}{dt}$ $\displaystyle=$
$\displaystyle\gamma I-\delta R,$ (3)

This model is an intermediate model of Eq. (1) and Eq. (2). The total
population $S+I+R$ is conserved similarly to Eq. (1). Since we study the case
of sufficiently small $\delta$, our model is close to the SIR model. However,
the small perturbation terms $\delta R$ change the system behavior
drastically. That is, only one infection wave appears in Eq. (1), but
repetitive infection waves appear in the model Eq. (3). We will show the
mechanism of the repetitive infection waves using an approximate map derived
from Eq. (3).

Since the total population is constant $S+I+R=N$, Eq. (3) is reduced to

$\displaystyle\frac{dS}{dt}$ $\displaystyle=$ $\displaystyle-\beta
SI+\delta(N-S-I),$ $\displaystyle\frac{dI}{dt}$ $\displaystyle=$
$\displaystyle\beta SI-\gamma I.$ (4)

There are two stationary solutions: $S=N,\;I=0$ and
$S=\gamma/\beta,\;I=\delta(\beta N-\gamma)/\\{\beta(\delta+\gamma)\\}$. If
$\beta N>\gamma$, the first solution of no infection becomes unstable and the
second solution becomes a stable solution. The stability of the second
solution is determined by eigenvalues of the Jacobian matrix:

$A=\left(\begin{array}[]{cc}-\delta(\beta
N-\gamma)/(\delta+\gamma)-\delta,&-\gamma-\delta\\\ \delta(\beta
N-\gamma)/(\delta+\gamma),&0\end{array}\right).$

The trace $t$ of the matrix is $t=-\delta(\beta
N-\gamma)/(\delta+\gamma)-\delta$ and the determinant $D$ is $D=\delta(\beta
N-\gamma)$. Because $\beta N>\gamma$, $D>0$, and $t<0$, the real part of the
eigenvalues is negative, and the solution is stable. If $D>t^{2}/4$, the
eigenvalues are complex. If $\delta$ is sufficiently small, $D>t^{2}/4$ is
satisfied and the damping oscillation occurs. The critical value of
$\delta_{c}$ between the damping oscillation and exponential decay is

$\delta_{c}=\frac{4(\beta N-\gamma)(\delta+\gamma)^{2}}{(\beta
B-\gamma+\delta+\gamma)^{2}}.$

Figure 1: (a) Time evolutions of $I(t)$ at $\delta=0.001$, $N=5$, $\beta=10$,
and $\gamma=20$. The initial condition is $S(0)=4.999$ and $I(0)=0.001$. (b)
Magnification of Fig. 2(a) for $983<t<986$. (c) Time evolution of $S(t)$.

For sufficiently small $\delta$, the damping oscillation of $I(t)$ seems to be
a spiking oscillation. Figure 1(a) shows the time evolutions of $I(t)$ at
$\delta=0.001$, $N=5$, $\beta=10$, and $\gamma=20$. The initial condition is
$S(0)=4.999$ and $I(0)=0.001$. The spike-like infection spreading repeats in
this long timescale. Figure 1(b) is the magnification of Fig. 1(a) for
$983<t<986$. The spreading and ending of infection occur continuously in this
short timescale. Figure 1(c) shows the time evolution of $S(t)$ in the long
timescale. After the infection spreads, $S(t)$ decreases rapidly, however,
$S(t)$ increases slowly by the term of $\delta R$ after the end of infection
spreading. When $\beta S(t)$ is larger than $\gamma$, the infection spreading
starts again.

Figure 2: (a) Time evolutions of $I(t)$ at $\delta=0.0005$, $N=5$, $\beta=2$,
and $\gamma=4$. The initial condition is $S(0)=4.999$ and $I(0)=0.001$. (b)
Time interval between the second and third peaks of $I(t)$ for several
$\delta$’s in the double-logarithmic scale. The dashed line is $0.505/\delta$.
(c) Decay constant $\kappa$ obtained by the fitting of the first seven peaks
of $I(t)$ for several $\delta$’s. The dashed line is $1.644\delta$.

Figure 2(a) shows the time evolution of $I(t)$ at $N=5$, $\beta=2$,
$\gamma=4$, and $\delta=0.0005$ in a semi-logarithmic scale. The time
intervals between the two successive peaks of $I(t)$ decrease with $t$. The
dashed line is $1.11e^{-\kappa t}$ with $\kappa=0.000839$. That is, the peak
amplitude of $I(t)$ decays exponentially initially, although $I(t)$ approaches
a constant value at $t=\infty$. Figure 2(b) shows the time interval between
the second and third peaks of $I(t)$ for several $\delta$’s in the double-
logarithmic scale. (The first peak appears at $t\simeq 0$.) The dashed line is
$0.505/\delta$. The time intervals increase as $\delta$ approaches zero.
Figure 2(c) shows the decay constant $\kappa$ obtained by the fitting of the
first seven peaks of $I(t)$ for several $\delta$’s. The dashed line is
$1.644\delta$. As $\delta$ approaches 0, the decay constant decreases linearly
with $\delta$.

We can study the repetitive infection spreading with a map if $\delta$ is
sufficiently small. In the original SIR model Eq. (1), a quantity
$Q=S+I-(\gamma/\beta)\log S$ takes a constant value in the time evolution.
After the end of infection, $I(\infty)=0$ and $S(\infty)$ satisfies

$S(\infty)-(\gamma/\beta)\log S(\infty)=S(0)-(\gamma/\beta)\log S(0),$ (5)

if $I(0)$ is sufficiently small. From this relation, the final susceptible
population $S(\infty)$ is numerically evaluated from the initial population
$S(0)$. Similarly, it can be shown that $I$ takes a peak value

$S(0)-(\gamma/\beta)\log S(0)-\gamma/\beta+(\gamma/\beta)\log(\gamma/\beta)$

when $\partial I/\partial S=\partial\\{Q-S+(\gamma/\beta)\log S\\}/\partial
S=0$ or $S=\gamma/\beta$ by the conservation law of $Q=S+I-(\gamma/\beta)\log
S$.

Next, we consider the dynamics in the long timescale. We assume that the $n$th
infection spreading starts at $t_{n0}$ and ends at $t_{n\infty}$. The
susceptible population $S(t_{n\infty})$ just after the $n$th infection
spreading is determined by the susceptible population $S(t_{n0})$ just before
the $n$th infection spreading by Eq. (5), where $S(t_{n0})$ and
$S(t_{n\infty})$ correspond respectively to $S(0)$ and $S(\infty)$ in Eq. (5).
If $S(t_{n\infty})$ is determined by Eq. (5) from $S(t_{n0})$,
$I(t_{n\infty})=I(t_{n0})$ is satisfied since $Q=S+I-(\gamma/\beta)\log S$
takes the same value at $t_{n0}$ and $t_{n\infty}$. Hereafter, we assume
approximately $t_{n0}=t_{n\infty}=t_{n}$, since $t_{n\infty}-t_{n0}$ is
considered to be negligible in the long-time dynamics as shown in Fig. 2(a).
After the end time $t_{n\infty}$ of the $n$th infection spreading, $S(t)$ is
assumed to obey the equation:

$\frac{dS}{dt}=\delta R=\delta(N-S),$ (6)

because $I$ is almost zero after the infection. The initial value is
$S(t_{n\infty})$ in the time evolution of Eq. (6) for $t_{n}<t<t_{n+1}$. The
linear differential equation Eq. (6) can be solved as
$S(t)=N-\\{N-S(t_{n\infty})\\}e^{-\delta(t-t_{n})}$. For $t_{n}<t<t_{n+1}$,
$I(t)$ changes as

$I(t)=I(t_{n})e^{\int_{t_{n}}^{t}(\beta S(t)-\gamma)dt}.$ (7)

$I(t)$ initially decreases since $S(t)$ satisfies $S(t)<\gamma/\beta$ at
$t\simeq t_{n}$. However, $S(t)$ increases with time, and then $I(t)$
increases after the time when $S(t)$ satisfies $\beta S(t)=\gamma$. The
$(n+1)$th infection spreading sets in at the time $t=t_{n+1}$ satisfying

$\int_{t_{n}}^{t}(\beta
S(t)-\gamma)dt=\beta\left\\{N(t-t_{n})-(N-S(t_{n\infty}))\frac{1-e^{-\delta(t-t_{n})}}{\delta}\right\\}-\gamma(t-t_{n})=0.$
(8)

At $t=t_{n+1}=t_{(n+1)0}$, the small but non-negligible infection population
$I(t)$ recovers to the initial value $I(t_{n\infty})$, i.e.,
$I(t_{(n+1)0})=I(t_{n\infty})$. The interval $\Delta t$ between the $n$th and
$(n+1)$th infection spreading is numerically evaluated from the relation

$\Delta t=\\{N-S(t_{n\infty})\\}\frac{1-\exp(-\delta\Delta
t)}{N\delta}-\frac{\gamma}{N\beta}\Delta t.$ (9)

Because that $\beta S(t_{(n+1)0})-\gamma$ takes a positive value at the start
time $t=t_{(n+1)0}$ of the $(n+1)$th infection spreading, the infection
spreads exponentially again.

Thus, the times of the infection spreading $t_{n}$, the time evolutions of
$S(t)$ for $t_{n}<t<t_{n+1}$, and the peak value of $I$ at the $n$th infection
spreading can be successively evaluated.

Figure 3: (a) Comparison of the direct numerical simulation (solid line) of
$S(t)$ with the theoretical approximation (dashed line) at $N=5$, $\beta=2$,
$\gamma=4$, and $\delta=0.001$. The initial condition is $S(0)=4.999$ and
$I(0)=0.001$. (b) Comparison of the direct numerical simulation (solid line)
of $I(t)$ with the theoretical approximation (dashed line) at the same
parameter set.

Figure 3(a) compares the direct numerical simulation (solid line) of $S(t)$
and the theoretical approximation (dashed line) at $N=5$, $\beta=2$,
$\gamma=4$, and $\delta=0.001$. The initial condition is $S(0)=4.999$ and
$I(0)=0.001$. Figure 3(b) compares the direct numerical simulation (solid
line) of $I(t)$ and the theoretical approximation (dashed line) at the same
parameter set. The dotted line shows the exponential decay $\exp(-0.00167t)$.
Good agreement is seen. However, the timing of the infection spreading
deviates for $t>2500$ partly because the time interval between the start and
end of each infection spreading is neglected in the above approximation.

##  III Nonlinear dynamics in the simplified SIRS model with periodic
infection rate

In general, the infection rate can change in time. The SIR model considering
the seasonality has been studied by many authors. Complicated dynamics
including chaos was found in the SIR model with the time-dependent infection
rate Aron ; Rand . In this section, we consider the simplified SIRS model with
the time-dependent infection rate. We will show that the generalized model can
be analyzed with the mapping method. The model equation is written as

$\displaystyle\frac{dS}{dt}$ $\displaystyle=$
$\displaystyle-(\beta_{0}+\beta_{1}\sin\omega t)SI+\delta(N-S-I),$
$\displaystyle\frac{dI}{dt}$ $\displaystyle=$
$\displaystyle(\beta_{0}+\beta_{1}\sin\omega t)SI-\gamma I.$ (10)

Here, the infection rate $\beta$ is assumed to change periodically in time as
$\beta(t)=\beta_{0}+\beta_{1}\sin\omega t$. Figures 4(a), (b), (c), and (d)
show the time evolutions of $I(t)$ at $\beta_{1}=0.1$, $0.2$, $0.21$, and
$0.2135$ for $\beta_{0}=2$, $\gamma=4$, $\omega=2\pi/300$, $\delta=0.001$, and
$N=5$. The initial condition is $S(0)=4.999$ and $I(0)=0.001$. Owing to the
periodic change of the infection rate, the infection spreading occurs
periodically at $\beta_{1}=0.1$. That is, the spike-like infection spreading
repeats infinitely, although $I(t)$ approaches a constant value at $t=\infty$
in case of $\beta_{1}=0$. As $\beta_{1}$ increases, more complicated nonlinear
phenomena such as the period-doubling bifurcation and chaos appear. A period-2
solution appears in Fig. 4(b), and a period-4 solution appears in Fig. 4(c).
The chaotic time evolution appears in Fig. 4(d). We have calculated the
Lyapunov exponent of Eq. (10) for $\beta_{0}=2$, $\gamma=4$,
$\omega=2\pi/300$, $\delta=0.001$, and $N=5$. Figure 5 shows the Lyapunov
exponent of Eq. (10) as a function of $\beta_{1}$. The period-doubling
bifurcation occurs at $\beta_{1}\simeq 0.2054$ and the period-2 solution
appears. The Lyapunov exponent takes 0 at the critical parameter of the
period-doubling bifurcation. The period-doubling bifurcation from the period-2
solution to the period-4 solution occurs at $\beta_{1}\simeq 0.2099$, and the
period-doubling bifurcation from the period-4 solution to the period-8
solution occurs at $\beta_{1}\simeq 0.21085$. The chaos appears for
$\beta_{1}>0.2111$ where the Lyapunov exponent takes a positive value.

Figure 4: Time evolutions of $I(t)$ at (a) $\beta_{1}=0.1$, (b) $0.2$, (c)
$0.21$, and (d) $0.2135$ for $\beta_{0}=2$, $\gamma=4$, $\omega=2\pi/300$,
$\delta=0.001$, and $N=5$ in Eq. (10). Figure 5: Lyapunov exponent of Eq. (10)
as a function of $\beta_{1}$ for $\beta_{0}=2$, $\gamma=4$, $\omega=2\pi/300$,
$\delta=0.001$, and $N=5$.

When $\delta$ is sufficiently small, the approximation in the previous section
can be applied. $S(t_{n\infty})$ at the final time of the $n$th infection
spreading is determined by the initial value $S(t_{n0})$ of the $n$th
infection spreading from the equation:

$S(t_{n\infty})-(\gamma_{0}+\gamma_{1}\sin\omega t_{n})/\beta\log
S(t_{n\infty})=S(t_{n0})-(\gamma_{0}+\gamma_{1}\sin\omega t_{n})/\beta\log
S(t_{n0}).$ (11)

$S(t)$ increases as $S(t)=N-\\{N-S(t_{n\infty})\\}e^{-\delta(t-t_{n})}$ for
$t_{n}<t<t_{n+1}$, and $I(t)$ changes as

$I(t)=I(t_{n})e^{f(t)},$ (12)

where

$f(t)=\int_{t_{n}}^{t}\\{(\beta_{0}+\beta_{1}\sin\omega t)S(t)-\gamma\\}dt.$

$f(t)$ can be explicitly expressed as

$\displaystyle f(t)$ $\displaystyle=$
$\displaystyle\beta_{0}N(t-t_{n})+\beta_{1}N\omega(-\cos\omega t+\cos\omega
t_{n})-\beta_{0}\\{N-S(t_{n\infty})\\}\frac{1-e^{-\delta(t-t_{n})}}{\delta}$
$\displaystyle+$
$\displaystyle\frac{e^{-\delta(t-t_{n})}}{\delta^{2}+\omega^{2}}(-\delta\sin\omega
t-\omega\cos\omega t)-\frac{1}{\delta^{2}+\omega^{2}}(-\delta\sin\omega
t_{n}-\omega\cos\omega t_{n})-\gamma(t-t_{n})$

The interval $\Delta t$ between the $n$th and $(n+1)$th infection spreading is
numerically evaluated from $f(t_{n+1})=0$. The initial value $S(t_{(n+1)0})$
at the $(n+1)th$ infection spreading is given by
$N-\\{N-S(t_{n\infty})\\}e^{-\delta(t_{n+1}-t_{n})}$, and the peak value of
$I$ at the $(n+1)$th infection spreading is evaluated at

$S(t_{(n+1)0})-(\gamma/\beta)\log
S(t_{(n+1)0})-\gamma/\beta+(\gamma/\beta)\log(\gamma/\beta),$

where $\beta=\beta_{0}+\beta_{1}\sin\omega t_{n+1}$. The map from $S(t_{n0})$
to $S(t_{(n+1)0})$ is determined by these relations. Figures 6(a), (b), (c),
and (d) show the time evolutions of $I(t_{n})$ at $\beta_{1}=0.1$, 0.2, 0.21,
and 0.1235 obtained by the map. The other parameters are set to be
$\beta_{0}=2$, $\gamma=4$, $\omega=2\pi/300$, $\delta=0.001$, and $N=5$. The
period-doubling bifurcation and chaos are reproduced in this map.

Figure 6: Time evolutions of $I(t)$ at (a) $\beta_{1}=0.1$, (b) $0.2$, (c)
$0.21$, and (d) $0.2135$ for $\beta_{0}=2$, $\gamma=4$, $\omega=2\pi/300$,
$\delta=0.001$, and $N=5$ by the approximate map.

##  IV Directed evolution in the coupled SIRS model with mutation on a one-
dimensional lattice

In the COVID-19 pandemic, new types of viruses such as $\beta$, $\delta$, and
$o$ appeared successively. It is known that infectivity tends to increase, and
virulence tends to decrease in the process of mutation and evolution. The
evolution of infectivity and virulence has been studied by many authors
Anderson ; May2 ; Sasaki . We propose a simple model equation to study the
mutation and evolution process using the simplified SIRS model. It has a form
of coupled equations of the simplified SIRS models. The population infected by
the $i$th type of virus is expressed as $I_{i}(t)$ and the corresponding
recovered population is denoted by $R_{i}(t)$. The susceptible population has
no type and is expressed as $S(t)$. The mutation occurs in the genetic
sequence, that is, the base substitution in DNA is a fundamental process of
the mutation. Although the diffusion in the hyper-space of the base sequence
is more realistic, we consider a one-dimensional lattice model for the sake of
simplicity. That is, we assume that the $(i+1)$th and $(i-1)$th types appear
from the $i$th type with mutation rate $D$. In other words, the population
$I_{i}$ increases by $D(I_{i+1}+I_{i-1})$ owing to the mutation from the
$(i+1)$th and $(i-1)$th types and decreases by $2DI_{i}$ due to the mutation
of the $i$th type. The infection rate $\beta$ depends on the type and is
expressed as $\beta_{i}$. As a simple model, $\beta_{i}$ is assumed to be
$\beta_{i}=\beta_{0}+\beta_{1}(i-L/2)$, where $L$ is the total number of
types. That is, the infection rate is assumed to be higher for larger $i$. The
initial type is $L/2$ and the infection rate is $\beta_{0}$. Thus, the model
equation is written as

$\displaystyle\frac{dS}{dt}$ $\displaystyle=$
$\displaystyle-\sum\beta_{i}SI_{i}+\delta\sum R_{i},$
$\displaystyle\frac{dI_{i}}{dt}$ $\displaystyle=$
$\displaystyle\beta_{i}SI_{i}-\gamma I_{i}+D(I_{i+1}-2I_{i}+I_{i-1}),$
$\displaystyle\frac{dR_{i}}{dt}$ $\displaystyle=$ $\displaystyle\gamma
I_{i}-\delta R_{i}.$ (14)

The total population $N=S+\sum I_{i}+\sum R_{i}$ is conserved. The boundary
conditions are $I_{0}=I_{1}$ and $I_{L+1}=I_{L}$. Equation (14) takes the form
of coupled equations of the simplified SIRS models which can also be
interpreted as a discrete-type reaction-diffusion equation.

Figure 7: (a) Time evolution of the profile $I_{i}$ by Eq. (14) for
$\beta_{0}=2$,$\beta_{1}=0.00005$, $\gamma=4$, $\delta=0.001$, $D=0.1$, and
$L=10000$. The initial condition is $S=4.999$ and $I_{L/2}=0.001$. (b) Time
evolution of the average displacement $X$. (c) Time evolution of the width
$W$.

Figure 7(a) is the time evolution of the profile $I_{i}$ for $\beta_{0}=2$,
$\beta_{1}=0.00005$, $\gamma=4$, $\delta=0.001$, $D=0.1$, and $L=10000$. The
initial condition is $S=4.999$, $I_{L/2}=0.001$, $I_{i}=0$ for $i\neq L/2$,
and $R_{i}=0$. A pulse-like wave propagates in the right direction. To
characterize the pulse-like wave, we calculated the average position $\langle
i(t)\rangle$ and the width $W(t)$:

$\langle i(t)\rangle=\sum iI_{i}/\sum I_{i},\;\;W(t)=\\{(\sum i^{2}I_{i}/\sum
I_{i})^{2}-\langle i(t)\rangle^{2}\\}^{1/2}.$

Figures 7(b) and 7(c) shows the time evolution of the average displacement
$X(t)=\langle i(t)\rangle-\langle i(0)\rangle$ and width $W(t)$ of the pulse-
like wave. The dashed line in Fig. 7(c) is $\sqrt{2Dt}$. That is, the width
$W$ grows as $t^{1/2}$ initially, $W$ decreases around $t=15000$, and tends to
saturate for $t>15000$. The displacement $X$ increases as $t^{2}$ initially,
but $X$ seems to change in proportion to $t$ for $t>15000$. For sufficiently
large $t$, $X$ takes a constant value close to $L/2$, that is, the evolution
stops at the boundary. This simple model reproduces the evolution toward the
type of higher infection rate. We do not understand the reason of the non-
smooth dynamical behavior of $W$ as shown in Fig. 7(b) yet. The rapid increase
of $X$ might cause the non-smooth behavior.

In the stationary state of Eq. (3), $S$ takes a value $\gamma/\beta$. Under
this approximation, $I_{i}$ obeys an equation:

$\frac{dI_{i}}{dt}=(\beta_{i}\gamma/\beta-\gamma)I_{i}+D(I_{i+1}-2I_{i}+I_{i-1}).$

However, the total sum $\sum I_{i}$ is not conserved in this model equation. A
simple model satisfying the conservation law is

$\frac{dI_{i}}{dt}=(\beta_{i}\gamma/\beta-\langle\beta_{i}\gamma/\beta\rangle)I_{i}+D(I_{i+1}-2I_{i}+I_{i-1}),$
(15)

where
$\langle\beta_{i}\gamma/\beta\rangle=\sum(\beta_{i}\gamma/\beta)I_{i}/\sum
I_{i}$ is the weighted average of the effective infection rate
$\beta_{i}\gamma/\beta$ by the weight $I_{i}$. The constant part of the
effective infection rate independent of $i$ such as $-\gamma$ does not appear
in this model equation because only the difference of the effective infection
rate from the weighted average determines the dynamics. The total sum of
$I_{i}$ is conserved in the time evolution by Eq. (15). The first term of Eq.
(15) implies that the population $I_{i}$ with the effective infection rate
larger than the weighted average grows in time. This type of evolution
equation was proposed by one of the authors Sakaguchi2 . A similar pulse-like
wave appears even in Eq. (15).

Figure 8: Time evolutions of (a) the average displacement $X$ and (b) width
$W$ by Eq. (15) for $\gamma/\beta=2$, $\beta_{i}=2+0.00005(i-L/2)$, $D=0.1$,
and $L=10000$.

Figures 8(a) and (b) are the time evolutions of the average displacement $X$
and width $W$ for $\gamma/\beta=2$, $\beta_{i}=2+0.00005(i-L/2)$, $D=0.1$, and
$L=10000$ for Eq. (15). The initial condition is $I_{L/10}=0.001$ and
$I_{i}=0$ for $i\neq L/10$. Time evolutions of $X$ and $W$ shown in Figs. 8(a)
and (b) are similar to those shown in Figs. 7(b) and (c). The dashed line in
Fig. 8(b) is $\sqrt{2Dt}$. The time evolution of $W$ begins to deviate from
the square root law near $t=5000$. The average displacement initially grows as
$t^{2}$ but changes in proportion to $t$ for $t>15000$.

The continuum approximation of Eq. (15) yields

$\frac{\partial I}{\partial
t}=\\{\beta(x)\gamma/\beta-\langle\beta(x)\gamma/\beta\rangle\\}I+D\frac{\partial^{2}I}{\partial
x^{2}}.$ (16)

Here, $\beta(x)=\beta_{0}+\beta_{1}(x-L/2)$. There is a Gaussian-type solution
to Eq. (16): $I(x)=A\sqrt{\alpha(t)/\pi}\exp[-\alpha(t)\\{x-X(t)\\}^{2}]$,
where $\alpha(t)$ and $X(t)$ satisfy

$\frac{d\alpha}{dt}=-4D\alpha^{2},\;\;\;\frac{dX}{dt}=\frac{\gamma\beta_{1}}{2\beta\alpha}.$
(17)

The solution of Eq. (17) is

$\alpha(t)=\frac{\alpha(0)}{1+4\alpha(0)Dt},X(t)=X(0)+\frac{\gamma\beta_{1}}{2\beta\alpha(0)}t+(\gamma/\beta)\beta_{1}Dt^{2},$

where $\alpha(0)$ and $X(0)$ are initial values of $\alpha(t)$ and $X(t)$. The
width $W$ of the Gaussian function is expressed with $\alpha$ as
$W=1/(2\alpha(t))$. If the initial condition is $I(x)=\delta(x)$,
$W=\sqrt{2Dt}$ and $X=(\gamma/\beta)\beta_{1}Dt^{2}$. The dashed line in Fig.
8(a) is $(\gamma/\beta)\beta_{1}Dt^{2}$ and the dashed line in Fig. 8(b) is
$W=\sqrt{2Dt}$. The dashed line in Fig. 7(b) is
$X=\gamma/(\beta_{0}+\beta_{1}X_{0})\beta_{1}Dt^{2}$ where $X_{0}$ is the
first approximation of $X$: $X_{0}=(\gamma/\beta_{0})\beta_{1}Dt^{2}$ and
$\beta$ is approximated at $\beta_{0}+\beta_{1}X_{0}$. The continuum
approximation is initially good, but the deviation is observed for $t>5000$.
Similar behavior was observed in Eq. (14) as shown in Figs. 7(b) and (c).

##  V Summary

We studied the nonlinear dynamics of two simple epidemic models based on the
SIRS model. In Sec. 2 and 3, we studied the repetitive spreading of infection
in a simplified SIRS model. If $\delta$ is zero, the model equation is reduced
to the SIR model, and the infection spreading occurs only once. Even if
$\delta$ is sufficiently small, the infection spreading occurs many times,
that is, the dynamical behavior changes drastically owing to the terms $\delta
R$. We showed that the time evolution can be reduced to a map when the decay
rate $\delta$ of acquired immune is sufficiently small. The repetitive
spreading of infection was reproduced by the map. The approximation can also
be applied to the SIRS model with the periodic variation of infection rate.
The map reproduced the period-doubling bifurcation and chaos. In Sec. 4, we
studied a simple evolution model which has a form of the coupled SIRS
equations where the mutation is represented by a discrete-type diffusion
equation. We confirmed that the strengthening of infection rate and weakening
of virulence occur in this type of simple model of evolution. By the continuum
approximation, a kind of reaction-diffusion equation for $I$ is derived. The
reaction-diffusion equation with the infection rate depending linearly on $x$
has an exact solution. The initial time evolutions in the discrete models Eqs.
(14) and (15) were well approximated at the solution of Eq. (16).

The main results in this paper are that the nonlinear dynamics and the
evolution can be reproduced respectively by the map derived from the
differential equation and the exact solution of the reaction-diffusion
equation. We studied the SIRS model and the coupled SIRS model for several
parameter sets. We would like to study the model equations further in the
future.

## References

  * (1) W. O. Kermack and A. G. McKendrick, Proc. Roy. Soc. A 115, 700 (1927).
  * (2) R. M. Anderson and R. M. May, Science 215, 1053 (1982).
  * (3) H.W. Hethcote, Math. Biosciences 28, 335 (1976), H. W. Hethcote, SIAM Rev. 42, 599 (2000).
  * (4) J. D. Murray, Mathematical Biology (Springer-Verlag, Berlin Heidelberg, 2003).
  * (5) O. Diekmann and J. A. P. Heesterbeek, Mathematical Epidemiology of Infectious Diseases (John Wiley and Sons, Chichester, 2000).
  * (6) V. Cappaso and G. Serio, Math. Bio. 42, 43 (1978).
  * (7) J. V. Nobel, Nature 250, 726 (1974).
  * (8) M. Katori and N. Konno, J. Phys. Soc. Jpn. 59, 877 (1990).
  * (9) R. Pastor-Satorras and A. Vespignani, Phys. Rev. Lett. 86, 3200 (2001).
  * (10) G. Odor, Phys. Rev. E 103, 062112 (2021).
  * (11) H. Sakaguchi and Y. Nakao, Phys. Rev. E 103, 012301 (2021).
  * (12) B. Shulgin, L. Stone, and Z. Agur, Bull. Math. Biol. 60, 1123 (1998).
  * (13) L. Laguzet and G. Turinici, Bull. Math. Biol. 77, 1955 (2015).
  * (14) L. Ferretti, C. Wymant, M. Kendall, L. Zhao, A. Nurtay, L. Abeler-Dörner, M. Parker, D. Bonsall, and C. Fraser, Science 10.1126/science.abb6936 (2020).
  * (15) S. M. Kissler, C. Tedijanto, E. Goldstein., Y. H. Grad, and M. Lipsitch, Science 369, 860 (2020).
  * (16) E. B. Are, Y. Song, J. E. Stockdale, P. Tupper and C. Colijn, J. Theor. Biol. 559, 111368 (2023).
  * (17) J. L. Aron and I. B. Schwartz, J. Theor. Biol. 110, 665 (1984).
  * (18) D. A. Rand and H. B. Wilson, Proc. Roy. Soc. Lond. B 246, 179 (1991).
  * (19) A. Dutta, Scientific Reports 12, 9332 (2022).
  * (20) R. M. May, Nature 261, 459 (1976).
  * (21) R. M. May and R. M. Anderson, Nature 280, 455 (1979).
  * (22) A. Sasaki and Y. Haraguchi, J. Mol. Biol. 51, 245 (2000).
  * (23) H. Sakaguchi, Europhys. Lett. 37, 571 (1997).