# Spatiotemporal Besov Priors for Bayesian Inverse Problems

Shiwei Lan<EMAIL_ADDRESS>School of Mathematical & Statistical Sciences, Arizona
State University, 901 S Palm Walk, Tempe, AZ 85287  Mirjeta Pasha  Laboratory
for Decision and Information System, Massachusetts Institute of Technology, 77
Massachusetts Avenue, Cambridge, MA 02139  Shuyi Li  School of Mathematical &
Statistical Sciences, Arizona State University, 901 S Palm Walk, Tempe, AZ
85287  Weining Shen  Department of Statistics, University of California -
Irvine, Donald Bren Hall, Irvine, CA 92697

###### Abstract

Fast development in science and technology has driven the need for proper
statistical tools to capture special data features such as abrupt changes or
sharp contrast. Many inverse problems in data science require spatiotemporal
solutions derived from a sequence of time-dependent objects with these spatial
features, e.g. dynamic reconstruction of computerized tomography (CT) images
with edges. Conventional methods based on Gaussian processes (GP) often fall
short in providing satisfactory solutions since they tend to offer over-smooth
priors. Recently, the Besov process (BP), defined by wavelet expansions with
random coefficients, has emerged as a more suitable prior for Bayesian inverse
problems of this nature. While BP excels in handling spatial inhomogeneity, it
does not automatically incorporate temporal correlation inherited in the
dynamically changing objects. In this paper, we generalize BP to a novel
spatiotemporal Besov process (STBP) by replacing the random coefficients in
the series expansion with stochastic time functions as Q-exponential process
(Q-EP) which governs the temporal correlation structure. We thoroughly
investigate the mathematical and statistical properties of STBP. A white-noise
representation of STBP is also proposed to facilitate the inference.
Simulations, two limited-angle CT reconstruction examples and a highly non-
linear inverse problem involving Navier-Stokes equation are used to
demonstrate the advantage of the proposed STBP in preserving spatial features
while accounting for temporal changes compared with the classic STGP and a
time-uncorrelated approach.

Keywords: Spatiotemporal functional data analysis, Inhomogeneous data, $L_{q}$
regularization, Q-Exponential process, Edge-preserving priors

##  1 Introduction

Many modern science and engineering applications are presented as inverse
problems whose main goal is to recover parameters of interest from observed
data. These data may possess inhomogeneity in the sense that certain portions
differs from others significantly. For example, some medical images exhibit
sharp edges where properties undergo dramatic changes. Furthermore, these
complex datasets may be spatiotemporal, with solutions extending across both
space and time. One of the key challenges in these types of inverse problems
is to effectively capture the distinctive characteristics of high-dimensional
(potentially infinite-dimensional) objects using limited data. Surging need
has been posed for statistical methodology to appropriately impose
regularization or fill in prior information for these ill-posed inverse
problems in order to construct meaningful solutions [17].

In the literature of optimization, deterministic regularization methods for
inverse problems date back to 1943 with the seminal work by Tikhonov [69],
followed by Ivanov [32]. A regularization term reflecting properties of the
target solution is typically added to a pre-determined energy function that
depends on the forward model and the statistical assumptions of the
observational noise. More recent methods and algorithms developed in this
class have been applied in signal and image processing [36], geophysics and
seismic monitoring [66], tomography reconstruction [56], etc. In statistics,
there is also a long history of developing models with penalty based on the
prior knowledge of unknown parameters including but not limited to Bayesian
Lasso [55, 12] and bridge [59], spike-and-slab [52, 49, 31, 3], and horseshoe
[11, 57] priors. If we view the Bayesian regression model as an inverse
problem, there is a natural analogy between the likelihood (prior) and the
energy (regularization) function.

In nonparametric modeling, GP [61] has been widely used as an $L_{2}$ penalty
(negative log-density as a quadratic form) or a prior on the function space.
However, despite the flexibility, random functions generated from GPs often
exhibit excessive smoothing, which is not ideal for modeling inhomogeneous
objects such as images with sharp edges (See Figure 10 for illustration). To
address this issue, researchers have proposed a class of $L_{1}$ penalty based
priors including Laplace random field [58, 39, 2] and Besov process (BP) [42,
14, 33, 15]. There are also many heavy-tailed priors based on Markov random
fields such as Cauchy [50, 68] and total variation (TV) [45, 78], and data-
informed priors based on level set functions [16, 30] proposed for handling
inhomogeneity. These approaches have found significant applications in spatial
modeling [58], signal processing [39], imaging analysis [56], and inverse
problems [14].

Most of these sparsity-promoting and edge-preserving priors in the literature
work well in characterizing spatial inhomogeneity, but few of them are
tailored to specially address spatiotemporal targets and their temporal
correlations. In this paper, we focus on BP proposed by [42] for imaging
analysis and generalize it to the spatiotemporal domain. Historically, [43]
discovered that the TV prior degenerates to a Gaussian prior as the
discretization mesh becomes denser and loses the edge-preserving properties in
high-dimensional applications. Therefore, [42] proposed the BP prior defined
using wavelet basis and random coefficients following a (univariate)
$q$-exponential distribution and proved its discretization-invariant property.
Recently, [47] introduced a stochastic process based on a consistent
multivariate generalization of the $q$-exponential distribution, hence named
$Q$-exponential process (Q-EP). The Q-EP can be viewed as an explicit
probabilistic definition of BP with direct control on the correlation
structure and tractable prediction formula. In this paper, we propose a novel
_spatiotemporal Besov process (STBP)_ by replacing the (univariate)
$q$-exponential random coefficients in the series definition of BP with
stochastic time functions as Q-EP. The proposed STBP offers a flexible prior
in modeling functional data with spatial features while controlling the
temporal correlations explicitly through a covariance kernel. Similarly as BP,
STBP also includes spatiotemporal GP (STGP) as a special case for $q=2$ (See
Figure 1 for their relationship). Since our motivation is to model
inhomogeneous data with priors imposing sharp $L_{q}$ regularization, we focus
on $1\leq q\leq 2$ in this paper (See Figure 10 for the regularization effect
of parameter $q$). To the best of our knowledge, this is by far the first
spatiotemporal generalization of BP.

To justify STBP as a working prior for spatiotemporal inverse problems, Bayes
theorem in this setting is re-examined based on [14, 15]. With proper
assumptions on the likelihood, the posterior contraction theorems for Bayesian
inverse models with STBP priors are established based on [21, 72, 1].
Posterior contraction properties with Gaussian priors have been extensively
studied by [72, 73, 74, 22] for regression and classification, density
estimation, white noise models, and Bayesian linear inverse problems [37] and
nonlinear inverse problems [76]. [23] also studied the posterior contraction
of density estimation for a class of $L^{r}$-metrics ($1\leq r\leq\infty$)
based priors including GP, wavelet series and normal mixture. [1] studied the
posterior contraction theorems for a re-branded BP named $p$-exponential
process (which only differs from BP by a constant in the univariate
$q$-exponential distribution) for density estimation and white noise model.
Other works on the posterior contraction of Besov-type priors include [62, 20,
24]. Compared with the existing literature, our theoretic results are novel in
terms of: 1) generalization to spatiotemporal models for inverse problems; 2)
simplified contraction rates given in a class of Besov-type spaces contained
in $L^{q}$ spaces.

To facilitate Bayesian inference for models with STBP priors, we introduce a
novel white-noise representation [13] of the random function drawn from STBP
and take advantage of the dimension-independent MCMC algorithms [8] for
efficient implementation. The numerical advantages of the proposed STBP over
STGP have been supported by multiple experiments in dynamic CT reconstruction
and inverse problems constrained by highly nonlinear differential equation.
Our proposed work on STBP has multiple contributions to the literature of
spatiotemporal inverse problems:

  1. 1.

It generalizes BP to the spatiotemporal domain to simultaneously model the
spatial inhomogeneity and the temporal correlations.

  2. 2.

It provides theoretic characterization on the posterior contraction in the
infinite data limit, justifying its validity as a nonparametric learning tool.

  3. 3.

It demonstrates utility in dynamic CT reconstruction and indicates potential
impact on medical imaging analysis.

The rest of the paper is organized as follows. Section 2 provides a background
review on the Bayesian inverse problems and BP used as a flexible edge-
preserving prior. Section 3 introduces Q-EP as random coefficient functions on
the time domain. We then formally define STBP and study its theoretic
properties in Section 4. In Section 5 we describe a white-noise representation
of STBP that facilitates the inference for models with STBP prior. In Section
6 we demonstrate the advantage of the proposed STBP prior in retaining spatial
features and capturing temporal correlations for the spatiotemporal inverse
problems using a simulated regression, two dynamic CT reconstruction examples
and a nonlinear inverse problem involving Navier-Stokes equation. Finally we
conclude with some discussion on future research in Section 7.

##  2 Background on Besov priors for Bayesian Inverse Problems

The Bayesian approach to inverse problems [15] has recently gained increasing
popularity because it provides a natural framework for model calibration and
uncertainty quantification (UQ) as well. In this section, we review some
background about the set up of Bayesian inverse problems and the definition of
BP as a flexible prior for modeling objects with spatial features.

###  2.1 Bayesian Inverse Problems

We consider the inverse problem of recovering an unknown parameter
$u\in{\mathbb{X}}$ from a noisy observation $y\in{\mathbb{Y}}$ based on the
following Bayesian model

$\displaystyle y$
$\displaystyle={\mathcal{G}}(u)+\eta,\quad\eta\sim{\mathbb{Q}}_{0},$ (1)
$\displaystyle u$ $\displaystyle\sim\Pi,$

where both ${\mathbb{X}}$ and ${\mathbb{Y}}$ are separable Banach spaces,
${\mathcal{G}}:{\mathbb{X}}\to{\mathbb{Y}}$ is a forward mapping from the
parameter space ${\mathbb{X}}$ to the data space ${\mathbb{Y}}$, and
$\eta\in{\mathbb{Y}}$ denotes the random noise whose distribution
${\mathbb{Q}}_{0}$ is independent from the prior $\Pi$. We assume the
conditional $y|u$ is distributed according to the measure
${\mathbb{Q}}_{u}\ll{\mathbb{Q}}_{0}$ for $u$ $\Pi$-almost surely (a.s.), and
hence define the potential (negative log-likelihood) function
$\Phi:{\mathbb{X}}\times{\mathbb{Y}}\to{\mathbb{R}}$:

$\frac{d{\mathbb{Q}}_{u}}{d{\mathbb{Q}}_{0}}(y)=\exp(-\Phi(u;y)).$ (2)

The objective of Bayesian inverse problem is to seek the posterior solution of
$u|y$ whose distribution, denoted as $\Pi^{y}$, according to the Bayes’
theorem [67, 15], satisfies the requirement that if
$0<Z:=\int_{{\mathbb{X}}}\exp(-\Phi(u;y))\Pi(du)<+\infty$ for $y$
${\mathbb{Q}}_{0}$-a.s., then

$\frac{d\Pi^{y}}{d\Pi}(u)=\frac{1}{Z}\,\exp(-\Phi(u;y)).$ (3)

The forward operator ${\mathcal{G}}$ could be linear or nonlinear, possibly
encoding physical information represented by a system of ordinary or partial
differential equations (ODE/PDE). The resulted posterior $\Pi^{y}$ is usually
non-Gaussian with a complicated geometric structure even if a Gaussian prior
$\Pi=\mathcal{GP}(0,{\mathcal{C}})$ is adopted.

When there are sparse data but the targets are high-dimensional, the inverse
problems are ill-posed. Proper prior information is crucial to induce well-
defined solutions [17]. In a plethora of applications such as image deblurring
[63] and CT reconstruction [56], methods are developed with priors that are
sparse [9], and/or have sparse gradients [42, 15] to preserve the sharp edge
features. In this work we focus on Besov-type priors [42].

###  2.2 Besov Process

Let ${\mathcal{X}}\subset{\mathbb{R}}^{d}$ be the spatial domain, e.g.
$d$-dimensional torus, ${\mathcal{X}}={\mathbb{T}}^{d}=(0,1]^{d}$ for $d\leq
3$. Consider a separable Banach space $({\mathbb{X}},\|\cdot\|)$ with a
Schauder basis $\\{\phi_{\ell}\\}_{\ell=1}^{\infty}$, e.g. square integrable
function
$L^{2}({\mathcal{X}}):=\\{u:{\mathcal{X}}\to{\mathbb{R}}|\int_{\mathcal{X}}|u({\bf
x})|^{2}d{\bf x}<\infty\\}$ with Fourier basis. Any function
$u\in{\mathbb{X}}$ can then be represented by the following series:

$u({\bf x})=\sum_{\ell=1}^{\infty}u_{\ell}\phi_{\ell}({\bf x})$ (4)

Based on (4), we consider a norm $\|\cdot\|_{s,q}$ defined with a smoothness
parameter $s>0$ and an integrability parameter $q\geq 1$ [42, 14]:

$\|u(\cdot)\|_{s,q}=\left(\sum_{\ell=1}^{\infty}\ell^{\tau_{q}(s)q}|u_{\ell}|^{q}\right)^{\frac{1}{q}},\quad\tau_{q}(s)=\frac{s}{d}+\frac{1}{2}-\frac{1}{q}.$
(5)

We define the Banach space
$B^{s,q}({\mathcal{X}}):=\\{u:{\mathcal{X}}\to{\mathbb{R}}\,|\,\|u(\cdot)\|_{s,q}<\infty\\}$.
If $\\{\phi_{\ell}\\}_{\ell=1}^{\infty}$ is an $r$-regular wavelet basis for
$r>s$, then $B^{s,q}({\mathcal{X}})$ becomes the Besov space $B^{s}_{qq}$
[70]. In particular, if $q=2$ and $\\{\phi_{\ell}\\}_{\ell=1}^{\infty}$ form
the Fourier basis, then $B^{s,2}({\mathcal{X}})$ reduces to the Sobolev space
$H^{s}({\mathcal{X}})$ with the special case
$B^{0,2}({\mathcal{X}})=L^{2}({\mathcal{X}})$ assuming $s=0$.

For a given basis $\\{\phi_{\ell}\\}_{\ell=1}^{\infty}$, based on the the
series expansion (4), there is a one-to-one correspondence between the
$B^{s,q}({\mathcal{X}})$ function $u(\cdot)$ and the infinite sequence
$u:=\\{u_{\ell}\\}_{\ell=1}^{\infty}$ in a weighted $\ell^{q}$ space,
$\ell^{q,\tau}:=\\{u\in{\mathbb{R}}^{\infty}\,|\,\|u\|_{\tau,q}=\left(\sum_{\ell=1}^{\infty}\ell^{\tau
q}|u_{\ell}|^{q}\right)^{\frac{1}{q}}<\infty\\}$, which reduces to the regular
$\ell^{q}$ space when $\tau=0$. Hence
$\|u(\cdot)\|_{s,q}=\|u\|_{\tau_{q}(s),q}$. In the following, we will use $u$
to refer to both notations under such correspondence when there is no
confusion.

Now we define _Besov process (BP)_ $u(\cdot)$ in (4) by randomizing the
coefficients $\\{u_{\ell}\\}_{\ell=1}^{\infty}$. More specifically, we set for
$\ell\in{\mathbb{N}}$

$u_{\ell}:=\gamma_{\ell}\xi_{\ell},\quad\gamma_{\ell}=\kappa^{-\frac{1}{q}}\ell^{-\tau_{q}(s)},\quad\xi_{\ell}\overset{i.i.d.}{\sim}\pi_{\xi}(\cdot)\propto\exp\left(-\frac{1}{2}|\cdot|^{q}\right),$
(6)

where $\kappa>0$ is a scaling factor, and $\pi_{\xi}$ denotes the probability
density function of the _$q$ -exponential distribution_ [14, 15]. Though not
spelled out, such $q$-exponential distribution is actually a special case of
the following exponential power (EP, a.k.a generalized normal) distribution
$\mathrm{EP}(\mu,\sigma,q)$ with $\mu=0$, $\sigma=1$:

$p(\xi|\mu,\sigma,q)=\frac{q}{2^{1+1/q}\sigma\Gamma(1/q)}\exp\left\\{-\frac{1}{2}\left|\frac{\xi-\mu}{\sigma}\right|^{q}\right\\}.$
(7)

When $q=2$, this is just a normal distribution
${\mathcal{N}}(\mu,\sigma^{2})$. When $q=1$, it becomes a Laplace distribution
$L(\mu,b)$ with $\sigma=2^{-1/q}b$.

Denote infinite sequences $\gamma=\\{\gamma_{\ell}\\}_{\ell=1}^{\infty}$ and
$\xi=\\{\xi_{\ell}\\}_{\ell=1}^{\infty}$. Then $\xi$ is a random element of
the probability space $(\Omega,{\mathcal{B}}(\Omega),\mathbb{P})$ with
$\Omega={\mathbb{R}}^{\infty}$, product $\sigma$-algebra
${\mathcal{B}}(\Omega)$ and probability measure $\mathbb{P}$ defined by
extending the finite product of $\pi_{\xi}$ to the infinite product by the
Kolmogorov extension theorem [c.f. Theorem 29 in section A.2.1 of 15]. Then we
define _Besov measure_ as the pushforward of $\mathbb{P}$ as follows.

######  Definition 1 (Besov Measure).

Let $\mathbb{P}$ be the measure of random sequences $\xi\in\Omega$. Suppose we
have the following map

$f_{\gamma}:\Omega\to B^{s,q}({\mathcal{X}}),\quad\xi\mapsto
u=\sum_{\ell=1}^{\infty}u_{\ell}\phi_{\ell}=\sum_{\ell=1}^{\infty}\gamma_{\ell}\xi_{\ell}\phi_{\ell},$
(8)

where $\gamma_{\ell}$ and $\xi_{\ell}$ are defined in (6). Then the
pushforward $f_{\gamma}^{\sharp}\mathbb{P}$ is Besov measure on
$B^{s,q}({\mathcal{X}})$, denoted as
${\mathcal{B}}(\kappa,B^{s,q}({\mathcal{X}}))$, and we say BP $u$ follows the
Besov measure, i.e., $u\sim{\mathcal{B}}(\kappa,B^{s,q}({\mathcal{X}}))$.

######  Remark 1.

If $q=2$ and $\\{\phi_{\ell}\\}_{\ell=1}^{\infty}$ is either a wavelet or
Fourier basis, we obtain a Gaussian measure with the Cameron-Martin space
$B_{22}^{s}$ [70], which is the Hilbert space $H^{s}({\mathcal{X}})$. Indeed,
we have (8) reduced to a zero mean Gaussian random element based on Karhunen-
Lovève representation (refer to Figure 1):

$u({\bf
x})=\kappa^{-\frac{1}{2}}\sum_{\ell=1}^{\infty}\ell^{-\frac{s}{d}}\xi_{\ell}\phi_{\ell}({\bf
x}),\quad\xi_{\ell}\overset{iid}{\sim}{\mathcal{N}}(0,1).$

Define
$\|\xi\|_{q}:=\left(\sum_{\ell=1}^{\infty}|\xi_{\ell}|^{q}\right)^{\frac{1}{q}}$
for $\xi\in\Omega$. Then we have
$\|u\|_{s,q}=\kappa^{-\frac{1}{q}}\|\xi\|_{q}$. The following formal Lebesgue
density can be made rigorous by Fernique theorem [14] (Refer to Theorem 4.2):

$\mathbb{P}(d\xi)=p(\xi)d\xi,\quad
p(\xi)=\prod_{\ell=1}^{\infty}\pi_{\xi}(\xi_{\ell})\propto\exp\left\\{-\frac{1}{2}\|\xi\|_{q}^{q}\right\\}=\exp\left\\{-\frac{\kappa}{2}\|u\|_{s,q}^{q}\right\\}.$
(9)

When used in the optimization to obtain a parameter estimate, the logarithm of
Besov prior density (9) serves as an $L_{q}$ regularization term. Larger
regularization parameter ($q>0$) promotes smoother solutions, as demonstrated
in a simulated regression example in Figure 10. In practice, Besov prior is
often adopted for $q=1$ to preserve edges in imaging analysis.

Define the Banach space of $q$-integrable functions in
$(\Omega,{\mathcal{B}}(\Omega),\mathbb{P})$ as
$L^{q}_{\mathbb{P}}(\Omega;B^{s,q}({\mathcal{X}}))=\\{u:{\mathcal{X}}\times\Omega\to{\mathbb{R}}|\mathbb{E}(\|u\|_{s,q}^{q})<\infty\\}$.
We notice that for $u\sim{\mathcal{B}}(\kappa,B^{s,q}({\mathcal{X}}))$,
$u\notin L^{q}_{\mathbb{P}}(\Omega;B^{s,q}({\mathcal{X}}))$ because
$\mathbb{E}(\|u\|_{s,q}^{q})=\kappa^{-1}\mathbb{E}(\|\xi\|_{q}^{q})=\infty$
due to the iid assumption on $\xi$ in (6). However, the following theorem
states that Besov random draw as in (8) has limit in a proper $q$-integrable
function space [Thorem 4 of 15] and hence in a subset of
$L^{q}({\mathcal{X}})$ in probability (almost surely).

