# On the Convergence of Stochastic Gradient Descent for Linear Inverse
Problems in Banach Spaces

B. Jin and Z. Kereta

###### Abstract

In this work we consider stochastic gradient descent (SGD) for solving linear
inverse problems in Banach spaces. SGD and its variants have been established
as one of the most successful optimisation methods in machine learning,
imaging and signal processing, etc. At each iteration SGD uses a single datum,
or a small subset of data, resulting in highly scalable methods that are very
attractive for large-scale inverse problems. Nonetheless, the theoretical
analysis of SGD-based approaches for inverse problems has thus far been
largely limited to Euclidean and Hilbert spaces. In this work we present a
novel convergence analysis of SGD for linear inverse problems in general
Banach spaces: we show the almost sure convergence of the iterates to the
minimum norm solution and establish the regularising property for suitable a
priori stopping criteria. Numerical results are also presented to illustrate
features of the approach.

##  1 Introduction

This work considers (stochastic) iterative solutions for linear operator
equations of the form

$\displaystyle{\mathbf{A}}{\boldsymbol{\mathsf{x}}}={\boldsymbol{\mathsf{y}}}$
(1)

where ${\mathbf{A}}:{\cal X}\rightarrow{\cal Y}$ is a bounded linear operator
between Banach spaces ${\cal X}$ and ${\cal Y}$ (equipped with the norms
$\|\cdot\|_{{\cal X}}$ and $\|\cdot\|_{{\cal Y}}$, respectively), and
${\boldsymbol{\mathsf{y}}}\in\operatorname{range}{({\mathbf{A}})}$ is the
exact data. In practice, we only have access to noisy data
${\boldsymbol{\mathsf{y}}}^{\delta}={\boldsymbol{\mathsf{y}}}+{\boldsymbol{\xi}}$,
where $\boldsymbol{\xi}$ denotes the measurement noise with a noise level
$\delta\geq 0$ such that
$\|{\boldsymbol{\mathsf{y}}}^{\delta}-{\boldsymbol{\mathsf{y}}}\|_{{\cal
Y}}{\leq}\delta$. Linear inverse problems arise naturally in many applications
in science and engineering, and also form the basis for studying nonlinear
inverse problems. Hence, design and analysis of stable reconstruction methods
for linear inverse problems have received much attention.

Iterative regularisation is a powerful algorithmic paradigm that has been
successfully employed for many inverse problems [14, Chapters 6 and 7] [30].
Classical iterative methods for inverse problems include (accelerated)
Landweber method, conjugate gradient method, Levenberg-Marquardt method, and
Gauss-Newton method, to name a few. The per-iteration computational bottleneck
of many iterative methods lies in utilising all the data at each iteration,
which can be of a prohibitively large size. For example, this occurs while
computing the derivative of an objective. One promising strategy to overcome
this challenge is stochastic gradient descent (SGD), due to Robbins and Monro
[40]. SGD decomposes the original problem into (finitely many) sub-problems,
and then at each iteration uses only a single datum, or a mini-batch of data,
typically selected uniformly at random. This greatly reduces the computational
complexity per-iteration, and enjoys excellent scalability with respect to
data size. In the standard, and best studied setting, ${\cal X}$ and ${\cal
Y}$ are finite dimensional Euclidean spaces and the corresponding data fitting
objective is the (rescaled) least squares
$\Psi({\boldsymbol{\mathsf{x}}})=\frac{1}{2N}\|{\mathbf{A}}{\boldsymbol{\mathsf{x}}}-{\boldsymbol{\mathsf{y}}}\|_{{\cal
Y}}^{2}=\frac{1}{N}\sum_{i=1}^{N}\frac{1}{2}\|{\mathbf{A}}_{i}{\boldsymbol{\mathsf{x}}}-{\boldsymbol{\mathsf{y}}}_{i}\|_{{\cal
Y}}^{2}$. In this setting SGD takes the form

$\displaystyle{\boldsymbol{\mathsf{x}}}_{k+1}={\boldsymbol{\mathsf{x}}}_{k}-\mu_{k+1}{\mathbf{A}}_{i_{k+1}}^{\ast}({\mathbf{A}}_{i_{k+1}}{\boldsymbol{\mathsf{x}}}_{k}-{\boldsymbol{\mathsf{y}}}_{i_{k+1}}),\quad
k=0,1,\ldots,$

where $\mu_{k}>$ is the step-size, ${i_{k+1}}$ a randomly selected index,
${\mathbf{A}}_{i}$ the $i$-th row of a matrix ${\mathbf{A}}$, and
${\boldsymbol{\mathsf{y}}}_{i}$ the $i$-th entry of
${\boldsymbol{\mathsf{y}}}$. In the seminal work [40], Robbins and Monro
presented SGD as a Markov chain, laying the groundwork for the field of
stochastic approximation [32]. SGD has since had a major impact on statistical
inference and machine learning, especially for the training of neural
networks. SGD has been extensively studied in the Euclidean setting; see [2]
for an overview of the convergence theory from the viewpoint of optimisation.

SGD has also been a popular method for image reconstruction, especially in
medical imaging. For example, the (randomised) Kaczmarz method is a reweighted
version of SGD that has been extensively used in computed tomography [18, 37].
Other applications of SGD and its variants include optical tomography [6],
phonon transmission coefficient recovery [15], positron emission tomography
[31], as well as general sparse recovery [42, 43]. For linear inverse problems
in Euclidean spaces, Jin and Lu [24] gave a first proof of convergence of SGD
iterates towards the minimum norm solution, and analysed the regularising
behaviour in the presence of noise; see [22, 25, 26, 34, 39] for further
convergence results, a posteriori stopping rules (discrepancy principle),
nonlinear problems, and general step-size schedules, etc.

Iterative methods in Euclidean and Hilbert spaces are effective for
reconstructing smooth solutions but fail to capture special features of the
solutions, such as sparsity and piecewise constancy. In practice, many imaging
inverse problems are more adequately described in non-Hilbert settings,
including sequence spaces $\ell^{p}(\mathbb{R})$ and Lebesgue spaces ${\cal
L}^{p}(\Omega)$, with $p\in[1,\infty]\setminus\\{2\\}$, which requires
changing either the solution, the data space, or both. For example, inverse
problems with impulse noise are better modelled by setting the data space
${\cal Y}$ to a Lebesgue space ${\cal L}^{p}(\Omega)$ with $p\approx 1$ [11],
whereas the recovery of sparse solutions is modelled by doing the same to the
solution space ${\cal X}$ [4]. Thus, it is of great importance to develop and
analyse algorithms for inverse problems in Banach spaces, and this has
received much attention [41, 46]. For the Landweber method for linear inverse
problems in Banach spaces, Schöpfer et al [44] were the first to prove strong
convergence of the iterates under a suitable step-size schedule for a smooth
and uniformly convex Banach space ${\cal X}$ and an arbitrary Banach space
${\cal Y}$. This has since been extended and refined in various aspects, e.g.
regarding acceleration [45, 49, 17, 51], nonlinear forward models [12, 35],
and Gauss-Newton methods [29].

In this work, we investigate SGD for inverse problems in Banach spaces, which
has thus far lagged behind due to outstanding challenges in extending the
analysis of standard Hilbert space approaches to the Banach space setting. The
main challenges in analysing SGD-like gradient-based methods in Banach spaces
are two-fold:

  1. 1.

The use of duality maps results in non-linear update rules, which greatly
complicates the convergence analysis. For example, the (expected) difference
between successive updates can no longer be identified as the (sub-)gradient
of the objective.

  2. 2.

Due to geometric characteristics of Banach spaces, it is more common to use
the Bregman distance for the convergence analysis, which results in the loss
of useful algebraic tools, e.g. triangle inequality and bias-variance
decomposition, that are typically needed for the analysis.

In this work, we develop an SGD approach for the numerical solution of linear
inverse problems in Banach spaces, using the sub-gradient approach based on
duality maps, and present a novel convergence analysis. We first consider the
case of exact data, and show that SGD iterates converge to a minimising
solution (first almost surely and then in expectation) under standard
assumptions on summability of step-sizes, and geometric properties of the
space ${\cal X}$, cf. Theorems 3.8 and 3.10. This solution is identified as
the minimum norm solution if the initial guess ${\boldsymbol{\mathsf{x}}}_{0}$
satisfies the range condition ${\cal J}_{p}^{{\cal
X}}({\boldsymbol{\mathsf{x}}}_{0})\in\overline{\rm{range}({\mathbf{A}}^{\ast})}$.
Further, we give a convergence rate in Theorem 3.14 when the forward operator
${\mathbf{A}}$ satisfies a conditional stability estimate. In case of noisy
observations, we show the regularising property of SGD, for properly chosen
stopping indices, cf. Theorem 4.3. The analysis rests on a descent property in
Lemma 3.6 and Robbins-Siegmund theorem for almost super-martingales. In
addition, we perform extensive numerical experiments on a model inverse
problem (linear integral equation) and computed tomography (with parallel beam
geometry) to illustrate distinct features of the proposed Banach space SGD,
and examine the influence of various factors, such as the choice of the spaces
${\cal X}$ and ${\cal Y}$, mini-batch size and noise characteristics (Gaussian
or impulse).

When finalising the paper, we became aware of the independent and simultaneous
work [27] on a stochastic mirror descent method for linear inverse problems
between a Banach space $\mathcal{X}$ and a Hilbert space $\mathcal{Y}$. The
method is a randomised version of the well-known Landweber-Kaczmarz method.
The authors prove convergence results under a priori stopping rules, and also
establish an order-optimal convergence rate when the exact solution
${\boldsymbol{\mathsf{x}}}^{\dagger}$ satisfies a benchmark source condition,
by interpreting the method as a randomised block gradient method applied to
the dual problem. Thus, the current work differs significantly from [27] in
terms of problem setting, main results and analysis techniques.

The rest of the paper is organised as follows. In Section 2, we recall
background materials on the geometry of Banach spaces, e.g. duality maps and
Bregman distance. In Section 3, we present the convergence of SGD for exact
data and in Section 4, we discuss the regularising property of SGD for noisy
observations. Finally, in Section 5, we provide some experimental results on a
model inverse problem and computed tomography. In the Appendix we collect
several useful inequalities and auxiliary estimates.

Throughout, let ${\cal X}$ and ${\cal Y}$ be two real Banach spaces, with
their norms denoted by $\|\cdot\|_{{\cal X}}$ and $\|\cdot\|_{{\cal Y}}$,
respectively. ${\cal X}^{\ast}$ and ${\cal Y}^{\ast}$ are their respective
dual spaces, with their norms denoted by $\|\cdot\|_{{\cal X}^{\ast}}$ and
$\|\cdot\|_{{\cal Y}^{\ast}}$, respectively. For
${\boldsymbol{\mathsf{x}}}\in{\cal X}$ and
${\boldsymbol{\mathsf{x}}}^{\ast}\in{\cal X}^{\ast}$, we denote the
corresponding duality pairing by
$\left<{{\boldsymbol{\mathsf{x}}}^{\ast}},{{\boldsymbol{\mathsf{x}}}}\right>=\left<{{\boldsymbol{\mathsf{x}}}^{\ast}},{{\boldsymbol{\mathsf{x}}}}\right>_{{\cal
X}^{\ast}\times{\cal
X}}={\boldsymbol{\mathsf{x}}}^{\ast}({\boldsymbol{\mathsf{x}}})$. For a
continuous linear operator ${\mathbf{A}}:{\cal X}\rightarrow{\cal Y}$, we use
$\|{\mathbf{A}}\|_{{\cal X}\to{\cal Y}}$ to denote the operator norm (often
with the subscript omitted). The adjoint of ${\mathbf{A}}$ is denoted as
${\mathbf{A}}^{\ast}:{\cal Y}^{\ast}\rightarrow{\cal X}^{\ast}$ and it is a
continuous linear operator, with $\|{\mathbf{A}}\|_{{\cal X}\to{\cal
Y}}=\|{\mathbf{A}}^{\ast}\|_{{\cal Y}^{\ast}\to{\cal X}^{\ast}}$. The
conjugate exponent of $p\in(1,\infty)$ is denoted by $p^{\ast}$, such that
$1/p+1/p^{\ast}=1$ holds. The Cauchy-Schwarz inequality of the following form
holds for any ${\boldsymbol{\mathsf{x}}}\in{\cal X}$ and
${\boldsymbol{\mathsf{x}}}^{\ast}\in{\cal X}^{\ast}$

$|\left<{{\boldsymbol{\mathsf{x}}}^{\ast}},{{\boldsymbol{\mathsf{x}}}}\right>|\leq\|{\boldsymbol{\mathsf{x}}}^{\ast}\|_{{\cal
X}^{\ast}}\|{\boldsymbol{\mathsf{x}}}\|_{{\cal X}}.$ (2)

For reals $a,b$ we write $a\wedge b=\min\\{a,b\\}$ and $a\vee
b=\max\\{a,b\\}$. By $({\cal F}_{k})_{k\in\mathbb{N}}$, we denote the natural
filtration, i.e. a growing sequence of $\sigma$-algebras such that ${\cal
F}_{k}\subset{\cal F}_{k+1}\subset{\cal F}$, for all $k\in\mathbb{N}$ and a
$\sigma$-algebra ${\cal F}$, and ${\cal F}_{k}$ is generated by random indices
$i_{j}$, for $j\leq k$. In the context of SGD, $k\in\mathbb{N}$ is the
iteration number and ${\cal F}_{k}$ denotes the iteration history, that is,
information available at time $k$, and for a given initialisation
${\boldsymbol{\mathsf{x}}}_{0}$, we can identify ${\cal
F}_{k}=\sigma({\boldsymbol{\mathsf{x}}}_{1},\ldots,{\boldsymbol{\mathsf{x}}}_{k})$.
For a filtration $({\cal F}_{k})_{k\in\mathbb{N}}$ we denote by
$\mathbb{E}_{k}[\cdot]=\mathbb{E}[\cdot\mid{\boldsymbol{\mathsf{x}}}_{1},\ldots{\boldsymbol{\mathsf{x}}}_{k}]$
the conditional expectation with respect to ${\cal F}_{k}$. A sequence of
random variables $(x_{k})_{k\in\mathbb{N}}$ (adapted to the filtration
$(\mathcal{F}_{k})_{k\in\mathbb{N}}$) is a called super-martingale if
$\mathbb{E}_{k}[x_{k+1}]\leq x_{k}.$ Throughout, the notation a.s. denotes
almost sure events.

##  2 Preliminaries on Banach spaces

In this section we recall relevant concepts from Banach space theory and the
geometry of Banach spaces.

###  2.1 Duality map

In a Hilbert space ${\cal H}$, for every ${\boldsymbol{\mathsf{x}}}\in{\cal
H}$, there exists a unique ${\boldsymbol{\mathsf{x}}}^{\ast}\in{\cal
H}^{\ast}$ such that
$\left<{{\boldsymbol{\mathsf{x}}}^{\ast}},{{\boldsymbol{\mathsf{x}}}}\right>=\|{\boldsymbol{\mathsf{x}}}\|_{{\cal
H}}\|{\boldsymbol{\mathsf{x}}}^{\ast}\|_{{\cal H}^{\ast}}$ and
$\|{\boldsymbol{\mathsf{x}}}^{\ast}\|_{{\cal
H}^{\ast}}=\|{\boldsymbol{\mathsf{x}}}\|_{{\cal H}}$, by the Riesz
representation theorem. For Banach spaces, however, such an
${\boldsymbol{\mathsf{x}}}^{\ast}$ is not necessarily unique, motivating the
notion of duality maps.

######  Definition 2.1 (Duality map).

For any $p>1$, a _duality map_ ${\cal J}_{p}^{{\cal X}}:{\cal X}\rightarrow
2^{{\cal X}^{\ast}}$ is the sub-differential of the (convex) functional
$\frac{1}{p}\|{\boldsymbol{\mathsf{x}}}\|_{{\cal X}}^{p}$

$\displaystyle{\cal J}_{p}^{{\cal
X}}({\boldsymbol{\mathsf{x}}})=\left\\{{\boldsymbol{\mathsf{x}}}^{\ast}\in{\cal
X}^{\ast}:\left<{{\boldsymbol{\mathsf{x}}}^{\ast}},{{\boldsymbol{\mathsf{x}}}}\right>=\|{\boldsymbol{\mathsf{x}}}\|_{{\cal
X}}\|{\boldsymbol{\mathsf{x}}}^{\ast}\|_{{\cal X}^{\ast}},\text{ and
}\|{\boldsymbol{\mathsf{x}}}\|_{{\cal
X}}^{p-1}=\|{\boldsymbol{\mathsf{x}}}^{\ast}\|_{{\cal X}^{\ast}}\right\\},$
(3)

with gauge function $t\mapsto t^{p-1}$. A single-valued selection of ${\cal
J}_{p}^{{\cal X}}$ is denoted by $\jmath_{p}^{{\cal X}}$.

In practice, the choice of the power parameter $p$ depends on geometric
properties of the space ${\cal X}$. For single-valued duality maps, we use
${\cal J}_{p}^{{\cal X}}$ and $\jmath_{p}^{{\cal X}}$ interchangeably. Next we
recall standard notions of smoothness and convexity of Banach spaces. For an
overview of Banach space geometry, we refer an interested reader to the
monographs [9, 10, 46].

######  Definition 2.2.

Let ${\cal X}$ be a Banach space. ${\cal X}$ is said to be reflexive if the
canonical map
${\boldsymbol{\mathsf{x}}}\mapsto\widehat{\boldsymbol{\mathsf{x}}}$ between
${\cal X}$ and the bidual ${\cal X}^{\ast\ast}$, defined by
$\widehat{\boldsymbol{\mathsf{x}}}({\boldsymbol{\mathsf{x}}}^{\ast})={\boldsymbol{\mathsf{x}}}^{\ast}({\boldsymbol{\mathsf{x}}})$,
is surjective. ${\cal X}$ is smooth if for every
$0\neq{\boldsymbol{\mathsf{x}}}\in{\cal X}$ there is a unique
${\boldsymbol{\mathsf{x}}}^{\ast}\in{\cal X}^{\ast}$ such that
$\left<{{\boldsymbol{\mathsf{x}}}^{\ast}},{{\boldsymbol{\mathsf{x}}}}\right>=\|{\boldsymbol{\mathsf{x}}}\|_{{\cal
X}}$ and $\|{\boldsymbol{\mathsf{x}}}^{\ast}\|_{{\cal X}^{\ast}}=1$. The
function $\delta_{\cal X}:(0,2]\rightarrow\mathbb{R}$ defined as

$\displaystyle\delta_{\cal
X}(\tau)=\inf\Big{\\{}1-\tfrac{1}{2}\|{\boldsymbol{\mathsf{z}}}+{\boldsymbol{\mathsf{w}}}\|_{{\cal
X}}:\|{\boldsymbol{\mathsf{z}}}\|_{{\cal
X}}=\|{\boldsymbol{\mathsf{w}}}\|_{{\cal
X}}=1,\|{\boldsymbol{\mathsf{z}}}-{\boldsymbol{\mathsf{w}}}\|_{{\cal
X}}\geq\tau\Big{\\}}$

is the _modulus of convexity_ of ${\cal X}$. ${\cal X}$ is said to be
_uniformly convex_ if $\delta_{\cal X}(\tau)>0$ for all $\tau\in(0,2]$, and
_$p$ -convex_, for $p>1$, if $\delta_{\cal X}(\tau)\geq K_{p}\tau^{p}$ for
some $K_{p}>0$ and all $\tau\in(0,2]$. The function $\rho_{\cal
X}:{[0,\infty)\to[0,\infty)}$ defined as

$\displaystyle\rho_{\cal
X}(\tau)=\sup\Big{\\{}\frac{\|{\boldsymbol{\mathsf{z}}}+\tau{\boldsymbol{\mathsf{w}}}\|_{{\cal
X}}+\|{\boldsymbol{\mathsf{z}}}-\tau{\boldsymbol{\mathsf{w}}}\|_{{\cal
X}}}{2}-1:\|{\boldsymbol{\mathsf{z}}}\|_{{\cal
X}}=\|{\boldsymbol{\mathsf{w}}}\|_{{\cal X}}=1\Big{\\}}$

is the _modulus of smoothness_ of ${\cal X}$, and is a convex and continuous
function such that $\frac{\rho_{\cal X}(\tau)}{\tau}$ is a non-decreasing
function with $\rho_{\cal X}(\tau)\leq\tau$. ${\cal X}$ is said to be
_uniformly smooth_ if $\lim_{\tau\searrow 0}\frac{\rho_{\cal
X}(\tau)}{\tau}=0$, and _$p$ -smooth_, for $p>1$, if $\rho_{\cal X}(\tau)\leq
K_{p}\tau^{p}$ for some $K_{p}>0$ and all $\tau\in(0,\infty)$.

The following relationships between Banach spaces and duality maps will be
used extensively.

######  Theorem 2.3 ([46, Theorems 2.52 and 2.53, and Lemma 5.16]).

  1. (i)

For every ${\boldsymbol{\mathsf{x}}}\in{\cal X}$, the set ${\cal J}_{p}^{{\cal
X}}({\boldsymbol{\mathsf{x}}})$ is non-empty, convex, and weakly-$\star$
closed in ${\cal X}^{\ast}$.

  2. (ii)

${\cal X}$ is $p$-smooth if and only if ${\cal X}^{\ast}$ is
$p^{\ast}$-convex. ${\cal X}$ is $p$-convex if and only if ${\cal X}^{\ast}$
is $p^{\ast}$-smooth.

  3. (iii)

${\cal X}$ is smooth if and only if ${\cal J}_{p}^{{\cal X}}$ is single
valued. If ${\cal X}$ is convex of power type and smooth, then ${\cal
J}_{p}^{{\cal X}}$ is invertible and $\big{(}{\cal J}_{p}^{{\cal
X}}\big{)}^{-1}={\cal J}_{p^{\ast}}^{{\cal X}^{\ast}}$. If ${\cal X}$ is
uniformly smooth and uniformly convex, then ${\cal J}_{p}^{{\cal X}}$ and
${\cal J}_{p^{\ast}}^{{\cal X}^{\ast}}$ are both uniformly continuous.

  4. (iv)

Let ${\cal X}$ be a uniformly smooth Banach space with duality map ${\cal
J}_{p}^{{\cal X}}$ with $p\geq 2$. Then, for all
${\boldsymbol{\mathsf{x}}},\widetilde{\boldsymbol{\mathsf{x}}}\in{\cal X}$,
there holds

$\displaystyle\|{\cal J}_{p}^{{\cal X}}({\boldsymbol{\mathsf{x}}})-{\cal
J}_{p}^{{\cal X}}(\widetilde{\boldsymbol{\mathsf{x}}})\|_{{\cal
X}^{\ast}}^{p^{\ast}}\leq C\max\\{1,\|{\boldsymbol{\mathsf{x}}}\|_{{\cal
X}},\|\widetilde{\boldsymbol{\mathsf{x}}}\|_{{\cal
X}}\\}^{p}\,{\overline{\rho}_{\cal
X}}(\|{\boldsymbol{\mathsf{x}}}-\widetilde{\boldsymbol{\mathsf{x}}}\|_{{\cal
X}})^{p^{\ast}},$

where $\overline{\rho}_{\cal X}(\tau)=\rho_{\cal X}(\tau)/\tau$ is a modulus
of smoothness function such that $\overline{\rho}(\tau)\leq 1$.

Next we list some common Banach spaces, the corresponding duality maps and
convexity and smoothness properties.

######  Example 2.4.

  1. (i)

