# Dual gradient method for ill-posed problems using multiple repeated
measurement data

Qinian Jin  Mathematical Sciences Institute, Australian National University,
Canberra, ACT 2601, Australia<EMAIL_ADDRESS>and  Wei Wang  College of
Data Science, Jiaxing University, Jiaxing, Zhejiang Province, 314001, China
<EMAIL_ADDRESS>
(Date: 24 June 2022)

###### Abstract.

We consider determining ${\mathcal{R}}$-minimizing solutions of linear ill-
posed problems $Ax=y$, where $A:{\mathscr{X}}\to{\mathscr{Y}}$ is a bounded
linear operator from a Banach space ${\mathscr{X}}$ to a Hilbert space
${\mathscr{Y}}$ and ${\mathcal{R}}:{\mathscr{X}}\to[0,\infty]$ is a proper
strongly convex penalty function. Assuming that multiple repeated independent
identically distributed unbiased data of $y$ are available, we consider a dual
gradient method to reconstruct the ${\mathcal{R}}$-minimizing solution using
the average of these data. By terminating the method by either an a priori
stopping rule or a statistical variant of the discrepancy principle, we
provide the convergence analysis and derive convergence rates when the sought
solution satisfies certain variational source conditions. Various numerical
results are reported to test the performance of the method.

##  1\. Introduction

Inverse problems arise from many practical applications in science and
engineering when one aims at determining unobservable causes from observed
effects. Due to their significance in wide range of applications, inverse
problems have received tremendous attention in the literature.

In this paper we will consider linear inverse problems of the form

$\displaystyle Ax=y,$ (1.1)

where $A:{\mathscr{X}}\to{\mathscr{Y}}$ is a bounded linear operator from a
Banach space ${\mathscr{X}}$ to a Hilbert space ${\mathscr{Y}}$. Throughout
the paper we will assume that $y\in\mbox{Ran}(A)$, the range of $A$, which
guarantees (1.1) has a solution. However, (1.1) may have many solutions. In
order to find a solution with desired feature, by incorporating a priori
available information we take a proper, lower semi-continuous, convex function
${\mathcal{R}}:{\mathscr{X}}\to[0,\infty]$ and search for a solution
$x^{\dagger}$ of (1.1) satisfying

$\displaystyle{\mathcal{R}}(x^{\dagger})=\min\\{{\mathcal{R}}(x):x\in{\mathscr{X}}\mbox{
and }Ax=y\\}$ (1.2)

which is called a ${\mathcal{R}}$-minimizing solution of (1.1).

In practical applications, the exact data $y$ is generally unknown. Very
often, we only have corrupted measurement data at hand. Due to the inherent
ill-posedness property of inverse problems, the ${\mathcal{R}}$-minimizing
solution of (1.1) may not depend continuously on the data. How to use the
noisy data to reconstruct the sought solution therefore becomes a central
topic in computational inverse problems.

To conquer ill-posedness, many iterative regularization methods have been
developed to solve inverse problems using noisy data; see [4, 8, 10, 23, 24,
32] for instance. The performance of these methods depends crucially on a
proper termination of the iterations. In case an accurate upper bound of noise
level is available, various a posteriori rules, including the prominent
discrepancy principle, have been proposed to choose the stopping index of
iteration. In case the information on noise level is difficult to infer, one
has to resort to the so-called heuristic rules ([12, 20, 21, 25, 26]) which
use only the noisy data to select a termination index. Although Bakushiski’s
veto ([1]) states that heuristic rules can not lead to convergence in the
sense of worst case scenario for any regularisation method, such rules often
work as well as, or even better than, the rules using information on noise
level. The drawback of a heuristic rule is that it chooses a stopping index of
iteration based on minimizing a function over ${\mathbb{N}}:=\\{1,2,\cdots\\}$
for which the existence of a minimizer can not be guaranteed unless the noisy
data satisfies certain restrictive conditions.

In this paper we assume that the data can be measured multiple times
repeatedly and the sought solution does not change during the measurement.
This is a common practice in applications where the experiments are set up to
allow acquiring multiple observation data by repeating the experiments.
Therefore, instead of the unknown exact data $y$, we are given multiple
unbiased measurements $y_{1},y_{2},\cdots$ of $y$ which are independent
identically distributed ${\mathscr{Y}}$-valued random variables. Based on
these multiple measurements, for each $n\geq 1$ we use

$\displaystyle\bar{y}^{(n)}:=\frac{1}{n}\sum_{i=1}^{n}y_{i}$ (1.3)

as an estimator of $y$. We then use $\bar{y}^{(n)}$ to construct regularized
solutions and consider their behavior as $n\to\infty$.

Using the average of multiple measurements to decrease the data error is a
standard noise reduction technique in engineering which is called the signal
averaging method ([14, 29, 35]). Assume $(y_{i})$ are ${\mathscr{Y}}$-valued
random variables defined on a probability space
$(\Omega,{\mathcal{F}},{\mathbb{P}})$ with ${\mathbb{E}}$ denoting the
expectation. If $(y_{i})$ are unbiased, independent, identically distributed
estimators of $y$ with a finite variance
$\sigma^{2}:={\mathbb{E}}[\|y_{1}-y\|^{2}]$, we have

$\displaystyle{\mathbb{E}}\left[\|\bar{y}^{(n)}-y\|^{2}\right]=\frac{1}{n^{2}}\sum_{i,j=1}^{n}{\mathbb{E}}\left[\langle
y_{i}-y,y_{j}-y\rangle\right]=\frac{1}{n^{2}}\sum_{i=1}^{n}{\mathbb{E}}\left[\|y_{i}-y\|^{2}\right]=\frac{\sigma^{2}}{n}$

which demonstrates that using $\bar{y}^{(n)}$ as an estimator of $y$ can
reduce the variance by a factor of $n$. Therefore, it is reasonable to expect
more accurate approximate solution can be constructed using the averaged data
$\bar{y}^{(n)}$ than using an individual data $y_{i}$.

In [13] a general class of linear spectral regularization methods have been
considered for solving linear ill-posed problems in Hilbert spaces using the
average of multiple measurement data. Without using any knowledge of the noise
distribution, the convergence and rates of convergence of the regularized
solutions, as the number of measurements increases, have been provided under
either an a priori parameter choice rule or the discrepancy principle. The
results obtained in [13] can be applied mainly for (1.2) with ${\mathscr{X}}$
being a Hilbert space and ${\mathcal{R}}(x)=\|x\|^{2}/2$. In many
applications, however, the sought solution may sit in a Banach space instead
of a Hilbert space, and the sought solution may have a priori known special
features, such as nonnegativity, sparsity and piecewise constancy. Therefore,
it is necessary to develop regularization methods using multiple repeated
measurement data for finding a ${\mathcal{R}}$-minimizing solution of (1.1) in
a general Banach space ${\mathscr{X}}$ with a general convex penalty function
${\mathcal{R}}$.

In this paper we will consider a dual gradient method for solving (1.2) using
multiple repeated measurement data. When the exact data $y$ is available, this
method can be derived by applying the gradient method to the dual problem of
(1.2) and it takes the form (see [22])

$\displaystyle\begin{split}&x_{t}=\arg\min_{x\in{\mathscr{X}}}\left\\{{\mathcal{R}}(x)-\langle\lambda_{t},Ax\rangle\right\\},\\\
&\lambda_{t+1}=\lambda_{t}-\gamma(Ax_{t}-y)\end{split}$ (1.4)

with the initial guess $\lambda_{0}=0$, where
$A^{*}:{\mathscr{Y}}\to{\mathscr{X}}^{*}$ denotes the adjoint of $A$ and
$\gamma>0$ is a step-size. This method actually belongs to the class of Uzawa
methods, see [2, 33, 34]. This method also has connection to the mirror
descent method, see [4, 22, 24, 30]. In case the exact data $y$ is unavailable
and we are given multiple repeated random measurement data
$y_{1},y_{2},\cdots$, we may define $\bar{y}^{(n)}$ by (1.3) for each $n\geq
1$ and replace $y$ in (1.4) by $\bar{y}^{(n)}$ to obtain

$\displaystyle\begin{split}x_{t}^{(n)}&=\arg\min_{x\in{\mathscr{X}}}\left\\{{\mathcal{R}}(x)-\langle\lambda_{t}^{(n)},Ax\rangle\right\\},\\\
\lambda_{t+1}^{(n)}&=\lambda_{t}^{(n)}-\gamma(Ax_{t}^{(n)}-\bar{y}^{(n)})\end{split}$
(1.5)

with the initial guess $\lambda_{0}^{(n)}=0$ and a step size $\gamma>0$. This
is the dual gradient method that we will consider for solving (1.2) using the
average of multiple random measurement data. Since $y_{1},y_{2},\cdots$ are
random variables, we need to consider the convergence of (1.5) based on a
stochastic analysis.

When $X$ is a Hilbert space and ${\mathcal{R}}(x)=\|x\|^{2}/2$, the method
(1.5) becomes the Landweber iteration in Hilbert spaces using the average of
multiple repeated measurement data which has been analyzed in [13] as a
special case of a class of linear spectral regularization methods. The
analysis in [13] is based on the spectral theory and singular value
decomposition of bounded linear self-adjoint operators in Hilbert spaces.
These tools however are no longer applicable to the method (1.5) due to the
possible non-Hilbertian structure of $X$ and non-quadraticity of the
regularization function ${\mathcal{R}}$. In this paper we use tools of convex
analysis in Banach spaces to analyze the method (1.5) with general Banach
space $X$ and general strongly convex function ${\mathcal{R}}$. Without
assuming any knowledge of the noise distribution, we obtain the convergence
and rates of convergence for the method (1.5) as $n\to\infty$ when the method
is terminated by either an a priori stopping rule or a modification of the
discrepancy principle.

The paper is organized as follows, In Section 2 we give a brief review of some
basic facts from convex analysis in Banach spaces. In Section 3 we provide the
convergence analysis of the dual gradient method (1.5) when the method is
terminated by either an a priori stopping rule or a modification of the
discrepancy principle; in particular we derive the convergence rates when the
sought solution satisfies certain variational source conditions. Finally in
Section 4, we provide various numerical results to test the performance of the
proposed method.

##  2\. Preliminaries

In this section, we collect some basic facts on convex analysis in Banach
spaces which will be used in the analysis of the dual gradient method (1.5);
for more details please refer to [36].

Let ${\mathscr{X}}$ be a Banach space, we use ${\mathscr{X}}^{*}$ to denote
its dual space; in case ${\mathscr{X}}$ is a Hilbert space,
${\mathscr{X}}^{*}$ is identified with ${\mathscr{X}}$. Given
$x\in{\mathscr{X}}$ and $\xi\in{\mathscr{X}}^{*}$ we write
$\langle\xi,x\rangle=\xi(x)$ for the duality pairing. For a convex function
$f:{\mathscr{X}}\to(-\infty,\infty]$, its effective domain is

$\mbox{dom}(f):=\\{x\in{\mathscr{X}}:f(x)<\infty\\}.$

If $\mbox{dom}(f)\neq\emptyset$, $f$ is called proper. Given
$x\in\mbox{dom}(f)$, an element $\xi\in{\mathscr{X}}^{*}$ is called a
subgradient of $f$ at $x$ if

$f(\bar{x})\geq
f(x)+\langle\xi,\bar{x}-x\rangle,\quad\forall\bar{x}\in{\mathscr{X}}.$

The collection of all subgradients of $f$ at $x$ is denoted as $\partial f(x)$
and is called the subdifferential of $f$ at $x$. If $\partial
f(x)\neq\emptyset$, then $f$ is called subdifferentiable at $x$. By definition
we have

$\displaystyle x^{*}\in\arg\min_{x\in{\mathscr{X}}}f(x)\Longleftrightarrow
0\in\partial f(x^{*}).$ (2.1)

Note that $x\to\partial f(x)$ defines a set-valued mapping $\partial f$. Let

$\mbox{dom}(\partial f):=\\{x\in\mbox{dom}(f):\partial f(x)\neq\emptyset\\}.$

which is the domain of definition of $\partial f$. Given
$x\in\mbox{dom}(\partial f)$ and $\xi\in\partial f(x)$, the Bregman distance
induced by $f$ at $x$ in the direction $\xi$ is defined by

$D_{f}^{\xi}(\bar{x},x):=f(\bar{x})-f(x)-\langle\xi,\bar{x}-x\rangle,\quad\forall\bar{x}\in{\mathscr{X}}$

which is always nonnegative.

For a proper convex function $f:{\mathscr{X}}\to(-\infty,\infty]$, its convex
conjugate is defined by

$f^{*}(\xi):=\sup_{x\in
X}\\{\langle\xi,x\rangle-f(x)\\},\quad\xi\in{\mathscr{X}}^{*}$

which is a convex function taking values in $(-\infty,\infty]$. By definition
there holds the Fenchel-Young inequality

$\displaystyle f^{*}(\xi)+f(x)\geq\langle\xi,x\rangle,\quad\forall
x\in{\mathscr{X}},\ \xi\in{\mathscr{X}}^{*}.$ (2.2)

If $f:X\to(-\infty,\infty]$ is proper, lower semi-continuous and convex,
$f^{*}$ is also proper and

$\displaystyle\xi\in\partial f(x)\Longleftrightarrow x\in\partial
f^{*}(\xi)\Longleftrightarrow f(x)+f^{*}(\xi)=\langle\xi,x\rangle.$ (2.3)

In applications, optimization problems of the form

$\inf_{x\in{\mathscr{X}}}\left\\{f(x)+g(Ax)\right\\}$

are frequently encountered, where $A:{\mathscr{X}}\to{\mathscr{Y}}$ is a
bounded linear operator between two Banach spaces ${\mathscr{X}}$ and
${\mathscr{Y}}$, and $f:{\mathscr{X}}\to(-\infty,\infty]$ and
$g:{\mathscr{Y}}\to(-\infty,\infty]$ are proper convex functions. Let
$A^{*}:{\mathscr{Y}}^{*}\to{\mathscr{X}}^{*}$ denote the adjoint of $A$. By
the Fenchel-Young inequality we have

$f(x)\geq\langle\xi,x\rangle-f^{*}(\xi),\qquad
g(Ax)\geq\langle-\eta,Ax\rangle-g^{*}(-\eta)$

and therefore

$f(x)+g(Ax)\geq\langle\xi-A^{*}\eta,x\rangle-f^{*}(\xi)-g^{*}(-\eta)$

for all $x\in{\mathscr{X}}$, $\xi\in{\mathscr{X}}^{*}$ and
$\eta\in{\mathscr{Y}}^{*}$. Thus, taking $\xi=A^{*}\eta$ gives

$\inf_{x\in{\mathscr{X}}}\left\\{f(x)+g(Ax)\right\\}\geq\sup_{\eta\in{\mathscr{Y}}^{*}}\left\\{-f^{*}(A^{*}\eta)-g^{*}(-\eta)\right\\}.$

It is natural to ask under what conditions on $f$ and $g$ we have equality in
the above equation. This is answered by the Fenchel-Rockafellar duality
formula given below (see [36, Corollary 2.8.5]).

######  Proposition 2.1.

Let ${\mathscr{X}}$ and ${\mathscr{Y}}$ be Banach spaces, let
$f:{\mathscr{X}}\to(-\infty,\infty]$ and $g:{\mathscr{Y}}\to(-\infty,\infty]$
be proper, convex functions, and let $A:{\mathscr{X}}\to{\mathscr{Y}}$ be a
bounded linear operator. If there is $x_{0}\in\emph{dom}(f)$ such that
$Ax_{0}\in\emph{dom}(g)$ and $g$ is continuous at $Ax_{0}$, then

$\displaystyle\inf_{x\in{\mathscr{X}}}\\{f(x)+g(Ax)\\}=\sup_{\eta\in{\mathscr{Y}}^{*}}\\{-f^{*}(A^{*}\eta)-g^{*}(-\eta)\\}.$
(2.4)

