SGD with memory: fundamental properties and stochastic acceleration

## Abstract

Abstract An important open problem is the theoretically feasible acceleration of mini-batch SGD-type algorithms on quadratic problems with power-law spectrum. In the non-stochastic setting, the optimal exponent ξ 𝜉 \xi italic_ξ in the loss convergence L t ∼ C L ⁢ t − ξ similar-to subscript 𝐿 𝑡 subscript 𝐶 𝐿 superscript 𝑡 𝜉 L_{t}\sim C_{L}t^{-\xi} italic_L start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT ∼ italic_C start_POSTSUBSCRIPT italic_L end_POSTSUBSCRIPT italic_t start_POSTSUPERSCRIPT - italic_ξ end_POSTSUPERSCRIPT is double that in plain GD and is achievable using Heavy Ball (HB) with a suitable schedule; this no longer works in the presence of mini-batch noise. We address this challenge by considering first-order methods with an arbitrary fixed number M 𝑀 M italic_M of auxiliary velocity vectors ( memory- M 𝑀 M italic_M algorithms ). We first prove an equivalence between two forms of such algorithms and describe them in terms of suitable characteristic polynomials. Then we develop a general expansion of the loss in terms of signal and noise propagators . Using it, we show that losses of stationary stable memory- M 𝑀 M italic_M algorithms always retain the exponent ξ 𝜉 \xi italic_ξ of plain GD, but can have different constants C L subscript 𝐶 𝐿 C_{L} italic_C start_POSTSUBSCRIPT italic_L end_POSTSUBSCRIPT depending on their effective learning rate that generalizes that of HB. We prove that in memory-1 algorithms we can make C L subscript 𝐶 𝐿 C_{L} italic_C start_POSTSUBSCRIPT italic_L end_POSTSUBSCRIPT arbitrarily small while maintaining stability. As a consequence, we propose a memory-1 algorithm with a time-dependent schedule that we show heuristically and experimentally to improve the exponent ξ 𝜉 \xi italic_ξ of plain SGD.

## 1 Introduction

Optimization is a central component of deep learning, with accelerated training algorithms providing both better final results and reducing the large amounts of compute used to train modern large models. The key aspects of this optimization is a huge number of degrees of freedom and the batch-wise learning process. The theoretical understanding of even simple gradient-based algorithms such as Stochastic Gradient Descent (SGD) is currently limited, and their acceleration challenging.

Deep neural networks can be approximated with linear models when close to convergence or in wide network regimes (Lee et al., ; Fort et al., ). For regression problems with MSE loss, this leads to quadratic optimization problem. The resulting spectral distributions – the eigenvalues $\lambda_{k}$ of the Hessian and the coefficients $c_{k}$ of the expansion of the optimal solution $\mathbf{w}_{*}$ over the Hessian eigenvectors – are often well described by power-laws when numerically computed on real-world data (Cui et al. (); Bahri et al. (); Wei et al. (), see more details and related work in Sec. A). With this motivation, we focus on infinite-dimensional quadratic problems, and especially (but not exclusively) those with a power-law asymptotic form of eigenvalues $\lambda_{k}$ and coefficients $c_{k}$:

$$ $\lambda_{k}=\Lambda k^{-\nu}(1+o(1)),\quad k\to\infty;\qquad\sum_{k:\lambda_{k }<\lambda}\lambda_{k}c_{k}^{2}=Q\lambda^{\zeta}(1+o(1)),\quad\lambda\to 0.$ (1) $$

These power laws can be thought of as stricter versions of the classical capacity and source conditions (Caponnetto & De Vito, ). Convergence rates of different first-order algorithms are determined by the eigenvalue and target exponents $\nu,\zeta>0$.

Theoretically more tractable are *stationary* algorithms, i.e. with iteration-independent parameters. In the non-stochastic setting, both stationary GD and Heavy Ball (HB) exhibit the same loss rate $O(t^{-\zeta})$ (Nemirovskiy & Polyak, ). In the presence of mini-batch noise, analysis of convergence becomes more challenging, with only the recent works Berthier et al. (); Varre et al. () showing that the rate of stationary GD remains $O(t^{-\zeta})$ for sufficiently hard targets with $\zeta<2-\tfrac{1}{\nu}$, while slowing down by noise to $O(t^{-(2-\frac{1}{\nu})})$ for easier targets with $\zeta>2-\tfrac{1}{\nu}$. In the rest of the paper, we will refer to these two regimes as *signal-* and *noise-dominated*, respectively.

With non-stationary algorithms, the rate in non-stochastic setting can be accelerated to $O(t^{-2\zeta})$ using HB with a schedule of learning rate and momentum based on Jacobi polynomials (Brakhage, ). However, in the stochastic mini-batch setting we are not aware of similar acceleration results for general $\zeta,\nu$ in the literature. Moreover, in Fig. 1 we show that direct application of HB with Jacobi schedule to noisy problems consistently leads to divergence of the optimization trajectory.

Figure: Figure 1: *Divergence of Jacobi accelerated HB vs. stability of our accelerated memory-1 (AM1) method.* For both synthetic Gaussian data (left) and MNIST classification with shallow ReLU network (right), Jacobi HB enjoys the accelerated rate $L_{t}=O(t^{-2\zeta})$ in full-batch setting, but eventually diverges in stochastic one. Increasing batch size only delays the loss explosion. In contrast, our AM1 algorithm, although having weaker acceleration than full-batch Jacobi HB, is stable for small batch sizes and a long training duration. For Gaussian data with an ideal power-law spectrum, dashed lines show theoretical loss power laws with predicted exponents ($t^{-\zeta}$ for SGD and $t^{-\zeta(1+\overline{\alpha})}$ for AM1) and asymptotically match the experimental loss trajectories. Here $\overline{\alpha}$ is the schedule exponent of the effective learning rate $\alpha_{\mathrm{eff},t}\approx t^{-\overline{\alpha}}$ used by AM1 algorithm. See precise definition of the AM1 schedule and parameters in sec. 6. Also, see sec. L for figure details and extra experiments.
Refer to caption: x1.png

Previous works point to a connection between the observed divergence of Jacobi scheduled HB and a weaker (constant-level) acceleration of the stationary algorithms. Under power law spectrum (1), stationary Heavy Ball in both non-stochastic (Velikanov & Yarotsky, ) and signal-dominated stochastic (Velikanov et al., ) settings has the loss asymptotic $L_{t}\sim C_{L}t^{-\zeta}$. While the exponent $\zeta$ is fixed, the constant $C_{L}$ can vary greatly with learning rate $\alpha$ and momentum $\beta$ parameters of HB as $C_{L}\propto(\alpha_{\mathrm{eff}})^{-\zeta}$, where $\alpha_{\mathrm{eff}}=\tfrac{\alpha}{1-\beta}$ is the *effective learning rate*. On the one hand, non-stochastic setting allows to take $\beta\to 1$ without destroying the stability of the algorithm, making the constant $C_{L}$ arbitrarily small. In a sense, this mechanism is behind the Jacobi scheduled HB with the accelerated rate $O(t^{-2\zeta})$ which slowly takes $\beta_{t}\to 1$ along the optimization trajectory. On the other hand, Velikanov et al. () shows that maintaining stability under a mini-batch sampling noise requires $\alpha_{\mathrm{eff}}$ to be bounded, thus putting a lower bound on the loss constant $C_{L}$.

#### Our contributions.

In this work, we overcome the above limitation by developing a general framework to analyze loss trajectories of a broad family of first-order methods. Then, we identify a region in the space of update rules that achieves acceleration while being stable. Specifically:

- 1.
Memory-$M$ algorithms. We define a broad family of first-order algorithms that allow to store an arbitrary fixed amount $M$ of auxiliary vectors (i.e. generalized velocities), and use any linear combinations of them in gradient computation and parameter update. First, we connect these memory-$M$ algorithms with other multistep first-order methods. Then, in the mini-batch setting, we provide a unified description of their loss trajectories in terms of expansion in various combinations of elementary signal and noise *propagators*.
- 2.
Combinatorial approach to loss asymptotics. We use the obtained propagator expansion to enrich in several directions the phase diagram of mini-batch SGD converges rates, for stationary algorithms. First, we show that the same signal-dominated $O(t^{-\zeta})$ and noise-dominated $O(t^{-(2-\frac{1}{\nu})})$ rates remain valid for any memory-$M$ algorithm. Second, we give a simple interpretation of the two phases in terms of our propagator expansion: the loss in the signal- (resp, noise-) dominated phase is asymptotically mostly determined by configurations including one long signal (resp, noise) propagator. Finally, we generalize *effective learning rate* $\alpha_{\mathrm{eff}}$ from HB to all memory-$M$ algorithms, and show that it largely determines the constant $C_{L}$ in the loss asymptotic.
- 3.
Accelerated update rule. We then focus on general memory-1 algorithms. First, we find explicit expressions for loss asymptotic and full stability conditions of general stationary memory-$1$ algorithms. Then, in this space of stable algorithms, we find a region with arbitrary large effective learning rate $\alpha_{\mathrm{eff}}$, which translates into arbitrary small constant $C_{L}$ in the loss asymptotic. The update rule in this region can be obtained by adding an extra learning rate to HB. Finally, we design a schedule that travels deep within the region in a power-law manner, and demonstrate that it achieves accelerated convergence rates up to $O(t^{-\zeta(2-\frac{1}{\nu})})$ in the signal-dominated phase.

## 2 The setting

We consider the regression problem of fitting a target function $y(\mathbf{x})$ with a linear model $f(\mathbf{w},\mathbf{x})=\langle\mathbf{x},\mathbf{w}\rangle$. Inputs $\mathbf{x}$ and parameters $\mathbf{w}$ are elements of a (infinite-dimensional) Hilbert space $\mathcal{H}$, and $\langle\cdot,\cdot\rangle$ denotes the scalar product in $\mathcal{H}$. Inputs $\mathbf{x}$ are drawn from a distribution $\rho$ on $\mathcal{H}$, and we measure the loss of the model as $L(\mathbf{w})=\tfrac{1}{2}\mathbb{E}_{\mathbf{x}\sim\rho}(\langle\mathbf{x},
\mathbf{w}\rangle-y(\mathbf{x}))^{2}$. Importantly, we assume that the target function is representable in $\mathcal{H}$ as $y(\mathbf{x})=\langle\mathbf{x},\mathbf{w}_{*}\rangle$, and require only the norm $E_{\mathbf{x}\sim\rho}[y(\mathbf{x})^{2}]$ of the target but not of the optimal parameters $\mathbf{w}_{*}$ to be finite. This setting is relevant to kernel methods or overparametrized neural networks, with $\mathbf{x}$ being the mapping of original inputs to either kernel feature space or gradients of the model, respectively; see Sec. A for more details. This brings us to the quadratic objective

$$ $L(\mathbf{w})=\frac{1}{2}\langle\Delta\mathbf{w},\mathbf{H}\Delta\mathbf{w} \rangle,\qquad\mathbf{H}=\mathbb{E}_{\mathbf{x}\sim\rho}[\mathbf{x}\otimes \mathbf{x}],$ (2) $$

where $\mathbf{H}$ is the model’s Hessian and $\Delta\mathbf{w}=\mathbf{w}-\mathbf{w}_{*}$.

We will always assume that the Hessian $\mathbf{H}$ has a discrete spectrum with eigenvalues $\lambda_{k}\searrow 0,\lambda_{k}\neq 0$. We will also often (but not always) assume that the eigenvalues $\lambda_{k}$ and the respective eigenexpansion coefficients $c_{k}=\langle\mathbf{e}_{k},\mathbf{w}_{*}\rangle$ of $\mathbf{w}_{*}$ obey power laws (1) with some $\nu,\zeta>0$. We allow the “infeasible” case $0<\zeta<1$ when $\langle\mathbf{w}_{*},\mathbf{H}\mathbf{w}_{*}\rangle<\infty$ but $\|\mathbf{w}_{*}\|^{2}=\langle\mathbf{w}_{*},\mathbf{w}_{*}\rangle=\infty.$

For stochastic optimization, we consider the random mini-batch setting. At each iteration $t$ we choose a batch $B=\{\mathbf{x}_{i}\}_{i=1}^{b}$ of $|B|$ i.i.d. samples from $\rho$, and compute the gradient on the respective empirical loss $L_{B}(\mathbf{w})=\tfrac{1}{2|B|}\sum_{i=1}^{|B|}(\langle\mathbf{x}_{i},\Delta
\mathbf{w}\rangle)^{2}.$
In this work, we focus on the evolution of the mean loss $L_{t}\equiv\mathbb{E}[L(\mathbf{w}_{t})],$ averaged w.r.t. the choice of all the batches $B_{1},\ldots,B_{t}$ along the optimization trajectory.

## 3 Gradient descent with memory

#### General form of the algorithm.

We start with describing our gradient descent with memory in the non-stochastic setting. While later we will specialize to quadratic objectives $L$ as in Eq. (2), GD with memory can be defined as the first-order minimization method applicable to any loss $L$:

$$ $\begin{pmatrix}{\mathbf{w}}_{t+1}-{\mathbf{w}}_{t}\\ \mathbf{u}_{t+1}\end{pmatrix}=\begin{pmatrix}-\alpha_{t}&\mathbf{b}_{t}^{T}\\ \mathbf{c}_{t}&D_{t}\end{pmatrix}\begin{pmatrix}\nabla L({\mathbf{w}}_{t}+ \mathbf{a}_{t}^{T}\mathbf{u}_{t})\\ \mathbf{u}_{t}\end{pmatrix},\quad t=0,1,2,\ldots$ (3) $$

Here, $\mathbf{w}_{t}$ is the step-$t$ approximation to an optimal vector $\mathbf{w}_{*}$, and $\mathbf{u}_{t}$ is an auxiliary vector representing the “memory” of the optimizer. We assume that these auxiliary vectors have the form $\mathbf{u}=(\mathbf{u}^{(1)},\ldots\mathbf{u}^{(M)})^{T}$ with $\mathbf{u}^{(m)}\in\mathcal{H};$ in other words, $\mathbf{u}_{t}$ is a size-$M$ column with each component belonging to the same space $\mathcal{H}$ as the vector $\mathbf{w}_{t}$. We will refer to $M$ as the *memory size*. The parameters $\alpha_{t}$ (learning rates) are scalar, the parameters $\mathbf{a}_{t},\mathbf{b}_{t},\mathbf{c}_{t}$ are $M$-dimensional column vectors, and $D_{t}$ are $M\times M$ scalar matrices. By $\mathbf{a}_{t}^{T}\mathbf{u}_{t}$ we mean $\sum_{m=1}^{M}a^{(m)}_{t}\mathbf{u}^{(m)}_{t},$ where $a_{t}^{(m)}$ are the scalar components of the vector $\mathbf{a}_{t}$. The algorithm can be viewed as a sequence of transformations of size-$(M+1)$ column vectors $(\begin{smallmatrix}\mathbf{w}_{t}\\
\mathbf{u}_{t}\end{smallmatrix})$ with $\mathcal{H}$-valued components.

Note a few important special cases:

- 1.
The basic gradient descent (GD) $\mathbf{w}_{t+1}=\mathbf{w}_{t}-\alpha_{t}\nabla L({\mathbf{w}}_{t})$ is the memory-0 case.
- 2.
The Heavy Ball (a.k.a. GD with momentum, Polyak ()) can be defined by the two-step recurrence $\mathbf{w}_{t+1}=\mathbf{w}_{t}-\alpha_{t}\nabla L({\mathbf{w}}_{t})+\beta_{t}
(\mathbf{w}_{t}-\mathbf{w}_{t-1})$. It can be cast in the form (3) with memory $M=1$ by introducing the auxiliary vector $\mathbf{u}_{t}=\mathbf{w}_{t}-\mathbf{w}_{t-1}$:
$\begin{pmatrix}\mathbf{w}_{t+1}-\mathbf{w}_{t}\\
\mathbf{u}_{t+1}\end{pmatrix}=\begin{pmatrix}-\alpha_{t}&\beta_{t}\\
-\alpha_{t}&\beta_{t}\end{pmatrix}\begin{pmatrix}\nabla L(\mathbf{w}_{t})\\
\mathbf{u}_{t}\end{pmatrix}.$
(4)
- 3.
Averaged GD (Polyak & Juditsky, ) can be defined by $\mathbf{w}_{t}=\tfrac{1}{t}\sum_{s=1}^{t}\widetilde{\mathbf{w}}_{s}$, where $\widetilde{\mathbf{w}}_{t+1}=\widetilde{\mathbf{w}}_{t}-\alpha_{t}\nabla L({
\widetilde{\mathbf{w}}}_{t})$ (i.e., the original GD trajectory $\widetilde{\mathbf{w}}_{t}$ is subsequently averaged, usually with the purpose of noise suppression). It can be cast in the form (3) with memory $M=1$ by introducing the auxiliary vector $\mathbf{u}_{t}=\widetilde{\mathbf{w}}_{t}-\mathbf{w}_{t}$:
$\begin{pmatrix}\mathbf{w}_{t+1}-\mathbf{w}_{t}\\
\mathbf{u}_{t+1}\end{pmatrix}=\begin{pmatrix}-\frac{\alpha_{t}}{t+1}&\frac{1}{
t+1}\\
-\frac{\alpha_{t}t}{t+1}&\frac{t}{t+1}\end{pmatrix}\begin{pmatrix}\nabla L(
\mathbf{w}_{t}+\mathbf{u}_{t})\\
\mathbf{u}_{t}\end{pmatrix}.$
(5)

Flammarion & Bach () consider a wider family of memory-1 algorithms that also includes Nesterov’s accelerated gradient (Nesterov, ). Some algorithms in the literature use memory size larger than 1 (Varre & Flammarion, ).

Iterations (3) thus represent a general framework encompassing various standard algorithms. The first entry in the l.h.s. is written as $\mathbf{w}_{t+1}-\mathbf{w}_{t}$ because it is natural to have all matrix elements in Eq. (3) *shift-invariant* while the sequence $\mathbf{w}_{t}$ is naturally *shift-equivariant* (see B.1).

Clearly, there is some redundancy in the definition of the algorithms (3). In particular, any linear transformation $\mathbf{u}^{\prime}_{t}=C_{t}\mathbf{u}_{t}$ of the auxiliary vectors with some non-degenerate $M\times M$ matrices $C_{t}$ produces an equivalent algorithm, having new parameters $\mathbf{a}^{\prime}_{t},\mathbf{b}^{\prime}_{t},\mathbf{c}^{\prime}_{t},D^{
\prime}_{t}$, but preserving the trajectory $\mathbf{w}_{t}$. A large part of our work will be focused on the stationary setting, i.e. when $\alpha,\mathbf{a},\mathbf{b},\mathbf{c},D$ are $t$-independent. In this setting we will be able to largely remove redundancy by going to an equivalent form of the dynamics.

#### Quadratic losses.

We assume now a quadratic objective $L$ as in Eq. (2). In this case $\nabla L(\mathbf{w})=\mathbf{H}\Delta\mathbf{w}$ and iterations (3) acquire the form

$$ $\displaystyle\begin{pmatrix}\Delta\mathbf{w}_{t+1}\\ \mathbf{u}_{t+1}\end{pmatrix}={}$ $\displaystyle\bigg{[}\begin{pmatrix}1&\mathbf{b}_{t}^{T}\\ 0&D_{t}\end{pmatrix}+\mathbf{H}\otimes\begin{pmatrix}-\alpha_{t}\\ \mathbf{c}_{t}\end{pmatrix}(1,\mathbf{a}_{t}^{T})\bigg{]}\begin{pmatrix}\Delta \mathbf{w}_{t}\\ \mathbf{u}_{t}\end{pmatrix}.$ (6) $$

Here, if $(\begin{smallmatrix}\Delta\mathbf{w}_{t}\\
\mathbf{u}_{t}\end{smallmatrix})$ is viewed as an element of the tensor product $\mathcal{H}\otimes\mathbb{R}^{M+1},$ the operator $\mathbf{H}$ acts in the first factor, and the rank-1 matrix $(\begin{smallmatrix}-\alpha_{t}\\
\mathbf{c}_{t}\end{smallmatrix})(\begin{smallmatrix}1&\mathbf{a}_{t}^{T}\end{
smallmatrix})$ in the second. In particular, if $\Delta\mathbf{w}_{t}$ and $\mathbf{u}_{t}$ are eigenvectors of $\mathbf{H}$ with eigenvalue $\lambda,$ then

$$ $\displaystyle\begin{pmatrix}\Delta\mathbf{w}_{t+1}\\ \mathbf{u}_{t+1}\end{pmatrix}={}$ $\displaystyle S_{\lambda,t}\begin{pmatrix}\Delta\mathbf{w}_{t}\\ \mathbf{u}_{t}\end{pmatrix},\quad S_{\lambda,t}=\bigg{[}\begin{pmatrix}1& \mathbf{b}_{t}^{T}\\ 0&D_{t}\end{pmatrix}+\lambda\begin{pmatrix}-\alpha_{t}\\ \mathbf{c}_{t}\end{pmatrix}(1,\mathbf{a}_{t}^{T})\bigg{]},$ (7) $$

and the new vectors $\Delta\mathbf{w}_{t+1},\mathbf{u}_{t+1}$ are again eigenvectors of $\mathbf{H}$ with eigenvalue $\lambda$. As a result, performing the spectral decomposition of $\Delta\mathbf{w}_{t},\mathbf{u}_{t}$ reduces the original dynamics (3) acting in $\mathcal{\mathbf{missing}}{H}\otimes\mathbb{R}^{M+1}$ to a $\lambda$-indexed collection of independent dynamics each acting in $\mathbb{R}^{M+1}$.

#### Equivalence with multistep recurrences.

As noted earlier, Heavy Ball can be alternatively written either in the single-step matrix form (4) with auxiliary vectors $\mathbf{u}_{t}$, or as the two-step recurrence $\mathbf{w}_{t+1}=\mathbf{w}_{t}-\alpha_{t}\nabla L({\mathbf{w}}_{t})+\beta_{t}
(\mathbf{w}_{t}-\mathbf{w}_{t-1})$ not involving auxiliary vectors, but involving two previous vectors $\mathbf{w}_{t},\mathbf{w}_{t-1}$. A natural question is whether a similar equivalence holds for general algorithms (3). It appears that there is no such equivalence for a general loss function $L$. However, if $L$ is quadratic as in Eq. (2) and, moreover, the parameters $\alpha,\mathbf{a,b,c},D$ are $t$-independent, we can prove an equivalence of (3) to multistep iterations

$$ $\mathbf{w}_{t+M+1}=\sum_{m=0}^{M}p_{m}\mathbf{w}_{t+m}+\sum_{m=0}^{M}q_{m} \nabla L(\mathbf{w}_{t+m}),\quad t=0,1,\ldots,$ (8) $$

where $\sum_{m=0}^{M}p_{m}=1$. Such iterations can be viewed as general stationary shift-invariant algorithms linearly depending on 0’th and 1’st order information from previous $M$ steps.

###### Theorem 1 ( B.2 ) .

Suppose that the loss is quadratic as in Eq. (2).

- 1.
Given any $t$-independent parameters $\alpha\in\mathbb{R},\mathbf{a,b,c}\in\mathbb{R}^{M},D\in\mathbb{R}^{M\times M}$, define constants $(p_{m})_{m=0}^{M},(q_{m})_{m=0}^{M}$ to be the coefficients in the polynomial expansions
$\displaystyle\sum_{m=0}^{M}p_{m}x^{m}={}$
$\displaystyle x^{M+1}-(x-1)\det(x-D),$
(9)
$\displaystyle\sum_{m=0}^{M}q_{m}x^{m}={}$
$\displaystyle\det\begin{pmatrix}\mathbf{a}^{T}\mathbf{c}-\alpha&\mathbf{a}^{T}
(1-D)-\mathbf{b}^{T}\\
\mathbf{c}&x-D\end{pmatrix}.$
(10)
Then $\sum_{m=0}^{M}p_{m}=1$ and any solution $(\begin{smallmatrix}\mathbf{w}_{t}\\
\mathbf{u}_{t}\end{smallmatrix})$ of matrix dynamics (3) obeys recurrence (8).
- 2.
Conversely, given constants $(p_{m})_{m=0}^{M},(q_{m})_{m=0}^{M}$ such that $\sum_{m=0}^{M}p_{m}=1$, there exist parameters $\alpha\in\mathbb{R},\mathbf{a,b,c}\in\mathbb{R}^{M},D\in\mathbb{R}^{M\times M}$ such that any vector sequence $\mathbf{w}_{t}$ generated by recurrence (8) satisfies the respective $t$-independent memory-$M$ matrix dynamics (3) with suitable auxiliary vectors $\mathbf{u}_{t}$. Moreover, we can set $\mathbf{a}=\mathbf{0}$ and $\mathbf{b}=(1,0.\ldots,0)^{T}$.

Theorem 1 along with its proof provide a constructive correspondence between the two forms of iterations. The fact that we can set $\mathbf{a}=0$ when converting recurrence (8) to matrix iterations (3) (statement 2) shows that for quadratic losses and stationary dynamics the possibility of shifting the argument in $\nabla L$ with $\mathbf{a}\neq 0$ does not increase the expressiveness of the algorithm.

The natural non-stationary analog of Theorem 1 fails to hold in general:

###### Proposition 1 ( B.3 ) .

There exist a non-stationary dynamics (3) with $M=1$ that cannot be written in the form $\mathbf{w}_{t+2}=p_{t,1}\mathbf{w}_{t+1}+p_{t,0}\mathbf{w}_{t}+q_{t,1}\nabla L
(\mathbf{w}_{t+1})+q_{t,0}\nabla L(\mathbf{w}_{t})$, even for quadratic problems.

###### Theorem 1 ( B.2 ) .

###### Proposition 1 ( B.3 ) .

## 4 Batch SGD with memory

#### Evolution of second moments.

SGD with memory is defined similarly to GD with memory (3), but with the batch loss $L_{B}$ instead of $L.$

For quadratic loss functions, we saw the deterministic GD with memory reduced to a simple evolution (6) of the state $(\begin{smallmatrix}\Delta\mathbf{w}_{t}\\
\mathbf{u}_{t}\end{smallmatrix})$. Such a description is no longer feasible for SGD, but we can describe deterministically the evolution of the second moments $\mathbf{M}_{t}$ of the state $(\begin{smallmatrix}\Delta\mathbf{w}_{t}\\
\mathbf{u}_{t}\end{smallmatrix})$

