text:
\section{Introduction and problem setting}

Random block-coordinate descent algorithms are nowadays among the methods of choice for solving large scale optimization problems \cite{Nes12,Ric16,Wri15}. Indeed, they have low complexity and low memory requirements and, additionally, they are amenable for distributed and parallel implementations \cite{Ric16c,Ric16}.
In the last decade a number of works have appeared on the topic which address
several aspects, that is: the way the block sampling is performed, the composite structure, the partial separability, and the smoothness/geometrical properties of the objective function, accelerations, and iteration complexity  \cite{Com15,Com18,Fer15,LuX15,Nec16,Nec17,Nes12,Qu16,Qu16b,Ric15,Ric16,Tap18}.

In this work we consider the following optimization problem
\begin{equation}
\label{eq:20190110b}
\minimize{\ensuremath \xx \in \ensuremath {\mathsf{H}}}{\ensuremath{{\mathsf{f}}}(\ensuremath \xx)+\ensuremath{{\gs}}(\ensuremath \xx)},
\qquad\ensuremath{{\gs}}(\ensuremath \xx)= \sum_{i=1}^m \ensuremath{{\mathsf{h}}}_i(\ensuremath \mathsf{x}_i),
\end{equation}
where $\ensuremath {\mathsf{H}}$ is the direct sum of $m$ separable real Hilbert spaces
 $(\ensuremath \mathsf{H}_i)_{1 \leq i \leq m}$, that is,
\begin{equation*}
\ensuremath {\mathsf{H}} = \bigoplus_{i=1}^m \ensuremath \mathsf{H}_i,\qquad
(\forall\, \ensuremath \xx = (\ensuremath \mathsf{x}_i)_{1 \leq i \leq m}, \ensuremath \yy = (\ensuremath \mathsf{y}_i)_{1 \leq i \leq m} \in \ensuremath {\mathsf{H}})\quad
\scalarp{\ensuremath \xx,\ensuremath \yy} = \sum_{i=1}^m \scalarp{\ensuremath \mathsf{x}_i,\ensuremath \mathsf{y}_i},
\end{equation*}
and the following assumptions hold:
\begin{enumerate}[{\rm H1}]
\item\label{eq:A1}
$\ensuremath{{\mathsf{f}}}\colon \ensuremath {\mathsf{H}} \to \ensuremath \mathbb{R}$ is convex and differentiable,
\item\label{eq:A2}
for every $i=1, \dots, m$, $\ensuremath{{\mathsf{h}}}_i\colon \ensuremath \mathsf{H}_i \to \left]-\infty, +\infty\right]$ is  proper, convex, and lower semicontinuous.
\end{enumerate}
The objective of this study is a stochastic algorithm, called
\emph{parallel random block-coordinate forward-backward algorithm}, that depends
on a random variable $\ensuremath \varepsilon$ satisfying the following hypothesis
\begin{enumerate}[{\rm H3}]
\item\label{eq:A3} $\ensuremath \varepsilon = (\varepsilon_1,\dots, \varepsilon_m)$ is a random variable with values in $\{0,1\}^m$
such that, for every $i \in \{1,\dots,m\}$, $\ensuremath \mathsf{p}_i:=\ensuremath{\mathsf P}(\varepsilon_i = 1)>0$
and $\ensuremath{\mathsf P}\big(\ensuremath \varepsilon = (0,\dots, 0)\big) =0$.
\end{enumerate}
\begin{algorithm}
\label{algoRCD}
Let $(\ensuremath \varepsilon^\ensuremath n)_{\ensuremath n \in \ensuremath \mathbb{N}}
=  (\varepsilon^\ensuremath n_1, \dots, \varepsilon^\ensuremath n_m)_{\ensuremath n \in \ensuremath \mathbb{N}}$ be a sequence of
independent copies of $\ensuremath \varepsilon$.
Let $(\gamma_i)_{1 \leq i \leq m} \in \ensuremath \mathbb{R}_{++}^m$ and
 $\ensuremath x^0 = (\ensuremath x^0_1,\dots, \ensuremath x^0_m)  \equiv \ensuremath \mathsf{x}^0 \in \dom \ensuremath{{\gs}}$ be a constant random variable. Iterate

\vspace{-2.5ex}
\begin{equation}
\label{eq:algoPRCD2}
\begin{array}{l}
\text{for}\;n=0,1,\ldots\\
\left\lfloor
\begin{array}{l}
\text{for}\;i=1,\dots, m\\[0.7ex]
\left\lfloor
\begin{array}{l}
x^{\ensuremath n+1}_i =
x^\ensuremath n_i + \varepsilon^\ensuremath n_i \big(\ensuremath{\text{\sf prox}}_{\gamma_i \ensuremath{{\mathsf{h}}}_{i}}
\big(x^\ensuremath n_i - \gamma_{i} \nabla_i \ensuremath{{\mathsf{f}}} (\ensuremath x^{\ensuremath n})\big) - x^\ensuremath n_i\big).
\end{array}
\right.
\end{array}
\right.
\end{array}
\end{equation}
For every $\ensuremath n \in \ensuremath \mathbb{N}$,
we denote by $\ensuremath{\mathfrak{E}}_{\ensuremath n}$ the sigma-algebra generated by $\ensuremath \varepsilon^{0},\dots, \ensuremath \varepsilon^{\ensuremath n}$.
\end{algorithm}

In Algorithm~\ref{algoRCD}, the role of the random variable
$\ensuremath \varepsilon^\ensuremath n$ is to select, at iteration $\ensuremath n$,
the blocks to update in parallel (those indexed in
$\{i \in \{1,\dots, m\} \,\,\vert\, \varepsilon_i^\ensuremath n = 1\}$).
When all block-coordinates are simultaneously updated at each iteration, Algorithm~\ref{algoRCD} reduces to
the (deterministic) forward-backward algorithm, which converges
 only if the stepsizes are appropriately set.
More specifically, if $\nabla \ensuremath{{\mathsf{f}}}$ is  $L$-Lipschitz continuous,
then convergence is ensured if the stepsizes $\gamma_i$ are all equal and \emph{strictly less} than $2/L$ \cite{Com05,Dav16}. This fact is proved by using
the so called \emph{descent lemma}, i.e.,
\begin{equation}
\label{eq:20191012b}
(\forall\, \ensuremath \xx \in \ensuremath {\mathsf{H}})(\forall\, \ensuremath \vv\in \ensuremath {\mathsf{H}})\qquad \ensuremath{{\mathsf{f}}}(\ensuremath \xx + \ensuremath \vv) \leq \ensuremath{{\mathsf{f}}}(\ensuremath \xx) + \scalarp{\nabla \ensuremath{{\mathsf{f}}}(\ensuremath \xx), \ensuremath \vv}
+ \frac{L}{2} \norm{\ensuremath \vv}^2.
\end{equation}
Indeed, \eqref{eq:20191012b} is itself an assumption concerning the smoothness of $\ensuremath{{\mathsf{f}}}$,
since it is well-known to be equivalent to the Lipschitz continuity of the gradient of $\ensuremath{{\mathsf{f}}}$ \cite[Theorem~18.15]{book1}.
By contrast, when
the block-coordinates are updated one by one in a serial manner,
it is desirable to
allow moving along the block-coordinates with different
 stepsizes, depending on the
Lipschitz constants of the partial gradients of $\ensuremath{{\mathsf{f}}}$ across
the block-coordinates
\cite{Bec13,Nes12}. So, in this case
it is more appropriate to assume that a descent lemma holds on each block-coordinate subspace individually, that is,
\begin{equation}
\label{eq:S00}
(\forall\, i=1,\dots,m)(\forall\, \ensuremath \xx \in \ensuremath {\mathsf{H}})(\forall\,\ensuremath \mathsf{v}_i \in \ensuremath {\mathsf{H}}_i)\quad
\ensuremath{{\mathsf{f}}}(\ensuremath \xx + \ensuremath{\mathsf{J}}_i\ensuremath \vv_i)
\leq \ensuremath{{\mathsf{f}}}(\ensuremath \xx) +  \scalarp{\nabla_i \ensuremath{{\mathsf{f}}}(\ensuremath \xx), \ensuremath \vv_i}
+ \frac{L_i}{2} \norm{\ensuremath \mathsf{v}_i}^2,
\end{equation}
where $\ensuremath{\mathsf{J}}_i \ensuremath \vv_i = (0, \dots, 0, \ensuremath \vv_i,0, \dots, 0)$  ($\ensuremath \vv_i$ occurring at the $i$-th position),
for some positive constants $L_i$'s.
In the setting of Algorithm~\ref{algoRCD}, multiple block-coordinates may be updated in parallel at each iteration,
according to the random sampling $\varepsilon$. Therefore,
it is reasonable to assume that there exists $(\nu_i)_{1 \leq i \leq m} \in \ensuremath \mathbb{R}_{++}^\ensuremath \mathbb{N}$
so that one of the \emph{generalized smoothness conditions} below holds
\begin{enumerate}[{\rm S1}]
\item
\label{eq:S1}
$(\forall\, \ensuremath \xx,\ensuremath \vv \in \ensuremath {\mathsf{H}})\quad
\ensuremath{\mathsf E}[\ensuremath{{\mathsf{f}}}(\ensuremath \xx + \ensuremath \varepsilon \odot \ensuremath \vv)]
\leq \ensuremath{{\mathsf{f}}}(\ensuremath \xx) +  \ensuremath{\mathsf E}[\scalarp{\nabla \ensuremath{{\mathsf{f}}}(\ensuremath \xx), \ensuremath \varepsilon \odot \ensuremath \vv}]
+\displaystyle\frac{1}{2} \sum_{i=1}^m \ensuremath \mathsf{p}_i \nu_i \norm{\ensuremath \mathsf{v}_i}^2
$,
\item\label{eq:S2}
$(\forall\, \ensuremath \xx,\ensuremath \vv \in \ensuremath {\mathsf{H}})\quad
\ensuremath{{\mathsf{f}}}(\ensuremath \xx + \ensuremath \varepsilon \odot \ensuremath \vv)
\leq \ensuremath{{\mathsf{f}}}(\ensuremath \xx) +  \scalarp{\nabla \ensuremath{{\mathsf{f}}}(\ensuremath \xx), \ensuremath \varepsilon \odot \ensuremath \vv}
+ \displaystyle\frac{1}{2} \sum_{i=1}^m \nu_i \varepsilon_i \norm{\ensuremath \mathsf{v}_i}^2
\quad \ensuremath{\mathsf P} \text{ a.s.}$,
\end{enumerate}
where $\ensuremath \varepsilon \odot \ensuremath \vv = (\varepsilon_i \ensuremath \mathsf{v}_i)_{1 \leq i \leq m} \in \ensuremath {\mathsf{H}}$.
Conditions \ref{eq:S1} and \ref{eq:S2} can be interpreted
as descent lemmas on random block-coordinate subspaces,
 depending on the chosen random sampling of the block-coordinates.
They reduce to \eqref{eq:S00}, with $\nu_i=L_i$, if the sampling $\varepsilon$
selects only one block at a time almost surely (see Section~\ref{subsec:beta}).
We call $(\nu_i)_{1 \leq i \leq m}$ the
 \emph{smoothness parameters} of $\ensuremath{{\mathsf{f}}}$.
 Then, similarly to the deterministic case,
we will adopt the following stepsize rule
\begin{equation}
\label{eq:stepsizerule}
(\forall\, i \in \{1,\dots, m\})\qquad \gamma_i < \frac{2}{\nu_i}.
\end{equation}
Another smoothness condition suitable for Algorithm~\ref{algoRCD}, which was considered in \cite{Nec16}, is
\begin{enumerate}[{\rm S3}]
\item
\label{eq:S0}
$(\forall\, \ensuremath \xx,\ensuremath \vv \in \ensuremath {\mathsf{H}})\quad
\ensuremath{{\mathsf{f}}}(\ensuremath \xx + \ensuremath \vv)
\leq \ensuremath{{\mathsf{f}}}(\ensuremath \xx) +  \scalarp{\nabla \ensuremath{{\mathsf{f}}}(\ensuremath \xx), \ensuremath \vv}
+\displaystyle\frac{1}{2} \sum_{i=1}^m \nu_i \norm{\ensuremath \mathsf{v}_i}^2$.
\end{enumerate}
Note that, \ref{eq:S0} $\ \Rightarrow\ $ \ref{eq:S2} $\ \Rightarrow\ $ \ref{eq:S1},
which in turn implies the Lipschitz continuity of
the gradient of $\ensuremath{{\mathsf{f}}}$
(see Theorem~\ref{thm:stepsizes}\ref{thm:stepsizes_iv}). So, possibly with different values of the $\nu_i$'s,
the above conditions are all equivalent. The point is that in the parallel setting
(where multiple blocks are updated in parallel at each iteration),
 \ref{eq:S1} may be fulfilled with values of $\nu_i$ that are much smaller than
those related to the other two conditions, ultimately allowing to significantly increase the stepsizes
and hence speeding up the convergence. Moreover,
\ref{eq:S1} makes parallelization particularly effective
on problems with a sparse structure and superior to the serial strategy (which updates a single block per iteration).
See the discussion after Theorem~\ref{thm:20171207a}.
The critical role played by assumption \ref{eq:S1} in the analysis of
parallel randomized block-coordinate descent methods was pointed out
in \cite{Qu16b,Ric16,Ric16b,Tap18}. There, it was called \emph{expected separable overapproximation} (ESO) inequality.
 Condition
\ref{eq:S2} is new and serves to guarantee that Algorithm~\ref{algoRCD}
is almost surely descending (Proposition~\ref{p:20181219c}), which is a property
that is especially relevant when error bound conditions hold (see Section~\ref{sec:errorbounds}).
Note that in \cite{Ric16} the issue of monotonicity of the algorithm was addressed for each sampling separately without any general guidance.
Finally, we stress that, except for \cite{Com15,Com18}
(which study the convergence of the iterates only), in all previous works
the stepsizes $\gamma_i$'s are set equal to $1/\nu_i$.
This is an unnecessary limitation that we remove,
so to match the standard stepsize rule of the forward-backward algorithm \cite{Com05,Dav16}.

\begin{remark}
\label{rmk:20190126c}
For every $i=1,\dots, m$, the canonical embedding of $\ensuremath \mathsf{H}_i$ into $\ensuremath {\mathsf{H}}$ is
the operator $\ensuremath{\mathsf{J}}_i\colon \ensuremath \mathsf{H}_i \to \ensuremath {\mathsf{H}}$, $\ensuremath \mathsf{x} \mapsto (0, \dots, 0, \ensuremath \mathsf{x},0, \dots, 0)$,
where $\ensuremath \mathsf{x}$ occurs in the $i$-th position.
Then Algorithm~\ref{algoRCD} can be written as
\begin{equation*}
\ensuremath x^{\ensuremath n+1} = \ensuremath x^\ensuremath n + \sum_{i=1}^m \varepsilon_i^\ensuremath n \ensuremath{\mathsf{J}}_i
\big(\ensuremath{\text{\sf prox}}_{\gamma_i \ensuremath{{\mathsf{h}}}_{i}} \big(x^\ensuremath n_i - \gamma_{i} \nabla_i \ensuremath{{\mathsf{f}}} (\ensuremath x^{\ensuremath n})\big) - x^\ensuremath n_i\big).
\end{equation*}
\end{remark}

\subsection{Main contributions and comparison to previous work}

In the following we summarize the main contributions of this paper, where,
for the sake of brevity,  we set $\ensuremath{{\mathsf{F}}}=\ensuremath{{\mathsf{f}}}+\ensuremath{{\gs}}$.
We assume that \ref{eq:A1}--\ref{eq:A3} are satisfied and that \ref{eq:S1} is met.
Then, the following hold.
\begin{itemize}
\item
Algorithm~\ref{algoRCD} is descending in expectation and
$\ensuremath{\mathsf E}[\ensuremath{{\mathsf{F}}}(\ensuremath x^{\ensuremath n})] - \inf \ensuremath{{\mathsf{F}}}  \to 0$, even if the infimum is not attained.
If $\argmin \ensuremath{{\mathsf{F}}} \neq \varnothing$, then
$\ensuremath{\mathsf E}[\ensuremath{{\mathsf{F}}}(\ensuremath x^{\ensuremath n})] -  \inf \ensuremath{{\mathsf{F}}}= o(1/\ensuremath n)$. In addition,
a nonasymptotic bound for $\ensuremath{\mathsf E}[\ensuremath{{\mathsf{F}}}(\ensuremath x^{\ensuremath n})] -  \inf \ensuremath{{\mathsf{F}}}$ of order $O(1/\ensuremath n)$
holds.
Finally,
there exists a random variable $\ensuremath x_*$ with values in $\argmin \ensuremath{{\mathsf{F}}}$ such that
$\ensuremath x^\ensuremath n \rightharpoonup \ensuremath x_*$ $\ensuremath{\mathsf P}$-a.s. See Theorem~\ref{thm:20171207a}.
\item
If $\ensuremath{{\mathsf{F}}}$ satisfies an error bound condition of Luo-Tseng type
(see condition~\ref{eq:EB1}),
then
 the  iterates
 as well as the corresponding
function values generated by Algorithm~\ref{algoRCD}, converge linearly in expectation.
See Theorem~\ref{thm:20181214a} and
Theorem~\ref{thm20200218a}.
\end{itemize}

Our results advance the state-of-the-art in the study of random block-coordinate
descent methods under several aspects. We comment on this below.
1) While convergence of the function values has been intensively studied in the related literature
(see e.g., \cite{Kar18,LuX15,Nec16,Nes12,Qu15,Ric16,Ric16b,Tap18}), surprisingly, in a convex setting,
convergence of the iterates has been investigated only recently in \cite{Com15}, but with stepsizes set
according to the global Lipschitz constant of $\nabla \ensuremath{{\mathsf{f}}}$.
See also \cite{Fer19} which addresses the convergence of the iterates in the framework
of primal-dual algorithms with a serial and uniform block sampling.
We improve the existing results, since we show \emph{convergence of the iterates}
for Algorithm~\ref{algoRCD} in an \emph{infinite dimensional setting}
even when the stepsizes are chosen according to the condition \ref{eq:S1},
which can incorporate the block Lipschitz constants of the gradient of $\ensuremath{{\mathsf{f}}}$ and is at the basis of the effectiveness of the
parallel block-coordinatewise approach.
2) The worst case \emph{asymptotic rate $o(1/\ensuremath n)$} for the mean of the function values is new in the setting of stochastic algorithms.
3) Our analysis spotlights an abstract convergence principle for stochastic descent algorithms (Theorem~\ref{p:SMFejer}) which is essentially a special form
 of the stochastic quasi-Fej\'er monotonicity property, involving also the values of the objective functions.
This principle, previously investigated in
a deterministic setting in \cite{Sal17}, allows to prove in
a unified way both the almost sure convergence of the iterates and rates of convergence for the mean of the function values.
4) As a by-product of the above analysis we single out an inequality (Proposition~\ref{p:20190610a})
which is pivotal for studying
\emph{the convergence under error
bound conditions}, improving the results and simplifying the analysis in \cite{Nec16}.
5) We allow for  {\em parallel and arbitrary sampling
of the blocks} in a composite setting.
The benefit of
such sampling in terms of convergence rate have been first investigated in \cite{Ric16b}
for a strongly convex and smooth objective function.
In \cite{Qu16} a composite objective optimization problem was analyzed but for a slightly different algorithm. The rest of the studies deal either with parallel uniform sampling of the blocks \cite{Ric16}, or with the case where a single block is updated at each iteration \cite{LuX15,Nes12}.
6) We also allow for {\em stepsizes larger than those considered in literature}
\cite{Kar18,LuX15,Nec16,Nes12,Qu15,Ric16,Ric16b,Tap18}, since we can let the stepsizes be arbitrarily close to
$2/\nu_i$, matching the standard rule for the forward-backward algorithm.

The rest of the paper is organized as follows.
In Section~\ref{sec:notation} we give notation and basic facts.
Section~\ref{sec:stepsizes} shows how to determine the smoothness parameters
$\nu_i$ when $\ensuremath{{\mathsf{f}}}$ features a partially separable structure.
In Section~\ref{sec:conana} we carry out the convergence analysis and give the related theorems.
Finally, Section~\ref{sec:app} shows three applications
and Section~\ref{sec:experiments} provides some numerical experiments.

\section{Notation and background}
\label{sec:notation}

\paragraph{Notation.}
We define $\ensuremath \mathbb{R}_+ = \left[0,+\infty\right[$, $\ensuremath \mathbb{R}_{++} = \left]0,+\infty\right[$,
for every integer $s \geq 1$, $[s] = \{1, \dots, s\}$, and for every $a \in \ensuremath \mathbb{R}^s$,
 $\mathrm{spt}(a) = \{i \in [s] \,\vert\, a_i\neq 0\}$.
 Scalar products and norms in Hilbert spaces are denoted by $\scalarp{\cdot,\cdot}$
 and $\norm{\cdot}$ respectively.
If $\ensuremath{\mathsf{U}}\colon \ensuremath \mathsf{H} \to \ensuremath \mathsf{G}$ is a bounded linear operator between
real Hilbert spaces,
$\ensuremath{\mathsf{U}}^\top \colon \ensuremath \mathsf{G} \to \ensuremath \mathsf{H}$ is its transpose operator, that is, the one satisfying $\scalarp{\ensuremath{\mathsf{U}} \ensuremath \mathsf{x}, \ensuremath \mathsf{y}} = \scalarp{\ensuremath \mathsf{x}, \ensuremath{\mathsf{U}}^\top \ensuremath \mathsf{y}}$,
for every $(\ensuremath \mathsf{x},\ensuremath \mathsf{y}) \in \ensuremath \mathsf{H}\times \ensuremath \mathsf{G}$.
Let $(\ensuremath \mathsf{H}_i)_{1 \leq i \leq m}$ be $m$ separable real Hilbert spaces and let
$\ensuremath {\mathsf{H}}=\bigoplus_{i=1}^m \ensuremath \mathsf{H}_i$ be their direct sum.
For every $\ensuremath \vv \in \ensuremath {\mathsf{H}} $ and $\ensuremath \epsilon \in \{0,1\}^m$
we set $\ensuremath \epsilon \odot\ensuremath \vv = (\epsilon_i \ensuremath \mathsf{v}_i)_{1 \leq i \leq m} \in \ensuremath {\mathsf{H}}$.
We will consider random variables with underlying probability space $(\Omega, \mathfrak{A}, \ensuremath{\mathsf P})$ taking values in $\ensuremath \mathsf{H}_i$ or $\ensuremath {\mathsf{H}}$.  We use
the default font for random variables and sans serif font
for their realizations. The expected value operator is denoted by $\ensuremath{\mathsf E}$.
 A copy of a random variable is  random variable having the same distribution of the given one.
Let $(w_i)_{1 \leq i \leq m} \in \ensuremath \mathbb{R}^m_{++}$.
The direct sum operator
$\ensuremath{{\WW}} = \bigoplus_{i=1}^m w_i \ensuremath{\mathsf{Id}}_i$,
where $\ensuremath{\mathsf{Id}}_i$
is the identity operator on $\ensuremath \mathsf{H}_i$, is the positive bounded linear operator on $\ensuremath {\mathsf{H}}$ acting as
$\ensuremath \xx=(\ensuremath \mathsf{x}_i)_{1 \leq i \leq m} \mapsto (w_i \ensuremath \mathsf{x}_i)_{1 \leq i \leq m}$.
$\ensuremath{{\WW}}$ defines an equivalent inner product on $\ensuremath {\mathsf{H}}$
\begin{equation*}
(\forall\, \ensuremath \xx \in \ensuremath {\mathsf{H}})(\forall\, \ensuremath \yy \in \ensuremath {\mathsf{H}})\qquad\scalarp{\ensuremath \xx, \ensuremath \yy}_{\ensuremath{{\WW}}} = \scalarp{\ensuremath{{\WW}} \ensuremath \xx, \ensuremath \yy}
= \sum_{i=1}^m w_i \scalarp{\ensuremath \mathsf{x}_i, \ensuremath \mathsf{y}_i},
\end{equation*}
which gives the norm $\norm{\ensuremath \xx}_{\ensuremath{{\WW}}}^2
= \sum_{i=1}^m w_i \norm{\ensuremath \mathsf{x}_i}^2$.
If $\ensuremath {\mathsf{S}} \subset \ensuremath {\mathsf{H}}$ and $\ensuremath \xx \in \ensuremath {\mathsf{H}}$, we set
 $\mathrm{dist}_{\ensuremath{{\WW}}}(\ensuremath \xx,\ensuremath {\mathsf{S}})
= \inf_{\bm{\ensuremath \mathsf{z}} \in \ensuremath {\mathsf{S}}} \norm{\ensuremath \xx - \bm{\ensuremath \mathsf{z}}}_{\ensuremath{{\WW}}}$.
Let $\ensuremath{{\varphi}}\colon \ensuremath {\mathsf{H}} \to \left]-\infty,+\infty\right]$ be proper, convex, and lower semicontinuous. The domain of $\ensuremath{{\varphi}}$ is $\dom \ensuremath{{\varphi}}
= \{\ensuremath \xx \in \ensuremath \mathsf{H} \,\vert\, \ensuremath{{\varphi}}(\ensuremath \xx)<+\infty\}$ and the set
of minimizers of $\ensuremath{{\varphi}}$ is $\argmin \ensuremath{{\varphi}}
= \{\ensuremath \xx \in \ensuremath {\mathsf{H}} \,\vert\, \ensuremath{{\varphi}}(\ensuremath \xx) = \inf \ensuremath{{\varphi}}\}$.
The subdifferential of $\ensuremath{{\varphi}}$ in the metric $\scalarp{\cdot,\cdot}_{\ensuremath{{\WW}}}$
is the multivalued operator
\begin{equation*}
\partial^{\ensuremath{{\WW}}} \ensuremath{{\varphi}} \colon \ensuremath {\mathsf{H}} \to 2^{\ensuremath {\mathsf{H}}},\ \
\ensuremath \xx \mapsto \partial^{\ensuremath{{\WW}}} \ensuremath{{\varphi}}(\ensuremath \xx)= \big\{\ensuremath \uu \in \ensuremath {\mathsf{H}} \,\vert\,
(\forall\, \ensuremath \yy \in \ensuremath {\mathsf{H}})\ \ensuremath{{\varphi}}(\ensuremath \yy) \geq \ensuremath{{\varphi}}(\ensuremath \xx)
+ \scalarp{\ensuremath \uu, \ensuremath \yy - \ensuremath \xx}_{\ensuremath{{\WW}}}\big\}.
\end{equation*}
In case $\ensuremath{{\WW}} = \ensuremath{\mathsf{Id}}$, it is simply denoted by $\partial \ensuremath{{\varphi}}$.
Clearly $\partial^{\ensuremath{{\WW}}} \ensuremath{{\varphi}} = \ensuremath{{\WW}}^{-1} \partial \ensuremath{{\varphi}}$.
If the function $\ensuremath{{\varphi}}\colon \ensuremath {\mathsf{H}}\to \ensuremath \mathbb{R}$ is differentiable,
then, for every $\ensuremath \xx \in \ensuremath {\mathsf{H}}$, $\partial^{\ensuremath{{\WW}}} \ensuremath{{\varphi}}(\ensuremath \xx) = \{\nabla^{\ensuremath{{\WW}}} \ensuremath{{\varphi}}(\ensuremath \xx)\}$ and for all $\ensuremath \vv \in \ensuremath {\mathsf{H}}$, $\scalarp{\nabla^{\ensuremath{{\WW}}} \ensuremath{{\varphi}}(\ensuremath \xx), \ensuremath \vv}_{\ensuremath{{\WW}}} = \scalarp{\nabla \varphi(\ensuremath \xx), \ensuremath \vv}$.
The proximity operator of $\ensuremath{{\varphi}}$ in the metric $\scalarp{\cdot,\cdot}_{\ensuremath{{\WW}}}$  is defined as
\begin{equation*}
\ensuremath{\text{\sf prox}}_{\ensuremath{{\varphi}}}^{\ensuremath{{\WW}}} \colon \ensuremath {\mathsf{H}}\to \ensuremath {\mathsf{H}},\quad
\ensuremath{\text{\sf prox}}_{\ensuremath{{\varphi}}}^{\ensuremath{{\WW}}}(\ensuremath \xx) = \argmin_{\ensuremath \zz \in \ensuremath {\mathsf{H}}} \ensuremath{{\varphi}}(\ensuremath \zz) + \frac 1 2 \norm{\ensuremath \xx - \ensuremath \zz}^2_{\ensuremath{{\WW}}}.
\end{equation*}
Referring to the functions in \eqref{eq:20190110b},
we denote by $\mu_{\ensuremath \mathsf{\Gamma}^{-1}}$ and $\sigma_{\ensuremath \mathsf{\Gamma}^{-1}}$ the moduli of strong convexity
of $\ensuremath{{\mathsf{f}}}$ and $\ensuremath{{\gs}}$ respectively,
 in the norm $\norm{\cdot}_{\ensuremath{{\Gammas}}^{-1}}$, where
 $\ensuremath{{\Gammas}} = \bigoplus_{i=1}^m \gamma_i \ensuremath{\mathsf{Id}}_i $
 and the $\gamma_i$'s
 are the stepsizes occurring in Algorithm~\ref{algoRCD}.
 This means that $\mu_{\ensuremath \mathsf{\Gamma}^{-1}}, \sigma_{\ensuremath \mathsf{\Gamma}^{-1}} \in \ensuremath \mathbb{R}_+$ and that,
 for every $\ensuremath \xx,\ensuremath \yy \in \ensuremath {\mathsf{H}}$,
