# Convergence and concentration properties of constant step-size SGD through
Markov chains

Ibrahim Merad
LPSM, UMR 8001, Université Paris Cité, Paris, France Stéphane Gaïffas LPSM,
UMR 8001, Université Paris Cité, Paris, France and DMA, École normale
supérieure

###### Abstract

We consider the optimization of a smooth and strongly convex objective using
constant step-size stochastic gradient descent (SGD) and study its properties
through the prism of Markov chains. We show that, for unbiased gradient
estimates with mildly controlled variance, the iteration converges to an
invariant distribution in total variation distance. We also establish this
convergence in Wasserstein-2 distance under a relaxed assumption on the
gradient noise distribution compared to previous work. Thanks to the
invariance property of the limit distribution, our analysis shows that the
latter _inherits_ sub-Gaussian or sub-exponential concentration properties
when these hold true for the gradient. This allows the derivation of high-
confidence bounds for the final estimate. Finally, under such conditions in
the linear case, we obtain a dimension-free deviation bound for the Polyak-
Ruppert average of a tail sequence. All our results are non-asymptotic and
their consequences are discussed through a few applications.

_Keywords._ stochastic approximation; Markov chains; Polyak-Ruppert averaging;
generalization error

##  1 Introduction

We consider the following stochastic optimization problem

$\min_{\theta\in\mathbb{R}^{d}}\mathcal{L}(\theta):=\mathbb{E}_{\zeta}[\ell(\theta,\zeta)],$
(1)

where $\mathcal{L}$ is a smooth strongly convex objective only accessible
through unbiased random gradient samples
$G(\theta,\zeta)=\nabla\ell(\theta,\zeta)$ which may be queried at any
parameter value $\theta\in\mathbb{R}^{d}$. Given an initial point $\theta_{0}$
and a step-size $\beta,$ problem (1) is commonly solved using the well-known
stochastic gradient descent (SGD) algorithm defined by the iteration

$\theta_{t+1}=\theta_{t}-\beta G(\theta_{t},\zeta_{t}),\quad\text{ for }\quad
t\geq 0.$ (2)

We study the convergence properties of the Markov chain $(\theta_{t})_{t\geq
0}$ generated by the above iteration as well as the concentration properties
satisfied by a derived estimator $\widehat{\theta}$ of the global optimum
$\theta^{\star}=\operatorname*{argmin}_{\theta}\mathcal{L}(\theta)$ based on
the concentration of the gradient samples $G(\theta_{t},\zeta_{t}).$

Problem (1) is the common formulation for a large fraction of statistical
learning problems where the objective $\mathcal{L}(\theta)$ is defined as the
expectation of a loss function $\ell$ over a random variable $\zeta$ following
an unknown distribution of samples. In a practical setting, the random
gradients $G(\theta_{t},\zeta_{t})$ are computed using a dataset of
independent and identically distributed samples $(\zeta_{i})_{i=1}^{n}.$ The
SGD algorithm is employed to solve (1) in two situations. Either the samples
$(\zeta_{i})_{i=1}^{n}$ are available offline but in such a great amount that
using the whole dataset at each gradient step incurs an excessive
computational load, therefore, individual samples or small batches are used at
each iteration instead. Or, the samples $\zeta_{i}$ are received individually
in an online fashion and optimization must be run using one instance at a
time. Our framework covers both cases provided that each iteration uses new
data which is independent from the past.

Thanks to its simplicity and efficiency, the SGD algorithm is widely adopted
as the go-to approach for stochastic optimization problems in general. Since
its first appearance in the seminal work of [73] the theoretical properties of
SGD have been investigated in a series of pioneering works [13, 80, 22]. A
notable milestone in these theoretical developments was the discovery of
Polyak-Ruppert averaging [79, 69] which allows to reduce the impact of noise
and improve the convergence rate for certain cases of interest. The subject
benefited from a growing attention with the advent of complex machine learning
models such as neural networks and a rich literature has appeared to address
the surfacing questions about SGD and its numerous variants and use cases [84,
1, 61, 9, 64].

Although the basic definition of the SGD iteration (2) is quite simple, a
great number of variations are possible by playing on various aspects among
which the choice of step-size is critical. Early work [73] suggested a
decaying step-size of order $t^{-1}$ but this leads to poor dependence on
problem conditioning [1, 22] while other step-size schedules with slower decay
of order $t^{-\alpha}$ with $\alpha\in(1/2,1]$ combined with averaging achieve
better practical and theoretical performance [79, 69, 61]. In this work, we
consider constant step-size SGD which is also a commonly adopted choice due to
its usually fast convergence [82, 87, 54].

###  1.1 Main contributions

This paper studies constant step-size SGD as a Markov chain and makes the
following contributions.

  * •

We state two convergence results of the Markov chain to an invariant
distribution. The first ergodicity theorem states convergence in total
variation distance and the second one in terms of the Wasserstein-$2$
distance. While similar results exist in the literature, our version for the
Wasserstein convergence mode uses a significantly weaker assumption than
previous results.

  * •

We show that sub-Gaussian and sub-exponential concentration of the gradient
samples imply the same property for the invariant limit distribution with
$\Psi_{1}/\Psi_{2}$ constant proportional to the step-size. Thanks to this
property, we obtain high-confidence deviation bounds on the final SGD iterate.

  * •

Provide a slightly stronger concentration assumption on the gradient samples,
we show similar but dimension-free high-confidence bounds on the last SGD
iterate.

  * •

Finally, for the special case of a linear gradient, we obtain a high-
confidence dimension-free bound for the Polyak-Ruppert average of a tail
sequence of the SGD iterates. This is achieved, in part, thanks to a more
generic concentration result which holds for any Lipschitz function applied to
a stationary sequence.

All our results are non-asymptotic.

###  1.2 Related works

#### SGD as a Markov chain.

A fairly limited portion of the SGD literature adopts the Markov chain
approach. Among the earliest, [68] studied the iteration in question for
vanishing step-size, while [2] considers constant step-size averaged SGD for
non-strongly-convex smooth objectives and shows $L_{p}$ convergence of the
excess risk for all $p\geq 1.$ Although their analysis does not use Markov
chain theory, they discuss properties of the invariant distribution which the
iteration converges to, including a few properties we state in this paper.
However, they do not derive high-confidence estimation bounds as we do. More
recently, convergence in Wasserstein distance was established by [18] for
constant step-size SGD applied to a strongly convex and smooth objective,
albeit under a fairly strong co-coercivity condition on each gradient sample.
Further, an expansion of the asymptotic moments of averaged SGD is provided in
[18] and the Richardson-Romberg extrapolation strategy is studied which allows
to reduce the estimation error on the global optimum. Most recently, [91]
studied SGD run on a non-convex, non-smooth but quadratically growing
objective. Under such weakened conditions, they show that the generated Markov
chain is geometrically ergodic (see [58]) and proceed to establish a CLT for
the generated Markov iterates. They also state results controlling the bias of
the limit distribution under additional assumptions such as convexity, $L_{4}$
control of the gradient noise and a generalized Łojasiewicz condition [39].

#### High probability bounds.

In addition to establishing the convergence of SGD in expectation, the works
of [72, 2, 1] make the further step of stating high-confidence bounds on the
final error. Still, sub-Gaussian concentration only holds under strong bounded
gradient assumptions. High-confidence deviation results are also stated in
[26] where an accelerated stochastic optimization method for strongly convex
composite objectives is studied. However, the bounds are sub-exponential while
the gradient is assumed to be sub-Gaussian.

In [38], high probability bounds are proved for the PEGASOS algorithm using
Freedman’s inequality for martingales [23]. A generalization of the said
inequality was used by [32, 31] to prove such bounds for SGD in the non-smooth
strongly convex case. Most recently, for a careful choice of step-size, [36]
obtained high-confidence results on the last SGD iterate. Unfortunately, both
previous works require a deterministic bound to hold over the gradient or its
noise which strongly constrains their applicability. Finally, a high
probability analysis of Delayed AdaGrad with momentum was presented by [48] in
the smooth non-convex setting.

Note that certain recent works design _robust_ variants of SGD achieving sub-
Gaussian deviation bounds on the last iterate with only a second moment
assumption on the gradient [85, 27]. Similar results were later obtained under
even weaker gradient moment assumptions [81, 67]. However, in this work, we
focus on the _classical_ SGD algorithm and the properties inherited by its
iterates from the gradient samples.

#### Polyak-Ruppert averaging.

The averaging procedure introduced by [69, 79] was also studied by [30, 14]
who proved asymptotic convergence properties. Non-asymptotic results and
additional developments appeared in the works of [61, 17, 19, 35, 34, 46] with
particular attention to least-squares, logistic regression and kernel-based
methods. In particular, non-asymptotic results of convergence in expectation
were obtained for averaged SGD in [19, 46, 66, 72]. Among such results, some
demonstrate the advantages of special averaging schemes [84, 45]. Finally,
some relatively recent works obtained high probability concentration bounds
for Polyak Ruppert averaging with and without sub-Gaussian assumptions on the
data [60, 49].

###  1.3 Paper organization

Section 2 lays out the basic setting and assumptions necessary for SGD
convergence. Section 3 states our first SGD ergodicity result. In Section 4,
we first state a basic result on the invariant measure’s expectation, bias and
variance and proceed to derive concentration properties based on analogous
assumptions on the gradient. Section 5 presents an additional convergence
result in Wasserstein distance. In Section 6, we give deviation bounds on the
final SGD iterate which follow from preceding results. We also formulate our
high-confidence bound on a tail Polyak-Ruppert average for the linear case.
Finally, we discuss a few applications in Section 7 and conclude.

##  2 Setting and notations

Let $\Theta$ denote either a convex subset of $\mathbb{R}^{d}$ or
$\mathbb{R}^{d}$ itself depending on context. We refer to the Borel
$\sigma$-algebra of $\mathbb{R}^{d}$ as $\mathcal{B}(\mathbb{R}^{d}).$ For any
random variable $X,$ we denote $\mathcal{D}(X)$ its distribution. We refer to
the space of square-integrable measures on $\mathbb{R}^{d}$ as
$\mathcal{P}_{2}(\mathbb{R}^{d}).$ We denote $\mathcal{M}_{1}(\mathbb{R}^{d})$
the set of probability measures over $\mathbb{R}^{d}.$ For real numbers $a$
and $b,$ we denote $\min(a,b)=a\wedge b$ and $\max(a,b)=a\vee b.$ We denote
$\operatorname{Lip}(\mathcal{X})$ the set of $1$-Lipschitz functions
$h:\mathcal{X}\to\mathbb{R}.$ For $p\in\mathbb{N}^{*},$ we denote
$\|X\|_{L_{p}}=(\mathbb{E}|X|^{p})^{1/p}$ the $L_{p}$ norm of a random
variable $X.$

In the entirety of this work, we assume that $\mathcal{L}$ satisfies

######  Assumption 1.

There exist positive constants $0<\mu\leq L<+\infty$ such that

$\frac{\mu}{2}\|\theta-\theta^{\prime}\|^{2}\leq\mathcal{L}(\theta)-\mathcal{L}(\theta^{\prime})-\langle\nabla\mathcal{L}(\theta^{\prime}),\theta-\theta^{\prime}\rangle\leq\frac{L}{2}\|\theta-\theta^{\prime}\|^{2}$

for all $\theta,\theta^{\prime}\in\mathbb{R}^{d},$ i.e. $\mathcal{L}$ is $L$
gradient-Lipschitz and $\mu$-strongly convex.

As an immediate consequence, $\mathcal{L}$ admits a unique minimum
$\theta^{\star}$ which is a critical point:

$\theta^{\star}=\operatorname*{argmin}_{\theta\in\mathbb{R}^{d}}\mathcal{L}(\theta)\quad\text{such
that}\quad\nabla\mathcal{L}(\theta^{\star})=0.$ (3)

For an initial $\theta_{0}\in\Theta,$ step-size $\beta>0$ and all $t\geq 0,$
we recall the basic SGD iteration

$\theta_{t+1}=\theta_{t}-\beta G(\theta_{t},\zeta_{t}).$ (4)

In this work, we consider constant step-size SGD so that $\beta$ is fixed
along the iteration. We require some basic assumptions on the samples
$G(\theta_{t},\zeta_{t})$ in order to prove the convergence of SGD. Namely,
$G(\theta_{t},\zeta_{t})$ needs to be an unbiased estimator of the true
gradient $\nabla\mathcal{L}(\theta_{t})$ with controlled variance as we
formally state in

######  Assumption 2.

Given any parameter $\theta\in\Theta$, the random gradient sample
$G(\theta,\zeta)$ can be written as

$G(\theta,\zeta)=\nabla\mathcal{L}(\theta)+\varepsilon_{\zeta}(\theta),$ (5)

where $\varepsilon_{\zeta}(\theta)$ is a centered noise i.e.
$\mathbb{E}\big{[}\varepsilon_{\zeta}(\theta)|\theta\big{]}=0$ whose
distribution can be written as
$\mathcal{D}(\varepsilon_{\zeta}(\theta))=\delta\nu_{\theta,1}+(1-\delta)\nu_{\theta,2}$
with $\delta>0$ and $\nu_{\theta,1},\nu_{\theta,2}$ two probability
distributions over $\mathbb{R}^{d}$ such that $\nu_{\theta,1}$ admits a
density $h(\theta,\cdot)$ w.r.t. Lebesgue’s measure satisfying:

$\inf_{\omega\in S}h(\theta,\omega)>0\quad\text{for all $\theta$ and compact
$S\subset\mathbb{R}^{d}$}.$

Moreover, there are positive constants $L_{\sigma}$ and $\sigma^{2}$ such that
for all $\theta$ we have

$\mathbb{E}\big{[}\|\varepsilon_{\zeta}(\theta)\|^{2}\>|\>\theta\big{]}=\mathbb{E}\big{[}\|G(\theta,\zeta)-\nabla\mathcal{L}(\theta)\|^{2}\>|\>\theta\big{]}\leq
L_{\sigma}\|\theta-\theta^{\star}\|^{2}+\sigma^{2}.$ (6)

The additional assumptions on the distribution of the noise
$\varepsilon_{\zeta}(\theta)$ are needed in order to establish the ergodicity
of the resulting Markov chain $(\theta_{t})_{t\geq 0}$ (Theorem 1 below) by
ensuring that the associated transition kernel satisfies a _minorization_
property implying that the chain will sufficiently explore the state space,
see [58] for more details. Note also that these requirements are fairly mild:
they only require that the noise distribution admits a diffuse component.

##  3 Markov Chain and Geometric Ergodicity

Before stating the convergence result for the SGD Markov chain we introduce
some further useful notation. For a given step-size $\beta>0,$ we will denote
$P_{\beta}$ the Markov transition kernel governing the Markov chain
$(\theta_{t})_{t\geq 0}$ generated by iteration (4), so that for any $t\geq 0$
and $A\in\mathcal{B}(\mathbb{R}^{d})$ we have:

$\mathbb{P}(\theta_{t+1}\in A\>|\>\theta_{t})=P_{\beta}(\theta_{t},A).$

The transition kernel $P_{\beta}$ acts on probability distributions
$\nu\in\mathcal{M}_{1}(\mathbb{R}^{d})$ through the mapping $\nu\to\nu
P_{\beta}$ which is defined, for all $A\in\mathcal{B}(\mathbb{R}^{d})$, by
$\nu
P_{\beta}(A)=\int_{A}P_{\beta}(\theta,A)d\nu(\theta)=\mathbb{P}(\theta_{t+1}\in
A\>|\>\theta_{t}\sim\nu).$ For $n\geq 1,$ we similarly define the multi-step
transition kernel $P_{\beta}^{n}$ which is such that
$P_{\beta}^{n}(\theta_{t},A)=\mathbb{P}(\theta_{t+n}\in A\>|\>\theta_{t})$ and
acts on probability distributions $\nu\in\mathcal{M}_{1}(\mathbb{R}^{d})$
through $\nu P_{\beta}^{n}=(\nu P_{\beta})P_{\beta}^{n-1}.$ Finally, we define
the total-variation norm of a signed measure $\nu$ as

$2\|\nu\|_{\operatorname{TV}}=\sup_{f:|f|\leq 1}\int
f(\theta)\nu(d\theta)=\sup_{A\in\mathcal{B}(\mathbb{R}^{d})}\nu(A)-\inf_{A\in\mathcal{B}(\mathbb{R}^{d})}\nu(A).$

In particular, we recover the total-variation _distance_ between two
probability distributions $\nu_{1},\nu_{2}\in\mathcal{M}_{1}(\mathbb{R}^{d})$
as
$d_{\operatorname{TV}}(\nu_{1},\nu_{2})=\|\nu_{1}-\nu_{2}\|_{\operatorname{TV}}.$
We are now ready to state the geometric ergodicity result for the SGD Markov
chain.

######  Theorem 1.

Under Assumptions 1 and 2, the Markov chain $(\theta_{t})_{t\geq 0}$ defined
by iteration (4) with step-size

$\beta<\frac{2\mu}{\mu^{2}+(\mu L\vee L_{\sigma})}$

converges geometrically to a unique invariant measure $\pi_{\beta}$. Namely,
for any initial $\theta_{0}\in\mathbb{R}^{d},$ there exist $\rho<1$ and
$M<+\infty$ such that

$\big{\|}\delta_{\theta_{0}}P^{n}-\pi_{\beta}\big{\|}_{\operatorname{TV}}\leq
M\rho^{n}\big{(}1+\|\theta_{0}-\theta^{\star}\|^{2}\big{)},$ (7)

where $\delta_{\theta_{0}}$ is the Dirac measure located at $\theta_{0}.$

The proof of Theorem 1 is given in Section B.1 and is based on [58, Theorem
15.0.1] and a _drift_ condition in terms of a Lyapunov function. A similar
method was previously used in [91] to establish the convergence of SGD for
non-convex, non-smooth objectives with quadratic growth. However, the focus in
[91] is on proving a central limit theorem for the Markov sequence
$(\theta_{t})_{t\geq 0}$ and bounding the invariant distribution’s bias under
certain regularity conditions. In contrast, we focus on proving concentration
properties of the invariant distribution $\pi_{\beta}$ which allow us to
obtain non-asymptotic deviation bounds on the estimation of the optimum
$\theta^{\star}.$

One of the main questions which arise after the statement of Theorem 1 is
regarding the speed of the claimed convergence. It is clear that the
contraction factor $\rho$ and the constant $M$ mainly depend on the initial
state $\theta_{0}$ and the step-size $\beta,$ however, a precise
quantification is lacking.

The issue is closely related to renewal theory and can be traced back to
Kendall’s theorem [40] and more generally concerns Markov chains satisfying a
drift property. A rich literature [78, 74, 52, 77, 83, 75, 20, 4, 59, 37, 5].
Near optimal results were obtained for stochastically ordered Markov Processes
[52, 75, 53, 24]. Other examples especially amenable to such analysis include
reversible Markov chains [16, 15, 37, 76] and chains satisfying special
assumptions on their renewal distribution [6, 51, 41]. However, the SGD Markov
chain does not satisfy such criteria. An estimation of $\rho$ may be obtained
using results based on renewal theory and Kendall’s theorem [40, 4, 5].
However, the resulting estimations are notoriously pessimistic [70]. Indeed,
let $\alpha:=1-\beta\mu$ be the contraction factor in the absence of gradient
noise (i.e. simple gradient descent) so that we have

$\|\theta_{t}-\beta\nabla\mathcal{L}(\theta_{t})-\theta^{\star}\|\leq\alpha\|\theta_{t}-\theta^{\star}\|\quad\text{for
all }\quad t\geq 0,$

then the worst-case bound on $\rho$ obtained thanks to [4, 5] is such that
$1-\rho\sim(\beta\mu)^{3}$ which is far worse than the intuitive expectation
that $\rho\approx\alpha.$ It is unclear whether the previous estimation can be
improved through a careful study of the renewal properties of the SGD Markov
chain or if a different approach such as the study of the spectral properties
of the transition kernel $P_{\beta}$ is more appropriate. Nevertheless, we
will see in Section 5 below that $\rho$ can be estimated close to $\alpha$
under additional conditions by leveraging Wasserstein convergence.

##  4 Invariant Distribution Properties

This section states that concentration properties of the random gradient
samples used for SGD _transfer_ to the invariant distribution $\pi_{\beta}$
towards which iteration (4) converges as stated in Theorem 1. We begin with a
basic statement which holds without additional assumptions and bounds the
variance of $\pi_{\beta}$ and its bias w.r.t. the true optimum
$\theta^{\star}.$

######  Proposition 1.

In the setting of Theorem 1, let $\pi_{\beta}$ be the invariant measure and
$\bar{\theta}_{\beta}:=\mathbb{E}_{\theta\sim\pi_{\beta}}[\theta]$ be its
expectation. For $\theta\sim\pi_{\beta}$, the invariant distribution
$\pi_{\beta}$ satisfies the following properties.

  1. (a)

