11institutetext:  Correspondence to Tianjiao Li. Coauthors are listed
according to the alphabetic order.
1 LJK, University Grenoble Alpes, 38401 Domaine Universitaire de Saint-
Martin-d’Hères, France
Email<EMAIL_ADDRESS>anatoli.juditsky@univ-
grenoble-alpes.fr
2 H. Milton Stewart School of Industrial and Systems Engineering, Georgia
Institute of Technology, Atlanta, GA, USA
Email<EMAIL_ADDRESS><EMAIL_ADDRESS>
# Accelerated stochastic approximation with state-dependent noise

Sasila Ilandarideva1 Anatoli Juditsky1 Guanghui Lan2 Tianjiao Li2 SI and AJ
were partially supported by Multidisciplinary Institute in Artificial
intelligence MIAI @ Grenoble Alpes ANR-19-P3IA-0003. GL and TL were partially
supported by Division of Mathematical Sciences grant DMS-1953199 and Air Force
Office of Scientific Research grant FA9550-22-1-0447.

(the date of receipt and acceptance should be inserted later)

###### Abstract

We consider a class of stochastic smooth convex optimization problems under
rather general assumptions on the noise in the stochastic gradient
observation. As opposed to the classical problem setting in which the variance
of noise is assumed to be uniformly bounded, herein we assume that the
variance of stochastic gradients is related to the “sub-optimality” of the
approximate solutions delivered by the algorithm. Such problems naturally
arise in a variety of applications, in particular, in the well-known
generalized linear regression problem in statistics. However, to the best of
our knowledge, none of the existing stochastic approximation algorithms for
solving this class of problems attain optimality in terms of the dependence on
accuracy, problem parameters, and mini-batch size.

We discuss two non-Euclidean accelerated stochastic approximation
routines—stochastic accelerated gradient descent (SAGD) and stochastic
gradient extrapolation (SGE)—which carry a particular duality relationship. We
show that both SAGD and SGE, under appropriate conditions, achieve the optimal
convergence rate, attaining the optimal iteration and sample complexities
simultaneously. However, corresponding assumptions for the SGE algorithm are
more general; they allow, for instance, for efficient application of the SGE
to statistical estimation problems under heavy tail noises and discontinuous
score functions. We also discuss the application of the SGE to problems
satisfying quadratic growth conditions, and show how it can be used to recover
sparse solutions. Finally, we report on some simulation experiments to
illustrate numerical performance of our proposed algorithms in high-
dimensional settings.

###### Keywords:

stochastic optimization state-dependent noise accelerated stochastic
approximation stochastic gradient extrapolation sparse recovery

###### MSC:

90C15 90C25 62L20 68Q25 62J12

##  1 Introduction

This paper focuses on the stochastic optimization problem given by

$f^{*}:=\min\limits_{x\in X}f(x)$ (1.1)

where $X$ is a closed convex subset of a Euclidean space $E$ and
$f:X\to\mathbb{R}$ is a smooth convex function with Lipschitz continuous
gradient, i.e., for some $L\geq 0$,

$0\leq f(y)-f(x)-\langle\nabla f(x),y-x\rangle\leq\tfrac{L}{2}\|y-x\|^{2},\ \
\ \forall x,y\in X.$ (1.2)

We assume throughout the paper that the set of optimal solutions $X^{*}$ is
nonempty.

We consider the stochastic setting where only stochastic first-order
information about $f$ is available for solving problem (1.1). Specifically, at
the current search point $x_{t}\in X$, a stochastic oracle (SO) generates the
stochastic operator ${\cal G}(x_{t},\xi_{t})$, where $\xi_{t}\in\Xi$ denotes a
random variable, whose probability distribution is supported on a set $\Xi$.
We assume $\xi_{t}$ is independent of $x_{0},...,x_{t}$, and
$\\{\xi_{t}\\}_{t\geq 0}$ are mutually independent. We assume that ${\cal
G}(x_{t},\xi_{t})$ is an unbiased estimator of $g(x_{t})=\nabla f(x_{t})$
satisfying

$\displaystyle{\mathbb{E}}_{\xi_{t}}[{\cal G}(x_{t},\xi_{t})]=g(x_{t})$ (1.3)

(expectation w.r.t. to the distribution $\xi_{t}$).

Stochastic approximation (SA) and stochastic mirror descent (SMD) methods are
routinely used to solve stochastic optimization problems; see, e.g.,
nemirovskii1979complexity ; polyak1990new ; polyak1992acceleration ;
nemirovski2009robust . More specifically, it was shown in nemirovski2009robust
that SMD can achieve the optimal sample complexity for general nonsmooth
optimization and saddle point problems. For smooth stochastic optimization
problems, Lan lan2012optimal ; LanBook2020  introduced an accelerated
stochastic approximation (AC-SA), also known as stochastic accelerated
gradient descent (SAGD), which was obtained by replacing exact gradients with
their unbiased estimators in the celebrated accelerated gradient methods
nesterov1983method  (see also YudinNemirovski77 ; NemirovskiYudin83-smooth ;
nemyud:83 ; NemirovskiNester85  for early developments). It was shown in
lan2012optimal  that AC-SA achieves the optimal sample complexity for smooth,
nonsmooth and stochastic convex optimization (see ghadimi2012optimal ;
ghadimi2013optimal  for generalization to the strongly convex settings). It
should be noted that the original analysis of AC-SA in lan2012optimal  was
carried out under uniformly bounded variance condition

${\mathbb{E}}_{\xi_{t}}[\|{\cal
G}(x_{t},\xi_{t})-g(x_{t})\|_{*}^{2}]\leq\sigma^{2}$ (1.4)

where $\|\cdot\|_{*}$ is the norm conjugate to $\|\cdot\|$. However, it has
been observed recently (see, e.g., srebro2010smoothness ; juditsky2023sparse )
that this uniformly bounded variance condition is not necessarily satisfied in
some important applications in which the variance of the stochastic gradient
depends on the search point $x_{t}$. As a motivation, consider the fundamental
to Statistical Learning problem of parameter estimation in the generalized
linear regression (GLR) model in which one aims to estimate an unknown
parameter vector $x^{*}\in X\subset\mathbb{R}^{n}$ given observations
$(\phi_{t},\eta_{t})$,

$\displaystyle\eta_{t}=u(\phi^{T}_{t}x^{*})+\zeta_{t},\quad t=1,2,...,$ (1.5)

where, in generalized linear models terminology, $\eta_{t}\in\mathbb{R}$ are
responses, $\phi_{t}\in\mathbb{R}^{n}$ are random regressors,
$\zeta_{t}\in\mathbb{R}$ are zero-mean random noises which are assumed to be
mutually independent and independent of $\phi_{t}$, and
$u:\mathbb{R}\rightarrow\mathbb{R}$ is the (generally nonlinear) “activation
function”. Then it follows directly that

${\mathbb{E}}[\phi_{t}(u(\phi^{T}_{t}x^{*})-\eta_{t})]={\mathbb{E}}[\phi_{t}\,\zeta_{t}]=0.$
(1.6)

Thus, the problem of recovery of $x^{*}$ from observations $\eta_{t}$ and
$\phi_{t}$ may be formulated as a stochastic optimization problem.
Specifically, when denoting $v:\mathbb{R}\to\mathbb{R}$ the primitive of $u$,
i.e., $v^{\prime}(t)=u(t)$ and assuming that $x^{*}\in\hbox{\rm int\,}X$,
(1.6) may be seen as as the optimality condition for the problem

$\min_{x\in
X}\left\\{f(x):={\mathbb{E}}[v(\phi^{T}x)-\phi^{T}x\eta]\right\\}.$ (1.7)

Clearly, the gradient of $f$ is given by
$g(x)={\mathbb{E}}[\phi(u(\phi^{T}x)-\eta)]$, and one of its unbiased
stochastic estimator is ${\cal
G}(x,\overbrace{(\phi,\zeta)}^{=:\xi})=\phi(u(\phi^{T}x)-\eta)$. Under mild
assumptions, one can show (see Section 6.1) that the noise

${\cal
G}(x_{t},(\phi_{t},\zeta_{t}))-g(x_{t})=\phi_{t}[u(\phi_{t}^{T}x_{t})-u(\phi_{t}^{T}x^{*})]-\phi_{t}\zeta_{t}+g(x^{*})-g(x_{t})$

of the gradient observation ${\cal G}(x_{t},(\phi_{t},\zeta_{t}))$ at $x_{t}$
satisfies the condition

${\mathbb{E}}_{\xi_{t}}\big{[}\|{\cal
G}(x_{t},\xi_{t})-g(x_{t})\|_{*}^{2}\big{]}\leq\sigma_{t}^{2}(x_{t})={\cal
L}[f(x_{t})-f^{*}-\langle g(x^{*}),x_{t}-x^{*}\rangle]+\sigma_{*}^{2}$ (SN)

for some ${\cal L},\sigma_{*}>0$. In what follows, with some terminology
abuse, we refer to (SN) and similar conditions as state-dependent noise
assumptions. More generally, if compared to the “uniform noise” condition
(1.4), (SN) may be seen as a refined assumption on the structure of the
stochastic oracle. Furthermore, one can easily verify that in various settings
of the GLR problem, condition (SN) holds, while condition (1.4) of uniformly
bounded noise variance is violated (e.g., in the case of unbounded $X$ and
linear $u$). Other stochastic optimization problems in which similar state-
dependent noise assumptions apply can also be found in recent literature on
machine learning woodworth2021even ; juditsky2023sparse ;
ilandarideva2023stochastic  and reinforcement learning, see, e.g.,
kotsalis2022simple1 ; kotsalis2022simple ; li2021accelerated ;
li2022stochastic .

Motivated by these aforementioned statistical applications, Juditsky et al.
juditsky2023sparse  proposed an SMD algorithm that exploits state-dependent
noise assumption to attain optimal convergence rates in the situation of
“dominating stochastic error”, i.e., when the amplitude of the error of the
gradient observation is comparable to the amplitude of the gradient $\nabla f$
of the problem objective. In the similar setting, the authors of
woodworth2021even  have recently established sharp lower complexity bounds for
stochastic optimization under state-dependent noise in the Euclidean setting
for both general and strongly convex situations. However, it is well-known
that SMD is suboptimal in the so-call mini-batch setting, where the noise of
the gradient estimator is reduced by using a batch of samples. This setting
has been widely used for applications of stochastic optimization algorithms
especially under a distributed computing environment. To achieve the
accelerated convergence, the application of the classical SAGD algorithm of
lan2012optimal  to the state-dependent noise setting was the subject of
cotter2011better ; liu2018mass ; woodworth2021even . The authors of
woodworth2021even  proved (expected) optimal accuracy bound
$\mathcal{O}\big{(}{\cal H}/k^{2}\big{)}$ after $k$ iterations for the
“standard” SAGD under condition of uniform (for all $\xi$) ${\cal
H}$-Lipschitz continuity on the stochastic operator ${\cal G}(\cdot,\xi)$.
However, this assumption impose significant limitations on the form of the
stochastic operator in statistical learning applications and is violated in
the simple case of unbounded (e.g., Gaussian) regressors $\phi_{t}$, etc. To
conclude, in spite of these efforts, to the best of our knowledge, the
question of building an optimal stochastic approximation routine of general
smooth convex optimization under state-dependent noise assumption (SN) (when
only $\nabla f$ rather than $\cal G(\cdot,\xi)$ is Lipschitz continuous) has
not received a complete answer.

###  1.1 Contributions and organization

Given the state of affairs, this paper focuses on designing accelerated
algorithms and providing sharp analysis for the general stochastic
optimization problem (1.1) with state-dependent noise. Our contribution is
threefold.

  1. 1.

We analyze the convergence rates of the generic (non-Euclidean) SAGD for
solving stochastic optimization with state-dependent noise. We show that under
condition (SN), SAGD attains a convergence rate

$\displaystyle\mathcal{O}\left(\tfrac{LR^{2}}{k^{2}}+\tfrac{{\cal
L}R^{2}}{km}+\tfrac{\sqrt{L{\cal
L}}R^{2}}{k\sqrt{m}}+\sqrt{\tfrac{\sigma_{*}^{2}R^{2}}{km}}\right)$

where $k$ is the number of iterations, $R$ is the initial distance to the
optimal solution, and $m$ is the batch size. The terms in the above bound are
optimal, except for the third term which has a sub-optimal dependence on the
batch size $m$. As a consequence, in order to achieve the optimal iteration
complexity $\mathcal{O}\big{(}\sqrt{LR^{2}/\epsilon}\big{)}$, SAGD requires a
larger batch size, resulting in a sub-optimal sample complexity111In what
follows we refer to the total number $N=N(\epsilon)$ of calls to the
stochastic oracle which are necessary for the approximate solution
$\widehat{x}_{k}$ after $k=k(\epsilon)$ iterations and $N$ oracle calls to
attain the expected (in)accuracy $\epsilon$, i.e.,
${\mathbb{E}}[f(\widehat{x}_{k})]-f^{*}\leq\epsilon$ (1.8) as sample (or
information) $\epsilon$-complexity of the method. We also call iteration
$\epsilon$-complexity the minimal iteration count $k$ such that (1.8) holds.

$\mathcal{O}\left(\sqrt{\tfrac{LR^{2}}{\epsilon}}+{\tfrac{\sqrt{L}{\cal
L}R^{3}}{\epsilon^{3/2}}}+{\tfrac{R^{2}\sigma_{*}^{2}}{\epsilon^{2}}}\right).$

However, imposing the condition of boundness of the second moment of the
Lipschitz constant of the gradient observation ${\cal G}(\cdot,\xi)$ allows to
improve the second term in the above sample complexity bound to
$\mathcal{O}\left(1/\epsilon\right)$. The corresponding iteration complexity
bound $\mathcal{O}\big{(}\sqrt{LR^{2}/\epsilon}\big{)}$ is an improvement
w.r.t. the bound $\mathcal{O}\big{(}\sqrt{{\cal H}R^{2}/\epsilon}\big{)}$ of
woodworth2021even , and our assumption is more general than the “uniform”
Lipschitz continuity assumption in woodworth2021even . Moreover, unlike
woodworth2021even , our analysis does not require the feasible region $X$ to
be a bounded set.

  2. 2.

Under state-dependent noise assumption (SN) we analyze an alternative
accelerated SA algorithm—stochastic gradient extrapolation method (SGE). The
gradient extrapolation method was introduced in lan2018random  by exchanging
the primal and the dual variables in a game interpretation of Nesterov’s
accelerated gradient method (see lan2018random  and Chapter 3 and 4 of
LanBook2020 ). SGE uses the same sequence of points for both gradient
estimations and output solutions. This appears to be a significant advantage
of the SGE over the SAGD in the present setting and allows for direct
“compensation” of the state-dependent noise term by the suboptimality gap of
approximate solutions. As a result, SGE achieves the optimal convergence rate
after $k$ iterations

$\displaystyle\mathcal{O}\left(\tfrac{LR^{2}}{k^{2}}+\tfrac{{\cal
L}R^{2}}{km}+\sqrt{\tfrac{\sigma_{*}^{2}R^{2}}{km}}\right).$

Consequently, it attains the optimal iteration complexity
$\mathcal{O}\big{(}\sqrt{LR^{2}/\epsilon}\big{)}$ along with the optimal
sample complexity $\mathcal{O}\big{(}\sqrt{LR^{2}/\epsilon}+{\cal
L}R^{2}/{\epsilon}+R^{2}\sigma_{*}^{2}/\epsilon^{2}\big{)}$, as supported by
lower bounds in nemyud:83 ; woodworth2021even .

  3. 3.

We propose a multi-stage algorithm with restarts for solving problems
satisfying the quadratic growth condition stating that for some $\mu>0$ and
$x^{*}\in X$,222We suppose for convenience that in this case the optimal
solution $x^{*}$ is unique. Note that (1.9) can be seen as a relaxation of the
strong convexity assumption, i.e., for any $x,y\in X$,
$f(y)-f(x)-\langle\nabla f(x),y-x\rangle\geq\mu\|y-x\|^{2}/2.$

$\displaystyle f(x)-f^{*}\geq\tfrac{\mu}{2}\|x-x^{*}\|^{2},~{}~{}\forall x\in
X.$ (1.9)

We show that the proposed procedure attains the optimal iteration complexity
$\mathcal{O}\big{(}\sqrt{L/\mu}\log(1/{\epsilon})\big{)}$ and the optimal
sample complexity $\mathcal{O}\big{(}\sqrt{L/\mu}\log(1/{\epsilon})+{\cal
L}/\mu\log(1/\epsilon)+\sigma_{*}^{2}/(\mu\epsilon)\big{)}$ simultaneously.
Furthermore, we specify the multi-stage SGE to solve the sparse recovery
problem. This is done by incorporating hard-thresholding of the approximate
solution at the end of each algorithm stage to enforce the sparsity. The
convergence results match the quadratic growth setting up to a multiplicative
factor of the sparsity level $s$, with only logarithmic dependence on the
dimension $n$. To the best of our knowledge, the corresponding convergence
guarantees are new to the optimization literature and lead to extremely
efficient algorithms for sparse recovery problems in the distributed setting.

Remaining sections of the paper are organized as follows. In Section 2, we
formalize the general problem statement and introduce the mini-batch setup. In
Section 3, we study the SAGD method for solving the general convex problem
with state-dependent noise. Section 4 introduces SGE and provides convergence
guarantees in the general convex problem. We introduce the multi-stage SGE for
solving problems satisfying quadratic growth condition in Section 5. Section 6
further extends the multi-stage SGE to be applied to the sparse recovery
problem. Finally, in Section 7, we present some results of a preliminary
simulation study illustrating the numerical performance of the proposed
algorithms in the high-dimensional setting of sparse recovery. Proofs of the
statements are postponed till the Appendix.

###  1.2 Notation

For any $n\geq 1$, we use $[n]$ to denote the set of integers $\\{1,...,n\\}$.
For $x\in\mathbb{R}$, we let $(x)_{+}=\max(x,0)$ and $(x)_{-}=\max(-x,0)$. In
what follows, $E$ is a finite-dimensional real-vector (Euclidean) space. Given
a norm $\|\cdot\|$ on $E$, the associated dual norm $\|\cdot\|_{*}$ is defined
as $\|z\|_{*}:=\sup\\{\langle x,z\rangle:\|x\|\leq 1\\}$. We define
$\omega:E\rightarrow\mathbb{R}$, the distance generating function, which is a
continuously differentiable strongly convex function with modulus $1$. Without
loss of generality, we assume that $\omega(x)\geq\omega(0)=0$ and for some
$\Omega\geq 1$,

$\displaystyle\omega(x)\leq\tfrac{\Omega}{2}\|x\|^{2},\quad\forall x\in E.$
(1.10)

Ideally, we want $\Omega$ to be “not too large”. Meanwhile, a desired distance
generating function should be “prox-friendly”, i.e., for any $a\in E$, the
minimization problem

