Neural Networks and Optimization · Study tools

InteractiveGoodfellow–Bengio–Courville, Ch. 6: Deep Feedforward Networks

Neural networks as functions

A neural network is just a family of functions with adjustable weights, and machine learning, at its core, is picking the member of the family that fits the data. The family is built from one simple principle: an affine map followed by a nonlinearity, repeated. (In reinforcement learning, this is how value functions and policies are represented once states are continuous.)

The idea

Learning is function fitting. Given data $(x_i,y_i)$ , find a function $f$ with $f(x_i)\approx y_i$ , then predict $f(x)$ at inputs you haven't seen. The family of functions you search matters: if it only contains straight lines and the data curve, every prediction will be off.

Feedforward network. Each node takes a weighted sum of every node in the previous layer plus a bias, and applies an activation

\sigma

. In matrix form, one layer is

h\mapsto\sigma(Wh+b)

, and a network with two hidden layers is

f_\theta(x)=W_3\,\sigma\big(W_2\,\sigma(W_1x+b_1)+b_2\big)+b_3,

where

\theta

collects every weight matrix and bias vector. The input layer's size is the input dimension (two nodes means

x\in\R^2

); the width is the number of nodes in a hidden layer and the depth is the number of hidden layers.

Without an activation, depth buys nothing. If

\sigma

is the identity,

W_3(W_2(W_1x+b_1)+b_2)+b_3=Mx+c

: a composition of affine maps is affine. However many layers, the network draws a plane.

Universal approximation (Cybenko 1989; Leshno et al. 1993; Goodfellow et al. §6.4.1). With any continuous activation that is not a polynomial, for example

\mathrm{ReLU}(z)=\max(0,z)

, one hidden layer suffices: for every continuous

g

on a compact set

K

and every

\varepsilon>0

there is a network with

\sup_{x\in K}|f_\theta(x)-g(x)|<\varepsilon

.

The theorem says a good network exists; it does not say how wide it must be, or how to find its weights. Smaller $\varepsilon$ generally needs more parameters, and under smoothness assumptions on $g$ one can say how many. In the figure, a ReLU network with $n$ hidden units is the piecewise-linear interpolant at $n+1$ points; for twice-differentiable $g$ its error falls like $1/n^2$ .

Explore

Target

Activation

Hidden units n

target

g

network

f_\theta

. The faint lines are the hidden units' contributions, each a scaled

\max(0,x-k)

with its kink at a knot

k

.

Worked example

A network $\R^2\to\R$ with one hidden layer of two ReLU units: $h_1=\mathrm{ReLU}(x_1-2x_2+0.1)$ , $h_2=\mathrm{ReLU}(-x_1+x_2)$ , output $y=h_1+3h_2-1$ . Evaluate it at $x=(1,-1)$ and $x=(-1,1)$ .

Forward pass at $(1,-1)$ .

$h_1=\mathrm{ReLU}(1+2+0.1)=3.1$ , $h_2=\mathrm{ReLU}(-2)=0$ . $y=3.1+0-1=2.1$ .
Forward pass at $(-1,1)$ .

$h_1=\mathrm{ReLU}(-1-2+0.1)=0$ , $h_2=\mathrm{ReLU}(2)=2$ . $y=0+6-1=5$ .
Count the parameters, and say what the network looks like as a function.

First layer: a $2\times2$ matrix and 2 biases (6); output: 2 weights and 1 bias (3). Total 9. Each ReLU is zero on one side of a line and linear on the other, so $f$ is piecewise linear with at most four pieces, cut by the lines $x_1-2x_2+0.1=0$ and $x_2=x_1$ . Large networks are the same thing with far more pieces.

Check your understanding

A network has ten hidden layers of width 100, but no activation functions. It can represent:

Any continuous function, given enough layers.Composing affine maps never leaves the affine maps.
Only affine functions $x\mapsto Mx+c$ .

The nonlinearity is what lets layers build on each other. Even the very simple

\max(0,z)

, two linear pieces glued together, is enough for universal approximation.

Universal approximation guarantees that:

Training by gradient descent will find a network within $\varepsilon$ of the target.It says nothing about training.
Some network within $\varepsilon$ of the target exists, of some size.
A network with a fixed number of units can approximate every continuous function.The size depends on $\varepsilon$ and on the target.

Existence, not construction and not size. Rates (how many parameters for a given

\varepsilon

) need extra assumptions on the target, and they are what make bigger models more expressive.

In the figure, the steep step needs many more units than sin for the same error. Why?

Piecewise-linear interpolation error is controlled by the target's slope and curvature, and the step is very steep in one place.
ReLU networks cannot represent increasing functions.Each unit is increasing; sums of them certainly can be.

With knots spaced

h

apart, the error is at most

\tfrac{h^2}8\max|g''|

, or

\tfrac h2\max|g'|

. Spending knots evenly wastes them on the flat parts; putting them where the function bends is what training (or adaptive knots) does.

Practice

The hat functionThree ReLUs make a bump

Write the hat $g(x)=\max(0,1-|x|)$ as a combination of ReLUs $\max(0,x-k)$ .

Solution

