In From derivatives to backpropagation, we followed individual numbers through a computation graph. Each operation supplied local derivatives. During the backward pass, we multiplied those derivatives by incoming gradients and added contributions wherever paths met.

A neural-network layer usually processes many numbers together. One input can affect several outputs, and each output can depend on several inputs. How do we organize all those derivatives, and what do they tell us about the layer?

The Jacobian is a matrix of first-order partial derivatives. We’ll build one from a two-input, two-output function, use it to predict small changes, and see why its rank can vary from point to point. Then we’ll connect that geometric picture to the matrix chain rule and explain why backpropagation usually needs products with Jacobians rather than the full matrices.

From partial derivatives to a matrix

A partial derivative measures how one output responds to one input while the other independent inputs stay fixed. With several outputs, we repeat that question for each output–input pair.

Consider a function with two inputs and two outputs:

$$f(\mathbf x)= \begin{bmatrix} f_1(x_1,x_2)\\f_2(x_1,x_2) \end{bmatrix} = \begin{bmatrix} x_1^2+x_2\\x_1+x_2^2 \end{bmatrix}.$$

For the first output, changing $x_1$ acts through the square, while changing $x_2$ acts through addition. For the second output, those roles reverse:

Output	With respect to $x_1$	With respect to $x_2$
$f_1=x_1^2+x_2$	$2x_1$	$1$
$f_2=x_1+x_2^2$	$1$	$2x_2$

Arrange this table into a matrix, keeping outputs as rows and inputs as columns:

$$J_f(\mathbf x)= \begin{bmatrix} \frac{\partial f_1}{\partial x_1}&\frac{\partial f_1}{\partial x_2}\\[6pt] \frac{\partial f_2}{\partial x_1}&\frac{\partial f_2}{\partial x_2} \end{bmatrix} = \begin{bmatrix}2x_1&1\\1&2x_2\end{bmatrix}.$$

This is the Jacobian matrix. For a function $f:\mathbb R^n\to\mathbb R^m$, it has $m$ rows and $n$ columns. The entry in row $j$, column $i$ is $\partial f_j/\partial x_i$.

A derivative rule and its value at a point

The expression above is a rule for constructing the Jacobian at any input. At $\mathbf x=(1,2)$, it gives:

$$f(1,2)=\begin{bmatrix}3\\5\end{bmatrix}, \qquad J_f(1,2)=\begin{bmatrix}2&1\\1&4\end{bmatrix}.$$

The function gives the current outputs. The Jacobian tells us how those outputs respond to changes around the current inputs. For example, its first column says that a small increase in $x_1$ changes $f_1$ at rate $2$ and $f_2$ at rate $1$, holding $x_2$ fixed.

Its second column answers the corresponding question for $x_2$: the rates are $1$ and $4$. A column therefore collects one input’s effect on all outputs. A row collects all inputs’ effects on one output.

Where the gradient fits

These familiar objects organize the same kind of information:

Function	Derivative information
One input, one output	A scalar derivative
Several inputs, one scalar output	A gradient vector
Several inputs, several outputs	A Jacobian matrix

The Jacobian also covers the first two cases. For a scalar loss $L$, its Jacobian has one row. If we write gradients as column vectors, then:

$$J_L=\begin{bmatrix}\frac{\partial L}{\partial x_1}&\cdots&\frac{\partial L}{\partial x_n}\end{bmatrix}, \qquad \nabla_{\mathbf x}L=J_L^{\mathsf T}.$$

The superscript $\mathsf T$ means transpose: swap rows and columns. We’ll consistently use column vectors for gradients below.

The Jacobian predicts local changes

For a scalar function, we used $\Delta y\approx f'(x)\Delta x$. The vector version is:

$$\Delta\mathbf f\approx J_f(\mathbf x)\Delta\mathbf x.$$

Matrix multiplication combines the effects of the input changes. At $(1,2)$, suppose we change the inputs by $\Delta\mathbf x=(0.01,-0.02)$:

$$\Delta\mathbf f\approx \begin{bmatrix}2&1\\1&4\end{bmatrix} \begin{bmatrix}0.01\\-0.02\end{bmatrix} = \begin{bmatrix}0\\-0.07\end{bmatrix}.$$

For the first output, the two first-order effects cancel: $2(0.01)+1(-0.02)=0$. For the second, they combine to $1(0.01)+4(-0.02)=-0.07$.

Evaluating the original function at the new inputs $(1.01,1.98)$ gives the exact change:

$$f(1.01,1.98)-f(1,2) =\begin{bmatrix}0.0001\\-0.0696\end{bmatrix}.$$

The small discrepancy comes from the squared changes. For this function, expanding the two squares gives:

$$\Delta\mathbf f =J_f(\mathbf x)\Delta\mathbf x +\begin{bmatrix}(\Delta x_1)^2\\(\Delta x_2)^2\end{bmatrix}.$$

The Jacobian captures the first-order terms. It does not describe every finite change exactly.

Why diagrams zoom into a small region

A nonlinear transformation can bend an entire grid. Near a chosen point, its behavior can be approximated by a linear transformation of small displacements. The full approximation is:

$$f(\mathbf x+\Delta\mathbf x) \approx f(\mathbf x)+J_f(\mathbf x)\Delta\mathbf x.$$

Notice the two parts: $f(\mathbf x)$ locates the transformed center, and the Jacobian transforms the displacement from that center. Multiplying an arbitrary point by the Jacobian alone does not reproduce the original nonlinear function.

The explorer below compares the actual transformed grid with the Jacobian’s prediction. Both plots use coordinates relative to their respective centers. Reduce Neighborhood radius to see the prediction improve; the axes zoom with the radius so the shapes remain visible.

A nonlinear map and its local linear approximation

Polynomial transformationSine transformationRank 1 point

Current x₁

1

Current x₂

2

Neighborhood radius

0.5

f(x₁, x₂) = (x₁² + x₂, x₁ + x₂²)

Near the input point

Changes from the output center

Blue: actual gridDashed orange: Jacobian prediction

J = 2114

det J = 7
Rank = 2

A curved grid and its local matrix

The explorer’s Sine transformation uses:

$$g(x,y)=\begin{bmatrix}x+\sin y\\y+\sin x\end{bmatrix}, \qquad J_g(x,y)=\begin{bmatrix}1&\cos y\\\cos x&1\end{bmatrix}.$$

Here $x$ and $y$ name the two input coordinates. At $(-2,1)$:

$$J_g(-2,1)\approx\begin{bmatrix}1&0.5403\\-0.4161&1\end{bmatrix}.$$

A small horizontal displacement $(0.1,0)$ therefore produces an approximate output displacement:

$$J_g(-2,1)\begin{bmatrix}0.1\\0\end{bmatrix} \approx\begin{bmatrix}0.1\\-0.0416\end{bmatrix}.$$

Moving right in the input space also moves downward in the output space near this point. The first column of the Jacobian describes that response. The second column describes the response to moving vertically.

This is why a diagram may highlight just one small rectangle within a curved grid: it is comparing the nonlinear map with its local linear approximation at that rectangle’s center.

Shape counts variables; rank counts independent responses

The number of inputs and outputs determines the Jacobian’s shape. Its rank measures how many independent directions of output change its linear approximation can produce:

$$\operatorname{rank}(J_f)\leq\min(m,n).$$

More precisely, these first-order output changes lie in the Jacobian’s column space, whose dimension is its rank. Two inputs do not necessarily provide two independent output directions: their effects might point along the same line.

Return to our polynomial example:

$$J_f(\mathbf x)=\begin{bmatrix}2x_1&1\\1&2x_2\end{bmatrix}, \qquad \det J_f=4x_1x_2-1.$$

For this $2\times2$ matrix, a nonzero determinant means rank $2$. When $4x_1x_2=1$, its rank is $1$. It cannot have rank $0$ because its off-diagonal entries are always $1$.

Current input	Jacobian	Rank
$(1,2)$	$\begin{bmatrix}2&1\\1&4\end{bmatrix}$	$2$
$(0.5,0.5)$	$\begin{bmatrix}1&1\\1&1\end{bmatrix}$	$1$

At $(0.5,0.5)$, both columns are identical. Changing either input alone produces the same first-order output direction. A change $(h,-h)$ cancels in the linear approximation:

$$\begin{bmatrix}1&1\\1&1\end{bmatrix} \begin{bmatrix}h\\-h\end{bmatrix} =\begin{bmatrix}0\\0\end{bmatrix}.$$

This direction lies in the Jacobian’s null space. But the exact output change is $(h^2,h^2)$, so zero first-order change does not imply zero finite change.

Select Rank 1 point in the explorer. The Jacobian maps the small square to a line, while the actual transformed grid still has curvature. The rank describes the linear approximation at this point; it does not establish that an entire neighborhood of the nonlinear function collapses to a line.

This connects Jacobians to matrix rank: the linear-algebra definition stays the same, but the matrix itself can change as we move through the input space.