A proper function $f:{\mathscr{X}}\to(-\infty,\infty]$ is called
$c_{0}$-strongly convex for some constant $c_{0}>0$ if

$\displaystyle f(s\bar{x}+(1-s)x)+c_{0}s(1-s)\|\bar{x}-x\|^{2}\leq
sf(\bar{x})+(1-s)f(x)$ (2.5)

for all $\bar{x},x\in\mbox{dom}(f)$ and $s\in[0,1]$. It is easy to show that
for a proper $c_{0}$-strongly convex function
$f:{\mathscr{X}}\to(-\infty,\infty]$ there holds

$\displaystyle D_{f}^{\xi}(\bar{x},x)\geq c_{0}\|x-\bar{x}\|^{2}$ (2.6)

for all $\bar{x}\in\mbox{dom}(f)$, $x\in\mbox{dom}(\partial f)$ and
$\xi\in\partial f(x)$. Furthermore, [36, Corollary 3.5.11] contains the
following important result concerning the differentiability of $f^{*}$ for
strongly convex $f$.

######  Proposition 2.2.

Let ${\mathscr{X}}$ be a Banach space and let
$f:{\mathscr{X}}\to(-\infty,\infty]$ be a proper, lower semi-continuous,
$c_{0}$-strongly convex function for some constant $c_{0}>0$. Then
$\emph{dom}(f^{*})={\mathscr{X}}^{*}$, $f^{*}$ is Fréchet differentiable and
its gradient $\nabla f^{*}$ maps ${\mathscr{X}}^{*}$ into ${\mathscr{X}}$ with

$\|\nabla f^{*}(\xi)-\nabla f^{*}(\eta)\|\leq\frac{\|\xi-\eta\|}{2c_{0}}$

for all $\xi,\eta\in{\mathscr{X}}^{*}$.

It should be emphasized that ${\mathscr{X}}$ in Proposition 2.2 can be an
arbitrary Banach space; it can even be a general normed vector spaces. The
gradient $\nabla f^{*}$ of $f^{*}$ is in general a mapping from
${\mathscr{X}}^{*}\to{\mathscr{X}}^{**}$, the second dual space of
${\mathscr{X}}$. Proposition 2.2 actually concludes that, for each
$\xi\in{\mathscr{X}}^{*}$, $\nabla f^{*}(\xi)$ is an element in
${\mathscr{X}}$, and thus $\nabla f^{*}$ is a mapping from ${\mathscr{X}}^{*}$
to ${\mathscr{X}}$.

Given a proper strongly convex function $f:{\mathscr{X}}\to(-\infty,\infty]$,
we may consider for each $\xi\in{\mathscr{X}}^{*}$ the convex minimization
problem

$\displaystyle\min_{x\in{\mathscr{X}}}\left\\{f(x)-\langle\xi,x\rangle\right\\}$
(2.7)

which is involved in the dual gradient method (1.5). According to [36, Theorem
3.5.8], (2.1), (2.3) and Proposition 2.2 we have

######  Proposition 2.3.

If $f:{\mathscr{X}}\to(-\infty,\infty]$ is a proper, lower semi-continuous,
strongly convex function, then for any $\xi\in{\mathscr{X}}^{*}$ the
minimization problem (2.7) has a unique minimizer given by $\nabla
f^{*}(\xi)$.

##  3\. Convergence analysis

In this section we will provide convergence analysis of the method (1.5) when
it is terminated by either an a priori stopping rule or a statistical variant
of the discrepancy principle. We will make use of the following assumption.

######  Assumption 3.1.

  1. _(i)_

${\mathscr{X}}$ is a Banach space, ${\mathscr{Y}}$ is a Hilbert space, and
$A:{\mathscr{X}}\to{\mathscr{Y}}$ is a bounded linear operator.

  2. _(ii)_

${\mathcal{R}}:{\mathscr{X}}\to[0,\infty]$ is proper, lower semi-continuous
and $c_{0}$-strongly convex for some constant $c_{0}>0$.

  3. _(iii)_

The equation $Ax=y$ has a solution in $\emph{dom}({\mathcal{R}})$.

  4. _(iv)_

$y_{1},y_{2},\cdots$ is a sequence of independent identically distributed
${\mathscr{Y}}$-valued random variables defined on a probability space
$(\Omega,{\mathcal{F}},{\mathbb{P}})$ with ${\mathbb{E}}[y_{1}]=y$ and
$0<\sigma^{2}:={\mathbb{E}}[\|y_{1}-y\|^{2}]<\infty$.

For a sequence of noisy data $\\{y_{i}\\}$ satisfying (iv) in Assumption 3.1,
we define $\bar{y}^{(n)}$ by (1.3) for each $n\geq 1$. There holds

$\displaystyle{\mathbb{E}}\left[\|\bar{y}^{(n)}-y\|^{2}\right]=\frac{\sigma^{2}}{n}.$
(3.1)

According to (i)-(iii) in Assumption 3.1 and Proposition 2.3, the equation
(1.1) has a unique ${\mathcal{R}}$-minimizing solution, denoted as
$x^{\dagger}$, and the method (1.5) is well-defined with

$\displaystyle
A^{*}\lambda_{t}^{(n)}\in\partial{\mathcal{R}}(x_{t}^{(n)})\quad\mbox{ and
}\quad x_{t}^{(n)}=\nabla{\mathcal{R}}^{*}(A^{*}\lambda_{t}^{(n)})$ (3.2)

for each integer $t\geq 0$.

Our convergence analysis of the method (1.5), in particular the derivation of
convergence rates, is based on the following result which has been established
in [22].

######  Proposition 3.1.

Let Assumption 3.1 hold and let
$d_{\bar{y}^{(n)}}(\lambda):={\mathcal{R}}^{*}(A^{*}\lambda)-\langle\lambda,\bar{y}^{(n)}\rangle$.
Let $L:=\|A\|^{2}/(2c_{0})$. Consider the dual gradient method (1.5) with
$\lambda_{0}^{(n)}=0$. If $0<\gamma<1/L$ then for any
$\lambda\in{\mathscr{Y}}$ there holds

$\displaystyle
d_{\bar{y}^{(n)}}(\lambda)-d_{\bar{y}^{(n)}}(\lambda_{t+1}^{(n)})$
$\displaystyle\geq\frac{1}{2\gamma(t+1)}\left(\|\lambda-\lambda_{t+1}^{(n)}\|^{2}-\|\lambda\|^{2}\right)$
$\displaystyle\quad\,+\left\\{\left(\frac{1}{2}-\frac{L\gamma}{4}\right)t+\left(\frac{1}{2}-\frac{L\gamma}{2}\right)\right\\}\gamma\|Ax_{t}^{(n)}-\bar{y}^{(n)}\|^{2}$

for all integers $t\geq 0$.

Based on Proposition 3.1, we can obtain the following estimates which will be
used frequently in the forthcoming analysis.

######  Lemma 3.2.

Let Assumption 3.1 hold and consider the method (1.5) with
$\lambda_{0}^{(n)}=0$. Let $L:=\|A\|^{2}/(2c_{0})$ and assume $0<\gamma<1/L$.
Then

$c\gamma(t+1)\|Ax_{t}^{(n)}-\bar{y}^{(n)}\|^{2}+\frac{1}{8\gamma(t+1)}\|\lambda_{t+1}^{(n)}\|^{2}\leq\eta_{t}+\gamma(t+1)\|\bar{y}^{(n)}-y\|^{2}$

for all $t\geq 0$, where $c:=1/2-L\gamma/2>0$ and

$\eta_{t}:=\sup_{x\in{\mathscr{X}}}\left\\{{\mathcal{R}}(x^{\dagger})-{\mathcal{R}}(x)-\frac{1}{3}\gamma(t+1)\|Ax-y\|^{2}\right\\}.$

###### Proof.

From Proposition 3.1 it follows for all $\lambda\in{\mathscr{Y}}$ that

$\displaystyle c\gamma(t+1)\|Ax_{t}^{(n)}-\bar{y}^{(n)}\|^{2}$
$\displaystyle\leq
d_{\bar{y}^{(n)}}(\lambda)-d_{\bar{y}^{(n)}}(\lambda_{t+1}^{(n)})-\frac{1}{2\gamma(t+1)}\left(\|\lambda_{t+1}^{(n)}-\lambda\|^{2}-\|\lambda\|^{2}\right).$

Let $d_{y}(\lambda):={\mathcal{R}}^{*}(A^{*}\lambda)-\langle\lambda,y\rangle$.
We then have

$\displaystyle c\gamma(t+1)\|Ax_{t}^{(n)}-\bar{y}^{(n)}\|^{2}$
$\displaystyle\leq
d_{y}(\lambda)-d_{y}(\lambda_{t+1}^{(n)})+\langle\lambda_{t+1}^{(n)}-\lambda,\bar{y}^{(n)}-y\rangle$
$\displaystyle\quad\,-\frac{1}{2\gamma(t+1)}\left(\|\lambda_{t+1}^{(n)}-\lambda\|^{2}-\|\lambda\|^{2}\right).$
(3.3)

By the Cauchy-Schwarz inequality we have

$\displaystyle\langle\lambda_{t+1}^{(n)}-\lambda,\bar{y}^{(n)}-y\rangle$
$\displaystyle\leq\|\lambda_{t+1}^{(n)}-\lambda\|\|\bar{y}^{(n)}-y\|\leq\frac{1}{4\gamma(t+1)}\|\lambda_{t+1}^{(n)}-\lambda\|^{2}+\gamma(t+1)\|\bar{y}^{(n)}-y\|^{2}.$

Combining this with (3) gives

$\displaystyle
c\gamma(t+1)\|Ax_{t}^{(n)}-\bar{y}^{(n)}\|^{2}+\frac{1}{4\gamma(t+1)}\|\lambda_{t+1}^{(n)}-\lambda\|^{2}$
$\displaystyle\leq
d_{y}(\lambda)-d_{y}(\lambda_{t+1}^{(n)})+\frac{1}{2\gamma(t+1)}\|\lambda\|^{2}+\gamma(t+1)\|\bar{y}^{(n)}-y\|^{2}$

which together with the inequality $\|\lambda_{t+1}^{(n)}\|^{2}\leq
2(\|\lambda_{t+1}^{(n)}-\lambda\|^{2}+\|\lambda\|^{2})$ then shows that

$\displaystyle
c\gamma(t+1)\|Ax_{t}^{(n)}-\bar{y}^{(n)}\|^{2}+\frac{1}{8\gamma(t+1)}\|\lambda_{t+1}^{(n)}\|^{2}$
$\displaystyle\leq
d_{y}(\lambda)-d_{y}(\lambda_{t+1}^{(n)})+\frac{3\|\lambda\|^{2}}{4\gamma(t+1)}+\gamma(t+1)\|\bar{y}^{(n)}-y\|^{2}.$
(3.4)

By using the fact $y=Ax^{\dagger}$ and the Fenchel-Young inequality (2.2), we
have

$d_{y}(\lambda_{t+1}^{(n)})={\mathcal{R}}^{*}(A^{*}\lambda_{t+1}^{(n)})-\langle
A^{*}\lambda_{t+1}^{(n)},x^{\dagger}\rangle\geq-{\mathcal{R}}(x^{\dagger}).$

Therefore

$\displaystyle
c\gamma(t+1)\|Ax_{t}^{(n)}-\bar{y}^{(n)}\|^{2}+\frac{1}{8\gamma(t+1)}\|\lambda_{t+1}^{(n)}\|^{2}$
$\displaystyle\leq{\mathcal{R}}^{*}(A^{*}\lambda)-\langle\lambda,y\rangle+{\mathcal{R}}(x^{\dagger})+\frac{3\|\lambda\|^{2}}{4\gamma(t+1)}+\gamma(t+1)\|\bar{y}^{(n)}-y\|^{2}$

for all $\lambda\in{\mathscr{Y}}$ and hence

$\displaystyle
c\gamma(t+1)\|Ax_{t}^{(n)}-\bar{y}^{(n)}\|^{2}+\frac{1}{8\gamma(t+1)}\|\lambda_{t+1}^{(n)}\|^{2}$
$\displaystyle\leq\inf_{\lambda\in{\mathscr{Y}}}\left\\{{\mathcal{R}}^{*}(A^{*}\lambda)-\langle\lambda,y\rangle+{\mathcal{R}}(x^{\dagger})+\frac{3\|\lambda\|^{2}}{4\gamma(t+1)}\right\\}+\gamma(t+1)\|\bar{y}^{(n)}-y\|^{2}$
$\displaystyle={\mathcal{R}}(x^{\dagger})-\sup_{\lambda\in{\mathscr{Y}}}\left\\{-{\mathcal{R}}^{*}(A^{*}\lambda)+\langle\lambda,y\rangle-\frac{3\|\lambda\|^{2}}{4\gamma(t+1)}\right\\}+\gamma(t+1)\|\bar{y}^{(n)}-y\|^{2}$

By the Fenchel-Rockafellar duality formula given in Proposition 2.1 we finally
obtain

$\displaystyle
c\gamma(t+1)\|Ax_{t}^{(n)}-\bar{y}^{(n)}\|^{2}+\frac{1}{8\gamma(t+1)}\|\lambda_{t+1}^{(n)}\|^{2}$
$\displaystyle\leq{\mathcal{R}}(x^{\dagger})-\inf_{x\in{\mathscr{X}}}\left\\{{\mathcal{R}}(x)+\frac{1}{3}\gamma(t+1)\|Ax-y\|^{2}\right\\}+\gamma(t+1)\|\bar{y}^{(n)}-y\|^{2}$
$\displaystyle=\eta_{t}+\gamma(t+1)\|\bar{y}^{(n)}-y\|^{2}.$

The proof is therefore complete. ∎

###  3.1. Convergence analysis under a priori stopping rule

In this subsection we consider the dual gradient method (1.5) and show that
$x_{t_{n}}^{(n)}$ converges to the unique ${\mathcal{R}}$-minimizing solution
$x^{\dagger}$ of (1.1) if $t_{n}$ is chosen such that $t_{n}\to\infty$ and
$t_{n}/n\to 0$ as $n\to\infty$. Furthermore we derive the convergence rates
under suitable a priori choice of $t_{n}$ when $x^{\dagger}$ satisfies certain
variational source conditions.

By using [24, Lemma 3.7] we immediately obtain the following convergence
result for the sequence $\\{x_{t},\lambda_{t}\\}$ defined by the method (1.4)
with exact data.

######  Lemma 3.3.

Let (i)-(iii) in Assumption 3.1 hold. Consider the method (1.4) with exact
data. If $\lambda_{0}=0$ and $0<\gamma<4c_{0}/\|A\|^{2}$, then
$D_{\mathcal{R}}^{A^{*}\lambda_{t}}(x^{\dagger},x_{t})\to 0$ as $t\to\infty$.

Although $\\{x_{t},\lambda_{t}\\}$ are deterministic,
$\\{x_{t}^{(n)},\lambda_{t}^{(n)}\\}$ are random variables. Therefore, we need
to consider $\\{x_{t}^{(m)},\lambda_{t}^{(n)}\\}$ by a stochastic analysis.

######  Lemma 3.4.

Let Assumption 3.1 hold. Let $\\{x_{t},\lambda_{t}\\}$ and
$\\{x_{t}^{(n)},\lambda_{t}^{(n)}\\}$ be defined by (1.4) and (1.5)
respectively with the initial guess $\lambda_{0}=\lambda_{0}^{(n)}=0$. Then
for all integers $t\geq 0$ there hold

${\mathbb{E}}\left[\|x_{t}^{(n)}-x_{t}\|^{2}\right]\to 0\quad\mbox{ and
}\quad{\mathbb{E}}\left[\|\lambda_{t}^{(n)}-\lambda_{t}\|^{2}\right]\to 0$

as $n\to\infty$.

###### Proof.