$$ $\mathbf{M}_{t}\equiv\big{(}\begin{smallmatrix}\mathbf{\mathbf{C}}_{t}&\mathbf{ J}_{t}\\ \mathbf{J}_{t}^{T}&\mathbf{V}_{t}\end{smallmatrix}\big{)},\quad\mathbf{C}_{t} \equiv\mathbb{E}[\Delta\mathbf{w}_{t}\otimes\Delta\mathbf{w}_{t}],\quad\mathbf {J}_{t}\equiv\mathbb{E}[\Delta\mathbf{w}_{t}\otimes\mathbf{u}^{T}_{t}],\quad \mathbf{V}_{t}\equiv\mathbb{E}[\mathbf{u}_{t}\otimes\mathbf{u}_{t}^{T}].$ (11) $$

The moments $\mathbf{C}_{t}$ determine the loss as $L_{t}=\tfrac{1}{2}\operatorname{Tr}(\mathbf{H}\mathbf{C}_{t})$. A direct computation (see sec. C.1) yields

$$ $\displaystyle\mathbf{M}_{t+1}={}$ $\displaystyle\bigg{[}\begin{pmatrix}1&\mathbf{b}_{t}^{T}\\ 0&D_{t}\end{pmatrix}+\mathbf{H}\begin{pmatrix}-\alpha_{t}\\ \mathbf{c}_{t}\end{pmatrix}(1,\mathbf{a}_{t}^{T})\bigg{]}\mathbf{M}_{t}\bigg{[ }\begin{pmatrix}1&\mathbf{b}_{t}^{T}\\ 0&D_{t}\end{pmatrix}+\mathbf{H}\begin{pmatrix}-\alpha_{t}\\ \mathbf{c}_{t}\end{pmatrix}(1,\mathbf{a}_{t}^{T})\bigg{]}^{T}$ (12) $\displaystyle+\frac{1}{|B|}{\Sigma_{\rho}}\Big{(}(1,\mathbf{a}_{t}^{T})\mathbf {M}_{t}(1,\mathbf{a}_{t}^{T})^{T}\Big{)}\otimes\begin{pmatrix}-\alpha_{t}\\ \mathbf{c}_{t}\end{pmatrix}\begin{pmatrix}-\alpha_{t}&\mathbf{c}_{t}\end{ pmatrix},$ (13) $$

where

$$ ${\Sigma_{\rho}}({\mathbf{C}})=\mathbb{E}_{\mathbf{x}\sim\rho}[\langle\mathbf{x },{\mathbf{C}}\mathbf{x}\rangle\mathbf{x}\otimes\mathbf{x}]-\mathbf{H}\mathbf{ C}\mathbf{H}.$ (14) $$

The term in line (12) corresponds to the noiseless dynamics, while the term in line (13) reflects the stochastic noise vanishing in the limit of batch size $|B|\to\infty$. The map $\Sigma_{\rho}(\mathbf{C})$ preserves the positive semi-definiteness of $\mathbf{C}$ (see Lemma 1), so the term in line (13) is also positive semi-definite, and $L_{t}$ is always not less than its counterpart for the deterministic GD with memory.

For general $\rho$, the first term in $\Sigma_{\rho}$ cannot be expressed via $\mathbf{H}$, thus preventing us from describing the loss evolution in terms of the spectral information $\lambda_{k},c_{k}$. Following Velikanov et al. (), we overcome this by replacing $\Sigma_{\rho}$ with the *Spectrally Expressible (SE)* approximation

$$ $\Sigma_{\mathrm{SE}}^{(\tau_{1},\tau_{2})}(\mathbf{C})\equiv\tau_{1} \operatorname{Tr}[\mathbf{H}\mathbf{C}]\mathbf{H}-\tau_{2}\mathbf{H}\mathbf{C} \mathbf{H},$ (15) $$

with some parameters $\tau_{1},\tau_{2}$. We list basic properties of $\Sigma_{\mathrm{SE}}^{(\tau_{1},\tau_{2})}$ in Lemma 1. In particular, $\Sigma_{\mathrm{SE}}^{(\tau_{1},\tau_{2})}$ preserves semi-definiteness of $\mathbf{C}$ whenever $0\leq\tau_{1}\geq\tau_{2}$.

We can look at this replacement in two ways. First, Velikanov et al. () show that, for certain choices of $\tau_{1},\tau_{2}$, the equality $\Sigma_{\rho}=\Sigma_{\mathrm{SE}}^{(\tau_{1},\tau_{2})}$ holds *exactly* in several natural special scenarios (gaussian and translation-invariant data) and *approximately* for a number of real-world problems such as MNIST classification. Second, generalizing earlier ideas of Bordelon & Pehlevan () and “uniform kurtosis” assumptions (Varre & Flammarion, ), we show that the loss obtained with $\Sigma_{\mathrm{SE}}^{(\tau_{1},\tau_{2})}$ can serve as an upper or lower bound for the loss obtained with $\Sigma_{\rho}$:

###### Proposition 2 ( C.2 ) .

Denote by $L_{t}^{(\rho)}$ and $L_{t}^{(\tau_{1},\tau_{2})}$ the losses obtained using the same initial condition $\mathbf{M}_{0}\geq 0$ and either the map $\Sigma_{\rho}$ or the map $\Sigma_{\mathrm{SE}}^{(\tau_{1},\tau_{2})}$ in iterations (13), respectively. Suppose that $0\leq\tau_{1}\geq\tau_{2}$. Then, if $\Sigma_{\rho}(\mathbf{C})\leq\Sigma_{\mathrm{SE}}^{(\tau_{1},\tau_{2})}(
\mathbf{C})$ for all $\mathbf{C}\geq 0$, then $L_{t}^{(\rho)}\leq L_{t}^{(\tau_{1},\tau_{2})}$ for all $t$. Conversely, if $\Sigma_{\rho}(\mathbf{C})\geq\Sigma_{\mathrm{SE}}^{(\tau_{1},\tau_{2})}(
\mathbf{C})$ for all $\mathbf{C}\geq 0$, then $L_{t}^{(\rho)}\geq L_{t}^{(\tau_{1},\tau_{2})}$ for all $t$.

###### Proposition 2 ( C.2 ) .

#### Loss expansion.

We now expand $L_{t}$ by exploiting the eigenstructure of the quadratic problem. The derived $\mathbf{M}$ matrix iterations (12)–(13) can be compactly written as $\mathbf{M}_{t+1}=F_{t}\mathbf{M}_{t}$ with a linear operator $F_{t}$. Consider the matrix representations of the blocks $\mathbf{C}_{t},\mathbf{J}_{t},$ of $\mathbf{M}_{t}$ w.r.t. the eigenbasis $\mathbf{e}_{k}$ of $\mathbf{H}$. Under the SE approximation (15), the diagonal components of these block evolve independently of the off-diagonal components, which allows us to only consider the diagonal components in the sequel. Assuming $\mathbf{w}_{0}=0,\mathbf{u}_{0}=0$, we have $\Delta\mathbf{w}_{0}=-\mathbf{w}_{*}$ and accordingly can take $\mathbf{M}_{0}=(\begin{smallmatrix}\mathbf{C}_{0}&0\\
0&0\end{smallmatrix}),$ where $\mathbf{C}_{0}=\operatorname{diag}(c_{k}^{2}),c_{k}^{2}=\langle\mathbf{e}_{k},
\mathbf{w}_{*}\rangle^{2}$. Then,

$$ $L_{t}=\tfrac{1}{2}\operatorname{Tr}[\mathbf{H}\mathbf{C}_{t}]=\tfrac{1}{2} \langle\boldsymbol{\lambda}\otimes(\begin{smallmatrix}1\\ 0\end{smallmatrix})(\begin{smallmatrix}1&0\end{smallmatrix}),F_{t}F_{t-1} \cdots F_{1}(\begin{smallmatrix}\mathbf{C}_{0}&0\\ 0&0\end{smallmatrix})\rangle,$ (16) $$

where $\boldsymbol{\lambda}=\operatorname{diag}(\lambda_{1},\lambda_{2},\ldots)$ and the inner product is between matrices, $\langle A,B\rangle\equiv\operatorname{Tr}(AB^{*}).$
The evolution operator $F_{t}$ can be represented as the sum $F_{t}=P_{t}+A_{t}$ of a rank-one operator $P_{t}$ (noise term) and an operator $A_{t}$ (main term) that acts independently at each eigenvalue $\lambda_{k}$ as a linear transformation of $(1+M)\times(1+M)$ matrices (denoted below by $Z$):

$$ $\displaystyle P_{t}={}$ $\displaystyle\tfrac{\tau_{1}}{|B|}\boldsymbol{\lambda}\otimes(\begin{ smallmatrix}-\alpha_{t}\\ \mathbf{c}_{t}\end{smallmatrix})(\begin{smallmatrix}-\alpha_{t}&\mathbf{c}_{t} ^{T}\end{smallmatrix})\langle\boldsymbol{\lambda}\otimes(\begin{smallmatrix}1 \\ \mathbf{a}_{t}\end{smallmatrix})(\begin{smallmatrix}1&\mathbf{a}_{t}^{T}\end{ smallmatrix}),\cdot\rangle,$ (17) $\displaystyle A_{t}={}$ $\displaystyle(A_{t,\lambda}),\quad A_{t,\lambda}Z=S_{\lambda,t}ZS_{\lambda,t}^ {T}-\tfrac{\tau_{2}}{|B|}\lambda^{2}(\begin{smallmatrix}-\alpha_{t}\\ \mathbf{c}_{t}\end{smallmatrix})(\begin{smallmatrix}1\\ \mathbf{a}_{t}\end{smallmatrix})^{T}Z(\begin{smallmatrix}1\\ \mathbf{a}_{t}\end{smallmatrix})(\begin{smallmatrix}-\alpha_{t}\\ \mathbf{c}_{t}\end{smallmatrix})^{T},$ (18) $$

with $S_{\lambda,t}$ defined in Eq. (7). Substituting $F_{t}=P_{t}+A_{t}$ for each $t$ in Eq. (16) yields our desired expansion of $L_{t}$. In the stationary setting this is an expansion in terms of *propagators* $U_{t},U^{\prime}_{t},V_{t},V^{\prime}_{t}$:

$$ $\displaystyle L_{t}={}$ $\displaystyle\frac{1}{2}\Big{(}V_{t+1}+\sum_{m=1}^{t}\sum_{0<t_{1}<\ldots<t_{m }<t+1}U_{t+1-t_{m}}U^{\prime}_{t_{m}-t_{m-1}}U^{\prime}_{t_{m-1}-t_{m-2}} \cdots U^{\prime}_{t_{2}-t_{1}}V^{\prime}_{t_{1}}\Big{)},$ (19) $\displaystyle U_{t}^{\prime}={}$ $\displaystyle\tfrac{\tau_{1}}{|B|}\sum_{k}\lambda_{k}^{2}\langle(\begin{ smallmatrix}1\\ \mathbf{a}\end{smallmatrix})(\begin{smallmatrix}1\\ \mathbf{a}\end{smallmatrix})^{T},A^{t-1}_{\lambda_{k}}(\begin{smallmatrix}- \alpha\\ \mathbf{c}\end{smallmatrix})(\begin{smallmatrix}-\alpha\\ \mathbf{c}\end{smallmatrix})^{T}\rangle,\quad V^{\prime}_{t}=\sum_{k}\lambda_{ k}c_{k}^{2}\langle(\begin{smallmatrix}1\\ \mathbf{a}\end{smallmatrix})(\begin{smallmatrix}1\\ \mathbf{a}\end{smallmatrix})^{T},A^{t-1}_{\lambda_{k}}(\begin{smallmatrix}1&0 \\ 0&0\end{smallmatrix})\rangle,$ (20) $$

and $U_{t},V_{t}$ are obtained from $U^{\prime}_{t},V^{\prime}_{t}$ by replacing $(\begin{smallmatrix}1\\
\mathbf{a}\end{smallmatrix})(\begin{smallmatrix}1\\
\mathbf{a}\end{smallmatrix})^{T}$ with $(\begin{smallmatrix}1&0\\
0&0\end{smallmatrix})$ in the left arguments of the scalar products. The more general non-stationary expansion is similar and involves propagators $U_{t,s},U^{\prime}_{t,s},V_{t},V^{\prime}_{t}$, see Section D.1.
The *signal propagators* $V_{t}$ and $V_{t}^{\prime}$ represent the noiseless evolution of the original signal encapsulated in $\mathbf{C}_{0}$, while the *noise propagators* $U_{t},U^{\prime}_{t}$ reflect the injected sampling noise and its evolution.

## 5 Stability and convergence of stationary memory- M 𝑀 M italic_M SGD

#### Preliminaries.

We will now use expansions (19), (20) to establish various general stability and convergence properties of stationary SGD with memory. It will be convenient to assume henceforth that all the propagators $U^{\prime}_{t},U_{t},V^{\prime}_{t},V_{t}$ are nonnegative. A particular sufficient condition for this is $\tau_{1}\geq 0,\tau_{2}\leq 0,$ since for $\tau_{2}\leq 0$ the maps $A_{\lambda}$ preserve positive semi-definiteness (see Eq. (18)).

###### Proposition 3 .

If $\tau_{1}\geq 0$ and $\tau_{2}\leq 0$, then all $U^{\prime}_{t},U_{t},V^{\prime}_{t},V_{t}\geq 0.$

###### Proposition 3 .

#### Immediate divergence.

Recall that $\dim\mathcal{H}=\infty$ and so the Hessian $\mathbf{H}$ has infinitely many eigenvalues. In this case, even single propagators $U^{\prime}_{t},U_{t},V^{\prime}_{t},V_{t}$ may diverge. Note that $U^{\prime}_{1}=\tfrac{\tau_{1}}{|B|}(\mathbf{a}^{T}\mathbf{c}-\alpha)^{2}\sum_
{k}\lambda_{k}^{2}$ and $V^{\prime}_{1}=\sum_{k}\lambda_{k}c_{k}^{2}$. Accordingly, necessary conditions to avoid divergence in one step for generic $\mathbf{a},\mathbf{c},\alpha$ are $\sum_{k}\lambda_{k}^{2}<\infty$ and $\sum_{k}\lambda_{k}c_{k}^{2}<\infty$. These conditions are clearly also sufficient to ensure $U_{t},U_{t}^{\prime},V_{t},V_{t}^{\prime}<\infty$ for all $t$. Under power law spectral assumptions (1) we automatically have $\sum_{k}\lambda_{k}c_{k}^{2}<\infty$, while $\sum_{k}\lambda_{k}^{2}<\infty$ is fulfilled whenever $\nu>\tfrac{1}{2}.$

#### Convergence of the propagator expansion.

We analyze now the loss expansion (19) and relate convergence of $L_{t}$ to properties of the propagators $U_{t},U_{t}^{\prime},V_{t},V_{t}^{\prime}$. To make the theorem less cumbersome, we state it for the case when $U_{t}=U_{t}^{\prime}$ and $V_{t}=V_{t}^{\prime}$ for all $t$, which holds if $\mathbf{a}=0$. See Section F for the general version. Denote $U_{\Sigma}=\sum_{t=1}^{\infty}U_{t},V_{\Sigma}=\sum_{t=1}^{\infty}V_{t}$.

###### Theorem 2 ( E ) .

Let $L_{t}$ be given by expansion (19) with some $U_{t}=U_{t}^{\prime}\geq 0,V_{t}=V_{t}^{\prime}\geq 0.$

- 1.
[Convergence] Suppose that $U_{\Sigma}<1$. At $t\to\infty,$ if $V_{t}=O(1)$ (respectively, $V_{t}=o(1)$), then also $L_{t}=O(1)$ (respectively, $L_{t}=o(1)$).
- 2.
[Divergence] If $U_{\Sigma}>1$ and
$V_{t}>0$ for at least one $t$,
then $\sup_{t=1,2,\ldots}L_{t}=\infty.$
- 3.
[Signal-dominated regime] Suppose that there exist constants $\xi_{V},C_{V}>0$ such that $V_{t}=C_{V}t^{-\xi_{V}}(1+o(1))$ as $t\to\infty$. Suppose also that $U_{\Sigma}<1$ and $U_{t}=O(t^{-\xi_{U}})$ with some $\xi_{U}>\max(\xi_{V},1).$ Then
$L_{t}=\frac{C_{V}}{2(1-U_{\Sigma})}t^{-\xi_{V}}(1+o(1)).$
(21)
- 4.
[Noise-dominated regime] Suppose that there exist constants $\xi_{V}>\xi_{U}>1,C_{U}>0$ such that $U_{t}=C_{U}t^{-\xi_{U}}(1+o(1))$ and $V_{t}=O(t^{-\xi_{V}})$ as $t\to\infty$. Suppose also that $U_{\Sigma}<1$. Then
$L_{t}=\frac{V_{\Sigma}C_{U}}{2(1-U_{\Sigma})^{2}}t^{-\xi_{U}}(1+o(1)).$
(22)

The first two statement characterize convergence in general, while the last two describe the special regimes occurring under power-law propagator decay assumptions. Recall that in the non-stochastic GD setting the noise propagators $U_{t}$ vanish so that $L_{t}=\tfrac{1}{2}V_{t+1}$. Then, statements 1 and 2 essentially say that convergence of SGD with memory is equivalent to that of the respective non-stochastic GD plus the condition $U_{\Sigma}<1$.

Formulas (21), (22) have a simple meaning (see Fig. 2, right). In the signal-dominated case the leading contribution to $L_{t}$ in expansion (19) comes from terms having one long signal propagator and several short noise propagators. The factor $C_{V}t^{-\xi_{V}}$ in Eq. (21) comes from the long signal propagator, while the coefficient $(1-U_{\Sigma})^{-1}$ comes from the accumulated contribution of the products of the noise propagators, summed over all possible partitions. In the noise-dominated case, the leading terms have one long noise propagator contributing the factor $C_{U}t^{-\xi_{U}}$ in Eq. (22). The remaining coefficient $\tfrac{V_{\Sigma}}{(1-U_{\Sigma})^{2}}$ is the accumulated effect of remaining short propagators before and after the long noise propagator.

###### Theorem 2 ( E ) .

#### Stability of the maps A λ subscript 𝐴 𝜆 A_{\lambda} italic_A start_POSTSUBSCRIPT italic_λ end_POSTSUBSCRIPT .

Now we start analyzing the particular form (20) of our propagators $U_{t},U_{t}^{\prime},V_{t},V_{t}^{\prime}$. Note first that the linear matrix maps $A=(A_{\lambda})$ defined in Eq. (18) may not be stable in the sense of having eigenvalues larger than 1 in absolute value. As a result, this would generally lead to an exponential growth of the propagators and hence an eventual divergence of $L_{t}$.

At $\lambda=0,$ $A_{\lambda=0}$ acts as the tensor square $S_{\lambda=0}\otimes S_{\lambda=0}$, and its eigenvalues are the products $\mu_{1}\mu_{2}$, where $\mu_{1},\mu_{2}$ are the eigenvalues of $S_{\lambda=0}$. The eigenvalues of $S_{\lambda=0}$ are just 1 and the eigenvalues of $D$. We will assume henceforth that all the eigenvalues of $D$ are strictly less than 1 in absolute value. We will refer to matrices with this property as *strictly stable*. When $D$ is strictly stable, $A_{\lambda=0}$ has the simple leading eigenvalue 1 with the eigenvector $\mathbf{M}_{0}=(\begin{smallmatrix}1&0\\
0&0\end{smallmatrix}).$ All the other eigenvalues of $A_{\lambda=0}$ are less than 1 in absolute value. This condition generalizes the standard condition $|\beta|<1$ for the momentum (Tugay & Tanik, ; Roy & Shynk, ) (note also that it holds for averaging (5) at all $t$).

As $\lambda$ increases, $A_{\lambda}$ acquires a more complicated form, and its eigenvalues continuously deform from the eigenvalues of $A_{\lambda=0}$.
It is hard to give an exact condition of stability for general $\lambda$ and $\tfrac{\tau_{2}}{|B|}$, but we prove the following:

###### Theorem 3 ( G ) .

Suppose that $D$ is strictly stable. Define effective learning rate $\alpha_{\mathrm{eff}}$ by

$$ $\alpha_{\mathrm{eff}}=\alpha-\mathbf{b}^{T}(1-D)^{-1}\mathbf{c}=\frac{\sum_{m= 0}^{M}q_{m}}{\sum_{m=0}^{M}mp_{m}-M-1},$ (23) $$

and suppose that $\alpha_{\mathrm{eff}}>0$. Then for sufficiently small $\lambda$ the linear maps $A_{\lambda}$ are also strictly stable, and their leading eigenvalue is $\mu_{A,\lambda}=1-2\alpha_{\mathrm{eff}}\lambda+O(\lambda^{2}).$

###### Theorem 3 ( G ) .

#### The power-law phase diagram.

We verify now the hypotheses of Theorem 2 under the power-law spectral assumptions (1).

###### Theorem 4 ( H ) .

Keep assumptions of Theorem 3 and suppose that $A_{\lambda}$ are strictly stable for all $0<\lambda\leq\lambda_{\max}\equiv\max_{k}\lambda_{k}$.
Then, under spectral assumptions (1) with $\nu>\tfrac{1}{2}$, the values $V_{t},V_{t}^{\prime},U_{t},U_{t}^{\prime}$ given by Eqs. (20) obey, as $t\to\infty$,

$$ $\displaystyle V_{t}^{\prime},V_{t}={}$ $\displaystyle(1+o(1))Q\Gamma(\zeta+1)(2\alpha_{\mathrm{eff}}t)^{-\zeta},$ (24) $\displaystyle U^{\prime}_{t},U_{t}={}$ $\displaystyle(1+o(1))\frac{(\alpha_{\mathrm{eff}}\Lambda)^{1/\nu}\tau_{1} \Gamma(2-1/\nu)}{|B|\nu}(2t)^{1/\nu-2}.$ (25) $$

We see that the phase diagram for stationary SGD Berthier et al. () remains valid for arbitrary memory size (see Fig. 2). Moreover, the algorithm parameters affect the leading terms in $V_{t},V_{t}^{\prime},U_{t},U_{t}^{\prime}$ only though the effective learning rate $\alpha_{\mathrm{eff}}$ and the batch size $|B|$. However, the loss asymptotics (21), (22) depend not only on the leading terms in $V_{t},U_{t}$, but also on the sums $V_{\Sigma},U_{\Sigma}$. In particular, in the signal-dominated scenario, if we wish to accelerate convergence of $L_{t},$ we need to increase $\alpha_{\mathrm{eff}}$ while keeping $U_{\Sigma}<1.$ In principle, if $U_{\Sigma}<\infty,$ we can always ensure $U_{\Sigma}<1$ by taking $|B|$ large enough. One can ask, however, if $\alpha_{\mathrm{eff}}$ can be indefinitely increased while keeping $|B|$ and preserving the condition $U_{\Sigma}<1$. It was shown in Velikanov et al. () that this is impossible for SGD with momentum. We will prove, however, that this can be done even with memory $M=1$, by slightly generalizing the model.

Figure: Figure 2: Left: The phase diagram of stationary SGD with momentum from Velikanov et al. (). Theorems 2, 4 show that it remains valid for general memory size $M$. Right: Leading terms in the loss expansion (19). In the signal-dominated phase the leading terms have one long signal propagator $V^{\prime}$ and many short noise propagators $U^{\prime}$. In the noise-dominated phase the leading terms have one long noise propagator $U^{\prime}$, while the other propagators $U^{\prime}$ and the signal propagators $V^{\prime}$ are short.

###### Theorem 4 ( H ) .

## 6 Detailed investigation of memory- 1 1 1 1 algorithm

We aim to check whether memory-$1$ algorithm can achieve accelerated convergence in the signal-dominated phase $\zeta<2-\tfrac{1}{\nu}$ with a suitable choice of its five scalar hyperparameters $\alpha_{t},a_{t},b_{t},c_{t},D_{t}$. In this section, we will assume that $\tau_{2}=0$ to ensure factorization $A_{\lambda}=S_{\lambda}\otimes S_{\lambda}$, and strict stability of $A_{\lambda}$ reduces to that of $S_{\lambda}$.

Let us start with the stationary case. According to theorems 2 and 4, stable loss trajectories have the loss asymptotic $L_{t}\sim C_{L}t^{-\zeta}$ with $C_{L}\propto\alpha_{\mathrm{eff}}^{-\zeta}$. Hence, we want to maximize $\alpha_{\mathrm{eff}}(\alpha,a,b,c,D)$ under the stability constraints. The latter consist of the noiseless stability condition $|S_{\lambda}|<1$ followed by the noisy condition $U_{\Sigma}<1$. To proceed, we find it convenient to denote $\delta=1-D$ and use multi-step recurrence parametrization of Theorem 1 with $p_{1}=2-\delta$ and $q_{0},q_{1}$ given by (10). In particular, applying (23) to memory-1 case gives $\alpha_{\mathrm{eff}}=-\frac{q_{0}+q_{1}}{\delta}$. Then, we characterize the stability with

###### Theorem 5 ( I ) .

Stability of all memory-1 algorithms is fully determined by the triplet $(\delta,\alpha_{\mathrm{eff}},q_{0})$.

- 1.
Given $\lambda_{\max}>0$, the matrix $D$ and the maps $S_{\lambda}$ are strictly stable for all $0<\lambda\leq\lambda_{\max}$ if and only if $0<\delta<2$, $-q_{0}<\tfrac{\delta}{\lambda_{\max}},$ and $0<\delta\alpha_{\mathrm{eff}}<\tfrac{4-2\delta}{\lambda_{\max}}-2q_{0}.$
- 2.
Assume $0<\delta<2$ and $\alpha_{\mathrm{eff}}>0$. If $q_{1}\leq-\alpha_{\mathrm{eff}},$ then the eigenvalues of $S_{\lambda}$ are real for all $\lambda>0.$ If $q_{1}>-\alpha_{\mathrm{eff}},$ then they are either real or lie on the circle $\{x\in\mathbb{C}:|q_{1}x+q_{0}|^{2}=\delta^{2}\alpha_{\mathrm{eff}}(\alpha_{
\mathrm{eff}}+q_{1})\}$ (see Fig. 3).
- 3.
Denote $R_{\lambda}=\sum_{t=0}^{\infty}\langle(\begin{smallmatrix}1\\
\ a\end{smallmatrix})(\begin{smallmatrix}1\\
\ a\end{smallmatrix})^{T},A^{t}_{\lambda}(\begin{smallmatrix}-\alpha\\
\ c\end{smallmatrix})(\begin{smallmatrix}-\alpha\\
\ c\end{smallmatrix})^{T}\rangle$ so that $U_{\Sigma}=\tfrac{\tau_{1}}{|B|}\sum_{k}\lambda_{k}^{2}R_{\lambda_{k}}$. Then for strictly stable $S_{\lambda}$ and $\tau_{2}=0$,
$R_{\lambda}=\frac{2\lambda q_{0}^{2}+(\lambda q_{0}+2-\delta)\delta\alpha_{
\mathrm{eff}}}{\lambda(\delta+\lambda q_{0})(4-2\delta-\lambda(2q_{0}+\delta
\alpha_{\mathrm{eff}}))}.$
(26)

