25 - Matrix Calculus and Automatic Differentiation: From the Jacobian to an Autograd Written From Scratch

Machine Learning Engineer course · Module 02 · Mathematics and statistics · Lesson 10

Goal. Differentiate a function with respect to a vector and with respect to a matrix without getting the shapes wrong, and understand how a library manages to compute the gradients of millions of parameters on its own. By the end you can declare a layout convention and stick to it, you can derive for yourself the five identities that recur in every network, you can tell forward mode from reverse mode on the basis of cost, and you can write a hundred-and-twenty-line autograd whose gradients agree with centred finite differences to within a relative error of 10⁻⁷.

Myth to dispel. "Differentiating with respect to a matrix is something exotic, a job for mathematicians." No: it is the same partial derivative as always, repeated for every cell and arranged back into a table of the same shape as the matrix. The only real difficulty is deciding how to stack the results, and 90% of sign and shape errors come from never having declared that decision.

Introduction

In the lesson on the gradient and gradient descent you derived the chain rule on a single neuron: one number in, one number out, one derivative. In the lesson on vectors, matrices and systems you learnt to multiply matrices and to read shapes. These two things put together are still not enough: in the lesson on backpropagation you will need ∂L/∂W with W of shape (128, 64), that is, the derivative of a single number with respect to 8192 numbers at once. This lesson is the bridge.

It is also the lesson that settles a debt: so far we have used ŷ − y without deriving it and have called "autograd" a black box. By the end there is nothing black left inside that box.

Final result of this lesson. A hundred-and-twenty-line file that builds a computational graph, walks it backwards exactly once and produces the gradients of a 2→16→16→2 network, verified against two independent oracles: the closed-form formulas derived by hand, and centred finite differences.

Roadmap of the lesson.

1. Differentiating with respect to a vector, and why the layout convention must be declared.
2. The five identities that actually matter, derived in full, component by component.
3. The Jacobian, the vector-Jacobian product and the memory arithmetic.
4. The Hessian, curvature, Newton's step and why deep learning does not use it.
5. Forward mode and reverse mode: the same result, two different costs.
6. A reverse-mode autograd written from scratch.
7. Gradient checking done properly: centred differences, the choice of h, relative error.
8. The experiment: a two-armed spiral, trained with our own autograd.

The code in this lesson was run with Python 3.14.4 and NumPy 2.3.5, and uses only NumPy: no deep learning library, because the whole point is to build its central piece. Where random data appear, the generator is initialised with a fixed seed, so the outputs are reproducible line by line; with different versions of NumPy the last digits may change, the conclusions do not.

Differentiating with respect to a vector

The analogy: the control panel with many knobs

Imagine a control unit with 8192 knobs and a single dial showing the loss. The partial derivative with respect to one knob is how far the needle moves per millimetre of that knob, holding all the others still.

The gradient is not one of these readings: it is the simultaneous reading of all 8192, noted next to each knob. And the point is that this sheet has the shape of the control panel: if the knobs sit in a 128 by 64 grid, so do the notes. This is not cosmetics: it is what lets you write the update as W ← W − η · ∂L/∂W, a subtraction between two tables of the same shape.

Definition. Given a real-valued function f(x) with x ∈ Rⁿ, the gradient is the vector of partial derivatives

∇f(x) = ∂f/∂x = (∂f/∂x₁, ∂f/∂x₂, …, ∂f/∂xₙ).

Read it as a sentence: the i-th component says how much f grows if you move coordinate i alone by a little, leaving the others still.

The layout convention, declared once and for all

Here begins the part that books take for granted and that costs hours of debugging. When you differentiate a scalar with respect to a vector you get n numbers, and you must decide whether to arrange them in a column or in a row. The two choices have names.

In this lesson, and throughout the course, we adopt the denominator convention for gradients: ∂f/∂x always has the shape of x. The reason is practical: the update is a subtraction between a parameter and its gradient, and if the two shapes do not match you have to transpose at every step — and sooner or later you forget.

For Jacobians we use the numerator convention instead: J has shape (m, n), one row per output and one column per input, the shape in which the chain rule is written as a matrix product in the natural order. The two conventions coexist because they concern different objects: the gradient of a scalar and the Jacobian of a vector-valued function. Declaring it at the top of the file is worth more than any later comment.

The transposed gradient that almost works

Here is the costliest mistake of the lesson, and it is worth seeing before the formulas. You compute the gradient with respect to a weight matrix and write it transposed. If the matrix is rectangular, NumPy raises a shape error at the first subtraction and in thirty seconds you have found the bug. If it is square — and in hidden layers it almost always is, 512→512 — nobody complains: the code runs, the loss may even go down a little, and you spend two days wondering why the network learns badly.

The rule to take away. A transposed gradient on a square matrix produces no error at all: it produces a wrong number. This is why gradient checking is not a luxury for perfectionists but the first thing to write, before the training loop itself.

The five identities that actually matter

There are manuals with hundreds of matrix calculus identities — The Matrix Cookbook collects very many. In machine learning practice five of them are needed, and we derive all five component by component, never writing "it can be shown that". The method is always the same: write the function with explicit indices and sums, differentiate with respect to the k-th component as you would with a single-variable function, then recognise a matrix product in the result.

The first two: the dot product and the quadratic form

Identity 1. f(x) = aᵀx = Σᵢ aᵢ xᵢ. Differentiate with respect to x_k: in the sum there is a single term containing x_k, namely a_k x_k; all the others are constants and vanish. Hence ∂f/∂x_k = a_k. Stacking the n components:

∂(aᵀx)/∂x = a.

In words: the derivative of a linear function is its vector of coefficients, exactly as the derivative of ax is a in the scalar case. ∎

Identity 2. f(x) = xᵀAx = Σᵢ Σⱼ xᵢ Aᵢⱼ xⱼ. Here x_k appears in two families of terms: those with i = k (that is, Σⱼ x_k A_kⱼ xⱼ) and those with j = k (that is, Σᵢ xᵢ Aᵢ_k x_k), plus the term A_kk x_k² which belongs to both families. Differentiate with the product rule, treating the two roles separately:

∂f/∂x_k = Σⱼ A_kⱼ xⱼ + Σᵢ xᵢ Aᵢ_k = (Ax)_k + (Aᵀx)_k

The double counting of the diagonal term sorts itself out: ∂(A_kk x_k²)/∂x_k = 2A_kk x_k, which is precisely the contribution the two sums give together. Stacking:

∂(xᵀAx)/∂x = (A + Aᵀ)x, which becomes 2Ax only if A is symmetric. ∎

That "only if" is the second trap of the lesson. The formula 2Ax is so common that it looks universal, but it holds for symmetric matrices — covariances, AᵀA, Gram matrices. On an arbitrary matrix it is wrong, and wrong silently, because the shape of the result is correct.

import numpy as np

def riga(v, cifre=4):
    """Print an array with a fixed number of digits: two runs give the same text."""
    return " ".join(f"{val:+.{cifre}f}" for val in np.ravel(v))

def gradiente_numerico(f, x, h=1e-5):
    """Centred differences, one component at a time. Works on matrices too."""
    g = np.zeros_like(x, dtype=float)
    for i in range(x.size):
        passo = np.zeros_like(x, dtype=float)
        passo.flat[i] = h
        g.flat[i] = (f(x + passo) - f(x - passo)) / (2 * h)
    return g

a = np.array([1.0, 2.0, 3.0])
x = np.array([4.0, 5.0, 6.0])
print("f(x) = aᵀx")
print("  analytic a         :", riga(a))
print("  numeric            :", riga(gradiente_numerico(lambda v: a @ v, x)))

A = np.array([[2.0, 1.0],
              [0.0, 3.0]])                       # NOT symmetric, on purpose
x2 = np.array([1.0, 2.0])
print("f(x) = xᵀAx, with A non-symmetric")
print("  (A + Aᵀ)x          :", riga((A + A.T) @ x2))
print("  2Ax (wrong here)   :", riga(2 * (A @ x2)))
print("  numeric            :", riga(gradiente_numerico(lambda v: v @ A @ v, x2)))

Output obtained:

f(x) = aᵀx
  analytic a         : +1.0000 +2.0000 +3.0000
  numeric            : +1.0000 +2.0000 +3.0000
f(x) = xᵀAx, with A non-symmetric
  (A + Aᵀ)x          : +6.0000 +13.0000
  2Ax (wrong here)   : +8.0000 +12.0000
  numeric            : +6.0000 +13.0000