We will use an induction argument. Since $\lambda_{0}^{(n)}=\lambda_{0}=0$, we
have
$x_{0}^{(n)}=\nabla{\mathcal{R}}^{*}(A^{*}\lambda_{0}^{(n)})=\nabla{\mathcal{R}}^{*}(A^{*}\lambda_{0})=x_{0}$.
Thus the result is true for $t=0$. Now we assume that the result is true for
some $t\geq 0$. Then, by the definition of $\lambda_{t+1}^{(n)}$ and
$\lambda_{t+1}$, we have

$\|\lambda_{t+1}^{(n)}-\lambda_{t+1}\|^{2}=\|(\lambda_{t}^{(n)}-\lambda_{t})-\gamma
A(x_{t}^{(n)}-x_{t})+\gamma(\bar{y}^{(n)}-y)\|^{2}.$

By using the inequality $(a+b+c)^{2}\leq 3(a^{2}+b^{2}+c^{2})$, taking the
expectation, and using (3.1), we can obtain

$\displaystyle{\mathbb{E}}\left[\|\lambda_{t+1}^{(n)}-\lambda_{t+1}\|^{2}\right]$
$\displaystyle\leq
3{\mathbb{E}}\left[\|\lambda_{t}^{(n)}-\lambda_{t}\|^{2}\right]+3\gamma^{2}{\mathbb{E}}\left[\|A(x_{t}^{(n)}-x_{t})\|^{2}\right]+3\gamma^{2}{\mathbb{E}}\left[\|\bar{y}^{(n)}-y\|^{2}\right]$
$\displaystyle\leq
3{\mathbb{E}}\left[\|\lambda_{t}^{(n)}-\lambda_{t}\|^{2}\right]+3\gamma^{2}\|A\|^{2}{\mathbb{E}}\left[\|x_{t}^{(n)}-x_{t}\|^{2}\right]+\frac{3\gamma^{2}\sigma^{2}}{n}.$

Thus, by the induction hypothesis we can conclude that
${\mathbb{E}}[\|\lambda_{t+1}^{(n)}-\lambda_{t+1}\|^{2}]\to 0$ as
$n\to\infty$. Note that

$x_{t+1}^{(n)}=\nabla{\mathcal{R}}^{*}(A^{*}\lambda_{t+1}^{(n)})\quad\mbox{
and }\quad x_{t+1}=\nabla{\mathcal{R}}^{*}(A^{*}\lambda_{t+1}).$

By Proposition 2.2 we have

$\displaystyle{\mathbb{E}}\left[\|x_{t+1}^{(n)}-x_{t+1}\|^{2}\right]$
$\displaystyle={\mathbb{E}}\left[\|\nabla{\mathcal{R}}^{*}(A^{*}\lambda_{t+1}^{(n)})-\nabla{\mathcal{R}}^{*}(A^{*}\lambda_{t+1})\|^{2}\right]$
$\displaystyle\leq\frac{1}{4c_{0}^{2}}{\mathbb{E}}\left[\|A^{*}(\lambda_{t+1}^{(n)}-\lambda_{t+1})\|^{2}\right]$
$\displaystyle\leq\frac{\|A\|^{2}}{4c_{0}^{2}}{\mathbb{E}}\left[\|\lambda_{t+1}^{(n)}-\lambda_{t+1}\|^{2}\right]$
$\displaystyle\to 0$

as $n\to\infty$. The proof is complete. ∎

######  Lemma 3.5.

Let Assumption 3.1 hold. Consider the method (1.5) with $\lambda_{0}^{(n)}=0$
and assume that $0<\gamma<4c_{0}/\|A\|^{2}$. Let
$\Delta_{t}^{(n)}:=D_{\mathcal{R}}^{A^{*}\lambda_{t}^{(n)}}(x^{\dagger},x_{t}^{(n)})$.
Then

${\mathbb{E}}\left[\Delta_{t+1}^{(n)}\right]-{\mathbb{E}}\left[\Delta_{t}^{(n)}\right]\leq
C_{1}\frac{\sigma^{2}}{n}$

for all integers $t\geq 0$, where
$C_{1}:=c_{0}\gamma/(4c_{0}-\gamma\|A\|^{2})$.

###### Proof.

Note that

$\displaystyle\Delta_{t+1}^{(n)}-\Delta_{t}^{(n)}=\left(\langle
A^{*}\lambda_{t+1}^{(n)},x_{t+1}^{(n)}-x^{\dagger}\rangle-{\mathcal{R}}(x_{t+1}^{(n)})\right)+\left({\mathcal{R}}(x_{t}^{(n)})+\langle
A^{*}\lambda_{t}^{(n)},x^{\dagger}-x_{t}^{(n)}\rangle\right).$

By using the fact
$A^{*}\lambda_{t}^{(n)}\in\partial{\mathcal{R}}(x_{t}^{(n)})$ we have from
(2.3) that

${\mathcal{R}}(x_{t}^{(n)})+{\mathcal{R}}^{*}(A^{*}\lambda_{t}^{(n)})=\langle
A^{*}\lambda_{t}^{(n)},x_{t}^{(n)}\rangle$

for all $t\geq 0$. Therefore

$\displaystyle\Delta_{t+1}^{(n)}-\Delta_{t}^{(n)}=\left({\mathcal{R}}^{*}(A^{*}\lambda_{t+1}^{(n)})-\langle
A^{*}\lambda_{t+1}^{(n)},x^{\dagger}\rangle\right)-\left({\mathcal{R}}^{*}(A^{*}\lambda_{t}^{(n)})-\langle
A^{*}\lambda_{t}^{(n)},x^{\dagger}\rangle\right).$

Since $x_{t}^{(n)}=\nabla{\mathcal{R}}^{*}(A^{*}\lambda_{t}^{(n)})$, by using
Proposition 2.2 we can obtain

$\displaystyle\Delta_{t+1}^{(n)}-\Delta_{t}^{(n)}$
$\displaystyle=\left({\mathcal{R}}^{*}(A^{*}\lambda_{t+1}^{(n)})-{\mathcal{R}}^{*}(A^{*}\lambda_{t}^{(n)})-\langle
A^{*}\lambda_{t+1}^{(n)}-A^{*}\lambda_{t}^{(n)},\nabla{\mathcal{R}}^{*}(A^{*}\lambda_{t}^{(n)})\rangle\right)$
$\displaystyle\quad\,+\langle
A^{*}\lambda_{t+1}^{(n)}-A^{*}\lambda_{t}^{(n)},x_{t}^{(n)}-x^{\dagger}\rangle$
$\displaystyle\leq\frac{1}{4c_{0}}\|A^{*}\lambda_{t+1}^{(n)}-A^{*}\lambda_{t}^{(n)}\|^{2}+\langle
A^{*}\lambda_{t+1}^{(n)}-A^{*}\lambda_{t}^{(n)},x_{t}^{(n)}-x^{\dagger}\rangle.$

According to the definition of $\lambda_{t+1}^{(n)}$ it is easy to see that

$A^{*}\lambda_{t+1}^{(n)}-A^{*}\lambda_{t}^{(n)}=-\gamma
A^{*}(Ax_{t}^{(n)}-\bar{y}^{(n)}).$

Therefore

$\displaystyle\Delta_{t+1}^{(n)}-\Delta_{t}^{(n)}$
$\displaystyle\leq\frac{\gamma^{2}}{4c_{0}}\|A^{*}(Ax_{t}^{(n)}-\bar{y}^{(n)})\|^{2}-\gamma\langle
Ax_{t}^{(n)}-\bar{y}^{(n)},Ax_{t}^{(n)}-y\rangle$
$\displaystyle\leq-\left(1-\frac{\gamma\|A\|^{2}}{4c_{0}}\right)\gamma\|Ax_{t}^{(n)}-\bar{y}^{(n)}\|^{2}+\gamma\|\bar{y}^{(n)}-y\|\|Ax_{t}^{(n)}-\bar{y}^{(n)}\|.$

By using the Young’s inequality we have

$\displaystyle\gamma\|\bar{y}^{(n)}-y\|\|Ax_{t}^{(n)}-\bar{y}^{(n)}\|\leq\left(1-\frac{\gamma\|A\|^{2}}{4c_{0}}\right)\gamma\|Ax_{t}^{(n)}-\bar{y}^{(n)}\|^{2}+C_{1}\|\bar{y}^{(n)}-y\|^{2}.$

with $C_{1}$ defined above. Therefore

$\displaystyle\Delta_{t+1}^{(n)}-\Delta_{t}^{(n)}\leq
C_{1}\|\bar{y}^{(n)}-y\|^{2}.$

By taking the expectation and using (3.1) we then obtain

$\displaystyle{\mathbb{E}}\left[\Delta_{t+1}^{(n)}\right]-{\mathbb{E}}\left[\Delta_{t}^{(n)}\right]\leq
C_{1}{\mathbb{E}}\left[\|\bar{y}^{(n)}-y\|^{2}\right]=C_{1}\frac{\sigma^{2}}{n}$

which shows the desired inequality. ∎

######  Theorem 3.6.

Let Assumption 3.1 hold. Consider the method (1.5) with $\lambda_{0}^{(n)}=0$
and assume that $0<\gamma<4c_{0}/\|A\|^{2}$. If $t_{n}$ is chosen such that
$t_{n}\to\infty$ and $t_{n}/n\to 0$ as $n\to\infty$, then

${\mathbb{E}}\left[D_{\mathcal{R}}^{A^{*}\lambda_{t_{n}}^{(n)}}(x^{\dagger},x_{t_{n}}^{(n)})\right]\to
0$

and hence ${\mathbb{E}}\left[\|x_{t_{n}}^{(n)}-x^{\dagger}\|^{2}\right]\to 0$
as $n\to\infty$.

###### Proof.

By the strong convexity of ${\mathcal{R}}$ it suffices to show
${\mathbb{E}}[\Delta_{t_{n}}^{(n)}]\to 0$ as $n\to\infty$. Let $t\geq 0$ be
any fixed integer. Since $t_{n}\to\infty$ as $n\to\infty$, we have $t_{n}>t$
for large $n$. Thus, we may repeatedly use Lemma 3.5 to obtain

${\mathbb{E}}\left[\Delta_{t_{n}}^{(n)}\right]\leq{\mathbb{E}}\left[\Delta_{t}^{(n)}\right]+C_{1}(t_{n}-t)\frac{\sigma^{2}}{n}.$

Since $t_{n}/n\to 0$ as $n\to\infty$, we thus have

$\displaystyle\limsup_{n\to\infty}{\mathbb{E}}\left[\Delta_{t_{n}}^{(n)}\right]\leq\limsup_{n\to\infty}{\mathbb{E}}\left[\Delta_{t}^{(n)}\right]$
(3.5)

for any $t\geq 0$. Note that

$\displaystyle{\mathbb{E}}\left[\Delta_{t}^{(n)}\right]-D_{{\mathcal{R}}}^{A^{*}\lambda_{t}}(x^{\dagger},x_{t})$
$\displaystyle=\left(\langle
A^{*}\lambda_{t},x^{\dagger}-x_{t}\rangle-{\mathbb{E}}\left[\langle
A^{*}\lambda_{t}^{(n)},x^{\dagger}-x_{t}^{(n)}\rangle\right]\right)$
$\displaystyle\quad\,+\left({\mathcal{R}}(x_{t})-{\mathbb{E}}\left[{\mathcal{R}}(x_{t}^{(n)})\right]\right).$
(3.6)

With the help of the Cauchy-Schwarz inequality we have

$\displaystyle\left|{\mathbb{E}}\left[\langle
A^{*}\lambda_{t}^{(n)},x^{\dagger}-x_{t}^{(n)}\rangle\right]-\langle
A^{*}\lambda_{t},x^{\dagger}-x_{t}\rangle\right|$
$\displaystyle\leq\left|{\mathbb{E}}\left[\langle
A^{*}(\lambda_{t}^{(n)}-\lambda_{t}),x^{\dagger}-x_{t}^{(n)}\rangle\right]\right|+\left|{\mathbb{E}}\left[\langle
A^{*}\lambda_{t},x_{t}-x_{t}^{(n)}\rangle\right]\right|$
$\displaystyle\leq\|A\|\left({\mathbb{E}}\left[\|\lambda_{t}^{(n)}-\lambda_{t}\|^{2}\right]\right)^{1/2}\left({\mathbb{E}}\left[\|x^{\dagger}-x_{t}^{(n)}\|^{2}\right]\right)^{1/2}$
$\displaystyle\quad\,+\|A^{*}\lambda_{t}\|\left({\mathbb{E}}\left[\|x_{t}-x_{t}^{(n)}\|^{2}\right]\right)^{1/2}.$

Thus we may use Lemma 3.4 to conclude

$\displaystyle\lim_{n\to\infty}{\mathbb{E}}\left[\langle
A^{*}\lambda_{t}^{(n)},x^{\dagger}-x_{t}^{(n)}\rangle\right]=\langle
A^{*}\lambda_{t},x^{\dagger}-x_{t}\rangle.$ (3.7)

Next we will show that

$\displaystyle{\mathcal{R}}(x_{t})\leq\liminf_{n\to\infty}{\mathbb{E}}\left[{\mathcal{R}}(x_{t}^{(n)})\right].$
(3.8)

To see this, we take a subsequence $\\{n_{k}\\}$ with $n_{k}\to\infty$ as
$k\to\infty$ such that

$\lim_{k\to\infty}{\mathbb{E}}\left[{\mathcal{R}}(x_{t}^{(n_{k})})\right]=\liminf_{n\to\infty}{\mathbb{E}}\left[{\mathcal{R}}(x_{t}^{(n)})\right].$

According to Lemma 3.4 we have ${\mathbb{E}}[\|x_{t}^{(n_{k})}-x_{t}\|^{2}]\to
0$ as $k\to\infty$. By taking a subsequence of $\\{n_{k}\\}$ if necessary, we
can guarantee $\|x_{t}^{(n_{k})}-x_{t}\|\to 0$ as $k\to\infty$ almost surely.
Thus, from the lower semi-continuity of ${\mathcal{R}}$ and Fatou’s lemma it
follows

$\displaystyle{\mathcal{R}}(x_{t})$
$\displaystyle={\mathbb{E}}[{\mathcal{R}}(x_{t})]\leq{\mathbb{E}}\left[\liminf_{k\to\infty}{\mathcal{R}}(x_{t}^{(n_{k})})\right]\leq\liminf_{k\to\infty}{\mathbb{E}}\left[{\mathcal{R}}(x_{t}^{(n_{k})})\right]$
$\displaystyle=\lim_{k\to\infty}{\mathbb{E}}\left[{\mathcal{R}}(x_{t}^{(n_{k})})\right]=\liminf_{n\to\infty}{\mathbb{E}}\left[{\mathcal{R}}(x_{t}^{(n)})\right]$

which shows (3.8). Consequently it follows from (3.1), (3.7) and (3.8) that

$\limsup_{n\to\infty}{\mathbb{E}}\left[\Delta_{t}^{(n)}\right]\leq
D_{{\mathcal{R}}}^{A^{*}\lambda_{t}}(x^{\dagger},x_{t}).$

Therefore we may use (3.5) to obtain

$\displaystyle\limsup_{n\to\infty}{\mathbb{E}}\left[\Delta_{t_{n}}^{(n)}\right]\leq
D_{{\mathcal{R}}}^{A^{*}\lambda_{t}}(x^{\dagger},x_{t})$

for all $t\geq 0$. Letting $t\to\infty$ and using Lemma 3.3 we thus obtain
${\mathbb{E}}[\Delta_{t_{n}}^{(n)}]\to 0$ as $n\to\infty$. ∎

We next consider deriving convergence rates of the method (1.5) under an a
priori stopping rule when the sought solution $x^{\dagger}$ satisfies the
variational source conditions specified in the following assumption.

######  Assumption 3.2.

For the unique solution $x^{\dagger}$ of (1.1) there is an error measure
function ${\mathcal{E}}^{\dagger}:\emph{dom}({\mathcal{R}})\to[0,\infty)$ with
${\mathcal{E}}^{\dagger}(x^{\dagger})=0$ such that