Let us discuss implications of the above result. First, taking into account $-q_{0}<\tfrac{\delta}{\lambda_{\mathrm{max}}}$, the last stability condition in statement 1 gives the effective learning rate bound $\alpha_{\mathrm{eff}}<\delta^{-1}\tfrac{4}{\lambda_{\max}}$. Therefore, the only possibility to have $\alpha_{\mathrm{eff}}\to\infty$ is to take $\delta\to 0$. Now consider the noise stability condition $U_{\Sigma}<1$.
We show that $R_{\lambda}$ simplifies in the limit $\delta\to 0$, leading to three spectral regions separated by $\lambda^{\mathrm{cr}}_{1}=\tfrac{\delta}{q_{0}}$ and $\lambda^{\mathrm{cr}}_{2}=\tfrac{\delta\alpha_{\mathrm{eff}}}{q_{0}^{2}}$ and having a distinct impact on the stability sum $\sum_{k=1}^{\infty}\lambda_{k}^{2}R_{\lambda_{k}}$:

$$ $\lambda^{2}R_{\lambda}\asymp\lambda\frac{\delta\alpha_{\mathrm{eff}}+\lambda q _{0}^{2}}{\delta+\lambda q_{0}}\asymp\begin{cases}\alpha_{\mathrm{eff}}\lambda ,\quad&0<\lambda<\lambda_{1}^{\mathrm{cr}}\\ \delta\alpha_{\mathrm{eff}}/q_{0},\quad&\lambda_{1}^{\mathrm{cr}}<\lambda< \lambda_{2}^{\mathrm{cr}}\\ q_{0}\lambda,\quad&\lambda_{2}^{\mathrm{cr}}<\lambda<\lambda_{\mathrm{max}}- \varepsilon,\end{cases}$ (27) $$

where $\epsilon$ is a small constant and $f\asymp g$ denotes $cg<f<c^{\prime}g$ for some $c,c^{\prime}>0$. This observation leads to a simple criterion of stable acceleration, in particular for a power-law spectrum:

###### Theorem 6 ( J ) .

Consider a spectrum $\lambda_{k}$ s.t. $\sum_{k}\lambda_{k}<\infty,\lambda_{k}<\lambda_{\mathrm{max}}=1$, and assume $\tau_{2}=0$.

- 1.
Suppose we are adjusting the parameters $\delta,\alpha_{\mathrm{eff}},q_{0}$ so that $\alpha_{\mathrm{eff}}\to\infty$ while maintaining the strict stability of $S_{\lambda}$ as described in Theorem 5, part 1. Then,
the necessary and sufficient conditions for having $\sum_{k=1}^{\infty}\lambda_{k}^{2}R_{\lambda_{k}}=O(1)$ are
$\delta\to 0;\quad q_{0}>0;\quad\sum_{k:\lambda_{k}<\frac{\delta}{q_{0}}}\!\!
\lambda_{k}=O(\tfrac{1}{\alpha_{\mathrm{eff}}});\quad\big{|}\big{\{}k:\tfrac{
\delta}{q_{0}}<\lambda_{k}<\tfrac{\delta\alpha_{\mathrm{eff}}}{q_{0}^{2}}\big{
\}}\big{|}=O\big{(}\tfrac{q_{0}}{\delta\alpha_{\mathrm{eff}}}\big{)}.$
- 2.
Conditions from item 1 can always be satisfied. In the particular case of power law $\lambda_{k}=(1+o(1))\Lambda k^{-\nu},\nu<1,$ this can be achieved by choosing $\alpha_{\mathrm{eff}}=\delta^{-h},\;q_{0}=\delta^{g}$ with any $0\leq g<1$ and $h\leq h_{\mathrm{max}}=(1-\tfrac{1}{\nu})(1-g).$ Moreover, for $g>0$ and $h<h_{\mathrm{max}}$, the noisy stability condition $U_{\Sigma}<1$ can be satisfied with any batch size $|B|$ as $\delta\to 0$.

This result fully settles our initial question of having arbitrarily small constant $C_{L}$ in the loss asymptotic, while maintaining stability.

Figure: Figure 3: The geometric position and movement of complex eigenvalues of a memory-1 matrix $S_{\lambda}$ as $\lambda$ increases from 0 to 1 (see Theorem 5, part 2). Left: Classical Heavy Ball ($q_{0}=0$, Polyak (); Tugay & Tanik ()) with $\delta=0.25,\alpha_{\mathrm{eff}}=14$. The circle of non-real eigenvalues is centered at the origin and has radius $r=\sqrt{1-\delta}$. Acceleration ($\alpha_{\mathrm{eff}}\gg 1$) requires $\delta\ll 1$ and hence a large circle, $r\approx 1$. Right: Our generalized memory-1 algorithm with $\delta=0.15,\alpha_{\mathrm{eff}}=4,q_{0}=1.3$. The accelerated regime of Theorem 6 corresponds to small circles close to point 1.

Now, let us look closer at the parameters of the obtained accelerated algorithm. First, we see from the second item of theorem 6 that, at given $\delta$, the maximal $\alpha_{\mathrm{eff}}=\delta^{-h}$ is reached with $h=1-\tfrac{1}{\nu}$ by choosing $g=0$, corresponding to $q_{0}=const$. Although these values are specific to power-law spectrum with exponent $\nu$, they suggest that in the space of algorithm parameters $(\delta,\alpha_{\mathrm{eff}},q_{0})\in\mathbb{R}^{3}$ the region $\delta^{-1}\gg\alpha_{\mathrm{eff}}\gg q_{0}=const$ is the most promising for stable acceleration.

###### Theorem 5 ( I ) .

###### Theorem 6 ( J ) .

#### Power-law schedule with accelerated rate.

Jacobi scheduled HB uses momentum $1-\beta_{t}=\delta_{t}\propto t^{-1}$ to achieve accelerated rate $O(t^{-2\zeta})$ in the noiseless setting (Brakhage, ). Heuristically, we can view this non-stationary algorithm as a sequence of stationary algorithms gradually changing so as to increase the effective learning rate $\alpha_{\mathrm{eff}}=\tfrac{\alpha}{\delta}\to\infty$ while maintaining stability.
Mimicking this connection, we take the accelerated stationary algorithm obtained above and construct a schedule $(\delta_{t},\alpha_{\mathrm{eff},t},q_{0,t})$ that also achieves improved convergence rate. We call the result *accelerated memory-1 algorithm* (AM1), heuristically derive its convergence rate below, and verify it empirically in figures 1 and  4. We additionally discuss intuition behind AM1 update rule in sec. L.

We construct AM1 using the power-law ansatz $(\delta_{t},\alpha_{\mathrm{eff},t},q_{0,t})\sim(t^{-\overline{\delta}},t^{
\overline{\alpha}},const)$. Assuming an approximate momentary stationarity as above and substituting this ansatz in Theorem 6 part 2, we get the stability condition $\overline{\alpha}\leq\overline{\delta}(1-\tfrac{1}{\nu})$. (We emphasize, however, that, strictly speaking, Theorem 6 only holds for stationary algorithms.)

Now we examine the loss convergence rate for these stable $\overline{\delta},\overline{\alpha}$. By Theorem 2, the rate of $L_{t}$ in the signal-dominated phase is the same as that of the noiseless loss $V_{t}$. We assume that this is also true for our non-stationary algorithm, so we only need to estimate $V_{t}$. In sec.  K, we replace evolution matrix $A_{\lambda}$ with its maximal eigenvalue $|\mu_{\mathrm{max}}(A_{\lambda})|$ in eq. (77), and obtain $V_{t}=O(t^{-\zeta(1+\overline{\alpha})})$ for $\overline{\delta}<1$. Summarizing, we arrive at the following convergence rate of AM1 algorithm:

$$ $L_{t}=O(t^{-\zeta(1+\overline{\alpha})}),\qquad 0<\overline{\alpha}\leq \overline{\delta}(1-\tfrac{1}{\nu}),\quad 0<\overline{\delta}<1.$ (28) $$

The maximal limiting convergence rate $O(t^{-\zeta(2-\frac{1}{\nu})})$ is obtained by setting $\overline{\delta}=1,\;\overline{\alpha}=1-\tfrac{1}{\nu}$.

## 7 Discussion

We have developed a general framework of SGD algorithms with size-$M$ memory and analyzed their convergence properties. We developed a loss expansion in terms of signal and noise propagators, and generalized the existing SGD phase diagram of power-law regimes to arbitrary $M$. In the signal-dominated regime, we have found the asymptotic loss convergence to depend on the algorithm mostly through the batch size and a scalar parameter $\alpha_{\mathrm{eff}}$ that we call the effective learning rate. We have proved that already within memory-1 algorithms the effective learning rate can be arbitrarily increased while maintaining stability. As a consequence, we proposed a time-dependent schedule for memory-1 SGD which is heuristically expected to improve the standard SGD convergence $L_{t}=O(t^{-\zeta})$ up to $O(t^{-\zeta(2-\frac{1}{\nu})})$. This conjecture was confirmed by experiments with MNIST and synthetic problems.

Interestingly, while our framework allows shifts $\mathbf{w}_{t}+\mathbf{a}^{T}\mathbf{u}_{t}$ in the arguments of computed gradients $\nabla L$, and these shifts are important for Nesterov momentum (Nesterov, ) and averaging (Polyak & Juditsky, ), we have not found $\mathbf{a}\neq 0$ to be useful in our setting of quadratic problems and mini-batch noise. The special regime of stably increasing $\alpha_{\mathrm{eff}}$ that we identified in Theorem 6
requires separation of scales between the three main memory-1 parameters unrelated to $\mathbf{a}$, and does not seem to be easily interpretable in terms of averaging or related general ideas.

Our results leave two obvious open problems: rigorously proving the accelerated convergence $L_{t}=O(t^{-\zeta(2-\frac{1}{\nu})})$ and further improving convergence by considering larger memory sizes $M.$

## Appendix A Related work

#### Linearization of neural networks.

The relevance of quadratic optimization problems to deep learning is given by various scenarios where a neural network function $f(\mathbf{w},\mathbf{x})$ can be linearized with respect to parameters $\mathbf{w}$

$$ $f_{\mathrm{lin}}(\mathbf{w},\mathbf{x})=f(\mathbf{w}_{0},\mathbf{x})+\langle \mathbf{w}-\mathbf{w}_{0},\nabla_{\mathbf{w}}f(\mathbf{w}_{0},\mathbf{x})\rangle,$ (29) $$

often referred to as kernel regime since $f_{\mathrm{lin}}(\mathbf{w},\mathbf{x})$ can be described in terms of respective neural tangent kernel (NTK) $K(\mathbf{x},\mathbf{x^{\prime}})=\langle\nabla_{\mathbf{w}}f(\mathbf{w}_{0},
\mathbf{x}),\nabla_{\mathbf{w}}f(\mathbf{w}_{0},\mathbf{x}^{\prime})\rangle$. One of such scenario is infinite width limit in NTK regime Jacot et al. (); Lee et al. (); Chizat et al. (), where linearization becomes exact. Another scenario, is an approximate linearization of practical neural network at certain stage of its training Fort et al. (); Maddox et al. (); Ortiz-Jimenez et al. (); Lee et al. (). We cover this scenario by linearized network features $\nabla_{\mathbf{w}}f(\mathbf{w}_{0},\mathbf{x})$ as inputs $\mathbf{x}$ in (2).

#### Power laws.

Both classical kernels as well as NTK often display power-law shapes of their spectral distributions when applied to real world data. Typical examples include MNIST and CIFAR on either classical kernels or NTK of different architectures such as MLP or resnet (Cui et al., ; Bahri et al., ; Lee et al., ; Canatar et al., ; Kopitkov & Indelman, ; Dou & Liang, ; Atanasov et al., ; Bordelon & Pehlevan, ; Basri et al., ; Bietti, ; Wei et al., ; Atanasov et al., ). Power-law spectrum also appears theoretically in settings involing highly oscilating Fourier harmonics (Yang & Salman, ; Spigler et al., ; Basri et al., ; Velikanov & Yarotsky, ). Finally, power-law spectrum is often used as an assumption (Caponnetto & De Vito, ) for different optimization and generalization results on kernel methods.

Importantly, the exponent $\zeta$ in the power law spectra of real world problems is typically $\zeta<1$, for example for MNIST different works report $\zeta\in[0.2,0.4]$. This implies that the respective infinite dimensional problem with $\zeta<1$ will have divergent series of coefficients $\|\mathbf{w}_{*}\|^{2}=\sum_{k}c_{k}^{2}=\infty$, so we can refer to respective targets as *infeasible*. Such setting makes vacuous many classical results relying on norm $\|\mathbf{w}_{*}\|^{2}$ (Wei et al., ; Velikanov & Yarotsky, ). Results in the present paper were developed to include the case $\zeta<1$.

#### Quadratic optimization under mini-batch noise.

The properties of stochastic optimization can greatly vary with the structure of the noise. In particular, a common setting of additive noise (Polyak & Juditsky, ; Wu et al., ; Zhu et al., ) differs substantially from multiplicative noise (Moulines & Bach, ; Bach & Moulines, ) considered in the current paper. The mini-batch setting for power-law spectrum and its phase diagram with signal-dominated $\zeta<2-\tfrac{1}{\nu}$ and noise-dominated $\zeta>2-\tfrac{1}{\nu}$ regimes was identified in (Berthier et al., ; Varre et al., ) on the level of rigorous upper bounds, and in (Velikanov et al., ) on the level of loss symptotics under spectrally expressible (SE) approximation. The latter work relies on a generating function technique and provides the asymptotic form for the partial sums $\sum_{t\leq T}tL_{t}$ rather than for the losses $L_{t}$. In contrast, the combinatorial approach of propagator expansions used in the present work allows to directly obtain the asymptotics of $L_{t}$.

As we mentioned in the main paper, SE approximation, or its $\tau_{2}=0$ version, is often used as an upper bound assumption on the noise covariance (Bordelon & Pehlevan, ; Flammarion & Bach, ; Dieuleveut et al., ; Zou et al., ; Varre & Flammarion, ; Wu et al., ). Another series of work (Paquette & Paquette, ; Paquette et al., ; Lee et al., ; Paquette et al., ; Collins-Woodfin et al., ) consider a different high-dimensional setting where SGD converges deterministic dynamics described by convolution Volterra equation. This setting differes from ours in that eigenvalues are described by a limiting density, leading to the regime of immediate divergence $\sum_{k}\lambda_{k}^{2}=\infty$, which has to be overcome by scaling the learning rate with dimension as $d^{-1}$. Apart from this difference, algebraic structure of continuous time Volterra equation is close to our discrete time propagator expansion (19), and we expect some of our results to be extendable to the aforementioned setting.

#### Accelerated algorithms.

It is a standard textbook material that Heavy Ball accelerates convergence of plain non-stochastic gradient descent on quadratic problems (Polyak, ). In the setting of non-stochastic optimization of non strictly convex quadratic problems with power law spectral data, convergence rates of different GD-related methods were first analized by Nemirovskiy & Polyak (); Nemirovsky & Polyak (). In particular, they showed that under the source condition corresponding to the second equation in (1) with some $\zeta$, the loss convergence rate of standard stationary GD was $L_{t}=O(t^{-\zeta})$, and also gave a number of other upper and lower bounds for different algorithms. Brakhage () introduced Heavy Ball with a time-dependent schedule based on Jacobi polynomials and showed that such schedules provide acceleration of the standard GD rate $L_{t}=O(t^{-\zeta})$ to $L_{t}=O(t^{-2\zeta})$ for any spectral exponent $\zeta$. This rate is optimal among algorithms with predefined schedule, though, as shown already by Nemirovskiy & Polyak (), Conjugate Gradients can provide an even faster convergence $L_{t}=O(t^{-(2+\nu)\zeta})$ (but in this algorithm the update coefficients are not predefined). See Hanke (, ); Velikanov & Yarotsky () for several other results concerning acceleration of non-stochastic gradient descent on non-strongly convex quadratic problems under power law spectral assumptions.

A series of works (Flammarion & Bach, ; Varre & Flammarion, ) consider a special schedule of memory-1 SGD parameters such that after rescaling the iterates $\Delta\widetilde{\mathbf{w}}_{t}=t\Delta\mathbf{w}_{t}$ the update rule become stationary. Then, the stability of $\widetilde{\mathbf{w}}_{t}$ leads to $O(t^{-2})$ convergence of original iterates $\Delta\mathbf{w}_{t}$ as long as $\|\Delta\mathbf{w}_{0}\|^{2}=\|\mathbf{w}_{*}\|^{2}<\infty$. In particular, for mini-batch SGD under uniform kurtosis assumption (i.e. noise is bound by our SE term), Varre & Flammarion () gives an upper bound $O(\frac{d}{t^{2}})$, where problem dimension $d$ appears from downscaling learning rate as $\alpha=O(d^{-1})$ to ensure stability. Potentially, the necessity to downscale learning rate of the algorithm of Varre & Flammarion () can be traced to third stability condition of theorem 6, with the set $\big{\{}k:\tfrac{\delta}{q_{0}}<\lambda_{k}<\tfrac{\delta\alpha_{\mathrm{eff}}
}{q_{0}^{2}}\big{\}}$ covering the whole spectrum, thus requiring small learning rate inversely proportional to the size of the set.

## Appendix B GD with memory

### B.1 Shift-invariant and equivariant vectors

The vectors $\mathbf{w}_{t}$ in the evolution (3) are *shift-equivariant* while $\mathbf{u}_{t}$ are *shift-invariant*, in the following sense. Given a loss function $L(\cdot)$ and an initial approximation $\mathbf{w}_{0},$ consider their shifts by some vector $\mathbf{v}:$ $L^{\prime}(\cdot)=L(\cdot-\mathbf{v}),\mathbf{w}_{0}^{\prime}=\mathbf{w}_{0}+
\mathbf{v}.$ Then running the algorithm (3) with the same parameters $\alpha_{t},\mathbf{a}_{t},\mathbf{b}_{t},\mathbf{c}_{t},D_{t}$ yields the shifted trajectory $\mathbf{w}_{t}^{\prime}=\mathbf{w}_{t}+\mathbf{v}$, as is commonly expected, while preserving the auxiliary trajectory $\mathbf{u}_{t}^{\prime}=\mathbf{u}_{t}$. This can be seen as a consequence of using the difference $\mathbf{w}_{t+1}-\mathbf{w}_{t}$ in the l.h.s., since $\mathbf{w}_{t}$ being shift-equivariant for all $t$ is equivalent to $\mathbf{w}_{t+1}-\mathbf{w}_{t}$ being shift-invariant for all $t$. Since the gradients $\nabla L(\mathbf{w}_{t}+\mathbf{a}_{t}^{T}\mathbf{u}_{t})$ are also shift-invariant, all the matrix elements appearing in Eq. (4) are shift-invariant, which is consistent with the mentioned in- and equi-variance of the trajectories $\mathbf{u}_{t}$ and $\mathbf{w}_{t}$.

### B.2 Proof of Theorem 1

###### Proof.

1. Consider the characteristic polynomial

$$ $\chi(x,\lambda)=\det(x-S_{\lambda})=\det\bigg{[}x-\bigg{(}\begin{pmatrix}1& \mathbf{b}^{T}\\ 0&D\end{pmatrix}+\lambda\begin{pmatrix}-\alpha\\ \mathbf{c}\end{pmatrix}(1,\mathbf{a}^{T})\bigg{)}\bigg{]}=\sum_{m,k}r_{mk}x^{m }\lambda^{k}.$ (30) $$

We can express $\chi$ in terms of $\alpha,\mathbf{a,b,c},D$ by conjugating $S_{\lambda}$ with the matrices

$$ $U=\begin{pmatrix}1&\mathbf{a}^{T}\\ 0&\mathbf{1}\end{pmatrix},\quad U^{-1}=\begin{pmatrix}1&-\mathbf{a}^{T}\\ 0&\mathbf{1}\end{pmatrix},$ (31) $$

namely

$$ $\displaystyle\chi(x,\lambda)={}$ $\displaystyle\det(x-US_{\lambda}U^{-1})$ (32) $\displaystyle={}$ $\displaystyle\det\bigg{[}x-\bigg{(}\begin{pmatrix}1&\mathbf{b}^{T}+\mathbf{a}^ {T}(D-1)\\ 0&D\end{pmatrix}+\lambda\begin{pmatrix}\mathbf{a}^{T}\mathbf{c}-\alpha\\ \mathbf{c}\end{pmatrix}(1,\mathbf{0})\bigg{)}\bigg{]}$ (33) $\displaystyle={}$ $\displaystyle\det\begin{pmatrix}x-1-\lambda(\mathbf{a}^{T}\mathbf{c}-\alpha)&- (\mathbf{b}^{T}+\mathbf{a}^{T}(D-1))\\ -\lambda\mathbf{c}&x-D\end{pmatrix}$ (34) $\displaystyle={}$ $\displaystyle\det\begin{pmatrix}x-1&-(\mathbf{b}^{T}+\mathbf{a}^{T}(D-1))\\ 0&x-D\end{pmatrix}-\lambda\det\begin{pmatrix}\mathbf{a}^{T}\mathbf{c}-\alpha&- (\mathbf{b}^{T}+\mathbf{a}^{T}(D-1))\\ \mathbf{c}&x-D\end{pmatrix}$ (35) $\displaystyle={}$ $\displaystyle(x-1)\det(x-D)-\lambda\det\begin{pmatrix}\mathbf{a}^{T}\mathbf{c} -\alpha&\mathbf{a}^{T}(1-D)-\mathbf{b}^{T}\\ \mathbf{c}&x-D\end{pmatrix}.$ (36) $$

This shows, in particular, that $\chi$ can be written in the form

$$ $\chi(x,\lambda)=x^{M+1}-\sum_{m=0}^{M}p_{m}x^{m}-\lambda\sum_{m=0}^{M}q_{m}x^{m}$ (37) $$

with some coefficients $(p_{m})_{m=0}^{M},(q_{m})_{m=0}^{M}$ such that $\sum_{m=0}^{M}p_{m}=1.$

On the other hand, by the Cayley-Hamilton theorem, $\chi(S_{\lambda},\lambda)=0.$ In the time-independent scenario $\lambda$-eigenvectors $(\begin{smallmatrix}\Delta\mathbf{w}_{t}\\
\mathbf{u}_{t}\end{smallmatrix})$ are evolved by

$$ $\begin{pmatrix}\Delta\mathbf{w}_{t+m}\\ \mathbf{u}_{t+m}\end{pmatrix}=S_{\lambda}^{m}\begin{pmatrix}\Delta\mathbf{w}_{ t}\\ \mathbf{u}_{t}\end{pmatrix}.$ (38) $$

Recalling the expansion $\chi(x,\lambda)=\sum_{m,k}r_{mk}x^{m}\lambda^{k}$ of the characteristic polynomial, it follows that

$$ $0=\Big{(}\sum_{m,k}r_{mk}S_{\lambda}^{m}\lambda^{k}\Big{)}\begin{pmatrix} \Delta\mathbf{w}_{t}\\ \mathbf{u}_{t}\end{pmatrix}=\sum_{m,k}r_{mk}\lambda^{k}\begin{pmatrix}\Delta \mathbf{w}_{t+m}\\ \mathbf{u}_{t+m}\end{pmatrix}=\sum_{m,k}r_{mk}\mathbf{H}^{k}\begin{pmatrix} \Delta\mathbf{w}_{t+m}\\ \mathbf{u}_{t+m}\end{pmatrix}.$ (39) $$

In particular, taking the top row,

$$ $\sum_{m,k}r_{mk}\mathbf{H}^{k}\Delta\mathbf{w}_{t+m}=0.$ (40) $$

Since this equation is $\lambda$-independent, it holds not only for eigenvectors $\Delta\mathbf{w}_{t+m}$, but, by linearity, for any vectors $\Delta\mathbf{w}_{t+m}$ connected by the matrix dynamics (3).

Finally, we obtain the desired identity (8) by substituting the expression (37) for $\chi$, recalling that $\sum_{m=0}^{1}p_{m}=1$, and noting that $\mathbf{H}\Delta\mathbf{w}=\nabla L(\mathbf{w}):$

$$ $\displaystyle 0={}$ $\displaystyle\Delta\mathbf{w}_{t+M+1}-\sum_{m=0}^{M}p_{m}\Delta\mathbf{w}_{t+m }-\sum_{m=0}^{M}q_{m}\mathbf{H}\Delta\mathbf{w}_{t+m}$ (41) $\displaystyle={}$ $\displaystyle\mathbf{w}_{t+M+1}-\sum_{m=0}^{M}p_{m}\mathbf{w}_{t+m}-\sum_{m=0} ^{M}q_{m}\nabla L(\mathbf{w}_{t+m}).$ (42) $$

2. The proof will consist of two components:

- a)
Construct $\alpha,\mathbf{a,b,c},D$ so as to satify Eqs. (9),(10). By part 1) of the theorem, this will ensure that the sequence $\mathbf{w}_{t}$ satisfies the recurrence (8).
- b)
Show that for any sequence $\mathbf{w}_{0},\mathbf{w}_{1},\ldots,\mathbf{w}_{M}$ required to initialize recurrence (8) one can find an initial condition $(\begin{smallmatrix}\mathbf{w}_{0}\\
\mathbf{u}_{0}\end{smallmatrix})$ for the matrix evolution (3) to produce this initial sequence.

For a), we first satisfy Eq. (9), i.e. construct a matrix $D$ such that

$$ $\det(x-D)=\frac{x^{M+1}-\sum_{m=0}^{M}p_{m}x^{m}}{x-1}=x^{M}+\sum_{m=0}^{M-1}p ^{\prime}_{m}x^{m},$ (43) $$

where

$$ $p_{m}^{\prime}=\sum_{k=0}^{m}p_{k},\quad m=0,\ldots,M-1.$ (44) $$

A $M\times M$ matrix $D$ with such $\det(x-D)$ can be given as

