# Stochastic asymptotical regularization
for linear inverse problems

Ye Zhang  School of Mathematics and Statistics, Beijing Institute of
Technology, 100081 Beijing, China  Shenzhen MSU-BIT University, 518172
Shenzhen, China  Chuchu Chen111Corresponding author Academy of Mathematics and
Systems Science, Chinese Academy of Sciences, 100080 Beijing, China

###### Abstract

We introduce Stochastic Asymptotical Regularization (SAR) methods for the
uncertainty quantification of the stable approximate solution of ill-posed
linear-operator equations, which are deterministic models for numerous inverse
problems in science and engineering. We prove the regularizing properties of
SAR with regard to mean-square convergence. We also show that SAR is an
optimal-order regularization method for linear ill-posed problems provided
that the terminating time of SAR is chosen according to the smoothness of the
solution. This result is proven for both _a priori_ and _a posteriori_
stopping rules under general range-type source conditions. Furthermore, some
converse results of SAR are verified. Two iterative schemes are developed for
the numerical realization of SAR, and the convergence analyses of these two
numerical schemes are also provided. A toy example and a real-world problem of
biosensor tomography are studied to show the accuracy and the advantages of
SAR: compared with the conventional deterministic regularization approaches
for deterministic inverse problems, SAR can provide the uncertainty
quantification of the quantity of interest, which can in turn be used to
reveal and explicate the hidden information about real-world problems, usually
obscured by the incomplete mathematical modeling and the ascendence of
complex-structured noise.

##  1 Introduction

In this paper, we are interested in solving deterministic linear inverse
problems of the form

$Ax=y,$ (1)

where $A$ is a compact linear operator acting between two infinite dimensional
Hilbert spaces $\mathcal{X}$ and $\mathcal{Y}$. For simplicity, we use
$\langle\cdot,\cdot\rangle$ and $\|\cdot\|$ to denote the inner products and
norms, respectively, for both $\mathcal{X}$ and $\mathcal{Y}$. Suppose that,
instead of the exact data $y$, we are given noisy measurement
$y^{\delta}\in\mathcal{Y}$, which also obeys the deterministic noise model
with noise level $\delta>0$:

$\|y^{\delta}-y\|\leq\delta.$ (2)

In addition, we assume that $\mathcal{R}(A)$ is an infinite dimensional
subspace of $\mathcal{Y}$. Then, since $A$ is compact and ${\rm
dim}(\mathcal{R}(A))=\infty$, we have
$\mathcal{R}(A)\neq\overline{\mathcal{R}(A)}$, and the problem (1) is ill-
posed and of type II, as defined by Nashed [21]. Hence, to obtain the stable
approximate solution of problem (1)-(2), an appropriate regularization method
should be designed. The formulation (1)-(2) is a typical deterministic
mathematical model for many linear inverse problems with applications in the
natural sciences if it is supposed that unknown physical characteristics $x$
cannot be measured directly. From experiments, it is possible to obtain only
noisy data $y^{\delta}$ connected with $x$ with the help of a forward
mathematical model $A$. Some important areas in which (1) can be applied are
inverse source problems in partial differential equations and the Fredholm
integral equation of the first kind, which are the basis of fundamental
mathematical models in geophysics, imaging, and many other fields.

In classical regularization theory, the deterministic regularization methods
are usually designed for deterministic inverse problems (1). These
deterministic regularization methods can be roughly classified into two
categories: variational regularization methods and iterative regularization
methods. Tikhonov regularization is probably the most prominent variational
regularization method, while the Landweber iteration is certainly the most
popular iterative regularization approach. The development of our new
regularization approach in this paper was motivated by the latter, as, from
the viewpoint of computation, the iterative approach seems more attractive,
especially for large-scale problems. In view of this, we begin with the
Landweber iteration for (1), which is given by

$\displaystyle x^{\delta}_{k+1}=x^{\delta}_{k}+\Delta
tA^{*}(y^{\delta}-Ax^{\delta}_{k}),\quad x^{\delta}_{0}=x_{0},\quad\Delta
t\in(0,2/\|A\|^{2}),$ (3)

where $A^{*}$ denotes the adjoint operator of $A$. It is well known that the
scheme (3) is turned into regularization algorithms, i.e.
$x^{\delta}_{k^{*}}\to x^{\dagger}$ as $\delta\to 0$, by stopping the
iteration after an adequate number, $k^{*}=k^{*}(\delta)$, of steps. Here and
later on, $x^{\dagger}$ denotes the $x_{0}$-minimum norm solution of the
operator equation (1), which is assumed solvable with exact data.

The continuous analog to (3) as $\Delta t$ tends to zero is known as
asymptotic regularization or Showalter’s method (see, e.g., [25, 26]). It is
written as a first-order evolution equation of the form

$\dot{x}^{\delta}(t)=A^{*}(y^{\delta}-Ax^{\delta}(t))$ (4)

with some initial condition, where an artificial scalar time $t$ is
introduced. It has been shown that, by using Runge-Kutta integrators, all the
properties of asymptotic regularization (4) are carried over to its numerical
realization [22, 31]. This continuous framework has recently received much
attention in the field of computational mathematics. For instance, in the
field of inverse problems, the authors in [19] studied a modified evolution
equation, (4), for linear inverse problems (1) in the presence of white noise.
In particular, they connected this formulation, i.e. (4), to some classical
methods from data assimilation, namely the Kalman-Bucy filter and 3DVAR.
Moving in a different direction, the authors in [29, 8, 28, 2] extended (4) to
second-order and fractional-order gradient flows. They proved that the
developed high-order flows are accelerated optimal regularization methods,
i.e. the optimal convergence rates can be obtained with approximately the
square root of the iteration number that is needed for the first-order flow
(4). Besides the community of inverse problems, this idea has also recently
been adopted in the field of machine learning for the theoretical
understanding of machine learning: the continuous dynamical flows can be used
to recover existing machine-learning models and construct new ones (see [5,
6]). Building on this idea, in this paper we are interested in the following
stochastic differential equation, called the _Stochastic Asymptotical
Regularization_ (SAR),

$dx^{\delta}=A^{*}(y^{\delta}-Ax^{\delta})dt+f(t)dB_{t},\quad
x^{\delta}(0)=x_{0},$ (5)

where the initial point $x_{0}\in\mathcal{X}$ is non-random, the auxiliary
function $f(t)\in L^{\infty}({\mathbb{R}}_{+})$ will be defined later, and
$B_{t}$ is an ${\mathcal{X}}$-valued $Q$-Wiener process, i.e.

$B_{t}=\sum_{j=1}^{\infty}\ \sqrt{q_{j}}u_{j}\beta_{j}(t),$ (6)

where $\\{u_{j}\\}$ is the orthogonal basis of
$N(A)^{\bot}\subset{\mathcal{X}}$, $\\{q_{j}\\}$ denotes the eigenvalues of
the covariance operator $Q$ under the orthorgomal basis $\\{u_{j}\\}$, and
$\\{\beta_{j}\\}$ is a family of independent ${\mathbb{R}}$-valued Brownian
motions. Throughout this paper, we assume that $Q$ commutes with $A^{*}A$ and
${\rm tr}(Q(A^{*}A)^{-1})<\infty$. We also note that, compared with
conventional regularization methods for deterministic inverse problems (1),
SAR offers a stochastic solution, which can provide the uncertainty
quantification for the solution of the considered inverse problem.

It should be noted that the formulation (5) is different from the conventional
_Stochastic Gradient Descent_ (SGD) method [23], which employs an unbiased
estimator of the full gradient computed from one single randomly selected data
point at each iteration. Recently, SGD and its variants have been established
as the workhorses behind many challenging training tasks in deep learning [3].
The regularizing properties of SGD for linear inverse problems with an _a
priori_ and an _a posteriori_ stopping rule can be found in [13] and [11],
respectively. An extended error estimation of SGD, taking into account the
discretization levels, the decay of the step-size, and the general source
conditions, can be found in [18]. SGD for non-linear ill-posed problems was
investigated in [14].

The remainder of this paper is structured as follows: in Section 2, we perform
the convergence analysis of SAR. Section 3 is devoted to the numerical
realization and the convergence analysis of the corresponding numerical scheme
of SAR. An abstract mathematical example and a real-world problem of biosensor
tomography are presented in Section 4. Finally, concluding remarks are made in
Section 5.

##  2 Convergence analysis of SAR

This section presents some basis properties of our new regularization approach
based on the stochastic partial differential equation – SAR in (5), including
the regularization property, convergence rates under general source
conditions, some converse results, and the convergence-rate results of the
best worst-case mean-square error.

###  2.1 Preliminary results of SAR

We start with the well-posedness of the dynamical flow (5).

######  Proposition 1.

For any $f\in L^{\infty}({\mathbb{R}}_{+})$, the stochastic differential
equation (5) has a unique mild solution $x^{\delta}(t)\in\mathcal{X}$, given
by

$x^{\delta}(t)=e^{-A^{*}At}x_{0}+\int_{0}^{t}e^{-A^{*}A(t-s)}A^{*}y^{\delta}ds+\int_{0}^{t}e^{-A^{*}A(t-s)}f(s)dB_{s}.$
(7)

The random variable $x^{\delta}(t)$ is Gaussian on $\mathcal{X}$ with mean

$\mathbb{E}x^{\delta}(t)=e^{-A^{*}At}x_{0}+\int_{0}^{t}e^{-A^{*}A(t-s)}A^{*}y^{\delta}ds$
(8)

and variance operator given by

$\textup{Var}(x^{\delta}(t))=\int_{0}^{t}e^{-A^{*}A(t-s)}Qe^{-A^{*}A(t-s)}[f(s)]^{2}ds.$
(9)

###### Proof.

It is sufficient to deal with the stochastic term. From the Itô isometry, we
have

$\displaystyle{\mathbb{E}}\Big{\|}\int_{0}^{t}e^{-A^{*}A(t-s)}f(s)dB(s)\Big{\|}^{2}=$
$\displaystyle\int_{0}^{t}\|e^{-A^{*}A(t-s)}f(s)Q^{\frac{1}{2}}\|^{2}_{HS(\mathcal{X};\mathcal{X})}ds$
$\displaystyle\leq\|f\|_{\infty}^{2}\int_{0}^{t}\|e^{-A^{*}A(t-s)}Q^{\frac{1}{2}}\|^{2}_{HS(\mathcal{X};\mathcal{X})}ds,$

where $HS(\mathcal{X};\mathcal{X})$ denotes the space of Hilbert-Schmidt
operators from $\mathcal{X}$ to $\mathcal{X}$ endowed with the norm
$\|\cdot\|_{HS(\mathcal{X};\mathcal{X})}=\Big{(}\sum_{j=1}^{\infty}\|\cdot
u_{j}\|^{2}\Big{)}^{\frac{1}{2}}$. Note that

$\displaystyle\int_{0}^{t}\|e^{-A^{*}A(t-s)}Q^{\frac{1}{2}}\|^{2}_{HS(\mathcal{X};\mathcal{X})}ds$
$\displaystyle\leq{\rm
tr}(Q(A^{*}A)^{-1})\int_{0}^{t}\|(A^{*}A)^{\frac{1}{2}}e^{-A^{*}A(t-s)}\|^{2}ds$
$\displaystyle\leq{\rm
tr}(Q(A^{*}A)^{-1})\int_{0}^{t}\sup_{\lambda>0}\big{\\{}\lambda
e^{-2\lambda(t-s)}\big{\\}}ds$ $\displaystyle\leq\frac{1}{2}{\rm
tr}(Q(A^{*}A)^{-1})<\infty.$

Therefore, according to [4, Theorem 5.4], the proposition holds true. ∎

We now introduce some definitions that will be frequently used for convergence
analysis in this paper.

######  Definition 1.

A real function $\varphi:(0,\infty)\to(0,\infty)$ is called an index function
if it is continuous and strictly increasing, and satisfies the condition
$\lim_{\lambda\to 0+}\varphi(\lambda)=0$. Let $\mathcal{I}$ denote the set of
all index functions.

######  Definition 2.

An index function $\varphi$ is called $g$-subhomogeneous if there exists a
decreasing function $g$ such that $\varphi(\gamma\lambda)\leq
g(\gamma)\varphi(\lambda)$ for all $\gamma>0,\ \lambda>0$. Let
$\mathcal{I}_{g}$ denote the set of all $g$-subhomogeneous index functions.

In order to study the optimality results of convergence rates for SAR (5), we
need to consider the following classes of decreasing functions:

