Hessian-Vector Products & the Pearlmutter Trick
Getting Hv without ever forming the (enormous) Hessian H.
On this page
The one-sentence idea: you almost never actually need the full Hessian matrix (section 1.16) itself — you usually only need to know what it does to a specific vector, and there's a way to get that answer directly, without ever building the (potentially enormous) matrix in between.
Beginner: for a neural network with a million parameters, the Hessian would be a million-by-million matrix — roughly a trillion numbers, far too large to store, let alone compute directly. But many algorithms only ever need the result of multiplying the Hessian by one particular vector, Hv, not the Hessian itself.
Intermediate: the Pearlmutter trick (also called the "R-op") computes exactly this Hessian-vector product using two ordinary backward passes (or a forward and backward pass), reusing the same automatic differentiation machinery from section 1.17 — at a cost comparable to just a couple of regular gradient computations, regardless of how many parameters the model has.
Advanced: the key algebraic observation is that Hv equals the gradient of the scalar quantity (∇f · v) — since that's a dot product of the gradient with a fixed vector, differentiating it again is just one more ordinary backpropagation pass, not a fundamentally different (and far more expensive) second-order operation.
The inner gradient is one ordinary backward pass; taking the gradient of the resulting scalar (the dot product) with respect to x is a second ordinary backward pass — two passes total, no matrix ever explicitly formed.
Write the gradient as a vector-valued function g(x) = ∇f(x), with components g_i(x) = ∂f/∂x_i. The Hessian is by definition the Jacobian of g: H_{ij} = ∂g_i/∂x_j. Now consider the scalar function φ(x) = g(x)·v = Σᵢ g_i(x)vᵢ for a fixed constant vector v, and differentiate it with respect to xⱼ, using that v doesn't depend on x:
(using that H is symmetric for a twice-differentiable f, so H_{ji} = H_{ij}). So the j-th component of ∇φ is exactly the j-th component of Hv — meaning ∇φ = Hv for every component simultaneously, which is the identity used above. Nothing here required ever writing down H itself; the derivation only ever manipulated the scalar function φ = ∇f · v, which is exactly why two ordinary backward passes suffice.
Where this is used: any second-order optimizer, curvature diagnostic, or influence-function computation that needs "Hessian times a vector" rather than the full Hessian relies on exactly this identity — it's the mathematical fact, not just an implementation trick, that makes second-order information tractable for million-parameter models.
For f(x, y) = x²y + y³, the true Hessian is H = [[2y, 2x], [2x, 6y]]. Pick v = [1, 0] and evaluate at (x,y) = (2,1): directly, Hv = [2(1), 2(2)] = [2, 4]. Via the trick: the gradient is ∇f = [2xy, x²+3y²] = [4, 7] at that point; the dot product with v is 4·1 + 7·0 = 4, a function of (x, y); differentiating that scalar expression, 2xy, with respect to x and y gives [2y, 2x] = [2, 4] — the same answer, reached without ever writing down the full 2×2 Hessian matrix.
create_graph=True on the first call is the crucial detail — it keeps the computation graph alive so the second grad call can differentiate through the first gradient itself.
import torch
x = torch.tensor([2.0, 1.0], requires_grad=True)
def f(x):
return x[0]**2 * x[1] + x[1]**3
v = torch.tensor([1.0, 0.0])
grad = torch.autograd.grad(f(x), x, create_graph=True)[0] # first backward pass
hvp = torch.autograd.grad(grad @ v, x)[0] # second backward pass -> H @ v
print(hvp) # tensor([2., 4.]) -- matches the hand-derivation above- K-FAC and other second-order optimizers use Hessian-vector products (or close approximations) to take smarter, curvature-aware update steps than plain gradient descent, without the prohibitive cost of forming a full Hessian.
- Influence functions (estimating how much a single training example affected a trained model's predictions) rely on efficiently solving linear systems involving the Hessian, which is only tractable via Hessian-vector products plus an iterative solver like conjugate gradient (section 1.34).
- Sharpness-aware training methods use the largest eigenvalue of the Hessian (estimated via power iteration, section 1.24, using nothing but repeated Hessian-vector products) as a proxy for how "flat" or "sharp" a solution is.
- Forgetting
create_graph=Trueon the first gradient call — without it, the graph needed for the second differentiation pass is discarded and the Hessian-vector product silently cannot be computed. - Attempting to form the full Hessian explicitly "just to be safe" on a large model — this defeats the entire point and can exhaust memory instantly; always reach for Hessian-vector products when only directional curvature information is actually needed.
Going deeper
The same trick generalizes to Jacobian-vector products for any vector-valued function, not just gradients of scalar losses — this is exactly what "forward-mode automatic differentiation" (mentioned in section 1.17's expert note) computes directly and efficiently.
At the master level: combining Hessian-vector products with the power iteration algorithm of section 1.24 (using Hv repeatedly instead of a fixed matrix multiply) is precisely how practitioners estimate a neural network loss landscape's sharpest curvature directions without ever materializing a Hessian that would otherwise be far too large to store.