${\mathcal{E}}^{\dagger}(x)\leq{\mathcal{R}}(x)-{\mathcal{R}}(x^{\dagger})+M\|Ax-y\|^{q},\quad\forall
x\in\emph{dom}({\mathcal{R}})$

for some $0<q\leq 1$ and some constant $M>0$.

######  Remark 3.1.

Variational source conditions were first introduced in [15], as a
generalization of the spectral source conditions in Hilbert spaces, to derive
convergence rates of Tikhonov regularization in Banach spaces. This kind of
source conditions was further generalized and refined subsequently, see [9,
16, 17, 19] for instance. The error measure function ${\mathcal{E}}^{\dagger}$
in Assumption 3.2 is used to measure the speed of convergence; the usual
choice of ${\mathcal{E}}^{\dagger}$ is the Bregman distance induced by
${\mathcal{R}}$.

######  Remark 3.2.

When $x^{\dagger}$ satisfies the benchmark source condition
$A^{*}\lambda^{\dagger}\in\partial{\mathcal{R}}(x^{\dagger})$ for some
$\lambda^{\dagger}\in{\mathscr{Y}}$, it is straightforward to see that

$D_{\mathcal{R}}^{A^{*}\lambda^{\dagger}}(x,x^{\dagger})\leq{\mathcal{R}}(x)-{\mathcal{R}}(x^{\dagger})+\|\lambda^{\dagger}\|\|Ax-y\|,\quad\forall
x\in\mbox{dom}({\mathcal{R}})$

which shows the variational source condition is satisfied with
${\mathcal{E}}^{\dagger}(x)=D_{\mathcal{R}}^{A^{*}\lambda^{\dagger}}(x,x^{\dagger})$,
$M=\|\lambda^{\dagger}\|$ and $q=1$; this is a well-known fact, see [15].

######  Remark 3.3.

When both $X$ and $Y$ are Hilbert spaces, ${\mathcal{R}}$ is $c_{0}$-strongly
convex for some constant $c_{0}>0$, and $x^{\dagger}$ satisfies the source
condition

$\displaystyle\xi^{\dagger}:=(A^{*}A)^{\nu/2}\omega\in\partial{\mathcal{R}}(x^{\dagger})$
(3.9)

for some $0<\nu\leq 1$ and $\omega\in X$, then the variational source
condition in Assumption 3.2 holds with

${\mathcal{E}}^{\dagger}(x)=\frac{1}{2}D_{\mathcal{R}}^{\xi^{\dagger}}(x,x^{\dagger}),\quad
M=C_{\nu}\|\omega\|^{\frac{2}{1+\nu}}\quad\mbox{and}\quad
q=\frac{2\nu}{1+\nu},$

where
$C_{\nu}:=\frac{1+\nu}{2}\left(\frac{1-\nu}{c_{0}}\right)^{\frac{1-\nu}{1+\nu}}$.
Indeed, by the given condition (3.9), we can obtain

$\displaystyle D_{\mathcal{R}}^{\xi^{\dagger}}(x,x^{\dagger})$
$\displaystyle={\mathcal{R}}(x)-{\mathcal{R}}(x^{\dagger})-\langle\omega,(A^{*}A)^{\nu/2}(x-x^{\dagger})\rangle$
$\displaystyle\leq{\mathcal{R}}(x)-{\mathcal{R}}(x^{\dagger})+\|\omega\|\|(A^{*}A)^{\nu/2}(x-x^{\dagger})\|.$

By the interpolation inequality ([8]) we then have

$D_{\mathcal{R}}^{\xi^{\dagger}}(x,x^{\dagger})\leq{\mathcal{R}}(x)-{\mathcal{R}}(x^{\dagger})+\|\omega\|\|x-x^{\dagger}\|^{1-\nu}\|A(x-x^{\dagger})\|^{\nu}.$

Thus, an application of the Young’s inequality gives

$\displaystyle D_{\mathcal{R}}^{\xi^{\dagger}}(x,x^{\dagger})$
$\displaystyle\leq{\mathcal{R}}(x)-{\mathcal{R}}(x^{\dagger})+\frac{1}{2}c_{0}\|x-x^{\dagger}\|^{2}+C_{\nu}\|\omega\|^{\frac{2}{1+\nu}}\|Ax-y\|^{\frac{2\nu}{1+\nu}},$

Finally, by invoking the strong convexity of ${\mathcal{R}}$ and (2.6) we can
obtain

$\displaystyle
D_{\mathcal{R}}^{\xi^{\dagger}}(x,x^{\dagger})\leq{\mathcal{R}}(x)-{\mathcal{R}}(x^{\dagger})+\frac{1}{2}D_{\mathcal{R}}^{\xi^{\dagger}}(x,x^{\dagger})+C_{\nu}\|\omega\|^{\frac{2}{1+\nu}}\|Ax-y\|^{\frac{2\nu}{1+\nu}}$

which implies the assertion.

It should be pointed out that the source condition (3.9) has been used in [11]
to derive convergence rates of Tikhonov regularization in Hilbert spaces with
non-quadratic penalty terms. We would also like to mention that for the
special case ${\mathcal{R}}(x)=\frac{1}{2}\|x\|^{2}+\iota_{\mathscr{C}}(x)$,
where $\iota_{\mathscr{C}}$ denotes the indicator function of a closed convex
set ${\mathscr{C}}\subset{\mathscr{X}}$, i.e. $\iota_{\mathscr{C}}(x)=0$ if
$x\in{\mathscr{C}}$ and $\infty$ otherwise, the source condition (3.9) becomes
the projected spectral source condition

$x^{\dagger}=P_{\mathscr{C}}((A^{*}A)^{\nu/2}\omega)$

for which it has been shown in [22] that the projected source condition
implies the variational source condition. Here $P_{\mathscr{C}}$ denotes the
metric projection of ${\mathscr{X}}$ onto ${\mathscr{C}}$.

######  Remark 3.4.

The variational source conditions have been verified for various concrete
inverse problems, see [6, 7, 17, 18] for instance.

######  Theorem 3.7.

Let Assumption 3.1 hold and consider the method (1.5) with
$\lambda_{0}^{(n)}=0$. Assume that $0<\gamma<2c_{0}/\|A\|^{2}$ and that
$x^{\dagger}$ satisfies Assumption 3.2. If $t_{n}$ is chosen such that
$t_{n}\sim(n/\sigma^{2})^{\frac{2-q}{2}}$, then

${\mathbb{E}}\left[{\mathcal{E}}^{\dagger}(x_{t_{n}}^{(n)})\right]=O\left(\left(\frac{\sigma^{2}}{n}\right)^{q/2}\right).$

###### Proof.

From Lemma 3.2 and (3.1) it follows that

$\displaystyle
c\gamma(t+1){\mathbb{E}}\left[\|Ax_{t}^{(n)}-\bar{y}^{(n)}\|^{2}\right]+\frac{1}{8\gamma(t+1)}{\mathbb{E}}\left[\|\lambda_{t+1}^{(n)}\|^{2}\right]$
$\displaystyle\leq\eta_{t}+\gamma(t+1){\mathbb{E}}\left[\|\bar{y}^{(n)}-y\|^{2}\right]=\eta_{t}+\gamma(t+1)\frac{\sigma^{2}}{n}.$

By using the variational source condition and the nonnegativity of
${\mathcal{E}}^{\dagger}$, we have

$0\leq{\mathcal{R}}(x)-{\mathcal{R}}(x^{\dagger})+M\|Ax-y\|^{q},\quad\forall
x\in{\mathscr{X}}.$

Therefore

$\displaystyle\eta_{t}$
$\displaystyle\leq\sup_{x\in{\mathscr{X}}}\left\\{M\|Ax-y\|^{q}-\frac{1}{3}\gamma(t+1)\|Ax-y\|^{2}\right\\}$
$\displaystyle\leq\sup_{s\geq
0}\left\\{Ms^{q}-\frac{1}{3}\gamma(t+1)s^{2}\right\\}=c_{2}(t+1)^{-\frac{q}{2-q}},$
(3.10)

where
$c_{2}:=(1-\frac{q}{2})\left(\frac{3qM}{2\gamma}\right)^{\frac{q}{2-q}}M$.
Consequently

$\displaystyle
c\gamma(t+1){\mathbb{E}}\left[\|Ax_{t}^{(n)}-\bar{y}^{(n)}\|^{2}\right]+\frac{1}{8\gamma(t+1)}{\mathbb{E}}\left[\|\lambda_{t+1}^{(n)}\|^{2}\right]$
$\displaystyle\leq
c_{2}(t+1)^{-\frac{q}{2-q}}+\gamma(t+1)\frac{\sigma^{2}}{n}.$

which shows that

$\displaystyle{\mathbb{E}}\left[\|Ax_{t}^{(n)}-\bar{y}^{(n)}\|^{2}\right]\leq
C\left((t+1)^{-\frac{2}{2-q}}+\frac{\sigma^{2}}{n}\right),$
$\displaystyle{\mathbb{E}}\left[\|\lambda_{t+1}^{(n)}\|^{2}\right]\leq
C\left((t+1)^{\frac{2(1-q)}{2-q}}+(t+1)^{2}\frac{\sigma^{2}}{n}\right),$

where here and below we use $C$ to denote a generic constant independent of
$t$, $n$ and $\sigma$. Since $t_{n}$ is chosen such that
$t_{n}\sim(n/\sigma^{2})^{\frac{2-q}{2}}$, we have

$\displaystyle{\mathbb{E}}\left[\|Ax_{t_{n}}^{(n)}-\bar{y}^{(n)}\|^{2}\right]\leq
C\frac{\sigma^{2}}{n}\quad\mbox{and}\quad{\mathbb{E}}\left[\|\lambda_{t_{n}}^{(n)}\|^{2}\right]\leq
C\left(\frac{n}{\sigma^{2}}\right)^{1-q}.$

The first estimate and (3.1) in particular imply

${\mathbb{E}}\left[\|Ax_{t_{n}}^{(n)}-y\|^{2}\right]\leq
2{\mathbb{E}}\left[\|Ax_{t_{n}}^{(n)}-\bar{y}^{(n)}\|^{2}\right]+2{\mathbb{E}}\left[\|\bar{y}^{(n)}-y\|^{2}\right]\leq
C\frac{\sigma^{2}}{n}.$

Now we are ready to complete the proof. By using the variational source
condition on $x^{\dagger}$, the convexity of ${\mathcal{R}}$, and the fact
$A^{*}\lambda_{t_{n}}^{(n)}\in\partial{\mathcal{R}}(x_{t_{n}}^{(n)})$ we have

$\displaystyle{\mathcal{E}}^{\dagger}(x_{t_{n}}^{(n)})$
$\displaystyle\leq{\mathcal{R}}(x_{t_{n}}^{(n)})-{\mathcal{R}}(x^{\dagger})+M\|Ax_{t_{n}}^{(n)}-y\|^{q}$
$\displaystyle\leq\langle
A^{*}\lambda_{t_{n}}^{(n)},x_{t_{n}}^{(n)}-x^{\dagger}\rangle+M\|Ax_{t_{n}}^{(n)}-y\|^{q}$
$\displaystyle=\langle\lambda_{t_{n}}^{(n)},Ax_{t_{n}}^{(n)}-y\rangle+M\|Ax_{t_{n}}^{(n)}-y\|^{q}$
$\displaystyle\leq\|\lambda_{t_{n}}^{(n)}\|\|Ax_{t_{n}}^{(n)}-y\|+M\|Ax_{t_{n}}^{(n)}-y\|^{q}.$
(3.11)

Therefore

$\displaystyle{\mathbb{E}}\left[{\mathcal{E}}^{\dagger}(x_{t_{n}}^{(n)})\right]$
$\displaystyle\leq{\mathbb{E}}\left[\|\lambda_{t_{n}}^{(n)}\|\|Ax_{t_{n}}^{(n)}-y\|\right]+M{\mathbb{E}}\left[\|Ax_{t_{n}}^{(n)}-y\|^{q}\right]$
$\displaystyle\leq\left({\mathbb{E}}\left[\|\lambda_{t_{n}}^{(n)}\|^{2}\right]\right)^{1/2}\left({\mathbb{E}}\left[\|Ax_{t_{n}}^{(n)}-y\|^{2}\right]\right)^{1/2}$
$\displaystyle\quad\,+M\left({\mathbb{E}}\left[\|Ax_{t_{n}}^{(n)}-y\|^{2}\right]\right)^{q/2}$
$\displaystyle\leq C\left(\frac{\sigma^{2}}{n}\right)^{q/2}.$

The proof is thus complete. ∎

###  3.2. Convergence analysis under a posteriori stopping rule

We next consider the dual gradient method (1.5) terminated by an a posteriori
stopping rule. The discrepancy principle is one of the most prominent rule
that has been studied extensively. In this stopping rule, one needs the
information on the noise level $\|\bar{y}^{(n)}-y\|$. Since $y$ is unknown, we
can not use this quantity directly. Recall that
${\mathbb{E}}\left[\|\bar{y}^{(n)}-y\|^{2}\right]=\sigma^{2}/n.$ To get an
estimate on $\sigma^{2}$, we consider the square root of the sample variance

$s_{n}:=\sqrt{\frac{1}{n-1}\sum_{i=1}^{n}\|y_{i}-\bar{y}^{(n)}\|^{2}}$

for which it is known that ${\mathbb{E}}[s_{n}^{2}]=\sigma^{2}$. Therefore, we
may use $s_{n}/\sqrt{n}$ as an estimator of $\|\bar{y}^{(n)}-y\|$. This leads
us to propose the following stopping rule which is a statistical variant of
the discrepancy principle.

######  Rule 3.8.

Let $\beta_{0}>0$ be a given number. We define $t_{n}$ to be the first integer
such that $t_{n}\leq\beta_{0}n$ and

$\displaystyle\|Ax_{t_{n}}^{(n)}-\bar{y}^{(n)}\|\leq\tau_{n}\frac{s_{n}}{\sqrt{n}},$
(3.12)

where $\tau_{n}>1$ is a number that may depend on $n$; if there is no such an
integer $t_{n}$ satisfying (3.12), we take $t_{n}=\beta_{0}n$.

Recall that we have used $s_{n}/\sqrt{n}$ as an estimator of
$\|\bar{y}^{(n)}-y\|$. In case $s_{n}/\sqrt{n}$ is an underestimated
estimator, (3.12) may not be satisfied at a right number of iterations;
continuing the iterations until (3.12) holds may result in a bad
reconstruction result. Therefore, the requirement $t_{n}\leq\beta_{0}n$ serves
as an emergency stop. It should be emphasized that $t_{n}$ is a random
variable instead of a deterministic quantity; this is a key difference between
Rule 3.8 and the classical discrepancy principle.

According to the spirit of the discrepancy principle, it is suggestive to take
$\tau_{n}>1$ in Rule 3.8 to be a number close to $1$. Due to the technical
reasons, however, in the forthcoming convergence analysis, the number
$\tau_{n}$ in Rule 3.8 is required to tend to $\infty$ as $n\to\infty$ but
does not go to $\infty$ too fast; specifically we require

$\displaystyle\tau_{n}\to\infty\quad\mbox{ and }\quad\tau_{n}/\sqrt{n}\to
0\quad\mbox{ as }n\to\infty.$ (3.13)

For instance, we may take

$\displaystyle\tau_{n}:=\max\\{\tau_{0},\log|\log|\log n||\\},$ (3.14)

where $\tau_{0}>1$ is a fixed number. This requirement on $\tau_{n}$ is
imposed to guarantee the existence of an event $\Omega_{n}$ with
${\mathbb{P}}(\Omega_{n})\to 1$ as $n\to\infty$ on which we can perform the
convergence analysis. Note that

$s_{n}^{2}=\frac{n}{n-1}\left(\frac{1}{n}\sum_{i=1}^{n}\|y_{i}\|^{2}-\|\bar{y}^{(n)}\|^{2}\right).$

