# Bayesian Spatiotemporal Modeling for Inverse Problems

Shiwei Lan Shuyi Li Mirjeta Pasha<EMAIL_ADDRESS>School of Mathematical &
Statistical Sciences, Arizona State University,
Tempe, AZ 85287

###### Abstract

Inverse problems with spatiotemporal observations are ubiquitous in scientific
studies and engineering applications. In these spatiotemporal inverse
problems, observed multivariate time series are used to infer parameters of
physical or biological interests. Traditional solutions for these problems
often ignore the spatial or temporal correlations in the data (static model),
or simply model the data summarized over time (time-averaged model). In either
case, the data information that contains the spatiotemporal interactions is
not fully utilized for parameter learning, which leads to insufficient
modeling in these problems. In this paper, we apply Bayesian models based on
spatiotemporal Gaussian processes (STGP) to the inverse problems with
spatiotemporal data and show that the spatial and temporal information
provides more effective parameter estimation and uncertainty quantification
(UQ). We demonstrate the merit of Bayesian spatiotemporal modeling for inverse
problems compared with traditional static and time-averaged approaches using a
time-dependent advection-diffusion partial different equation (PDE) and three
chaotic ordinary differential equations (ODE). We also provide theoretic
justification for the superiority of spatiotemporal modeling to fit the
trajectories even it appears cumbersome (e.g. for chaotic dynamics).

Keywords: Spatiotemporal Inverse Problems, Spatiotemporal Gaussian Process,
Chaotic Dynamics, Trajectory Fitting, Uncertainty Quantification

##  1 Introduction

Many inverse problems in science and engineering involve large scale
spatiotemporal data, typically recorded as multivariate time series. There are
examples in fluid dynamics that describes the flow of liquid (e.g. petroleum)
or gas (e.g. flame jet) [31]. Other examples include dynamical systems with
chaotic behavior prevalent in weather prediction [37], biology [35], economics
[8] etc. where small perturbation of the initial condition could lead to large
deviation from what is observed/calculated in time. The goal of such inverse
problems is to recover the parameters from given observations and knowledge of
the underlying physics. The spatiotemporal information is crucial and should
be respected when considering proper statistical models for parameter
learning. This is not only of interest in statistics, but also beneficial for
practical applications of physics and biology to obtain inverse solutions and
UQ more effectively.

Traditional methods for these spatiotemporal inverse problems often ignore the
time dependence in the data for a simplified solution [54, 11, 34]. They
either treat the observed time series statically as independent identically
distributed (i.i.d.) observations across times [54, 34] (hence we refer to it
as “static” model), or summarize them by taking time average or higher order
moments [40, 11, 27] (referred as “time-averaged” approach). The former is
prevalent in Bayesian inverse problems with time series observations [34]. The
latter is especially common in parameter learning of chaotic dynamics, e.g.
Lorenz systems [37, 11], due to their sensitivity to the initial conditions
and the system parameters, which in turn causes a rough landscape of the
objective function. In both scenarios, the spatiotemporal information is not
fully integrated into the statistical modeling.

In this paper, we propose to apply Bayesian methods based on GP to the inverse
problems with sptiotemporal data to account for the space-time inter-
dependence. This leads to fitting the whole trajectories of the observed data,
rather than their statistical summaries, with elaborated models. More
specifically, we use the STGP model [14] to fit the observed multivariate time
series in comparison with the static or the time-averaged (for summarized
data) models. Theoretically, we justify why the STGP model should be preferred
to by investigating their Fisher information, which can be used as a
measurement of convexity: STGP renders a more convex likelihood than the other
two models and leads to an easier learning of the parameters. We also
demonstrate in numerical experiments (Section 4) that the STGP model yields
parameter estimates closer to the truth with smaller observation window
required, and also provides more reasonable UQ results. Note this implies
faster convergence (future work) by the STGP model, which is computationally
important because complex ODE/PDE systems are usually expensive to solve.

Spatiotemporal reasoning/modeling was introduced to inverse problems. However,
it was either qualitatively applied to specific domains such as functional
magnetic resonance imaging (fMRI) [56], electroencephalography (EEG) [52] and
electrocardiography (ECG) [51], or to a simplified Gauss-linear problem [36,
42, 12, 58]. Spatiotemporal information was also used to construct prior [62]
and regularization [59, 44], or to reduce the number of parameters [17].
However, none of them formulates the spatiotemporal modeling in the general
framework of Bayesian inverse problems with spatiotemporal observations. We
summarize the main contributions of this work as follows:

  * •

It formulates a Bayesian modeling framework for inverse problems with
spatiotemporal data that includes traditional static and time-averaged
methods;

  * •

It provides a theoretical justification on why the STGP model is preferable in
the spatiotemporal inverse problems;

  * •

It numerically demonstrates the advantage of the STGP model in parameter
learning and UQ.

The rest of the paper is organized as follows: Section 2 reviews the
background of Bayesian UQ for inverse problems, with a particular framework
named _Calibration-Emulation-Sampling (CES)_ [11, 34]. In Section 3 we
generalize the problem setup to include spatiotemporal observations and
compare the STGP model (Section 3.3) with the static model (Section 3.1) and
the time-averaged model (Section 3.2). We prove in theorems 3.1 and 3.2 that
the STGP model can have more convex likelihood than the static and the time-
averaged models. Then in Section 4 we demonstrate the advantage of the STGP
model over the other two traditional approaches with inverse problems
involving an advection-diffusion equation and three chaotic dynamics. Finally
we conclude with some discussions on future directions in Section 5.

##  2 Background: Bayesian UQ for Inverse Problems

In many inverse problems, we are interested in finding an unknown parameter,
$u$ (which could be a function or a vector), given the observed data, $y$. The
parameter $u$ usually appears as a quantity of interest in the inverse
problem, e.g. the initial condition of a time-dependent advection-diffusion
problem (Section 4.1) or the coefficient vector in the chaotic dynamics
(Section 4.2). Let ${\mathbb{X}}$ and ${\mathbb{Y}}$ be two separable Hilbert
spaces. A forward mapping ${\mathcal{G}}:{\mathbb{X}}\rightarrow{\mathbb{Y}}$
from the parameter space ${\mathbb{X}}$ to the data space ${\mathbb{Y}}$ (e.g.
${\mathbb{Y}}=\mathbb{R}^{m}$ for $m\geq 1$) connects $u\in{\mathbb{X}}$ to
$y\in{\mathbb{Y}}$ as follows:

$y={\mathcal{G}}(u)+\eta,\qquad\eta\sim{\mathcal{N}}(0,\Gamma)$ (1)

We can define the following potential function (negative log-likelihood),
$\Phi:{\mathbb{X}}\times{\mathbb{Y}}\to\mathbb{R}$, often with
$\Gamma=\sigma^{2}I$:

$\Phi(u;y)=\frac{1}{2}\|y-{\mathcal{G}}(u)\|^{2}_{\Gamma}=\frac{1}{2}\langle
y-{\mathcal{G}}(u),\Gamma^{-1}(y-{\mathcal{G}}(u))\rangle$ (2)

The forward mapping ${\mathcal{G}}$ represents physical laws usually expressed
as large and complex ODE/PDE systems that could be highly non-linear.
Therefore repeated evaluations of $\Phi(u;y)$ (and hence ${\mathcal{G}}(u)$)
are expensive for different $u$’s.

In the Bayesian setting, a prior measure $\mu_{0}$ is imposed on $u$,
independent of $\eta$. For example, we could assume a Gaussian prior
$\mu_{0}=\mathcal{N}(0,\mathcal{C})$ with the covariance $\mathcal{C}$ being a
positive, self-adjoint and trace-class operator on ${\mathbb{X}}$. Then we can
obtain the posterior of $u$, denoted as $\mu(u|y)$, using Bayes’ theorem [53,
15]:

$\frac{d\mu}{d\mu_{0}}(u)=\frac{1}{Z}\,\exp(-\Phi(u;y))\ ,\quad\textrm{if}\
0<Z:=\int_{{\mathbb{X}}}\exp(-\Phi(u;y))\mu_{0}(du)<+\infty\ .$ (3)

Bayesian UQ for inverse problems involves learning the posterior distribution
$\mu(du)$ which often exhibits strongly non-Gaussian behavior, posing
significant challenges for efficient inference methods such as Markov Chain
Monte Carlo (MCMC).

There are three urging computational challenges in the Bayesian UQ for inverse
problems: 1) intensive computation for likelihood evaluations, which requires
expensive solving of forward problems; 2) complex (non-Gaussian) posterior
distributions; and 3) high dimensionality of the discretized parameter (still
denoted as $u$ when there is no confusion from the context). The latter makes
the first two more difficult in the sense that high dimensionality not only
makes the forward solutions more expensive, but also challenges the robustness
of sampling algorithms. To address these challenges, an approximate inference
framework named Calibration-Emulation-Sampling (CES) has recently been
proposed by [11] and developed by [34]. It consists of the following three
stages:

  1. 1.

Calibration: using optimization-based (ensemble Kalman) algorithms to obtain
parameter estimation and collect expensive forward evaluations for the
emulation step;

  2. 2.

Emulation: recycling forward evaluations from the calibration stage to build
an emulator for sampling;

  3. 3.

Sampling: sampling the posterior approximately based on the emulator, which is
much cheaper than the original forward mapping.

CES calibrates the model with ensemble Kalman (EnK) methods [19, 20]. Two
algorithms, ensemble Kalman inversion (EKI) [48, 22] and ensemble Kalman
sampler (EKS) [22, 23], evolve $J$ ensemble particles
$\\{u^{(j)}\\}_{j=1}^{J}$ according to the following equations respectively:

$\displaystyle\textrm{EKI}:$ $\displaystyle\frac{du^{(j)}}{dt}$
$\displaystyle=\frac{1}{J}\sum_{k=1}^{J}\left\langle{\mathcal{G}}(u^{(k)})-\overline{\mathcal{G}},y-{\mathcal{G}}(u^{(j)})+\sqrt{\varSigma}\frac{dW^{(j)}}{dt}\right\rangle_{\Gamma}(u^{(k)}-\overline{u})$
(4a) $\displaystyle\textrm{EKS}:$ $\displaystyle\frac{du^{(j)}}{dt}$
$\displaystyle=\frac{1}{J}\sum_{k=1}^{J}\left\langle{\mathcal{G}}(u^{(k)})-\overline{\mathcal{G}},y-{\mathcal{G}}(u^{(j)})\right\rangle_{\Gamma}(u^{(k)}-\overline{u})-C(u){\mathcal{C}}^{-1}u^{(j)}+\sqrt{2C(u)}\frac{dW^{(j)}}{dt}$
(4b)

where $\overline{u}:=\frac{1}{J}\sum_{j=1}^{J}u^{(j)}$,
$\overline{\mathcal{G}}:=\frac{1}{J}\sum_{j=1}^{J}{\mathcal{G}}(u^{(j)})$,
$\varSigma=0$ or $\Gamma$, $\\{W^{(j)}\\}$ are independent cylindrical
Brownian motions on ${\mathbb{Y}}$, and
$C(u):=\frac{1}{J}\sum_{j=1}^{J}(u^{(j)}-\overline{u})\otimes(u^{(j)}-\overline{u})$.
Implemented in parallel, EnK algorithms converge quickly to the optimal
parameter with a few (usually hundreds of) ensembles without explicit
calculation of gradients. However, due to the collapse of ensembles [48, 49,
16, 9], the sample variance given by $\\{u^{(j)}\\}_{j=1}^{J}$ tends to
underestimate the actual uncertainty [see Figure 1 in 34].

CES recovers the proper uncertainty by running sampling algorithms based on
emulators ${\mathcal{G}}^{e}:{\mathbb{X}}\rightarrow{\mathbb{Y}}$ trained on
data $\\{u_{n}^{(j)},{\mathcal{G}}(u_{n}^{(j)})\\}_{j=1,n=0}^{J,N}$ that have
been collected in the calibration stage. The emulator can be GP [11] or neural
network (NN), e.g. convolutional NN (CNN) [34], with the latter being
preferred to the former for its computational efficiency and no need to design
an optimal training set with controlled size.

Once the emulator is built, CES approximately samples from the posterior with
dimension-independent MCMC algorithms based on the emulated likelihood
$\Phi^{e}$ and its gradient $D\Phi^{e}$ (defined by substituting
${\mathcal{G}}$ with ${\mathcal{G}}^{e}$ in (2)) at much lower computational
cost. A class of dimension-independent algorithms – including _preconditioned
Crank-Nicolson (pCN)_ [13], _infinite-dimensional MALA ( $\infty$-MALA)_ [6],
_infinite-dimensional HMC ( $\infty$-HMC)_ [3], and _infinite-dimensional
manifold MALA ( $\infty$-mMALA) [4] and HMC ($\infty$-mHMC) [5]_ – are used to
overcome the deteriorating mixing time of traditional Metropolis-Hastings
algorithms as the dimension of parameter space increases. These algorithms can
all be derived from the following Hamiltonian dynamics on manifold
$({\mathbb{X}},\mathcal{K}(u))$:

$\frac{d^{2}u}{dt^{2}}+\mathcal{K}(u)\,\big{\\{}\,\mathcal{C}^{-1}u+\alpha\nabla\Phi(u)\big{\\}}=0,\quad\left.\left(v:=\frac{du}{dt}\right)\right|_{t=0}\sim\mathcal{N}(0,\mathcal{K}(u))\
.$ (5)

where $\nabla\Phi(u)$ denotes the Fréchet derivative of $\Phi$, and
$\mathcal{K}(u)^{-1}={\mathcal{C}}^{-1}+\beta\mathcal{H}(u)$ with
$\mathcal{H}(u)$ being chosen as Hessian. They are implemented by numerically
simulating (5) for $I$ steps to generate a proposal $u^{\prime}$ that is
accepted with certain probability. We have $\infty$-mHMC for $\alpha\equiv
1,\beta\equiv 1$. With $\alpha\equiv 1,\beta\equiv 0$, it reduces to
$\infty$-HMC. If we let $I=1$, then $\infty$-mHMC reduces to $\infty$-mMALA,
which becomes pCN further with $\alpha=0$.

##  3 Spatiotemporal Inverse Problems (STIP)

When the observations are taken from a spatiotemporal process, $y({\bf x},t)$,
simple Gaussian likelihood function as (2) with $\Gamma=\sigma^{2}I$, for
example, may not be sufficient to describe the space-time interactions. To
address this issue, we propose to rewrite the data model (1) in terms of a GP
with spatiotemporal kernel $\Gamma({\bf x},t)$:

$y({\bf x},t)={\mathcal{G}}(u)({\bf x},t)+\eta({\bf x},t),\qquad\eta({\bf
x},t)\sim\mathcal{GP}(0,\Gamma({\bf x},t))$ (6)

In practice, the forward model often involves time-dependent PDE, e.g. heat
equation and Navier–Stokes equations. Therefore, it is crucial to allow for
the spatiotemporal correlations in the statistical analysis of such inverse
problems. Compared to (1), model (6) offers a more appropriate definition of
the likelihood by incorporating the spatiotemporal structures in the data.

Note, the proposed general framework (6) also includes many existing
statistical models as special cases. For example, if we define the forward map
based on some covariates, ${\bf X}({\bf x},t)$, ${\mathcal{G}}(\beta)({\bf
x},t)=f(\beta,{\bf X}({\bf x},t))$, e.g. $f(\beta,{\bf X}({\bf x},t))={\bf
X}({\bf x},t)\beta({\bf x},t)$, then (6) is simply a regression model. If we
set ${\mathcal{G}}(u)({\bf x},t)={\bf L}({\bf x},t)u({\bf x},t)$ with loading
matrix ${\bf L}({\bf x},t)$, then (6) becomes a latent factor model.

In the following, we will introduce the static (Section 3.1) and the time-
averaged (Section 3.2) models and unify them in the framework of STGP model
(Section 3.3). For the convenience of exposition, we fix some notations in the
following. Denote ${\bf X}:=\\{{\bf x}_{i}\\}_{i=1}^{I}$, ${\bf
t}:=\\{t_{j}\\}_{j=0}^{J-1}$, and ${\bf Y}:=y({\bf X},{\bf t})=\\{y({\bf
x}_{i},t_{j})\\}_{i=1,j=0}^{I,J-1}$. ${\bf C}_{*}$ is the covariance matrix of
the covariance kernel ${\mathcal{C}}_{*}$ restricted on the finite-dimensional
discrete space.

###  3.1 Static model

In the literature of Bayesian inverse problems, the noise $\eta$ is often
assumed i.i.d. over time in (6), i.e. $\eta({\bf
x},t_{j})\overset{iid}{\sim}{\mathcal{N}}(0,{\mathcal{C}}_{\bf x})$. This
leads to the following static model where the temporal correlation is ignored:

$\displaystyle y({\bf x},t)|u,\Gamma$
$\displaystyle\sim\mathcal{GP}({\mathcal{G}}(u)({\bf x},t),\Gamma({\bf x},t))$
(7) $\displaystyle\textrm{static}:\qquad\Gamma({\bf x},t)$
$\displaystyle={\mathcal{C}}_{\bf x}\otimes{\mathcal{I}}_{t}$

where ${\mathcal{I}}_{t}$ is the Dirac operator such that
${\mathcal{I}}_{t}(t,t^{\prime})=1$ only if $t=t^{\prime}$. When the spatial
dependence is also suppressed (as in the advection-diffusion example of
Section 4.1 and in [54, 34]), we have ${\mathcal{C}}_{\bf
x}=\sigma^{2}_{\varepsilon}{\mathcal{I}}_{\bf x}$.

Figure 1: Advection-diffusion inverse problem: comparing maximum _a
posteriori_ (MAP) estimates of parameter $u_{0}=u({\bf x},0)$ by the static
model (middle) and the STGP model (right) with the truth $u_{0}^{\dagger}$
(left).

Temporal correlation is disregarded in the static model (7). When there is
(spatio-)temporal effect in the residual $\eta$, the static model (7) may be
insufficient to account for the spatiotemporal relationships contained in the
data. For illustration, we consider an inverse problem involving advection-
diffusion (Section 4.1) equation [54, 34] of an evolving concentration field
$u({\bf x},t)$, e.g., temperature for heat transfer, and seek the solution to
the initial condition, $u_{0}=u({\bf x},0)$, based on spatiotemporal solutions
observed (through an observation operator ${\mathcal{O}}$) on the boundaries
of two boxes (Figure 1, left panel) for a given time period, i.e.
$y={\mathcal{O}}u({\bf x},t)+\eta,\;\eta\sim N(0,\sigma^{2}_{\eta})$. As shown
in Figure 1, the simple static model (7) used in [34] does not account for
space-time interactions hence yields the result underestimating the true
function $u_{0}^{\dagger}$ (left panel). On the contrary, the estimate by the
spatiotemporal model (16) (right panel) is much closer to the truth.

###  3.2 Time-averaged model

In many chaotic dynamics, we observe the trajectories as multivariate time
series that are very sensitive to the initial condition and the parameters.
This usually results in a complex objective function with multiple local
minima [1]. They in turn form a rough landscape of the objective and pose
extreme difficulties on parameter learning [11] (See also Figure 6). The time-
averaged approach is commonly used in the spirit of extracting sufficient
statistics from the raw data [21].

We consider the same data model as in (6) with ${\mathcal{G}}(u)$ being the
observed solution ${\bf x}(t;u,{\bf x}_{0})$ of the following chaotic dynamics
($r$-th order ODE) for a given parameter $u\in{\mathbb{R}}^{p}$:

$\dot{\bf x}:=\frac{d{\bf x}}{dt}=f(t,{\bf x},{\bf x}^{(1)},\cdots,{\bf
x}^{(r)};u),\qquad{\bf x}(0)={\bf x}_{0}\in{\mathbb{R}}^{I}$ (8)

That is, ${\mathcal{G}}(u)={\mathcal{O}}{\bf x}(t;u,{\bf x}_{0})$ with an
observation operator ${\mathcal{O}}$. At each time $t$, the observed vector
could include components of ${\bf x}$ and up to their $k$-th order
interactions for $k\geq 1$. For example, if ${\bf x}=[x_{1},\cdots,x_{I}]$, we
could include all the first and second order terms in the observation vector,
${\mathcal{O}}{\bf
x}=[x_{1},\cdots,x_{I},x_{1}^{2},x_{1}x_{2},\cdots,x_{i}x_{j},\cdots,x_{I}^{2}]$.
Because the trajectories of ${\mathcal{G}}(u)$ are usually complex, it is
often to average them over time and consider the following forward mapping
instead:

${\mathcal{G}}_{T}(u;{\bf
x}_{0}):=\frac{1}{T}\int_{t_{0}}^{t_{0}+T}{\mathcal{O}}{\bf x}(t;u,{\bf
x}_{0})dt$ (9)

where $t_{0}$ is the spin-up time and $T$ is the window length for averaging
the observed trajectories of the dynamics.

Following [11], we make the same assumption regarding the dynamical system
(8):

######  Assumption 1.

  1. 1.

For $u\in{\mathbb{X}}$, (8) has a compact attractor $\mathcal{A}$, supporting
an invariant measure $\mu(d{\bf x};u)$. The system is ergodic, and the
following limit of Law of Large Numbers (LLN) is satisfied: for ${\bf
x}_{0}\sim\mu(\cdot;u)$ fixed, with probability one,

$\lim_{T\to\infty}{\mathcal{G}}_{T}(u;{\bf
x}_{0})={\mathcal{G}}_{\infty}(u):=\int_{\mathcal{A}}{\mathcal{O}}{\bf
x}(t;u,{\bf x}_{0})\mu(d{\bf x};u)$ (10)

  2. 2.

The Central Limit Theorem (CLT) holds quantifying the ergodicity: for ${\bf
x}_{0}\sim\mu(\cdot;u)$,

${\mathcal{G}}_{T}(u;{\bf
x}_{0})\overset{\cdot}{\sim}{\mathcal{N}}\left({\mathcal{G}}_{\infty}(u),T^{-1}\varSigma(u)\right)$
(11)

The limit ${\mathcal{G}}_{\infty}(u)$ becomes independent of the initial
condition ${\bf x}_{0}$. However, the finite-time truncation in
${\mathcal{G}}_{T}(u;{\bf x}_{0})$, with different random initializations
${\bf x}_{0}$, generates random errors from the limit
${\mathcal{G}}_{\infty}(u)$, which are assumed approximately Gaussian. Assume
the data $y$ can be observed with a true parameter $u^{\dagger}$, i.e.
$y={\mathcal{G}}_{T}(u^{\dagger};{\bf x}_{0})$. The following time-averaged
model is usually adopted for the inverse problems involving chaotic dynamics
[11]:

$\displaystyle y|u,\varSigma(u)$
$\displaystyle\sim{\mathcal{N}}({\mathcal{G}}_{\infty}(u),T^{-1}\varSigma(u))$
(12) $\displaystyle\text{time-average}:\qquad T^{-1}\varSigma(u)$
$\displaystyle\approx\Gamma_{\textrm{obs}}$

where the empirical covariance $\Gamma_{\textrm{obs}}$ can be estimated with
${\mathcal{G}}_{\tau}(u;{\bf x}_{0})$ for $\tau\gg T$.

In practice, we replace ${\mathcal{G}}_{\infty}(u)$ with
${\mathcal{G}}_{T}(u;{\bf x}_{0})$ in (12) and define the potential
$\Phi_{\textrm{\tiny T}}(u)$ of parameter $u$ for the time-averaged model (12)
as follows:

$\Phi_{\textrm{\tiny T}}(u)=\frac{1}{2}\|y-{\mathcal{G}}_{T}(u;{\bf
x}_{0})\|_{\Gamma_{\textrm{obs}}}^{2}$ (13)

If we observe the trajectories (without component interaction terms, i.e.
${\mathcal{O}}{\bf x}={\bf x}$) at discrete time points ${\bf t}$ with
$t_{J-1}=t_{0}+T$, then ${\mathcal{O}}{\bf x}(t;u)$ yields multivariate time
series, denoted as ${\bf X}(u)_{I\times J}={\bf x}({\bf t};u)=[{\bf
x}(t_{0};u),\cdots,{\bf x}(t_{J-1};u)]$. Then we have

${\mathcal{G}}_{T}(u;{\bf x}_{0})=\overline{{\bf X}}(u):={\bf
X}(u)\frac{\bm{1}_{J}}{J},\quad y={\bf
X}(u^{\dagger})\frac{\bm{1}_{J}}{J},\quad\Gamma_{\textrm{obs}}={\bf
X}(u^{\dagger})\left[{\bf
I}_{J}-\frac{\bm{1}_{J}{\bm{1}}^{\mathsf{T}}_{J}}{J}\right]{{\bf
X}(u^{\dagger})}^{\mathsf{T}}$ (14)

Denote ${\bf X}_{0}={\bf X}(u)-{\bf X}(u^{\dagger})$. Therefore the potential
$\Phi_{\textrm{\tiny T}}$ becomes

$\Phi_{\textrm{\tiny
T}}(u)=\frac{1}{2}\frac{{\bm{1}}^{\mathsf{T}}_{J}}{J}{{\bf
X}}^{\mathsf{T}}_{0}\Gamma_{\textrm{obs}}^{-1}{\bf
X}_{0}\frac{\bm{1}_{J}}{J}=\frac{1}{2}\mathrm{tr}\left[\frac{\bm{1}_{J}{\bm{1}}^{\mathsf{T}}_{J}}{J^{2}}{{\bf
X}}^{\mathsf{T}}_{0}\Gamma_{\textrm{obs}}^{-1}{\bf X}_{0}\right]$ (15)

Note, averaging the trajectories over time does not ease the difficulty of
rough landscape, see for instance Figure 6 for an illustration. However, the
potential function for the following STGP model (16) is more convex around the
true values $u^{\dagger}$ compared with the time-averaged approach (12).

The aforementioned two approaches, the static model (7) and the time-averaged
model (12), can be recognized as special cases of a more general framework of
spatiotemporal modeling based on STGP, to be discussed in the following
section.