What the code does. gradiente_numerico perturbs one component at a time by ±h and computes the centred difference quotient; it uses passo.flat[i] so the same function works for vectors and matrices, while riga formats to a fixed number of digits because NumPy's printing changes alignment depending on the values. For aᵀx the numerical result matches a. For the quadratic form with non-symmetric A it gives (6, 13), which is (A + Aᵀ)x, not 2Ax = (8, 12): the shortcut is wrong on both components, with the right shape. Centred differences act here as an impartial judge, and they are exact on second-degree polynomials — a property we shall prove when we come to gradient checking.

Least squares and the outer product

Identity 3. f(x) = ‖Ax − b‖². Write the residual with indices: rᵢ = Σⱼ Aᵢⱼ xⱼ − bᵢ, and f = Σᵢ rᵢ². By the chain rule, ∂f/∂x_k = Σᵢ 2 rᵢ · ∂rᵢ/∂x_k, and ∂rᵢ/∂x_k = Aᵢ_k because bᵢ is constant. Hence ∂f/∂x_k = 2 Σᵢ Aᵢ_k rᵢ, which is the k-th component of Aᵀr:

∂‖Ax − b‖²/∂x = 2Aᵀ(Ax − b). ∎

This is the derivative of the least squares seen in the statistics lesson. Setting it to zero gives AᵀAx = Aᵀb, that is, the normal equation from the lesson on vectors and matrices: the minimum point is the solution of a linear system, and the whole of linear regression sits in these two lines.

Identity 4. Now we differentiate with respect to a matrix. Let z = Wx with W of shape (m, n), and let L be a scalar loss that depends on W only through z. Call δ = ∂L/∂z (shape (m,), the gradient arriving from above). By the chain rule, ∂L/∂W_pq = Σᵢ δᵢ · ∂zᵢ/∂W_pq. But zᵢ = Σⱼ Wᵢⱼ xⱼ, so ∂zᵢ/∂W_pq is zero unless i = p, in which case it equals x_q. The sum collapses onto a single term:

∂L/∂W_pq = δ_p x_q, that is, ∂L/∂W = δ xᵀ, the outer product between the gradient coming down and the input that had gone up. ∎

It is worth pausing, because it is the heart of backpropagation: the gradient of a weight is the product of just two numbers, what entered that weight and what came out of the error. No giant Jacobians: an outer product of cost m·n, as much as the forward product.

A = np.array([[1.0, 0.0],
              [1.0, 1.0],
              [0.0, 2.0]])
b = np.array([1.0, 1.0, 1.0])
x = np.array([1.0, 2.0])

r = A @ x - b
print("residual r = Ax − b  :", riga(r))
print("‖Ax − b‖²            :", f"{float(r @ r):.4f}")
print("analytic 2Aᵀ(Ax − b) :", riga(2 * A.T @ r))
print("numeric              :", riga(gradiente_numerico(lambda v: float((A @ v - b) @ (A @ v - b)), x)))

x_stella = np.linalg.solve(A.T @ A, A.T @ b)      # normal equation AᵀA x = Aᵀb
print("solution of AᵀA x = Aᵀb:", riga(x_stella, 6))
print("the gradient at the solution is zero:", np.allclose(2 * A.T @ (A @ x_stella - b), 0.0))

W = np.array([[1.0, 2.0, 3.0],
              [4.0, 5.0, 6.0]])                   # shape (2, 3)
x3 = np.array([1.0, 0.0, -1.0])                   # shape (3,)
delta = np.array([10.0, 20.0])                    # ∂L/∂(Wx), shape (2,)
print("shape of W:", W.shape, "| shape of δxᵀ:", np.outer(delta, x3).shape)
print("analytic δxᵀ :", riga(np.outer(delta, x3)))
print("numeric      :", riga(gradiente_numerico(lambda M: delta @ (M @ x3), W)))
print("xδᵀ would have shape", np.outer(x3, delta).shape, "→ with a square W the error would not show")

Output obtained:

residual r = Ax − b  : +0.0000 +2.0000 +3.0000
‖Ax − b‖²            : 13.0000
analytic 2Aᵀ(Ax − b) : +4.0000 +16.0000
numeric              : +4.0000 +16.0000
solution of AᵀA x = Aᵀb: +0.777778 +0.444444
the gradient at the solution is zero: True
shape of W: (2, 3) | shape of δxᵀ: (2, 3)
analytic δxᵀ : +10.0000 +0.0000 -10.0000 +20.0000 +0.0000 -20.0000
numeric      : +10.0000 +0.0000 -10.0000 +20.0000 +0.0000 -20.0000
xδᵀ would have shape (3, 2) → with a square W the error would not show

What the code does. The first part verifies identity 3: residual (0, 2, 3), loss 13, analytical gradient (4, 16) equal to the numerical one; then it solves the normal equation, and the point (7/9, 4/9) makes the gradient vanish — the very definition of a least-squares solution. The second verifies identity 4 by differentiating with respect to a matrix: gradiente_numerico perturbs one cell of W at a time and rebuilds a (2, 3) table identical to δxᵀ. The last line prints the shape of the transposed version, (3, 2): incompatible here, perfectly compatible if W were square.

Softmax and cross entropy: where ŷ − y comes from

The last identity is the one that in every course appears as a gift from heaven. The model produces K raw scores z (the logits); the softmax turns them into probabilities ŷᵢ = e^{zᵢ} / Σ_j e^{z_j}; the cross entropy against the true label y (a vector with one 1 and all 0s) is L = −Σᵢ yᵢ log ŷᵢ. Let us differentiate L with respect to z_k.

First step: the logarithm of the softmax simplifies. log ŷᵢ = zᵢ − log(Σ_j e^{z_j}), hence L = −Σᵢ yᵢ zᵢ + (Σᵢ yᵢ) · log(Σ_j e^{z_j}). Since Σᵢ yᵢ = 1, what remains is L = −Σᵢ yᵢ zᵢ + log(Σ_j e^{z_j}).

Second step: differentiate the two terms. The first is linear and gives −y_k. The second is the logarithm of a sum: ∂/∂z_k log(Σ_j e^{z_j}) = e^{z_k} / Σ_j e^{z_j} = ŷ_k. Adding up:

∂L/∂z = ŷ − y. ∎

In words: the gradient on the logits is the difference between what the model said and what it should have said. Two consequences: there is never any need to build the K × K Jacobian of the softmax, because its composition with cross entropy collapses into a subtraction; and the components sum to zero, because both ŷ and y sum to 1 — if one logit rises, another falls.

def softmax(z):
    z = z - z.max()                              # numerical stability: see the IEEE 754 lesson
    e = np.exp(z)
    return e / e.sum()

z = np.array([1.0, 2.0, 0.0])
y = np.array([0.0, 1.0, 0.0])                    # the true class is the second one
p = softmax(z)
perdita = lambda v: float(-np.log(softmax(v) @ y))

print("ŷ = softmax(z)  :", riga(p))
print("loss −log ŷ_y   :", f"{perdita(z):.4f}")
print("analytic ŷ − y  :", riga(p - y))
print("numeric         :", riga(gradiente_numerico(perdita, z)))
print("the sum of the gradient is zero:", np.allclose((p - y).sum(), 0.0))

Output obtained:

ŷ = softmax(z)  : +0.2447 +0.6652 +0.0900
loss −log ŷ_y   : 0.4076
analytic ŷ − y  : +0.2447 -0.3348 +0.0900
numeric         : +0.2447 -0.3348 +0.0900
the sum of the gradient is zero: True

What the code does. Subtracting the maximum inside softmax does not change the mathematical result but it stops np.exp overflowing on large logits: the attention to floating point of the lesson on floating-point numbers. The analytical gradient ŷ − y matches the numerical one to four digits and sums to zero; the correct class has a negative gradient (its logit must rise), the other two a positive one.

Trade-off. The fused "softmax plus cross entropy" version is more stable and faster than the two pieces separately, which is why libraries expose a single function (in PyTorch it is cross_entropy, which expects logits and not probabilities). The price is not being able to slip anything in between: for a calibration or a per-class weighted loss you must go back to the two pieces and accept the worse stability.

The Jacobian and the vector-Jacobian product

So far the function at the top has always been scalar: a loss. But a network layer does not produce a number, it produces a vector. We need the object that describes the derivative of a vector-valued function.

Definition. Given F: Rⁿ → Rᵐ, the Jacobian at x is the matrix J of shape (m, n) with

Jᵢⱼ = ∂Fᵢ/∂xⱼ.

Read it by rows and by columns. Row i is the gradient of the i-th output: it says which inputs it depends on. Column j says what effect the j-th input has on all the outputs together. For m = 1 the Jacobian has a single row and is the transposed gradient: the gradient is a special case of the Jacobian, not a different object.

The chain rule as a product of Jacobians

