g_{n}^{T}}{\langle\Delta g_{n},\Delta
w_{n}\rangle}\right)\hat{H}_{n-1}^{-1}\left(v-\tfrac{\langle\Delta
g_{n},v\rangle}{\langle\Delta g_{n},\Delta w_{n}\rangle}\Delta
w_{n}\right)+\tfrac{\langle\Delta w_{n},v\rangle}{\langle\Delta g_{n},\Delta
w_{n}\rangle}\Delta w_{n}.$

Assuming $\langle\Delta g_{n},\Delta w_{n}\rangle$ is precalculated, we only
need to calculate a scalar product and subtract two vectors until we have
reached a point where we need to only apply $\hat{H}_{n-1}^{-1}$ to a vector.
Going through this recursively and assuming $\hat{H}_{0}=\mathbb{I}$, we only
need $O(nd)$ operations to calculate $\hat{H}_{n}^{-1}v$ this way. For $n\leq
d$ this can be cheaper than matrix vector multiplication requiring $O(d^{2})$
operations. Unfortunately after some time it becomes more and more expensive.

So limited memory BFGS simply throws away information older than $k$
iterations, effectively restarting the estimation of the inverse Hessian with
$\hat{H}_{n-k}^{-1}=\mathbb{I}$, guaranteeing a cost of $O(kd)=O(d)$.

###  4.5 Global Newton Methods

So far we have only motivated the Newton-Raphson method on quadratic
functions. Assuming we are in some local area around the minimum, it can be
shown that successive Newton-Raphson steps result in convergence. But these
proofs rely on the fact that the Hessian is locally accurate enough. Therefore
even these proofs implicitly bound the rate of change of the curvature
(Hessian). For global Newton methods we just go about this more explicitly,
like we did with the rate of change of gradient descent (Section 1.3).

####  4.5.1 Cubic Regularization of Newton’s Method

In Section 1.3 we used Lipschitz continuity of the gradient to get an upper
bound on our Loss (1.17)

$\displaystyle\mathcal{L}(\theta)\leq\mathcal{L}(w)+\langle\nabla\mathcal{L}(w),\theta-w\rangle+\tfrac{L}{2}\|\theta-w\|^{2},$

which we then optimized over in Lemma 1.3.1 to find that gradient descent with
learning rate $\tfrac{1}{L}$ is optimal with regards to this bound. The
Lipschitz constant $L$ bounds the second derivative and thus the rate of
change of the gradient. This allows us to formulate for how long we can trust
our gradient evaluation in some sense.

In the case of the Newton-Raphson method, we need a similar rate of change
bound on the Hessian. In this case we want to bound the rate of change of the
second derivative so we assume a Lipschitz continuous second derivative (or in
other words bounded third derivative). This similarly results in the upper
bound

$\displaystyle\mathcal{L}(\theta)\leq\mathcal{L}(w)+\langle\nabla\mathcal{L}(w),\theta-w\rangle+\tfrac{1}{2}\langle\nabla^{2}\mathcal{L}(w)(\theta-w),\theta-w\rangle+\tfrac{M}{6}\|\theta-w\|^{3}.$

Minimizing this upper bound results in the cubic regularization of Newton’s
method [50, Section 4.1]. The solution of which can unfortunately not be
written explicitly.

Since this method is history agnostic again, Nesterov reuses estimating
sequences, to derive an “accelerated cubic Newton scheme” [50, Section 4.2]
similar to Nesterov momentum, i.e. the accelerated gradient descent.

####  4.5.2 Levenberg-Marquard Regularization

Instead of using a cubic bound motivated by a bounded third derivative, one
could still use a quadratic bound

$\displaystyle\mathcal{L}(\theta)$
$\displaystyle\leq\mathcal{L}(w)+\langle\nabla\mathcal{L}(w),\theta-w\rangle+\tfrac{1}{2}\langle\nabla^{2}\mathcal{L}(w)(\theta-w),\theta-w\rangle+\tfrac{M}{2}\|\theta-w\|^{2}$
$\displaystyle=\mathcal{L}(w)+\langle\nabla\mathcal{L}(w),\theta-w\rangle+\tfrac{1}{2}\langle[\nabla^{2}\mathcal{L}(w)+M\mathbb{I}](\theta-w),\theta-w\rangle$

which might perhaps be motivated by a bounded second derivative

$\displaystyle\mu\mathbb{I}\precsim\nabla^{2}\mathcal{L}(w)\precsim
L\mathbb{I}$

which would justify the bound $M:=L-\mu$, as we have

$\displaystyle\mathcal{L}(\theta)$
$\displaystyle\leq\mathcal{L}(w)+\langle\nabla\mathcal{L}(w),\theta-w\rangle+\underbrace{\tfrac{L}{2}\|\theta-w\|^{2}\mathrlap{.}}_{\leq\tfrac{1}{2}\langle[\nabla^{2}\mathcal{L}(w)+M\mathbb{I}]\mathrlap{(\theta-w),\theta-w\rangle}}$
(4.10)

Of course in practice one would probably use smaller $M$, hoping that the
changes are not so drastic that the smallest eigenvalues become the upper
bound. While minimizing the cubic bound was difficult, minimizing this bound
results in

$\displaystyle
w_{n+1}=w_{n}-[\nabla^{2}\mathcal{L}(w_{n})+M\mathbb{I}]^{-1}\nabla\mathcal{L}(w_{n}).$
(4.11)

For this result one simply replaces $\nabla^{2}\mathcal{L}(w)$ with
$[\nabla^{2}\mathcal{L}(w)+M\mathbb{I}]$ in the derivation of the Newton-
Raphson step.

For very small $M$ this method behaves like the Newton-Raphson method, while
for large $M$ the step in every eigenspace is roughly $\tfrac{1}{M}$ and the
method therefore behaves like gradient descent. Levenberg-Marquard
regularization (4.11) can therefore also be viewed as a way to interpolate
between the faster Newton-Raphson method and the more stable gradient descent.

It moves small eigenvalues further away from zero, making the Hessian easier
to invert (without numerical errors). And it can turn negative eigenvalues
into (small) positive eigenvalues, which results in large movements in these
eigenspaces in the correct direction (away from the vertex). On the other hand
it might also fail to make large negative eigenvalues positive, turning them
into small negative eigenvalues causing big jumps towards the vertex in this
eigenspace. So if we pick $M$ too small, this method can make things worse. If
it is too big, it will just slow down convergence and behave roughly like
gradient descent with a small learning rate.

##### Theoretical Guarantees

As we have used the Lipschitz continuity to justify its bound (4.10), we are
already starting out with a worse bound than the one we have used for the
convergence of gradient descent. So this approach does not result in better
convergence guarantees. It is therefore more of a heuristic than theoretically
motivated. In his motivation for the cubic regularization Nesterov states
“Unfortunately, none of these approaches seems to be useful in addressing the
global behavior of second-order schemes.” [50, p. 242]. This also includes
trust bounds (Subsection 4.5.3) and damped Newton (Subsection 4.5.4).

##### Hessian Free Optimization

Since (4.11) still requires the inversion of a matrix, which is infeasible in
high dimension, [46] suggests a Hessian free variant. He notes that the bound

$\displaystyle
B_{\mathcal{L}(w)}(\theta):=\mathcal{L}(w)+\langle\nabla\mathcal{L}(w),\theta-w\rangle+\tfrac{1}{2}\langle[\nabla^{2}\mathcal{L}(w)+M\mathbb{I}](\theta-w),\theta-w\rangle$
(4.12)

we are optimizing over is a quadratic, (and hopefully convex) function and
suggests using conjugate gradient descent (Section 4.3) on this bound.
Conjugate gradient descent does not require inversion of the hessian and it
does not even require the hessian itself. It only requires the ability to
evaluate $\nabla^{2}\mathcal{L}(w)v$ for any vector $v$. For this evaluation
[46] suggests using either the finite differences approximation

$\displaystyle\nabla^{2}\mathcal{L}(w)v\approx\frac{\nabla\mathcal{L}(w+\epsilon
v)-\nabla\mathcal{L}(w)}{\epsilon}.$

or adapting automatic differentiation frameworks to this task, as discussed in
[54]. This avoids storage cost $O(d^{2})$ of the Hessian. Stopping conjugate
gradient descent early, avoids most of the computational costs as well. This
is justified with the fact that the first few iterations of conjugate gradient
achieve most of the optimization.

##### Krylov Subspace Descent

Let us assume we are using Hesse Free Optimization and run conjugate gradient
descent for $m$ steps on $w_{n}$, then we know that conjugate gradient will
find the optimal $w_{n+1}$ in the Krylov subspace $\mathcal{K}_{m-1}$

$\displaystyle w_{n+1}=\arg\min\\{B_{\mathcal{L}(w_{n})}(w):w\in
w_{n}+\mathcal{K}_{m-1}\\},$ (Hessian free)

where $B_{\mathcal{L}(w_{n})}$ is the upper bound defined in (4.12) and the
Krylov subspace is accordingly

$\displaystyle\mathcal{K}_{m}$
$\displaystyle=\text{span}\\{[\nabla^{2}\mathcal{L}(w_{n})+M\mathbb{I}]^{0}\nabla\mathcal{L}(w_{n}),\dots,[\nabla^{2}\mathcal{L}(w_{n})+M\mathbb{I}]^{m}\nabla\mathcal{L}(w_{n})\\}$
$\displaystyle=\text{span}\\{[\nabla^{2}\mathcal{L}(w_{n})]^{0}\nabla\mathcal{L}(w_{n}),\dots,[\nabla^{2}\mathcal{L}(w_{n})]^{m}\nabla\mathcal{L}(w_{n})\\},$

where we get the second equation as $M\mathbb{I}\nabla\mathcal{L}(w_{n})$ is
just a multiple of the first element and therefore the linear hulls are the
same.

Noticing this fact, [77] argue that we can avoid choosing the “correct”
regularization parameter $M$, if we select

$\displaystyle w_{n+1}=\arg\min\\{\mathcal{L}(w):w\in
w_{n}+\mathcal{K}_{m-1}\\}$ (Krylov subspace descent)

instead of the (Hessian free) recursion, which requires
$B_{\mathcal{L}(w_{n})}$ and in turn the regularization parameter $M$ for this
upper bound. Of course at this point we can not use conjugate gradient descent
anymore. They then suggest using L-BFGS to solve this optimization problem.
Why not use L-BFGS in the first place? Well, the problem

$\displaystyle\min\\{\mathcal{L}(w_{n}+w):w\in\mathcal{K}_{m-1}\\}$

can be reparametrized using an (orthogonal) basis $d^{(1)},\dots,d^{(m)}$ of
$\mathcal{K}_{m-1}$ into

$\displaystyle\min_{\alpha}\mathcal{L}\left(w_{n}+\sum_{k=1}^{m}\alpha_{k}d^{(k)}\right)$