###  3.3 Spatiotemporal GP model

For the spatiotemporal data $y({\bf x},t)$ in the inverse problems, we
consider the following likelihood model based on STGP:

$\displaystyle y({\bf x},t)|u,\Gamma$
$\displaystyle\sim\mathcal{GP}({\mathcal{G}}(u)({\bf x},t),\Gamma({\bf x},t))$
(16) $\displaystyle\textrm{STGP}:\qquad\Gamma({\bf x},t)$
$\displaystyle={\mathcal{C}}_{\bf x}\otimes{\mathcal{C}}_{t}$

where ${\mathcal{C}}_{\bf x}$ and ${\mathcal{C}}_{t}$ are spatial and temporal
kernel respectively.

If we observe the process $y({\bf x},t)$ according to (16), the resulted data
matrix ${\bf Y}={\mathcal{G}}(u)({\bf X},{\bf t})$ follows the matrix normal
distribution (denoted as ‘$\mathcal{MN}$’) [24] for which we can also specify
the above-mentioned three models

$\displaystyle{\bf Y}\,|\,{\bf M},{\bf U},{\bf V}$
$\displaystyle\sim\mathcal{MN}({\bf M},{\bf U},{\bf V}),\quad{\bf
M}={\mathcal{G}}(u^{\dagger})({\bf X},{\bf t})$
$\displaystyle\textrm{static}:\qquad$ $\displaystyle{\bf U}_{\textrm{\tiny
S}}=\sigma^{2}_{\varepsilon}{\bf I}_{\bf x},\quad{\bf V}_{\textrm{\tiny
S}}={\bf I}_{t}$ (17a) $\displaystyle\textrm{time-average}:\qquad$
$\displaystyle{\bf U}_{\textrm{\tiny T}}=\Gamma_{\textrm{obs}},\quad{\bf
V}_{\textrm{\tiny T}}=J^{2}(\bm{1}_{J}{\bm{1}}^{\mathsf{T}}_{J})^{-}$ (17b)
$\displaystyle\textrm{STGP}:\qquad$ $\displaystyle{\bf U}_{\textrm{\tiny
ST}}={\bf C}_{\bf x},\quad{\bf V}_{\textrm{\tiny ST}}={\bf C}_{t}$ (17c)

where ${\bf Y}={\mathcal{O}}{\bf x}(t;u)={\bf X}(u)$ for the static model and
$M^{-}$ is the pseudo-inverse of $M$.

In all the above three models (17), we assume ${\bf Y}$ i.i.d. over $u$’s.
Denote $\Phi_{*}$ and ${\mathcal{I}}_{*}$ as potential function and Fisher
information matrix with $*$ being ‘S’ for the static model (17a), ‘T’ for the
time-averaged model (17b) and ‘ST’ for the STGP model (17c) respectively. The
following theorem compares the convexity of their likelihoods and indicates
that the STGP model (17c) with proper configuration has the advantage of
parameter learning with the most convex likelihood among the three models.

######  Theorem 3.1.

If we set the maximal eigenvalues of ${\bf C}_{\bf x}$ and ${\bf C}_{t}$ such
that $\lambda_{\max}({\bf C}_{\bf x})\lambda_{\max}({\bf
C}_{t})\leq\sigma^{2}_{\varepsilon}$, then the following inequality holds
regarding the Fisher information matrices, ${\mathcal{I}}_{\textrm{\tiny S}}$
and ${\mathcal{I}}_{\textrm{\tiny ST}}$, of the static model and the STGP
model respectively:

${\mathcal{I}}_{\textrm{\tiny ST}}(u)\geq{\mathcal{I}}_{\textrm{\tiny S}}(u)$
(18)

If we control the maximal eigenvalues of ${\bf C}_{\bf x}$ and ${\bf C}_{t}$
such that $\lambda_{\max}({\bf C}_{\bf x})\lambda_{\max}({\bf C}_{t})\leq
J\lambda_{\min}(\Gamma_{\textrm{obs}})$, then the following inequality holds
regarding the Fisher information matrices, ${\mathcal{I}}_{\textrm{\tiny T}}$
and ${\mathcal{I}}_{\textrm{\tiny ST}}$, of the time-averaged model and the
STGP model respectively:

${\mathcal{I}}_{\textrm{\tiny ST}}(u)\geq{\mathcal{I}}_{\textrm{\tiny T}}(u)$
(19)

###### Proof.

See Appendix A.

∎

The following theorem considers a special case, ${\bf C}_{\bf
x}=\Gamma_{\textrm{obs}}$, under milder condition in comparing the likelihood
convexity of the time-averaged model and the STGP model.

######  Theorem 3.2.

If we choose ${\bf C}_{\bf x}=\Gamma_{\textrm{obs}}$ and require the maximal
eigenvalue of ${\bf C}_{t}$, $\lambda_{\max}({\bf C}_{t})\leq J$, then the
following inequality holds regarding the Fisher information matrices,
${\mathcal{I}}_{\textrm{\tiny T}}$ and ${\mathcal{I}}_{\textrm{\tiny ST}}$, of
the time-averaged model and the STGP model respectively:

${\mathcal{I}}_{\textrm{\tiny ST}}(u)\geq{\mathcal{I}}_{\textrm{\tiny T}}(u)$
(20)

###### Proof.

See Appendix A.

∎

######  Remark 1.

In general, $\Phi_{*}(u)$ is not the potential of a Gaussian distribution
because of the possible non-linearity of ${\mathcal{G}}(u)$. Theorems 3.1 and
3.2 indicate that for each $u\in{\mathbb{X}}$, the STGP model can have a more
convex Gaussian proxy in the Laplace approximation.

######  Remark 2.

If we view Fisher information as a measurement of (statistical) convexity, the
above theorems 3.1 and 3.2 indicate that the STGP model can have a likelihood
more convex around the true parameter value than either the static model or
the time-averaged model does. This implies that parameter learning method
based on the STGP model could be more effective in the sense that it may
converge faster.

Often we are interested in predicting the underlying process $y({\bf x},t)$ at
future time $t_{*}$ given the spatiotemporal observations ${\bf Y}$. Based on
the STGP model (16), we could use the following posterior predicative
distribution

$p(y({\bf x},t_{*})|{\bf Y})=\int p(y({\bf x},t_{*})|u,{\bf Y})p(u|{\bf Y})du$
(21)

Denote the conditional prediction $\mathrm{E}[y({\bf x},t_{*})|u,{\bf Y}]$ as

${\mathcal{G}}^{*}(u)({\bf x},t_{*})=\underbrace{{\mathcal{G}}(u)({\bf
x},t_{*})}_{Physical}+\underbrace{\Gamma_{t_{*}{\bf t}}\Gamma_{{\bf t}{\bf
t}}^{-1}({\bf Y}-{\mathcal{G}}(u)({\bf X},{\bf t}))}_{Statistical}$ (22)

Then we predict $y({\bf x},t_{*})$ with the following predicative mean

$\mathrm{E}[y({\bf x},t_{*})|{\bf Y}]=\mathrm{E}_{u|{\bf
Y}}[\mathrm{E}_{y_{*}|u,{\bf Y}}[y({\bf x},t_{*})]]=\mathrm{E}_{u|{\bf
Y}}[{\mathcal{G}}^{*}(u)({\bf x},t_{*})]\approx\overline{{\mathcal{G}}}({\bf
x},t_{*})+\Gamma_{t_{*}{\bf t}}\Gamma_{{\bf t}{\bf t}}^{-1}({\bf
Y}-\overline{{\mathcal{G}}}({\bf X},{\bf t}))$ (23)

where $\overline{{\mathcal{G}}}({\bf
x},t_{*}):=\frac{1}{S}\sum_{s=1}^{S}{\mathcal{G}}(u^{(s)})({\bf x},t_{*})$
with $u^{(s)}\sim p(u|{\bf Y})$. And we can quantify the uncertainty using the
law of total conditional variance:

$\displaystyle\mathrm{Var}[y({\bf x},t_{*})|{\bf Y}]$
$\displaystyle=\mathrm{E}_{u|{\bf Y}}[\mathrm{Var}_{y_{*}|u,{\bf Y}}[y({\bf
x},t_{*})]]+\mathrm{Var}_{u|{\bf Y}}[\mathrm{E}_{y_{*}|u,{\bf Y}}[y({\bf
x},t_{*})]]$ (24) $\displaystyle=\Gamma_{t_{*}t_{*}}-\Gamma_{t_{*}{\bf
t}}\Gamma^{-1}_{{\bf t}{\bf t}}\Gamma_{{\bf t}t_{*}}+\mathrm{Var}_{u|{\bf
Y}}[{\mathcal{G}}^{*}(u)({\bf x},t_{*})]$
$\displaystyle\approx\Gamma_{t_{*}t_{*}}-\Gamma_{t_{*}{\bf
t}}\Gamma^{-1}_{{\bf t}{\bf t}}\Gamma_{{\bf
t}t_{*}}+s^{2}_{{\mathcal{G}}^{*}}({\bf x},t_{*})$

where $s^{2}_{{\mathcal{G}}^{*}}({\bf
x},t_{*}):=\frac{1}{S}\sum_{s=1}^{S}[{\mathcal{G}}^{*}(u^{(s)})({\bf
x},t_{*})-\overline{{\mathcal{G}}^{*}}({\bf x},t_{*})]^{2}$ with $u^{(s)}\sim
p(u|{\bf Y})$.

Assume $t_{*}\notin{\bf t}$. For static model (7), we have $\Gamma_{t_{*}{\bf
t}}=0$ thus ${\mathcal{G}}^{*}(u)({\bf x},t_{*})={\mathcal{G}}(u)({\bf
x},t_{*})$. Therefore we have the simplified results

$\mathrm{E}[y({\bf x},t_{*})|{\bf Y}]\approx\overline{{\mathcal{G}}}({\bf
x},t_{*}),\qquad\mathrm{Var}[y({\bf x},t_{*})|{\bf
Y}]\approx\sigma^{2}_{\varepsilon}+s^{2}_{{\mathcal{G}}}({\bf x},t_{*})$ (25)

This may underestimate the uncertainty compared with the more general STGP
model (16). If we are only interested in predicting the forward map
${\mathcal{G}}(u)$ to new time $t=t_{*}$, we actually have similar results

$\mathrm{E}[{\mathcal{G}}(u)({\bf x},t_{*})|{\bf
Y}]\approx\overline{{\mathcal{G}}}({\bf
x},t_{*}),\qquad\mathrm{Var}[{\mathcal{G}}(u)({\bf x},t_{*})|{\bf Y}]\approx
s^{2}_{{\mathcal{G}}}({\bf x},t_{*})$ (26)

Note all the above prediction is feasible only if we are able to solve ODE/PDE
systems to time $t_{*}$, i.e. we can evaluate ${\mathcal{G}}(u^{(s)})({\bf
x},t)$ at $t=t_{*}$. When we do not have the computer codes available for
doing so, we could model ${\mathcal{G}}(u)({\bf x},t)$ with another GP
$\mathcal{GP}(0,\Gamma^{\mathcal{G}})$ and further predict the forward
mapping:

${\mathcal{G}}(u)({\bf x},t_{*})|{\mathcal{G}}(u)({\bf X},{\bf
t})\sim{\mathcal{N}}(\Gamma^{\mathcal{G}}_{t_{*}{\bf
t}}(\Gamma^{\mathcal{G}}_{{\bf t}{\bf t}})^{-1}{\mathcal{G}}(u)({\bf X},{\bf
t}),\Gamma^{\mathcal{G}}_{t_{*}t_{*}}-\Gamma^{\mathcal{G}}_{t_{*}{\bf
t}}(\Gamma^{\mathcal{G}}_{{\bf t}{\bf t}})^{-1}\Gamma^{\mathcal{G}}_{{\bf
t}t_{*}})$ (27)

##  4 Numerical Experiments

In this section, we demonstrate the numerical advantage of spatiotemporal
modeling in parameter estimation and UQ. More specifically, we compare the
STGP model (16) with the static model (7) using an advection-diffusion inverse
problem (Section 4.1) previously considered in [54, 34] with the static
method. Then we compare the STGP model (16) with the time-averaged model (12)
using three chaotic dynamical inverse problems (Section 4.2) of which the
Lorenz problem (Section 4.2.1) was studied by [11] with the time-averaged
approach. Numerical evidences are presented to support that the STGP model
(16) is preferable to the other two models. All the computer codes are
publicly available at https://github.com/lanzithinking/Spatiotemporal-inverse-
problem.

###  4.1 Advection-diffusion inverse problem

In this section, we consider an inverse problem governed by a parabolic PDE
within the Bayesian inference framework. The underlying PDE is a time-
dependent advection-diffusion equation that can be applied to heat transfer,
air pollution, etc. The inverse problem involves inferring an unknown initial
condition $u_{0}\in L^{2}(\Omega)$ from spatiotemporal point measurements
$\\{y({\bf x}_{i},t_{j})\\}$.

(a) True initial condition (top left), and the solutions $u(s)$ at different
time points.