$g(x)=\mathrm{ReLU}(x+1)-2\,\mathrm{ReLU}(x)+\mathrm{ReLU}(x-1)$ . Check: for $x\le-1$ all are 0; on $[-1,0]$ it is $x+1$ ; on $[0,1]$ it is $x+1-2x=1-x$ ; for $x\ge1$ , $x+1-2x+x-1=0$ . Three hidden units. Sums of shifted hats are exactly the piecewise-linear interpolants, which proves universal approximation in one dimension.

Error versus widthA rate

Using $n$ equally spaced pieces on $[-1,1]$ , how many hidden units are needed to approximate $g(x)=\sin(3x)$ to within $10^{-3}$ ?

Solution

$h=2/n$ , $|g''|\le9$ , so the error is at most $\frac{(2/n)^2}{8}\cdot9=\frac{4.5}{n^2}$ . Need $n^2\ge4500$ : $n\ge68$ . Ten times more accuracy costs only about $\sqrt{10}\approx3.2$ times more units, because $g$ is smooth.

InteractiveGoodfellow–Bengio–Courville, §6.1–6.3

The forward pass, layer by layer

Every node does the same small job: take a weighted sum of all the nodes in the layer before, add a bias, and apply the activation. Its output becomes an input to every node in the next layer. Do that layer after layer and the input becomes the output.

The idea

Take two inputs $x_1,x_2$ and a first hidden layer of three nodes. Each node has its own weights and bias:

z^{(1)}_1=w_{11}x_1+w_{12}x_2+b_1,\quad z^{(1)}_2=w_{21}x_1+w_{22}x_2+b_2,\quad z^{(1)}_3=w_{31}x_1+w_{32}x_2+b_3,

and outputs $h^{(1)}_j=\sigma(z^{(1)}_j)$ . Stacking the weights in a matrix, the whole layer is $h^{(1)}=\sigma(W_1x+b_1)$ with $W_1$ of size $3\times2$ : row $j$ holds node $j$ 's weights. The next layer takes $h^{(1)}$ as its input in exactly the same way, and so on.

Reading the shapes. In a network with layer sizes

2\to3\to3\to2

:

x\in\R^2

;

W_1\in\R^{3\times2}

,

b_1\in\R^3

;

W_2\in\R^{3\times3}

,

b_2\in\R^3

;

W_3\in\R^{2\times3}

,

b_3\in\R^2

; output

y\in\R^2

. A weight matrix always has shape (nodes in this layer)

\times

(nodes in the previous layer). If the shapes don't fit, the formula is wrong.

Parameters. The weights and biases together, often written

\theta

. This network has

(3\cdot2+3)+(3\cdot3+3)+(2\cdot3+2)=29

. When a language model is described by its number of parameters, this is the count being reported (in far larger networks).

Training chooses $\theta$ . Once it is chosen, the network is just a fixed function $\R^2\to\R^2$ , and evaluating it is the forward pass you can step through below.

Explore

Input x₁ Input x₂

Activation in hidden layers

Click a node to see how its value is computed. Line thickness is the size of the weight; positive negative. A grey node is switched off: ReLU sent its negative sum to 0.

Worked example

The network in the figure has $W_1=\begin{pmatrix}1&-2\\0.5&1\\-1&0.5\end{pmatrix},\ b_1=\begin{pmatrix}0.1\\0.2\\-0.4\end{pmatrix}$ , $W_2=\begin{pmatrix}0&2&0.5\\1&-1&0.5\\-0.5&0&1\end{pmatrix},\ b_2=\begin{pmatrix}0\\0\\0.2\end{pmatrix}$ , $W_3=\begin{pmatrix}1&-1&0.5\\0.5&1&-1\end{pmatrix},\ b_3=\begin{pmatrix}0.1\\-0.2\end{pmatrix}$ , ReLU in the hidden layers and none at the output. Evaluate it at $x=(1,-1)$ .

First hidden layer: compute $z^{(1)}=W_1x+b_1$ , then $h^{(1)}$ .

$z^{(1)}=(1+2+0.1,\ 0.5-1+0.2,\ -1-0.5-0.4)=(3.1,\,-0.3,\,-1.9)$ . ReLU zeroes the negatives: $h^{(1)}=(3.1,0,0)$ .
Second hidden layer.

$z^{(2)}=W_2h^{(1)}+b_2=(0,\ 3.1,\ -1.55+0.2)=(0,3.1,-1.35)$ , so $h^{(2)}=(0,3.1,0)$ .
Output.

$y=W_3h^{(2)}+b_3=(-3.1+0.1,\ 3.1-0.2)=(-3.0,\ 2.9)$ . Set the sliders to $x=(1,-1)$ and check.
Now remove every activation. What single affine map is the network?

$y=W_3(W_2(W_1x+b_1)+b_2)+b_3=Mx+c$ with $M=W_3W_2W_1$ (a $2\times2$ matrix) and $c=W_3W_2b_1+W_3b_2+b_3$ . Three layers, 29 parameters, and the result is something 6 parameters already describe.

Check your understanding

A layer maps 128 nodes to 64 nodes. Its weight matrix has shape:

$128\times64$ It multiplies the previous layer's vector, which has 128 entries, so it needs 128 columns.
$64\times128$

