Jacobians, Hessians and Gauss–Newton
Most of the problems in this site's other guides — bundle adjustment, pose graphs, sensor calibration — are the same shape: minimise a sum of squared residuals that depends nonlinearly on parameters. You cannot solve them with the linear algebra of Part 10 directly, because the residual is not a line. What you can do is linearise the residual around the current guess, solve a linear least-squares problem for the correction, and repeat. This part is that loop, and the matrix that drives it.
The problem
A sum of squares with no closed-form answer
Take a model with a handful of parameters p and a stack of measurements. For each measurement i define a residual ri(p): how far the model's prediction is from what you observed. The cost is the sum of squares,
When each ri is linear in p, this is exactly the least-squares problem of Part 10, and the normal equations give the answer in one step. When the residuals are nonlinear — a camera projection, an exponential decay, a rotation — there is no formula, only a search. The search is built out of linear solves, which is why this part is linear algebra.
Linearise, then solve
One Taylor step, one least-squares problem
Suppose you have a current guess p and want a correction δ. A first-order Taylor expansion says r(p+δ) ≈ r(p) + Jδ, where J is the Jacobian: one row per residual, one column per parameter, entry Jij = ∂ri/∂pj. Minimising the linearised cost over δ is now a linear least-squares problem,
That is the Gauss–Newton step: form JᵀJ and Jᵀr, solve a linear system, add δ to p, and do it again. The matrix JᵀJ is square, symmetric and positive semidefinite, so it is the same kind of object as Part 15's quadratic forms. Watch a single step on a decaying-exponential fit: the model snaps toward the data, and the residuals shrink.
Blue points are the data; the curve is a·e−bx + c; the grey segments are the residuals ri. One step solves JᵀJ δ = −Jᵀr.
Each step is a linear solve, not a search. The nonlinearity is handled entirely by re-forming J and r at the new p and stepping again, which is why the method is called iterative: it walks a sequence of quadratic approximations toward a minimum.
The normal matrix JᵀJ
What the Jacobian's shape buys you
The Jacobian carries all the local information. Its rows say how sensitive each residual is to each parameter; inside JᵀJ those rows become inner products, exactly as in Part 7. A parameter direction that no residual responds to gives a nearly zero diagonal entry — a redundant parameter, and an ill-conditioned normal matrix. The heatmap and the eigenvalues below update as you step the fit above.
The Jacobian J: 12 residuals (rows) × 3 parameters (columns).
Because JᵀJ is symmetric, its eigenvalues are real and its eigenvectors are orthogonal — Part 15 again. They are the principal curvatures of the local quadratic, so their ratio is the condition number of the linear solve at the heart of every Gauss–Newton step. That ratio is the practical version of the conditioning story in Part 19, and it is why real solvers add a damping term before inverting.
The Hessian, and why JᵀJ is enough
Dropping the second-order term on purpose
The true second derivative of the cost is the Hessian. Differentiating C = ½‖r‖² twice gives two pieces,
Gauss–Newton throws away the second sum. That sounds reckless until you notice what it buys: JᵀJ is guaranteed symmetric positive semidefinite and costs nothing extra to form, while the discarded term needs every residual's second derivative and can make H indefinite. Near a good fit the residuals are small, so the dropped term is small too, and the approximation is excellent; far from a fit it can be poor, which is what trust-region and Levenberg–Marquardt machinery exists to manage.
Gauss–Newton
Solve JᵀJ δ = −Jᵀr. Quadratic convergence near a good fit, cheap, but can overshoot when the linearisation is bad.
Levenberg–Marquardt
Solve (JᵀJ + λI) δ = −Jᵀr. The λ term interpolates between Gauss–Newton and a small gradient step, keeping the step inside the region where the linearisation holds. This is the same ridge idea as the damped pseudoinverse in Part 18.
Sparsity is the whole game
Why pose graphs with a million variables are solvable
A bundle-adjustment problem has one small block of parameters per camera and per 3D point, and each measurement touches only two of them. So a Jacobian row has a handful of nonzeros, and JᵀJ is sparse in the same pattern: a block for every parameter, and an off-diagonal block only where two parameters appear in a common measurement. A chain of poses produces the block-tridiagonal pattern below. The matrix is enormous and nearly empty, which is what makes the linear solve affordable.
The nonzero pattern of JᵀJ for a chain of N poses. One node per parameter block; an edge means the two appear together in a measurement.
The picture generalises directly. Replace each scalar node with a 6-DOF pose block and you have the structure of a SLAM problem; the nonzero count stays proportional to the number of measurements rather than the square of the number of parameters, and sparse Cholesky or QR exploits that to solve systems that would be hopeless as dense matrices. Knowing the pattern is knowing which factorisation you can afford, which is the practical content of Part 19.
Where this shows up
The loop behind every geometry pipeline
Pose graphs and bundle adjustment
The pose-graph part of the optimization guide is exactly this loop with SE(3) blocks, and structure from motion runs it over camera poses and 3D points at once. The sparsity you just saw is the only reason those systems close.
MAP estimation and second-order training
Fitting a probabilistic model by maximum a posteriori estimation is the same least-squares problem, and second-order optimisers approximate the curvature with JᵀJ. The architecture chapter shows the forward pass these derivatives come from.
Cheat sheet
| Object | Size | Meaning |
|---|---|---|
| r(p) | m × 1 | Residual for each of the m measurements |
| J | m × n | Jij = ∂ri/∂pj, one row per measurement |
| JᵀJ | n × n | Symmetric PSD approximation of the Hessian |
| Jᵀr | n × 1 | Gradient of the cost, up to sign |
| (JᵀJ + λI)δ = −Jᵀr | n × 1 | Levenberg–Marquardt step; λ → 0 is Gauss–Newton |
Further reading
- This site's nonlinear optimization guide — the same material developed for robotics.
- K. Madsen, H. B. Nielsen and O. Tingleff, Methods for Non-Linear Least Squares Problems, IMM, 2004.
- Gilbert Strang, 18.06 Linear Algebra, Lecture 17 (least squares and the normal equations).
- Triggs et al., Bundle Adjustment — A Modern Synthesis, 2000 — where the sparsity pattern earns its keep.