Since $y_{1},y_{2},\cdots$ is a sequence of independent identically
distributed random variables, by the strong law of large numbers ([28,
Corollary 7.10]) we have

$\|\bar{y}^{(n)}\|^{2}\to\|y\|^{2}\quad\mbox{ and
}\quad\frac{1}{n}\sum_{i=1}^{n}\|y_{i}\|^{2}\to{\mathbb{E}}\left[\|y_{1}\|^{2}\right]$

as $n\to\infty$ almost surely. Therefore

$s_{n}^{2}\to{\mathbb{E}}[\|y_{1}\|^{2}]-\|y\|^{2}=\sigma^{2}$

and hence $s_{n}\to\sigma$ as $n\to\infty$ almost surely. Since almost sure
convergence implies convergence in probability, we have

$\displaystyle{\mathbb{P}}\left(|s_{n}-\sigma|\leq\frac{\sigma}{2}\right)\to
1\quad\mbox{as }n\to\infty.$ (3.15)

Now we define the event

$\Omega_{n}:=\left\\{|s_{n}-\sigma|\leq\frac{\sigma}{2},\,\|\bar{y}^{(n)}-y\|\leq\tau_{n}\sqrt{\frac{cs_{n}^{2}}{2n}}\right\\},$

where $c>0$ is the constant appearing in Lemma 3.2. We claim that

$\displaystyle{\mathbb{P}}(\Omega_{n})\to 1\quad\mbox{ as }n\to\infty.$ (3.16)

Indeed, by the Markov’s inequality, (3.1), and the assumption
$\tau_{n}\to\infty$, we have

$\displaystyle{\mathbb{P}}\left(|s_{n}-\sigma|\leq\frac{\sigma}{2},\,\|\bar{y}^{(n)}-y\|>\tau_{n}\sqrt{\frac{cs_{n}^{2}}{2n}}\right)$
$\displaystyle\leq{\mathbb{P}}\left(\|\bar{y}^{(n)}-y\|>\sigma\tau_{n}\sqrt{\frac{c}{8n}}\right)\leq\frac{8n}{c\sigma^{2}\tau_{n}^{2}}{\mathbb{E}}\left[\|\bar{y}^{(n)}-y\|^{2}\right]=\frac{8}{c\tau_{n}^{2}}\to
0$

as $n\to\infty$. Therefore, by using (3.15), we have

$\displaystyle{\mathbb{P}}(\Omega_{n})={\mathbb{P}}\left(|s_{n}-\sigma|\leq\frac{\sigma}{2}\right)-{\mathbb{P}}\left(|s_{n}-\sigma|\leq\frac{\sigma}{2},\,\|\bar{y}^{(n)}-y\|>\tau_{n}\sqrt{\frac{cs_{n}^{2}}{2n}}\right)\to
1$

as $n\to\infty$.

The following result gives an upper bound estimate of $t_{n}$ on $\Omega_{n}$
defined by Rule 3.8.

######  Lemma 3.9.

Let Assumption 3.1 hold and consider the method (1.5) with
$\lambda_{0}^{(n)}=0$. Assume that $0<\gamma<2c_{0}/\|A\|^{2}$. Let $t_{n}$ be
defined by Rule 3.8 with $\tau_{n}$ satisfying (3.13) and set
$\phi(n):=\frac{\sigma\tau_{n}}{2\sqrt{n}}$. Then for any $\varepsilon>0$
there is an integer $n_{\varepsilon}$ such that

$t_{n}\phi(n)^{2}\leq\varepsilon\quad\mbox{ on }\Omega_{n}$

for all $n\geq n_{\varepsilon}$.

###### Proof.

Note that on $\Omega_{n}$ we always have $\sigma/2\leq s_{n}\leq 3\sigma/2$.
If $t_{n}\leq 1/\phi(n)$, then $t_{n}\phi(n)^{2}\leq\phi(n)$. In the following
we will assume $t_{n}>1/\phi(n)$. By the definition of $t_{n}$ we have

$\|Ax_{t_{n}-1}^{(n)}-\bar{y}^{(n)}\|>\tau_{n}\frac{s_{n}}{\sqrt{n}}.$

Therefore, by using (3) we have on $\Omega_{n}$ that

$\displaystyle c\gamma t_{n}\frac{\tau_{n}^{2}s_{n}^{2}}{n}+\frac{1}{8\gamma
t_{n}}\|\lambda_{t_{n}}^{(n)}\|^{2}$ $\displaystyle\leq c\gamma
t_{n}\|Ax_{t_{n}-1}^{(n)}-\bar{y}^{(n)}\|^{2}+\frac{1}{8\gamma
t_{n}}\|\lambda_{t_{n}}^{(n)}\|^{2}$ $\displaystyle\leq
d_{y}(\lambda)-d_{y}(\lambda_{t_{n}}^{(n)})+\frac{3\|\lambda\|^{2}}{4\gamma
t_{n}}+\gamma t_{n}\|\bar{y}^{(n)}-y\|^{2}$ $\displaystyle\leq
d_{y}(\lambda)-d_{y}(\lambda_{t_{n}}^{(n)})+\frac{3\|\lambda\|^{2}}{4\gamma
t_{n}}+c\gamma t_{n}\frac{\tau_{n}^{2}s_{n}^{2}}{2n}$

which implies that

$\displaystyle\frac{c}{2}\gamma t_{n}\phi(n)^{2}+\frac{1}{8\gamma
t_{n}}\|\lambda_{t_{n}}^{(n)}\|^{2}$ $\displaystyle\leq c\gamma
t_{n}\frac{\tau_{n}^{2}s_{n}^{2}}{2n}+\frac{1}{8\gamma
t_{n}}\|\lambda_{t_{n}}^{(n)}\|^{2}$ $\displaystyle\leq
d_{y}(\lambda)-d_{y}(\lambda_{t_{n}}^{(n)})+\frac{3\|\lambda\|^{2}}{4\gamma
t_{n}}$ (3.17)

for all $\lambda\in{\mathscr{Y}}$. Note that

$d_{y}(\lambda)={\mathcal{R}}^{*}(A^{*}\lambda)-\langle
A^{*}\lambda,x^{\dagger}\rangle\geq-{\mathcal{R}}(x^{\dagger})$

which shows that

$\inf_{\lambda\in{\mathscr{Y}}}d_{y}(\lambda)\geq-{\mathcal{R}}(x^{\dagger})>-\infty.$

Therefore, for any $\varepsilon>0$ we can find
$\lambda_{\varepsilon}\in{\mathscr{Y}}$ such that

$d_{y}(\lambda_{\varepsilon})\leq\inf_{\lambda\in{\mathscr{Y}}}d_{y}(\lambda)+\frac{c\gamma}{4}\varepsilon\leq
d_{y}(\lambda_{t_{n}}^{(n)})+\frac{c\gamma}{4}\varepsilon.$

This and (3.2) with $\lambda=\lambda_{\varepsilon}$ show that

$\displaystyle\frac{c}{2}\gamma t_{n}\phi(n)^{2}+\frac{1}{8\gamma
t_{n}}\|\lambda_{t_{n}}^{(n)}\|^{2}\leq\frac{c\gamma}{4}\varepsilon+\frac{3\|\lambda_{\varepsilon}\|^{2}}{4\gamma
t_{n}}$ (3.18)

which together with $t_{n}>1/\phi(n)$ in particular implies

$\displaystyle
t_{n}\phi(n)^{2}\leq\frac{\varepsilon}{2}+\frac{3\|\lambda_{\varepsilon}\|^{2}}{2c\gamma^{2}t_{n}}\leq\frac{\varepsilon}{2}+\frac{3\|\lambda_{\varepsilon}\|^{2}}{2c\gamma^{2}}\phi(n).$

Therefore, we always have

$t_{n}\phi(n)^{2}\leq\max\left\\{\phi(n),\frac{\varepsilon}{2}+\frac{3\|\lambda_{\varepsilon}\|^{2}}{2c\gamma^{2}}\phi(n)\right\\}\quad\mbox{
on }\Omega_{n}.$

Since $\phi(n)\to 0$ as $n\to\infty$, we can find $n_{\varepsilon}$ such that
$\phi(n)\leq\varepsilon$ and
$\frac{3\|\lambda_{\varepsilon}\|^{2}}{2c\gamma^{2}}\phi(n)\leq\frac{\varepsilon}{2}$
for all $n\geq n_{\varepsilon}$. Thus $t_{n}\phi(n)^{2}\leq\varepsilon$ on
$\Omega_{n}$ for all $n\geq n_{\varepsilon}$. ∎

Based on Lemma 3.9, we are now ready to show a convergence result for
$x_{t_{n}}^{(n)}$ with $t_{n}$ determined by Rule 3.8.

######  Theorem 3.10.

Let Assumption 3.1 hold and consider the method (1.5) with
$\lambda_{0}^{(n)}=0$. Assume that $0<\gamma<2c_{0}/\|A\|^{2}$. Let $t_{n}$ be
determined by Rule 3.8 with $\tau_{n}$ satisfying (3.13). Assume that there
exists $\xi^{\dagger}\in\partial{\mathcal{R}}(x^{\dagger})$ such that
$\xi^{\dagger}\in\overline{A^{*}({\mathscr{Y}})}$. Then for any $0<\alpha<1$
there holds

$\displaystyle{\mathbb{E}}\left[\left(D_{\mathcal{R}}^{A^{*}\lambda_{t_{n}}^{(n)}}(x^{\dagger},x_{t_{n}}^{(n)})\right)^{\alpha}\right]\to
0\quad\mbox{ as }n\to\infty$ (3.19)

and consequently for any $0<p<2$ there holds

$\displaystyle{\mathbb{E}}\left[\|x_{t_{n}}^{(n)}-x^{\dagger}\|^{p}\right]\to
0\quad\mbox{ as }n\to\infty.$ (3.20)

###### Proof.

First we have

$\displaystyle
D_{\mathcal{R}}^{A^{*}\lambda_{t_{n}}^{(n)}}(x^{\dagger},x_{t_{n}}^{(n)})$
$\displaystyle\leq
D_{\mathcal{R}}^{A^{*}\lambda_{t_{n}}^{(n)}}(x^{\dagger},x_{t_{n}}^{(n)})+D_{\mathcal{R}}^{\xi^{\dagger}}(x_{t_{n}}^{(n)},x^{\dagger})$
$\displaystyle=\langle
A^{*}\lambda_{t_{n}}^{(n)}-\xi^{\dagger},x_{t_{n}}^{(n)}-x^{\dagger}\rangle.$
(3.21)

Since $\xi^{\dagger}\in\overline{A^{*}({\mathscr{Y}})}$, for any
$\varepsilon>0$ we can find $\lambda^{\dagger}_{\varepsilon}\in{\mathscr{Y}}$
and $v_{\varepsilon}\in{\mathscr{X}}^{*}$ such that

$\xi^{\dagger}=A^{*}\lambda^{\dagger}_{\varepsilon}+v_{\varepsilon}\quad\mbox{
and }\quad\|v_{\varepsilon}\|^{2}\leq\frac{1}{4}c_{0}\varepsilon.$

Therefore

$\displaystyle
D_{\mathcal{R}}^{A^{*}\lambda_{t_{n}}^{(n)}}(x^{\dagger},x_{t_{n}}^{(n)})$
$\displaystyle\leq\langle
A^{*}\lambda_{t_{n}}^{(n)}-A^{*}\lambda^{\dagger}_{\varepsilon}-v_{\varepsilon},x_{t_{n}}^{(n)}-x^{\dagger}\rangle$
$\displaystyle\leq\langle\lambda_{t_{n}}^{(n)}-\lambda^{\dagger}_{\varepsilon},Ax_{t_{n}}^{(n)}-y\rangle+\|v_{\varepsilon}\|\|x_{t_{n}}^{(n)}-x^{\dagger}\|.$

By using the strong convexity of ${\mathcal{R}}$ and the estimate on
$\|v_{\varepsilon}\|$ we have

$\displaystyle\|v_{\varepsilon}\|\|x_{t_{n}}^{(n)}-x^{\dagger}\|\leq\frac{1}{2}c_{0}\|x_{t_{n}}^{(n)}-x^{\dagger}\|^{2}+\frac{\|v_{\varepsilon}\|^{2}}{2c_{0}}\leq\frac{1}{2}D_{\mathcal{R}}^{A^{*}\lambda_{t_{n}}^{(n)}}(x^{\dagger},x_{t_{n}}^{(n)})+\frac{\varepsilon}{8}.$

Consequently

$\displaystyle
D_{\mathcal{R}}^{A^{*}\lambda_{t_{n}}^{(n)}}(x^{\dagger},x_{t_{n}}^{(n)})$
$\displaystyle\leq
2\langle\lambda_{t_{n}}^{(n)}-\lambda^{\dagger}_{\varepsilon},Ax_{t_{n}}^{(n)}-y\rangle+\frac{\varepsilon}{4}$
$\displaystyle\leq
2\left(\|\lambda_{t_{n}}^{(n)}\|+\|\lambda^{\dagger}_{\varepsilon}\|\right)\|Ax_{t_{n}}^{(n)}-y\|+\frac{\varepsilon}{4}.$
(3.22)

We now show that there exists $n_{\varepsilon}$ such that

$\displaystyle
D_{\mathcal{R}}^{A^{*}\lambda_{t_{n}}^{(n)}}(x^{\dagger},x_{t_{n}}^{(n)})\leq\varepsilon\quad\mbox{
on }\Omega_{n}$ (3.23)

for all $n\geq n_{\varepsilon}$. According to Lemma 3.9 and (3.13) we have
$t_{n}<1/\phi(n)^{2}\leq\beta_{0}n$ on $\Omega_{n}$ for sufficiently large
$n$. Therefore

$\|Ax_{t_{n}}^{(n)}-\bar{y}^{(n)}\|\leq\tau_{n}\frac{s_{n}}{\sqrt{n}}$

and consequently

$\displaystyle\|Ax_{t_{n}}^{(n)}-y\|$
$\displaystyle\leq\|Ax_{t_{n}}^{(n)}-\bar{y}^{(n)}\|+\|\bar{y}^{(n)}-y\|\leq\tau_{n}\frac{s_{n}}{\sqrt{n}}+\tau_{n}\sqrt{\frac{cs_{n}^{2}}{2n}}$
$\displaystyle\leq\frac{3\sigma\tau_{n}}{2\sqrt{n}}+\frac{3\sigma\tau_{n}}{2}\sqrt{\frac{c}{2n}}=3\left(1+\sqrt{\frac{c}{2}}\right)\phi(n)$
$\displaystyle<\frac{9}{2}\phi(n).$ (3.24)

By using (3.18) we have

$\|\lambda_{t_{n}}^{(n)}\|^{2}\leq 2c\gamma^{2}\varepsilon
t_{n}+6\|\lambda_{\varepsilon}\|^{2},$

where $\lambda_{\varepsilon}\in{\mathscr{Y}}$ is an element chosen in the
proof of Lemma 3.9. Therefore, it follows from (3.2) that

$\displaystyle
D_{\mathcal{R}}^{A^{*}\lambda_{t_{n}}^{(n)}}(x^{\dagger},x_{t_{n}}^{(n)})\leq
9\left(\sqrt{2c\gamma^{2}\varepsilon
t_{n}+6\|\lambda_{\varepsilon}\|^{2}}+\|\lambda^{\dagger}_{\varepsilon}\|\right)\phi(n)+\frac{\varepsilon}{4}$
(3.25)

on $\Omega_{n}$. By using Lemma 3.9 and the property $\phi(n)\to 0$ as
$n\to\infty$, we may find a sufficiently large $n_{\varepsilon}$ such that

$\left(\sqrt{2c\gamma^{2}\varepsilon
t_{n}+6\|\lambda_{\varepsilon}\|^{2}}+\|\lambda^{\dagger}_{\varepsilon}\|\right)\phi(n)\leq\frac{\varepsilon}{12}$