$\displaystyle\min_{x\in X}\\{\langle a,x\rangle+\omega(x)\\}$

can be easily solved. Note that when $\|\cdot\|$ is the Euclidean norm, the
corresponding $\Omega=1$. Another “standard” choice is
$\|\cdot\|=\|\cdot\|_{1}$, $\|\cdot\|_{*}=\|\cdot\|_{\infty}$, and one can
choose the distance generating function $\omega$ (cf. nesterov2013first )

$\displaystyle\omega(x)=\tfrac{1}{2}e\log n\cdot
n^{(p-1)(2-p)/p}\|x\|_{p}^{2},\quad p=1+\tfrac{1}{\log n};$

the corresponding $\Omega$ satisfies $\Omega\leq e^{2}\log n$ in this case.

Given an initialization $x_{0}\in X$, we define the $x_{0}$-associated
Bregman’s divergence of $x,y\in X$ as

$\displaystyle
V_{x_{0}}(x,y)=\omega(y-x_{0})-\omega(x-x_{0})-\langle\nabla\omega(x-x_{0}),y-x\rangle.$

Clearly, for any $y,x,x_{0}\in X$, we have

$\displaystyle
V_{x_{0}}(x_{0},y)\leq\tfrac{\Omega}{2}\|y-x_{0}\|^{2}\quad\text{and}\quad
V_{x_{0}}(x,y)\geq\tfrac{1}{2}\|x-y\|^{2}.$ (1.11)

We use the shorthand notation $V(x,y)$ for $V_{x_{0}}(x,y)$ when $x_{0}$ is
clear in the content.

Unless stated otherwise, all relations between random variables are assumed to
hold almost surely.

##  2 Problem statement

We summarize below the setting of the stochastic optimization problem and our
assumptions.

###  2.1 Assumptions

The problem under consideration is a stochastic optimization problem (1.1)
with a convex and smooth objective function as given in (1.2).

We assume that stochastic oracle ${\cal G}(\cdot,\cdot)$ is unbiased, i.e.,
satisfies (1.3). Furthermore, we consider the following “state-dependent”
oracle noise assumption.

  * •

[State-dependent variance] We assume that for some ${\cal L}<\infty$,

${\mathbb{E}}_{\xi_{t}}\big{[}\|{\cal
G}(x_{t},\xi_{t})-g(x_{t})\|_{*}^{2}\big{]}\leq{\cal L}[f(x_{t})-f^{*}-\langle
g(x^{*}),x_{t}-x^{*}\rangle]+\sigma_{*}^{2},\quad t\in\mathbb{Z}_{+},$ (SN)

where $x^{*}\in X^{*}$.

Assumption (SN) can be further weakened to $\sigma_{t}^{2}:={\cal
L}[f(x_{t})-f^{*}]+\sigma_{*}^{2}$. In the case of unconditional minimizer
$x^{*}\in\hbox{\rm int\,}X$, the latter condition is clearly equivalent to
(SN). Using the relaxed assumption in the case of $\nabla g(x^{*})\neq 0$
results in convergence guarantees which depend explicitly on $f(x_{0})-f^{*}$.
We utilize (SN) for the sake of convenience, the term $\langle
g(x^{*}),x_{t}-x^{*}\rangle$ in the right-hand side together with smoothness
of $f$ allow us to upper bound the variance of ${\cal G}(x_{0},\xi_{t})$ at
$x_{0}$ by the term proportional to $\|x_{0}-x^{*}\|^{2}$; see, e.g., (4.5).

When proving the convergence rates for the stochastic variant of Nesterov’s
accelerated gradient descent method (SAGD; see Section 3), we also consider
the following assumption:

  * •

[Lipschitz continuous stochastic gradient] For each $\xi\in\Xi$, there exists
a ${\cal K}(\xi)>0$, such that ${\mathbb{E}}_{\xi}[{\cal K}(\xi)^{2}]<\infty$
and

$\|{\cal G}(x,\xi)-{\cal G}(y,\xi)\|_{*}\leq{\cal K}(\xi)\|x-y\|,\quad\forall
x,y\in X.$ (LP)

This assumption relaxes the assumption in woodworth2021even  that assumes
${\cal G}(\cdot,\xi)$ is ${\cal H}$-Lipschitz continuous for all $\xi\in\Xi$.
For instance, for GLR model with Gaussian/sub-Gaussian regressors $\phi_{t}$,
(LP) holds, but the ${\cal H}$-uniform Lipschitz continuous condition is
violated. Nevertheless, (LP) is not a necessary condition of Assumption (SN)
since the latter one may hold in various situations of interest where ${\cal
G}(\cdot,\xi)$ is not Lipschitz (and even not continuous). In Figure 1 we
present the plot of the expectation in the left-hand side of (SN) as a
function of $x$ for two choices of scalar discontinuous gradient observation
${\cal G}_{1}(x,\phi)=\phi\cdot u(\phi x)$ and ${\cal
G}_{2}(x,[\phi,\zeta])=\phi\cdot u(\phi x+\zeta)$ where $\phi$ and $\zeta$ are
independent r.v. with Student $t_{4}$ distribution and
$u(t)=(\tfrac{1}{2}+\sqrt{|t|})\mathrm{sign}(t)$.

|
---|---
Figure 1: Left plot: variance of the stochastic oracle ${\cal G}_{1}(x,\xi)$
as function of $x$ (solid line) and upper bound $4\sigma^{2}(0)+3[f(x)-f(0)]$
(dashed line); right plot: variance of ${\cal G}_{2}(x,\xi)$ as function of
$x$ (solid line) and upper bound $2.3\sigma^{2}(0)+3[f(x)-f(0)]$ (dashed
line).

Given the limitations of Assumption (LP), we will only use it partially in
Section 3 in order to improve the convergence rates of SAGD. In the following
sections, we will propose an alternative accelerated algorithm called SGE that
does not rely on Assumption (LP) but attains stronger convergence guarantees;
see Section 4 for more details.

###  2.2 Mini-batch setup

We consider the mini-batch approach widely used in practice. Specifically, we
assume that at each search point $u_{t}$, the stochastic oracle is called
repeatedly, thus generating $m_{t}$ i.i.d. samples
$\\{\xi_{t,i}\\}_{i=1}^{m_{t}}$, $m_{t}$ being the number of oracle calls.
Next, we compute the unbiased estimator $G_{t}(u_{t})$ of $g(u_{t})$,

$\displaystyle G_{t}(u_{t})=\tfrac{1}{m_{t}}\textstyle\sum_{i=1}^{m_{t}}{\cal
G}(u_{t},\xi_{t,i}).$

We define filtration
$\mathcal{F}_{t}=\sigma\big{(}u_{0},\xi_{0,1},...\xi_{0,m_{0}},...,\xi_{t,1},...\xi_{t,m_{t}}\big{)}$,
so that

  * •

random variables in $\\{\xi_{t,i}\\}_{i=1}^{m_{t}}$ are
$\mathcal{F}_{t}$-measurable.

  * •

search points with index $t$ which are deterministic functions of
$G_{\tau}(u_{\tau}),\tau\leq t-1$, are $\mathcal{F}_{t-1}$-measurable.333Note
that $u_{t}$ here is a general place holder for a
$\mathcal{F}_{t-1}$-measurable search point. With a slight ambiguity of
notation, search points $z_{t}$ and $x_{t}$ in Algorithms 1 and 2 are
$\mathcal{F}_{t}$-measurable.

We use the shorthand notation ${\mathbb{E}}_{\lceil t\rceil}$ to denote the
conditional expectation with respect to the filtration $\mathcal{F}_{t}$.

Based on the state-dependent noise Assumption (SN), we have the following
characterization of the properties of the mini-batch estimator:

###### Lemma 1

The mini-batch estimator $G_{t}$ satisfies

${\mathbb{E}}_{\lceil
t-1\rceil}[\|G_{t}(u_{t})-g(u_{t})\|_{*}^{2}]\leq\tfrac{\Omega}{m_{t}}\cdot{\mathbb{E}}[\|{\cal
G}_{t}(u_{t},\xi_{t,1})-g(u_{t})\|_{*}^{2}],$ (2.1)

where $\Omega$ is defined in (1.10). Consequently, under Assumption (SN), we
have

${\mathbb{E}}_{\lceil
t-1\rceil}[\|G_{t}(u_{t})-g(u_{t})\|_{*}^{2}]\leq\tfrac{\Omega}{m_{t}}\left\\{{\cal
L}[f(u_{t})-f^{*}-\langle
g(x^{*}),u_{t}-x^{*}\rangle]+\sigma_{*}^{2}\right\\}.$ (2.2)

In the bound (2.2), one has $\Omega=1$ in the case of Euclidean setup (when
$\|\cdot\|=\|\cdot\|_{*}=\|\cdot\|_{2})$; in the $\ell_{1}$-setup, one has
$\Omega=\mathcal{O}(\log n)$. One can easily see that, in general, the
logarithmic factor is unavoidable in this case. Indeed, when ${\cal
G}(x,\xi_{i})\sim{\cal R}(n)$ (are $n$-dimensional Rademacher vectors), one
has $\|{\cal G}(x,\xi_{i})\|_{\infty}\leq 1$, while

${\mathbb{E}}\left[\left\|\tfrac{1}{m}\textstyle\sum_{i=1}^{m}{\cal
G}(x,\xi_{i})\right\|_{\infty}^{2}\right]\asymp{\tfrac{\log n}{m}}.$

On the other hand, for ${\cal G}(x,\xi_{i})\sim{\cal N}(0,I_{n})$, we have
that ${\mathbb{E}}[\|{\cal G}(x,\xi_{i})\|_{\infty}^{2}]=2\log n$, and
$\tfrac{1}{m}\sum_{i=1}^{m}{\cal G}(x,\xi_{i})\sim{\cal N}(0,I_{n}/m)$ with

${\mathbb{E}}\left[\left\|\tfrac{1}{m}\textstyle\sum_{i=1}^{m}{\cal
G}(x,\xi_{i})\right\|_{\infty}^{2}\right]={\tfrac{2\log n}{m}},$

and there is no extra logarithmic factor in the bound (2.2) for the norm of
the sum in this situation (in a sense, it is “already paid” in the bound for
the expectation of the norm of a summand).

In the following sections, we present accelerated algorithms equipped with
mini-batches in order to achieve the optimal iteration complexity for
stochastic optimization with state-dependent noise.

##  3 SAGD for general convex problem with state-dependent noise

Algorithm 1 describes a mini-batch variant of the standard stochastic
accelerated gradient descent (SAGD) method. The SAGD method, also called
accelerated stochastic approximation (AC-SA), maintains three sequences of
points. Specifically, $\\{z_{t}\\}$ is the sequence of “prox-centers” for the
prox-mapping updates (3.1b), and $y_{t}$ and $x_{t}$ are weighted averages of
the past $z_{t}$’s; $\\{y_{t}\\}$ is the sequence of search points where
stochastic gradients are estimated using mini-batches, and points
$\\{x_{t}\\}$ is the trajectory of approximate solutions (outputs) at each
iteration.

Algorithm 1 Stochastic accelerated gradient descent method (SAGD)

Input: initial point $z_{0}=x_{0}$, nonnegative nonrandom parameters
$\\{\beta_{t}\\}$ and $\\{\eta_{t}\\}$, and batch size $\\{m_{t}\\}$.

for $t=1,2,\ldots,$ do $\displaystyle y_{t}$
$\displaystyle=(1-\beta_{t})x_{t-1}+\beta_{t}z_{t-1}$ (3.1a) $\displaystyle
G_{t}$ $\displaystyle=\tfrac{1}{m_{t}}\textstyle\sum_{i=1}^{m_{t}}{\cal
G}(y_{t},\xi_{t,i}).$ $\displaystyle z_{t}$ $\displaystyle={\rm argmin}_{z\in
X}\\{\langle G_{t},z\rangle+\eta_{t}V(z_{t-1},z)\\},$ (3.1b) $\displaystyle
x_{t}$ $\displaystyle=(1-\beta_{t})x_{t-1}+\beta_{t}z_{t},$ (3.1c)

end for

We start with the following characterization of the approximate solution
$x_{k}$ of Algorithm 1.

###### Theorem 3.1

Suppose Assumption (SN) is satisfied. Let the algorithmic parameters
$\beta_{t}$ and $\eta_{t}$ satisfy for some $\theta_{t}\geq 0$,

$\displaystyle\theta_{t}\beta_{t}\eta_{t}$
$\displaystyle\leq\theta_{t-1}\beta_{t-1}\eta_{t-1},~{}t=2,...,k$ (3.2a)
$\displaystyle\eta_{t}$ $\displaystyle>L\beta_{t},~{}t=1,...,k.$ (3.2b)

Furthermore, suppose that $\beta_{1}=1$ and

$\displaystyle\theta_{t}(1-\beta_{t})(1+r_{t}{\cal
L})\leq\theta_{t-1},~{}t=2,...,k.$ (3.3)

Then

$\displaystyle\theta_{k}{\mathbb{E}}[f(x_{k})$
$\displaystyle-f^{*}]+\theta_{k}\beta_{k}\eta_{k}{\mathbb{E}}[V(z_{k},x^{*})]$
$\displaystyle\leq\theta_{1}\beta_{1}\eta_{1}V(z_{0},x^{*})+\textstyle\sum_{t=1}^{k}\theta_{t}r_{t}\beta_{t}L{\cal
L}{\mathbb{E}}[V(z_{t-1},x^{*})]+\textstyle\sum_{t=1}^{k}\theta_{t}r_{t}\sigma_{*}^{2}$
(3.4)

where $r_{t}:=\tfrac{\beta_{t}\Omega}{2(\eta_{t}-L\beta_{t})m_{t}}$.

###### Corollary 1

In the premise of Theorem 3.1, suppose that $k\geq 2$, $V(z_{0},x^{*})\leq
D^{2}$, and let $\theta_{t}=(t+1)(t+2)$, $\beta_{t}=\tfrac{3}{t+2}$,
$m_{t}=m$, and $\eta_{t}=\tfrac{\eta}{t+1}$ with

$\eta=\max\left\\{4L,\tfrac{6\Omega(k-1){\cal
L}}{m},\sqrt{\tfrac{9(k+{\color[rgb]{0,0,0}1})^{2}\Omega L{\cal
L}}{m}},\sqrt{\tfrac{2(k+2)^{3}\Omega\sigma_{*}^{2}}{3D^{2}m}}\right\\}.$

Then ${\mathbb{E}}[V(z_{t},x^{*})]\leq 3D^{2}$ for all $1\leq t\leq k$, and

$\displaystyle{\mathbb{E}}[f(x_{k})-f^{*}]\leq\tfrac{12LD^{2}}{(k+1)(k+2)}+\tfrac{6\Omega{\cal
L}D^{2}}{(k+2)m}+\tfrac{18\sqrt{\Omega L{\cal
L}}D^{2}}{(k+1)\sqrt{m}}+\tfrac{4\sqrt{2\Omega\sigma_{*}^{2}D^{2}}}{\sqrt{(k+1)m}}.$
(3.5)

#### Remarks.

Bound (3.5) of Corollary 1 allows us to establish the complexity bounds for
the SAGD algorithm when solving a general convex problem with state-dependent
noise. Let us first consider the case when $m=1$. In this case, we can take
$\Omega/m=1$, and the total number of iterations/oracle calls used by SAGD to
find an $\epsilon$-optimal solution, i.e., $\widehat{x}\in X$ such that
${\mathbb{E}}[f(\widehat{x})-f^{*}]\leq\epsilon$, is bounded by

$\displaystyle\mathcal{O}\left(\sqrt{\tfrac{LD^{2}}{\epsilon}}+\tfrac{{\cal
L}D^{2}}{\epsilon}+\tfrac{\sqrt{L{\cal
L}}D^{2}}{\epsilon}+\tfrac{D^{2}\sigma_{*}^{2}}{\epsilon^{2}}\right).$ (3.6)

In the usual setting with $L=\mathcal{O}({\cal L})$, this upper bound matches
the optimal sample complexity under Assumption (SN), supported by the lower
bound in Theorem 4 of woodworth2021even .

In the mini-batch setting, the bound in (3.5) means that in order to achieve
the optimal iteration complexity of
$\mathcal{O}\left(\sqrt{LD^{2}/\epsilon}\right)$, the batch size
$m\geq\max\left\\{1,\tfrac{k\Omega{\cal L}}{L},\tfrac{k^{2}\Omega{\cal
L}}{L},\tfrac{k^{3}\Omega\sigma_{*}^{2}}{D^{2}L^{2}}\right\\}=\max\left\\{1,\tfrac{k^{2}\Omega{\cal
L}}{L},\tfrac{k^{3}\Omega\sigma_{*}^{2}}{D^{2}L^{2}}\right\\}$ is needed.
Consequently, the total sample complexity is bounded by

$\displaystyle\mathcal{O}\left(\sqrt{\tfrac{LD^{2}}{\epsilon}}+\tfrac{\sqrt{L}{\cal
L}\Omega D^{3}}{\epsilon^{3/2}}+\tfrac{\Omega
D^{2}\sigma_{*}^{2}}{\epsilon^{2}}\right)$ (3.7)

in this situation. When comparing (3.7) to (3.6), we observe that the second
term in (3.7) is sub-optimal. This implies that, based on our analysis, SAGD
does not simultaneously achieve the optimal iteration complexity and sample
complexity under Assumption (SN).

The analysis on the convergence of SAGD in the state-dependent noise setting
reveals the “bottleneck”: in the recursion (3.1), different points $y_{t}$ and
$x_{t}$ are used for gradient estimations and output solutions. Improving the
convergence rates in Corollary 1 requires better control of the objective
value at the points of stochastic gradient estimation. This can be achived by
imposing the Lipschitz regularity assumption in (LP) on stochastic gradients.

###### Lemma 2

Suppose Assumptions (SN) and (LP) hold. Let $\\{x_{t}\\}$, $\\{y_{t}\\}$ and
$\\{z_{t}\\}$ be generated by Algorithm 1. We have

$\displaystyle{\mathbb{E}}_{\lceil
t-1\rceil}\big{[}\|G_{t}-g(y_{t})\|_{*}^{2}\big{]}\leq\tfrac{\Omega}{m_{t}}\left\\{{\cal
K}^{2}\beta_{t}^{2}\|x_{t-1}-z_{t-1}\|^{2}+3{\cal L}[f(x_{t-1})-f^{*}-\langle
g(x^{*}),x_{t-1}-x^{*}\rangle]+3\sigma_{*}^{2}\right\\},$

where

$\displaystyle{\cal K}:=\left(3{\mathbb{E}}_{\xi_{t,1}}[{\cal
K}(\xi_{t,1})^{2}]+3L^{2}\right)^{1/2}.$ (3.8)

We have the following analog of Theorem 3.1 under Assumptions (SN) and (LP).

###### Theorem 3.2