(b) Spatiotemporal observations at $80$ selected locations (color dots) across
different time points.

Figure 2: Advection-diffusion inverse problem.

The parameter-to-observable forward mapping
${\mathcal{G}}:u_{0}\to{\mathcal{O}}u$ maps the initial condition $u_{0}$ to
pointwise spatiotemporal observations of the concentration field $u({\bf
x},t)$ through the solution of the following advection-diffusion equation [45,
54]:

$\displaystyle u_{t}-\kappa\Delta u+{\bf v}\cdot\nabla u$
$\displaystyle=0\quad in\;\Omega\times(0,T)$ (28) $\displaystyle u(\cdot,0)$
$\displaystyle=u_{0}\quad in\;\Omega$ $\displaystyle\kappa\nabla
u\cdot\vec{n}$ $\displaystyle=0,\quad on\;\partial\Omega\times(0,T)$

where $\Omega\subset[0,1]^{2}$ is a bounded domain shown in Figure 2(a),
$\kappa=10^{-3}$ is the diffusion coefficient, and $T>0$ is the final time.
The velocity field ${\bf v}$ is computed by solving the following steady-state
Navier-Stokes equation with the side walls driving the flow [45]:

$\displaystyle-\frac{1}{\mathrm{Re}}\Delta{\bf v}+\nabla q+{\bf
v}\cdot\nabla{\bf v}$ $\displaystyle=0\quad in\;\Omega$ (29)
$\displaystyle\nabla\cdot{\bf v}$ $\displaystyle=0\quad in\;\Omega$
$\displaystyle{\bf v}$ $\displaystyle={\bf g},\quad on\;\partial\Omega$

Here, $q$ is the pressure, and $\mathrm{Re}$ is the Reynolds number, which is
set to 100 in this example. The Dirichlet boundary data ${\bf
g}\in{\mathbb{R}}^{2}$ is given by ${\bf g}={\bf e}_{2}=(0,1)$ on the left
wall of the domain, ${\bf g}=-{\bf e}_{2}$ on the right wall, and ${\bf
g}={\bf 0}$ everywhere else.

(a) Posterior mean estimates of the initial concentration field $u_{0}({\bf
x})$.

(b) Posterior standard deviation estimates of the initial concentration field
$u_{0}({\bf x})$.

Figure 3: Advection-diffusion inverse problem: comparing posterior estimates
of parameter $u_{0}$ in the static model (upper row) and the STGP model (lower
row) based on $5000$ samples by various MCMC algorithms.

We set the true initial condition
$u_{0}^{\dagger}=0.5\wedge\exp\\{-100[(x-0.35)^{2}+(y-0.7)^{2}]\\}$,
illustrated in the top left panel of Figure 2(a), which also shows a few
snapshots of solutions $u({\bf x},t)$ at other time points on a regular grid
mesh of size $61\times 61$. To obtain spatiotemporal observations $\\{y({\bf
x}_{i},t_{j})\\}$, we collect solutions $u({\bf x},t)$ solved on a refined
mesh at $I=80$ selected locations $\\{{\bf x}_{i}\\}_{i=1}^{I}$ across $J=16$
time points $\\{t_{j}\\}_{j=1}^{J}$ evenly distributed between $1$ and $4$
seconds (thus denoted as ${\mathcal{O}}u$) and inject some Gaussian noise
${\mathcal{N}}(0,\sigma^{2}_{\eta})$ such that the relative noise standard
deviation is $\sigma_{\eta}/\max{\mathcal{O}}u=0.5$, i.e.,

$y={\mathcal{G}}(u_{0}^{\dagger})={\mathcal{O}}u({\bf
x},t;u_{0}^{\dagger})+\eta,\quad\eta\sim\mathcal{N}(0,\sigma_{\eta}^{2}I_{1280})$

Figure 2(b) plots 4 snapshots of these observations at 80 locations along the
inner boundary. In the Bayesian setting, we adopt a GP prior for
$u_{0}\sim\mu_{0}=\mathcal{GP}(0,{\mathcal{C}})$ with the covariance kernel
${\mathcal{C}}=(\delta{\mathcal{I}}-\gamma\Delta)^{-2}$ defined through the
Laplace operator $\Delta$, where $\delta$ governs the variance of the prior
and $\gamma/\delta$ controls the correlation length. We set $\gamma=2$ and
$\delta=10$ in this example.

The Bayesian inverse problem involves obtaining an estimate of the initial
condition $u_{0}$ and quantifying its uncertainty based on the $80\times 16$
spatiotemporal observations. The Bayesian UQ in this example is especially
challenging not only because of its large dimensionality (3413) of spatially
discretized $u$ (Lagrange degree 1) at each time $t$, but also due to the
spatiotemporal correlations in these observations.

| Estimation | Prediction
---|---|---
Models | pCN |  $\infty$-MALA |  $\infty$-HMC | pCN |  $\infty$-MALA |  $\infty$-HMC
static | 0.83 (0.023) | 0.81 (0.011) | 0.79 (0.005) | 0.43 (0.013) | 0.4 (0.006) | 0.4 (0.003)
STGP | 0.74 (0.021) | 0.73 (0.012) | 0.73 (0.003) | 0.44 (0.068) | 0.32 (0.016) | 0.31 (0.005)
Table 1: Advection-diffusion inverse problem: comparing (i) posterior
estimates of parameter $u_{0}$ in terms of relative error of mean
$\mathrm{REM}=\frac{\|\hat{u}_{0}-u_{0}^{\dagger}\|}{\|u_{0}^{\dagger}\|}$ and
(ii) the forward predictions ${\mathcal{G}}(u)({\bf x},t_{*})$ in terms of
relative error $\frac{\|\overline{{\mathcal{G}}}({\bf
x},t_{*})-{\mathcal{G}}(u_{0}^{\dagger})({\bf
x},t_{*})\|}{\|{\mathcal{G}}(u_{0}^{\dagger})({\bf x},t_{*})\|}$ by two
likelihood models (static and STGP). Each experiment is repeated for 10 runs
of MCMC (pCN, $\infty$-MALA, and $\infty$-HMC respectively) and the numbers in
the bracket are standard deviations of these repeated experiments.

We compare two likelihood models (7) and (16). The static model (7) is
commonly used in the literature of Bayesian inverse problems [32, 54, 34].
Here the STGP model (16) is considered to better account for the
spatiotemporal relationships in the data. We estimate the variance parameter
of the joint kernel from data. The correlation length parameters are
determined ($\ell_{\bf x}=0.5$ and $\ell_{t}=0.2$) by investigating their
autocorrelations as in Figure 18. Figure 1 compares the maximum a posterior
(MAP) of the parameter $u_{0}$ by the two likelihood models (right two panels)
with the true parameter $u_{0}^{\dagger}$ (left panel). The STGP model yields
a better MAP estimate closer to the truth compared with the static model.

We also run MCMC algorithms (pCN, $\infty$-MALA, and $\infty$-HMC) to estimate
$u_{0}$. For each algorithm, we run 6000 iterations and burn in the first
1000. The remaining 5000 samples are used to obtain the posterior estimate
$\hat{u}_{0}$ (Figure 3(a)) and posterior standard deviation (Figure 3(b)).
The STGP model (16) consistently generates estimates closer to the true values
(refer to Figure 1) with smaller posterior standard deviation than the static
model (7) using various MCMC algorithms. Such improvement of parameter
estimation by the STGP model (16) is also verified by smaller relative error
of mean estimates
$\mathrm{REM}=\frac{\|\hat{u}_{0}-u_{0}^{\dagger}\|}{\|u_{0}^{\dagger}\|}$
reported in Table 1, which summarizes the results of 10 repeated experiments
with their standard deviations in the brackets.

Figure 4: Advection-diffusion inverse problem: comparing forward predictions,
$\overline{{\mathcal{G}}}({\bf x},t_{*})$, based on the static model and the
STGP model. The left panel plots the curves representing the percentage of 80
(corresponding to the selected locations) credible bands that cover the true
solution ${\mathcal{G}}(u_{0}^{\dagger})({\bf x},t_{*})$ at each time
$t_{*}\in[0,5]$. The right two panels show the predicted time series (blue
dashed and orange dot-dashed lines) along with the credible bands (shaded
regions) by the two models compared with the truth (red solid line) at two
selective locations ${\bf x}=(0.375,0.401)$ and ${\bf x}=(0.249,0.250)$. Blues
dots are observations.

Finally, we consider the forward prediction (26) over the time interval
$[0,5]$. We substitute each of the 5000 samples $\\{u^{(s)}\\}_{s=1}^{5000}$
generated by $\infty$-HMC into ${\mathcal{G}}(u^{(s)})({\bf x},t_{*})$ to
solve the advection-diffusion equation (28) for $t_{*}\in[0,5]$. We observe
each of these 5000 solutions at the 80 locations (Figure 2(b)) for $50$ points
equally spaced in $[0,5]$. Then we obtain the prediction by
$\overline{{\mathcal{G}}}({\bf x},t_{*})_{80\times
50}=\frac{1}{5000}\sum_{s=1}^{5000}{\mathcal{G}}(u^{(s)})({\bf x},t_{*})$, and
compute the relative errors in terms of the Frobenius norm of the difference
between the prediction and the true solution
${\mathcal{G}}(u_{0}^{\dagger})({\bf x},t_{*})$:
$\frac{\|\overline{{\mathcal{G}}}({\bf
x},t_{*})-{\mathcal{G}}(u_{0}^{\dagger})({\bf
x},t_{*})\|}{\|{\mathcal{G}}(u_{0}^{\dagger})({\bf x},t_{*})\|}$. Table 1
shows the STGP model (16) provides more accurate predictions with smaller
errors compared with the static model (7). Figure 4 depicts the predicted time
series $\overline{{\mathcal{G}}}({\bf x},t_{*})$ at two selective locations
based on the static (blue dashed line) and the STGP (orange dot-dashed line)
models along with their credible bands (shaded regions) compared with the
truth (red solid lines) in the two right panels. Note that with smaller
credible bands, the static model is more certain about its prediction that is
further away from the truth. While the STGP model provides wider credible
bands that cover more of the true trajectories, indicating a more appropriate
uncertainty being quantified. Therefore, on the left panel of Figure 4, the
STGP model has higher truth covering rate for its credible intervals among
these $80$ locations on most of $t_{*}\in[0,5]$. Note these models are trained
on $t\in[1,4]$, so the STGP model does not show much advantage at the
beginning but quickly outperforms the static model after $t_{*}=1$.

Figure 5: Lorenz63 dynamics: two-lobe orbits (left), chaotic solutions
(middle) and coordinates’ distributions (right).

###  4.2 Chaotic dynamical inverse problems

Chaos, refers to the behavior of a dynamical system that appears to be random
in long term even its evolution is fully determined by the initial condition.
Many physical systems are characterized by the presence of chaos that has been
extensively demonstrated [37, 30, 7]. The main challenges of analyzing chaotic
dynamical systems include the stability, the transitivity, and the sensitivity
to the initial conditions (which contributes to the seeming randomness) [18].
In the study of chaotic dynamical systems, one of the interests is determining
the essential system parameters given the observed data. In this section, we
will investigate three chaotic dynamical systems, Lorenz63 [37], Rössler [2]
and Chen [60], that can be summarized as the first-order ODE: $\dot{\bf
x}=f({\bf x};u)$. We will apply the CES framework (Section 2) to learn the
system parameter $u$ and quantify its associated uncertainty based on the
observed trajectories. We find the spatiotemporal models numerically more
advantageous by fitting the whole trajectories than the common approach by
averaging the trajectories over time [50, 11, 27].

####  4.2.1 Lorenz system

The most popular example of chaotic dynamics is the Lorenz63 system [named
after the author and the year it was proposed in 37] that represents a
simplified model of atmospheric convection for the chaotic behavior of the
weather. The governing equations of the Lorenz system are given by the
following ODE

$\begin{cases}\dot{x}&=\sigma(y-x),\\\ \dot{y}&=x(\rho-z)-y,\\\
\dot{z}&=xy-\beta z,\end{cases}$ (30)

where $x$, $y$, and $z$ denote variables proportional to convective intensity,
horizontal and vertical temperature differences and $u:=(\sigma,\rho,\beta)$
represents the model parameters known as Prandtl number ($\sigma$), Rayleigh
number ($\rho$), and an unnamed parameter ($\beta$) used for physical
proportions of the regions [43].

Figure 6: Lorenz inverse problem: marginal (diagonal) and pairwise (lower
triangle) sections of the joint density $p(u)$ by the time-averaged model
(left) and the STGP model (right) respectively.