which is easier to solve as the dimensionality of the problem is greatly
reduced for $m\ll d$. So in some sense it is similar to line search, but for
$m$ dimensions instead of one. Krylov subspace descent essentially bets on the
conjecture that the Krylov subspaces represent the “important” directions for
optimization. Just like the gradient direction is the optimal direction for
line search.

####  4.5.3 Trust Bounds

Another way to make sure that we can trust the Hessian is to have literal
trust bounds, i.e.

$\displaystyle\min_{\theta\in
B(w)}\underbrace{\mathcal{L}(w)+\langle\nabla\mathcal{L}(w),\theta-w\rangle+\langle\nabla^{2}\mathcal{L}(w)(\theta-w),\theta-w\rangle}_{=T_{2}\mathcal{L}(\theta)}$

where $B(w)$ is some trust set around $w$. For the euclidean ball

$\displaystyle B(w)=\\{\theta:\|\theta-w\|\leq\epsilon\\}$

this restriction can be converted into the secondary condition

$\displaystyle g(\theta):=\frac{\|\theta-w\|^{2}-\epsilon^{2}}{2}\leq 0.$

The Lagrangian function to this constrained optimization problem is

$\displaystyle L(\theta,M)$
$\displaystyle=T_{2}\mathcal{L}(\theta)+Mg(\theta)$
$\displaystyle=\mathcal{L}(w)-\tfrac{\epsilon^{2}}{2}+\langle\nabla\mathcal{L}(w),\theta-w\rangle+\langle[\nabla^{2}\mathcal{L}(w)+M\mathbb{I}](\theta-w),\theta-w\rangle.$

As the constant $\frac{\epsilon}{2}$ does not change this optimization
problem, its solution is the same as in Levenberg-Marquard regularization. So
trust bounds are the more general concept.

In particular [16] define the bound

$\displaystyle B(w)=\\{\theta:\|\theta-w\|_{|H|}\leq\epsilon\\}$

where $|H|$ is the hessian with its eigenvalues replaced by their absolute
values and $\|x\|_{|H|}:=x^{T}|H|x$. They show that this trust bound results
in

$\displaystyle
w_{n+1}=w_{n}-|\nabla^{2}\mathcal{L}(w_{n})|^{-1}\nabla\mathcal{L}(w_{n}).$
(saddle-free Newton)

This method avoids moving towards saddle points and maxima as it does not flip
the direction with negative eigenvalues in these eigenspaces. Since the
calculation of the eigenvalues of the Hessian is generally intractable for
high dimensional problems, [16] use the dimension reduction idea from Krylov
subspace descent to calculate and take the absolute value of the eigenvalues
of the reduced Hessian.

####  4.5.4 Damped Newton

“Not moving too far away” could also be phrased as a line search along the
direction given by the Newton-Raphson method, i.e.

$\displaystyle
h_{n}:=\arg\min_{h}\mathcal{L}\left(w_{n}-h[\nabla^{2}\mathcal{L}(w_{n})]^{-1}\nabla\mathcal{L}(w_{n})\right)$

resulting in

$\displaystyle
w_{n+1}=w_{n}-h_{n}[\nabla^{2}\mathcal{L}(w_{n})]^{-1}\nabla\mathcal{L}(w_{n}).$
(damped Newton)

###  4.6 Gauß-Newton Algorithm

The Gauß-Newton Algorithm opens up the black box of our loss function, to
create a relatively cheap approximation of the Hessian
$\nabla^{2}\mathcal{L}$. To motivate the (classical) Gauß-Newton algorithm, we
are going to assume a squared loss

$\displaystyle l(w,z):=\tfrac{1}{2}(f_{w}(x)-y)^{2},\quad z=(x,y).$

Assuming differentiation and integration can be freely swapped, we have

$\displaystyle\nabla\mathcal{L}(w)$
$\displaystyle=\mathbb{E}[\nabla\tfrac{1}{2}(f_{w}(X)-Y)^{2}]=\mathbb{E}[(f_{w}(X)-Y)\nabla
f_{w}(X)],$ (4.13) $\displaystyle\nabla^{2}\mathcal{L}(w)$
$\displaystyle=\mathbb{E}[\nabla[(f_{w}(X)-Y)\nabla f_{w}(X)]]$ (4.14)
$\displaystyle=\mathbb{E}[\nabla f_{w}(X)\nabla
f_{w}(X)^{T}+(f_{w}(X)-Y)\nabla^{2}f_{w}(X)].$

If our model is close to the truth, i.e. we are close to the minimum, then
$f_{w}(X)-Y$ should be small. Therefore we should approximately have

$\displaystyle\nabla^{2}\mathcal{L}(w)\approx\mathbb{E}[\nabla f_{w}(X)\nabla
f_{w}(X)^{T}].$

A different way to motivate this approximation, is to use the first taylor
approximation of $f_{w}$ inside the loss function, i.e.

$\displaystyle\mathcal{L}(\theta)$
$\displaystyle=\mathbb{E}[l(\theta,Z)]\approx\mathbb{E}[\tfrac{1}{2}(f_{w}(X)+\langle\nabla
f_{w}(X),\theta-w\rangle-Y)^{2}]$ $\displaystyle=\begin{aligned}
&\mathbb{E}[\tfrac{1}{2}(f_{w}(X)-Y)^{2}]\\\
&+\mathbb{E}[(f_{w}(X)-Y)\langle\nabla f_{w}(X),\theta-w\rangle]\\\
&+\mathbb{E}[\tfrac{1}{2}\langle\nabla
f_{w}(X),\theta-w\rangle^{2}]\end{aligned}$
$\displaystyle=\mathcal{L}(w)+\langle\nabla\mathcal{L}(w),\theta-w\rangle+\tfrac{1}{2}(\theta-w)^{T}\underbrace{\mathbb{E}[\nabla
f_{w}(X)\nabla f_{w}(X)^{T}]}_{\approx\nabla^{2}\mathcal{L}(w)}(\theta-w)$

The nice thing about this approximation is, that it is always positive
definite since we have

$\displaystyle v^{T}\mathbb{E}[\nabla f_{w}(X)\nabla f_{w}(X)^{T}]v$
$\displaystyle=\mathbb{E}[v^{T}\nabla f_{w}(X)\nabla f_{w}(X)^{T}v]$
$\displaystyle=\mathbb{E}[\langle\nabla f_{w}(X),v\rangle^{2}]\geq 0.$

####  4.6.1 Levenberg-Marquard Algorithm

Unfortunately, it is not necessarily strictly positive definite. That is, if
there is a vector $v$, which is orthogonal to $\nabla f_{w}(X)$ for all but a
zero probability set of $X$, then this vector is in the eigenspace with
eigenvalue zero. For the expectation one might still be able to discuss that
away. But if one optimizes with regard to a finite sample
$\mathcal{S}=(Z_{1},\dots,Z_{S})$, then we have

$\displaystyle\nabla^{2}\mathcal{L}_{\mathcal{S}}(w)\approx\mathbb{E}_{Z\sim\mathcal{D}_{\mathcal{S}}}[\langle\nabla
f_{w}(X),v\rangle^{2}]=\frac{1}{S}\sum_{k=1}^{S}[\langle\nabla
f_{w}(X_{k}),v\rangle^{2}]$

and it is much easier to find a vector orthogonal to a finite set of
gradients. Especially if we are talking about small mini-batches and do not
use full batch gradient descent.

For this reason the Gauß-Newton algorithm is usually combined with Levenberg-
Marquard regularization (Subsection 4.5.2) to allow for inversion. With this
regularization it is also called Levenberg-Marquard algorithm going back to
[42] and rediscovered by [45], since that regularization technique was
developed with the Gauß-Newton algorithm in mind.

####  4.6.2 Generalized Gauß-Newton

While the classical Gauß-Newton algorithm seems to be tailored to squared
losses only, this approach can actually be generalized to more general loss
functions [65]. Assuming the loss function is of the form $l(f,y)$, i.e with
some abuse of notation

$\displaystyle l(w,z):=l(f_{w}(x),y),\quad z=(x,y),$

then we obtain similarly to (4.13) and (4.14)

$\displaystyle\tfrac{d}{dw}\mathcal{L}(w)$
$\displaystyle=\mathbb{E}[\tfrac{d}{dw}(f_{w}(X),Y)]=\mathbb{E}[\tfrac{d}{df}l(f_{w}(X),Y)\tfrac{d}{dw}f_{w}(X)],$
(4.15) $\displaystyle\tfrac{d^{2}}{dw^{2}}\mathcal{L}(w)$
$\displaystyle=\mathbb{E}[\tfrac{d}{dw}[\tfrac{d}{df}l(f_{w}(X),Y)\tfrac{d}{dw}f_{w}(X)]]$
(4.16) $\displaystyle=\begin{aligned}
&\mathbb{E}[\tfrac{d}{dw}f_{w}(X)^{T}\tfrac{d^{2}}{df^{2}}l(f_{w}(X),Y)\tfrac{d}{dw}f_{w}(X)]\\\
&+\mathbb{E}[\tfrac{d}{df}l(f_{w}(X),Y)\tfrac{d^{2}}{dw^{2}}f_{w}(X)].\end{aligned}$

Again we throw away the part with the Hessian of the model $f_{w}$, to obtain
the approximation

$\displaystyle\nabla^{2}\mathcal{L}(w)$
$\displaystyle\approx\mathbb{E}[\tfrac{d}{dw}f_{w}(X)\tfrac{d^{2}}{df^{2}}^{T}l(f_{w}(X),Y)\tfrac{d}{dw}f_{w}(X)].$

The result is quite similar to the classical Gauß-Newton method, except for
the Hessian of the loss function in the middle. The calculation of which might
be a huge or tiny problem. If the model prediction is one dimensional it is
likely not an issue. Then it would only be a constant.

One should still come up with some justification, why the second part of the
Hessian part can be thrown away though, i.e. why is

$\displaystyle\mathbb{E}[\tfrac{d}{df}l(f_{w}(X),Y)\tfrac{d^{2}}{dw^{2}}f_{w}(X)]$

small? In the sqared loss case, we have used the fact that
$\tfrac{d}{df}l(f_{w}(X),Y)$ was small for good models.

###  4.7 Natural Gradient Descent

“Natural gradient descent” would require a deep dive into “information
geometry” for which we do not have time. But not mentioning this algorithm at
all would do it injustice. So to avoid playing a game of telephone, here is a
direct quote of the description by [5]:

“We have seen that Newton’s method is invariant to linear transformations of
the parameter vector w. By contrast, the natural gradient method [5, 6] aims
to be invariant with respect to all differentiable and invertible
transformations. The essential idea consists of formulating the gradient
descent algorithm in the space of prediction functions rather than specific
parameters. Of course, the actual computation takes place with respect to the
parameters, but accounts for the anisotropic relation between the parameters
and the decision function. That is, in parameter space, the natural gradient
algorithm will move the parameters more quickly along directions that have a
small impact on the decision function, and more cautiously along directions
that have a large impact on the decision function.