$$ $D=\begin{pmatrix}-p_{M-1}^{\prime}&1&0&\cdots&0\\ -p_{M-2}^{\prime}&0&1&\cdots&0\\ &&\cdots&&\\ -p_{1}^{\prime}&0&0&\cdots&1\\ -p_{0}^{\prime}&0&0&\cdots&0\end{pmatrix}.$ (45) $$

Now to satisfy Eq. (10) we choose $\mathbf{a}=\mathbf{0},\mathbf{b}=(1,0,\ldots,0)^{T}$, so that this condition becomes

$$ $\displaystyle\sum_{m=0}^{M}q_{m}x^{m}={}$ $\displaystyle\det\begin{pmatrix}-\alpha&-\mathbf{b}^{T}\\ \mathbf{c}&x-D\end{pmatrix}$ (46) $\displaystyle={}$ $\displaystyle-\alpha\det(x-D)+\begin{pmatrix}c_{1}&-1&0&\cdots&0\\ c_{2}&x&-1&\cdots&0\\ &&\cdots&&\\ c_{M-1}&0&0&\cdots&-1\\ c_{M}&0&0&\cdots&x\end{pmatrix}$ (47) $\displaystyle={}$ $\displaystyle-\alpha\det(x-D)+\sum_{m=0}^{M-1}c_{M-1-m}x^{m}.$ (48) $$

Observe that there is a unique assignment of $\alpha$ and $\mathbf{c}$ that satisfies this identity: first we assign $\alpha=-q_{M}$ so that $\sum_{m=0}^{M}q_{m}x^{m}+\alpha\det(x-D)$ is a polynomial of degree not greater than $M-1$, and its coefficients are then exactly the components of the vector $\mathbf{c}$.

Having constructed $\alpha,\mathbf{a,b,c},D$, we check now the completeness condition b). Note that since we have set $\mathbf{a}=0,$ the gradients $\nabla L(\mathbf{w}_{t}+\mathbf{a}^{T}\mathbf{u}_{t})=\nabla L(\mathbf{w}_{t})$ appearing in Eq. (3) only depend on the target sequence $\mathbf{w}_{t}$ and not on the auxiliary sequence $\mathbf{u}_{t}$. In the first step of matrix iterations we have

$$ $\mathbf{w}_{1}-\mathbf{w}_{0}=-\alpha\nabla L(\mathbf{w}_{0})+\mathbf{b}^{T} \mathbf{u}_{0}=\mathbf{v}_{1}+\mathbf{u}_{0}^{(1)}$ (49) $$

with some $\mathbf{v}_{1}$ that depends only on $\mathbf{w}_{0}$ and not on $\mathbf{u}_{0}$. In the second step, using the explicit form (45) of the matrix $D,$

$$ $\displaystyle\mathbf{w}_{2}-\mathbf{w}_{1}={}$ $\displaystyle-\alpha\nabla L(\mathbf{w}_{1})+\mathbf{b}^{T}\mathbf{u}_{1}$ (50) $\displaystyle={}$ $\displaystyle-\alpha\nabla L(\mathbf{w}_{1})+\mathbf{b}^{T}(\nabla L(\mathbf{w }_{0})\mathbf{c}+D\mathbf{u}_{0})$ (51) $\displaystyle={}$ $\displaystyle-\alpha\nabla L(\mathbf{w}_{1})+\nabla L(\mathbf{w}_{0})\mathbf{c }_{1}-p^{\prime}_{M-1}\mathbf{u}_{0}^{(1)}+\mathbf{u}^{(2)}_{0}$ (52) $\displaystyle={}$ $\displaystyle\mathbf{v}_{2}+\mathbf{u}_{0}^{(2)},$ (53) $$

where $\mathbf{v}_{2}$ is a vector depending on $\mathbf{w}_{0},\mathbf{w}_{1}$ and the first component $\mathbf{u}_{0}^{(1)}$ of $\mathbf{u}_{0},$ but not on the second component $\mathbf{u}_{0}^{(2)}.$

By repeating this process, we see that for all $m=1,2,\ldots M$ we can write

$$ $\mathbf{w}_{m}-\mathbf{w}_{m-1}=\mathbf{v}_{m}+\mathbf{u}_{0}^{(m)}$ (54) $$

with some vector $\mathbf{v}_{m}$ that depends on $\mathbf{w}_{0},\ldots,\mathbf{w}_{m-1}$ and on $\mathbf{u}_{0}^{(1)},\ldots,\mathbf{u}_{0}^{(m-1)},$ but not on $\mathbf{u}_{0}^{(m)}.$ This shows that by sequentially adjusting the components $\mathbf{u}_{0}^{(m)}$ we can ensure the desired values for all the differences $\mathbf{w}_{m}-\mathbf{w}_{m-1},m=1,\ldots,M,$ i.e. we can find the initial condition $(\begin{smallmatrix}\mathbf{w}_{0}\\
\mathbf{u}_{0}\end{smallmatrix})$ producing the desired sequence $\mathbf{w}_{0},\ldots,\mathbf{w}_{M}.$

∎

###### Proof.

### B.3 Proof of Proposition 1

Suppose that $L$ is quadratic and consider iterations (7) in the eigenspace representation, starting from some $(\begin{smallmatrix}\Delta\mathbf{w}_{0}\\
\mathbf{u}_{0}\end{smallmatrix})$. Observe that $(\begin{smallmatrix}\Delta\mathbf{w}_{t}\\
\mathbf{u}_{t}\end{smallmatrix})$ is polynomial in $\lambda$, with degree at most $t$ and the associated leading term resulting from applying the rank-1 $\lambda$-term in each step:

$$ $\begin{pmatrix}\Delta\mathbf{w}_{t}\\ \mathbf{u}_{t}\end{pmatrix}=\lambda^{t}\begin{pmatrix}-\alpha_{t}\\ c_{t}\end{pmatrix}\Big{[}\prod_{s=1}^{t-1}(c_{s}a_{s+1}-\alpha_{s})\Big{]}( \Delta\mathbf{w}_{0}+a\mathbf{u}_{0})+O(\lambda^{t-1}).$ (55) $$

In particular,

$$ $\displaystyle\Delta\mathbf{w}_{1}={}$ $\displaystyle\lambda(-\alpha_{1})(\Delta\mathbf{w}_{0}+a_{1}\mathbf{u}_{0})+O(1)$ (56) $\displaystyle\Delta\mathbf{w}_{2}={}$ $\displaystyle\lambda^{2}(-\alpha_{2})(c_{1}a_{2}-\alpha_{1})(\Delta\mathbf{w}_ {0}+a_{1}\mathbf{u}_{0})+O(\lambda)$ (57) $$

Suppose now that $\alpha_{1}=0$ while $\alpha_{2}\neq 0$ and $c_{1}a_{2}\neq 0$. Then the leading term in $\Delta\mathbf{w}_{1}$ vanishes and the degree of $\Delta\mathbf{w}_{1}$ drops to 0, while the degree of $\Delta\mathbf{w}_{2}$ is exactly 2. Then, it is impossible for generic $\Delta\mathbf{w}_{0},\mathbf{u}_{0}$ to represent $\Delta\mathbf{w}_{2}$ as

$$ $\Delta\mathbf{w}_{2}=\lambda q_{1}\Delta\mathbf{w}_{1}+\lambda q_{0}\Delta \mathbf{w}_{0}+p_{1}\Delta\mathbf{w}_{1}+p_{0}\Delta\mathbf{w}_{0},$ (58) $$

since the r.h.s. is a polynomial in $\lambda$ of degree at most 1. Representation (58) is equivalent to

$$ $\mathbf{w}_{2}=p_{1}\mathbf{w}_{1}+p_{0}\mathbf{w}_{0}+q_{1}\nabla L(\mathbf{w }_{1})+q_{0}\nabla L(\mathbf{w}_{0})-(p_{1}+p_{0}-1)\mathbf{w}_{*}.$ (59) $$

We can assume $\mathbf{w}_{*}=0$, then this has exactly the form as in the statement of the proposition.

## Appendix C SGD with memory

### C.1 Evolution of second moments

In a single step of SGD, the parameters and memory vectors are updated with stochastic batch Hessian $\mathbf{H}_{t}=\frac{1}{|B_{t}|}\sum_{\mathbf{x}_{i}\in B_{t}}\mathbf{x}_{i}
\otimes\mathbf{x}_{i}$ as

$$ $\begin{pmatrix}\Delta\mathbf{w}_{t+1}\\ \mathbf{u}_{t+1}\end{pmatrix}=\mathbf{S}_{t}\begin{pmatrix}\Delta\mathbf{w}_{t }\\ \mathbf{u}_{t}\end{pmatrix},\quad\mathbf{S}_{t}\equiv\begin{pmatrix}1&\mathbf{ b}_{t}^{T}\\ 0&D_{t}\end{pmatrix}+\mathbf{H}_{t}\otimes\begin{pmatrix}-\alpha_{t}\\ \mathbf{c}_{t}\end{pmatrix}(1,\mathbf{a}_{t}^{T})$ (60) $$

To lighten notation, we will sometimes drop the tensor product sign between Hessian and the algorithm parameter matrix, writing simply as $\mathbf{H}_{t}\big{(}\begin{smallmatrix}-\alpha_{t}\\
\mathbf{c}_{t}\end{smallmatrix}\big{)}(1,\mathbf{a}_{t}^{T}).$

To obtain the update of second moments (11), we simply substitute the update of parameters above, and take the expectation with respect to draw of batches $B_{1},B_{2},\ldots B_{t}$

$$ $\begin{split}\mathbf{M}_{t+1}&={}\mathbb{E}_{B_{1},\ldots,B_{t-1},B_{t}}\left[ \begin{pmatrix}\Delta\mathbf{w}_{t+1}\\ \mathbf{u}_{t+1}\end{pmatrix}\begin{pmatrix}\Delta\mathbf{w}_{t+1}\\ \mathbf{u}_{t+1}\end{pmatrix}^{T}\right]\\ &=\mathbb{E}_{B_{1},\ldots,B_{t-1},B_{t}}\left[\mathbf{S}_{t}\begin{pmatrix} \Delta\mathbf{w}_{t}\\ \mathbf{u}_{t}\end{pmatrix}\begin{pmatrix}\Delta\mathbf{w}_{t}\\ \mathbf{u}_{t}\end{pmatrix}^{T}\mathbf{S}_{t}^{T}\right]=\mathbb{E}_{B_{t}} \left[\mathbf{S}_{t}\mathbf{M}_{t}\mathbf{S}_{t}^{T}\right].\end{split}$ (61) $$

To calculate the last expectation, we use that samples $\mathbf{x}_{i}\in B_{t}$ are drawn i.i.d. from distribution $\rho$. Then, for any fixed matrix $\mathbf{G}$ in $\mathcal{H}$, the expectations involving two batch Hessians are computed as

$$ $\begin{split}\mathbb{E}_{B_{t}}[\mathbf{H}_{t}\mathbf{G}\mathbf{H}_{t}]&= \mathbb{E}_{B_{t}}\left[\frac{1}{|B_{t}|^{2}}\sum_{i,j=1}^{|B_{t}|}\langle \mathbf{x}_{i},\mathbf{G}\mathbf{x}_{j}\rangle\mathbf{x}_{i}\otimes\mathbf{x}_ {j}\right]\\ &=\frac{|B_{t}|-1}{|B_{t}|}\mathbf{H}\mathbf{G}\mathbf{H}+\frac{1}{|B_{t}|} \mathbb{E}_{\mathbf{x}\sim\rho}[\langle\mathbf{x},\mathbf{G}\mathbf{x}\rangle \mathbf{x}\otimes\mathbf{x}]\\ &=\mathbf{H}\mathbf{G}\mathbf{H}+\frac{1}{|B_{t}|}\Sigma_{\rho}(\mathbf{G}), \end{split}$ (62) $$

where the sampling noise variance $\Sigma_{\rho}(\mathbf{G})$ is defined by eq. (14).

Using the above, we finish the computation of second moments update

$$ $\begin{split}\mathbf{M}_{t+1}=&\mathbb{E}_{B_{t}}\Big{(}\begin{pmatrix}1& \mathbf{b}_{t}^{T}\\ 0&D_{t}\end{pmatrix}+\mathbf{H}_{t}\begin{pmatrix}-\alpha_{t}\\ \mathbf{c}_{t}\end{pmatrix}(1,\mathbf{a}_{t}^{T})\Big{)}\mathbf{M}_{t}\Big{(} \begin{pmatrix}1&\mathbf{b}_{t}^{T}\\ 0&D_{t}\end{pmatrix}+\mathbf{H}_{t}\begin{pmatrix}-\alpha_{t}\\ \mathbf{c}_{t}\end{pmatrix}(1,\mathbf{a}_{t}^{T})\Big{)}^{T}\\ =&\Big{(}\begin{pmatrix}1&\mathbf{b}_{t}^{T}\\ 0&D_{t}\end{pmatrix}+\mathbf{H}\begin{pmatrix}-\alpha_{t}\\ \mathbf{c}_{t}\end{pmatrix}(1,\mathbf{a}_{t}^{T})\Big{)}\mathbf{M}_{t}\Big{(} \begin{pmatrix}1&\mathbf{b}_{t}^{T}\\ 0&D_{t}\end{pmatrix}+\mathbf{H}\begin{pmatrix}-\alpha_{t}\\ \mathbf{c}_{t}\end{pmatrix}(1,\mathbf{a}_{t}^{T})\Big{)}^{T}\\ &+\frac{1}{|B_{t}|}\Sigma_{\rho}\left(\mathbf{C}_{t}+2\mathbf{J}_{t}\mathbf{a} _{t}+\mathbf{a}_{t}^{T}\mathbf{V}\mathbf{a}_{t}^{T}\right)\otimes\begin{ pmatrix}-\alpha_{t}\\ \mathbf{c}_{t}\end{pmatrix}(-\alpha_{t},\mathbf{c}_{t}).\end{split}$ (63) $$

### C.2 Spectrally-expressible approximation

As usual, for two Hermitian operators $A,B$ we write $A\geq B$ (or $A\leq B$) meaning that $A-B$ is positive (negative) semi-definite.

Recall the maps $\Sigma_{\rho}$ and $\Sigma_{\mathrm{SE}}^{(\tau_{1},\tau_{2})}$ defined in Eqs. (14), (15):

$$ $\displaystyle{\Sigma_{\rho}}({\mathbf{C}})={}$ $\displaystyle\mathbb{E}_{\mathbf{x}\sim\rho}[\langle\mathbf{x},{\mathbf{C}} \mathbf{x}\rangle\mathbf{x}\otimes\mathbf{x}]-\mathbf{H}\mathbf{C}\mathbf{H},$ (64) $\displaystyle\Sigma_{\mathrm{SE}}^{(\tau_{1},\tau_{2})}(\mathbf{C})={}$ $\displaystyle\tau_{1}\operatorname{Tr}[\mathbf{H}\mathbf{C}]\mathbf{H}-\tau_{2 }\mathbf{H}\mathbf{C}\mathbf{H}.$ (65) $$

In the following lemma we collect a few useful properties of these maps.

###### Lemma 1 .

Suppose that $\mathbf{H}\geq 0.$

- 1.
For any $\mathbf{C}\geq 0$ we have $\mathbf{H}\mathbf{C}\mathbf{H}\leq\operatorname{Tr}(\mathbf{H}\mathbf{C})
\mathbf{H}$.
- 2.
For any rank-1 Hermitian $\mathbf{H}$ we have $\mathbf{H}\mathbf{C}\mathbf{H}=\operatorname{Tr}(\mathbf{H}\mathbf{C})\mathbf{
H}.$
- 3.
If $0\leq\tau_{1}\geq\tau_{2},$ then $\Sigma_{\mathrm{SE}}^{(\tau_{1},\tau_{2})}(\mathbf{C})\geq 0$ for any $\mathbf{C}\geq 0,$ and, more generally, $\Sigma_{\mathrm{SE}}^{(\tau_{1},\tau_{2})}(\mathbf{C})\geq\Sigma_{\mathrm{SE}}
^{(\tau_{1},\tau_{2})}(\mathbf{C}^{\prime})$ for any $\mathbf{C}\geq\mathbf{C}^{\prime}.$
- 4.
If $\tau_{1}<0$ or $\tau_{2}>\tau_{1}$, then it is no longer true in general that $\Sigma_{\mathrm{SE}}^{(\tau_{1},\tau_{2})}(\mathbf{C})\geq 0$ for any $\mathbf{C}\geq 0.$
- 5.
$\Sigma_{\mathrm{SE}}^{(\tau_{1}+\tau,\tau_{2}+\tau)}(\mathbf{C})\geq\Sigma_{
\mathrm{SE}}^{(\tau_{1},\tau_{2})}(\mathbf{C})$ for any $\mathbf{C}\geq 0$ and $\tau\geq 0$.
- 6.
$\Sigma_{\rho}(\mathbf{C})\geq 0$ for any $\mathbf{C}\geq 0$ and, more generally, $\Sigma_{\rho}(\mathbf{C})\geq\Sigma_{\rho}(\mathbf{C}^{\prime})$ for any $\mathbf{C}\geq\mathbf{C}^{\prime}.$

###### Proof.

- 1.
We need to show that for any vector $\mathbf{u}$
$\langle\mathbf{u},\mathbf{H}\mathbf{C}\mathbf{H}\mathbf{u}\rangle\leq
\operatorname{Tr}(\mathbf{H}\mathbf{C})\langle\mathbf{u},\mathbf{H}\mathbf{u}\rangle.$
(66)
By considering $\mathbf{u}=\mathbf{H}^{\frac{1}{2}}\mathbf{u}$ and rescaling $\mathbf{v}=\mathbf{u}/\|\mathbf{u}\|$, it is sufficient to show that $\langle\mathbf{v},\mathbf{H}^{\frac{1}{2}}\mathbf{C}\mathbf{H}^{\frac{1}{2}}
\mathbf{v}\rangle\leq\operatorname{Tr}(\mathbf{H}\mathbf{C})$ for any unit vector $\mathbf{v}$. But this just follows from the positive semi-definiteness of $\mathbf{H}^{\frac{1}{2}}\mathbf{C}\mathbf{H}^{\frac{1}{2}}$ and the cyclic property of the trace:
$\langle\mathbf{v},\mathbf{H}^{\frac{1}{2}}\mathbf{C}\mathbf{H}^{\frac{1}{2}}
\mathbf{v}\rangle\leq\operatorname{Tr}(\mathbf{H}^{\frac{1}{2}}\mathbf{C}
\mathbf{H}^{\frac{1}{2}})=\operatorname{Tr}(\mathbf{H}\mathbf{C}).$
(67)
- 2.
Immediate.
- 3.
For $\mathbf{C}\geq 0$ we have $\Sigma_{\mathrm{SE}}^{(\tau_{1},\tau_{2})}(\mathbf{C})=\tau_{1}[\operatorname{
Tr}(\mathbf{H}\mathbf{C})\mathbf{H}-\mathbf{H}\mathbf{C}\mathbf{H}]+(\tau_{2}-
\tau_{1})\mathbf{H}\mathbf{C}\mathbf{H}\geq 0$
by statement 1 and positive semi-definiteness of $\mathbf{H}\mathbf{C}\mathbf{H}.$
Monotonicity follows from the linearity of $\Sigma_{\mathrm{SE}}^{(\tau_{1},\tau_{2})}$.
- 4.
If $\tau_{1}<0,$ consider $\mathbf{H}=\mathbf{1}$ and a rank-1 $\mathbf{C}$. Then $\Sigma_{\mathrm{SE}}^{(\tau_{1},\tau_{2})}(\mathbf{C})=\tau_{1}(\operatorname{
Tr}\mathbf{C})\mathbf{1}-\tau_{2}\mathbf{C},$ which is not positive semi-definite unless the space is 1-dimensional.
If $\tau_{2}>\tau_{1}$, consider a rank-1 $\mathbf{H}$. By Statement 2, $\Sigma_{\mathrm{SE}}^{(\tau_{1},\tau_{2})}(\mathbf{C})=(\tau_{1}-\tau_{2})
\mathbf{H}\mathbf{C}\mathbf{H},$ which is negative semi-definite.
- 5.
Follows from statement 1.
- 6.
It suffices to show that for any two vectors $\mathbf{u},\mathbf{v}$ we have $\langle\mathbf{u},\Sigma_{\rho}(\mathbf{v}\otimes\mathbf{v})\mathbf{u}\rangle
\geq 0,$ i.e.
$\mathbb{E}_{\mathbf{x}\sim\rho}[\langle\mathbf{u},\mathbf{x}\rangle^{2}\langle
\mathbf{x},\mathbf{v}\rangle^{2}]\geq\mathbb{E}_{\mathbf{x}\sim\rho}^{2}[
\langle\mathbf{u},\mathbf{x}\rangle\langle\mathbf{x},\mathbf{v}\rangle].$
(68)
But this is just a special case of the Cauchy inequality $\mathbb{E}^{2}_{\mathbf{x}\sim\rho}[f]\leq\mathbb{E}_{\mathbf{x}\sim\rho}[f^{2}]$ with $f(\mathbf{x})=\langle\mathbf{u},\mathbf{x}\rangle\langle\mathbf{x},\mathbf{v}\rangle.$
Monotonicity again follows from the linearity of $\Sigma_{\rho}$.

∎

We can now prove Proposition 2.

###### Proof of Proposition 2 .

By statements 3 and 6 of Lemma, both $\Sigma_{\mathrm{SE}}^{(\tau_{1},\tau_{2})}$ and $\Sigma_{\rho}$ preserve positive semi-definiteness and are monotone in the operator sense.

Let us prove that $\Sigma_{\rho}(\mathbf{C})\leq\Sigma_{\mathrm{SE}}^{(\tau_{1},\tau_{2})}(
\mathbf{C})$ for all $\mathbf{C}\geq 0$ implies $L_{t}^{(\rho)}\leq L_{t}^{(\tau_{1},\tau_{2})}$ for all $t$.
It suffices to show that for all $t=0,1\ldots$ we have $0\leq\mathbf{C}_{t}\leq\mathbf{C}_{t}^{\prime},$ where $\mathbf{C}_{t}$ and $\mathbf{C}_{t}^{\prime}$ denote the $\mathbf{C}$-components of the matrices $\mathbf{M}_{t},\mathbf{M}_{t}^{\prime}$ corresponding to the dynamics with $\Sigma_{\rho}$ and $\Sigma_{\mathrm{SE}}^{(\tau_{1},\tau_{2})},$ respectively. In turn, it is sufficient to show that $0\leq\mathbf{M}_{t}\leq\mathbf{M}_{t}^{\prime}$ for all $t$. By induction, suppose that $0\leq\mathbf{M}_{t}\leq\mathbf{M}_{t}^{\prime}$ holds for some $t.$ Then, in particular,

$$ $0\leq(1,\mathbf{a}_{t}^{T})\mathbf{M}_{t}(1,\mathbf{a}_{t}^{T})^{T}\leq(1, \mathbf{a}_{t}^{T})\mathbf{M}_{t}^{\prime}(1,\mathbf{a}_{t}^{T})^{T}.$ (69) $$

Then, using the hypothesis of the proposition and the monotonicity of $\Sigma_{\mathrm{SE}}^{(\tau_{1},\tau_{2})}$,

$$ $\displaystyle 0\leq{}$ $\displaystyle{\Sigma_{\rho}}\Big{(}(1,\mathbf{a}_{t}^{T})\mathbf{M}_{t}(1, \mathbf{a}_{t}^{T})^{T}\Big{)}$ (70) $\displaystyle\leq{}$ $\displaystyle{\Sigma_{\mathrm{SE}}^{(\tau_{1},\tau_{2})}}\Big{(}(1,\mathbf{a}_ {t}^{T})\mathbf{M}_{t}(1,\mathbf{a}_{t}^{T})^{T}\Big{)}$ (71) $\displaystyle\leq{}$ $\displaystyle{\Sigma_{\mathrm{SE}}^{(\tau_{1},\tau_{2})}}\Big{(}(1,\mathbf{a}_ {t}^{T})\mathbf{M}^{\prime}_{t}(1,\mathbf{a}_{t}^{T})^{T}\Big{)}.$ (72) $$

It follows that the noise term (13) of the $\Sigma_{\rho}$-dynamics is positive semi-definite and dominated by the respective term of the ${\Sigma_{\mathrm{SE}}^{(\tau_{1},\tau_{2})}}$-dynamics. The same holds for the signal terms (12). It follows that $0\leq\mathbf{M}_{t+1}\leq\mathbf{M}_{t+1}^{\prime},$ as desired.

The proof of the converse statement is completely analogous.
∎

###### Lemma 1 .

###### Proof.

###### Proof of Proposition 2 .

## Appendix D Propagators and related formulas

In this section we collect auxiliary results that will be used in the proof of Theorem 2.

### D.1 General non-stationary case

In the general non-stationary case the loss expansion obtained by substituting $F_{t}=P_{t}+A_{t}$ in Eq. (16) has the form

$$ $L_{t}=\frac{1}{2}\Big{(}V_{t+1}+\sum_{m=1}^{t}\sum_{0<t_{1}<\ldots<t_{m}<t+1}U _{t+1,t_{m}}U^{\prime}_{t_{m},t_{m-1}}U^{\prime}_{t_{m-1},t_{m-2}}\cdots U^{ \prime}_{t_{2},t_{1}}V^{\prime}_{t_{1}}\Big{)}.$ (73) $$

Here the *propagators* $U^{\prime}_{t,s},U_{t,s},V^{\prime}_{t},V_{t}$ have the form