The chain rule, in vector form, is a matrix multiplication. If z = G(y) and y = F(x), then

J_{G∘F}(x) = J_G(y) · J_F(x), with shapes (p, m) · (m, n) = (p, n).

It is the formula that justifies the whole of backpropagation: a network of L layers is a composition of L functions, and the derivative of the loss is the product of L Jacobians. The shapes match automatically, and that is a free check: if they do not fit together, you have the order of the factors wrong. The product is not commutative: swapping two Jacobians does not give "the same result but transposed", it gives a shape error or a wrong number.

The VJP, and why nobody materialises a Jacobian

Here the mathematics parts company with the engineering: writing the chain rule as a product of Jacobians is correct, but computing them one by one and then multiplying them would be memory suicide — and it is never necessary. In backpropagation you always need the Jacobian left-multiplied by a vector, the one arriving from the next layer. The operation u ↦ uᵀJ is called the vector-Jacobian product (VJP) and it is the only thing a library actually implements. For a linear layer z = Wx the Jacobian with respect to x is W itself, and the VJP is uᵀW: a matrix-vector product of cost m·n, the same as the forward pass, not its square. The dual is the JVP (Jacobian-vector product), w ↦ Jw, and it is what forward mode needs.

def G(v):
    """G: R² → R³. Three outputs depending on the two inputs in different ways."""
    x, y = v
    return np.array([x**2 * y, x + y**3, x * y])

def jacobiano_analitico(v):
    x, y = v
    return np.array([[2*x*y, x**2],
                     [1.0,   3*y**2],
                     [y,     x]])

def jacobiano_numerico(f, v, h=1e-5):
    colonne = []
    for i in range(v.size):
        passo = np.zeros_like(v)
        passo[i] = h
        colonne.append((f(v + passo) - f(v - passo)) / (2 * h))
    return np.stack(colonne, axis=1)              # one column per input

v = np.array([2.0, 3.0])
J = jacobiano_analitico(v)
print("shape of J:", J.shape, "(rows = outputs, columns = inputs)")
print("J analytic:", riga(J))
print("J numeric :", riga(jacobiano_numerico(G, v)))

u = np.array([1.0, 1.0, 1.0])
print("VJP uᵀJ, one backward pass         :", riga(u @ J))
print("JVP Jw with w = (1, 0), one forward pass:", riga(J @ np.array([1.0, 0.0])))

for n in (256, 1024, 4096, 16384):
    print(f"layer {n} → {n}: Jacobian of {n*n} numbers = {n*n*4/2**20:.2f} MiB in float32")

Output obtained:

shape of J: (3, 2) (rows = outputs, columns = inputs)
J analytic: +12.0000 +4.0000 +1.0000 +27.0000 +3.0000 +2.0000
J numeric : +12.0000 +4.0000 +1.0000 +27.0000 +3.0000 +2.0000
VJP uᵀJ, one backward pass         : +16.0000 +33.0000
JVP Jw with w = (1, 0), one forward pass: +12.0000 +1.0000 +3.0000
layer 256 → 256: Jacobian of 65536 numbers = 0.25 MiB in float32
layer 1024 → 1024: Jacobian of 1048576 numbers = 4.00 MiB in float32
layer 4096 → 4096: Jacobian of 16777216 numbers = 64.00 MiB in float32
layer 16384 → 16384: Jacobian of 268435456 numbers = 1024.00 MiB in float32

What the code does, and what it tells us. jacobiano_numerico builds the matrix by columns: it perturbs one input at a time and watches how all the outputs move together — forward mode implemented naively, at a cost of one pass per input. The analytical Jacobian matches; the VJP contracts the three outputs and returns two numbers, the JVP contracts the two inputs and returns three. The last loop counts bytes, and that is the part that weighs: a 4096 → 4096 layer has a 64 MiB Jacobian in single precision, a 16384 → 16384 one reaches a gigabyte, for one example and one layer; with a minibatch of 32 and ten layers you multiply by 320. No library does this, and not out of laziness: the Jacobian of a linear layer is the weight matrix, which you already have in memory, and the VJP uses it without copying it. The rule to take away: never build a Jacobian that you only need multiplied.

Hessian, curvature and Newton's step

The gradient says in which direction the function rises, but says nothing about how fast that direction changes. That information lives in the second derivatives.

Definition. For f: Rⁿ → R, the Hessian is the n × n matrix of second derivatives,

Hᵢⱼ = ∂²f/∂xᵢ∂xⱼ.

It is the Jacobian of the gradient: H = J_{∇f}. With continuous second derivatives Schwarz's theorem guarantees ∂²f/∂xᵢ∂xⱼ = ∂²f/∂xⱼ∂xᵢ, so H is symmetric. No detail: the spectral theorem from the lesson on eigenvalues, SVD and PCA applies to the Hessian — real eigenvalues, perpendicular eigenvectors, H = QΛQᵀ.

The eigenvalue test: minimum, maximum or saddle

At a point where the gradient vanishes, the nature of the point can be read off the eigenvalues of the Hessian. The reason is the second-order Taylor expansion: moving by d from the critical point, f(x + d) ≈ f(x) + ½ dᵀHd, and the sign of dᵀHd decides whether you go up or down. Writing d in the eigenvector basis, dᵀHd = Σᵢ λᵢ cᵢ²: each eigenvalue weights the square of the component along its eigenvector.

• All eigenvalues positive: dᵀHd > 0 in every direction, so you always go up. It is a local minimum.
• All negative: you always go down. It is a maximum.
• Mixed signs: you go up along some directions and down along others. It is a saddle, and in deep networks saddles are the vast majority of critical points.
• Some eigenvalue zero: second order does not decide, higher derivatives are needed. In practice it is a flat direction, and gradient descent drags itself along it.

The ratio λ_max/λ_min is the condition number of the Hessian, and it measures how squashed the bowl is. It is the same "narrow valley" as in the lesson on gradient descent: the larger the ratio, the more descent bounces between the walls instead of running down the floor.

Newton's step, derived from Taylor

If you know the curvature you can do better than moving a little in the direction opposite to the gradient: you can jump to the minimum of the quadratic approximation. Expand f to second order around x:

f(x + d) ≈ f(x) + ∇f(x)ᵀd + ½ dᵀHd

This is a function of d whose gradient we know how to take with identities 1 and 2: ∂/∂d = ∇f(x) + Hd, because H is symmetric and therefore ∂(½dᵀHd)/∂d = Hd. Setting it to zero gives Hd = −∇f(x), that is,

x ← x − H⁻¹∇f(x), Newton's step. ∎

In words: instead of choosing a step size by hand, you deduce it from the curvature. On an exactly quadratic function a single step lands in the minimum, from whatever point you start. And H⁻¹ need not be inverted: you solve Hd = −∇f, which costs less and is more stable — the lesson of the normal equation, once again.

def f2(v):
    """f(x₀, x₁) = 2x₀² + 3x₀x₁ + 4x₁² − 5x₀ − 6x₁"""
    return 2*v[0]**2 + 3*v[0]*v[1] + 4*v[1]**2 - 5*v[0] - 6*v[1]

def gradiente_f2(v):
    return np.array([4*v[0] + 3*v[1] - 5, 3*v[0] + 8*v[1] - 6])

def hessiano_numerico(f, x, h=1e-4):
    n = x.size
    H = np.zeros((n, n))
    for i in range(n):
        for j in range(n):
            ei = np.zeros(n); ei[i] = h
            ej = np.zeros(n); ej[j] = h
            H[i, j] = (f(x+ei+ej) - f(x+ei-ej) - f(x-ei+ej) + f(x-ei-ej)) / (4*h*h)
    return H

H = np.array([[4.0, 3.0],
              [3.0, 8.0]])                        # ∂²f/∂xᵢ∂xⱼ, worked out by hand
x0 = np.array([1.0, 1.0])
print("H analytic:", riga(H))
print("H numeric :", riga(hessiano_numerico(f2, x0)))
print("symmetric, as Schwarz requires:", np.allclose(H, H.T))

autov = np.linalg.eigvalsh(H)
print("eigenvalues of H:", riga(autov), "→ all positive:", bool(np.all(autov > 0)), "→ minimum")
print("condition number λ_max/λ_min:", f"{autov[-1] / autov[0]:.4f}")
sella = lambda v: v[0]**2 - v[1]**2
print("eigenvalues of x₀² − x₁²:", riga(np.linalg.eigvalsh(hessiano_numerico(sella, x0))), "→ saddle")