We remark at the outset that many authors [119, 99] propose quasi-natural-
gradient methods that are strikingly similar to the quasi-Newton methods
described in §6.2. The natural gradient approach therefore offers a different
justification for these algorithms, one that involves qualitatively different
approximations. It should also be noted that research on the design of methods
inspired by the natural gradient is ongoing and may lead to markedly different
algorithms.”

###  4.8 Coordinate Descent

Coordinate descent only ever updates one coordinate. “[I]n situations in which
$d$ coordinate updates can be performed at cost similar to the evaluation of
one full gradient, the method is competitive with a full gradient method both
theoretically and in practice.” [5, 7.3, p. 72]. This might be interesting as
the Hessian is of course much easier to handle in one dimension. A line search
might not be the right thing to though, as [58] provides an example of a three
dimensional function on which coordinate descent jumps between the corners of
the cube, never converging. This behavior can be fixed with a stochastic
selection of the coordinates, instead of cycling through them. In fact [10]
argues that coordinate descent _is_ just SGD with a random filter, filtering
out all but one coordinate and disappearing in expectation. I.e. for a random
index $I$ and standard basis vectors $e_{1},\dots,e_{d}$ we can define

$\displaystyle\nabla
l(w,i):=\tfrac{1}{P(I=i)}\langle\nabla\mathcal{L}(w),e_{i}\rangle e_{i}$

to get the original loss in expectation

$\displaystyle\mathbb{E}[\nabla l(w,I)]=\sum_{k=1}^{d}P(I=k)\nabla
l(w,k)=\sum_{k=1}^{d}\langle\nabla\mathcal{L}(w),e_{k}\rangle
e_{k}=\nabla\mathcal{L}(w).$

With conditional independence the same argument can be applied to already
stochastic losses. So we can apply our SGD convergence statements to prove
convergence of random coordinate descent.

##  Appendix A Technical Proofs

###  A.1 Non-Trivial Basics

####  A.1.1 Convex Analysis

######  Lemma A.1.1.

If $f$ is differentiable, then the derivative $\nabla f$ is bounded (w.r.t.
the operator norm) by constant $K$ if and only if the function $f$ is
$K$-Lipschitz continuous.

###### Proof.

Lipschitz continuity is implied by the mean value theorem

$\displaystyle\|f(x_{1})-f(x_{0})\|$ $\displaystyle\leq\|\nabla
f(x_{0}+\xi(x_{1}-x_{0}))\|\|x_{1}-x_{0}\|$ $\displaystyle\leq
K\|x_{1}-x_{0}\|.$

The opposite direction is implied by

$\displaystyle\|\nabla f(x_{0})\|$ $\displaystyle\equiv\sup_{v}\frac{\|\nabla
f(x_{0})v\|}{\|v\|}=\sup_{v}\lim_{\lambda\to 0}\frac{|\lambda|\|\nabla
f(x_{0})v\|}{|\lambda|\|v\|}$
$\displaystyle\mathrel{\smashoperator[l]{\mathop{\leq}_{}^{\Delta}}}\sup_{v}\lim_{\lambda\to
0}(\underbrace{\tfrac{\|\nabla f(x_{0})\lambda v+f(x_{0})-f(x_{0}+\lambda
v)\|}{\|\lambda v\|}}_{\to 0\text{ (derivative
definition)}}+\underbrace{\tfrac{\|f(x_{0}+\lambda v)-f(x_{0})\|}{\|\lambda
v\|}}_{\leq K})$

where we have used the scalability of the norm to multiply $v$ with a
decreasing factor both in the numerator and denominator in order to introduce
the limit. ∎

######  Lemma A.1.2.

If $\nabla f$ is $L$-Lipschitz continuous, then

$\displaystyle|f(y)-f(x)-\langle\nabla
f(x),y-x\rangle|\leq\tfrac{L}{2}\|y-x\|^{2}$

If $f$ is convex, then the opposite direction is also true.

###### Proof.

The first direction is taken from [50, Lemma 1.2.3].

$\displaystyle f(y)=f(x)+\int_{0}^{1}\langle\nabla f(x+\tau(y-x)),y-x\rangle
d\tau$

implies (using the Cauchy-Schwarz inequality)

$\displaystyle|f(y)-f(x)-\langle\nabla f(x),y-x\rangle|$
$\displaystyle\leq\int_{0}^{1}|\langle\nabla f(x+\tau(y-x))-\nabla
f(x),y-x\rangle|d\tau$
$\displaystyle\mathrel{\smashoperator[l]{\mathop{\leq}_{}^{\text{C.S.}}}}\int_{0}^{1}\|\langle\nabla
f(x+\tau(y-x))-\nabla f(x)\|\cdot\|y-x\|d\tau$
$\displaystyle\leq\int_{0}^{1}L\|\tau(y-x)\|\cdot\|y-x\|d\tau=\tfrac{L}{2}\|y-x\|^{2}.$

The opposite direction is taken from [50, Lemma 2.1.5]. As we have only used

$\displaystyle 0\;\mkern-5.0mu plus
-5.0mu\mathop{\leq}\limits^{\text{Convexity}}\;\mkern-5.0mu plus
-5.0muf(y)-f(x)-\langle\nabla f(x),y-x\rangle\leq\tfrac{L}{2}\|y-x\|^{2}$

in Lemma 1.4.2 (and Lemma 1.3.1 which was used in the proof), we know that
equation (1.27b) follows from it and after applying Cauchy-Schwarz to this
equation

$\displaystyle\tfrac{1}{L}\|\nabla f(x)-\nabla f(y)\|^{2}$
$\displaystyle\leq\langle\nabla f(x)-\nabla f(y),x-y\rangle$
$\displaystyle\mathrel{\smashoperator[l]{\mathop{\leq}_{}^{\text{C.S.}}}}\|\nabla
f(x)-\nabla f(y)\|\|x-y\|,$

we only have to divide both sides by $\tfrac{1}{L}\|\nabla f(x)-\nabla f(y)\|$
to get $L$-Lipschitz continuity back. ∎

####  A.1.2 Sequences

######  Theorem A.1.3 (Cesàro-Stolz).

Let $(x_{n},n\in\mathbb{N}),(y_{n},n\in\mathbb{N})$ be real sequences with
$\lim_{n\to\infty}y_{n}\to\infty$, then for $\Delta x_{n}:=x_{n+1}-x_{n}$

$\displaystyle\liminf_{n}\frac{\Delta x_{n}}{\Delta
y_{n}}\leq\liminf_{n}\frac{x_{n}}{y_{n}}\leq\limsup_{n}\frac{x_{n}}{y_{n}}\leq\limsup_{n}\frac{\Delta
x_{n}}{\Delta y_{n}}$ (A.1)

In particular if the limit $\lim_{n}\frac{\Delta x_{n}}{\Delta y_{n}}$ exists
we have equality. Thereby it becomes a discrete version of L’Hôpital’s rule.

###### Proof [34].

To show the first inequality let us assume we have

$\displaystyle\gamma<\liminf_{n}\frac{\Delta x_{n}}{\Delta y_{n}}.$

Then by definition of the limes inferior there exists $N$ such that

$\displaystyle\gamma\Delta y_{n}<\Delta x_{n}\qquad\forall n\geq N$

Using telescoping sums this implies

$\displaystyle\gamma(y_{m}-y_{N})=\gamma\sum_{n=N}^{m-1}\Delta
y_{n}<\sum_{n=N}^{m-1}\Delta x_{n}=x_{m}-x_{N}.$

Dividing both sides by $y_{m}$ and taking the limes inferior keeping in mind
that $y_{m}\to\infty$ results in

$\displaystyle\gamma=\liminf_{m}\gamma\left(\frac{y_{m}}{y_{m}}-\frac{y_{N}}{y_{m}}\right)<\liminf_{m}\frac{x_{m}}{y_{m}}-\frac{x_{N}}{y_{m}}=\liminf_{m}\frac{x_{m}}{y_{m}}.$

Since $\gamma$ was arbitrary this proves the first inequality. The second
inequality is true by definition of the limes inferior (and superior) and the
third inequality can be proven in the same way as the first starting with
$\limsup_{n}\frac{\Delta x_{n}}{\Delta y_{n}}<\gamma$. ∎

######  Lemma A.1.4.

Let $a_{0}\in[0,1/q)$ for $q>0$ and assume

$\displaystyle a_{n+1}=(1-qa_{n})a_{n}\quad\forall n\geq 0,$

then we have

$\displaystyle\lim_{n\to\infty}na_{n}=1/q$ (A.2)

###### Proof.

The proof follows [35]. By induction we have

$\displaystyle 0\leq a_{n}\leq a_{n-1}<1/q.$

Therefore $a_{n}$ is a falling bounded sequences and thus converges to some
$a\in[0,1/q)$ which necessarily fulfills $a(1-qa)=a$ which implies

$\displaystyle\lim_{n\to\infty}a_{n}=a=0$

Thus for $y_{n}:=1/a_{n}$ and $x_{n}:=n$ we have

$\displaystyle\frac{\Delta x_{n}}{\Delta
y_{n}}=\frac{1}{\frac{1}{a_{n+1}}-\frac{1}{a_{n}}}=\frac{a_{n+1}a_{n}}{a_{n}-a_{n+1}}=\frac{(1-qa_{n})a_{n}}{1-(1-qa_{n})}=\frac{1-qa_{n}}{q}\to
1/q$

and since $y_{n}\to\infty$ we can apply Theorem A.1.3 to get

$\displaystyle\lim_{n\to\infty}na_{n}$
$\displaystyle=\lim_{n\to\infty}\frac{x_{n}}{y_{n}}=\lim_{n\to\infty}\frac{\Delta
x_{n}}{\Delta y_{n}}=1/q.\qed$

######  Lemma A.1.5 (Discrete Grönwall).

Let $g_{n},a_{n},b_{n}\geq 0$, $g_{0}=0$ and

$\displaystyle g_{n+1}\leq a_{n}g_{n}+b_{n}$

then we have

$\displaystyle g_{n}\leq\sum_{k=0}^{n-1}b_{k}\prod_{i=k+1}^{n-1}a_{i}$

###### Proof.

By induction where $n=0$ is trivial interpreting the empty sum as zero, we
have

$\displaystyle g_{n+1}$ $\displaystyle\leq a_{n}g_{n}+b_{n}\;\mkern-5.0mu plus
-5.0mu\mathop{\leq}\limits^{\text{ind.}}\;\mkern-5.0mu plus
-5.0mua_{n}\left(\sum_{k=0}^{n-1}b_{k}\prod_{i=k+1}^{n-1}a_{i}\right)+b_{n}\leq\sum_{k=0}^{n}b_{k}\prod_{i=k+1}^{n}a_{i}\qed$

####  A.1.3 Operators

######  Lemma A.1.6 (Absolute Values Inside Operator Norms).

For complex matrices $A^{(1)},\dots,A^{(n)}\in\mathbb{C}^{d\times d}$ we
define $B^{(1)},\dots,B^{(n)}\in\mathbb{C}^{d\times d}$ with

$\displaystyle B^{(k)}_{ij}:=|A^{(k)}_{ij}|\qquad\forall k,i,j.$