A Hilbert space ${\cal X}$ is $2$-smooth and $2$-convex, and ${\cal
J}_{2}^{{\cal X}}$ is the identity.

  2. (ii)

If ${\cal X}$ is smooth, then ${\cal J}_{p}^{{\cal X}}$ is the Gateaux
derivative of the functional
${\boldsymbol{\mathsf{x}}}\mapsto\frac{1}{p}\|{\boldsymbol{\mathsf{x}}}\|_{{\cal
X}}^{p}$.

  3. (iii)

If ${\cal X}=\ell^{r}(\mathbb{R})$ with $1<r<\infty$, then ${\cal
J}_{p}^{{\cal X}}$ is single-valued, and the duality map is given by ${\cal
J}_{p}^{{\cal
X}}({\boldsymbol{\mathsf{x}}})=\|{\boldsymbol{\mathsf{x}}}\|_{r}^{p-r}|{\boldsymbol{\mathsf{x}}}|^{r-1}\operatorname{sign}({\boldsymbol{\mathsf{x}}}).$
Moreover, ${\cal J}_{p}^{{\cal X}}=\nabla({\frac{1}{p}\|\cdot\|_{{\cal
X}}^{p}})$ since ${\cal X}$ is smooth.

  4. (iv)

Lebesgue spaces ${\cal L}^{p}(\Omega)$, Sobolev spaces $W^{s,p}(\Omega)$, with
$s>0$, (for an open bounded domain $\Omega$), and sequence spaces
$\ell^{p}(\mathbb{R})$ are $p\wedge 2$-smooth and $p\vee 2$-convex, for
$1<p<\infty$. For $p\in\\{1,\infty\\}$, they are neither smooth nor strictly
convex.

###  2.2 Bregman distance

Due to the geometry of Banach spaces, it is often more convenient to use the
Bregman distance than the standard Banach space norm $\|\cdot\|_{{\cal X}}$ in
the convergence analysis.

######  Definition 2.5 (Bregman distance).

For a smooth Banach space ${\cal X}$, the functional

$\displaystyle{\mathbf{B}}_{p}({\boldsymbol{\mathsf{z}}},{\boldsymbol{\mathsf{w}}})$
$\displaystyle=\frac{1}{p^{\ast}}\|{\boldsymbol{\mathsf{z}}}\|_{{\cal
X}}^{p}+\frac{1}{p}\|{\boldsymbol{\mathsf{w}}}\|_{{\cal X}}^{p}-\left<{{\cal
J}_{p}^{{\cal
X}}({\boldsymbol{\mathsf{z}}})},{{\boldsymbol{\mathsf{w}}}}\right>,$

is called the _Bregman distance_ , where $1/p+1/p^{\ast}=1$.

Note that the dependence of the Bregman distance
${\mathbf{B}}_{p}({\boldsymbol{\mathsf{z}}},{\boldsymbol{\mathsf{w}}})$ on the
space ${\cal X}$ is omitted, which is often clear from the context. The
Bregman distance does not satisfy the triangle inequality, and is generally
non-symmetric. Thus it is not a distance. The next theorem lists useful
properties of the Bregman distance, which show the relationship between the
geometry of the underlying Banach space and duality maps.

######  Theorem 2.6 ([46, Theorem 2.60, Lemmas 2.62 and 2.63]).

The following properties hold.

  1. (i)

If ${\cal X}$ is smooth and reflexive, then
${\mathbf{B}}_{p}({\boldsymbol{\mathsf{z}}},{\boldsymbol{\mathsf{w}}})={{\mathbf{B}}_{p^{\ast}}\\!\Big{(}\\!{\cal
J}_{p}^{{\cal X}}({\boldsymbol{\mathsf{w}}}),{\cal J}_{p}^{{\cal
X}}({\boldsymbol{\mathsf{z}}})\\!\Big{)}}.$

  2. (ii)

Bregman distance satisfies the three-point identity

$\displaystyle{\mathbf{B}}_{p}({\boldsymbol{\mathsf{z}}},{\boldsymbol{\mathsf{w}}})={\mathbf{B}}_{p}({\boldsymbol{\mathsf{z}}},{\boldsymbol{\mathsf{v}}})+{\mathbf{B}}_{p}({\boldsymbol{\mathsf{v}}},{\boldsymbol{\mathsf{w}}})+\left<{{\cal
J}_{p}^{{\cal X}}({\boldsymbol{\mathsf{v}}})-{\cal J}_{p}^{{\cal
X}}({\boldsymbol{\mathsf{z}}})},{{\boldsymbol{\mathsf{w}}}-{\boldsymbol{\mathsf{v}}}}\right>.$
(4)

  3. (iii)

If ${\cal X}$ is $p$-convex, then it is reflexive, $p\geq 2$ and there exists
$C_{p}>0$ such that

$\displaystyle{\mathbf{B}}_{p}({\boldsymbol{\mathsf{z}}},{\boldsymbol{\mathsf{w}}})\geq
p^{-1}C_{p}\|{\boldsymbol{\mathsf{w}}}-{\boldsymbol{\mathsf{z}}}\|_{{\cal
X}}^{p}.$ (5)

  4. (iv)

If ${\cal X}^{\ast}$ is $p^{\ast}$-smooth, then it is reflexive, $p^{\ast}\leq
2$ and there exists $G_{p^{\ast}}>0$ such that

$\displaystyle{{\mathbf{B}}_{p^{\ast}}\big{(}{\boldsymbol{\mathsf{z}}}^{\ast},{\boldsymbol{\mathsf{w}}}^{\ast}\big{)}\leq(p^{\ast})^{-1}G_{p^{\ast}}\|{\boldsymbol{\mathsf{w}}}^{\ast}-{\boldsymbol{\mathsf{z}}}^{\ast}\|_{{\cal
X}^{\ast}}^{p^{\ast}}}.$ (6)

  5. (v)

${\mathbf{B}}_{p}({\boldsymbol{\mathsf{z}}},{\boldsymbol{\mathsf{w}}})\geq 0$,
and if ${\cal X}$ is uniformly convex,
${\mathbf{B}}_{p}({\boldsymbol{\mathsf{z}}},{\boldsymbol{\mathsf{w}}})=0$ if
and only if ${\boldsymbol{\mathsf{z}}}={\boldsymbol{\mathsf{w}}}$.

  6. (vi)

${\mathbf{B}}_{p}({\boldsymbol{\mathsf{z}}},{\boldsymbol{\mathsf{w}}})$ is
continuous in the second argument. If ${\cal X}$ is smooth and uniformly
convex, then ${\cal J}_{p}^{{\cal X}}$ is continuous on bounded subsets and
${\mathbf{B}}_{p}({\boldsymbol{\mathsf{z}}},{\boldsymbol{\mathsf{w}}})$ is
continuous in its first argument.

##  3 Convergence analysis for exact data

Now we develop an SGD type approach for problem (1) and analyse its
convergence. Throughout, we make the following assumption on the Banach spaces
${\cal X}$ and ${\cal Y}$, unless indicated otherwise.

######  Assumption 3.1.

The Banach space ${\cal X}$ is $p$-convex and smooth, and ${\cal Y}$ is
arbitrary.

To recover the solution ${\boldsymbol{\mathsf{x}}}^{\dagger}$, we minimise a
least-squares type objective
$\operatorname*{argmin}_{{\boldsymbol{\mathsf{x}}}\in{\cal
X}}\tfrac{1}{p}\|{\mathbf{A}}{\boldsymbol{\mathsf{x}}}-{\boldsymbol{\mathsf{y}}}\|_{{\cal
Y}}^{p},$ for some $p>1$. By ${\cal X}_{\min}$, we denote the (non-empty) set
of minimisers over ${\cal X}$. Among the elements of ${\cal X}_{\min}$, the
regularisation theory focuses on the so-called minimum norm solution.

######  Definition 3.2.

An element ${\boldsymbol{\mathsf{x}}}^{\dagger}\in{\cal X}$ is called a
_minimum norm solution_ (MNS ) of (1) if

${\mathbf{A}}{\boldsymbol{\mathsf{x}}}^{\dagger}={\boldsymbol{\mathsf{y}}}\quad\text{
and }\quad\|{\boldsymbol{\mathsf{x}}}^{\dagger}\|_{{\cal
X}}=\inf\\{\|{\boldsymbol{\mathsf{x}}}\|_{{\cal
X}}:{\boldsymbol{\mathsf{x}}}\in{\cal
X},\,{\mathbf{A}}{\boldsymbol{\mathsf{x}}}={\boldsymbol{\mathsf{y}}}\\}.$

The MNS ${\boldsymbol{\mathsf{x}}}^{\dagger}$ is not unique for general Banach
spaces. The following lemma states sufficient geometric assumptions on ${\cal
X}$ for uniqueness.

######  Lemma 3.3 ([46, Lemma 3.3]).

Let Assumption 3.1 hold. Then there exists a unique MNS
${\boldsymbol{\mathsf{x}}}^{\dagger}$. Furthermore, ${\cal J}_{p}^{{\cal
X}}({\boldsymbol{\mathsf{x}}}^{\dagger})\in\overline{\operatorname{range}({\mathbf{A}}^{\ast})}$,
for $1<p<\infty$. If some $\widehat{{\boldsymbol{\mathsf{x}}}}\in{\cal X}$
satisfies ${\cal J}_{p}^{{\cal
X}}(\widehat{{\boldsymbol{\mathsf{x}}}})\in\overline{\operatorname{range}({\mathbf{A}}^{\ast})}$
and
$\widehat{{\boldsymbol{\mathsf{x}}}}-{\boldsymbol{\mathsf{x}}}^{\dagger}\in\operatorname{null}({\mathbf{A}})$,
then
$\widehat{{\boldsymbol{\mathsf{x}}}}={\boldsymbol{\mathsf{x}}}^{\dagger}$.

By Lemma 3.3, the MNS ${\boldsymbol{\mathsf{x}}}^{\dagger}$ is unique modulo
the null space of ${\mathbf{A}}$, under certain smoothness and convexity
assumptions on ${\cal X}$. These conditions exclude Lebesgue and sequence
spaces ${\cal L}^{1}(\Omega)$ and $\ell^{1}(\mathbb{R})$, cf. Example 2.4(iv).
The standard Landweber method [33, 44] constructs an approximation to the MNS
${\boldsymbol{\mathsf{x}}}^{\dagger}$ by running the iterations

$\displaystyle{\boldsymbol{\mathsf{x}}}_{k+1}={\cal J}_{p^{\ast}}^{{\cal
X}^{\ast}}\left({\cal J}_{p}^{{\cal
X}}({\boldsymbol{\mathsf{x}}}_{k})-\mu_{k+1}{\mathbf{A}}^{\ast}\jmath_{p}^{{\cal
Y}}({\mathbf{A}}{\boldsymbol{\mathsf{x}}}_{k}-{\boldsymbol{\mathsf{y}}})\right),\quad
k=0,1,\ldots,$ (7)

where $\mu_{k+1}>0$ is the step-size. Asplund’s theorem [46, Theorem 2.28]
allows characterising the duality map as the sub-differential, ${\cal
J}_{p}^{{\cal X}}=\partial({\frac{1}{p}\|\cdot\|_{{\cal X}}^{p}})$ for $p>1$.
This identifies the descent direction ${\mathbf{A}}^{\ast}\jmath_{p}^{{\cal
Y}}({\mathbf{A}}{\boldsymbol{\mathsf{x}}}_{k}-{\boldsymbol{\mathsf{y}}})$ as
the sub-gradient: ${\mathbf{A}}^{\ast}\jmath_{p}^{{\cal
Y}}({\mathbf{A}}{\boldsymbol{\mathsf{x}}}_{k}-{\boldsymbol{\mathsf{y}}})=\partial(\frac{1}{p}\|{\mathbf{A}}\cdot-{\boldsymbol{\mathsf{y}}}\|_{{\cal
Y}})({\boldsymbol{\mathsf{x}}}_{k})$. Note that ${\cal J}_{p}^{{\cal X}}$ is
single valued by Assumption 3.1 and Theorem 2.3, though ${\cal J}_{p}^{{\cal
Y}}$ is not. For well selected step-sizes, Landweber iterations (7) converge
to an MNS of (1) [44, Theorem 3.3].

The evaluation of the sub-gradient ${\mathbf{A}}^{\ast}\jmath_{p}^{{\cal
Y}}({\mathbf{A}}{\boldsymbol{\mathsf{x}}}_{k}-{\boldsymbol{\mathsf{y}}})$
represents the main per-iteration cost of the iteration (7). In this work, we
consider the following Kaczmarz type setting:

$\displaystyle{\mathbf{A}}=\begin{pmatrix}{\mathbf{A}}_{1}\\\ \vdots\\\
{\mathbf{A}}_{N}\end{pmatrix}\quad\text{and}\quad{\mathbf{A}}{\boldsymbol{\mathsf{x}}}=\begin{pmatrix}{\mathbf{A}}_{1}{\boldsymbol{\mathsf{x}}}\\\
\vdots\\\
{\mathbf{A}}_{N}{\boldsymbol{\mathsf{x}}}\end{pmatrix}=\begin{pmatrix}{\boldsymbol{\mathsf{y}}}_{1}\\\
\vdots\\\ {\boldsymbol{\mathsf{y}}}_{N}\end{pmatrix},$ (8)

where ${\mathbf{A}}_{i}:{\cal X}\rightarrow{\cal Y}_{i}$,
${\boldsymbol{\mathsf{y}}}_{i}\in{\cal Y}_{i}$, for
$i\in[N]=\\{1,\ldots,N\\}$. Problem (8) is defined on the direct product
$(\otimes_{i=1}^{N}{\cal Y}_{i},\ell^{r})$, equipped with the $\ell^{r}$ norm,
for $r\geq 1$

$\displaystyle\|{\boldsymbol{\mathsf{y}}}\|_{\mathcal{Y}}:=\|({\boldsymbol{\mathsf{y}}}_{1},\ldots,{\boldsymbol{\mathsf{y}}}_{N})\|_{\mathcal{Y}}=\|(\|{\boldsymbol{\mathsf{y}}}_{1}\|_{{\cal
Y}_{1}},\ldots,\|{\boldsymbol{\mathsf{y}}}_{N}\|_{{\cal
Y}_{N}})\|_{\ell^{r}}=\Big{(}\sum_{i=1}^{N}\|{\boldsymbol{\mathsf{y}}}_{i}\|_{\mathcal{Y}_{i}}^{r}\Big{)}^{1/r}.$
(9)

Below we identify ${\cal Y}_{i}={\cal Y}$ for notational brevity, and use
$\|\cdot\|_{{\cal Y}}$ to denote both the norm of the direct product space and
the component spaces, though all the relevant proofs and concepts easily
extend to the general case. Then the objective
$\Psi({\boldsymbol{\mathsf{x}}})$ is given by

$\Psi({\boldsymbol{\mathsf{x}}})=\frac{1}{N}\sum_{i=1}^{N}\Psi_{i}({\boldsymbol{\mathsf{x}}}),\quad\mbox{with
}\Psi_{i}({\boldsymbol{\mathsf{x}}})=\frac{1}{p}\|{\mathbf{A}}_{i}{\boldsymbol{\mathsf{x}}}-{\boldsymbol{\mathsf{y}}}_{i}\|_{{\cal
Y}}^{p}.$

Note that for many common imaging problems we use
$\mathcal{Y}=\ell^{p}(\mathbb{R})$, which then naturally gives
$\Psi({\boldsymbol{\mathsf{x}}})=\frac{1}{pN}\|{\mathbf{A}}{\boldsymbol{\mathsf{x}}}-{\boldsymbol{\mathsf{y}}}\|_{{\cal
Y}}^{p}$. To reduce the computational cost per-iteration, we exploit the
finite-sum structure of the objective $\Psi({\boldsymbol{\mathsf{x}}})$ and
adopt SGD iterations of the form

$\displaystyle{\boldsymbol{\mathsf{x}}}_{k+1}={\cal J}_{p^{\ast}}^{{\cal
X}^{\ast}}\left({\cal J}_{p}^{{\cal
X}}({\boldsymbol{\mathsf{x}}}_{k})-\mu_{k+1}{\boldsymbol{\mathsf{g}}}_{k+1}\right),$
(10)

where
${\boldsymbol{\mathsf{g}}}_{k+1}=g({\boldsymbol{\mathsf{x}}}_{k},{\boldsymbol{\mathsf{y}}},i_{k+1})$
is the stochastic update direction given by

$\displaystyle
g({\boldsymbol{\mathsf{x}}},{\boldsymbol{\mathsf{y}}},i)={\mathbf{A}}_{i}^{\ast}\jmath_{p}^{{\cal
Y}}({\mathbf{A}}_{i}{\boldsymbol{\mathsf{x}}}-{\boldsymbol{\mathsf{y}}}_{i})=\partial\Big{(}\tfrac{1}{p}\|{\mathbf{A}}_{i}\cdot-{\boldsymbol{\mathsf{y}}}_{i}\|_{{\cal
Y}}^{p}\Big{)}({\boldsymbol{\mathsf{x}}}),$ (11)

and the random index $i_{k}$ is sampled uniformly over the index set $[N]$,
independent of ${\boldsymbol{\mathsf{x}}}_{k}$. Clearly, it is an unbiased
estimator of the sub-gradient $\partial\Psi({\boldsymbol{\mathsf{x}}})$, i.e.
$\mathbb{E}[g({\boldsymbol{\mathsf{x}}},{\boldsymbol{\mathsf{y}}},i)]=\partial\Psi({\boldsymbol{\mathsf{x}}})$,
and the per-iteration cost is reduced by a factor of $N$.

######  Remark 3.4.

In the model (8), if ${\cal Y}$ admits a complemented sum ${\cal
Y}=\sum_{i=1}^{N}{\cal Y}_{i}$, we can take the (internal ) direct sum
$(\oplus_{i=1}^{N}{\cal Y}_{i},\ell^{r})$, so that
${\boldsymbol{\mathsf{y}}}={\boldsymbol{\mathsf{y}}}_{1}+\ldots+{\boldsymbol{\mathsf{y}}}_{N}$
and the corresponding norm $\|{\boldsymbol{\mathsf{y}}}\|=\|(\|{\rm
Proj}_{{\cal Y}_{1}}({\boldsymbol{\mathsf{y}}})\|_{{\cal Y}_{1}},\ldots,\|{\rm
Proj}_{{\cal Y}_{N}}({\boldsymbol{\mathsf{y}}})\|_{{\cal Y}_{N}})\|_{r}$. With
this identification the spaces $(\otimes_{i=1}^{N}{\cal Y}_{i},\ell^{r})$ and
$(\oplus_{i=1}^{N}{\cal Y}_{i},\ell^{r})$ are isometrically isomorphic [48]
and the norms are equivalent for all $r\geq 1$.

We now collect some useful properties about the objective $\Psi$ and the
Bregman divergence. Throughout, $L_{\max}=\max_{i\in[N]}\|{\mathbf{A}}_{i}\|.$
Note that $c_{N}=1/N$ if ${\cal Y}={\cal L}^{p}(\Omega)$.

######  Lemma 3.5.

For all $i\in[N]$, ${\boldsymbol{\mathsf{x}}}\in{\cal X}$, and any
$\widehat{\boldsymbol{\mathsf{x}}}\in\mathcal{X}_{\min}$ (such that
${\mathbf{A}}\widehat{\boldsymbol{\mathsf{x}}}={\boldsymbol{\mathsf{y}}}$), we
have

$\displaystyle\left<{\partial\Psi_{i}({\boldsymbol{\mathsf{x}}})},{{\boldsymbol{\mathsf{x}}}-\widehat{\boldsymbol{\mathsf{x}}}}\right>=p\Psi_{i}({\boldsymbol{\mathsf{x}}})\quad\text{and}\quad\left<{\partial\Psi({\boldsymbol{\mathsf{x}}})},{{\boldsymbol{\mathsf{x}}}-\widehat{\boldsymbol{\mathsf{x}}}}\right>=p\Psi({\boldsymbol{\mathsf{x}}}).$
(12)

Moreover,
$\Psi_{i}({\boldsymbol{\mathsf{x}}})\leq\frac{\|{\mathbf{A}}_{i}\|^{p}}{C_{p}}{\mathbf{B}}_{p}({\boldsymbol{\mathsf{x}}},\widehat{\boldsymbol{\mathsf{x}}})$,
$\Psi({\boldsymbol{\mathsf{x}}})\leq\frac{L_{\max}^{p}}{C_{p}}{\mathbf{B}}_{p}({\boldsymbol{\mathsf{x}}},\widehat{\boldsymbol{\mathsf{x}}})$,
and for some $C_{N}>0$ we have
$\Psi({\boldsymbol{\mathsf{x}}})\geq\frac{C_{N}}{p}\|{\mathbf{A}}{\boldsymbol{\mathsf{x}}}-{\boldsymbol{\mathsf{y}}}\|_{{\cal
Y}}^{p}.$

###### Proof.

It follows from the identity
${\mathbf{A}}\widehat{\boldsymbol{\mathsf{x}}}={\boldsymbol{\mathsf{y}}}$ that

$\displaystyle\left<{\partial\Psi_{i}({\boldsymbol{\mathsf{x}}})},{{\boldsymbol{\mathsf{x}}}-\widehat{\boldsymbol{\mathsf{x}}}}\right>=\left<{{\mathbf{A}}_{i}^{\ast}\jmath_{p}^{{\cal
Y}}({\mathbf{A}}_{i}{\boldsymbol{\mathsf{x}}}-{\boldsymbol{\mathsf{y}}}_{i})},{{\boldsymbol{\mathsf{x}}}-\widehat{\boldsymbol{\mathsf{x}}}}\right>=\left<{\jmath_{p}^{{\cal
Y}}({\mathbf{A}}_{i}{\boldsymbol{\mathsf{x}}}-{\boldsymbol{\mathsf{y}}}_{i})},{{\mathbf{A}}_{i}{\boldsymbol{\mathsf{x}}}-{\boldsymbol{\mathsf{y}}}_{i}}\right>=p\Psi_{i}({\boldsymbol{\mathsf{x}}}).$

Since
$\partial\Psi({\boldsymbol{\mathsf{x}}})=\frac{1}{N}\sum_{i=1}^{N}\partial\Psi_{i}({\boldsymbol{\mathsf{x}}})$,
the second identity in (12) follows from the linearity of the dual product. By
the $p$-convexity of the space ${\cal X}$ and Theorem 2.6(iii), we get

$\displaystyle\Psi_{i}({\boldsymbol{\mathsf{x}}})=\frac{1}{p}\|{\mathbf{A}}_{i}{\boldsymbol{\mathsf{x}}}-{\boldsymbol{\mathsf{y}}}_{i}\|_{{\cal
Y}}^{p}=\frac{1}{p}\|{\mathbf{A}}_{i}({\boldsymbol{\mathsf{x}}}-\widehat{\boldsymbol{\mathsf{x}}})\|_{{\cal
Y}}^{p}\leq\frac{\|{\mathbf{A}}_{i}\|^{p}}{p}\|{\boldsymbol{\mathsf{x}}}-\widehat{\boldsymbol{\mathsf{x}}}\|_{{\cal
X}}^{p}\leq\frac{\|{\mathbf{A}}_{i}\|^{p}}{C_{p}}{\mathbf{B}}_{p}({\boldsymbol{\mathsf{x}}},\widehat{\boldsymbol{\mathsf{x}}}).$

The second claim follows since
$\Psi({\boldsymbol{\mathsf{x}}})=\frac{1}{N}\sum_{i=1}^{N}\Psi_{i}({\boldsymbol{\mathsf{x}}})$.
Lastly, by the norm equivalence (9) for $1<r<\infty$, there exists $C_{N}>0$
such that

