Stochastic mirror descent method for linear ill-posed problems in Banach spaces

Qinian Jin Mathematical Sciences Institute, Australian National University, Canberra, ACT 2601, Australia qinian.jin@anu.edu.au , Xiliang Lu School of Mathematics and Statistics, Wuhan University, Wuhan 430072, China & Hubei Computational Science Key Laboratory, Wuhan University, Wuhan 430072, China xllv.math@whu.edu.cn and Liuying Zhang School of Mathematics and Statistics, Wuhan University, Wuhan 430072, China lyzhang.math@whu.edu.cn

Abstract.

Consider linear ill-posed problems governed by the system $A_{i}x=y_{i}$ for $i=1,\cdots,p$ , where each $A_{i}$ is a bounded linear operator from a Banach space $X$ to a Hilbert space $Y_{i}$ . In case $p$ is huge, solving the problem by an iterative regularization method using the whole information at each iteration step can be very expensive, due to the huge amount of memory and excessive computational load per iteration. To solve such large-scale ill-posed systems efficiently, we develop in this paper a stochastic mirror descent method which uses only a small portion of equations randomly selected at each iteration steps and incorporates convex regularization terms into the algorithm design. Therefore, our method scales very well with the problem size and has the capability of capturing features of sought solutions. The convergence property of the method depends crucially on the choice of step-sizes. We consider various rules for choosing step-sizes and obtain convergence results under a priori early stopping rules. In particular, by incorporating the spirit of the discrepancy principle we propose a choice rule of step-sizes which can efficiently suppress the oscillations in iterates and reduce the effect of semi-convergence. Furthermore, we establish an order optimal convergence rate result when the sought solution satisfies a benchmark source condition. Various numerical simulations are reported to test the performance of the method.

1. Introduction

Due to rapid growth of data sizes in practical applications, in recent years stochastic optimization methods have received tremendous attention and have been proved to be efficient in various applications of science and technology including in particular the machine learning researches ([5, 13]). In this paper we will develop a stochastic mirror descent method for solving linear ill-posed inverse problems in Banach spaces.

We will consider linear ill-posed inverse problems governed by the system

\displaystyle A_{i}x=y_{i},\qquad i=1,\cdots,p

(1.1)

consisting of $p$ linear equations, where, for each $i=1,\cdots,p$ , $A_{i}:X\to Y_{i}$ is a bounded linear operator from a Banach space $X$ to a Hilbert space $Y_{i}$ . Such systems arise naturally in many practical applications. For instance, many linear ill-posed inverse problems can be described by an integral equation of the first kind ([10, 14])

(Ax)(s):=\int_{\mathscr{D}}k(s,t)x(t)dt=y(s),\quad s\in\mathscr{D}^{\prime},

where $\mathscr{D}$ and $\mathscr{D}^{\prime}$ are bounded domains in Euclidean spaces and the kernel $k$ is a bounded continuous function on $\mathscr{D}^{\prime}\times\mathscr{D}$ . Clearly $A$ is a bounded linear operator from $L^{r}(\mathscr{D})$ to $C(\mathscr{D}^{\prime})$ for any $1\leq r\leq\infty$ . By taking $p$ sample points $s_{1},\cdots,s_{p}$ in $\mathscr{D}^{\prime}$ , then the problem of determining a solution using only the knowledge of $y_{i}:=y(s_{i})$ for $i=1,\cdots,p$ reduces to solving a linear system of the form (1.1), where $A_{i}:L^{r}(\mathscr{D})\to{\mathbb{R}}$ is given by

A_{i}x:=\int_{\mathscr{D}}k(s_{i},t)x(t)dt

for each $i$ . Further examples of (1.1) can be found in various tomographic techniques using multiple measurements [27].

Throughout the paper we always assume (1.1) has a solution. The system (1.1) may have many solutions. By taking into account of a priori information about the sought solution, we may use a proper, lower semi-continuous, convex function ${\mathcal{R}}:X\to(-\infty,\infty]$ to select a solution $x^{\dagger}$ of (1.1) such that

\displaystyle{\mathcal{R}}(x^{\dagger})=\min\left\{{\mathcal{R}}(x):A_{i}x=y_{i}\mbox{ for }i=1,\cdots,p\right\}

(1.2)

which, if exists, is called a ${\mathcal{R}}$ -minimizing solution of (1.1). In practical applications, the exact data $y:=(y_{1},\cdots,y_{p})$ is in general not available, instead we only have noisy data $y^{\delta}:=(y_{1}^{\delta},\cdots,y_{p}^{\delta})$ satisfying

\displaystyle\|y_{i}^{\delta}-y_{i}\|\leq\delta_{i},\quad i=1,\cdots,p,

(1.3)

where $\delta_{i}>0$ denotes the noise level corresponding to data in the space $Y_{i}$ . How to use the noisy data $y^{\delta}$ to construct an approximate solution of the ${\mathcal{R}}$ -minimizing solution of (1.1) is an important topic.

Let $Y:=Y_{1}\times\cdots\times Y_{p}$ and define $A:X\to Y$ by

Ax=(A_{1}x,\cdots,A_{p}x),\quad x\in X.

Then (1.2) can be equivalently stated as

\displaystyle{\mathcal{R}}(x^{\dagger})=\min\left\{{\mathcal{R}}(x):Ax=y\right\}.

(1.4)

which has been considered by various variational and iterative regularization methods, see [4, 7, 21, 23, 35] for instance. In particular, the Landweber iteration in Hilbert spaces has been extended for solving (1.4), leading to the iterative method of the form

\displaystyle\begin{split}&x_{n}^{\delta}=\arg\min_{x\in X}\left\{{\mathcal{R}}(x)-\langle\xi_{n}^{\delta},x\rangle\right\},\\ &\xi_{n+1}^{\delta}=\xi_{n}^{\delta}-t_{n}^{\delta}A^{*}(Ax_{n}^{\delta}-y^{\delta}),\end{split}

(1.5)

where $A^{*}:Y\to X^{*}$ denotes the adjoint of $A$ and $t_{n}^{\delta}$ is the step-size. This method can be derived as a special case of the mirror descent method; see Section 2 for a brief account. The method (1.5) has been investigated in a number of references, see [4, 12, 21, 22, 23]. In particular, the convergence and convergence rates have been derived in [21] when the method is terminated by either an a priori stopping rule or the discrepancy principle

\|Ax_{n_{\delta}}^{\delta}-y^{\delta}\|\leq\tau\delta<\|Ax_{n}^{\delta}-y^{\delta}\|,\quad 0\leq n<n_{\delta},

where

\delta:=\sqrt{\delta_{1}^{2}+\cdots+\delta_{p}^{2}}

denotes the total noise level of the noisy data. We remark that the minimization problem in (1.5) for defining $x_{n}^{\delta}$ can be solved easily in general as it does not depend on $A$ ; in fact $x_{n}^{\delta}$ can be given by an explicit formula in many interesting cases; even if $x_{n}^{\delta}$ does not have an explicit formula, there exist fast algorithms for determining $x_{n}^{\delta}$ efficiently. However, note that

A^{*}(Ax_{n}^{\delta}-y^{\delta})=\sum_{i=1}^{p}A_{i}^{*}(A_{i}x_{n}^{\delta}-y_{i}^{\delta}).

Therefore, updating $\xi_{n}^{\delta}$ to $\xi_{n+1}^{\delta}$ requires calculating $A_{i}^{*}(A_{i}x_{n}^{\delta}-y_{i}^{\delta})$ for all $i=1,\cdots,p$ . In case $p$ is huge, using the method (1.5) to solve (1.4) can be inefficient because it requires a huge amount of memory and excessive computational work per iteration.

In order to relieve the drawback of the method (1.5), by extending the Kaczmarz-type method [15] in Hilbert spaces, a Landweber-Kaczmarz method in Banach spaces has been proposed in [20, 23] to solve (1.2) which cyclically considers each equation in (1.2) in a Gauss-Seidel manner and the iteration scheme takes the form

\displaystyle\begin{split}&x_{n}^{\delta}=\arg\min_{x\in X}\left\{{\mathcal{R}}(x)-\langle\xi_{n}^{\delta},x\rangle\right\},\\ &\xi_{n+1}^{\delta}=\xi_{n}^{\delta}-t_{n}^{\delta}A_{i_{n}}^{*}(A_{i_{n}}x_{n}^{\delta}-y_{i_{n}}^{\delta}),\end{split}

(1.6)

where $i_{n}=(n\mod p)+1$ . The convergence of this method has been shown in [23] in which the numerical results demonstrate its nice performance. However, it should be pointed out that the efficiency of the method (1.6) depends crucially on the order of the equations and its convergence speed is difficult to be quantified. In order to resolve these issues, in this paper we will consider a stochastic version of (1.6), namely, instead of taking $i_{n}$ cyclically, we will choose $i_{n}$ from $\{1,\cdots,p\}$ randomly at each iteration step. The corresponding method will be called the stochastic mirror descent method and more details will be presented in Section 2 where we also propose the mini-batch version of the stochastic mirror descent method.

The stochastic mirror descent method, that we will consider in this paper, includes the stochastic gradient descent as a special case. Indeed, when $X$ is a Hilbert space and ${\mathcal{R}}(x)=\|x\|^{2}/2$ , the method (1.6) becomes

\displaystyle x_{n+1}^{\delta}=x_{n}^{\delta}-t_{n}^{\delta}A_{i_{n}}^{*}(A_{i_{n}}x_{n}^{\delta}-y_{i_{n}}^{\delta})

(1.7)

which is exactly the stochastic gradient descent method studied in [18, 19, 25, 33] for solving linear ill-posed problems in Hilbert spaces. 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. Unfortunately, the stochastic gradient descent method does not have the capability to incorporate these information into the algorithm. However, this can be handled by the stochastic mirror descent method with careful choices of a suitable Banach space $X$ and a strongly convex penalty functional ${\mathcal{R}}$ .

In this paper we will use tools from convex analysis in Banach spaces to analyze the stochastic mirror descent method. The choice of the step-size plays a crucial role on the convergence of the method. We consider several rules for choosing the step-sizes and provide criteria for terminating the iterations in order to guarantee a convergence when the noise level tends to zero. The iterates produced by the stochastic mirror descent method exhibit salient oscillations and, due to the ill-posedness of the underlying problems, the method using noisy data demonstrates the semi-convergence property, i.e. the iterate tends to the sought solution at the beginning and, after a critical number of iterations, the iterates diverges. The oscillations and semi-convergence make it difficult to determine an output with good approximation property, in particular when the noise level is relatively large. When the information on noise level is available, by incorporating the spirit of the discrepancy principle we propose a rule for choosing step-size. This rule enables us to efficiently suppress the oscillations of iterates and remove the semi-convergence of the method as indicated by the extensive numerical simulations. Furthermore, we obtain an order optimal convergence rate result for the stochastic mirror descent method with constant step-size when the sought solution satisfies a benchmark source condition. We achieve this by interpreting the stochastic mirror descent method equivalently as a randomized block gradient method applied to the dual problem of (1.1). Even for the stochastic gradient descent method, our convergence rate result supplements the existing results since only sub-optimal convergence rates have been derived under diminishing step-sizes, see [18, 25].

This paper is organized as follows. In Section 2 we first collect some basic facts on convex analysis in Banach spaces and then give an account on the stochastic mirror descent method. In Section 3 we prove some convergence results on the stochastic mirror descent method under various choices of the step-sizes. When the sought solution satisfies a benchmark source condition, in Section 4 we establish an order optimal convergence rate result. Finally, in Section 5 we present extensive numerical simulations to test the performance of the stochastic mirror descent method.

2. The method

2.1. Preliminaries

In this section, we will collect some basic facts on convex analysis in Banach spaces which will be used in the analysis of the stochastic mirror descent method; for more details one may refer to [34, 38] for instance.

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

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

If $\mbox{dom}(f)\neq\emptyset$ , $f$ is called proper. Given $x\in\mbox{dom}(f)$ , an element $\xi\in 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 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$ . Thus $x\to\partial f(x)$ defines a set-valued mapping $\partial f$ whose domain of definition is defined as

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

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 X

which is always nonnegative.

A proper function $f:X\to(-\infty,\infty]$ is called strongly convex if there exists a constant $\sigma>0$ such that

\displaystyle f(t\bar{x}+(1-t)x)+\sigma t(1-t)\|\bar{x}-x\|^{2}\leq tf(\bar{x})+(1-t)f(x)

(2.1)

for all $\bar{x},x\in\mbox{dom}(f)$ and $t\in[0,1]$ . The largest number $\sigma>0$ such that (2.1) holds true is called the modulus of convexity of $f$ . It is easy to see that for a proper, strongly convex function $f:X\to(-\infty,\infty]$ with modulus of convexity $\sigma>0$ there holds

\displaystyle D_{f}^{\xi}(\bar{x},x)\geq\sigma\|x-\bar{x}\|^{2}

(2.2)

for all $\bar{x}\in\mbox{dom}(f)$ , $x\in\mbox{dom}(\partial f)$ and $\xi\in\partial f(x)$ .

For a proper function $f:X\to(-\infty,\infty]$ , its Legendre-Fenchel conjugate is defined by

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

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

\displaystyle f^{*}(\xi)+f(x)\geq\langle\xi,x\rangle

(2.3)

for all $x\in X$ and $\xi\in X^{*}$ . 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.4)

The following important result gives further properties of $f^{*}$ which in particular shows that the strong convexity of $f$ implies the continuous differentiability of $f^{*}$ with gradient $\nabla f^{*}$ mapping $X^{*}$ to $X$ ; see [38, Corollary 3.5.11]

Proposition 2.1.

Let $X$ be a Banach space and let $f:X\to(-\infty,\infty]$ be a proper, lower semi-continuous, strongly convex function with modulus of convexity $\sigma>0$ . Then $\mbox{dom}(f^{*})=X^{*}$ , $f^{*}$ is Fréchet differentiable and its gradient $\nabla f^{*}$ maps $X^{*}$ into $X$ satisfying

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

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

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

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

(2.5)

which is involved in the formulation of the stochastic mirror descent method below. According to [38, Theorem 3.5.8], (2.4) and Proposition 2.1 we have

Proposition 2.2.

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

2.2. Description of the method

In order to motivate our method, we first briefly review how to extend the gradient method in Hilbert spaces to solve the minimization problem

\displaystyle\min_{x\in X}f(x)

(2.6)

in Banach spaces, where $f:X\to{\mathbb{R}}$ is a Fréchet differentiable function defined on a Banach space $X$ . When $X$ is a Hilbert space, the gradient method for solving (2.6) takes the form $x_{n+1}=x_{n}-t_{n}\nabla f(x_{n})$ which can be equivalently stated as

	$\displaystyle x_{n+1}$	$\displaystyle=\arg\min_{x\in X}\left\{\frac{1}{2}\\|x-(x_{n}-t_{n}\nabla f(x_{n}))\\|^{2}\right\}$
		$\displaystyle=\arg\min_{x\in X}\left\{\frac{1}{2}\\|x-x_{n}\\|^{2}+t_{n}\langle\nabla f(x_{n}),x\rangle\right\},$		(2.7)

where $\langle\cdot,\cdot\rangle$ denotes the inner product on $X$ and $\|\cdot\|$ denotes the induced norm. When $X$ is a general Banach space, inner product is no longer available and $\nabla f(x_{n})$ is not necessarily an element in $X$ ; instead $\nabla f(x_{n})$ is an element in $X^{*}$ . Therefore, in order to guarantee $\langle\nabla f(x_{n}),x\rangle$ to be meaningful, this expression should be understood as the dual pairing between $X^{*}$ and $X$ . On the other hand, since no Hilbert space norm is available, $\frac{1}{2}\|x-x_{n}\|^{2}$ in (2.2) should be replaced by other suitable distance-like functionals. In order to capture the feature of the sought solution, a suitable convex penalty function ${\mathcal{R}}:X\to(-\infty,\infty]$ is usually chosen to enhance the feature. In such a situation, one may use the Bregman distance induced by ${\mathcal{R}}$ to fulfill the purpose. To be more precise, if $\partial{\mathcal{R}}(x_{n})\neq\emptyset$ , we may take $\xi_{n}\in\partial{\mathcal{R}}(x_{n})$ suitably and then use the Bregman distance

D_{\mathcal{R}}^{\xi_{n}}(x,x_{n})={\mathcal{R}}(x)-{\mathcal{R}}(x_{n})-\langle\xi_{n},x-x_{n}\rangle

to replace $\frac{1}{2}\|x-x_{n}\|^{2}$ in (2.2). This leads to the new updating formula

x_{n+1}\in\arg\min_{x\in X}\left\{D_{\mathcal{R}}^{\xi_{n}}(x,x_{n})+t_{n}\langle\nabla f(x_{n}),x\rangle\right\}

for $x_{n+1}$ . By the expression of $D_{\mathcal{R}}^{\xi_{n}}(x,x_{n})$ , we have

x_{n+1}\in\arg\min_{x\in X}\left\{{\mathcal{R}}(x)-\langle\xi_{n}-t_{n}\nabla f(x_{n}),x\rangle\right\}.

Let $\xi_{n+1}=\xi_{n}-t_{n}\nabla f(x_{n})$ , then

x_{n+1}\in\arg\min_{x\in X}\left\{{\mathcal{R}}(x)-\langle\xi_{n+1},x\rangle\right\}.

If $x_{n+1}$ is well-defined, then by the optimality condition on $x_{n+1}$ we have $\xi_{n+1}\in\partial{\mathcal{R}}(x_{n+1})$ . Therefore, we can repeat the above procedure, leading to the algorithm

\displaystyle\begin{split}&x_{n}=\arg\min_{x\in X}\left\{{\mathcal{R}}(x)-\langle\xi_{n},x\rangle\right\},\\ &\xi_{n+1}=\xi_{n}-t_{n}\nabla f(x_{n})\end{split}

(2.8)

for solving (2.6) in Banach spaces, which is called the mirror descent method in optimization community; see [1, 6, 30, 31].

In order to apply the mirror descent method to solve (1.2) when only noisy data $y_{i}^{\delta}$ are available, we may consider a problem of the form (2.6) with

f(x)=\frac{1}{2}\sum_{i=1}^{p}\|A_{i}x-y_{i}^{\delta}\|^{2}.

Note that $\nabla f(x)=\sum_{i=1}^{p}A_{i}^{*}(A_{i}x-y_{i}^{\delta})$ . Therefore, an application of (2.8) gives

\displaystyle\begin{split}&x_{n}=\arg\min_{x\in X}\left\{{\mathcal{R}}(x)-\langle\xi_{n},x\rangle\right\},\\ &\xi_{n+1}=\xi_{n}-t_{n}\sum_{i=1}^{p}A_{i}^{*}(A_{i}x_{n}-y_{i}^{\delta}).\end{split}

(2.9)

This is exactly the method (1.5) which has been considered in [4, 21, 23]. Clearly, the implementation of (2.9) requires to calculate $A_{i}^{*}(A_{i}x_{n}-y_{i}^{\delta})$ for all $i=1,\cdots,p$ at each iteration. In case $p$ is huge, each iteration step in (2.9) requires a huge amount of computational work. In order to reduce the computational load per iteration, one may randomly choose a subset $I_{n}$ from $\{1,\cdots,p\}$ with small size to form the partial term

f_{I_{n}}(x):=\frac{1}{2}\sum_{i\in I_{n}}\|A_{i}x-y_{i}^{\delta}\|^{2}

of $f$ and use its gradient $\nabla f_{I_{n}}(x_{n})=\sum_{i\in I_{n}}A_{i}^{*}(A_{i}x_{n}-y_{i}^{\delta})$ at $x_{n}$ as a replacement of $\nabla f(x_{n})$ in (2.9). This leads to the following mini-batch stochastic mirror descent method for solving (1.2) with noisy data.

Algorithm 1.

