Bundle Adjustment & Structure from Motion
Every part so far used two or three cameras. Real reconstructions — a phone scan, a drone survey, a SLAM map — use dozens to thousands of views of hundreds of thousands of points. Making every camera and every point agree with every observation, all at once, is bundle adjustment: the payoff this whole series has been building toward, and the finale of this guide.
From pairs and triples to everything at once
Setup
Structure from Motion (SfM) is the general problem: given a pile of overlapping photos with unknown camera poses, recover both the 3D structure of the scene and where every camera was standing. A typical pipeline chains together everything from this series: match features and estimate F/E between pairs (Part 6), recover relative poses (Part 10), triangulate an initial point cloud (Part 10), extend to more views via point transfer (Part 12) — and then, because every step so far only used two or three views at a time, small errors accumulate. Bundle adjustment is the cleanup step that fixes that: refine every camera pose and every 3D point together, so the whole reconstruction agrees with every single observation as well as possible.
One cost function, all cameras, all points
Foundations
Stack every camera's pose and every point's position into one big set of unknowns. For every observation — camera i seeing point j at pixel x⁅₠ — define a reprojection residual, and sum the squared residuals over the whole dataset:
That's it — that's the entire bundle adjustment problem. It looks just like the robot localization cost function from the Nonlinear Optimization series — because it is the same kind of problem: a sum of squared nonlinear residuals, minimized with exactly the same machinery. Gauss-Newton and, in practice, Levenberg-Marquardt (see Part 5 of that series and its treatment of 3D rotations, since camera orientations live on the same curved rotation manifold) are the standard solvers — this is literally the textbook motivating example for why that series' machinery exists. And when the same optimization runs incrementally, one new camera pose at a time, as a robot moves — instead of once, in a batch, over a fixed photo set — it's SLAM: the online sibling of this page's offline problem.
Working example: a toy structure-from-motion pipeline
Interactive · the finale
Five cameras (fixed, known poses — as if already calibrated by the earlier parts of this series) each observe the cube's 8 corners, with pixel noise added. Corners are first initialized by linear triangulation (Part 10) from the first two cameras that see them — a fast but rough estimate. Click Run bundle adjustment to jointly refine every point's position against every camera that observed it, via Gauss-Newton, and watch the cost fall and the reconstructed cube tighten up around the truth. For a non-iterative weak-perspective alternative to this iterative refinement, see affine factorization (Part 14).
Drag to orbit, scroll to zoom. Diamonds = the 5 fixed cameras.
Total reprojection cost vs. iteration.
With the checkbox off, this demo only refines the points, holding the 5 camera poses fixed, exactly as before. Check it and every camera's pose is also perturbed away from truth at scene start (simulating imperfect initial calibration) and refined every step — jointly with the points, using the block-coordinate version of the same Gauss-Newton loop above (points refined with poses held fixed, then poses refined with points held fixed, repeated). Camera 0 is always left untouched — that pins the gauge freedom explained below.
Jacobian sparsity: one row per observation, columns are the 5×6 camera-pose block ■ followed by the 8×3 point block ■. Each row only touches one camera-column-group and one point-column-group.
Ucc (30×30)
W (30×24)
Vpp (24×24)
|JᵀJ| for this exact toy problem, split into camera-camera, camera-point, and point-point blocks. Only the block-diagonal structure of Ucc and Vpp is exploited by the Schur complement below — W is left dense.
The Jacobian is almost all zeros
Why naive bundle adjustment doesn't scale
Stack every camera's 6 pose unknowns and every point's 3 position unknowns into one long vector, and stack every residual (one observation, camera i seeing point j) into one long vector r. Gauss-Newton needs the Jacobian J = ∂r/∂(cameras, points) of that whole residual vector. Here's the fact that makes bundle adjustment special: the residual for observation (i,j) is r⁅₠ = π(camera⁅, point₠) − x⁅₠ — it is a function of exactly two blocks of unknowns, camera i's 6 params and point j's 3 params, and literally nothing else. Its derivative with respect to every other camera's params and every other point's params is exactly, structurally zero — not small, zero.
So each row-pair of J (one observation, 2 rows for u,v) has nonzero entries in only 9 of the (6·#cameras + 3·#points) columns. Grouping columns as "all camera blocks" then "all point blocks", J splits cleanly:
row (i,j): zeros … ↑ ∂r⁅₠/∂camera⁅ ↑ … zeros | zeros … ↑ ∂r⁅₠/∂point₠ ↑ … zeros
Rendered as a grid (rows = observations, columns = parameters), the camera half of J looks block-diagonal-ish — each row lights up exactly one 6-wide camera stripe — and the point half looks the same with 3-wide point stripes. The demo above renders this pattern for the toy 5-camera, 8-point cube directly: scroll up and toggle "optimize camera poses" to see it change shape live as visibility changes.
Eliminate the points first: the Schur complement
The one trick that makes BA practical
The Gauss-Newton normal equations are JᵀJ·Δ = −Jᵀr. Split the unknown update Δ into a camera part Δc and a point part Δp, and — because J = [Jc|Jp] — the normal equations split into a matching 2×2 block system:
where Ucc = JcᵀJc, Vpp = JpᵀJp, W = JcᵀJp, and ec, ep are the matching pieces of −Jᵀr. The sparsity from the last section says something very precise about these three blocks: an entry of Ucc connecting camera i to camera i′ (i≠i′) sums over residuals that depend on both cameras' params at once — but no residual ever does, since one observation only ever involves one camera. So Ucc is exactly block-diagonal (one 6×6 block per camera), and by the identical argument so is Vpp (one 3×3 block per point). Only W, whose (i,j) block sums over the single residual that does involve both camera i and point j, is nonzero off the block-diagonal — dense-ish, patterned by whichever cameras actually see whichever points.
A block-diagonal matrix inverts block-by-block, for free (invert each tiny 3×3 point block independently — trivial compared to inverting the whole thing). That's the opening Vpp gives you. Take the bottom block row of the system, WᵀΔc + VppΔp = ep, and solve it for Δp:
Substitute that into the top block row, UccΔc + WΔp = ec, and the Δp terms rearrange into a system purely in Δc — the Schur complement, a.k.a. the reduced camera system:
How dense the reduced system actually is depends on the visibility graph: the bipartite graph with cameras on one side, points on the other, and an edge for every observation. W's nonzero pattern is this graph, and W Vpp⁻¹Wᵀ fills in a camera-camera entry exactly when two cameras share at least one point — so a well-connected scene (lots of shared points between many camera pairs) gives a denser but still far smaller reduced system, while a weakly-connected scene (long chains of cameras with thin overlap, or loops that barely close) gives a sparser reduced system that's harder to solve well — the same weak connectivity that lets small errors accumulate and drift, the failure mode that motivates incremental SfM's loop-closure problem later in this series.
Robust losses: the outlier problem lives inside the optimizer too
IRLS: same machinery, one extra step
RANSAC filters out most bad correspondences before bundle adjustment ever runs — but "most" isn't "all." A handful of bad matches always survive: a repeated texture, a moving object, a mismatch just inside the inlier threshold. Plain squared loss ρ(r) = r² gives every residual influence proportional to 2r — so one leftover outlier residual, ten times larger than the rest, pulls on the solution ten times harder than any good one, and can visibly drag a converged reconstruction off course. This is the exact same failure mode that motivated RANSAC in the first place, just relocated: now it's happening inside the optimizer, on data RANSAC already mostly cleaned.
The fix is the same idea as RANSAC's inlier threshold, made smooth and differentiable: a robust loss ρ(r) that grows like r² for small residuals (so it still behaves like ordinary least squares near the optimum) but grows much more slowly — or not at all — for large ones.
— quadratic near zero, linear beyond δ: caps an outlier's influence at a constant rather than letting it grow with r.
— grows only logarithmically for large r: even more aggressive down-weighting than Huber, but non-convex, so it can (rarely) trap Gauss-Newton in a bad local minimum.
— fully saturates: past the cutoff c the loss is flat, its gradient is exactly zero, and the residual is hard-rejected — it stops influencing the solve at all.
None of these need a new solver. Each robust loss can be implemented as plain weighted least squares, where the weight is recomputed from the current residual at every iteration — Iteratively Reweighted Least Squares (IRLS):
normal equations become: JᵀWJ·Δ = −JᵀWr, W = diag(w(r₁), w(r₂), …)
At the current estimate, compute every residual, turn it into a weight w(r) (1 for a small residual, shrinking toward 0 for a large one), solve the ordinary weighted Gauss-Newton step — the exact same JᵀJ, same block-sparse structure, same Schur complement from the last section, just with each observation's contribution to Ucc, Vpp, W, ec, ep scaled by its current weight — then repeat. Robust bundle adjustment isn't a different algorithm; it's this exact same machinery with one extra reweighting pass bolted onto the top of every iteration.
One 3D point is observed by 4 fixed cameras. Camera 0's observation can be dragged arbitrarily far from its true pixel (simulating one bad correspondence that survived RANSAC); the other three stay clean. Each checked loss is solved to convergence via IRLS from the same initial triangulation, and its result is plotted against ground truth.
Black = ground truth, diamonds = cameras (camera 0 is the one with the draggable outlier).
Gauge freedom: the whole reconstruction can float
A hidden rank deficiency
Reprojection error depends only on the geometry between cameras and points — how each camera's ray relates to each point, not where the whole configuration happens to sit in an abstract world frame. Concretely: take every camera pose and every 3D point and apply the exact same rigid transform (a rotation and translation — or, since scale is separately unobservable from monocular two-view geometry too, as Part 10's four-candidate-pose decomposition of E only fixes t up to scale, a full similarity transform including a uniform scale factor) to all of them at once. Every camera-to-point ray is exactly preserved — same direction, same relative depth ratios — so every reprojected pixel is bit-for-bit unchanged and the cost C doesn't move by a single floating-point unit.
A direction that changes every unknown but leaves the cost function exactly flat is, by definition, a direction of zero curvature — a null vector of the Hessian JᵀJ. So JᵀJ is exactly rank-deficient by the dimension of that invariance group: 6 if only a rigid transform is free (3 rotation + 3 translation), 7 if scale is unconstrained too (the monocular case, by default — see Part 10). A rank-deficient Hessian can't be inverted; naive Gauss-Newton on an un-anchored bundle-adjustment problem is solving a system with infinitely many equally-good solutions differing only by where the whole reconstruction sits, and the linear solve either fails outright or wanders arbitrarily along that flat direction.
Drag the sliders — with the datum free, the cube and cameras visibly move but the cost stays ≈0. Pin camera 0 and moving them stops being free.
Updating rotations without breaking them
SO(3), SE(3), and why the demo above froze poses in the first place
A rotation matrix R has 9 numbers but only 3 real degrees of freedom — the other 6 are eaten by the orthogonality constraint RᵀR = I. An ordinary Gauss-Newton step treats every unknown as free: R + ΔR for some 3×3 ΔR straight out of the normal equations. Do that even once and R + ΔR is, generically, no longer orthogonal — not a rotation at all anymore, and every camera that "drifts off SO(3)" this way corrupts the whole reconstruction it feeds into.
The fix reuses the skew-symmetric matrix [·]× from the essential-matrix derivation (E = [t]×R), where [v]× is the matrix with [v]×w = v×w. Instead of updating R additively, update it multiplicatively by a small rotation built from a free 3-vector δ via the matrix exponential:
Because exp([δ]×) is a genuine rotation matrix for any δ (it's the Rodrigues formula — the same exponential map that turns an axis-angle vector into a rotation), Rnew is automatically, exactly orthogonal again — no re-normalization hack required. Gauss-Newton's actual unknown for the rotation part is the free 3-vector δ, not the 9 entries of R itself; the Jacobian column for "camera i's rotation" is ∂r/∂δ, evaluated at δ=0. Paired with 3 more unknowns for translation, that's exactly the 6 pose DOF per camera claimed all the way back at the start of this section.
This is precisely the missing piece the toy demo skipped by freezing camera poses: naively perturbing R's 9 numbers and adding the result back is exactly the trap described above. The "optimize camera poses" toggle in the demo implements the rule on this page — every pose update is R ← R·exp([δ]×), computed by finite-differencing the reprojection residual with respect to δ at δ=0, never by touching R's entries directly.
Orthogonality error per applied increment (log scale feel: the additive curve is visually off the chart within ~10 steps). Green = exp map (Rodrigues), blue = additive + polar repair, red = raw additive.
Series recap
The whole pipeline, end to end
| Part | Question answered |
|---|---|
| 1 · Pinhole camera | How does one camera turn a 3D point into a pixel? |
| 2 · Epipolar geometry | Given a match in one image, where can it be in another? |
| 3 · Homography | When can one image predict another exactly? |
| 4 · Triangulation & pose | Given two views, where exactly is the 3D point? |
| 5 · Trifocal tensor | What does a third view add? |
| 6 · Bundle adjustment | How do all cameras and all points agree at once? |
| Schur complement | (Ucc − WVpp⁻¹Wᵀ)Δc = ec − WVpp⁻¹ep — eliminate the (many, cheap) points, solve the much smaller (few, dense) camera-only system, back-substitute |
| Robust losses | Huber: quadratic/linear split at δ; Cauchy: (c²/2)log(1+r²/c²); Tukey: saturates to zero gradient past c |
| IRLS weight | w(r) = ρ′(r)/r, recomputed from the current residual every iteration — robust BA is ordinary BA plus this one reweighting step |
| Gauge freedom | rigid transform of the whole reconstruction is invisible to cost ⇒ JᵀJ is rank-deficient by 6 (rigid) or 7 (+scale); fix by pinning camera 0 (and one baseline length) |
| SE(3) pose update | Rnew = Rold·exp([δ]×) — update via the exponential map, never additively on R's 9 entries |
That's the field, start to finish: perspective projection sets up the problem, epipolar geometry and homographies constrain matches between views, triangulation and the trifocal tensor turn constraints into 3D structure, and bundle adjustment ties the whole thing together into one globally consistent reconstruction — the same nonlinear least-squares machinery from the Nonlinear Optimization series, aimed at cameras and 3D points instead of a single robot.