######  Theorem 2.1.

If $u\sim{\mathcal{B}}(\kappa,B^{s,q}({\mathcal{X}}))$ as in (8), then $u\in
L^{q}_{\mathbb{P}}(\Omega;B^{s^{\prime},q}({\mathcal{X}}))$ for all
$s^{\prime}<s-\frac{d}{q}$.

###### Proof.

See Appendix A.1. ∎

######  Remark 2.

This theorem implies that for any random draw from
${\mathcal{B}}(\kappa,B^{s,q}({\mathcal{X}}))$, we need to consider its
$q$-integrability in a larger ambient space $B^{s^{\prime},q}({\mathcal{X}})$
for some $s^{\prime}<s-\frac{d}{q}$ ($B^{s,q}({\mathcal{X}})\subset
B^{s^{\prime},q}({\mathcal{X}})$, refer to Proposition 4.1).

Next, we generalize the series representation of a Besov random function (8)
to a representation for STBP by replacing the random variable $\xi_{\ell}$
with a stochastic process $\xi_{\ell}(\cdot)$ on the temporal domain
${\mathcal{T}}\subset{\mathbb{R}}_{+}$. For this purpose, in the following we
will first introduce a properly defined process Q-EP for $\xi_{\ell}(\cdot)$
that generalizes the $q$-exponential random variable $\xi_{\ell}$ and has the
capability of capturing the temporal dependence in data.

##  3 $Q$-exponential Process Valued Random Coefficients

###  3.1 Multivariate Generalization of $Q$-exponential Distribution

To generalize aforementioned $q$-exponential (6) (or univariate EP (7)) random
variable, $\xi_{\ell}$, to a multivariate random vector, $\bm{\xi}_{\ell}$,
and further a stochastic process, $\xi_{\ell}(\cdot)$, we have two important
requirements by the Kolmogorov’ extension theorem [54]: i) exchangeability of
the joint distribution, i.e. $p(\bm{\xi}_{1:J})=p(\bm{\xi}_{\tau(1:J)})$ for
any finite permutation $\tau$; and ii) consistency of the marginalization,
i.e. $p(\bm{\xi}_{1})=\int p(\bm{\xi}_{1},\bm{\xi}_{2})d\bm{\xi}_{2}$.

Consider the process $\xi_{\ell}(\cdot)$ defined in a finite temporal domain
${\mathcal{T}}\subset{\mathbb{R}}_{+}$, i.e., there exists $T<\infty$ such
that ${\mathcal{T}}\subset[0,T]$. Suppose we observe $\xi_{\ell}(t)$ at $J$
time points, $t_{1},\cdots,t_{J}\in{\mathcal{T}}$, then we need to define the
distribution of
$\bm{\xi}_{\ell}=(\xi_{\ell}(t_{1}),\cdots,\xi_{\ell}(t_{J}))$. [47]
investigate the family of elliptic contour distributions [34, 19] and propose
the following consistent _multivariate $q$-exponential distribution_ for
$\bm{\xi}_{\ell}$.

######  Definition 2.

A multivariate $q$-exponential distribution, denoted as
$\mathrm{q}\\!-\\!\mathrm{ED}_{J}(\bm{\mu},{\bf C})$, has the following
density

$p(\bm{\xi}|\bm{\mu},{\bf C},q)=\frac{q}{2}(2\pi)^{-\frac{J}{2}}|{\bf
C}|^{-\frac{1}{2}}r^{(\frac{q}{2}-1)\frac{J}{2}}\exp\left\\{-\frac{r^{\frac{q}{2}}}{2}\right\\},\quad
r(\bm{\xi})={(\bm{\xi}-\bm{\mu})}^{\mathsf{T}}{\bf
C}^{-1}(\bm{\xi}-\bm{\mu}).$ (10)

######  Remark 3.

When $q=2$, $\mathrm{q}\\!-\\!\mathrm{ED}_{J}(\bm{\mu},{\bf C})$ reduces to
MVN ${\mathcal{N}}_{J}(\bm{\mu},{\bf C})$. When $J=1$, if we let $C=1$, then
we have the density for $\xi$ as
$p(\xi)\propto|\xi|^{\frac{q}{2}-1}\exp\left\\{-\frac{1}{2}|\xi|^{q}\right\\}$,
differing from the original setup (6) by a term $|\xi|^{\frac{q}{2}-1}$. This
is needed for the consistency of temporal process generalization.

Based on [Theorem 1 of 35], [47] prove the following theorem [Theorem 3.3 of
47] which states that the multivariate $q$-exponential random vector following
distribution (10) satisfies the conditions of Kolmogorov’s extension theorem
hence can be generalized to a stochastic process.

######  Theorem 3.1.

The multivariate $q$-exponential distribution (10) is both exchangeable and
consistent.

To generate random vectors
$\bm{\xi}\sim\mathrm{q}\\!-\\!\mathrm{ED}_{J}(\bm{\mu},{\bf C})$, one can take
advantage of the stochastic representation [64, 10, 34], as detailed in the
following proposition [c.f. Theorem 2.1 and Proposition A.1 47]. This will be
needed for the Bayesian inference in Section 5.

######  Proposition 3.1.

If $\bm{\xi}\sim\mathrm{q}\\!-\\!\mathrm{ED}_{J}(\bm{\mu},{\bf C})$, then we
have

$\bm{\xi}=\bm{\mu}+R{\bf L}S$ (11)

where $S\sim\mathrm{Unif}({\mathcal{S}}^{J+1})$ uniformly distributed on the
unit-sphere ${\mathcal{S}}^{J+1}$, ${\bf L}$ is the Cholesky factor of ${\bf
C}$ such that ${\bf C}={\bf L}{{\bf L}}^{\mathsf{T}}$, $R\perp S$ and
$R^{q}\overset{d}{=}r^{\frac{q}{2}}\sim\Gamma\left(\alpha=\frac{J}{2},\beta=\frac{1}{2}\right)=\chi^{2}(J)$.

###  3.2 $Q$-exponential Process

With a covariance (symmetric and positive-definite) kernel
${\mathcal{C}}:{\mathcal{T}}\times{\mathcal{T}}\to{\mathbb{R}}$, we define
_$q$ -exponential process (Q-EP)_ based on the multivariate $q$-exponential
distribution (10).

######  Definition 3.

A (centered) $q$-exponential process $\xi(t)$ with kernel ${\mathcal{C}}$,
$\mathrm{q}\\!-\\!\mathcal{EP}(0,{\mathcal{C}})$, is a collection of random
variables such that any finite set,
$\bm{\xi}:=(\xi(t_{1}),\cdots,\xi(t_{J}))$, follows a multivariate
$q$-exponential distribution $\mathrm{q}\\!-\\!\mathrm{ED}_{J}({\bf 0},{\bf
C})$, where ${\bf C}=[{\mathcal{C}}(t_{j},t_{j^{\prime}})]{J\times J}$.

######  Remark 4.

When $q=2$, $\mathrm{q}\\!-\\!\mathcal{EP}(0,{\mathcal{C}})$ becomes
$\mathcal{GP}(0,{\mathcal{C}})$. When $q\in[1,2)$,
$\mathrm{q}\\!-\\!\mathcal{EP}(0,{\mathcal{C}})$ lends flexibility to modeling
functional data with more regularization than GP. See Figures 1 and 10, and
more details in Section 6.

In the square-integrable function space over ${\mathcal{T}}$, i.e.
$L^{2}({\mathcal{T}})$, the covariance kernel ${\mathcal{C}}$ is associated
with a Hilbert-Schmidt (HS) integral operator
$T_{\mathcal{C}}:L^{2}({\mathcal{T}})\to
L^{2}({\mathcal{T}}),\xi(\cdot)\mapsto\int_{\mathcal{T}}{\mathcal{C}}(\cdot,t^{\prime})\xi(t^{\prime})\mu(dt^{\prime})$
which has eigen-pairs
$\\{\lambda_{\ell},\psi_{\ell}(\cdot)\\}_{\ell=1}^{\infty}$ such that for
$\forall\ell\in{\mathbb{N}}$,
$T_{\mathcal{C}}\psi_{\ell}(t)=\psi_{\ell}(t)\lambda_{\ell}$ and
$\|\psi_{\ell}\|_{2}=1$. Then $\\{\psi_{\ell}\\}_{\ell=1}^{\infty}$ serves a
basis of $L^{2}({\mathcal{T}})$. Denote
$\lambda:=\\{\lambda_{\ell}\\}_{\ell=1}^{\infty}$. We assume $T_{\mathcal{C}}$
is a trace-class operator, i.e.,
$\mathrm{tr}(T_{\mathcal{C}}):=\|\lambda\|_{1}=\sum_{\ell=1}^{\infty}\lambda_{\ell}<\infty$.
[Theorem 3.4 of 47] shows that we have a Karhunen-Loéve type of theorem on the
series representation of random function
$\xi(\cdot)\sim\mathrm{q}\\!-\\!\mathcal{EP}(0,{\mathcal{C}})$ drawn from
Q-EP.

######  Theorem 3.2 (Karhunen-Loéve).

If $\xi(\cdot)\sim\mathrm{q}\\!-\\!\mathcal{EP}(0,{\mathcal{C}})$ with a
trace-class HS operator $T_{\mathcal{C}}$ having eigen-pairs
$\\{\lambda_{\ell},\psi_{\ell}(\cdot)\\}_{\ell=1}^{\infty}$, then we have the
following series representation for $\xi(t)$:

$\xi(t)=\sum_{\ell=1}^{\infty}\xi_{\ell}\psi_{\ell}(t),\quad\xi_{\ell}:=\int_{\mathcal{T}}\xi(t)\psi_{\ell}(t)\mu(dt)\sim\mathrm{q}\\!-\\!\mathrm{ED}(0,\lambda_{\ell}),$
(12)

where $\mathrm{E}[\xi_{\ell}]=0$ and
$\mathrm{Cov}(\xi_{\ell},\xi_{\ell^{\prime}})=\lambda_{\ell}\delta_{\ell\ell^{\prime}}$
with Dirac function $\delta_{\ell\ell^{\prime}}=1$ if $\ell=\ell^{\prime}$ and
$0$ otherwise. Moreover, we have

$\mathrm{E}[\|\xi(\cdot)\|_{2}^{2}]=\sum_{\ell=1}^{\infty}\mathrm{E}[\xi_{\ell}^{2}]=\mathrm{tr}(T_{\mathcal{C}})<\infty.$

######  Remark 5.

By re-scaling $\xi_{\ell}$ in (12), we have the series representation of Q-EP
$\xi(\cdot)$ in the same format as BP in (6):

$\xi_{\ell}:=\gamma_{\ell}\xi_{\ell}^{\star},\quad\gamma_{\ell}=\sqrt{\lambda_{\ell}},\quad\xi_{\ell}^{\star}\overset{iid}{\sim}\mathrm{q}\\!-\\!\mathrm{ED}(0,1)\sim\pi_{\xi}(\cdot).$
(13)

If we choose $\sqrt{\lambda_{\ell}}=\ell^{-\tau_{q}(s)}$ as in (13), then
$\mathrm{q}\\!-\\!\mathcal{EP}(0,{\mathcal{C}})$ process becomes equivalent to
${\mathcal{B}}(1,B^{s,q}({\mathcal{T}}))$ process. From this perspective, we
can view Q-EP as a probabilistic definition of BP (See Figure 1).

The above Theorem 3.2 states that for the given basis
$\\{\psi_{\ell}\\}_{\ell=1}^{\infty}$ on the time domain ${\mathcal{T}}$, we
can identify any random draw
$\xi(\cdot)\sim\mathrm{q}\\!-\\!\mathcal{EP}(0,{\mathcal{C}})$ with the
associated infinite sequence $\xi=\\{\xi_{\ell}\\}_{\ell=1}^{\infty}$ as in
(12). Similarly as in Section 2.2, we can define
$\|\xi(\cdot)\|_{s,q}=\|\xi\|_{\tau_{q}(s),q}=\left(\sum_{\ell=1}^{\infty}\ell^{\tau_{q}(s)q}|\xi_{\ell}|^{q}\right)^{\frac{1}{q}}$
and
$B^{s,q}({\mathcal{T}})=\\{\xi:{\mathcal{T}}\to{\mathbb{R}}\,|\,\|\xi(\cdot)\|_{s,q}<\infty\\}$.

In the probability space $(\Omega,{\mathcal{B}}(\Omega),\mathbb{P})$ with
$\Omega={\mathbb{R}}$ and $\mathbb{P}$ defined by
$\mathrm{q}\\!-\\!\mathcal{EP}(0,{\mathcal{C}})$, we consider a Banach space
defined as
$L^{p}_{\mathbb{P}}(\Omega,L^{p}({\mathcal{T}}))=\\{\xi:{\mathcal{T}}\times\Omega\to{\mathbb{R}}\,|\,\mathbb{E}(\|\xi\|_{p}^{p})<\infty\\}$.
From Theorem 3.2, we immediately have that if
$\xi(\cdot)\sim\mathrm{q}\\!-\\!\mathcal{EP}(0,{\mathcal{C}})$ with a trace-
class HS operator $T_{\mathcal{C}}$, then $\xi(\cdot)\in
L^{2}_{\mathbb{P}}(\Omega,L^{2}({\mathcal{T}}))$. More general integrability
of $\xi(\cdot)$ relates to the summability of eigenvalues $\lambda$ of the HS
operator $T_{\mathcal{C}}$, as expressed in the following assumption.

######  Assumption 1.

Suppose $\lambda=\\{\lambda_{\ell}\\}_{\ell=1}^{\infty}$ are eigenvalues of HS
operator $T_{\mathcal{C}}$ for the kernel ${\mathcal{C}}$ in Definition 3.

  1. (i)

We assume for $s^{\prime}<s-\frac{d}{q}$,

$\sqrt{\lambda}\in\ell^{q,\tau_{q}(s^{\prime})},\quad
i.e.,\quad\|\sqrt{\lambda}\|_{\tau_{q}(s^{\prime}),q}^{q}=\sum_{\ell=1}^{\infty}\ell^{\tau_{q}(s^{\prime})q}\lambda_{\ell}^{\frac{q}{2}}<\infty.$

  2. (ii)

In particular when $\tau_{q}(s^{\prime})=0$, i.e.
$s^{\prime}=\frac{d}{2}-\frac{d}{q}$, we assume

$\lambda\in\ell^{\frac{q}{2}},\quad
i.e.,\quad\|\lambda\|_{\frac{q}{2}}^{\frac{q}{2}}=\sum_{\ell=1}^{\infty}\lambda_{\ell}^{\frac{q}{2}}<\infty.$

######  Theorem 3.3.

If $\xi(\cdot)\sim\mathrm{q}\\!-\\!\mathcal{EP}(0,{\mathcal{C}})$ with a
trace-class HS operator $T_{\mathcal{C}}$ satisfying Assumption 1-(i), then
$\xi(\cdot)\in L^{q}_{\mathbb{P}}(\Omega,B^{s^{\prime},q}({\mathcal{T}}))$. If
Assumption 1-(ii) holds instead, then $\xi(\cdot)\in
L^{q}_{\mathbb{P}}(\Omega,L^{q}({\mathcal{T}}))$ and in particular
$\mathbb{E}[\|\xi(\cdot)\|_{q}^{q}]=\|\lambda\|_{\frac{q}{2}}^{\frac{q}{2}}<\infty$.

###### Proof.

See Appendix A.1 ∎

Under Assumption 1-(i), e.g. $\sqrt{\lambda_{\ell}}=\ell^{-\tau_{q}(s)}$,
$\mathrm{q}\\!-\\!\mathcal{EP}(0,{\mathcal{C}})$ process becomes equivalent to
the BP, ${\mathcal{B}}(1,B^{s,q}({\mathcal{T}}))$, which only differs from the
_$p$ -exponential process_ [1] by a constant in the definition
$\pi_{\xi}(\cdot)\propto\exp\left(-\frac{1}{p}|\cdot|^{p}\right)$. See Figure
1 for more illustration of their relationship. Therefore, the posterior
concentration theory developed by [1] for the $p$-exponential process applies
to the $\mathrm{q}\\!-\\!\mathcal{EP}(0,{\mathcal{C}})$ process. For
$\xi(\cdot)\in B^{s^{\prime},q}({\mathcal{T}})$ with
$s^{\prime}<s-\frac{d}{q}$, we define the concentration function of the
$q$-exponential measure $\mu$ (hence
${\mathcal{B}}(1,B^{s,q}({\mathcal{T}}))$) at $\xi=\xi^{\dagger}$ as follows.

$\varphi_{\xi^{\dagger}}(\varepsilon)=\inf_{h\in
B^{s,q}({\mathcal{T}}):\|h-\xi^{\dagger}\|_{s^{\prime},q}\leq\varepsilon}\frac{1}{2}\|h\|_{s,q}^{q}-\log\mu(\|\xi\|_{s^{\prime},q}\leq\varepsilon).$
(14)

Quoted from [Theorem 3.1 and Lemma 5.14 of 1], the following general
contraction theorem for $p$-exponential process (also BP, and hence Q-EP)
priors will be used in the proof of posterior contraction Theorem 4.7 for STBP
priors.

######  Theorem 3.4.

Let $\mu$ be a $\mathrm{q}\\!-\\!\mathcal{EP}(0,{\mathcal{C}})$ measure
satisfying Assumption 1-(i) in the separable Banach space
$(B^{s^{\prime},q}({\mathcal{T}}),\|\cdot\|_{s^{\prime},q})$, where
$q\in[1,2]$. Let $\xi\sim\mu$ and the true parameter $\xi^{\dagger}\in
B^{s^{\prime},q}({\mathcal{T}})$. Assume $\varepsilon_{n}>0$ such that
$\varphi_{\xi^{\dagger}}(\varepsilon_{n})\leq n\varepsilon_{n}^{2}$, where
$n\varepsilon_{n}^{2}\gtrsim 1$. Then for any $C>1$, there exists a measurable
set $B_{n}\subset B^{s^{\prime},q}({\mathcal{T}})$ and a constant $R>0$
depending on $C$ and $q$, such that

$\displaystyle\log N(4\varepsilon_{n},B_{n},\|\cdot\|_{s^{\prime},q})$
$\displaystyle\leq$ $\displaystyle Rn\varepsilon_{n}^{2},$ (15)
$\displaystyle\mu(\xi\notin B_{n})$ $\displaystyle\leq$
$\displaystyle\exp(-Cn\varepsilon_{n}^{2}),$ (16)
$\displaystyle\mu(\|\xi-\xi^{\dagger}\|_{s^{\prime},q}<2\varepsilon_{n})$
$\displaystyle\geq$ $\displaystyle\exp(-n\varepsilon_{n}^{2}).$ (17)

where $N(4\varepsilon_{n},B_{n},\|\cdot\|_{s^{\prime},q})$ is the minimal
number of $\|\cdot\|_{s^{\prime},q}$-balls of radius $4\varepsilon_{n}$ to
cover $B_{n}$.

##  4 Spatiotemporal Besov Process

Now we generalize the Banach space $B^{s,q}({\mathcal{X}})$ to include the
temporal domain. Let the coefficients $\\{u_{\ell}\\}_{\ell=1}^{\infty}$ in
(4) be $L^{p}({\mathcal{T}})$ functions over some bounded temporal domain
${\mathcal{T}}\subset{\mathbb{R}}_{+}$. Denote
${\mathcal{Z}}={\mathcal{X}}\times{\mathcal{T}}$ and ${\bf z}=({\bf x},t)$.
Then we obtain a spatiotemporal function $u({\bf z})=u({\bf x},t)$ on
${\mathcal{Z}}$ by the following series expansion with an infinite sequence of
$L^{p}({\mathcal{T}})$ functions:

$u({\bf z})=\sum_{\ell=1}^{\infty}u_{\ell}(t)\phi_{\ell}({\bf x}),\quad
u_{\ell}(\cdot)\in L^{p}({\mathcal{T}}),\;\forall\ell\in{\mathbb{N}}.$ (18)

Denote the infinite sequence
$u_{\mathcal{T}}:=\\{u_{\ell}(\cdot)\\}_{\ell=1}^{\infty}$. We define the
following norm $\|\cdot\|_{\tau,q,p}$ for $u_{\mathcal{T}}$ with a spatial
(BP) index $q\geq 1$ and a temporal (Q-EP) index $p\geq 1$:

$\|u_{\mathcal{T}}\|_{\tau,q,p}=\left(\sum_{\ell=1}^{\infty}\ell^{\tau
q}\|u_{\ell}(\cdot)\|^{q}_{p}\right)^{\frac{1}{q}}.$ (19)

