Calculus in one variable, plus linear algebra, equals the calculus that machine learning actually runs on. The gradient points uphill, the Jacobian is a derivative that is a matrix, the Hessian is curvature — and the multivariable chain rule is backpropagation.
A loss depends on millions of parameters at once. The trick is to change one at a time.
For \(f:\RR^d\to\RR\), the partial derivative with respect to \(x_i\) is
\[ \frac{\partial f}{\partial x_i}(\vx) = \lim_{h\to0}\frac{f(\vx+h\ve_i)-f(\vx)}{h}, \]where \(\ve_i\) is the \(i\)-th standard basis vector. In practice: differentiate with respect to \(x_i\), treating every other variable as a constant.
It is Lecture 4's definition verbatim, with the increment restricted to one axis. Every rule from Lecture 5 applies unchanged. For \(f(x,y)=x^2y+3y\):
\[\frac{\partial f}{\partial x} = 2xy, \qquad \frac{\partial f}{\partial y} = x^2+3.\]
The gradient collects all the partials into one vector:
\[ \grad f(\vx) = \begin{bmatrix}\partial f/\partial x_1 \\ \vdots \\ \partial f/\partial x_d\end{bmatrix} \in \RR^{d}. \]Shapes matter: \(f\) maps \(\RR^d\to\RR\), and \(\grad f\) maps \(\RR^d\to\RR^d\). The gradient of a scalar function has the same shape as its input — which is what makes \(\vw \leftarrow \vw - \eta\grad L\) type-check.
Partials tell you the rate of change along the axes. But you can move in any direction. Let \(\vu\) be a unit vector; the rate of change of \(f\) along \(\vu\) is
\[D_{\vu}f(\vx) = \lim_{h\to0}\frac{f(\vx+h\vu)-f(\vx)}{h}.\]
For a vector-valued \(\vf:\RR^n\to\RR^m\), the Jacobian is the \(m\times n\) matrix of all partials:
\[ \mJ_{\vf}(\vx) = \frac{\partial \vf}{\partial \vx} = \begin{bmatrix} \partial f_1/\partial x_1 & \cdots & \partial f_1/\partial x_n\\ \vdots & \ddots & \vdots\\ \partial f_m/\partial x_1 & \cdots & \partial f_m/\partial x_n \end{bmatrix}, \qquad (\mJ)_{ij} = \frac{\partial f_i}{\partial x_j}. \]Row \(i\) is \(\grad f_i^\top\). When \(m=1\) the Jacobian is a single row — the gradient, transposed. That row-versus-column distinction is a notational nuisance and the source of endless sign-and-shape confusion; the rule to hold on to is the Jacobian has as many rows as outputs and as many columns as inputs.
The Jacobian is the best linear approximation to \(\vf\) near \(\vx\): \(\;\vf(\vx+\vh) \approx \vf(\vx) + \mJ\vh\). Compare Lecture 4's \(f(x+h)\approx f(x)+f'(x)h\) — the derivative was a number that multiplied the nudge; now it is a matrix that multiplies the nudge vector. Same idea, richer type.
If \(\vf:\RR^n\to\RR^m\) and \(\vg:\RR^m\to\RR^p\), then for \(\vh = \vg\circ\vf\),
\[ \mJ_{\vh}(\vx) = \mJ_{\vg}\big(\vf(\vx)\big)\;\mJ_{\vf}(\vx). \]Shapes: \((p\times m)(m\times n) = (p\times n)\) ✓. The chain rule is matrix multiplication. In coordinates,
\[ \frac{\partial h_i}{\partial x_k} = \sum_{j=1}^{m}\frac{\partial g_i}{\partial u_j}\,\frac{\partial f_j}{\partial x_k}, \]summing over every path from input \(k\) to output \(i\) through the intermediate layer.
Now look at a network. With \(\vz^{(\ell)} = f_\ell(\vz^{(\ell-1)})\) and a scalar loss \(L\) at the end, the chain rule gives
\[ \frac{\partial L}{\partial \vz^{(0)}} = \frac{\partial L}{\partial \vz^{(L)}}\,\mJ_{L}\,\mJ_{L-1}\cdots\mJ_{1}. \]A product of Jacobians, exactly as Lecture 5's scalar chain rule was a product of derivatives. Two observations that explain the entire design of a deep learning framework:
Dual class from Lecture 5. Reverse mode wins whenever there is one output and
many inputs, which is exactly the shape of a loss function.Lecture 10 used two gradient rules on credit. Here they are, proved by writing out coordinates — which is always the safe way to settle a matrix-calculus question you are unsure about.
Using the two identities above, differentiate
\[ J(\vw) = \|\mX\vw-\vy\|^2 \]with respect to \(\vw\), and confirm you recover the normal equations. Check the shape of every term as you go.
Expand (as in Lecture 10):
\[J(\vw) = \vw^\top\underbrace{(\mX^\top\mX)}_{d\times d,\ \text{symmetric}}\vw \;-\;2\,\vw^\top\underbrace{(\mX^\top\vy)}_{d\times 1}\;+\;\|\vy\|^2 .\]Differentiate term by term, using linearity:
Alternative route, worth knowing. Treat it as a composition: \(\vr(\vw)=\mX\vw-\vy\) with Jacobian \(\mJ_{\vr}=\mX\) (\(n\times d\)), and \(J=\vr^\top\vr\) with \(\grad_{\vr}J = 2\vr\) (\(n\times1\)). The chain rule for a scalar output reads \(\grad_{\vw}J = \mJ_{\vr}^\top\,\grad_{\vr}J = \mX^\top(2\vr) = 2\mX^\top(\mX\vw-\vy)\) ✓. Same answer, and note the transpose: pulling a gradient backwards through a linear map multiplies by \(\mJ^\top\), which is exactly what a backward pass does.
For \(f:\RR^d\to\RR\), the Hessian is the \(d\times d\) matrix of second partials:
\[ \mH = \grad^2 f, \qquad H_{ij} = \frac{\partial^2 f}{\partial x_i \partial x_j}. \]It is the Jacobian of the gradient. Clairaut's theorem: if the second partials are continuous then \(\frac{\partial^2 f}{\partial x_i\partial x_j} = \frac{\partial^2 f}{\partial x_j\partial x_i}\), so \(\mH\) is symmetric — and therefore everything from Lecture 11 applies to it: real eigenvalues, orthogonal eigenvectors, a definiteness classification.
The second-order Taylor expansion in \(d\) dimensions is exactly what you would guess:
\[f(\vx+\vh) \;\approx\; f(\vx) + \grad f(\vx)^\top\vh + \tfrac12\,\vh^\top\mH(\vx)\,\vh.\]
Three terms: a number, a linear form, a quadratic form. Compare Lecture 6's \(f+f'h+\tfrac12 f''h^2\) — identical, with the scalars promoted to a vector and a matrix. And now the second-derivative test writes itself, using Lecture 11's classification.
Let \(\grad f(\vx^\star)=\vzero\) and let \(\mH=\grad^2f(\vx^\star)\), with eigenvalues \(\lambda_1,\dots,\lambda_d\).
Why: at a critical point the linear term vanishes, so \(f(\vx^\star+\vh)-f(\vx^\star) \approx \tfrac12\vh^\top\mH\vh\), and Lecture 11 showed the sign of that quadratic form is decided entirely by the signs of the eigenvalues.
Saddle points dominate in high dimensions. Here is a counting argument that is worth internalising. At a critical point of a random smooth function in \(\RR^d\), think of each of the \(d\) Hessian eigenvalues as being positive or negative roughly independently. For a local minimum you need all \(d\) to be positive — probability about \(2^{-d}\). For \(d=10^6\) that number is indistinguishable from zero.
So essentially every critical point of a large network's loss is a saddle. Consequences:
Compute \(\grad f\) and \(\mH\) for each, and classify every critical point:
(a) \(f(x,y)=x^2+3y^2-4x+6y\) (b) \(f(x,y)=x^3-3xy+y^3\) (c) \(f(x,y)=e^{-(x^2+y^2)}\)
(a) \(\grad f = (2x-4,\ 6y+6)\), zero at \((2,-1)\). \(\mH=\begin{bmatrix}2&0\\0&6\end{bmatrix}\), constant. Eigenvalues \(2,6\), both positive → strict local (and, since \(\mH\) is constant and positive definite, global) minimum, value \(f(2,-1)=4+3-8-6=-7\).
(b) \(\grad f = (3x^2-3y,\ -3x+3y^2)\). Setting both to zero: \(y=x^2\) and \(x=y^2\).
Substituting, \(x=x^4\), so \(x(x^3-1)=0\), giving \(x=0\) or \(x=1\). Critical points \((0,0)\) and
\((1,1)\).
\(\mH=\begin{bmatrix}6x&-3\\-3&6y\end{bmatrix}\).
At \((0,0)\): \(\begin{bmatrix}0&-3\\-3&0\end{bmatrix}\), \(\det=-9<0\). A negative determinant in 2D means
\(\lambda_1\lambda_2<0\), hence mixed signs → saddle. (Eigenvalues \(\pm3\).)
At \((1,1)\): \(\begin{bmatrix}6&-3\\-3&6\end{bmatrix}\), \(\tr=12>0\) and \(\det=36-9=27>0\), so both
eigenvalues are positive (\(3\) and \(9\)) → local minimum, value \(1-3+1=-1\).
Useful 2×2 shortcut: \(\det<0\) ⟹ saddle; \(\det>0\) with \(\tr>0\) ⟹ minimum; \(\det>0\) with
\(\tr<0\) ⟹ maximum.
(c) \(\grad f = -2e^{-(x^2+y^2)}(x,y)\), zero only at \((0,0)\) since the exponential
never vanishes.
At the origin, \(\mH = \begin{bmatrix}-2&0\\0&-2\end{bmatrix}\) (differentiate again and evaluate; the
cross terms vanish by symmetry). Both eigenvalues negative → local maximum, value 1. It is
also global, since \(0<f\le1\) everywhere.
Prove that \(\grad f(\vx^\star)=\vzero\) is necessary for \(\vx^\star\) to be a local minimum of a differentiable \(f:\RR^d\to\RR\).
Hint: reduce to the one-dimensional case using \(g(t)=f(\vx^\star+t\vu)\) and Lecture 7's Fermat rule.
Then explain why the converse fails in a way that is worse in high dimensions than in one.
Proof. Suppose \(\vx^\star\) is a local minimum, so \(f(\vx^\star)\le f(\vx)\) for all \(\vx\) within some \(\delta\). Fix any unit direction \(\vu\) and define \(g(t)=f(\vx^\star+t\vu)\), a function of one real variable, differentiable at 0 with \(g'(0)=D_{\vu}f(\vx^\star)=\grad f(\vx^\star)^\top\vu\) by the directional-derivative theorem.
For \(|t|<\delta\) the point \(\vx^\star+t\vu\) is within \(\delta\) of \(\vx^\star\), so \(g(t)\ge g(0)\): the point \(t=0\) is a local minimum of \(g\), and it is interior. By Fermat's rule (Lecture 7), \(g'(0)=0\), i.e.
\[\grad f(\vx^\star)^\top\vu = 0.\]This holds for every unit \(\vu\). Take \(\vu=\ve_i\) for each \(i\): every component of the gradient is zero. Hence \(\grad f(\vx^\star)=\vzero\). ∎
(Equivalently: a vector orthogonal to every vector must be the zero vector — take \(\vu\) parallel to the gradient itself and you get \(\|\grad f\|^2=0\).)
Why the converse fails worse in high dimensions. In one dimension a critical point is a minimum, a maximum, or an inflection — and generically (with \(f''\neq0\)) it is one of the first two, so a critical point has roughly a 50% chance of being a minimum.
In \(d\) dimensions the classification is by the sign pattern of \(d\) eigenvalues. A minimum needs all of them positive. If the signs were independent and equally likely, minima would have probability \(2^{-d}\), while saddles — anything with at least one sign of each — would have probability \(1-2^{-d+1}\).
The independence assumption is crude, and random-matrix theory sharpens it considerably (the semicircle law, and results of Bray–Dean and Dauphin et al. on spin glasses and neural losses), but the qualitative conclusion survives and is well supported empirically: in high dimensions almost all critical points are saddles, and the ones with low loss are overwhelmingly the ones that are minima. This reframes what optimisation is fighting. It is not searching a landscape littered with bad local minima; it is navigating long flat plateaus around saddles.
Derive the gradient of the logistic regression loss
\[ L(\vw) = -\sum_{i=1}^{n}\Big[y_i\log\sigma(\vw^\top\vx_i) + (1-y_i)\log\big(1-\sigma(\vw^\top\vx_i)\big)\Big] \]and show it simplifies to \(\;\grad L = \mX^\top(\hat{\vp}-\vy)\), where \(\hat p_i = \sigma(\vw^\top\vx_i)\).
Comment on the resemblance to the least-squares gradient.
Work with one term. Write \(z=\vw^\top\vx\), \(p=\sigma(z)\), and \(\ell = -[y\log p + (1-y)\log(1-p)]\).
Step 1 — differentiate the loss with respect to \(p\):
\[\frac{\partial \ell}{\partial p} = -\frac{y}{p} + \frac{1-y}{1-p} = \frac{-y(1-p)+p(1-y)}{p(1-p)} = \frac{p-y}{p(1-p)}.\]Step 2 — the sigmoid derivative (Lecture 5): \(\dfrac{dp}{dz} = p(1-p)\).
Step 3 — chain rule, and the beautiful cancellation:
\[\frac{\partial \ell}{\partial z} = \frac{p-y}{p(1-p)}\cdot p(1-p) = p - y.\]Step 4 — chain rule again, with \(\dfrac{\partial z}{\partial \vw}=\vx\):
\[\grad_{\vw}\ell = (p-y)\,\vx.\]Summing over the dataset and stacking, \(\grad L = \sum_i (\hat p_i-y_i)\vx_i = \mX^\top(\hat{\vp}-\vy)\). ∎
The resemblance. Least squares gave \(\grad J = 2\mX^\top(\mX\vw-\vy) = 2\mX^\top (\hat{\vy}-\vy)\). Both are design matrix transposed, times prediction error. This is not a coincidence: both losses are negative log-likelihoods of exponential-family models with a canonical link function (identity for the Gaussian, logit for the Bernoulli), and for that whole family the gradient always takes this form. It is why the same optimiser code works for both, and why the update rule \(\vw \leftarrow \vw - \eta\mX^\top(\hat{\vy}-\vy)\) feels like the natural thing to write.
Note also that the messy \(p(1-p)\) cancelled exactly. Had we used squared error on top of a sigmoid instead, it would not have — the gradient would carry a factor of \(\sigma'\), which vanishes when the model is confidently wrong. That is a concrete reason cross-entropy trains better than squared error for classification: it removes the saturation from the gradient.
Coding. Implement reverse-mode automatic differentiation for a two-layer network and verify it against numerical gradients.
num_grad to at least \(10^{-7}\) relative error.2. The backward pass. Work from the loss backwards, carrying \(\bar{\vv} := \partial L/\partial \vv\) at each node. Each step is a vector-Jacobian product.
\[\bar{\vz}_2 = \vz_2-\vy, \qquad \bar{\mW}_2 = \bar{\vz}_2\va_1^\top, \qquad \bar{\vb}_2 = \bar{\vz}_2,\] \[\bar{\va}_1 = \mW_2^\top\bar{\vz}_2, \qquad \bar{\vz}_1 = \bar{\va}_1 \odot (1-\va_1^2), \] \[\bar{\mW}_1 = \bar{\vz}_1\vx^\top, \qquad \bar{\vb}_1 = \bar{\vz}_1.\]Two patterns to notice, and they are the whole of backpropagation. Pulling a gradient back through a linear map multiplies by the transpose of the weight matrix. Pulling it back through an elementwise map multiplies elementwise by the local derivative. The gradient with respect to a weight matrix is always an outer product of (incoming gradient) and (stored input).
import numpy as np
rng = np.random.default_rng(0)
din, dh, dout = 5, 7, 3
W1 = rng.normal(size=(dh, din))*0.5; b1 = rng.normal(size=dh)*0.1
W2 = rng.normal(size=(dout, dh))*0.5; b2 = rng.normal(size=dout)*0.1
x = rng.normal(size=din); y = rng.normal(size=dout)
def forward(W1, b1, W2, b2):
z1 = W1 @ x + b1; a1 = np.tanh(z1)
z2 = W2 @ a1 + b2
return 0.5*np.sum((z2 - y)**2), (z1, a1, z2)
def backward(W1, b1, W2, b2):
_, (z1, a1, z2) = forward(W1, b1, W2, b2)
gz2 = z2 - y # dL/dz2
gW2 = np.outer(gz2, a1); gb2 = gz2
ga1 = W2.T @ gz2 # transpose: linear layer
gz1 = ga1 * (1 - a1**2) # elementwise: tanh
gW1 = np.outer(gz1, x); gb1 = gz1
return gW1, gb1, gW2, gb2
def num_grad(f, p, h=1e-6):
g = np.zeros_like(p); it = np.nditer(p, flags=['multi_index'])
while not it.finished:
i = it.multi_index; old = p[i]
p[i] = old + h; fp = f()
p[i] = old - h; fm = f()
p[i] = old; g[i] = (fp - fm)/(2*h)
it.iternext()
return g
gW1, gb1, gW2, gb2 = backward(W1, b1, W2, b2)
for name, param, ana in [("W1",W1,gW1), ("b1",b1,gb1), ("W2",W2,gW2), ("b2",b2,gb2)]:
num = num_grad(lambda: forward(W1,b1,W2,b2)[0], param)
rel = np.abs(num-ana).max() / max(np.abs(ana).max(), 1e-12)
print(f"{name}: max abs diff {np.abs(num-ana).max():.3e} rel {rel:.3e}"
f" {'OK' if rel < 1e-7 else 'MISMATCH'}")3. All four should report relative error around \(10^{-9}\) or better. If one is wrong and the others are right, the error is almost always a missing transpose or a missing elementwise factor — the two patterns above.
4. Cost. The forward pass does two matrix-vector products: \(d_h d_{in} + d_{out}d_h\) multiply-adds. The backward pass does two matrix-vector products (\(\mW_2^\top\bar{\vz}_2\), and the analogous one if there were more layers) plus two outer products of the same total size. So the backward pass costs roughly two forward passes, independent of the number of parameters. This is the central fact that makes deep learning possible: you obtain all \(P\) partial derivatives for the price of a small constant number of forward evaluations.
Forward mode (the Dual class from Lecture 5) propagates one perturbation
direction at a time, so it needs \(P\) sweeps to get \(P\) partials — here \(P = d_hd_{in}+d_h+d_{out}d_h+
d_{out} = 35+7+21+3 = 66\), so about \(66\times\) the work. For a model with \(10^9\) parameters the ratio
becomes \(10^9\). Reverse mode wins whenever outputs are few and inputs are many; forward mode wins in the
opposite case (for example, sensitivity of many outputs to one hyperparameter). Both are the same chain rule,
associated in opposite orders.
Next: Part V. Everything so far has been deterministic. Real data is noisy, losses are expectations, and generalisation is a probabilistic claim. Lecture 13 builds random variables, expectation and variance from the ground up.