\begin{align}
\label{eq:20181112b}\ensuremath{{\mathsf{f}}}(\ensuremath \yy) &\geq \ensuremath{{\mathsf{f}}}(\ensuremath \xx) + \scalarp{\nabla \ensuremath{{\mathsf{f}}}(\ensuremath \xx), \ensuremath \yy - \ensuremath \xx}
+ \frac{\mu_{\ensuremath \mathsf{\Gamma}^{-1}}}{2} \sum_{i=1}^m \frac{1}{\gamma_i}\norm{\ensuremath \mathsf{y}_i - \ensuremath \mathsf{x}_i}^2,\\
\label{eq:20181112b1}(\forall\, \ensuremath \vv \in \partial \ensuremath{{\gs}}(\ensuremath \xx))\qquad
\ensuremath{{\gs}}(\ensuremath \yy) &\geq \ensuremath{{\gs}}(\ensuremath \xx) + \scalarp{\ensuremath \vv, \ensuremath \yy - \ensuremath \xx}
+ \frac{\sigma_{\ensuremath \mathsf{\Gamma}^{-1}}}{2} \sum_{i=1}^m \frac{1}{\gamma_i}\norm{\ensuremath \mathsf{y}_i - \ensuremath \mathsf{x}_i}^2.
\end{align}
Note that, since $\ensuremath{{\gs}}$ is separable, by taking $\ensuremath \yy= \ensuremath \xx + \ensuremath{\mathsf{J}}_i(\ensuremath \mathsf{y}_i - \ensuremath \mathsf{x}_i)$ in \eqref{eq:20181112b1}, we have
\begin{equation}
\label{eq:20181112a}
(\forall\, i \in [m])\quad
\ensuremath{{\mathsf{h}}}_i(\ensuremath \mathsf{y}_i) \geq \ensuremath{{\mathsf{h}}}_i(\ensuremath \mathsf{x}_i) + \scalarp{\ensuremath \mathsf{v}_i, \ensuremath \mathsf{y}_i - \ensuremath \mathsf{x}_i} + \frac{\sigma_{\ensuremath \mathsf{\Gamma}^{-1}}}{2} \frac{1}{\gamma_i} \norm{\ensuremath \mathsf{y}_i - \ensuremath \mathsf{x}_i}^2.
\end{equation}

\begin{remark}
If \ref{eq:S1} is satisfied, the $\gamma_i$'s are chosen
as in \eqref{eq:stepsizerule}, and $\delta = \max_{1 \leq i \leq m} \gamma_i \nu_i$ (according to the convergence theorems), then we have
\begin{equation}
\label{eq:20190126e}
\mu_{\ensuremath \mathsf{\Gamma}^{-1}} \leq \min_{1 \leq i \leq m}\gamma_i \nu_i \leq \delta <2.
\end{equation}
Indeed,
let $\ensuremath \xx\in \ensuremath {\mathsf{H}}$,  $i \in [m]$ and $\ensuremath \mathsf{v}_i \in \ensuremath \mathsf{H}_i$, $\ensuremath \mathsf{v}_i \neq 0$.
It follows from \eqref{eq:20181112b} with
$\ensuremath \yy = \ensuremath \xx + \ensuremath \varepsilon\odot \ensuremath{\mathsf{J}}_i \ensuremath \mathsf{v}_i$ (where $\ensuremath{\mathsf{J}}_i$ is defined in
Remark~\ref{rmk:20190126c}) and  \ref{eq:S1} that
\begin{align*}
\ensuremath{{\mathsf{f}}}(\ensuremath \xx) + \ensuremath{\mathsf E}[\scalarp{\nabla \ensuremath{{\mathsf{f}}}(\ensuremath \xx), \ensuremath \varepsilon \odot \ensuremath{\mathsf{J}}_i \ensuremath \mathsf{v}_i}] + \frac{\mu_{\ensuremath \mathsf{\Gamma}^{-1}}}{2} \frac{\ensuremath \mathsf{p}_i}{\gamma_i}\norm{\ensuremath \mathsf{v}_i}^2
&\leq \ensuremath{\mathsf E}[\ensuremath{{\mathsf{f}}}(\ensuremath \xx + \ensuremath \varepsilon\odot \ensuremath{\mathsf{J}}_i \ensuremath \mathsf{v}_i)] \\
&\leq \ensuremath{{\mathsf{f}}}(\ensuremath \xx) +  \ensuremath{\mathsf E}[\scalarp{\nabla \ensuremath{{\mathsf{f}}}(\ensuremath \xx), \ensuremath \varepsilon \odot \ensuremath{\mathsf{J}}_i \ensuremath \mathsf{v}_i}]
+\displaystyle\frac{1}{2}  \ensuremath \mathsf{p}_i \nu_i \norm{\ensuremath \mathsf{v}_i}^2.
\end{align*}
Thus, \eqref{eq:20190126e} follows.
\end{remark}

\begin{fact}[{\cite[Example~5.1.5]{book2}}]
\label{f:20190112d}
Let $\zeta_1$ and $\zeta_2$ be independent random variables with values in
the measurable spaces $\mathcal{Z}_1$ and $\mathcal{Z}_2$ respectively.
Let $\varphi\colon \mathcal{Z}_1\times \mathcal{Z}_2 \to \ensuremath \mathbb{R}$ be measurable
and suppose that $\ensuremath{\mathsf E}[\abs{\varphi(\zeta_1,\zeta_2)}]<+\infty$.
Then $\ensuremath{\mathsf E}[\varphi(\zeta_1,\zeta_2) \,\vert\, \zeta_1] = \psi(\zeta_1)$,
where for all $z_1 \in \mathcal{Z}_1$, $\psi(z_1) = \ensuremath{\mathsf E}[\varphi(z_1, \zeta_2)]$.
\end{fact}

\begin{fact}
\label{lem:20181024a}
Let $\ensuremath \varepsilon$ be a random variable with values in  $\{0,1\}^m$ and, for all $i \in [m]$, $\ensuremath \mathsf{p}_i = \ensuremath{\mathsf P}(\varepsilon_i = 1)$.
Then $\ensuremath{\mathsf E}[\varepsilon_i] = \ensuremath \mathsf{p}_i$ and, for every
 $\ensuremath \vv = (\ensuremath \mathsf{v}_i)_{1 \leq i\leq m} \in \ensuremath \mathbb{R}^m$,  $\ensuremath{\mathsf E}[\scalarp{\ensuremath \varepsilon, \ensuremath \vv} ]
 = \sum_{i=1}^m \ensuremath \mathsf{p}_i \ensuremath \mathsf{v}_i$.
\end{fact}

\begin{fact}[\cite{book3}]
\label{lem:monosequence}
Let $(a_\ensuremath n)_{\ensuremath n \in \ensuremath \mathbb{N}}$ be a decreasing sequence in $\ensuremath \mathbb{R}_+$.
If $\sum_{\ensuremath n=0}^{+\infty} a_\ensuremath n <+\infty$, then,
for every $\ensuremath n \in \ensuremath \mathbb{N}$, $a_\ensuremath n \leq (1/(\ensuremath n+1)) \sum_{\ensuremath n=0}^{+\infty} a_\ensuremath n$
and $a_\ensuremath n = o\big( 1/(\ensuremath n+1)\big)$.
\end{fact}

\section{Determining the smoothness parameters}
\label{sec:stepsizes}
In this section we provide few scenarios for which the relaxed
smoothness conditions \ref{eq:S1} and \ref{eq:S2} can be fully exploited,
attaining tight values for the $\nu_i$'s.
This ultimately allows to take larger stepsizes and improves rates of convergence.
In \cite{Qu16b,Tap18} an extensive analysis of cases in which \ref{eq:S1}
is satisfied is presented.

\subsection{General estimates.}
We consider the following setting.
\begin{enumerate}[{\rm H4}]
\item\label{eq:B0}
The function $\ensuremath{{\mathsf{f}}}\colon \ensuremath {\mathsf{H}} \to \ensuremath \mathbb{R}$ is such that
\begin{equation}
\label{eq:20190605b}
(\forall\, \ensuremath \xx \in \ensuremath {\mathsf{H}})\qquad
\ensuremath{{\mathsf{f}}}(\ensuremath \xx) = \sum_{k=1}^p \ensuremath{{\mathsf{g}}}_k \bigg(\sum_{i=1}^m \ensuremath{\mathsf{U}}_{k,i} \ensuremath \mathsf{x}_i\bigg),
\end{equation}
where, for every $k =1,\dots, p$,
 $\ensuremath{{\mathsf{g}}}_k \colon \ensuremath \mathsf{G}_k \to \ensuremath \mathbb{R}$ is a convex differentiable function
 defined on a real Hilbert space $\ensuremath \mathsf{G}_k$
 and, for every $i \in [m]$, $\ensuremath{\mathsf{U}}_{k,i}\colon \ensuremath \mathsf{H}_i \to \ensuremath \mathsf{G}_k$
 is a bounded linear operator. Moreover,
 $\bigcup_{k = 1}^p I_k \neq \varnothing$, where, for all
 $k =1,\dots, p$, $I_k = \big\{i \in [m] \,\vert\, \ensuremath{\mathsf{U}}_{k,i} \neq 0\big\}$,
and $\eta = \max_{1 \leq k \leq p} \mathrm{card}(I_k)$.
\end{enumerate}
We will also consider one of the following conditions.
 \begin{enumerate}[{\rm L1}]
\item\label{eq:B1} For every $i=1, \dots, m$ there exists $L_i>0$ such that, for every
$\ensuremath \xx \in \ensuremath {\mathsf{H}}$,  the function
$\nabla_i \ensuremath{{\mathsf{f}}} (\ensuremath \mathsf{x}_1, \dots, \ensuremath \mathsf{x}_{i-1}, \cdot, \ensuremath \mathsf{x}_{i+1}, \dots, \ensuremath \mathsf{x}_m) \colon \ensuremath \mathsf{H}_i \to \ensuremath \mathsf{H}_i$
 is $L_i$-Lipschitz continuous.
 \item\label{eq:B2}
For every $k=1,\dots, p$,
 $\nabla \ensuremath{{\mathsf{g}}}_k \colon \ensuremath \mathsf{G}_k \to \ensuremath \mathsf{G}_k$ is $L^{(k)}$-Lipschitz continuous and for
 every $i,j \in [m]$, $i\neq j$, the ranges of   $\ensuremath{\mathsf{U}}_{k,i}$ and $\ensuremath{\mathsf{U}}_{k,j}$
 are orthogonal.
\end{enumerate}

Assumption \ref{eq:B0} concerns the \emph{partial separability} of the function $\ensuremath{{\mathsf{f}}}$.
Depending on the number of the nonzero operators $\ensuremath{\mathsf{U}}_{k,i}$,
$\ensuremath{{\mathsf{g}}}_k$ might depend only on few block-variables $\ensuremath \mathsf{x}_i$'s:
if $\eta = 1$, $\ensuremath{{\mathsf{f}}}$ is fully separable, whereas if $\eta = m$, $\ensuremath{{\mathsf{f}}}$ is not separable.
Note that \ref{eq:B1} is equivalent to \eqref{eq:S00} and,
since $\ensuremath{{\mathsf{f}}}$ is convex, implies the global Lipschitz continuity of the
gradient of $\ensuremath{{\mathsf{f}}}$ (Corollary~\ref{p:20181014a}). So either \ref{eq:B1} or \ref{eq:B2}
implies the global Lipschitz smoothness of $\ensuremath{{\mathsf{f}}}$.
However, considering the constants $L_i$'s or $L^{(k)}$'s leads in general to
a finer analysis of the smoothness properties of $\ensuremath{{\mathsf{f}}}$,
eventually determining parameters $\nu_i$ that are smaller than the global Lipschitz constant of $\nabla \ensuremath{{\mathsf{f}}}$.
Instances of problem~\eqref{eq:20190110b} where $\ensuremath{{\mathsf{f}}}$
has the structure shown in \ref{eq:B0},
occur very often in applications. In particular,
a prominent example is that of the \emph{Lasso problem} which
will be discussed in Section~\ref{subsec:lasso}.
The following theorem, which is proved in Appendix~\ref{sec:appB},
relates the smoothness parameters $(\nu_i)_{1 \leq i \leq m}$  to the block Lipschitz constants
of the partial gradients of $\ensuremath{{\mathsf{f}}}$ and to the Lipschitz constants of the gradients
of its components $\ensuremath{{\mathsf{g}}}_k$'s in
\eqref{eq:20190605b}, as well as to the distribution of the random variable $\ensuremath \varepsilon$.

\begin{theorem}
\label{thm:stepsizes}
Assume  \ref{eq:A3} and \ref{eq:B0} and let
 $(\nu_i)_{1 \leq i \leq m} \in \ensuremath \mathbb{R}_{++}^\ensuremath \mathbb{N}$.
 Then the following hold.
 \begin{enumerate}[{\rm (i)}]
 \item\label{thm:stepsizes_i} \ref{eq:B1}$\ \Rightarrow\  $ \ref{eq:S1} provided that
\begin{equation*}
(\forall\, i \in [m])\quad
\nu_i \geq \beta_{1,i} L_i,\quad\text{where}\quad
\beta_{1,i} := \ensuremath{\mathsf E} \Big[
\max_{1 \leq k \leq p} \Big(\sum_{j\in I_k} \varepsilon_j \Big)\,\Big\vert\, \varepsilon_i=1
\Big].
\end{equation*}
\vspace{-4ex}
\item \label{thm:stepsizes_ii}
\ref{eq:B1}$\ \Rightarrow\  $ \ref{eq:S2} provided that
\begin{equation*}
(\forall\, i \in [m])\quad \nu_i \geq \beta_2 L_i,\quad\text{where}\quad
\beta_2 :=\esssup \Big(\max_{1 \leq k \leq p}\Big(\sum_{j \in I_k} \varepsilon_j \Big) \Big).
\end{equation*}
\vspace{-4ex}
\item\label{thm:stepsizes_iii}
\ref{eq:B2}$\ \Rightarrow\  $ \ref{eq:S0}
provided that
\begin{equation}
(\forall\, i \in [m])\quad
\nu_i \geq \tilde{L}_i := \Big\lVert \sum_{k=1}^p L^{(k)} \ensuremath{\mathsf{U}}_{k,i}^\top \ensuremath{\mathsf{U}}_{k,i}\Big\rVert.
\end{equation}
\vspace{-3ex}
\item\label{thm:stepsizes_iv}
\ref{eq:S1}$\ \Rightarrow\  $\ref{eq:B1} with, for all $i \in [m]$, $L_i = \nu_i$.
In particular, \ref{eq:S1} implies that $\ensuremath{{\mathsf{f}}}$ is Lipschitz smooth.
 \end{enumerate}
\end{theorem}

\begin{remark}\
\label{rmk:20190128e}
\begin{enumerate}[{\rm (i)}]
\item\label{rmk:20190128e_i} Suppose that in \ref{eq:B0},
for all $k \in [p]$,
$\ensuremath \mathsf{G}_k = \ensuremath {\mathsf{H}}$, $\ensuremath{{\mathsf{g}}}_k$ is
$L^{(k)}$-Lipschitz smooth, and, for all $i \in I_k$,
$\ensuremath{\mathsf{U}}_{k,i} = \ensuremath{\mathsf{J}}_i$, the canonical embedding of $\ensuremath \mathsf{H}_i$ into $\ensuremath {\mathsf{H}}$ (see Remark~\ref{rmk:20190126c}).
Then, \ref{eq:B2} holds and, for every $i \in [m]$,
$\tilde{L}_i = \sum_{k\,\vert i \in I_k} L^{(k)}$.
Hence, in view of Theorem~\ref{thm:stepsizes}\ref{thm:stepsizes_iii},
\ref{eq:S0}
is met with $\nu_i = \tilde{L}_i$.
This setting was studied in \cite{Nec16}.
\item\label{rmk:20190128e_ii} If $\nabla \ensuremath{{\mathsf{f}}}$ is $L$-Lipschitz continuous,
then \ref{eq:S0}
is satisfied with,
for every $i \in [m]$, $\nu_i = L$.
Therefore, we cover the analysis of the random
block-coordinate forward-backward algorithm given in
 \cite{Com15,Com18} which set the stepsizes as $\gamma_i < 2/L$.
\item\label{rmk:20190128e_iii} Let, for every $k \in [p]$, $\ensuremath{{\mathsf{f}}}_k(\ensuremath \mathsf{x}) =  \ensuremath{{\mathsf{g}}}_k \big(\sum_{i=1}^m \ensuremath{\mathsf{U}}_{k,i} \ensuremath \mathsf{x}_i\big)$.
If, for every $k \in [p]$, \ref{eq:S1} (resp. \ref{eq:S2}) holds for $\ensuremath{{\mathsf{f}}}_k$ with $(\nu_i^{(k)})_{1 \leq i \leq m}$ ,
then \ref{eq:S1} (resp. \ref{eq:S2}) holds for $\ensuremath{{\mathsf{f}}}$ with $\nu_i = \sum_{k=1}^p \nu_i^{(k)}$.
\item\label{rmk:20190128e_iv} Using similar ideas as in the proof of \cite[Theorem~12]{Ric16} we show in Appendix~\ref{sec:appB} that
item \ref{thm:stepsizes_i} in Theorem~\ref{thm:stepsizes} remains true with
\begin{equation}
\label{eq:20200316a}
\beta_{1,i} := \sum_{t=1}^\eta t \max_{\substack{1 \leq k \leq p \\[0.3ex] i \in I_k}}
\ensuremath{\mathsf P}\Big(\sum_{j \in I_k} \varepsilon_j = t \,\big\vert\, \varepsilon_i=1 \Big).
\end{equation}
\end{enumerate}
\end{remark}
\begin{remark}
\label{rmk:20181219ee}
Referring to Theorem~\ref{thm:stepsizes},
for all $i \in [m]$, we have $1 \leq \beta_{1,i} \leq \beta_2
\leq \min\{\eta,\tau_{\max}\}$,
where
\begin{align*}
\tau_{\max} :=\esssup \Big(\sum_{i =1}^m \varepsilon_i \Big)
&=\min\bigg\{ \tau^\prime \in \ensuremath \mathbb{N} \,\Big\vert\, \ensuremath{\mathsf P}\bigg( \sum_{i =1}^m
\varepsilon_i \leq \tau^\prime \bigg)=1 \bigg\}
\end{align*}
is the maximum number of blocks processed in parallel.
Indeed, since
$\ensuremath{\mathsf P}(\ensuremath \varepsilon \equiv 0) = 0$
we have $\ensuremath{\mathsf P}(\max_{1 \leq i \leq  m} \varepsilon_i \geq 1) = 1$. Moreover,
since $\max_{1 \leq k \leq p} \big(\sum_{i \in I_k} \varepsilon_i \big)
\geq \max_{1 \leq i \leq m} \varepsilon_i$, we have $1 \leq \beta_{1,i}$.
The inequality $\beta_{1,i} \leq \beta_2$ is immediate, while the last one derives from the following
\begin{equation*}
(\forall\, k \in [p])
\quad\sum_{i \in I_k} \varepsilon_i \leq \min\Big\{\mathrm{card}(I_k), \sum_{i=1}^m \varepsilon_i \Big\}
\leq \min\Big\{\eta, \sum_{i=1}^m \varepsilon_i \Big\} \leq \min\{\eta, \tau_{\max}\}.
\end{equation*}
\end{remark}

\subsection{The smoothness parameters for some special block samplings.}
\label{subsec:beta}
Here
we show how to compute (or estimate)
the constants $(\beta_{1,i})_{1 \leq i \leq m}$ and $\beta_2$ in
Theorem~\ref{thm:stepsizes} and Remark~\ref{rmk:20190128e}\ref{rmk:20190128e_iv},
and the related $(\nu_i)_{1 \leq i \leq m}$, in some relevant scenarios, when \ref{eq:B0} and \ref{eq:B1}  are satisfied.
\setlist[description]{font=\normalfont\itshape}
\begin{description}[leftmargin=0cm]
\item[Arbitrary parallel sampling.]
It follows from Theorem~\ref{thm:stepsizes}\ref{thm:stepsizes_ii} and Remark~\ref{rmk:20181219ee}
that for an arbitrary (possibly nonuniform) block sampling $\ensuremath \varepsilon$, \ref{eq:S2} is satisfied provided that
$\nu_i = \min\{\eta,\tau_{\max}\} L_i$, for every $i \in [m]$.
Additionally, if we denote by $L^{(k)}_i$ the blockwise Lipschitz constants of the gradient of the function
 $\ensuremath \mathsf{x} \mapsto \ensuremath{{\mathsf{g}}}_k \big(\sum_{i=1}^m \ensuremath{\mathsf{U}}_{k,i} \ensuremath \mathsf{x}_i\big)$, then we derive from
 Remark~\ref{rmk:20190128e}\ref{rmk:20190128e_iii}
 and the above discussion that \ref{eq:S2} holds
 with\footnote{If $i \notin I_k$, then $L^{(k)}_i = 0$.}
 $\nu_i = \sum_{k \vert i \in I_k} \min\{\mathrm{card}(I_k), \tau_{\max}\} L^{(k)}_i$.
However, the above estimates are rather conservative
and can be improved for special choices of the block sampling as we will show below.
We refer to \cite{Qu16b,Ric16b} for further results on nonuniform samplings.
\item[Serial sampling or full separability.]
Suppose that $\tau_{\max}=1$  or $\eta=1$.
Then, Remark~\ref{rmk:20181219ee} yields
$\beta_{1,i}=\beta_2= \min\{\eta,\tau_{\max}\}=1$. Moreover, recalling
Theorem~\ref{thm:stepsizes}\ref{thm:stepsizes_ii}-\ref{thm:stepsizes_iv}, this also shows that
 \ref{eq:S1}, \ref{eq:S2},  and \ref{eq:B1} are indeed equivalent with
the same smoothness parameters $\nu_i = L_i$.
So, conditions \ref{eq:S1} or \ref{eq:S2} find their justification
 only in the parallel case ($\tau_{\max}>1$) and when $\ensuremath{{\mathsf{f}}}$ is not
fully separable ($\eta>1$).
\item[Fully Parallel.] If $\ensuremath{\mathsf P}(\sum_{i=1}^m \varepsilon_i = m) = 1$, then
 for every $i \in [m]$
$\ensuremath \mathsf{p}_i =1$. This yields a fully parallel (deterministic) algorithm.
Moreover, since $\ensuremath{\mathsf P}\big(\ensuremath \varepsilon = (1,\dots, 1)\big)=1$, we have
$\beta_{1,i} = \eta =\beta_2$ and hence \ref{eq:S2} holds with $\nu_i = \eta L_i$. Actually,
also \ref{eq:S0} holds with $\nu_i = \eta L_i$ (see Corollary~\ref{p:20181014a}\ref{p:20181014a_3}).
\item[Uniform samplings.]
Suppose that $m>1$.
The sampling is \emph{uniform} if $\ensuremath \mathsf{p}_i = \ensuremath \mathsf{p}_j$, with $i \neq j$.
In this case if we denote by $\bar{\tau}$ the average number of block updates per iteration, we have
$\bar{\tau}= \ensuremath{\mathsf E} [\sum_{i=1}^m \varepsilon_i] = \sum_{i=1}^m \ensuremath \mathsf{p}_i$ and hence $\ensuremath \mathsf{p}_i = \bar{\tau}/m$,
for every $i \in [m]$.
In \cite{Ric16} several types of uniform samplings are studied.
In the following we single out two of them.
The sampling is said to be \emph{doubly uniform} if any
 two sets of blocks with the same number of blocks have the same probability to be chosen.
 In formula, this means that for every $J_1, J_2 \subset [m]$ such that $\mathrm{card}(J_1) = \mathrm{card}(J_2)$,  $\ensuremath{\mathsf P}(\cap_{i \in J_1} \{\varepsilon_i = 1\}) = \ensuremath{\mathsf P}(\cap_{i \in J_2} \{\varepsilon_i = 1\})$.
 For such sampling one directly derives from \eqref{eq:20200316a}
 in Remark~\ref{rmk:20190128e}\ref{rmk:20190128e_iv} (see Appendix~\ref{sec:appB}) that
 \begin{equation}
 \label{eq:20200316b}
 \beta_{1,i} =  \beta_{1} :=  1 + \frac{\eta - 1}{m-1} \bigg(
 \frac{\ensuremath{\mathsf E} \big[\big(\sum_{i=1}^m \varepsilon_i \big)^2\big]}{\ensuremath{\mathsf E} [\sum_{i=1}^m \varepsilon_i]} - 1 \bigg).
 \end{equation}
A special type of doubly uniform sampling is the \emph{$\tau$-nice} sampling in which
$\sum_{i=1}^m \varepsilon_i = \tau$ $\ensuremath{\mathsf P}$-a.s. for some $\tau \in [m]$. In this case \eqref{eq:20200316b}
reduces to
 \begin{equation}
 \label{eq:20200316c}
 \beta_{1,i} =  \beta_{1} := 1 + \frac{(\eta - 1)(\tau-1)}{m-1}.
 \end{equation}
 Now, according to Remark~\ref{rmk:20190128e}\ref{rmk:20190128e_iv},
 if we set, for every $i \in [m]$, $\nu_i = \beta_{1} L_i$, then condition \ref{eq:S1} holds.
Additionally,
 if we denote by $L^{(k)}_i$ the blockwise Lipschitz constants of the gradient of the function
 $\ensuremath \mathsf{x} \mapsto \ensuremath{{\mathsf{g}}}_k \big(\sum_{i=1}^m \ensuremath{\mathsf{U}}_{k,i} \ensuremath \mathsf{x}_i\big)$, then we derive from
 Remark~\ref{rmk:20190128e}\ref{rmk:20190128e_iii}
 and \eqref{eq:20200316c} that \ref{eq:S1} is satisfied with $\nu_i = \sum_{k \vert i \in I_k} (1 + (\tau-1)(\mathrm{card}(I_k)-1)/(m-1)) L^{(k)}_i$. This result provides possibly even smaller values
 for the parameters $(\nu_i)_{1 \leq i \leq m}$ and was given,
in the special setting of Remark~\ref{rmk:20190128e}\ref{rmk:20190128e_i},
 in \cite{Fer15,Tap18}.
\end{description}

\section{Convergence analysis}
\label{sec:conana}