The behavior of Lorenz63 system (30) strongly relies on the parameters. In
many studies, the parameter $\rho$ varies in $(0,\infty)$ and the other
parameters $\sigma$ and $\beta$ are held constant. In particular, (30) has a
stable equilibrium point at the origin for $\rho\in(0,1)$. For
$\rho\in(1,\gamma)$ with $\gamma=\sigma\frac{\sigma+\beta+3}{\sigma-\beta-1}$,
(30) has three equilibrium points, one unstable equilibrium point at the
origin and two stable equilibrium points at
${(\sqrt{\beta(\rho-1)},\sqrt{\beta(\rho-1)},\rho-1)}^{\mathsf{T}}$ and
${(-\sqrt{\beta(\rho-1)},-\sqrt{\beta(\rho-1)},\rho-1)}^{\mathsf{T}}$. When
$\rho>\gamma$, the equilibrium points become unstable and it results in
unstable spiral shaped trajectories. One classical configuration for the
parameter in (30) is $\sigma=10$, $\beta=\frac{8}{3}$, $\rho=28$ when the
system exhibits two-lobe orbits, also known as the butterfly effect [57] (See
the left panel of Figure 5). In this example, we seek to infer such parameter
$u^{\dagger}=(\sigma^{\dagger},\beta^{\dagger},\rho^{\dagger})=(10,8/3,28)$
based on the observed chaotic trajectories demonstrated in the middle panel of
Figure 5. Note the solutions $(x(t),y(t),z(t))$ highly depend on the initial
conditions $(x(0),y(0),z(0))$, we hence fix $(x(0),y(0),z(0))$ in the
following.

Due to the chaotic nature of the states $\\{(x(t),y(t),z(t)):t\in[0,\tau]\\}$,
we can treat these coordinates as random variables. In right panel of Figure
5, we demonstrate their marginal and pairwise distributions (diagonal and
lower triangle) estimated by a collection of states (upper triangle) along a
long-time trajectory solved with $u^{\dagger}$. For a given parameter
$u=(\sigma,\beta,\rho)$, we have the trajectory ${\mathcal{G}}(u)$ as the
following map:

${\mathcal{G}}(u):\;{\mathbb{R}}_{+}\rightarrow{\mathbb{R}}^{3},\quad
t\mapsto(x(t;u),y(t;u),z(t;u))$ (31)

where $(x(t;u),y(t;u),z(t;u))$ is the solution of (30) for given parameter
$u$. We generate spatiotemporal data from the chaotic dynamics (30) with
$u^{\dagger}=(\sigma^{\dagger},\beta^{\dagger},\rho^{\dagger})$ by observing
its trajectory on $J=100$ equally spaced time points
$t_{j}\in[t_{0},t_{0}+T]$: ${\bf X}(u^{\dagger})_{3\times
100}:=\\{{\mathcal{G}}(u^{\dagger})(t_{j})=(x(t_{j};u^{\dagger}),y(t_{j};u^{\dagger}),z(t_{j};u^{\dagger}))\\}_{j=1}^{J}$.
These observations can be viewed as a 3-dimensional time series that estimate
the empirical covariance $\Gamma_{\textrm{obs}}$ as in [11]. The inverse
problem involves learning the parameter $u$ given these observations, also
known as parameter identification [41].

Figure 7: Lorenz inverse problem: comparing posterior estimates of parameter
$u$ for two models (time-average and STGP) in terms of relative error of mean
$\mathrm{REM}=\frac{\|\hat{u}-u^{\dagger}\|}{\|u^{\dagger}\|}$. The upper row
shows the results by varying the spin-up $t_{0}$ and fixing $T=10$. The lower
row shows the results by varying the observation window size $T$ and fixing
$t_{0}=100$. Each experiment is repeated for 10 runs of EnK (EKI and EKS
respectively) with $J=500$ ensembles and the shaded regions indicate standard
deviations of such repeated experiments.

Following [11], we endow a log-Normal prior on $u$: $\log
u\sim{\mathcal{N}}(\mu_{0},\sigma_{0}^{2})$ with $\mu_{0}=(2.0,1.2,3.3)$ and
$\sigma_{0}=(0.2,0.5,0.15)$. We compare the two likelihood models (12) and
(16) for this dynamical inverse problem. For the time-averaged model (12),
instead of the 3-dimensional time series from the trajectory (31), we
substitute $\overline{{\bf X}}(u)_{3\times 1}$ with $\overline{{\bf
X}^{\star}}(u)_{9\times 1}={\mathcal{O}}{\mathcal{G}}^{\star}(u)(t)$ by
averaging the following augmented trajectory ${\mathcal{G}}^{\star}(u)(t)$ in
time [11]:

${\mathcal{G}}^{\star}(u)(t)=(x(t),y(t),z(t),x^{2}(t),y^{2}(t),z^{2}(t),x(t)y(t),x(z)z(t),y(t)z(t))$

For the spatiotemporal likelihood model STGP (16), we set the correlation
length $\ell_{\bf x}=0.4$ and $\ell_{t}=0.1$ for the spatial kernel
${\mathcal{C}}_{\bf x}$ and the temporal kernel ${\mathcal{C}}_{t}$
respectively. They are chosen to reflect the spatial and temporal resolutions.

We first notice that the spatiotemporal modeling facilitates the learning of
the true parameter $u^{\dagger}$. As illustrated in Figure 6, despite of the
rough landscape, the marginal (e.g.
$p(\sigma,\beta^{\dagger},\rho^{\dagger})$) and pairwise (e.g.
$p(\sigma,\beta,\rho^{\dagger})$) sections of the joint density $p(u)$ by the
STGP model (16) are more convex in the neighbourhood of $u^{\dagger}$ compared
with the time-averaged model (12). This verifies the implication of Theorem
3.2 on their difference in convexity. Therefore, particle based algorithms
such as EnK methods have higher chance of concentrating their ensemble
particles around the true parameter value $u^{\dagger}$, leading to better
estimates. Here, the roughness of the posterior creates barrier for the direct
application of MCMC algorithms. Therefore, we apply more robust EnK methods
for parameter estimation.

We run each of the EnK algorithms for $N=50$ iterations and choose the
ensembles (of size $J$) when its ensemble mean attains the minimal error in
estimating the parameter $u$ with reference to its true value $u^{\dagger}$.
In practice, EnK algorithms usually converge quickly within a few iterations
so $N=50$ suffices the need for most applications.

To investigate the roles of spin-up length $t_{0}$ and observation window size
$T$, we run EnK multiple times while varying each of the two quantities one at
a time. Seen from Figure 7, we observe consistently smaller errors by the STGP
model (16) compared with the time-averaged model (12). More specifically, the
upper row indicates that the estimation errors, measured by
$\mathrm{REM}=\frac{\|\hat{u}-u^{\dagger}\|}{\|u^{\dagger}\|}$, are not every
sensitive to the spin-up $t_{0}$ given sufficient window size $T=10$. On the
other hand, for fixed spin-up $t_{0}=100$, both models decrease errors with
increasing window size $T$ as they aggregate more information. However, the
STGP model requires only about $\frac{1}{4}$ time length as the time-averaged
model to attain accuracy at the same level ($T=1$ vs $T=4$). This supports
that the STGP is preferable to the time-average approach as the former may add
a small overhead for the statistical inference but could save much more in
resolving the physics (solving ODE/PDE), which is usually more expensive.

Figure 8: Lorenz inverse problem: marginal (diagonal) and pairwise (lower triangle) distributions estimated with 10000 samples (upper triangle) by the pCN algorithm based on NN emulators for the time-averaged model (left) and the STGP model (right) respectively. Red dots (lower triangle) are selective 10000 ensemble particles from running the EKS algorithm. Model-Algorithms | J=50 | J=100 | J=200 | J=500 | J=1000
---|---|---|---|---|---
Tavg-EKI | 0.06 (0.03) | 0.09 (0.03) | 0.09 (0.01) | 0.06 (0.04) | 0.07 (0.02)
Tavg-EKS | 0.10 (0.02) | 0.07 (4.62e-03) | 0.05 (2.60e-03) | 0.03 (3.04e-03) | 0.03 (8.56e-04)
STGP-EKI | 0.07 (0.03) | 0.04 (0.03) | 0.03 (0.02) | 0.02 (0.03) | 0.02 (0.01)
STGP-EKS | 0.09 (0.03) | 0.05 (0.03) | 0.03 (0.02) | 3.97e-04 (1.06e-03) | 5.52e-04 (6.37e-04)
Table 2: Lorenz inverse problem: comparing posterior estimates of parameter
$u$ for two models, time-average (Tavg) and STGP, in terms of relative error
of median $\mathrm{REM}=\frac{\|\hat{u}-u^{\dagger}\|}{\|u^{\dagger}\|}$. Each
experiment is repeated for 10 runs of EnK (EKI and EKS respectively) and the
numbers in the bracket are standard deviations of such repeated experiments.

Now we set spin-up $t_{0}=100$ long enough to ignore effect of initial
condition in the dynamics and choose the observation window size $T=10$. We
compare the two models (12) (16) using EnK algorithms with different ensemble
sizes ($J$) to obtain an estimate $\hat{u}$ of the parameter $u$. Figure 19
shows that the STGP model performs better than the time-averaged model in
generating smaller errors (REM) for almost all cases. In general more
ensembles help reduce the errors except for the time-averaged model using EKI
algorithm. Note, the STGP model with EKS algorithm yields parameter estimates
with the lowest errors. Table 2 summarizes the REM’s by different combinations
of the two likelihood models (time-averaged and STGP) and two EnK algorithms
(EKI and EKS). Again we can see consistent advantage of the spatiotemporal
likelihood model STGP (16) over the simple time-averaged model (12) in
producing more accurate parameter estimation.

Figure 9: Lorenz inverse problem: comparing forward predictions
$\overline{{\mathcal{G}}}({\bf x},t_{*})$ based on the time-averaged model and
the STGP model.

Next, we apply CES (Section 2) [11, 34] to quantify the uncertainty of the
estimate $\hat{u}$. Direct application of MCMC suffers from the extremely low
acceptance rate because of the rough density landscape (Figure 6). Ensemble
particles from EnK algorithm cannot provide rigorous systematic UQ due to the
ensemble collapse [48, 49, 16, 9] (See red dots in Figure 8). Therefore, we
run approximate MCMC based on NN emulators built from EnK outputs
$\\{u_{n}^{(j)},{\mathcal{G}}(u_{n}^{(j)})\\}_{j=1,n=0}^{J,N}$. Note, we have
different structures for the observed data in the two models (12) (16):
9-dimensional summary of time series for the time-averaged model (12) and
$3\times 100$ time series for the STGP model (16). Therefore we build densely
connected NN (DNN) ${\mathcal{G}}^{e}:{\mathbb{R}}^{3}\to{\mathbb{R}}^{9}$ for
the former and DNN-RNN (recurrent NN) type of network
${\mathcal{G}}^{e}:{\mathbb{R}}^{3}\to{\mathbb{R}}^{3\times 100}$ for the
latter to account for their different data structures in the forward output.
Figure 8 compares the marginal (diagonal) and pairwise (lower triangle)
posterior densities of $u$ estimated by 10000 samples (upper triangle) of the
pCN algorithm based on the corresponding NN emulators for the two models. The
spatiotemporal model STGP (16) yields more reasonable UQ results compared with
the time-averaged model (12).

Finally, we consider the forward prediction $\overline{{\mathcal{G}}}({\bf
x},t_{*})$ (26) for $t_{*}\in[t_{0},t_{0}+1.5T]$ with $J=500$ EKS ensembles
corresponding to the lowest error. Figure 9 compares the prediction results
given by these two models. The result by the STGP model is very close to the
truth till $t=113$ while the prediction by the time-averaged model quickly
departs from the truth only after $t=102$. The STGP model predicts the future
of this challenging chaotic dynamics significantly better than the time-
averaged model.

####  4.2.2 Rössler system

Figure 10: Rössler dynamics: single-lobe orbits (left), chaotic solutions
(middle) and coordinates’ distributions (right).

Next we consider the following Rössler system [25] governed by the system of
autonomous differential equations:

$\begin{cases}\dot{x}&=-y-z,\\\ \dot{y}&=x+ay,\\\
\dot{z}&=b+z(x-c).\end{cases}$ (32)

where $a,b,c>0$ are parameters determining the behavior of the system. The
Rössler attractor was originally discovered by German biochemist Otto Eberhard
Rössler [46, 47]. When $c^{2}>4ab$, the system (32) exhibits continuous-time
chaos and has two unstable equilibrium points
$(a\gamma_{-},-\gamma_{-},\gamma_{-})$ and
$(a\gamma_{+},-\gamma_{+},\gamma_{+})$ with
$\gamma_{+}=\frac{c+\sqrt{c^{2}-4ab}}{2a},\,\gamma_{-}=\frac{c-\sqrt{c^{2}-4ab}}{2a}$.
Note that the Rössler attractor has similarities to the Lorenz attractor,
nevertheless it has a single lobe and offers more flexibility in the
qualitative analysis. The true parameter that we try to infer is
$u^{\dagger}=(a^{\dagger},b^{\dagger},c^{\dagger})=(0.2,0.2,5.7)$. Figure 10
illustrates the single-lobe orbits (left), the chaotic solutions (middle) and
their marginal and pairwise distributions (right) of their coordinates viewed
as random variables.