$\displaystyle\Psi({\boldsymbol{\mathsf{x}}})=\frac{1}{N}\sum_{i=1}^{N}\frac{1}{p}\|{\mathbf{A}}_{i}{\boldsymbol{\mathsf{x}}}-{\boldsymbol{\mathsf{y}}}_{i}\|_{{\cal
Y}}^{p}\geq\frac{C_{N}}{p}\|{\mathbf{A}}{\boldsymbol{\mathsf{x}}}-{\boldsymbol{\mathsf{y}}}\|_{{\cal
Y}}^{p}.$

∎

We now focus on the convergence study of the iterations (10), without and with
noise in the data, and discuss convergence rates under conditional stability.

###  3.1 Convergence for the Kaczmarz model

Below the notation $\mathbb{E}[\cdot]$ denotes taking expectation with respect
to the sampling of the random indices $i_{k}$ and $\mathbb{E}_{k}[\cdot]$
denotes taking conditional expectation with respect to $\mathcal{F}_{k}$. The
remaining variables, e.g. ${\boldsymbol{\mathsf{x}}}$ and
${\boldsymbol{\mathsf{y}}}$, are measurable with respect to the underlying
probability measure. To study the convergence of SGD (10), we first establish
a descent property in terms of the Bregman distance.

######  Lemma 3.6.

Let Assumption 3.1 hold. For any $\widehat{\boldsymbol{\mathsf{x}}}\in{\cal
X}$, the iterates in (10) satisfy

$\displaystyle{\mathbf{B}}_{p}({\boldsymbol{\mathsf{x}}}_{k+1},\widehat{\boldsymbol{\mathsf{x}}})\leq{\mathbf{B}}_{p}({\boldsymbol{\mathsf{x}}}_{k},\widehat{\boldsymbol{\mathsf{x}}})-\mu_{k+1}\left<{{\boldsymbol{\mathsf{g}}}_{k+1}},{{\boldsymbol{\mathsf{x}}}_{k}-\widehat{\boldsymbol{\mathsf{x}}}}\right>+\frac{G_{p^{\ast}}}{p^{\ast}}\mu_{k+1}^{p^{\ast}}\|{\boldsymbol{\mathsf{g}}}_{k+1}\|_{{\cal
X}^{\ast}}^{p^{\ast}}.$ (13)

###### Proof.

Let
$\Delta_{k}:={\mathbf{B}}_{p}({\boldsymbol{\mathsf{x}}}_{k},\widehat{\boldsymbol{\mathsf{x}}})$.
By Definition 2.5 and expression (10), we have

$\displaystyle\Delta_{k+1}$
$\displaystyle=\frac{1}{p}\|\widehat{\boldsymbol{\mathsf{x}}}\|_{{\cal
X}}^{p}+\frac{1}{p^{\ast}}\|{\boldsymbol{\mathsf{x}}}_{k+1}\|_{{\cal
X}}^{p}-\left<{{\cal J}_{p}^{{\cal
X}}({\boldsymbol{\mathsf{x}}}_{k+1})},{\widehat{\boldsymbol{\mathsf{x}}}}\right>$
$\displaystyle{=\frac{1}{p}\|\widehat{\boldsymbol{\mathsf{x}}}\|_{{\cal
X}}^{p}+\frac{1}{p^{\ast}}\|{\cal J}_{p^{\ast}}^{{\cal X}^{\ast}}\left({\cal
J}_{p}^{{\cal
X}}({\boldsymbol{\mathsf{x}}}_{k})-\mu_{k+1}{\boldsymbol{\mathsf{g}}}_{k+1}\right)\|_{{\cal
X}}^{p}-\left<{{\cal J}_{p}^{{\cal
X}}({\boldsymbol{\mathsf{x}}}_{k+1})},{\widehat{\boldsymbol{\mathsf{x}}}}\right>}.$

Using Definition 2.1, the identity $p(p^{\ast}-1)=p^{\ast}$ and Theorem
2.3(iii), we deduce

$\displaystyle\Delta_{k+1}$
$\displaystyle=\frac{1}{p}\|\widehat{\boldsymbol{\mathsf{x}}}\|_{{\cal
X}}^{p}+\frac{1}{p^{\ast}}\|{\cal J}_{p}^{{\cal
X}}({\boldsymbol{\mathsf{x}}}_{k})-\mu_{k+1}{\boldsymbol{\mathsf{g}}}_{k+1}\|_{{\cal
X}^{\ast}}^{p(p^{\ast}-1)}-\left<{{\cal J}_{p}^{{\cal
X}}({\boldsymbol{\mathsf{x}}}_{k})-\mu_{k+1}{\boldsymbol{\mathsf{g}}}_{k+1}},{\widehat{\boldsymbol{\mathsf{x}}}}\right>$
$\displaystyle=\frac{1}{p}\|\widehat{\boldsymbol{\mathsf{x}}}\|_{{\cal
X}}^{p}+\frac{1}{p^{\ast}}\|{\cal J}_{p}^{{\cal
X}}({\boldsymbol{\mathsf{x}}}_{k})-\mu_{k+1}{\boldsymbol{\mathsf{g}}}_{k+1}\|_{{\cal
X}^{\ast}}^{p^{\ast}}-\left<{{\cal J}_{p}^{{\cal
X}}({\boldsymbol{\mathsf{x}}}_{k})-\mu_{k+1}{\boldsymbol{\mathsf{g}}}_{k+1}},{\widehat{\boldsymbol{\mathsf{x}}}}\right>.$

Since ${\cal X}$ is $p$-convex, ${\cal X}^{\ast}$ is $p^{\ast}$-smooth, cf.
Theorem 2.3(i). By [9, Corollary 5.8], this implies

$\frac{1}{p^{\ast}}\|{\boldsymbol{\mathsf{x}}}^{\ast}-\tilde{\boldsymbol{\mathsf{x}}}^{\ast}\|_{{\cal
X}^{\ast}}^{p^{\ast}}\leq\frac{1}{p^{\ast}}\|{\boldsymbol{\mathsf{x}}}^{\ast}\|_{{\cal
X}^{\ast}}^{p^{\ast}}+\frac{G_{p^{\ast}}}{p^{\ast}}\|\tilde{\boldsymbol{\mathsf{x}}}^{\ast}\|_{{\cal
X}^{\ast}}^{p^{\ast}}-\left<{{\cal J}_{p^{\ast}}^{{\cal
X}^{\ast}}({\boldsymbol{\mathsf{x}}}^{\ast})},{\tilde{\boldsymbol{\mathsf{x}}}^{\ast}}\right>,\quad\forall{\boldsymbol{\mathsf{x}}}^{\ast},\tilde{\boldsymbol{\mathsf{x}}}^{\ast}\in{\cal
X}^{\ast}.$

Using the identities $p^{\ast}(p-1)=p$ and $({\cal J}_{p}^{{\cal
X}})^{-1}={\cal J}_{{p}^{\ast}}^{{\cal X}^{\ast}}$, cf. Theorem 2.3(iii), we
get

$\displaystyle\frac{1}{p^{\ast}}\|{\cal J}_{p}^{{\cal
X}}({\boldsymbol{\mathsf{x}}}_{k})-\mu_{k+1}{\boldsymbol{\mathsf{g}}}_{k+1}\|_{{\cal
X}^{\ast}}^{p^{\ast}}$ $\displaystyle\leq\frac{1}{p^{\ast}}\|{\cal
J}_{p}^{{\cal X}}({\boldsymbol{\mathsf{x}}}_{k})\|_{{\cal
X}^{\ast}}^{p^{\ast}}+\frac{G_{p^{\ast}}}{p^{\ast}}\|\mu_{k+1}{\boldsymbol{\mathsf{g}}}_{k+1}\|_{{\cal
X}^{\ast}}^{p^{\ast}}-\left<{\mu_{k+1}{\boldsymbol{\mathsf{g}}}_{k+1}},{{\boldsymbol{\mathsf{x}}}_{k}}\right>$
$\displaystyle=\frac{1}{p^{\ast}}\|{\boldsymbol{\mathsf{x}}}_{k}\|_{{\cal
X}}^{p}+\frac{G_{p^{\ast}}}{p^{\ast}}\mu_{k+1}^{p^{\ast}}\|{\boldsymbol{\mathsf{g}}}_{k+1}\|_{{\cal
X}^{\ast}}^{p^{\ast}}-\mu_{k+1}\left<{{\boldsymbol{\mathsf{g}}}_{k+1}},{{\boldsymbol{\mathsf{x}}}_{k}}\right>.$

Combining the preceding estimates gives the desired assertion through

$\displaystyle\Delta_{k+1}$
$\displaystyle\leq\frac{1}{p}\|\widehat{\boldsymbol{\mathsf{x}}}\|_{{\cal
X}}^{p}+\frac{1}{p^{\ast}}\|{\boldsymbol{\mathsf{x}}}_{k}\|_{{\cal
X}^{\ast}}^{p}-\left<{{\cal J}_{p}^{{\cal
X}}({\boldsymbol{\mathsf{x}}}_{k})},{\widehat{\boldsymbol{\mathsf{x}}}}\right>+\frac{G_{p^{\ast}}}{p^{\ast}}\mu_{k+1}^{p^{\ast}}\|{\boldsymbol{\mathsf{g}}}_{k+1}\|_{{\cal
X}^{\ast}}^{p^{\ast}}-\mu_{k+1}\left<{{\boldsymbol{\mathsf{g}}}_{k+1}},{{\boldsymbol{\mathsf{x}}}_{k}-\widehat{\boldsymbol{\mathsf{x}}}}\right>$
$\displaystyle=\Delta_{k}-\mu_{k+1}\left<{{\boldsymbol{\mathsf{g}}}_{k+1}},{{\boldsymbol{\mathsf{x}}}_{k}-\widehat{\boldsymbol{\mathsf{x}}}}\right>+\frac{G_{p^{\ast}}}{p^{\ast}}\mu_{k+1}^{p^{\ast}}\|{\boldsymbol{\mathsf{g}}}_{k+1}\|_{{\cal
X}^{\ast}}^{p^{\ast}}.$

∎

Lemma 3.6 allows showing that the sequence of Bregman distances
$({\mathbf{B}}_{p}({\boldsymbol{\mathsf{x}}}_{k},\widehat{\boldsymbol{\mathsf{x}}}))_{k\in\mathbb{N}}$
forms an almost super-martingale (in the Robbins-Siegmund sense defined below)
for $\widehat{\boldsymbol{\mathsf{x}}}\in{\cal X}_{\min}$ and well chosen
step-sizes $(\mu_{k})_{k\in\mathbb{N}}$. We will show the almost sure
convergence of the iterates using Robbins-Siegmund theorem.

######  Theorem 3.7 (Robbins-Siegmund theorem on the convergence of almost
super-martingales, [38, Lemma 11]).

Consider a filtration $({\cal F}_{k})_{k\in\mathbb{N}}$ and four non-negative,
$({\cal F}_{k})_{k\in\mathbb{N}}$ adapted processes
$(\alpha_{k})_{k\in\mathbb{N}}$, $(\beta_{k})_{k\in\mathbb{N}}$,
$(\gamma_{k})_{k\in\mathbb{N}}$, and $(\delta_{k})_{k\in\mathbb{N}}$. Let the
sequence $(\alpha_{k})_{k\in\mathbb{N}}$ be an _almost super-martingale_ ,
i.e. for all $k$ we have
$\mathbb{E}_{k}[\alpha_{k+1}]\leq(1+\beta_{k})\alpha_{k}+\gamma_{k}-\delta_{k}.$
Then the sequence $(\alpha_{k})_{k\in\mathbb{N}}$ converges a.s. to a random
variable $\alpha_{\infty}$, and $\sum_{k=1}^{\infty}\delta_{k}<\infty$ a.s. on
the set
$\\{\sum_{k=1}^{\infty}\beta_{k}<\infty,\,\sum_{k=1}^{\infty}\gamma_{k}<\infty\\}$.

Under certain conditions on ${\boldsymbol{\mathsf{x}}}_{0}$, the limit is the
MNS ${\boldsymbol{\mathsf{x}}}^{\dagger}$. Below $\mathbb{E}_{k}$ denotes the
conditional expectation with respect to the filtration ${\cal F}_{k}$.

######  Theorem 3.8.

Let $(\mu_{k})_{k\in\mathbb{N}}$ satisfy $\sum_{k=1}^{\infty}\mu_{k}=\infty$
and $\sum_{k=1}^{\infty}\mu_{k}^{p^{\ast}}<\infty,$ Assumption 3.1 hold, and
${\boldsymbol{\mathsf{x}}}^{\dagger}$ be the MNS. Then the sequence
$({\boldsymbol{\mathsf{x}}}_{k})_{k\in\mathbb{N}}$ converges a.s. to a
solution of (1):

$\displaystyle\mathbb{P}\Big{(}\lim_{k\rightarrow\infty}\inf_{\widetilde{\boldsymbol{\mathsf{x}}}\in{\cal
X}_{\min}}\|{\boldsymbol{\mathsf{x}}}_{k}-\widetilde{\boldsymbol{\mathsf{x}}}\|_{{\cal
X}}=0\Big{)}=1.$

Moreover, if ${\cal J}_{p}^{{\cal
X}}({\boldsymbol{\mathsf{x}}}_{0})\in\overline{\operatorname{range}({\mathbf{A}}^{\ast})}$,
we have
$\lim_{k\rightarrow\infty}{\mathbf{B}}_{p}({\boldsymbol{\mathsf{x}}}_{k},{\boldsymbol{\mathsf{x}}}^{\dagger})=0$
a.s.

###### Proof.

By Lemma 3.5, we have
$\left<{\partial\Psi({\boldsymbol{\mathsf{x}}}_{k})},{{\boldsymbol{\mathsf{x}}}_{k}-{\boldsymbol{\mathsf{x}}}^{\dagger}}\right>=p\Psi({\boldsymbol{\mathsf{x}}}_{k}).$
Moreover,

$\displaystyle\|g({\boldsymbol{\mathsf{x}}},{\boldsymbol{\mathsf{y}}},i)\|_{{\cal
X}^{\ast}}=\|{\mathbf{A}}_{i}^{\ast}\jmath_{p}^{{\cal
Y}}({\mathbf{A}}_{i}{\boldsymbol{\mathsf{x}}}-{\boldsymbol{\mathsf{y}}}_{i})\|_{{\cal
X}^{\ast}}\leq\|{\mathbf{A}}_{i}\|\|\jmath_{p}^{{\cal
Y}}({\mathbf{A}}_{i}{\boldsymbol{\mathsf{x}}}-{\boldsymbol{\mathsf{y}}}_{i})\|_{{\cal
Y}^{\ast}}\leq
L_{\max}\|{\mathbf{A}}_{i}{\boldsymbol{\mathsf{x}}}-{\boldsymbol{\mathsf{y}}}_{i}\|_{{\cal
Y}}^{p-1},$

with $L_{\max}=\max_{i\in[N]}\|{\mathbf{A}}_{i}\|$. Thus, since
$p^{\ast}(p-1)=p$, we have

$\displaystyle\mathbb{E}\big{[}\|g({\boldsymbol{\mathsf{x}}},{\boldsymbol{\mathsf{y}}},i)\|_{{\cal
X}^{\ast}}^{p^{\ast}}\big{]}\leq
pL_{\max}^{p^{\ast}}\frac{1}{N}\sum_{i=1}^{N}\frac{1}{p}\|{\mathbf{A}}_{i}{\boldsymbol{\mathsf{x}}}-{\boldsymbol{\mathsf{y}}}_{i}\|_{{\cal
Y}}^{p}=pL_{\max}^{p^{\ast}}\Psi({\boldsymbol{\mathsf{x}}}).$

Upon taking the conditional expectation $\mathbb{E}_{k}[\cdot]$ of the descent
property (13) (with
$\widehat{\boldsymbol{\mathsf{x}}}={\boldsymbol{\mathsf{x}}}^{\dagger}$), and
using the measurability of ${\boldsymbol{\mathsf{x}}}_{k}$ with respect to
${\cal F}_{k}$, we deduce

$\displaystyle\mathbb{E}_{k}[\Delta_{k+1}]$
$\displaystyle\leq\Delta_{k}-p\mu_{k+1}\Psi({\boldsymbol{\mathsf{x}}}_{k})+pL_{\max}^{p^{\ast}}\frac{G_{p^{\ast}}}{p^{\ast}}\mu_{k+1}^{p^{\ast}}\Psi({\boldsymbol{\mathsf{x}}}_{k}).$

Using Lemma 3.5 again we have
$\Psi({\boldsymbol{\mathsf{x}}}_{k})\leq\frac{L_{\max}^{p}}{C_{p}}\Delta_{k}$,
which yields

$\displaystyle\mathbb{E}_{k}[\Delta_{k+1}]$
$\displaystyle\leq\left(1+L_{\max}^{p^{\ast}+p}\frac{p}{C_{p}}\frac{G_{p^{\ast}}}{p^{\ast}}\mu_{k+1}^{p^{\ast}}\right)\Delta_{k}-p\mu_{k+1}\Psi({\boldsymbol{\mathsf{x}}}_{k}).$

Since $\sum_{k=1}^{\infty}\mu_{k}^{p^{\ast}}<\infty$ by assumption, we can
apply Theorem 3.7 and deduce that the sequence $(\Delta_{k})_{k\in\mathbb{N}}$
converges a.s. to a random variable $\Delta_{\infty}$ and
$\sum_{k=0}^{\infty}\mu_{k+1}\Psi({\boldsymbol{\mathsf{x}}}_{k})<\infty$ a.s.
Let $\Omega$ be the measurable set on which $(\Delta_{k})_{k\in\mathbb{N}}$
converges,
$\sum_{k=0}^{\infty}\mu_{k+1}\Psi({\boldsymbol{\mathsf{x}}}_{k})<\infty$, and
$\mathbb{P}(\Omega)=1$. Next we show
$\liminf_{k}\Psi({\boldsymbol{\mathsf{x}}}_{k})=0$ a.s. Consider an event
$\omega$ on which this is not the case, i.e. where
$\liminf_{k}\Psi({\boldsymbol{\mathsf{x}}}_{k})>0$. Then there exist
$\epsilon>0$ and $k_{\epsilon}\in\mathbb{N}$ such that for all $k\geq
k_{\epsilon}$, $\Psi({\boldsymbol{\mathsf{x}}}_{k})\geq\epsilon$, giving
$\sum_{k\geq
k_{\epsilon}}\mu_{k+1}\Psi({\boldsymbol{\mathsf{x}}}_{k})\geq\epsilon\sum_{k\geq
k_{\epsilon}}\mu_{k+1}$. Since for all events in $\Omega$ this would lead to a
contradiction: the right hand side diverges
($\sum_{k=1}^{\infty}\mu_{k}=\infty$ by assumption), whereas the left hand
side is the remainder of a convergent series, we conclude
$\omega\not\in\Omega$. Since $\mathbb{P}(\Omega^{c})=0$, we have
$\liminf_{k}\Psi({\boldsymbol{\mathsf{x}}}_{k})=0$ a.s. For every event in the
set where $\liminf_{k}\Psi({\boldsymbol{\mathsf{x}}}_{k})=0$ holds we can then
find a sub-sequence $({\boldsymbol{\mathsf{x}}}_{n_{k}})_{k\in\mathbb{N}}$
such that
$\lim_{k\rightarrow\infty}\Psi({\boldsymbol{\mathsf{x}}}_{n_{k}})=0$. Define
also
$\widehat{\Psi}({\boldsymbol{\mathsf{x}}})=\sum_{i=1}^{N}\widehat{\Psi}_{i}({\boldsymbol{\mathsf{x}}})$,
with
$\widehat{\Psi}_{i}({\boldsymbol{\mathsf{x}}})=\|{\mathbf{A}}_{i}{\boldsymbol{\mathsf{x}}}-{\boldsymbol{\mathsf{y}}}_{i}\|_{{\cal
Y}}$. We have $\liminf_{k}\widehat{\Psi}({\boldsymbol{\mathsf{x}}}_{k})=0$ and
$\lim_{j\rightarrow\infty}\widehat{\Psi}({\boldsymbol{\mathsf{x}}}_{n_{j}})=0$
(on the same subsequence), since by Young’s inequality,

$\displaystyle\Big{(}\sum_{i=1}^{N}\|{\mathbf{A}}_{i}{\boldsymbol{\mathsf{x}}}-{\boldsymbol{\mathsf{y}}}_{i}\|_{{\cal
Y}}^{p}\Big{)}^{1/p}\leq\sum_{i=1}^{N}\|{\mathbf{A}}_{i}{\boldsymbol{\mathsf{x}}}-{\boldsymbol{\mathsf{y}}}_{i}\|_{{\cal
Y}}\leq
N\Big{(}\frac{1}{N}\sum_{i=1}^{N}\|{\mathbf{A}}_{i}{\boldsymbol{\mathsf{x}}}-{\boldsymbol{\mathsf{y}}}_{i}\|_{{\cal
Y}}^{p}\Big{)}^{1/p}.$

Moreover, $\widehat{\Psi}({\boldsymbol{\mathsf{x}}})^{p}\leq
pN^{p}\Psi({\boldsymbol{\mathsf{x}}})$. The following argument is understood
pointwise on the a.s. set $\Omega$ where $(\Delta_{k})_{k\in\mathbb{N}}$
converges,
$\sum_{k=0}^{\infty}\mu_{k+1}\Psi({\boldsymbol{\mathsf{x}}}_{k})<\infty$, and
$\liminf_{k}\widehat{\Psi}({\boldsymbol{\mathsf{x}}}_{k})=0$. Since
$(\Delta_{k})_{k\in\mathbb{N}}$ converges it is bounded. By the coercivity of
the Bregman distance (see Lemma A.3) so are
$({\boldsymbol{\mathsf{x}}}_{k})_{k\in\mathbb{N}}$ and $({\cal J}_{p}^{{\cal
X}}({\boldsymbol{\mathsf{x}}}_{k}))_{k\in\mathbb{N}}$. By further passing to a
subsequence, we can find a subsequence of
$({\boldsymbol{\mathsf{x}}}_{n_{k}})_{k\in\mathbb{N}}$, that we denote the
same, such that $(\|{\boldsymbol{\mathsf{x}}}_{n_{k}}\|_{{\cal
X}})_{k\in\mathbb{N}}$ is convergent, $({\cal J}_{p}^{{\cal
X}}({\boldsymbol{\mathsf{x}}}_{n_{k}}))_{k\in\mathbb{N}}$ is weakly
convergent, and

$\displaystyle\lim_{k\rightarrow\infty}\widehat{\Psi}({\boldsymbol{\mathsf{x}}}_{n_{k}})=0\quad\text{and}\quad\widehat{\Psi}({\boldsymbol{\mathsf{x}}}_{n_{k}})\leq\widehat{\Psi}({\boldsymbol{\mathsf{x}}}_{n})\text{
for all }n<n_{k}.$ (14)

The latter can be obtained by setting $n_{1}=1$, and then recursively defining
$n_{k+1}=\min\\{k>n_{k}:\Psi({\boldsymbol{\mathsf{x}}}_{k})\leq\Psi({\boldsymbol{\mathsf{x}}}_{n_{k}})/2\\}$,
for $k\in\mathbb{N}$. Any following subsequence satisfies the same property.
Using Theorem 2.6(ii), we have for $k>\ell$