x = np.array([0.0, 0.0])
print("start:", riga(x, 6), "| f =", f"{f2(x):.6f}")
x_newton = x - np.linalg.solve(H, gradiente_f2(x))
print("after ONE Newton step:", riga(x_newton, 6), "| f =", f"{f2(x_newton):.6f}")
print("the gradient at the new point is zero:", np.allclose(gradiente_f2(x_newton), 0.0))

eta = 0.1
contrazione = max(abs(1 - eta*autov[0]), abs(1 - eta*autov[-1]))
print(f"gradient descent with η = {eta}: the error shrinks by {contrazione:.4f} per step")
print("steps to cut the error by a factor of 1000:", int(np.ceil(np.log(1e-3) / np.log(contrazione))))
for n in (8192, 10**6):
    print(f"n = {n}: Hessian of {n*n} numbers ({n*n*4/2**30:.2f} GiB in float32), n³/3 = {n**3/3:.2e} operations")

Output obtained:

H analytic: +4.0000 +3.0000 +3.0000 +8.0000
H numeric : +4.0000 +3.0000 +3.0000 +8.0000
symmetric, as Schwarz requires: True
eigenvalues of H: +2.3944 +9.6056 → all positive: True → minimum
condition number λ_max/λ_min: 4.0116
eigenvalues of x₀² − x₁²: -2.0000 +2.0000 → saddle
start: +0.000000 +0.000000 | f = 0.000000
after ONE Newton step: +0.956522 +0.391304 | f = -3.565217
the gradient at the new point is zero: True
gradient descent with η = 0.1: the error shrinks by 0.7606 per step
steps to cut the error by a factor of 1000: 26
n = 8192: Hessian of 67108864 numbers (0.25 GiB in float32), n³/3 = 1.83e+11 operations
n = 1000000: Hessian of 1000000000000 numbers (3725.29 GiB in float32), n³/3 = 3.33e+17 operations

What the code does. hessiano_numerico uses the four-point formula — two coordinates perturbed at a time, all the sign combinations, divided by 4h² — and confirms the Hessian worked out by hand. eigvalsh gives the eigenvalues, both positive: we are in a bowl, and the condition number says how squashed it is; x₀² − x₁² provides the counter-check with one eigenvalue per sign, that is, a saddle. Finally Newton's step starts from the origin and lands on the exact minimum (22/23, 9/23) in one step, with a vanishing gradient, while gradient descent with η = 0.1 contracts the error by max|1 − ηλᵢ| per step and needs 26 of them.

Why deep learning does not use Newton's step

After a derivation like that, the question is unavoidable: if Newton converges in one step and gradient descent in dozens, why does everyone use the second? There are three reasons, and all three are quantitative.

Memory. The Hessian has n² cells. For the single matrix W of shape (128, 64) from the next lesson, n = 8192 and the Hessian takes 0.25 GiB in single precision; for a model with a million parameters — a small model — it is a thousand billion numbers, almost four terabytes. It is not a matter of better hardware: it is quadratic growth against linear growth.

Time. Solving Hd = −∇f costs of the order of n³/3 operations: with n around a million that is 3·10¹⁷ operations per step, against the few tens of millions of one descent step. One better step is not worth a hundred thousand steps.

Conditioning. Even if you could, the Hessian of a deep network is ill-conditioned and often indefinite: at a saddle it has negative eigenvalues, and there −H⁻¹∇f does not point downhill but upwards. Newton, taken literally, runs towards the nearest critical point, which in a deep network is almost always a saddle.

Trade-off. What survives of second order is the idea, not the matrix: quasi-Newton methods (L-BFGS) keep a low-rank approximation of it with a few dozen vectors, and adaptive optimisers such as Adam use a diagonal estimate of the curvature, that is, n numbers instead of n². They cost almost as much as gradient descent and capture the part of the curvature that matters. The price is that the approximation can be wrong and it introduces further hyperparameters to choose.

Forward mode and reverse mode

The analogy: the bill of materials against the quotation

A company assembles a finished product out of many raw materials. News arrives that copper has gone up by a cent a kilo, and the director wants to know how much the cost of the finished product rises. You take the bill of materials and propagate the increase forwards: copper goes into the cables, the cables into the module, the module into the product. One pass from raw material to product, and you have the answer — for copper. If aluminium goes up tomorrow, you redo the whole pass.

The next day the director changes the question: "the finished product has gone up by a euro; whose fault is it?" Now the bill of materials is walked backwards: the euro is split among the sub-assemblies, each sub-assembly splits its share among its components, and in a single pass you have the responsibility of all the raw materials together. It is the reversed quotation, and in calculus it is called reverse mode.

The analogy holds on a less comfortable point too: to split things backwards you must remember how the assembly was put together. The forward bill of materials can be walked while forgetting what lies behind; the backward one requires having kept the structure. That is exactly the memory cost of reverse mode, the famous "tape".

The exact arithmetic, and where the choice comes from

Let us formalise. Let f: Rⁿ → Rᵐ be computable with C elementary operations.

• Forward mode: one pass propagates forwards, alongside each value, its derivative with respect to one chosen input. One pass gives one column of the Jacobian, that is, Jw for w a unit vector. For the whole Jacobian you need n passes: cost ≈ n · C.
• Reverse mode: one forward pass records the values, one backward pass propagates the derivative of one chosen output with respect to everything. One pass gives one row of the Jacobian, that is, uᵀJ. For the whole Jacobian you need m passes: cost ≈ m · C, with a constant factor between 2 and 4 and memory proportional to the length of the computation.

Here the mathematics ends and the engineering choice begins. The function to differentiate is the loss: a single number. So m = 1 and n is the number of parameters, between millions and billions. Reverse mode costs about two or three evaluations of the loss and returns all the gradients; forward mode would cost a billion of them. This is not a matter of style: it is a ratio of a billion to one.

The mirror case holds just as well: with one parameter and a thousand outputs (n = 1, m = 1000) forward mode would be a thousand times cheaper and would need no tape. Reverse mode is not "better": it is better when the outputs are few and the inputs many.

Now let us measure it. Take f(x) = Σᵢ xᵢ², with n inputs and one output, implement it in both modes and count the elementary operations. Counting is better than timing: a stopwatch also measures the computer, a counter measures only the algorithm and is reproducible on any machine.

class Duale:
    """Dual number a + b·ε with ε² = 0: carries the value and ONE directional derivative."""

    def __init__(self, valore, derivata, conta):
        self.valore, self.derivata, self.conta = valore, derivata, conta

    def __add__(self, altro):
        self.conta[0] += 2                        # one for the value, one for the derivative
        return Duale(self.valore + altro.valore, self.derivata + altro.derivata, self.conta)

    def __mul__(self, altro):
        self.conta[0] += 2
        return Duale(self.valore * altro.valore,
                     self.valore * altro.derivata + self.derivata * altro.valore, self.conta)

def modo_diretto(valori, i_seme, conta):
    """f(x) = Σ xᵢ²: one pass gives ONE partial derivative, the seeded one."""
    s = Duale(0.0, 0.0, conta)
    for i, v in enumerate(valori):
        xi = Duale(v, 1.0 if i == i_seme else 0.0, conta)
        s = s + xi * xi
    return s.derivata

def modo_inverso(valori, conta):
    """The same f: one forward pass and one backward pass give ALL the derivatives."""
    s = 0.0
    for v in valori:                              # forward pass
        q = v * v;  conta[0] += 1
        s = s + q;  conta[0] += 1
    g = np.zeros(len(valori))
    aggiunto_s = 1.0                              # the seed: ∂f/∂f = 1
    for i, v in enumerate(valori):                # backward pass
        aggiunto_q = aggiunto_s;      conta[0] += 1          # ∂s/∂qᵢ = 1
        g[i] = aggiunto_q * (2 * v);  conta[0] += 1          # ∂qᵢ/∂xᵢ = 2xᵢ
    return g

rng = np.random.default_rng(7)
print(f"{'n':>4} | {'forward':>8} | {'reverse':>8} | {'ratio':>8} | same gradient")
for n in (2, 10, 100):
    valori = rng.normal(size=n)
    c_d = [0]; g_d = np.array([modo_diretto(valori, i, c_d) for i in range(n)])
    c_i = [0]; g_i = modo_inverso(valori, c_i)
    print(f"{n:>4} | {c_d[0]:>8} | {c_i[0]:>8} | {c_d[0]/c_i[0]:>8.1f} | {np.allclose(g_d, g_i)}")

Output obtained:

   n |  forward |  reverse |    ratio | same gradient
   2 |       16 |        8 |      2.0 | True
  10 |      400 |       40 |     10.0 | True
 100 |    40000 |      400 |    100.0 | True

