Nonasymptotic Analysis of Stochastic Gradient Descent with the Richardson–Romberg Extrapolation

## Abstract

Abstract We address the problem of solving strongly convex and smooth minimization problems using stochastic gradient descent (SGD) algorithm with a constant step size. Previous works suggested to combine the Polyak-Ruppert averaging procedure with the Richardson-Romberg extrapolation technique to reduce the asymptotic bias of SGD at the expense of a mild increase of the variance. We significantly extend previous results by providing an expansion of the mean-squared error of the resulting estimator with respect to the number of iterations n 𝑛 n italic_n . More precisely, we show that the mean-squared error can be decomposed into the sum of two terms: a leading one of order 𝒪 ⁢ ( n − 1 / 2 ) 𝒪 superscript 𝑛 1 2 \mathcal{O}(n^{-1/2}) caligraphic_O ( italic_n start_POSTSUPERSCRIPT - 1 / 2 end_POSTSUPERSCRIPT ) with explicit dependence on a minimax-optimal asymptotic covariance matrix, and a second-order term of order 𝒪 ⁢ ( n − 3 / 4 ) 𝒪 superscript 𝑛 3 4 \mathcal{O}(n^{-3/4}) caligraphic_O ( italic_n start_POSTSUPERSCRIPT - 3 / 4 end_POSTSUPERSCRIPT ) where the power 3 / 4 3 4 3/4 3 / 4 can not be improved in general. We also extend this result to the p 𝑝 p italic_p -th moment bound keeping optimal scaling of the remainders with respect to n 𝑛 n italic_n . Our analysis relies on the properties of the SGD iterates viewed as a time-homogeneous Markov chain. In particular, we establish that this chain is geometrically ergodic with respect to a suitably defined weighted Wasserstein semimetric.

## 1 Introduction

Stochastic gradient methods are a fundamental approach for solving a wide range of optimization problems, with a broad range of applications including generative modeling (Goodfellow et al., ; ), empirical risk minimization (Van der Vaart, ), and reinforcement learning (Sutton & Barto, ; Schulman et al., ; Mnih et al., ). These methods are devoted to solving the stochastic minimization problem

$$ $\textstyle\min_{\theta\in\mathbb{R}^{d}}f(\theta)\,,\qquad\nabla f(\theta)= \mathsf{E}_{\xi\sim\mathbb{P}_{\xi}}[\nabla F(\theta,\xi)]\,,$ (1) $$

where $\xi$ is a random variable on $(\mathsf{Z},\mathcal{Z})$ and we can access the gradient $\nabla f$ of the function $f$ only through (unbiased) noisy estimates $\nabla F$. Throughout this paper, we consider strongly convex minimization problems admitting a unique solution $\theta^{\star}$. Arguably the simplest and one of the most widely used approaches to solve (1) is the stochastic gradient descent (SGD), defining the sequence of updates

$$ $\textstyle\theta_{k+1}=\theta_{k}-\gamma_{k+1}\nabla F(\theta_{k},\xi_{k+1})\, ,\quad\theta_{0}\in\mathbb{R}^{d}\,,$ (2) $$

where $\{\gamma_{k}\}_{k\in\mathbb{N}}$ are step sizes, either diminishing or constant, and $\{\xi_{k}\}_{k\in\mathbb{N}}$ is an i.i.d.  sequence with distribution $\mathbb{P}_{\xi}$. The algorithm (2) can be viewed as a special instance of the Robbins-Monro procedure (Robbins & Monro, ). While the SGD algorithm remains one of the core algorithms in statistical inference, its performance can be enhanced by means of additional techniques that use e.g., momentum (Qian, ), averaging (Polyak & Juditsky, ), or variance reduction (Defazio et al., ; Nguyen et al., ). In particular, the celebrated Polyak-Ruppert algorithm proceeds with a trajectory-wise averaging of the estimates

$$ $\textstyle\bar{\theta}_{n_{0},n}=\frac{1}{n}\sum_{k=n_{0}+1}^{n+n_{0}}\theta_{ k}\,$ (3) $$

for some $n_{0}>0.$ It is known (Polyak & Juditsky, ; Fort, ), that under appropriate assumptions on $f$ and $\gamma_{k}$, the sequence of estimates $\{\bar{\theta}_{n_{0},n}\}_{n\in\mathbb{N}}$ is asymptotically normal, that is,

$$ $\sqrt{n}(\bar{\theta}_{n_{0},n}-\theta^{\star})\overset{d}{\to}\mathrm{N}(0, \Sigma_{\infty})\,,\quad n\to\infty$ (4) $$

where $\overset{d}{\to}$ denotes the convergence in distribution and $\mathrm{N}(0,\Sigma_{\infty})$ denotes the zero-mean Gaussian distribution with covariance matrix $\Sigma_{\infty}$, which is asymptotically optimal from the Rao-Cramer lower bound, see Fort () for a discussion. On the other hand, quantitative counterparts of (4) rely on the mean-squared error bounds of the form

$$ $\mathsf{E}^{1/2}[\|\bar{\theta}_{n_{0},n}-\theta^{\star}\|^{2}]\leq\frac{\sqrt {\operatorname{Tr}{\Sigma_{\infty}}}}{n^{1/2}}+\frac{C(f,d)}{n^{1/2+\delta}}+ \mathcal{R}(\|\theta_{0}-\theta^{\star}\|,n)\,.$ (5) $$

Here $\mathcal{R}(\|\theta_{0}-\theta^{\star}\|,n)$ is a remainder term which reflects the dependence upon initial condition, $C(f,d)$ is some instance-dependent constant and $\delta>0$. There are many studies establishing (5) for Polyak-Ruppert averaged SGD under various model assumptions, including Bach & Moulines (), Gadat & Panloup (). In particular, Li et al. () derived the bound (5) with the optimal rate $\delta=1/4$ and proved the corresponding lower bound. However, their results apply a modified algorithm with control variates and multiple restarts. In our work, we show that the same non-asymptotic upper bound is achieved by a simple modification of the estimate $\bar{\theta}_{n_{0},n}$ based on Richardson-Romberg extrapolation. The main contributions of the current paper are as follows:

- •
We show that a version of SGD algorithm with constant step size, Polyak-Ruppert averaging, and Richardson-Romberg extrapolation lead to the optimal root-MSE bound (5) with $\delta=1/4$ when applied to strongly convex minimization problems. We obtain this result by leveraging the analysis of iterates generated by the constant step-size SGD as a Markov chain. This process turns out to be geometrically ergodic with respect to a carefully designed Wasserstein semi-metric (see detailed discussion in Section 3). It is important to note that this result is obtained for a fixed step size $\gamma$ of order $1/\sqrt{n}$ with $n$ being a total number of iterations.
- •
We generalize the above result to the $p$-th moment error bounds. Under a similar step size $\gamma\asymp 1/\sqrt{n},$ we obtain the error bound of the order
$\mathsf{E}^{1/p}[\|\bar{\theta}_{n}^{(RR)}-\theta^{\star}\|^{p}]\leq\frac{Cp^{
1/2}\sqrt{\operatorname{Tr}{\Sigma_{\infty}}}}{n^{1/2}}+\frac{C(f,d,p)}{n^{3/4
}}+\mathcal{R}(\|\theta_{0}-\theta^{\star}\|,n,p)\,,$
(6)
where $C$ is a universal constant, and $\bar{\theta}_{n}^{(RR)}$ is a counterpart of the quantity $\bar{\theta}_{n_{0},n}$ when using Richardson-Romberg extrapolation, see related definitions at Section 4. Our proof is based on a novel version of the Rosenthal inequality, which might be of independent interest.

The rest of the paper is organized as follows. First, we provide a literature review on the non-asymptotic analysis of the first order optimization methods, with an emphasis on the constant step-size algorithms and Richardson-Romberg procedure in Section 2. Then, we provide analysis of the constant step size SGD viewed as a Markov chain together with the properties of the Polyak-Ruppert averaged estimator (3) in Section 3. In Section 4, we discuss the Richardson-Romberg extrapolation applied to the Polyak-Ruppert averaged SGD and derive the respective $2$-nd and $p$-th moment error bounds.

#### Notations and definitions.

For $\theta_{1},\ldots,\theta_{k}$ being the iterates of stochastic first-order method, we denote ${\mathcal{F}}_{k}=\sigma(\theta_{0},\theta_{1},\ldots,\theta_{k})$ and $\mathsf{E}_{k}$ be an alias for $\mathsf{E}\bigl{[}\cdot|{\mathcal{F}}_{k}\bigr{]}$. We call a function $c:\mathsf{Z}\times\mathsf{Z}\rightarrow\mathbb{R}_{+}$ a *distance-like* function, if it is symmetric, lower semi-continuous and
$c(x,y)=0$ if and only if $x=y$, and there exists $q\in\mathbb{N}$ such that
$(d(x,y)\wedge 1)^{q}\leq c(x,y)$. For two probability measures $\xi$ and $\xi^{\prime}$ we denote by $\mathscr{C}(\xi,\xi^{\prime})$ the set of couplings of two probability measures, that is, for any $\mathcal{C}\in\mathscr{C}(\xi,\xi^{\prime})$ and any $A\in\mathcal{Z}$ it holds $\mathcal{C}(\mathsf{Z}\times A)=\xi^{\prime}(A)$ and $\mathcal{C}(A\times\mathsf{Z})=\xi(A)$. We define the Wasserstein semi-metrics associated to the distance-like function $c(\cdot,\cdot)$, as

$$ $\textstyle{\mathbf{W}_{{c}}(\xi,\xi^{\prime})=\underset{\mathcal{C}\in\mathscr {C}(\xi,\xi^{\prime})}{\inf}\int_{\mathsf{Z}\times\mathsf{Z}}c(z,z^{\prime}) \mathcal{C}(\mathrm{d}z,\mathrm{d}z^{\prime})}\,.$ (7) $$

Note that $\mathbf{W}_{{c}}(\xi,\xi^{\prime})$ is not necessarily a distance, as it may fail to satisfy the triangle inequality. In the particular case of $\mathsf{Z}=\mathbb{R}^{d}$, and $c_{p}(x,y)=\|x-y\|^{p}$, $x,y\in\mathbb{R}^{d}$, $p\geq 1$, we denote the corresponding Wasserstein metrics by $\mathbf{W}_{{p}}(\xi,\xi^{\prime})$. Let ${\rm Q}(z,A)$ be a Markov kernel on $(\mathsf{Z},\mathcal{Z})$. We say that ${\rm K}$ is a Markov coupling of ${\rm Q}$ if for all $(z,z^{\prime})\in\mathsf{Z}^{2}$ and $A\in\mathcal{Z}$, ${\rm K}((z,z^{\prime}),A\times\mathsf{Z})={\rm Q}(z,A)$ and ${\rm K}((z,z^{\prime}),\mathsf{Z}\times A)={\rm Q}(z^{\prime},A)$.
If ${\rm K}$ is a kernel coupling of ${\rm Q}$, then for all $n\in\mathbb{N}$, ${\rm K}^{n}$ is a kernel coupling of ${\rm Q}^{n}$ and for any $\mathcal{C}\in\mathscr{C}(\xi,\xi^{\prime})$, $\mathcal{C}{\rm K}^{n}$ is a coupling of $(\xi{\rm Q}^{n},\xi^{\prime}{\rm Q}^{n})$ and it holds

$$ $\textstyle\mathbf{W}_{{c}}(\xi{\rm Q}^{n},\xi^{\prime}{\rm Q}^{n})\leq\int_{ \mathsf{Z}\times\mathsf{Z}}{\rm K}^{n}c(z,z^{\prime})\mathcal{C}(\mathrm{d}z \mathrm{d}z^{\prime})\,,$ $$

see (Douc et al., , Corollary 20.1.4). For any probability measure $\mathcal{C}$ on $(\mathsf{Z}^{2},\mathcal{Z}^{\otimes 2})$, we denote by $\mathbb{P}_{\zeta}^{\rm K}$ and $\mathsf{E}_{\zeta}^{\rm K}$ the probability and the
expectation on the canonical space $((\mathsf{Z}^{2})^{\mathbb{N}},(\mathcal{Z}^{\otimes 2})^{\otimes\mathbb{N}})$ such that the canonical process
$\{(Z_{n},Z_{n}^{\prime}),n\in\mathbb{N}\}$ is a Markov chain with initial
probability $\mathcal{C}$ and Markov kernel ${\rm K}$. We write
$\mathsf{E}_{z,z^{\prime}}^{{\rm K}}$ instead $\mathsf{E}_{\delta_{z,z^{\prime}}}^{{\rm K}}$. For all $x,y\in\mathbb{R}^{d}$ denote by $x\otimes y$ the tensor product of $x$ and $y$ and by $x^{\otimes k}$ the $k$-th tensor power of $x$. In addition, for a function $f:\mathbb{R}^{d}\rightarrow\mathbb{R}$ we denote by $\nabla^{k}f(\theta)$ the $k$-th differential of $f$, that is $\nabla^{k}f(\theta)_{i_{1},\cdots,i_{k}}=\frac{\partial^{k}f}{\prod_{j=1}^{k}
\partial x_{i_{j}}}$. For any tensor $M\in(\mathbb{R}^{d})^{\otimes(k-1)}$, we define $\nabla^{k}f(\theta)M\in\mathbb{R}^{d}$ by the relation $(\nabla^{k}f(\theta)M)_{l}=\sum_{i_{1},\ldots,i_{k-1}}M_{i_{1},\cdots,i_{k-1}}
\nabla^{k}f(\theta)_{i_{1},\cdots,i_{k-1},l}$, where $l\in\{1,\ldots,d\}$.

For two sequences $\{a_{n}\}_{n\in\mathbb{N}}$ and $\{b_{n}\}_{n\in\mathbb{N}}$ we write $a_{n}\lesssim b_{n}$, if there is an absolute constant $c_{0}$, such that $a_{n}\leq c_{0}b_{n}$.

## 2 Literature review

Gadat & Panloup () derive (5) with $\delta=1/8$ for a certain class of functions $f$, including strongly convex functions, improving Moulines & Bach () which obtain this result for $\delta=1/12$. Li et al. () suggested the Root-SGD algorithm combining the ideas of the Polyak-Ruppert averaged SGD with control variates and established (5) with $\delta=1/4$. The recent series of papers (Huo et al., ; Zhang & Xie, ; Zhang et al., ) investigate stochastic approximation algorithms with both i.i.d. and Markovian data and constant step sizes. The authors consider both linear SA problems and $Q$-learning, quantify bias, and propose precise characterization of the bias together with a Richardson-Romberg extrapolation procedure. However, these results only consider $2$-nd moment of the error and provide MSE bounds of order $\mathcal{O}(1/n)+\mathcal{O}(\gamma)$ with no explicit expression for the leading term.

#### Richardson-Romberg extrapolation.

Richardson-Romberg extrapolation is a technique used to improve the accuracy of numerical approximations (Hildebrand, ), such as those from numerical differentiation or integration. It involves using approximations with different step sizes and then extrapolating to reduce the error, typically by removing the leading term in the error expansion. The one-step Richardson-Romberg was introduced to reduce the discretization error induced by an Euler scheme to simulate stochastic differential equation in Talay & Tubaro (), and later generalized for non-smooth functions in Bally & Talay (). This technique was extended using multistep discretizations in Pagès (). Finally, Richardson-Romberg extrapolation have been applied to Stochastic Gradient Descent (SGD) methods in Durmus et al. (), Merad & Gaïffas () and Huo et al. (), to improve convergence and reduce error in optimization problems, particularly when dealing with noisy or high-variance gradient estimates.

## 3 Finite-time analysis of the SGD dynamics for strongly convex minimization problems

### 3.1 Geometric ergodicity of SGD iterates

We consider the following assumption on the function $f$ in the minimization problem (1).

###### A 1 .

The function $f$ is $\mu$-strongly convex on $\mathbb{R}^{d}$, that is, it is continuously differentiable and there exists a constant $\mu>0$, such that for any $\theta,\theta^{\prime}\in\mathbb{R}^{d}$, it holds that

$$ $\frac{\mu}{2}\|\theta-\theta^{\prime}\|^{2}\leq f(\theta)-f(\theta^{\prime})- \langle\nabla f(\theta^{\prime}),\theta-\theta^{\prime}\rangle\,.$ (8) $$

###### A 2 .

The function $f$ is $4$ times continuously differentiable and $\operatorname{L}_{2}$-smooth on $\mathbb{R}^{d}$, i.e., it is continuously differentiable and there is a constant $\operatorname{L}_{2}\geq 0$, such that for any $\theta,\theta^{\prime}\in\mathbb{R}^{d}$,

$$ $\|\nabla f(\theta)-\nabla f(\theta^{\prime})\|\leq\operatorname{L}_{2}\|\theta -\theta^{\prime}\|\,.$ (9) $$

Moreover, $f$ has uniformly bounded $3$-rd and $4$-th derivatives, there exist $\operatorname{L}_{3},\operatorname{L}_{4}\geq 0$ such that

$$ $\|\nabla^{i}f(\theta)\|\leq\operatorname{L}_{i}\,\text{ for }i\in\{3,4\}\,.$ (10) $$

We aim to solve the problem (1) using SGD with a constant step size, starting from initial distribution $\nu$. That is, for $k\geq 0$ and a step size $\gamma\geq 0,$ we consider the following recurrent scheme

$$ $\theta^{(\gamma)}_{k+1}=\theta^{(\gamma)}_{k}-\gamma\nabla F(\theta^{(\gamma)} _{k},\xi_{k+1})\,,\quad\theta^{(\gamma)}_{0}=\theta_{0}\sim\nu\,,$ (11) $$

where $\{\xi_{k}\}_{k\in\mathbb{N}}$ is a sequence satisfying the following condition.

###### A 3 (p) .

$\{\xi_{k}\}_{k\in\mathbb{N}}$ is a sequence of independent and identically distributed (i.i.d.) random variables with distribution $\mathbb{P}_{\xi}$, such that $\xi_{i}$ and $\theta_{0}$ are independent and for any $\theta\in\mathbb{R}^{d}$ it holds that

$$ $\mathsf{E}_{\xi\sim\mathbb{P}_{\xi}}[\nabla F(\theta,\xi)]=\nabla f(\theta)\,.$ $$

Moreover, there exists $\tau_{p}$, such that $\mathsf{E}^{1/p}[\|\nabla F(\theta^{\star},\xi)\|^{p}]\leq\tau_{p}$, and for any $q=2,\ldots,p$ it holds with some $\operatorname{L}_{1}>0$ that for any $\theta_{1},\theta_{2}\in\mathbb{R}^{d}$,

$$ $\operatorname{L}_{1}^{q-1}\|\theta_{1}-\theta_{2}\|^{q-2}\langle\nabla f( \theta_{1})-\nabla f(\theta_{2}),\theta_{1}-\theta_{2}\rangle\geq\mathsf{E}_{ \xi\sim\mathbb{P}_{\xi}}[\|\nabla F(\theta_{1},\xi)-\nabla F(\theta_{2},\xi)\| ^{q}]\,.$ (12) $$

Assumption A 3($p$) generalizes the well-known $\operatorname{L}_{1}$-co-coercivity assumption, see Dieuleveut et al. (). A sufficient condition which allows for A 3($p$) is to assume that $F(\theta,\xi)$ is $\mathbb{P}_{\xi}$-a.s. convex with respect to $\theta\in\mathbb{R}^{d}$. For ease of notation, we set

$$ $\textstyle\operatorname{L}=\max(\operatorname{L}_{1},\operatorname{L}_{2}, \operatorname{L}_{3},\operatorname{L}_{4})\,,$ (13) $$

and trace only $\operatorname{L}$ in our subsequent bounds. In this paper we focus on the convergence to $\theta^{\star}$ of the Polyak-Ruppert averaging estimator defined for any $n\geq 0$,

$$ $\textstyle\bar{\theta}_{n}^{(\gamma)}=\frac{1}{n}\sum_{k=n+1}^{2n}\theta^{( \gamma)}_{k}\,.$ (14) $$

Many previous studies instead consider $\bar{\vartheta}_{n}^{(\gamma)}=\frac{1}{n-n_{0}}\sum_{k=n_{0}+1}^{n}\theta_{k}
^{(\gamma)}$ rather than $\bar{\theta}_{n}^{(\gamma)}$, where $n\geq n_{0}+1$ and $n_{0}$ denotes a burn-in period. However, when the sample size $n$ is sufficiently large, the choice of the optimal burn-in size $n_{0}$ affects the leading terms in the MSE bound of $\bar{\theta}_{n}^{(\gamma)}-\theta^{\star}$ only by a constant factor. Therefore, we focus on (14), or equivalently, use $2n$ observations and set $n_{0}=n$.

###### A 1 .

###### A 2 .

###### A 3 (p) .

#### Properties of { θ k ( γ ) } k ∈ ℕ subscript superscript subscript 𝜃 𝑘 𝛾 𝑘 ℕ \{\theta_{k}^{(\gamma)}\}_{k\in\mathbb{N}} { italic_θ start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT start_POSTSUPERSCRIPT ( italic_γ ) end_POSTSUPERSCRIPT } start_POSTSUBSCRIPT italic_k ∈ blackboard_N end_POSTSUBSCRIPT viewed as a Markov chain.