$$ $\displaystyle U^{\prime}_{t,s}={}$ $\displaystyle\tfrac{\tau_{1}}{|B|}\langle\boldsymbol{\lambda}\otimes(\begin{ smallmatrix}1\\ \mathbf{a}_{t}\end{smallmatrix})(\begin{smallmatrix}1&\mathbf{a}_{t}^{T}\end{ smallmatrix}),A_{t-1}A_{t-2}\cdots A_{s+1}[\boldsymbol{\lambda}\otimes(\begin{ smallmatrix}-\alpha_{s}\\ \mathbf{c}_{s}\end{smallmatrix})(\begin{smallmatrix}-\alpha_{s}&\mathbf{c}_{s} ^{T}\end{smallmatrix})]\rangle$ (74) $\displaystyle={}$ $\displaystyle\tfrac{\tau_{1}}{|B|}\sum_{k}\lambda_{k}^{2}\langle(\begin{ smallmatrix}1\\ \mathbf{a}_{t}\end{smallmatrix})(\begin{smallmatrix}1&\mathbf{a}_{t}^{T}\end{ smallmatrix}),A_{t-1,\lambda_{k}}\cdots A_{s+1,\lambda_{k}}(\begin{smallmatrix }-\alpha_{s}\\ \mathbf{c}_{s}\end{smallmatrix})(\begin{smallmatrix}-\alpha_{s}&\mathbf{c}_{s} ^{T}\end{smallmatrix})\rangle,$ (75) $\displaystyle V^{\prime}_{t}={}$ $\displaystyle\langle\boldsymbol{\lambda}\otimes(\begin{smallmatrix}1\\ \mathbf{a}_{t}\end{smallmatrix})(\begin{smallmatrix}1&\mathbf{a}_{t}^{T}\end{ smallmatrix}),A_{t-1}A_{t-2}\cdots A_{1}(\begin{smallmatrix}\mathbf{C}_{0}&0\\ 0&0\end{smallmatrix})\rangle$ (76) $\displaystyle={}$ $\displaystyle\sum_{k}\lambda_{k}c_{k}^{2}\langle(\begin{smallmatrix}1\\ \mathbf{a}_{t}\end{smallmatrix})(\begin{smallmatrix}1&\mathbf{a}_{t}^{T}\end{ smallmatrix}),A_{t-1,\lambda_{k}}\cdots A_{1,\lambda_{k}}(\begin{smallmatrix}1 &0\\ 0&0\end{smallmatrix})\rangle,$ (77) $$

and $U_{t,s},V_{t}$ are obtained from $U^{\prime}_{t,s},V^{\prime}_{t}$ by replacing $(\begin{smallmatrix}1\\
\mathbf{a}_{t}\end{smallmatrix})(\begin{smallmatrix}1&\mathbf{a}_{t}^{T}\end{
smallmatrix})$ with $(\begin{smallmatrix}1&0\\
0&0\end{smallmatrix})$ in the left arguments of scalar products in these formulas. In the stationary case the above formulas simplify as given in Eqs. (19), (20): we can use a single time index $U^{\prime}_{t,s}\equiv U^{\prime}_{t-s},U_{t,s}\equiv U_{t-s}$, and replace the products $A_{t-1}\ldots A_{s+1}$ by the power $A^{t-s-1}$.

We will derive a convenient representation of the vector of loss values,

$$ $\mathbf{L}=(L_{0},L_{1},\ldots),$ (78) $$

in terms of suitable linear operators. First, consider the values

$$ $L^{\prime}_{t}=\frac{1}{2}\Big{(}V^{\prime}_{t+1}+\sum_{m=1}^{t}\sum_{0<t_{1}< \ldots<t_{m}<t+1}U^{\prime}_{t+1,t_{m}}U^{\prime}_{t_{m},t_{m-1}}U^{\prime}_{t _{m-1},t_{m-2}}\cdots U^{\prime}_{t_{2},t_{1}}V^{\prime}_{t_{1}}\Big{)}$ (79) $$

and the respective vector

$$ $\mathbf{L}^{\prime}=(L^{\prime}_{0},L^{\prime}_{1},\ldots)^{T}$ (80) $$

Observe that

$$ $\mathbf{L}=\tfrac{1}{2}\mathbf{V}+\mathbb{U}\mathbf{L}^{\prime},$ (81) $$

where

$$ $\mathbf{V}=(V_{1},V_{2},\ldots)^{T}$ (82) $$

and $\mathbb{U}$ is the infinite matrix

$$ $(\mathbb{U})_{ts}=\begin{cases}U_{t+1,s+1},&t>s,\\ 0,&t\leq s.\end{cases}$ (83) $$

Note now that $\mathbf{L}^{\prime}$ satisfies the equation

$$ $\mathbf{L}^{\prime}=\tfrac{1}{2}\mathbf{V}^{\prime}+\mathbb{U}^{\prime}\mathbf {L}^{\prime},$ (84) $$

where

$$ $\displaystyle\mathbf{V}^{\prime}={}$ $\displaystyle(V^{\prime}_{1},V^{\prime}_{2},\ldots)^{T},$ (85) $\displaystyle(\mathbb{U}^{\prime})_{ts}={}$ $\displaystyle\begin{cases}U^{\prime}_{t+1,s+1},&t>s,\\ 0,&t\leq s.\end{cases}$ (86) $$

The representation (79) corresponds to the formal solution of Eq. (84):

$$ $\mathbf{L}^{\prime}=\tfrac{1}{2}(1-\mathbb{U}^{\prime})^{-1}\mathbf{V}^{\prime }=\tfrac{1}{2}(\mathbf{V}^{\prime}+\mathbb{U}^{\prime}\mathbf{V}^{\prime}+ \mathbb{U}^{\prime 2}\mathbf{V}^{\prime}+\ldots).$ (87) $$

When viewed in suitable Banach spaces, these operator representations can be conveniently used to derive upper bounds on $L_{t}$. Specifically, let $a_{t}>0$ be a positive sequence and consider the Banach space of sequences with the norm

$$ $\|\mathbf{f}\|=\sup_{t}|a_{t}f_{t}|,\quad\mathbf{f}=(f_{0},f_{1},\ldots).$ (88) $$

Then

$$ $\displaystyle\|\mathbb{U}\|={}$ $\displaystyle\sup_{t}\sum_{s}|U_{t+1,s+1}|\frac{a_{t}}{a_{s}},$ (89) $\displaystyle\|\mathbb{U}^{\prime}\|={}$ $\displaystyle\sup_{t}\sum_{s}|U^{\prime}_{t+1,s+1}|\frac{a_{t}}{a_{s}}.$ (90) $$

Assuming $\|\mathbb{U}^{\prime}\|<1$, we get

$$ $\|\mathbf{L}^{\prime}\|=\tfrac{1}{2}\|(1-\mathbb{U}^{\prime})\mathbf{V}^{ \prime}\|\leq\frac{\|\mathbf{V}^{\prime}\|}{2(1-\|\mathbb{U}^{\prime}\|)}$ (91) $$

and then

$$ $\|\mathbf{L}\|=\|\tfrac{1}{2}\mathbf{V}+\mathbb{U}\mathbf{L}^{\prime}\|\leq \tfrac{1}{2}\|\mathbf{V}\|+\|\mathbb{U}\|\|\mathbf{L}^{\prime}\|.$ (92) $$

For example, choosing $a_{t}=(t+1)^{\xi},t=0,1,\ldots,$ with a constant exponent $\xi\geq 0$, we obtain a simple sufficient condition for the bound $L_{t}=O(t^{-\xi})$:

###### Proposition 4 .

Suppose that

$$ $\displaystyle|V_{t}|\leq{}$ $\displaystyle C_{V}t^{-\xi},\quad t=1,2,\ldots$ (93) $\displaystyle|V^{\prime}_{t}|\leq{}$ $\displaystyle C_{V^{\prime}}t^{-\xi},\quad t=1,2,\ldots$ (94) $\displaystyle C_{U}={}$ $\displaystyle\sup_{t\geq 1}\sum_{s=1}^{t-1}|U_{t,s}|\Big{(}\frac{t}{s}\Big{)}^ {\xi}<\infty$ (95) $\displaystyle C_{U^{\prime}}={}$ $\displaystyle\sup_{t\geq 1}\sum_{s=1}^{t-1}|U^{\prime}_{t,s}|\Big{(}\frac{t}{s }\Big{)}^{\xi}<1$ (96) $$

with some constants $C_{V},C_{V^{\prime}},C_{U},C_{U^{\prime}}$. Then

$$ $L_{t}\leq C_{L}(t+1)^{-\xi},\quad C_{L}=\tfrac{1}{2}C_{V}+C_{U}\frac{C_{V^{ \prime}}}{2(1-C_{U^{\prime}})}.$ (97) $$

In particular, by taking $\xi=0$ we see that if $\sup_{t\geq 1}\sum_{s=1}^{t-1}|U^{\prime}_{t,s}|<1$ and $\sup_{t\geq 1}\sum_{s=1}^{t-1}|U_{t,s}|<\infty$, then uniform boundedness of $V_{t},V_{t}^{\prime}$ implies uniform boundedness of $L_{t}.$

###### Proposition 4 .

### D.2 Stationary case

The stationary version of loss expansion (73) reads

$$ $\displaystyle L_{t}={}$ $\displaystyle\frac{1}{2}\Big{(}V_{t+1}+\sum_{m=1}^{t}\sum_{0<t_{1}<\ldots<t_{m }<t+1}U_{t+1-t_{m}}U^{\prime}_{t_{m}-t_{m-1}}U^{\prime}_{t_{m-1}-t_{m-2}} \cdots U^{\prime}_{t_{2}-t_{1}}V^{\prime}_{t_{1}}\Big{)}.$ (98) $$

It is will be occasionally convenient to slightly adjust the operator formalism of previous section by considering two-sided sequences. We adjust the definitions (78) and (82) of the vectors $\mathbf{L}$ and $\mathbf{V}$ by padding them by 0 to the left:

$$ $\displaystyle\mathbf{L}={}$ $\displaystyle(\ldots,0,L_{0},L_{1},\ldots),$ (99) $\displaystyle\mathbf{V}={}$ $\displaystyle(\ldots,0,V_{1},V_{2},\ldots).$ (100) $$

Similar convention applies to $\mathbf{L}^{\prime},\mathbf{V}^{\prime}$.

We extend the definition (83) of operator $\mathbb{U}$ to double-sided sequences, also taking into account that due to the stationarity $\mathbb{U}$ is now translation equivariant:

$$ $(\mathbb{U})_{ts}=\begin{cases}U_{t-s},&t>s\\ 0,&t\leq s\end{cases}\equiv(\mathbb{U})_{t-s}.$ (101) $$

Similar convention applies to $\mathbb{U}^{\prime}.$

We now establish a useful general result showing that $U_{t}=O(t^{-\xi_{U}})$ with $\xi_{U}>1$ implies a similar bound for the matrix elements of $(1-\mathbb{U})^{-1}$. This can be easily done if we assume not just $\Sigma_{U}<1$ and $U_{t}=O(t^{-\xi_{U}})$ as in Theorem 2, but rather $U_{t}\leq C^{\prime}t^{-\xi_{U}}$ with a sufficiently small $C^{\prime}$:

###### Lemma 2 .

Suppose that $U_{t}=0$ for $t\leq 0$ and $0\leq U_{t}\leq C^{\prime}t^{-\xi_{U}}$ for $t>0$ with some $C^{\prime}$ and $\xi_{U}>1$. Then, if $C^{\prime}$ is sufficiently small, we have $(\mathbb{U}^{2})_{t}\leq\tfrac{1}{2}C^{\prime}t^{-\xi_{U}}$ for all $t>0$. By iterating this bound, we get

$$ $((1-\mathbb{U})^{-1})_{t}\begin{cases}\leq 2C^{\prime}t^{-\xi_{U}},&t>0\\ =1,&t=0\\ =0,&t<0\end{cases}$ (102) $$

###### Proof.

We have for $t>0$

$$ $\displaystyle(\mathbb{U}^{2})_{t}={}$ $\displaystyle\sum_{r=1}^{t-1}U_{t-r}U_{r}$ (103) $\displaystyle\leq{}$ $\displaystyle 2\sum_{r=1}^{t/2}C^{\prime}(t/2)^{-\xi_{U}}U_{r}$ (104) $\displaystyle\leq{}$ $\displaystyle C^{\prime}t^{-\xi_{U}}2^{1+\xi_{U}}\sum_{r=1}^{t/2}U_{r}$ (105) $\displaystyle\leq{}$ $\displaystyle\Big{[}2^{1+\xi_{U}}C^{\prime}\sum_{r=1}^{\infty}r^{-\xi_{U}}\Big {]}C^{\prime}t^{-\xi_{U}}.$ (106) $$

The factor $2^{1+\xi_{U}}C^{\prime}\sum_{r=1}^{\infty}r^{-\xi_{U}}$ can be made less than $1/2$ by taking $C^{\prime}$ small enough.

The relation (102) follows by expanding $(1-\mathbb{U})^{-1}=1+\mathbb{U}+\mathbb{U}^{2}+\ldots$ and repeatedly applying $\sum_{r=1}^{t-1}[C^{\prime}(t-r)^{-\xi_{U}}][C^{\prime}r^{-\xi_{U}}]\leq\tfrac
{1}{2}C^{\prime}(t-r)^{-\xi_{U}}.$
∎

Now we prove a more subtle statement that only uses the weaker assumptions made in Theorem 2.

###### Lemma 3 .

Suppose that $\Sigma_{U}<1$ and $U_{t}=O(t^{-\xi})$ for some $\xi>1$. Then $((1-\mathbb{U})^{-1})_{t}=O(t^{-\xi}).$

###### Proof.

The operator $(1-\mathbb{U})^{-1}$ is the operator of convolution with the kernel $((1-U)^{-1})_{t}$ which can be expanded as

$$ $((1-U)^{-1})_{t}=\delta_{t}+U_{t}+(U*U)_{t}+\ldots,$ (107) $$

where $*$ denotes the convolution of two sequences, i.e.

$$ $(A*B)_{t}=\sum_{r=0}^{t}A_{t-r}B_{r},\quad r=0,1,\ldots$ (108) $$

for any one-sided sequences $A,B$.
Let us split the sequence $U_{t}$ into a finite initial segment $A_{t}$ and a remainder $B_{t}$ with a small sum:

$$ $\displaystyle U={}$ $\displaystyle A+B,$ (109) $\displaystyle A_{t}={}$ $\displaystyle U_{t}\mathbf{1}_{t\leq t_{0}},$ (110) $\displaystyle B_{t}={}$ $\displaystyle U_{t}\mathbf{1}_{t>t_{0}},$ (111) $$

where $t_{0}$ is such that

$$ $\sum_{t}B_{t}<\epsilon$ (113) $$

with a sufficiently small $\epsilon$ to be chosen later. We then have

$$ $(1-U)^{-1}=\sum_{n=0}^{\infty}U^{*n}=\sum_{n=0}^{\infty}\sum_{k=0}^{n}\binom{n }{k}A^{*k}*B^{*(n-k)},$ (114) $$

where the binomial expansion is legitimate due to the commutativity and associativity of convolution.
Consider the convolutional power $B^{*r}$. Recall that $B_{t}=O(t^{-\xi})$. We can now argue as in Lemma 2 and show by induction that for any $r=1,2,\ldots$

$$ $(B^{*r})_{t}\leq C\epsilon_{1}^{r-1}t^{-\xi}$ (115) $$

with $\epsilon_{1}$ that can be made arbitrarily small by making $\epsilon$ in Eq. (113) small enough. Indeed, for $r\geq 2$

$$ $\displaystyle(B^{*r})_{t}={}$ $\displaystyle\sum_{s=1}^{t-1}(B^{*(r-1)})_{s}B_{t-s}$ (116) $\displaystyle={}$ $\displaystyle\sum_{s=1}^{t/2}(B^{*(r-1)})_{s}B_{t-s}+\sum_{s=t/2+1}^{t-1}(B^{* (r-1)})_{s}B_{t-s}$ (117) $\displaystyle\leq{}$ $\displaystyle\Big{(}\sum_{s=1}^{\infty}(B^{*(r-1)})_{s}\Big{)}C(t/2)^{-\xi}+C \epsilon_{1}^{r-2}(t/2)^{-\xi}\sum_{s=1}^{\infty}B_{s}$ (118) $\displaystyle\leq{}$ $\displaystyle\epsilon^{r-1}C(t/2)^{-\xi}+C\epsilon_{1}^{r-2}(t/2)^{-\xi}\epsilon$ (119) $\displaystyle={}$ $\displaystyle C\epsilon_{1}^{r-1}t^{-\xi}((\epsilon/\epsilon_{1})^{r-1}2^{\xi} +(\epsilon/\epsilon_{1})2^{\xi})$ (120) $\displaystyle\leq{}$ $\displaystyle C\epsilon_{1}^{r-1}t^{-\xi}$ (121) $$

by choosing $\epsilon$ so that $(\epsilon/\epsilon_{1})^{r-1}2^{\xi}+(\epsilon/\epsilon_{1})2^{\xi}<1$ for $r\geq 2$. This proves (115).

Now let

$$ $a=\sum_{t}A_{t}.$ (122) $$

By assumption, $a\leq\sum_{t}U_{t}<1$. Observe that $(A^{*k})_{t}=0$ for $t>kt_{0}$. It follows that

$$ $\displaystyle(A^{*k}*B^{*(n-k)})_{t}\leq{}$ $\displaystyle\Big{(}\sum_{s=1}^{kt_{0}}(A^{*k})_{s}\Big{)}C\epsilon_{1}^{n-k-1 }\max(1,t-kt_{0})^{-\xi}$ (123) $\displaystyle\leq{}$ $\displaystyle a^{k}C\epsilon_{1}^{n-k-1}\max(1,t-kt_{0})^{-\xi}.$ (124) $$

Then from the binomial expansion (114)

$$ $\displaystyle((1-\mathbb{U})^{-1}-1)_{t}\leq{}$ $\displaystyle\sum_{n=1}^{\infty}\Big{(}\sum_{k=0}^{n-1}\binom{n}{k}a^{k}C \epsilon_{1}^{n-1-k}\max(1,t-kt_{0})^{-\xi}+a^{n}\mathbf{1}_{t\leq nt_{0}}\Big {)}$ (125) $\displaystyle={}$ $\displaystyle\sum_{n=1}^{t/(2t_{0})}+\sum_{n=t/(2t_{0})+1}^{\infty}$ (126) $\displaystyle\leq{}$ $\displaystyle(t/2)^{-\xi}\sum_{n=1}^{t/(2t_{0})}\sum_{k=0}^{n}\binom{n}{k}a^{k }C\epsilon_{1}^{n-1-k}+\sum_{n=t/(2t_{0})+1}^{\infty}\sum_{k=0}^{n}\binom{n}{k }a^{k}C\epsilon_{1}^{n-1-k}$ (127) $\displaystyle={}$ $\displaystyle(t/2)^{-\xi}C\sum_{n=1}^{t/(2t_{0})}(a+\epsilon_{1})^{n}+C\sum_{n =t/(2t_{0})+1}^{\infty}(a+\epsilon_{1})^{n}$ (128) $\displaystyle\leq{}$ $\displaystyle C(1-(a+\epsilon_{1}))^{-1}((t/2)^{-\xi}+(a+\epsilon_{1})^{t/(2t_ {0})+1}).$ (129) $$

By choosing $\epsilon_{1}$ sufficiently small so that $a+\epsilon_{1}<1$ we get the desired property $((1-\mathbb{U})^{-1})_{t}=O(t^{-\xi}).$
∎

###### Lemma 2 .

###### Proof.

###### Lemma 3 .

###### Proof.

## Appendix E Proof of Theorem 2

Throughout this section we will use the stationary-case representation of $L_{t}$ from Eq. (19), with the identified $U_{t}=U_{t}^{\prime},V_{t}=V_{t}^{\prime}$:

$$ $\displaystyle L_{t}={}$ $\displaystyle\frac{1}{2}\Big{(}V_{t+1}+\sum_{m=1}^{t}\sum_{0<t_{1}<\ldots<t_{m }<t+1}U_{t+1-t_{m}}U_{t_{m}-t_{m-1}}U_{t_{m-1}-t_{m-2}}\cdots U_{t_{2}-t_{1}}V _{t_{1}}\Big{)}.$ (130) $$

We use the operator formalism and results from Section D.

### E.1 Part 1 (convergence)

If $U_{\Sigma}<1,$ then $V_{t}=O(1)$ implies $L_{t}=O(1)$ by the remark after Proposition 4.

To conclude $L_{t}=o(1)$ from $V_{t}=o(1),$ write

$$ $L_{t}=\frac{1}{2}\sum_{s=1}^{t}((1-\mathbb{U})^{-1})_{s}V_{t-s+1}.$ (131) $$

Here $((1-\mathbb{U})^{-1})_{s}\geq 0$ and

$$ $\sum_{s=1}^{\infty}((1-\mathbb{U})^{-1})_{s}=\frac{1}{1-U_{\Sigma}}<\infty.$ (132) $$

For any $\epsilon>0$ we can choose $r$ such that $\sum_{s=r}^{\infty}((1-\mathbb{U})^{-1})_{s}<\epsilon.$ Then by Eq. (131)

$$ $L_{t}\leq\frac{\max_{s=1,\ldots,r-1}V_{t-s+1}}{2(1-U_{\Sigma})}+\frac{\epsilon \max_{s}V_{s}}{2}\stackrel{{\scriptstyle t\to\infty}}{{\longrightarrow}}\frac{ \epsilon\max_{s}V_{s}}{2}.$ (133) $$

It follows that $L_{t}\to 0.$

### E.2 Part 2 (divergence)

Since $U_{t},V_{t}\geq 0$, we can lower bound, for any $k\geq 0$ and $s\geq 0$,

$$ $L_{t+s}\geq(\mathbb{U}^{k})_{t}V_{s+1}.$ (134) $$

We choose $s$ so that $V_{s+1}>0$. Since $U_{\Sigma}>1,$ we can choose $t_{0}$ such that

$$ $\sum_{t=1}^{t_{0}}U_{t}=a>1.$ (135) $$

Consider the respective truncated sequence $A_{t}$:

$$ $A_{t}=\begin{cases}U_{t}.&t\leq t_{0}\\ 0,&t>t_{0}.\end{cases}$ (136) $$

It follows from Eq. (134) that

$$ $L_{t+s}\geq(A^{*k})_{t}V_{s+1}.$ (137) $$

Note that $(A^{*k})_{t}=0$ for $t>kt_{0}$ and $\sum_{t}(A^{*k})_{t}=(\sum_{t}A_{t})^{k}=a^{k}.$ It follows that

$$ $\sum_{t=s}^{kt_{0}+s}L_{t}\geq a^{k}V_{s+1}.$ (138) $$

Hence

$$ $\sup_{t=0,1,\ldots}L_{t}\geq\sup_{k\geq 0}\frac{a^{k}V_{s+1}}{kt_{0}+s}=\infty.$ (139) $$

### E.3 Part 3 (signal-dominated regime)

The idea of the proof is that in the signal dominated case we expect the main contribution to $L_{t}$ to come from the terms with large time $t_{1}$ of the $V$-factor and small time increments of the other, $U$-factors. We then define the potential leading term $L_{t}^{(V)}$ by first replacing $V_{t_{1}}$ by $V_{t+1}$ in Eq. (130),

$$ $\displaystyle L_{t}^{(V,0)}={}$ $\displaystyle\frac{1}{2}\Big{(}\sum_{m=0}^{t}\sum_{0<t_{1}<\ldots<t_{m}<t+1}U_ {T+1-t_{m}}U_{t_{m}-t_{m-1}}U_{t_{m-1}-t_{m-2}}\cdots U_{t_{2}-t_{1}}\Big{)}V_ {t+1},$ (140) $$

and then extending summation to negative $t$:

$$ $\displaystyle L_{t}^{(V)}={}$ $\displaystyle\frac{1}{2}\Big{(}\sum_{m=0}^{\infty}\sum_{-\infty<t_{1}<\ldots<t _{m}<t+1}U_{t+1-t_{m}}U_{t_{m}-t_{m-1}}U_{t_{m-1}-t_{m-2}}\cdots U_{t_{2}-t_{1 }}\Big{)}V_{t+1}$ (141) $\displaystyle={}$ $\displaystyle\frac{1}{2}\Big{(}\sum_{m=0}^{\infty}\sum_{s_{1},\ldots,s_{m}=1}^ {\infty}U_{s_{m}}U_{s_{m-1}}\cdots U_{s_{1}}\Big{)}V_{t+1}$ (142) $\displaystyle={}$ $\displaystyle\frac{V_{t+1}}{2(1-U_{\Sigma})}$ (143) $\displaystyle={}$ $\displaystyle\frac{C_{V}}{2(1-U_{\Sigma})}(t+1)^{-\xi_{V}}(1+o(1)).$ (144) $$

Thus, to establish the desired asymptotics , it is sufficient to show that the differences $L_{t}^{(V)}-L_{t}^{(V,0)}$ and $L_{t}^{(V,0)}-L_{t}$ are sub-leading compared to the leading term $L_{t}^{(V)}$.

#### 

#### The difference L t ( V ) − L t ( V , 0 ) superscript subscript 𝐿 𝑡 𝑉 superscript subscript 𝐿 𝑡 𝑉 0 L_{t}^{(V)}-L_{t}^{(V,0)} italic_L start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT start_POSTSUPERSCRIPT ( italic_V ) end_POSTSUPERSCRIPT - italic_L start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT start_POSTSUPERSCRIPT ( italic_V , 0 ) end_POSTSUPERSCRIPT .

In terms of the operator $\mathbb{U}$ we can write the vector representing the difference $L_{t}^{(V)}-L_{t}^{(V,0)}$ as

$$ $\mathbf{L}^{(V)}-\mathbf{L}^{(V,0)}=\tfrac{1}{2}\mathbf{V}\odot(1-\mathbb{U})^ {-1}\mathbf{1}_{t<0},$ (145) $$

where $\mathbf{1}_{t<0}(t)=0$ if $t\geq 0$ and 1 otherwise, and $\odot$ is the pointwise product. Consider the norm (88) defined with

$$ $a_{t}=\begin{cases}1,&t\leq 0\\ (t+1)^{\xi},&t\geq 0\end{cases}$ (146) $$

with some $\xi>0$ to be chosen later. Then if we check that $\|\mathbb{U}\|<1,$ then $((1-\mathbb{U})^{-1}\mathbf{1}_{t<0})_{t}=O(t^{-\xi})$ as $t\to+\infty$, in particular implying that

$$ $L_{t}^{(V)}-L_{t}^{(V,0)}=V_{t}\cdot o(1)=o(L_{t}^{(V)}),$ (147) $$

as desired. We have

$$ $\displaystyle\|\mathbb{U}\|={}$ $\displaystyle\sup_{t}\sum_{s}|U_{t-s}|\frac{|a_{t}|}{|a_{s}|}$ (148) $\displaystyle\leq{}$ $\displaystyle\sum_{r>0}U_{r}\sup_{t}\frac{a_{t}}{a_{t-r}}$ (149) $\displaystyle\leq{}$ $\displaystyle\sum_{r>0}U_{r}(r+1)^{\xi}.$ (150) $$

By assumption, $U_{r}=O(r^{-\xi_{U}})$ with some $\xi_{U}>1.$ Then, by the dominated convergence theorem,

$$ $\lim_{\xi\searrow 0}\sum_{r>0}U_{r}(r+1)^{\xi}=\sum_{r>0}U_{r}<1,$ (151) $$