Then we have

$\displaystyle\Bigg{\|}\prod_{k=1}^{n}A^{(k)}\Bigg{\|}\leq\Bigg{\|}\prod_{k=1}^{n}B^{(k)}\Bigg{\|}.$

###### Proof.

Let $x\in\mathbb{C}^{d}$ be an arbitrary vector and define
$y\in\mathbb{C}^{d}$ with

$\displaystyle y_{i}=|x_{i}|\qquad\forall i=1,\dots,d.$

Then we have

$\displaystyle\Bigg{\|}\prod_{k=1}^{n}A^{(k)}x\Bigg{\|}$
$\displaystyle=\sum_{j=1}^{d}\Big{|}\sum_{1\leq i_{1},\dots,i_{n}\leq
d}A^{(1)}_{ji_{1}}A^{(2)}_{i_{1}i_{2}}\dots
A^{(n)}_{i_{n-1}i_{n}}x_{i_{n}}\Big{|}^{2}$
$\displaystyle\leq\sum_{j=1}^{d}\Big{|}\sum_{1\leq i_{1},\dots,i_{n}\leq
d}B^{(1)}_{ji_{1}}B^{(2)}_{i_{1}i_{2}}\dots
B^{(n)}_{i_{n-1}i_{n}}y_{i_{n}}\Big{|}^{2}$
$\displaystyle=\Bigg{\|}\prod_{k=1}^{n}B^{(k)}y\Bigg{\|}$

Since we have $\|x\|=\|y\|$ as the norm takes the absolute value anyway, we
get our claim

$\displaystyle\Bigg{\|}\prod_{k=1}^{n}A^{(k)}\Bigg{\|}$
$\displaystyle=\sup_{\|x\|=1}\Bigg{\|}\prod_{k=1}^{n}A^{(k)}x\Bigg{\|}\leq\sup_{\|y\|=1}\Bigg{\|}\prod_{k=1}^{n}B^{(k)}y\Bigg{\|}=\Bigg{\|}\prod_{k=1}^{n}B^{(k)}\Bigg{\|}.\qed$

###  A.2 Momentum Convergence

####  A.2.1 Eigenvalue Analysis

######  Theorem A.2.1 ([59]).

Let

$\displaystyle
r_{1/2}=\tfrac{1}{2}\left(1+\beta-\alpha\lambda\pm\sqrt{(1+\beta-\alpha\lambda)^{2}-4\beta}\right)$

then

  1. 1.

$\max\\{|r_{1}|,|r_{2}|\\}<1$ if and only if

$\displaystyle 0<\alpha\lambda<2(1+\beta)\qquad\text{and}\qquad|\beta|<1$

  2. 2.

The complex case can be characterized by either

$\displaystyle 0<(1-\sqrt{\alpha\lambda})^{2}<\beta<1$

or alternatively $\beta>0$ and

$\displaystyle(1-\sqrt{\beta})^{2}<\alpha\lambda<(1+\sqrt{\beta})^{2},$

for which we have $|r_{1}|=|r_{2}|=\sqrt{\beta}$.

  3. 3.

In the real case we have $r_{1}>r_{2}$ and

$\displaystyle\max\\{|r_{1}|,|r_{2}|\\}=\begin{cases}|r_{1}|=r_{1}&\alpha\lambda<1+\beta\\\
|r_{2}|=-r_{2}&\alpha\lambda\geq 1+\beta.\end{cases}$ (A.3)

Restricted to $1>\beta>0$ this results in two different behaviors. For

$\displaystyle 0<\alpha\lambda\leq(1-\sqrt{\beta})^{2}<1+\beta$

we have $1>r_{1}>r_{2}>0$. For

$\displaystyle 1+\beta<(1+\sqrt{\beta})^{2}\leq\alpha\lambda<2(1+\beta)$

on the other hand, we get $-1<r_{2}<r_{1}<0$.

###### Proof.

Define

$\displaystyle\Delta(\beta)$
$\displaystyle:=(1+\beta-\alpha\lambda)^{2}-4\beta$
$\displaystyle=\beta^{2}-2(1+\alpha\lambda)\beta+(1-\alpha\lambda)^{2}$

then $r_{1/2}$ is complex iff $\Delta(\beta)<0$. Since it is a convex parabola
this implies that $\beta$ needs to be between the roots

$\displaystyle\beta_{1/2}=(1+\alpha\lambda)\pm\underbrace{\sqrt{(1+\alpha\lambda)^{2}-(1-\alpha\lambda)^{2}}}_{=\sqrt{4\alpha\lambda}=\mathrlap{2\sqrt{\alpha\lambda}}}=(1\pm\sqrt{\alpha\lambda})^{2}$

of $\Delta$. Assuming those roots are real, which requires $\alpha\lambda>0$.
But if $\alpha\lambda\leq 0$ then we would have

$\displaystyle
r_{1}=\tfrac{1}{2}\Big{(}\underbrace{1+\beta+|\alpha\lambda|}_{\geq
1+\beta}+\sqrt{\underbrace{(1+\beta+|\alpha\lambda|)^{2}-4\beta}_{\smash{\geq(1-\beta)^{2}}}}\Big{)}\geq\frac{1+\beta+|1-\beta|}{2}\geq
1$

which means that $\max\\{|r_{1}|,|r_{2}|\\}<1$ implies $\alpha\lambda>0$. The
roots $\beta_{1/2}$ are therefore real no matter the direction we wish to
prove.

Complex Case:

$r_{1/2}$ are complex iff

$\displaystyle
0<(1-\sqrt{\alpha\lambda})^{2}<\beta<(1+\sqrt{\alpha\lambda})^{2}.$ (A.4)

In that case we have

$\displaystyle|r_{1}|=|r_{2}|$
$\displaystyle=\sqrt{|\Re(r_{1})|^{2}+|\Im(r_{1})|^{2}}$
$\displaystyle=\tfrac{1}{2}\sqrt{(1+\beta-\alpha\lambda)^{2}+4\beta-(1+\beta-\alpha\lambda)^{2}}=\sqrt{\beta}.$

So the condition $\beta<1$ is necessary and sufficient for
$\max\\{|r_{1}|,|r_{2}|\\}<1$ in the complex case. And since we have seen that
$\alpha\lambda>0$ is necessary and $0<\beta$ from (A.4) implies the condition
$|\beta|<1$, we only need to show that $\alpha\lambda<2(1+\beta)$ is necessary
in the complex case to have this case covered. Using
$(\sqrt{\alpha\lambda}-1)^{2}<\beta$ from (A.4) and $ab\leq a^{2}+b^{2}$ we
get

$\displaystyle\sqrt{\alpha\lambda}-1<\beta\implies\alpha\lambda<(1+\sqrt{\beta})^{2}\leq
2(1+\beta).$

The remaining characterization of the complex case follows from restricting
(A.4) to $\beta<1$ and reordering.

Real Case:

Since we have $r_{2}\leq r_{1}$, $\max\\{|r_{1}|,|r_{2}|\\}<1$ is equivalent
to

$\displaystyle-1<r_{2}\qquad\text{and}\qquad r_{1}<1$

By subtracting the part before the ”$\pm$” from the equation $r_{1}<1$ after
multiplying it by two, we get

$\displaystyle 0\;\mkern-5.0mu plus -5.0mu\mathop{\leq}\limits^{\text{Real
Case}}\;\mkern-5.0mu plus
-5.0mu\sqrt{(1+\beta-\alpha\lambda)^{2}-4\beta}<1-\beta+\alpha\lambda.$ (A.5)

And since $(1+\sqrt{\alpha\lambda})^{2}\leq\beta$ leads to a contradiction

$\displaystyle 0\leq 1-\beta+\alpha\lambda\leq
1-(1+\sqrt{\alpha\lambda})^{2}+\alpha\lambda=-2\sqrt{\alpha\lambda}\geq 0,$

the real case is restricted to the case $\beta<(1-\sqrt{\alpha\lambda})^{2}$
as we are otherwise in the complex case (A.4). The inequality (A.5) is (in the
real case) equivalent to

$\displaystyle(1+\beta-\alpha\lambda)^{2}-4\beta<(1-\beta+\alpha\lambda)^{2}$
$\displaystyle\iff\cancel{1^{2}}+2(\beta-\alpha\lambda)+\cancel{(\beta-\alpha\lambda)^{2}}-4\beta<\cancel{1^{2}}-2(\beta-\alpha\lambda)+\cancel{(\beta-\alpha\lambda)^{2}}$
$\displaystyle\iff 0<4\alpha\lambda.$

Which is a given for positive learning rates in the convex case. This is not
surprising if one remembers the proof from gradient descent. The real issue is
not overshooting $-1$ with large learning rates. So multiplying $-1<r_{2}$ by
two, adding 2 and moving the root to the other side we get

$\displaystyle 0\;\mkern-5.0mu plus -5.0mu\mathop{\leq}\limits^{\text{Real
Case}}\;\mkern-5.0mu plus
-5.0mu\sqrt{(1+\beta-\alpha\lambda)^{2}-4\beta}<3+\beta-\alpha\lambda.$ (A.6)

From this we get $\alpha\lambda<3+\beta$. We also get

$\displaystyle(1+\beta-\alpha\lambda)^{2}-4\beta<(3+\beta-\alpha\lambda)^{2}$
$\displaystyle\iff
1^{2}+2(\beta-\alpha\lambda)+\cancel{(\beta-\alpha\lambda)^{2}}-4\beta<3^{2}+6(\beta-\alpha\lambda)+\cancel{(\beta-\alpha\lambda)^{2}}$
$\displaystyle\iff 4\alpha\lambda<8+8\beta\iff\alpha\lambda<2(1+\beta)$ (A.7)

which is the important requirement on the learning rate we are looking for.
Now for $\beta<1$ this requirement is actually stronger than
$\alpha\lambda<3+\beta$ which means that we have proven that our requirements
are sufficient for $\max\\{|r_{1}|,|r_{2}|\\}<1$. What is left to show is that
$|\beta|<1$ is also necessary in the real case. From (A.7) we immediately get

$\displaystyle 0\leq\tfrac{\alpha\lambda}{2}<1+\beta$

which implies $-1<\beta$. At first glance it might look like
$\beta<(1-\sqrt{\alpha\lambda})^{2}$ is already sufficient for the upper
bound, but for $\alpha\lambda>4$ this bound is greater than one again. So let
us assume $\alpha\lambda>2$ and $\beta>1$, then
$\beta<(1-\sqrt{\alpha\lambda})^{2}$ implies a contradiction:

$\displaystyle\sqrt{\beta}<\sqrt{\alpha\lambda}-1$
$\displaystyle\implies\underbrace{(1+\sqrt{\beta})^{2}}_{1+2\sqrt{\beta}+\beta}<\alpha\lambda\;\mkern-5.0mu
plus -5.0mu\mathop{<}\limits^{(\ref{eq: from -1<sigma_2
followed})}\;\mkern-5.0mu plus -5.0mu3+\beta$
$\displaystyle\implies-2+2\sqrt{\beta}<0\implies\sqrt{\beta}<1$