Let Assumptions (SN) and (LP) hold. Let $\beta_{t},\eta_{t}$ satisfy (3.2) for
some $\theta_{t}\geq 0$, and let also $\beta_{1}=1$ and

$\displaystyle\theta_{t}(1-\beta_{t}+3r_{t}{\cal
L})\leq\theta_{t-1},~{}t=2,...,k.$ (3.9)

Then

$\displaystyle\theta_{k}$
$\displaystyle{\mathbb{E}}[f(x_{k})-f^{*}]+\theta_{k}\beta_{k}\eta_{k}{\mathbb{E}}[V(z_{k},x^{*})]$
$\displaystyle\leq\theta_{1}\beta_{1}\eta_{1}V(z_{0},x^{*})+\tfrac{3\theta_{1}r_{1}{\cal
L}L}{2}\|z_{0}-x^{*}\|^{2}+\textstyle\sum_{t=1}^{k}\theta_{t}r_{t}\beta_{t}^{2}{\cal
K}^{2}{\mathbb{E}}\|\bar{z}_{t-1}-z_{t-1}\|^{2}+\textstyle\sum_{t=1}^{k}3\theta_{t}r_{t}\sigma_{*}^{2}$
(3.10)

where ${\cal K}$ is defined in (3.8) and
$r_{t}:=\tfrac{\beta_{t}\Omega}{2(\eta_{t}-L\beta_{t})m_{t}}$.

###### Corollary 2

In the premise of Theorem 3.2, suppose that $k\geq 2$, $V(z_{0},x^{*})\leq
D^{2}$, and let $\beta_{t}=\tfrac{3}{t+2}$, $m_{t}=m$ and
$\eta_{t}=\tfrac{\eta}{t+1}$ with

$\eta=\max\left\\{4L,\tfrac{18\Omega(k+1){\cal
L}}{m},12\sqrt{\tfrac{k\Omega{\cal
K}^{2}}{m}},\sqrt{\tfrac{2(k+2)^{3}\Omega\sigma_{*}^{2}}{D^{2}m}}\right\\}.$

Then ${\mathbb{E}}[V(z_{t},z_{0})]\leq 3D^{2}$, $t=1,...,k$, and

$\displaystyle{\mathbb{E}}[f(x_{k})-f^{*}]\leq\tfrac{13LD^{2}}{(k+1)(k+2)}+\tfrac{54\Omega{\cal
L}D^{2}}{(k+2)m}+\tfrac{72{\cal
K}D^{2}}{(k+2)}\cdot\sqrt{\tfrac{\Omega}{m(k+1)}}+4\sqrt{\tfrac{6\Omega\sigma_{*}^{2}D^{2}}{(k+1)m}}.$
(3.11)

#### Remarks.

Let $m=1$ (and $\Omega/m=1$ in Lemmas 1 and 2). By (3.11), the
iteration/sample complexity of the SAGD is bounded by

$\displaystyle\mathcal{O}\left\\{\sqrt{\tfrac{LD^{2}}{\epsilon}}+\tfrac{{\cal
L}D^{2}}{\epsilon}+\left(\tfrac{{\cal
K}D^{2}}{\epsilon}\right)^{\frac{2}{3}}+\tfrac{D^{2}\sigma_{*}^{2}}{\epsilon^{2}}\right\\}.$

When $\epsilon=\mathcal{O}(L^{3}D^{2}/{\cal K}^{2})$, this sample complexity
is optimal.

In the mini-batch setting, to obtain the optimal iteration complexity
$\mathcal{O}\left(\sqrt{LD^{2}/\epsilon}\right)$, SAGD batch size $m$ should
be in the order of

$\displaystyle\max\left\\{1,\tfrac{k\Omega{\cal L}}{L},\tfrac{k\Omega{\cal
K}^{2}}{L^{2}},\tfrac{k^{3}\Omega\sigma_{*}^{2}}{L^{2}D^{2}}\right\\}.$

As a result, the corresponding sample complexity becomes

$\displaystyle\mathcal{O}\left\\{\sqrt{\tfrac{LD^{2}}{\epsilon}}+\tfrac{\Omega{\cal
L}D^{2}}{\epsilon}+\tfrac{\Omega{\cal K}^{2}D^{2}}{L\epsilon}+\tfrac{\Omega
D^{2}\sigma_{*}^{2}}{\epsilon^{2}}\right\\}.$ (3.12)

When ${\cal K}^{2}=\mathcal{O}(L{\cal L})$, this complexity bound is optimal
(cf. Theorem 4 of woodworth2021even ). The result of Corollary 2 refines the
corresponding statement of woodworth2021even  in three aspects. First, the
corresponding iteration complexity bound
$\mathcal{O}\left(\sqrt{LD^{2}/\epsilon}\right)$ is stated in terms of the
Lipschitz constant of the expected gradient. Second, it relies upon Assumption
(LP) which is significantly weaker than the assumption of uniform Lipschitz
continuity of the stochastic gradient observation ${\cal G}(x,\cdot)$ used in
woodworth2021even . Third, unlike woodworth2021even , our analysis does not
require the feasible region $X$ to be a bounded set.

##  4 SGE for general convex problem with state-dependent noise

In this section, we discuss an alternative acceleration scheme to solve the
general smooth convex problem with state-dependent noise which we refer to as
stochastic gradient extrapolation (SGE). SGE (Algorithm 2) is a variant of the
gradient extrapolation method proposed in lan2018random . We consider here the
mini-batch version of the routine.

Algorithm 2 Stochastic gradient extrapolation method (SGE)

Input: initial point $x_{0}=z_{0}$, nonnegative parameters $\\{\alpha_{t}\\}$,
$\\{\eta_{t}\\}$ and $\\{\beta_{t}\\}$, and batch size $\\{m_{t}\\}$.

$G_{0}=G_{-1}=\tfrac{1}{m_{0}}\textstyle\sum_{i=1}^{m_{0}}{\cal
G}(x_{0},\xi_{0,i})$.

for $t=1,2,\ldots,$ do $\displaystyle\widetilde{G}_{t}$
$\displaystyle=G_{t-1}+\alpha_{t}(G_{t-1}-G_{t-2}),$ (4.1a) $\displaystyle
z_{t}$ $\displaystyle={\rm argmin}_{x\in
X}\\{\langle\widetilde{G}_{t},x\rangle+\eta_{t}V_{x_{0}}(z_{t-1},x)\\},$
(4.1b) $\displaystyle x_{t}$
$\displaystyle={\color[rgb]{0,0,0}(1-\beta_{t})x_{t-1}+\beta_{t}z_{t}},$
(4.1c) $\displaystyle G_{t}$
$\displaystyle=\tfrac{1}{m_{t}}\textstyle\sum_{i=1}^{m_{t}}{\cal
G}(x_{t},\xi_{t,i}).$

end for

As shown in Algorithm 2, the basic algorithmic scheme is conceptually simple.
It involves two sequences of search points $\\{z_{t}\\}$ and $\\{x_{t}\\}$,
the latter being weighted averages of the former. Note that both $x_{t}$ and
$z_{t}$ are $\mathcal{F}_{t-1}$-measurable. Notably, the stochastic gradients
are estimated at the search points $\\{x_{t}\\}$, which are also the
approximate solutions output by the algorithm at each iteration. The key
feature which distinguishes the SGE from the accelerated gradient descent
(SAGD) is the gradient extrapolation step (4.1a), which is essential for
achieving accelerated convergence. The gradient extrapolation approach can be
viewed as the dual counterpart of Nesterov’s accelerated gradient descent.
This connection has been extensively investigated in Section 3 of
lan2018random .

Next, we establish the convergence guarantees of SGE method in expectation,
i.e., ${\mathbb{E}}[f(x_{k})-f^{*}]$. We should stress here that the
convergence analysis of SGE under Assumption (SN) is much more involved than
the ones for its basic scheme in lan2018random , and the SAGD method in
Section 3. Therefore, the details are deferred to the Appendix.

###### Theorem 4.1

Suppose Assumption (SN) holds. Assume that the parameters $\\{\alpha_{t}\\}$,
$\\{\eta_{t}\\}$ and $\\{\beta_{t}\\}$ of Algorithm 2 satisfy

$\displaystyle\theta_{t-1}$
$\displaystyle=\alpha_{t}\theta_{t},~{}~{}\eta_{t}\leq\alpha_{t}\eta_{t-1},~{}~{}t=2,\ldots,k$
(4.2a) $\displaystyle\tfrac{\eta_{t}{(1-\beta_{t})}}{\alpha_{t}}$
$\displaystyle\geq 5L\beta_{t},\quad t=3,\ldots,k$ (4.2b)
$\displaystyle\tfrac{\eta_{1}\eta_{2}}{\alpha_{2}}$ $\displaystyle\geq
25L^{2},~{}~{}\eta_{k}(1-\beta_{k})\geq L\beta_{k},$ (4.2c)

for some $\theta_{t}\geq 0$. Denote

$q_{t}:=\tfrac{\theta_{t+1}(1+\alpha_{t+1}^{2})}{\eta_{t+1}}+\tfrac{\theta_{t+2}\alpha_{t+2}^{2}}{\eta_{t+2}}+\tfrac{\theta_{t+2}^{2}\alpha_{t+2}^{2}}{\theta_{t+1}\eta_{t+1}}\quad\text{and}\quad\epsilon_{t}:=\tfrac{5q_{t}\Omega}{2m_{t}}.$
(4.3)

If $\beta_{1}=1$ and the parameters satisfy

$\theta_{t}\tau_{t}+{\cal
L}\epsilon_{t-1}\leq\theta_{t-1}(1+\tau_{t-1}),~{}t\geq 2,$ (4.4)

then

$\displaystyle\theta_{k}(1+\tau_{k}){\mathbb{E}}[f(x_{k})-f^{*}]\leq\theta_{1}\eta_{1}V(x_{0},x^{*})+\tfrac{{\cal
L}\epsilon_{0}L}{2}\|x_{0}-x^{*}\|^{2}+\sigma_{*}^{2}\cdot\textstyle\sum_{t=0}^{k-1}\epsilon_{t}.$
(4.5)

We now specify a particular stepsize policy in order to establish the
convergence guarantee of SGE.

###### Corollary 3

Let $\theta_{t}=t,~{}\alpha_{t}=\tfrac{t-1}{t},~{}m_{t}=m,$
$\beta_{t}=\tfrac{3}{t+2}$, and $\eta_{t}=\tfrac{\eta}{t}$ for $\eta>0$.
Suppose that $V(x_{0},x^{*})\leq D^{2}$ and that

$\displaystyle\eta=\max\left\\{30L,\tfrac{30\Omega(k+2){\cal
L}}{m},\sqrt{\tfrac{10\Omega(k+1)^{3}\sigma_{*}^{2}}{3mD^{2}}}\right\\}.$

Then

$\displaystyle{\mathbb{E}}[f(x_{k})-f^{*}]\leq\tfrac{91LD^{2}}{k(k+2)}+\tfrac{90\Omega{\cal
L}D^{2}}{mk}+\sqrt{\tfrac{120\Omega\sigma_{*}^{2}D^{2}}{mk}}.$ (4.6)

Furthermore, when $\|x_{0}-x^{*}\|\leq R$ and

$\displaystyle\eta=\max\left\\{30L,\tfrac{30\Omega(k+2){\cal
L}}{m},\sqrt{\tfrac{20(k+1)^{3}\sigma_{*}^{2}}{3mR^{2}}}\right\\},$ (4.7)

we have

$\displaystyle{\mathbb{E}}[f(x_{k})-f^{*}]\leq\tfrac{91L\Omega
R^{2}}{2k(k+2)}+\tfrac{45\Omega^{2}{\cal
L}R^{2}}{mk}+\sqrt{\tfrac{60\Omega^{2}\sigma_{*}^{2}R^{2}}{mk}}.$ (4.8)

#### Remarks.

The bounds of Corollary 3 merit some comments. Observe first that SGE achieves
the optimal sample complexity

$\displaystyle\mathcal{O}\bigg{\\{}\sqrt{\tfrac{LD^{2}}{\epsilon}}+\tfrac{{\cal
L}D^{2}}{\epsilon}+\tfrac{D^{2}\sigma_{*}^{2}}{\epsilon^{2}}\bigg{\\}}$ (4.9)

in the case of $m=1$. In the mini-batch setting, by setting the batch size of
$m\geq\max\left\\{1,\tfrac{\Omega{\cal L}k}{L},\tfrac{\Omega
k^{3}\sigma_{*}^{2}}{D^{2}L^{2}}\right\\}$, the iteration complexity of the
algorithm is bounded by $\mathcal{O}\left(\sqrt{LD^{2}/{\epsilon}}\right)$ and
the overall sample complexity is

$\displaystyle\mathcal{O}\bigg{\\{}\sqrt{\tfrac{LD^{2}}{\epsilon}}+\tfrac{{\cal
L}\Omega D^{2}}{\epsilon}+\tfrac{\Omega
D^{2}\sigma_{*}^{2}}{\epsilon^{2}}\bigg{\\}}.$ (4.10)

Notably, the SGE attains the optimal iteration and sample complexity bounds.

Similarly, Ineq. (4.8) states a convergence rate with condition
$\|x_{0}-x^{*}\|\leq R$. This bound will be be further used in the proof of
the multi-stage SGE method in the next section.

##  5 SGE for convex problem with quadratic growth condition

In this section, we consider the problem setting in which the smooth objective
$f$ in (1.1) satisfies the quadratic growth condition (cf. (1.9)), i.e., when
for some $\mu>0$ and $x_{*}\in X$

$\displaystyle f(x)-f^{*}\geq\tfrac{\mu}{2}\|x-x^{*}\|^{2},~{}~{}\forall x\in
X.$ (5.1)

We propose a multi-stage routine with restarts that utilizes Algorithm 2 as a
working horse.

Algorithm 3 Multi-stage stochastic gradient extrapolation method

Input: initial point ${y}^{0}\in X$. Let $R_{0}$ be a positive real.

for $k=1,2,\ldots,K$ do

(a) Set $R_{k}=R_{0}2^{-k/2}$, $N=\left\lceil 10\sqrt{\tfrac{2\Omega
L}{\mu}}\,\right\rceil$. Run $N$ iterations of SGE (Algorithm 2) with
$x_{0}=z_{0}={y}^{k-1}$ and  $\displaystyle\theta_{t}$
$\displaystyle=t,~{}~{}\alpha_{t}=\tfrac{t-1}{t},~{}~{}{\beta_{t}=\tfrac{3}{t+2}},~{}~{}\eta_{t}=\tfrac{\eta}{t},~{}~{}t=1,...,N$
$\displaystyle\eta$ $\displaystyle=\max\left\\{30L,\tfrac{30\Omega({N}+2){\cal
L}}{{m^{k}}},\sqrt{\tfrac{20({N}+1)^{3}\sigma_{*}^{2}}{{6m^{k}R_{k}^{2}}}}\right\\},$
$\displaystyle m^{k}$
$\displaystyle=\max\left\\{1,\left\lceil\tfrac{18\Omega{\cal
L}({N}+2)}{L}\right\rceil,\left\lceil\tfrac{15{N}({N}+2)^{2}\sigma_{*}^{2}}{{2L^{2}R_{k}^{2}}}\right\rceil\right\\},~{}~{}t=0,...,{N}.$
(5.2)

(b) Set ${y}^{k}=x_{N}$, where $x_{N}$ is the SGE solution obtained in Step
(a).

end for

Algorithm 3 has a simple structure. Each stage $k$ of the routine, consists of
$N_{k}=N$ iterations of the SGE with initial condition $x_{0}={y}^{k-1}$ being
the approximate solution at the end of the stage $k-1$. Method parameters are
selected in such a way that the upper bound $R_{k}^{2}$ for the expected
squared distance ${\mathbb{E}}\|{y}^{k}-x^{*}\|^{2}$ between the approximate
solution ${y}^{k}$ at the end of the $k$-th stage and the optimal solutions
reduces by factor 2.

###### Corollary 4

Let $\\{{y}^{K}\\}$ be approximate solution by Algorithm 3 after $K\geq 1$
stages. Assume that $\|{y}^{0}-x^{*}\|^{2}\leq R_{0}^{2}$. Then

$\displaystyle{\mathbb{E}}[f({y}^{K})-f^{*}]\leq\mu R_{0}^{2}\cdot
2^{-K-1}\quad\text{and}\quad{\mathbb{E}}[\|{y}^{K}-x^{*}\|^{2}]\leq
R_{0}^{2}\cdot 2^{-K}.$

#### Remarks.

By Corollary 4, the number of stages of Algorithm 3 required to attain the
expected inaccuracy $\epsilon$ is bounded with $\mathcal{O}\left(\log(\mu
R_{0}^{2}/{\epsilon})\right)$. When recalling what is the total number $N$ of
iterations at each stage, we conclude that the “total” iteration complexity of
the method is $\mathcal{O}\left(\sqrt{L\Omega/{\mu}}\cdot\log(\mu
R_{0}^{2}/{\epsilon})\right)$. Consequently, the corresponding sample
complexity $\sum_{k=1}^{K}\sum_{t=1}^{N}m^{k}$ is order of

$\displaystyle\mathcal{O}\left\\{\sqrt{\tfrac{L\Omega}{\mu}}\log(\tfrac{\mu
R_{0}^{2}}{\epsilon})+\tfrac{{\cal L}\Omega^{2}}{\mu}\log(\tfrac{\mu
R_{0}^{2}}{\epsilon})+\tfrac{\Omega^{2}\sigma_{*}^{2}}{\mu\epsilon}\right\\}.$

Similarly, the iteration complexity of solution ${y}^{K}$ satisfying
$\mathbb{E}[\|{y}^{K}-x^{*}\|^{2}]\leq\epsilon^{2}$ does not exceed

$\mathcal{O}\left(\sqrt{L\Omega/{\mu}}\cdot\log(R_{0}/{\epsilon})\right);$

the corresponding sample complexity is then bounded by

$\displaystyle\mathcal{O}\bigg{\\{}\sqrt{\tfrac{L\Omega}{\mu}}\log(\tfrac{R_{0}}{\epsilon})+\tfrac{{\cal
L}\Omega^{2}}{\mu}\log(\tfrac{R_{0}}{\epsilon})+\tfrac{\Omega^{2}\sigma_{*}^{2}}{\mu^{2}\epsilon^{2}}\bigg{\\}}.$

Similar to Corollary 3, the three terms in the sample complexity bounds
represent the deterministic error, the state-dependent stochastic error, and
the state-independent stochastic error, respectively. Regardless of the
dependence on $\Omega$, the multi-stage SGE achieves the optimal iteration
complexity and sample complexity simultaneously, supported by the lower bound
in Theorem 5 of woodworth2021even .

##  6 SGE for sparse recovery

