3.8 Matrix calculus and backpropagation

Training a neural network means computing how a single number (the loss) changes with millions of weights. The derivative of a composition is a product of matrices, and evaluating that product from the output backwards is backpropagation.

A neural network is a long composition of simple maps: multiply by a weight matrix, add a vector, squash each entry, repeat. Training it means nudging every weight to make the error smaller, which requires the derivative of one number (the loss) with respect to every weight. With millions of weights, computing each derivative separately would take millions of passes through the network. Backpropagation gets all of them in one backward pass, for about the cost of one forward pass. It is nothing but the chain rule, with the matrices multiplied in the right order.

Before you start

Chapter 1.4 (the derivative as best linear approximation, the chain rule), 2.5 (composition is matrix multiplication) and 2.10 (least squares).

Quick checkIn one variable, the chain rule says

(f \circ g)'(a) = f'(g(a))\,g'(a)

. What happens to the 'slopes'?

Quick checkIf

A

2 \times 3

and

B

3 \times 4

, what is the shape of

AB

The picture: derivatives as matrices

For a function $g: \mathbb R^n \to \mathbb R^m$ , the derivative at a point is the best linear approximation (1.4’s idea, one dimension up): $g(x + h) \approx g(x) + J_g(x)\,h$ , where $J_g(x)$ is the $m \times n$ Jacobian matrix whose entry $(i, j)$ is $\frac{\partial g_i}{\partial x_j}$ , the rate at which output $i$ changes when input $j$ is nudged. For a scalar function $f: \mathbb R^n \to \mathbb R$ the Jacobian is a single row, and its transpose, a column, is the gradient $\nabla f$ , pointing uphill. Chapters 5.2 and 5.3 make these definitions rigorous; here we use them to compute.

Where this came from· Backpropagation was discovered several times

The chain rule is from Leibniz (1676). Running it backwards through a computation (reverse-mode differentiation) was described by Seppo Linnainmaa in his 1970 master’s thesis. Paul Werbos applied it to neural networks in his 1974 PhD thesis. It became famous when David Rumelhart, Geoffrey Hinton and Ronald Williams showed in 1986 (Nature) that it lets multi-layer networks learn useful internal representations. Hinton later shared the 2018 Turing Award with Yoshua Bengio and Yann LeCun for work on deep learning. In 2024 he shared the Nobel Prize in Physics with John Hopfield for foundational discoveries on machine learning with artificial neural networks.

The rules you need

Three gradient formulas cover most of linear models. In each, $x \in \mathbb R^n$ is the variable and everything else is constant.

Function $f(x)$	Gradient $\nabla f(x)$	Why
$a^{\mathsf T}x = \sum_i a_i x_i$	$a$	$\frac{\partial}{\partial x_j}\sum_i a_ix_i = a_j$
$x^{\mathsf T}Ax$	$(A + A^{\mathsf T})x$ , which is $2Ax$ if $A$ is symmetric	expand $\sum_{i,j}a_{ij}x_ix_j$ and differentiate in $x_k$
$\lVert Ax - b \rVert^2$	$2A^{\mathsf T}(Ax - b)$	expand, then use the first two rows

Worked example

Derive the gradient of

f(x) = x^{\mathsf T}Ax

by differentiating in one coordinate.

Solution.

Expand into a double sum. $f(x) = \sum_{i=1}^n\sum_{j=1}^n a_{ij}x_ix_j$ . Row-times-column written out entry by entry.
Differentiate in $x_k$ . Only terms containing $x_k$ survive: those with $i = k$ give $\sum_j a_{kj}x_j$ , and those with $j = k$ give $\sum_i a_{ik}x_i$ . Product rule: $x_ix_j$ depends on $x_k$ through either factor; when $i = j = k$ both contributions appear, which matches $\frac{d}{dx}x^2 = 2x$ .
Recognise the sums. $\sum_j a_{kj}x_j = (Ax)_k$ and $\sum_i a_{ik}x_i = (A^{\mathsf T}x)_k$ . Entry $k$ of a matrix-vector product.