for all $n\geq n_{\varepsilon}$. Combining this with (3.25) we thus obtain
(3.23).

Next we will show (3.19). Let $\chi_{\Omega_{n}}$ denote the characteristic
function of $\Omega_{n}$, i.e. $\chi_{\Omega_{n}}(\omega)=1$ if
$\omega\in\Omega_{n}$ and $\chi_{\Omega_{n}}(\omega)=0$ otherwise. By using
(3.23) and the Hölder inequality we have for sufficiently large $n$ that

$\displaystyle{\mathbb{E}}\left[\left(D_{\mathcal{R}}^{A^{*}\lambda_{t_{n}}^{(n)}}(x^{\dagger},x_{t_{n}}^{(n)})\right)^{\alpha}\right]$
$\displaystyle={\mathbb{E}}\left[\left(D_{\mathcal{R}}^{A^{*}\lambda_{t_{n}}^{(n)}}(x^{\dagger},x_{t_{n}}^{(n)})\right)^{\alpha}\chi_{\Omega_{n}}\right]+{\mathbb{E}}\left[\left(D_{\mathcal{R}}^{A^{*}\lambda_{t_{n}}^{(n)}}(x^{\dagger},x_{t_{n}}^{(n)})\right)^{\alpha}\chi_{\Omega_{n}^{c}}\right]$
$\displaystyle\leq\varepsilon^{\alpha}+\left({\mathbb{E}}\left[D_{\mathcal{R}}^{A^{*}\lambda_{t_{n}}^{(n)}}(x^{\dagger},x_{t_{n}}^{(n)})\right]\right)^{\alpha}{\mathbb{P}}(\Omega_{n}^{c})^{1-\alpha}.$
(3.26)

If we are able to show

$\displaystyle{\mathbb{E}}\left[D_{\mathcal{R}}^{A^{*}\lambda_{t_{n}}^{(n)}}(x^{\dagger},x_{t_{n}}^{(n)})\right]\leq
C$ (3.27)

for some constant $C$ independent of $n$, then, by using the fact that
${\mathbb{P}}(\Omega_{n}^{c})\to 0$ as $n\to\infty$, we can conclude

$\limsup_{n\to\infty}{\mathbb{E}}\left[\left(D_{\mathcal{R}}^{A^{*}\lambda_{t_{n}}^{(n)}}(x^{\dagger},x_{t_{n}}^{(n)})\right)^{\alpha}\right]\leq\varepsilon^{\alpha}$

and thus obtain (3.19), due to the arbitrariness of $\varepsilon$.

It remains only to show (3.27). From (3.2), the Cauchy-Schwarz inequality and
the strong convexity of ${\mathcal{R}}$ it follows that

$\displaystyle
D_{\mathcal{R}}^{A^{*}\lambda_{t_{n}}^{(n)}}(x^{\dagger},x_{t_{n}}^{(n)})$
$\displaystyle\leq\langle\lambda_{t_{n}}^{(n)},Ax_{t_{n}}^{(n)}-y\rangle+\|\xi^{\dagger}\|\|x_{t_{n}}^{(n)}-x^{\dagger}\|$
$\displaystyle\leq\|\lambda_{t_{n}}^{(n)}\|\|Ax_{t_{n}}^{(n)}-y\|+\frac{1}{2c_{0}}\|\xi^{\dagger}\|^{2}+\frac{1}{2}D_{\mathcal{R}}^{A^{*}\lambda_{t_{n}}^{(n)}}(x^{\dagger},x_{t_{n}}^{(n)})$

which implies

$\displaystyle
D_{\mathcal{R}}^{A^{*}\lambda_{t_{n}}^{(n)}}(x^{\dagger},x_{t_{n}}^{(n)})$
$\displaystyle\leq
2\|\lambda_{t_{n}}^{(n)}\|\|Ax_{t_{n}}^{(n)}-y\|+\frac{1}{c_{0}}\|\xi^{\dagger}\|^{2}$

We next use Lemma 3.2. By the definition of $\eta_{t}$ and the nonnegativity
of ${\mathcal{R}}$, we have
$0\leq\eta_{t}\leq{\mathcal{R}}(x^{\dagger})<\infty$ for all $t\geq 0$. Thus,
it follows from Lemma 3.2 that

$\displaystyle\|Ax_{t_{n}}^{(n)}-\bar{y}^{(n)}\|^{2}\leq
C\left(\frac{1}{t_{n}}+\|\bar{y}^{(n)}-y\|^{2}\right),\quad\|\lambda_{t_{n}}^{(n)}\|^{2}\leq
C\left(t_{n}+t_{n}^{2}\|\bar{y}^{(n)}-y\|^{2}\right).$

Therefore, by using $t_{n}\leq\beta_{0}n$, we have

$\displaystyle
D_{\mathcal{R}}^{A^{*}\lambda_{t_{n}}^{(n)}}(x^{\dagger},x_{t_{n}}^{(n)})$
$\displaystyle\leq
C\left(\sqrt{t_{n}}+t_{n}\|\bar{y}^{(n)}-y\|\right)\left(\frac{1}{\sqrt{t_{n}}}+\|\bar{y}^{(n)}-y\|\right)+\frac{1}{c_{1}}\|\xi^{\dagger}\|^{2}$
$\displaystyle\leq
C\left(1+\sqrt{t_{n}}\|\bar{y}^{(n)}-y\|+t_{n}\|\bar{y}^{(n)}-y\|^{2}\right)+\frac{1}{c_{1}}\|\xi^{\dagger}\|^{2}$
$\displaystyle\leq
C\left(1+\sqrt{n}\|\bar{y}^{(n)}-y\|+n\|\bar{y}^{(n)}-y\|^{2}\right)+\frac{1}{c_{1}}\|\xi^{\dagger}\|^{2}.$
(3.28)

By taking the expectation and using (3.1), we can obtain

$\displaystyle{\mathbb{E}}\left[D_{\mathcal{R}}^{A^{*}\lambda_{t_{n}}^{(n)}}(x^{\dagger},x_{t_{n}}^{(n)})\right]$
$\displaystyle\leq
C\left(1+\sqrt{n}({\mathbb{E}}[\|\bar{y}^{(n)}-y\|^{2}])^{1/2}+n{\mathbb{E}}[\|\bar{y}^{(n)}-y\|^{2}]\right)+\frac{1}{c_{1}}\|\xi^{\dagger}\|^{2}$
$\displaystyle\leq
C\left(1+\sigma+\sigma^{2}\right)+\frac{1}{c_{1}}\|\xi^{\dagger}\|^{2}$

which shows (3.27). The proof is therefore complete. ∎

######  Remark 3.5.

In Theorem 3.10 we have obtained the convergence result (3.19) with
$0<\alpha<1$. If in addition

$\sigma_{4}:={\mathbb{E}}[\|y_{1}-y\|^{4}]<\infty,$

we can improve the convergence result (3.19) to include $\alpha=1$ and hence
${\mathbb{E}}[\|x_{t_{n}}^{(n)}-x^{\dagger}\|^{2}]\to 0$ as $n\to\infty$. To
see this, we may use the similar argument for deriving (3.2) to obtain

$\displaystyle{\mathbb{E}}\left[D_{\mathcal{R}}^{A^{*}\lambda_{t_{n}}^{(n)}}(x^{\dagger},x_{t_{n}}^{(n)})\right]\leq\varepsilon+\left({\mathbb{E}}\left[\left(D_{\mathcal{R}}^{A^{*}\lambda_{t_{n}}^{(n)}}(x^{\dagger},x_{t_{n}}^{(n)})\right)^{2}\right]\right)^{1/2}{\mathbb{P}}(\Omega_{n}^{c})^{1/2}.$

Therefore, it suffices to show

$\displaystyle{\mathbb{E}}\left[\left(D_{\mathcal{R}}^{A^{*}\lambda_{t_{n}}^{(n)}}(x^{\dagger},x_{t_{n}}^{(n)})\right)^{2}\right]\leq
C$ (3.29)

for some constant $C$ independent of $n$. By virtue of (3.2) it is easy to
derive

$\displaystyle{\mathbb{E}}\left[\left(D_{\mathcal{R}}^{A^{*}\lambda_{t_{n}}^{(n)}}(x^{\dagger},x_{t_{n}}^{(n)})\right)^{2}\right]$
$\displaystyle\leq
C\left(1+n{\mathbb{E}}\left[\|\bar{y}^{(n)}-y\|^{2}\right]+n^{2}{\mathbb{E}}\left[\|\bar{y}^{(n)}-y\|^{4}\right]\right)+\frac{2}{c_{1}}\|\xi^{\dagger}\|^{2}.$
(3.30)

With the similar argument in the proof of [13, Corollary 3], we can derive

$\displaystyle{\mathbb{E}}[\|\bar{y}^{(n)}-y\|^{4}]$
$\displaystyle\leq\frac{1}{n^{3}}{\mathbb{E}}\left[\|y_{1}-y\|^{4}\right]+\frac{3(n-1)}{n^{3}}\left({\mathbb{E}}\left[\|y_{1}-y\|^{2}\right]\right)^{2}$
$\displaystyle=\frac{\sigma_{4}}{n^{3}}+\frac{3(n-1)}{n^{3}}\sigma^{4}.$

By invoking this estimate and (3.1), we can obtain (3.29) from (3.5)
immediately.

######  Theorem 3.11.

Let Assumption 3.1 hold and consider the method (1.5) with
$\lambda_{0}^{(n)}=0$. Assume that $0<\gamma<2c_{0}/\|A\|^{2}$. Let $t_{n}$ be
determined by Rule 3.8 with $\tau_{n}$ satisfying (3.13). If $x^{\dagger}$
satisfies the variational source condition given in Assumption 3.2, then there
is a constant $C$ such that

${\mathbb{P}}\left({\mathcal{E}}^{\dagger}(x_{t_{n}}^{(n)})\leq
C\phi(n)^{q}\right)\to 1$

as $n\to\infty$, where $\phi(n):=\frac{\sigma\tau_{n}}{2\sqrt{n}}$ is defined
in Lemma 3.9.

###### Proof.

Since ${\mathbb{P}}(\Omega_{n})\to 1$ as $n\to\infty$, it suffices to
establish

${\mathcal{E}}^{\dagger}(x_{t_{n}}^{(n)})\leq C\phi(n)^{q}\quad\mbox{ on
}\Omega_{n}$

for some constant $C$ independent of $n$. Under the given variational source
condition on $x^{\dagger}$, we have the estimate (3.1) on $\eta_{t}$.
Combining this with Lemma 3.2 shows that

$\displaystyle c\gamma\|Ax_{t_{n}-1}^{(n)}-\bar{y}^{(n)}\|^{2}\leq
c_{2}t_{n}^{-\frac{2}{2-q}}+\gamma\|\bar{y}^{(n)}-y\|^{2}\leq
c_{2}t_{n}^{-\frac{2}{2-q}}+c\gamma\frac{\tau_{n}^{2}s_{n}^{2}}{2n}$ (3.31)

and

$\displaystyle\|\lambda_{t_{n}}^{(n)}\|^{2}\leq 8\gamma
t_{n}^{\frac{2(1-q)}{2-q}}+8\gamma^{2}t_{n}^{2}\|\bar{y}^{(n)}-y\|^{2}\leq
8\gamma
t_{n}^{\frac{2(1-q)}{2-q}}+4\gamma^{2}t_{n}^{2}\frac{\tau_{n}^{2}s_{n}^{2}}{n}$
(3.32)

on $\Omega_{n}$. Since
$\|Ax_{t_{n}-1}^{(n)}-\bar{y}^{(n)}\|>\tau_{n}\frac{s_{n}}{\sqrt{n}}$, we have
from (3.31) that

$c\gamma\frac{\tau_{n}^{2}s_{n}^{2}}{n}\leq
c_{2}t_{n}^{-\frac{2}{2-q}}+c\gamma\frac{\tau_{n}^{2}s_{n}^{2}}{2n}$

which together with $s_{n}\geq\sigma/2$ on $\Omega_{n}$ implies

$\frac{1}{2}c\gamma\phi(n)^{2}\leq c_{2}t_{n}^{-\frac{2}{2-q}}.$

Therefore

$\displaystyle t_{n}\leq c_{3}\phi(n)^{q-2}$ (3.33)

where $c_{3}>0$ is a constant independent of $n$. This in particular shows
that $t_{n}<\beta_{0}n$ and hence
$\|Ax_{t_{n}}^{(n)}-\bar{y}^{(n)}\|\leq\tau_{n}\frac{s_{n}}{\sqrt{n}}$ on
$\Omega_{n}$ for sufficiently large $n$. Consequently

$\|Ax_{t_{n}}^{(n)}-y\|\leq\frac{9}{2}\phi(n)$

as argued in (3.2). By using (3.32), (3.33) and $s_{n}\leq 3\sigma/2$ we also
have

$\displaystyle\|\lambda_{t_{n}}^{(n)}\|^{2}\leq 8\gamma
t_{n}^{\frac{2(1-q)}{2-q}}+4c\gamma^{2}t_{n}^{2}\phi(n)^{2}\leq
c_{4}\phi(n)^{2(q-1)}$

for some constant $c_{4}$ independent of $n$. Finally we may use (3.1) to
obtain

$\displaystyle{\mathcal{E}}^{\dagger}(x_{t_{n}}^{(n)})$
$\displaystyle\leq\|\lambda_{t_{n}}^{(n)}\|\|Ax_{t_{n}}^{(n)}-y\|+C\|Ax_{t_{n}}^{(n)}-y\|^{q}\leq
c_{5}\phi(n)^{q}$

for some constant $c_{5}$ independent of $n$. The proof is complete. ∎

######  Remark 3.6.

According to the definition of $\phi(n)$, Theorem 3.11 gives the convergence
rate

${\mathbb{P}}\left({\mathcal{E}}(x_{t_{n}}^{(n)})\leq
C\tau_{n}^{q}\left(\frac{\sigma^{2}}{n}\right)^{q/2}\right)\to 1\quad\mbox{ as
}n\to\infty.$

Because $\tau_{n}\to\infty$, this rate is worse than the one derived in
Theorem 3.7 under an a priori stopping rule. The requirement
$\tau_{n}\to\infty$ is used to construct $\Omega_{n}$ with
${\mathbb{P}}(\Omega_{n})\to 1$ as $n\to\infty$. If such an even $\Omega_{n}$
could be constructed without using the requirement on $\tau_{n}$, the
convergence rate in Theorem 3.11 could be upgraded to
${\mathbb{P}}({\mathcal{E}}(x_{t_{n}}^{(n)})\leq C(\sigma^{2}/n)^{q/2})\to 1$
as $n\to\infty$. This question however remains open.

##  4\. Numerical results

In this section we will report various numerical results to test the
performance of the method (1.5).

######  Example 4.1.

We first consider the application of the method (1.5) to solve linear ill-
posed problems in Hilbert spaces with convex constraint. Let
$A:{\mathscr{X}}\to{\mathscr{Y}}$ be a bounded linear operator between two
Hilbert spaces ${\mathscr{X}}$ and ${\mathscr{Y}}$ and let
${\mathscr{C}}\subset{\mathscr{X}}$ be a closed convex set. Given $y\in
A({\mathscr{C}})$, we consider finding the unique solution $x^{\dagger}$ of
$Ax=y$ in ${\mathscr{C}}$ with minimal norm which can be stated as (1.2) with

${\mathcal{R}}(x):=\frac{1}{2}\|x\|^{2}+\iota_{\mathscr{C}}(x),$

where $\iota_{\mathscr{C}}$ denotes the indicator function of ${\mathscr{C}}$.
Clearly ${\mathcal{R}}$ satisfies Assumption 3.1 (ii) with $c_{0}=1/2$. It is
easy to see that the method (1.5) takes the form