In the rest of the paper, referring to Algorithm~\ref{algoRCD}, we set
\begin{equation}
\label{eq:W}
\ensuremath{{\Gammas}}^{-1} = \bigoplus_{i=1}^m \frac{1}{\gamma_i} \ensuremath{\mathsf{Id}}_i,
\quad(w_i)_{1 \leq i \leq m} = \Big( \frac{1}{\gamma_i \ensuremath \mathsf{p}_i}\Big)_{1 \leq i \leq m},
\quad\ensuremath{{\WW}} = \bigoplus_{i=1}^m w_i \ensuremath{\mathsf{Id}}_i,
\end{equation}
where $\ensuremath{\mathsf{Id}}_i$ is the identity operator on $\ensuremath \mathsf{H}_i$, and
\begin{equation}
\label{eq:20171206a}
\bar{\ensuremath x}^{\ensuremath n+1} = \big(\ensuremath{\text{\sf prox}}_{\gamma_{i} \ensuremath{{\mathsf{h}}}_i} ( x_i^{\ensuremath n}
- \gamma_i \nabla_i \ensuremath{{\mathsf{f}}}(\ensuremath x^\ensuremath n) ) \big)_{1 \leq i \leq m},
\quad
\ensuremath{{\Delta}}^\ensuremath n =  \ensuremath x^\ensuremath n - \bar{\ensuremath x}^{\ensuremath n+1}.
\end{equation}
Then, we have
\begin{equation}
\label{eq:20181114c}
\bar{\ensuremath x}^{\ensuremath n+1} = \ensuremath{\text{\sf prox}}_{\ensuremath{{\gs}}}^{\ensuremath{{\Gammas}}^{-1}}\big( \ensuremath x^\ensuremath n - \nabla^{\ensuremath{{\Gammas}}^{-1}} \ensuremath{{\mathsf{f}}}(\ensuremath x^\ensuremath n) \big),
\quad\ensuremath x^{\ensuremath n+1} = \ensuremath x^\ensuremath n + \ensuremath \varepsilon^\ensuremath n \odot (\bar{\ensuremath x}^{\ensuremath n+1} - \ensuremath x^\ensuremath n),
\end{equation}
and, recalling \eqref{eq:algoPRCD2}, that for every $i \in [m]$ such that $\varepsilon^\ensuremath n_i = 1$,
\begin{equation}
\bar{x}^{\ensuremath n+1}_{i} = \ensuremath{\text{\sf prox}}_{\gamma_{i} \ensuremath{{\mathsf{h}}}_{i}} \big( x^{\ensuremath n}_{i}
- \gamma_{i} \nabla_{i} \ensuremath{{\mathsf{f}}}(\ensuremath x^\ensuremath n) \big) =x^{\ensuremath n+1}_i,
\quad \Delta^\ensuremath n_{i} = x^\ensuremath n_{i} - x^{\ensuremath n+1}_{i}.
\end{equation}
Note that $\ensuremath x^\ensuremath n$ and $\bar{\ensuremath x}^{\ensuremath n+1}$ are functions of the random
variables $\ensuremath \varepsilon^0, \dots, \ensuremath \varepsilon^{\ensuremath n-1}$ only,
hence they are both discrete random variables, which are measurable with respect
to $\ensuremath{\mathfrak{E}}_{\ensuremath n-1}$.

\subsection{An abstract principle for stochastic convergence}

We provide an abstract convergence principle for stochastic descent algorithms
 in the same spirit of \cite[Theorem~3.10]{Sal17}. It simultaneously addresses
 the convergence of the iterates and that of the function values.

\begin{theorem}
\label{p:SMFejer}
Let $\ensuremath \mathsf{H}$ be a separable real Hilbert space with norm $\norm{\cdot}$.
Let $\ensuremath{{\mathsf{\Phi}}}\colon \ensuremath \mathsf{H} \to \left]-\infty,+\infty\right]$ be a proper, lower semicontinuous, and convex function and set $\ensuremath {\mathsf{S}}_* = \argmin \ensuremath{{\mathsf{\Phi}}}$ and $\ensuremath{{\mathsf{\Phi}}}_* = \inf \ensuremath{{\mathsf{\Phi}}}$.
Let $(\ensuremath x^{\ensuremath n})_{\ensuremath n \in \ensuremath \mathbb{N}}$ be a sequence of $\ensuremath {\mathsf{H}}$-valued random variables
such that $\ensuremath x^0 \equiv \ensuremath \xx^0 \in \dom \ensuremath{{\mathsf{\Phi}}}$ and, for every $\ensuremath n \in \ensuremath \mathbb{N}$, $\ensuremath{{\mathsf{\Phi}}}(\ensuremath x^{\ensuremath n})$ is $\ensuremath{\mathsf P}$-summable.
Consider the following conditions
\begin{enumerate}[{\rm P1}]
\item\label{SMFejer_a} $(\ensuremath{\mathsf E}[\ensuremath{{\mathsf{\Phi}}}(\ensuremath x^{\ensuremath n})])_{\ensuremath n \in \ensuremath \mathbb{N}}$ is decreasing.
\item\label{SMFejer_b}
There exist a sequence $(\mathfrak{X}_{\ensuremath n})_{\ensuremath n \in \ensuremath \mathbb{N}}$ of sub-sigma algebras of
$\mathfrak{A}$ such that, $(\forall\, \ensuremath n \in \ensuremath \mathbb{N})$
$\mathfrak{X}_{\ensuremath n} \subset \mathfrak{X}_{\ensuremath n+1}$ and
$\ensuremath x^{\ensuremath n}$ is $\mathfrak{X}_{\ensuremath n}$-measurable,
 a sequence $(\xi_\ensuremath n)_{\ensuremath n \in \ensuremath \mathbb{N}}$
of $\mathfrak{X}_{\ensuremath n}$-measurable real-valued positive random variables such that
$\sum_{\ensuremath n \in \ensuremath \mathbb{N}} \ensuremath{\mathsf E}[\xi_\ensuremath n] \leq b<+\infty$, and
 $a>0$ such that, for every $\ensuremath \xx \in \dom \ensuremath{{\mathsf{\Phi}}}$ and $\ensuremath n \in \ensuremath \mathbb{N}$,
\begin{equation}
\label{eq:SMFejer}
\ensuremath{\mathsf E}[\norm{\ensuremath x^{\ensuremath n+1} - \ensuremath \xx}^2 \,\vert\, \mathfrak{X}_{\ensuremath n} ]
\leq   \norm{\ensuremath x^{\ensuremath n} - \ensuremath \xx}^2
 + a \ensuremath{\mathsf E}[ \ensuremath{{\mathsf{\Phi}}}(\ensuremath \xx) - \ensuremath{{\mathsf{\Phi}}}(\ensuremath x^{\ensuremath n+1}) \,\vert\, \mathfrak{X}_{\ensuremath n}] +  \xi_n
\quad \ensuremath{\mathsf P}\text{-a.s.}
\end{equation}
\item\label{SMFejer_c} There exist
$(\ensuremath y^\ensuremath n)_{\ensuremath n \in \ensuremath \mathbb{N}}$ and $(\ensuremath v^\ensuremath n)_{\ensuremath n \in \ensuremath \mathbb{N}}$,
sequences of $\ensuremath {\mathsf{H}}$-valued random variables, such that
 $(\forall\, \ensuremath n \in \ensuremath \mathbb{N})$
$\ensuremath v^\ensuremath n\in \partial \ensuremath{{\mathsf{\Phi}}}(\ensuremath y^{\ensuremath n})$,
  $\ensuremath y^{\ensuremath n} - \ensuremath x^{\ensuremath n} \rightharpoonup 0$, and $\ensuremath v^{\ensuremath n} \to 0$ $\ensuremath{\mathsf P}$-a.s.
\end{enumerate}
Assume \ref{SMFejer_a} and that
$(\inf_{\ensuremath n \in \ensuremath \mathbb{N}} \ensuremath{\mathsf E}[\ensuremath{{\mathsf{\Phi}}}(\ensuremath x^\ensuremath n)]>-\infty ) \Rightarrow\ $\ref{SMFejer_b}.
 Then, the following hold.
\begin{enumerate}[{\rm (i)}]
\item\label{SMFejer_i} $\ensuremath{\mathsf E}[\ensuremath{{\mathsf{\Phi}}}(\ensuremath x^{\ensuremath n})] \to \ensuremath{{\mathsf{\Phi}}}_*$.
\item\label{SMFejer_ii} Suppose that $\ensuremath {\mathsf{S}}_*\neq \varnothing$.
Then $\ensuremath{\mathsf E}[\ensuremath{{\mathsf{\Phi}}}(\ensuremath x^{\ensuremath n})] -  \ensuremath{{\mathsf{\Phi}}}_* = o(1/\ensuremath n)$ and,
\begin{equation*}
(\forall\, \ensuremath n \in \ensuremath \mathbb{N},\ensuremath n \geq 1)\quad
\ensuremath{\mathsf E}[\ensuremath{{\mathsf{\Phi}}}(\ensuremath x^{\ensuremath n})] - \ensuremath{{\mathsf{\Phi}}}_*  \leq
\bigg[ \frac{ \mathrm{dist}^2(\ensuremath \xx^0, \ensuremath {\mathsf{S}}_*) }{a}
+ \frac{b}{a}\bigg] \frac{1}{n}.
\end{equation*}
\item\label{SMFejer_iii}
Suppose that \ref{SMFejer_c} holds and $\ensuremath {\mathsf{S}}_*\neq \varnothing$.
Then, there exists a random variable $\ensuremath x_*$ taking values in $\ensuremath {\mathsf{S}}_*$ such that
$\ensuremath x^\ensuremath n \rightharpoonup \ensuremath x_*$ $\ensuremath{\mathsf P}$-a.s.
\end{enumerate}
\end{theorem}
\begin{proof}
Taking the expectation in \eqref{eq:SMFejer}, we obtain
\begin{equation}
\label{eq:SMFejer2}
a ( \ensuremath{\mathsf E}[\ensuremath{{\mathsf{\Phi}}}(\ensuremath x^{\ensuremath n+1})] -\ensuremath{{\mathsf{\Phi}}}(\ensuremath \xx)) \leq
 \ensuremath{\mathsf E}[\norm{\ensuremath x^{\ensuremath n} - \ensuremath \xx}^2 ] - \ensuremath{\mathsf E}[\norm{\ensuremath x^{\ensuremath n+1} - \ensuremath \xx}^2 ]
 +  \ensuremath{\mathsf E}[\xi_n].
\end{equation}

\ref{SMFejer_i}:
Since $(\ensuremath{\mathsf E}[\ensuremath{{\mathsf{\Phi}}}(\ensuremath x^\ensuremath n)])_{n \in \ensuremath \mathbb{N}}$ is decreasing,
$\ensuremath{\mathsf E}[\ensuremath{{\mathsf{\Phi}}}(\ensuremath x^\ensuremath n)]\to\inf_{n \in \ensuremath \mathbb{N}} \ensuremath{\mathsf E}[\ensuremath{{\mathsf{\Phi}}}(\ensuremath x^\ensuremath n)] \geq \ensuremath{{\mathsf{\Phi}}}_*$. Thus,
the statement is true if $\inf_{n \in \ensuremath \mathbb{N}} \ensuremath{\mathsf E}[\ensuremath{{\mathsf{\Phi}}}(\ensuremath x^\ensuremath n)] = -\infty$.
Suppose that  $\inf_{n \in \ensuremath \mathbb{N}} \ensuremath{\mathsf E}[\ensuremath{{\mathsf{\Phi}}}(\ensuremath x^\ensuremath n)]>- \infty$ and let $\ensuremath \xx \in \dom \ensuremath{{\mathsf{\Phi}}}$. Then, \ref{SMFejer_b} holds and
the right hand side of \eqref{eq:SMFejer2}, being summable, converges to zero.
Therefore,
$\ensuremath{{\mathsf{\Phi}}}_* \leq \lim_{n \to +\infty} \ensuremath{\mathsf E}[\ensuremath{{\mathsf{\Phi}}}(\ensuremath x^{\ensuremath n+1})] \leq \ensuremath{{\mathsf{\Phi}}}(\ensuremath \xx)$. Since $\ensuremath \xx$
is arbitrary in $\dom \ensuremath{{\mathsf{\Phi}}}$, $\ensuremath{\mathsf E}[\ensuremath{{\mathsf{\Phi}}}(\ensuremath x^{\ensuremath n})]\to \ensuremath{{\mathsf{\Phi}}}_*$.

\ref{SMFejer_ii}:
Let $\ensuremath \xx \in \ensuremath {\mathsf{S}}_*$. Then,  $\inf_{n \in \ensuremath \mathbb{N}} \ensuremath{\mathsf E}[\ensuremath{{\mathsf{\Phi}}}(\ensuremath x^\ensuremath n)]\geq \ensuremath{{\mathsf{\Phi}}}(\ensuremath \xx)>- \infty$.
Hence \ref{SMFejer_b} holds and \eqref{eq:SMFejer2} yields
\begin{align*}
a \sum_{\ensuremath n \in \ensuremath \mathbb{N}} \big(\ensuremath{\mathsf E}[\ensuremath{{\mathsf{\Phi}}}(\ensuremath x^{\ensuremath n+1})] - \ensuremath{{\mathsf{\Phi}}}_* \big)
&\leq \ensuremath{\mathsf E} \big[ \norm{\ensuremath x^0 -  \ensuremath \xx}^2\big]
+ \sum_{\ensuremath n \in \ensuremath \mathbb{N}} \ensuremath{\mathsf E}[\xi_{\ensuremath n}]
\leq \norm{\ensuremath \xx^0 -  \ensuremath \xx}^2
+ b.
\end{align*}
Therefore, $\sum_{\ensuremath n \in \ensuremath \mathbb{N}} (\ensuremath{\mathsf E}[\ensuremath{{\mathsf{\Phi}}}(\ensuremath x^{\ensuremath n+1})] - \ensuremath{{\mathsf{\Phi}}}_* )
\leq (1/a) \norm{\ensuremath \xx^0 -  \ensuremath \xx}^2  + b/a$.
Since $(\ensuremath{\mathsf E}[\ensuremath{{\mathsf{\Phi}}}(\ensuremath x^{\ensuremath n+1})] - \ensuremath{{\mathsf{\Phi}}}_*)_{\ensuremath n \in \ensuremath \mathbb{N}}$ is decreasing,
the statement follows from Fact~\ref{lem:monosequence}.

\ref{SMFejer_iii}:
Let $\ensuremath \xx \in \ensuremath {\mathsf{S}}_*$. Then \ref{SMFejer_b} holds and,
since $\ensuremath{{\mathsf{\Phi}}}(\ensuremath \xx) \leq \ensuremath{{\mathsf{\Phi}}}(\ensuremath x^{\ensuremath n+1})$,
we derive from \eqref{eq:SMFejer} that,
\begin{equation}
(\forall\, \ensuremath n \in \ensuremath \mathbb{N})\quad
\ensuremath{\mathsf E}\big[\norm{\ensuremath x^{\ensuremath n+1} - \ensuremath \xx}^2\,\vert\,
\mathfrak{X}_{\ensuremath n}\big]
\leq
 \norm{\ensuremath x^\ensuremath n -  \ensuremath \xx}^2 + \xi_\ensuremath n \quad \ensuremath{\mathsf P}\text{-a.s.}
\end{equation}
Note that  $\xi_\ensuremath n$ and $ \norm{\ensuremath x^\ensuremath n -  \ensuremath \xx}^2$
are $\mathfrak{X}_{\ensuremath n}$-measurable.
Moreover $\ensuremath{\mathsf E}[\sum_{\ensuremath n \in \ensuremath \mathbb{N}} \xi_\ensuremath n] = \sum_{\ensuremath n \in \ensuremath \mathbb{N}} \ensuremath{\mathsf E}[\xi_\ensuremath n]<+\infty$
and hence $\sum_{\ensuremath n \in \ensuremath \mathbb{N}} \xi_\ensuremath n<+\infty$ $\ensuremath{\mathsf P}$-a.s.
Therefore $(\ensuremath x^\ensuremath n)_{\ensuremath n \in \ensuremath \mathbb{N}}$ is a stochastic quasi-Fej\'er sequence with respect to $\ensuremath {\mathsf{S}}_*$
\cite{Erm69}.
Then, in view of \cite[Proposition~2.3(iv)]{Com15} it is sufficient to prove that
 the weak limit points of
$(\ensuremath x^\ensuremath n)_{\ensuremath n \in \ensuremath \mathbb{N}}$ are contained in $\ensuremath {\mathsf{S}}_*$ $\ensuremath{\mathsf P}$-a.s.
By assumption \ref{SMFejer_c} there exist two sequences of
$\ensuremath {\mathsf{H}}$-valued random variables
 $(\ensuremath y^\ensuremath n)_{\ensuremath n \in \ensuremath \mathbb{N}}$
and $(\ensuremath v^\ensuremath n)_{\ensuremath n \in \ensuremath \mathbb{N}}$  and $\tilde{\Omega} \subset \Omega$, $\ensuremath{\mathsf P}(\tilde{\Omega})=1$
such that, for every $\omega \in \tilde{\Omega}$,
 $\ensuremath v^\ensuremath n(\omega)\in \partial \ensuremath{{\mathsf{\Phi}}}(\ensuremath y^{\ensuremath n}(\omega))$,
$\ensuremath y^{\ensuremath n}(\omega) - \ensuremath x^{\ensuremath n}(\omega) \rightharpoonup 0$, $\ensuremath v^{\ensuremath n}(\omega) \to 0$.
Let $\omega \in \tilde{\Omega}$ and let $(\ensuremath x^{n_k}(\omega))_{\ensuremath n \in \ensuremath \mathbb{N}}$
be a subsequence of $(\ensuremath x^\ensuremath n(\omega))_{\ensuremath n \in \ensuremath \mathbb{N}}$ such that
 $\ensuremath x^{n_k}(\omega) \rightharpoonup \bar{\ensuremath \xx}$, for some $\bar{\ensuremath \xx} \in \ensuremath {\mathsf{H}}$.
Then,
\begin{equation*}
\ensuremath y^{n_k}(\omega) \rightharpoonup \bar{\ensuremath \xx}, \ensuremath v^{n_k}(\omega) \to 0, \
\ensuremath v^{n_k}(\omega) \in \partial \ensuremath{{\mathsf{\Phi}}}(\ensuremath y^{n_k}(\omega)).
\end{equation*}
Since $\partial \ensuremath{{\mathsf{\Phi}}}$ is weakly-strongly closed \cite{book1}, we have $0 \in \partial \ensuremath{{\mathsf{\Phi}}}(\bar{\ensuremath \xx})$,
so $\bar{\ensuremath \xx} \in \ensuremath {\mathsf{S}}_*$.
\end{proof}

\begin{remark}
 Inequalities similar to \eqref{eq:SMFejer} appear implicitly in
the analysis of several deterministic and stochastic
 algorithms \cite{Ber11,Kiw06,Nem09}, to get rate of convergence for the function values. Moreover, \eqref{eq:SMFejer} is related also to the concept introduced in \cite{Lin18},
 in a deterministic setting.
\end{remark}

\subsection{Convergence under convexity and strong convexity assumptions}

In this section we address the convergence of Algorithm~\ref{algoRCD} in the convex and strongly convex case. The main results consist in the $o(1/\ensuremath n)$ rate of convergence
for the mean of the function values and in the
almost sure weak convergence of the iterates.
We start by recalling a standard result
(see \cite[Lemma~3.12(iii)]{Sal17}). Here we give a slightly
more general version, including the moduli of strong convexity.
The proof is given in Appendix~\ref{sec:appB} for reader's convenience.
\begin{lemma}
\label{lem:20190313b}
Let $\ensuremath \mathsf{H}$ be a real Hilbert space.
Let $\varphi\colon \ensuremath \mathsf{H} \to \ensuremath \mathbb{R}$ be differentiable and convex with modulus of strong convexity $\mu_\varphi \geq 0$
and $\psi\colon \ensuremath \mathsf{H}\to \left]-\infty,+\infty\right]$ be proper, lower semicontinuous,
and convex with modulus of strong convexity $\mu_\psi \geq 0$.
Let $\ensuremath \mathsf{x} \in \ensuremath \mathsf{H}$ and set $\ensuremath \mathsf{x}^+ = \ensuremath{\text{\sf prox}}_{\psi}(\ensuremath \mathsf{x} - \nabla \varphi(\ensuremath \mathsf{x}))$.
Then,
for every $\ensuremath \mathsf{z} \in \ensuremath \mathsf{H}$,
\begin{align*}
\nonumber(1 + \mu_\psi) \scalarp{\ensuremath \mathsf{x} - \ensuremath \mathsf{x}^+ , \ensuremath \mathsf{z} - \ensuremath \mathsf{x}}
&\leq  \Big((\varphi + \psi)(\ensuremath \mathsf{z}) - (\varphi + \psi)(\ensuremath \mathsf{x})  - \frac{\mu_\varphi + \mu_\psi}{2} \norm{\ensuremath \mathsf{z} - \ensuremath \mathsf{x}}^2 \Big)\\[1ex]
&\ + \big(\psi(\ensuremath \mathsf{x}) - \psi(\ensuremath \mathsf{x}^+) +  \scalarp{\nabla \varphi(\ensuremath \mathsf{x}), \ensuremath \mathsf{x} - \ensuremath \mathsf{x}^+}\big)
-  \Big(1 + \frac{\mu_\psi}{2}\Big) \norm{\ensuremath \mathsf{x} - \ensuremath \mathsf{x}^+}^2.
\end{align*}
\end{lemma}

\begin{proposition}
\label{p:20190610a}
 Let \ref{eq:A1}--\ref{eq:A3} be satisfied.
Let $(\nu_{i})_{1 \leq i \leq m} \in \ensuremath \mathbb{R}^m_{++}$ and suppose that
\ref{eq:S1} holds.
Let $(\ensuremath x^\ensuremath n)_{\ensuremath n \in \ensuremath \mathbb{N}}$ be generated by  Algorithm~\ref{algoRCD}
with, for every $i \in [m]$, $\gamma_i < 2/\nu_i$.
Set $\delta = \max_{1 \leq i \leq m} \gamma_i  \nu_i$
 and
$\ensuremath \mathsf{p}_{\min} = \min_{1 \leq i \leq m} \ensuremath \mathsf{p}_i$.
Let $\ensuremath \mathsf{\Gamma}^{-1}$ be as in \eqref{eq:W} and
 $\mu_{\ensuremath \mathsf{\Gamma}^{-1}}$ and $\sigma_{\ensuremath \mathsf{\Gamma}^{-1}}$ be the moduli of strong convexity of
$\ensuremath{{\mathsf{f}}}$ and $\ensuremath{{\gs}}$ respectively, in the norm $\norm{\cdot}_{\ensuremath \mathsf{\Gamma}^{-1}}$.
Set $\ensuremath{{\mathsf{F}}} = \ensuremath{{\mathsf{f}}}+\ensuremath{{\gs}}$.
Then,
\begin{align}
\nonumber
(1 + \sigma_{\ensuremath \mathsf{\Gamma}^{-1}})) \scalarp{\ensuremath x^{\ensuremath n} - \bar{\ensuremath x}^{\ensuremath n+1}, \ensuremath \xx - \ensuremath x^{\ensuremath n}}_{\ensuremath \mathsf{\Gamma}^{-1}}
&\leq  \frac{1}{\ensuremath \mathsf{p}_{\min}}\ensuremath{\mathsf E} \big[ \ensuremath{{\mathsf{F}}}(\ensuremath x^n) - \ensuremath{{\mathsf{F}}}(\ensuremath x^{n+1}) \,\vert\, \ensuremath{\mathfrak{E}}_{\ensuremath n-1}  \big]  \\[1ex]
\nonumber&\qquad + \bigg( \ensuremath{{\mathsf{F}}}(\ensuremath \xx) -\ensuremath{{\mathsf{F}}}(\ensuremath x^\ensuremath n)
- \frac{\mu_{\ensuremath \mathsf{\Gamma}^{-1}} + \sigma_{\ensuremath \mathsf{\Gamma}^{-1}}}{2} \norm{\ensuremath x^\ensuremath n - \ensuremath \xx}_{\ensuremath{{\Gammas}}^{-1}}^2 \bigg)\\[1ex]
\label{eq:20190610a}&\qquad +
\frac{\delta - 2 - \sigma_{\ensuremath \mathsf{\Gamma}^{-1}}}{2} \norm{\ensuremath x^{\ensuremath n} - \bar{\ensuremath x}^{\ensuremath n+1}}_{\ensuremath{{\Gammas}}^{-1}}^2.
\end{align}
\end{proposition}
\begin{proof}
Let $\ensuremath \xx \in \dom \ensuremath{{\mathsf{F}}}$ and $\ensuremath n \in \ensuremath \mathbb{N}$. Since for all $\ensuremath \vv \in \ensuremath {\mathsf{H}}$, $\scalarp{\nabla^{\ensuremath{{\Gammas}}^{-1}} \ensuremath{{\mathsf{f}}} (\ensuremath x^n), \ensuremath \vv}_{\ensuremath{{\Gammas}}^{-1}} =
 \scalarp{\nabla \ensuremath{{\mathsf{f}}}(\ensuremath x^n), \ensuremath \vv}$, we
derive from Lemma~\ref{lem:20190313b},
 written in the norm $\norm{\cdot}_{\ensuremath{{\Gammas}}^{-1}}$, and
 \eqref{eq:20181114c} that
\begin{align}
\nonumber
(1 + \sigma_{\ensuremath \mathsf{\Gamma}^{-1}})) \scalarp{\ensuremath x^{\ensuremath n} - \bar{\ensuremath x}^{\ensuremath n+1}, \ensuremath \xx - \ensuremath x^{\ensuremath n}}_{\ensuremath \mathsf{\Gamma}^{-1}}
&\leq  \big(\ensuremath{{\gs}}(\ensuremath x^n) - \ensuremath{{\gs}}(\bar{\ensuremath x}^{n+1}) +  \scalarp{\nabla \ensuremath{{\mathsf{f}}}(\ensuremath x^n), \ensuremath x^\ensuremath n - \bar{\ensuremath x}^{\ensuremath n+1}}\big)  \\[0.8ex]
\nonumber&\qquad + \Big( \ensuremath{{\mathsf{F}}}(\ensuremath \xx) -\ensuremath{{\mathsf{F}}}(\ensuremath x^\ensuremath n)
- \frac{\mu_{\ensuremath \mathsf{\Gamma}^{-1}} + \sigma_{\ensuremath \mathsf{\Gamma}^{-1}}}{2} \norm{\ensuremath x^\ensuremath n - \ensuremath \xx}_{\ensuremath{{\Gammas}}^{-1}}^2 \Big)\\[0.8ex]
\label{eq:20190313f}&\qquad - \Big(1 + \frac{\sigma_{\ensuremath \mathsf{\Gamma}^{-1}}}{2} \Big)\norm{\ensuremath x^{\ensuremath n} - \bar{\ensuremath x}^{\ensuremath n+1}}_{\ensuremath{{\Gammas}}^{-1}}^2.
\end{align}
Now, we majorize $\ensuremath{{\gs}}(\ensuremath x^n) - \ensuremath{{\gs}}(\bar{\ensuremath x}^{n+1}) +  \scalarp{\nabla \ensuremath{{\mathsf{f}}}(\ensuremath x^n), \ensuremath x^n - \bar{\ensuremath x}^{n+1}}$. By
Fact~\ref{f:20190112d} and  Fact~\ref{lem:20181024a}, we have
\begin{align*}
\ensuremath{{\gs}}(\ensuremath x^n) - \ensuremath{{\gs}}(\bar{\ensuremath x}^{n+1}) &+  \scalarp{\nabla \ensuremath{{\mathsf{f}}}(\ensuremath x^n), \ensuremath x^n - \bar{\ensuremath x}^{n+1}}\\
& = \ensuremath{\mathsf E} \bigg[ \sum_{i=1}^m \frac{\varepsilon^n_i}{\ensuremath \mathsf{p}_i}
\Big( \ensuremath{{\mathsf{h}}}_i(x_i^n) - \ensuremath{{\mathsf{h}}}_i(\bar{x}_i^{n+1}) + \scalarp{\nabla_i \ensuremath{{\mathsf{f}}}(\ensuremath x^n), x_i^n - \bar{x}_i^{n+1}} \Big) \,\Big\vert\, \ensuremath{\mathfrak{E}}_{\ensuremath n-1} \bigg].
\end{align*}
Moreover,
\begin{align*}
\sum_{i=1}^m \frac{\varepsilon^n_i}{\ensuremath \mathsf{p}_i}
&\Big( \ensuremath{{\mathsf{h}}}_i(x_i^n) - \ensuremath{{\mathsf{h}}}_i(\bar{x}_i^{n+1}) + \scalarp{\nabla_i \ensuremath{{\mathsf{f}}}(\ensuremath x^n), x_i^n - \bar{x}_i^{n+1}} \Big)\\
& =  \sum_{i=1}^m \frac{1}{\ensuremath \mathsf{p}_i}
\Big( \ensuremath{{\mathsf{h}}}_i(x_i^n) - \ensuremath{{\mathsf{h}}}_i(x_i^{n+1})
+ \scalarp{\nabla_i \ensuremath{{\mathsf{f}}}(\ensuremath x^n), x_i^n - x_i^{n+1}} \Big)\\[0.8ex]
& = \frac{1}{\ensuremath \mathsf{p}_{\min}}
\Big( \ensuremath{{\gs}}(\ensuremath x^n) - \ensuremath{{\gs}}(\ensuremath x^{n+1}) +  \scalarp{\nabla \ensuremath{{\mathsf{f}}}(\ensuremath x^n), \ensuremath x^\ensuremath n - \ensuremath x^{\ensuremath n+1}} \Big)\\[0.8ex]
 &\qquad - \sum_{i=1}^m  \bigg(\underbrace{\frac{1}{\ensuremath \mathsf{p}_{\min}} - \frac{1}{\ensuremath \mathsf{p}_{i}}}_{\geq 0} \bigg)