$\displaystyle{\mathbf{B}}_{p}({\boldsymbol{\mathsf{x}}}_{n_{\ell}},{\boldsymbol{\mathsf{x}}}_{n_{k}})$
$\displaystyle\\!=\\!\frac{1}{p^{\ast}}\\!\Big{(}\\!\|{\boldsymbol{\mathsf{x}}}_{n_{\ell}}\|_{{\cal
X}}^{p}\\!-\\!\|{\boldsymbol{\mathsf{x}}}_{n_{k}}\|_{{\cal
X}}^{p}\\!\Big{)}\\!+\\!\left<{\\!{\cal J}_{p}^{{\cal
X}}({\boldsymbol{\mathsf{x}}}_{n_{k}}\\!)\\!-\\!{\cal J}_{p}^{{\cal
X}}({\boldsymbol{\mathsf{x}}}_{n_{\ell}}\\!)},{{\boldsymbol{\mathsf{x}}}^{\dagger}\\!}\right>\\!+\\!\left<{\\!{\cal
J}_{p}^{{\cal X}}({\boldsymbol{\mathsf{x}}}_{n_{k}}\\!)\\!-\\!{\cal
J}_{p}^{{\cal
X}}({\boldsymbol{\mathsf{x}}}_{n_{\ell}}\\!)},{{\boldsymbol{\mathsf{x}}}_{n_{k}}\\!\\!-\\!{\boldsymbol{\mathsf{x}}}^{\dagger}\\!}\right>.$

Since the first two terms involve Cauchy sequences, it suffices to treat the
last term, denoted by ${\rm I}_{k,\ell}$. Using telescopic sum and applying
the iterate update rule, we have

$\displaystyle{\rm I}_{k,\ell}$
$\displaystyle=\sum_{n=n_{\ell}}^{n_{k}-1}\\!\left<{{\cal J}_{p}^{{\cal
X}}({\boldsymbol{\mathsf{x}}}_{n+1})\\!-\\!{\cal J}_{p}^{{\cal
X}}({\boldsymbol{\mathsf{x}}}_{n})},{{\boldsymbol{\mathsf{x}}}_{n_{k}}\\!-\\!{\boldsymbol{\mathsf{x}}}^{\dagger}}\right>=\\!\sum_{n=n_{\ell}}^{n_{k}-1}\\!\mu_{n+1}\left<{{\mathbf{A}}_{i_{n+1}}^{\ast}\jmath_{p}^{{\cal
Y}}({\mathbf{A}}_{i_{k+1}}{\boldsymbol{\mathsf{x}}}_{n}-{\boldsymbol{\mathsf{y}}}_{i_{n+1}})},{{\boldsymbol{\mathsf{x}}}_{n_{k}}-{\boldsymbol{\mathsf{x}}}^{\dagger}}\right>$
$\displaystyle=\sum_{n=n_{\ell}}^{n_{k}-1}\mu_{n+1}\left<{\jmath_{p}^{{\cal
Y}}({\mathbf{A}}_{i_{n+1}}{\boldsymbol{\mathsf{x}}}_{n}-{\boldsymbol{\mathsf{y}}}_{i_{k+1}})},{{\mathbf{A}}_{i_{n+1}}{\boldsymbol{\mathsf{x}}}_{n_{k}}-{\boldsymbol{\mathsf{y}}}_{i_{n+1}}}\right>.$

By the Cauchy-Schwarz inequality and properties of the duality map, we get

$\displaystyle|{\rm I}_{k,\ell}|$
$\displaystyle\leq\sum_{n=n_{\ell}}^{n_{k}-1}\\!\\!\mu_{n+1}\\!\|{\mathbf{A}}_{i_{n+1}}\\!{\boldsymbol{\mathsf{x}}}_{n}\\!-\\!{\boldsymbol{\mathsf{y}}}_{i_{n+1}}\|_{{\cal
Y}}^{p-1}\|{\mathbf{A}}_{i_{n+1}}{\boldsymbol{\mathsf{x}}}_{n_{k}}\\!-\\!{\boldsymbol{\mathsf{y}}}_{i_{n+1}}\\!\|_{{\cal
Y}}\leq\\!\sum_{n=n_{\ell}}^{n_{k}-1}\\!\mu_{n+1}\widehat{\Psi}_{i_{n+1}}({\boldsymbol{\mathsf{x}}}_{n})^{p-1}\widehat{\Psi}_{i_{n+1}}({\boldsymbol{\mathsf{x}}}_{n_{k}}).$

Since
$\widehat{\Psi}_{i}({\boldsymbol{\mathsf{x}}})\leq\widehat{\Psi}({\boldsymbol{\mathsf{x}}})$,
for all $i\in[N]$, we use (14) and get

$\displaystyle|{\rm I}_{k,\ell}|$
$\displaystyle\leq\sum_{n=n_{\ell}}^{n_{k}-1}\mu_{n+1}\widehat{\Psi}({\boldsymbol{\mathsf{x}}}_{n})^{p-1}\widehat{\Psi}({\boldsymbol{\mathsf{x}}}_{n_{k}})\leq\sum_{n=n_{\ell}}^{n_{k}-1}\mu_{n+1}\widehat{\Psi}({\boldsymbol{\mathsf{x}}}_{n})^{p}.$

Since $\widehat{\Psi}({\boldsymbol{\mathsf{x}}})^{p}\leq
pN^{p}\Psi({\boldsymbol{\mathsf{x}}})$, the right hand side of the inequality
converges to $0$ as $n_{\ell}\to\infty$. Therefore, by [44, Theorem 2.12(e)],
it follows that $({\boldsymbol{\mathsf{x}}}_{n_{k}})_{k\in\mathbb{N}}$, is a
Cauchy sequence, and thus converges strongly to an
$\widehat{\boldsymbol{\mathsf{x}}}$ such that
$\Psi(\widehat{\boldsymbol{\mathsf{x}}})=0$.

The above argument showing the a.s. convergence of
$(\Delta_{k})_{k\in\mathbb{N}}$ can be applied pointwise to any solution.
Namely, on the event where
$({\boldsymbol{\mathsf{x}}}_{n_{k}})_{k\in\mathbb{N}}$ converges strongly to
an $\widehat{\boldsymbol{\mathsf{x}}}\in\mathcal{X}_{\min}$ (i.e.
${\mathbf{A}}\widehat{\boldsymbol{\mathsf{x}}}={\boldsymbol{\mathsf{y}}}$),
define
$\widehat{\Delta}_{k}:={\mathbf{B}}_{p}({\boldsymbol{\mathsf{x}}}_{k},\widehat{\boldsymbol{\mathsf{x}}})$.
By repeating the argument using Lemma 3.5, we deduce

$\displaystyle\widehat{\Delta}_{k+1}\leq\left(1+L_{\max}^{p^{\ast}+p}\frac{p}{C_{p}}\frac{G_{p^{\ast}}}{p^{\ast}}\mu_{k+1}^{p^{\ast}}\right)\widehat{\Delta}_{k}-p\mu_{k+1}\Psi_{i_{k}}({\boldsymbol{\mathsf{x}}}_{k}).$

Since $\sum_{k=1}^{\infty}\mu_{k}^{p^{\ast}}<\infty$, it follows that the
(deterministic) sequence $(\widehat{\Delta}_{k})_{k\in\mathbb{N}}$ converges
to a $\widehat{\Delta}_{\infty}\geq 0$. The continuity of the Bregman distance
in the first argument (Theorem 2.6(vi)) gives
$\lim_{j\rightarrow\infty}{\mathbf{B}}_{p}({\boldsymbol{\mathsf{x}}}_{n_{j}},\widehat{\boldsymbol{\mathsf{x}}})={\mathbf{B}}_{p}(\widehat{\boldsymbol{\mathsf{x}}},\widehat{\boldsymbol{\mathsf{x}}})=0$,
and thus $\widehat{\Delta}_{\infty}=0$. Moreover, by the $p$-convexity of
${\cal X}$ (Theorem 2.6(iii)), we have
$0\leq\|{\boldsymbol{\mathsf{x}}}_{k}-\widehat{\boldsymbol{\mathsf{x}}}\|_{{\cal
X}}^{p}\leq\frac{p}{C_{p}}\widehat{\Delta}_{k}.$ From the squeeze theorem it
follows that
$\lim_{k\rightarrow\infty}\|{\boldsymbol{\mathsf{x}}}_{k}-\widehat{\boldsymbol{\mathsf{x}}}\|_{{\cal
X}}=0$. Thus, for every event in an a.s. set $\Omega$, the sequence
$({\boldsymbol{\mathsf{x}}}_{k})_{k\in\mathbb{N}}$ strongly converge to some
minimising solution, that is

$\displaystyle\mathbb{P}\Big{(}\lim_{k\rightarrow\infty}\inf_{\widetilde{\boldsymbol{\mathsf{x}}}\in{\cal
X}_{\min}}\|{\boldsymbol{\mathsf{x}}}_{k}-\widetilde{\boldsymbol{\mathsf{x}}}\|_{{\cal
X}}=0\Big{)}=1.$

Next assume ${\cal J}_{p}^{{\cal
X}}({\boldsymbol{\mathsf{x}}}_{0})\in\overline{\mathrm{range}({\mathbf{A}}^{\ast})}$.
From (10), it follows that ${\cal J}_{p}^{{\cal
X}}({\boldsymbol{\mathsf{x}}}_{k})\in\overline{\operatorname{range}({\mathbf{A}}^{\ast})}$
holds for all $k\geq 1$. By the continuity of ${\cal J}_{p}^{{\cal X}}$, we
have ${\cal J}_{p}^{{\cal
X}}(\widehat{\boldsymbol{\mathsf{x}}})\in\overline{\operatorname{range}({\mathbf{A}}^{\ast})}$.
Thus, from
${\mathbf{A}}(\widehat{\boldsymbol{\mathsf{x}}}-{\boldsymbol{\mathsf{x}}}^{\dagger})=0$
and Lemma 3.3 it follows
$\widehat{\boldsymbol{\mathsf{x}}}={\boldsymbol{\mathsf{x}}}^{\dagger}$. ∎

The assumptions and conclusions of Theorem 3.8 can be broken down into two
parts. The step-size conditions $\sum_{k=1}^{\infty}\mu_{k}=\infty$ and
$\sum_{k=1}^{\infty}\mu_{k}^{p^{\ast}}<\infty$ are required to show the a.s.
convergence of
$({\mathbf{B}}_{p}({\boldsymbol{\mathsf{x}}}_{k},\widehat{\boldsymbol{\mathsf{x}}}))_{k\in\mathbb{N}}$
to $0$, for some non-deterministic
$\widehat{\boldsymbol{\mathsf{x}}}\in\mathcal{X}_{\min}$. The remaining
assumption ${\cal J}_{p}^{{\cal
X}}({\boldsymbol{\mathsf{x}}}_{0})\in\overline{\mathrm{range}({\mathbf{A}}^{\ast})}$
is needed to identify this limit as the MNS
${\boldsymbol{\mathsf{x}}}^{\dagger}$, as the Landweber method [44, Remark
3.12]. If ${\cal J}_{p}^{{\cal
X}}({\boldsymbol{\mathsf{x}}}_{0})\not\in\overline{\mathrm{range}({\mathbf{A}}^{\ast})}$,
we commonly establish convergence to an MNS relative to
${\boldsymbol{\mathsf{x}}}_{0}$, i.e. a solution which minimises
$\|{\boldsymbol{\mathsf{x}}}-{\boldsymbol{\mathsf{x}}}_{0}\|_{{\cal X}}$,
analogous to the Euclidean case [24].

######  Remark 3.9.

The stepsize conditions $\sum_{k=1}^{\infty}\mu_{k}=\infty$ and
$\sum_{k=1}^{\infty}\mu_{k}^{p^{\ast}}<\infty$ are satisfied by a polynomially
decaying step-size schedule
$(\mu_{k})_{k\in\mathbb{N}}=(\mu_{0}k^{-\beta})_{k\in\mathbb{N}}$, with
$\frac{1}{p^{\ast}}<\beta\leq 1$.

Theorem 3.8 states sufficient conditions ensuring the a.s. convergence of
$({\mathbf{B}}_{p}({\boldsymbol{\mathsf{x}}}_{k},{\boldsymbol{\mathsf{x}}}^{\dagger}))_{k\in\mathbb{N}}$
to $0$. To strengthen this to the convergence in expectation, we require an
additional assumption to ensure that
$({\mathbf{B}}_{p}({\boldsymbol{\mathsf{x}}}_{k},{\boldsymbol{\mathsf{x}}}^{\dagger}))_{k\in\mathbb{N}}$
is a uniformly integrable super-martingale and the space ${\cal X}$ being
uniformly smooth. Note that removing the assumptions of Theorem 3.8 from
Theorem 3.10 would still result in convergence in expectation to some non-
negative random variable, but not necessarily to $0$. Recall that a family
$(X_{t})_{t}$ of random variables is uniformly integrable provided
$\lim_{k\rightarrow\infty}\sup_{t}\mathbb{E}[\|X_{t}\|\mid\boldsymbol{1}_{\|X_{t}\|\geq
k}]=0$, where $\boldsymbol{1}(\cdot)$ is the indicator function.

######  Theorem 3.10.

Let the conditions of Theorem 3.8 hold with ${\cal J}_{p}^{{\cal
X}}({\boldsymbol{\mathsf{x}}}_{0})\in\overline{\operatorname{range}({\mathbf{A}}^{\ast})}$
and let
$\mu_{k}^{p^{\ast}-1}\leq\frac{p^{\ast}}{G_{p^{\ast}}L_{\max}^{p^{\ast}}}$ for
all $k\in\mathbb{N}$. Then there holds
$\lim_{k\rightarrow\infty}\mathbb{E}[{\mathbf{B}}_{p}({\boldsymbol{\mathsf{x}}}_{k},{\boldsymbol{\mathsf{x}}}^{\dagger})]=0.$
Moreover, for $1\leq r\leq p$, we have
$\lim_{k\rightarrow\infty}\mathbb{E}[\|{\boldsymbol{\mathsf{x}}}_{k}-{\boldsymbol{\mathsf{x}}}^{\dagger}\|_{{\cal
X}}^{r}]=0$ and if ${\cal X}$ is additionally uniformly smooth, then
$\lim_{k\rightarrow\infty}\mathbb{E}[\|{\cal J}_{p}^{{\cal
X}}({\boldsymbol{\mathsf{x}}}_{k})-{\cal J}_{p}^{{\cal
X}}({\boldsymbol{\mathsf{x}}}^{\dagger})\|_{{\cal X}^{\ast}}^{p^{\ast}}]=0$.

###### Proof.

The step-size conditions allow applying Lemma A.2, which yields
${\mathbf{B}}_{p}({\boldsymbol{\mathsf{x}}}_{k},{\boldsymbol{\mathsf{x}}}^{\dagger})\leq{\mathbf{B}}_{p}({\boldsymbol{\mathsf{x}}}_{0},{\boldsymbol{\mathsf{x}}}^{\dagger})$
for all $k$. It follows that
$({\mathbf{B}}_{p}({\boldsymbol{\mathsf{x}}}_{k},{\boldsymbol{\mathsf{x}}}^{\dagger}))_{k\in\mathbb{N}}$
is bounded, and is thus uniformly integrable, and by Theorem 3.8, it converges
a.s. to $0$. Then, by Vitali’s convergence theorem [1, Theorem 4.5.4], we
deduce that $(\Delta_{k})_{k\in\mathbb{N}}$ converges to $0$ in expectation as
well. Using now the $p$-convexity of ${\cal X}$ and the monotonicity of
expectation, we have

$\displaystyle
0\leq\frac{C_{p}}{p}\lim_{k\rightarrow\infty}\mathbb{E}[\|{\boldsymbol{\mathsf{x}}}_{k}-{\boldsymbol{\mathsf{x}}}^{\dagger}\|_{{\cal
X}}^{p}]\leq\lim_{k\rightarrow\infty}\mathbb{E}[{\mathbf{B}}_{p}({\boldsymbol{\mathsf{x}}}_{k},{\boldsymbol{\mathsf{x}}}^{\dagger})]=0.$

By the continuity of the power function and the Lyapunov inequality for $1\leq
r{\leq p}$, we have

$\displaystyle
0\leq\lim_{k\rightarrow\infty}\mathbb{E}[\|{\boldsymbol{\mathsf{x}}}_{k}-{\boldsymbol{\mathsf{x}}}^{\dagger}\|_{{\cal
X}}^{r}]\leq\lim_{k\rightarrow\infty}(\mathbb{E}[\|{\boldsymbol{\mathsf{x}}}_{k}-{\boldsymbol{\mathsf{x}}}^{\dagger}\|_{{\cal
X}}^{p}])^{r/p}=0.$

To prove the last claim we use uniform smoothness of ${\cal X}$ and Theorem
2.3(iv), to deduce

$\displaystyle\|{\cal J}_{p}^{{\cal X}}({\boldsymbol{\mathsf{x}}}_{k})-{\cal
J}_{p}^{{\cal X}}({\boldsymbol{\mathsf{x}}}^{\dagger})\|_{{\cal
X}^{\ast}}^{p^{\ast}}\leq C\max\\{1,\|{\boldsymbol{\mathsf{x}}}_{k}\|_{{\cal
X}},\|{\boldsymbol{\mathsf{x}}}^{\dagger}\|_{{\cal
X}}\\}^{p}\,{\overline{\rho}_{\cal
X}}(\|{\boldsymbol{\mathsf{x}}}_{k}-{\boldsymbol{\mathsf{x}}}^{\dagger}\|_{{\cal
X}})^{p^{\ast}},$

where $\overline{\rho}_{\cal X}(\tau)=\rho_{\cal X}(\tau)/\tau$ is a modulus
of smoothness function such that $\overline{\rho}(\tau)\leq 1$ and
$\lim_{\tau\rightarrow 0}\overline{\rho}(\tau)=0$, cf. Definition 2.2. By
Lemmas A.2 and A.3 $(\|{\boldsymbol{\mathsf{x}}}_{k}\|_{{\cal
X}}^{p})_{k\in\mathbb{N}}$ is (uniformly) bounded, giving that the sequence
$(\|{\cal J}_{p}^{{\cal X}}({\boldsymbol{\mathsf{x}}}_{k})-{\cal J}_{p}^{{\cal
X}}({\boldsymbol{\mathsf{x}}}^{\dagger})\|_{{\cal
X}^{\ast}}^{p^{\ast}})_{k\in\mathbb{N}}$ is bounded and thus uniformly
integrable. Since
$\lim_{k\rightarrow\infty}\mathbb{E}[\|{\boldsymbol{\mathsf{x}}}_{k}-{\boldsymbol{\mathsf{x}}}^{\dagger}\|_{{\cal
X}}]=0$, it follows that
$\|{\boldsymbol{\mathsf{x}}}_{k}-{\boldsymbol{\mathsf{x}}}^{\dagger}\|_{{\cal
X}}$ converges to $0$ in probability, and thus by the continuous mapping
theorem ${\overline{\rho}_{\cal
X}}(\|{\boldsymbol{\mathsf{x}}}_{k}-{\boldsymbol{\mathsf{x}}}^{\dagger}\|_{{\cal
X}})^{p^{\ast}}$ also converges to $0$ in probability. Applying Vitaly’s
theorem to the uniformly integrable sequence $(\|{\cal J}_{p}^{{\cal
X}}({\boldsymbol{\mathsf{x}}}_{k})-{\cal J}_{p}^{{\cal
X}}({\boldsymbol{\mathsf{x}}}^{\dagger})\|_{{\cal
X}^{\ast}}^{p^{\ast}})_{k\in\mathbb{N}}$ yields that it converges to $0$ in
measure, and the claim follows. ∎

######  Remark 3.11.

Note that the condition ${\cal J}_{p}^{{\cal
X}}({\boldsymbol{\mathsf{x}}}_{0})\in\overline{\operatorname{range}({\mathbf{A}}^{\ast})}$
on ${\boldsymbol{\mathsf{x}}}_{0}$ is crucial for ensuring that all the limits
are the same. Landweber iterations converge for uniformly convex and smooth
${\cal X}$, and any Banach space ${\cal Y}$ [44, Theorem 3.3]. In our
analysis, we have assumed that ${\cal X}$ is $p$-convex to simplify the
analysis. First, $p$-convexity is used in the proof of Lemma 3.6. If ${\cal
X}$ were only uniformly convex (and ${\cal X}^{\ast}$ only uniformly smooth),
then we may use the modulus of smoothness function $\rho_{\mathcal{X}}$, cf.
(2.2) and [46, Theorem 2.41], to establish a suitable analogue of the descent
property (13). Second, $p$-convexity is used in the proof of Theorem 3.8,
allowing a more direct application of Robbins-Siegmund theorem by relating the
objective values to Bregman distances. Meanwhile, the Landweber method in [44]
requires step-sizes that depend on the modulus of smoothness, the current
iterate and objective value, which is more restrictive than that in this work.

###  3.2 Convergence analysis for the generalised Kaczmarz model

Schöpfer et al [44] studied general powers of the Banach space norm and sub-
gradients of the form
$\partial(\frac{1}{q}\|{\mathbf{A}}\cdot-{\boldsymbol{\mathsf{y}}}\|_{{\cal
Y}}^{q})({\boldsymbol{\mathsf{x}}})$. Now we take an analogous perspective for
the objective

$\Psi({\boldsymbol{\mathsf{x}}})=\frac{1}{N}\sum_{i=1}^{N}\Psi_{i}({\boldsymbol{\mathsf{x}}}),\quad\mbox{with
}\Psi_{i}({\boldsymbol{\mathsf{x}}}):=\frac{1}{q}\|{\mathbf{A}}_{i}{\boldsymbol{\mathsf{x}}}-{\boldsymbol{\mathsf{y}}}_{i}\|_{{\cal
Y}}^{q},$

with $1<q\leq 2$. This model is herein called the generalised Kaczmarz model.
(Note that this is different from the randomised extended Kaczmarz method
[53].) We shall show the convergence of SGD with stochastic directions

$\displaystyle
g({\boldsymbol{\mathsf{x}}},{\boldsymbol{\mathsf{y}}},i)={\mathbf{A}}_{i}^{\ast}\jmath_{q}^{{\cal
Y}}({\mathbf{A}}_{i}{\boldsymbol{\mathsf{x}}}-{\boldsymbol{\mathsf{y}}}_{i})=\partial(\tfrac{1}{q}\|{\mathbf{A}}_{i}\cdot-{\boldsymbol{\mathsf{y}}}_{i}\|_{{\cal
Y}}^{q})({\boldsymbol{\mathsf{x}}}).$ (15)

The descent property (13) is unaffected, and a direct computation again yields

${\mathbf{B}}_{p}({\boldsymbol{\mathsf{x}}}_{k+1},{\boldsymbol{\mathsf{x}}}^{\dagger})\leq{\mathbf{B}}_{p}({\boldsymbol{\mathsf{x}}}_{k},{\boldsymbol{\mathsf{x}}}^{\dagger})-\mu_{k+1}\left<{{\boldsymbol{\mathsf{g}}}_{k+1}},{{\boldsymbol{\mathsf{x}}}_{k}-{\boldsymbol{\mathsf{x}}}^{\dagger}}\right>+\frac{G_{p^{\ast}}}{p^{\ast}}\mu_{k+1}^{p^{\ast}}\|{\boldsymbol{\mathsf{g}}}_{k+1}\|_{{\cal
X}^{\ast}}^{p^{\ast}}.$ (16)

However, Robbins-Siegmund theorem cannot be applied directly. Instead, we
pursue a different proof strategy by first establishing the uniform
boundedness of iterates.

######  Lemma 3.12.