Answer. $\nabla f = Ax + A^{\mathsf T}x = (A + A^{\mathsf T})x$ . For symmetric $A$ this is $2Ax$ , the matrix version of $\frac{d}{dx}(ax^2) = 2ax$ .

Worked example

Recover the normal equations of least squares by setting a gradient to zero.

Solution.

Write the loss. $L(x) = \lVert Ax - b \rVert^2 = x^{\mathsf T}A^{\mathsf T}Ax - 2b^{\mathsf T}Ax + b^{\mathsf T}b$ . Expand $(Ax - b)^{\mathsf T}(Ax - b)$ ; the two middle terms are equal numbers.
Gradient term by term. $\nabla(x^{\mathsf T}A^{\mathsf T}Ax) = 2A^{\mathsf T}Ax$ (symmetric matrix $A^{\mathsf T}A$ ), $\nabla(-2b^{\mathsf T}Ax) = -2A^{\mathsf T}b$ (first rule with $a = A^{\mathsf T}b$ ), and the constant contributes 0. The rules in the table.
Set to zero. $2A^{\mathsf T}Ax - 2A^{\mathsf T}b = O$ , so $A^{\mathsf T}Ax = A^{\mathsf T}b$ . At a minimum of a smooth function, every directional slope is zero (Fermat’s theorem in many variables, 5.6).

Answer. The normal equations of 2.10, reached by calculus instead of geometry. Both routes agree: “the error is perpendicular to every column” is “the gradient is zero”.

Worked example

Find

\nabla f

for

f(x) = \lVert x \rVert^2 = x^{\mathsf T}x

, and the shape and value of the Jacobian of

g(x) = Wx

for

W

of size

3 \times 5

Solution.

Match a rule. $x^{\mathsf T}x = x^{\mathsf T}Ix$ , which is the second row of the table with $A = I$ . Writing the function in the form of a known rule is the first move for every gradient.
Apply it. $\nabla f = (I + I^{\mathsf T})x = 2x$ . $I$ is symmetric, so $(A + A^{\mathsf T})x = 2Ax$ .
Sanity check with one variable. For $n = 1$ , $f(x) = x^2$ and $f'(x) = 2x$ . The new rule must reduce to the old one.
Jacobian shape. $g: \mathbb R^5 \to \mathbb R^3$ , so $J_g$ is $3 \times 5$ . Rows are outputs, columns are inputs.
Jacobian value. Entry $(i, j)$ of $J_g$ is $\frac{\partial}{\partial x_j}\sum_k W_{ik}x_k = W_{ij}$ , so $J_g = W$ . A linear map is its own best linear approximation.

Answer. $\nabla\lVert x \rVert^2 = 2x$ , and $J_{Wx} = W$ (size $3 \times 5$ ).

Worked example

Find the minimiser of the ridge-regression loss

L(x) = \lVert Ax - b \rVert^2 + \lambda\lVert x \rVert^2

for

A = \begin{bmatrix}1&0\\0&1\\1&1\end{bmatrix}

b = (1, 2, 3)

\lambda = 1

Solution.

Gradient. $\nabla L = 2A^{\mathsf T}(Ax - b) + 2\lambda x$ . Sum of the third table row and the gradient of $\lambda x^{\mathsf T}x$ from the previous example.
Set to zero. $(A^{\mathsf T}A + \lambda I)x = A^{\mathsf T}b$ . Divide by 2 and collect the $x$ terms.
Compute the pieces. $A^{\mathsf T}A = \begin{bmatrix}2&1\\1&2\end{bmatrix}$ , so $A^{\mathsf T}A + I = \begin{bmatrix}3&1\\1&3\end{bmatrix}$ , and $A^{\mathsf T}b = (1 + 3,\ 2 + 3) = (4, 5)$ . Column dot products, then add $\lambda$ on the diagonal.
Solve. The determinant is $8$ , so $x = \frac18\begin{bmatrix}3&-1\\-1&3\end{bmatrix}\begin{bmatrix}4\\5\end{bmatrix} = \frac18(7, 11)$ . The inverse of a $2 \times 2$ matrix.