Fix a batch size $b$ , pick the initial guess $\xi_{0}=0$ in $X^{*}$ and set $\xi_{0}^{\delta}:=\xi_{0}$ . For $n\geq 0$ do the following:

(i)

Calculate $x_{n}^{\delta}\in X$ by solving

$\displaystyle x_{n}^{\delta}=\arg\min_{x\in X}\left\{{\mathcal{R}}(x)-\langle\xi_{n}^{\delta},x\rangle\right\};$ (2.10)
(ii)

Randomly select a subset $I_{n}\subset\{1,\cdots,p\}$ with $|I_{n}|=b$ via the uniform distribution;
(iii)

Choose a step-size $t_{n}^{\delta}\geq 0$ ;

(iv)

Define $\xi_{n+1}^{\delta}\in X^{*}$ by

\displaystyle\xi_{n+1}^{\delta}=\xi_{n}^{\delta}-t_{n}^{\delta}\sum_{i\in I_{n}}A_{i}^{*}\left(A_{i}x_{n}^{\delta}-y_{i}^{\delta}\right).

(2.11)

It should be pointed out that the choice of the step-size $t_{n}^{\delta}$ plays a crucial role on the convergence of Algorithm 1. We will consider in Section 3 several rules for choosing the step-size. Note also that, at each iteration step of Algorithm 1, $x_{n}^{\delta}$ is defined by a minimization problem (2.10). When ${\mathcal{R}}$ is proper, lower semi-continuous and strongly convex, it follows from Proposition 2.2 that $x_{n}^{\delta}$ is well-defined. Moreover, for many important choices of ${\mathcal{R}}$ , $x_{n}^{\delta}$ can be given by an explicit formula, see Section 5 for instance, and thus the calculation of $x_{n}^{\delta}$ does not take much time.

There exist extensive studies on the mirror descent method and its stochastic variants in optimization, see [6, 8, 28, 29, 39] for instance. The existing works either depend crucially on the finite-dimensionality of the underlying spaces or establish only error estimates in terms of objective function values, and therefore they are not applicable to our Algorithm 1 for ill-posed problems. We need to develop new analysis. Our analysis of Algorithm 1 is based on the following assumption which is assumed throughout the paper.

Assumption 1.

(i)

$X$ is a Banach space, $Y_{i}$ is a Hilbert space and $A_{i}:X\to Y_{i}$ is a bounded linear operator for each $i=1,\cdots,p$ ;
(ii)

${\mathcal{R}}:X\to(-\infty,\infty]$ is a proper, lower semi-continuous, strongly convex function with modulus of convexity $\sigma>0$ ;
(iii)

The system $A_{i}x=y_{i}$ , $i=1,\cdots,p$ , has a solution in $\emph{dom}({\mathcal{R}})$ .

According to Assumption 1 (ii) and Proposition 2.1, the Legendre-Fenchel conjugate ${\mathcal{R}}^{*}$ of ${\mathcal{R}}$ is continuous differentiable with Lipschitz continuous gradient, i.e.

\displaystyle\|\nabla{\mathcal{R}}^{*}(\xi)-\nabla{\mathcal{R}}^{*}(\eta)\|\leq\frac{\|\xi-\eta\|}{2\sigma},\quad\forall\xi,\eta\in X^{*}.

(2.12)

Moreover, by virtue of Assumption 1 and Proposition 2.2, the linear system (1.1) has a unique ${\mathcal{R}}$ -minimizing solution $x^{\dagger}$ and Algorithm 1 is well-defined. By the optimality condition on $x_{n}^{\delta}$ and Proposition 2.2 we have

\displaystyle\xi_{n}^{\delta}\in\partial{\mathcal{R}}(x_{n}^{\delta})\quad\mbox{ and }\quad x_{n}^{\delta}=\nabla{\mathcal{R}}^{*}(\xi_{n}^{\delta})

(2.13)

which is the starting point of our convergence analysis.

3. Convergence

In order to establish a convergence result on Algorithm 1, we need to specify a probability space on which the analysis will be carried out. Let $\Lambda_{p}:=\{1,\cdots,p\}$ and let ${\mathcal{S}}_{p}$ denote the $\sigma$ -algebra consisting of all subsets of $\Lambda_{p}$ . Recall that, at each iteration step, a subset of indices is randomly chosen from $\Lambda_{p}$ via the uniform distribution. Therefore, for each $n\geq 1$ , it is natural to consider $x_{n}^{\delta}$ and $\xi_{n}^{\delta}$ on the sample space

\Lambda_{p}^{n}:=\underbrace{\Lambda_{p}\times\cdots\times\Lambda_{p}}_{n\mbox{ copies}}

equipped with the $\sigma$ -algebra ${\mathcal{S}}_{p}^{\otimes n}$ and the uniform distributed probability, denoted as ${\mathbb{P}}_{n}$ . According to the Kolmogorov extension theorem ([2]), there exists a unique probability ${\mathbb{P}}$ defined on the measurable space $(\Omega,{\mathcal{F}}):=(\Lambda_{p}^{\infty},{\mathcal{S}}_{p}^{\otimes\infty})$ such that each ${\mathbb{P}}_{n}$ is consistent with ${\mathbb{P}}$ . Let ${\mathbb{E}}$ denote the expectation on the probability space $(\Omega,{\mathcal{F}},{\mathbb{P}})$ . Given a Banach space $X$ , we use $L^{2}(\Omega,X)$ to denote the space consisting of all random variables $x$ with values in $X$ such that $\mathbb{E}[\|x\|^{2}]$ is finite; this is a Banach space under the norm

\left(\mathbb{E}[\|x\|^{2}]\right)^{1/2}.

Concerning Algorithm 1 we will use $\{{\mathcal{F}}_{n}\}$ to denote the natural filtration, where ${\mathcal{F}}_{n}:=\sigma(I_{0},\cdots,I_{n-1})$ for each $n\geq 1$ . We will frequently use the identity

\displaystyle\mathbb{E}[\phi]=\mathbb{E}[\mathbb{E}[\phi|{\mathcal{F}}_{n}]]

(3.1)

for any random variable $\phi$ on $(\Omega,{\mathcal{F}},{\mathbb{P}})$ . where $\mathbb{E}[\phi|{\mathcal{F}}_{n}]$ denotes the expectation of $\phi$ conditioned on ${\mathcal{F}}_{n}$ .

In this section we will prove some convergence results on Algorithm 1 under suitable choices of the step-sizes. The convergence of Algorithm 1 will be established by investigating the convergence property of its counterpart for exact data together with its stability property. For simplicity of exposition, for each index set $I=\{i_{1},\cdots,i_{b}\}\subset\{1,\cdots,p\}$ of size $b$ we set

Y_{I}:=Y_{i_{1}}\times\cdots\times Y_{i_{b}},\quad y_{I}^{\delta}:=(y_{i_{1}}^{\delta},\cdots,y_{i_{b}}^{\delta}),\quad\delta_{I}^{2}=\delta_{i_{1}}^{2}+\cdots+\delta_{i_{b}}^{2}

and define $A_{I}:X\to Y_{I}$ by

A_{I}x:=(A_{i_{1}}x,\cdots,A_{i_{b}}x),\quad\forall x\in X.

Let $A_{I}^{*}$ denote the adjoint of $A_{I}$ . Then the updating formula (2.11) of $\xi_{n+1}^{\delta}$ from $\xi_{n}^{\delta}$ can be rephrased as

\xi_{n+1}^{\delta}=\xi_{n}^{\delta}-t_{n}^{\delta}A_{I_{n}}^{*}(A_{I_{n}}x_{n}^{\delta}-y_{I_{n}}^{\delta}).

We start with the following result.

Lemma 3.1.

Let Assumption 1 hold. Consider Algorithm 1 and assume that

0\leq t_{n}^{\delta}\leq\min\left\{\frac{\mu_{0}\|A_{I_{n}}x_{n}^{\delta}-y_{I_{n}}^{\delta}\|^{2}}{\|A_{I_{n}}^{*}(A_{I_{n}}x_{n}^{\delta}-y_{I_{n}}^{\delta})\|^{2}},\mu_{1}\right\}

for all $n\geq 0$ , where $\mu_{0}$ and $\mu_{1}$ are two positive constants with $c_{0}:=1-\mu_{0}/(4\sigma)>0$ . Let $\hat{x}$ be any solution of (1.1) in $\mbox{dom}({\mathcal{R}})$ and let

\Delta_{n}^{\delta}:=D_{\mathcal{R}}^{\xi_{n}^{\delta}}(\hat{x},x_{n}^{\delta}),\quad n=0,1,\cdots.

Then

\displaystyle\mathbb{E}\left[\Delta_{n+1}^{\delta}\right]-\mathbb{E}\left[\Delta_{n}^{\delta}\right]\leq\frac{\mu_{1}b}{4c_{0}p}\delta^{2}.

Proof.

Note that

\Delta_{n+1}^{\delta}-\Delta_{n}^{\delta}=\left(\langle\xi_{n+1}^{\delta},x_{n+1}^{\delta}-\hat{x}\rangle-{\mathcal{R}}(x_{n+1}^{\delta})\right)+\left({\mathcal{R}}(x_{n}^{\delta})+\langle\xi_{n}^{\delta},\hat{x}-x_{n}^{\delta}\rangle\right).

By using (2.13) and (2.4) we have

{\mathcal{R}}(x_{n}^{\delta})+{\mathcal{R}}^{*}(\xi_{n}^{\delta})=\langle\xi_{n}^{\delta},x_{n}^{\delta}\rangle

for all $n\geq 0$ . Therefore

\displaystyle\Delta_{n+1}^{\delta}-\Delta_{n}^{\delta}

\displaystyle=\left({\mathcal{R}}^{*}(\xi_{n+1}^{\delta})-\langle\xi_{n+1}^{\delta},\hat{x}\rangle\right)-\left({\mathcal{R}}^{*}(\xi_{n}^{\delta})-\langle\xi_{n}^{\delta},\hat{x}\rangle\right).

Using $x_{n}^{\delta}=\nabla{\mathcal{R}}^{*}(\xi_{n}^{\delta})$ in (2.13) and the inequality (2.12) on $\nabla{\mathcal{R}}^{*}$ , we can obtain

	$\displaystyle\Delta_{n+1}^{\delta}-\Delta_{n}^{\delta}$	$\displaystyle=\left({\mathcal{R}}^{}(\xi_{n+1}^{\delta})-{\mathcal{R}}^{}(\xi_{n}^{\delta})-\langle\xi_{n+1}^{\delta}-\xi_{n}^{\delta},\nabla{\mathcal{R}}^{*}(\xi_{n}^{\delta})\rangle\right)$
		$\displaystyle\quad\,+\langle\xi_{n+1}^{\delta}-\xi_{n}^{\delta},x_{n}^{\delta}-\hat{x}\rangle$
		$\displaystyle\leq\frac{1}{4\sigma}\\|\xi_{n+1}^{\delta}-\xi_{n}^{\delta}\\|^{2}+\langle\xi_{n+1}^{\delta}-\xi_{n}^{\delta},x_{n}^{\delta}-\hat{x}\rangle.$

According to the definition of $\xi_{n+1}^{\delta}$ and $A_{I_{n}}\hat{x}=y_{I_{n}}$ , we can further obtain

\displaystyle\Delta_{n+1}^{\delta}-\Delta_{n}^{\delta}

\displaystyle\leq\frac{1}{4\sigma}\left(t_{n}^{\delta}\right)^{2}\|A_{I_{n}}^{*}(A_{I_{n}}x_{n}^{\delta}-y_{I_{n}}^{\delta})\|^{2}-t_{n}^{\delta}\langle A_{I_{n}}x_{n}^{\delta}-y_{I_{n}}^{\delta},A_{I_{n}}x_{n}^{\delta}-y_{I_{n}}\rangle.

By the given condition on $t_{n}^{\delta}$ and the Cauchy-Schwarz inequality, we then obtain

	$\displaystyle\Delta_{n+1}^{\delta}-\Delta_{n}^{\delta}$	$\displaystyle\leq\frac{\mu_{0}}{4\sigma}t_{n}^{\delta}\\|A_{I_{n}}x_{n}^{\delta}-y_{I_{n}}^{\delta}\\|^{2}-t_{n}^{\delta}\langle A_{I_{n}}x_{n}^{\delta}-y_{I_{n}}^{\delta},A_{I_{n}}x_{n}^{\delta}-y_{I_{n}}\rangle$
		$\displaystyle\leq-c_{0}t_{n}^{\delta}\\|A_{I_{n}}x_{n}^{\delta}-y_{I_{n}}^{\delta}\\|^{2}+t_{n}^{\delta}\delta_{I_{n}}\\|A_{I_{n}}x_{n}^{\delta}-y_{I_{n}}^{\delta}\\|$
		$\displaystyle\leq\frac{1}{4c_{0}}t_{n}^{\delta}\delta_{I_{n}}^{2}\leq\frac{\mu_{1}}{4c_{0}}\delta_{I_{n}}^{2}.$

By taking the expectation and using $\sum_{I:|I|=b}$ to denote a sum over all subsets $I\subset\{1,\cdots,n\}$ with $|I|=b$ , we have

	$\displaystyle\mathbb{E}[\Delta_{n+1}^{\delta}]-\mathbb{E}[\Delta_{n}^{\delta}]$	$\displaystyle\leq\frac{\mu_{1}}{4c_{0}}\mathbb{E}[\delta_{I_{n}}^{2}]=\frac{\mu_{1}}{4c_{0}}\frac{1}{{p\choose b}}\sum_{I:\|I\|=b}\delta_{I}^{2}=\frac{\mu_{1}}{4c_{0}}\frac{1}{{p\choose b}}\sum_{I:\|I\|=b}\sum_{i\in I}\delta_{i}^{2}$
		$\displaystyle=\frac{\mu_{1}}{4c_{0}}\frac{1}{{p\choose b}}\sum_{i=1}^{p}\sum_{I:\|I\|=b\,\&\,i\in I}\delta_{i}^{2}=\frac{\mu_{1}}{4c_{0}}\frac{{p-1\choose b-1}}{{p\choose b}}\sum_{i=1}^{p}\delta_{i}^{2}$
		$\displaystyle=\frac{\mu_{1}b}{4c_{0}p}\delta^{2}$

which shows the desired inequality. ∎

Next we consider Algorithm 1 with exact data and drop the superscript $\delta$ for every quantity defined by the algorithm, Thus $x_{n},\xi_{n}$ denote the corresponding iterative sequences and $t_{n}$ denotes the step-size. We now show a convergence result for Algorithm 1 with exact data by demonstrating that $\{x_{n}\}$ is a Cauchy sequence in $L^{2}(\Omega,X)$ .

Theorem 3.2.

Let Assumption 1 hold. Consider Algorithm 1 with exact data and assume that

\displaystyle\mu_{2}\leq t_{n}\leq\min\left\{\frac{\mu_{0}\|A_{I_{n}}x_{n}-y_{I_{n}}\|^{2}}{\|A_{I_{n}}^{*}(A_{I_{n}}x_{n}-y_{I_{n}})\|^{2}},\mu_{1}\right\}\quad\mbox{ when }A_{I_{n}}x_{n}\neq y_{I_{n}},

(3.2)

where $\mu_{0}$ , $\mu_{1}$ and $\mu_{2}$ are positive numbers with $c_{0}:=1-\mu_{0}/(4\sigma)>0$ . Then

\mathbb{E}[\|x_{n}-x^{\dagger}\|^{2}]\to 0\quad\mbox{ and }\quad\mathbb{E}\left[D_{\mathcal{R}}^{\xi_{n}}(x^{\dagger},x_{n})\right]\to 0

as $n\to\infty$ , where $x^{\dagger}$ denotes the unique ${\mathcal{R}}$ -minimizing solution of (1.1).

Proof.

Let $\hat{x}$ be any solution of (1.1) in $\mbox{dom}({\mathcal{R}})$ and define

\Delta_{n}:=D_{\mathcal{R}}^{\xi_{n}}(\hat{x},x_{n}).

By the similar argument in the proof of Lemma 3.1 we can obtain

\Delta_{n+1}-\Delta_{n}\leq-c_{0}t_{n}\|A_{I_{n}}x_{n}-y_{I_{n}}\|^{2}\leq-c_{0}\mu_{2}\|A_{I_{n}}x_{n}-y_{I_{n}}\|^{2}.

Consequently

		$\displaystyle\mathbb{E}[\Delta_{n+1}]-\mathbb{E}[\Delta_{n}]$
		$\displaystyle\leq-c_{0}\mu_{2}\mathbb{E}\left[\\|A_{I_{n}}x_{n}-y_{I_{n}}\\|^{2}\right]=-c_{0}\mu_{2}\mathbb{E}\left[\mathbb{E}\left[\\|A_{I_{n}}x_{n}-y_{I_{n}}\\|^{2}\|{\mathcal{F}_{n}}\right]\right]$
		$\displaystyle=-c_{0}\mu_{2}\mathbb{E}\left[\frac{1}{{p\choose b}}\sum_{I:\|I\|=b}\\|A_{I}x_{n}-y_{I}\\|^{2}\right].$		(3.3)

This shows that $\{\mathbb{E}[\Delta_{n}]\}$ is monotonically decreasing and therefore

\displaystyle\lim_{n\to\infty}\mathbb{E}[\Delta_{n}]\mbox{ exists}\quad\mbox{and}\quad\lim_{n\to\infty}\Phi_{n}=0

(3.4)

where

\Phi_{n}:=\mathbb{E}\left[\frac{1}{{p\choose b}}\sum_{I:|I|=b}\|A_{I}x_{n}-y_{I}\|^{2}\right].

If $\Phi_{n}=0$ for some $n$ , we must have

\sum_{I:|I|=b}\|A_{I}x_{n}-y_{I}\|^{2}=0

along any sample path $(I_{0},\cdots,I_{n-1})$ since there exist only finite many such sample paths each with a positive probability; consequently $\xi_{m}=\xi_{n}$ and $x_{m}=x_{n}$ and hence $\Phi_{m}=0$ for all $m\geq n$ . Based on these properties of $\Phi_{n}$ , it is possible to choose a strictly increasing sequence $\{n_{l}\}$ of integers by setting $n_{0}:=0$ and, for each $l\geq 1$ , by letting $n_{l}$ be the first integer satisfying

n_{l}\geq n_{l-1}+1\quad\mbox{ and }\quad\Phi_{n_{l}}\leq\Phi_{n_{l-1}}.

It is easy to see that for this sequence $\{n_{l}\}$ there holds

\displaystyle\Phi_{n_{l}}\leq\Phi_{n},\quad\forall 0\leq n\leq n_{l}.

(3.5)

For the above chosen sequence $\{n_{l}\}$ of integers, we are now going to show that

\displaystyle\sup_{l\geq k}\mathbb{E}\left[D_{\mathcal{R}}^{\xi_{n_{k}}}(x_{n_{l}},x_{n_{k}})\right]\to 0\quad\mbox{as }k\to\infty.

(3.6)

By the definition of the Bregman distance we have for any $l>k$ that

D_{\mathcal{R}}^{\xi_{n_{k}}}(x_{n_{l}},x_{n_{k}})=\Delta_{n_{k}}-\Delta_{n_{l}}+\langle\xi_{n_{l}}-\xi_{n_{k}},x_{n_{l}}-\hat{x}\rangle.

Taking the expectation gives

\displaystyle\mathbb{E}\left[D_{\mathcal{R}}^{\xi_{n_{k}}}(x_{n_{l}},x_{n_{k}})\right]=\mathbb{E}[\Delta_{n_{k}}]-\mathbb{E}[\Delta_{n_{l}}]+\mathbb{E}\left[\langle\xi_{n_{l}}-\xi_{n_{k}},x_{n_{l}}-\hat{x}\rangle\right].

(3.7)

