13. Jacobians and Hessians: derivatives as matrices

📖 Reading · 10 min
💡 Every code box below is live — edit it and hit Run.

§10.5's gradient packaged a scalar function's partial derivatives into a vector. This lesson packages derivatives into full matrices — the Jacobian, for functions with several outputs, and the Hessian, for a scalar function's second derivatives — and shows both are already familiar objects from earlier in the course, just organized and named for their role in optimization.

Prerequisite note. Like §13.8, this lesson uses linear algebra the calculus series does not develop: matrix–vector multiplication, matrix symmetry, the 2\times2 determinant, and the vocabulary of positive/negative definite, indefinite and positive semidefinite — including the characterization "all eigenvalues \ge0". Each term is defined where it is first used and the lesson stands on its own, but a prior acquaintance with matrices will make it read much more smoothly — the companion linear-algebra course's §16.2–§16.4 (matrix operations and determinants) and §21.2–§21.3 (positive definite matrices and quadratic forms) cover this in full.

The Jacobian matrix

For a function \vec F:\mathbb R^n\to\mathbb R^m with m output components, each depending on n input variables, \vec F(\vec x)=\langle f_1(\vec x),\ldots,f_m(\vec x)\rangle, the Jacobian matrix collects every partial derivative into one m\times n array:

J=\begin{pmatrix}\dfrac{\partial f_1}{\partial x_1}&\cdots&\dfrac{\partial f_1}{\partial x_n}\\\vdots&\ddots&\vdots\\\dfrac{\partial f_m}{\partial x_1}&\cdots&\dfrac{\partial f_m}{\partial x_n}\end{pmatrix}

— row i holds \nabla f_i, so the Jacobian is simply every output's gradient (§10.5), stacked into rows. Two special cases already appeared earlier in this course: when m=1 (a single scalar output), the Jacobian is the gradient, just written as a row instead of a column; when m=n and only the determinant of J is used, this is exactly §11.5's Jacobian determinant from the change-of-variables formula — the same matrix, put to a different use there (measuring area/volume distortion) than here (representing a linear approximation to \vec F itself).

The Jacobian is the best linear approximation to \vec F near a point — the direct multi-output generalization of §10.3's tangent plane, since $\vec F(\vec x)\approx\vec F(\vec x_0)+J(\vec x_0)(\vec x-\vec x_0)$ using ordinary matrix-vector multiplication.

The Hessian matrix

For a scalar function f:\mathbb R^n\to\mathbb R, the Hessian collects every second-order partial derivative:

H=\begin{pmatrix}f_{x_1x_1}&f_{x_1x_2}&\cdots\\f_{x_2x_1}&f_{x_2x_2}&\cdots\\\vdots&\vdots&\ddots\end{pmatrix}

The Hessian is symmetricH_{ij}=H_{ji} — a direct consequence of Clairaut's theorem (§10.2), the identical fact that made mixed partials order-independent back in Module 10.

For two variables, \det(H) is exactly §10.6's discriminant:

H=\begin{pmatrix}f_{xx}&f_{xy}\\f_{xy}&f_{yy}\end{pmatrix},\qquad\det(H)=f_{xx}f_{yy}-(f_{xy})^2=D

— §10.6's second-derivative test, restated: \det(H)>0 with f_{xx}>0 is a local minimum, \det(H)>0 with f_{xx}<0 is a local maximum, \det(H)<0 is a saddle. In the language of linear algebra, these conditions are usually phrased as H being positive definite (minimum), negative definite (maximum), or indefinite (saddle) — the same classification, in terminology that generalizes cleanly to more than two variables, where a 3\times3 or larger determinant test alone no longer suffices.

Connecting back to convexity

§14.1's convexity condition, f''(x)\ge0 for one variable, generalizes directly: a multivariable function f is convex exactly when its Hessian is positive semidefinite at every point in its domain (all eigenvalues \ge0) — the Hessian is quite literally §14.1's second derivative, promoted to a matrix, and the same payoff (every local minimum is global) transfers unchanged.

Doing it in Python

The Jacobian of a vector-valued function with two outputs and two inputs, \vec F(x,y)=\langle x^2y,\,x+y^2\rangle:

import sympy as sp

x, y = sp.symbols('x y')
F = sp.Matrix([x**2*y, x + y**2])

J = F.jacobian([x, y])
print(f"Jacobian:\n{J}")

The Hessian of f(x,y)=x^3-3xy+y^3 at (1,1) — the exact function and point from §10.6's worked example — confirming \det(H) reproduces that lesson's discriminant D=27 precisely:

import sympy as sp