Answer. $x = (\frac78, \frac{11}8)$ . The matrix $A^{\mathsf T}A + \lambda I$ is positive definite for $\lambda > 0$ (3.2), so this is always solvable, even when $A$ has dependent columns.

Your answerFor

f(x, y) = 3x + 5y

(that is,

a^{\mathsf T}x

with

a = (3, 5)

), what is

\frac{\partial f}{\partial y}

Worked example

Prove

\nabla(a^{\mathsf T}x) = a

, then use the product rule to prove

\nabla(x^{\mathsf T}Ax) = (A + A^{\mathsf T})x

without expanding into sums.

Solution.

The linear rule. $a^{\mathsf T}x = \sum_ia_ix_i$ , so $\frac{\partial}{\partial x_k}(a^{\mathsf T}x) = a_k$ . Only the term with $i = k$ contains $x_k$ .
Split the quadratic. Write $x^{\mathsf T}Ax = u^{\mathsf T}v$ with $u = x$ and $v = Ax$ . Two vector-valued functions of $x$ multiplied by a dot product.
Product rule for a dot product. $\nabla(u^{\mathsf T}v) = J_u^{\mathsf T}v + J_v^{\mathsf T}u$ . Differentiate $\sum_iu_iv_i$ in $x_k$ : $\sum_i\big(\frac{\partial u_i}{\partial x_k}v_i + u_i\frac{\partial v_i}{\partial x_k}\big)$ , which is entry $k$ of $J_u^{\mathsf T}v + J_v^{\mathsf T}u$ .
Jacobians. $J_u = I$ and $J_v = A$ . $x$ is its own best linear approximation, and $Ax$ is linear with matrix $A$ .
Substitute. $\nabla(x^{\mathsf T}Ax) = I^{\mathsf T}Ax + A^{\mathsf T}x = Ax + A^{\mathsf T}x$ .

Answer. $\nabla(x^{\mathsf T}Ax) = (A + A^{\mathsf T})x$ . $\blacksquare$

Proof workshopRebuild the gradient of a quadratic form

Tap the lines in the order a proof would write them. Each line can only use what is already above it.

Your answerFor

f(x) = x^{\mathsf T}Ax

with

A = \begin{bmatrix}2&1\\1&3\end{bmatrix}

x = (1, 1)

, what is the first entry of

\nabla f

The chain rule as a product of matrices

If $x \mapsto u = g(x) \mapsto y = f(u)$ , then near a point

J_{f \circ g}(x) = J_f(u)\,J_g(x),

the product of the Jacobians, in the same order as the composition (Apostol 8.20, proved in 5.4). A network with layers $x \to h_1 \to h_2 \to \dots \to h_L \to \ell$ (a scalar loss) therefore has

\nabla_x \ell^{\mathsf T} = J_\ell\,J_{h_L}\cdots J_{h_2}\,J_{h_1},

a long product of matrices whose first factor is a single row.

Why· Why the order of multiplication matters so much

Matrix multiplication is associative, so any bracketing gives the same answer, but the cost differs enormously. Multiplying from the left starts with a row vector: (row) × (matrix) is a row, so every step is a cheap vector-times-matrix product. Multiplying from the right starts with two full matrices: every step is an expensive matrix-times-matrix product. For one output (the loss) and many inputs (the weights), going from the output backwards is far cheaper. That is the entire idea of backpropagation, also called reverse-mode differentiation.

Definition B

Computational graph

A computational graph records a computation as nodes (the inputs, every intermediate value, and the output) joined by arrows from the values each node uses. Evaluating it from inputs to output is the forward pass.

In one lineA diagram of a calculation: nodes are intermediate values, arrows say which values each one is computed from.

Why it existsIt turns any program into a composition of simple steps, so the chain rule can be applied mechanically, node by node.

Definition B

Backpropagation