W\in\R^{64\times128}

times

h\in\R^{128}

gives

\R^{64}

. Rows are the new nodes, columns the old ones. With the bias, the layer has

64\cdot128+64=8256

parameters.

In the figure, a hidden node is grey for the current input. If you nudge the input slightly, how does that node's output change?

Not at all: it stays at 0 while its sum is negative.
Linearly with the input.Only on the positive side of its kink.

ReLU is flat on the negative side. Which nodes are on or off splits the input space into regions, and on each region the network is a single affine map: a ReLU network is piecewise linear.

Does the bias matter? Without biases, a ReLU network $f$ always satisfies:

$f(0)=0$ and $f(cx)=c\,f(x)$ for $c>0$ .
Nothing special; biases only shift the output.Try $x=0$ in every layer.

With no biases each layer is

\sigma(Wh)

, and

\mathrm{ReLU}(cz)=c\,\mathrm{ReLU}(z)

for

c>0

. So the network can't represent even a constant function unless biases are allowed: the kinks of every unit are stuck at the origin.

Practice

Counting parametersA wider, deeper network

Count the parameters of a network $784\to256\to256\to10$ (an image classifier for 28-by-28 images). Which layer has the most?

Solution

$(784\cdot256+256)+(256\cdot256+256)+(256\cdot10+10)=200{,}960+65{,}792+2{,}570=269{,}322$ . The first layer has three quarters of them: parameter counts are dominated by the widest pair of adjacent layers.

Regions of a ReLU networkOne hidden layer in the plane

A network $\R^2\to\R$ with one hidden layer of $n$ ReLU units. Each unit is off on one side of a line. At most how many regions can the $n$ lines cut the plane into, and what does the network look like on each one?

Solution

$n$ lines in general position cut the plane into $1+n+\binom n2$ regions. On each region every unit is either off or exactly linear, so the network is affine there: a piecewise-linear surface with up to $1+n+n(n-1)/2$ flat pieces. Depth multiplies the count further, which is one explanation of why deep networks are more expressive than shallow ones of the same size.

InteractiveChen et al. (2018), Neural Ordinary Differential Equations

Neural ODEs: learning a flow

Usually you are given a differential equation and asked for its solutions. Turn it around: you are given solutions (observed trajectories of some system) and asked for the equation. A neural ODE answers with a neural network for the velocity, one network for every starting point.

The idea

The system. An ordinary differential equation

\dot z(t)=F\big(z(t),t\big)

,

z(0)=z_0

, with

z(t)\in\R^d

. Written out, it is a system of

d

equations

\dot z_i=F_i(z_1,\dots,z_d,t)

, each with a single derivative, in time. (A partial differential equation, like the heat equation

\partial_tu=\partial_x^2u

, mixes derivatives in several variables; that is the difference.)

F

is the vector field: at each point and time, the velocity.

Examples. $\dot z=t$ integrates directly: $z(t)=z_0+t^2/2$ . For $\dot z=-z$ , separate the variables: $\frac{dz}z=-dt$ , so $\log|z|=-t+c$ and $z(t)=z_0e^{-t}$ . In general one can write $z(t)=z_0+\int_0^tF(z(s),s)\,ds$ , but that is circular as a recipe: the unknown $z$ is inside the integral.

The flow map. Fix a time

T

. The map

\varphi_T:\R^d\to\R^d

,

z_0\mapsto z(T)

, sends each starting point to where its solution is at time

T

. Check what kind of object it is: a vector goes in and a vector comes out. It is not a curve; the curves appear when you let

T

vary, and the family of all

\varphi_t

is the flow.

Neural ODE. Replace the unknown

F

by a network

F_\theta(z,t)

, solve

\dot x=F_\theta(x,t)

,

x(0)=z_0

, and choose

\theta

so that the solutions

x_{z_0}(t)

match the observed ones. The same

\theta

is used for every starting point: the goal is the law of motion, not a separate fit per trajectory.

What “approximate the flow” means. For a compact (closed and bounded) set

K

of starting points and a finite horizon

T

, universal-approximation results for neural ODEs give, for every

\varepsilon>0

, parameters

\theta

with

\sup_{z_0\in K}\ \sup_{t\in[0,T]}\ \big|z_{z_0}(t)-x_{z_0}(t)\big|<\varepsilon.

One network, every starting point in

K

, every time up to

T

. The error goes to 0 as the number of parameters grows; with more smoothness of

F

, one can say how fast.

Sup, not max. $\sup$ is the least upper bound, which always exists for a bounded set of real numbers; a maximum need not. On $(0,1)$ , $f(x)=x$ has supremum 1 but no maximum, since 1 is never reached. When a maximum exists the two agree, so writing $\sup$ is simply the safe choice.

Variants. The network can take time as an input in different ways: weights that are functions of time (or piecewise constant in time), or, in semi-autonomous neural ODEs, one fixed set of weights with time entering only linearly inside the activation, $F_\theta(z,t)=W\sigma(A_1z+a_2t+b)$ . Check the shapes: $A_1$ is square, $z$ , $a_2$ , $b$ are column vectors and $t$ is a scalar, so $A_1z+a_2t+b$ is a column vector. Setting $a_2=0$ gives an autonomous field, one that does not depend on time.