$\mathbb{E}_{\theta\sim\pi_{\beta}}\big{[}\nabla\mathcal{L}(\theta)\big{]}=0.$
In particular, if the gradient $\nabla\mathcal{L}$ is linear (see Assumption 6
below) then we have $\bar{\theta}_{\beta}=\theta^{\star}$.

  2. (b)

The variance and the bias are bounded as follows:

$\operatorname{Var}_{\pi_{\beta}}(\theta)\leq\frac{\beta\sigma^{2}}{2\mu-\beta(\mu^{2}+L_{\sigma})}\quad\text{
and
}\quad\big{\|}\bar{\theta}_{\beta}-\theta^{\star}\big{\|}\leq\sqrt{\operatorname{Var}_{\pi_{\beta}}(\theta)}.$

Proposition 1 is proven in Section B.2 and relies on the invariance property
of $\pi_{\beta}$ and the contraction property of the optimisation iteration.
The second part of this statement reflects the known-fact that the iterates
have an asymptotic magnitude of $\sqrt{\beta}$ [68, 63].

Before stating further results, we need to define sub-Gaussian and sub-
exponential concentration properties for real random variables. Among the many
known equivalent characterizations, we only introduce those required for the
proofs of our results, see [89, Chapter 2] for other characterizations.

######  Definition 1.

Let $X$ be a real random variable. We say that $X$ is $K$-sub-Gaussian for
some $K>0$ whenever

  1. (i)

we have

$\mathbb{E}\exp(\lambda^{2}X^{2})\leq e^{\lambda^{2}K^{2}}\quad\text{for}\quad
0\leq\lambda\leq 1/K,$ (8)

which we will denote $X\in\widetilde{\Psi}_{2}(K),$

  2. (ii)

or we have

$\mathbb{E}\exp(\lambda X)\leq\exp(\lambda^{2}K^{2})\quad\text{ for
all}\quad\lambda\in\mathbb{R},$ (9)

which we will denote $X\in\Psi_{2}(K)$.

######  Definition 2.

Let $X$ be a real random variable. We say that $X$ is sub-exponential if one
of the two following conditions holds.

  1. (i)

There exists $K_{1}>0$ such that

$\|X\|_{L_{p}}\leq K_{1}p\quad\text{for all}\quad p\geq 1,$ (10)

in which case we write $X\in\widetilde{\Psi}_{1}(K_{1}).$

  2. (ii)

There exists $K_{2}$ such that

$\mathbb{E}\exp(\lambda X)\leq\exp(\lambda^{2}K_{2}^{2})\quad\text{ for
all}\quad|\lambda|\leq 1/K_{2},$ (11)

in which case we write $X\in\Psi_{1}(K_{2}).$

Note that, for a centered variable $X,$ the two characterizations of
Definition 1 are equivalent with the same constant $K.$ Analogously, for
centered $X,$ we have that $X\in\Psi_{1}(K)$ entails
$X\in\widetilde{\Psi}_{1}(2eK)$ and $X\in\widetilde{\Psi}_{1}(K)$ entails
$X\in\Psi_{1}(2eK).$ Namely, the two characterizations of Definition 2 imply
each other but with worse constants (see [89, Section 2] for a proof). Since
the constants in Definition 2 degrade by switching between the two properties,
we will specify which property is meant in each subsequent statement in order
to minimize these degradations.

We first formulate a sub-Gaussian/sub-exponential concentration assumption on
the norms of the gradient errors.

######  Assumption 3.

There exists $\overline{K}<+\infty$ such that one of the following holds:

  1. (a)

For all $\theta\in\Theta,$ the gradient error satisfies
$\|\varepsilon_{\zeta}(\theta)\|\in\widetilde{\Psi}_{2}(\overline{K}).$

  2. (b)

For all $\theta\in\Theta,$ the gradient error satisfies
$\|\varepsilon_{\zeta}(\theta)\|\in\widetilde{\Psi}_{1}(\overline{K}).$

In combination with Assumptions 1 and 2, the above pair of assumptions imply
the following concentration properties for $\pi_{\beta}.$

######  Proposition 2.

In the setting of Theorem 1, for $\theta\sim\pi_{\beta},$ the invariant
distribution $\pi_{\beta}$ satisfies the following properties:

  1. (a)

If Assumption 3 (a) holds then
$\|\theta-\theta^{\star}\|\in\widetilde{\Psi}_{2}\big{(}\overline{K}\sqrt{8\beta/\mu}\big{)}.$

  2. (b)

If Assumption 3 (b) holds and if $\beta\leq(2\mu)^{-1}$ then
$\|\theta-\theta^{\star}\|\in\widetilde{\Psi}_{1}\big{(}2\overline{K}\sqrt{\beta/\mu}\big{)}.$

The proof of Proposition 2 is given in Section B.3. The most important aspect
of this statement is that the quantity $\|\theta-\theta^{\star}\|$ is sub-
Gaussian/sub-exponential with a constant _depending_ on the step-size $\beta.$
Indeed, it is fairly easy to show, for example, that
$\|\theta-\theta^{\star}\|\in\widetilde{\Psi}_{2}(\overline{K}/\mu)$ under
Assumption 3 (a), however, this constant is too pessimistic since it fails to
take advantage of a small step-size which leads to stronger concentration. The
improved constants above are obtained by carefully leveraging the unbiased
property of the gradient error (see Assumption 2). Previous characterizations
of $\pi_{\beta}$ obtained bounds on the bias w.r.t. $\theta^{\star}$ [91] and
moment expansions of $\widehat{\theta}-\theta^{\star}$ for
$\widehat{\theta}\sim\pi_{\beta}$ or $\widehat{\theta}$ equal to a Polyak-
Ruppert average [18], however, the sub-exponential and sub-Gaussian
characterizations of Proposition 2 appear to be new.

Note that the constant $\overline{K}$ obtained from Assumption 3 and appearing
in Proposition 2 may hide a dependence on the dimension in $\sqrt{d}$ since it
is related to the Euclidean norm $\|\varepsilon_{\zeta}(\theta)\|$ of the
gradient noise. In this respect, Proposition 2 resembles the results of [56]
where similar hypotheses to Assumption 3 were used entailing the same
dimension dependence. In order to avoid this shortcoming, one needs a stronger
assumption which is stated along with the associated results further below.
Note also that Assumption 3 considerably strengthens Assumption 2 by requiring
that $\|\varepsilon_{\zeta}(\theta)\|$ admits a finite exponential moment. In
addition, the involved bound is uniform w.r.t. $\theta.$ However, under a non-
uniform finite $p$-moment assumption, it is still possible to show the
following.

######  Lemma 1.

Grant Assumptions 1 and 2 and assume that there is $K,\underline{K}>0$ and
$p\in\mathbb{N}^{*}$ such that conditionally on any $\theta$ we have

$\big{\|}\|\varepsilon_{\zeta}(\theta)\|\big{\|}_{L_{p}}\leq
K\|\theta-\theta^{\star}\|+\underline{K},$ (12)

then for step-size $\beta$ as in Theorem 1 and satisfying the additional
condition $\beta\leq\frac{\mu}{j(\mu^{2}+K^{2})}$ with $j\leq p,$ the Markov
chain $(\theta_{t})_{t\geq 0}$ converges to an invariant distribution
$\pi_{\beta}$ with at least $j$ finite moments.

Lemma 1 is proved in Section B.4 and shows that $\pi_{\beta}$ can have as many
finite moments as the gradient, provided that the step-size is small enough.
This implies that even weaker concentration properties transfer to the
invariant distribution. Note that a non-uniform sub-exponential (resp. sub-
Gaussian) assumption would correspond to condition (12) with $K$ replaced by
$Kp$ (resp. $K\sqrt{p}$) in which case the condition on $\beta$ becomes at
least $\beta\leq\mathcal{O}\big{(}\mu/(jK^{2}p)\big{)}$. This suggests that,
for arbitrary $p$, $\pi_{\beta}\in L_{p}$ may only hold in the limit $\beta\to
0$. Lemma 1 shares a similarity with [2, Theorem 2] which states an $L_{p}$
convergence result under a resembling condition on the step-size. However, [2,
Theorem 2] uses a uniform $L_{p}$ condition on the gradient error whereas we
allow the upperbound to depend on $\|\theta-\theta^{\star}\|$ in Inequality
(12).

We now introduce a stronger analog to Assumption 3, which will enable later
the proof of dimension-free deviation bounds..

######  Assumption 4.

There is $K<+\infty$ such that one of the following holds:

  1. (a)

For all $\theta\in\Theta$ and all $f\in\operatorname{Lip}(\mathbb{R}^{d}),$ we
have $f(G(\theta,\zeta))-\mathbb{E}f(G(\theta,\zeta))\in\Psi_{2}(K).$

  2. (b)

For all $\theta\in\Theta$ and all $f\in\operatorname{Lip}(\mathbb{R}^{d}),$ we
have $f(G(\theta,\zeta))-\mathbb{E}f(G(\theta,\zeta))\in\Psi_{1}(K).$

As announced, the subtle difference with Assumption 3 is that the involved
constants are, a priori, independent from the dimension. The so-called Bobkov-
Götze theorem [7] states that Assumption 4 (a) is equivalent to the fact that
$\nu_{\theta}:=\mathcal{D}(G(\theta,\zeta))$ satisfies the following
_Transportation-Information_ inequality

$\mathcal{W}_{1}\big{(}\nu,\nu_{\theta}\big{)}\leq\sqrt{2K^{2}D(\nu\|\nu_{\theta})}\quad\text{
for all }\quad\nu\in\mathcal{M}_{1}(\mathbb{R}^{d}),$ (13)

where $\mathcal{W}_{1}$ and $D(\cdot\|\cdot)$ are the Wasserstein-$1$ distance
[90] (see definition below) and the Kullback-Leibler divergence [44] between
probability measures respectively. An analogous equivalence may be established
for the sub-exponential case of Assumption 4 (b) (for instance, by adapting
the proof given in [86, Theorem 4.8]).

By restricting the functions $f$ in Assumption 4 to be linear, we recover the
assumption that the vector $G(\theta,\zeta)$ is sub-Gaussian/sub-exponential.
An interesting question is then whether this weaker property implies
Assumption 4 with a dimension independent constant. To the best of our
knowledge of the current literature, this is only known to hold for Gaussian
vectors (see for instance [86, Theorem 3.25]). In fact, Talagrand’s well-known
transport inequality states that Gaussian vectors satisfy Inequality (13) for
the $\mathcal{W}_{2}$ distance rather than $\mathcal{W}_{1}$, which is an even
stronger property. Since Inequality (13) involves two very different forms of
distance between probability measures, a direct intuitive understanding of its
meaning is elusive. However, the above inequality is related to a host of
properties used to describe the concentration of measure phenomenon including
Poincaré inequalities [8, 28], logarithmic Sobolev inequalities [7, 47] and
modified logarithmic Sobolev inequalities [25, 3] to mention only a few. A
broad and comprehensive survey on transport inequalities and their
consequences on concentration and deviation inequalities is available in [29].

Using the previous assumption, we can show that the invariant distribution
inherits similar properties.

######  Proposition 3.

Grant Assumptions 1 and 2 and let $\pi_{\beta}$ be the invariant distribution
obtained for step-size $\beta$ as in Theorem 1. For $\theta\sim\pi_{\beta}$,
the invariant distribution $\pi_{\beta}$ satisfies the following properties.

  1. (a)

If Assumption 4 (a) holds then
$f(\theta)-\mathbb{E}f(\theta)\in\Psi_{2}\big{(}K\sqrt{\beta/\mu}\big{)}$ for
all $f\in\operatorname{Lip}(\mathbb{R}^{d}).$

  2. (b)

If Assumption 4 (b) holds then
$f(\theta)-\mathbb{E}f(\theta)\in\Psi_{1}\big{(}K\sqrt{\beta/\mu}\big{)}$ for
all $f\in\operatorname{Lip}(\mathbb{R}^{d}).$

Proposition 3 is proven in Section B.5 and will be used in Section 6 to derive
dimension-free deviation bounds. Note that the $\Psi_{1}$/$\Psi_{2}$ constants
in Proposition 3 also display the crucial $\sqrt{\beta/\mu}$ dependence as in
Proposition 2 and without further degradation. Before proceeding to the
statement of such results we explore another convergence mode of the SGD
Markov chain.

##  5 Wasserstein Convergence

This section complements Theorem 1 with an additional convergence result
w.r.t. the Wasserstein metric. We recall that, for $p\geq 1$ and two
distributions $\varpi,\nu\in\mathcal{M}_{1}(\mathbb{R}^{d}),$ the
Wasserstein-$p$ distance is defined by

$\mathcal{W}_{p}^{p}(\varpi,\nu)=\inf_{\xi\in\Pi(\varpi,\nu)}\mathbb{E}_{X,Y\sim\xi}\|X-Y\|^{p},$

where $\Pi(\varpi,\nu)$ is the set of all couplings of $\varpi$ and $\nu$ i.e.
distributions over $\mathbb{R}^{d}\times\mathbb{R}^{d}$ with first and second
marginals equal to $\varpi$ and $\nu$ respectively.

In order to show that the SGD iteration converges w.r.t. the Wasserstein-$2$
distance, we require the following assumption.

######  Assumption 5.

There is $L_{\mathcal{W}}<+\infty$ such that for all $\theta,\theta^{\prime},$
the gradient noise distributions $\mathcal{D}(\varepsilon_{\zeta}(\theta))$
and $\mathcal{D}(\varepsilon_{\zeta}(\theta^{\prime}))$ at $\theta$ and
$\theta^{\prime}$ satisfy

$\mathcal{W}_{2}^{2}\big{(}\mathcal{D}(\varepsilon_{\zeta}(\theta)),\mathcal{D}(\varepsilon_{\zeta}(\theta^{\prime}))\big{)}\leq
L_{\mathcal{W}}\|\theta-\theta^{\prime}\|^{2}.$

In words, we assume that the change in the gradient noise distribution
measured with the $\mathcal{W}_{2}$ metric is controlled by the change in the
parameter $\theta.$ This assumption is discussed below and allows to obtain
the following result.

######  Proposition 4.

Grant Assumptions 1, 2 and 5. Let
$\nu_{1},\nu_{2}\in\mathcal{P}_{2}(\mathbb{R}^{d})$ be two initial
distributions and let $\beta$ be a step-size such that

$\beta<\frac{2\mu}{\mu^{2}+(L_{\mathcal{W}}\vee\mu L)},$

then we have the contraction

$\mathcal{W}_{2}^{2}(\nu_{1}P_{\beta},\nu_{2}P_{\beta})\leq\big{(}(1-\beta\mu)^{2}+\beta^{2}L_{\mathcal{W}}\big{)}\mathcal{W}_{2}^{2}(\nu_{1},\nu_{2}).$

Consequently, for such a $\beta$ and any initial
$\theta_{0}\sim\nu\in\mathcal{P}_{2}(\mathbb{R}^{d}),$ the Markov chain
generated by iteration (4) converges to a unique stationary measure
$\pi_{\beta}$ in $\mathcal{W}_{2}$ distance.

The proof of Proposition 4 is given in Section B.6. The intuition behind it is
that, if the Markov chain evolves according to a locally similar dynamic when
started from different points then, for small enough step-size, the
contraction phenomenon coming from the optimization will prevail so that
trajectories associated to different initializations join even before
convergence. A similar result was previously stated in [18, Proposition 2] for
smooth and strongly convex functions also. The main difference with
Proposition 4 is that we substitute Assumption 5 for Assumption A4 from [18]
which requires that every random gradient $G(\theta,\zeta)$ be almost surely
co-coercive with fixed constant. We propose to compare these assumptions for
the simple example of least-squares linear regression in which, given a sample
$\zeta=(X,\xi)\in\mathbb{R}^{d}\times\mathbb{R}$, a random gradient is
computed as:

$G(\theta,\zeta)=XX^{\top}\theta-XY\quad\text{with}\quad
Y=X^{\top}\theta^{\star}+\xi,$

where $\xi$ is an independent centered noise. The co-coercivity assumption
requires that, for all $\theta,\theta^{\prime},$ the following inequality
holds almost surely

$L\langle
G(\theta,\zeta)-G(\theta^{\prime},\zeta),\theta-\theta^{\prime}\rangle\geq\|G(\theta,\zeta)-G(\theta^{\prime},\zeta)\|^{2}.$

In this particular example, this implies that $\|X\|^{2}\leq L$ almost surely,
which is quite restrictive. On the other hand we have
$\varepsilon_{\zeta}(\theta)=G(\theta,\zeta)-\nabla\mathcal{L}(\theta)=(XX^{\top}-\mathbb{E}XX^{\top})(\theta-\theta^{\star})-X\xi$
and it is easy to couple the distributions of
$\varepsilon_{\zeta}(\theta^{\prime})$ and
$\varepsilon_{\zeta^{\prime}}(\theta^{\prime})$ by defining them with the same
variables $\zeta=\zeta^{\prime}=(X,\xi)$ so that we find

$\mathcal{W}_{2}^{2}\big{(}\mathcal{D}(\varepsilon_{\zeta}(\theta)),\mathcal{D}(\varepsilon_{\zeta^{\prime}}(\theta^{\prime}))\big{)}\leq\mathbb{E}\|\varepsilon_{\zeta}(\theta)-\varepsilon_{\zeta^{\prime}}(\theta^{\prime})\|^{2}\leq\mathbb{E}\|XX^{\top}-\mathbb{E}XX^{\top}\|_{2}^{2}\|\theta-\theta^{\prime}\|^{2},$

where $\|\cdot\|_{2}$ is the operator norm. Assumption 5 is then satisfied as
soon as $X$ has a bounded fourth moment which is a significantly more lenient
condition since it is even weaker than a finite variance assumption on the
gradient in this particular case.

More generally, one may consider an objective $\mathcal{L}$ defined by a
linear learning task such that
$\mathcal{L}(\theta)=\mathbb{E}_{(X,Y)}[\ell(X^{\top}\theta,Y)]$ for a convex
loss $\ell$ and random samples and labels $(X,Y)$ so that the gradient samples
are $G(\theta,(X,Y))=X\ell^{\prime}(X^{\top}\theta,Y)$ with $\ell^{\prime}$
the derivative in the first argument. By a similar computation as above, it is
possible to verify Assumption 5 as soon as $\ell$ is smooth in its first
argument and $X$ has a finite fourth moment.

In the same vein as Assumption 5, it is possible to introduce a regularity
condition on the transition kernel $P_{\beta}$ in terms of the
$\operatorname{TV}$ distance which allows to obtain the following result.

######  Proposition 5.

Let the assumptions of Proposition 4 hold and further assume that$:$

  * •

For all $\theta\in\mathbb{R}^{d}$ the probability measure
$P_{\beta}(\theta,\cdot)$ admits a density $p_{\beta}(\theta,\omega)$ w.r.t.
Lebesgue’s measure.

  * •

There exists $A<\infty$ such that for all
$\theta,\theta^{\prime}\in\mathbb{R}^{d}$

$\|P_{\beta}(\theta,\cdot)-P_{\beta}(\theta^{\prime},\cdot)\|_{\operatorname{TV}}=\frac{1}{2}\int_{\mathbb{R}^{d}}|p_{\beta}(\theta,\omega)-p_{\beta}(\theta^{\prime},\omega)|d\omega\leq
A\|\theta-\theta^{\prime}\|.$ (14)

Then the ergodicity result of Theorem 1 holds with
$\rho\leq\sqrt{(1-\beta\mu)^{2}+\beta^{2}L_{\mathcal{W}}}.$

###### Proof.

Using [55, Theorem 12] (see also [55, Lemma 13]), the assumptions imply that
for all $\varpi,\nu\in\mathcal{M}_{1}(\mathbb{R}^{d})$ we have:

$\|\varpi P_{\beta}-\nu P_{\beta}\|_{\operatorname{TV}}\leq
A\mathcal{W}_{1}(\varpi,\nu).$

It then only remains to use Proposition 4 with
$\varpi=\delta_{\theta_{0}}P_{\beta}^{n-1}$ and
$\nu=\pi_{\beta}=\pi_{\beta}P_{\beta}$ along with the inequality:

$\mathcal{W}_{1}(\varpi,\nu)\leq\sqrt{\mathcal{W}_{2}^{2}(\varpi,\nu)}$

valid for all $\varpi,\nu\in\mathcal{M}_{1}(\mathbb{R}^{d}).$ ∎

Proposition 5 uses the “Wasserstein-to-TV” method [71, 55] in order to derive
convergence in $\operatorname{TV}$ distance from Proposition 4 which leads to
an explicit estimate of $\rho$ in Theorem 1. Note however that although the
necessary condition (14) is quite intuitive, its verification is not
straightforward even for a toy example.

