Gradients, Jacobians, and Hessians

The derivative shape follows the dimensions of the input and output spaces: scalar to scalar gives a scalar, vector to scalar a gradient vector under a declared convention, vector to vector a Jacobian matrix, and a second derivative of a scalar function a square Hessian. The last row of the table is different in kind. Frameworks frequently never build the Jacobian at all. They evaluate Jacobian–vector and vector–Jacobian products instead, and the reason is a published cost law rather than an implementation habit.

Building the whole m×n Jacobian costs one sweep per input coordinate in forward mode, and one sweep per output coordinate in reverse mode. The same small constant sits in front of both. A 2018 survey of automatic differentiation in the Journal of Machine Learning Research sets the arithmetic out: “if we denote the operation count to evaluate the original function by ops(f), the time it takes to calculate the m×n Jacobian by the forward mode is nc ops(f), whereas the same computation can be done via reverse mode in mc ops(f), where c is a constant guaranteed to be c<6 and typically c∼[2,3] (Griewank and Walther, 2008)”. That constant is below 6, usually between 2 and 3. The survey draws the conclusion in one sentence: “That is to say, reverse mode AD performs better when m≪n.”

That inequality is why a scalar training loss over millions of parameters is differentiated backwards: m is 1 and n is the parameter count. JAX ships the asymmetry as an API rule rather than advice. “jacfwd uses forward-mode automatic differentiation, which is more efficient for 'tall' Jacobian matrices”, while “jacrev uses reverse-mode, which is more efficient for 'wide' Jacobian matrices”. The shape of the map does not only tell you which derivative object exists. It tells you what computing it will cost.

For f:ℝⁿ→ℝ, the gradient collects partial derivatives into a vector. For g:ℝⁿ→ℝᵐ, the derivative is a Jacobian matrix that maps small input changes to output changes. A Hessian stores second partial derivatives of a scalar function and describes how the gradient changes locally. These objects are related, but they are not interchangeable. Their shapes follow directly from the function's input and output dimensions.

Both matrices carry a nineteenth-century name, and both names can be traced to a printed page. The functional determinant has been Jacobi's since a memoir he published in Crelle's Journal in 1841. The determinant that became the Hessian arrived the following year, in an 1842 paper by Ludwig Otto Hesse on cubic and quadratic curves; the MacTutor history at St Andrews says he “is remembered particularly for introducing the Hessian determinant”.

The name itself came later, and from someone else: “The Hessian matrix was developed by Ludwig Otto Hesse (1811 – 1874), a German mathematician, though the term was first used by James Joseph Sylvester (1814 – 1897)”. The objects are older than the vocabulary now used to teach them. The vocabulary is what obscures that they answer different questions.

Both are matrices in common cases, which makes them easy to confuse. A matrix of the right dimensions will multiply without complaint even when it is the wrong object.

The Jacobian J is first-order and applies to vector-valued outputs. It is the local linear map g(x+Δ)≈g(x)+JΔ, with rows and columns following a stated convention, and it is what sensitivity analysis, coordinate changes and backpropagation are made of. The Hessian H is second-order and applies to scalar-valued outputs: an n×n square matrix for n inputs, symmetric when the mixed partials are regular. It enters the local quadratic model through the term ½ΔᵀHΔ, and the signs of its eigenvalues sort directions into upward, flat and downward. The gradient ∇f is first-order for a scalar output: one entry per input, producing directional derivatives through an inner product.

The distinction that matters in practice is the number of derivatives, not the number of indices. A Jacobian tells you where a perturbation goes. A Hessian tells you how the answer to that question is itself changing as you move. Reading one as the other is the error that shape checking exists to catch.

Let g(x,y)=(x²+y, xy). Its Jacobian is [[2x,1],[y,x]] — two rows because the output has two components, two columns because the input has two coordinates. The shape is fixed before any entry is computed.

Along a unit direction u, the directional derivative of f is ∇f(x)·u. The gradient therefore encodes local change along any direction through a single dot product. Under Euclidean geometry it points toward steepest local increase, and its negative points toward steepest decrease.

That last sentence carries a condition that teaching usually drops. Martens attaches it: “The negative gradient −∇h can be interpreted as the steepest descent direction for h in the sense that it yields the greatest instantaneous rate of reduction in h per unit of change in θ, where change in θ is measured using the standard Euclidean norm”. Steepest is defined relative to a way of measuring change in the parameters. Rescale the parameters and you have changed the norm. That changes the gradient coordinates, and it can change the apparent steepest direction even when the underlying function is the same function.

The directional-derivative view is also how gradients are tested in practice, because probing every coordinate is unaffordable. JAX's jax.test_util.check_grads compares automatic derivatives against finite differences, and documents the compromise plainly: “Gradients are only checked in a single randomly chosen direction”. One direction, one inner product, one number to compare. The gradient's defining property, used as a test procedure.

The output dimension determines whether local sensitivity is naturally a scalar, a vector, a matrix, or a higher-order object. Scalar to scalar gives an ordinary derivative: one local rate. Vector to scalar gives a gradient: one sensitivity coordinate per input dimension. Vector to vector gives a Jacobian: a local linear transformation between two spaces, m×n for f:ℝⁿ→ℝᵐ. Vector to scalar at second order gives a Hessian: local curvature across input directions, n×n and symmetric under regularity.

Writing the shapes first is not bookkeeping. It is the cheapest available test. A claimed derivative of the wrong shape cannot be repaired by checking the algebra inside it, because the object itself is answering a different question.

Use dimensions before symbols to test a derivative claim. First, write the input and output shapes: state f:ℝⁿ→ℝᵐ explicitly, before writing any partial derivative. Second, choose the derivative object the map actually calls for — gradient, Jacobian or Hessian. Third, test a perturbation: verify that the derivative maps Δx into the expected output shape. A product that will not conform is a claim that cannot be right. Fourth, check symmetry and conventions: confirm Hessian regularity and the row/column orientation your source uses for the Jacobian.

Fifth, validate numerically. Here the audit stops being a discipline you maintain by hand, because both major frameworks ship it as a tested utility. PyTorch's torch.autograd.gradcheck is documented as follows: “Check gradients computed via small finite differences against analytical gradients wrt tensors in inputs that are of floating point or complex type and with requires_grad=True.” Its published signature carries the tolerances it will judge you by — gradcheck(func, inputs, *, eps=1e-06, atol=1e-05, rtol=0.001, ...) — and one warning decides whether a failure is your bug or your dtype: “The default values are designed for input of double precision. This check will likely fail if input is of less precision, e.g., FloatTensor.”

JAX's equivalent is jax.test_util.check_grads, with signature check_grads(f, args, order, modes=('fwd', 'rev'), atol=None, rtol=None, eps=None). It compares automatic derivatives to finite differences, checks gradients in a single randomly chosen direction to stay affordable in high dimensions, and raises AssertionError on mismatch. Step five is not an exhortation to compare local changes at several scales. It is two function calls whose default tolerances are written down and can be read before you trust the result.