Bundle Adjustment: Polishing the Reconstruction

"Every camera is a little wrong and every point is a little wrong, and at two in the morning I redistribute the blame across nineteen thousand parameters until nobody can point fingers."

A Sleep-Deprived Bundle Adjuster

The incremental loop of Section 14.2 makes locally reasonable decisions: each pose from PnP is the best fit to the map as it stood, each triangulation the best fit to the cameras as they stood. But "as it stood" is the problem: the map used to register camera 41 already contained the small errors of cameras 1 through 40, so errors do not merely persist, they compound. What the reconstruction needs is a stage where no estimate is privileged: cameras and points all become variables simultaneously, and the data, the observed pixel locations of the tracks from Section 14.1, arbitrates among them. That stage is bundle adjustment, named for the bundles of light rays connecting each 3D point to the cameras that observe it; the name predates computer vision by decades, coming from the photogrammetrists who adjusted aerial-survey bundles by hand calculation in the 1950s.

1. The Objective: Total Reprojection Error Intermediate

Let camera $i$ carry parameters $\boldsymbol{\theta}_i$ (rotation and translation; optionally intrinsics and distortion, as calibrated in Chapter 12), let $\mathbf{X}_j$ be the position of track $j$, and let $\mathcal{O}$ be the set of observations: $(i, j) \in \mathcal{O}$ when camera $i$ observes track $j$ at pixel $\mathbf{x}_{ij}$. With $\pi$ the projection function from Chapter 12, bundle adjustment solves

$$ \min_{\{\boldsymbol{\theta}_i\}, \{\mathbf{X}_j\}} \; \sum_{(i,j) \in \mathcal{O}} \rho\!\left( \left\lVert \pi(\boldsymbol{\theta}_i, \mathbf{X}_j) - \mathbf{x}_{ij} \right\rVert^2 \right), $$

where $\rho$ is a robust loss we return to in Section 5 (for now, read it as the identity). The objective is exactly the reprojection error that Chapter 13 used to evaluate triangulations, summed over the entire model. Two properties shape everything that follows. First, the problem is nonlinear: $\pi$ involves a rotation and a perspective division, so no closed-form solution exists and we must iterate from the initialization that incremental SfM supplies. Second, it is not fully determined: moving every camera and every point by a common rigid transform, or scaling the whole scene, changes nothing about reprojection. This seven-dimensional family of equivalent solutions (3 rotation, 3 translation, 1 scale) is the gauge freedom; solvers handle it by fixing one camera and one distance, or by letting the normal equations' damping absorb it.

2. The Solver: Gauss-Newton and Levenberg-Marquardt Intermediate

Stack all residuals $\mathbf{r}(\mathbf{p}) = \left( \pi(\boldsymbol{\theta}_i, \mathbf{X}_j) - \mathbf{x}_{ij} \right)_{(i,j)\in\mathcal{O}}$ into one vector, where $\mathbf{p}$ collects every parameter. Gauss-Newton linearizes the residuals around the current estimate via the Jacobian $J = \partial \mathbf{r} / \partial \mathbf{p}$ and solves the normal equations for an update $\boldsymbol{\delta}$:

$$ (J^\top J)\, \boldsymbol{\delta} = -J^\top \mathbf{r}. $$

Far from the solution the linearization can overshoot, so Levenberg-Marquardt (LM) adds adaptive damping,

$$ \left(J^\top J + \lambda \,\mathrm{diag}(J^\top J)\right) \boldsymbol{\delta} = -J^\top \mathbf{r}, $$

shrinking $\lambda$ when a step reduces the cost (behaving like fast Gauss-Newton) and growing it when a step fails (retreating toward cautious gradient descent). LM is the default optimizer of every serious bundle adjuster, and its character matters for intuition: unlike the first-order stochastic optimizers you will meet training networks in Chapter 18, LM uses curvature, converging in typically 10 to 50 iterations to sub-pixel consistency, and it is deterministic: the same input yields the same polished model.

3. The Sparsity Miracle and the Schur Complement Advanced

Count parameters for a modest reconstruction: 500 cameras at 6 parameters and 200,000 points at 3 parameters give $n = 603{,}000$ unknowns. The matrix $J^\top J$ is $n \times n$; storing it densely needs about 2.9 terabytes, and solving it densely costs $O(n^3)$. Bundle adjustment is feasible only because of the structure visible in the objective: each residual involves exactly one camera and exactly one point. The Jacobian row for observation $(i,j)$ is zero everywhere except in the 6 columns of camera $i$ and the 3 columns of point $j$, so $J^\top J$ assembles into the famous arrowhead form

$$ J^\top J = \begin{bmatrix} B & E \\ E^\top & C \end{bmatrix}, $$

where $B$ (cameras by cameras) is small, $C$ (points by points) is enormous but block diagonal with independent $3 \times 3$ blocks, and $E$ encodes which camera sees which point. Block-diagonal matrices invert block by block, which enables the Schur complement trick: eliminate all point parameters analytically and solve the small reduced camera system