An interesting application of stochastic optimization with state-dependent
noise arises in the relation to the problem of sparse recovery when it is
assumed that (1.1) has a sparse or low-rank solution $x^{*}$. This problem is
motivated by applications in high-dimensional statistics, where stochastic
estimation under sparsity or low-rank constraints has garnered significant
attention. To solve this problem, one usually builds a sample average
approximation (SAA) of the expected risk and solves the resulting minimization
problem. To enhance sparsity, a classical approach consists of incorporating a
regularization term, e.g., $\ell_{1}$\- or trace-norm penalty (as in Lasso),
and minimizing the norm of the solution under constraint (as in Dantzig
Selector), see, e.g., bickel2009simultaneous ; candes2007dantzig ;
candes2011probabilistic ; candes2011tight ; candes2012exact ; candes2006stable
; candes2009near ; fazel2008compressed ; juditsky2011accuracy ;
tsybakov2011nuclear  among many others.

Stochastic approximation also serves as a standard approach to deal with the
sparse recovery problem. However, utilizing traditional Euclidean stochastic
approximation usually leads to sub-optimal complexity bounds: in this setting,
the expected squared $\ell_{2}$-error of the stochastic operator ${\cal
G}(x,\xi)$ is usually proportional to the problem dimension $n$. Therefore,
non-Euclidean stochastic mirror descent (SMD) methods have been applied to
address this issue. In particular, the SMD algorithm in shalev2009stochastic ;
srebro2010smoothness  attains the high probability complexity bound
$f(\widehat{x})-f^{*}=\mathcal{O}(\sigma\sqrt{s/N})$ (up to some “logarithmic
factors”) under the sub-Gaussian noise assumption, referred to as “slow rates”
for sparse recovery. To improve the convergence rate to
$\mathcal{O}(\sigma^{2}s/N)$, multi-stage routines exploiting the properties
similar to strong/uniform convexity could be used, cf. ghadimi2012optimal ;
juditsky2011first ; juditsky2014deterministic . In agarwal2012stochastic ;
gaillard2017sparse , the authors utilize the “restricted” strong convexity
condition to establish $\mathcal{O}(\tfrac{\sigma^{2}s}{N})$ complexity bounds
when assuming that $N\gg s^{2}$. The latter assumption means that the optimal
rates are only valid in the range $s\ll\sqrt{N}$ of sparsity parameter.
Recently, in juditsky2023sparse ; ilandarideva2023stochastic , new multi-stage
stochastic mirror descent algorithms were proposed which rely on the idea of
variance reduction. The proposed method improved the required number of
iterations in each stage to $\mathcal{O}(s)$, thus attaining the best-known
state-dependent stochastic error. However, the iteration complexity of those
routines is still sub-optimal in the mini-batch setting, leaving room for
further acceleration.

In this section, we suppose that stochastic optimization problem (1.1) admits
a sparse solution $x^{*}\in X$. We consider three types of sparsity structures
on $E$:

  * •

“Vanilla” sparsity: we assume that an optimal solution $x^{*}\in X$ has at
most $s\ll n$ nonvanishing entries. We put $\|\cdot\|=\|\cdot\|_{1}$ and
$\|\cdot\|_{*}=\|\cdot\|_{\infty}$.

  * •

Group sparsity: let us partition the set $[n]$ into $b$ subsets
$\\{I_{1},...,I_{B}\\}$, and let $x_{b}$ the $b$th block of $x$, meaning that
$[x_{b}]_{i}=0$ for all $i\notin I_{b}$. We assume that the optimal $x^{*}\in
X$ is a block vector with at most $s\leq B$ nonvanishing blocks $x_{b}$. We
define $\|x\|=\textstyle\sum_{b=1}^{B}\|x_{b}\|_{2}$ (block
$\ell_{1}/\ell_{2}$-norm) and $\|x\|_{*}=\max_{b\leq B}\|x_{b}\|_{2}$ (block
$\ell_{\infty}/\ell_{2}$-norm).

  * •

Low rank sparsity: consider the matrix space $\mathbb{R}^{p\times q}$ where
$p\geq q$ eqipped with Frobenius inner product. We assume that the optimal
$x^{*}\in X$ satisfies $\text{rank}(x^{*})\leq s$. We consider the nuclear
norm $\|x\|=\textstyle\sum_{i=1}^{q}\sigma_{i}(x)$ where $\sigma_{i}(\cdot)$
are the singular values of $x$, so that
$\|x\|_{*}=\max_{i\in[q]}\sigma_{i}(x)$ is the spectral norm.

In what follows we assume that problem (1.1) and norms $\|\cdot\|$ and
$\|\cdot\|_{*}$ satisfy conditions (1.2), (1.3), and state-dependent noise
conditions in Assumptions (SN) and (LP) of Section 2. Also, instead of
$\mu$-quadratic growth (with $\mu>0$) condition with respect to $\|\cdot\|$ as
in the previous section, we suppose $f$ and $x^{*}$ verify the quadratic
growth condition with respect to the Euclidean norm, i.e.,

$\displaystyle
f(x)-f^{*}\geq\tfrac{1}{2}\underline{\kappa}\|x-x^{*}\|_{2}^{2},~{}~{}\forall
x\in X,$ (6.1)

for some $\underline{\kappa}>0$.

###  6.1 Application: sparse generalized linear regression

Let us check that problem assumptions of this section hold in the case of
generalized linear regression problem (GLR), as described in the introduction.
Let us assume that

  * •

regressors $\phi_{t}$ satisfy
$\mathbb{E}[\phi_{1}\phi_{1}^{T}]=\Sigma\succeq\kappa I$ with $\kappa>0$ and
$\|\Sigma\|_{\infty}\leq\nu$;

  * •

noises $\zeta_{t}$ are zero mean with bounded variance, i.e.,
$\mathbb{E}[\zeta_{1}^{2}]\leq 1$ without loss of generality;

  * •

activation function $u$ is strongly monotone and Lipschitz continuous, i.e.,

$\displaystyle\underline{r}(t-t^{\prime})\leq
u(t)-u(t^{\prime})\leq\bar{r}(t-t^{\prime}).$ (6.2)

As already explained in the introduction, estimation of $x^{*}\in\hbox{\rm
int\,}X$ may be addressed through solving the stochastic optimization problem

$\displaystyle\min_{x\in
X}\left\\{f(x):={\mathbb{E}}[v(\phi^{T}x)-\phi^{T}x\eta]\right\\}$ (6.3)

where $v^{\prime}(t)=u(t)$. The gradient of the problem objective and its
stochastic estimate are given by

$\displaystyle
g(x)={\mathbb{E}}[\phi\big{(}u(\phi^{T}x)-\eta\big{)}]\quad\text{and}\quad{\cal
G}(x,(\phi,\zeta)):=\phi\big{(}u(\phi^{T}x)-\eta\big{)}=\phi\big{(}u(\phi^{T}x)-u(\phi^{T}x^{*})\big{)}-\phi\zeta.$

It is easy to see that condition (1.3) is verified in this case, and invoking
${\mathbb{E}}[\eta]={\mathbb{E}}[u(\phi^{T}x)]$, we conclude that
$g(x^{*})=0$. To check the quadratic growth condition (6.1), we write

$\displaystyle f(x)-f^{*}$
$\displaystyle=\int_{0}^{1}g(x^{*}+t(x-x^{*}))^{T}(x-x^{*})dt$
$\displaystyle=\int_{0}^{1}{\mathbb{E}}\\{\phi[u(\phi^{T}x^{*}+t(x-x^{*}))-u(\phi^{T}x^{*})]\\}^{T}(x-x^{*})dt$
[by (6.2)]
$\displaystyle\geq\int_{0}^{1}\underline{r}{\mathbb{E}}\\{[\phi^{T}(x-x^{*})]^{2}\\}tdt=\tfrac{\underline{r}}{2}\|x-x^{*}\|_{\Sigma}^{2}\geq\tfrac{\underline{r}\kappa}{2}\|x-x^{*}\|_{2}^{2}.$
(6.4)

Therefore, condition (6.1) holds with
$\underline{\kappa}=\underline{r}\cdot\kappa$ which is independent of problem
dimension $n$.

Since we are interested in the high-dimensional setting, the desired recovery
error should have, at most, logarithmic dependence in the problem dimension
$n$. However, the $\ell_{2}$ variance of the stochastic first-order
information ${\mathbb{E}}\|\mathcal{G}-g\|_{2}^{2}$ is proportional to the
problem dimension $n$, making the standard Euclidean SA methods not
applicable. To address this issue, we work in the non-Euclidean setting with
$\|\cdot\|=\|\cdot\|_{1}$ and $\|\cdot\|_{*}=\|\cdot\|_{\infty}$. Next, let us
examine the smoothness condition (1.2) and state-dependent variance condition
(SN). Note that for all $x,x^{\prime}\in\mathbb{R}^{n}$,

$\displaystyle\|g(x)-g(x^{\prime})\|_{\infty}$
$\displaystyle=\sup_{\|z\|_{1}\leq 1}\langle
g(x)-g(x^{\prime}),z\rangle=\sup_{\|z\|_{1}\leq
1}{\mathbb{E}}\\{\phi^{T}z[u(\phi^{T}x)-u(\phi^{T}x^{\prime})]\\}$
$\displaystyle\overset{(i)}{\leq}\sup_{\|z\|_{1}\leq
1}\bar{r}{\mathbb{E}}\\{|\phi^{T}z||\phi^{T}(x-x^{\prime})|\\}\overset{(ii)}{\leq}\bar{r}\sup_{\|z\|_{1}\leq
1}\sqrt{{\mathbb{E}}\\{(\phi^{T}z)^{2}\\}}\cdot\|x-x^{\prime}\|_{\Sigma}$
$\displaystyle\leq\bar{r}\sqrt{\nu}\|x-x^{\prime}\|_{\Sigma},$ (6.5)

where (i) is a consequence of (6.2) and (ii) follows from the Cauchy
inequality. Consequently, we have

$\displaystyle{\mathbb{E}}\|{\cal G}(x,(\phi,\eta))-g(x)\|_{\infty}^{2}$
$\displaystyle={\mathbb{E}}\|\phi\big{(}u(\phi^{T}x)-u(\phi^{T}x^{*})\big{)}+\phi\zeta+[g(x^{*})-g(x)]\|_{\infty}^{2}$
[by (6.1)] $\displaystyle\leq
3{\mathbb{E}}\|\phi\big{(}u(\phi^{T}x)-u(\phi^{T}x^{*})\big{)}\|_{\infty}^{2}+3{\mathbb{E}}\\{\|\phi\|_{\infty}^{2}\\}\sigma^{2}+3\bar{r}^{2}\nu\|x-x^{*}\|_{\Sigma}^{2}$
[by (6.2)] $\displaystyle\leq
3\bar{r}^{2}{\mathbb{E}}\\{\|\phi\|_{\infty}^{2}\big{(}\phi^{T}x-\phi^{T}x^{*}\big{)}^{2}\\}+3{\mathbb{E}}\\{\|\phi\|_{\infty}^{2}\\}\sigma^{2}+3\bar{r}^{2}\nu\|x-x^{*}\|_{\Sigma}^{2}.$

By (6.1), we conclude that the condition (SN) holds whenever

$\displaystyle{\mathbb{E}}[\|\phi\|_{\infty}^{2}\big{(}\phi^{T}x-\phi^{T}x^{*}\big{)}^{2}]\lesssim{\mathbb{E}}[\big{(}\phi^{T}(x-x^{*})\big{)}^{2}]\lesssim\|x-x^{*}\|_{\Sigma}^{2};$

e.g., when the regressor $\phi$ is bounded or sub-Gaussian. Finally, under
similar assumptions on the regressors and sub-Gaussian assumption on the
additive noise $\zeta$, condition (LP) naturally follows.

###  6.2 SGE-SR: stochastic gradient extrapolation for sparse recovery

We extend SGE to solve the sparse recovery problem. Similarly to Algorithm 3,
sparse recovery routine is organized in stages; each stage represents a run of
SGE (Algorithm 2). The principal difference with Algorithm 3, apart from the
different choice of algorithm parameters, is the sparsity enforcing step (see,
e.g., blumensath2009iterative ; jain2014iterative ; liu2020between )
implemented at the end of each stage.

Observe that for $x\in X$ one can efficiently compute a sparse approximation
of $x$, specifically, $x_{s}=\mathsf{sparse}(x)$, an optimal solution to

$\min\|x-z\|_{2}~{}\text{over $s$-sparse $z\in X$}.$

When set $X$ is positive monotone (respectively, block-monotone)444Let for
$x\in\mathbb{R}^{n}$ $|x|$ be a vector in $\mathbb{R}^{n}_{+}$ with entries
which are absolute values of corresponding entries of $x$. We say that $X$ is
positive monotone if whenever $x\in X$ and $|y|\leq|x|$ (the inequality is
understood coordinate-wise), one also has $y\in X$. A typical example of a
monotone convex set $X$ is a ball of an absolute norm in $\mathbb{R}^{n}$.
$x_{s}$ is obtained by “sparsification”—by zeroing all but $s$ largest in
amplitude entries (respectively, $s$ largest in $\ell_{2}$-norm blocks) of
$x$. Moreover, for any $s$-sparse $z$ we have

$\displaystyle\|z\|\leq\sqrt{s}\|z\|_{2}.$ (6.6)

Note also that if $x,y\in E$, $x_{s}:=\mathsf{sparse}(x)$ and $y$ is
$s$-sparse, then

$\|x_{s}-y\|\leq\sqrt{2s}\|x_{s}-y\|_{2}\leq 2\sqrt{2s}\|x-y\|_{2}.$ (6.7)

Indeed, given that both $x_{s}$ and $y$ are $s$-sparse, we have that $x_{s}-y$
is $2s$-sparse, thus $\|x_{s}-y\|\leq\sqrt{2s}\|x_{s}-y\|_{2}.$ On the other
hand, by optimality of $x_{s}$,

$\|x_{s}-y\|_{2}\leq\|x_{s}-x\|_{2}+\|x-y\|_{2}\leq 2\|x-y\|_{2}.$

Algorithm 4 Stochastic gradient extrapolation method for sparse recovery (SGE-
SR)

Input: initial point $\overline{y}^{0}\in X$.

for $k=1,2,\ldots,K$ do

(a) Set $N=40\sqrt{sL\Omega/{\underline{\kappa}}}$ and $R_{k}=2^{-k/2}R_{0}$.
Run $N$ iterations of SGE (Algorithm 2) with $x_{0}=z_{0}=\overline{y}^{k-1}$
and  $\displaystyle\theta_{t}$
$\displaystyle=t,~{}~{}\alpha_{t}=\tfrac{t-1}{t},~{}~{}{\beta_{t}=\tfrac{3}{t+2}},~{}~{}\eta_{t}=\tfrac{\eta}{t},~{}~{}t=1,...,N$
$\displaystyle\eta$ $\displaystyle=\max\left\\{30L,\tfrac{30\Omega(N+2){\cal
L}}{m^{k}},\sqrt{\tfrac{20(N+1)^{3}\sigma_{*}^{2}}{6m^{k}R_{k}^{2}}}\right\\},$
$\displaystyle m^{k}$
$\displaystyle=\max\left\\{1,\left\lceil\tfrac{18\Omega{\cal
L}(N+2)}{L}\right\rceil,\left\lceil\tfrac{15N(N+2)^{2}\sigma_{*}^{2}}{2L^{2}R_{k}^{2}}\right\rceil\right\\},~{}~{}t=0,...,N$

(b) Set ${y}^{k}=x_{N}$, where $x_{N}$ is the solution obtained in Step (a).
Calculate  $\overline{y}^{k}=\mathsf{sparse}({y}^{k}).$

end for

The following corollary characterizes the convergence rate of SGE-SR for
solving the sparse recovery problem.

###### Corollary 5

Let $\\{{y}^{k},\overline{y}^{k}\\}$ be computed by Algorithm 4. Assume
$\|\overline{y}^{0}-x^{*}\|^{2}\leq R_{0}^{2}$. Then we have for $k\geq 1$

$\displaystyle{\mathbb{E}}[f({y}^{k})-f^{*}]\leq\underline{\kappa}s^{-1}R_{0}^{2}\cdot
2^{-k+4}\quad\text{and}\quad{\mathbb{E}}[\|\overline{y}^{k}-x^{*}\|^{2}]\leq
R_{0}^{2}\cdot 2^{-k}.$

#### Remarks.

From the result of Corollary 5 we conclude that the SGE-SR algorithm finds an
$s$-sparse $\overline{y}^{k}\in X$ such that
${\mathbb{E}}[\|y-x^{*}\|^{2}]\leq\epsilon^{2}$ for any $\epsilon\in(0,R_{0})$
in at most $k=\mathcal{O}\left(\log(R_{0}/\epsilon)\right)$ stages. The
corresponding iteration complexity of the SGE-SR is
$\mathcal{O}\left(\sqrt{\tfrac{sL\Omega}{\underline{\kappa}}}\log(\tfrac{R_{0}}{\epsilon})\right)$,
and the overall sample complexity is

$\displaystyle\mathcal{O}\bigg{\\{}\sqrt{\tfrac{sL\Omega}{\underline{\kappa}}}\log(\tfrac{R_{0}}{\epsilon})+\tfrac{s{\cal
L}\Omega^{2}}{\underline{\kappa}}\log(\tfrac{R_{0}}{\epsilon})+\tfrac{\Omega^{2}s^{2}\sigma_{*}^{2}}{\underline{\kappa}^{2}\epsilon^{2}}\bigg{\\}}.$

Similarly, the iteration complexity of the solution ${y}^{k}\in X$ (which is
not $s$-sparse in general) such that
${\mathbb{E}}[f({y}^{k})-f^{*}]\leq\epsilon$ is
$\mathcal{O}\left(\sqrt{\tfrac{sL\Omega}{\underline{\kappa}}}\log(\tfrac{\underline{\kappa}R_{0}^{2}}{s\epsilon})\right)$,
while the total sample complexity is

$\displaystyle\mathcal{O}\bigg{\\{}\sqrt{\tfrac{sL\Omega}{\underline{\kappa}}}\log(\tfrac{\underline{\kappa}R_{0}^{2}}{s\epsilon})+\tfrac{s{\cal
L}\Omega^{2}}{\underline{\kappa}}\log(\tfrac{\underline{\kappa}R_{0}^{2}}{s\epsilon})+\tfrac{\Omega^{2}s\sigma_{*}^{2}}{\underline{\kappa}\epsilon}\bigg{\\}}.$

The above results may be compared to the convergence guarantees obtained in
juditsky2023sparse  in the similar setting of the sparse recovery problem. The
iteration complexity of the SGE-SR algorithm attains the optimal dependence on
the problem’s condition number, $\sqrt{sL/\underline{\kappa}}$, which improves
over the corresponding result in juditsky2023sparse  by a factor of
$\mathcal{O}(\sqrt{sL/\underline{\kappa}})$. Moreover, except for the
dependence on $\Omega$, the proposed solution matches the best-known sample
complexity bounds for the stochastic error; the extra factor $\Omega$ being
due to the mini-batch use (cf. Lemma 1 and subsequent remark) and could
theoretically hinder the method precision. Note that in the problems of
interest, $\Omega$ is logarithmic in problem dimension. It should be mentioned
that in our numerical experiments we did not observe any accuracy degradation
when using the mini-batch algorithm.

