# All-at-once formulation meets the Bayesian approach: A study of two
prototypical linear inverse problems

Anna Schlintl and Barbara Kaltenbacher
Alpen-Adria-Universität Klagenfurt, Austria
<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

In this work, the Bayesian approach to inverse problems is formulated in an
all-at-once setting. The advantages of the all-at-once formulation are known
to include the avoidance of a parameter-to-state map as well as numerical
improvements, especially when considering nonlinear problems. In the Bayesian
approach, prior knowledge is taken into account with the help of a prior
distribution. In addition, the error in the observation equation is formulated
by means of a distribution. This method naturally results in a whole posterior
distribution for the unknown target, not just point estimates. This allows for
further statistical analysis including the computation of credible intervals.
We combine the Bayesian setting with the all-at-once formulation, resulting in
a novel approach for investigating inverse problems. With this combination we
are able to chose a prior not only for the parameter, but also for the state
variable, which directly influences the parameter. Furthermore, errors not
only in the observation equation, but additionally, in the model can be taken
into account. We analyze this approach with the help of two linear standard
examples, namely the inverse source problem for the Poisson equation and the
backwards heat equation, i.e. a stationary and a time dependent problem.
Appropriate function spaces and derivation of adjoint operators are
investigated. To assess the degree of ill-posedness, we analyze the singular
values of the corresponding all-at-once forward operators. Finally, joint
priors are designed and numerically tested.

Keywords: inverse problems, all-at-once formulations, Bayesian inverse
problems

##  1 Introduction.

In a general setting, the model to work with in the Bayesian approach is given
by

$y^{\delta}=Gx+\delta\eta,$ (1)

where $G$ is the (linear) forward operator, mapping between Hilbert spaces $X$
and $Y$, $\delta\geq 0$ describes the noise level, $x$ is the unknown target
and $y^{\delta}$ denotes the data, when $\eta$ is a random variable describing
the noise. We restrict our considerations to Gaussian noise, i.e.
$\eta\sim\mathcal{N}(0,\Sigma)$. This equation can be transformed such that it
results in a problem under Gaussian white noise given by

$z^{\delta}=\Sigma^{-1}Gx+\delta\Sigma^{-1}\eta.$ (2)

The prior distribution for $x$ is chosen to be normal, i.e.
$\eta\sim\mathcal{N}(0,\frac{\delta^{2}}{\alpha}C_{0})$, with the noise level
$\delta$ and a scaling parameter $\alpha$. In that setting, much is known,
especially that the posterior is also normally distributed with mean and
covariance given by

$\displaystyle x_{\alpha}^{\delta}$ $\displaystyle=C_{0}^{1/2}(\alpha
I+H)^{-1}B^{*}z^{\delta},$ (3) $\displaystyle C_{\alpha}^{\delta}$
$\displaystyle=\delta^{2}C_{0}^{1/2}(\alpha I+H)^{-1}C_{0}^{1/2},$ (4)