by the assumption $U_{\Sigma}<1$. Therefore, by choosing $\xi$ small enough, we can ensure $\|\mathbb{U}\|<1,$ as desired.

#### The difference L t ( V , 0 ) − L t superscript subscript 𝐿 𝑡 𝑉 0 subscript 𝐿 𝑡 L_{t}^{(V,0)}-L_{t} italic_L start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT start_POSTSUPERSCRIPT ( italic_V , 0 ) end_POSTSUPERSCRIPT - italic_L start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT .

We start by writing

$$ $L_{t}-L_{t}^{(V,0)}=\frac{1}{2}\sum_{m=1}^{t}\sum_{0<t_{1}<\ldots<t_{m}<t+1}U_ {t+1-t_{m}}U_{t_{m}-t_{m-1}}\cdots U_{t_{2}-t_{1}}(V_{t_{1}}-V_{t+1})$ (152) $$

or equivalently

$$ $\mathbf{L}-\mathbf{L}_{V}^{(V,0)}=\tfrac{1}{2}(1-\mathbb{U})^{-1}\mathbf{V}- \tfrac{1}{2}\mathbf{V}\odot(1-\mathbb{U})^{-1}\mathbf{1}_{t\geq 0}.$ (153) $$

Let $k_{t}$ be any integer sequence such that $k_{t}\to\infty$ and $k_{t}=o(t)$. Using $V_{t}=C_{V}t^{-\xi_{V}}(1+o(1))$ and Lemma 3, we have

$$ $\displaystyle 2|L_{t}^{(V,0)}-L_{t}|\leq{}$ $\displaystyle\sum_{s=1}^{t}((1-\mathbb{U})^{-1})_{s}|V_{t+1-s}-V_{t+1}|$ (154) $\displaystyle={}$ $\displaystyle\sum_{s=1}^{k_{t}}+\sum_{s=k_{t}+1}^{t}$ (155) $\displaystyle\leq{}$ $\displaystyle\sum_{s=1}^{k_{t}}((1-\mathbb{U})^{-1})_{s}o(t^{-\xi_{V}})+\sum_{ s=k_{t}+1}^{t}O(s^{-\xi_{U}})O((1+t-s)^{-\xi_{V}})$ (156) $\displaystyle={}$ $\displaystyle o(t^{-\xi_{V}})+\sum_{s=k_{t}}^{t/2}O(s^{-\xi_{U}})O((t-s)^{-\xi _{V}})+\sum_{s=t/2+1}^{t-1}O(s^{-\xi_{U}})O((t-s)^{-\xi_{V}})$ (157) $\displaystyle={}$ $\displaystyle o(t^{-\xi_{V}})+o(1)O(t^{-\xi_{V}})+O(t^{-\xi_{U}})\sum_{s=t/2+1 }^{t-1}O((t-s)^{-\xi_{V}})$ (158) $\displaystyle={}$ $\displaystyle o(t^{-\xi_{V}})+O(t^{-\xi_{U}})\begin{cases}O(1),&\xi_{V}>1\\ O(\ln t),&\xi_{V}=1\\ O(t^{1-\xi_{V}}),&\xi_{V}<1\end{cases}$ (159) $\displaystyle={}$ $\displaystyle o(t^{-\xi_{V}})+\begin{cases}O(t^{-\xi_{U}}),&\xi_{V}>1\\ O(t^{-\xi_{U}}\ln t),&\xi_{V}=1\\ O(t^{(1-\xi_{U})-\xi_{V}}),&\xi_{V}<1\end{cases}$ (160) $\displaystyle={}$ $\displaystyle o(t^{-\xi_{V}}),$ (161) $$

where we used the assumptions $\xi_{V}<\xi_{U}$ and $\xi_{U}>1.$

Having established Eqs. (147), (161), we have completed the proof that

$$ $L_{t}=L_{t}^{(V)}(1+o(1))=\frac{C_{V}}{2(1-U_{\Sigma})}t^{-\xi_{V}}(1+o(1)).$ (162) $$

### E.4 Part 4 (noise-dominated regime)

We transform the original loss $L_{T}$ to an approximation $L_{T}^{(U)}$ through the following three steps:

1. Isolate in $L_{T}$ a long factor $U_{t_{m_{1}}+1-t_{m_{1}}}$ of length $T(1-o(1))$:

$$ $\displaystyle L_{T}^{(U,0)}={}$ $\displaystyle\frac{1}{2}\sum_{\begin{subarray}{c}m_{1}=1\\ m_{2}=0\end{subarray}}^{\infty}\sum_{\begin{subarray}{c}0<t_{1}<\ldots<t_{m_{1 }}<T_{1}(T)\\ T-T_{1}(T)<t_{m_{1}+1}<\ldots<t_{m_{1}+m_{2}}<T+1\end{subarray}}(U_{T+1-t_{m_{ 1}+m_{2}}}\cdots U_{t_{m_{1}+2}-t_{m_{1}+1}})$ (163) $\displaystyle\times U_{t_{m_{1}+1}-t_{m_{1}}}(U_{t_{m_{1}}-t_{m_{1}-1}}\cdots U _{t_{2}-t_{1}})V_{t_{1}}$ (164) $\displaystyle={}$ $\displaystyle\frac{1}{2}\sum_{\begin{subarray}{c}0<t_{1}<t<T_{1}(T)\\ T-T_{1}(T)<s<T+1\end{subarray}}((1-\mathbb{U})^{-1})_{T+1-s}U_{s-t}((1-\mathbb {U})^{-1})_{t-t_{1}}V_{t_{1}},$ (165) $$

where $T_{1}(T)=T^{\tfrac{\xi_{U}+\xi_{V}}{2\xi_{V}}}$ so that in particular $T_{1}(T)\to\infty$ and $T_{1}(T)=o(T),$ due to $\xi_{V}>\xi_{U}$.

2. Replace in $L_{T}^{(U,0)}$ the long factor $U_{t_{m_{1}+1}-t_{m_{1}}}$ by $U_{T}$:

$$ $\displaystyle L_{T}^{(U,1)}={}$ $\displaystyle\frac{1}{2}\sum_{\begin{subarray}{c}0<t_{1}<t<T_{1}(T)\\ T-T_{1}(T)<s<T+1\end{subarray}}((1-\mathbb{U})^{-1})_{T+1-s}U_{T}((1-\mathbb{U })^{-1})_{t-t_{1}}V_{t_{1}}.$ (166) $$

3. Finally, extend summation in $L_{T}^{(U,1)}$, dropping the upper constraints on $t_{1},t$ and the lower constraints on $s$:

$$ $\displaystyle L_{T}^{(U)}={}$ $\displaystyle\frac{1}{2}\sum_{\begin{subarray}{c}0<t_{1}<t<\infty\\ -\infty<s<T+1\end{subarray}}((1-\mathbb{U})^{-1})_{T+1-s}U_{T}((1-\mathbb{U})^ {-1})_{t-t_{1}}V_{t_{1}}$ (167) $\displaystyle={}$ $\displaystyle\frac{V_{\Sigma}U_{T}}{2(1-U_{\Sigma})^{2}}$ (168) $\displaystyle={}$ $\displaystyle\frac{V_{\Sigma}C_{U}}{2(1-U_{\Sigma})^{2}}(T+1)^{-\xi_{U}}(1+o(1 )).$ (169) $$

We argue now that each of these approximations introduces an error of order $o(T^{-\xi_{U}})$.

- 1.
$[L\to L^{(U,0)}].$ Note that the difference $L_{T}-L_{T}^{(U,0)}$ consists of the terms where either $t_{1}>T_{1}(T)$, or there exists $k>1$ for which $T_{1}(T)<t_{k}<T-T_{1}(T)$. Let us bound the contribution of such terms.
In the first case, the contribution is bounded by $O((T_{1}(T))^{-\xi_{V}})=O(T^{-\tfrac{\xi_{U}+\xi_{V}}{2}})=o(T^{-\xi_{U}})$, since $\xi_{V}>\xi_{U}$.
In the second case using $((1-\mathbb{U})^{-1})_{t}=O(t^{-\xi_{U}})$ and $V_{t}=O(t^{-\xi_{U}})$, the contribution of these terms can be upper bounded by
$\displaystyle\sum_{t_{k}=T_{1}(T)}^{T-T_{1}(T)}O((T-t_{k})^{-\xi_{U}})O(t_{k}^
{-\xi_{U}})={}$
$\displaystyle O(T^{-\xi_{U}})\sum_{t_{k}=T_{1}(T)}^{T/2}O(t_{k}^{-\xi_{U}})$
(170)
$\displaystyle={}$
$\displaystyle O(T^{-\xi_{U}})O((T_{1}(T))^{1-\xi_{U}}-T^{1-\xi_{U}})$
(171)
$\displaystyle={}$
$\displaystyle o(T^{-\xi_{U}}),$
(172)
since $\xi_{U}>1$.
- 2.
$[L^{(U,0)}\to L^{(U,1)}].$ Since $s-t=T(1+o(1)),$ we have $U_{s-t}-U_{T}=o(T^{-\xi_{U}})$ and then also $|L^{(U,0)}_{T}-L_{T}^{(U,1)}|=o(T^{-\xi_{U}})$ by the boundedness of the summed products of the remaining factors.
- 3.
$[L^{(U,1)}\to L^{(U)}].$ Here one can use similar arguments, taking into account the relation $U_{T}=O(T^{-\xi_{U}})$ and the tail bound
$\sum_{t>T_{1}}((1-\mathbb{U})^{-1})_{t}=O(T_{1}^{1-\xi_{U}}),$
(173)
following from Lemma 3.

This completes the proof of this part of Theorem 2.

## Appendix F Generalization of Theorem 2 to V t ≠ V t ′ , U t ≠ U t ′ formulae-sequence subscript 𝑉 𝑡 superscript subscript 𝑉 𝑡 ′ subscript 𝑈 𝑡 superscript subscript 𝑈 𝑡 ′ V_{t}\neq V_{t}^{\prime},U_{t}\neq U_{t}^{\prime} italic_V start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT ≠ italic_V start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT , italic_U start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT ≠ italic_U start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT

In this section we drop the assumption $V_{t}=V_{t}^{\prime},U_{t}=U_{t}^{\prime}$ adopted for Theorem 2. Accordingly, $V_{t},U_{t},L_{t}$ are now distinct from $V^{\prime}_{t},U^{\prime}_{t},L^{\prime}_{t}$. We define the sums $U^{\prime}_{\Sigma}=\sum tU_{t}^{\prime},V^{\prime}_{\Sigma}=\sum_{t}V_{t}^{\prime}$ similarly to the sums $U_{\Sigma},V_{\Sigma}$. The respective generalized version of Theorem 2 reads

###### Theorem 7 .

Suppose that $L_{t}$ is an expansion (19) with some $U_{t},U_{t}^{\prime},V_{t},V_{t}^{\prime}\geq 0.$

- 1.
[Convergence] Let $U^{\prime}_{\Sigma}<1$ and $U_{\Sigma}<\infty$. At $t\to\infty$, if $V_{t},V^{\prime}_{t}=O(1)$ (respectively, $V_{t},V^{\prime}_{t}=o(1)$), then also $L_{t}=O(1)$ (respectively, $L_{t}=o(1)$).
- 2.
[Divergence] If $U^{\prime}_{\Sigma}>1$, $V^{\prime}_{t}>0$ at least for one $t$ and $U_{t}>0$ at least for one $t$, then $\sup_{t=1,2,\ldots}L_{t}=\infty.$
- 3.
[Signal-dominated regime] Suppose that there exist constants $\xi_{V},C_{V}>0$ such that $V_{t},V^{\prime}_{t}=C_{V}t^{-\xi_{V}}(1+o(1))$ as $t\to\infty$. Suppose also that $U^{\prime}_{\Sigma}<1$ and $U_{t},U^{\prime}_{t}=O(t^{-\xi_{U}})$ with some $\xi_{U}>\max(\xi_{V},1).$ Then
$L_{t}=\frac{(1-U^{\prime}_{\Sigma}+U_{\Sigma})C_{V}}{2(1-U^{\prime}_{\Sigma})}
t^{-\xi_{V}}(1+o(1)).$
(174)
- 4.
[Noise-dominated regime] Suppose that there exist constants $\xi_{V}>\xi_{U}>1,C_{U}>0$ such that $U_{t},U_{t}^{\prime}=C_{U}t^{-\xi_{U}}(1+o(1))$ and $V_{t},V_{t}^{\prime}=O(t^{-\xi_{V}})$ as $t\to\infty$. Suppose also that $U^{\prime}_{\Sigma}<1$. Then
$L_{t}=\frac{(1-U^{\prime}_{\Sigma}+U_{\Sigma})V^{\prime}_{\Sigma}C_{U}}{2(1-U^
{\prime}_{\Sigma})^{2}}t^{-\xi_{U}}(1+o(1)).$
(175)

Note that the only difference in the conclusions of Theorems 2 and 7 is the factor $1-U^{\prime}_{\Sigma}+U_{\Sigma}$ appearing in the leading terms. In the case $U^{\prime}_{\Sigma}=U_{\Sigma}$ this factor evaluates to 1.

We remark also that theorem 7 assumes in item 3 that the quantities $V_{t}$ and $V^{\prime}_{t}$, though generally different, have the same leading terms $C_{V}t^{-\xi_{V}};$ in particular, with the same constant $C_{V}$ rather than two different constants $C_{V},C_{V^{\prime}}$. This assumption is not important and is made simply because it is fulfilled by the actual $V_{t},V^{\prime}_{t}$ as analyzed in Theorem 4. The same remark applies to $U_{t},U_{t}$ in item 4.

###### Proof of Theorem 7 .

1. [Convergence] The $O(1)$ part of the claim follows from the remark after Proposition 4.

For the $o(1)$ part, we already know from Theorem 2 that $U^{\prime}_{\Sigma}<1$ and $V_{t}^{\prime}=o(1)$ imply $L_{t}^{\prime}=o(1)$. Using the relation

$$ $L_{t}=\frac{1}{2}\Big{(}V_{t+1}+\sum_{s=1}^{t}U_{s}L^{\prime}_{t-s}\Big{)}$ (176) $$

and the conditions $V_{t}=o(1),$ $U_{\Sigma}<\infty$, we then get $L_{t}=o(1)$.

2. [Divergence] We already know from Theorem 2 that $\sup_{t}L_{t}^{\prime}=\infty$. Since, for any $s>0$, $L_{t}\geq U_{s}L_{t-s}^{\prime}$, and $U_{s}>0$ at least for one $s$, we also have $\sup_{t}L_{t}=\infty$.

3. [Signal-dominated regime] We already know by Theorem 2 that

$$ $L_{t}^{\prime}=\frac{C_{V}}{2(1-U^{\prime}_{\Sigma})}t^{-\xi_{V}}(1+o(1)).$ (177) $$

Let $\gamma=\frac{\xi_{U}+\max(\xi_{V},1)}{2\xi_{U}}\in(0,1)$ and observe that

$$ $\displaystyle L_{t}={}$ $\displaystyle\frac{1}{2}V_{t+1}+\sum_{s=0}^{t-1}U_{t-s}L^{\prime}_{s}$ (178) $\displaystyle={}$ $\displaystyle\frac{C_{V}}{2}t^{-\xi_{V}}(1+o(1))+\sum_{s=0}^{t-t^{\gamma}}U_{t -s}L^{\prime}_{s}+\sum_{s=t-t^{\gamma}+1}^{t}U_{t-s}L^{\prime}_{s}$ (179) $\displaystyle={}$ $\displaystyle\frac{C_{V}}{2}t^{-\xi_{V}}(1+o(1))+O\Big{(}t^{-\gamma\xi_{U}} \sum_{s=1}^{t}s^{-\xi_{V}}\Big{)}+\frac{C_{V}}{2(1-U^{\prime}_{\Sigma})}t^{- \xi_{V}}(1+o(1))\sum_{s=t-t^{\gamma}+1}^{t}U_{t-s}$ (180) $\displaystyle={}$ $\displaystyle\frac{1}{2}\Big{(}C_{V}+\frac{C_{V}U_{\Sigma}}{1-U_{\Sigma}^{ \prime}}\Big{)}t^{-\xi_{V}}(1+o(1))+O(t^{-\frac{\xi_{U}+\max(\xi_{V},1)}{2}}) \begin{cases}O(1),&\xi_{V}>1\\ O(\ln t),&\xi_{V}=1\\ O(t^{1-\xi_{V}}),&\xi_{V}<1\end{cases}$ (181) $\displaystyle={}$ $\displaystyle\frac{1}{2}\Big{(}C_{V}+\frac{C_{V}U_{\Sigma}}{1-U_{\Sigma}^{ \prime}}\Big{)}t^{-\xi_{V}}(1+o(1))$ (182) $$

as claimed.

4. [Noise-dominated regime] We know that, by Theorem 2,

$$ $L_{t}^{\prime}=\frac{V^{\prime}_{\Sigma}C_{U}}{2(1-U^{\prime}_{\Sigma})^{2}}t^ {-\xi_{U}}(1+o(1)).$ (183) $$

We have

$$ $\displaystyle L_{t}={}$ $\displaystyle\frac{1}{2}V_{t+1}+\sum_{s=0}^{t-1}U_{t-s}L^{\prime}_{s}=\sum_{s= 0}^{t-1}U_{t-s}L^{\prime}_{s}+o(t^{-\xi_{U}}).$ (184) $$

Arguing as earlier, the leading contribution to the sum $\sum_{s=0}^{t-1}U_{t-s}L^{\prime}_{s}$ comes from terms with $s$ close to 0 or $t$, since the terms in the middle are $O(t^{-2\xi_{U}})$. It follows that

$$ $L_{t}=\Big{[}C_{U}\sum_{s=0}^{\infty}L_{s}^{\prime}+\frac{V^{\prime}_{\Sigma}C _{U}}{2(1-U^{\prime}_{\Sigma})^{2}}U_{\Sigma}\Big{]}t^{-\xi_{U}}(1+o(1)).$ (185) $$

We can compute

$$ $\sum_{s=0}^{\infty}L_{s}^{\prime}=\frac{1}{2}\sum_{s=0}^{\infty}\sum_{t=1}^{s} (1-\mathbb{U}^{\prime})^{-1}_{s-t}V^{\prime}_{t}=\frac{1}{2}\sum_{s=0}^{\infty }(1-\mathbb{U}^{\prime})^{-1}_{s}\sum_{t=1}^{\infty}V^{\prime}_{t}=\frac{V^{ \prime}_{\Sigma}}{2(1-U^{\prime}_{\Sigma})}.$ (186) $$

Substituting in (185),

$$ $\displaystyle L_{t}={}$ $\displaystyle\Big{[}C_{U}\frac{V^{\prime}_{\Sigma}}{2(1-U^{\prime}_{\Sigma})}+ \frac{V^{\prime}_{\Sigma}C_{U}}{2(1-U^{\prime}_{\Sigma})^{2}}U_{\Sigma}\Big{]} t^{-\xi_{U}}(1+o(1))$ (187) $\displaystyle={}$ $\displaystyle\frac{(1-U^{\prime}_{\Sigma}+U_{\Sigma})V^{\prime}_{\Sigma}C_{U}} {2(1-U^{\prime}_{\Sigma})^{2}}t^{-\xi_{U}}(1+o(1)).$ (188) $$

∎

###### Theorem 7 .

###### Proof of Theorem 7 .

## Appendix G Proof of Theorem 3

Recall that if $B,\Delta B$ are linear operators and $B$ has a simple eigenvalue $\mu_{0}$ associated with a one-dimensional eigenprojector $P$, then for small $\Delta B$ the operator $B+\Delta B$ has a simple eigenvalue $\mu$ given by

$$ $\mu=\mu_{0}+\operatorname{Tr}[P\Delta B]+O(\|\Delta B\|^{2}).$ (189) $$

We apply this to $B=A_{\lambda=0}$ and $\Delta B=A_{\lambda}-A_{\lambda=0}$ with $A_{\lambda}$ given by (18). Observe that the noisy term in Eq. (18) is $O(\lambda^{2}),$ so it will only contribute a correction $O(\lambda^{2})$ to $\mu_{\lambda}$. On the other hand, the main term in $A_{\lambda}$, given by $\mathbf{M}\mapsto S_{\lambda}\mathbf{M}S_{\lambda}^{T}$, acts as a tensor square $S_{\lambda}\otimes S_{\lambda}$. Accordingly, the eigenvalue $\mu_{A,\lambda}$ of $A_{\lambda}$ can be written as

$$ $\mu_{A,\lambda}=z_{\lambda}^{2}+O(\lambda^{2}),$ (190) $$

where $z_{\lambda}$ is the leading eigenvalue of the matrix $S_{\lambda}$. The leading linear term in $z_{\lambda}$ can be found, for example, from the condition

$$ $\chi(z_{\lambda},\lambda)\equiv 0,$ (191) $$

where $\chi(x,\lambda)$ is the characteristic polynomial of $S_{\lambda}$ defined in Eq. (30). Differentiating condition (191) over $\lambda$, we get

$$ $\frac{\partial\chi}{\partial x}\frac{dz_{\lambda}}{d\lambda}+\frac{\partial \chi}{\partial\lambda}=0.$ (192) $$

Using representation (37) of $\chi$ and taking into account that $z_{\lambda=0}=1$, this implies

$$ $\displaystyle\frac{dz_{\lambda}}{d\lambda}(\lambda=0)={}$ $\displaystyle-\Big{(}\frac{\partial\chi}{\partial\lambda}\Big{/}\frac{\partial \chi}{\partial x}\Big{)}(\lambda=0)=\frac{\sum_{m=0}^{M}q_{m}}{\sum_{m=0}^{M} mp_{m}-M-1}=-\alpha_{\mathrm{eff}},$ (193) $$

where $\alpha_{\mathrm{eff}}$ is given by the second equality in (23). To see that it is equivalent to the first equality $\alpha_{\mathrm{eff}}=\alpha-\mathbf{b}^{T}(1-D)^{-1}\mathbf{c}$, note from Eq. (36) that

$$ $\displaystyle\frac{\partial\chi(x,\lambda)}{\partial x}(x=1,\lambda=0)={}$ $\displaystyle\det(1-D),$ (194) $\displaystyle\frac{\partial\chi(x,\lambda)}{\partial\lambda}(x=1,\lambda=0)={}$ $\displaystyle-\det\begin{pmatrix}\mathbf{a}^{T}\mathbf{c}-\alpha&\mathbf{a}^{T }(1-D)-\mathbf{b}^{T}\\ \mathbf{c}&1-D\end{pmatrix}$ (195) $\displaystyle={}$ $\displaystyle\det\begin{pmatrix}\alpha&\mathbf{b}^{T}\\ \mathbf{c}&1-D\end{pmatrix}$ (196) $\displaystyle={}$ $\displaystyle\det(1-D)\det\begin{pmatrix}\alpha&\mathbf{b}^{T}\\ (1-D)^{-1}\mathbf{c}&I\end{pmatrix}$ (197) $\displaystyle={}$ $\displaystyle\det(1-D)(\alpha-\mathbf{b}^{T}(1-D)^{-1}\mathbf{c}).$ (198) $$

We conclude that $z_{\lambda}=1-\alpha_{\mathrm{eff}}\lambda+O(\lambda^{2})$ and hence
$\mu_{A,\lambda}=1-2\alpha_{\mathrm{eff}}\lambda+O(\lambda^{2})$, as claimed.

## Appendix H Proof of Theorem 4

#### Asymptotics of V t , V t ′ subscript 𝑉 𝑡 subscript superscript 𝑉 ′ 𝑡 V_{t},V^{\prime}_{t} italic_V start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT , italic_V start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT .

Recall that

$$ $V^{\prime}_{t}=\sum_{k}\lambda_{k}c_{k}^{2}\langle(\begin{smallmatrix}1\\ \mathbf{a}\end{smallmatrix})(\begin{smallmatrix}1\\ \mathbf{a}\end{smallmatrix})^{T},A^{t-1}_{\lambda_{k}}(\begin{smallmatrix}1&0 \\ 0&0\end{smallmatrix})\rangle.$ (199) $$

Recall by Theorem 3 that for sufficiently small $\lambda>0$, say for $\lambda<\lambda_{*}$, the maps $A_{\lambda}$ have a simple leading eigenvalue

$$ $\mu_{A,\lambda}=1-2\alpha_{\mathrm{eff}}\lambda+O(\lambda^{2}).$ (200) $$

Let $P_{\lambda}:\mathbb{R}^{(M+1)\times(M+1)}\to\mathbb{R}^{(M+1)\times(M+1)}$ denote the corresponding one-dimensional spectral projector. Then, by the assumption of strict stability of $A_{\lambda},$

$$ $V^{\prime}_{t}=\sum_{k:\lambda_{k}<\lambda_{*}}\lambda_{k}c_{k}^{2}\langle( \begin{smallmatrix}1\\ \mathbf{a}\end{smallmatrix})(\begin{smallmatrix}1\\ \mathbf{a}\end{smallmatrix})^{T},P_{\lambda_{k}}(\begin{smallmatrix}1&0\\ 0&0\end{smallmatrix})\rangle\mu^{t-1}_{A,\lambda_{k}}+O(r^{t}),\quad t\to\infty,$ (201) $$

for some $0<r<1.$ The second term will be asymptotically negligible, so we can focus on the first one.

Next, observe that the projector $P_{\lambda}$ continuously depends on $\lambda$. At $\lambda=0$ the matrix $\mathbf{M}_{0}=(\begin{smallmatrix}1&0\\
0&0\end{smallmatrix})$ is precisely the leading eigenvector of $A_{\lambda=0}$, so

$$ $\langle(\begin{smallmatrix}1\\ \mathbf{a}\end{smallmatrix})(\begin{smallmatrix}1\\ \mathbf{a}\end{smallmatrix})^{T},P_{\lambda}(\begin{smallmatrix}1&0\\ 0&0\end{smallmatrix})\rangle=1+o(1),\quad\lambda\to 0.$ (202) $$

Since for any fixed $\widetilde{\lambda}>0$ the contribution of the terms with $\lambda_{k}>\widetilde{\lambda}$ to the expansion (201) is exponentially small, we get

$$ $V^{\prime}_{t}=(1+o(1))\sum_{k:\lambda_{k}<\lambda_{*}}\lambda_{k}c_{k}^{2}\mu ^{t-1}_{A,\lambda_{k}}+O(r^{t}),\quad t\to\infty.$ (203) $$

