https://distill.pub/2017/momentum/

Distill
About Prize Submit
[?]\star[?] +++ --- === a\alphaa l\lambdal b\betab RRR a=\alpha=a= b=\
beta=b= b=0\beta = 0b=0 b=1\beta=1b=1 a=1/li\alpha = 1/\lambda_ia=1/l
 i   model\text{model}model 0p10 p_10p 1   0p-10 \bar{p}_10 p -   1  
2b2\sqrt{\beta}2[?] b    li\lambda_il i   li=0\lambda_i = 0l i  =0 a>1/
li\alpha > 1/\lambda_ia>1/l i   max{|s1|,|s2|}>1\max\{|\sigma_1|,|\
sigma_2|\} > 1max{|s 1  |,|s 2  |}>1 xik-xi*x_i^k - x_i^*x i k  -x i 
*   xi\xi_ix i   b=(1-ali)2\beta = (1 - \sqrt{\alpha \lambda_i})^2b=(
1-[?] al i     ) 2  

Why Momentum Really Works

Step-size a = 0.02
Momentum b = 0.99
We often think of Momentum as a means of dampening oscillations and
speeding up the iterations, leading to faster convergence. But it has
other interesting behavior. It allows a larger range of step-sizes to
be used, and creates its own oscillations. What is going on?
Gabriel Goh UC Davis
April. 4
2017
 
Citation:
Goh, 2017

Here's a popular story about momentum [1, 2, 3]: gradient descent is
a man walking down a hill. He follows the steepest path downwards;
his progress is slow, but steady. Momentum is a heavy ball rolling
down the same hill. The added inertia acts both as a smoother and an
accelerator, dampening oscillations and causing us to barrel through
narrow valleys, small humps and local minima.

This standard story isn't wrong, but it fails to explain many
important behaviors of momentum. In fact, momentum can be understood
far more precisely if we study it on the right model.

One nice model is the convex quadratic. This model is rich enough to
reproduce momentum's local dynamics in real problems, and yet simple
enough to be understood in closed form. This balance gives us
powerful traction for understanding this algorithm.

---------------------------------------------------------------------

We begin with gradient descent. The algorithm has many virtues, but
speed is not one of them. It is simple -- when optimizing a smooth
function fff, we make a small step in the gradient wk+1=wk-a[?]f(wk).w^
{k+1} = w^k-\alpha\nabla f(w^k).w k+1  =w k  -a[?]f(w k  ). For a
step-size small enough, gradient descent makes a monotonic
improvement at every iteration. It always converges, albeit to a
local minimum. And under a few weak curvature conditions it can even
get there at an exponential rate.

But the exponential decrease, though appealing in theory, can often
be infuriatingly small. Things often begin quite well -- with an
impressive, almost immediate decrease in the loss. But as the
iterations progress, things start to slow down. You start to get a
nagging feeling you're not making as much progress as you should be.
What has gone wrong?

The problem could be the optimizer's old nemesis, pathological
curvature. Pathological curvature is, simply put, regions of fff
which aren't scaled properly. The landscapes are often described as
valleys, trenches, canals and ravines. The iterates either jump
between valleys, or approach the optimum in small, timid steps.
Progress along certain directions grind to a halt. In these
unfortunate regions, gradient descent fumbles.

Momentum proposes the following tweak to gradient descent. We give
gradient descent a short-term memory: zk+1=bzk+[?]f(wk)wk+1=wk-azk+1 \
begin{aligned} z^{k+1}&=\beta z^{k}+\nabla f(w^{k})\\[0.4em] w^{k+1}&
=w^{k}-\alpha z^{k+1} \end{aligned}  z k+1   w k+1     =bz k  +[?]f(w k
  ) =w k  -az k+1     The change is innocent, and costs almost
nothing. When b=0\beta = 0b=0 , we recover gradient descent. But for 
b=0.99\beta = 0.99b=0.99 (sometimes 0.9990.9990.999, if things are
really bad), this appears to be the boost we need. Our iterations
regain that speed and boldness it lost, speeding to the optimum with
a renewed energy.

Optimizers call this minor miracle "acceleration".

The new algorithm may seem at first glance like a cheap hack. A
simple trick to get around gradient descent's more aberrant
behavior -- a smoother for oscillations between steep canyons. But the
truth, if anything, is the other way round. It is gradient descent
which is the hack. First, momentum gives up to a quadratic speedup on
many functions. ^1 This is no small matter -- this is similar to the
speedup you get from the Fast Fourier Transform, Quicksort, and
Grover's Algorithm. When the universe gives you quadratic speedups,
you should start to pay attention.

But there's more. A lower bound, courtesy of Nesterov [5], states
that momentum is, in a certain very narrow and technical sense,
optimal. Now, this doesn't mean it is the best algorithm for all
functions in all circumstances. But it does satisfy some curiously
beautiful mathematical properties which scratch a very human itch for
perfection and closure. But more on that later. Let's say this for
now -- momentum is an algorithm for the book.

---------------------------------------------------------------------

First Steps: Gradient Descent

We begin by studying gradient descent on the simplest model possible
which isn't trivial -- the convex quadratic, f(w)=12wTAw-bTw,w[?]Rn. f
(w) = \tfrac{1}{2}w^TAw - b^Tw, \qquad w \in \mathbf{R}^n. f(w)= 2  1
  w T  Aw-b T  w,w[?]R n  . Assume AAA is symmetric and invertible,
then the optimal solution w[?]w^{\star}w [?]   occurs at w[?]=A-1b. w^{\
star} = A^{-1}b.w [?]  =A -1  b. Simple as this model may be, it is
rich enough to approximate many functions (think of AAA as your
favorite model of curvature -- the Hessian, Fisher Information Matrix 
[6], etc) and captures all the key features of pathological
curvature. And more importantly, we can write an exact closed formula
for gradient descent on this function.

This is how it goes. Since [?]f(w)=Aw-b\nabla f(w)=Aw - b[?]f(w)=Aw-b,
the iterates are wk+1=wk-a(Awk-b). w^{k+1}=w^{k}- \alpha (Aw^{k} -
b). w k+1  =w k  -a(Aw k  -b). Here's the trick. There is a very
natural space to view gradient descent where all the dimensions act
independently -- the eigenvectors of AAA.

Every symmetric matrix AAA has an eigenvalue decomposition A=Q diag
(l1,...,ln) QT,Q=[q1,...,qn], A=Q\ \text{diag}(\lambda_{1},\ldots,\
lambda_{n})\ Q^{T},\qquad Q = [q_1,\ldots,q_n], A=Q diag(l 1  ,...,l n 
 ) Q T  ,Q=[q 1  ,...,q n  ], and, as per convention, we will assume
that the li\lambda_il i  's are sorted, from smallest l1\lambda_1l 1 
  to biggest ln\lambda_nl n  . If we perform a change of basis, xk=QT
(wk-w[?])x^{k} = Q^T(w^{k} - w^\star)x k  =Q T  (w k  -w [?]  ), the
iterations break apart, becoming: xik+1=xik-alixik=(1-ali)xik=(1-ali)
k+1xi0 \begin{aligned} x_{i}^{k+1} & =x_{i}^{k}-\alpha \lambda_ix_{i}
^{k} \\[0.4em] &= (1-\alpha\lambda_i)x^k_i=(1-\alpha \lambda_i)^{k+1}
x^0_i \end{aligned}  x i k+1      =x i k  -al i  x i k   =(1-al i  )x
 i k  =(1-al i  ) k+1  x i 0     Moving back to our original space ww
w, we can see that wk-w[?]=Qxk=[?]inxi0(1-ali)kqi w^k - w^\star = Qx^k=\
sum_i^n x^0_i(1-\alpha\lambda_i)^k q_i w k  -w [?]  =Qx k  = i [?] n  x i
 0  (1-al i  ) k  q i   and there we have it -- gradient descent in
closed form.

Decomposing the Error

The above equation admits a simple interpretation. Each element of 
x0x^0x 0   is the component of the error in the initial guess in the 
QQQ-basis. There are nnn such errors, and each of these errors
follows its own, solitary path to the minimum, decreasing
exponentially with a compounding rate of 1-ali1-\alpha\lambda_i1-al i
  . The closer that number is to 111, the slower it converges.

For most step-sizes, the eigenvectors with largest eigenvalues
converge the fastest. This triggers an explosion of progress in the
first few iterations, before things slow down as the smaller
eigenvectors' struggles are revealed. By writing the contributions of
each eigenspace's error to the loss f(wk)-f(w[?])=[?](1-ali)2kli[xi0]2 f
(w^{k})-f(w^{\star})=\sum(1-\alpha\lambda_{i})^{2k}\lambda_{i}[x_{i}^
{0}]^2 f(w k  )-f(w [?]  )=[?](1-al i  ) 2k  l i  [x i 0  ] 2   we can
visualize the contributions of each error component to the loss.

Optimization can be seen as combination of several component
problems, shown here as 1 2 3 with eigenvalues l1=0.01\lambda_1=0.01l
 1  =0.01, l2=0.1\lambda_2=0.1l 2  =0.1, and l3=1\lambda_3=1l 3  =1
respectively.
Step-size
Optimal Step-size

Choosing A Step-size

The above analysis gives us immediate guidance as to how to set a
step-size a\alphaa. In order to converge, each |1-ali||1-\alpha \
lambda_i||1-al i  | must be strictly less than 1. All workable
step-sizes, therefore, fall in the interval 0<ali<2.0<\alpha\lambda_i
<2.0<al i  <2. The overall convergence rate is determined by the
slowest error component, which must be either l1\lambda_1l 1   or ln\
lambda_nl n  : rate(a) = maxi|1-ali| = max{|1-al1|, |1-aln|} \begin
{aligned}\text{rate}(\alpha) & ~=~ \max_{i}\left|1-\alpha\lambda_{i}\
right|\\[0.9em] & ~=~ \max\left\{|1-\alpha\lambda_{1}|,~ |1-\alpha\
lambda_{n}|\right\} \end{aligned}  rate(a)     =  i max  |1-al i  |  
= max{|1-al 1  |, |1-al n  |}  