$\displaystyle
x_{t}^{(n)}=P_{\mathscr{C}}(A^{*}\lambda_{t}^{(n)}),\qquad\lambda_{t+1}^{(n)}=\lambda_{t}^{(n)}-\gamma(Ax_{t}^{(n)}-\bar{y}^{(n)}),$
(4.1)

where $P_{\mathscr{C}}$ denotes the metric projection of ${\mathscr{X}}$ onto
${\mathscr{C}}$. In case ${\mathscr{X}}=L^{2}(\mathscr{D})$ for some domain
$\mathscr{D}\subset{\mathbb{R}}^{d}$ and
${\mathscr{C}}=\\{x\in{\mathscr{X}}:x\geq 0\mbox{ a.e. on }\mathscr{D}\\}$,
the iteration scheme (4.1) becomes

$\displaystyle
x_{t}^{(n)}=\max\\{A^{*}\lambda_{t}^{(n)},0\\},\qquad\lambda_{t+1}^{(n)}=\lambda_{t}^{(n)}-\gamma(Ax_{t}^{(n)}-\bar{y}^{(n)})$
(4.2)

with initial guess $\lambda_{0}^{(n)}=0$. This method can be viewed as an
application of the Landweber iteration to the dual variable with nonnegative
constraint on the primal variable.

Figure 1. The noisy data (10 samples, green circle), the exact data (blue one)
and the reconstruction results with different sample sizes of Example 4.1,
where the data are corrupted by Gaussian noise.

We now test the performance of the method (4.2) by reconstructing nonnegative
solutions of linear ill-posed problems and compare it with the Landweber
iteration

$\displaystyle x_{t+1}^{(n)}=x_{t}^{(n)}-\gamma
A^{*}(Ax_{t}^{(n)}-\bar{y}^{(n)})$ (4.3)

that is considered in [13] which does not incorporate the nonnegative
constraint. Let us consider the first kind Fredholm integral equation of the
form

$Ax(s)=\int_{0}^{1}k(s,t)x(t)dt=y(s)\quad\mbox{ on }[0,1]$ (4.4)

with the kernel