For the _characterization of the real case_ one only needs to observe that for
any $b>0$

$\displaystyle\max|a\pm b|=\begin{cases}a+b&a\geq 0\\\ -(a-b)&a<0\end{cases}.$
(A.8)

In our case $b$ is a square root in the real case, so by definition positive
and $a>0$ is equivalent to $1+\beta-\alpha\lambda>0$ and thus implies (A.3).

For $\beta>0$ one only needs to keep in mind that for
$\alpha\lambda<(1+\sqrt{\beta})^{2}<1+\beta$

$\displaystyle
1+\beta-\alpha\lambda=\sqrt{(1+\beta-\alpha\lambda)^{2}}>\sqrt{(1+\beta-\alpha\lambda)^{2}-4\beta}$

ensures that $r_{2}>0$. On the other hand for $\beta>0$ and
$\alpha\lambda>(1-\sqrt{\beta})^{2}>1+\beta$

$\displaystyle-(1+\beta-\alpha\lambda)=\sqrt{(1+\beta-\alpha\lambda)^{2}}>\sqrt{(1+\beta-\alpha\lambda)^{2}-4\beta}$

implies $r_{1}<0$. ∎

######  Lemma A.2.2.

$\displaystyle\frac{d}{d\beta}\max\\{|r_{1}|,|r_{2}|\\}<0\qquad\forall-1<\beta<(1-\sqrt{\alpha\lambda})^{2}$

###### Proof.

Using the characterizations of the real case from Theorem 3.1.1, we get

$\displaystyle\frac{d}{d\beta}\max\\{|r_{1}|,|r_{2}|\\}=\begin{cases}\frac{dr_{1}}{d\beta}&\beta\geq\alpha\lambda-1\\\
-\frac{dr_{2}}{d\beta}&\beta\leq\alpha\lambda-1\end{cases}$

This allows us to treat those cases separately. For the first case we use the
same definition for $\Delta$ as in Theorem 3.1.1. Calculating its derivative

$\displaystyle\Delta^{\prime}(\beta)=2(1+\beta-\alpha\lambda)-4$

allows us to represent $r_{1}$ as

$\displaystyle
1>r_{1}=1+\frac{\Delta^{\prime}(\beta)}{4}+\tfrac{1}{2}\sqrt{\Delta(\beta)},$

which results in the inequality

$\displaystyle
0>\tfrac{1}{4}\Delta^{\prime}(\beta)+\tfrac{1}{2}\sqrt{\Delta(\beta)}\implies-\tfrac{1}{2}>\tfrac{1}{4}\frac{\Delta^{\prime}(\beta)}{\sqrt{\Delta{\beta}}}$

We can now use this inequality to show that111 I want to apologize for the
lack of insight in this proof upon which I stumbled by accident playing around
with the $\Delta(\beta)$ representation of $r_{1/2}$. Its only redeeming merit
is its brevity.

$\displaystyle\frac{dr_{1}}{d\beta}=\frac{1}{2}+\frac{1}{4}\frac{\Delta^{\prime}(\beta)}{\sqrt{\Delta(\beta)}}<0.$

The second case on the other hand requires us to show that

$\displaystyle\frac{dr_{2}}{d\beta}=\frac{1}{2}-\frac{1}{4}\frac{\Delta^{\prime}(\beta)}{\sqrt{\Delta(\beta)}}>0\iff
2>\frac{\Delta^{\prime}(\beta)}{\sqrt{\Delta(\beta)}}.$

Which is true due to

$\displaystyle 2\sqrt{\Delta(\beta)}>0>-4\sqrt{\alpha\lambda}$
$\displaystyle=2\overbrace{(1-\alpha\lambda+(1-\sqrt{\alpha\lambda})^{2})}^{\smash{=2-2\sqrt{\alpha\lambda}}}-4$
$\displaystyle=\Delta^{\prime}((1-\sqrt{\alpha\lambda})^{2})>\Delta^{\prime}(\beta),$

since $\beta<(1-\sqrt{\alpha\lambda})^{2}$ and $\Delta^{\prime}$ is
monotonously increasing in $\beta$. ∎

####  A.2.2 Nesterov’s Momentum Convergence

######  Lemma A.2.3 (Dumpster Dive).

The conditions

$\displaystyle L\gamma_{n}^{2}\leq\mu_{n}$

and

$\displaystyle y_{n}$
$\displaystyle=\tfrac{\mu_{n}w_{n-1}+\mu_{n-1}\gamma_{n}z_{n-1}}{\gamma_{n}\mu+\mu_{n-1}}$

are sufficient for lower bounding

$\displaystyle
J:=(1-\gamma_{n})\left(\langle\mathcal{L}(y_{n}),w_{n-1}-y_{n}\rangle+\tfrac{\mu_{n-1}}{2}\|y_{n}-z_{n-1}\|^{2}\right)-\tfrac{\mu_{n}}{2}\|y_{n}-z_{n}\|^{2}$

by $-\tfrac{1}{2L}\|\nabla\mathcal{L}(y_{n})\|^{2}$.

###### Proof.

By definition $z_{n}$ is the unique minimum of $\Phi_{n}$ so we can find it
with

$\displaystyle 0\;\mkern-5.0mu plus -5.0mu\mathop{=}\limits^{!}\;\mkern-5.0mu
plus
-5.0mu\Phi^{\prime}_{n}(z_{n})=\gamma_{n}\underbrace{\phi_{n}^{\prime}(z_{n})}_{=\nabla\mathcal{L}(y_{n})\mathrlap{+\mu(z_{n}-y_{n})}}+(1-\gamma_{n})\underbrace{\Phi^{\prime}_{n-1}(z_{n})}_{=\mu_{n-1}\mathrlap{(z_{n}-z_{n-1})}},$

where we have used the representation (3.20) for $\Phi_{n-1}$. It is also
quite reassuring that the gradient at $y_{n}$ which we ultimately need makes
its first appearance here. Collecting all the $z_{n}$ we get

$\displaystyle\overbrace{\mu_{n}}^{=\mathrlap{\gamma_{n}\mu+(1-\gamma_{n})\mu_{n-1}}}z_{n}=\gamma_{n}\mu
y_{n}+(1-\gamma_{n})\mu_{n}z_{n-1}-\gamma_{n}\nabla\mathcal{L}(y_{n}).$ (A.9)

Subtracting $\mu_{n}y_{n}$ results in

$\displaystyle\mu_{n}(z_{n}-y_{n})=(1-\gamma_{n})\mu_{n-1}(z_{n-1}-y_{n})-\gamma_{n}\nabla\mathcal{L}(y_{n}).$

Plugging this into the last part of our “junk”, we get

$\displaystyle\tfrac{\mu_{n}}{2}\|y_{n}-z_{n}\|^{2}$
$\displaystyle=\tfrac{1}{2\mu_{n}}\underbrace{\|{\scriptstyle(1-\gamma_{n})\mu_{n-1}}(z_{n-1}-y_{n})-{\scriptstyle\gamma_{n}}\nabla\mathcal{L}(y_{n})\|^{2}}_{\begin{aligned}
=&{\scriptstyle(1-\gamma_{n})^{2}\mu_{n-1}^{2}}\|z_{n-1}-y_{n}\|^{2}\\\
&-{\scriptstyle
2(1-\gamma_{n})\mu_{n-1}\gamma_{n}}\langle\nabla\mathcal{L}(y_{n}),z_{n-1}-y_{n}\rangle\\\
&+{\scriptstyle\gamma_{n}^{2}}\|\nabla\mathcal{L}(y_{n})\|^{2}\end{aligned}}$

Therefore our “junk” is equal to

$\displaystyle-\tfrac{\gamma_{n}^{2}}{2\mu_{n}}\|\nabla\mathcal{L}(y_{n})\|^{2}\begin{aligned}
&+{\scriptstyle(1-\gamma_{n})}\left\langle\mathcal{L}(y_{n}),(w_{n}-y_{n})+\tfrac{\mu_{n-1}\gamma_{n}}{\mu_{n}}(z_{n-1}-y_{n})\right\rangle\\\
&+\tfrac{(1-\gamma_{n})\mu_{n-1}}{2}\underbrace{\left(1-\tfrac{(1-\gamma_{n})\mu_{n-1}}{\mu_{n}}\right)}_{=\tfrac{\gamma_{n}\mu}{\mu_{n}}\geq
0}\|y_{n}-z_{n-1}\|^{2}.\end{aligned}$

As the last summand is positive we can discard it and end up with two
equations sufficient to achieve our lower bound

$\displaystyle-\tfrac{\gamma_{n}^{2}}{\mu_{n}}$
$\displaystyle\geq-\tfrac{1}{L}\iff L\gamma_{n}^{2}\leq\mu_{n}$ $\displaystyle
0$
$\displaystyle=(w_{n-1}-y_{n})+\tfrac{\mu_{n}\gamma_{n}}{\mu_{n}}(z_{n-1}-y_{n})\;\mkern-5.0mu
plus -5.0mu\mathop{\iff}\limits^{\text{reorder}}\;\mkern-5.0mu plus
-5.0muy_{n}=\tfrac{\mu_{n}w_{n-1}+\mu_{n-1}\gamma_{n}z_{n-1}}{\gamma_{n}\mu+\mu_{n-1}}\qed$

######  Lemma A.2.4.

If we use the recursion

$\displaystyle w_{n}=y_{n}-\tfrac{1}{L}\mathcal{L}(y_{n})$ (A.10)

in Algorithm 1 line 1, then we can simplify $z_{n}$ and $y_{n}$ to

$\displaystyle z_{n}$
$\displaystyle=w_{n-1}+\tfrac{1}{\gamma_{n}}(w_{n}-w_{n-1})\qquad\forall n\geq
1$ (A.11) $\displaystyle y_{n+1}$
$\displaystyle=w_{n}+\beta_{n}(w_{n}-w_{n-1})\qquad\forall n\geq 0$ (A.12)

for $w_{-1}:=w_{0}=z_{0}$ and

$\displaystyle\beta_{n}=\frac{\mu_{n}\gamma_{n+1}(1-\gamma_{n})}{\gamma_{n}(\gamma_{n+1}\mu+\mu_{n})}=\frac{\gamma_{n}(1-\gamma_{n})}{\gamma_{n}^{2}+\gamma_{n+1}}$
(A.13)

###### Proof.

After the assignment in line 1 we can reorder it as an equation to get

$\displaystyle\mu_{n-1}z_{n-1}=\tfrac{1}{\gamma_{n}}[(\gamma_{n}\mu+\mu_{n-1})y_{n}-\mu_{n}w_{n-1}]$

plugging this into our definition of $z_{n}$ results in our first claim