##  6 Confidence bounds

Using the convergence and concentration results formulated in the previous
sections for the invariant distribution of the SGD Markov chain, we are ready
to state confidence bounds on the estimation of the optimal $\theta^{\star}.$
Recall that by Proposition 1, the invariant distribution $\pi_{\beta}$ may not
be centered around $\theta^{\star}$ unless the gradient is linear, which is a
particular case. In general, the expectation of $\pi_{\beta}$ may not be equal
to $\theta^{\star}$ but the bias is controlled by the step-size $\beta.$
Therefore, two possibilities are available for the final estimator:

  * •

The last iterate $\theta_{T}$: with $T$ the optimization horizon. In which
case a small step size is appropriate.

  * •

A tail average $\frac{1}{n}\sum_{j=n_{0}+1}^{n_{0}+n}\theta_{j}$: in which
case the step size may be chosen reasonably large within the convergence
conditions.

###  6.1 Final iterate concentration bounds

When the expectation of the invariant measure $\overline{\theta}_{\beta}$
differs from the true optimum $\theta^{\star}$, one may choose a small step-
size $\beta$ in order to obtain a precise estimator of $\theta^{\star}$
through the final iterate $\theta_{T}$ of iteration (4). In the event where
Assumption 3 holds, a quick application of Proposition 2 combined with Theorem
1 leads to the following first deviation bounds.

######  Corollary 1.

In the setting of Proposition 2, for a step-size $\beta$ and horizon $T,$ we
have the following high-confidence bounds:

  1. (a)

Under Assumption 3 (a), for $\delta>0,$ we have

$\mathbb{P}\Big{(}\big{\|}\theta_{T}-\theta^{\star}\big{\|}>\overline{K}\sqrt{8\beta\log(e/\delta)/\mu}\Big{)}\leq\delta+\rho^{T}M\big{(}1+\|\theta_{0}-\theta^{\star}\|^{2}\big{)}.$

  2. (b)

Under Assumption 3 (b), for $\delta>0,$ we have

$\mathbb{P}\Big{(}\big{\|}\theta_{T}-\theta^{\star}\big{\|}>4e\overline{K}\log(2/\delta){\sqrt{\beta/\mu}}\Big{)}\leq\delta+\rho^{T}M\big{(}1+\|\theta_{0}-\theta^{\star}\|^{2}\big{)},$

where $\rho$ and $M$ are the constants obtained from Theorem 1.

The proof of Corollary 1 is given in Section B.7. As discussed earlier, the
constants $\overline{K}$ drawn from Assumption 3 may have a poor dependence on
the dimension in $\sqrt{d}$ which leaves room for improvement in the above
bounds. This can be achieved when the requirements of Assumption 4 are met
leading to the following _dimension-free_ deviation bounds.

######  Corollary 2.

In the setting of Proposition 3, for a step-size $\beta$ such that
$\beta\leq\frac{\mu}{\mu^{2}+L_{\sigma}},$ we have the following high-
confidence bounds:

  1. (a)

Under Assumption 4 (a), for $\delta>0,$ we have

$\mathbb{P}\bigg{(}\big{\|}\theta_{n}-\theta^{\star}\big{\|}>\sqrt{\frac{\beta\sigma^{2}}{\mu}}+2K\sqrt{\frac{\beta\log(1/\delta)}{\mu}}\bigg{)}\leq\delta+\rho^{n}M\big{(}1+\|\theta_{0}-\theta^{\star}\|^{2}\big{)}.$
(15)

  2. (b)

Under Assumption 4 (b), for $\delta>0,$ we have

$\mathbb{P}\bigg{(}\big{\|}\theta_{n}-\theta^{\star}\big{\|}>\sqrt{\frac{\beta\sigma^{2}}{\mu}}+2K\bigg{(}\sqrt{\frac{\beta\log(1/\delta)}{\mu}}\>\>\vee\>\>\beta\log(1/\delta)\bigg{)}\bigg{)}\leq\delta+\rho^{n}M\big{(}1+\|\theta_{0}-\theta^{\star}\|^{2}\big{)},$
(16)

where $\rho$ and $M$ are the constants obtained from Theorem 1.

Corollary 2 is proven in Section B.8 and is a consequence of Propositions 1, 3
and Theorem 1. Inequality (15) is an example of a sub-Gaussian deviation bound
since the uncertainty term has no dependence on the dimension [50].

The bounds given in Corollaries 1 and 2 can be made more explicit by granting
the conditions of Proposition 5 which yields $\rho\approx 1-\beta\mu.$ In this
case, it is possible to set $\beta=\frac{\log(T)}{\mu T}$ leading to a nearly
optimal statistical rate in $T$ albeit with a sub-optimal dependence in the
problem conditioning [1, 2].

###  6.2 Polyak-Ruppert averaging

In this part, we consider the case where the step-size $\beta$ is chosen as a
constant order value satisfying the convergence criteria required in our
previous results. Our goal is to obtain a high-confidence bound for the
Polyak-Ruppert average $\frac{1}{n}\sum_{j=n_{0}+1}^{n_{0}+n}\theta_{j}$
computed after a burn-in period of $n_{0}$ iterations. This raises two
challenges, the first of which is that, even for a very long burn-in period
$n_{0},$ the stationary regime is never reached in theory so that one cannot
immediately use the concentration properties of $\pi_{\beta}.$ The second
challenge comes from the lack of independence of the Markov chain iterates.
This prevents the adoption of certain approaches such as the entropy method as
done in [56] for example.

Notice that, unless the gradient is linear, there is little hope to estimate
$\theta^{\star}$ using the Polyak-Ruppert average since it is bound to
approach
$\mathbb{E}_{\theta\sim\pi_{\beta}}[\theta]=\overline{\theta}_{\beta}$ which
may differ from $\theta^{\star}$ by up to $\sigma\sqrt{\beta/\mu}$ in the non
linear case. Nevertheless, the following initial statement holds without this
assumption.

######  Theorem 2.

Grant Assumptions 1, 2, 4 (a) and 5. Let $f:\Theta^{n}\to\mathbb{R}$ be a
$1$-Lipschitz function in each of its parameters and
$\vec{\theta}:=(\theta_{1},\dots,\theta_{n})$ be a sequence of SGD iterates
with step size $\beta<\frac{2\mu}{\mu^{2}+(\mu L\vee L_{\mathcal{W}})}$
started from stationarity i.e. such that $\theta_{1}\sim\pi_{\beta}$. Then we
have

$f(\vec{\theta}\>)-\mathbb{E}f(\vec{\theta})\in\Psi_{2}\big{(}KC_{\mathcal{W}}\sqrt{\beta/\mu+(n-1)\beta^{2}}\big{)},$

where
$C_{\mathcal{W}}=\big{(}1-\sqrt{(1-\beta\mu)^{2}+\beta^{2}L_{\mathcal{W}}}\big{)}^{-1}$.
If Assumption 4 (a) is replaced by Assumption 4 (b) then

$f(\vec{\theta}\>)-\mathbb{E}f(\vec{\theta})\in\Psi_{1}\big{(}KC_{\mathcal{W}}\sqrt{\beta/\mu+(n-1)\beta^{2}}\big{)}.$

The proof of Theorem 2 is given in Section B.9 and employs a hybrid martingale
transportation method (see [10, 57, 21, 12] for a reference) leveraging the
$\mathcal{W}_{2}$ convergence established in Proposition 4 in combination with
[42, Theorem 4.3].

Theorem 2 may be used in a variety of ways by plugging different choices of
the function $f.$ For instance, one may choose
$f(\vec{\theta})=\sum_{i}g(\theta_{i})$ for any
$g\in\operatorname{Lip}(\mathbb{R}^{d}).$ Instead, in what follows, we focus
on the choice

$f(\vec{\theta})=\Big{\|}\sum_{i=n_{0}+1}^{n_{0}+n}\theta_{i}-n\theta^{\star}\Big{\|}.$

Before we proceed, we formalize the gradient linearity assumption.

######  Assumption 6.

The gradient $\nabla\mathcal{L}$ is linear i.e. for all $\theta\in\Theta$ it
is equal to $\nabla\mathcal{L}(\theta)=\Sigma(\theta-\theta^{\star})$ for some
symmetric positive definite matrix $\Sigma\in\mathbb{R}^{d\times d}.$

Note that the positive definiteness of $\Sigma$ in Assumption 6 is a
consequence of strong convexity while its symmetry is a result of the Hessian
$\nabla^{2}\mathcal{L}$ being constant in this case and therefore continuous.
We are now ready to state our non-asymptotic deviation bound for Polyak-
Ruppert averaging.

######  Proposition 6.

Grant Assumptions 1, 2, 4 (a), 5 and 6. Let $(\theta_{t})_{t\geq 0}$ be the
Markov sequence obtained by running SGD with step-size

$\beta<\frac{2\mu}{\mu^{2}+(\mu L\vee
L_{\mathcal{W}})}\wedge\frac{\mu}{\mu^{2}+L_{\sigma}}$

and initial distribution $\theta_{0}\sim\nu.$ Then there exist $\rho<1$ and
$M<\infty$ such that

$\begin{split}\Big{\|}\frac{1}{n}\sum_{t=n_{0}+1}^{n_{0}+n}\theta_{t}-\theta^{\star}\Big{\|}&\leq\sqrt{\frac{2}{n}\frac{1+\alpha}{1-\alpha}\Big{(}\alpha_{\mathcal{W}}^{n_{0}}\mathcal{W}_{2}^{2}(\nu,\pi_{\beta})+\frac{\beta\sigma^{2}}{\mu}\Big{)}}\\\
&\quad+\frac{2K\sqrt{\beta/\mu}}{1-\alpha_{\mathcal{W}}}\sqrt{\beta\mu+\frac{1}{n}}\sqrt{\frac{\log(1/\delta)}{n}}\end{split}$
(17)

for $\delta>0$ and $n,n_{0}>0$ with probability at least
$1-\Upsilon(\nu,n_{0})\delta,$ where

$\displaystyle\alpha=1-\beta\mu,\quad\alpha_{\mathcal{W}}=\sqrt{\alpha^{2}+\beta^{2}L_{\mathcal{W}}}\quad\text{and}\quad\Upsilon(\nu,n_{0})=1+M\rho^{n_{0}}\Big{\|}\frac{d\nu}{d\pi_{\beta}}\Big{\|}_{\infty}.$

If Assumption 4 (a) is replaced by Assumption 4 (b) then

$\begin{split}\Big{\|}\frac{1}{n}\sum_{t=n_{0}+1}^{n_{0}+n}\theta_{t}-\theta^{\star}\Big{\|}&\leq\sqrt{\frac{2}{n}\frac{1+\alpha}{1-\alpha}\Big{(}\alpha_{\mathcal{W}}^{n_{0}}\mathcal{W}_{2}^{2}(\nu,\pi_{\beta})+\frac{\beta\sigma^{2}}{\mu}\Big{)}}\\\
&\quad+\frac{2K\sqrt{\beta/\mu}}{1-\alpha_{\mathcal{W}}}\bigg{(}\sqrt{\beta\mu+\frac{1}{n}}\sqrt{\frac{\log(1/\delta)}{n}}\vee\frac{\log(1/\delta)}{n}\bigg{)}\end{split}$
(18)

holds with the same probability.

The proof of Proposition 6 is given in Section B.10 and takes advantage of the
convergence both in total-variation distance and in the $\mathcal{W}_{2}$
metric. Note that the given bounds are also dimension-free thanks to
Assumption 4. It is possible to derive a weaker result using only Assumption 3
but we omit it to avoid repetition. The variance terms in the upperbounds of
(17) and (18) (those independent of $\delta$) are controlled thanks to a
geometric decorrelation phenomenon which can be shown for the Markov chain
iterates under Assumption 6 (see Lemma 6 in the Appendix). This phenomenon
becomes weaker for smaller step-size $\beta,$ therefore, it only makes sense
to apply Proposition 6 with $\beta$ of constant order to avoid excessive
correlation between the averaged samples. Finally, the lack of stationarity of
the involved Markov samples is tackled by taking advantage of a spectral gap
property satisfied by the transition kernel $P_{\beta}$ under the conditions
of Theorem 1 (see [43]).

Proposition 6 may be compared to the works of [60] and [49]. The former
derives a similar high probability bound for linear stochastic approximation
under a generalized sub-Gaussianity assumption and uncorrelated noise. The
latter considers a weaker finite $L_{p}$ moment assumption on the SGD data and
uses mini-batching to obtain Nagaev type concentration bounds with provably
optimal dependence in the confidence level. However, the results of [60, 49]
both lack the dimension-free property of Proposition 6.

##  7 Applications

We discuss the consequences of our results for two common use-cases of SGD.

###  7.1 Linear regression

Linear regression is one of the most popular and most used standard models.
The aim is to predict a real variable $Y$ based on a random vector
$X\in\mathbb{R}^{d}$ according to the linear model

$Y=X^{\top}\theta^{\star}+\epsilon$

where $\theta^{\star}$ is an unknown parameter and $\epsilon$ a centered
noise. The estimation of $\theta^{\star}$ may be carried out by minimizing the
least-squares objective
$\mathcal{L}(\theta):=\frac{1}{2}\mathbb{E}\big{(}X^{\top}\theta-Y\big{)}^{2}$
with respect to $\theta\in\mathbb{R}^{d}.$ This may be done by running SGD
with the random gradient $G(\theta,(X,Y))=X(X^{\top}\theta-Y).$

Provided the previous gradient admits a finite second moment, Theorem 1 and
Proposition 4 apply and guarantee the convergence of the SGD Markov chain in
total-variation and $\mathcal{W}_{2}$ distance. If the covariates $X$ and the
noise $\epsilon$ are both Gaussian then the gradient $G(\theta,(X,Y))$ is sub-
exponential. However, note that Assumption 3 or (b) are not immediately
satisfied since the associated $\Psi_{1}$ constant may be unbounded for
arbitrarily high values of $\|\theta-\theta^{\star}\|.$ This problem can be
remedied thanks to the following lemma.

######  Lemma 2.

Let Assumption 1 hold and assume that gradient errors write
$\varepsilon_{\zeta_{t}}(\theta_{t})=\Xi_{t}(\theta_{t}-\theta^{\star})+\xi_{t}$
where the pairs $(\Xi_{t},\xi_{t})_{t\geq 0}$ are i.i.d in
$\mathbb{R}^{d\times d}\times\mathbb{R}^{d}$ with $\Xi_{t}$ symmetric and such
that for all $u\in\mathbb{R}^{d},\|u\|=1$ we have $\langle
u,\Xi_{t}u\rangle\in\Psi_{1}(K_{\Xi})$ and $\langle
u,\xi_{t}\rangle\in\Psi_{1}(K_{\xi})$ for $K_{\Xi},K_{\xi}>0.$ Assume the
following minibatch SGD iteration is run starting from $\theta_{0}$ such that
$\|\theta_{0}-\theta^{\star}\|\leq C$ for some $C>0$ for a finite horizon $T$

$\theta_{t+1}=\theta_{t}-\beta\overline{G}_{N}(\theta_{t})\quad\text{with}\quad\overline{G}_{N}(\theta_{t})=\frac{1}{N}\sum_{i=1}^{N}G(\theta_{t},\zeta_{tN+i})$

with $N$ the minibatch size. For a confidence level $\delta>0,$ assume that
$N$ and $\beta$ satisfy

$\frac{N}{\log(4T/\delta)\\!+\\!3d}\geq
1\vee\Big{(}\frac{6}{\mu}\big{(}3K_{\Xi}\vee
4K_{\xi}/C\big{)}\Big{)}^{2}\quad\text{and}\quad\beta\leq\frac{\mu
N}{54K_{\Xi}^{2}(\log(4T/\delta)\\!+\\!3d)}\wedge\frac{2}{\mu\\!+\\!L}.$

Then, with probability at least $1-\delta,$ we have $\max_{0\leq s\leq
T}\|\theta_{s}-\theta^{\star}\|\leq C$.

Lemma 2 is proven in Section B.11 and guarantees that, using a small step size
and minibatching to reduce the gradient variance, with high probability, the
iteration does not stray from the vicinity of the optimum during a finite
horizon. This shows that the uniform aspect of Assumptions 3 and 4 does not
prevent the application of the results given in the previous sections. Note
that although Lemma 2 requires that $N=\Omega(d),$ the constant $C$ is
arbitrary and may be taken dimension-free, for instance, by starting the
iteration from a preliminary estimator $\theta_{0}=\widehat{\theta}.$

For the example of linear regression with sub-Gaussian samples
$(X_{t},Y_{t})_{t},$ Lemma 2 applies with
$\Xi_{t}=X_{t}X_{t}^{\top}-\mathbb{E}[X_{t}X_{t}^{\top}]$ and
$\xi_{t}=-\epsilon_{t}X_{t}.$ Thus, for finite horizon, one may consider the
event where the bound of Lemma 2 holds to apply results from Sections 4 and 6.

Alternatively, it is also possible to restrict the optimization to a convex
and bounded subset $\Theta\subset\mathbb{R}^{d}$ such that
$\theta^{\star}\in\Theta.$ By letting $\Pi_{\Theta}(\cdot)$ be the projection
onto $\Theta$ and replacing iteration (4) with

$\theta_{t+1}=\Pi_{\Theta}\big{(}\theta_{t}-\beta
G(\theta_{t},\zeta_{t})\big{)},$ (19)

we obtain a Markov chain to which Proposition 2 (b) applies and leads to the
deviation bound (16). Indeed, it is easy to verify that these results still
hold for iteration (19) thanks to the inequality

$\big{\|}\Pi_{\Theta}\big{(}\theta-\beta
G(\theta,\zeta)\big{)}-\theta^{\star}\big{\|}\leq\big{\|}\theta-\beta
G(\theta,\zeta)-\theta^{\star}\big{\|},$

valid for all $\theta\in\mathbb{R}^{d}$ since $\theta^{\star}\in\Theta$ which
is convex. However, by considering the projected iteration (19), Proposition 1
(a) may no longer hold so that $\overline{\theta}_{\beta}\neq\theta^{\star}$
making Proposition 6 no longer applicable.

###  7.2 Logistic regression

Logistic regression corresponds to the model

$1-\mathbb{P}(Y=-1|X)=\mathbb{P}(Y=+1|X)=\sigma(X^{\top}\theta^{\star}),$

where $\sigma$ is the sigmoid function $\sigma(x)=1/(1+e^{-x}).$ For a
parameter $\theta$ and a sample $X,$ the predicted probability is
$\mathbb{P}(Y=+1|X)=\sigma(X^{\top}\theta)$ and the model is trained using the
log-loss $\ell(\theta,(X,Y))=-\log(\sigma(YX^{\top}\theta))$ which yields the
objective $\mathcal{L}(\theta)=\mathbb{E}_{(X,Y)}\ell(\theta,(X,Y)).$

In order to ensure the objective is strongly-convex, it is necessary to
restrict the parameter $\theta$ to a bounded convex set $\Theta.$ This is
commonly done by setting $\Theta=\\{\theta\in\mathbb{R}^{d}\>,\>\|\theta\|\leq
R\\}$ for some radius $R>0$ [33, 1, 62].

In this case, the projected iteration (19) may be used. In this setting, one
may easily check that the gradient is sub-Gaussian/sub-exponential as soon as
the covariates $X$ satisfy one or the other of these properties. Therefore,
the results of Propositions 2 and 3 apply in this context as well.

##  8 Conclusion and Discussion

The Markov chain point of view for SGD is very useful since it allows to draw
conclusions and establish a number of characterizations for the invariant
limit distribution. Convergence of the SGD Markov chain holds under fairly
weak conditions [58]. As evidenced by our results, this opens doors for a
better characterization of the limit distribution when the associated
optimization iteration progresses at _geometric_ speed, for instance, when
strong convexity holds. The precise determination of the speed of convergence
in distribution constitutes a particular difficulty which more generally
concerns Markov chains with a geometric drift property. However, this
difficulty may be circumvented for SGD by leveraging Wasserstein convergence
provided a regularity condition on the noise distribution and transition
kernel. Obtaining such properties from generic assumptions on the gradient
distribution represents an interesting perspective.

Finally, despite being quite productive, the Markov chain study of SGD remains
limited to the constant step-size setting. This excludes the combination of a
decreasing step-size with averaging which is known for its better dependence
on problem conditioning [1, 2].

Acknowledgments

This research is supported by the Agence Nationale de la Recherche as part of
the “Investissements d’avenir” program (reference ANR-19-P3IA-0001; PRAIRIE
3IA Institute).