By using the definition of $\{\xi_{n}\}$ we have

	$\displaystyle\langle\xi_{n_{l}}-\xi_{n_{k}},x_{n_{l}}-\hat{x}\rangle$	$\displaystyle=\sum_{n=n_{k}}^{n_{l}-1}\langle\xi_{n+1}-\xi_{n},x_{n_{l}}-\hat{x}\rangle$
		$\displaystyle=-\sum_{n=n_{k}}^{n_{l}-1}t_{n}\langle A_{I_{n}}^{*}(A_{I_{n}}x_{n}-y_{I_{n}}),x_{n_{l}}-\hat{x}\rangle$
		$\displaystyle=-\sum_{n=n_{k}}^{n_{l}-1}t_{n}\langle A_{I_{n}}x_{n}-y_{I_{n}},A_{I_{n}}x_{n_{l}}-y_{I_{n}}\rangle.$

Therefore, by taking the expectation and using the Cauchy-Schwarz inequality, we can obtain

	$\displaystyle\left\|\mathbb{E}\left[\langle\xi_{n_{l}}-\xi_{n_{k}},x_{n_{l}}-\hat{x}\rangle\right]\right\|\leq\mu_{1}\sum_{n=n_{k}}^{n_{l}-1}\left\|\mathbb{E}\left[\langle A_{I_{n}}x_{n}-y_{I_{n}},A_{I_{n}}x_{n_{l}}-y_{I_{n}}\rangle\right]\right\|$
	$\displaystyle\qquad\qquad\qquad\leq\mu_{1}\sum_{n=n_{k}}^{n_{l}-1}\left(\mathbb{E}\left[\\|A_{I_{n}}x_{n}-y_{I_{n}}\\|^{2}\right]\right)^{1/2}\left(\mathbb{E}\left[\\|A_{I_{n}}x_{n_{l}}-y_{I_{n}}\\|^{2}\right]\right)^{1/2}.$

Since $x_{n}$ is ${\mathcal{F}}_{n}$ -measurable, we have

\displaystyle\mathbb{E}\left[\|A_{I_{n}}x_{n}-y_{I_{n}}\|^{2}\right]=\mathbb{E}\left[\mathbb{E}\left[\|A_{I_{n}}x_{n}-y_{I_{n}}\|^{2}|{\mathcal{F}}_{n}\right]\right]=\Phi_{n}.

We can not treat the term $\mathbb{E}\left[\|A_{I_{n}}x_{n_{l}}-y_{I_{n}}\|^{2}\right]$ in the same way because $x_{n_{l}}$ is not necessarily ${\mathcal{F}}_{n}$ -measurable. However, by noting that

\|A_{I_{n}}x_{n_{l}}-y_{I_{n}}\|^{2}\leq\sum_{I:|I|=b}\|A_{I}x_{n_{l}}-y_{I}\|^{2},

we have

\mathbb{E}\left[\|A_{I_{n}}x_{n_{l}}-y_{I_{n}}\|^{2}\right]\leq\mathbb{E}\left[\sum_{I:|I|=b}\|A_{I}x_{n_{l}}-y_{I}\|^{2}\right]\leq{p\choose b}\Phi_{n_{l}}.

Therefore, with $c_{p,b}:={p\choose b}^{1/2}$ , we obtain

\displaystyle\left|\mathbb{E}\left[\langle\xi_{n_{l}}-\xi_{n_{k}},x_{n_{l}}-\hat{x}\rangle\right]\right|\leq\mu_{1}c_{p,b}\sum_{n=n_{k}}^{n_{l}-1}\Phi_{n}^{1/2}\Phi_{n_{l}}^{1/2}.

(3.8)

By virtue of (3.5) and (3) we then have

\displaystyle\left|\mathbb{E}\left[\langle\xi_{n_{l}}-\xi_{n_{k}},x_{n_{l}}-\hat{x}\rangle\right]\right|\leq\mu_{1}c_{p,b}\sum_{n=n_{k}}^{n_{l}-1}\Phi_{n}\leq\frac{\mu_{1}c_{p,b}}{c_{0}\mu_{2}}\left(\mathbb{E}[\Delta_{n_{k}}]-\mathbb{E}[\Delta_{n_{l}}]\right).

(3.9)

Combining this with (3.7) yields

\displaystyle\mathbb{E}\left[D_{\mathcal{R}}^{\xi_{n_{k}}}(x_{n_{l}},x_{n_{k}})\right]\leq\left(1+\frac{\mu_{1}c_{p,b}}{c_{0}\mu_{2}}\right)\left(\mathbb{E}[\Delta_{n_{k}}]-\mathbb{E}[\Delta_{n_{l}}]\right)

which together with the first equation in (3.4) shows (3.6) immediately.

By using the strong convexity of ${\mathcal{R}}$ , we can obtain from (3.6) that

\sup_{l\geq k}\mathbb{E}[\|x_{n_{l}}-x_{n_{k}}\|^{2}]\to 0\quad\mbox{as }k\to\infty

which means that $\{x_{n_{l}}\}$ is a Cauchy sequence in $L^{2}(\Omega,X)$ . Thus there exists a random vector $x^{*}\in L^{2}(\Omega,X)$ such that

\displaystyle\mathbb{E}[\|x_{n_{l}}-x^{*}\|^{2}]\to 0\quad\mbox{ as }l\to\infty.

(3.10)

By taking a subsequence of $\{n_{l}\}$ if necessary, we can obtain from (3.10) and the second equation in (3.4) that

\displaystyle\lim_{l\to\infty}\|x_{n_{l}}-x^{*}\|=0\quad\mbox{and}\quad\lim_{l\to\infty}\sum_{I:|I|=b}\|A_{I}x_{n_{l}}-y_{I}\|^{2}=0

(3.11)

almost surely. Consequently

\sum_{I:|I|=b}\|A_{I}x^{*}-y_{I}\|^{2}=0\ \ \mbox{almost surely},

i.e. $x^{*}$ is a solution of (1.1) almost surely.

We next show that $x^{*}\in\mbox{dom}({\mathcal{R}})$ almost surely. It suffices to show $\mathbb{E}[{\mathcal{R}}(x^{*})]<\infty$ . Recall that $\xi_{n_{l}}\in\partial{\mathcal{R}}(x_{n_{l}})$ , we have

\displaystyle{\mathcal{R}}(x_{n_{l}})\leq{\mathcal{R}}(x)+\langle\xi_{n_{l}},x_{n_{l}}-x\rangle,\quad\forall x\in X.

(3.12)

By using (3.12) with $x=x^{\dagger}$ , taking the expectation, and noting $\xi_{n_{0}}=\xi_{0}=0$ , we can obtain from (3.9) that

	$\displaystyle\mathbb{E}[{\mathcal{R}}(x_{n_{l}})]$	$\displaystyle\leq{\mathcal{R}}(x^{\dagger})+\mathbb{E}\left[\langle\xi_{n_{l}}-\xi_{n_{0}},x_{n_{l}}-x^{\dagger}\rangle\right]$
		$\displaystyle\leq{\mathcal{R}}(x^{\dagger})+\frac{\mu_{1}c_{p,b}}{c_{0}\mu_{2}}\mathbb{E}\left[D_{{\mathcal{R}}}^{\xi_{0}}(x^{\dagger},x_{0})\right]$
		$\displaystyle=:C<\infty.$

Therefore, it follows from the first equation in (3.11), the lower semi-continuity of ${\mathcal{R}}$ and Fatou’s Lemma that

\displaystyle\mathbb{E}[{\mathcal{R}}(x^{*})]\leq\mathbb{E}\left[\liminf_{l\to\infty}{\mathcal{R}}(x_{n_{l}})\right]\leq\liminf_{l\to\infty}\mathbb{E}\left[{\mathcal{R}}(x_{n_{l}})\right]\leq C<\infty.

(3.13)

We have thus shown that $x^{*}\in L^{2}(\Omega,X)$ is a solution of (1.1) in $\mbox{dom}({\mathcal{R}})$ almost surely.

In order to proceed further, we will show that, for any $\hat{x}\in L^{2}(\Omega,X)$ that is a solution of (1.1) in $\mbox{dom}({\mathcal{R}})$ almost surely, there holds

\displaystyle\lim_{l\to\infty}\mathbb{E}\left[\langle\xi_{n_{l}},x_{n_{l}}-\hat{x}\rangle\right]=0.

(3.14)

To see this, for any $k<l$ we write

\displaystyle\mathbb{E}\left[\langle\xi_{n_{l}},x_{n_{l}}-\hat{x}\rangle\right]=\mathbb{E}\left[\langle\xi_{n_{l}}-\xi_{n_{k}},x_{n_{l}}-\hat{x}\rangle\right]+\mathbb{E}\left[\langle\xi_{n_{k}},x_{n_{l}}-\hat{x}\rangle\right].

Since $\hat{x}\in L^{2}(\Omega,X)$ is a solution of (1.1) in $\mbox{dom}({\mathcal{R}})$ almost surely, by the definition of $\{\xi_{n}\}$ and $\xi_{0}=0$ we can use the similar argument for deriving (3.8) to obtain

\displaystyle\left|\mathbb{E}\left[\langle\xi_{n_{k}},x_{n_{l}}-\hat{x}\rangle\right]\right|\leq\mu_{1}c_{p,b}\left(\sum_{n=0}^{n_{k}-1}\Phi_{n}^{1/2}\right)\Phi_{n_{l}}^{1/2}.

This and the second equation in (3.4) imply for any fixed $k$ that $\mathbb{E}[\langle\xi_{n_{k}},x_{n_{l}}-\hat{x}\rangle]\to 0$ as $l\to\infty$ . Hence

\displaystyle\limsup_{l\to\infty}\left|\mathbb{E}\left[\langle\xi_{n_{l}},x_{n_{l}}-\hat{x}\rangle\right]\right|

\displaystyle\leq\sup_{l\geq k}\left|\mathbb{E}\left[\langle\xi_{n_{l}}-\xi_{n_{k}},x_{n_{l}}-\hat{x}\rangle\right]\right|

for all $k$ . By virtue of (3.9) and letting $k\to\infty$ we thus obtain (3.14).

Based on (3.12) with $x=x^{*}$ and (3.14), we can obtain

\displaystyle\limsup_{l\to\infty}\mathbb{E}[{\mathcal{R}}(x_{n_{l}})]\leq\mathbb{E}\left[{\mathcal{R}}(x^{*})\right]+\lim_{l\to\infty}\mathbb{E}\left[\langle\xi_{n_{l}},x_{n_{l}}-x^{*}\rangle\right]=\mathbb{E}[{\mathcal{R}}(x^{*})]

which together with (3.13) then implies

\displaystyle\lim_{l\to\infty}\mathbb{E}[{\mathcal{R}}(x_{n_{l}})]=\mathbb{E}[{\mathcal{R}}(x^{*})]

(3.15)

By using (3.15), (3.12) with $x=x^{\dagger}$ , and (3.14) with $\hat{x}=x^{\dagger}$ , we have

\displaystyle\mathbb{E}[{\mathcal{R}}(x^{*})]=\lim_{l\to\infty}\mathbb{E}[{\mathcal{R}}(x_{n_{l}})]\leq{\mathcal{R}}(x^{\dagger})+\lim_{l\to\infty}\mathbb{E}[\langle\xi_{n_{l}},x_{n_{l}}-x^{\dagger}\rangle]={\mathcal{R}}(x^{\dagger}).

Since $x^{*}$ is a solution of (1.1) in $\mbox{dom}({\mathcal{R}})$ almost surely, we have ${\mathcal{R}}(x^{*})\geq{\mathcal{R}}(x^{\dagger})$ almost surely which implies $\mathbb{E}[{\mathcal{R}}(x^{*})]\geq{\mathcal{R}}(x^{\dagger})$ . Consequently $\mathbb{E}[{\mathcal{R}}(x^{*})]={\mathcal{R}}(x^{\dagger})$ and thus it follows from (3.15) that

\displaystyle\lim_{l\to\infty}\mathbb{E}[{\mathcal{R}}(x_{n_{l}})]={\mathcal{R}}(x^{\dagger}).

(3.16)

Finally, from (3.16) and (3.14) with $\hat{x}=x^{\dagger}$ it follows that

\displaystyle\lim_{l\to\infty}\mathbb{E}\left[D_{\mathcal{R}}^{\xi_{n_{l}}}(x^{\dagger},x_{n_{l}})\right]=\lim_{l\to\infty}\left(\mathbb{E}[{\mathcal{R}}(x^{\dagger})-{\mathcal{R}}(x_{n_{l}})-\langle\xi_{n_{l}},x^{\dagger}-x_{n_{l}}\rangle]\right)=0.

By the monotonicity of $\{\mathbb{E}[D_{\mathcal{R}}^{\xi_{n}}(x^{\dagger},x_{n})]\}$ we can conclude

\lim_{n\to\infty}\mathbb{E}\left[D_{\mathcal{R}}^{\xi_{n}}(x^{\dagger},x_{n})\right]=0

and hence $\lim_{n\to\infty}\mathbb{E}[\|x_{n}-x^{\dagger}\|^{2}]=0$ by the strong convexity of ${\mathcal{R}}$ . The proof is therefore complete. ∎

In order to use Theorem 3.2 to establish the convergence of Algorithm 1 with noisy data, we need to investigate the stability of the algorithm, i.e. the behavior of $x_{n}^{\delta}$ and $\xi_{n}^{\delta}$ as $\delta\to 0$ for each fixed $n$ . In the following we will provide stability results of Algorithm 1 under the following three choices of step-sizes:

(s1)

$t_{n}^{\delta}$ depends only on $I_{n}$ , i.e. $t_{n}^{\delta}=t_{I_{n}}$ and $0<t_{I}<4\sigma/\|A_{I}\|^{2}$ for all $I\subset\{1,\cdots,p\}$ with $|I|=b$ .

(s2)

$t_{n}^{\delta}$ is chosen according to the formula

t_{n}^{\delta}=\min\left\{\frac{\mu_{0}\|A_{I_{n}}x_{n}^{\delta}-y_{I_{n}}^{\delta}\|^{2}}{\|A_{I_{n}}^{*}(A_{I_{n}}x_{n}^{\delta}-y_{I_{n}}^{\delta})\|^{2}},\tilde{\mu}_{1}\right\},

where $\mu_{0}$ and $\tilde{\mu}_{1}$ are two positive constants with $0<\mu_{0}<4\sigma$ .

(s3)

In case the information of $\delta_{1},\cdots,\delta_{p}$ is available, choose the step-size $t_{n}^{\delta}$ according to the rule