##  7 Numerical experiments

In this section we present a simulation study illustrating numerical
performance of the proposed routines. We consider the sparse recovery problem
in generalized linear regression (GLR) model with random design as discussed
in the previous section. Recall that we are looking to recover the $s$-sparse
vector $x^{*}\in\mathbb{R}^{n}$ from i.i.d observations

$\eta_{i}=u(\phi_{i}^{T}x^{*})+\sigma\zeta_{i},~{}i=1,2,\dots,N.$

In experiments we report on below, the activation function $u(\cdot)$ is of
the form

$u_{\alpha}(x)=x\mathbf{1}\\{|x|\leq\alpha\\}+\mathrm{sign}(x)[\alpha^{-1}(|x|^{\alpha}-1)+1]\mathbf{1}\\{|x|>\alpha\\},\quad\alpha>0,\;x\in\mathbb{R}.$

We consider three different activations, namely, the linear link function
$u_{1}(\cdot)$, $u_{1/2}(\cdot)$, and $u_{1/10}(\cdot)$ (cf. Figure 2).

Figure 2: Activation functions

In our simulations, $s$ nonvanishing components of the signal $x^{*}$ are
sampled from the $s$-dimensional standard Gaussian distribution. We explore
two setups—light tailed and heavy tailed—for generating regressors and
additive noises. In the light-tail setup, regressors $\phi_{i}$ are
independently drawn from a multivariate Gaussian distribution
$\phi_{i}\sim\mathcal{N}(0,\Sigma)$, where $\Sigma$ is a diagonal covariance
matrix with diagonal entries $0<\Sigma_{1,1}\leq\dots\leq\Sigma_{n,n}$. In the
heavy-tail setup, regressors are independently drawn from a multivariate
Student distribution $\phi_{i}\sim t_{n}(\nu,0,\Sigma)$, $\nu$ being the
corresponding degree of freedom kotz2004multivariate . The condition number
$\kappa$ of the problem is defined as the ratio of the largest and the
smallest eigenvalues of $\Sigma$. The additive noise of the model in the
light-tail setup is the zero-mean Gaussian noise with variance $\sigma^{2}$;
in the heavy-tail setup, the additive noises have (scaled) univariate Student
distribution $\eta_{i}\sim\lambda t(\nu)$, $\nu\geq 3$, with scale parameter
$\lambda=\sqrt{(\nu-2)/\nu}$ with unit variance. Because of the memory
limitations, observations $(\eta_{i},\phi_{i})$ are generated on the fly at
each oracle call.

In all our experiments, we run 50 simulation trials (with randomly generated
regressors and noises); then we trace in the plots the median and the first
and the last deciles of the error $\|x_{t}-x^{*}\|_{2}$.

The aim of the first series of experiments is to compare the procedure
described in Section 6 to the SMD-SR algorithm of juditsky2023sparse ;555SMD-
SR is a stochastic approximation algorithm for sparse recovery utilizing hard
thresholding which relies upon “vanilla” non-Euclidean mirror descent; both
algorithms use the same distance generating function
$\omega(x)=c(n)\|x\|^{2}_{p}$. in the light-tail noise setting, the
corresponding results are presented in Figure 3. In Figure 4, we present
results of simulations of the accelerated algorithm in the light-tail and
heavy-tail setup. We used the same algorithmic parameters in both simulation
setups. In the above experiments, $n=500\,000$, the maximal number of calls to
the stochastic oracle (estimation sample size) $N=250\,000$, and sparsity
level $s=250$; unless stated otherwise, the problem condition number is set to
$\kappa=1$.

Figure 3: Estimation error $\|x_{t}-x^{*}\|_{2}$ against the number of
stochastic oracle calls for SGE-SR and SMD-SR algorithms. In the left, middle,
and right columns of the plot we show results for the linear activation
$u_{1}$, and nonlinear $u_{1/2}$ and $u_{1/10}$, respectively. Two figure rows
correspond to two different noise levels, $\sigma=0.1$ (the upper row) and
$\sigma=0.001$ (the bottom row). The legend specifies the value $m_{0}$ of the
batch size of the preliminary phase of the algorithm for both routines.

Figure 4: Estimation error $\|x_{t}-x^{*}\|_{2}$ against the number of
stochastic oracle calls for SGE-SR in Gaussian (light-tail) and Student
$t_{5}$ (heavy-tail) regressor and noise generation setups. In the left,
middle, and right columns of the plot we show results for the linear
activation $u_{1}$, and nonlinear $u_{1/2}$ and $u_{1/10}$, respectively. Two
figure rows correspond to two different noise levels, $\sigma=0.1$ (the upper
row) and $\sigma=0.001$ (the bottom row). The legend specifies the value
$m_{0}$ of the batch size of the preliminary phase of the algorithm for both
routines.

In the second series of experiments we put $n=100\,000$, $N=200\,000$, and
$s=50$. Experiments reported in Figure 5 illustrate the impact of the
condition number on the convergence of the SGE-SR algorithm.

Figure 5: Error $\|x_{t}-x^{*}\|_{2}$ of the SGE-SR algorithm for three values
of the problem condition number. First row: algorithm error against the number
of oracle calls; second row: error against the number of algorithm iterations.
Figure columns correspond to the results for $u_{1}$, $u_{1/2}$, and
$u_{1/10}$ activation functions and $\sigma=0.001$.

Finally, we illustrate the performance of SGE-SR and SMD-SR algorithms which
share the same size of mini-batch. Recall that both algorithms if “normally
set”—SGE-SR with mini-batch of optimal size and SMD-SR with “trivial” mini-
batch (of size $m_{0}=1$)—converge linearly during the preliminary phase. In
Figure 6 we report on the simulation of algorithms with the same (optimal for
the accelerated algorithm) size of the mini-batch juditsky2023sparse .

The series of experiments conducted indicates that the SGE-SR algorithm
outperforms its non-accelerated counterpart. Despite both algorithms
exhibiting linear rates of convergence, SGE-SR reaches a better precision for
a fixed number of samples in every setting. This advantage is also observed
when the two algorithms are compared in terms of the number of iterations,
where the accelerated algorithm clearly outperforms its non-accelerated
counterpart by a significant margin. The improved iteration complexities make
SGE-SR a viable solution for distributed sparse recovery problems, where it is
crucial to reduce the number of communication rounds between the servers and
the clients while keeping high precision.

Figure 6: SGE-SR compared to “vanilla” SMD-SR and its mini-batch variant.
First row: error $\|x_{t}-x^{*}\|_{2}$ against the number of oracle calls;
second row: the error against the number of algorithm iterations. Figure
columns correspond to the results for $u_{1}$, $u_{1/2}$, and $u_{1/10}$
activation functions and $\sigma=0.001$.

##  8 Concluding remarks

In this paper, we investigate the problem of stochastic smooth convex
optimization with “state-dependent” variance of stochastic gradients. We study
two non-Euclidean accelerated stochastic approximation algorithms, stochastic
accelerated gradient descent (SAGD) and stochastic gradient extrapolation
(SGE), and provide optimal iteration and sample complexities for both
algorithms under appropriate conditions. However, the optimal convergence
guarantees for SGE require less restrictive assumptions, thus leading to wider
applications such as statistical estimation problems with heavy tail noises.
In addition, we propose a multi-stage routine of SGE to solve problems that
satisfy the quadratic growth condition and further extend it to the sparse
recovery problem. Our theoretical guarantees are corroborated by numerical
experiments in high-dimensional settings. Further research will be directed to
proving large deviation bounds to ensure the reliability and robustness of the
solutions.

##  Appendix A Appendix: proofs

###  A.1 Proof of Lemmas 1 and 2

#### Proof of Lemma 1.

For notational simplicity, put

$\displaystyle\delta_{t}:=G_{t}(u_{t})-g(u_{t}),\quad\delta_{t,i}:={\cal
G}_{t}(u_{t},\xi_{t,i})-g(u_{t}).$

Let $\omega^{*}$ denote the convex conjugate of the distance generating
function $\omega$. We write

$\displaystyle\|\delta_{t}\|_{*}^{2}$
$\displaystyle=\Omega^{2}\|\tfrac{\delta_{t}}{\Omega}\|_{*}^{2}=2\Omega^{2}\max_{z}[\langle\tfrac{\delta_{t}}{\Omega},z\rangle-\tfrac{1}{2}\|z\|^{2}]\underset{\text{by
(\ref{V_to_norm_general})}}{\leq}2\Omega\max_{z}[\langle\delta_{t},z\rangle-\omega(z)]=2\Omega\,\omega^{*}\big{(}\delta_{t}\big{)}.$
(A.1)

By the strong convexity of $\omega$, $\omega^{*}$ is smooth with $1$-Lipschitz
continuous gradient. Thus

$\displaystyle\omega^{*}\big{(}\delta_{t}\big{)}\leq\omega^{*}\big{(}\tfrac{1}{m_{t}}\textstyle\sum_{i=1}^{m_{t}-1}\delta_{t,i}\big{)}+\left\langle\nabla\big{(}\tfrac{1}{m_{t}}\textstyle\sum_{i=1}^{m_{t}-1}\delta_{t,i}\big{)},\tfrac{\delta_{t,m_{t}}}{m_{t}}\right\rangle+\tfrac{1}{2}\Big{\|}\tfrac{\delta_{t,m_{t}}}{m_{t}}\Big{\|}_{*}^{2}.$

Now, recursively using the above relationship and the independence of
$\delta_{t,i},~{}i\in[m_{t}]$, we get

${\mathbb{E}}_{\lceil
t-1\rceil}[\omega^{*}\big{(}\delta_{t}\big{)}]\leq\tfrac{1}{2m_{t}}\cdot{\mathbb{E}}\|\delta_{t,1}\|_{*}^{2},$

and we obtain the bound in (2.2) by combing the above inequality with (A.1).
$\square$

#### Proof of Lemma 2.

Note that $x_{t-1}$ and $y_{t}$ are $\mathcal{F}_{t-1}$-measurable. By Lemma
1,

$\displaystyle{\mathbb{E}}_{\lceil
t-1\rceil}\big{[}\|G_{t}-g(y_{t})\|_{*}^{2}\big{]}\leq\tfrac{\Omega}{m_{t}}{\mathbb{E}}_{\lceil
t-1\rceil}\big{[}\|{\cal G}(y_{t},\xi_{t,1})-g(y_{t})\|_{*}^{2}\big{]}$
$\displaystyle\leq\tfrac{3\Omega}{m_{t}}\Big{\\{}{\mathbb{E}}_{\lceil
t-1\rceil}\big{[}\|{\cal G}(y_{t},\xi_{t,1})-{\cal
G}(x_{t-1},\xi_{t,1})\|_{*}^{2}\big{]}+{\mathbb{E}}_{\lceil
t-1\rceil}\big{[}\|{\cal G}(x_{t-1},\xi_{t,1})-g(x_{t-1})\|_{*}^{2}\big{]}$
(A.2) $\displaystyle\qquad+{\mathbb{E}}_{\lceil
t-1\rceil}\big{[}\|g(x_{t-1})-g(y_{t})\|_{*}^{2}\big{]}\Big{\\}}$
$\displaystyle\leq\tfrac{3\Omega}{m_{t}}\left\\{\left({\mathbb{E}}_{\xi_{t,1}}[{\cal
K}(\xi_{t,1})^{2}]+L^{2}\right)\cdot\|y_{t}-x_{t-1}\|^{2}+\big{[}{\cal
L}[f(x_{t-1})-f^{*}-\langle
g(x^{*}),x_{t-1}-x^{*}\rangle]+\sigma_{*}^{2}\big{]}\right\\},$ (A.3)

the last inequality being the consequence of Assumptions (LP) and (SN). Now,
utilizing
$x_{t-1}-y_{t}=x_{t-1}-(1-\beta_{t})x_{t-1}-\beta_{t}z_{t-1}=\beta_{t}(x_{t-1}-z_{t-1})$,
we obtain the desired result. $\square$

###  A.2 Proof of Theorem 3.1

For sack of completeness, we start with the following statement similar to the
previous results for ASGD (e.g., Proposition 3.1 in LanBook2020 ).

###### Proposition 1

Let $\\{z_{t}\\},\\{y_{t}\\}$ and $\\{x_{t}\\}$ be generated by Algorithm 1.
Suppose that $\\{\beta_{t}\\}$ and $\\{\eta_{t}\\}$ satisfy (3.2) for some
$\theta_{t}\geq 0$. Then for any $x\in X$ one has

$\displaystyle\textstyle\sum_{t=1}^{k}\theta_{t}[f(x_{t})-f(x)]+\theta_{k}\beta_{k}\eta_{k}V(z_{k},x)$
$\displaystyle\leq\textstyle\sum_{t=1}^{k}\theta_{t}(1-\beta_{t})[f(x_{t-1})-f(x)]+\theta_{1}\beta_{1}\eta_{1}V(z_{0},x)+\textstyle\sum_{t=1}^{k}\theta_{t}\beta_{t}\langle\delta_{t},x-z_{t-1}\rangle+\textstyle\sum_{t=1}^{k}\tfrac{\theta_{t}\beta_{t}\|\delta_{t}\|_{*}^{2}}{2(\eta_{t}-L\beta_{t})}$
(A.4)

where $\delta_{t}:=G_{t}-g(y_{t})$.

#### Proof of the proposition.

First, by convexity of $f$ and due to the definition of $x_{t}$, we have for
all $z\in X$,

$\displaystyle f(z)+\langle g(z),x_{t}-z\rangle$
$\displaystyle=(1-\beta_{t})[f(z)+\langle
g(z),x_{t-1}-z\rangle]+\beta_{t}[f(z)+\langle g(z),z_{t}-z\rangle]$
$\displaystyle\leq(1-\beta_{t})f(x_{t-1})+\beta_{t}[f(z)+\langle
g(z),z_{t}-z\rangle].$

Now, by the smoothness of $f$,

$\displaystyle f(x_{t})$ $\displaystyle\leq f(z)+\langle
g(z),x_{t}-z\rangle+\tfrac{L}{2}\|x_{t}-z\|^{2}$
$\displaystyle\leq(1-\beta_{t})f(x_{t-1})+\beta_{t}[f(z)+\langle
g(z),z_{t}-z\rangle]+\tfrac{L}{2}\|x_{t}-z\|^{2},$

so that for $z=y_{t}$ we have

$\displaystyle f(x_{t})\leq(1-\beta_{t})f(x_{t-1})+\beta_{t}[f(y_{t})+\langle
g(y_{t}),z_{t}-y_{t}\rangle]+\tfrac{L}{2}\|x_{t}-y_{t}\|^{2}.$ (A.5)

On the other hand, the optimality condition for (3.1b) in Algorithm 1 implies
the following relationship (see, e.g., Lemma 3.5 of LanBook2020 ):

$\displaystyle\langle
G_{t},z_{t}-y_{t}\rangle+\eta_{t}V(z_{t-1},z_{t})+\eta_{t}V(z_{t},x)\leq\langle
G_{t},x-y_{t}\rangle+\eta_{t}V(z_{t-1},x),\quad\forall x\in X.$ (A.6)

By combining (A.5) and (A.6) we obtain for all $x\in X$,

$\displaystyle
f(x_{t})+\beta_{t}\eta_{t}V(z_{t-1},z_{t})+\beta_{t}\eta_{t}V(z_{t},x)$
$\displaystyle\qquad\leq(1-\beta_{t})f(x_{t-1})+\beta_{t}[f(y_{t})+\langle
g(y_{t}),x-y_{t}\rangle]$
$\displaystyle\qquad\quad+\beta_{t}\langle\delta_{t},x-z_{t}\rangle+\beta_{t}\eta_{t}V(z_{t-1},x)+\tfrac{L}{2}\|x_{t}-y_{t}\|^{2}$
$\displaystyle\qquad\leq(1-\beta_{t})f(x_{t-1})+\beta_{t}[f(y_{t})+\langle
g(y_{t}),x-y_{t}\rangle]$
$\displaystyle\qquad\quad+\beta_{t}\langle\delta_{t},z_{t-1}-z_{t}\rangle+\beta_{t}\langle\delta_{t},x-z_{t-1}\rangle+\beta_{t}\eta_{t}V(z_{t-1},x)+\tfrac{L}{2}\|x_{t}-y_{t}\|^{2}.$
(A.7)

From (3.1a) and (3.1c) we have $x_{t}-y_{t}=\beta_{t}(z_{t}-z_{t-1})$; after
substituting into (A.2) and taking into account (1.11), we obtain

$\displaystyle f(x_{t})+\beta_{t}\eta_{t}V(z_{t},x)$
$\displaystyle\quad\leq(1-\beta_{t})f(x_{t-1})+\beta_{t}[f(y_{t})+\langle
g(y_{t}),x-y_{t}\rangle]$
$\displaystyle\quad\quad+\beta_{t}\langle\delta_{t},x-z_{t-1}\rangle+\beta_{t}\eta_{t}V(z_{t-1},x)+\beta_{t}\|\delta_{t}\|_{*}\|z_{t-1}-z_{t}\|+(\tfrac{L\beta_{t}^{2}}{2}-\tfrac{\eta_{t}\beta_{t}}{2})\|z_{t}-z_{t-1}\|^{2}$
$\displaystyle\quad\leq(1-\beta_{t})f(x_{t-1})+\beta_{t}[f(y_{t})+\langle
g(y_{t}),x-y_{t}\rangle]+\beta_{t}\langle\delta_{t},x-z_{t-1}\rangle+\beta_{t}\eta_{t}V(z_{t-1},x)+\tfrac{\beta_{t}\|\delta_{t}\|_{*}^{2}}{2(\eta_{t}-L\beta_{t})},$

the last inequality being a consequence of (3.1c) and Young’s inequality. When
subtracting $f(x)$ on both sides we get

$\displaystyle[f(x_{t})-f(x)]+\beta_{t}\eta_{t}V(z_{t},x)$
$\displaystyle\leq(1-\beta_{t})[f(x_{t-1})-f(x)]+\beta_{t}{\color[rgb]{0,0,0}}[\underbrace{f(y_{t})+\langle
g(y_{t}),x-y_{t}\rangle-f(x)}_{\leq
0}]+\beta_{t}\eta_{t}V(z_{t-1},x)+\beta_{t}\langle\delta_{t},x-z_{t-1}\rangle+\tfrac{\beta_{t}\|\delta_{t}\|_{*}^{2}}{2(\eta_{t}-L\beta_{t})}$
$\displaystyle\leq(1-\beta_{t})[f(x_{t-1})-f(x)]+\beta_{t}\eta_{t}V(z_{t-1},x)+\beta_{t}\langle\delta_{t},x-z_{t-1}\rangle+\tfrac{\beta_{t}\|\delta_{t}\|_{*}^{2}}{2(\eta_{t}-L\beta_{t})},$