\Big( \ensuremath{{\mathsf{h}}}_i(x_i^n) - \ensuremath{{\mathsf{h}}}_i(x_i^{n+1})
+ \scalarp{\nabla_i \ensuremath{{\mathsf{f}}}(\ensuremath x^n), x_i^n - x_i^{n+1}} \Big)\\[0.8ex]
&\leq \frac{1}{\ensuremath \mathsf{p}_{\min}}
\Big( \ensuremath{{\gs}}(\ensuremath x^n) - \ensuremath{{\gs}}(\ensuremath x^{n+1})
+  \scalarp{\nabla \ensuremath{{\mathsf{f}}}(\ensuremath x^n), \ensuremath x^\ensuremath n - \ensuremath x^{\ensuremath n+1}} \Big)\\[0.8ex]
 &\qquad - \bigg( 1 + \frac{\sigma_{\ensuremath \mathsf{\Gamma}^{-1}}}{2} \bigg)\sum_{i=1}^m \bigg(\frac{1}{\ensuremath \mathsf{p}_{\min}} - \frac{1}{\ensuremath \mathsf{p}_{i}} \bigg)
 \frac{\varepsilon^\ensuremath n_i}{\gamma_i} \norm{\Delta_i^\ensuremath n}^2,
\end{align*}
where in the last inequality we used that
\begin{equation}
\label{eq:20181108b}
-\big(\ensuremath{{\mathsf{h}}}_i(x_i^{\ensuremath n}) - \ensuremath{{\mathsf{h}}}_i(x_i^{\ensuremath n+1}) + \scalarp{\nabla_i \ensuremath{{\mathsf{f}}}(\ensuremath x^\ensuremath n), x_i^\ensuremath n - x_i^{\ensuremath n+1}}
 \big)
\leq - \frac{\varepsilon^\ensuremath n_i}{\gamma_i}
 \bigg( 1 + \frac{\sigma_{\ensuremath \mathsf{\Gamma}^{-1}}}{2} \bigg)\norm{\Delta_i^\ensuremath n}^2,
\end{equation}
which was obtained from
\eqref{eq:20181112a} with
\begin{equation*}
\ensuremath \mathsf{x}_i = x_i^\ensuremath n, \quad
\ensuremath \mathsf{y}_i = x_i^{\ensuremath n+1},\quad
\ensuremath \mathsf{v}_i = \frac{x^\ensuremath n_{i} - x^{\ensuremath n+1}_{i}}{\gamma_{i}} - \nabla_{i} \ensuremath{{\mathsf{f}}} (\ensuremath x^\ensuremath n)
\in \partial \ensuremath{{\mathsf{h}}}_{i}(x^{\ensuremath n+1}_{i}),\quad\text{for }\varepsilon^\ensuremath n_i=1.
\end{equation*}
Therefore,
\begin{align}
\nonumber \ensuremath{{\gs}}(\ensuremath x^n) - \ensuremath{{\gs}}(\bar{\ensuremath x}^{n+1})
&+  \scalarp{\nabla \ensuremath{{\mathsf{f}}}(\ensuremath x^n), \ensuremath x^n - \bar{\ensuremath x}^{n+1}}\\
\nonumber&\leq \frac{1}{\ensuremath \mathsf{p}_{\min}}\ensuremath{\mathsf E} \big[ \ensuremath{{\gs}}(\ensuremath x^n) - \ensuremath{{\gs}}(\ensuremath x^{n+1}) +  \scalarp{\nabla \ensuremath{{\mathsf{f}}}(\ensuremath x^n), \ensuremath x^n - \ensuremath x^{n+1}} \,\vert\, \ensuremath{\mathfrak{E}}_{\ensuremath n-1}  \big]\\[1ex]
\label{eq:20190313i}&\qquad - \frac{1}{\ensuremath \mathsf{p}_{\min}}\bigg( 1 + \frac{\sigma_{\ensuremath \mathsf{\Gamma}^{-1}}}{2} \bigg) \sum_{i=1}^m
 \frac{\ensuremath \mathsf{p}_i}{\gamma_i} \norm{\Delta_i^\ensuremath n}^2
 + \bigg( 1 + \frac{\sigma_{\ensuremath \mathsf{\Gamma}^{-1}}}{2} \bigg) \norm{\bar{\ensuremath x}^{n+1} - \ensuremath x^n}^2_{\ensuremath{{\Gammas}}^{-1}}.
\end{align}
Next, it follows from \eqref{eq:20181114c}, \ref{eq:S1}, and Fact~\ref{f:20190112d} that
\begin{equation*}
\ensuremath{\mathsf E}[\scalarp{\nabla \ensuremath{{\mathsf{f}}}(\ensuremath x^\ensuremath n), \ensuremath x^{\ensuremath n} - \ensuremath x^{\ensuremath n+1}} \,\vert\, \ensuremath{\mathfrak{E}}_{\ensuremath n-1}]
 \leq \ensuremath{\mathsf E}[\ensuremath{{\mathsf{f}}}(\ensuremath x^\ensuremath n) -  \ensuremath{{\mathsf{f}}}(\ensuremath x^{\ensuremath n+1}) \,\vert\, \ensuremath{\mathfrak{E}}_{\ensuremath n-1}]
+ \frac{1}{2}\sum_{i=1}^m \ensuremath \mathsf{p}_i \nu_i \norm{\Delta_i^\ensuremath n}^2.
\end{equation*}
Then, we derive from \eqref{eq:20190313i} that
\begin{align*}
 \ensuremath{{\gs}}(\ensuremath x^n) - \ensuremath{{\gs}}(\bar{\ensuremath x}^{n+1})
&+  \scalarp{\nabla \ensuremath{{\mathsf{f}}}(\ensuremath x^n), \ensuremath x^n - \bar{\ensuremath x}^{n+1}}\\
&\leq \frac{1}{\ensuremath \mathsf{p}_{\min}}\ensuremath{\mathsf E} \big[ \ensuremath{{\mathsf{F}}}(\ensuremath x^n) - \ensuremath{{\mathsf{F}}}(\ensuremath x^{n+1}) \,\vert\, \ensuremath{\mathfrak{E}}_{\ensuremath n-1}  \big]\\[1ex]
&\qquad - \frac{1}{2\ensuremath \mathsf{p}_{\min}}\sum_{i=1}^m  \bigg( 2 + \sigma_{\ensuremath \mathsf{\Gamma}^{-1}}
- \gamma_i \nu_i \bigg)
 \frac{\ensuremath \mathsf{p}_i}{\gamma_i} \norm{\Delta_i^\ensuremath n}^2
 + \bigg( 1 + \frac{\sigma_{\ensuremath \mathsf{\Gamma}^{-1}}}{2} \bigg) \norm{\bar{\ensuremath x}^{n+1} - \ensuremath x^n}^2_{\ensuremath{{\Gammas}}^{-1}}.
\end{align*}
The statement follows from \eqref{eq:20190313f}, considering that
\begin{align*}
\frac{1}{\ensuremath \mathsf{p}_{\min}} \sum_{i=1}^m \bigg(\gamma_i\nu_i - 2  - \sigma_{\ensuremath \mathsf{\Gamma}^{-1}}\bigg)
\frac{\ensuremath \mathsf{p}_i}{\gamma_i}\norm{\Delta_i^\ensuremath n}^2
&\leq
\underbrace{\frac{\delta - 2 - \sigma_{\ensuremath \mathsf{\Gamma}^{-1}}}{\ensuremath \mathsf{p}_{\min}}}_{\leq0} \sum_{i=1}^m
\frac{\ensuremath \mathsf{p}_i}{\gamma_i}\norm{\Delta_i^\ensuremath n}^2\\
&\leq
\frac{\delta - 2 - \sigma_{\ensuremath \mathsf{\Gamma}^{-1}}}{\ensuremath \mathsf{p}_{\min}} \sum_{i=1}^m
\frac{\ensuremath \mathsf{p}_{\min}}{\gamma_i}\norm{\Delta_i^\ensuremath n}^2\\
&=(\delta - 2 - \sigma_{\ensuremath \mathsf{\Gamma}^{-1}}) \norm{\ensuremath x^n - \bar{\ensuremath x}^{n+1}}^2_{\ensuremath{{\Gammas}}^{-1}}.
\qedhere
\end{align*}
\end{proof}

\begin{proposition}
\label{p:20190313c}
 Let \ref{eq:A1}--\ref{eq:A3} be satisfied.
 Let $\ensuremath{{\Gammas}}^{-1}$ and $\ensuremath{{\WW}}$ be as in \eqref{eq:W} and
 $(\ensuremath x^\ensuremath n)_{\ensuremath n \in \ensuremath \mathbb{N}}$ be generated by Algorithm~\ref{algoRCD}.
Let $n \in \ensuremath \mathbb{N}$ and $\ensuremath x$ be an $\ensuremath {\mathsf{H}}$-valued  random variable
which is measurable w.r.t. $\ensuremath{\mathfrak{E}}_{\ensuremath n-1}$. Then
\begin{equation}
\label{eq:20190313c}
\ensuremath{\mathsf E}[\norm{\ensuremath x^{\ensuremath n+1} - \ensuremath x}_\ensuremath{{\WW}}^2 \,\vert\, \ensuremath{\mathfrak{E}}_{\ensuremath n-1} ] -  \norm{\ensuremath x^\ensuremath n - \ensuremath x}^2_\ensuremath{{\WW}}=
\norm{\bar{\ensuremath x}^{\ensuremath n+1} - \ensuremath x}^2_{\ensuremath{{\Gammas}}^{-1}}
 - \norm{\ensuremath x^\ensuremath n - \ensuremath x}^2_{\ensuremath{{\Gammas}}^{-1}}
\end{equation}
and
$\ensuremath{\mathsf E}[\norm{\ensuremath x^{\ensuremath n+1} - \ensuremath x^\ensuremath n}_\ensuremath{{\WW}}^2 \,\vert\, \ensuremath{\mathfrak{E}}_{\ensuremath n-1} ]
= \norm{\bar{\ensuremath x}^{\ensuremath n+1} - \ensuremath x^\ensuremath n}^2_{\ensuremath{{\Gammas}}^{-1}}$.
\end{proposition}
\begin{proof}
If follows from \eqref{eq:20181114c}, Fact~\ref{f:20190112d},
and Fact~\ref{lem:20181024a} that
\begin{align*}
\ensuremath{\mathsf E}[\norm{\ensuremath x^{\ensuremath n+1} - \ensuremath x}_\ensuremath{{\WW}}^2 \,\vert\, \ensuremath{\mathfrak{E}}_{\ensuremath n-1} ]
&= \ensuremath{\mathsf E}\bigg[ \sum_{i=1}^m \frac{1}{\gamma_i \ensuremath \mathsf{p}_i} \norm{x^{\ensuremath n+1}_i - x_i}^2
\,\big\vert\, \ensuremath{\mathfrak{E}}_{\ensuremath n-1}  \bigg]\\
& = \ensuremath{\mathsf E}\bigg[ \sum_{i=1}^m \frac{\varepsilon_i^\ensuremath n}{\gamma_i \ensuremath \mathsf{p}_i} \norm{\bar{x}^{\ensuremath n+1}_i - x_i}^2
\,\big\vert\, \ensuremath{\mathfrak{E}}_{\ensuremath n-1}  \bigg]
+ \ensuremath{\mathsf E}\bigg[ \sum_{i=1}^m \frac{1 -\varepsilon_i^\ensuremath n}{\gamma_i \ensuremath \mathsf{p}_i} \norm{x^{\ensuremath n}_i - x_i}^2
\,\big\vert\, \ensuremath{\mathfrak{E}}_{\ensuremath n-1}  \bigg]\\
& = \norm{\bar{\ensuremath x}^{\ensuremath n+1} - \ensuremath x}^2_{\ensuremath{{\Gammas}}^{-1}}
+ \norm{\ensuremath x^\ensuremath n - \ensuremath x}^2_W - \norm{\ensuremath x^\ensuremath n - \ensuremath x}^2_{\ensuremath{{\Gammas}}^{-1}}.
\end{align*}
The second equation follows from \eqref{eq:20190313c}, by choosing $\ensuremath x = \ensuremath x^\ensuremath n$.
\end{proof}

The following result is a stochastic version of \cite[Proposition~3.15]{Sal17}.

\begin{proposition}
 \label{p:20180927a}
 Let \ref{eq:A1}--\ref{eq:A3} be satisfied.
Let $(\nu_{i})_{1 \leq i \leq m} \in \ensuremath \mathbb{R}^m_{++}$ and suppose that
\ref{eq:S1} holds.
Let $(\ensuremath x^\ensuremath n)_{\ensuremath n \in \ensuremath \mathbb{N}}$ be generated by  Algorithm~\ref{algoRCD}
with, for every $i \in [m]$, $\gamma_i < 2/\nu_i$.
Set $\delta = \max_{1 \leq i \leq m} \gamma_i  \nu_i$
 and
$\ensuremath \mathsf{p}_{\min} = \min_{1 \leq i \leq m} \ensuremath \mathsf{p}_i$.
Let $\ensuremath \mathsf{\Gamma}^{-1}$ and $\ensuremath{{\WW}}$ be as in \eqref{eq:W} and
 $\mu_{\ensuremath \mathsf{\Gamma}^{-1}}$ and $\sigma_{\ensuremath \mathsf{\Gamma}^{-1}}$ be the moduli of strong convexity of
$\ensuremath{{\mathsf{f}}}$ and $\ensuremath{{\gs}}$ respectively, in the norm $\norm{\cdot}_{\ensuremath \mathsf{\Gamma}^{-1}}$.
Set $\ensuremath{{\mathsf{F}}} = \ensuremath{{\mathsf{f}}}+\ensuremath{{\gs}}$.
Then, the following hold.
\begin{enumerate}[{\rm (i)}]
\item\label{p:20180927a_0} $( \ensuremath{\mathsf E}[\ensuremath{{\mathsf{F}}}(\ensuremath x^\ensuremath n)])_{\ensuremath n \in \ensuremath \mathbb{N}}$ is descreasing.
\item\label{p:20180927a_i} Suppose that $\inf_{\ensuremath n \in \ensuremath \mathbb{N}} \ensuremath{\mathsf E}[\ensuremath{{\mathsf{F}}}(\ensuremath x^\ensuremath n)]>0$.
Then,
\begin{equation*}
\sum_{\ensuremath n \in \ensuremath \mathbb{N}} \norm{\bar{\ensuremath x}^{\ensuremath n+1} - \ensuremath x^\ensuremath n}_{\ensuremath \mathsf{\Gamma}^{-1}}^2
= \sum_{\ensuremath n \in \ensuremath \mathbb{N}} \ensuremath{\mathsf E} \big[ \norm{\ensuremath x^\ensuremath n - \ensuremath x^{\ensuremath n+1}}_\ensuremath{{\WW}}^2 \big\vert
\ensuremath{\mathfrak{E}}_{\ensuremath n-1} \big]
<+\infty\quad\ensuremath{\mathsf P}\text{ a.s.}
\end{equation*}
\item\label{p:20180927a_ii} For every $\ensuremath n \in \ensuremath \mathbb{N}$ and every $\ensuremath \xx \in \dom \ensuremath{{\mathsf{F}}}$
\begin{align}
\nonumber
\nonumber(1 + \sigma_{\ensuremath \mathsf{\Gamma}^{-1}})\ensuremath{\mathsf E}[\norm{\ensuremath x^{\ensuremath n+1} - \ensuremath \xx}_\ensuremath{{\WW}}^2 \,\vert\, \ensuremath{\mathfrak{E}}_{\ensuremath n-1} ] &
\leq (1 + \sigma_{\ensuremath \mathsf{\Gamma}^{-1}})  \norm{\ensuremath x^{\ensuremath n} - \ensuremath \xx}_{\ensuremath{{\WW}}}^2
\\[1ex]
\nonumber&\ + 2 \bigg(  \ensuremath{{\mathsf{F}}}(\ensuremath \xx) -\ensuremath{{\mathsf{F}}}(\ensuremath x^\ensuremath n)
+ \frac{\mu_{\ensuremath \mathsf{\Gamma}^{-1}} + \sigma_{\ensuremath \mathsf{\Gamma}^{-1}}}{2} \norm{\ensuremath x^\ensuremath n - \ensuremath \xx}_{\ensuremath{{\Gammas}}^{-1}}^2 \bigg)\\[1ex]
\nonumber&\ +\frac{2}{\ensuremath \mathsf{p}_{\min}}
\Big(\frac{(\delta-1)_+}{2+\sigma_{\ensuremath \mathsf{\Gamma}^{-1}}-\delta}+1\Big)
\ensuremath{\mathsf E}[  \ensuremath{{\mathsf{F}}}(\ensuremath x^\ensuremath n) - \ensuremath{{\mathsf{F}}}(\ensuremath x^{\ensuremath n+1}) \,\vert\, \ensuremath{\mathfrak{E}}_{\ensuremath n-1} ].
\end{align}
\end{enumerate}
\end{proposition}
\begin{proof}
Let $\ensuremath n \in \ensuremath \mathbb{N}$ and $\ensuremath \xx \in \dom \ensuremath{{\mathsf{F}}}$.
Since
\begin{equation*}
\norm{\ensuremath x^\ensuremath n - \ensuremath \xx}_{\ensuremath{{\Gammas}}^{-1}}^2 - \norm{\bar{\ensuremath x}^{\ensuremath n+1} - \ensuremath \xx}_{\ensuremath{{\Gammas}}^{-1}}^2
= - \norm{ \ensuremath x^\ensuremath n - \bar{\ensuremath x}^{\ensuremath n+1}}_{\ensuremath{{\Gammas}}^{-1}}^2
+ 2 \scalarp{\ensuremath x^{\ensuremath n} - \bar{\ensuremath x}^{\ensuremath n+1}, \ensuremath x^\ensuremath n - \ensuremath \xx}_{\ensuremath{{\Gammas}}^{-1}},
\end{equation*}
we derive from \eqref{eq:20190610a}, multiplied by $2$, that
\begin{align}
\nonumber(1 + \sigma_{\ensuremath \mathsf{\Gamma}^{-1}}) \norm{\bar{\ensuremath x}^{\ensuremath n+1} - \ensuremath \xx}_{\ensuremath{{\Gammas}}^{-1}}^2
&\leq (1 + \sigma_{\ensuremath \mathsf{\Gamma}^{-1}})  \norm{\ensuremath x^{\ensuremath n} - \ensuremath \xx}_{\ensuremath{{\Gammas}}^{-1}}^2
+ (\delta - 1)  \norm{\bar{\ensuremath x}^{\ensuremath n+1} - \ensuremath x^{\ensuremath n}}_{\ensuremath{{\Gammas}}^{-1}}^2 \\[1ex]
\nonumber&\quad+\frac{2}{\ensuremath \mathsf{p}_{\min}}
\ensuremath{\mathsf E}[  \ensuremath{{\mathsf{F}}}(\ensuremath x^\ensuremath n) - \ensuremath{{\mathsf{F}}}(\ensuremath x^{\ensuremath n+1}) \,\vert\, \ensuremath{\mathfrak{E}}_{\ensuremath n-1} ]  \\[1ex]
\label{eq:20181124a}&\quad - 2 \bigg( \ensuremath{{\mathsf{F}}}(\ensuremath x^\ensuremath n) - \ensuremath{{\mathsf{F}}}(\ensuremath \xx)
+ \frac{\mu_{\ensuremath \mathsf{\Gamma}^{-1}} + \sigma_{\ensuremath \mathsf{\Gamma}^{-1}}}{2} \norm{\ensuremath x^\ensuremath n - \ensuremath \xx}_{\ensuremath{{\Gammas}}^{-1}}^2 \bigg).
\end{align}
Then for an $\ensuremath {\mathsf{H}}$-valued $\ensuremath{\mathfrak{E}}_{\ensuremath n-1}$-measurable random variable $\ensuremath x$, Proposition~\ref{p:20190313c} yields
\begin{align}
\nonumber
(1 + \sigma_{\ensuremath \mathsf{\Gamma}^{-1}})&\ensuremath{\mathsf E}[\norm{\ensuremath x^{\ensuremath n+1} - \ensuremath x}_\ensuremath{{\WW}}^2 \,\vert\, \ensuremath{\mathfrak{E}}_{\ensuremath n-1} ] \\[1ex]
\nonumber&\leq (1 + \sigma_{\ensuremath \mathsf{\Gamma}^{-1}})  \norm{\ensuremath x^{\ensuremath n} - \ensuremath x}_{\ensuremath{{\WW}}}^2
+ (\delta - 1) \ensuremath{\mathsf E}[ \norm{\ensuremath x^{\ensuremath n+1} - \ensuremath x^{\ensuremath n}}_\ensuremath{{\WW}}^2
\,\vert\, \ensuremath{\mathfrak{E}}_{\ensuremath n-1} ] \\[1ex]
\nonumber&\qquad\qquad+\frac{2}{\ensuremath \mathsf{p}_{\min}}
\ensuremath{\mathsf E}[  \ensuremath{{\mathsf{F}}}(\ensuremath x^\ensuremath n) - \ensuremath{{\mathsf{F}}}(\ensuremath x^{\ensuremath n+1}) \,\vert\, \ensuremath{\mathfrak{E}}_{\ensuremath n-1} ]  \\[1ex]
\label{eq:20180925h}&\qquad\qquad - 2 \bigg( \ensuremath{{\mathsf{F}}}(\ensuremath x^\ensuremath n) - \ensuremath{{\mathsf{F}}}(\ensuremath x)
+ \frac{\mu_{\ensuremath \mathsf{\Gamma}^{-1}} + \sigma_{\ensuremath \mathsf{\Gamma}^{-1}}}{2} \norm{\ensuremath x^\ensuremath n - \ensuremath x}_{\ensuremath{{\Gammas}}^{-1}}^2 \bigg).
\end{align}
Taking $\ensuremath x=\ensuremath x^\ensuremath n$ in \eqref{eq:20180925h}, we have
\begin{equation}
\label{eq:20180925g}
\frac{\ensuremath \mathsf{p}_{\min}}{2}(2 +\sigma_{\ensuremath \mathsf{\Gamma}^{-1}} - \delta)\ensuremath{\mathsf E}[\norm{\ensuremath x^{\ensuremath n+1} - \ensuremath x^\ensuremath n}_\ensuremath{{\WW}}^2
\,\vert\, \ensuremath{\mathfrak{E}}_{\ensuremath n-1} ]
\leq \ensuremath{\mathsf E}\big[  \ensuremath{{\mathsf{F}}}(\ensuremath x^\ensuremath n) - \ensuremath{{\mathsf{F}}}(\ensuremath x^{\ensuremath n+1})
\,\vert\, \ensuremath{\mathfrak{E}}_{\ensuremath n-1} \big],
\end{equation}
which plugged into \eqref{eq:20180925h},
with $\ensuremath x\equiv \ensuremath \xx \in \dom \ensuremath{{\mathsf{F}}}$, gives \ref{p:20180927a_ii}.
Moreover, taking the expectation in  \eqref{eq:20180925g}, we obtain
\begin{equation}
\label{eq:20190128a}
\frac{\ensuremath \mathsf{p}_{\min}}{2}(2 + \sigma_{\ensuremath \mathsf{\Gamma}^{-1}} - \delta)\ensuremath{\mathsf E}\big[\norm{\ensuremath x^{\ensuremath n+1} - \ensuremath x^\ensuremath n}_\ensuremath{{\WW}}^2 \big]
\leq \ensuremath{\mathsf E}[  \ensuremath{{\mathsf{F}}}(\ensuremath x^\ensuremath n)]  -  \ensuremath{\mathsf E}[\ensuremath{{\mathsf{F}}}(\ensuremath x^{\ensuremath n+1})],
\end{equation}
which gives \ref{p:20180927a_0}.
Finally, set for all $\ensuremath n \in \ensuremath \mathbb{N}$,
 $\xi_\ensuremath n = \ensuremath{\mathsf E}\big[ \ensuremath{{\mathsf{F}}}(\ensuremath x^\ensuremath n) - \ensuremath{{\mathsf{F}}}(\ensuremath x^{\ensuremath n+1}) \big\vert
\ensuremath{\mathfrak{E}}_{\ensuremath n-1} \big] \geq 0$. Then
\begin{equation*}
\ensuremath{\mathsf E} \bigg[\sum_{\ensuremath n=0}^{+\infty} \xi_\ensuremath n \bigg]
= \sum_{\ensuremath n=0}^{+\infty} \ensuremath{\mathsf E}[\xi_\ensuremath n] = \sum_{\ensuremath n=0}^{+\infty}
\ensuremath{\mathsf E}[\ensuremath{{\mathsf{F}}}(\ensuremath x^\ensuremath n)] - \ensuremath{\mathsf E}[\ensuremath{{\mathsf{F}}}(\ensuremath x^{\ensuremath n+1})] \leq \ensuremath{\mathsf E}[\ensuremath{{\mathsf{F}}}(\ensuremath x^0)] - \inf_{\ensuremath n \in \ensuremath \mathbb{N}} \ensuremath{\mathsf E}[\ensuremath{{\mathsf{F}}}(\ensuremath x^\ensuremath n)].
\end{equation*}
This shows that if $\inf_{\ensuremath n \in \ensuremath \mathbb{N}} \ensuremath{\mathsf E}[\ensuremath{{\mathsf{F}}}(\ensuremath x^\ensuremath n)]>0$, then
$\sum_{\ensuremath n=0}^{+\infty} \xi_\ensuremath n$ is $\ensuremath{\mathsf P}$-integrable and hence
it is $\ensuremath{\mathsf P}$-a.s.~finite.
Then \ref{p:20180927a_i} follows from \eqref{eq:20180925g} and Proposition~\ref{p:20190313c}.
\end{proof}