Under assumptions A 1, A 2 and A 3($2$), the sequence $\{\theta_{k}^{(\gamma)}\}_{k\in\mathbb{N}}$ defined by the relation (11) is a time-homogeneous Markov chain with the Markov kernel

$$ ${\rm Q}_{\gamma}(\theta,\mathsf{A})=\int_{\mathbb{R}^{d}}\mathbbm{1}_{\mathsf{ A}}(\theta-\gamma\nabla F(\theta,z))\mathsf{P}_{\xi}(\mathrm{d}z)\,,\quad \theta\in\mathbb{R}^{d}\,,\,\mathsf{A}\in\mathsf{B}(\mathbb{R}^{d})\,,$ (15) $$

where $\mathsf{B}(\mathbb{R}^{d})$ denoted the Borel $\sigma$-field of $\mathbb{R}^{d}$. In Dieuleveut et al. () it has been established that, under the stated assumptions, ${\rm Q}_{\gamma}$ admits a unique invariant distribution $\pi_{\gamma}$, if $\gamma$ is small enough. Previous studies, such as Dieuleveut et al. () or Merad & Gaïffas (), studied the convergence of the distributions of $\{\theta_{k}^{(\gamma)}\}_{k\in\mathbb{N}}$ to $\pi_{\gamma}$ in the $2$-Wasserstein distance $\mathbf{W}_{{2}}$, associated with the Euclidean distance in $\mathbb{R}^{d}$. However, our main results would require to switch to the non-standard distance-like function, which is defined under A 1 and A 3($2$) as follows:

$$ $\textstyle c(\theta,\theta^{\prime})=\|\theta-\theta^{\prime}\|\left(\|\theta- \theta^{*}\|+\|\theta^{\prime}-\theta^{*}\|+\frac{2\sqrt{2}\tau_{2}\sqrt{ \gamma}}{\sqrt{\mu}}\right)\,,\qquad\theta,\theta^{\star}\in\mathbb{R}^{d}\,.$ (16) $$

Here the constants $\tau_{2}$ and $\mu$ are given in A 3($2$) and A 1, respectively. Note that this distance-like function is specifically designed to analyze $\{\theta_{k}^{(\gamma)}\}_{k\in\mathbb{N}}$ under A 1 and A 3($2$). In particular, it depends on the step size $\gamma$ and $\theta^{\star}$. Our first main result establishes geometric ergodicity of the sequence $\{\theta_{k}^{(\gamma)}\}_{k\in\mathbb{N}}$ with respect to the distance-like function $c$ from (16).

###### Proposition \theproposition .

Assume A 1, A 2, and A 3($2$). Then for any $\gamma\in(0;1/(2\operatorname{L})],$ the Markov kernel ${\rm Q}_{\gamma}$ defined in (15) admits a unique invariant distribution $\pi_{\gamma}$. Moreover, for any initial distribution $\nu$ on $\mathbb{R}^{d}$ and $k\in\mathbb{N}$,

$$ $\mathbf{W}_{{c}}(\nu{\rm Q}_{\gamma}^{k},\pi_{\gamma})\leq 4(1/2)^{k/m(\gamma) }\mathbf{W}_{{c}}(\nu,\pi_{\gamma})\,,$ (17) $$

where $m(\gamma)=\lceil 2\log 4/(\gamma\mu)\rceil$.

###### Proposition \theproposition .

#### Discussion.

The proof of Section 3.1 is provided in Section A.1. Properties of the invariant distribution $\pi_{\gamma}$ were previously studied in literature, see e.g. Dieuleveut et al. (). It particular, it is known (Dieuleveut et al., , Lemma 13), that the $2$-nd moment of $\theta_{\infty}^{(\gamma)}$, where $\theta_{\infty}^{(\gamma)}$ is distributed according to the stationary distribution $\pi_{\gamma}$, scales linearly with $\gamma$:

$$ $\textstyle\int_{\mathbb{R}^{d}}\|\theta-\theta^{\star}\|^{2}\pi_{\gamma}( \mathrm{d}\theta)\leq\frac{\mathsf{D}_{\operatorname{last},2}\gamma\tau_{2}}{ \mu}\,.$ (18) $$

This property yields, using Lyapunov’s inequality, that

$$ $\textstyle\int_{\mathbb{R}^{d}\times\mathbb{R}^{d}}\|\theta-\theta^{\prime}\| \pi_{\gamma}(\mathrm{d}\theta)\pi_{\gamma}(\mathrm{d}\theta^{\prime})\lesssim \sqrt{\frac{\mathsf{D}_{\operatorname{last},2}\gamma\tau_{2}}{\mu}}\,.$ $$

At the same time, expectation of the cost function $c(\theta,\theta^{\prime})$ scale linearly with the step size $\gamma$:

$$ $\textstyle\int_{\mathbb{R}^{d}}c(\theta,\theta^{\prime})\pi(\mathrm{d}\theta) \pi(\mathrm{d}\theta^{\prime})\lesssim\frac{\mathsf{D}_{\operatorname{last},2} \gamma\tau_{2}}{\mu}\,.$ (19) $$

The property (19) is crucial to obtain tighter (with respect to the step size $\gamma$) error bounds for the Richardson-Romberg estimator, as well as in the Rosenthal inequality for additive functional of $\{\theta_{k}^{(\gamma)}\}_{k\in\mathbb{N}}$ derived in Section 4.

Next, we analyze the error $\theta_{\infty}^{(\gamma)}-\theta^{\star}$ where $\theta_{\infty}^{(\gamma)}$ is distributed according to the stationary distribution $\pi_{\gamma}$. To this end, we consider the following condition.

###### C 1 (p) .

There exist constants $\mathsf{D}_{\operatorname{last},p},\mathsf{C}_{\operatorname{step},p}\geq 2$ depending only on $p$, such that for any step size $\gamma\in(0,1/(\operatorname{L}\mathsf{C}_{\operatorname{step},p})]$, and any initial distribution $\nu$ it holds that

$$ $\textstyle\mathsf{E}_{\nu}^{2/p}\bigl{[}\|\theta_{k}^{(\gamma)}-\theta^{\star} \|^{p}\bigr{]}\leq(1-\gamma\mu)^{k}\mathsf{E}_{\nu}^{2/p}\bigl{[}\|\theta_{0}- \theta^{\star}\|^{p}\bigr{]}+\mathsf{D}_{\operatorname{last},p}\gamma\tau_{p}^ {2}/\mu\,.$ (20) $$

Moreover, for the stationary distribution $\pi_{\gamma}$ it holds that

$$ $\mathsf{E}_{\pi_{\gamma}}^{2/p}\bigl{[}\|\theta_{\infty}^{(\gamma)}-\theta^{ \star}\|^{p}\bigr{]}\leq\mathsf{D}_{\operatorname{last},p}\gamma\tau_{p}^{2}/ \mu\,.$ (21) $$

Note that, since $\operatorname{L}\geq\mu$, it holds that $\gamma\mu\leq 1/2$. It is important to recognize that C 1 is not independent from the previous assumptions A 1 - A 3($p$). In particular, Dieuleveut et al. (, Lemma $13$) implies that, under A 1, A 3($p$) with $p\geq 2$, and A 2, the bound (21) holds for $\gamma\in(0,1/(\operatorname{L}\mathsf{C}_{\operatorname{step},p})]$ with some constants $\mathsf{D}_{\operatorname{last},p}$ and $\mathsf{C}_{\operatorname{step},p}$, which depends only upon $p$. Unfortunately, it is complicated to obtain precise dependence of $\mathsf{C}_{\operatorname{step},p}$ and $\mathsf{D}_{\operatorname{last},p}$ upon $p$, as well as to obtain the bound (21) with tight numerical constants. The results available in the literature (Gadat & Panloup, ; Li et al., ; Merad & Gaïffas, ) either are obtained under alternative set of assumptions, or are not explicit with respect to their dependence upon $p$. That is why we prefer to state C 1($p$) as a separate assumption. In the subsequent bounds we use C 1($p$) together with A 1, A 3($p$) with $p\geq 2$, and A 2, tracking the dependence of our bounds upon $\mathsf{C}_{\operatorname{step},p}$ and $\mathsf{D}_{\operatorname{last},p}$. We leave the problem of deriving C 1($p$) with sharp constants $\mathsf{D}_{\operatorname{last},p},\mathsf{C}_{\operatorname{step},p}$ as an interesting direction for the future research.

Under the assumption C 1, we control the fluctuations of $\{\theta_{k}^{(\gamma)}\}$ around the solution $\theta^{\star}$ of (1). However, unless the function $f$ is quadratic, it is known that $\int_{\mathbb{R}^{d}}\theta\pi_{\gamma}(\mathrm{d}\theta)\neq\theta^{\star}$.
In the following proposition, we quantify this bias under milder assumptions compared to the ones from Dieuleveut et al. (, Theorem 4). Namely, the following result holds:

###### Proposition \theproposition .

Assume A 1, A 2, A 3($6$), and C 1($6$). Then there exist such $\Delta_{1}\in\mathbb{R}^{d},\Delta_{2}\in\mathbb{R}^{d\times d}$, not depending upon $\gamma$, that for any $\gamma\in(0,1/(\operatorname{L}\mathsf{C}_{\operatorname{step},6})]$, it holds

$$ $\bar{\theta}_{\gamma}:=\int_{\mathbb{R}^{d}}\theta\pi_{\gamma}(\mathrm{d} \theta)=\theta^{\star}+\gamma\Delta_{1}+B_{1}\gamma^{3/2}\,,$ (22) $$

$$ $\bar{\Sigma}_{\gamma}:=\int_{\mathbb{R}^{d}}(\theta-\theta^{\star})^{\otimes 2 }\pi_{\gamma}(\mathrm{d}\theta)=\gamma\Delta_{2}+B_{2}\gamma^{3/2}\,.$ (23) $$

Here $B_{1}\in\mathbb{R}^{d}$ and $B_{2}\in\mathbb{R}^{d\times d}$ satisfy $\|B_{i}\|\leq\mathsf{C}_{1}$, $i=1,2$, where $\mathsf{C}_{1}$ defined in (50) is a constant independent of $\gamma$. Moreover, for any initial distribution $\nu$ on $\mathbb{R}^{d}$, it holds that

$$ $\mathsf{E}_{\nu}[\bar{\theta}_{n}^{(\gamma)}]=\theta^{\star}+\gamma\Delta_{1}+ B_{1}\gamma^{3/2}+\mathcal{R}_{1}(\theta_{0}-\theta^{\star},\gamma,n)\,,$ (24) $$

where

$$ $\textstyle\|\mathcal{R}_{1}(\theta_{0}-\theta^{\star},\gamma,n)\|\lesssim\frac {\mathrm{e}^{-\gamma\mu(n+1)/2}}{n\gamma\mu}\left(\mathsf{E}_{\nu}^{1/2}\bigl{ [}\|\theta_{0}-\theta^{\star}\|^{2}\bigr{]}+\frac{\sqrt{\gamma}\tau_{2}}{\sqrt {\mu}}\right)\,.$ (25) $$

The proof is postponed to Appendix A. Results of this type were already obtained in the literature for stochastic approximation algorithms, see e.g. Huo et al. () and Allmeier & Gast (). As already highlighted, the additive term $\Delta_{1}$ vanishes in the case of minimizing the quadratic function $f$, see Bach & Moulines ().

###### C 1 (p) .

###### Proposition \theproposition .

### 3.2 Analysis of the Polyak-Ruppert averaged estimator θ ¯ n ( γ ) superscript subscript ¯ 𝜃 𝑛 𝛾 \bar{\theta}_{n}^{(\gamma)} over¯ start_ARG italic_θ end_ARG start_POSTSUBSCRIPT italic_n end_POSTSUBSCRIPT start_POSTSUPERSCRIPT ( italic_γ ) end_POSTSUPERSCRIPT .

In this section, we analyze the finite-sample properties of the estimator $\bar{\theta}_{n}^{(\gamma)}$ from (14). The analysis of this section is based on techniques previously used in Moulines & Bach (), as well as in the analysis of the Polyak-Ruppert averaged LSA (Linear Stochastic Approximation) algorithms, see Mou et al. (); Durmus et al. (). Next, we define the $k$-th step noise level at the point $\theta\in\mathbb{R}^{d}$ by:

$$ ${\varepsilon_{k}}(\theta)=\nabla F(\theta,\xi_{k})-\nabla f(\theta)\,,$ (26) $$

Note that ${\varepsilon_{k+1}}(\theta_{k}^{(\gamma)})$ is a martingale-difference sequence w.r.t. the filtration $({\mathcal{F}}_{k})_{k\in\mathbb{N}}$. With these notations, (11) writes

$$ $\theta_{k+1}^{(\gamma)}-\theta^{\star}=\theta_{k}^{(\gamma)}-\theta^{\star}- \gamma\bigl{(}\nabla f(\theta_{k}^{(\gamma)})+{\varepsilon_{k+1}}(\theta_{k}^{ (\gamma)})\bigr{)}\,.$ (27) $$

Setting

$$ $\eta(\theta)=\nabla f(\theta)-\operatorname{H}^{\star}(\theta-\theta^{\star})\,,$ (28) $$

with

$$ $\operatorname{H}^{\star}=\nabla^{2}f(\theta^{\star})\in\mathbb{R}^{d\times d}\,.$ (29) $$

we get from (27) that

$$ $\theta_{k+1}^{(\gamma)}-\theta^{\star}=(\mathrm{I}-\gamma\operatorname{H}^{ \star})(\theta_{k}^{(\gamma)}-\theta^{\star})-\gamma{\varepsilon_{k+1}}(\theta _{k}^{(\gamma)})-\gamma\eta(\theta_{k}^{(\gamma)})\,,$ $$

or, rearranging the terms,

$$ $\operatorname{H}^{\star}(\theta_{k}^{(\gamma)}-\theta^{\star})=\frac{\theta_{k }^{(\gamma)}-\theta_{k+1}^{(\gamma)}}{\gamma}-{\varepsilon_{k+1}}(\theta_{k}^{ (\gamma)})-\eta(\theta_{k}^{(\gamma)})\,.$ (30) $$

Taking average of (30) with $k=n+1$ to $2n$, we arrive at the final representation:

$$ $\operatorname{H}^{\star}(\bar{\theta}_{n}^{(\gamma)}-\theta^{\star})=\frac{ \theta_{n+1}^{(\gamma)}-\theta^{\star}}{\gamma n}-\frac{\theta_{2n}^{(\gamma)} -\theta^{\star}}{\gamma n}-\frac{1}{n}\sum_{k=n+1}^{2n}{\varepsilon_{k+1}}( \theta_{k}^{\gamma})-\frac{1}{n}\sum_{k=n+1}^{2n}\eta(\theta_{k}^{\gamma})\,.$ (31) $$

Representation (31) is the key ingredient in the proof of the next result where the variance of noise ${\varepsilon_{k}}(\theta^{\star})$ measured at the optimal point $\theta^{\star}$ naturally appears, that is,

$$ $\Sigma_{\varepsilon}^{\star}=\mathsf{E}_{\xi\sim\mathbb{P}_{\xi}}[\nabla F( \theta^{\star},\xi)^{\otimes 2}]\,.$ (32) $$

Note that $\Sigma_{\varepsilon}^{\star}$ does not depend on the step size $\gamma$ of (11) and is related to the ”optimal” covariance matrix of the Polyak-Ruppert averaged iterates $\bar{\theta}_{n}^{(\gamma)}$, see Fort (). In our first main result below, we establish non-asymptotic properties of the averaged Polyak-Ruppert estimator (14).

###### Theorem 1 .

Assume A 1, A 2, A 3($6$), and C 1($6$). Then for any $\gamma\in(0,1/(\operatorname{L}\mathsf{C}_{\operatorname{step},6})]$, $n\in\mathbb{N}$, and initial distribution $\nu$ on $\mathbb{R}^{d}$, the sequence of Polyak-Ruppert estimates (14) satisfies

$$ $\mathsf{E}_{\nu}^{1/2}[\|\operatorname{H}^{\star}(\bar{\theta}_{n}^{(\gamma)}- \theta^{\star})\|^{2}]\lesssim\frac{\sqrt{\operatorname{Tr}{\Sigma_{ \varepsilon}^{\star}}}}{\sqrt{n}}+\frac{\mathsf{C}_{2}}{\gamma^{1/2}n}+\mathsf {C}_{3}\gamma+\frac{\mathsf{C}_{4}\gamma^{1/2}}{n^{1/2}}+\mathcal{R}_{2}(n, \gamma,\|\theta_{0}-\theta^{\star}\|)\,,$ (33) $$

where the constants $\mathsf{C}_{2}$ to $\mathsf{C}_{4}$ are defined in Appendix B (see equation (60)), and

$$ $\mathcal{R}_{2}(n,\gamma,\|\theta_{0}-\theta^{\star}\|)=\frac{(1-\gamma\mu)^{( n+1)/2}}{\gamma n}\left(\frac{\operatorname{L}}{\mu}+1\right)\mathsf{E}_{\nu}^ {1/2}\bigl{[}\|\theta_{0}-\theta^{\star}\|^{2}\bigr{]}\\ +\frac{\operatorname{L}(1-\gamma\mu)^{n+1}}{2n\gamma\mu}\mathsf{E}^{1/2}_{\nu} \bigl{[}\|\theta_{0}-\theta^{\star}\|^{4}\bigr{]}\,.$ $$

The version of Theorem 1 with explicit constants together with the proof is provided in Appendix B, see Theorem 4. Note that the result of Theorem 1 is valid for arbitrary $\gamma\in(0,1/(\operatorname{L}\mathsf{C}_{\operatorname{step},6})]$. At the same time, this bound can be optimized over step size of the form $\gamma=n^{-\beta}$, $\beta\in(0,1)$.

###### Corollary \thecorollary .

Under the assumptions of Theorem 1, by setting $\gamma=n^{-2/3}$ with $n\geq\left(\operatorname{L}\mathsf{C}_{\operatorname{step},6}\right)^{3/2}$ it holds that

$$ $\textstyle\mathsf{E}^{1/2}_{\nu}[\|\operatorname{H}^{\star}(\bar{\theta}_{n}^{ (\gamma)}-\theta^{\star})\|^{2}]\lesssim\frac{\sqrt{\operatorname{Tr}{\Sigma_{ \varepsilon}^{\star}}}}{n^{1/2}}+\mathcal{O}\bigl{(}n^{-2/3}\bigr{)}+\mathcal{ R}_{2}(n,1/n^{2/3},\|\theta_{0}-\theta^{\star}\|)\,.$ (34) $$

Section 3.2 implies that, if $n$ is known in advance and $\gamma\asymp n^{-2/3}$, then $\bar{\theta}_{n}^{(\gamma)}$ satisfies (5) with $\delta=1/6$. A closer inspection of the sum (31) reveals that $\mathsf{E}_{\pi_{\gamma}}[\eta(\theta_{k}^{\gamma})]\asymp\gamma$ meaning that we can not hope to provide a better bound for the term $\frac{1}{n}\sum_{k=n+1}^{2n}\eta(\theta_{k}^{\gamma})$ compared to the one coming from the Minkowski’s inequality. Thus, this is the *bias* of the stationary distribution, which prevents us from gaining the optimal second-order term wrt the sample size $n$ from Section 3.2.

The decay rate of the second-order (w.r.t. $n$) term in (34) is well studied in the literature. In particular, Moulines & Bach () obtained a second order term of order $n^{-5/8}$ for the SGD algorithm with Polyak-Ruppert averaging. A similar rate under more general assumptions was reported in (Gadat & Panloup, , Theorem 2) for $p$-th moment bounds. However, all these results are known to be suboptimal for first-order methods. The recent work by Li et al. () shows that an optimal second-order error term in the bounds (34) is of order $\mathcal{O}\bigl{(}n^{-3/4}\bigr{)}$ and can be achieved by the Root-SGD algorithm. In fact, the Root-SGD algorithm requires the same smoothness assumptions (bounded $3$rd derivative) as Theorem 1. In the next section we can mirror this bound using the constant step-size SGD algorithm combined with the Richardson-Romberg extrapolation technique (Dieuleveut et al., ).

###### Theorem 1 .

###### Corollary \thecorollary .

## 4 Richardson-Romberg extrapolation

Our analysis presented in Theorem 1 was based on the summation by parts formula (31) and Taylor expansion of the gradient $\nabla f(\theta)$ in the vicinity of $\theta^{\star}$, yielding the remainder quantity $\eta(\theta)$. It is important to notice that

$$ $\textstyle\int_{\mathbb{R}^{d}}\eta(\theta)\pi_{\gamma}(\mathrm{d}\theta)\neq 0\,,$ (35) $$

which prevents us from using more aggressive (larger) step sizes $\gamma$ in the optimized bound (34). In this section we show that Richardson-Romberg extrapolation technique is sufficient to significantly reduce the bias associated with $\eta(\theta)$ and improve the second-order term in the MSE bound (34). Instead of considering a single SGD trajectory $\{\theta_{k}^{(\gamma)}\}_{k\in\mathbb{N}}$, and then relying on the tail-averaged estimator $\bar{\theta}_{n}^{(\gamma)}$, we construct two parallel chains based on the same sequence of noise variables $\{\xi_{k}\}_{k\in\mathbb{N}}$:

$$ $\begin{split}\textstyle\theta_{k+1}^{(\gamma)}&=\textstyle\theta_{k}^{(\gamma) }-\gamma\nabla F(\theta_{k}^{(\gamma)},\xi_{k+1})\,,\quad\bar{\theta}_{n}^{( \gamma)}=\frac{1}{n}\sum_{k=n+1}^{2n}\theta_{k}^{(\gamma)}\,,\\ \textstyle\theta_{k+1}^{(2\gamma)}&=\textstyle\theta_{k}^{(2\gamma)}-2\gamma \nabla F(\theta_{k}^{(2\gamma)},\xi_{k+1})\,,\quad\bar{\theta}_{n}^{(2\gamma)} =\frac{1}{n}\sum_{k=n+1}^{2n}\theta_{k}^{(2\gamma)}\,.\end{split}$ (36) $$

Based on $\bar{\theta}_{n}^{(\gamma)}$ and $\bar{\theta}_{n}^{(2\gamma)}$ defined above, we construct a Richardson-Romberg estimator as

$$ $\bar{\theta}_{n}^{(RR)}:=2\bar{\theta}_{n}^{(\gamma)}-\bar{\theta}_{n}^{(2 \gamma)}\,.$ (37) $$

Note that it is possible to use different sources of randomness $\{\xi_{k}\}_{k\in\mathbb{N}}$ and $\{\xi^{\prime}_{k}\}_{k\in\mathbb{N}}$ when constructing the sequences $\{\theta_{k}^{(\gamma)}\}_{k\in\mathbb{N}}$ and $\{\theta_{k}^{(2\gamma)}\}_{k\in\mathbb{N}}$, respectively. At the same time, it is possible to show the benefits of using the same sequence of random variables $\{\xi_{k}\}_{k\in\mathbb{N}}$ in (36). Indeed, consider the decomposition (31) and further expand the term $\eta(\theta)$ defined in (28) as

$$ $\eta(\theta)=\psi(\theta)+G(\theta)\,,$ $$

where we have defined, for $\theta\in\mathbb{R}^{d}$, the following vector-valued functions:

$$ $\textstyle\begin{split}\psi(\theta)&=(1/2)\nabla^{3}f(\theta^{*})(\theta- \theta^{\star})^{\otimes 2}\,,\\ G(\theta)&=\textstyle(1/6)\left(\int_{0}^{1}\nabla^{4}f(t\theta^{\star}+(1-t) \theta)\,\mathrm{d}t\right)(\theta-\theta^{\star})^{\otimes 3}\,.\end{split}$ (38) $$

We further rewrite the decomposition (31) as

$$ $\operatorname{H}^{\star}(\bar{\theta}_{n}^{(\gamma)}-\theta^{\star})=\frac{ \theta_{n+1}^{(\gamma)}-\theta^{\star}}{\gamma n}-\frac{\theta_{2n}^{(\gamma)} -\theta^{\star}}{\gamma n}-\frac{1}{n}\sum_{k=n+1}^{2n}{\varepsilon_{k+1}}( \theta^{\star})\\ -\frac{1}{n}\sum_{k=n+1}^{2n}\{{\varepsilon_{k+1}}(\theta_{k}^{(\gamma)})-{ \varepsilon_{k+1}}(\theta^{\star})\}-\frac{1}{n}\sum_{k=n+1}^{2n}\psi(\theta_{ k}^{(\gamma)})-\frac{1}{n}\sum_{k=n+1}^{2n}G(\theta_{k}^{(\gamma)})\,.$ (39) $$

Note that in the decomposition (39), the linear statistics $W=n^{-1}\sum_{k=n+1}^{2n}{\varepsilon_{k+1}}(\theta^{\star})$ does not depend upon $\gamma$. Moreover, when setting the step size $\gamma\simeq n^{-\beta}$ with an appropriate $\beta\in(0,1)$, we can show that the moments of all other terms except for $W$ in the r.h.s. of (39) are small (see Theorem 3 for more details). Hence, using the same sequence $\{\xi_{k}\}_{k\in\mathbb{N}}$ of noise variables in (36) yields an estimator $\bar{\theta}_{n}^{(RR)}$, such that its leading component of the variance still equals $W$. Hence, using the Richardson-Romberg procedure will increase only the second-order (w.r.t. $n$) components of the variance. At the same time, using different random sequences $\{\xi_{k}\}_{k\in\mathbb{N}}$ and $\{\xi^{\prime}_{k}\}_{k\in\mathbb{N}}$ for $\bar{\theta}_{n}^{(\gamma)}$ and $\bar{\theta}_{n}^{(2\gamma)}$ increase the leading component of the MSE by a constant factor. Hence, it is preferable to use synchronous noise construction as introduced in (36).

Section 3.1 implies the following improved bound on the bias of $\bar{\theta}_{n}^{(RR)}$:

###### Proposition \theproposition .

Assume A 1, A 2, A 3($6$), and C 1($6$). Then, for any $\gamma\in(0,1/(\operatorname{L}\mathsf{C}_{\operatorname{step},6})]$, and any initial distribution $\nu$ on $\mathbb{R}^{d}$, it holds that

$$ $\mathsf{E}_{\nu}[\bar{\theta}_{n}^{(RR)}]=\theta^{\star}+B_{3}\gamma^{3/2}+ \mathcal{R}_{3}(\theta_{0}-\theta^{\star},\gamma,n)\,,$ (40) $$

where $B_{3}\in\mathbb{R}^{d}$ is a vector such that $\|B_{3}\|\leq\mathsf{C}_{1}$, and

$$ $\|\mathcal{R}_{3}(\theta_{0}-\theta^{\star},\gamma,n)\|\lesssim\frac{\mathrm{e }^{-\gamma\mu(n+1)/2}}{n\gamma\mu}\left(\mathsf{E}_{\nu}^{1/2}\bigl{[}\|\theta _{0}-\theta^{\star}\|^{2}\bigr{]}+\frac{\sqrt{\gamma}\tau_{2}}{\sqrt{\mu}} \right)\,.$ $$

The proof of Section 4 is provided in Appendix A. This result is a simple consequence of Section 4, since the linear in $\gamma$ component of the bias $\gamma\Delta_{1}$ from (24) cancels out when computing $\bar{\theta}_{n}^{(RR)}$. We are now ready to formulate the main result for the Richardson-Romberg estimate $\bar{\theta}_{n}^{(RR)}$.

###### Theorem 2 .

Assume A 1, A 2, A 3($6$), and C 1($6$). Then for any $\gamma\in(0,1/(\operatorname{L}\mathsf{C}_{\operatorname{step},6})]$, initial distribution $\nu$ and $n\in\mathbb{N}$, the Richardson-Romberg estimator $\bar{\theta}_{n}^{(RR)}$ defined in (37) satisfies

$$ $\displaystyle\mathsf{E}_{\nu}^{1/2}[\|\operatorname{H}^{\star}(\bar{\theta}_{n }^{(RR)}-\theta^{\star})\|^{2}]$ $\displaystyle\lesssim\frac{\sqrt{\operatorname{Tr}{\Sigma_{\varepsilon}^{\star }}}}{n^{1/2}}+\frac{\mathsf{C}_{\operatorname{RR},1}\gamma^{1/2}}{n^{1/2}}+ \frac{\mathsf{C}_{\operatorname{RR},2}}{\gamma^{1/2}n}+\mathsf{C}_{ \operatorname{RR},3}\gamma^{3/2}+\frac{\mathsf{C}_{\operatorname{RR},4}\gamma} {n^{1/2}}$ $\displaystyle\qquad\qquad+\mathcal{R}_{4}(n,\gamma,\|\theta_{0}-\theta^{\star} \|)\,,$ $$

where the constants $\mathsf{C}_{\operatorname{RR},1}$ to $\mathsf{C}_{\operatorname{RR},4}$ are defined in (68) and

$$ $\mathcal{R}_{4}(n,\gamma,\|\theta_{0}-\theta^{\star}\|)=\frac{\operatorname{L} (1-\gamma\mu)^{(n+1)/2}}{n\gamma\mu}\\ \times\left(\mathsf{E}_{\nu}^{1/2}[\|\theta_{0}-\theta^{\star}\|^{6}]+\mathsf{ E}_{\nu}^{1/2}[\|\theta_{0}-\theta^{\star}\|^{4}]+\mathsf{E}_{\nu}^{1/2}[\| \theta_{0}-\theta^{\star}\|^{2}]+\frac{\mathsf{D}_{\operatorname{last},4} \gamma\tau_{4}^{2}}{\mu}\right)\,.$ $$

Proof of Theorem 2 is provided in Appendix C. Similarly to Theorem 1, we can optimize the above bound setting $\gamma$ depending upon $n$.

###### Corollary \thecorollary .

Under the assumptions of Theorem 2, by setting $\gamma=n^{-1/2}$ with $n\geq\left(\operatorname{L}\mathsf{C}_{\operatorname{step},6}\right)^{2}$, it holds that

$$ ${\textstyle\mathsf{E}^{1/2}_{\nu}[\|\operatorname{H}^{\star}(\bar{\theta}_{n}^ {(RR)}-\theta^{\star})\|^{2}]}\lesssim\frac{\sqrt{\operatorname{Tr}{\Sigma_{ \varepsilon}^{\star}}}}{n^{1/2}}+\mathcal{O}\bigl{(}n^{-3/4}\bigr{)}+\mathcal{ R}_{4}(n,1/\sqrt{n},\|\theta_{0}-\theta^{\star}\|)\,.$ (41) $$

###### Proposition \theproposition .

###### Theorem 2 .

###### Corollary \thecorollary .

#### Discussion.

Note that the result of Section 4 is a counterpart of (5) with $\delta=1/4$. This decay rate of the second order term is the same as for the Root-SGD algorithm, and in general can not be improved, see Li et al. (). However, the assumptions of Theorem 2 are stronger compared to the ones imposed by Li et al. (). In particular, in A 2 we require that $f$ is $4$ times continuously differentiable and uniformly bounded. At the same time, Li et al. () impose Lipschitz continuity of the Hessian of $f$, which is essentially equivalent to bounded $3$-rd derivative of $f$. Our proof of Theorem 2 essentially relies on the $4$-th order Taylor expansion, and it is not clear, if this assumption can be relaxed. We leave further investigations of this question for future research.

Now we aim to generalize the previous result for the $p$-th moment bounds with $p\geq 2$. The key technical element of our proof for the $p$-th moment bound is the following statement, which can be viewed as a version of Rosenthal’s inequality (Rosenthal, ; Pinelis, ).

###### Proposition \theproposition .

Let $p\geq 2$ and assume A 1, A 2, A 3($2p$), and C 1($2p$). Then for any $\gamma\in(0,1/(\operatorname{L}\mathsf{C}_{\operatorname{step},2p})]$, it holds that