and after multiplying by $\theta_{t}$ and summing up from $t=1$ to $k$ we
arrive at

$\displaystyle\textstyle\sum_{t=1}^{k}$
$\displaystyle\theta_{t}[f(x_{t})-f(x)]+\theta_{k}\beta_{k}\eta_{k}V(z_{k},x)$
$\displaystyle\leq\textstyle\sum_{t=1}^{k}\theta_{t}(1-\beta_{t})[f(x_{t-1})-f(x)]+\theta_{1}\beta_{1}\eta_{1}V(z_{0},x)+\textstyle\sum_{t=1}^{k}\theta_{t}\beta_{t}\langle\delta_{t},x-z_{t-1}\rangle+\textstyle\sum_{t=1}^{k}\tfrac{\theta_{t}\beta_{t}\|\delta_{t}\|_{*}^{2}}{2(\eta_{t}-L\beta_{t})},$

which is (A.4). $\square$

#### Proof of the theorem.

When setting $x=x^{*}$, using (2.2) in Lemma 1, and taking expectation on both
sides of (A.4), we get

$\displaystyle\textstyle\sum_{t=1}^{k}\theta_{t}{\mathbb{E}}[f(x_{t})-f^{*}]+\theta_{k}\beta_{k}\eta_{k}{\mathbb{E}}[V(z_{k},x^{*})]$
$\displaystyle\quad\leq\textstyle\sum_{t=1}^{k}\theta_{t}(1-\beta_{t}){\mathbb{E}}[f(x_{t-1})-f^{*}]+\theta_{1}\beta_{1}\eta_{1}V(z_{0},x^{*})+\textstyle\sum_{t=1}^{k}\theta_{t}r_{t}{\mathbb{E}}\big{[}{\cal
L}\big{(}f(y_{t})-f^{*}-\langle
g(x^{*}),y_{t}-x^{*}\rangle\big{)}+\sigma_{*}^{2}\big{]}.$

Recall that $y_{t}=(1-\beta_{t})x_{t-1}+\beta_{t}z_{t-1}$. On the other hand,
by smoothness of $f$, when recalling that $\langle
g(x^{*}),x_{t-1}-x^{*}\rangle\geq 0$,

$\displaystyle f(y_{t})-f^{*}-\langle g(x^{*}),y_{t}-x^{*}\rangle$
$\displaystyle\quad\leq(1-\beta_{t})\big{[}f(x_{t-1})-f^{*}-\langle
g(x^{*}),x_{t-1}-x^{*}\rangle\big{]}+\beta_{t}\big{[}f(z_{t-1})-f^{*}-\langle
g(x^{*}),z_{t-1}-x^{*}\rangle\big{]}$
$\displaystyle\quad\leq(1-\beta_{t})\big{[}f(x_{t-1})-f^{*}-\langle
g(x^{*}),x_{t-1}-x^{*}\rangle\big{]}+\tfrac{L\beta_{t}}{2}\|z_{t-1}-x^{*}\|^{2}$
$\displaystyle\quad\leq(1-\beta_{t})\big{[}f(x_{t-1})-f^{*}\big{]}+L\beta_{t}V(z_{t-1},x^{*}).$

When combining the above inequalities we obtain

$\displaystyle\textstyle\sum_{t=1}^{k}\theta_{t}{\mathbb{E}}[f(x_{t})-f^{*}]+\theta_{k}\beta_{k}\eta_{k}{\mathbb{E}}[V(z_{k},x^{*})]$
$\displaystyle\leq\textstyle\sum_{t=1}^{k}\theta_{t}(1-\beta_{t})(1+r_{t}{\cal
L}){\mathbb{E}}[f(y_{t-1})-f^{*}]+\theta_{1}\beta_{1}\eta_{1}V(z_{0},x^{*})+\textstyle\sum_{t=1}^{k}\left\\{\theta_{t}r_{t}\beta_{t}L{\cal
L}{\mathbb{E}}[V(z_{t-1},x^{*})]+\theta_{t}r_{t}\sigma_{*}^{2}\right\\}.$

Due to (3.3), and taking into account that $\beta_{1}=1$, we conclude that

$\displaystyle\theta_{k}{\mathbb{E}}[f(x_{k})$
$\displaystyle-f^{*}]+\theta_{k}\beta_{k}\eta_{k}{\mathbb{E}}[V(z_{k},x^{*})]\leq\theta_{1}\beta_{1}\eta_{1}V(z_{0},x^{*})+\textstyle\sum_{t=1}^{k}\theta_{t}r_{t}\beta_{t}L{\cal
L}{\mathbb{E}}[V(z_{t-1},x^{*})]+\textstyle\sum_{t=1}^{k}\theta_{t}r_{t}\sigma_{*}^{2},$

what is (3.4). $\square$

###  A.3 Proof of Corollary 1

#### 1o.

Observe first that for $\eta\geq 4L$ and $m_{t}=m$ one has

$\displaystyle
r_{t}=\frac{\beta_{t}\Omega}{2[\eta/(t+1)-3L/(t+2)]m}\leq\frac{2\beta_{t}\Omega(t+1)}{\eta
m}=\frac{2\beta_{t}\Omega}{\eta_{t}m}.$ (A.8)

As a result, when $\eta\geq\tfrac{6\Omega(k-1){\cal L}}{m}$ and given the
definition of $\beta_{t}$ and $\theta_{t}$, we have

$\displaystyle\theta_{t}(1-\beta_{t})(1+r_{t}{\cal
L})\leq(t+1)(t-1)\left(1+\tfrac{6\Omega(t+1){\cal L}}{(t+2)\eta
m}\right)\leq(t+1)(t-1)\left(1+\tfrac{(t+1)}{(t+2)(k-1)}\right)\leq
t(t+1)=\theta_{t-1},$

so that (3.3) holds. On the other hand,

$\displaystyle\theta_{t}\beta_{t}\eta_{t}=(t+1)(t+2)\cdot\tfrac{3\eta}{(t+2)(t+1)}=3\eta=\theta_{t-1}\beta_{t-1}\eta_{t-1},$
(A.9)

so condition (3.2a) is satisfied, so Theorem 3.1 applies. By (A.8), we also
have

$\displaystyle\theta_{t}r_{t}\beta_{t}L{\cal
L}\leq\tfrac{2\theta_{t}\beta_{t}^{2}\Omega L{\cal
L}}{\eta_{t}m}\leq\tfrac{18\Omega L{\cal L}}{\eta_{t}m}.$ (A.10)

#### 2o.

Next, let us check that ${\mathbb{E}}[V(z_{t},x^{*})]\leq 3D^{2}$ for all
$t\geq 1$. Indeed, for $t=1$ we have by (3.4), the definition of
$\theta_{t},\beta_{t},\eta_{t}$ and due to (A.8)–(A.10)

$\displaystyle{\mathbb{E}}[V(z_{1},x^{*})]$ $\displaystyle\leq
V(z_{0},x^{*})+\tfrac{\theta_{1}\beta_{1}r_{1}L{\cal
L}}{\theta_{1}\beta_{1}\eta_{1}}V(z_{0},x^{*})+\tfrac{\theta_{1}r_{1}}{\theta_{1}\beta_{1}\eta_{1}}\sigma_{*}^{2}\leq
V(z_{0},x^{*})+\tfrac{{\color[rgb]{0,0,0}8}\Omega L{\cal
L}}{m\eta^{2}}V(z_{0},x^{*})+\tfrac{8\Omega}{\eta^{2}m}\sigma_{*}^{2}$
$\displaystyle\overset{(ii)}{\leq}V(z_{0},x^{*})+V(z_{0},x^{*})+D^{2}\leq
3D^{2}.$

where the last inequality follows from
$\eta^{2}\geq\max\left\\{{\tfrac{9(k+{\color[rgb]{0,0,0}1})^{2}\Omega L{\cal
L}}{m}},{\tfrac{2(k+2)^{3}\Omega\sigma_{*}^{2}}{3D^{2}m}}\right\\}$.

Now we assume that for $s=1,...,t-1$, ${\mathbb{E}}[V(z_{s},x^{*})]\leq
3D^{2}$. Then

$\displaystyle{\mathbb{E}}[V(z_{t},x^{*})]$
$\displaystyle\leq\tfrac{\theta_{1}\beta_{1}\eta_{1}}{\theta_{t}\beta_{t}\eta_{t}}V(z_{0},x^{*})+\textstyle\sum_{s=1}^{t}\tfrac{\theta_{s}r_{s}\beta_{s}L{\cal
L}}{\theta_{t}\beta_{t}\eta_{t}}{\mathbb{E}}[V(z_{s-1},x^{*})]+\textstyle\sum_{s=1}^{t}\tfrac{\theta_{s}r_{s}}{\theta_{t}\beta_{t}\eta_{t}}\sigma_{*}^{2}$
$\displaystyle\leq
V(z_{0},x^{*})+\textstyle\sum_{s=1}^{t}\tfrac{r_{s}}{\eta_{s}}L{\cal
L}{\mathbb{E}}[V(z_{s-1},x^{*})]+\textstyle\sum_{s=1}^{t}\tfrac{r_{s}}{\beta_{t}\beta_{s}\eta_{s}}\sigma_{*}^{2}$
$\displaystyle\leq V(z_{0},x^{*})+\tfrac{3(t+{\color[rgb]{0,0,0}1})^{2}\Omega
L{\cal L}}{\eta^{2}m}\cdot
3D^{2}+\tfrac{2(t+2)^{3}\Omega}{3\eta^{2}m}\cdot\sigma_{*}^{2}$
$\displaystyle\leq V(z_{0},x^{*})+D^{2}+D^{2}\leq 3D^{2}$

where the last inequality is, again, a consequence of
$\eta^{2}\geq\max\left\\{{\tfrac{9(k+{\color[rgb]{0,0,0}1})^{2}\Omega L{\cal
L}}{m}},{\tfrac{2(k+2)^{3}\Omega\sigma_{*}^{2}}{3D^{2}m}}\right\\}$.

#### 3o.

Finally, when substituting the above estimates into (3.4), we obtain

$\displaystyle{\mathbb{E}}[f(x_{k})-f^{*}]$
$\displaystyle\leq\tfrac{1}{\theta_{k}}\cdot\left(\theta_{1}\beta_{1}\eta_{1}V(z_{0},x^{*})+\textstyle\sum_{t=1}^{k}3\theta_{t}r_{t}\beta_{t}L{\cal
L}D^{2}+\textstyle\sum_{t=1}^{k}\theta_{t}r_{t}\sigma_{*}^{2}\right)$
$\displaystyle\overset{(i)}{\leq}\tfrac{1}{(k+1)(k+2)}\cdot\left(3\eta
D^{2}+\tfrac{27(k+{\color[rgb]{0,0,0}1})^{2}\Omega L{\cal L}D^{2}}{\eta
m}+\tfrac{2(k+2)^{3}\Omega\sigma_{*}^{2}}{\eta m}\right)$
$\displaystyle\overset{(ii)}{\leq}\tfrac{12LD^{2}}{(k+1)(k+2)}+\tfrac{18\Omega{\cal
L}D^{2}}{(k+2)m}+\tfrac{18\sqrt{\Omega L{\cal
L}}D^{2}}{(k+1)\sqrt{m}}+\tfrac{2\sqrt{6(k+2)\Omega\sigma_{*}^{2}D^{2}}}{(k+1)\sqrt{m}}$

where (i) follows from (A.8)–(A.10), and (ii) is a consequence of the
definition of $\eta$. This implies (3.5) due to
$\tfrac{k+2}{k+1}\leq\tfrac{4}{3}$ for $k\geq 2$. $\square$

###  A.4 Proof of Theorem 3.2 and Corollary 2

The subsequent proofs follow those of Theorem 3.1 and Corollary 1. We present
them here for reader’s convenience.

#### Proof of Theorem 3.2.

Note that assumptions of Proposition 1 hold. When taking expectation on both
sides of (A.4) and using the bound Lemma 2, we obtain for $x=x_{*}$,

$\displaystyle\textstyle\sum_{t=1}^{k}\theta_{t}{\mathbb{E}}[f(x_{t})-f^{*}]+\theta_{k}\beta_{k}\eta_{k}{\mathbb{E}}[V(z_{k},x^{*})]$
$\displaystyle\quad\leq\textstyle\sum_{t=1}^{k}\theta_{t}(1-\beta_{t}){\mathbb{E}}[f(x_{t-1})-f^{*}]+\theta_{1}\beta_{1}\eta_{1}V(z_{0},x^{*})$
$\displaystyle\quad+\textstyle\sum_{t=1}^{k}\theta_{t}r_{t}{\mathbb{E}}\big{[}{\cal
K}^{2}\beta_{t}^{2}\|x_{t-1}-z_{t-1}\|^{2}+3{\cal
L}\big{(}f(x_{t-1})-f^{*}-\langle
g(x^{*}),x_{t-1}-x^{*}\rangle\big{)}+3\sigma_{*}^{2}\big{]}.$

By rearranging the terms and utilizing (3.9) and $\beta_{1}=1$, we get

$\displaystyle\theta_{k}{\mathbb{E}}[f(x_{k})-f^{*}]+\theta_{k}\beta_{k}\eta_{k}{\mathbb{E}}[V(z_{k},x^{*})]$
$\displaystyle\leq\theta_{1}\beta_{1}\eta_{1}V(z_{0},x^{*})+3\theta_{1}r_{1}{\cal
L}\big{(}f(z_{0})-f^{*}-\langle
g(x^{*}),z_{0}-x^{*}\rangle\big{)}+\textstyle\sum_{t=1}^{k}\theta_{t}r_{t}{\mathbb{E}}\big{[}{\cal
K}^{2}\beta_{t}^{2}{\mathbb{E}}\|x_{t-1}-z_{t-1}\|^{2}+3\sigma_{*}^{2}\big{]}$
$\displaystyle\leq\theta_{1}\beta_{1}\eta_{1}V(z_{0},x^{*})+\tfrac{3\theta_{1}r_{1}{\cal
L}L}{2}\|z_{0}-x^{*}\|^{2}+\textstyle\sum_{t=1}^{k}\theta_{t}r_{t}{\mathbb{E}}\big{[}{\cal
K}^{2}\beta_{t}^{2}{\mathbb{E}}\|x_{t-1}-z_{t-1}\|^{2}+3\sigma_{*}^{2}\big{]},$

what is (3.10). $\square$

#### Proof of Corollary 2.

Observe that for $\eta\geq 4L$ and $m_{t}=m$ (cf. (A.8))

$\displaystyle r_{t}\leq\tfrac{2\beta_{t}\Omega(t+1)}{\eta m}.$ (A.11)

Then, due to $\eta\geq\tfrac{18\Omega(k+1){\cal L}}{m}$,

$\displaystyle\theta_{t}(1-\beta_{t}+3r_{t}{\cal
L})\leq(t+1)(t+2)\left(\tfrac{t-1}{t+2}+\tfrac{18\Omega(t+1){\cal
L}}{(t+2)\eta m}\right)\leq(t+1)\left(t-1+\tfrac{t+1}{k+1}\right)\leq
t(t+1)=\theta_{t-1},$

thus (3.9) holds. We have (cf. (A.9))

$\displaystyle\theta_{t}\beta_{t}\eta_{t}=(t+1)(t+2)\cdot\tfrac{3\eta}{(t+2)(t+1)}=3\eta=\theta_{t-1}\beta_{t-1}\eta_{t-1},$
(A.12)

so (3.2a) is verified and the conclusion of Theorem 3.2 applies. We also have

$\displaystyle\theta_{t}r_{t}\beta_{t}^{2}\leq\tfrac{2\theta_{t}\beta_{t}^{3}\Omega(t+1)}{\eta
m}\leq\tfrac{54\Omega}{\eta m}.$ (A.13)

Let us now check that ${\mathbb{E}}[V(z_{t},x^{*})]\leq 3D^{2}$ for all $t\geq
1$. First, for $t=1$, applying (3.10) and taking into account (A.11)–(A.13),
we get

$\displaystyle{\mathbb{E}}[V(z_{1},x^{*})]$ $\displaystyle\leq
V(z_{0},x^{*})+\tfrac{3\theta_{1}r_{1}{\cal
L}L}{2\theta_{1}\beta_{1}\eta_{1}}\|z_{0}-x^{*}\|^{2}+\tfrac{3\theta_{1}r_{1}}{\theta_{1}\beta_{1}\eta_{1}}\sigma_{*}^{2}$
$\displaystyle\leq V(z_{0},x^{*})+\tfrac{12\Omega{\cal
L}L}{m\eta^{2}}\|z_{0}-x^{*}\|^{2}+\tfrac{24\Omega}{\eta^{2}m}\sigma_{*}^{2}\leq
V(z_{0},x^{*})+\tfrac{1}{24}V(z_{0},x^{*})+D^{2}\leq 3D^{2},$

the second inequality being due to (1.11) and
$\eta\geq\max\left\\{4L,\tfrac{18\Omega(k+1){\cal
L}}{m},\sqrt{\tfrac{2(k+2)^{3}\Omega\sigma_{*}^{2}}{D^{2}m}}\right\\}$.

Now, assume that for $s=1,...,t-1$, ${\mathbb{E}}[V(z_{s-1},x^{*})]\leq
3D^{2}$. Because $x_{s}$ is a weighted average of the previous iterates, by
convexity of $\|\cdot\|^{2}$, we have for $s<t$,

$\displaystyle{\mathbb{E}}\|x_{s}-z_{s}\|^{2}$ $\displaystyle\leq
2{\mathbb{E}}\|x_{s}-x^{*}\|^{2}+2{\mathbb{E}}\|z_{s}-x^{*}\|^{2}$
$\displaystyle\leq 2\max_{0\leq i\leq
s}{\mathbb{E}}\|z_{i}-x^{*}\|^{2}+2{\mathbb{E}}\|z_{s}-x^{*}\|^{2}$
$\displaystyle\leq 4\max_{0\leq i\leq
s}{\mathbb{E}}[V(z_{i},x^{*})]+4{\mathbb{E}}[V(z_{s},x^{*})]\leq 24D^{2}.$
(A.14)

As a consequence, when substituting into (3.10) the bounds of (A.11)–(A.13)