This overall rate is minimized when the rates for l1\lambda_1l 1  
and ln\lambda_nl n   are the same -- this mirrors our informal
observation in the previous section that the optimal step-size causes
the first and last eigenvectors to converge at the same rate. If we
work this through we get: optimal a = argmina rate(a) =
 2l1+lnoptimal rate = mina rate(a) = ln/l1-1ln/l1+1 \begin{aligned} \
text{optimal }\alpha ~=~{\mathop{\text{argmin}}\limits_\alpha} ~\text
{rate}(\alpha) & ~=~\frac{2}{\lambda_{1}+\lambda_{n}}\\[1.4em] \text
{optimal rate} ~=~{\min_\alpha} ~\text{rate}(\alpha) & ~=~\frac{\
lambda_{n}/\lambda_{1}-1}{\lambda_{n}/\lambda_{1}+1} \end{aligned}  
optimal a =  a argmin   rate(a) optimal rate =  a min   rate(a)    = 
 l 1  +l n    2    =  l n  /l 1  +1  l n  /l 1  -1    

Notice the ratio ln/l1\lambda_n/\lambda_1l n  /l 1   determines the
convergence rate of the problem. In fact, this ratio appears often
enough that we give it a name, and a symbol -- the condition number. 
condition number:=k:=lnl1 \text{condition number} := \kappa :=\frac{\
lambda_n}{\lambda_1} condition number:=k:= l 1    l n     The
condition number means many things. It is a measure of how close to
singular a matrix is. It is a measure of how robust A-1bA^{-1}bA -1  
b is to perturbations in bbb. And, in this context, the condition
number gives us a measure of how poorly gradient descent will
perform. A ratio of k=1\kappa = 1k=1 is ideal, giving convergence in
one step (of course, the function is trivial). Unfortunately the
larger the ratio, the slower gradient descent will be. The condition
number is therefore a direct measure of pathological curvature.

---------------------------------------------------------------------

Example: Polynomial Regression

The above analysis reveals an insight: all errors are not made equal.
Indeed, there are different kinds of errors, nnn to be exact, one for
each of the eigenvectors of AAA. And gradient descent is better at
correcting some kinds of errors than others. But what do the
eigenvectors of AAA mean? Surprisingly, in many applications they
admit a very concrete interpretation.

Lets see how this plays out in polynomial regression. Given 1D data, 
xi\xi_ix i  , our problem is to fit the model model(x)=w1p1(x)+[?]+wnpn
(x)pi=x-xi-1 \text{model}(\xi)=w_{1}p_{1}(\xi)+\cdots+w_{n}p_{n}(\xi)
\qquad p_{i}=\xi\mapsto\xi^{i-1} model(x)=w 1  p 1  (x)+[?]+w n  p n  (
x)p i  =x-x i-1   to our observations, did_id i  . This model, though
nonlinear in the input x\xix, is linear in the weights, and therefore
we can write the model as a linear combination of monomials, like:

Because of the linearity, we can fit this model to our data xi\xi_ix 
i   using linear regression on the model mismatch minimizew12[?]i(model
(xi)-di)2  =  12[?]Zw-d[?]2 \text{minimize}_w \qquad\tfrac{1}{2}\sum_i (\
text{model}(\xi_{i})-d_{i})^{2} ~~=~~ \tfrac{1}{2}\|Zw - d\|^2 
minimize w   2  1   i [?]  (model(x i  )-d i  ) 2    =   2  1  [?]Zw-d[?] 2
   where Z=(1x1x12...x1n-11x2x22...x2n-1[?][?][?][?][?]1xmxm2...xmn-1). Z=\left(\
begin{array}{ccccc} 1 & \xi_{1} & \xi_{1}^{2} & \ldots & \xi_{1}^
{n-1}\\ 1 & \xi_{2} & \xi_{2}^{2} & \ldots & \xi_{2}^{n-1}\\ \vdots &
\vdots & \vdots & \ddots & \vdots\\ 1 & \xi_{m} & \xi_{m}^{2} & \
ldots & \xi_{m}^{n-1} \end{array}\right). Z= [?] [?] [?] [?]   1 1 [?] 1   x 1 
  x 2   [?] x m     x 1 2   x 2 2   [?] x m 2     ... ... [?] ...   x 1 n-1   x 2
 n-1   [?] x m n-1     [?] [?] [?] [?]  .

The path of convergence, as we know, is elucidated when we view the
iterates in the space of QQQ (the eigenvectors of ZTZZ^T ZZ T  Z). So
let's recast our regression problem in the basis of QQQ. First, we do
a change of basis, by rotating www into QwQwQw, and counter-rotating
our feature maps ppp into eigenspace, p-\bar{p} p -  . We can now
conceptualize the same regression as one over a different polynomial
basis, with the model model(x) = x1p-1(x) + [?] + xnp-n(x)p-i=[?]qijpj. \
text{model}(\xi)~=~x_{1}\bar{p}_{1}(\xi)~+~\cdots~+~x_{n}\bar{p}_{n}
(\xi)\qquad \bar{p}_{i}=\sum q_{ij}p_j. model(x) = x 1   p -   1  (x)
 + [?] + x n   p -   n  (x) p -   i  =[?]q ij  p j  . This model is
identical to the old one. But these new features p-\bar{p} p -  
(which I call "eigenfeatures") and weights have the pleasing property
that each coordinate acts independently of the others. Now our
optimization problem breaks down, really, into nnn small 1D
optimization problems. And each coordinate can be optimized greedily
and independently, one at a time in any order, to produce the final,
global, optimum. The eigenfeatures are also much more informative:

The observations in the above diagram can be justified
mathematically. From a statistical point of view, we would like a
model which is, in some sense, robust to noise. Our model cannot
possibly be meaningful if the slightest perturbation to the
observations changes the entire model dramatically. And the
eigenfeatures, the principal components of the data, give us exactly
the decomposition we need to sort the features by its sensitivity to
perturbations in did_id i  's. The most robust components appear in
the front (with the largest eigenvalues), and the most sensitive
components in the back (with the smallest eigenvalues).