where $B=\Sigma^{-1/2}GC_{0}^{1/2}$ and $H=BB^{*}$, see, e.g., ([1, 2, 3, 4,
5, 6, 7, 8, 9, 10]. Here and below the superscript ∗ denotes the Hilbert space
adjoint.

In equation (1) the inverse problem is given in its reduced form, i.e. the
model

$\mathbb{M}(u,x)=0$ (5)

and the observation equation

$\mathbb{O}(u)=y$ (6)

with operators $\mathbb{M}:U\times X\to W^{\prime}\times Y$, $\mathbb{O}:U\to
Y$ mapping between Hilbert spaces $U,X,W^{\prime},Y$, are combined through the
parameter-to-state map $S$, which maps the parameter $x$ to the state $u$ and
is determined by the identity

$\mathbb{M}(S(x),x)=0.$

The problem then transforms to

$G(x):=\mathbb{O}(S(x))=y.$

Contrary to the reduced approch, the all-at-once formulation combines the
model (5) and the observation equation (6) in one system given by

$\mathcal{G}(\mathbf{x})=\mathcal{G}(u,x)=\left(\begin{array}[]{c}\mathbb{M}(u,x)\\\
\mathbb{O}(u)\end{array}\right)=\left(\begin{array}[]{c}0\\\
y\end{array}\right)=\mathbf{y}.$ (7)

This approach has recently gained attention, see, e.g., ([11, 12, 13, 14, 15,
16, 17]). While these paper remain in a purely deterministic setting, it is
the aim of this paper to apply the Bayesian approach to the all-at once
formulation. Incorporating a priori information not only on the parameter but
also on the state is expected to potentially improve reconstructions.
Moreover, considering possible perturbations in both model and observations,
i.e., assuming
$\mathbf{y}^{\delta}=\mathbf{y}+\delta\left(\begin{array}[]{c}\eta_{1}\\\
\eta_{2}\end{array}\right)$, allows to take into account not only noise in the
data but also uncertainty in the model, which is relevant in many
applications. Combinations to a posterior distribution obtained by Bayes’
formula allows to quantify the resulting uncertainty in both state and
parameter.

###  1.1 Examples.

To illustrate the all-at-once setting for inverse problems, two examples of
linear inverse problems for PDEs are given in this section. These examples
will serve as prototypical showcases in the analytic and numerical
considerations of the following sections. Here we will only sketch these
examples; the function spaces and operator definitions will be made clear in
the next section.

To avoid confusion with the space variable, we will rename the parameter
$\theta$ here.

###### Example 1

The first example is a linear inverse source problem where we aim to recover
the source $u$ from the model

$\begin{array}[]{rcll}-\Delta u&=&f+\theta&\quad\mbox{ in }\Omega\\\
u&=&0&\quad\mbox{ on }\partial\Omega\end{array}$ (8)

with noisy observations $y^{\delta}$ given through the equation

$y^{\delta}=\mathbb{O}u+\delta\eta\quad\in\Omega.$ (9)

The problem can be reformulated in an all-at-once fashion resulting in

$\displaystyle\mathbf{y}^{\delta}=\mathcal{G}\mathbf{x}+\delta\boldsymbol{\eta}\Longleftrightarrow\left(\begin{array}[]{c}y^{\delta}_{1}(x)\\\
y^{\delta}_{2}(x)\end{array}\right)=\left(\begin{array}[]{cc}-\Delta&-I\\\
\mathbb{O}&0\end{array}\right)\left(\begin{array}[]{c}u(x)\\\
\theta(x)\end{array}\right)+\delta\left(\begin{array}[]{c}\eta_{1}\\\
\eta_{2}\end{array}\right),$ (18)

where $\mathbf{y}$ is now a vector consisting of the known inhomogeneity term
$f$ and of the observation $y$, $\mathcal{G}$ is a block matrix operator, the
target $\mathbf{x}$ is a vector of the state variable $u$ and the source
$\theta$. Additionally, the noise consists of two components, where the
possible perturbation in the system is modeled through $\eta_{1}$ and the
error in the observation equation by $\eta_{2}$.

Inverse source problems with possibly more complicated elliptic operators in
place of $-\Delta$ as well as restricted measurements $\mathbb{O}u$, arise in
numerous applications. For example, with the Helmholtz equation in place of
the Poisson equation as a model, the problem corresponds to the frequency
domain formulation of reconstructing sound sources from measurements taken by
a microphone array, see, e.g. [18, 19].

###### Example 2

The second example of interest in this paper is the backward heat equation.
The ambition is to recover the initial temperature when measurements are only
available at some point later in time. This means, contrary to the first
example this problem is time dependent. We consider the model as

$\begin{array}[]{rcll}\partial_{t}u-\Delta u&=&f&\mbox{in
}(0,T)\times\Omega\\\ u&=&0&\mbox{on }(0,T)\times\partial\Omega\\\
u(0,x)&=&\theta(x)&x\in\Omega.\end{array}$

with given observation equation

$y^{\delta}(x)=u(T,x)+\delta\eta\quad\mbox{ in }\Omega.$

To formulate the model and the observation equation in one system, some
rewriting has to be carried out first. To symbolize the observation of
$v(t,x)$ in point $T$, we use the Dirac operator defined by
$\delta_{T}v(t,x):=v(T,x)$, which maps the state to its value at the final
time $T$. Further, we want to incorporate the initial condition into the state
space, which we do by making the ansatz $u(t,x):=\hat{u}(t,x)+\theta(x)$ for
some $\hat{u}$ contained in a linear space of functions vanishing at initial
time. This (after skipping the hat) leads to the all-at-once formulation

$\displaystyle\mathbf{y}^{\delta}=\mathcal{G}\mathbf{x}+\delta\boldsymbol{\eta}\Longleftrightarrow\left(\begin{array}[]{c}y^{\delta}_{1}\\\
y^{\delta}_{2}\end{array}\right)=\left(\begin{array}[]{cc}\partial_{t}-\Delta&-\underline{I}\Delta\\\
\delta_{T}&\underline{I}\end{array}\right)\left(\begin{array}[]{c}u\\\
\theta\end{array}\right)+\delta\left(\begin{array}[]{c}\eta_{1}\\\
\eta_{2}\end{array}\right).$ (27)

Here $y_{1}(t,x)=f(t,x)$, $y_{2}(x)=y(x)$ and the operator $\underline{I}$
maps an only space-dependent function to a formally space and time dependent
function by assigning $v$ to the time-constant function $t\mapsto v$.

A classical application that can be modeled by backwards diffusion is
deconvolution of images. Further related application examples are the
identification of airborne contaminants [20] and imaging with acoustic or
elastic waves in the presence of strong attenuation, arising, e.g., in
photoacoustic tomography [21].

The remainder of this paper is organized as follows. In Section 2 we provide
more details on appropriate function space settings for the two prototypical
examples above. In particular, since the adjoint of $\mathcal{G}$ is required
for the computation of the posterior mean and covariance (3), (4), but also
for many other reconstruction methods, we will provide details on this. For
both of these examples, the degree of ill-posedness, i.e., the decay reate of
the eigenvalues of $G^{*}G$ is well known in the reduced setting and the
question arises whether this behaviour may change in the all-at-once
formulation. We therefore investigate the eigenvalues of the operators
$\mathcal{G}^{*}\mathcal{G}$ both analytically and numerically in Section 3.
Section 4 is devoted to a convergence analysis of the Bayesian all-at-once
approach, where we can heavily rely on existing literature, especially the
results from [9, 10], that largely carry over to the all-at-once formulation;
however, certain conditions need a different interpretation than in the
reduced setting, which we do by means of the prototypical examples above. The
important question of how to choose priors not only for the parameter but also
for the state is discussed in Section 5. Finally, in Section 6, we provide
some numerical results.

##  2 Function space setting and computation of adjoints.

For the following theoretical and numerical considerations, the two examples
given in Section 1.1 are stated more precisely in terms of functions spaces
and their adjoint operators are computed in this section. To this end, we will
restrict ourselves to full observations $(\mathbb{O}u)(x)=u(x)$, $x\in\Omega$,
and vanishing inhomogeneity $f$, and employ the following notations. The
superscript ∗ denotes the Hilbert space adjoint $B^{*}:Z\to A$ of a linear
operator $B:A\to Z$, whereas $B^{\star}:Z^{\prime}\to A^{\prime}$ denotes its
Banach space adjoint, mapping beween the dual spaces $A^{\prime}$,
$Z^{\prime}$. The Laplacian $-\Delta$ with homogeneous Dirichlet boundary
conditions will be denoted by $\mathcal{A}$ both when acting from
$H_{0}^{1}(\Omega)$ to its dual $H^{-1}(\Omega)$ and when acting from
$H_{0}^{1}(\Omega)\cap H^{2}(\Omega)$ to $L^{2}(\Omega)$. In the context of
the time dependent model of the backwards diffusion problem, we will denote by
$C(0,T;H)$, $L^{2}(0,T;H)$, $H^{1}(0,T;H)$ or simply shorthand $C(H)$,
$L^{2}(H)$, $H^{1}(H)$, the Bochner spaces of time dependent functions with
values in the space $H$ (which will typically be a space of $x$ dependent
functions). Moreover, to map elements of such a space $H$ into formally time
dependent functions, we will use the operator $\underline{I}$ defined by

$\underline{I}:H\to L^{2}(0,T;H)\,\quad v\mapsto(t\mapsto v)$

(which can as well be considered as an operator mapping into $C(0,T;H)$, or
$H^{1}(0,T;H)$ or even $C^{\infty}(0,T;H)$), whose $H$-$L^{2}(0,T;H)$ Hilbert
space adjoint is the averaging operator defined by

$\underline{I}^{*}:L^{2}(0,T;H)\to H\,\quad z\mapsto\int_{0}^{T}z(t)\,dt\,.$

Finally, by $(e^{-\mathcal{A}t})_{t>0}$ we denote the semigroup associated
with the heat equation, cf., e.g., [22] and Subsection 5.4 below.

###  2.1 Inverse source problem.

The all-at-once formulation of the inverse source problem, as given in (18)
consists of the vector of functions $\mathbf{x}(x)=(u(x),\theta(x))^{T}$ and
the block operator matrix

$\mathcal{G}=\left(\begin{array}[]{cc}\mathcal{A}&-I\\\
I&0\end{array}\right),$ (28)

which acts from the space $U\times X:=H_{0}^{1}(\Omega)\cap
H^{2}(\Omega)\times L^{2}(\Omega)$ to $W^{\prime}\times Y:=L^{2}(\Omega)\times
L^{2}(\Omega)$. Here $I$ denotes both the embedding of $H_{0}^{1}(\Omega)\cap
H^{2}(\Omega)$ into $L^{2}(\Omega)$ and the identity on $L^{2}(\Omega)$. The
scalar product on $U$ for elements $u,v\in U$ is defined as

$\langle
v(x),u(x)\rangle_{U}:=\int_{\Omega}\mathcal{A}v(x)\mathcal{A}u(x)\mbox{ dx}.$
(29)

For functions $\mathbf{x_{1}}(x)=(u_{1}(x),\theta_{1}(x))$ and
$\mathbf{x_{2}}(x)=(u_{2}(x),\theta_{2}(x))\in U\times L^{2}(\Omega)$ it holds
that

$\langle\mathbf{x_{1}},\mathbf{x_{2}}\rangle_{U\times
L^{2}(\Omega)}=\int_{\Omega}\left[\mathcal{A}u_{1}(x)\mathcal{A}u_{2}(x)+\theta_{1}(x)\theta_{2}(x)\right]\mbox{
dx }.$ (30)

Therefore, the adjoint $\mathcal{G}^{*}$ of $\mathcal{G}$ can be computed from

$\displaystyle\langle\mathcal{G}\mathbf{x_{1}},\mathbf{x_{2}}\rangle_{L^{2}\times
L^{2}}$
$\displaystyle=\int_{\Omega}\left[(\mathcal{A}u_{1}-\theta_{1})u_{2}+u_{1}\theta_{2}\right]\mbox{
dx }$ (31)
$\displaystyle=\int_{\Omega}\left[\mathcal{A}u_{1}\mathcal{A}(\mathcal{A}^{-1}u_{2}+\mathcal{A}^{-2}\theta_{2})-\theta_{1}u_{2}\right]\mbox{
dx }=(\mathbf{x}_{1},\mathcal{G}^{*}\mathbf{x}_{2})_{U\times L^{2}},$ (32)

resulting in

$\mathcal{G}^{*}=\left(\begin{array}[]{cc}\mathcal{A}^{-1}&\mathcal{A}^{-2}\\\
-I&0\end{array}\right),\quad\mathcal{G}^{*}\mathcal{G}=\left(\begin{array}[]{cc}I+\mathcal{A}^{-2}&-\mathcal{A}^{-1}\\\
-\mathcal{A}&I\end{array}\right).$ (33)

###  2.2 Backwards heat problem.

The same analysis can be done for the backwards heat problem, although it is a
bit trickier, due to time dependence. The operator of interest in the all-at-
once formulation of the backwards heat equation is given by

$\mathcal{G}=\left(\begin{array}[]{cc}\partial_{t}+\mathcal{A}&\underline{I}\mathcal{A}\\\
\delta_{T}&I\end{array}\right),$ (34)

which is an operator from the space $U_{0}\times X$ to $W^{\prime}\times
Y:=L^{2}(H^{-1}(\Omega)\times L^{2}(\Omega)$ where $X=H_{0}^{1}(\Omega)$,

$U_{0}:=\\{w\in L^{2}(H_{0}^{1}(\Omega))\cap
H^{1}(H^{-1}(\Omega)):w(0,x)=0\\}\,,$

which means, that the initial condition $(u+\underline{I}\theta)(0)=\theta$ is
implicitly enforced through the function space we are using for $u$. The
operator $\underline{I}\mathcal{A}$ maps the static parameter into a time-
dependent function space, i.e. $\theta\mapsto(t\mapsto\mathcal{A}\theta)$, and
$I$ here denotes the embedding of $H_{0}^{1}(\Omega)$ into $L^{2}(\Omega)$.
The scalar product on $U_{0}$ for all $v(t,x),u(t,x)\in U_{0}$ is given by

$\displaystyle\langle v,u\rangle_{U_{0}}$
$\displaystyle=\int_{0}^{T}\int_{\Omega}\left[\nabla v\nabla
u+\nabla\mathcal{A}^{-1}v_{t}\nabla\mathcal{A}^{-1}u_{t}\right]\mbox{dx
dt}+\int_{\Omega}v(T)u(T)\mbox{ dx}$ (35)
$\displaystyle=\int_{0}^{T}\int_{\Omega}\mathcal{A}^{-1/2}(\partial_{t}+\mathcal{A})v\,\mathcal{A}^{-1/2}(\partial_{t}+\mathcal{A})u\mbox{
dx dt}$ (36)
$\displaystyle=\int_{0}^{T}\int_{\Omega}v\left[\mathcal{A}u-\mathcal{A}^{-1}u_{tt}\right]\mbox{
dx dt}+\int_{\Omega}v(T)\left[u(T)+\mathcal{A}^{-1}u_{t}(T)\right]\mbox{ dx}.$
(37)

The other scalar products of interest are those of $L^{2}(H^{-1})$ and of
$H_{0}^{1}(\Omega)$:

$\langle
v,u\rangle_{L^{2}(H^{-1})}=\int_{0}^{T}\int_{\Omega}\nabla\mathcal{A}^{-1}v\cdot\nabla\mathcal{A}^{-1}u\mbox{
dx dt}\,,\quad\langle v,u\rangle_{H_{0}^{1}}=\int_{\Omega}\nabla v\cdot\nabla
u\mbox{ dx}$ (38)

With the help of these scalar products the operator $\mathcal{G}^{*}$ will
form as

$\mathcal{G}^{*}=\left(\begin{array}[]{cc}(\partial_{t}+\mathcal{A})^{*}&\delta_{T}^{*}\\\
(\underline{I}\mathcal{A})^{*}&I^{*}\end{array}\right),$

where for some components of the operator a bit more investigation is needed.
We start with computing the first operator $(\partial_{t}+\mathcal{A})^{*}$
with the help of the scalar product equation

$\langle(\partial_{t}+\mathcal{A})u,v\rangle_{L^{2}(H^{-1})}=\langle
u,(\partial_{t}+\mathcal{A})^{*}v\rangle_{U_{0}},$ (39)

where the left hand side computes as

$\displaystyle\langle(\partial_{t}+\mathcal{A})u,v\rangle_{L^{2}(H^{-1})}$
$\displaystyle=\int_{0}^{T}\int_{\Omega}u(t,x)\left[\mathcal{A}^{-1}v_{t}+v\right](t,x)\mbox{
dx dt}$ $\displaystyle+\int_{\Omega}u(T,x)\mathcal{A}^{-1}v(T,x)\mbox{ dx}$

and the right hand side, with $(\partial_{t}+\mathcal{A})^{*}v=:z$ according
to the identity (35) works out as

$\displaystyle\langle
u,(\partial_{t}+\mathcal{A})^{*}v\rangle_{U_{0}}=\int_{0}^{T}\int_{\Omega}u(t,x)\left[\mathcal{A}z-\mathcal{A}^{-1}z_{tt}\right](t,x)\mbox{
dx dt}$
$\displaystyle+\int_{\Omega}u(T,x)[z(T,x)+\mathcal{A}^{-1}z_{t}(T,x)]\mbox{
dx}.$

Then, equation (39) leads to the system

$\begin{array}[]{rcll}\mathcal{A}^{2}z-z_{tt}&=&\mathcal{A}v-v_{t}&\mbox{in
}(0,T)\\\ \mathcal{A}z(T)+z_{t}(T)&=&v(T)\end{array}$ (40)

By factorizing
$\mathcal{A}^{2}-\partial_{tt}=(\mathcal{A}-\partial_{t})(\mathcal{A}+\partial_{t})$
and from the fact that $z\in U_{0}$ the system in (40) results in

$\begin{array}[]{rcll}(\partial_{t}+\mathcal{A})z&=&v&\mbox{in }(0,T)\\\
z(0)&=&0\,.\end{array}$

Therefore the adjoint operator $(\partial_{t}+\mathcal{A})^{*}$ applied to $v$
can be written as the solution of the heat equation using the variation of
constants formula for the heat semigroup

$((\partial_{t}+\mathcal{A})^{*}v)(t,x)=\int_{0}^{t}e^{-\mathcal{A}(t-s)}v(s,x)\mbox{
ds}.$

The next operator computed is $\delta_{T}^{*}$. We proceed as before with the
scalar product equation

$(\delta_{T}u,v)_{L^{2}(\Omega)}=(u,\delta_{T}^{*}v)_{U_{0}},$

which, with rewriting $z=\delta_{T}^{*}v(x)$, due to (35) is equivalent to

$\displaystyle\int_{\Omega}u(T,x)v(x)\mbox{
dx}=\int_{0}^{T}\int_{\Omega}u(t,x)[\mathcal{A}z-\mathcal{A}^{-1}z_{tt}](t,x)\mbox{
dx dt}$ (41)
$\displaystyle+\int_{\Omega}u(T,x)[z(T,x)+\mathcal{A}^{-1}z_{t}(T,x)]\mbox{
dx}.$ (42)

Then, (41) leads to the system

$\begin{array}[]{rcll}(\partial_{t}-\mathcal{A})(\partial_{t}+\mathcal{A})z&=&0&\mbox{in
}(0,T)\\\ (\partial_{t}+\mathcal{A})z(T)&=&\mathcal{A}v\,.\end{array}$

Again using the heat semigroup, the solution can be given as

$(\delta_{T}^{*}v)(t,x)=\mathcal{A}\int_{0}^{t}e^{-\mathcal{A}(T+t-2s)}v(x)\mbox{
ds}={\textstyle\frac{1}{2}}[e^{-\mathcal{A}(T-t)}-e^{-\mathcal{A}(T+t)}]v(x).$

The other two adjoint operators compute as

$((\mathcal{A}\underline{I})^{*}u)(x)=\mathcal{A}^{-1}\int_{0}^{T}u(s,x)\mbox{
ds},\quad(I^{*}v)(x)=(\mathcal{A}^{-1}v)(x).$

Altogether, the adjoint operator $\mathcal{G}^{*}$ is given by

$\mathcal{G}^{*}=\left(\begin{array}[]{cc}\int_{0}^{\cdot}e^{-\mathcal{A}(\cdot-s)}.(s)\mbox{
ds}&{\textstyle\frac{1}{2}}[e^{-\mathcal{A}(T-\cdot)}-e^{-\mathcal{A}(T+\cdot)}]\\\
\mathcal{A}^{-1}\underline{I}^{*}&\mathcal{A}^{-1}\end{array}\right),$

and $\mathcal{G}^{*}\mathcal{G}$ is given by

$\mathcal{G}^{*}\mathcal{G}=\left(\begin{array}[]{cc}I+{\textstyle\frac{1}{2}}[e^{-\mathcal{A}(T-\cdot)}-e^{-\mathcal{A}(T+\cdot)}]\delta_{T}&\underline{I}-e^{-\mathcal{A}\cdot}+{\textstyle\frac{1}{2}}[e^{-\mathcal{A}(T-\cdot)}-e^{-\mathcal{A}(T+\cdot)}]\\\
\underline{I}^{*}+2\mathcal{A}^{-1}\delta_{T}&TI+\mathcal{A}^{-1}\end{array}\right).$

##  3 Analysis of the eigenvalues.

In this section the eigenvalues of the all-at-once operators are analyzed in
terms of the two prototypical examples. For the analysis, the adjoint
operators from Section 2 will be used. Especially the operator
$\mathcal{G}^{*}\mathcal{G}$ is of interest. As one might have already
noticed, these operators do not look symmetric. However, they are indeed
symmetric, but with respect to the specific inner product on $U\times
L^{2}(\Omega)$ and $U_{0}\times H_{0}^{1}(\Omega)$, respectively. Therefore, a
transformation is applied first, to find a representation of the operator in
an $L^{2}$-related inner product. It can be shown that this transformed
operator leads to approximate eigenvalues of the true operator.

###### Lemma 3.1

Let $V$ and $H$ be Hilbert spaces with $C:V\to V$ self-adjoint and compact and
$\mathcal{T}:V\to H$ boundedly invertible with $\mathcal{T}\in L(V,H)$ and
$\mathcal{T}^{-1}\in L(H,V)$. Then the operator
$\tilde{C}:=(\mathcal{T}^{-1})^{*}C\mathcal{T}^{-1}:H\to H$ is self-adjoint
and compact and the eigenvalues $\lambda_{k}$ of $C$ and $\mu_{k}$ of
$\tilde{C}$ decay at the same rate, more precisely it holds

$\frac{1}{\|\mathcal{T}^{-1}||^{2}}\mu_{k}\leq\lambda_{k}\leq\|\mathcal{T}\|^{2}\mu_{k},$

with $\lambda_{1}\geq\lambda_{2}\geq\cdots\geq
0,\mu_{1}\geq\mu_{2}\geq\cdots\geq 0$.

Proof. The proof is based on the Courant-Fischer Theorem, which we quote here
for the convenience of the reader

> Theorem (Courant-Fischer). Let $C:V\to V$ be a selfadjoint and compact
> operator. Then the eigenvalues of $C$ fulfill
>
> $\lambda_{k}=\max\\{\min\\{(Cx,x)_{V}:x\in
> S_{k},\|x\|=1\\}:\dim(S_{k})=k,S_{k}\mbox{ subspace of V}\\}$

Let

$\displaystyle\lambda_{k}$ $\displaystyle=\max\\{\min\\{(Cx,x)_{V}:x\in
S_{k},\|x\|=1\\}:\dim(S_{k})=k,S_{k}\mbox{ subspace of }V\\}=$
$\displaystyle=\max_{\dim(S_{k})=k}\min_{x\in
S_{k},\|x\|=1}(Cx,x)_{V}=\max_{\dim(S_{k})=k}\min_{x\in
S_{k},\|x\|=1}((\mathcal{T}^{-1})^{*}C\mathcal{T}^{-1}\,\mathcal{T}x,\mathcal{T}x)_{V}=$
$\displaystyle=\max_{\dim(S_{k})=k}\min_{x\in
S_{k},\hat{x}=\mathcal{T}x/\|\mathcal{T}x\|,\|x\|=1}(\tilde{C}\hat{x},\hat{x})_{H}\|\mathcal{T}x\|^{2}\quad(\star).$

Due to the fact that

$\hat{x}=\frac{\mathcal{T}x}{\|\mathcal{T}x\|}\in\hat{S}_{k}=\mathcal{T}S_{k},$

and the dimension of $S_{k}$ being $k$, due to regularity of $\mathcal{T}$,
$\hat{S}_{k}$ is of dimension $k$ as well. Therefore, taking the minimum over
a superset and using $\|\mathcal{T}x\|\geq\frac{1}{\|\mathcal{T}^{-1}\|}\|x\|$
results in

$(\star)\geq\frac{1}{\|\mathcal{T}^{-1}\|^{2}}\max_{\dim(\hat{S}_{k})=k}\min_{\hat{x}\in\hat{S}_{k},\|x\|=1}(\tilde{C}x,x)_{H}=\frac{1}{\|\mathcal{T}^{-1}\|^{2}}\mu_{k}.$

Analogously it holds

$\mu_{k}\geq\frac{1}{\|(\mathcal{T}^{-1})^{-1}\|^{2}}\lambda_{k}=\frac{1}{\|\mathcal{T}\|^{2}}\lambda_{k}.$

$\diamondsuit$

With the help of Lemma 3.1 we will transform the operator
$\mathcal{G}^{*}\mathcal{G}$ both for the inverse source and backwards heat
problem and then investigate in the computation of the eigenvalues of the
resulting operators.

###  3.1 Inverse source problem.

The operator $\mathcal{G}^{*}\mathcal{G}$ as stated in (33), will be
transformed according to Lemma 3.1 with the operator $\mathcal{T}:U\times
L^{2}(\Omega)\to L^{2}(\Omega)\times L^{2}(\Omega)$ given by

$\mathcal{T}=\left(\begin{array}[]{cc}\mathcal{A}&0\\\ 0&I\end{array}\right),$

which can be easily seen to be unitary in this setting of spaces, i.e.,
$(\mathcal{T}^{-1})^{*}=\mathcal{T}$, therefore $\mathcal{G}^{*}\mathcal{G}$
transforms to

$\displaystyle\widetilde{\mathcal{G}^{*}\mathcal{G}}=(\mathcal{T}^{-1})^{*}\mathcal{G}^{*}\mathcal{G}\mathcal{T}^{-1}=\left(\begin{array}[]{cc}\mathcal{A}&0\\\
0&I\end{array}\right)\left(\begin{array}[]{cc}I+\mathcal{A}^{-2}&-\mathcal{A}^{-1}\\\
-\mathcal{A}&I\end{array}\right)\left(\begin{array}[]{cc}\mathcal{A}^{-1}&0\\\
0&I\end{array}\right)$
$\displaystyle=\left(\begin{array}[]{cc}I+\mathcal{A}^{-2}&-I\\\
-I&I\end{array}\right),$

which is obviously a compact perturbation of
$\overline{\mathcal{G}^{*}\mathcal{G}}=\left(\begin{array}[]{cc}I&-I\\\
-I&I\end{array}\right)$ whose eigenvalues are $0$ and $2$ with eigenspaces
$E_{0}=\\{(g,g)^{T}\,:\,g\in L^{2}(\Omega)\\}$, $E_{2}=\\{(-g,g)^{T}\,:\,g\in
L^{2}(\Omega)\\}$, that actually span all of $L^{2}(\Omega)$, $E_{0}\oplus
E_{2}=L^{2}(\Omega)$.

####  3.1.1 Analytic computation of the eigenvalues.

The eigenvalues $\lambda$ are computed both analytically and numerically. For
the analytic computation the characteristic eigenvalue equation is considered

$\displaystyle\left(\begin{array}[]{cc}I+\mathcal{A}^{-2}&-I\\\
-I&I\end{array}\right)\left(\begin{array}[]{c}f\\\
g\end{array}\right)=\lambda\left(\begin{array}[]{c}f\\\ g\end{array}\right),$

which leads to the system

$\displaystyle(1-\lambda)f+\mathcal{A}^{-2}f-g=0,\quad(1-\lambda)g-f=0,$

resulting in

$\displaystyle\mathcal{A}^{-2}g=\frac{\lambda(2-\lambda)}{(1-\lambda)}g$
$\displaystyle f=(1-\lambda)g.$

Therefore, let $\mu=\lambda\frac{2-\lambda}{1-\lambda}$ be an eigenvalue of
$\mathcal{A}^{-2}$. Then $\lambda$ solves the quadratic equation

$\lambda^{2}-\lambda(2+\mu)+\mu=0,$

resulting in the two solutions

$\displaystyle\lambda_{1,2}=\frac{2+\mu}{2}\pm\sqrt{\frac{\mu^{2}}{4}+1}=1+\frac{\mu}{2}+\sqrt{1+\frac{\mu^{2}}{4}}\,,\
\frac{\mu}{1+\frac{\mu}{2}+\sqrt{1+\frac{\mu^{2}}{4}}}.$

As $\mu\to 0$, the solutions for $\lambda_{1,2}$ tend to $2$ and $0$,
respectively. Therefore, the eigenvalues of
$\widetilde{\mathcal{G}^{*}\mathcal{G}}$ accumulate at 0 and 2, at the same
(basically linear) speed as $\mu\to 0$.

Thus, besides the singular values tending to zero at a linear rate, known from
the reduced setting and reflecting the mild ill-posedness of the inverse
source problem, we have another sequence tending to a positive value.

####  3.1.2 Numerical computation of the eigenvalues.

The computation is done in python with the help of the finite element
discretization in FEniCS and the eigensolver from SLEPc. The discretzation is
done on a unit square mesh with degree 1 Lagrange elements. To overcome the
computation of $\mathcal{A}^{-2}$ the inverse operator
$\widetilde{\mathcal{G}^{*}\mathcal{G}}^{-1}=\left(\begin{array}[]{cc}\mathcal{A}^{2}&\mathcal{A}^{2}\\\
\mathcal{A}^{2}&\mathcal{A}^{2}+I\end{array}\right)$ is used. The matrices
resulting from the finite element discretization
$u(x)\approx\sum_{i=1}^{n_{el}}u_{i}\phi_{i}(x)$ with the FE basis functions
$\phi_{i}$, and the coefficient vector
$\underline{u}=(u_{1},\ldots,u_{n_{el}})$, according to the identities

$\mathcal{A}u=\lambda u\ \Leftrightarrow\ K\underline{u}=\lambda
M\underline{u}\,,\quad\mathcal{A}^{2}u=\lambda u\ \Leftrightarrow\
KM^{-1}K\underline{u}=\lambda M\underline{u}\,,$

are

$\widetilde{\mathcal{G}^{*}\mathcal{G}}^{-1}_{h}=\left(\begin{array}[]{cc}KM^{-1}K&KM^{-1}K\\\
KM^{-1}K&KM^{-1}K+M\end{array}\right),\quad\mathcal{M}=\left(\begin{array}[]{cc}M&0\\\
0&M\end{array}\right),$

with $M$ and $K$ being the assembled mass and stiffness matrices according to

$\displaystyle M_{i,j}$ $\displaystyle=$
$\displaystyle\int_{\Omega}\phi_{i}\phi_{j},$ (46) $\displaystyle K_{i,j}$
$\displaystyle=$
$\displaystyle\int_{\Omega}\nabla\phi_{i}\cdot\nabla\phi_{j},$ (47)

to state the eigenvalue equation as

$\mathcal{M}\mathbf{x}=\lambda\widetilde{\mathcal{G}^{*}\mathcal{G}}^{-1}_{h}\mathbf{x},$

with $\lambda$ and $\mathbf{x}$ denoting the eigenvalue and eigenvector,
respectively. Sampling of 500 eigenvalues leads then to the following
visualized output as seen in Figure 1(a), the square root scaled plot in
Figure 1(b) suggests quadratic decay, corresponding to the fact that
$\lambda_{2,n}\sim\frac{\mu_{n}}{2}\sim\frac{C}{n^{2}}$ as the eigenvalues of
the Laplacian on the 2-d unit square are given by
$\mu_{j,k}=\frac{j^{2}+k^{2}}{4}$, $j,k\in\mathbb{N}$ and therefore, upon
proper renumbering, decay linearly in $n$.

(a) Eigenvalues

(b) Detailed plot

Figure 1: Numerical results for the eigenvalues of
$\widetilde{\mathcal{G}^{*}\mathcal{G}}$ for the inverse source problem.

###  3.2 Backwards heat equation.

The same analysis is done with the backwards heat equation. Here, the operator
of interest is

$\mathcal{G}^{*}\mathcal{G}=\left(\begin{array}[]{cc}I+\mathcal{A}\int_{0}^{t}e^{-\mathcal{A}(T+t-2s)}\mbox{
.(s)}\delta_{T}&I-e^{-\mathcal{A}t}+\mathcal{A}\int_{0}^{t}e^{-\mathcal{A}(T+t-2s)}\mbox{
.(s)}\\\
\underline{I}^{*}+2\mathcal{A}^{-1}\delta_{T}&TI+\mathcal{A}^{-1}\end{array}\right).$
(48)

The transformation from the space $U_{0}\times H_{0}^{1}(\Omega)$ to the space
$L^{2}(L^{2}(\Omega))\times L^{2}(\Omega)$ is computed with Lemma 3.1 where
the operator $\mathcal{T}$ for the transformation is given by

$\mathcal{T}=\left(\begin{array}[]{cc}\mathcal{A}^{-1/2}(\partial_{t}+\mathcal{A})&0\\\
0&\mathcal{A}^{1/2}\end{array}\right)$

and therefore,

$\widetilde{\mathcal{G}^{*}\mathcal{G}}=(\mathcal{G}\mathcal{T}^{-1})^{*}\mathcal{G}\mathcal{T}^{-1}.$

Here

$\mathcal{T}^{-1}\left(\begin{array}[]{c}f\\\
g\end{array}\right)=\left(\begin{array}[]{c}u\\\ \theta\end{array}\right)\
\Leftrightarrow\
\left\\{\begin{array}[]{c}(\partial_{t}+\mathcal{A})u=\mathcal{A}^{1/2}f\\\
\theta=\mathcal{A}^{-1/2}g\end{array}\right.$

and

$\mathcal{G}\left(\begin{array}[]{c}u\\\
\theta\end{array}\right)=\left(\begin{array}[]{c}a\\\ b\end{array}\right)\
\Leftrightarrow\
\left\\{\begin{array}[]{c}(\partial_{t}+\mathcal{A})u+\mathcal{A}\theta=a\\\
u(T)+\theta=b\end{array}\right.,$

thus

$\left\\{\begin{array}[]{c}\mathcal{A}^{1/2}f(t)+\mathcal{A}^{1/2}g=a(t),\
t\in(0,T)\\\
\int_{0}^{T}e^{-\mathcal{A}(T-t)}\mathcal{A}^{1/2}f(t)\,dt+\mathcal{A}^{-1/2}g=b\end{array}\right.\,,$

i.e.,

$\mathcal{G}\mathcal{T}^{-1}=\left(\begin{array}[]{cc}\mathcal{A}^{1/2}&\mathcal{A}^{1/2}\underline{I}\\\
\int_{0}^{T}e^{-\mathcal{A}(T-s)}\mathcal{A}^{1/2}.(s)\mbox{
ds}&\mathcal{A}^{-1/2}\end{array}\right).$

For computing the adjoint $(\mathcal{G}\mathcal{T}^{-1})^{*}$ we consider the
identity

$\displaystyle\left\langle\mathcal{G}\mathcal{T}^{-1}\left(\begin{array}[]{c}f\\\
g\end{array}\right),\left(\begin{array}[]{c}a\\\
b\end{array}\right)\right\rangle_{L^{2}(H^{-1})\times L^{2}}$ $\displaystyle=$
$\displaystyle\int_{0}^{T}\int_{\Omega}(f(t,x)+g(x))\mathcal{A}^{-1/2}a(t,x)\mbox{
dx}\mbox{ dt}$
$\displaystyle+\int_{\Omega}\left(\int_{0}^{T}e^{-\mathcal{A}(T-s)}\mathcal{A}^{1/2}f(s,x)\mbox{
ds}+\mathcal{A}^{-1/2}g(x)\right)b(x)\mbox{ dx}$ $\displaystyle=$
$\displaystyle\int_{0}^{T}\int_{\Omega}f(t,x)\left(\mathcal{A}^{-1/2}a(t,x)+e^{-\mathcal{A}(T-t)}\mathcal{A}^{1/2}b\right)\mbox{
dx}\mbox{ dt}$
$\displaystyle+\int_{\Omega}g(x)\mathcal{A}^{-1/2}\left(\int_{0}^{T}a(t,x)\mbox{
dt}+b\right)\mbox{ dx},$

which yields

$(\mathcal{G}\mathcal{T}^{-1})^{*}=\left(\begin{array}[]{cc}\mathcal{A}^{-1/2}&e^{-\mathcal{A}(T-\cdot)}\mathcal{A}^{1/2}\\\
\mathcal{A}^{-1/2}\underline{I}^{*}&\mathcal{A}^{-1/2}\end{array}\right).$

Therefore,

$\displaystyle\widetilde{\mathcal{G}^{*}\mathcal{G}}=(\mathcal{G}\mathcal{T}^{-1})^{*}\mathcal{G}\mathcal{T}^{-1}=$
$\displaystyle=\left(\begin{array}[]{cc}I+e^{-\mathcal{A}(T-\cdot)}\mathcal{A}^{1/2}\int_{0}^{T}e^{-\mathcal{A}(T-s)}\mathcal{A}^{1/2}.(s)\mbox{
ds}&\underline{I}+e^{-\mathcal{A}(T-\cdot)}\\\
\underline{I}^{*}+\int_{0}^{T}e^{-\mathcal{A}(T-s)}.(s)\mbox{
ds}&TI+\mathcal{A}^{-1}\end{array}\right).$ (51)

This shows that $\widetilde{\mathcal{G}^{*}\mathcal{G}}$ is a compact
perturbation of the operator
$\overline{\mathcal{G}^{*}\mathcal{G}}=\left(\begin{array}[]{cc}I+D&\underline{I}\\\
\underline{I}^{*}&TI\end{array}\right)$, where
$D=e^{-\mathcal{A}(T-\cdot)}\mathcal{A}^{1/2}\int_{0}^{T}e^{-\mathcal{A}(T-s)}\mathcal{A}^{1/2}.(s)\mbox{
ds}$. To see this, consider a system of eigenvalues $\gamma_{n}$ and
eigenfunctions $\phi_{n}$ of $\mathcal{A}$ and the estimate

$\displaystyle\|\int_{0}^{T}e^{-\mathcal{A}(T-s)}u(s)\mbox{
ds}\|_{\dot{H}^{1}(\Omega)}=\|\mathcal{A}^{1/2}\int_{0}^{T}e^{-\mathcal{A}(T-s)}u(s)\mbox{
ds}\|_{\dot{H}^{1}(\Omega)}$
$\displaystyle=\left(\sum_{n=1}^{\infty}\gamma_{n}\Bigl{(}\int_{0}^{T}e^{-\gamma_{n}(T-s)}\langle
u(s),\phi_{n}\rangle_{L^{2}}\mbox{ ds}\Bigr{)}^{2}\right)^{1/2}$
$\displaystyle\leq\Bigl{(}\sum_{n=1}^{\infty}\int_{0}^{T}\langle
u(s),\phi_{n}\rangle_{L^{2}}^{2}\mbox{
ds}\underbrace{\gamma_{n}\int_{0}^{T}e^{-2\gamma_{n}(T-s)}\mbox{
ds}}_{=\frac{1}{2}(1-e^{-2\gamma_{n}T})}\Bigr{)}^{1/2}\
\leq\frac{1}{\sqrt{2}}\|u\|_{L^{2}(L^{2})}\,,$

which together with compactness of the embedding
$\dot{H}^{1}(\Omega):=\mathcal{D}(\mathcal{A}^{1/2})\to L^{2}(\Omega)$, shows
that the operator
$\int_{0}^{T}e^{-\mathcal{A}(T-s)}.(s):L^{2}(0,T;L^{2}(\Omega))\to
L^{2}(\Omega)$ is compact and so is its adjoint
$e^{-\mathcal{A}(T-\cdot)}:L^{2}(\Omega)\to L^{2}(0,T;L^{2}(\Omega))$. However
$D:L^{2}(0,T;L^{2}(\Omega))\to L^{2}(0,T;L^{2}(\Omega))$ is only bounded, as
the following computation shows

$\displaystyle\|e^{-\mathcal{A}(T-\cdot)}\mathcal{A}^{1/2}\int_{0}^{T}e^{-\mathcal{A}(T-s)}\mathcal{A}^{1/2}u(s)\mbox{
ds}\|_{L^{2}}$
$\displaystyle=\left(\int_{0}^{T}\sum_{n=1}^{\infty}\Bigl{(}\gamma_{n}e^{-\gamma_{n}(T-t)}\int_{0}^{T}e^{-\gamma_{n}(T-s)}\langle
u(s),\phi_{n}\rangle_{L^{2}}\mbox{ ds}\Bigr{)}^{2}\mbox{ dt}\right)^{1/2}$
$\displaystyle\leq\left(\int_{0}^{T}\sum_{n=1}^{\infty}\gamma_{n}^{2}e^{-2\gamma_{n}(T-t)}\int_{0}^{T}e^{-2\gamma_{n}(T-s)}\,ds\
\int_{0}^{T}\langle u(s),\phi_{n}\rangle_{L^{2}}^{2}\mbox{ ds}\mbox{
dt}\right)^{1/2}$ $\displaystyle=\left(\sum_{n=1}^{\infty}\int_{0}^{T}\langle
u(s),\phi_{n}\rangle_{L^{2}}^{2}\mbox{
ds}\Bigl{(}\gamma_{n}e^{-2\gamma_{n}(T-t)}\mbox{ dt}\Bigr{)}^{2}\right)^{1/2}\
\leq\frac{1}{2}\|u\|_{L^{2}(L^{2})}\,.$

####  3.2.1 Analytic computation of the eigenvalues.

With that operator we now again state the eigenvalue equation

$\displaystyle\widetilde{\mathcal{G}^{*}\mathcal{G}}\left(\begin{array}[]{c}f\\\
g\end{array}\right)=\lambda\left(\begin{array}[]{c}f\\\ g\end{array}\right),$
(56)

where we write $f(t)$ and $g$ in terms of their generalized Fourier series
with respect to the eigensystem
$(\mu_{n},\phi_{n})_{n\in\mathbb{N}}\subseteq\mathbb{R}\times L^{2}(\Omega)$
of $\mathcal{A}^{-1}$

$f(t,x)=\sum_{n=1}^{\infty}f_{n}(t)\phi_{n}(x)\,,\quad
g(x)=\sum_{n=1}^{\infty}g_{n}\phi_{n}(x),$

with $f_{n}(t)=\langle f(t),\phi_{n}\rangle_{L^{2}(\Omega)}$, $g_{n}=\langle
g,\phi_{n}\rangle_{L^{2}(\Omega)}$. Then, with (51), upon taking inner
products with $\phi_{n}$, equation (56) leads to the system

$\begin{array}[]{l}(1-\lambda)f_{n}(t)+e^{-\frac{1}{\mu_{n}}(T-t)}\frac{1}{\mu_{n}}\int_{0}^{T}e^{-\frac{1}{\mu_{n}}(T-s)}f_{n}(s)\mbox{
ds}+(1+e^{-\frac{1}{\mu_{n}}(T-t)})g_{n}=0,\\\ \int_{0}^{T}f_{n}(s)\mbox{
ds}+\int_{0}^{T}e^{-\frac{1}{\mu_{n}}(T-s)}f_{n}(s)\mbox{
ds}+(T+\mu_{n}-\lambda)g_{n}=0,\end{array}$ (57)

for all $n\in\mathbb{N}$ and $t\in(0,T)$. This has a nontrivial solution
$(f_{n},g_{n})_{n\in\mathbb{N}}\in L^{2}(0,T,\ell^{2})\times\ell^{2}$, if and
only if there exists $m\in\mathbb{N}$ such that (57) holds for $n=m$ and
$f_{m}\not\equiv 0$ or $g_{m}\not=0$ (to see sufficiency, set all other
components of $(f_{n},g_{n})_{n\in\mathbb{N}}$ to zero and they will trivially
satisfy the linear system (57)). With the particular time dependent functions
$b^{m}_{0}:t\mapsto 1$ and $b^{m}_{1}:t\mapsto e^{-\frac{1}{\mu_{n}}(T-t)}$
playing a role here, this is the case iff either

$\mbox{(a) }\lambda=1\mbox{ and
}\left\\{\begin{array}[]{l}g_{m}b^{m}_{0}(t)+\Bigl{(}{\textstyle\frac{1}{\mu_{m}}}\int_{0}^{T}b^{m}_{1}(s)f_{m}(s)\mbox{
ds}+g_{m}\Bigr{)}b^{m}_{1}(t)=0\,,\ t\in(0,T)\\\
\int_{0}^{T}b^{m}_{0}(s)f_{m}(s)\mbox{
ds}+\int_{0}^{T}b^{m}_{1}(s)f_{m}(s)\mbox{
ds}+(T+\mu_{m}-\lambda)g_{m}=0\end{array}\right.$

or

$\mbox{(b) }\lambda\not=1\mbox{ and
}\left\\{\begin{array}[]{l}f_{m}=a_{m}b^{m}_{0}+c_{m}b^{m}_{1}\mbox{ where}\\\
a_{m}={\textstyle\frac{1}{\lambda-1}}g_{m}\,,\quad
c_{m}={\textstyle\frac{1}{\lambda-1}}\Bigl{(}{\textstyle\frac{1}{\mu_{m}}}\int_{0}^{T}b^{m}_{1}(s)f_{m}(s)\mbox{
ds}+g_{m}\Bigr{)}\\\ \int_{0}^{T}b^{m}_{0}(s)f_{m}(s)\mbox{
ds}+\int_{0}^{T}b^{m}_{1}(s)f_{m}(s)\mbox{
ds}+(T+\mu_{m}-\lambda)g_{m}=0\end{array}\right..$

The first equation of case (a) $\lambda=1$ due to linear independence of the
functions $b^{m}_{0}$, $b^{m}_{1}$ leads to the conditions $g_{m}=0$ and
${\textstyle\frac{1}{\mu_{m}}}\int_{0}^{T}b^{m}_{1}(s)f_{m}(s)\mbox{
ds}+g_{m}=0$, that, combined with the second equation, yield the necessary and
sufficient conditions

$g_{m}=0,\int_{0}^{T}b^{m}_{1}(s)f_{m}(s)\mbox{ ds}=0\mbox{ and
}\int_{0}^{T}b^{m}_{0}(s)f_{m}(s)\mbox{ ds}=0\,.$

This corresponds to a unit eigenvalue with the infinite dimensional eigenspace
$E_{1}$

$\lambda_{0}=1\,,\quad E_{1}=\mbox{span}\\{(x,t)\mapsto\psi(t)\phi_{m}(x)\ :\
\psi\in\\{b^{m}_{0},\,b^{m}_{1}\\}^{\bot}\,,\ m\in\mathbb{N}\\}$ (58)

where ${.}^{\bot}$ denotes the $L^{2}(0,T)$ orthogonal complement.
In case (b) $\lambda_{m}\not=1$ inserting the representation
$f_{m}=a_{m}b^{m}_{0}+c_{m}b^{m}_{1}$ into the two equations involving $f_{m}$
and $g_{m}=(\lambda-1)a_{m}$, we get the following system for $a_{m}$ and
$b_{m}$

$\displaystyle
c_{m}={\textstyle\frac{1}{\lambda-1}}{\textstyle\frac{1}{\mu_{m}}}\int_{0}^{T}b^{m}_{1}(s)(a_{m}b^{m}_{0}(s)+c_{m}b^{m}_{1}(s))\mbox{
ds}+a_{m},$
$\displaystyle\int_{0}^{T}b^{m}_{0}(s)(a_{m}b^{m}_{0}(s)+c_{m}b^{m}_{1}(s))\mbox{
ds}+\int_{0}^{T}b^{m}_{1}(s)(a_{m}b^{m}_{0}(s)+c_{m}b^{m}_{1}(s))\mbox{ ds}$
$\displaystyle\quad+(T+\mu_{m}-\lambda)(\lambda-1)a_{m}=0\,,$

which with the integrals $\int_{0}^{T}b^{m}_{0}(s)^{2}\mbox{ ds}=T$,
$\int_{0}^{T}b^{m}_{1}(s)^{2}\mbox{
ds}=\frac{\mu_{m}}{2}(1-e^{-\frac{1}{\mu_{m}}2T})$,
$\int_{0}^{T}b^{m}_{0}(s)b^{m}_{1}(s)\mbox{
ds}=\mu_{m}(1-e^{-\frac{1}{\mu_{m}}T})$ reads as

$\displaystyle\Bigl{(}(\lambda-1)+1-e^{-\frac{1}{\mu_{m}}T}\Bigr{)}a_{m}+\Bigl{(}-(\lambda-1)+{\textstyle\frac{1}{2}}(1-e^{-\frac{1}{\mu_{m}}2T})\Bigr{)}c_{m}=0,$
$\displaystyle\Bigl{(}T+(T+\mu_{m}-\lambda)(\lambda-1)+\mu_{m}(1-e^{-\frac{1}{\mu_{m}}T})\Bigr{)}a_{m}$
$\displaystyle\quad+\mu_{m}({\textstyle\frac{3}{2}}-e^{-\frac{1}{\mu_{m}}T}-e^{-\frac{1}{\mu_{m}}2T})c_{m}=0\,,$

i.e.,

$\displaystyle(\lambda-e^{-\frac{1}{\mu_{m}}T})a_{m}+({\textstyle\frac{3}{2}}-\lambda-{\textstyle\frac{1}{2}}e^{-\frac{1}{\mu_{m}}2T}))c_{m}=0,$
$\displaystyle\Bigl{(}-\lambda^{2}+\lambda(T+1+\mu_{m})-\mu_{m}e^{-\frac{1}{\mu_{m}}T}\Bigr{)}a_{m}+\mu_{m}({\textstyle\frac{3}{2}}-e^{-\frac{1}{\mu_{m}}T}-e^{-\frac{1}{\mu_{m}}2T})c_{m}=0\,.$

Existence of a nontrivial solution $a_{m}$, $c_{m}$ by setting the determinant
of this system to zero is equivalent to the following cubic equation for
$\lambda$

$\lambda^{3}-\lambda^{2}(T+{\textstyle\frac{5}{2}}+\alpha_{m})+\lambda\left({\textstyle\frac{3}{2}}(T+1)+\beta_{m}\right)=\mu_{m}e^{-\frac{1}{\mu_{m}}2T}\,,$

with $\alpha_{m}=\mu_{m}-{\textstyle\frac{1}{2}}e^{-\frac{1}{\mu_{m}}T}$,
$\beta_{m}=2\mu_{m}e^{-\frac{1}{\mu_{m}}T}-{\textstyle\frac{1}{2}}e^{-\frac{1}{\mu_{m}}2T}(T+1)$
whose solutions $\lambda_{\mu_{m}}^{1}$, $\lambda_{\mu_{m}}^{2}$,
$\lambda_{\mu_{m}}^{3}$ will be the remaining (besides $\lambda_{0}=1$)
eigenvalues of $\widetilde{\mathcal{G}^{*}\mathcal{G}}$. To investigate the
asymptotics of the eigenvalues, in particular of $\lambda^{1}_{\mu_{n}}$, we
consider another sequence of values, namely

$\displaystyle\bar{\lambda}_{1,m}$ $\displaystyle=$ $\displaystyle
0=:\bar{\lambda}_{1},$ $\displaystyle\bar{\lambda}_{2,m}$ $\displaystyle=$
$\displaystyle{\textstyle\frac{T+{\textstyle\frac{5}{2}}+\alpha_{m}}{2}}+{\textstyle\frac{\sqrt{(T-{\textstyle\frac{1}{2}})^{2}+\alpha_{m}^{2}+(2T+5)\alpha_{m}-4\beta_{m}}}{2}}\to\left\\{\begin{array}[]{l}T+1\mbox{
if }T\geq{\textstyle\frac{1}{2}}\\\ {\textstyle\frac{3}{2}}\mbox{ if
}T<{\textstyle\frac{1}{2}}\end{array}\right.=:\bar{\lambda}_{2}\,,$
$\displaystyle\bar{\lambda}_{3,m}$ $\displaystyle=$
$\displaystyle{\textstyle\frac{{\textstyle\frac{3}{2}}(T+1)+\beta_{m}}{\bar{\lambda}_{2,m}}}\to\left\\{\begin{array}[]{l}{\textstyle\frac{3}{2}}\mbox{
if }T\geq{\textstyle\frac{1}{2}}\\\ T+1\mbox{ if
}T<{\textstyle\frac{1}{2}}\end{array}\right.=:\bar{\lambda}_{3}\,.$

(where the limits are to be understood as $m\to\infty$ and therewith
$\mu_{m}\to 0$), so that the cubic equation for $\lambda$ above can be written
as

$\lambda(\lambda-\bar{\lambda}_{2,m})(\lambda-\bar{\lambda}_{3,m})=\mu_{m}e^{-\frac{1}{\mu_{m}}2T}\,.$
(61)

Since $\lambda^{i}_{\mu_{m}}$ are roots of cubic polynomial whose coefficients
converge to those of a cubic polynomial with the three single roots
$\bar{\lambda}_{1}$, $\bar{\lambda}_{2}$, $\bar{\lambda}_{3}$, we also have
the convergence

$\lambda^{i}_{\mu_{m}}\to\bar{\lambda}_{i}\mbox{ as }m\to\infty\,,\
i\in\\{1,2,3\\}\,.$

In particular this means that since
$\bar{\lambda}_{1}=0<\bar{\lambda}_{2},\bar{\lambda}_{3}$ and from (61) we
have

$\lambda^{1}_{\mu_{m}}=\frac{\mu_{m}e^{-\frac{1}{\mu_{m}}2T}}{(\lambda^{1}_{\mu_{m}}-\bar{\lambda}_{2,m})(\lambda^{1}_{\mu_{m}}-\bar{\lambda}_{3,m})}\,,$

where the denominator on the right hand side is positive and bounded away from
zero for $m$ sufficiently large. Hence, there exists $m_{0}\in\mathbb{N}$ such
that for all $m\geq n_{0}$ we have $\lambda^{1}_{\mu_{m}}>0$ (compatibly with
the fact that $\lambda^{1}_{\mu_{m}}$ is an eigenvalue of the positive
semidefinite operator $\widetilde{\mathcal{G}^{*}\mathcal{G}}$) and
$\lambda^{1}_{\mu_{m}}=O(\mu_{m}e^{-\frac{1}{\mu_{m}}2T})$.
Thus, we here have, in addition to $\lambda_{0}=1$, three sequences of
eigenvalues: One tending to zero at an exponential rate (like in the reduced
setting) in according to the severe ill-posedness of the backwards diffusion
problem, and two further sequences accumulating at the positive values
${\textstyle\frac{3}{2}}$ and $T+1$.

####  3.2.2 Numerical computation of the eigenvalues.

The computation is done in python similarly to the computation of the
eigenvalues of the inverse source problem, relying on a finite element
discretization of the Laplacian. To obtain the right matrices for representing
the semigroup expressions appearing in the definition of
$\widetilde{\mathcal{G}^{*}\mathcal{G}}$, we exemplarily consider the term
$w(t)=\int_{0}^{t}e^{-(t-s)\mathcal{A}}f(s)\,ds$ for some $f\in
L^{2}(L^{2}(\Omega))$, and make a semidiscretization in space with a finite
element ansatz $f(t,x)\approx\sum_{i=1}^{n_{el}}f_{i}(t)\phi_{i}(x)$,
$w(t,x)\approx\sum_{i=1}^{n_{el}}w_{i}(t)\phi_{i}(x)$ with basis functions
$\phi_{1},\ldots,\phi_{n}$. Taking into account the fact that $w$ solves

$\begin{array}[]{rcll}\partial_{t}w+\mathcal{A}w&=&f&\mbox{ in
}(0,T)\times\Omega\\\ w&=&0&\mbox{ on }(0,T)\times\partial\Omega\\\
w(0,x)&=&0&x\in\Omega,\end{array}$

i.e., inserting the above ansatz, and testing with FE shape functions
$\phi_{j}$, $j\in\\{1,\ldots,n_{el}\\}$, (a procedure known as Faedo-Galerkin
approximation,) we end up with the system of ODEs

$M\dot{\underline{w}}(t)+K\underline{w}(t)=M\underline{f}(t)\,,t\in(0,T)\,,\quad\underline{w}(0)=0\,,$

where $M$ and $K$ are FE mass and stiffness matrices according to (46), (47),
which, in order to obtain symmetry with respect to the Euclidean inner
product, we write as
$M^{1/2}\dot{\underline{w}}(t)+M^{-1/2}KM^{-1/2}\,M^{1/2}\underline{w}(t)=M^{1/2}\underline{f}(t)$.
Thus the coefficient vector function
$\underline{w}(t)=(w_{1}(t),\ldots,w_{n_{el}}(t))^{T}$ is determined by the
identity

$M^{1/2}\underline{w}(t)=\int_{0}^{t}e^{-(t-s)A_{h}}M^{1/2}\underline{f}(s)\,,ds$

with $A_{h}=M^{-1/2}KM^{-1/2}$.
Additional discretization in time is done by piecewise constant basis
functions on a uniform partition $0=t_{0}<t_{1}<\cdots<t_{N}=T$,
$t_{k}=\frac{T}{N}j=\tau k$, of the time interval, i.e,
$f_{i}(t)=\sum_{k=0}^{N}f_{i}^{k}\psi^{k}(t)$ with
$\psi^{0}(t)=\chi_{[0,\tau/2)}$, $\psi^{N}(t)=\chi_{[T-\tau/2,T)}$,
$\psi^{k}(t)=\chi_{[t_{k}-\tau/2,t_{k}+\tau/2)}$, $k=1,\ldots,N-1$, so that
with

$f_{h}(t,x)=\sum_{k=0}^{N}\sum_{i=1}^{n_{el}}f_{i}^{k}\psi^{k}(t)\phi_{i}(x)\,,\quad
g_{h}(x)=\sum_{i=1}^{n_{el}}g_{i}\phi_{i}(x)\,,$

the discretized eigenvalue equation in variational form reads as

$\displaystyle\left\langle\widetilde{\mathcal{G}^{*}\mathcal{G}}\left(\begin{array}[]{c}f_{h}\\\
g_{h}\end{array}\right),\left(\begin{array}[]{c}\psi^{\ell}\phi_{j}\\\
\phi_{j}\end{array}\right)\right\rangle_{L^{2}(L^{2})\times
L^{2}}=\left\langle\left(\begin{array}[]{c}f_{h}\\\
g_{h}\end{array}\right),\left(\begin{array}[]{c}\psi^{\ell}\phi_{j}\\\
\phi_{j}\end{array}\right)\right\rangle_{L^{2}(L^{2})\times L^{2}},$ (70)
$\displaystyle\qquad\forall\ell\in\\{0,\ldots,N\\},j\in\\{1,\ldots,n_{el}\\},$
(71)

with $\widetilde{\mathcal{G}^{*}\mathcal{G}}$ as in (51), i.e., in matrix-
vector form

$\left(\begin{array}[]{cc}\mathcal{M}+\mathcal{E}&(\mathcal{M}^{1}+\mathcal{E}^{1})^{T}\\\
\mathcal{M}^{1}+\mathcal{E}^{1}&T\,M+M^{1/2}A_{h}^{-1}M^{1/2}\end{array}\right)\left(\begin{array}[]{c}\underline{\underline{f}}\\\
\underline{g}\end{array}\right)=\lambda\left(\begin{array}[]{cc}\mathcal{M}&0\\\
0&M\end{array}\right)\left(\begin{array}[]{c}\underline{\underline{f}}\\\
\underline{g}\end{array}\right)\,,$

with the matrices

$\displaystyle\mathcal{M}$ $\displaystyle=$
$\displaystyle\tau\mbox{diag}({\textstyle\frac{1}{2}}M,M,\cdots,M,{\textstyle\frac{1}{2}}M)\,,$
$\displaystyle\mathcal{M}^{1}$ $\displaystyle=$
$\displaystyle\tau({\textstyle\frac{1}{2}}M,M,\cdots,M,{\textstyle\frac{1}{2}}M)\,,$
$\displaystyle\mathcal{E}_{\ell,k}$ $\displaystyle=$
$\displaystyle\tau^{2}M^{1/2}A_{h}e^{-A_{h}(2T-t_{\ell}-t_{k})}M^{1/2},$
$\displaystyle\mathcal{E}_{0,k}$ $\displaystyle=$
$\displaystyle\mathcal{E}_{k,0}={\textstyle\frac{\tau^{2}}{2}}M^{1/2}A_{h}e^{-A_{h}(2T-t_{k})}M^{1/2},\quad\mathcal{E}_{0,0}={\textstyle\frac{\tau^{2}}{4}}M^{1/2}A_{h}e^{-A_{h}2T}M^{1/2},$
$\displaystyle\mathcal{E}_{N,k}$ $\displaystyle=$
$\displaystyle\mathcal{E}_{k,N}={\textstyle\frac{\tau^{2}}{2}}M^{1/2}A_{h}e^{-A_{h}(T-t_{k})}M^{1/2},\quad\mathcal{E}_{N,N}={\textstyle\frac{\tau^{2}}{4}}M^{1/2}A_{h}M^{1/2},$
$\displaystyle\mathcal{E}^{1}$ $\displaystyle=$
$\displaystyle\tau({\textstyle\frac{1}{2}}M^{1/2}e^{-A_{h}T}M^{1/2},\ldots,M^{1/2}e^{-A_{h}(T-t_{k})}M^{1/2},\ldots,{\textstyle\frac{1}{2}}M)$

for $\ell,k\in\\{1,\ldots,N-1\\}$. Here we have computed and approximated the
time integrals in (70) as follows

$\displaystyle\int_{t_{k}-\tau/2}^{t_{k}+\tau/2}\mbox{ ds}=\tau,$
$\displaystyle\int_{t_{k}-\tau/2}^{t_{k}+\tau/2}e^{-A_{h}(T-s)}\mbox{
ds}=A_{h}^{-1}\left(e^{A_{h}\tau/2}-e^{-A_{h}\tau/2}\right)e^{-A_{h}(T-t_{k})}\approx\tau
e^{-A_{h}(T-t_{k})},$
$\displaystyle\int_{t_{\ell}-\tau/2}^{t_{\ell}+\tau/2}e^{-A_{h}(T-t)}A_{h}\int_{t_{k}-\tau/2}^{t_{k}+\tau/2}e^{-A_{h}(T-s)}\mbox{
ds}\mbox{
dt}=A_{h}^{-1}\left(e^{A_{h}\tau}-2+e^{-A_{h}\tau}\right)e^{-A_{h}(2T-t_{\ell}-t_{k})}$
$\displaystyle\quad\approx\tau^{2}A_{h}e^{-A_{h}(2T-t_{\ell}-t_{k})}\,.$

To simplify the implementation we transform the eigenvalue equation by pre-
multiplication with $\left(\begin{array}[]{cc}\mathcal{M}&0\\\
0&M\end{array}\right)^{-\frac{1}{2}}$ and setting
$\left(\begin{array}[]{c}\underline{\underline{\tilde{f}}}\\\
\underline{\tilde{g}}\end{array}\right)=\left(\begin{array}[]{cc}\mathcal{M}&0\\\
0&M\end{array}\right)^{\frac{1}{2}}\left(\begin{array}[]{c}\underline{\underline{f}}\\\
\underline{g}\end{array}\right)$ which results in the following equation

$\displaystyle\left(\begin{array}[]{cc}\mathcal{I}+\tilde{\mathcal{E}}&(\mathcal{I}^{1}+\tilde{\mathcal{E}}^{1})^{T}\\\
\mathcal{I}^{1}+\tilde{\mathcal{E}}^{1}&TI+A_{h}^{-1}\end{array}\right)\left(\begin{array}[]{c}\underline{\underline{\tilde{f}}}\\\
\underline{\tilde{g}}\end{array}\right)=\lambda\left(\begin{array}[]{c}\underline{\underline{\tilde{f}}}\\\
\underline{\tilde{g}}\end{array}\right)\,,$ (78)

with the matrices

$\displaystyle\mathcal{I}$ $\displaystyle=$
$\displaystyle\mbox{diag}(I,\cdots,I)\,,$ $\displaystyle\mathcal{I}^{1}$
$\displaystyle=$
$\displaystyle\sqrt{\tau}({\textstyle\frac{1}{\sqrt{2}}}I,I,\cdots,I,{\textstyle\frac{1}{\sqrt{2}}}I)\,,$
$\displaystyle\tilde{\mathcal{E}}_{\ell,k}$ $\displaystyle=$
$\displaystyle\tau A_{h}e^{-A_{h}(2T-t_{\ell}-t_{k})},$
$\displaystyle\tilde{\mathcal{E}}_{0,k}$ $\displaystyle=$
$\displaystyle\tilde{\mathcal{E}}_{k,0}={\textstyle\frac{\tau}{\sqrt{2}}}A_{h}e^{-A_{h}(2T-t_{k})},\quad\tilde{\mathcal{E}}_{0,0}={\textstyle\frac{\tau}{2}}A_{h}e^{-A_{h}2T},$
$\displaystyle\tilde{\mathcal{E}}_{N,k}$ $\displaystyle=$
$\displaystyle\tilde{\mathcal{E}}_{k,N}={\textstyle\frac{\tau}{\sqrt{2}}}A_{h}e^{-A_{h}(T-t_{k})},\quad\tilde{\mathcal{E}}_{N,N}={\textstyle\frac{\tau}{2}}A_{h},$
$\displaystyle\tilde{\mathcal{E}}_{0,N}$ $\displaystyle=$
$\displaystyle\tilde{\mathcal{E}}_{N,0}=\tau A_{h}e^{-A_{h}},$
$\displaystyle\tilde{\mathcal{E}}^{1}$ $\displaystyle=$
$\displaystyle\sqrt{\tau}({\textstyle\frac{1}{\sqrt{2}}}e^{-A_{h}T},\ldots,e^{-A_{h}(T-t_{k})},\ldots,{\textstyle\frac{1}{\sqrt{2}}}I).$

With that setting the eigenvalues are computed with the time interval $[0,T]$
with $T=1$ discretized in the points
$t_{0}=0,t_{1}=0.25,t_{2}=0.5,t_{3}=0.75,t_{4}=T=1$. Sampling of 700
eigenvalues leads to the following visualized output, which can be seen in
Figure 2. The eigenvalues tend to zero, in three steps, which was also found
in the analytic considerations.

Figure 2: Eigenvalues of $\widetilde{\mathcal{G}^{*}\mathcal{G}}$ for
backwards heat problem.

##  4 Convergence analysis.

For the convergence analysis, the results in [9, 10] are extended to problems
formulated in an all-at-once fashion, which can be done in a very
straightforward fashion, so we keep this section short by more or less
recalling the essential results from [9, 10].

The aim is to find out about the convergence of the posterior towards the true
element $\mathbf{x}^{*}$ which generates the data $\mathbf{y}^{\delta}$. This
is typically done by analyzing the squared posterior contraction, given by

$\mbox{SPC
}:=\mathbb{E}^{\mathbf{x}^{*}}\mathbb{E}_{\alpha}^{\delta}\|\mathbf{x}^{*}-\mathbf{x}\|,\quad\alpha,\delta>0,$
(79)

where the outward expectation is taken with respect to the data generating
function, which means with given $\mathbf{x}^{*}$, the distribution which
generates the data $\mathbf{y}^{\delta}$. The inward expectation is taken with
respect to the posterior distribution with given data $\mathbf{y}^{\delta}$
and chosen scaling parameter $\alpha$. The squared posterior contraction can
be decomposed into the squared bias, the estimation variance and the posterior
spread, i.e.

$\mbox{SPC
}=\|\mathbf{x}^{*}-\mathbb{E}^{\mathbf{x}^{*}}\mathbf{x}_{\alpha}^{\delta}\|^{2}+\mathbb{E}^{\mathbf{x}^{*}}\|\mathbf{x}_{\alpha}^{\delta}-\mathbb{E}^{\mathbf{x}^{*}}\mathbf{x}_{\alpha}^{\delta}\|^{2}+\mbox{tr
}[\mathcal{C}_{\alpha}^{\delta}],$ (80)

where $\mathbf{x}_{\alpha}^{\delta}$ and $\mathcal{C}_{\alpha}^{\delta}$ are
denoting the posterior mean and covariance, respectively. It suffices to bound
the bias
$b_{\mathbf{x}^{*}}(\alpha)=\|\mathbf{x}^{*}-\mathbb{E}^{\mathbf{x}^{*}}\mathbf{x}_{\alpha}^{\delta}\|$
and the posterior spread $\mbox{tr }[\mathcal{C}_{\alpha}^{\delta}]$, as the
estimation variance
$\mathbb{E}^{\mathbf{x}^{*}}\|\mathbf{x}_{\alpha}^{\delta}-\mathbb{E}^{\mathbf{x}^{*}}\mathbf{x}_{\alpha}^{\delta}\|^{2}$
is always bounded by the posterior spread. To bound the posterior spread, it
is shown that the necessary methods given in [9, 10] can be extended to block
matrix operators. First, the concept of index functions is discussed.

###### Definition 4.1

(index function) A function $\psi:(0,\infty)\to\mathbb{R}^{+}$ is called an
index function if it is a continuous non-decreasing function with $\psi(0)=0$.

To be able to compare index functions the following partial ordering is
defined.

###### Definition 4.2

(partial ordering for index functions) Given two index functions $g,h$ we
write $g\prec h$ if the function $h(t)/g(t)$ is an index function (which means
$h$ tends to zero faster than $g$).

In addition to ordering index function, also self-adjoint operators need to be
partially ordered.

###### Definition 4.3

(partial ordering for self-adjoint operators) Let $\mathcal{G}$ and
$\mathcal{G}^{\prime}$ be bounded self-adjoint opertors in some Hilbert space
$X\times Y$. We say that $\mathcal{G}\leq\mathcal{G}^{\prime}$ if for all
$\mathbf{x}\in X\times Y$ the inequality
$\langle\mathcal{G}\mathbf{x},\mathbf{x}\rangle\leq\langle\mathcal{G}^{\prime}\mathbf{x},\mathbf{x}\rangle$
holds.

The last definition of interest states operator concavity for self-adjoint
operators.

###### Definition 4.4

Let $f:[0,a]\to\mathbb{R}^{+}$ be a continuous function. It is called operator
concave if we have for any pair $\mathcal{G},\mathcal{H}\geq 0$ of self-
adjoint operators with spectra in $[0,a]$

$f\left(\frac{\mathcal{G}+\mathcal{H}}{2}\right)\geq\frac{f(\mathcal{G})+f(\mathcal{H})}{2}.$

The theorem stated below is the range inclusion theorem, which also holds for
block operator matrices.

###### Theorem 4.5

(Douglas’ Range Inclusion Theorem; Theorem 1 in [23]) Let the operators
$\mathcal{S},\mathcal{T}:X\times Y\to X\times Y$ be bounded and act between
Hilbert spaces. Then the following statements are equivalent:

  1. 1.

$\mathcal{R}(\mathcal{S})\subset\mathcal{R}(\mathcal{T})$

  2. 2.

$\mathcal{S}\mathcal{S}^{*}\leq C^{2}\mathcal{T}\mathcal{T}^{*}$, for some
$C\geq 0$

  3. 3.

there exists an bounded operator $\mathcal{R}:X\times Y\to X\times Y$ with
$||\mathcal{R}||\leq C$, such that $\mathcal{S}=\mathcal{T}\mathcal{R}$.

Before continuing we introduce the functions

$\Theta(t)=\Theta_{\psi}(t):=\sqrt{t}\psi(t),\quad t>0,\quad
s_{\alpha}(t):=\frac{\alpha}{\alpha+t},\quad\alpha>0,$ (81)

which will be needed in the further analysis. The main concept of the
convergence analysis in [10] is the fulfillment of two assumptions, namely a
link and a source condition, which will be given next. With these two
assumptions one is able to bound the squared bias and the posterior spread,
and therefore, the squared posterior contraction.

###### Assumption 4.6

(link condition) There is an index function $\psi$, and there are constants
$0<\underline{m}\leq\overline{m}<\infty$ such that

$\underline{m}\|\psi(\mathcal{C}_{0})\mathbf{x}\|_{W^{\prime}\times
Y}\leq\|\Sigma^{1/2}\mathcal{G}\mathbf{x}\|_{U\times
X}\leq\overline{m}\|\psi(\mathcal{C}_{0})\mathbf{x}\|_{W^{\prime}\times
Y},\quad\mathbf{x}\in U\times X.$ (82)

In addition, with the function $\Theta$ from (81) the function

$f_{0}(s):=\left((\Theta^{2})^{-1}(s)\right)^{-1/2},~{}~{}s>0,$

has an operator concave square $f_{0}^{2}$.

For the following analysis we define

$\varphi_{0}(t):=\sqrt{t},~{}~{}t>0$ (83)

and state the following proposition, which will allow us to bound the bias.

###### Proposition 4.7

(Proposition 1 in [10]) Under Assumption 4.6 we have that
$\mathcal{R}(\mathcal{C}_{0}^{1/2})=\mathcal{R}(f_{0}(\mathcal{H}))$. Mainly
we have that the operator
$f_{0}(\mathcal{H})\varphi_{0}(\mathcal{C}_{0})^{-1}$ is norm bounded by
$\overline{m}$.

As mentioned before, the second assumption of importance is a source
condition, which will be stated next.

###### Assumption 4.8

(source set) There is an index function $\varphi$ such that

$\mathbf{x}^{*}\in\mathcal{S}_{\varphi}:=\\{\mathbf{x},\quad\mathbf{x}=\varphi(\mathcal{C}_{0})\mathbf{v},~{}||\mathbf{v}||\leq
1\\}.$ (84)

By using Proposition 4.7 we can bound the bias as

$\displaystyle b_{\mathbf{x}^{*}}(\alpha)$
$\displaystyle\leq\frac{1}{\underline{m}}||f_{0}(\mathcal{H})s_{\alpha}(\mathcal{H})\varphi_{0}(\mathcal{C}_{0})^{-1}\varphi(\mathcal{C}_{0})||$
(85)
$\displaystyle=\frac{1}{\underline{m}}||s_{\alpha}(\mathcal{H})f_{0}(\mathcal{H})\varphi_{0}(\mathcal{C}_{0})^{-1}\varphi(\mathcal{C}_{0})||.$

Now, with the assumptions made, one can give an upper bound for the squared
bias.

###### Proposition 4.9

(Proposition 3 in [10]) Suppose that either
$\varphi_{0}\prec\varphi\prec\theta$, and the function

$g^{2}(t):=\left(\frac{\varphi}{\varphi_{0}}\right)^{2}\left(\left(\theta^{2}\right)^{-1}(t)\right),\quad
t>0,$

is operator concave, or $1\prec\varphi\prec\varphi_{0}$ and $\varphi$ is
operator concave. Under Assumptions 4.6 and 4.8 we have that

$b_{\mathbf{x}^{*}}(\alpha)\leq\frac{\overline{m}}{\underline{m}}||s_{\alpha}(\mathcal{H})\varphi(f_{0}^{2}(\mathcal{H}))||.$
(86)

With the help of Assumption 4.6, one can bound the posterior spread as stated
in the next proposition.

###### Proposition 4.10

(Proposition 6 in [10]) Under Assumption 4.6 we have that

$\mbox{tr}[\mathcal{C}_{\alpha}^{\delta}]\leq\frac{\delta^{2}}{\underline{m}^{2}}\mbox{tr}[(\alpha+\mathcal{H})^{-1}f_{0}^{2}(\mathcal{H})].$
(87)

Under certain circumstances, see [10] and the references therein, we can
rewrite the upper bound for the bias as

$b_{\mathbf{x}^{*}}(\alpha)\leq\frac{\overline{m}}{\underline{m}}\varphi(f_{0}^{2}(\alpha)),\quad\alpha>0,$
(88)

using the fact that $s_{\alpha}(t)t\leq\alpha$, for $t,\alpha>0$.

Combination of the results from Proposition 4.9 and 4.10 leads to a bound on
the squared posterior contraction as stated in the following theorem.

###### Theorem 4.11

(Theorem 3 in [10]) Supposed that Assumption 4.6 and Assumption 4.8 hold for
some index functions $\psi$ and $\varphi$ and additionally, the assumptions in
Proposition 4.9 and holds true, then

$\mbox{SPC}(\alpha,\delta)\leq\frac{1}{\underline{m}^{2}}\Bigl{(}\overline{m}^{2}||s_{\alpha}(\mathcal{H})\varphi(f_{0}^{2}(\mathcal{H}))||+\delta^{2}\mbox{
tr
}[(\alpha+\mathcal{H})^{-1}f_{0}^{2}(\mathcal{H})]\Bigr{)},\quad\alpha,\delta>0.$
(89)

###### Remark 1

If only the lower bound in (82) holds, with (85) and under a source condition
one still obtains the estimate

$\mbox{SPC}(\alpha,\delta)=\frac{1}{\underline{m}^{2}}\Bigl{(}||s_{\alpha}(\mathcal{H})f_{0}(\mathcal{H})\varphi_{0}(\mathcal{C}_{0})^{-1}\varphi(\mathcal{C}_{0})||+\delta^{2}\mbox{
tr }[(\alpha+\mathcal{H})^{-1}f_{0}^{2}(\mathcal{H})]\Bigr{)},$ (90)

for $\alpha,\delta>0.$

###  4.1 Fulfillment of the link condition for the all-at-once-formulation.

A straightforward example that always allows to fulfill the link condition, is
as follows. Choosing the prior covariance to be given by

$\mathcal{C}_{0}:=\mathcal{G}^{*}\mathcal{G},$ (91)

and the link condition is always fulfilled by choosing the link function to be
given by

$\psi(t):=t^{1/2}.$ (92)

Assume $\mathcal{G}:X\times Y\to S\times T$, then $\mathcal{G}^{*}:S\times
T\to X\times Y$ and $\mathcal{G}^{*}\mathcal{G}:X\times Y\to X\times Y$ which
yields

$\|\psi(\mathcal{G}^{*}\mathcal{G})\mathbf{x}\|_{X\times
Y}=\|\mathcal{G}\mathbf{x}\|_{S\times T},\quad\mbox{ for all }\mathbf{x}\in
X\times Y,$ (93)

and the link condition is fulfilled. However, since
$\mathcal{G}^{*}\mathcal{G}$ is ill-posed, (91) will not be a good
regularizer, so we look at further possible choices for our two examples in
the following Section 5.

In particular, for the backwards heat equation, the severe ill-posedness is
provibitive for (91). In the Section 5.2 below, we therefore derive a prior
that is motivated by the link condition but mach easier to handle than (91).

##  5 Choice of joint priors.

###  5.1 Block diagonal priors satisfying unilateral link estimates.

In this section, we investigate to which extent the link condition (82) can be
satisfied by block diagonal operators $\psi(\mathcal{C}_{0})$, as these are
convenient for applying $\psi^{-1}$ in order to obtain the prior covariance
$\mathcal{C}_{0}$ itself again in block diagonal form, which in its turn is
useful for numerical computations. For this purpose we again focus on our two
prototypical examples from Sections 2.1 and 2.2. We will consider the upper
and lower bound in (82) separately and, assuming that $\Sigma$ is an
isomorphism on $Y$, without loss of generality, set it to the identity, i.e.,
we will investigate fulfillment of

$\|\psi(\mathcal{C}_{0})\mathbf{x}\|_{U\times
X}\leq\frac{1}{\underline{m}}\|\mathcal{G}\mathbf{x}\|_{W^{\prime}\times
Y}\mbox{ for all }\mathbf{x}=(u,\theta)\in U\times X,$ (94)
$\|\mathcal{G}\mathbf{x}\|_{W^{\prime}\times
Y}\leq\overline{m}\|\psi(\mathcal{C}_{0})\mathbf{x}\|_{U\times X}\mbox{ for
all }\mathbf{x}=(u,\theta)\in U\times X,$ (95)

for $\psi(\mathcal{C}_{0})$ of the form

$\psi(\mathcal{C}_{0})=\left(\begin{array}[]{cc}\mathcal{B}&0\\\
0&\mathcal{D}\end{array}\right).$ (96)

For the inverse source problem from Section 2.1 we have

$\|\mathcal{G}(u,\theta)\|_{W^{\prime}\times
Y}=\|\mathcal{A}u-\theta\|_{L^{2}}+\|u\|_{L^{2}}$

and

$\|\psi(\mathcal{C}_{0})(u,\theta)\|_{U\times
X}=\|\mathcal{A}\mathcal{B}u\|_{L^{2}}+\|\mathcal{D}\theta\|_{L^{2}}$

so setting $\theta$ and $u$ to zero separately we get the following two
necessary conditions

$\displaystyle\|\mathcal{A}u\|_{L^{2}}\leq\overline{m}\|\mathcal{A}\mathcal{B}u\|_{L^{2}}\mbox{
for all }u\in U=H_{0}^{1}(\Omega)\cap H^{2}(\Omega)\,,$ i.e.,
$\mathcal{B}:U\to U$ bijective with bounded inverse $\mathcal{B}^{-1}$ (97)

and

$\displaystyle\|\theta\|_{L^{2}}\leq\overline{m}\|\mathcal{D}\theta\|_{L^{2}}\mbox{
for all }\theta\in X=L^{2}(\Omega)\,,$ i.e., $\mathcal{D}:L^{2}(\Omega)\to
L^{2}(\Omega)$ bijective with bounded inverse $\mathcal{D}^{-1}$ (98)

for (95). On the other hand, an easy application of the triangle inequality
together with boundedness of $\mathcal{A}^{-1}:L^{2}(\Omega)\to L^{2}(\Omega)$
implies sufficiency of (5.1), (5.1) for (95).
The lower bound (94), with $\theta=\mathcal{A}u$ implies

$\|\mathcal{A}\mathcal{B}u\|_{L^{2}}+\|\mathcal{D}\mathcal{A}u\|_{L^{2}}\leq\frac{1}{\underline{m}}\|u\|_{L^{2}}\mbox{
for all }u\in U=H_{0}^{1}(\Omega)\cap H^{2}(\Omega)\,,$

thus since $\|\mathcal{D}\mathcal{A}\|_{L^{2}\to
L^{2}}=\|\mathcal{A}\mathcal{D}^{*}\|_{L^{2}\to L^{2}}$ we get the necessary
conditions

$\mathcal{B}:L^{2}(\Omega)\to H_{0}^{1}(\Omega)\cap H^{2}(\Omega)$ and
$\mathcal{D}^{*}:L^{2}(\Omega)\to H_{0}^{1}(\Omega)\cap H^{2}(\Omega)$ bounded
(99)

for (94). To see sufficiency of (99) together with boundedness of
$\mathcal{D}:L^{2}(\Omega)\to L^{2}(\Omega)$ for (94), we set
$f=\mathcal{A}u-\theta$ to obtain

$\displaystyle\|\psi(\mathcal{C}_{0})(u,\theta)\|_{U\times X}$
$\displaystyle=$
$\displaystyle\|\mathcal{A}\mathcal{B}u\|_{L^{2}}+\|\mathcal{D}(\mathcal{A}u-f)\|_{L^{2}}$
$\displaystyle\leq$ $\displaystyle(\|\mathcal{B}\|_{L^{2}\to
H^{2}}+\|\mathcal{D}^{*}\|_{L^{2}\to
H^{2}})\|u\|_{L^{2}}+\|\mathcal{D}\|_{L^{2}\to L^{2}})\|f\|_{L^{2}}$
$\displaystyle\leq$
$\displaystyle\frac{1}{\underline{m}}(\|f\|_{L^{2}}+\|u\|_{L^{2}})=\frac{1}{\underline{m}}\|\mathcal{G}(u,\theta)\|_{W^{\prime}\times
Y}$

with $\underline{m}:=\frac{1}{\max\\{\|\mathcal{B}\|_{L^{2}\to
H^{2}}+\|\mathcal{D}^{*}\|_{L^{2}\to H^{2}},\,\|\mathcal{D}\|_{L^{2}\to
L^{2}}\\}}$.
Unfortunately, conditions (5.1) and (99) are contradictory, so the full
squared posterior contraction estimate (89) cannot be applied. Still, with,
e.g. $\psi(\mathcal{C}_{0})=\left(\begin{array}[]{cc}\mathcal{A}^{-s}&0\\\
0&\mathcal{A}^{-p}\end{array}\right)$ and $s,p\geq 1$, which satisfies (94),
we obtain (90).
For the backwards heat problem from Section 2.1 with

$\|\mathcal{G}(u,\theta)\|_{W^{\prime}\times
Y}=\|\mathcal{A}^{-1/2}(\partial_{t}+\mathcal{A})u+\mathcal{A}^{1/2}\theta\|_{L^{2}(L^{2})}+\|u(T)+\theta\|_{L^{2}}$

and

$\|\psi(\mathcal{C}_{0})(u,\theta)\|_{U\times
X}=\|\mathcal{A}^{-1/2}(\partial_{t}+\mathcal{A})\mathcal{B}u\|_{L^{2}(L^{2})}+\|\mathcal{A}^{1/2}\mathcal{D}\theta\|_{L^{2}}$

we see, by considering the two special cases $\theta=0$ and $u=0$, that a
necessary condition for (95) is

$\mathcal{B}^{-1}:U_{0}\to U_{0}$ and $\mathcal{D}^{-1}:H_{0}^{1}(\Omega)\to
H_{0}^{1}(\Omega)$ bounded. (100)

Indeed, as can be seen by the triangle inequality and fact that that
$\|u(T)\|_{L^{2}}\leq\|\mathcal{A}^{-1/2}(\partial_{t}+\mathcal{A})u\|_{L^{2}(L^{2})}$
(cf. (35)) condition (100) is also sufficient for (95).
The lower bound (94) is again more challenging to obtain. With the particular
choices $\theta=0$ on one hand and
$(\partial_{t}+\mathcal{A})u=-\mathcal{A}\theta$, i.e.,
$u(t)=-(I-e^{-t\mathcal{A}})\theta$ implying
$u(T)+\theta=e^{-T\mathcal{A}}\theta$ on the other hand, it yields the
necesssary conditions

$\displaystyle\|\mathcal{B}u\|_{U_{0}}\leq\frac{1}{\underline{m}}\Bigl{(}\|u\|_{U_{0}}+\|u(T)\|_{L^{2}}\Bigr{)}\leq\frac{2}{\underline{m}}\|u\|_{U_{0}}\mbox{
for all }u\in U_{0},$ i.e., $\mathcal{B}:U_{0}\to U_{0}$ bounded

and

$\displaystyle\|\mathcal{B}(\underline{I}-e^{-\cdot\mathcal{A}})\theta\|_{U_{0}}+\|\mathcal{D}\theta\|_{H_{0}^{1}}\leq\frac{1}{\underline{m}}\|e^{-T\mathcal{A}}\theta\|_{L^{2}}\mbox{
for all }\theta\in X,$ i.e., $\mathcal{D}e^{T\mathcal{A}}:L^{2}(\Omega)\to
H_{0}^{1}(\Omega)$ bounded.

To obtain sufficient conditions for (94), with
$f=(\partial_{t}+\mathcal{A})u+\mathcal{A}\theta$, $g=u(T)+\theta$, so that
$\|\mathcal{G}(u,\theta)\|_{W^{\prime}\times
Y}=\|f\|_{L^{2}(H^{-1})}+\|g\|_{L^{2}}$,
$u(t)=\int_{0}^{t}e^{-(t-s)\mathcal{A}}f(s)\,ds-(I-e^{-t\mathcal{A}})\theta$,
hence we can express $u$ and $\theta$ via $f$ and $g$ as

$\displaystyle u(t)$ $\displaystyle=$
$\displaystyle\int_{0}^{t}e^{-(t-s)\mathcal{A}}f(s)\,ds-(I-e^{-t\mathcal{A}})(e^{T\mathcal{A}}g-\int_{0}^{T}e^{s\mathcal{A}}f(s)\,ds),$
$\displaystyle\theta$ $\displaystyle=$ $\displaystyle
e^{T\mathcal{A}}g-\int_{0}^{T}e^{s\mathcal{A}}f(s)\,ds,$

and, therefore ,

$\displaystyle\|\psi(\mathcal{C}_{0})(u,\theta)\|_{U\times X}$
$\displaystyle=\|\mathcal{B}u\|_{U_{0}}+\|\mathcal{D}\theta\|_{H_{0}^{1}}$
$\displaystyle\leq\|\mathcal{B}\|_{L^{\infty}(L^{2})\to
U_{0}}\|\int_{0}^{\cdot}e^{-(\cdot-s)\mathcal{A}}f(s)\,ds\|_{L^{\infty}(L^{2})}$
$\displaystyle\quad+\Bigl{(}\|\mathcal{B}e^{T\mathcal{A}}\|_{L^{\infty}(L^{2})\to
U_{0}}\|\underline{I}-e^{-\cdot\mathcal{A}}\|_{L^{2}\to
L^{\infty}(L^{2})}+\|\mathcal{D}e^{T\mathcal{A}}\|_{L^{2}\to
H_{0}^{1}}\Bigr{)}$
$\displaystyle\qquad\qquad\cdot\Bigl{(}\|g\|_{L^{2}}+\|\int_{0}^{T}e^{-(T-s)\mathcal{A}}f(s)\,ds\|_{L^{2}}\Bigr{)}\,,$

where $\|\underline{I}-e^{-\cdot\mathcal{A}}\|_{L^{2}\to
L^{\infty}(L^{2})}\leq 1$ and by self-adjoinedness of $\mathcal{A}$ as well as
the Cauchy-Scharz inequality

$\displaystyle\|\int_{0}^{\cdot}e^{-(\cdot-s)\mathcal{A}}f(s)\,ds\|_{L^{\infty}(L^{2})}$
$\displaystyle=\sup_{t\in[0,T]\,v\in
L^{2}\setminus\\{0\\}}{\textstyle\frac{1}{\|v\|_{L^{2}}}}\int_{0}^{t}\int_{\Omega}e^{-(t-s)\mathcal{A}}\mathcal{A}^{1/2}v\,\mathcal{A}^{-1/2}f(s)\,dx\,ds$
$\displaystyle\leq\sup_{t\in[0,T]\,v\in
L^{2}\setminus\\{0\\}}{\textstyle\frac{1}{\|v\|_{L^{2}}}}\left(\int_{0}^{t}\sum_{n=1}^{\infty}e^{-2(t-s)\frac{1}{\mu_{n}}}\frac{1}{\mu_{n}}\langle
v,\phi_{n}\rangle_{L^{2}}^{2}\,ds\right)^{1/2}\
\|\mathcal{A}^{-1/2}f\|_{L^{2}(L^{2})}$
$\displaystyle\leq\frac{1}{2}\|f\|_{L^{2}(H^{-1})}$

with an eigensystem
$(\mu_{n},\phi_{n})_{n\in\mathbb{N}}\subseteq\mathbb{R}\times L^{2}(\Omega)$
of $\mathcal{A}^{-1}$; likewise
$\|\int_{0}^{T}e^{-(T-s)\mathcal{A}}f(s)\,ds\|_{L^{2}}\leq\frac{1}{2}\|f\|_{L^{2}(H^{-1})}$.
Since these estimates are sharp when having to hold for arbitrary $f\in
L^{2}(H^{-1})$, we can conclude that validity of (94) is equivalent to

$\mathcal{B}e^{T\mathcal{A}}:L^{\infty}(L^{2})\to U_{0}$ and
$\mathcal{D}e^{T\mathcal{A}}:L^{2}(\Omega)\to H_{0}^{1}(\Omega)$ bounded.
(101)

Again the conditions for (94) and (95) are unfortunately contradictory. Still,
for example the choice
$\psi(\mathcal{C}_{0})=\left(\begin{array}[]{cc}\mathcal{A}^{1/2}e^{-T\mathcal{A}}\int_{0}^{\cdot}e^{-(\cdot-s)\mathcal{A}}\cdot(s)\,ds&0\\\
0&\mathcal{A}^{-1/2}e^{-T\mathcal{A}}\end{array}\right)$ satisfies (94) and
therewith (90).

###  5.2 Heuristic choice of $\mathcal{C}_{0}$ for the backwards heat problem.

To find out more about suitable operators fulfilling the link condition
besides the trivial choice one attempt is made in the setting of the backwards
heat problem. To do so, the norms are compared and so, the operator
$\psi(\mathcal{C}_{0})$ can be given by

$\psi(\mathcal{C}_{0})=\left(\begin{array}[]{cc}I&\underline{I}-e^{-\mathcal{A}t}\\\
\mathcal{A}^{-1}\delta_{T}+B&\mathcal{A}^{-1/2}\end{array}\right),$

where $\mathcal{A}:=-\Delta$, $Bu:=\int_{0}^{T}u(t)\mbox{ dt}$ and
$(\underline{I}\theta)(t)=\theta$. This is a symmetric operator from
$U_{0}\times H_{0}^{1}(\Omega)\to U_{0}\times H_{0}^{1}(\Omega)$ as

$\displaystyle((I-e^{-\mathcal{A}t})u,v)_{U_{0}}$
$\displaystyle=\int_{0}^{T}\int_{\Omega}\mathcal{A}^{-1/2}(\partial_{t}+\mathcal{A})(\underline{I}-e^{-\mathcal{A}t})u\mathcal{A}^{-1/2}(\partial_{t}+\mathcal{A})v\mbox{
dx}\mbox{ dt}$
$\displaystyle=\int_{0}^{T}\int_{\Omega}\mathcal{A}^{1/2}u\mathcal{A}^{-1/2}(\partial_{t}+\mathcal{A})v\mbox{
dx}\mbox{ dt}$
$\displaystyle=\int_{0}^{T}\int_{\Omega}\mathcal{A}^{1/2}u\mathcal{A}^{-1/2}\partial_{t}v\mbox{
dx}\mbox{
dt}+\int_{0}^{T}\int_{\Omega}\mathcal{A}^{1/2}u\mathcal{A}^{1/2}v\mbox{
dx}\mbox{ dt}$
$\displaystyle=\int_{\Omega}\mathcal{A}^{1/2}u\mathcal{A}^{-1/2}v(T)\mbox{
dx}+\int_{\Omega}\mathcal{A}^{1/2}u\mathcal{A}^{1/2}\int_{0}^{T}v(t)\mbox{
dt}\mbox{ dx},$
$\displaystyle(u,(\mathcal{A}^{-1}\delta_{T}+B)v)_{H_{0}^{1}(\Omega)}$
$\displaystyle=\int_{\Omega}\mathcal{A}^{1/2}u\mathcal{A}^{1/2}(\mathcal{A}^{-1}\delta_{T}+B)v\mbox{
dx}$
$\displaystyle=\int_{\Omega}\mathcal{A}^{1/2}u\mathcal{A}^{-1/2}v(T)\mbox{
dx}+\int_{\Omega}\mathcal{A}^{1/2}u\mathcal{A}^{1/2}\int_{0}^{T}v(t)\mbox{
dt}.$

Now the norm of $\mathcal{G}$ and $\psi(\mathcal{C})$ are computed with the
help of a rewriting of $\psi(\mathcal{C})$ as

$\displaystyle\psi(\mathcal{C}_{0})$ $\displaystyle=$
$\displaystyle\left(\begin{array}[]{cc}I&\underline{I}-e^{-\mathcal{A}t}\\\
\mathcal{A}^{-1}\delta_{T}+B&\mathcal{A}^{-1/2}\end{array}\right)$ (104)
$\displaystyle=$ $\displaystyle\underbrace{\left(\begin{array}[]{cc}I&0\\\
0&\mathcal{A}^{-1/2}\end{array}\right)}_{=:\tilde{\mathcal{A}}}\underbrace{\left(\begin{array}[]{cc}I&\underline{I}-e^{-\mathcal{A}t}\\\
\mathcal{A}^{-1/2}\delta_{T}&\mathcal{A}^{-1/2}\end{array}\right)}_{\tilde{\mathcal{C}}}+\underbrace{\left(\begin{array}[]{cc}0&0\\\
B&\mathcal{A}^{-1/2}-\mathcal{A}^{-1}\end{array}\right)}_{=:\mathcal{R}}.$
(111)

Instead of $\psi(\mathcal{C}_{0})$ now we use $\tilde{\mathcal{C}}$ under the
norm. The norms compute as follows for $x=(u,\theta)\in U_{0}\times H_{0}^{1}$

$\displaystyle\|\mathcal{G}\textbf{x}\|_{L^{2}(H^{-1})\times L^{2}}^{2}$
$\displaystyle=$
$\displaystyle\|(\partial_{t}+\mathcal{A})u+\mathcal{A}\theta||_{L^{2}(H^{-1})}^{2}+\|\delta_{T}u+I\theta\|_{L^{2}}^{2}$
$\displaystyle=$
$\displaystyle\|(\partial_{t}+\mathcal{A})u\|_{L^{2}(H^{-1})}^{2}+\|\mathcal{A}\theta\|_{L^{2}(H^{-1})}^{2}+2((\partial_{t}+\mathcal{A})u,\mathcal{A}\theta)_{L^{2}(H^{-1})}$
$\displaystyle+\|\delta_{T}u\|_{L^{2}}^{2}+\|I\theta\|_{L^{2}}^{2}+2(\delta_{T}u,I\theta)_{L^{2}}$
$\displaystyle=$
$\displaystyle\int_{0}^{T}\int_{\Omega}|\mathcal{A}^{-1/2}(\partial_{t}+\mathcal{A})u|^{2}\mbox{
dx}\mbox{ dt}+\int_{0}^{T}\int_{\Omega}|\mathcal{A}^{1/2}\theta|^{2}\mbox{
dx}\mbox{ dt}$
$\displaystyle+2\int_{0}^{T}\int_{\Omega}\mathcal{A}^{-1/2}(\partial_{t}+\mathcal{A})u\mathcal{A}^{1/2}\theta\mbox{
dx}\mbox{ dt}+\int_{\Omega}|u(T)|^{2}\mbox{ dx}$
$\displaystyle+\int_{\Omega}|\theta|^{2}\mbox{
dx}+2\int_{\Omega}u(T)I\theta\mbox{ dx},$

$\displaystyle\|\tilde{\mathcal{C}}\textbf{x}\|_{U_{0}\times H_{0}^{1}}^{2}$
$\displaystyle=$
$\displaystyle\|u+(\underline{I}-e^{-\mathcal{A}t}\theta)||_{U_{0}}^{2}+\|\mathcal{A}^{-1/2}\delta_{T}u+\mathcal{A}^{-1/2}\theta\|_{H_{0}^{1}}^{2}$
$\displaystyle=$
$\displaystyle\|u\|_{U_{0}}^{2}+\|(I-e^{-\mathcal{A}t})\theta\|_{U_{0}}^{2}+2(u,(\underline{I}-e^{-\mathcal{A}t})\theta)_{U_{0}}$
$\displaystyle+\|\mathcal{A}^{-1/2}\delta_{T}u\|_{H_{0}^{1}}^{2}+\|\mathcal{A}^{-1/2}\theta\|_{H_{0}^{1}}^{2}+2(\mathcal{A}^{-1/2}\delta_{T}u,\mathcal{A}^{-1/2}\theta)_{H_{0}^{1}}$
$\displaystyle=$
$\displaystyle\int_{0}^{T}\int_{\Omega}|\mathcal{A}^{-1/2}(\partial_{t}+\mathcal{A})u|^{2}\mbox{
dx}\mbox{ dt}$
$\displaystyle+\int_{0}^{T}\int_{\Omega}|\mathcal{A}^{-1/2}(\partial_{t}+\mathcal{A})(I-e^{-\mathcal{A}t})\theta|^{2}\mbox{
dx}\mbox{ dt}$
$\displaystyle+2\int_{0}^{T}\int_{\Omega}\mathcal{A}^{-1/2}(\partial_{t}+\mathcal{A})u\mathcal{A}^{-1/2}(\partial_{t}+\mathcal{A})(\underline{I}-e^{-\mathcal{A}t})\theta\mbox{
dx}\mbox{ dt}$
$\displaystyle+\int_{\Omega}|\mathcal{A}^{1/2}\mathcal{A}^{-1/2}\delta_{T}u|^{2}\mbox{
dx}$
$\displaystyle+\int_{\Omega}|\mathcal{A}^{1/2}\mathcal{A}^{-1/2}\theta|^{2}\mbox{
dx}+2\int_{\Omega}\mathcal{A}^{1/2}\mathcal{A}^{-1/2}\delta_{T}u\mathcal{A}^{1/2}\mathcal{A}^{-1/2}\theta\mbox{
dx}$ $\displaystyle=$
$\displaystyle\int_{0}^{T}\int_{\Omega}|\mathcal{A}^{-1/2}(\partial_{t}+\mathcal{A})u|^{2}\mbox{
dx}\mbox{ dt}+\int_{0}^{T}\int_{\Omega}|\mathcal{A}^{1/2}\theta|^{2}\mbox{
dx}\mbox{ dt}$
$\displaystyle+2\int_{0}^{T}\int_{\Omega}\mathcal{A}^{-1/2}(\partial_{t}+\mathcal{A})u\mathcal{A}^{1/2}\theta\mbox{
dx}\mbox{ dt}+\int_{\Omega}|u(T)|^{2}\mbox{ dx}$
$\displaystyle+\int_{\Omega}|\theta|^{2}\mbox{ dx}+2\int_{\Omega}v(T)u\mbox{
dx}.$

This shows that
$\|\tilde{\mathcal{C}}\mathbf{x}\|^{2}=\|\mathcal{G}\mathbf{x}\|^{2}$.
Motivated by this, we use $\tilde{\mathcal{C}}$ in place of $\mathcal{C}$ as a
prior covariance matrix. Note that estimating the remainder
$\|\mathcal{R}\mathbf{x}\|$ by $\|\mathcal{G}\mathbf{x}\|$ unfortunately does
not seem to be possible, since in $\mathcal{G}\mathbf{x}$ terms containing $u$
and terms containing $\theta$ may possibly cancel. Nevertheless, since
$\mathcal{R}$ contains negative order differential operators as compared to
$\tilde{\mathcal{C}}$, we can regard it as a perturbation of the latter and
therefore skip it for computational purposes.

In Subsection 6.2.4 below we demonstrate numerically that block diagonal of
$\tilde{\mathcal{C}}$ as a prior yields reasonable results in the
reconstruction.

###  5.3 Priors for the inverse source problem.

The choice of joint priors in this work relies on priors which are already
well known in the literature and also implemented in hippylib, see, e.g.
[8][25]. We restrict ourselves to normal priors with the covariance operator
given by

$C_{p}=(-\gamma\Delta+\delta I)^{-n}$ (112)

where $n$ most often take the values $1$ or $2$. For joint priors for the
parameter and the state, priors of the same type have shown to be useful in
the numerical experiments in this study.

###  5.4 Prior for the state variable of the backwards heat problem.

For the backwards heat problem we require a prior which is suitable for the
initial condition, but in addition, we also need to find a prior for the state
variable for every $t\in[0,T]$. To do so, the heat equation is analyzed. At
first the solution to the heat equation is considered and the homogeneous case
$f=0$ is examined. In case of bounded $\Omega$ with, e.g., homogeneous
Dirichlet boundary conditions and denoting by $\mathcal{A}=-\Delta$ the
Laplace operator equipped with these boundary conditions, we can use semigroup
theory to express $u$. To use the semigroup theory some preparatory work has
to be done.

###### Definition 5.1

(Section 7.4.1 in [26]) Let $X$ be a real Banach space. A family of linear,
bounded operators $\\{S(t)\\}_{t\geq 0}$, $S(t):X\to X$ is called a semigroup
iff

$S(0)=I\quad\mbox{and}\quad S(s+t)=S(s)S(t)\quad\forall t,s\geq 0.$

As already assumed the solution of the heat equation with initial data
$\theta$ can be written in terms of semigroups using the so called heat
semigroup

$S(t)\theta:=u(t)\mbox{ where }\left\\{\begin{array}[]{rcl}u_{t}&=&\Delta
u\mbox{ in }\Omega\times(0,\infty)\\\ u&=&0\mbox{ on
}\partial\Omega\times(0,\infty)\\\ \quad u&=&\theta\mbox{ on
}\Omega\times\\{0\\}\end{array}\right.$ (113)

with $X=L^{2}(\Omega)$. Clearly, $S(t)$ fulfills the conditions to be a
semigroup. To investigate this equation the definition of the infinitesimal
generator is needed.

###### Definition 5.2

(Section 7.4.1 in [26]) Let $\\{S(t)\\}_{t\geq 0}$ be a semigroup on a Banach
space $X$. Write

$D(A):=\left\\{u\in X:\lim_{t\to 0^{+}}\frac{S(t)u-u}{t}\mbox{ exists in
}X\right\\}$

and

$Au:=\lim_{t\to 0+}\frac{S(t)u-u}{t}\quad(u\in D(A)).$

We call $A:D(A)\to X$ the infinitesimal generator of the semigroup
$\\{S(t)\\}_{t\geq 0}$, $D(A)$ is the domain of $A$.

It can be shown, that the generator of the heat semigroup is defined by the
Laplace operator with homogeneous Dirichet boundary conditions $\mathcal{A}$.
Therewith, we can write $u$ as

$u(t,x)=(S(t)\theta)(x)=(e^{-t\mathcal{A}}\theta)(x).$ (114)

In terms of finding a prior for the state $u(t,x)$ at any time $t$, this is a
convenient setting as we can make use of the prior measure for $\theta(x)$. We
denote the measure for the parameter $\theta$ by $\mu_{0}$ and given by

$\mu_{0}=\mathcal{N}(m_{0},C_{p})$ (115)

Using the properties of the normal distribution one can calculate the
resulting prior for $u(t,x)$. In addition the fact that $e^{-t\mathcal{A}}$ is
a self-adjoint operator is used. Then the prior for $u(t,x)$ is given by

$\mu_{s}=\mathcal{N}(e^{-t\mathcal{A}}m_{0},e^{-t\mathcal{A}}C_{p}e^{-t\mathcal{A}}).$
(116)

In case $f\neq 0$, one can show that the solution to the heat equation changes
to

$u(t,x)=S(t)\theta(x)+\int_{0}^{t}S(t-s)f(s,x)\mbox{
ds}=e^{-t\mathcal{A}}\theta(x)+\int_{0}^{t}e^{-(t-s)\mathcal{A}}f(s,x)\mbox{
ds}.$ (117)

In the corresponding prior measure only the expectation value has changed, as
the covariance operator is invariant under translation. This leads to the new
prior

$\mu_{s}=\mathcal{N}(e^{-t\mathcal{A}}m_{0}+\int_{0}^{t}e^{-(t-s)\mathcal{A}}f(s,x)\mbox{
ds},e^{-t\mathcal{A}}C_{p}e^{-t\mathcal{A}}).$ (118)

##  6 Numerical Experiments.

In this section the prototypical inverse problems discussed in the previous
section will be solved numerically. To do so, first the theoretical
considerations for solving the problems are described, while later the
implementation is discussed. In this work, we used the software FEniCS [27]
with the extension of hIPPYlib [25] which we extended further for the needs of
our all-at-once formulations. In this work the method for discretizing and
solving the inverse problem is based on the methods described in [6]. However,
we change the computation of the Hessian including all the prior information
we have in the all-at-once setting. The equation describing the relation
between the data is given via

$\mathbb{O}u=y+\eta,$ (119)

where $\mathbb{O}:U\to Y$ is the operator which maps the state variable into
the observation space $Y$ and $\eta$ is the random variable describing the
noise in this equation. We assume that $\eta$ is random white noise.

###  6.1 Lagrangian method for computing the adjoint based Hessian and
gradient.

In this section we compute the adjoint based Hessian and gradient for a cost
functional formulated in the all-at-once setting. This functional is then
analyzed twice, first for the parameter and then also for the state variable.

For the computation of the Hessian with respect to the parameter we consider
the following minimization problem

$\displaystyle\min_{\theta}$ $\displaystyle
J_{\alpha}(\theta)=\|\mathbb{O}u-y\|_{Y}^{2}+\frac{\alpha}{2}\mathcal{R}(u,\theta)$
(120) $\displaystyle s.t.~{}~{}$ $\displaystyle\mathbb{M}(u,\theta)=0$ (121)

with the regularization given by the prior distribution and written as

$\mathcal{R}(u,\theta)=\langle u,C_{1}u\rangle_{U}+2\langle
u,C_{2}\theta\rangle_{U}+\langle\theta,C_{3}\theta\rangle_{X},$ (122)

with selfadjoint operators $C_{1}:U\to U$, $C_{3}:X\to X$. Together, this
leads to the following Lagrange functional of $u\in U$, $\theta\in X$, $p\in
W$,

$\displaystyle\mathcal{L}_{\alpha}(u,\theta,p):=$
$\displaystyle\frac{1}{2}\|\mathbb{O}u-y\|_{Y}^{2}+\frac{\alpha}{2}\big{(}\langle
u,C_{1}u\rangle_{U}+2\langle
u,C_{2}\theta\rangle_{U}+\langle\theta,C_{3}\theta\rangle_{X}\big{)}$ (123)
$\displaystyle+\mathbb{M}(u,\theta)p,$ (124)

which is used in order to derive adjoint based methods for gradient and
Hessian computation. We start with computing the first derivatives and for
this purpose use the Riesz isomorphisms $I_{R}^{U}:U\to U^{\prime}$,
$u\mapsto(v\mapsto\langle u,v\rangle_{U})$, $I_{R}^{X}:X\to X^{\prime}$

$\displaystyle\frac{\partial\mathcal{L}_{\alpha}}{\partial p}(u,\theta,p)$
$\displaystyle=\mathbb{M}(u,\theta)\in W^{\prime}$ (125)
$\displaystyle\frac{\partial\mathcal{L}_{\alpha}}{\partial u}(u,\theta,p)$
$\displaystyle=I_{R}^{U}\Bigl{(}\mathbb{O}^{*}(\mathbb{O}u-y)+\alpha\left(C_{1}u+C_{2}\theta\right)\Bigl{)}+\frac{\partial\mathbb{M}}{\partial
u}(u,\theta)^{\star}p\in U^{\prime}$ (126)
$\displaystyle\frac{\partial\mathcal{L}_{\alpha}}{\partial\theta}(u,\theta,p)$
$\displaystyle=\alpha
I_{R}^{X}\left(C_{2}^{*}u+C_{3}\theta\right)+\frac{\partial\mathbb{M}}{\partial\theta}(u,\theta)^{\star}p\in
X^{\prime}\,,$ (127)

where the derivative with respect to $p$ set to zero gives the forward
problem. The following operators are now used as abbreviations of the
derivatives of $\mathbb{M}(u,\theta)$

$-\frac{\partial\mathbb{M}}{\partial u}(u,\theta)=:K:U\to
W^{\prime},\qquad\frac{\partial\mathbb{M}}{\partial\theta}(u,\theta)=:L:X\to
W^{\prime},$ (128)

along with their Banach space adjoints $K^{\star}:W\to U^{\prime}$,
$L^{\star}:W\to X^{\prime}$ and their liftings to $U$ and $X$, respectively
$K^{*}:=(I_{R}^{U})^{-1}K^{\star}:W\to U$,
$L^{*}:=(I_{R}^{X})^{-1}L^{\star}:W\to X$. Since we deal with linear models
here, $K$ and $L$ are in fact independent of $u$ and $\theta$ and we can write
$\mathbb{M}(u,\theta)=-Ku+L\theta+f$. Then, the derivative with respect to $u$
when set to zero leads to the adjoint problem with the adjoint linearization
of the model defining the differential operator and the residual plus some
additional term from the prior define the right hand side

$K^{*}p=\mathbb{O}^{*}(\mathbb{O}u-y)+\alpha(C_{1}u+C_{2}\theta)\,.$ (129)

With $p=p(\theta)$ satisfying (129), and $u=u(\theta)$ (also inserted into
(129)) satisfying the state equation (121), the derivative with respect to
$\theta$ gives us the gradient of the reduced cost function
$j_{\alpha}(\theta)=J_{\alpha}(u(\theta),\theta))$, which results in

$\displaystyle j_{\alpha}^{\prime}(\theta)=$
$\displaystyle\frac{d}{d\theta}J_{\alpha}(u(\theta),\theta))=\frac{d}{d\theta}\mathcal{L}_{\alpha}(u(\theta),\theta,p(\theta)))$
$\displaystyle=\frac{\partial\mathcal{L}_{\alpha}}{\partial\theta}(u(\theta),\theta,p(\theta))=\alpha
I_{R}^{X}\left(C_{2}^{*}u+C_{3}\theta\right)+L^{\star}p(\theta)\,.$

Now we can write down a new Lagrange functional for the computation of the
Hessian matrix using the concept of second order adjoints. This new functional
denoted by $\mathcal{L}^{\mathcal{H}}$ is given by

$\displaystyle\mathcal{L}^{\mathcal{H}}(u,\theta,p;\hat{u},\hat{\theta},\hat{p}):=\frac{\partial\mathcal{L}_{\alpha}}{\partial
p}(u,\theta,p)\hat{p}+\frac{\partial\mathcal{L}_{\alpha}}{\partial
u}(u,\theta,p)\hat{u}+\frac{\partial\mathcal{L}_{\alpha}}{\partial\theta}(u,\theta,p)\hat{\theta}$
$\displaystyle=\mathbb{M}(u,\theta)\hat{p}+\langle\mathbb{O}u-y,\mathbb{O}\hat{u}\rangle_{Y}+\alpha\langle
C_{1}u+C_{2}\theta,\hat{u}\rangle_{U}-\langle
K\hat{u},p\rangle_{W^{\prime},W}$ $\displaystyle\qquad+\alpha\langle
C_{2}^{*}u+C_{3}\theta,\hat{\theta}\rangle_{X}+\langle
L\hat{\theta},p\rangle_{W^{\prime},W}.$

Again, the derivatives are computed and finally, the Hessian of the reduced
cost functional can be given.

$\displaystyle j_{\alpha}^{\prime\prime}(\theta)(\hat{\theta},h)=$
$\displaystyle\frac{d^{2}}{d\theta^{2}}J_{\alpha}(u(\theta),\theta))(\hat{\theta},h)=\frac{d}{d\theta}\mathcal{L}^{H}(u(\theta),\theta,p(\theta);\hat{u},\hat{\theta},\hat{p})h$
$\displaystyle=$ $\displaystyle\frac{\partial\mathcal{L}^{H}}{\partial
u}(u(\theta),\theta,p(\theta);\hat{u},\hat{\theta},\hat{p})\frac{\partial
u}{\partial\theta}(\theta)h+\frac{\partial\mathcal{L}^{H}}{\partial\theta}(u(\theta),\theta,p(\theta);\hat{u},\hat{\theta},\hat{p})h$
$\displaystyle+\frac{\partial\mathcal{L}^{H}}{\partial
p}(u(\theta),\theta,p(\theta);\hat{u},\hat{\theta},\hat{p})\frac{\partial
p}{\partial\theta}(\theta)h$

Choosing $\hat{u}=\hat{u}(\theta,\hat{\theta})$,
$\hat{p}=\hat{p}(\theta,\hat{\theta})$ such that the first and the last term
vanish for all $h$, i.e.,

$\displaystyle 0$ $\displaystyle=\frac{\partial\mathcal{L}^{H}}{\partial
u}(u(\theta),\theta,p(\theta);\hat{u},\hat{\theta},\hat{p})$
$\displaystyle=I_{R}^{U}\left(-K^{*}\hat{p}(\theta,\hat{\theta})+\mathbb{O}^{*}\mathbb{O}\hat{u}(\theta,\hat{\theta})+\alpha(C_{1}^{*}\hat{u}(\theta,\hat{\theta})+C_{2}\hat{\theta})\right),$
$\displaystyle 0$ $\displaystyle=\frac{\partial\mathcal{L}^{H}}{\partial
p}(u(\theta),\theta,p(\theta);\hat{u},\hat{\theta},\hat{p})=-K\hat{u}(\theta,\hat{\theta})+L\hat{\theta}=0,$

where the lower equation is solved first and the resulting $\hat{u}(\theta)$
is inserted into the upper equation, which is the resolved for
$\hat{p}(\theta)$. This together with

$\frac{\partial\mathcal{L}^{H}}{\partial\theta}(u,\theta,p;\hat{u},\hat{\theta},\hat{p})h=\langle
Lh,\hat{p}\rangle_{W^{\prime},W}+\alpha\langle
C_{2}h,\hat{u}\rangle_{U}+\alpha\langle C_{3}h,\hat{\theta}\rangle_{U}$

yields, for any $h\in X$,

$j_{\alpha}^{\prime\prime}(\theta)(\hat{\theta},h)=\langle
L^{*}\hat{p}(\theta,\hat{\theta})+\alpha\left(C_{2}^{*}\hat{u}(\theta,\hat{\theta})+C_{3}\hat{\theta}\right),h\rangle_{X},$

thus, abbreviating $K^{-*}=(K^{*})^{-1}$,

$lj_{\alpha}^{\prime\prime}(\theta)(\hat{\theta},\check{\theta})=\langle\left(L^{*}K^{-*}[\mathbb{O}^{*}\mathbb{O}+\alpha
C_{1}]K^{-1}L+\alpha C_{3}+\alpha L^{*}K^{-*}C_{2}+\alpha
C_{2}^{*}K^{-1}L\right)\hat{\theta},\check{\theta}\rangle_{X}.$

Thus, altogether the gradient and Hessian of the reduced cost function $j$
read as follows

$\displaystyle j_{\alpha}^{\prime}(\theta)=\alpha
I_{R}^{X}\left(C_{2}^{*}K^{-1}L+C_{3}\right)\theta+L^{\star}K^{-*}\bigl{[}\mathbb{O}^{*}(\mathbb{O}K^{-1}L\theta-y)+\alpha(C_{1}K^{-1}L+C_{2})\theta\bigr{]}\,,$
$\displaystyle
j_{\alpha}^{\prime\prime}(\theta)=I_{R}^{X}H_{\theta}=I_{R}^{X}\Bigl{(}L^{*}K^{-*}[\mathbb{O}^{*}\mathbb{O}+\alpha
C_{1}]K^{-1}L+\alpha C_{3}+\alpha L^{*}K^{-*}C_{2}+\alpha
C_{2}^{*}K^{-1}L\Bigr{)}.$

If $L$ is invertible, then the role of $u$ and $\theta$ in the above procedure
can be exchanged and for the reduced cost function
$\ell_{\alpha}(u)=J_{\alpha}(u,\theta(u))$ with $\theta(u)=L^{-1}(Ku-f)$ and
$p(u)=-\alpha L^{-*}(C_{2}^{*}u+C_{3}\theta(u))$ we get

$\displaystyle\ell_{\alpha}^{\prime}(u)=\frac{d}{du}J_{\alpha}(u,\theta(u)))=\frac{d}{du}\mathcal{L}_{\alpha}(u,\theta(u),p(u))=\frac{\partial\mathcal{L}_{\alpha}}{\partial
u}(u,\theta(u),p(u))$
$\displaystyle=I_{R}^{U}\Bigl{(}\mathbb{O}^{*}(\mathbb{O}u-y)+\alpha\Bigl{(}C_{1}u+C_{2}L^{-1}(Ku-f)+K^{*}L^{-*}(C_{2}^{*}u+C_{3}L^{-1}(Ku-f))\Bigr{)}\Bigr{)}$

and, with $\hat{\theta}(u,\hat{u})=L^{-1}K\hat{u}$,
$\hat{p}(u,\hat{u})=-\alpha L^{-*}(C_{2}^{*}+C_{3}L^{-1}K)\hat{u}$,

$\displaystyle\ell_{\alpha}^{\prime\prime}(u)(\hat{u},\check{u})$
$\displaystyle=\frac{d}{du}\mathcal{L}^{H}(u,\theta(u),p(u);\hat{u},\hat{\theta}(u,\hat{u}),\hat{p}(u,\hat{u}))\check{u}$
$\displaystyle=\frac{\partial\mathcal{L}^{H}}{\partial
u}(u,\theta(u),p(u);\hat{u},\hat{\theta}(u,\hat{u}),\hat{p}(u,\hat{u}))\check{u}$
$\displaystyle=\langle-K^{*}\hat{p}(u,\hat{u})+\mathbb{O}^{*}\mathbb{O}\hat{u}+\alpha(C_{1}\hat{u}+C_{2}\hat{\theta}(u,\hat{u})),\check{u}\rangle_{U}$
$\displaystyle=\langle\left[\mathbb{O}^{*}\mathbb{O}+\alpha(C_{1}+C_{2}L^{-1}K+K^{*}L^{-*}C_{2}^{*}+K^{*}L^{-*}C_{3}L^{-1}K)\right]\hat{u},\check{u}\rangle_{U}$

that is, $\ell_{\alpha}$ is the quadratic functional
$\ell_{\alpha}(u)=\frac{1}{2}\langle Hu,u\rangle+\langle g,u\rangle$ with

$\displaystyle
H=(I_{R}^{U})^{-1}\ell_{\alpha}^{\prime\prime}(u)=\mathbb{O}^{*}\mathbb{O}+\alpha(C_{1}+C_{2}L^{-1}K+K^{*}L^{-*}C_{2}^{*}+K^{*}L^{-*}C_{3}L^{-1}K)\,,$
$\displaystyle
g=(I_{R}^{U})^{-1}\ell_{\alpha}^{\prime}(0)=-\mathbb{O}^{*}y-\alpha(C_{2}+K^{*}L^{-*}C_{3})L^{-1}f.$

The computed Hessians can now be used to evaluate the MAP estimator for the
posterior mean and the posterior covariance.

####  6.1.1 Inverse source problem.

$U=H_{0}^{1}(\Omega)\cap H^{2}(\Omega)$, $X=Y=W=W^{\prime}=L^{2}(\Omega)$,
$\mathbb{M}(u,\theta)=-\mathcal{A}u+\theta$, $\mathbb{O}u=u$,
$K=\mathcal{A}:H_{0}^{1}(\Omega)\cap H^{2}(\Omega)\to L^{2}(\Omega)$,
$K^{*}=\mathcal{A}^{-1}:L^{2}(\Omega)\to H_{0}^{1}(\Omega)\cap H^{2}(\Omega)$,
$L=\mbox{id}:L^{2}(\Omega)\to L^{2}(\Omega)$,
$L^{*}=\mbox{id}:L^{2}(\Omega)\to L^{2}(\Omega)$,
$\mathbb{O}:H_{0}^{1}(\Omega)\cap H^{2}(\Omega)\to L^{2}(\Omega)$,
$\mathbb{O}^{*}=\mathcal{A}^{-2}:L^{2}(\Omega)\to H_{0}^{1}(\Omega)\cap
H^{2}(\Omega)$, which results in

$\displaystyle g_{\theta}(\theta)$
$\displaystyle:=(I_{R}^{X})^{-1}j_{\alpha}^{\prime}(\theta)$
$\displaystyle=\alpha(C_{2}^{*}\mathcal{A}^{-1}+C_{3})\theta+\mathcal{A}\bigl{[}\mathcal{A}^{-2}(\mathcal{A}^{-1}\theta-y)+\alpha(C_{1}\mathcal{A}^{-1}+C_{2})\theta\bigr{]},$
$\displaystyle H_{\theta}$
$\displaystyle:=(I_{R}^{X})^{-1}j_{\alpha}^{\prime\prime}=\mathcal{A}(\mathcal{A}^{-2}+\alpha
C_{1})\mathcal{A}^{-1}+\alpha C_{3}+\alpha\mathcal{A}C_{2}+\alpha
C_{2}^{*}\mathcal{A}^{-1},$

for the parameter $\theta$, where the computation of the adjoint is as
described in Section 2. For the state $u$ we have

$\displaystyle g_{u}(u)$
$\displaystyle:=(I_{R}^{U})^{-1}\ell_{\alpha}^{\prime}(u)$
$\displaystyle=\mathcal{A}^{-2}(u-y)+\alpha\Bigl{(}C_{1}u+C_{2}(\mathcal{A}u-f)+\mathcal{A}^{-1}(C_{2}^{*}u+C_{3}(\mathcal{A}u-f)),$
$\displaystyle H_{u}$
$\displaystyle:=(I_{R}^{U})^{-1}\ell_{\alpha}^{\prime\prime}=\mathcal{A}^{-2}+\alpha(C_{1}+C_{2}\mathcal{A}+\mathcal{A}^{-1}C_{2}^{*}+\mathcal{A}^{-1}C_{3}\mathcal{A})).$

####  6.1.2 Backwards heat problem.

$U=U_{0}=\\{w\in L^{2}(H_{0}^{1}(\Omega))\cap
H^{1}(H^{-1}(\Omega)):w(0,x)=0\\}$, $W=L^{2}(H_{0}^{1}(\Omega))$,
$W^{\prime}=L^{2}(H^{-1}(\Omega))$, $X=H_{0}^{1}(\Omega)$, $Y=L^{2}(\Omega)$,
$\mathbb{M}(u,\theta)=(\partial_{t}+\mathcal{A})u+\underline{I}\mathcal{A}\theta$,
$\mathbb{O}=\delta_{T}$, $K=(\partial_{t}+\mathcal{A}):U_{0}\to W^{\prime}$,
$K^{-1}=\int_{0}^{.}e^{-\mathcal{A}(.-s)}\cdot(s)\mbox{ d}s:W^{\prime}\to
U_{0}$, $K^{*}=\int_{0}^{.}e^{-\mathcal{A}(.-s)}\cdot(s)\mbox{
d}s:W^{\prime}\to U_{0}$, $K^{-*}=K$,
$L=\underline{I}\mathcal{A}:H_{0}^{1}(\Omega)\to W^{\prime}$,
$L^{*}=\mathcal{A}^{-1}\int_{0}^{T}\cdot(s)\mbox{ d}s:W^{\prime}\to
H_{0}^{1}(\Omega)$,
$K^{-1}L=\underline{I}-e^{\mathcal{A}\cdot}:H_{0}^{1}(\Omega)\to U_{0}$,
$L^{*}K^{-*}=\mathcal{A}^{-1}\delta_{T}-\underline{I}^{*}:U_{0}\to
H_{0}^{1}(\Omega)$, $\mathbb{O}=\delta_{T}:U_{0}\to L^{2}(\Omega)$,
$\mathbb{O}^{*}=\frac{1}{2}(e^{-\mathcal{A}(T-.)}-e^{-\mathcal{A}(T+.)}):L^{2}(\Omega)\to
U_{0}$, resulting in the following instances

$\displaystyle
g_{\theta}(\theta):=(I_{R}^{X})^{-1}j_{\alpha}^{\prime}(\theta)=\alpha(C_{2}^{*}(\underline{I}-e^{\mathcal{A}\cdot})+C_{3})\theta+(\mathcal{A}^{-1}\delta_{T}-\underline{I}^{*})$
$\displaystyle\qquad\cdot\Bigl{[}{\textstyle\frac{1}{2}}(e^{-\mathcal{A}(T-.)}-e^{-\mathcal{A}(T+.)})(\delta_{T}[\underline{I}-e^{\mathcal{A}\cdot}]\theta-y)+\alpha(C_{1}[\underline{I}-e^{\mathcal{A}\cdot}]+C_{2})\theta\Bigr{]},$
$\displaystyle
H_{\theta}:=(I_{R}^{X})^{-1}j_{\alpha}^{\prime\prime}=(\mathcal{A}^{-1}\delta_{T}-\underline{I}^{*}){\textstyle\frac{1}{2}}(e^{-\mathcal{A}(T-.)}-e^{-\mathcal{A}(T+.)})(\delta_{T}[\underline{I}-e^{\mathcal{A}\cdot}]$
$\displaystyle\qquad+\alpha\Bigl{(}(\mathcal{A}^{-1}\delta_{T}-\underline{I}^{*})C_{1}[\underline{I}-e^{\mathcal{A}\cdot}]+(\mathcal{A}^{-1}\delta_{T}-\underline{I}^{*})C_{2}+C_{2}^{*}(\underline{I}-e^{\mathcal{A}\cdot})+C_{3},\Bigr{)}$

for the parameter.

Since $L$ is not invertible here, we substitute it by its right inverse
$L^{\sharp}=\frac{1}{T}L^{*}$ to obtain, for the state,

$\displaystyle
g_{u}(u):=(I_{R}^{U})^{-1}\ell_{\alpha}^{\prime}(u)={\textstyle\frac{1}{2}}(e^{-\mathcal{A}(T-.)}-e^{-\mathcal{A}(T+.)})(\delta_{T}u-y)$
$\displaystyle\qquad+\alpha\Bigl{(}C_{1}u+C_{2}\mathcal{A}^{-1}{\textstyle\frac{1}{T}}{\textstyle\int}_{0}^{T}((\partial_{t}+\mathcal{A})u-f)(s)\mbox{
d}s$
$\displaystyle\qquad+{\textstyle\frac{1}{T}}(\mbox{id}-e^{-\mathcal{A}t})\Bigl{(}C_{2}^{*}u+C_{3}{\textstyle\frac{1}{T}}{\textstyle\int}_{0}^{T}((\partial_{t}+\mathcal{A})u-f)(s)\mbox{
d}s\Bigr{)}\Bigr{)},$ $\displaystyle
H_{u}:=(I_{R}^{U})^{-1}\ell_{\alpha}^{\prime\prime}={\textstyle\frac{1}{2}}(e^{-\mathcal{A}(T-.)}-e^{-\mathcal{A}(T+.)})\delta_{T}$
$\displaystyle\qquad+\alpha\Bigl{(}C_{1}+C_{2}\mathcal{A}^{-1}{\textstyle\frac{1}{T}}{\textstyle\int}_{0}^{T}(\partial_{t}+\mathcal{A})\cdot(s)\mbox{
d}s$
$\displaystyle\qquad+{\textstyle\frac{1}{T}}(\mbox{id}-e^{-\mathcal{A}t})\Bigl{(}C_{2}^{*}+C_{3}{\textstyle\frac{1}{T}}{\textstyle\int}_{0}^{T}(\partial_{t}+\mathcal{A})\cdot(s)\mbox{
d}s\Bigr{)}.$

###  6.2 Implementation.

As already mentioned in the beginning of this section the analysis is done in
python, with the help of Fenics and hippylib. With hippylib, the machinery of
Bayesian inversion can be applied almost automatically to all stationary
problems. However, concerning time-dependent problems a bit more work is
necessary to get the results which are needed. In our case this means that the
backwards heat problem had to be formulated in hippylib fashion to hand over
to the algortihms for reconstruction. These algorithms were adapted to the
all-at-once setting. In the following two sections the implementation of the
problems as well as the results of the reconstruction are described in detail.

From now on, we fix the spatial domain $\Omega$ on which both inverse problems
are considered to be the unit square and space discretization is done by means
of Lagrange finite elements. For time discretization in the backwards heat
problem we use zero order discontinuos Galerkin method as already described
previously, which leads to a backward Euler time stepping scheme. The
discretized versions of paramter and state will be denoted by $u_{h}$
$\theta_{h}$, the corresponding finite dimensional function spaces by $U_{h}$
and $X_{h}$.

To generate the synthetic observations the forward model is solved with the
true parameter. Here, the dimension of the spatial discretization space for
the state and the parameter is both 1681. In order to avoid an inverse crime,
this forward simulation is done on a finer grid as compared to the
reconstruction, where the dimension is chosen to be 961. Discrete observation
points are generated as 100 random points on the computational mesh, and the
data at these points is perturbed with random noise. The measurement noise is
assumed to follow a Gaussian distribution with zero mean and covariance
$\Sigma=\delta^{2}I$, where $\delta=0.01$, and $I\in\mathbb{R}^{100\times
100}$. These synthetic observations are used (together with the prior) as
input for the reconstructions.

####  6.2.1 Inverse source problem.

The first test example considered for the implementation of the all-at-once
version of Bayesian reconstruction is the inverse source problem. For this
experiment the source is given by the function visualized in Figure 3(a). As
described above, the Poisson equation with this source is solved and random
discrete observations are constructed from the resulting state, see Figure
3(b) and 3(c). To reconstruct the source these observations are used together
with the prior distribution. In this case the prior distribution is defined by

$\mu_{\mbox{prior}}=\mathcal{N}(\mathbf{m},\mathcal{C}_{0}),\mbox{ with
}\mathbf{m}=\left(\begin{array}[]{c}m_{u}\\\
m_{\theta}\end{array}\right),\mathcal{C}_{0}=\left(\begin{array}[]{cc}\kappa_{s}M+\gamma_{s}K&0\\\
0&\kappa_{p}M+\gamma_{p}K\end{array}\right)^{-1},$ (130)

with $M$ and $K$ denoting the mass and the stiffness matrix, respectively, cf.
(46), (47), and chosen parameters $m_{\theta}=0,\kappa_{p}=\kappa_{s}=10^{-2}$
and $\gamma_{p}=\gamma_{s}=35$. The prior matrix is chosen to be block
diagonal, as operators on the off diagonal have not improved the
reconstruction for this problem. The mean $m_{\theta}$ is chosen to be zero
and the mean $m_{u}$ is the output of the forward problem when the source is
$m_{\theta}$. The reconstruction using the synthetic observations and prior
with the help of the CG-solver for the Hessian shown in the above computations
leads to the results displayed in Figure 4(a) for the parameter and Figure
4(b) for the state. The posterior distribution is visualized with the help of
samples for both the parameter in Figure 5 and the state in Figure 6. Let us
comment on the choice of the parameters $\kappa$ and $\gamma$ in the prior
covariance. The reconstruction of the parameter appears to be strongly
influenced by the choice of $\kappa_{p}$, $\gamma_{p}$ for the parameter
prior, while the choice of $\kappa_{s}$, $\gamma_{s}$ in the state prior do
not really influence the reconstruction neither for the parameter nor for the
state. Some comparison of different choices of parameters can be seen in
Figure 7 for the parameter and Figure 8. Nonzero choices for the covariance
operators $C_{2}$ or $C_{2}^{*}$ other than zero have not improved the
reconstructions for this example.

In Figures 4 – 8, the noise was set to 1% as described in the introduction of
this section. However, the reconstructions are quite good up to a noise level
of 3%, as it can be seen in Figure 9 and Figure 10.

(a) True source.

(b) True state.

(c) Discrete observations.

Figure 3: Simulated observations for the inverse source problem.

(a) Reconstructed source.

(b) Reconstructed state.

Figure 4: Reconstructions for the inverse source problem.

Figure 5: Samples from the posterior for the parameter variable of the inverse
source problem.

Figure 6: Samples from the posterior for the state variable of the inverse
source problem.

(a) Source: $\gamma_{p}=1$ and $\kappa_{p}=10^{-2}$.

(b) State: $\gamma_{p}=1$ and $\kappa_{p}=10^{-2}$.

(c) Source: $\gamma_{p}=10^{-2}$ and $\kappa_{p}=70$.

(d) State: $\gamma_{p}=10^{-2}$ and $\kappa_{p}=70$.

Figure 7: Reconstructions for the inverse source problem for different values
for $\kappa_{p}$ and $\gamma_{p}$ and $\kappa_{s}=10^{-2}$ and
$\gamma_{s}=35$.

(a) Source: $\gamma_{s}=10^{-2}$ and $\kappa_{s}=70$.

(b) State: $\gamma_{s}=10^{-2}$ and $\kappa_{s}=70$.

(c) Source: $\gamma_{s}=10^{-5}$ and $\kappa_{s}=100$.

(d) State: $\gamma_{s}=10^{-5}$ and $\kappa_{s}=100$.

Figure 8: Reconstructions for the inverse source problem for different values
for $\kappa_{s}$ and $\gamma_{s}$ and $\kappa_{p}=10^{-2}$ and
$\gamma_{p}=35$.

Figure 9: Reconstructions for the inverse source problem with different noise
levels.

Figure 10: Reconstructions for the inverse source problem with different noise
levels.

####  6.2.2 Backwards heat equation, sampled initial condition.

For the backwards heat equation the prior as described in Section 5 is used,
with the covariance operator for the initial condition chosen to be

$C_{p}=(\kappa M+\gamma K)^{-1},$ (131)

where $M$ and $K$ are the mass and stiffness matrices resulting from the
finite element discretization, with $\kappa=1.5$ and $\gamma=0.5$. The time
variable is discretized in the same way as in Section 3 resulting in the
covariance operator

$C_{s}=\mbox{diag}((e^{-t_{i}K}C_{p}e^{-t_{i}K})_{i=0}^{N})$ (132)

for the state variable, so that altogether we arrive at

$\mathcal{C}_{0}=\left(\begin{array}[]{cc}C_{s}&0\\\
0&C_{p}\end{array}\right).$

Also here, nonzero covariance operators $C_{2}$ or $C_{2}^{*}$ did not improve
the reconstructions but would only complicate computation of powers of
$\mathcal{C}_{0}$.

Here, we have the setting $T=0.1,N=4$ and $t_{i}=i\frac{T}{N}$. The prior mean
for the inital condition is chosen to be zero, therefore, also the mean for
the state is zero for each $t\in[0,T]$. Due to the diagonal structure of the
operator matrix it is not necessary to save the whole matrix, instead for
every time step only the application of the matrix vector product is needed.
This overcomes the problem of block matrix usage in python. For the
experiment, the initial condition is sampled from the prior distribution as
seen in Figure 11(b). Then, the forward problem is solved (Figure 11(b) and
discrete synthetic observations of the state at time $T=0.1$ are constructed,
see Figure 11(c). Finally, the MAP estimators can be computed as described in
the beginning of this section for both the parameter and the state variable,
see Figure 12. In this example the posterior samples are only shown for the
parameter in Figure 13.

(a) True initial condition.

(b) True state at time $T$.

(c) Discrete observations.

Figure 11: Simulated observations for the backwards heat problem.

Figure 12: Reconstructions for the backwards heat problem.

Figure 13: Samples from the posterior for the initial condition of the
backwards heat problem.

####  6.2.3 Backwards heat equation, chosen initial condition.

We now demonstrate performance of the method with a fixed initial condition as
visualized in Figure 14(a). Again, discrete synthetic observations are made
which are used for the reconstruction, see Figure 14(c). The prior is chosen
as in Subsection 6.2.2. The reconstructions of the parameter and the state are
shown in Figure 15. Samples are shown in Figure 16.

(a) True initial condition.

(b) True state at time $T$.

(c) Discrete observations.

Figure 14: Simulated observations for the backwards heat problem.

Figure 15: Reconstructions for the backwards heat problem.

Figure 16: Samples from the posterior for the initial condition of the
backwards heat problem.

####  6.2.4 Backwards heat equation, chosen initial condition, prior motivated
by the link condition.

Referring to Section 5.2 we do the same experiments as in the last section
with the prior given by $\tilde{\mathcal{C}}$ as in (104) but skipping the
off-diagonal blocks for ease of implementation

$\mu_{\mbox{prior}}=\mathcal{N}(\mathbf{m},\mathcal{C}_{0}),\mbox{ with
}\mathbf{m}=\left(\begin{array}[]{c}m_{u}\\\
m_{\theta}\end{array}\right),\quad\mathcal{C}_{0}=\left(\begin{array}[]{cc}\mathcal{M}&0\\\
0&K^{1/2}\end{array}\right).$ (133)

The operator $\mathcal{M}$ denotes an operator consisting of the mass matrix
$M$ for every timestep $t\in[0,T]$. The setup of the numerical experiment is
the same as in the last example concerning discretization, initial condition
and observations. The reconstructions can be seen in Figure 17, which suggests
quite a good fit also for this prior.

Figure 17: Reconstructions for the backwards heat problem with the heuristic
prior.

##  7 Conclusions and Remarks.

In this paper we have combined the Bayesian approach with an all-at-once
formulation of the inverse problem. We have done so for linear problems, in
particular focusing on two prototypical examples, namely the inverse source
problem for the Poisson equation and the backwards heat equation.
Our next step will be to extend this approach to nonlinear problems such as
the identification of coefficients in time dependent and stationary PDEs.
Moreover we will further investigate the use of joint state and parameter
priors, possibly also taking covariances between them into account.

## Acknowledgment

This work was supported by the Austrian Science Fund FWF under the grants
P30054 and DOC 78.

## References

  * [1] Alen Alexanderian, Noemi Petra, Georg Stadler, and Omar Ghattas.  Optimal design of experiments for infinite-dimensional bayesian linear inverse problems with regularized l-sparsification.  SIAM Journal on Scientific Computing, 36(5):A2122–A2148, 2014.
  * [2] A.W. van der Vaart B.T. Knapik and J.H. van Zanten.  Bayesian recovery of the initial condition for the heat equation.  Communications in Statistics, Theory and Methods, 42, 2013.
  * [3] S. Gugushvili, A. van der Vaart, and D. Yan.  Bayesian inverse problems with partial observations.  ArXiv e-prints, February 2018.
  * [4] H.P.Flath, L.C. Wilcox, V. Akcelik, J.Hill, B. Van Bloemen Waanders, and O. Ghattas.  Fast algorithms for bayesian uncertainty quantification in large-scale linear inverse problems based on low-rank partial hessian approximations.  SIAM J. Scientific Computing, 33, 2011.
  * [5] James Martin, Lucas C. Wilcox, Carsten Burstedde, and Omar Ghattas.  A stochastic newton mcmc method for large-scale statistical inverse problems with application to seismic inversion.  SIAM Journal on Scientific Computing, 34(3):A1460–A1487, 2012\.
  * [6] James Martin Tan Bui-Thanh, Omar Ghattas and Georg Stadler.  A computational framework for infinite-dimensional bayesian inverse problems part i: The linearized case, with application to global seismic inversion.  SIAM J. Scientific Computing, 35, 2013.
  * [7] S. L. Cotter, M. Dashti, and A. M. Stuart.  Approximation of bayesian inverse problems for pdes.  SIAM Journal on Numerical Analysis, 48(1):322–345, 2010.
  * [8] Masoumeh Dashti and Andrew M. Stuart.  The bayesian approach to inverse problems.  In Roger Ghanem, David Higdon, and Houman Owhadi, editors, Handbook of Uncertainty Quantification, pages 311–428. Springer International Publishing, Cham, 2017.
  * [9] Sergios Agapiou and Peter Mathé.  Posterior contraction in bayesian inverse problems under gaussian priors.  In Bernd Hofmann, Antonio Leitão, and Jorge P. Zubelli, editors, New Trends in Parameter Identification for Mathematical Models, pages 1–29. Springer International Publishing, Cham, 2018.
  * [10] Peter Mathé.  Bayesian inverse problems with non-commuting operators.  Mathematics of Computation, 01 2018.
  * [11] M. Burger and W. Mühlhuber.  Iterative regularization of parameter identification problems by sequential quadratic programming methods.  Inverse Problems, 18:943–969, 2002.
  * [12] M. Burger and W. Mühlhuber.  Numerical approximation of an sqp-type method for parameter identification.  SIAM J. Numer. Anal., 40:1775–1797, 2002.
  * [13] E Haber and U M Ascher.  Preconditioned all-at-once methods for large, sparse parameter estimation problems.  Inverse Problems, 17(6):1847–1864, nov 2001.
  * [14] B Kaltenbacher.  All-at-once versus reduced iterative methods for time dependent inverse problems.  Inverse Problems, 33(6):064002, 2017.
  * [15] B. Kaltenbacher.  Regularization based on all-at-once formulations for inverse problems.  SIAM Journal on Optimization, 28:620–645, 2018.
  * [16] B Kaltenbacher, A. Kirchner, and B. Vexler.  Goal oriented adaptivity in the irgnm for parameter identification in pdes ii: all-at once formulations.  Inverse Problems, 30:045002, 2014.
  * [17] Tram Thi Ngoc Nguyen.  Landweber–kaczmarz for parameter identification in time-dependent inverse problems: all-at-once versus reduced version.  Inverse Problems, 35(3):035009, feb 2019.
  * [18] B. Kaltenbacher, M. Kaltenbacher, and S. Gombots.  Inverse scheme for acoustic source localization using microphone measurements and finite element simulations.  Acta Acustica united with Acustica, pages 647–656, 2018.
  * [19] A. Schuhmacher, K. Rasmussen, and C. Hansen.  Sound source reconstruction using inverse boundary element calculations.  J. Acoust. Soc. Am., 113:114–127, 2003.
  * [20] Volkan Akcelik, George Biros, Andrei Draganescu, Omar Ghattas, Judith Hill, and Bart G. van Bloemen Waanders.  Inversion of airborne contaminants in a regional model.  In Computational Science - ICCS 2006, 6th International Conference, Reading, UK, May 28-31, 2006, Proceedings, Part III, pages 481–488, 2006.
  * [21] Richard Kowar and Otmar Scherzer.  Attenuation models in photoacoustics.  In H. Ammari, editor, Mathematical Modeling in Biomedical Imaging II: Optical, Ultrasound, and Opto-Acoustic Tomographies, volume 2035 of Lecture Notes in Mathematics, pages 85–130. Springer Verlag, Berlin Heidelberg, 2012.
  * [22] K.J. Engel and R. Nagel.  One-Parameter Semigroups for Linear Evolution Equations.  Graduate Texts in Mathematics. Springer New York, 1999.
  * [23] R.G.Douglas.  On majorization, factorization, and range inclusion of operators on hilbert space.  Proc. Amer. Math. Soc., 17:413–415, 1966.
  * [24] Peter Mathé and Ulrich Tautenhahn.  Interpolation in variable hilbert scales with application to inverse problems.  Inverse Problems, 22(6):2271–2297, oct 2006.
  * [25] U. Villa, N. Petra, and O. Ghattas.  hIPPYlib: an Extensible Software Framework for Large-scale Deterministic and Bayesian Inverse Problems.  Journal of Open Source Software, 3(30), 2018.
  * [26] Lawrence C. Evans.  Partial Differential Equations, volume 19 of Graduate Studies in Mathematics.  AMS, Providence, Rhode Island, 1998.
  * [27] Christian Clason.  Numerical partial differential equations.  Lecture notes, University of Duisburg-Essen, 2017.