Denote the space of such infinite sequences as
$\ell^{q,\tau}(L^{p}({\mathcal{T}})):=\\{u_{\mathcal{T}}\,|\,\|u_{\mathcal{T}}\|_{\tau,q,p}<\infty\\}$.
For a fixed spatial basis $\\{\phi_{\ell}({\bf x})\\}_{\ell=1}^{\infty}$, we
can identify $u$ with $u_{\mathcal{T}}$ based on the series representation
(18). Let $\|u\|_{s,q,p}=\|u_{\mathcal{T}}\|_{\tau_{q}(s),q,p}$ with
$\tau_{q}(s)=\frac{s}{d}+\frac{1}{2}-\frac{1}{q}$ as in (5). Then we define
the Banach space of spatiotemporal functions
$B^{s,q,p}({\mathcal{Z}}):=\\{u:{\mathcal{Z}}\to{\mathbb{R}}\,|\,\|u\|_{s,q,p}<\infty\\}$.

Next we generalize BP $u({\bf
x})\sim{\mathcal{B}}(\kappa,B^{s,q}({\mathcal{X}}))$ as in (8) to be
spatiotemporal by letting random coefficients
$\\{\xi_{\ell}\\}_{\ell=1}^{\infty}$ vary in time domain according to
$\mathrm{q}\\!-\\!\mathcal{EP}(0,{\mathcal{C}})$. For this purpose, we make
the following assumption.

######  Assumption 2.

In (18) we let

$u_{\ell}(t)=\gamma_{\ell}\xi_{\ell}(t),\quad\gamma_{\ell}=\kappa\ell^{-\tau_{q}(s)},\quad\xi_{\ell}(\cdot)\overset{i.i.d.}{\sim}\mathrm{q}\\!-\\!\mathcal{EP}(0,{\mathcal{C}}).$
(20)

Compared with (6), we absorb the scaling factor $\kappa>0$ into the covariance
kernel ${\mathcal{C}}$ and set $\kappa=1$ for the most part. We only vary
$\kappa$ in the posterior contraction with adaptive priors in Section 4.2.2.
Under Assumption 2, we have $u$ in (18) as a stochastic process termed
_spatiotemporal Besov process (STBP)_ , denoted as
$\mathcal{STBP}({\mathcal{C}},B^{s,q,p}({\mathcal{Z}}))$. Similarly as in
Section 2.2, the infinite random sequence
$\xi_{\mathcal{T}}:=\\{\xi_{\ell}(\cdot)\\}_{\ell=1}^{\infty}$ is a random
element of the probability space $(\Omega,{\mathcal{B}}(\Omega),\mathbb{P})$
with $\Omega=(L^{p}({\mathcal{T}}))^{\infty}$, product $\sigma$-algebra
${\mathcal{B}}(\Omega)$ and probability measure $\mathbb{P}$ defined by the
infinite product of $\mathrm{q}\\!-\\!\mathcal{EP}(0,{\mathcal{C}})$ measures.
Then we can define a spatiotemporal Besov measure on
$B^{s,q,p}({\mathcal{Z}})$ as the law of STBP.

######  Definition 4 (Spatiotemporal Besov Measure).

Let $\mathbb{P}$ be the measure of random sequences
$\xi_{\mathcal{T}}\in\Omega$. Suppose we have the following map

$f_{\gamma}:\Omega\to B^{s,q,p}({\mathcal{Z}}),\quad\xi_{\mathcal{T}}\mapsto
u({\bf z})=\sum_{\ell=1}^{\infty}u_{\ell}(t)\phi_{\ell}({\bf
x})=\sum_{\ell=1}^{\infty}\gamma_{\ell}\xi_{\ell}(t)\phi_{\ell}({\bf x}),$
(21)

where $\gamma_{\ell}$ and $\xi_{\ell}(\cdot)$ are defined in (20). Then the
pushforward $f_{\gamma}^{\sharp}\mathbb{P}$ is spatiotemporal Besov measure
$\Pi$ on $B^{s,q,p}({\mathcal{Z}})$.

Based on (21), we need to bound $\|u_{\ell}(\cdot)\|_{p}^{q}$ (or
$\|\xi_{\ell}(\cdot)\|_{p}^{q}$) so the norm (19) is well defined. By Theorem
3.3, $\|\xi_{\ell}(\cdot)\|_{p}^{q}$ has bounded mean for $1\leq p\leq q$. For
the convenience of exposition and theoretic development, in the following we
only consider $p=q\in[1,2]$. This is also when most interesting applications
happen (See more details in Section 6). Denote
$\|u\|_{s,q}:=\|u\|_{s,q,q}=\|\xi_{\mathcal{T}}\|_{q,q}=\left(\sum_{\ell=1}^{\infty}\|\xi_{\ell}(\cdot)\|^{q}_{q}\right)^{\frac{1}{q}}$
and $B^{s,q,q}({\mathcal{Z}}):=B^{s,q}({\mathcal{Z}})$. Figure 1 summarizes
the relationship between GP, BP, Q-EP and their spatiotemporal variants.

${p\textrm{-exponential
\cite[citep]{[\@@bibref{Number}{Agapiou_2021}{}{}]}}}$${\mathcal{GP}(0,{\mathcal{C}}_{\bf
x})}$${{\mathcal{B}}(\kappa,B^{s,q}({\mathcal{X}}))\text{\cite[citep]{[\@@bibref{Number}{Lassas_2009}{}{}]}}}$${\mathcal{STBP}({\mathcal{C}}_{t},B^{s,q}({\mathcal{Z}}))}$${\mathcal{STGP}(0,{\mathcal{C}}_{\bf
z})}$${\mathrm{q}\\!-\\!\mathcal{EP}(0,{\mathcal{C}}_{t})\text{\cite[citep]{[\@@bibref{Number}{Li_2023}{}{}]}}}$time-
uncorrelatedmodified
$\pi_{\xi}(\cdot)\propto\exp\left(-\frac{1}{p}|\cdot|^{p}\right)$$q=2$K-L
expansion$\xi_{\ell}\to\xi_{\ell}(\cdot)$$q=2$${\mathcal{C}}_{t}={\mathcal{I}}_{t}$K-L
expansion$q=2$K-L
expansion$\xi_{\ell}(\cdot)\sim\mathrm{q}\\!-\\!\mathcal{EP}(0,{\mathcal{C}}_{t})$
Figure 1: Relationship between GP, BP, Q-EP and their spatiotemporal variants.

###  4.1 STBP as A Prior

In this section, we study multiple properties of STBP used as a prior.
Similarly as in Section 2.2, for a given random draw
$u\sim\mathcal{STBP}({\mathcal{C}},B^{s,q}({\mathcal{Z}}))$ as in (21), we
have
$\mathbb{E}[\|u\|_{s,q}^{q}]=\mathbb{E}[\|\xi_{\mathcal{T}}\|_{q,q}^{q}]=\infty$
due to the iid assumption on $\xi_{\mathcal{T}}$ in (20). However, we have the
following integrability of STBP random sample in the ambient space
$B^{s^{\prime},q}({\mathcal{Z}})$ similar to Theorem 2.1.

######  Theorem 4.1.

If $u\sim\mathcal{STBP}({\mathcal{C}},B^{s,q}({\mathcal{Z}}))$ as in (21)
satisfies both Assumptions 1-(ii) and 2, then $u\in
L^{q}_{\mathbb{P}}(\Omega;B^{s^{\prime},q}({\mathcal{Z}}))$ for all
$s^{\prime}<s-\frac{d}{q}$.

The following (Fernique type) theorem enhances such result of well-definedness
in the context of almost sure convergence for STBP random draws [42, 14, 15].

######  Theorem 4.2.

Let $u\sim\mathcal{STBP}({\mathcal{C}},B^{s,q}({\mathcal{Z}}))$ as in (21)
satisfy both Assumptions 1-(ii) and 2. The following statements are
equivalent:

  1. (i)

$u\in B^{s^{\prime},q}({\mathcal{Z}})$ $\Pi$-a.s.

  2. (ii)

$\mathbb{E}[\exp(\alpha\|u\|_{s^{\prime},q}^{q})]<\infty\;for\;any\;\alpha\in(0,(\sup_{\ell}\lambda_{\ell})^{-\frac{q}{2}}/2)$.

  3. (iii)

$s^{\prime}<s-\frac{d}{q}$.

###### Proof.

See Appendix A.1.

∎

######  Remark 6.

Fernique type result (ii) is important to make rigorous the formal Lebesgue
density (9) and the conditions of similar format in Assumption 4 as well. The
results $(i)$ and $(iii)$ immediately imply that the
$\mathcal{STBP}({\mathcal{C}},B^{s,q}({\mathcal{Z}}))$ measure $\Pi$ is
supported on the ambient space $B^{s^{\prime},q}({\mathcal{Z}})$ for
$s^{\prime}<s-\frac{d}{q}$, as stated in the following corollary.

######  Corollary 4.1.

Let $\Pi$ be an $\mathcal{STBP}({\mathcal{C}},B^{s,q}({\mathcal{Z}}))$ measure
satisfying both Assumptions 1-(ii) and 2. Then
$\Pi(B^{s^{\prime},q}({\mathcal{Z}}))=1$ for any $s^{\prime}<s-\frac{d}{q}$
and $\Pi(B^{s^{\prime},q}({\mathcal{Z}}))=0$ for any $s^{\prime}\geq
s-\frac{d}{q}$.

It is known that $B^{s,q}({\mathcal{Z}})\subset
B^{s^{\prime},q}({\mathcal{Z}})$ for $s^{\prime}<s$. More generally, we have
the following continuous embeddings.

######  Proposition 4.1.

For $q,q^{\dagger}\in[1,2]$ and
$s^{\prime}<s^{\dagger}-\left(\frac{d}{q^{\dagger}}-\frac{d}{q}\right)_{+}$
with $x_{+}:=\max\\{x,0\\}$, we have $B^{s^{\dagger},q^{\dagger}\wedge
q,q}({\mathcal{Z}})\hookrightarrow B^{s^{\prime},q}$, where $q^{\dagger}\wedge
q:=\min\\{q^{\dagger},q\\}$.

###### Proof.

See Appendix A.1. ∎

Similarly to Theorem 3.2, we have the Karhunen-Loéve theorem for an STBP
$u(\cdot)$ as in (21) represented completely in the spatial
($\\{\phi_{\ell}\\}_{\ell=1}^{\infty}$) and temporal
($\\{\psi_{\ell^{\prime}}\\}_{\ell^{\prime}=1}^{\infty}$) bases.

######  Theorem 4.3.

[Karhunen-Loéve] If
$u\sim\mathcal{STBP}({\mathcal{C}},B^{s,q}({\mathcal{Z}}))$ as in (21) with a
trace-class HS operator $T_{\mathcal{C}}$ having eigen-pairs
$\\{\lambda_{\ell},\psi_{\ell}(\cdot)\\}_{\ell=1}^{\infty}$, then we have the
following series representation for $u({\bf z})$:

$u({\bf
z})=\sum_{\ell=1}^{\infty}\sum_{\ell^{\prime}=1}^{\infty}u_{\ell\ell^{\prime}}\phi_{\ell}({\bf
x})\psi_{\ell^{\prime}}(t),\quad
u_{\ell\ell^{\prime}}:=\int_{\mathcal{T}}u_{\ell}(t)\psi_{\ell^{\prime}}(t)dt\sim\mathrm{q}\\!-\\!\mathrm{ED}(0,\gamma_{\ell}^{2}\lambda_{\ell^{\prime}})\,.$
(22)

Moreover, the spatiotemporal covariance of STBP bears a separable structure,
i.e.,

$\mathrm{Cov}(u({\bf z}),u({\bf
z}^{\prime}))=\sum_{\ell=1}^{\infty}\gamma_{\ell}^{2}\phi_{\ell}({\bf
x})\phi_{\ell}({\bf x}^{\prime}){\mathcal{C}}(t,t^{\prime})\,.$ (23)

###### Proof.

See Appendix A.1 ∎

######  Remark 7.

In practice (Section 6), we do not model with such pure series based
representation (22) because it is not straighforward to specify correlations
among sample functions which requires tweaking basis functions according to
(23). Instead, it is our motivation to include a kernel ${\mathcal{C}}$ in
Definition 4 of STBP to directly model the temporal correlations through
covariance functions, e.g. squared exponential, Matérn (See Section 5.3). In
Section 6, we will also investigate the importance of temporal kernel
${\mathcal{C}}$ by comparing against a time-uncorrelated method which sets it
to be identity ${\mathcal{C}}={\mathcal{I}}$ (Refer to Figure 1).

The regularity of an STBP random draw $u({\bf x},t)$ as in (22) also depends
on the properties of spatial ($\\{\phi_{\ell}\\}_{\ell=1}^{\infty}$) and
temporal ($\\{\psi_{\ell^{\prime}}\\}_{\ell^{\prime}=1}^{\infty}$) bases. To
study its Hölder continuity, we make the following assumption that mainly
states that the bases are Hölder continuous with summable Lipschitz constants.

######  Assumption 3.

In the series expansion (22), we assume the spatial
($\\{\phi_{\ell}\\}_{\ell=1}^{\infty}$) and temporal
($\\{\psi_{\ell^{\prime}}\\}_{\ell^{\prime}=1}^{\infty}$) bases satisfy for
$\alpha>0$ and $s^{\prime}<s-\frac{d}{q}$:

  1. (i)

for $\forall{\bf x},{\bf x}^{\prime}\in{\mathcal{X}}$, the following holds

$|\phi_{\ell}({\bf x})-\phi_{\ell}({\bf x}^{\prime})|\leq L(\phi_{\ell})|{\bf
x}-{\bf
x}^{\prime}|^{\alpha+\frac{d}{q}},\quad\sum_{\ell=1}^{\infty}\ell^{\frac{s^{\prime}-s}{d}q}L(\phi_{\ell})<\infty,\quad\sum_{\ell=1}^{\infty}\ell^{\frac{s^{\prime}-s}{d}q}\|\phi_{\ell}(\cdot)\|_{\infty}^{q}<\infty;$

  2. (ii)

for $\forall t,t^{\prime}\in{\mathcal{T}}$, the following holds

$|\psi_{\ell^{\prime}}(t)-\psi_{\ell^{\prime}}(t^{\prime})|\leq
L(\psi_{\ell})|t-t^{\prime}|^{\alpha+\frac{d}{q}},\quad\sum_{\ell^{\prime}=1}^{\infty}\lambda_{\ell^{\prime}}^{\frac{q}{2}}L(\psi_{\ell^{\prime}})<\infty.$

The following theorem regarding the Hölder continuity of STBP random functions
can be proved by the Kolmogorov continuity test [Theorem 30 in Section A.2.5
of 15].

######  Theorem 4.4.

[Hölder Continuity] Let
$u\sim\mathcal{STBP}({\mathcal{C}},B^{s,q}({\mathcal{Z}}))$ as in (21) satisfy
both Assumptions 1-(ii) and 2. Suppose the spatial
($\\{\phi_{\ell}\\}_{\ell=1}^{\infty}$) and temporal
($\\{\psi_{\ell^{\prime}}\\}_{\ell^{\prime}=1}^{\infty}$) bases in the series
representation (22) satisfy Assumption 3. Then for any $\beta<\alpha$, there
exists a version 111A version of stochastic process $u({\bf z})$ is
$\tilde{u}({\bf z})$ such that $\Pi[\tilde{u}({\bf z})=u({\bf z})]=1$ for
$\forall{\bf z}\in{\mathcal{Z}}$. $\tilde{u}({\bf z})$ of $u({\bf z})$ in
$C^{0,\beta}({\mathcal{Z}},B^{s^{\prime},q}({\mathcal{Z}}))$.

###### Proof.

See Appendix A.1. ∎

###  4.2 Posterior Theorems for Inverse Problems with STBP Priors

In this section, we study the posterior properties of the Bayesian inverse
model (1) with STBP priors. In particular, we set the separable Banach space
${\mathbb{X}}=B^{s^{\prime},q}({\mathcal{Z}})$ for some
$s^{\prime}<s-\frac{d}{q}$, and let ${\mathbb{Y}}$ be another separable Banach
space, e.g. ${\mathbb{Y}}=H^{s}({\mathcal{X}})$ or
${\mathbb{Y}}={\mathbb{R}}^{m}$, depending on the applications. Following
[14], we also impose some conditions on the potential (negative log-
likelihood) function
$\Phi:B^{s^{\prime},q}({\mathcal{Z}})\times{\mathbb{Y}}\to{\mathbb{R}}$ as in
(2) regarding its lower (i) and upper (ii) bounds, and the Lipschitz
continuity in $u$ (iii) and $y$ (iv) in the following assumption.

######  Assumption 4.

The potential function
$\Phi:B^{s^{\prime},q}({\mathcal{Z}})\times{\mathbb{Y}}\to{\mathbb{R}}$
satisfies:

  1. (i)

there is an $\alpha_{1}>0$ and for every $r>0$, an $M\in{\mathbb{R}}$, such
that for all $u\in B^{s^{\prime},q}({\mathcal{Z}})$, and for all
$y\in{\mathbb{Y}}$ such that $\|y\|_{\mathbb{Y}}<r$,

$\Phi(u,y)\geq M-\alpha_{1}\|u\|_{s^{\prime},q}\;;$

  2. (ii)

for every $r>0$ there exists $K=K(r)>0$ such that for all $u\in
B^{s^{\prime},q}({\mathcal{Z}}),y\in{\mathbb{Y}}$ with
$\max\\{\|u\|_{s^{\prime},q},\|y\|_{\mathbb{Y}}\\}<r$,

$\Phi(u,y)\leq K\;;$

  3. (iii)

for every $r>0$ there exists $L=L(r)>0$ such that for every $y\in{\mathbb{Y}}$
and for all $u_{1},u_{2}\in B^{s^{\prime},q}({\mathcal{Z}})$ with
$\max\\{\|u_{1}\|_{s^{\prime},q},\|u_{2}\|_{s^{\prime},q}\\}<r$,

$|\Phi(u_{1},y)-\Phi(u_{2},y)|\leq L\|u_{1}-u_{2}\|_{s^{\prime},q}\;;$

  4. (iv)

there is an $\alpha_{2}>0$ and for every $r>0$ a $C\in{\mathbb{R}}$ such that
for all $y_{1},y_{2}\in{\mathbb{Y}}$ with
$\max\\{\|y_{1}\|_{\mathbb{Y}},\|y_{2}\|_{\mathbb{Y}}\\}<r$ and for every
$u\in B^{s^{\prime},q}({\mathcal{Z}})$,

$|\Phi(u,y_{1})-\Phi(u,y_{2})|\leq\exp(\alpha_{2}\|u\|_{s^{\prime},q}+C)\|y_{1}-y_{2}\|_{\mathbb{Y}}\;.$

####  4.2.1 Well-definedness and Well-posedness of Posteriors

Recall that $\Pi$ is the STBP prior defined by (21) and $\Pi^{y}$ is the
resulting posterior measure in the model (1). We re-examine the well-
definedness and well-posedness of the posterior measure based on Assumption 4.

######  Theorem 4.5.

[Well-definedness of Posterior] Let the potential function $\Phi$ in (2)
satisfy Assumption 4 (i)-(iii). If $\Pi$ is an
$\mathcal{STBP}({\mathcal{C}},B^{s,q}({\mathcal{Z}}))$ measure with
$s>s^{\prime}+\frac{d}{q}$, then $\Pi^{y}\ll\Pi$ and satisfies

$\frac{d\Pi^{y}}{d\Pi}(u)=\frac{1}{Z(y)}\exp(-\Phi(u;y)),$ (24)

with the normalizing factor
$0<Z(y)=\int_{B^{s^{\prime},q}({\mathcal{Z}})}\exp(-\Phi(u;y))\Pi(du)<\infty$.

###### Proof.

See Appendix A.2.

∎

Now we show the well-posedness of the posterior measure $\Pi^{y}$ with respect
to the data $y$. Define the Hellinger metric as
$d_{H}(\mu,\mu^{\prime})=\sqrt{\frac{1}{2}\int\left(\sqrt{\frac{d\mu}{d\nu}}-\sqrt{\frac{d\mu^{\prime}}{d\nu}}\right)^{2}d\nu}$.
Note we require $\mu,\mu^{\prime}\ll\nu$, but this definition is independent
of the choice of the measure $\nu$. The following theorem states that the
posterior measure is Lipschitz with respect to data $y$, in the Hellinger
metric.

######  Theorem 4.6.

[Well-posedness of Posterior] Let the potential function $\Phi$ in (2) satisfy
Assumption 4 (i)-(iv). If $\Pi$ is an
$\mathcal{STBP}({\mathcal{C}},B^{s,q}({\mathcal{Z}}))$ measure with
$s>s^{\prime}+\frac{d}{q}$, then

$d_{H}(\Pi^{y},\Pi^{y^{\prime}})\leq C\|y-y^{\prime}\|_{\mathbb{Y}},$

where $C=C(r)$ is a constant depending on $r$ such that
$\max\\{\|y\|_{\mathbb{Y}},\|y^{\prime}\|_{\mathbb{Y}}\\}\leq r$.

###### Proof.

See Appendix A.2.

∎

####  4.2.2 Posterior Contraction of Bayesian Inverse Problems

