Gradients, Jacobians, and Hessians: Calculus in Many Dimensions

A loss depends on millions of parameters at once, but the derivative machinery only knows how to handle one. So freeze all but one and differentiate normally. That is a partial derivative, written $\partial f / \partial x_i$, and it answers a deliberately narrow question: how does $f$ respond to this coordinate alone, with everything else held still?

Collect all of them into a vector and you have the gradient:

$\nabla f(x) = \left[\frac{\partial f}{\partial x_1},\ \frac{\partial f}{\partial x_2},\ \dots,\ \frac{\partial f}{\partial x_n}\right]^{\mathsf{T}}$

Definition: The gradient of a scalar function is the vector of its partial derivatives. It has exactly as many entries as the function has inputs, which is why a gradient in a language model is the same size as the model.

Partial derivatives only describe moves along the axes. To move in an arbitrary direction $u$, with $\|u\| = 1$, take the directional derivative:

$D_u f(x) = \nabla f(x) \cdot u$

Now ask which unit $u$ makes that largest. It is a dot product with a fixed vector, so the Cauchy-Schwarz inequality answers immediately: it is maximised when $u$ points the same way as $\nabla f(x)$, and the maximum value is $\|\nabla f(x)\|$.

Two facts come out of that one line. The gradient's direction is the steepest ascent, which is why descent moves the other way. Its magnitude is the steepest possible rate of change at that point, which is the number people mean by "gradient norm" when they watch it explode during training.

That derivation quietly assumed something. "Unit vector" means $\|u\| = 1$ in the ordinary Euclidean norm, so steepest ascent means steepest per unit of Euclidean distance in parameter space.

Change the units of one parameter and the geometry changes with them. Suppose $f$ depends on a weight measured in metres and you switch to kilometres. The function is the same function; the same physical move is now a different-sized step; and the direction the gradient calls steepest has rotated.

Gotcha: The gradient is not an intrinsic "downhill" arrow. It is downhill relative to a particular notion of distance, and that notion is an arbitrary choice inherited from how you happened to parameterise the model. Adam, RMSProp and every other per-parameter learning rate scheme is in effect a rescaling of that geometry, which is why they can outrun plain gradient descent without ever computing a better gradient.

A neural network layer is not a scalar function. It maps $n$ numbers to $m$ numbers, so there is a derivative for every input-output pair. Stack them into a matrix and you have the Jacobian $J$, with $J_{ij} = \partial f_i / \partial x_j$: one row per output, one column per input.

The gradient is just the Jacobian of a function with one output, transposed. Everything in the table is the same idea counted differently, and getting the count right is most of the work in practice.

Take a linear layer $y = Wx$ with $W$ of shape $(4, 3)$. The Jacobian $\partial y / \partial x$ is $W$ itself, shape $(4, 3)$. But the backward pass multiplies by its transpose:

import numpy as np

W = np.random.randn(4, 3)      # layer R^3 -> R^4
x = np.random.randn(3)
y = W @ x                      # forward

dL_dy = np.random.randn(4)     # gradient arriving from above
dL_dx = W.T @ dL_dy            # (3,)   : Jacobian transpose times a vector
dL_dW = np.outer(dL_dy, x)     # (4, 3) : same shape as W

Gotcha: Textbooks disagree on whether a derivative matrix is laid out as the Jacobian or its transpose, the numerator and denominator conventions. Frameworks sidestep the argument with a rule you can always fall back on: a gradient has the same shape as the thing it is a gradient of. dL_dW is shaped like W, always. If your shapes do not line up, you have the convention backwards.

Differentiate every entry of the gradient with respect to every input and you get the Hessian, the $n \times n$ matrix of second partial derivatives, $H_{ij} = \partial^2 f / \partial x_i \partial x_j$. For any function whose second partials are continuous, mixed partials commute, so $H$ is symmetric.

It is the coefficient of the quadratic term in the multivariable Taylor expansion, the direct analogue of $f''$ from the previous lesson:

$f(x + h) \approx f(x) + \nabla f(x)^{\mathsf{T}} h + \tfrac{1}{2} h^{\mathsf{T}} H h$

The gradient tells you which way the surface tilts. The Hessian tells you how fast that tilt is changing, which is the same as asking how far you can trust the gradient before it goes stale.

The Hessian is symmetric, so it has real eigenvalues and orthogonal eigenvectors. Each eigenvector is a direction, and its eigenvalue is the curvature along that direction: positive means the surface bends upward, negative means downward.

This is why "the optimiser got stuck in a local minimum" is usually the wrong diagnosis for a large model. Flat regions around saddle points, where the gradient is small but not zero, slow training far more often than genuine basins do.

The ratio of largest to smallest Hessian eigenvalue is the condition number $\kappa$. Geometrically it is how elongated the bowl is: $\kappa = 1$ is a perfect circle, large $\kappa$ is a narrow ravine where the steepest direction points across the valley rather than along it.

Kantorovich's classical analysis of steepest descent on a quadratic gives the price directly. Each step shrinks the gap to the optimum by a factor of $\left(\frac{\kappa - 1}{\kappa + 1}\right)^2$, so the iterations needed grow linearly with $\kappa$:

A badly conditioned problem is not harder to differentiate. It is harder to follow, and the cost is paid entirely in wall-clock time.

Curvature is genuinely useful, and the classic use of it, Newton's method, steps to $x - H^{-1}\nabla f$, which is condition-number-blind and converges in far fewer iterations. It is also unusable at scale, for a reason that is pure arithmetic.

A model with $10^9$ parameters has a Hessian with $10^{18}$ entries. In 32-bit floats that is 4 exabytes, and inverting it is worse than storing it.

In practice: What survives is the Hessian-vector product $Hv$. Pearlmutter showed in 1994 that $Hv$ can be computed exactly at roughly the cost of one gradient, without ever forming $H$, by differentiating the gradient in a single direction. Everything second-order that runs on real models, Newton-Krylov methods, curvature-aware optimisers, sharpness measurements, is built on that product rather than on the matrix.

Everything so far assumed you may move anywhere. Add a constraint $g(x) = 0$ and $\nabla f = 0$ stops being the right stopping condition, because the best point on the constraint surface is usually not flat.

The correct condition is still pure calculus. At a constrained optimum you cannot improve $f$ without leaving the surface, which happens exactly when the two gradients are parallel:

$\nabla f(x) = \lambda \nabla g(x)$

The multiplier $\lambda$ measures how much the optimum would improve if the constraint were relaxed slightly, which makes it a price rather than a bookkeeping device. Generalising this equation to inequality constraints produces the KKT conditions and Lagrangian duality, the subject of the separate Lagrangian Duality lesson in this catalogue. The calculus you need for it is the calculus in this lesson.