$\displaystyle{\mathbb{E}}[V(z_{k},x^{*})]$ $\displaystyle\leq
V(z_{0},x^{*})+\tfrac{3\theta_{1}r_{1}{\cal
L}L}{2\theta_{k}\beta_{k}\eta_{k}}\|z_{0}-x^{*}\|^{2}+\textstyle\sum_{t=2}^{k}\tfrac{24\theta_{t}r_{t}\beta_{t}^{2}{\cal
K}^{2}D^{2}}{\theta_{k}\beta_{k}\eta_{k}}+\textstyle\sum_{t=1}^{k}\tfrac{3\theta_{t}r_{t}\sigma_{*}^{2}}{\theta_{k}\beta_{k}\eta_{k}}$
$\displaystyle\overset{(ii)}{\leq}V(z_{0},x^{*})+\tfrac{12\Omega{\cal
L}L}{m\eta^{2}}\|z_{0}-x^{*}\|^{2}+\tfrac{144\Omega(k-1){\cal
K}^{2}D^{2}}{\eta^{2}m}+\tfrac{2(k+2)^{3}\Omega\sigma_{*}^{2}}{\eta^{2}m}\overset{(i)}{\leq}3D^{2},$

due to

$\eta\geq\max\left\\{4L,\tfrac{18\Omega(k+1){\cal
L}}{m},12\sqrt{\tfrac{k\Omega{\cal
K}^{2}}{m}},\sqrt{\tfrac{2(k+2)^{3}\Omega\sigma_{*}^{2}}{D^{2}m}}\right\\}.$

Finally, from (3.10) and the definition of $\eta$ we conclude that

$\displaystyle{\mathbb{E}}[f(x_{k})-f^{*}]$
$\displaystyle\leq\tfrac{1}{\theta_{k}}\cdot\left(\theta_{1}\beta_{1}\eta_{1}V(z_{0},x^{*})+\tfrac{3\theta_{1}r_{1}{\cal
L}L}{2}\|z_{0}-x^{*}\|^{2}+\textstyle\sum_{t=2}^{k}\theta_{t}r_{t}\beta_{t}^{2}{\cal
K}^{2}{\mathbb{E}}\|\bar{z}_{t-1}-z_{t-1}\|^{2}+\textstyle\sum_{t=1}^{k}3\theta_{t}r_{t}\sigma_{*}^{2}\right)$
$\displaystyle\leq\tfrac{1}{(k+1)(k+2)}\cdot\left(3\eta
D^{2}+\tfrac{36\Omega{\cal L}L}{m\eta}+\tfrac{432(k-1)\Omega{\cal
K}^{2}D^{2}}{\eta m}+\tfrac{6(k+2)^{3}\Omega\sigma_{*}^{2}}{\eta m}\right)$
$\displaystyle\leq\tfrac{13LD^{2}}{(k+1)(k+2)}+\tfrac{54\Omega{\cal
L}D^{2}}{(k+2)m}+\tfrac{72{\cal
K}D^{2}}{(k+2)}\cdot\sqrt{\tfrac{\Omega}{m(k+1)}}+\tfrac{6}{k+1}\cdot\sqrt{\tfrac{2(k+2)\Omega\sigma_{*}^{2}D^{2}}{m}}.$

This completes the proof due to $\tfrac{k+2}{k+1}\leq\tfrac{4}{3}$ for $k\geq
2$. $\square$

###  A.5 Proof of Theorem 4.1

For $0\leq\beta_{t}\leq 1$, denote $\tau_{t}=\tfrac{1-\beta_{t}}{\beta_{t}}$.
Relationships (4.2) in variables $\theta_{t},\alpha_{t},\eta_{t}$, and
$\tau_{t}$ become

$\displaystyle\theta_{t-1}$
$\displaystyle=\alpha_{t}\theta_{t},~{}~{}\eta_{t}\leq\alpha_{t}\eta_{t-1},~{}~{}t=2,\ldots,k$
(A.15a) $\displaystyle\tfrac{\eta_{t}\tau_{t-1}}{\alpha_{t}}$
$\displaystyle\geq 5L,\quad t=3,\ldots,k$ (A.15b)
$\displaystyle\tfrac{\eta_{1}\eta_{2}}{\alpha_{2}}$ $\displaystyle\geq
25L^{2},~{}~{}\eta_{k}\tau_{k}\geq L,$ (A.15c)

with (4.1c) of Algorithm 2 replaced with

$x_{t}=(z_{t}+\tau_{t}x_{t-1})/(1+\tau_{t}).$ (A.16)

We need the following technical statement.

###### Proposition 2

If the algorithm parameters $\\{\alpha_{t}\\}$, $\\{\eta_{t}\\}$ and
$\\{\tau_{t}\\}$ satisfy $\tau_{1}\geq 0,$ $\tau_{t}>0$ for all $t\geq 2$,
along with relations (A.15a)–(A.15c) for some $\theta_{t}\geq 0$. Then for all
$x\in X$

$\displaystyle\textstyle\sum_{t=1}^{k}\theta_{t}\left\\{\tau_{t}[f(x_{t})-f(x_{t-1})]-\langle
g(x_{t}),x-x_{t}\rangle\right\\}\leq\theta_{1}\eta_{1}V(x_{0},x)+\textstyle\sum_{t=0}^{k-1}\left[p_{t}\langle\delta_{t},x-z_{t}\rangle+\tfrac{5q_{t}}{2}\|\delta_{t}\|_{*}^{2}\right]$
(A.17)

where $~{}p_{t}:=\begin{cases}~{}\theta_{t},&t\leq k-2,\\\
~{}\theta_{t+1}+\theta_{t},&t=k-1.\end{cases}$

#### Proof of the proposition.

#### 1o.

Denote $g_{t}:=g(x_{t})$. By the smoothness of $f$,

$\tfrac{1}{2L}\|g_{t}-g_{t-1}\|_{*}^{2}\leq f(x_{t-1})-[f(x_{t})+\langle
g_{t},x_{t-1}-x_{t}\rangle].$

Then for any $x\in X$, using the above inequality and the fact that
$x-x_{t}=x-z_{t}-\tau_{t}(x_{t-1}-x_{t})$ due to (A.16), we get

$\displaystyle\tau_{t}f(x_{t})-\langle g_{t},x-x_{t}\rangle$
$\displaystyle=\tau_{t}[f(x_{t})+\langle g_{t},x_{t-1}-x_{t}\rangle]+\langle
g_{t},z_{t}-x\rangle$
$\displaystyle\leq\tau_{t}[f(x_{t-1})-\tfrac{1}{2L}\|g_{t}-g_{t-1}\|_{*}^{2}]+\langle
g_{t},z_{t}-x\rangle.$

By the optimality condition of (4.1b), we have

$\displaystyle\langle\widetilde{G}_{t},z_{t}-x\rangle\leq\eta_{t}V(z_{t-1},x)-\eta_{t}V(z_{t},x)-\eta_{t}V(z_{t-1},z_{t}).$

Combining two previous inequalities, we obtain

$\displaystyle\tau_{t}f(x_{t})-\langle
g_{t},x-x_{t}\rangle-\tau_{t}f(x_{t-1})$
$\displaystyle\leq\eta_{t}V(z_{t-1},x)-\eta_{t}V(z_{t},x)+\langle
g_{t}-\widetilde{G}_{t},z_{t}-x\rangle-\tfrac{\tau_{t}}{2L}\|g_{t}-g_{t-1}\|_{*}^{2}-\eta_{t}V(z_{t-1},z_{t}).$

Note that

$\displaystyle\langle g_{t}-\widetilde{G}_{t},z_{t}-x\rangle$
$\displaystyle=\langle
g_{t}-g_{t-1}-\alpha_{t}(g_{t-1}-g_{t-2}),z_{t}-x\rangle-\langle\delta_{t-1}+\alpha_{t}(\delta_{t-1}-\delta_{t-2}),z_{t}-x\rangle$
$\displaystyle=\langle g_{t}-g_{t-1},z_{t}-x\rangle-\alpha_{t}\langle
g_{t-1}-g_{t-2},z_{t-1}-x\rangle$ $\displaystyle\quad+\alpha_{t}\langle
g_{t-1}-g_{t-2},z_{t}-z_{t-1}\rangle-\langle\delta_{t-1}+\alpha_{t}(\delta_{t-1}-\delta_{t-2}),z_{t}-x\rangle.$

When taking the $\theta_{t}$-weighted sum of the above inequalities for
$t=1,\ldots,k$, noting that $x_{0}=z_{0}$, and using (A.15a), we obtain

$\displaystyle\textstyle\sum_{t=1}^{k}\theta_{t}\left\\{\tau_{t}[f(x_{t})-f(x_{t-1})]-\langle
g_{t},x-x_{t}\rangle\right\\}$
$\displaystyle\leq\theta_{1}\eta_{1}V(x_{0},x)-\theta_{k}\eta_{k}V(z_{k},x)+\theta_{k}\langle
g_{k}-g_{k-1},z_{k}-x\rangle+\Delta_{k},$ (A.18)

where

$\displaystyle\Delta_{k}$
$\displaystyle:=\textstyle\sum_{t=1}^{k}\theta_{t}\big{[}\alpha_{t}\langle
g_{t-1}-g_{t-2},z_{t}-z_{t-1}\rangle-\tfrac{\tau_{t}}{2L}\|g_{t}-g_{t-1}\|_{*}^{2}$
$\displaystyle~{}\qquad\quad\qquad-\langle\delta_{t-1}+\alpha_{t}(\delta_{t-1}-\delta_{t-2}),z_{t}-x\rangle-\eta_{t}V(z_{t-1},z_{t})\big{]}.$

Observe that

$\theta_{t}\alpha_{t}\langle
g_{t-1}-g_{t-2},z_{t}-z_{t-1}\rangle\leq{\theta_{t}^{2}\alpha_{t}^{2}L\over
2\theta_{t-1}\tau_{t-1}}\|z_{t}-z_{t-1}\|^{2}+{\theta_{t-1}\tau_{t-1}\over
2L}\|g_{t-1}-g_{t-2}\|_{*}^{2}$

so that, after rearranging the terms,

$\displaystyle\Delta_{k}$
$\displaystyle\leq\underbrace{\theta_{2}\alpha_{2}\langle
g_{1}-g_{0},z_{2}-z_{1}\rangle+\textstyle\sum_{t=3}^{k}\tfrac{L\theta_{t}^{2}\alpha_{t}^{2}}{2\theta_{t-1}\tau_{t-1}}\|z_{t}-z_{t-1}\|^{2}-\tfrac{\theta_{k}\tau_{k}}{2L}\|g_{k}-g_{k-1}\|_{*}^{2}-\tfrac{1}{5}\textstyle\sum_{t=1}^{k}\theta_{t}\eta_{t}V(z_{t-1},z_{t})}_{=:\Delta_{k,1}}$
$\displaystyle\quad\underbrace{-\textstyle\sum_{t=1}^{k}\theta_{t}\langle\delta_{t-1}+\alpha_{t}(\delta_{t-1}-\delta_{t-2}),z_{t}-x\rangle-\tfrac{4}{5}\textstyle\sum_{t=1}^{k}\theta_{t}\eta_{t}V(z_{t-1},z_{t})}_{=:\Delta_{k,2}}.$

#### 2o.

Let us bound $\Delta_{k,1}$. By the Young inequality,

$\alpha_{2}\langle
g_{1}-g_{0},z_{2}-z_{1}\rangle\leq\frac{5\alpha_{2}^{2}}{2\eta_{2}}\|g_{1}-g_{0}\|_{*}^{2}+\frac{\eta_{2}}{10}\|z_{2}-z_{1}\|^{2}\leq\frac{5\alpha_{2}^{2}}{2\eta_{2}}\|g_{1}-g_{0}\|_{*}^{2}+\frac{\eta_{2}}{5}V(z_{1},z_{2}),$

thus

$\displaystyle\Delta_{k,1}$
$\displaystyle\leq\tfrac{5\theta_{2}\alpha_{2}^{2}\|g_{1}-g_{0}\|_{*}^{2}}{2\eta_{2}}-\tfrac{1}{5}\theta_{1}\eta_{1}V(z_{1},z_{0})+\textstyle\sum_{t=3}^{k}\left(\tfrac{L\theta_{t}^{2}\alpha_{t}^{2}}{2\theta_{t-1}\tau_{t-1}}-\tfrac{\theta_{t}\eta_{t}}{10}\right)\|z_{t}-z_{t-1}\|^{2}-\tfrac{\theta_{k}\tau_{k}}{2L}\|g_{k}-g_{k-1}\|_{*}^{2}$
$\displaystyle\leq-\tfrac{\theta_{k}\tau_{k}}{2L}\|g_{k}-g_{k-1}\|_{*}^{2}$
(A.19)

where the last inequality follows from (A.15b) and the bound

$\displaystyle\tfrac{5\theta_{2}\alpha_{2}^{2}\|g_{1}-g_{0}\|_{*}^{2}}{2\eta_{2}}\leq\tfrac{5\theta_{2}\alpha_{2}^{2}L^{2}\|x_{1}-x_{0}\|^{2}}{2\eta_{2}}=\tfrac{5\theta_{2}\alpha_{2}^{2}L^{2}\|z_{1}-z_{0}\|^{2}}{2\eta_{2}}\leq\tfrac{5\theta_{2}\alpha_{2}^{2}L^{2}V(z_{1},z_{0})}{\eta_{2}}\leq\tfrac{\theta_{1}\eta_{1}V(z_{1},z_{0})}{5}.$

Note that

$\displaystyle-\eta_{k}V(z_{k},x)+\langle
g_{k}-g_{k-1},z_{k}-x\rangle-\tfrac{\tau_{k}}{2L}\|g_{k}-g_{k-1}\|_{*}^{2}$
$\displaystyle\qquad\leq-\tfrac{\eta_{k}}{2}\|z_{k}-x\|^{2}+\langle
g_{k}-g_{k-1},z_{k}-x\rangle-\tfrac{\tau_{k}}{2L}\|g_{k}-g_{k-1}\|_{*}^{2}$
$\displaystyle\qquad\leq-(\tfrac{\eta_{k}}{2}-\tfrac{L}{2\tau_{k}})\|z_{k}-x\|^{2}\leq
0$

due to the second relationship in (A.15c). Now, substituting the bound (A.5)
into (A.18) results in

$\displaystyle\textstyle\sum_{t=1}^{k}\theta_{t}\left\\{\tau_{t}[f(x_{t})-f(x_{t-1})]-\langle
g_{t},x-x_{t}\rangle\right\\}$
$\displaystyle\qquad\leq\theta_{1}\eta_{1}V(x_{0},x)-\theta_{k}\eta_{k}V(z_{k},x)+\theta_{k}\langle
g_{k}-g_{k-1},z_{k}-x\rangle-\tfrac{\theta_{k}\tau_{k}}{2L}\|g_{k}-g_{k-1}\|_{*}^{2}+\Delta_{k,2}$
$\displaystyle\qquad\leq\theta_{1}\eta_{1}V(x_{0},x)+\Delta_{k,2},$ (A.20)

#### 3o.

To bound $\Delta_{k,2}$ we act as follows. Observe that

$\displaystyle\Delta_{k,2}$
$\displaystyle\leq-\textstyle\sum_{t=1}^{k}\theta_{t}[\langle\delta_{t-1},z_{t}-z_{t-1}\rangle+\langle\delta_{t-1},z_{t-1}-x\rangle]$
$\displaystyle\quad-\textstyle\sum_{t=1}^{k}\theta_{t}\alpha_{t}[\langle\delta_{t-1},z_{t}-z_{t-1}\rangle+\langle\delta_{t-1},z_{t-1}-x\rangle]$
$\displaystyle\quad+\textstyle\sum_{t=2}^{k}\theta_{t}\alpha_{t}[\langle\delta_{t-2},z_{t}-z_{t-1}\rangle+\langle\delta_{t-2},z_{t-1}-z_{t-2}\rangle+\langle\delta_{t-2},z_{t-2}-x\rangle]$
$\displaystyle\quad-\tfrac{2}{5}\textstyle\sum_{t=1}^{k}\theta_{t}\eta_{t}\|z_{t-1}-z_{t}\|^{2}.$

Recall that $z_{t}$ is ${\cal F}_{t-1}$-measurable. When using the bound

$\langle\delta_{s},z_{r}-z_{r-1}\rangle-a\|z_{t-1}-z_{t}\|^{2}\leq\|\delta_{s}\|_{*}^{2}/(4a)$

for $s<r$ and $a>0$, we get

$\displaystyle\Delta_{k,2}$
$\displaystyle\leq\textstyle\sum_{t=1}^{k}[\tfrac{5\theta_{t}}{2\eta_{t}}\|\delta_{t-1}\|_{*}^{2}-\theta_{t}\langle\delta_{t-1},z_{t-1}-x\rangle]$
$\displaystyle\quad+\textstyle\sum_{t=1}^{k}[\tfrac{5\theta_{t}\alpha_{t}^{2}}{2\eta_{t}}\|\delta_{t-1}\|_{*}^{2}-\theta_{t}\alpha_{t}\langle\delta_{t-1},z_{t-1}-x\rangle]$
$\displaystyle\quad+\textstyle\sum_{t=2}^{k}[\tfrac{5\theta_{t}\alpha_{t}^{2}}{2\eta_{t}}\|\delta_{t-2}\|_{*}^{2}+\tfrac{5\theta_{t}^{2}\alpha_{t}^{2}}{2\theta_{t-1}\eta_{t-1}}\|\delta_{t-2}\|_{*}^{2}+\theta_{t}\alpha_{t}\langle\delta_{t-2},z_{t-2}-x\rangle]$
$\displaystyle\leq\textstyle\sum_{t=1}^{k}\tfrac{5}{2}[\tfrac{\theta_{t}(1+\alpha_{t}^{2})}{\eta_{t}}+\tfrac{\theta_{t+1}\alpha_{t+1}^{2}}{\eta_{t+1}}+\tfrac{\theta_{t+1}^{2}\alpha_{t+1}^{2}}{\theta_{t}\eta_{t}}]\|\delta_{t-1}\|_{*}^{2}$
$\displaystyle\quad-\textstyle\sum_{t=1}^{k}\theta_{t}(1+\alpha_{t})\langle\delta_{t-1},z_{t-1}-x\rangle+\textstyle\sum_{t=1}^{k-1}\theta_{t}\langle\delta_{t-1},z_{t-1}-x\rangle$
$\displaystyle\leq\textstyle\sum_{t=1}^{k}\tfrac{5}{2}[\tfrac{\theta_{t}(1+\alpha_{t}^{2})}{\eta_{t}}+\tfrac{\theta_{t+1}\alpha_{t+1}^{2}}{\eta_{t+1}}+\tfrac{\theta_{t+1}^{2}\alpha_{t+1}^{2}}{\theta_{t}\eta_{t}}]\|\delta_{t-1}\|_{*}^{2}$
$\displaystyle\quad-\textstyle\sum_{t=0}^{k-1}\theta_{t}\langle\delta_{t},z_{t}-x\rangle-\theta_{k}\langle\delta_{k-1},z_{k-1}-x\rangle.$
(A.21)

When substituting the bound (A.21) for $\Delta_{k,2}$ into (A.20) we obtain
(A.17). $\square$

#### Proof of the theorem.

By taking expectation on both sides of (A.17), and using (SN), (2.2), and the
fact that