Let Assumption 3.1 hold. Consider SGD with descent directions (15) for
$1<q\leq 2$, and assume that
$\mu_{k}^{p^{\ast}-1}<\frac{p^{\ast}}{G_{p^{\ast}}L_{\max}^{p^{\ast}}}$ holds
for all $k\in\mathbb{N}$ and
$\sum_{k=1}^{\infty}\mu_{k}^{p^{\ast}}=:\Gamma<\infty$. Then
$({\mathbf{B}}_{p}({\boldsymbol{\mathsf{x}}}_{k},{\boldsymbol{\mathsf{x}}}^{\dagger}))_{k\in\mathbb{N}}$
and $({\boldsymbol{\mathsf{x}}}_{k})_{k\in\mathbb{N}}$ are uniformly bounded.

###### Proof.

Let
$\overline{\Psi}_{i}({\boldsymbol{\mathsf{x}}})=\|{\mathbf{A}}_{i}{\boldsymbol{\mathsf{x}}}-{\boldsymbol{\mathsf{y}}}_{i}\|_{{\cal
Y}}^{q}$, and
$\Delta_{k}={\mathbf{B}}_{p}({\boldsymbol{\mathsf{x}}}_{k},{\boldsymbol{\mathsf{x}}}^{\dagger})$.
Then we have
$\left<{{\boldsymbol{\mathsf{g}}}_{k+1}},{{\boldsymbol{\mathsf{x}}}_{k}-{\boldsymbol{\mathsf{x}}}^{\dagger}}\right>=\overline{\Psi}_{i_{k+1}}({\boldsymbol{\mathsf{x}}}_{k})$
and

$\displaystyle\|{\boldsymbol{\mathsf{g}}}_{k+1}\|_{{\cal
X}^{\ast}}^{p^{\ast}}$
$\displaystyle=\|{\mathbf{A}}_{i_{k+1}}^{\ast}\jmath_{q}^{{\cal
Y}}({\mathbf{A}}_{i_{k+1}}{\boldsymbol{\mathsf{x}}}-{\boldsymbol{\mathsf{y}}}_{i_{k+1}})\|_{{\cal
X}^{\ast}}^{p^{\ast}}\leq
L_{\max}^{p^{\ast}}\|{\mathbf{A}}_{i_{k+1}}{\boldsymbol{\mathsf{x}}}-{\boldsymbol{\mathsf{y}}}_{i_{k+1}}\|_{{\cal
Y}}^{p^{\ast}(q-1)}$ $\displaystyle\leq
L_{\max}^{p^{\ast}}\overline{\Psi}_{i_{k+1}}({\boldsymbol{\mathsf{x}}}_{k})^{p^{\ast}\frac{q-1}{q}}=L_{\max}^{p^{\ast}}\overline{\Psi}_{i_{k+1}}({\boldsymbol{\mathsf{x}}}_{k})^{\frac{p^{\ast}}{q^{\ast}}},$

where $q^{\ast}\geq 2$ is the conjugate exponent of $q$. Plugging this into
(16) gives

$\displaystyle\Delta_{k+1}\leq\Delta_{k}-\mu_{k+1}\overline{\Psi}_{i_{k+1}}({\boldsymbol{\mathsf{x}}}_{k})+L_{\max}^{p^{\ast}}\frac{G_{p^{\ast}}}{p^{\ast}}\mu_{k+1}^{p^{\ast}}\overline{\Psi}_{i_{k+1}}({\boldsymbol{\mathsf{x}}}_{k})^{\frac{p^{\ast}}{q^{\ast}}}.$
(17)

Since $1<p^{\ast}\leq 2$ by Theorem 2.6(iii), and $q^{\ast}\geq 2$, we have
$\frac{p^{\ast}}{q^{\ast}}\leq 1$. Now we define two sets of indices

${\cal I}=\\{j\leq
k:\overline{\Psi}_{i_{j+1}}({\boldsymbol{\mathsf{x}}}_{j})\geq 1\\}\text{ and
}{\cal J}=\\{j\leq
k:\overline{\Psi}_{i_{j+1}}({\boldsymbol{\mathsf{x}}}_{j})<1\\},$

so that ${\cal I}\cap{\cal J}=\emptyset$, and ${\cal I}\cup{\cal J}=[k]$. Note
that ${\cal I}$ and ${\cal J}$ actually depend on the current iterate index
$k$. Applying the inductive argument to (17) gives

$\displaystyle\Delta_{k+1}$
$\displaystyle\leq\Delta_{0}-\sum_{j=0}^{k}\mu_{j+1}\overline{\Psi}_{i_{j+1}}({\boldsymbol{\mathsf{x}}}_{j})+L_{\max}^{p^{\ast}}\frac{G_{p^{\ast}}}{p^{\ast}}\sum_{j=0}^{k}\mu_{j+1}^{p^{\ast}}\overline{\Psi}_{i_{j+1}}({\boldsymbol{\mathsf{x}}}_{j})^{\frac{p^{\ast}}{q^{\ast}}}$
$\displaystyle=\Delta_{0}\underbrace{-\sum_{j\in{\cal
I}}\mu_{j+1}\overline{\Psi}_{i_{j+1}}({\boldsymbol{\mathsf{x}}}_{j})+L_{\max}^{p^{\ast}}\frac{G_{p^{\ast}}}{p^{\ast}}\sum_{j\in{\cal
I}}\mu_{j+1}^{p^{\ast}}\overline{\Psi}_{i_{j+1}}({\boldsymbol{\mathsf{x}}}_{j})^{\frac{p^{\ast}}{q^{\ast}}}}_{(\star)}\underbrace{-\sum_{j\in{\cal
J}}\mu_{j+1}\overline{\Psi}_{i_{j+1}}({\boldsymbol{\mathsf{x}}}_{j})}_{(\star\star)}$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\quad+\underbrace{L_{\max}^{p^{\ast}}\frac{G_{p^{\ast}}}{p^{\ast}}\sum_{j\in{\cal
J}}\mu_{j+1}^{p^{\ast}}\overline{\Psi}_{i_{j+1}}({\boldsymbol{\mathsf{x}}}_{j})^{\frac{p^{\ast}}{q^{\ast}}}}_{(\star\star\star)}.$

Next we analyse these three terms separately. First, for $j\in{\cal I}$, we
have $\overline{\Psi}_{i_{j+1}}({\boldsymbol{\mathsf{x}}}_{j})\geq 1$ and
since $\frac{p^{\ast}}{q^{\ast}}<1$, we have
$\overline{\Psi}_{i_{j+1}}({\boldsymbol{\mathsf{x}}}_{j})^{\frac{p^{\ast}}{q^{\ast}}}\leq\overline{\Psi}_{i_{j+1}}({\boldsymbol{\mathsf{x}}}_{j})$,
giving

$\displaystyle(\star)$ $\displaystyle\\!\leq\\!-\\!\sum_{j\in{\cal
I}}\mu_{j+1}\\!\overline{\Psi}_{i_{j+1}}({\boldsymbol{\mathsf{x}}}_{j})\\!+\\!L_{\max}^{p^{\ast}}\\!\frac{G_{p^{\ast}}}{p^{\ast}}\\!\sum_{j\in{\cal
I}}\mu_{j+1}^{p^{\ast}}\overline{\Psi}_{i_{j+1}}({\boldsymbol{\mathsf{x}}}_{j})$
$\displaystyle=\\!-\\!\sum_{j\in{\cal
I}}\\!\Big{(}\\!1\\!-\\!L_{\max}^{p^{\ast}}\frac{G_{p^{\ast}}}{p^{\ast}}\mu_{j+1}^{p^{\ast}-1}\\!\Big{)}\\!\mu_{j+1}\\!\overline{\Psi}_{i_{j+1}}({\boldsymbol{\mathsf{x}}}_{j}).$

Since
$\mu_{j+1}^{p^{\ast}-1}<\frac{p^{\ast}}{G_{p^{\ast}}L_{\max}^{p^{\ast}}}$
holds by assumption, the term $(\star)$ is non-positive. Moreover,
$(\star\star)$ is trivially non-positive. Since
$\overline{\Psi}_{i_{j+1}}({\boldsymbol{\mathsf{x}}}_{j})<1$ for $j\in{\cal
J}$, the last term $(\star\star\star)$ can be bounded as

$\displaystyle L_{\max}^{p^{\ast}}\frac{G_{p^{\ast}}}{p^{\ast}}\sum_{j\in{\cal
J}}\mu_{j+1}^{p^{\ast}}\overline{\Psi}_{i_{j+1}}({\boldsymbol{\mathsf{x}}}_{j})^{\frac{p^{\ast}}{q^{\ast}}}\leq
L_{\max}^{p^{\ast}}\frac{G_{p^{\ast}}}{p^{\ast}}\sum_{j\in{\cal
J}}\mu_{j+1}^{p^{\ast}}\leq
L_{\max}^{p^{\ast}}\frac{G_{p^{\ast}}}{p^{\ast}}\sum_{j=1}^{\infty}\mu_{j}^{p^{\ast}}=L_{\max}^{p^{\ast}}\frac{G_{p^{\ast}}}{p^{\ast}}\Gamma.$

By combining the last three bounds on $(\star)$, $(\star\star)$ and
$(\star\star\star)$, we get

$\Delta_{k+1}\leq\Delta_{0}+L_{\max}^{p^{\ast}}\frac{G_{p^{\ast}}}{p^{\ast}}\Gamma,\text{
for all }k\geq 0.$

Thus, $(\Delta_{k})_{k\in\mathbb{N}}$ is uniformly bounded and by Lemma A.3,
so is $({\boldsymbol{\mathsf{x}}}_{k})_{k\in\mathbb{N}}$. ∎

The proof of Lemma 3.12 exposes the challenge in extending the convergence
results to general stochastic directions. Namely, in the proof of Theorem 3.8,
we showed the convergence by taking conditional expectation of (13), recasting
the resulting expression as an almost super-martingale, and then relating
objective values to Bregman distances via
$\Psi({\boldsymbol{\mathsf{x}}}_{k})\leq C\Delta_{k}$, for some $C>0$. Here,
using $\frac{q}{q^{\ast}}=q-1$ and $\frac{p^{\ast}}{p}=p^{\ast}-1$, we instead
have

$\Psi({\boldsymbol{\mathsf{x}}}_{k})^{\frac{p^{\ast}}{q^{\ast}}}\leq
C\Delta_{k}^{(p^{\ast}-1)(q-1)},\quad\mbox{with
}C=q^{-\frac{p^{\ast}}{q^{\ast}}}L_{\max}^{p^{\ast}(q-1)}\Big{(}\frac{p}{C_{p}}\Big{)}^{(p^{\ast}-1)(q-1)},$

which gives

$\displaystyle\mathbb{E}_{k}[\Delta_{k+1}]\leq\Delta_{k}+CL_{\max}^{p^{\ast}}q^{\frac{p^{\ast}}{q^{\ast}}}\frac{G_{p^{\ast}}}{p^{\ast}}\mu_{k+1}^{p^{\ast}}\Delta_{k}^{(p^{\ast}-1)(q-1)}-q\mu_{k+1}\Psi({\boldsymbol{\mathsf{x}}}_{k}).$

Here $0<(p^{\ast}-1)(q-1)<1$, provided $p^{\ast}\neq 2$ and $q\neq 2$.
Therefore, Robbins-Siegmund theorem cannot be applied directly. Nonetheless,
we still have the following analogue of Theorem 3.10.

######  Theorem 3.13.

Consider iterations (10) with descent directions (15) for $1<q\leq 2$ and let
Assumption 3.1 hold and ${\boldsymbol{\mathsf{x}}}^{\dagger}$ be the MNS. Let
the step-sizes $(\mu_{k})_{k\in\mathbb{N}}$ satisfy
$\sum_{k=1}^{\infty}\mu_{k}=\infty$,
$\sum_{k=1}^{\infty}\mu_{k}^{p^{\ast}}<\infty$, and
$\mu_{k}^{p^{\ast}-1}<\frac{p^{\ast}}{G_{p^{\ast}}L_{\max}^{p^{\ast}}}$ for
all $k\in\mathbb{N}$. Then the sequence
$({\boldsymbol{\mathsf{x}}}_{k})_{k\in\mathbb{N}}$ converges a.s. to a
solution of (1):

$\displaystyle\mathbb{P}\Big{(}\lim_{k\rightarrow\infty}\inf_{\widetilde{\boldsymbol{\mathsf{x}}}\in{\cal
X}_{\min}}\|{\boldsymbol{\mathsf{x}}}_{k}-\widetilde{\boldsymbol{\mathsf{x}}}\|_{{\cal
X}}=0\Big{)}=1.$

Moreover, if ${\cal J}_{p}^{{\cal
X}}({\boldsymbol{\mathsf{x}}}_{0})\in\overline{\mathrm{range}({\mathbf{A}}^{\ast})}$,
we have

$\lim_{k\rightarrow\infty}{\mathbf{B}}_{p}({\boldsymbol{\mathsf{x}}}_{k},{\boldsymbol{\mathsf{x}}}^{\dagger})=0\
\mbox{
a.s.}\quad\mbox{and}\quad\lim_{k\rightarrow\infty}\mathbb{E}[{\mathbf{B}}_{p}({\boldsymbol{\mathsf{x}}}_{k},{\boldsymbol{\mathsf{x}}}^{\dagger})]=0.$

###### Proof.

To establish the a.s. convergence of iterates, we first take the conditional
expectation of the descent property (16) and obtain

$\mathbb{E}_{k}[\Delta_{k+1}]\leq\Delta_{k}-\mu_{k+1}\left<{\mathbb{E}_{k}[{\boldsymbol{\mathsf{g}}}_{k+1}]},{{\boldsymbol{\mathsf{x}}}_{k}-{\boldsymbol{\mathsf{x}}}^{\dagger}}\right>+\frac{G_{p^{\ast}}}{p^{\ast}}\mu_{k+1}^{p^{\ast}}\mathbb{E}_{k}\big{[}\|{\boldsymbol{\mathsf{g}}}_{k+1}\|_{{\cal
X}^{\ast}}^{p^{\ast}}\big{]}.$ (18)

We now have
$\left<{\mathbb{E}_{k}[{\boldsymbol{\mathsf{g}}}_{k+1}]},{{\boldsymbol{\mathsf{x}}}_{k}-{\boldsymbol{\mathsf{x}}}^{\dagger}}\right>=\left<{\partial\Psi({\boldsymbol{\mathsf{x}}}_{k})},{{\boldsymbol{\mathsf{x}}}_{k}-{\boldsymbol{\mathsf{x}}}^{\dagger}}\right>=q\Psi({\boldsymbol{\mathsf{x}}}_{k})$,
and

$\|g({\boldsymbol{\mathsf{x}}},{\boldsymbol{\mathsf{y}}},i)\|_{{\cal
X}^{\ast}}\leq
L_{\max}\|{\mathbf{A}}_{i}{\boldsymbol{\mathsf{x}}}-{\boldsymbol{\mathsf{y}}}_{i}\|_{{\cal
Y}}^{q-1}.$

Then taking the conditional expectation of
$\|g({\boldsymbol{\mathsf{x}}},{\boldsymbol{\mathsf{y}}},i)\|_{{\cal
X}^{\ast}}^{p^{\ast}}$ yields

$\displaystyle\mathbb{E}\big{[}\|g({\boldsymbol{\mathsf{x}}},{\boldsymbol{\mathsf{y}}},i)\|_{{\cal
X}^{\ast}}^{p^{\ast}}\big{]}$ $\displaystyle\leq
L_{\max}^{p^{\ast}}\mathbb{E}\Big{[}\|{\mathbf{A}}_{i}{\boldsymbol{\mathsf{x}}}-{\boldsymbol{\mathsf{y}}}_{i}\|_{{\cal
Y}}^{p^{\ast}(q-1)}\Big{]}=L_{\max}^{p^{\ast}}\mathbb{E}\Big{[}(\|{\mathbf{A}}_{i}{\boldsymbol{\mathsf{x}}}-{\boldsymbol{\mathsf{y}}}_{i}\|_{{\cal
Y}}^{q})^{\frac{p^{\ast}}{q^{\ast}}}\Big{]}.$

We have $0<\frac{p^{\ast}}{q^{\ast}}\leq 1$, with the equality achieved only
if $p^{\ast}=q^{\ast}=2$. In the latter case, it trivially follows that
$\mathbb{E}[\|g({\boldsymbol{\mathsf{x}}},{\boldsymbol{\mathsf{y}}},i)\|_{{\cal
X}^{\ast}}^{p^{\ast}}]\leq
qL_{\max}^{p^{\ast}}\Psi({\boldsymbol{\mathsf{x}}})$. If
$0<\frac{p^{\ast}}{q^{\ast}}<1$, by Jensen’s inequality, we have

$\displaystyle\mathbb{E}\big{[}\|g({\boldsymbol{\mathsf{x}}},{\boldsymbol{\mathsf{y}}},i)\|_{{\cal
X}^{\ast}}^{p^{\ast}}\big{]}$ $\displaystyle\leq
L_{\max}^{p^{\ast}}\mathbb{E}\Big{[}(\|{\mathbf{A}}_{i}{\boldsymbol{\mathsf{x}}}-{\boldsymbol{\mathsf{y}}}_{i}\|_{{\cal
Y}}^{q})^{\frac{p^{\ast}}{q^{\ast}}}\Big{]}\leq
L_{\max}^{p^{\ast}}(\mathbb{E}[\|{\mathbf{A}}_{i}{\boldsymbol{\mathsf{x}}}-{\boldsymbol{\mathsf{y}}}_{i}\|_{{\cal
Y}}^{q}])^{\frac{p^{\ast}}{q^{\ast}}}\leq
L_{\max}^{p^{\ast}}q^{\frac{p^{\ast}}{q^{\ast}}}\Psi({\boldsymbol{\mathsf{x}}})^{\frac{p^{\ast}}{q^{\ast}}}.$

Plugging this estimate into the conditional descent property (18) yields

$\displaystyle\mathbb{E}_{k}[\Delta_{k+1}]\leq\Delta_{k}-q\mu_{k+1}\Psi({\boldsymbol{\mathsf{x}}}_{k})+L_{\max}^{p^{\ast}}q^{\frac{p^{\ast}}{q^{\ast}}}\frac{G_{p^{\ast}}}{p^{\ast}}\mu_{k+1}^{p^{\ast}}\Psi({\boldsymbol{\mathsf{x}}}_{k})^{\frac{p^{\ast}}{q^{\ast}}}.$

Since the sequence $({\boldsymbol{\mathsf{x}}}_{k})_{k\in\mathbb{N}}$ is
uniformly bounded by Lemma 3.12, so is
$(\Psi({\boldsymbol{\mathsf{x}}}_{k}))_{k\in\mathbb{N}}$, and we thus have

$\sum_{k=0}^{\infty}\mu_{k+1}^{p^{\ast}}\Psi({\boldsymbol{\mathsf{x}}}_{k})^{\frac{p^{\ast}}{q^{\ast}}}\leq
C\sum_{k=0}^{\infty}\mu_{k+1}^{p^{\ast}}<\infty.$

Thus, we can apply Robbins-Siegmund theorem for almost super-martingales, and
deduce that $(\Delta_{k})_{k\in\mathbb{N}}$ converges a.s. to a non-negative
random variable $\Delta_{\infty}$. Moreover,
$\sum_{k=0}^{\infty}\mu_{k+1}\Psi({\boldsymbol{\mathsf{x}}}_{k})<\infty$ holds
a.s. By repeating the argument for Theorem 3.8, there exists a subsequence
$({\boldsymbol{\mathsf{x}}}_{k_{j}})_{j\in\mathbb{N}}$ that a.s. converges to
some $\widehat{\boldsymbol{\mathsf{x}}}\in{\cal X}_{\min}$, and hence
$\Delta_{\infty}=0$, as desired. Moreover, by Lemma 3.12, the sequence
$(\Delta_{k})_{k\in\mathbb{N}}$ is bounded, and thus uniformly integrable.
Since it converges to $0$ a.s., from Vitali’s theorem it follows that
$\lim_{k\rightarrow\infty}\mathbb{E}[{\mathbf{B}}_{p}({\boldsymbol{\mathsf{x}}}_{k},{\boldsymbol{\mathsf{x}}}^{\dagger})]=0$.
∎

The results in Theorem 3.13 are similar to that of Theorem 3.10, but the
generality of the latter is compensated for by an additional step-size
assumption ensuring boundedness of iterates
$({\boldsymbol{\mathsf{x}}}_{k})_{k\in\mathbb{N}}$.

###  3.3 Convergence rates for conditionally stable operators

Theorem 3.10 states the conditions needed for the convergence of Bregman
distances in expectation. However, it does not provide convergence rates. In
order to obtain convergence rates, one needs additional conditions on the MNS
${\boldsymbol{\mathsf{x}}}^{\dagger}$, which are collectively known as source
conditions. One approach is via conditional stability: for a locally
conditionally stable operator, we can extract convergence in expectation and
quantify the convergence speed. Conditional stability is known for many
inverse problems for PDEs, and has been used extensively to investigate
regularised solutions [8, 13]. It is useful for analysing ill-posed problems
that are locally well-posed, and in case of a (possibly) non-linear forward
operator $F$ it is of the form

$\displaystyle\|{\boldsymbol{\mathsf{x}}}_{1}-{\boldsymbol{\mathsf{x}}}_{2}\|_{{\cal
X}}\leq\Phi(\|F({\boldsymbol{\mathsf{x}}}_{1})-F({\boldsymbol{\mathsf{x}}}_{2})\|_{{\cal
Y}}),\quad\forall{\boldsymbol{\mathsf{x}}}_{1},{\boldsymbol{\mathsf{x}}}_{2}\in{\cal
M}\subset{\cal X},$ (19)

where $\Phi:[0,\infty)\rightarrow[0,\infty)$ with $\Phi(0)=0$ is a continuous,
non-decreasing function, and ${\cal M}$ is typically a ball in the ambient
norm [19]. In Banach space settings, the conditional stability needs to be
adjusted, by replacing the left hand side of (19) with a non-negative error
measure [7]. Since the most relevant error measure for Banach space analysis
is the Bregman distance
${\mathbf{B}}_{p}({\boldsymbol{\mathsf{x}}}_{1},{\boldsymbol{\mathsf{x}}}_{2})$,
a Hölder type stability estimate then reads: for some $\alpha\geq 1$ and
$C_{\alpha}>0$

${\mathbf{B}}_{p}({\boldsymbol{\mathsf{x}}},{\boldsymbol{\mathsf{x}}}^{\dagger})^{\alpha}\leq
C_{\alpha}^{-1}\|{\mathbf{A}}{\boldsymbol{\mathsf{x}}}-{\mathbf{A}}{\boldsymbol{\mathsf{x}}}^{\dagger}\|_{{\cal
Y}}^{p}.$ (20)

Now we give a convergence rate under conditional stability bound (20). The
constant $C_{N}$ appears in Lemma 3.5 and denotes the norm equivalence
constant.

######  Theorem 3.14.

Let the forward operator ${\mathbf{A}}$ satisfy the conditional stability
bound (20) for some $\alpha\geq 1$ and $C_{\alpha}>0$. Let ${\cal
J}_{p}^{{\cal
X}}({\boldsymbol{\mathsf{x}}}_{0})\in\overline{\operatorname{range}({\mathbf{A}}^{\ast})}$,
and for
$C_{k}=C_{N}C_{\alpha}(1-L_{\max}^{p^{\ast}}\frac{G_{p^{\ast}}}{p^{\ast}}\mu_{k}^{p^{\ast}-1})>0$,
the step-sizes satisfy $\sum_{k=1}^{\infty}\mu_{k}C_{k}=\infty$. Then there
holds