What the code does, and what it tells us. The dual number is forward mode in object form: every value carries its own derivative with respect to the seeded input, and __add__ and __mul__ are the sum rule and the product rule written once and for all. For the full gradient the pass must be redone n times; reverse mode does one forward and one backward pass and produces all the derivatives in one go. The two gradients agree, being the same object computed in two different orders. The ratio between the two counts is exactly n, by construction and not by chance. In exchange reverse mode has had to keep something — here s and the value of each input, in a deep network every intermediate activation. The rule to take away: reverse mode buys time with memory.

An autograd from scratch

Now we have everything we need to write the real thing. A reverse-mode automatic differentiation system needs four ingredients, and none of the four is hard: an object that remembers where it comes from, a local rule for every operation, an order in which to walk the graph backwards, and a seed to start from.

The Nodo class, and why a closure

Every node keeps four things: the dato (the value computed forwards), the grad (the derivative of the loss with respect to that value, initially zero), the genitori, and _indietro, a function with no arguments that distributes its own gradient to the parents. The last is the trick that keeps the code short: the differentiation rule is written while the operation is executed, when the operands are still at hand, and kept as a closure to be called later. No table of derivatives, no dispatch on strings: the derivative is already packed inside the node.

def somma_su_broadcast(grad, forma):
    """Put the gradient back into the operand's shape: sum the axes added by broadcasting."""
    while grad.ndim > len(forma):
        grad = grad.sum(axis=0)
    for asse, dim in enumerate(forma):
        if dim == 1 and grad.shape[asse] != 1:
            grad = grad.sum(axis=asse, keepdims=True)
    return grad

class Nodo:
    """An array that remembers where it comes from: value, gradient, parents, backward rule."""

    def __init__(self, dato, genitori=(), operazione="foglia"):
        self.dato = np.asarray(dato, dtype=float)
        self.grad = np.zeros_like(self.dato)
        self.genitori = genitori
        self.operazione = operazione
        self._indietro = lambda: None             # leaves have nothing to distribute

    def __add__(self, altro):
        altro = altro if isinstance(altro, Nodo) else Nodo(altro)
        fuori = Nodo(self.dato + altro.dato, (self, altro), "+")
        def propaga():                            # a sum passes the gradient through untouched
            self.grad += somma_su_broadcast(fuori.grad, self.dato.shape)
            altro.grad += somma_su_broadcast(fuori.grad, altro.dato.shape)
        fuori._indietro = propaga
        return fuori

    def __mul__(self, altro):
        altro = altro if isinstance(altro, Nodo) else Nodo(altro)
        fuori = Nodo(self.dato * altro.dato, (self, altro), "*")
        def propaga():                            # product rule, element by element
            self.grad += somma_su_broadcast(altro.dato * fuori.grad, self.dato.shape)
            altro.grad += somma_su_broadcast(self.dato * fuori.grad, altro.dato.shape)
        fuori._indietro = propaga
        return fuori

    def __matmul__(self, altro):
        fuori = Nodo(self.dato @ altro.dato, (self, altro), "@")
        def propaga():
            self.grad += fuori.grad @ altro.dato.T        # δWᵀ: the VJP of the linear layer
            altro.grad += self.dato.T @ fuori.grad        # Xᵀδ: the outer product, identity 4
        fuori._indietro = propaga
        return fuori

How to read the code. The three operations are the three lines of mathematics already derived. The sum has derivative 1 with respect to each addend and passes the gradient through untouched; the only complication is broadcasting, because adding a (16,) bias to a (600, 16) matrix makes the gradient come back with 600 rows, and it must be compacted back by summing the added axis. The matrix product is identity 4 and its twin: gradient on the input δWᵀ, on the weights Xᵀδ, the outer product generalised to the minibatch. And gradients are accumulated (+=), not assigned: if a node is used twice the two contributions must add up.

The operators, the topological order and the seed

What is missing is the non-linearities, the final reduction and the engine that walks the graph. The order is not free: a node can distribute its own gradient only after receiving the contributions of all its children. The right order is the topological order reversed, which is what a depth-first visit produces when it appends a node after its parents.

    def relu(self):
        fuori = Nodo(np.maximum(self.dato, 0.0), (self,), "relu")
        def propaga():
            self.grad += (self.dato > 0) * fuori.grad     # at 0 the derivative is CHOSEN: 0 here
        fuori._indietro = propaga
        return fuori

    def tanh(self):
        t = np.tanh(self.dato)
        fuori = Nodo(t, (self,), "tanh")
        def propaga():
            self.grad += (1.0 - t**2) * fuori.grad        # tanh' = 1 − tanh²
        fuori._indietro = propaga
        return fuori

    def log(self):
        fuori = Nodo(np.log(self.dato), (self,), "log")
        def propaga():
            self.grad += fuori.grad / self.dato
        fuori._indietro = propaga
        return fuori

    def sum(self):
        fuori = Nodo(self.dato.sum(), (self,), "sum")
        def propaga():
            self.grad += np.ones_like(self.dato) * fuori.grad
        fuori._indietro = propaga
        return fuori

    def indietro(self):
        """Topological ordering, then the graph walked in reverse starting from the seed."""
        ordine, visti = [], set()
        def visita(nodo):
            if id(nodo) in visti:
                return
            visti.add(id(nodo))
            for genitore in nodo.genitori:
                visita(genitore)
            ordine.append(nodo)                   # after the parents: topological order
        visita(self)
        self.grad = np.ones_like(self.dato)       # the seed: ∂L/∂L = 1
        for nodo in reversed(ordine):
            nodo._indietro()

def entropia_incrociata(logit, y_onehot):
    """Softmax and cross entropy in a single node: the derivative is (ŷ − y)/m, identity 5."""
    z = logit.dato - logit.dato.max(axis=1, keepdims=True)
    e = np.exp(z)
    p = e / e.sum(axis=1, keepdims=True)
    m = p.shape[0]
    fuori = Nodo(-np.log((p * y_onehot).sum(axis=1)).mean(), (logit,), "entropia_incrociata")
    def propaga():
        logit.grad += (p - y_onehot) / m * fuori.grad
    fuori._indietro = propaga
    return fuori

How to read the code. The four non-linearities follow the same pattern: compute the forward value, build the node, pack the local derivative into a closure. relu deserves a note: at 0 it is not differentiable, and the code chooses 0 by writing self.dato > 0 instead of >= — a convention, not a theorem, and it is precisely there that the numerical check will fail. The indietro method does three things in ten lines: it visits the graph appending every node after its parents, puts the seed 1 on the root, calls the closures in reverse order. entropia_incrociata is a fused node: it uses identity 5 directly instead of differentiating softmax and logarithm separately.

Let us try it on an example checkable by hand: X = (1, −2, 3) as a row, W of shape (3, 2), loss L = sum(relu(XW)). The product gives Z = (4, 0): first channel active, second exactly on the ReLU's knee. The upstream gradient is therefore δ = (1, 0), and the two identities give ∂L/∂W = Xᵀδ and ∂L/∂X = δWᵀ: numbers you can count on your fingers.

X = Nodo([[1.0, -2.0, 3.0]])                      # (1, 3)
W = Nodo([[1.0, 2.0],
          [0.0, 1.0],
          [1.0, 0.0]])                            # (3, 2)
L = (X @ W).relu().sum()
L.indietro()
print("Z = X @ W :", riga(X.dato @ W.dato))
print("L         :", f"{float(L.dato):.4f}")
print("∂L/∂W     :", riga(W.grad), "  (expected Xᵀδ with δ = (1, 0))")
print("∂L/∂X     :", riga(X.grad), "  (expected δWᵀ)")

a = Nodo(3.0)
b = a * a + a                                     # a appears twice: contributions add up
b.indietro()
print("d(a² + a)/da at a = 3:", f"{float(a.grad):.4f}", "(expected 2a + 1 = 7)")

Output obtained:

Z = X @ W : +4.0000 +0.0000
L         : 4.0000
∂L/∂W     : +1.0000 +0.0000 -2.0000 +0.0000 +3.0000 +0.0000   (expected Xᵀδ with δ = (1, 0))
∂L/∂X     : +1.0000 +0.0000 +1.0000   (expected δWᵀ)
d(a² + a)/da at a = 3: 7.0000 (expected 2a + 1 = 7)

What the code does. It builds the graph of the figure and walks it backwards once. The gradient on W is the outer product between input and δ: first column (1, −2, 3), second all zeros because the ReLU closed the second channel. The one on X is δWᵀ. The second example checks accumulation: a enters a² + a twice, and the += adds the two contributions, giving 2a + 1 = 7. With = instead of += one of the two would be lost and the result would be 6 or 1 depending on the visiting order: a bug that depends on ordering is the worst kind.

