Jacobians & Linearization

Chapter 1: The Best Linear Approximation

Let's build the one dimensional version of the fix properly, with real numbers, before we generalize to many inputs. Take a smooth function f(x). Pick a point a where you know the value f(a) and the slope f′(a). The claim is that near a, the straight line through (a, f(a)) with slope f′(a) is the best possible linear stand in for f.

Read this formula as a sentence. Start at the known height f(a). Walk a distance of (x − a) away from a. For each unit you walk, climb by f′(a), the slope at the starting point. That's the whole idea: start where you know the answer, then extrapolate using the one thing you know about how it's changing.

A hand computed worked example

Take f(x) = √x, and linearize it at a = 4. First, f(4) = √4 = 2. Next, the derivative of √x is 1 / (2√x), so f′(4) = 1 / (2·2) = 0.25. The linear approximation is:

Now push this line past its comfort zone and watch it drift from the truth. The table below computes the true value, the linear guess, and the error at four increasing distances from a = 4.

Look at the error column. Doubling δ from 0.5 to 1.0 roughly quadruples the error, from about 0.0037 to about 0.0139. Doubling again to δ = 2.0 roughly quadruples it again, to about 0.0505. This isn't a coincidence, and it isn't specific to square roots. It's the signature of a squared term.

Code: the same computation

python
import math

def linearize(f, fprime, a, x):
    # L(x) = f(a) + f'(a) * (x - a)
    return f(a) + fprime(a) * (x - a)

f      = math.sqrt
fprime = lambda x: 1 / (2 * math.sqrt(x))

for x in [4.5, 5.0, 6.0, 8.0]:
    true_val = f(x)
    approx   = linearize(f, fprime, 4.0, x)
    print(f"x={x}  true={true_val:.6f}  approx={approx:.6f}  err={abs(true_val-approx):.6f}")

Chapter 2: Partials and the Gradient

The tangent line trick from Chapter 1 handles a function with one input and one output. Our range sensor has two inputs, x and y. We need to generalize "the slope at a point" to a setting where there's more than one direction to walk in.

The tool is the partial derivative, and the rule for building one is refreshingly simple: freeze every variable except the one you're differentiating with respect to, and treat the frozen ones as plain constants. For f(x, y), the partial with respect to x, written ∂f/∂x or f_x, holds y fixed and asks how f changes as x alone wiggles. The partial with respect to y does the mirror image.

Stack the two partials into a vector and you get the gradient, written ∇f:

The gradient is a compass. It points in the direction that increases f the fastest, and its length tells you how steep that fastest direction is. Just as one derivative gave us a tangent line in Chapter 1, the gradient gives us a tangent plane in two dimensions: the flattest possible sheet that still touches the surface at one point and leans in exactly the direction the surface leans.

A hand computed worked example

Take f(x, y) = x² + xy, and linearize it at the point (1, 2).

Step 1, the partials. Freeze y, differentiate in x:

Freeze x, differentiate in y:

Step 2, evaluate at (1, 2).

So ∇f(1, 2) = (4, 1), and the tangent plane is:

Step 3, test it nearby. At (1.1, 2.1):

An error of just 0.02 for a point that moved a full 0.1 units in each direction. That's the tangent plane earning its keep: close to (1, 2), it's an excellent stand in for the true curved surface.

Drag the sliders around the bowl. Notice that wherever you land, the warm gradient arrow always points straight out, perpendicular to the ring you're standing on, toward higher ground. The rings are lines of equal height, so walking along a ring changes nothing, all the change is in the perpendicular direction. That's not a coincidence of this particular bowl, it's true of every smooth surface, and it's why gradient descent, the workhorse of machine learning training, always steps perpendicular to the current loss contour.

Chapter 3: The Jacobian

The gradient handled a function with many inputs and one output. Our range and bearing sensor has two inputs and, once we add bearing to the mix, two outputs at once. We need one more generalization: a function f: ℝⁿ → ℝ^m, n numbers in, m numbers out. The tool for linearizing it is the Jacobian matrix, and it's built from a very natural idea: stack a gradient for every output.

Row i of J is exactly the gradient of output f_i, the i-th output's own private compass. So J is an m by n matrix: m rows, one per output, n columns, one per input. Entry J_ij = ∂f_i/∂x_j answers a precise question: if I wiggle input j by a tiny amount, how much does output i move? That's it, that's the entire object. Everything else in this lesson is putting that one idea to work.

Hand computing the polar to Cartesian Jacobian

Here is the single most useful worked example in robotics: converting a polar coordinate, a distance r and an angle θ, into a Cartesian point (x, y). The map is:

Two outputs, x and y. Two inputs, r and θ. The Jacobian is 2 by 2. Let's fill in all four entries with the freeze the other variable rule from Chapter 2, one at a time, showing every step.

Now plug in real numbers. Take r = 4, θ = 60° (that's π/3 radians), so cosθ = 0.5 and sinθ ≈ 0.866025.

Read the first column, [0.5, 0.866], as an instruction: increasing r by one unit moves the point 0.5 units in x and 0.866 units in y, exactly the unit direction the point is already sitting in. Read the second column, [−3.4641, 2.0], as a different instruction: increasing θ by one radian sweeps the point sideways, perpendicular to the first column, and four times farther, because sweeping a lever of length r = 4 through an angle covers 4 times as much arc length as sweeping a lever of length 1.

Slide r up. The ellipse stretches wider in the direction perpendicular to the line to the origin, exactly as the "sweeping a longer lever covers more distance" argument predicts. Slide θ around a full circle instead, and the ellipse just spins in place, since neither cosθ nor sinθ changes the size of the stretch, only its direction. This picture, a small circle of uncertainty becoming an ellipse under a Jacobian, is exactly the picture we'll use in Chapter 8 to compare true and linearized sensor noise.

Chapter 4: The Chain Rule, in Matrix Form

Real systems are rarely one function. A lidar return goes through a polar to Cartesian conversion, then a rotation into the robot's body frame, then a translation into a world map frame, three functions composed end to end. A neural network is dozens of layers composed the same way. How does linearization survive being composed?

Beautifully, it turns out. If y = g(f(x)), meaning you feed x into f, then feed the result into g, the single variable chain rule you already know says the derivative of the composition is the product of the two derivatives: dy/dx = g′(f(x))·f′(x). The matrix version says exactly the same thing, with matrix multiplication standing in for ordinary multiplication:

Read this carefully, because the evaluation points matter. J_f is evaluated at x, the original input. J_g is evaluated at f(x), the output of the first stage, since that's the input g actually sees. Multiply those two matrices together, in that order, and you get the Jacobian of the whole pipeline, no matter how curved each individual stage is.

A hand computed two stage example

Reuse Chapter 3's polar to Cartesian map as stage one: u = f(r, θ) = (r·cosθ, r·sinθ), evaluated at r = 4, θ = 60°, giving u = (2, 3.4641) and J_f = [[0.5, −3.4641], [0.8660, 2.0]] exactly as before.

Now bolt on a second, deliberately nonlinear stage: g(u₁, u₂) = (u₁² + u₂,   u₁·u₂). Its Jacobian, by the same freeze and differentiate rule as always:

Evaluated at u = (2, 3.4641): J_g = [[4, 1], [3.4641, 2]]. Now multiply, row by column, tracking every number:

Four numbers, and they cost nothing more than two smaller Jacobians and one matrix multiply, no matter how the two stages were defined. This is exactly the computation that makes deep learning practical: backpropagation is nothing more than this chain rule, applied layer by layer, from the output back to the input, multiplying local Jacobians as it goes. See the Backpropagation lesson for that story told end to end.

python
import numpy as np

def Jf(r, theta):
    return np.array([[np.cos(theta), -r*np.sin(theta)],
                     [np.sin(theta),  r*np.cos(theta)]])

def Jg(u1, u2):
    return np.array([[2*u1, 1],
                     [u2,      u1]])

r, theta = 4.0, np.pi/3
u1, u2   = r*np.cos(theta), r*np.sin(theta)     # = f(r, theta)
J_total  = Jg(u1, u2) @ Jf(r, theta)             # chain rule: Jg(f(x)) . Jf(x)
print(np.round(J_total, 3))

Chapter 5: Numerical vs Analytic Jacobians

Every Jacobian so far came from algebra: freeze a variable, differentiate, write down the formula. That's the analytic route, and it's exact. But analytic Jacobians are also where real EKF code most often breaks: a dropped minus sign, a forgotten chain rule term buried three functions deep, a copy paste mistake between the x row and the y row. You need a second, independent way to compute the same numbers, one that never touches your formula at all.

That's a finite difference. Instead of differentiating symbolically, nudge one input by a tiny step h, see how much the output moved, and divide.

The central version costs one extra function call per input but is dramatically more accurate for the same h, because it cancels the leading error term that the forward version leaves behind. This is worth understanding precisely, since the size of h is not a free choice, it's a genuine engineering trade off.

The trade off: too big, too small

Make h too large, and the straight line secant you're drawing cuts across real curvature in f, so you're not really measuring the slope at x anymore, that's truncation error, and it shrinks as h shrinks. Make h too small, and you run into a completely different problem: computers store numbers with finite precision, roughly 16 significant digits for a standard float. Subtracting two nearly identical numbers, f(x+h) and f(x), throws away almost all of those digits, and dividing by a tiny h then blows up whatever tiny error was left. That's round off error, and it grows as h shrinks. Somewhere in between sits a sweet spot.

Drag h from the right side of the plot toward the left. At first, both curves fall as h shrinks, truncation error dominates and smaller steps help. Keep going and the central difference curve, which starts with much smaller error thanks to canceling that leading term, bottoms out and starts climbing again as round off takes over. The forward difference curve does the same thing, but its floor sits many orders of magnitude higher. Somewhere around h = 10⁻⁵ to 10⁻⁴ is usually a safe, general purpose choice for central differences in double precision arithmetic.

The same computation, three ways

Let's put this into practice on the polar to Cartesian Jacobian from Chapter 3, computed three independent ways: by hand coded analytic formulas, by a from scratch central finite difference in NumPy, and by automatic differentiation in PyTorch, which computes exact derivatives without you writing any formula at all. All three should land on the same numbers we computed by hand: [[0.5, −3.4641], [0.8660, 2.0]].

python · plain, analytic formulas
import math

def polar_to_cart(r, theta):
    return (r * math.cos(theta), r * math.sin(theta))

def analytic_jacobian(r, theta):
    c, s = math.cos(theta), math.sin(theta)
    return [[c,  -r * s],
            [s,   r * c]]

r, theta = 4.0, math.pi / 3
print("analytic:", analytic_jacobian(r, theta))
# [[0.5, -3.4641], [0.866, 2.0]]

numpy · central finite differences
import numpy as np

def polar_to_cart(rt):
    r, theta = rt
    return np.array([r * np.cos(theta), r * np.sin(theta)])

def jacobian_fd(f, x, h=1e-5):
    x = np.asarray(x, dtype=float)
    n = len(x)
    J = np.zeros((len(f(x)), n))
    for j in range(n):
        step = np.zeros(n); step[j] = h
        J[:, j] = (f(x + step) - f(x - step)) / (2 * h)   # central diff
    return J

x0 = np.array([4.0, np.pi / 3])
print("finite-diff:", np.round(jacobian_fd(polar_to_cart, x0), 4))

pytorch · automatic differentiation
import torch

def polar_to_cart(rt):
    r, theta = rt[0], rt[1]
    return torch.stack([r * torch.cos(theta), r * torch.sin(theta)])

x0 = torch.tensor([4.0, math.pi / 3])
J  = torch.autograd.functional.jacobian(polar_to_cart, x0)
print("autograd:", J)

Three completely different code paths, none of which knows about the other two, and all three land on [[0.5, −3.4641], [0.8660, 2.0]] to within rounding. That agreement is exactly the check a real filter's test suite runs before the code ever ships.

Chapter 6: Linearizing a Measurement Model

Time to solve the actual problem from Chapter 0. Give the tracking station a real sensor: it reports range and bearing to the drone, not just range. If the drone sits at (x, y) with the dish at the origin, the measurement function is:

The first output is range, the second is bearing, the angle from the dish to the drone. This is the Cartesian to polar map, the inverse direction from Chapter 3's worked example. In an Extended Kalman Filter this h is what plugs into the update step, and the matrix we need for the covariance math is its Jacobian, usually written H instead of J when it's specifically the measurement Jacobian:

Computing H by hand

Freeze y, differentiate range with respect to x: range = (x²+y²)^1/2, so ∂range/∂x = x·(x²+y²)^−1/2 = x / range. By the same logic, ∂range/∂y = y / range. The bearing partials take one more step of chain rule through arctangent, but the result is clean:

Plug in a concrete drone position: (x, y) = (3, 4), so r = √(9+16) = √25 = 5.

A remarkable hidden structure

Look closely at the top row, (0.6, 0.8). Its length is √(0.6² + 0.8²) = √(0.36+0.64) = √1 = 1, always, for any (x, y), because that row is just (cosθ, sinθ) for the bearing angle θ, a unit vector by definition. Multiply H by its own transpose, H·H^T, and the cross terms between the two rows cancel:

Check it with the numbers above: row one dotted with itself is 0.6² + 0.8² = 1. Row one dotted with row two is 0.6·(−0.16) + 0.8·0.12 = −0.096 + 0.096 = 0, they're exactly perpendicular. Row two dotted with itself is (−0.16)² + 0.12² = 0.0256 + 0.0144 = 0.04 = 1/25 = 1/r². This diagonal, orthogonal structure means range sensitivity and bearing sensitivity point in completely independent directions, and it's the exact formula we'll reuse in Chapter 8 to build the linearized uncertainty ellipse without any extra machinery.

Chapter 7: Linearizing Dynamics

Sensors aren't the only nonlinear piece. A robot's own motion model usually is too. Consider a simple ground robot with state (x, y, θ), position plus heading, driven forward at speed v with turn rate ω for a short time step dt:

Position updates depend on the heading through sine and cosine, curved functions of θ, so this motion model needs the same linearization treatment as any sensor. The Jacobian here is usually called F, the process Jacobian, and because there are three state variables, F is 3 by 3.

Computing F by hand

Two of the three columns are almost free. Differentiating the x update with respect to x itself gives 1, and with respect to y gives 0, since x doesn't appear in the y update at all, and vice versa. The only genuinely nonlinear column is the θ column, since that's the only place sine and cosine show up:

Plug in v = 2 m/s, θ = 30° (sinθ = 0.5, cosθ ≈ 0.8660), dt = 0.1 s: −v·sinθ·dt = −2·0.5·0.1 = −0.1, and v·cosθ·dt = 2·0.8660·0.1 ≈ 0.1732.

Why we need F: propagating the covariance

The state estimate itself gets pushed through the true, curved f. But the uncertainty around that estimate, the covariance matrix P, gets pushed through the linearized F instead, using the sandwich formula:

Let's do this by hand too, with a starting covariance P = diag(0.04, 0.04, 0.01), meaning position uncertainty of about 0.2 units in x and y and heading uncertainty of about 0.1 radians, all currently uncorrelated with each other. First, F·P, easy since P is diagonal, it just scales F's columns:

Now multiply that result by F^T, whose columns are F's rows:

Chapter 8: The Linearization Lab

Everything in this lesson lands here. On the left, a small tracking dish sits at the origin and a drone sits somewhere nearby, the operating point where we'll build the tangent plane. Drag the drone anywhere you like. Around it hangs a cloud of sample points, real position uncertainty, its spread set by the sigma slider.

On the right, each of those sample points is pushed through the true, curved range and bearing function h(x, y) from Chapter 6, plotted as small dots in range bearing space, the true distribution. On top of that cloud sits an ellipse, drawn from the closed form linearized covariance H·Σ·H^T we derived in that same chapter, the linearized approximation.

Start with the drone far from the dish, several units out. The cloud and the ellipse sit almost exactly on top of each other; the linear approximation is excellent here, exactly what Chapter 1's error grows with δ² result predicts when your uncertainty σ is small relative to how fast the curvature changes. Now click "Very close" and watch the mismatch appear. The true cloud bends into a curved, banana shaped smear as bearing wraps around near the dish, while the linearized ellipse, forced by construction to be a perfect axis aligned oval, cannot bend at all. It captures the right overall size but completely misses the true shape.

Try "Straight ahead," drone directly in front of the dish along the y axis. Even at a moderate range, notice the true cloud is slightly skewed relative to the symmetric ellipse, a small foretaste of the same nonlinearity that becomes dramatic up close. There is no quiz for this chapter. If you can predict, before dragging, whether the ellipse will hug the cloud or badly miss it, purely from how close the drone is to the dish, you understand linearization.

Chapter 9: Choose Your Tool

Four objects showed up across this lesson, and it's easy to blur them together. They're related, but each answers a different question and each has a different shape. Here's the map.

Quick cheat sheet

Where to go next

• Extended Kalman Filter: the F and H you built here plug directly into its predict and update equations. If you understood Chapters 6 and 7, the EKF's linearization step will feel obvious rather than magical.
• Unscented Kalman Filter: the direct answer to Chapter 8's showcase. Instead of one linearized ellipse, it propagates a small deterministic set of sample points through the true nonlinear function and rebuilds the covariance from where they land.
• The Hessian & Curvature: the natural next step once linear approximations aren't enough. A Hessian is a Jacobian of a gradient, second derivatives instead of first.
• Backpropagation: Chapter 4's chain rule, applied to a network with dozens of layers instead of two, computed efficiently by multiplying local Jacobians from the output back toward the input.

"What I cannot create, I do not understand." Richard Feynman.
You hand computed a Jacobian, checked it two independent ways, linearized a real sensor, watched it fail exactly where the geometry predicted, and built the lab that compares truth to approximation. You can now create it.