$\mathcal{S}_{C_{\sigma}}:=\left\\{\varphi\in\mathcal{I}:\forall\lambda\in(0,\|A\|^{2}],~{}\varphi(\lambda)/\varphi(t^{-1})\leq
C_{\sigma}e^{\lambda t\sigma}\right\\},$ (10)

where positive numbers $C_{\sigma}$ and $\sigma$ are independent on $\lambda$
and $t$.

$\mathcal{S}^{g}_{\zeta}=\left\\{\begin{array}[]{ll}\varphi\in\mathcal{I}_{g}:~{}\text{There
exists a monotonically decreasing, integrable function~{}}\zeta\\\
\qquad\qquad\text{~{}such that~{}}e^{2\lambda
t}\zeta\left(\frac{\varphi(\lambda)}{\varphi(t^{-1})}\right)\geq
1\text{~{}for~{}}\lambda t>1.\end{array}\right\\}.$ (11)

In order to investigate the properties of sets $\mathcal{S}_{C_{\sigma}}$ and
$\mathcal{S}^{g}_{\zeta}$, we introduce the following concept, which is
generalized from [20, Definition 2]:

######  Definition 3.

A function $\varphi$ is said to be _covered_ by another function $\psi$ if
there are $\underline{c},\bar{T}>0$ such that, for all $t\geq\bar{T}$,

$\underline{c}\,\frac{\psi(t^{-1})}{\varphi(t^{-1})}\leq\min\limits_{t^{-1}\leq\lambda\leq\|A\|^{2}}\frac{\psi(\lambda)}{\varphi(\lambda)}.$
(12)

The following assertions are straightforward:

######  Proposition 2.

(i) $\varphi\in\mathcal{S}_{C_{\sigma}}$ if $\psi\in\mathcal{S}_{C_{\sigma}}$
and $\varphi$ is covered by $\psi$.

(ii) $\varphi\in\mathcal{S}^{g}_{\zeta}$ if $\psi\in\mathcal{S}^{g}_{\zeta}$
and $\varphi$ is covered by $\psi$.

(iii) If the quotient function
$\lambda\mapsto\frac{\psi(\lambda)}{\varphi(\lambda)}$ is increasing for
$0<\lambda\leq\bar{\lambda}$ and some $\bar{\lambda}>0$, $\varphi$ is covered
by $\psi$. If, in particular, $\varphi(\lambda)\in\mathcal{S}_{C_{\sigma}}$ is
concave for $0<\lambda\leq\bar{\lambda}$, $\varphi$ is covered by
$\psi(\lambda)=\lambda$.

Now, we study the elements in the sets $\mathcal{S}_{C_{\sigma}}$ and
$\mathcal{S}^{g}_{\zeta}$, which are related to the range-type source
conditions in the regularization theory of inverse problems:

$x_{0}-x^{\dagger}=\varphi(A^{*}A)v,\quad\|v\|\leq\rho,\quad\varphi\in\mathcal{S}_{C_{\sigma}}.$
(13)

In this paper, we are mainly interested in the following two examples of
source conditions:

######  Example 1.

(i) Hölder-type source conditions $\varphi_{p}$ ($p>0$):
$\varphi_{p}(\lambda)=\lambda^{p}$.

(ii) Logarithmic source conditions $\varphi_{\mu}$ ($\mu>0$):

