Neural networks do not have one parameter — they have millions, arranged in matrices. To derive gradients efficiently we need calculus that speaks the language of vectors and matrices. Matrix calculus looks intimidating, but it reduces to a handful of identities and one golden rule: check the shapes.
Layout convention#
In this course we use the denominator layout for gradients of scalars: if $L$ is a scalar and $\mathbf{W}$ is an $m \times n$ matrix, then $\partial L/\partial \mathbf{W}$ is also $m \times n$, with entry $(i,j)$ equal to $\partial L/\partial W_{ij}$. This matches what PyTorch's .grad stores: the gradient has the same shape as the parameter.
The Jacobian#
For a vector function $\mathbf{f}: \mathbb{R}^n \to \mathbb{R}^m$, the Jacobian is the $m \times n$ matrix of all first partial derivatives:
It is the best linear approximation of $\mathbf{f}$: $\mathbf{f}(\mathbf{x} + \boldsymbol{\delta}) \approx \mathbf{f}(\mathbf{x}) + \mathbf{J}\boldsymbol{\delta}$. The chain rule for vector functions becomes matrix multiplication of Jacobians:
Important Jacobians:
- Linear map $\mathbf{y} = \mathbf{W}\mathbf{x}$: $\partial\mathbf{y}/\partial\mathbf{x} = \mathbf{W}$.
- Element-wise activation $\mathbf{y} = \phi(\mathbf{x})$: $\partial\mathbf{y}/\partial\mathbf{x} = \text{diag}(\phi'(\mathbf{x}))$ — diagonal, so we implement it as an element-wise product, never as a full matrix.
- Softmax $\mathbf{s} = \text{softmax}(\mathbf{z})$: $\partial s_i/\partial z_j = s_i(\delta_{ij} - s_j)$, i.e. $\text{diag}(\mathbf{s}) - \mathbf{s}\mathbf{s}^\top$.
The Hessian#
For scalar $f: \mathbb{R}^n \to \mathbb{R}$, the Hessian $\mathbf{H}_{ij} = \partial^2 f/\partial x_i\partial x_j$ is symmetric (for smooth $f$). It describes curvature: positive definite at a point with zero gradient means a local minimum; indefinite means a saddle point. In high-dimensional deep learning loss surfaces, saddle points vastly outnumber poor local minima.
Identities to memorise#
| Scalar function $f$ | Gradient $\nabla_{\mathbf{x}} f$ |
|---|---|
| $\mathbf{a}^\top\mathbf{x}$ | $\mathbf{a}$ |
| $\mathbf{x}^\top\mathbf{x}$ | $2\mathbf{x}$ |
| $\mathbf{x}^\top\mathbf{A}\mathbf{x}$ | $(\mathbf{A} + \mathbf{A}^\top)\mathbf{x}$ (equals $2\mathbf{A}\mathbf{x}$ if symmetric) |
| $\|\mathbf{A}\mathbf{x} - \mathbf{b}\|^2$ | $2\mathbf{A}^\top(\mathbf{A}\mathbf{x} - \mathbf{b})$ |
| Scalar function of a matrix | Gradient w.r.t. $\mathbf{W}$ |
|---|---|
| $\mathbf{a}^\top\mathbf{W}\mathbf{b}$ | $\mathbf{a}\mathbf{b}^\top$ |
| $\text{tr}(\mathbf{A}\mathbf{W})$ | $\mathbf{A}^\top$ |
| $\|\mathbf{W}\|_F^2$ | $2\mathbf{W}$ |
Deriving linear regression in one line#
Minimise $L(\mathbf{w}) = \|\mathbf{X}\mathbf{w} - \mathbf{y}\|^2$. Using the table:
These are the normal equations. The Hessian is $2\mathbf{X}^\top\mathbf{X}$, which is positive semi-definite, so the loss is convex and any solution is a global minimum.
The backprop rule for a linear layer#
Consider a layer $\mathbf{Y} = \mathbf{X}\mathbf{W}$ with a batch $\mathbf{X} \in \mathbb{R}^{B \times d}$, weights $\mathbf{W} \in \mathbb{R}^{d \times k}$, output $\mathbf{Y} \in \mathbb{R}^{B \times k}$. Suppose we already know the upstream gradient $\mathbf{G} = \partial L/\partial \mathbf{Y} \in \mathbb{R}^{B \times k}$. Then
import numpy as np
rng = np.random.default_rng(0)
B, d, k = 4, 3, 2
X, W = rng.normal(size=(B, d)), rng.normal(size=(d, k))
T = rng.normal(size=(B, k)) # targets
def L(W): # L = 0.5 * ||XW - T||_F^2
return 0.5 * np.sum((X @ W - T) ** 2)
G = X @ W - T # dL/dY
grad_W = X.T @ G # our formula
num = np.zeros_like(W); eps = 1e-6
for i in range(d):
for j in range(k):
E = np.zeros_like(W); E[i, j] = eps
num[i, j] = (L(W + E) - L(W - E)) / (2 * eps)
print(np.allclose(grad_W, num, atol=1e-6)) # TrueSoftmax + cross-entropy#
For logits $\mathbf{z}$, probabilities $\mathbf{p} = \text{softmax}(\mathbf{z})$ and one-hot target $\mathbf{y}$, the loss $L = -\sum_i y_i\ln p_i$ has the elegant gradient
— the same "prediction minus target" pattern we found for the sigmoid. Frameworks fuse softmax and cross-entropy into one operation precisely to use this simple, numerically stable gradient.
Vector–Jacobian products#
Frameworks never build full Jacobians — for a layer with a million inputs and outputs that would be $10^{12}$ numbers. Reverse-mode automatic differentiation computes vector–Jacobian products $\mathbf{g}^\top\mathbf{J}$ directly, which cost about as much as the forward computation. This is why training costs only a small constant factor more than inference.