The chain rule becomes matrix multiplication

Now connect several transformations:

$$\mathbf a=f_1(x),\qquad \mathbf b=f_2(\mathbf a),\qquad p=f_3(\mathbf b).$$

Suppose $x$ and $p$ are scalars, while $\mathbf a$ and $\mathbf b$ each contain two values. The diagram below connects each local partial derivative to the input–output connection it describes. Select an entry in a Jacobian to highlight that connection.

One connection, one local partial derivative

J₁: 2 × 1x → a₁, a₂

out ↓ / in →	x
a₁	∂a₁∂x
a₂	∂a₂∂x

J₂: 2 × 2a₁, a₂ → b₁, b₂

out ↓ / in →	a₁	a₂
b₁	∂b₁∂a₁	∂b₁∂a₂
b₂	∂b₂∂a₁	∂b₂∂a₂

J₃: 1 × 2b₁, b₂ → p

out ↓ / in →	b₁	b₂
p	∂p∂b₁	∂p∂b₂

J₂, row 1, column 2: a₂ → b₁How quickly b₁ changes as we change a₂, holding a₁ fixed.

dp/dx = J₃ J₂ J₁
(1 × 2)(2 × 2)(2 × 1) = (1 × 1)

The number of outputs determines the rows, and the number of inputs determines the columns:

Transformation	Jacobian	Shape
$x\to\mathbf a$	$J_1$	$2\times1$
$\mathbf a\to\mathbf b$	$J_2$	$2\times2$
$\mathbf b\to p$	$J_3$	$1\times2$

A small change passes through the layers as:

$$\Delta\mathbf a\approx J_1\Delta x,\qquad \Delta\mathbf b\approx J_2\Delta\mathbf a,\qquad \Delta p\approx J_3\Delta\mathbf b.$$

Combining these gives the matrix chain rule:

$$\frac{dp}{dx}=J_3J_2J_1, \qquad (1\times2)(2\times2)(2\times1)=(1\times1).$$

Each Jacobian is evaluated at the input to its own transformation: $J_1(x)$, $J_2(\mathbf a)$, and $J_3(\mathbf b)$. The intermediate values come from the forward pass. The chain rule gives an exact derivative where these functions are differentiable, even though predicting a finite change with it is generally an approximation.

Multiplication adds the paths automatically

Use a concrete example:

$$\mathbf a=\begin{bmatrix}x^2\\x\end{bmatrix},\qquad \mathbf b=\begin{bmatrix}a_1+a_2\\2a_1\end{bmatrix},\qquad p=b_1+3b_2.$$

At $x=2$, the forward pass gives $\mathbf a=(4,2)$, $\mathbf b=(6,8)$, and $p=30$. The local Jacobians are:

$$J_1=\begin{bmatrix}4\\1\end{bmatrix},\qquad J_2=\begin{bmatrix}1&1\\2&0\end{bmatrix},\qquad J_3=\begin{bmatrix}1&3\end{bmatrix}.$$

Therefore:

$$\frac{dp}{dx} =\begin{bmatrix}1&3\end{bmatrix} \begin{bmatrix}1&1\\2&0\end{bmatrix} \begin{bmatrix}4\\1\end{bmatrix} =\begin{bmatrix}7&1\end{bmatrix} \begin{bmatrix}4\\1\end{bmatrix} =29.$$

The $7$ collects the two paths from $a_1$ to $p$: one contributes $1$, and the other contributes $3\times2=6$. Matrix multiplication performs the same multiply-along-paths, add-between-paths calculation as our scalar engine. Expanding $p=7x^2+x$ confirms that $dp/dx=14x+1=29$ at $x=2$.

Backpropagation needs Jacobian products

During training, we usually have one scalar loss and many parameters. For an operation $\mathbf y=F(\mathbf x)$, suppose the backward pass has already calculated how the final loss responds to each output $y_j$.

A layer’s input can affect several outputs, each of which can affect the final loss. To calculate that input’s gradient, backpropagation must account for every one of these routes. The Jacobian organizes the local derivatives needed for this calculation: one column records how a particular input affects all the outputs. We multiply those entries by the corresponding output gradients and add the contributions:

$$\frac{\partial L}{\partial x_i} =\sum_j\frac{\partial L}{\partial y_j} \frac{\partial y_j}{\partial x_i}.$$

Writing the gradients as columns turns this into:

$$\boxed{\nabla_{\mathbf x}L=J_F(\mathbf x)^{\mathsf T}\nabla_{\mathbf y}L.}$$

