Backpropagation
Last modified: 2026-10-25
Source
Built from the training and back-propagation half of Brian Hutchinson’s Fall 2025 lecture notes, which the neural networks chapter left as an outtake. The source states the notation — pre-activations \(\va\) and \(\v{b}\), post-activations \(\vhidval\) and \(\vy\) — and stops there, so we wrote the derivations, the cost accounting and the depth material here. The section on gradients with respect to the input comes from Logan Sizemore’s Spring 2025 lecture on backpropagation. We converted the sections on exploding gradients, choosing the hidden non-linearity and residual connections, and the remark on shattered gradients, from the “Training deep networks” material the neural networks chapter keeps as an outtake, and added their measurements and He initialization.
We computed every number in the solutions rather than quoting it. We also checked the three gradient formulas — single example, \(K\) layers, and minibatched — against central differences on random networks, and they agreed to about one part in \(10^{9}\).
1 Why this needs an algorithm
The neural networks chapter built a model and gradient descent built a way to fit one. Putting them together needs one thing neither chapter supplied: the gradient of a network’s loss with respect to its weights.
This looks like a detail. A gradient is a list of partial derivatives, each partial derivative is a limit, and we can approximate a limit by moving the parameter a little and dividing. We need no new mathematics at all.
The trouble is arithmetic rather than mathematics, so we first measure how large it is before looking for a way to avoid it.
Exercise 1 (Price the obvious way of getting a gradient)
Exercise 1 (Price the obvious way of getting a gradient) You already know one way to compute a derivative without doing any calculus: move a parameter a little, see how much the loss changes, and divide. For a single weight \(w\), the central difference estimate is
\[\frac{\partial \mathcal{L}}{\partial w} \approx \frac{\mathcal{L}(w + \smallconst) - \mathcal{L}(w - \smallconst)}{2\smallconst}.\]
This needs no derivation, works for any loss, and is easy to get right. Before building anything more clever, find out what it would cost.
Take the network from the neural networks chapter: \(\ninputvars = 784\) inputs, \(L = 100\) hidden units, \(C = 10\) outputs, so \(W\) is \(785 \times 100\) and \(U\) is \(101 \times 10\).
How many parameters does the network have, and how many forward passes does the central-difference estimate need to produce the whole gradient?
Solution 1. \(W\) holds \(785 \times 100 = 78{,}500\) parameters and \(U\) holds \(101 \times 10 = 1{,}010\), so \(P = 79{,}510\). Almost \(99\%\) of them are in \(W\).
Each parameter needs two losses, at \(w + \smallconst\) and \(w - \smallconst\), and each loss is a full forward pass. That is \(2P = 159{,}020\) forward passes for one gradient. (A one-sided difference halves this to \(P + 1\) and is considerably less accurate; Solution 4 says why.)
Exercise 2 (Price the obvious way of getting a gradient (part 2))
Exercise 2 (Price the obvious way of getting a gradient (part 2)) Use the setup of Exercise 1.
Count the multiply-accumulate operations in one forward pass. (Count only the matrix products; the activation functions are one operation per unit and do not change the order of magnitude.) Multiply through: what does one numerical gradient cost, for one training example?
Solution 2. The forward pass is two matrix products: \(\tp{W} \vx\) costs \(785 \times 100 = 78{,}500\) multiply-accumulates and \(\tp{U} [1; \vhidval]\) costs \(101 \times 10 = 1{,}010\), for \(79{,}510\) in total — one per parameter, which is not a coincidence. Every weight is used exactly once in a forward pass.
So one numerical gradient, for one training example, costs
\[159{,}020 \times 79{,}510 \approx 1.26 \times 10^{10} \text{ multiply-accumulates}.\]
The cost is quadratic in the number of parameters, which is the real problem: doubling the width of the hidden layer multiplies the cost by four.
Exercise 3 (Price the obvious way of getting a gradient (part 3))
Exercise 3 (Price the obvious way of getting a gradient (part 3)) Use the setup of Exercise 1.
MNIST has \(60{,}000\) training images. What does one epoch of full-batch gradient descent cost this way? Compare it with a method that computes the whole gradient in the time of about three forward passes.
Solution 3. Multiply by \(60{,}000\):
\[60{,}000 \times 1.26 \times 10^{10} \approx 7.6 \times 10^{14}\]
multiply-accumulates for a single step. A method costing three forward passes needs \(60{,}000 \times 3 \times 79{,}510 \approx 1.4 \times 10^{10}\), a factor of about \(53{,}000\) less. Put differently: one step of the numerical method costs what fifty thousand steps of the other one costs. Training is thousands of steps, so this is not a constant factor to optimize later — it decides whether we can fit the model at all.
Exercise 4 (Price the obvious way of getting a gradient (part 4))
Exercise 4 (Price the obvious way of getting a gradient (part 4)) Use the setup of Exercise 1.
The method has an accuracy problem as well as a speed problem. Taking \(\smallconst\) too large leaves a truncation error, and taking it too small leaves a rounding error, because the two losses being subtracted agree in their leading digits. Which way does each error move as \(\smallconst\) shrinks, and roughly where is the best compromise in double precision?
Solution 4. Write \(\mathcal{L}\) as a smooth function of the one weight. Taylor expansion gives
\[\frac{\mathcal{L}(w+\smallconst) - \mathcal{L}(w-\smallconst)}{2\smallconst} = \mathcal{L}'(w) + \frac{\smallconst^2}{6}\mathcal{L}'''(w) + \cdots,\]
so the truncation error falls like \(\smallconst^2\). Against that, the computer evaluates both losses to a relative accuracy of about \(\machineps_{\text{mach}} \approx 2.2 \times 10^{-16}\), and the subtraction cancels their leading digits, so the rounding error in the quotient grows like \(\machineps_{\text{mach}} / \smallconst\).
The total is smallest where the two are comparable, at \(\smallconst \sim \machineps_{\text{mach}}^{1/3} \approx 6 \times 10^{-6}\), leaving an error of about \(\machineps_{\text{mach}}^{2/3} \approx 4 \times 10^{-11}\). A one-sided difference has truncation error \(O(\smallconst)\) instead, so its best compromise is at \(\smallconst \sim \sqrt{\machineps_{\text{mach}}} \approx 1.5 \times 10^{-8}\), where the truncation and rounding terms are both that size, leaving an error worse by a factor of \(\sqrt{\machineps_{\text{mach}}} / \machineps_{\text{mach}}^{2/3} \approx 400\) for half the work.
Neither is catastrophic on its own. You pay \(10^{10}\) operations for an approximate gradient when an exact one is available for \(10^5\).
Exercise 5 (Price the obvious way of getting a gradient (part 5))
Exercise 5 (Price the obvious way of getting a gradient (part 5)) Use the setup of Exercise 1.
The gradient is a list of \(P\) numbers, so no method can produce it in less than \(P\) operations. Say what the best imaginable algorithm would cost here, in units of forward passes, and what that target implies about how much work each parameter may receive.
Solution 5. The output has \(P\) numbers in it, so \(P\) operations is a lower bound, and one forward pass already costs \(P\) here. The best imaginable algorithm therefore runs in a small constant number of forward passes — and the constant has to be small enough that each parameter receives only \(O(1)\) work per training example.
That is exactly what back-propagation achieves, and the rest of this chapter builds it. Solution 2 is where to look for the reason: a forward pass touches each weight once, so a single pass already carries past each weight everything that is local to it. It carries nothing about the loss, which sits at the far end of the network. The next few exercises ask what has to come back the other way, and how little of it there turns out to be.
Watch
- Sanderson explains backpropagation with pictures and little algebra, from 3:07 [1, Ch. 3].
- Starmer explains the main ideas of backpropagation [2].
Watch
Hutchinson covers the idea of backpropagation (15 min) [3]. The video site requires a login, which enrolled students find in the syllabus.
Read
GeeksforGeeks explains backpropagation [4].
2 Naming the pieces
The obvious method fails because it treats the network as a function it can only evaluate, and evaluates it \(159{,}020\) separate times. Everything that follows comes from using the network’s structure: a network is a composition of simple steps, the derivative of each step is easy, and the chain rule assembles them.
So we give each step a name on both sides of it.
Exercise 6 (Give the intermediate quantities names)
Exercise 6 (Give the intermediate quantities names) The network chapter wrote the whole model in two lines. To differentiate it, we break the two lines into four, so that every place a function is applied has a name on each side of it:
\[\va = \tp{W} \vx, \qquad \vhidval = f(\va), \qquad \v{b} = \tp{U} \begin{bmatrix} 1 \\ \vhidval \end{bmatrix}, \qquad \vy = g(\v{b}).\]
\(\va\) and \(\v{b}\) are the pre-activations, \(\vhidval\) and \(\vy\) the post-activations.
Give the shape of each of the four vectors, for \(\ninputvars\) inputs, \(L\) hidden units and \(C\) outputs.
Solution 6. \(\vx \in \reals^{\ninputvars+1}\) (the \(+1\) is the constant input), \(\va \in \reals^{L}\), \(\vhidval \in \reals^{L}\), \(\v{b} \in \reals^{C}\), \(\vy \in \reals^{C}\). \(W\) is \((\ninputvars+1) \times L\) and \(U\) is \((L+1) \times C\); the bracket \([1; \vhidval]\) supplies \(U\)’s extra row.
Exercise 7 (Give the intermediate quantities names (part 2))
Exercise 7 (Give the intermediate quantities names (part 2)) Use the setup of Exercise 6.
The equation \(\vhidval = f(\tp{W} \vx)\) says the same thing as the first two equations together. Give a reason to keep \(\va\) as a separate name anyway, in terms of what the chain rule will need.
Solution 7. Because we compute the derivative with respect to the pre-activation once and use it twice.
Writing \(\backerr^a_j = \partial \mathcal{L} / \partial a_j\), every weight in column \(j\) of \(W\) has
\[\frac{\partial \mathcal{L}}{\partial W_{ij}} = \backerr^a_j \, x_i,\]
so the \(L\) numbers in \(\vbackerr^a\) give all \((\ninputvars+1)L\) derivatives in \(W\). We also pass the same \(\vbackerr^a\) further back in a deeper network. If you never name \(\va\), you have no name for the quantity the algorithm is actually propagating — which is why the method is named after it.
Exercise 8 (Give the intermediate quantities names (part 3))
Exercise 8 (Give the intermediate quantities names (part 3)) Use the setup of Exercise 6.
The backward pass will reuse quantities from the forward pass rather than recompute them. Before deriving it, predict which of \(\vx, \va, \vhidval, \v{b}, \vy\) have to survive until the backward pass and which you can discard as soon as you compute the next line, and say how much memory your answer costs for a minibatch of \(B\) examples. Exercise 11 and Exercise 18 will settle it, and predicting first gives you something to check them against.
Solution 8. Deriving Exercise 11 and Exercise 18 settles it. The backward pass needs:
- \(\vx\), to form \(\partial \mathcal{L} / \partial W\);
- \(\vhidval\), to form \(\partial \mathcal{L} / \partial U\);
- \(\vy\), to start the recursion;
- and \(f'(\va)\), which for the logistic sigmoid is recoverable from \(\vhidval\) alone, since \(f'(a) = \hidval(1-\hidval)\).
We can discard \(\v{b}\) the moment we compute \(\vy\), and \(\va\) too whenever \(f'\) is expressible in \(\vhidval\) — which is why implementations usually store the post-activations and not the pre-activations.
For a minibatch of \(B\) examples that is \(B(\ninputvars + L + C + 1)\) numbers, against \(P = (\ninputvars+1)L + (L+1)C\) for the parameters themselves. With \(\ninputvars = 784\), \(L = 100\), \(C = 10\) and \(B = 128\), the stored activations come to \(114{,}560\) numbers and the parameters to \(79{,}510\) — the activations are the larger of the two, by a factor of \(1.44\). Holding both at once therefore costs \(2.44\) times what the model alone does, and in a deep network it is worse, because we hold every layer’s activations until the backward pass reaches them.
Exercise 9 (Give the intermediate quantities names (part 4))
Exercise 9 (Give the intermediate quantities names (part 4)) Use the setup of Exercise 6.
Separate the quantities that depend on the training example from those that do not, and say what that separation implies about which of them a minibatch has \(B\) copies of.
Solution 9. \(W\) and \(U\) do not depend on the example; \(\vx, \va, \vhidval, \v{b}, \vy\) all do. A minibatch therefore carries \(B\) copies of each activation and one copy of each parameter.
That asymmetry is the one Exercise 70 turns into a number. Whether reading each weight once for \(B\) examples gives any real benefit depends on figures that exercise supplies, so note the asymmetry here and leave the conclusion until then.
Exercise 10 (Give the intermediate quantities names (part 5))
Exercise 10 (Give the intermediate quantities names (part 5)) Use the setup of Exercise 6.
Suppose \(f\) were the identity. Which of the four names would still be worth keeping, and what does your answer say about why the nonlinearity is what makes the bookkeeping necessary?
Solution 10. If \(f\) were the identity then \(\vhidval = \va\) and one of the two names is redundant — and more than that, the whole two-layer structure collapses, since \(\v{b} = \tp{U} [1; \tp{W} \vx]\) is affine in \(\vx\) and could be replaced by a single matrix. There would be nothing to propagate through.
The bookkeeping exists because \(f\) lies between two linear maps and does not commute with them. Every complication in this chapter — the elementwise factor \(f'(\va)\), the stored activations, the vanishing gradients of Exercise 34 — traces back to that one function, and so does every gain in representational power from the neural networks chapter.
3 The last layer
Differentiation runs in the opposite direction to the computation: the loss is at the end, so the derivative that is easiest to write is the one nearest to it.
We carry over one result rather than rederive it. The neural networks chapter showed that a softmax output with a cross-entropy loss — and equally a sigmoid output with a binary cross-entropy loss, and a linear output with the halved squared error \(\tfrac12\norm{ \vy - \vt}^2\) — all give the same derivative at the output pre-activation:
\[ \begin{aligned} \vbackerr^{\,b} &= \frac{\partial \mathcal{L}}{\partial \v{b}} \\ &= \vy - \vt \end{aligned} \]
Prediction minus target. Each of the three matched pairs produces it, which is why those pairs are the standard ones. It is the only fact about the loss that the rest of this chapter uses. The \(\tfrac12\) cancels the \(2\) that differentiating a square produces, which is why the squared-error convention includes it. Exercise 17 drops it and shows the \(2\) reappear.
Exercise 11 (Differentiate the last layer first)
Exercise 11 (Differentiate the last layer first) The neural networks chapter established the one derivative this chapter takes as given: with a softmax output and the cross-entropy loss, the derivative of the loss with respect to the output pre-activation is
\[ \begin{aligned} \vbackerr^{\,b} &= \frac{\partial \mathcal{L}}{\partial \v{b}} \\ &= \vy - \vt, \end{aligned} \]
where \(\vt\) is the one-hot target. Prediction minus target, with the softmax and the logarithm cancelling each other out. Everything else in this chapter is the chain rule applied to that.
Write \(b_c = \sum_{j} U_{jc} [1; \vhidval]_j\) and differentiate it with respect to a single entry \(U_{jc'}\). Be explicit about when the answer is zero.
Solution 11. \(b_c\) involves only column \(c\) of \(U\), so
\[\frac{\partial b_c}{\partial U_{jc'}} = \begin{cases} [1; \vhidval]_j & c' = c \\ 0 & c' \neq c. \end{cases}\]
The zero is the useful half: changing a weight in column \(3\) cannot affect output \(5\) at all.
Exercise 12 (Differentiate the last layer first (part 2))
Exercise 12 (Differentiate the last layer first (part 2)) Use the setup of Exercise 11.
Combine Exercise 11 with \(\vbackerr^{\,b}\) to get \(\partial \mathcal{L} / \partial U_{jc}\).
Solution 12. The chain rule sums over the outputs, and Solution 11 makes every term zero except one:
\[ \begin{aligned} \frac{\partial \mathcal{L}}{\partial U_{jc}} &= \sum_{c'} \backerr^{\,b}_{c'} \frac{\partial b_{c'}}{\partial U_{jc}} \\ &= \backerr^{\,b}_{c} \, [1; \vhidval]_j \\ &= (y_c - t_c)\, [1; \vhidval]_j. \end{aligned} \]
Read it as a product of two things the forward pass already has: how wrong output \(c\) was, and how active unit \(j\) was. The gradient moves a weight in proportion to both. That is the form of the learning rule Donald Hebb proposed in 1949 — strengthen a connection between two units that are active together — with the second factor replaced by an error rather than an activity, so the connection that changes is the one between an active input and a wrong output. The resemblance is structural, and we should not take it further. Hebb’s rule is local and unsupervised, and this one needs a loss.
Exercise 13 (Differentiate the last layer first (part 3))
Exercise 13 (Differentiate the last layer first (part 3)) Use the setup of Exercise 11.
Write the whole matrix \(\partial \mathcal{L} / \partial U\) in one expression, with no indices, and check that its shape matches \(U\).
Solution 13. That is an outer product:
\[\frac{\partial \mathcal{L}}{\partial U} = \begin{bmatrix} 1 \\ \vhidval \end{bmatrix} \pt{\vbackerr^{\,b}}.\]
The shapes work out: \((L+1) \times 1\) times \(1 \times C\) gives \((L+1) \times C\), which is \(U\)’s shape.
Exercise 14 (Differentiate the last layer first (part 4))
Exercise 14 (Differentiate the last layer first (part 4)) Use the setup of Exercise 11.
What is the rank of that matrix, for a single training example?
Solution 14. An outer product of two non-zero vectors has rank 1.
Exercise 15 (Differentiate the last layer first (part 5))
Exercise 15 (Differentiate the last layer first (part 5)) Use the setup of Exercise 11.
What does that rank say about how much a single example can change \(U\)?
Solution 15. So a single example can only push \(U\) along one direction in a \((L+1)C\)-dimensional space, no matter how large the step. Fitting \(U\) properly needs at least \(\min(L+1, C)\) independent examples before the accumulated gradient can even have full rank. This is one concrete reason the stochastic steps of gradient descent have to be numerous and small: each one is a rank-one correction, and the useful direction only emerges from their average.
Exercise 16 (Differentiate the last layer first (part 6))
Exercise 16 (Differentiate the last layer first (part 6)) Use the setup of Exercise 11.
Row \(0\) of \(U\) multiplies the constant \(1\). Write its gradient, and say in one sentence what quantity it is the gradient of.
Solution 16. Setting \(j = 0\), where \([1; \vhidval]_0 = 1\):
\[\frac{\partial \mathcal{L}}{\partial U_{0c}} = y_c - t_c.\]
That row is the output bias, and its gradient is the prediction error itself, unweighted by any activation. It is the one parameter that moves on every example regardless of what the hidden layer did.
Exercise 17 (Differentiate the last layer first (part 7))
Exercise 17 (Differentiate the last layer first (part 7)) Use the setup of Exercise 11.
Suppose the task were regression instead: \(g\) the identity and \(\mathcal{L} =\norm{ \vy - \vt}^2\). Recompute \(\vbackerr^{\,b}\) and say which parts of your answers to Exercise 12 and Exercise 13 change.
Solution 17. With \(g\) the identity, \(\vy = \v{b}\), and
\[ \begin{aligned} \vbackerr^{\,b} &= \frac{\partial}{\partial \v{b}}\norm{ \v{b} - \vt}^2 \\ &= 2(\vy - \vt). \end{aligned} \]
Parts 2 and 3 do not change at all beyond the factor of \(2\) — the outer product is the same expression with a different \(\vbackerr^{\,b}\) in it.
This is the first sign of the structure the whole chapter has: the layer does not know what loss it is serving. The single vector \(\vbackerr^{\,b}\) summarizes everything above the pre-activation \(\v{b}\), and the layer’s own derivative is the same formula whatever produced it. That makes the recursion in Exercise 29 possible, and lets a library offer losses and layers as separate, interchangeable pieces.
Watch
- Sanderson applies the chain rule to a network, from 0:38 [1, Ch. 4].
- Starmer explains the chain rule [5], optimizes the last layer [6] of a small network, and differentiates cross-entropy [7] and the softmax [8].
4 One layer back
\(U\) was the easy matrix: it sits directly beneath the loss. \(W\) is underneath \(U\) and underneath \(f\), and it holds almost all of the parameters.
Passing the derivative back through those two obstacles is the whole of back-propagation [9]. The first is a matrix, and passing back through it costs a sum. The second is an elementwise function, and passing back through it costs a multiplication.
Exercise 18 (Push the error back through the hidden layer)
Exercise 18 (Push the error back through the hidden layer) Exercise 11 stopped at \(U\). The other matrix, \(W\), is the one holding \(99\%\) of the parameters, and \(W\) is separated from the loss by two things: the multiplication by \(U\), and the activation function \(f\). This exercise handles the multiplication and the activation function, one at a time.
Throughout, \(\vbackerr^{\,b} = \vy - \vt\) is known from Exercise 11.
Unlike a weight in \(U\), the hidden value \(\hidvalof{j}\) affects every output. Write \(\partial \mathcal{L} / \partial \hidvalof{j}\) as a sum over \(c\), and simplify that sum using \(b_c = \sum_{j'} U_{j'c}[1; \vhidval]_{j'}\).
Solution 18. \(\hidvalof{j}\) reaches the loss through all \(C\) outputs, so the chain rule sums over them:
\[ \begin{aligned} \frac{\partial \mathcal{L}}{\partial \hidvalof{j}} &= \sumto{c}{C} \backerr^{\,b}_c \frac{\partial b_c}{\partial \hidvalof{j}} \\ &= \sumto{c}{C} \backerr^{\,b}_c \, U_{jc}. \end{aligned} \]
Exercise 19 (Push the error back through the hidden layer (part 2))
Exercise 19 (Push the error back through the hidden layer (part 2)) Use the setup of Exercise 18.
Write your answer to Exercise 18 as a matrix-vector product, \(\vbackerr^{\,\hidval} = \;?\). The matrix is not quite \(U\): say which part of \(U\) it leaves out and why.
Solution 19. That sum is row \(j\) of \(U\) against \(\vbackerr^{\,b}\), so stacking over \(j\),
\[\vbackerr^{\,\hidval} = U_{1:L,\,:} \; \vbackerr^{\,b},\]
where \(U_{1:L,\,:}\) is \(U\) with its row \(0\) deleted. That row multiplies the constant \(1\), not any \(\hidvalof{j}\), so it contributes to no \(\hidvalof{j}\)’s derivative. Keeping it would produce an \((L+1)\)-vector whose first entry is the derivative of the loss with respect to a constant, which is not a quantity that exists.
(The constant sits first throughout this chapter, in position \(0\), as the neural networks chapter sets it up and as linear regression pads its inputs, so row \(j\) of \(U\) multiplies \(\hidvalof{j}\) and row \(0\) is the bias.)
This is the step that gives the method its name: the same matrix \(U\) that carried activations forward carries errors backward, transposed.
Exercise 20 (Push the error back through the hidden layer (part 3))
Exercise 20 (Push the error back through the hidden layer (part 3)) Use the setup of Exercise 18.
Now pass back through \(f\), using \(\hidvalof{j} = f(a_j)\). Give \(\backerr^{\,a}_j\), and explain why this step produces no sum, where Exercise 18 did. Express the result in vector form with the elementwise product \(\odot\).
Solution 20. \(f\) is applied elementwise, so \(a_j\) affects \(\hidvalof{j}\) and nothing else:
\[\backerr^{\,a}_j = \backerr^{\,\hidval}_j \, f'(a_j), \qquad \vbackerr^{\,a} = \vbackerr^{\,\hidval} \odot f'(\va).\]
There is no sum because the Jacobian \(\partial \vhidval / \partial \va\) is diagonal. Solution 18 needed a sum because \(\partial \v{b} / \partial \vhidval = \pt{U_{1:L,\,:}}\) is a full \(C \times L\) matrix. Elementwise functions are cheap to differentiate for exactly this reason, which is why we apply activation functions coordinate by coordinate rather than as a general \(\reals^L \to \reals^L\) map.
Exercise 21 (Push the error back through the hidden layer (part 4))
Exercise 21 (Push the error back through the hidden layer (part 4)) Use the setup of Exercise 18.
Finally, \(\va = \tp{W} \vx\) has exactly the shape Exercise 11 already handled. Give \(\partial \mathcal{L}/\partial W_{ij}\) and then the whole matrix without indices.
Solution 21. \(a_j = \sum_i W_{ij} \predvar_i\), so \(\partial a_j / \partial W_{ij} = \predvar_i\) and
\[\frac{\partial \mathcal{L}}{\partial W_{ij}} = \backerr^{\,a}_j \, x_i, \qquad \frac{\partial \mathcal{L}}{\partial W} = \vx \pt{\vbackerr^{\,a}}.\]
The same outer-product shape as \(U\)’s gradient, with \(\vx\) in place of \([1; \vhidval]\) and \(\vbackerr^{\,a}\) in place of \(\vbackerr^{\,b}\). The input to a layer and the error at its output are all a layer’s gradient ever needs.
Exercise 22 (Push the error back through the hidden layer (part 5))
Exercise 22 (Push the error back through the hidden layer (part 5)) Use the setup of Exercise 18.
Lay the four results out in order and check that the shapes compose: \((C) \to (L) \to (L) \to (\ninputvars+1) \times L\).
Solution 22. In order:
\[\vbackerr^{\,b} = \vy - \vt \;\in \reals^C \;\longrightarrow\; \vbackerr^{\,\hidval} = U_{1:L,:}\vbackerr^{\,b} \;\in \reals^L \;\longrightarrow\; \vbackerr^{\,a} = \vbackerr^{\,\hidval} \odot f'(\va) \;\in \reals^L \;\longrightarrow\; \vx \pt{\vbackerr^{\,a}} \in \reals^{(\ninputvars+1) \times L}.\]
\(U_{1:L,:}\) is \(L \times C\) against a \(C\)-vector, giving an \(L\)-vector; the elementwise product preserves that; and \((\ninputvars+1) \times 1\) against \(1 \times L\) gives \(W\)’s shape.
Exercise 23 (Push the error back through the hidden layer (part 6))
Exercise 23 (Push the error back through the hidden layer (part 6)) Use the setup of Exercise 18.
For \(f = \sigmoid\), write \(f'(a_j)\) in terms of \(\hidvalof{j}\) alone, give its largest possible value, and say what happens to \(\backerr^{\,a}_j\) when unit \(j\) is strongly on or strongly off.
Solution 23. \(\sigmoid' (a) = \sigmoid(a)\paren{1 - \sigmoid(a)} = \hidvalof{j}(1 - \hidvalof{j})\), so the derivative is recoverable from the stored activation and \(\va\) need not be kept.
It is largest at \(\hidvalof{j} = 1/2\), where it equals \(1/4\), and it goes to zero at both ends. A unit that is strongly on (\(\hidvalof{j} \approx 1\)) or strongly off (\(\hidvalof{j} \approx 0\)) has \(\backerr^{\,a}_j \approx 0\), so none of the weights feeding it move, however wrong the network’s output was.
That is the saturation problem. The same \(\sigmoid' \leq 1/4\) made the logic gates of the neural networks chapter work so cleanly — a unit driven to \(\pm 10\) gives a clear answer precisely because it has stopped responding to its input — and Exercise 34 shows what it costs once such units are stacked.
Watch
- Hutchinson covers backprop for a single hidden layer (28 min) [3].
- Starmer carries the chain rule back a layer [10].
5 The algorithm
Collecting the results, for one training example and one hidden layer:
Backpropagation algorithm (one hidden layer)
Forward. \(\va = \tp{W} \vx\), \(\;\vhidval = f(\va)\), \(\;\v{b} = \tp{U} [1; \vhidval]\), \(\;\vy = g(\v{b})\).
Backward.
\(\vbackerr^{\,b} = \vy - \vt\)
\(\dfrac{\partial \mathcal{L}}{\partial U} = [1; \vhidval]\,\pt{\vbackerr^{\,b}}\)
\(\vbackerr^{\,a} = \paren{U_{1:L,:}\,\vbackerr^{\,b}} \odot f'(\va)\)
\(\dfrac{\partial \mathcal{L}}{\partial W} = \vx\,\pt{\vbackerr^{\,a}}\)
Each of the four lines is a matrix-vector product, an outer product, or an elementwise multiplication. They share a shape: every parameter gradient is the layer’s input times the error at the layer’s output. That is the pattern the next section turns into a loop.
Exercise 24 (Count what the backward pass actually costs)
Exercise 24 (Count what the backward pass actually costs) Exercise 1 measured the cost of the obvious method and set a target: a gradient in a small constant number of forward passes. Exercise 11 and Exercise 18 produced an algorithm. Now check the algorithm against the target.
Use the same network: \(\ninputvars = 784\), \(L = 100\), \(C = 10\), so \(W\) is \(785 \times 100\), \(U\) is \(101 \times 10\), and one forward pass is \(79{,}510\) multiply-accumulates.
Count the multiply-accumulates in each of the four backward steps separately: \(\vbackerr^{\,\hidval} = U_{1:L,:}\vbackerr^{\,b}\), the outer product for \(\partial \mathcal{L}/\partial U\), the elementwise product giving \(\vbackerr^{\,a}\), and the outer product for \(\partial \mathcal{L}/\partial W\).
Solution 24. Each step is a product whose cost is the number of entries it touches:
| step | shape | multiply-accumulates |
|---|---|---|
| \(\vbackerr^{\,\hidval} = U_{1:L,:}\vbackerr^{\,b}\) | \((100 \times 10)(10)\) | \(1{,}000\) |
| \(\partial \mathcal{L}/\partial U = [1; \vhidval]\pt{\vbackerr^{\,b}}\) | \((101)(10)\) | \(1{,}010\) |
| \(\vbackerr^{\,a} = \vbackerr^{\,\hidval} \odot f'(\va)\) | \((100)\) | \(100\) |
| \(\partial \mathcal{L}/\partial W = \vx \pt{\vbackerr^{\,a}}\) | \((785)(100)\) | \(78{,}500\) |
Exercise 25 (Count what the backward pass actually costs (part 2))
Exercise 25 (Count what the backward pass actually costs (part 2)) Use the setup of Exercise 24.
Total them, and give the ratio of backward cost to forward cost. Which single step dominates, and why is that the same step that dominated the forward pass?
Solution 25. The total is \(80{,}610\), against \(79{,}510\) forward, a ratio of \(1.014\). The backward pass costs essentially the same as the forward pass.
\(\partial \mathcal{L}/\partial W\) dominates both, at \(78{,}500\) of each total, and for the same reason in both directions: it has one entry per parameter in \(W\), and \(W\) holds \(98.7\%\) of the parameters. Forward or backward, the work is proportional to the number of weights, because each pass touches each weight once.
Exercise 26 (Count what the backward pass actually costs (part 3))
Exercise 26 (Count what the backward pass actually costs (part 3)) Use the setup of Exercise 24.
Give the cost of forward-plus-backward in units of forward passes, and compare it with the \(2P\) forward passes of Exercise 1. By what factor has the gradient got cheaper?
Solution 26. Forward plus backward is \(160{,}120\) multiply-accumulates, or \(2.01\) forward passes. Common estimates usually say “about three”, and that figure is a count of matrix products rather than an allowance for overheads: one product forward, two back, as Solution 28 works out. This network comes to two rather than three because layer \(1\) skips one of its two backward products. \(\vbackerr^{\,\predvar}\) would be the error on the input, and the input is data, not a parameter — we fit nothing below it, so there is nothing to propagate to. Layer \(1\) holds \(98.7\%\) of the weights, so skipping its input-error product removes almost exactly one forward pass’ worth of work: restoring it gives \(238{,}620\) multiply-accumulates, or \(3.0011\) forward passes. A deep network needs the full three, because only its first layer skips that product.
Exercise 1 needed \(2P = 159{,}020\) forward passes, or \(1.26 \times 10^{10}\) multiply-accumulates. The gradient is cheaper by a factor of
\[\frac{1.26 \times 10^{10}}{160{,}120} \approx 79{,}000.\]
Exercise 27 (Count what the backward pass actually costs (part 4))
Exercise 27 (Count what the backward pass actually costs (part 4)) Use the setup of Exercise 24.
The quadratic in Exercise 1 has become linear. Say where it went: what does back-propagation compute once that the numerical method recomputed for every parameter?
Solution 27. It computes \(\vbackerr^{\,b}\) and \(\vbackerr^{\,a}\) once, and every parameter’s derivative is a product of one of those numbers with one stored activation.
The numerical method does not share any work. To learn about weight \(W_{ij}\) it runs the entire network twice, and the \(79{,}509\) other weights take part in that computation without teaching us anything about them. Back-propagation does the shared part once — the \(L\) numbers in \(\vbackerr^{\,a}\) — and then spends exactly one multiplication per parameter.
That is the whole idea, stated as a rule: the cost is quadratic when the work per parameter is a function evaluation, and linear when the work per parameter is a multiplication.
Exercise 28 (Count what the backward pass actually costs (part 5))
Exercise 28 (Count what the backward pass actually costs (part 5)) Use the setup of Exercise 24.
None of these counts depended on \(f\), \(g\), or the loss. State the general rule for the cost of a backward pass through any layer whose forward pass is a matrix product, and say what this implies for a network of \(K\) layers.
Solution 28. A layer computing \(\vv = \tp{M} \vu\) costs \(\abs{ M}\) multiply-accumulates forward. Backward it costs \(\abs{ M}\) for the input error \(\vbackerr^{\,u} = M \vbackerr^{\,v}\) and another \(\abs{ M}\) for the parameter gradient \(\vu \pt{\vbackerr^{\,v}}\) — so a backward pass is about twice a forward pass per layer, and the whole of training is a small constant times the cost of prediction.
For \(K\) layers, both passes are a sum over layers of \(\abs{ M^{(k)}}\), so everything stays linear in the total parameter count however deep the network gets. Depth costs memory, as Exercise 8 showed, and it costs numerical trouble, as Exercise 34 shows. It does not cost arithmetic.
Watch
The CS231N lecture on backpropagation builds the algorithm on a computational graph, from 35:52 and then the backpropagation algorithm, from 40:36 [11], Lecture 4.
Watch
The CS229 lecture covers backpropagation and improving neural networks [12], Lecture 11.
6 Depth
Nothing in the derivation used the fact that there was one hidden layer. The backward step took a derivative with respect to one pre-activation and produced a derivative with respect to the pre-activation below it — the same kind of object, one layer down.
A procedure that turns an object into another object of the same kind is a recursion, and running it \(K\) times makes a network deep.
Exercise 29 (Do it again, \(K\) times)
Exercise 29 (Do it again, \(K\) times) Nothing in Exercise 18 used the fact that there was only one hidden layer. The derivation passed back through a matrix, then through an activation function, and arrived at a quantity of exactly the same kind it started with. That is a recursion, and we can now write it down.
Take \(K\) hidden layers, with widths \(L_0 = \ninputvars\), \(L_1, \ldots, L_K\), and \(L_{K+1} = C\). Write \(\vhidval^{(0)} = \vx\) and, for \(k = 1, \ldots, K+1\),
\[\va^{(k)} = \pt{W^{(k)}} \begin{bmatrix} 1 \\ \vhidval^{(k-1)} \end{bmatrix}, \qquad \vhidval^{(k)} = f\paren{\va^{(k)}},\]
with \(g\) in place of \(f\) at the last layer, so \(\vy = \vhidval^{(K+1)}\).
Give the shape of \(W^{(k)}\) and the total parameter count.
Exercise 30 (Do it again, \(K\) times (part 2))
Exercise 30 (Do it again, \(K\) times (part 2)) Use the setup of Exercise 29.
Write the recursion that takes \(\vbackerr^{(k+1)} = \partial \mathcal{L} / \partial \va^{(k+1)}\) to \(\vbackerr^{(k)}\), and say which parts of Exercise 18 each piece comes from.
Solution 30. Both steps of Exercise 18, in the same order:
\[\vbackerr^{(k)} = \underbrace{\paren{W^{(k+1)}_{1:L_k,\,:} \; \vbackerr^{(k+1)}}} _{\text{across the matrix, @sol-deep-recursion-part2}} \odot \underbrace{f'\paren{\va^{(k)}}} _{\text{across the activation, @sol-deep-recursion-part3}}.\]
The row deletion is the same one: \(W^{(k+1)}\)’s row \(0\) multiplies the constant \(1\), not any unit of layer \(k\).
The recursion starts at the top with \(\vbackerr^{(K+1)} = \vy - \vt\) from Exercise 11. That first value is special only because the loss enters there.
Exercise 31 (Do it again, \(K\) times (part 3))
Exercise 31 (Do it again, \(K\) times (part 3)) Use the setup of Exercise 29.
Give \(\partial \mathcal{L} / \partial W^{(k)}\) in terms of \(\vbackerr^{(k)}\) and the stored activations.
Solution 31. Exactly the outer product of Exercise 21, with the layer’s own input in front:
\[\frac{\partial \mathcal{L}}{\partial W^{(k)}} = \begin{bmatrix} 1 \\ \vhidval^{(k-1)} \end{bmatrix} \pt{\vbackerr^{(k)}}.\]
Every layer’s gradient has the same form: the input to the layer, times the error at its output.
Exercise 32 (Do it again, \(K\) times (part 4))
Exercise 32 (Do it again, \(K\) times (part 4)) Use the setup of Exercise 29.
Write the whole algorithm as two loops, forward then backward, and say what has to be kept between them.
Solution 32. Forward, storing as you go:
for \(k = 1 \ldots K+1\): compute \(\va^{(k)}\), then \(\vhidval^{(k)}\), and keep \(\vhidval^{(k-1)}\).
Backward, starting from \(\vbackerr^{(K+1)} = \vy - \vt\):
for \(k = K+1 \ldots 1\): record \(\partial \mathcal{L}/\partial W^{(k)} = [1; \vhidval^{(k-1)}]\pt{\vbackerr^{(k)}}\), then, if \(k > 1\), form \(\vbackerr^{(k-1)} = \paren{W^{(k)}_{1:L_{k-1},:}\vbackerr^{(k)}} \odot f'(\va^{(k-1)})\).
Between the two loops we keep every layer’s post-activation \(\vhidval^{(k)}\) — which is also enough to recover \(f'(\va^{(k)})\) when \(f\) is the logistic sigmoid, since \(\sigmoid'(a) = \hidval(1-\hidval)\).
The order inside the backward loop matters. The parameter gradient uses \(\vbackerr^{(k)}\), and only then do we use \(W^{(k)}\) to produce \(\vbackerr^{(k-1)}\). Using the old \(W^{(k)}\) for that second step is essential. An implementation that updated \(W^{(k)}\) before propagating through it would be differentiating a network that never existed.
Exercise 33 (Do it again, \(K\) times (part 5))
Exercise 33 (Do it again, \(K\) times (part 5)) Use the setup of Exercise 29.
There is an alternative that never runs backwards: start at the inputs and propagate derivatives with respect to one chosen parameter forward through the network, the way the values themselves travel. Work out what that costs for all \(P\) parameters, and use the answer to explain why the algorithm is called back-propagation.
Solution 33. Propagating forward, you would carry \(\partial \vhidval^{(k)} / \partial \param\) for one chosen parameter \(\param\) alongside the values, at roughly the cost of the forward pass itself — and you would have to do it again for each of the \(P\) parameters. That is \(O(P)\) forward passes — the same order as Exercise 1, which needed \(2P\) of them, and it is no simpler to write.
The asymmetry is the reason for the direction. A forward sweep computes the derivatives of many outputs with respect to one input; a backward sweep computes the derivatives of one output with respect to many inputs. Training has one output — the scalar loss — and \(P\) inputs — the parameters. That is the case where going backwards wins, and it wins by a factor of \(P\).
Both sweeps are standard and both are used: automatic-differentiation libraries call them forward mode and reverse mode, and choose between them by comparing the number of inputs with the number of outputs. Back-propagation is reverse mode, specialized to the one function form that machine learning cares about.
7 What depth costs
The arithmetic stays linear in the parameter count however deep the network is, which makes depth look free. It is not.
The recursion multiplies \(\vbackerr\) by a derivative of the activation function at every layer it passes, and the logistic sigmoid’s derivative is never larger than \(1/4\). Over ten layers that is a factor of \(10^{-6}\). This is a large part of why networks deeper than two or three layers resisted training from the 1980s until the late 2000s, when better initialization, unsupervised pre-training and then the rectified linear unit (ReLU) each addressed one part of it.
8 Measuring it
The bound of Section 7 is for the worst case. Whether it matters in practice is a question about real initialized networks, so we measure it on \(200\) sigmoid networks with \(10\) hidden layers of \(20\) units, \(20\) inputs and \(3\) classes, each weight drawn from \(\mathcal{N}(0, 1/(\text{fan-in} + 1))\), each fed one standard normal input against one uniformly chosen one-hot target under softmax and cross-entropy. We report the Frobenius norm of \(\partial \mathcal{L}/\partial W^{(k)}\), averaged over the \(200\) networks, together with the geometric mean of the ratios between consecutive layers.
Each version draws its own random numbers, so the three agree to about two significant figures rather than exactly. Replacing the sigmoid with \(\max(a, 0)\) gives the ReLU figures of Exercise 40.
import numpy as np
K, L, D, C, N = 10, 20, 20, 3, 200
rng = np.random.default_rng(0)
widths = [D] + [L] * K + [C]
norms = np.zeros(K + 1)
for _ in range(N):
W = [rng.normal(0, 1 / np.sqrt(widths[k] + 1), (widths[k] + 1, widths[k + 1]))
for k in range(K + 1)]
post = [np.concatenate(([1.0], rng.normal(0, 1, D)))]
pre = []
for k in range(K + 1):
a = W[k].T @ post[k]
pre.append(a)
if k < K:
post.append(np.concatenate(([1.0], 1 / (1 + np.exp(-a)))))
e = np.exp(pre[K] - pre[K].max())
t = np.zeros(C)
t[rng.integers(C)] = 1.0
d = e / e.sum() - t
for k in range(K, -1, -1):
norms[k] += np.linalg.norm(np.outer(post[k], d))
if k > 0:
h = post[k][1:]
d = (W[k][1:, :] @ d) * h * (1 - h)
norms /= N
r = norms[1:] / norms[:-1]
print(f"layer 1 {norms[0]:.2e}")layer 1 1.11e-06output 2.12ratio 5.2e-07geomean 4.25# nolint start: object_name_linter. K, L, D, C, N and W match the prose.
K <- 10
L <- 20
D <- 20
C <- 3
N <- 200
widths <- c(D, rep(L, K), C)
set.seed(0)
norms <- numeric(K + 1)
for (n in seq_len(N)) {
W <- lapply(seq_len(K + 1), function(k) {
matrix(rnorm((widths[k] + 1) * widths[k + 1], 0, 1 / sqrt(widths[k] + 1)),
widths[k] + 1, widths[k + 1])
})
post <- list(c(1, rnorm(D)))
pre <- list()
for (k in seq_len(K + 1)) {
a <- as.vector(t(W[[k]]) %*% post[[k]])
pre[[k]] <- a
if (k <= K) post[[k + 1]] <- c(1, 1 / (1 + exp(-a)))
}
e <- exp(pre[[K + 1]] - max(pre[[K + 1]]))
t_vec <- numeric(C)
t_vec[sample.int(C, 1)] <- 1
d <- e / sum(e) - t_vec
for (k in seq(K + 1, 1)) {
norms[k] <- norms[k] + sqrt(sum(outer(post[[k]], d)^2))
if (k > 1) {
h <- post[[k]][-1]
d <- as.vector(W[[k]][-1, , drop = FALSE] %*% d) * h * (1 - h)
}
}
}
norms <- norms / N
r <- norms[-1] / head(norms, -1)
cat(sprintf("layer 1 %.2e\n", norms[1]))layer 1 1.13e-06output 2.06ratio 5.5e-07geomean 4.23using Random, LinearAlgebra, Statistics, Printf
K, L, D, C, N = 10, 20, 20, 3, 200;
widths = [D; fill(L, K); C];
rng = MersenneTwister(0);
norms = zeros(K + 1);
for _ in 1:N
W = [randn(rng, widths[k] + 1, widths[k+1]) / sqrt(widths[k] + 1) for k in 1:K+1]
post = [[1.0; randn(rng, D)]]
pre = Vector{Vector{Float64}}()
for k in 1:K+1
a = W[k]' * post[k]
push!(pre, a)
k <= K && push!(post, [1.0; 1 ./ (1 .+ exp.(-a))])
end
e = exp.(pre[K+1] .- maximum(pre[K+1]))
t = zeros(C)
t[rand(rng, 1:C)] = 1.0
d = e ./ sum(e) .- t
for k in K+1:-1:1
norms[k] += norm(post[k] * d')
if k > 1
h = post[k][2:end]
d = (W[k][2:end, :] * d) .* h .* (1 .- h)
end
end
end
norms ./= N;
r = norms[2:end] ./ norms[1:end-1];
@printf("layer 1 %.2e\n", norms[1])layer 1 1.06e-06output 2.06ratio 5.2e-07geomean 4.25Show code
# The Python measurement above as a function of the activation and its
# derivative, so the exercises below can rerun it with another activation.
# Run with the sigmoid it reproduces the Python tab's norms (checked below).
def measure(act, dact, K=10, L=20, D=20, C=3, N=200, seed=0):
rng = np.random.default_rng(seed)
widths = [D] + [L] * K + [C]
norms = np.zeros(K + 1)
for _ in range(N):
W = [rng.normal(0, 1 / np.sqrt(widths[k] + 1),
(widths[k] + 1, widths[k + 1]))
for k in range(K + 1)]
post = [np.concatenate(([1.0], rng.normal(0, 1, D)))]
pre = []
for k in range(K + 1):
a = W[k].T @ post[k]
pre.append(a)
if k < K:
post.append(np.concatenate(([1.0], act(a))))
e = np.exp(pre[K] - pre[K].max())
t = np.zeros(C)
t[rng.integers(C)] = 1.0
d = e / e.sum() - t
for k in range(K, -1, -1):
norms[k] += np.linalg.norm(np.outer(post[k], d))
if k > 0:
h = post[k][1:]
d = (W[k][1:, :] @ d) * dact(h)
return norms / N
sigmoid_norms = measure(lambda a: 1 / (1 + np.exp(-a)), lambda h: h * (1 - h))Show code
# Figures quoted by the two exercises below and their solutions,
# from the Python measurements.
stopifnot(isTRUE(all.equal(py$sigmoid_norms, py$norms)))
sci <- function(v, digits) {
parts <- strsplit(formatC(v, format = "e", digits = digits), "e")[[1]]
sprintf("%s \\times 10^{%d}", parts[1], as.integer(parts[2]))
}
summarize_norms <- function(norms) {
k <- length(norms)
list(
first = norms[1],
last = norms[k],
ratio = norms[1] / norms[k],
geomean = exp(mean(log(norms[-1] / norms[-k])))
)
}
english_count <- function(n) {
c("one", "two", "three", "four", "five", "six", "seven", "eight", "nine")[n]
}
sig <- summarize_norms(py$norms)
relu <- summarize_norms(py$relu_norms)Show code
# HTML formats get an interactive plotly figure;
# the handout and other formats get the same line drawn with ggplot2.
# The norms are the Python measurement above,
# so the figure shows the numbers the prose quotes.
grad_norms <- data.frame(layer = seq_along(py$norms), norm = py$norms)
# Plain-number tick labels on the log axis, not scientific notation.
tick_values <- 10^seq(
floor(log10(min(grad_norms$norm))),
ceiling(log10(max(grad_norms$norm)))
)
tick_labels <- format(
tick_values,
scientific = FALSE,
trim = TRUE,
drop0trailing = TRUE
)
if (knitr::is_html_output()) {
library(plotly)
widget <- plot_ly(
grad_norms,
x = ~layer,
y = ~norm,
type = "scatter",
mode = "lines+markers",
line = list(color = "#1f77b4"),
marker = list(color = "#1f77b4"),
hovertemplate = "layer %{x}<br>norm %{y:.3g}<extra></extra>"
) |>
layout(
xaxis = list(title = "layer (1 = nearest the input)", dtick = 1),
yaxis = list(
title = "average gradient norm (log scale)",
type = "log",
tickvals = tick_values,
ticktext = tick_labels
),
margin = list(t = 10)
) |>
config(displaylogo = FALSE, responsive = TRUE)
# The alt text: a plotly widget has none of its own.
htmltools::div(
role = "img",
`aria-label` = paste(
"Interactive plot of average gradient norm against layer number",
"on a logarithmic axis.",
"The norm rises from the layer nearest the input to the output layer,",
"by several orders of magnitude."
),
widget
)
} else {
library(ggplot2)
ggplot(grad_norms, aes(x = layer, y = norm)) +
geom_line(colour = "#1f77b4") +
geom_point(colour = "#1f77b4") +
scale_x_continuous(breaks = grad_norms$layer) +
scale_y_log10(breaks = tick_values, labels = tick_labels) +
labs(
x = "layer (1 = nearest the input)",
y = "average gradient norm (log scale)"
)
}Exercise 34 (Find out what the recursion does to the size of a gradient)
Exercise 34 (Find out what the recursion does to the size of a gradient) Exercise 29 says depth costs no extra arithmetic. It says nothing about whether the numbers coming back are useful when they arrive.
Unroll the recursion from the output down to layer \(1\):
\[\vbackerr^{(1)} = \sb{ \prod_{k=2}^{K+1} \text{(a matrix, then an elementwise } f' \text{)} } \vbackerr^{(K+1)}.\]
Write that product out in full, as an alternating chain of matrices \(W^{(k)}_{1:L_{k-1},:}\) and diagonal matrices holding \(f'(\va^{(k)})\). How many factors of \(f'\) does \(\vbackerr^{(1)}\) carry?
Solution 34. Writing \(F^{(k)} = \diagf{f'(\va^{(k)})}\),
\[\vbackerr^{(1)} = F^{(1)} \paren{W^{(2)}_{1:L,:}} F^{(2)} \paren{W^{(3)}_{1:L,:}} \cdots F^{(K)} \paren{W^{(K+1)}_{1:L,:}} \vbackerr^{(K+1)},\]
so \(\vbackerr^{(1)}\) carries \(K\) factors of \(f'\) — one per hidden layer it has passed through. The count is the whole story: every one of them is a number no larger than \(\sigmoid\)’s maximum slope.
Exercise 35 (Find out what the recursion does to the size of a gradient (part 2))
Exercise 35 (Find out what the recursion does to the size of a gradient (part 2)) Use the setup of Exercise 34.
For \(f = \sigmoid\), Exercise 23 gives \(\sigmoid'(a) = \hidval(1-\hidval)\). Give its maximum over all \(a\), and deduce a bound on how much each activation factor can multiply the size of \(\vbackerr\) by.
Solution 35. \(\sigmoid'(a) = \hidval(1-\hidval)\) with \(\hidval \in (0,1)\), maximized at \(\hidval = 1/2\):
\[ \begin{aligned} \max_a \sigmoid'(a) &= \tfrac12 \times \tfrac12 \\ &= \tfrac14. \end{aligned} \]
Each activation factor therefore multiplies the size of \(\vbackerr\) by at most \(1/4\), and by less than that whenever any unit is away from the middle of its range.
Exercise 36 (Find out what the recursion does to the size of a gradient (part 3))
Exercise 36 (Find out what the recursion does to the size of a gradient (part 3)) Use the setup of Exercise 34.
Ignoring the weight matrices for a moment, evaluate that bound for \(K = 2, 5, 10\) and \(20\) hidden layers. State plainly what it says about the gradient reaching layer \(1\) of a twenty-layer sigmoid network.
Solution 36. The best case is \(4^{-K}\):
| \(K\) | \(4^{-K}\) |
|---|---|
| \(2\) | \(6.25 \times 10^{-2}\) |
| \(5\) | \(9.8 \times 10^{-4}\) |
| \(10\) | \(9.5 \times 10^{-7}\) |
| \(20\) | \(9.1 \times 10^{-13}\) |
At twenty layers the gradient reaching layer \(1\) is smaller than the one at the output by twelve orders of magnitude in the most favorable case available — every unit sitting exactly at \(\hidval = 1/2\), which is also the case where the units are carrying the least information. In double precision it is not yet zero; in the single precision used for training, and against the noise of a stochastic gradient, it is effectively zero.
Exercise 37 (Find out what the recursion does to the size of a gradient (part 4))
Exercise 37 (Find out what the recursion does to the size of a gradient (part 4)) Use the setup of Exercise 34.
The bound is worst-case only if the weights do nothing. Say what a weight matrix could do to push in the opposite direction, and why that remedy conflicts with Exercise 23 in the other direction. Then say what goes wrong at the other extreme, when the weights are small: follow a signal forward through the chain rather than a gradient backward, and say what reaches the output.
Solution 37. Large weights. \(W^{(k+1)}\) can multiply \(\vbackerr\) up, and with entries large enough it can exceed the factor of \(1/4\).
But the same large weights drive the pre-activations \(\va^{(k)}\) to large magnitudes on the forward pass, which is precisely where \(\sigmoid'\) collapses: \(\sigmoid'(4) = 0.018\), \(\sigmoid'(6) = 0.0025\). So the remedy makes the problem worse, and by more than it helps. Small weights fail the other way, and the failure is on the forward pass rather than the backward one. A matrix with small entries shrinks \(\va^{(k)}\) towards \(0\), where \(\sigmoid\) is close to affine and every unit reports very nearly \(1/2\) whatever its input was. After ten such layers the network has compressed the inputs into a small neighborhood of one point, so the output carries almost no information about \(\vx\) — and there is nothing useful for the gradient to be large about. The signal shrinks to nothing going forward; with large weights the gradient shrinks to nothing coming back.
The two failure modes are the same phenomenon seen from each end, which is why no choice of scale for a sigmoid network escapes it — there is a best range of scales, it is narrow, and depth narrows it.
Exercise 38 (Find out what the recursion does to the size of a gradient (part 5))
Exercise 38 (Find out what the recursion does to the size of a gradient (part 5)) Use the setup of Exercise 34.
Section 8 measures exactly this, on \(200\) freshly initialized sigmoid networks with \(10\) hidden layers of \(20\) units, \(20\) inputs and \(3\) classes. Averaged over the \(200\), the Frobenius norm of \(\partial \mathcal{L}/\partial W^{(k)}\) runs from \(1.11 \times 10^{-6}\) at layer \(1\) to \(2.12\) at the output — a ratio of \(5.2 \times 10^{-7}\) — and the ratio between consecutive layers has geometric mean \(4.25\). Say which of those numbers is the \(1/4\) of Exercise 35 reappearing, and then describe what training such a network looks like: which parameters move, which do not, and what the network effectively becomes.
Solution 38. The geometric mean of \(4.25\) is the \(1/4\) of Solution 35 reappearing, almost exactly: the bound said each layer costs at most a factor of four, and freshly initialized networks lose very nearly that full factor. The end-to-end ratio \(5.2 \times 10^{-7}\) adds no evidence of its own: the ten consecutive ratios multiply to it, so it is exactly the geometric mean to the power \(-10\). The evidence is the geometric mean itself, which is within 6% of the factor of four the bound allows, so nothing else has a significant effect — the weight matrices are neither restoring the gradient nor making it much worse.
Training such a network moves the last layer or two and leaves the rest where initialization put them. The network is then, in effect, a random feature map with a trained classifier on top: whatever the early layers compute, they compute it by accident, and the only thing being fitted is the last matrix. That is a real model. Fitting a linear classifier on top of a fixed random map is a method in its own right, and one with theory behind it — Rahimi and Recht’s random features, published in 2007, are the best-known version. But it is not the deep network anyone intended to train, and ten sigmoid layers are a very expensive way to draw random features.
Exercise 39 (Find out what the recursion does to the size of a gradient (part 6))
Exercise 39 (Find out what the recursion does to the size of a gradient (part 6)) Use the setup of Exercise 34.
Does making the learning rate larger fix it? Answer for the early layers and the late layers separately.
Solution 39. No, and it fails in both directions at once.
A learning rate large enough to move the early layers is \(10^{6}\) times too large for the last layer, which will diverge immediately. A learning rate that keeps the last layer stable leaves the early layers unchanged. There is no single number that serves both. Fixing that needs either a per-layer step size — the direction the adaptive optimizers of gradient descent take — or an activation function that does not shrink the gradient at all. Exercise 40 takes the second route.
The diagnosis names the cause precisely enough to suggest a remedy: if \(\max_a f'(a)\) is small, use a function whose derivative is not.
Exercise 40 (Change the activation function and see what moves)
Exercise 40 (Change the activation function and see what moves) Exercise 34 traced the problem to a single number: the largest value \(\sigmoid'\) can take. That suggests a simple fix — use an activation function whose derivative is not small.
The rectified linear unit (ReLU) activation function is
\[f(a) = \max(a, 0).\]
Give \(f'(a)\), and say what happens at \(a = 0\).
Solution 40. \[f'(a) = \begin{cases} 1 & a > 0 \\ 0 & a < 0. \end{cases}\]
At \(a = 0\) the function is not differentiable: the left slope is \(0\) and the right slope is \(1\). Implementations pick one by convention, usually \(0\). This is less troubling than it sounds — the set of parameters for which any pre-activation falls exactly on \(0\) has measure zero — but it is worth knowing that the “gradient” your library returns is a subgradient, and that the convergence guarantees of gradient descent were not written for it.
Exercise 41 (Change the activation function and see what moves (part 2))
Exercise 41 (Change the activation function and see what moves (part 2)) Use the setup of Exercise 40.
Substitute your answer into the unrolled product of Exercise 34. What does each activation factor do to the size of \(\vbackerr\) now, and what replaces the \(4^{-K}\) bound?
Solution 41. Each factor is now \(1\) or \(0\): the derivative either passes through unchanged or is switched off entirely.
The product of \(K\) factors is therefore \(1\) along any path whose units are all active, rather than \(4^{-K}\). There is no shrinkage bound to replace it, because there is no shrinkage: the activation function selects which paths carry gradient and does not attenuate the gradient on the paths it selects.
Exercise 42 (Change the activation function and see what moves (part 3))
Exercise 42 (Change the activation function and see what moves (part 3)) Use the setup of Exercise 40.
Repeating Section 8 with ReLU in place of \(\sigmoid\) gives a first-to-last gradient ratio of \(5.6 \times 10^{-2}\) where the sigmoid gave \(5.2 \times 10^{-7}\), and a geometric-mean per-layer ratio of \(1.33\) against \(4.25\). Say which of those four numbers is the one the fix was aimed at, and whether the remaining decay is the same phenomenon.
Solution 42. The ReLU measurement is the Python one with the sigmoid’s two functions swapped for ReLU’s:
Show code
# ReLU in place of the sigmoid: max(a, 0), with derivative [h > 0].
# Where it runs, before "Change the activation function and see what moves",
# this line is hidden, because it gives away part 1 of that exercise, and the
# exercise quotes its result, so it cannot run later. Its solution shows it.
relu_norms = measure(lambda a: np.maximum(a, 0), lambda h: (h > 0) * 1.0)The fix was aimed at the per-layer ratio: \(4.25\) was the \(1/4\) of \(\sigmoid'\) made visible, and \(1.33\) is the ratio once that factor is gone.
The remaining decay is a different phenomenon. It comes from the weight matrices and from the units that happen to be inactive on a given input, not from an upper limit on \(f'\) — so it depends on the initialization scale and can be tuned, where the sigmoid’s \(1/4\) could not be. Over ten layers the difference is \(5.6 \times 10^{-2}\) against \(5.2 \times 10^{-7}\): five orders of magnitude, from one change to one function.
Exercise 43 (Change the activation function and see what moves (part 4))
Exercise 43 (Change the activation function and see what moves (part 4)) Use the setup of Exercise 40.
The fix has a cost. Consider a unit whose pre-activation is negative for every training example. What is its gradient, what will it be at the next step, and what have you lost?
Solution 43. Its gradient is zero, because \(f'(a) = 0\) for every example, so \(\backerr^{\,a}_j = 0\) and none of the weights feeding it move. Since they do not move, the pre-activation is negative for every example at the next step too, and at every step after that.
The unit is dead: it has no output and no prospect of recovering one, and it will consume memory and arithmetic for the rest of training without contributing anything. You have lost the guarantee, which the sigmoid did give, that every parameter receives some gradient from every example. The leaky variants — \(f(a) = \max(a, \leak a)\) for a small \(\leak > 0\), which Maas, Hannun and Ng introduced in 2013 — exist for exactly this reason: a negative pre-activation still passes a small gradient, so no unit can be permanently stuck. They give up the exact zeros that make ReLU fast.
Exercise 44 (Change the activation function and see what moves (part 5))
Exercise 44 (Change the activation function and see what moves (part 5)) Use the setup of Exercise 40.
ReLU is linear on each side of the origin, and the neural networks chapter showed that a linear \(f\) collapses the network to a single matrix. Explain why this argument does not apply, and check your explanation by building the exclusive-or (XOR) function from two ReLU units.
Solution 44. Because ReLU is piecewise linear rather than linear. The neural networks chapter’s collapse needed one affine map valid everywhere; ReLU gives a different affine map on each side of the hyperplane \(\dprod{\vw_j}{\vx} = 0\), and with \(L\) units the weights can cut the plane into as many as \(1 + L + \binom{L}{2}\) pieces, with a different affine function on each.
XOR in two units, on binary inputs:
\[\hidvalof{1} = \max(x_1 + x_2,\; 0), \qquad \hidvalof{2} = \max(x_1 + x_2 - 1,\; 0), \qquad y = \hidvalof{1} - 2\hidvalof{2}.\]
| \(x_1\) | \(x_2\) | \(\hidvalof{1}\) | \(\hidvalof{2}\) | \(y\) | XOR |
|---|---|---|---|---|---|
| \(0\) | \(0\) | \(0\) | \(0\) | \(0\) | \(0\) |
| \(0\) | \(1\) | \(1\) | \(0\) | \(1\) | \(1\) |
| \(1\) | \(0\) | \(1\) | \(0\) | \(1\) | \(1\) |
| \(1\) | \(1\) | \(2\) | \(1\) | \(0\) | \(0\) |
Both hidden units see the same quantity \(x_1 + x_2\) and differ only in their bias; the second one switches on only when both inputs are, and subtracting twice it turns the ramp back down.
Exercise 45 (Change the activation function and see what moves (part 6))
Exercise 45 (Change the activation function and see what moves (part 6)) Use the setup of Exercise 40.
The construction in Exercise 44 is exact, where the sigmoid version in the neural networks chapter was only accurate to about \(4.5 \times 10^{-5}\). Say which property of ReLU gives that, and what ReLU gives up in exchange.
Solution 45. Exactness comes from ReLU being unbounded above and exactly zero below. The sigmoid can only approach \(0\) and \(1\) asymptotically, so its logic gates are always off by about \(\sigmoid(-10) = 4.5 \times 10^{-5}\) at each unit, and no choice of finite weights removes that fixed error. The error does not grow with depth — weights of \(\pm 20\) push each pre-activation about \(10\) away from the boundary, so every stage re-sharpens its input and the final output of the sigmoid network is off by \(4.5440 \times 10^{-5}\) against the \(4.5398 \times 10^{-5}\) of a single unit. The error has a fixed size, and the ReLU construction has none.
ReLU gives up the probabilistic reading. A sigmoid output is in \((0,1)\) and we can report it as a probability. A ReLU output is any non-negative number and means nothing on its own. So we use ReLU for hidden layers and the output layer keeps \(g\) = softmax or \(g\) = \(\sigmoid\): the network computes in unbounded units and converts to probabilities once, at the end.
Watch
Hutchinson covers exploding and vanishing gradients in Deep Neural Networks (21 min) [3].
Read
GeeksforGeeks explains:
9 Exploding gradients
ReLU took away the \(1/4\) cap on each factor of \(f'\). The weights now decide each factor’s size.
Exercise 46 (Turn the weights up)
Exercise 46 (Turn the weights up) Exercise 40 replaced a factor of \(f'\) that was at most \(1/4\) with one that is \(0\) or \(1\). That leaves the weights alone to set how much each layer of the chain in the vanishing-gradient exercise (Exercise 34) multiplies a vector by, and nothing stops that from being more than one.
Draw every weight of layer \(k\) from \(\mathcal{N}(0, s^2 / (L_{k-1} + 1))\), with ReLU units throughout. The measurement of Section 8 used \(s = 1\).
Take a layer with input \(\vhidval\), ignore its bias, and write \(n\) for its fan-in. Give the mean and variance of one pre-activation \(a_j\), and say what ReLU does to it. Deduce that, when the layer is as wide as its input, it multiplies the typical size of \(\vhidval\) by about \(s / \sqrt 2\). Argue that the backward pass multiplies \(\vbackerr\) by the same factor, and say which \(s\) makes that factor \(1\).
Solution 46. Without the bias, \(a_j = \sumin w_{ij} \hidvalof{i}\), a sum of independent terms with mean \(0\). Taking the weight variance as \(s^2/n\) (the \(+1\) added to the fan-in is for the bias, set aside here),
\[ \Ep[a_j] = 0, \qquad \begin{aligned} \Var{a_j} &= \frac{s^2}{n} \sumin \hidvalof{i}^2 \\ &= \frac{s^2\norm{\vhidval}^2}{n}. \end{aligned} \]
ReLU keeps \(a_j\) when it is positive and replaces it with \(0\) otherwise. Since \(a_j\) is symmetric about \(0\), that keeps half its mean square: \(\Ep\sb{\max(a_j, 0)^2} = s^2\norm{\vhidval}^2 / (2n)\). Adding up the \(n\) units of a layer as wide as its input, the output has mean square \(s^2\norm{\vhidval}^2 / 2\), so the layer multiplies the typical size of \(\vhidval\) by \(s/\sqrt 2\).
The backward pass has the same two ingredients: \(\vbackerr^{(k)} = F^{(k)} W^{(k+1)}_{1:L,:} \vbackerr^{(k+1)}\) multiplies by a matrix of the same independent weights, then by \(F^{(k)}\), whose diagonal is \(1\) for the units that are on and \(0\) for the rest. Treating \(\vbackerr^{(k+1)}\) as independent of \(W^{(k+1)}\), and about half the units as on, as He and colleagues do [15, Sec. 2.2], the same count gives the same factor, \(s/\sqrt 2\). Going backward the bias row is dropped, so the factor there is in fact \(\frac{s}{\sqrt 2}\sqrt{L/(L+1)}\), just under \(1\) at \(s = \sqrt 2\): part of why the gradient at \(s = \sqrt 2\) still falls a little towards layer \(1\).
So \(s = \sqrt 2\) makes it \(1\): weights with variance \(2/n\).
Exercise 47 (Turn the weights up (part 2))
Exercise 47 (Turn the weights up (part 2)) Use the setup of Exercise 46.
With \(K = 10\) hidden layers, by what factor does Exercise 46 predict \(\norm{\vhidval^{(K)}}\) to differ from \(\norm{\vx}\) at \(s = 1\), \(\sqrt 2\) and \(2\)?
Solution 47. The factor \((s/\sqrt 2)^K\) is \(2^{-5} = 1/32\) at \(s = 1\), exactly \(1\) at \(s = \sqrt 2\), and \(2^{5} = 32\) at \(s = 2\).
Exercise 48 (Turn the weights up (part 3))
Exercise 48 (Turn the weights up (part 3)) Use the setup of Exercise 46.
Repeating Section 8 with ReLU units at those three values of \(s\) gives these averages over the same \(200\) networks:
| \(s = 1\) | \(s = \sqrt 2\) | \(s = 2\) | |
|---|---|---|---|
| \(\norm{\vx}\) | \(4.44\) | \(4.44\) | \(4.44\) |
| \(\norm{\vhidval^{(K)}}\) | \(0.925\) | \(3.80\) | \(87.1\) |
| \(\norm{\partial \mathcal{L} / \partial W^{(1)}}\) | \(0.0629\) | \(2.21\) | \(91.7\) |
| \(\norm{\partial \mathcal{L} / \partial W^{(K+1)}}\) | \(1.13\) | \(3.59\) | \(91.6\) |
Compare the second row with Exercise 47, and say what keeps the \(s = 1\) figure from falling as far as predicted. Then explain why, at \(s = 2\), the gradient is about the same size at layer \(1\) as at the output, although everything in the network is growing.
Solution 48. The measurement uses the same draws as Section 8. At \(s = 1\) it reproduces the ReLU norms of Exercise 40, which the page checks before quoting any of these figures.
# The measurement with ReLU units, with every weight's standard
# deviation multiplied by `scale`. Besides the per-layer gradient norms it
# returns the average norm and mean square of each layer's output, the norm
# of each network's whole gradient, and the norm of each network's W^(1).
def measure_scaled(scale, K=10, L=20, D=20, C=3, N=200, seed=0):
rng = np.random.default_rng(seed)
widths = [D] + [L] * K + [C]
norms = np.zeros(K + 1)
forward = np.zeros(K + 1)
forward_sq = np.zeros(K + 1)
total = np.zeros(N)
first_weights = np.zeros(N)
for i in range(N):
W = [rng.normal(0, scale / np.sqrt(widths[k] + 1),
(widths[k] + 1, widths[k + 1]))
for k in range(K + 1)]
post = [np.concatenate(([1.0], rng.normal(0, 1, D)))]
pre = []
for k in range(K + 1):
a = W[k].T @ post[k]
pre.append(a)
if k < K:
post.append(np.concatenate(([1.0], np.maximum(a, 0))))
sizes = np.array([np.linalg.norm(p[1:]) for p in post])
forward += sizes
forward_sq += sizes ** 2
e = np.exp(pre[K] - pre[K].max())
t = np.zeros(C)
t[rng.integers(C)] = 1.0
d = e / e.sum() - t
sq = 0.0
for k in range(K, -1, -1):
g = np.linalg.norm(np.outer(post[k], d))
norms[k] += g
sq += g ** 2
if k > 0:
h = post[k][1:]
d = (W[k][1:, :] @ d) * (h > 0)
total[i] = np.sqrt(sq)
first_weights[i] = np.linalg.norm(W[0])
return {"norms": norms / N, "forward": forward / N,
"forward_sq": forward_sq / N, "total": total,
"first_weights": first_weights}
scaled = {name: measure_scaled(scale) for name, scale in
[("one", 1.0), ("root2", np.sqrt(2)), ("two", 2.0)]}The ordering is right, with \(s = \sqrt 2\) holding the size steady and \(s = 2\) growing it by a factor of \(19.6\). But \(s = 1\) shrinks it by a factor of only \(4.80\), where Solution 47 predicted \(32\).
Solution 46 left out the bias. Every layer’s input is \([1;\, \vhidval]\), so every pre-activation keeps a term with variance \(s^2/(L+1)\) however small \(\vhidval\) gets. With it, the mean square \(m_k = \Ep\norm{\vhidval^{(k)}}^2\) obeys
\[m_k = r\,(m_{k-1} + 1), \qquad r = \frac{L s^2}{2(L+1)}.\]
At \(s = 1\), \(r = \tfrac{10}{21}\), and \(m_k\) settles at the fixed point \(r/(1-r) = \tfrac{10}{11}\), a norm of about \(0.953\), which stops the \(s = 1\) network.
Solution 46 also predicts mean squares, where the table averages norms. Run from the measured mean square of \(\vx\), the recursion predicts root-mean-square norms for \(\vhidval^{(K)}\) of \(0.960\), \(4.49\) and \(119\); the networks give \(0.963\), \(4.55\) and \(115\). The average norm at \(s = 2\), \(87.1\), is smaller than its root mean square because the growth varies from network to network, and the networks that grow fastest raise a mean square more than a mean.
The gradient at layer \(k\) is \(\partial \mathcal{L} / \partial W^{(k)} = [1;\, \vhidval^{(k-1)}] \pt{\vbackerr^{(k)}}\). Going forward, \(\vhidval\) gains a factor of \(g = s/\sqrt 2\) per layer; going backward, \(\vbackerr\) gains the same \(g\) per layer. So \(\vhidval^{(k-1)}\) carries \(g^{k-1}\) and \(\vbackerr^{(k)}\) carries \(g^{K+1-k}\), and their product carries \(g^K\) at every layer. So the \(s = 2\) gradient is \(91.7\) at layer \(1\) and \(91.6\) at the output: the growth is the same everywhere, and large. At \(s = 1\) the bias stops \(\vhidval\) from shrinking but nothing stops \(\vbackerr\), so the gradient there still falls towards layer \(1\).
Exercise 49 (Turn the weights up (part 4))
Exercise 49 (Turn the weights up (part 4)) Use the setup of Exercise 46.
At \(s = 2\), \(\norm{W^{(1)}}\) averages \(8.96\). Gradient descent with learning rate \(\learnrate = 0.1\) would move \(W^{(1)}\) by \(\learnrate\norm{\partial \mathcal{L} / \partial W^{(1)}}\) on its first step. Compare the two, and say why the gradient that chose that step says little about where the step lands.
Solution 49. The step has length about \(0.1 \times 91.7 = 9.17\), against the norm \(\norm{W^{(1)}} \approx 8.96\): one step changes \(W^{(1)}\) by about as much as \(W^{(1)}\) is. The gradient describes the loss only near the current weights, as the slope of its linear approximation there. A step that size lands where the approximation no longer holds: units switch on and off, the loss can rise instead of fall, and a few such steps in a row can overflow to Inf and then NaN.
Definition 1 (Exploding gradient)
Definition 1 (Exploding gradient) An exploding gradient is one whose size grows geometrically with the number of layers it passes back through, because the factors of the chain in Exercise 34 are larger than one. A gradient step that large overshoots the region where the gradient describes the loss [16].
Definition 2 (He initialization)
Definition 2 (He initialization) He initialization draws each weight feeding a layer of ReLU units from \(\mathcal{N}(0, 2/n)\), where \(n\) is the layer’s fan-in. It keeps the mean square of the vector passing through a layer as wide as its input unchanged, on the forward pass and on the backward pass [15, Sec. 2.2]. He and colleagues set the biases to zero. The experiment in Exercise 46 draws them too, and its \(s = \sqrt 2\) is the same rule with the bias counted in the fan-in.
Definition 3 (Glorot initialization)
Definition 3 (Glorot initialization) Glorot initialization takes the scale of a layer’s random weights from the layer’s shape rather than from experiment [17]. For a layer with \(n_{\text{in}}\) inputs and \(n_{\text{out}}\) outputs, draw each weight with mean zero and
\[\Var{w} = \frac{2}{n_{\text{in}} + n_{\text{out}}},\]
which is a standard deviation of \(\sqrt{2/(n_{\text{in}} + n_{\text{out}})}\) for a normal draw, or bounds \(\pm\sqrt{6/(n_{\text{in}} + n_{\text{out}})}\) for the uniform draw the original paper states.
Example 1 (A layer from 784 inputs to 256 units)
Example 1 (A layer from 784 inputs to 256 units) Take a layer with \(n_{\text{in}} = 784\) inputs, the pixels of a \(28 \times 28\) image, and \(n_{\text{out}} = 256\) units.
Show code
# The Glorot variance, the standard deviation of the normal version,
# and the half-width of the uniform version, which has the same variance.
n_in <- 784
n_out <- 256
glorot_var <- 2 / (n_in + n_out)
glorot_sd <- sqrt(glorot_var)
glorot_bound <- sqrt(6 / (n_in + n_out))
# The variance He initialization would use, which counts only the inputs.
he_var <- 2 / n_in
# A uniform distribution on (-a, a) has variance a^2 / 3,
# so this equals glorot_var; it confirms the two versions match.
stopifnot(isTRUE(all.equal(glorot_bound^2 / 3, glorot_var)))Glorot initialization draws each of the \(784 \times 256\) weights with variance \(2/(784 + 256) =\) 0.00192, so a normal draw has standard deviation 0.0439, and a uniform draw is on \((-a, a)\) with \(a =\) 0.076. He initialization, which counts only the inputs, would use variance \(2/784 =\) 0.00255. Glorot’s is smaller, because its denominator counts the \(256\) outputs as well as the inputs.
The denominator is a compromise between two requirements that cannot both be met unless the layer is square. Holding the variance of the activations steady from one layer to the next asks for \(1/n_{\text{in}}\); holding the variance of the back-propagated error steady asks for \(1/n_{\text{out}}\). The Glorot variance is the harmonic mean of those two. He initialization makes the same trade for ReLU units, but keeps only the forward requirement, with the factor of \(2\) that the rectifier’s discarded half calls for.
The two schemes set the scale of the weights. Biases are a separate and less important choice, because a bias belongs to one unit and so breaks no symmetry by itself. Zero is the usual value. A common exception starts the biases of rectified linear units at a small positive number, \(0.1\) say, so that the units begin on the active side of the rectifier where the derivative is not zero.
Watch
Hutchinson covers weight initialization (16 min) [3].
Read
GeeksforGeeks explains weight initialization techniques [18].
Exercise 50 (Cap the length of a step)
Exercise 50 (Cap the length of a step) He initialization (Definition 2) chooses a good starting point, but training moves the weights away from it, and a gradient can still grow too large later on. A safeguard should touch only the steps that would do damage.
Stack the gradient of every parameter into one vector \(\vg\), and pick a threshold \(\clipthresh > 0\).
Write a rule that returns a vector pointing the same way as \(\vg\), returns \(\vg\) itself when \(\norm{\vg} \le \clipthresh\), and otherwise returns a vector of length \(\clipthresh\).
Solution 50. \[\vg \leftarrow \minf{1, \frac{\clipthresh}{\norm{\vg}}} \vg.\]
The multiplier is positive, so the direction does not change. It is \(1\) when \(\norm{\vg} \le \clipthresh\), and otherwise it is \(\clipthresh /\norm{\vg}\), which makes the length exactly \(\clipthresh\).
Exercise 51 (Cap the length of a step (part 2))
Exercise 51 (Cap the length of a step (part 2)) Use the setup of Exercise 50.
On the networks of Exercise 46, the norm of \(\vg\) has median \(1.50\), \(6.64\) and \(216\) at \(s = 1\), \(\sqrt 2\) and \(2\), and with \(\clipthresh = 10\) your rule would change \(0\)%, \(34\)% and \(76\)% of them. With \(\learnrate = 0.1\), how long can a step now be?
Solution 51. A step is now at most \(\learnrate \clipthresh = 1\) long.
Exercise 52 (Cap the length of a step (part 3))
Exercise 52 (Cap the length of a step (part 3)) Use the setup of Exercise 50.
How does that length compare with the step length of Exercise 49?
Solution 52. That compares with the \(9.17\) of Solution 49 of Exercise 49, which was about the size of \(W^{(1)}\) itself.
Exercise 53 (Cap the length of a step (part 4))
Exercise 53 (Cap the length of a step (part 4)) Use the setup of Exercise 50.
What do the three percentages of Exercise 51 say about when the rule is the right remedy?
Solution 53. At \(s = 1\) the rule never acts, and at \(s = 2\) it acts on most networks. At \(s = \sqrt 2\) it acts on \(34\)%, which says \(\clipthresh = 10\) is low for those networks. We choose \(\clipthresh\) from the gradient norms that a healthy network actually produces, for example their average over many updates [19, Sec. 3.2]. The rule is a guard against the occasional spike. A network that needs it on most steps is badly scaled, and the remedy is to fix the scale, starting with Definition 2.
Exercise 54 (Cap the length of a step (part 5))
Exercise 54 (Cap the length of a step (part 5)) Use the setup of Exercise 50.
Explain why the rule cannot help with the vanishing gradients of Exercise 34.
Solution 54. The rule only ever shrinks a gradient, and it multiplies every layer’s part of \(\vg\) by the same factor. A vanishing gradient is too small, and its trouble is the ratio between layer \(1\) and the output, which the rule leaves exactly as it was.
Definition 4 (Gradient clipping)
Definition 4 (Gradient clipping) Gradient clipping by norm replaces the gradient \(\vg\) with \(\minf{1, \clipthresh /\norm{\vg}} \vg\) before each update, for a chosen threshold \(\clipthresh\) [16], [19, Sec. 3.2].
Example 2 (Clipping a large gradient)
Example 2 (Clipping a large gradient) Fit mpg on standardized horsepower in the Auto data [20] by gradient descent from \(\vcoef = \vzero\). At the start the gradient of the sample mean squared error (sMSE) is large, because the fitted line is far from the data. Compute it, and clip it to norm \(\clipthresh = 5\):
import numpy as np
import pandas as pd
auto_gc = pd.read_csv("data/auto.csv")[["mpg", "horsepower"]].dropna()
y_gc = auto_gc["mpg"].to_numpy(float)
hp_gc = auto_gc["horsepower"]
X_gc = np.column_stack([np.ones(len(y_gc)), ((hp_gc - hp_gc.mean()) / hp_gc.std()).to_numpy()])
# Gradient of the mean squared error at beta = 0.
g_gc = 2 * X_gc.T @ (X_gc @ np.zeros(2) - y_gc) / len(y_gc)
# Clip by norm: shrink the gradient to length clip_gc if it is longer, keep its direction.
clip_gc = 5
g_clipped_gc = min(1, clip_gc / np.linalg.norm(g_gc)) * g_gc
print(g_gc, np.linalg.norm(g_gc))[-46.89183673 12.12025566] 48.4328911950955[-4.84090827 1.25124222] 5.0auto_gc <- na.omit(read.csv("data/auto.csv")[, c("mpg", "horsepower")])
y_gc <- auto_gc$mpg
X_gc <- cbind(1, as.vector(scale(auto_gc$horsepower)))
# Gradient of the mean squared error at beta = 0.
g_gc <- as.vector(2 * t(X_gc) %*% (X_gc %*% c(0, 0) - y_gc) / length(y_gc))
# Clip by norm: shrink the gradient to length clip_gc if it is longer, keep its direction.
clip_gc <- 5
g_clipped_gc <- min(1, clip_gc / sqrt(sum(g_gc^2))) * g_gc
rbind(gradient = g_gc, clipped = g_clipped_gc) [,1] [,2]
gradient -46.891837 12.120256
clipped -4.840908 1.251242 norm clipped_norm
48.43289 5.00000 The gradient has norm 48.4, far above the threshold, so clipping multiplies it by 0.103. The clipped gradient has norm 5.0 and points the same way. Clipping changes how far a step goes and not which way it goes. Once the fit improves and the gradient norm falls below 5, clipping leaves it alone.
10 Choosing the hidden non-linearity
\(\sigmoid\) and ReLU are two of several functions in use for hidden units. The measurement can compare them too.
Exercise 55 (Compare four non-linearities)
Exercise 55 (Compare four non-linearities) This chapter has used \(\sigmoid\) and ReLU as the non-linearity \(f\), also called the activation function. Two more are common in hidden layers: \(\tanh\), which part 5 of What the logistic function does showed is \(\tanh(a) = 2\sigmoid(2a) - 1\), and
\[f(a) = \max(a, \leak a) \quad\text{for a small } \leak > 0,\]
called the leaky ReLU; take \(\leak = 0.01\).
Use \(\tanh(a) = 2\sigmoid(2a) - 1\) to show that
\[ \begin{aligned} \tanh'(a) &= 4\sigmoid'(2a) \\ &= 1 - \tanh(a)^2. \end{aligned} \]
Give the largest value \(\tanh'\) takes, and compare it with the \(1/4\) that Exercise 35 found for \(\sigmoid'\).
Solution 55. Differentiating \(\tanh(a) = 2\sigmoid(2a) - 1\) gives
\[ \begin{aligned} \tanh'(a) &= 2 \cdot 2\sigmoid'(2a) \\ &= 4\sigmoid'(2a). \end{aligned} \]
Since \(\sigmoid' = \sigmoid(1-\sigmoid)\) and \(\sigmoid(2a) = (1 + \tanh a)/2\),
\[ \begin{aligned} \tanh'(a) &= 4 \cdot \frac{1 + \tanh a}{2} \cdot \frac{1 - \tanh a}{2} \\ &= 1 - \tanh(a)^2. \end{aligned} \]
Its largest value is \(4 \times \tfrac14 = 1\), at \(a = 0\): four times \(\sigmoid\)’s. A factor of the chain in Exercise 34 can now be as large as \(1\), so the bound no longer forces the gradient to shrink.
Exercise 56 (Compare four non-linearities (part 2))
Exercise 56 (Compare four non-linearities (part 2)) Use the setup of Exercise 55.
A unit is saturated when its derivative is close to \(0\), so that almost no gradient passes through it. The derivatives at \(a = 0\), \(2\) and \(5\) are:
| \(a = 0\) | \(a = 2\) | \(a = 5\) | |
|---|---|---|---|
| \(\sigmoid'(a)\) | \(0.250\) | \(0.105\) | \(0.00665\) |
| \(\tanh'(a)\) | \(1.00\) | \(0.0707\) | \(0.000182\) |
For each of the four functions, say for which \(a\) a unit is saturated, and what happens to a unit whose pre-activation is negative for every training example.
Solution 56. \(\sigmoid\) and \(\tanh\) are saturated when \(|a|\) is large, in either direction, and \(\tanh\) more sharply: at \(a = 5\) its derivative is \(0.000182\) against \(\sigmoid\)’s \(0.00665\). A unit of either whose pre-activation is negative on every example still has a nonzero derivative there, so it keeps learning, slowly if it sits far out in the tail.
ReLU’s derivative is exactly \(1\) for \(a > 0\), so it is never saturated there, and exactly \(0\) for \(a < 0\): by the definition, saturated completely. The leaky ReLU’s derivative is \(\leak\) for \(a < 0\): nearly saturated, but not completely. That difference matters for a unit negative on every example. The ReLU unit is dead, with a gradient of \(0\) that stays \(0\) (Solution 58 of Exercise 40); the leaky one still receives \(\leak\) of the gradient and can recover.
Exercise 57 (Compare four non-linearities (part 3))
Exercise 57 (Compare four non-linearities (part 3)) Use the setup of Exercise 55.
Repeating Section 8 with each function in turn gives the ratio of the gradient at layer \(1\) to the gradient at the output, the root mean square (RMS) of the hidden units’ derivatives, and the share of hidden units whose derivative is below \(0.1\):
| \(\sigmoid\) | \(\tanh\) | ReLU | leaky ReLU | |
|---|---|---|---|---|
| layer \(1\) / output | \(5.2 \times 10^{-7}\) | \(0.269\) | \(0.0557\) | \(0.0551\) |
| RMS derivative | \(0.232\) | \(0.846\) | \(0.708\) | \(0.708\) |
| below \(0.1\) | \(0.4\)% | \(0.7\)% | \(49.8\)% | \(49.9\)% |
Explain why \(\tanh\) keeps far more of the gradient than \(\sigmoid\), although both can saturate. Explain ReLU’s RMS derivative, given that about half its units pass nothing, and say whether this comparison is fair to ReLU. Then explain why the leaky ReLU’s figures are almost those of ReLU, and when you would expect the two to differ.
Solution 57. These figures come from the measurement of Section 8, run once for each function on the same draws.
Show code
# The forward pass of the measurement above, on the same draws, recording
# every hidden unit's derivative. It returns their root mean square and the
# fraction of them below flat_cut; measure() supplies the gradient norms.
flat_cut = 0.1
def hidden_slopes(act, dact, K=10, L=20, D=20, C=3, N=200, seed=0):
rng = np.random.default_rng(seed)
widths = [D] + [L] * K + [C]
slopes = []
for _ in range(N):
W = [rng.normal(0, 1 / np.sqrt(widths[k] + 1),
(widths[k] + 1, widths[k + 1]))
for k in range(K + 1)]
post = [np.concatenate(([1.0], rng.normal(0, 1, D)))]
for k in range(K):
h = act(W[k].T @ post[k])
slopes.append(dact(h))
post.append(np.concatenate(([1.0], h)))
rng.integers(C)
slopes = np.concatenate(slopes)
return np.sqrt(np.mean(slopes ** 2)), np.mean(slopes < flat_cut)
leak = 0.01
nonlinearities = {
"sigmoid": (lambda a: 1 / (1 + np.exp(-a)), lambda h: h * (1 - h)),
"tanh": (np.tanh, lambda h: 1 - h ** 2),
"relu": (lambda a: np.maximum(a, 0), lambda h: (h > 0) * 1.0),
"leaky": (lambda a: np.maximum(a, leak * a),
lambda h: np.where(h > 0, 1.0, leak)),
}
activation_results = {}
for name, (act, dact) in nonlinearities.items():
rms, flat = hidden_slopes(act, dact)
activation_results[name] = {"norms": measure(act, dact),
"rms": rms, "flat": flat}Each layer multiplies \(\vbackerr\) by \(F^{(k)}\), the diagonal matrix of its units’ derivatives, so \(\norm{F^{(k)} \vbackerr}^2 = \sum_j f'(a_j)^2 \backerr_j^2\). With the derivatives independent of \(\vbackerr\), that multiplies the size of \(\vbackerr\) by about the root mean square of the derivatives, not their mean.
At this initialization the pre-activations are small, and almost no unit of \(\sigmoid\) or \(\tanh\) is saturated (\(0.4\)% and \(0.7\)%). Near \(a = 0\), though, \(\sigmoid'\) is about \(1/4\) and \(\tanh'\) is about \(1\): \(\sigmoid\) loses most of the gradient at every layer even unsaturated, and \(\tanh\) keeps most of it.
On each input about half of ReLU’s units are on, with derivative \(1\), and half are off, with derivative \(0\). The mean square of the derivatives is \(1/2\), so the root mean square is \(1/\sqrt 2 \approx 0.707\): the \(s/\sqrt 2\) of Solution 55 of Exercise 46 at \(s = 1\). The comparison is not fair to ReLU. The experiment draws weights with variance \(1/(L+1) = 1/21\), close to Glorot and Bengio’s \(2/(L_{k-1} + L_k) = 1/20\) for these equal widths, which is the scale recommended for \(\tanh\) [17, Sec. 4.2]. For ReLU it is half of He initialization’s variance (Definition 2). At He’s scale, \(s = \sqrt 2\), the ReLU networks of Exercise 46 keep \(0.614\) of the gradient between the output and layer \(1\), against \(\tanh\)’s \(0.269\) here.
The leaky ReLU differs from ReLU only through the units that are off. At initialization each of those is off for one input, not dead, and \(\leak\) adds only \(1\)% of the gradient through it, so the figures barely move. The two differ during training, when a unit can end up negative on every example: the ReLU unit stops learning, and the leaky one does not.
Exercise 58 (Compare four non-linearities (part 4))
Exercise 58 (Compare four non-linearities (part 4)) Use the setup of Exercise 55.
Which would you use in the hidden layers of a deep network, and with which initialization?
Solution 58. ReLU is the usual default for the hidden layers of a deep network [21, Sec. 6.3], with He initialization [15, Sec. 2.2]: its derivative is exactly \(1\) wherever a unit is on, it is never saturated on that side, and it is cheap to compute. The leaky ReLU is the fix when too many units die. Of the two functions that saturate on both sides, \(\tanh\) is the better hidden unit [21, Sec. 6.3.2], with Glorot and Bengio’s initialization. \(\sigmoid\) is now rarely used in hidden layers, and belongs at the output of a binary classifier.
Definition 5 (Saturated unit)
Definition 5 (Saturated unit) A unit is saturated when its pre-activation lies where the derivative of its non-linearity is close to \(0\), so that almost no gradient passes back through it. \(\sigmoid\) and \(\tanh\) saturate when \(|a|\) is large; ReLU is saturated completely for \(a < 0\), where its derivative is \(0\).
Example 3 (How much gradient passes through a unit)
Example 3 (How much gradient passes through a unit) The derivative of a non-linearity is the factor that multiplies the gradient on its way back through the unit. Compute it for the sigmoid, \(\tanh\) and ReLU at several pre-activations \(a\):
import numpy as np
a_sat = np.array([-10.0, -5.0, -2.0, 0.0, 2.0, 5.0, 10.0])
sigmoid_sat = 1 / (1 + np.exp(-a_sat))
# Closed forms: sigmoid' = s(1 - s), tanh' = 1 - tanh^2, ReLU' = 1 for a > 0 and 0 otherwise.
d_sigmoid_sat = sigmoid_sat * (1 - sigmoid_sat)
d_tanh_sat = 1 - np.tanh(a_sat) ** 2
d_relu_sat = (a_sat > 0).astype(float)
print(np.column_stack([a_sat, d_sigmoid_sat, d_tanh_sat, d_relu_sat]).round(5))[[-1.0000e+01 5.0000e-05 0.0000e+00 0.0000e+00]
[-5.0000e+00 6.6500e-03 1.8000e-04 0.0000e+00]
[-2.0000e+00 1.0499e-01 7.0650e-02 0.0000e+00]
[ 0.0000e+00 2.5000e-01 1.0000e+00 0.0000e+00]
[ 2.0000e+00 1.0499e-01 7.0650e-02 1.0000e+00]
[ 5.0000e+00 6.6500e-03 1.8000e-04 1.0000e+00]
[ 1.0000e+01 5.0000e-05 0.0000e+00 1.0000e+00]] a sigmoid tanh relu
[1,] -10 0.00005 0.00000 0
[2,] -5 0.00665 0.00018 0
[3,] -2 0.10499 0.07065 0
[4,] 0 0.25000 1.00000 0
[5,] 2 0.10499 0.07065 1
[6,] 5 0.00665 0.00018 1
[7,] 10 0.00005 0.00000 1At \(a = 0\) the sigmoid passes a quarter of the gradient and \(\tanh\) passes all of it. By \(|a| = 5\) the sigmoid passes 0.0066 and \(\tanh\) passes 0.00018, so a unit whose pre-activation sits there is saturated: almost nothing flows back, and the weights feeding it barely learn. ReLU never shrinks the gradient for \(a > 0\), but for every \(a < 0\) the factor is exactly \(0\).
Definition 6 (Leaky ReLU (leaky rectifier))
Definition 6 (Leaky ReLU (leaky rectifier)) The leaky ReLU, or leaky rectifier, is \(f(a) = \max(a, \leak a)\) for a small \(\leak > 0\), often \(0.01\). Its derivative is \(1\) for \(a > 0\) and \(\leak\) for \(a < 0\), so no unit can die [21, Sec. 6.3.1], [22].
Example 4 (A unit that would be dead on half the cars)
Example 4 (A unit that would be dead on half the cars) Let one hidden unit read standardized horsepower in the Auto data [20], with weight \(w = -1\) and bias \(0\), so its pre-activation is \(a = -z\). Count the cars that give the unit a zero gradient under ReLU, and the gradient each car passes under a leaky ReLU with \(\leak = 0.01\):
import numpy as np
import pandas as pd
hp_lr = pd.read_csv("data/auto.csv")["horsepower"].dropna()
a_lr = -((hp_lr - hp_lr.mean()) / hp_lr.std()).to_numpy() # pre-activation: weight -1, bias 0
leak_lr = 0.01
d_relu_lr = (a_lr > 0).astype(float)
d_leaky_lr = np.where(a_lr > 0, 1.0, leak_lr)
print(len(a_lr), (d_relu_lr == 0).sum(), np.unique(d_leaky_lr))392 148 [0.01 1. ]hp_lr <- na.omit(read.csv("data/auto.csv")$horsepower)
a_lr <- -as.vector(scale(hp_lr)) # pre-activation: weight -1, bias 0
leak_lr <- 0.01
d_relu_lr <- as.numeric(a_lr > 0)
d_leaky_lr <- ifelse(a_lr > 0, 1, leak_lr)
c(cars = length(a_lr), relu_zero = sum(d_relu_lr == 0), leaky_values = length(unique(d_leaky_lr))) cars relu_zero leaky_values
392 148 2 For 392 cars, 148 have \(a \le 0\). Under ReLU the unit passes no gradient for any of them, so training on those cars cannot move its weight. Under the leaky ReLU each of them still passes a gradient of 0.01, which is small but not zero, so the weight can recover.
11 Residual connections
ReLU and He initialization keep each factor of the chain near one. The other way out is to give the gradient a route that avoids the chain.
Definition 7 (Residual connection (skip connection, shortcut connection))
Definition 7 (Residual connection (skip connection, shortcut connection)) A residual connection around hidden layer \(k\) adds the layer’s input to its output, so that the layer computes
\[\vhidval^{(k)} = \vhidval^{(k-1)} + f\paren{\va^{(k)}} \qquad\text{instead of}\qquad \vhidval^{(k)} = f\paren{\va^{(k)}},\]
with \(\va^{(k)}\) as before. The layer then learns a change to its input, the residual, rather than a replacement for it [23], [24]. The sum needs the layer’s input and output to have the same width, \(L_{k-1} = L_k\).
Other names: skip connection, shortcut connection.
Example 5 (Twenty layers with and without the skip)
Example 5 (Twenty layers with and without the skip) Stack 20 hidden layers of width 4 with tanh units and small random weights, drawn once in R, and pass one input vector through them. Compare the layers \(\vhidval^{(k)} = f\paren{\va^{(k)}}\) with the residual layers \(\vhidval^{(k)} = \vhidval^{(k-1)} + f\paren{\va^{(k)}}\), using the same weights for both:
set.seed(3)
width_rc <- 4
# Small weights (sd 0.2), so each layer's change to its input is small.
Ws_rc <- lapply(1:20, function(k) matrix(rnorm(width_rc^2, sd = 0.2), width_rc, width_rc))
x_rc <- c(1, -0.5, 0.25, 2)
h_plain_rc <- h_resid_rc <- x_rc
for (W in Ws_rc) {
h_plain_rc <- tanh(W %*% h_plain_rc) # replaces the input
h_resid_rc <- h_resid_rc + tanh(W %*% h_resid_rc) # adds a change to the input
}
# Hand the weights and the input to the Python chunk below, which reads them as ordinary variables.
py$Ws_from_r_rc <- Ws_rc
py$x_from_r_rc <- x_rc
c(input = sqrt(sum(x_rc^2)), plain = sqrt(sum(h_plain_rc^2)), residual = sqrt(sum(h_resid_rc^2))) input plain residual
2.304886e+00 1.639422e-08 2.443157e+00 import numpy as np
# The weights and the input come from the R chunk of this example,
# so that both languages use the same numbers.
Ws_py_rc = [np.asarray(W) for W in Ws_from_r_rc]
x_py_rc = np.asarray(x_from_r_rc)
h_plain_py_rc = h_resid_py_rc = x_py_rc
for W in Ws_py_rc:
h_plain_py_rc = np.tanh(W @ h_plain_py_rc) # replaces the input
h_resid_py_rc = h_resid_py_rc + np.tanh(W @ h_resid_py_rc) # adds a change to the input
print(np.linalg.norm(x_py_rc), np.linalg.norm(h_plain_py_rc), np.linalg.norm(h_resid_py_rc))2.3048861143232218 1.6394217748097892e-08 2.443157472227823The input has norm 2.30. After 20 plain layers it has norm 1.6e-08: each layer shrank the signal, and little of the input survives. The residual stack still has norm 2.44, because every layer starts from its input and only adds a correction. The same holds for gradients going backward, since the identity path contributes a factor of \(1\) at each layer.
Exercise 59 (Give the gradient a path around the chain)
Exercise 59 (Give the gradient a path around the chain) Exercise 34 found that \(\vbackerr^{(1)}\) carries one factor of \(f'\) for every hidden layer, and Exercise 40 made each factor larger. A residual connection (Definition 7) attacks the count instead.
Put a residual connection around every hidden layer from \(k = 2\) to \(K\), with all hidden widths equal to \(L\), and write \(\vhidgrad^{(k)} = \partial \mathcal{L} / \partial \vhidval^{(k)}\) for the gradient with respect to a layer’s output.
Differentiate \(\vhidval^{(k)} = \vhidval^{(k-1)} + f(\va^{(k)})\) to get the recursion that takes \(\vhidgrad^{(k)}\) to \(\vhidgrad^{(k-1)}\), by way of \(\vbackerr^{(k)} = \partial \mathcal{L} / \partial \va^{(k)}\). Compare it with the recursion for the same network without the residual connection.
Solution 59. Write \(W^{(k)}_{1:L,:}\) for \(W^{(k)}\) without its bias row, as in Solution 34, and \(F^{(k)} = \diag(f'(\va^{(k)}))\). The output of the block depends on its input twice: directly, and through \(\va^{(k)}\). So
\[\frac{\partial \vhidval^{(k)}}{\partial \vhidval^{(k-1)}} = I + F^{(k)} \pt{W^{(k)}_{1:L,:}},\]
and the chain rule gives the recursion
\[\vbackerr^{(k)} = f'\paren{\va^{(k)}} \odot \vhidgrad^{(k)}, \qquad \vhidgrad^{(k-1)} = \vhidgrad^{(k)} + W^{(k)}_{1:L,:} \, \vbackerr^{(k)}.\]
Without the residual connection the second equation loses its first term: \(\vhidgrad^{(k-1)} = W^{(k)}_{1:L,:} \, \vbackerr^{(k)}\). The residual connection hands the layer below a copy of the gradient it received, and adds the usual matrix-and-\(f'\) term on top.
Exercise 60 (Give the gradient a path around the chain (part 2))
Exercise 60 (Give the gradient a path around the chain (part 2)) Use the setup of Exercise 59.
Unroll your recursion from \(\vhidgrad^{(K)}\) down to \(\vhidgrad^{(1)}\) as a product of \(K - 1\) matrix factors. Expand the product, and find the term that carries no factor of \(f'\) at all. How many factors of \(f'\) does \(\vbackerr^{(1)}\) then carry along that term, and how many did it carry in Exercise 34?
Solution 60. Substituting the first equation into the second, \(\vhidgrad^{(k-1)} = \paren{I + W^{(k)}_{1:L,:} F^{(k)}} \vhidgrad^{(k)}\), so
\[\vhidgrad^{(1)} = \paren{I + W^{(2)}_{1:L,:} F^{(2)}} \paren{I + W^{(3)}_{1:L,:} F^{(3)}} \cdots \paren{I + W^{(K)}_{1:L,:} F^{(K)}} \vhidgrad^{(K)}.\]
Expanding the product gives \(2^{K-1}\) terms, one for each choice of which factors contribute their \(W F\) part and which contribute \(I\). Choosing \(I\) every time gives \(\vhidgrad^{(K)}\) itself, with no factor of \(f'\) and no weight matrix: the gradient arriving at the top of the stack reaches layer \(1\) at full size. The other terms add to it; the shortest of them passes through one block.
Then \(\vbackerr^{(1)} = F^{(1)} \vhidgrad^{(1)}\) carries one factor of \(f'\) along that term, where Exercise 34’s \(\vbackerr^{(1)}\) carried \(K\). Along that term the bound \(4^{-K}\) becomes \(4^{-1}\), whatever the depth.
Exercise 61 (Give the gradient a path around the chain (part 3))
Exercise 61 (Give the gradient a path around the chain (part 3)) Use the setup of Exercise 59.
With \(f = \mathrm{ReLU}\), set every weight in \(W^{(2)}, \ldots, W^{(K)}\) to zero. What does the network compute, and what is \(\vhidgrad^{(1)}\)? Now do the same with \(f = \sigmoid\), and say what each block adds to its input instead.
Solution 61. With \(f = \mathrm{ReLU}\) and \(W^{(k)} = 0\), each block has \(\va^{(k)} = \vzero\) and \(f(\vzero) = \vzero\), so \(\vhidval^{(k)} = \vhidval^{(k-1)}\): every block is the identity, and the network computes exactly what its first hidden layer and output layer compute on their own. The recursion of Solution 59 gives \(\vhidgrad^{(k-1)} = \vhidgrad^{(k)}\) at every block, so \(\vhidgrad^{(1)} = \vhidgrad^{(K)}\). A residual network can start as a shallow one and use its extra layers only where they help. A plain stack has no such default: to pass its input through unchanged, each layer would have to learn the identity from its weights.
With \(f = \sigmoid\) the same weights give \(\sigmoid(0) = 1/2\), so each block adds \(1/2\) to every unit, and the \(K - 1 = 9\) blocks of the measurement add \(4.5\) to each entry of \(\vhidval^{(1)}\) between them.
Exercise 62 (Give the gradient a path around the chain (part 4))
Exercise 62 (Give the gradient a path around the chain (part 4)) Use the setup of Exercise 59.
Repeating Section 8 on the same \(200\) sigmoid networks, with a residual connection around each of hidden layers \(2\) to \(10\), moves the gradient norm at layer \(1\) from \(1.11 \times 10^{-6}\) to \(0.95\). The geometric-mean ratio between consecutive hidden layers, leaving out the output layer, falls from \(4.22\) to \(1.10\), and the gradient norm at the output layer rises from \(2.12\) to \(21.1\). Say which of these numbers Exercise 60 predicted, and why the hidden layers’ gradients are still not all equal. Then use Exercise 61 to explain the rise at the output layer, and say what a transformer pairs with each residual connection to stop it.
Solution 62. The measurement uses the same draws as Section 8, so the two runs see the same \(200\) networks, inputs and targets. Without the residual connections it reproduces the norms printed there, which the page checks before quoting either set.
Show code
# The measurement of the vanishing-gradient experiment, on the same draws, with an
# optional residual connection around hidden layers 2 to K. Besides the
# per-layer gradient norms it returns the average norm of [1; h^(K)], the
# input to the output layer.
def measure_residual(act, dact, residual, K=10, L=20, D=20, C=3, N=200,
seed=0):
rng = np.random.default_rng(seed)
widths = [D] + [L] * K + [C]
norms = np.zeros(K + 1)
stream = 0.0
for _ in range(N):
W = [rng.normal(0, 1 / np.sqrt(widths[k] + 1),
(widths[k] + 1, widths[k + 1]))
for k in range(K + 1)]
post = [np.concatenate(([1.0], rng.normal(0, 1, D)))]
pre, outs = [], []
for k in range(K + 1):
a = W[k].T @ post[k]
pre.append(a)
if k < K:
out = act(a)
outs.append(out)
h = out + post[k][1:] if residual and k >= 1 else out
post.append(np.concatenate(([1.0], h)))
stream += np.linalg.norm(post[K])
e = np.exp(pre[K] - pre[K].max())
t = np.zeros(C)
t[rng.integers(C)] = 1.0
d = e / e.sum() - t
norms[K] += np.linalg.norm(np.outer(post[K], d))
eps = W[K][1:, :] @ d
for k in range(K - 1, -1, -1):
d = eps * dact(outs[k])
norms[k] += np.linalg.norm(np.outer(post[k], d))
if k > 0:
eps = (eps if residual else 0) + W[k][1:, :] @ d
return norms / N, stream / N
sigmoid = lambda a: 1 / (1 + np.exp(-a))
dsigmoid = lambda h: h * (1 - h)
plain_norms, plain_stream = measure_residual(sigmoid, dsigmoid, False)
res_norms, res_stream = measure_residual(sigmoid, dsigmoid, True)Solution 60 predicted the first number. The gradient at layer \(1\) grows by a factor of \(8.6 \times 10^{5}\), from \(1.11 \times 10^{-6}\) to \(0.95\), and the ratio between consecutive hidden layers falls from \(4.22\), the \(1/4\) of Exercise 34 made visible, to \(1.10\). These ratios leave out the output layer, so the first differs a little from the \(4.25\) of Section 8; with the residual connections the output layer’s gradient rises for a separate reason, below.
The hidden layers’ gradients still differ, for two reasons. Each \(\partial \mathcal{L} / \partial W^{(k)} = [1;\, \vhidval^{(k-1)}] \pt{\vbackerr^{(k)}}\) still has its own factor of \(f'\) in \(\vbackerr^{(k)}\), which Solution 60 cannot remove because it sits inside the block, not on the path around it. And its other factor, \(\vhidval^{(k-1)}\), gets larger towards the output. Layer \(1\) is the exception: its input is \(\vx\) itself, whose entries are not squashed into \((0, 1)\).
That growth is Solution 61’s sigmoid case. Every block adds an output in \((0, 1)\) to the stream, so the input to the output layer has average norm \(22.8\) with the residual connections, against \(2.5\) without them. The output layer’s gradient is that input times \(\vbackerr^{(K+1)}\), so it rises with it, from \(2.12\) to \(21.1\). A deeper stack would push it further. A transformer follows each residual connection with a layer normalization, which rescales each position’s vector to mean \(0\) and variance \(1\) before a learned scale and shift; this is the “Add & Norm” of Attention and Transformers.
Exercise 63 (Give the gradient a path around the chain (part 5))
Exercise 63 (Give the gradient a path around the chain (part 5)) Use the setup of Exercise 59.
The sum in Definition 7 needs \(L_{k-1} = L_k\). Suggest how to keep a path around a layer that changes the width.
Solution 63. Replace the identity on the path around the layer with a matrix:
\[\vhidval^{(k)} = P^{(k)} \vhidval^{(k-1)} + f\paren{\va^{(k)}},\]
where \(P^{(k)}\) is \(L_k \times L_{k-1}\), either learned or fixed. He and colleagues describe a learned linear projection of this kind and a fixed one that pads the shorter vector with zeros, keeping the identity wherever the width does not change [23, Secs. 3.2–3.3], and find the learned projection only slightly better [23, Sec. 4.1]. The path around the layer then multiplies the gradient by \(\pt{P^{(k)}}\) once, rather than by a chain of matrices and \(f'\) factors.
Remark 1 (Shattered gradients)
Remark 1 (Shattered gradients). A gradient can be the right size and still be useless. Balduzzi and colleagues compared the gradients of deep rectifier networks at nearby inputs. In a plain network the correlation between them falls like \(2^{-L}\) with depth \(L\), until the gradient resembles white noise. Residual connections like those of Definition 7 slow this to \((3/4)^{L}\), and with batch normalization added the correlation falls only like \(1/\sqrt{L}\) [25, Secs. 3.1–3.3]. Descent on such a gradient makes no consistent progress, and checking the gradient’s norm, which catches vanishing and exploding gradients, does not catch this.
12 A whole minibatch at once
Everything so far handled one training example. Gradient descent already argued for minibatches on statistical grounds — a gradient averaged over \(B\) examples is a better estimate than a gradient from one.
A second argument decides the batch size in practice. The same formulas, with vectors replaced by matrices, run far faster per example than a loop over the same examples would.
Exercise 64 (Do a whole minibatch at once)
Exercise 64 (Do a whole minibatch at once) Gradient descent established why training uses minibatches rather than single examples. This exercise does the bookkeeping: every formula so far handled one example, and a minibatch of \(B\) of them should not mean running those formulas \(B\) times in a loop.
Stack the examples as columns. Write \(\mpredvar \in \reals^{(\ninputvars+1) \times B}\) for the inputs and \(T \in \reals^{C \times B}\) for the one-hot targets, and take the loss to be the mean over the batch, \(\mathcal{L} = \frac{1}{B}\sumto{i}{B} \mathcal{L}_i\).
Capital letters denote the batched version of the matching lower-case vector: \(A\), \(H\) and \(Y\) for \(\va\), \(\vhidval\) and \(\vy\). The output pre-activation \(\v{b}\) is the one exception — its capital is already the batch size, so write it \(Z \in \reals^{C \times B}\).
Rewrite the four forward equations of Exercise 6 with \(\mpredvar\) in place of \(\vx\), and give the shape of each result. Say what happens to the bracket \([1; \vhidval]\).
Solution 64. \[A = \tp{W} \mpredvar \;\in \reals^{L \times B}, \qquad H = f(A) \;\in \reals^{L \times B},\] \[Z = \tp{U} \begin{bmatrix} \tp{\vone} \\ H \end{bmatrix} \;\in \reals^{C \times B}, \qquad Y = g(Z) \;\in \reals^{C \times B}.\]
The bracket becomes a row of ones stacked on top of \(H\), giving an \((L+1) \times B\) matrix: each column still carries its own constant \(1\). \(f\) is elementwise as before; \(g\) is softmax applied down each column separately, since each column is one example’s own distribution over the \(C\) classes.
Exercise 65 (Do a whole minibatch at once (part 2))
Exercise 65 (Do a whole minibatch at once (part 2)) Use the setup of Exercise 64.
\(\mbackerr^{b} = \partial \mathcal{L} / \partial Z\) collects the output errors of the whole batch. Give it, including whatever the \(\frac{1}{B}\) contributes.
Solution 65. Differentiating the mean rather than the sum divides everything by \(B\):
\[\mbackerr^{b} = \frac{1}{B}\paren{Y - T} \;\in \reals^{C \times B}.\]
Column \(i\) is that example’s \(\vbackerr^{\,b}\), scaled by \(1/B\). Carrying the \(1/B\) here rather than at the end is the usual choice, because it means every gradient downstream comes out already averaged.
Exercise 66 (Do a whole minibatch at once (part 3))
Exercise 66 (Do a whole minibatch at once (part 3)) Use the setup of Exercise 64.
Give \(\partial \mathcal{L} / \partial U\) and \(\partial \mathcal{L} / \partial W\) as matrix products. Confirm that they are the same shape as \(U\) and \(W\) — not \(B\) times larger.
Solution 66. \[\frac{\partial \mathcal{L}}{\partial U} = \begin{bmatrix} \tp{\vone} \\ H \end{bmatrix} \pt{\mbackerr^{b}} \;\in \reals^{(L+1) \times C}, \qquad \frac{\partial \mathcal{L}}{\partial W} = \mpredvar \pt{\mbackerr^{a}} \;\in \reals^{(\ninputvars+1) \times L},\]
with \(\mbackerr^{a}\) from Solution 69. \((L+1) \times B\) against \(B \times C\) gives \((L+1) \times C\), and \((\ninputvars+1) \times B\) against \(B \times L\) gives \((\ninputvars+1) \times L\): both match their parameter matrix exactly.
They must. There is one gradient per parameter no matter how many examples produced it.
Exercise 67 (Do a whole minibatch at once (part 4))
Exercise 67 (Do a whole minibatch at once (part 4)) Use the setup of Exercise 64.
In the single-example version those two gradients were outer products. Say what happened to the sum over the batch, and where the batched expression now performs it.
Solution 67. In Exercise 67, the sum over the batch became the inner dimension of the matrix product.
\([\tp{\vone}; H]\pt{\mbackerr^{b}}\) computes \(\frac{1}{B}\sumto{i}{B} [1; \vhidval_i]\pt{\vbackerr^{\,b}_i}\), because a matrix product sums over the shared index — and here the shared index is the example. The \(\frac{1}{B}\) is the one Solution 65 included in \(\mbackerr^{b}\), so what comes out is the average over the batch rather than the sum. The loop is still there. We have given it to the routine that is best at it. Every optimization anyone has made to dense matrix multiplication, over several decades, now applies to averaging gradients over a batch.
Exercise 68 (Do a whole minibatch at once (part 5))
Exercise 68 (Do a whole minibatch at once (part 5)) Use the setup of Exercise 64.
Exercise 14 found that one example gives a rank-one gradient for \(U\). Give the largest possible rank of \(\partial \mathcal{L}/\partial U\) for a batch of \(B\), and say what batch size is needed before that limit stops applying.
Solution 68. A product of an \((L+1) \times B\) matrix with a \(B \times C\) one has rank at most \(\min(L+1, C, B)\).
With \(C = 10\) outputs, the batch stops being the limiting factor once \(B \geq 10\). Below that, the batch itself limits how many directions in \(U\) a single step can move, regardless of how informative the examples are — which is a second reason, alongside the arithmetic of Exercise 70, that batches of one are a poor default.
Exercise 69 (Do a whole minibatch at once (part 6))
Exercise 69 (Do a whole minibatch at once (part 6)) Use the setup of Exercise 64.
Write the whole backward pass out, and check that every line is a matrix product or an elementwise product, with no explicit loop over examples anywhere in it.
Solution 69. \[\mbackerr^{b} = \tfrac{1}{B}(Y - T), \qquad \frac{\partial \mathcal{L}}{\partial U} = \begin{bmatrix} \tp{\vone} \\ H \end{bmatrix}\pt{\mbackerr^{b}},\] \[\mbackerr^{\hidval} = U_{1:L,:}\,\mbackerr^{b}, \qquad \mbackerr^{a} = \mbackerr^{\hidval} \odot f'(A), \qquad \frac{\partial \mathcal{L}}{\partial W} = \mpredvar \pt{\mbackerr^{a}}.\]
Three matrix products and one elementwise product. This is the whole of back-propagation for a one-hidden-layer network, and it is the form an implementation actually takes: no loop over examples appears anywhere, and the only thing that changed from the single-example version is that vectors became matrices with \(B\) columns.
Exercise 70 (Work out why a batch is faster per example)
Exercise 70 (Work out why a batch is faster per example) Exercise 64 shows that a batch changes no formula and no gradient — \(B\) examples produce the average of the \(B\) single-example gradients, as a loop would have. So the batch gives no mathematical benefit. It gives an arithmetic benefit, and this exercise says what.
Keep \(\ninputvars = 784\), \(L = 100\), \(C = 10\), so \(P = 79{,}510\), and suppose the parameters are stored in \(32\)-bit floats.
Give the number of multiply-accumulates in the forward pass for a batch of \(B\), and the number of bytes of weights that have to be read to perform it. Say which of the two depends on \(B\).
Solution 70. Every weight is used once per example, so the forward pass is \(PB = 79{,}510B\) multiply-accumulates.
The weights are read once regardless of \(B\): \(4P = 318{,}040\) bytes. The arithmetic scales with the batch; the memory traffic for the weights does not.
Exercise 71 (Work out why a batch is faster per example (part 2))
Exercise 71 (Work out why a batch is faster per example (part 2)) Use the setup of Exercise 70.
Form the ratio: how many floating-point operations does the hardware perform per byte of weights it reads? (Count a multiply-accumulate as two operations.) Give the formula and evaluate it at \(B = 1, 8, 32, 128\).
Solution 71. Two operations per multiply-accumulate gives
\[\frac{2PB}{4P} = \frac{B}{2} \quad\text{operations per byte},\]
which is independent of the network’s size entirely — the \(P\) cancels.
| \(B\) | operations per byte of weights |
|---|---|
| \(1\) | \(0.5\) |
| \(8\) | \(4\) |
| \(32\) | \(16\) |
| \(128\) | \(64\) |
Exercise 72 (Work out why a batch is faster per example (part 3))
Exercise 72 (Work out why a batch is faster per example (part 3)) Use the setup of Exercise 70.
An accelerator has some fixed ratio of arithmetic throughput to memory bandwidth, in operations per byte. Explain, in terms of your Exercise 71 answer, why small batches leave it idle and what changes as \(B\) grows.
Solution 72. With \(B = 1\) the hardware performs half an operation for every byte it fetches, so it spends nearly all its time waiting on memory: the arithmetic units are idle, and the step takes as long as reading the weights takes.
As \(B\) grows, the same fetched weights serve more examples, and the ratio rises linearly. Once \(B/2\) exceeds the machine’s own operations-per-byte ratio, the weights arrive faster than the arithmetic consumes them and the step is bounded by compute — which is the regime the hardware was designed for. Below that point, increasing \(B\) costs almost no extra wall-clock time, because the extra arithmetic fills gaps that were previously idle.
Solution 70 counted only the weights, though, and Solution 73 is about to show that the activations outweigh them. The activations scale with \(B\), so they cancel out of the ratio the way the weights do not, and the true intensity saturates rather than rising without limit: with \(895\) floats of activation per example, it is
\[\frac{2PB}{4P + 4 \cdot 895 B} \;\xrightarrow{\;B \to \infty\;}\; \frac{2P}{4 \cdot 895} \approx 44 \quad\text{operations per byte},\]
reaching \(26\) at \(B = 128\) rather than the \(64\) the table gives. A network of this shape therefore cannot be driven compute-bound on an accelerator whose own ratio exceeds about \(44\), however large the batch. The direction of Solution 71’s argument survives — raising \(B\) does raise the intensity, steeply at first — but the upper limit is set by the shape of the network, not by the batch size.
This is the real content of “minibatching is faster”: a batch is more arithmetic per step and less wall-clock time per example, because a matrix-matrix product reuses data that a matrix-vector product fetches and throws away.
Exercise 73 (Work out why a batch is faster per example (part 4))
Exercise 73 (Work out why a batch is faster per example (part 4)) Use the setup of Exercise 70.
Exercise 8 counted the stored activations. At \(B = 128\), compare the bytes of activations with the bytes of weights, and say which one bounds the batch size in practice.
Solution 73. At \(B = 128\) the stored activations are \(128 \times (784 + 100 + 10 + 1) \times 4 = 458{,}240\) bytes, against \(318{,}040\) bytes of weights.
The activations are already larger, and unlike the weights they grow with \(B\) — and, in a deep network, with the number of layers as well, since every layer’s activations are held until the backward pass reaches them. So activation memory, not the model, sets the largest batch that fits, and a practitioner meets that limit first.
Exercise 74 (Work out why a batch is faster per example (part 5))
Exercise 74 (Work out why a batch is faster per example (part 5)) Use the setup of Exercise 70.
Against these benefits, a larger batch means fewer updates for the same amount of arithmetic. The minibatch gradient is an average of \(B\) independent estimates, so say how its standard error scales with \(B\), and how the work per step scales. What does the comparison say about doubling the batch size repeatedly?
Solution 74. The batch gradient is a mean of \(B\) independent unbiased estimates (gradient descent), so its standard error falls like \(1/\sqrt{B}\) while the work per step rises like \(B\):
| \(B\) | noise | work per step |
|---|---|---|
| \(1\) | \(1\) | \(1\times\) |
| \(4\) | \(0.5\) | \(4\times\) |
| \(16\) | \(0.25\) | \(16\times\) |
| \(64\) | \(0.125\) | \(64\times\) |
Quadrupling the batch halves the noise and quadruples the cost. Doubling repeatedly therefore has sharply diminishing returns: past some point you are paying linearly for a square-root improvement in an estimate that was never the bottleneck — gradient descent tolerates a noisy direction, as gradient descent shows.
Exercise 75 (Work out why a batch is faster per example (part 6))
Exercise 75 (Work out why a batch is faster per example (part 6)) Use the setup of Exercise 70.
Put Exercise 72 and Exercise 74 together into a statement about how to choose \(B\), and say which part of it is a fact about the estimator and which is a fact about the machine.
Solution 75. Raise \(B\) until the arithmetic intensity is near its upper limit, and no further than the memory allows; past that point extra examples give \(\sqrt{\cdot}\) and cost \(B\).
The \(1/\sqrt{B}\) is a fact about sampling: the variance of a mean of \(B\) independent estimates is \(\varpar/B\), which is an exact identity at every \(B\) and needs only that the variance is finite. It would hold on any machine ever built. The operations-per-byte ratio and the activation-memory limit are facts about the machine, and they change with every hardware generation. So batch size is the hyperparameter that moves when the hardware changes — and it is a useful habit to ask of any training recipe which of its numbers are facts about the problem and which are facts about the machine. The two are not independent: the common rules for changing \(B\) scale the learning rate with it, linearly or as \(\sqrt{B}\), precisely because a less noisy gradient can support a longer step.
Exercise 76 (Fit a big batch in a small memory)
Exercise 76 (Fit a big batch in a small memory) Solution 70 found that the activations of a batch, which grow with the batch size \(B\), are what limit \(B\) in practice. Suppose the activations of only \(m < B\) examples fit in memory at once, yet you want the gradient step of a batch of \(B\).
Describe a way to take that step using \(m\) examples at a time, and say what you add up and when you take the step.
Solution 76. Split the batch into groups of at most \(m\) examples. For each group in turn, run forward and backward propagation, and add the group’s summed per-example gradients to a running total. Then free that group’s activations and move to the next group. Take one gradient step only after the last group, with the total divided by \(B\). This is gradient accumulation. Memory holds the activations of \(m\) examples, plus one extra array the size of the parameters for the total.
Exercise 77 (Fit a big batch in a small memory (part 2))
Exercise 77 (Fit a big batch in a small memory (part 2)) Use the setup of Exercise 76.
Show that your procedure gives exactly the gradient of the full batch.
Solution 77. The batch gradient is an average of per-example gradients, which is a sum divided by \(B\):
\[ \grad J(\vcoef) = \frac{1}{B} \sumto{i}{B} \grad \obsloss{i}(\vcoef). \]
A sum can be split into groups, \(\mathcal{G}_1, \dots, \mathcal{G}_K\), and added in any order:
\[ \sumto{i}{B} \grad \obsloss{i}(\vcoef) = \sumto{k}{K} \sum_{i \in \mathcal{G}_k} \grad \obsloss{i}(\vcoef). \]
Each group’s inner sum needs only that group’s activations, and the parameters do not change until the step, so the accumulated total equals the full-batch sum exactly, up to floating-point rounding. It is the same step a memory large enough for \(B\) examples would take.
Exercise 78 (Fit a big batch in a small memory (part 3))
Exercise 78 (Fit a big batch in a small memory (part 3)) Use the setup of Exercise 76.
Suppose \(B\) is not a multiple of \(m\), so the last group is smaller. What goes wrong if you average the group means instead of what you added up in Exercise 76?
Solution 78. Each group’s mean divides by its own size \(n_k\). Averaging the \(K\) group means gives every group the weight \(1/K\), so an example in the small last group counts for more than one in a full group. The result is not the batch gradient. Dividing the accumulated sum by \(B\), as in Solution 76, weights every example by \(1/B\).
Exercise 79 (Fit a big batch in a small memory (part 4))
Exercise 79 (Fit a big batch in a small memory (part 4)) Use the setup of Exercise 76.
Check Exercise 77 and Exercise 78 for the sample mean squared error (sMSE) of a linear model on \(B = 100\) simulated examples with \(m = 32\).
Solution 79.
# Simulate B = 100 examples of a linear model with 3 coefficients.
set.seed(571)
B <- 100
m <- 32
X <- matrix(rnorm(B * 3), nrow = B)
y <- X %*% c(1, -2, 0.5) + rnorm(B)
coef <- c(0, 0, 0)
# The summed gradient of the squared error over some rows, at the current
# coefficients: the sum over examples of 2 x (x'coef - y).
grad_sum <- function(rows) {
2 * crossprod(X[rows, , drop = FALSE], X[rows, , drop = FALSE] %*% coef - y[rows])
}
# The full-batch gradient, all rows at once.
full <- grad_sum(1:B) / B
# Accumulate over groups of at most m rows, then divide by B.
total <- numeric(3)
means <- list()
starts <- seq(1, B, by = m)
for (s in starts) {
rows <- s:min(s + m - 1, B)
total <- total + grad_sum(rows)
means[[length(means) + 1]] <- grad_sum(rows) / length(rows)
}
accumulated <- total / B
last_n <- B - m * (length(starts) - 1)
# Averaging the group means gives every group the same weight.
mean_of_means <- Reduce(`+`, means) / length(means)
isTRUE(all.equal(c(full), accumulated))[1] FALSE[1] FALSEimport numpy as np
# Simulate B = 100 examples of a linear model with 3 coefficients.
rng = np.random.default_rng(571)
B, m = 100, 32
X = rng.normal(0, 1, size=(B, 3))
y = X @ np.array([1, -2, 0.5]) + rng.normal(0, 1, size=B)
coef = np.zeros(3)
# The summed gradient of the squared error over some rows, at the current
# coefficients: the sum over examples of 2 x (x'coef - y).
def grad_sum(rows):
return 2 * X[rows].T @ (X[rows] @ coef - y[rows])
# The full-batch gradient, all rows at once.
full = grad_sum(slice(0, B)) / B
# Accumulate over groups of at most m rows, then divide by B.
total = np.zeros(3)
means = []
for s in range(0, B, m):
rows = slice(s, min(s + m, B))
total += grad_sum(rows)
means.append(grad_sum(rows) / len(y[rows]))
accumulated = total / B
# Averaging the group means gives every group the same weight.
mean_of_means = np.mean(means, axis=0)
print(np.allclose(full, accumulated), np.allclose(full, mean_of_means))True FalseThe accumulated gradient matches the full-batch gradient, and the average of the group means does not. Here the batch splits into 3 full groups of 32 and a last group of 4, so with the group means averaged, each example in the last group counts 6.25 times as heavily as in the batch.
A caveat: this equivalence needs the loss of each example to depend only on that example. A layer that computes statistics across the batch, such as batch normalization, couples the examples, so accumulating small groups does not reproduce the full-batch step.
Watch
- The CS231N lecture shows vectorized backpropagation, from 1:05:07 [11], Lecture 4.
- Starmer explains the matrix algebra of a network [26].
Watch
Hutchinson covers minibatched backprop (11 min) [3].
13 Gradients with respect to the input
Back-propagation differentiates the loss with respect to the weights, but the chain rule works for any variable it is applied to. Run the same backward sweep one step further, across \(W\), and it returns the gradient with respect to the input:
\[\frac{\partial \mathcal{L}}{\partial \vx_{1:\ninputvars}} = W_{1:\ninputvars,\,:}\;\vbackerr^{\,a}.\]
This is the crossing Exercise 18 made across \(U\), with the row for the constant dropped for the same reason: the constant is not an input anyone can change.
That derivative answers a different question — not “how should the model change?” but “which parts of this input is the prediction sensitive to?” Starting the sweep from the indicator of one class in place of \(\vy - \vt\) gives the gradient of that class’s score instead of the loss’s.
Definition 8 (Saliency map)
Definition 8 (Saliency map) For an image classifier, a saliency map is the gradient of the predicted class’s score with respect to each pixel of one input image, displayed as a picture. A large magnitude at a pixel means a small change to that pixel would move the prediction noticeably.
Example 6 (Which pixels decide a 3 versus an 8)
Example 6 (Which pixels decide a 3 versus an 8) The digits data of scikit-learn [27] are \(8 \times 8\) grey-scale images of handwritten digits. Fit a logistic regression that tells 3s from 8s, and take one 8. The gradient of its predicted probability of “8” with respect to the 64 pixels is the saliency map of Definition 8. For a logistic regression, \(p = \sigmoidf{\dprod{\vw}{\vx} + b}\), so the gradient is \(p(1 - p)\,\vw\):
Show code
import matplotlib.pyplot as plt
import numpy as np
from sklearn.datasets import load_digits
from sklearn.linear_model import LogisticRegression
digits_sm = load_digits()
keep_sm = np.isin(digits_sm.target, [3, 8])
X_sm = digits_sm.data[keep_sm] / 16.0 # pixel values scaled to [0, 1]
y_sm = (digits_sm.target[keep_sm] == 8).astype(int)
clf_sm = LogisticRegression(max_iter=1000).fit(X_sm, y_sm)
# Take the first image of an 8, and compute p(1 - p) * w, the gradient of p with respect to the pixels.
i_sm = int(np.flatnonzero(y_sm == 1)[0])
p_sm = clf_sm.predict_proba(X_sm[i_sm:i_sm + 1])[0, 1]
saliency_sm = p_sm * (1 - p_sm) * clf_sm.coef_[0]
fig, axes = plt.subplots(1, 2, figsize=(6, 3), layout="constrained")
axes[0].imshow(X_sm[i_sm].reshape(8, 8), cmap="Greys");
axes[0].set_title("an 8");
bound_sm = np.abs(saliency_sm).max()
axes[1].imshow(saliency_sm.reshape(8, 8), cmap="coolwarm", vmin=-bound_sm, vmax=bound_sm);
axes[1].set_title("saliency map");
for ax in axes:
ax.set_xticks([]);
ax.set_yticks([]);
plt.show()
The model gives this image a probability of 0.985 of being an 8. Because the model is linear in the pixels, the map is the weight vector times the constant \(p(1 - p)\), so its pattern is the same for every image and only its size changes with the confidence: it is largest for images the model is unsure about, and nearly zero for images it is sure of. For a deep network, where the gradient depends on the image, the map differs from image to image.
Read
GeeksforGeeks explains saliency maps [28].
A saliency map costs one backward sweep and one extra matrix-vector product, because the computation already exists for training.
14 Checking the gradient in code
The code below implements the formulas of this chapter and then checks them. The code nudges each weight by \(\pm\smallconst\) and re-evaluates the loss, which gives a central finite-difference estimate of that weight’s partial derivative. The analytic gradient should agree with it to several decimal places. The check earns its place because a wrong gradient is not a loud failure: as long as it points roughly downhill it still descends, so the loss falls, the run looks healthy, and the only symptom is that it converges slower or to a worse optimum than it should. Nothing in the training curve distinguishes that from a hard problem.
The problem is XOR, which Neural Networks uses to motivate the hidden layer: four points, not linearly separable, and so unlearnable by anything without one.
import numpy as np
def forward(W, U, X):
a = X @ W # a = W'x, one row per example
h = 1 / (1 + np.exp(-a)) # h = f(a), logistic
hb = np.column_stack([np.ones(len(h)), h])
b = hb @ U # b = U'[1; h]
b = b - b.max(axis=1, keepdims=True)
p = np.exp(b)
return h, hb, p / p.sum(axis=1, keepdims=True) # y = g(b), softmax
def loss(W, U, X, t):
*_, y = forward(W, U, X)
return -np.mean(np.log(y[np.arange(len(t)), t]))
def gradients(W, U, X, t):
h, hb, y = forward(W, U, X)
d_b = y.copy()
d_b[np.arange(len(t)), t] -= 1 # softmax + cross-entropy
d_b /= len(t)
d_h = d_b @ U[1:].T # drop the bias row
d_a = d_h * h * (1 - h) # f'(a) = h(1-h)
return X.T @ d_a, hb.T @ d_b
rng = np.random.default_rng(0)
X = np.array([[1., 0, 0], [1., 0, 1], [1., 1, 0], [1., 1, 1]]) # bias column
t = np.array([0, 1, 1, 0]) # XOR
W = rng.normal(scale=1.0, size=(3, 4))
U = rng.normal(scale=1.0, size=(5, 2))
# Does the derivation match the function it claims to differentiate?
gW, gU = gradients(W, U, X, t)
eps, worst = 1e-6, 0.0
for M, G in ((W, gW), (U, gU)):
for i in range(M.shape[0]):
for j in range(M.shape[1]):
keep = M[i, j]
M[i, j] = keep + eps; hi = loss(W, U, X, t)
M[i, j] = keep - eps; lo = loss(W, U, X, t)
M[i, j] = keep
worst = max(worst, abs((hi - lo) / (2 * eps) - G[i, j]))
print("largest gap, analytic vs finite difference:", f"{worst:.2e}")largest gap, analytic vs finite difference: 1.57e-100 1.4194
5000 0.0019
10000 0.0009
15000 0.0006
20000 0.0004predictions: [0 1 1 0] targets: [0 1 1 0]# nolint start: object_name_linter. W, U, X and M match the prose.
forward <- function(W, U, X) {
a <- X %*% W # a = W'x, one row per example
h <- 1 / (1 + exp(-a)) # h = f(a), logistic
hb <- cbind(1, h)
b <- hb %*% U # b = U'[1; h]
b <- b - apply(b, 1, max)
p <- exp(b)
list(h = h, hb = hb, y = p / rowSums(p)) # y = g(b), softmax
}
loss <- function(W, U, X, t) {
y <- forward(W, U, X)$y
-mean(log(y[cbind(seq_along(t), t)]))
}
gradients <- function(W, U, X, t) {
f <- forward(W, U, X)
d_b <- f$y
d_b[cbind(seq_along(t), t)] <- d_b[cbind(seq_along(t), t)] - 1
d_b <- d_b / length(t)
d_h <- d_b %*% t(U[-1, , drop = FALSE]) # drop the bias row
d_a <- d_h * f$h * (1 - f$h) # f'(a) = h(1-h)
list(W = t(X) %*% d_a, U = t(f$hb) %*% d_b)
}
set.seed(0)
X <- cbind(1, c(0, 0, 1, 1), c(0, 1, 0, 1)) # bias column
t <- c(1L, 2L, 2L, 1L) # XOR, as 1-based class labels
W <- matrix(rnorm(12), 3, 4)
U <- matrix(rnorm(10), 5, 2)
# Does the derivation match the function it claims to differentiate?
g <- gradients(W, U, X, t)
eps <- 1e-6
worst <- 0
for (nm in c("W", "U")) {
M <- get(nm)
for (i in seq_len(nrow(M))) for (j in seq_len(ncol(M))) {
keep <- M[i, j]
M[i, j] <- keep + eps
assign(nm, M)
hi <- loss(W, U, X, t)
M[i, j] <- keep - eps
assign(nm, M)
lo <- loss(W, U, X, t)
M[i, j] <- keep
assign(nm, M)
worst <- max(worst, abs((hi - lo) / (2 * eps) - g[[nm]][i, j]))
}
}
cat("largest gap, analytic vs finite difference:",
format(worst, digits = 3), "\n")largest gap, analytic vs finite difference: 6.44e-11 0 0.6982
5000 0.0019
10000 9e-04
15000 6e-04
20000 4e-04 predictions: 1 2 2 1 targets: 1 2 2 1 using Random, Statistics
function forward(W, U, X)
a = X * W # a = W'x, one row per example
h = 1 ./ (1 .+ exp.(-a)) # h = f(a), logistic
hb = hcat(ones(size(h, 1)), h)
b = hb * U # b = U'[1; h]
b = b .- maximum(b, dims = 2)
p = exp.(b)
return h, hb, p ./ sum(p, dims = 2) # y = g(b), softmax
end;
function loss(W, U, X, t)
_, _, y = forward(W, U, X)
-mean(log(y[i, t[i]]) for i in eachindex(t))
end;
function gradients(W, U, X, t)
h, hb, y = forward(W, U, X)
d_b = copy(y)
for i in eachindex(t)
d_b[i, t[i]] -= 1 # softmax + cross-entropy
end
d_b ./= length(t)
d_h = d_b * U[2:end, :]' # drop the bias row
d_a = d_h .* h .* (1 .- h) # f'(a) = h(1-h)
return X' * d_a, hb' * d_b
end;
function gradient_check(W, U, X, t; eps = 1e-6)
gW, gU = gradients(W, U, X, t)
worst = 0.0
for (M, G) in ((W, gW), (U, gU)), i in axes(M, 1), j in axes(M, 2)
keep = M[i, j]
M[i, j] = keep + eps; hi = loss(W, U, X, t)
M[i, j] = keep - eps; lo = loss(W, U, X, t)
M[i, j] = keep
worst = max(worst, abs((hi - lo) / (2 * eps) - G[i, j]))
end
return worst
end;
function train!(W, U, X, t; h = 0.5, steps = 20_000)
for step in 0:steps
gW, gU = gradients(W, U, X, t)
W .-= h .* gW
U .-= h .* gU
step % 5_000 == 0 && println(step, " ", round(loss(W, U, X, t), digits = 4))
end
end;
Random.seed!(0);
# bias column
X = [1.0 0 0; 1.0 0 1; 1.0 1 0; 1.0 1 1];
# XOR, as 1-based class labels
t = [1, 2, 2, 1];
W = randn(3, 4);
U = randn(5, 2);
# Does the derivation match the function it claims to differentiate?
println("largest gap, analytic vs finite difference: ", gradient_check(W, U, X, t))largest gap, analytic vs finite difference: 1.0225707780531224e-100 0.6997
5000 0.0021
10000 0.001
15000 0.0006
20000 0.0005predictions: [1, 2, 2, 1] targets: [1, 2, 2, 1]In all three, the analytic gradient agrees with the finite difference to about ten decimal places. That is roughly the best a central difference can do, and the reason is worth knowing, because shrinking \(\smallconst\) — the obvious way to sharpen the agreement — makes the agreement worse. Two errors pull in opposite directions as \(\smallconst\) shrinks. Truncation error — the \(O(\smallconst^2)\) remainder in the Taylor expansion the difference quotient comes from — falls. Cancellation error rises: the two losses agree in their leading digits, so subtracting them throws those digits away, and dividing the remainder by a smaller \(\smallconst\) magnifies what is left. Sweeping \(\smallconst\) over each tab’s check — that block with eps replaced and nothing else changed — gives
| worst gap | \(\smallconst = 10^{-2}\) | \(10^{-4}\) | \(10^{-6}\) | \(10^{-8}\) | \(10^{-10}\) | \(10^{-12}\) |
|---|---|---|---|---|---|---|
| Python | \(9.1\times10^{-7}\) | \(9.2\times10^{-11}\) | \(1.6\times10^{-10}\) | \(2.0\times10^{-8}\) | \(1.2\times10^{-6}\) | \(1.5\times10^{-4}\) |
| R | \(7.8\times10^{-7}\) | \(7.8\times10^{-11}\) | \(6.4\times10^{-11}\) | \(6.9\times10^{-9}\) | \(7.6\times10^{-7}\) | \(7.2\times10^{-5}\) |
| Julia | \(4.1\times10^{-7}\) | \(4.2\times10^{-11}\) | \(1.0\times10^{-10}\) | \(1.0\times10^{-8}\) | \(7.6\times10^{-7}\) | \(9.7\times10^{-5}\) |
The digits differ across the three, because each draws its own weights. The shape does not, and the shape is the point: every row is a valley rather than a slope, bottoming out at \(\smallconst =\) \(10^{-4}\) for Python, \(10^{-6}\) for R and \(10^{-4}\) for Julia, and climbing from there, until at \(10^{-12}\) the check is less accurate than it was at \(10^{-2}\).
A gradient check that comes back disappointing is therefore not evidence of a wrong gradient until we have varied \(\smallconst\). And \(\smallconst = 10^{-6}\), the value the three tabs above use, sits near the floor. Its worst gap is within a factor of \(2.5\) of the best in every row. That is a small margin, and that is the useful part: at \(10^{-6}\) the two error terms are still close to balanced, so the check is near its best and the penalty for the habitual value is mild. Cancellation runs away from \(10^{-8}\) onward, and that costs orders of magnitude.
The three tabs print different losses at step 0 for one reason only: each language’s generator produces different initial weights from the same scale-\(1\) normal. The size of the gradient gap and the final predictions do reproduce. All three recover XOR exactly, though they spell the answer differently. Python labels its classes from \(0\), and R and Julia label from \(1\). The label base is an indexing convention and nothing more: both conventions select the same two softmax columns, so it cannot move the loss.
Two limits of this example are worth naming. Four points is not a data set, so there is nothing here to overfit and no held-out set to check against — success means memorizing the training data, which is the one setting where that is the goal. And full-batch descent over four points is not how we train a network in practice. See A whole minibatch at once above.
Versions
Rendered with NumPy 2.5.4 on Python 3.13.18, R 4.6.1, and Julia 1.13.3. The gradient-check magnitude is a property of the arithmetic and should reproduce anywhere; the per-step losses depend on each language’s random number generator and will not match across the three tabs.
15 Where this goes next
This chapter leads to three topics.
Nobody writes these formulas any more. Solution 33 names the general procedure: reverse-mode automatic differentiation, which computes the derivative of one scalar with respect to many inputs in a constant multiple of the cost of evaluating it. A library records the operations of the forward pass and replays them backwards, so a new layer needs only its own local derivative and everything above composes automatically. The hand derivation here is worth doing once, because it tells you which of a library’s error messages are about shapes and which are about your model.
Initialization stopped being arbitrary. The vanishing-gradient exercise (Exercise 34) and Section 9 show that the scale of the initial weights decides whether a deep network trains at all: too small and the signal shrinks going forward, too large and sigmoid units saturate or ReLU gradients explode. The standard schemes — Glorot and Bengio’s for sigmoid and tanh units [17, Sec. 4.2] and He initialization (Definition 2) for ReLU — set the variance of each layer’s weights from its width, so that the size of the activations is roughly preserved from layer to layer. The experiment in Section 8 scales its weights by \(1/\sqrt{\text{fan-in}+1}\) for the same reason, so that it measures depth rather than one fixed scale. For ReLU that scale still shrinks the signal, and the factor of \(\sqrt 2\) in He initialization holds it steady.
Depth is where this chapter stops. Every network here is a stack of fully connected layers. The architectures that made deep learning work — convolutional layers for images, attention for sequences, residual connections that give the gradient a path around the activation functions entirely — are all answers to questions this chapter raises. Section 11 works through residual connections, Convolutional and Recurrent Networks builds convolutional and recurrent layers, and Attention and Transformers covers attention.
16 Further reading
Goodfellow et al. [21], Section 6.5
Glorot and Bengio [17] — why deep sigmoid networks were hard to train
Starmer [29], on backpropagation, with the chain rule worked through in pictures
Watch
Andrew Ng’s lectures for DeepLearning.AI work through backpropagation from computation graphs: the computation graph [30]; derivatives on the graph [31]; the intuition for a network [32]; the forward and backward passes of a deep network [33]; vanishing and exploding gradients [34]; and checking a gradient numerically [35].