Backpropagation computes $\frac{\partial \ell}{\partial v}$ for every node $v$ of a computational graph whose output $\ell$ is a scalar. It starts with $\frac{\partial \ell}{\partial \ell} = 1$ and moves backwards: each node’s derivative is the sum, over the nodes that use it, of (that node’s derivative) × (the local derivative of that node with respect to this one).

In one lineCompute the derivative of a scalar output with respect to every input by applying the chain rule backwards through the computational graph, reusing each intermediate result once.

Why it existsIt gets all n partial derivatives for about the cost of one forward pass, instead of n separate passes.

Worked example

Backpropagate through the network

z = w_1x

h = \max(0, z)

y = w_2h

\ell = \frac12(y - t)^2

x = 2

w_1 = 3

w_2 = -1

t = 1

, and take one gradient step.

Solution. Do a forward pass and keep every value. Then walk backwards, multiplying by one local derivative per arrow. The step-through shows each stage.

Answer. $\frac{\partial \ell}{\partial w_1} = 14$ and $\frac{\partial \ell}{\partial w_2} = -42$ . After one step with learning rate $0.01$ the weights are $w_1 = 2.86$ , $w_2 = -0.58$ and the loss falls from $24.5$ to about $9.3$ .

Worked example

Same network, but with input

x = -2

. What is

\frac{\partial \ell}{\partial w_1}

, and what has happened?

Solution.

Forward pass. $z = w_1x = -6$ , $h = \max(0, -6) = 0$ , $y = w_2h = 0$ , $\ell = \frac12(0 - 1)^2 = 0.5$ . The ReLU clips the negative value to 0.
Backward from the loss. $\frac{\partial \ell}{\partial y} = y - t = -1$ , and $\frac{\partial \ell}{\partial h} = -1\cdot w_2 = 1$ . Same local derivatives as before.
Through the ReLU. $\frac{\partial \ell}{\partial z} = \frac{\partial \ell}{\partial h}\cdot\mathbf 1[z > 0] = 1\cdot 0 = 0$ . $z < 0$ , so the ReLU has derivative 0 here.
To the weight. $\frac{\partial \ell}{\partial w_1} = 0\cdot x = 0$ . Everything upstream of a zero is zero.

Answer. $\frac{\partial \ell}{\partial w_1} = 0$ : the unit is “dead” for this input, so gradient descent cannot move $w_1$ using it. With the identity in place of ReLU the derivative would be $1\cdot(-2) = -2 \ne 0$ .

Worked example

A dense layer has

W = \begin{bmatrix}1&0&2\\3&-1&1\end{bmatrix}

, input

h = (1, 2, 3)

and incoming gradient

g = \frac{\partial \ell}{\partial z} = (1, -1)

. Find

\frac{\partial \ell}{\partial W}

and

\frac{\partial \ell}{\partial h}

Solution.

Shapes first. $z = Wh$ has 2 entries, so $g$ has 2. $\frac{\partial \ell}{\partial W}$ must be $2 \times 3$ and $\frac{\partial \ell}{\partial h}$ must have 3 entries. Each gradient has the shape of the thing it differentiates by.
Weight gradient. $\frac{\partial \ell}{\partial W} = gh^{\mathsf T} = \begin{bmatrix}1\\-1\end{bmatrix}\begin{bmatrix}1&2&3\end{bmatrix} = \begin{bmatrix}1&2&3\\-1&-2&-3\end{bmatrix}$ . Entry $(i, j)$ is $g_ih_j$ : weight $W_{ij}$ only affects output $i$ , by a factor $h_j$ .
Input gradient. $\frac{\partial \ell}{\partial h} = W^{\mathsf T}g = \begin{bmatrix}1&3\\0&-1\\2&1\end{bmatrix}\begin{bmatrix}1\\-1\end{bmatrix} = \begin{bmatrix}-2\\1\\1\end{bmatrix}$ . Input $h_j$ affects every output $i$ through $W_{ij}$ , so we add $g_iW_{ij}$ over $i$ , which is a multiplication by $W^{\mathsf T}$ .