$\varphi_{\mu}(\lambda)=\left\\{\begin{array}[]{ll}\log^{-\mu}(1/\lambda)\qquad\textrm{~{}for}\quad
0<\lambda\leq e^{-\mu-1},\\\ \mbox{arbitrarily extended as index function
for}\;\;\lambda>e^{-\mu-1}.\end{array}\right.$ (14)

The relationship between source condition functions $\varphi_{p}$,
$\varphi_{\mu}$ and the classes $\mathcal{S}_{C_{\sigma}}$ and
$\mathcal{S}^{g}_{\zeta}$ is investigated in the following lemma.

######  Lemma 1.

(i) $\varphi_{p}\in\mathcal{S}_{C_{\sigma}}\cap\mathcal{S}^{g}_{\zeta}$ with
$C_{\sigma}=p^{p}/(e^{p}\sigma^{p})$, arbitrary positive number $\sigma$,
$g(\lambda)=\lambda$ and $\zeta(\lambda)=(\ln\lambda/2)^{-p}$.

(ii) $\varphi_{\mu}\in\mathcal{S}_{C_{\sigma}}\cap\mathcal{S}^{g}_{\zeta}$
with arbitrary positive number $\sigma$.

###### Proof.

(i) From

$\sup_{t>0}e^{-\lambda\sigma
t}t^{p}=_{t=p/(\lambda\sigma)}p^{p}/(e^{p}\sigma^{p}\lambda^{p})$

we conclude that, for any $\sigma>0$, $\varphi_{p}\in\mathcal{S}_{C_{\sigma}}$
with $C_{\sigma}=p^{p}/(e^{p}\sigma^{p})$. The relation
$\varphi_{p}\in\mathcal{S}^{g}_{\zeta}$ holds, with $g(\lambda)=\lambda$ and
$\zeta(\lambda)=e(-2\lambda^{1/p})$.

(ii) Note that for arbitrary $\mu>0$ the index function
$\varphi_{\mu}(\lambda)$ is concave for all $0<\lambda\leq e^{-\mu-1}$, and
hence, because of (iii) of Proposition 2, is covered by
$\varphi^{1}(\lambda)=\lambda$. Consequently,
$\varphi_{\mu}\in\mathcal{S}_{C_{\sigma}}\cap\mathcal{S}^{g}_{\zeta}$. ∎

###  2.2 Regularization property of SAR

In this subsection, we perform the convergence analysis of the method (5) with
regard to regularization theory. Throughout this paper, we assume $x_{0}\in
N(A)^{\bot}$. Let $\\{\lambda_{j};u_{j},v_{j}\\}_{j=1}^{\infty}$ be the well-
defined singular system for the compact linear operator $A$, i.e. we have
$Au_{j}=\lambda_{j}v_{j}$ and $A^{*}v_{j}=\lambda_{j}u_{j}$ with ordered
singular values
$\|A\|=\lambda_{1}\geq\lambda_{2}\geq\cdot\cdot\cdot\geq\lambda_{j}\geq\lambda_{j+1}\geq\cdot\cdot\cdot\to
0$ as $j\to\infty$. Since the eigenelement $\\{u_{j}\\}_{j=1}^{\infty}$ forms
an orthogonal basis in $N(A)^{\bot}\subset\mathcal{X}$, according to the
construction of the approximate solution $x^{\delta}(t)$ by SAR in (5), it has
the decomposition $x^{\delta}(t)=\sum_{j}\xi_{j}(t)u_{j}$. Hence, by using

$\quad\langle dx^{\delta},u_{j}\rangle=\langle
A^{*}(y^{\delta}-Ax^{\delta})dt,u_{j}\rangle+\langle
f(t)dB_{t},u_{j}\rangle,~{}j=1,2,\cdots,$ (15)

we obtain the stochastic differential equations for coefficients
$\\{\xi_{j}(t)\\}^{\infty}_{j=1}$:

$d\xi_{j}(t)=\left(\lambda_{j}\langle
y^{\delta},v_{j}\rangle-\lambda^{2}_{j}\xi_{j}(t)\right)dt+f(t)\sqrt{q_{j}}d\beta_{j}(t),~{}\xi_{j}(0)=\langle
x_{0},u_{j}\rangle,~{}j=1,2,\cdots.$ (16)

######  Proposition 3.

The stochastic differential equation (16) has a unique solution,

$\displaystyle\xi_{j}(t)=e^{-\lambda^{2}_{j}t}\langle
x_{0},u_{j}\rangle+\frac{1-e^{-\lambda^{2}_{j}t}}{\lambda_{j}}\langle
y^{\delta},v_{j}\rangle+\int^{t}_{0}\sqrt{q_{j}}e^{-\lambda^{2}_{j}(t-s)}f(s)d\beta_{j}(s),$

where the stochastic integral
$\int^{t}_{0}\sqrt{q_{j}}e^{-\lambda^{2}_{j}(t-s)}f(s)d\beta_{j}(s)$ is
Gaussian with distribution

$\mathcal{N}\left(0,\int^{t}_{0}q_{j}e^{-2\lambda^{2}_{j}(t-s)}[f(s)]^{2}ds\right).$

Moreover, $\xi_{j}(t)$ is also Gaussian, with mean

$\mathbb{E}\xi_{j}(t)=e^{-\lambda^{2}_{j}t}\langle
x_{0},u_{j}\rangle+\frac{1-e^{-\lambda^{2}_{j}t}}{\lambda_{j}}\langle
y^{\delta},v_{j}\rangle$ (17)

and variance

$\textup{Var}\left(\xi_{j}(t)\right)=\int^{t}_{0}q_{j}e^{-2\lambda^{2}_{j}(t-s)}[f(s)]^{2}ds.$
(18)

Consequently, if $f\in\mathcal{I}$,
$\xi_{j}(t)\sim\mathcal{N}\left(\frac{\langle
y^{\delta},v_{j}\rangle}{\lambda_{j}},0\right)$ as $t\to\infty$.

###### Proof.

The existence and uniqueness of (16) is a standard result for linear
stochastic differential equations; see, e.g., [4]. (17) is from the mean-zero
property of stochastic integrals, and (18) is from the Itô isometry. ∎

From Proposition 3 and the decomposition
$x^{\delta}(t)=\sum_{j}\xi_{j}(t)u_{j}$ we obtain the explicit formula for the
mild solution of (5):

$\begin{array}[]{ll}x^{\delta}(t)=\sum\limits_{j}e^{-\lambda^{2}_{j}t}\langle
x_{0},u_{j}\rangle
u_{j}+\sum\limits_{j}\frac{1-e^{-\lambda^{2}_{j}t}}{\lambda_{j}}\langle
y^{\delta},v_{j}\rangle
u_{j}+\sum\limits_{j}\int^{t}_{0}e^{-\lambda^{2}_{j}(t-s)}f(s)\langle
dB_{s},u_{j}\rangle u_{j},\\\
=:(1-A^{*}Ag(t,A^{*}A))x_{0}+g(t,A^{*}A)A^{*}y^{\delta}+\int^{t}_{0}e^{-A^{*}A(t-s)}f(s)dB_{s},\end{array}$
(19)

where (we identify $\lambda^{2}_{j}$ as $\lambda$)

$\displaystyle g(t,\lambda)=\frac{1-e^{-\lambda t}}{\lambda}.$ (20)

######  Theorem 1.

Let $x^{\delta}(t)$ be the stochastic dynamic solution of (5) with
$f(t)\in\mathcal{I}$. Then, if the terminating time
$t^{*}=t^{*}(\delta,y^{\delta})$ is chosen so that

$\lim_{\delta\to 0}t^{*}=\infty\textrm{~{}and~{}}\lim_{\delta\to 0}\delta\cdot
t^{*}=0,$ (21)

the $x^{\delta}(t^{*})$ converges to $x^{\dagger}$ in the sense of mean
square, i.e. $\mathbb{E}\|x^{\delta}(t)-x^{\dagger}\|^{2}\to 0$, as $\delta\to
0$.

###### Proof.

According to the bias-variance decomposition for the mean-squared error,
namely,

$\mathbb{E}\|x^{\delta}(t)-x^{\dagger}\|^{2}=\|\mathbb{E}x^{\delta}(t)-x^{\dagger}\|^{2}+\mathbb{E}\|x^{\delta}(t)-\mathbb{E}x^{\delta}(t)\|^{2},$
(22)

it is sufficient to find estimates for both bias error
$\|\mathbb{E}x^{\delta}(t)-x^{\dagger}\|^{2}$ and variance error
$\mathbb{E}\|x^{\delta}(t)-\mathbb{E}x^{\delta}(t)\|^{2}$, respectively.

The bias error is deterministic, and thus can be estimated by standard
techniques in regularization theory, i.e. we have (let $r(t,\lambda)=1-\lambda
g(t,\lambda)=e^{-\lambda t}$)

$\begin{array}[]{ll}\|\mathbb{E}x^{\delta}(t)-x^{\dagger}\|=\|r(t,A^{*}A)(x_{0}-x^{\dagger})+g(t,A^{*}A)A^{*}(y^{\delta}-y)\|\\\
\qquad\leq\|e^{-tA^{*}A}(x_{0}-x^{\dagger})\|+\vartheta
t^{1/2}\delta,\end{array}$ (23)

where we have used the inequalities

$\|g(t,A^{*}A)A^{*}(y^{\delta}-y)\|\leq\sup_{\lambda\in(0,\|A\|^{2})}\sqrt{\lambda}g(t,\lambda)\|y^{\delta}-y\|\leq\delta\sup_{\lambda\in(0,\|A\|^{2})}\frac{1-e^{-\lambda
t}}{\sqrt{\lambda}}\leq\vartheta t^{1/2}\delta,$ (24)

where
$\vartheta=\sup_{\lambda\in\mathbf{R}_{+}}\sqrt{\lambda}(\lambda-e^{-\lambda})\approx
0.6382$ [26]. Hence, with the choice of terminating time in (21), we conclude
that

$\|\mathbb{E}x^{\delta}(t)-x^{\dagger}\|\to 0\textrm{~{}as~{}}\delta\to 0.$
(25)

On the other hand, according to L’Hospital’s rule, we get, for any
$\lambda>0$,

$\lim_{t\to\infty}\frac{\int_{0}^{t}e^{\lambda
s}[f(s)]^{2}ds}{\int^{t}_{0}e^{\lambda s}ds}=\lim_{t\to\infty}\frac{e^{\lambda
t}[f(t)]^{2}}{e^{\lambda t}}=\lim_{t\to\infty}[f(t)]^{2}.$ (26)

Hence, $T$ exists such that, for all $t\geq T$,

$\frac{1}{2}[f(t)]^{2}\int^{t}_{0}e^{\lambda s}ds\leq\int_{0}^{t}e^{\lambda
s}[f(s)]^{2}ds\leq 2[f(t)]^{2}\int^{t}_{0}e^{\lambda s}ds.$ (27)

By using the Itô isometry, we derive, together with the assumption
$f(t)\in\mathcal{I}$, the assumption $tr(Q(A^{*}A)^{-1})<\infty$ and the
inequality (27) such that, for $t\geq T$,

$\begin{array}[]{ll}\mathbb{E}\|x^{\delta}(t)-\mathbb{E}x^{\delta}(t)\|^{2}=\mathbb{E}\left\|\int^{t}_{0}e^{-A^{*}A(t-s)}f(s)dB_{s}\right\|^{2}=\sum\limits^{\infty}\limits_{j=1}\int^{t}_{0}q_{j}e^{-2\lambda^{2}_{j}(t-s)}[f(s)]^{2}ds\\\
\qquad=\sum\limits^{\infty}\limits_{j=1}q_{j}e^{-2\lambda^{2}_{j}t}\int^{t}_{0}e^{2\lambda^{2}_{j}s}[f(s)]^{2}ds\leq
2\sum\limits^{\infty}\limits_{j=1}q_{j}e^{-2\lambda^{2}_{j}t}[f(t)]^{2}\int^{t}_{0}e^{2\lambda^{2}_{j}s}ds\\\
\qquad=[f(t)]^{2}\sum\limits^{\infty}\limits_{j=1}\frac{q_{j}}{\lambda^{2}_{j}}\left(1-e^{-2\lambda^{2}_{j}t}\right)\leq
tr(Q(A^{*}A)^{-1})[f(t)]^{2}\to 0\end{array}$ (28)

as $t\to\infty$. We combine (25) and (28) to obtain the convergence of the
full error of our regularization method (5). ∎

###  2.3 Convergence rates with noisy data under _a priori_ and _a posteriori_
stopping rules

According to the standard argument in regularization theory [24], under the
general assumptions of the previous sections, the rate of convergence of mean-
square error $\mathbb{E}\|x^{\delta}(t(\delta))-x^{\dagger}\|^{2}\to 0$ as
$\delta\to 0$ can be arbitrarily slow for solutions $x^{\dagger}$ which are
not smooth enough. In order to prove convergence rates, some smoothness
assumptions imposed on the exact solution must be employed. This subsection
focuses on the range-type source conditions (13), while in the next subsection
we adopt the variational inequalities, which are more natural in the study of
converse results. Therefore, we first consider the convergence-rate results
under an _a priori_ stopping rule.

######  Theorem 2.

Let $x^{\delta}(t)$ be the stochastic dynamic solution of (5) with
$f(t)=\mathcal{O}(\varphi(1/t))$. Then, under the source condition (13), if
the terminating time is chosen as $t^{*}=\Theta^{-1}(\delta)$ with
$\Theta(t)=t^{-1/2}\varphi(t^{-1})$, we have the convergence rate

$\mathbb{E}\|x^{\delta}(t^{*})-x^{\dagger}\|^{2}=\mathcal{O}([\varphi([\Theta^{-1}(\delta)]^{-1})]^{2})\text{~{}as~{}}\delta\to
0.$

Consequently, if $\varphi=\varphi_{p}$, we have
$\mathbb{E}\|x^{\delta}(t^{*})-x^{\dagger}\|^{2}=\mathcal{O}(\delta^{\frac{4p}{2p+1}})$;
if $\varphi=\varphi_{\mu}$, we have
$\mathbb{E}\|x^{\delta}(t^{*})-x^{\dagger}\|^{2}=\mathcal{O}(\log^{-2\mu}(\delta^{-1}))$.

###### Proof.

According to the proof of Theorem 1 (see (23) and (28)), we find, together
with the assumption $f(t)=\mathcal{O}(\varphi(1/t))$ and the relation
$t^{*}(\delta)\to\infty$ as $\delta\to 0$, that for a small-enough $\delta$
there exists a constant $C>0$ such that

$\begin{array}[]{ll}\mathbb{E}\|x^{\delta}(t)-x^{\dagger}\|^{2}=\|\mathbb{E}x^{\delta}(t)-x^{\dagger}\|^{2}+\mathbb{E}\|x^{\delta}(t)-\mathbb{E}x^{\delta}(t)\|^{2}\\\
\qquad\leq\left(\|e^{-tA^{*}A}(x_{0}-x^{\dagger})\|+\vartheta
t^{1/2}\delta\right)^{2}+tr(Q(A^{*}A)^{-1})[f(t)]^{2}\\\ \qquad\leq
2C^{2}_{1}\rho^{2}[\varphi(1/t)]^{2}+2\vartheta^{2}t\delta^{2}+C^{2}[\varphi(1/t)]^{2},\end{array}$
(29)

which yields the required result through the selection method of $t^{*}$. ∎

For the _a posteriori_ stopping rule, we consider the following stochastic
version of Morozov’s discrepancy principle, i.e. $t^{*}$ is chosen such that

$t^{*}_{i}:=\inf\left\\{t>0:\chi_{i}(t)<0\right\\},\quad i=1,2,$ (30)

where

$\chi_{1}(t):=\|A\mathbb{E}x^{\delta}(t)-y^{\delta}\|-\tau\delta,$ (31)

and

$\chi_{2}(t):=\mathbb{E}\|Ax^{\delta}(t)-y^{\delta}\|^{2}-\tau\delta^{2}\quad\text{with~{}}f(t)\in\mathcal{I},$
(32)

where we assume $\tau>1$ for the occurring factor of the noise level $\delta$.

######  Remark 1.

Clearly, $t^{*}_{i}$ is the root of $\chi_{i}(t)$, i.e.
$\chi_{i}(t^{*}_{i})=0$. The formulation (30) is more appropriate for a
practical realization of SAR when $t^{*}$ is the first time point for which
the size of the sample residual $\|Ax^{\delta}(t)-y^{\delta}\|$ is
approximately of the order of the data error.

######  Proposition 4.

If $\|Ax_{0}-y^{\delta}\|>\tau\delta$, there always exists a unique
$t^{*}_{i}$ ($i=1,2$), defined by (30).

###### Proof.

From the explicit formula of $x^{\delta}(t)$ (see (19)), we derive

$\begin{array}[]{ll}\chi_{1}(t)&=\|A\mathbb{E}x^{\delta}(t)-y^{\delta}\|-\tau\delta=\left\|Ae^{-tA^{*}A}(x_{0}-x^{\dagger})+e^{-tAA^{*}}(y-y^{\delta})\right\|-\tau\delta\\\
&\leq\left\|e^{-tAA^{*}}A(x_{0}-x^{\dagger})\right\|-(\tau-1)\delta\to-(\tau-1)\delta<0\end{array}$

as $t\to\infty$. The continuity of $\chi_{1}(t)$ is clear as we are dealing
with the linear problem. Since $\chi(0)=\|Ax_{0}-y^{\delta}\|-\tau\delta>0$,
the existence of the root of $\chi(t)$ follows from Bolzano’s theorem. Using
similar reasoning, we can obtain the existence and uniqueness of $t^{*}_{2}$.
∎

We are now able to provide the convergence-rate results of SAR under the _a
posteriori_ stopping rule (30).

######  Theorem 3.

Let the terminating time $t^{*}=t^{*}(\delta,y^{\delta})$ of SAR (5) with
$f(t)=\mathcal{O}(\varphi(1/t))$ be chosen according to the stopping rule
(30). Then, we have the following convergence-rate results:

(i) Under the Hölder-type source conditions $\varphi_{p}$:

$t^{*}=\mathcal{O}\left(\delta^{-\frac{2}{2p+1}}\right)\qquad\mbox{and}\qquad\mathbb{E}\|x^{\delta}(t^{*})-x^{\dagger}\|^{2}=\mathcal{O}\left(\delta^{\frac{4p}{2p+1}}\right).$
(33)

(ii) Under the logarithmic source conditions $\varphi_{\mu}$:

$t^{*}=o\left(\delta^{-\frac{2}{2p+1}}\log^{-\frac{2}{2p+1}}(\delta^{-1})\right)\quad\mbox{and}\quad\mathbb{E}\|x^{\delta}(t^{*})-x^{\dagger}\|^{2}=\mathcal{O}\left(\log^{-2\mu}(\delta^{-1})\right).$
(34)

###### Proof.

We show only the convergence-rate results under the Hölder-type source
conditions and the stopping rule (30) with discrepancy function $\chi_{1}(t)$.
Other cases (i.e. the case of the stopping rule (30) with $\chi_{2}(t)$ and
the case with the logarithmic source condition) can be proved in a similar
way. To that end, using the interpolation inequality
$\|B^{p}u\|\leq\|B^{q}u\|^{p/q}\|u\|^{1-p/q}$ and the source conditions
$x_{0}-x^{\dagger}=(A^{*}A)^{p}v$, we deduce that

$\displaystyle\begin{array}[]{ll}\|e^{-tA^{*}A}(x_{0}-x^{\dagger})\|=\|(A^{*}A)^{p}e^{-tA^{*}A}v\|\\\
\quad\leq\|(A^{*}A)^{(p+1/2)}e^{-tA^{*}A}v\|^{2p/(2p+1)}\cdot\|e^{-tA^{*}A}v\|^{1/(2p+1)}\\\
\quad=\|Ae^{-tA^{*}A}(x_{0}-x^{\dagger})\|^{2p/(2p+1)}\cdot\|e^{-tA^{*}A}v\|^{1/(2p+1)}.\end{array}$
(38)

Since $t^{*}_{1}$ is chosen according to the equation $\chi_{1}(t)=0$, we
derive

$\displaystyle\begin{array}[]{ll}\tau\delta=\|A\mathbb{E}x^{\delta}(t^{*}_{1})-y^{\delta}\|\\\
\qquad=\left\|Ae^{-t^{*}_{1}A^{*}A}(x_{0}-x^{\dagger})+e^{-t^{*}_{1}A^{*}A}(y-y^{\delta})\right\|\\\
\qquad\geq\|Ae^{-t^{*}_{1}A^{*}A}(x_{0}-x^{\dagger})\|-\|e^{-t^{*}_{1}A^{*}A}(y-y^{\delta})\|.\end{array}$
(42)

Now we combine the estimates (38) and (42) to obtain, with the source
conditions,

$\displaystyle\begin{array}[]{ll}\|e^{-t^{*}_{1}A^{*}A}(x_{0}-x^{\dagger})\|\leq\|Ae^{-t^{*}_{1}A^{*}A}(x_{0}-x^{\dagger})\|^{2p/(2p+1)}\cdot\|e^{-t^{*}_{1}A^{*}A}v\|^{1/(2p+1)}\\\
\quad\leq\left(\tau\delta+\|e^{-t^{*}_{1}A^{*}A}(y^{\delta}-y)\|\right)^{2p/(2p+1)}\rho^{1/(2p+1)}\\\
\quad\leq c_{1}\rho^{1/(2p+1)}\delta^{2p/(2p+1)}\end{array}$ (46)

where $c_{1}:=\left(\tau+1\right)^{2p/(2p+1)}$.

On the other hand, in a similar fashion to (42), it is easy to show that

$\begin{array}[]{ll}\tau\delta&=\|A\mathbb{E}x^{\delta}(t^{*}_{1})-y^{\delta}\|=\|e^{-t^{*}_{1}AA^{*}}A(x_{0}-x^{\dagger})+e^{-t^{*}_{1}A^{*}A}(y-y^{\delta})\|\\\
&\leq\|Ae^{-t^{*}_{1}A^{*}A}(x_{0}-x^{\dagger})\|+\|e^{-t^{*}_{1}A^{*}A}(y-y^{\delta})\|\\\
&\leq\|Ae^{-t^{*}_{1}A^{*}A}(x_{0}-x^{\dagger})\|+\delta.\end{array}$ (47)

If we combine the above inequality with the source conditions (13) with
$\varphi=\varphi_{p}$, we obtain

$\displaystyle(\tau-1)\delta\leq\|Ae^{-t^{*}_{1}A^{*}A}(y-y^{\delta})\|\leq\|(A^{*}A)^{p+1/2}e^{-t^{*}_{1}A^{*}A}v\|\leq\rho
C_{1}(t^{*}_{1})^{-(2p+1)/2},$ (48)

which yields the estimate for $t^{*}_{1}$ in (33). Finally, according to the
proof of Theorem 1 (cf. (23)), the estimate for $t^{*}$, and (46), we conclude
that a constant $C>0$ exists such that, for a small-enough $\delta$,

$\displaystyle\begin{array}[]{ll}\mathbb{E}\|x^{\delta}(t^{*})-x^{\dagger}\|^{2}\leq\|e^{-t^{*}_{1}A^{*}A}(x_{0}-x^{\dagger})\|^{2}+\vartheta^{2}t^{*}_{1}\delta^{2}+tr(Q(A^{*}A)^{-1})[f(t^{*}_{1})]^{2}\\\
\qquad\leq
c^{2}_{1}\rho^{2/(2p+1)}\delta^{4p/(2p+1)}+\vartheta^{2}\left(\frac{\rho
C_{1}}{\tau-1}\right)^{2/(2p+1)}\delta^{4p/(2p+1)}+C\left[\varphi_{p}\left(\delta^{\frac{2}{2p+1}}\right)\right]^{2},\end{array}$

which completes the proof. ∎

###  2.4 Convergence rates with exact data and converse results

LetDefine the _spectral tail_ be defined by

$\omega(\lambda)=\|E_{[0,\lambda]}(x_{0}-x^{\dagger})\|^{2},$ (50)

where the mapping $\Theta\mapsto E_{\Theta}$ means the spectral measure of
$A^{*}A$ on Borel set $\Theta\subset[0,\infty)$. Then, we have the
representation

$\|\mathbb{E}x(t)-x^{\dagger}\|^{2}=\int^{\|A\|^{2}}_{0}e^{-\lambda
t}d\omega(\lambda),$ (51)

where $x(t)$ represents the regularized stochastic solution $x(t)$ of (5) with
exact data $y$. In this subsection, we establish an equivalent relation
between the spectral tail $\omega(\cdot)$ and the convergence rate of $x(t)$
to the $x_{0}$-minimum norm solution $x^{\dagger}$.

######  Theorem 4.

Let $\varphi\in\mathcal{S}_{C_{\sigma}}$ with $\sigma\in(0,1)$. Then, the
following two statements are equivalent:

  * (i)

There exists a constant $C_{3}>0$ with

$\left\|\mathbb{E}x(t)-x^{\dagger}\right\|^{2}\leq
C_{3}\varphi(1/t)\quad\text{for all }t>0.$ (52)

  * (ii)

There exists a constant $C_{4}>0$ with

$\omega(\lambda)\leq C_{4}\varphi(\lambda)\quad\text{for all }\lambda>0.$ (53)

Moreover, if $f(t)=C_{\varphi}\sqrt{\varphi(1/t)}$ ($C_{\varphi}$ is a
positive number independent of $t$), every one of the above statements is also
equivalent to the following one:

  * (iii)

There exists a constant $C_{5}>0$ with

$\mathbb{E}\left\|x(t)-x^{\dagger}\right\|^{2}\leq
C_{5}\varphi(1/t)\quad\text{for all }t>0.$ (54)

###### Proof.

According to the definitions of $x(t)$ and $\omega(\lambda)$ and the relation
(51), for all $t>0$,

$\omega(1/t)=\int_{0}^{1/t}\mathrm{d}\omega(\lambda)\leq
e\int_{0}^{1/t}e^{-\lambda t}d\omega(\lambda)\leq
e\left\|\mathbb{E}x(t)-x^{\dagger}\right\|^{2}.$ (55)

Let fFirst, let (52) hold. Then, (53) holds with $C_{4}=eC_{3}$, through theby
combining of (52) and (55).

Conversely, let (53) hold. Since
$\left\|\mathbb{E}x(t)-x^{\dagger}\right\|^{2}\leq\left\|x_{0}-x^{\dagger}\right\|^{2}$
(which follows from the representation (51)), it is sufficientenough to check
the condition (52) for all $t>1/\|A\|^{2}$.

Integrating the right- hand side by parts in (51), we obtain that

$\left\|\mathbb{E}x(t)-x^{\dagger}\right\|^{2}=e^{-\|A\|^{2}t}\omega\left(\|A\|^{2}\right)+t\int_{0}^{\|A\|^{2}}\omega(\lambda)e^{-\lambda
t}d\lambda$ (56)

Note that

$t\int_{0}^{1/t}\omega(\lambda)e^{-\lambda t}d\lambda\leq
t\omega(1/t)\int_{0}^{1/t}e^{-\lambda t}d\lambda=(1-e^{-1})\omega(1/t)\leq
C_{4}(1-e^{-1})\varphi(1/t),$ (57)

where we used the assumption (53) and the monotonic increasing of $\omega$. On
the other hand, the assumptions $\varphi\in\mathcal{S}_{C_{\sigma}}$ and (53)
imply that

$\displaystyle t$
$\displaystyle\int_{1/t}^{\|A\|^{2}}\omega(\lambda)e^{-\lambda t}d\lambda\leq
C_{4}t\int_{1/t}^{\|A\|^{2}}\varphi(\lambda)e^{-\lambda t}d\lambda$ (58)
$\displaystyle=C_{4}t\int_{1/t}^{\|A\|^{2}}\varphi(\lambda)e^{-\lambda
t\sigma}e^{-\lambda t(1-\sigma)}d\lambda\leq
C_{\sigma}C_{4}\varphi(1/t)t\int_{1/t}^{\|A\|^{2}}e^{-\lambda
t(1-\sigma)}d\lambda$
$\displaystyle\leq\frac{C_{\sigma}C_{4}}{(1-\sigma)e^{1-\sigma}}\varphi(1/t).$

Inserting (57) and (58) into (56), we find that, with
$\omega\left(\|A\|^{2}\right)=\left\|x_{0}-x^{\dagger}\right\|^{2}$, that

$\left\|\mathbb{E}x(t)-x^{\dagger}\right\|^{2}\leq
e^{-\|A\|^{2}t}\left\|x_{0}-x^{\dagger}\right\|^{2}+C_{4}(1-e^{-1})\varphi(1/t)+\frac{C_{\sigma}C_{4}}{(1-\sigma)e^{1-\sigma}}\varphi(1/t)$
(59)

From the definition of $\mathcal{S}_{C_{\sigma}}$, we deduce further that

$e^{-\|A\|^{2}t}\leq\frac{C_{\sigma}^{1/\sigma}}{\varphi^{1/\sigma}\left(\|A\|^{2}\right)}\varphi^{1/\sigma}(1/t)=\frac{C_{\sigma}^{1/\sigma}\varphi^{1/\sigma-1}(1/t)}{\varphi^{1/\sigma}\left(\|A\|^{2}\right)}\varphi(1/t)\leq\frac{C_{\sigma}^{1/\sigma}}{\varphi\left(\|A\|^{2}\right)}\varphi(1/t),$

since $\varphi$ is increasing and $\sigma<1$.

Thus, we obtain,get from (59), that

$\left\|\mathbb{E}x(t)-x^{\dagger}\right\|^{2}\leq C_{3}\varphi(1/t)$

with

$C_{3}=\frac{C_{\sigma}^{1/\sigma}}{\varphi\left(\|A\|^{2}\right)}\left\|x_{0}-x^{\dagger}\right\|^{2}+C_{4}(1-e^{-1})+\frac{C_{\sigma}C_{4}}{(1-\sigma)e^{1-\sigma}}.$

Finally, the equivalence of inequalities (52) and (54) follows from the
following two inequalities for a large- enough $t$:

$\displaystyle\mathbb{E}\|x^{\delta}(t)-x^{\dagger}\|^{2}=\|\mathbb{E}x^{\delta}(t)-x^{\dagger}\|^{2}+\mathbb{E}\|x^{\delta}(t)-\mathbb{E}x^{\delta}(t)\|^{2}$
$\displaystyle\leq\|\mathbb{E}x^{\delta}(t)-x^{\dagger}\|^{2}+tr(Q(A^{*}A)^{-1})[f(t)]^{2}\leq\|\mathbb{E}x^{\delta}(t)-x^{\dagger}\|^{2}+C\varphi(1/t),$
$\displaystyle\mathbb{E}\|x^{\delta}(t)-x^{\dagger}\|^{2}\geq\|\mathbb{E}x^{\delta}(t)-x^{\dagger}\|^{2}+\frac{1}{2}[f(t)]^{2}\sum\limits^{\infty}\limits_{j=1}\frac{q_{j}}{\lambda^{2}_{j}}\left(1-e^{-2\lambda^{2}_{j}t}\right)$
$\displaystyle\geq\|\mathbb{E}x^{\delta}(t)-x^{\dagger}\|^{2}+\frac{1}{2}tr(Q(A^{*}A)^{-1})[f(t)]^{2}\geq\|\mathbb{E}x^{\delta}(t)-x^{\dagger}\|^{2}+C^{\prime}\varphi(1/t),$

which is obtained from the bias-variance decomposition (22) and the assumption
$f(t)=C_{\varphi}\sqrt{\varphi(1/t)}$. Here, $C$ and $C^{\prime}$ are two
fixed numbers independent of $t$. ∎

######  Remark 2.

The inequality (53) has close connections to variational inequalities and
range-type source conditions, which are frequently used in regularization
theory, (see, e.g., [7, 9, 10, 20]). To be more precise, let
$\varphi:[0,\infty)\to[0,\infty)$ be an increasing, continuous function and
$\nu\in(0,1)$. Then, according to [1], the following two statements are
equivalent:

(i) There exists a constant $C>0$ with

$\omega(\lambda)\leq C_{a}\varphi^{2\nu}(\lambda)\quad\text{for all}\
\lambda>0.$

(ii) There exists a constant $C_{b}>0$ such that

$\left|\langle x_{0}-x^{\dagger},x\rangle\right|\leq
C_{b}\|\varphi(L^{*}L)x\|^{\nu}\|x\|^{1-\nu}\quad\text{for all}\ x\in X.$ (60)

Moreover, the range-type source conditions, (i.e.
$x_{0}-x^{\dagger}\in\mathcal{R}(\psi^{\nu}(A^{*}A))$) implyies the
variational inequality (60). Conversely, the variational inequality (60)
implies that the relation
$x_{0}-x^{\dagger}\in\mathcal{R}(\psi^{\nu}(L^{*}L))$ holds for every
continuous function $\varphi:[0,\infty)\to[0,\infty)$ with $\psi\geq
c\varphi^{\mu}$ for some constant $c>0$ and some $\mu\in(0,\nu)$.

###  2.5 The best worst-case mean-square error

Let
$B_{\delta}(y):=\left\\{\tilde{y}\in\mathcal{Y}:~{}\|\tilde{y}-y\|\leq\delta\right\\}$
denote the ball containing all possible measurement data $\tilde{y}$ with
given accuracy level $\delta$, and let $x(t;\tilde{y})$ be the stochastic
dynamic solution of (5), with $y^{\delta}$ replaced with $\tilde{y}\in
B_{\delta}(y)$. In the case of $\tilde{y}=y$, we have $x(t)\equiv x(t;y)$. In
this subsection, we are interested in the convergence-rate results for the
_best worst-case mean-square error_ $\sup_{\tilde{y}\in
B_{\delta}(y)}\inf_{t>0}\mathbb{E}\|x(t,\tilde{y})-x^{\dagger}\|^{2}$, which
represents the distance between the $x_{0}$-minimum norm solution
$x^{\dagger}$ and the regularized stochastic solution $x(t,\tilde{y})$ for
some data $\tilde{y}$ belonging to the ball $B_{\delta}(y)$ under the optimal
choice of the regularization parameter $t$. The proof technique is similar to
that in [1], but some proof details were simplified.

We first investigate the convergence-rate results of best worst-case mean-
square error in some special situations.

######  Lemma 2.

Assume that there exists a constant $C_{6}>0$ such that

$\mathbb{E}\|x(t)-x^{\dagger}\|=0\emph{ for all }t\geq C_{6}.$

Then, we have

$\sup_{\tilde{y}\in
B_{\delta}(y)}\inf_{t>0}\mathbb{E}\|x(t,\tilde{y})-x^{\dagger}\|^{2}\leq
2C_{6}\vartheta^{2}\delta^{2}.$ (61)

###### Proof.

Let $\tilde{y}\in B_{\delta}(y)$ be fixed. With the use of the inequality
(24), it follows that, together with the identity $Ag(t,A^{*}A)=g(t,AA^{*})A$,

$\mathbb{E}\left\|x(t;\tilde{y})-x(t)\right\|^{2}=\left\langle\tilde{y}-y,g_{\alpha}^{2}(AA^{*})AA^{*}(\tilde{y}-y)\right\rangle\leq\delta^{2}\max_{\lambda>0}\lambda
g^{2}(t,\lambda)\leq\vartheta^{2}\delta^{2}t.$ (62)

Since the right-hand side is uniform for all $\tilde{{y}}\in B_{\delta}(y)$,
picking $t=C_{6}$, we get

$\sup_{\tilde{y}\in
B_{\delta}(y)}\inf_{t>0}\mathbb{E}\|x(t;\tilde{y})-x^{\dagger}\|^{2}\leq\inf_{t>0}\left\\{2\mathbb{E}\|x(t)-x^{\dagger}\|^{2}+2\vartheta^{2}\delta^{2}t\right\\}\leq
2C_{6}\vartheta^{2}\delta^{2},$

which is (61). ∎

######  Lemma 3.

Suppose that $\mathbb{E}\|x(t)-x^{\dagger}\|>0$ for all $t>0$. If we choose
for every $\delta>0$ the parameter $t_{\delta}>0$ such that

$\mathbb{E}\|x(t_{\delta})-x^{\dagger}\|^{2}=\delta^{2}t_{\delta},$ (63)

there exist two positive numbers $C_{7}$ and $C_{8}$ such that

$C_{8}\delta^{2}t_{\delta}\leq\sup_{\tilde{y}\in
B_{\delta}(y)}\inf_{t>0}\mathbb{E}\|x(t,\tilde{y})-x^{\dagger}\|^{2}\leq
C_{7}\delta^{2}t_{\delta}\,\,\text{~{}for all~{}}\,\delta>0.$ (64)

###### Proof.

By combining (19), (22), (28), and (51), we have

$\mathbb{E}\|x(t)-x^{\dagger}\|^{2}=\int^{\|A\|^{2}}_{0}e^{-\lambda
t}d\omega(\lambda)+\sum\limits^{\infty}\limits_{j=1}\int^{t}_{0}q_{j}e^{-2\lambda^{2}_{j}(t-s)}[f(s)]^{2}ds.$
(65)

Note that the function

$\xi(t)=t^{-1}\mathbb{E}\|x(t)-x^{\dagger}\|^{2}=\int^{\|A\|^{2}}_{0}t^{-1}e^{-\lambda
t}d\omega(\lambda)+t^{-1}\sum\limits^{\infty}\limits_{j=1}\int^{t}_{0}q_{j}e^{-2\lambda^{2}_{j}(t-s)}[f(s)]^{2}ds$

is, according to the assumption that $\mathbb{E}\|x(t)-x^{\dagger}\|>0$ for
all $t>0$, continuous and strictly decreasing, and satisfies
$\lim_{t\to\infty}\xi(t)=0$ and $\lim_{t\to 0}\xi(t)=\infty$. Therefore, we
find for every $\delta>0$ a unique value $t_{\delta}=\xi^{-1}(\delta^{2})$.

Let $\tilde{y}\in B_{\delta}(y)$. Then, as in the proof of Lemma 2 (see (62)),
we find that

$\mathbb{E}\|x(t;\tilde{y})-x(t)\|^{2}\leq\vartheta^{2}\delta^{2}t.$

From this estimate, we obtain, with the triangular inequality and the
definition (63) of $t_{\delta}$,

$\sup_{\tilde{y}\in
B_{\delta}(y)}\inf_{t>0}\mathbb{E}\|x(t;\tilde{y})-x^{\dagger}\|^{2}\leq\inf_{t>0}\left\\{2\mathbb{E}\|x(t)-x^{\dagger}\|^{2}+2\vartheta^{2}\delta^{2}t\right\\}^{2}\leq
2(1+\vartheta^{2})\delta^{2}t_{\delta},$

which is the upper bound in (64) with the constant $C_{1}=2(1+\vartheta^{2})$.

For the lower bound in (64), we write, similarly,

$\displaystyle\mathbb{E}\|x(t;\tilde{y})-x^{\dagger}\|^{2}$
$\displaystyle=\mathbb{E}\|x(t)-x^{\dagger}\|^{2}+\mathbb{E}\|x(t;\tilde{y})-x(t)\|^{2}+2\mathbb{E}\langle
x(t;\tilde{y})-x(t),x(t)-x^{\dagger}\rangle$ (66)
$\displaystyle=\mathbb{E}\|x(t)-x^{\dagger}\|^{2}+\mathbb{E}\langle\tilde{y}-y,\,g^{2}(t,AA^{*})AA^{*}(\tilde{y}-y)\rangle$
$\displaystyle\qquad+2\mathbb{E}\langle
g(t,AA^{*})(\tilde{y}-y),\,g(t,AA^{*})AA^{*}y-y\rangle.$

Note that for every $t>0$ there exists a large-enough number
$T_{\delta}>t_{\delta}$ such that $[T^{-1}_{\delta},t^{-1}_{\delta}]$ contains
at least one eigenvalue of $AA^{*}$ and the inequality
$1-e^{-t_{\delta}/T_{\delta}}\geq C_{A}$ holds with a fixed small positive
number $C_{A}$, dependent on only the distribution of eigenvalues
$\\{\sigma^{2}_{j}\\}^{\infty}_{j=1}$ of $AA^{*}$. For instance, if
$\max\limits_{j}\frac{\sigma^{2}_{j}}{\sigma^{2}_{j+2}}=:\kappa_{A}<\infty$,
we set $T_{\delta}=\kappa_{A}t_{\delta}$, and consequently can choose
$C_{A}=1-e^{-1/\kappa_{A}}$. In view of this, we let

$z_{\delta}=F_{[T^{-1}_{\delta},t^{-1}_{\delta}]}(g(t_{\delta},AA^{*})AA^{*}y-y),$
(67)

where $F$ represents the spectral measure of the operator $AA^{*}$. Through
the choice of $T_{\delta}$, $F_{[T^{-1}_{\delta},t^{-1}_{\delta}]}\neq 0$.

We first consider the case when $z_{\delta}\neq 0$. Then, if we set
$\tilde{y}=y+\delta\frac{z_{\delta}}{\|z_{\delta}\|}$ , the equation (66)
becomes

$\mathbb{E}\|x(t,\tilde{y})-x^{\dagger}\|^{2}=\mathbb{E}\|x(t)-x^{\dagger}\|^{2}+\frac{\delta^{2}}{\|z_{\delta}\|^{2}}\mathbb{E}\langle
z_{\delta},g^{2}(t,AA^{*})AA^{*}z_{\delta}\rangle+\frac{2\delta}{\|z_{\delta}\|}\mathbb{E}\langle
g(t,AA^{*})z_{\delta},z_{\delta}\rangle.$

We may drop the last term as it is non-negative, which gives us the lower
bound

$\sup_{\tilde{y}\in
B_{\delta}(y)}\inf_{t>0}\mathbb{E}\|x(t,\tilde{y})-x^{\dagger}\|^{2}\geq\inf_{t>0}\left(\mathbb{E}\|x(t)-x^{\dagger}\|^{2}+\delta^{2}\min_{\lambda\in[T^{-1}_{\delta},t^{-1}_{\delta}]}\lambda
g^{2}(t,\lambda)\right).$

Now, from

$\lambda g^{2}(t,\lambda)=\frac{(1-e^{-\lambda
t})^{2}}{\lambda}\geq\frac{(1-e^{-T^{-1}_{\delta}t})^{2}}{t^{-1}_{\delta}}\,\,\text{for
all}\,\lambda\in[T^{-1}_{\delta},t^{-1}_{\delta}],$

we can estimate further:

$\sup_{\tilde{y}\in
B_{\delta}(y)}\inf_{t>0}\mathbb{E}\|x(t,\tilde{y})-x^{\dagger}\|^{2}\geq\inf_{t>0}\left(\mathbb{E}\|x(t)-x^{\dagger}\|^{2}+\delta^{2}t_{\delta}\left(1-e^{-T^{-1}_{\delta}t}\right)^{2}\right).$

Now, since the first term is decreasing in $t$ (see (6)) and the second term
is increasing in $t$, we can estimate the expression for $t>t_{\delta}$ from
below using the second term at $t=t_{\delta}$, and for $t\leq t_{\delta}$
using the first term at $t=t_{\delta}$:

$\sup_{\tilde{y}\in
B_{\delta}(y)}\inf_{t>0}\mathbb{E}\|x(t,\tilde{y})-x^{\dagger}\|^{2}\geq\min\left\\{\mathbb{E}\|x(t)-x^{\dagger}\|^{2},\delta^{2}t_{\delta}(1-e^{-T^{-1}_{\delta}t_{\delta}})^{2}\right\\}\geq(1-e^{-1})^{2}\delta^{2}t_{\delta},$

which is the lower bound in (64) with $C_{8}=(1-e^{-1})^{2}$.

If $z_{\delta}$, as defined by (67), happens to vanish, the same argument
works with an arbitrary non-zero element
$z_{\delta}\in\mathcal{R}(F_{[T^{-1}_{\delta},t^{-1}_{\delta}]})$ since the
last term in (66) is zero for
$\tilde{y}=y+\delta\frac{z_{\delta}}{\|z_{\delta}\|}$. ∎

From Lemma 2 and Lemma 3, we now get an equivalence relation between the noisy
and noise-free convergence rates.

######  Theorem 5.

Let $\phi(1/\cdot)\equiv\varphi(\cdot)\in\mathcal{S}^{g}_{\zeta}$ and, denote
by

$\tilde{\phi}(t)=\sqrt{t^{-1}\phi(t)}\ \text{and }\
\psi(\delta)=\delta^{2}\tilde{\phi}^{-1}(\delta).$ (68)

Then, the following two statements are equivalent:

  * (a)

There exists a constant $c>0$ such that

$\sup_{\tilde{y}\in\tilde{B}_{\delta}(y)}\inf_{t>0}\mathbb{E}\|x(t,\tilde{y})-x^{\dagger}\|^{2}\leq
c\psi(\delta)\quad\text{for all}\ \delta>0.$ (69)

  * (b)

There exists a constant $\tilde{c}>0$ such that

$\mathbb{E}\|x(t)-x^{\dagger}\|^{2}\leq\tilde{c}\phi(t)\quad\text{for all}\
t>0.$ (70)

###### Proof.

From $\varphi\in\mathcal{S}^{g}_{\zeta}$, we have $\phi(\gamma t)\leq
g(1/\gamma)\phi(t)$, which implies that $\tilde{\phi}(\gamma
t)\leq\sqrt{\gamma^{-1}g(1/\gamma)}\tilde{\phi}(t)$, and so, by setting
$\tilde{g}(\gamma)=\sqrt{\gamma^{-1}g(1/\gamma)},\delta=\tilde{\phi}(t)$ and
$\tilde{\gamma}=\tilde{g}(\gamma)$, we get

$\tilde{g}^{-1}(\tilde{\gamma})\tilde{\phi}^{-1}(\delta)\geq\tilde{\phi}^{-1}(\tilde{\gamma}\delta).$

Thus, we have

$\psi(\tilde{\gamma}\delta)=\tilde{\gamma}^{2}\delta^{2}\tilde{\phi}^{-1}(\tilde{\gamma}\delta)\leq\tilde{\gamma}^{2}\delta^{2}\tilde{g}^{-1}(\tilde{\gamma})\tilde{\phi}^{-1}(\delta)=h(\tilde{\gamma})\psi(\delta),$
(71)

where $h(\tilde{\gamma})=\tilde{\gamma}^{2}\tilde{g}^{-1}(\tilde{\gamma})$.

In the case where $\mathbb{E}\|x(t)-x^{\dagger}\|=0$ for all $t>C_{6}$ for
some $C_{6}>0$, the inequality (70) is trivially fulfilled for some
$\tilde{c}>0$. Moreover, we know from Lemma 2 that then the inequality (61)
holds, which implies the inequality (69) for some constant $c>0$, since we
have, according to the definition of the function $\psi$ and the decreasing
property of $\tilde{\phi}(t)$, $\psi(\delta)\geq a\delta^{2}$ for all
$\delta\in(0,\delta_{0})$ for some constants $a>0$ and $\delta_{0}>0$.

Thus, we may assume that $\mathbb{E}\|x(t)-x^{\dagger}\|>0$ for all $t>0$.

Let (70) hold. For arbitrary $\delta>0$, we use the parameter $t_{\delta}$
defined in (63). Then, the inequality (70) implies that

$\delta^{2}t_{\delta}\leq\tilde{c}\phi(t_{\delta}).$

Consequently,

$\tilde{\phi}^{-1}\Big{(}\frac{\delta}{\sqrt{\tilde{c}}}\Big{)}\geq
t_{\delta},$

and therefore, using the upper bound in (64) obtained in Lemma 3, we find with
(71) that

$\sup_{\tilde{y}\in\tilde{B}_{\delta}(y)}\inf_{t>0}\mathbb{E}\|x(t;\tilde{y})-x^{\dagger}\|^{2}\leq
C_{7}\delta^{2}t_{\delta}\leq
C_{7}\delta^{2}\tilde{\phi}^{-1}\left(\frac{\delta}{\sqrt{\tilde{c}}}\right)=C_{7}\tilde{c}\psi\left(\frac{\delta}{\sqrt{\tilde{c}}}\right)\leq
C_{7}\tilde{c}h(\tilde{c}^{-1/2})\psi(\delta),$

which is the estimate (69) with $c=C_{7}\tilde{c}h(\tilde{c}^{-1/2})$.

Conversely, if (69) holds, we choose an arbitrary $\delta>0$ such that
$t_{\delta}$ is defined by (63). Then, we can use the lower bound in (64) of
Lemma 3 to obtain, from the condition (69),

$C_{8}\delta^{2}t_{\delta}\leq c\psi(\delta).$

Thus, from the definition of $\psi$, we have
$\frac{C_{8}}{c}t_{\delta}\leq\tilde{\phi}^{-1}(\delta)$, and consequently

$\left(\frac{C_{8}}{c}t_{\delta}\right)^{-1}\phi\left(\frac{C_{8}}{c}t_{\delta}\right)=\left[\tilde{\phi}\left(\frac{C_{8}}{c}t_{\delta}\right)\right]^{2}\geq\delta^{2},$

So, finally, we obtain, with the inequality $\phi(\gamma t)\leq
g(1/\gamma)\phi(t)$,

$\mathbb{E}\|x(t_{\delta})-x^{\dagger}\|^{2}=\delta^{2}t_{\delta}\leq\left(\frac{C_{8}}{c}t_{\delta}\right)^{-1}\phi\left(\frac{C_{8}}{c}t_{\delta}\right)t_{\delta}=\frac{c}{C_{8}}\phi\Big{(}\frac{C_{8}}{c}t_{\delta}\Big{)}\leq\frac{c}{C_{8}}g\left(\frac{c}{C_{8}}\right)\phi(t_{\delta}),$

and, since this holds for every $\delta$, (70) holds with
$\hat{c}=\frac{c}{C_{8}}g(\frac{c}{C_{8}})$.

∎

##  3 Numerical realization of SAR

In order to use SAR in practice, we need to discretize the stochastic flow (5)
for the introduced artificial time variable. Numerous algorithms have been
proposed (see, e.g., [12, 17]) for accurate numerical approximation of
stochastic differential equations. Here, we propose two of the simplest
numerical approaches: the Euler method,

$\displaystyle x^{\delta}_{k+1}=x^{\delta}_{k}+\Delta
tA^{*}(y^{\delta}-Ax^{\delta}_{k})+f_{k}\Delta B_{k},\quad
x^{\delta}_{0}=x_{0},\quad\Delta t\in(0,2/\|A\|^{2}),$ (72)

and the exponential Euler method,

$x^{\delta}_{k+1}=e^{-A^{*}A\Delta
t}\big{[}x_{k}^{\delta}+A^{*}y^{\delta}\Delta t+f_{k}\Delta B_{k}\big{]},$
(73)

Where $f_{k}=f(t_{k})$, $\Delta t$ is the size of the uniform time step, and

$\Delta
B_{k}=B(t_{k+1})-B(t_{k})=\sum_{j=1}^{\infty}\sqrt{q_{j}}u_{j}\big{(}\beta_{j}(t_{k+1})-\beta_{j}(t_{k})\big{)}=\sqrt{\Delta
t}\sum_{j=1}^{\infty}\sqrt{q_{j}}u_{j}\xi_{j}^{k},$

with $\\{\xi_{j}^{k}\\}_{j,k}$ being a sequence of independent
${\mathcal{N}}(0,1)$-random variables.

######  Proposition 5.

Let $f\in L^{\infty}({\mathbb{R}}_{+})$ be globally Lipschitz continuous with
Lipschitz constant $\|f\|_{Lip}$. The numerical solutions of (72) and (73)
converge to the exact solution (7) of (5) with order one in the mean-square
sense, i.e. for $k=1,2,\cdots$,

${\mathbb{E}}\|x^{\delta}(t_{k})-x^{\delta}_{k}\|^{2}\leq C\Delta t^{2},$ (74)

where $C:=C(t_{k},{\rm tr}(Q(A^{*}A)^{-1}),\|A\|,\|A^{*}y^{\delta}\|,f)$.

###### Proof.

It follows from (72) and (73) that

$x^{\delta}_{k}=(I-A^{*}A\Delta t)^{k}x_{0}+\sum_{j=0}^{k-1}(I-A^{*}A\Delta
t)^{k-1-j}A^{*}y^{\delta}\Delta t+\sum_{j=0}^{k-1}(I-A^{*}A\Delta
t)^{k-1-j}f_{j}\Delta B_{j}$

and

$x^{\delta}_{k}=e^{-A^{*}At_{k}}x_{0}+\sum_{j=0}^{k-1}e^{-A^{*}A(t_{k}-t_{j})}A^{*}y^{\delta}\Delta
t+\sum_{j=0}^{k-1}e^{-A^{*}A(t_{k}-t_{j})}f_{j}\Delta B_{j}.$

We recall that the exact solution (7) of (5) at time $t_{k}$ is

$x^{\delta}(t_{k})=e^{-A^{*}At_{k}}x_{0}+\int_{0}^{t_{k}}e^{-A^{*}A(t_{k}-s)}A^{*}y^{\delta}ds+\int_{0}^{t_{k}}e^{-A^{*}A(t_{k}-s)}f(s)dB_{s}.$

Subtracting the above two equations and taking the
${\mathbb{E}}\|\cdot\|^{2}$-norm, we get the estimate of the error for the
Euler method (72),

$\displaystyle{\mathbb{E}}\|x^{\delta}(t_{k})-x^{\delta}_{k}\|^{2}$
$\displaystyle\leq 3\|(e^{-A^{*}At_{k}}-(I-A^{*}A\Delta t)^{k})x_{0}\|^{2}$
$\displaystyle\quad+3{\mathbb{E}}\Big{\|}\sum_{j=0}^{k-1}\int_{t_{j}}^{t_{j+1}}\big{(}e^{-A^{*}A(t_{k}-s)}-(I-A^{*}A\Delta
t)^{k-1-j}\big{)}A^{*}y^{\delta}ds\Big{\|}^{2}$
$\displaystyle\quad+3{\mathbb{E}}\Big{\|}\sum_{j=0}^{k-1}\int_{t_{j}}^{t_{j+1}}\big{(}e^{-A^{*}A(t_{k}-s)}f(s)-(I-A^{*}A\Delta
t)^{k-1-j}f_{j}\big{)}dB(s)\Big{\|}^{2}$ $\displaystyle\leq
3\|(e^{-A^{*}At_{k}}-(I-A^{*}A\Delta t)^{k})x_{0}\|^{2}$
$\displaystyle\quad+3k\Delta
t\sum_{j=0}^{k-1}\int_{t_{j}}^{t_{j+1}}\big{\|}\big{(}e^{-A^{*}A(t_{k}-s)}-(I-A^{*}A\Delta
t)^{k-1-j}\big{)}A^{*}y^{\delta}\big{\|}^{2}ds$
$\displaystyle\quad+3{\mathbb{E}}\sum_{j=0}^{k-1}\int_{t_{j}}^{t_{j+1}}\big{\|}\big{(}e^{-A^{*}A(t_{k}-s)}f(s)-(I-A^{*}A\Delta
t)^{k-1-j}f_{j}\big{)}Q^{\frac{1}{2}}\big{\|}_{HS(\mathcal{X};\mathcal{X})}^{2}ds,$

and the estimate of the error for the exponential Euler method (73),

$\displaystyle{\mathbb{E}}\|x^{\delta}(t_{k})-x^{\delta}_{k}\|^{2}$
$\displaystyle\leq
2{\mathbb{E}}\Big{\|}\sum_{j=0}^{k-1}\int_{t_{j}}^{t_{j+1}}\big{(}e^{-A^{*}A(t_{k}-s)}-e^{-A^{*}A(t_{k}-t_{j})}\big{)}A^{*}y^{\delta}ds\Big{\|}^{2}$
$\displaystyle\quad+2{\mathbb{E}}\Big{\|}\sum_{j=0}^{k-1}\int_{t_{j}}^{t_{j+1}}\big{(}e^{-A^{*}A(t_{k}-s)}f(s)-e^{-A^{*}A(t_{k}-t_{j})}f_{j}\big{)}dB(s)\Big{\|}^{2}$
$\displaystyle\leq 2k\Delta
t\sum_{j=0}^{k-1}\int_{t_{j}}^{t_{j+1}}\big{\|}\big{(}e^{-A^{*}A(t_{k}-s)}-e^{-A^{*}A(t_{k}-t_{j})}\big{)}A^{*}y^{\delta}\big{\|}^{2}ds$
$\displaystyle\quad+2{\mathbb{E}}\sum_{j=0}^{k-1}\int_{t_{j}}^{t_{j+1}}\big{\|}\big{(}e^{-A^{*}A(t_{k}-s)}f(s)-e^{-A^{*}A(t_{k}-t_{j})}f_{j}\big{)}Q^{\frac{1}{2}}\big{\|}_{HS(\mathcal{X};\mathcal{X})}^{2}ds.$

Note that

$\displaystyle\|(e^{-A^{*}At_{k}}-(I-A^{*}A\Delta
t)^{k})\|_{{\mathcal{L}}({\mathcal{X}};{\mathcal{X}})}$
$\displaystyle=\Big{\|}\sum_{\ell=0}^{k-1}e^{-A^{*}At_{k-1-\ell}}(e^{-A^{*}A\Delta
t}-(I-A^{*}A\Delta t))(I-A^{*}A\Delta
t)^{\ell}\Big{\|}_{{\mathcal{L}}({\mathcal{X}};{\mathcal{X}})}$
$\displaystyle\leq k\|e^{-A^{*}A\Delta t}-(I-A^{*}A\Delta
t)\|_{{\mathcal{L}}({\mathcal{X}};{\mathcal{X}})}$ $\displaystyle\leq
k\sup_{0<\lambda\leq\|A\|^{2}}|e^{-\lambda\Delta t}-(1-\lambda\Delta t)|\leq
t_{k}\|A\|^{2}\Delta t,$

and, for $0<s-t\leq\Delta t$,

$\displaystyle\big{\|}e^{-A^{*}At}-e^{-A^{*}As}\big{\|}_{{\mathcal{L}}({\mathcal{X}};{\mathcal{X}})}\leq\sup_{\lambda\in(0,\|A\|^{2}]}e^{-\lambda
t}\big{(}1-e^{-\lambda(s-t)}\big{)}\leq\|A\|^{2}\Delta t.$

Hence, for the Euler method (72),

$\displaystyle\big{\|}\big{(}e^{-A^{*}A(t_{k}-s)}-(I-A^{*}A\Delta
t)^{k-1-j}\big{)}A^{*}y^{\delta}\big{\|}$
$\displaystyle\leq\big{\|}\big{(}e^{-A^{*}A(t_{k}-s)}-e^{-A^{*}A(t_{k}-t_{j+1})}\big{)}A^{*}y^{\delta}\big{\|}$
$\displaystyle+\big{\|}\big{(}e^{-A^{*}A(t_{k}-t_{j+1})}-(I-A^{*}A\Delta
t)^{k-1-j}\big{)}A^{*}y^{\delta}\big{\|}$
$\displaystyle\leq\|A^{*}y^{\delta}\|\|A\|^{2}(1+t_{k})\Delta t,$

and

$\displaystyle\big{\|}\big{(}e^{-A^{*}A(t_{k}-s)}f(s)-(I-A^{*}A\Delta
t)^{k-1-j}f_{j}\big{)}Q^{\frac{1}{2}}\big{\|}_{HS(\mathcal{X};\mathcal{X})}^{2}$
$\displaystyle\leq{\rm
tr}(Q(A^{*}A)^{-1})\|A\|^{2}\Big{[}\|(e^{-A^{*}A(t_{k}-s)}-(I-A^{*}A\Delta
t)^{k-1-j})f(s)\|^{2}$ $\displaystyle\quad\quad+\|(I-A^{*}A\Delta
t)^{k-1-j}(f(s)-f(t_{j}))\|^{2}\Big{]}$ $\displaystyle\leq{\rm
tr}(Q(A^{*}A)^{-1})\|A\|^{2}\Big{[}t_{k-1-j}\|f\|^{2}_{L^{\infty}}\|A\|^{4}+\|f\|_{Lip}^{2}\Big{]}(\Delta
t)^{2},$

which lead to the assertion (74). Similarly, for the exponential Euler method,
we have (73),

$\displaystyle\big{\|}\big{(}e^{-A^{*}A(t_{k}-s)}-e^{-A^{*}A(t_{k}-t_{j})}\big{)}A^{*}y^{\delta}\big{\|}$
$\displaystyle\leq\big{\|}e^{-A^{*}A(t_{k}-s)}-e^{-A^{*}A(t_{k}-t_{j})}\big{\|}_{{\mathcal{L}}({\mathcal{X}};{\mathcal{X}})}\|A^{*}y^{\delta}\|$
$\displaystyle\leq\|A^{*}y^{\delta}\|\sup_{\lambda\in(0,\|A\|^{2}]}e^{-\lambda(t_{k}-s)}\big{(}1-e^{-\lambda(s-t_{j})}\big{)}$
$\displaystyle\leq\|A^{*}y^{\delta}\|\|A\|^{2}\Delta t,$

and

$\displaystyle\big{\|}\big{(}e^{-A^{*}A(t_{k}-s)}f(s)-e^{-A^{*}A(t_{k}-t_{j})}f_{j}\big{)}Q^{\frac{1}{2}}\big{\|}_{HS(\mathcal{X};\mathcal{X})}^{2}$
$\displaystyle\leq{\rm
tr}(Q(A^{*}A)^{-1})\|A\|^{2}\Big{[}\|(e^{-A^{*}A(t_{k}-s)}-e^{-A^{*}A(t_{k}-t_{j})})f(s)\|^{2}+\|e^{-A^{*}A(t_{k}-t_{j})}(f(s)-f(t_{j}))\|^{2}\Big{]}$
$\displaystyle\leq{\rm
tr}(Q(A^{*}A)^{-1})\|A\|^{2}\Big{[}\|f\|^{2}_{L^{\infty}}\|A\|^{4}+\|f\|_{Lip}^{2}\Big{]}(\Delta
t)^{2},$

which yield to the assertion (74) for the exponential Euler method.

∎

##  4 Numerical experiments

This section presents some numerical experiments to illustrate the numerical
behavior of SAR. In the first example of abstract integral equations, we
demonstrate that, in contrast to the conventional deterministic regularization
methods, SAR provides the uncertainty quantification of the regularized
solution. In the second example, we further show that the uncertainty
quantification of the estimated physical quantity can reveal and explicate the
hidden information about the real-world inverse problem.

###  4.1 A toy example showing the uncertainty capability of SAR

Our first group of examples are based on the following integral equation:

$Ax(s):=\int^{1}_{0}K(s,t)x(t)dt=y(s),\quad K(s,t)=s(1-t)\chi_{s\leq
t}+t(1-s)\chi_{s>t}.$ (75)

If we choose $\mathcal{X}=\mathcal{Y}=L^{2}[0,1]$, the operator $A$ is
compact, selfadjoint, and injective. It is well known that the integral
equation (79) has a solution $x=-y^{\prime\prime}$ if $y\in H^{2}[0,1]\cap
H^{1}_{0}[0,1]$. Furthermore, using the interpolation theory (see, e.g.,
[16]), it is not difficult to show that for $4p-1/2\not\in\mathbb{N}$

$\displaystyle\mathcal{R}((A^{*}A)^{p})=\left\\{x\in
H^{4p}[0,1]:~{}x^{2l}(0)=x^{2l}(1)=0,~{}l=0,1,\cdots,\lfloor
2p-1/4\rfloor\right\\}.$

If we choose $y(s)=s^{4}(1-s)^{3}$, $x^{\dagger}=-6t^{2}(1-t)(2-8t+7t^{2})$,
and $x^{\dagger}\in R((A^{*}A)^{p})$ for all $p<5/8$. The numerical results
are displayed in Figures 1 and 2, where, in addition to a deterministic
approximation solution (the expectation of SAR), the 85% and 70% confidence
intervals are also provided. In the numerical application, we recommend the
use of a small confidence interval when the data contain large noise, since in
such a setting the exact solution function can be located in the corresponding
confidence interval with high probability.

Figure 1: The expectation of SAR and the 85% confidence interval for problem
(79) with noise level $\delta=1\%$. Other parameters: the size of
discretization is 100, $\tau=1.1$, $\Delta t=0.1$. Figure 2: The expectation
of SAR and the 70% confidence interval for problem (79) with noise level
$\delta=5\%$. Other parameters: the size of discretization is 100, $\tau=1.1$,
$\Delta t=0.1$.

###  4.2 Biosensor tomography with real data

In this section, the developed SAR is applied to a biosensor tomography
problem, i.e. finding the interaction information from biosensor data. We
demonstrate that, in contrast to the conventional deterministic regularization
methods, SAR provides the uncertainty quantification of the regularized
solution, which can be further used to reveal and explicate the hidden
information about the physical quantity of interest. To this end, we first
briefly review the mechanism of biosensor tomography.

####  4.2.1 The mathematical model

For the considered biosensors, the analyte binds to an antibody that is
attached to a surface. Then, a physical method (e.g. evanescent wave, surface
plasmon resonance) is used to measure the surface concentration of the
antibody-analyte complex. To better understand the construction of the
mathematical model, we consider the three steps of the binding process during
the modeling.

Step 1: “1-to-1” meta kinetic model. For each interaction we have

$[A]+[L]\autoleftrightharpoons{k_{a}}{k_{d}}[AL],$

where $k_{a}$ and $k_{d}$ denote the association and dissociation rate
constants, and $[A]$, $[L]$, and $[AL]$ represent the concentrations of the
analyte, ligand, and complex, respectively. In our simplified model, the
analyte $A$ is injected and flushed over the surface in such a way that the
concentration $[A]$ can be assumed to be constant during the study. The amount
of free ligand decreases with time according to $[L](t)=[L](0)-[AL](t)$. We
assume that the sensor response $R$ is proportional to the complex
concentration $[AL](t)$, i.e. $R(t)=C_{b}\cdot[AL](t)$, where $C_{b}$ is a
constant. Let $R_{max}=C_{b}\cdot[L](0)$. Then, if the mass-transfer kinetics
are extremely fast, the rate of complex formation satisfies the following
dynamical equation:

$\frac{dR(t)}{dt}=k_{a}\cdot[A](t)\cdot\left(R_{max}-R(t)\right)-k_{d}\cdot
R(t).$ (76)

Setting $C=[A]$ (as mentioned above, $C$ is a fixed number) and $R(t_{0})=0$,
the solution to (76) is

$R(t)=R_{max}\cdot\frac{k_{a}C}{k_{d}+k_{a}C}\cdot\left(1-e^{(k_{d}+k_{a}C)(t-t_{0})}\right).$
(77)

Step 2: “m-to-n” kinetic model. We now assume that there are $m$ analytes and
$n$ binding sites on the biosensor surface and first-order kinetics. Let
$(k_{a,i},k_{d,j})$ denote the pair of association and dissociation constants
for the interaction between the $i$th analyte and $j$th binding site. Let
$R_{i,j}(t)$ be the response at time $t$ of a complex with association
constant $k_{a,i}$ and dissociation constant $k_{d,j}$. Then, according to
(77), we have

$R_{i,j}(t)=\left\\{\begin{array}[]{>{\displaystyle}l>{\displaystyle}l}&0,\qquad
t\leq t_{0}+\Delta t,\\\
&R^{max}_{i,j}\frac{k_{a,i}C}{k_{d,j}+k_{a,i}C}\left(1-e^{-(k_{d,j}+k_{a,i}C)(t-t_{0})}\right),\\\
&\qquad\qquad\qquad\qquad\qquad t_{0}+\Delta t<t\leq t_{0}+t_{inj}+\Delta
t,\\\
&R^{max}_{i,j}\frac{k_{a,i}C}{k_{d,j}+k_{a,i}C}\left(1-e^{-(k_{d,j}+k_{a,i}C)t_{inj}}\right)e^{-k_{d,j}(t-t_{0}-t_{inj})},\\\
&\qquad\qquad\qquad\qquad\qquad\qquad\qquad t>t_{0}+t_{inj}+\Delta
t,\end{array}\right.$ (78)

where constant $C$ represents the concentration of the analyte, $t_{0}$ is the
time at which the injection of the analyte begins, and $t_{inj}$ is the
injection time. The adjustment parameter $\Delta t$ is a time delay that
accounts for the fact that it usually takes some time for the detector to
respond to the injection. In our software, $\Delta t$ can be automatically
determined; it can be understood as a preconditioning of a real piece of data.
Constant $R^{max}_{i,j}$ is the total surface-binding capacity, corresponding
to association and dissociation constants $k_{a,i}$ and $k_{d,j}$, i.e. the
detector response when every binding site on the biosensor surface has formed
a complex with the analyte.

Step 3: Continuous kinetic model. We now use the functions $R_{i,j}$ to make
an approximation of the measured sensorgrams $R_{obs}$. By employing the
principle of superposition, the total measured response, $R_{obs}$, of a
system can be written as a linear combination of some individual responses,
namely $R_{obs}=\sum^{m,n}_{i,j=1}R_{i,j}$. If we let $m,n\to+\infty$ in the
above equation, we obtain the final mathematical model of biosensor
tomography:

$R_{obs}(t;C)=\int_{\Omega}K(t,C;k_{a},k_{d})x(k_{a},k_{d})dk_{a}dk_{d},\quad(k_{a},k_{d})\in\Omega,$
(79)

where $\Omega\subset\mathbb{R}^{2}$ is the domain of rate constants that is of
interest, and the kernel function $K(\cdot)$ is defined as

$K(t,C;k_{a},k_{d})=\left\\{\begin{array}[]{>{\displaystyle}l>{\displaystyle}l}&0,\qquad
t\leq t_{0}+\Delta t,\\\
&\frac{k_{a}C}{k_{d}+k_{a}C}\left(1-e^{-(k_{d}+k_{a}C)(t-t_{0})}\right),\\\
&\qquad\qquad\qquad\qquad\qquad t_{0}+\Delta t<t\leq t_{0}+t_{inj}+\Delta
t,\\\
&\frac{k_{a}C}{k_{d}+k_{a}C}\left(1-e^{-(k_{d}+k_{a}C)t_{inj}}\right)e^{-k_{d}(t-t_{0}-t_{inj})},\\\
&\qquad\qquad\qquad\qquad\qquad\qquad\qquad t>t_{0}+t_{inj}+\Delta
t,\end{array}\right.$ (80)

The function $x(k_{a},k_{d})$, which is the generalization of the total
surface-binding capacity $\\{R^{max}_{i,j}\\}$, is known as the (continuous)
rate-constant map. In the _rate-constant-map theory_ of chemical reactions,
the local and global maximums of a rate-constant map reflect the interaction
information about a chemical system consisting of various molecules:

  * •

The global maximum of a rate-constant map reflects the principal interaction
among molecules, and the coordinate of the global maximum represents the value
of association and dissociation rate constants of the principle interaction.

  * •

The local maximums of a rate-constant map reflect the minor interactions among
molecules. The number of local maximums shows the number of minor interactions
in the considered environment. The corresponding coordinates of the local
maximums represent the value of association and dissociation rate constants of
the corresponding minor interactions.

In sum, according to the rate-constant-map theory, the biosensor tomography
can be solved via the following two steps:

  * •

Given noisy biosensor data $R_{obs}(t;C)$, find the approximate rate-constant
map $x(k_{a},k_{d})$.

  * •

Determine the interaction information (i.e. the number of interactions and the
corresponding association and dissociation rate constants) from the estimated
rate-constant map $x(k_{a},k_{d})$.

In the next subsection, we show that SAR can substantially improve the second
step of biosensor tomography through the use of the statistical property of
SAR.

####  4.2.2 Simulation study

We tested SAR for biosensor tomography with real experimental data – the
parathyroid hormone (PTH). In the experiment, the human PTH1R receptor was
immobilized on a LNB-carboxyl biosensor chip using amine coupling according to
the manufacturer’s instructions. Using the flow rate 25 $\mu L/min$ at 20.0∘C,
we did 35 $\mu L$ injections of the peptide PTH(1-34) at 9 concentration
levels from 1214 nM to 14571 nM; see the solid lines in Figure 3. The
sensorgrams were measured with a QCM biosensor Attana Cell 200 (Attana AB,
Stockholm, Sweden) instrument.

A standard adaptive linear finite-element method was adopted to discretize the
two-dimensional integral equation (79). The initial triangulation was
uniformly distributed in the log-scale domain
$(\log_{10}(k_{d}),\log_{10}(k_{a}))\in[-4,0]\times[3,7]$ with $10\times
10=200$ node points. Then, a standard adaptive approach (see, e.g., [15, 27])
was employed to achieve better accuracy with minimum degrees of freedom. More
precisely, we first solved integral equation (79) to obtain the solution on
the current triangulation. The error was then estimated using the solution,
and used to mark a set of triangles that were to be refined. The triangles
were refined in such a way as to keep two of the most important properties of
the triangulations: shape regularity and conformity. Our algorithm stopped at
the 6-th iteration with 600 nodes and 1149 triangles, as shown in Figure 3.
The total running time was 13.6 seconds 222All the computations were performed
on a dual-core personal computer with 1.00 GB RAM and MatLab version R2019b..
The estimated rate-constant-mean map and corresponding intensity map and
contour map are shown in Figures 4 and 5, respectively. Figure 3 also displays
a comparison between the experimental data (solid line) and simulated response
curves (dashed line), which was obtained by solving the forward problem with
the expectation of SAR. The results, with approximately 95.7% overlap,
indicate that the estimated rate-constant-mean map can be used as the real
rate-constant map in a deterministic model.

Figure 3: The experimental data (solid line) and simulated response (dashed
line) with different concentrations of PTH at a temperature of 20
${}^{\circ}C$. The relative residual error equals 0.013. Figure 4: The
estimated rate-constant-mean map.

Figure 5: (a) The intensity map. (b) The contour map.

The results obtained so far (as shown in Figures 3-5) are similar to those
obtained from conventional deterministic regularization approaches, e.g. the
Landweber iteration (3) and the Tikhonov regularization proposed in [6]. Note
that the approximate solution from FAR is a random variable, which provides
uncertainty quantification in the estimation that would be useful for a real
problem. More precisely, for the considered biosensor tomography, we can
consider interaction-information maps from the high-order moments of FAR;
Figures 6 and 7 show the results of maps based on second moments. From these
figures, we can easily conclude that three interactions of the considered
chemical systems exist, which theoretically confirms the empirical assumption
of chemists and biologists that at least two interactions of the PTH system
exist. Furthermore, from the maps based on second moments, we can easily
determine the corresponding detailed information about three interactions: a
principal interaction of the considered biosensor system exists. The
corresponding dissociation and association rate constants are
$(k_{d},k_{a})=(10^{-0.9},10^{4.4})$. Moreover, two minor interactions exist,
and the corresponding dissociation and association rate constants are
$(k_{d},k_{a})=(10^{-3.5},10^{3.9})$ and $(k_{d},k_{a})=(10^{-2.9},10^{2.6})$,
respectively.

Figure 6: The estimated second-moment map.

Figure 7: (a) The second-moment intensity map. (b) The second-moment contour
map.

It should be noted that these results are similar to the results in [30],
which uses the Bayesian model333The Bayesian model is a statistical model,
while our model (operator equation (79)) is deterministic. with Markov Chain
Monte Carlo and variational Bayesian approach. The computational cost of the
Bayesian-based approach is clearly very high (to obtain a convergent result
takes hours for a problem of similar size to those above). Finally, for the
considered real data problem, the maps based on high-order moments, i.e. the
third and fourth moment maps, present very similar sharps to the second-order
maps, and hence we did not show the corresponding numerical results. However,
for other biosensor systems, we remain inclined to suggest that this analysis
be performed, since it is possible that the sharps of the maps based on high-
order moments behave differently, which will be interesting to study when
analyzing the advanced properties of interaction information about biosensor
systems.

##  5 Conclusions

We develop a statistical approach – SAR for deterministic linear inverse
problems. In this paper, a complete regularization theory of SAR has been
presented and analyzed. In particular, it has been shown that SAR yields an
optimal regularization method with regard to mean-square convergence. Compared
with the conventional deterministic regularization method, SAR provides the
uncertainty quantification of the estimated quantity, which is described by a
deterministic forward model. A real data application of SAR for biosensor
tomography showed that the uncertainty quantification of SAR can reveal and
explicate the interaction information about a chemical system, which is hidden
behind the original deterministic forward model. We therefore believe that SAR
will be a useful tool for studying biosensor interactions and other real-world
inverse problems with deterministic ill-posed forward models.

##  6 Acknowledgement

This work of Y. Zhang is supported by the Beijing Natural Science Foundation
(Key project No. Z210001), the National Natural Science Foundation of China
(No. 12171036), the Shenzhen Stable Support Fund for College Researchers
20200827173701001 and the Guangdong Fundamental and Applied Research Fund (No.
2019A1515110971). The author C. Chen is supported by the National key R&D
Program of China (No. 2020YFA0713701), the National Natural Science Foundation
of China (No. 11871068, No. 12022118) and by Youth Innovation Promotion
Association CAS.

## References

  * [1] V. Albani, P. Elbau, M.V. de Hoop, and O. Scherzer.  Optimal convergence rates results for linear inverse problems in hilbert spaces.  Numerical Functional Analysis and Optimization, 37:521–540, 2016\.
  * [2] R. Boţ, G. Dong, P. Elbau, and O. Scherzer.  Convergence rates of first- and higher-order dynamics for solving linear ill-posed problems.  Foundations of Computational Mathematics, pages https://doi.org/10.1007/s10208–021–09536–6, 2021.
  * [3] L. Bottou, F.E. Curtis, and J. Nocedal.  Optimization methods for large-scale machine learning.  SIAM Review, 60:223–311, 2018.
  * [4] G. Da Prato and J. Zabczyk.  Stochastic equations in infinite dimensions, volume 152 of Encyclopedia of Mathematics and its Applications.  Cambridge University Press, Cambridge, second edition, 2014.
  * [5] W. E.  A proposal on machine learning via dynamical systems.  Communications in Mathematics and Statistics, 5:1–11, 2017.
  * [6] W. E, C. Ma, and L. Wu.  Machine learning from a continuous viewpoint, i.  Science China Mathematics, 63:2233–2266, 2020.
  * [7] J. Flemming.  Variational smoothness assumptions in convergence rate theory – an overview.  Journal of Inverse and Ill-posed Problems, 21:395–409, 2013.
  * [8] R. Gong, B. Hofmann, and Y. Zhang.  A new class of accelerated regularization methods, with application to bioluminescence tomography.  Inverse Problems, 36:055013, 2020.
  * [9] B. Hofmann, B. Kaltenbacher, C. Pöschl, and O. Scherzer.  A convergence rates result for tikhonov regularization in banach spaces with non-smooth operators.  Inverse Problems, 23:987–1010, 2007.
  * [10] B. Hofmann and M. Yamamoto.  On the interplay of source conditions and variational inequalities for nonlinear ill-posed problems.  Applicable Analysis, 89:1705–1727, 2010.
  * [11] T. Jahn and B. Jin.  On the discrepancy principle for stochastic gradient descent.  Inverse Problems, 36:095009, 2020.
  * [12] Arnulf Jentzen and Peter E. Kloeden.  Overcoming the order barrier in the numerical approximation of stochastic partial differential equations with additive space-time noise.  Proceedings of the Royal Society A: Mathematical,Physical and Engineering Sciences, 465(2102):649–667, 2009.
  * [13] B. Jin and X. Lu.  On the regularizing property of stochastic gradient descent.  Inverse Problems, 35:015004, 2019.
  * [14] B. Jin, Z Zhou, and J. Zou.  On the convergence of stochastic gradient descent for nonlinear ill-posed problems.  SIAM Journal on Optimization, 30:1421–1450, 2020.
  * [15] N. Koshev and L. Beilina.  An adaptive finite element method for fredholm integral equations of the first kind and its verificationon experimental data.  Central European Journal of Mathematics, 11:1489–1509, 2013.
  * [16] J. Lions and E. Magenes.  Non-Homogeneous Boundary Value Problems and Applications, Volumes I.  Berlin: Springer, 1972.
  * [17] Gabriel J. Lord, Catherine E. Powell, and Tony Shardlow.  An introduction to computational stochastic PDEs.  Cambridge Texts in Applied Mathematics. Cambridge University Press, New York, 2014.
  * [18] S. Lu and P. Mathé.  Stochastic gradient descent for linear inverse problems in Hilbert spaces. doi: https://doi.org/10.1090/mcom/3714.  Mathematics of Computation, 2021.
  * [19] S. Lu, P. Niu, and F. Werner.  On the asymptotical regularization for linear inverse problems in presence of white noise.  SIAM/ASA J. Uncertainty Quantification, 9:1–28, 2021.
  * [20] P. Mathé and S. V. Pereverzev.  Geometry of linear ill-posed problems in variable Hilbert scales.  Inverse Problems, 19(3):789–803, 2003.
  * [21] M.Z. Nashed.  A new approach to classification and regularization of ill-posed operator equations.  In H.W. Engl and C.W. Groetsch, editors, Inverse and Ill-Posed Problems, pages 53–75. Academic Press, 1987.
  * [22] A. Rieder.  Runge-Kutta integrators yield optimal regularization schemes.  Inverse Problems, 21:453–471, 2005.
  * [23] H. Robbins and S. Monro.  A stochastic approximation method.  The Annals of Mathematical Statistics, 22:400–407, 1951.
  * [24] E. Schock.  Approximate solution of ill-posed equations: arbitrarily slow convergence vs. superconvergence.  Constructive methods for the practical treatment of integral equations, 73:234–243, 1985.
  * [25] U. Tautenhahn.  On the asymptotical regularization of nonlinear ill-posed problems.  Inverse Problems, 10:1405–1418, 1994.
  * [26] G. Vainikko and A. Veretennikov.  Iteration Procedures in Ill-Posed Problems.  Moscow: Nauka (In Russian), 1986.
  * [27] Y. Zhang, P. Forssen, T. Fornstedt, M. Gulliksson, and X. Dai.  An adaptive regularization algorithm for recovering the rate constant distribution from biosensor data.  Inverse Problems in Science and Engineering, 26:1464–1489, 2018\.
  * [28] Y. Zhang and B. Hofmann.  On fractional asymptotical regularization of linear ill-posed problems in Hilbert spaces.  Fractional Calculus and Applied Analysis, 22:699–721, 2019.
  * [29] Y. Zhang and B. Hofmann.  On the second order asymptotical regularization of linear ill-posed inverse problems.  Applicable Analysis, 99:1000–1025, 2020.
  * [30] Y. Zhang, Z. Yao, P. Forssen, and T. Fornstedt.  Estimating the rate constant from biosensor data via an adaptive variational Bayesian approach.  Annals of Applied Statistics, 13:2011–2042, 2019.
  * [31] Y. Zhao, P. Mathé, and S. Lu.  Convergence analysis of asymptotical regularization and Runge-Kutta integrators for linear inverse problems under variational source conditions.  The CSIAM Transactions on Applied Mathematics, 1:693–714, 2020\.