$\displaystyle\langle g(x_{t}),x-x_{t}\rangle\leq f(x)-f(x_{t})$

we have for all $x\in X$,

$\displaystyle\textstyle\sum_{t=1}^{k}\theta_{t}{\mathbb{E}}\left\\{\tau_{t}[f(x_{t})-f(x_{t-1})]-[f(x)-f(x_{t})]\right\\}$
$\displaystyle\leq\theta_{1}\eta_{1}V(x_{0},x)+\textstyle\sum_{t=1}^{k}\epsilon_{t-1}\left\\{{\cal
L}{\mathbb{E}}[f(x_{t-1})-f^{*}-\langle
g(x^{*}),x_{t-1}-x^{*}\rangle]+\sigma_{*}^{2}\right\\}.$

Note that (4.4) implies that

$\theta_{t}\tau_{t}+{\cal
L}\epsilon_{t-1}\leq\theta_{t-1}(1+\tau_{t-1}),~{}t\geq 2.$ (A.22)

Thus, using $\langle g(x^{*}),x_{t}-x^{*}\rangle\geq 0$ and
$f(x_{0})-f^{*}-\langle
g(x^{*}),x_{0}-x^{*}\rangle\leq\tfrac{L}{2}\|x_{0}-x^{*}\|^{2}$, for $x=x^{*}$
we obtain

$\displaystyle\theta_{k}(1+\tau_{t}){\mathbb{E}}[f(x_{k})-f^{*}]$
$\displaystyle\leq\textstyle\sum_{t=1}^{k}\theta_{t}(1+\tau_{t}){\mathbb{E}}[f(x_{t})-f^{*}]-\textstyle\sum_{t=1}^{k-1}\theta_{t}(1+\tau_{t}){\mathbb{E}}[f(x_{t})-f^{*}]$
$\displaystyle\leq\textstyle\sum_{t=1}^{k}\theta_{t}(1+\tau_{t}){\mathbb{E}}[f(x_{t})-f^{*}]-\textstyle\sum_{t=2}^{k}(\theta_{t}\tau_{t}+\epsilon_{t-1}{\cal
L}){\mathbb{E}}[f(x_{t-1})-f^{*}]$
$\displaystyle\qquad\leq\theta_{1}\eta_{1}V(x_{0},x^{*})+\tfrac{\epsilon_{0}{\cal
L}L}{2}\|x_{0}-x^{*}\|^{2}+\textstyle\sum_{t=1}^{k}\epsilon_{t-1}\sigma_{*}^{2}$

which is (4.5). $\square$

###  A.6 Proof of Corollary 3

#### 1o.

Note that in the premise of the corollary one has $\tau_{t}={t-1\over 3}$. Let
us check that with the present choice of stepsize parameters conditions
(A.15a)–(A.15c) and (A.22) (which are equivalents (4.2) and (4.4)) are
satisfied. It is easy to see that (A.15a) holds. Because $\eta\geq 30L$, we
have

$\displaystyle\eta_{1}\eta_{2}$
$\displaystyle{\color[rgb]{0,0,0}=}\tfrac{1}{2}\cdot(30L)^{2}\geq
25\alpha_{2}L^{2},$ $\displaystyle\eta_{t}\tau_{t-1}$
$\displaystyle{\color[rgb]{0,0,0}=}\tfrac{\eta}{t}\cdot\tfrac{t-2}{3}\geq\tfrac{10L(t-2)}{t}\geq
5L\alpha_{t},~{}t=3,...,k$ $\displaystyle\eta_{k}\tau_{k}$
$\displaystyle\geq\tfrac{\eta}{k}\cdot\tfrac{k-1}{3}\geq\tfrac{10L(k-1)}{k}\geq
L.$

To check (A.22), notice that for $t\geq 1$,

$\displaystyle
q_{t}=\tfrac{\theta_{t+1}(1+\alpha_{t+1}^{2})}{\eta_{t+1}}+\tfrac{\theta_{t+2}\alpha_{t+2}^{2}}{\eta_{t+2}}+\tfrac{\theta_{t+2}^{2}\alpha_{t+2}^{2}}{\theta_{t+1}\eta_{t+1}}=\tfrac{3(t+1)^{2}+t^{2}}{\eta}\leq\tfrac{4(t+1)^{2}}{\eta},$

thus

$\displaystyle\epsilon_{t}=\tfrac{5\Omega
q_{t}}{2m_{t}}\leq\tfrac{10\Omega(t+1)^{2}}{\eta m}.$

Given the above inequality, for $t\geq 3$ we have

$\displaystyle\theta_{t}\tau_{t}+{\cal
L}\epsilon_{t-1}=t\cdot\tfrac{t-1}{3}+{\cal
L}\epsilon_{t-1}\leq\tfrac{t(t-1)}{3}+\tfrac{10\Omega t^{2}{\cal L}}{\eta
m}\leq\tfrac{t(t-1)}{3}+\tfrac{t-1}{3}=\tfrac{(t+1)(t-1)}{3},$

where the second inequality follows from $\eta\geq\tfrac{30\Omega(k+2){\cal
L}}{m}$. Combining the above bound with
$\theta_{t-1}(1+\tau_{t-1})=\tfrac{(t-1)(t+1)}{3}$, we arrive at

$\displaystyle\theta_{t-1}(1+\tau_{t-1})-(\theta_{t}\tau_{t}+{\cal
L}\epsilon_{t-1})\geq 0,\qquad t\geq 2,$

which is (A.22). We conclude that the bound (4.5) of Theorem (4.1) holds.

On the other hand, when
$\eta\geq\sqrt{\tfrac{10\Omega(k+1)^{3}\sigma_{*}^{2}}{3mD^{2}}}$ we have

$\displaystyle\textstyle\sum_{t=0}^{k-1}\epsilon_{t}\sigma_{*}^{2}\leq\textstyle\sum_{t=0}^{k-1}\tfrac{10\Omega(t+1)^{2}\sigma_{*}^{2}}{\eta
m}\leq\tfrac{10\Omega(k+1)^{3}\sigma_{*}^{2}}{3\eta
m}\leq\sqrt{\tfrac{10\Omega(k+1)^{3}\sigma_{*}^{2}D^{2}}{3m}}.$

When substituting the above bound into (4.5) and noticing that
$\theta_{k}(1+\tau_{k})=\tfrac{k(k+2)}{3}$ and ${\cal
L}\epsilon_{0}\leq\tfrac{1}{3}$, we conclude that

$\displaystyle{\mathbb{E}}[f(x_{k})-f^{*}]\leq\tfrac{91LD^{2}}{k(k+2)}+\tfrac{90\Omega{\cal
L}D^{2}}{mk}+\sqrt{\tfrac{120\Omega\sigma_{*}^{2}D^{2}}{mk}},$

which completes the proof of (4.6).

#### 2o.

The “Furthermore” part of the statement immediately follows from (1.11) and
the fact that when $\eta\geq\sqrt{\tfrac{20(k+1)^{3}\sigma_{*}^{2}}{3mR^{2}}}$
one has

$\textstyle\sum_{t=0}^{k-1}\epsilon_{t}\sigma_{*}^{2}\leq\tfrac{10\Omega(k+1)^{3}\sigma_{*}^{2}}{3\eta
m}\leq\sqrt{\tfrac{5\Omega^{2}(k+1)^{3}\sigma_{*}^{2}R^{2}}{3m}}.$ $\square$

###  A.7 Proofs of Corollaries 4 and Corollary 5

#### Proof of Corollary 4.

Note that bound (4.8) of Corollary 3 implies that whenever size $m$ of the
mini-batch satisfies

$m\geq\max\left\\{1,\tfrac{18(k+2)\Omega{\cal
L}}{L},\tfrac{15N(N+2)^{2}\sigma_{*}^{2}}{L^{2}R^{2}}\right\\},$

we have for the approximate solution $x_{N}$ by SGE after $N$ iterations,

$\displaystyle{\mathbb{E}}[f(x_{N})-f^{*}]\leq\tfrac{91L\Omega
R^{2}}{2N(N+2)}+\tfrac{45\Omega^{2}{\cal
L}R^{2}}{mk}+\sqrt{\tfrac{60\Omega^{2}\sigma_{*}^{2}R^{2}}{mk}}\leq\tfrac{91L\Omega
R^{2}}{2N(N+2)}+\tfrac{5L\Omega R^{2}}{2N(N+2)}+\tfrac{2L\Omega
R^{2}}{N(N+2)}=\tfrac{50L\Omega R^{2}}{N(N+2)}$ (A.23)

where $R$ is an upper bound for the “initial distance to $x^{*}$.”

Note that $\|{y}^{0}-x^{*}\|\leq R_{0}$. Let us now assume that for $1\leq
k\leq K$, $\|{y}^{k-1}-x^{*}\|\leq R_{k}$, so that at the beginning of the
$k$th stage $\|x_{0}-x^{*}\|\leq R_{k}$. Based on the above inequality, by the
choice of $N$,

$\displaystyle{\mathbb{E}}[f(x_{N})-f^{*}]\leq\tfrac{50L\Omega
R_{k-1}^{2}}{N(N+2)}\leq\tfrac{\mu R_{k-1}^{2}}{4}=2^{-k-1}\mu
R_{0}^{2}=\tfrac{1}{2}\mu R_{k}^{2}.$

Due to (5.1) this also means that

${\mathbb{E}}[\|x_{N}-x^{*}\|^{2}]\leq{R_{k}^{2}}=2^{-k}R_{0}^{2}.$ $\square$

#### Proof of Corollary 5.

We have $\|\overline{y}^{0}-x^{*}\|\leq R_{0}$. Let us assume that
$\|\overline{y}^{k-1}-x^{*}\|\leq R_{k-1}$ for some $1\leq k\leq K$. From
(A.23) we conclude that

$\displaystyle{\mathbb{E}}[f({y}^{k})-f^{*}]\leq\tfrac{50L\Omega
R_{k-1}^{2}}{N(N+2)}\leq\tfrac{\underline{\kappa}R_{k}^{2}}{16s}={2^{-k-4}s^{-1}\underline{\kappa}R_{0}^{2}}$

by the choice of $N$. Thus, recalling that
$\overline{y}^{k}=\mathsf{sparse}({y}^{k})$, by (6.7) we conclude that

$\displaystyle{\mathbb{E}}[\|\overline{y}^{k}-x^{*}\|^{2}]\leq
8s\|{y}^{k}-x^{*}\|^{2}_{2}\leq
16s\underline{\kappa}^{-1}{\mathbb{E}}[f({y}^{k})-f^{*}]\leq
R_{k}^{2}=2^{-k}R_{0}^{2}$

which completes the proof. $\square$

## Reference

  * [1] A. Agarwal, S. Negahban, and M. J. Wainwright.  Stochastic optimization and sparse statistical recovery: Optimal algorithms for high dimensions.  Advances in Neural Information Processing Systems, 25, 2012.
  * [2] P. J. Bickel, Y. Ritov, and A. B. Tsybakov.  Simultaneous analysis of lasso and dantzig selector.  2009\.
  * [3] T. Blumensath and M. E. Davies.  Iterative hard thresholding for compressed sensing.  Applied and computational harmonic analysis, 27(3):265–274, 2009\.
  * [4] E. Candes and B. Recht.  Exact matrix completion via convex optimization.  Communications of the ACM, 55(6):111–119, 2012.
  * [5] E. Candes and T. Tao.  The dantzig selector: Statistical estimation when p is much larger than n.  2007\.
  * [6] E. J. Candes and Y. Plan.  Near-ideal model selection by $\ell_{1}$ minimization.  Annals of Statistics, 37(5A):2145–2177, 2009.
  * [7] E. J. Candes and Y. Plan.  A probabilistic and ripless theory of compressed sensing.  IEEE transactions on information theory, 57(11):7235–7254, 2011\.
  * [8] E. J. Candes and Y. Plan.  Tight oracle inequalities for low-rank matrix recovery from a minimal number of noisy random measurements.  IEEE Transactions on Information Theory, 57(4):2342–2359, 2011\.
  * [9] E. J. Candes, J. K. Romberg, and T. Tao.  Stable signal recovery from incomplete and inaccurate measurements.  Communications on Pure and Applied Mathematics: A Journal Issued by the Courant Institute of Mathematical Sciences, 59(8):1207–1223, 2006.
  * [10] A. Cotter, O. Shamir, N. Srebro, and K. Sridharan.  Better mini-batch algorithms via accelerated gradient methods.  Advances in neural information processing systems, 24, 2011.
  * [11] M. Fazel, E. Candes, B. Recht, and P. Parrilo.  Compressed sensing and robust recovery of low rank matrices.  In 2008 42nd Asilomar Conference on Signals, Systems and Computers, pages 1043–1047. IEEE, 2008.
  * [12] P. Gaillard and O. Wintenberger.  Sparse accelerated exponential weights.  In Artificial Intelligence and Statistics, pages 75–82. PMLR, 2017\.
  * [13] S. Ghadimi and G. Lan.  Optimal stochastic approximation algorithms for strongly convex stochastic composite optimization i: A generic algorithmic framework.  SIAM Journal on Optimization, 22(4):1469–1492, 2012.
  * [14] S. Ghadimi and G. Lan.  Optimal stochastic approximation algorithms for strongly convex stochastic composite optimization, ii: shrinking procedures and optimal algorithms.  SIAM Journal on Optimization, 23(4):2061–2089, 2013.
  * [15] S. Ilandarideva, Y. Bekri, A. Iouditski, and V. Perchet.  Stochastic mirror descent for large-scale sparse recovery.  In International Conference on Artificial Intelligence and Statistics, pages 5931–5957. PMLR, 2023.
  * [16] P. Jain, A. Tewari, and P. Kar.  On iterative hard thresholding methods for high-dimensional m-estimation.  Advances in neural information processing systems, 27, 2014.
  * [17] A. Juditsky, A. Kulunchakov, and H. Tsyntseus.  Sparse recovery by reduced variance stochastic approximation.  Information and Inference: A Journal of the IMA, 12(2):851–896, 2023.
  * [18] A. Juditsky and A. Nemirovski.  Accuracy guarantees for $\ell_{1}$-recovery.  IEEE Transactions on Information Theory, 57(12):7818–7839, 2011\.
  * [19] A. Juditsky, A. Nemirovski, et al.  First order methods for nonsmooth convex large-scale optimization, i: general purpose methods.  Optimization for Machine Learning, 30(9):121–148, 2011.
  * [20] A. Juditsky and Y. Nesterov.  Deterministic and stochastic primal-dual subgradient algorithms for uniformly convex minimization.  Stochastic Systems, 4(1):44–80, 2014.
  * [21] G. Kotsalis, G. Lan, and T. Li.  Simple and optimal methods for stochastic variational inequalities, i: Operator extrapolation.  SIAM Journal on Optimization, 32(3):2041–2073, 2022.
  * [22] G. Kotsalis, G. Lan, and T. Li.  Simple and optimal methods for stochastic variational inequalities, ii: Markovian noise and policy evaluation in reinforcement learning.  SIAM Journal on Optimization, 32(2):1120–1155, 2022.
  * [23] S. Kotz and S. Nadarajah.  Multivariate t-distributions and their applications.  Cambridge University Press, 2004.
  * [24] G. Lan.  An optimal method for stochastic composite optimization.  Mathematical Programming, 133(1-2):365–397, 2012.
  * [25] G. Lan.  First-order and stochastic optimization methods for machine learning, volume 1.  Springer, 2020.
  * [26] G. Lan and Y. Zhou.  Random gradient extrapolation for distributed and stochastic optimization.  SIAM Journal on Optimization, 28(4):2753–2782, 2018.
  * [27] T. Li, G. Lan, and A. Pananjady.  Accelerated and instance-optimal policy evaluation with linear function approximation.  arXiv preprint arXiv:2112.13109, 2021.
  * [28] T. Li, F. Wu, and G. Lan.  Stochastic first-order methods for average-reward markov decision processes.  arXiv preprint arXiv:2205.05800, 2022.
  * [29] C. Liu and M. Belkin.  Mass: an accelerated stochastic method for over-parametrized learning.  arXiv preprint arXiv:1810.13395, 2018.
  * [30] H. Liu and R. Foygel Barber.  Between hard and soft thresholding: optimal iterative thresholding algorithms.  Information and Inference: A Journal of the IMA, 9(4):899–933, 2020\.
  * [31] A. Nemirovski, A. Juditsky, G. Lan, and A. Shapiro.  Robust stochastic approximation approach to stochastic programming.  SIAM Journal on Optimization, 19(4):1574–1609, 2009.
  * [32] A. Nemirovski and D. Yudin.  Complexity of problems and effectiveness of methods of optimization(russian book).  Moscow, Izdatel’stvo Nauka, 1979. 384, 1979.
  * [33] A. S. Nemirovski and D. Yudin.  Problem complexity and method efficiency in optimization.  Wiley-Interscience Series in Discrete Mathematics. John Wiley, XV, 1983\.
  * [34] A. Nemirovskii and Y. Nesterov.  Optimal methods of smooth convex minimization.  USSR Computational Mathematics and Mathematical Physics, 25:21–30, 1985.
  * [35] A. Nemirovskii and D. Yudin.  Information-based complexity of mathematical programming.  Engineering Cybernetics, 1:76–100, 1983.
  * [36] Y. Nesterov.  A method for unconstrained convex minimization problem with the rate of convergence $O(1/k^{2})$.  In Doklady an ussr, volume 269, pages 543–547, 1983.
  * [37] Y. Nesterov and A. Nemirovski.  On first-order algorithms for $\ell_{1}$/nuclear norm minimization.  Acta Numerica, 22:509–575, 2013.
  * [38] B. T. Polyak.  New stochastic approximation type procedures.  Automat. i Telemekh, 7(98-107):2, 1990.
  * [39] B. T. Polyak and A. B. Juditsky.  Acceleration of stochastic approximation by averaging.  SIAM journal on control and optimization, 30(4):838–855, 1992.
  * [40] S. Shalev-Shwartz and A. Tewari.  Stochastic methods for $\ell_{1}$ regularized loss minimization.  In Proceedings of the 26th Annual International Conference on Machine Learning, pages 929–936, 2009.
  * [41] N. Srebro, K. Sridharan, and A. Tewari.  Smoothness, low noise and fast rates.  Advances in neural information processing systems, 23, 2010.
  * [42] A. Tsybakov, V. Koltchinskii, and K. Lounici.  Nuclear-norm penalization and optimal rates for noisy low-rank matrix completion.  Annals of Statistics, 39(5):2302–2329, 2011.
  * [43] B. E. Woodworth and N. Srebro.  An even more optimal stochastic optimization algorithm: minibatching and interpolation learning.  Advances in Neural Information Processing Systems, 34:7333–7345, 2021.
  * [44] D. Yudin and A. Nemirovski.  Computational complexity of strictly convex programming.  Eknomika i Matematicheskie Metody, 3:550–569, 1977.