In this section, we consider the concentration property of the posterior
$\Pi^{y}$ in the infinitely-informative data-limit. Unlike the Gaussian
measure with reproducible kernel Hilbert space (RKHS), the lack of inner
product structure on $B^{s,q}({\mathcal{Z}})\subset
B^{s^{\prime},q}({\mathcal{Z}})$ makes the posterior contraction theories more
challenging [1]. We consider the separable Banach space
$(B^{s^{\prime},q}({\mathcal{Z}}),\|\cdot\|_{s^{\prime},q})$. Define the
concentration function of STBP measure $\Pi$ at $u=u^{\dagger}$ as

$\varphi_{u^{\dagger}}(\varepsilon)=\inf_{h\in
B^{s,q}({\mathcal{Z}}):\|h-u^{\dagger}\|_{s^{\prime},q}\leq\varepsilon}\frac{1}{2}\|h\|_{s,q}^{q}-\log\Pi(\|u\|_{s^{\prime},q}\leq\varepsilon).$
(25)

Let $n=I\wedge J=\min\\{I,J\\}$ with $I,J$ being discrete spatial and temporal
dimensions, respectively. We consider the inverse model (1) with observations
$Y^{(n)}:=\\{Y_{j}=y({\bf X},t_{j})\\}_{j=1}^{n}$ independent but not
identically distributed (the mean ${\mathcal{G}}(u)({\bf X},t_{j})$ changes
for different $j$). Denote $P_{u}^{(n)}:=\bigotimes_{j=1}^{n}P_{u,j}$ as the
product measure on
$\bigotimes_{j=1}^{n}({\mathbb{Y}}_{j},{\mathcal{B}}_{j},\mu_{j})$. Each
$P_{u,j}$ has a density $p_{u_{j}}$ with respect to the $\sigma$-finite
reference measure $\mu_{j}$, which could be
$p_{{\mathcal{N}}({\mathcal{G}}(u_{j}),\sigma^{2}_{\varepsilon})}$. Define the
average Hellinger distance as
$d^{2}_{n,H}(u,u^{\prime})=\frac{1}{n}\sum_{j=1}^{n}\int(\sqrt{p_{u_{j}}}-\sqrt{p_{u^{\prime}_{j}}})^{2}d\mu_{j}$.
Then we have the following posterior contraction theorem.

######  Theorem 4.7.

[Posterior Contraction] Let $u$ be an
$\mathcal{STBP}({\mathcal{C}},B^{s,q}({\mathcal{Z}}))$ random element as in
(21) satisfying both Assumptions 1-(ii) and 2 in
$\Theta:=B^{s^{\prime},q}({\mathcal{Z}})$ with $s^{\prime}<s-\frac{d}{q}$ and
$P_{u}^{(n)}:=\bigotimes_{j=1}^{n}P_{u,j}$ be the product measure of $Y^{(n)}$
parameterized by $u$ with the potential function $\Phi$ (2) satisfying
Assumption 4 (i)-(iii). If the true value $u^{\dagger}\in\Theta$ is in the
support of $u$, and $\varepsilon_{n}$ satisfies the rate equation
$\varphi_{u^{\dagger}}(\varepsilon_{n})\leq n\varepsilon_{n}^{2}$ with
$\varepsilon_{n}\geq n^{-\frac{1}{2}}$, then there exists
$\Theta_{n}\subset\Theta$ such that
$\Pi_{n}(u\in\Theta_{n}:d_{n,H}(u,u^{\dagger})\geq
M_{n}\varepsilon|Y^{(n)})\to 0$ in $P_{u^{\dagger}}^{(n)}$-probability for
every $M_{n}\to\infty$.

###### Proof.

See Appendix A.2.

∎

Denote $a\wedge b=\min\\{a,b\\}$, $a\vee b=\max\\{a,b\\}$, and $x_{+}=x\vee
0$. By solving the inequality $\varphi_{u^{\dagger}}(\varepsilon_{n})\leq
n\varepsilon_{n}^{2}$ for the minimal $\varepsilon_{n}$, we obtain the
posterior contraction rate as follows.

######  Theorem 4.8.

[Posterior Contraction Rate] Let $u$ be an
$\mathcal{STBP}({\mathcal{C}},B^{s,q}({\mathcal{Z}}))$ random element in
$\Theta:=B^{s^{\prime},q}({\mathcal{Z}})$ with $s^{\prime}<s-\frac{d}{q}$. The
rest of the settings are the same as in Theorem 4.7. If the true value
$u^{\dagger}\in B^{s^{\dagger},q^{\dagger}}({\mathcal{Z}})$ with
$s^{\dagger}>s^{\prime}+\left(\frac{d}{q^{\dagger}}-\frac{d}{q}\right)_{+}$
and $q^{\dagger},q\in[1,2]$, then we have the rate of the posterior
contraction as
$\varepsilon_{n}=n^{-\frac{\sigma(s,q,s^{\dagger},q^{\dagger})-s^{\prime}}{2(\sigma(s,q,s^{\dagger},q^{\dagger})-s^{\prime})+q(s-\sigma(s,q,s^{\dagger},q^{\dagger}))}}$,
where
$\sigma(s,q,s^{\dagger},q^{\dagger})=\left(s-\frac{d}{q}\right)\bigwedge\left(s^{\dagger}-\left(\frac{d}{q^{\dagger}}-\frac{d}{q}\right)_{+}\right)$.

###### Proof.

See Appendix A.2.

∎

######  Remark 8.

The contraction rate $\varepsilon_{n}$ becomes optimal
$\varepsilon_{n}^{*}=n^{-\frac{1}{2+\frac{d}{s^{\dagger}-s^{\prime}-\left(\frac{d}{q^{\dagger}}-\frac{d}{q}\right)_{+}}}}$
if
$s=s^{\dagger}+\frac{d}{q}-\left(\frac{d}{q^{\dagger}}-\frac{d}{q}\right)_{+}$,
which is further maximized as
$\varepsilon_{n}^{\dagger}=n^{-\frac{1}{2+\frac{d}{s^{\dagger}-s^{\prime}}}}$
when $q\leq q^{\dagger}$. Note, such optimal rate is achieved regardless of
the value of modeling regularization parameter $q$ as long as $q\leq
q^{\dagger}$. This implies that when modeling inhomogeneous data, under-
smoothing (with smaller regularization parameter $q$) is preferred to over-
smoothing (with larger $q$, refer to Figure 10). This is also the reason why
$q=1$ is often adopted – the posterior can contract the fastest if the true
integrability $q^{\dagger}$ is at least $L_{1}$.

######  Remark 9.

Another observation is that the ambient space
$B^{s^{\prime},q}({\mathcal{Z}})$ can be chosen for the smoothness parameter
$s^{\prime}<\sigma(s,q,s^{\dagger},q^{\dagger})$. In particular, we consider
two cases:

  1. (i)

If we set $\tau_{q}(s^{\prime})=0$, i.e.
$s^{\prime}=\frac{d}{q}-\frac{d}{2}\geq 0$, then
$B^{s^{\prime},q}({\mathcal{Z}})\cong\ell^{q}(L^{q}({\mathcal{T}}))$. For the
Gaussian case ($q^{\dagger}=2$), if we adopt $q=2$ and hence $s^{\prime}=0$
and $B^{0,2}({\mathcal{Z}})\cong\ell^{2}(L^{2}({\mathcal{T}}))\cong
L^{2}({\mathcal{Z}})$, then the optimal rate
$\varepsilon_{n}^{\dagger}=n^{-\frac{s^{\dagger}}{2s^{\dagger}+d}}$ is minimax
[72]. For other sub-Gaussian cases ($q^{\dagger}<2$), such optimal rate
$\varepsilon_{n}^{\dagger}=n^{-\frac{s^{\dagger}-s^{\prime}}{2(s^{\dagger}-s^{\prime})+d}}$
is not minimax [1] regardless of the choice of $q\in[1,2]$ because either
$s^{\prime}>0$ ($q\leq q^{\dagger}<2$) or such optimal rate is not attained
($q>q^{\dagger}$).

  2. (ii)

On the other hand, if we allow $s^{\prime}=0$ and consider a larger ambient
space
$B^{0,q}({\mathcal{Z}})\cong\ell^{q,\frac{1}{2}-\frac{1}{q}}(L^{q}({\mathcal{T}}))\supset\ell^{q}(L^{q}({\mathcal{T}}))$,
then the minimax rate
$\varepsilon_{n}^{\dagger}=n^{-\frac{s^{\dagger}}{2s^{\dagger}+d}}$ can be
obtained for $q\leq q^{\dagger}$.

In general, the minimax posterior contraction rate is not achieved when
$q>q^{\dagger}$. Therefore we typically rescale the prior to infuse additional
regularity [73, 22, 1]. That is, we vary the scaling factor $\kappa>0$ as in
(20) and rescale the Banach space $(\kappa
B^{s,q}({\mathcal{Z}}),\|\cdot\|_{s,q})\cong(B^{s,q}({\mathcal{Z}}),\kappa^{-1}\|\cdot\|_{s,q})$.
Denote
$(B^{s,q}_{\kappa}({\mathcal{Z}}),\|\cdot\|):=(B^{s,q}({\mathcal{Z}}),\kappa^{-1}\|\cdot\|_{s,q})$
and the corresponding (rescaled) STBP measure as $\Pi_{\kappa}$ in Definition
4. Now we redefine the posterior concentration function (25) to be

$\varphi_{u^{\dagger},\kappa}(\varepsilon)=\inf_{h\in
B^{s,q}_{\kappa}({\mathcal{Z}}):\|h-u^{\dagger}\|_{s^{\prime},q}\leq\varepsilon}\frac{\kappa^{-q}}{2}\|h\|_{s,q}^{q}-\log\Pi_{\kappa}(\|u\|_{s^{\prime},q}\leq\varepsilon).$
(26)

The following theorem regards the posterior contraction rate with rescaled
STBP prior.

######  Theorem 4.9.

[Adaptive Posterior Contraction Rate] Let $u$ be an
$\mathcal{STBP}({\mathcal{C}},B^{s,q}_{\kappa}({\mathcal{Z}}))$ random element
in $\Theta:=B^{s^{\prime},q}_{\kappa}({\mathcal{Z}})$ with
$s^{\prime}<s-\frac{d}{q}$. Suppose $\varepsilon_{n}$ satisfies the rate
equation $\varphi_{u^{\dagger},\kappa}(\varepsilon_{n})\leq
n\varepsilon_{n}^{2}$ with $\varepsilon_{n}\geq n^{-\frac{1}{2}}$. The rest of
the settings are the same as in Theorem 4.7. If the true value $u^{\dagger}\in
B^{s^{\dagger},q^{\dagger}}_{\kappa}({\mathcal{Z}})$ with
$s^{\dagger}>s^{\prime}+\left(\frac{d}{q^{\dagger}}-\frac{d}{q}\right)_{+}$
and $1\leq q^{\dagger}<q\leq 2$, then the minimax posterior contraction rate
$\varepsilon_{n}^{\dagger}=n^{-\frac{s^{\dagger}}{2s^{\dagger}+d}}$ can be
attained at
$s=\frac{s^{\dagger}(s^{\dagger}-s^{\prime})}{s^{\prime}+\left(\frac{d}{q^{\dagger}}-\frac{d}{q}\right)}+s^{\prime}$
with the scaling factor $\kappa_{n}\asymp
n^{-\frac{1}{2s^{\dagger}+d}\left[-\frac{s^{\dagger}(s^{\dagger}-s^{\prime})}{s^{\prime}+\left(\frac{d}{q^{\dagger}}-\frac{d}{q}\right)}+\frac{d}{q}+s^{\dagger}\right]}$.

###### Proof.

See Appendix A.2.

∎

##  5 Bayesian Inference

In this section, we describe the inference of the Bayesian inverse problem (1)
with spatiotemporal observations using an STBP prior. Assume the unknown
function $u$ is evaluated at $I$ locations ${\bf X}:=\\{{\bf
x}_{i}\\}_{i=1}^{I}$ and $J$ time points ${\bf t}:=\\{t_{j}\\}_{j=1}^{J}$,
that is $u({\bf X},{\bf t}):=\\{u({\bf x}_{i},t_{j})\\}_{i,j=1}^{I,J}$. The
data ${\bf Y}=\\{{\bf y}_{j}\\}_{j=1}^{J}$ with ${\bf
y}_{j}\in{\mathbb{R}}^{m}$ is observed through the forward operator
${\mathcal{G}}$, which could be a linear mapping or a nonlinear one governed
by a PDE. Here we consider Gaussian noise and recap the model as follows:

$\displaystyle{\bf y}_{j}$ $\displaystyle={\mathcal{G}}(u)({\bf
X},t_{j})+\bm{\varepsilon}_{j},\quad\bm{\varepsilon}_{j}\overset{iid}{\sim}{\mathcal{N}}_{I}(0,\Gamma_{\mathrm{noise}}),\quad
j=1,2,\dots,J,$ (27) $\displaystyle u$
$\displaystyle\sim\mathcal{STBP}({\mathcal{C}},B^{s,q}({\mathcal{Z}})).$

In our applications of inverse problems, the spatial dimension $I$ is usually
tremendously higher than the temporal dimension ($I\gg J$). Therefore we
truncate $u$ in (21) for the first $L>0$ terms: $u({\bf x},t)\approx
u^{L}({\bf x},t)=\sum_{\ell=1}^{L}\gamma_{\ell}\xi_{\ell}(t)\phi_{\ell}({\bf
x})$, and choose $L=2000$ in the numerical experiments (Section 6). Denote
${\bf u}_{j}=u({\bf X},t_{j})\in{\mathbb{R}}^{I}$, and ${\bf U}=[{\bf
u}_{1},\cdots,{\bf u}_{J}]_{I\times J}=u^{L}({\bf X},{\bf
t})={\bm{\Phi}}\operatorname{diag}({\bm{\gamma}}){{\bm{\Xi}}}^{\mathsf{T}}$
where ${\bm{\Phi}}=[\phi_{1}({\bf X}),\cdots,\phi_{L}({\bf X})]_{I\times L}$,
${\bm{\gamma}}=(\gamma_{1},\cdots,\gamma_{L})$, and
${\bm{\Xi}}=[\bm{\xi}_{1},\cdots,\bm{\xi}_{L}]_{J\times L}$ with
$\bm{\xi}_{\ell}=\xi_{\ell}({\bf t})$. Instead of the large dimensional matrix
${\bf U}$, we work with ${\bm{\Xi}}$ of much smaller size. Let
$r_{\ell}={\bm{\xi}_{\ell}}^{\mathsf{T}}{\bf C}_{J}^{-1}\bm{\xi}_{\ell}$. The
log-posterior for ${\bm{\Xi}}$ is computed directly as

$\displaystyle\log p({\bm{\Xi}},\theta|{\bf Y})=$
$\displaystyle-\frac{J}{2}\log|\Gamma_{\mathrm{noise}}|-\frac{1}{2}\sum_{j=1}^{J}\|{\bf
y}_{j}-{\mathcal{G}}({\bf u}_{j})\|^{2}_{\Gamma_{\mathrm{noise}}}$ (28)
$\displaystyle-\frac{L}{2}\log|{\bf
C}_{J}|+\frac{J}{2}(\frac{q}{2}-1)\sum_{\ell=1}^{L}\log
r_{\ell}-\frac{1}{2}\sum_{\ell=1}^{L}r_{\ell}^{\frac{q}{2}}.$

We optimize (28) to obtain the maximum a posterior (MAP) estimate. To quantify
the uncertainty efficiently, we need effective inference algorithms for high-
dimensional models with non-Gaussian priors. We refer to the work of
dimension-robust MCMC proposed by [13] based on the pushforward of Gaussian
white noise which in turn takes advantage of the dimension-independent
sampling algorithms for Gaussian priors [8]. For the convenience of
application, in the following, we introduce a new white noise representation
for STBP which is different from the one used in [13] for series based priors.

###  5.1 White Noise Representation

Recall we have the stochastic representation (11) of
$\bm{\xi}\sim\mathrm{q}\\!-\\!\mathrm{ED}_{J}({\bf 0},{\bf C})$:
$\bm{\xi}=R{\bf L}S$ with $R^{q}\sim\chi^{2}(J)$ and
$S\sim\mathrm{Unif}({\mathcal{S}}^{J+1})$. We can write

$S=\frac{\bm{\zeta}}{\|\bm{\zeta}\|_{2}},\quad
R^{q}=\|\bm{\zeta}\|_{2}^{2},\quad\textrm{for}\;\bm{\zeta}\sim{\mathcal{N}}_{J}({\bf
0},{\bf I}_{J}).$

Therefore, $\bm{\xi}$ can be represented in terms of the white noise
$\bm{\zeta}$ by a pushforward mapping
$\Lambda:{\mathbb{R}}^{J}\to{\mathbb{R}}^{J}$ and vice versa with its inverse
$\Lambda^{-1}$:

$\bm{\xi}=\Lambda(\bm{\zeta})={\bf
L}\bm{\zeta}\|\bm{\zeta}\|_{2}^{\frac{2}{q}-1},\quad\bm{\zeta}=\Lambda^{-1}(\bm{\xi})={\bf
L}^{-1}\bm{\xi}\|{\bf L}^{-1}\bm{\xi}\|_{2}^{\frac{q}{2}-1}.$ (29)

For
$\zeta(\cdot)=\sum_{\ell^{\prime}=1}^{\infty}\zeta_{\ell^{\prime}}\psi_{\ell^{\prime}}(\cdot),\;\xi(\cdot)=\sum_{\ell^{\prime}=1}^{\infty}\xi_{\ell^{\prime}}\psi_{\ell^{\prime}}(\cdot)\in
L^{2}({\mathcal{T}})$, we can extend $\Lambda$ and its inverse $\Lambda^{-1}$
to $L^{2}({\mathcal{T}})$ and have

$\xi(\cdot)=\Lambda(\zeta(\cdot))=\sum_{\ell^{\prime}=1}^{\infty}\lambda_{\ell^{\prime}}^{\frac{1}{2}}\zeta_{\ell^{\prime}}\psi_{\ell^{\prime}}(\cdot)\|\zeta(\cdot)\|_{2}^{\frac{2}{q}-1},\quad\zeta(\cdot)=\Lambda^{-1}(\xi(\cdot))=\sum_{\ell^{\prime}=1}^{\infty}\lambda_{\ell^{\prime}}^{-\frac{1}{2}}\xi_{\ell^{\prime}}\psi_{\ell^{\prime}}(\cdot)\|\xi(\cdot)\|_{2,{\mathcal{C}}}^{\frac{2}{q}-1}\,,$

where
$\|\zeta(\cdot)\|_{2}^{2}=\sum_{\ell^{\prime}=1}^{\infty}\zeta_{\ell^{\prime}}^{2}$
and
$\|\xi(\cdot)\|_{2,{\mathcal{C}}}^{2}=\sum_{\ell^{\prime}=1}^{\infty}\lambda_{\ell^{\prime}}^{-1}\xi_{\ell^{\prime}}^{2}$.
We propose the following representation of $u({\bf z})$ in terms of an
infinite sequence of white noises, i.e.
$\zeta:=\\{\zeta_{\ell}(\cdot)\\}_{\ell=1}^{\infty}$:

$u({\bf
z})=T(\zeta)=\sum_{\ell=1}^{\infty}\gamma_{\ell}\Lambda(\zeta_{\ell}(t))\phi_{\ell}({\bf
x}),\quad\zeta_{\ell}(\cdot)\overset{i.i.d.}{\sim}\mathcal{GP}(0,{\mathcal{I}}).$
(30)

Denote $\bm{\mathrm{Z}}=[\bm{\zeta}_{1},\cdots,\bm{\zeta}_{L}]_{J\times L}$
with $\bm{\zeta}_{\ell}=\zeta_{\ell}({\bf t})$. From (30), we have ${\bf
U}=T(\bm{\mathrm{Z}})={\bm{\Phi}}\operatorname{diag}({\bm{\gamma}}){\Lambda(\bm{\mathrm{Z}})}^{\mathsf{T}}$.
Then the log-posterior (28) can be rewritten in terms of $\bm{\mathrm{Z}}$:

$\displaystyle\log p(\bm{\mathrm{Z}},\theta|{\bf Y})=$
$\displaystyle-\frac{J}{2}\log|\Gamma_{\mathrm{noise}}|-\frac{1}{2}\sum_{j=1}^{J}\mathrm{tr}({({\bf
Y}-{\mathcal{G}}(T(\bm{\mathrm{Z}})))}^{\mathsf{T}}\Gamma_{\mathrm{noise}}^{-1}({\bf
Y}-{\mathcal{G}}(T(\bm{\mathrm{Z}}))))$ (31)
$\displaystyle-\frac{L}{2}\log|{\bf
C}_{J}|+\frac{J}{2}(\frac{q}{2}-1)\sum_{\ell=1}^{L}\log
r_{\ell}-\frac{1}{2}\sum_{\ell=1}^{L}r_{\ell}^{\frac{q}{2}},$

where $r_{\ell}={\Lambda(\bm{\zeta}_{\ell})}^{\mathsf{T}}{\bf
C}_{J}^{-1}\Lambda(\bm{\zeta}_{\ell})$. Once the MAP
$\bm{\mathrm{Z}}_{\text{\tiny MAP}}$ is obtained by maximizing the above log-
posterior, we can obtain ${\bf U}_{\text{\tiny
MAP}}=T(\bm{\mathrm{Z}}_{\text{\tiny MAP}})$. We refer to this process as
“optimization in the whitened space”. The objective function can be explored
more efficiently in the whitened space with variables de-correlated (See
Figure 5).