$\lim_{k\rightarrow\infty}\mathbb{E}[{\mathbf{B}}_{p}({\boldsymbol{\mathsf{x}}}_{k},{\boldsymbol{\mathsf{x}}}^{\dagger})]=0.$

Moreover,

$\mathbb{E}[{\mathbf{B}}_{p}({\boldsymbol{\mathsf{x}}}_{k},{\boldsymbol{\mathsf{x}}}^{\dagger})]\leq\left\\{\begin{aligned}
\frac{{\mathbf{B}}_{p}({\boldsymbol{\mathsf{x}}}_{0},{\boldsymbol{\mathsf{x}}}^{\dagger})}{\Big{(}1+(\alpha-1){\mathbf{B}}_{p}({\boldsymbol{\mathsf{x}}}_{0},{\boldsymbol{\mathsf{x}}}^{\dagger})^{\alpha-1}\sum_{j=1}^{k}\mu_{j}C_{j}\Big{)}^{\frac{1}{\alpha-1}}},&\quad\text{
if }\alpha>1,\\\
{\exp\Big{(}-\sum_{j=1}^{k}\mu_{j}C_{j}\Big{)}{\mathbf{B}}_{p}({\boldsymbol{\mathsf{x}}}_{0},{\boldsymbol{\mathsf{x}}}^{\dagger}),}&\quad\text{
if }\alpha=1.\end{aligned}\right.$

###### Proof.

Let
$\Delta_{k}:={\mathbf{B}}_{p}({\boldsymbol{\mathsf{x}}}_{k},{\boldsymbol{\mathsf{x}}}^{\dagger})$.
The proof of Theorem 3.8 and the conditional stability bound (20) imply

$\displaystyle\mathbb{E}_{k}[\Delta_{k+1}]$
$\displaystyle\leq\Delta_{k}-p\mu_{k+1}\left(1-L_{\max}^{p^{\ast}}\frac{G_{p^{\ast}}}{p^{\ast}}\mu_{k+1}^{p^{\ast}-1}\right)\Psi({\boldsymbol{\mathsf{x}}}_{k})$
(21)
$\displaystyle\leq\Delta_{k}-p\mu_{k+1}\frac{C_{N}C_{\alpha}}{p}\Big{(}1-L_{\max}^{p^{\ast}}\frac{G_{p^{\ast}}}{p^{\ast}}\mu_{k+1}^{p^{\ast}-1}\Big{)}\Delta_{k}^{\alpha},$

since by Lemma 3.5, there exists a $C_{N}>0$ such that
$\Psi({\boldsymbol{\mathsf{x}}})\geq\frac{C_{N}}{p}\|{\mathbf{A}}{\boldsymbol{\mathsf{x}}}-{\boldsymbol{\mathsf{y}}}\|_{{\cal
Y}}^{p}$. Taking the full expectation and using Jensen’s inequality lead to

$\mathbb{E}[\Delta_{k+1}]\leq\mathbb{E}[\Delta_{k}]-\mu_{k+1}C_{k+1}\mathbb{E}[\Delta_{k}]^{\alpha}.$

Since $C_{k+1}>0$ by assumption, $(\mathbb{E}[\Delta_{k}])_{k\in\mathbb{N}}$
is a monotonically decreasing sequence. By the convexity of the function
$x\mapsto x^{\alpha}$ (for $\alpha\geq 1$), for any $\epsilon>0$ and
$x\geq\epsilon$, we have $\epsilon^{\alpha}\geq\frac{\epsilon}{x}x^{\alpha}$.
We claim that for every $\epsilon>0$, there exists a
$k_{\epsilon}\in\mathbb{N}$ such that $\mathbb{E}[\Delta_{k}]\leq\epsilon$ for
all $k\geq k_{\epsilon}$. Assuming the contrary,
$\mathbb{E}[\Delta_{k}]\geq\epsilon$ for all $k$, gives

$\displaystyle\mathbb{E}[\Delta_{k+1}]\leq\mathbb{E}[\Delta_{k}]-\mu_{k+1}C_{k+1}\mathbb{E}[\Delta_{k}]^{\alpha}\leq\mathbb{E}[\Delta_{k}]-\mu_{k+1}C_{k+1}\epsilon^{\alpha}\leq\Delta_{0}-\epsilon^{\alpha}\sum_{j=1}^{k+1}\mu_{j}C_{j}\rightarrow-\infty,$

since $\sum_{j=1}^{\infty}\mu_{j}C_{j}=\infty$ by assumption, which is a
contradiction. Therefore, $\lim_{k\rightarrow\infty}\mathbb{E}[\Delta_{k}]=0$.
For $\alpha>1$, by Polyak’s inequality (cf. Lemma A.1), we have

$\displaystyle\mathbb{E}[\Delta_{k+1}]\leq\frac{\Delta_{0}}{\Big{(}1+(\alpha-1)\Delta_{0}^{\alpha-1}\sum_{j=1}^{k+1}\mu_{j}C_{j}\Big{)}^{\frac{1}{\alpha-1}}}.$

Meanwhile, for $\alpha=1$, using the inequality $1-x\leq e^{-x}$ for $x\geq
0$, a direct computation yields

$\displaystyle\mathbb{E}[\Delta_{k+1}]\leq(1-\mu_{k+1}C_{k+1})\mathbb{E}[\Delta_{k}]\leq\prod_{j=1}^{k+1}(1-\mu_{j}C_{j})\Delta_{0}\leq{\exp\Big{(}-\sum_{j=1}^{k+1}\mu_{j}C_{j}\Big{)}\Delta_{0},}$

completing the proof of the theorem. ∎

######  Remark 3.15.

We have the following comments on Theorem 3.14.

  * (i)

The estimates for $\alpha>1$ and $\alpha=1$ in Theorem 3.14 are consistent in
the sense that

$\lim_{\alpha\searrow
1}\frac{{\mathbf{B}}_{p}({\boldsymbol{\mathsf{x}}}_{0},{\boldsymbol{\mathsf{x}}}^{\dagger})}{\Big{(}1+(\alpha-1){\mathbf{B}}_{p}({\boldsymbol{\mathsf{x}}}_{0},{\boldsymbol{\mathsf{x}}}^{\dagger})^{\alpha-1}\sum_{j=1}^{k}\mu_{j}C_{j}\Big{)}^{\frac{1}{\alpha-1}}}=\exp\Big{(}-\sum_{j=1}^{k}\mu_{j}C_{j}\Big{)}{\mathbf{B}}_{p}({\boldsymbol{\mathsf{x}}}_{0},{\boldsymbol{\mathsf{x}}}^{\dagger}).$

  * (ii)

While it might seem counter-intuitive, $\alpha=1$ gives a better convergence
rate than $\alpha>1$, because of the following

$\displaystyle{\mathbf{B}}_{p}({\boldsymbol{\mathsf{x}}},{\boldsymbol{\mathsf{x}}}^{\dagger})^{\alpha}\geq{\mathbf{B}}_{p}({\boldsymbol{\mathsf{x}}},{\boldsymbol{\mathsf{x}}}^{\dagger})^{\tilde{\alpha}}\text{
if and only if
}\alpha\log{\mathbf{B}}_{p}({\boldsymbol{\mathsf{x}}},{\boldsymbol{\mathsf{x}}}^{\dagger})\geq\tilde{\alpha}\log{\mathbf{B}}_{p}({\boldsymbol{\mathsf{x}}},{\boldsymbol{\mathsf{x}}}^{\dagger}).$

Hence, whenever
${\mathbf{B}}_{p}({\boldsymbol{\mathsf{x}}},{\boldsymbol{\mathsf{x}}}^{\dagger})<1$,
we have
${\mathbf{B}}_{p}({\boldsymbol{\mathsf{x}}},{\boldsymbol{\mathsf{x}}}^{\dagger})\geq{\mathbf{B}}_{p}({\boldsymbol{\mathsf{x}}},{\boldsymbol{\mathsf{x}}}^{\dagger})^{\alpha}$
for $\alpha>1$. Plugging this into the conditional stability bound (20) yields

$\displaystyle{\mathbf{B}}_{p}({\boldsymbol{\mathsf{x}}},{\boldsymbol{\mathsf{x}}}^{\dagger})^{\alpha}\leq{\mathbf{B}}_{p}({\boldsymbol{\mathsf{x}}},{\boldsymbol{\mathsf{x}}}^{\dagger})\leq
C_{1}^{-1}\|{\mathbf{A}}{\boldsymbol{\mathsf{x}}}-{\mathbf{A}}{\boldsymbol{\mathsf{x}}}^{\dagger}\|_{{\cal
Y}}^{p}=C_{1}^{-1}pN\Psi({\boldsymbol{\mathsf{x}}}).$

Meanwhile, the proof of Theorem 3.14 uses the conditional stability bound to
establish a relationship between the objective value and the Bregman distance
to the MNS ${\boldsymbol{\mathsf{x}}}^{\dagger}$, cf. (21). Putting these
together gives that $\alpha=1$ provides a greater decrease of the expected
Bregman distance, once we are close enough to the solution.

The conditional stability estimate (20) for a linear operator ${\mathbf{A}}$
implies its injectivity. Then the objective $\Psi({\boldsymbol{\mathsf{x}}})$
is strongly convex. Under condition (20), there can indeed be only one
solution: if
${\mathbf{A}}\tilde{\boldsymbol{\mathsf{x}}}={\mathbf{A}}{\boldsymbol{\mathsf{x}}}$,
then
${\mathbf{B}}_{p}(\tilde{\boldsymbol{\mathsf{x}}},{\boldsymbol{\mathsf{x}}})=0$
follows from (20). The step-size condition
$\sum_{k=1}^{\infty}\mu_{k}C_{k}=\infty$ is weaker than that in Theorem 3.10.
Namely, it follows from step-size conditions in Theorem 3.8, since

$\sum_{k=1}^{\infty}\mu_{k}C_{k}=C_{N}C_{\alpha}\Big{(}\sum_{k=1}^{\infty}\mu_{k}-L_{\max}^{p^{\ast}}\frac{G_{p^{\ast}}}{p^{\ast}}\sum_{k=1}^{\infty}\mu_{k}^{p^{\ast}}\Big{)}=\infty$

holds if $\sum_{k=1}^{\infty}\mu_{k}=\infty$ and
$\sum_{k=1}^{\infty}\mu_{k}^{p^{\ast}}<\infty$. Further, if there exists a
$C>0$ such that
$1-L_{\max}^{p^{\ast}}\frac{G_{p^{\ast}}}{p^{\ast}}\mu_{k}^{p^{\ast}-1}>C$
holds for all $k\in\mathbb{N}$, e.g. if $\mu_{k}$ is a constant satisfying
this condition, then $\sum_{k=1}^{\infty}\mu_{k}C_{k}=\infty$ is weaker than
the conditions in Theorem 3.8, since the condition
$\sum_{k=1}^{\infty}\mu_{k}^{p^{\ast}}<\infty$ is no longer needed for
convergence, and $\sum_{k=1}^{\infty}\mu_{k}=\infty$ suffices. Moreover, we
can choose constant step-sizes. Indeed, setting $\mu_{k}=\mu_{0}$, with
$1-L_{\max}^{p^{\ast}}\frac{G_{p^{\ast}}}{p^{\ast}}\mu_{0}^{p^{\ast}-1}=\frac{1}{2}$,
we get an exponential convergence rate for $\alpha=1$, since
$C_{k}=\frac{C_{N}C_{\alpha}}{2}$, we have

$\displaystyle\mathbb{E}[\Delta_{k+1}]$
$\displaystyle\leq(1-\mu_{0}C_{k+1})\mathbb{E}[\Delta_{k}]\leq\bigg{(}1-2^{-1-1/(p^{\ast}-1)}L_{\max}^{-p^{\ast}/p^{\ast}-1}\Big{(}\frac{p^{\ast}}{G_{p^{\ast}}}C_{N}C_{\alpha}\Big{)}^{1/p^{\ast}-1}\bigg{)}^{k}\mathbb{E}[\Delta_{0}]$
$\displaystyle\leq\bigg{(}1-2^{-p}L_{\max}^{-p}\Big{(}\frac{p^{\ast}}{G_{p^{\ast}}}\Big{)}^{p^{\ast}/p}C_{N}C_{\alpha}\bigg{)}^{k}\Delta_{0}.$

Note that this convergence rate is largely comparable with that in the Hilbert
case: the conditional stability bound implies the strict convexity of the
quadratic objective $\Psi({\boldsymbol{\mathsf{x}}})$, and the SGD is known to
converge exponentially fast (see e.g. [16, Theorem 3.1]), with the rate
determined by a variant of the condition number.

######  Remark 3.16.

The conditional stability bound (20) is stated globally. However, such
conditions are often valid only locally. A local definition could have been
employed in (20), with minor modifications of the argument. Indeed, by the
argument of Theorem 3.10, we appeal to Lemma A.2, showing that the Bregman
distances of the iterates are non-increasing. Thus, it suffices to assume that
the initial point ${\boldsymbol{\mathsf{x}}}_{0}$ is sufficiently close to the
MNS ${\boldsymbol{\mathsf{x}}}^{\dagger}$.

######  Remark 3.17.

Conditional stability is intimately tied with classical source conditions. For
example, as shown in [41], assuming $\alpha=1$ in (20) allows to show a
variational inequality

$\left<{{\cal J}_{p}^{{\cal
X}}({\boldsymbol{\mathsf{x}}}^{\dagger})},{{\boldsymbol{\mathsf{x}}}-{\boldsymbol{\mathsf{x}}}^{\dagger}}\right>\leq\|{\boldsymbol{\mathsf{x}}}^{\dagger}\|_{{\cal
X}}^{p-1}C_{\alpha}^{-1}(pC_{p}^{-1})^{1/p}\|{\mathbf{A}}({\boldsymbol{\mathsf{x}}}-{\boldsymbol{\mathsf{x}}}^{\dagger})\|_{{\cal
Y}}.$

Then Hahn-Banach theorem and [41, Lemma 8.21] give the canonical range type
condition ${\cal J}_{p}^{{\cal
X}}({\boldsymbol{\mathsf{x}}}^{\dagger})={\mathbf{A}}^{\ast}{\boldsymbol{\mathsf{w}}}$,
for ${\boldsymbol{\mathsf{w}}}\in{\cal X}$ such that
$\|{\boldsymbol{\mathsf{w}}}\|_{{\cal X}}\leq 1$. Connections between source
conditions and conditional stability estimates have been studied, e.g. for
linear operators in Hilbert spaces [47] and in ${\cal L}^{p}$ spaces [5].
Moreover, variational source conditions often imply conditional stability
estimates [20], and in case of bijective and continuous operators they are
trivially inferred by a standard source condition (albeit only in a possibly
small neighbourhood around the solution). See the book [50] about the
connections between source conditions and conditional stability estimates, and
[21] for inverse problems for differential equations.

##  4 Regularising property

In practice, we often do not have access to the exact data
${\boldsymbol{\mathsf{y}}}$ but only to noisy observations
${\boldsymbol{\mathsf{y}}}^{\delta}$, such that
$\|{\boldsymbol{\mathsf{y}}}^{\delta}-{\boldsymbol{\mathsf{y}}}\|_{{\cal
Y}}{\leq}\delta$. The convergence study in the presence of observational noise
requires a different approach, since the sequence of objective values
$(\|{\mathbf{A}}{\boldsymbol{\mathsf{x}}}_{k}^{\delta}-{\boldsymbol{\mathsf{y}}}^{\delta}\|_{{\cal
Y}}^{p})_{k\in\mathbb{N}}$ generally will not converge to $0$. In this section
we show that SGD has a regularising effect, in the sense that the expected
error
$\mathbb{E}[{\mathbf{B}}_{p}({\boldsymbol{\mathsf{x}}}_{k(\delta)}^{\delta},{\boldsymbol{\mathsf{x}}}^{\dagger})]$
converges to $0$ as the noise level $\delta$ decays to $0$, for properly
selected stopping indices $k(\delta)$.

Let $({\boldsymbol{\mathsf{x}}}_{k})_{k\in\mathbb{N}}$ and
$({\boldsymbol{\mathsf{x}}}_{k}^{\delta})_{k\in\mathbb{N}}$ be the noiseless
and noisy iterates, defined respectively by

$\displaystyle{\boldsymbol{\mathsf{x}}}_{k+1}$ $\displaystyle={\cal
J}_{p^{\ast}}^{{\cal X}^{\ast}}\left({\cal J}_{p}^{{\cal
X}}({\boldsymbol{\mathsf{x}}}_{k})-\mu_{k+1}{\boldsymbol{\mathsf{g}}}_{k+1}\right),\quad\mbox{with
}{\boldsymbol{\mathsf{g}}}_{k+1}=g({\boldsymbol{\mathsf{x}}}_{k},{\boldsymbol{\mathsf{y}}},i_{k+1}),$
(22) $\displaystyle{\boldsymbol{\mathsf{x}}}_{k+1}^{\delta}$
$\displaystyle={\cal J}_{p^{\ast}}^{{\cal X}^{\ast}}\left({\cal J}_{p}^{{\cal
X}}({\boldsymbol{\mathsf{x}}}_{k}^{\delta})-\mu_{k+1}{\boldsymbol{\mathsf{g}}}^{\delta}_{k+1}\right),\quad\mbox{with
}{\boldsymbol{\mathsf{g}}}_{k+1}^{\delta}=g({\boldsymbol{\mathsf{x}}}_{k}^{\delta},{\boldsymbol{\mathsf{y}}}^{\delta},i_{k+1}).$
(23)

The key step in proving the regularising property is to show the stability of
SGD iterates with respect to noise. The noise enters into the iterations
through the update directions ${\boldsymbol{\mathsf{g}}}^{\delta}_{k+1}$ and
thus, the stability of the iterates requires that of update directions. This
however requires imposing suitable assumptions on the observation space ${\cal
Y}$ since in general, the single valued duality maps $\jmath_{p}^{{\cal Y}}$
are continuous only at $0$. If ${\cal Y}$ is uniformly smooth, the
corresponding duality maps are also smooth. This assumption is also needed for
deterministic iterates, cf. [46, Proposition 6.17] or [35, Lemma 9]. Thus we
make the following assumption.

######  Assumption 4.1.

The Banach space ${\cal X}$ is $p$-convex and uniformly smooth, and ${\cal Y}$
is uniformly smooth.

We then have the following stability result on the iterates with respect to
noise, whose elementary but lengthy proof is deferred to the appendix.

######  Lemma 4.2.

Let Assumption 4.1 hold. Consider the iterations (22) and (23) with the same
initialisation
${\boldsymbol{\mathsf{x}}}_{0}^{\delta}={\boldsymbol{\mathsf{x}}}_{0}$, and
following the same path (i.e. using same random indices $i_{k}$). Then, for
any fixed $k\in\mathbb{N}$, we have

$\displaystyle\lim_{\delta\searrow
0}\mathbb{E}[{\mathbf{B}}_{p}({\boldsymbol{\mathsf{x}}}^{\delta}_{k},{\boldsymbol{\mathsf{x}}}_{k})]=\lim_{\delta\searrow
0}\mathbb{E}[\|{\boldsymbol{\mathsf{x}}}^{\delta}_{k}-{\boldsymbol{\mathsf{x}}}_{k}\|_{{\cal
X}}]=\lim_{\delta\searrow 0}\mathbb{E}[\|{\cal J}_{p}^{{\cal
X}}({\boldsymbol{\mathsf{x}}}_{k}^{\delta})-{\cal J}_{p}^{{\cal
X}}({\boldsymbol{\mathsf{x}}}_{k})\|_{{\cal X}^{\ast}}]=0.$

Now we show the regularising property of SGD for suitable stopping indices
$k(\delta)$.

######  Theorem 4.3.

Let Assumption 4.1 hold, and the step-sizes $(\mu_{k})_{k\in\mathbb{N}}$
satisfy $\sum_{k=1}^{\infty}\mu_{k}=\infty$,
$\sum_{k=1}^{\infty}\mu_{k}^{p^{\ast}}<\infty$ and
$1-L_{\max}^{p^{\ast}}\frac{G_{p^{\ast}}}{p^{\ast}}\mu_{k}^{p^{\ast}-1}>C>0$.
If $\lim_{\delta\searrow 0}k(\delta)=\infty$ and $\lim_{\delta\searrow
0}\delta^{p}\sum_{\ell=1}^{k(\delta)}\mu_{\ell}=0$, then

$\displaystyle\lim_{\delta\searrow
0}\mathbb{E}[{\mathbf{B}}_{p}({\boldsymbol{\mathsf{x}}}_{k(\delta)}^{\delta},{\boldsymbol{\mathsf{x}}}^{\dagger})]=0.$

###### Proof.

Let
$\Delta_{k}={\mathbf{B}}_{p}({\boldsymbol{\mathsf{x}}}_{k},{\boldsymbol{\mathsf{x}}}^{\dagger})$
and
$\Delta_{k}^{\delta}={\mathbf{B}}_{p}({\boldsymbol{\mathsf{x}}}_{k}^{\delta},{\boldsymbol{\mathsf{x}}}^{\dagger})$.
Take any $\delta>0$ and $k\in\mathbb{N}$. By the three point identity (4), we
have

$\displaystyle\Delta_{k}^{\delta}$
$\displaystyle={\mathbf{B}}_{p}({\boldsymbol{\mathsf{x}}}_{k}^{\delta},{\boldsymbol{\mathsf{x}}}_{k})+\Delta_{k}+\left<{{\cal
J}_{p}^{{\cal X}}({\boldsymbol{\mathsf{x}}}_{k})-{\cal J}_{p}^{{\cal
X}}({\boldsymbol{\mathsf{x}}}_{k}^{\delta})},{{\boldsymbol{\mathsf{x}}}_{k}-{\boldsymbol{\mathsf{x}}}^{\dagger}}\right>$
$\displaystyle\leq{\mathbf{B}}_{p}({\boldsymbol{\mathsf{x}}}_{k}^{\delta},{\boldsymbol{\mathsf{x}}}_{k})+\Delta_{k}+\|{\cal
J}_{p}^{{\cal X}}({\boldsymbol{\mathsf{x}}}_{k})-{\cal J}_{p}^{{\cal
X}}({\boldsymbol{\mathsf{x}}}_{k}^{\delta})\|_{{\cal
X}^{\ast}}\|{\boldsymbol{\mathsf{x}}}_{k}-{\boldsymbol{\mathsf{x}}}^{\dagger}\|_{{\cal
X}}.$ (24)

Consider a sequence $(\delta_{j})_{j\in\mathbb{N}}$ decaying to zero. Taking
any $\epsilon>0$, it suffices to find a $j_{\epsilon}\in\mathbb{N}$ such that
for all $j\geq j_{\epsilon}$ we have
$\mathbb{E}[\Delta_{k(\delta_{j})}^{\delta_{j}}]\leq 4\epsilon$. By Theorem
3.10, there exists a $k_{\epsilon}\in\mathbb{N}$ such that for all $k\geq
k_{\epsilon}$ we have

$\displaystyle\mathbb{E}[\Delta_{k}]<\epsilon\quad\text{and}\quad\mathbb{E}[\|{\boldsymbol{\mathsf{x}}}_{k}-{\boldsymbol{\mathsf{x}}}^{\dagger}\|_{{\cal
X}}]<\epsilon^{1/2}.$ (25)

Moreover, for any fixed $k_{\epsilon}$, by Lemma 4.2, there exists
$j_{1}\in\mathbb{N}$ such that for all $j\geq j_{1}$ we have

$\displaystyle\mathbb{E}[{\mathbf{B}}_{p}({\boldsymbol{\mathsf{x}}}_{k_{\epsilon}}^{\delta_{j}},{\boldsymbol{\mathsf{x}}}_{k_{\epsilon}})]<\epsilon\quad\text{and}\quad\mathbb{E}[\|{\cal
J}_{p}^{{\cal X}}({\boldsymbol{\mathsf{x}}}_{k_{\epsilon}})-{\cal
J}_{p}^{{\cal
X}}({\boldsymbol{\mathsf{x}}}_{k_{\epsilon}}^{\delta_{j}})\|_{{\cal
X}^{\ast}}]<\epsilon^{1/2}.$ (26)

Thus, plugging the estimates (25) and (26) into (4), we have
$\mathbb{E}[\Delta_{k_{\epsilon}}^{\delta_{j}}]<3\epsilon$, for all $j\geq
j_{1}$. Note, however, that the same does not necessarily hold for all $k\geq
k_{\epsilon}$, and thus for a monotonically increasing sequence of stopping
indices $k(\delta_{j})$, since
$\mathbb{E}[\Delta_{k(\delta_{j})}^{\delta_{j}}]$ are not necessarily
monotone. Instead, taking the expectation of the descent property (13) with
respect to ${\cal F}_{k}$ yields

$\displaystyle\mathbb{E}_{k}[\Delta_{k+1}^{\delta}]\leq\Delta_{k}^{\delta}-\mu_{k+1}\left<{\mathbb{E}_{k}[{\boldsymbol{\mathsf{g}}}_{k+1}^{\delta}]},{{\boldsymbol{\mathsf{x}}}_{k}^{\delta}-{\boldsymbol{\mathsf{x}}}^{\dagger}}\right>+pL_{\max}^{p^{\ast}}\frac{G_{p^{\ast}}}{p^{\ast}}\mu_{k+1}^{p^{\ast}}\Psi({\boldsymbol{\mathsf{x}}}_{k}^{\delta}).$

Then we decompose the middle term into

$\displaystyle\left<{\mathbb{E}_{k}[{\boldsymbol{\mathsf{g}}}_{k+1}^{\delta}]},{{\boldsymbol{\mathsf{x}}}^{\dagger}-{\boldsymbol{\mathsf{x}}}_{k}^{\delta}}\right>$
$\displaystyle=\frac{1}{N}\sum_{i=1}^{N}\left<{\jmath_{p}^{{\cal
Y}}({\mathbf{A}}_{i}{\boldsymbol{\mathsf{x}}}_{k}^{\delta}-{\boldsymbol{\mathsf{y}}}_{i}^{\delta})},{-({\mathbf{A}}_{i}{\boldsymbol{\mathsf{x}}}_{k}^{\delta}-{\boldsymbol{\mathsf{y}}}_{i}^{\delta})+{\boldsymbol{\mathsf{y}}}_{i}-{\boldsymbol{\mathsf{y}}}_{i}^{\delta}}\right>$
$\displaystyle=-p\Psi({\boldsymbol{\mathsf{x}}}_{k}^{\delta})+\frac{1}{N}\sum_{i=1}^{N}\left<{\jmath_{p}^{{\cal
Y}}({\mathbf{A}}_{i}{\boldsymbol{\mathsf{x}}}_{k}^{\delta}-{\boldsymbol{\mathsf{y}}}_{i}^{\delta})},{{\boldsymbol{\mathsf{y}}}_{i}-{\boldsymbol{\mathsf{y}}}_{i}^{\delta}}\right>$
$\displaystyle\leq-p\Psi({\boldsymbol{\mathsf{x}}}_{k}^{\delta})+\frac{1}{N}\sum_{i=1}^{N}\|{\mathbf{A}}_{i}{\boldsymbol{\mathsf{x}}}_{k}^{\delta}-{\boldsymbol{\mathsf{y}}}_{i}^{\delta}\|_{{\cal
Y}}^{p-1}\|{\boldsymbol{\mathsf{y}}}_{i}-{\boldsymbol{\mathsf{y}}}_{i}^{\delta}\|_{{\cal
Y}}$
$\displaystyle\leq-p\Psi({\boldsymbol{\mathsf{x}}}_{k}^{\delta})+\delta\frac{1}{N}\sum_{i=1}^{N}\|{\mathbf{A}}_{i}{\boldsymbol{\mathsf{x}}}_{k}^{\delta}-{\boldsymbol{\mathsf{y}}}_{i}^{\delta}\|_{{\cal
Y}}^{p-1},$

where we have used (3) and the Cauchy-Schwarz inequality. Taking the full
expectation gives

$\displaystyle\mathbb{E}[\Delta^{\delta}_{k+1}]$
$\displaystyle\leq\mathbb{E}[\Delta^{\delta}_{k}]\\!-\\!p\mu_{k+1}\mathbb{E}[\Psi({\boldsymbol{\mathsf{x}}}_{k}^{\delta})]\\!+\\!pL_{\max}^{p^{\ast}}\frac{G_{p^{\ast}}}{p^{\ast}}\mu_{k+1}^{p^{\ast}}\mathbb{E}[\Psi({\boldsymbol{\mathsf{x}}}_{k}^{\delta})]\\!+\\!\delta\mu_{k+1}\frac{1}{N}\sum_{i=1}^{N}\mathbb{E}[\|{\mathbf{A}}_{i}{\boldsymbol{\mathsf{x}}}_{k}^{\delta}-{\boldsymbol{\mathsf{y}}}_{i}^{\delta}\|_{{\cal
Y}}^{p-1}]$
$\displaystyle=\mathbb{E}[\Delta^{\delta}_{k}]-p\mu_{k+1}C_{k+1}\mathbb{E}[\Psi({\boldsymbol{\mathsf{x}}}_{k}^{\delta})]+\delta\mu_{k+1}\frac{1}{N}\sum_{i=1}^{N}\mathbb{E}[\|{\mathbf{A}}_{i}{\boldsymbol{\mathsf{x}}}_{k}^{\delta}-{\boldsymbol{\mathsf{y}}}_{i}^{\delta}\|_{{\cal
Y}}^{p-1}],$

where
$C_{k}\\!=\\!1-L_{\max}^{p^{\ast}}\frac{G_{p^{\ast}}}{p^{\ast}}\mu_{k}^{p^{\ast}-1}>C>0$.
Now using the Lyapunov inequality

$\displaystyle\frac{1}{N}\sum_{i=1}^{N}\mathbb{E}[\|{\mathbf{A}}_{i}{\boldsymbol{\mathsf{x}}}_{k}^{\delta}-{\boldsymbol{\mathsf{y}}}_{i}^{\delta}\|_{{\cal
Y}}^{p-1}]\leq\frac{1}{N}\sum_{i=1}^{N}\Big{(}\mathbb{E}[\|{\mathbf{A}}_{i}{\boldsymbol{\mathsf{x}}}_{k}^{\delta}-{\boldsymbol{\mathsf{y}}}_{i}^{\delta}\|_{{\cal
Y}}^{p}]\Big{)}^{(p-1)/p}=p^{1/p^{\ast}}\frac{1}{N}\sum_{i=1}^{N}\Big{(}\mathbb{E}[\Psi_{i}({\boldsymbol{\mathsf{x}}}_{k}^{\delta})]\Big{)}^{1/p^{\ast}},$

we deduce

$\displaystyle\mathbb{E}[\Delta^{\delta}_{k+1}]$
$\displaystyle\leq\mathbb{E}[\Delta^{\delta}_{k}]-p\mu_{k+1}C_{k+1}\mathbb{E}[\Psi({\boldsymbol{\mathsf{x}}}_{k}^{\delta})]+\delta\mu_{k+1}p^{1/p^{\ast}}\frac{1}{N}\sum_{i=1}^{N}\Big{(}\mathbb{E}[\Psi_{i}({\boldsymbol{\mathsf{x}}}_{k}^{\delta})]\Big{)}^{1/p^{\ast}}.$
(27)

Next we remove the exponent in the last term. Using Young’s inequality
$ab\leq\frac{a^{p}}{p}\omega^{-p}+\frac{b^{p^{\ast}}}{p^{\ast}}\omega^{p^{\ast}}$,
with $a=\delta$ and
$b=\mathbb{E}[\Psi_{i}({\boldsymbol{\mathsf{x}}}_{k}^{\delta})]^{1/p^{\ast}}$,
we have

$\displaystyle\frac{1}{N}\sum_{i=1}^{N}\delta\Big{(}\mathbb{E}[\Psi_{i}({\boldsymbol{\mathsf{x}}}_{k}^{\delta})]\Big{)}^{1/p^{\ast}}\leq\delta^{p}\frac{\omega^{-p}}{p}+\mathbb{E}\Big{[}\frac{1}{N}\sum_{i=1}^{N}\Psi_{i}({\boldsymbol{\mathsf{x}}}_{k}^{\delta})\Big{]}\frac{\omega^{p^{\ast}}}{p^{\ast}}\leq\delta^{p}\frac{\omega^{-p}}{p}+\mathbb{E}[\Psi({\boldsymbol{\mathsf{x}}}_{k}^{\delta})]\frac{\omega^{p^{\ast}}}{p^{\ast}}.$

Plugging this back in (27) gives

$\displaystyle\mathbb{E}[\Delta^{\delta}_{k+1}]$
$\displaystyle\leq\mathbb{E}[\Delta^{\delta}_{k}]-p\mu_{k+1}C_{k+1}\mathbb{E}[\Psi({\boldsymbol{\mathsf{x}}}_{k}^{\delta})]+p^{1/p^{\ast}}(p^{\ast})^{-1}\omega^{p^{\ast}}\mu_{k+1}\mathbb{E}[\Psi({\boldsymbol{\mathsf{x}}}_{k}^{\delta})]+p^{-1/p}\delta^{p}\omega^{-p}\mu_{k+1}.$

Taking $\omega>0$ small enough so that $\omega^{p^{\ast}}\leq
p^{\ast}p^{1/p}C_{k}$ (which can be made uniformly on $k$, thanks to the
positive lower bound on $C_{k}$), replacing $k+1$ with $k(\delta)$ and using
the inductive argument, we have

$\displaystyle\mathbb{E}[\Delta^{\delta}_{k(\delta)}]$
$\displaystyle\leq\mathbb{E}[\Delta^{\delta}_{k(\delta)-1}]+p^{-1/p}\omega^{-p}\delta^{p}\mu_{k(\delta)}\leq\mathbb{E}[\Delta^{\delta}_{k_{\epsilon}}]+p^{-1/p}\omega^{-p}\delta^{p}\sum_{\ell=1}^{k(\delta)}\mu_{\ell}.$

Since $\lim_{\delta\searrow 0}\delta^{p}\sum_{\ell=1}^{k(\delta)}\mu_{\ell}=0$
and $\lim_{\delta\searrow 0}k(\delta)=\infty$, there exists
$j_{2}\in\mathbb{N}$ such that for all $j\geq j_{2}$ we have
$k(\delta_{j})\geq k_{\epsilon}$ and
$p^{-1/p}\omega^{-p}\delta_{j}^{p}\sum_{\ell=1}^{k(\delta_{j})}\mu_{\ell}<\epsilon$.
Taking ${j_{\epsilon}}=j_{1}\vee j_{2}$ shows
$\mathbb{E}[\Delta_{k(\delta_{j})}^{\delta_{j}}]<4\epsilon$ for all $j\geq
j_{\epsilon}$, and hence the desired claim follows. ∎

######  Remark 4.4.

In the constant step-size regime, such as in the case of conditionally stable
operators, the correspondence between the noise level and the step-size regime
takes a more standard form. Namely, the condition in Theorem 4.3 reduces to
$\lim_{\delta\searrow 0}\delta^{p}k(\delta)=0$. In other words, we have
$k(\delta)={\cal O}(\delta^{-p})$, mirroring the traditional conditions in
Euclidean spaces. Note that the condition on $k(\delta)$ is fairly broad, and
does not give useful concrete stopping rules directly. Generally, the issue of
a posterior stopping rules for stochastic iterative methods is completely
open, even for the Hilbert setting [22]. For a polynomially decaying step-
sizes $\mu_{k}=c_{0}k^{-\beta}$, the conditions $\frac{1}{p^{*}}<\beta\leq 1$
and
$c_{0}<(\frac{P^{\ast}}{L_{\max}^{p^{\ast}}G_{p^{\ast}}})^{\frac{1}{p^{\ast}-1}}$
give a valid step-size choice, and the stopping index $k(\delta)$ should
satisfy $\lim_{\delta\searrow 0}k(\delta)=\infty$ and $\lim_{\delta\searrow
0}k(\delta)\delta^{\frac{p}{1-\beta}}=0$.

######  Remark 4.5.

It is of much interest to derive a convergence rate for noisy data under a
conditional stability condition as in Theorem 3.14, as a natural extension of
the regularising property. However, this is still unavailable. Within the
current analysis strategy, deriving the rate would require quantitative
versions of stability estimates in Lemma 4.2 in terms of $\delta$ and $k$.
Generally the convergence rate analysis for iterative regularisation methods
in Banach space remains a very challenging task, and much more work is still
needed.

##  5 Numerical experiments

We present numerical results on two sets of experiments to illustrate distinct
features of the SGD (10). The first set of experiments deals with an integral
operator and the reconstruction of a sparse signal in the presence of either
Gaussian or impulse noise. On this model example, we investigate the impact of
the number of batches and the choice of the spaces ${\cal X}$ and ${\cal Y}$
on the performance of the algorithm. To simplify the study we investigate
spaces ${\cal X}$ and ${\cal Y}$ that are smooth and convex of power type, and
thus the corresponding duality maps are singletons. To facilitate a direct
comparison of the SGD with the Landweber method, we count the computational
complexity with respect to the number of epochs, i.e. the size $N_{b}$ of
partition defined below. Note moreover that our implementation of the
Landweber method does not use the stepsizes described in [44, Method 3.1],
since the latter requires knowledge of quantities that are inconvenient to
compute in practice. The second set of experiments is about tomographic
reconstruction, with respect to different types of noise. All the shown
reconstructions are obtained with a single stochastic run, as is often done in
practice, and the stopping index is determined in a trial and error manner so
that the corresponding reconstruction yields small errors.

###  5.1 Model linear inverse problem

First we consider the following model inverse problem studied in [28]. Let
$\kappa:\overline{\Omega}\times\overline{\Omega}\rightarrow\mathbb{R}^{+}$,
with $\Omega=(0,1)$, be a continuous function, and define an integral operator
${\cal T}_{\kappa}:{\cal L}^{r_{{\cal X}}}(\Omega)\rightarrow{\cal
L}^{r_{{\cal Y}}}(\Omega)$, for $1<r_{{\cal X}},r_{{\cal Y}}<\infty$, by

$\displaystyle({\cal T}_{\kappa}x)(t)=\int_{\Omega}\kappa(t,s)x(s)ds.$ (28)

This is a compact linear operator between ${\cal L}^{r_{{\cal X}}}(\Omega)$
and ${\cal L}^{r_{{\cal Y}}}(\Omega)$, with the adjoint ${\cal
T}_{\kappa}^{\ast}\\!:\\!{\cal L}^{r^{\ast}_{{\cal
Y}}}(\Omega)\\!\rightarrow\\!{\cal L}^{r^{\ast}_{{\cal X}}}(\Omega)$ given by
$({\cal T}_{\kappa}^{\ast}y)(s)\\!=\\!\int_{\Omega}\kappa(t,s)y(t)dt$. To
approximate the integrals, we subdivide the interval $\overline{\Omega}$ into
$N\\!=\\!1000$ subintervals $[\frac{k}{N},\frac{k+1}{N}]$, for
$k\\!=\\!0,\\!\ldots\\!,N\\!-\\!1$, and then use quadrature, giving a finite-
dimensional model
${\mathbf{A}}{\boldsymbol{\mathsf{x}}}\\!=\\!{\boldsymbol{\mathsf{y}}}$, with
${\mathbf{A}}\\!=\\!\frac{1}{N}\\!\left(\kappa\\!\left(\frac{j-1}{N},\frac{2k-1}{N}\right)\\!\right)_{j,k=1}^{N}$
and
${\boldsymbol{\mathsf{x}}}\\!=\\!\left(x\left(\frac{2j-1}{2N}\right)\right)_{j=1}^{N}$.
For SGD we use $N_{b}\in[N]$ mini-batches. To obtain equisized batches, we
assume that $N_{b}$ divides $N$. The mini-batch matrices ${\mathbf{A}}_{j}$
are then constructed by taking every $N_{b}$-th row of ${\mathbf{A}}$, shifted
by $j$, resulting in well-balanced mini-batches, in the sense that the norm
$\|{\mathbf{A}}_{j}\|$ is (nearly) independent of $j$.

The kernel function $k(t,s)$ and the exact signal $x^{\dagger}$ are defined
respectively by

$\displaystyle\kappa(t,s)=\begin{cases}40t(1-s),&\text{ if }t\leq s,\\\
40s(1-t),&\text{ otherwise},\end{cases}\quad\mbox{and}\quad
x^{\dagger}(s)=\begin{cases}1,&\text{ if
}s\in[\frac{9}{40},\frac{11}{40}]\cup[\frac{29}{40},\frac{31}{40}],\\\
2,&\text{ if }s\in[\frac{19}{40},\frac{21}{40}],\\\ 0,&\text{
otherwise}.\end{cases}$

This is a sparse signal and we expect sparsity promoting norms to perform
well. To illustrate this, we compare the following four settings: (a) ${\cal
X}={\cal Y}={\cal L}^{2}(\Omega)$; (b) ${\cal X}={\cal L}^{2}(\Omega)$ and
${\cal Y}={\cal L}^{1.1}(\Omega)$; (c) ${\cal X}={\cal L}^{1.5}(\Omega)$ and
${\cal Y}={\cal L}^{2}(\Omega)$; (d) ${\cal X}={\cal L}^{1.1}(\Omega)$ and
${\cal Y}={\cal L}^{2}(\Omega)$. Setting (a) is the standard Hilbert space
setting, suitable for recovering smooth solutions from measurement data with
i.i.d. Gaussian noise, whereas settings (b)-(d) use Banach spaces. Settings
(c) and (d) both aim at sparse solutions, and we expect the latter to yield
sparser solutions, since spaces ${\cal L}^{r}(\Omega)$ progressively enforce
sparser solutions as the exponent $r$ gets closer to $1$. In the experiments,
we employ the step-size schedule
$\mu_{k}=\frac{L_{\max}}{1+0.05(k/N_{b})^{1/p^{\ast}+0.01}},$ with
$L_{\max}=\max_{j\in[N_{b}]}\|{\mathbf{A}}_{j}\|$. This satisfies the
summability conditions $\sum_{k=1}^{\infty}\mu_{k}=\infty$ and
$\sum_{k=1}^{\infty}\mu_{k}^{p^{*}}<\infty$ required by Theorem 3.8. The
operator norm $\|{\mathbf{A}}_{j}\|=\|{\mathbf{A}}_{j}\|_{{\cal L}^{r_{{\cal
X}}}\rightarrow{\cal L}^{r_{{\cal Y}}}}=\max_{{\boldsymbol{\mathsf{x}}}\neq
0}\frac{\|{\mathbf{A}}_{j}{\boldsymbol{\mathsf{x}}}\|_{{\cal L}^{r_{{\cal
Y}}}}}{\|{\boldsymbol{\mathsf{x}}}\|_{{\cal L}^{r_{{\cal X}}}}}$ is estimated
using Boyd’s power method [3]. All the reconstruction algorithms are
initialised with a zero vector.

In Fig. 1, we compare the reconstructions with settings (a)-(d) for exact
data. We observe from Fig. 1(a) that settings (a) and (b), with ${\cal
X}={\cal L}^{2}(\Omega)$, result in smooth solutions that fail to capture the
sparsity structure of the true signal ${\boldsymbol{\mathsf{x}}}^{\dagger}$.
In contrast, the choice ${\cal X}={\cal L}^{1.5}(\Omega)$ recovers a sparser
solution, and the choice ${\cal X}={\cal L}^{1.1}(\Omega)$ gives a truly
sparse reconstruction, but with peaks that overshoot the magnitude of
${\boldsymbol{\mathsf{x}}}^{\dagger}$. This might be related to the fact
${\boldsymbol{\mathsf{x}}}^{\dagger}$ exhibits a cluster structure in addition
to sparsity, which is not accounted for in the choice of the space ${\cal
X}={\cal L}^{1.1}(\Omega)$ [52, 23]. Fig. 1(b) indicates that early stopping
would result in lower peaks and significantly reduce the overshooting, but a
more explicit form of regularisation [52, 23] might allow faster convergence.

(a) Changing ${\cal X}$ and ${\cal Y}$ for $N_{b}=100$ | (b) Progression of iterates for ${\cal X}={\cal L}^{1.1}(\Omega)$
---|---
Figure 1: Comparison of reconstructed solutions after $500$ epochs.

In Fig. 2, we investigate the convergence of the objective value with respect
to the number of batches $N_{b}$ and the choice of the solution space ${\cal
X}$. As expected, having a larger number of batches results in a faster
initial convergence, but also in increased variance, as shown by the
oscillations. Moreover, the variance is lower in the case of a smoother space
${\cal X}$ (promoting smoother solutions), where the variance existing in
early epochs is dramatically reduced later on. This observation can be
explained by the gradient expression
$g({\boldsymbol{\mathsf{x}}},{\boldsymbol{\mathsf{y}}},i)={\mathbf{A}}_{i}^{\ast}\jmath_{p}^{{\cal
Y}}({\mathbf{A}}_{i}{\boldsymbol{\mathsf{x}}}-{\boldsymbol{\mathsf{y}}}_{i})$,
which tends to zero as SGD iterates converge to the true solution
${\boldsymbol{\mathsf{x}}}^{\dagger}$ and so does its variance, and the larger
is the exponent $p$, the faster is the convergence.

(a) ${\cal X}={\cal L}^{1.1}(\Omega)$ and ${\cal Y}={\cal L}^{2}(\Omega)$ | (b) ${\cal X}={\cal L}^{1.5}(\Omega)$ and ${\cal Y}={\cal L}^{2}(\Omega)$
---|---
Figure 2: The variation of
$\frac{1}{p}\sum_{i=1}^{N_{b}}\|{\mathbf{A}}_{i}{\boldsymbol{\mathsf{x}}}_{k}-{\boldsymbol{\mathsf{y}}}_{i}\|_{{\cal
Y}}^{p}$ with respect to the number of batches $N_{b}$.

Next we examine the performance of the algorithm when the observational data
${\boldsymbol{\mathsf{y}}}^{\delta}$ contains (random-valued) impulse noise,
cf. Fig. 3, which is generated by

$\displaystyle y^{\delta}_{i}=\left\\{\begin{aligned}
y_{i}^{\dagger},&\quad\mbox{with probability }1-p,\\\
(1-\xi)y_{i}^{\dagger},&\quad\mbox{with probability }p/2,\\\
1.4\xi+(1-\xi)y_{i}^{\dagger},&\quad\mbox{with probability
}p/2,\end{aligned}\right.$

where $p\in(0,1)$ denotes the percentage of corruption (which is set to $0.05$
in the experiment) and $\xi\sim{\rm Uni}(0.1,0.4)$ follows a uniform
distribution over the interval $(0.1,0.4)$. It is known that
$\mathcal{L}^{r}(\Omega)$ fittings with $r$ close to 1 is suitable for
impulsive noise. This allows investigating the role of not only the space
${\cal X}$ but also ${\cal Y}$. The results in Fig. 3(b) show that the choice
${\cal Y}={\cal L}^{r_{{\cal Y}}}(\Omega)$, with $r_{{\cal Y}}$ close to $1$,
performs significantly better. Indeed, the Hilbert setting ${\cal X}={\cal
Y}={\cal L}^{2}(\Omega)$ produces overly smooth, non-sparse solutions with
pronounced artefacts. In sharp contrast, setting ${\cal X}={\cal Y}={\cal
L}^{1.1}(\Omega)$ yields solutions that can correctly identify the sparsity
structure of the true solution, and have no artefacts. Similar as before, the
reconstruction in this setting overestimates the signal magnitude on its
support, which is exacerbated as the exponent $r_{{\cal Y}}$ gets closer to
$1$.

(a) Data with impulse noise | (b) Reconstructions with respect to ${\cal X}$ and ${\cal Y}$
---|---
Figure 3: The reconstruction performance in case of impulse noise. The
algorithms utilised $N_{b}=100$ batches and were run for $250$ epochs.

Lastly, we investigate the convergence behaviour of the method for the
generalised model (15) in Section 3.2, where stochastic directions
$g({\boldsymbol{\mathsf{x}}},{\boldsymbol{\mathsf{y}}},i)$ are defined as
$g({\boldsymbol{\mathsf{x}}},{\boldsymbol{\mathsf{y}}},i)={\mathbf{A}}_{i}^{\ast}\jmath_{q}^{{\cal
Y}}({\mathbf{A}}_{i}{\boldsymbol{\mathsf{x}}}-{\boldsymbol{\mathsf{y}}}_{i})$,
with $q=r_{{\cal Y}}$ different from the convexity parameter $p$ of the space
${\cal X}$. The results in Fig. 4 show that this can indeed be beneficial for
the performance of the method: the reconstructions are more accurate not only
in terms of the solution support, but also in terms of the magnitudes of the
non-zero entries. However, the precise mechanism of the excellent performance
remains largely elusive.

(a) Standard vs generalised Kaczmarz | (b) Changing ${\cal X}$ in generalised Kaczmarz
---|---
Figure 4: The dependence of the reconstructions in the case of impulse noise
on the choice of $q$ parameter in the generalised model (15). The results are
obtained using $N_{b}=100$ batches, after $250$ epochs.

###  5.2 Computed Tomography

Now we numerically investigate the behaviour of SGD on computed tomography
(CT), with respect to the model spaces ${\cal X}$ and ${\cal Y}$ and data
noise. In CT reconstruction, we aim at determining the density of cross
sections of an object by measuring the attenuation of X-rays as they propagate
through the object [36]. Mathematically, the forward map is given by the Radon
transform. In the experiments, the discrete forward operator ${\mathbf{A}}$ is
defined by a $2D$ parallel beam geometry, with $180$ projection angles on a
$1$ angle separation, $256$ detector elements, and pixel size of $0.1$. The
sought-for signal ${\boldsymbol{\mathsf{x}}}^{\dagger}$ is a (sparse) phantom,
cf. Fig. 5(a). After applying the forward operator ${\mathbf{A}}$, either
Gaussian (with mean zero and variance $0.01$) or salt-and-pepper noise is
added. In the latter setting we consider low (with $5\%$ of values changed to
either salt or pepper values) and high ($10\%$ of values changed) noise
regimes. The resulting sinograms (i.e. measurement data) are shown in Fig.
5(b)-(d). Note that standard quality metrics in image assessment, such as peak
signal to noise ratio or mean squared error, are computed using the distance
between images in the $\ell^{2}$-norm, which have an implicit bias towards
Hilbert spaces and smooth signals, whereas using a metric that emphasises
sparsity is more pertinent to sparsity promoting spaces. To provide a balanced
comparison, we report the following two metrics based on normalised
$\ell^{1}$\- and $\ell^{2}$-norms:
$\delta_{1}({\boldsymbol{\mathsf{x}}})=\|{\boldsymbol{\mathsf{x}}}^{\dagger}-{\boldsymbol{\mathsf{x}}}\|_{\ell^{1}}/\|{\boldsymbol{\mathsf{x}}}^{\dagger}\|_{\ell^{1}}$
and
$\delta_{2}({\boldsymbol{\mathsf{x}}})=\|{\boldsymbol{\mathsf{x}}}^{\dagger}-{\boldsymbol{\mathsf{x}}}\|_{\ell^{2}}/\|{\boldsymbol{\mathsf{x}}}^{\dagger}\|_{\ell^{2}}.$

(a) Original phantom | (b) Gaussian noise measurement
---|---
(c) Low noise salt-and-pepper measurement | (d) High noise salt-and-pepper measurement
Figure 5: The plot in (a) shows the phantom to be recovered and (b)-(d) show
noisy measurements used in the recovery: in (b), random Gaussian noise was
added, and (c)-(d) are sinogram data degraded by salt-and-pepper noise in the
low ($5\%$) and high ($10\%$) noise regimes.

First, we show the performance on Gaussian noise, where we compare the Hilbert
setting (${\cal X}={\cal Y}={\cal L}^{2}$) with two Banach settings (${\cal
X}={\cal L}^{1.1}$, ${\cal Y}={\cal L}^{2}$, and ${\cal X}={\cal Y}={\cal
L}^{1.1}$). In the reconstruction, we employ step-sizes
$\mu_{k}=\frac{L_{\max}/2}{1+0.05(k/N_{b})^{1/p^{\ast}+0.01}}$, with
$L_{\max}=\max_{j\in[N_{b}]}\|{\mathbf{A}}_{j}\|$. Fig. 6 shows exemplary
reconstructions. In all three settings much of the noise is retained in the
reconstruction, and whereas the Hilbert setting is better at recovering the
magnitude of non-zero entries, the Banach settings are better at recovering
the support. Moreover, we observe that the Banach setting with a sparse signal
space ${\cal X}={\cal L}^{1.1}$, and a smooth observation space ${\cal
Y}={\cal L}^{2}$, has the best performance in terms of $\delta_{1}$ and
$\delta_{2}$ metrics. The Hilbert model performs better than the fully sparse
model ${\cal X}={\cal Y}={\cal L}^{1.1}$ in terms of the smooth metric
$\delta_{2}$, but worse in the sparsity promoting metric ($\delta_{1}$). We
also consider the Banach setting for the generalised model (15), with ${\cal
X}={\cal Y}={\cal L}^{1.1}$ and $p_{{\cal Y}}=1.1$, where we study the effects
of early stopping. Fig. 7 shows that this setting recovers the support more
accurately (and actually does so very early on) and recovers the magnitudes
better, but that a form of regularisation (through e.g. early stopping) can be
beneficial, since in the later epochs SGD iterates again tend to overshoot on
the support. A similar behaviour can observed for other studied Banach space
settings, but not for the Hilbert space setting, which does not recover the
support.

(a) ${\cal X}={\cal Y}={\cal L}^{2}$ | (b)${\cal X}={\cal L}^{1.1}$, ${\cal Y}={\cal L}^{2}$ | (c) ${\cal X}={\cal Y}={\cal L}^{1.1}$
---|---|---
$\delta_{1}({\boldsymbol{\mathsf{x}}})/\delta_{2}({\boldsymbol{\mathsf{x}}})$: $2.643/0.528$ |  $\delta_{1}({\boldsymbol{\mathsf{x}}})/\delta_{2}({\boldsymbol{\mathsf{x}}})$: $0.711/0.341$ |  $\delta_{1}({\boldsymbol{\mathsf{x}}})/\delta_{2}({\boldsymbol{\mathsf{x}}})$: $2.195/0.620$
Figure 6: The reconstruction of the phantom from the observed sinograms degraded by Gaussian noise, cf. Fig. 5(b). The algorithms use $N_{b}=60$ batches and were run for $200$ epochs. (a) $5$ epochs | (b) $50$ epochs | (c) $200$ epochs
---|---|---
$\delta_{1}({\boldsymbol{\mathsf{x}}})/\delta_{2}({\boldsymbol{\mathsf{x}}})$: $0.702/0.627$ |  $\delta_{1}({\boldsymbol{\mathsf{x}}})/\delta_{2}({\boldsymbol{\mathsf{x}}})$: $0.263/0.235$ |  $\delta_{1}({\boldsymbol{\mathsf{x}}})/\delta_{2}({\boldsymbol{\mathsf{x}}})$: $0.604/0.308$
Figure 7: The evolution of the quality of reconstruction from sinograms
degraded by Gaussian noise with respect to the number of epochs. The algorithm
uses ${\cal X}={\cal Y}={\cal L}^{1.1}$ and $p_{{\cal Y}}=1.1$, with
$N_{b}=60$ batches.

We next investigate the performance for low and high salt-and-pepper noise. We
compare the Hilbert setting with two Banach settings: the standard SGD with
${\cal X}={\cal Y}={\cal L}^{1.1}$ and the generalised model (15) with ${\cal
X}={\cal Y}={\cal L}^{1.1}$ and $p_{{\cal Y}}=1.1$. For the reconstruction, we
employ step-sizes $\mu_{k}=\frac{0.5}{1+0.05(k/N_{b})^{1/p^{\ast}+0.01}}$. The
results in Fig. 8 show the reconstructions after $200$ epochs with $N_{b}=60$
batches. In the low noise regime, the Hilbert setting can reconstruct the
general shape of the phantom, but retains a lot of the noise and exhibits
streaking artefacts in the background. The reconstruction in the high noise
regime is of much poorer quality. The standard Banach SGD shows good behaviour
in the low-noise setting, reconstructing well both the sparsity structure and
the magnitudes, but its performance degrades in the high noise setting. In
sharp contrast, the model (15) shows a nearly perfect reconstruction
performance - the phantom is well recovered, with intensities on the correct
scale, for both low and high noise regimes. Similar as before, we observe that
Banach methods tend to slightly overestimate the overall intensities, though
the recovered values are comparable to the true solution. Overall, the Hilbert
setting shows a qualitatively worst performance, in both $\ell^{1}$\- and
$\ell^{2}$-norm sense, and the model (15) shows the best performance.

(a) ${\cal X}={\cal Y}={\cal L}^{2}$ in low noise | (b)${\cal X}={\cal Y}={\cal L}^{1.1}$ in low noise | (c) ${\cal X}={\cal Y}={\cal L}^{1.1}$, $p_{{\cal Y}}=1.1$ in low noise
---|---|---
$\delta_{1}({\boldsymbol{\mathsf{x}}})/\delta_{2}({\boldsymbol{\mathsf{x}}})$: $18.67/3.71$ | $\delta_{1}({\boldsymbol{\mathsf{x}}})/\delta_{2}({\boldsymbol{\mathsf{x}}})$: $1.80/0.544$ | $\delta_{1}({\boldsymbol{\mathsf{x}}})/\delta_{2}({\boldsymbol{\mathsf{x}}})$: $2.43/3.68\cdot\text{e-}3$
(a) ${\cal X}={\cal Y}={\cal L}^{2}$ in high noise | (b) ${\cal X}={\cal Y}={\cal L}^{1.1}$ in high noise | (c) ${\cal X}={\cal Y}={\cal L}^{1.1}$, $p_{{\cal Y}}=1.1$ in high noise
$\delta_{1}({\boldsymbol{\mathsf{x}}})/\delta_{2}({\boldsymbol{\mathsf{x}}})$: $26.61/5.19$ |  $\delta_{1}({\boldsymbol{\mathsf{x}}})/\delta_{2}({\boldsymbol{\mathsf{x}}})$: $5.37/1.54$ |  $\delta_{1}({\boldsymbol{\mathsf{x}}})/\delta_{2}({\boldsymbol{\mathsf{x}}})$: $3.72/6.03\cdot\text{e-}3$
Figure 8: The reconstruction of the phantom from the observed sinograms,
degraded with low (top) and high (bottom) salt-and-pepper noise, respectively,
obtained using the Hilbert space model (${\cal X}={\cal Y}={\cal L}^{2}$)
(left), the Banach model (${\cal X}={\cal Y}={\cal L}^{1.1}$) (middle) and the
Banach model (${\cal X}={\cal Y}={\cal L}^{1.1}$) with the generalised
Kaczmarz scheme ($p_{{\cal Y}}=1.1$) (right). The algorithms use $N_{b}=60$
batches and were run for $200$ epochs.

Lastly, we investigate a more challenging setting with noise affecting not
only the sinograms, but also the original phantoms. Then the ground-truth
image is only approximately sparse. The phantom is degraded with Gaussian
noise (zero mean and variance $0.01$) after which we apply the forward
operator to the resulting noisy phantom. We then add either Gaussian (zero
mean and variance $0.01$) or salt-and-pepper noise (affecting $3\%$ of
measurements); see Fig. 9 for representative images. The reconstruction
algorithms use SGD with a decaying step-size schedule,
$\mu_{k}=\frac{0.2}{1+0.05(k/N_{b})^{1/p^{\ast}+0.01}}$.

(a) Noisy Phantom | (b) Gaussian measurement noise | (c) Salt-and-pepper measurement noise
---|---|---
Figure 9: The phantoms and sinograms for the forward model with both pre and
post measurement noise. The phantom on the left is degraded by Gaussian noise.
After applying the forward operator, either Gaussian (middle) or salt-and-
pepper noise (right) is added to the sinogram.

The reconstructions for data with Gaussian noise in both phantom and sinogram
are shown in Fig. 10. As before, reconstructions in the Hilbert setting are
comparable, but slightly worse than that with the Banach ones. Banach methods
are better at recovering the sparsity structure of the solution, and have
better reconstruction quality metrics, though they do not completely remove
the noise. In the second setting, with the Gaussian noise affecting the
phantom and salt-and-pepper noise affecting the sinogram, the difference in
reconstruction quality in the Hilbert space and Banach space settings is
significantly more pronounced, cf. Fig. 11. In both settings, the choice of
spaces ${\cal X}$ and ${\cal Y}$ can have a big impact on the reconstruction
quality, especially on the amount of noise retained in the background.
Moreover, further improvements can be achieved by explicitly penalising the
objective function.

(a) ${\cal X}={\cal Y}={\cal L}^{2}$ | (b) ${\cal X}={\cal Y}={\cal L}^{1.1}$, $p_{{\cal Y}}=1.1$ | (c) ${\cal X}={\cal L}^{1.1}$, ${\cal Y}={\cal L}^{1.9}$, $p_{{\cal Y}}=1.9$
---|---|---
$\delta_{1}({\boldsymbol{\mathsf{x}}})/\delta_{2}({\boldsymbol{\mathsf{x}}})$: $5.65/1.12$ |  $\delta_{1}({\boldsymbol{\mathsf{x}}})/\delta_{2}({\boldsymbol{\mathsf{x}}})$: $3.16/0.632$ |  $\delta_{1}({\boldsymbol{\mathsf{x}}})/\delta_{2}({\boldsymbol{\mathsf{x}}})$: $2.99/0.561$
Figure 10: The reconstruction of the phantom from the observed sinograms with pre- and post-measurement Gaussian noise. The algorithms use $N_{b}=60$ batches and were run for $200$ epochs. (a) ${\cal X}={\cal Y}={\cal L}^{2}$ |  (b)${\cal X}={\cal Y}={\cal L}^{1.3}$, $p_{{\cal Y}}=1.3$ |  (c) ${\cal X}={\cal Y}={\cal L}^{1.1}$, $p_{{\cal Y}}=1.1$
---|---|---
$\delta_{1}({\boldsymbol{\mathsf{x}}})/\delta_{2}({\boldsymbol{\mathsf{x}}})$: $17.48/3.52$ |  $\delta_{1}({\boldsymbol{\mathsf{x}}})/\delta_{2}({\boldsymbol{\mathsf{x}}})$: $3.46/0.86$ |  $\delta_{1}({\boldsymbol{\mathsf{x}}})/\delta_{2}({\boldsymbol{\mathsf{x}}})$: $3.14/0.62$
Figure 11: The reconstructed phantom from the sinograms with a Gaussian pre-
measurement and a salt-and-pepper (post-)measurement noise. The algorithms use
$N_{b}=60$ batches and were run for $400$ epochs.

## Acknowledgements

We are very grateful to three anonymous referees for their constructive
comments which have led to a significant improvement of the quality of the
paper.

##  Appendix A Technical results and proofs

######  Lemma A.1 ([38, Lemma 6]).

Let $(\delta_{n})_{n}$ be a sequence of non-negative scalars, $(\mu_{n})_{n}$
a sequence of positive scalars, and $\alpha>0$. If

$\delta_{n+1}\leq\delta_{n}-\mu_{n+1}\delta_{n}^{1+\alpha},\text{ for all
}n=0,\ldots,N,$

then

$\delta_{N}\leq{\delta_{0}}{\Big{(}1+\alpha\delta_{0}^{\alpha}\sum_{n=1}^{N}\mu_{n}\Big{)}^{-1/\alpha}}.$

###  A.1 Two elementary estimates

In this section, we present two elementary estimates on the SGD iterates for
exact data that are useful in establishing the regularising property.

######  Lemma A.2.

Let the sequence $({\boldsymbol{\mathsf{x}}}_{k})_{k\in\mathbb{N}}$ be
generated by iterations (10), and let the step-sizes
$(\mu_{k})_{k\in\mathbb{N}}$ satisfy
$\mu_{k}^{p^{\ast}-1}\leq\frac{p^{\ast}}{G_{p^{\ast}}L_{\max}^{p^{\ast}}}$ for
all $k\in\mathbb{N}$, and stochastic update directions
${\boldsymbol{\mathsf{g}}}_{k}$ be of the form (11). Then for any
$\widehat{\boldsymbol{\mathsf{x}}}\in{\cal X}_{\min}$, the sequence
$({\mathbf{B}}_{p}({\boldsymbol{\mathsf{x}}}_{k},\widehat{\boldsymbol{\mathsf{x}}}))_{k\in\mathbb{N}}$
is non-increasing. In particular, if
${\mathbf{B}}_{p}({\boldsymbol{\mathsf{x}}}_{0},\widehat{\boldsymbol{\mathsf{x}}})\leq\rho$,
then
${\mathbf{B}}_{p}({\boldsymbol{\mathsf{x}}}_{k},\widehat{\boldsymbol{\mathsf{x}}})\leq\rho$
for all $k$.

###### Proof.

Let
$\Delta_{k}={\mathbf{B}}_{p}({\boldsymbol{\mathsf{x}}}_{k},\widehat{\boldsymbol{\mathsf{x}}})$.
By Lemma 3.6, we have

$\displaystyle\Delta_{k+1}$
$\displaystyle\leq\Delta_{k}-\mu_{k+1}\left<{{\boldsymbol{\mathsf{g}}}_{k+1}},{{\boldsymbol{\mathsf{x}}}_{k}-\widehat{\boldsymbol{\mathsf{x}}}}\right>+\frac{G_{p^{\ast}}}{p^{\ast}}\mu_{k+1}^{p^{\ast}}\|{\boldsymbol{\mathsf{g}}}_{k+1}\|_{{\cal
X}^{\ast}}^{p^{\ast}}.$

By the definition of duality map and the choice of the update directions
${\boldsymbol{\mathsf{g}}}_{k}$, we have

$\displaystyle\left<{{\boldsymbol{\mathsf{g}}}_{k+1}},{{\boldsymbol{\mathsf{x}}}_{k}-\widehat{\boldsymbol{\mathsf{x}}}}\right>$
$\displaystyle=\left<{\jmath_{p}^{{\cal
Y}}({\mathbf{A}}_{i_{k+1}}{\boldsymbol{\mathsf{x}}}_{k}-{\boldsymbol{\mathsf{y}}}_{i_{k+1}})},{{\mathbf{A}}_{i_{k+1}}{\boldsymbol{\mathsf{x}}}_{k}-{\boldsymbol{\mathsf{y}}}_{i_{k+1}}}\right>=\|{\mathbf{A}}_{i_{k+1}}{\boldsymbol{\mathsf{x}}}_{k}-{\boldsymbol{\mathsf{y}}}_{i_{k+1}}\|_{{\cal
Y}}^{p},$ $\displaystyle\|{\boldsymbol{\mathsf{g}}}_{k+1}\|_{{\cal
X}^{\ast}}^{p^{\ast}}$
$\displaystyle=\|{\mathbf{A}}_{i_{k+1}}^{\ast}\jmath_{p}^{{\cal
Y}}({\mathbf{A}}_{i_{k+1}}{\boldsymbol{\mathsf{x}}}_{k}-{\boldsymbol{\mathsf{y}}}_{i_{k+1}})\|_{{\cal
X}^{\ast}}^{p^{\ast}}\leq\|{\mathbf{A}}_{i_{k+1}}^{\ast}\|^{p^{\ast}}\|\jmath_{p}^{{\cal
Y}}({\mathbf{A}}_{i_{k+1}}{\boldsymbol{\mathsf{x}}}_{k}-{\boldsymbol{\mathsf{y}}}_{i_{k+1}})\|_{{\cal
Y}^{\ast}}^{p^{\ast}}$ $\displaystyle\leq
L_{\max}^{p^{\ast}}\|{\mathbf{A}}_{i_{k+1}}{\boldsymbol{\mathsf{x}}}_{k}-{\boldsymbol{\mathsf{y}}}_{i_{k+1}}\|_{{\cal
Y}}^{(p-1)p^{\ast}}=L_{\max}^{p^{\ast}}\|{\mathbf{A}}_{i_{k+1}}{\boldsymbol{\mathsf{x}}}_{k}-{\boldsymbol{\mathsf{y}}}_{i_{k+1}}\|_{{\cal
Y}}^{p}.$

Consequently,

$\displaystyle\Delta_{k+1}$
$\displaystyle\leq\Delta_{k}-\mu_{k+1}\left<{{\boldsymbol{\mathsf{g}}}_{k+1}},{{\boldsymbol{\mathsf{x}}}_{k}-\widehat{\boldsymbol{\mathsf{x}}}}\right>+\frac{G_{p^{\ast}}}{p^{\ast}}\mu_{k+1}^{p^{\ast}}\|{\boldsymbol{\mathsf{g}}}_{k+1}\|_{{\cal
X}^{\ast}}^{p^{\ast}}$
$\displaystyle\leq\Delta_{k}-\Big{(}1-L_{\max}^{p^{\ast}}\frac{G_{p^{\ast}}}{p^{\ast}}\mu_{k+1}^{p^{\ast}-1}\Big{)}\mu_{k+1}\|{\mathbf{A}}_{i_{k+1}}{\boldsymbol{\mathsf{x}}}_{k}-{\boldsymbol{\mathsf{y}}}_{i_{k+1}}\|_{{\cal
Y}}^{p}.$

Since
$\mu_{k}^{p^{\ast}-1}\leq\frac{p^{\ast}}{G_{p^{\ast}}L_{\max}^{p^{\ast}}}$ by
assumption, $\Delta_{k+1}\leq\Delta_{k}\leq\Delta_{0}$, completing the proof.
∎

######  Lemma A.3 (Coercivity of the Bregman distance).

If
$\Delta_{k}={\mathbf{B}}_{p}({\boldsymbol{\mathsf{x}}}_{k},{\boldsymbol{\mathsf{x}}}^{\dagger})\leq
C<\infty$ for all $k$, then $\|{\boldsymbol{\mathsf{x}}}_{k}\|_{{\cal
X}}^{p}\leq(2p^{\ast})^{p}(\|{\boldsymbol{\mathsf{x}}}^{\dagger}\|_{{\cal
X}}^{p}\vee C)$, for all $k\in\mathbb{N}$.

###### Proof.

By the definition of $\Delta_{k}$ and the Cauchy-Schwarz inequality, we have

$\Delta_{k}\geq\frac{1}{p^{\ast}}\|{\boldsymbol{\mathsf{x}}}_{k}\|_{{\cal