This measure of robustness, by a rather convenient coincidence, is
also a measure of how easily an eigenspace converges. And thus, the
"pathological directions" -- the eigenspaces which converge the
slowest -- are also those which are most sensitive to noise! So
starting at a simple initial point like 000 (by a gross abuse of
language, let's think of this as a prior), we track the iterates till
a desired level of complexity is reached. Let's see how this plays
out in gradient descent.

This effect is harnessed with the heuristic of early stopping : by
stopping the optimization early, you can often get better
generalizing results. Indeed, the effect of early stopping is very
similar to that of more conventional methods of regularization, such
as Tikhonov Regression. Both methods try to suppress the components
of the smallest eigenvalues directly, though they employ different
methods of spectral decay.^2 But early stopping has a distinct
advantage. Once the step-size is chosen, there are no regularization
parameters to fiddle with. Indeed, in the course of a single
optimization, we have the entire family of models, from underfitted
to overfitted, at our disposal. This gift, it seems, doesn't come at
a price. A beautiful free lunch [7] indeed.

---------------------------------------------------------------------

The Dynamics of Momentum

Let's turn our attention back to momentum. Recall that the momentum
update is zk+1=bzk+[?]f(wk)wk+1=wk-azk+1. \begin{aligned} z^{k+1}&=\
beta z^{k}+\nabla f(w^{k})\\[0.4em] w^{k+1}&=w^{k}-\alpha z^{k+1}. \
end{aligned}  z k+1   w k+1     =bz k  +[?]f(w k  ) =w k  -az k+1  .  
Since [?]f(wk)=Awk-b\nabla f(w^k) = Aw^k - b[?]f(w k  )=Aw k  -b, the
update on the quadratic is zk+1=bzk+(Awk-b)wk+1=wk-azk+1. \begin
{aligned} z^{k+1}&=\beta z^{k}+ (Aw^{k}-b)\\[0.4em] w^{k+1}&=w^{k}-\
alpha z^{k+1}. \end{aligned}  z k+1   w k+1     =bz k  +(Aw k  -b) =w
 k  -az k+1  .   Following [8], we go through the same motions, with
the change of basis xk=Q(wk-w[?]) x^{k} = Q(w^{k} - w^\star)x k  =Q(w k
  -w [?]  ) and yk=Qzk y^{k} = Qz^{k}y k  =Qz k  , to yield the update
rule yik+1=byik+lixikxik+1=xik-ayik+1. \begin{aligned} y_{i}^{k+1}&=\
beta y_{i}^{k}+\lambda_{i}x_{i}^{k}\\[0.4em] x_{i}^{k+1}&=x_{i}^{k}-\
alpha y_{i}^{k+1}. \end{aligned}  y i k+1   x i k+1     =by i k  +l i
  x i k   =x i k  -ay i k+1  .   in which each component acts
independently of the other components (though xikx^k_ix i k   and 
yiky^k_iy i k   are coupled). This lets us rewrite our iterates as ^3
(yikxik)=Rk(yi0xi0)R=(bli-ab1-ali). \left(\!\!\begin{array}{c} y_{i}^
{k}\\ x_{i}^{k} \end{array}\!\!\right)=R^k\left(\!\!\begin{array}{c}
y_{i}^{0}\\ x_{i}^{0} \end{array}\!\!\right) \qquad R = \left(\!\!\
begin{array}{cc} \beta & \lambda_{i}\\ -\alpha\beta & 1-\alpha\
lambda_{i} \end{array}\!\!\right). ( y i k   x i k    )=R k  ( y i 0 
  x i 0    )R=( b -ab   l i   1-al i    ). There are many ways of
taking a matrix to the kthk^{th}k th   power. But for the 2x22 \times
22x2 case there is an elegant and little known formula [9] in terms
of the eigenvalues of RRR, s1\sigma_1s 1   and s2\sigma_2s 2  . Rk=
{s1kR1-s2kR2s1[?]s2s1k(kR/s1-(k-1)I)s1=s2,Rj=R-sjIs1-s2 \color{#AAA}{\
color{black}{R^{k}}=\begin{cases} \color{black}{\sigma_{1}^{k}}R_{1}-
\color{black}{\sigma_{2}^{k}}R_{2} & \sigma_{1}\neq\sigma_{2}\\ \
sigma_{1}^{k}(kR/\sigma_1-(k-1)I) & \sigma_{1}=\sigma_{2} \end
{cases},\qquad R_{j}=\frac{R-\sigma_{j}I}{\sigma_{1}-\sigma_{2}}} R k
  ={ s 1 k  R 1  -s 2 k  R 2   s 1 k  (kR/s 1  -(k-1)I)   s 1  [?]s 2  
 s 1  =s 2    ,R j  = s 1  -s 2    R-s j  I   This formula is rather
complicated, but the takeaway here is that it plays the exact same
role the individual convergence rates, 1-ali1-\alpha\lambda_i1-al i  
do in gradient descent. But instead of one geometric series, we have
two coupled series, which may have real or complex values. The
convergence rate is therefore the slowest of the two rates, max
{|s1|,|s2|}\max\{|\sigma_{1}|,|\sigma_{2}|\}max{|s 1  |,|s 2  |} ^4.
By plotting this out, we see there are distinct regions of the
parameter space which reveal a rich taxonomy of convergence behavior 
[10]:

Convergence Rate
A plot of max{|s1|,|s2|}\max\{|\sigma_1|, |\sigma_2|\}max{|s 1  |,|s 
2  |} reveals distinct regions, each with its own style of
convergence.

For what values of a\alphaa and b\betab does momentum converge? Since
we need both s1\sigma_1s 1   and s2\sigma_2s 2   to converge, our
convergence criterion is now max{|s1|,|s2|}<1\max\{|\sigma_{1}|,|\
sigma_{2}|\} < 1max{|s 1  |,|s 2  |}<1. The range of available
step-sizes work out ^5 to be 0<ali<2+2bfor0<=b<10<\alpha\lambda_{i}
<2+2\beta \qquad \text{for} \qquad 0 \leq \beta < 10<al i  <2+2bfor0<=
b<1 We recover the previous result for gradient descent when b=0\beta
= 0b=0. But notice an immediate boon we get. Momentum allows us to
crank up the step-size up by a factor of 2 before diverging.

---------------------------------------------------------------------

The Critical Damping Coefficient

The true magic happens, however, when we find the sweet spot of a\
alphaa and b\betab. Let us try to first optimize over b\betab.

Momentum admits an interesting physical interpretation when a\alphaa
is [11] small: it is a discretization of a damped harmonic
oscillator. Consider a physical simulation operating in discrete time
(like a video game).

yik+1y_{i}^{k+1}y i k+1   === +++ lixik\lambda_{i}x_{i}^{k}l i  x i k
   and perturbed by an external force field We can think of -yik-y_i^
k-y i k   as velocity byik\beta y_{i}^{k}by i k   which is dampened
at each step xik+1x_i^{k+1}x i k+1   === xik-ayik+1x_i^k - \alpha y_i
^{k+1}x i k  -ay i k+1   And xxx is our particle's position which is
moved at each step by a small amount in the direction of the velocity
yik+1y^{k+1}_iy i k+1  .

We can break this equation apart to see how each component affects
the dynamics of the system. Here we plot, for 150150150 iterates, the
particle's velocity (the horizontal axis) against its position (the
vertical axis), in a phase diagram.

This system is best imagined as a weight suspended on a spring. We
pull the weight down by one unit, and we study the path it follows as
it returns to equilibrium. In the analogy, the spring is the source
of our external force lixik\lambda_ix^k_il i  x i k  , and
equilibrium is the state when both the position xikx^k_ix i k   and
the speed yiky^k_iy i k   are 0. The choice of b\betab crucially
affects the rate of return to equilibrium.

The critical value of b=(1-ali)2\beta = (1 - \sqrt{\alpha \lambda_i})
^2b=(1-[?] al i     ) 2   gives us a convergence rate (in eigenspace ii
i) of 1-ali.1 - \sqrt{\alpha\lambda_i}.1-[?] al i     . A square root
improvement over gradient descent, 1-ali1-\alpha\lambda_i1-al i  !
Alas, this only applies to the error in the ithi^{th}i th  
eigenspace, with a\alphaa fixed.

Optimal parameters

To get a global convergence rate, we must optimize over both a\alphaa
and b\betab. This is a more complicated affair,^6 but they work out
to be a=(2l1+ln)2b=(ln-l1ln+l1)2 \alpha = \left(\frac{2}{\sqrt{\
lambda_{1}}+\sqrt{\lambda_{n}}}\right)^{2} \quad \beta = \left(\frac
{\sqrt{\lambda_{n}}-\sqrt{\lambda_{1}}}{\sqrt{\lambda_{n}}+\sqrt{\
lambda_{1}}}\right)^{2} a=( [?] l 1     +[?] l n       2  ) 2  b=( [?] l n 
    +[?] l 1       [?] l n     -[?] l 1       ) 2   Plug this into the
convergence rate, and you get

k-1k+1\frac{\sqrt{\kappa}-1}{\sqrt{\kappa}+1} [?] k   +1  [?] k   -1   
Convergence rate, Momentum k-1k+1 \frac{\kappa-1}{\kappa+1} k+1  k-1 
  Convergence rate, Gradient Descent

With barely a modicum of extra effort, we have essentially square
rooted the condition number! These gains, in principle, require
explicit knowledge of l1\lambda_1l 1   and ln\lambda_nl n  . But the
formulas reveal a simple guideline. When the problem's conditioning
is poor, the optimal a\alphaa is approximately twice that of gradient
descent, and the momentum term is close to 111. So set b\betab as
close to 111 as you can, and then find the highest a\alphaa which
still converges. Being at the knife's edge of divergence, like in
gradient descent, is a good place to be.

We can do the same decomposition here with momentum, with eigenvalues
l1=0.01\lambda_1=0.01l 1  =0.01, l2=0.1\lambda_2=0.1l 2  =0.1, and l3
=1\lambda_3=1l 3  =1. Though the decrease is no longer monotonic, but
significantly faster. f(wk)-f(w[?])f(w^k) - f(w^\star)f(w k  )-f(w [?]  )
Note that the optimal parameters do not necessarily imply the fastest
convergence, though, only the fastest asymptotic convergence rate.
Step-size a =
Momentum b =

While the loss function of gradient descent had a graceful, monotonic
curve, optimization with momentum displays clear oscillations. These
ripples are not restricted to quadratics, and occur in all kinds of
functions in practice. They are not cause for alarm, but are an
indication that extra tuning of the hyperparameters is required.

---------------------------------------------------------------------

Example: The Colorization Problem

Let's look at how momentum accelerates convergence with a concrete
example. On a grid of pixels let GGG be the graph with vertices as
pixels, EEE be the set of edges connecting each pixel to its four
neighboring pixels, and DDD be a small set of a few distinguished
vertices. Consider the problem of minimizing

minimize\text{minimize} minimize 12[?]i[?]D(wi-1)2\qquad \frac{1}{2} \
sum_{i\in D} (w_i - 1)^2  2  1   i[?]D [?]  (w i  -1) 2   The colorizer
pulls distinguished pixels towards 1 +++ 12[?]i,j[?]E(wi-wj)2.\frac{1}{2}
\sum_{i,j\in E} (w_i - w_j)^2. 2  1   i,j[?]E [?]  (w i  -w j  ) 2  . The
smoother spreads out the color

The optimal solution to this problem is a vector of all 111's ^7. An
inspection of the gradient iteration reveals why we take a long time
to get there. The gradient step, for each component, is some form of
weighted average of the current value and its neighbors: wik+1=
wik-a[?]j[?]N(wik-wjk)-{a(wik-1)i[?]D0i[?]D w_{i}^{k+1}=w_{i}^{k}-\alpha\sum_
{j\in N}(w_{i}^{k}-w_{j}^{k})-\begin{cases} \alpha(w_{i}^{k}-1) & i\
in D\\ 0 & i\notin D \end{cases} w i k+1  =w i k  -a j[?]N [?]  (w i k  -
w j k  )-{ a(w i k  -1) 0   i[?]D i[?]D   This kind of local averaging is
effective at smoothing out local variations in the pixels, but poor
at taking advantage of global structure. The updates are akin to a
drop of ink, diffusing through water. Movement towards equilibrium is
made only through local corrections and so, left undisturbed, its
march towards the solution is slow and laborious. Fortunately,
momentum speeds things up significantly.

The eigenvectors of the colorization problem form a generalized
Fourier basis for RnR^nR n  . The smallest eigenvalues have low
frequencies, hence gradient descent corrects high frequency errors
well but not low frequency ones.

In vectorized form, the colorization problem is

minimize\text{minimize}minimize The smoother's quadratic form is the
Graph Laplacian 12[?]i[?]D(xTeieiTx-eiTx)\frac{1}{2}\sum_{i\in D}\left(x^
{T}e_{i}e_{i}^{T}x-e_{i}^{T}x\right) 2  1   i[?]D [?]  (x T  e i  e i T  
x-e i T  x) +++ 12xTLGx\frac{1}{2}x^{T}L_{G}x 2  1  x T  L G  x And
the colorizer is a small low rank correction with a linear term. 
eie_ie i   is the ithi^{th}i th   unit vector.

The Laplacian matrix, LGL_GL G   ^8, which dominates the behavior of
the optimization problem, is a valuable bridge between linear algebra
and graph theory. This is a rich field of study, but one fact is
pertinent to our discussion here. The conditioning of LGL_GL G  ,
here defined as the ratio of the second eigenvector to the last (the
first eigenvalue is always 0, with eigenvector equal to the matrix of
all 1's), is directly connected to the connectivity of the graph.

Small world graphs, like expanders and dense graphs, have excellent
conditioning The conditioning of grids improves with its
dimensionality. And long, wiry graphs, like paths, condition poorly.

These observations carry through to the colorization problem, and the
intuition behind it should be clear. Well connected graphs allow
rapid diffusion of information through the edges, while graphs with
poor connectivity do not. And this principle, taken to the extreme,
furnishes a class of functions so hard to optimize they reveal the
limits of first order optimization.

---------------------------------------------------------------------

The Limits of Descent

Let's take a step back. We have, with a clever trick, improved the
convergence of gradient descent by a quadratic factor with the
introduction of a single auxiliary sequence. But is this the best we
can do? Could we improve convergence even more with two sequences?
Could one perhaps choose the a\alphaa's and b\betab's intelligently
and adaptively? It is tempting to ride this wave of optimism - to the
cube root and beyond!

Unfortunately, while improvements to the momentum algorithm do exist,
they all run into a certain, critical, almost inescapable lower
bound.

Adventures in Algorithmic Space

To understand the limits of what we can do, we must first formally
define the algorithmic space in which we are searching. Here's one
possible definition. The observation we will make is that both
gradient descent and momentum can be "unrolled". Indeed, since w1=
w0 - a[?]f(w0)w2=w1 - a[?]f(w1)=w0 - a[?]f(w0) - a[?]f(w1) [?]wk+1=w0 - a[?]f(w0)
 -    [?][?]    - a[?]f(wk) \begin{array}{lll} w^{1} & \!= & \!w^{0} ~-~ \
alpha\nabla f(w^{0})\\[0.35em] w^{2} & \!= & \!w^{1} ~-~ \alpha\nabla
f(w^{1})\\[0.35em] & \!= & \!w^{0} ~-~ \alpha\nabla f(w^{0}) ~-~ \
alpha\nabla f(w^{1})\\[0.35em] & ~ \!\vdots \\ w^{k+1} & \!= & \!w^
{0} ~-~ \alpha\nabla f(w^{0}) ~-~~~~ \cdots\cdots ~~~~-~ \alpha\nabla
f(w^{k}) \end{array}  w 1   w 2     w k+1     = = =  [?] =   w 0   - a[?]
f(w 0  ) w 1   - a[?]f(w 1  ) w 0   - a[?]f(w 0  ) - a[?]f(w 1  ) w 0   - a
[?]f(w 0  ) -    [?][?]    - a[?]f(w k  )   we can write gradient descent as 
wk+1  =  w0 - a[?]ik[?]f(wi). w^{k+1} ~~=~~ w^{0} ~-~ \alpha\sum_i^k\
nabla f(w^{i}). w k+1    =  w 0   - a i [?] k  [?]f(w i  ). A similar
trick can be done with momentum: wk+1  =  w0 + a[?]ik(1-bk+1-i)1-b[?]f
(wi). w^{k+1} ~~=~~ w^{0} ~+~ \alpha\sum_i^k\frac{(1-\beta^{k+1-i})}
{1-\beta}\nabla f(w^i). w k+1    =  w 0   + a i [?] k   1-b  (1-b k+1-i
  )  [?]f(w i  ). In fact, all manner of first order algorithms,
including the Conjugate Gradient algorithm, AdaMax, Averaged Gradient
and more, can be written (though not quite so neatly) in this
unrolled form. Therefore the class of algorithms for which wk+1  =
  w0 + [?]ikgik[?]f(wi) for some gik w^{k+1} ~~=~~ w^{0} ~+~ \sum_{i}^{k}
\gamma_{i}^{k}\nabla f(w^{i}) \qquad \text{ for some } \gamma_{i}^{k}
w k+1    =  w 0   +  i [?] k  g i k  [?]f(w i  ) for some g i k  
contains momentum, gradient descent and a whole bunch of other
algorithms you might dream up. This is what is assumed in Assumption
2.1.4 [5] of Nesterov. But let's push this even further, and expand
this class to allow different step-sizes for different directions. 
wk+1  =  w0 + [?]ikGik[?]f(wi) for some diagonal matrix Gik. w^{k+1} ~~=
~~ w^{0} ~+~ \sum_{i}^{k}\Gamma_{i}^{k}\nabla f(w^{i}) \quad \text{
for some diagonal matrix } \Gamma_{i}^{k} . w k+1    =  w 0   +  i [?] 
k  G i k  [?]f(w i  ) for some diagonal matrix G i k  . This class of
methods covers most of the popular algorithms for training neural
networks, including ADAM and AdaGrad. We shall refer to this class of
methods as "Linear First Order Methods", and we will show a single
function all these methods ultimately fail on.

The Resisting Oracle

Earlier, when we talked about the colorizer problem, we observed that
wiry graphs cause bad conditioning in our optimization problem.
Taking this to its extreme, we can look at a graph consisting of a
single path -- a function so badly conditioned that Nesterov called a
variant of it "the worst function in the world". The function follows
the same structure as the colorizer problem, and we shall call this
the Convex Rosenbrock,

fn(w)f^n(w)f n  (w) === with a colorizer of one node 12(w1-1)2\frac
{1}{2}\left(w_{1}-1\right)^{2} 2  1  (w 1  -1) 2   +++ 12[?]i=1n
(wi-wi+1)2\frac{1}{2}\sum_{i=1}^{n}(w_{i}-w_{i+1})^{2} 2  1   i=1 [?] n
  (w i  -w i+1  ) 2   strong couplings of adjacent nodes in the path,
+++ 2k-1[?]w[?]2.\frac{2}{\kappa-1}\|w\|^{2}. k-1  2  [?]w[?] 2  . and a
small regularization term.

The optimal solution of this problem is wi[?]=(k-1k+1)i w_{i}^{\star}=\
left(\frac{\sqrt{\kappa}-1}{\sqrt{\kappa}+1}\right)^{i} w i [?]  =( [?] k
   +1  [?] k   -1  ) i   and the condition number of the problem fnf^nf
 n   approaches k\kappak as nnn goes to infinity. Now observe the
behavior of the momentum algorithm on this function, starting from w0
=0w^0 = 0w 0  =0.

Step-size a =
Momentum b =
Here we see the first 50 iterates of momentum on the Convex
Rosenbrock for n=25n=25n=25. The behavior here is similar to that of
any Linear First Order Algorithm.
This triangle is a "dead zone" of our iterates. The iterates are
always 0, no matter what the parameters. The remaining expanding
space is the "light cone" of our iterate's influence. Momentum does
very well here with the optimal parameters.
Error
Weights

The observations made in the above diagram are true for any Linear
First Order algorithm. Let us prove this. First observe that each
component of the gradient depends only on the values directly before
and after it: [?]f(x)i=2wi-wi-1-wi+1+4k-1wi,i[?]1. \nabla f(x)_{i}=2w_{i}
-w_{i-1}-w_{i+1} +\frac{4}{\kappa-1} w_{i}, \qquad i \neq 1. [?]f(x) i 
 =2w i  -w i-1  -w i+1  + k-1  4  w i  ,i[?]1. Therefore the fact we
start at 0 guarantees that that component must remain stoically there
till an element either before or after it turns nonzero. And
therefore, by induction, for any linear first order algorithm,

w0=[  0,0,0,...0,0,...0 ]w1=[ w11,0,0,...0,0,...0 ]w2=[ w12,w22,0,...0,0,...0 ]
 [?]wk=[ w1k,w2k,w3k,...wkk,0,...0 ]. \begin{array}{lllllllll} w^{0} & = &
[~~0, & 0, & 0, & \ldots & 0, & 0, & \ldots & 0~]\\[0.35em] w^{1} & =
& [~w_{1}^{1}, & 0, & 0, & \ldots & 0, & 0, & \ldots & 0~]\\[0.35em]
w^{2} & = & [~w_{1}^{2}, & w_{2}^{2}, & 0, & \ldots & 0, & 0, & \
ldots & 0~]\\[0.35em] & ~ \vdots \\ w^{k} & = & [~w_{1}^{k}, & w_{2}^
{k}, & w_{3}^{k}, & \ldots & w_{k}^{k}, & 0, & \ldots & 0~].\\ \end
{array}  w 0   w 1   w 2    w k      = = =  [?] =   [  0, [ w 1 1  , [ 
w 1 2  , [ w 1 k  ,   0, 0, w 2 2  , w 2 k  ,   0, 0, 0, w 3 k  ,   ...
 ... ... ...   0, 0, 0, w k k  ,   0, 0, 0, 0,   ... ... ... ...   0 ] 0 ] 0 ] 0 ].
  

Think of this restriction as a "speed of light" of information
transfer. Error signals will take at least kkk steps to move from 
w0w_0w 0   to wkw_kw k  . We can therefore sum up the errors which
cannot have changed yet^9: [?]wk-w[?][?][?]>=maxi>=k+1{|wi[?]|}=(k-1k+1)k+1=
(k-1k+1)k[?]w0-w[?][?][?]. \begin{aligned} \|w^{k}-w^{\star}\|_{\infty}&\geq\
max_{i\geq k+1}\{|w_{i}^{\star}|\}\\[0.9em]&=\left(\frac{\sqrt{\
kappa}-1}{\sqrt{\kappa}+1}\right)^{k+1}\\[0.9em]&=\left(\frac{\sqrt{\
kappa}-1}{\sqrt{\kappa}+1}\right)^{k}\|w^{0}-w^{\star}\|_{\infty}. \
end{aligned}  [?]w k  -w [?]  [?] [?]       >= i>=k+1 max  {|w i [?]  |} =( [?] k  
 +1  [?] k   -1  ) k+1   =( [?] k   +1  [?] k   -1  ) k  [?]w 0  -w [?]  [?] [?]  .
   As nnn gets large, the condition number of fnf^nf n   approaches k
\kappak. And the gap therefore closes; the convergence rate that
momentum promises matches the best any linear first order algorithm
can do. And we arrive at the disappointing conclusion that on this
problem, we cannot do better.

Like many such lower bounds, this result must not be taken literally,
but spiritually. It, perhaps, gives a sense of closure and finality
to our investigation. But this is not the final word on first order
optimization. This lower bound does not preclude the possibility, for
example, of reformulating the problem to change the condition number
itself! There is still much room for speedups, if you understand the
right places to look.

Momentum with Stochastic Gradients

There is a final point worth addressing. All the discussion above
assumes access to the true gradient -- a luxury seldom afforded in
modern machine learning. Computing the exact gradient requires a full
pass over all the data, the cost of which can be prohibitively
expensive. Instead, randomized approximations of the gradient, like
minibatch sampling, are often used as a plug-in replacement of [?]f(w)\
nabla f(w)[?]f(w). We can write the approximation in two parts,

[?]f(w)\nabla f(w)[?]f(w) the true gradient +++ error(w).\text{error}(w).
error(w). and an approximation error.
If the estimator is unbiased e.g. E[error(w)]=0\mathbf{E}[\text
{error}(w)] = 0E[error(w)]=0

It is helpful to think of our approximate gradient as the injection
of a special kind of noise into our iteration. And using the
machinery developed in the previous sections, we can deal with this
extra term directly. On a quadratic, the error term cleaves cleanly
into a separate term, where ^10

(yikxik) \left(\begin{array}{c} y_{i}^{k}\\ x_{i}^{k} \end{array}\
right)( y i k   x i k    ) the noisy iterates are a sum of === Rk
(yi0xi0)R^{k}\left(\begin{array}{c} y_{i}^{0}\\ x_{i}^{0} \end{array}
\right)R k  ( y i 0   x i 0    ) the noiseless, deterministic
iterates and +++ [?]ik[?]j=1kRk-j(1-a)\epsilon^k_i \sum_{j=1}^{k}R^{k-j}\
left(\begin{array}{c} 1\\ -\alpha \end{array}\right)[?] i k   j=1 [?] k  
R k-j  ( 1 -a  ) a decaying sum of the errors, where [?]k=Q[?]error(wk)\
epsilon^k = Q \cdot \text{error}(w^k)[?] k  =Q[?]error(w k  ).

The error term, [?]k\epsilon^k[?] k  , with its dependence on the wkw^kw 
k  , is a fairly hairy object. Following [10], we model this as
independent 0-mean Gaussian noise. In this simplified model, the
objective also breaks into two separable components, a sum of a
deterministic error and a stochastic error ^11, visualized here.

We decompose the expected value of the objective value Ef(w)-f(w[?])\
mathbf{E} f(w) - f(w^\star)Ef(w)-f(w [?]  ) into a deterministic part
and a stochastic part . Ef(w)-f(w[?])\mathbf{E} f(w) - f(w^\star) Ef(w)
-f(w [?]  ) The small black dots are a single run of stochastic
gradient
Step-size a =
Momentum b =
As [1] observes, the optimization has two phases. In the initial
transient phase the magnitude of the noise is smaller than the
magnitude of the gradient, and Momentum still makes good progress. In
the second, stochastic phase, the noise overwhelms the gradient, and
momentum is less effective.

Note that there are a set of unfortunate tradeoffs which seem to pit
the two components of error against each other. Lowering the
step-size, for example, decreases the stochastic error, but also
slows down the rate of convergence. And increasing momentum, contrary
to popular belief, causes the errors to compound. Despite these
undesirable properties, stochastic gradient descent with momentum has
still been shown to have competitive performance on neural networks.
As [1] has observed, the transient phase seems to matter more than
the fine-tuning phase in machine learning. And in fact, it has been
recently suggested [12] that this noise is a good thing -- it acts as
a implicit regularizer, which, like early stopping, prevents
overfitting in the fine-tuning phase of optimization.

---------------------------------------------------------------------

Onwards and Downwards

The study of acceleration is seeing a small revival within the
optimization community. If the ideas in this article excite you, you
may wish to read [13], which fully explores the idea of momentum as
the discretization of a certain differential equation. But other,
less physical, interpretations exist. There is an algebraic
interpretation of momentum in terms of approximating polynomials [3,
14]. Geometric interpretations are emerging [15, 16], connecting
momentum to older methods, like the Ellipsoid method. And finally,
there are interpretations relating momentum to duality [17], perhaps
providing a clue as how to accelerate second order methods and Quasi
Newton (for a first step, see [18]). But like the proverbial blind
men feeling an elephant, momentum seems like something bigger than
the sum of its parts. One day, hopefully soon, the many perspectives
will converge into a satisfying whole.

Acknowledgments

I am deeply indebted to the editorial contributions of Shan Carter
and Chris Olah, without which this article would be greatly
impoverished. Shan Carter provided complete redesigns of many of my
original interactive widgets, a visual coherence for all the figures,
and valuable optimizations to the page's performance. Chris Olah
provided impeccable editorial feedback at all levels of detail and
abstraction - from the structure of the content, to the alignment of
equations.

I am also grateful to Michael Nielsen for providing the title of this
article, which really tied the article together. Marcos Ginestra
provided editorial input for the earliest drafts of this article, and
spiritual encouragement when I needed it the most. And my gratitude
extends to my reviewers, Matt Hoffman and Anonymous Reviewer B for
their astute observations and criticism. I would like to thank
Reviewer B, in particular, for pointing out two non-trivial errors in
the original manuscript (discussion here). The contour plotting
library for the hero visualization is the joint work of Ben
Frederickson, Jeff Heer and Mike Bostock.

Many thanks to the numerous pull requests and issues filed on github.
Thanks in particular, to Osemwaro Pedro for spotting an off by one
error in one of the equations. And also to Dan Schmidt who did an
editing pass over the whole project, correcting numerous
typographical and grammatical errors.

Discussion and Review

Reviewer A - Matt Hoffman
Reviewer B - Anonymous
Discussion with User derifatives

Footnotes

 1. It is possible, however, to construct very specific
    counterexamples where momentum does not converge, even on convex
    functions. See [4] for a counterexample.
 2. In Tikhonov Regression we add a quadratic penalty to the
    regression, minimizing minimize12[?]Zw-d[?]2+e2[?]w[?]2=12wT(ZTZ+eI)w-
    (Zd)Tw \text{minimize}\qquad\tfrac{1}{2}\|Zw-d\|^{2}+\frac{\eta}
    {2}\|w\|^{2}=\tfrac{1}{2}w^{T}(Z^{T}Z+\eta I)w-(Zd)^{T}w minimize
     2  1  [?]Zw-d[?] 2  + 2  e  [?]w[?] 2  = 2  1  w T  (Z T  Z+eI)w-(Zd) T 
     w Recall that ZTZ=Q diag(L1,...,Ln) QTZ^{T}Z=Q\ \text{diag}(\
    Lambda_{1},\ldots,\Lambda_{n})\ Q^TZ T  Z=Q diag(L 1  ,...,L n  ) Q
     T  . The solution to Tikhonov Regression is therefore (ZTZ+eI)-1
    (Zd)=Q diag(1l1+e,[?],1ln+e)QT(Zd) (Z^{T}Z+\eta I)^{-1}(Zd)=Q\ \
    text{diag}\left(\frac{1}{\lambda_{1}+\eta},\cdots,\frac{1}{\
    lambda_{n}+\eta}\right)Q^T(Zd) (Z T  Z+eI) -1  (Zd)=Q diag( l 1  
    +e  1  ,[?], l n  +e  1  )Q T  (Zd) We can think of regularization
    as a function which decays the largest eigenvalues, as follows: 
    Tikhonov Regularized li=1li+e=1li(1-(1+li/e)-1). \text{Tikhonov
    Regularized } \lambda_i = \frac{1}{\lambda_{i}+\eta}=\frac{1}{\
    lambda_{i}}\left(1-\left(1+\lambda_{i}/\eta\right)^{-1}\right). 
    Tikhonov Regularized l i  = l i  +e  1  = l i    1  (1-(1+l i  /e
    ) -1  ). Gradient descent can be seen as employing a similar
    decay, but with the decay rate  Gradient Descent Regularized li=
    1li(1-(1-ali)k) \text{ Gradient Descent Regularized } \lambda_i =
    \frac{1}{\lambda_i} \left( 1-\left(1-\alpha\lambda_{i}\right)^{k}
    \right) Gradient Descent Regularized l i  = l i    1  (1-(1-al i 
     ) k  ) instead. Note that this decay is dependent on the
    step-size.
 3. This is true as we can write updates in matrix form as (10a1)
    (yik+1xik+1)=(bli01)(yikxik) \left(\!\!\begin{array}{cc} 1 & 0\\
    \alpha & 1 \end{array}\!\!\right)\Bigg(\!\!\begin{array}{c} y_{i}
    ^{k+1}\\ x_{i}^{k+1} \end{array}\!\!\Bigg)=\left(\!\!\begin
    {array}{cc} \beta & \lambda_{i}\\ 0 & 1 \end{array}\!\!\right)\
    left(\!\!\begin{array}{c} y_{i}^{k}\\ x_{i}^{k} \end{array}\!\!\
    right) ( 1 a   0 1  )( y i k+1   x i k+1    )=( b 0   l i   1  )(
     y i k   x i k    ) which implies, by inverting the matrix on the
    left, (yik+1xik+1)=(bli-ab1-ali)(yikxik)=Rk+1(xi0yi0) \Bigg(\!\!\
    begin{array}{c} y_{i}^{k+1}\\ x_{i}^{k+1} \end{array}\!\!\Bigg)=\
    left(\!\!\begin{array}{cc} \beta & \lambda_{i}\\ -\alpha\beta &
    1-\alpha\lambda_{i} \end{array}\!\!\right)\left(\!\!\begin{array}
    {c} y_{i}^{k}\\ x_{i}^{k} \end{array}\!\!\right)=R^{k+1}\left(\!
    \!\begin{array}{c} x_{i}^{0}\\ y_{i}^{0} \end{array}\!\!\right) (
     y i k+1   x i k+1    )=( b -ab   l i   1-al i    )( y i k   x i 
    k    )=R k+1  ( x i 0   y i 0    )
 4. We can write out the convergence rates explicitly. The
    eigenvalues are s1=12(1-al+b+(-al+b+1)2-4b)s2=12(1-al+b-(-al+b+1)
    2-4b) \begin{aligned} \sigma_{1} & =\frac{1}{2}\left(1-\alpha\
    lambda+\beta+\sqrt{(-\alpha\lambda+\beta+1)^{2}-4\beta}\right)\\
    [0.6em] \sigma_{2} & =\frac{1}{2}\left(1-\alpha\lambda+\beta-\
    sqrt{(-\alpha\lambda+\beta+1)^{2}-4\beta}\right) \end{aligned}  s
     1   s 2     = 2  1  (1-al+b+[?] (-al+b+1) 2  -4b   ) = 2  1  (1-al
    +b-[?] (-al+b+1) 2  -4b   )   When the (-al+b+1)2-4b<0(-\alpha\
    lambda+\beta+1)^{2}-4\beta<0(-al+b+1) 2  -4b<0 is less than zero,
    then the roots are complex and the convergence rate is |s1|=|s2|=
    (1-al+b)2+|(-al+b+1)2-4b|=2b \begin{aligned} |\sigma_{1}|=|\
    sigma_{2}| & =\sqrt{(1-\alpha\lambda+\beta)^{2}+|(-\alpha\lambda+
    \beta+1)^{2}-4\beta|}=2\sqrt{\beta} \end{aligned}  |s 1  |=|s 2  
    |   =[?] (1-al+b) 2  +|(-al+b+1) 2  -4b|   =2[?] b      Which is,
    surprisingly, independent of the step-size or the eigenvalue al\
    alpha\lambdaal. When the roots are real, the convergence rate is 
    max{|s1|,|s2|}=12max{|1-ali+b+-(1-ali+b)2-4b|} \max\{|\sigma_{1}|,
    |\sigma_{2}|\}=\tfrac{1}{2}\max\left\{ |1-\alpha\lambda_{i}+\beta
    \pm\sqrt{(1-\alpha\lambda_{i}+\beta)^{2}-4\beta}|\right\} max{|s 
    1  |,|s 2  |}= 2  1  max{|1-al i  +b+-[?] (1-al i  +b) 2  -4b   |}
 5. This can be derived by reducing the inequalities for all 4 + 1
    cases in the explicit form of the convergence rate above.
 6. We must optimize over mina,bmax{[?](bli-ab1-ali)[?],...,[?](bln-ab1-aln)
    [?]}. \min_{\alpha,\beta}\max\left\{ \bigg\| \! \left(\begin{array}
    {cc} \beta & \lambda_{i}\\ -\alpha\beta & 1-\alpha\lambda_{i} \
    end{array}\right) \! \bigg\|,\ldots,\bigg\| \! \left(\begin
    {array}{cc} \beta & \lambda_{n}\\ -\alpha\beta & 1-\alpha\lambda_
    {n} \end{array}\right)\! \bigg\|\right\}.  a,b min  max{ [?] [?] [?] [?] 
     ( b -ab   l i   1-al i    ) [?] [?] [?] [?]  ,..., [?] [?] [?] [?]  ( b -ab   l n 
      1-al n    ) [?] [?] [?] [?]  }. ( [?][?][?]\|\cdot \|[?][?][?] here denotes the
    magnitude of the maximum eigenvalue), and occurs when the roots
    of the characteristic polynomial are repeated for the matrices
    corresponding to the extremal eigenvalues.
 7. The above optimization problem is bounded from below by 000, and
    vector of all 111's achieve this.
 8. This can be written explicitly as [LG]ij={degree of vertex ii=
    j-1i[?]j,(i,j) or (j,i)[?]E0otherwise [L_{G}]_{ij}=\begin{cases} \
    text{degree of vertex }i & i=j\\ -1 & i\neq j,(i,j)\text{ or }
    (j,i)\in E\\ 0 & \text{otherwise} \end{cases} [L G  ] ij  = [?] [?] [?]
     [?] [?]   degree of vertex i -1 0   i=j i[?]j,(i,j) or (j,i)[?]E 
    otherwise  
 9. We use the infinity norm to measure our error, similar results
    can be derived for the 1 and 2 norms.
10. The momentum iterations are zk+1=bzk+Awk+error(wk)wk+1=wk-azk+1.
    \begin{aligned} z^{k+1}&=\beta z^{k}+ A w^{k} + \text{error}(w^k)
    \\[0.4em] w^{k+1}&=w^{k}-\alpha z^{k+1}. \end{aligned}  z k+1   w
     k+1     =bz k  +Aw k  +error(w k  ) =w k  -az k+1  .   which,
    after a change of variables, become (10a1)(yik+1xik+1)=(bli01)
    (yikxik)+([?]ik0) \left(\!\!\begin{array}{cc} 1 & 0\\ \alpha & 1 \
    end{array}\!\!\right)\Bigg(\!\!\begin{array}{c} y_{i}^{k+1}\\ x_
    {i}^{k+1} \end{array}\!\!\Bigg)=\left(\!\!\begin{array}{cc} \beta
    & \lambda_{i}\\ 0 & 1 \end{array}\!\!\right)\left(\!\!\begin
    {array}{c} y_{i}^{k}\\ x_{i}^{k} \end{array}\!\!\right)+\left(\!
    \!\begin{array}{c} \epsilon_{i}^{k}\\ 0 \end{array}\!\!\right) ( 
    1 a   0 1  )( y i k+1   x i k+1    )=( b 0   l i   1  )( y i k   
    x i k    )+( [?] i k   0  ) Inverting the 2x22 \times 22x2 matrix
    on the left, and applying the formula recursively yields the
    final solution.
11. On the 1D function f(x)=l2x2f(x)=\frac{\lambda}{2}x^{2}f(x)= 2  l
      x 2  , the objective value is Ef(xk)=l2E[(xk)2]=l2E(e2TRk(y0x0)
    +[?]ke2T[?]i=1kRk-i(1-a))2=l2e2TRk(y0x0)+l2E([?]ke2T[?]i=1kRk-i(1-a))2=
    l2e2TRk(y0x0)+l2E[[?]k][?][?]i=1k(e2TRk-i(1-a))2=l2e2TRk(y0x0)+lE
    [[?]k2[?][?]i=1kgi2,gi=e2TRk-i(1-a) \begin{aligned} \mathbf{E}f(x^{k})&
    =\frac{\lambda}{2}\mathbf{E}[(x^{k})^{2}]\\&=\frac{\lambda}{2}\
    mathbf{E}\left(e_{2}^{T}R^{k}\left(\begin{array}{c} y^{0}\\ x^{0}
    \end{array}\right)+\epsilon^{k}e_{2}^{T}\sum_{i=1}^{k}R^{k-i}\
    left(\begin{array}{c} 1\\ -\alpha \end{array}\right)\right)^{2}\\
    &=\frac{\lambda}{2}e_{2}^{T}R^{k}\left(\begin{array}{c} y^{0}\\ x
    ^{0} \end{array}\right)+\frac{\lambda}{2}\mathbf{E}\left(\epsilon
    ^{k}e_{2}^{T}\sum_{i=1}^{k}R^{k-i}\left(\begin{array}{c} 1\\ -\
    alpha \end{array}\right)\right)^{2}\\&=\frac{\lambda}{2}e_{2}^{T}
    R^{k}\left(\begin{array}{c} y^{0}\\ x^{0} \end{array}\right)+\
    frac{\lambda}{2}\mathbf{E}[\epsilon^{k}]\,\cdot\,\sum_{i=1}^{k}\
    left(e_{2}^{T}R^{k-i}\left(\begin{array}{c} 1\\ -\alpha \end
    {array}\right)\right)^{2}\\&=\frac{\lambda}{2}e_{2}^{T}R^{k}\left
    (\begin{array}{c} y^{0}\\ x^{0} \end{array}\right)+\frac{\lambda\
    mathbf{E}[\epsilon^{k}}{2}\cdot\sum_{i=1}^{k}\gamma_{i}^{2}, \
    qquad \gamma_i = e_{2}^{T}R^{k-i}\left(\begin{array}{c} 1\\ -\
    alpha \end{array}\right) \end{aligned}  Ef(x k  )       = 2  l  E
    [(x k  ) 2  ] = 2  l  E(e 2 T  R k  ( y 0   x 0    )+[?] k  e 2 T  
     i=1 [?] k  R k-i  ( 1 -a  )) 2   = 2  l  e 2 T  R k  ( y 0   x 0  
      )+ 2  l  E([?] k  e 2 T   i=1 [?] k  R k-i  ( 1 -a  )) 2   = 2  l  
    e 2 T  R k  ( y 0   x 0    )+ 2  l  E[[?] k  ][?] i=1 [?] k  (e 2 T  R 
    k-i  ( 1 -a  )) 2   = 2  l  e 2 T  R k  ( y 0   x 0    )+ 2  lE[[?]
     k    [?] i=1 [?] k  g i 2  ,g i  =e 2 T  R k-i  ( 1 -a  )   The
    third inequality uses the fact that E[?]k=0\mathbf{E} \epsilon^k =
    0E[?] k  =0 and the fourth uses the fact they are uncorrelated.

References

 1. On the importance of initialization and momentum in deep
    learning.  [PDF]
    Sutskever, I., Martens, J., Dahl, G.E. and Hinton, G.E., 2013.
    ICML (3), Vol 28, pp. 1139--1147.
 2. Some methods of speeding up the convergence of iteration methods
     [PDF]
    Polyak, B.T., 1964. USSR Computational Mathematics and
    Mathematical Physics, Vol 4(5), pp. 1--17. Elsevier. DOI: 10.1016/
    0041-5553(64)90137-5
 3. Theory of gradient methods
    Rutishauser, H., 1959. Refined iterative methods for computation
    of the solution and the eigenvalues of self-adjoint boundary
    value problems, pp. 24--49. Springer. DOI: 10.1007/
    978-3-0348-7224-9_2
 4. Analysis and design of optimization algorithms via integral
    quadratic constraints  [PDF]
    Lessard, L., Recht, B. and Packard, A., 2016. SIAM Journal on
    Optimization, Vol 26(1), pp. 57--95. SIAM.
 5. Introductory lectures on convex optimization: A basic course
    Nesterov, Y., 2013. , Vol 87. Springer Science \& Business Media.
    DOI: 10.1007/978-1-4419-8853-9
 6. Natural gradient works efficiently in learning  [link]
    Amari, S., 1998. Neural computation, Vol 10(2), pp. 251--276. MIT
    Press. DOI: 10.1162/089976698300017746
 7. Deep Learning, NIPS'2015 Tutorial  [PDF]
    Hinton, G., Bengio, Y. and LeCun, Y., 2015.
 8. Adaptive restart for accelerated gradient schemes  [PDF]
    O'Donoghue, B. and Candes, E., 2015. Foundations of computational
    mathematics, Vol 15(3), pp. 715--732. Springer. DOI: 10.1007/
    s10208-013-9150-3
 9. The Nth Power of a 2x2 Matrix.  [PDF]
    Williams, K., 1992. Mathematics Magazine, Vol 65(5), pp. 336.
    MAA. DOI: 10.2307/2691246
10. From Averaging to Acceleration, There is Only a Step-size.  [PDF]
    Flammarion, N. and Bach, F.R., 2015. COLT, pp. 658--695.
11. On the momentum term in gradient descent learning algorithms  
    [PDF]
    Qian, N., 1999. Neural networks, Vol 12(1), pp. 145--151.
    Elsevier. DOI: 10.1016/s0893-6080(98)00116-6
12. Understanding deep learning requires rethinking generalization  
    [PDF]
    Zhang, C., Bengio, S., Hardt, M., Recht, B. and Vinyals, O.,
    2016. arXiv preprint arXiv:1611.03530.
13. A differential equation for modeling Nesterov's accelerated
    gradient method: Theory and insights  [PDF]
    Su, W., Boyd, S. and Candes, E., 2014. Advances in Neural
    Information Processing Systems, pp. 2510--2518.
14. The Zen of Gradient Descent  [HTML]
    Hardt, M., 2013.
15. A geometric alternative to Nesterov's accelerated gradient
    descent  [PDF]
    Bubeck, S., Lee, Y.T. and Singh, M., 2015. arXiv preprint
    arXiv:1506.08187.
16. An optimal first order method based on optimal quadratic
    averaging  [PDF]
    Drusvyatskiy, D., Fazel, M. and Roy, S., 2016. arXiv preprint
    arXiv:1604.06543.
17. Linear coupling: An ultimate unification of gradient and mirror
    descent  [PDF]
    Allen-Zhu, Z. and Orecchia, L., 2014. arXiv preprint
    arXiv:1407.1537.
18. Accelerating the cubic regularization of Newton's method on
    convex problems  [PDF]
    Nesterov, Y., 2008. Mathematical Programming, Vol 112(1), pp.
    159--181. Springer. DOI: 10.1007/s10107-006-0089-x

Updates and Corrections

View all changes to this article since it was first published. If you
see a mistake or want to suggest a change, please create an issue on
GitHub.

Citations and Reuse

Diagrams and text are licensed under Creative Commons Attribution
CC-BY 2.0, unless noted otherwise, with the source available on
GitHub. The figures that have been reused from other sources don't
fall under this license and can be recognized by a note in their
caption: "Figure from ...".

For attribution in academic contexts, please cite this work as

Goh, "Why Momentum Really Works", Distill, 2017. http://doi.org/10.23915/distill.00006

BibTeX citation

@article{goh2017why,
  author = {Goh, Gabriel},
  title = {Why Momentum Really Works},
  journal = {Distill},
  year = {2017},
  url = {http://distill.pub/2017/momentum},
  doi = {10.23915/distill.00006}
}

On the importance of initialization and momentum in deep learning.  
[PDF]
I. Sutskever, J. Martens, G.E. Dahl, G.E. Hinton.
ICML (3), Vol 28, pp. 1139--1147. 2013.

Some methods of speeding up the convergence of iteration methods  
[PDF]
B.T. Polyak.
USSR Computational Mathematics and Mathematical Physics, Vol 4(5),
pp. 1--17. Elsevier. 1964.
DOI: 10.1016/0041-5553(64)90137-5

Theory of gradient methods
H. Rutishauser.
Refined iterative methods for computation of the solution and the
eigenvalues of self-adjoint boundary value problems, pp. 24--49.
Springer. 1959.
DOI: 10.1007/978-3-0348-7224-9_2
Analysis and design of optimization algorithms via integral quadratic
constraints  [PDF]
L. Lessard, B. Recht, A. Packard.
SIAM Journal on Optimization, Vol 26(1), pp. 57--95. SIAM. 2016.
Introductory lectures on convex optimization: A basic course
Y. Nesterov.
, Vol 87. Springer Science \& Business Media. 2013.
DOI: 10.1007/978-1-4419-8853-9
Natural gradient works efficiently in learning  [link]
S. Amari.
Neural computation, Vol 10(2), pp. 251--276. MIT Press. 1998.
DOI: 10.1162/089976698300017746
Deep Learning, NIPS'2015 Tutorial  [PDF]
G. Hinton, Y. Bengio, Y. LeCun.
2015.
Adaptive restart for accelerated gradient schemes  [PDF]
B. O'Donoghue, E. Candes.
Foundations of computational mathematics, Vol 15(3), pp. 715--732.
Springer. 2015.
DOI: 10.1007/s10208-013-9150-3
The Nth Power of a 2x2 Matrix.  [PDF]
K. Williams.
Mathematics Magazine, Vol 65(5), pp. 336. MAA. 1992.
DOI: 10.2307/2691246
From Averaging to Acceleration, There is Only a Step-size.  [PDF]
N. Flammarion, F.R. Bach.
COLT, pp. 658--695. 2015.
On the momentum term in gradient descent learning algorithms  [PDF]
N. Qian.
Neural networks, Vol 12(1), pp. 145--151. Elsevier. 1999.
DOI: 10.1016/s0893-6080(98)00116-6
Introductory lectures on convex optimization: A basic course
Y. Nesterov.
, Vol 87. Springer Science \& Business Media. 2013.
DOI: 10.1007/978-1-4419-8853-9
From Averaging to Acceleration, There is Only a Step-size.  [PDF]
N. Flammarion, F.R. Bach.
COLT, pp. 658--695. 2015.
On the importance of initialization and momentum in deep learning.  
[PDF]
I. Sutskever, J. Martens, G.E. Dahl, G.E. Hinton.
ICML (3), Vol 28, pp. 1139--1147. 2013.
On the importance of initialization and momentum in deep learning.  
[PDF]
I. Sutskever, J. Martens, G.E. Dahl, G.E. Hinton.
ICML (3), Vol 28, pp. 1139--1147. 2013.
Understanding deep learning requires rethinking generalization  [PDF]
C. Zhang, S. Bengio, M. Hardt, B. Recht, O. Vinyals.
arXiv preprint arXiv:1611.03530. 2016.
A differential equation for modeling Nesterov's accelerated gradient
method: Theory and insights  [PDF]
W. Su, S. Boyd, E. Candes.
Advances in Neural Information Processing Systems, pp. 2510--2518.
2014.
Theory of gradient methods
H. Rutishauser.
Refined iterative methods for computation of the solution and the
eigenvalues of self-adjoint boundary value problems, pp. 24--49.
Springer. 1959.
DOI: 10.1007/978-3-0348-7224-9_2

The Zen of Gradient Descent  [HTML]
M. Hardt. 2013.
A geometric alternative to Nesterov's accelerated gradient descent  
[PDF]
S. Bubeck, Y.T. Lee, M. Singh.
arXiv preprint arXiv:1506.08187. 2015.

An optimal first order method based on optimal quadratic averaging  
[PDF]
D. Drusvyatskiy, M. Fazel, S. Roy.
arXiv preprint arXiv:1604.06543. 2016.
Linear coupling: An ultimate unification of gradient and mirror
descent  [PDF]
Z. Allen-Zhu, L. Orecchia.
arXiv preprint arXiv:1407.1537. 2014.
Accelerating the cubic regularization of Newton's method on convex
problems  [PDF]
Y. Nesterov.
Mathematical Programming, Vol 112(1), pp. 159--181. Springer. 2008.
DOI: 10.1007/s10107-006-0089-x
It is possible, however, to construct very specific counterexamples
where momentum does not converge, even on convex functions. See [4]
for a counterexample.
In Tikhonov Regression we add a quadratic penalty to the regression,
minimizing minimize12[?]Zw-d[?]2+e2[?]w[?]2=12wT(ZTZ+eI)w-(Zd)Tw \text
{minimize}\qquad\tfrac{1}{2}\|Zw-d\|^{2}+\frac{\eta}{2}\|w\|^{2}=\
tfrac{1}{2}w^{T}(Z^{T}Z+\eta I)w-(Zd)^{T}w minimize 2  1  [?]Zw-d[?] 2  +
 2  e  [?]w[?] 2  = 2  1  w T  (Z T  Z+eI)w-(Zd) T  w Recall that ZTZ=
Q diag(L1,...,Ln) QTZ^{T}Z=Q\ \text{diag}(\Lambda_{1},\ldots,\Lambda_
{n})\ Q^TZ T  Z=Q diag(L 1  ,...,L n  ) Q T  . The solution to Tikhonov
Regression is therefore (ZTZ+eI)-1(Zd)=Q diag(1l1+e,[?],1ln+e)QT(Zd) (Z
^{T}Z+\eta I)^{-1}(Zd)=Q\ \text{diag}\left(\frac{1}{\lambda_{1}+\
eta},\cdots,\frac{1}{\lambda_{n}+\eta}\right)Q^T(Zd) (Z T  Z+eI) -1  
(Zd)=Q diag( l 1  +e  1  ,[?], l n  +e  1  )Q T  (Zd) We can think of
regularization as a function which decays the largest eigenvalues, as
follows: Tikhonov Regularized li=1li+e=1li(1-(1+li/e)-1). \text
{Tikhonov Regularized } \lambda_i = \frac{1}{\lambda_{i}+\eta}=\frac
{1}{\lambda_{i}}\left(1-\left(1+\lambda_{i}/\eta\right)^{-1}\right). 
Tikhonov Regularized l i  = l i  +e  1  = l i    1  (1-(1+l i  /e) -1
  ). Gradient descent can be seen as employing a similar decay, but
with the decay rate  Gradient Descent Regularized li=1li(1-(1-ali)k)
\text{ Gradient Descent Regularized } \lambda_i = \frac{1}{\lambda_i}
\left( 1-\left(1-\alpha\lambda_{i}\right)^{k} \right)
 Gradient Descent Regularized l i  = l i    1  (1-(1-al i  ) k  )
instead. Note that this decay is dependent on the step-size.
This is true as we can write updates in matrix form as (10a1)
(yik+1xik+1)=(bli01)(yikxik) \left(\!\!\begin{array}{cc} 1 & 0\\ \
alpha & 1 \end{array}\!\!\right)\Bigg(\!\!\begin{array}{c} y_{i}^
{k+1}\\ x_{i}^{k+1} \end{array}\!\!\Bigg)=\left(\!\!\begin{array}{cc}
\beta & \lambda_{i}\\ 0 & 1 \end{array}\!\!\right)\left(\!\!\begin
{array}{c} y_{i}^{k}\\ x_{i}^{k} \end{array}\!\!\right) ( 1 a   0 1  
)( y i k+1   x i k+1    )=( b 0   l i   1  )( y i k   x i k    )
which implies, by inverting the matrix on the left, (yik+1xik+1)=
(bli-ab1-ali)(yikxik)=Rk+1(xi0yi0) \Bigg(\!\!\begin{array}{c} y_{i}^
{k+1}\\ x_{i}^{k+1} \end{array}\!\!\Bigg)=\left(\!\!\begin{array}{cc}
\beta & \lambda_{i}\\ -\alpha\beta & 1-\alpha\lambda_{i} \end{array}
\!\!\right)\left(\!\!\begin{array}{c} y_{i}^{k}\\ x_{i}^{k} \end
{array}\!\!\right)=R^{k+1}\left(\!\!\begin{array}{c} x_{i}^{0}\\ y_
{i}^{0} \end{array}\!\!\right) ( y i k+1   x i k+1    )=( b -ab   l i
   1-al i    )( y i k   x i k    )=R k+1  ( x i 0   y i 0    )
We can write out the convergence rates explicitly. The eigenvalues
are s1=12(1-al+b+(-al+b+1)2-4b)s2=12(1-al+b-(-al+b+1)2-4b) \begin
{aligned} \sigma_{1} & =\frac{1}{2}\left(1-\alpha\lambda+\beta+\sqrt
{(-\alpha\lambda+\beta+1)^{2}-4\beta}\right)\\[0.6em] \sigma_{2} & =\
frac{1}{2}\left(1-\alpha\lambda+\beta-\sqrt{(-\alpha\lambda+\beta+1)^
{2}-4\beta}\right) \end{aligned}  s 1   s 2     = 2  1  (1-al+b+[?] (-a
l+b+1) 2  -4b   ) = 2  1  (1-al+b-[?] (-al+b+1) 2  -4b   )   When the 
(-al+b+1)2-4b<0(-\alpha\lambda+\beta+1)^{2}-4\beta<0(-al+b+1) 2  -4b<
0 is less than zero, then the roots are complex and the convergence
rate is |s1|=|s2|=(1-al+b)2+|(-al+b+1)2-4b|=2b \begin{aligned} |\
sigma_{1}|=|\sigma_{2}| & =\sqrt{(1-\alpha\lambda+\beta)^{2}+|(-\
alpha\lambda+\beta+1)^{2}-4\beta|}=2\sqrt{\beta} \end{aligned}  |s 1 
 |=|s 2  |   =[?] (1-al+b) 2  +|(-al+b+1) 2  -4b|   =2[?] b      Which
is, surprisingly, independent of the step-size or the eigenvalue al\
alpha\lambdaal. When the roots are real, the convergence rate is max
{|s1|,|s2|}=12max{|1-ali+b+-(1-ali+b)2-4b|} \max\{|\sigma_{1}|,|\
sigma_{2}|\}=\tfrac{1}{2}\max\left\{ |1-\alpha\lambda_{i}+\beta\pm\
sqrt{(1-\alpha\lambda_{i}+\beta)^{2}-4\beta}|\right\} max{|s 1  |,|s 
2  |}= 2  1  max{|1-al i  +b+-[?] (1-al i  +b) 2  -4b   |}
This can be derived by reducing the inequalities for all 4 + 1 cases
in the explicit form of the convergence rate above.
We must optimize over mina,bmax{[?](bli-ab1-ali)[?],...,[?](bln-ab1-aln)[?]}. \
min_{\alpha,\beta}\max\left\{ \bigg\| \! \left(\begin{array}{cc} \
beta & \lambda_{i}\\ -\alpha\beta & 1-\alpha\lambda_{i} \end{array}\
right) \! \bigg\|,\ldots,\bigg\| \! \left(\begin{array}{cc} \beta & \
lambda_{n}\\ -\alpha\beta & 1-\alpha\lambda_{n} \end{array}\right)\!
\bigg\|\right\}.  a,b min  max{ [?] [?] [?] [?]  ( b -ab   l i   1-al i    ) 
[?] [?] [?] [?]  ,..., [?] [?] [?] [?]  ( b -ab   l n   1-al n    ) [?] [?] [?] [?]  }. ( [?][?][?]\|
\cdot \|[?][?][?] here denotes the magnitude of the maximum eigenvalue),
and occurs when the roots of the characteristic polynomial are
repeated for the matrices corresponding to the extremal eigenvalues.
The above optimization problem is bounded from below by 000, and
vector of all 111's achieve this.
This can be written explicitly as [LG]ij={degree of vertex ii=j-1i[?]j,
(i,j) or (j,i)[?]E0otherwise [L_{G}]_{ij}=\begin{cases} \text{degree of
vertex }i & i=j\\ -1 & i\neq j,(i,j)\text{ or }(j,i)\in E\\ 0 & \text
{otherwise} \end{cases} [L G  ] ij  = [?] [?] [?] [?] [?]   degree of vertex i 
-1 0   i=j i[?]j,(i,j) or (j,i)[?]E otherwise  
We use the infinity norm to measure our error, similar results can be
derived for the 1 and 2 norms.
The momentum iterations are zk+1=bzk+Awk+error(wk)wk+1=wk-azk+1. \
begin{aligned} z^{k+1}&=\beta z^{k}+ A w^{k} + \text{error}(w^k) \\
[0.4em] w^{k+1}&=w^{k}-\alpha z^{k+1}. \end{aligned}  z k+1   w k+1  
   =bz k  +Aw k  +error(w k  ) =w k  -az k+1  .   which, after a
change of variables, become (10a1)(yik+1xik+1)=(bli01)(yikxik)+([?]ik0)
\left(\!\!\begin{array}{cc} 1 & 0\\ \alpha & 1 \end{array}\!\!\right)
\Bigg(\!\!\begin{array}{c} y_{i}^{k+1}\\ x_{i}^{k+1} \end{array}\!\!\
Bigg)=\left(\!\!\begin{array}{cc} \beta & \lambda_{i}\\ 0 & 1 \end
{array}\!\!\right)\left(\!\!\begin{array}{c} y_{i}^{k}\\ x_{i}^{k} \
end{array}\!\!\right)+\left(\!\!\begin{array}{c} \epsilon_{i}^{k}\\ 0
\end{array}\!\!\right) ( 1 a   0 1  )( y i k+1   x i k+1    )=( b 0  
 l i   1  )( y i k   x i k    )+( [?] i k   0  ) Inverting the 2x22 \
times 22x2 matrix on the left, and applying the formula recursively
yields the final solution.
On the 1D function f(x)=l2x2f(x)=\frac{\lambda}{2}x^{2}f(x)= 2  l  x 
2  , the objective value is Ef(xk)=l2E[(xk)2]=l2E(e2TRk(y0x0)+[?]ke2T[?]i
=1kRk-i(1-a))2=l2e2TRk(y0x0)+l2E([?]ke2T[?]i=1kRk-i(1-a))2=l2e2TRk(y0x0)
+l2E[[?]k][?][?]i=1k(e2TRk-i(1-a))2=l2e2TRk(y0x0)+lE[[?]k2[?][?]i=1kgi2,gi=
e2TRk-i(1-a) \begin{aligned} \mathbf{E}f(x^{k})&=\frac{\lambda}{2}\
mathbf{E}[(x^{k})^{2}]\\&=\frac{\lambda}{2}\mathbf{E}\left(e_{2}^{T}R
^{k}\left(\begin{array}{c} y^{0}\\ x^{0} \end{array}\right)+\epsilon^
{k}e_{2}^{T}\sum_{i=1}^{k}R^{k-i}\left(\begin{array}{c} 1\\ -\alpha \
end{array}\right)\right)^{2}\\&=\frac{\lambda}{2}e_{2}^{T}R^{k}\left
(\begin{array}{c} y^{0}\\ x^{0} \end{array}\right)+\frac{\lambda}{2}\
mathbf{E}\left(\epsilon^{k}e_{2}^{T}\sum_{i=1}^{k}R^{k-i}\left(\begin
{array}{c} 1\\ -\alpha \end{array}\right)\right)^{2}\\&=\frac{\
lambda}{2}e_{2}^{T}R^{k}\left(\begin{array}{c} y^{0}\\ x^{0} \end
{array}\right)+\frac{\lambda}{2}\mathbf{E}[\epsilon^{k}]\,\cdot\,\
sum_{i=1}^{k}\left(e_{2}^{T}R^{k-i}\left(\begin{array}{c} 1\\ -\alpha
\end{array}\right)\right)^{2}\\&=\frac{\lambda}{2}e_{2}^{T}R^{k}\left
(\begin{array}{c} y^{0}\\ x^{0} \end{array}\right)+\frac{\lambda\
mathbf{E}[\epsilon^{k}}{2}\cdot\sum_{i=1}^{k}\gamma_{i}^{2}, \qquad \
gamma_i = e_{2}^{T}R^{k-i}\left(\begin{array}{c} 1\\ -\alpha \end
{array}\right) \end{aligned}  Ef(x k  )       = 2  l  E[(x k  ) 2  ] 
= 2  l  E(e 2 T  R k  ( y 0   x 0    )+[?] k  e 2 T   i=1 [?] k  R k-i  (
 1 -a  )) 2   = 2  l  e 2 T  R k  ( y 0   x 0    )+ 2  l  E([?] k  e 2 
T   i=1 [?] k  R k-i  ( 1 -a  )) 2   = 2  l  e 2 T  R k  ( y 0   x 0   
 )+ 2  l  E[[?] k  ][?] i=1 [?] k  (e 2 T  R k-i  ( 1 -a  )) 2   = 2  l  e 
2 T  R k  ( y 0   x 0    )+ 2  lE[[?] k    [?] i=1 [?] k  g i 2  ,g i  =e 2
 T  R k-i  ( 1 -a  )   The third inequality uses the fact that E[?]k=0\
mathbf{E} \epsilon^k = 0E[?] k  =0 and the fourth uses the fact they
are uncorrelated.
Distill is dedicated to clear explanations of machine learning
About Submit Prize Archive RSS GitHub Twitter      ISSN 2476-0757