Answer. $\frac{\partial \ell}{\partial W} = \begin{bmatrix}1&2&3\\-1&-2&-3\end{bmatrix}$ (rank 1, an outer product), $\frac{\partial \ell}{\partial h} = (-2, 1, 1)$ .

Worked example algorithms

A network has

L = 10

layers, each with an

n \times n

weight matrix and

n = 1000

. Compare the cost of the products

r^{\mathsf T}M_1\cdots M_L

evaluated left to right and right to left.

Solution.

Left to right. Each step multiplies a row by an $n \times n$ matrix, about $n^2 = 10^6$ operations. There are $L = 10$ steps: about $10^7$ . A row times a matrix needs $n$ multiplications for each of $n$ outputs.
Right to left. The first $L - 1 = 9$ steps multiply two $n \times n$ matrices, about $n^3 = 10^9$ operations each, then a final row product: about $9\times10^9$ . A matrix times a matrix needs $n$ multiplications for each of $n^2$ outputs.
Ratio. $\frac{9\times10^9}{10^7} = 900$ . Roughly $n$ times cheaper, as the general count $\frac{Ln^3}{Ln^2} = n$ says.

Answer. Left to right is about 900 times cheaper here, and about $n$ times cheaper in general. Reverse-mode differentiation is the left-to-right order, because it starts from the single-row end (the scalar loss).

Python labCheck backpropagation against finite differences, then train

import numpy as np

# the tiny network of the step-through: z = w1*x, h = max(0, z), y = w2*h, loss = (y - t)^2 / 2
def loss(w1, w2, x=2.0, t=1.0):
    h = max(0.0, w1 * x)
    return 0.5 * (w2 * h - t) ** 2

w1, w2, x, t = 3.0, -1.0, 2.0, 1.0
# backpropagation, by hand
z = w1 * x; h = max(0.0, z); y = w2 * h
dy = y - t
dw2 = dy * h
dh = dy * w2
dz = dh * (1.0 if z > 0 else 0.0)
dw1 = dz * x
print("backprop      :", dw1, dw2)

# finite differences: one extra forward pass per weight
eps = 1e-6
print("finite differ.:", (loss(w1 + eps, w2) - loss(w1, w2)) / eps,
                         (loss(w1, w2 + eps) - loss(w1, w2)) / eps)

# gradient descent: watch the loss fall
for step in range(5):
    print(step, round(loss(w1, w2), 4))
    z = w1 * x; h = max(0.0, z); y = w2 * h
    dy = y - t
    g1 = dy * w2 * (1.0 if z > 0 else 0.0) * x
    g2 = dy * h
    w1 -= 0.01 * g1; w2 -= 0.01 * g2

Try this. Run it. The backprop gradients $(14, -42)$ and the finite differences agree to about five digits, and the loss falls with every step. Change the input x in the first line of the loss function and in the hand computation to -2.0: the gradient with respect to w1 becomes 0, the dead-ReLU case from the example. Change both learning rates 0.01 to 0.1 in the training loop: the first step makes the loss rise from 24.5 to about 42.7 (the step overshoots), and then the unit goes dead and the loss freezes at 0.5.

Worked example

Check a backpropagated derivative numerically.

Solution.

Pick the derivative to check. Backprop gave $\frac{\partial \ell}{\partial w_1} = 14$ in the step-through. Every framework’s tests compare backprop against a finite difference.
Nudge the weight. With $w_1 = 3.001$ : $z = 6.002$ , $h = 6.002$ , $y = -6.002$ , $\ell = \frac12(7.002)^2 = 24.514002$ . One forward pass with the weight changed by $10^{-3}$ .
Form the difference quotient. $\frac{24.514002 - 24.5}{0.001} = 14.002$ . The definition of the derivative (1.4), with a small $h$ instead of a limit.

Answer. $14.002 \approx 14$ , matching backprop up to the $O(h)$ error of a one-sided difference. A numerical check costs one forward pass per weight, which is why it is only used for testing, never for training.