$\displaystyle z_{n}$ $\displaystyle=\tfrac{1}{\mu_{n}}\left[\gamma_{n}\mu
y_{n}+\tfrac{(1-\gamma_{n})}{\gamma_{n}}[(\gamma_{n}\mu+\mu_{n-1})y_{n}-\mu_{n}w_{n-1}]-\gamma_{n}\nabla\mathcal{L}(y_{n})\right]$
$\displaystyle=\tfrac{1}{\mu_{n}}\underbrace{\Big{[}\underbrace{\gamma_{n}\mu+(1-\gamma_{n})\mu}_{=\mu}+\tfrac{1-\gamma_{n}}{\gamma_{n}}\mu_{n-1}\Big{]}}_{=\mu_{n}/\gamma_{n}}y_{n}-\tfrac{1-\gamma_{n}}{\gamma_{n}}w_{n-1}-\underbrace{\tfrac{\gamma_{n}}{\mu_{n}}}_{=\frac{1}{L\gamma_{n}}\mathrlap{\
(L\gamma_{n}^{2}=\mu_{n})}}\nabla\mathcal{L}(y_{n})$
$\displaystyle=w_{n-1}+\tfrac{1}{\gamma_{n}}(y_{n}-w_{n-1})-\tfrac{1}{L\gamma_{n}}\nabla\mathcal{L}(y_{n})$
$\displaystyle\mathrel{\smashoperator[l]{\mathop{=}_{}^{(\ref{eq-lem-assmpt:
weigth recursion})}}}w_{n-1}+\tfrac{1}{\gamma_{n}}(w_{n}-w_{n-1})$

Now we will rewrite $y_{n+1}$:

$\displaystyle y_{n+1}$
$\displaystyle=\frac{\overbrace{[\gamma_{n+1}\mu+(1-\gamma_{n+1})\mu_{n}]}^{\mu_{n+1}}w_{n}+\mu_{n}\gamma_{n+1}z_{n}}{\gamma_{n+1}\mu+\mu_{n}}=w_{n}+\frac{\mu_{n}\gamma_{n+1}(z_{n}-w_{n})}{\gamma_{n+1}\mu+\mu_{n}}$

For $n\geq 1$ we can apply our previous result (A.11)

$\displaystyle
z_{n}-w_{n}=\left(\tfrac{1}{\gamma_{n}}-1\right)(w_{n}-w_{n-1}),$

for $n=0$ both sides are zero and we therefore obtain our second recursion

$\displaystyle
y_{n+1}=w_{n}+\underbrace{\frac{\mu_{n}\gamma_{n+1}(1-\gamma_{n})}{\gamma_{n}(\gamma_{n+1}\mu+\mu_{n})}}_{=\beta_{n}}(w_{n}-w_{n-1}).$

Now we just use the identity
$\gamma_{n+1}\mu=L\gamma_{n+1}^{2}-(1-\gamma_{n+1})\mu_{n}$ from the
definition of $\gamma_{n+1}$ to simplify our momentum coefficient

$\displaystyle\beta_{n}$
$\displaystyle=\frac{\mu_{n}\gamma_{n+1}(1-\gamma_{n})}{\gamma_{n}([L\gamma_{n+1}^{2}-(\cancel{1}-\gamma_{n+1})\mu_{n}]+\cancel{\mu_{n}})}=\frac{\mu_{n}\cancel{\gamma_{n+1}}(1-\gamma_{n})}{\gamma_{n}\cancel{\gamma_{n+1}}(L\gamma_{n+1}+\underbrace{\mu_{n}}_{=L\gamma_{n}^{2}})}$
$\displaystyle\mathrel{\smashoperator[l]{\mathop{=}_{}^{\bcancel{L\gamma_{n}}}}}\frac{\gamma_{n}(1-\gamma_{n})}{\gamma_{n+1}+\gamma_{n}^{2}}\qed$

######  Lemma A.2.5 ([50, Lemma 2.2.4]).

For $\kappa_{0}^{-1}\in(\kappa^{-1},3+\kappa^{-1}]$ we have

$\displaystyle\Gamma^{n}\leq\frac{4\kappa^{-1}}{(\kappa_{0}^{-1}-\kappa^{-1})\left[\exp\left(\frac{n+1}{2\sqrt{\kappa}}\right)-\exp\left(-\frac{n+1}{2\sqrt{\kappa}}\right)\right]^{2}}\leq\frac{4}{(\kappa_{0}^{-1}-\kappa^{-1})(n+1)^{2}}$

###### Proof.

We are interested in finding a function fulfilling

$\displaystyle\Gamma^{n}\leq\frac{1}{f(n)^{2}}\iff
f(n)\leq\frac{1}{\sqrt{\Gamma^{n}}}$ (A.14)

Since the weaker inequality would require $f(n)\sim n$ it might be a good idea
to look at the increments. In particular it would be sufficient if we had
$f(0)\leq\frac{1}{\sqrt{\Gamma^{0}}}=1$ and

$\displaystyle
f(k)-f(k-1)\leq\frac{1}{\sqrt{\Gamma^{k}}}-\frac{1}{\sqrt{\Gamma^{k-1}}}\qquad\forall
k>0.$

For our weaker inequality this requires the right side to be lower bounded by
a constant. Now if we had not started out looking for a squared $f$, we would
have instead ended up with the increment

$\displaystyle\frac{1}{\Gamma^{k}}-\frac{1}{\Gamma^{k-1}}.$

This provides some intuition how Nesterov might have come up with the
following idea. Let us consider the increment above and multiply it by
$\Gamma^{k}$. Then we get

$\displaystyle 1-\frac{\Gamma^{k}}{\Gamma^{k-1}}$
$\displaystyle\mathrel{\smashoperator[l]{\mathop{=}_{}^{\text{def.}}}}\gamma_{k}\;\mkern-5.0mu
plus -5.0mu\mathop{=}\limits^{(\ref{eq: morphing speed
equation})}\;\mkern-5.0mu plus -5.0mu\sqrt{\frac{\mu_{k}}{L}}\;\mkern-5.0mu
plus -5.0mu\mathop{=}\limits^{(\ref{eq: lbound recursion})}\;\mkern-5.0mu plus
-5.0mu\sqrt{\frac{\Gamma^{k}\mu_{0}+(1-\Gamma^{k})\mu}{L}}$
$\displaystyle\mathrel{\smashoperator[l]{\mathop{=}_{}^{\text{def.}}}}\sqrt{\kappa^{-1}+\Gamma^{k}(\kappa_{0}^{-1}-\kappa^{-1})}.$

Now we not only need to return to the increments by dividing by $\Gamma^{k}$
again, we also need to somehow add roots. To do that we use the third binomial
formula

$\displaystyle\tfrac{1}{\sqrt{\Gamma^{k}}}\sqrt{\tfrac{\kappa^{-1}}{\Gamma^{k}}+\kappa_{0}^{-1}-\kappa^{-1}}$
$\displaystyle=\frac{1}{\Gamma^{k}}-\frac{1}{\Gamma^{k-1}}$
$\displaystyle\mathrel{\smashoperator[l]{\mathop{=}_{}^{\text{bin.}}}}\left(\frac{1}{\sqrt{\Gamma^{k}}}+\frac{1}{\sqrt{\Gamma^{k-1}}}\right)\left(\frac{1}{\sqrt{\Gamma^{k}}}-\frac{1}{\sqrt{\Gamma^{k-1}}}\right)$
$\displaystyle\leq\frac{2}{\sqrt{\Gamma^{k}}}\left(\frac{1}{\sqrt{\Gamma^{k}}}-\frac{1}{\sqrt{\Gamma^{k-1}}}\right)$
(A.15)

where we have used $\Gamma^{k}\leq\Gamma^{k-1}$ for the inequality. From this
we get the sufficient condition

$\displaystyle
f(k)-f(k-1)\leq\tfrac{1}{2}\sqrt{\tfrac{\kappa^{-1}}{\Gamma^{k}}+\kappa_{0}^{-1}-\kappa^{-1}}=\tfrac{1}{2}\underbrace{\sqrt{\kappa_{0}^{-1}-\kappa^{-1}}}_{=:c}\sqrt{\tfrac{\kappa^{-1}}{\Gamma^{k}}c^{-2}+1}$

In particular we could throw away the complicated positive part in the root to
get the sufficient condition

$\displaystyle f(k)-f(k-1):=\tfrac{c}{2}\quad\forall k>0\quad\text{and}\quad
f(0):=\tfrac{c}{2}=\tfrac{\sqrt{\kappa_{0}^{-1}-\kappa^{-1}}}{2}\leq\tfrac{\sqrt{3}}{2}<1.$

And this immediately leads to $f(n)=\tfrac{c}{2}(n+1)$ which provides us with
the weaker second bound. Note that the “complicated” part we have thrown away
is zero in case of $\kappa^{-1}$ i.e. when we do not have strong convexity
this is the best we can do.

But for strong convexity we need to find a tighter upper bound $f$. If we
shuffle around (A.15) a bit we get

$\displaystyle\frac{1}{\sqrt{\Gamma^{k-1}}}\leq\frac{1}{\sqrt{\Gamma^{k}}}-\frac{c}{2}\sqrt{\tfrac{\kappa^{-1}}{c^{2}}\left(\tfrac{1}{\sqrt{\Gamma^{k}}}\right)^{2}+1}=:g\left(\tfrac{1}{\sqrt{\Gamma^{k}}}\right)$

Now here is the trick: $g$ is strictly monotonously increasing (its derivative
is positive). So if $n$ was the first number violating our requirement (A.14),
(i.e. $\tfrac{1}{\sqrt{\Gamma^{n}}}<f(n)$), then due to

$\displaystyle f(n-1)\;\mkern-5.0mu plus
-5.0mu\mathop{\leq}\limits^{\text{(\ref{eq: bounding function
properties})}}\;\mkern-5.0mu plus -5.0mu\frac{1}{\sqrt{\Gamma^{n-1}}}\leq
g\left(\tfrac{1}{\sqrt{\Gamma^{n}}}\right)<g(f(n))$

it would be enough for

$\displaystyle g(f(n))\leq f(n-1)$ (A.16)

to be true, to obtain a contradiction. So if $f(0)\leq 1$ is a given, $f(n)$
would then fulfill (A.14) by induction222Unfortunately, I have no idea how
Nesterov came up with this trick, nor how he eventually found the function $f$
that satisfies the requirements. As we know the result we can define

$\displaystyle
f(n):=\tfrac{c}{2\sqrt{\kappa^{-1}}}\left[\exp\left(\tfrac{n+1}{2\sqrt{\kappa}}\right)-\exp\left(-\tfrac{n+1}{2\sqrt{\kappa}}\right)\right]=\frac{c}{4\delta}\left[e^{(n+1)\delta}-e^{-(n+1)\delta}\right]$

for $\delta:=\sqrt{\kappa^{-1}}/2\leq 1/2$. Since $f(0)$ is increasing in
$\delta$ we have the induction start

$\displaystyle
f(0)=\tfrac{c}{4\delta}\left[e^{\delta}-e^{-\delta}\right]\leq\tfrac{c}{2}\left[e^{1/2}-e^{-1/2}\right]\leq\tfrac{c}{\sqrt{3}}\leq
1.$