Note, the Rössler dynamics evolve at a lower rate compared with the Lorenz63
dynamics (compare the middle panels of Figure 10 and Figure 5). Therefore, we
adopt a longer spin-up length ($t_{0}=1000$) and a larger window size
($T=100$). For the STGP model (16), spatiotemporal data are generated by
observing the trajectory (31) of the chaotic dynamics (32) with
$u^{\dagger}=(a^{\dagger},b^{\dagger},c^{\dagger})$ for $J=100$ time points in
$[t_{0},t_{0}+T]$. We also augment the time-averaged data with second-order
moments for the time-averaged model (12). In this Bayesian inverse problem, we
adopt a log-Normal prior on $u$: $\log
u\sim{\mathcal{N}}(\mu_{0},\sigma_{0}^{2})$ with $\mu_{0}=(-1.5,-1.5,2.0)$ and
$\sigma_{0}=(0.15,0.15,0.2)$. Once again, with spatiotemporal likelihood model
STGP (16), learning the true parameter value $u^{\dagger}$ becomes easier
because the posterior density $p(u)$ concentrates more on $u^{\dagger}$
compared with the time-averaged model (12), as indicated by Theorem 3.2. See
Figure 20 for the comparison on their marginal and pairwise sections of the
joint density $p(u)$.

Figure 11: Rössler inverse problem: comparing posterior estimates of parameter
$u$ for two models (time-average and STGP) in terms of relative error of mean
$\mathrm{REM}=\frac{\|\hat{u}-u^{\dagger}\|}{\|u^{\dagger}\|}$. The upper row
shows the results by varying the spin-up $t_{0}$ and fixing $T=100$. The lower
row shows the results by varying the window size $T$ and fixing $t_{0}=100$.
Each experiment is repeated for 10 runs of EnK (EKI and EKS respectively) with
$J=500$ ensembles and the shaded regions indicate standard deviations of such
repeated experiments.

We also compare the two models (12) (16) when investigating the roles of spin-
up length $t_{0}$ and observation window size $T$ in Figure 11. Despite of the
consistent smaller errors (expressed in terms of REM) by the STGP model, REM
is not every sensitive to the spin-up $t_{0}$ given sufficient window size
$T=100$. However, for fixed spin-up $t_{0}=100$, the STGP model (16) is
superior than the time-averaged approach (12) in reducing the estimation error
using smaller observation time window $T$: the former requires only half time
length as the latter to attain the same level of accuracy ($T=30$ vs $T=60$
with EKI and $T=20$ vs $T=40$ with EKS).

Figure 12: Rössler inverse problem: marginal (diagonal) and pairwise (lower triangle) distributions estimated with 10000 samples (upper triangle) by the pCN algorithm based on NN emulators for the time-averaged model (left) and the STGP model (right) respectively. Red dots (lower triangle) are selective 10000 ensemble particles from running the EKS algorithm. Model-Algo | J=50 | J=100 | J=200 | J=500 | J=1000
---|---|---|---|---|---
Tavg-EKI | 0.16 (0.09) | 0.11 (0.06) | 0.10 (0.07) | 0.07 (0.04) | 0.11 (0.07)
Tavg-EKS | 0.06 (0.02) | 0.06 (7.61e-03) | 0.06 (6.20e-03) | 0.06 (5.37e-03) | 0.06 (2.53e-03)
STGP-EKI | 0.02 (0.02) | 0.01 (0.01) | 0.02 (0.02) | 0.01 (9.09e-03) | 0.01 (0.02)
STGP-EKS | 0.02 (0.01) | 2.47e-03 (0.02) | 7.63e-04 (2.86e-03) | 4.23e-04 (2.45e-04) | 3.62e-04 (1.19e-04)
Table 3: Rössler inverse problem: comparing posterior estimates of parameter
$u$ for two models (time-average and STGP) in terms of relative error of
median $\mathrm{REM}=\frac{\|\hat{u}-u^{\dagger}\|}{\|u^{\dagger}\|}$. Each
experiment is repeated for 10 runs of EnK (EKI and EKS respectively) and the
numbers in the bracket are standard deviations of such repeated experiments.

Now we fix $t_{0}=1000$ and $T=100$. Figure 21 compares these two models (12)
(16) in terms of REM’s of the parameter estimation by EnK algorithms with
different ensemble sizes ($J$). The STGP model (16) shows universal advantage
over the time-averaged model (12) in generating smaller REM’s. Note, the time-
averaged model becomes over-fitting if running EKS more than 10 iterations, a
phenomenon also reported in [29, 28]. Table 3 summarizes the REM’s by
different combinations of likelihood models and EnK algorithms and confirms
the consistent advantage of the STGP model over the time-averaged model in
rendering more accurate parameter estimation.

We apply CES (Section 2) [11, 34] for the UQ. Based on the EKS ($J=500$)
outputs, we build DNN ${\mathcal{G}}^{e}:{\mathbb{R}}^{3}\to{\mathbb{R}}^{9}$
for the time-averaged model (12) and DNN-RNN
${\mathcal{G}}^{e}:{\mathbb{R}}^{3}\to{\mathbb{R}}^{3\times 100}$ for the STGP
model (16) to account for their different data structures. Figure 12 compares
the marginal and pairwise posterior densities of $u$ estimated by 10000
samples of the pCN algorithm based on the corresponding NN emulators for the
two models. The STGP model (16) generates more appropriate UQ results than the
time-averaged model (12) does. Finally, we consider the forward prediction
$\overline{{\mathcal{G}}}({\bf x},t_{*})$ (26) for
$t_{*}\in[t_{0},t_{0}+1.5T]$ with $J=500$ EKS ensembles corresponding to the
lowest error. Figure 13 shows that the STGP model provides better prediction
consistent with the truth throughout the whole time window while the result by
the time-averaged model deviates from the truth quickly after $t=1020$.

Figure 13: Rössler inverse problem: comparing forward predictions
$\overline{{\mathcal{G}}}({\bf x},t_{*})$ based on the time-averaged model and
the STGP model.

####  4.2.3 Chen system

Figure 14: Chen dynamics: double-scroll attractor (left), chaotic solutions
(middle) and coordinates’ distributions (right).

Yet another chaotic dynamical system we consider is the Chen system [10]
described by the following ODE:

$\begin{cases}\dot{x}&=a(y-x),\\\ \dot{y}&=(c-a)x-xz+cy,\\\ \dot{z}&=xy-
bz.\end{cases}$ (33)

where $a,b,c>0$ are parameters. When $a=35,b=3,c=28$, the system (33) has a
double-scroll chaotic attractor often observed from a physical, electronic
chaotic circuit. The true parameter that we will infer is
$u^{\dagger}=(a^{\dagger},b^{\dagger},c^{\dagger})=(35,3,28)$. With
$u^{\dagger}$, the system has three unstable equilibrium states given by
$(0,0,0)$, $(\gamma,\gamma,2c-a)$, and $(-\gamma,-\gamma,2c-a)$ where
$\gamma=\sqrt{b(2c-a)}$ [60]. Figure 14 illustrates the two-scroll attractor
(left), the chaotic trajectories (middle) and their marginal and pairwise
distributions (right) of their coordinates viewed as random variables.

The Chen dynamics has trajectories changing rapidly as the Lorenz63 dynamics
(compare the middle panels of Figure 14 and Figure 5). Therefore we adopt the
same spin-up length ($t_{0}=100$) and observation window size ($T=10$) as in
the Lorenz inverse problem (Section 4.2.1). We generate the spatiotemporal
data and the augmented time-averaged summary data by observing the trajectory
of (33) over $[t_{0},t_{0}+T]$ solved with $u^{\dagger}$ similarly as in the
previous sections. A log-Nomral prior is adopted for $u$: $\log
u\sim{\mathcal{N}}(\mu_{0},\sigma_{0}^{2})$ with $\mu_{0}=(3.5,1.2,3.3)$ and
$\sigma_{0}=(0.35,0.5,0.15)$. The STGP model (16) still posses more convex
posterior density $p(u)$ than the time-averaged model (12) as illustrated by
its marginal and pairwise sections plotted in Figure 22.

Figure 15: Chen inverse problem: comparing posterior estimates of parameter
$u$ for two models (time-average and STGP) in terms of relative error of mean
$\mathrm{REM}=\frac{\|\hat{u}-u^{\dagger}\|}{\|u^{\dagger}\|}$. The upper row
shows the results by varying the spin-up $t_{0}$ and fixing $T=10$. The lower
row shows the results by varying the observation window size $T$ and fixing
$t_{0}=100$. Each experiment is repeated for 10 runs of EnK (EKI and EKS
respectively) with $J=500$ ensembles and the shaded regions indicate standard
deviations of such repeated experiments.

Varying the spin-up length $t_{0}$ and the observation window size $T$ one at
a time in Figure 15, we observe similar advantage of the STGP model compared
with the time-averaged model regardless of the insensitivity of errors with
respect to $t_{0}$. Similarly, the STGP model demands a smaller observation
window than the time-averaged model ($T=2$ vs $T=6$ with EKI and $T=2$ vs
$T=3$ with EKS) to reach the same level of accuracy.

Model-Algo | J=50 | J=100 | J=200 | J=500 | J=1000
---|---|---|---|---|---
Tavg-EKI | 0.07 (0.03) | 0.04 (0.04) | 0.04 (0.04) | 0.05 (0.04) | 0.04 (0.04)
Tavg-EKS | 0.12 (0.03) | 0.10 (0.02) | 0.09 (0.02) | 0.09 (0.01) | 0.09 (0.01)
STGP-EKI | 0.14 (0.09) | 0.09 (0.08) | 0.09 (0.08) | 0.03 (0.03) | 0.01 (9.87e-03)
STGP-EKS | 0.07 (0.04) | 0.05 (0.04) | 0.01 (0.01) | 2.89e-03 (6.07e-03) | 3.32e-04 (4.66e-04)
Table 4: Chen inverse problem: comparing posterior estimates of parameter $u$
for two models (time-average and STGP) in terms of relative error of median
$\mathrm{REM}=\frac{\|\hat{u}-u^{\dagger}\|}{\|u^{\dagger}\|}$. Each
experiment is repeated for 10 runs of EnK (EKI and EKS respectively) and the
numbers in the bracket are standard deviations of such repeated experiments.

Figure 16: Chen inverse problem: marginal (diagonal) and pairwise (lower
triangle) distributions estimated with 10000 samples (upper triangle) by the
pCN algorithm based on NN emulators for the time-averaged model (left) and the
STGP model (right) respectively. Red dots (lower triangle) are selective 10000
ensemble particles from running the EKS algorithm.

Again we see the merit of the STGP model (16) in reducing the error (REM) of
parameter estimation compared with the time-averaged model (12) in various
combinations of EnK algorithms with different ensemble sizes ($J$) in Figure
23 and Table 4. As in the previous problem (Section 4.2.2), similar over-
fitting (bottom left of Figure 23) by the time-averaged model occurs if
running EKS algorithms more than 5 iterations (or earlier).

UQ results (Figure 16) by CES show the STGP model estimates the uncertainty of
parameter $u$ more appropriately than the time-averaged model. Finally, though
the prediction is challenging to the Chen dynamics (33), the STGP model still
performs much better than the time-averaged model by predicting more accurate
trajectory for longer time ($t=111$ vs $t=101$) as shown in Figure 17.

Figure 17: Chen inverse problem: comparing forward predictions
$\overline{{\mathcal{G}}}({\bf x},t_{*})$ based on the time-averaged model and
the STGP model.

##  5 Conclusion

In this paper, we investigate the inverse problems with spatiotemporal data.
We compare the Bayesian models based on STGP with traditional static and time-
averaged models that do not fully integrate the spatiotemporal information. By
fitting the trajectories of the observed data, the STGP model provides more
effective parameter estimation and more appropriate UQ. We explain the
superiority of the STGP model in theorems showing that it renders more convex
likelihood that facilitates the parameter learning. We demonstrate the
advantage of the spatiotemporal modeling using an inverse problem constrained
by an advection-diffusion PDE and three inverse problems involving chaotic
dynamics.

Theorems 3.1 and 3.2 compare the STGP model with the static and the time-
averaged models regarding their statistical convexity. These novel qualitative
results imply that the parameter learning (based on EnK methods) with the STGP
model converges faster than the other two traditional methods. In the future
work, we will explore a quantitative characterization on their convergence
rates particularly in terms of covariance properties.

The STGP model (16) considered in this paper has a classical separation
structure in their joint kernel. This may not be sufficient to characterize
complex spatiotemporal relationships, e.g. the temporal evoluation of spatial
dependence (TESD) [33]. We will expand this work by considering non-stationary
non-separable STGP models [14, 61, 55] to account for more complicated space-
time interactions in these spatiotemporal inverse problems.

## Acknowledgement

SL is supported by NSF grant DMS-2134256.

## Appendix

##  Appendix A Proofs

######  Theorem (3.1).

If we set the maximal eigenvalues of ${\bf C}_{\bf x}$ and ${\bf C}_{t}$ such
that $\lambda_{\max}({\bf C}_{\bf x})\lambda_{\max}({\bf
C}_{t})\leq\sigma^{2}_{\varepsilon}$, then the following inequality holds
regarding the Fisher information matrices, ${\mathcal{I}}_{\textrm{\tiny S}}$
and ${\mathcal{I}}_{\textrm{\tiny ST}}$, of the static model and the STGP
model respectively:

${\mathcal{I}}_{\textrm{\tiny ST}}(u)\geq{\mathcal{I}}_{\textrm{\tiny S}}(u)$
(34)

If we control the maximal eigenvalues of ${\bf C}_{\bf x}$ and ${\bf C}_{t}$
such that $\lambda_{\max}({\bf C}_{\bf x})\lambda_{\max}({\bf C}_{t})\leq
J\lambda_{\min}(\Gamma_{\textrm{obs}})$, then the following inequality holds
regarding the Fisher information matrices, ${\mathcal{I}}_{\textrm{\tiny T}}$
and ${\mathcal{I}}_{\textrm{\tiny ST}}$, of the time-averaged model and the
STGP model respectively:

${\mathcal{I}}_{\textrm{\tiny ST}}(u)\geq{\mathcal{I}}_{\textrm{\tiny T}}(u)$
(35)

###### Proof.

Denote ${\bf Y}_{0}={\bf Y}-{\bf M}$. We have
$\Phi_{*}(u)=\frac{1}{2}\mathrm{tr}\left[{\bf V}_{*}^{-1}{{\bf
Y}}^{\mathsf{T}}_{0}{\bf U}_{*}^{-1}{\bf Y}_{0}\right]$ with $*$ being S or
ST. ${\bf U}_{\textrm{\tiny S}}$, ${\bf V}_{\textrm{\tiny S}}$, ${\bf
U}_{\textrm{\tiny ST}}$ and ${\bf V}_{\textrm{\tiny ST}}$ are specified in
(17). We notice that both ${\bf U}_{*}$ and ${\bf V}_{*}$ are symmetric, then
we have

$\displaystyle\frac{\partial\Phi_{*}}{\partial u_{i}}$
$\displaystyle=\frac{1}{2}\left\\{\mathrm{tr}\left[{\bf
V}_{*}^{-1}\frac{\partial{{\bf Y}}^{\mathsf{T}}_{0}}{\partial u_{i}}{\bf
U}_{*}^{-1}{\bf Y}_{0}\right]+\mathrm{tr}\left[{\bf V}_{*}^{-1}{{\bf
Y}}^{\mathsf{T}}_{0}{\bf U}_{*}^{-1}\frac{\partial{\bf Y}_{0}}{\partial
u_{i}}\right]\right\\}=\mathrm{tr}\left[{\bf V}_{*}^{-1}{{\bf
Y}}^{\mathsf{T}}_{0}{\bf U}_{*}^{-1}\frac{\partial{\bf Y}_{0}}{\partial
u_{i}}\right]$ (36) $\displaystyle\frac{\partial^{2}\Phi_{*}}{\partial
u_{i}\partial u_{j}}$ $\displaystyle=\mathrm{tr}\left[{\bf V}_{*}^{-1}{{\bf
Y}}^{\mathsf{T}}_{0}{\bf U}_{*}^{-1}\frac{\partial^{2}{\bf Y}_{0}}{\partial
u_{i}\partial u_{j}}\right]+\mathrm{tr}\left[{\bf
V}_{*}^{-1}\frac{\partial{{\bf Y}}^{\mathsf{T}}_{0}}{\partial u_{i}}{\bf
U}_{*}^{-1}\frac{\partial{\bf Y}_{0}}{\partial u_{j}}\right]$

Due to the i.i.d. assumption in both models, ${\bf Y}_{0}$ is independent of
either $\frac{\partial{\bf Y}_{0}}{\partial u_{i}}$ or $\frac{\partial^{2}{\bf
Y}_{0}}{\partial u_{i}\partial u_{j}}$. Therefore

$\displaystyle({\mathcal{I}}_{*})_{ij}$
$\displaystyle=\mathrm{E}\left[\frac{\partial^{2}\Phi_{*}}{\partial
u_{i}\partial u_{j}}\right]=\mathrm{E}\left[\mathrm{tr}\left({\bf
V}_{*}^{-1}\frac{\partial{{\bf Y}}^{\mathsf{T}}_{0}}{\partial u_{i}}{\bf
U}_{*}^{-1}\frac{\partial{\bf Y}_{0}}{\partial u_{j}}\right)\right]$ (37)
$\displaystyle=\mathrm{E}\left[\mathrm{vec}{\left(\frac{\partial{\bf
Y}_{0}}{\partial u_{i}}\right)}^{\mathsf{T}}({\bf V}_{*}^{-1}\otimes{\bf
U}_{*}^{-1})\mathrm{vec}\left(\frac{\partial{\bf Y}_{0}}{\partial
u_{j}}\right)\right]$

For any ${\bf w}=(w_{1},\cdots,w_{p})\in{\mathbb{R}}^{p}$ and ${\bf w}\neq{\bf
0}$, denote $\tilde{\bf
w}:=\sum_{i=1}^{p}w_{i}\mathrm{vec}\left(\frac{\partial{\bf Y}_{0}}{\partial
u_{i}}\right)$. To prove ${\mathcal{I}}_{\textrm{\tiny
ST}}(u)\geq{\mathcal{I}}_{\textrm{\tiny S}}(u)$, it suffices to show
${\tilde{\bf w}}^{\mathsf{T}}({\bf V}_{\textrm{\tiny ST}}\otimes{\bf
U}_{\textrm{\tiny ST}})^{-1}\tilde{\bf w}\geq{\tilde{\bf w}}^{\mathsf{T}}({\bf
V}_{\textrm{\tiny S}}\otimes{\bf U}_{\textrm{\tiny S}})^{-1}\tilde{\bf w}$.

By [Theorem 4.2.12 in 26], we know that any eigenvalue of ${\bf
V}_{*}\otimes{\bf U}_{*}$ has the format as a product of eigenvalues of ${\bf
V}_{*}$ and ${\bf U}_{*}$ respectively, i.e. $\lambda_{k}({\bf
V}_{*}\otimes{\bf U}_{*})=\lambda_{i}({\bf V}_{*})\lambda_{j}({\bf U}_{*})$,
where where $\\{\lambda_{j}(M)\\}$ are the ordered eigenvalues of $M$, i.e.
$\lambda_{1}(M)\geq\cdots\geq\lambda_{d}(M)$. By the given condition we have

$\lambda_{IJ}(({\bf V}_{\textrm{\tiny ST}}\otimes{\bf U}_{\textrm{\tiny
ST}})^{-1})=\lambda_{1}^{-1}({\bf V}_{\textrm{\tiny ST}}\otimes{\bf
U}_{\textrm{\tiny ST}})=\lambda_{1}^{-1}({\bf C}_{t})\lambda_{1}^{-1}({\bf
C}_{\bf x})\geq\sigma^{-2}_{\varepsilon}=\lambda_{1}(({\bf V}_{\textrm{\tiny
S}}\otimes{\bf U}_{\textrm{\tiny S}})^{-1})$ (38)

Thus it completes the proof of the first inequality.

Similarly by the second condition, we have

$\lambda_{IJ}(({\bf V}_{\textrm{\tiny ST}}\otimes{\bf U}_{\textrm{\tiny
ST}})^{-1})=\lambda_{1}^{-1}({\bf C}_{t})\lambda_{1}^{-1}({\bf C}_{\bf x})\geq
J^{-1}\lambda_{\min}^{-1}(\Gamma_{\textrm{obs}})=\lambda_{1}({\bf
V}_{\textrm{\tiny S}}^{-}\otimes{\bf U}_{\textrm{\tiny S}}^{-1})$ (39)

and complete the proof of the second inequality. ∎

######  Theorem (3.2).

If we choose ${\bf C}_{\bf x}=\Gamma_{\textrm{obs}}$ and require the maximal
eigenvalue of ${\bf C}_{t}$, $\lambda_{\max}({\bf C}_{t})\leq J$, then the
following inequality holds regarding the Fisher information matrices,
${\mathcal{I}}_{\textrm{\tiny T}}$ and ${\mathcal{I}}_{\textrm{\tiny ST}}$, of
the time-averaged model and the STGP model respectively:

${\mathcal{I}}_{\textrm{\tiny ST}}(u)\geq{\mathcal{I}}_{\textrm{\tiny T}}(u)$
(40)

###### Proof.

Denote ${\bf Y}_{0}={\bf Y}-{\bf M}$. We have
$\Phi_{*}(u)=\frac{1}{2}\mathrm{tr}\left[{\bf V}_{*}^{-1}{{\bf
Y}}^{\mathsf{T}}_{0}{\bf U}_{*}^{-1}{\bf Y}_{0}\right]$ with $*$ being T or
ST. ${\bf U}_{\textrm{\tiny T}}$, ${\bf V}_{\textrm{\tiny T}}$, ${\bf
U}_{\textrm{\tiny ST}}$ and ${\bf V}_{\textrm{\tiny ST}}$ are specified in
(17).

By the similar argument of the proof in Theorem 3.1, we have

$({\mathcal{I}}_{*})_{ij}=\mathrm{E}\left[\frac{\partial^{2}\Phi_{*}}{\partial
u_{i}\partial u_{j}}\right]=\mathrm{tr}\left[{\bf
V}_{*}^{-1}\mathrm{E}\left(\frac{\partial{{\bf Y}}^{\mathsf{T}}_{0}}{\partial
u_{i}}{\bf U}_{*}^{-1}\frac{\partial{\bf Y}_{0}}{\partial
u_{j}}\right)\right]$ (41)

For any ${\bf w}=(w_{1},\cdots,w_{p})\in{\mathbb{R}}^{p}$ and ${\bf w}\neq{\bf
0}$, denote ${\bf W}:=\sum_{i,j=1}^{p}w_{i}\mathrm{E}\left(\frac{\partial{{\bf
Y}}^{\mathsf{T}}_{0}}{\partial u_{i}}{\bf U}_{*}^{-1}\frac{\partial{\bf
Y}_{0}}{\partial u_{j}}\right)w_{j}$. We know ${\bf W}\geq{\bf 0}_{J\times
J}$. It suffices to show $\mathrm{tr}[{\bf V}_{\textrm{\tiny ST}}^{-1}{\bf
W}]\geq\mathrm{tr}[{\bf V}_{\textrm{\tiny T}}^{-1}{\bf W}]$.

By the corollary [38] of Von Neumann’s trace inequality [39], we have

$\sum_{j=1}^{J}\lambda_{j}({\bf V}_{*}^{-1})\lambda_{J-j+1}({\bf
W})\leq\mathrm{tr}({\bf V}_{*}^{-1}{\bf W})\leq\sum_{j=1}^{J}\lambda_{j}({\bf
V}_{*}^{-1})\lambda_{j}({\bf W})$ (42)

where $\\{\lambda_{j}(M)\\}$ are the ordered eigenvalues of $M$, i.e.
$\lambda_{1}(M)\geq\cdots\geq\lambda_{d}(M)$. The only non-zero eigenvalue of
${\bf V}_{\tiny T}^{-}=J^{-2}(\bm{1}_{J}{\bm{1}}^{\mathsf{T}}_{J})$ is
$\lambda_{1}({\bf V}_{\tiny T}^{-})=J^{-1}$. Therefore, we have

$\mathrm{tr}[{\bf V}_{\tiny T}^{-}{\bf W}]\leq J^{-1}\lambda_{1}({\bf
W})\leq\lambda_{J}({\bf V}_{\textrm{\tiny ST}}^{-1})\lambda_{1}({\bf
W})+\sum_{j=1}^{J-1}\lambda_{j}({\bf V}_{\textrm{\tiny
ST}}^{-1})\lambda_{J-j+1}({\bf W})\leq\mathrm{tr}[{\bf V}_{\textrm{\tiny
ST}}^{-1}{\bf W}]$ (43)

where $\lambda_{J}({\bf V}_{\textrm{\tiny ST}}^{-1})=\lambda_{1}^{-1}({\bf
C}_{t})\geq J^{-1}$ and $\lambda_{j}({\bf V}_{\textrm{\tiny
ST}}^{-1}),\lambda_{j}({\bf W})\geq 0$. ∎

##  Appendix B More Numerical Results

Figure 18: Advection-diffusion inverse problem: auto-correlations of
observations in space (left) and time (right) respectively. Figure 19: Lorenz
inverse problem: comparing posterior estimates of parameter $u$ for two models
(time-average and STGP) in terms of relative error of median
$\mathrm{REM}=\frac{\|\hat{u}-u^{\dagger}\|}{\|u^{\dagger}\|}$. Each
experiment is repeated for 10 runs of EnK (EKI and EKS respectively) and
shaded regions indicate $5\sim 95\%$ quantiles of such repeated results.

Figure 20: Rössler inverse problem: marginal (diagonal) and pairwise (lower
triangle) sections of the joint density $p(u)$ by the time-averaged model
(left) and the STGP model (right) respectively. Figure 21: Rössler inverse
problem: comparing posterior estimates of parameter $u$ for two models (time-
average and STGP) in terms of relative error of median
$\mathrm{REM}=\frac{\|\hat{u}-u^{\dagger}\|}{\|u^{\dagger}\|}$. Each
experiment is repeated for 10 runs of EnK (EKI and EKS respectively) and
shaded regions indicate $5\sim 95\%$ quantiles of such repeated results.