Your answerIn the step-through, change the target to

t = -6

(so the prediction

y = -6

is perfect). What is

\frac{\partial \ell}{\partial w_1}

Quick checkA layer computes

z = Wh

with

W

of size

4 \times 7

. What is the shape of

\frac{\partial \ell}{\partial W}

Quick checkFor a scalar function

f: \mathbb R^n \to \mathbb R

, the Jacobian

J_f

and the gradient

\nabla f

are related by:

Quick checkIn a network, a ReLU unit has

z < 0

for the current input. What gradient reaches the weights that feed it?

Quick checkWhy is backpropagation (reverse mode) used for training instead of computing each weight's derivative separately?

Proof workshopGive the reason for each backward step

Each line of the proof is shown in order. Choose the reason that justifies it, then press Check.

$\partial \ell/\partial y = y - t$
$\partial \ell/\partial h = (\partial \ell/\partial y)\, w_2$
$\partial \ell/\partial z = (\partial \ell/\partial h)\cdot 1[z > 0]$
$\partial \ell/\partial w_1 = (\partial \ell/\partial z)\, x$

Layers of vectors and matrices

For a dense layer $z = Wh + b$ with $W$ an $m \times n$ weight matrix, the same bookkeeping in matrix form gives, for the incoming gradient $g = \frac{\partial \ell}{\partial z}$ (an $m$ -vector):

\frac{\partial \ell}{\partial W} = g\,h^{\mathsf T}\ (m \times n),\qquad \frac{\partial \ell}{\partial b} = g,\qquad \frac{\partial \ell}{\partial h} = W^{\mathsf T}g.

The gradient flows backwards through a layer by multiplying by the transpose $W^{\mathsf T}$ . That is why the conditioning of the weight matrices (3.5) controls whether gradients explode (singular values above 1, multiplied many times) or vanish (below 1).

Your answerA layer has

W

of size

512 \times 2048

. How many numbers does

\frac{\partial \ell}{\partial W}

contain?

Where it is used

Used inside maths· Inside maths: optimisation in many variables

Gradient descent (5.10) moves against $\nabla f$ ; Newton’s method uses the Hessian (5.6); constrained problems add Lagrange multipliers (5.8). Each needs the derivatives this chapter computes. Chapter 5.4 proves the chain rule used here.

Used in the real world· AI: every deep learning framework

When you call loss.backward() in PyTorch, it walks the computational graph it recorded during the forward pass, in reverse, applying the rules above to every operation, and leaves $\frac{\partial \ell}{\partial w}$ in each weight’s .grad. Training a large language model is this, repeated over trillions of tokens. The cost of a backward pass is about twice that of the forward pass, which is why “training costs about three times inference” is a standard estimate.

Traps

Watch out· Shapes first

Before writing any gradient formula, write the shape of every quantity. $\frac{\partial \ell}{\partial W}$ always has the shape of $W$ . If your formula gives another shape, a transpose is missing.

Watch out· Gradient versus Jacobian

For a scalar $f$ , the gradient is a column and the Jacobian is a row: $\nabla f = J_f^{\mathsf T}$ . Mixing them up flips every product.

Problems

check

Find $\nabla f$ for $f(x) = \lVert x \rVert^2 = x^{\mathsf T}x$ .

Solution

$x^{\mathsf T}Ix$ with $I$ symmetric, so $\nabla f = 2x$ .

check

What is the shape of the Jacobian of $g(x) = Wx$ for $W$ of size $3 \times 5$ , and what is it?

Solution

$3 \times 5$ , and it equals $W$ : a linear map is its own best linear approximation.

use

Compute $\nabla L$ for the ridge-regression loss $L(x) = \lVert Ax - b \rVert^2 + \lambda\lVert x \rVert^2$ and solve $\nabla L = O$ .

Solution

$\nabla L = 2A^{\mathsf T}(Ax - b) + 2\lambda x$ . Setting it to zero: $(A^{\mathsf T}A + \lambda I)x = A^{\mathsf T}b$ . For $\lambda > 0$ the matrix is positive definite (3.2), hence invertible, even if $A$ has dependent columns.