$\displaystyle k(s,t)=\left\\{\begin{array}[]{lll}40s(1-t),&s\leq t\\\
40t(1-s),&s\geq t.\end{array}\right.$ (4.7)

It is easy to see that $A$ is a compact linear operator from
${\mathscr{X}}=L^{2}[0,1]$ to ${\mathscr{Y}}=L^{2}[0,1]$. We assume the sought
solution is given by

$\displaystyle
x^{\dagger}(s)=\left\\{\begin{array}[]{lll}20s(s-0.2)(0.7-s),&0.2\leq s\leq
0.7,\\\ 0,&\mbox{elsewhere}\end{array}\right.$

which is nonnegative and the exact data is $y:=Ax^{\dagger}$. By adding
independent Gaussian noise $\epsilon_{i}$ in $N(0,\sigma^{2})$ with
$\sigma=0.2$ to $y$ we produce the independent identically distributed noisy
data $y_{i}=y+\epsilon_{i}$, $i=1,\cdots,n$. For implementing (4.2) and (4.3)
with the noisy data $y_{i}$ we use the step-size $\gamma=2/\|A\|^{2}$ and
terminate the iterations by Rule 3.8 with $\beta_{0}=10$ and $\tau_{n}$ given
by (3.14), where $\tau_{0}=1.1$. In our numerical simulations, we divide
$[0,1]$ into $m=400$ subintervals of equal length and approximate integrals by
the trapezoidal rule. In Figure 1 we plot the exact data (blue one) and the
noisy data (10 samples, green circles). We consider several different sample
sizes $n=10,10^{2},10^{3},10^{4},10^{5}$, each is run by 200 simulations. In
order to visualize the performance, we plot in Figure 1 the reconstructed
solutions (the mean of 200 simulations), where “dgm-NN” and “Landweber”
represent the results obtained by the methods (4.2) and (4.3) respectively.
The results indicate that, as the sample size increases, more accurate
reconstruction results can be obtained. Since the method (4.2) incorporates
the nonnegativity constraint, it produces satisfactory results; while the
reconstruction by (4.3) includes undesired negative values.

Table 1. Numerical results of Example 4.1, where the data are corrupted by Gaussian noise. n | Method | Iteration numbers | Relative error | Emergency stops
---|---|---|---|---
|  | (in average) |  |
$10$ | dgm-NN | 17 | 3.9138e-01 | 0
| Landweber | 17 | 3.9254e-01 | 0
$10^{2}$ | dgm-NN | 115 | 1.2506e-01 | 0
| Landweber | 101 | 1.5405e-01 | 0
$10^{3}$ | dgm-NN | 252 | 7.4015e-02 | 0
| Landweber | 356 | 9.7710e-02 | 0
$10^{4}$ | dgm-NN | 887 | 4.3395e-02 | 0
| Landweber | 1174 | 6.1918e-02 | 0
$10^{5}$ | dgm-NN | 2686 | 2.2675e-02 | 0
| Landweber | 4373 | 3.7267e-02 | 0

In Table 1 we report the numerical results including the required average
number of iterations and the relative error (in average) which is calculated
by $\sqrt{\frac{1}{200}\sum_{i=1}^{200}e_{i}^{2}}$ as an approximation of
$({\mathbb{E}}[\|x_{t_{n}}^{(n)}-x^{\dagger}\|^{2}/\|x^{\dagger}\|^{2}])^{1/2}$,
where $e_{i}$ denotes the relative error between the exact solution and the
reconstructed solution of the $i$-th run respectively. We also record the
number of the times that $t_{n}$ reaches $\beta_{0}n$ to see if Rule 3.8
requires the emergency stop.

Figure 2.  Boxplots of the relative errors for 200 simulations for Example 4.1
with different sample sizes for noisy data corrupted by Gaussan noise
$N(0,\sigma^{2})$ with $\sigma=0.2$; the left one is for the method (4.3) and
the right one is for the method (4.2).

In Figure 2 we present the boxplots of the relative errors given by the
methods (4.2) and (4.3) with different sample sizes. On each box, the central
mark is the median, the bottle and top edges of the box indicate the 25th and
75th percentiles, the whiskers extend to the most extreme data points the
algorithm considers to be not outliers, and the outliers (red crosses) are
plotted individually. It is visible that the proposed method is convergent.
The red crosses below the blue box implies that the real noise levels have
been underestimated. In this case the upper bound $\beta_{0}n$ of iteration
numbers plays the important role of emergency stop. On the other hand, the red
crosses above the blue box implies the real noise levels have been
overestimated or the semi-convergence phenomenon has been happened already.

######  Example 4.2.

Consider the equation $Ax=y$, where $A:L^{1}(\mathscr{D})\to{\mathscr{Y}}$ is
a bounded linear operator, ${\mathscr{Y}}$ is a Hilbert space, and
$\mathscr{D}\subset{\mathbb{R}}^{d}$ is a bounded domain. Assuming the sought
solution is a probability density function, we may find such a solution by
considering the convex minimization problem (1.2) with

$\displaystyle{\mathcal{R}}(x):=f(x)+\iota_{\Delta}(x),$ (4.8)

where $\iota_{\Delta}$ denotes the indicator function of the closed convex set

$\displaystyle\Delta:=\left\\{x\in L^{1}(\mathscr{D}):x\geq 0\mbox{ a.e. on
}\mathscr{D}\mbox{ and }\int_{\mathscr{D}}x=1\right\\}$

in $L^{1}(\mathscr{D})$ and $f$ denotes the negative of the Boltzmann-Shannon
entropy, i.e.

$f(x):=\left\\{\begin{array}[]{lll}\int_{\mathscr{D}}x\log x&\mbox{ if }x\in
L_{+}^{1}(\mathscr{D})\mbox{ and }x\log x\in
L^{1}(\mathscr{D}),\\\\[5.16663pt] \infty&\mbox{
otherwise},\end{array}\right.$

where $L_{+}^{1}(\mathscr{D}):=\\{x\in L^{1}(\mathscr{D}):x\geq 0\mbox{ a.e.
on }\mathscr{D}\\}$. According to [3], ${\mathcal{R}}$ satisfies Assumption
3.1 (ii) with $c_{0}=1/2$. By the Karush-Kuhn-Tucker theory, for any $\ell\in
L^{\infty}(\mathscr{D})$ the unique minimizer of

$\displaystyle\min_{x\in
L^{1}(\mathscr{D})}\left\\{{\mathcal{R}}(x)-\int_{\mathscr{D}}\ell x\right\\}$

is given by $\hat{x}:=e^{\ell}/\int_{\mathscr{D}}e^{\ell}$. Therefore the dual
gradient method (1.5) using multiple repeated measurement data takes the form

$\displaystyle
x_{t}^{(n)}=\frac{1}{\int_{\Omega}e^{A^{*}\lambda_{t}^{(n)}}}e^{A^{*}\lambda_{t}^{(n)}},\qquad\lambda_{t+1}^{(n)}=\lambda_{t}^{(n)}-\gamma(Ax_{t}^{(n)}-\bar{y}^{(n)})$
(4.9)

with initial guess $\lambda_{0}^{(n)}=0$, which is an entropic dual gradient
method ([22]). This method is actually equivalent to the exponentiated gradent
method ([5, 27])

$x_{t+1}^{(n)}=\frac{x_{t}^{(n)}e^{\gamma
A^{*}(\bar{y}^{(n)}-Ax_{t}^{(n)})}}{\int_{\Omega}x_{t}^{(n)}e^{\gamma
A^{*}(\bar{y}^{(n)}-Ax_{t}^{(n)})}}.$

Figure 3. The noisy data (10 samples, green circle), the exact data (blue one)
and the reconstruction results with different sample sizes for Example 4.2.

In the numerical experiment, we consider again ill-posed integral equation
(4.4) with kernel

$k(s,t)=3e^{-\frac{(s-t)^{2}}{0.04}}$

and exact solution

$x^{\dagger}(s)=c\left(e^{-\frac{(s-0.3)^{2}}{0.01}}+3e^{-\frac{(s-0.7)^{2}}{0.005}}\right),$

where $c$ is chosen to ensure $\int_{0}^{1}x(s)ds=1$ so that $x^{\dagger}$ is
a probability density function. Clearly the integral operator is a compact
linear operator from $L^{1}[0,1]$ to $L^{2}[0,1]$. To perform numerical
simulations, we divide $[0,1]$ into $m=400$ subintervals of equal length and
approximate integrals by the trapezoidal rule. In Figure 3, we plot the exact
data (the blue one) and the noisy data (10 samples, green circle); these noisy
data are produced from the exact data $y=Ax^{\dagger}$ by adding independent
Gaussian noise in $N(0,\sigma^{2})$ with $\sigma=0.4$.

Table 2. Numerical results for Example 4.2, where the sought solution is a probability density function. n | Iteration numbers | Relative error | Emergency stops
---|---|---|---
| (in average) |  |
$10$ | 51 | 4.2691e-01 | 2
$10^{2}$ | 161 | 1.7945e-01 | 0
$10^{3}$ | 370 | 1.0920e-01 | 0
$10^{4}$ | 1142 | 9.7181e-02 | 0
$10^{5}$ | 3418 | 9.6206e-02 | 0

When implementing the method (4.9), we use the step-size $\gamma=2/\|A\|^{2}$
and terminate the iteration by Rule 3.8 with $\beta_{0}=10$ and $\tau_{n}$
given by (3.14) with $\tau_{0}=1.1$. We perform the numerical computation for
several different sample sizes $n=10,10^{2},10^{3},10^{4},10^{5}$ each is ran
by 200 simulations. The reconstruction results in average are plotted in
Figure 3. In Table 2 we report further computational results, including the
averaged number of iterations, the $L^{1}$ relative error in average, and the
utilized number of emergence stops. We also present the boxplots of the
$L^{1}$ relative errors for 200 simulations in Figure 4. The results indicate
that, when sample size increases, the relative error is reduced and more
accurate reconstruction result can be produced.

Figure 4.  Boxplots of the relative errors for 200 simulations with different
sample sizes for Example 4.2.

######  Example 4.3.

In this example we consider using the method (1.5) to reconstruct sparse
solutions for ill-posed problems $Af=u$, where $A:L^{2}(\mathscr{D})\to
L^{2}(\mathscr{D})$ is a bounded linear operator and
$\mathscr{D}\subset{\mathbb{R}}^{d}$ is a bounded domain. For this purpose, we
take ${\mathcal{R}}(f)$ to be a strongly convex perturbation of
$\|f\|_{L^{1}(\mathscr{D})}$, i.e.

${\mathcal{R}}(f)=\|f\|_{L^{1}(\mathscr{D})}+\frac{1}{2\beta}\|f\|_{L^{2}(\mathscr{D})}^{2},$

where $\beta>0$ is a large number. The method (1.5) then becomes

$\displaystyle\begin{split}f_{t}^{(n)}&=\beta\mbox{sign}(A^{*}\lambda_{t}^{(n)})\max\\{|A^{*}\lambda_{t}^{(n)}|-1\\},\\\
\lambda_{t+1}^{(n)}&=\lambda_{t}^{(n)}-\gamma(Af_{t}^{(n)}-\bar{u}^{(n)}),\end{split}$
(4.10)

where $\bar{u}^{(n)}$ is the average of a sequence of unbiased independent
identically distributed noisy data $\\{u_{i}\\}$, $i=1,\cdots,n$, of the exact
data $u$.

For numerical simulations we consider determining the initial data $f$ in the
time fractional diffusion equation

$\displaystyle\left\\{\begin{array}[]{lll}\partial_{t}^{\alpha}u(x,t)-\triangle
u(x,t)=0,&x\in\mathscr{D},\,t>0,\\\
u(x,t)=0,&x\in\partial\mathscr{D},\,t>0,\\\
u(x,0)=f(x),&x\in\mathscr{D}\end{array}\right.$

from the measurement of $u(x,T)$ at a fixed later time $T>0$, where
$\mathscr{D}=[0,1]\times[0,1]$, $0<\alpha<1$, and $\partial_{t}^{\alpha}u$
denotes the Caputo fractional derivative

$\partial_{t}^{\alpha}u(x,t)=\frac{1}{\Gamma(1-\alpha)}\int_{0}^{t}\frac{1}{(t-s)^{\alpha}}\frac{\partial
u}{\partial s}(x,s)dx$

with $\Gamma$ denoting the Gamma function. The mapping from $f$ to $u(x,T)$ is
a compact linear operator $A:L^{2}(\mathscr{D})\to L^{2}(\mathscr{D})$. Time
fractional diffusion equations occur naturally in anomalous diffusion in which
the variance of the process behaves like a non-integer power of time ([31]),
sharply contrast to the classical normal diffusion which is governed by the
heat equation.

Figure 5. The exact solution, the exact data and the noisy data (random
sample) for Example 4.3.

To solve this inverse problem numerically, we discretize $\mathscr{D}$ by
taking $(N+1)\times(N+1)$ grid points $(x_{i},y_{j}):=(i/N,j/N)$,
$i,j=0,1,\cdots,N$ and write $u_{i,j}(t)$ for $u(x_{i},y_{j},t)$ and $f_{i,j}$
for $f(x_{i},y_{j})$. Let $h=1/N$. Then, by the finite difference
approximation of $-\triangle u$, the diffusion equation becomes

$\displaystyle\begin{split}&\partial_{t}^{\alpha}u_{i,j}+\frac{1}{h^{2}}\left(4u_{i,j}-u_{i+1,j}-u_{i-1,j}-u_{i,j+1}-u_{i,j-1}\right)=0,\\\
&\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad i,j=1,\cdots,N-1,\\\
&u_{0,j}=u_{N,j}=u_{i,0}=u_{i,N}=0,\,\,\,\,\qquad i,j=0,\cdots,N,\\\
&u_{i,j}(0)=f_{i,j},\qquad i,j=1,\cdots,N-1.\end{split}$ (4.11)

It turns out that the solution of (4.11) has the form

$u_{i,j}(t)=\sum_{p,q=1}^{N-1}c_{p,q}(t)\sin(iph\pi)\sin(jqh\pi),$

where each $c_{p,q}(t)$ satisfies the fractional ordinary differential
equation

$\partial_{t}^{\alpha}c_{p,q}+\mu_{p,q}c_{p,q}=0,\quad\mbox{with
}\mu_{p,q}=\frac{1}{h^{2}}\left(4-2\cos(ph\pi)-2\cos(qh\pi)\right).$

Let $S$ and $S^{-1}$ denote the discrete sine transform and the inverse
disrete sine transform defined respectively by

$\displaystyle(Sv)_{p,q}$
$\displaystyle:=4h^{2}\sum_{i,j=1}^{N-1}v_{i,j}\sin(iph\pi)\sin(jqh\pi),$
$\displaystyle(S^{-1}v)_{i,j}$
$\displaystyle:=\sum_{p,q=1}^{N-1}v_{p,q}\sin(iph\pi)\sin(jqh\pi)$

for any $(N-1)\times(N-1)$ matrix $(v_{i,j})$. Then we have
$c_{p,q}(0)=(Sf)_{p,q}$ and thus
$c_{p,q}(t)=(Sf)_{p,q}E_{\alpha}(-\mu_{p,q}t^{\alpha})$, where $E_{\alpha}$
denotes the Mittag-Leffler function

$E_{\alpha}(z)=\sum_{k=0}^{\infty}\frac{z^{k}}{\Gamma(\alpha k+1)}.$

Consequently

$u_{i,j}(t)=\sum_{p,q=1}^{N-1}(Sf)_{p,q}E_{\alpha}(-\mu_{p,q}t^{\alpha})\sin(iph\pi)\sin(jqh\pi)$

for $i,j=1,\cdots,N-1$. If we define the transform $W_{t}$ by
$(W_{t}v)_{p,q}=E_{\alpha}(-\mu_{p,q}t^{\alpha})v_{p,q}$, then
$u(t)=S^{-1}W_{t}Sf$. Let $A=S^{-1}W_{T}S$. Then $f$ can be determined by
solving $Af=u$, where $u:=(u_{i,j}(T))$.

Figure 6. The reconstruction results of the method (4.10) with different
sample sizes for Example 4.3.

In our numerical simulation, we assume the sought solution is sparse and
consider $\mathscr{D}=[0,1]\times[0,1]$ on an equidistant grid of $256\times
256$ points. The sought solution $f^{\dagger}$ and the corresponding exact
data $u$ are plotted in Figure 5 (the left and middle ones). To reconstruct
$f^{\dagger}$, we use multiple repeated measurement data of different sample
sizes which are generated from $u$ by adding independent Gaussian noise
$N(0,\sigma^{2})$ with $\sigma=0.2\max|u|$; one sample of measurement data is
plotted in Figure 5 (the right one).

Table 3. Numerical results for Example 4.3, where the sought solution is sparse. n | Iteration numbers | Relative error | Emergency stops
---|---|---|---
| (in average) |  |
$10$ | 174 | 5.9865e-01 | 0
$10^{2}$ | 458 | 3.8077e-01 | 0
$10^{3}$ | 933 | 1.7786e-01 | 0
$10^{4}$ | 1821 | 1.1831e-01 | 0

When applying the method (4.10), we take $\beta=300$ and use the step-size
$\gamma=2/(\beta\|A\|^{2})$. The iteration is terminated by Rule 3.8 with
$\beta_{0}=50$ and $\tau_{n}$ given by (3.14); where $\tau_{0}=1.1$. The
reconstruction results for several different sample sizes
$n=10,10^{2},10^{3},10^{4}$ are plotted in Figure 6, each is based on the
average of 200 runs. In Table 2 we report the average number of iterations,
the average relative errors, and the number of emergency stops. Furthermore,
we present the boxplots of the relative errors in Figure 7. All these result
demonstrate that the proposed method can give satisfactory reconstruction
results when the sample size increases and the sparsity of the sought solution
can be captured well.

Figure 7. Boxplots of the relative errors for 200 simulations with different
sample sizes for Example 4.3.

######  Example 4.4.

In this final example we consider using the method (1.5) to reconstruct
piecewise constant solutions. We consider again the equation (4.4) with the
kernel given by (4.7) and assume that the sought solution $x^{\dagger}$ is
piecewise constant. By dividing $[0,1]$ into $m=400$ subintervals of equal
length and approximating integrals by the trapezoidal rule, we have a discrete
ill-posed problem ${\bf A}{\bf x}={\bf y}$, where ${\bf A}$ is a matrix. We
use the model

$\displaystyle\min\left\\{{\mathcal{R}}({\bf x}):=\|{\bf D}{\bf
x}\|_{1}+\frac{1}{2\beta}\|{\bf D}{\bf x}\|_{2}^{2}+\frac{1}{2\beta}\|{\bf
x}\|_{2}^{2}:{\bf A}{\bf x}={\bf y}\right\\},$ (4.12)

where ${\bf D}$ denotes the discrete gradient operator and $\beta$ is a large
positive number. Thus, ${\mathcal{R}}$ is a strongly convex perturbation of
the total variation $\|{\bf D}{\bf x}\|_{1}$. If we apply the method (1.5) to
(4.12) directly, we need to solve a minimization problem related to
${\mathcal{R}}$ to obtain ${\bf x}_{t}^{(n)}$ at each iteration. This can make
the algorithm time-consuming since those minimization problems can not be
solved explicitly.

Figure 8. The noisy data (10 samples, green circle), the exact data (blue one)
and the reconstruction results with different sample sizes for Example 4.4.

To circumvent this difficulty, by introducing ${\bf z}={\bf D}{\bf x}$ we
reformulate (4.12) as

$\displaystyle\min\left\\{f({\bf z})+\phi({\bf x}):{\bf A}{\bf x}={\bf
y}\mbox{ and }{\bf D}{\bf x}-{\bf z}=0\right\\},$ (4.13)

where

$f({\bf z})=\|{\bf z}\|_{1}+\frac{1}{2\beta}\|{\bf
z}\|_{2}^{2},\qquad\phi({\bf x})=\frac{1}{2\beta}\|{\bf x}\|_{2}^{2}.$

Since $f({\bf z})+\phi({\bf x})$ is strongly convex, we may apply the method
(1.5) to (4.13) to obtain the iteration scheme

$\displaystyle({\bf x}_{t}^{(n)},{\bf z}_{t}^{(n)})$
$\displaystyle=\arg\min_{{\bf x},{\bf z}}\left\\{f({\bf z})+\phi({\bf
x})-\langle\lambda_{t}^{(n)},{\bf A}{\bf x}\rangle-\langle\mu_{t}^{(n)},{\bf
D}{\bf x}-{\bf z}\rangle\right\\},$
$\displaystyle(\lambda_{t+1}^{(n)},\mu_{t+1}^{(n)})$
$\displaystyle=(\lambda_{t}^{(n)},\mu_{t}^{(n)})-\gamma({\bf A}{\bf
x}_{t}^{(n)}-\bar{\bf y}^{(n)},{\bf D}{\bf x}_{t}^{(n)}-{\bf z}_{t}^{(n)}).$

Since ${\bf x}_{t}^{(n)}$ and ${\bf z}_{t}^{(n)}$ can be given explicitly,
this leads to the following algorithm

$\displaystyle\begin{split}{\bf x}_{t}^{(n)}&=\beta({\bf
A}^{T}\lambda_{t}^{(n)}+{\bf D}^{T}\mu_{t}^{(n)}),\\\ {\bf
z}_{t}^{(n)}&=-\beta\mbox{sign}(\mu_{t}^{(n)})\max\\{|\mu_{t}^{(n)}|-1,0\\},\\\
\lambda_{t+1}^{(n)}&=\lambda_{t}^{(n)}-\gamma({\bf A}{\bf
x}_{t}^{(n)}-\bar{\bf y}^{(n)}),\\\ \mu_{t+1}^{(n)}&=\mu_{t}^{(n)}-\gamma({\bf
D}{\bf x}_{t}^{(n)}-{\bf z}_{t}^{(n)}).\end{split}$ (4.14)

Table 4. Numerical results for Example 4.4, where the sought solution is piecewise constant. n | Iteration numbers | Relative error | Emergency stops
---|---|---|---
| (in average) |  |
$10$ | 1109 | 3.5812e-01 | 10
$10^{2}$ | 17877 | 2.2729e-01 | 70
$10^{3}$ | 64078 | 1.7208e-01 | 0
$10^{4}$ | 242623 | 1.3263e-01 | 0
$10^{5}$ | 675068 | 1.0761e-01 | 0

In our numerical experiments, the sought solution is piecewise constant, whose
graph together with the exact data is plotted in Figure 8 in blue. In order to
apply the method (4.14), we use multiple repeated measurement data of several
different sample sizes $n=10,10^{2},10^{3},10^{4},10^{5}$ corrupted by
independent Gaussian noise $N(0,\sigma^{2})$ with $\sigma=0.2$ and let
$\bar{\bf y}^{(n)}$ be the average of these data; 10 samples of these
measurement data are plotted in Figure 8 in green. The step-size in (4.14) is
chosen to be $\gamma=2/[\beta(\|{\bf A}\|^{2}+4)]$. In order to enhance the
effect of total variation in reconstruction, we take $\beta=100$. The method
is terminated by Rule 3.8 adapted to (4.13) with $\beta_{0}=200$ and
$\tau_{n}$ given by (3.14), where $\tau_{0}=1.1$. For each sample size, we run
200 simulations and the reconstruction results in average are plotted in
Figure 8 in red. In Table 4 we report further numerical results including the
average number of iterations, the relative error in average, and the number of
emergency stops used. The boxplot of the relative error is presented in Figure
9. These results clearly demonstrate that the proposed method converges and
more and more satisfactory reconstruction results can be obtained when the
sample size increases.

Figure 9. Boxplots of the relative errors for 200 simulations with different
sample sizes for Example 4.4.

## References

  * [1] A. B. Bakushinskii, Remarks on choosing a regularization parameter using the quasioptimality and ratio criterion, Comp. Math. Math. Phys., 24 (1984), 181–182.
  * [2] A. Beck and M. Teboulle, A fast dual proximal gradient algorithm for convex minimization and applications, Oper. Res. Lett., 42 (2014), no. 1, 1–6.
  * [3] J. M. Borwein and C. S. Lewis, Convergence of best entropy estimates, SIAM J. Optim., 1 (1991), pp. 191–205.
  * [4] R. Boţ and T. Hein, Iterative regularization with a geeral penalty term – theory and applications to $L^{1}$ and TV regularization, Inverse Problems, 28 (2012), 104010.
  * [5] M. Burger, E. Resmerita and M. Benning, An entropic Landweber method for linear ill-posed problems, Inverse Problems, 36 (2019), 015009.
  * [6] D. H. Chen, D. Jiang, I. Yousept and J. Zou, Variational source conditions for inverse Robin and flux problems by partial measurements, Inverse Probl. Imaging, 16 (2022), no. 2, 283–304.
  * [7] D. H. Chen and I. Yousept, Variational source condition for ill-posed backward nonlinear Maxwell’s equations, Inverse Problem, 35 (2019), 025001.
  * [8] H. W. Engl, M. Hanke and A. Neubauer, Regularization of Inverse Problems, Kluwer, Dordrecht, 1996.
  * [9] J. Flemming,Existence of variational source conditions for nonlinear inverse problems in Banach spaces, J. Inverse Ill-posed Problems, 26 (2018), no. 2, 277–286.
  * [10] K. Frick and O. Scherzer,Regularization of ill-posed linear equations by the non-stationary augmented Lagrangian method, J. Integral Equ. Appl., 22 (2010), no. 2, 217–257.
  * [11] M. Grasmair, Source conditions for non-quadratic Tikhonov regularization, Numer. Funct. Anal. Optim., 41 (2020), no. 11, 1352–1372.
  * [12] P. C. Hansen and D. P. O. Leary, The use of the L-curve in the regularization of discrete ill-posed problems, SIAM J. Sci. Comput., 14 1487–503.
  * [13] B. Harrach, T. Jahn and R. Potthast, Beyond the Bakushinkii veto: regularising linear inverse problems without knowing the noise distribution, Numer. Math., 145 (2020), no. 3, 581–603.
  * [14] U. Hassan and M. S. Anwar, Reducing noise by repetition: introduction to signal averaging. Eur. J. Phys., 31 (2010), no. 3, 453–465.
  * [15] 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 (2007), 987–1010.
  * [16] B. Hofmann and P. Mathé, Parameter choice under variational inequalities, Inverse Problems, 28 (2012), 104006.
  * [17] 24\. T. Hohage and F. Weidling, Verification of a variational source condition for acoustic inverse medium scattering problems, Inverse Problems 31 (2015), no. 7, 075006.
  * [18] T. Hohage and F. Weilding, Variational source condition and stability estimates for inverse electromagnetic medium scattering problems, Inverse Probl. Imaging, 11 (2017), 203–220.
  * [19] T. Hohage and F. Weidling, Characterizations of variational source conditions, converse results, and maxisets of spectral regularization methods, SIAM J. Numer. Anal., 55 (2017), no. 2, 598–620.
  * [20] Q. Jin, Hanke-Raus heuristic rule for variational regularization in Banach spaces, Inverse Problems, 32 (2016), no. 8, 085008.
  * [21] Q. Jin, On a heuristic stopping rule for the regularization of inverse problems by the augmented Lagrangian method, Numer. Math., 136 (2017), no. 4, 973–992.
  * [22] Q. Jin, Convergence rates of a dual gradient method for constrained linear ill-posed problems, Numer. Math., 151 (2022), no. 4, 841–871.
  * [23] Q. Jin and L. Stals, Nonstationary iterated Tikhonov regularization for ill-posed problems in Banach spaces, Inverse Problems, 28 (2012), no. 10, 104011.
  * [24] Q. Jin and W. Wang, Landweber iteration of Kaczmarz type with general non-smooth convex penalty functionals, Inverse Problems, 29 (2013), no. 8, 085011, 22 pp.
  * [25] S. Kindermann and A. Neubauer, On the convergence of the quasioptimality criterion for (iterated) Tikhonov regularization, Inverse Probl. Imaging, 2 (2008), 291–299.
  * [26] S. Kindermann and K. Raik, Convergence of heuristic parameter choice rules for convex Tikhonov regularization, SIAM J. Numer. Anal., 58 (2020), no. 3, 1773–1800.
  * [27] J. Kivinen and M. K. Warmuth, Exponentiated gradient versus gradient descent for linear predictors, Information and Computation, 132 (1997), no. 1, 1–63.
  * [28] M. Ledoux and M. Talagrand, Probability in Banach Spaces: Isoperimetry and Processes, vol. 23. Springer, Berlin, 1991.
  * [29] R. G. Lyons, Understanding Digital Signal Processing, 3/E. Pearson Education India, Chennai, 2004.
  * [30] S. Nemirovski and D. B. Yudin, Problem Complexity and Method Efficiency in Optimization, Wiley, New York, 1983.
  * [31] I. Podlubny, Fractional Differential Equations, San Diego, CA: Academic, 1999.
  * [32] T. Schuster, B. Kaltenbacher, B. Hofmann and K. S. Kazimierski,Regularization Methods in Banach Spaces, Radon Series on Computational and Applied Mathematics 10 Walter de Gruyter, Berlin 2012.
  * [33] P. Tseng, Applications of a splitting algorithm to decomposition in convex programming and variational inequalities, SIAM J. Control Optim., 29 (1991), no. 1, 119–138.
  * [34] H. Uzawa, Iterative methods for concave programming, in: Studies in Linear and Nonlinear Programming, 1958, pp. 154–165.
  * [35] W. van Drongelen, Signal Processing for Neuroscientists, 2nd edition, Academic Press, 2018.
  * [36] C. Zălinscu, Convex Analysis in General Vector Spaces, World Scientific Publishing Co., Inc., River Edge, New Jersey, 2002.