###  5.2 White Noise MCMC

Denote the measure formed by infinite product of
$\mathcal{GP}(0,{\mathcal{I}})$ as $\Pi_{0}$. Then our STBP prior measure
$\Pi$ can be regained by the pushforward using $T$, i.e.
$\Pi=T^{\sharp}\Pi_{0}$. A class of dimension-independent MCMC algorithms [8]
for models with Gaussian prior $\Pi_{0}$ can be reintroduced to posterior
sampling with STBP prior $\Pi$.

Let $u=T(\zeta)$ with $\zeta\sim\Pi_{0}$. Consider the following continuous-
time Hamiltonian dynamics:

$\frac{d^{2}\zeta}{dt^{2}}+\mathcal{K}(\zeta)\,[\zeta+\nabla\Phi(\zeta)]=0,\quad\left.\left(\eta:=\frac{d\zeta}{dt}\right)\right|_{t=0}\sim{\mathcal{N}}(0,\mathcal{K}(\zeta))\
,$ (32)

where $\Phi(\zeta):=\Phi(T(\zeta))-\log|dT(\zeta)|$. More generally, we could
set $\mathcal{K}(\zeta)^{-1}={\mathcal{I}}+\beta\mathcal{H}(\zeta)$ where
$\mathcal{H}(\zeta)$ can be chosen as Hessian, Gauss-Newton Hessian, or Fisher
information operator [40]. For example, we can choose the Gauss-Newton Hessian
computed as $\mathcal{H}(\zeta)=dT^{*}\mathcal{H}(u)dT$ with $dT$ being the
Jacobian. Let
$g(\zeta):=-\mathcal{K}(\zeta)\\{\alpha\nabla\Phi(\zeta)-\beta\mathcal{H}(\zeta)\zeta\\}$
where
$\nabla_{\zeta}\Phi(\zeta)=dT^{*}\nabla_{u}\Phi(u)-\nabla_{\zeta}\log|dT(\zeta)|$.
The Hamiltonian Monte Carlo (HMC) algorithm [53] solves the dynamics (32)
using the Störmer-Verlet symplectic (leapfrog) integrator with step size
$\varepsilon$:

$\displaystyle\eta^{-}$
$\displaystyle=\eta_{0}+\tfrac{\varepsilon}{2}\,g(\zeta_{0})\ ;$ (33)
$\displaystyle\begin{bmatrix}\zeta_{\varepsilon}\\\ \eta^{+}\end{bmatrix}$
$\displaystyle=\begin{bmatrix}\cos\varepsilon&\sin\varepsilon\\\
-\sin\varepsilon&\cos\varepsilon\end{bmatrix}\begin{bmatrix}\zeta_{0}\\\
\eta^{-}\end{bmatrix}\ ;$ $\displaystyle\eta_{\varepsilon}$
$\displaystyle=\eta^{+}+\tfrac{\varepsilon}{2}\,g(\zeta_{\varepsilon})\ .$

Equation (33) gives rise to the leapfrog map
$\Psi_{\varepsilon}:(\zeta_{0},\eta_{0})\mapsto(\zeta_{\varepsilon},\eta_{\varepsilon})$.
Given a time horizon $\tau$ and current position $\zeta$, the MCMC mechanism
proceeds by concatenating $I=\lfloor\tau/\varepsilon\rfloor$ steps of leapfrog
map consecutively,

$\zeta^{\prime}=\mathcal{P}_{\zeta}\big{\\{}\Psi_{\varepsilon}^{I}(\zeta,\eta)\big{\\}}\
,\quad\eta\sim{\mathcal{N}}(0,\mathcal{K}(\zeta))\ ,$

where $\mathcal{P}_{\zeta}$ denotes the projection onto the $\zeta$-argument.
Then, the proposal $\zeta^{\prime}$ is accepted with probability
$a(\zeta,\zeta^{\prime})=1\wedge\exp(-\Delta E(\zeta,\eta))$, where

$\displaystyle\Delta E(\zeta,\eta)=$ $\displaystyle
E(\Psi_{\varepsilon}^{I}(\zeta,\eta))-E(\zeta,\eta)$ (34) $\displaystyle=$
$\displaystyle\Phi(\zeta_{I})-\Phi(\zeta_{0})+\frac{\beta}{2}\langle\eta_{I},\mathcal{H}(\zeta_{I})\eta_{I}\rangle-\frac{\beta}{2}\langle\eta_{0},\mathcal{H}(\zeta_{0})\eta_{0}\rangle-\log|\mathcal{K}^{-\frac{1}{2}}(\zeta_{I})|+\log|\mathcal{K}^{-\frac{1}{2}}(\zeta_{0})|$
$\displaystyle-\frac{\varepsilon^{2}}{8}\left\\{\|g(\zeta_{I})\|^{2}-\|g(\zeta_{0})\|^{2}\right\\}-\frac{\varepsilon}{2}\sum_{i=0}^{I-1}(\langle
g(\zeta_{i}),\eta_{i}\rangle+\langle g(\zeta_{i+1}),\eta_{i+1}\rangle).$

At last we convert the sample $\zeta$ back to $u=T(\zeta)$. This yields white-
noise infinite-dimensional manifold HMC (wn-$\infty$-mHMC) [8] which reduces
to white-noise infinite-dimensional manifold Metropolis adjusted Langevin
algorithm (wn-$\infty$-mMALA) [7] when $I=1$, and white-noise infinite-
dimensional HMC (wn-$\infty$-HMC) [6] when $\beta=0$. We set $\alpha=1$ for
both scenarios and summarize all the these methods in Algorithm 1 of Appendix
B named as _white-noise dimension-independent MCMC (wn- $\infty$-MCMC)_.

###  5.3 Hyper-parameter Tuning

In Definitions 3 and 4, there is a temporal kernel ${\mathcal{C}}$ that has
not been specified. This is the key component to capture the temporal
correlation which is absent in a pure series based approach (See Remark 7 and
Section 6). There are hyper-parameters, denoted as $\theta$, in the covariance
kernel ${\mathcal{C}}$, e.g. variance magnitude ($\kappa$) and correlation
length ($\rho$), i.e. $\theta=(\kappa,\rho)$, that require careful adjustment
and fine tuning as in, e.g., Matérn kernel:

${\mathcal{C}}(t,t^{\prime})=\kappa\frac{2^{1-\nu}}{\Gamma(\nu)}w^{\nu}K_{\nu}(w),\quad
w=\sqrt{2\nu}(\|t-t^{\prime}\|/\rho)^{s}.$ (35)

Unless we assume the likelihood in the Bayesian inverse model (1) to be
another $\mathrm{q}\\!-\\!\mathrm{ED}$ and the forward mapping is linear [c.f.
Theorem 3.5 of 47], we do not have a tractable marginal likelihood to optimize
for these hyper-parameters [61]. In general settings, e.g., in the model (27)
with a Gaussian likelihood, we need to jointly update $({\bm{\Xi}},\theta)$
based on (28) or $(\bm{\mathrm{Z}},\theta)$ according to (31). Denote ${\bf
C}_{0}=\kappa^{-1}{\bf C}$ and $r_{0,\ell}={\bm{\xi}_{\ell}}^{\mathsf{T}}{\bf
C}_{0}^{-1}\bm{\xi}_{\ell}$. The following proposition permits conditional
conjugacy for the variance magnitude ($\kappa$) given an appropriate hyper-
prior.

######  Proposition 5.1.

If we assume a inverse-gamma hyper-prior for the variance magnitude
$\kappa^{\frac{q}{2}}\sim\Gamma^{-1}(\alpha,\beta)$ such that
$\bm{\xi}_{\ell}|\kappa\overset{iid}{\sim}\mathrm{q}\\!-\\!\mathrm{ED}_{J}({\bf
0},{\bf C})$ in (27), then we have

$\kappa^{\frac{q}{2}}|{\bf
u}\sim\Gamma^{-1}(\alpha^{\prime},\beta^{\prime}),\quad\alpha^{\prime}=\alpha+\frac{JL}{2},\quad\beta^{\prime}=\beta+\frac{1}{2}\sum_{\ell=1}^{L}r_{0,\ell}^{\frac{q}{2}}$
(36)

###### Proof.

See Appendix B. ∎

Therefore, we could either update
$\kappa\leftarrow\left(\frac{\beta^{\prime}}{\alpha^{\prime}+1}\right)^{\frac{2}{q}}$
or sample $\kappa$ according to (36). In general, there is no such conditional
conjugacy for the correlation length ($\rho$). We impose a hyper-prior for
$\rho$ and sample from $p(\rho|{\bm{\Xi}})$. Alternatively, one could also use
Bayesian optimization methods [46, 18, 5] for tuning $\rho$.

##  6 Numerical Experiments

In this section, we compare the proposed STBP
($\mathcal{STBP}({\mathcal{C}},B^{s,q}({\mathcal{Z}}))$) with STGP (equivalent
to $\mathcal{STBP}({\mathcal{C}},B^{s,2}({\mathcal{Z}}))$) and a time-
uncorrelated prior ($\mathcal{STBP}({\mathcal{I}},B^{s,q}({\mathcal{Z}}))$)
using a simulated regression, two dynamic tomography imaging examples and an
inverse problem of recovering a spatiotemporal function. Since the main focus
is to model inhomogeneous data such as images with edges, we tend to adopt
sharper regularization and set $q=1$ for STBP throughout this section (See
also Remark 8). Our numerical results demonstrate the advantage of Besov
($L_{1}$) type priors over Gaussian ($L_{2}$) type priors in modeling
inhomogeneous data. Moreover, these examples highlight the importance of
appropriately modeling temporal correlations in spatiotemporal inverse
problems. All the computer codes are publicly available at
https://github.com/lanzithinking/Spatiotemporal-Besov-prior.

In all these applications, $u({\bf x}_{i},t_{j})$ refers to the image pixel
value of point ${\bf x}_{i}$ at time $t_{j}$ with resolution $I=n_{x}\times
n_{y}$. To asses the quality of reconstructed images (viewed as inverse
solutions $u({\bf x},t)$ defined on 2d space, i.e.
${\mathcal{X}}\subset{\mathbb{R}}^{2}$), we refer to several quantitative
measures including the relative error,
$\text{RLE}=\frac{\|u^{*}-u^{\dagger}\|}{\|u^{\dagger}\|}$, where
$u^{\dagger}$ denotes the reference/true image and $u^{*}$ is its
reconstruction. Additionally, we adopt the peak signal-to-noise ratio,
$\text{PSNR}=10*\log_{10}(\frac{\|u^{\dagger}\|_{\infty}^{2}}{\|u^{*}-u^{\dagger}\|_{2}^{2}})$,
by using the maximum possible pixel value as a reference point to normalise
the MSE. We could also consider the structured similarity index [77],
$\text{SSIM($u^{*},u^{\dagger}$)}=\frac{(2\overline{u^{*}}\overline{u^{\dagger}}+c_{1})(2s_{u^{*}u^{\dagger}}+c_{2})}{(\overline{u^{*}}^{2}+\overline{u^{\dagger}}^{2}+c_{1})(s^{2}_{u^{*}}+s^{2}_{u^{\dagger}}+c_{2})}$,
where $\overline{u}$, $s^{2}_{u}$ and $s_{u_{1}u_{2}}$ denote the sample mean,
sample variance, and sample covariance respectively, $c_{i}=(k_{i}L)^{2}$ for
$i=1,2$, $k_{1}=0.01$, $k_{2}=0.03$ and $L$ is the dynamic range of the pixel
values of the reference images. These image quality metrics are extensively
used to benchmark algorithms in alignment with human perception, ensuring
greater consistency.

Truth |  Observations |  STBP |  STGP |  time-uncorrelated
---|---|---|---|---

---

---

---

---

Figure 2: MAP reconstruction of simulated annulus with $I=256\times 256$,
$J=100$. Columns from left to right: true images, observations, MAP estimates
by STBP, STGP and time-uncorrelated models respectively. Rows from top to
bottom: time step $t_{j}=0.1,0.3,0.6,0.9$.

###  6.1 Simulation

First, we consider a simulated regression problem of a rising and shrinking 2d
annulus:

$u({\bf x},t)=t\delta(\sin(\pi\|{\bf x}\|_{2})\geq t),\quad{\bf
x}\in{\mathbb{R}}^{2}\;such\;that\;\|{\bf x}\|_{\infty}\leq 1,\quad
t\in(0,1],$ (37)

where $\delta(\cdot)$ is the Dirac function returning 1 when the condition
(${\bf x}$ is inside the 2d annulus) is met and 0 otherwise. The first column
of Figure 2 plots this function at a few time points illustrating a 2d annulus
forming and shrinking as time goes by. We simulate the data by discretizing
the function $u({\bf x},t)$ in (37) on a $I=n_{x}\times n_{y}$ mesh (denoted
as ${\bf X}$) in the boxed spatial domain ${\mathcal{X}}=[-1,1]^{2}$ over a
grid of $J$ time points (denoted as ${\bf t}$) in ${\mathcal{T}}=(0,1]$, and
adding some Gaussian noise with $\sigma_{\varepsilon}=0.1$, i.e.

${\bf y}_{j}=u({\bf
X},t_{j})+\bm{\varepsilon}_{j},\quad\bm{\varepsilon}_{j}\overset{iid}{\sim}{\mathcal{N}}_{I}(0,\sigma^{2}_{\varepsilon}{\bf
I}).$

The noisy spatiotemporal data are demonstrated in the second column of Figure
2. Based on the observed data, the goal of this Bayesian inverse problem to
recover the ground truth (37) using STBP, STGP and time-uncorrelated priors.
We use this example to numerically investigate the posterior contraction
studied in Section 4.2.2 as $n=I\wedge J\to\infty$. In particular, we will
consider the problem with data observed at various spatiotemporal resolutions
by combining $I=16\times 16,32\times 32,128\times 128,256\times 256$ and
$J=10,20,50,100$ respectively.

Table 1: Comparison of MAP estimates for simulated annulus generated by STBP, STGP and time-uncorrelated prior models in terms of RLE with increasing data. Standard deviations (in bracket) are obtained by repeating the experiments for 10 times with different random seeds for initialization. $I=n_{x}\times n_{y}$ | $J$ | time-uncorrelated | STGP | STBP
---|---|---|---|---
| 10 |  0.2768 (5.32e-6) |  0.2632 (7.91e-6) |  0.2770 (4.14e-6)
$16\times 16$ | 20 |  0.2399 (3.52e-6) |  0.2063 (5.97e-6) |  0.2421 (6.10e-6)
| 50 |  0.2067 (1.02e-6) |  0.1594 (1.08e-5) |  0.1553 (1.07e-5)
| 100 |  0.1850 (8.12e-7) |  0.1217 (8.60e-6) |  0.1211 (7.79e-6)
| 10 |  0.3677 (3.00e-6) |  0.2465 (7.08e-6) |  0.3672 (6.81e-6)
$32\times 32$ | 20 |  0.2460 (1.49e-6) |  0.1911 (4.88e-6) |  0.2359 (5.02e-6)
| 50 |  0.2236 (1.30e-6) |  0.1464 (7.27e-6) |  0.1662 (7.36e-6)
| 100 |  0.1918 (9.90e-7) |  0.1203 (7.73e-6) |  0.1241 (9.31e-6)
| 10 |  0.1943 (1.57e-5) |  0.2052 (4.51e-6) |  0.1937 (1.36e-5)
$128\times 128$ | 20 |  0.1474 (7.23e-6) |  0.1497 (4.16e-6) |  0.1399 (1.77e-5)
| 50 |  0.1227 (4.95e-6) |  0.1083 (1.03e-5) |  0.1030 (1.51e-5)
| 100 |  0.1146 (2.31e-6) |  0.0934 (1.51e-5) |  0.0909 (1.24e-5)
| 10 |  0.1635 (1.47e-5) |  0.1970 (2.70e-6) |  0.1633 (1.03e-5)
$256\times 256$ | 20 |  0.1190 (7.24e-6) |  0.1442 (7.18e-6) |  0.1088 (1.24e-5)
| 50 |  0.0966 (5.77e-6) |  0.1012 (1.46e-5) |  0.0858 (1.31e-5)
| 100 |  0.0910 (3.89e-6) |  0.0864 (6.68e-6) |  0.0808 (1.36e-5)

Note, the spatial image of the function at each time point, when viewed as a
picture, has clear edges. This imposes challenges for GP as it tends to over-
smooth when modeling inhomogeneous objects while BP is more amenable. On the
other hand, these sequential images are not isolated from each other in time,
and the temporal kernel ${\mathcal{C}}$ in STBP (STGP) can be well-used to
capture such dependence. More specifically, we adopt the Matérn kernel (35)
with $\nu=\frac{1}{2}$, $\sigma^{2}=1$, $\rho=0.1$ and $s=1$. The MAP estimate
for ${\bf U}=u({\bf X},{\bf t})$ is obtained by minimizing the negative log-
posterior (31) in the whitened space of $\bm{\mathrm{Z}}$ and converting
$\bm{\mathrm{Z}}_{\text{\tiny MAP}}$ back to ${\bf U}_{\text{\tiny
MAP}}=T(\bm{\mathrm{Z}}_{\text{\tiny MAP}})$. The last three columns of Figure
2 compare the MAP estimates by STBP, STGP and time-uncorrelated models at
$I=256\times 256$, $J=100$. The STGP model indeed returns an over-smoothed
result; while the time-uncorrelated model yields a more noisy estimate due to
the negligence of temporal correlation. Figure 10 also demonstrates different
degrees of regularization interpolating with parameter $q$ in the range of
$(0,2]$ with $q=2$ (STGP) yielding the most blurry solution.

Next, we vary the spatiotemporal resolution by changing the mesh and the time
interval for observations. Figure 3 investigates the MAP estimates by the STBP
model with increasing data. They gradually approximate the true function (37)
as the spatiotemporal resolution is refined. This verifies the posterior
consistency described in Theorem 4.7 in terms of point estimation. Table 1
also shows an error reducing phenomenon with increasing data for all three
models. To make a fair comparison across different resolutions, we adopt
$\|{\bf U}\|_{\infty,1}=\max_{1\leq i\leq I}\sum_{j=1}^{J}|u({\bf
x}_{i},t_{j})|$ in the RLE to focus on the pixel differences while averaging
over the time domain. The STBP model outperforms the other two in most cases.
Though not a direct verification of the posterior contraction rate Theorem
4.8, it shows that STBP reduces error with increasing data at rates not slower
than STGP.

Truth |  $\overset{I=16\times 16}{J=10}$ |  $\overset{I=32\times 32}{J=20}$ |  $\overset{I=128\times 128}{J=50}$ |  $\overset{I=256\times 256}{J=100}$
---|---|---|---|---

---

---

---

---

Figure 3: MAP reconstruction of simulated annulus by STBP model with
increasing data. Columns from left to right: true images, MAP estimates
obtained at different spatiotemporal resolutions. Rows from top to bottom:
time step $t_{j}=0.1,0.3,0.6,0.9$.

###  6.2 Dynamic Tomography Reconstruction

In this section, we investigate the dynamic reconstruction of a simulated
(STEMPO) and a real (emoji) tomography problems. Computed tomography (CT) is a
medical imaging technique used to non-intrusively obtain detailed internal
images of a subject such as human body [65, 28]. CT scanners project (Radon
transformation) X-ray over the subject at different angles and measure the
attenuated signals by an array of sensors recorded as sinograms.
Reconstructing the (high-resolution) internal images from the limited observed
sinograms is usually regarded as a linear inverse problem which is severely
underdetermined (ill-posed). Therefore, appropriate (edge-preserving) priors
are especially important for these tomography reconstruction problems.

Truth |  Observations |  STBP |  STGP |  time-uncorrelated
---|---|---|---|---

---

---

---

---

Figure 4: MAP reconstruction of dynamic STEMPO tomography in the whitened
space. Columns from left to right: true images, sinograms, MAP estimates by
STBP, STGP and time-uncorrelated models respectively. Rows from top to bottom:
time step $j=0,6,13,19$.

####  6.2.1 STEMPO Tomography