Gradient checking done properly

An autograd written by hand is code like any other, and its output is a number nobody can eyeball. An independent oracle is needed, and the only one available without installing anything is finite differences. But they must be used with three precautions, and each of the three can be derived.

First precaution: centred differences, not forward ones. Expand f in a Taylor series at the two points x ± h:

f(x + h) = f(x) + h f'(x) + (h²/2) f''(x) + (h³/6) f'''(x) + …

f(x − h) = f(x) − h f'(x) + (h²/2) f''(x) − (h³/6) f'''(x) + …

Subtracting, the even terms cancel: f(x + h) − f(x − h) = 2h f'(x) + (h³/3) f'''(x) + …, whence

(f(x + h) − f(x − h)) / 2h = f'(x) + (h²/6) f'''(x) + …

The error is O(h²). The forward difference (f(x + h) − f(x))/h, by contrast, does not cancel the f'' term and is off by O(h). For the same h the centred version costs one extra evaluation and gains five orders of magnitude with h = 1e-5. It is the best bargain in the lesson. ∎

How large h should be

Here two errors that pull in opposite directions collide. Truncation is the one just computed, ≈ (h²/6)·|f'''|, and it decreases with h. Rounding is the one from the lesson on floating-point numbers: f(x + h) and f(x − h) are two extremely close numbers, and subtracting them is catastrophic cancellation; the relative uncertainty on each is the machine epsilon ε ≈ 2.2·10⁻¹⁶, so the absolute error on the difference is ≈ ε·|f| and, divided by 2h, it grows like ε·|f|/h.

How to read the table, row by row. The orders of magnitude come from the two formulas, with |f| ≈ 1 and |f'''| ≈ 6. The minimum of the sum is found by differentiating: d/dh (h²/3 + ε/h) = 0 gives h ≈ ∛ε ≈ 5·10⁻⁶. Hence the recipe h = 1e-5, not superstition but the cube root of the machine epsilon — and hence the apparent nonsense: shrinking h beyond that point makes the estimate worse.

Second precaution: the criterion is the relative error. Comparing |numerical − analytical| against an absolute threshold makes no sense, because the gradients of different layers differ by orders of magnitude. One uses |a − b| / max(|a|, |b|, η) with a small η protecting the division by zero, and accepts anything below 10⁻⁶.

Third precaution: do not check where the function is not differentiable. This is the mistake we shall see in the experiment, and it is worth anticipating: on a ReLU evaluated exactly at 0, finite differences return a sensible number that is different from any analytical convention.

def gradient_check(perdita_di, parametri, h=1e-5, eta=1e-12):
    """Compare the gradients already in .grad with centred differences. Criterion: RELATIVE error."""
    peggiore = 0.0
    for p in parametri:
        for i in range(p.dato.size):
            originale = p.dato.flat[i]
            p.dato.flat[i] = originale + h; piu = perdita_di()
            p.dato.flat[i] = originale - h; meno = perdita_di()
            p.dato.flat[i] = originale                       # ALWAYS restore
            numerico = (piu - meno) / (2 * h)
            analitico = float(p.grad.flat[i])
            scala = max(abs(numerico), abs(analitico), eta)
            peggiore = max(peggiore, abs(numerico - analitico) / scala)
    return peggiore

# (a) smooth loss: L = Σ (XW)², a second-degree polynomial in the weights
X = Nodo([[1.0, -2.0, 3.0]])
W = Nodo([[1.0, 2.0], [0.0, 1.0], [1.0, 0.0]])
Z = X @ W
L = (Z * Z).sum()
L.indietro()
perdita_liscia = lambda: float(((X.dato @ W.dato) ** 2).sum())
print("(a) smooth loss L = Σ (XW)²")
print("    analytic:", riga(W.grad))
print("    numeric :", riga(gradiente_numerico(lambda M: float(((X.dato @ M) ** 2).sum()), W.dato)))
print("    maximum relative error:", f"{gradient_check(perdita_liscia, [W]):.2e}")

# (b) the same network with the ReLU, evaluated EXACTLY on the knee
X = Nodo([[1.0, -2.0, 3.0]])
W = Nodo([[1.0, 2.0], [0.0, 1.0], [1.0, 0.0]])
L = (X @ W).relu().sum()
L.indietro()
perdita_relu = lambda: float(np.maximum(X.dato @ W.dato, 0.0).sum())
print("(b) loss with ReLU, at the point of non-differentiability")
print("    analytic:", riga(W.grad))
print("    numeric :", riga(gradiente_numerico(lambda M: float(np.maximum(X.dato @ M, 0.0).sum()), W.dato)))
print("    maximum relative error:", f"{gradient_check(perdita_relu, [W]):.4f}")

# (c) the U of the step: truncation on the left, cancellation on the right, on f(t) = t³ at t = 1
print("(c) choice of h on f(t) = t³ at t = 1, exact derivative 3")
for h in (1e-1, 1e-2, 1e-3, 1e-5, 1e-7, 1e-9, 1e-11, 1e-13):
    stima = ((1.0 + h) ** 3 - (1.0 - h) ** 3) / (2 * h)
    print(f"    h = {h:.0e}: estimate {stima:.12f} | relative error {abs(stima - 3.0) / 3.0:.3e}")

Output obtained:

(a) smooth loss L = Σ (XW)²
    analytic: +8.0000 +0.0000 -16.0000 +0.0000 +24.0000 +0.0000
    numeric : +8.0000 +0.0000 -16.0000 +0.0000 +24.0000 +0.0000
    maximum relative error: 2.68e-11
(b) loss with ReLU, at the point of non-differentiability
    analytic: +1.0000 +0.0000 -2.0000 +0.0000 +3.0000 +0.0000
    numeric : +1.0000 +0.5000 -2.0000 -1.0000 +3.0000 +1.5000
    maximum relative error: 1.0000
(c) choice of h on f(t) = t³ at t = 1, exact derivative 3
    h = 1e-01: estimate 3.010000000000 | relative error 3.333e-03
    h = 1e-02: estimate 3.000100000000 | relative error 3.333e-05
    h = 1e-03: estimate 3.000001000000 | relative error 3.333e-07
    h = 1e-05: estimate 3.000000000097 | relative error 3.246e-11
    h = 1e-07: estimate 3.000000000086 | relative error 2.876e-11
    h = 1e-09: estimate 3.000000081688 | relative error 2.723e-08
    h = 1e-11: estimate 3.000000248221 | relative error 8.274e-08
    h = 1e-13: estimate 2.999267501025 | relative error 2.442e-04

What the code does. Block (a) checks the autograd on a smooth loss: analytical 2·Xᵀ(XW) and numerical agree, with an error of the order of rounding. Block (b) repeats the check with the ReLU on an input that sends the second channel exactly to 0: the autograd answers 0 by convention, centred differences answer with the average between the slope on the left (zero) and the one on the right, and the relative error reaches 1. This is not a bug in the autograd, it is a badly posed question. Block (c) shows the U of the step.

What the sweep over h tells us. From 1e-1 to 1e-3 the relative error falls exactly like h²: 3.3·10⁻³, then 3.3·10⁻⁵, then 3.3·10⁻⁷, a hundred times smaller at every step — Taylor's prediction, confirmed digit by digit. Then the descent stops: between 1e-5 and 1e-7 the error is flat around 3·10⁻¹¹, because truncation is by now below rounding and no longer in charge. From there on it climbs back: 2.7·10⁻⁸ at 1e-9, 8.3·10⁻⁸ at 1e-11, and 2.4·10⁻⁴ at 1e-13, that is, almost a thousand times worse than the minimum. Note that at h = 1e-13 the estimate is 2.999267…: it has lost four significant digits. The rule to take away: the bottom of the U is three orders of magnitude wide, and any h between 10⁻⁶ and 10⁻⁴ will do; below 10⁻⁸ you are measuring noise.

The experiment: a two-armed spiral

We need a problem a straight line cannot solve, otherwise the gradient is never put to the test. Let us build two interleaved spirals: two arms starting at the origin half a turn out of phase. The data are simulated by us with numpy.random.default_rng, and we say so: what we are measuring is our own code, not nature.

The network is 2→16→16→2 with tanh in the two hidden layers and cross entropy at the output: 354 parameters, zero deep learning libraries. The last layer starts at zero, and the effect can be predicted exactly: with zero logits the softmax gives (0.5, 0.5) for every example and the initial loss is log 2 = 0.693147, whatever the first two layers do. It is the first test: if that number does not come out, there is an error in the forward pass and there is no point looking at the gradients.

