Jacobian-vector product

In the vast and complex systems that define modern science—from the climate models spanning the globe to the neural networks powering artificial intelligence—a central challenge is understanding change. How does tweaking one parameter out of millions affect the final outcome? The classical answer lies in the Jacobian matrix, a comprehensive map of all possible sensitivities. However, for systems with millions of variables, this map is often too colossal to compute or even store, creating a formidable barrier to analysis and optimization.

This article introduces a powerful and elegant solution: the ​​Jacobian-vector product (JVP)​​. It is a computational technique that allows us to determine the rate of change in a specific, chosen direction without ever constructing the full, impossibly large Jacobian matrix. By focusing on the action of the derivative rather than the derivative object itself, the JVP unlocks the ability to solve problems at a scale previously thought to be intractable.

We will embark on a journey to understand this pivotal concept. In the first chapter, ​​Principles and Mechanisms​​, we will delve into the heart of the JVP, exploring the clever techniques like automatic differentiation and finite differences that make it possible and examining its role as the linchpin of powerful Newton-Krylov solvers. Following this, the chapter on ​​Applications and Interdisciplinary Connections​​ will reveal the astonishing breadth of the JVP's impact, showcasing how this single idea serves as a unifying thread across fields as diverse as computational fluid dynamics, machine learning, quantum chemistry, and even evolutionary biology.

Imagine you are standing on a rolling hillside, and you want to know the slope. But you don't want to know the slope in every possible direction—north, east, southeast, and so on. You only care about the slope in the exact direction you are about to take your next step. Would you need to first create an exhaustive topographical map of the entire hill, detailing the slope in every conceivable direction, just to find the one value you're interested in? Of course not. You'd find a way to measure the change in elevation along your specific path.

This simple idea is the heart of the Jacobian-vector product. In the world of multivariable functions, which describe everything from weather patterns to the behavior of neural networks, the ​​Jacobian matrix​​, often written as $J$, is that complete topographical map. For a function $f$ that takes $n$ inputs and produces $m$ outputs, the Jacobian is an $m \times n$ matrix of all the possible partial derivatives. It tells you how every output changes in response to a small change in every input. The action of this matrix on a vector $v$, the product $Jv$, gives you the directional derivative: the rate of change of the function's output if you "move" the input in the direction of $v$. This product is the answer to our question: "what is the slope along my chosen path?"

The profound insight, which has revolutionized computational science, is this: if all we need is the slope along one path, we don't need to build the whole map. We can compute the Jacobian-vector product, or ​​JVP​​, directly.

How can we possibly compute the product $Jv$ without first knowing $J$? It feels a bit like magic, but it rests on the solid foundation of calculus. There are two primary techniques, each with its own flavor of elegance.

Automatic Differentiation: The Exact Path

The first method, known as ​​Automatic Differentiation (AD)​​, is a clever computational technique that treats the function not as a black box, but as a sequence of elementary operations (like addition, multiplication, sine, cosine). AD calculates the derivative by applying the chain rule to this sequence, step by step.

In its ​​forward mode​​, AD provides a beautiful way to compute a JVP. Imagine we want to compute the JVP for a function $f$ at a point $a$ in the direction $v$. We can use a special kind of number, a "dual number," of the form $x_{real} + \epsilon x_{dual}$, where $\epsilon$ is a curious symbol with the property that $\epsilon^2 = 0$. We set our input to be $a + \epsilon v$. Then, we simply evaluate the function $f(a + \epsilon v)$. As we propagate this dual number through each step of the function, the rules of calculus (encoded in how we define operations on dual numbers) do all the work for us. Because any term with $\epsilon^2$ vanishes, a Taylor expansion tells us the final result will look like $f(a) + \epsilon (J(a)v)$. The JVP we were looking for appears automatically as the "dual" part of the output!.

This is why forward-mode AD is sometimes called "tangent mode." It computes the function's value and its directional derivative (the tangent vector) simultaneously. Crucially, this calculation is not an approximation; it is exact up to the limits of computer floating-point precision, a stark contrast to other methods.

Finite Differences: The Intuitive Approximation

The second method is even more direct and intuitive. It takes us back to the very definition of a derivative. The derivative of a function is the limit of the change in the function divided by the change in the input. We can use this idea to approximate a JVP:

Here, we simply take a tiny step $h$ in the direction $\mathbf{v}$, evaluate our function $F$, see how much it changed from its value at $\mathbf{x}$, and divide by the step size $h$. This gives us an approximation of the directional derivative.

But this elegant simplicity hides a subtle and beautiful trade-off. How small should our step $h$ be? If we make $h$ too large, our linear approximation becomes poor, and we suffer from a large ​​truncation error​​ that grows with $h$. If we make $h$ too small, we fall victim to the limitations of our digital world. The subtraction in the numerator, $F(\mathbf{x}+h\mathbf{v}) - F(\mathbf{x})$, becomes a subtraction of two nearly identical numbers, a recipe for ​​catastrophic cancellation​​ in floating-point arithmetic. This round-off error gets amplified when we divide by the tiny $h$. This effect is even more pronounced if the function evaluations themselves are noisy, perhaps coming from a simulation or a physical experiment with inherent uncertainty $\delta$. In that case, the noise error scales like $\delta/h$.

So, there is a "sweet spot" for $h$—not too big, not too small—that optimally balances the truncation error and the noise or round-off error. This tells us something profound: the accuracy of our computed derivative is fundamentally limited by the precision of our tools, whether it's the finite precision of a computer or the noise in an experiment. Better accuracy can sometimes be squeezed out by using a more sophisticated formula, like a central difference, but the fundamental trade-off always remains.

Why go to all this trouble to avoid forming the Jacobian matrix? The answer is scale. In modern science and engineering, we often deal with functions that have millions or even billions of inputs. A neural network for image recognition or a finite element model of a car crash can easily have millions of parameters. The Jacobian matrix for such a system would have millions of rows and millions of columns, containing trillions ($10^{12}$) of numbers. Simply storing such a matrix is impossible on any current or foreseeable computer.

This is where the JVP becomes a superhero. The matrix is impossibly large, but the cost of computing a single JVP, using either AD or finite differences, is typically only a small constant multiple of the cost of evaluating the function itself. We can find the slope along any one path without ever needing the map.

This capability is the key that unlocks a powerful class of algorithms known as ​​Newton-Krylov methods​​. Many of the hardest problems in science boil down to solving a huge system of nonlinear equations, written as $F(x) = 0$. Newton's method is the classic way to do this: start with a guess, and iteratively improve it by solving a linear system based on the Jacobian, $J(x_k)s_k = -F(x_k)$, to find the next step $s_k$.

For large systems, we can't solve this linear system directly. Instead, we use iterative linear solvers, most famously the family of ​​Krylov subspace methods​​ (like GMRES). And here is the miracle: these solvers don't need to know the matrix $J(x_k)$ explicitly. All they require is a "black box" function that, given any vector $v$, returns the product $J(x_k)v$.

This is a perfect marriage. The Newton step needs a linear system solved. The Krylov solver can do it, provided it gets JVPs. And we have matrix-free ways to provide those JVPs! This synergy, known as the ​​Jacobian-Free Newton-Krylov (JFNK)​​ method, allows us to apply the power of Newton's method to problems of a scale that would have been unimaginable a few decades ago.

Practical JFNK methods include even more cleverness. They solve the linear system only approximately (an "inexact" Newton step), just enough to make progress on the outer nonlinear problem, saving immense computational effort. They also use ​​preconditioners​​—cheap, approximate versions of the Jacobian—to guide the Krylov solver and dramatically accelerate its convergence, even while the "true" JVP is still computed matrix-free.

The story doesn't end there. The JVP, $Jv$, asks how a change in the inputs affects the outputs. But we can also ask the "adjoint" question: how does a change in the outputs trace back to a change in the inputs? This corresponds to the ​​vector-Jacobian product (VJP)​​, written as $v^T J$.

This is the domain of ​​reverse-mode automatic differentiation​​, which is the engine behind the deep learning revolution, where it is famously known as ​​backpropagation​​. If you have a function with many inputs and only one output (like a loss function that measures the error of a neural network), reverse mode is astonishingly efficient. In a single "backward pass" through the function's operations, it can compute the VJP, which in this case gives you the entire gradient—the derivative of the scalar output with respect to all inputs. This is exactly what's needed to train a neural network or, in computational chemistry, to compute the forces acting on every atom in a molecule from a neural network potential energy surface.