Explore

Network width n Horizon T

Show

A damped pendulum,

\dot z_1=z_2,\ \dot z_2=-\sin z_1-0.3z_2

. Drag the ring

K

of starting points. true trajectories trajectories of the fitted network

F_\theta(z)=W\tanh(Az+b)

, with

n

hidden units fitted to the true field on the grey box.

Worked example

Why the theorem needs a compact set and a finite horizon.

Solve $\dot z=z^2$ , $z(0)=z_0>0$ . What happens?

Separate: $dz/z^2=dt$ , so $-1/z=t-1/z_0$ and $z(t)=\dfrac{z_0}{1-z_0t}$ . It blows up at $t=1/z_0$ . On $[0,T]$ with $T\ge1/z_0$ there is no solution to approximate.
Take $K=[0.1,0.5]$ and $T=1$ . Does every solution exist on $[0,T]$ , and how large does it get?

Yes: blow-up happens at $1/z_0\ge2>1$ . The largest value is $z(1)=0.5/(1-0.5)=1$ from $z_0=0.5$ . Everything stays in the compact interval $[0.1,1]$ , which is where the network needs to be accurate.
Why can't one network work for all $z_0\in(0,\infty)$ and all $t$ ?

Solutions from large $z_0$ leave every bounded set almost immediately, and a network with finitely many units is approximately linear far from the data (or bounded, with tanh), so it can't track $z^2$ everywhere. Restricting to a compact $K$ and finite $T$ keeps all the trajectories inside one bounded region. In the figure, drag $K$ far out or raise $T$ and watch the error grow.

Check your understanding

The flow map $\varphi_T$ for a system in $\R^d$ is a function from:

$\R^d$ to curves in $\R^d$ .That is what you get when $T$ is allowed to vary.
$\R^d$ to $\R^d$ : starting point to position at time $T$ .
$\R$ to $\R^d$ .That describes one trajectory $t\mapsto z(t)$ , with $z_0$ fixed.

Fix

T

and vary

z_0

: a map from vectors to vectors. Fix

z_0

and vary

t

: a curve. Asking “what kind of object is this?” sorts out most confusion in a paper.

Why train one network for all starting points rather than one network per trajectory?

One network approximates the law of motion $F$ , which predicts trajectories from starting points never observed.
It is cheaper to train.Maybe, but the real reason is what you learn.

Three separate fits would memorize three curves and say nothing about a fourth starting point. A single vector field is a model of the system.

In the figure, the trajectories of the fitted field drift away from the true ones as $T$ grows, even with a good fit. Why?

Small errors in the velocity accumulate along the trajectory and can grow over time.
The numerical solver loses accuracy after a while.The same solver is used for both fields.

By Grönwall's inequality, if

|F-F_\theta|\le\delta

and

F

is

L

-Lipschitz, trajectories differ by at most

\delta\,(e^{LT}-1)/L

: the bound grows exponentially in

T

. That is why approximation results fix a finite horizon.

Practice

A linear flowFlow maps compose

For $\dot z=Az$ with $A\in\R^{d\times d}$ constant, find $\varphi_t$ and show $\varphi_{s+t}=\varphi_s\circ\varphi_t$ .

Solution

$z(t)=e^{tA}z_0$ , so $\varphi_t$ is multiplication by the matrix $e^{tA}$ , and $e^{(s+t)A}=e^{sA}e^{tA}$ . Flowing for $t$ and then for $s$ is flowing for $s+t$ . This holds for any autonomous field; with a time-dependent $F$ the maps depend on the start time too.

Grönwall in one lineHow errors grow

Let $\dot z=F(z)$ and $\dot x=G(x)$ with $x(0)=z(0)$ , $|F-G|\le\delta$ everywhere and $F$ $L$ -Lipschitz. Show $|z(t)-x(t)|\le\frac\delta L(e^{Lt}-1)$ .

Solution

$e(t)=|z(t)-x(t)|$ satisfies $e(t)\le\int_0^t|F(z)-F(x)|+|F(x)-G(x)|\,ds\le\int_0^t(Le(s)+\delta)\,ds$ . Grönwall's lemma gives $e(t)\le\delta t$ when $L=0$ and $\frac\delta L(e^{Lt}-1)$ in general. So a velocity error of $\delta$ costs at most $\delta\,T\,e^{LT}$ on $[0,T]$ , which is why controlling the sup error of the field controls the sup error of the flow on a finite horizon.

InteractiveGoodfellow–Bengio–Courville, Ch. 8: Optimization for Training Deep Models

Gradient descent, momentum and mini-batches

Fitting a network means minimizing a loss over its weights. Gradient descent does it by walking downhill: from wherever you are, take a small step in the direction in which the loss decreases fastest.

The idea

Why the negative gradient? The rate of change of $L$ in a unit direction $v$ is the directional derivative $\nabla L\cdot v=\|\nabla L\|\,\|v\|\cos\phi$ , where $\phi$ is the angle between them. It is most negative when $\cos\phi=-1$ , that is, when $v$ points opposite to $\nabla L$ . So