def spirale_a_due_bracci(per_braccio=300, rumore=0.18, seme=42):
    """Two interleaved arms. Data SIMULATED by us, not a real dataset."""
    rng_dati = np.random.default_rng(seme)
    t = np.linspace(0.3, 3.2, per_braccio)
    blocchi, etichette = [], []
    for classe in (0, 1):
        angolo = 2.5 * t + classe * np.pi         # the second arm is half a turn out of phase
        px = t * np.cos(angolo) + rng_dati.normal(0.0, rumore, per_braccio)
        py = t * np.sin(angolo) + rng_dati.normal(0.0, rumore, per_braccio)
        blocchi.append(np.stack([px, py], axis=1))
        etichette.append(np.full(per_braccio, classe))
    return np.concatenate(blocchi), np.concatenate(etichette)

X_dati, y_dati = spirale_a_due_bracci()
Y = np.eye(2)[y_dati]                             # labels in one-hot form, (600, 2)
print("simulated data: X", X_dati.shape, "| y", y_dati.shape,
      "| class 0:", int((y_dati == 0).sum()), "| class 1:", int((y_dati == 1).sum()))

from sklearn.linear_model import LogisticRegression          # only to measure the linear wall
lineare = LogisticRegression().fit(X_dati, y_dati)
print("accuracy of a LINEAR model on the same data:", f"{lineare.score(X_dati, y_dati):.4f}")

rng = np.random.default_rng(0)
def strato(n_in, n_out, nullo=False):
    pesi = np.zeros((n_in, n_out)) if nullo else rng.normal(0.0, np.sqrt(2.0 / n_in), (n_in, n_out))
    return Nodo(pesi), Nodo(np.zeros(n_out))

W1, b1 = strato(2, 16)
W2, b2 = strato(16, 16)
W3, b3 = strato(16, 2, nullo=True)                # last layer at zero: logits start at zero
parametri = [W1, b1, W2, b2, W3, b3]
print("parameters of the 2→16→16→2 network:", sum(p.dato.size for p in parametri))

def avanti(X):
    H1 = (Nodo(X) @ W1 + b1).tanh()
    H2 = (H1 @ W2 + b2).tanh()
    return H2 @ W3 + b3

perdita = entropia_incrociata(avanti(X_dati), Y)
print("initial loss:", f"{float(perdita.dato):.6f}", "| log 2 =", f"{float(np.log(2.0)):.6f}")
print("initial accuracy:", f"{float((avanti(X_dati).dato.argmax(axis=1) == y_dati).mean()):.4f}")
perdita.indietro()
print("gradient shapes:", [p.grad.shape for p in parametri])
print("with the last layer zeroed, ∂L/∂W1 is exactly zero:", bool(np.abs(W1.grad).max() == 0.0))

Output obtained:

simulated data: X (600, 2) | y (600,) | class 0: 300 | class 1: 300
accuracy of a LINEAR model on the same data: 0.5567
parameters of the 2→16→16→2 network: 354
initial loss: 0.693147 | log 2 = 0.693147
initial accuracy: 0.5000
gradient shapes: [(2, 16), (16,), (16, 16), (16,), (16, 2), (2,)]
with the last layer zeroed, ∂L/∂W1 is exactly zero: True

What the code does, and what it tells us. The function generates 300 points per arm along a spiral of growing radius, with noise of standard deviation 0.18 — enough for the arms to brush past each other without overlapping; the logistic regression measures the wall, that is, how little can be achieved without non-linearity. Then it builds the network and checks the initial loss against log 2. Here there is a lesson inside the lesson: with the last layer at zero the gradient on W1 is exactly zero, not approximately, because the gradient coming down passes through W3 and if W3 = 0 nothing passes. Hence only the last layer is initialised at zero, and only for this test: keeping it as the definitive initialisation would freeze the two hidden layers at the first step. An all-zero initialisation is not "a neutral choice": it is a network that never learns.

The wrong transposition, on purpose

Now let us give the last layer real weights, check all 354 gradients with gradient checking, and deliberately introduce the mistake from the first section: swapping Xᵀδ with its transpose. On W2, which is 16×16, the program will not complain.

W3.dato = rng.normal(0.0, np.sqrt(2.0 / 16), (16, 2))        # the last layer is no longer zero
for p in parametri:
    p.grad = np.zeros_like(p.dato)                           # zero the gradients BEFORE every pass
perdita = entropia_incrociata(avanti(X_dati), Y)
perdita.indietro()

perdita_corrente = lambda: float(entropia_incrociata(avanti(X_dati), Y).dato)
errore = gradient_check(perdita_corrente, parametri)
print("gradient check on all 354 parameters: maximum relative error", f"{errore:.3e}")
print("below the 1e-6 threshold:", errore < 1e-6)

giusto = W2.grad.copy()
print("the transposed version has the same shape:", giusto.T.shape == giusto.shape)
print("but is it the same gradient?", bool(np.allclose(giusto, giusto.T)))
print("maximum gap between the gradient and its transpose:", f"{np.abs(giusto - giusto.T).max():.3e}")
print("on W1, which is (2, 16), the transposition immediately gives:", end=" ")
try:
    W1.dato - W1.grad.T
except ValueError as errore_di_forma:
    print(type(errore_di_forma).__name__, "— the shapes do not match")

for epoca in range(1, 3001):                                 # training with our hand-written autograd only
    perdita = entropia_incrociata(avanti(X_dati), Y)
    for p in parametri:
        p.grad = np.zeros_like(p.dato)
    perdita.indietro()
    for p in parametri:
        p.dato -= 0.5 * p.grad
    if epoca % 500 == 0:
        previste = avanti(X_dati).dato.argmax(axis=1)
        print(f"epoch {epoca:>4}: loss {float(perdita.dato):.4f} | accuracy {(previste == y_dati).mean():.4f}")

Output obtained:

gradient check on all 354 parameters: maximum relative error 1.078e-07
below the 1e-6 threshold: True
the transposed version has the same shape: True
but is it the same gradient? False
maximum gap between the gradient and its transpose: 1.074e-01
on W1, which is (2, 16), the transposition immediately gives: ValueError — the shapes do not match
epoch  500: loss 0.0582 | accuracy 0.9917
epoch 1000: loss 0.0117 | accuracy 0.9967
epoch 1500: loss 0.0069 | accuracy 0.9983
epoch 2000: loss 0.0047 | accuracy 0.9983
epoch 2500: loss 0.0035 | accuracy 1.0000
epoch 3000: loss 0.0027 | accuracy 1.0000

What the code does, and what it tells us. The block puts random weights back into the last layer, zeroes the gradients, redoes a full pass and then checks all 354 parameters against centred differences. The largest relative error is 1.078·10⁻⁷: autograd and numerical oracle agree to seven significant digits, and that is the proof the code is right — not "it seems to work", it is verified. The middle part documents the transposition: on W2 the shape matches and nobody complains, but the gradient differs from its transpose by up to 0.107, enormous compared with the gradients themselves, whose largest magnitude is of the order of 10⁻⁴; on W1, rectangular, NumPy raises a ValueError straight away. The last part trains for 3000 epochs on our autograd's gradients alone: from an initial loss of 0.693147 and an accuracy of 0.5000 it reaches 1.0000 on the 600 points, against the 0.5567 of the best linear separator. That difference between 0.56 and 1.00 is the value added by two non-linear layers, obtained with a hundred and twenty lines of NumPy. The accuracy is measured on the same data used for training and says nothing about generalisation: it only shows that the optimisation works.

Trade-off. Gradient checking on every parameter costs 2 × 354 forward passes: here a few seconds, on a real network impractical. In practice one checks a random sample (twenty or thirty parameters are enough), on a tiny minibatch, and once only, right after writing the code: it is not a regression test to leave running, it is a commissioning check.

Typical mistakes (symptom, cause, remedy)

1. The transposed gradient that converges anyway. Symptom: the network learns worse than expected, and only on square layers. Cause: xδᵀ instead of δxᵀ. Remedy: gradient checking on a small minibatch, not a shape check.
2. h too small in gradient checking. Symptom: the relative error grows as you shrink the step. Cause: catastrophic cancellation in the subtraction. Remedy: h = 1e-5, and if the test does not pass look for the bug, not for a smaller h.
3. In-place mutation that breaks the graph. Symptom: zero or wrong gradients, with no error at all. Cause: the backward closure uses the current value, not the one at the time. Remedy: update the parameters only after indietro(), and rebuild the graph on every pass.
4. Gradients not zeroed between two steps. Symptom: the loss diverges after a few epochs. Cause: the += inside the closures accumulates on the previous pass. Remedy: p.grad = np.zeros_like(p.dato) before every indietro() — it is the libraries' zero_grad(), and it is not a formality.
5. Derivative of the ReLU at zero. Symptom: gradient checking fails on a few weights and passes on all the others. Cause: at 0 finite differences return the average of the two slopes. Remedy: perturb the data by an epsilon, or check with tanh; do not "fix" the autograd.
6. Materialising a Jacobian. Symptom: running out of memory on a network that was perfectly happy in the forward pass. Cause: having built J instead of implementing uᵀJ. Remedy: always think in terms of VJPs.

