Matrix Calculus: Gradients and Jacobians

Calculus on scalars asks how one output moves when one input moves. Vector calculus asks the same of many outputs and many inputs at once, and packs the answer into a matrix. That matrix, the Jacobian, is the derivative: it is the best linear map that copies the behavior of a nonlinear function near a point. Every dynamics linearization, policy gradient, and differentiable rollout is an instance of assembling and multiplying these matrices.

Gradients and Jacobians

For a scalar field $f:\mathbb{R}^n\to\mathbb{R}$ the derivative is the gradient, the column of partials $\nabla f = (\partial f/\partial x_1,\dots,\partial f/\partial x_n)^{\top}$. For a vector field $f:\mathbb{R}^n\to\mathbb{R}^m$ each output has its own gradient, and stacking them by row gives the Jacobian

Row $i$ is $\nabla f_i^{\top}$; column $j$ is $\partial f/\partial x_j$, the direction the whole output moves when input $j$ moves. Two conventions coexist. The numerator layout above sends outputs to rows, so $J\in\mathbb{R}^{m\times n}$. The denominator layout transposes this, putting inputs on rows, so a scalar gradient comes out as a column. The math is identical; only the shapes and the order of factors in the chain rule flip. Pick one and hold it fixed.¹

The Jacobian is the best local linear map

Near a point $x$ the first-order Taylor expansion is

So $J\delta$ is the linear map that best reproduces $f$ for small displacements: the derivative is not a number attached to a point, it is a matrix. Feed it the basis vector $e_j$ and you recover column $j$, confirming that columns are the responses to unit input moves.

Because $J$ is a linear map, its geometry is the SVD story. Write $J = U\Sigma V^{\top}$. A small sphere of perturbations $\lVert\delta\rVert = \varepsilon$ around $x$ is sent to an ellipsoid around $f(x)$ whose semi-axes are $\varepsilon\,\sigma_i u_i$: the singular values are the local stretch factors and $\det J = \prod_i \sigma_i$ is the local volume scale.¹ The right singular vectors $v_i$ mark the input directions that stretch most and least; each maps to $\sigma_i u_i$.

The chain rule is matrix multiplication

Compose $h = f\circ g$. The Jacobian of the composition is the product of the Jacobians, in output-to-input order,

This one identity is the whole of backpropagation and of gradient flow through a simulator: a deep network or a $T$-step rollout is a long composition, and its Jacobian is the product of the per-layer or per-step Jacobians.² Forward mode multiplies this product right to left onto an input perturbation; reverse mode multiplies left to right onto an output cotangent, which is why reverse mode is cheap when there are few outputs, such as a scalar loss.

A worked 2×2 example

Take $f(x,y) = (x^2 - y,\; xy)$. Differentiating each component,

Verify column one against a central finite difference with $h = 10^{-3}$. The column is $\partial f/\partial x$, so perturb only $x$ at the point $(1,2)$:

The finite difference reproduces the analytic column to the digits shown, with error $O(h^2)$. This is exactly the check the demo runs live for every point.

Where these Jacobians show up

Linearizing continuous dynamics $\dot{x} = f(x,u)$ gives the pair $A = \partial\dot{x}/\partial x$ and $B = \partial\dot{x}/\partial u$; LQR and iLQR consume exactly these Jacobians at each nominal point to build their local quadratic models. In reinforcement learning the policy gradient is a Jacobian of sampled returns with respect to parameters, and value-function fitting rides the Jacobian of the critic. Differentiable simulation multiplies per-step Jacobians through a rollout so a loss can be pushed back onto controls, masses, or contact parameters.² In every case the object being assembled is the same matrix of partials, and the chain rule is the same matrix product.