w_{t+1}=w_t-\alpha\,\nabla L(w_t).

The gradient says which way to go, not how far: the learning rate $\alpha$ is a choice. In machine learning $L(\theta)=\frac1N\sum_{i=1}^N\ell\big(f_\theta(x_i),y_i\big)$ is the average loss over the data, and computing its gradient through the network's layers is backpropagation.

Quadratics, one direction at a time. For

L(w)=\tfrac12(w-w^\star)^\top H(w-w^\star)

with

H

symmetric positive definite,

\nabla L=H(w-w^\star)

, so

w_{t+1}-w^\star=(I-\alpha H)(w_t-w^\star).

Along an eigenvector of

H

with eigenvalue

\lambda

, the error is multiplied by

1-\alpha\lambda

each step. Steep directions (large

\lambda

) need

\alpha<2/\lambda_{\max}

or the iterates blow up; flat directions (small

\lambda

) then crawl at rate

1-\alpha\lambda_{\min}

. The ratio

\kappa=\lambda_{\max}/\lambda_{\min}

sets the speed: about

\kappa

steps per factor of

e

.

Momentum (heavy ball). Keep part of the previous step:

w_{t+1}=w_t-\alpha\nabla L(w_t)+\beta\,(w_t-w_{t-1}).

In a flat valley the steps all point the same way and add up; across a narrow valley they alternate and cancel. Well tuned, it needs about

\sqrt\kappa

steps instead of

\kappa

.

Mini-batches. With millions of data points, use the average gradient over a random batch

B_t

of

|B|

of them:

g_t=\frac1{|B|}\sum_{i\in B_t}\nabla\ell_i

. It is unbiased,

\mathbb E[g_t]=\nabla L

, but noisy, with variance shrinking like

1/|B|

. The noise is the price of cheap steps, and it keeps the iterates jittering near the minimum unless

\alpha

is decreased.

Explore

Steepness ratio κ Learning rate α Momentum β Batch noise

Contours of

L(w)=\tfrac12(w_1^2+\kappa w_2^2)

; the minimum is at the centre. Drag the start. the first 60 steps.

Worked example

$L(w)=\tfrac12(w_1^2+10w_2^2)$ , so $H=\mathrm{diag}(1,10)$ . Start at $w_0=(1,1)$ .

With $\alpha=0.1$ : the per-step factors in each coordinate, and $w_t$ .

$1-0.1\cdot1=0.9$ and $1-0.1\cdot10=0$ . So $w_1=(0.9,0)$ : the steep coordinate is solved in one step, and $w_t=(0.9^t,0)$ after that. Reaching $|w_1|<10^{-3}$ takes $t\ge66$ steps.
With $\alpha=0.25$ ?

Factors $0.75$ and $1-2.5=-1.5$ . The second coordinate flips sign and grows by 1.5 each step: divergence, since $\alpha>2/\lambda_{\max}=0.2$ .
What is the best fixed $\alpha$ , and how fast is it?

Balance the two factors: $1-\alpha=-(1-10\alpha)$ , so $\alpha=2/11$ and both have size $9/11=\frac{\kappa-1}{\kappa+1}\approx0.82$ . About 35 steps for three digits. Momentum improves this to roughly $\frac{\sqrt\kappa-1}{\sqrt\kappa+1}\approx0.52$ per step.

Check your understanding

Why step against the gradient rather than along it?

$\nabla L\cdot v=\|\nabla L\|\|v\|\cos\phi$ is smallest at $\phi=\pi$ : the loss falls fastest opposite to the gradient.
Along the gradient the step would be too long.Length is set by $\alpha$ ; direction is the issue.

The gradient points uphill, in the direction of fastest increase. Gradient ascent (

+\alpha\nabla L

) maximizes, as in policy-gradient methods that maximize expected return.

On an ill-conditioned quadratic ( $\kappa=100$ ) plain gradient descent is slow because:

$\alpha$ must be below $2/\lambda_{\max}$ for the steep direction, which makes progress in the flat direction tiny.
The gradient points away from the minimum.It points downhill, just not at the minimum: it zigzags across the valley.

One learning rate serves all directions. Momentum, preconditioning and adaptive methods (Adam) are different ways around this.

Doubling the batch size $|B|$ changes the mini-batch gradient's standard deviation by a factor of:

$1/2$ That is the variance, not the standard deviation.
$1/\sqrt2$
$1$ : the noise comes from the data, not the batch.Averaging more independent samples reduces it.

The average of

|B|

independent gradients has variance

\sigma^2/|B|

. Twice the work per step buys only

\sqrt2

less noise, which is one reason small batches are popular.

Practice

The gradient of a quadraticThree lines

For $L(w)=\tfrac12w^\top Hw-q^\top w$ with $H$ symmetric, show $\nabla L=Hw-q$ , find the minimizer, and rewrite gradient descent as $w_{t+1}-w^\star=(I-\alpha H)(w_t-w^\star)$ .

Solution

$\partial_k\tfrac12\sum_{ij}H_{ij}w_iw_j=\tfrac12\sum_j(H_{kj}+H_{jk})w_j=(Hw)_k$ by symmetry, and $\partial_kq^\top w=q_k$ . The minimizer solves $Hw^\star=q$ . Then $w_{t+1}=w_t-\alpha(Hw_t-q)=w_t-\alpha H(w_t-w^\star)$ ; subtract $w^\star$ .