And we can go even one level higher. What about second derivatives? The matrix of second derivatives is called the ​​Hessian matrix​​, $H$. It tells us about the curvature of a function. The Hessian is simply the Jacobian of the gradient function ($H = J(\nabla F)$). This means a ​​Hessian-vector product​​, $Hv$, is just a JVP of the gradient map! We can use our same bag of tricks—forward- or reverse-mode AD—to compute HVPs efficiently without ever forming the massive Hessian matrix. This is indispensable for advanced optimization algorithms and for quantifying the uncertainty in our models. The JVP is a unifying principle that extends from first to second derivatives and beyond.

What happens if we try to use these methods on a function that isn't smooth, one with sharp "kinks" or corners? Consider the absolute value function, $|x|$, or the ReLU function, $\max(0, x)$, which is the cornerstone of modern neural networks.

Away from the kink (e.g., for $|x|$ where $x \neq 0$), the function is perfectly smooth, and a finite-difference JVP will compute the exact derivative. A Newton step will send you straight to the solution.

But if your iterate lands exactly on the kink ($x=0$), something strange happens. The finite-difference formula $\frac{F(x_k+hv) - F(x_k)}{h}$ still produces a well-defined value. For $F(x)=|x|$ at $x_k=0$, it simply returns $|v|$. The problem is that the resulting operator, the map $v \mapsto |v|$, is no longer linear. This is a disaster for Krylov solvers like GMRES, which are built on the fundamental assumption of linearity. The method breaks down not because the derivative is undefined, but because the very structure of the problem has changed.

This failure is deeply instructive. It reveals the hidden machinery we rely on. It teaches us that the power of matrix-free Newton-Krylov methods comes not just from a clever computational trick, but from the beautiful interplay between the local linear structure of smooth functions and the algebraic properties of Krylov subspaces. When that local linearity disappears, so does the magic. This boundary case pushes us toward even more advanced ideas, like "semismooth" Newton methods, but it also solidifies our appreciation for the elegant and powerful principles that govern the smooth world where Jacobian-vector products reign supreme.

Now that we have explored the beautiful machinery of the Jacobian-vector product (JVP), we are ready for a grand tour. We have seen what the JVP is—the directional derivative of a function, the response of a system to a small push in a specific direction. But its true power and elegance are revealed only when we see it in action. Our expedition will take us across the vast landscape of modern science and engineering, from the swirling of galaxies to the intricate dance of electrons, from the design of an airplane wing to the very logic of the tree of life.

You will find that this single, beautiful idea is a kind of universal key, unlocking problems that at first seem impossibly large and unrelated. It is a testament to the profound unity of computational science. So, let us begin.

Much of our understanding of the universe is written in the language of differential equations. They describe how things change, evolve, and interact. But when we bring these elegant equations to a computer, discretizing space and time, they transform into something monstrous: a system of millions, or even billions, of coupled nonlinear algebraic equations.

Imagine trying to solve such a system head-on. The standard approach, Newton's method, requires calculating the system's Jacobian matrix—a matrix that could have trillions of entries, far too many to store, let alone invert. It is like trying to understand a city's traffic flow by creating a master chart of how every car's movement affects every other car, simultaneously. The task is hopeless.

This is where the JVP makes its grand entrance. The brilliant insight of so-called ​​Newton-Krylov methods​​ is that we don't need the entire blueprint of the Jacobian matrix. To solve the system, iterative methods like GMRES only need to know the action of the Jacobian on a given vector. They need to ask: if we make a small change $v$ to the state of our system, what is the first-order change in the system's governing equations? This is precisely what the JVP, $Jv$, tells us.

What is truly remarkable is how simply we can compute this action. We don't need to derive the complex analytical form of the Jacobian. We can approximate the JVP with a finite difference:

for some tiny step $h$. All we need is the ability to evaluate our system's equations, $F(\mathbf{x})$, which we must have anyway! This "matrix-free" approach is the secret behind many of the most powerful simulations of the physical world. It allows us to tackle problems of enormous scale in computational fluid dynamics, modeling the flow of air over a wing or the churning of a chemical reactor. It is the engine that drives the stable integration of "stiff" differential equations, which are notorious in combustion modeling and circuit simulation for having wildly different timescales. It is also the cornerstone of modeling complex reaction-diffusion systems, like those that create the beautiful Turing patterns seen in nature. In all these fields, the JVP allows us to probe the dynamics of a system without being crushed by its immense complexity.