\begin{proposition}
\label{p:20181219c}
Under the same assumptions of Proposition~\ref{p:20180927a},
suppose that
condition \ref{eq:S1} is replaced by condition \ref{eq:S2}.
Then
\begin{equation*}
(\forall\, \ensuremath n \in \ensuremath \mathbb{N})\qquad\frac{2 + \sigma_{\ensuremath \mathsf{\Gamma}^{-1}} - \delta}{2}
\norm{\ensuremath x^{\ensuremath n+1} - \ensuremath x^\ensuremath n}_{\ensuremath{{\Gammas}}^{-1}}^2 \leq
 \ensuremath{{\mathsf{F}}}(\ensuremath x^{\ensuremath n}) - \ensuremath{{\mathsf{F}}}(\ensuremath x^{\ensuremath n+1})
 \quad \ensuremath{\mathsf P} \text{ a.s.}
\end{equation*}
\end{proposition}
\begin{proof}
We derive from \ref{eq:S2} (since $\ensuremath \varepsilon^n$ has the same distribution of
$\ensuremath \varepsilon$) and \eqref{eq:20181114c} that
\begin{equation}
\label{eq:20181219b}
\scalarp{\nabla \ensuremath{{\mathsf{f}}}(\ensuremath x^\ensuremath n), \ensuremath x^{\ensuremath n+1} - \ensuremath x^{\ensuremath n}}
 \leq \ensuremath{{\mathsf{f}}}(\ensuremath x^{\ensuremath n}) -  \ensuremath{{\mathsf{f}}}(\ensuremath x^{\ensuremath n+1})
+ \sum_{i=1}^m
\frac 1 2 \varepsilon^\ensuremath n_i \nu_i \norm{\Delta^\ensuremath n_{i}}^2  \quad \ensuremath{\mathsf P} \text{ a.s.}
\end{equation}
Therefore, summing \eqref{eq:20181108b}, from $i=1$ to $m$, we have
\begin{align*}
\frac{2 + \sigma_{\ensuremath \mathsf{\Gamma}^{-1}}}{2}\sum_{i=1}^m \varepsilon_i^\ensuremath n
 \frac{1}{\gamma_i}\norm{\Delta_i^\ensuremath n}^2 &\leq
 \ensuremath{{\gs}}(\ensuremath x^{\ensuremath n}) - \ensuremath{{\gs}}(\ensuremath x^{\ensuremath n+1}) + \scalarp{\nabla \ensuremath{{\mathsf{f}}}(\ensuremath x^\ensuremath n), \ensuremath x^\ensuremath n - \ensuremath \xx^{\ensuremath n+1}}\\
 & \leq \ensuremath{{\gs}}(\ensuremath x^{\ensuremath n}) - \ensuremath{{\gs}}(\ensuremath x^{\ensuremath n+1}) +  \ensuremath{{\mathsf{f}}}(\ensuremath x^{\ensuremath n}) -  \ensuremath{{\mathsf{f}}}(\ensuremath x^{\ensuremath n+1})
 + \frac 1 2 \sum_{i=1}^m
\varepsilon^\ensuremath n_i  \nu_i \norm{\Delta^\ensuremath n_{i}}^2\quad \ensuremath{\mathsf P} \text{ a.s.}
\end{align*}
Hence
$(1/2) \sum_{i=1}^m
(2 + \sigma_{\ensuremath \mathsf{\Gamma}^{-1}} - \gamma_i \nu_i )\gamma_i^{-1}\varepsilon_i^\ensuremath n\norm{\Delta_i^\ensuremath n}^2
\leq  \ensuremath{{\mathsf{F}}}(\ensuremath x^{\ensuremath n}) - \ensuremath{{\mathsf{F}}}(\ensuremath x^{\ensuremath n+1})$ $\ensuremath{\mathsf P}$-a.s.
\end{proof}

\begin{proposition}
\label{prop:20180927a}
Under the assumptions of Proposition~\ref{p:20180927a},
suppose in addition that
  $\ensuremath{{\mathsf{F}}}$ is bounded from below.
Then, there exist $(\ensuremath y^\ensuremath n)_{\ensuremath n \in \ensuremath \mathbb{N}}$ and $(\ensuremath v^\ensuremath n)_{\ensuremath n \in \ensuremath \mathbb{N}}$,
sequences of $\ensuremath {\mathsf{H}}$-valued random variables, such that
the following hold.
\begin{enumerate}[{\rm (i)}]
\item\label{eq:20190113a_i}
$\ensuremath v^\ensuremath n\in \partial \ensuremath{{\mathsf{F}}}(\ensuremath y^{\ensuremath n})$ $\ensuremath{\mathsf P}$-a.s.
\item\label{eq:20190113a_ii}  $\ensuremath y^{\ensuremath n} - \ensuremath x^{\ensuremath n} \to 0$
and $\ensuremath v^{\ensuremath n} \to 0$ $\ensuremath{\mathsf P}$-a.s.
\end{enumerate}
\end{proposition}
\begin{proof}
It follows from \eqref{eq:20171206a} that,
$(x_i^\ensuremath n(\omega) - \bar{x}_i^{\ensuremath n+1}(\omega))/\gamma_{i} - \nabla_i \ensuremath{{\mathsf{f}}} (\ensuremath x^\ensuremath n(\omega))
\in \partial \ensuremath{{\mathsf{h}}}_i(\bar{x}_i^{\ensuremath n+1}(\omega))$,
for all $i \in [m]$ and $\omega \in \Omega$. Hence
\begin{equation*}
\Big( \frac{x_i^\ensuremath n(\omega) - \bar{x}_i^{\ensuremath n+1}(\omega)}{\gamma_i} \Big)_{1 \leq i \leq m} - \nabla \ensuremath{{\mathsf{f}}}(\ensuremath x^\ensuremath n(\omega))
\in \partial \ensuremath{{\gs}}(\bar{\ensuremath x}^{\ensuremath n+1}).
\end{equation*}
 Set $\ensuremath y^\ensuremath n = \bar{\ensuremath x}^{\ensuremath n+1}$ and let
 $\ensuremath v^\ensuremath n \colon \Omega \to \ensuremath {\mathsf{H}}$ be such that, for every $\omega \in \Omega$,
\begin{multline*}
\ensuremath v^\ensuremath n(\omega) = \Big(\frac{x_i^\ensuremath n(\omega) - y_i^\ensuremath n(\omega)}{\gamma_i} \Big)_{1 \leq i \leq m}
 + \nabla \ensuremath{{\mathsf{f}}} (\ensuremath y^\ensuremath n(\omega)) - \nabla \ensuremath{{\mathsf{f}}} (\ensuremath x^\ensuremath n(\omega))  \\
 \in \partial \ensuremath{{\gs}}(\ensuremath y^{\ensuremath n}(\omega)) + \nabla \ensuremath{{\mathsf{f}}}(\ensuremath y^\ensuremath n(\omega)) = \partial \ensuremath{{\mathsf{F}}}(\ensuremath y^\ensuremath n(\omega)).
\end{multline*}
Clearly $\ensuremath v^\ensuremath n$ is measurable and hence it is a random variable.
Moreover,
for every $\omega \in \Omega$,
\begin{align*}
\norm{\ensuremath v^\ensuremath n(\omega)} \leq \frac{1}{\gamma_{\min}} \norm{\ensuremath x^\ensuremath n(\omega) - \ensuremath y^\ensuremath n(\omega)}
+ \norm{\nabla \ensuremath{{\mathsf{f}}}(\ensuremath y^\ensuremath n(\omega)) - \nabla \ensuremath{{\mathsf{f}}}(\ensuremath x^\ensuremath n(\omega))}.
\end{align*}
Now, since $\ensuremath{{\mathsf{F}}}$ is bounded from below, Proposition~\ref{p:20180927a}\ref{p:20180927a_i} yields that
$(\norm{\ensuremath y^\ensuremath n - \ensuremath x^\ensuremath n}_{\ensuremath \mathsf{\Gamma}^{-1}}^2)_{\ensuremath n \in \ensuremath \mathbb{N}}$ is summable $\ensuremath{\mathsf P}$-a.s.~and hence
$\ensuremath y^\ensuremath n - \ensuremath x^\ensuremath n\to 0$ $\ensuremath{\mathsf P}$-a.s. The statement follows from
the fact  that $\nabla \ensuremath{{\mathsf{f}}}$ is Lipschitz continuous
(see Theorem~\ref{thm:stepsizes}\ref{thm:stepsizes_iv}).
\end{proof}

Now we are ready to state one of the main convergence results of this paper.
From one hand, it
 extends to the stochastic setting
a well-known convergence rate of the (deterministic) forward-backward algorithm
\cite{Dav16,Gar17,Sal17}. On the other hand, it proves the almost sure weak convergence of the iterates of Algorithm~\ref{algoRCD} in the convex case.
We stress that none of the works \cite{LuX15,Nec16,Qu16b,Ric15,Ric16,Ric16b,Tap18} addresses this latter aspect.
To the best of our knowledge, \cite{Com15} is the only work that
proves almost sure weak convergence of the iterates.
However, in \cite[Corollary~5.11]{Com15} the stepsize is set according to the (global)
Lipschitz constant of $\nabla \ensuremath{{\mathsf{f}}}$ which,
in general, leads to smaller stepsizes and worse upper bounds on convergence rates.
See the subsequent discussion.

\begin{theorem}
\label{thm:20171207a}
 Let \ref{eq:A1}--\ref{eq:A3} be satisfied.
Let $(\nu_{i})_{1 \leq i \leq m} \in \ensuremath \mathbb{R}^m_{++}$ and suppose that
 \ref{eq:S1} holds.
Let $(\ensuremath x^\ensuremath n)_{\ensuremath n \in \ensuremath \mathbb{N}}$ be generated by  Algorithm~\ref{algoRCD}
with, for every $i \in [m]$, $\gamma_i < 2/\nu_i$.
Set $\delta = \max_{1 \leq i \leq m} \gamma_i  \nu_i$
 and
$\ensuremath \mathsf{p}_{\min} = \min_{1 \leq i \leq m} \ensuremath \mathsf{p}_i$.
Let $\ensuremath{{\WW}}$ be as in \eqref{eq:W} and
set  $\ensuremath{{\mathsf{F}}}= \ensuremath{{\mathsf{f}}}+\ensuremath{{\gs}}$, $\ensuremath{{\mathsf{F}}}_* = \inf \ensuremath{{\mathsf{F}}}$, and $\ensuremath {\mathsf{S}}_*
= \argmin \ensuremath{{\mathsf{F}}} \subset \ensuremath {\mathsf{H}}$.
Then, the following hold.
\begin{enumerate}[{\rm (i)}]
\item\label{thm:20171207a_0} $\ensuremath{\mathsf E}[\ensuremath{{\mathsf{F}}}(\ensuremath x^{\ensuremath n})]\to \ensuremath{{\mathsf{F}}}_*$.
\item\label{thm:20171207a_i} Suppose that $\ensuremath {\mathsf{S}}_*\neq \varnothing$.
Then $\ensuremath{\mathsf E}[\ensuremath{{\mathsf{F}}}(\ensuremath x^{\ensuremath n})] -  \ensuremath{{\mathsf{F}}}_* = o(1/\ensuremath n)$ and, for every integer $\ensuremath n \geq 1$,
\begin{equation}
\label{eq:20171207a}
\ensuremath{\mathsf E}[\ensuremath{{\mathsf{F}}}(\ensuremath x^{\ensuremath n})] - \ensuremath{{\mathsf{F}}}_*  \leq
\bigg[ \frac{ \mathrm{dist}_\ensuremath{{\WW}}^2(\ensuremath \xx^0, \ensuremath {\mathsf{S}}_*) }{2}
+ \Big(\frac{\max\big\{1, (2-\delta)^{-1}\big\}}{\ensuremath \mathsf{p}_{\min}} -1\Big)
 (\ensuremath{{\mathsf{F}}}(\ensuremath \xx^{0})  - \ensuremath{{\mathsf{F}}}_*)\bigg] \frac{1}{n}.
\end{equation}
Moreover,
there exists a random variable $\ensuremath x_*$ taking values in $\ensuremath {\mathsf{S}}_*$ such that
$\ensuremath x^\ensuremath n \rightharpoonup \ensuremath x_*$ $\ensuremath{\mathsf P}$-a.s.
\end{enumerate}
\end{theorem}
\begin{proof}
Proposition~\ref{p:20180927a}\ref{p:20180927a_ii}
with $\mu_{\ensuremath \mathsf{\Gamma}^{-1}}=\sigma_{\ensuremath \mathsf{\Gamma}^{-1}}=0$ gives,
for all $\ensuremath \xx \in \dom \ensuremath{{\mathsf{F}}}$ and $\ensuremath n \in \ensuremath \mathbb{N}$,
\begin{equation*}
\ensuremath{\mathsf E}[\norm{\ensuremath x^{\ensuremath n+1} - \ensuremath \xx}_\ensuremath{{\WW}}^2 \,\vert\, \ensuremath{\mathfrak{E}}_{\ensuremath n-1} ]
\leq   \norm{\ensuremath x^{\ensuremath n} - \ensuremath \xx}_{\ensuremath{{\WW}}}^2  + 2\ensuremath{\mathsf E}[ \ensuremath{{\mathsf{F}}}(\ensuremath \xx) -\ensuremath{{\mathsf{F}}}(\ensuremath x^{\ensuremath n+1})\,\vert\,\ensuremath{\mathfrak{E}}_{\ensuremath n-1}] + \xi_\ensuremath n,
\end{equation*}
where
\begin{equation*}
\xi_\ensuremath n = b_1  \ensuremath{\mathsf E}[\ensuremath{{\mathsf{F}}}(\ensuremath x^\ensuremath n) - \ensuremath{{\mathsf{F}}}(\ensuremath x^{\ensuremath n+1})\,\vert\,
\ensuremath{\mathfrak{E}}_{\ensuremath n-1}],\qquad b_1 = 2  \bigg(\frac{\max\{1,1/(2-\delta)\}}{\ensuremath \mathsf{p}_{\min}} - 1\bigg).
\end{equation*}
Note that the random variables $\ensuremath x^\ensuremath n$'s are discrete with finite range
and $(\ensuremath{\mathsf E}[\ensuremath{{\mathsf{F}}}(\ensuremath x^\ensuremath n)])_{\ensuremath n \in \ensuremath \mathbb{N}}$ is decreasing.
Moreover,
$\sum_{\ensuremath n \in \ensuremath \mathbb{N}} \ensuremath{\mathsf E}[\xi_\ensuremath n] \leq b_1 (\ensuremath{{\mathsf{F}}}(\ensuremath \xx^0) - \inf_{\ensuremath n \in \ensuremath \mathbb{N}} \ensuremath{\mathsf E}[\ensuremath{{\mathsf{F}}}(\ensuremath x^\ensuremath n)])$.
Therefore,
 the statement follows from
 Theorem~\ref{p:SMFejer} and
Proposition~\ref{prop:20180927a}.
\end{proof}

