text:
\section{Supplementary}

\subsection{Consequences of our Assumptions}  \label{sec:examples}
We repeatedly use the following result on input distributions of assumption~\ref{assume:symmetric_inputs}.
\begin{lemma} \label{lemma:symmetric_inputs}
Suppose matrix ${\bf M}$ can be diagonalised as
\begin{align}
{\bf M} = {\bf U} {\textbf{diag}}(\boldsymbol{ \lambda}) {\bf U}^\top
\end{align}
Under Assumption~\ref{assume:symmetric_inputs} holds, the following holds:
\begin{align}
 {\mathbf E} \left[ {\bf x} {\bf x}^\top {\bf M} {\bf x} {\bf x}^\top \right] = {\bf U}^\top {\textbf{diag}}(\boldsymbol{ \lambda}') {\bf U}, \quad \boldsymbol{ \lambda}'_i
 = \left({\textbf{diag}}\left(\k- (\boldsymbol{\sigma})^2\right)  + \boldsymbol{\sigma} \boldsymbol{\sigma}^\top\right) \boldsymbol{ \lambda}
\end{align}
\end{lemma}
\begin{proof}
Using the structure of ${\bf M}$, we first simplify the expression of ${\bf x}^\top {\bf M} {\bf x}$ :
\begin{align}
 {\bf x}^\top {\bf M} {\bf x} & = \v^\top {\bf U}^\top {\bf U} {\textbf{diag}}(\boldsymbol{ \lambda}) {\bf U}^\top {\bf U} \v \\
 & = \v^\top {\textbf{diag}}(\boldsymbol{ \lambda}) \v = \sum_{i=1}^d \v_i^2 \boldsymbol{ \lambda}_i
\end{align}
Replacing this into the desired expectation yields
\begin{align} \label{eq:symmetric_xMx_expansion}
 {\mathbf E} \left[ {\bf x} {\bf x}^\top {\bf M} {\bf x} {\bf x}^\top \right]  = {\bf U} {\mathbf E} \left[ \left( \sum_{i=1}^d \v_i^2 \boldsymbol{ \lambda}_i
 \right) \v \v^\top \right] {\bf U}.
\end{align}
Since distribution of $\v_i$ is symmetric, we conclude
\begin{align} \label{eq:off_diagonal_symmetric}
 {\mathbf E} \left[ \left( \sum_{i=1}^d \v_i^2 \boldsymbol{ \lambda}_i
 \right) \v_p \v_q \right] =  {\mathbf E} \left[ \left( \sum_{i=1}^d \v_i^2 \boldsymbol{ \lambda}_i
 \right) -\v_p \v_q \right] = 0
\end{align}
holds for all $\v_p \neq \v_q$. For $p=q$,
\begin{align}
{\mathbf E} \left[  \left( \sum_{i=1}^d \v_i^2 \boldsymbol{ \lambda}_i
 \right) \v_p^2 \right] = {\mathbf E} \left[ \v_p^4 \right] - ({\mathbf E} \left[\v_p^2\right])^2 + \langle {\mathbf E}\left[ \v^2 \right], \boldsymbol{ \lambda} \rangle {\mathbf E}\left[ \v_p^2 \right]
\end{align}
holds. Replacing the above result together with Eq.~\eqref{eq:off_diagonal_symmetric} into Eq.~\eqref{eq:symmetric_xMx_expansion} concludes the proof.

\end{proof}

A straightforward application of the above result concludes the following corollary on Gaussian inputs.
\begin{corollary}
\label{cor:gaussianity_consequence}
Suppose that ${\bf x}$ is drawn from a zero-mean multivariate normal distribution, i.e. ${{\bf x} \sim {\mathcal N}(0,\cov)}$, then the following holds for every matrix ${\bf M}$:
\begin{align} \label{eq:gaussianity_consequence}
    {\mathbf E} \left[ {\bf x} {\bf x}^\top {\bf M} {\bf x} {\bf x}^\top \right] = 2 \cov {\bf M} \cov + \langle {\bf M}, \cov \rangle \cov
\end{align}
where $ \langle {\bf M}, \cov \rangle$ is the inner product of vectorized ${\bf M}$ and $\cov$.
\end{corollary}
Using the result of last corollary, we can establish a lower-bound on strong grow constant in Eq.~\eqref{assum:sgc}.
\begin{lemma} \label{lemma:example1support}
    Consider the regression objective $f$ on Gaussian input ${\bf x} \sim {\mathcal N}(0,\cov)$ and suppose that $y=0$. There exists a vector ${\bf w}_0 \in {\mathbb{R}}^d$ such that
    \begin{align}
        {\mathbf E}_{\bf z} \left[\| \nabla f_{\bf z}({\bf w}_0) \|^2 \right] \geq L/\mu \| \nabla f({\bf w}_0) \|^2.
    \end{align}
\end{lemma}
\begin{proof}
Let ${\bf w}_0$ be the eigenvector of $\cov$ associated with the smallest eigenvalue of $\cov$, i.e. ${\cov {\bf w}_0 = \mu {\bf w}_0}$. Then we use the result of corollary~\ref{cor:gaussianity_consequence} to compute expected norm of stochastic gradients as
\begin{align}
    {\mathbf E} \| \nabla f_{\bf z}({\bf w}_0) \|^2 & = {\bf w}^\top_0 {\mathbf E} \left[ {\bf x} {\bf x}^\top {\bf x} {\bf x}^\top \right] {\bf w}_0 \\
    & \stackrel{\eqref{eq:gaussianity_consequence}}{=} {\bf w}_0^\top \cov^2 {\bf w}_0 + \text{Tr}[\cov]{\bf w}_0^\top \cov {\bf w}_0 \\
    & = \mu^2 + \text{Tr}[\cov] \mu \geq \mu (\mu+L)
\end{align}
The norm of gradient evaluated at ${\bf w}_0$ is
\begin{align}
    \|\nabla f({\bf w}) \|^2 = {\bf w}^\top_0 \cov^2 {\bf w}_0 = \mu^2
\end{align}
To satisfy assumption of Eq.~\eqref{assum:sgc}, we need to choose $\rho$ such that
\begin{align}
    \rho \geq 1 + L/\mu \implies   \mu (\mu+L) \leq \rho \mu^2
\end{align}
\end{proof}
Correspondingly we have a similar lemma for Example~\ref{exam:uniform_example}.
\begin{lemma}\label{lemma:uniform_consequence}
    Consider the regression objective $f$ on the distribution defined in Example~\ref{exam:uniform_example} and suppose $y=0$. There exists a vector ${\bf w}_0$ such that
    \begin{align}
        {\mathbf E}\|\nabla f_{\bf z}({\bf w}_0)\|^2 = (L/2\mu+1)\|\nabla f({\bf w}_0)\|^2.
    \end{align}
\end{lemma}
\begin{proof}
The covariance matrix of the variable ${\bf x}$ is
\begin{align}
    \cov :=\begin{bmatrix} 1 & 0 \\ 0 & \kappa^{-1}/3 \end{bmatrix}.
\end{align}
Let ${\bf w}_0$ be the eigenvector of $\cov$ associated with the smallest eigenvalue of $\cov$, i.e. $\cov{\bf w}_0={\bf w}_0$. Similar to Corollary~\ref{cor:gaussianity_consequence} we can calculate
\begin{align}
    {\mathbf E}\big[{\bf x}\x^\top{\bf x}\x^\top\big]=\begin{bmatrix}\kappa^{-1}/3+1 & 0 \\ 0 & \kappa^{-2}/5\end{bmatrix}.
\end{align}
Then the expected norm of stochastic gradients is
\begin{align}
    {\mathbf E}\|\nabla f_{\bf z}({\bf w}_0)\|^2 & = {\bf w}_0^\top{\mathbf E}\big[{\bf x}\x^\top{\bf x}\x^\top\big]{\bf w}_0 \\
    & ={\bf w}_0^\top\begin{bmatrix}\kappa^{-1}/3+1 & 0 \\ 0 & \kappa^{-2}/5\end{bmatrix}{\bf w}_0 \\
    & = \kappa^{-1}/3+1.
\end{align}
Writing in term of $\mu$ and $L$ leads to conclusion
\begin{align}
    {\mathbf E}\|\nabla f_{\bf z}({\bf w}_0)\|^2 = (L/2\mu+1)\|\nabla f({\bf w}_0)\|^2.
\end{align}
\end{proof}

\subsection{Proof of Theorem~\ref{thm:convergence_char}} \label{sec:proof_of_theorem_appendix}
\paragraph{Deterministic--stochastic decomposition.}
Recall matrix ${\bf A}_n$ in recurrence of Eq.~\eqref{eq:recurrrence_wns}:
\begin{align}
    {\bf A}_n := \begin{bmatrix}
    (1+\beta) {\mathbf I}- (1+\alpha) \gamma {\bf x}_n{\bf x}^\top_n & \alpha \gamma {\bf x}_n {\bf x}^\top_n - \beta {\mathbf I} \\
    {\mathbf I} & 0
    \end{bmatrix}.
\end{align}
${\bf A}_{n+1}$ decomposes as
 \begin{align} \label{eq:An_decomposition}
     {\bf A}_{n+1} = {\bf A} -\gamma
     \epsilon, \quad \epsilon := \begin{bmatrix}
     (1+\alpha) \Delta & -\alpha\Delta \\
     0 & 0
     \end{bmatrix}, \Delta :=(\cov - {\bf x} {\bf x}^\top)
 \end{align}
 where ${\mathbf E} \left[ \epsilon \right] = 0$ and matrix ${\bf A}$ is
 \begin{align} \label{eq:A}
     {\bf A} = \begin{bmatrix}
     {\bf A}_1 & {\bf A}_2 \\
     {\mathbf I} & 0
     \end{bmatrix}, {\bf A}_1 := (1+\beta){\mathbf I} - \gamma (1+\alpha) \cov, {\bf A}_2 := \alpha \gamma \cov - \beta {\mathbf I}
 \end{align}
 The above decomposition allows us to analyze the covariance term induced by the noise separately. Recall that matrices ${\bf M}_n$ and ${\bf B}_n$ were defined as
 \begin{align}
 {\bf M}_n = {\mathbf E} \left[ {\bf B}_n^\top {\bf B}_n \right], \quad {\bf B}_n = {\bf A}_n {\bf A}_{n-1} \dots {\bf A}_1
 \end{align}

 \begin{lemma}
Suppose assumption~\ref{assume:symmetric_inputs} holds on ${\bf x}$. We further assume that ${\bf M}$ has the following block structure:
\begin{align}
      {\bf M}_n = \begin{bmatrix}
      {\bf M}_1 & {\bf M}_2 \\
      {\bf M}_2 & {\bf M}_3
     \end{bmatrix}, \quad {\bf M}_{i} = {\bf U} {\textbf{diag}}(\boldsymbol{ \lambda}_n^{(i)}) {\bf U}^\top
\end{align}

Then matrix ${\bf M}_{n+1}$ can be decomposed as
\begin{align} \label{eq:E_nplus_expansion_result}
     {\bf M}_{n+1} = {\bf A}^\top  {\bf M}_n {\bf A} + \gamma^2 \underbrace{\begin{bmatrix} (1+\alpha)^2 {\bf K} &  -\alpha(1+\alpha) {\bf K} \\
 -\alpha (1+\alpha) {\bf K}  & \alpha^2 {\bf K}
 \end{bmatrix}}_{\text{due to the noise}}
\end{align}
where  \begin{align}
    {\bf K} := {\bf U} {\textbf{diag}}\left( \left({\textbf{diag}}(\k - 2(\boldsymbol{\sigma})^2) + \boldsymbol{\sigma} \boldsymbol{\sigma}^\top \right)\boldsymbol{ \lambda}_n^{(1)} \right) {\bf U}^\top
\end{align}
\end{lemma}
\begin{proof}
 Using the decomposition of matrix ${\bf A}_n$ in Eq.~\eqref{eq:An_decomposition}, we expand the expectation as
 \begin{align}\label{eq:M_nplus}
     {\bf M}_{n+1}  & = {\mathbf E} \left[{\bf A}_{n+1}^\top  {\bf M}_n {\bf A}_{n+1}  \right] \\
     & = {\bf A}^\top  {\bf M}_n {\bf A} + \gamma^2  {\mathbf E} \left[ \epsilon^\top  {\bf M}_n \epsilon \right]
 \end{align}

We compute 2nd term (induced by the noise)
\begin{align} \label{eq:variance_noise}
    {\mathbf E} \left[ \epsilon^\top  {\bf M}_n \epsilon \right]  = \begin{bmatrix}
    (1+\alpha)^2 {\bf K} & -\alpha (1+\alpha) {\bf K} \\
    -\alpha (1+\alpha) {\bf K}  & \alpha^2 {\bf K}
    \end{bmatrix}
\end{align}
The result of Lemma~\ref{lemma:symmetric_inputs} on symmetric inputs concludes the expression of ${\bf K}$.
\end{proof}

\paragraph{The block structure of matrix ${\bf M}_n$.}
Lemmas~\ref{lemma:Mn_blockstructure} and \ref{lemma:blambdan} are consequences of the last lemma.
\begin{corollary} [Combined lemmas \ref{lemma:Mn_blockstructure} and \ref{lemma:blambdan}]
Suppose assumption~\ref{assume:symmetric_inputs} holds on ${\bf x}$; then
\begin{align} \label{eq:Mn_condition}
      {\bf M}_n  = \begin{bmatrix}
      {\bf M}_1 & {\bf M}_2 \\
      {\bf M}_2 & {\bf M}_3
     \end{bmatrix}, \quad {\bf M}_{i} = {\bf U} {\textbf{diag}}(\boldsymbol{ \lambda}_n^{(i)}) {\bf U}^\top
\end{align}
where
\begin{align} \label{eq:an_appendix}
    \a_n=
     {\bf C}^n_\Theta
    {\bf 1}, \quad \a_n := \begin{bmatrix}
    \boldsymbol{ \lambda}_n^{(1)}
    &
    \boldsymbol{ \lambda}_n^{(2)}
    &
    \boldsymbol{ \lambda}_n^{(3)}
     \end{bmatrix}
\end{align}
and matrix ${\bf C}_\Theta$ is
\begin{align}
    {\bf C}_\Theta & = \begin{bmatrix} {\bf D}_1^2 + (1+\alpha)^2 {\bf K}' & 2 {\bf D}_1 & {\mathbf I} \\
     {\bf D}_1 {\bf D}_2 - \alpha(1+\alpha) {\bf K}'& {\bf D}_2 & 0 \\
     {\bf D}_2^2 +\alpha^2 {\bf K}' & 0 & 0
    \end{bmatrix}  \label{eq:C_matrix}\\
    {\bf D}_1 &= (1+\beta) {\mathbf I} - \gamma (1+\alpha) {\textbf{diag}}(\boldsymbol{\sigma}),\\
    {\bf D}_2 &=  \alpha \gamma {\textbf{diag}}(\boldsymbol{\sigma}) - \beta {\mathbf I} \\
    {\bf K}' &= \gamma^2 \left( {\textbf{diag}}(\k - 2(\boldsymbol{\sigma})^2) + \boldsymbol{\sigma} \boldsymbol{\sigma}^\top \right) \label{eq:cmat_appendix}
\end{align}
\end{corollary}
\begin{proof}
  We prove the result by induction.
  Recall matrix ${\bf A}$ with the following block structure:
\begin{align}
     {\bf A} = \begin{bmatrix}
     {\bf A}_1 & {\bf A}_2 \\
     {\mathbf I} & 0
     \end{bmatrix}, \quad  {\bf A}_1 := (1+\beta){\mathbf I} - \gamma (1+\alpha) \cov, \quad  {\bf A}_2 := \alpha \gamma \cov - \beta {\mathbf I}
 \end{align}
Given the decomposition $\cov = {\bf U} {\textbf{diag}}(\boldsymbol{\sigma}) {\bf U}^\top$ (in Eq.~\eqref{eq:cov}), submatrices ${\bf A}_i$ decomposes as:
\begin{align}
    {\bf A}_i = {\bf U} {\bf D}_i {\bf U}^\top
\end{align}
where ${\bf U}$ is eigenvectors of the covariance  and ${\bf D}_i$ are diagonal matrices of eigenvalues in Eq.~\eqref{eq:cmat_appendix}.  Using the assumption of  Eq.~\eqref{eq:Mn_condition},
  we expand the first term in Eq.~\eqref{eq:E_nplus_expansion_result}
\begin{align}
    {\bf A}^\top  {\bf M}_n {\bf A} & = \begin{bmatrix}
    {\bf A}_1 {\bf M}_1 {\bf A}_1 + 2{\bf M}_2 {\bf A}_1 + {\bf M}_3 & {\bf A}_1 {\bf M}_1 {\bf A}_2 + {\bf M}_2 {\bf A}_2 \\
    {\bf A}_2 {\bf M}_1 {\bf A}_1 + {\bf M}_2 {\bf A}_2 & {\bf A}_2 {\bf M}_1 {\bf A}_2
    \end{bmatrix} \\
    & =
    \begin{bmatrix}
     {\bf B}_1 & {\bf B}_2 \\
     {\bf B}_2 & {\bf B}_3
    \end{bmatrix}
\end{align}
where
\begin{align}
    {\bf B}_1 &= {\bf U} \left( {\bf D}_1^2 {\textbf{diag}}(\boldsymbol{ \lambda}_n^{(1)}) + 2 {\bf D}_1{\textbf{diag}}(\boldsymbol{ \lambda}_n^{(2)}) + {\textbf{diag}}(\boldsymbol{ \lambda}_n^{(3)})\right) {\bf U}^\top \\
    {\bf B}_2 &= {\bf U} \left( {\bf D}_1 {\bf D}_2 {\textbf{diag}}(\boldsymbol{ \lambda}_n^{(1)}) + {\bf D}_2 {\textbf{diag}}(\boldsymbol{ \lambda}_n^{(2)}) \right) {\bf U}^\top \\
    {\bf B}_3 & = {\bf U} \left(  {\bf D}_2^2 {\textbf{diag}}(\boldsymbol{ \lambda}_n^{(1)}) \right) {\bf U}^\top
\end{align}
Using eigendecomposition of $\cov$ and ${\bf M}_1$, the matrix ${\bf K}$ in the noise-induced term can be written as
\begin{align}
     {\bf K}  = {\bf U} \left(  \left( {\textbf{diag}}(\k - 2(\boldsymbol{\sigma})^2) + \boldsymbol{\sigma} \boldsymbol{\sigma}^\top \right) {\textbf{diag}}(\boldsymbol{ \lambda}_n^{(1)}) \right) {\bf U}^\top
\end{align}
Putting all together, we complete the proof:
\begin{align}
     {\bf M}_{n+1}  = \begin{bmatrix}
    {\bf M}_{1}^+ & {\bf M}_2^+ \\
    {\bf M}_2^+ &  {\bf M}_3^+
    \end{bmatrix}, \quad  {\bf M}_i^+ = {\bf U} {\textbf{diag}}(\boldsymbol{ \lambda}_{n+1}^{(i)}) {\bf U}^\top
\end{align}
where
\begin{align}
    \a_{n+1} = {\bf C}_{\Theta}\a_n
\end{align}
where matrix ${\bf C}_{\Theta}$ is those of Eq.~\eqref{eq:C_matrix}.
The above results concludes our induction. \end{proof}
\paragraph{Bound spectral norm of ${\bf M}_n$.}
Given bounds on $\a_n$, how we can establish lowerbound and upperbound on eigenvalues of ${\bf M}_n$? Next lemma address this result. The result will be used in the proof of Theorem~\ref{thm:convergence_char} (very last step).
\begin{lemma} \label{lemma:v1_times_Mn_bound}
Recall matrix ${\bf M}_n$ in the last corollary:
\begin{align}
      {\bf M}_n = \begin{bmatrix}
      {\bf M}_1 & {\bf M}_2 \\
      {\bf M}_2 & {\bf M}_3
     \end{bmatrix}, \quad {\bf M}_{i} = {\bf U} {\textbf{diag}}(\boldsymbol{ \lambda}_n^{(i)}) {\bf U}^\top, \quad \a_n = \begin{bmatrix}
    \boldsymbol{ \lambda}_n^{(1)} &
    \boldsymbol{ \lambda}_n^{(2)} &
    \boldsymbol{ \lambda}_n^{(3)}
    \end{bmatrix}
\end{align}
For ${\bf M}_n$, the following holds
\begin{align}
   \| \a_n \|_\infty \| \v_1 \|^2 \leq  \begin{bmatrix}
     \v_1^\top & \v_1^\top
    \end{bmatrix} {\bf M}_n \begin{bmatrix}
    \v_1 \\
    \v_1
    \end{bmatrix}
    \leq  6 \sqrt{d} \|\a_n\| \| \v_1 \|^2.
\end{align}
\end{lemma}
\begin{proof}
The proof is straight-forward
 \begin{multline} \label{eq:v1_Mn_v1_expansion}
      \begin{bmatrix}
    \v_1^\top & \v_1^\top
    \end{bmatrix} \begin{bmatrix}
    {\bf U} {\textbf{diag}}(\boldsymbol{ \lambda}_n^{(1)}) {\bf U}^\top & {\bf U} {\textbf{diag}}(\boldsymbol{ \lambda}_n^{(2)}) {\bf U}^\top \\
    {\bf U} {\textbf{diag}}(\boldsymbol{ \lambda}_n^{(2)}) {\bf U}^\top & {\bf U} {\textbf{diag}}(\boldsymbol{ \lambda}_n^{(3)}) {\bf U}^\top
    \end{bmatrix}
    \begin{bmatrix}
    \v_1 \\
    \v_1
    \end{bmatrix} \\
    = \v_1^\top {\bf U}^\top {\textbf{diag}}\left( \boldsymbol{ \lambda}_n^{(1)} + 2 \boldsymbol{ \lambda}_n^{(2)} + \boldsymbol{ \lambda}_n^{(3)} \right) {\bf U} \v_1
 \end{multline}
Hence
\begin{align}
    \begin{bmatrix}
     \v_1^\top & \v_1^\top
    \end{bmatrix} {\bf M}_n \begin{bmatrix}
    \v_1 \\
    \v_1
    \end{bmatrix} & \leq 2 \left( \sum_{i=1}^3 \|\boldsymbol{ \lambda}_n^{(i)}\|_1\right)  \| {\bf U} \v_1 \|^2 \\
    & \leq 2 \sqrt{3} \left( \sum_{i=1}^d \| \boldsymbol{ \lambda}_n^{(i)}\| \right) \| \v_1 \|^2  \\
    & \leq 6 \sqrt{d} \|\a_n\|_2 \| \v_1 \|^2
\end{align}
which concludes the desired upper-bound. For the lower-bound, we use the fact that all $\boldsymbol{ \lambda}_n^{(i)}$ have positive coordinates (as they are eigenvalues of symmetric matrices). Plugging this into the Eq.~\eqref{eq:v1_Mn_v1_expansion} together with orthogonality of ${\bf U}$ concludes the proof:
\begin{align}
    \| \a_n \|_\infty \| \v_1 \|^2 \leq  \begin{bmatrix}
     \v_1 & \v_1^\top
    \end{bmatrix} {\bf M}_n \begin{bmatrix} \v_1^\top \\ \v_1 \end{bmatrix}
\end{align}

\end{proof}
\subsection{The spectral analysis} \label{sec:spectral_analysis}
\begin{lemma}[Restated Lemma~\ref{lemma:spectral_radius_bound}]
The spectral radius of matrix ${\bf C}_{\gamma, \alpha,\beta}$ is bounded as
\begin{align}
    \| {\bf C}_{\gamma,\alpha,\beta} \|_{\rho} \leq \max_{i=1,\dots,d} \| {\mathbf J}_i \|_\rho + \epsilon + 3(1+\alpha)^2 \gamma^2   \| {\textbf{diag}}(\boldsymbol{\sigma})^2 - \boldsymbol{\sigma} \boldsymbol{\sigma}^\top \|
\end{align}
where ${\mathbf J}_j$ is a $3 \times 3$ matrix as
\begin{align}
    {\mathbf J}_{i} := \begin{bmatrix}
      [{\bf D}_1]_{ii}^2 + (1+\alpha)^2 \gamma^2 (k_{i} - \sigma^2_i)   &  2[{\bf D}_1]_{ii}  & 1 \\
    [{\bf D}_1]_{ii}[{\bf D}_2]_{ii} -  \alpha (1+\alpha)  \gamma^2 (k_i - \sigma_i^2)  & [{\bf D}_{2}]_{ii} & 0 \\
    [{\bf D}_{2}]_{ii}^2 + \gamma^2 \alpha^2 (k_i - \sigma_i^2) & 0 & 0
    \end{bmatrix}
\end{align}
\end{lemma}
\begin{proof}
Recall matrix ${\bf C} := {\bf C}_{\gamma,\alpha,\beta}$ in Eq.~\eqref{eq:cmatrix} formulation:

\begin{align}
    {\bf C} : = \begin{bmatrix} {\bf C}_1 & 2 {\bf D}_1 & {\mathbf I} \\
     {\bf C}_2 & {\bf D}_2 & 0 \\
     {\bf C}_3 & 0 & 0
    \end{bmatrix}
\end{align}
where
\begin{align}
    {\bf C}_1 & = {\bf D}_1^2 + (1+\alpha)^2 {\bf K}', \quad {\bf D}_1 := (1+\beta) {\mathbf I} - \gamma (1+\alpha) {\textbf{diag}}(\boldsymbol{\sigma})\\
    {\bf C}_2 & = {\bf D}_1 {\bf D}_2 - \alpha(1+\alpha) {\bf K}', \quad {\bf D}_2 =  \alpha \gamma {\textbf{diag}}(\boldsymbol{\sigma}) - \beta  {\bf K}'\\ {\bf C}_3 & = {\bf D}_2^2 +\alpha^2 {\bf K}', \quad {\bf K}' = \gamma^2 \left( {\textbf{diag}}(\k - 2(\boldsymbol{\sigma})^2) + \boldsymbol{\sigma} \boldsymbol{\sigma}^\top \right)
\end{align}
\paragraph{A perturbation form.}
We can decompose ${\bf K}'$ (the induced matrix by noise) to sum of diagonal matrices and a non-diagonal matrix:
\begin{align}
    {\bf K}' := \gamma^2 \left( {\textbf{diag}}(\k- \boldsymbol{\sigma}^2) + \xi \right), \quad \xi = \boldsymbol{\sigma} \boldsymbol{\sigma}^\top  - {\textbf{diag}}(\boldsymbol{\sigma})^2
\end{align}
Using the decomposition of ${\bf K}'$, we decompose ${\bf C}$ as a perturbation of a block diagonal matrix:
\begin{align}
    {\bf C} = \bar{{\bf C}} + \gamma^2 \boldsymbol{\xi}, \quad \boldsymbol{\xi} := \begin{bmatrix}
    (1+\alpha)^2 \xi & 0 & 0 \\
    -\alpha(1+\alpha) \xi & 0 & 0 \\
    \alpha^2 \xi & 0 & 0
    \end{bmatrix}
\end{align}
where matrix $\bar{{\bf C}}$ is
\begin{align}
    \bar{{\bf C}} := \begin{bmatrix} {\bf D}_1^2 + (1+\alpha)^2 \bar{{\bf K}}   & 2 {\bf D}_1 & {\mathbf I} \\
     {\bf D}_1 {\bf D}_2 - \alpha(1+\alpha) \bar{{\bf K}}& {\bf D}_2 & 0 \\
     {\bf D}_2^2 +\alpha^2 \bar{{\bf K}} & 0 & 0
    \end{bmatrix}, \quad \bar{{\bf K}} :=  \gamma^2 \left( {\textbf{diag}}(\k)-{\textbf{diag}}(\boldsymbol{\sigma})^2  \right)
\end{align}
\paragraph{Exploiting the diagonal structure.} Matrix $\bar{{\bf C}}$ has a particular structure: its $d\times d$ blocks are diagonal. Similar to the classical analysis of Heavy ball method \cite{stephen2017behavior}, we permute rows and columns to compute eigenvalues of $\bar{{\bf C}}$. Let $\Pi$ be a permutation matrix that swaps column(and row) $i+d$ and $i+2d$ with columns $i+1$ and $i+2$, respectively. Then,
\begin{align}
    \Pi^\top \bar{{\bf C}} \Pi  = \begin{bmatrix}
     {\mathbf J}_{1} & 0&  0 & \dots & 0 \\
     0 & {\mathbf J}_{2} & 0 & \dots & 0 \\
     \vdots & \vdots & \vdots & \vdots & \vdots \\
     0 & \dots & \dots & 0 & {\mathbf J}_{d}
    \end{bmatrix}
\end{align}
where where ${\mathbf J}_j$ is a $3 \times 3$ matrix as
\begin{align}
    {\mathbf J}_{i} := \begin{bmatrix}
      [{\bf D}_1]_{ii}^2 + (1+\alpha)^2 \gamma^2 (k_{i} - \sigma^2_i)   &  2 [{\bf D}_1]_{ii}  & 1 \\
    [{\bf D}_1]_{ii}[{\bf D}_2]_{ii} -  \alpha (1+\alpha)  \gamma^2 (k_i - \sigma_i^2)  & [{\bf D}_{2}]_{ii} & 0 \\
    [{\bf D}_{2}]_{ii}^2 + \gamma^2 \alpha^2 (k_i - \sigma_i^2) & 0 & 0
    \end{bmatrix}
\end{align}
Given spectral decomposition of ${\mathbf J}_i = {\bf U}_i {\textbf{diag}}(\d_i) {\bf U}_i^\top$, we decompose $\bar{{\bf C}}$:
\begin{align}
   {\bf U}_J^\top \Pi^\top  \bar{{\bf C}} \Pi {\bf U}_J = \begin{bmatrix}
   \d_1 & 0 & \dots & 0 \\
    0 & \d_2 & \dots & 0 \\
    \vdots & \vdots & \vdots & \vdots \\
    0 & 0 & \dots & \d_d
   \end{bmatrix}
   , \quad {\bf U}_J :=  \begin{bmatrix}
    {\bf U}_1 & 0 & \dots & 0 \\
    0 & {\bf U}_2 & \dots & 0 \\
    \vdots & \vdots & \vdots & \vdots \\
    0 & 0 & \dots & {\bf U}_d
    \end{bmatrix}
\end{align}
Since matrices ${\bf U}_j$ and perturbation matrix $\Pi$ are orthogonal, we conclude that eigenvalues of $\bar{{\bf C}}$ are those of ${\mathbf J}_i$s (i.e. $\d_i$s). Hence
\begin{align} \label{eq:Cbar_bound}
    \| \bar{{\bf C}} \|_\rho \leq \max_{i=1,\dots, d} \| {\mathbf J}_i \|_\rho
\end{align}
\paragraph{A bound on pseudospectrum.}
So far, we have established a bound on spectral radius of $\bar{{\bf C}}$. Yet, we need to bound the spectral radius of $\bar{{\bf C}} + \gamma^2 \boldsymbol{\xi}$. To this end, we use results of results Pseudospectrum in section \ref{sec:pseudospectrum}.
\begin{align}
    \| {\bf C} \|_{\rho_\epsilon} & = \| \bar{{\bf C}} + \gamma^2 \boldsymbol{\xi} \|_{\rho_\epsilon} \\
    & \stackrel{\text{Lemma}~\ref{sec:robust_pseudo}}{\leq } \| \bar{{\bf C}} \|_{\rho_{\epsilon + \gamma^2 \| \boldsymbol{\xi}\|}} \\
    & \stackrel{\text{Lemma}~\ref{lemma:baur_fike}}{\leq}  \| \bar{{\bf C}} \|_\rho  + (\epsilon + \gamma^2 \| \boldsymbol{\xi} \|)
\end{align}
Note that in the last step, we have used the fact that the conditioning of $\Pi {\bf U}_J$ is bounded by 1 (which is a consequence of orthogonality of matrices $\Pi$ and ${\bf U}_J$). Replacing the above result into Eq.~\eqref{eq:Cbar_bound} concludes the proof:
\begin{align}
    \| {\bf C} \|_{\rho_\epsilon} \leq  \max_{i=1,\dots, d} \| {\mathbf J}_i \|_\rho + \epsilon + \gamma^2 \| \boldsymbol{\xi} \| \leq  \max_{i=1,\dots, d} \| {\mathbf J}_i \|_\rho + \epsilon + 3(1+\alpha)^2 \gamma^2 \| \xi \|
\end{align}
\end{proof}

\subsection{Examples of improved rates} \label{sec:converge_examples_app}
\begin{lemma}[Restated Lemma \ref{lemma:gaussian_gaurantees}]
Suppose input and label distributions are those of example~\ref{exam:gaussian_example}.
For $\mu\leq0.02$, consider stochastic accelerated gradient descent with parameters: $\alpha = 2$, $\beta=1-10^{-1/2}\sqrt{\mu}$ and $\gamma = 0.1$. Then,
\begin{align}
    {\mathbf E} \| {\bf w}_n - {\bf w}_* \|^2 \leq \frac{1200}{\sqrt{\mu}} \left( 1- \sqrt{\mu}/5  \right)^n \| {\bf w}_0 - {\bf w}_* \|^2.
\end{align}
\end{lemma}

\begin{proof}
Recall ${\bf x}\sim\mathcal{N}(0,{\textbf{diag}}([\mu,1]))$, as defined in Example \ref{exam:gaussian_example}, and $\boldsymbol{\sigma}=[\sigma_1,\sigma_2]=[\mu,1]$. Therefore their corresponding fourth moments are
\begin{align}
    k_1=3\mu^2,\quad k_2=3.
\end{align}
Plug $k_i$ and $\sigma_i$ into ${\mathbf J}_i$
\begin{align}
    {\mathbf J}_{i} = \begin{bmatrix}
      [{\bf D}_1]_{ii}^2 + 2(1+\alpha)^2 \gamma^2 \sigma^2_i   &  2 [{\bf D}_1]_{ii}  & 1 \\
    [{\bf D}_1]_{ii}[{\bf D}_2]_{ii} -  2\alpha (1+\alpha)  \gamma^2 \sigma_i^2  & [{\bf D}_{2}]_{ii} & 0 \\
    [{\bf D}_{2}]_{ii}^2 + 2\gamma^2 \alpha^2 \sigma_i^2 & 0 & 0
    \end{bmatrix}.
\end{align}
Employing MATLAB symbolic tools, we can check that
\begin{align}
 \| {\mathbf J}_2 \|_\rho \leq 0.966
\end{align}
holds for our choice of parameters. Furthermore, the result of
Lemma~\ref{lemma:spectral_radius_mu} guarantees
\begin{align}
    \|{\mathbf J}_1\|_\rho\leq 1-10^{-1/2}\sqrt{\mu}.
\end{align}
Moreover, we calculate $\boldsymbol{\xi}$ as
\begin{align}
    \boldsymbol{\xi} & =\boldsymbol{\sigma}\bsigma^\top-{\textbf{diag}}(\boldsymbol{\sigma})^2 =\begin{bmatrix} \mu^2 & \mu \\ \mu & 1 \end{bmatrix} - \begin{bmatrix} \mu^2 & 0 \\ 0 & 1 \end{bmatrix}=\begin{bmatrix} 0 & \mu \\ \mu & 0 \end{bmatrix}
\end{align}
whose nor is bounded as
\begin{align}
    \|\boldsymbol{\xi}\|\leq\mu.
\end{align}
$\mu\leq0.02$ concludes the proof:
\begin{align}
    \|{\bf C}\|_{\rho_\epsilon} & \leq \max\{1-10^{-1/2}\sqrt{\mu},0.966\}+\epsilon+0.27\mu \\
    & \leq 1-\sqrt{\mu}/4+\epsilon.
\end{align}
Choosing $\epsilon=0.05\sqrt{\mu}$ concludes the proof.
\end{proof}

\begin{lemma}[Restated Lemma \ref{lemma:uniform_gaurantees}] \label{lemma:uniform_convergence_app}
For $\kappa\leq0.02$, consider running stochastic acceleration method on example~\eqref{exam:uniform_example}. If $\alpha=2$, $\beta=1-10^{-1/2}\sqrt{\kappa}$ and $\gamma=\kappa/10$, then
\begin{align}
    {\mathbf E} \| {\bf w}_n - {\bf w}_* \|^2 \leq \frac{1200}{\sqrt{\kappa}} \left( 1- \sqrt{\kappa}/5 \right)^n \| {\bf w}_0 - {\bf w}_* \|^2.
\end{align}
holds.
\end{lemma}

\begin{proof}
The first coordinate of ${\bf x}$ is a Rademacher random variable with fourth moment $k_1=1/4$ and $\sigma_1=1/2$. For the second coordinate, the moments of uniform distribution on range $[-\kappa^{-1/2},\kappa^{-1/2}]$ are
\begin{align}
    k_2=\int_{-\kappa^{-1/2}}^{\kappa^{-1/2}}\frac{x^4}{2\kappa^{-1/2}}dx=\kappa^{-2}/5
\end{align}
and
\begin{align}
    \sigma_2=\int_{-\kappa^{-1/2}}^{\kappa^{-1/2}}\frac{x^2}{2\kappa^{-1/2}}dx=\kappa^{-1}/3.
\end{align}
Plug $k_i$'s and $\sigma_i$'s into ${\mathbf J}_i$ as defined in Lemma \ref{lemma:spectral_radius_bound}
\begin{align}
    {\mathbf J}_{1} := \begin{bmatrix}
      [{\bf D}_1]_{11}^2   &  2[{\bf D}_1]_{11}  & 1 \\
    [{\bf D}_1]_{11}[{\bf D}_2]_{11}   & [{\bf D}_{2}]_{11} & 0 \\
    [{\bf D}_{2}]_{11}^2  & 0 & 0
    \end{bmatrix}
\end{align}
and
\begin{align}
    {\mathbf J}_{2} := \begin{bmatrix}
      [{\bf D}_1]_{22}^2 + 4(1+\alpha)^2 \gamma^2 \kappa^{-2}/45   &  2[{\bf D}_1]_{22}  & 1 \\
    [{\bf D}_1]_{22}[{\bf D}_2]_{22} -  4\alpha (1+\alpha)  \gamma^2  \kappa^{-2}/45   & [{\bf D}_{2}]_{22} & 0 \\
    [{\bf D}_{2}]_{2}^2 + 4\gamma^2 \alpha^2  \kappa^{-2}/45  & 0 & 0
    \end{bmatrix}.
\end{align}

For the Rademacher coordinate, the additional noise term cancels out. Clearly $-[{\bf D}_2]_{11}$ is an eigenvalue and a simple deduction from the proof for Lemma \ref{lemma:spectral_radius_mu} yields $\|{\mathbf J}_1\|_\rho \leq 1-10^{-1/2}\sqrt{\kappa}$. For ${\mathbf J}_2$, we employ MATLAB to verify following fact.
\paragraph{Fact 3.}
\textit{
For $\kappa\leq0.02,\alpha=2,\beta=1-10^{-1/2}\sqrt{\kappa}$ and $\gamma=\kappa/10$, $\|{\mathbf J}_2\|_\rho\leq0.965.$}

Next we calculate $\boldsymbol{\xi}$,
\begin{align}
    \boldsymbol{\xi} & =\boldsymbol{\sigma}\bsigma^\top-{\textbf{diag}}(\boldsymbol{\sigma})^2 \\
    & =\begin{bmatrix} 1 & \kappa^{-1}/3 \\ \kappa^{-1}/3 & \kappa^{-2}/9 \end{bmatrix} - \begin{bmatrix} 1 & 0 \\ 0 & \kappa^{-2}/9 \end{bmatrix}=\begin{bmatrix} 0 & \kappa^{-1}/3 \\ \kappa^{-1}/3 & 0 \end{bmatrix}
\end{align}
whose norm is
$\|\boldsymbol{\xi}\| = \kappa^{-1}/3$.
And $\|{\bf C}\|_{\rho_\epsilon}$ is upper bounded as
\begin{align}
    \|{\bf C}\|_{\rho_\epsilon} & \leq \max\{1-10^{-1/2}\sqrt{\kappa},0.95\}+\epsilon+0.09\kappa \\
    & \leq 1-\sqrt{\kappa}/4+\epsilon.
\end{align}
Choosing $\epsilon = 0.05\sqrt{\kappa}$ concludes the proof.
\end{proof}

\subsubsection{Spectral radius at $\mu$}\label{sec:spectral_radius_mu}
 Here we establish an upper bound for the spectral radius of $3\times3$ matrix ${\mathbf J}_1$ when $\sigma_1=\mu$. Subscript $i$ are temporarily omitted without ambiguity.

\paragraph{Root of characteristic equation.}
In the rest part of this subsection, we will fix $\sigma=\mu$ and set $\alpha,\gamma$ and $\beta$ as defined in Lemma \ref{lemma:gaussian_gaurantees}, on which all claims and lemmas mentioned below are based without special statements. We denote the entries of ${\mathbf J}$ as
\begin{align}
    {\mathbf J} = \begin{bmatrix}
      D_1^2 + 2(1+\alpha)^2 \gamma^2 \mu^2   &  2 D_1  & 1 \\
    D_1D_2 -  2\alpha (1+\alpha)  \gamma^2 \mu^2  & D_{2} & 0 \\
    D_{2} + 2\gamma^2 \alpha^2 \mu^2 & 0 & 0
    \end{bmatrix}.
\end{align}
where
\begin{align}
    D_1=(1+\beta)-(1+\alpha)\gamma\mu,\quad D_2=\alpha\gamma\mu-\beta.
\end{align}
Now, consider the characteristic equation of ${\mathbf J}$
\begin{align}\label{eq:J_charact}
    x^3+bx^2+cx+d=0
\end{align}
where
\begin{align}
    b &{} :=-D_1^2 + 2(1+\alpha)^2 \gamma^2 \mu^2-D_{2} \\
    c & := 2(1+\alpha)^2 \gamma^2 \mu^2D_{2} - D_{2} - 2\gamma^2 \alpha^2 \mu^2 - D_1^2D_2 +  4\alpha (1+\alpha)  \gamma^2 \mu^2 D_1 \\
    d & :=D_{2}\big(D_{2} + 2\gamma^2 \alpha^2 \mu^2\big).
\end{align}
 The following fact about the roots of Eq. (\ref{eq:J_charact}) can be verified by MATLAB symbolic tools.
\paragraph{Fact 1.}
\textit{For $\mu\leq0.1$, the discriminant of Eq. (\ref{eq:J_charact}), i.e. $\Delta=b^2c^2-4c^3-4b^3d-27d^2+18bcd$, is positive. Therefore the characteristic equation of ${\mathbf J}$ has one real root $x_1$ and two conjugate complex roots $x_2,x_3$.}

We turn our focus into the real root $x_1$. In Lemma \ref{lemma:spectral_radius_mu}, we extend our bounds to absolut values of complex roots. Due to the complexity of closed-form roots of the above cubic equation, we  approximate them.
\paragraph{Main idea.}
The key idea is based on viewing ${\mathbf J}$ as a perturbation of matrix ${\mathbf J}_s$ which reads as
\begin{align}
    {\mathbf J}_s = \begin{bmatrix}
      D_1^2   &  2 D_1  & 1 \\
    D_1D_2  & D_{2} & 0 \\
    D_{2}  & 0 & 0
    \end{bmatrix}.
\end{align}
For small $\gamma$, we expect that ${\mathbf J}_s-{\mathbf J}$ be close to zero. The characteristic equation of ${\mathbf J}_s$ is
\begin{align}\label{eq:Js_charact}
    y^3+b_sy^2+c_sy+d_s=0
\end{align}
where
\begin{align}
    b  :=-D_1^2 -D_{2}, \quad  c  := - D_{2} - D_1^2D_2, \quad d :=D_{2}^2 .
\end{align}

Notice that $y=-D_2$ is a root of Eq. (\ref{eq:Js_charact}). These simple solutions provide us proper estimates of the real root of Eq. (\ref{eq:J_charact}).

Let $z=x-y$. Plugging $x=y+z$ and $y=-D_2$ into Eq. (\ref{eq:J_charact}) and Eq. (\ref{eq:Js_charact}) respectively and subtracting from both sides leads to a cubic equation about $z$
\begin{align}\label{eq:z_cubic}
    z^3+pz^2+qz+r=0
\end{align}
where
\begin{align}
    p := b-3D_2, \quad q :=-2bD_2+3D_2^2+c, r :=(b-b_s)D_2^2-(c-c_s)D_2+d-d_s.
\end{align}
Let  $z_1,z_2$ and $z_3$ be roots of above cubic equation. These roots relate to those of Eq. (\ref{eq:J_charact}) as
$
 x_i =-D_2+z_1
$
for $i=1,2,3$. A natural consequence of this is $z_1$ is real and $z_2$, $z_3$ are conjugate complex (considering that $D_2$ is real).
\paragraph{Properties of the real root.}
Let's focus on the real root $z_1$. MATLAB verification indicates coefficients $p$, $q$ and $r$ obey the following property.

\paragraph{Fact 2.}
\textit{
For $\mu\leq0.1$, coefficients in Eq. (\ref{eq:z_cubic}) satisfy $p,\ q> 0$ and $r<0$.
}

Given the above fact, we prove that the real root $z_1$ is positive.
\begin{lemma}
For $\mu\leq0.1$, the real root $z_1$ of Eq.  (\ref{eq:z_cubic}) is positive.
\end{lemma}

\begin{proof}
Consider three roots $z_1,\ z_2$ and $z_3$ of Eq. (\ref{eq:z_cubic}). We calculate
\begin{align}
    (z-z_1)(z-z_2)(z-z_3) & =z^3-(z_1+z_2+z_3)z^2+(z_1z_2+z_2z_3+z_3z_1)z-z_1z_2z_3 \\
    & =z^3+pz^2+qz+r
\end{align}
which boils down to $z_1z_2z_3=-r$. The existence of zero roots is ruled out since $r\ne0$. The last fact implies that
\begin{align}
    z_1=\frac{z_1z_2z_3}{|z_2|^2}=-\frac{r}{|z_2|^2}>0
\end{align}
since $z_2$ and $z_3$ are conjugate complex roots.
\end{proof}

\paragraph{Cubic root approximation.}
Now, we establish an upper-bound on $z_1$.
\begin{lemma}
For $\mu\leq0.1$, $0\leq z_1\leq-r/q$ holds.
\end{lemma}
\begin{proof}
The last fact implies $p>0,q>0,r>0$ and $z_1>0$. By rearranging of terms in Eq. (\ref{eq:z_cubic}), we have
\begin{align}\label{eq:z_quad}
    pz^2+q^z=-(r+z_1^3)
\end{align}
of which $z_1$ is still a root. Since $z_1>0$ then $d+z_1^3=-bz_1^2-cz_1<0$ holds.
\begin{align}
    z_1\overset{\textrm{(A)}}{\leq} \frac{-q+\sqrt{q^2-4p(d+z_1^3)}}{2p}\overset{\textrm{(B)}}{<}\frac{-q+\sqrt{q^2-4pd}}{2p}\leq-\frac{r}{q}
\end{align}
where (A) is due to $0<-4p(d+z_1^3)<-4pd$ and (B) comes from inequality $\sqrt{1+x}\leq1+x/2$ for $x\geq 0$.
\end{proof}

Since both $q$ and $r$ are polynomials in $\sigma=\mu$ and other parameters, a rational-form upper bound on $z_1$ can be established. This allows establishing a bound on $x_1=z_1-D_2$ by analyzing $m$ and $n$.

\begin{lemma}\label{lemma:upper_bound_x}
  For $\mu\leq0.1$, $-D_2\leq x_1\leq1-10^{-1/2}\sqrt{\mu}$ holds.
\end{lemma}
\begin{proof}
The lower bound is immediate. Let's consider the compact notation $\theta=1-\beta$. Our choice of parameters leads to $\theta^2=0.1\mu\leq0.01$. Then, we write $q,r$ as polynomials in $\theta$
\begin{align}
    q & = 54\theta^6 + 51\theta^5 - 23\theta^4 - 9\theta^3 + 3\theta^2\geq2\theta^4>0
\end{align}
and
\begin{align}
    r & = - 24\theta^7 - 4\theta^6 + 16\theta^5 - 4\theta^4\geq- 4\theta^4<0.
\end{align}
Now we calculate
\begin{align}
    z_1\leq-r/q=\frac{4\theta^4}{2\theta^2}\leq2\theta^2.
\end{align}
This leads to an upper bound on $x_1$
\begin{align}
    x_1=-D_2+z_1\leq 1-\theta-2\theta^2+2\theta^2\leq1-\theta=1-10^{-1/2}\sqrt{\mu}.
\end{align}
\end{proof}

\paragraph{Bound for spectral radius.}
So far, we have proven the real eigenvalue is bounded by $1-\mathcal{O}(\sqrt{\mu})$ It remains to bound complex eigenvalues.
\begin{lemma} \label{lemma:spectral_radius_mu}
For $\mu\leq0.1$, $\alpha=2$, $\beta=1-10^{-1/2}\sqrt{\mu}$ and $\gamma=0.1$, $\|{\bf M}\|_\rho\leq1-10^{-1/2}\sqrt{\mu}$.
\end{lemma}
\begin{proof}
We reuse the notation $\theta=1-\beta$. The determinant of ${\mathbf J}$ read as
\begin{align}
    \det{\mathbf J} & =1-3\theta-3\theta^2+11\theta^3+14\theta^4-20\theta^5-24\theta^6
\end{align}
is upper bounded by $|D_2|(1-\theta)^2$, which can be verified by MATLAB or manual calculation.
On the other hand, we have $x_2x_3=|x_2|^2>0$ since $x_2$ and $x_3$ are conjugate. Consider $x_1>-D_2>0$ and calculate the spectral radius
\begin{align}
    \|{\mathbf J}\|_\rho & =\max\{x_1,|x_2|\}\max\{x_1,\sqrt{\frac{\det{\mathbf J}}{x_1}}\} \\
    &\leq\max\{x_1,\sqrt{\frac{|D_2|(1-\theta)^2}{|D_2|}}\}\overset{\textrm{(A)}}{\leq}1-\theta=1-10^{-1/2}\sqrt{\mu}
\end{align}
where (A) comes from $x_1\leq1-\theta$ as shown in Lemma \ref{lemma:upper_bound_x}.
\end{proof}

\section{Clarifying collaborations}
\section{Notes on the spectral radius}
Let consider the following 2d matrix
\begin{align}
    {\bf M} = \begin{bmatrix}
    1+\beta - \gamma \lambda_i & - \beta  \\
    1 & 0
    \end{bmatrix}
\end{align}
To compute eigenvalues of the above matrix, we need to compute roots of the charat. function as
\begin{align}
((1+\beta) - \gamma \lambda_i - x)x + \beta  = 0
\end{align}
After a rearrangement of the terms, we have
\begin{align}
    x^2 - (1+\beta - \gamma \lambda_i) x - \beta = 0
\end{align}
The roots are
\begin{align}
    {\bf x}_{\mp} = \frac{(1+\beta - \gamma \lambda_i) \mp \sqrt{(1+\beta - \gamma \lambda_i)^2 - 4 \beta}}{2}
\end{align}
Suppose that $\beta$ is large enough such that the square root term gets negative, then solutions would be complex in the following form: \begin{align}
     {\bf x}_{\mp} = \frac{(1+\beta - \gamma \lambda_i) \mp i \sqrt{4 \beta-(1+\beta - \gamma \lambda_i)^2 }}{2}
\end{align}
Finally, the absolute value of the the above complex numbers is
\begin{align}
    4(|{\bf x}_{\mp}|)^2 & = (1+\beta - \gamma \lambda_i)^2 + (4 \beta - (1+\beta - \gamma \lambda_i)^2)  = 4 \beta
\end{align}

\section*{References}
\bibliographystyle{plain}

\section{Introduction}

\section{Preliminaries}
\subsection{Assumptions on the data}
\subsection{Gaussian inputs}
\subsection{Rotation of symmetric random variables*}

\section{Acceleration in realizable case}
\subsection{Results on 2d Gaussian}
\subsection{Results on d Gaussian}
\subsection{Results on d symmetric*}
\subsection{Comparison: why Gaussian is not easy?}
\subsection{Proof sketch}

\section{Acceleration in Markov chain}
\subsection{Stochastic optimization and Markov-Chain}
\subsection{Invariant distribution}
\subsection{Acceleration of discrete Markov Chains}
\subsection{Acceleration of continuous Markov Chains}

\section{Experiments}
\subsection{Least-squares: Gaussian input, non-Gaussian inputs, real datasets}
\subsection{Logistic-regression: Gaussian input, non-Gaussian inputs, real datasets }

\subsection{Least squares regression}
Let's consider least-squares regression, namely the following minimization
\begin{align}
   \arg \min_{\bf w} \left( \hat{f}({\bf w}) = \frac 12 {\mathbf E}_{{\bf z}} \|y-{\bf w}^\top {\bf x} \|^2 \right), \quad {\bf z}:= ({\bf x},y) \sim \P,
\end{align}
where ${\bf x} \in {\mathbb{R}}^d$ is the input variable the  and $y \in {\mathbb{R}}$ is the response variable. We can equivalently minimize
\begin{align} \label{eq:f}
    f({\bf w}) = \frac{1}{2} {\bf w}^\top \cov {\bf w} - {\mathbf E} \left[ y {\bf x}^\top {\bf w} \right], \quad \cov := {\mathbf E}\left[ \cov \right]
\end{align}
The gradient and stochastic gradient of the above objective are
\begin{align}
    \nabla f({\bf w}) = \cov {\bf w} - {\mathbf E} \left[ y {\bf x} \right], \nabla f_{\bf z}({\bf w}) = {\bf x} {\bf x}^\top {\bf w} - y{\bf x}.
\end{align}
Let $\mu = \lambda_{\min}(\cov)$ and $L=\lambda_{\max}(\cov)$. The condition number of the matrix $\cov$, denoted by $\kappa=L/\mu$, determines the convergence of GD on $f$. GD enjoys a linear convergence of $(1-c/\kappa)$ on $f$.
\subsection{Acceleration for stochastic optimization}
 \paragraph{Heavy ball (HB) method} HB method is an accelerated method designed particularly for \textbf{quadratic} convex objectives. This method optimizes $f$ through the following recurrence: \begin{align}
      {\bf w}_{n+1} = {\bf w}_n - \gamma \nabla f({\bf w}_n) + \beta ({\bf w}_n - {\bf w}_{n-1}), \quad {\bf w}_1 = {\bf w}_0.
 \end{align}
 For $\beta = (\kappa-1)/(\kappa+1)$ and $\gamma=\min\{4/(L+\mu)^2,1/L\}$, HB converges linearly to the minimizer ${\bf w}^*$ (of $f$) \cite{jung2017fixed}:
 \begin{align}
     \| {\bf w}_n - {\bf w}^* \|^2 \leq C(\kappa)n \left(\frac{\sqrt{\kappa}-1}{\sqrt{\kappa}+1}\right)^n \| {\bf w}_0 - {\bf w}^*\|^2,
 \end{align}
 where $C(\kappa) = 4 (2+2\beta^2+4\gamma)/(\beta^{-1}) $
 (See \cite{recht2010}.).
 \paragraph{Nesterov Acceleration}
 Although HB successfully accelerates the convergence of GD on convex quadratics, there is no global convergence guarantee for this method on non-quadratic functions \footnote{There are global convergence guarantees for HB but these guarantees does not show any acceleration compared to GD.}. Nesterov adapted HB to reach a globally accelerated rate for convex and strongly convex objectives. Nesterov's proposal is based on two sequence of iterations as
 \begin{align}
     {\bf w}_{n+1} & = \v_n - \gamma \nabla f(\v_n) \\
     \v_{n+1} & = {\bf w}_{n+1} + \beta ({\bf w}_{n+1}- {\bf w}_n)
 \end{align}
 We can show that ${\bf w}_{n+1}$ converges to ${\bf w}^*$ as
 \begin{align}
     f({\bf w}_{n}) - f({\bf w}^*) \leq (L+\mu) \exp(-(n-1)/\sqrt{\kappa})\| {\bf w}_0 - {\bf w}^*\|^2
 \end{align}
 where $\mu$ is strong convexity constant and $L$ is the smoothness parameter \footnote{There is summarized proof for the above result in  \href{https://blogs.princeton.edu/imabandit/2014/03/06/nesterovs-accelerated-gradient-descent-for-smooth-and-strongly-convex-optimization/}{this blog}.)}. The key element in the Nesterov's modification is extrapolation. To observe that, we swap updates as
 \begin{align}
 \v_{n} & = {\bf w}_{n} + \beta ({\bf w}_{n}- {\bf w}_{n-1}) \\
 {\bf w}_{n+1} & = \v_n - \gamma \nabla f(\v_n)
 \end{align}
 Then we combine the above two sequences in one:
 \begin{align}
     {\bf w}_{n+1}= {\bf w}_{n+1/2} - \gamma \nabla f({\bf w}_{n+1/2}), \quad  {\bf w}_{n+1/2} := {\bf w}_{n} + \beta ({\bf w}_{n}- {\bf w}_{n-1})
 \end{align}
 As we see, the gradient is evaluated at ${\bf w}_{n+1/2}$ obtained by the extrapolation ${\bf w}_{n+1/2}$.

\paragraph{Stochastic acceleration schemes}
Replacing the gradient with a stochastic approximation of gradient in HB and Nesterov's acceleration method leads to stochastic HB (SHB) and stochastic Nesterov's acceleration (SNA):
\begin{align}
\text{SHB:} & {\bf w}_{n+1} = {\bf w}_n - \gamma \nabla f_{\bf z}({\bf w}_n) + \beta ({\bf w}_{n} - {\bf w}_{n-1}) \\
\text{SNA:} & {\bf w}_{n+1} = {\bf w}_{n+1/2} - \gamma \nabla f_{\bf z}({\bf w}_{n+1/2}), \quad {\bf w}_{n+1/2}:={\bf w}_{n} + \beta ({\bf w}_{n}- {\bf w}_{n-1})
\end{align}
What is the difference between SGD and the above accelerated stochastic methods? How the noise of stochastic approximation interacts with the momentum term? Can we robustify these method against the noise of stochastic approximation? These are questions that we want to investigate here.
\paragraph{Gaussian input}
 We repeatedly use the following result on Gaussian inputs.
\begin{lemma} \label{lemma:gaussianity_consequence}
Suppose that ${\bf x}$ is drawn from a zero-mean multivariate normal distribution, i.e. ${{\bf x} \sim {\mathcal N}(0,\cov)}$, then the following holds for every matrix ${\bf M}$:
\begin{align} \label{eq:gaussianity_consequence}
    {\mathbf E} \left[ {\bf x} {\bf x}^\top {\bf M} {\bf x} {\bf x}^\top \right] = 2 \cov {\bf M} \cov + \langle {\bf M}, \cov \rangle \cov
\end{align}
where $ \langle {\bf M}, \cov \rangle$ is the inner product of vectorized ${\bf M}$ and $\cov$.
\end{lemma}
\begin{proof}
Let ${\bf A}$ be the desired matrix as
\begin{align} \label{eq:a_matrix}
{\bf A} := {\mathbf E} \left[ {\bf x} {\bf x}^\top  {\bf M} {\bf x} {\bf x}^\top  \right].
\end{align}
Element ${\bf A}_{ij}$ can be computed as
\begin{align}
    {\bf A}_{ij} & = {\mathbf E} \left[ {\bf x}_i {\bf x}_j \sum_{k,l}  {\bf x}_l {\bf M}_{lk} {\bf x}_{k}  \right] \\
    & = \sum_{lk} {\mathbf E} \left[ {\bf x}_i {\bf x}_j {\bf x}_l {\bf x}_k \right]  {\bf M}_{lk}.
\end{align}
Since ${\bf x}$ is Gaussian, we can compute its 4th order statistics using the 2nd order statistics:
\begin{align}
    {\bf A}_{ij} = \cov_{ij} \left( \sum_{kl} \cov_{kl} {\bf M}_{kl} \right) + 2 \sum_{lk} \cov_{ik}  {\bf M}_{kl} \cov_{lj}
\end{align}
Writing the above equation in the matrix equivalent form concludes the proof:
\begin{align} \label{eq:E_y_mn_y}
    {\bf A} := \langle \cov, {\bf M} \rangle \cov + 2 \cov   {\bf M} \cov
\end{align}

\end{proof}

\section{Markov chain and stochastic optimization}
\subsection{SGD and Markov chains}
\paragraph{SGD and Markov chains} Recall SGD iterates on objective $f$ as
\begin{align}
    {\bf w}_{n} = {\bf w}_n - \gamma \nabla f_{\bf z}({\bf w}), \quad \nabla f_{\bf z}({\bf w}) = {\bf x}\x^\top {\bf w} - {\bf y} {\bf x}
\end{align}
Sequence $\{{\bf w}_k\}_{k=1}^n$ makes a homogeneous Markov chain on ${\mathbb{R}}^d$ convergent to a unique stationary distribution. Bridging SGD with Markov chain allows us not only understand the behaviour of SGD better but also it is useful in designing new optimization techniques~\cite{dieuleveut2017bridging}.  For example, stochastic averaging of iterates can be very well-analyzed by tools from Ergodicity theory (See \cite{babichev2018constant} for an overview). Here, we particularly focus on the convergence of constant stepsize SGD to a stationary distribution. Later, we will show the applications of the convergence of SGD to its stationary distribution.
\paragraph{Convergence to the stationary distribution for SGD: A coupling approach}
Here, we prove the convergence of SGD to an stationary distribution --independent from its initialization scheme. Our prove is based on the result of ~\cite{dieuleveut2017bridging}.
Let's consider the space of probability measures on ${\mathbb{R}}^d$ endowed with Wasserstein 2 metric $W$ defined as
\begin{align}
    (W_2(\mu,\nu))^2 = \inf_{(\v,{\bf w}) \in L(\mu,\nu)} {\mathbf E} \left[ \| \v - {\bf w} \|^2 \right]
\end{align}
We know that the space of probability measure with Wasserstein 2 metric is complete. Consider two different sequences of SGD starting at random, namely
\begin{align}
    {\bf w}_{n+1}^{(1)} &= {\bf w}_{n}^{(1)} - \gamma \nabla f_{{\bf z}_n}({\bf w}_{n}^{(1)}), \quad  {\bf w}_0^{(1)} \sim \mu_0 \\
    {\bf w}_{n+1}^{(2)} &= {\bf w}_{n}^{(2)} - \gamma \nabla f_{{\bf z}_n}({\bf w}_{n}^{(2)}), \quad  {\bf w}_0^{(2)} \sim  \nu_0
\end{align}
Remarkably, we the used the same ${\bf z}_n$ in both of updates. Let consider the compact notation $\Delta {\bf w}_{n} := {\bf w}_n^{(1)} - {\bf w}_n^{(2)}$.
This notation and SGD updates yields
\begin{align}
\Delta {\bf w}_{n+1} = \left( {\mathbf I} - \gamma {\bf x}_n {\bf x}^\top_n \right) \Delta {\bf w}_{n} = {\bf B}_n \Delta {\bf w}_{n}
\end{align}
where ${\bf B}_n$ is the product of random matrices in Eq.~\eqref{eq:SGD_iterates}. According to Theorem 4.1. of \cite{villani2008optimal}, there are a couple of random variables ${\bf w}_0$ and $\v_0$, independent from $\{{\bf z}_i\}_{i=1}^n$s, such that
\begin{align} \label{eq:w_match_0}
    (W_2(\nu_0,\mu_0))^2 = {\mathbf E} \| {\bf w}_0 - \v_0 \|^2
\end{align}
Let $\mu_n$ and $\nu_n$ denote the probability laws of random variables ${\bf w}_{n}^{(1)}$ and ${\bf w}_{n}^{(2)}$ (starting from ${\bf w}_0$ and $\v_0$). Then according to the definition of the Wasserstein distance, we have
\begin{align}
    (W_2(\mu_n,\nu_n))^2  & \leq {\mathbf E} \left[ \| {\bf w}_n^{(1)} - {\bf w}_n^{(2)} \|^2  \right]  \\
    & = {\mathbf E} \left[ \| \Delta {\bf w}_n \|^2 \right] \\
    & = (\Delta {\bf w}_0)^\top {\mathbf E} \left[ {\bf B}_n^\top {\bf B}_n \right] (\Delta {\bf w}_0) \\
    & \leq  {\mathbf E} \| \Delta {\bf w}_0 \|^2 \| {\mathbf E} \left[ {\bf B}_n^\top {\bf B}_n \right] \|  \\
    & \stackrel{\eqref{eq:w_match_0}}{=} (W_2(\mu_0,\nu_0))^2 \| {\mathbf E} \left[ {\bf B}_n^\top {\bf B}_n \right] \| \\
    & \leq (W_2(\mu_0,\nu_0))^2 \text{Tr}( {\mathbf E} \left[ {\bf B}_n^\top {\bf B}_n \right]) \\
    & \stackrel{\eqref{eq:eigenvalues_sgd}}{\leq } (W_2(\mu_0,\nu_0))^2 \| {\bf C}^n {\bf 1} \|, \quad \text{Lemma~\ref{lemma:sgd}}\\
    & \leq (W_2(\mu_0,\nu_0))^2 \| {\bf C} \|^n d
\end{align}
A reasonable choice of the stepsize gaurantees $\| {\bf C} \| < 1$ holds in the above inequality concludes an exponential convergence rate in $W_2$. This result is sufficient to conclude that $\mu_n$ is Cauchy. The completeness of the space concludes that $\mu_n$ converges to a unique stationary distribution~\footnote{This result can be extended to an arbitrary strongly convex function under weak assumptions on stochastic gradient (see Proposition 2. of \cite{dieuleveut2017bridging}).}. Therefore, the convergence analysis of the realizable case in Lemma~\ref{lemma:sgd} is equivalent to establishing a rate on the convergence of SGD to its stationary distribution. We can extend this result to stochastic accelerated methods.
\subsection{Acceleration of discrete-state Markov chains}
Consider a discrete state Markov chain with $d$ statates. Transition matrix

\section{Acceleration in Realizable settings}
We turn our focus into realizable least-squares where there is a ${\bf w}^*$ generating the response as $y = {\bf x}^\top {\bf w}^*$.  This setting is recently motivated by over-parameterization models of neural networks \cite{schmidt2018fast,jain2018accelerating}. We will further motivate this setting from Markov chain perspective.

\subsection{Extra extrapolation}
Throughout the rest of the manuscript, we assume that the input is drawn from a zero-mean multivariate normal distribution, i.e. ${\bf x} \sim {\mathcal N}(0,\cov)$. Let's define extra-extrapolated stochastic gradient descent (EESGD) as
\begin{align}
    {\bf w}_{n+1} = {\bf w}_n + \beta ({\bf w}_n - {\bf w}_{n-1}) - \gamma \nabla f_{\bf z}({\bf w}_n + \alpha({\bf w}_{n}-{\bf w}_{n-1})
\end{align}
where $\beta<1$ but $\alpha$ can be greater than 1. We observe that  stochastic HB and stochastic Nesterov's acceleration  requires a small momentum constant $\beta$ to avoid divergence caused by stochasticity of the gradient estimate. A large $\alpha$ allows us to use relatively high $\beta$ in that reach an acceleration. Results in Figure~\ref{fig:extrapolation_result} compare the convergence of SHB, SNA and the proposed extra-extrapolated stochastic gradient descent. Note that we search over $\gamma$ and $\beta$ for this result and it is not clear what is the best choice for stochastic optimization. A clear comparison between these method requires a theoretical justification on the convergence rate obtained by optimal parameters for each method. Indeed, we need to find the optimal $\alpha$, $\beta$, and $\gamma$ to achieve the best rate.
\begin{figure} \label{fig:extrapolation_result}
    \centering
    \includegraphics[width=0.8\textwidth]{images/extra_extrapolation.pdf}
    \caption{Extra-extrapolation is useful in a stochastic setting. In this experiement, the horizontal axis is $\log(\| {\bf w} - {\bf w}^*\|^2)$. The input is ${\bf x} \sim N(0,{\textbf{diag}}(\boldsymbol{\sigma}))$ where $\boldsymbol{\sigma}_i = i^{-4}$ for $i=1, \dots, 16$.}
    \label{fig:my_label}
\end{figure}
\subsection{Warm-up: convergence analysis of SGD}
To warm-up, we prove the convergence rate of SGD on Gaussian inputs. This result give us some intuitions about the noise of SGD and how it can effect the convergence in a realizable setting. Without loss of generality, we assume that ${\bf w}^* = 0$, hence
\begin{align}
    \nabla f_{\bf x}({\bf w}) = {\bf x} {\bf x}^\top {\bf w}.
\end{align}
Using the stochastic gradient formulation, we derive SGD updates by product of random matrices
\begin{align} \label{eq:SGD_iterates}
    {\bf w}_{n} = ({\mathbf I} - \gamma {\bf x}_n {\bf x}^\top_n) {\bf w}_{n-1} = {\bf B}_n {\bf w}_0, \quad {\bf B}_n :=\prod_{k=1}^{n} ({\mathbf I} - \gamma {\bf x}_k {\bf x}^\top_k)
\end{align}
Our goal is to show that ${\bf w}_n$ converges to zero, i.e. the following holds
\begin{align}
    \lim_{n\to \infty}  \| {\bf w}_n \|^2  = 0
\end{align}

\paragraph{Why we can assume that ${\bf w}^* =0$?}
Suppose that ${\bf w}^* \neq 0$, then stochastic gradient can be written as
\begin{align}
     \nabla f_{\bf z}({\bf w}) & = {\bf x}\x^\top {\bf w} - y {\bf x} \\
     & \stackrel{(\eta)}{=} {\bf x} {\bf x}^\top {\bf w} - {\bf x} {\bf x}^\top {\bf w}^* \\
     & = {\bf x} {\bf x}^\top ({\bf w} - {\bf w}^*)
\end{align}
 In $\eta$, we do exploit the realizability as $y = {\bf x}^\top {\bf w}^*$. Now, we can write SGD updates as
 \begin{align}
     {\bf w}_{n+1} - {\bf w}^* = (I - \gamma {\bf x}_n{\bf x}^\top_n)({\bf w}_n - {\bf w}^*)
 \end{align}
 Therefore, we can replace ${\bf w}_n$ (and ${\bf w}_{n-1}$) in Eq.~\eqref{eq:SGD_iterates} by ${\bf w}_n - {\bf w}^*$ (and ${\bf w}_{n-1} - {\bf w}^*$) and all the rest of analysis still holds.

We can prove the above result in two ways: (i) using existing result on product of random matrices or (ii) compute $E\left[ {\bf B}_n^\top {\bf B}_n\right] $ and prove its eigenvalues tends to zero for large $n$. Product of random matrices studies the asymptotic grow rate of $\| {\bf B}_n \|$ by showing that
\begin{align}
    \lim_{n\to \infty} \frac 1n \log(\|{\bf B}_n \|) = \lambda_1
\end{align}
holds with probability one. $\lambda_1$ is the top Lyapunouv exponent of the product of the random matrices (You can find a very short proof of the above result in Example 7.5.1. of \cite{durrett2010probability} which is based on subadditive Ergodic theorem. \cite{bougerol2012products} is an introductory lecture on product of random matrices.).
The later approach (ii) allows us to reach a non-asymptotic result (for finite $n$), while the existing result on product of random matrices all are asymptotic. To show the desired norm decay, we need to bound the top eigenvalue of ${\bf B}_n^\top {\bf B}_n$.
\begin{align}
    {\mathbf E} \| {\bf w}_n \|^2  = {\bf w}_0^\top {\mathbf E} \left[ {\bf B}_n^\top {\bf B}_n \right] {\bf w}_0
\end{align}
To compute the above expectation, we first reverse the indexing, namely we assume that ${\bf B}_n = ({\mathbf I}-\gamma {\bf x}_1 {\bf x}_1^\top) \dots ({\mathbf I} - \gamma {\bf x}_n {\bf x}_n^\top)$. Then, we make the following recurrence
\begin{align}
    {\mathbf E} \left[ {\bf B}_{n+1}^\top {\bf B}_{n+1} \right] &= {\mathbf E} \left[ ({\mathbf I} - \gamma {\bf x} {\bf x}^\top) {\bf B}_n^\top {\bf B}_n ({\mathbf I} - \gamma {\bf x} {\bf x}^\top) \right] \\
    &= {\mathbf E} \left[ {\bf B}_n {\bf B}_n^\top \right]  - \gamma \cov {\mathbf E} \left[ {\bf B}_n^\top {\bf B}_n \right] - \gamma {\mathbf E} \left[ {\bf B}_n^\top {\bf B}_n \right] \cov + \gamma^2 {\mathbf E} \left[ {\bf x} {\bf x}^\top {\bf B}_n^\top {\bf B}_n {\bf x} {\bf x}^\top \right]\nonumber \\
    & \stackrel{\eqref{eq:gaussianity_consequence}}{=} ({\mathbf I} - \gamma \cov) {\mathbf E} \left[ {\bf B}_n {\bf B}_n^\top \right] ({\mathbf I} - \gamma \cov) + \gamma^2 \cov {\mathbf E} \left[ {\bf B}_n^\top {\bf B}_n \right] + \langle \cov, {\mathbf E} \left[ {\bf B}_n^\top {\bf B}_n \right] \rangle \cov \label{eq:E_BB_recursive}
\end{align}
The above recurrence allows us to track the expectation ${\mathbf E}\left[ {\bf B}_n^\top {\bf B}_n \right]$ in a closed-form. The next lemma tracks eigenvalues and eigenvectors of  this matrix separately \footnote{We established a similar result for PCA problem during my visit at INRIA, but we later found a better way to solve our question. The analysis of the Gaussian case gave us very helpful intitutions there. }.
\begin{lemma} \label{lemma:sgd}
    Consider the following eigenvalue decomposition for the covariance matrix of the input $\cov = {\bf U} {\textbf{diag}}(\boldsymbol{\sigma}){\bf U}^\top$. If ${\bf x} \sim {\mathcal N}(0,\cov)$ then
    \begin{align}
        {\mathbf E} \left[ {\bf B}_n^\top {\bf B}_n \right] = {\bf U} {\textbf{diag}}(\boldsymbol{ \lambda}_n) {\bf U}^\top
    \end{align}
    where
    \begin{align} \label{eq:eigenvalues_sgd}
        \boldsymbol{ \lambda}_n = {\bf C}^n {\bf 1}, \quad {\bf C} := ({\mathbf I} - \gamma {\textbf{diag}}(\boldsymbol{\sigma}))^2 + \gamma^2 {\textbf{diag}}(\boldsymbol{\sigma})^2 + \gamma^2 \boldsymbol{\sigma} \boldsymbol{\sigma}^\top
    \end{align}
\end{lemma}
\begin{proof}
Our proof is based on induction: supposing the above result holds for $n$, we prove that it holds for $n+1$.
\begin{multline}
    {\mathbf E} \left[ {\bf B}_{n+1} {\bf B}_{n+1}^\top \right] \stackrel{\eqref{eq:E_BB_recursive}}{=}  {\bf U}  \left(  ({\mathbf I}-\gamma {\textbf{diag}}(\boldsymbol{\sigma})) {\textbf{diag}}(\boldsymbol{ \lambda}_n) ({\mathbf I} - \gamma {\textbf{diag}}(\boldsymbol{\sigma})) \right) {\bf U}^\top
    \\ +  \gamma^2 {\bf U} \left( {\textbf{diag}}(\boldsymbol{\sigma}) {\textbf{diag}}(\boldsymbol{ \lambda}_n) {\textbf{diag}}(\boldsymbol{\sigma}) + \langle\boldsymbol{\sigma},\boldsymbol{ \lambda}_n\rangle {\textbf{diag}}(\boldsymbol{\sigma}) \right) {\bf U}^\top
\end{multline}
From the above equation, we can prove that $\boldsymbol{ \lambda}_{n+1}$ is obtained by the following affine transformation of $\boldsymbol{ \lambda}_n$:
\begin{align}
     \boldsymbol{ \lambda}_{n+1} = {\bf C} \boldsymbol{ \lambda}_{n}
\end{align}
where ${\bf C}$ is the matrix in Eq.~\eqref{eq:eigenvalues_sgd}. This result conclude our induction.
\end{proof}
According to the result of last lemma, SGD with stepsize $0.5/(L+\text{Tr}[\cov])$ linearly converge to ${\bf w}^* =0$ with rate $1-0.5\mu/(L+\text{Tr}[\cov])$. The rate of convergence is slightly worst than the convergence rate of gradient descent, which is $1-\mu/L$.

\subsection{Convergence analysis for accelerated stochastic optimization}
Recall our general stochastic accelerated method as
\begin{align}
    {\bf w}_{n+1} = {\bf w}_n + \beta ({\bf w}_n - {\bf w}_{n-1}) - \gamma \nabla f_{\bf z}({\bf w}_n + \alpha({\bf w}_{n}-{\bf w}_{n-1}))
\end{align}
Setting $\alpha =0$ in the above equation recovers SHB method and $\alpha = \beta$ reduces to Stochastic Nestrov's acceleration method. Our question is: is there any choice for parameters $\alpha$, $\beta$ and $\gamma$ such that the above method improves upon the established convergence rate of SGD in the last section. To address this question, we first write the update in the following linear recurrence:
\begin{align}
    \begin{bmatrix}
    {\bf w}_{n+1} \\
    {\bf w}_{n}
    \end{bmatrix} = {\bf A}_n \begin{bmatrix}
    {\bf w}_{n} \\
    {\bf w}_{n-1}
    \end{bmatrix}, \quad {\bf A}_n := \begin{bmatrix}
    (1+\beta) {\mathbf I}- (1+\alpha) \gamma {\bf x}_n{\bf x}^\top_n & \alpha \gamma {\bf x}_n {\bf x}^\top_n - \beta {\mathbf I} \\
    {\mathbf I} & 0
    \end{bmatrix}
\end{align}
We write matrix ${\bf A}_{n+1}$ as
 \begin{align} \label{eq:An_decomposition}
     {\bf A}_{n+1} = {\bf A} -\gamma
     \epsilon, \quad \epsilon := \begin{bmatrix}
     (1+\alpha) \Delta & -\alpha\Delta \\
     0 & 0
     \end{bmatrix}, \Delta :=(\cov - {\bf x} {\bf x}^\top)
 \end{align}
 where ${\mathbf E} \left[ \epsilon \right] = 0$ and matrix ${\bf A}$ is
 \begin{align}
     {\bf A} = \begin{bmatrix}
     {\bf A}_1 & {\bf A}_2 \\
     {\mathbf I} & 0
     \end{bmatrix}, {\bf A}_1 := (1+\beta){\mathbf I} - \gamma (1+\alpha) \cov, {\bf A}_2 := \alpha \gamma \cov - \beta {\mathbf I}
 \end{align}
 The above decomposition allows us to analyze the covariance term induced by the noise separately.
Similar to our analysis for SGD, we aim at tracking matrix $ {\bf M}_n := {\mathbf E} \left[{\bf A}_n^\top \dots {\bf A}_1^\top {\bf A}_1 \dots {\bf A}_n \right] $ in a closed form. The next lemma provides a recurrence that later will be used to track the desired matrix.
\begin{lemma}
Suppose ${\bf x} \sim {\mathcal N}(0,\cov)$. Consider the following block structure:
\begin{align}
     {\mathbf E} \left[ {\bf M}_n \right] = \begin{bmatrix}
      {\bf M}_1 & {\bf M}_2 \\
      {\bf M}_2 & {\bf M}_3
     \end{bmatrix}.
\end{align}
Then matrix ${\mathbf E}\left[ {\bf M}_{n+1} \right]$ can be decomposed as
\begin{align} \label{eq:E_nplus_expansion_result}
    {\mathbf E} \left[ {\bf M}_{n+1} \right] = {\bf A}^\top  {\mathbf E} \left[{\bf M}_n \right]{\bf A} + \gamma^2 \underbrace{\begin{bmatrix} (1+\alpha)^2 {\bf K} &  -\alpha(1+\alpha) {\bf K} \\
 -\alpha (1+\alpha) {\bf K}  & \alpha^2 {\bf K}
 \end{bmatrix}}_{\text{due to the noise}}
\end{align}
where  \begin{align}
    {\bf K} := \cov {\bf M}_1 \cov + \langle \cov,{\bf M}_1 \rangle \cov.
\end{align}
\end{lemma}
\begin{proof}
 Using the decomposition of matrix ${\bf A}_n$ in Eq.~\eqref{eq:An_decomposition}, we expand the expectation as
 \begin{align}\label{eq:M_nplus}
     {\mathbf E} \left[ {\bf M}_{n+1} \right] & = {\mathbf E} \left[{\bf A}_{n+1}^\top  {\bf M}_n {\bf A}_{n+1}  \right] \\
     & = {\bf A}^\top  {\mathbf E} \left[{\bf M}_n \right]{\bf A} + \gamma^2  {\mathbf E} \left[ \epsilon^\top  {\bf M}_n \epsilon \right]
 \end{align}

We compute 2nd term (induced by the noise)
\begin{align} \label{eq:variance_noise}
    {\mathbf E} \left[ \epsilon^\top  {\bf M}_n \epsilon \right]  = \begin{bmatrix}
    (1+\alpha)^2 \Delta {\bf M}_1 \Delta & -\alpha (1+\alpha) \Delta{\bf M}_1 \Delta \\
    -\alpha (1+\alpha) \Delta{\bf M}_1 \Delta & \alpha^2 \Delta {\bf M}_1 \Delta
    \end{bmatrix}
\end{align}
Using lemma~\ref{lemma:gaussianity_consequence} on Gaussian inputs, we have
\begin{align}
    {\mathbf E} \left[ \Delta {\bf M}_1 \Delta \right] =  \cov {\bf M}_1 \cov + \langle {\bf M}_1,\cov \rangle \cov
\end{align}
Replacing this result into Eq.~\eqref{eq:M_nplus} concludes the proof.
\end{proof}
The first term in Eq.~\eqref{eq:E_nplus_expansion_result} is the update with exact gradient and the 2nd term is caused by the noise of the stochastic gradient.
Using the above result, next corollary proves that ${\mathbf E}\left[ {\bf M}_n \right]$ can be computed through a linear fixed point iteration.
\begin{corollary}
Suppose that ${\bf x} \sim N(0,\cov) $ where $\cov = {\bf U} {\textbf{diag}}(\boldsymbol{\sigma}) {\bf U}^\top$; then
\begin{align}
     {\mathbf E} \left[ {\bf M}_n \right] = \begin{bmatrix}
      {\bf M}_1 & {\bf M}_2 \\
      {\bf M}_2 & {\bf M}_3
     \end{bmatrix}, \quad {\bf M}_{i} = {\bf U} {\textbf{diag}}(\boldsymbol{ \lambda}_n^{(i)}) {\bf U}^\top
\end{align}
where
\begin{align}
    \begin{bmatrix}
    \boldsymbol{ \lambda}_n^{(1)}
     \\
    \boldsymbol{ \lambda}_n^{(2)}
    \\
    \boldsymbol{ \lambda}_n^{(3)}
     \end{bmatrix} =
     {\bf C}^n
    {\bf 1}
\end{align}
and matrix ${\bf C}$ is
\begin{align}
    {\bf C} & = \begin{bmatrix} {\bf D}_1^2 + (1+\alpha)^2 {\bf K}_\epsilon & 2 {\bf D}_1 & {\mathbf I} \\
     {\bf D}_1 {\bf D}_2 - \alpha(1+\alpha) {\bf K}_\epsilon& {\bf D}_2 & 0 \\
     {\bf D}_2^2 +\alpha^2 {\bf K}_\epsilon & 0 & 0
    \end{bmatrix}  \label{eq:C_matrix}\\
    {\bf D}_1 &= (1+\beta) {\mathbf I} - \gamma (1+\alpha) {\textbf{diag}}(\boldsymbol{\sigma}),\\
    {\bf D}_2 &=  \alpha \gamma {\textbf{diag}}(\boldsymbol{\sigma}) - \beta {\mathbf I} \\
    {\bf K}_\epsilon &= \gamma^2 \left( {\textbf{diag}}(\boldsymbol{\sigma})^2 +\boldsymbol{\sigma} \boldsymbol{\sigma}^\top \right)
\end{align}
\end{corollary}
\begin{proof}
  We prove the result by induction.
  Recall matrix ${\bf A}$ with the following block structure:
\begin{align}
     {\bf A} = \begin{bmatrix}
     {\bf A}_1 & {\bf A}_2 \\
     {\mathbf I} & 0
     \end{bmatrix}, {\bf A}_1 := (1+\beta){\mathbf I} - \gamma (1+\alpha) \cov, {\bf A}_2 := \alpha \gamma \cov - \beta {\mathbf I}
 \end{align}
Indeed, ${\bf D}_i$ are the diagonal matrix of eigenvalues of matrix ${\bf A}_i$. We decompose ${\mathbf E}\left[ {\bf M}_n \right]$ as
 \begin{align}
      {\mathbf E} \left[ {\bf M}_n \right] = \begin{bmatrix}
      {\bf M}_1 & {\bf M}_2 \\
      {\bf M}_2 & {\bf M}_3
     \end{bmatrix}, {\bf M}_i = {\bf U} {\textbf{diag}}(\boldsymbol{ \lambda}_n^{(i)}) {\bf U}^\top
 \end{align}
  Using the above expression,
  we expand the first term in Eq.~\eqref{eq:E_nplus_expansion_result}
\begin{align}
    {\bf A}^\top {\mathbf E} \left[ {\bf M}_n \right] {\bf A} & = \begin{bmatrix}
    {\bf A}_1 {\bf M}_1 {\bf A}_1 + 2{\bf M}_2 {\bf A}_1 + {\bf M}_3 & {\bf A}_1 {\bf M}_1 {\bf A}_2 + {\bf M}_2 {\bf A}_2 \\
    {\bf A}_2 {\bf M}_1 {\bf A}_1 + {\bf M}_2 {\bf A}_2 & {\bf A}_2 {\bf M}_1 {\bf A}_2
    \end{bmatrix} \\
    & =
    \begin{bmatrix}
     {\bf B}_1 & {\bf B}_2 \\
     {\bf B}_2 & {\bf B}_3
    \end{bmatrix}
\end{align}
where
\begin{align}
    {\bf B}_1 &= {\bf U} \left( {\bf D}_1^2 {\textbf{diag}}(\boldsymbol{ \lambda}_n^{(1)}) + 2 {\bf D}_1{\textbf{diag}}(\boldsymbol{ \lambda}_n^{(2)}) + {\textbf{diag}}(\boldsymbol{ \lambda}_n^{(3)})\right) {\bf U}^\top \\
    {\bf B}_2 &= {\bf U} \left( {\bf D}_1 {\bf D}_2 {\textbf{diag}}(\boldsymbol{ \lambda}_n^{(1)}) + {\bf D}_2 {\textbf{diag}}(\boldsymbol{ \lambda}_n^{(2)}) \right) {\bf U}^\top \\
    {\bf B}_3 & = {\bf U} \left(  {\bf D}_2^2 {\textbf{diag}}(\boldsymbol{ \lambda}_n^{(1)}) \right) {\bf U}^\top
\end{align}
Using eigendecomposition of $\cov$ and ${\bf M}_1$, the matrix ${\bf K}$ in the noise-induced term can be written as
\begin{align}
     {\bf K}  = {\bf U} \left(  {\textbf{diag}}(\boldsymbol{\sigma})^2 {\textbf{diag}}(\boldsymbol{ \lambda}_n^{(1)}) + \langle \boldsymbol{ \lambda}_n^{(1)}, \boldsymbol{\sigma} \rangle {\textbf{diag}}(\boldsymbol{\sigma}) \right) {\bf U}^\top
\end{align}
Putting all together, we complete the proof:
\begin{align}
    {\mathbf E} \left[ {\bf M}_{n+1} \right] = \begin{bmatrix}
    {\bf M}_{1}^+ & {\bf M}_2^+ \\
    {\bf M}_2^+ &  {\bf M}_3^+
    \end{bmatrix}, \quad  {\bf M}_i^+ = {\bf U} {\textbf{diag}}(\boldsymbol{ \lambda}_{n+1}^{(i)}) {\bf U}^\top
\end{align}
where
\begin{align}
    \begin{bmatrix}
    \boldsymbol{ \lambda}_{n+1}^{(1)}
     \\
    \boldsymbol{ \lambda}_{n+1}^{(2)}
    \\
    \boldsymbol{ \lambda}_{n+1}^{(3)}
     \end{bmatrix} =
     {\bf C}
     \begin{bmatrix}
    \boldsymbol{ \lambda}_n^{(1)}
     \\
    \boldsymbol{ \lambda}_n^{(2)}
    \\
    \boldsymbol{ \lambda}_n^{(3)}
     \end{bmatrix}
\end{align}
where matrix ${\bf C}$ is those of Eq.~\eqref{eq:C_matrix}.
The above results concludes our induction.
\end{proof}
Using the result of last corollary, we only need to analyze eigenvalues of matrix ${\bf C}$ to prove the convergence of stochastic accelerated methods. We highlighted the noise effect with the subscript $\epsilon$. We clearly see the effect of changing extrapolation constant $\alpha$ in matrix ${\bf C}$. If $\alpha=0$, then the noise will be applied just to $\boldsymbol{ \lambda}_{n+1}^{(1)}$.  For a larger $\alpha$ the noise will be distributed through all eigenvalues $\boldsymbol{ \lambda}_{n+1}^{(i)}$. The complex eigenvalues of matrix ${\bf C}$ can be analyze using a similar proof technique that is used in Heavy ball analysis. Here is our plan for computing eigenvalues of ${\bf C}$.
\begin{itemize}
    \item We first establish a bound on eigenvalues of matrix ${\bf K}_\epsilon$. Note that this matrix is obtained by a rank one modification of a diagonal matrix, hence we can easily bound its spectrum.
    \item Then we use the rearrangement of row and columns of ${\bf C}$ to make it block diagonal. In this way, we need to analysis eigenvalues of $3\times 3$ matrices to derive eigenvalues of ${\bf C}$.
\end{itemize}

\subsection{Spectral radius analysis: 2D case}
In the rest of the document, we establish a bound on the spectral radius of ${\bf C}$. Recall matrix ${\bf C}$ formulation:
\begin{align}
    {\bf C} : = \begin{bmatrix} {\bf C}_1 & 2 {\bf D}_1 & {\mathbf I} \\
     {\bf C}_2 & {\bf D}_2 & 0 \\
     {\bf C}_3 & 0 & 0
    \end{bmatrix}
\end{align}
where
\begin{align}
    {\bf C}_1 & = {\bf D}_1^2 + (1+\alpha)^2 {\bf K}_\epsilon, \quad {\bf D}_1 := (1+\beta) {\mathbf I} - \gamma (1+\alpha) {\textbf{diag}}(\boldsymbol{\sigma})\\
    {\bf C}_2 & = {\bf D}_1 {\bf D}_2 - \alpha(1+\alpha) {\bf K}_\epsilon, \quad {\bf D}_2 =  \alpha \gamma {\textbf{diag}}(\boldsymbol{\sigma}) - \beta  {\bf K}_\epsilon\\ {\bf C}_3 & = {\bf D}_2^2 +\alpha^2 {\bf K}_\epsilon, \quad {\bf K}_\epsilon := \gamma^2 \left( {\textbf{diag}}(\boldsymbol{\sigma})^2 +\boldsymbol{\sigma} \boldsymbol{\sigma}^\top \right)
\end{align}

To warm-up, we start our analysis with a 2D Gaussian input such that
\begin{align}
    {\mathbf E} \left[ {\bf x} {\bf x}^\top \right] = {\textbf{diag}}( \boldsymbol{\sigma}), \quad \boldsymbol{\sigma} :=[L, \mu]
\end{align}
where $L$ is the maximum eigenvalue of input covariance and $\mu$ is the minimum eigenvalue of the input. Our focus is on cases where $L$ is significantly larger than $\mu$, i.e. $L \gg \mu$. In this case, we can decompose ${\bf K}_\epsilon$ (the induced matrix by noise) as sum of diagonal matrices and a non-diagonal matrix:
\begin{align}
    {\bf K}_\epsilon := \gamma^2 \left( {\textbf{diag}}(\boldsymbol{\sigma})^2 + {\textbf{diag}}\left([L^2,\mu^2]\right) + \xi\right), \quad \xi = \begin{bmatrix}
    0 &  \mu L \\
    \mu L & 0
    \end{bmatrix}
\end{align}
The non-diagonal matrix $\xi$ has a bounded norm since
\begin{align}
    \| \xi \| \leq \| \xi \|_F  \leq 2 \mu L.
\end{align}
Using the decomposition of ${\bf K}_\epsilon$, we decompose ${\bf C}$ as a perturbation of a block diagonal matrix:
\begin{align}
    {\bf C} = \bar{{\bf C}} + \gamma^2
    \begin{bmatrix}
    (1+\alpha)^2 \xi & 0 & 0 \\
    -\alpha(1+\alpha) \xi & 0 & 0 \\
    \alpha^2 \xi & 0 & 0
    \end{bmatrix}
\end{align}
where matrix $\bar{C}$ is
\begin{align}
    \bar{{\bf C}} := \begin{bmatrix} {\bf D}_1^2 + (1+\alpha)^2 \bar{{\bf K}}   & 2 {\bf D}_1 & {\mathbf I} \\
     {\bf D}_1 {\bf D}_2 - \alpha(1+\alpha) \bar{{\bf K}}& {\bf D}_2 & 0 \\
     {\bf D}_2^2 +\alpha^2 \bar{{\bf K}} & 0 & 0
    \end{bmatrix}, \quad \bar{{\bf K}} := 2 \gamma^2 \left( {\textbf{diag}}(\boldsymbol{\sigma})^2  \right)
\end{align}
Our goal is using pseudospectrum of ${\bf C}$ to bound  $\| {\bf C}^n \|$ (see~\cite{jensen09}).
\paragraph{Background on pseudospectrum}
Pseudospectrum is mainly developed for perturbation analysis of non-hermitian matrix. Let $\sigma({\bf A})$ be the set of complex eigenvalues of the non-symmetric matrix ${\bf A}$; then the spectral radius of ${\bf A}$ is $\rho({\bf A}) = \sup\{ |{\bf z}|\; | \; {\bf z} \in \sigma({\bf A}) \}$. $\epsilon$-pseudospectrum of ${\bf A}$, denoted by  $\sigma_{\epsilon}({\bf A})$, is defined as
\begin{align}
    \sigma_{\epsilon}({\bf A}) := \sigma({\bf A}) \cup \{ {\bf z} \in C \; | \; \| ({\bf A} - {\bf z} {\mathbf I})^{-1} \| \geq 1/\epsilon \}.
\end{align}
Pseudospectrum of ${\bf A}$ is a genealization of spectral radius for $\epsilon$-pseudospectrum as
\begin{align}
    \rho_\epsilon({\bf A}) = \sup\{ |{\bf z}| \; | \; {\bf z} \in \sigma({\bf A}) \}
\end{align}
Three key results on $\rho_\epsilon({\bf A})$ are presented in next three lemmas.
\begin{lemma}[Matrix power and Pseudospectrum (Theorem 9.2\cite{jensen09})]
The following holds for any $\epsilon$ and all $n$:
\begin{align}
    \| {\bf A}^n \| \leq \frac{(\rho_{\epsilon}({\bf A}))^{n+1}}{\epsilon}
\end{align}
\end{lemma}
\begin{lemma}[Robustness of Pseudospectrum (Theorem 5.12.\cite{jensen09})]
For all matrices ${\bf M}$, the following holds
\begin{align}
    \rho_{\epsilon}({\bf A} + {\bf M}) \leq \rho_{\epsilon+\|{\bf M}\|}({\bf A})
\end{align}
\end{lemma}
\begin{lemma}[Bauer–Fike; see Thm. 5.11 of \cite{jensen09}]
Let ${\bf A}$ be a diagonalizable square matrix such that ${\bf A} = {\cal V} \Lambda {\cal V}^{-1}$. Then for $\epsilon>0$, the following holds:
\begin{align}
    \rho_\epsilon({\bf A}) \leq \kappa \epsilon
\end{align}
where $\kappa$ is the condition number of ${\cal V}$.
\end{lemma}
\paragraph{Diagonalization of $\bar{{\bf C}}$}

Since matrix $\bar{{\bf C}}$ has a block diagonal structure, we can permute its rows and columns to make a block diagonal matrix as
\begin{align}
    \Pi \bar{{\bf C}} \Pi^\top = \begin{bmatrix}
     {\bf M}_{1} & 0&  0 & \dots & 0 \\
     0 & {\bf M}_{2} & 0 & \dots & 0 \\
     \vdots & \vdots & \vdots & \vdots & \vdots \\
     0 & \dots & \dots & 0 & {\bf M}_{d}
    \end{bmatrix}
\end{align}
where ${\bf M}_{i}$ is a $3 \times 3$ matrix as
\begin{align}
    {\bf M}_{i} := \begin{bmatrix}
     m_1  & m_2 & 1 \\
    m_3 & m_4 & 0 \\
    m_5 & 0 & 0
    \end{bmatrix}
\end{align}
where
\begin{align}
    m_1 & := D_{1i}^2 + 2 (1+\alpha)^2 \gamma^2 \sigma_i^2  \\
    m_2 & := 2 D_{1i} \\
    m_3 & := D_{1i}D_{2i} - 2 \alpha (1+\alpha)  \gamma^2 \sigma_i^2 \\
    m_4 & := D_{2i} \\
    m_5 & := D_{2i}^2 + 2 \gamma^2 \alpha^2 \sigma_i^2
\end{align}
Given diagonalization ${\bf M}_i = {\bf U}_i \Lambda_i {\bf U}_i^\top $ , then we have
\begin{align}
     \bar{{\bf C}}
    =
    {\cal V}
    \begin{bmatrix}
    \Lambda_{1} & 0&  0 & \dots & 0 \\
     0 & \Lambda_{2} & 0 & \dots & 0 \\
     \vdots & \vdots & \vdots & \vdots & \vdots \\
     0 & \dots & \dots & 0 & \Lambda_{d}
     \end{bmatrix}
    {\cal V}^\top
\end{align}
where matrix ${\cal V}$ is an orthogonal matrix as
\begin{align}
    {\cal V} := \Pi^\top
    \begin{bmatrix}
        {\bf U}_1  & 0 & \dots \\
        0 & {\bf U}_2  & \dots \\
        \vdots & \vdots & \vdots
    \end{bmatrix}.
\end{align}
Since ${\cal V}$ is orthogonal, its condition number is $1$.

\subsection{Generalization}
Here we generalize our result to a structural input distribution as
\begin{align}
    {\bf x} = {\bf U} \v
\end{align}
where ${\bf U}$ is an orthogonal matrix and $\v$ is from a symmetric distribution, namely $\v_i$ is distributed as $-\v_i$. We will repeatly use the result of

\section{Related works}
\cite{dieuleveut2017harder} has proved an improved convergence rate for averaged iterates of stochastic accelerated Nesterov when the covariance matrix is degenerated. In this case, they show that acceleration can improve $1/n$ convergence rate of SGD to $1/n^2$. When the matrix is not degenerated SGD can reach an exponential convergence rate $(1-c\mu/(L+\text{Tr}{\cov}))$ \cite{pillaud2017exponential,dieuleveut2017bridging} (we will prove this result later). \cite{jain2018accelerating} has shown that Nesterov's acceleration can improve the rate of SGD to $(1-\sqrt{\mu/(L}K)$. However, the constant $K$ in the rate can be larger than $\sqrt{\mu/L}$ in worst-cases. \cite{schmidt2018fast} has shown that accelerated stochastic gradient descent can obtain the square root improvement in the convergence rate for learning halfspace (which is a non-convex program). However, this result strongly rely on strong growth condition. The stochastic gradients satisfies strong grow condition if the following holds uniformly in ${\bf w}$:
\begin{align} \label{assum:sgc}
    {\mathbf E}_{\bf z}  \left[ \| \nabla f_{\bf z}({\bf w}) \|^2 \right] \leq \rho \| \nabla f({\bf w})\|^2
\end{align}
Under the above assumption,~\cite{schmidt2018fast} show that a complicated update scheme over  3 sequences of iterates can reach to $1-\sqrt{\mu/(\rho^2L)}$ convergence rate. In the next lemma, we show that the constant $\rho$ can be as larger as $L/\mu$ for Gaussian inputs, hence their rate can is $1-(\mu/L)^{3/2}$ in worst case. This rate is worst than the rate for SGD.
\begin{lemma}
    Consider the regression objective of Eq.~\ref{eq:f} on Gaussian input ${\bf x} \sim {\mathcal N}(0,\cov)$ and suppose that $y=0$. There exists a vector ${\bf w}_0 \in {\mathbb{R}}^d$ such that assumption~\ref{assum:sgc} does not hold for $\rho<L/\mu$.
\end{lemma}
\begin{proof}
Let ${\bf w}_0$ be the eigenvector of $\cov$ associated with the smallest eigenvalue of $\cov$, i.e. ${\cov {\bf w}_0 = \mu {\bf w}_0}$. Then we use the result of lemma~\ref{lemma:gaussianity_consequence} to compute expected norm of stochastic gradients as
\begin{align}
    {\mathbf E} \| \nabla f_{\bf z}({\bf w}_0) \|^2 & = {\bf w}^\top_0 {\mathbf E} \left[ {\bf x} {\bf x}^\top {\bf x} {\bf x}^\top \right] {\bf w}_0 \\
    & \stackrel{\eqref{eq:gaussianity_consequence}}{=} {\bf w}_0^\top \cov^2 {\bf w}_0 + \text{Tr}[\cov]{\bf w}_0^\top \cov {\bf w}_0 \\
    & = \mu^2 + \text{Tr}[\cov] \mu \geq \mu (\mu+L)
\end{align}
The norm of gradient evaluated at ${\bf w}_0$ is
\begin{align}
    \|\nabla f({\bf w}) \|^2 = {\bf w}^\top_0 \cov^2 {\bf w}_0 = \mu^2
\end{align}
To satisfy assumption of Eq.~\eqref{assum:sgc}, we need to choose $\rho$ such that
\begin{align}
    \mu (\mu+L) \leq \rho \mu^2 \longleftarrow \rho \geq 1 + L/\mu.
\end{align}
\end{proof}

\section{Experiments}
We empirically validate the established result in Theorem~\ref{thm:convergence_char}, and Lemmas~\ref{lemma:gaussian_gaurantees} and~\ref{lemma:uniform_gaurantees}. Then, we empirically show that the result of Theorem~\ref{thm:convergence_char} may hold even on real data sets on which the Assumption~\ref{assume:symmetric_inputs} does not necessarily hold.
\paragraph{Experiments on examples.}
Through an experiment, we check whether the accelerated rate established in Lemmas~\ref{lemma:gaussian_gaurantees} and~\ref{lemma:uniform_gaurantees} are achievable. Recall that these results hold for examples~\ref{exam:gaussian_example} and \ref{exam:uniform_example}. Our experimental results, presented in Figure~\ref{fig:accelerated_rates_real_world_theory}, confirm that SAGD-chain enjoys the accelerated mixing rate $\mathcal{O}((1-\sqrt{\mu/L})^n)$ using the choice of parameters in corresponding lemmas.
\begin{figure}[h!]
    \centering
    \begin{tabular}{c}
        \includegraphics[width=0.4\textwidth]{images/accelerated_gaussian.pdf}  \\
        Example~\ref{exam:gaussian_example}  \\ \includegraphics[width=0.4\textwidth]{images/accelerated_uniform.pdf}  \\
        Example~\ref{exam:uniform_example}
    \end{tabular}
    \caption{Mixing of SAGD on examples: horizontal axis shows number of iterations and vertical axis is $\log_{10}{\mathbf E}\| \u_{n}^{(0)} - \u_n^{(1)} \|^2 $ where the expectation is taken over 5 independent runs. Remarkably, $W_2^2(\nu(\u_n^{(0)}),\nu(\u_n^{(1)})) \leq {\mathbf E}\| \u_{n}^{(0)} - \u_n^{(1)} \|^2$. See Lemma~\ref{lemma:gaussian_gaurantees} and ~\ref{lemma:uniform_gaurantees} for details on the parameter choice.  }
    \label{fig:accelerated_rates_examples}
\end{figure}

\paragraph{Validation of the establish mixing rate.}
Next, we compare the theoretical mixing rate established in Theorem~\ref{thm:convergence_char}, with empirical ones under different parameter configurations $\Theta$. Experiments run on  Example~\ref{exam:gaussian_example} and \ref{exam:uniform_example} and the comparison is represented in Table~\ref{tbl:gaussian} and \ref{tbl:symmetric}, respectively. A total iteration of $n=1000$ is employed for both two examples. This experiment confirms empirical rates are generally consistent with established theoretical mixing rates up to constant factors. We note that we estimate the Wasserstein distance with ${\mathbf E} \left[ \| \u_n^{(0)}- \u_n^{(1)} \|^2 \right]$ for the sake of simplicity. The expectation is taken over 10 independent runs.

\begin{table}[h]
    \caption{${\mathbf E} \left[ \| \u_n^{(0)}- \u_n^{(1)} \|^2 \right]$ for Example 8, $\mu=0.05$} \label{tbl:gaussian}
    \begin{center}
        \begin{tabular}{lll}
            \textbf{$\Theta=(\gamma,\beta,\alpha)$}  & Empirical rate &  Theoretical rate \\
            \hline \\
            $(10^{-1},0.95,2)$  & $e^{-0.0605n}$  & $ e^{-0.0568n}$ \\
            $(10^{-1},0.99,2)$  & $e^{-0.0164n}$  & $ e^{-0.0156n}$ \\
            $(10^{-1},0.95,3)$  & $ e^{-0.0628n}$  & $e^{-0.0623n}$ \\
            $(10^{-2},0.95,2)$  & $e^{-0.0260n}$  & $e^{-0.0257n}$ \\
        \end{tabular}
    \end{center}
\end{table}

\begin{table}[h]
    \caption{${\mathbf E} \left[ \| \u_n^{(0)}- \u_n^{(1)} \|^2 \right]$ for Example 11, $\kappa^{-1}=20$} \label{tbl:symmetric}
    \begin{center}
        \begin{tabular}{lll}
            \textbf{$\Theta=(\gamma,\beta,\alpha)$}  & Empirical rate & Theoretical rate \\
            \hline \\
            $(2\times10^{-3},0.95,2)$  & $e^{-0.0626n}$   & $e^{-0.0527n}$ \\
            $(2\times10^{-3},0.99,2)$  & $e^{-0.0188n}$  & $ e^{-0.0118n}$ \\
            $(2\times10^{-3},0.95,3)$  & $e^{-0.0572n}$  & $ e^{-0.0537n}$ \\
            $(4\times10^{-4},0.95,2)$  & $e^{-0.0206n}$  & $e^{-0.0086n}$ \\
        \end{tabular}
    \end{center}
\end{table}

\paragraph{Experiments on real-world data.}
We substantiate our results on two real-world datasets: Boston Housing and California Housing. These datasets are 13 and 8 dimensional respectively. As a reprocessing step, we normalized inputs (this guarantees $L<1$). We further used $\ell_2$ regularization with the penalty factor $10^{-3}$ (this guarantees $\mu>10^{-3}$). The experiments compare contraction rate, established in Theorem~\ref{thm:convergence_char}, on these datasets. Although Assumption~\ref{assume:symmetric_inputs} does not hold for these datasets, we can reuse the choice of parameters in Lemma~\ref{lemma:gaussian_gaurantees} to achieve an accelerated rate.

 \begin{figure}[h!]
    \centering
    \begin{tabular}{c c}
        \includegraphics[width=0.4\textwidth]{images/accelerated_boston_theory_2.pdf}  \\
        Boston Housing  \\
        \includegraphics[width=0.4\textwidth]{images/accelerated_california_theory_2.pdf}  \\
        California Housing
    \end{tabular}
    \caption{Accelerated rates achieved on real-world data for theoretical parameters. The horizontal axis shows the number of stochastic iterations $n$ and the vertical axis is $\log_{10}\| \u_{n}^{(0)} - \u_n^{(1)} \|^2 $. }
    \label{fig:accelerated_rates_real_world_theory}
\end{figure}

\section{A quest for understanding SAGD} \label{sec:intution}

\section{Introduction}\label{sec:intro}
Stochastic variants of gradient based optimization methods have become the \textit{de facto} standard optimization technique for large scale learning problems. Trading off statistical and computational aspects, stochastic approximation methods attempt to obtain high statistical accuracy while keeping the computational per-iteration  costs low~\cite{bottou2008tradeoffs}. The vast empirical success of such methods has motivated a growing body of theoretical studies on stochastic approximation in both convex (e.g.\cite{moulines2011non,schmidt2017minimizing,rakhlin2011making}) and non-convex optimization (e.g. \cite{ge2015escaping,daneshmand2018escaping,zhu2018anisotropic}). Most remarkably, the analysis in~\cite{moulines2011non} establishes fast, non-asymptotic convergence rates for Stochastic Gradient Descent (SGD). Despite this growing understanding of SGD, the inner workings of stochastic \textit{accelerated} first-order methods are still not very well-understood.  Inspired by the success of accelerated schemes (such as stochastic momentum and Adam~\cite{kingma2014adam}) in optimization of deep neural networks, interesting recent works are starting to improve the current theoretical understanding of this empirical observation \cite{jain2018accelerating,schmidt2018fast,dieuleveut2017harder}.

We here contribute to this line of research, starting from the simplest possible setting relevant for machine learning: ordinary least-squares regression. Namely we consider following problem set-up throughout this paper:
\begin{align}\label{eq:main_leastsquares}
    \begin{split}
       \min_{{\bf w} \in \mathbb{R}^d} \Big( f({\bf w}) & = \frac 12 {\mathbf E}_{{\bf z}} \left[ f_{\bf z} \right] \Big), \\
       \quad f_{\bf z}  = \|y-{\bf w}^\top {\bf x} \|^2 &, \quad {\bf z}:= ({\bf x},y) \sim \P,
    \end{split}
\end{align}
where ${\bf x} \in {\mathbb{R}}^d$ is the input variable and $y \in {\mathbb{R}}$ is the response variable. We assume that ${\bf x}$ is zero-mean with covariance matrix $\mu {\mathbf I} \preceq \cov \preceq L {\mathbf I}$, where $\mu$ is positive.

\paragraph{Mixing time of SGD.}
Stochastic gradient optimizes $f$ through the following iterative scheme:
\begin{align} \label{eq:sgdupdate}
    {\bf w}_{n+1} = {\bf w}_n - \gamma \nabla f_{{\bf z}_n}({\bf w}_n)
\end{align}
where ${\bf z}_n \stackrel{\text{i.i.d.}}{\sim} \P$ and $\gamma$ is a constant stepsize. $\{ {\bf w}_k \}_{k=1}^n$ makes a time-homogeneous Markov chain in ${\mathbb{R}}^d$. Under regularity assumptions on $\cov$, the chain admits a unique invariant distribution denoted by $\pi_\gamma$ \cite{dieuleveut2017bridging}. \cite{dieuleveut2017bridging} proves that this Markov chain enjoys an exponential mixing time with rate $\mu/L$, namely
\begin{align}
    W_2^2(\nu({\bf w}_n), \pi_\gamma) \leq C (1-c\mu/L)^n
\end{align}
holds, where $W_2$ is Wasserstein-2 distance and $\nu({\bf w}_n)$ denotes the probability measure induced by random variable ${\bf w}_n$.

\paragraph{Stochastic Accelerated Gradient Descent (SAGD).}
Starting with ${\bf w}'_{1} = {\bf w}'_{0}$, stochastic accelerated gradient descent optimizes $f$ through the following recurrence
\begin{align} \label{eq:our_accelerated_sgd}
    \begin{split}
         {\bf w}'_{n+1} = {\bf w}'_n & + \beta ({\bf w}'_{n}-{\bf w}'_{n-1}) \\
         & - \gamma \nabla f_{{\bf z}_n}({\bf w}'_{n} + \alpha ({\bf w}'_n - {\bf w}'_{n-1})),
    \end{split}
\end{align}
where ${\bf z}_n \stackrel{\text{i.i.d.}}{\sim} \P$. Although the above sequence is not Markovian, $\{ \u_k := ({\bf w}'_{k},{\bf w}'_{k-1}) \}_{k=1}^n$ is a Markov chain running on ${\mathbb{R}}^{2d}$, since
\begin{align} \label{eq:recurrrence_wns}
    \u_{n+1} = {\bf A}_n \u_n+ \begin{bmatrix}
    \gamma y_n {\bf x}_n \\
    0
    \end{bmatrix}
\end{align}
where
\begin{align*}
    {\bf A}_n := \begin{bmatrix}
    (1+\beta) {\mathbf I}- (1+\alpha) \gamma {\bf x}_n{\bf x}^\top_n & \alpha \gamma {\bf x}_n {\bf x}^\top_n - \beta {\mathbf I} \\
    {\mathbf I} & 0
    \end{bmatrix}
\end{align*}

\paragraph{A mixing rate for SAGD.}
 The invariant distribution of $\{ \u_k \}_{k=1}^n$ is $(\pi_\gamma, \pi_\gamma)$, where $\pi_\gamma$ is the invariant distribution of the SGD-chain (see Lemma~\ref{lemma:invariance_sagd}). In this regard, it is natural to ask:
\begin{center}
\textit{
    When does $\{ \u_k \}_{k=1}^n$ exhibit better mixing properties compared to $\{ {\bf w}_{k} \}_{k=1}^n$? }
\end{center}
Considering that both of chains are simulating the same invariant distribution $\pi_\gamma$, a chain with better mixing properties is more computationally efficient. Notably, better mixing properties often leads to a better convergence rate for the mean and variance of ergodic average. Also, it may lead to better properties for solutions Poisson equation associated the chain. The solution of Poisson equation, in turn, plays an important role in the establishment of Central Limit Theorem for the ergodic average. In this paper, we further highlight a novel application for mixing-analysis: applications in the convergence of SAGD in realizable cases.  Motivated by these applications of mixing properties, this paper aims at characterizing mixing properties of $\{\u_k\}_{k=1}^n$. Specifically, we will prove that there exits a constant $c \in {\mathbb{R}}_+$ and a $3d \times 3d$-matrix ${\bf C}_\theta(\P)$ determined by first 4 moments of $\P$ such that
\begin{align*}
    W_2^2(\nu(\u_n), (\pi_\gamma,\pi_\gamma)) \leq c \| {\bf C}_\Theta(\P) \|^n_{\rho_\epsilon}/\epsilon
\end{align*}
holds where $\| {\bf M} \|_{\rho_\epsilon}$ denotes pseudospectrum of matrix ${\bf M}$, $\P$ belongs to a broad class of distributions, and $\Theta := \{ \gamma, \alpha, \beta \}$. By means of illustrative examples, we show how our results can be employed to derive the accelerated mixing rate $\mathcal{O}((1-\sqrt{\mu/L})^n)$. Then, we show that the mixing rate equates the convergence rate of SAGD in realizable cases. The proposed analysis is based on a novel \textit{non-asymptotic} analysis of products of random matrices. Although  the asymptotic analysis of products of i.i.d. random matrices is an old and rich literature~\cite{furstenberg1960products}, non-asymptotic analyses are rare.

\section{Invariant distribution}
Leveraging the invariance property, next lemma proves that the invariant distribution of SAGD-iterates simulates the same distribution as SGD-iterates.
\begin{lemma} \label{lemma:invariance_sagd}
 If chain $\{ \u_k \}_{k=1}^n$ obtained by the recurrence of Eq.~\eqref{eq:recurrrence_wns} admits a unique invariant distribution, then the invariant distribution is $(\pi_\gamma, \pi_\gamma)$ where $\pi_\gamma$ is the invariant distribution of SGD-iterates.
\end{lemma}
\begin{proof}
Proof is based on a simple application of the invarince property. Suppose that $\u_1 =[{\bf w}'_1,{\bf w}'_0]$ is drawn from the invariant distribution associated with $\{\u_k\}_{k=1}^n$. Since $[{\bf w}'_{2},{\bf w}'_1]$ is distributed as $[{\bf w}'_1, {\bf w}'_{0}]$, ${\bf w}'_0$ and ${\bf w}'_1$ are identically distributed. Since the invariant distribution is assume to be unique, we need to show that $(\pi_\gamma, \pi_\gamma)$, i.e. ${\bf w}'_{0,1} \sim \pi_\gamma$, is invariant with respect to the SAGD-update in Eq.~\eqref{eq:our_accelerated_sgd}. The SAGD update for the particular case of ridge-regression can be written alternatively as
\begin{align*}
    {\bf w}'_{2} & = \widehat{{\bf w}} + \left( \beta {\mathbf I} - \gamma{\bf x} {\bf x}^\top\right) \left( {\bf w}'_{1} - {\bf w}'_{0} \right),
\end{align*}
where $\widehat{{\bf w}} = {\bf w}'_1 - \gamma \nabla f_{{\bf z}} ({\bf w}'_1)$. Since $\pi_\gamma$ is the invariant with respect to SGD-update, $\widehat{{\bf w}} \sim \pi_\gamma$. It remains to prove that ${\bf w}'_2 \sim \pi_\gamma$.
Recall the definition of ${\bf \gamma}(\nu({\bf w}'_2),\nu(\widehat{{\bf w}}))$ used in $W_2$ notation at Eq.~\eqref{eq:w2_distance}. Suppose that ${\bf w}'_0 = {\bf w}'_1 \sim \pi_\gamma$. For this particular case, the joint distribution $(\nu({\bf w}'_2),\nu(\widehat{{\bf w}}))$ belongs ${\bf \gamma}(\nu({\bf w}'_2),\nu(\widehat{{\bf w}}))$, hence $W_2(\nu({\bf w}'_2),\nu(\widehat{{\bf w}}))=0$. This concludes the proof: ${\bf w}'_2 \sim \pi_\gamma$.
\end{proof}
\section{Mixing analysis}
In the last section, we prove that SAGD and SGD are simulating the same invariant distribution. Yet, the mixing time for SAGD is unknown to the best of our knowledge. In this section, we prove that SAGD-chain is geometrically ergodic.
\paragraph{A coupling analysis.} Similar to the analysis of SGD in \cite{dieuleveut2017bridging}, we propose a coupling analysis for SAGD. Consider two sequences $\{\u_{k}^{(0)}\}_{k=1}^n $ and $\{\u_{k}^{(1)}\}_{k=1}^n$ of SAGD-iterates starting from two different initial random vectors $\u^{(0)}_0$ and $\u^{(1)}_{0}$, respectively. These sequences are assumed to be coupled by sharing the same sequence of random variables $\{{\bf z}_{k}:=({\bf x}_k,y_k)\}_{k=1}^n$ in the recurrence of Eq.~\eqref{eq:our_accelerated_sgd}. More precisely, these sequences are obtained by following iterative schemes:
\begin{align*}
    \u_{n+1}^{(i)} = {\bf A}_n \u_{n}^{(i)} + \begin{bmatrix}
    \gamma y_n {\bf x}_n \\
    0
    \end{bmatrix}, \quad i \in \{1, 2 \}.
\end{align*}
Next Theorem proves that probability measures $\nu(\u_n^{(0)})$ converges to $\nu(\u^{(1)}_n)$ in an exponential rate.
\begin{theorem} \label{thm:convergence_char}
Suppose that Assumption~\ref{assume:symmetric_inputs} holds. Let $\{\u_{k}^{(0)}\}_{k=1}^n $ and $\{\u_{k}^{(1)}\}_{k=1}^n$ be two sequences of SAGD-iterates coupled with $\{{\bf z}_{k}\}_{k=1}^n$; then,
\begin{align*}
     W_2^2(\nu(\u_n^{(0)}),\nu(\u_n^{(1)}))\leq  c \| {\bf C}_{\Theta}(\P) \|_{\rho_\epsilon}^{n+1}  /\epsilon,
\end{align*}
holds for all $\epsilon>0$ where \footnote{Vectors $\boldsymbol{\sigma} \in {\mathbb{R}}^d $, $c:= 18 d^{3/2} {\mathbf E} \| \u_0^{(1)}- \u_0^{(0)}\|^2 $ and $\k \in {\mathbb{R}}^d $ are defined in Eq.~\eqref{eq:cov} and \eqref{eq:symmetric_inputs}, respectively.}
\begin{align}
    {\bf C}_{\Theta}(\P) & = \begin{bmatrix} {\bf D}_1^2 + (1+\alpha)^2 {\bf K} & 2 {\bf D}_1 & {\mathbf I} \\
     {\bf D}_1 {\bf D}_2 - \alpha(1+\alpha) {\bf K}& {\bf D}_2 & 0 \\
     {\bf D}_2^2 +\alpha^2 {\bf K} & 0 & 0
    \end{bmatrix}  \label{eq:cmatrix}
\end{align}
and
\begin{align*}\nonumber
    {\bf D}_1 & = (1+\beta) {\mathbf I} - \gamma (1+\alpha) {\textbf{diag}}(\boldsymbol{\sigma}), \\
    {\bf D}_2 & =  \alpha \gamma {\textbf{diag}}(\boldsymbol{\sigma}) - \beta {\mathbf I} \\
    {\bf K} & = \gamma^2 \left( {\textbf{diag}}(\k - 2(\boldsymbol{\sigma})^2) + \boldsymbol{\sigma} \boldsymbol{\sigma}^\top \right).
\end{align*}
\end{theorem}
Section~\ref{sec:convergence_char_proof} outlines the proof of the last Theorem. An immediate consequence of the above result is the uniqueness of the invariant distribution when $\| {\bf C}_{\Theta}(\P)\|_{\rho_\epsilon}<1$. Later, we will show how we can choose parameters to achieve a fast mixing.
\paragraph{An exponential rate for Ergodicity}
Replacing $\u^{(0)}_0 \sim (\pi_\gamma, \pi_\gamma)$ into the result of last Theorem leads to a mixing rate for SAGD. Next Corollary states this mixing result.
\begin{corollary} [Mixing of SAGD] \label{cor:mixing_SAGD}
Suppose that Assumption~\ref{assume:symmetric_inputs} holds and $\{ \u_k \}_{k=1}^n$ are obtained from Eq.~\eqref{eq:our_accelerated_sgd}; then,
\begin{align}
     W_2^2(\nu(\u_k),(\pi_\gamma,\pi_\gamma))\leq  c' \| {\bf C}_{\Theta}(\P) \|_{\rho_\epsilon}^{n+1}  /\epsilon
\end{align}
holds for all $\epsilon>0$ as long as $\| {\bf C}_{\Theta}(\P) \|_{\rho_\epsilon}<1$, where ${\bf C}_\Theta(\P)$ is defined in Eq.~\eqref{eq:cmatrix} and constant \[ c':=18d^{3/2}{\mathbf E} \left[ \| \u_0 - {\mathbf E}_{\u \sim (\pi_\gamma,\pi_\gamma)} \left[ \u \right] \| \right]. \]
\end{corollary}
  As a result, the mixing of SAGD depends on the 4th moment of the input due to the dependency of ${\bf C}_\Theta(\P)$ on vector $\k$ that goes into the Matrix ${\bf K}$, which in turn arises from the stochastic gradient estimates. This is in contrast to the convergence of (deterministic) accelerated gradient descent which depends only on the covariance matrix of the input.

\section{Spectral analysis and parameter tuning}

 A closer look at Eq.~\eqref{eq:cmatrix} conveys that a small stepsize choice reduces the contribution of $\gamma^2 {\bf K}$ in ${\bf C}_\Theta(\P)$ which thus reduces the noise effect in the convergence rate. Yet, the optimization process slows down for a small $\gamma$. But how can we find the proper choice of the stepsize to balance this trade-off?
  Given the eigenvalues of the covariance matrix $\cov$ as well as the vector $\k$, one can construct the matrix ${\bf C}_{\Theta}(\P)$ and minimize $\| {\bf C}_{\Theta} (\P)\|_{\rho_\epsilon}$ in $\Theta$. Since there are only 3-parameters to estimate, this problem can be solved using a simple grid search. We further simply this optimization problem based on the following key observation presented in the next Lemma: ${\bf C}_{\Theta}(\P)$ has a diagonal structure that can be employed to bound $\| {\bf C}_{\Theta}(\P)\|_{\rho_\epsilon}$ by spectral-radiuses of $3\times 3$-matrices.
\begin{lemma}[A bound on the mixing rate of SAGD] \label{lemma:spectral_radius_bound}
The spectral radius of matrix ${\bf C}_{\Theta}(\P)$ is bounded as
\begin{align}
    \begin{split}
        \| {\bf C}_{\Theta}(\P) \|_{\rho_\epsilon} & \leq \max_{i=1,\dots,d} \| {\mathbf J}_i(\Theta)\|_\rho + \epsilon \\
        & + 3(1+\alpha)^2 \gamma^2  \| {\textbf{diag}}(\boldsymbol{\sigma})^2 - \boldsymbol{\sigma} \boldsymbol{\sigma}^\top \|
    \end{split}
\end{align}
where ${\mathbf J}_j(\Theta)$ is a $3 \times 3$ matrix:
\begin{align*}
        {\mathbf J}_{i}(\Theta) :=
         \begin{bmatrix}
        {\mathbf J}_i^{(1)}   &  2 [{\bf D}_1]_{ii}  & 1 \\
        {\mathbf J}_i^{(2)}  & [{\bf D}_{2}]_{ii} & 0 \\
        {\mathbf J}_i^{(3)} & 0 & 0
        \end{bmatrix},
\end{align*}
and
\begin{align*}
    {\mathbf J}_i^{(1)} & := [{\bf D}_1]_{ii}^2 + (1+\alpha)^2 \gamma^2 (k_{i} - \sigma^2_i) \\
    {\mathbf J}_i^{(2)} & := [{\bf D}_1]_{ii}[{\bf D}_2]_{ii} -  \alpha (1+\alpha)  \gamma^2 (k_i - \sigma_i^2) \\
    {\mathbf J}_i^{(3)} & := [{\bf D}_{2}]_{ii}^2 + \gamma^2 \alpha^2 (k_i - \sigma_i^2).
\end{align*}
\end{lemma}
For the proof of the last lemma, we refer readers to Appendix \ref{sec:spectral_analysis}.
Assuming that stepsize $\gamma$ is sufficiently small,  $ \|{\mathbf J}_1(\Theta)\|_{\rho} \geq \|{\mathbf J}_i(\Theta)\|_{\rho} \geq \| {\mathbf J}_d(\Theta) \|_{\rho}$ holds. Hence, the stochastic acceleration ties to 4 parameters: (i) smoothness parameter $\sigma_{d}= L$, (ii) strong convexity $\sigma_{1} = \mu$ and (iii,iv) 4th order statistics $k_d$ and $k_1$. For the choice of $\Theta=\{\alpha,\beta,\gamma\}$, one can solve the following 3 dimensional problem
\begin{align}
    \begin{split}
        \min_{\Theta} \quad & \| {\mathbf J}_1(\Theta) \|_\rho \\
        \text{subject to} \quad & \| {\mathbf J}_d(\Theta) \|_\rho \leq 1-c\sqrt{\mu/L}.
    \end{split}
\end{align}
This provides us a practical method for the acceleration of mixing time based on minimal statistics from the input, including $(\sigma_1,\sigma_d)$ and $(k_1,k_d)$.

\section{Examples}
But, does the result of Theorem~\ref{thm:convergence_char} lead to an accelerated mixing time for SAGD, faster than the mixing rate of SGD? By means of two examples, we show that accelerated mixing rate $\mathcal{O}((1-\sqrt{\mu/L})^n)$ is achievable.

\begin{example}[Gaussian inputs] \label{exam:gaussian_example} ${\bf x}$ is a 2-dimensional Gaussian random vector with zero mean, i.e. ${{\bf x} \sim {\mathcal N}(0, {\textbf{diag}}([\mu,1]))}$.
\end{example}

By combining Lemma~\ref{lemma:spectral_radius_bound} and Theorem~\ref{thm:convergence_char}, next Lemma established the accelerated rate.
\begin{lemma}[The acceleration on Example~\ref{exam:gaussian_example}]  \label{lemma:gaussian_gaurantees}
Suppose input and label distributions are those of Example~\ref{exam:gaussian_example}.
Consider SAGD with parameters: $\alpha = 2$, $\beta=1-10^{-1/2}\sqrt{\mu}$ and $\gamma = 0.1$. Then,
\begin{multline*}
    W_2^2(\nu(\u_n),(\pi_\gamma,\pi_\gamma)) \\ \leq \frac{1200}{\sqrt{\mu}} \left( 1- \sqrt{\mu}/5  \right)^n {\mathbf E}\| {\bf w}_0 - {\bf w}_* \|^2.
\end{multline*}
holds as long as $\mu\leq0.02$.
\end{lemma}
We postpone the proof to Section \ref{sec:converge_examples_app} in the appendix.
 The result of the last lemma highlights that the sequence of SAGD-iterates enjoys better mixing properties compared to those of SGD -- if the parameters are chosen properly. Remarkably, our parameter choice implies that more extrapolation (i.e. $\alpha> \beta$) is needed in stochastic settings. In our experiments, we observe that this choice of $\alpha$ is indeed very important for the convergence rate. Therefore, it is very importance to tune parameters using the proposed spectral analysis in the last section. We note that above guarantee readily extends to non-Gaussian data with the same 4-order statistics, since the convergence only depends on the first four moments. Next corollary states this extension.
\begin{corollary}
Suppose ${\bf x} \in {\mathbb{R}}^2$ is generated from an orthogonal transformation of a random vector ${\bf z}$, i.e. ${\bf x} = {\bf U} {\bf z}$. If the coordinates of ${\bf z}$ are drawn from a symmetric zero-kurtosis distribution, then the result of Lemma~\ref{lemma:gaussian_gaurantees} holds (with the same rate using the same parameters).
\end{corollary}
Note that the zero-kurtosis property in the last corollary guarantees that the first 4 moments of the input distribution match those of a Gaussian distribution. Yet, this condition is not necessary for the accelerated mixing. Next example presents an other input distribution on which SAGD enjoys the accelerated mixing rate $\mathcal{O}((1-\sqrt{\mu/L})^n)$ in terms of Wasserstein-2 distance.
\begin{example}[Uniform-Rademacher input]
\label{exam:uniform_example}
    ${\bf x}$ is a two dimensional random variable. The first coordinate of ${\bf x}$ is a Rademacher random variable. The second coordinate is uniform on $[- \kappa^{-1/2}, \kappa^{-1/2}]$ for $\kappa<1$. In this case $\mu=1/2$ and $L=\kappa^{-1}/3$ (see Lemma~\ref{lemma:uniform_convergence_app} in Appendix).
\end{example}
The next lemma proposes an accelerated mixing rate for SAGD on the above example.
\begin{lemma}[Results on Example~\ref{exam:uniform_example}]  \label{lemma:uniform_gaurantees}
Suppose the sequence $\{\u_k\}_{k=1}^n$ is obtained by running SAGD on Example~\eqref{exam:uniform_example}. If $\alpha=2$, $\beta=1-10^{-1/2}\sqrt{\kappa}$ and $\gamma=\kappa/10$, then
\begin{align*}
    W_2^2(\nu(\u_n),(\pi_\gamma,\pi_\gamma)) & \leq \frac{1200}{\sqrt{\kappa}} \left( 1- \sqrt{\kappa}/5 \right)^n \| {\bf w}_0 - {\bf w}_* \|^2
\end{align*}
holds as long as $\kappa = 2\mu/3L \leq0.02$.
\end{lemma}
 Finally, we stress the fact that the above guarantees are only exemplary. As a matter of fact, our approach can be employed for all datasets obeying Assumption~\ref{assume:symmetric_inputs}.
 \section{Applications in realizable least-squares}
As mentioned in the introduction, mixing properties play roles in the convergence of ergodic average, central limit theorem, and even optimization in over-parameterized settings. In this section, we particularly highlight applications in over-parameterized settings. This setting has attracted attentions due to recent observations in optimization of deep neural networks. Deep nets are almost perfectly optimizable in the sense that simple gradient methods achieve a zero-training error on these networks. This is often attributed to over-parameterized weight-spaces of neural networks that may contain billions parameters. Inspired by this, recent optimization-studies  have focused on the particular case of realizable models, where the minimal objective error zero is achievable ~\cite{schmidt2018fast,pillaud2017exponential}. Next lemma proves that the invariant measure of SAGA is a Dirac measure on the minimizer in realizable cases.
\begin{lemma}
Suppose that there exits ${\bf w}_* \in {\mathbb{R}}^d$ such that $f({\bf w}_*)=0$. If the invariant measure of SAGD is unique, then it equates $(\delta({\bf w}_*),\delta({\bf w}_*))$ where $\delta({\bf w}_*)$ is the Dirac measure concentrated on ${\bf w}_*$.
\end{lemma}
\begin{proof}
Since $f({\bf w}_*)=0$, $f_{\bf z}({\bf w}_*)=0$ holds almost surely. Hence $\nabla f_{\bf z}({\bf w}_*)=0$, in that $(\delta({\bf w}_*),\delta({\bf w}_*))$ is invariant with respect to SAGD-update in Eq.~\eqref{eq:our_accelerated_sgd}. Since the invariant distribution is assumed to be unique, $(\delta({\bf w}_*),\delta({\bf w}_*))$ is THE invariant distribution.
\end{proof}
Combining the result of last lemma and Corollary~\ref{thm:convergence_char} yeilds a convergence rate for SAGD in realizable cases:
\begin{align*}
    W_2^2(\nu(\u_n),(\delta({\bf w}_*),\delta({\bf w}_*))) \leq c \| {\bf C}_{\theta}(\P)\|^{n+1}_{\rho_\epsilon}/\epsilon.
\end{align*}
Hence, our established mixing rate equates the convergence rate of SAGD in realizable cases.
More interestingly, the established rate is $\mathcal{O}((1-\sqrt{\mu/L})^n)$ for examples~\ref{exam:gaussian_example} and \ref{exam:uniform_example} under the realizability assumption. Let us compare this rate with the existing established convergence rate of \cite{schmidt2018fast} for a modified stochastic accelerated scheme. As mentioned in \textit{Related Works} section, this rate relies on constant $\rho$ in strong growth condition in Eq.~\eqref{assum:sgc}. Lemma~\ref{lemma:example1support} and \ref{lemma:uniform_consequence} prove that $\rho>\mathcal{O}(\mu/L)$ for these examples, hence the established convergence of \cite{schmidt2018fast} is not better than $\mathcal{O}((1-(\mu/L)^{\sfrac{3}{2}})^n)$ on these examples. This comparison highlights the novelty of the established mixing rate as well as the sharpness of our theoretical guarantees.
\section{Proof outline for Theorem~\ref{thm:convergence_char}}\label{sec:convergence_char_proof}

In this section, we outline the proof of Theorem~\ref{thm:convergence_char}. Let $\v_k = \u_k^{(0)}- \u_k^{(1)}$. Eq.~\eqref{eq:recurrrence_wns} yields
\begin{align} \label{eq:vn}
    \v_{n} = {\bf B}_n \v_0, \quad {\bf B}_n = {\bf A}_n {\bf A}_{n-1} \dots {\bf A}_1.
\end{align}
According to the definition of $W_2$,
\begin{align*}
    W_2^2( \nu(\u^{(0)}_n), \nu(\u^{(1)}_n)) \leq {\mathbf E} \left[ \| \v_n \|^2 \right] = {\mathbf E} \left[ \| {\bf B}_n \v_0 \|^2 \right].
\end{align*}
Hence, we need to bound the spectral norm of matrix ${\bf B}_n$ to establish the desired convergence in $W_2$. Notably, ${\bf B}_n$ is obtained by products of random non-symmetric matrices.
\paragraph{Asymptotic analyses of products of random matrices.}
We can now leverage results of the well-studied field of products of random matrices (see e.g. \cite{bougerol2012products,furstenberg1960products}), which gives interesting asymptotic characterizations of $\|{\bf B}_n\|$. For example, one can show that there is a constant $\lambda_1$ (called Lyapunov exponent) such that
\begin{align*}
   \lim_{n\to \infty} \frac{1}{n} \log(\|{\bf B}_n\|) = \lambda_1
\end{align*}
holds. A straight-forward implication of this is that $\|{\bf B}_n\|$ --depending on the sign of $\lambda_1$-- grows or decays in an exponential rate. More interestingly, the asymptotic convergence rate is a constant for all random samples ${\bf B}_n$. Yet, this result is asymptotic and it is not easy to estimate the exponent $\lambda_1$. In our setting, we need an implicit formulation (or at least a bound) for $\lambda_1$ in terms of the parameters $\Theta$ such that we can tune $\Theta$ and derive a rate.
\paragraph{A non-asymptotic analysis.}
Here, we establish a non-asymptotic expression for the Lyapunov exponent $\lambda_1$ which clearly reflects the dependency to parameters $\Theta$.
The expectation of the squared norm of iterate $\v_n$ (derived in Eq.~\eqref{eq:vn}) reads as
\begin{align} \label{eq:vn_Mn}
   {\mathbf E} \left[  \|\v_n \|^2 \right] \stackrel{\eqref{eq:vn}}{=} \v_0^\top {\bf M}_n \v_0, \quad {\bf M}_n := {\mathbf E} \left[ {\bf B}_n^\top {\bf B}_n \right].
\end{align}
To bound ${\mathbf E} \|\v_n\|^2$, we thus need to bound $\|{\bf M}_n\|$. Exploiting Assumption~\ref{assume:symmetric_inputs}, one can prove that ${\bf M}_n$ obeys an interesting block structure, which is deriven in the next lemma.
\begin{lemma} \label{lemma:Mn_blockstructure}
Assuming~\ref{assume:symmetric_inputs} holds, ${\bf M}_n$ decomposes as
\begin{align} \label{eq:Mn_decomp}
    {\bf M}_{n} = \begin{bmatrix}
    {\bf U}^\top {\textbf{diag}}(\boldsymbol{ \lambda}^{(1)}_n) {\bf U} & {\bf U}^\top {\textbf{diag}}(\boldsymbol{ \lambda}^{(2)}_n) {\bf U} \\
     {\bf U}^\top {\textbf{diag}}(\boldsymbol{ \lambda}^{(2)}_n) {\bf U} & {\bf U}^\top {\textbf{diag}}(\boldsymbol{ \lambda}^{(3)}_n) {\bf U},
    \end{bmatrix}
\end{align}
where the orthogonal matrix ${\bf U}$ contains the eigenvectors of the covariance matrix $\cov$.
\end{lemma}
 Remarkably, ${\bf U}$ is independent of the number iterations. Hence, we only need to track diagonal matrices $\boldsymbol{ \lambda}_n^{(1-3)}$ in the derived expression for ${\bf M}_n$ in Eq.~\eqref{eq:Mn_decomp}. The next lemma establishes a closed-form expression for $\boldsymbol{ \lambda}_n^{(1-3)}$.
\begin{lemma} \label{lemma:blambdan}
Consider $\boldsymbol{ \lambda}_{n}^{(1)}$--$\boldsymbol{ \lambda}_n^{(3)}$ in Eq.~\eqref{eq:Mn_decomp}.  Then
\begin{align*}
   \a_{n} = {\bf C}_{\Theta}^n(\P) {\bf 1}, \quad \a_{n} =  \begin{bmatrix}
    \boldsymbol{ \lambda}_{n}^{(1)} & \boldsymbol{ \lambda}_{n}^{(2)} & \boldsymbol{ \lambda}_{n}^{(3)}
    \end{bmatrix}
\end{align*}
holds (see Eq.~\eqref{eq:cmatrix} for the exact expression of ${\bf C}_{\Theta}(\P)$).
\end{lemma}
For the proof of the last two Lemmas, we refer the reader to Corollary \ref{sec:proof_of_theorem_appendix} in the appendix. Remarkably, the result of the last Lemma allows to compute the eigenvalues of ${\bf M}_n= {\mathbf E}\left[ {\bf B}_n^\top {\bf B}_n \right] $ in a closed form\footnote{Exploiting the block diagonal structure of ${\bf M}_n$, we can extract eigenvalues values of ${\bf M}_n$ from $\a_n$.}. Furthermore, it obtains a closed form for the Lyapunov exponent: $\lambda_1^2 = \| \a_n \|_{\max}$. Combining the result of the Lemma~\ref{lemma:Mn_blockstructure} and \ref{lemma:blambdan} concludes the proof of Theorem~\ref{thm:convergence_char}: A straight-forward application of the pseudospectrum properties (developed in Lemma~\ref{lemma:pseudospectrum_power}) yields
\begin{align*} \label{eq:an_bound}
    \| \a_n \| \leq 3d \| {\bf C}_\Theta \|^{n+1}_{\rho_\epsilon}/\epsilon.
\end{align*}
The result of Lemma~\ref{lemma:Mn_blockstructure} concludes the desired bound:
\begin{align*}
    {\mathbf E} \| \v_n \|^2 \leq \v_0^\top {\bf M}_n  \v_0 \leq 18 d^{3/2} \| {\bf C}_\Theta \|^{n+1}_{\rho_\epsilon}{\mathbf E} \|\v_0\|^2 /\epsilon.
\end{align*}

\subsubsection*{\bibname}}

\usepackage[utf8]{inputenc}
\usepackage[T1]{fontenc}
\usepackage{hyperref}
\usepackage{url}
\usepackage{booktabs}
\usepackage{amsfonts}
\usepackage{nicefrac}
\usepackage{microtype}
\include{defs}

\begin{document}

\twocolumn[

\aistatstitle{Mixing of Stochastic Accelerated Gradient Descent}

\aistatsauthor{  Peiyuan Zhang* \And Hadi Daneshmand* \And  Thomas Hofmann }

\aistatsaddress{ ETH Zurich \And ETH Zurich \And ETH Zurich}
]

\begin{abstract}
  We study the mixing properties for stochastic accelerated gradient descent (SAGD) on least-squares regression. First, we show that stochastic gradient descent (SGD) and SAGD are simulating the same invariant distribution. Motivated by this, we then establish mixing rate for SAGD-iterates and compare it with those of SGD-iterates. Theoretically, we prove that the chain of SAGD iterates is geometrically ergodic --using a proper choice of parameters and under regularity assumptions on the input distribution. More specifically, we derive an explicit mixing rate depending on the first 4 moments of the data distribution. By means of illustrative examples, we prove that SAGD-iterate chain mixes faster than the chain of iterates obtained by SGD. Furthermore, we highlight applications of the established mixing rate in the convergence analysis of SAGD on realizable objectives. The proposed analysis is based on a novel \textit{non-asymptotic} analysis of products of random matrices. This theoretical result is substantiated and validated by experiments.
\end{abstract}

\input{introduction.tex}
\input{related-works.tex}
\input{preliminaries.tex}
\input{idea.tex}
\input{main-results.tex}

\input{experiments.tex}
\section{Discussions} \label{sec:discussions}
We have established the mixing rate for stochastic accelerated gradient descent on least-squares. Using examples, we have shown than that these iterates can mix faster than SGD-iterates depending on the first 4 moments of the input distribution. This result inspires two important follow-up topics: (i) mixing analysis of SAGD on more general optimization problems, and (ii) relaxing the regularity assumption~\ref{assume:symmetric_inputs} on the input distribution.

\bibliographystyle{plain}

\section{Non-realizable cases} \label{sec:non_realizable}
 We shall now consider the question to what extent our results generalize to non-realizable cases where $f({\bf w}_*) \neq 0$. Let vector $\v_k \in {\mathbb{R}}^{2d}$ be the concatenation of two consecutive iterates,  $\v_k:=[{\bf w}_{k+1},{\bf w}_k]$. Then the sequence $\{\v_k\}_{k=1}^n$ constitutes a Markov chain. The connection between stochastic optimization and Markov chains has been explored in \cite{dieuleveut2017bridging}. Theorem~\ref{thm:convergence_char} allows us to make a similar argument on this sequence of pairs of iterates. We follow the coupling approach proposed in \cite{dieuleveut2017bridging} by considering two coupling two sequences of iterates starting from different initial parameters:
\begin{equation} \label{eq:coupling}
\begin{aligned}
    {\bf w}_{n+1} &= {\bf w}_n + \beta ({\bf w}_{n}-{\bf w}_{n-1}) - \gamma \nabla f_{{\bf z}_n}({\bf w}_{n} + \alpha ({\bf w}_n - {\bf w}_{n-1})), {\bf w}_1 ={\bf w}_0 = {\bf w} \\
    {\bf w}_{n+1}'&= {\bf w}_n' + \beta ({\bf w}_{n}'-{\bf w}_{n-1}') - \gamma \nabla f_{{\bf z}_n}({\bf w}_{n}' + \alpha ({\bf w}_n' - {\bf w}_{n-1}')), {\bf w}_1'= {\bf w}_0' = {\bf w}'
    \end{aligned}
\end{equation}
The two sequences are coupled via the shared random samples $\{{\bf z}_k\}_{k=1}^n$ in the above recurrences. For least-squares, it is easy to show that these sequences converges in a rate that is same as the convergence rate in realizable settings. The next lemma summarizes this results.
\begin{lemma} [The forgetting speed in non-realizable cases] \label{lemma:forgetting}
    Consider the coupled sequences in Eq.~\eqref{eq:coupling} running on inputs following Assumption~\ref{assume:symmetric_inputs}. These sequence converges as
    \begin{align}
        {\mathbf E} \| {\bf w}'_{n} - {\bf w}_{n} \|^2 \leq 18 d^{3/2} \left(\| {\bf C}_{\Theta} \|_{\rho_\epsilon}\right)^{n+1}  \| {\bf w}' - {\bf w} \|^2/\epsilon
    \end{align}
    for all $\epsilon>0$ (${\bf C}_{\Theta}$ is introduced in Eq.~\eqref{eq:cmatrix}).
\end{lemma}

Appendix~\ref{sec:nonrealizable_app} presents the proof. The above results implies that the Markov chain $\{\v_k\}_{k=1}^n$ is stationary and forgetting its initial condition in rate $\mathcal{O}(\| {\bf C}_{\Theta} \|_{\rho_\epsilon}^n)$. This rate can be used to derive a composite rate for the convergence of tail-averaged iterates (i.e. bias-variance decomposition in \cite{dieuleveut2017harder}). For more details on the composite rate, we refer readers to Appendix~\ref{sec:nonrealizable_app}.

\section{Preliminaries}
\paragraph{Notations.}
We will repeatedly use eigenvalue decomposition of the covariance matrix $\cov$ as
\begin{align} \label{eq:cov}
    \begin{split}
        \cov = & {\bf U}^\top {\textbf{diag}}(\boldsymbol{\sigma}) {\bf U}, \quad \boldsymbol{\sigma} :=[\sigma_1, \dots, \sigma_d] , \\
        & 0 < \sigma_{1}  = \mu \leq \dots  \leq \sigma_d = L.
    \end{split}
\end{align}
Using $\cov$, we can rewrite the gradient and stochastic gradient of $f$ as follows
\begin{align} \label{eq:grads}
    \nabla f({\bf w}) & = \cov {\bf w} - {\mathbf E} \left[ y {\bf x} \right],\\ \nabla f_{\bf z}({\bf w}) &= {\bf x} {\bf x}^\top {\bf w} - y{\bf x}.
\end{align}
Let $[{\bf M}]_{ij}$ be the element $(i,j)$ of matrix ${\bf M}$.
  We will repeatedly use the compact notation $\Theta:=\{\alpha,\beta,\gamma\}$ for a set containing the hyperparameters of SAGD. $\|\cdot\|_p$ denotes $p$-norm. For the sake of simplicity we define $\|\cdot\| := \|\cdot\|_2$.

   Our mixing rate is established in terms of Wasserstein-2 metric defined on the set of probability measures on $({\mathbb{R}}^{d},{\mathcal B}({\mathbb{R}}^d))$ with bounded second moment, denoted by $\P_2({\mathbb{R}}^d)$ \footnote{Notations are borrowed from \cite{dieuleveut2017bridging}}. More precisely,
  \begin{align} \label{eq:w2_distance}
      W_2^2(\nu,\mu) = \inf_{p \in {\bf \gamma}(\nu,\mu)}  \left( \int \| \v - {\bf w} \|^2 p(d\v,d{\bf w}) \right)
  \end{align}
  where for all $p \in {\bf \gamma}(\nu,\mu)$, $\nu= \int p(.,{\bf w}) \mu(d{\bf w})$ and $\mu = \int p(\v,.)\nu(d\v)$
  . Let $\nu(Z)$ be the probability measure induced by the random variable $Z$. Notation $Z \sim \mu$ is equivalent to $\nu(Z) = \mu$.

\paragraph{Pseudospectrum} \label{sec:pseudospectrum}
Let $\sigma({\bf M})$ be the set of complex eigenvalues of the non-symmetric matrix ${\bf M}$; then the spectral radius of ${\bf M}$ is $\rho({\bf M}) = \sup\{ |{\bf z}|\; | \; {\bf z} \in \sigma({\bf M}) \}$. $\epsilon$-pseudospectrum of ${\bf M}$, denoted by  $\sigma_{\epsilon}({\bf M})$, is defined as
\begin{align*}
    \sigma_{\epsilon}({\bf M}) := \sigma({\bf M}) \cup \{ {\bf z} \in C \; | \; \| ({\bf M} - {\bf z} {\mathbf I})^{-1} \| \geq 1/\epsilon \}.
\end{align*}
Pseudospectrum of ${\bf M}$ is a genealization of spectral radius for $\epsilon$-pseudospectrum:
\begin{align*}
    \rho_\epsilon({\bf M}) = \sup\{ |{\bf z}| \; | \; {\bf z} \in \sigma_\epsilon({\bf M}) \}.
\end{align*}
As the next Lemma states, pseudospectrum bounds the spectral norm of power of non-symmetric matrices.
\begin{lemma}[Matrix power and Pseudospectrum (Theorem 9.2\cite{jensen09})] \label{lemma:pseudospectrum_power}
The following holds for any $\epsilon$ and all $n$:
\begin{align*}
    \| {\bf M}^n \| \leq \frac{(\rho_{\epsilon}({\bf M}))^{n+1}}{\epsilon}
\end{align*}
\end{lemma}
 Pseudospectrum is mainly developed for perturbation analysis of non-hermitian matrix~\cite{jensen09}. The result of next lemma shows how a pertubation of a matrix reflects in its pseudospectrum.
\begin{lemma}[Robustness of Pseudospectrum (Theorem 5.12.\cite{jensen09})] \label{sec:robust_pseudo}
For all matrices ${\bf M}$, the following holds
\begin{align*}
    \rho_{\epsilon}({\bf A} + {\bf M}) \leq \rho_{\epsilon+\|{\bf M}\|}({\bf A})
\end{align*}
\end{lemma}
Next lemma establishes the connection between pseudospectrum and spectral radius of a matrix.
\begin{lemma}[Bauer--Fike (Theorem 5.11 of \cite{jensen09})] \label{lemma:baur_fike}
Let ${\bf M}$ be a diagonalizable square matrix such that ${\bf M} = {\cal V} \Lambda {\cal V}^{-1}$. Then for $\epsilon>0$, the following holds:
\begin{align}
    \rho_\epsilon({\bf M}) \leq \rho({\bf M}) + \kappa \epsilon
\end{align}
where $\kappa$ is the condition number of ${\cal V}$.
\end{lemma}

\paragraph{Input distribution.}
This paper focuses on a structured input distribution.
\begin{assumption}[Symmetric input] \label{assume:symmetric_inputs}
${\bf x}$ is generated by an orthogonal transformation of a random vector $\v$ whose coordinates have symmetric distribution. More precisely,
\begin{align} \label{eq:symmetric_inputs}
    {\bf x} & = {\bf U} \v, \quad {\mathbf E} \left[ \v_i^2 \right] = \sigma_i, \\
    {\mathbf E} \left[ \v_i^4 \right] & = k_i, \quad \k := [k_1, \dots, k_d],
\end{align}
where ${\bf U} \in {\mathbb{R}}^{d\times d}$ is an orthogonal matrix  and $\v \in {\mathbb{R}}^d $ is from a symmetric distribution, i.e. $v_i$ is distributed as $-v_i$\footnote{Exploiting symmetricity of $\v$, one can check that ${\bf U}$ is equal to those of Eq.~\eqref{eq:cov}.}.
\end{assumption}
The above assumption simplifies our theoretical analysis. Notably, all results can be extended to the case where coordinates of $\v$, defined in the last assumption, are independent random variable. We further remark that the above assumption naturally holds in some practical applications, such as speech recognition.

\section{Related Works}

Recent results show that -- despite the potential problem of noise instability and error accumulation (see e.g. \cite{devolder2014first}) -- stochastic accelerated methods can indeed be provably faster than SGD in certain settings~\cite{dieuleveut2017harder,jain2018accelerating,schmidt2018fast}. Among these results, \cite{dieuleveut2017harder} and \cite{jain2018accelerating} focus on least-squares regression (the same setting considered in this paper), but the focus of \cite{schmidt2018fast} is on a more general setting of learning halfspaces.  \cite{dieuleveut2017harder} has shown that stochastic Nesterov's acceleration combined with stochastic averaging accelerates the convergence of stochastic gradient descent on quadratic objectives when $\mu=0$. This combination improves the convergence of stochastic gradient descent from $1/n$ to $1/n^2$ for realizable cases, where $\exists {\bf w}_*$ such that $f({\bf w}_*)=0$. \cite{schmidt2018fast} shows that accelerated stochastic gradient descent can obtain an accelerated $\mathcal{O}\left(1-\sqrt{\mu/(\rho^2L)}\right)$ rate, if the following strong growth condition over $f_{\bf z}$ holds uniformly in ${\bf w}$:
\begin{equation} \label{assum:sgc}
    {\mathbf E}_{\bf z}  \left[ \| \nabla f_{\bf z}({\bf w}) \|^2 \right] \leq \rho \| \nabla f({\bf w})\|^2.
\end{equation}
 \cite{jain2018accelerating} proves a modified version of SAGD improves the convergence of SGD. Their results rely on a statistical condition number defined as minimum number $\Tilde{\kappa}$ such that
 \begin{align} \label{eq:statistical}
     {\mathbf E} \left[ \| {\bf x} \|_{\cov^{-1}}^2 {\bf x} {\bf x}^\top  \right]  \preceq \Tilde{\kappa} \cov
 \end{align}
 holds, which allows to prove convergence with rate $\mathcal{O}(1-(\sqrt{\mu/(\Tilde{\kappa} L)}))^n$. By proving that $\Tilde{\kappa} \leq L/\mu$, \cite{jain2018accelerating} shows that their method always outperforms SGD. For Gaussian inputs, their method enjoys the the accelerated $\mathcal{O}\left(1-\sqrt{\mu/(d L)}\right)$ rate.

 Yet, the goal of this research is different from three valuable piece of works listed above, i.e. \cite{dieuleveut2017harder}, \cite{jain2018accelerating}, and \cite{schmidt2018fast}. Here, we analyze SAGD through the framework of Markov chains. Our goal is extending the established connection between Markov chain and stochastic optimization in \cite{dieuleveut2017bridging}. This paper formulates the connection religiously for SGD. We aim at extending their result to SAGD.

redpajama_set_name: RedPajamaArXiv