In the first example, we consider a simulated dynamic tomography named Spatio-
TEmporal Motor-POwered (STEMPO) ground truth phantom from [27]. The dataset
contains 360 snapshots of CT images each of size $I=560\times 560$ but we
choose $J=20$ at equal time intervals. For each of $J=20$ snapshot, the ASTRA
toolbox [See more details in 71] is used to generate the forward operators
${\mathcal{G}}_{j}$, $j=1,2,\dots,J$. Each operator
${\mathcal{G}}_{j}\in{\mathbb{R}}^{8701\times 313600}$ then projects the true
image $u^{\dagger}({\bf X},t_{j})$ to a sinogram
${\mathcal{G}}_{j}(u^{\dagger}({\bf X},t_{j}))\in{\mathbb{R}}^{791\times 11}$
at $n_{a}=11$ angles ($x$-axis) with $n_{s}=791$ equally spaced X-ray
detectors ($y$-axis) as shown in the second column of Figure 4. Finally, as in
the model (27), Gaussian white noise
$\bm{\varepsilon}_{j}\sim{\mathcal{N}}_{n_{a}n_{s}}(0,\sigma_{j}^{2}{\bf
I}_{n_{a}n_{s}})$ is added to the signogram to obtain the observation ${\bf
Y}_{j}={\mathcal{G}}_{j}(u^{\dagger}({\bf X},t_{j}))+\bm{\varepsilon}_{j}$
such that the noise level
$\|\sigma_{j}\|_{2}/\|{\mathcal{G}}_{j}\left(u^{\dagger}({\bf
X},t_{j}))\right)\|_{2}=0.01$. The true images $\\{u^{\dagger}({\bf
X},t_{j})\\}_{j=1}^{J}$ and the noisy observations $\\{{\bf
Y}_{j}\\}_{j=1}^{J}$ at time $j=0,6,13,19$ are shown in the first two columns
of Figure 4, respectively.

Figure 5: Dynamic STEMPO tomography reconstruction: negative posterior densities and relative errors for the optimization in the original space (left two) and in the whitened space (right two) as functions of iterations in the BFGS algorithm used to obtain MAP estimates. Early termination is implemented if the error falls below some threshold or the maximal iteration (1000) is reached. Table 2: Comparison of MAP estimates for STEMPO tomography generated by STBP, STGP and time-uncorrelated prior models in terms of RLE, log-likelihood, PSNR, and SSIM measures. Standard deviations (in bracket) are obtained by repeating the experiments for 10 times with different random seeds for initialization. | time-uncorrelated | STGP | STBP
---|---|---|---
RLE | 0.4354 (2.91e-5) | 0.3512(1.42e-4) | 0.3217 (2.72e-5)
log-likelihood | -39190.72 (0.65) | -39085.37 (5.49) | -39697.93 (0.71)
PSNR | 16.6235 (5.80e-4) | 18.4896 (3.50e-3) | 19.2532 (7.33e-4)
SSIM | 0.1469 (3.50e-5) | 0.3486 (3.47e-4) | 0.2318 (7.10e-5)

We minimize the negative log-posterior densities (31) in terms of the whitened
coordinates $\bm{\mathrm{Z}}$ for the three models with STBP, STGP and time-
uncorrelated priors to obtain the MAP estimates. The rightmost three columns
of Figure 4 compare these MAP estimates obtained in the whitened space and
mapped to the original space. With STBP, we have the sharpest reconstruction
on the third column that is the closest to the truth. However, the results by
the other two models are either blurry (by STGP on the forth column) or noisy
(by the time-uncorrelated model on the last column). Table 2 confirms that the
STBP model yields the reconstruction of the highest quality. Though their log-
likelihood values are not comparable in the regularized optimization, the STBP
model achieves the lowest $\textrm{RLE}=32.17\%$ on average in 10 experiments
repeated with different random seeds. The same advantage is supported by high
values in other image quality measures such as PSNR, and SSIM.

Truth |  STBP (mean) |  STGP (mean) |  STBP (std) |  STGP (std)
---|---|---|---|---

---

---

---

---

Figure 6: MCMC reconstruction of dynamic STEMPO tomography in the whitened
space. Columns from left to right: true images, posterior mean estimates by
STBP and STGP, posterior standard deviation estimates by STBP and STGP models
respectively. Rows from top to bottom: time step $j=0,6,13,19$.

On the other hand, the MAP estimates generated by minimizing the negative log-
posterior (28) in terms of the original parameters ${\bm{\Xi}}$ are compared
in Figure 11. They have more than $40\%$ RLE’s and are generally more blurry
than those obtained in the whitened space (See Figure 4). Such difference is
also observed in Figure 5 where the objective functions and RLE’s are compared
between optimization in the original space (w.r.t. ${\bm{\Xi}}$, left two
panels) and optimization in the whitened space (w.r.t. ${\bf Z}$, right two
panels) for these three models: the latter yields better results within fewer
iterations bearing lower errors, possibly due to faster exploration in the
whitened space with variables de-correlated. In general, STBP converges faster
to the lowest error state among the three models.

Lastly, we apply the white-noise manifold infinite-dimensional MALA
(wn-$\infty$-mMALA) (Algorithm 1) to sample from the posterior distributions
of the two models with STBP and STGP priors respectively (the result for time-
uncorrelated prior is far worse and hence omitted) and compare their posterior
estimates in Figure 6. We generate 3000 samples and discard the first 1000
samples. The remaining 2000 samples are used to estimate the posterior means
(the second and the third columns) and posterior standard deviations (the last
two columns). Due to the large dimensionality ($560\times 560\times 20$) and
limited number of samples, these posterior estimates tend to be noisy. The
posterior mean estimates are not as good reconstructions as their MAP
estimates. Yet the posterior standard deviations by STBP (the forth column)
provides uncertainty information with more clear spatial features than those
by STGP model (the last column).

Truth |  Observations |  STBP |  STGP |  time-uncorrelated
---|---|---|---|---

---

---

---

---

Figure 7: MAP reconstruction for the dynamic emoji tomography in the whitened
space. Columns from left to right: true images, sinograms, MAP estimates by
STBP, STGP and time-uncorrelated models respectively. Rows from top to bottom:
time step $j=6,14,22,30$.

####  6.2.2 Emoji Tomography

Next, we test our methods on a real data of dynamic “emoji” phantom measured
at the University of Helsinki [See more details in 51, about the machine
(forward operator) set-up and data collection]. The available spatiotemporal
data represent $J=33$ time steps of a series of the X-ray sinogram of emojis
made of small ceramic stones obtained by shining $n_{s}=217$ X-ray projections
from $n_{a}=10$ angles.

The inverse problem involves reconstructing a sequence of images $u({\bf
X},t_{j})$, $t=1,2,\dots,J$, each of size $I=128\times 128$, from low-dose
observations measured at a limited number of $n_{a}$ angles. Hence, the
unknown images are collected in ${\bf U}=u({\bf X},{\bf
t})\in{\mathbb{R}}^{16,384\times 33}$, representing the dynamic sequence of
the emoji images changing from an expressionless face with closed eyes and a
straight mouth to a face with smiling eyes and mouth, where the outmost
circular shape does not change. We refer to Figure 7 for a sample of 4 setup
images (first column) and sinograms (second column) at time steps
$t=6,14,22,30$. The low-dose observations are modeled as in the model (27):
${\bf Y}_{j}\sim{\mathcal{N}}_{2170}({\mathcal{G}}_{j}(u^{\dagger}({\bf
X},t_{j})),\Gamma_{\mathrm{noise}})$ with $\Gamma_{\mathrm{noise}}$ being the
empirical covariance obtained from $J=33$ images, and measurement matrix
${\mathcal{G}}_{j}$ being the result of the same Radon transform as above that
represents line integrals [51]. Although the ground truth is not available, we
can qualitatively compare the visual outputs from STBP, STGP and time-
uncorrelated models.

Figure 7 compares the MAP estimates by STBP (the third column), STGP (the
forth row) and the time-uncorrelated (the last column) prior models in the
whitened space. Again we observe similar advantage in reconstructing a
sequence of sharper tomography images by STBP compared with those more blurry
results by STGP. Note, due to the absence of temporal correlation, the time-
uncorrelated prior model reconstructs images that are noisy and difficult to
recognize. We also compare the UQ results generated by wn-$\infty$-mMALA (See
Algorithm 1) for the two models, STBP and STGP, respectively in Figure 13.
Again we observe noisy posterior mean estimates (the second and the third
columns) for both models. However, the posterior standard deviation estimates
(the forth column) by STBP are slightly clearer than those (the last column)
by STGP in characterizing the uncertainty field representing the changing
smiling faces.

Truth |  Observations |  STBP |  STGP |  time-uncorrelated
---|---|---|---|---

---

---

---

---

Figure 8: MAP inverse solutions of Navier-Stokes equation in the whitened
space. Columns from left to right: true vorticity, observations, MAP estimates
by STBP, STGP and time-uncorrelated models respectively. Rows from top to
bottom: time step $j=0,3,6,9$.

###  6.3 Navier-Stokes Inverse Problem

Lastly, we consider a complex non-linear inverse problem involving the
following 2-d Navier-Stokes equation (NSE) for a viscous, incompressible fluid
in vorticity form on the unit torus ${\mathbb{T}}^{2}=(0,1)^{2}$:

$\displaystyle\begin{split}\partial_{t}w(x,t)+u(x,t)\cdot\nabla
w(x,t)&=\nu\Delta w(x,t)+f(x),\qquad x\in(0,1)^{2},t\in(0,T],\\\ \nabla\cdot
u(x,t)&=0,\qquad\qquad\qquad\qquad\quad x\in(0,1)^{2},t\in[0,T],\\\
w(x,0)&=w_{0}(x),\qquad\qquad\qquad\quad x\in(0,1)^{2}.\end{split}$

where $u\in C([0,T];H^{r}({\mathbb{T}}^{2};{\mathbb{R}}^{2}))$ for any $r>0$
is the velocity field, $w=\nabla\times u$ is the vorticity, $w_{0}\in
L^{2}({\mathbb{T}}^{2};{\mathbb{R}})$ is the initial vorticity,
$\nu\in{\mathbb{R}}_{+}$ is the viscosity coefficient, and $f\in
L^{2}({\mathbb{T}}^{2};{\mathbb{R}})$ is the forcing function.

Because NSE is computationally intensive to solve, we build an emulator based
on the Fourier operator neural network (FNO) [48, 38] that maps the vorticity
up to time $T_{0}=10$ to the vorticity up to some later time $T>10$:

${\mathcal{G}}:C([0,T_{0}];H^{r}({\mathbb{T}}^{2};{\mathbb{R}}^{2}))\to
C((T_{0},T];H^{r}({\mathbb{T}}^{2};{\mathbb{R}}^{2})),\quad
w|_{(0,1)^{2}\times[0,10]}\mapsto w|_{(0,1)^{2}\times(10,T]}$

One of the attractive features of FNO is that the neural network is built to
learn operators defined on function spaces. Compared with traditional neural
networks for simulating PDE solutions including CNN [29, 75], PINNs [60], and
ResNet [26], FNO is mesh-independent and very efficient for the inference of
Bayesian inverse problem constrained by NSE.

In this example, we choose the viscosity $\nu=1\mathrm{e}-3$ and set
$T-T_{0}=30$. Since the target operator, ${\mathcal{G}}^{\dagger}$, is time-
dependent, we train a 3-d FNO (FNO-3d) based on 5000 pairs of input vorticity
(for the first 10 unit time) and output vorticity (for the following 30 unit
time) solved on $I=64\times 64$ spatial mesh (denoted as ${\bf X}$) using the
network configuration the same as in [48]. We initialize the vorticity $w_{0}$
with a (star-convex) polygon shown as in the top left of Figure 8 which also
demonstrates a few snapshots of true vorticity trajectory,
$w^{\dagger}|_{(0,1)^{2}\times[0,10]}({\bf X},t_{j})$, at $j=0,3,6,9$ in the
first column. We then observe data of vorticity
$w|_{(0,1)^{2}\times(10,40]}({\bf X},t_{j})$ with $t_{j}\in(T_{0},T]$ for
$j=0,\cdots,29$ based on the true initial inputs
$w^{\dagger}|_{(0,1)^{2}\times[0,10]}$, with Gaussian noise contamination,
i.e., ${\bf y}_{j}={\mathcal{G}}(w^{\dagger}|_{(0,1)^{2}\times[0,10]}({\bf
X},t_{j}))+\bm{\eta}_{j}$ with $\bm{\eta}_{j}\sim
N(0,\Gamma_{\mathrm{noise}})$, and $\Gamma_{\mathrm{noise}}$ empirically
estimated as in the previous example. A few time snapshots of the observed
vorticity are illustrated in the second column of Figure 8. Figure 14 compares
the emulated NSE trajectory by the FNO network (lower row) against the true
trajectory solved by the classical PDE solver (upper row) in the observation
time window $(T_{0},T]$. Their visual difference is hardly discernible.

Figure 9: Navier-Stokes inverse problem: negative posterior densities (left) and relative errors (right) for the optimization in the whitened space as functions of iterations in the BFGS algorithm used to obtain MAP estimates. Early termination is implemented if the error falls below some threshold or the maximal iteration (1000) is reached. Table 3: Comparison of MAP estimates of NSE trajectory generated by STBP, STGP and time-uncorrelated prior models in terms of RLE, log-likelihood, PSNR, and SSIM measures. Standard deviations (in bracket) are obtained by repeating the experiments for 10 times with different random seeds for initialization. | time-uncorrelated | STGP | STBP
---|---|---|---
RLE | 0.7656 (7.60e-5) | 0.7457 (3.041e-5) | 0.6618 (1.07e-4)
log-lik | -229.18 (0.21) | -1586.11 (0.47) | -173.33 (0.10)
PSNR | 15.7267 (8.62e-4) | 15.9555 (3.54e-4) | 16.9921 (1.40e-3)
SSIM | 0.1842 (7.84e-5) | 0.2213 (5.68e-5) | 0.3416 (1.67e-4)

Unlike the traditional time-dependent inverse problems seeking the solution of
the initial condition $w_{0}$ alone, we are interested in the inverse solution
of vorticity for an initial period, i.e., $w|_{(0,1)^{2}\times[0,10]}$. What
is more, we want to obtain UQ for such spatiotemporal object in addition to
its point estimate (MAP) using STBP, STGP and time-uncorrelated priors. Their
MAP estimates are compared in the last three columns of Figure 8. Note, this
inverse problem for a spatiotemporal solution is much more challenging than
the traditional inverse problem for just the initial condition based on the
same amount of downstream observations. STBP still yields the inverse solution
(the third column) closest to the true trajectory (the first column) among the
three models, especially the initial condition at $t=0$ (the first row) which
is the most difficult because it is the farthest from the observation window
$(T_{0},T]$. Note, due to the lack of temporal correlation, the solution
trajectory from the time-uncorrelated prior model appears excessively erratic.
Table 3 further confirms that STBP prior model yields the best inverse
solution with the lowest RLE, $66.18\%$, compared with the true trajectory,
almost $10\%$ lower than the other two methods. The high values of image
reconstruction metrics also support the superiority of STBP model compared
with the other two. Figure 9 compares the optimization objective (the negative
log-posterior) and the relative error as functions of iterations. STBP
converges to lower RLE value, while STGP terminates earlier at a higher RLE
value.

Lastly, because of the computational cost, we apply white-noise infinite-
dimensional HMC (wn-$\infty$-HMC) (Algorithm 1) instead of wn-$\infty$-mMALA
for the UQ. We run 20,000 iterations, discard the first 5,000, and sub-sample
one of every three. The remaining 5,000 samples are used to obtain posterior
estimates illustrated in Figure 15 comparing STBP model (the second and the
forth columns) against STGP model (the third and the last columns). The
posterior mean by STBP (the second column) is more noisy compared with that by
STGP (the third row). However the posterior standard deviation by STBP (the
forth column) is more informative than that of STGP (the last column).

##  7 Conclusion

In this paper we propose a nonparametric Bayesian framework for spatiotemporal
inverse problems with inhomogeneous data, such as sequential images with
edges. Our proposed STBPs are generalizations of BP from spatial to
spatiotemporal domain. The key idea is to replace the random coefficients
(following a $q$-exponential distribution) in the series definition of BP with
the recently proposed $Q$-EP [47] to account for the temporal correlations
among spatial function images through a covariance kernel, similar to GP.
Moreover, STBP controls the regularization of the posterior solutions through
a parameter $q\in[1,2]$ and includes STGP as a special case ($q=2$).

We conduct a thorough theoretical investigation regarding well-definedness,
Karhunen-Loéve series representation of STBP priors and their posterior
contraction property to justify the appropriateness and superiority of the
proposed methodology. To address the challenges of posterior inference, we
propose dimension-independent MCMC algorithms based on a new white-noise
representation for series based priors [13]. Through extensive numerical
experiments from various spatiotemporal inverse problems we demonstrate that
STBP ($q=1$) has advantage in handling spatial inhomogeity over STGP (which
tends to be over-smooth) and in capturing temporal correlations over a time-
uncorrelated approach. These promising results from limited-angle CT
reconstruction suggest potential applications in medical imaging analysis.

Several directions remain open for future research. Currently, the constraint
of the STBP prior on a common spatial ($q$) and temporal ($p$) regularization
parameter, i.e., $q=p$ (Section 4) can be relaxed to independently control the
spatial and temporal regularities of the posterior solution – this might help
with slow changing spatial data while $L^{2}({\mathcal{T}})$ ($p=2$) suffices
the need for the temporal space. We also aim to extend the current work to the
non-convex regime for $q\in(0,1)$. It imposes more regularization than the
adopted case ($q=1$) (refer to Figure 10) however the development is expected
to be considerably more complex.

## Acknowledgments

SL is supported by NSF grant DMS-2134256. MP gratefully acknowledges support
from the NSF award 2202846.