Given power-law spectral assumptions (1) and relation (200), the asymptotic of $\sum_{k:\lambda_{k}<\lambda_{*}}\lambda_{k}c_{k}^{2}\mu^{t-1}_{A,\lambda_{k}}$ is derived by standard methods (see e.g. Theorem 1 in Appendix B of Velikanov & Yarotsky ()):

$$ $\sum_{k:\lambda_{k}<\lambda_{*}}\lambda_{k}c_{k}^{2}\mu^{t-1}_{A,\lambda_{k}}= (1+o(1))Q\Gamma(\zeta+1)(2\alpha_{\mathrm{eff}}t)^{-\zeta}.$ (204) $$

We conclude that

$$ $V_{t}^{\prime},V_{t}=(1+o(1))Q\Gamma(\zeta+1)(2\alpha_{\mathrm{eff}}t)^{-\zeta},$ (205) $$

as claimed.

#### Asymptotics of U t , U t ′ subscript 𝑈 𝑡 subscript superscript 𝑈 ′ 𝑡 U_{t},U^{\prime}_{t} italic_U start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT , italic_U start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT .

We have

$$ $U_{t}^{\prime}=\frac{\tau_{1}}{|B|}\sum_{k}\lambda_{k}^{2}\langle(\begin{ smallmatrix}1\\ \mathbf{a}\end{smallmatrix})(\begin{smallmatrix}1\\ \mathbf{a}\end{smallmatrix})^{T},A^{t-1}_{\lambda_{k}}(\begin{smallmatrix}- \alpha\\ \mathbf{c}\end{smallmatrix})(\begin{smallmatrix}-\alpha\\ \mathbf{c}\end{smallmatrix})^{T}\rangle.$ (206) $$

We argue as in the previous paragraph, but the difference now is that the matrix $(\begin{smallmatrix}-\alpha\\
\mathbf{c}\end{smallmatrix})(\begin{smallmatrix}-\alpha\\
\mathbf{c}\end{smallmatrix})^{T}$ is not an eigenvector of $A_{\lambda=0},$ so we need to find the projection $P_{\lambda=0}\big{[}(\begin{smallmatrix}-\alpha\\
\mathbf{c}\end{smallmatrix})(\begin{smallmatrix}-\alpha\\
\mathbf{c}\end{smallmatrix})^{T}\big{]}$. At $\lambda=0$, this can be factored as

$$ $P_{\lambda=0}\big{[}(\begin{smallmatrix}-\alpha\\ \mathbf{c}\end{smallmatrix})(\begin{smallmatrix}-\alpha\\ \mathbf{c}\end{smallmatrix})^{T}\big{]}=\big{[}P_{S_{0}}(\begin{smallmatrix}- \alpha\\ \mathbf{c}\end{smallmatrix})\big{]}\big{[}P_{S_{0}}(\begin{smallmatrix}-\alpha \\ \mathbf{c}\end{smallmatrix})\big{]}^{T},$ (207) $$

where $P_{S_{0}}$ is the spectral projector for the matrix $S_{\lambda=0}$ and its eigenvalue 1. The corresponding eigenvector of $S_{\lambda=0}$ is $(\begin{smallmatrix}1\\
0\end{smallmatrix}),$ so

$$ $P_{S_{0}}(\begin{smallmatrix}-\alpha\\ \mathbf{c}\end{smallmatrix})=x(\begin{smallmatrix}1\\ 0\end{smallmatrix})$ (208) $$

with a coefficient $x$ that can be found, for example, from the relation

$$ $A_{\lambda=0}^{t}(\begin{smallmatrix}-\alpha\\ \mathbf{c}\end{smallmatrix})=x(\begin{smallmatrix}1\\ 0\end{smallmatrix})+o(1),\quad t\to\infty.$ (209) $$

We have

$$ $A_{\lambda=0}^{t}=\begin{pmatrix}1&\mathbf{b}^{T}\sum_{s=0}^{t-1}D^{s}\\ 0&D^{t}\end{pmatrix},$ (210) $$

so

$$ $x=\lim_{t\to\infty}\Big{(}-\alpha+\mathbf{b}^{T}\sum_{s=0}^{t-1}D^{s}\mathbf{c }\Big{)}=-\alpha_{\mathrm{eff}}.$ (211) $$

It follows that

$$ $\langle(\begin{smallmatrix}1\\ \mathbf{a}\end{smallmatrix})(\begin{smallmatrix}1\\ \mathbf{a}\end{smallmatrix})^{T},P_{\lambda}\big{[}(\begin{smallmatrix}-\alpha \\ \mathbf{c}\end{smallmatrix})(\begin{smallmatrix}-\alpha\\ \mathbf{c}\end{smallmatrix})^{T}\big{]}\rangle=\alpha^{2}_{\mathrm{eff}}+o(1), \quad\lambda\to 0$ (212) $$

and hence

$$ $U^{\prime}_{t}=(1+o(1))\frac{\tau_{1}}{|B|}\alpha^{2}_{\mathrm{eff}}\sum_{k: \lambda_{k}<\lambda_{*}}\lambda^{2}_{k}\mu^{t-1}_{A,\lambda_{k}}+O(r^{t}), \quad t\to\infty.$ (213) $$

Applying the power law spectral assumption on $\lambda_{k}$ and computing the standard asymptotic, we finally get

$$ $U^{\prime}_{t},U_{t}=(1+o(1))\frac{(\alpha_{\mathrm{eff}}\Lambda)^{1/\nu}\tau_ {1}\Gamma(2-1/\nu)}{|B|\nu}(2t)^{1/\nu-2}.$ (214) $$

## Appendix I Proof of Theorem 5

1. Strict stability of $D$ is equivalent to

$$ $0<\delta<2.$ (215) $$

By Theorem 3, $\mu_{A,\lambda}=1-2\alpha_{\mathrm{eff}}\lambda+O(\lambda),$ so

$$ $\alpha_{\mathrm{eff}}\equiv-\frac{q_{0}+q_{1}}{\delta}\geq 0$ (216) $$

is a necessary condition of strict stability of $S_{\lambda}$ at small $\lambda>0$. We can exclude the case of equality $\alpha_{\mathrm{eff}}=0:$ in this case $q_{0}=-q_{1}$ and so $\chi(1,\lambda)=\det(1-S_{\lambda})=0$ for all $\lambda$, i.e. 1 is an eigenvalue of $S_{\lambda}$ for all $\lambda$.

Another necessary condition of strict stability is that the product of both eigenvalues of $S_{\lambda}$ is less than 1 in absolute value, i.e. $|p_{0}+\lambda q_{0}|<1.$ This holds for all $0\leq\lambda\leq\Lambda$ iff this holds for $\lambda=0$ and $\lambda=\lambda_{\max}$. The condition for $\lambda=0$ is the already established condition $0<\delta<2$, while the condition for $\lambda_{\max}$ means $-1<p_{0}+\lambda q_{0}<1,$ i.e.

$$ $-\frac{\delta}{\lambda_{\max}}<q_{0}<\frac{2-\delta}{\lambda_{\max}}.$ (217) $$

Under condition (217), stability can only fail to hold if $S_{\lambda}$ has real eigenvalues. At the lowest point $\lambda$ of this failure, one of the eigenvalues of $S_{\lambda}$ must be 1 or $-1$. Since $\chi(1,\lambda)=-\lambda(q_{0}+q_{1}),$ 1 cannot be an eigenvalue of $S_{\lambda}$ for $\lambda>0$ if $\alpha_{\mathrm{eff}}>0.$ On the other hand the condition $\chi(-1,\lambda)=0$ reads

$$ $4-2\delta+\lambda(q_{1}-q_{0})=0.$ (218) $$

At $\lambda=0$ the l.h.s. is positive, so $S_{\lambda}$ does not have the eigenvalue $-1$ for $0\leq\lambda\leq\lambda_{\max}$ iff the l.h.s. is also positive at $\lambda=\lambda_{\max}$, i.e.

$$ $\delta\alpha_{\mathrm{eff}}-2q_{0}=q_{0}-q_{1}<\frac{4-2\delta}{\lambda_{\max}}.$ (219) $$

Note that adding this inequality to $\alpha_{\mathrm{eff}}>0$ gives the right inequality (217), so we can drop the latter.

Summarizing, the full set of necessary and sufficient conditions of strict stability is

$$ $0<\delta<2,\quad\alpha_{\mathrm{eff}}>0,\quad-q_{0}<\frac{\delta}{\lambda_{ \max}},\quad\delta\alpha_{\mathrm{eff}}<\frac{4-2\delta}{\lambda_{\max}}-2q_{0}.$ (220) $$

2. The characteristic equation $\det(x-S_{\lambda})=0$ is

$$ $(x-1)(x-1+\delta)-\lambda(q_{1}x+q_{0})=0.$ (221) $$

a) Suppose that $q_{1}\leq-\alpha_{\mathrm{eff}}.$ From the condition $\alpha_{\mathrm{eff}}>0$ we have $q_{0}+q_{1}<0$. We also have $q_{1}\leq-\alpha_{\mathrm{eff}}<0$. It follows that

$$ $-\frac{q_{0}}{q_{1}}<1.$ (222) $$

On the other hand, $q_{1}\leq-\alpha_{\mathrm{eff}}$ means $q_{1}\leq\tfrac{q_{0}+q_{1}}{\delta}$ or, equivalently,

$$ $1-\delta\leq-\frac{q_{0}}{q_{1}}.$ (223) $$

Inequalities (222), (223) show that the root $-\tfrac{q_{0}}{q_{1}}$ of the linear $\lambda$-term $\lambda(q_{1}x+q_{0})$ lies between the roots $1,1-\delta$ of the quadratic term $(x-1)(x-1+\delta)$. It follows that at all $\lambda\in\mathbb{R}$ the matrix $S_{\lambda}$ has two real eigenvalues (possibly coinciding if $-\tfrac{q_{0}}{q_{1}}=1-\delta$).

b) Now suppose that $q_{1}>-\alpha_{\mathrm{eff}}$. Suppose first that $q_{1}>0.$ Then

$$ $-\frac{q_{0}}{q_{1}}>1,$ (224) $$

i.e. the root of the term $\lambda(q_{1}x+q_{0})$ lies to the right of both roots $1,1-\delta$. Then, since $q_{1}+q_{0}<0,$ the polynomial $\chi(x,\lambda)$ has two different real $x$ roots at sufficiently small $\lambda>0.$ Then, as $\lambda$ is increased, the roots become not real and complex conjugate, and then again become real, at even larger $\lambda$.

To find the geometric position of the roots, note that they satisfy the condition

$$ $\frac{(x-1)(x-1+\delta)}{q_{1}x+q_{0}}\in\mathbb{R}.$ (225) $$

This can be written as

$$ $hJ\Big{(}\frac{q_{1}x+q_{0}}{h}\Big{)}\in\mathbb{R},$ (226) $$

where

$$ $h=\sqrt{(q_{0}+q_{1})(q_{0}+q_{1}-\delta q_{1})}=\delta\sqrt{\alpha_{\mathrm{ eff}}(\alpha_{\mathrm{eff}}+q_{1})}$ (227) $$

and $J(z)=z+\tfrac{1}{z}$ is the Joukowsky function. Since $J^{-1}(\mathbb{R})=\mathbb{R}\setminus\{0\}\cup\{z\in\mathbb{C}:|z|=1\},$ the roots $x$ either are real or belong to the circle $|q_{1}x+q_{0}|=h$, as claimed.

If $q_{1}>-\alpha_{\mathrm{eff}}$ but $q_{1}<0,$ we have

$$ $-\frac{q_{0}}{q_{1}}<1-\delta,$ (228) $$

i.e. the root of $q_{1}x+q_{0}$ is to the left of both roots $1,1-\delta.$ The same pattern as above then applies – at small $\lambda>0$ there are two real eigenvalues, then complex eigenvalues, then again real, and the non-real eigenvalues again lie on the circle $|q_{1}x+q_{0}|=h$.

In the special case $q_{1}=0$ the root $-\tfrac{q_{0}}{q_{1}}$ becomes infinite. For sufficiently small $\lambda=0$ the eigenvalues are real, and then for all larger $\lambda$ become complex. The circle $|q_{1}x+q_{0}|=h$ degenerates into the straight line $\Re x=1-\tfrac{\delta}{2}.$

3. We can compute the value $R_{\lambda}$ using the identity

$$ $R_{\lambda}\equiv\sum_{t=0}^{\infty}\langle(\begin{smallmatrix}1\\ \ a\end{smallmatrix})(\begin{smallmatrix}1\\ \ a\end{smallmatrix})^{T},A^{t}_{\lambda}(\begin{smallmatrix}-\alpha\\ \ c\end{smallmatrix})(\begin{smallmatrix}-\alpha\\ \ c\end{smallmatrix})^{T}\rangle=\langle(\begin{smallmatrix}1\\ \ a\end{smallmatrix})(\begin{smallmatrix}1\\ \ a\end{smallmatrix})^{T},(1-A_{\lambda})^{-1}(\begin{smallmatrix}-\alpha\\ \ c\end{smallmatrix})(\begin{smallmatrix}-\alpha\\ \ c\end{smallmatrix})^{T}\rangle.$ (229) $$

Here, $A_{\lambda}$ can be viewed as a linear operator on $\mathbb{R}^{2\times 2}$ or, simpler, on the space of symmetric matrices $\mathbb{R}_{\mathrm{sym}}^{2\times 2}\cong\mathbb{R}^{3}.$ The expression $(1-A_{\lambda})^{-1}(\begin{smallmatrix}-\alpha\\
\ c\end{smallmatrix})(\begin{smallmatrix}-\alpha\\
\ c\end{smallmatrix})^{T}$ can be found by symbolically solving the $3\times 3$ linear system $(1-A_{\lambda})Z=(\begin{smallmatrix}-\alpha\\
\ c\end{smallmatrix})(\begin{smallmatrix}-\alpha\\
\ c\end{smallmatrix})^{T}$ for $Z\in\mathbb{R}_{\mathrm{sym}}^{2\times 2}$. We perform this computation in Sympy (Meurer et al., ) and obtain the reported result.

## Appendix J Proof of theorem 6

#### Part 1.

Note that the theorem assumption $\lambda_{k}<1$ implies that for some $\varepsilon>0$ all the eigenvalues $\lambda_{k}\leq 1-\varepsilon$. This means that the algorithm parameters never approach the boundary of strict stability of $S_{\lambda}$, given by theorem 3 at $\lambda_{\mathrm{max}}=1$. This observation will be useful throughout the proof.

We start with showing the necessity of the first two conditions $\delta\to 0$ and $q_{0}>0$ that are agnostic to spectrum $\lambda_{k}$. As mentioned in the main text the necessity of $\delta\to 0$ to have $\alpha_{\mathrm{eff}}\to\infty$ follows directly from strict stability bound $\alpha_{\mathrm{eff}}<\frac{4}{\delta}$. Next, we observe that at fixed $\delta$, pairs of strictly stable $(\alpha_{\mathrm{eff}},q_{0})$ are given by $q_{0}>-\delta$ and $0<\alpha_{\mathrm{eff}}<\frac{4-2\delta}{\delta}-\frac{2q_{0}}{\delta}$. These inequalities define a right-angle triangle on $(\alpha_{\mathrm{eff}},q_{0})$ plane with vertices at $\{(\tfrac{4}{\delta},-\delta),\;(0,-\delta),\;(0,2-\delta)\}$, which can serve as a useful geometric picture of relevant $(\alpha_{\mathrm{eff}},q_{0})$ values.

The condition $q_{0}>0$ will follow from noisy stability condition $\sum_{k}\lambda_{k}^{2}R_{\lambda_{k}}=O(1)$. Hence, we bound different terms appearing in $R_{\lambda}$, taking into account $\lambda_{k}\leq 1-\varepsilon$, as

$$ $\begin{split}2-2\delta&<\Big{[}\lambda q_{0}+2-\delta\Big{]}<4-2\delta,\\ \varepsilon(4-2\delta)&<\Big{[}4-2\delta-\lambda(2q_{0}+\delta\alpha_{\mathrm{ eff}})\Big{]}<4.\end{split}$ (230) $$

As $\delta\to 0$, all the bounds above can be replaced by some global $\delta$-independent values. Then, we can simplify $R_{\lambda}$ up to a multiplicative constant as

$$ $R_{\lambda}\asymp\frac{\lambda q_{0}^{2}+\delta\alpha_{\mathrm{eff}}}{\lambda( \delta+\lambda q_{0})},$ (231) $$

reproducing the first part of eq. (27). Observe that for $q_{0}\leq 0$ we have $\delta+\lambda q_{0}\geq\delta$, and, therefore, $R_{\lambda}\geq c\tfrac{\alpha_{\mathrm{eff}}}{\lambda}$ with some $c>0$. Then, for the stability sum we would have $\sum_{k}\lambda_{k}^{2}R_{\lambda_{k}}\geq c\alpha_{\mathrm{eff}}\sum_{k}
\lambda_{k}\to\infty,$ since $\alpha_{\mathrm{eff}}\to\infty$. This settles impossibility of $q_{0}\leq 0$ for noisy stability.

Now, we proceed to the last two stability conditions condition, assuming $q_{0}>0$. Then, both nominator and denominator in (231) have only positive terms, and we can write them $\delta\asymp\max(\delta,\lambda q_{0})$ and $\lambda q_{0}^{2}+\delta\alpha_{\mathrm{eff}}\asymp\max(\lambda q_{0}^{2},
\delta\alpha_{\mathrm{eff}})$. The characteristic spectral scales that separate two choices in these maximums are $\lambda^{\mathrm{cr}}_{1}=\tfrac{\delta}{q_{0}}$ and $\lambda^{\mathrm{cr}}_{2}=\tfrac{\delta\alpha_{\mathrm{eff}}}{q_{0}^{2}}$, as given in the main text. Note that we asymptotically have $\lambda^{\mathrm{cr}}_{1}<\lambda^{\mathrm{cr}}_{2}$ since $\alpha_{\mathrm{eff}}\to\infty$ while $q_{0}$ is bounded. This gives the second part of (27)

$$ $R_{\lambda}\asymp\begin{cases}\frac{\alpha_{\mathrm{eff}}}{\lambda},\quad&0< \lambda<\lambda_{1}^{\mathrm{cr}}\\ \frac{\delta\alpha_{\mathrm{eff}}}{\lambda^{2}q_{0}},\quad&\lambda_{1}^{ \mathrm{cr}}<\lambda<\lambda_{2}^{\mathrm{cr}}\\ \frac{q_{0}}{\lambda},\quad&\lambda_{2}^{\mathrm{cr}}<\lambda<\lambda_{\mathrm {max}}-\varepsilon.\end{cases}$ (232) $$

Substituting the above into the stability sum leads to the condition

$$ $\sum_{k}\lambda_{k}^{2}R_{\lambda_{k}}\asymp\alpha_{\mathrm{eff}}\Big{[}\sum_{ k:\lambda_{k}\leq\frac{\delta}{q_{0}}}\lambda_{k}\Big{]}+\frac{\delta\alpha_{ \mathrm{eff}}}{q_{0}}\Big{[}\sum_{k:\frac{\delta}{q_{0}}<\lambda_{k}<\frac{ \delta\alpha_{\mathrm{eff}}}{q_{0}^{2}}}1\Big{]}+q_{0}\Big{[}\sum_{k:\lambda_{ k}\geq\frac{\delta\alpha_{\mathrm{eff}}}{q_{0}^{2}}}\lambda_{k}\Big{]}=O(1)$ (233) $$

Here, the third term is always bounded since $\sum_{k}\lambda_{k}<\infty$ and $q_{0}$ is also bounded. The boundness of the first and second term corresponds to the third and forth stability conditions of theorem 6, respectively. From the argumentation presented above, it is clear that those conditions are both necessary and sufficient.

#### Part 2.

In fact, we can consider even weaker power-law spectrum assumption $\lambda_{k}\asymp k^{-\nu}$, while in the theorem statement we opted for spectral condition (1) in the sake of main text consistency. Assuming $\frac{\delta}{q_{0}}\to 0$, the sums in (233) can be estimated with a standard power-law arguments in terms of characteristic spectral index $k^{\mathrm{cr}}_{1}\asymp(\tfrac{\delta}{q_{0}})^{-\frac{1}{\nu}}$ defined by $\lambda_{k^{\mathrm{cr}}_{1}}=\frac{\delta}{q_{0}}$

$$ $\Big{[}\sum_{k:\lambda_{k}\leq\frac{\delta}{q_{0}}}\lambda_{k}\Big{]}\asymp \left(\frac{\delta}{q_{0}}\right)^{-\frac{1-\nu}{\nu}},\qquad\Big{[}\sum_{k: \frac{\delta}{q_{0}}<\lambda_{k}<\frac{\delta\alpha_{\mathrm{eff}}}{q_{0}^{2}} }1\Big{]}\asymp\left(\frac{\delta}{q_{0}}\right)^{-\frac{1}{\nu}}.$ (234) $$

Then, both the third and fourth stability conditions become $\alpha_{\mathrm{eff}}\delta^{-\frac{1-\nu}{\nu}}q_{0}^{\frac{1-\nu}{\nu}}=O(1)$. Substituting the ansatz $\alpha_{\mathrm{eff}}=\delta^{-h}$ and $q_{0}=\delta^{g}$, we get $\delta^{-h-(1-g)\frac{1-\nu}{\nu}}=O(1)$, which is equivalent to $h\leq(1-g)(1-\tfrac{1}{\nu})$ given in the statement of the theorem. As for the case $g\geq 1$, it breaks the previous assumption $\frac{\delta}{q_{0}}\to 0$. This brings us back to $q_{0}\leq 0$ case, where the first sum in (233) is estimated as $\Omega(\alpha_{\mathrm{eff}})$, leading to the divergence of $\sum_{k}\lambda_{k}^{2}R_{\lambda_{k}}$. Finally, when $h<h_{\mathrm{max}}$ and $g>0$, we can see that all three sums in (233) converge to zero, which translates to $U_{\Sigma}$ at any constant batch size $|B|$.

## Appendix K Convergence of non-stationary algorithm with power-law schedule

In this section, we derive the convergence rate of signal propagator $V_{t}$, given by (77) (with $\mathbf{a}_{t}=0$ in the left scalar product factor), for the power-law ansatz schedule $(\delta_{t},\alpha_{\mathrm{eff},t},q_{0,t})\sim(t^{-\overline{\delta}},t^{
\overline{\alpha}},const)$.

For a non-stationary algorithm, evolution operators $A_{t,\lambda}$ essentially depend on $t$, making the analysis of the respective products $A_{\lambda,t}A_{\lambda,t-1}\ldots A_{\lambda,0}$ very challenging. To overcome this, we rely on the fact that in the considered schedule $\delta_{t}=t^{-\overline{\delta}},\alpha_{\mathrm{eff},t}=t^{\overline{\alpha}}$, the parameters are changing very slowly at large $t$. Specifically, we assume that a significant change of eigenvectors of $A_{\lambda,t}$ also requires a significant relative change of algorithm parameters (e.g. 2 times increase/decrease). To achieve this change on a time window $(t,t+\Delta t)$ for any power-law relation $f_{t}=t^{-F}$, a large window size $\Delta t\sim t$ is required. Let $\mu_{1},\mu_{2},\mu_{3},\mu_{4}$ be the eigenvalues of $A_{\lambda_{t}}$, ordered as $|\mu_{1}|>|\mu_{2}|\geq|\mu_{3}|\geq|\mu_{4}|$, and $\mathbf{S}=(\mathbf{s_{1}}\ldots\mathbf{s}_{4})$ be the matrix of respective eigenvectors. Then, if we manage to pick window size such that $\Delta t\ll t$ and $|\mu_{2}|^{\Delta t}\ll|\mu_{1}|^{\Delta t}\lesssim 1$, for any starting vector $\mathbf{r}\in\mathbb{R}^{2}\otimes\mathbb{R}^{2}$ we have

$$ $\prod_{i=t}^{t+\Delta t-1}A_{\lambda,i}\mathbf{r}\approx\big{(}A_{\lambda,t} \big{)}^{\Delta t}\mathbf{r}=\sum_{e=1}^{4}\mathbf{s}_{e}(\mu_{e})^{\Delta t} \big{(}\mathbf{S}^{-1}\mathbf{r}\big{)}_{e}\approx\mathbf{s}_{1}(\mu_{1})^{ \Delta t}\big{(}\mathbf{S}^{-1}\mathbf{r}\big{)}_{1}.$ (235) $$