\paragraph{Discussion.}
In the following we examine some crucial aspects related to Algorithm~\ref{algoRCD}.
We suppose that \ref{eq:B0} and \ref{eq:B1} hold and that
for every $i \in [m]$, $\gamma_i = \delta/\nu_i$ with $\delta \in \left]0,2\right[$.
\setlist[description]{font=\normalfont\itshape}
\begin{description}[leftmargin=0cm]
\item[The benefit of a parallel block update.]
Here we discuss the advantage of updating multiple blocks in parallel instead of just a single block.
We consider the setting of a $\tau$-nice uniform block sampling, which was described in
Section~\ref{subsec:beta}.
In this case, for ever $i,j \in [m]$, with $i\neq j$,
$\ensuremath \mathsf{p}_i = \ensuremath \mathsf{p}_j$ and, since $\tau = \sum_{i=1}^m \ensuremath{\mathsf E}[\varepsilon_i]
= \sum_{i=1}^m \ensuremath \mathsf{p}_i$,  we have, for every $i \in [m]$, $\ensuremath \mathsf{p}_i = \ensuremath \mathsf{p}: = \tau/m$.
Moreover, we can set for every $i \in [m]$,
$\nu_i = \beta_1 L_i$ with $\beta_1$ defined as in \eqref{eq:20200316c}.
In order to compare different choices of $\tau$, we normalize the iterations so to match
the same computational cost per iteration of the standard (full parallel) forward-backward algorithm (FB).
It follows from \eqref{eq:20171207a} that after $\ensuremath n_\tau =\lceil m \bar{\ensuremath n}/\tau \rceil$ iterations of Algorithm~\ref{algoRCD},
which have the same total computational cost of  $\bar{\ensuremath n}$ iterations of FB, we have
\begin{equation}
\label{eq:20200330a}
\ensuremath{\mathsf E}[\ensuremath{{\mathsf{F}}}(\ensuremath x^{\lceil m \bar{\ensuremath n}/\tau \rceil})] -  \ensuremath{{\mathsf{F}}}_*
\leq \bigg[ \beta_1\frac{\mathrm{dist}_{\ensuremath \mathsf{\Lambda}}^2(\ensuremath \xx^0, \ensuremath {\mathsf{S}}_*)}{2\delta}
+ \Big(\max\Big\{1,\frac{1}{2-\delta}\Big\} -\ensuremath \mathsf{p}\Big)  ( \ensuremath{{\mathsf{F}}}(\ensuremath \xx^{0}) - \ensuremath{{\mathsf{F}}}_*)\bigg] \frac{1}{\bar{\ensuremath n}},
\end{equation}
where $\ensuremath \mathsf{\Lambda} = \bigoplus_{i=1}^m L_i \ensuremath{\mathsf{Id}}_i$.
Now, since, $\beta_1 = 1 + (\tau-1)(\eta-1)/(m-1)$, we see that
if $\eta \ll m$ and $\tau\ll m$, then $\beta_1$
is close to 1 (and $\ensuremath \mathsf{p}$ is close to zero), so that $\beta_1$ nearly does not depend on $\tau$, as long as $\tau$
 remains sufficiently small.
For instance, in Section~\ref{sec:experiments} we consider the setting where $m = 10^5$ and
$\eta = 148$. In such case, if we let $\tau=1, 5, 10, 50$, the corresponding $\beta_1$'s
and $\ensuremath \mathsf{p}$'s are essentially the same so that the right hand side of \eqref{eq:20200330a}
does not change much. Therefore, the above options for $\tau$
require the same total amount of computations (i.e., $m \bar{n}$ block-coordinate updates) and lead essentially
to the same improvement in the objective function.
However, a parallel implementation,
say with $\tau = 50$ on a CPU with $50$ cores,
will be $50$ times faster than a serial implementation ($\tau=1$)
which uses only one core per iteration.
In summary, in the large scale ($m$ large) and sparse ($\eta\ll m$) setting,
the parallel strategy ($\tau>1$ and $\tau$ equal to the number of CPU cores)
is definitely advantageous provided that $\tau$
is sufficiently small compared to $m$.
\item[Comparison with \cite{Com15}.]
The almost sure weak convergence of the iterates of Algorithm~\ref{algoRCD} is also obtained in \cite{Com15},
but with stepsizes set according to the global Lipschitz constant of the gradient of $\ensuremath{{\mathsf{f}}}$.
Let $L$ be the Lipschitz constant of $\nabla \ensuremath{{\mathsf{f}}}$ and note that
$\ensuremath{{\mathsf{f}}}$ is also Lipschitz smooth in the norm $\norm{\cdot}_{\ensuremath \mathsf{\Lambda}}$,
defined by the operator
$\ensuremath \mathsf{\Lambda} = \bigoplus_{i=1}^m L_i \ensuremath{\mathsf{Id}}_i$, with constant $\eta$ (see Corollary~\ref{p:20181014a}\ref{p:20181014a_3}).
Therefore, the results in \cite{Com15} can be applied in the original norm $\norm{\cdot}$ or in the
norm $\norm{\cdot}_{\ensuremath \mathsf{\Lambda}}$.
In this respect we note that
since
\begin{equation*}
\ensuremath{\text{\sf prox}}_{\alpha \ensuremath{{\gs}}}^\ensuremath \mathsf{\Lambda}(\ensuremath x^\ensuremath n - \alpha \nabla^\ensuremath \mathsf{\Lambda} \ensuremath{{\mathsf{f}}}(\ensuremath x^\ensuremath n)) = (\ensuremath{\text{\sf prox}}_{(\alpha/L_i) \ensuremath{{\mathsf{h}}}_i} (x_i^\ensuremath n - (\alpha/L_i) \nabla_i \ensuremath{{\mathsf{f}}}(\ensuremath x^\ensuremath n)))_{1 \leq i \leq m}
\qquad (\alpha=\delta/\eta, 0<\delta<2),
\end{equation*}
the implementation in the norm $\norm{\cdot}_{\ensuremath \mathsf{\Lambda}}$ is nothing but Algorithm~\ref{algoRCD}
with stepsizes $\gamma_i = \delta/(\eta L_i)$.
In both cases Corollary~5.11 in \cite{Com15}, applied in the corresponding norms,
proves weak convergence of the iterates for
Algorithm~\ref{algoRCD} with stepsizes $\gamma_i \equiv \delta/L$
and $\gamma_i = \delta/(\eta L_i)$ respectively.
However, Theorem~\ref{thm:20171207a}, together with Theorem~\ref{thm:stepsizes},
allows to set the stepsizes as $\gamma_i = \delta/(\beta_{1,i} L_i)$.
Since $\beta_{1,i} L_i$ may be much smaller than $L$ and, in view of Remark~\ref{rmk:20181219ee},
it is always smaller than $\eta L_i$, Theorem~\ref{thm:20171207a}
provides a significant improvement over \cite{Com15} in terms of flexibility in the stepsizes.
\item[The advantage over the standard FB.]
We consider the forward-backward algorithm (FB)
in the original norm of $\ensuremath \mathsf{H}$ and in the norm $\norm{\cdot}_{\ensuremath \mathsf{\Lambda}}$.
All the remarks about the stepsizes discussed in the previous paragraph apply also here.
Moreover,
in the case $\delta = 1$, standard convergence rate for FB (see e.g., \cite{Dav16}) yields that after $\bar{\ensuremath n}$ iterations, we have
\begin{equation}
\label{eq:20200330b}
\ensuremath{{\mathsf{F}}}(\ensuremath \mathsf{x}^{\bar{\ensuremath n}}) - \ensuremath{{\mathsf{F}}}_*
\leq \frac{L}{2} \frac{\mathrm{dist}^2(\ensuremath \mathsf{x}^0, \ensuremath {\mathsf{S}}_*)}{\bar{\ensuremath n}}
\quad\text{or}\quad
 \ensuremath{{\mathsf{F}}}(\ensuremath \mathsf{x}^{\bar{\ensuremath n}}) - \ensuremath{{\mathsf{F}}}_*
\leq \frac{\eta}{2} \frac{\mathrm{dist}_{\ensuremath \mathsf{\Lambda}}^2(\ensuremath \mathsf{x}^0, \ensuremath {\mathsf{S}}_*)}{\bar{\ensuremath n}},
\end{equation}
depending on which of the two above implementations of FB we consider.
In order to appropriately compare the rates \eqref{eq:20200330b} with that
of Algorithm~\ref{algoRCD} given in Theorem~\ref{thm:20171207a}
in the following we set $\delta=1$ and analyze two choices of the block sampling.
\begin{enumerate}[(i)]
\item Assume that we perform a $\tau$-nice block sampling.
Then we saw that the (normalized) convergence rate of Algorithm~\ref{algoRCD} is
\eqref{eq:20200330a}.
We first note that \eqref{eq:20200330a} reduces to
the second inequality in \eqref{eq:20200330b} when $\delta=1$ and $\tau = m$.
Comparing the bounds in \eqref{eq:20200330b} with
\eqref{eq:20200330a} (with $\delta=1$) we see that, if we assume that the terms
$\mathrm{dist}_{\ensuremath \mathsf{\Lambda}}^2(\ensuremath \mathsf{x}^0, \ensuremath {\mathsf{S}}_*)$ and $\ensuremath{{\mathsf{F}}}(\ensuremath \xx^{0}) - \ensuremath{{\mathsf{F}}}_*$ are about of the same magnitude, then
Algorithm~\ref{algoRCD} features always a better rate than FB if
implemented in the norm $\norm{\cdot}_{\ensuremath \mathsf{\Lambda}}$ (since $\beta_1 \leq \eta$), whereas if FB
is implemented in the original norm of $\ensuremath {\mathsf{H}}$, Algorithm~\ref{algoRCD} is still a better choice provided that $\beta_1\max_{1 \leq i \leq m} L_i \leq L$.
\item Suppose that the block sampling performs on average $\tau$
updates per iteration and that, for every $i \in [m]$,  $\ensuremath \mathsf{p}_i$ is
proportional to the Lipschitz constant $L_i$, that is,
$\ensuremath \mathsf{p}_i = \tau L_i/(\sum_{j=1}^m L_j)$ (provided that $\tau \leq (\sum_{j=1}^m L_j)/\max_{1 \leq j \leq m} L_j$).
In this case,
as stated at the beginning of Section~\ref{subsec:beta},  we can let
  $\nu_i = \beta_2 L_i$ and $\gamma_i = 1/\nu_i$.
Then,
\eqref{eq:20171207a}
becomes for $\delta=1$ and $\ensuremath n =\lceil m \bar{\ensuremath n}/\tau \rceil$,
\begin{equation*}
\ensuremath{\mathsf E}[\ensuremath{{\mathsf{F}}}(\ensuremath x^{\lceil m \bar{\ensuremath n}/\tau \rceil})] - \ensuremath{{\mathsf{F}}}_*  \leq \bigg[ \beta_2
\bar{L}\frac{ \mathrm{dist}^2(\ensuremath \xx^0, \ensuremath {\mathsf{S}}_*) }{2}
+ \bigg( \frac{\bar{L}}{L_{\mathrm{min}}} - \frac{\tau}{m}\bigg)
(\ensuremath{{\mathsf{F}}}(\ensuremath \xx^{0}) - \ensuremath{{\mathsf{F}}}_*)\bigg] \frac{1}{\bar{n}},
\end{equation*}
where  $\bar{L} = \sum_{i=1}^m L_i/m$ and $L_{\min} = \min_{1 \leq i \leq m} L_i$.
Here we see that Algorithm~\ref{algoRCD} can be superior to FB if
$(\beta_2 +2/L_{\min})\bar{L} \leq L$.
\end{enumerate}
\end{description}

We now provide an additional convergence theorem, including the strongly convex case,
which extends \cite[Theorem~1]{LuX15} and
\cite[Theorem~3]{Tap18} to an arbitrary (not necessarily uniform) sampling and
to the more general stepsize rule \eqref{eq:stepsizerule}.
The proof is still based on Proposition~\ref{p:20180927a}\ref{p:20180927a_ii},
but will closely follow that of the above cited works,
so it will be postponed in Appendix~\ref{sec:appB}.

\begin{theorem}
\label{p:20181130a}
Under the same assumptions of Theorem~\ref{thm:20171207a},
let $\mu_{\ensuremath \mathsf{\Gamma}^{-1}}$ and $\sigma_{\ensuremath \mathsf{\Gamma}^{-1}}$ be the moduli of strong convexity of
$\ensuremath{{\mathsf{f}}}$ and $\ensuremath{{\gs}}$ respectively, in the norm $\norm{\cdot}_{\ensuremath{{\Gammas}}^{-1}}$,
and suppose that $\ensuremath {\mathsf{S}}_* \neq \varnothing$.
Then, the following hold.
\begin{enumerate}[{\rm (i)}]
\item\label{p:20181130a_ii}
For every integer $n \geq 1$,
\begin{equation*}
\ensuremath{\mathsf E} [\ensuremath{{\mathsf{F}}}(\ensuremath x^{\ensuremath n})] - \ensuremath{{\mathsf{F}}}_*
\leq \frac{1}{1 + \ensuremath \mathsf{p}_{\min} n}
\bigg( \ensuremath \mathsf{p}_{\min}\frac{\mathrm{dist}_\ensuremath{{\WW}}^2(\ensuremath \xx^0, \ensuremath {\mathsf{S}}_*)}{2}
+ \max\Big\{1, \frac{1}{2-\delta}\Big\}(\ensuremath{{\mathsf{F}}}(\ensuremath \xx^{0}) - \ensuremath{{\mathsf{F}}}_*) \bigg).
\end{equation*}
\item\label{p:20181130a_iii} Suppose
 that $\mu_{\ensuremath \mathsf{\Gamma}^{-1}}+\sigma_{\ensuremath \mathsf{\Gamma}^{-1}}>0$ and $\delta \leq 1$.
 Set $\rho = (1 - 2 \ensuremath \mathsf{p}_{\min} (\mu_{\ensuremath \mathsf{\Gamma}^{-1}} + \sigma_{\ensuremath \mathsf{\Gamma}^{-1}})/(1 + \mu_{\ensuremath \mathsf{\Gamma}^{-1}} + 2 \sigma_{\ensuremath \mathsf{\Gamma}^{-1}}))$.
 Then,
for every $n \in \ensuremath \mathbb{N}$,
\begin{equation*}
\ensuremath{\mathsf E}[\ensuremath{{\mathsf{F}}}(\ensuremath x^{\ensuremath n})] - \ensuremath{{\mathsf{F}}}_*
\leq
\rho^\ensuremath n
\bigg( \ensuremath \mathsf{p}_{\min} (1+\sigma_{\ensuremath \mathsf{\Gamma}^{-1}}) \frac{\mathrm{dist}_\ensuremath{{\WW}}^2(\ensuremath \xx^0, \ensuremath {\mathsf{S}}_*)}{2}
+ \ensuremath{{\mathsf{F}}}(\ensuremath \xx^0) - \ensuremath{{\mathsf{F}}}_* \bigg).
\end{equation*}
\end{enumerate}
\end{theorem}

\begin{remark}
\label{rmk:20190130a}
Let $\mu$ and $\sigma$ be the moduli of strong convexity
of $\ensuremath{{\mathsf{f}}}$ and $\ensuremath{{\gs}}$ respectively, in the original norm $\norm{\cdot}$.
Let, for every $i \in [m]$, $\gamma_i = 1/\nu_i$ and set $\nu_{\max} = \max_{1 \leq i \leq m} \nu_i$.
Then $\mu/\nu_{\max}$ and $\sigma/\nu_{\max}$
are the moduli of strong convexity of $\ensuremath{{\mathsf{f}}}$ and $\ensuremath{{\gs}}$ respectively, w.r.t.~$\norm{\cdot}_{\ensuremath \mathsf{\Gamma}^{-1}}$.
Therefore the rate in
Theorem~\ref{p:20181130a}\ref{p:20181130a_iii} becomes
\begin{equation*}
\ensuremath{\mathsf E}[\ensuremath{{\mathsf{F}}}(\ensuremath x^{\ensuremath n})] - \ensuremath{{\mathsf{F}}}_* \leq \bigg( 1 - \ensuremath \mathsf{p}_{\min}\frac{2 (\mu + \sigma)}{\nu_{\max} + \mu + 2\sigma} \bigg)^\ensuremath n \mathrm{const}.
\end{equation*}
If $\sigma =0$ and the algorithm is deterministic and fully parallel (i.e., $\ensuremath{\mathsf P}(\sum_{i=1}^m \varepsilon_i =m)=1$),
then
\begin{equation*}
\ensuremath{{\mathsf{F}}}(\ensuremath x^{\ensuremath n}) - \ensuremath{{\mathsf{F}}}_* \leq \bigg(\frac{\nu_{\max}-\mu}{\nu_{\max} + \mu} \bigg)^\ensuremath n \mathrm{const},
\end{equation*}
which agrees with the classic result \cite[Theorem~2.1.15]{book4}
when the stepsizes are all equal to
$1/\nu_{\max}$.
\end{remark}

\begin{remark}
Suppose that the block sampling is uniform, that is,
$\ensuremath \mathsf{p}_i = \ensuremath \mathsf{p}$ for all $i \in [m]$ and let, for every $i \in [m]$, $\gamma_i = 1/\nu_i$.
Then $\delta = 1$ and
Theorem~\ref{p:20181130a}\ref{p:20181130a_ii}-\ref{p:20181130a_iii} reduce to
\begin{align*}
\ensuremath{\mathsf E} [\ensuremath{{\mathsf{F}}}(\ensuremath x^{\ensuremath n})] - \ensuremath{{\mathsf{F}}}_*
&=\frac{1}{1 + \ensuremath \mathsf{p} \ensuremath n}
\bigg(\frac{\mathrm{dist}_{\ensuremath{{\Gammas}}^{-1}}^2(\ensuremath \xx^0, \ensuremath {\mathsf{S}}_*)}{2}
+ \ensuremath{{\mathsf{F}}}(\ensuremath \xx^{0}) - \ensuremath{{\mathsf{F}}}_* \bigg)\qquad\text{and}\\[1ex]
\ensuremath{\mathsf E}[\ensuremath{{\mathsf{F}}}(\ensuremath x^{\ensuremath n})] - \ensuremath{{\mathsf{F}}}_*
&\leq \bigg(1 - \ensuremath \mathsf{p}\frac{2(\mu_{\ensuremath \mathsf{\Gamma}^{-1}} + \sigma_{\ensuremath \mathsf{\Gamma}^{-1}})}{1 + \mu_{\ensuremath \mathsf{\Gamma}^{-1}} + 2 \sigma_{\ensuremath \mathsf{\Gamma}^{-1}}}\bigg)^\ensuremath n
\bigg(  (1+\sigma_{\ensuremath \mathsf{\Gamma}^{-1}}) \frac{\mathrm{dist}_{\ensuremath{{\Gammas}}^{-1}}^2(\ensuremath \xx^0, \ensuremath {\mathsf{S}}_*)}{2}
+ \ensuremath{{\mathsf{F}}}(\ensuremath \xx^0) - \ensuremath{{\mathsf{F}}}_* \bigg)
\end{align*}
respectively. These results were obtained in \cite[Theorem~3]{Tap18},
which are in turn generalizations of \cite[Theorem~1]{LuX15},
treating the serial case ($\ensuremath{\mathsf P}(\sum_{i=1}^m \varepsilon_i = 1)=1$).
\end{remark}

\subsection{Linear convergence under error bound conditions}
\label{sec:errorbounds}

In this section we analyze the convergence of Algorithm~\ref{algoRCD}
under error bound conditions. We improve and simplify
the results given in \cite{Nec16}.
In the rest of the section we assume \ref{eq:A1} and \ref{eq:A2}.
Moreover, we let $\ensuremath{{\XX}} \subset \ensuremath {\mathsf{H}}$,  $\ensuremath{{\mathsf{F}}} = \ensuremath{{\mathsf{f}}} + \ensuremath{{\gs}}$, $\ensuremath{{\mathsf{F}}}_* = \inf \ensuremath{{\mathsf{F}}}$,
and suppose $\ensuremath {\mathsf{S}}_*:=\argmin \ensuremath{{\mathsf{F}}} \neq \varnothing$.

We consider the following
condition,  which
was studied in \cite{Dru18} in connection with
the proximal gradient method
and is known as \emph{Luo-Tseng error bound} condition \cite{Luo93}.

\begin{enumerate}[{\rm EB}]
\item\label{eq:EB1} For some $c_{\ensuremath{{\XX}},\ensuremath{{\Gammas}}^{-1}}>0$, we have
\begin{equation}
(\forall\,\ensuremath \xx \in \ensuremath{{\XX}})\quad\mathrm{dist}_{\ensuremath{{\Gammas}}^{-1}}(\ensuremath \xx, \ensuremath {\mathsf{S}}_*)
\leq c_{\ensuremath{{\XX}},\ensuremath{{\Gammas}}^{-1}}
\norm{\ensuremath \xx - \ensuremath{\text{\sf prox}}^{\ensuremath{{\Gammas}}^{-1}}_{\ensuremath{{\gs}}} (\ensuremath \xx - \nabla^{\ensuremath{{\Gammas}}^{-1}}\ensuremath{{\mathsf{f}}}(\ensuremath \xx) )}_{\ensuremath{{\Gammas}}^{-1}}.
\end{equation}
\end{enumerate}

 \begin{remark}\
 \label{rmk:20190603a}
 \begin{enumerate}[{\rm (i)}]
\item Another popular error bound condition is that of the metric subregularity of the subdifferential.
More precisely, $\partial^{\ensuremath{{\Gammas}}^{-1}} \ensuremath{{\mathsf{F}}}$ is $2$-\emph{metrically subregular}
on $\ensuremath{{\XX}}$ with respect to the metric $\norm{\cdot}_{\ensuremath{{\Gammas}}^{-1}}$\cite{Dru18,Gar17}
if for some
 $\zeta_{\ensuremath{{\XX}},\ensuremath{{\Gammas}}^{-1}}>0$ the following holds
\begin{equation}
\label{eq:20181114h}
(\forall\, \ensuremath \xx \in \ensuremath{{\XX}})\quad
\mathrm{dist}_{\ensuremath{{\Gammas}}^{-1}}(\ensuremath \xx, \ensuremath {\mathsf{S}}_*) \leq \frac{1}{\zeta_{\ensuremath{{\XX}},\ensuremath{{\Gammas}}^{-1}}}
\mathrm{dist}_{\ensuremath{{\Gammas}}^{-1}}(\ensuremath{{0}}, \partial^{\ensuremath{{\Gammas}}^{-1}} \ensuremath{{\mathsf{F}}}(\ensuremath \xx)).
\end{equation}
 \item\label{rmk:20190129e_00} \ref{eq:EB1} and \eqref{eq:20181114h} are equivalent if $\ensuremath{{\gs}}=0$,
 since in that case $\ensuremath{\text{\sf prox}}_{\ensuremath{{\gs}}}^{\ensuremath{{\Gammas}}^{-1}} = \ensuremath{\mathsf{Id}}$
 and $c_{\ensuremath{{\XX}},\ensuremath{{\Gammas}}^{-1}} = \zeta_{\ensuremath{{\XX}}, \ensuremath{{\Gammas}}^{-1}}^{-1}$.
 \item\label{rmk:20190129e_0} Since $\partial^{\ensuremath{{\Gammas}}^{-1}} \ensuremath{{\mathsf{F}}} (\ensuremath x) = \ensuremath{{\Gammas}} \partial \ensuremath{{\mathsf{F}}}(\ensuremath x)$
 and $\norm{\cdot} \geq \gamma_{\min}^{1/2}\norm{\cdot}_{\ensuremath{{\Gammas}}^{-1}}$,
it follows that if for every $\ensuremath x \in \ensuremath{{\XX}}$, $
\mathrm{dist}(\ensuremath \xx, \ensuremath {\mathsf{S}}_*) \leq \zeta^{-1}_{\ensuremath{{\XX}},\ensuremath{\mathsf{Id}}}
\mathrm{dist}(\ensuremath{{0}}, \partial \ensuremath{{\mathsf{F}}}(\ensuremath \xx))$,
then \eqref{eq:20181114h} holds with constant $\zeta_{\ensuremath{{\XX}}, \ensuremath{{\Gammas}}^{-1}} = \gamma_{\min} \zeta_{\ensuremath{{\XX}},\ensuremath{\mathsf{Id}}}$.
\item\label{rmk:20190603a_iv} \cite[Theorem~3.5]{Dru18} yields that
for any Hilbert norm $\norm{\cdot}_\ensuremath{{\WW}}$,
 $\norm{\ensuremath \xx - \ensuremath{\text{\sf prox}}^{\ensuremath{{\WW}}}_{\ensuremath{{\gs}}} (\ensuremath \xx - \nabla^{\ensuremath{{\WW}}}\ensuremath{{\mathsf{f}}}(\ensuremath \xx) )}_{\ensuremath{{\WW}}}
\leq \mathrm{dist}_{\ensuremath{{\WW}}}(\ensuremath{{0}}, \partial^{\ensuremath{{\WW}}}\ensuremath{{\mathsf{F}}}(\ensuremath \xx))$. So,
 if \ref{eq:EB1} holds on $\ensuremath{{\XX}}$, then \eqref{eq:20181114h} holds on $\ensuremath{{\XX}}$
 with $\zeta_{\ensuremath{{\XX}},\ensuremath{{\Gammas}}^{-1}} = c_{\ensuremath{{\XX}},\ensuremath{{\Gammas}}^{-1}}^{-1}$.
 In \cite[Theorem~3.4-3.5]{Dru18} also the reverse implication was shown when
 $\ensuremath{{\mathsf{f}}}$ is Lipschitz smooth and
 $\ensuremath{{\XX}}$ is a sublevel set of $\ensuremath{{\mathsf{F}}}$.
 \end{enumerate}
 \end{remark}

\begin{remark}
\label{rmk:20190129e}
In \cite[Corollary~3.6]{Dru18} condition \ref{eq:EB1} was shown to be equivalent to the
following \emph{quadratic growth condition} (also called $2$-conditioning in \cite{Gar17})
\begin{equation}
\label{eq:20181218b}
(\exists \alpha_{\ensuremath{{\XX}},\ensuremath{{\Gammas}}^{-1}} >0)(\forall\, \ensuremath \xx \in \ensuremath{{\XX}})\quad \ensuremath{{\mathsf{F}}}(\ensuremath \xx) - \inf \ensuremath{{\mathsf{F}}}
\geq \frac{\alpha_{\ensuremath{{\XX}},\ensuremath{{\Gammas}}^{-1}}}{2} \mathrm{dist}_{\ensuremath{{\Gammas}}^{-1}}^2(\ensuremath \xx, \ensuremath {\mathsf{S}}_*),
\end{equation}
on every sublevel set $\ensuremath{{\XX}} = \{\ensuremath \xx \in \ensuremath {\mathsf{H}} \,\vert\, \ensuremath{{\mathsf{F}}}(\ensuremath \xx) - \ensuremath{{\mathsf{F}}}_* \leq r\}$.
Moreover, the relationships between the constants are $c_{\ensuremath{{\XX}},\ensuremath{{\Gammas}}^{-1}} = (1 + 2/\alpha_{\ensuremath{{\XX}},\ensuremath{{\Gammas}}^{-1}})(1+ L_{\norm{\cdot}_{\ensuremath{{\Gammas}}^{-1}}})$
and $\alpha_{\ensuremath{{\XX}},\ensuremath{{\Gammas}}^{-1}} <1/c_{\ensuremath{{\XX}},\ensuremath{{\Gammas}}^{-1}}$.
Finally,  if the quadratic growth condition \eqref{eq:20181218b}
holds, then \eqref{eq:20181114h} holds on $\ensuremath{{\XX}}$,
with $\zeta_{\ensuremath{{\XX}},\ensuremath{{\Gammas}}^{-1}} = \alpha_{\ensuremath{{\XX}},\ensuremath{{\Gammas}}^{-1}}/2$.
\end{remark}

We now analyze the convergence of Algorithm~\ref{algoRCD}
under condition \ref{eq:EB1}.

\begin{theorem}
\label{thm:20181214a}
Under the assumptions of Theorem~\ref{thm:20171207a},
suppose that $\ensuremath {\mathsf{S}}_* \neq \varnothing$ and that
\ref{eq:EB1} holds on a set $\ensuremath{{\XX}}$ such that $\ensuremath{{\XX}} \supset \big\{ \ensuremath x^\ensuremath n \,\vert\, \ensuremath n \in \ensuremath \mathbb{N}\big\}$ $\ensuremath{\mathsf P}$-a.s.~with $c_{\ensuremath{{\XX}},\ensuremath{{\Gammas}}^{-1}}>0$.
Then,
\begin{equation}
\label{eq:20190131e}
(\forall\, \ensuremath n \in \ensuremath \mathbb{N})\quad
\ensuremath{\mathsf E}[\ensuremath{{\mathsf{F}}}(\ensuremath x^{\ensuremath n})] - \ensuremath{{\mathsf{F}}}_* \leq \bigg( 1-  \ensuremath \mathsf{p}_{\min}\min\bigg\{1,\frac{2 - \delta}{2 c_{\ensuremath{{\XX}},\ensuremath{{\Gammas}}^{-1}}} \bigg\} \bigg)^\ensuremath n
\big(\ensuremath{\mathsf E}[\ensuremath{{\mathsf{F}}}(\ensuremath x^{0})] - \ensuremath{{\mathsf{F}}}_* \big).
\end{equation}
Moreover,
there exists a random variable $\ensuremath x_*$ which takes values in $\ensuremath {\mathsf{S}}_*$ such that
$\ensuremath x^\ensuremath n \to \ensuremath x_*$ $\ensuremath{\mathsf P}$-a.s.~and
$\ensuremath{\mathsf E}[\norm{\ensuremath x^{n} - \ensuremath x_*}_{\ensuremath{{\WW}}}] = O\big(
( 1-  \ensuremath \mathsf{p}_{\min}\min\{1,(2 - \delta)/(2 c_{\ensuremath{{\XX}},\ensuremath{{\Gammas}}^{-1}})\})^{n/2} \big)$.
\end{theorem}
\begin{proof}
Let $\ensuremath n \in \ensuremath \mathbb{N}$ and $\ensuremath \xx \in \ensuremath {\mathsf{S}}_*$. Then,
 \eqref{eq:20190610a}
with $\mu_{\ensuremath \mathsf{\Gamma}^{-1}}=\sigma_{\ensuremath \mathsf{\Gamma}^{-1}}=0$, yields
\begin{align}
\nonumber
\frac{1}{\ensuremath \mathsf{p}_{\min}}
\ensuremath{\mathsf E}[  &\ensuremath{{\mathsf{F}}}(\ensuremath x^{\ensuremath n+1}) - \ensuremath{{\mathsf{F}}}(\ensuremath x^{\ensuremath n}) \,\vert\, \ensuremath{\mathfrak{E}}_{\ensuremath n-1} ]\\
\label{eq:20181213a}
&\leq  \frac{\delta - 2}{2} \norm{ \ensuremath x^\ensuremath n - \bar{\ensuremath x}^{\ensuremath n+1}}_{\ensuremath{{\Gammas}}^{-1}}^2
+ \norm{\ensuremath x^{\ensuremath n} - \bar{\ensuremath x}^{\ensuremath n+1}}_{\ensuremath{{\Gammas}}^{-1}} \norm{\ensuremath x^\ensuremath n - \ensuremath \xx}_{\ensuremath{{\Gammas}}^{-1}}-  \big( \ensuremath{{\mathsf{F}}}(\ensuremath x^\ensuremath n) - \ensuremath{{\mathsf{F}}}_* \big).
\end{align}
Since \eqref{eq:20181213a} holds for all $\ensuremath \xx \in \ensuremath {\mathsf{S}}_*$,
using  \ref{eq:EB1}, \eqref{eq:20180925g}, and Proposition~\ref{p:20190313c},
we have
\begin{align}
\nonumber
\frac{2}{\ensuremath \mathsf{p}_{\min}}&
\ensuremath{\mathsf E}[  \ensuremath{{\mathsf{F}}}(\ensuremath x^{\ensuremath n+1}) - \ensuremath{{\mathsf{F}}}(\ensuremath x^\ensuremath n) \,\vert\, \ensuremath{\mathfrak{E}}_{\ensuremath n-1} ] \\[1ex]
\nonumber&\leq  (\delta - 2) \norm{ \ensuremath x^\ensuremath n - \bar{\ensuremath x}^{\ensuremath n+1}}_{\ensuremath{{\Gammas}}^{-1}}^2
+ 2 \norm{\ensuremath x^{\ensuremath n} - \bar{\ensuremath x}^{\ensuremath n+1}}_{\ensuremath{{\Gammas}}^{-1}} \mathrm{dist}_{\ensuremath{{\Gammas}}^{-1}}(\ensuremath x^\ensuremath n, \ensuremath {\mathsf{S}}_*)
- 2 \big( \ensuremath{{\mathsf{F}}}(\ensuremath x^\ensuremath n) - \ensuremath{{\mathsf{F}}}_* \big)\\[1ex]
\label{eq:20181214g}&\leq \frac{2(2 c_{\ensuremath{{\XX}},\ensuremath{{\Gammas}}^{-1}} + \delta - 2)_+}{\ensuremath \mathsf{p}_{\min}(2 - \delta)}
\ensuremath{\mathsf E}[  \ensuremath{{\mathsf{F}}}(\ensuremath x^{\ensuremath n}) - \ensuremath{{\mathsf{F}}}(\ensuremath x^{\ensuremath n+1}) \,\vert\, \ensuremath{\mathfrak{E}}_{\ensuremath n-1} ]
- 2 ( \ensuremath{{\mathsf{F}}}(\ensuremath x^\ensuremath n) - \ensuremath{{\mathsf{F}}}_*)\quad \ensuremath{\mathsf P} \text{ a.s.},
\end{align}
which can be equivalently written as
\begin{multline*}
\ensuremath{\mathsf E}[  \ensuremath{{\mathsf{F}}}(\ensuremath x^{\ensuremath n+1}) - \ensuremath{{\mathsf{F}}}_*  \,\vert\, \ensuremath{\mathfrak{E}}_{\ensuremath n-1} ]
- \big(\ensuremath{{\mathsf{F}}}(\ensuremath x^\ensuremath n) - \ensuremath{{\mathsf{F}}}_* \big) \\[1ex]
\leq \bigg(\frac{2 c_{\ensuremath{{\XX}},\ensuremath{{\Gammas}}^{-1}}}{2 - \delta} -1\bigg)_+\Big(
\ensuremath{{\mathsf{F}}}(\ensuremath x^{\ensuremath n}) - \ensuremath{{\mathsf{F}}}_*  - \ensuremath{\mathsf E}[  \ensuremath{{\mathsf{F}}}(\ensuremath x^{\ensuremath n+1}) - \ensuremath{{\mathsf{F}}}_* \,\vert\, \ensuremath{\mathfrak{E}}_{\ensuremath n-1} ] \Big)
- \ensuremath \mathsf{p}_{\min} ( \ensuremath{{\mathsf{F}}}(\ensuremath x^\ensuremath n) - \ensuremath{{\mathsf{F}}}_* ).
\end{multline*}
Therefore,
\begin{equation*}
\max\bigg\{1,\frac{2 c_{\ensuremath{{\XX}},\ensuremath{{\Gammas}}^{-1}}}{2 - \delta}\bigg\}
\ensuremath{\mathsf E}[  \ensuremath{{\mathsf{F}}}(\ensuremath x^{\ensuremath n+1}) - \ensuremath{{\mathsf{F}}}_* \,\vert\, \ensuremath{\mathfrak{E}}_{\ensuremath n-1} ]
\leq \bigg(\max\bigg\{1,\frac{2 c_{\ensuremath{{\XX}},\ensuremath{{\Gammas}}^{-1}}}{2 - \delta}\bigg\} -\ensuremath \mathsf{p}_{\min}\bigg)
\big( \ensuremath{{\mathsf{F}}}(\ensuremath x^{\ensuremath n}) - \ensuremath{{\mathsf{F}}}_* \big),
\end{equation*}
which gives \eqref{eq:20190131e}.
Now we set $\rho = 1-  \ensuremath \mathsf{p}_{\min}\min\big\{1,(2 - \delta)/(2 c_{\ensuremath{{\XX}},\ensuremath{{\Gammas}}^{-1}}) \big\}$ and $\theta = \ensuremath \mathsf{p}_{\min}(2 - \delta)/2$.
Then, Jensen inequality, \eqref{eq:20190128a}, and \eqref{eq:20190131e} yield
\begin{equation}
\label{eq:20190131g}
\ensuremath{\mathsf E}[\norm{\ensuremath x^n - \ensuremath x^{\ensuremath n+1}}_{\ensuremath{{\WW}}}]
\leq \theta^{-1/2}
\sqrt{\ensuremath{\mathsf E}[\ensuremath{{\mathsf{F}}}(\ensuremath x^\ensuremath n)] - \ensuremath{{\mathsf{F}}}_*}
 \leq \theta^{-1/2}
\rho^{\ensuremath n/2}\sqrt{\ensuremath{\mathsf E}[\ensuremath{{\mathsf{F}}}(\ensuremath x^0)] - \ensuremath{{\mathsf{F}}}_*}.
\end{equation}
Therefore, since $\rho^{1/2}<1$, we have
$\ensuremath{\mathsf E}[\sum_{\ensuremath n \in \ensuremath \mathbb{N}} \norm{\ensuremath x^\ensuremath n - \ensuremath x^{\ensuremath n+1}}_{\ensuremath{{\WW}}}]
= \sum_{\ensuremath n \in \ensuremath \mathbb{N}} \ensuremath{\mathsf E}[\norm{\ensuremath x^\ensuremath n - \ensuremath x^{\ensuremath n+1}}_{\ensuremath{{\WW}}}] <+\infty$.
Hence $\sum_{\ensuremath n \in \ensuremath \mathbb{N}} \norm{\ensuremath x^\ensuremath n - \ensuremath x^{\ensuremath n+1}}_{\ensuremath{{\WW}}}<+\infty$ $\ensuremath{\mathsf P}$-a.s.,
which means that $(\ensuremath x^\ensuremath n)_{n \in \ensuremath \mathbb{N}}$ is a Cauchy sequence $\ensuremath{\mathsf P}$-a.s.
Now, Theorem~\ref{thm:20171207a}\ref{thm:20171207a_i} yields that there exists a random variable $\ensuremath x_*$
with values in $\ensuremath {\mathsf{S}}_*$ such that $\ensuremath x^n \rightharpoonup \ensuremath x_*$
$\ensuremath{\mathsf P}$-a.s. Therefore, $\ensuremath x^n \to \ensuremath x_*$ $\ensuremath{\mathsf P}$-a.s.
Finally, let $\ensuremath n \in \ensuremath \mathbb{N}$. Then, for every $p \in \ensuremath \mathbb{N}$,
\begin{equation*}
\norm{\ensuremath x^{\ensuremath n} - \ensuremath x^{\ensuremath n+p}}_{\ensuremath{{\WW}}} \leq
\sum_{i=0}^{p-1} \norm{\ensuremath x^{n + i} - \ensuremath x^{\ensuremath n + i + 1}}_{\ensuremath{{\WW}}}
\leq \sum_{i=0}^{+\infty} \norm{\ensuremath x^{\ensuremath n + i} - \ensuremath x^{\ensuremath n + i + 1}}_\ensuremath{{\WW}}.
\end{equation*}
Hence, letting $p\to +\infty$, we have
$\ensuremath{\mathsf E}[\norm{\ensuremath x^{\ensuremath n} - \ensuremath x_*}_\ensuremath{{\WW}}] \leq \sum_{i=0}^{+\infty} \ensuremath{\mathsf E}[\norm{\ensuremath x^{\ensuremath n + i} - \ensuremath x^{\ensuremath n + i + 1}}_\ensuremath{{\WW}}]$. Therefore,  it follows from \eqref{eq:20190131g} that
\begin{align*}
\ensuremath{\mathsf E}[\norm{\ensuremath x^{n} - \ensuremath x_*}_\ensuremath{{\WW}}]
& \leq \theta^{-1/2}
\sqrt{\ensuremath{\mathsf E}[\ensuremath{{\mathsf{F}}}(\ensuremath x^0)] - \ensuremath{{\mathsf{F}}}_*}
\sum_{i=0}^{+\infty}  \rho^{(n+i)/2}
 = \theta^{-1/2}
\sqrt{\ensuremath{\mathsf E}[\ensuremath{{\mathsf{F}}}(\ensuremath x^0)] - \ensuremath{{\mathsf{F}}}_*}
\frac{\rho^{n/2}}{1-\rho^{1/2}}.
\qedhere
\end{align*}
\end{proof}

\begin{remark}\
\label{rmk:20181219e}
\begin{enumerate}[(i)]
\item The rate given in Theorem~\ref{thm:20181214a} matches the one given in
\cite[Theorem~3.2]{Dru18}
for the deterministic case ($\ensuremath \mathsf{p}_{\min}=1$).
\item In \cite[Definition~5.2]{Nec16}, in relation to Algorithm~\ref{algoRCD} but with uniform block
sampling and assuming \ref{eq:S0} and $R_{\ensuremath{{\Gammas}}^{-1}}(\ensuremath \xx^0):=\sup_{\ensuremath \xx \in \{\ensuremath{{\mathsf{F}}} \leq \ensuremath{{\mathsf{F}}}(\ensuremath \xx^0)\} } \mathrm{dist}_{\ensuremath{{\Gammas}}^{-1}}(\ensuremath \xx, \ensuremath {\mathsf{S}}_*)<+\infty$, the following error bound condition
is considered
\begin{equation}
\label{eq:20181114b}
\mathrm{dist}_{\ensuremath{{\Gammas}}^{-1}}(\ensuremath \xx, \ensuremath {\mathsf{S}}_*)
\leq \big(\kappa_{1,\ensuremath{{\XX}},\ensuremath{{\Gammas}}^{-1}} + \kappa_{2,\ensuremath{{\XX}},\ensuremath{{\Gammas}}^{-1}}
\mathrm{dist}_{\ensuremath{{\Gammas}}^{-1}}^2(\ensuremath \xx, \ensuremath {\mathsf{S}}_*)\big)
\norm{\ensuremath \xx - \ensuremath{\text{\sf prox}}^{\ensuremath{{\Gammas}}^{-1}}_{\ensuremath{{\gs}}} (\ensuremath \xx - \nabla^{\ensuremath{{\Gammas}}^{-1}}\ensuremath{{\mathsf{f}}}(\ensuremath \xx) )}_{\ensuremath{{\Gammas}}^{-1}},
\end{equation}
for some constants $\kappa_{1,\ensuremath{{\XX}},\ensuremath{{\Gammas}}^{-1}}>0$ and $\kappa_{2,\ensuremath{{\XX}},\ensuremath{{\Gammas}}^{-1}} \geq 0$.
The authors show several examples in which such condition is satisfied
with $\ensuremath{{\XX}}= \dom \ensuremath{{\gs}}$ and possibly $\kappa_{2,\ensuremath{{\XX}},\ensuremath{{\Gammas}}^{-1}}>0$.
The above error bound looks
more general then \ref{eq:EB1}. However, for the purpose of analyzing Algorithm~\ref{algoRCD}
and under the assumptions considered in \cite{Nec16} this is not the case. Indeed,
in \cite[equation~(3.11)]{Nec16} it was shown
that the algorithm is descending almost surely\footnote{Alternatively,
note that \ref{eq:S0} implies \ref{eq:S2} which in turn,
in view of Proposition~\ref{p:20181219c}, ensures the descending property.},
so  $\{\ensuremath x^\ensuremath n \,\vert\, \ensuremath n \in \ensuremath \mathbb{N}\} \subset \{\ensuremath{{\mathsf{F}}} \leq \ensuremath{{\mathsf{F}}}(\ensuremath \xx^0)\}$ $\ensuremath{\mathsf P}$-a.s.
Therefore,
since $\sup_{\ensuremath \xx \in \{\ensuremath{{\mathsf{F}}} \leq \ensuremath{{\mathsf{F}}}(\ensuremath \xx^0)\} } \mathrm{dist}_{\ensuremath{{\Gammas}}^{-1}}(\ensuremath \xx, \ensuremath {\mathsf{S}}_*)
= R_{\ensuremath{{\Gammas}}^{-1}}(\ensuremath \xx^0)<+\infty$,
if  \eqref{eq:20181114b} holds
on a set $\ensuremath{{\XX}}$ containing $\ensuremath{\mathsf P}$-a.s.~the set $\{\ensuremath x^\ensuremath n \,\vert\, \ensuremath n \in \ensuremath \mathbb{N}\}$,
then \ref{eq:EB1} holds on $\ensuremath{{\XX}}^\prime := \ensuremath{{\XX}} \cap \{\ensuremath{{\mathsf{F}}} \leq \ensuremath{{\mathsf{F}}}(\ensuremath \xx^0)\} \supset
\{\ensuremath x^\ensuremath n \,\vert\, \ensuremath n \in \ensuremath \mathbb{N}\}$ $\ensuremath{\mathsf P}$-a.s.~with
$c_{\ensuremath{{\XX}}^\prime,\ensuremath{{\Gammas}}^{-1}} := \kappa_{1,\ensuremath{{\XX}},\ensuremath{{\Gammas}}^{-1}} + \kappa_{2,\ensuremath{{\XX}},\ensuremath{{\Gammas}}^{-1}} R_{\ensuremath{{\Gammas}}^{-1}}(\ensuremath \xx^0)^2$.
Thus, Theorem~\ref{thm:20181214a}
applies accordingly. Moreover,
\cite[Theorem~5.5]{Nec16} gives the linear rate
\begin{multline*}
\ensuremath{\mathsf E}[\ensuremath{{\mathsf{F}}}(\ensuremath x^{\ensuremath n})] - \ensuremath{{\mathsf{F}}}_* \leq \bigg(1 - \frac{1}{1+\bar{c}} \bigg)^\ensuremath n (\ensuremath{\mathsf E}[\ensuremath{{\mathsf{F}}}(\ensuremath \xx^{0})] - \ensuremath{{\mathsf{F}}}_*),
\quad\text{where}\\
\bar{c} = \frac{1}{\ensuremath \mathsf{p}} \bigg( 2 + \frac{2 c_{\ensuremath{{\XX}}^\prime, \ensuremath{{\Gammas}}^{-1}}}{\sqrt{\ensuremath \mathsf{p}}} + (1 - \ensuremath \mathsf{p})
\frac{c_{\ensuremath{{\XX}}^\prime, \ensuremath{{\Gammas}}^{-1}}^2}{\ensuremath \mathsf{p}} + 2 c_{\ensuremath{{\XX}}^\prime, \ensuremath{{\Gammas}}^{-1}} + 1 - \ensuremath \mathsf{p} \bigg).
\end{multline*}
Then, we have $1+ \bar{c} \geq (3 + 4 c_{\ensuremath{{\XX}}^\prime,\ensuremath{{\Gammas}}^{-1}})/\ensuremath \mathsf{p}$ and hence
\begin{equation*}
\frac{1}{1 + \bar{c}} \leq \frac{\ensuremath \mathsf{p}}{3 + 4 c_{\ensuremath{{\XX}}^\prime, \ensuremath{{\Gammas}}^{-1}}} < \frac{\ensuremath \mathsf{p}}{\max\{1,2c_{\ensuremath{{\XX}}^\prime, \ensuremath{{\Gammas}}^{-1}}\}}.
\end{equation*}
This shows that Theorem~\ref{thm:20181214a}  improves the rate in \cite[Theorem~5.5]{Nec16}. Moreover, the analysis given here,
relying on Proposition~\ref{p:20180927a}, relaxes the assumptions and is significantly simpler.
\item\label{rmk:20181219e_iii}
It follows from \cite[Theorem~6.8]{Nec16} that if $\ensuremath{{\mathsf{f}}}$ is a quadratic function and
$\ensuremath{{\gs}}$ is an indicator function of a polyhedral set, then
\eqref{eq:20181114b} is satisfied on $\dom \ensuremath{{\gs}}$. Therefore, if $\dom \ensuremath{{\gs}}$ is bounded,
then \ref{eq:EB1} holds on $\ensuremath{{\XX}} =\dom \ensuremath{{\gs}}$ with
$c_{\ensuremath{{\XX}},\ensuremath{{\Gammas}}^{-1}} := \kappa_{1,\ensuremath{{\XX}},\ensuremath{{\Gammas}}^{-1}} + \kappa_{2,\ensuremath{{\XX}},\ensuremath{{\Gammas}}^{-1}} \diam^2_{\ensuremath{{\Gammas}}^{-1}}(\dom \ensuremath{{\gs}})$ and Theorem~\ref{thm:20181214a} can be applied, since
$\{x^\ensuremath n \,\vert\, \ensuremath n \in \ensuremath \mathbb{N}\} \subset \dom \ensuremath{{\gs}}$.
\item Several works address the convergence of random coordinate descent methods
under error bound conditions. We mention \cite{Kar18} which considers a serial sampling and stepsizes
 set according to the global Lipschitz constant of $\nabla \ensuremath{{\mathsf{f}}}$ and \cite{Fer20},
which analyzes restarting procedures for accelerated and parallel coordinate descent methods
using assumptions \ref{eq:S1} and \eqref{eq:20181218b}.
\end{enumerate}
\end{remark}

\begin{remark}
\label{rmk:20200331a}
Often error bound conditions or quadratic growth conditions are satisfied when $\ensuremath{{\XX}}$ is a sublevel set (see Remark~\ref{rmk:20190129e}). So, in such scenarios,
apart when $\dom \ensuremath{{\gs}}$ is a sublevel set of $\ensuremath{{\mathsf{F}}}$,
in order to fulfill the assumption
$\ensuremath{{\XX}} \supset \{\ensuremath x^\ensuremath n \,\vert\, \ensuremath n \in \ensuremath \mathbb{N}\}$ $\ensuremath{\mathsf P}$-a.s.~in Theorem~\ref{thm:20181214a},
it is desirable for Algorithm~\ref{algoRCD} to be (a.s.) descending.
This occurs if condition \ref{eq:S2} holds (Proposition~\ref{p:20181219c}),
whereas, in general, \ref{eq:S1} does not guarantee any such
descending property.
However, especially when $\eta \ll m$, condition \ref{eq:S2} may be much more restrictive
than \ref{eq:S1}, thus leading to a significant reduction of the stepsizes, which ultimately
slows down the convergence.
The next result shows that
Algorithm~\ref{algoRCD} can be slightly modified so to ensure the descending property
while keeping the validity of Theorem~\ref{thm:20181214a}.
\end{remark}

\begin{theorem}
\label{thm20200218a}
 Let \ref{eq:A1}--\ref{eq:A3} be satisfied.
Let $(\nu_{i})_{1 \leq i \leq m} \in \ensuremath \mathbb{R}^m_{++}$ and suppose that
 \ref{eq:S1} holds. Suppose in addition that
 $\ensuremath {\mathsf{S}}_* \neq \varnothing$, and that
\ref{eq:EB1} holds on the set $\ensuremath{{\XX}} = \{\ensuremath{{\mathsf{F}}} \leq\ensuremath{{\mathsf{F}}}(\ensuremath \mathsf{x}^0)\}$
with $c_{\ensuremath{{\XX}},\ensuremath{{\Gammas}}^{-1}}>0$.
Let $(\ensuremath x^\ensuremath n)_{\ensuremath n \in \ensuremath \mathbb{N}}$ be generated by the following variation of Algorithm~\ref{algoRCD}
\begin{equation}
\label{eq:algoPRCD3}
\begin{array}{l}
\text{for}\;n=0,1,\ldots\\
\left\lfloor
\begin{array}{l}
\text{for}\;i=1,\dots, m\\[0.7ex]
\left\lfloor
\begin{array}{l}
\tilde{x}^{\ensuremath n+1}_i =
x^\ensuremath n_i + \varepsilon^\ensuremath n_i \big(\ensuremath{\text{\sf prox}}_{\gamma_i \ensuremath{{\mathsf{h}}}_{i}}
\big(x^\ensuremath n_i - \gamma_{i} \nabla_i \ensuremath{{\mathsf{f}}} (\ensuremath x^{\ensuremath n})\big) - x^\ensuremath n_i\big)
\end{array}
\right.\\[1ex]
\text{if } \ensuremath{{\mathsf{F}}}(\tilde{x}^{\ensuremath n+1}) \leq \ensuremath{{\mathsf{F}}}(x^\ensuremath n)\\[0.8ex]
\quad x^{\ensuremath n+1}=\tilde{x}^{\ensuremath n+1}\\
\text{else}\\
\quad x^{\ensuremath n+1} = x^{\ensuremath n}.
\end{array}
\right.
\end{array}
\end{equation}
Then the conclusions of Theorem~\ref{thm:20181214a} still hold.
\end{theorem}
\begin{proof}
It follows from the definition of $\ensuremath x^{n+1}$ that
$\ensuremath{{\mathsf{F}}}(\ensuremath x^{n+1}) \leq \ensuremath{{\mathsf{F}}}(\ensuremath x^n)$. Therefore
$\ensuremath{{\XX}} \supset \{\ensuremath x^\ensuremath n \,\vert\, \ensuremath n \in \ensuremath \mathbb{N}\}$.
Recalling \eqref{eq:20171206a} and \eqref{eq:20181114c},
algorithm \eqref{eq:algoPRCD3} can be alternatively written as
\begin{equation}
\label{eq:algoPRCD4}
\begin{array}{l}
\text{for}\;n=0,1,\ldots\\
\left\lfloor
\begin{array}{l}
\text{for}\;i=1,\dots, m\\[0.7ex]
\left\lfloor
\begin{array}{l}
\bar{\ensuremath x}^{\ensuremath n+1}_i = \ensuremath{\text{\sf prox}}_{\gamma_{i} \ensuremath{{\mathsf{h}}}_i} ( x_i^{\ensuremath n}
- \gamma_i \nabla_i \ensuremath{{\mathsf{f}}}(\ensuremath x^\ensuremath n) ),\\
\end{array}
\right.\\[1ex]
\tilde{\ensuremath x}^{\ensuremath n+1} = \ensuremath x^{\ensuremath n} + \ensuremath \varepsilon^{\ensuremath n} \odot (\bar{\ensuremath x}^{\ensuremath n+1} - \ensuremath x^{\ensuremath n})\\
x^{\ensuremath n+1} = \tilde{\ensuremath x}^{\ensuremath n+1} \mathbf{1}_{\{\ensuremath{{\mathsf{F}}}(\tilde{x}^{\ensuremath n+1}) \leq \ensuremath{{\mathsf{F}}}(x^\ensuremath n)\}}
+ \ensuremath x^{\ensuremath n}\mathbf{1}_{\{\ensuremath{{\mathsf{F}}}(\tilde{x}^{\ensuremath n+1}) > \ensuremath{{\mathsf{F}}}(x^\ensuremath n)\}}
\end{array}
\right.
\end{array}
\end{equation}
and we have $\ensuremath{{\mathsf{F}}}(x^{\ensuremath n+1}) \leq \ensuremath{{\mathsf{F}}}(\tilde{x}^{\ensuremath n+1})$.
Then we can essentially repeat the argument in the proof of Theorem~\ref{thm:20181214a}.
First we note that \eqref{eq:20190610a} and hence
\eqref{eq:20181213a} holds with $x^{\ensuremath n+1}$ replaced by $\tilde{x}^{\ensuremath n+1}$.
This follows from the definition of $\tilde{\ensuremath x}^{\ensuremath n+1}$.
Moreover, also \eqref{eq:20181124a} holds with $x^{\ensuremath n+1}$ replaced by $\tilde{x}^{\ensuremath n+1}$
and hence we derive (with $\ensuremath \xx = \ensuremath x^{\ensuremath n}$ and $\sigma_{\ensuremath \mathsf{\Gamma}^{-1}}=0$) that
\begin{equation}
\label{20200218b}
(2 -\delta) \norm{\bar{\ensuremath x}^{\ensuremath n+1} - \ensuremath x^{\ensuremath n}}_{\ensuremath{{\Gammas}}^{-1}}^2 \leq
\frac{2}{\ensuremath \mathsf{p}_{\min}}
\ensuremath{\mathsf E}[  \ensuremath{{\mathsf{F}}}(\ensuremath x^\ensuremath n) - \ensuremath{{\mathsf{F}}}(\tilde{\ensuremath x}^{\ensuremath n+1}) \,\vert\, \ensuremath{\mathfrak{E}}_{\ensuremath n-1} ].
\end{equation}
Then, we have
\begin{align}
\nonumber
\frac{2}{\ensuremath \mathsf{p}_{\min}}&
\ensuremath{\mathsf E}[  \ensuremath{{\mathsf{F}}}(\tilde{\ensuremath x}^{\ensuremath n+1}) - \ensuremath{{\mathsf{F}}}(\ensuremath x^\ensuremath n) \,\vert\, \ensuremath{\mathfrak{E}}_{\ensuremath n-1} ] \\[1ex]
\nonumber&\leq  (\delta - 2) \norm{ \ensuremath x^\ensuremath n - \bar{\ensuremath x}^{\ensuremath n+1}}_{\ensuremath{{\Gammas}}^{-1}}^2
+ 2 \norm{\ensuremath x^{\ensuremath n} - \bar{\ensuremath x}^{\ensuremath n+1}}_{\ensuremath{{\Gammas}}^{-1}} \mathrm{dist}_{\ensuremath{{\Gammas}}^{-1}}(\ensuremath x^\ensuremath n, \ensuremath {\mathsf{S}}_*)
- 2 \big( \ensuremath{{\mathsf{F}}}(\ensuremath x^\ensuremath n) - \ensuremath{{\mathsf{F}}}_* \big)\\[1ex]
\label{eq:20200218a}&\leq \frac{2(2 c_{\ensuremath{{\XX}},\ensuremath{{\Gammas}}^{-1}} + \delta - 2)_+}{\ensuremath \mathsf{p}_{\min}(2 - \delta)}
\ensuremath{\mathsf E}[  \ensuremath{{\mathsf{F}}}(\ensuremath x^{\ensuremath n}) - \ensuremath{{\mathsf{F}}}(\tilde{\ensuremath x}^{\ensuremath n+1}) \,\vert\, \ensuremath{\mathfrak{E}}_{\ensuremath n-1} ]
- 2 ( \ensuremath{{\mathsf{F}}}(\ensuremath x^\ensuremath n) - \ensuremath{{\mathsf{F}}}_*)\quad \ensuremath{\mathsf P} \text{ a.s.}
\end{align}
and hence, since $\ensuremath{{\mathsf{F}}}(\ensuremath x^{\ensuremath n+1}) \leq \ensuremath{{\mathsf{F}}}(\tilde{x}^{\ensuremath n+1})$,  \eqref{eq:20181214g}
still holds (for the new definition of $\ensuremath x^{\ensuremath n+1}$). Thus,
 \eqref{eq:20190131e} follows. As for the second part of the statement,
we note that, since $\ensuremath{{\mathsf{F}}}(\ensuremath x^{\ensuremath n+1}) \leq \ensuremath{{\mathsf{F}}}(\tilde{x}^{\ensuremath n+1})$, by \eqref{20200218b}, we have $(2 -\delta) \norm{\bar{\ensuremath x}^{\ensuremath n+1} - \ensuremath x^{\ensuremath n}}_{\ensuremath{{\Gammas}}^{-1}}^2 \leq
(2/\ensuremath \mathsf{p}_{\min}) \ensuremath{\mathsf E}[  \ensuremath{{\mathsf{F}}}(\ensuremath x^\ensuremath n) - \ensuremath{{\mathsf{F}}}(\ensuremath x^{\ensuremath n+1}) \,\vert\, \ensuremath{\mathfrak{E}}_{\ensuremath n-1} ]$.
Moreover, it follows from the definitions of
$\ensuremath x^{\ensuremath n+1}$ and $\tilde{\ensuremath x}^{\ensuremath n+1}$ in algorithm~\eqref{eq:algoPRCD4} that
\begin{equation*}
\ensuremath x^{\ensuremath n+1} - \ensuremath x^\ensuremath n = (\tilde{\ensuremath x}^{\ensuremath n+1} - \ensuremath x^\ensuremath n) \mathbf{1}_{\{\ensuremath{{\mathsf{F}}}(\tilde{\ensuremath x}^{\ensuremath n+1}) \leq \ensuremath{{\mathsf{F}}}(\ensuremath x^{\ensuremath n})\}}
\end{equation*}
and hence, by Proposition~\ref{p:20190313c}, we have
\begin{equation*}
\ensuremath{\mathsf E}[\norm{\ensuremath x^{\ensuremath n+1} - \ensuremath x^\ensuremath n}_\ensuremath{{\WW}}^2 \,\vert\, \ensuremath{\mathfrak{E}}_{\ensuremath n-1} ] \leq
\ensuremath{\mathsf E}[\norm{\tilde{\ensuremath x}^{\ensuremath n+1} - \ensuremath x^\ensuremath n}_\ensuremath{{\WW}}^2 \,\vert\, \ensuremath{\mathfrak{E}}_{\ensuremath n-1} ]
= \norm{\bar{\ensuremath x}^{\ensuremath n+1} - \ensuremath x^\ensuremath n}^2_{\ensuremath{{\Gammas}}^{-1}}.
\end{equation*}
In the end \eqref{eq:20190128a} with $\sigma_{\ensuremath{{\Gammas}}^{-1}}=0$ still holds
and the proof can continue as in that of Theorem~\ref{thm:20181214a}.
\end{proof}

\vspace{-2ex}
\section{Applications}
\label{sec:app}

In this section we show some relevant optimization problems for which
the theoretical analysis of Algorithm~\ref{algoRCD}
can be particularly useful.

\label{sec:appl}

\subsection{The Lasso problem}
\label{subsec:lasso}
Many papers study the convergence of coordinate descent methods for the Lasso problem and recent works prove linear convergence (see e.g., \cite{Kar18,Nec16,Nec19}). In \cite{Kar18} a random serial update of blocks is considered while in \cite{Nec19} the general framework of feasible descend methods is analyzed which include
(nonrandom) cyclic coordinate methods. In the following we discuss our contribution comparing with \cite{Nec16}.
Let $\ensuremath{{\As}} \in \ensuremath \mathbb{R}^{p\times m}$ and $\ensuremath \bs \in \ensuremath \mathbb{R}^p$.
We consider the problem
\begin{equation}
\label{eq:20181015a}
\min_{\ensuremath \xx \in \ensuremath \mathbb{R}^m} \frac 1 2 \norm{\ensuremath{{\As}} \ensuremath \xx - \ensuremath \bs}_2^2 + \lambda \norm{\ensuremath \xx}_1 \quad(\lambda>0).
\end{equation}
We denote by
$\ensuremath \as^i$ and $\ensuremath \as_k$ the $i$-th column and $k$-th row of $\ensuremath{{\As}}$ respectively.
Since
\begin{equation*}
\frac 1 2 \norm{\ensuremath{{\As}} \ensuremath \xx - \ensuremath \bs}^2
= \frac1 2 \sum_{k=1}^p ( \scalarp{\ensuremath \as_k, \ensuremath \xx} - \ensuremath \mathsf{b}_k)^2
= \frac 1 2 \sum_{k=1}^p \Big(\sum_{i=1}^m \ensuremath \mathsf{a}_k^i \ensuremath \mathsf{x}_i - \ensuremath \mathsf{b}_k \Big)^2,
\end{equation*}
\ref{eq:B0} holds and $\eta = \max_{1 \leq k \leq p} \mathrm{card}(\mathrm{spt}(\ensuremath \as_k))$.
Moreover, since
 $\nabla_i \ensuremath{{\mathsf{f}}}(\ensuremath \xx) = \scalarp{\ensuremath \as^i, \ensuremath{{\As}} \ensuremath \xx - \ensuremath \bs}$
and recalling Remark~\ref{rmk:20190128e}\ref{rmk:20190128e_i},
 conditions \ref{eq:B1} and  \ref{eq:B2} are satisfied with  $L_i = \norm{\ensuremath \as^i}^2$
 and $L^{(k)} = \norm{\ensuremath \as_k}^2$ respectively.
Then, Algorithm~\ref{algoRCD} (assuming that each block is made of
one coordinate only)
writes as
\begin{equation}
\label{eq:20170915a}
\begin{aligned}
\ensuremath x^{\ensuremath n+1}
&= \ensuremath x^\ensuremath n
+  \sum_{i=1}^m \varepsilon_i^\ensuremath n \big[ \mathrm{soft}_{\gamma_{i} \lambda}
\big( x_{i}^\ensuremath n - \gamma_{i} {\ensuremath \as^i}^\top(\ensuremath{{\As}} \ensuremath x^\ensuremath n - \ensuremath \bs) \big) - x_{i}^\ensuremath n \big]\ensuremath \ee_{i},
\end{aligned}
\end{equation}
where the \emph{soft thresholding operator} $\mathrm{soft}_{\gamma_{i} \lambda}$
 is defined as $\mathrm{soft}_{\gamma_{i} \lambda} (t)
 = \mathrm{sign}(t)\max\{0, \abs{t}- \gamma_{i} \lambda\}$ and $(\ensuremath \ee_{i})_{1 \leq i \leq m}$
 is the canonical basis of $\ensuremath \mathbb{R}^m$.
Now, define $u^\ensuremath n = \ensuremath{{\As}} \ensuremath x^\ensuremath n - \ensuremath \bs$. Then multiplying the equation in \eqref{eq:20170915a}
by $\ensuremath{{\As}}$ and subtracting $\ensuremath \bs$ by both terms, the algorithm is equivalently written as
\begin{equation}
\label{eq:algoLasso}
\begin{array}{l}
\left\lfloor
\begin{array}{l}
\text{for}\;i \in \mathrm{spt}(\varepsilon^\ensuremath n)\\
\left\lfloor
\begin{array}{l}
\xi_i = \mathrm{soft}_{\gamma_{i} \lambda}
\big( x_{i}^\ensuremath n - \gamma_{i} {\ensuremath \as^i}^\top u^\ensuremath n \big) - x_{i}^\ensuremath n,
\end{array}
\right.\\[1ex]
\ensuremath x^{\ensuremath n+1} = \ensuremath x^{\ensuremath n}  + \sum_{i \in \mathrm{spt}(\varepsilon^\ensuremath n)}
\xi_i \ensuremath \ee_{i}\\[0.5ex]
u^{\ensuremath n+1}  = u^{\ensuremath n}  + \sum_{i \in \mathrm{spt}(\varepsilon^\ensuremath n)}
\xi_i \ensuremath \as^i,
\end{array}
\right.
\end{array}
\end{equation}
showing that each iteration costs $O(p\tau_{\max})$ multiplications,
where $\tau_{\max}$ is the maximum number of block updates per iteration.
We now address the determination of the smoothness parameters $(\nu_i)_{1 \leq i \leq m}$.
We first give a general rule which holds for any arbitrary sampling.
Recalling  Theorem~\ref{thm:stepsizes}\ref{thm:stepsizes_iii} and
Remark~\ref{rmk:20190128e}\ref{rmk:20190128e_i} and noting that
$\{k \vert i \in I_k\} = \{k \vert i \in \mathrm{spt}(\ensuremath \as_k)\} = \mathrm{spt}(\ensuremath \as^i)$, if we set
$\nu_i = \sum_{k\in \mathrm{spt}(\ensuremath \as^i)}
\norm{\ensuremath \as_k}^2$, then \ref{eq:S0} and hence \ref{eq:S2} holds.
This choice was considered in \cite{Nec16}.
Moreover, according to the discussion at the beginning of Section~\ref{subsec:beta}
other options for satisfying \ref{eq:S2} are
$\nu_i=\min\{\eta,\tau_{\max}\} \norm{\ensuremath \as^i}^2$ or
 $\nu_i = \sum_{k \in \mathrm{spt}(\ensuremath \as^i)} \min\{\mathrm{card}(\mathrm{spt}(\ensuremath \as_k)), \tau_{\max}\} (\ensuremath \as_k^i)^2$. This latter choice is better than the second one and,
if we assume that the nonzero entries of $A$ are about of the same magnitude,
it is also better than the first one.
Next, we face the special case of the $\tau$-nice sampling
which allows to reduce the $\nu_i$'s while satisfying \ref{eq:S1}.
Recalling the corresponding discussion in Section~\ref{subsec:beta},
we have the following alternatives: (1) set, for every $i \in [m]$,
$\nu_i = (1 + (\tau-1)(\eta-1)/(m-1)) \norm{\ensuremath \as^i}^2$;
(2) set for every $i \in [m]$, $\nu_i =
\sum_{k\in \mathrm{spt}(\ensuremath \as^i)} (1 + (\tau-1)(\mathrm{card}(\mathrm{spt}(\ensuremath \as_k))-1)/(m-1)) (\ensuremath \as_k^i)^2$.
Finally, we make few remarks on the convergence properties of algorithm~\eqref{eq:algoLasso}.
Since the objective function in \eqref{eq:20181015a} satisfies a quadratic growth condition on its sublevel sets \cite[Example~3.8]{Gar17}, then Remark~\ref{rmk:20190129e}, Remark~\ref{rmk:20200331a},
and Theorem~\ref{thm:20181214a}, yield linear convergence of algorithm
\eqref{eq:algoLasso} provided that \ref{eq:S2} holds.
Whereas Theorem~\ref{thm20200218a}
ensures that if we modify algorithm \eqref{eq:algoLasso}
so that we accept the next iterate $\ensuremath x^{\ensuremath n+1}$ only if $\norm{u^{\ensuremath n+1}}^2
+ 2\lambda \norm{\ensuremath x^{\ensuremath n+1}}_1\leq \norm{u^\ensuremath n}^2 + 2\lambda \norm{\ensuremath x^{\ensuremath n}}_1$,
then the resulting algorithm
converges linearly under condition \ref{eq:S1}.
If the violation of the monotonicity condition above
occurs few times along all the iterations, this modification does not increase much
the computational cost of the algorithm (see also Section~\ref{sec:experiments}).
We stress that both
Theorem~\ref{thm:20181214a} and Theorem~\ref{thm20200218a}
ensure also almost sure and linear convergence in mean of the iterates
of \eqref{eq:20170915a}.
This latter result is new and is especially relevant in this context, since
 the iterates carry sparsity information.

\subsection{Computing the minimal norm solution of a linear system}
Let $\ensuremath{{\As}} \in \ensuremath \mathbb{R}^{m\times p}$ and $\ensuremath \bs \in R(\ensuremath{{\As}})$ (the range of $\ensuremath{{\As}}$).
Let us consider the problem
\begin{equation}
\label{eq:20170920a}
\minimize{\substack{\ensuremath \xx \in \ensuremath \mathbb{R}^p \\ \ensuremath{{\As}} \ensuremath \xx = \ensuremath \bs}}{\frac 1 2 \norm{\ensuremath \xx}^2}.
\end{equation}
Here, we denote by $\ensuremath \as_i \in \ensuremath \mathbb{R}^p$ and $\ensuremath \as^k \in \ensuremath \mathbb{R}^m$
the $i$-th row and the $k$-th column of $\ensuremath{{\As}}$.
The dual problem is
\begin{equation}
\label{eq:20190118h}
\minimize{\ensuremath \uu \in \ensuremath \mathbb{R}^m}{\frac 1 2 \norm{ \ensuremath{{\As}}^\top\ensuremath \uu}^2 - \scalarp{\ensuremath \uu,\ensuremath \bs}}:=\mathcal{D}(\ensuremath \uu),
\end{equation}
which is a smooth convex optimization problem.
Moreover, if $\ensuremath \xx_*$ is the solution of \eqref{eq:20170920a} and
$\ensuremath \xx = \ensuremath{{\As}}^\top \ensuremath \uu$ (the primal-dual relationship), then we have
\begin{equation}
\label{eq:20190118j}
\frac 1 2 \norm{\ensuremath \xx - \ensuremath \xx_*}^2
\leq \mathcal{D}(\ensuremath \uu) - \inf \mathcal{D}.
\end{equation}
Then, the dual problem is clearly of the form \eqref{eq:20190110b},
with $\ensuremath{{\gs}}=0$,
and \ref{eq:B0} and \ref{eq:B1}
are satisfied, assuming that each block is made of
one coordinate only, with $L_i = \norm{\ensuremath \as_i}^2$ and $\eta = \max_{1 \leq k \leq p} \mathrm{card}(\mathrm{spt}(\ensuremath \as^k))$. So, Algorithm~\ref{algoRCD} applied to \eqref{eq:20190118h}, turns into
\begin{equation*}
\ensuremath u^{\ensuremath n+1} = \ensuremath u^\ensuremath n - \sum_{i=1}^m \varepsilon_i^\ensuremath n\gamma_{i}
\big(\scalarp{\ensuremath \as_{i},
\ensuremath{{\As}}^\top \ensuremath u_\ensuremath n} - \ensuremath \mathsf{b}_{i}\big)
\ensuremath \ee_{i}.
\end{equation*}
Now, setting $\ensuremath x^\ensuremath n = \ensuremath{{\As}}^\top \ensuremath u^\ensuremath n$ and multiplying the above equality by $\ensuremath{{\As}}^\top$, we have
\begin{equation}
\label{eq:20190129f}
\ensuremath x^{\ensuremath n+1} = \ensuremath x^\ensuremath n - \sum_{i=1}^m
\varepsilon_i^\ensuremath n\gamma_{i} \big(\scalarp{\ensuremath \as_{i}, \ensuremath x^\ensuremath n} - \ensuremath \mathsf{b}_{i}\big)
\ensuremath \as_{i}.
\end{equation}
Since, $\ensuremath \bs \in R(\ensuremath{{\As}})$, it is easy to see, through a singular value decomposition of $\ensuremath{{\As}}$, that, for every $\ensuremath \uu \in \ensuremath \mathbb{R}^m$,
$ \sigma^2_{\min}(\ensuremath{{\As}})\,\mathrm{dist}(\ensuremath \uu, \argmin \mathcal{D})
\leq \norm{\nabla \mathcal{D}(\ensuremath \uu)}$ (where $\sigma_{\min}(\ensuremath{{\As}})$ is the minimum  singular value of $\ensuremath{{\As}}$) \cite[Example 3.6]{Gar17}.
So, in view of Remark~\ref{rmk:20190603a}\ref{rmk:20190129e_00}-\ref{rmk:20190129e_0},
\ref{eq:EB1} is satisfied on the entire space
with constant $c_{\ensuremath \mathbb{R}^m, \ensuremath{{\Gammas}}^{-1}} = (\gamma_{\min}
\sigma^2_{\min}(\ensuremath{{\As}}))^{-1}$. Therefore, if, for every $i \in [m]$, $\gamma_i <2/(\beta_{1,i}\norm{\ensuremath \as_i}^2)$, Theorem~\ref{thm:20181214a} and \eqref{eq:20190118j}
ensure the linear convergence of the iterates $\ensuremath x^n$'s towards the solution of \eqref{eq:20170920a} with rate
$
\big(1 - \ensuremath \mathsf{p}_{\min}\min\big\{1, \gamma_{\min}\sigma^2_{\min}(\ensuremath{{\As}})(2-\delta)/2\big\} \big)^{1/2}
$.
We remark that \eqref{eq:20190129f} is nothing but a stochastic
gradient descent algorithm on the problem
\begin{equation*}
\minimize{\ensuremath \xx \in \ensuremath \mathbb{R}^p}{\frac 1 2 \norm{\ensuremath{{\As}} \ensuremath \xx - \ensuremath \bs}^2}=
\frac 1 2 \sum_{i=1}^m (\scalarp{\ensuremath \as_i, \ensuremath \xx} - \ensuremath \mathsf{b}_i)^2.
\end{equation*}
Since $\norm{\ensuremath{{\As}} \ensuremath \xx - \ensuremath \bs}^2 \leq \norm{\ensuremath{{\As}}}^2 \norm{\ensuremath \xx - \ensuremath \xx_*}^2$,
we have then showed the linear convergence rate
\begin{equation*}
\frac1 2 \norm{\ensuremath{{\As}} \ensuremath x^n - \ensuremath \bs}^2 - \frac 1 2 \norm{\ensuremath{{\As}} \ensuremath \xx_* - \ensuremath \bs}^2 = O\bigg(\sigma^2_{\max}(\ensuremath{{\As}})
\bigg(1 - \ensuremath \mathsf{p}_{\min}\min\bigg\{1, \frac{\gamma_{\min} \sigma^2_{\min}(\ensuremath{{\As}})(2-\delta)}{2}\bigg\} \bigg)^n \bigg),
\end{equation*}
of the stochastic gradient descent with arbitrary
and possibly variable batch size for least squares problems.
We finally note that in the serial case, that is, if for every $\ensuremath n \in \ensuremath \mathbb{N}$
$\mathrm{spt}(\ensuremath \varepsilon^\ensuremath n) = \{ i_\ensuremath n \}$, then multiplying equation \eqref{eq:20190129f} by $\ensuremath \as_{i_\ensuremath n}^\top$,
we have
\begin{equation*}
\scalarp{\ensuremath \as_{i_\ensuremath n}, \ensuremath x^{\ensuremath n+1}} =
\scalarp{\ensuremath \as_{i_\ensuremath n}, \ensuremath x^\ensuremath n} - \gamma_{i_\ensuremath n}
\big(\scalarp{\ensuremath \as_{i_\ensuremath n}, \ensuremath x^\ensuremath n} - \ensuremath \mathsf{b}_{i_\ensuremath n}\big)
\norm{\ensuremath \as_{i_\ensuremath n}}^2.
\end{equation*}
Therefore, since in this case $\beta_{1,i}= \beta_2=1$, we can chose the stepsizes such that $\gamma_{i}\norm{\ensuremath \as_{i}}^2 = 1$ (so that $\delta=1$) and hence
 $\ensuremath x^{\ensuremath n+1}$ is a solution of the $i_\ensuremath n$-th equation of the linear system
 $\ensuremath{{\As}} \ensuremath \xx = \ensuremath \bs$.
Moreover, $\ensuremath x^{\ensuremath n+1}$ is the projection of $\ensuremath x^\ensuremath n$ onto the affine space defined
by the equation $\ensuremath \as_{i_\ensuremath n} \ensuremath \xx = \ensuremath \mathsf{b}_{i_\ensuremath n}$ \cite{Wri15}.
Thus, this method is nothing but the \emph{randomized Kaczmarz method} \cite{Stro09}  and we proved linear convergence for general probabilities $\ensuremath \mathsf{p}_i$'s, although the
constants we derive are not optimal (see \cite{Lev10,Stro09,Wri15}).

\subsection{Regularized empirical risk minimization}

Let $\ensuremath \mathsf{H}$ be a separable real Hilbert space.
Regularized empirical risk estimation solves the following optimization problem
\begin{equation}
\label{eq:20170217q}
\minimize{\ensuremath \mathsf{w} \in \ensuremath \mathsf{H}}{\frac{1}{\lambda m}
\sum_{i=1}^m \ell(\ensuremath \mathsf{y}_i, \scalarp{\ensuremath \mathsf{w},\ensuremath \mathsf{x}_i}) + \frac 1 2 \norm{\ensuremath \mathsf{w}}^2}
:= \mathcal{P}(\ensuremath \mathsf{w}),
\end{equation}
where $(\ensuremath \mathsf{x}_i, \ensuremath \mathsf{y}_i)_{1 \leq i \leq m}$ is the training set (input-output pairs),
 $\ell\colon \ensuremath \mathcal{Y}\times\ensuremath \mathbb{R} \to \ensuremath \mathbb{R}_+$, $\ensuremath \mathcal{Y} \subset \ensuremath \mathbb{R}$, is the \emph{loss} function, which is convex in the second variable,
 and $\lambda>0$
is a regularization parameter.
The dual problem of \eqref{eq:20170217q} is
\begin{equation}
\label{eq:gendual}
\minimize{\ensuremath \uu \in \ensuremath \mathbb{R}^m}{\frac 1 2 \ensuremath \uu^\top \mathsf{K} \ensuremath \uu
+ \frac{1}{\lambda m} \sum_{i=1}^m \ell^* (\ensuremath \mathsf{y}_i, - \ensuremath \mathsf{u}_i \lambda m)}
:=\mathcal{D}(\ensuremath \uu),
\end{equation}
where $\ell^*(\ensuremath \mathsf{y}_i, \cdot)$ is the Fenchel conjugate of $\ell(\ensuremath \mathsf{y}_i, \cdot)$ and
$\mathsf{K} = \ensuremath{{\XX}} \ensuremath{{\XX}}^\top \in \ensuremath \mathbb{R}^{m\times m}$  is the Gram matrix
of $(\ensuremath \mathsf{x}_i)_{1 \leq i \leq m}$.
Moreover, the solutions $(\bar{\ensuremath \mathsf{w}}, \bar{\ensuremath \uu})$
of the primal and dual problems are characterized by the following KKT conditions
\begin{equation}
\label{eq:kkt}
\begin{cases}
\displaystyle\bar{\ensuremath \mathsf{w}} = \ensuremath{{\XX}}^\top \bar{\ensuremath \uu}
= \sum\limits_{i=1}^m \bar{\ensuremath \mathsf{u}}_i \ensuremath \mathsf{x}_i,\\[2.5ex]
\forall\, i \in \{1,\dots, m\}\ \
- \bar{\ensuremath \mathsf{u}}_i m \lambda \in \partial \ell(\ensuremath \mathsf{y}_i, \scalarp{\ensuremath \mathsf{x}_i, \bar{\ensuremath \mathsf{w}}}),
\end{cases}
\end{equation}
where $\partial \ell(\ensuremath \mathsf{y}_i,\cdot)$ is the subdifferential of $\ell(\ensuremath \mathsf{y}_i, \cdot)$.
Note also that the first of \eqref{eq:kkt} gives the link between the dual and the primal variable
and, if $\ensuremath \mathsf{w} = \ensuremath{{\XX}}^\top \ensuremath \uu$, then it holds $(1/2) \norm{\ensuremath \mathsf{w} - \bar{\ensuremath \mathsf{w}}}^2 \leq \mathcal{D}(\ensuremath \uu) - \inf \mathcal{D}$.
Now, the dual problem \eqref{eq:gendual} is of the form \eqref{eq:20190110b}
and hence Algorithm~\ref{algoRCD} can be applied. The following examples give implementation details for two specific losses.

\begin{example}[Ridge regression]\
The least squares loss is  $\ell(s,t) = (1/2) \abs{s - t}^2$.
Then $\ell^*(s,r) = (1/2) r^2 + r s$ and,
in this case, \eqref{eq:gendual} reduces to
\begin{equation*}
\minimize{\ensuremath \uu \in \ensuremath \mathbb{R}^m}
{\frac 1 2 \ensuremath \uu^\top ( \mathsf{K} + \lambda m \ensuremath{\mathsf{Id}} )\ensuremath \uu
- \ensuremath \yy^\top \ensuremath \uu} :=\mathcal{D}(\ensuremath \uu)
\end{equation*}
which  is strongly convex with modulus $\lambda m$ and has  solution $\bar{\ensuremath \uu} = (\mathsf{K} + m \lambda \ensuremath{\mathsf{Id}})^{-1} \ensuremath \yy$.
Since $\mathcal{D}$ is smooth and
$\nabla_i \mathcal{D} (\ensuremath \uu) = \ensuremath \ee_i^\top (\mathsf{K} + \lambda m\ensuremath{\mathsf{Id}}) \ensuremath \uu - \ensuremath \mathsf{y}_i$,  conditions \ref{eq:B0} and \ref{eq:B1} hold with $L_i = \norm{\ensuremath \mathsf{x}_i}^2 + \lambda m = \mathsf{K}_{i,i} + \lambda m$ and
 Algorithm~\ref{algoRCD} (with $\ensuremath{{\gs}}=0$) becomes
 \begin{equation}
\label{eq:20190123b}\ensuremath u^{\ensuremath n+1}
  = \ensuremath u^\ensuremath n - \sum_{i=1}^m \varepsilon^\ensuremath n_i \gamma_i
\big( \ensuremath \ee_i^\top\mathsf{K}\ensuremath u^\ensuremath n + \lambda m u^\ensuremath n_i
- \ensuremath \mathsf{y}_i \big) \ensuremath \ee_i.
\end{equation}
Moreover, multiplying
\eqref{eq:20190123b} by $\ensuremath{{\XX}}^\top$, defining
$w^\ensuremath n = \ensuremath{{\XX}}^\top \ensuremath u^\ensuremath n$,
and recalling that
$\mathsf{K} = \ensuremath{{\XX}} \ensuremath{{\XX}}^\top$, we have
\begin{align}
\label{eq:20190131a}
w^{\ensuremath n+1}
&= w^\ensuremath n - \sum_{i=1}^m \varepsilon^\ensuremath n_i \gamma_i
\big(\scalarp{w^\ensuremath n, \ensuremath \mathsf{x}_i}  - \ensuremath \mathsf{y}_i\big) \ensuremath \mathsf{x}_i
-  \lambda m \sum_{i=1}^m \varepsilon^\ensuremath n_i \gamma_i  \ensuremath u^\ensuremath n_i \ensuremath \mathsf{x}_i.
\end{align}
Note that, since the dual problem is strongly convex with modulus $\lambda m$,
then it follows from Theorem~\ref{p:20181130a}\ref{p:20181130a_iii}, Remark~\ref{rmk:20190130a}, and Theorem~\ref{thm:dualgap}\ref{thm:dualgap_i}
that, setting, for every $i \in [m]$,
$\nu_i \geq \beta_{1,i} (\mathsf{K}_{i,i} + \lambda m)$
and $\gamma_i = 1/\nu_i$, we have
\begin{equation*}
\ensuremath{\mathsf E}[\mathcal{P}(w^\ensuremath n)] - \inf \mathcal{P} \leq
\bigg(1 + \frac{\norm{\mathsf{K}}}{\lambda m } \bigg)
\bigg( 1 - \ensuremath \mathsf{p}_{\min} \frac{2 \lambda m}{\nu_{\max} + \lambda m} \bigg)^\ensuremath n
\mathrm{const.}
\end{equation*}
Now, we compare algorithm \eqref{eq:20190131a} with the stochastic gradient
descent on problem \eqref{eq:20170217q}.
Assume that $\ensuremath{\mathsf P}(\sum_{i=1}^m \varepsilon^\ensuremath n_i = \tau) = 1$ for some
$\tau \in [m]$. Then,
$\ensuremath \mathsf{p}_{\min} = \tau/m$ and
we can take
$\zeta_i \leq (\mathsf{K}_{i,i} + \lambda m)^{-1}$,
and set
$\nu_i = \tau/\zeta_i$ and $\gamma_i = 1/\nu_i$, so that
 algorithm \eqref{eq:20190131a} turns into
\begin{equation}
\label{eq:20190123d}
w^{\ensuremath n+1} =w^\ensuremath n - \sum_{i=1}^m \varepsilon^\ensuremath n_i \frac{\zeta_i}{\tau}
\big(\scalarp{w^\ensuremath n, \ensuremath \mathsf{x}_i}  - \ensuremath \mathsf{y}_i\big) \ensuremath \mathsf{x}_i
-  \lambda m \sum_{i=1}^m \varepsilon^\ensuremath n_i  \frac{\zeta_i}{\tau} u^\ensuremath n_i \ensuremath \mathsf{x}_i.
\end{equation}
If we apply stochastic gradient descent with batch size $\tau \in [m]$
and stepsize $\zeta>0$ directly on the primal problem \eqref{eq:20170217q} (multiplied by $\lambda m$), and recalling that $w^\ensuremath n = \sum_{i=1}^m u^\ensuremath n_i \ensuremath \mathsf{x}_i$,
 we have
\begin{align}
 w^{\ensuremath n + 1}
\label{eq:20190123e}= w^\ensuremath n - \frac{\zeta}{\tau} \sum_{i=1}^m \varepsilon^\ensuremath n_i
( \scalarp{w^\ensuremath n, \ensuremath \mathsf{x}_i} - \ensuremath \mathsf{y}_i ) \ensuremath \mathsf{x}_i -  \lambda m
\frac{\zeta}{m} \sum_{i=1}^m u_i^\ensuremath n \ensuremath \mathsf{x}_i.
\end{align}
Then, comparing \eqref{eq:20190123d} and \eqref{eq:20190123e} we see that,
provided that $\zeta_i = \zeta$ for every $i \in [m]$,
they only differ for the replacement $(1/m) \sum_{i=1}^m u^\ensuremath n_i \ensuremath \mathsf{x}_i\ \leftrightarrow\
(1/\tau) \sum_{i=1}^m \varepsilon^\ensuremath n_i u^\ensuremath n_i \ensuremath \mathsf{x}_i$.
We stress that the stepsize $\zeta$ in the stochastic gradient descent
algorithm \eqref{eq:20190123e} is normally
set according to the spectral norm of
$\mathsf{K} + \lambda m \ensuremath{\mathsf{Id}}$, which may be difficult to compute.
On the contrary in algorithm \eqref{eq:20190123d} the stepsizes $\zeta_i$'s
are simply set as $\zeta_i \leq 1/(\mathsf{K}_{i,i} + \lambda m)$,
so they allow possibly much longer steps and also do not require any SVD computation.
\end{example}

\begin{example}[Support vector machines]
The \emph{hinge loss} is
$\ell(s,t) = (1-s t)_+$. Then we have $\ell^*(s,r) = r + \iota_{[0,1]}(s r)$
and the dual problem \eqref{eq:gendual}
is
\begin{equation}
\label{eq:20190131b}
\minimize{\ensuremath \uu \in \ensuremath \mathbb{R}^m}
{\frac 1 2 \ensuremath \uu^\top \mathsf{K} \ensuremath \uu
 - \ensuremath \yy^\top\ensuremath \uu + \iota_{(\lambda m)^{-1}[0,1]^m}(\ensuremath \yy \odot \ensuremath \uu)}.
\end{equation}
Then Algorithm~\ref{algoRCD} on the dual turns into a parallel random block-coordinate
projected gradient descent method.
Moreover, it follows from Remark~\ref{rmk:20181219e}\ref{rmk:20181219e_iii}
that the objective in \eqref{eq:20190131b}
satisfies \ref{eq:EB1} on its domain.
Therefore, it follows from
Theorem~\ref{thm:20181214a}, Theorem~\ref{thm:stepsizes}\ref{thm:stepsizes_i}, and
Theorem~\ref{thm:dualgap}\ref{thm:dualgap_ii}
that
$\ensuremath{\mathsf E}[\mathcal{P}(w^\ensuremath n)]- \inf \mathcal{P}$ converges linearly to zero,
provided that, for all $i \in [m]$, $\nu_i \geq \beta_{1,i} \mathsf{K}_{ii}$ and $\gamma_i < 2/\nu_i$.
\end{example}

\section{Numerical Experiments}
\label{sec:experiments}
In this section we consider a Lasso problem,
that is,
\begin{equation}
\label{eq:20200312a}
\min_{\ensuremath \xx \in \ensuremath \mathbb{R}^m} \frac 1 2 \norm{\ensuremath{{\As}} \ensuremath \xx - \ensuremath \bs}_2^2 + \lambda \norm{\ensuremath \xx}_1 \quad(\lambda>0),
\end{equation}
where $\ensuremath{{\As}} \in \ensuremath \mathbb{R}^{p\times m}$ is generated with random entries uniformly distributed in $[-1,1]$ so that
each row is sparse and $\ensuremath \bs = \ensuremath{{\As}} \bar{\ensuremath \xx} + 0.06\cdot \alpha$ with $\alpha \sim N(0,\ensuremath{\mathsf{Id}}_p)$ and $\bar{x}$
a sparse vector in $\ensuremath \mathbb{R}^m$.
We implement Algorithm~\ref{algoRCD} with $\gamma_i = \delta/\nu_i$ ($0<\delta<2$) and a $\tau$-nice uniform sampling, as described
in Section~\ref{subsec:lasso}.
We present two experiments. The first compares conditions
\ref{eq:S1} and \ref{eq:S2} for the determination of the stepsizes $\gamma_i$. The second one investigates the role played by $\delta$.
In all the experiments we empirically checked that the algorithm is essentially descending
in the sense that during the iterations there are very few violations of the descent property and with low magnitude.
So, since the objective function in \eqref{eq:20200312a} satisfies \ref{eq:EB1} on the sublevel sets,
in virtue of Theorem~\ref{thm20200218a}, linear convergence holds.

\subsection*{Condition~\ref{eq:S1} vs \ref{eq:S2} and the effectiveness of the parallel strategy.}
We compare the conditions \ref{eq:S1} and \ref{eq:S2} for the stepsizes selection
and we checked the critical role played by \ref{eq:S1}
for the effectiveness of the parallel strategy on problems with sparse structure.
Here we set $\delta=1$. In Figure~\ref{lasso_stepsizes}, Algo1 uses smoothness parameters specifically
designed for the $\tau$-nice sampling, that is, $\nu_i=\beta_1 \norm{\ensuremath \as^i}^2$ with
$\beta_1 = 1 + (\tau-1)(\eta-1)/(m-1)$ (making \ref{eq:S1} satisfied), while Algo2 uses a more conservative choice for the smoothness parameters which is valid for any sampling updating a maximum of $\tau$ blocks per iteration, that is $\nu_i = \beta_2 \norm{\ensuremath \as^i}^2$ with $\beta_2 = \min\{\tau, \eta\}$ (making also \ref{eq:S2} satisfied).
\begin{figure}[t]
\begin{center}
\vspace{-2ex}
\includegraphics[width=0.5\textwidth]{hugelasso1_S1_S2.eps}
\hspace{-3ex}
\includegraphics[width=0.5\textwidth]{lasso_S1_S2.eps}
\end{center}
\vspace{-2em}
\caption{\small Comparison between \ref{eq:S1} and \ref{eq:S2} for the the stepsizes selection in a Lasso problem.
Left: $5\cdot 10^4$ equations in $10^5$ unknowns; degree of partial separability $\eta=148$.
Right: $10^3$ equations in $5\cdot 10^3$ unknowns;
degree of partial separability $\eta = 563$.}
\label{lasso_stepsizes}
\end{figure}
\begin{figure}[h!]
\begin{center}
\vspace{-1ex}
\includegraphics[width=0.5\textwidth]{lasso_delta_eta71.eps}
\hspace{-3ex}
\includegraphics[width=0.5\textwidth]{lasso_delta_eta563_1.eps}
\includegraphics[width=0.5\textwidth]{lasso_delta_eta563_2.eps}
\hspace{-3ex}
\includegraphics[width=0.5\textwidth]{lasso_delta4.eps}
\end{center}
\vspace{-2em}
\caption{\small Comparison between $\delta=1$ and $\delta>1$.
Lasso problem with $n=10^3$ equations in $m= 5\cdot 10^3$ unknowns.
Degree of partial separability $\eta=71$ (top left), $\eta= 563$ (top right and bottom left), and $\eta = 2594$ (bottom right). The choice $\delta=1$ is better than $\delta>1$ only for $\eta=71$ and $\tau=5,10$.}
\label{lasso_overrelaxed2}
\end{figure}
In the left diagram, we considered a large scale setting with $\eta\ll m$. In that case
$\beta_1$ may be much smaller than $\beta_2$ leading to significantly larger stepsizes. Moreover, and more importantly, we note that as long as $\tau$ is small enough the behavior of Algo1 does not depend on $\tau$ (indeed $\tau=10,50,100$ perform equally well), whereas this is not true for Algo2.
This feature of \ref{eq:S1}, first noted in \cite{Ric16}, has been already discussed after
Theorem~\ref{thm:20171207a} and  is at the basis
of the effectiveness of the parallel strategy.
Indeed in the small-$\tau$ regime described above,
the various versions of Algo1 depicted in Figure~\ref{lasso_stepsizes}(left)
have the same total computation cost ($\bar{n} m$ block-coordinate updates),
but the parallel implementation on $\tau$ cores is $\tau$ times faster than the serial one ($\tau=1$).
Finally, in the right diagram of Figure~\ref{lasso_stepsizes} we show a scenario in which
$\eta/m$ is larger (the problem is less sparse).
In such situation we see that the difference between the two
stepsize selection criteria is less evident for $\tau \geq 50$.
Moreover, Algo1 is more sensitive to $\tau$ (compare $\tau=10$ and $\tau= 50$),
so that the benefit of the parallel scheme is reduced.
\subsection*{The effect of $\delta>1$}
Here we study the effect of over-relaxing the stepsizes, meaning choosing $\delta>1$.
We compare Algorithm~\ref{algoRCD} with $\delta=1$ and several choices of $\delta>1$.
Figure~\ref{lasso_overrelaxed2} considers different scenarios for the degree of separability $\eta$ of $f$.
In those cases we see that choosing $\delta>1$ usually speeds up the convergence,
depending on the parameter $\eta$ of partial separability of $f$
and the number $\tau$ of parallel block updates.
This fact seems not to occur when both $\eta/m$ and $\tau/m$ are very small.