$$ \left( B - E\, C^{-1} E^\top \right) \boldsymbol{\delta}_{\text{cam}} = -\mathbf{g}_{\text{cam}} + E\, C^{-1} \mathbf{g}_{\text{pt}}, $$

then back-substitute each point's update independently (and in parallel). The expensive linear solve shrinks from 603,000 unknowns to 3,000, a workload a laptop handles in milliseconds per iteration. Figure 14.3.1 makes the structure visual; once you have seen it, you cannot unsee why bundle adjustment scales.

4. Bundle Adjustment in Code Intermediate

SciPy's least_squares will not win speed records against Ceres, but it exposes every concept of this section honestly, and on problems up to a few hundred cameras it is perfectly serviceable. The implementation needs three ingredients: a residual function (project every observed point through its camera, subtract the observation), a parameter vector packing camera 6-vectors (axis-angle rotation plus translation) and point 3-vectors, and, crucially, the Jacobian sparsity pattern, without which SciPy would attempt dense finite differences over the full parameter count and never finish. The arrays cam_idx, pt_idx, and obs are exactly the observation list of the objective: row $k$ says camera cam_idx[k] saw point pt_idx[k] at pixel obs[k], data that falls straight out of the track table built in Section 14.1.

import numpy as np
from scipy.optimize import least_squares
from scipy.sparse import lil_matrix

def rotate(pts, rvecs):
    """Apply batched axis-angle rotations (Rodrigues formula)."""
    theta = np.linalg.norm(rvecs, axis=1, keepdims=True)
    with np.errstate(invalid="ignore"):
        k = np.nan_to_num(rvecs / theta)              # unit axes (0 where theta=0)
    cos, sin = np.cos(theta), np.sin(theta)
    return (cos * pts + sin * np.cross(k, pts)
            + (pts * k).sum(1, keepdims=True) * (1 - cos) * k)

def residuals(params, n_cams, n_pts, cam_idx, pt_idx, obs, K):
    cams = params[:n_cams * 6].reshape(n_cams, 6)     # [rvec | tvec] per camera
    pts = params[n_cams * 6:].reshape(n_pts, 3)
    Xc = rotate(pts[pt_idx], cams[cam_idx, :3]) + cams[cam_idx, 3:6]
    proj = Xc @ K.T
    proj = proj[:, :2] / proj[:, 2:]                  # perspective division
    return (proj - obs).ravel()

def jac_sparsity(n_cams, n_pts, cam_idx, pt_idx):
    """Mark which parameters each residual can possibly touch."""
    m = cam_idx.size * 2                              # two residuals per observation
    A = lil_matrix((m, n_cams * 6 + n_pts * 3), dtype=int)
    row = np.arange(cam_idx.size)
    for s in range(6):                                # this residual's camera...
        A[2 * row, cam_idx * 6 + s] = 1
        A[2 * row + 1, cam_idx * 6 + s] = 1
    for s in range(3):                                # ...and its point, nothing else
        A[2 * row, n_cams * 6 + pt_idx * 3 + s] = 1
        A[2 * row + 1, n_cams * 6 + pt_idx * 3 + s] = 1
    return A

x0 = np.hstack([cams_init.ravel(), pts_init.ravel()])
A = jac_sparsity(n_cams, n_pts, cam_idx, pt_idx)
res = least_squares(residuals, x0, jac_sparsity=A, method="trf",
                    loss="huber", f_scale=1.0, x_scale="jac",
                    args=(n_cams, n_pts, cam_idx, pt_idx, obs, K))

r0 = residuals(x0, n_cams, n_pts, cam_idx, pt_idx, obs, K)
print(f"reprojection RMSE: {np.sqrt((r0**2).mean()):.2f}px "
      f"-> {np.sqrt((res.fun**2).mean()):.2f}px in {res.nfev} evaluations")
# reprojection RMSE: 3.84px -> 0.47px in 23 evaluations

Watch what the printed before-and-after line says: the optimizer did not move any single camera dramatically, yet the collective consistency improved by nearly an order of magnitude. That is the signature of bundle adjustment: it harvests hundreds of small, correlated corrections that no per-camera or per-point estimator could see.

5. Robustness in Practice: Losses, Outliers, Scheduling Intermediate

A plain squared loss hands control of the model to its worst data: one contaminated track observation, 40 pixels off, contributes as much gradient as 1,600 observations that are one pixel off. Hence the robust loss $\rho$. The Huber loss is quadratic inside a threshold and linear beyond, so outliers retain influence but no longer dominate; the Cauchy loss goes further, asymptotically ignoring gross errors. Production systems pair a robust kernel with an eviction schedule: adjust, delete observations whose residuals remain large, adjust again, exactly the audit rhythm of Section 14.2 applied at global scale. Scheduling is the other practical lever. Running global bundle adjustment after every new camera would be ruinously slow, so incremental mappers interleave local BA (the new camera plus its covisible neighborhood, a few dozen cameras) after every registration with global BA at exponentially spaced milestones, when the model has grown by some percentage. And when many images share one physical camera, refining shared intrinsics inside the final global BA (self-calibration) recovers focal length and distortion better than any single image could; this is why COLMAP often reports better intrinsics than the EXIF data claimed.