If $J_F$ has shape $m\times n$, its transpose has shape $n\times m$. Multiplying it by the $m$ output gradients produces the $n$ input gradients. The transpose does not reverse the forward function or invert its matrix; it organizes how output sensitivities combine for each input.

For the numerical example above, treat $p$ as the scalar objective and start with $dp/dp=1$. Because $p=b_1+3b_2$, the incoming gradient for the middle layer is:

$$\nabla_{\mathbf b}p=J_3^{\mathsf T}\cdot1=\begin{bmatrix}1\\3\end{bmatrix}.$$

That layer computes $b_1=a_1+a_2$ and $b_2=2a_1$. The input $a_1$ affects both outputs, so its gradient must include both routes. The input $a_2$ affects only $b_1$:

Input	Contribution through $b_1$	Contribution through $b_2$	Total input gradient
$a_1$	$1\times1=1$	$3\times2=6$	$7$
$a_2$	$1\times1=1$	$3\times0=0$	$1$

Each contribution multiplies an incoming output gradient by a local derivative. The Jacobian-transpose product performs both rows of this calculation together:

$$\nabla_{\mathbf a}p =J_2^{\mathsf T}\nabla_{\mathbf b}p =\begin{bmatrix}1&2\\1&0\end{bmatrix} \begin{bmatrix}1\\3\end{bmatrix} =\begin{bmatrix}7\\1\end{bmatrix}.$$

Finally, pass these gradients through $\mathbf a=(x^2,x)$ at $x=2$:

$$\frac{dp}{dx}=J_1^{\mathsf T}\begin{bmatrix}7\\1\end{bmatrix} =4\times7+1\times1=29.$$

This is why Jacobians appear in explanations of backpropagation: they organize all the local derivative contributions needed to move a gradient through a vector operation. Our scalar engine accumulated the same contributions one edge at a time.

Why engines don’t usually build the matrices

An operation can compute the required product from its derivative rule and saved values. It need not materialize every entry of its Jacobian. PyTorch describes this as a vector–Jacobian product. In row-vector notation it is $\mathbf v^{\mathsf T}J$; our column-vector notation writes its transpose, $J^{\mathsf T}\mathbf v$.

Consider an elementwise operation $y_i=\tanh(x_i)$. Its Jacobian is diagonal, with entries $1-y_i^2$. Given an incoming gradient $g_i=\partial L/\partial y_i$, the backward rule is simply:

$$\frac{\partial L}{\partial x_i}=g_i(1-y_i^2).$$

A vector of length $n$ needs $n$ such multiplications; there is no reason to allocate a dense $n\times n$ matrix full of off-diagonal zeros.

For a dense layer $\mathbf y=W\mathbf x+\mathbf b$, the input Jacobian is $W$. With incoming column gradient $\mathbf g$, its backward rules are:

$$\nabla_{\mathbf x}L=W^{\mathsf T}\mathbf g,\qquad \nabla_W L=\mathbf g\mathbf x^{\mathsf T},\qquad \nabla_{\mathbf b}L=\mathbf g.$$

These formulas describe one example; contributions from multiple examples are combined according to the loss’s reduction. They also show why input gradients and parameter gradients are distinct: the same operation has derivatives with respect to each of its inputs, including its weights and biases.

Input changes forward, loss sensitivities backward

The two directions use the same local derivative information for different questions:

Calculation	Question
$J\,\Delta\mathbf x$	How does a small input displacement affect the outputs?
$J^{\mathsf T}\nabla_{\mathbf y}L$	How do output sensitivities to the loss combine into input sensitivities?

Propagating a tangent vector forward computes a Jacobian–vector product. Propagating a loss gradient backward computes the corresponding vector–Jacobian product. Neither is the ordinary forward pass, which evaluates the function’s values.

Our earlier engine implements the backward calculation through scalar edges. A tensor engine groups many of those contributions into operations on arrays. The Jacobian connects the two descriptions: its columns explain local changes geometrically, its rank describes independent first-order responses, and its products express the chain-rule calculation used in backpropagation.

Follow-up: using Jacobians to interpret a language model

Sida Liu’s My Attempt to Explain the Jacobian Lens: A Single Jacobian Already Works applies these ideas to model interpretability. It follows how $Jv$ maps a small change along an intermediate activation direction to a change in the final activation, and how $J^{\mathsf T}w$ traces an output direction’s sensitivity back to the intermediate layer. The note also tests whether a Jacobian measured on one passage can help read activations from unrelated passages—a useful next step from understanding local derivatives to investigating what a trained model represents.