## References

  * [1] Francis Bach.  Adaptivity of averaged stochastic gradient descent to local strong convexity for logistic regression.  The Journal of Machine Learning Research, 15(1):595–627, 2014.
  * [2] Francis Bach and Eric Moulines.  Non-strongly-convex smooth stochastic approximation with convergence rate $O(1/n)$.  Advances in Neural Information Processing Systems, 26:773–781, 2013\.
  * [3] Franck Barthe and Cyril Roberto.  Modified logarithmic Sobolev inequalities on $\mathbb{R}$.  Potential Analysis, 29(2):167–193, 2008.
  * [4] Peter H Baxendale.  Renewal theory and computable convergence rates for geometrically ergodic Markov chains.  The Annals of Applied Probability, 15(1B):700–738, 2005.
  * [5] Witold Bednorz.  The Kendall theorem and its application to the geometric ergodicity of Markov chains.  Applicationes Mathematicae, 40(2):129–165, 2013.
  * [6] Kenneth S Berenhaut and Robert Lund.  Geometric renewal convergence rates from hazard rates.  Journal of applied probability, 38(1):180–194, 2001.
  * [7] Sergej G Bobkov and Friedrich Götze.  Exponential integrability and transportation cost related to logarithmic Sobolev inequalities.  Journal of Functional Analysis, 163(1):1–28, 1999.
  * [8] Sergey Bobkov and Michel Ledoux.  Poincaré’s inequalities and Talagrand’s concentration phenomenon for the exponential distribution.  Probability Theory and Related Fields, 107:383–400, 1997.
  * [9] Léon Bottou and Olivier Bousquet.  The tradeoffs of large scale learning.  Advances in Neural Information Processing Systems, 20:161–168, 2007\.
  * [10] Stéphane Boucheron, Gábor Lugosi, and Pascal Massart.  Concentration Inequalities: A Nonasymptotic Theory of Independence.  Oxford university press, 2013.
  * [11] Sébastien Bubeck.  Convex optimization: Algorithms and complexity.  Foundations and Trends® in Machine Learning, 8(3-4):231–357, 2015.
  * [12] Fan Chung and Linyuan Lu.  Concentration inequalities and martingale inequalities: a survey.  Internet mathematics, 3(1):79–127, 2006.
  * [13] Kai Lai Chung.  On a stochastic approximation method.  The Annals of Mathematical Statistics, pages 463–483, 1954.
  * [14] Alexandre Défossez and Francis Bach.  Averaged least-mean-squares: Bias-variance trade-offs and optimal sampling distributions.  In Artificial Intelligence and Statistics, pages 205–213. PMLR, 2015.
  * [15] Persi Diaconis and Laurent Saloff-Coste.  Comparison theorems for reversible Markov chains.  The Annals of Applied Probability, 3(3):696–730, 1993.
  * [16] Persi Diaconis and Daniel Stroock.  Geometric Bounds for Eigenvalues of Markov Chains.  The Annals of Applied Probability, 1(1):36–61, 1991.
  * [17] Aymeric Dieuleveut and Francis Bach.  Nonparametric stochastic approximation with large step-sizes.  The Annals of Statistics, 44(4):1363–1399, 2016.
  * [18] Aymeric Dieuleveut, Alain Durmus, and Francis Bach.  Bridging the gap between constant step size stochastic gradient descent and Markov chains.  The Annals of Statistics, 48(3):1348–1382, 2020.
  * [19] Aymeric Dieuleveut, Nicolas Flammarion, and Francis Bach.  Harder, better, faster, stronger convergence rates for least-squares regression.  The Journal of Machine Learning Research, 18(1):3520–3570, 2017\.
  * [20] Randal Douc, Eric Moulines, and Jeff Rosenthal.  Quantitative bounds for geometric convergence rates of Markov chains.  Ann. Appl. Probab, 14(4):1643–1665, 2004.
  * [21] Devdatt P Dubhashi and Alessandro Panconesi.  Concentration of Measure for the Analysis of Randomized Algorithms.  Cambridge University Press, 2009.
  * [22] Vaclav Fabian.  On asymptotic normality in stochastic approximation.  The Annals of Mathematical Statistics, pages 1327–1332, 1968.
  * [23] David A Freedman.  On tail probabilities for martingales.  the Annals of Probability, pages 100–118, 1975.
  * [24] Julia Gaudio, Saurabh Amin, and Patrick Jaillet.  Exponential convergence rates for stochastically ordered Markov processes under perturbation.  Systems & Control Letters, 133:104515, 2019.
  * [25] Ivan Gentil, Arnaud Guillin, and Laurent Miclo.  Modified logarithmic Sobolev inequalities and transportation inequalities.  Probability Theory and Related Fields, 133:409–436, 2005.
  * [26] Saeed Ghadimi and Guanghui 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.
  * [27] Eduard Gorbunov, Marina Danilova, and Alexander Gasnikov.  Stochastic optimization with heavy-tailed noise via accelerated gradient clipping.  Advances in Neural Information Processing Systems, 33:15042–15053, 2020.
  * [28] Nathael Gozlan.  A characterization of dimension free concentration in terms of transportation inequalities.  The Annals of Probability, 37(6):2480–2498, 2009.
  * [29] Nathael Gozlan and Christian Léonard.  Transport Inequalities. A Survey.  Markov Processes and Related Fields, 16:635–736, 2010.
  * [30] László Györfi and Harro Walk.  On the averaged stochastic approximation for linear regression.  SIAM Journal on Control and Optimization, 34(1):31–61, 1996.
  * [31] Nicholas JA Harvey, Christopher Liaw, Yaniv Plan, and Sikander Randhawa.  Tight analyses for non-smooth stochastic gradient descent.  In Conference on Learning Theory, pages 1579–1613. PMLR, 2019.
  * [32] Nicholas JA Harvey, Christopher Liaw, and Sikander Randhawa.  Simple and optimal high-probability bounds for strongly-convex stochastic gradient descent.  arXiv preprint arXiv:1909.00843, 2019.
  * [33] Elad Hazan, Tomer Koren, and Kfir Y Levy.  Logistic regression: Tight bounds for stochastic and online optimization.  In Conference on Learning Theory, pages 197–209. PMLR, 2014.
  * [34] Prateek Jain, Sham M Kakade, Rahul Kidambi, Praneeth Netrapalli, Venkata Krishna Pillutla, and Aaron Sidford.  A Markov chain theory approach to characterizing the minimax optimality of stochastic gradient descent (for least squares).  arXiv preprint arXiv:1710.09430, 2017.
  * [35] Prateek Jain, Sham M. Kakade, Rahul Kidambi, Praneeth Netrapalli, and Aaron Sidford.  Parallelizing stochastic gradient descent for least squares regression: Mini-batching, averaging, and model misspecification.  Journal of Machine Learning Research, 18(223):1–42, 2018.
  * [36] Prateek Jain, Dheeraj Nagaraj, and Praneeth Netrapalli.  Making the last iterate of SGD information theoretically optimal.  In Conference on Learning Theory, pages 1752–1755. PMLR, 2019.
  * [37] Daniel C Jerison.  Quantitative convergence rates for reversible Markov chains via strong random times.  arXiv preprint arXiv:1908.06459, 2019.
  * [38] Sham M Kakade and Ambuj Tewari.  On the generalization ability of online strongly convex programming algorithms.  Advances in Neural Information Processing Systems, 21, 2008.
  * [39] Hamed Karimi, Julie Nutini, and Mark Schmidt.  Linear convergence of gradient and proximal-gradient methods under the Polyak-Łojasiewicz condition.  In Machine Learning and Knowledge Discovery in Databases: European Conference, ECML PKDD 2016, Riva del Garda, Italy, September 19-23, 2016, Proceedings, Part I 16, pages 795–811. Springer, 2016.
  * [40] David G Kendall.  Unitary dilations of Markov transition operators, and the corresponding integral representations for transition-probability matrices.  Probability and statistics, pages 139–161, 1959.
  * [41] Masaaki Kijima.  Markov Processes for Stochastic Modeling, volume 6.  CRC Press, 1997.
  * [42] Aryeh Kontorovich and Maxim Raginsky.  Concentration of measure without independence: a unified approach via the martingale method.  In Convexity and Concentration, pages 183–210. Springer, 2017.
  * [43] Ioannis Kontoyiannis and Sean P Meyn.  Geometric ergodicity and the spectral gap of non-reversible Markov chains.  Probability Theory and Related Fields, 154(1-2):327–339, 2012.
  * [44] Solomon Kullback and Richard A Leibler.  On information and sufficiency.  The Annals of Mathematical Statistics, 22(1):79–86, 1951.
  * [45] Simon Lacoste-Julien, Mark Schmidt, and Francis Bach.  A simpler approach to obtaining an $O(1/t)$ convergence rate for the projected stochastic subgradient method.  arXiv preprint arXiv:1212.2002, 2012.
  * [46] Chandrashekar Lakshminarayanan and Csaba Szepesvari.  Linear stochastic approximation: How far does constant step-size and iterate averaging go?  In International Conference on Artificial Intelligence and Statistics, pages 1347–1355. PMLR, 2018.
  * [47] Michel Ledoux.  On Talagrand’s deviation inequalities for product measures.  ESAIM: Probability and Statistics, 1:63–87, 1997.
  * [48] Xiaoyu Li and Francesco Orabona.  A high probability analysis of adaptive SGD with momentum.  arXiv preprint arXiv:2007.14294, 2020.
  * [49] Zhipeng Lou, Wanrong Zhu, and Wei Biao Wu.  Beyond sub-Gaussian noises: Sharp concentration analysis for stochastic gradient descent.  Journal of Machine Learning Research, 23:1–22, 2022.
  * [50] Gábor Lugosi and Shahar Mendelson.  Mean estimation and regression under heavy-tailed distributions: A survey.  Foundations of Computational Mathematics, 19(5):1145–1190, 2019\.
  * [51] Robert Lund, Ying Zhao, and Peter C Kiessler.  A monotonicity in reversible Markov chains.  Journal of Applied Probability, 43(2):486–499, 2006.
  * [52] Robert B Lund, Sean P Meyn, and Richard L Tweedie.  Computable exponential convergence rates for stochastically ordered Markov processes.  The Annals of Applied Probability, 6(1):218–237, 1996.
  * [53] Robert B Lund and Richard L Tweedie.  Geometric convergence rates for stochastically ordered Markov chains.  Mathematics of Operations Research, 21(1):182–194, 1996.
  * [54] Siyuan Ma, Raef Bassily, and Mikhail Belkin.  The Power of Interpolation: Understanding the Effectiveness of SGD in Modern Over-parametrized Learning.  In International Conference on Machine Learning, 2017.
  * [55] Neal Madras and Deniz Sezer.  Quantitative bounds for Markov chain convergence: Wasserstein and total variation distances.  Bernoulli, 16(3):882–908, 2010.
  * [56] Andreas Maurer and Massimiliano Pontil.  Concentration inequalities under sub-Gaussian and sub-exponential conditions.  Advances in Neural Information Processing Systems, 34:7588–7597, 2021.
  * [57] Colin McDiarmid.  Concentration.  Probabilistic Methods for Algorithmic Discrete Mathematics, pages 195–248, 1998.
  * [58] Sean P Meyn and Richard L Tweedie.  Markov Chains and Stochastic Stability.  Springer London, 1993.
  * [59] Sean P Meyn and Robert L Tweedie.  Computable bounds for geometric convergence rates of Markov chains.  The Annals of Applied Probability, pages 981–1011, 1994.
  * [60] Wenlong Mou, Chris Junchi Li, Martin J Wainwright, Peter L Bartlett, and Michael I Jordan.  On linear stochastic approximation: Fine-grained Polyak-Ruppert and non-asymptotic concentration.  In Conference on Learning Theory, pages 2947–2997. PMLR, 2020.
  * [61] Eric Moulines and Francis Bach.  Non-asymptotic analysis of stochastic approximation algorithms for machine learning.  Advances in Neural Information Processing Systems, 24:451–459, 2011\.
  * [62] Jaouad Mourtada and Stéphane Gaïffas.  An improper estimator with optimal excess risk in misspecified density estimation and logistic regression.  Journal of Machine Learning Research, 23(31):1–49, 2022.
  * [63] Angelia Nedić and Dimitri Bertsekas.  Convergence rate of incremental subgradient algorithms.  Stochastic Optimization: Algorithms and Applications, pages 223–264, 2001.
  * [64] Deanna Needell, Rachel Ward, and Nati Srebro.  Stochastic gradient descent, weighted sampling, and the randomized Kaczmarz algorithm.  Advances in Neural Information Processing Systems, 27, 2014.
  * [65] Yurii Nesterov.  Introductory Lectures on Convex Optimization: A Basic Course.  Springer Publishing Company, Incorporated, 1 edition, 2014.
  * [66] Gergely Neu and Lorenzo Rosasco.  Iterate averaging as regularization for stochastic gradient descent.  In Conference On Learning Theory, pages 3222–3242. PMLR, 2018.
  * [67] Ta Duy Nguyen, Thien Hang Nguyen, Alina Ene, and Huy Le Nguyen.  High probability convergence of clipped-sgd under heavy-tailed noise.  arXiv preprint arXiv:2302.05437, 2023.
  * [68] Georg Ch Pflug.  Stochastic minimization with constant step-size: asymptotic laws.  SIAM Journal on Control and Optimization, 24(4):655–666, 1986.
  * [69] Boris T Polyak and Anatoli B Juditsky.  Acceleration of stochastic approximation by averaging.  SIAM journal on control and optimization, 30(4):838–855, 1992.
  * [70] Qian Qin and James P Hobert.  On the limitations of single-step drift and minorization in Markov chain convergence analysis.  The Annals of Applied Probability, 31(4):1633–1659, 2021.
  * [71] Qian Qin and James P Hobert.  Wasserstein-based methods for convergence complexity analysis of MCMC with applications.  The Annals of Applied Probability, 32(1):124–166, 2022.
  * [72] Alexander Rakhlin, Ohad Shamir, and Karthik Sridharan.  Making gradient descent optimal for strongly convex stochastic optimization.  In Proceedings of the 29th International Coference on International Conference on Machine Learning, page 1571–1578. Omnipress, 2012\.
  * [73] Herbert Robbins and Sutton Monro.  A stochastic approximation method.  The Annals of Mathematical Statistics, pages 400–407, 1951.
  * [74] Gareth O Roberts and Richard L Tweedie.  Bounds on regeneration times and convergence rates for Markov chains.  Stochastic Processes and their applications, 80(2):211–229, 1999\.
  * [75] Gareth O Roberts and Richard L Tweedie.  Rates of convergence of stochastically monotone and continuous time Markov models.  Journal of Applied Probability, 37(2):359–373, 2000.
  * [76] Gareth O Roberts and Richard L Tweedie.  Geometric $L_{2}$ and $L_{1}$ convergence are equivalent for reversible Markov chains.  Journal of Applied Probability, 38(A):37–41, 2001.
  * [77] Jeffrey Rosenthal.  Quantitative convergence rates of Markov chains: A simple account.  Electronic Communications in Probability, 7:123–128, 2002.
  * [78] Jeffrey S Rosenthal.  Convergence rates for Markov chains.  Siam Review, 37(3):387–405, 1995.
  * [79] David Ruppert.  Efficient estimations from a slowly convergent Robbins-Monro process.  Technical report, Cornell University Operations Research and Industrial Engineering, 1988.
  * [80] Jerome Sacks.  Asymptotic distribution of stochastic approximation procedures.  The Annals of Mathematical Statistics, 29(2):373–405, 1958.
  * [81] Abdurakhmon Sadiev, Marina Danilova, Eduard Gorbunov, Samuel Horváth, Gauthier Gidel, Pavel Dvurechensky, Alexander Gasnikov, and Peter Richtárik.  High-probability bounds for stochastic optimization and variational inequalities: the case of unbounded variance.  arXiv preprint arXiv:2302.00999, 2023.
  * [82] Mark Schmidt and Nicolas Le Roux.  Fast convergence of stochastic gradient descent under a strong growth condition.  arXiv preprint arXiv:1308.6370, 2013.
  * [83] DJ Scott and RL Tweedie.  Explicit rates of convergence of stochastically ordered Markov chains.  In Athens Conference on Applied Probability and Time Series Analysis, pages 176–191. Springer, 1996.
  * [84] Ohad Shamir and Tong Zhang.  Stochastic gradient descent for non-smooth optimization: Convergence results and optimal averaging schemes.  In International Conference on Machine Learning, pages 71–79. PMLR, 2013.
  * [85] Che-Ping Tsai, Adarsh Prasad, Sivaraman Balakrishnan, and Pradeep Ravikumar.  Heavy-tailed streaming statistical estimation.  In International Conference on Artificial Intelligence and Statistics, pages 1251–1282. PMLR, 2022.
  * [86] Ramon Van Handel.  Probability in high dimension.  Technical report, Pinceton University, 2014.
  * [87] Sharan Vaswani, Francis Bach, and Mark Schmidt.  Fast and faster convergence of sgd for over-parameterized models and an accelerated perceptron.  In Proceedings of the Twenty-Second International Conference on Artificial Intelligence and Statistics, volume 89, pages 1195–1204. PMLR, 2019\.
  * [88] Roman Vershynin.  Introduction to the non-asymptotic analysis of random matrices.  arXiv preprint arXiv:1011.3027, 2010.
  * [89] Roman Vershynin.  High-dimensional probability: An introduction with applications in data science, volume 47.  Cambridge university press, 2018.
  * [90] Cédric Villani et al.  Optimal Transport: Old and New.  Grundlehren der mathematischen Wissenschaften. Springer Berlin, Heidelberg, 2009.
  * [91] Lu Yu, Krishnakumar Balasubramanian, Stanislav Volgushev, and Murat A Erdogdu.  An analysis of constant step size SGD in the non-convex regime: Asymptotic normality and bias.  Advances in Neural Information Processing Systems, 34:4234–4248, 2021.

##  Appendix A Preliminary lemmas

######  Lemma 3.

Let $X$ be a real random variable such that $X\in\widetilde{\Psi}_{2}(K)$
then, with probability at least $\delta,$ we have

$|X|\leq K\sqrt{\log(e/\delta)}.$

###### Proof.

Using Chernoff’s method, we find for $t>0$ and $\lambda>0$

$\displaystyle\mathbb{P}\big{(}|X|>t\big{)}$
$\displaystyle=\mathbb{P}\big{(}\lambda^{2}X^{2}>\lambda^{2}t^{2}\big{)}=\mathbb{P}\big{(}\exp(\lambda^{2}X^{2})>\exp(\lambda^{2}t^{2})\big{)}$
$\displaystyle\leq\mathbb{E}\exp(\lambda^{2}X^{2})e^{-\lambda^{2}t^{2}}\leq\exp\big{(}\lambda^{2}(K^{2}-t^{2})\big{)}.$

Choosing $\lambda\\!=\\!1/K,$ we have
$\exp\big{(}1\\!-\\!(t/K)^{2}\big{)}\\!\leq\\!\delta\\!\\!\iff\\!\\!t\\!\geq\\!K\sqrt{\log(e/\delta)}$
and the result follows. ∎

######  Lemma 4.

Let $X$ be a real random variable such that $X\in\widetilde{\Psi}_{1}(K)$
then, with probability at least $\delta,$ we have

$|X|\leq 2eK\log(2/\delta).$

###### Proof.

Using Stirling’s approximation, we find for $|\lambda|<(eK)^{-1}:$

$\displaystyle\mathbb{E}\exp(\lambda|X|)$ $\displaystyle=\sum_{p\geq
0}\frac{\lambda^{p}\mathbb{E}|X|^{p}}{p!}\leq 1+\sum_{p\geq 1}\frac{(\lambda
Kp)^{p}}{p!}$ $\displaystyle\leq 1+\sum_{p\geq 1}\frac{(\lambda
eK)^{p}}{\sqrt{2\pi p}}\leq 1+\frac{1}{\sqrt{2\pi}}\frac{\lambda eK}{1-\lambda
eK}\leq\exp\Big{(}\frac{1}{\sqrt{2\pi}}\frac{\lambda eK}{1-\lambda
eK}\Big{)},$

where we used the inequality $1+x\leq e^{x}$ in the last step. For $t>0,$
using Chernoff’s method and choosing $\lambda=(2eK)^{-1}$, we find:

$\displaystyle\mathbb{P}\big{(}|X|>t\big{)}$
$\displaystyle=\mathbb{P}\big{(}\lambda|X|>\lambda
t\big{)}=\mathbb{P}\big{(}\exp(\lambda|X|)>\exp(\lambda t)\big{)}$
$\displaystyle\leq\mathbb{E}\exp(\lambda|X|)e^{-\lambda
t}\leq\exp\Big{(}\frac{1}{\sqrt{2\pi}}-\frac{t}{2eK}\Big{)}.$

It only remains to choose $t=2eK\log(2/\delta)$ to obtain the desired bound. ∎