Figure 22: Chen inverse problem: marginal (diagonal) and pairwise (lower
triangle) sections of the joint density $p(u)$ by the time-averaged model
(left) and the STGP model (right) respectively. Figure 23: Chen inverse
problem: comparing posterior estimates of parameter $u$ for two models (time-
average and STGP) in terms of relative error of median
$\mathrm{REM}=\frac{\|\hat{u}-u^{\dagger}\|}{\|u^{\dagger}\|}$. Each
experiment is repeated for 10 runs of EnK (EKI and EKS respectively) and
shaded regions indicate $5\sim 95\%$ quantiles of such repeated results.

## References

  * [1] Henry Abarbanel.  Predicting the Future: completing models of observed complex systems, volume 1.  Springer New York, 2013.
  * [2] HN Agiza and MT Yassen.  Synchronization of Rossler and Chen chaotic dynamical systems using active control.  Physics Letters A, 278(4):191–197, 2001.
  * [3] A. Beskos, F. J. Pinski, J. M. Sanz-Serna, and A. M. Stuart.  Hybrid Monte-Carlo on Hilbert spaces.  Stochastic Processes and their Applications, 121:2201–2230, 2011\.
  * [4] Alexandros Beskos.  A stable manifold MCMC method for high dimensions.  Statistics & Probability Letters, 90:46–52, 2014.
  * [5] Alexandros Beskos, Mark Girolami, Shiwei Lan, Patrick E. Farrell, and Andrew M. Stuart.  Geometric MCMC for infinite-dimensional inverse problems.  Journal of Computational Physics, 335, 2017.
  * [6] Alexandros Beskos, Gareth Roberts, Andrew Stuart, and Jochen Voss.  MCMC methods for diffusion bridges.  Stochastics and Dynamics, 8(03):319–350, 2008.
  * [7] Robert Bishop.  Chaos.  In Edward N. Zalta, editor, The Stanford Encyclopedia of Philosophy. Metaphysics Research Lab, Stanford University, Spring 2017 edition, 2017.
  * [8] Chris Brooks.  Chaos in foreign exchange markets: a sceptical view.  Computational Economics, 11(3):265–281, 1998.
  * [9] Neil K. Chada, Andrew M. Stuart, and Xin T. Tong.  Tikhonov regularization within ensemble kalman inversion, 2019.
  * [10] Guanrong Chen and Tetsushi Ueta.  Yet another chaotic attractor.  International Journal of Bifurcation and chaos, 9(07):1465–1466, 1999.
  * [11] Emmet Cleary, Alfredo Garbuno-Inigo, Shiwei Lan, Tapio Schneider, and Andrew M Stuart.  Calibrate, emulate, sample, 2020.
  * [12] Maxime Conjard and Henning Omre.  Spatio-temporal inversion using the selection kalman model.  Frontiers in Applied Mathematics and Statistics, 7, apr 2021.
  * [13] Simon L Cotter, Gareth O Roberts, AM Stuart, and David White.  MCMC methods for functions: modifying old algorithms to make them faster.  Statistical Science, 28(3):424–446, 2013.
  * [14] N. Cressie and C.K. Wikle.  Statistics for Spatio-Temporal Data.  CourseSmart Series. Wiley, 2011.
  * [15] Masoumeh Dashti and Andrew M. Stuart.  The Bayesian Approach to Inverse Problems, pages 311–428.  Springer International Publishing, Cham, 2017.
  * [16] J. de Wiljes, S. Reich, and W. Stannat.  Long-time stability and accuracy of the ensemble kalman–bucy filter for fully observed processes and small measurement noise.  SIAM Journal on Applied Dynamical Systems, 17(2):1152–1181, 2018\.
  * [17] David Echeverría Ciaurri and T. Mukerji.  A robust scheme for spatio-temporal inverse modeling of oil reservoirs.  01 2009.
  * [18] S Effah-Poku, William Obeng-Denteh, and IK Dontwi.  A study of chaos in dynamical systems.  Journal of Mathematics, 2018, 2018.
  * [19] Geir Evensen.  Sequential data assimilation with a nonlinear quasi-geostrophic model using monte carlo methods to forecast error statistics.  Journal of Geophysical Research, 99(C5):10143, 1994.
  * [20] Geir Evensen and Peter Jan van Leeuwen.  Assimilation of geosat altimeter data for the agulhas current using the ensemble kalman filter with a quasigeostrophic model.  Monthly Weather Review, 124(1):85–96, jan 1996.
  * [21] R. A. Fisher and Edward John Russell.  On the mathematical foundations of theoretical statistics.  Philosophical Transactions of the Royal Society of London. Series A, Containing Papers of a Mathematical or Physical Character, 222(594-604):309–368, 1922.
  * [22] Alfredo Garbuno-Inigo, Franca Hoffmann, Wuchen Li, and Andrew M. Stuart.  Interacting langevin diffusions: Gradient structure and ensemble kalman sampler.  SIAM Journal on Applied Dynamical Systems, 19(1):412–441, 2020\.
  * [23] Alfredo Garbuno-Inigo, Nikolas Nüsken, and Sebastian Reich.  Affine invariant interacting langevin dynamics for bayesian inference.  SIAM Journal on Applied Dynamical Systems, 19(3):1633–1658, jan 2020.
  * [24] A.K. Gupta and D.K. Nagar.  Matrix Variate Distributions, chapter Chapter 2: MATRIX VARIATE NORMAL DISTRIBUTION.  Chapman and Hall/CRC, may 2018.
  * [25] A Hegazi, HN Agiza, and MM El-Dessoky.  Controlling chaotic behaviour for spin generator and Rossler dynamical systems with feedback control.  Chaos, Solitons & Fractals, 12(4):631–658, 2001.
  * [26] Roger A. Horn and Charles R. Johnson.  Topics in Matrix Analysis.  Cambridge University Press, apr 1991.
  * [27] Daniel Zhengyu Huang, Jiaoyang Huang, Sebastian Reich, and Andrew M. Stuart.  Efficient derivative-free bayesian inference for large-scale inverse problems, 2022.
  * [28] Marco A Iglesias.  A regularizing iterative ensemble kalman method for PDE-constrained inverse problems.  Inverse Problems, 32(2):025002, Jan 2016.
  * [29] Marco A Iglesias, Kody J H Law, and Andrew M Stuart.  Ensemble kalman methods for inverse problems.  Inverse Problems, 29(4):045001, Mar 2013.
  * [30] Vladimir G. Ivancevic and Tijana T. Ivancevic.  Complex Nonlinearity.  Springer Berlin Heidelberg, 2008.
  * [31] Charles E. Baukal Jr., Vladimir Gershtein, and Xianming Jimmy Li, editors.  Computational Fluid Dynamics in Industrial Combustion.  CRC Press, oct 2000.
  * [32] Shiwei Lan.  Adaptive dimension reduction to accelerate infinite-dimensional geometric markov chain monte carlo.  Journal of Computational Physics, 392:71 – 95, September 2019.
  * [33] Shiwei Lan.  Learning temporal evolution of spatial dependence with generalized spatiotemporal gaussian process models.  arXiv:1901.04030, 08 2021.
  * [34] Shiwei Lan, Shuyi Li, and Babak Shahbaba.  Scaling up bayesian uncertainty quantification for inverse problems using deep neural networks.  SIAM/ASA Journal on Uncertainty Quantification, to appear, 2022\.
  * [35] Eduardo Liz and Alfonso Ruiz-Herrera.  Chaos in discrete structured population models.  SIAM Journal on Applied Dynamical Systems, 11(4):1200–1214, jan 2012.
  * [36] Christopher J. Long, Patrick L. Purdon, Simona Temereanca, Neil U. Desai, Matti S. Hämäläinen, and Emery N. Brown.  State-space solutions to the dynamic magnetoencephalography inverse problem using high performance computing.  The Annals of Applied Statistics, 5(2B), jun 2011.
  * [37] Edward N. Lorenz.  Deterministic nonperiodic flow.  Journal of the Atmospheric Sciences, 20(2):130–141, mar 1963.
  * [38] Albert W. Marshall, Ingram Olkin, and Barry C. Arnold.  Inequalities: Theory of Majorization and Its Applications.  Springer New York, 2nd edition, 2011.
  * [39] L. Mirsky.  A trace inequality of john von neumann.  Monatshefte für Mathematik, 79:303–306, 1975.
  * [40] M. Morzfeld, J. Adams, S. Lunderman, and R. Orozco.  Feature-based data assimilation in geophysics.  Nonlinear Processes in Geophysics, 25(2):355–374, 2018.
  * [41] Elisa Negrini, Giovanna Citti, and Luca Capogna.  System identification through lipschitz regularized deep neural networks.  Journal of Computational Physics, 444:110549, nov 2021.
  * [42] +Alejandro Ojeda, +Marius Klug, +Kenneth Kreutz-Delgado, +Klaus Gramann, and +Jyoti Mishra.  A bayesian framework for unifying data cleaning, source separation and imaging of electroencephalographic signals.  bioRxiv, 2019.
  * [43] Edward Ott.  Strange attractors and chaotic motions of dynamical systems.  Reviews of Modern Physics, 53(4):655, 1981.
  * [44] Mirjeta Pasha, Arvind K. Saibaba, Silvia Gazzola, Malena I. Espanol, and Eric de Sturler.  Efficient edge-preserving methods for dynamic inverse problems, 2021.
  * [45] N. Petra and G. Stadler.  Model variational inverse problems governed by partial differential equations.  Technical report, The Institute for Computational Engineering and Sciences, The University of Texas at Austin., 2011.
  * [46] O.E. Rössler.  An equation for continuous chaos.  Physics Letters A, 57(5):397–398, 1976.
  * [47] O.E. Rossler.  An equation for hyperchaos.  Physics Letters A, 71(2):155–157, 1979.
  * [48] C. Schillings and A. Stuart.  Analysis of the ensemble kalman filter for inverse problems.  SIAM Journal on Numerical Analysis, 55(3):1264–1290, 2017.
  * [49] C. Schillings and A. M. Stuart.  Convergence analysis of ensemble kalman inversion: the linear, noisy case.  Applicable Analysis, 97(1):107–123, Oct 2017.
  * [50] Tapio Schneider, Shiwei Lan, Andrew Stuart, and João Teixeira.  Earth system modeling 2.0: A blueprint for models that learn from observations and targeted high-resolution simulations.  Geophysical Research Letters, 44(24):12,396–12,417, 2017.
  * [51] A I Shcherbakova, Y A Kupriyanova, and G V Zhikhareva.  Spatio-temporal analysis the results of solving the inverse problem of electrocardiography.  Journal of Physics: Conference Series, 2091(1):012028, nov 2021\.
  * [52] Pridi Siregar and Jean-Paul Sinteff.  Introducing spatio-temporal reasoning into the inverse problem in electroencephalography.  Artificial Intelligence in Medicine, 8(2):97–122, 1996.
  * [53] Andrew M Stuart.  Inverse problems: a Bayesian perspective.  Acta Numerica, 19:451–559, 2010.
  * [54] Umberto Villa, Noemi Petra, and Omar Ghattas.  hIPPYlib: An extensible software framework for large-scale inverse problems governed by PDEs; part i: Deterministic Inversion and Linearized Bayesian Inference, 2020.
  * [55] Kangrui Wang, Oliver Hamelijnck, Theodoros Damoulas, and Mark Steel.  Non-separable non-stationary random fields.  In Hal Daumé III and Aarti Singh, editors, Proceedings of the 37th International Conference on Machine Learning, volume 119 of Proceedings of Machine Learning Research, pages 9887–9897. PMLR, 13–18 Jul 2020\.
  * [56] M.W. Woolrich, M. Jenkinson, J.M. Brady, and S.M. Smith.  Fully bayesian spatio-temporal modeling of FMRI data.  IEEE Transactions on Medical Imaging, 23(2):213–231, feb 2004\.
  * [57] Shyi-Kae Yang, Chieh-Li Chen, and Her-Terng Yau.  Control of chaos in Lorenz system.  Chaos, Solitons & Fractals, 13(4):767–780, 2002.
  * [58] Ying Yang.  Source-space analyses in meg/eeg and applications to explore spatio-temporal neural dynamics in human vision.  2017\.
  * [59] Bing Yao and Hui Yang.  Physics-driven spatiotemporal regularization for high-dimensional predictive modeling: A novel approach to solve the inverse ECG problem.  Scientific Reports, 6(1), dec 2016.
  * [60] MT Yassen.  Chaos control of Chen chaotic dynamical system.  Chaos, Solitons & Fractals, 15(2):271–283, 2003.
  * [61] Bohai Zhang and Noel Cressie.  Bayesian inference of spatio-temporal changes of arctic sea ice.  Bayesian Analysis, 15(2):605–631, jun 2020.
  * [62] Yiheng Zhang, Alireza Ghodrati, and Dana H Brooks.  An analytical comparison of three spatio-temporal regularization methods for dynamic linear inverse problems in a common statistical framework.  Inverse Problems, 21(1):357–382, jan 2005.