Momentum on a lineThe characteristic equation

Apply heavy ball to $L(w)=\tfrac\lambda2w^2$ . Show that $e_t=w_t$ satisfies $e_{t+1}=(1+\beta-\alpha\lambda)e_t-\beta e_{t-1}$ , and find when it converges for $\beta=0.5$ , $\lambda=1$ .

Solution

Substitute $\nabla L=\lambda w$ . The roots of $z^2-(1+\beta-\alpha\lambda)z+\beta=0$ have product $\beta=0.5$ ; the iteration converges iff both have modulus below 1, which (Jury's test) means $|1+\beta-\alpha\lambda|<1+\beta$ : $0<\alpha<3$ . Momentum enlarges the stable range from $\alpha<2$ to $\alpha<2(1+\beta)$ .

InteractiveGoodfellow–Bengio–Courville, §8.1.3 and §8.3.1

Mini-batches and stochastic gradients

The loss is an average over the whole dataset, and so is its gradient. With millions of data points, computing it exactly for every step is too expensive. Instead, pick a small random batch, average its gradients, and step. The step is now random, and that randomness is the whole story of this module.

The idea

With data $(x_i,y_i)$ , $i=1,\dots,N$ , and a per-example loss $\ell_i(\theta)=\ell(f_\theta(x_i),y_i)$ ,

L(\theta)=\frac1N\sum_{i=1}^N\ell_i(\theta),\qquad\nabla L(\theta)=\frac1N\sum_{i=1}^N\nabla\ell_i(\theta).

One full gradient step costs $N$ per-example gradients.

Mini-batch SGD. At step

t

, draw a batch

B_t

of

|B|

indices uniformly at random and step along

g_t=\frac1{|B|}\sum_{i\in B_t}\nabla\ell_i(\theta_t),\qquad\theta_{t+1}=\theta_t-\alpha\,g_t.

Each step costs

|B|

gradients instead of

N

.

Why it works, and what it costs. Where does the randomness come from? Only from which examples were picked; a different batch gives a different step. Averaged over the possible batches,

\mathbb E[g_t]=\nabla L(\theta_t)

: the step is right on average (unbiased). Its variance is

\frac{\sigma^2}{|B|}

when sampling with replacement, where

\sigma^2

is the variance of a single example's gradient (without replacement, multiply by

\frac{N-|B|}{N-1}

). So the noise shrinks like

1/\sqrt{|B|}

while the cost grows like

|B|

.

The consequence: with a fixed step size $\alpha$ , SGD does not converge to the minimizer. It reaches the neighbourhood quickly and then jitters around it, at a scale set by $\alpha$ and the noise. Decreasing $\alpha$ over time (a schedule) or increasing $|B|$ shrinks the jitter.

Explore

Batch size |B| Step size α

100 data points. this step's batch and the current line

y=ax+c

.

Loss over

(a,c)

. full gradient descent mini-batch SGD minimizer.

Worked example

Fit $f_\theta(x)=\theta x$ to four points $x=(1,2,3,4)$ , $y=(2,3,7,8)$ with $\ell_i=\tfrac12(\theta x_i-y_i)^2$ . Take $\theta=1$ .

Compute each $\nabla\ell_i=(\theta x_i-y_i)x_i$ and the full gradient.

Residuals $\theta x_i-y_i=(-1,-1,-4,-4)$ ; gradients $(-1,-2,-12,-16)$ . Full gradient $-31/4=-7.75$ .
List the mini-batch gradients for all six batches of size 2.

$\{1,2\}:-1.5$ , $\{1,3\}:-6.5$ , $\{1,4\}:-8.5$ , $\{2,3\}:-7$ , $\{2,4\}:-9$ , $\{3,4\}:-14$ . They range from $-1.5$ to $-14$ .
Average them.

$-46.5/6=-7.75$ , the full gradient exactly. Unbiased, but any single step can be five times too short or nearly twice too long.

Check your understanding

The randomness in a mini-batch gradient comes from:

Noise in the labels $y_i$ .The data are fixed once collected; full-batch descent on them is deterministic.
Which examples happened to be drawn into the batch.
Rounding errors.Those are tiny and not the point.

Whenever a quantity is random, ask what it depends on. Here the only random input is the choice of indices, made to avoid computing all

N

gradients.

You quadruple the batch size. Per step, the gradient's standard deviation:

halves, and each step costs four times as much.
quarters, at four times the cost.That's the variance.

Standard deviation

\sigma/\sqrt{|B|}

. Diminishing returns: past some batch size, more examples per step buy little extra accuracy per unit of compute.

In the figure, with $|B|=1$ and a large step, the SGD path ends in a cloud around the minimizer rather than at it. To land closer you could:

Decrease the step size over time, or use bigger batches.
Run more steps with the same settings.The cloud has a fixed size for fixed α and |B|: more steps just wander inside it.

With step

\alpha

and noise variance

\sigma^2/|B|

, the stationary jitter has size of order