The following fundamental lemma will be often used in our proofs.

######  Lemma 5.

Grant Assumption 1. For any $\theta,\theta^{\prime}\in\mathbb{R}^{d}$ and
$\beta\leq\frac{2}{\mu+L}$ we have

$\|\theta-\beta\nabla\mathcal{L}(\theta)-(\theta^{\prime}-\beta\nabla\mathcal{L}(\theta^{\prime}))\|^{2}\leq(1-\beta\mu)^{2}\|\theta-\theta^{\prime}\|^{2}.$
(20)

###### Proof.

For $\beta\leq\frac{2}{\mu+L},$ we have

$\displaystyle\|\theta-\beta\nabla\mathcal{L}(\theta)-$
$\displaystyle(\theta^{\prime}-\beta\nabla\mathcal{L}(\theta^{\prime}))\|^{2}$
$\displaystyle=\|\theta-\theta^{\prime}\|^{2}-2\beta\langle\theta-\theta^{\prime},\nabla\mathcal{L}(\theta)-\nabla\mathcal{L}(\theta^{\prime})\rangle+\beta^{2}\|\nabla\mathcal{L}(\theta)-\nabla\mathcal{L}(\theta^{\prime})\|^{2}$
$\displaystyle\leq(1-\beta^{2}\mu
L)\|\theta-\theta^{\prime}\|^{2}-\beta(2-\beta(\mu+L))\langle\nabla\mathcal{L}(\theta)-\nabla\mathcal{L}(\theta^{\prime}),\theta-\theta^{\prime}\rangle)$
$\displaystyle\leq(1-\beta^{2}\mu
L)\|\theta-\theta^{\prime}\|^{2}-\beta(2-\beta(\mu+L))\mu\|\theta-\theta^{\prime}\|^{2}$
$\displaystyle=(1-\beta^{2}\mu
L-2\beta\mu+\beta^{2}\mu(\mu+L))\|\theta-\theta^{\prime}\|^{2}$
$\displaystyle=(1-\beta\mu)^{2}\|\theta-\theta^{\prime}\|^{2},$

where we used the inequalities

$\displaystyle\|\nabla\mathcal{L}(\theta)-\nabla\mathcal{L}(\theta^{\prime})\|^{2}$
$\displaystyle\leq(\mu+L)\langle\nabla\mathcal{L}(\theta)-\nabla\mathcal{L}(\theta^{\prime}),\theta-\theta^{\prime}\rangle-\mu
L\|\theta-\theta^{\prime}\|^{2}$ (21)
$\displaystyle\mu\|\theta-\theta^{\prime}\|^{2}$
$\displaystyle\leq\langle\nabla\mathcal{L}(\theta)-\nabla\mathcal{L}(\theta^{\prime}),\theta-\theta^{\prime}\rangle,$
(22)