## References

  * [1] Sergios Agapiou, Masoumeh Dashti, and Tapio Helin.  Rates of contraction of posterior distributions based on p-exponential priors.  Bernoulli, 27(3):1616 – 1642, 2021.
  * [2] Yves A. Atchadaé.  On the contraction properties of some high-dimensional quasi-posterior distributions.  The Annals of Statistics, 45(5), October 2017.
  * [3] Yves Atchade and Anwesha Bhattacharyya.  An approach to large-scale quasi-bayesian inference with spike-and-slab priors.  arXiv:1803.10282, 03 2018.
  * [4] Frank Aurzada.  On the lower tail probabilities of some random sequences in lp.  Journal of Theoretical Probability, 20(4):843–858, Dec 2007.
  * [5] Noor Awad, Neeratyoy Mallik, and Frank Hutter.  DEHB: Evolutionary hyberband for scalable, robust and efficient hyperparameter optimization.  In Proceedings of the Thirtieth International Joint Conference on Artificial Intelligence. International Joint Conferences on Artificial Intelligence Organization, aug 2021.
  * [6] 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\.
  * [7] Alexandros Beskos.  A stable manifold MCMC method for high dimensions.  Statistics & Probability Letters, 90:46–52, 2014.
  * [8] 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.
  * [9] Daniela Calvetti, Monica Pragliola, Erkki Somersalo, and Alexander Strang.  Sparse reconstructions from few noisy data: analysis of hierarchical Bayesian models with generalized gamma hyperpriors.  Inverse Problems, 36(2):025010, 2020.
  * [10] Stamatis Cambanis, Steel Huang, and Gordon Simons.  On the theory of elliptically contoured distributions.  Journal of Multivariate Analysis, 11(3):368–385, 1981.
  * [11] Carlos M Carvalho, Nicholas G Polson, and James G Scott.  The horseshoe estimator for sparse signals.  Biometrika, 97(2):465–480, 2010.
  * [12] George Casella, Malay Ghosh, Jeff Gill, and Minjung Kyung.  Penalized regression, standard errors, and Bayesian Lassos.  Bayesian analysis, 5(2):369–411, 2010.
  * [13] Victor Chen, Matthew M Dunlop, Omiros Papaspiliopoulos, and Andrew M Stuart.  Dimension-robust mcmc in bayesian inverse problems.  arXiv preprint arXiv:1803.03344, 2018.
  * [14] Masoumeh Dashti, Stephen Harris, and Andrew Stuart.  Besov priors for bayesian inverse problems.  Inverse Problems and Imaging, 6(2):183–200, may 2012.
  * [15] Masoumeh Dashti and Andrew M. Stuart.  The Bayesian Approach to Inverse Problems, pages 311–428.  Springer International Publishing, Cham, 2017.
  * [16] Matthew M Dunlop and Andrew M Stuart.  Map estimators for piecewise continuous inversion.  Inverse Problems, 32(10):105003, 2016.
  * [17] Heinz Werner Engl, Martin Hanke, and Andreas Neubauer.  Regularization of inverse problems, volume 375.  Springer Science & Business Media, 1996.
  * [18] Stefan Falkner, Aaron Klein, and Frank Hutter.  BOHB: Robust and efficient hyperparameter optimization at scale.  In Jennifer Dy and Andreas Krause, editors, Proceedings of the 35th International Conference on Machine Learning, volume 80 of Proceedings of Machine Learning Research, pages 1437–1446. PMLR, 10–15 Jul 2018\.
  * [19] K. Fang and Y.T. Zhang.  Generalized Multivariate Analysis.  Science Press, 1990.
  * [20] Chao Gao, Aad W. van der Vaart, and Harrison H. Zhou.  A general framework for bayes structured linear models.  The Annals of Statistics, 48(5), October 2020.
  * [21] S. Ghosal and A.W. van der Vaart.  Convergence rates of posterior distributions for non-i.i.d. observations.  Annals of Statistics, 35(1):192–223, 2007.  MR2332274.
  * [22] Subhashis Ghosal and Aad van der Vaart.  Fundamentals of nonparametric bayesian inference.  2017\.
  * [23] Evarist Giné and Richard Nickl.  Rates of contraction for posterior distributions in $L^{r}$-metrics, $1\leq r\leq\infty$.  The Annals of Statistics, 39(6), December 2011.
  * [24] Matteo Giordano.  Besov-laplace priors in density estimation: optimal posterior contraction rates and adaptation.  Electronic Journal of Statistics, 17(2), January 2023.
  * [25] M. Hairer, A. M. Stuart, and J. Voss.  Analysis of spdes arising in path sampling part ii: The nonlinear case.  The Annals of Applied Probability, 17(5–6), October 2007.
  * [26] Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun.  Deep Residual Learning for Image Recognition.  2015\.
  * [27] Tommi Heikkilä.  Stempo–dynamic X-ray tomography phantom.  arXiv preprint arXiv:2209.12471, 2022.
  * [28] Young M. Hermena S.  Ct-scan image production procedures.  StatPearls Publishing, August 2023.
  * [29] Geoffrey Hinton.  Imagenet classification with deep convolutional neural networks.  NIPS’12: Proceedings of the 25th International Conference on Neural Information Processing Systems - Volume 1, 1:1097–1105, 2012.
  * [30] Marco A Iglesias, Yulong Lu, and Andrew M Stuart.  A bayesian level set method for geometric inverse problems.  Interfaces and free boundaries, 18(2):181–217, 2016.
  * [31] Hemant Ishwaran and J. Sunil Rao.  Spike and slab variable selection: Frequentist and bayesian strategies.  The Annals of Statistics, 33(2), April 2005.
  * [32] Valentin Konstantinovich Ivanov.  On linear problems which are not well-posed.  In Doklady akademii nauk, volume 145, pages 270–272. Russian Academy of Sciences, 1962.
  * [33] Junxiong Jia, Jigen Peng, and Jinghuai Gao.  Bayesian approach to inverse problems for functions with a variable-index Besov prior.  Inverse Problems, 32(8):085006, 2016.
  * [34] Mark E. Johnson.  Multivariate Statistical Simulation, chapter 6 Elliptically Contoured Distributions, pages 106–124.  Probability and Statistics. John Wiley & Sons, Ltd, 1987.
  * [35] Y. Kano.  Consistency property of elliptic probability density functions.  Journal of Multivariate Analysis, 51(1):139–147, 1994.
  * [36] W Clem Karl.  Regularization in image restoration and reconstruction.  In Handbook of Image and Video Processing, pages 183–V. Elsevier, 2005.
  * [37] Bartek T Knapik, Aad W Van Der Vaart, and J Harry van Zanten.  Bayesian inverse problems with gaussian priors.  The Annals of Statistics, 39(5):2626–2657, 2011.
  * [38] Nikola Kovachki, Zongyi Li, Burigede Liu, Kamyar Azizzadenesheli, Kaushik Bhattacharya, Andrew Stuart, and Anima Anandkumar.  Neural operator: Learning maps between function spaces with applications to pdes.  Journal of Machine Learning Research, 24(89):1–97, 2023.
  * [39] Tomasz J. Kozubowski, Krzysztof Podgórski, and Igor Rychlik.  Multivariate generalized laplace distribution and related random fields.  Journal of Multivariate Analysis, 113:59–72, 2013.  Special Issue on Multivariate Distribution Theory in Memory of Samuel Kotz.
  * [40] Shiwei Lan.  Adaptive dimension reduction to accelerate infinite-dimensional geometric markov chain monte carlo.  Journal of Computational Physics, 392:71 – 95, September 2019.
  * [41] Shiwei Lan.  Learning temporal evolution of spatial dependence with generalized spatiotemporal gaussian process models.  Journal of Machine Learning Research, 23(259):1–53, 2022.
  * [42] Matti Lassas, Eero Saksman, and Samuli Siltanen.  Discretization-invariant bayesian inversion and besov space priors.  Inverse Problems and Imaging, 3(1):87–122, 2009.
  * [43] Matti Lassas and Samuli Siltanen.  Can one use total variation prior for edge-preserving Bayesian inversion?  Inverse Problems, 20(5):1537, 2004.
  * [44] L. LeCam.  Convergence of estimates under dimensionality restrictions.  The Annals of Statistics, 1(1):38–53, Jan 1973.
  * [45] J Lee and PK Kitanidis.  Bayesian inversion with total variation prior for discrete geologic structure identification.  Water Resources Research, 49(11):7658–7669, 2013.
  * [46] Lisha Li, Kevin Jamieson, Giulia DeSalvo, Afshin Rostamizadeh, and Ameet Talwalkar.  Hyperband: A novel bandit-based approach to hyperparameter optimization.  J. Mach. Learn. Res., 18(1):6765–6816, jan 2017.
  * [47] Shuyi Li, Michael O’Connor, and Shiwei Lan.  Bayesian learning via q-exponential process.  In Proceedings of the 37th Conference on Neural Information Processing Systems. NeurIPS, 12 2023.
  * [48] Zongyi Li, Nikola Borislavov Kovachki, Kamyar Azizzadenesheli, Burigede liu, Kaushik Bhattacharya, Andrew Stuart, and Anima Anandkumar.  Fourier neural operator for parametric partial differential equations.  In International Conference on Learning Representations, 2021.
  * [49] David Madigan and Adrian E. Raftery.  Model selection and accounting for model uncertainty in graphical models using occam’s window.  Journal of the American Statistical Association, 89(428):1535–1546, December 1994.
  * [50] Markku Markkanen, Lassi Roininen, Janne MJ Huttunen, and Sari Lasanen.  Cauchy difference priors for edge-preserving Bayesian inversion.  Journal of Inverse and Ill-posed Problems, 27(2):225–240, 2019\.
  * [51] Alexander Meaney, Zenith Purisha, and Samuli Siltanen.  Tomographic x-ray data of 3d emoji.  arXiv preprint arXiv:1802.09397, 2018.
  * [52] T. J. Mitchell and J. J. Beauchamp.  Bayesian variable selection in linear regression.  Journal of the American Statistical Association, 83(404):1023–1032, December 1988.
  * [53] R. M. Neal.  MCMC using Hamiltonian dynamics.  In S. Brooks, A. Gelman, G. Jones, and X. L. Meng, editors, Handbook of Markov Chain Monte Carlo. Chapman and Hall/CRC, 2010.
  * [54] Bernt Øksendal.  Stochastic Differential Equations.  Springer Berlin Heidelberg, 2003.
  * [55] Trevor Park and George Casella.  The Bayesian Lasso.  Journal of the American Statistical Association, 103(482):681–686, 2008.
  * [56] Mirjeta Pasha, Arvind K. Saibaba, Silvia Gazzola, Malena I. Espa nol, and Eric de Sturler.  A computational framework for edge-preserving regularization in dynamic inverse problems.  Electron. Trans. Numer. Anal., 58:486–516, 2023.
  * [57] Juho Piironen and Aki Vehtari.  Sparsity information and regularization in the horseshoe and other shrinkage priors.  Electronic Journal of Statistics, 11(2), January 2017.
  * [58] Krzysztof Podgórski and Jörg Wegener.  Estimation for stochastic models driven by laplace motion.  Communications in Statistics - Theory and Methods, 40(18):3281–3302, sep 2011.
  * [59] Nicholas G. Polson, James G. Scott, and Jesse Windle.  The Bayesian Bridge.  Journal of the Royal Statistical Society Series B: Statistical Methodology, 76(4):713–733, 11 2013.
  * [60] M. Raissi, P. Perdikaris, and G.E. Karniadakis.  Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations.  Journal of Computational Physics, 378:686–707, 2019.
  * [61] Carl Edward Rasmussen and Christopher K. I. Williams.  Gaussian Processes for Machine Learning.  The MIT Press, 2005.
  * [62] Vincent Rivoirard and Judith Rousseau.  Posterior concentration rates for infinite dimensional exponential families.  Bayesian Analysis, 7(2), June 2012.
  * [63] Uwe Schmidt, Kevin Schelten, and Stefan Roth.  Bayesian deblurring with integrated noise estimation.  In CVPR 2011, pages 2625–2632. IEEE, 2011.
  * [64] I. J. Schoenberg.  Metric spaces and completely monotone functions.  Annals of Mathematics, 39:811–841, 1938.
  * [65] L. A. Shepp and B. F. Logan.  The fourier reconstruction of a head section.  IEEE Transactions on Nuclear Science, 21(3):21–43, 1974.
  * [66] Roel Snieder and Jeannot Trampert.  Inverse problems in geophysics.  In Wavefield inversion, pages 119–190. Springer, 1999.
  * [67] Andrew M Stuart.  Inverse problems: a Bayesian perspective.  Acta Numerica, 19:451–559, 2010.
  * [68] Jarkko Suuronen, Neil K. Chada, and Lassi Roininen.  Cauchy markov random field priors for bayesian inversion.  Statistics and Computing, 32(2):33, Mar 2022.
  * [69] Andrey Nikolayevich Tikhonov.  On the stability of inverse problems.  In Dokl. Akad. Nauk SSSR, volume 39, pages 195–198, 1943.
  * [70] Hans Triebel.  Theory of Function Spaces.  Springer Basel, 1983.
  * [71] Wim Van Aarle, Willem Jan Palenstijn, Jan De Beenhouwer, Thomas Altantzis, Sara Bals, K Joost Batenburg, and Jan Sijbers.  The ASTRA Toolbox: A platform for advanced algorithm development in electron tomography.  Ultramicroscopy, 157:35–47, 2015.
  * [72] A. W. van der Vaart and J. H. van Zanten.  Rates of contraction of posterior distributions based on gaussian process priors.  The Annals of Statistics, 36(3):1435–1463, 2008.
  * [73] A. W. van der Vaart and J. H. van Zanten.  Adaptive bayesian estimation using a gaussian random field with inverse gamma bandwidth.  Ann. Statist., 37(5B):2655–2675, 10 2009.
  * [74] Aad van der Vaart and Harry van Zanten.  Information rates of nonparametric gaussian process methods.  J. Mach. Learn. Res., 12:2095–2119, July 2011.
  * [75] Ragav Venkatesan and Baoxin Li.  Convolutional Neural Networks in Visual Computing: A Concise Guide.  CRC Press, 2017.
  * [76] Sebastian J Vollmer.  Posterior consistency for bayesian inverse problems through stability and regression results.  Inverse Problems, 29(12):125011, November 2013.
  * [77] Zhou Wang, Alan C Bovik, Hamid R Sheikh, and Eero P Simoncelli.  Image quality assessment: From error visibility to structural similarity.  IEEE transactions on image processing, 13(4):600–612, 2004.
  * [78] Zhewei Yao, Zixi Hu, and Jinglai Li.  A tv-gaussian prior for infinite-dimensional bayesian inverse problems and its numerical implementations.  Inverse Problems, 32(7):075006, May 2016.

## Appendix

##  Appendix A Proofs

###  A.1 Theorems of BP, STBP Priors

See 2.1

###### Proof of Theorem 2.1.

Based on (6), it is straightforward to verify

$\mathbb{E}[\|u\|_{s^{\prime},q}^{q}]=\sum_{\ell=1}^{\infty}\ell^{\tau_{q}(s^{\prime})q}\mathbb{E}|u_{\ell}|^{q}=\kappa^{-1}\mathbb{E}[|\xi_{1}|^{q}]\sum_{\ell=1}^{\infty}\ell^{(\tau_{q}(s^{\prime})-\tau_{q}(s))q}<\infty$

if $(\tau_{q}(s^{\prime})-\tau_{q}(s))q=\frac{s^{\prime}-s}{d}q<-1$, i.e.
$s^{\prime}<s-\frac{d}{q}$. ∎

See 3.3

###### Proof of Theorem 3.3.

Note
$r(\xi_{\ell})^{\frac{q}{2}}=\lambda_{\ell}^{-\frac{q}{2}}|\xi_{\ell}|^{q}\sim\chi^{2}(1)$
for all $\ell\in{\mathbb{N}}$ by Proposition 3.1. Denote
$\chi^{2}_{\ell}:=\lambda_{\ell}^{-\frac{q}{2}}|\xi_{\ell}|^{q}\overset{iid}{\sim}\chi^{2}(1)$.
Hence
$\|\xi\|_{s^{\prime},q}^{q}=\sum_{\ell=1}^{\infty}\ell^{\tau_{q}(s^{\prime})q}\lambda_{\ell}^{\frac{q}{2}}\chi^{2}_{\ell}$
becomes an infinite mixture of chi-squared random variables whose density is
analytically intractable. Yet we have

$\mathbb{E}[\|\xi(\cdot)\|_{s^{\prime},q}^{q}]=\sum_{\ell=1}^{\infty}\ell^{\tau_{q}(s^{\prime})q}\lambda_{\ell}^{\frac{q}{2}}\mathbb{E}[\chi^{2}_{\ell}]=\sum_{\ell=1}^{\infty}\ell^{\tau_{q}(s^{\prime})q}\lambda_{\ell}^{\frac{q}{2}}<\infty$

if Assumption 1-(i) holds. From above we have
$\mathbb{E}[\|\xi(\cdot)\|_{q}^{q}]=\|\lambda\|_{\frac{q}{2}}^{\frac{q}{2}}<\infty$
if Assumption 1-(ii) holds. ∎

See 4.2

###### Proof of Theorem 4.2.

First, by Assumption-2 and the proof of Theorem 3.3 we have

$\|u\|_{s^{\prime},q}^{q}=\sum_{\ell=1}^{\infty}\ell^{\tau_{q}(s^{\prime})q}\|u_{\ell}(\cdot)\|_{q}^{q}=\sum_{\ell=1}^{\infty}\ell^{(\tau_{q}(s^{\prime})-\tau_{q}(s))q}\|\xi_{\ell}(\cdot)\|_{q}^{q}=\sum_{\ell=1}^{\infty}\ell^{\frac{s^{\prime}-s}{d}q}\sum_{\ell^{\prime}=1}^{\infty}\lambda_{\ell^{\prime}}^{\frac{q}{2}}\chi^{2}_{\ell\ell^{\prime}}$

where $\chi^{2}_{\ell\ell^{\prime}}\overset{iid}{\sim}\chi^{2}(1)$. Denote
$\alpha_{\ell\ell^{\prime}}=\alpha\ell^{\frac{s^{\prime}-s}{d}q}\lambda_{\ell^{\prime}}^{\frac{q}{2}}$.
Then we have

$\mathbb{E}[\exp(\alpha_{\ell\ell^{\prime}}\chi^{2}_{\ell\ell^{\prime}})]=M_{\chi^{2}(1)}(\alpha_{\ell\ell^{\prime}})=\left[1-2\alpha_{\ell\ell^{\prime}}\right]^{-\frac{1}{2}},\quad
for\;\alpha_{\ell\ell^{\prime}}<\frac{1}{2}$

#####  $(iii)\implies(ii)$.

Now that

$\mathbb{E}[\exp(\alpha\|u\|_{s^{\prime},q}^{q})]=\mathbb{E}\left[\exp\left(\sum_{\ell=1}^{\infty}\sum_{\ell^{\prime}=1}^{\infty}\alpha_{\ell\ell^{\prime}}\chi^{2}_{\ell\ell^{\prime}})\right)\right]=\prod_{\ell=1}^{\infty}\prod_{\ell^{\prime}=1}^{\infty}\left[1-2\alpha_{\ell\ell^{\prime}}\right]^{-\frac{1}{2}}$

Assume $\alpha>0$, we have each item in the product bigger than 1. To bound
such infinite product, it suffices to bound the following infinite sum

$\sum_{\ell=1}^{\infty}\sum_{\ell^{\prime}=1}^{\infty}\left(\left[1-2\alpha_{\ell\ell^{\prime}}\right]^{-\frac{1}{2}}-1\right)=\sum_{\ell=1}^{\infty}\sum_{\ell^{\prime}=1}^{\infty}\frac{2\alpha_{\ell\ell^{\prime}}}{\sqrt{1-2\alpha_{\ell\ell^{\prime}}}+1-2\alpha_{\ell\ell^{\prime}}}\leq\frac{2\alpha\sum_{\ell=1}^{\infty}\ell^{\frac{s^{\prime}-s}{d}q}\|\lambda\|_{\frac{q}{2}}^{\frac{q}{2}}}{\sqrt{1-2\alpha\|\lambda\|_{\infty}^{\frac{q}{2}}}+1-2\alpha\|\lambda\|_{\infty}^{\frac{q}{2}}}$

By Assumptions 1-(ii), it is finite provided that $\frac{s^{\prime}-s}{d}q<-1$
and $1-2\alpha\|\lambda\|_{\infty}^{\frac{q}{2}}>0$, that is,
$s^{\prime}<s-\frac{d}{q}$ and
$\alpha<\frac{1}{2}\|\lambda\|_{\infty}^{-\frac{q}{2}}$ where
$\|\lambda\|_{\infty}=\sup_{\ell}\lambda_{\ell}$.

#####  $(ii)\implies(i)\implies(iii)$.

The proof remains the same as in [42, 15]. ∎

See 4.1

###### Proof of Proposition 4.1.

If $q^{\dagger}=q$, we have $B^{s^{\dagger},q^{\dagger}}({\mathcal{Z}})\subset
B^{s^{\prime},q}({\mathcal{Z}})$ when $s^{\prime}<s^{\dagger}$; if
$q^{\dagger}>q$, for $u\in B^{s^{\dagger},q^{\dagger}}({\mathcal{Z}})$ by
Hölder inequality with
$(\frac{q^{\dagger}}{q},\frac{q^{\dagger}}{q^{\dagger}-q})$,

$\displaystyle\|u\|_{s^{\prime},q}^{q}$
$\displaystyle=\sum_{\ell=1}^{\infty}\ell^{\tau_{q}(s^{\prime})q}\|u_{\ell}(\cdot)\|_{q}^{q}\leq\sum_{\ell=1}^{\infty}\ell^{\tau_{q^{\dagger}}(s^{\dagger})q}\|u_{\ell}(\cdot)\|_{q^{\dagger}}^{q}\ell^{(\tau_{q}(s^{\prime})-\tau_{q^{\dagger}}(s^{\dagger}))q}$
(38)
$\displaystyle\leq\|u\|_{s^{\dagger},q^{\dagger}}^{q}\left[\sum_{\ell=1}^{\infty}\ell^{(\tau_{q}(s^{\prime})-\tau_{q^{\dagger}}(s^{\dagger}))qq^{\dagger}/(q^{\dagger}-q)}\right]^{1-\frac{q}{q^{\dagger}}}<\infty$

holds when
$(\tau_{q}(s^{\prime})-\tau_{q^{\dagger}}(s^{\dagger}))qq^{\dagger}/(q^{\dagger}-q)=\frac{(s^{\prime}-s^{\dagger})/d}{\frac{1}{q}-\frac{1}{q^{\dagger}}}-1<-1$,
i.e., $s^{\prime}<s^{\dagger}$; if $q^{\dagger}<q$, for $u\in
B^{s^{\dagger},q^{\dagger},q}({\mathcal{Z}})$, we have

$\displaystyle\|u\|_{s^{\prime},q}^{q}$
$\displaystyle=\sum_{\ell}\ell^{\tau_{q}(s^{\prime})q}\|u_{\ell}(\cdot)\|_{q}^{q}=\sum_{\ell}\ell^{\tau_{q^{\dagger}}(s^{\dagger})q^{\dagger}}\|u_{\ell}(\cdot)\|_{q}^{q^{\dagger}}\ell^{(\tau_{q}(s^{\prime})q-\tau_{q^{\dagger}}(s^{\dagger})q^{\dagger})}\|u_{\ell}(\cdot)\|_{q}^{(q-q^{\dagger})}$
(39)
$\displaystyle\leq\sum_{\ell}\ell^{\tau_{q^{\dagger}}(s^{\dagger})q^{\dagger}}\|u_{\ell}(\cdot)\|_{q}^{q^{\dagger}}\ell^{(\tau_{q}(s^{\prime})q-\tau_{q^{\dagger}}(s^{\dagger})q^{\dagger})}\|u\|_{s^{\dagger},q^{\dagger},q}^{(q-q^{\dagger})}\ell^{-\tau_{q^{\dagger}}(s^{\dagger})(q-q^{\dagger})}$
$\displaystyle=\|u\|_{s^{\dagger},q^{\dagger},q}^{(q-q^{\dagger})}\sum_{\ell}\ell^{\tau_{q^{\dagger}}(s^{\dagger})q^{\dagger}}\|u_{\ell}(\cdot)\|_{q}^{q^{\dagger}}\ell^{(\tau_{q}(s^{\prime})-\tau_{q^{\dagger}}(s^{\dagger}))q}\lesssim\|u\|_{s^{\dagger},q^{\dagger},q}^{q}$

holds when
$(\tau_{q}(s^{\prime})-\tau_{q^{\dagger}}(s^{\dagger}))q=\left(\frac{s^{\prime}-s^{\dagger}}{d}-\frac{1}{q}+\frac{1}{q^{\dagger}}\right)q<0$,
i.e., $s^{\prime}<s^{\dagger}+\frac{d}{q}-\frac{d}{q^{\dagger}}$. ∎