\sqrt{\alpha\sigma^2/|B|}

. Robbins–Monro step sizes (

\sum\alpha_t=\infty

,

\sum\alpha_t^2<\infty

) drive it to zero.

Practice

UnbiasednessProve it

Let $B$ be $|B|$ indices drawn independently and uniformly from $\{1,\dots,N\}$ . Show $\mathbb E[g]=\nabla L$ and $\mathrm{Var}(g)=\sigma^2/|B|$ (coordinatewise), where $\sigma^2=\frac1N\sum_i|\nabla\ell_i-\nabla L|^2$ is computed per coordinate.

Solution

Each draw $I$ has $\mathbb E[\nabla\ell_I]=\frac1N\sum_i\nabla\ell_i=\nabla L$ and variance $\sigma^2$ . $g$ is the average of $|B|$ independent copies: same mean, variance divided by $|B|$ .

EpochsCounting passes

With $N=50{,}000$ examples and $|B|=100$ , how many SGD steps make one pass over the data (an epoch)? How many full-gradient steps cost the same?

Solution

500 SGD steps per epoch; the same compute buys one full-gradient step. Early in training, 500 noisy steps make far more progress than one exact one, which is why mini-batches are standard.

SimulationGoodfellow–Bengio–Courville, §8.2 and §8.3.2

Quadratics, spectra and power laws

On a quadratic loss, gradient descent splits into independent one-dimensional problems, one for each eigenvector of the Hessian. Each shrinks at its own rate, and the slow ones, with small eigenvalues, set the pace. When the eigenvalues follow a power law, so does the loss: $L_t\approx t^{-\zeta}$ .

The idea

The model. $L(w)=\tfrac12w^\top Hw-q^\top w$ . Check the objects: $w,q$ are column vectors, $H$ is a symmetric matrix, so $w^\top Hw$ and $q^\top w$ are numbers and $L$ is a number. The gradient with respect to $w$ (the variable we choose; $H$ and $q$ are given) is $Hw-q$ , so the minimizer has $Hw^\star=q$ .

Where $H$ comes from. For a linear model with squared loss, $\ell=\tfrac12(\langle w,x\rangle-y)^2$ averaged over data, the Hessian is $H=\mathbb E[xx^\top]$ . Note the order: $x^\top x$ is a number (inner product), $xx^\top$ is a matrix (outer product). $H$ is the covariance of the inputs, so its eigenvalues measure how much the data vary in each direction.

Spectral decomposition. Write

H=\sum_k\lambda_ku_ku_k^\top

with orthonormal eigenvectors

u_k

. “Decomposing” the error vector

e_t=w_t-w^\star

means writing it in this basis,

e_t=\sum_kc_k(t)u_k

. Since

e_{t+1}=(I-\alpha H)e_t

,

c_k(t)=(1-\alpha\lambda_k)^t\,c_k(0),\qquad L(w_t)-L(w^\star)=\tfrac12\sum_k\lambda_k(1-\alpha\lambda_k)^{2t}c_k(0)^2.

Directions with

\alpha\lambda_k

near 1 are fixed in a step or two; directions with tiny

\lambda_k

barely move.

Power laws. Suppose the error mass in small eigen-directions follows a power law:

\sum_{k:\lambda_k\le\lambda}\lambda_kc_k(0)^2\approx\lambda^{\zeta}

for small

\lambda

, with an exponent

\zeta>0

that is a property of the problem, not something you choose. Since

(1-\alpha\lambda)^{2t}\approx e^{-2\alpha\lambda t}

is near 0 when

\lambda\gg1/(\alpha t)

and near 1 when

\lambda\ll1/(\alpha t)

, after

t

steps only the directions with

\lambda\lesssim1/(\alpha t)

remain:

L(w_t)-L(w^\star)\approx\tfrac12\sum_{k:\lambda_k\lesssim1/(\alpha t)}\lambda_kc_k(0)^2\approx(\alpha t)^{-\zeta}.

Larger

\zeta

means less error hides in the flat directions, and faster learning. Power-law spectra like this are what is observed for real data and kernels, which is one reason training losses often fall like a power of the step count.

Momentum and acceleration. With constant $\alpha$ and $\beta$ , heavy ball behaves like gradient descent with step $\alpha/(1-\beta)$ on the slow directions: a constant-factor gain, the same exponent $\zeta$ . Schedules that vary $\alpha_t,\beta_t$ with $t$ (built from orthogonal polynomials such as Chebyshev or Jacobi polynomials) can roughly double the exponent to $t^{-2\zeta}$ on deterministic problems. With mini-batch noise these aggressive schedules can amplify the noise and diverge, which is an active research question.

Polynomials and contours. After $t$ steps of any such method the error is $e_t=p_t(H)e_0$ for a polynomial $p_t$ of degree $t$ with $p_t(0)=1$ (for gradient descent, $p_t(\lambda)=(1-\alpha\lambda)^t$ ). So the question “how fast can a method be?” becomes “how small can such a polynomial be on the spectrum?”. Functions of a matrix can be written with Cauchy's integral formula: for a closed curve $\Gamma$ in the complex plane enclosing all the eigenvalues,