use

In the step-through, replace ReLU by the identity ( $h = z$ ). Recompute $\frac{\partial \ell}{\partial w_1}$ .

Solution

The only change is the local derivative $\frac{\partial h}{\partial z} = 1$ , which was already 1 because $z > 0$ . So $\frac{\partial \ell}{\partial w_1} = 14$ again. (With $x = -2$ , ReLU would give 0 while the identity would not.)

prove

Prove $\nabla(a^{\mathsf T}x) = a$ and use it, with the product rule, to prove $\nabla(x^{\mathsf T}Ax) = (A + A^{\mathsf T})x$ without coordinates.

Solution

$\frac{\partial}{\partial x_k}\sum a_ix_i = a_k$ . For the quadratic: $x^{\mathsf T}Ax = x^{\mathsf T}(Ax)$ ; treat it as $u^{\mathsf T}v$ with $u = x$ , $v = Ax$ . The product rule gives $\nabla(u^{\mathsf T}v) = J_u^{\mathsf T}v + J_v^{\mathsf T}u = Ax + A^{\mathsf T}x$ , since $J_u = I$ and $J_v = A$ .

prove

For a dense layer $z = Wh$ , prove $\frac{\partial \ell}{\partial W_{ij}} = g_i h_j$ where $g = \frac{\partial \ell}{\partial z}$ , so $\frac{\partial \ell}{\partial W} = gh^{\mathsf T}$ .

Solution

$z_k = \sum_j W_{kj}h_j$ , so $\frac{\partial z_k}{\partial W_{ij}}$ is $h_j$ if $k = i$ and 0 otherwise. By the chain rule, $\frac{\partial \ell}{\partial W_{ij}} = \sum_k g_k\frac{\partial z_k}{\partial W_{ij}} = g_ih_j$ , which is entry $(i, j)$ of $gh^{\mathsf T}$ .

stretchalgorithms

Multiplying a $1 \times n$ row by an $n \times n$ matrix costs about $n^2$ operations; an $n \times n$ by $n \times n$ product costs about $n^3$ . For a chain $r^{\mathsf T}M_1M_2\cdots M_L$ of $L$ such matrices, compare the cost of evaluating left-to-right with right-to-left.

Solution

Left-to-right: each step is row × matrix, so $L$ steps of $n^2$ : about $Ln^2$ . Right-to-left: $L - 1$ matrix products of $n^3$ each, then one row product: about $Ln^3$ . For $n = 1000$ that is a factor of 1000. Reverse-mode differentiation is the left-to-right order, starting from the scalar loss.

Summary

The derivative of $g: \mathbb R^n \to \mathbb R^m$ is the $m \times n$ Jacobian; for a scalar function, the gradient is its transpose.
Rules: $\nabla(a^{\mathsf T}x) = a$ , $\nabla(x^{\mathsf T}Ax) = (A + A^{\mathsf T})x$ , $\nabla\lVert Ax - b \rVert^2 = 2A^{\mathsf T}(Ax - b)$ . Setting the last to zero gives the normal equations.
Chain rule: the Jacobian of a composition is the product of Jacobians.
Backpropagation multiplies that product starting from the scalar loss, so all gradients cost about one forward pass. Through a layer $z = Wh$ : $\frac{\partial \ell}{\partial W} = gh^{\mathsf T}$ , $\frac{\partial \ell}{\partial h} = W^{\mathsf T}g$ .
Check gradients numerically with a difference quotient; never train that way.

Watch and read

The same ideas explained by someone else. Use these when an explanation here does not land.

3.8 Matrix calculus and backpropagation

The picture: derivatives as matrices

The rules you need

The chain rule as a product of matrices

Computational graph

Backpropagation

Layers of vectors and matrices

Where it is used

Traps

Problems

Summary

Watch and read

Where this sits

Comes from

Leads to

Same idea elsewhere