Simulating the world as it is is one thing; changing it to be what we want is another. The JVP and its inseparable twin, the vector-Jacobian product (VJP), are the central tools for optimization, design, and learning.

Suppose we have a complex system—a bridge designed by the Finite Element Method (FEM), for instance—and we want to improve it. We can describe its shape with a set of parameters, $\boldsymbol{\theta}$. We also have a single objective we care about, like minimizing its weight while maintaining its strength. The question is: how does our objective change as we tweak all the design parameters?

Asking this question for each parameter one by one would be painfully slow. This is where the VJP, $\mathbf{w}^{\top}J$, comes to the rescue. By framing the problem in a special "adjoint" way, we can calculate the sensitivity of our single objective with respect to all parameters at once. This computation, which is the heart of ​​reverse-mode automatic differentiation (AD)​​, has a cost that is miraculously independent of the number of parameters. For the price of roughly one simulation, we get the complete gradient, telling us the most efficient way to improve our design.

This principle has fueled a revolution at the intersection of computational mechanics and machine learning. Imagine that the material properties of our bridge are not given by a simple textbook constant, but are instead described by a complex neural network. To "teach" this data-driven material model how to behave correctly, we need to adjust the millions of weights and biases of the network. The VJP provides the exact gradient of the simulation's error with respect to every single network parameter, forming a direct bridge between the physical simulation and the learning algorithm.

Perhaps the most impactful application of the VJP is one many of us use every day. The algorithm that trains almost all modern neural networks, ​​backpropagation​​, is nothing more than a clever, recursive algorithm for computing VJPs. The "loss function" is the single output, and the parameters are the millions of network weights. Backpropagation efficiently computes the gradient of this loss with respect to all weights by propagating sensitivities backward through the network's layers. This gradient is exactly what optimizers like Stochastic Gradient Descent use to make the network learn. Every time you speak to a voice assistant or use an image recognition app, you are witnessing the power of countless VJPs at work.

The reach of the JVP extends into the most fundamental domains of science, providing a computational microscope to explore complexity at its deepest levels.

Let's shrink down to the quantum realm. The "gold standard" for predicting the properties of molecules is a theory called Coupled Cluster (CC). Its equations are notoriously difficult, involving a dizzying number of tensor contractions. To calculate how a molecule interacts with light, one must solve an enormous eigenvalue problem for the CC Jacobian. As you might now guess, this matrix is never actually constructed. Instead, iterative solvers compute the excited states by repeatedly asking for the Jacobian's action on a vector—a perfect job for the JVP, computed via forward-mode AD.

Or consider the world of soft matter. The elegant, ordered patterns formed by self-assembling block copolymers are described by a sophisticated framework called Self-Consistent Field Theory (SCFT). The resulting system of integro-differential equations is formidable. Yet, here again, a matrix-free Newton-Krylov solver, powered by a JVP, tames this complexity. In a particularly beautiful twist, the JVP itself is calculated by solving a set of linearized diffusion equations, which describe how a perturbation to the field propagates along the polymer chains.

Finally, let us zoom out to the grandest scale of all: the tree of life. To understand how species and their traits evolved, biologists build statistical models of evolution on phylogenetic trees. The likelihood of the observed data (like DNA sequences from living species) is calculated using a recursive "pruning" algorithm that moves from the tips of the tree to its root. To fit the model's parameters, such as substitution rates or branch lengths, one needs the gradient of the log-likelihood. This entire calculation—a cascade of matrix-vector products and matrix exponentials across the tree—can be viewed as a single, giant computational graph. Applying reverse-mode AD to this graph allows for the efficient computation of the exact gradient of the likelihood with respect to every single parameter. This computation is, in its essence, one large VJP, and it is the engine powering modern statistical phylogenetics.

From the flow of air, to the training of AI, to the structure of molecules and the history of life, we have seen the same fundamental idea appear again and again. The Jacobian-vector product, in its various guises, is a unifying concept in computational science. It teaches us a profound lesson: to understand, predict, and optimize a complex system, we often do not need to know everything about it all at once. We just need to know how to ask the right question: "If I push here, what happens there?" The JVP is the embodiment of that question, and the elegant answer it provides is one of the great triumphs of our quest to understand the world.