Thus, evolving starting vector $\mathbf{r}$ with $A_{\lambda,t}$ aligns the result with eigenvector $\mathbf{s}_{1}$ corresponding to the slowest eigenvalue $\mu_{1}$ of $A_{\lambda,t}$. The magnitude of the evolution result is given, up to initial projection $\big{(}\mathbf{S}^{-1}\mathbf{r}\big{)}_{1}$, by the product of slowest eigenvalue $\mu_{1}$. Extending this argument to many consecutive time windows $(t_{j},t_{j+1})$, with $\Delta t_{i}=t_{i+j}-t_{j}$ growing with $t_{i}$, we naturally arrive to the following *adiabatic approximation*(^1^11our usage of this approximation has many parallels with adiabatic approximation in quantum mechanics [wiki/Adiabatic theorem](https://en.wikipedia.org/wiki/Adiabatic_theorem))

$$ $\prod_{i=1}^{t}A_{\lambda,i}\mathbf{r}\asymp\left(\prod_{i=1}^{t}\mu_{1,i} \right)\mathbf{s}_{1,t},\quad t\to\infty$ (236) $$

We will use this approximation to compute asymptotic of $V_{t}$. Applying (236) to a single term $V_{\lambda,t}$ in $V_{t}$ sum gives

$$ $\begin{split}V_{\lambda,T}\equiv\langle(\begin{smallmatrix}1\\ \mathbf{a}\end{smallmatrix})(\begin{smallmatrix}1\\ \mathbf{a}\end{smallmatrix})^{T},A_{T,\lambda}A_{T-1,\lambda}\ldots A_{1, \lambda}(\begin{smallmatrix}1&0\\ 0&0\end{smallmatrix})\rangle&\asymp\mu_{1,T}\mu_{1,T-1}\ldots\mu_{1,1}\\ &\asymp\operatorname{exp}\left[\int_{1}^{T}\log\mu_{1}(t)dt\right]\equiv \widetilde{V}_{\lambda,T}.\end{split}$ (237) $$

Here in the last transition, we used that our schedule can be smoothly interpolated to non-integer $t$. Then, replacing the sum by the integral can change the result by at most a constant, as soon as $\varepsilon<|\mu_{1,t}|<1$ for all $t$.

Before proceeding to the calculation of the signal propagator as $V_{T}\asymp\sum_{k}c_{k}^{2}\lambda_{k}\widetilde{V}_{\lambda_{k},T}$, let us write down the required properties of eigenvalues of $A_{\lambda}$. First, we remind that, as in sec. 6, we consider $\tau_{2}=0$, which leads to factorizable evolution operator $A_{\lambda}=S_{\lambda}\otimes S_{\lambda}$ with eigenvalues equal to products of eigenvalues of $S_{\lambda}$. With a slight abuse of notation, in the rest of the section we will denote $\mu_{1},\mu_{2}$ the eigenvalues of $S_{\lambda}$. These eigenvalues are given by zeros of characteristic equation $\chi(x,\lambda)=0$, which in variables $(\delta,\alpha_{\mathrm{eff}},q_{0})$ is written as

$$ $x^{2}-x\big{(}2-\delta-\lambda(\delta\alpha_{\mathrm{eff}}+q_{0})\big{)}+1- \delta-\lambda q_{0}=0.$ (238) $$

The respective roots are

$$ $\mu_{1,2}=1-\frac{\delta+\lambda(\delta\alpha_{\mathrm{eff}}+q_{0})\pm\sqrt{E} }{2},\quad E=(\delta+\lambda(\delta\alpha_{\mathrm{eff}}+q_{0}))^{2}-4\lambda \delta\alpha_{\mathrm{eff}}$ (239) $$

In the case $\delta\ll\alpha_{\mathrm{eff}}\ll q_{0}\sim 1$, the solutions of $E(\lambda)=0$ are given by $\lambda_{1}^{\mathrm{cr}}\approx\frac{\delta}{4\alpha_{\mathrm{eff}}}$ and $\lambda_{2}^{\mathrm{cr}}\approx\frac{4\delta\alpha\mathrm{eff}}{q_{0}^{2}}$. They separate values $\lambda$ with real and complex eigenvalues $\mu_{1,2}$. In particular, for $\lambda<\lambda_{1}^{\mathrm{cr}}$ the eigenvalues are real, can be approximated for $\lambda\ll\lambda_{1}^{\mathrm{cr}}$ as $\mu_{1}\approx 1-\alpha_{\mathrm{eff}}\lambda$ and $\mu_{2}\approx 1-\delta$, and $\mu_{1}(\lambda_{1}^{\mathrm{cr}})=\mu_{2}(\lambda_{1}^{\mathrm{cr}})\approx 1
-\frac{\delta}{2}$ (see also fig. 3 and theorem 3).

Coming back to signal propagator $V_{T}$ and its terms $\widetilde{V}_{\lambda,T}$, let us fix some $\lambda\ll 1$ and increase time $T$. At first, the chosen eigenvalue $\lambda$ will be very small compared to $\lambda_{1}^{\mathrm{cr}}(t)$. But as the schedule advances, it will enter the complex region at a transition time $t_{*}(\lambda):\lambda_{1}^{\mathrm{cr}}(t)=\lambda$. By substituting $\delta_{t},\alpha_{\mathrm{eff},t}$ from our schedule, we get $t_{*}(\lambda)\approx(4\lambda)^{-\frac{1}{\overline{\alpha}+\overline{\delta}}}$. Although for $T>t_{*}(\lambda)$ adiabatic approximation is no longer applicable, we can use it for $T\ll t_{*}(\lambda)$ to obtain (recall that $\mu_{1}(A_{\lambda})=\big{(}\mu_{1}(S_{\lambda})\big{)}^{2}$)

$$ $\widetilde{V}_{\lambda,T}=\operatorname{exp}\left[\int_{1}^{T}2\log\mu_{1}(t) dt\right]\approx\operatorname{exp}\left[-\int_{1}^{T}2\lambda t^{\overline{ \alpha}}dt\right]\approx\operatorname{exp}\left[-\frac{2}{1+\overline{\alpha}} \lambda T^{1+\overline{\alpha}}\right]$ (240) $$

Let us look at the convergence reached by the end of adiabatic approximation applicability time $t_{*}(\lambda)\approx(4\lambda)^{-\frac{1}{\overline{\alpha}+\overline{\delta}}}$. Substituting this value into $\widetilde{V}_{\lambda,T}$ leads to

$$ $\widetilde{V}_{\lambda,t_{*}(\lambda)}\approx\operatorname{exp}\left[-C_{ \widetilde{V}}\lambda^{\frac{1-\overline{\delta}}{\overline{\alpha}+\overline{ \delta}}}\right],\quad C_{\widetilde{V}}=\frac{2^{\frac{2+\overline{\alpha}- \overline{\delta}}{\overline{\alpha}+\overline{\delta}}}}{1+\overline{\alpha}}.$ (241) $$

We see that, for $\overline{\delta}<1$(^2^22Note that Jacobi scheduled HB asymptotically uses $\overline{\delta}=\overline{\alpha}=1$, which suggests that the requirement $\overline{\delta}<1$ is reasonable and approaches the optimal scaling of momentum term $\delta$. Yet, the current approach based on adiabatic approximation cannot treat $\overline{\delta}=1$ case.), by the end of adiabatic approximation applicability period, the signal at the chosen eigenvalue $\lambda$ has already converged exponentially as $\lambda\to 0$. Moreover, since at all $t$ we have stable $S_{\lambda}$ with $|\mu_{1,2}|<1$, $\widetilde{V}_{\lambda,T}$ can not increase back to non-negligible values at some $T>t_{*}(\lambda)$. Summarizing, in the limit $T\to\infty$ we can use the expression (240) at all $T$, and also at non-vanishing $\lambda$ since larger eigenvalues are converging even faster.

Now, we can combine all eigenvalues under power-law conditions (1), and obtain $T\to\infty$ asymptotic of the signal propagator of AM1 algorithm

$$ $\begin{split}V_{T}&\asymp\sum_{k}c_{k}^{2}\lambda_{k}\widetilde{V}_{\lambda_{k },T}\approx\int_{0}^{\lambda_{\mathrm{max}}}Q\zeta\lambda^{\zeta-1}\widetilde{ V}_{\lambda,T}d\lambda\\ &\approx\int_{0}^{\lambda_{\mathrm{max}}}Q\zeta\lambda^{\zeta-1}e^{-\frac{2}{1 +\overline{\alpha}}\lambda T^{1+\overline{\alpha}}}d\lambda\approx Q\Gamma( \zeta+1)\left(\frac{2}{1+\overline{\alpha}}\right)^{-\zeta}T^{-\zeta(1+ \overline{\alpha})}.\end{split}$ (242) $$

The integral here is computed by reducing it to Gamma function integral with linear change of variables $z=\frac{2}{1+\overline{\alpha}}\lambda T^{1+\overline{\alpha}}$.

## Appendix L Experiments

In all the experiments we consider Memory-1 algorithms with $a=0$. There are many equivalent ways to implement a Memory-1 algorithm. In practice, we found it the most convenient to use a minimal modification of standard *SGD with momentum* algorithm (referred to as HB in this work) implemented in modern deep learning frameworks(^3^33See, for example, [pytorch/SGD](https://pytorch.org/docs/stable/generated/torch.optim.SGD.html))

$$ $\mathbf{u}_{t+1}=\beta\mathbf{u}_{t}-\nabla L(\mathbf{w}_{t}),\qquad\begin{ cases}\text{HB}:\quad&\mathbf{w}_{t+1}=\mathbf{w}_{t}+\alpha\mathbf{u}_{t+1}\\ \text{Memory-1}:\quad&\mathbf{w}_{t+1}=\mathbf{w}_{t}+\alpha_{2}\mathbf{u}_{t+ 1}-\alpha_{1}\nabla L(\mathbf{w}_{t}),\\ \end{cases}$ (243) $$

This way, one can use an update step of standard SGD optimizer class with learning rate $\alpha=\alpha_{2}$, and then manually subtract parameter gradients with new learning rate $\alpha_{1}$

The new learning rates in parametrization (243) are expressed in terms of our main notations as $\beta=1-\delta$, $\alpha_{1}=\tfrac{q_{0}}{\beta},\alpha_{2}=\delta(\alpha_{\mathrm{eff}}-\alpha
_{1})$. Then, the accelerated region $\delta^{-1}\gg\alpha_{\mathrm{eff}}\gg q_{0}=const$ translates to $\delta\ll\alpha_{2}\ll\alpha_{1}=const$.

This gives another perspective on the accelerated AM1 algorithm we identified in sec. 6. To achieve acceleration, we still want, as in HB, to make the velocity vector more inertive (i.e. making the ball “heavier”). However, while the optimally accelerated HB adds this velocity with $\alpha=O(1)$ learning rate, we add it with vanishing $\alpha_{2}=O(\delta\alpha_{\mathrm{eff}})\to 0$ learning rate but compensate with an extra “kick” term $\nabla L(\mathbf{w}_{t})$ with constant $\alpha_{1}=O(1)$ learning rate. Apparently, this combination leads both acceleration and bounded effect of mini-batch noise.

Next, we discuss our two experimental setups.

### L.1 Gaussian data

For Gaussian data experiments, we generate inputs $\mathbf{x}\sim\mathcal{N}(0,\Lambda)$ from Gaussian distribution with diagonal covariance $\Lambda=\operatorname{diag}(\lambda_{1},\lambda_{2},\ldots\lambda_{M})$, leading to diagonal Hessian $\mathbf{H}=\Lambda$. For the optimal parameters, we simply take the vector of target coefficients $\mathbf{w}_{*}=(c_{1},c_{2},\ldots,c_{M})$. Then, we set ideal power-laws for $\lambda_{k}=k^{-\nu}$ and $c_{k}^{2}=k^{-\kappa-1}$ which satisfy our asymptotic power-law conditions (1) with $\zeta=\frac{\kappa}{\nu}$. This way, we create a simple setting where all the loss curves are expected to have a pronounced power-law shape, and we will be able to directly compare experimental convergence exponents with theoretically predicted ones. We set $M=10^{4}$ to ensure that truncation of the spectrum has negligible finite size effects and the experiments are well described by the asymptotic regime $t,k\to\infty$.

For the experiments in fig. 1 and fig. 4, we pick spectral exponents $\zeta=0.5,\nu=3$. For the AM1 schedule exponent we use $\overline{\delta}=0.95,\overline{\alpha}=\overline{\delta}(1-\tfrac{1}{\nu})$, at the edge of allowed AM1 exponents values. In both figures 1 and 4, we compare empirical loss trajectories with power-law curves (dashed lines) with predicted convergence rate, which show a good agreement between each other. In fig. 4 we try more $\overline{\delta},\overline{\alpha}$ values, and find good agreement of the experimental loss rates with predicted $L_{t}=O(t^{-\zeta(1+\overline{\alpha})})$.

### L.2 MNIST

We design the MNIST experiment to be the opposite of the Gaussian data setting to test our results in more general scenario beyond quadratic problems with asymptotically power-law spectrum. Specifically, we apply a single-hidden layer neural network to a flattened MNIST images $\mathbf{x}\in\mathbb{R}^{768}$:

$$ $f(\mathbf{W}_{1},\mathbf{W}_{2},\mathbf{b},\mathbf{x})=\frac{1}{\sqrt{n}} \mathbf{W}_{2}\operatorname{ReLU}(\mathbf{W}_{1}\mathbf{x}-\mathbf{b})\in \mathbb{R}^{10},$ (244) $$

where $\mathbf{W}_{1}\in\mathbb{R}^{n\times 768},\;\mathbf{b},\mathbf{x}\in\mathbb{R}
^{768},\;\mathbf{W}_{2}\in\mathbb{R}^{10\times n}$, with model width taken as $n=1000$. Then, we consider one hot encoded vectors $y=(0,\ldots,1,\ldots,0)\in\mathbb{R}^{1}0$ as quadratic loss $L_{B}=\frac{1}{2b}\sum_{i=1}^{b}\|f(\mathbf{x}_{i})-\mathbf{y}_{i}\|^{2}$ as an empirical loss function to be minimized.

The above setting, although much simpler than modern deep-learning problems, provides an example of non-convex optimization function. In particular, if one linearizes $f(\mathbf{W}_{1},\mathbf{W}_{2},\mathbf{b},\mathbf{x})$ with respect to parameters on the current optimization step, the resulting quadratic approximation to the original loss function will evolve during and the respective spectral distributions $\lambda_{k},c_{k}$ are expected to deform significantly from the roughly power-law shapes with $\nu\approx 1.33$ and $\kappa=\zeta\nu\approx 0.33$ reported in Velikanov & Yarotsky () for the limit $n\to\infty$.

Yet, on figure1 we can still see that both Jacobi-scheduled HB and AM1 display qualitatively similar behavior to that of quadratic problems with power-law spectrum, though with slightly different exponents than can be expected from $\nu\approx 1.33,\;\kappa\approx 0.33$ estimation. Most importantly, Jacobi-scheduled HB retains its instability in stochastic settings, while AM1 is stable at least for extremely long $t=10^{5}$-steps duration of the experiment.

AM1 algorithm in fig.1 uses $\overline{\delta}=1$ and $\overline{\alpha}=0.5$. This chosen effective learning rate exponent is, in fact, larger than the maximal stable exponent $\overline{\alpha}_{\mathrm{max}}=\overline{\delta}(1-\tfrac{1}{\nu})\approx 0.25$ that could be estimated with MNIST eigenvalue exponent $\nu\approx 1.33$. Yet, we don’t see any signs of divergence on the whole training duration, but which were observed on Gaussian data experiments in fig. 4. That might suggest that non-quadratic models could exhibit different stability conditions due to some recovery mechanism, with the catapult effect Lewkowycz et al. () being one example of such recovery. We provide MNIST loss trajectories for extra values of Batch size and algorithm parameters in fig. 5.

Figure: Figure 4: Comparison of plain GD and AM1 on a grid of $\overline{\delta},\overline{\alpha}$. We run both algorithms in full-batch setting and mini-batch with $|B|=1$, leading to 4 trajectories on each subplot. Also, we choose $\alpha=0.1$ for both algorithms, and $\alpha_{\mathrm{eff},t}=0.1t^{\overline{\alpha}}$ to avoid any non-asymptotic divergence. For each couple $\overline{\delta},\overline{\alpha}$, we color in green the subplot titles if stability condition $\overline{\alpha}\leq\overline{\delta}(1-\tfrac{1}{\nu})$ holds, and in red otherwise. We observe that all empirically diverged trajectories are predicted correctly by this stability condition, while a single trajectory $(\overline{\delta},\overline{\alpha})=(1.0,0.75)$ is expected to diverge but demonstrate only large but bounded noise fluctuations. Finally, along each experimental trajectory, we plot an exact power-law (dashed) line with rate: $t^{-\zeta}$ for GD and $t^{-\zeta(1+\overline{\alpha})}$ for AM1. Again, we observe very good agreement between theoretical prediction and empirical rates, validating adiabatic approximation of sec. K. Moreover, the rate $t^{-\zeta(1+\overline{\alpha})}$ remains valid even for the case $\overline{\delta}=1$ not covered by adiabatic approximation.
Refer to caption: x2.png

Figure: Figure 5: Trajectories of plain GD, Jacobi scheduled HB and AM1 for different learning rates and batch sizes. For each subplot, learning rate $\alpha$ of GD and Jacobi-scheduled HB is equal to “kick” learning rate $\alpha_{1}$ of AM1. Also, Jacobi-scheduled HB and AM1 have exactly the same values of momentum $\beta_{t}=1-\tfrac{2}{t}$ on each iteration. We observe that AM1 is indeed much more stable than Jacobi-scheduled HB for a range of learning rates and batch sizes. Moreover, it is stable for larger values of $\overline{\alpha}$ than estimated maximal value $\overline{\alpha}_{\mathrm{max}}\approx 0.25$. Another interesting observation is that the divergent trajectories of Jacobi scheduled HB, unlike for pure quadratic problems, only get stuck around the loss of random prediction instead of diverging to infinity.
Refer to caption: x3.png

## Peer Reviews

confidence: 2
contribution: 3
initial_rating: 8
presentation: 3
questions: As mentioned, I wonder if there are any ways of improving the readability of the technical portions with further addition of intuitions and heuristics. This wouldn’t change the technical contributions of the paper but would make the paper more accessible for the readers and thus facilitate the propagation of ideas in it.
soundness: 4
strengths: Overall, this is a very good paper. I think it is very methodical in its development of the memory-M framework - starting with existing issues of accelerated HB for SGD, moving to the description of the generalized framework for GD and SGD, loss decomposition and improved algorithm and convergence properties. The topic coverage is very comprehensive, covering the expected questions raising with the development of the generalized framework. There are clearly a lot of efforts made towards the balance between rigor and readability - with moving a lot of the technical details into the appendix (e.g. the phrasing of Theorem 2). The contributions are novel to the field of stochastic optimization
summary: This paper establishes a theoretical framework for the so-called memory-M algorithms - first order algorithms that are augmented with storage for M vectors in the gradient space, and uses a linear combination of them to compute the gradients and modify weight updates. The usefulness of this framework is that it is a natural generalization of Heavy Ball and a number of other first order algorithms. The work provides an equivalence of this formulation with recursive multi-step first-order algorithms. The work also provides a generalized loss decomposition into signal and noise propagations. The paper also provides a detailed description of convergence properties of memory-M (and focusing on memory-1) algorithms - including the main parameters controlling this convergence and an algorithm providing equivalent convergence rates for the SGD case.
weaknesses: An expected consequence of the comprehensiveness of the paper is its density - making it relatively hard to read. As mentioned before, there was clearly a lot of effort put towards reorganization and moving details into the appendix - but I wonder if more could be done in that direction. For example, the Loss expansion section with the introduction of propagators seemed too dense - and was lacking intuitions that would facilitate further reading of the paper, since those concepts were heavily relied on for further developments in the paper. Those intuitions were present in the earlier sections and were thus helpful for speeding up understanding of the paper

---

confidence: 3
contribution: 4
initial_rating: 8
presentation: 2
questions: 1. Could such a framework be extended to non-IID batching, that is, batching with replacement (recent work [1,2])?

2. Could you say anything about the generalization error using [3]?

3. Could we add label noise or regularization?

[1] Lok, Jackie, Rishi Sonthalia, and Elizaveta Rebrova. "Discrete error dynamics of mini-batch gradient descent for least squares regression." arXiv preprint arXiv:2406.03696 (2024).

[2] Beneventano, Pierfrancesco. "On the trajectories of sgd without replacement." arXiv preprint arXiv:2312.16143 (2023).

[3]  Wang, Y., Sonthalia, R. &amp; Hu, W.. (2024). Near-Interpolators: Rapid Norm Growth and the Trade-Off between Interpolation and Generalization. <i>Proceedings of The 27th International Conference on Artificial Intelligence and Statistics</i>, in <i>Proceedings of Machine Learning Research</i> 238:4483-4491 Available from https://proceedings.mlr.press/v238/wang24l.html.
soundness: 4
strengths: The paper's biggest strength is the establishment of strong convergence rates for iterative methods. Mini-batched algorithms are difficult to analyze; hence, establishing the convergence rate for regression in the noiseless regime is interesting. 

Additionally, the method allows for controlling the constant and shows that it can also be arbitrarily improved. 

I think both of these are important contributions. 

In addition, the paper introduces a general framework for analyzing such methods. i think this is an important contribution as well.
summary: This paper introduces a framework for analyzing iterative optimization methods, referred to as memory-$M$ algorithms, which allow the gradient to be evaluated at shifted points relative to $w_t$ additionally a shift term is added to the gradient when doing the update. The authors explore both stationary (with fixed update parameters) and non-stationary versions of these algorithms. They establish that memory-$M$ algorithms can achieve the optimal convergence rate for data with a power-law spectrum, even under mini-batch noise. Notably, they propose a memory-1 algorithm with a time-dependent power-law schedule, which provides a heuristic pathway to accelerated convergence, supported by empirical results.


**Justification for the score**

While I think the result is important in establishing convergence results for a large family of methods, I think the presentation of the paper could be significantly improved.
weaknesses: The writing of the paper is a major weakness. It is very dense notationally, not necessarily conceptually. This makes it quite challenging to build intuition as one reads the paper. One possible suggestion is that the paper develops the results initially in the $a = 0$ as the shifts are not needed for the main contributions and then provide extensions, possibly in the appendix. The other is to hide some of the intermediate steps to provide space for more intuition for the different objects. 

One thing that threw me off is using capital letters for the propagators as these are scalars.

---

confidence: 2
contribution: 3
initial_rating: 3
presentation: 1
questions: Can authors convince me that their work fits the spirit of this conference? This question is not obvious to me at this point.....
soundness: 1
strengths: The work is a definite interest to read, and the results seem to hold a lot of promise.
summary: This paper explores how to make stochastic gradient descent (SGD) faster on certain types of optimization problems, specifically those that involve a noisy mini-batch setting. The authors examine how adding "memory" (extra helper variables) to SGD can help. They introduce "memory-$M$" algorithms, where $M$ is the number of extra velocity vectors added. Through mathematical analysis, they show that these memory-$M$ algorithms can achieve better performance by improving certain constants that affect convergence speed, even though they don’t change the basic rate of convergence.

The authors also propose a new "memory-1" algorithm, which has only one extra velocity vector and uses a special time-based schedule for learning. This algorithm is shown to improve performance both theoretically and in experiments.
weaknesses: - Introduction: The reader who is not an expert in this field will find this introduction difficult to understand. I recommend that it be significantly reworked using various welcome tricks, such as tables. In general, the text is difficult to read....

- I don't fully understand about the experimental part. Is it just not there?

---

confidence: 2
contribution: 3
initial_rating: 6
presentation: 2
questions: 1. I do not fully understand the image on the right of Figure 2. What are $t_1, ..., t_4$ supposed to denote or perhaps more generally what does it mean that a signal/noise propagator is "long" or "short"?
2. I have another question regarding the stability of the maps $A_\lambda$: If I understood it correctly, $\lambda$ corresponds to the eigenvalues of the matrix $\mathbf{H}$ (line 201). Thus, it seems to me that we have to assume $\lambda > 0$ in general. Does Theorem 3 then imply that one can only have strict stability of $A_{\lambda}$ if $\lambda_{\max}$ is sufficiently small? 
3. What are the values of $\xi_V$ and $\xi_U$ respectively in Eq. (21) and (22)? 
4. Have the authors also considered non-quadratic settings and how does it perform here?
soundness: 4
strengths: - The authors are able to (heuristically) derive an accelerated and stable rate of convergence in their algorithm AM1 and verify the rate empirically. (However I am not too familiar with the field to judge the novelty of this contribution.) 
- The derivation of the stability of stationary stable memory-M algorithms seems rigorous and it is interesting to see that stationary stable algorithms cannot achieve a faster rate than plain GD even using arbitrary numbers $M$ auxiliary vectors. 
- The analysis of the convergence constant $C_L$ and the possibility of making it arbitrarily small for stationary stable memory-1 algorithms and how it relates to the effective learning rate $\alpha_{eff}$ was also insightful.
summary: In this work, the authors address the open problem of theoretical feasible acceleration of mini-batch SGD-type algorithms (i.e. in the stochastic setting) on infinite-dimensional quadratic problems, particularly with power-law spectrum. 
The authors introduce memory-M algorithms with M auxiliary velocity vectors and show the equivalence to the multistep recurrence formulation if the loss $L$ is quadratic and hyperparameters of the algorithm are stationary. 
By expanding the loss in terms of so-called signal and noise propagators, the authors show that the loss convergence $L_t \sim C_l t^{-\xi}$ of stationary stable memory-M algorithms retain the exponent $\xi$ of plain GD, but can have different constants $C_L$. Finally the authors prove that memory-1 algorithms can make $C_L$ arbitrarily small while maintaining stability. The theoretical results are empirically verified on synthetic Gaussian data and MNIST using a shallow ReLU network.
weaknesses: - The paper seemed quite convoluted with many different symbols, equations and indices due to its mathematic rigor and it could benefit from restructuring it such that it becomes more accessible. 
- The empirical validation of the results remains quite limited with only experiments on synthetic Gaussian data and MNIST. Particularly, the spectrally expressible approximation holds exactly for gaussian data and approximately for MNIST classification. It would be interesting to see also experiments where this approximation is less exact.
- It is seems like the optimal rate depends on the right choice of hyperparameters, such as $\alpha_0$, $q_0$ and $\delta$. However it is unclear to me how to select or estimate them in practice. This limits its application to real-world problems.

Meta-review: This paper establishes a theoretical framework for the so-called memory-M algorithms - first order algorithms that are augmented with storage for M vectors in the gradient space, and uses a linear combination of them to compute the gradients and modify weight updates. The usefulness of this framework is that it is a natural generalization of Heavy Ball and a number of other first order algorithms. The work provides an equivalence of this formulation with recursive multi-step first-order algorithms. The work also provides a generalized loss decomposition into signal and noise propagations. The paper also provides a detailed description of convergence properties of memory-M (and focusing on memory-1) algorithms - including the main parameters controlling this convergence and an algorithm providing equivalent convergence rates for the SGD case.

- Overall, the work is methodical in its development of the memory-M framework - starting with existing issues of accelerated HB for SGD, moving to the description of the generalized framework for GD and SGD, loss decomposition and improved algorithm and convergence properties. 
- The topic coverage is comprehensive. 
- The contributions are novel to the field of stochastic optimization

However, the reviewers agreed that the paper is not very well written, mainly due to it being overly dense and technical. The empirical validation of the results is limited with only experiments on synthetic Gaussian data and MNIST; however, since this is a theoretical work, I am putting lower weight on this.

Three reviewers rated the paper with "8"; and one rated it as "3". However, this negative review did not offer substantial reasons for the criticism raised, and I did not see how the score 3 was justified. Hence, I put a lower weight on this review.

Overall, I am happy to recommend acceptance. 

---

I wish to add that the authors should consider (I am not forcing them to) citing & discussing these works in relation to their results, which seem closely relevant, since they consider acceleration of SGD for quadratic problems:

1. Nicolas Loizou and Peter Richtarik. Momentum and stochastic momentum for stochastic gradient, Newton, proximal point and subspace descent methods, arXiv:1712.09677, 2017
2. Robert M. Gower, Filip Hanzely, Peter Richtarik and Sebastian Stich. Accelerated stochastic matrix inversion: general theory and speeding up BFGS rules for faster second-order optimization, NeurIPS 2018
3. Peter Richtarik and Martin Takac. Stochastic reformulations of linear systems: algorithms and convergence theory, arXiv:1706.01108, 2017

Decision: Accept (Poster)