See 4.3

###### Proof of Theorem 4.3.

The series expansion (22) is the result of applying Theorem 3.2 to each
$\xi_{\ell}(\cdot)\overset{iid}{\sim}\mathrm{q}\\!-\\!\mathcal{EP}(0,{\mathcal{C}})$
and the convergence is in $L^{2}_{\mathbb{P}}(\Omega,B^{s,2}({\mathcal{Z}}))$.
We then directly compute the spatiotemporal covariance

$\displaystyle~{}~{}~{}~{}\mathrm{Cov}(u({\bf z}),u({\bf z}^{\prime}))$
$\displaystyle=\mathrm{E}(u({\bf z})u({\bf
z}^{\prime}))=\mathrm{E}\left[\sum_{\ell=1}^{\infty}\gamma_{\ell}\xi_{\ell}(t)\phi_{\ell}({\bf
x})\sum_{\ell^{\prime}=1}^{\infty}\gamma_{\ell^{\prime}}\xi_{\ell^{\prime}}(t^{\prime})\phi_{\ell^{\prime}}({\bf
x}^{\prime})\right]$
$\displaystyle=\sum_{\ell,\ell^{\prime}=1}^{\infty}\gamma_{\ell}\gamma_{\ell^{\prime}}\phi_{\ell}({\bf
x})\phi_{\ell^{\prime}}({\bf
x}^{\prime})\mathrm{E}[\xi_{\ell}(t)\xi_{\ell^{\prime}}(t^{\prime})]=\sum_{\ell=1}^{\infty}\gamma_{\ell}^{2}\phi_{\ell}({\bf
x})\phi_{\ell}({\bf
x}^{\prime})\mathrm{E}[\xi_{\ell}(t)\xi_{\ell}(t^{\prime})]$
$\displaystyle=\sum_{\ell=1}^{\infty}\gamma_{\ell}^{2}\phi_{\ell}({\bf
x})\phi_{\ell}({\bf x}^{\prime}){\mathcal{C}}(t,t^{\prime})$

where we use the iid assumption of $\xi_{\ell}(\cdot)$ so that
$\mathrm{E}[\xi_{\ell}(t)\xi_{\ell^{\prime}}(t^{\prime})]=\mathrm{E}[\xi_{\ell}(t)\xi_{\ell}(t^{\prime})]\delta_{\ell\ell^{\prime}}$.
∎

See 4.4

###### Proof of Theorem 4.4.

Based on the series representation of $u(\cdot)$ in (21), we have by Jensen’s
inequality

$\displaystyle\mathbb{E}[\|u({\bf z})-u({\bf
z}^{\prime})\|_{\tau_{q}(s^{\prime}),q}^{q}]$
$\displaystyle\leq\sum_{\ell=1}^{\infty}\mathbb{E}[\ell^{(\tau_{q}(s^{\prime})-\tau_{q}(s))q}|\phi_{\ell}({\bf
x})\xi_{\ell}(t)-\phi_{\ell}({\bf x}^{\prime})\xi_{\ell}(t^{\prime})|^{q}]$
$\displaystyle\lesssim\sum_{\ell=1}^{\infty}\ell^{\frac{s^{\prime}-s}{d}q}\mathbb{E}[|\Delta\phi_{\ell}|^{q}|\xi_{\ell}(t)|^{q}+|\phi_{\ell}({\bf
x}^{\prime})|^{q}|\Delta\xi_{\ell}|^{q}]\quad(convexity\;of\;|\cdot|^{q})$
$\displaystyle\lesssim\sum_{\ell=1}^{\infty}\ell^{\frac{s^{\prime}-s}{d}q}[L(\phi_{\ell})|{\bf
x}-{\bf
x}^{\prime}|^{q\alpha+d}+\|\phi_{\ell}(\cdot)\|_{\infty}^{q}\mathbb{E}|\Delta\xi_{\ell}|^{q}]$

where $\Delta\phi=\phi({\bf x})-\phi({\bf x}^{\prime})$ and
$\Delta\xi=\xi(t)-\xi(t^{\prime})$. Now based on the series representation
$\xi(\cdot)$ in (12), by Jensen’s inequality and the proof of Theorem 3.3 we
have further for $\forall\ell\in{\mathbb{N}}$ (due to the iid assumption of
$\xi_{\ell}(\cdot)$ in Assumption 2)

$\mathbb{E}[|\Delta\xi_{\ell}|^{q}]\leq\sum_{\ell^{\prime}=1}^{\infty}\mathbb{E}[|\xi_{\ell^{\prime}}|^{q}|\psi_{\ell^{\prime}}(t)-\psi_{\ell^{\prime}}(t^{\prime})|^{q}]\leq\sum_{\ell^{\prime}=1}^{\infty}\lambda_{\ell^{\prime}}^{\frac{q}{2}}L(\psi_{\ell^{\prime}})|t-t^{\prime}|^{q\alpha+d}$

Therefore by Assumption 3 we have

$\displaystyle\mathbb{E}[\|u({\bf z})-u({\bf
z}^{\prime})\|_{\tau_{q}(s^{\prime}),q}^{q}]$ $\displaystyle\lesssim|{\bf
x}-{\bf
x}^{\prime}|^{q\alpha+d}\sum_{\ell=1}^{\infty}\ell^{\frac{s^{\prime}-s}{d}q}L(\phi_{\ell})+|t-t^{\prime}|^{q\alpha+d}\sum_{\ell=1}^{\infty}\ell^{\frac{s^{\prime}-s}{d}q}\|\phi_{\ell}(\cdot)\|_{\infty}^{q}\sum_{\ell^{\prime}=1}^{\infty}\lambda_{\ell^{\prime}}^{\frac{q}{2}}L(\psi_{\ell^{\prime}})\lesssim|{\bf
z}-{\bf z}^{\prime}|^{q\alpha+d}$

The conclusion follows by the Kolmogorov continuity theorem [15]. ∎

###  A.2 Posterior Theorems of Inverse Problems with STBP Priors

See 4.5

###### Proof of Theorem 4.5.

The proof is based on [Theorem 3.2 of 14] and [Theorem 14 15] and we include
it here for the completeness. Define
$\pi_{0}(du,dy)=\Pi(du)\otimes{\mathbb{Q}}_{0}(dy)$ and
$\pi(du,dy)=\Pi(du)\otimes{\mathbb{Q}}_{u}(dy)$. We assume
${\mathbb{Q}}_{0}\ll{\mathbb{Q}}$ and the Radon-Nikodym derivative (2) holds
for $\Pi$-a.s.. Thus, for fixed $u$,
$\Phi(u;\cdot):{\mathbb{Y}}\to{\mathbb{R}}$ is ${\mathbb{Q}}_{0}$-measurable
and $\int_{\mathbb{Y}}\exp(-\Phi(u;y)){\mathbb{Q}}_{0}(dy)=1$. On other hand,
by Assumption 4-(iii), we have
$\Phi(\cdot;y):B^{s^{\prime},q}({\mathcal{Z}})\to{\mathbb{R}}$ is continuous
on $B^{s^{\prime},q}({\mathcal{Z}})$. By Corollary 4.1
$\Pi(B^{s^{\prime},q}({\mathcal{Z}}))=1$. Hence $\Phi(\cdot;y)$ is
$\Pi$-measurable. Therefore, $\Phi$ is $\pi_{0}$-measurable and
$\pi\ll\pi_{0}$ with

$\frac{d\pi}{d\pi_{0}}(u;y)=\exp(-\Phi(u;y)),\quad\int_{B^{s^{\prime},q}({\mathcal{Z}})\times{\mathbb{Y}}}\exp(-\Phi(u;y))\pi_{0}(du,dy)=1$

Then by [Lemma 5.3 of 25] (also [Theorem 13 of 15]), the conditional
distribution
$\Pi^{y}(du):=\pi(du,dy^{\prime}|y^{\prime}=y)\ll\Pi(du)=\pi_{0}(du,dy^{\prime}|y^{\prime}=y)$
due to the definition of $\pi_{0}$. The same Lemma/Theorem implies (24) if the
normalizing constant $Z(y)>0$.

First, by Assumption 4-(i), we have

$Z(y)=\int_{B^{s^{\prime},q}({\mathcal{Z}})}\exp(-\Phi(u;y))\Pi(du)\leq\int_{B^{s^{\prime},q}({\mathcal{Z}})}\exp(-M+\alpha_{1}\|u\|_{s^{\prime},q})\Pi(du)<\infty$

where the boundedness is the result of Theorem 4.2-(ii). Now we show $Z(y)>0$.
By Theorem 4.1, we have $R=\mathbb{E}\|u\|_{s^{\prime},q}<\infty$. Since
$\|u\|_{s^{\prime},q}$ is nonnegative, we have
$\Pi(\|u\|_{s^{\prime},q}<R)>0$. Let $r=\max\\{\|y\|_{\mathbb{Y}},R\\}$. By
Assumption 4-(ii), we have

$Z(y)=\int_{B^{s^{\prime},q}({\mathcal{Z}})}\exp(-\Phi(u;y))\Pi(du)\geq\int_{\|u\|_{s^{\prime},q}<R}\exp(-K)\Pi(du)=\exp(-K)\Pi(\|u\|_{s^{\prime},q}<R)>0$

∎

See 4.6

###### Proof of Theorem 4.6.

The proof is based on [Theorem 3.3 of 14] and [Theorem 16 of 15] and included
for the completeness. As in Theorem 4.5, $Z(y),Z(y^{\prime})\in(0,\infty)$. We
directly compute

$\displaystyle 2d_{H}^{2}(\Pi^{y},\Pi^{y^{\prime}})=$
$\displaystyle\int_{B^{s^{\prime},q}({\mathcal{Z}})}\left[Z(y)^{-\frac{1}{2}}\exp\left(-\frac{1}{2}\Phi(u;y)\right)-Z(y^{\prime})^{-\frac{1}{2}}\exp\left(-\frac{1}{2}\Phi(u;y^{\prime})\right)\right]^{2}\Pi(du)$
$\displaystyle\leq$
$\displaystyle\frac{2}{Z(y)}\int_{B^{s^{\prime},q}({\mathcal{Z}})}\left[\exp\left(-\frac{1}{2}\Phi(u;y)\right)-\exp\left(-\frac{1}{2}\Phi(u;y^{\prime})\right)\right]^{2}\Pi(du)$
$\displaystyle+2|Z(y)^{-\frac{1}{2}}-Z(y^{\prime})^{-\frac{1}{2}}|^{2}Z(y^{\prime})$

By the mean value theorem and Assumptions 4 (i) and (iv), we have

$\displaystyle\int_{B^{s^{\prime},q}({\mathcal{Z}})}\left[\exp\left(-\frac{1}{2}\Phi(u;y)\right)-\exp\left(-\frac{1}{2}\Phi(u;y^{\prime})\right)\right]^{2}\Pi(du)$
$\displaystyle\leq\int_{B^{s^{\prime},q}({\mathcal{Z}})}\frac{1}{4}\exp(\alpha_{1}\|u\|_{s^{\prime},q}-M)\exp(2\alpha_{2}\|u\|_{s^{\prime},q}+2C)\|y_{1}-y_{2}\|_{\mathbb{Y}}^{2}\Pi(du)\leq
C\|y_{1}-y_{2}\|_{\mathbb{Y}}^{2}$

where the boundedness is the result of Theorem 4.2-(ii). By the mean value
theorem, Jesen’s inequality and Assumptions 4 (i) and (iv) we have

$\displaystyle
2|Z(y)^{-\frac{1}{2}}-Z(y^{\prime})^{-\frac{1}{2}}|^{2}Z(y^{\prime})\leq
C|Z(y)-Z(y^{\prime})|^{2}$ $\displaystyle\leq
C\left[\int_{B^{s^{\prime},q}({\mathcal{Z}})}|\exp(-\Phi(u;y))-\exp(-\Phi(u;y^{\prime}))|\Pi(du)\right]^{2}$
$\displaystyle\leq
C\int_{B^{s^{\prime},q}({\mathcal{Z}})}|\exp(-\Phi(u;y))-\exp(-\Phi(u;y^{\prime}))|^{2}\Pi(du)\leq
C\|y_{1}-y_{2}\|_{\mathbb{Y}}^{2}$

where we used the above result. ∎

See 4.7

Before we proceed with the proof, we need the following lemma to bound the
Hellinger distance, Kullback-Leibler (K-L) divergence and K-L variation.

######  Lemma A.1.

Suppose the inverse model (1) has the potential function $\Phi$ (2) satisfy
Assumption 4-(iii). Then we have

  * •

$d_{H}(p_{u},p_{u^{\prime}})\lesssim\|u-u^{\prime}\|_{s^{\prime},q}$

  * •

$K(p_{u},p_{u^{\prime}})\lesssim\|u-u^{\prime}\|_{s^{\prime},q}$

  * •

$V(p_{u},p_{u^{\prime}})\lesssim\|u-u^{\prime}\|^{2}_{s^{\prime},q}$

###### Proof.

First, we consider K-L divergence:

$K(p_{u},p_{u^{\prime}})=\int\log\frac{p_{u}}{p_{u^{\prime}}}p_{u}d\mu=\int(\Phi(u^{\prime};y)-\Phi(u;y))p_{u}d\mu(y)\leq
L\|u-u^{\prime}\|_{s^{\prime},q}$

by Assumption 4-(iii). Similarly, we have for K-L variation:

$V(p_{u},p_{u^{\prime}})=\int\left(\log\frac{p_{u}}{p_{u^{\prime}}}\right)^{2}p_{u}d\mu=\int|\Phi(u^{\prime};y)-\Phi(u;y)|^{2}p_{u}d\mu(y)\leq
L^{2}\|u-u^{\prime}\|_{s^{\prime},q}^{2}$

Lastly, we bound the Hellinger distance:

$\displaystyle 2d^{2}_{H}(p_{u},p_{u^{\prime}})$
$\displaystyle=\int(\sqrt{p_{u}}-\sqrt{p_{u^{\prime}}})^{2}d\mu=\int\left[1-\exp\left(\frac{1}{2}\Phi(u;y)-\frac{1}{2}\Phi(u^{\prime};y)\right)\right]^{2}p_{u}d\mu(y)$
$\displaystyle\leq\int\frac{C}{4}|\Phi(u^{\prime};y)-\Phi(u;y)|^{2}p_{u}d\mu(y)\leq\frac{CL^{2}}{4}\|u-u^{\prime}\|_{s^{\prime},q}^{2}$

where the inequality holds for $\|u-u^{\prime}\|_{s^{\prime},q}^{2}$ small
enough. ∎

###### Proof of Theorem 4.7.

Based on [Theorem 1 of 21], it suffices to verify the following two conditions
(the entropy condition (2.4), and the prior mass condition (2.5)) for some
universal constants $\eta,K>0$ and sufficiently large $k\in{\mathbb{N}}$,

$\displaystyle\sup_{\varepsilon>\varepsilon_{n}}\log
N(\eta\varepsilon/2,\\{u\in\Theta_{n}:d_{n,H}(u,u^{\dagger})<\varepsilon\\},d_{n,H})$
$\displaystyle\leq$ $\displaystyle n\varepsilon_{n}^{2}$ (40)
$\displaystyle\frac{\Pi_{n}(u\in\Theta_{n}:k\varepsilon_{n}<d_{n,H}(u,u^{\dagger})\leq
2k\varepsilon_{n})}{\Pi_{n}(B_{n}(u^{\dagger},\varepsilon_{n}))}$
$\displaystyle\leq$ $\displaystyle e^{Kn\varepsilon_{n}^{2}k^{2}/2}$ (41)

where the left side of (40) is called _Le Cam dimension_ [44], logarithm of
the minimal number of $d_{n,H}$-balls of radius $\xi\varepsilon/2$ needed to
cover a ball of radius $\varepsilon$ around the true value
$u^{\dagger}=\sum_{\ell=1}^{\infty}\gamma_{\ell}\xi_{\ell}^{\dagger}(t)\phi_{\ell}({\bf
x})$;
$B_{n}(u^{\dagger},\varepsilon_{n})=\\{u\in\Theta:\frac{1}{n}K_{j}(u^{\dagger},u)\leq\varepsilon^{2},\frac{1}{n}V_{j}(u^{\dagger},u)\leq\varepsilon^{2}\\}$
with $K_{j}(u^{\dagger},u)=K(P_{u^{\dagger},j},P_{u,j})$ and
$V_{j}(u^{\dagger},u)=V(P_{u^{\dagger},j},P_{u,j})$.

We adopt the argument for infinite sequence of functions in [Theorem 4 of 41].
For each $1\leq\ell\leq n$, $\xi_{\ell}(\cdot)\in
B^{s^{\prime},q}({\mathcal{T}})\subset L^{q}({\mathcal{T}})$ satisfy
conditions for Theorem 3.4. Therefore, there exists $B_{n,\ell}\subset
B^{s^{\prime},q}({\mathcal{T}})$ such that (15)-(17) hold for each $\ell$ with
$B_{n}$ replaced by $B_{n,\ell}$, and $\varepsilon_{n}$ replaced by
$\varepsilon_{n,\ell}=c2^{-\ell}\ell^{-\frac{s^{\prime}-s}{d}}\varepsilon_{n}^{2}$
for some constant $c>0$. Note, for given spatial basis $\\{\phi_{\ell}({\bf
x})\\}_{\ell=1}^{\infty}$ and $\gamma=\\{\gamma_{\ell}\\}_{\ell=1}^{\infty}$
in (21), $u\in\Theta=B^{s^{\prime},q}({\mathcal{Z}})$ can be identified with
$\xi_{\mathcal{T}}=\\{\xi_{\ell}(\cdot)\\}_{\ell=1}^{\infty}\in\ell^{q,\tau_{q}(s^{\prime})}(L^{q}({\mathcal{T}}))$
through $f_{\gamma}$ in Definition 4, i.e.
$\Theta\cong\ell^{q,\tau_{q}(s^{\prime})}(L^{q}({\mathcal{T}}))$. Now we set

$\Theta_{n}=\\{u=f_{\gamma}(\xi_{\mathcal{T}})\in\Theta|\xi_{\ell}(\cdot)\in
B_{n,\ell},\;for\;\ell=1,\cdots,n\\}\subset\Theta$

For $\forall u,u^{\prime}\in\Theta_{n}$ such that
$\|\xi_{\ell}(\cdot)-\xi^{\prime}_{\ell}(\cdot)\|_{s^{\prime},q}\leq\varepsilon_{n,\ell}$
for $\ell=1,\cdots,n$, we have
$d_{n,H}(u,u^{\prime})\lesssim\|u-u^{\prime}\|_{s^{\prime},q}\leq\varepsilon_{n}$
by Lemma A.1 and (15). Therefore
$N(\varepsilon_{n},\Theta_{n},d_{n,H})=\max_{1\leq\ell\leq
n}N(4\varepsilon_{n,\ell},B_{n,\ell},\|\cdot\|_{s^{\prime},q})$ and we have
the following global entropy bound holds

$\log N(\varepsilon_{n},\Theta_{n},d_{n,H})\leq
Rn(c2^{-\ell}\ell^{-\frac{s^{\prime}-s}{d}}\varepsilon_{n}^{2})^{2}\leq
n\varepsilon_{n}^{2}$

for some $1\leq\ell\leq n$ and $c>0$, which is stronger than the local entropy
condition (40).

Now by Lemma A.1 and (17), we have

$\displaystyle\Pi_{n}(B_{n}(u^{\dagger},\varepsilon_{n}))$
$\displaystyle\geq\Pi_{n}(\|u^{\dagger}-u\|_{s^{\prime},q}\leq\varepsilon_{n}^{2},\|u^{\dagger}-u\|_{s^{\prime},q}^{2}\leq\varepsilon_{n}^{2})$
$\displaystyle=\Pi_{n}(\|u^{\dagger}-u\|_{s^{\prime},q}^{q}\leq\varepsilon_{n}^{2q})\geq\exp\left\\{\sum_{\ell=1}^{n}\log\mu(\|\xi_{\ell}(\cdot)-\xi^{\prime}_{\ell}(\cdot)\|_{s^{\prime},q}<2\varepsilon_{n,\ell})\right\\}$
$\displaystyle\geq
e^{-n\sum_{\ell=1}^{n}\varepsilon_{n,\ell}^{2}}=e^{-Knk^{2}\varepsilon_{n}^{4}/2},\quad
withK=2,\;k^{2}=\sum_{\ell=1}^{n}2^{-2\ell}\ell^{-2\frac{s^{\prime}-s}{d}}$

Then the prior mass condition (41) is satisfied because the numerator is
bounded by 1. The proof is hence completed. ∎

See 4.8

The following lemma studies the small ball probability in the concentration
function (25).

######  Lemma A.2 (Small ball probability).

Let $\Pi$ be an $\mathcal{STBP}({\mathcal{C}},B^{s,q}({\mathcal{Z}}))$ prior
on $B^{s^{\prime},q}({\mathcal{Z}})$ with $s^{\prime}<s-\frac{d}{q}$. Then as