$$ ${\textstyle\mathsf{E}^{1/p}_{\pi_{\gamma}}\bigl{[}\|\sum_{k=0}^{n-1}\{\psi( \theta_{k}^{(\gamma)})-\pi_{\gamma}(\psi)\|^{p}\bigr{]}}\lesssim\frac{ \operatorname{L}\mathsf{D}_{\operatorname{last},2p}p\tau^{2}_{2p}\sqrt{n\gamma }}{\mu^{3/2}}+\frac{\operatorname{L}\mathsf{D}_{\operatorname{last},2p}\tau_{2 p}}{\mu^{2}}\,,$ (42) $$

where $\psi$ is defined in (38).

###### Proposition \theproposition .

#### Discussion.

The proof of Section 4 is provided in Section D.1. It is important to acknowledge that there are numerous Rosenthal-type inequalities for dependent sequences in the literature. Section 4 can be viewed as an analogue to the classical Rosenthal inequality for strongly mixing sequences, see (Rio, , Theorem 6.3). However, it should be emphasized that the Markov chain $\{\theta_{k}^{(\gamma)}\}_{k\in\mathbb{N}}$ is geometrically ergodic under the assumptions A 1-A 3($p$) only in sense of the weighted Wasserstein semi-metric $\mathbf{W}_{{c}}(\xi,\xi^{\prime})$ with respect to a cost function $c$ defined in (16). As a result, the sequence $\{\theta_{k}^{(\gamma)}\}_{k\in\mathbb{N}}$ does not necessarily satisfy strong mixing conditions. At the same time, $\{\theta_{k}^{(\gamma)}\}_{k\in\mathbb{N}}$ satisfies the $\tau$-mixing condition, see Merlevède et al. (). However, the considered function $\psi(\theta)$ is quadratic, which makes the respective result of (Merlevède et al., , Theorem 1) inapplicable. Similar Rosenthal-type inequalities have been explored in (Durmus et al., ), but in Section 4 we obtain the bound with tighter dependence of the right-hand side upon $\gamma$.

Below we provide the $p$-th moment bound together with its corollary for the step size $\gamma$ optimized with respect to $n$.

###### Theorem 3 .

Let $p\geq 2$ and assume A 1, A 2, A 3($3p$), and C 1($3p$). Then for any step size $\gamma\in(0,1/(\operatorname{L}\mathsf{C}_{\operatorname{step},3p})]$, initial distribution $\nu$, and $n\in\mathbb{N}$, the estimator $\bar{\theta}_{n}^{(RR)}$ defined in (37) satisfies

$$ $\begin{split}\mathsf{E}_{\nu}^{1/p}[\|\operatorname{H}^{\star}(\bar{\theta}_{n }^{(RR)}-\theta^{\star})\|^{p}]&\lesssim\frac{\sqrt{\operatorname{Tr}{\Sigma_{ \varepsilon}^{\star}}}p^{1/2}}{n^{1/2}}+\frac{p\tau_{p}}{n^{1-1/p}}+\frac{ \mathsf{C}_{\operatorname{RR},5}}{n\gamma^{1/2}}+\frac{\mathsf{C}_{ \operatorname{RR},6}\gamma^{1/2}}{n^{1/2}}+\mathsf{C}_{\operatorname{RR},7} \gamma^{3/2}\\ &\qquad\qquad+\frac{\mathsf{C}_{\operatorname{RR},8}}{n}+\mathcal{R}_{5}(n, \gamma,\|\theta_{0}-\theta^{\star}\|)\,,\end{split}$ (43) $$

where constants $\mathsf{C}_{\operatorname{RR},5}$ to $\mathsf{C}_{\operatorname{RR},8}$ are defined in (95) and

$$ $\textstyle\mathcal{R}_{5}(n,\gamma,\|\theta_{0}-\theta^{\star}\|)=(1-\gamma\mu )^{(n+1)/2}C_{f,p}\bigl{(}\mathsf{E}_{\nu}^{1/p}\bigl{[}\|\theta_{0}-\theta^{ \star}\|^{p}\bigr{]}+\mathsf{E}_{\nu}^{1/p}\bigl{[}\|\theta_{0}-\theta^{\star} \|^{2p}\bigr{]}+\mathsf{E}_{\nu}^{1/p}\bigl{[}\|\theta_{0}-\theta^{\star}\|^{3 p}\bigr{]}\bigr{)}\,,$ $$

where the constant $C_{f,p}$ can be traced from (96).

###### Corollary \thecorollary .

Under the assumptions of Theorem 3, by setting $\gamma=n^{-1/2}$ with $n\geq\left(\operatorname{L}\mathsf{C}_{\operatorname{step},3p}\right)^{2}$, it holds that

$$ $\textstyle\mathsf{E}^{1/p}_{\nu}[\|\operatorname{H}^{\star}(\bar{\theta}_{n}^{ (RR)}-\theta^{\star})\|^{p}]\lesssim\frac{\sqrt{\operatorname{Tr}{\Sigma_{ \varepsilon}^{\star}}}p^{1/2}}{n^{1/2}}+\mathcal{O}\bigl{(}n^{-3/4}\bigr{)}+ \mathcal{R}_{5}(n,1/\sqrt{n},\|\theta_{0}-\theta^{\star}\|)\,.$ (44) $$

###### Theorem 3 .

###### Corollary \thecorollary .

#### Discussion.

Proof of Theorem 3 is provided in Appendix D. Note that the result above is a direct generalization of Theorem 2, which reveals the same optimal scaling of the step size $\gamma$ with respect to $n$. To the best of our knowledge, this is the first analysis of a first-order method, which provides a bound for the second-order term of order $\mathcal{O}\bigl{(}n^{-3/4}\bigr{)}$ while keeping the precise leading term related to the minimax-optimal covariance matrix $\operatorname{H}^{\star}$. Such results were previously known only for the setting of linear stochastic approximation (LSA), which corresponds to the case of quadratic function $f$ in the initial minimization problem (1), see Durmus et al. (). In such a case, no bias occurs: $\int_{\mathbb{R}^{d}}\theta\pi_{\gamma}(\mathrm{d}\theta)=\theta^{\star}$, and Polyak-Ruppert averaging with the step size $\gamma\simeq n^{-1/2}$ allows for the same scaling of the remainder terms in $n$, as in (44). This result can be found, for example, in Durmus et al. (, Theorem 1). Hence, the main result of Theorem 3 is as follows: Richardson-Romberg extrapolation applied to strongly convex minimization problems allows to restore the $p$-th moment error moment bounds from the LSA setting.

## 5 Conclusion

In this paper, we have demonstrated that Polyak-Ruppert averaged SGD iterates with a constant step size achieve optimal root-MSE and $p$-th moment error bounds. More precisely, we have shown that these bounds admit both the sharp, optimal leading term, which aligns with the optimal covariance matrix $\Sigma_{\infty}$, and a second-order term of order $\mathcal{O}(n^{-3/4})$, which is known to be generally non-improvable. directions for future research include, firstly, generalizing the proposed algorithm to the setting of dependent noise sequences $\{\xi_{k}\}_{k\in\mathbb{N}}$ in the stochastic gradients (1). Another natural question is to study the properties of $\bar{\theta}_{n}^{(RR)}$ under relaxed assumptions on $f$. In particular, it would be interesting to remove additional smoothness assumptions on $f$ (bounded $3$-rd and $4$-th derivatives), and to relax the strong convexity condition A 1. One more research direction is to study a relation between the parameter $\delta$ in (5) and rates in the corresponding Berry-Esseen type results, following the technique of Shao & Zhang ().

## Appendix A Proof of Section 3.1 and Section 4

### A.1 Proof of Section 3.1

Consider the synchronous coupling construction defined by the recursions

$$ $\begin{split}\theta^{(\gamma)}_{k+1}&=\theta^{(\gamma)}_{k}-\gamma\nabla F( \theta^{(\gamma)}_{k},\xi_{k+1})\,,\quad\theta^{(\gamma)}_{0}=\theta\in\mathbb {R}^{d}\,,\\ \tilde{\theta}^{(\gamma)}_{k+1}&=\tilde{\theta}^{(\gamma)}_{k}-\gamma\nabla F( \tilde{\theta}^{(\gamma)}_{k},\xi_{k+1})\,,\quad\tilde{\theta}^{(\gamma)}_{0}= \tilde{\theta}\in\mathbb{R}^{d}\,.\end{split}$ (45) $$

The pair $(\theta^{(\gamma)}_{k},\tilde{\theta}^{(\gamma)}_{k})_{k\in\mathbb{N}}$ defines a Markov chain with the Markov kernel ${\rm K}_{\gamma}(\cdot,\cdot)$, which is a coupling kernel of $({\rm Q}_{\gamma},{\rm Q}_{\gamma})$. From now on we omit an upper index $(\gamma)$ and write simply $(\theta_{k},\tilde{\theta}_{k})_{k\in\mathbb{N}}$. Applying now A 3($2$), for $\gamma\leq 2/\operatorname{L}$, we get that

$$ $\displaystyle\mathsf{E}_{\nu}[\|\theta_{k+1}-\tilde{\theta}_{k+1}\|^{2}|{ \mathcal{F}}_{k}]=\mathsf{E}[\|\theta_{k}-\tilde{\theta}_{k}-\gamma(\nabla F( \theta_{k},\xi_{k+1})-\nabla F(\tilde{\theta}_{k},\xi_{k+1}))\|^{2}|{\mathcal{ F}}_{k}]$ $\displaystyle\qquad\qquad=\|\theta_{k}-\tilde{\theta}_{k}\|^{2}+\gamma^{2} \mathsf{E}[\|\nabla F(\theta_{k},\xi_{k+1})-\nabla F(\tilde{\theta}_{k},\xi_{k +1})\|^{2}|{\mathcal{F}}_{k}]$ $\displaystyle\qquad\qquad\qquad\qquad-2\gamma\langle\nabla f(\theta_{k})- \nabla f(\tilde{\theta}_{k}),\theta_{k}-\tilde{\theta}_{k}\rangle$ $\displaystyle\qquad\qquad\leq(1-\gamma\mu)\|\theta_{k}-\tilde{\theta}_{k}\|^{2},$ (46) $$

where in the last inequality we additionally used $(1-2\gamma\mu(1-\gamma\operatorname{L}/2))\leq 1-\gamma\mu$. Similarly, for a cost function $c$ defined in (16), we get using Hölder’s and Minkowski’s inequalities, that for any $r\in\mathbb{N}$

$$ $\mathsf{E}[c(\theta_{k+r},\tilde{\theta}_{k+r})|{\mathcal{F}}_{k}]\leq\mathsf{ E}^{1/2}[\|\theta_{k+r}-\tilde{\theta}_{k+r}\|^{2}|{\mathcal{F}}_{k}]\bigl{(} \mathsf{E}^{1/2}[\|\theta_{k+r}-\theta^{\star}\|^{2}|{\mathcal{F}}_{k}]\\ +\mathsf{E}^{1/2}[\|\tilde{\theta}_{k+r}-\theta^{\star}\|^{2}|{\mathcal{F}}_{k }]+\frac{2^{3/2}\gamma^{1/2}\tau_{2}}{\mu^{1/2}}\bigr{)}\,.$ $$

Combining the above inequalities and using (63), we obtain

$$ $\displaystyle\mathsf{E}[c(\theta_{k+r},\theta^{\prime}_{k+r})|{\mathcal{F}}_{k}]$ $\displaystyle\leq(1-\gamma\mu)^{r/2}\|\theta_{k}-\tilde{\theta}_{k}\|\bigl{(}( 1-\gamma\mu)^{r/2}(\|\theta_{k}-\theta^{\star}\|+\|\tilde{\theta}_{k}-\theta^{ \star}\|)+\frac{2^{5/2}\gamma^{1/2}\tau_{2}}{\mu^{1/2}}\bigr{)}$ $\displaystyle\leq 2(1-\gamma\mu)^{r/2}c(\theta_{k},\theta^{\prime}_{k})\,.$ $$

Note that $2(1-\gamma\mu)^{r/2}\leq 2$ for any $r\leq m(\gamma)-1$ and $2(1-\gamma\mu)^{m(\gamma)/2}\leq 1/2$. Hence, applying the theorem Douc et al. (, Theorem 20.3.4), we obtain that the Markov kernel ${\rm Q}_{\gamma}$ admits a unique invariant distribution $\pi_{\gamma}$. Moreover,

$$ $\mathbf{W}_{{c}}(\nu{\rm Q}_{\gamma}^{k},\pi_{\gamma})\leq 2(1/2)^{\lfloor k/m (\gamma)\rfloor}\mathbf{W}_{{c}}(\nu,\pi_{\gamma})\,.$ (47) $$

It remains to note that $(1/2)^{\lfloor k/m(\gamma)\rfloor}\leq 2(1/2)^{k/m(\gamma)}$.

### A.2 Proof of Section 3.1

We begin with proving (22) and (23). First we introduce some additional notations. Under assumptions A 1 – A 3($2$), we define a matrix-valued function $\mathcal{C}(\theta):\mathbb{R}^{d}\rightarrow\mathbb{R}^{d\times d}$ as

$$ $\mathcal{C}(\theta)=\mathsf{E}[{\varepsilon_{1}}(\theta)^{\otimes 2}]\,.$ (48) $$

The result below is essentially based on an appropriate modification of the bounds presented in Dieuleveut et al. (, Lemma 18). A careful inspection of the respective proof reveals that we do not need specific assumptions for $\mathcal{C}(\theta)$ defined in (48), instead we use Section A.2. For completeness, we present the respective result below:

###### Lemma \thelemma .

Assume A 1, A 2, A 3($6$), and C 1($6$). Then, for any $\gamma\in(0,1/(\operatorname{L}\mathsf{C}_{\operatorname{step},6})]$, it holds

$$ $\bar{\theta}_{\gamma}-\theta^{\star}=-(\gamma/2)\{\operatorname{H}^{\star}\}^{ -1}\{\nabla^{3}f(\theta^{\star})\}\mathbf{T}\mathcal{C}(\theta^{\star})+B_{1} \gamma^{3/2}\,,$ (49) $$

where $\bar{\theta}_{\gamma}$ is defined in (22), $\mathcal{C}(\theta)$ is defined in (48), and $B_{1}\in\mathbb{R}^{d}$ satisfies $\|B_{1}\|\leq\mathsf{C}_{1}$, where

$$ $\mathsf{C}_{1}=\left(\frac{\operatorname{L}^{2}\mathsf{D}_{\operatorname{last} ,2}}{\sqrt{\mu}}+\frac{\operatorname{L}\sqrt{\mathsf{D}_{\operatorname{last},2 }}}{\sqrt{\mu}}\right)\frac{\tau_{2}^{2}}{\mu}+\frac{\operatorname{L}}{\mu}\,.$ (50) $$

Moreover,

$$ $\bar{\Sigma}_{\gamma}=\gamma\mathbf{T}\mathcal{C}(\theta^{\star})+B_{2}\gamma^ {3/2}\,,$ (51) $$

where the operator $\mathbf{T}:\mathbb{R}^{d\times d}\to\mathbb{R}^{d\times d}$ is defined by the relation

$$ $\operatorname{vec}\left(\mathbf{T}A\right)=(\operatorname{H}^{\star}\otimes \mathrm{I}+\mathrm{I}\otimes\operatorname{H}^{\star})^{-1}\operatorname{vec} \left(A\right)$ $$

for any matrix $A\in\mathbb{R}^{d\times d}$, and $B_{2}\in\mathbb{R}^{d\times d}$ is a matrix, such that $\|B_{2}\|\leq\mathsf{C}_{1}$.

###### Proof.

Let $(\theta_{k}^{(\gamma)})_{k\in\mathbb{N}}$ be a recurrence defined in (11) with initial distribution $\theta_{0}\sim\pi_{\gamma}$. Moreover, we assume that $\theta_{0}$ is independent of the noise variables $(\xi_{k})_{k\geq 1}$. First, applying a third-order Taylor expansion of $\nabla f(\theta)$ around $\theta^{\star}$, for any $\theta\in\mathbb{R}^{d}$, we obtain

$$ $\nabla f(\theta)=\operatorname{H}^{\star}(\theta-\theta^{\star})+(1/2)\{\nabla ^{3}f(\theta^{\star})\}(\theta-\theta^{\star})^{\otimes 2}+G(\theta)\,,$ (52) $$

where $G(\theta)$ is defined in (38) and writes as

$$ $G(\theta)=\frac{1}{6}\left(\int_{0}^{1}\nabla^{4}f(t\theta^{\star}+(1-t)\theta )\,dt\right)(\theta-\theta^{\star})^{\otimes 3}\,.$ $$

Thus, using A 2,

$$ $\|G(\theta)\|\lesssim\operatorname{L}_{4}\|\theta-\theta^{\star}\|^{3}\,.$ $$

Thus, integrating (52) with respect to $\pi_{\gamma}$, we get from (52) that

$$ $\operatorname{H}^{\star}(\bar{\theta}_{\gamma}-\theta^{\star})+(1/2)\{\nabla^{ 3}f(\theta^{\star})\}\left[\int_{\mathbb{R}^{d}}(\theta-\theta^{\star})^{ \otimes 2}\pi_{\gamma}(\mathrm{d}\theta)\right]=-\int_{\mathbb{R}^{d}}G(\theta )\pi_{\gamma}(\mathrm{d}\theta)\,.$ (53) $$

Now we need to provide an explicit expression for the covariance matrix

$$ $\bar{\Sigma}_{\gamma}=\int_{\mathbb{R}^{d}}(\theta-\theta^{\star})^{\otimes 2} \pi_{\gamma}(\mathrm{d}\theta)\,.$ (54) $$

Using the recurrence (11), we obtain that

$$ $\theta_{1}-\theta^{\star}=(\mathrm{I}-\gamma\operatorname{H}^{\star})(\theta_{ 0}-\theta^{\star})-\gamma{\varepsilon_{1}}(\theta_{0})-\gamma\eta(\theta_{0})\,,$ $$

where the function $\eta(\cdot)$ is defined in (28). Hence, taking second moment w.r.t. $\pi_{\gamma}$ from both sides, we get that

$$ $\bar{\Sigma}_{\gamma}=(\mathrm{I}-\gamma\operatorname{H}^{\star})\bar{\Sigma}_ {\gamma}(\mathrm{I}-\gamma\operatorname{H}^{\star})+\gamma^{2}\int_{\mathbb{R} ^{d}}\mathcal{C}(\theta)\pi_{\gamma}(\mathrm{d}\theta)+\gamma^{2}\int_{\mathbb {R}^{d}}\{\eta(\theta)\}^{\otimes 2}\pi_{\gamma}(\mathrm{d}\theta)\\ -\gamma\int_{\mathbb{R}^{d}}\left[(\mathrm{I}-\gamma\operatorname{H}^{\star})( \theta-\theta^{\star})\{\eta(\theta)\}^{\top}+\eta(\theta)(\theta-\theta^{ \star})^{\top}(\mathrm{I}-\gamma\operatorname{H}^{\star})\right]\pi_{\gamma}( \mathrm{d}\theta)\,.$ (55) $$

In the above equation $\mathcal{C}(\theta)$ is defined in (48), and we additionally used that $\mathsf{E}\bigl{[}{\varepsilon_{1}}(\theta_{0})|{\mathcal{F}}_{0}\bigr{]}=0$. Moreover, (48) together with C 1($6$) implies that

$$ $\int_{\mathbb{R}^{d}}\mathcal{C}(\theta)\pi_{\gamma}(\mathrm{d}\theta)= \mathcal{C}(\theta^{\star})+B\gamma^{1/2}\,,$ $$

where $B\in\mathbb{R}^{d\times d}$ satisfies $\|B\|\leq\mathsf{C}_{2}$. Thus, from (55) together with C 1($6$) we obtain that $\bar{\Sigma}_{\gamma}$ is a solution to the matrix equation

$$ $\operatorname{H}^{\star}\bar{\Sigma}_{\gamma}+\bar{\Sigma}_{\gamma} \operatorname{H}^{\star}-\gamma\operatorname{H}^{\star}\bar{\Sigma}_{\gamma} \operatorname{H}^{\star}=\gamma\mathcal{C}(\theta^{\star})+B\gamma^{3/2}\,,$ $$

which can be written using the vectorization operation as

$$ $\operatorname{vec}\left(\bar{\Sigma}_{\gamma}\right)=\gamma(\operatorname{H}^{ \star}\otimes\mathrm{I}+\mathrm{I}\otimes\operatorname{H}^{\star}-\gamma \operatorname{H}^{\star}\otimes\operatorname{H}^{\star})^{-1}\operatorname{vec }\left(\mathcal{C}(\theta^{\star})\right)+\gamma^{3/2}(\operatorname{H}^{\star }\otimes\mathrm{I}+\mathrm{I}\otimes\operatorname{H}^{\star}-\gamma \operatorname{H}^{\star}\otimes\operatorname{H}^{\star})^{-1}\operatorname{vec }\left(B\right)\,.$ $$

Now we check that the latter operator $\operatorname{H}^{\star}\otimes\mathrm{I}+\mathrm{I}\otimes\operatorname{H}^{
\star}-\gamma\operatorname{H}^{\star}\otimes\operatorname{H}^{\star}$ is indeed invertible for $\gamma\in(0,2/\operatorname{L})$. Moreover, assumption A 1 guarantees that the symmetric matrix $\operatorname{H}^{\star}$ is non-degenerate and positive-definite. Let $u_{1},\ldots,u_{d}\in\mathbb{R}^{d}$ and $\lambda_{1}\geq\lambda_{2}\geq\ldots\geq\lambda_{d}\geq\mu>0$ be its eigenvectors and eigenvalues, respectively. Then we notice that

$$ $\operatorname{H}^{\star}\otimes\mathrm{I}+\mathrm{I}\otimes\operatorname{H}^{ \star}-\gamma\operatorname{H}^{\star}\otimes\operatorname{H}^{\star}= \operatorname{H}^{\star}\otimes(\mathrm{I}-(\gamma/2)\operatorname{H}^{\star}) +(\mathrm{I}-(\gamma/2)\operatorname{H}^{\star})\otimes\operatorname{H}^{\star }\,.$ $$

Hence, the latter operator is also diagonalizable in the orthogonal basis $u_{i}\otimes u_{j}$ in $\mathbb{R}^{d^{2}}$ with the respective eigenvalues being equal to $\lambda_{i}(1-(\gamma/2)\lambda_{j})+\lambda_{j}(1-(\gamma/2)\lambda_{i})$. Set now

$$ $\begin{split}S&=\operatorname{H}^{\star}\otimes\mathrm{I}+\mathrm{I}\otimes \operatorname{H}^{\star}\in\mathbb{R}^{d^{2}\times d^{2}}\\ R&=\operatorname{H}^{\star}\otimes\operatorname{H}^{\star}\in\mathbb{R}^{d^{2} \times d^{2}}\,.\end{split}$ (56) $$

Then it is easy to observe that

$$ $\displaystyle(S-\gamma R)^{-1}=S^{-1}+S^{-1}\sum_{k=1}^{\infty}\gamma^{k}(RS^{ -1})^{k}\,,$ $$

provided that $\gamma\|RS^{-1}\|\leq 1$. Sice $R$ and $S$ are diagonalizable in the same orthogonal basis $\{u_{i}\otimes u_{j}\}_{1\leq i,j\leq d}$ with the eigenvalues $\lambda_{i}\lambda_{j}$ and $\lambda_{i}+\lambda_{j}$, respectively, the condition $\gamma\|RS^{-1}\|\leq 1$ holds provided that $\gamma\leq 2/\operatorname{L}$. Hence, for $\gamma\leq 1/\operatorname{L}$, it holds that

$$ $(\operatorname{H}^{\star}\otimes\mathrm{I}+\mathrm{I}\otimes\operatorname{H}^{ \star}-\gamma\operatorname{H}^{\star}\otimes\operatorname{H}^{\star})^{-1}=( \operatorname{H}^{\star}\otimes\mathrm{I}+\mathrm{I}\otimes\operatorname{H}^{ \star})^{-1}+D\,,$ $$

where $D\in\mathbb{R}^{d^{2}\times d^{2}}$ satisfies

$$ $\|D\|\lesssim\gamma\|S\|^{-1}\|RS^{-1}\|\lesssim\frac{\gamma\operatorname{L}}{ \mu}\,.$ $$

Combining the above bounds in (53), we arrive at the expansion formula (49).
∎

We now state an auxiliary lemma about the function $\mathcal{C}(\theta)$ from (48).

###### Lemma \thelemma .

Assume A 1, A 2, A 3($2$), and C 1($2$). Then, for any $\gamma\in(0,1/(\operatorname{L}\mathsf{C}_{\operatorname{step},2})]$, it holds

$$ $\|\int_{\mathbb{R}^{d}}\mathcal{C}(\theta)\pi_{\gamma}(\mathrm{d}\theta)- \mathcal{C}(\theta^{\star})\|\leq\mathsf{C}_{2}\gamma^{1/2}\,,$ $$

where the constant $\mathsf{C}_{2}$ is given by

$$ $\mathsf{C}_{2}=\left(\frac{\operatorname{L}^{2}\mathsf{D}_{\operatorname{last} ,2}}{\sqrt{\mu}}+\frac{\operatorname{L}\sqrt{\mathsf{D}_{\operatorname{last},2 }}}{\sqrt{\mu}}\right)\tau_{2}^{2}\,.$ (57) $$

###### Proof.

Recall that

$$ $\displaystyle{\varepsilon_{1}}(\theta)=\nabla F(\theta,\xi_{1})-\nabla f( \theta)\,.$ $$

Hence, using the definition of $\mathcal{C}(\theta)$ in (48), we get

$$ $\displaystyle\mathcal{C}(\theta)-\mathcal{C}(\theta^{\star})=\mathsf{E}_{\xi_{ 1}\sim\mathbb{P}_{\xi}}[({\varepsilon_{1}}(\theta)-{\varepsilon_{1}}(\theta^{ \star}))({\varepsilon_{1}}(\theta)-{\varepsilon_{1}}(\theta^{\star}))^{T}]+ \mathsf{E}_{\xi_{1}\sim\mathbb{P}_{\xi}}[{\varepsilon_{1}}(\theta^{\star})({ \varepsilon_{1}}(\theta)-{\varepsilon_{1}}(\theta^{\star}))^{T}]$ $\displaystyle+\mathsf{E}_{\xi_{1}\sim\mathbb{P}_{\xi}}[({\varepsilon_{1}}( \theta)-{\varepsilon_{1}}(\theta^{\star})){\varepsilon_{1}}(\theta^{\star})^{T }].$ $$

Using A 3($2$), we obtain

$$ $\mathsf{E}_{\xi}[\|{\varepsilon_{1}}(\theta)-{\varepsilon_{1}}(\theta^{\star}) \|^{2}]\lesssim\operatorname{L}\langle\nabla f(\theta)-\nabla f(\theta^{\star} ),\theta-\theta^{\star}\rangle-\|\nabla f(\theta)-\nabla f(\theta^{\star})\|^{ 2}\lesssim\operatorname{L}^{2}\|\theta-\theta^{\star}\|^{2}.$ $$

Hence, combining the previous inequalities and using Hölder’s inequality, we obtain for any $\theta\in\mathbb{R}^{d}$, that

$$ $\|\mathcal{C}(\theta)-\mathcal{C}(\theta^{\star})\|\lesssim\operatorname{L}^{2 }\|\theta-\theta^{\star}\|^{2}+\tau_{2}\operatorname{L}\|\theta-\theta^{\star}\|.$ $$

Applying now C 1($2$), we obtain

$$ $\|\int_{\mathbb{R}^{d}}\mathcal{C}(\theta)\pi_{\gamma}(\mathrm{d}\theta)- \mathcal{C}(\theta^{\star})\|\leq\int_{\mathbb{R}^{d}}\|\mathcal{C}(\theta)- \mathcal{C}(\theta^{\star})\|\pi_{\gamma}(\mathrm{d}\theta)\lesssim \operatorname{L}^{2}\frac{\mathsf{D}_{\operatorname{last},2}\gamma\tau_{2}^{2} }{\mu}+\tau_{2}\operatorname{L}\sqrt{\frac{\mathsf{D}_{\operatorname{last},2} \gamma\tau_{2}^{2}}{\mu}}.$ $$

We conclude the proof by noting that $\gamma\mu\leq 1$.
∎

Now we prove (24). We use synchronous coupling construction defined by the pair of recursions:

$$ $\displaystyle\theta_{k+1}=\theta_{k}-\gamma\nabla F(\theta_{k},\xi_{k+1}), \quad\theta_{0}\sim\nu$ $\displaystyle\tilde{\theta}_{k+1}=\tilde{\theta}_{k}-\gamma\nabla F(\tilde{ \theta}_{k},\xi_{k+1}),\quad\tilde{\theta}_{0}\sim\pi_{\gamma}\,.$ $$

Recall that the corresponding coupling kernel is denoted as ${\rm K}_{\gamma}(\cdot,\cdot)$. Then we obtain

$$ $\displaystyle\mathsf{E}_{\nu}[\bar{\theta}_{n}]-\theta^{\star}$ $\displaystyle=n^{-1}\sum_{k=n+1}^{2n}\mathsf{E}^{{\rm K}_{\gamma}}_{\nu,\pi_{ \gamma}}[\theta_{k}-\tilde{\theta}_{k}]+n^{-1}\sum_{k=n+1}^{2n}\mathsf{E}_{\pi _{\gamma}}[\tilde{\theta}_{k}-\theta^{\star}]$ $\displaystyle=n^{-1}\sum_{k=n+1}^{2n}\mathsf{E}^{{\rm K}_{\gamma}}_{\nu,\pi_{ \gamma}}[\theta_{k}-\tilde{\theta}_{k}]+(\bar{\theta}_{\gamma}-\theta^{\star})\,.$ $$

Using (46) and C 1($2$), we obtain

$$ $\displaystyle\|\mathsf{E}^{K}_{\nu,\pi_{\gamma}}[\{\theta_{k}-\tilde{\theta}_{ k}\}]\|$ $\displaystyle\leq(1-\gamma\mu)^{k/2}\{\mathsf{E}^{{\rm K}_{\gamma}}_{\nu,\pi_{ \gamma}}\|\theta_{0}-\tilde{\theta}_{0}\|^{2}\}^{1/2}$ $\displaystyle\leq(1-\gamma\mu)^{k/2}(\mathsf{E}_{\nu}^{1/2}\bigl{[}\|\theta_{0 }-\theta^{\star}\|^{2}\bigr{]}+\frac{\sqrt{2\gamma}\tau_{2}}{\sqrt{\mu}})\,.$ $$

Summing the above bounds for $k$ from $n+1$ to $2n$, we obtain (24).

###### Lemma \thelemma .

###### Proof.

###### Lemma \thelemma .

###### Proof.

### A.3 Proof of Section 4

Note that

$$ $\mathsf{E}_{\nu}[\bar{\theta}_{n}^{(RR)}-\theta^{\star}]=2\mathsf{E}_{\nu}[ \bar{\theta}_{n}^{\gamma}-\theta^{\star}]-\mathsf{E}_{\nu}[\bar{\theta}_{n}^{2 \gamma}-\theta^{\star}].$ $$

Applying (24), we obtain

$$ $\|\mathsf{E}_{\nu}[\bar{\theta}_{n}^{(RR)}-\theta^{\star}]\|\lesssim\mathsf{C} _{1}\gamma^{3/2}+\mathcal{R}_{3}(\theta_{0}-\theta^{\star},\gamma,n),$ (58) $$

where

$$ $\|\mathcal{R}_{3}(\theta_{0}-\theta^{\star},\gamma,n)\|\lesssim\frac{(1-\gamma \mu)^{(n+1)/2}}{n\gamma\mu}(\mathsf{E}_{\nu}^{1/2}\bigl{[}\|\theta_{0}-\theta^ {\star}\|^{2}\bigr{]}+\frac{\sqrt{\gamma}\tau_{2}}{\sqrt{\mu}})\,.$ (59) $$

## Appendix B Proof of Theorem 1

###### Theorem 4 (Version of Theorem 1 with explicit constants) .

Assume A 1, A 2, A 3($6$), and C 1($6$). Then for any $\gamma\in(0,1/(\operatorname{L}\mathsf{C}_{\operatorname{step},6})]$, $n\in\mathbb{N}$, and initial distribution $\nu$ on $\mathbb{R}^{d}$, the sequence of Polyak-Ruppert estimates (14) satisfies

$$ $\mathsf{E}_{\nu}^{1/2}[\|\operatorname{H}^{\star}(\bar{\theta}_{n}^{(\gamma)}- \theta^{\star})\|^{2}]\lesssim\frac{\sqrt{\operatorname{Tr}{\Sigma_{ \varepsilon}^{\star}}}}{\sqrt{n}}+\frac{\mathsf{C}_{2}}{\gamma^{1/2}n}+\mathsf {C}_{3}\gamma+\frac{\mathsf{C}_{4}\gamma^{1/2}}{n^{1/2}}+\mathcal{R}_{2}(n, \gamma,\|\theta_{0}-\theta^{\star}\|)\,,$ $$

where we have set

$$ $\begin{split}\mathsf{C}_{2}&=\mathsf{D}_{\operatorname{last},2}^{1/2}\tau_{2} \,,\quad\mathsf{C}_{3}=\frac{\operatorname{L}\mathsf{D}_{\operatorname{last},4 }\tau_{4}^{2}}{2\mu}\,,\quad\mathsf{C}_{4}=\operatorname{L}\mathsf{D}_{ \operatorname{last},2}^{1/2}\tau_{2}\,.\end{split}$ (60) $$

and the remainder term $\mathcal{R}_{2}(n,\gamma,\|\theta_{0}-\theta^{\star}\|)$ is given by

$$ $\mathcal{R}_{2}(n,\gamma,\|\theta_{0}-\theta^{\star}\|)=\frac{(1-\gamma\mu)^{( n+1)/2}}{\gamma n}\left(\frac{\operatorname{L}}{\mu}+1\right)\mathsf{E}_{\nu}^ {1/2}\bigl{[}\|\theta_{0}-\theta^{\star}\|^{2}\bigr{]}\\ +\frac{\operatorname{L}(1-\gamma\mu)^{n+1}}{2n\gamma\mu}\mathsf{E}^{1/2}_{\nu} \bigl{[}\|\theta_{0}-\theta^{\star}\|^{4}\bigr{]}\,.$ (61) $$

###### Proof.

Throughout the proof we omit upper index $(\gamma)$ both for the elements of the sequence $\{\theta_{k}^{(\gamma)}\}_{k\in\mathbb{N}}$ and Polyak-Ruppert averaged estimates $\bar{\theta}_{n}^{(\gamma)}$. Instead, we write simply $\theta_{k}$ and $\bar{\theta}_{n}$, respectively. Summing the recurrence (31), we obtain that

$$ $\operatorname{H}^{\star}(\bar{\theta}_{n}-\theta^{\star})=\frac{\theta_{n+1}- \theta^{\star}}{\gamma n}-\frac{\theta_{2n}-\theta^{\star}}{\gamma n}-\frac{1} {n}\sum_{k=n+1}^{2n}{\varepsilon_{k+1}}(\theta_{k})-\frac{1}{n}\sum_{k=n+1}^{2 n}\eta(\theta_{k})\,.$ (62) $$

Applying the $3$-rd order Taylor expansion with integral remainder, we get that

$$ $\nabla f(\theta_{k})=\operatorname{H}^{\star}(\theta_{k}-\theta^{\star})+\frac {1}{2}\left(\int_{0}^{1}\nabla^{3}f(t\theta^{\star}+(1-t)\theta_{k})\,dt\right )(\theta_{k}-\theta^{\star})^{\otimes 2}\,,$ $$

where $\nabla^{3}f(\cdot)\in\mathbb{R}^{d\times d\times d}$. Using A 2, we thus obtain that

$$ $\|\eta(\theta_{k})\|\leq\frac{1}{2}\operatorname{L}_{3}\|\theta_{k}-\theta^{ \star}\|^{2}\,.$ $$

Applying Minkowski’s inequality to the decomposition (33) and to the last term therein, we get

$$ $\displaystyle\mathsf{E}_{\nu}^{1/2}[\|\operatorname{H}^{\star}(\bar{\theta}_{n }-\theta^{\star})\|^{2}]\leq\frac{\mathsf{E}_{\nu}^{1/2}[\|\theta_{n+1}-\theta ^{\star}\|^{2}]}{\gamma n}+\frac{\mathsf{E}_{\nu}^{1/2}[\|\theta_{2n}-\theta^{ \star}\|^{2}]}{\gamma n}+\frac{1}{n}\mathsf{E}_{\nu}^{1/2}\bigl{[}\|\sum_{k=n+ 1}^{2n}{\varepsilon_{k+1}}(\theta_{k})\|^{2}\bigr{]}$ $\displaystyle+\frac{\operatorname{L}_{3}}{2n}\sum_{k=n+1}^{2n}\mathsf{E}_{\nu} ^{1/2}\bigl{[}\|\theta_{k}-\theta^{\star}\|^{4}\bigr{]}\,.$ $$

Applying C 1($2$), we obtain that for $\gamma\in(0;1/(\operatorname{L}\mathsf{C}_{\operatorname{step},2})]$ it holds that

$$ $\mathsf{E}_{\nu}{\|\theta_{k}-\theta^{\star}\|^{2}}\lesssim(1-\gamma\mu)^{k} \mathsf{E}_{\nu}\bigl{[}\|\theta_{0}-\theta^{\star}\|^{2}\bigr{]}+\frac{ \mathsf{D}_{\operatorname{last},2}\gamma\tau_{2}^{2}}{\mu}\,.$ (63) $$

Moreover, from $\gamma\in(0;1/(\operatorname{L}\mathsf{C}_{\operatorname{step},4})]$ it holds that

$$ $\mathsf{E}_{\nu}^{1/2}{\|\theta_{k}-\theta^{\star}\|^{4}}\lesssim(1-\gamma\mu) ^{k}\mathsf{E}^{1/2}_{\nu}\bigl{[}\|\theta_{0}-\theta^{\star}\|^{4}\bigr{]}+ \frac{\mathsf{D}_{\operatorname{last},4}\gamma\tau_{4}^{2}}{\mu}\,.$ (64) $$

Combining Appendix B with previous inequalities, we obtain

$$ $\mathsf{E}_{\nu}^{1/2}[\|\operatorname{H}^{\star}(\bar{\theta}_{n}-\theta^{ \star})\|^{2}]\lesssim\frac{\sqrt{\operatorname{Tr}{\Sigma_{\varepsilon}^{ \star}}}}{\sqrt{n}}+\frac{\mathsf{D}_{\operatorname{last},2}^{1/2}\tau_{2}}{ \gamma^{1/2}n}+\frac{\operatorname{L}\mathsf{D}_{\operatorname{last},4}\gamma \tau_{4}^{2}}{2\mu}+\frac{\operatorname{L}\mathsf{D}_{\operatorname{last},2}^{ 1/2}\gamma^{1/2}\tau_{2}}{\mu^{1/2}n^{1/2}}\\ +\frac{(1-\gamma\mu)^{(n+1)/2}}{\gamma n}\left(\frac{\operatorname{L}}{\mu}+1 \right)\mathsf{E}_{\nu}^{1/2}\bigl{[}\|\theta_{0}-\theta^{\star}\|^{2}\bigr{]} +\frac{\operatorname{L}(1-\gamma\mu)^{n+1}}{n\gamma\mu}\mathsf{E}_{\nu}^{1/2} \bigl{[}\|\theta_{0}-\theta^{\star}\|^{4}\bigr{]}\,,$ $$

and the result follows.
∎

###### Lemma \thelemma .

Assume A 1, A 3($2$), A 2, and C 1($2$). Then for any $\gamma\in(0;1/(\operatorname{L}\mathsf{C}_{\operatorname{step},2})]$ and any $n\in\mathbb{N}$, it holds

$$ $\mathsf{E}_{\nu}^{1/2}[\|\sum_{k=n+1}^{2n}\{{\varepsilon_{k+1}}(\theta_{k})-{ \varepsilon_{k+1}}(\theta^{\star})\}\|^{2}]\lesssim\frac{\operatorname{L} \mathsf{D}_{\operatorname{last},2}^{1/2}\sqrt{\gamma n}\tau_{2}}{\mu^{1/2}}+ \frac{\operatorname{L}(1-\gamma\mu)^{(n+1)/2}}{\gamma\mu}\mathsf{E}_{\nu}^{1/2 }\bigl{[}\|\theta_{0}-\theta^{\star}\|^{2}\bigr{]}\,.$ (65) $$

Moreover, under A 1, A 3($p$), A 2, and C 1($p$), for any $\gamma\in(0;1/(\operatorname{L}\mathsf{C}_{\operatorname{step},p})]$ and $n\in\mathbb{N}$ it holds that

$$ $\begin{split}\mathsf{E}_{\nu}^{1/p}[\|\sum_{k=n+1}^{2n}\{{\varepsilon_{k+1}}( \theta_{k})-{\varepsilon_{k+1}}(\theta^{\star})\}\|^{p}]&\lesssim\frac{ \operatorname{L}\mathsf{D}_{\operatorname{last},p}^{1/2}\sqrt{\gamma n}p\tau_{ p}}{\mu^{1/2}}\\ &\qquad+\frac{\operatorname{L}p(1-\gamma\mu)^{(n+1)/2}}{\mu^{1/2}\gamma^{1/2}} \mathsf{E}^{1/p}_{\nu}[\|\theta_{0}-\theta^{\star}\|^{p}]\,.\end{split}$ (66) $$

###### Proof.

Since $\{{\varepsilon_{k+1}}(\theta_{k})-{\varepsilon_{k+1}}(\theta^{\star})\}$ is a martingale-difference sequence with respect to ${\mathcal{F}}_{k}$, we have

$$ $\displaystyle\mathsf{E}_{\nu}\bigl{[}\|\sum_{k=n+1}^{2n}\{{\varepsilon_{k+1}}( \theta_{k})-{\varepsilon_{k+1}}(\theta^{\star})\}\|^{2}\bigr{]}$ $\displaystyle=\sum_{k=n+1}^{2n}\mathsf{E}_{\nu}\bigl{[}\|\{{\varepsilon_{k+1}} (\theta_{k})-{\varepsilon_{k+1}}(\theta^{\star})\}\|^{2}\bigr{]}\,.$ $$

where ${\varepsilon_{k+1}}(\theta^{\star})=\nabla F(\theta^{\star},\xi_{k+1})$ uses the same noise variable $\xi_{k+1}$ as $F(\theta_{k},\xi_{k+1})$.
Note that

$$ $\mathsf{E}_{\nu}[\|{\varepsilon_{k+1}}(\theta_{k})-{\varepsilon_{k+1}}(\theta^ {\star})\|^{2}]=\mathsf{E}_{\nu}[\|\nabla F(\theta_{k},\xi_{k+1})-\nabla F( \theta^{\star},\xi_{k+1})\|^{2}\\ \qquad-2\mathsf{E}_{\nu}\bigl{[}\langle\nabla F(\theta_{k},\xi_{k+1})-\nabla F (\theta^{\star},\xi_{k+1}),\nabla f(\theta_{k})-\nabla f(\theta^{\star}) \rangle\bigr{]}+\|\nabla f(\theta_{k})-\nabla f(\theta^{\star})\|^{2}]\,.$ $$

Using A 2, A 3($2$), and taking conditional expectation with respect to $\mathcal{F}_{k}$, we obtain

$$ $\displaystyle\mathsf{E}_{\nu}[\|{\varepsilon_{k+1}}(\theta_{k})-{\varepsilon_{ k+1}}(\theta^{\star})\|^{2}]$ $\displaystyle\leq\mathsf{E}_{\nu}[\operatorname{L}\langle\nabla f(\theta_{k})- \nabla f(\theta^{\star}),\theta_{k}-\theta^{\star}\rangle-\|\nabla f(\theta_{k })-\nabla f(\theta^{\star})\|^{2}]$ $\displaystyle\leq\operatorname{L}^{2}\mathsf{E}_{\nu}[\|\theta_{k}-\theta^{ \star}\|^{2}].$ $$

Thus, we obtain that

$$ $\displaystyle\mathsf{E}_{\nu}[\|\sum_{k=n+1}^{2n}\{{\varepsilon_{k+1}}(\theta_ {k})-{\varepsilon_{k+1}}(\theta^{\star})\}\|^{2}]\leq\operatorname{L}^{2}\sum_ {k=n+1}^{2n}\mathsf{E}_{\nu}[\|\theta_{k}-\theta^{\star}\|^{2}]\,,$ $$

and the statement (65) follows from the assumption C 1($2$). In order to prove (66), we apply Burkholder’s inequality Osekowski (, Theorem 8.6) and obtain

$$ $\displaystyle\mathsf{E}_{\nu}^{1/p}[\|\sum_{k=n+1}^{2n}\{{\varepsilon_{k+1}}( \theta_{k})-{\varepsilon_{k+1}}(\theta^{\star})\}\|^{p}]\leq p\mathsf{E}_{\nu} ^{1/p}\bigl{[}\bigl{(}\sum_{k=n+1}^{2n}\|{\varepsilon_{k+1}}(\theta_{k})-{ \varepsilon_{k+1}}(\theta^{\star})\|^{2}\bigr{)}^{p/2}\bigr{]}$ $\displaystyle\qquad\leq p\bigl{(}\sum_{k=n+1}^{2n}\mathsf{E}_{\nu}^{2/p}\bigl{ [}\|{\varepsilon_{k+1}}(\theta_{k})-{\varepsilon_{k+1}}(\theta^{\star})\|^{p} \bigr{]}\bigr{)}^{1/2}$ $\displaystyle\qquad\lesssim p\operatorname{L}\bigl{(}\sum_{k=n+1}^{2n}\mathsf{ E}^{2/p}_{\nu}\bigl{[}\|\theta_{k}-\theta^{\star}\|^{p}\bigr{]}\bigr{)}^{1/2}$ $\displaystyle\qquad\overset{(a)}{\lesssim}\frac{\operatorname{L}\mathsf{D}_{ \operatorname{last},p}^{1/2}\sqrt{\gamma n}p\tau_{p}}{\mu^{1/2}}+\frac{ \operatorname{L}p(1-\gamma\mu)^{(n+1)/2}}{\mu^{1/2}\gamma^{1/2}}\mathsf{E}^{1/ p}_{\nu}[\|\theta_{0}-\theta^{\star}\|^{p}]\,,$ $$

where in (a) we have additionally used C 1($p$).
∎

###### Theorem 4 (Version of Theorem 1 with explicit constants) .

###### Proof.

###### Lemma \thelemma .

###### Proof.

## Appendix C Proof of Theorem 2

Within this section we often use the definition of the function $\psi:\mathbb{R}^{d}\to\mathbb{R}^{d}$ from (38):

$$ $\psi(\theta)=(1/2)\nabla^{3}f(\theta^{*})(\theta-\theta^{\star})^{\otimes 2}$ (67) $$

###### Theorem 5 (Version of Theorem 2 with explicit constants) .

Assume A 1, A 2, A 3($6$), and C 1($6$). Then for any $\gamma\in(0,1/(\operatorname{L}\mathsf{C}_{\operatorname{step},6})]$, initial distribution $\nu$ and $n\in\mathbb{N}$, the Richardson-Romberg estimator $\bar{\theta}_{n}^{(RR)}$ defined in (37) satisfies

$$ $\displaystyle\mathsf{E}_{\nu}^{1/2}[\|\operatorname{H}^{\star}(\bar{\theta}_{n }^{(RR)}-\theta^{\star})\|^{2}]$ $\displaystyle\lesssim\frac{\sqrt{\operatorname{Tr}{\Sigma_{\varepsilon}^{\star }}}}{n^{1/2}}+\frac{\mathsf{C}_{\operatorname{RR},1}\gamma^{1/2}}{n^{1/2}}+ \frac{\mathsf{C}_{\operatorname{RR},2}}{\gamma^{1/2}n}+\mathsf{C}_{ \operatorname{RR},3}\gamma^{3/2}+\frac{\mathsf{C}_{\operatorname{RR},4}\gamma} {n^{1/2}}$ $\displaystyle\qquad\qquad+\mathcal{R}_{4}(n,\gamma,\|\theta_{0}-\theta^{\star} \|)\,,$ $$

where we have set

$$ $\begin{split}\mathsf{C}_{\operatorname{RR},1}&=\frac{\mathsf{D}_{\operatorname {last},4}\operatorname{L}\tau_{4}^{2}}{\mu^{3/2}}+\frac{\operatorname{L} \mathsf{D}_{\operatorname{last},2}^{1/2}\tau_{2}}{\mu^{1/2}}\,,\quad\mathsf{C} _{\operatorname{RR},2}=\frac{\mathsf{D}_{\operatorname{last},2}^{1/2}\tau_{2}} {\mu^{1/2}}\\ \mathsf{C}_{\operatorname{RR},3}&=\frac{\operatorname{L}\mathsf{D}_{ \operatorname{last},6}^{3/2}\tau_{6}^{3}}{\mu^{3/2}}+\mathsf{C}_{1}\,,\quad \mathsf{C}_{\operatorname{RR},4}=\frac{\mathsf{D}_{\operatorname{last},4} \operatorname{L}\tau_{4}^{2}}{\mu}\,,\end{split}$ (68) $$

$\mathsf{C}_{1}$ is defined in (50), and the remainder term $\mathcal{R}_{4}(n,\gamma,\|\theta_{0}-\theta^{\star}\|)$ is given by

$$ $\mathcal{R}_{4}(n,\gamma,\|\theta_{0}-\theta^{\star}\|)=\frac{\operatorname{L} (1-\gamma\mu)^{(n+1)/2}}{n\gamma\mu}\\ \times\left(\mathsf{E}_{\nu}^{1/2}[\|\theta_{0}-\theta^{\star}\|^{6}]+\mathsf{ E}_{\nu}^{1/2}[\|\theta_{0}-\theta^{\star}\|^{4}]+\mathsf{E}_{\nu}^{1/2}[\| \theta_{0}-\theta^{\star}\|^{2}]+\frac{\mathsf{D}_{\operatorname{last},4} \gamma\tau_{4}^{2}}{\mu}\right)\,.$ (69) $$

###### Proof.

Using the recursion (31), we obtain that

$$ $\displaystyle\operatorname{H}^{\star}(\bar{\theta}_{n}^{(RR)}-\theta^{\star})= \frac{2(\theta_{n+1}^{(\gamma)}-\theta^{\star})}{\gamma n}-\frac{2(\theta_{2n} ^{(\gamma)}-\theta^{\star})}{\gamma n}-\frac{\theta_{n+1}^{(2\gamma)}-\theta^{ \star}}{2\gamma n}+\frac{\theta_{2n}^{(2\gamma)}-\theta^{\star}}{2\gamma n}$ $\displaystyle-\frac{1}{n}\sum_{k=n+1}^{2n}[2{\varepsilon_{k+1}}(\theta_{k}^{( \gamma)})-{\varepsilon_{k+1}}(\theta_{k}^{(2\gamma)})]-\frac{1}{n}\sum_{k=n+1} ^{2n}[2\eta(\theta_{k}^{(\gamma)})-\eta(\theta_{k}^{(2\gamma)})]\,.$ (70) $$

Therefore, applying Minkowski’s inequality to the decomposition (C), we obtain for any initial distribution $\nu$ that

$$ $\displaystyle\mathsf{E}_{\nu}^{1/2}[\|\operatorname{H}^{\star}(\bar{\theta}_{n }^{(RR)}-\theta^{\star})\|^{2}]$ $\displaystyle\lesssim\underbrace{\frac{1}{n}\mathsf{E}_{\nu}^{1/2}[\|\sum_{k=n +1}^{2n}{\varepsilon_{k+1}}(\theta^{\star})\|^{2}]}_{T_{1}}+\underbrace{\frac{ 1}{\gamma n}\mathsf{E}_{\nu}^{1/2}[\|\theta_{n+1}^{(\gamma)}-\theta^{\star}\|^ {2}]+\frac{1}{\gamma n}\mathsf{E}_{\nu}^{1/2}[\|\theta_{2n}^{(\gamma)}-\theta^ {\star}\|^{2}]}_{T_{2}}$ $\displaystyle+\underbrace{\frac{1}{\gamma n}\mathsf{E}_{\nu}^{1/2}[\|\theta_{n +1}^{(2\gamma)}-\theta^{\star}\|^{2}]+\frac{1}{\gamma n}\mathsf{E}_{\nu}^{1/2} [\|\theta_{2n}^{(2\gamma)}-\theta^{\star}\|^{2}]}_{T_{3}}$ $\displaystyle+\underbrace{\frac{1}{n}\mathsf{E}_{\nu}^{1/2}[\|\sum_{k=n+1}^{2n }{\varepsilon_{k+1}}(\theta_{k}^{(\gamma)})-{\varepsilon_{k+1}}(\theta^{\star} )\|^{2}]}_{T_{4}}$ $\displaystyle+\underbrace{\frac{1}{n}\mathsf{E}_{\nu}^{1/2}[\|\sum_{k=n+1}^{2n }{\varepsilon_{k+1}}(\theta_{k}^{(2\gamma)})-{\varepsilon_{k+1}}(\theta^{\star })\|^{2}]}_{T_{5}}+\underbrace{\|2\pi_{\gamma}(\psi)-\pi_{2\gamma}(\psi)\|}_{T _{6}}$ $\displaystyle+\underbrace{\frac{1}{n}\mathsf{E}_{\nu}^{1/2}[\|\sum_{k=n+1}^{2n }\eta(\theta_{k}^{(\gamma)})-\pi_{\gamma}(\psi)\|^{2}]+\frac{1}{n}\mathsf{E}_{ \nu}^{1/2}[\|\sum_{k=n+1}^{2n}\eta(\theta_{k}^{(2\gamma)})-\pi_{2\gamma}(\psi) \|^{2}]}_{T_{7}}\,.$ $$

Now we upper bound the terms in the right-hand side of the above bound separately. First, we note that

$$ $\displaystyle T_{1}=\frac{\sqrt{\operatorname{Tr}{\Sigma_{\varepsilon}^{\star} }}}{\sqrt{n}}\,.$ $$

Using C 1($2$), we get

$$ $\displaystyle T_{2}+T_{3}\lesssim\frac{(1-\gamma\mu)^{{n+1}/2}}{\gamma n} \mathsf{E}_{\nu}^{1/2}[\|\theta_{0}-\theta^{\star}\|^{2}]+\frac{\mathsf{D}_{ \operatorname{last},2}^{1/2}\tau_{2}}{\mu^{1/2}\gamma^{1/2}n}\,.$ $$

Applying Appendix B, we get

$$ $\displaystyle T_{4}+T_{5}\lesssim\frac{\operatorname{L}\mathsf{D}_{ \operatorname{last},2}^{1/2}\gamma^{1/2}\tau_{2}}{\mu^{1/2}n^{1/2}}+\frac{ \operatorname{L}(1-\gamma\mu)^{(n+1)/2}}{\mu\gamma n}\mathsf{E}_{\nu}^{1/2} \bigl{[}\|\theta_{0}-\theta^{\star}\|^{2}\bigr{]}\,.$ $$

Now we proceed with the term $T_{6}$. Applying the recurrence (11), we obtain that

$$ $\theta_{1}^{(\gamma)}-\theta^{\star}=(\mathrm{I}-\gamma\operatorname{H}^{\star })(\theta_{0}^{(\gamma)}-\theta^{\star})-\gamma{\varepsilon_{1}}(\theta_{0}^{( \gamma)})-\gamma\eta(\theta_{0}^{(\gamma)})\,.$ (71) $$

Thus, taking expectation w.r.t. $\pi_{\gamma}$ in both sides above, we get

$$ $\operatorname{H}^{\star}(\bar{\theta}_{\gamma}-\theta^{\star})=\mathsf{E}_{\pi _{\gamma}}\bigl{[}\eta(\theta_{0}^{(\gamma)})\bigr{]}=\pi_{\gamma}(\psi)+\pi_{ \gamma}(G)\,,$ $$

where $G(\theta)$ is defined in (38) and writes as

$$ $G(\theta)=\frac{1}{6}\left(\int_{0}^{1}\nabla^{4}f(t\theta^{\star}+(1-t)\theta )\,dt\right)(\theta-\theta^{\star})^{\otimes 3}\,.$ $$

Hence, applying A 2 together with Section 3.1, we obtain that

$$ $\displaystyle T_{6}=\|2\pi_{\gamma}(\psi)-\pi_{2\gamma}(\psi)\|\lesssim\mathsf {C}_{1}\gamma^{3/2}\,.$ (72) $$

Finally, using Appendix C, Appendix C, and Appendix C, we obtain that

$$ $\displaystyle T_{7}$ $\displaystyle\lesssim\frac{\mathsf{D}_{\operatorname{last},4}\operatorname{L} \gamma\tau_{4}^{2}}{\mu n^{1/2}}+\frac{\mathsf{D}_{\operatorname{last},4} \operatorname{L}\gamma^{1/2}\tau_{4}^{2}}{\mu^{3/2}n^{1/2}}+\frac{ \operatorname{L}\mathsf{D}_{\operatorname{last},6}^{3/2}\gamma^{3/2}\tau_{6}^{ 3}}{\mu^{3/2}}$ $\displaystyle\qquad+\frac{\operatorname{L}(1-\gamma\mu)^{(n+1)/2}}{n\gamma\mu} \left(\mathsf{E}_{\nu}^{1/2}[\|\theta_{0}-\theta^{\star}\|^{6}]+\mathsf{E}^{1/ 2}_{\nu}[\|\theta_{0}-\theta^{\star}\|^{4}]+\frac{\mathsf{D}_{\operatorname{ last},4}\gamma\tau_{4}^{2}}{\mu}\right)\,.$ $$

Combining the bounds above completes the proof.
∎

Below we provide some auxiliary technical lemmas.

###### Lemma \thelemma .

Assume A 1, A 2, A 3($4$), and C 1($4$). Then for any $\gamma\in(0;1/(\operatorname{L}\mathsf{C}_{\operatorname{step},4})]$ and any $n\in\mathbb{N}$ it holds

$$ $n^{-1}\mathsf{E}^{1/2}_{\pi_{\gamma}}\bigl{[}\|\sum_{k=n+1}^{2n}\{\psi(\theta_ {k})-\pi_{\gamma}(\psi)\}\|^{2}\bigr{]}\lesssim\frac{\mathsf{D}_{\operatorname {last},4}\operatorname{L}_{3}\gamma\tau_{4}^{2}}{\mu n^{1/2}}+\frac{\mathsf{D} _{\operatorname{last},4}\operatorname{L}_{3}\gamma^{1/2}\tau_{4}^{2}}{\mu^{3/2 }n^{1/2}}\,.$ (73) $$

###### Proof.

Using the fact that $\pi_{\gamma}$ is a stationary distribution, we obtain that

$$ $\mathsf{E}_{\pi_{\gamma}}\bigl{[}\|\sum_{k=n+1}^{2n}\{\psi(\theta_{k})-\pi_{ \gamma}(\psi)\}\|^{2}\bigr{]}=n\mathsf{E}_{\pi_{\gamma}}[\|\psi(\theta_{0})- \pi_{\gamma}(\psi)\|^{2}]\\ +\sum_{k=1}^{n-1}(n-k)\mathsf{E}_{\pi_{\gamma}}[(\psi(\theta_{0})-\pi_{\gamma} (\psi))^{T}(\psi(\theta_{k})-\pi_{\gamma}(\psi))]$ $$

Using the Markov property, Cauchy–Schwartz inequality, Section 3.1, and Appendix C, we obtain

$$ $\displaystyle\mathsf{E}_{\pi_{\gamma}}[(\psi(\theta_{0})-\pi_{\gamma}(\psi))^{ T}(\psi(\theta_{k})-\pi_{\gamma}(\psi))]$ (74) $\displaystyle\qquad=\mathsf{E}_{\pi_{\gamma}}[(\psi(\theta_{0})-\pi_{\gamma}( \psi))^{T}({\rm Q}_{\gamma}^{k}\psi(\theta_{0})-\pi_{\gamma}(\psi))]$ (75) $\displaystyle\qquad\overset{(a)}{\lesssim}(1/2)^{k/m(\gamma)}\operatorname{L}_ {3}\mathsf{E}_{\pi_{\gamma}}\bigl{[}\|\psi(\theta_{0})-\pi_{\gamma}(\psi)\| \int c(\theta_{0},\vartheta)d\pi_{\gamma}(\vartheta)\bigr{]}\,,$ (76) $$

where in (a) we additionally used the fact that

$$ $\mathbf{W}_{{c}}(\delta_{\theta_{0}},\pi_{\gamma})=\int c(\theta_{0},\vartheta )d\pi_{\gamma}(\vartheta)\,.$ $$

Using C 1($4$), we get

$$ $\mathsf{E}_{\pi_{\gamma}}[\|\psi(\theta_{0})-\pi_{\gamma}\|^{2}]\leq\mathsf{E} _{\pi_{\gamma}}[\|\psi(\theta_{0})\|^{2}]\leq\operatorname{L}_{3}^{2}\mathsf{E }_{\pi_{\gamma}}[\|\theta_{0}-\theta^{\star}\|^{4}]\leq\frac{\operatorname{L}_ {3}^{2}\mathsf{D}_{\operatorname{last},4}\gamma^{2}\tau_{4}^{4}}{\mu^{2}}\,,$ (77) $$

and, using C 1($2$) and C 1($4$), we get

$$ $\displaystyle\int\int c^{2}(\theta_{0},\vartheta)d\pi_{\gamma}(\vartheta)d\pi_ {\gamma}(\theta_{0})$ (78) $\displaystyle\leq\int\int\|\theta_{0}-\vartheta\|^{2}\left(\|\theta_{0}-\theta ^{\star}\|+\|\vartheta-\theta^{\star}\|+\frac{2^{3/2}\gamma^{1/2}\tau_{2}}{\mu ^{1/2}}\right)^{2}d\pi_{\gamma}(\vartheta)d\pi_{\gamma}(\theta_{0})$ (79) $\displaystyle\lesssim\int\int(\|\theta_{0}-\theta^{\star}\|^{4}+\|\vartheta- \theta^{\star}\|^{4})+\frac{\gamma\tau_{2}^{2}}{\mu}(\|\theta_{0}-\theta^{ \star}\|^{2}+\|\vartheta-\theta^{\star}\|^{2})d\pi_{\gamma}(\vartheta)d\pi_{ \gamma}(\theta_{0})$ (80) $\displaystyle\lesssim\frac{\mathsf{D}_{\operatorname{last},4}\gamma^{2}\tau_{4 }^{2}}{\mu^{2}}+\frac{\mathsf{D}_{\operatorname{last},2}\gamma^{2}\tau_{2}^{4} }{\mu^{2}}\lesssim\frac{\mathsf{D}_{\operatorname{last},4}\gamma^{2}\tau_{4}^{ 4}}{\mu^{2}}\,.$ (81) $$

Using (77), (78), and Cauchy–Schwartz inequality for (74), we obtain

$$ $\displaystyle\mathsf{E}_{\pi_{\gamma}}[(\psi(\theta_{0})-\pi_{\gamma}(\psi))^{ T}(\psi(\theta_{k})-\pi_{\gamma}(\psi))]\lesssim(1/2)^{k/m(\gamma)}\frac{ \operatorname{L}_{3}\mathsf{D}_{\operatorname{last},4}\gamma^{2}\tau_{4}^{4}}{ \mu^{2}}\,.$ $$

Combining the inequalities above and using that $m(\gamma)=\lceil 2\frac{\log 4}{\gamma\mu}\rceil\leq\frac{2\log 4+1}{\gamma\mu}$, we get

$$ $\displaystyle n^{-1}\mathsf{E}^{1/2}_{\pi_{\gamma}}\bigl{[}\|\sum_{k=n+1}^{2n} \{\psi(\theta_{k})-\pi_{\gamma}(\psi)\}\|^{2}\bigr{]}$ $\displaystyle\leq\left(\frac{\mathsf{D}_{\operatorname{last},4}\operatorname{L }_{3}^{2}\gamma^{2}\tau_{4}^{4}}{\mu^{2}n}+\frac{\mathsf{D}_{\operatorname{ last},4}m(\gamma)\operatorname{L}_{3}^{2}\gamma^{2}\tau_{4}^{4}}{\mu^{2}n} \right)^{1/2}$ $\displaystyle\lesssim\frac{\mathsf{D}_{\operatorname{last},4}\operatorname{L}_ {3}\gamma\tau_{4}^{2}}{\mu n^{1/2}}+\frac{\mathsf{D}_{\operatorname{last},4} \operatorname{L}_{3}\gamma^{1/2}\tau_{4}^{2}}{\mu^{3/2}n^{1/2}}\,.$ $$

∎

###### Lemma \thelemma .

Assume A 1, A 2, A 3($4$). Then for any $\gamma\in(0;\frac{2}{11L}]$, and any $k\in\mathbb{N}$ it holds that

$$ $\mathsf{E}[\|\theta_{k+1}-\tilde{\theta}_{k+1}\|^{4}|{\mathcal{F}}_{k}]\leq(1- \gamma\mu)^{2}\|\theta_{k}-\tilde{\theta}_{k}\|^{4}.$ (82) $$

###### Proof.

Recall that the sequences $\{\theta_{k}\}_{k\in\mathbb{N}}$ and $\{\tilde{\theta}_{k}\}_{k\in\mathbb{N}}$ are defined by the recurrences

$$ $\displaystyle\theta_{k+1}$ $\displaystyle=\theta_{k}-\gamma\nabla F(\theta_{k},\xi_{k+1})\,,\quad\theta_{0 }=\theta\in\mathbb{R}^{d}\,,$ (83) $\displaystyle\tilde{\theta}_{k+1}$ $\displaystyle=\tilde{\theta}_{k}-\gamma\nabla F(\tilde{\theta}_{k},\xi_{k+1}) \,,\quad\tilde{\theta}_{0}=\tilde{\theta}\in\mathbb{R}^{d}\,.$ (84) $$

Expanding the brackets, we obtain that

$$ $\displaystyle\|\theta_{k+1}-\tilde{\theta}_{k+1}\|^{4}=\|\theta_{k}-\tilde{ \theta}_{k}\|^{4}+\gamma^{4}\|\nabla F(\theta_{k},\xi_{k+1})-\nabla F(\tilde{ \theta}_{k},\xi_{k+1})\|^{4}$ $\displaystyle+4\gamma^{2}\langle\nabla F(\theta_{k},\xi_{k+1})-\nabla F(\tilde {\theta}_{k},\xi_{k+1}),\theta_{k}-\tilde{\theta}_{k}\rangle^{2}$ $\displaystyle+2\gamma^{2}\|\nabla F(\theta_{k},\xi_{k+1})-\nabla F(\tilde{ \theta}_{k},\xi_{k+1})\|^{2}\|\theta_{k}-\tilde{\theta}_{k}\|^{2}$ $\displaystyle-4\gamma\langle\nabla F(\theta_{k},\xi_{k+1})-\nabla F(\tilde{ \theta}_{k},\xi_{k+1}),\theta_{k}-\tilde{\theta}_{k}\rangle\|\theta_{k}-\tilde {\theta}_{k}\|^{2}$ $\displaystyle-4\gamma^{3}\langle\nabla F(\theta_{k},\xi_{k+1})-\nabla F(\tilde {\theta}_{k},\xi_{k+1}),\theta_{k}-\tilde{\theta}_{k}\rangle\|\nabla F(\theta_ {k},\xi_{k+1})-\nabla F(\tilde{\theta}_{k},\xi_{k+1})\|^{2}$ $$

Using A 3($4$) and Cauchy–Schwartz inequality, we get

$$ $\displaystyle\mathsf{E}[\|\nabla F(\theta_{k},\xi_{k+1})-\nabla F(\tilde{ \theta}_{k},\xi_{k+1})\|^{4}|{\mathcal{F}}_{k}]$ $\displaystyle\leq\operatorname{L}^{3}\langle\nabla f(\theta_{k})-\nabla f( \tilde{\theta}_{k}),\theta_{k}-\tilde{\theta}_{k}\rangle\|\theta_{k}-\tilde{ \theta}_{k}\|^{2},$ $\displaystyle\mathsf{E}[\langle\nabla F(\theta_{k},\xi_{k+1})-\nabla F(\tilde{ \theta}_{k},\xi_{k+1}),\theta_{k}-\tilde{\theta}_{k}\rangle^{2}|{\mathcal{F}}_ {k}]$ $\displaystyle\leq\operatorname{L}\langle\nabla f(\theta_{k})-\nabla f(\tilde{ \theta}_{k}),\theta_{k}-\tilde{\theta}_{k}\rangle\|\theta_{k}-\tilde{\theta}_{ k}\|^{2},$ $\displaystyle\mathsf{E}[\|\nabla F(\theta_{k},\xi_{k+1})-\nabla F(\tilde{ \theta}_{k},\xi_{k+1})\|^{2}\|\theta_{k}-\theta^{\prime}_{k}\|^{2}|{\mathcal{F }}_{k}]$ $\displaystyle\leq\operatorname{L}\langle\nabla f(\theta_{k})-\nabla f(\tilde{ \theta}_{k}),\theta_{k}-\tilde{\theta}_{k}\rangle\|\theta_{k}-\tilde{\theta}_{ k}\|^{2}$ $\displaystyle\mathsf{E}[\langle\nabla F(\theta_{k},\xi_{k+1})-\nabla F(\tilde{ \theta}_{k},\xi_{k+1}),\theta_{k}-\tilde{\theta}_{k}\rangle\|\theta_{k}-\tilde {\theta}_{k}\|^{2}|{\mathcal{F}}_{k}]$ $\displaystyle=\langle\nabla f(\theta_{k})-\nabla f(\tilde{\theta}_{k}),\theta_ {k}-\tilde{\theta}_{k}\rangle\|\theta_{k}-\tilde{\theta}_{k}\|^{2}\,.$ $$

Similarly,

$$ $\displaystyle\mathsf{E}[\langle\nabla F(\theta_{k},\xi_{k+1})-\nabla F(\tilde{ \theta}_{k},\xi_{k+1}),\theta_{k}-\tilde{\theta}_{k}\rangle\|\nabla F(\theta_{ k},\xi_{k+1})-\nabla F(\tilde{\theta}_{k},\xi_{k+1})\|^{2}|{\mathcal{F}}_{k}]$ $\displaystyle\qquad\qquad\qquad\leq\operatorname{L}^{2}\langle\nabla f(\theta_ {k})-\nabla f(\tilde{\theta}_{k}),\theta_{k}-\tilde{\theta}_{k}\rangle\|\theta _{k}-\tilde{\theta}_{k}\|^{2}$ $$

Combining all inequalities above, we obtain

$$ $\mathsf{E}[\|\theta_{k+1}-\theta^{\prime}_{k+1}\|^{4}|{\mathcal{F}}_{k}]\leq\| \theta_{k}-\tilde{\theta}_{k}\|^{4}\\ -(4\gamma-\gamma^{4}\operatorname{L}^{3}-4\gamma^{2}\operatorname{L}-2\gamma^{ 2}\operatorname{L}-4\gamma^{3}\operatorname{L}^{2})\langle\nabla f(\theta_{k}) -\nabla f(\tilde{\theta}_{k}),\theta_{k}-\tilde{\theta}_{k}\rangle\|\theta_{k} -\tilde{\theta}_{k}\|^{2}$ $$

Using A 1 and since $1-\gamma^{3}\operatorname{L}^{3}/4-3\gamma\operatorname{L}/2-\gamma^{2}
\operatorname{L}^{2}\geq 1-11\gamma\operatorname{L}/4$, we get

$$ $\displaystyle\mathsf{E}[\|\theta_{k+1}-\tilde{\theta}_{k+1}\|^{4}|{\mathcal{F} }_{k}]$ $\displaystyle\leq(1-4\gamma\mu(1-11\gamma\operatorname{L}/4))\|\theta_{k}- \tilde{\theta}_{k}\|^{4}$ $\displaystyle\leq(1-2\gamma\mu(1-11\gamma\operatorname{L}/4))^{2}\|\theta_{k}- \tilde{\theta}_{k}\|^{4}\,.$ $$

Since $1-11\gamma\operatorname{L}/4\geq 1/2$ for $\gamma\leq 2/(11\operatorname{L})$, we complete the proof.
∎

###### Lemma \thelemma .

Assume A 1, A 2, A 3($4$), and C 1($4$). Then for any $\gamma\in(0;1/(\operatorname{L}\mathsf{C}_{\operatorname{step},4})]$, any $n\in\mathbb{N}$ and initial distribution $\nu$ it holds

$$ $\displaystyle n^{-1}\mathsf{E}^{1/2}_{\nu}[\|\sum_{k=n+1}^{2n}\{\psi(\theta_{k })-\pi_{\gamma}(\psi)\}\|^{2}]$ $\displaystyle\lesssim n^{-1}\mathsf{E}^{1/2}_{\pi_{\gamma}}[\|\sum_{k=n+1}^{2n }\{\psi(\theta_{k})-\pi_{\gamma}(\psi)\}\|^{2}]$ $\displaystyle\qquad+\frac{\operatorname{L}_{3}(1-\gamma\mu)^{(n+1)/2}}{n\gamma \mu}\left(\mathsf{E}^{1/2}_{\nu}[\|\theta_{0}-\theta^{\star}\|^{4}]+\frac{ \mathsf{D}_{\operatorname{last},4}\gamma\tau_{4}^{2}}{\mu}\right)\,.$ $$

###### Proof.

Using the synchronous coupling construction defined in (45) and the corresponding coupling kernel ${\rm K}_{\gamma}$, we obtain that

$$ $\displaystyle\mathsf{E}_{\nu}^{1/2}[\|\sum_{k=n+1}^{2n}\{\psi(\theta_{k})-\pi_ {\gamma}(\psi)\}\|^{2}]=(\mathsf{E}^{{\rm K}_{\gamma}}_{\nu,\pi_{\gamma}}[\| \sum_{k=n+1}^{2n}\{\psi(\theta_{k})-\pi_{\gamma}(\psi)\}\|^{2}])^{1/2}$ (85) $\displaystyle\quad\quad\leq\mathsf{E}^{1/2}_{\pi_{\gamma}}[\|\sum_{k=n+1}^{2n} \{\psi(\tilde{\theta}_{k})-\pi_{\gamma}(\psi)\}\|^{2}]+(\mathsf{E}^{{\rm K}_{ \gamma}}_{\nu,\pi_{\gamma}}[\|\sum_{k=n+1}^{2n}\{\psi(\theta_{k})-\psi(\tilde{ \theta}_{k})\}\|^{2}])^{1/2}$ $$

Applying Minkowski’s inequality to the last term and using Appendix C, we get

$$ $\displaystyle(\mathsf{E}^{{\rm K}_{\gamma}}_{\nu,\pi_{\gamma}}[\|\sum_{k=n+1}^ {2n}\{\psi(\theta_{k})-\psi(\tilde{\theta}_{k})\}\|^{2}])^{1/2}$ $\displaystyle\leq\sum_{k=n+1}^{2n}(\mathsf{E}^{K}_{\nu,\pi_{\gamma}}[\|\{\psi( \theta_{k})-\psi(\tilde{\theta}_{k})\}\|^{2}])^{1/2}$ $\displaystyle\leq\frac{\operatorname{L}_{3}}{2}\sum_{k=n+1}^{2n}(\mathsf{E}^{{ \rm K}_{\gamma}}_{\nu,\pi_{\gamma}}[c^{2}(\theta_{k},\tilde{\theta}_{k})])^{1/ 2}\,.$ $$

Using Hölder’s and Minkowski’s inequality and applying Appendix C , (63) and (64), we obtain

$$ $\displaystyle(\mathsf{E}^{{\rm K}_{\gamma}}_{\nu,\pi_{\gamma}}[c^{2}(\theta_{k },\tilde{\theta}_{k})])^{1/2}$ $\displaystyle\qquad\leq(\mathsf{E}^{{\rm K}_{\gamma}}_{\nu,\pi_{\gamma}}[\| \theta_{k}-\tilde{\theta}_{k}\|^{4}])^{1/4}\bigl{(}\mathsf{E}_{\pi_{\gamma}}^{ 1/4}[\|\tilde{\theta}_{k}-\theta^{\star}\|^{4}]+\mathsf{E}_{\eta}^{1/4}[\| \theta_{k}-\theta^{\star}\|^{4}+\frac{\gamma^{1/2}\tau_{2}}{\mu^{1/2}}]\bigr{)}$ $\displaystyle\qquad\leq(1-\gamma\mu)^{k/2}(\mathsf{E}^{{\rm K}_{\gamma}}_{\nu, \pi_{\gamma}}[\|\theta_{0}-\tilde{\theta}_{0}\|^{4}])^{1/4}(\mathsf{E}_{\eta}^ {1/4}[\|\theta_{0}-\theta^{\star}\|^{4}]+\frac{\mathsf{D}_{\operatorname{last} ,4}^{1/2}\gamma^{1/2}\tau_{4}}{\mu^{1/2}}+\frac{\gamma^{1/2}\tau_{2}}{\mu^{1/2 }})$ $\displaystyle\qquad\lesssim(1-\gamma\mu)^{k/2}\left(\frac{\mathsf{D}_{ \operatorname{last},4}\gamma\tau_{4}^{2}}{\mu}+\mathsf{E}_{\nu}^{1/2}\|\theta_ {0}-\theta^{\star}\|^{4}\right)$ $$

Combining all inequalities above, we get

$$ $(\mathsf{E}^{{\rm K}_{\gamma}}_{\nu,\pi_{\gamma}}[\|\sum_{k=n+1}^{2n}\{\psi( \theta_{k})-\psi(\theta_{k}^{\prime})\}\|^{2}])^{1/2}\lesssim\frac{ \operatorname{L}_{3}(1-\gamma\mu)^{(n+1)/2}}{\gamma\mu}\left(\mathsf{E}^{1/2}_ {\nu}[\|\theta_{0}-\theta^{\star}\|^{4}]+\frac{\mathsf{D}_{\operatorname{last} ,4}\gamma\tau_{4}^{2}}{\mu}\right)\,.$ $$

Substituting the last inequality into (LABEL:eq:_coupling_arbitrary_eta) we complete the proof.
∎

###### Lemma \thelemma .

Assume A 1, A 2, A 3($6$), and C 1($6$). Then for any $\gamma\in(0;1/(\operatorname{L}\mathsf{C}_{\operatorname{step},6})]$, $n\in\mathbb{N}$, and initial distribution $\nu$, it holds that

$$ $\begin{split}n^{-1}\mathsf{E}_{\nu}^{1/2}[\sum_{k=n+1}^{2n}\|\eta(\theta_{k})- \pi_{\gamma}(\psi)\|^{2}]&\leq n^{-1}\mathsf{E}_{\nu}^{1/2}\bigl{[}\sum_{k=n+1 }^{2n}\|\psi(\theta_{k})-\pi_{\gamma}(\psi)\|^{2}\bigr{]}\\ &+\frac{\operatorname{L}_{4}(1-\gamma\mu)^{(n+1)/2}}{n\gamma\mu}\mathsf{E}_{ \nu}^{1/2}[\|\theta_{0}-\theta^{\star}\|^{6}]+\frac{\operatorname{L}_{4} \mathsf{D}_{\operatorname{last},6}^{3/2}\gamma^{3/2}\tau_{6}^{3}}{3\mu^{3/2}} \,.\end{split}$ (86) $$

###### Proof.

Applying the 4-rd order Taylor expansion with integral remainder, we get that

$$ $\eta(\theta)=\psi(\theta)+\frac{1}{6}\left(\int_{0}^{1}\nabla^{4}f(t\theta^{ \star}+(1-t)\theta)\,dt\right)(\theta-\theta^{\star})^{\otimes 3}\,,$ (87) $$

and using A 2, we obtain

$$ $\|\left(\int_{0}^{1}\nabla^{4}f(t\theta^{\star}+(1-t)\theta)\,dt\right)(\theta -\theta^{\star})^{\otimes 3}\|\leq\operatorname{L}_{4}\|\theta-\theta^{\star} \|^{3}\,.$ (88) $$

Therefore, combining (87), A 2, and applying Minkowski’s inequality, we get

$$ $\mathsf{E}_{\nu}^{1/2}\bigl{[}\sum_{k=n+1}^{2n}\|\eta(\theta_{k})-\pi_{\gamma} (\psi)\|^{2}\bigr{]}\leq\mathsf{E}_{\nu}^{1/2}\bigl{[}\sum_{k=n+1}^{2n}\|\psi( \theta_{k})-\pi_{\gamma}(\psi)\|^{2}\bigr{]}\\ +\frac{\operatorname{L}_{4}}{6}\sum_{k=n+1}^{2n}\mathsf{E}_{\nu}^{1/2}[\| \theta_{k}-\theta^{\star}\|^{6}]$ (89) $$

Applying C 1($6$) for the last term of (89), we get

$$ $\displaystyle\mathsf{E}_{\nu}^{1/2}\bigl{[}\sum_{k=n+1}^{2n}\|\eta(\theta_{k}) -\pi_{\gamma}(\psi)\|^{2}\bigr{]}$ $\displaystyle\lesssim\mathsf{E}_{\nu}^{1/2}\bigl{[}\sum_{k=n+1}^{2n}\|\psi( \theta_{k})-\pi_{\gamma}(\psi)\|^{2}\bigr{]}+\frac{\operatorname{L}_{4}n \mathsf{D}_{\operatorname{last},6}^{3/2}\gamma^{3/2}\tau_{6}^{3}}{\mu^{3/2}}$ (90) $\displaystyle\qquad+\frac{\operatorname{L}_{4}(1-\gamma\mu)^{3(n+1)/2}}{1-(1- \gamma\mu)^{3/2}}\mathsf{E}_{\nu}^{1/2}[\|\theta_{0}-\theta^{\star}\|^{6}]\,.$ $$

In remains to notice that $(1-\gamma\mu)^{3/2}\leq(1-\gamma\mu)$, and the statement follows.
∎

We conclude this section with a technical statement on the properties of the function $\psi$ from (67).

###### Lemma \thelemma .

Let $\psi(\cdot)$ be a function defined in (67). Then for any $\theta,\theta^{\prime}\in\mathbb{R}^{d}$, it holds that

$$ $\|\psi(\theta)-\psi(\theta^{\prime})\|\leq\frac{1}{2}\operatorname{L}_{3}c( \theta,\theta^{\prime}).$ $$

###### Proof.

For simplicity, let us denote $T=\nabla^{3}f(\theta^{*})$.
Hence,

$$ $\|\psi(\theta)-\psi(\theta^{\prime})\|\leq\frac{1}{2}\|T(\theta-\theta^{\star} )^{\otimes 2}-T(\theta^{\prime}-\theta^{\star})^{\otimes 2}\|.$ (91) $$

Note that

$$ $\|T\|=\underset{x\neq 0,y\neq 0,z\neq 0}{\sup}\frac{\underset{i,j,k}{\sum}{T_{ ijk}x_{i}y_{j}z_{k}}}{\|x\|\|y\|\|z\|}\geq\underset{x\neq 0,y\neq 0}{\sup} \underset{z\neq 0}{\sup}\frac{\underset{k}{\sum}z_{k}\underset{i,j}{\sum}T_{ ijk}x_{i}y_{j}}{\|z\|\|y\|\|x\|}=\underset{x\neq 0,y\neq 0}{\sup}\frac{\|t(x,y )\|}{\|y\|\|x\|},$ (92) $$

where $t(x,y)_{k}=\underset{i,j}{\sum}T_{ijk}x_{i}y_{j}$.
Therefore, for any $x,y\in\mathbb{R}^{d}$, it holds that

$$ $\|t(x,y)\|\leq\|x\|\|y\|\|T\|$ (93) $$

We denote $v=Tx^{\otimes 2}-Ty^{\otimes 2}$. Then

$$ $\displaystyle v_{k}=\underset{i,j}{\sum}T_{ijk}(x_{i}x_{j}-y_{i}y_{j})= \underset{i,j}{\sum}T_{ijk}((x_{i}-y_{i})x_{j}+(x_{i}-y_{i})y_{j})=$ $\displaystyle\underset{i,j}{\sum}T_{ijk}(x_{i}-y_{i})x_{j}+\underset{i,j}{\sum }T_{ijk}(x_{i}-y_{i})y_{j},$ (94) $$

where the first inequality is true since $T_{ijk}=T_{jik}$ by definition of $T$.
Combining (93) and (C) and using triangle inequality, we obtain

$$ $\|v\|\leq\|T\|\|x-y\|(\|x\|+\|y\|)\leq\|T\|\|x-y\|(\|x\|+\|y\|+\frac{2\sqrt{2} \tau_{2}\sqrt{\gamma}}{\sqrt{\mu}}).$ $$

We complete the proof setting $x=\theta-\theta^{\star},y=\theta^{\prime}-\theta^{\star}$
∎

###### Theorem 5 (Version of Theorem 2 with explicit constants) .

###### Proof.

###### Lemma \thelemma .

###### Proof.

###### Lemma \thelemma .

###### Proof.

###### Lemma \thelemma .

###### Proof.

###### Lemma \thelemma .

###### Proof.

###### Lemma \thelemma .

###### Proof.

## Appendix D Proof of Theorem 3

###### Theorem 6 (Version of Theorem 3 with explicit constants) .

Let $p\geq 2$ and assume A 1, A 2, A 3($3p$), and C 1($3p$). Then for any $\gamma\in(0;1/(\operatorname{L}\mathsf{C}_{\operatorname{step},3p})]$, initial distribution $\nu$, and $n\in\mathbb{N}$, the estimator $\bar{\theta}_{n}^{(RR)}$ defined in (37) satisfies

$$ $\displaystyle\mathsf{E}_{\nu}^{1/p}[\|\operatorname{H}^{\star}(\bar{\theta}_{n }^{(RR)}-\theta^{\star})\|^{p}]$ $\displaystyle\lesssim\frac{\sqrt{\operatorname{Tr}{\Sigma_{\varepsilon}^{\star }}}p^{1/2}}{n^{1/2}}+\frac{p\tau_{p}}{n^{1-1/p}}+\frac{\mathsf{C}_{ \operatorname{RR},5}}{n\gamma^{1/2}}+\frac{\mathsf{C}_{\operatorname{RR},6} \gamma^{1/2}}{n^{1/2}}+\mathsf{C}_{\operatorname{RR},7}\gamma^{3/2}$ $\displaystyle\qquad\qquad+\frac{\mathsf{C}_{\operatorname{RR},8}}{n}+\mathcal{ R}_{5}(n,\gamma,\|\theta_{0}-\theta^{\star}\|)\,,$ $$

where we have set

$$ $\begin{split}\mathsf{C}_{\operatorname{RR},5}&=\frac{\mathsf{D}_{\operatorname {last},p}^{1/2}\tau_{p}}{\mu^{1/2}}\,,\quad\mathsf{C}_{\operatorname{RR},6}= \frac{\operatorname{L}\mathsf{D}_{\operatorname{last},p}^{1/2}p\tau_{p}}{\mu^{ 1/2}}+\frac{\operatorname{L}\mathsf{D}_{\operatorname{last},2p}p\tau^{2}_{2p}} {\mu^{3/2}}\,,\\ \mathsf{C}_{\operatorname{RR},7}&=\mathsf{C}_{1}+\frac{\operatorname{L}\mathsf {D}_{\operatorname{last},3p}^{3/2}\tau_{3p}^{3}}{\mu^{3/2}}\,,\quad\mathsf{C}_ {\operatorname{RR},8}=\frac{\operatorname{L}\mathsf{D}_{\operatorname{last},2p }\tau_{2p}}{\mu^{2}}\,,\end{split}$ (95) $$

$\mathsf{C}_{1}$ is defined in (50), and the remainder term $\mathcal{R}_{5}(n,\gamma,\|\theta_{0}-\theta^{\star}\|)$ is given by

$$ $\mathcal{R}_{5}(n,\gamma,\|\theta_{0}-\theta^{\star}\|)=\frac{(1-\gamma\mu)^{( n+1)/2}}{\gamma n}\mathsf{E}_{\nu}^{1/p}\bigl{[}\|\theta_{0}-\theta^{\star}\|^ {p}\bigr{]}+\frac{\operatorname{L}p(1-\gamma\mu)^{(n+1)/2}}{\mu^{1/2}\gamma^{1 /2}n}\mathsf{E}^{1/p}_{\nu}[\|\theta_{0}-\theta^{\star}\|^{p}]\\ +\frac{\operatorname{L}(1-\gamma\mu)^{(n+1)/2}p^{2}}{\gamma\mu^{2}}\mathsf{E}_ {\nu}^{1/p}\bigl{[}\|\theta_{0}-\theta^{\star}\|^{2p}\bigr{]}+\frac{ \operatorname{L}(1-\gamma\mu)^{(3/2)n}}{\gamma\mu}\mathsf{E}_{\nu}^{1/p}\bigl{ [}\|\theta_{0}-\theta^{\star}\|^{3p}\bigr{]}$ (96) $$

###### Proof.

Using the decomposition (C), we obtain that for any $p\geq 2$, it holds that

$$ $\displaystyle\mathsf{E}_{\nu}^{1/p}[\|\operatorname{H}^{\star}(\bar{\theta}_{n }^{(RR)}-\theta^{\star})\|^{p}]$ $\displaystyle\lesssim\underbrace{\frac{1}{n}\mathsf{E}_{\nu}^{1/p}[\|\sum_{k=n +1}^{2n}{\varepsilon_{k+1}}(\theta^{\star})\|^{p}]}_{T_{1}}+\underbrace{\frac{ 1}{\gamma n}\mathsf{E}_{\nu}^{1/p}[\|\theta_{n+1}^{(\gamma)}-\theta^{\star}\|^ {p}]+\frac{1}{\gamma n}\mathsf{E}_{\nu}^{1/p}[\|\theta_{2n}^{(\gamma)}-\theta^ {\star}\|^{p}]}_{T_{2}}$ $\displaystyle+\underbrace{\frac{1}{\gamma n}\mathsf{E}_{\nu}^{1/p}[\|\theta_{n +1}^{(2\gamma)}-\theta^{\star}\|^{p}]+\frac{1}{\gamma n}\mathsf{E}_{\nu}^{1/p} [\|\theta_{2n}^{(2\gamma)}-\theta^{\star}\|^{p}]}_{T_{3}}$ $\displaystyle+\underbrace{\frac{1}{n}\mathsf{E}_{\nu}^{1/p}[\|\sum_{k=n+1}^{2n }{\varepsilon_{k+1}}(\theta_{k}^{(\gamma)})-{\varepsilon_{k+1}}(\theta^{\star} )\|^{p}]}_{T_{4}}$ $\displaystyle+\underbrace{\frac{1}{n}\mathsf{E}_{\nu}^{1/p}[\|\sum_{k=n+1}^{2n }{\varepsilon_{k+1}}(\theta_{k}^{(2\gamma)})-{\varepsilon_{k+1}}(\theta^{\star })\|^{p}]}_{T_{5}}+\underbrace{\|2\pi_{\gamma}(\psi)-\pi_{2\gamma}(\psi)\|}_{T _{6}}$ $\displaystyle+\underbrace{\frac{1}{n}\mathsf{E}_{\nu}^{1/p}[\|\sum_{k=n+1}^{2n }\psi(\theta_{k}^{(\gamma)})-\pi_{\gamma}(\psi)\|^{p}]+\frac{1}{n}\mathsf{E}_{ \nu}^{1/p}[\|\sum_{k=n+1}^{2n}\psi(\theta_{k}^{(2\gamma)})-\pi_{2\gamma}(\psi) \|^{p}]}_{T_{7}}$ $\displaystyle+\underbrace{\frac{1}{n}\sum_{k=n+1}^{2n}\mathsf{E}^{1/p}_{\nu} \bigl{[}\|G(\theta_{k}^{(\gamma)})\|^{p}\bigr{]}+\frac{1}{n}\sum_{k=n+1}^{2n} \mathsf{E}^{1/p}_{\nu}\bigl{[}\|G(\theta_{k}^{(2\gamma)})\|^{p}\bigr{]}}_{T_{8 }}\,.$ $$

Now we upper bounds the terms above separately. Applying first the Pinelis version of Rosenthal inequality (Pinelis, ) together with A 3($p$), we obtain that

$$ $\displaystyle T_{1}\lesssim\frac{\sqrt{\operatorname{Tr}{\Sigma_{\varepsilon}^ {\star}}}p^{1/2}}{n^{1/2}}+\frac{p\tau_{p}}{n^{1-1/p}}\,.$ $$

Applying C 1($p$) (which is implied by C 1($3p$)), we obtain that

$$ $\displaystyle T_{2}+T_{3}\lesssim\frac{\mathsf{D}_{\operatorname{last},p}^{1/2 }\tau_{p}}{\mu^{1/2}n\gamma^{1/2}}+\frac{(1-\gamma\mu)^{(n+1)/2}}{\gamma n} \mathsf{E}_{\nu}^{1/p}\bigl{[}\|\theta_{0}-\theta^{\star}\|^{p}\bigr{]}\,.$ $$

Applying Appendix B (see the bound (66)), we get that

$$ $\displaystyle T_{4}+T_{5}\lesssim\frac{\operatorname{L}\mathsf{D}_{ \operatorname{last},p}^{1/2}\gamma^{1/2}p\tau_{p}}{\mu^{1/2}n^{1/2}}+\frac{ \operatorname{L}p(1-\gamma\mu)^{(n+1)/2}}{\mu^{1/2}\gamma^{1/2}n}\mathsf{E}^{1 /p}_{\nu}[\|\theta_{0}-\theta^{\star}\|^{p}]\,.$ $$

Using the bounds (71) and (72), we obtain

$$ $\displaystyle T_{6}\lesssim\mathsf{C}_{1}\gamma^{3/2}\,.$ $$

Applying Section 4, we get

$$ $\displaystyle\frac{1}{n}\mathsf{E}_{\nu}^{1/p}[\|\sum_{k=n+1}^{2n}\psi(\theta_ {k}^{(\gamma)})-\pi_{\gamma}(\psi)\|^{p}]\lesssim\frac{\operatorname{L}\mathsf {D}_{\operatorname{last},2p}p\tau^{2}_{2p}\gamma^{1/2}}{\mu^{3/2}n^{1/2}}+ \frac{\operatorname{L}\mathsf{D}_{\operatorname{last},2p}\tau_{2p}}{\mu^{2}n}\,.$ $$

Using this bound and adopting the result of (Durmus et al., , Theorem 4), we obtain that

$$ $\displaystyle T_{7}\lesssim\frac{\operatorname{L}\mathsf{D}_{\operatorname{ last},2p}p\tau^{2}_{2p}\gamma^{1/2}}{\mu^{3/2}n^{1/2}}+\frac{\operatorname{L} \mathsf{D}_{\operatorname{last},2p}\tau_{2p}}{\mu^{2}n}+\frac{\operatorname{L} (1-\gamma\mu)^{(n+1)/2}p^{2}}{\gamma\mu^{2}}\mathsf{E}_{\nu}^{1/p}\bigl{[}\| \theta_{0}-\theta^{\star}\|^{2p}\bigr{]}\,.$ $$

Finally, applying the definition of $G(\theta)$ in (38) together with C 1($3p$), we obtain that

$$ $\displaystyle T_{8}$ $\displaystyle\lesssim\frac{\operatorname{L}\mathsf{D}_{\operatorname{last},3p} ^{3/2}\gamma^{3/2}\tau_{3p}^{3}}{\mu^{3/2}}+\frac{\operatorname{L}}{n}\sum_{k= n+1}^{2n}(1-\gamma\mu)^{(3/2)k}\mathsf{E}_{\nu}^{1/p}\bigl{[}\|\theta_{0}- \theta^{\star}\|^{3p}\bigr{]}$ $\displaystyle\lesssim\frac{\operatorname{L}\mathsf{D}_{\operatorname{last},3p} ^{3/2}\gamma^{3/2}\tau_{3p}^{3}}{\mu^{3/2}}+\frac{\operatorname{L}(1-\gamma\mu )^{(3/2)n}}{\gamma\mu}\mathsf{E}_{\nu}^{1/p}\bigl{[}\|\theta_{0}-\theta^{\star }\|^{3p}\bigr{]}\,.$ $$

To complete the proof it remains to combine the bounds for $T_{1}$ to $T_{8}$.
∎

###### Theorem 6 (Version of Theorem 3 with explicit constants) .

###### Proof.

### D.1 Proof of Section 4

In the proof below we use the notation

$$ $\bar{\psi}(\theta)=\psi(\theta)-\pi_{\gamma}(\psi)\,.$ $$

We proceed with the blocking technique. Indeed, let us set the parameter

$$ $m=m(\gamma)=\biggl{\lceil}\frac{2\log 4}{\gamma\mu}\biggr{\rceil}\,.$ (97) $$

Our choice of parameter $m(\gamma)$ is due to Section 3.1. For notation conciseness we write it simply as $m$, dropping its dependence upon $\gamma$. Using Minkowski’s inequality, we obtain that

$$ $\mathsf{E}_{\pi_{\gamma}}^{1/p}\bigl{[}\|\sum_{k=0}^{n-1}\bar{\psi}(\theta_{k} )\|^{p}\bigr{]}\leq\mathsf{E}_{\pi_{\gamma}}^{1/p}\bigl{[}\|\sum_{k=0}^{ \lfloor n/m\rfloor m-1}\bar{\psi}(\theta_{k})\|^{p}\bigr{]}+m\mathsf{E}_{\pi_{ \gamma}}^{1/p}\bigl{[}\|\bar{\psi}(\theta_{0})\|^{p}\bigr{]}\,.$ (98) $$

Now we consider the Poisson equation, associated with $\mathrm{Q}_{\gamma}^{m}$ and function $\bar{\psi}$, that is,

$$ $g_{m}(\theta)-\mathrm{Q}_{\gamma}^{m}g_{m}(\theta)=\bar{\psi}(\theta)\,.$ (99) $$

The function

$$ $g_{m}(\theta)=\sum_{k=0}^{\infty}\mathrm{Q}_{\gamma}^{km}\bar{\psi}(\theta)$ (100) $$

is well-defined under the assumptions A 1, A 2, A 3($2p$), and C 1($2p$). Moreover, $g_{m}$ is a solution of the Poisson equation (99). Define $q:=\lfloor n/m\rfloor$, then we have

$$ $\sum_{k=0}^{qm-1}\bar{\psi}(\theta_{k})=\sum_{r=0}^{m-1}B_{m,r}\,,\quad\text{ with}\quad B_{m,r}=\sum_{k=0}^{q-1}\bigl{\{}g_{m}(\theta_{km+r})-\mathrm{Q}_{ \gamma}^{m}g_{m}(\theta_{km+r})\bigr{\}}\,.$ (101) $$

Using Minkowski’s inequality, we get from (98), that

$$ $\mathsf{E}_{\pi_{\gamma}}^{1/p}\bigl{[}\|\sum_{k=0}^{n-1}\bar{\psi}(\theta_{k} )\|^{p}\bigr{]}\leq m\mathsf{E}_{\pi_{\gamma}}^{1/p}\bigl{[}\|\sum_{k=1}^{q} \left\{g_{m}(\theta_{km})-\mathrm{Q}_{\gamma}^{m}g_{m}(\theta_{(k-1)m})\right \}\|^{p}\bigr{]}+2m\mathsf{E}_{\pi_{\gamma}}^{1/p}\bigl{[}\|\psi(\theta_{0}\|^ {p}\bigr{]}$ (102) $$

Now we upper bound both terms of (102) separately. Under assumption A 2, and applying C 1($2p$), we get

$$ $\mathsf{E}_{\pi_{\gamma}}^{1/p}\bigl{[}\|\psi(\theta_{0}\|^{p}\bigr{]}\leq \frac{\operatorname{L}}{2}\mathsf{E}_{\pi_{\gamma}}^{1/p}\bigl{[}\|\theta_{0}- \theta^{\star}\|^{2p}\bigr{]}\leq\frac{\operatorname{L}\mathsf{D}_{ \operatorname{last},2p}\gamma\tau_{2p}^{2}}{2\mu}\,.$ (103) $$

To proceed with the first term, we apply Burkholder’s inequality (Osekowski, , Theorem 8.6), and obtain that

$$ $\mathsf{E}_{\pi_{\gamma}}^{1/p}\bigl{[}\|\sum_{k=1}^{q}\left\{g_{m}(\theta_{km })-\mathrm{Q}_{\gamma}^{m}g_{m}(\theta_{(k-1)m})\right\}\|^{p}\bigr{]}\\ \leq p\mathsf{E}_{\pi_{\gamma}}^{1/p}\bigr{[}\bigl{(}\sum_{k=1}^{q}\|\left\{g_ {m}(\theta_{km})-\mathrm{Q}_{\gamma}^{m}g_{m}(\theta_{(k-1)m})\right\}\|^{2} \bigr{)}^{p/2}\bigr{]}\,.$ (104) $$

Applying now Minkowski’s inequality again, we get

$$ $\displaystyle\mathsf{E}_{\pi_{\gamma}}^{2/p}\bigl{[}\bigl{(}\sum_{k=1}^{q}\| \left\{g_{m}(\theta_{km})-\mathrm{Q}_{\gamma}^{m}g_{m}(\theta_{(k-1)m})\right \}\|^{2}\bigr{)}^{p/2}\bigr{]}\leq q\mathsf{E}_{\pi_{\gamma}}^{2/p}\bigl{[}\| \left\{g_{m}(\theta_{km})-\mathrm{Q}_{\gamma}^{m}g_{m}(\theta_{(k-1)m})\right \}\|^{p}\bigr{]}$ $\displaystyle\qquad\qquad\qquad\qquad\lesssim q\left(\mathsf{E}^{2/p}_{\pi_{ \gamma}}[\|g_{m}(\theta_{0})\|^{p}]+\mathsf{E}^{2/p}_{\pi_{\gamma}}[\|\mathrm{ Q}_{\gamma}^{m}g_{m}(\theta_{0})\|^{p}]\right)$ $\displaystyle\qquad\qquad\qquad\qquad\lesssim q\mathsf{E}^{2/p}_{\pi_{\gamma}} [\|g_{m}(\theta_{0})\|^{p}]\,.$ $$

It remains to upper bound the moment $\mathsf{E}^{2/p}_{\pi_{\gamma}}[\|g_{m}(\theta_{0})\|^{p}]$. In order to do this, we first note that due to the duality theorem (Douc et al., , Theorem 20.1.2.), we get that for any $k\in\mathbb{N}$,

$$ $\|\mathrm{Q}^{mk}\psi(\theta)-\pi_{\gamma}(\psi)\|\leq\frac{1}{2}\operatorname {L}_{3}\mathbf{W}_{{c}}(\delta_{\theta}{\rm Q}_{\gamma}^{km},\pi_{\gamma})\leq 2 \operatorname{L}_{3}(1/2)^{k}\mathbf{W}_{{c}}(\delta_{\theta},\pi_{\gamma})\,,$ $$

where the last inequality is due to Section 3.1. Hence, applying the definition of $g_{m}(\theta)$ in (100), we obtain that

$$ $\displaystyle\mathsf{E}^{1/p}_{\pi_{\gamma}}[\|g_{m}(\theta_{0})\|^{p}]\leq \sum_{k=0}^{\infty}\mathsf{E}^{1/p}_{\pi_{\gamma}}\bigl{[}\|{\rm Q}_{\gamma}^{ km}\bar{\psi}(\theta)\|^{p}\bigr{]}\leq 2\operatorname{L}_{3}\sum_{k=0}^{ \infty}(1/2)^{k}\mathsf{E}^{1/p}_{\pi_{\gamma}}\bigl{[}\{\mathbf{W}_{{c}}( \delta_{\theta},\pi_{\gamma})\}^{p}\bigr{]}\,.$ $$

To control the latter term, we simply apply the definition of $\mathbf{W}_{{c}}(\delta_{\theta},\pi_{\gamma})$ and a cost function $c(\theta,\theta^{\prime})$ together with C 1($2p$), we get

$$ $\displaystyle\mathsf{E}^{1/p}_{\pi_{\gamma}}\bigl{[}\{\mathbf{W}_{{c}}(\delta_ {\theta},\pi_{\gamma})\}^{p}\bigr{]}\lesssim\left(\int_{\mathbb{R}^{d}\times \mathbb{R}^{d}}\|\theta-\theta^{\prime}\|^{p}\left(\|\theta-\theta^{\star}\|+ \|\theta^{\prime}-\theta^{\star}\|+\frac{\tau_{2}\sqrt{\gamma}}{\sqrt{\mu}} \right)^{p}\pi_{\gamma}(\mathrm{d}\theta)\pi_{\gamma}(\mathrm{d}\theta^{\prime })\right)^{1/p}$ $\displaystyle\qquad\leq\bigl{(}\int\|\theta-\theta^{\prime}\|^{2p}\pi_{\gamma} (\mathrm{d}\theta)\pi_{\gamma}(\mathrm{d}\theta^{\prime})\bigr{)}^{1/2p}\bigl{ (}\int\left(\|\theta-\theta^{\star}\|+\|\theta^{\prime}-\theta^{\star}\|+\frac {\tau_{2}\sqrt{\gamma}}{\sqrt{\mu}}\right)^{2p}\pi_{\gamma}(\mathrm{d}\theta) \pi_{\gamma}(\mathrm{d}\theta^{\prime})\bigr{)}^{1/2p}$ $\displaystyle\qquad\lesssim\frac{\mathsf{D}_{\operatorname{last},2p}\tau^{2}_{ 2p}\gamma}{\mu}\,.$ $$

Combining now the bounds above in (104), we get that

$$ $\mathsf{E}_{\pi_{\gamma}}^{1/p}\bigl{[}\|\sum_{k=1}^{q}\left\{g_{m}(\theta_{km })-\mathrm{Q}_{\gamma}^{m}g_{m}(\theta_{(k-1)m})\right\}\|^{p}\bigr{]}\lesssim \frac{\mathsf{D}_{\operatorname{last},2p}p\operatorname{L}_{3}\tau^{2}_{2p} \gamma\sqrt{q}}{\mu}\,,$ (105) $$

and, hence, substituting into (98), we get

$$ $\mathsf{E}_{\pi_{\gamma}}^{1/p}\bigl{[}\|\sum_{k=0}^{n-1}\bar{\psi}(\theta_{k} )\|^{p}\bigr{]}\lesssim\frac{\mathsf{D}_{\operatorname{last},2p}p\operatorname {L}_{3}\tau^{2}_{2p}\gamma\sqrt{q}m}{\mu}+\frac{\operatorname{L}\mathsf{D}_{ \operatorname{last},2p}\tau_{2p}^{2}\gamma m}{2\mu}\,.$ (106) $$

Now the statement follows from the definition of $m=m(\gamma)$ in (97) and $q=\lfloor n/m\rfloor\leq n/m$.

## Peer Reviews

confidence: 2
contribution: 3
initial_rating: 6
presentation: 3
questions: See **Weaknesses** part.
soundness: 3
strengths: 1. The authors provide interesting results on SGD with Richardson-Romberg extrapolation.
2. I also appreciate the authors' honest discussion of their stronger assumption compared to the previous works, e.g., the higher-order smoothness. 
3. Moreover, the writing is easy to follow.
summary: The paper investigates the convergence of stochastic gradient descent (SGD) for strongly convex and smooth problems. Concretely, using the mean squared distance as the performance metric, the author proves that the Polyak-Ruppert averaging procedure with the Richardson-Romberg extrapolation technique can achieve the optimal convergence rate for both the leading and higher-order terms.
weaknesses: Since I am not very familiar with related works of the Richardson-Romberg extrapolation, the author could feel free to point out anything wrong in my following question.

1. Could the authors provide more discussion on C1? Especially, when can (20) hold?

2. The paper studies the convergence in expectation, is it possible to establish the high-probability convergence result with the same rate for both leading and higher-order terms?

3. One point I found is that when the stochastic problem degenerates into the deterministic problem, the current rate cannot recover the well-known linear convergence result. Could the authors make some statement on this point?

4. Another thing I am interested in is that the leading term $O(n^{-1/2})$ seems not to be infected by the stepsize (even if without the Richardson-Romberg extrapolation), which seems very different from the existing analysis of SGD for the same problem. Could the authors elaborate more on this fact?

---

confidence: 3
contribution: 3
initial_rating: 6
presentation: 2
questions: Could the authors provide a geometric interpretation for the definition of $c(\theta, \theta')$, and the ensuing geometric ergodicity of the SGD iterates with respect to $W_c$? Also, are the authors aware of other contexts were such non-standard distance-like functions have been used, which could be cited in the paper?
soundness: 3
strengths: The paper interprets SGD updates with constant step sizes as an homogenous Markov chain, and uses a trick from numerical analysis to show that the difference of the average processes obtained for step sizes $\gamma$ and $2\gamma$ (defined in eq. (37)) will converge to a point that is closer to the optimum (as shown by Dieuleveut et al. 2018).
summary: The paper considers the minimization of an objective function through unbiased estimates of its gradients. It brings tools from the theory of Markov chains to the study of stochastic gradient descent, showing that, under suitable assumptions, the Richardson-Romberg extrapolation achieves the optimal non-asymptotic risk bound (proved by Li et al. 2022 for the Root-SGD algorithm). Extensions to the p-th moment bound are also provided.
weaknesses: The paper could use some more polishing in its presentation:
* The geometric ergodicity of the sequence of updates is a crucial aspect, yet is not defined in the preliminaries. It would also be valuable to add an intuitive explanations for how this fact would be used to derive the risk bounds. 
* $D_{{\rm last}, 2}$ used eq. (19) in line 222 is not defined.
* In line 480, the beginning of the sentence should be "Directions ..."
* The authors should also make explicit the fact that the number of samples, $n$, needs to be known a priori to optimize the step size $\gamma$.

---

confidence: 4
contribution: 2
initial_rating: 6
presentation: 3
questions: 1. The extensive use of $\lesssim$ is rather confusing and potentially undermining in the paper. For example, in Corollaries 4 and 7, are there multiplicative constants suppressed on the leading $n^{-1/2} \sqrt{ \operatorname{tr} \Sigma }$ terms? 

* If there are: (1) these corollaries do not match in sharpness the existing results, and provide little insight into how RR-SGD compares to PR-SGD as the latter is favorable exactly because of its Cramer-Rao-type optimality in the leading term. 
(2) is it possible to spell out and compare the explicit leading terms?

* If there are not, I wonder why $ \lesssim$ instead of $\le$ is used.

Given how the notation $\lesssim$ is defined in Line 118 I'm afraid the former is the case, right?

Either way, the authors might want to say a few words on what terms/constants are being suppressed at least for presentation purposes.


----


2. Along the same line, we know that (considering decreasing step sizes for the moment --- you might be able to work back and forth between these two regimes) with PR-SGD the asymptotic normal distribution one can prove for
$\sqrt{t} (\widehat{\theta}_{t}^{\mathsf{avg}} - \theta^*)$
has a variance that matches the Cramer-Rao lower bound, 
which happens for any step sizes $ \gamma_t \propto t^{-w}$ where $w\in[1/2, 1]$.  Whereas for SGD without averaging, the asymptotic normality takes the form of 

$t^{w - 1/2} ( \widehat{\theta}_{t} - \theta^* ) $ where $w$ is the decay rate of step sizes defined above. Even in the $ w= 1 $ case, one gets a variance that's strictly larger than Cramer-Rao (hence than that of the iterate average). When $ w<1 $ this is a strictly slower rate of convergence. 

Long story short, PR averaging accelerates the optimization via a "removal of the step size effect". Is it the case that RR has a similar effect without averaging?
soundness: 4
strengths: - The paper is well-motivated and well-written. In particular, the last-iterate bound for RR-SGD, despite some issues raised below, is novel.
summary: The authors provide convergence analysis of both the Polyak-Ruppert averaged SGD (PR-SGD) and the Richardson-Romberg extrapolated SGD (RR-SGD) by a two-term expansion of the quantity $ \mathbb E^{1/2}[ || \hat \theta - \theta^* ||^2  ]$, a quantitative refinement of asymptotic normality. They show that RR-SGD achieves the same rate as PR-SGD which is known to be optimal.
weaknesses: See Questions

---

confidence: 3
contribution: 3
initial_rating: 6
presentation: 3
questions: 1) As stated: can the authors show that they can realize this n^{3/4} correction term in an experiment, (the setup should be such that I can reproduce it and, the 3/4 fit cleanly on a log-log plot)?  

2) The optimality of the 3/4 is asserted multiple places, but it is nowhere displayed.  What construction demonstrates the non-improvability of the 3/4?
soundness: 3
strengths: This improves an existing bound in various ways, both in the sense of
(1) achieving moment bounds and
(2) improving the rate of the correction term from $n^{1/2+\delta}$ to $n^{3/4}$, which is claimed as optimal.

It does it using the method of the `Richardson-Romberg' extrapolation.  

For more common Polyak-Ruppert averaging, it also achieves the same rate provided the minimization is strongly convex.
summary: The paper contains an analysis of one-pass SGD on stochastic objective functions:

$$\min\_{\theta \in \mathbb{R}^d} f(\theta), \quad \nabla f(\theta)=\mathrm{E}\_{\xi \sim \mathbb{P}_{\xi}}[\nabla F(\theta, \xi)],$$

This is run with constant step-size $\gamma \asymp 1/\sqrt{n}$ with $n$ the number of samples,  and it is shown that $p$-th moments satisfy a bound

$$\mathrm{E}^{1 / p}\left[\left\|\bar{\theta}\_n^{(R R)}-\theta^{\star}\right\|^p\right] \leq \frac{C p^{1 / 2} \sqrt{\operatorname{Tr} \Sigma_{\infty}}}{n^{1 / 2}}+\frac{C(f, d, p)}{n^{3 / 4}}+\mathcal{R}\left(\left\|\theta_0-\theta^{\star}\right\|, n, p\right),$$

Here 
* $\bar{\theta}\_n^{(R R)}$ is the ``Richardson-Romburg'' estimator for reducing lower order terms (specifically achieving the $n^{3/4}$ correction.
* $\mathcal{R}$ is an error term that controls the dependence on initial condition.

This improves the correction term in existing work from $n^{1/2+\delta}$ to $n^{3/4}$.
weaknesses: 1) This paper performs a detailed mathematical analysis of a correction term to the behavior of SGD on fixed problems.  It is hard to imagine important practical conclusions of this -- we are talking about improving subleading correction terms in an _asymptotic_ setting when the number of samples overwhelms all the important problem-dependent constants in a stochastic optimization context.  And yes, I use the term _asymptotic_ intentionally: the $n$ in this problem has to overwhelm constants which are _quite clearly_ not effective, and only once that happens will the subleading terms will have improved dependence.   The authors have _not_ attempted to demonstrate that there are settings where this theorem is experimentally realizable, nor honestly, do I think it would be reasonable to expect them to do so (if they can show their 3/4-power in an experiment, I'd retract this and happily ``eat crow'').

2) The 3/4 power seems to be the essential point of this story.  It is asserted at multiple places that the 3/4 is optimal, but there is no clear construction.  There is a citation to another paper (Li et al), which considers a different algorithm.  If the example from that other paper works, it should be either displayed or precisely cited here (it was not possible to find it in a quick perusal).

Meta-review: This paper considers constant step size SGD on strongly convex and smooth objectives and provides a Markov chain based analysis. Authors derive a bias-variance decomposition of the MSE with a dominant term and a higher term of respective orders 1/n^.5 and 1/n^.75. 


This paper was reviewed by four expert reviewers and received the following Scores/Confidence: 6/3, 6/3, 6/4, 6/2. I think the paper is studying an interesting topic and the results are relevant to ICLR community. The following concerns were brought up by the reviewers:

- There are some important and relevant references on the analysis of "constant step size SGD" missing. Authors should include references to existing works on this algorithm.

- The optimality of power 3/4 is not explained clearly in the paper. Examples would help.
  
- Several typos and ambiguous statements/notations were pointed by the reviewers. These should be carefully addressed.

- Certain technical terminology are not clearly defined in the paper. While I admit not all terms can be included in the paper, I suggest authors to make their paper more friendly to non-experts on MCMC.


Authors should carefully go over reviewers' suggestions and address any remaining concerns in their final revision.  Based on the reviewers' suggestion, as well as my own assessment of the paper, I recommend including this paper to the ICLR 2025 program.

Decision: Accept (Poster)