Now we only need to show (A.16) to be true to finish the proof:

$\displaystyle g(f(n))$
$\displaystyle=f(n)-\frac{c}{2}\sqrt{\tfrac{\kappa^{-1}}{c^{2}}f(n)^{2}+1}$
$\displaystyle=f(n)-\frac{c}{2}\sqrt{\tfrac{1}{2^{2}}\left[e^{(n+1)\delta}-e^{-(n+1)\delta}\right]^{2}+1}$
$\displaystyle=f(n)-\frac{c}{2}\sqrt{\tfrac{1}{2}\left[\tfrac{1}{2}e^{2(n+1)\delta}-\smash{\underbrace{e^{(n+1)\delta}e^{-(n+1)\delta}}_{=1}}+\tfrac{1}{2}e^{-2(n+1)\delta}\right]+1}$
$\displaystyle=f(n)-\underbrace{\frac{c}{4}\left[e^{(n+1)\delta}+e^{-(n+1)\delta}\right]}_{=f^{\prime}(n)}$
$\displaystyle=f(n)+\langle f^{\prime}(n),n-(n-1)\rangle\;\mkern-5.0mu plus
-5.0mu\mathop{\leq}\limits^{f\text{ convex}}\;\mkern-5.0mu plus
-5.0muf(n-1)\qed$

######  Lemma A.2.6.

$\displaystyle
c(\gamma_{1})=2L\frac{1-\gamma_{1}(1+\kappa^{-1})+\gamma_{1}^{2}}{\gamma_{1}^{2}-\kappa^{-1}}\qquad\gamma_{1}\in(\sqrt{\kappa^{-1}},1).$

is monotonously decreasing.

###### Proof.

We discard $2L$ without loss of generality to get

$\displaystyle\frac{dc}{d\gamma_{1}}$
$\displaystyle=\frac{(-(1+\kappa^{-1})+2\gamma_{1})(\gamma_{1}^{2}-\kappa^{-1})-(1-\gamma_{1}(1+\kappa^{-1})+\gamma_{1}^{2})2\gamma_{1}}{(\gamma_{1}^{2}-\kappa^{-1})^{2}}.$

Since we are interested in the sign of the derivative we can discard the
denominator and look at just the enumerator

$\displaystyle
2\gamma_{1}(\cancel{\gamma_{1}^{2}}-\kappa^{-1})-(1+\kappa^{-1})(\gamma_{1}^{2}-\kappa^{-1})-(1-\gamma_{1}(1+\kappa^{-1})+\cancel{\gamma_{1}^{2}})2\gamma_{1}$
$\displaystyle=-2\gamma_{1}\kappa^{-1}-(1+\kappa^{-1})(\gamma_{1}^{2}+\kappa^{-1})-2\gamma_{1}+2\gamma_{1}^{2}(1+\kappa^{-1})$
$\displaystyle=(1+\kappa^{-1})\underbrace{(\kappa^{-1}-2\gamma_{1}+2\gamma_{1}^{2})}_{=:f(\gamma_{1})}.$

In which we only need to consider the parabola $f$. Now since the parabola is
convex, it is negative for all points between two ends where it is negative.
And we have

$\displaystyle f(\sqrt{\kappa^{-1}})$
$\displaystyle=2(\kappa^{-1}-\sqrt{\kappa^{-1}})<0$ $\displaystyle f(1)$
$\displaystyle=\kappa^{-1}-1<0\qed$

##  Appendix B Code for Visualizations

The Programming Language Julia [4] was used to generate all the visualizations
without citation in this work. The code can be found in various .jl files on
GitHub at https://github.com/FelixBenning/masterthesis/tree/main/media. While
these are normal Julia files which can be run as-is, they were written for the
Pluto.jl notebook, which is required for the interactive features of these
plots.