The six rows have one thing in common: none of them produces an explicit error. The transposed gradient, the step that is too small, the in-place mutation and the gradients that were not zeroed give a wrong but plausible number, and that is why gradient checking is worth more than any review of the code by eye. The last two, the ReLU at zero and the materialised Jacobian, are instead false alarms or hidden costs: the first makes a correct check fail, the second exhausts memory on a network that ran fine going forward.

Best practices

• Declare the layout convention at the top of the file, in a comment, and never change it halfway through a project.
• Write the gradient check before the training loop, not after the first odd result.
• Check shapes with an assert at the entry of gradient functions: assert grad_W.shape == W.shape costs nothing and catches the transposition.
• Verify with two independent oracles: the closed-form formula derived by hand and centred differences. If both agree with your code, the code is right.
• Zero the gradients before every backward pass, and make it part of the loop, not an occasional gesture.
• Declare seed and versions (default_rng(42), NumPy 2.3.5): a gradient check that cannot be reproduced is not a verification.

The six rules read in two groups. The first three (declared convention, check written before training, assert on the shapes) prevent the most expensive silent error of the lesson; the last three (two oracles, zeroed gradients, seed and versions declared) make the verification repeatable by anyone. If you can adopt only one today, choose the check on the shapes: it costs one line.

Minimal glossary

• Gradient (gradiente): the vector of partial derivatives of a scalar; it has the shape of the object it is differentiated against.
• Jacobian (Jacobiano): the (m, n) matrix of the derivatives of a vector-valued function; one row per output.
• VJP (vector-Jacobian product): the operation u ↦ uᵀJ, the only one reverse mode really needs; its dual w ↦ Jw is the JVP.
• Hessian (Hessiano): the symmetric matrix of second derivatives; its eigenvalues say minimum, maximum or saddle.
• Forward mode (modo diretto): propagates derivatives alongside values; costs n passes, needs no memory.
• Reverse mode (modo inverso): one forward pass and one backward pass; costs m passes and the memory of the tape.
• Tape (nastro): the record of the operations performed, needed to walk the graph backwards.
• Seed (of the differentiation): the value 1 placed on the root before going down; it is ∂L/∂L.

The eight terms hold together with a single rule: every object has a shape, and every operation declares which shape it returns. Gradient and Jacobian differ by that shape, VJP and JVP by the side from which you multiply, forward and reverse mode by the direction in which you walk the graph. If you remember the shapes, the rest of the glossary follows from there.

Mini-project

Extend the lesson's autograd and verify it. All in NumPy, no deep learning libraries.

1. Add to Nodo the operators __sub__, __truediv__ and __pow__ with a constant exponent, writing the local derivative of each one in the closure; add sigmoide, and derive by hand that σ' = σ(1 − σ);
2. write controlla(nodo_perdita, parametri, quanti=20), which samples quanti parameters at random (with a declared seed) instead of checking them all and returns the worst relative error with the index of the guilty parameter;
3. deliberately break one local derivative — for instance 1 + t**2 instead of 1 − t**2 in tanh — and check that it finds the mistake and points to the right layer. Document the edge case: what happens if the sampled parameter has an exactly zero gradient?

Deliverable: a script that prints the worst relative error with the correct code, the one with the broken tanh, and the parameter on which the check lit up.

Challenge

• Easy. Pen and paper: derive ∂(xᵀAx)/∂x component by component for n = 2 with A = [[2, 1], [0, 3]], and check that the result is (A + Aᵀ)x and not 2Ax. Solution: f = 2x₀² + x₀x₁ + 3x₁², so ∂f/∂x₀ = 4x₀ + x₁ and ∂f/∂x₁ = x₀ + 6x₁, that is, [[4, 1], [1, 6]]x = (A + Aᵀ)x. With 2Ax you would get [[4, 2], [0, 6]]x: wrong on both components, with the right shape.
• Medium. Pen and paper: prove ∂‖Ax − b‖²/∂x = 2Aᵀ(Ax − b) by expanding (Ax − b)ᵀ(Ax − b), then show that setting the gradient to zero gives AᵀAx = Aᵀb. Solution: f = xᵀAᵀAx − 2bᵀAx + bᵀb; the first term is a quadratic form with AᵀA symmetric, so it differentiates to 2AᵀAx; the second is linear and gives −2Aᵀb; the third is constant. Adding up, 2Aᵀ(Ax − b). Setting it to zero gives the normal equation, with a unique solution when the columns of A are independent.
• Hard. Implement forward mode with dual numbers on a function f: Rⁿ → R of your choice and measure that the cost of the full gradient grows linearly in n, counting the elementary operations for n = 2, 5, 10, 20, 50, 100. Compare with reverse mode and check that the ratio is exactly n.
• Expert. Compute one row of the Hessian of the spiral loss with respect to the weights of the last layer by differentiating twice: obtain ∂L/∂W3 with the autograd, then apply centred differences to the gradient rather than to the loss, (∇L(W + h e_k) − ∇L(W − h e_k)) / 2h. Compare with the same row computed by the four-point formula on the loss, and explain why the first is more accurate for the same h.
• Enterprise. Take a cost function from your own work — a price that depends on a table of parameters, a risk score — and differentiate it by hand with respect to all its parameters, declaring the layout convention. Verify it with this lesson's gradient check and report the largest relative error. If it fails on some parameter, look for the point where the function is not differentiable: a max, a rounding, a threshold.

The five levels rise in difficulty by the kind of skill required: the first two are solved with pen and paper and consolidate the identities, the third and the fourth require you to implement and to compare two ways of differentiating, the last brings the method to a case from your own work. Do them in order, because each uses the result of the previous one as an oracle.

The coding exercises for this lesson are in a single file, esercizi.py (about 20 KB): the five identities verified with centred differences, the operation count in the two modes, the autograd to complete and gradient checking on a whole network. It requires Python 3.12 or later and NumPy 2.3. Run it with python3 esercizi.py: it works even with the parts to complete left empty and tells you which ones are missing.

Recap

Mental map. Declare the layout: ∂f/∂x has the shape of x. Five identities: ∂(aᵀx)/∂x = a; ∂(xᵀAx)/∂x = (A + Aᵀ)x; ∂‖Ax − b‖²/∂x = 2Aᵀ(Ax − b); ∂L/∂W = δxᵀ; softmax plus cross entropy gives ŷ − y. Jacobian (m, n), chain rule as a product of Jacobians, but in practice only VJPs. Symmetric Hessian, eigenvalues for minimum/maximum/saddle, Newton x ← x − H⁻¹∇f out of scale at O(n³). Forward mode costs n passes, reverse mode m: in ML m = 1. Autograd: value, grad, parents, closure; reversed topological order; seed 1. Gradient check: centred differences, h = ∛ε ≈ 1e-5, relative error.

Cheat sheet.

3 concepts: the gradient has the shape of the object it is differentiated against; the VJP is the only thing actually computed; reverse mode wins because the output is a single one. 3 mistakes: the transposed gradient on a square matrix; h too small in gradient checking; gradients not zeroed between two steps. 3 best practices: declare the layout convention; write the gradient check before the training loop; verify with two independent oracles. 1 exercise: the mini-project of this lesson.

Next lesson: optimisation with constraints, that is, where the multipliers the SVM uses without explaining them come from.

Questions to move forward

• How do you maximise a function when the variables are not free but must satisfy a constraint, such as ‖u‖ = 1 in PCA? (Lesson 11 of module 02: Lagrange multipliers, duality and the KKT conditions.)
• Why is the classification loss cross entropy in particular, and not any old measure of distance between distributions? (Lesson 12 of module 02: entropy, cross entropy and KL divergence.)
• How do you cut the memory of the tape when the network is too deep to fit? (Recomputing the intermediate values, gradient checkpointing: you trade time for memory, in the opposite direction to this lesson.)

The three questions are the hooks to the following lessons: the first two open lesson 11 and lesson 12 of module 02, the third is the other face of the time-memory trade-off met in the comparison of forward mode and reverse mode. They are not for revision: they are for deciding what to read next.

The same series also includes: 26 - Constrained Optimisation: Lagrange Multipliers, Duality and the KKT Conditions (the next lesson).