x, y = sp.symbols('x y')
f = x**3 - 3*x*y + y**3

H = sp.hessian(f, [x, y])
print(f"Hessian:\n{H}")

H_at_point = H.subs({x: 1, y: 1})
print(f"Hessian at (1,1):\n{H_at_point}")
print(f"det(H) at (1,1) = {H_at_point.det()}")
print("matches D=27 from section 10.6's second-derivative test exactly")

Checking whether a Hessian is positive definite (confirming a convex function, or a genuine local minimum) using its eigenvalues:

import sympy as sp

x, y = sp.symbols('x y')
f = x**2 + y**2   # a simple bowl

H = sp.hessian(f, [x, y])
eigenvalues = H.eigenvals()
print(f"Hessian:\n{H}")
print(f"eigenvalues: {eigenvalues}")
print("both positive everywhere -- positive definite -- f is convex")

Worked example

Compute the Hessian of f(x,y)=x^3-3xy+y^3 at (1,1), and use it to confirm §10.6's classification of that point.

f_x=3x^2-3y,\qquad f_y=-3x+3y^2

f_{xx}=6x,\qquad f_{yy}=6y,\qquad f_{xy}=-3

H=\begin{pmatrix}6x&-3\\-3&6y\end{pmatrix}

At (1,1):

H(1,1)=\begin{pmatrix}6&-3\\-3&6\end{pmatrix}

\det(H)=6\cdot6-(-3)(-3)=36-9=\boxed{27}

Since \det(H)=27>0 and H_{11}=f_{xx}(1,1)=6>0:

\boxed{(1,1)\text{ is a local minimum — matching §10.6's classification exactly}}

Sanity check. This is the identical conclusion §10.6 reached using D=f_{xx}f_{yy}-(f_{xy})^2 directly, computed here as \det(H) instead — the same arithmetic, organized into matrix form. The Hessian's eigenvalues (computable directly, though not needed for this 2\times2 shortcut) would both be positive here, which is the genuinely general form of "positive definite" that extends cleanly to three or more variables, where a single 2\times2-style discriminant no longer applies. ✓

Your turn

1. Find the Jacobian of \vec F(x,y)=\langle e^{xy},\,x^2-y^2\rangle.

2. Find the Hessian of f(x,y)=x^2+4xy+y^2, and use \det(H) to classify the critical point at the origin (find the critical point first by confirming \nabla f(0,0)=\vec0).

3. True or false: the Hessian of a function of three variables is a 3\times3 matrix.

Solutions

1.

J=\begin{pmatrix}ye^{xy}&xe^{xy}\\2x&-2y\end{pmatrix}

2. f_x=2x+4y, f_y=4x+2y. At (0,0): both are 0 ✓, a critical point.

f_{xx}=2,\quad f_{yy}=2,\quad f_{xy}=4

H=\begin{pmatrix}2&4\\4&2\end{pmatrix},\qquad\det(H)=4-16=-12<0

\boxed{\text{saddle point}}

3. True. For n variables, the Hessian collects every pairwise second partial derivative f_{x_ix_j} into an n\times n array — three variables give 3\times3=9 entries (though only 6 are independent, since Clairaut's theorem forces the matrix to be symmetric, H_{ij}=H_{ji}).

Check yourself in code

Compute the Hessian of f(x,y)=x^3-3xy+y^3 at (1,1) and its determinant.

Print exactly this:

Hessian at (1,1) = Matrix([[6, -3], [-3, 6]])
det(H) = 27
import sympy as sp

x, y = sp.symbols('x y')
f = x**3 - 3*x*y + y**3

H = sp.hessian(f, [x, y])
H_at_point = H.subs({x: 1, y: 1})
print("Hessian at (1,1) = ...")
print("det(H) = ...")
import sympy as sp

x, y = sp.symbols('x y')
f = x**3 - 3*x*y + y**3

H = sp.hessian(f, [x, y])
H_at_point = H.subs({x: 1, y: 1})
print(f"Hessian at (1,1) = {H_at_point}")
print(f"det(H) = {H_at_point.det()}")

The Jacobian stacks every output's gradient (§10.5) into a matrix, supplying the best linear approximation to a vector-valued function exactly the way §10.3's tangent plane did for a scalar one, while the Hessian collects a scalar function's second partials into a symmetric matrix (via Clairaut's theorem, §10.2) whose determinant is precisely §10.6's discriminant D for two variables, and whose positive- or negative-definiteness generalizes both that test and §14.1's convexity condition to any number of variables at once.

Next: the chain rule, rebuilt in matrix form — the exact operation behind backpropagation, the algorithm that trains every neural network.