f(H)=\frac1{2\pi i}\oint_\Gamma f(z)\,(zI-H)^{-1}\,dz,

which for diagonalizable

H

is the scalar formula

f(\lambda)=\frac1{2\pi i}\oint_\Gamma\frac{f(z)}{z-\lambda}dz

applied to each eigenvalue. Some recent work designs optimizers by choosing the contour

\Gamma

instead of the polynomial, and complex numbers are natural there: the per-direction iteration of a momentum method is a

2\times2

matrix whose eigenvalues are complex when it oscillates.

Explore

Spectrum decay ν (λₖ = k^−ν) Source exponent ζ Momentum β

A quadratic with 20,000 eigen-directions,

\lambda_k=k^{-\nu}

, step

\alpha=1

. Log–log axes. gradient descent heavy ball (constant β) slope

-\zeta

for reference.

Worked example

$H=\mathrm{diag}(1,\,0.01)$ , $w^\star=0$ , $w_0=(1,1)$ , step $\alpha=1$ .

Per-step factors and $w_t$ .

$1-\alpha\lambda_1=0$ and $1-\alpha\lambda_2=0.99$ : $w_1=(0,0.99)$ , then $w_t=(0,0.99^t)$ . The first direction is solved immediately.
The loss, and how many steps to reduce it by a factor of 1000.

$L_t=\tfrac12\cdot0.01\cdot0.99^{2t}$ . Need $0.99^{2t}\le10^{-3}$ : $t\ge\frac{\ln1000}{2\ln(1/0.99)}\approx344$ . All of that time is spent on the flat direction.
Could a larger step help?

Only up to $\alpha<2/\lambda_{\max}=2$ , where the steep factor is $|1-\alpha|$ . At $\alpha=1.98$ the flat factor is $0.9802$ , about 175 steps, but the steep direction now oscillates with factor $0.98$ . The best fixed step balances the two at $\alpha=2/1.01$ , giving $\frac{\kappa-1}{\kappa+1}\approx0.980$ for $\kappa=100$ . With thousands of directions spread over many scales, there is no balancing point, and the power law takes over.

Check your understanding

For a column vector $x\in\R^d$ , the object $xx^\top$ is:

The number $|x|^2$ .That is $x^\top x$ .
A $d\times d$ matrix of rank one.

A column times a row is a matrix (outer product); a row times a column is a number (inner product). Averaging the rank-one matrices

x_ix_i^\top

over the data gives the Hessian of the squared loss for a linear model.

“Spectral decomposition of the error vector $e_t$ ” means:

Computing the eigenvalues of $e_t$ .A vector has no eigenvalues; the spectrum is $H$ 's.
Writing $e_t$ as a combination of the eigenvectors of $H$ .

The spectrum belongs to the matrix. Decomposing a vector along it gives coordinates

c_k

, each of which gradient descent multiplies by its own factor

1-\alpha\lambda_k

.

In the figure, raising the momentum $\beta$ from 0 to 0.9 shifts the heavy-ball curve but leaves its slope at about $-\zeta$ . Why?

With constant $\alpha,\beta$ , slow directions behave like gradient descent with step $\alpha/(1-\beta)$ : a constant-factor speed-up, same exponent.
Momentum has no effect on quadratics.It moves the curve by a factor of about ten.

On a log-log plot, multiplying time by a constant shifts the curve sideways without changing its slope. Changing the exponent needs step sizes and momenta that change with

t

.

Practice

From ν and κ to ζBuild a power law

Let $\lambda_k=k^{-\nu}$ and $\lambda_kc_k(0)^2=k^{-1-\kappa}$ . Show that $\sum_{\lambda_k\le\lambda}\lambda_kc_k(0)^2\approx\lambda^{\kappa/\nu}$ , so $\zeta=\kappa/\nu$ .

Solution

$\lambda_k\le\lambda$ iff $k\ge\lambda^{-1/\nu}$ . The tail $\sum_{k\ge k_0}k^{-1-\kappa}\approx k_0^{-\kappa}/\kappa$ , and $k_0^{-\kappa}=\lambda^{\kappa/\nu}$ . This is how the figure builds its problem: for given $\nu$ and $\zeta$ it sets $\kappa=\zeta\nu$ .

Momentum per directionComplex eigenvalues

For one eigen-direction, heavy ball is $c_{t+1}=(1+\beta-\alpha\lambda)c_t-\beta c_{t-1}$ . Write it as $\begin{pmatrix}c_{t+1}\\c_t\end{pmatrix}=M\begin{pmatrix}c_t\\c_{t-1}\end{pmatrix}$ , and find when $M$ 's eigenvalues are complex. What is their modulus then?

Solution

$M=\begin{pmatrix}1+\beta-\alpha\lambda&-\beta\\1&0\end{pmatrix}$ with characteristic polynomial $\mu^2-(1+\beta-\alpha\lambda)\mu+\beta$ . The roots are complex when $(1+\beta-\alpha\lambda)^2<4\beta$ , and then both have modulus $\sqrt\beta$ (their product is $\beta$ ). In that regime the coordinate spirals in at rate $\sqrt\beta$ per step, whatever $\lambda$ is: this is how momentum speeds up the flat directions.