\displaystyle t_{n}^{\delta}=\left\{\begin{array}[]{lll}\min\left\{\frac{\mu_{0}\|A_{I_{n}}x_{n}^{\delta}-y_{I_{n}}^{\delta}\|^{2}}{\|A_{I_{n}}^{*}(A_{I_{n}}x_{n}^{\delta}-y_{I_{n}}^{\delta})\|^{2}},\tilde{\mu}_{1}\right\}&\mbox{ if }\|A_{I_{n}}x_{n}^{\delta}-y_{I_{n}}^{\delta}\|>\tau\delta_{I_{n}},\\[5.16663pt] 0&\mbox{ otherwise},\end{array}\right.

where $\tau\geq 1$ , $\mu_{0}>0$ and $\tilde{\mu}_{1}>0$ are preassigned numbers with $0<\mu_{0}<4\sigma$ .

The step-size chosen in (s1) is motivated by the Landweber iteration which uses constant step-size. Along the iterations in Algorithm 1, when the same subset $I$ is repeatedly chosen, (s1) suggests to use the same step-size in computation. The step-sizes chosen by (s1) could be small and thus slows down the computation. The rules given in (s2) and (s3) use the adaptive strategy which may produce large step-sizes. The step-size chosen by (s2) is motivated by the deterministic minimal error method, see [24]. The step-size given in (s3) is motivated by the work on deterministic Landweber-Kaczmarz method, see [15, 20], and it incorporates the spirit of the discrepancy principle into the selection.

Let $\{x_{n},\xi_{n}\}$ be defined by Algorithm 1 with exact data, where the step-size $t_{n}$ is chosen by (s1), (s2) or (s3) with the superscript $\delta$ dropped and with $\delta_{I_{n}}$ replaced by $0$ . It is easy to see that $t_{n}$ satisfies (3.2) in Theorem 3.2. Therefore, Theorem 3.2 is applicable to $\{x_{n},\xi_{n}\}$ .

The following result gives a stability property of Algorithm 1 with the step-sizes chosen by (s1), (s2) or (s3).

Lemma 3.3.

Let Assumption 1 hold. Consider Algorithm 1 with the step-size chosen by (s1), (s2) or (s3). Then for each fixed integer $n\geq 0$ there hold

\mathbb{E}[\|x_{n}^{\delta}-x_{n}\|^{2}]\to 0\quad\mbox{and}\quad\mathbb{E}[\|\xi_{n}^{\delta}-\xi_{n}\|^{2}]\to 0

as $\delta\to 0$ .

Proof.

We show the result by induction on $n$ . The result is trivial for $n=0$ because $\xi_{0}^{\delta}=\xi_{0}$ . Assuming the result holds for some $n\geq 0$ , we will show it also holds for $n+1$ . Since there are only finitely many sample paths of the form $(I_{0},\cdots,I_{n-1})$ each with positive probability, we can conclude from the induction hypothesis that

\displaystyle\|x_{n}^{\delta}-x_{n}\|\to 0\quad\mbox{ and }\quad\|\xi_{n}^{\delta}-\xi_{n}\|\to 0\quad\mbox{as }\delta\to 0

(3.17)

along every sample path $(I_{0},\cdots,I_{n-1})$ . We now show that along each sample path $(I_{0},\cdots,I_{n-1},I_{n})$ there hold

\displaystyle\|x_{n+1}^{\delta}-x_{n+1}\|\to 0\quad\mbox{ and }\quad\|\xi_{n+1}^{\delta}-\xi_{n+1}\|\to 0\quad\mbox{as }\delta\to 0.

(3.18)

To see this, we first show $\|\xi_{n+1}^{\delta}-\xi_{n+1}\|\to 0$ as $\delta\to 0$ by considering two cases.

We first consider the case $A_{I_{n}}x_{n}-y_{I_{n}}=0$ . For this case we have $\xi_{n+1}=\xi_{n}$ and hence

\displaystyle\xi_{n+1}^{\delta}-\xi_{n+1}=\xi_{n}^{\delta}-\xi_{n}-t_{n}^{\delta}A_{I_{n}}^{*}\left(A_{I_{n}}(x_{n}^{\delta}-x_{n})-(y_{I_{n}}^{\delta}-y_{I_{n}})\right).

Note that for all the choices of the step-sizes given in (s1), (s2) and (s3) we have $0\leq t_{n}^{\delta}\leq\hat{\mu}_{1}$ , where $\hat{\mu}_{1}$ is a constant independent of $n$ and $\delta$ . Note also that $\delta_{I_{n}}\leq\delta$ . Therefore

\displaystyle\|\xi_{n+1}^{\delta}-\xi_{n+1}\|\leq\|\xi_{n}^{\delta}-\xi_{n}\|+\hat{\mu}_{1}\|A_{I_{n}}\|\left(\|A_{I_{n}}\|\|x_{n}^{\delta}-x_{n}\|+\delta\right).

Consequently, by using (3.17) we can obtain $\|\xi_{n+1}^{\delta}-\xi_{n+1}\|\to 0$ as $\delta\to 0$ .

Next we consider the case $A_{I_{n}}x_{n}-y_{I_{n}}\neq 0$ . By using (3.17) we have $\|A_{I_{n}}x_{n}^{\delta}-y_{I_{n}}^{\delta}\|>\tau\delta_{I_{n}}$ when $\delta>0$ is sufficiently small. Note that

\displaystyle\langle A_{I_{n}}^{*}(A_{I_{n}}x_{n}-y_{I_{n}}),x_{n}-x^{\dagger}\rangle=\|A_{I_{n}}x_{n}-y_{I_{n}}\|^{2}>0

which implies that $A_{I_{n}}^{*}(A_{I_{n}}x_{n}-y_{I_{n}})\neq 0$ . Thus the step-size chosen by (s2) or (s3) is given by

t_{n}^{\delta}=\min\left\{\frac{\mu_{0}\|A_{I_{n}}x_{n}^{\delta}-y_{I_{n}}^{\delta}\|^{2}}{\|A_{I_{n}}^{*}(A_{I_{n}}x_{n}^{\delta}-y_{I_{n}}^{\delta})\|^{2}},\tilde{\mu}_{1}\right\}.

By using again (3.17) we have

\displaystyle t_{n}^{\delta}\to t_{n}\quad\mbox{ as }\delta\to 0.

(3.19)

For the step-size chosen by (s1) the equation (3.19) holds trivially as $t_{n}^{\delta}=t_{n}=t_{I_{n}}$ . Therefore, by noting that

	$\displaystyle\xi_{n+1}^{\delta}-\xi_{n+1}$	$\displaystyle=\xi_{n}^{\delta}-\xi_{n}-t_{n}^{\delta}A_{I_{n}}^{}(A_{I_{n}}x_{n}^{\delta}-y_{I_{n}}^{\delta})+t_{n}A_{I_{n}}^{}(A_{I_{n}}x_{n}-y_{I_{n}})$
		$\displaystyle=\xi_{n}^{\delta}-\xi_{n}+(t_{n}-t_{n}^{\delta})A_{I_{n}}^{*}(A_{I_{n}}x_{n}-y_{I_{n}})$
		$\displaystyle\quad\,+t_{n}^{\delta}A_{I_{n}}^{*}\left(A_{I_{n}}(x_{n}-x_{n}^{\delta})+(y_{I_{n}}^{\delta}-y_{I_{n}})\right),$

we can use (3.17) and (3.19) to conclude again $\|\xi_{n+1}^{\delta}-\xi_{n+1}\|\to 0$ as $\delta\to 0$ .

Now by virtue of the formula $x_{n+1}^{\delta}=\nabla{\mathcal{R}}^{*}(\xi_{n+1}^{\delta})$ , $x_{n+1}=\nabla{\mathcal{R}}^{*}(\xi_{n+1})$ and the continuity of $\nabla{\mathcal{R}}^{*}$ , we have $\|x_{n+1}^{\delta}-x_{n+1}\|\to 0$ as $\delta\to 0$ . We thus obtain (3.18).

Finally, since there are only finitely many sample paths of the form $(I_{0},\cdots,I_{n})$ , we can use (3.18) to conclude

\mathbb{E}\left[\|x_{n+1}^{\delta}-x_{n+1}\|^{2}\right]\to 0\quad\mbox{and}\quad\mathbb{E}\left[\|\xi_{n+1}^{\delta}-\xi_{n+1}\|^{2}\right]\to 0

as $\delta\to 0$ . The proof is thus complete by the induction principle. ∎

Theorem 3.4.

Let Assumption 1 hold. Consider Algorithm 1 with the step-size chosen by (s1), (s2) or (s3). If the integer $n_{\delta}$ is chosen such that $n_{\delta}\to\infty$ and $\delta^{2}n_{\delta}\to 0$ as $\delta\to 0$ , then

\mathbb{E}\left[\|x_{n_{\delta}}^{\delta}-x^{\dagger}\|^{2}\right]\to 0\quad\mbox{ and }\quad\mathbb{E}\left[D_{\mathcal{R}}^{\xi_{n_{\delta}}^{\delta}}(x^{\dagger},x_{n_{\delta}}^{\delta})\right]\to 0

as $\delta\to 0$ .

Proof.

Let $\Delta_{n}^{\delta}:=D_{\mathcal{R}}^{\xi_{n}^{\delta}}(x^{\dagger},x_{n}^{\delta})$ for all $n\geq 0$ . By the strong convexity of ${\mathcal{R}}$ it suffices to show $\mathbb{E}[\Delta_{n_{\delta}}^{\delta}]\to 0$ as $\delta\to 0$ . According to Lemma 3.1 we have

\mathbb{E}\left[\Delta_{n+1}^{\delta}\right]-\mathbb{E}\left[\Delta_{n}^{\delta}\right]\leq\frac{\mu_{1}b}{4c_{0}p}\delta^{2}

for all $n\geq 0$ . Let $n$ be any fixed integer. Since $n_{\delta}\to\infty$ as $\delta\to 0$ , we have $n_{\delta}>n$ for small $\delta>0$ . Thus, we may repeatedly use the above inequality to obtain

\displaystyle\mathbb{E}\left[\Delta_{n_{\delta}}^{\delta}\right]\leq\mathbb{E}\left[\Delta_{n}^{\delta}\right]+\frac{\mu_{1}b}{4c_{0}p}(n_{\delta}-n)\delta^{2}.

Since $\delta^{2}n_{\delta}\to 0$ as $\delta\to 0$ , we thus have

	$\displaystyle\limsup_{\delta\to 0}\mathbb{E}\left[\Delta_{n_{\delta}}^{\delta}\right]$
	$\displaystyle\leq\limsup_{\delta\to 0}\mathbb{E}\left[\Delta_{n}^{\delta}\right]=\limsup_{\delta\to 0}\left({\mathcal{R}}(x^{\dagger})-\mathbb{E}\left[{\mathcal{R}}(x_{n}^{\delta})\right]-\mathbb{E}\left[\langle\xi_{n}^{\delta},x^{\dagger}-x_{n}^{\delta}\rangle\right]\right)$

for any $n\geq 0$ . With the help of Cauchy-Schwarz inequality we have

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

Thus, we may use Lemma 3.3 to conclude

\displaystyle\lim_{\delta\to 0}\mathbb{E}\left[\langle\xi_{n}^{\delta},x^{\dagger}-x_{n}^{\delta}\rangle\right]=\mathbb{E}\left[\langle\xi_{n},x^{\dagger}-x_{n}\rangle\right].

According to the proof of Lemma 3.3 we have $\|x_{n}^{\delta}-x_{n}\|\to 0$ as $\delta\to 0$ along each sample path $(I_{0},\cdots,I_{n-1})$ . Thus, from the lower semi-continuity of ${\mathcal{R}}$ and Fatou’s lemma it follows

\mathbb{E}[{\mathcal{R}}(x_{n})]\leq\mathbb{E}\left[\liminf_{\delta\to 0}{\mathcal{R}}(x_{n}^{\delta})\right]\leq\liminf_{\delta\to 0}\mathbb{E}\left[{\mathcal{R}}(x_{n}^{\delta})\right].

Consequently

	$\displaystyle\limsup_{\delta\to 0}\mathbb{E}\left[\Delta_{n_{\delta}}^{\delta}\right]$	$\displaystyle\leq{\mathcal{R}}(x^{\dagger})-\liminf_{\delta\to 0}\mathbb{E}\left[{\mathcal{R}}(x_{n}^{\delta})\right]-\lim_{\delta\to 0}\mathbb{E}\left[\langle\xi_{n}^{\delta},x^{\dagger}-x_{n}^{\delta}\rangle\right]$
		$\displaystyle\leq{\mathcal{R}}(x^{\dagger})-\mathbb{E}[{\mathcal{R}}(x_{n})]-\mathbb{E}\left[\langle\xi_{n},x^{\dagger}-x_{n}\rangle\right]$
		$\displaystyle=\mathbb{E}[D_{\mathcal{R}}^{\xi_{n}}(x^{\dagger},x_{n})]$

for all $n\geq 0$ . Letting $n\to\infty$ and using Theorem 3.2 we therefore obtain $\mathbb{E}[\Delta_{n_{\delta}}^{\delta}]\to 0$ as $\delta\to 0$ . ∎

Theorem 3.4 provides convergence results on Algorithm 1 under various choices of the step-sizes. For the step-size chosen by (s3), if further conditions are imposed and $\tau$ and $\mu_{0}$ , we can show that the iteration error of Algorithm 1 always decreases.

Lemma 3.5.

Let Assumption 1 hold. Consider Algorithm 1 with the step-size chosen by (s3). If $\mu_{0}>0$ and $\tau>1$ are chosen such that

\displaystyle 1-\frac{1}{\tau}-\frac{\mu_{0}}{4\sigma}\geq 0,

(3.20)

then $\Delta_{n+1}^{\delta}\leq\Delta_{n}^{\delta}$ and consequently $\mathbb{E}\left[\Delta_{n+1}^{\delta}\right]\leq\mathbb{E}\left[\Delta_{n}^{\delta}\right]$ for all $n\geq 0$ , where $\Delta_{n}^{\delta}:=D_{\mathcal{R}}^{\xi_{n}^{\delta}}(x^{\dagger},x_{n}^{\delta})$ .

Proof.

According to the proof of Lemma 3.1, we have

	$\displaystyle\Delta_{n+1}^{\delta}-\Delta_{n}^{\delta}$	$\displaystyle\leq\frac{1}{4\sigma}\left(t_{n}^{\delta}\right)^{2}\\|A_{I_{n}}^{*}(A_{I_{n}}x_{n}^{\delta}-y_{I_{n}}^{\delta})\\|^{2}-t_{n}^{\delta}\langle A_{I_{n}}x_{n}^{\delta}-y_{I_{n}}^{\delta},A_{I_{n}}x_{n}^{\delta}-y_{I_{n}}\rangle$
		$\displaystyle\leq\frac{1}{4\sigma}\left(t_{n}^{\delta}\right)^{2}\\|A_{I_{n}}^{*}(A_{I_{n}}x_{n}^{\delta}-y_{I_{n}}^{\delta})\\|^{2}-t_{n}^{\delta}\\|A_{I_{n}}x_{n}^{\delta}-y_{I_{n}}^{\delta}\\|^{2}$
		$\displaystyle\quad\,+\delta_{I_{n}}t_{n}^{\delta}\\|A_{I_{n}}x_{n}^{\delta}-y_{I_{n}}^{\delta}\\|.$

By the definition of $t_{n}^{\delta}$ we have

	$\displaystyle t_{n}^{\delta}\\|A_{I_{n}}^{*}(A_{I_{n}}x_{n}^{\delta}-y_{I_{n}}^{\delta})\\|^{2}\leq\mu_{0}\\|A_{I_{n}}x_{n}^{\delta}-y_{I_{n}}^{\delta}\\|^{2},$
	$\displaystyle\delta_{I_{n}}t_{n}^{\delta}\\|A_{I_{n}}x_{n}^{\delta}-y_{I_{n}}^{\delta}\\|\leq\frac{1}{\tau}t_{n}^{\delta}\\|A_{I_{n}}x_{n}^{\delta}-y_{I_{n}}^{\delta}\\|^{2}.$

Therefore

\displaystyle\Delta_{n+1}^{\delta}-\Delta_{n}^{\delta}\leq-\left(1-\frac{1}{\tau}-\frac{\mu_{0}}{4\sigma}\right)t_{n}^{\delta}\|A_{I_{n}}x_{n}^{\delta}-y_{I_{n}}^{\delta}\|^{2}.

Since $\mu_{0}$ and $\tau$ satisfy (3.20), this implies $\Delta_{n+1}^{\delta}\leq\Delta_{n}^{\delta}$ . By taking the expectation we then obtain $\mathbb{E}[\Delta_{n+1}^{\delta}]\leq\mathbb{E}[\Delta_{n}^{\delta}]$ . ∎

Theorem 3.6.

Let Assumption 1 hold. Consider Algorithm 1 with step-size chosen by (s3), where $\tau>1$ and $\mu_{0}>0$ are chosen to satisfy (3.20). If the integer $n_{\delta}$ is chosen such that $n_{\delta}\to\infty$ as $\delta\to 0$ , then

\mathbb{E}\left[\|x_{n_{\delta}}^{\delta}-x^{\dagger}\|^{2}\right]\to 0\quad\mbox{ and }\quad\mathbb{E}\left[D_{\mathcal{R}}^{\xi_{n_{\delta}}^{\delta}}(x^{\dagger},x_{n_{\delta}}^{\delta})\right]\to 0

as $\delta\to 0$ .

Proof.

Let $n\geq 0$ be any fixed integer. Since $n_{\delta}\to\infty$ as $\delta\to 0$ , we have $n_{\delta}>n$ for sufficiently small $\delta>0$ . It then follows from Lemma 3.5 that

\mathbb{E}\left[\Delta_{n_{\delta}}^{\delta}\right]\leq\mathbb{E}\left[\Delta_{n}^{\delta}\right]

By using Lemma 3.3 and following the proof of Theorem 3.4 we can further conclude

\limsup_{\delta\to 0}\mathbb{E}\left[\Delta_{n_{\delta}}^{\delta}\right]\leq\mathbb{E}[\Delta_{n}]

for all integers $n\geq 0$ . Finally, by letting $n\to\infty$ , we can use Theorem 3.2 to obtain $\mathbb{E}[\Delta_{n_{\delta}}^{\delta}]\to 0$ as $\delta\to 0$ . ∎

We remark that, unlike Theorem 3.4, the convergence result given in Theorem 3.6 requires only $n_{\delta}$ to satisfy $n_{\delta}\to\infty$ as $\delta\to 0$ . This is not surprising because the discrepancy principle has been incorporated into the choice of the step-size by (s3). Theorem 3.6 can be viewed as a stochastic extension of the corresponding result for the deterministic Landweber-Kaczmarz method ([20]) for solving ill-posed problems. It should be pointed out that, due to (3.20), the application of Theorem 3.6 requires either $\tau$ to be sufficiently large or $\mu_{0}$ to be sufficiently small which may lead the corresponding algorithm to lose accuracy or converge slowly. However, our numerical simulations in Section 5 demonstrate that using step-sizes chosen by (s3) without satisfying (3.20) still has the effect of decreasing errors as the iteration proceeds.

4. Rate of convergence

In Section 3 we have proved the convergence of Algorithm 1 under various choices of the step-sizes. In this section we will focus on deriving convergence rates. For ill-posed problems, this requires the sought solution to satisfy suitable source conditions. We will consider the benchmark source condition of the form ([7])

\displaystyle A^{*}\lambda^{\dagger}\in\partial{\mathcal{R}}(x^{\dagger})

(4.1)

for some $\lambda^{\dagger}:=(\lambda_{,1}^{\dagger},\cdots,\lambda_{,p}^{\dagger})\in Y:=Y_{1}\times\cdots\times Y_{p}$ , where, here and below, for any $\lambda\in Y$ we use $\lambda_{,i}$ to denote its $i$ -th component in $Y_{i}$ . Deriving convergence rates under general choices of step-sizes can be very challenging. Therefore, in this section we will restrict ourselves to the step-sizes chosen by (s1) that includes the constant step-sizes. The corresponding algorithm can be reformulated as follows.

Algorithm 2.

Fix a batch size $b$ , pick the initial guess $\xi_{0}=0$ in $X^{*}$ and step-sizes $t_{I}$ for all $I\subset\{1,\cdots,p\}$ with $|I|=b$ . Set $\xi_{0}^{\delta}:=\xi_{0}$ . For $n\geq 0$ do the following:

(i)

Calculate $x_{n}^{\delta}\in X$ by solving (2.10);
(ii)

Randomly select a subset $I_{n}\subset\{1,\cdots,p\}$ with $|I_{n}|=b$ via the uniform distribution;
(iii)

Define $\xi_{n+1}^{\delta}\in X^{*}$ by (2.11).

For Algorithm 2 we have the following convergence rate result.

Theorem 4.1.

Let Assumption 1 hold and consider Algorithm 2. Assume that $0<t_{I}<4\sigma/\|A_{I}\|^{2}$ for all subsets $I\subset\{1,\cdots,p\}$ with $|I|=b$ ; furthermore $t_{I}=t$ with a constant $t$ for all such $I$ in case $b>1$ . If the sought solution $x^{\dagger}$ satisfies the source condition (4.1) and the integer $n_{\delta}$ is chosen such that $1+\frac{b}{p}n_{\delta}$ is of the magnitude of $\delta^{-1}$ , then

{\mathbb{E}}\left[\|x_{n_{\delta}}^{\delta}-x^{\dagger}\|^{2}\right]\leq C\delta,

where $C$ is a positive constant independent of $\delta$ , $p$ and $b$ .

When $X$ is a Hilbert space and ${\mathcal{R}}(x)=\|x\|^{2}/2$ , the stochastic mirror descent method becomes the stochastic gradient method. In the existing literature on stochastic gradient method for solving ill-posed problems, sub-optimal convergence rates have been derived for diminishing step-sizes under general source conditions on the sought solution, see [18, 25]. Our result given in Theorem 4.1 supplements these results by demonstrating that the stochastic gradient method can achieve the order optimal convergence rate $O(\delta)$ under constant step-sizes if the sought solution satisfies the source condition (4.1).

The result in Theorem 4.1 also demonstrates the role played by the batch size $b$ : to achieve the same convergence rate, less number of iterations is required if a larger batch size $b$ is used. However, this does not mean the corresponding algorithm with larger batch size $b$ is faster because the computational time at each iteration can increase as $b$ increases. How to determine a batch size with the best performance is an outstanding challenging question.

The proof of Theorem 4.1 is based on considering an equivalent formulation of Algorithm 2 which we will derive in the following. The derivation is based on applying the randomized block gradient method ([26, 32, 36]) to the dual problem of (1.4) with $y$ replaced by $y^{\delta}$ . The associated Lagrangian function is

L(x;\lambda):={\mathcal{R}}(x)-\langle\lambda,Ax-y^{\delta}\rangle,

where $x\in X$ and $\lambda\in Y$ . Thus the dual function is

\displaystyle\inf_{x\in X}L(x;\lambda)

\displaystyle=-\sup_{x\in X}\left\{\left\langle A^{*}\lambda,x\right\rangle-{\mathcal{R}}(x)\right\}+\langle\lambda,y^{\delta}\rangle=-{\mathcal{R}}^{*}\left(A^{*}\lambda\right)+\langle\lambda,y^{\delta}\rangle

which means the dual problem is

\displaystyle\min_{\lambda\in Y}\left\{d_{y^{\delta}}(\lambda):={\mathcal{R}}^{*}\left(A^{*}\lambda\right)-\langle\lambda,y^{\delta}\rangle\right\}.

(4.2)

Note that ${\mathcal{R}}^{*}$ is continuous differentiable and

d_{y^{\delta}}(\lambda)={\mathcal{R}}^{*}\left(\sum_{i=1}^{p}A_{i}^{*}\lambda_{,i}\right)-\sum_{i=1}^{p}\langle\lambda_{,i},y_{i}^{\delta}\rangle.

Therefore we may solve (4.2) by the randomized block gradient method which iteratively selects partial components of $\lambda$ at random to be updated by a partial gradient of $d_{y^{\delta}}$ and leave other components unchanged. To be more precisely, for any $\lambda\in Y$ and any subset $I\subset\{1,\cdots,p\}$ , let $\lambda_{I}$ denote the group of components $\lambda_{,i}$ of $\lambda$ with $i\in I$ . Assume $\lambda_{n}^{\delta}:=(\lambda_{n,1}^{\delta},\cdots,\lambda_{n,p}^{\delta})\in Y$ is a current iterate. We then choose a subset of indices $I_{n}$ from $\{1,\cdots,p\}$ with $|I_{n}|=b$ randomly with uniform distribution and define $\lambda_{n+1}^{\delta}:=(\lambda_{n+1,1}^{\delta},\cdots,\lambda_{n+1,p}^{\delta})$ by setting $\lambda_{n+1,i}^{\delta}=\lambda_{n,i}^{\delta}$ if $i\not\in I_{n}$ and

\displaystyle\lambda_{n+1,I_{n}}^{\delta}=\lambda_{n,I_{n}}^{\delta}-t_{I_{n}}\nabla_{I_{n}}d_{y^{\delta}}(\lambda_{n}^{\delta}),

(4.3)

with step-sizes $t_{I_{n}}>0$ depending only on $I_{n}$ , where, for each $I\subset\{1,\cdots,p\}$ , $\nabla_{I}d_{y^{\delta}}$ denotes the partial gradient of $d_{y^{\delta}}$ with respect to $\lambda_{I}$ . Note that

\nabla_{I}d_{y^{\delta}}(\lambda)=A_{I}\nabla{\mathcal{R}}^{*}\left(A^{*}\lambda\right)-y_{I}^{\delta}.

Therefore, by setting $x_{n}^{\delta}:=\nabla{\mathcal{R}}^{*}\left(A^{*}\lambda_{n}^{\delta}\right)$ , we can rewrite (4.3) as

\lambda_{n+1,I_{n}}^{\delta}=\lambda_{n,I_{n}}^{\delta}-t_{I_{n}}\left(A_{I_{n}}x_{n}^{\delta}-y_{I_{n}}^{\delta}\right).

By using the definition of $x_{n}^{\delta}$ and (2.4) we have $A^{*}\lambda_{n}^{\delta}\in\partial{\mathcal{R}}(x_{n}^{\delta})$ which implies

0\in\partial\left({\mathcal{R}}-\langle\lambda_{n}^{\delta},A\cdot\rangle\right)(x_{n}^{\delta})

and thus

\displaystyle x_{n}^{\delta}=\arg\min_{x\in X}\left\{{\mathcal{R}}(x)-\langle\lambda_{n}^{\delta},Ax\rangle\right\}.

(4.4)

Combining the above analysis we thus obtain the following mini-batch randomized dual block gradient method for solving (1.2) with noisy data.

Algorithm 3.

Fix a batch size $b$ , pick the initial guess $\lambda_{0}:=(0,\cdots,0)\in Y:=Y_{1}\times\cdots\times Y_{p}$ and step-sizes $t_{I}$ for all $I\subset\{1,\cdots,p\}$ with $|I|=b$ . Set $\lambda_{0}^{\delta}:=\lambda_{0}$ . For $n\geq 0$ do the following:

(i)

Calculate $x_{n}^{\delta}\in X$ by using (4.4);
(ii)

Randomly select a subset $I_{n}\subset\{1,\cdots,p\}$ with $|I_{n}|=b$ via the uniform distribution;
(iii)

Define $\lambda_{n+1}^{\delta}:=(\lambda_{n+1,1}^{\delta},\cdots,\lambda_{n+1,p}^{\delta})\in Y$ by setting $\lambda_{n+1,I_{n}^{c}}^{\delta}=\lambda_{n,I_{n}^{c}}^{\delta}$ and

$\displaystyle\lambda_{n+1,I_{n}}^{\delta}=\lambda_{n,I_{n}}^{\delta}-t_{I_{n}}\left(A_{I_{n}}x_{n}^{\delta}-y_{I_{n}}^{\delta}\right),$ (4.5)

where $I_{n}^{c}$ denotes the complement of $I_{n}$ in $\{1,\cdots,p\}$ .

Let us demonstrate the equivalence between Algorithm 2 and Algorithm 3. If $\{x_{n}^{\delta},\lambda_{n}^{\delta}\}$ is defined by Algorithm 3, by defining $\xi_{n}^{\delta}:=A^{*}\lambda_{n}^{\delta}$ we have

	$\displaystyle\xi_{n+1}^{\delta}$	$\displaystyle=A^{}\lambda_{n+1}^{\delta}=A_{I_{n}^{c}}^{}\lambda_{n+1,I_{n}^{c}}^{\delta}+A_{I_{n}}^{*}\lambda_{n+1,I_{n}}^{\delta}$
		$\displaystyle=A_{I_{n}^{c}}^{}\lambda_{n,I_{n}^{c}}^{\delta}+A_{I_{n}}^{}\left(\lambda_{n,I_{n}}^{\delta}-t_{I_{n}}(A_{I_{n}}x_{n}^{\delta}-y_{I_{n}}^{\delta})\right)$
		$\displaystyle=\left(A_{I_{n}^{c}}^{}\lambda_{n,I_{n}^{c}}^{\delta}+A_{I_{n}}^{}\lambda_{n,I_{n}}^{\delta}\right)-t_{I_{n}}A_{I_{n}}^{*}(A_{I_{n}}x_{n}^{\delta}-y_{I_{n}}^{\delta})$
		$\displaystyle=\xi_{n}^{\delta}-t_{I_{n}}A_{I_{n}}^{*}(A_{I_{n}}x_{n}^{\delta}-y_{I_{n}}^{\delta}).$

Therefore $\{x_{n}^{\delta},\xi_{n}^{\delta}\}$ can be produced by Algorithm 1. Conversely, if $\{x_{n}^{\delta},\xi_{n}^{\delta}\}$ is defined by Algorithm 1, by using $\xi_{0}^{\delta}=0$ one can easily see that $\xi_{n}^{\delta}\in\mbox{Ran}(A)$ , the range of $A$ . By writing $\xi_{n}^{\delta}=A^{*}\lambda_{n}^{\delta}$ for some $\lambda_{n}^{\delta}=(\lambda_{n,1}^{\delta},\cdots,\lambda_{n,p}^{\delta})\in Y$ , we have

	$\displaystyle\xi_{n+1}^{\delta}$	$\displaystyle=\xi_{n}^{\delta}-t_{I_{n}}A_{I_{n}}^{}(A_{I_{n}}x_{n}^{\delta}-y_{I_{n}}^{\delta})=A^{}\lambda_{n}^{\delta}-t_{I_{n}}A_{I_{n}}^{*}(A_{I_{n}}x_{n}^{\delta}-y_{I_{n}}^{\delta})$
		$\displaystyle=A_{I_{n}^{c}}^{}\lambda_{n,I_{n}^{c}}^{\delta}+A_{I_{n}}^{}\left(\lambda_{n,I_{n}}^{\delta}-t_{I_{n}}(A_{I_{n}}x_{n}^{\delta}-y_{I_{n}}^{\delta})\right)$

Define $\lambda_{n+1}^{\delta}\in Y$ by setting $\lambda_{n+1,I_{n}^{c}}^{\delta}=\lambda_{n,I_{n}^{c}}^{\delta}$ and $\lambda_{n+1,I_{n}}^{\delta}=\lambda_{n,I_{n}}^{\delta}-t_{I_{n}}(A_{I_{n}}x_{n}^{\delta}-y_{I_{n}}^{\delta})$ . We can see that $\xi_{n+1}^{\delta}=A^{*}\lambda_{n+1}^{\delta}$ and $\{x_{n}^{\delta},\lambda_{n}^{\delta}\}$ can be generated by Algorithm 3. Therefore we obtain

Lemma 4.2.

Algorithm 2 and Algorithm 3 are equivalent.

Based on Algorithm 3, we will derive the estimate on ${\mathbb{E}}[\|x_{n}^{\delta}-x^{\dagger}\|^{2}]$ under the benchmark source condition (4.1). Let

d_{y}(\lambda):={\mathcal{R}}^{*}\left(A^{*}\lambda\right)-\langle\lambda,y\rangle

which is obtained from $d_{y^{\delta}}(\lambda)$ with $y^{\delta}$ replaced by $y$ . From (4.1) and (2.4) it follows that $x^{\dagger}=\nabla{\mathcal{R}}^{*}(A^{*}\lambda^{\dagger})$ . Thus, by using $Ax^{\dagger}=y$ we have

\displaystyle\nabla d_{y}(\lambda^{\dagger})

\displaystyle=A\nabla{\mathcal{R}}^{*}\left(A^{*}\lambda^{\dagger}\right)-y=Ax^{\dagger}-y=0.

Since $d_{y}(\lambda)$ is convex, this shows that $\lambda^{\dagger}$ is a global minimizer of $d_{y}$ on $Y$ .

According to (4.4) we have $A^{*}\lambda_{n}^{\delta}\in\partial{\mathcal{R}}(x_{n}^{\delta})$ . Thus, we may consider the Bregman distance

\Delta_{n}^{\delta}:=D_{\mathcal{R}}^{A^{*}\lambda_{n}^{\delta}}(x^{\dagger},x_{n}^{\delta}).

By using (4.1), $A^{*}\lambda_{n}^{\delta}\in\partial{\mathcal{R}}(x_{n}^{\delta})$ , and (2.4) we have

\displaystyle{\mathcal{R}}(x^{\dagger})+{\mathcal{R}}^{*}\left(A^{*}\lambda^{\dagger}\right)=\left\langle A^{*}\lambda^{\dagger},x^{\dagger}\right\rangle,\quad{\mathcal{R}}(x_{n}^{\delta})+{\mathcal{R}}^{*}\left(A^{*}\lambda_{n}^{\delta}\right)=\left\langle A^{*}\lambda_{n}^{\delta},x_{n}^{\delta}\right\rangle.

Therefore

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

Consequently, by virtue of the strong convexity of ${\mathcal{R}}$ we have

\sigma\|x_{n}^{\delta}-x^{\dagger}\|^{2}\leq\Delta_{n}^{\delta}\leq d_{y}(\lambda_{n}^{\delta})-d_{y}(\lambda^{\dagger}).

Taking the expectation gives

\displaystyle\sigma\mathbb{E}\left[\|x_{n}^{\delta}-x^{\dagger}\|^{2}\right]\leq\mathbb{E}\left[d_{y}(\lambda_{n}^{\delta})-d_{y}(\lambda^{\dagger})\right].

(4.6)

This demonstrates that we can achieve our goal by bounding $\mathbb{E}[d_{y}(\lambda_{n}^{\delta})-d_{y}(\lambda^{\dagger})]$ .

Note that the sequence $\{\lambda_{n}^{\delta}\}$ is defined via the function $d_{y^{\delta}}(\lambda)$ . Although $\lambda^{\dagger}$ is a minimizer of $d_{y}$ , it may not be a minimizer of $d_{y^{\delta}}$ . In order to overcome this gap, we will make use of the relation

\displaystyle d_{y^{\delta}}(\lambda)-d_{y}(\lambda)=\langle\lambda,y-y^{\delta}\rangle

(4.7)

between $d_{y^{\delta}}(\lambda)$ and $d_{y}(\lambda)$ for all $\lambda\in Y$ .

Lemma 4.3.

Let Assumption 1 hold and consider Algorithm 3. Assume that $0<t_{I}<4\sigma/\|A_{I}\|^{2}$ for all $I\subset\{1,\dots,p\}$ with $|I|=b$ . Then

\mathbb{E}\left[d_{y}(\lambda_{n+1}^{\delta})\right]\leq\mathbb{E}\left[d_{y}(\lambda_{n}^{\delta})\right]+\frac{t_{\max}b}{2c_{1}p}\delta^{2}-\frac{c_{1}}{2{p\choose b}}\sum_{I:|I|=b}t_{I}\mathbb{E}\left[\|\nabla_{I}d_{y^{\delta}}(\lambda_{n}^{\delta})\|^{2}\right]

for all $n\geq 0$ , where

t_{\max}:=\max_{I:|I|=b}t_{I}\quad\mbox{ and }\quad c_{1}:=\min_{I:|I|=b}\left\{1-\frac{t_{I}\|A_{I}\|^{2}}{4\sigma}\right\}>0.

Proof.

For any $\lambda\in Y$ and any $I\subset\{1,\cdots,p\}$ with $|I|=b$ we may write $\lambda=(\lambda_{I},\lambda_{I^{c}})$ up to a permutation. Since the partial gradient of $d_{y^{\delta}}$ with respect to $\lambda_{I}$ is given by $\nabla_{I}d_{y^{\delta}}(\lambda)=A_{I}\nabla{\mathcal{R}}^{*}(A^{*}\lambda)-y_{I}^{\delta}$ , it follows from (2.12) that

\|\nabla_{I}d_{y^{\delta}}(\lambda_{I},\lambda_{I^{c}})-\nabla_{I}d_{y^{\delta}}(\lambda_{I}+h,\lambda_{I^{c}})\|\leq L_{I}\|h\|,\quad\forall h\in Y_{I},

where $L_{I}:=\|A_{I}\|^{2}/(2\sigma)$ . Consequently, by recalling the relation between $\lambda_{n+1}^{\delta}$ and $\lambda_{n}^{\delta}$ , we can obtain

	$\displaystyle d_{y^{\delta}}(\lambda_{n+1}^{\delta})$	$\displaystyle\leq d_{y^{\delta}}(\lambda_{n}^{\delta})+\langle\nabla_{I_{n}}d_{y^{\delta}}(\lambda_{n}^{\delta}),\lambda_{n+1,I_{n}}^{\delta}-\lambda_{n,I_{n}}^{\delta}\rangle+\frac{L_{I_{n}}}{2}\\|\lambda_{n+1,I_{n}}^{\delta}-\lambda_{n,I_{n}}^{\delta}\\|^{2}$
		$\displaystyle=d_{y^{\delta}}(\lambda_{n}^{\delta})-\left(1-\frac{L_{I_{n}}t_{I_{n}}}{2}\right)t_{I_{n}}\\|\nabla_{I_{n}}d_{y^{\delta}}(\lambda_{n}^{\delta})\\|^{2}.$

By using (4.7) and the definition of $\lambda_{n+1}^{\delta}$ , we can see that

	$\displaystyle d_{y^{\delta}}(\lambda_{n+1}^{\delta})-d_{y^{\delta}}(\lambda_{n}^{\delta})$	$\displaystyle=d_{y}(\lambda_{n+1}^{\delta})-d_{y}(\lambda_{n}^{\delta})+\langle\lambda_{n+1,I_{n}}^{\delta}-\lambda_{n,I_{n}}^{\delta},y_{I_{n}}-y_{I_{n}}^{\delta}\rangle$
		$\displaystyle=d_{y}(\lambda_{n+1}^{\delta})-d_{y}(\lambda_{n}^{\delta})-t_{I_{n}}\langle\nabla_{I_{n}}d_{y^{\delta}}(\lambda_{n}^{\delta}),y_{I_{n}}-y_{I_{n}}^{\delta}\rangle.$

Combining the above two equations and using the definition of $c_{1}$ , we thus have

	$\displaystyle d_{y}(\lambda_{n+1}^{\delta})-d_{y}(\lambda_{n}^{\delta})$
	$\displaystyle\leq t_{I_{n}}\langle\nabla_{I_{n}}d_{y^{\delta}}(\lambda_{n}^{\delta}),y_{I_{n}}-y_{I_{n}}^{\delta}\rangle-\left(1-\frac{L_{I_{n}}t_{I_{n}}}{2}\right)t_{I_{n}}\\|\nabla_{I_{n}}d_{y^{\delta}}(\lambda_{n}^{\delta})\\|^{2}$
	$\displaystyle\leq t_{I_{n}}\delta_{I_{n}}\\|\nabla_{I_{n}}d_{y^{\delta}}(\lambda_{n}^{\delta})\\|-c_{1}t_{I_{n}}\\|\nabla_{I_{n}}d_{y^{\delta}}(\lambda_{n}^{\delta})\\|^{2}$
	$\displaystyle\leq\frac{1}{2c_{1}}t_{I_{n}}\delta_{I_{n}}^{2}-\frac{1}{2}c_{1}t_{I_{n}}\\|\nabla_{I_{n}}d_{y^{\delta}}(\lambda_{n}^{\delta})\\|^{2}.$

By taking the expectation, using (3.1), and noting that $\mathbb{E}[\delta_{I_{n}}^{2}]=\frac{b}{p}\delta^{2}$ , we finally obtain

	$\displaystyle\mathbb{E}\left[d_{y}(\lambda_{n+1}^{\delta})-d_{y}(\lambda_{n}^{\delta})\right]$	$\displaystyle\leq\frac{t_{\max}}{2c_{1}}\mathbb{E}\left[\delta_{I_{n}}^{2}\right]-\frac{1}{2}c_{1}\mathbb{E}\left[\mathbb{E}\left[t_{I_{n}}\\|\nabla_{I_{n}}d_{y^{\delta}}(\lambda_{n}^{\delta})\\|^{2}\|{\mathcal{F}}_{n}\right]\right]$
		$\displaystyle=\frac{t_{\max}b}{2c_{1}p}\delta^{2}-\frac{1}{2}c_{1}\mathbb{E}\left[\frac{1}{{p\choose b}}\sum_{I:\|I\|=b}t_{I}\\|\nabla_{I}d_{y^{\delta}}(\lambda_{n}^{\delta})\\|^{2}\right]$

which completes the proof. ∎

Lemma 4.4.

Consider Algorithm 3 in which $0<t_{I}<4\sigma/\|A_{I}\|^{2}$ for all $I\subset\{1,\cdots,p\}$ with $|I|=b$ , furthermore $t_{I}=t$ with a constant $t$ for all such $I$ in case $b>1$ . Assume (4.1) holds and let

r_{n}:=\left\{\begin{array}[]{lll}\displaystyle{\sum_{i=1}^{p}\frac{1}{t_{i}}\|\lambda_{n,i}^{\delta}-\lambda_{,i}^{\dagger}\|^{2}}&\mbox{ if }b=1,\\[12.05553pt] \displaystyle{\frac{1}{t}\|\lambda_{n}^{\delta}-\lambda^{\dagger}\|^{2}}&\mbox{ if }b>1\end{array}\right.

for all $n\geq 0$ . Then

	$\displaystyle\mathbb{E}[r_{n+1}]$	$\displaystyle\leq\mathbb{E}[r_{n}]+\frac{2b}{p}\mathbb{E}\left[d_{y}(\lambda^{\dagger})-d_{y}(\lambda_{n}^{\delta})\right]+\frac{2t_{\max}^{1/2}b}{p}\delta\left(\mathbb{E}[r_{n}]\right)^{1/2}$
		$\displaystyle\quad\,+\frac{1}{{p\choose b}}\sum_{I:\|I\|=b}t_{I}\mathbb{E}\left[\\|\nabla_{I}d_{y^{\delta}}(\lambda_{n}^{\delta})\\|^{2}\right]$

for all integers $n\geq 0$ .

Proof.

We will only prove the result for $b>1$ ; the result for $b=1$ can be proved similarly. By using the definition of $\lambda_{n+1}^{\delta}$ and $t_{I_{n}}=t$ we have

	$\displaystyle r_{n+1}$	$\displaystyle=\frac{1}{t}\\|\lambda_{n+1,I_{n}^{c}}^{\delta}-\lambda_{I_{n}^{c}}^{\dagger}\\|^{2}+\frac{1}{t}\\|\lambda_{n+1,I_{n}}^{\delta}-\lambda_{I_{n}}^{\dagger}\\|^{2}$
		$\displaystyle=\frac{1}{t}\\|\lambda_{n,I_{n}^{c}}^{\delta}-\lambda_{I_{n}^{c}}^{\dagger}\\|^{2}+\frac{1}{t}\\|\lambda_{n,I_{n}}^{\delta}-\lambda_{I_{n}}^{\dagger}-t\nabla_{I_{n}}d_{y^{\delta}}(\lambda_{n}^{\delta})\\|^{2}$
		$\displaystyle=\frac{1}{t}\\|\lambda_{n,I_{n}^{c}}^{\delta}-\lambda_{I_{n}^{c}}^{\dagger}\\|^{2}+\frac{1}{t}\\|\lambda_{n,I_{n}}^{\delta}-\lambda_{I_{n}}^{\dagger}\\|^{2}-2\langle\nabla_{I_{n}}d_{y^{\delta}}(\lambda_{n}^{\delta}),\lambda_{n,I_{n}}^{\delta}-\lambda_{I_{n}}^{\dagger}\rangle$
		$\displaystyle\quad\,+t\\|\nabla_{I_{n}}d_{y^{\delta}}(\lambda_{n}^{\delta})\\|^{2}$
		$\displaystyle=r_{n}-2\langle\nabla_{I_{n}}d_{y^{\delta}}(\lambda_{n}^{\delta}),\lambda_{n,I_{n}}^{\delta}-\lambda_{I_{n}}^{\dagger}\rangle+t\\|\nabla_{I_{n}}d_{y^{\delta}}(\lambda_{n}^{\delta})\\|^{2}.$

Taking the expectation and using (3.1), it gives

	$\displaystyle\mathbb{E}[r_{n+1}]$	$\displaystyle=\mathbb{E}[r_{n}]-2\mathbb{E}\left[\mathbb{E}\left[\langle\nabla_{I_{n}}d_{y^{\delta}}(\lambda_{n}^{\delta}),\lambda_{n,I_{n}}^{\delta}-\lambda_{I_{n}}^{\dagger}\rangle\|{\mathcal{F}}_{n}\right]\right]$
		$\displaystyle\quad\,+t\mathbb{E}\left[\mathbb{E}\left[\\|\nabla_{I_{n}}d_{y^{\delta}}(\lambda_{n}^{\delta})\\|^{2}\|{\mathcal{F}}_{n}\right]\right]$
		$\displaystyle=\mathbb{E}[r_{n}]-2\mathbb{E}\left[\frac{1}{{p\choose b}}\sum_{I:\|I\|=b}\langle\nabla_{I}d_{y^{\delta}}(\lambda_{n}^{\delta}),\lambda_{n,I}^{\delta}-\lambda_{I}^{\dagger}\rangle\right]$
		$\displaystyle\quad\,+t\mathbb{E}\left[\frac{1}{{p\choose b}}\sum_{I:\|I\|=b}\\|\nabla_{I}d_{y^{\delta}}(\lambda_{n}^{\delta})\\|^{2}\right].$

Noting that

	$\displaystyle\sum_{I:\|I\|=b}\langle\nabla_{I}d_{y^{\delta}}(\lambda_{n}^{\delta}),\lambda_{n,I}^{\delta}-\lambda_{I}^{\dagger}\rangle$	$\displaystyle=\sum_{I:\|I\|=b}\sum_{i\in I}\langle\nabla_{i}d_{y^{\delta}}(\lambda_{n}^{\delta}),\lambda_{n,i}^{\delta}-\lambda_{,i}^{\dagger}\rangle$
		$\displaystyle=\sum_{i=1}^{p}\sum_{I:\|I\|=b\mbox{ and }i\in I}\langle\nabla_{i}d_{y^{\delta}}(\lambda_{n}^{\delta}),\lambda_{n,i}^{\delta}-\lambda_{,i}^{\dagger}\rangle$
		$\displaystyle={p-1\choose b-1}\sum_{i=1}^{p}\langle\nabla_{i}d_{y^{\delta}}(\lambda_{n}^{\delta}),\lambda_{n,i}^{\delta}-\lambda_{,i}^{\dagger}\rangle$
		$\displaystyle={p-1\choose b-1}\langle\nabla d_{y^{\delta}}(\lambda_{n}^{\delta}),\lambda_{n}^{\delta}-\lambda^{\dagger}\rangle.$

Therefore

\displaystyle\mathbb{E}[r_{n+1}]

\displaystyle=\mathbb{E}[r_{n}]+\frac{2b}{p}\mathbb{E}\left[\langle\nabla d_{y^{\delta}}(\lambda_{n}^{\delta}),\lambda^{\dagger}-\lambda_{n}^{\delta}\rangle\right]+\frac{t}{{p\choose b}}\sum_{I:|I|=b}\mathbb{E}\left[\|\nabla_{I}d_{y^{\delta}}(\lambda_{n}^{\delta})\|^{2}\right].

By the convexity of $d_{y^{\delta}}(\lambda)$ , the Cauchy-Schwarz inequality and the definition of $r_{n}$ , we can obtain

	$\displaystyle\langle\nabla d_{y^{\delta}}(\lambda_{n}^{\delta}),\lambda^{\dagger}-\lambda_{n}^{\delta}\rangle$	$\displaystyle\leq d_{y^{\delta}}(\lambda^{\dagger})-d_{y^{\delta}}(\lambda_{n}^{\delta})$
		$\displaystyle=d_{y}(\lambda^{\dagger})-d_{y}(\lambda_{n}^{\delta})+\langle\lambda_{n}^{\delta}-\lambda^{\dagger},y^{\delta}-y\rangle$
		$\displaystyle\leq d_{y}(\lambda^{\dagger})-d_{y}(\lambda_{n}^{\delta})+\delta\\|\lambda_{n}^{\delta}-\lambda^{\dagger}\\|$
		$\displaystyle=d_{y}(\lambda^{\dagger})-d_{y}(\lambda_{n}^{\delta})+\sqrt{t}\delta r_{n}^{1/2}.$

Therefore

	$\displaystyle\mathbb{E}[r_{n+1}]$	$\displaystyle\leq\mathbb{E}[r_{n}]+\frac{2b}{p}\mathbb{E}\left[d_{y}(\lambda^{\dagger})-d_{y}(\lambda_{n}^{\delta})\right]+\frac{2\sqrt{t}b}{p}\delta\mathbb{E}[r_{n}^{1/2}]$
		$\displaystyle\quad\,+\frac{t}{{p\choose b}}\sum_{I:\|I\|=b}\mathbb{E}\left[\\|\nabla_{I}d_{y^{\delta}}(\lambda_{n}^{\delta})\\|^{2}\right].$

Since $\mathbb{E}[r_{n}^{1/2}]\leq(\mathbb{E}[r_{n}])^{1/2}$ , we thus obtain the desired inequality. ∎

Lemma 4.5.

Consider Algorithm 3 with the step-sizes chosen as in Lemma 4.4, there holds

	$\displaystyle\mathbb{E}\left[d_{y}(\lambda_{n+1}^{\delta})-d_{y}(\lambda^{\dagger})+\frac{1}{2}c_{1}r_{n+1}\right]+\frac{bc_{1}}{p}\sum_{k=0}^{n}\mathbb{E}\left[d_{y}(\lambda_{k}^{\delta})-d_{y}(\lambda^{\dagger})\right]$
	$\displaystyle\leq M_{0}+\frac{t_{\max}b}{2c_{1}p}(n+1)\delta^{2}+\frac{t_{\max}^{1/2}bc_{1}}{p}\delta\sum_{k=0}^{n}\left(\mathbb{E}[r_{k}]\right)^{1/2}$

for all $n\geq 0$ , where $M_{0}:=\frac{1}{2}c_{1}r_{0}+\left(d_{y}(\lambda_{0})-d_{y}(\lambda^{\dagger})\right)$ .

Proof.

From Lemma 4.3 and Lemma 4.4 it follows that

	$\displaystyle\mathbb{E}\left[d_{y}(\lambda_{n+1}^{\delta})-d_{y}(\lambda^{\dagger})+\frac{1}{2}c_{1}r_{n+1}\right]$
	$\displaystyle\leq\mathbb{E}\left[d_{y}(\lambda_{n}^{\delta})-d_{y}(\lambda^{\dagger})+\frac{1}{2}c_{1}r_{n}\right]-\frac{bc_{1}}{p}\mathbb{E}\left[d_{y}(\lambda_{n}^{\delta})-d_{y}(\lambda^{\dagger})\right]$
	$\displaystyle\quad\,+\frac{t_{\max}b}{2c_{1}p}\delta^{2}+\frac{t_{\max}^{1/2}bc_{1}}{p}\delta\left(\mathbb{E}[r_{n}]\right)^{1/2}.$

Recursively using this inequality then completes the proof. ∎

In order to proceed further, we need an estimate on ${\mathbb{E}}[r_{n}]$ . We will use the following elementary result.

Lemma 4.6.

Let $\{a_{n}\}$ and $\{b_{n}\}$ be two sequences of nonnegative numbers such that

a_{n}^{2}\leq b_{n}^{2}+c\sum_{j=0}^{n-1}a_{j},\quad n=0,1,\cdots,

where $c\geq 0$ is a constant. If $\{b_{n}\}$ is non-decreasing, then

a_{n}\leq b_{n}+cn,\quad n=0,1,\cdots.

Proof.

We show the result by induction. The result is trivial for $n=0$ . Assume that the result is valid for all $0\leq n\leq k$ for some $k\geq 0$ . We show it is also true for $n=k+1$ . If $a_{k+1}\leq\max\{a_{0},\cdots,a_{k}\}$ , then $a_{k+1}\leq a_{l}$ for some $0\leq l\leq k$ . Thus, by the induction hypothesis and the monotonicity of $\{b_{n}\}$ we have

a_{k+1}\leq a_{l}\leq b_{l}+cl\leq b_{k+1}+c(k+1).

If $a_{k+1}>\max\{a_{0},\cdots,a_{k}\}$ , then

\displaystyle a_{k+1}^{2}\leq b_{k+1}^{2}+c\sum_{j=0}^{k}a_{j}\leq b_{k+1}^{2}+c(k+1)a_{k+1}

which implies that

	$\displaystyle\left(a_{k+1}-\frac{1}{2}c(k+1)\right)^{2}$	$\displaystyle=a_{k+1}^{2}-c(k+1)a_{k+1}+\frac{1}{4}c^{2}(k+1)^{2}$
		$\displaystyle\leq b_{k+1}^{2}+\frac{1}{4}c^{2}(k+1)^{2}$
		$\displaystyle\leq\left(b_{k+1}+\frac{1}{2}c(k+1)\right)^{2}.$

Taking square roots shows $a_{k+1}\leq b_{k+1}+c(k+1)$ again. ∎

Now we are ready to show the following result which together with Lemma 4.2 implies Theorem 4.1 immediately.

Theorem 4.7.

Let Assumption 1 hold and consider Algorithm 3 with the step-sizes chosen as in Theorem 4.1. If the sought solution $x^{\dagger}$ satisfies the source condition (4.1), then for all $n\geq 0$ there holds

\displaystyle\mathbb{E}\left[d_{y}(\lambda_{n}^{\delta})-d_{y}(\lambda^{\dagger})\right]\leq\frac{2M_{0}}{1+\frac{b}{p}c_{1}n}+\frac{5t_{\max}b}{2c_{1}p}n\delta^{2}

(4.8)

and consequently

\displaystyle\mathbb{E}\left[\|x_{n}^{\delta}-x^{\dagger}\|^{2}\right]\leq\frac{1}{\sigma}\left(\frac{2M_{0}}{1+\frac{b}{p}c_{1}n}+\frac{5t_{\max}b}{2c_{1}p}n\delta^{2}\right),

where $t_{\max}$ and $c_{1}$ are defined in Lemma 4.3.

Proof.

According to (4.6), it suffices to show (4.8). Since $\lambda^{\dagger}$ is a minimizer of $d_{y}(\lambda)$ over $Y$ , we have $d_{y}(\lambda_{n}^{\delta})-d_{y}(\lambda^{\dagger})\geq 0$ for all $n\geq 0$ . Thus, by noting that $\mathbb{E}[r_{0}]=r_{0}\leq 2M_{0}/c_{1}$ , it follows from Lemma 4.5 that

\displaystyle\mathbb{E}\left[r_{n}\right]

\displaystyle\leq\frac{2M_{0}}{c_{1}}+\frac{t_{\max}b}{c_{1}^{2}p}n\delta^{2}+\frac{2t_{\max}^{1/2}b}{p}\delta\sum_{k=0}^{n-1}\left(\mathbb{E}[r_{k}]\right)^{1/2}

for all $n\geq 0$ . Applying Lemma 4.6 and using the inequality $\sqrt{a+b}\leq\sqrt{a}+\sqrt{b}$ for any $a,b\geq 0$ , we can obtain

\displaystyle\left(\mathbb{E}[r_{n}]\right)^{1/2}\leq\sqrt{\frac{2M_{0}}{c_{1}}+\frac{t_{\max}b}{c_{1}^{2}p}n\delta^{2}}+\frac{2t_{\max}^{1/2}b}{p}n\delta\leq\sqrt{\frac{2M_{0}}{c_{1}}}+\left(\sqrt{\frac{bn}{c_{1}^{2}p}}+\frac{2bn}{p}\right)t_{\max}^{1/2}\delta.

Consequently

	$\displaystyle\sum_{k=0}^{n}\left(\mathbb{E}[r_{k}]\right)^{1/2}$	$\displaystyle\leq\sqrt{\frac{2M_{0}}{c_{1}}}(n+1)+\sum_{k=0}^{n}\left(\sqrt{\frac{bk}{c_{1}^{2}p}}+\frac{2bk}{p}\right)t_{\max}^{1/2}\delta$
		$\displaystyle\leq\sqrt{\frac{2M_{0}}{c_{1}}}(n+1)+\left(\frac{2b^{1/2}(n+1)^{3/2}}{3c_{1}p^{1/2}}+\frac{bn(n+1)}{p}\right)t_{\max}^{1/2}\delta.$

Combining this with the estimate in Lemma 4.5 we obtain

	$\displaystyle\mathbb{E}\left[d_{y}(\lambda_{n+1}^{\delta})-d_{y}(\lambda^{\dagger})\right]+\frac{bc_{1}}{p}\sum_{k=0}^{n}\mathbb{E}\left[d_{y}(\lambda_{k}^{\delta})-d_{y}(\lambda^{\dagger})\right]$
	$\displaystyle\leq M_{0}+\frac{t_{\max}b}{2c_{1}p}(n+1)\delta^{2}+\frac{b\sqrt{2t_{\max}c_{1}M_{0}}}{p}(n+1)\delta$
	$\displaystyle\quad\,+\left(\frac{2b^{3/2}(n+1)^{3/2}}{3p^{3/2}}+\frac{b^{2}c_{1}(n+1)^{2}}{p^{2}}\right)t_{\max}\delta^{2}.$

By using the inequality

\frac{b\sqrt{2t_{\max}c_{1}M_{0}}}{p}(n+1)\delta\leq M_{0}+\frac{b^{2}c_{1}(n+1)^{2}}{2p^{2}}t_{\max}\delta^{2}

we can further obtain

		$\displaystyle\mathbb{E}\left[d_{y}(\lambda_{n+1}^{\delta})-d_{y}(\lambda^{\dagger})\right]+\frac{bc_{1}}{p}\sum_{k=0}^{n}\mathbb{E}\left[d_{y}(\lambda_{k}^{\delta})-d_{y}(\lambda^{\dagger})\right]$
		$\displaystyle\leq 2M_{0}+\frac{t_{\max}b}{p}\left(\frac{1}{2c_{1}}+\frac{2b^{1/2}(n+1)^{1/2}}{3p^{1/2}}+\frac{3bc_{1}(n+1)}{2p}\right)(n+1)\delta^{2}.$		(4.9)

According to Lemma 4.3 we have for any $0\leq k\leq n$ that

\mathbb{E}\left[d_{y}(\lambda_{k}^{\delta})\right]\geq\mathbb{E}\left[d_{y}(\lambda_{n+1}^{\delta})\right]-\frac{t_{\max}b}{2c_{1}p}(n+1-k)\delta^{2}.

Therefore

	$\displaystyle\sum_{k=0}^{n}\mathbb{E}\left[d_{y}(\lambda_{k}^{\delta})-d_{y}(\lambda^{\dagger})\right]$
	$\displaystyle\geq(n+1)\mathbb{E}\left[d_{y}(\lambda_{n+1}^{\delta})-d_{y}(\lambda^{\dagger})\right]-\frac{t_{\max}b}{4c_{1}p}(n+1)(n+2)\delta^{2}.$

Combining this with (4) and noting that $n+2\leq 2(n+1)$ and

\frac{2b^{1/2}(n+1)^{1/2}}{3p^{1/2}}\leq\frac{2}{c_{1}}+\frac{bc_{1}(n+1)}{18p},

we obtain

	$\displaystyle\left(1+\frac{bc_{1}(n+1)}{p}\right)\mathbb{E}\left[d_{y}(\lambda_{n+1}^{\delta})-d_{y}(\lambda^{\dagger})\right]$
	$\displaystyle\leq 2M_{0}+\frac{t_{\max}b}{p}\left(\frac{1}{2c_{1}}+\frac{2b^{1/2}(n+1)^{1/2}}{3p^{1/2}}+\left(\frac{1}{2}+\frac{3c_{1}}{2}\right)\frac{b(n+1)}{p}\right)(n+1)\delta^{2}$
	$\displaystyle\leq 2M_{0}+\frac{t_{\max}b}{p}\left(\alpha+\frac{\beta b(n+1)}{p}\right)(n+1)\delta^{2},$

where

\alpha:=\frac{5}{2c_{1}}\quad\mbox{and}\quad\beta:=\frac{1}{2}+\frac{14}{9}c_{1}.

Note that $\beta/\alpha<c_{1}$ as $c_{1}<1$ . Dividing the both sides of the above equation by $1+bc_{1}(n+1)/p$ shows

\displaystyle\mathbb{E}\left[d_{y}(\lambda_{n+1}^{\delta})-d_{y}(\lambda^{\dagger})\right]

\displaystyle\leq\frac{2M_{0}}{1+\frac{b}{p}c_{1}(n+1)}+\frac{5t_{\max}b}{2c_{1}p}(n+1)\delta^{2}

which completes the proof. ∎

5. Numerical results

In this section we will present numerical simulations to test the performance of the stochastic mirror descent method for solving linear ill-posed problems in which the sought solutions have various special features that require to make particular choices of the regularization functional ${\mathcal{R}}$ and the solution space $X$ . Except Example 5.2, all the numerical simulations are performed on the linear ill-posed system

\displaystyle A_{i}x:=\int_{a}^{b}k(s_{i},t)x(t)dt=y_{i},\quad i=1,\cdots,p

(5.1)

obtained from linear integral equations of first kind on $[a,b]$ by sampling at $s_{i}\in[a,b]$ with $i=1,\cdots,p$ , where the kernel $k(\cdot,\cdot)$ is continuous on $[a,b]\times[a,b]$ and $s_{i}=a+(i-1)(b-a)/(p-1)$ for $i=1,\cdots,p$ .

Example 5.1.

Consider the linear system (1.1) where $X$ and $Y_{i}$ are all Hilbert spaces. In case the minimal norm solution is of interest, we may take ${\mathcal{R}}(x)=\|x\|^{2}/2$ in Algorithm 1, then the definition of $x_{n}^{\delta}$ shows $x_{n}^{\delta}=\xi_{n}^{\delta}$ . Consequently, Algorithm 1 becomes the form

\displaystyle x_{n+1}^{\delta}=x_{n}^{\delta}-t_{n}^{\delta}A_{I_{n}}^{*}(A_{I_{n}}x_{n}^{\delta}-y_{I_{n}}^{\delta}),

(5.2)

where $I_{n}\subset\{1,\cdots,p\}$ is a randomly selected subset with $|I_{n}|=b$ via the uniform distribution for a preassigned batch size $b$ . This is exactly the mini-batch stochastic gradient descent method for solving linear ill-posed problems in Hilbert spaces which has been studied recently in [17, 18, 25]. In these works the convergence analysis has been performed under diminishing step-sizes. Our work supplements the existing results by providing convergence results under new choices of step-sizes, in particular, for the method (5.2) with batch size $b=1$ , we obtain convergence and convergence rate results for the step-size

\displaystyle t_{n}^{\delta}=\frac{\mu_{0}}{\|A_{i_{n}}\|^{2}}\quad\mbox{ with }0<\mu_{0}<2

(5.3)

which obeys (s1) and (s2), where $i_{n}$ is the index chosen at the $n$ -th iteration step.

We now test the performance of the method (5.2) by considering the system (5.1) with $[a,b]=[-6,6]$ , $p=1000$ and $k(s,t):=\varphi(s-t)$ , where $\varphi(s)=(1+\cos(\pi s/3))\chi_{\{|s|<3\}}$ . We assume the sought solution is

x^{\dagger}(t)=\sin(\pi t/12)+\sin(\pi t/3)+\frac{1}{200}t^{2}(1-t).

Instead of the exact data $y_{i}:=A_{i}x^{\dagger}$ , we use the noisy data

\displaystyle y_{i}^{\delta}=y_{i}+\delta_{rel}|y_{i}|\varepsilon_{i},\quad i=1,\cdots,p,

(5.4)

where $\delta_{rel}$ is the relative noise level and $\varepsilon_{i}$ are random noises obeying the standard Gaussian distribution, i.e. $\varepsilon_{i}\sim N(0,1)$ . We execute the method (5.2) with the batch size $b=1$ and the initial guess $x_{0}^{\delta}=0$ together with the step-sizes given by (5.3); the integrals involved in the method are approximated by the trapezoidal rule based on the partition of $[-6,6]$ into $p-1$ subintervals of equal length. To illustrate the dependence of convergence on the magnitude of step-size, we consider the three distinct values $\mu_{0}=1.0,0.5$ and $0.1$ . We also use noisy data with three distinct relative noise levels $\delta_{rel}=10^{-1},10^{-2}$ and $10^{-3}$ . In Figure 1 we plot the corresponding reconstruction errors; the first row plots the relative mean square errors $\mathbb{E}[\|x_{n}^{\delta}-x^{\dagger}\|_{L^{2}}^{2}/\|x^{\dagger}\|_{L^{2}}^{2}]$ which are calculated approximately by the average of $100$ independent runs and the second row plots $\|x_{n}^{\delta}-x^{\dagger}\|_{L^{2}}^{2}/\|x^{\dagger}\|_{L^{2}}^{2}$ for a particular individual run. From these plots we can observe that, for each individual run, the iteration errors exhibit oscillations which are in particular dramatic when the noisy data have relative large noise levels. We can also observe that the method (5.2) demonstrates the semi-convergence property, i.e. the iterates converges to the sought solution at the beginning and, after a critical number of iterations, the iterates begin to diverge. This is typical for any iterative methods for solving ill-posed problems. Furthermore, for a fixed noise level, the semi-convergence occurs earlier if a larger step-size is used. We also note that, for noisy data with relatively large noise level, if a large step-size is used, the iterates can quickly produce a reconstruction result with minimal error and then diverge quickly from the sought solution; this makes it difficult to decide how to stop the iteration to produce a good reconstruction result. The use of small step-sizes has the advantage of suppressing oscillations and deferring semi-convergence. However, it slows down the convergence and hence makes the method time-consuming.

Refer to caption — Figure 1. Reconstruction errors of the stochastic gradient descent (5.2) using step-size (5.3) with various values of $\mu_{0}$ .

In order to efficiently suppress the oscillations and reduce the effect of semi-convergence, we next consider the method (5.2) using the step-size chosen by (s3) whose realization relies on the information of noise level. We now assume that the noisy data have the form (5.4), where $\varepsilon_{i}$ are noise uniformly distributed on $[-1,1]$ . Note that $\delta_{i}:=\delta_{rel}|y_{i}|$ are the noise levels with $|y_{i}^{\delta}-y_{i}|\leq\delta_{i}$ . Assuming the information on $\delta_{i}$ is known, the step-size chosen by (s3) takes the form

\displaystyle t_{n}^{\delta}=\left\{\begin{array}[]{lll}\mu_{0}/\|A_{i_{n}}\|^{2}&\mbox{ if }|A_{i_{n}}x_{n}^{\delta}-y_{i_{n}}^{\delta}|>\tau\delta_{i_{n}},\\ 0&\mbox{ otherwise}\end{array}\right.

(5.7)

which incorporates the spirit of the discrepancy principle. To illustrate the advantage of this choice of step-size, we compare the computed results with the ones obtained by the step-size chosen by (5.3). In Figure 2 we present the numerical results of reconstruction errors by the method (5.2) with batch size $b=1$ for various noise levels using the step-sizes chosen by (5.3) and (5.7) with $\mu_{0}=1$ and $\tau=1$ . where “DP” and “No DP” represent the results corresponding to the step-sizes chosen by (5.7) and (5.3) respectively. The first row in Figure 2 plots the mean square errors which are calculate approximately by the average of 100 independent runs and the second row plots the reconstruction errors for a typical individual run. The results show clearly that using the step-size (5.7) can significantly suppress the oscillations in reconstruction errors in particular when the data are corrupted by noise with relatively large noise level. We have performed extensive simulations which indicate that, due to the regularization effect of the discrepancy principle, using the step-size (5.7) has the tendency to decrease the iteration error as the number of iterations increases and has the ability to reduce the effect of semi-convergence.

Example 5.2.

We consider the linear system (1.1) in which $X$ and $Y_{i}$ are all Hilbert spaces and the sought solution satisfies the constraint $x\in\mathcal{C}$ , where $\mathcal{C}\subset X$ is a closed convex set. Finding the unique solution $x^{\dagger}$ of (1.1) in $\mathcal{C}$ with minimal norm can be stated as (1.2) with ${\mathcal{R}}(x):=\frac{1}{2}\|x\|^{2}+\iota_{\mathcal{C}}(x)$ , where $\iota_{\mathcal{C}}$ denotes the indicator function of $\mathcal{C}$ . Let $P_{\mathcal{C}}$ denote the metric projection of $X$ onto $\mathcal{C}$ . Then the mini-batch stochastic mirror descent method for determining $x^{\dagger}$ takes the form

\displaystyle x_{n}^{\delta}=P_{\mathcal{C}}(\xi_{n}^{\delta}),\qquad\xi_{n+1}^{\delta}=\xi_{n}^{\delta}-t_{n}^{\delta}A_{I_{n}}^{*}(A_{I_{n}}x_{n}^{\delta}-y_{I_{n}}^{\delta}).

(5.8)

In case $X=L^{2}(\mathscr{D})$ for some domain $\mathscr{D}\subset{\mathbb{R}}^{d}$ and $\mathcal{C}=\{x\in X:x\geq 0\mbox{ a.e. on }\mathscr{D}\}$ , the iteration scheme (5.8) becomes

\displaystyle x_{n}^{\delta}=\max\{\xi_{n}^{\delta},0\},\qquad\xi_{n+1}^{\delta}=\xi_{n}^{\delta}-t_{n}^{\delta}A_{I_{n}}^{*}(A_{I_{n}}x_{n}^{\delta}-y_{In}^{\delta})

(5.9)

with an initial guess $\xi_{0}^{\delta}=0$ , where $I_{n}$ is a randomly selected subset of $\{1,\cdots,p\}$ via the uniform distribution with $|I_{n}|=b$ for a preassigned batch size $b$ .

As a testing example, we consider the computed tomography which consists in determining the density of cross sections of a human body by measuring the attenuation of X-rays as they propagate through the biological tissues [27]. Mathematically, it requires to determine a function supported on a bounded domain from its line integrals. We consider here only the standard 2D parallel-beam tomography; tomography with other scan geometries can be considered similarly. We consider a full angle problem using $90$ projection angles evenly distributed between $1$ and $180$ degrees, with $367$ lines per projection. Assuming the sought image is discretized on a $256\times 256$ pixel grid, we may use the function paralleltomo in the MATLAB package AIR TOOLS [16] to discretize the problem. After deleting those rows with zero entries, it leads to a linear system $Ax=y$ , where $A$ is a matrix with size $29658\times 65536$ . This gives a problem of the form (1.1) with $p=29658$ , where $A_{i}$ corresponds to the $i$ -th row of $A$ and $y_{i}$ is the $i$ -th component of $y$ .

We assume that the true image is the Shepp-Logan phantom shown in figure 3 (a) discretized on a $256\times 256$ pixel grid with nonnegative pixel values. This phantom is widely used in evaluating tomographic reconstruction algorithms. Let $x^{\dagger}$ denote the vector formed by stacking all the columns of the true image and let $y=Ax^{\dagger}$ be the true data. By adding Gaussian noise to $y_{i}$ we get the noisy data $y_{i}^{\delta}$ of the form (5.4), where $\varepsilon_{i}\sim N(0,1)$ and $\delta_{rel}=0.01$ . Since the sought solution is nonnegative, we may reconstruct it by using the method (5.9). As a comparison, we also consider reconstructing the sought phantom by the method (5.2) for which the nonnegative constraint is not incorporated. For the both methods, we use the batch size $b=400$ and the step-size chosen by (s2) with $\mu_{0}=1$ . We execute the methods for 600 iterations and plot the reconstruction results in Figure 3 (b) and (c). In Figure 3 (d) we also plot the squares of the relative errors $\|x_{n}^{\delta}-x^{\dagger}\|^{2}/\|x^{\dagger}\|^{2}$ which indicate that incorporating the nonnegative constraint can produce more accurate reconstruction results. This example demonstrates that available a priori information on sought solutions should be incorporated into algorithm design to assist with better reconstruction results.

Example 5.3.

Let $\mathscr{D}\subset{\mathbb{R}}^{d}$ be a bounded domain. Consider the linear system (1.1), where $A_{i}:L^{1}(\mathscr{D})\to Y_{i}$ is a bounded linear operator and $Y_{i}$ is a Hilbert space for each $i$ . We assume the sought solution $x^{\dagger}$ is a probability density function, i.e. $x^{\dagger}\geq 0$ a.e. on $\mathscr{D}$ and $\int_{\mathscr{D}}x^{\dagger}=1$ . To determine such a solution, we take ${\mathcal{R}}(x):=f(x)+\iota_{\Delta}(x),$ where $\iota_{\Delta}$ denotes the indicator function of

\Delta:=\left\{x\in L_{+}^{1}(\mathscr{D}):\int_{\mathscr{D}}x^{\dagger}=1\right\}

and

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}),\\ +\infty,&\mbox{ otherwise}\end{array}\right.

is the negative of the Boltzmann-Shannon entropy. Here $L_{+}^{1}(\mathscr{D}):=\{x\in L^{1}(\mathscr{D}):x\geq 0\mbox{ a.e. on }\mathscr{D}\}$ . The Boltzmann-Shannon entropy has been used in Tikhonov regularization as a regularization functional to determine nonnegative solutions, see [9, 11] for instance. According to [3, 9], ${\mathcal{R}}$ is strongly convex on $L^{1}(\mathscr{D})$ with modulus of convexity $\sigma=1/2$ . By the Karush-Kuhn-Tucker theory, for any $\xi\in L^{\infty}(\mathscr{D})$ the unique minimizer of

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

is given by $\hat{x}=e^{\xi}/\int_{\mathscr{D}}e^{\xi}$ . Therefore the corresponding stochastic mirror descent method takes the form

\displaystyle x_{n}^{\delta}=\frac{1}{\int_{\mathscr{D}}e^{\xi_{n}^{\delta}}}e^{\xi_{n}^{\delta}},\qquad\xi_{n+1}^{\delta}=\xi_{n}^{\delta}-t_{n}^{\delta}A_{I_{n}}^{*}(A_{I_{n}}x_{n}^{\delta}-y_{I_{n}}^{\delta}),

(5.10)

where $I_{n}\subset\{1,\cdots,p\}$ is a randomly selected subset via the uniform distribution with a preassigned batch size $b$ and $t_{n}^{\delta}\geq 0$ is the step-size.

For numerical simulations we consider the linear system (5.1) with $[a,b]=[0,1]$ , $p=1000$ and $k(s,t):=4e^{-(s-t)^{2}/0.0064}$ . We assume the sought solution is

x^{\dagger}(t):=c\left(e^{-60(t-0.3)^{2}}+0.3e^{-40(t-0.8)^{2}}\right),

where $c>0$ is a constant to ensure that $\int_{0}^{1}x^{\dagger}(t)dt=1$ so that $x^{\dagger}$ is a probability density function. By adding random noise to the exact data $y_{i}:=A_{i}x^{\dagger}$ we get the noisy data $y_{i}^{\delta}$ ; we then use these noisy data and the method (5.10) with batch size $b=1$ to reconstruct the sought solution $x^{\dagger}$ ; the integrals involved in the method are approximated by the trapezoidal rule based on the partition of $[0,1]$ into $p-1$ subintervals of equal length.

We first test the performance of the method (5.10) using the noisy data given by (5.4) corrupted by Gaussian noise, where $\varepsilon_{i}\sim N(0,1)$ and $\delta_{rel}>0$ is the relative noise level. We execute the method (5.10) using the batch size $b=1$ , the initial guess $\xi_{0}^{\delta}=0$ and the step-size $t_{n}^{\delta}$ chosen by (s2) with three distinct values $\mu_{0}=2.0,1.0$ and $0.4$ . In Figure 4 we plot the reconstruction errors for three distinct relative noise levels $\delta_{rel}=0.5,0.1$ and $0.01$ . The first row plots the mean square errors $\mathbb{E}[\|x_{n}^{\delta}-x^{\dagger}\|_{L^{1}}^{2}]$ which are calculated approximately by the average of 100 independent runs. The second row plots $\|x_{n}^{\delta}-x^{\dagger}\|_{L^{1}}^{2}$ for a particular individual run. These plots demonstrate that the method (5.10) admits the semi-convergence property and the iterates exhibit dramatic oscillations. Furthermore, using large step-sizes can make convergence fast at the beginning but the iterates can diverge quickly; while using small step-sizes can suppress the oscillations and delay the occurrence of semi-convergence, it however makes the method converge slowly.

In order to efficiently remove oscillations and semi-convergence, we next consider the method (5.10) using the step-size chosen by (s3) whose realization relies on the information of noise level. We consider noisy data of the form (5.4), where $\varepsilon_{i}$ are uniform noise distributed on $[-1,1]$ . Note that $\delta_{i}:=\delta_{rel}|y_{i}|$ are the noise levels such that $|y_{i}^{\delta}-y_{i}|\leq\delta_{i}$ for all $i$ ; we assume the information on $\delta_{i}$ is known. To illustrate the advantage of the step-size chosen by (s3), we compare the computed results with the ones obtained by the step-size chosen by (s2). In Figure 5 we plot the reconstruction errors by the method (5.10) for various relative noise levels using the step-sizes chosen by (s2) and (s3) with $\mu_{0}=2$ and $\tau=1$ , where “DP” and “No DP” represent the results corresponding to the step-size chosen by (s3) and (s2) respectively. The first row in Figure 5 plots the mean square errors calculated approximately by the average of 100 independent runs and the second row plots the reconstruction errors for a particular individual run. The results show clearly that using step-size by (s3) can significantly suppress the oscillations in iterates and relieve the method from semi-convergence.

Example 5.4.

Let $\mathscr{D}\subset{\mathbb{R}}^{d}$ be a bounded domain. We consider reconstructing a sparse solution in the linear system (1.1), where $A_{i}:L^{2}(\mathscr{D})\to Y_{i}$ is a bounded linear operator and $Y_{i}$ is a Hilbert space for each $i$ . To determine such a solution, we take ${\mathcal{R}}(x):=\beta\|x\|_{L^{1}(\mathscr{D})}+\frac{1}{2}\|x\|_{L^{2}(\mathscr{D})}^{2}$ , where $\beta>0$ is a sufficiently large number whose choice reflects the role of the term $\|x\|_{L^{1}(\mathscr{D})}$ in the reconstruction of sparse solutions. Note that for any $\xi\in L^{2}(\mathscr{D})$ we have

\arg\min_{x\in L^{2}(\mathscr{D})}\left\{{\mathcal{R}}(x)-\langle\xi,x\rangle\right\}=\mbox{sign}(\xi)\max\{|\xi|-\beta,0\}.

Therefore the corresponding stochastic mirror descent method takes the form

\displaystyle x_{n}^{\delta}=\mbox{sign}(\xi_{n}^{\delta})\max\{|\xi_{n}^{\delta}|-\beta,0\},\quad\xi_{n+1}^{\delta}=\xi_{n}^{\delta}-t_{n}^{\delta}A_{I_{n}}^{*}(A_{I_{n}}x_{n}^{\delta}-y_{I_{n}}^{\delta}),

(5.11)

where $I_{n}\subset\{1,\cdots,p\}$ is a randomly selected subset satisfying $|I_{n}|=b$ via the uniform distribution with the preassigned batch size $b$ and $t_{n}^{\delta}\geq 0$ is the step-size. We remark that, when $|I_{n}|=1$ and $I_{n}=n\ (\mbox{mod }p)+1$ , the corresponding method is the sparse-Kaczmarz method which was first proposed in [23] to reconstruct sparse solutions of ill-posed inverse problems. For finite-dimensional linear systems with well-conditioned matrices, a randomized sparse-Kaczmarz method of the form (5.11) was introduced in [37] in which $I_{n}$ is chosen with a probability proportional to $\|A_{I_{n}}\|^{2}$ and a linear convergence rate is derived. The theoretical result in [37] however is not applicable to ill-posed problems because the analysis depends heavily on the finite-dimensionality of the underlying spaces and the well-conditioning property of the matrices.

For numerical simulations we consider the linear system (5.1) with $[a,b]=[0,1]$ and $k(s,t):=(0.1^{2}+(s-t)^{2})^{-3/2}$ . We assume the sought solution is

x^{\dagger}(t):=\chi_{\{0.19\leq t\leq 0.22\}}-\chi_{\{0.50\leq t\leq 0.52\}}+0.5\chi_{\{0.78\leq t\leq 0.80\}}

which is sparse on $[0,1]$ . By adding random noise to the exact data $y_{i}:=A_{i}x^{\dagger}$ we get the noisy data $y_{i}^{\delta}$ ; we then use these noisy data and the method (5.11) with batch size $b=1$ to reconstruct the sought solution $x^{\dagger}$ . In our numerical computation we take $p=1000$ and $\beta=80$ and the integrals involved in the method are approximated by the trapezoidal rule based on the partition of $[0,1]$ into $p-1$ subintervals of equal length.

We first test the performance of the method (5.11) using the noisy data given by (5.4) corrupted by Gaussian noise, where $\varepsilon_{i}\sim N(0,1)$ . We execute the method (5.11) using the batch size $b=1$ , the initial guess $\xi_{0}^{\delta}=0$ and the step-size $t_{n}^{\delta}$ chosen by (s2) with three distinct values $\mu_{0}=2.0,1.0$ and $0.4$ . In Figure 6 we plot the reconstruction errors for three distinct relative noise levels $\delta_{rel}=0.5,0.1$ and $0.01$ . The first row plots the relative mean square errors $\mathbb{E}[\|x_{n}^{\delta}-x^{\dagger}\|_{L^{2}}^{2}/\|x^{\dagger}\|_{L^{2}}^{2}]$ calculated approximately by the average of 100 independent runs, and the second row plots the relative errors $\|x_{n}^{\delta}-x^{\dagger}\|_{L^{2}}^{2}/\|x^{\dagger}\|_{L^{2}}^{2}$ for a particular individual run. From these plots we can observe the oscillations in iterates, the semi-convergence of the method, and the influence of the magnitude of step-sizes.

We next consider the method (5.11) using the step-size chosen by (s3). We use noisy data of the form (5.4) with $\varepsilon_{i}$ being uniformly distributed over $[-1,1]$ and assume that the noise levels $\delta_{i}=\delta_{rel}|y_{i}|$ are known. In Figure 7 we plot the reconstruction errors of the method (5.11) using the step-sizes chosen by (s3) with $\mu_{0}=2$ and $\tau=1.01$ which are labelled as “DP”; as comparisons we also plot the corresponding results using step-sizes chosen by (s2) with $\mu_{0}=2$ which are labelled as “No DP”. Sharply contrast to (s2), these results illustrate that using step-size by (s3) can significantly suppress the oscillations in iterates and reduce the effect of semi-convergence.

Example 5.5.

In this final example we consider using the stochastic mirror descent method to reconstruct piecewise constant solutions. We consider again the linear system (5.1) with $[a,b]=[0,1]$ and $k(s,t)=4e^{-(s-t)^{2}/0.0064}$ and assume that the sought solution $x^{\dagger}$ is piecewise constant. By dividing $[0,1]$ into $p-1$ subintervals of equal length and approximating integrals by the trapezoidal rule, we have a discrete ill-posed system ${\bf A}_{i}{\bf x}=y_{i}$ , where ${\bf A}_{i}$ is a row vector for each $i$ . To reconstruct a piecewise constant solution, we use the model

\displaystyle\min\left\{{\mathcal{R}}({\bf x}):=\beta\|{\bf D}{\bf x}\|_{1}+\frac{1}{2}\|{\bf D}{\bf x}\|_{2}^{2}+\frac{1}{2}\|{\bf x}\|_{2}^{2}:{\bf A}_{i}{\bf x}=y_{i},i=1,\cdots,p\right\},

(5.12)

where ${\bf D}$ denotes the discrete gradient operator and $\beta$ is a large positive number. If we apply the stochastic mirror descent method to solve (5.12) directly, we need to solve a minimization problem related to ${\mathcal{R}}$ to obtain ${\bf x}_{n}^{\delta}$ at each iteration step. This can make the algorithm time-consuming since those minimization problems can not be solved explicitly. To circumvent this difficulty, we introduce ${\bf z}:={\bf D}{\bf x}$ and rephrase (5.12) as

\displaystyle\min\left\{\beta\|{\bf z}\|_{1}+\frac{1}{2}\|{\bf z}\|_{2}^{2}+\frac{1}{2}\|{\bf x}\|_{2}^{2}:{\bf B}_{i}\begin{pmatrix}{\bf x}\\ {\bf z}\end{pmatrix}=\begin{pmatrix}y_{i}\\ {\bf 0}\end{pmatrix},i=1,\cdots,p\right\},

(5.13)

where ${\bf B}_{i}:=\begin{pmatrix}{\bf A}_{i}&{\bf 0}\\ {\bf D}&-{\bf I}\end{pmatrix}$ . Assuming, instead of $y_{i}$ , we have the noisy data $y_{i}^{\delta}$ . By applying the stochastic mirror descent method with batch size $b=1$ to (5.13) we can obtain

	$\displaystyle\begin{pmatrix}{\bf x}_{n}^{\delta}\\[4.30554pt] {\bf z}_{n}^{\delta}\end{pmatrix}=\arg\min_{{\bf x},{\bf z}}\left\{\beta\\|{\bf z}\\|_{1}+\frac{1}{2}\\|{\bf z}\\|_{2}^{2}+\frac{1}{2}\\|x\\|_{2}^{2}-\langle\xi_{n}^{\delta},{\bf x}\rangle-\langle\eta_{n}^{\delta},{\bf z}\rangle\right\},$
	$\displaystyle\begin{pmatrix}\xi_{n+1}^{\delta}\\[4.30554pt] \eta_{n+1}^{\delta}\end{pmatrix}=\begin{pmatrix}\xi_{n}^{\delta}\\[4.30554pt] \eta_{n}^{\delta}\end{pmatrix}-t_{n}^{\delta}{\bf B}_{i_{n}}^{T}\left({\bf B}_{i_{n}}\begin{pmatrix}{\bf x}_{n}^{\delta}\\[4.30554pt] {\bf z}_{n}^{\delta}\end{pmatrix}-\begin{pmatrix}y_{i_{n}}^{\delta}\\[4.30554pt] {\bf 0}\end{pmatrix}\right).$

By calculating ${\bf x}_{n}^{\delta},{\bf z}_{n}^{\delta}$ and noting ${\bf x}_{n}^{\delta}=\xi_{n}^{\delta}$ , we therefore obtain the following iteration scheme

\displaystyle\begin{split}{\bf z}_{n}^{\delta}&=\mbox{sign}(\eta_{n}^{\delta})\max\{|\eta_{n}^{\delta}|-\beta,0\},\\ {\bf x}_{n+1}^{\delta}&={\bf x}_{n}^{\delta}-t_{n}^{\delta}\left({\bf A}_{i_{n}}^{T}({\bf A}_{i_{n}}{\bf x}_{n}^{\delta}-y_{i_{n}}^{\delta})+{\bf D}^{T}({\bf D}{\bf x}_{n}^{\delta}-{\bf z}_{n}^{\delta})\right),\\ \eta_{n+1}^{\delta}&=\eta_{n}^{\delta}-t_{n}^{\delta}({\bf z}_{n}^{\delta}-{\bf D}{\bf x}_{n}^{\delta})\end{split}

(5.14)

with the initial guess $\eta_{0}^{\delta}={\bf 0}$ and ${\bf x}_{0}^{\delta}={\bf 0}$ , where $i_{n}\in\{1,\cdots,p\}$ is a randomly selected index via the uniform distribution and $t_{n}^{\delta}\geq 0$ denotes the step-size.

Assuming the sought solution ${\bf x}^{\dagger}$ is piecewise constant whose graph is plotted in Figure 8. By adding noise to the exact data $y_{i}:={\bf A}_{i}{\bf x}^{\dagger}$ , we produce the noisy data $y_{i}^{\delta}$ . We then use these noisy data in the method (5.14) to reconstruct ${\bf x}^{\dagger}$ . In the numerical computation we use $p=1000$ and $\beta=400$ . For the method (5.14) using step-size chosen by (s2), i.e.

\displaystyle t_{n}^{\delta}=\frac{\mu_{0}(|{\bf A}_{i_{n}}{\bf x}_{n}^{\delta}-y_{i_{n}}^{\delta}|^{2}+\|{\bf D}{\bf x}_{n}^{\delta}-{\bf z}_{n}^{\delta}\|_{2}^{2})}{\|{\bf A}_{i_{n}}^{T}({\bf A}_{i_{n}}{\bf x}_{n}^{\delta}-y_{i_{n}}^{\delta})\|_{2}^{2}+\|{\bf D}^{T}({\bf D}{\bf x}_{n}^{\delta}-{\bf z}_{n}^{\delta})\|_{2}^{2}}

(5.15)

we have performed the numerical computation for several distinct values of $\mu_{0}$ and observed semi-convergence property, oscillations of iterates and the effect of the magnitude of the step-sizes. Since these observations are very similar to the previous examples, we will not report them here. Instead we will focus on the computational effect of the step-size chosen by (s3). As before, the noisy data are generated by (5.4) with $\varepsilon_{i}$ being uniform noise on $[-1,1]$ and we assume the noise levels $\delta_{i}:=\delta_{rel}|y_{i}|$ are known. The step-size chosen by (s3) then takes the formula

\displaystyle t_{n}^{\delta}=\left\{\begin{array}[]{lll}\frac{\mu_{0}(|{\bf A}_{i_{n}}{\bf x}_{n}^{\delta}-y_{i_{n}}^{\delta}|^{2}+\|{\bf D}{\bf x}_{n}^{\delta}-{\bf z}_{n}^{\delta}\|_{2}^{2})}{\|{\bf A}_{i_{n}}^{T}({\bf A}_{i_{n}}{\bf x}_{n}^{\delta}-y_{i_{n}}^{\delta})\|_{2}^{2}+\|{\bf D}^{T}({\bf D}{\bf x}_{n}^{\delta}-{\bf z}_{n}^{\delta})\|_{2}^{2}}&\mbox{ if }|{\bf A}_{i_{n}}{\bf x}_{n}^{\delta}-y_{i_{n}}^{\delta}|>\tau\delta_{i_{n}},\\ 0&\mbox{ otherwise}.\end{array}\right.

(5.18)

In Figure 8 we plot the reconstruction results by the method (5.14) after $5\times 10^{5}$ iterations using noise data for three distinct relative noise levels $\delta_{rel}=0.5,0.1$ and $0.01$ , where “DP” and “No DP” represent the reconstruction results by the step-size chosen by (5.18) and (5.15) respectively; we use $\mu_{0}=1.0$ and $\tau=1.0$ . These results demonstrate that the method (5.14) can capture the piecewise constant feature very well. In order to give further comparison on the reconstruction results by the step-sizes (5.15) and (5.18), we present in Figure 9 the reconstruction errors, the first row plots approximations of the relative mean square errors $\mathbb{E}[\|{\bf x}_{n}^{\delta}-{\bf x}^{\dagger}\|_{2}^{2}/\|{\bf x}^{\dagger}\|_{2}^{2}]$ by an average of 100 independent runs and the second row plots the reconstruction errors $\|{\bf x}_{n}^{\delta}-{\bf x}^{\dagger}\|_{2}^{2}/\|{\bf x}^{\dagger}\|_{2}^{2}$ for a typical individual run. The results show clearly that using the step-size (5.18) can significantly suppress the oscillations in the iterates and reduce the effect of semi-convergence of the method.

Acknowledgement

Q. Jin would like to thank Peter Mathé from Weierstrass Institute for discussions on stochastic gradient descent. The work of Q. Jin is partially supported by the Future Fellowship of the Australian Research Council (FT170100231). The work of X. Lu is partially supported by the National Key Research and Development Program of China (No. 2020YFA0714200) and the National Science Foundation of China (No. 11871385).

References

[1] A. Beck and M. Teboulle, Mirror descent and nonlinear projected subgradient methods for convex optimization, Oper. Res. Lett., 31 (2003), no. 3, 167–175.
[2] R. Bhattacharya and E. C. Waymire, A basic course in probability theory, Second edition. Universitext. Springer, Cham, 2016.
[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] L. Bottou, F. E. Curtis and J. Nocedal, Optimization methods for large-scale machine learning, SIAM Rev., 60 (2018), no. 2, 223–311.
[6] S. Bubeck, Convex Optimization: Algorithms and Complexity, Foundations and Trends in Machine Learning: 8 (2015), no. 3-4, pp 231–357.
[7] M. Burger and S. Osher, Convergence rates of convex variational regularization, Inverse Problems 20 (2004), no. 5, 1411–1421.
[8] J. C. Duchi, Introductory lectures on stochastic convex optimization, In The Mathematics of Data, IAS/Park City Mathematics Series. American Mathematical Society, 2018.
[9] P. P. B. Eggermont, Maximum entropy regularization for Fredholm integral equations of the first kind, SIAM J. Math. Anal., 24 (1993), no. 6, 1557–1576.
[10] H. W. Engl, M. Hanke and A. Neubauer, Regularization of Inverse Problems, Kluwer, Dordrecht, 1996.
[11] H. W. Engl and G. Landl, Convergence rates for maximum entropy regularization, SIAM J. Numer. Anal., 30 (1993), 1509–1536.
[12] 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.
[13] I. Goodfellow, Y. Bengio and A. Courville, Deep Learning, Adaptive Computation and Machine Learning series, Cambridge, United States, 2017.
[14] C. W. Groetsch, The theory of Tikhonov regularization for Fredholm equations of the first kind, Research Notes in Mathematics, 105. Pitman (Advanced Publishing Program), Boston, MA, 1984.
[15] M. Haltmaier, A. Leitao and O. Scherzer, Regularization of systems of nonlinear ill-posed equations: I. Convergence Analysis, Inverse Problems and Imaging, 1 (2007), no. 2, 289–298.
[16] P. C. Hansen and M. Saxild-Hansen, AIR Tools–a MATLAB package of algebraic iterative reconstruction methods, J. Comput. Appl. Math. 236 (2012), 2167–2178.
[17] Y. Jiao, B. Jin and X. Lu, Preasymptotic convergence of randomized Kaczmarz method, Inverse Problems, 33 (2017), 125012, 21 pp.
[18] B. Jin and X. Lu, On the regularizing property of stochastic gradient descent, Inverse Problems, 35 (2019), no. 1, 015004, 27 pp.
[19] B. Jin, Z. Zhou and J. Zou, On the convergence of stochastic gradient descent for nonlinear ill-posed problems, SIAM J. Optim., 30 (2020), no. 2, 1421–1450.
[20] Q. Jin, Landweber-Kaczmarz method in Banach spaces with inexact inner solvers, Inverse Problems, 32 (2016), no. 10, 104005, 26 pp.
[21] Q. Jin, Convergence rates of a dual gradient method for linear ill-posed problems, Numer. Math., accepted for publication, 2022.
[22] Q. Jin and X. Lu, A fast nonstationary iterative method with convex penalty for inverse problems in Hilbert spaces, Inverse Problems, 30 (2014), 04501, 21 pp.
[23] 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.
[24] B. Kaltenbacher, A. Neubauer and O. Scherzer, Iterative Regularization Methods for Nonlinear Ill-Posed Problems, Berlin: de Gruyter, 2008.
[25] S. Lu and P. Mathé, Stochastic gradient descent for linear inverse problems in Hilbert spaces, Math. Comput., 91 (2022), 1763–1788.
[26] Z. Lu and L. Xiao, On the complexity analysis of randomized block-coordinate descent methods, Math. Program., 152 (2015), no. 1–2, Ser. A, 615–642.
[27] F. Natterer,The Mathematics of Computerized Tomography, SIAM, Philadelphia, 2001.
[28] A. Nedic and S. Lee, On stochastic subgradient mirror-descent algorithm with weighted averaging, SIAM J. Optim., 24 (2014), pp. 84–107.
[29] A. Nemirovski, A. Juditsky, G. Lan and A. Shapiro, Robust stochastic approximation approach to stochastic programming, SIAM J. Optim., 19 (2009), 1574–1609.
[30] S. Nemirovski and D. B. Yudin, Problem Complexity and Method Efficiency in Optimization, Wiley, New York, 1983.
[31] Y. Nesterov, Primal-dual subgradient methods for convex problems, Math. Program., Ser. B 120 (2009), 221–259.
[32] Y. Nesterov, Efficiency of coordinate descent methods on huge-scale optimization problems, SIAM J. Optim., 22 (2012), no. 2, 341–362.
[33] J. C. Rabelo, Y. F. Saporito and A. Leitão, On stochastic Kaczmarz type methods for solving large scale systems of ill-posed equations, Inverse Problems, 38 (2022), 025003.
[34] R. T. Rockafellar, Convex Analysis, Princeton University Press, Princeton, NJ, 1970.
[35] 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.
[36] P. Richtárik and M. Takác̆, Iteration complexity of randomized block-coordinate descent methods for minimizing a composite function, Math. Program., 144 (2014), no. 1–2, 1–38.
[37] F. Schöpfer and D. A. Lorenz, Linear convergence of the randomized sparse Kaczmarz method, Math. Program., 173 (2019), no. 1-2, Ser. A, 509–536.
[38] C. Zălinscu, Convex Analysis in General Vector Spaces, World Scientific Publishing Co., Inc., River Edge, New Jersey, 2002.
[39] Z. Zhou, P. Mertikopoulos, N. Bambos, S. P. Boyd and P. W. Glynn, On the convergence of mirror descent beyond stochastic convex programming, SIAM J. Optim., 30 (2020), no. 1, 687–716.

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

	$\displaystyle\Delta_{n+1}^{\delta}-\Delta_{n}^{\delta}$	$\displaystyle\leq\frac{1}{4\sigma}\left(t_{n}^{\delta}\right)^{2}\\|A_{I_{n}}^{*}(A_{I_{n}}x_{n}^{\delta}-y_{I_{n}}^{\delta})\\|^{2}-t_{n}^{\delta}\langle A_{I_{n}}x_{n}^{\delta}-y_{I_{n}}^{\delta},A_{I_{n}}x_{n}^{\delta}-y_{I_{n}}\rangle$
		$\displaystyle\leq\frac{1}{4\sigma}\left(t_{n}^{\delta}\right)^{2}\\|A_{I_{n}}^{*}(A_{I_{n}}x_{n}^{\delta}-y_{I_{n}}^{\delta})\\|^{2}-t_{n}^{\delta}\\|A_{I_{n}}x_{n}^{\delta}-y_{I_{n}}^{\delta}\\|^{2}$
		$\displaystyle\quad\,+\delta_{I_{n}}t_{n}^{\delta}\\|A_{I_{n}}x_{n}^{\delta}-y_{I_{n}}^{\delta}\\|.$

	$\displaystyle\xi_{n+1}^{\delta}$	$\displaystyle=A^{}\lambda_{n+1}^{\delta}=A_{I_{n}^{c}}^{}\lambda_{n+1,I_{n}^{c}}^{\delta}+A_{I_{n}}^{*}\lambda_{n+1,I_{n}}^{\delta}$
		$\displaystyle=A_{I_{n}^{c}}^{}\lambda_{n,I_{n}^{c}}^{\delta}+A_{I_{n}}^{}\left(\lambda_{n,I_{n}}^{\delta}-t_{I_{n}}(A_{I_{n}}x_{n}^{\delta}-y_{I_{n}}^{\delta})\right)$
		$\displaystyle=\left(A_{I_{n}^{c}}^{}\lambda_{n,I_{n}^{c}}^{\delta}+A_{I_{n}}^{}\lambda_{n,I_{n}}^{\delta}\right)-t_{I_{n}}A_{I_{n}}^{*}(A_{I_{n}}x_{n}^{\delta}-y_{I_{n}}^{\delta})$
		$\displaystyle=\xi_{n}^{\delta}-t_{I_{n}}A_{I_{n}}^{*}(A_{I_{n}}x_{n}^{\delta}-y_{I_{n}}^{\delta}).$

	$\displaystyle d_{y}(\lambda_{n+1}^{\delta})-d_{y}(\lambda_{n}^{\delta})$
	$\displaystyle\leq t_{I_{n}}\langle\nabla_{I_{n}}d_{y^{\delta}}(\lambda_{n}^{\delta}),y_{I_{n}}-y_{I_{n}}^{\delta}\rangle-\left(1-\frac{L_{I_{n}}t_{I_{n}}}{2}\right)t_{I_{n}}\\|\nabla_{I_{n}}d_{y^{\delta}}(\lambda_{n}^{\delta})\\|^{2}$
	$\displaystyle\leq t_{I_{n}}\delta_{I_{n}}\\|\nabla_{I_{n}}d_{y^{\delta}}(\lambda_{n}^{\delta})\\|-c_{1}t_{I_{n}}\\|\nabla_{I_{n}}d_{y^{\delta}}(\lambda_{n}^{\delta})\\|^{2}$
	$\displaystyle\leq\frac{1}{2c_{1}}t_{I_{n}}\delta_{I_{n}}^{2}-\frac{1}{2}c_{1}t_{I_{n}}\\|\nabla_{I_{n}}d_{y^{\delta}}(\lambda_{n}^{\delta})\\|^{2}.$

	$\displaystyle r_{n+1}$	$\displaystyle=\frac{1}{t}\\|\lambda_{n+1,I_{n}^{c}}^{\delta}-\lambda_{I_{n}^{c}}^{\dagger}\\|^{2}+\frac{1}{t}\\|\lambda_{n+1,I_{n}}^{\delta}-\lambda_{I_{n}}^{\dagger}\\|^{2}$
		$\displaystyle=\frac{1}{t}\\|\lambda_{n,I_{n}^{c}}^{\delta}-\lambda_{I_{n}^{c}}^{\dagger}\\|^{2}+\frac{1}{t}\\|\lambda_{n,I_{n}}^{\delta}-\lambda_{I_{n}}^{\dagger}-t\nabla_{I_{n}}d_{y^{\delta}}(\lambda_{n}^{\delta})\\|^{2}$
		$\displaystyle=\frac{1}{t}\\|\lambda_{n,I_{n}^{c}}^{\delta}-\lambda_{I_{n}^{c}}^{\dagger}\\|^{2}+\frac{1}{t}\\|\lambda_{n,I_{n}}^{\delta}-\lambda_{I_{n}}^{\dagger}\\|^{2}-2\langle\nabla_{I_{n}}d_{y^{\delta}}(\lambda_{n}^{\delta}),\lambda_{n,I_{n}}^{\delta}-\lambda_{I_{n}}^{\dagger}\rangle$
		$\displaystyle\quad\,+t\\|\nabla_{I_{n}}d_{y^{\delta}}(\lambda_{n}^{\delta})\\|^{2}$
		$\displaystyle=r_{n}-2\langle\nabla_{I_{n}}d_{y^{\delta}}(\lambda_{n}^{\delta}),\lambda_{n,I_{n}}^{\delta}-\lambda_{I_{n}}^{\dagger}\rangle+t\\|\nabla_{I_{n}}d_{y^{\delta}}(\lambda_{n}^{\delta})\\|^{2}.$