Multivariable calculus for training

Gradients, Jacobians and Hessians, the chain rule across layers, and what curvature tells you about why deep networks are hard to train.

Gradients and Jacobians

For a scalar loss L of many parameters, the gradient is a vector of the same shape as the parameters: each entry is the sensitivity of the loss to one parameter. For a vector-valued function f: R^n -> R^m, the derivative is the Jacobian, an m x n matrix whose row i is the gradient of output i.

import numpy as np

def f(x):
    return np.array([x[0] ** 2, x[0] * x[1], np.sin(x[1])])

def jacobian(fn, x, h=1e-6):
    x = np.asarray(x, dtype=float)
    out = np.zeros((len(fn(x)), len(x)))
    for i in range(len(x)):
        dx = np.zeros_like(x)
        dx[i] = h
        out[:, i] = (fn(x + dx) - fn(x - dx)) / (2 * h)
    return out

x = np.array([1.0, 2.0])
J = jacobian(f, x)
# analytic: [[2x0, 0], [x1, x0], [0, cos(x1)]]
print(np.round(J, 6))

• For a loss you almost never need the full Jacobian: reverse-mode autodiff computes a vector-Jacobian product, which is exactly one backward pass per output.
• Forward-mode autodiff computes Jacobian-vector products cheaply, which is why it wins when there are few inputs and many outputs.
• The gradient of a scalar loss has the same shape as the parameter tensor — a fast way to spot a wiring bug is a shape mismatch between loss gradient and parameter.

The chain rule through a network

A layer is a function; a network is a composition. The chain rule says the derivative of the composition is the product of the local derivatives, applied in order. Backpropagation is simply evaluating that product from the loss backwards, reusing shared intermediate results.

# two layers, hand-written backward pass
rng = np.random.default_rng(0)
X = rng.normal(size=(32, 4))
Y = rng.normal(size=(32, 1))

W1 = rng.normal(scale=0.1, size=(4, 8))
W2 = rng.normal(scale=0.1, size=(8, 1))

def relu(z):
    return np.maximum(z, 0.0)

for _ in range(200):
    H_pre = X @ W1                 # (32, 8)
    H = relu(H_pre)
    P = H @ W2                     # (32, 1)

    dP = 2 * (P - Y) / len(X)      # dL/dP for mean squared error
    dW2 = H.T @ dP                 # (8, 1)
    dH = dP @ W2.T                 # (32, 8)  gradient flows back through W2
    dH_pre = dH * (H_pre > 0)      # through the relu
    dW1 = X.T @ dH_pre             # (4, 8)

    for w, g in ((W1, dW1), (W2, dW2)):
        w -= 0.1 * g

Curvature, Hessians and fragility

# the Hessian is the Jacobian of the gradient: it describes curvature
def loss(w):
    return w[0] ** 2 + 20 * w[1] ** 2        # a long narrow valley

H = np.diag([2.0, 40.0])
np.linalg.eigvalsh(H)                       # [2, 40] -> condition number 20

# the largest safe learning rate for gradient descent is 2 / L
L = np.linalg.eigvalsh(H).max()
2 / L                                       # 0.05: above this it oscillates and diverges

• The condition number (largest eigenvalue / smallest) is how stretched the loss surface is. Gradient descent zigzags across a narrow valley instead of walking along it.
• Each Hessian eigenvalue is a curvature; the corresponding eigenvector is a direction. Ill-conditioning is the ratio between the worst and the best of them.
• Saddle points dominate in high dimensions, not local minima: at a random critical point the chance that every direction curves upward falls quickly with dimension.
• Computing a full Hessian is O(n^2) memory — for real networks you approximate it (diagonal, or a few eigenvalues via power iteration or Lanczos) or avoid it entirely with adaptive methods.

FAQ

Related

Gradients and calculus intuition Optimisation and gradient descent variants

Last refreshed 2026-09-18.