13 Matrix Calculus
Matrix calculus is not new mathematics. It is Chapter 12 done in bulk — the same partial derivatives, arranged so you can manipulate them without writing indices.
The payoff is that results which take a page of subscripts collapse to a line. The least squares solution, the gradient of a neural network layer, and the Fisher information matrix are all one or two steps once you have the notation. Getting there requires settling one tedious question first: what shape is a derivative?
13.1 Why matrix calculus
Suppose you want to minimize the least squares loss
\[ L(\boldsymbol{\beta}) = \|\mathbf{y} - \mathbf{X}\boldsymbol{\beta}\|^2 \]
Written out with indices, that is \(\sum_i (y_i - \sum_j x_{ij}\beta_j)^2\), and differentiating with respect to \(\beta_k\) means chasing two nested sums and a Kronecker delta. It works. It is miserable, and it does not scale to anything more complicated.
With matrix calculus:
\[ \nabla_{\boldsymbol{\beta}} L = -2\mathbf{X}^\top(\mathbf{y} - \mathbf{X}\boldsymbol{\beta}) \]
Set it to zero and the normal equations fall out (Section 13.7). Three lines instead of a page, and the result is in a form you can compute with directly.
13.2 Layout conventions
Here is the thing that trips everyone up. If \(f\) is a scalar and \(\mathbf{x}\) is a vector of length \(n\), then \(\partial f/\partial\mathbf{x}\) has \(n\) components — but is it a row or a column?
Both answers are in print. Numerator layout makes it a row; denominator layout makes it a column. They differ by a transpose, which means every identity you look up may or may not match your convention, and a mismatch produces shape errors that look like bugs.
This book uses the convention that dominates machine learning:
A gradient has the same shape as the variable it differentiates with respect to.
So for scalar \(f\) and \(\mathbf{x} \in \mathbb{R}^n\), the gradient \(\nabla_{\mathbf{x}}f\) is a column vector in \(\mathbb{R}^n\). For scalar \(f\) and matrix \(\mathbf{X} \in \mathbb{R}^{m \times n}\), the derivative \(\partial f/\partial\mathbf{X}\) is an \(m \times n\) matrix.
A Jacobian is \(m \times n\): for \(\mathbf{y} \in \mathbb{R}^m\) and \(\mathbf{x} \in \mathbb{R}^n\), the matrix \(\mathbf{J}\) has \(J_{ij} = \partial y_i/\partial x_j\).
This is a mixed convention — gradients follow the denominator, Jacobians the numerator — and it is chosen because it makes the two things you actually do work cleanly. Gradient descent needs \(\boldsymbol{\theta} \leftarrow \boldsymbol{\theta} - \eta\nabla L\), which requires the gradient to match \(\boldsymbol{\theta}\)’s shape. Backpropagation needs Jacobians to chain by ordinary matrix multiplication.
When you look up an identity, check the source’s convention before using it. If a result comes out with the wrong shape, a transpose is the most likely cause, not an error in the algebra.
The reliable sanity check: differentiate a small case numerically and compare. Every identity in this chapter is verified that way below.
13.3 Scalar with respect to vector
The gradient collects every partial derivative into a column:
\[ \nabla_{\mathbf{x}}f = \begin{bmatrix} \partial f/\partial x_1 \\ \vdots \\ \partial f/\partial x_n \end{bmatrix} \]
Three results carry most of the weight.
Linear form. For \(f(\mathbf{x}) = \mathbf{a}^\top\mathbf{x} = \sum_i a_ix_i\), the \(i\)-th partial is just \(a_i\):
\[ \nabla_{\mathbf{x}}(\mathbf{a}^\top\mathbf{x}) = \mathbf{a} \]
The matrix analogue of \(\frac{d}{dx}(ax) = a\).
Squared norm. For \(f(\mathbf{x}) = \mathbf{x}^\top\mathbf{x} = \sum_i x_i^2\):
\[ \nabla_{\mathbf{x}}(\mathbf{x}^\top\mathbf{x}) = 2\mathbf{x} \]
The analogue of \(\frac{d}{dx}(x^2) = 2x\).
Quadratic form. For \(f(\mathbf{x}) = \mathbf{x}^\top\mathbf{A}\mathbf{x}\):
\[ \nabla_{\mathbf{x}}(\mathbf{x}^\top\mathbf{A}\mathbf{x}) = (\mathbf{A} + \mathbf{A}^\top)\mathbf{x} \tag{13.1}\]
and when \(\mathbf{A}\) is symmetric this simplifies to \(2\mathbf{A}\mathbf{x}\) — but only then.
[,1] [,2]
numeric 3 -1
formula 3 -1 [,1] [,2]
numeric 2 4
formula 2 4Equation 13.1 is the identity people get wrong, because they remember the symmetric shortcut and apply it generally.
[,1] [,2]
numeric 6 13
correct 6 13
wrong 8 12\(2\mathbf{A}\mathbf{x}\) gives \((8, 12)\) where the truth is \((6, 13)\). In practice most matrices you differentiate through are symmetric — covariances, Gram matrices, Hessians — which is exactly why the mistake survives so long before it bites.
13.4 Vector with respect to vector
When both input and output are vectors, the derivative is the Jacobian (Section 12.7): an \(m \times n\) matrix whose \((i,j)\) entry is \(\partial y_i/\partial x_j\).
\[ \mathbf{J} = \frac{\partial\mathbf{y}}{\partial\mathbf{x}} = \begin{bmatrix} \partial y_1/\partial x_1 & \cdots & \partial y_1/\partial x_n \\ \vdots & \ddots & \vdots \\ \partial y_m/\partial x_1 & \cdots & \partial y_m/\partial x_n \end{bmatrix} \]
The fundamental case is a linear map:
\[ \mathbf{y} = \mathbf{A}\mathbf{x} \quad\Longrightarrow\quad \frac{\partial\mathbf{y}}{\partial\mathbf{x}} = \mathbf{A} \]
which is what you would hope: the derivative of a linear map is the map itself, exactly as \(\frac{d}{dx}(ax) = a\).
Row \(i\) of a Jacobian is the gradient of output \(i\). So a gradient is the special case \(m = 1\) — up to the transpose our convention introduces, which is the price of wanting gradients to match the variable’s shape.
An elementwise function has a diagonal Jacobian, since output \(i\) depends only on input \(i\). Applying \(\sigma\) componentwise:
\[ \frac{\partial\sigma(\mathbf{x})}{\partial\mathbf{x}} = \operatorname{diag}\bigl(\sigma'(x_1), \dots, \sigma'(x_n)\bigr) \]
That diagonal structure is why activation functions are cheap to backpropagate through: multiplying by a diagonal matrix is an elementwise product, not a matrix multiply.
13.5 Scalar with respect to matrix
When \(f\) is a scalar and \(\mathbf{X}\) a matrix, \(\partial f/\partial\mathbf{X}\) is a matrix of the same shape, with entry \((i,j)\) equal to \(\partial f/\partial x_{ij}\).
The identities worth knowing:
\[ \begin{aligned} \frac{\partial}{\partial\mathbf{X}}\operatorname{tr}(\mathbf{X}) &= \mathbf{I} \\ \frac{\partial}{\partial\mathbf{X}}\operatorname{tr}(\mathbf{A}\mathbf{X}) &= \mathbf{A}^\top \\ \frac{\partial}{\partial\mathbf{X}}(\mathbf{a}^\top\mathbf{X}\mathbf{b}) &= \mathbf{a}\mathbf{b}^\top \\ \frac{\partial}{\partial\mathbf{X}}\log|\mathbf{X}| &= \mathbf{X}^{-\top} \end{aligned} \]
The third is the one that appears constantly in neural networks: the gradient with respect to a weight matrix is an outer product of the incoming activation and the outgoing error. The last shows up in Gaussian likelihoods, where \(\log|\boldsymbol{\Sigma}|\) is the normalizing term.
13.6 An identities cheat sheet
Everything above in one table. \(\mathbf{a}\), \(\mathbf{b}\) are constant vectors; \(\mathbf{A}\) a constant matrix.
| Expression | Derivative | With respect to |
|---|---|---|
| \(\mathbf{a}^\top\mathbf{x}\) | \(\mathbf{a}\) | \(\mathbf{x}\) |
| \(\mathbf{x}^\top\mathbf{a}\) | \(\mathbf{a}\) | \(\mathbf{x}\) |
| \(\mathbf{x}^\top\mathbf{x}\) | \(2\mathbf{x}\) | \(\mathbf{x}\) |
| \(\mathbf{x}^\top\mathbf{A}\mathbf{x}\) | \((\mathbf{A}+\mathbf{A}^\top)\mathbf{x}\) | \(\mathbf{x}\) |
| \(\mathbf{x}^\top\mathbf{A}\mathbf{x}\), \(\mathbf{A}\) symmetric | \(2\mathbf{A}\mathbf{x}\) | \(\mathbf{x}\) |
| \(\mathbf{A}\mathbf{x}\) | \(\mathbf{A}\) | \(\mathbf{x}\) |
| \(\|\mathbf{A}\mathbf{x}-\mathbf{b}\|^2\) | \(2\mathbf{A}^\top(\mathbf{A}\mathbf{x}-\mathbf{b})\) | \(\mathbf{x}\) |
| \(\operatorname{tr}(\mathbf{A}\mathbf{X})\) | \(\mathbf{A}^\top\) | \(\mathbf{X}\) |
| \(\mathbf{a}^\top\mathbf{X}\mathbf{b}\) | \(\mathbf{a}\mathbf{b}^\top\) | \(\mathbf{X}\) |
| \(\log\lvert\mathbf{X}\rvert\) | \(\mathbf{X}^{-\top}\) | \(\mathbf{X}\) |
The seventh row is the one to memorize. It is least squares, ridge, and the linear part of every neural network layer.
13.7 Deriving the least squares solution
Now the payoff. Minimize
\[ L(\boldsymbol{\beta}) = \|\mathbf{y} - \mathbf{X}\boldsymbol{\beta}\|^2 \]
Step 1 — expand. Using \(\|\mathbf{v}\|^2 = \mathbf{v}^\top\mathbf{v}\):
\[ \begin{aligned} L &= (\mathbf{y} - \mathbf{X}\boldsymbol{\beta})^\top(\mathbf{y} - \mathbf{X}\boldsymbol{\beta}) \\ &= \mathbf{y}^\top\mathbf{y} - 2\boldsymbol{\beta}^\top\mathbf{X}^\top\mathbf{y} + \boldsymbol{\beta}^\top\mathbf{X}^\top\mathbf{X}\boldsymbol{\beta} \end{aligned} \]
The two cross terms combine because \(\boldsymbol{\beta}^\top\mathbf{X}^\top\mathbf{y}\) is a scalar, and a scalar equals its own transpose.
Step 2 — differentiate, term by term. The first term has no \(\boldsymbol{\beta}\), so it vanishes. The second is a linear form. The third is a quadratic form with the symmetric matrix \(\mathbf{X}^\top\mathbf{X}\) (symmetric for any \(\mathbf{X}\), as shown in Section 5.16), so the shortcut applies:
\[ \nabla_{\boldsymbol{\beta}}L = -2\mathbf{X}^\top\mathbf{y} + 2\mathbf{X}^\top\mathbf{X}\boldsymbol{\beta} \]
Step 3 — set to zero.
\[ \mathbf{X}^\top\mathbf{X}\hat{\boldsymbol{\beta}} = \mathbf{X}^\top\mathbf{y} \tag{13.2}\]
The normal equations (Section 6.9), derived in three steps rather than found geometrically. Both routes reach the same place, which is a good sign for both.
Step 4 — confirm it is a minimum. The Hessian is
\[ \nabla^2_{\boldsymbol{\beta}}L = 2\mathbf{X}^\top\mathbf{X} \]
which is positive definite whenever \(\mathbf{X}\) has independent columns (Section 8.9). Positive definite Hessian, so the critical point is a minimum — and since the loss is a convex quadratic, the only one.
[,1] [,2]
numeric 2.081668e-11 -1.249001e-10
formula 1.776357e-15 6.217249e-15Both zero, to rounding. The Hessian’s eigenvalues confirm the minimum:
Both positive, so positive definite.
bs <- seq(1.5, 2.5, length.out = 200)
Ls <- sapply(bs, function(b) loss(c(beta[1], b)))
# draw_line, not draw_plane: slope and squared loss have
# unrelated units, and horizontality survives unequal
# axis scaling anyway
draw_line(
x = bs,
y = list(
`L(beta)` = Ls,
`tangent at the minimum` = rep(loss(beta), length(bs))
),
points = FALSE,
xlab = "slope",
ylab = "loss"
)13.8 Backpropagation is the chain rule
A neural network is a composition (Section 3.7):
\[ L = \ell \circ f_K \circ f_{K-1} \circ \cdots \circ f_1 \]
The chain rule (Section 11.6) in matrix form says the derivative of a composition is the product of the Jacobians:
\[ \frac{\partial L}{\partial\mathbf{x}} = \frac{\partial L}{\partial\mathbf{h}_K} \frac{\partial\mathbf{h}_K}{\partial\mathbf{h}_{K-1}} \cdots \frac{\partial\mathbf{h}_1}{\partial\mathbf{x}} \tag{13.3}\]
That is the entire theory. Everything else is a question about the order in which you multiply.
Matrix multiplication is associative (Section 5.9), so you may evaluate Equation 13.3 from either end. The two choices have names:
| Order | Name | Good when |
|---|---|---|
| right to left | forward mode | few inputs, many outputs |
| left to right | reverse mode | many inputs, few outputs |
Reverse mode is backpropagation, and it is the right choice for machine learning because the loss is a single scalar. Starting from the left, \(\partial L/\partial \mathbf{h}_K\) is a row vector, and every subsequent product is vector-times-matrix — never matrix-times-matrix. Going the other way builds full Jacobians at every step.
This is exactly the lesson from Section 5.14, where grouping \(\mathbf{A}(\mathbf{B}\mathbf{C})\) instead of \((\mathbf{A}\mathbf{B})\mathbf{C}\) saved a factor of 500. Backpropagation is that same associativity choice, made once, for the whole of deep learning.
For a network with a million parameters, reverse mode costs roughly one forward pass. Forward mode would cost a million of them.
This asymmetry is why deep learning happened. Training needs \(\partial L/\partial \boldsymbol{\theta}\) for millions of parameters and one scalar loss — the exact case reverse mode is best at, by a factor equal to the parameter count. Had the useful case been many outputs and few inputs, the field would look entirely different.
13.9 Automatic differentiation
There are three ways to get a derivative on a computer, and they are often confused.
| Method | How | Cost | Error |
|---|---|---|---|
| Numerical | finite differences | one evaluation per input | truncation + rounding |
| Symbolic | manipulate the formula algebraically | expression can blow up | exact |
| Automatic | chain rule applied to the executed operations | ~1 evaluation (reverse mode) | exact to machine precision |
Automatic differentiation is neither of the other two. It does not approximate, and it does not produce a formula. It records the elementary operations a program actually performs — each with a known derivative — and applies Equation 13.3 to that record.
The consequences are worth being precise about:
-
No step size. Nothing to tune, and none of the U-shaped error curve from
- Accuracy is machine precision.
- Works on programs, not just formulas. Loops, branches and function calls are all fine, because it differentiates the trace of what ran.
- Costs about one function evaluation in reverse mode, regardless of parameter count.
R has autodiff through packages such as torch and RTMB; the calculations in this book are small enough to differentiate by hand, so we do that instead. But every gradient computed by a deep learning framework is Equation 13.3, evaluated left to right.
13.10 Summary
| Rule | Result |
|---|---|
| Convention | gradients match the variable’s shape; Jacobians are \(m \times n\) |
| \(\mathbf{a}^\top\mathbf{x}\) | \(\mathbf{a}\) |
| \(\mathbf{x}^\top\mathbf{A}\mathbf{x}\) | \((\mathbf{A}+\mathbf{A}^\top)\mathbf{x}\) |
| \(\|\mathbf{A}\mathbf{x}-\mathbf{b}\|^2\) | \(2\mathbf{A}^\top(\mathbf{A}\mathbf{x}-\mathbf{b})\) |
| Chain rule | product of Jacobians |
| Backpropagation | that product, evaluated left to right |
| Always | check a small case numerically |
13.11 Exercises
1. Find \(\nabla_{\mathbf{x}}(\mathbf{a}^\top\mathbf{x} + \mathbf{x}^\top\mathbf{x})\) and verify numerically.
2. For symmetric \(\mathbf{S}\), confirm that \(\nabla(\mathbf{x}^\top\mathbf{S}\mathbf{x}) = 2\mathbf{S}\mathbf{x}\), and check the shortcut fails for a non-symmetric matrix.
[,1] [,2]
numeric 8 10
shortcut 8 10They agree. For the non-symmetric A from earlier they do not, as Equation 13.1 requires — the general formula \((\mathbf{A}+\mathbf{A}^\top)
\mathbf{x}\) reduces to \(2\mathbf{A}\mathbf{x}\) only when \(\mathbf{A} =
\mathbf{A}^\top\).
3. Derive the gradient of the ridge objective \(\|\mathbf{y}-\mathbf{X}\boldsymbol{\beta}\|^2 + \lambda\|\boldsymbol{\beta}\|^2\) and solve for \(\hat{\boldsymbol{\beta}}\).
The penalty contributes \(2\lambda\boldsymbol{\beta}\), so
\[ \nabla L = -2\mathbf{X}^\top(\mathbf{y}-\mathbf{X}\boldsymbol{\beta}) + 2\lambda\boldsymbol{\beta} = \mathbf{0} \]
\[ (\mathbf{X}^\top\mathbf{X} + \lambda\mathbf{I})\hat{\boldsymbol{\beta}} = \mathbf{X}^\top\mathbf{y} \]
exactly Equation 6.3.
[,1] [,2]
numeric -2.220446e-10 4.440892e-10
should_be 0.000000e+00 0.000000e+00The gradient vanishes at the ridge solution, confirming the derivation. Note this version penalizes the intercept, which real implementations usually do not.
4. Show that \(\nabla_{\mathbf{x}}\|\mathbf{A}\mathbf{x}-\mathbf{b}\|^2 = 2\mathbf{A}^\top(\mathbf{A}\mathbf{x}-\mathbf{b})\) using the chain rule.
Let \(\mathbf{r} = \mathbf{A}\mathbf{x}-\mathbf{b}\), so \(L = \mathbf{r}^\top\mathbf{r}\). Then \(\partial L/\partial\mathbf{r} = 2\mathbf{r}\) and \(\partial\mathbf{r}/\partial\mathbf{x} = \mathbf{A}\). Chaining, and transposing to keep the gradient a column:
\[ \nabla_{\mathbf{x}}L = \mathbf{A}^\top(2\mathbf{r}) = 2\mathbf{A}^\top(\mathbf{A}\mathbf{x}-\mathbf{b}) \]
[,1] [,2]
numeric 56 198
formula 56 198The transpose is where the layout convention shows up: without it the result would be a row and would not match \(\mathbf{x}\).
5. A chain has Jacobians of shapes \(1\times100\), \(100\times100\) and \(100\times1000\). Count the multiply-adds for left-to-right versus right-to-left.
Multiplying \(m \times n\) by \(n \times p\) costs \(mnp\).
Left to right (reverse mode): \((1\times100)(100\times100)\) costs \(10{,}000\) and gives \(1\times100\); then \((1\times100)(100\times1000)\) costs \(100{,}000\). Total 110,000.
Right to left (forward mode): \((100\times100)(100\times1000)\) costs \(10{,}000{,}000\) and gives \(100\times1000\); then \((1\times100)(100\times1000)\) costs \(100{,}000\). Total 10,100,000.
c(
reverse = 1 * 100 * 100 + 1 * 100 * 1000,
forward = 100 * 100 * 1000 + 1 * 100 * 1000
) reverse forward
110000 10100000 Nearly a hundredfold, from the same product evaluated in a different order. With a scalar at the left end, reverse mode never forms a matrix.
6. Verify numerically that \(\partial(\mathbf{a}^\top\mathbf{X}\mathbf{b})/\partial\mathbf{X} = \mathbf{a}\mathbf{b}^\top\).
av <- c(1, 2)
bv <- c(3, -1)
Xm <- matrix(c(1, 2, 0, 4), nrow = 2, byrow = TRUE)
fmat <- function(M) drop(crossprod(av, M %*% bv))
hh <- 1e-6
numeric_grad <- outer(1:2, 1:2, Vectorize(function(i, j) {
E <- matrix(0, 2, 2)
E[i, j] <- hh
(fmat(Xm + E) - fmat(Xm - E)) / (2 * hh)
}))
round(numeric_grad, 6) [,1] [,2]
[1,] 3 -1
[2,] 6 -2outer(av, bv) [,1] [,2]
[1,] 3 -1
[2,] 6 -2The derivative is an outer product, and it does not depend on \(\mathbf{X}\) at all — \(\mathbf{a}^\top\mathbf{X}\mathbf{b}\) is linear in \(\mathbf{X}\). This is the identity behind every weight-matrix gradient in a neural network.