## References

  * [1] Stefan Ankirchner and Stefan Perko “Approximating Stochastic Gradient Descent with Diffusions: Error Expansions and Impact of Learning Rate Schedules”, 2021 URL: https://hal.archives-ouvertes.fr/hal-03262396
  * [2] Larry Armijo “Minimization of Functions Having Lipschitz Continuous First Partial Derivatives.” In _Pacific Journal of Mathematics_ 16.1 Pacific Journal of Mathematics, A Non-profit Corporation, 1966, pp. 1–3 URL: https://projecteuclid.org/journals/pacific-journal-of-mathematics/volume-16/issue-1/Minimization-of-functions-having-Lipschitz-continuous-first-partial-derivatives/pjm/1102995080.full
  * [3] Francis Bach and Eric Moulines “Non-Strongly-Convex Smooth Stochastic Approximation with Convergence Rate O(1/n)”, 2013 arXiv: http://arxiv.org/abs/1306.2119
  * [4] Jeff Bezanson, Alan Edelman, Stefan Karpinski and Viral B. Shah “Julia: A Fresh Approach to Numerical Computing” In _SIAM Review_ 59.1 Society for Industrial and Applied Mathematics, 2017, pp. 65–98 DOI: 10.1137/141000671
  * [5] Léon Bottou, Frank E. Curtis and Jorge Nocedal “Optimization Methods for Large-Scale Machine Learning”, 2018 arXiv: http://arxiv.org/abs/1606.04838
  * [6] Clément Bouttier and Ioana Gavra “Convergence Rate of a Simulated Annealing Algorithm with Noisy Observations” In _Journal of Machine Learning Research_ 20.4, 2019, pp. 1–45 URL: http://jmlr.org/papers/v20/16-588.html
  * [7] Anton Bovier and Frank Hollander “Metastability: A Potential-Theoretic Approach”, Grundlehren Der Mathematischen Wissenschaften Springer International Publishing, 2015 DOI: 10.1007/978-3-319-24777-9
  * [8] Alan J. Bray and David S. Dean “The Statistics of Critical Points of Gaussian Fields on Large-Dimensional Spaces” In _Physical Review Letters_ 98.15, 2007 DOI: 10.1103/PhysRevLett.98.150201
  * [9] C.. Broyden “The Convergence of a Class of Double-rank Minimization Algorithms 1. General Considerations” In _IMA Journal of Applied Mathematics_ 6.1, 1970, pp. 76–90 DOI: 10.1093/imamat/6.1.76
  * [10] Sébastien Bubeck “Convex Optimization: Algorithms and Complexity”, 2015 arXiv: http://arxiv.org/abs/1405.4980
  * [11] Sébastien Bubeck, Yin Tat Lee and Mohit Singh “A Geometric Alternative to Nesterov’s Accelerated Gradient Descent”, 2015 arXiv: http://arxiv.org/abs/1506.08187
  * [12] Augustin Cauchy “Méthode Générale Pour La Résolution Des Systemes d’équations Simultanées” In _Comptes rendus de l’Académie des Sciences_ 25, 1847, pp. 536–538
  * [13] Yuxin Chen “Large-Scale Optimization for Data Science: Mirror Descent”, 2019 URL: http://www.princeton.edu/~yc5/ele522_optimization/lectures/mirror_descent.pdf
  * [14] Anna Choromanska et al. “The Loss Surfaces of Multilayer Networks”, 2015 arXiv: http://arxiv.org/abs/1412.0233
  * [15] Christian Darken and John Moody “Towards Faster Stochastic Gradient Search” In _NIPs_ 91, 1991, pp. 1009–1016
  * [16] Yann Dauphin et al. “Identifying and Attacking the Saddle Point Problem in High-Dimensional Non-Convex Optimization”, 2014 arXiv:1406.2572
  * [17] Olivier Devolder, François Glineur and Yurii Nesterov “First-Order Methods of Smooth Convex Optimization with Inexact Oracle” In _Mathematical Programming_ 146.1, 2014, pp. 37–75 DOI: 10.1007/s10107-013-0677-5
  * [18] (https://stats.stackexchange.com/users/95569/dontloo) dontloo “What’s the Difference between Momentum Based Gradient Descent and Nesterov’s Accelerated Gradient Descent?”, Cross Validated, 2016 Cross Validated Stack Exchange URL: https://stats.stackexchange.com/q/191727
  * [19] Timothy Dozat “Incorporating Nesterov Momentum into Adam”, 2016 URL: https://openreview.net/forum?id=OM0jvwB8jIp57ZJjtNEZ
  * [20] John Duchi, Elad Hazan and Yoram Singer “Adaptive Subgradient Methods for Online Learning and Stochastic Optimization” In _The Journal of Machine Learning Research_ 12, 2011, pp. 2121–2159
  * [21] Nicolas Flammarion and Francis Bach “From Averaging to Acceleration, There Is Only a Step-size”, 2015 arXiv: http://arxiv.org/abs/1504.01577
  * [22] R. Fletcher “A New Approach to Variable Metric Algorithms” In _The Computer Journal_ 13.3, 1970, pp. 317–322 DOI: 10.1093/comjnl/13.3.317
  * [23] Timur Garipov et al. “Loss Surfaces, Mode Connectivity, and Fast Ensembling of DNNs”, 2018 arXiv: http://arxiv.org/abs/1802.10026
  * [24] Israel Gelfand “Normierte Ringe” In _Matematicheskii Sbornik_ 9(51).1, 1941, pp. 3–24 URL: http://www.mathnet.ru/links/78bac43c7b4d5544986886ed8be19701/sm6046.pdf
  * [25] Euhanna Ghadimi, Hamid Reza Feyzmahdavian and Mikael Johansson “Global Convergence of the Heavy-ball Method for Convex Optimization” In _2015 European Control Conference (ECC)_ , 2015, pp. 310–315 DOI: 10.1109/ECC.2015.7330562
  * [26] Gabriel Goh “Why Momentum Really Works” In _Distill_ 2.4, 2017, pp. e6 DOI: 10.23915/distill.00006
  * [27] Donald Goldfarb “A Family of Variable-Metric Methods Derived by Variational Means” In _Mathematics of Computation_ 24.109, 1970, pp. 23–26 DOI: 10.1090/S0025-5718-1970-0258249-6
  * [28] Gene H. Golub and Charles F. Van Loan “Matrix Computations”, Johns Hopkins Studies in the Mathematical Sciences Baltimore: The Johns Hopkins University Press, 2013 URL: http://math.ecnu.edu.cn/~jypan/Teaching/books/2013%20Matrix%20Computations%204th.pdf
  * [29] Anupam Gupta “Advanced Algorithms, Fall 2020, Lecture 19”, 2020 URL: http://www.cs.cmu.edu/~15850/notes/lec19.pdf
  * [30] Moritz Hardt, Benjamin Recht and Yoram Singer “Train Faster, Generalize Better: Stability of Stochastic Gradient Descent”, 2016 arXiv: http://arxiv.org/abs/1509.01240
  * [31] Jun He and Xinghuo Yu “Conditions for the Convergence of Evolutionary Algorithms” In _Journal of Systems Architecture_ 47.7, Evolutionary Computing, 2001, pp. 601–612 DOI: 10.1016/S1383-7621(01)00018-2
  * [32] Geoffrey Hinton “Neural Networks for Machine Learning”, 2012 URL: https://www.cs.toronto.edu/~hinton/coursera_lectures.html
  * [33] Sepp Hochreiter and Jürgen Schmidhuber “Flat Minima” In _Neural Computation_ 9.1, 1997, pp. 1–42 DOI: 10.1162/neco.1997.9.1.1
  * [34] “IMOmath: L’Hopital’s Theorem” URL: https://www.imomath.com/index.php?options=686
  * [35] Robert (https://math.stackexchange.com/users/8508) Israel “How to Work out This Problem Using Stolz’s Theorem?”, 2012 Mathematics Stack Exchange URL: https://math.stackexchange.com/q/207596
  * [36] Pavel Izmailov et al. “Averaging Weights Leads to Wider Optima and Better Generalization”, 2019 arXiv: http://arxiv.org/abs/1803.05407
  * [37] Hamed Karimi, Julie Nutini and Mark Schmidt “Linear Convergence of Gradient and Proximal-Gradient Methods Under the Polyak-Łojasiewicz Condition”, 2020 arXiv: http://arxiv.org/abs/1608.04636
  * [38] Diederik P. Kingma and Jimmy Ba “Adam: A Method for Stochastic Optimization”, 2017 arXiv: http://arxiv.org/abs/1412.6980
  * [39] Victor Kozyakin “On Accuracy of Approximation of the Spectral Radius by the Gelfand Formula” In _Linear Algebra and its Applications_ 431.11, 2009, pp. 2134–2141 DOI: 10.1016/j.laa.2009.07.008
  * [40] Harold J. Kushner “Stochastic Approximation Algorithms and Applications”, Applications of Mathematics; 35 Berlin, et al.: Springer, 1997
  * [41] Laurent Lessard, Benjamin Recht and Andrew Packard “Analysis and Design of Optimization Algorithms via Integral Quadratic Constraints” In _SIAM Journal on Optimization_ 26.1 Society for Industrial and Applied Mathematics, 2016, pp. 57–95 DOI: 10.1137/15M1009597
  * [42] Kenneth Levenberg “A Method for the Solution of Certain Non-Linear Problems in Least Squares” In _Quarterly of Applied Mathematics_ 2.2, 1944, pp. 164–168 DOI: 10.1090/qam/10666
  * [43] Hao Li et al. “Visualizing the Loss Landscape of Neural Nets”, 2018 arXiv: http://arxiv.org/abs/1712.09913
  * [44] Qianxiao Li, Cheng Tai and Weinan E “Stochastic Modified Equations and Adaptive Stochastic Gradient Algorithms” In _Proceedings of the 34th International Conference on Machine Learning_ PMLR, 2017, pp. 2101–2110 URL: https://proceedings.mlr.press/v70/li17f.html
  * [45] Donald W. Marquardt “An Algorithm for Least-Squares Estimation of Nonlinear Parameters” In _Journal of the Society for Industrial and Applied Mathematics_ 11.2 Society for Industrial and Applied Mathematics, 1963, pp. 431–441 DOI: 10.1137/0111030
  * [46] James Martens “Deep Learning via Hessian-free Optimization” In _ICML_ 27, 2010, pp. 8
  * [47] Luke Metz, C. Freeman, Samuel S. Schoenholz and Tal Kachman “Gradients Are Not All You Need”, 2021 arXiv: http://arxiv.org/abs/2111.05803
  * [48] A. Nemirovski, A. Juditsky, G. Lan and A. Shapiro “Robust Stochastic Approximation Approach to Stochastic Programming” In _SIAM Journal on Optimization_ 19.4 Society for Industrial and Applied Mathematics, 2009, pp. 1574–1609 DOI: 10.1137/070704277
  * [49] Yurii Evgen’evič Nesterov “A Method for Solving the Convex Programming Problem with Convergence Rate $O(1/k^{2})$” In _Dokl. Akad. Nauk SSSR_ 269, 1983, pp. 543–547 URL: http://m.mathnet.ru/php/archive.phtml?wshow=paper&jrnid=dan&paperid=46009&option_lang=eng
  * [50] Yurii Evgen’evič Nesterov “Lectures on Convex Optimization”, Springer Optimization and Its Applications; Volume 137 Cham: Springer, 2018 DOI: 10.1007/978-3-319-91578-4
  * [51] Thanh Huy Nguyen, Umut Şimşekli, Mert Gürbüzbalaban and Gaël Richard “First Exit Time Analysis of Stochastic Gradient Descent Under Heavy-Tailed Gradient Noise”, 2019 arXiv: http://arxiv.org/abs/1906.09069
  * [52] Jorge Nocedal and Stephen J. Wright “Quasi-Newton Methods” In _Numerical Optimization_ , Springer Series in Operations Research and Financial Engineering New York, NY: Springer, 2006, pp. 135–163 DOI: 10.1007/978-0-387-40065-5˙6
  * [53] Razvan Pascanu, Yann N. Dauphin, Surya Ganguli and Yoshua Bengio “On the Saddle Point Problem for Non-Convex Optimization”, 2014 arXiv: http://arxiv.org/abs/1405.4604
  * [54] Barak A. Pearlmutter “Fast Exact Multiplication by the Hessian” In _Neural Computation_ 6.1, 1994, pp. 147–160 DOI: 10.1162/neco.1994.6.1.147
  * [55] Scott Pesme, Aymeric Dieuleveut and Nicolas Flammarion “On Convergence-Diagnostic Based Step Sizes for Stochastic Gradient Descent” In _Proceedings of the 37th International Conference on Machine Learning_ PMLR, 2020, pp. 7641–7651 URL: https://proceedings.mlr.press/v119/pesme20a.html
  * [56] Boris T. Polyak “Some Methods of Speeding up the Convergence of Iteration Methods” In _USSR Computational Mathematics and Mathematical Physics_ 4.5, 1964, pp. 1–17 DOI: 10.1016/0041-5553(64)90137-5
  * [57] Boris T. Polyak “Introduction to Optimization” New York Optimization Software, 1987
  * [58] M… Powell “On Search Directions for Minimization Algorithms” In _Mathematical Programming_ 4.1, 1973, pp. 193–201 DOI: 10.1007/BF01584660
  * [59] Ning Qian “On the Momentum Term in Gradient Descent Learning Algorithms” In _Neural Networks_ 12.1 Pergamon, 1999, pp. 145–151 DOI: 10.1016/S0893-6080(98)00116-6
  * [60] Alec Radford “Visualizing Optimization Algos - GIFs”, 2014 Imgur URL: https://imgur.com/a/Hqolp
  * [61] Benjamin Recht “Optimization I”, 2013 URL: https://simons.berkeley.edu/talks/ben-recht-2013-09-04
  * [62] Sashank J. Reddi, Satyen Kale and Sanjiv Kumar “On the Convergence of Adam and Beyond”, 2019 arXiv: http://arxiv.org/abs/1904.09237
  * [63] Sebastian Ruder “An Overview of Gradient Descent Optimization Algorithms”, 2017 arXiv: http://arxiv.org/abs/1609.04747
  * [64] Robin M. Schmidt, Frank Schneider and Philipp Hennig “Descending through a Crowded Valley – Benchmarking Deep Learning Optimizers”, 2021 arXiv: http://arxiv.org/abs/2007.01547
  * [65] Nicol N. Schraudolph “Fast Curvature Matrix-Vector Products” In _Artificial Neural Networks — ICANN 2001_ , Lecture Notes in Computer Science Berlin, et al.: Springer, 2001, pp. 19–26 DOI: 10.1007/3-540-44668-0˙4
  * [66] Terrence J. Sejnowski “The Unreasonable Effectiveness of Deep Learning in Artificial Intelligence” In _Proceedings of the National Academy of Sciences_ 117.48 National Academy of Sciences, 2020, pp. 30033–30038 DOI: 10.1073/pnas.1907373117
  * [67] D.. Shanno “Conditioning of Quasi-Newton Methods for Function Minimization” In _Mathematics of Computation_ 24.111, 1970, pp. 647–656 DOI: 10.1090/S0025-5718-1970-0274029-X
  * [68] Jonathan Richard Shewchuk “An Introduction to the Conjugate Gradient Method without the Agonizing Pain” Carnegie-Mellon University. Department of Computer Science, 1994 URL: https://www.cs.cmu.edu/~quake-papers/painless-conjugate-gradient.pdf
  * [69] Umut Simsekli, Levent Sagun and Mert Gurbuzbalaban “A Tail-Index Analysis of Stochastic Gradient Noise in Deep Neural Networks”, 2019 arXiv: http://arxiv.org/abs/1901.06053
  * [70] Samuel L. Smith, Pieter-Jan Kindermans, Chris Ying and Quoc V. Le “Don’t Decay the Learning Rate, Increase the Batch Size”, 2018 arXiv: http://arxiv.org/abs/1711.00489
  * [71] Weijie Su, Stephen Boyd and Emmanuel J. Candes “A Differential Equation for Modeling Nesterov’s Accelerated Gradient Method: Theory and Insights”, 2015 arXiv: http://arxiv.org/abs/1503.01243
  * [72] Tao Sun et al. “Non-Ergodic Convergence Analysis of Heavy-Ball Algorithms” In _Proceedings of the AAAI Conference on Artificial Intelligence_ 33.01, 2019, pp. 5033–5040 DOI: 10.1609/aaai.v33i01.33015033
  * [73] Ilya Sutskever, James Martens, George Dahl and Geoffrey Hinton “On the Importance of Initialization and Momentum in Deep Learning” In _International Conference on Machine Learning_ PMLR, 2013, pp. 1139–1147 URL: http://proceedings.mlr.press/v28/sutskever13.html
  * [74] Richard S. Sutton and Andrew G. Barto “Reinforcement Learning: An Introduction”, Adaptive Computation and Machine Learning Series Cambridge, Massachusetts: The MIT Press, 2018
  * [75] Tuyen Trung Truong and Tuan Hang Nguyen “Backtracking Gradient Descent Method for General $C^{1}$ Functions, with Applications to Deep Learning”, 2019 arXiv: http://arxiv.org/abs/1808.05160
  * [76] V.N. Vapnik “An Overview of Statistical Learning Theory” In _IEEE Transactions on Neural Networks_ 10.5, 1999, pp. 988–999 DOI: 10.1109/72.788640
  * [77] Oriol Vinyals and Daniel Povey “Krylov Subspace Descent for Deep Learning” In _Proceedings of the Fifteenth International Conference on Artificial Intelligence and Statistics_ PMLR, 2012, pp. 1261–1268 URL: https://proceedings.mlr.press/v22/vinyals12.html
  * [78] Ashia C. Wilson et al. “The Marginal Value of Adaptive Gradient Methods in Machine Learning”, 2018 arXiv: http://arxiv.org/abs/1705.08292
  * [79] Matthew D. Zeiler “ADADELTA: An Adaptive Learning Rate Method”, 2012 arXiv: http://arxiv.org/abs/1212.5701