valid for all $\theta,\theta^{\prime}$. Equation (21) is stated, for example,
in [65, Theoerem 2.1.12] (see also [11, Lemma 3.11] and (22) is just a
characterization of strong convexity (see for instance [65, Theorem 2.1.9]). ∎

##  Appendix B Proof of Geometric Ergodicity

In the remaining part of this document, we make the dependences on $\zeta$ in
the gradient samples and errors implicit and write $G(\theta)$ and
$\varepsilon(\theta)$ instead of $G(\theta,\zeta)$ and
$\varepsilon_{\zeta}(\theta)$ respectively.

We show the geometric ergodicity of the SGD Markov chain $(\theta_{t})_{t\geq
0}$ by relying on [58, Theorem 15.0.1]. We will show that the following
function:

$V(\theta):=1+\|\theta-\theta^{\star}\|^{2},$

is a _drift_ function for this Markov chain. We define the action of the
transition kernel $P$ on integrable functions $f$ through

$\displaystyle Pf(\theta)=\mathbb{E}f(\theta-\beta G(\theta)).$

We also define the variation operator

$\Delta f(\theta):=Pf(\theta)-f(\theta).$

###  B.1 Proof of Theorem 1

First, we establish that the Markov chain is aperiodic. Indeed, by Assumption
2, for all $\theta,$ the gradient is distributed according to an everywhere
positive density, therefore, for all $\theta\in S\subset\mathbb{R}^{d}$ with
$S$ a set with non zero Lebesgue measure we have $P(\theta,S)>0.$ This implies
that the greatest possible period for the chain is $1$ which makes it
aperiodic.

We also show that the Markov chain is $\psi$-irreducible (see [58, Chapter
4]). For any initial $\theta_{0},$ its successor reads:

$\theta_{+1}=\theta_{0}-\beta(\nabla\mathcal{L}(\theta_{0})+\varepsilon(\theta_{0}))$

Given Assumption 2, the distribution of $\varepsilon(\theta_{0})$ is minorized
by $\delta\nu_{\theta_{0},1}$ where $\nu_{\theta_{0},1}$ is a probability
distribution which admits an everywhere positive density
$h(\theta_{0},\cdot).$ Consequently, for all $A\in\mathcal{B}(\mathbb{R}^{d})$
with non zero Lebesgue measure, we have the following minorization:

$\mathbb{P}(\theta_{+1}\in
A|\theta_{0})=P(\theta_{0},A)\geq\frac{\delta}{\beta^{d}}\int_{A}h\Big{(}\theta_{0},\frac{\theta-\theta_{0}}{\beta}-\nabla\mathcal{L}(\theta_{0})\Big{)}d\theta>0.$

It follows that the Markov chain is irreducible w.r.t. Lebesgue’s measure and
is thus $\psi$-irreducible.

For fixed $\theta,$ and step-size $\beta<\frac{2}{\mu+L},$ using Lemma 5 we
find:

$\displaystyle P\|\theta-\theta^{\star}\|^{2}$
$\displaystyle=\mathbb{E}\|\theta-\beta G(\theta)-\theta^{\star}\|^{2}$
$\displaystyle=\mathbb{E}\big{[}\|\theta-\beta\nabla\mathcal{L}(\theta)-\theta^{\star}\|^{2}-2\beta\langle\theta-\beta\nabla\mathcal{L}(\theta)-\theta^{\star},\varepsilon(\theta)\rangle+\beta^{2}\|\varepsilon(\theta)\|^{2}\big{]}$
$\displaystyle\leq(1-\beta\mu)^{2}\|\theta-\theta^{\star}\|^{2}+\beta^{2}\mathbb{E}\|\varepsilon(\theta)\|^{2}$
$\displaystyle\leq(1-\beta\mu)^{2}\|\theta-\theta^{\star}\|^{2}+\beta^{2}\big{(}L_{\sigma}\|\theta-\theta^{\star}\|^{2}+\sigma^{2}\big{)}$
$\displaystyle=\big{(}(1-\beta\mu)^{2}+\beta^{2}L_{\sigma}\big{)}\|\theta-\theta^{\star}\|^{2}+\beta^{2}\sigma^{2}$

The previous inequality yields a contraction for step-size satisfying
$0<\beta<\frac{2\mu}{\mu^{2}+L_{\sigma}}$ and, as a consequence we have:

$PV(\theta)\leq\underbrace{\big{(}(1-\beta\mu)^{2}+\beta^{2}L_{\sigma}\big{)}}_{=:\widetilde{\lambda}}V(\theta)+\underbrace{\beta^{2}\sigma^{2}+\big{(}1-\big{(}(1-\beta\mu)^{2}+\beta^{2}L_{\sigma}\big{)}\big{)}}_{=:\widetilde{b}}$

We now define the set
$\mathcal{C}=\big{\\{}\theta\in\mathbb{R}^{d}\>,\>V(\theta)\leq
2\widetilde{b}/(1-\widetilde{\lambda})\big{\\}}$ which satisfies:

$\Delta
V(\theta)\leq-\frac{1-\widetilde{\lambda}}{2}V(\theta)+\widetilde{b}\mathbf{1}_{\theta\in\mathcal{C}}.$

For such $\mathcal{C},$ let
$\underline{h}(\theta)=\inf_{\theta_{0}\in\mathcal{C}}h\Big{(}\theta_{0},\frac{\theta-\theta_{0}}{\beta}-\nabla\mathcal{L}(\theta_{0})\Big{)}$
and define the probability measure $\nu_{\mathcal{C}}$ by

$\nu_{\mathcal{C}}(A)=\frac{\int_{A\cap\mathcal{C}}\underline{h}(\theta)d\theta}{\int_{\mathcal{C}}\underline{h}(\theta^{\prime})d\theta^{\prime}}\quad\text{
for all }\quad A\in\mathcal{B}(\mathbb{R}^{d}).$

It follows that for all $\theta_{0}\in\mathcal{C},$ we have the following
minorization property:

$P(\theta_{0},A)\geq\xi\nu_{\mathcal{C}}(A)\quad\text{ for all }\quad
A\in\mathcal{B}(\mathbb{R}^{d}),$

where $\xi=\delta\int_{\mathcal{C}}\underline{h}(\theta)d\theta>0.$ In words,
the set $\mathcal{C}$ is a _small_ set and, thanks to [58, Proposition 5.5.3],
also a _petite_ set (see definitions in [58, Chapter 5]).

We have shown that the Markov chain $(\theta_{t})$ satisfies condition (iii)
of [58, Theorem 15.0.1]. By the latter result, it follows that it admits an
invariant probability measure $\pi_{\beta}$ and there exist $r>1$ and
$M<\infty$ such that:

$\sum_{t\geq 0}r^{t}\|P(\theta_{0},\cdot)-\pi_{\beta}\|_{V}\leq
MV(\theta_{0}).$ (23)

Taking $\rho=r^{-1}$ allows to conclude the proof.

###  B.2 Proof of Proposition 1

To prove (a), let $\theta\sim\pi_{\beta}$ and simply compute

$\displaystyle\mathbb{E}[\theta]=\mathbb{E}\big{[}\theta-\beta
G(\theta)\big{]}=\mathbb{E}\big{[}\theta-\beta\nabla\mathcal{L}(\theta)\big{]}=\mathbb{E}[\theta]-\beta\mathbb{E}[\nabla\mathcal{L}(\theta)]$

since we know that $\mathbb{E}[\theta]<\infty$ (this follows from (23) in the
proof of Theorem 1), this implies the first part of the claim. If we further
assume the gradient to be linear, we have in addition

$\mathbb{E}\nabla\mathcal{L}(\theta)=\nabla\mathcal{L}\big{(}\mathbb{E}\theta\big{)}=\nabla\mathcal{L}(\bar{\theta})=0,$

and the conclusion follows since $\theta^{\star}$ is the unique critical
point.

To prove (b), let $\theta\sim\pi_{\beta}$ which implies that $\theta-\beta
G(\theta)\sim\pi_{\beta}$ as well. We compute:

$\displaystyle\mathbb{E}\|\theta-\theta^{\star}\|^{2}$
$\displaystyle=\mathbb{E}\|\theta-\beta G(\theta)-\theta^{\star}\|^{2}$
$\displaystyle=\mathbb{E}\big{[}\|\theta-\beta\nabla\mathcal{L}(\theta)-\theta^{\star}\|^{2}+\beta^{2}\|\varepsilon(\theta)\|^{2}-2\beta\langle\theta-\beta\nabla\mathcal{L}(\theta)-\theta^{\star},\varepsilon(\theta)\rangle\big{]}$
$\displaystyle=\mathbb{E}\big{[}\|\theta-\beta\nabla\mathcal{L}(\theta)-\theta^{\star}\|^{2}+\beta^{2}\mathbb{E}\|\varepsilon(\theta)\|^{2}\big{]}$
$\displaystyle\leq\big{(}(1-\beta\mu)^{2}+\beta^{2}L_{\sigma}\big{)}\mathbb{E}\|\theta-\theta^{\star}\|^{2}+\beta^{2}\sigma^{2}$

which implies

$\mathbb{E}\|\theta-\theta^{\star}\|^{2}\leq\frac{\beta^{2}\sigma^{2}}{1-(1-\beta\mu)^{2}-\beta^{2}L_{\sigma}}=\frac{\beta\sigma^{2}}{2\mu-\beta(\mu^{2}+L_{\sigma})}.$

Moreover, by Jensen and Cauchy-Schwarz inequalities, we have

$\|\bar{\theta}_{\beta}-\theta^{\star}\|=\|\mathbb{E}\theta-\theta^{\star}\|\leq\mathbb{E}\|\theta-\theta^{\star}\|\leq\sqrt{\mathbb{E}\|\theta-\theta^{\star}\|^{2}}.$

###  B.3 Proof of Proposition 2

Let us prove (a), we again use the fact that $\theta$ and $\theta-\beta
G(\theta)$ have the same distribution when $\theta\sim\pi_{\beta}:$

$\displaystyle\mathbb{E}\exp(\lambda^{2}\|\theta-\theta^{\star}\|^{2})$
$\displaystyle=\mathbb{E}\exp(\lambda^{2}\|\theta-\beta
G(\theta)-\theta^{\star}\|^{2})$
$\displaystyle=\mathbb{E}\exp\big{(}\lambda^{2}(\|\theta\\!-\\!\beta\nabla\mathcal{L}(\theta)\\!-\\!\theta^{\star}\|^{2}\\!-\\!2\beta\langle\theta\\!-\\!\beta\nabla\mathcal{L}(\theta)\\!-\\!\theta^{\star},\varepsilon(\theta)\rangle\\!+\\!\beta^{2}\|\varepsilon(\theta)\|^{2})\big{)}$

Since we assume that
$\|\varepsilon(\theta)\|\in\widetilde{\Psi}_{2}(\overline{K}),$ it is easy to
check that for all $u\in\mathbb{R}^{d}$ with unit norm, $\langle
u,\varepsilon(\theta)\rangle\in\Psi_{2}(\overline{K})$ because
$\varepsilon(\theta)$ is centered. Therefore, conditioning on $\theta,$ we
have

$\displaystyle\mathbb{E}\big{[}\exp$
$\displaystyle(\lambda^{2}(-2\beta\langle\theta-\beta\nabla\mathcal{L}(\theta)-\theta^{\star},\varepsilon(\theta)\rangle+\beta^{2}\|\varepsilon(\theta)\|^{2}))|\theta\big{]}$
$\displaystyle\leq\mathbb{E}\big{[}\exp(-(2\lambda)^{2}\beta\langle\theta-\beta\nabla\mathcal{L}(\theta)-\theta^{\star},\varepsilon(\theta)\rangle|\theta\big{]}^{1/2}\mathbb{E}\big{[}\exp(2(\beta\lambda)^{2}\|\varepsilon(\theta)\|^{2})|\theta\big{]}^{1/2}$
$\displaystyle\leq\exp(8\lambda^{4}\beta^{2}\|\theta-\beta\nabla\mathcal{L}(\theta)-\theta^{\star}\|^{2}\overline{K}^{2}+\lambda^{2}\beta^{2}\overline{K}^{2})$
$\displaystyle\leq\exp(8\lambda^{4}\beta^{2}(1-\beta\mu)^{2}\|\theta-\theta^{\star}\|^{2}\overline{K}^{2}+\lambda^{2}\beta^{2}\overline{K}^{2}),$

where the last line uses that
$\|\theta-\beta\nabla\mathcal{L}(\theta)-\theta^{\star}\|\leq(1-\beta\mu)\|\theta-\theta^{\star}\|$.
The previous inequality holds for
$|\lambda|\leq\big{(}\sqrt{2}\beta\overline{K}\big{)}^{-1}$. We now restrict
$\lambda$ so that
$|\lambda|\leq\big{(}\overline{K}\sqrt{8\beta/\mu}\big{)}^{-1}$ and use
Jensen’s inequality to obtain

$\displaystyle\mathbb{E}\exp\big{(}\lambda^{2}\|\theta-\theta^{\star}\|^{2})$
$\displaystyle\leq\mathbb{E}\exp\big{(}\lambda^{2}(1-\beta\mu)^{2}\|\theta-\theta^{\star}\|^{2}(1+8\lambda^{2}\beta^{2}\overline{K}^{2})\big{)}\exp\big{(}\lambda^{2}\beta^{2}\overline{K}^{2}\big{)}$
$\displaystyle=\mathbb{E}\exp\big{(}\lambda^{2}(1-\beta\mu)^{2}\|\theta-\theta^{\star}\|^{2}(1+\beta\mu)\big{)}\exp\big{(}\lambda^{2}\beta^{2}\overline{K}^{2}\big{)}$
$\displaystyle=\mathbb{E}\exp\big{(}\lambda^{2}(1-\beta\mu)(1-(\beta\mu)^{2})\|\theta-\theta^{\star}\|^{2}\big{)}\exp\big{(}\lambda^{2}\beta^{2}\overline{K}^{2}\big{)}$
$\displaystyle\leq\Big{[}\mathbb{E}\exp\big{(}\lambda^{2}\|\theta-\theta^{\star}\|^{2}\big{)}\Big{]}^{(1-\beta\mu)(1-(\beta\mu)^{2})}\exp\big{(}\lambda^{2}\beta^{2}\overline{K}^{2}\big{)}$

and we arrive to the conclusion that

$\displaystyle\mathbb{E}\exp(\lambda^{2}\|\theta-\theta^{\star}\|^{2})$
$\displaystyle\leq\exp\Big{(}\frac{\lambda^{2}\beta^{2}\overline{K}^{2}}{\beta\mu+(\beta\mu)^{2}-(\beta\mu)^{3}}\Big{)}=\exp\Big{(}\frac{\lambda^{2}\beta\overline{K}^{2}}{\mu+\beta\mu^{2}-\beta^{2}\mu^{3}}\Big{)}$
$\displaystyle\leq\exp\Big{(}\frac{\lambda^{2}\beta\overline{K}^{2}}{\mu}\Big{)}$

where in the end we use that $\beta\mu^{2}\geq\beta^{2}\mu^{3}$ since
$\beta\leq 1/\mu$.

The task of proving (b) is more delicate. By assumption, we know that for all
$\theta$, the gradient error $\varepsilon(\theta)$ satisfies

$\big{\|}\|\varepsilon(\theta)\|\big{\|}_{L_{p}}\leq\overline{K}p\quad\text{for
}\quad p\geq 1.$ (24)

We denote $M_{p}^{p}=\mathbb{E}\|\theta-\theta^{\star}\|^{p}$. For $p=2,$ we
have through a similar computation to that for proving (b)

$M_{2}^{2}\leq\frac{\beta(2\overline{K})^{2}}{\mu(2-\beta\mu)}\leq\frac{\beta(2\overline{K})^{2}}{\mu},$

which immediately entails $M_{1}\leq 2\overline{K}\sqrt{\beta/\mu}.$ We will
show by induction that

$M_{p}\leq K_{\pi}p\quad\text{for all }\quad p\geq 1,$ (25)

with $K_{\pi}=C\overline{K}\sqrt{\beta/\mu}$ for some $C\geq 2.$ For $p\geq
2,$ we assume (25) holds up to $2p-2$ and consider $M_{2p}$, we compute

$\displaystyle M_{2p}^{2p}$
$\displaystyle=\mathbb{E}\|\theta-\theta^{\star}\|^{2p}=\mathbb{E}\|\theta-\beta
G(\theta)-\theta^{\star}\|^{2p}$
$\displaystyle=\mathbb{E}\big{(}\|\theta-\beta\nabla\mathcal{L}(\theta)-\theta^{\star}\|^{2}-2\beta\langle\theta-\beta\nabla\mathcal{L}(\theta)-\theta^{\star},\varepsilon(\theta)\rangle+\beta^{2}\|\varepsilon(\theta)\|^{2}\big{)}^{p}$
$\displaystyle\leq\mathbb{E}\sum_{k=0}^{p}\binom{p}{k}\|\theta-\beta\nabla\mathcal{L}(\theta)-\theta^{\star}\|^{2k}\big{(}\beta^{2}\|\varepsilon(\theta)\|^{2}-2\beta\langle\theta-\beta\nabla\mathcal{L}(\theta)-\theta^{\star},\varepsilon(\theta)\rangle\big{)}^{p-k}$
$\displaystyle\leq(1-\beta\mu)^{2p}M_{2p}^{2p}+p\mathbb{E}\|\theta-\beta\nabla\mathcal{L}(\theta)-\theta^{\star}\|^{2p-2}(\beta\|\varepsilon(\theta)\|)^{2}$
$\displaystyle\quad+\mathbb{E}\sum_{k=0}^{p-2}\binom{p}{k}\|\theta-\beta\nabla\mathcal{L}(\theta)-\theta^{\star}\|^{2k}\big{(}\beta^{2}\|\varepsilon(\theta)\|^{2}-2\beta\langle\theta-\beta\nabla\mathcal{L}(\theta)-\theta^{\star},\varepsilon(\theta)\rangle\big{)}^{p-k},$

where we isolated the two last terms of the sum in the last step and used that
$\varepsilon(\theta)$ is centered conditionally on $\theta.$ Further, we have

$\displaystyle\mathbb{E}$
$\displaystyle\sum_{k=0}^{p-2}\binom{p}{k}\|\theta-\beta\nabla\mathcal{L}(\theta)-\theta^{\star}\|^{2k}\big{(}\beta^{2}\|\varepsilon(\theta)\|^{2}-2\beta\langle\theta-\beta\nabla\mathcal{L}(\theta)-\theta^{\star},\varepsilon(\theta)\rangle\big{)}^{p-k}$
$\displaystyle\leq\mathbb{E}\sum_{k=0}^{p-2}\\!\binom{p}{k}\|\theta\\!-\\!\beta\nabla\mathcal{L}(\theta)\\!-\\!\theta^{\star}\|^{2k}\sum_{j=0}^{p-k}\\!\binom{p\\!-\\!k}{j}(\beta\|\varepsilon(\theta)\|)^{2j}(2\beta\|\theta\\!-\\!\beta\nabla\mathcal{L}(\theta)\\!-\\!\theta^{\star}\|\|\varepsilon(\theta)\|)^{p-k-j}$
$\displaystyle=\sum_{k=0}^{p-2}\sum_{j=0}^{p-k}\binom{p}{k}\binom{p-k}{j}\|\theta-\beta\nabla\mathcal{L}(\theta)-\theta^{\star}\|^{p+k-j}(\beta\|\varepsilon(\theta)\|)^{p-k+j}2^{p-k-j}.$

Now consider an index $l=p-k+j$, note that $2p-l=p+k-j$ and we have $2\leq
l\leq 2p.$ We compute the sum
$\sum_{k=0}^{p-2}\sum_{j=0}^{p-k}\binom{p}{k}\binom{p-k}{j}2^{p-k-j}$ for a
fixed value of $l:$

$\displaystyle\sum_{\begin{subarray}{c}k=0,j=0\\\
p-k+j=l\end{subarray}}^{k=p-2,j=p-k}$
$\displaystyle\binom{p}{k}\binom{p-k}{j}2^{p-k-j}=\sum_{k=0}^{p-2}\binom{p}{k}\binom{p-k}{l-(p-k)}2^{2(p-k)-l}\mathbf{1}_{p-k\leq
l\leq 2(p-k)}$ (26)
$\displaystyle=\\!\\!\\!\sum_{k=0\vee(p-l)}^{(p-2)\wedge(p-\lceil
l/2\rceil)}\\!\\!\binom{p}{k}\binom{p-k}{l\\!-\\!(p-k)}2^{2(p-k)-l}=\\!\\!\\!\sum_{k=0}^{(\lfloor
l/2\rfloor)\wedge(l-2)}\\!\\!\binom{p}{k,k\\!+\\!p\\!-\\!l,l\\!-\\!2k}2^{l-2k}$

where $\binom{p}{k,k+p-l,l-2k}=\frac{p!}{k!(k+p-l)!(l-2k)!}$ is the
_trinomial_ coefficient. Similarly we find that:

$\displaystyle\sum_{\begin{subarray}{c}k=0,j=0\\\
p-k+j=l\end{subarray}}^{k=p,j=p-k}\binom{p}{k}\binom{p-k}{j}2^{p-k-j}=\sum_{k=0}^{\lfloor
l/2\rfloor}\binom{p}{k,k\\!+\\!p\\!-\\!l,l\\!-\\!2k}2^{l-2k}.$

In what follows, we set the convention that $\binom{p}{k,k+p-l,l-2k}=0$
whenever $k\wedge(k+p-l)\wedge(l-2k)<0$ which allows us to sum over all
integer values without specifying the limits. For some variable $x,$ we
multiply by $x^{l},$ sum over $l$ and perform the change of variable $l\to
l+2k$ to find

$\displaystyle\sum_{l}$ $\displaystyle\sum_{k=0}^{\lfloor
l/2\rfloor}\binom{p}{k,k+p-l,l-2k}2^{l-2k}x^{l}=\sum_{l,k}\binom{p}{k,k+p-l,l-2k}2^{l-2k}x^{l}$
$\displaystyle=\sum_{l,k}\binom{p}{k,l,p-l-k}2^{l}x^{l+2k}=\sum_{l,k}\binom{p}{k,l,p-l-k}(2x)^{l}(x^{2})^{k}$
$\displaystyle=(x^{2}+2x+1)^{p}=(x+1)^{2p}=\sum_{l=0}^{2p}\binom{2p}{l}x^{l}.$

By identification of the terms in the sum over $l,$ this yields the equality

$\displaystyle\sum_{k=0}^{\lfloor
l/2\rfloor}\binom{p}{k,k\\!+\\!p\\!-\\!l,l\\!-\\!2k}2^{l-2k}=\binom{2p}{l}$

Plugging back into (26) and paying attention to the missing terms in the
original sum, we find, for $2\leq l\leq 2p:$

$\displaystyle\sum_{\begin{subarray}{c}k=0,j=0\\\
p-k+j=l\end{subarray}}^{k=p-2,j=p-k}\binom{p}{k}\binom{p-k}{j}2^{p-k-j}=\binom{2p}{l}-p\mathbf{1}_{l=2}.$

Plugging back in the original sum, we find

$\displaystyle(1-$ $\displaystyle(1-\beta\mu)^{2p})M_{2p}^{2p}\leq
p\mathbb{E}\|\theta-\beta\nabla\mathcal{L}(\theta)-\theta^{\star}\|^{2p-2}(\beta\|\varepsilon(\theta)\|)^{2}$
$\displaystyle\quad+\mathbb{E}\sum_{l=2}^{2p}\Big{(}\binom{2p}{l}-p\mathbf{1}_{l=2}\Big{)}\|\theta-\beta\nabla\mathcal{L}(\theta)-\theta^{\star}\|^{2p-l}(\beta\|\varepsilon(\theta)\|)^{l}$
$\displaystyle=\mathbb{E}\sum_{l=2}^{2p}\binom{2p}{l}\|\theta-\beta\nabla\mathcal{L}(\theta)-\theta^{\star}\|^{2p-l}(\beta\|\varepsilon(\theta)\|)^{l}$
(27)
$\displaystyle\stackrel{{\scriptstyle\Circled{1}}}{{\leq}}\sum_{l=2}^{2p}\binom{2p}{l}((1-\beta\mu)M_{2p-l})^{2p-l}(\beta
l\overline{K})$
$\displaystyle\stackrel{{\scriptstyle\Circled{2}}}{{\leq}}(\beta
2p\overline{K})^{2p}+\sum_{l=2}^{2p-1}\binom{2p}{l}((1-\beta\mu)(2p-l)K_{\pi})^{2p-l}(\beta
l\overline{K})^{l}$
$\displaystyle\stackrel{{\scriptstyle\Circled{3}}}{{\leq}}((2p\\!-\\!1)K_{\pi})^{2p}\Big{(}\frac{2p}{2p\\!-\\!1}\Big{)}^{2p}\Big{[}(\beta\overline{K}/K_{\pi})^{2p}+\frac{e^{\frac{1}{24p}}}{\sqrt{2\pi}}\sum_{l=2}^{2p-1}\sqrt{\frac{2p}{l(2p\\!-\\!l)}}(1\\!-\\!\beta\mu)^{2p-l}(\beta\overline{K}/K_{\pi})^{l}\Big{]}$
$\displaystyle\stackrel{{\scriptstyle\Circled{4}}}{{\leq}}((2p\\!-\\!1)K_{\pi})^{2p}\Big{(}\frac{2p}{2p\\!-\\!1}\Big{)}^{2p}\Big{[}(\beta\overline{K}/K_{\pi})^{2p}+\frac{e^{\frac{1}{24p}}}{\sqrt{4\pi}}(2p\\!-\\!2)\sqrt{\frac{2p}{2p\\!-\\!2}}(1\\!-\\!\beta\mu)^{2p-2}(\beta\overline{K}/K_{\pi})^{2}\Big{]},$

where $\Circled{1}$ uses (24) and Lemma 5, $\Circled{2}$ uses our induction
hypothesis, $\Circled{3}$ uses Stirling’s approximation and $\Circled{4}$ uses
that $1-\beta\mu>\beta\overline{K}/K_{\pi}.$ We now use the following
inequalities for $p\geq 2:$

$\displaystyle
1-(1-\beta\mu)^{2p}=\beta\mu\sum_{i=0}^{2p-1}(1-\beta\mu)^{i}\geq
2p\beta\mu(1-\beta\mu)^{2p-1}\quad\quad\beta\leq 1/(2\mu)$
$\displaystyle\Big{(}\frac{2p}{2p-1}\Big{)}^{2p-1}=\Big{(}1+\frac{1}{2p-1}\Big{)}^{2p-1}\leq
e\quad\text{and}\quad\frac{\sqrt{2p(2p-2)}}{2p-1}\leq 1,$

in addition to the choice $K_{\pi}=C\overline{K}\sqrt{\beta\mu}$ with $C=2$ to
find

$\displaystyle M_{2p}^{2p}$
$\displaystyle\leq((2p\\!-\\!1)K_{\pi})^{2p}\Big{(}\frac{e}{1-\beta\mu}\Big{)}\Big{[}\Big{(}\frac{\sqrt{\beta\mu}}{1-\beta\mu}\Big{)}^{2p-2}\frac{C^{-2p}}{2p-1}+\frac{e^{\frac{1}{48}}C^{-2}}{\sqrt{4\pi}}\Big{]}$
$\displaystyle\leq((2p\\!-\\!1)K_{\pi})^{2p}(2e)\Big{[}\frac{1}{6}\Big{(}\frac{1}{\sqrt{2}}\Big{)}^{2p}+\frac{e^{\frac{1}{48}}}{8\sqrt{\pi}}\Big{]}\leq((2p\\!-\\!1)K_{\pi})^{2p}$

which yields the desired bound (25) for $M_{2p}$ as well as $M_{2p-1}$ through
$M_{2p-1}\leq M_{2p}.$ This finishes the induction.

###  B.4 Proof of Lemma 1

Without loss of generality, we consider moments of even order. For $j\geq 1,$
denoting $M_{2j}^{2j}=\mathbb{E}\|\theta-\theta^{\star}\|^{2j}$ and starting
from Equation (27) which was obtained in the proof of Proposition 2 and using
Lemma 5 and our assumption on $\|\varepsilon(\theta)\|$ yields

$\displaystyle\big{(}1-(1-$
$\displaystyle\beta\mu)^{2j}\big{)}M_{2j}^{2j}\leq\mathbb{E}\sum_{l=2}^{2j}\binom{2j}{l}\|\theta-\beta\nabla\mathcal{L}(\theta)-\theta^{\star}\|^{2j-l}(\beta\|\varepsilon(\theta)\|)^{l}$
$\displaystyle\leq\mathbb{E}\sum_{l=2}^{2j}\binom{2j}{l}\big{(}(1-\beta\mu)\|\theta-\theta^{\star}\|\big{)}^{2j-l}\beta^{l}\big{(}K\|\theta-\theta^{\star}\|+\underline{K}\big{)}^{l}$
$\displaystyle\leq\mathbb{E}\sum_{l=2}^{2j}\binom{2j}{l}\big{(}(1-\beta\mu)\|\theta-\theta^{\star}\|\big{)}^{2j-l}\beta^{l}\sum_{k=0}^{l}\binom{l}{k}\big{(}K\|\theta-\theta^{\star}\|\big{)}^{l-k}\underline{K}^{k}$
$\displaystyle\leq\sum_{l=2}^{2j}\binom{2j}{l}(1-\beta\mu)^{2j-l}\beta^{l}\Big{(}K^{l}M_{2j}^{2j}+\sum_{k=1}^{l}\binom{l}{k}K^{l-k}\underline{K}^{k}M_{2j-k}^{2p-k}\Big{)}.$

By sorting out the factors of $M_{2j}^{2j}$ and rearranging the terms, we find

$\displaystyle\Big{(}1-(1-\beta\mu)^{2j}$
$\displaystyle-\sum_{l=2}^{2j}\binom{2j}{l}(1-\beta\mu)^{2j-l}(\beta
K)^{l}\Big{)}M_{2j}^{2j}\leq$
$\displaystyle\sum_{l=2}^{2j}\binom{2j}{l}(1-\beta\mu)^{2j-l}\beta^{l}\sum_{k=1}^{l}\binom{l}{k}K^{l-k}\underline{K}^{k}M_{2j-k}^{2j-k}.$

Assuming that $M_{i}<\infty$ for $i<2j,$ the above inequality would allow us
to show that $M_{2j}<\infty$ as well provided that the factor of $M_{2j}^{2j}$
on the LHS is positive. We now use the inequalities $0\leq(1-x)^{k}-(1-kx)\leq
k(k-1)x^{2}/2$ valid for $x\geq 0$ and $k\in\mathbb{N}^{*}$ to find

$\displaystyle 1-(1-\beta\mu)^{2j}$
$\displaystyle-\sum_{l=2}^{2j}\binom{2j}{l}(1-\beta\mu)^{2j-l}(\beta
K)^{l}=1-(1-\beta(\mu-K))^{2j}+2j\beta K(1-\beta\mu)^{2j-1}$
$\displaystyle\geq 2j\beta\mu-2j(2j-1)\beta^{2}((\mu-K)^{2}/2+\mu K)$
$\displaystyle=2j\beta(\mu-\beta(2j-1)(\mu^{2}+K^{2})/2)\geq 0,$

where the last inequality follows from the bound we imposed on $\beta.$

Therefore, we can deduce that $M_{2j}<\infty.$ Since a similar argument works
for $M_{i}<\infty$ with $i<2p$ with a weaker condition on $\beta,$ the result
follows.

###  B.5 Proof of Proposition 3

We now prove (a). Let $\theta,\theta^{\prime}\sim\pi_{\beta}$ be independent
and define the _gradient step_ function $g_{\beta}$ as

$g_{\beta}(\vartheta)=\vartheta-\beta\nabla\mathcal{L}(\vartheta)\quad\text{
for }\quad\vartheta\in\mathbb{R}^{d}.$

Note that, by Lemma 5, $g_{\beta}$ is $(1-\beta\mu)$-Lipschitz. Under
Assumption 4 (a) and using the invariance of $\pi_{\beta}$, we have for all
$\lambda\in\mathbb{R}:$

$\displaystyle\sup_{f\in\operatorname{Lip}(\mathbb{R}^{d})}$
$\displaystyle\mathbb{E}\exp\big{(}\lambda(f(\theta)\\!-\\!\mathbb{E}f(\theta^{\prime}))\big{)}$
$\displaystyle=\sup_{f\in\operatorname{Lip}(\mathbb{R}^{d})}\mathbb{E}\exp\big{(}\lambda(f(\theta\\!-\\!\beta
G(\theta))\\!-\\!\mathbb{E}f(\theta^{\prime}\\!-\\!\beta
G(\theta^{\prime})))\big{)}$
$\displaystyle=\sup_{f\in\operatorname{Lip}(\mathbb{R}^{d})}\mathbb{E}\exp\Big{(}\lambda\big{(}f((\theta-\beta\nabla\mathcal{L}(\theta))-\beta\varepsilon(\theta))-\mathbb{E}f((\theta^{\prime}-\beta\nabla\mathcal{L}(\theta^{\prime}))-\beta\varepsilon(\theta^{\prime}))\big{)}\Big{)}$
$\displaystyle=\sup_{f\in\operatorname{Lip}(\mathbb{R}^{d})}\mathbb{E}\exp\Big{(}\lambda\big{(}f(g_{\beta}(\theta))-\mathbb{E}f(g_{\beta}(\theta^{\prime}))\big{)}$
$\displaystyle\quad+\lambda\big{(}f(g_{\beta}(\theta)-\beta\varepsilon(\theta))-f(g_{\beta}(\theta))-\mathbb{E}[f(g_{\beta}(\theta^{\prime})-\beta\varepsilon(\theta^{\prime}))-f(g_{\beta}(\theta^{\prime}))]\big{)}\Big{)}.$

Conditioning on $\theta,$ it is clear that

$\phi\big{(}G(\theta)\big{)}:=f\big{(}g_{\beta}(\theta)-\beta(G(\theta)-\nabla\mathcal{L}(\theta))\big{)}-f\big{(}g_{\beta}(\theta)\big{)}=f\big{(}g_{\beta}(\theta)-\beta\varepsilon(\theta)\big{)}-f\big{(}g_{\beta}(\theta)\big{)}$

is a $\beta$-Lipschitz function of $G(\theta).$ In addition,
$f(g_{\beta}(\theta))$ is a $(1-\beta\mu)$-Lipschitz function of $\theta,$
therefore by reparametrizing the space of Lipschitz functions and using
Jensen’s inequality, we find

$\displaystyle\sup_{f\in\operatorname{Lip}(\mathbb{R}^{d})}$
$\displaystyle\mathbb{E}\exp\big{(}\lambda(f(\theta)\\!-\\!\mathbb{E}f(\theta^{\prime}))\big{)}$
$\displaystyle\leq\sup_{f\in\operatorname{Lip}(\mathbb{R}^{d})}\mathbb{E}\exp\Big{(}\lambda\big{(}f(g_{\beta}(\theta))-\mathbb{E}f(g_{\beta}(\theta^{\prime}))\big{)}\Big{)}\exp(\lambda^{2}\beta^{2}K^{2})$
$\displaystyle\leq\sup_{f\in\operatorname{Lip}(\mathbb{R}^{d})}\mathbb{E}\exp\Big{(}\lambda(1-\beta\mu)\big{(}f(\theta)-\mathbb{E}f(\theta^{\prime})\big{)}\Big{)}\exp(\lambda^{2}\beta^{2}K^{2})$
$\displaystyle\leq\Big{(}\sup_{f\in\operatorname{Lip}(\mathbb{R}^{d})}\mathbb{E}\exp\big{(}\lambda\big{(}f(\theta)-\mathbb{E}f(\theta^{\prime})\big{)}\big{)}\Big{)}^{1-\beta\mu}\exp(\lambda^{2}\beta^{2}K^{2})$
$\displaystyle\implies$
$\displaystyle\sup_{f\in\operatorname{Lip}(\mathbb{R}^{d})}\mathbb{E}\exp\big{(}\lambda(f(\theta)\\!-\\!\mathbb{E}f(\theta^{\prime}))\big{)}\leq\exp(\lambda^{2}\beta
K/\mu)$

The proof of (b) is analogous except for the fact that the above inequalities
only hold for $|\lambda|\leq(\beta K)^{-1}$ when $f(G(\theta))$ is $K$-sub-
exponential for all $f\in\operatorname{Lip}(\mathbb{R}^{d})$. The rest of the
proof is unchanged and since $\big{(}K\sqrt{\beta/\mu}\big{)}^{-1}<(\beta
K)^{-1}$, we similarly obtain the sub-exponential property.

###  B.6 Proof of Proposition 4

Let $\theta_{1}\sim\nu_{1}$ and $\theta_{2}\sim\nu_{2}$ be random variables
such that
$\mathcal{W}_{2}^{2}(\nu_{1},\nu_{2})=\mathbb{E}\big{[}\|\theta_{1}-\theta_{2}\|^{2}\big{]}.$
Such a pair of variables exists by [90, Theorem 4.1].

We consider the set of couplings of the distributions $\nu_{1}P$ and
$\nu_{2}P$ through the random variables $G(\theta_{1})$ and $G(\theta_{2})$
such that

$\displaystyle\theta_{1}-\beta G(\theta_{1})\sim\nu_{1}P\quad$
$\displaystyle\text{and}\quad\theta_{2}-\beta G(\theta_{2})\sim\nu_{2}P.$

Recall also that by Assumption 2, for $j=1,2,$ conditionally on $\theta_{j},$
we have

$G(\theta_{j})=\nabla\mathcal{L}(\theta_{j})+\varepsilon(\theta_{j})\quad\text{with}\quad\mathbb{E}[\varepsilon(\theta_{j})|\theta_{j}]=0.$
(28)

Taking the infimum over such variables $G(\theta_{j}),$ we compute

$\displaystyle\mathcal{W}_{2}^{2}(\nu_{1}P,\nu_{2}P)$
$\displaystyle=\inf_{G(\theta_{j})}\mathbb{E}\big{\|}\theta_{1}-\beta
G(\theta_{1})-(\theta_{2}-\beta G(\theta_{2}))\big{\|}^{2}$
$\displaystyle=\inf_{G(\theta_{j})}\mathbb{E}\big{[}\big{\|}\theta_{1}-\beta\nabla\mathcal{L}(\theta_{1})-(\theta_{2}-\beta\nabla\mathcal{L}(\theta_{2}))\big{\|}^{2}$
$\displaystyle\quad-2\beta\langle\theta_{1}-\beta\nabla\mathcal{L}(\theta_{1})-(\theta_{2}-\beta\nabla\mathcal{L}(\theta_{2})),\varepsilon(\theta_{1})-\varepsilon(\theta_{2})\rangle$
$\displaystyle\quad+\beta^{2}\|\varepsilon(\theta_{1})-\varepsilon(\theta_{2})\|^{2}\big{]}$
$\displaystyle\stackrel{{\scriptstyle\Circled{1}}}{{=}}\mathbb{E}\big{[}\big{\|}\theta_{1}-\beta\nabla\mathcal{L}(\theta_{1})-(\theta_{2}-\beta\nabla\mathcal{L}(\theta_{2}))\big{\|}^{2}$
$\displaystyle\quad+\beta^{2}\inf_{G(\theta_{j})}\mathbb{E}\big{[}\|\varepsilon(\theta_{1})-\varepsilon(\theta_{2})\|^{2}|\theta_{1},\theta_{2}\big{]}\big{]}$
$\displaystyle\stackrel{{\scriptstyle\Circled{2}}}{{\leq}}\mathbb{E}\big{[}(1-\beta\mu)^{2}\|\theta_{1}-\theta_{2}\|^{2}+\beta^{2}\mathcal{W}_{2}^{2}\big{(}\mathcal{D}(\varepsilon(\theta_{1})),\mathcal{D}(\varepsilon(\theta_{2}))\big{)}\big{]}$
$\displaystyle\stackrel{{\scriptstyle\Circled{3}}}{{\leq}}\mathbb{E}\big{[}\big{(}(1-\beta\mu)^{2}+\beta^{2}L_{\mathcal{W}}\big{)}\|\theta_{1}-\theta_{2}\|^{2}\big{]}$
$\displaystyle=\big{(}(1-\beta\mu)^{2}+\beta^{2}L_{\mathcal{W}}\big{)}\mathcal{W}_{2}^{2}(\nu_{1},\nu_{2}),$

where $\Circled{1}$ is obtained by conditioning on $\theta_{1},\theta_{2}$ and
using (28), $\Circled{2}$ uses Lemma 5 and $\Circled{3}$ uses Assumption 5.

Since $\beta<\frac{2\mu}{\mu^{2}+L_{\mathcal{W}}}$ by assumption, the obtained
inequality shows that the mapping $\nu\to\nu P$ is a contraction in the space
$\mathcal{P}_{2}(\mathbb{R}^{d})$ endowed with the $\mathcal{W}_{2}$ metric
which is complete and separable by [90, Theorem 6.18]. Consequently, by
Banach’s fixed-point theorem, the previous mapping admits a unique fixed point
$\pi_{\beta}\in\mathcal{P}_{2}(\mathbb{R}^{d})$ i.e. such that
$\pi_{\beta}P=\pi_{\beta}.$ Moreover, for any initial measure
$\xi_{0}\in\mathcal{P}_{2}(\mathbb{R}^{d}),$ the sequence
$(\xi_{n})_{n\in\mathbb{N}}$ defined by $\xi_{n}=\xi_{0}P^{n}$ converges to
$\pi_{\beta}$ w.r.t. the $\mathcal{W}_{2}$ metric.

###  B.7 Proof of Corollary 1

By Theorem 1, we know that the Markov chain $(\theta_{t})_{t\geq 0}$ is
geometrically converging to the invariant distribution $\pi_{\beta}$ w.r.t.
the total-variation distance so that for any event
$\mathcal{E}\in\mathcal{B}(\mathbb{R}^{d}),$ we have

$\big{|}\mathbb{P}\big{(}\theta_{n}\in\mathcal{E}\big{)}-\mathbb{P}_{\theta\sim\pi}\big{(}\theta\in\mathcal{E}\big{)}\big{|}\leq
M\rho^{n}V(\theta_{0}).$ (29)

By combining the above inequality with the conclusions of Proposition 2 and
using Lemmas 3 and 4 in the sub-Gaussian and sub-exponential cases
respectively we obtain the claimed bounds.

###  B.8 Proof of Corollary 2

Similarly to the proof of Corollary 1, inequality (29) holds by Theorem 1 for
any event $\mathcal{E}\in\mathcal{B}(\mathbb{R}^{d}).$ For
$\theta\sim\pi_{\beta}$ stationary, let $f(\theta)=\|\theta-\theta^{\star}\|$
and denote $\Delta=f(\theta)-\mathbb{E}f(\theta).$ Using Proposition 3 (a) and
Chernoff’s method for $t>0$ and $\lambda>0,$ we have

$\mathbb{P}(\Delta>t)=\mathbb{P}(e^{\lambda\Delta}>e^{\lambda
t})\leq\mathbb{E}\exp\big{(}\lambda\Delta-\lambda
t\big{)}\leq\exp\big{(}\lambda^{2}K^{2}\beta/\mu-\lambda t\big{)}.$ (30)

After minimizing over $\lambda,$ we get for $\delta>0,$ with probability at
least $1-\delta,$ the following inequality holds:

$\Delta\leq 2K\sqrt{\beta\log(1/\delta)/\mu}.$ (31)

Additionally, by Proposition 1 (b) we have

$\mathbb{E}f(\theta)=\mathbb{E}\|\theta-\theta^{\star}\|\leq\sqrt{\frac{\beta\sigma^{2}}{2\mu-\beta(\mu^{2}+L_{\sigma})}}.$
(32)

Taking $\beta$ such that $\beta\leq\frac{\mu}{\mu^{2}+L_{\sigma}}$ and
combining (31) with (32) yields (15).

To obtain (16), we proceed similarly using Proposition 3 (b) this time.
Applying the constraint $|\lambda|\leq(\beta K)^{-1}$ (see proof of
Proposition 3) into the optimization of (30) yields

$\mathbb{P}(\Delta>t)\leq\begin{cases}\exp\Big{(}\frac{-t^{2}}{4\beta
K^{2}/\mu}\Big{)}&\text{if}\quad t\leq 2K/\mu\\\ \exp\Big{(}\frac{-t}{2\beta
K}\Big{)}&\text{otherwise.}\end{cases}$ (33)

Expressing $t$ in terms of the failure probability $\delta$ and combining with
(32) as before finishes the proof.

######  Lemma 6.

Grant Assumption 1, 2, 5 and 6. Let the Markov chain $(\theta_{t})_{t\geq}$ be
initialized with $\theta_{0}\sim\nu$ and $\beta$ be chosen as in Proposition
4. The sequence of SGD iterates $\theta_{0},\dots,\theta_{n}$ satisfies for
$0\leq i,j\leq n:$

$\mathbb{E}\langle\theta_{i}-\theta^{\star},\theta_{j}-\theta^{\star}\rangle\leq
2(1-\beta\mu)^{|i-j|}\big{(}\big{(}(1-\beta\mu)+\beta^{2}L_{\mathcal{W}}\big{)}^{i}\mathcal{W}_{2}^{2}(\nu,\pi)+\operatorname{Var}_{\pi}(\theta)\big{)}.$

###### Proof.

We assume without loss of generality that $i\leq j.$ Since the gradient is
linear it commutes with the expectation. Therefore, by conditioning over
$\theta_{j-1}$ and later over $\theta_{j-1}$ we find

$\displaystyle\mathbb{E}\langle\theta_{j}-\theta^{\star},\theta_{i}-\theta^{\star}\rangle$
$\displaystyle=\mathbb{E}\langle\theta_{j-1}-\beta
G(\theta_{j-1})-\theta^{\star},\theta_{i}-\theta^{\star}\rangle$
$\displaystyle=\mathbb{E}\langle\theta_{j-1}-\beta\nabla\mathcal{L}(\theta_{j-1})-\theta^{\star},\theta_{i}-\theta^{\star}\rangle$
$\displaystyle=\mathbb{E}\langle\theta_{j-2}-\beta
G(\theta_{j-2})-\beta\nabla\mathcal{L}(\theta_{j-2}-\beta
G(\theta_{j-2}))-\theta^{\star},\theta_{i}-\theta^{\star}\rangle$
$\displaystyle=\mathbb{E}\langle\theta_{j-2}-\beta\nabla\mathcal{L}(\theta_{j-2})-\beta\nabla\mathcal{L}(\theta_{j-2}-\beta\nabla\mathcal{L}(\theta_{j-2}))-\theta^{\star},\theta_{i}-\theta^{\star}\rangle.$

It is clear that the previous steps can be repeated for the remaining indices
$j-3,j-4,\dots$ down to $i$ at which point the following identity is reached

$\mathbb{E}\langle\theta_{j}-\theta^{\star},\theta_{i}-\theta^{\star}\rangle=\mathbb{E}\langle\widecheck{\theta}_{j}-\theta^{\star},\theta_{i}-\theta^{\star}\rangle,$

where $\widecheck{\theta}_{j}$ is recursively defined by
$\widecheck{\theta}_{i}=\theta_{i}$ and
$\widecheck{\theta}_{k}=\widecheck{\theta}_{k-1}-\beta\nabla\mathcal{L}(\widecheck{\theta}_{k-1})$
for $i<k\leq j.$

Using Cauchy-Schwarz and iterating the inequality
$\|\widecheck{\theta}_{k}-\theta^{\star}\|\leq(1-\beta\mu)\|\widecheck{\theta}_{k-1}-\theta^{\star}\|$
yields that

$\mathbb{E}\langle\theta_{j}-\theta^{\star},\theta_{i}-\theta^{\star}\rangle\leq(1-\beta\mu)^{j-i}\mathbb{E}\|\theta_{i}-\theta^{\star}\|^{2}.$

Now, by [90, Theorem 4.1], there exists a random variable
$\widetilde{\theta}\sim\pi_{\beta}$ such that the coupling
$(\theta_{i},\widetilde{\theta})$ satisfies

$\displaystyle\mathbb{E}\|\theta_{i}-\widetilde{\theta}\big{\|}^{2}$
$\displaystyle=\mathcal{W}_{2}^{2}(\mathcal{D}(\theta_{i}),\pi)=\mathcal{W}_{2}^{2}(\nu
P^{i},\pi)\leq\big{(}(1-\beta\mu)+\beta^{2}L_{\mathcal{W}}\big{)}^{i}\mathcal{W}_{2}^{2}(\nu,\pi),$

where the inequality comes from Proposition 4. It then only remains to write

$\displaystyle\mathbb{E}\|\theta_{i}-\theta^{\star}\|^{2}$ $\displaystyle\leq
2\big{(}\mathbb{E}\|\theta_{i}-\widetilde{\theta}\|^{2}+\mathbb{E}\|\widetilde{\theta}-\theta^{\star}\|^{2}\big{)}\leq
2\big{(}(1-\beta\mu)+\beta^{2}L_{\mathcal{W}}\big{)}^{i}\mathcal{W}_{2}^{2}(\nu,\pi)+2\operatorname{Var}_{\pi}(\theta),$

which implies the result. The case $i\geq j$ is handled similarly. ∎

######  Lemma 7.

Let $A\in\mathbb{R}^{n\times n}$ be a matrix with positive entries such that
there exists $C>0$ and $0<\alpha<1$ such that

$A_{ij}\leq C\alpha^{|j-i|}\quad\text{ for }\quad 1\leq i,j\leq n,$

then we have

$\sum_{i,j}A_{ij}\leq
C\Big{(}n+\frac{2\alpha}{1-\alpha}\Big{(}n-\Big{(}\frac{1-\alpha^{n}}{1-\alpha}\Big{)}\Big{)}\Big{)}.$

###### Proof.

Straightforward computations yield

$\displaystyle\sum_{i,j}A_{ij}$
$\displaystyle=\sum_{i=1}^{n}A_{ii}+2\sum_{i<j}A_{ij}\leq
nC+2\sum_{i<j}A_{ij},$

and we have

$\displaystyle\sum_{i=1}^{n}\sum_{j=i+1}^{n}A_{ij}$ $\displaystyle\leq
C\sum_{i=1}^{n}\sum_{j=i+1}^{n}\alpha^{j-i}=C\alpha\sum_{i=1}^{n}\frac{1-\alpha^{n-i}}{1-\alpha}$
$\displaystyle=\frac{C\alpha}{1-\alpha}\Big{(}n-\sum_{i=1}^{n}\alpha^{n-i}\Big{)}=\frac{C\alpha}{1-\alpha}\Big{(}n-\frac{1-\alpha^{n}}{1-\alpha}\Big{)}.$

∎

###  B.9 Proof of Theorem 2

We introduce the notations
$\theta_{[i]}=(\theta_{1},\theta_{2},\dots,\theta_{i})$ and
$\theta_{[k,l]}=(\theta_{k},\theta_{k+1},\dots,\theta_{l})$ and define for
$1\leq i\leq n$ the variables

$M^{(i)}:=\mathbb{E}\big{[}f(\vec{\theta})|\theta_{[i]}\big{]}-\mathbb{E}\big{[}f(\vec{\theta})|\theta_{[i-1]}\big{]}\quad\text{so
that}\quad f(\vec{\theta})-\mathbb{E}[f(\vec{\theta})]=\sum_{i=1}^{n}M^{(i)}.$

Notice that, if we condition on $\theta_{[i-1]}$ then $M^{(i)}$ only depends
on $\theta_{i}$. We consider $M^{(i)}$ as a function of $\theta_{i}$ and
compute its Lipschitz constant. We temporarily consider $\theta_{i}$ and
$\theta_{i}^{\prime}$ as two fixed deterministic vectors of $\mathbb{R}^{d}$
and $\theta_{i+1},\theta_{i+2},\dots$ and
$\theta^{\prime}_{i+1},\theta^{\prime}_{i+2},\dots$ are the SGD trajectories
obtained from them i.e. for $j>i:$

$\theta_{j}=\theta_{j-1}-\beta G(\theta_{j-1})\quad\text{ and
}\quad\theta^{\prime}_{j}=\theta^{\prime}_{j-1}-\beta
G(\theta^{\prime}_{j-1}).$

In the following, we use the Lipschitz property of $f$ and the Kantorovich-
Rubinstein dual representation of the $\mathcal{W}_{1}$ metric

$\mathcal{W}_{1}(\nu_{1},\nu_{2})=\sup_{h\in\operatorname{Lip}(\mathbb{R}^{d})}\int
hd\nu_{1}-\int hd\nu_{2},$

in order to find

$\displaystyle\big{|}M^{(i)}(\theta_{i})-M^{(i)}(\theta_{i}^{\prime})\big{|}$
$\displaystyle=\bigg{|}\mathbb{E}\big{[}f(\vec{\theta})|\theta_{[i]}\big{]}-\mathbb{E}\big{[}f(\theta_{[i-1]},\theta_{[i,n]}^{\prime})|\theta_{i}^{\prime},\theta_{[i-1]}\big{]}\bigg{|}$
$\displaystyle=\bigg{|}\sum_{j=i}^{n}\mathbb{E}\Big{[}f\big{(}\theta_{[i-1]},\theta_{[i,j-1]}^{\prime},\theta_{[j,n]}\big{)}-f\big{(}\theta_{[i-1]},\theta_{[i,j]}^{\prime},\theta_{[j+1,n]}\big{)}|\theta_{[i-1]}\Big{]}\bigg{|}$
$\displaystyle\leq\sum_{j=i}^{n}\mathcal{W}_{1}\big{(}\mathcal{D}(\theta_{j}),\mathcal{D}(\theta_{j}^{\prime})\big{)}.$

Using Proposition 4 we have

$\displaystyle\mathcal{W}_{1}\big{(}\mathcal{D}(\theta_{j}),\mathcal{D}(\theta_{j}^{\prime})\big{)}$
$\displaystyle\leq\mathcal{W}_{2}\big{(}\mathcal{D}(\theta_{j}),\mathcal{D}(\theta_{j}^{\prime})\big{)}=\mathcal{W}_{2}\big{(}\mathcal{D}(\theta_{j-1})P,\mathcal{D}(\theta^{\prime}_{j-1})P\big{)}$
$\displaystyle\leq\underbrace{\sqrt{(1-\beta\mu)^{2}+\beta^{2}L_{\mathcal{W}}}}_{=:\alpha_{\mathcal{W}}(\beta,\mu)}\mathcal{W}_{2}\big{(}\mathcal{D}(\theta_{j-1}),\mathcal{D}(\theta^{\prime}_{j-1})\big{)}$
$\displaystyle\leq\dots$
$\displaystyle\leq\alpha_{\mathcal{W}}(\beta,\mu)^{j-i}\mathcal{W}_{2}\big{(}\mathcal{D}(\theta_{i}),\mathcal{D}(\theta^{\prime}_{i})\big{)}=\alpha_{\mathcal{W}}(\beta,\mu)^{j-i}\big{\|}\theta_{i}-\theta^{\prime}_{i}\big{\|},$

where the last equality follows from $\theta_{i}$ and $\theta_{i}^{\prime}$
being deterministic. Provided that
$\beta<\frac{2\mu}{\mu^{2}+L_{\mathcal{W}}}$ we have
$\alpha_{\mathcal{W}}(\beta,\mu)<1$ so that
$\mathcal{W}_{1}(\theta_{j},\theta_{j}^{\prime})\leq\alpha_{\mathcal{W}}(\beta,\mu)^{j-i}\|\theta_{i}-\theta_{i}^{\prime}\|$
for $i\leq j\leq n.$ By summing over $j,$ we find that the $M^{(i)}$s are
$(1-\alpha_{\mathcal{W}}(\beta,\mu))^{-1}$-Lipschitz

$\big{|}M^{(i)}(\theta_{i})-M^{(i)}(\theta_{i}^{\prime})\big{|}\leq\frac{\|\theta_{i}-\theta_{i}^{\prime}\|}{1-\alpha_{\mathcal{W}}(\beta,\mu)}.$

In what follows we denote
$\mathbb{E}_{k}[\cdot]=\mathbb{E}[\cdot|\theta_{[k]}]$ to lighten notation and
let $C_{\mathcal{W}}:=(1-\alpha_{\mathcal{W}}(\beta,\mu))^{-1}.$ Let
$\lambda\in\mathbb{R},$ by conditioning on $\theta_{[n-1]},$ we have

$\displaystyle\mathbb{E}\exp\big{(}\lambda(f(\vec{\theta})-\mathbb{E}f(\vec{\theta}))\big{)}$
$\displaystyle=\mathbb{E}\exp\Big{(}\lambda\sum_{i=1}^{n}M^{(i)}\Big{)}=\mathbb{E}\Big{[}\mathbb{E}\Big{[}\exp\Big{(}\lambda\sum_{i=1}^{n}M^{(i)}\Big{)}\big{|}\theta_{[n-1]}\Big{]}\Big{]}$
$\displaystyle=\mathbb{E}\Big{[}\exp\Big{(}\lambda\sum_{i=1}^{n-1}M^{(i)}\Big{)}\mathbb{E}\big{[}\exp(\lambda
M^{(n)})|\theta_{[n-1]}\big{]}\Big{]}.$

Recall that conditionally on $\theta_{[n-1]},$ we have that $M^{(n)}$ is a
function of $\theta_{n}=\theta_{n-1}-\beta G(\theta_{n-1})$ so that $M^{(n)}$
is a $\beta C_{\mathcal{W}}$-Lipschitz function of $G(\theta_{n-1})$ which
satisfies Assumption 4 (a) and thus

$\mathbb{E}\big{[}\exp(\lambda
M^{(n)})|\theta_{[n-1]}\big{]}\leq\exp\big{(}\lambda^{2}\beta^{2}C_{\mathcal{W}}^{2}K^{2}\big{)}.$

By repeating this argument $n-1$ times, we arrive at

$\displaystyle\mathbb{E}\exp\big{(}\lambda(f(\vec{\theta})-\mathbb{E}f(\vec{\theta}))\big{)}$
$\displaystyle\leq\mathbb{E}\big{[}\exp\big{(}\lambda
M^{(1)}\big{)}\big{]}\exp\big{(}(n-1)\lambda^{2}\beta^{2}C_{\mathcal{W}}^{2}K^{2}\big{)}$
$\displaystyle\leq\exp\big{(}\lambda^{2}C_{\mathcal{W}}^{2}K^{2}\beta/\mu+(n-1)\lambda^{2}\beta^{2}C_{\mathcal{W}}^{2}K^{2}\big{)},$

where the last inequality uses that $\theta_{1}\sim\pi_{\beta}$ which is
$K\sqrt{\beta/\mu}$-sub-Gaussian by Proposition 3 (a).

The proof in the sub-exponential case is completely analogous using Assumption
4 (b) and the result of Proposition 3 (b) with the main difference that the
obtained inequalities only hold for
$|\lambda|\leq\big{(}C_{\mathcal{W}}K\sqrt{\beta/\mu}\big{)}^{-1}\wedge\big{(}\beta
C_{\mathcal{W}}K\big{)}^{-1}=\big{(}C_{\mathcal{W}}K\sqrt{\beta/\mu}\big{)}^{-1}$
because $\beta<\mu^{-1}.$

###  B.10 Proof of Proposition 6

For $j\geq 0,$ we introduce the notation

$\Delta_{j}:=\Big{\|}\sum_{t=j+1}^{j+n}\theta_{t}-n\theta^{\star}\Big{\|}-\mathbb{E}\Big{\|}\sum_{t=j+1}^{j+n}\theta_{t}-n\theta^{\star}\Big{\|}.$

We are interested in obtaining a high probability bound on the quantity
$\Delta_{n_{0}}.$ We write $\mathbb{E}_{\nu}$ for the expectation when the
Markov chain is started with distribution $\nu$

$\displaystyle\mathbb{E}_{\nu}\big{[}\exp\big{(}\lambda\Delta_{n_{0}}\big{)}\big{]}$
$\displaystyle=\mathbb{E}_{\nu
P^{n_{0}}}\big{[}\exp\big{(}\lambda\Delta_{1}\big{)}\big{]}=\mathbb{E}_{\pi}\Big{[}\frac{d(\nu
P^{n_{0}})}{d\pi}\exp\big{(}\lambda\Delta_{1}\big{)}\Big{]}$
$\displaystyle\leq\Big{\|}\frac{d(\nu
P^{n_{0}})}{d\pi}\Big{\|}_{\pi,\infty}\mathbb{E}_{\pi}\big{[}\exp\big{(}\lambda\Delta_{1}\big{)}\big{]}.$

The last expectation can be bounded using Theorem 2. As for the factor coming
from the measure change, we write

$\Big{\|}\frac{d(\nu
P^{n_{0}})}{d\pi}\Big{\|}_{\pi,\infty}\leq\Big{\|}\frac{d(\nu
P^{n_{0}}-\pi)}{d\pi}\Big{\|}_{\pi,\infty}+1.$

For any function $F:\mathbb{R}^{d}\to\mathbb{R},$ we define the norm
$\|F\|_{V}=\sup_{\vartheta\in\mathbb{R}^{d}}\frac{|F(\vartheta)|}{V(\vartheta)}$
and its induced operator norm
${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|Q\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{V}=\sup_{F}\frac{\|QF\|_{V}}{\|F\|_{V}},$
where $V$ is the function defined in Section B. We also denote $\odot$ the
pointwise product between functions.

$\displaystyle\Big{\|}\frac{d(\nu P^{n_{0}}-\pi)}{d\pi}\Big{\|}_{\pi,\infty}$
$\displaystyle=\Big{\|}\frac{d(\nu(P^{n_{0}}-\mathbf{1}\otimes\pi))}{d\pi}\Big{\|}_{\pi,\infty}=\Big{\|}(P^{n_{0}}-\mathbf{1}\otimes\pi)^{*}\frac{d\nu}{d\pi}\Big{\|}_{\pi,\infty}$
$\displaystyle=\Big{\|}(P^{n_{0}}-\mathbf{1}\otimes\pi)^{*}\frac{d\nu}{d\pi}\odot
V\odot\frac{1}{V}\Big{\|}_{\pi,\infty}$
$\displaystyle\leq\Big{\|}(P^{n_{0}}-\mathbf{1}\otimes\pi)^{*}\frac{d\nu}{d\pi}\odot
V\Big{\|}_{V}$
$\displaystyle\leq{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|(P^{n_{0}}-\mathbf{1}\otimes\pi)^{*}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{V}\Big{\|}\frac{d\nu}{d\pi}\odot
V\Big{\|}_{V}={\left|\kern-1.07639pt\left|\kern-1.07639pt\left|P^{n_{0}}-\mathbf{1}\otimes\pi\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{V}\Big{\|}\frac{d\nu}{d\pi}\Big{\|}_{\infty}.$

The outer product $\mathbf{1}\otimes\pi_{\beta}$ denotes the kernel such that
$\mathbf{1}\otimes\pi(\vartheta,A)=\pi(A)$ for all $\vartheta$ and
$A\in\mathcal{B}(\mathbb{R}^{d}).$ By the proof of Theorem 1 and [43,
Proposition 1.1] (see also Equation (4)) the kernel $P$ has a spectral gap in
the Banach space $L_{\infty}^{V}$ of functions with finite norm
$\|\cdot\|_{V}$ and, therefore, there exist $\rho<1$ and $M<\infty$ such that

${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|P^{n_{0}}-\mathbf{1}\otimes\pi\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{V}\leq
M\rho^{n_{0}},$

which leads to

$\Big{\|}\frac{d(\nu P^{n_{0}})}{d\pi}\Big{\|}_{\pi,\infty}\leq
1+M\rho^{n_{0}}\Big{\|}\frac{d\nu}{d\pi}\Big{\|}_{\infty}=\Upsilon(\nu,n_{0}).$

Using Theorem 2 in the sub-Gaussian case, denoting
$\breve{K}=KC_{\mathcal{W}}\sqrt{\beta/\mu+(n-1)\beta^{2}},$ we find

$\mathbb{E}_{\nu}\Big{[}\exp\big{(}\lambda\Delta_{n_{0}}\big{)}\Big{]}\leq\Upsilon(\nu,n_{0})\exp(\lambda^{2}\breve{K}^{2}).$

Using Chernoff’s method for a random variable $X\in\Psi_{2}(\breve{K})$ and
$t>0$ and $\lambda>0,$ we have

$\mathbb{P}_{\nu}(\Delta_{n_{0}}>t)=\mathbb{P}_{\nu}(e^{\lambda\Delta_{n_{0}}}>e^{\lambda
t})\leq\mathbb{E}_{\nu}\exp\big{(}\lambda\Delta_{n_{0}}-\lambda
t\big{)}\leq\Upsilon(\nu,n_{0})\exp\big{(}\lambda^{2}\breve{K}^{2}-\lambda
t\big{)}.$

After minimizing over $\lambda,$ we get for $\delta>0,$ with probability at
least $1-\Upsilon(\nu,n_{0})\delta,$ the following inequality holds

$\Delta_{n_{0}}\leq 2\breve{K}\sqrt{\log(1/\delta)}.$ (34)

In the sub-exponential case (under Assumption 4 (b)), taking the constraint
$|\lambda|\leq\big{(}C_{\mathcal{W}}K\sqrt{\beta/\mu}\big{)}^{-1}$ into
account (see the proof of Theorem 2), we get that

$\mathbb{P}(\Delta_{n_{0}}>t)\leq\begin{cases}\Upsilon(\nu,n_{0})\exp\Big{(}\frac{-t^{2}}{4\breve{K}^{2}}\Big{)}&\text{if}\quad
t\leq\frac{2\breve{K}^{2}}{C_{\mathcal{W}}K\sqrt{\beta/\mu}}\\\
\Upsilon(\nu,n_{0})\exp\Big{(}\frac{-t}{2C_{\mathcal{W}}K\sqrt{\beta/\mu}}\Big{)}&\text{otherwise.}\end{cases}$

So that with probability at least $1-\Upsilon(\nu,n_{0})\delta:$

$\Delta_{n_{0}}\leq\max\big{(}2\breve{K}\sqrt{\log(1/\delta)},2C_{\mathcal{W}}K\sqrt{\beta/\mu}\log(1/\delta)\big{)}.$
(35)

It then only remains to bound the expectation
$\mathbb{E}\Big{\|}\sum_{t=n_{0}+1}^{n_{0}+n}\\!\theta_{t}-n\theta^{\star}\Big{\|},$
which can be done as follows

$\displaystyle\Big{(}\mathbb{E}\Big{\|}\sum_{t=n_{0}+1}^{n_{0}+n}\theta_{t}-n\theta^{\star}\Big{\|}\Big{)}^{2}$
$\displaystyle\leq\mathbb{E}\Big{\|}\sum_{t=n_{0}+1}^{n_{0}+n}(\theta_{t}-\theta^{\star})\Big{\|}^{2}=\sum_{i=n_{0}+1}^{n_{0}+n}\sum_{j=n_{0}+1}^{n_{0}+n}\mathbb{E}\langle\theta_{i}-\theta^{\star},\theta_{j}-\theta^{\star}\rangle.$

Using Lemmas 6 and 7 we find that

$\displaystyle\Big{(}\mathbb{E}\Big{\|}\sum_{t=n_{0}+1}^{n_{0}+n}\theta_{t}-\theta^{\star}\Big{\|}\Big{)}^{2}\leq
2n\frac{1+\alpha}{1-\alpha}\Big{(}\alpha_{\mathcal{W}}^{n_{0}}\mathcal{W}_{2}^{2}(\nu,\pi)+\operatorname{Var}_{\pi}(\theta)\Big{)},$

where $\alpha=1-\beta\mu$ and
$\alpha_{\mathcal{W}}=\sqrt{\alpha^{2}+\beta^{2}L_{\mathcal{W}}}.$ Moreover,
since $\beta<\frac{\mu}{\mu^{2}+L_{\sigma}},$ by Proposition 1 we have

$\operatorname{Var}_{\pi}(\theta)\leq\frac{\beta\sigma^{2}}{\mu}.$

Plugging into Inequalities (34) and (35) and dividing by $n$ finishes the
proof.

###  B.11 Proof of Lemma 2

Denote $\Xi_{t}^{(N)}=\frac{1}{N}\sum_{i=1}^{N}\Xi_{tN+i}$ and
$\xi_{t}^{(N)}=\frac{1}{N}\sum_{i=1}^{N}\xi_{tN+i}.$ By Lemma 8, we have the
following concentration inequalities for all $0\leq t<T:$

$\displaystyle\mathbb{P}\Big{(}\|\Xi_{t}^{(N)}\|_{2}>3K_{\Xi}\Big{(}\frac{\log(4T/\delta)+3d}{N}\vee\sqrt{\frac{\log(4T/\delta)+3d}{N}}\Big{)}\Big{)}\leq\delta/(2T)$
$\displaystyle\mathbb{P}\Big{(}\|\xi_{t}^{(N)}\|>4K_{\xi}\Big{(}\frac{\log(4T/\delta)+2d}{N}\vee\sqrt{\frac{\log(4T/\delta)+2d}{N}}\Big{)}\Big{)}\leq\delta/(2T).$

We will show by induction over $0\leq t\leq T$ that we have with probability
at least $1-t\delta/T$ that

$\max_{0\leq s\leq t}\|\theta_{s}-\theta^{\star}\|\leq C.$ (36)

The case $t=0$ holds by assumption. Further, assuming the property at rank $t$
and conditioning on $\theta_{t}$ we have with probability at least
$1-\delta/T:$

$\displaystyle\|\theta_{t+1}-\theta^{\star}\|^{2}$
$\displaystyle=\|\theta_{t}-\beta\nabla\mathcal{L}(\theta_{t})-\beta(\Xi_{t}^{(N)}(\theta_{t}-\theta^{\star})+\xi_{t}^{(N)})-\theta^{\star}\|^{2}$
$\displaystyle=\|\theta_{t}-\beta\nabla\mathcal{L}(\theta_{t})-\theta^{\star}\|^{2}-2\beta\langle\theta_{t}-\beta\nabla\mathcal{L}(\theta_{t})-\theta^{\star},\Xi_{t}^{(N)}(\theta_{t}-\theta^{\star})+\xi_{t}^{(N)}\rangle$
$\displaystyle\quad+\beta^{2}\|\Xi_{t}^{(N)}(\theta_{t}-\theta^{\star})+\xi_{t}^{(N)}\|^{2}$
$\displaystyle\stackrel{{\scriptstyle\Circled{1}}}{{\leq}}(1-\beta\mu)^{2}\|\theta_{t}-\theta^{\star}\|^{2}+2\beta(1-\beta\mu)\|\theta_{t}-\theta^{\star}\|(\|\Xi_{t}^{(N)}(\theta_{t}-\theta^{\star})\|+\|\xi_{t}^{(N)}\|)$
$\displaystyle\quad+2\beta^{2}\|\Xi_{t}^{(N)}(\theta_{t}-\theta^{\star})\|^{2}+2\beta^{2}\|\xi_{t}^{(N)}\|^{2}$
$\displaystyle\leq\big{[}(1-\beta\mu)^{2}+2\beta(1-\beta\mu)\|\Xi_{t}^{(N)}\|_{2}+2\beta^{2}\|\Xi_{t}^{(N)}\|_{2}^{2}\big{]}\|\theta_{t}-\theta^{\star}\|^{2}$
$\displaystyle\quad+2\beta(1-\beta\mu)\|\theta_{t}-\theta^{\star}\|\|\xi_{t}^{(N)}\|+2\beta^{2}\|\xi_{t}^{(N)}\|^{2}$
$\displaystyle\stackrel{{\scriptstyle\Circled{2}}}{{\leq}}\big{[}(1-\beta\mu)^{2}(1+\epsilon)+2\beta(1-\beta\mu)\|\Xi_{t}^{(N)}\|_{2}+2\beta^{2}\|\Xi_{t}^{(N)}\|_{2}^{2}\big{]}\|\theta_{t}-\theta^{\star}\|^{2}$
$\displaystyle\quad+\beta^{2}(2+1/\epsilon)\|\xi_{t}^{(N)}\|^{2}$
$\displaystyle\stackrel{{\scriptstyle\Circled{3}}}{{\leq}}\big{[}(1-\beta\mu)+2\beta\|\Xi_{t}^{(N)}\|_{2}+2\beta^{2}\|\Xi_{t}^{(N)}\|_{2}^{2}\big{]}\|\theta_{t}-\theta^{\star}\|^{2}+3\frac{\beta}{\mu}\|\xi_{t}^{(N)}\|^{2}$
$\displaystyle\stackrel{{\scriptstyle\Circled{4}}}{{\leq}}\big{[}(1-\beta\mu)+\beta\mu/3+\beta\mu/3\big{]}C^{2}+\beta\mu
C^{2}/3\leq C^{2},$

where $\Circled{1}$ uses Lemma 5 and the Cauchy-Schwarz inequality,
$\Circled{2}$ uses the inequality $2ab\leq a^{2}\epsilon+b^{2}/\epsilon$ valid
for all $\epsilon>0$ and $\Circled{3}$ sets the choice $\epsilon=\beta\mu$ and
uses that $\beta\leq 1/\mu.$ Finally $\Circled{4}$ uses the high probability
bounds stated above and the conditions on $N$ and $\beta$.

Using a union bound argument, we obtain (36) for $t+1$ with probability at
least $1-(t+1)\delta/T.$ The induction argument is completed and implies the
result for $t=T.$

######  Lemma 8.

Let $\Xi_{1},\dots,\Xi_{N}\in\mathbb{R}^{d\times d}$ be random matrices and
$\xi_{1},\dots,\xi_{N}\in\mathbb{R}^{d}$ random vectors as in Lemma 2. Then
denoting $\overline{\Xi}=\frac{1}{N}\sum_{i=1}^{N}\Xi_{i}$ and
$\overline{\xi}=\frac{1}{N}\sum_{i=1}^{N}\xi_{i},$ we have the high
probability bounds

$\displaystyle\mathbb{P}\Big{(}\|\overline{\Xi}\|_{2}>3K_{\Xi}\phi\Big{(}\frac{\log(2/\delta)+3d}{N}\Big{)}\Big{)}\leq\delta,$
(37)
$\displaystyle\mathbb{P}\Big{(}\|\overline{\xi}\|>4K_{\xi}\phi\Big{(}\frac{\log(2/\delta)+2d}{N}\Big{)}\Big{)}\leq\delta,$
(38)

where $\phi(x)=\max(x,\sqrt{x}).$

###### Proof.

We first prove (37). Denote $S^{d-1}=\\{u\in\mathbb{R}^{d}\>:\>\|u\|=1\\}$ and
let $u\in S^{d-1}$ and $|\lambda|\leq N/K_{\Xi},$ we have

$\displaystyle\mathbb{E}\exp(\lambda\langle
u,\overline{\Xi}u\rangle)=\prod_{i=1}^{N}\mathbb{E}\exp(\lambda\langle
u,\Xi_{i}u\rangle/N)\leq\prod_{i=1}^{N}\exp(\lambda^{2}K_{\Xi}^{2}/N^{2})=\exp(\lambda^{2}K_{\Xi}^{2}/N),$

so that for all $u\in S^{d-1}$ we have $\langle
u,\overline{\Xi}u\rangle\in\Psi_{1}(K_{\Xi}/\sqrt{N}).$

Let $\Omega_{\epsilon}$ be an $\epsilon$-net of $S^{d-1}.$ By [88, Lemma 5.2],
there exists an $\epsilon$-net such that
$|\Omega_{\epsilon}|\leq(1+2/\epsilon)^{d}$ and for all $u\in S^{d-1}$ there
exists $v\in\Omega_{\epsilon}$ such that $\|u-v\|\leq\epsilon.$ We write

$\langle u,\overline{\Xi}u\rangle=\langle v,\overline{\Xi}v\rangle+2\langle
u-v,\overline{\Xi}v\rangle+\langle u-v,\overline{\Xi}(u-v)\rangle,$

which allows us to deduce that

$\|\overline{\Xi}\|_{2}=\sup_{u\in S^{d-1}}|\langle
u,\overline{\Xi}u\rangle|\leq\sup_{v\in\Omega_{\epsilon}}|\langle
v,\overline{\Xi}v\rangle|+(2\epsilon+\epsilon^{2})\|\overline{\Xi}\|_{2}\implies\|\overline{\Xi}\|_{2}\leq\frac{\sup_{v\in\Omega_{\epsilon}}|\langle
v,\overline{\Xi}v\rangle|}{1-2\epsilon-\epsilon^{2}}.$

Let $v\in\Omega_{\epsilon},$ using Chernoff’s method and the sub-exponential
property of $\overline{\Xi}$ (see also the proof of Corollary 2), we find for
$t>0:$

$\mathbb{P}(|\langle
v,\overline{\Xi}v\rangle|>t)\leq\begin{cases}2\exp(-Nt^{2}/(4K_{\Xi}^{2}))&\text{
if }\quad t\leq 2K_{\Xi}\\\ 2\exp(-Nt/(2K_{\Xi}))&\text{
otherwise.}\end{cases}$

Reformulating in terms of a failure probability $\delta,$ we find that

$\mathbb{P}\Big{(}|\langle
v,\overline{\Xi}v\rangle|>2K_{\Xi}\phi\Big{(}\frac{\log(2/\delta)}{N}\Big{)}\Big{)}\leq\delta.$

Replacing $\delta$ with $\delta/(1+2/\epsilon)^{d}$ and using a union bound
argument over $\Omega_{\epsilon}$ we find

$\mathbb{P}\Big{(}\sup_{v\in\Omega_{\epsilon}}|\langle
v,\overline{\Xi}v\rangle|>2K_{\Xi}\phi\Big{(}\frac{\log(2/\delta)+d\log(1+2/\epsilon)}{N}\Big{)}\Big{)}\leq\delta.$

It only remains to set $\epsilon=1/8$ and plug back into the inequality
$\|\overline{\Xi}\|_{2}\leq\frac{\sup_{v\in\Omega_{\epsilon}}|\langle
v,\overline{\Xi}v\rangle|}{1-2\epsilon-\epsilon^{2}}$ in order to obtain (37).

To prove (38), we proceed similarly and first obtain for all $u\in S^{d-1}$
and $|\lambda|\leq N/K_{\xi}:$

$\mathbb{E}\exp(\langle
u,\overline{\xi}\rangle)\leq\exp(\lambda^{2}K_{\xi}^{2}/N).$

For $u\in S^{d-1}$ and $v\in\Omega_{\epsilon}$ such that
$\|u-v\|\leq\epsilon,$ we write $\langle u,\overline{\xi}\rangle=\langle
v,\overline{\xi}\rangle+\langle u-v,\overline{\xi}\rangle$ which yields the
inequality

$\|\overline{\xi}\|\leq\frac{\sup_{v\in\Omega_{\epsilon}}|\langle
v,\overline{\xi}\rangle|}{1-\epsilon}.$

As before, the sub-exponential property of $\overline{\xi}$ yields

$\mathbb{P}\Big{(}|\langle
v,\overline{\xi}\rangle|>2K_{\xi}\phi\Big{(}\frac{\log(2/\delta)}{N}\Big{)}\Big{)}\leq\delta,$

and using another union bound argument over $\Omega_{\epsilon}$ we find

$\mathbb{P}\Big{(}\sup_{v\in\Omega_{\epsilon}}|\langle
v,\overline{\xi}\rangle|>2K_{\xi}\phi\Big{(}\frac{\log(2/\delta)+d\log(1+2/\epsilon)}{N}\Big{)}\Big{)}\leq\delta.$

It only remains to set $\epsilon=1/2$ to finish the proof of (38). ∎