Bundle Adjustment

Bundle adjustment refines a whole reconstruction at once: every camera’s pose, often its intrinsics, and every 3D point are adjusted together to minimize the total reprojection error over all observations. It is the final step of most structure-from-motion pipelines and runs repeatedly inside visual SLAM systems, and points missing from some images need no special treatment.

The name, as Triggs, McLauchlan, Hartley and Fitzgibbon (2000) explain it, refers to the bundles of light rays that leave each 3D point and converge on each camera center, adjusted jointly with respect to both point and camera positions. The method comes from photogrammetry: their historical appendix credits Duane C. Brown and co-workers, working for the U.S. Air Force in 1957–59, with the basic bundle method for aerial mapping with calibrated cameras. Brown’s 1958 method already eliminated the point parameters to leave a small system in the cameras, the reduction derived below.

Problem

Given $m$ cameras, $n$ 3D points, and the measured image positions $\mathbf{x}_{ij}$ of point $j$ in each image $i$ that sees it, find camera parameters $\mathbf{a}_i$ (pose, and possibly focal length and distortion) and points $\mathbf{X}_j$ that minimize

$$E(\mathbf{a}, \mathbf{X}) = \sum_{(i,j) \in \mathcal{V}} \rho\Big(\big\lVert \mathbf{x}_{ij} - \pi(\mathbf{a}_i, \mathbf{X}_j) \big\rVert^2\Big),$$

where $\pi$ projects a point through a camera (see camera parameters), $\mathcal{V}$ is the set of observations that exist, and $\rho$ is the identity for plain least squares or a robust loss. It is a nonlinear least-squares problem, set apart by its size and structure.

Inputs and Outputs

Inputs: the observations $\mathbf{x}_{ij}$ with their visibility pattern (feature tracks), an initial estimate of every camera and point, and a camera model stating which parameters are free and which are shared, for example one intrinsic matrix for all images from the same camera. Optional: per-observation covariances, a robust loss, and fixed control points or poses.

Outputs: refined cameras and points, the final residuals, and optionally parameter covariances.

Intuition

In a reconstruction built step by step, errors compound: a slightly wrong camera gives slightly wrong triangulated points, which give a slightly wrong next camera. Bundle adjustment lets every quantity move at once until no small change to any camera or point reduces the image error.

The coupling is sparse: each observation involves one camera and one point, and points interact only through the cameras that see them. Given the cameras, each point could be solved on its own. The efficient solution therefore eliminates the points algebraically, solves a small system for the cameras, and then recovers each point independently.

Algorithm

Bundle adjustment is Levenberg–Marquardt, or a trust-region variant, with a structured linear solve in each iteration:

1. Initialize cameras and points, typically from incremental structure from motion: two-view geometry from the essential matrix, new cameras by perspective-n-point, new points by triangulation.
2. Fix the gauge, for example by holding one camera fixed and one distance constant.
3. Linearize: compute residuals $\mathbf{e}_{ij}$ and the Jacobian blocks $A_{ij} = \partial \mathbf{e}_{ij} / \partial \mathbf{a}_i$ and $B_{ij} = \partial \mathbf{e}_{ij} / \partial \mathbf{X}_j$.
4. Reduce and solve: form the damped normal equations, eliminate the points with the Schur complement, solve the reduced camera system, and back-substitute for the points.
5. Update the parameters (rotations through a local parametrization), adapt the damping, and repeat until the cost or the step stops changing.

Mathematical Formulation

Parametrization

Rotations are not a vector space, so each update is a small rotation vector $\boldsymbol{\omega}$ applied to the current estimate, $R_i \leftarrow R_i \exp([\boldsymbol{\omega}]_\times)$, or a quaternion that is renormalized after each step; Triggs et al. warn that global three-parameter forms such as Euler angles cause numerical trouble near their singularities. A calibrated camera then has $c = 6$ parameters. The Bundle Adjustment in the Large problems of Agarwal et al. (2010) use $c = 9$: pose, one focal length, and two radial distortion coefficients.

Points are usually Euclidean, which behaves badly for distant points: the cost becomes nearly flat in depth, and infinity lies at infinite parameter values. Triggs et al. recommend homogeneous coordinates $(X, Y, Z, W)^\top$ with unit norm for such points, which reach infinity at $W = 0$. In filter-based monocular SLAM, Civera, Davison and Montiel (2008) instead parametrized each point by its inverse depth along the ray from the camera that first observed it, which likewise stays finite at infinity.

Sparse structure

Order the unknowns as all camera parameters followed by all point parameters. Each observation contributes two rows to the Jacobian $J$ with exactly two nonzero blocks, $A_{ij}$ ($2 \times c$) and $B_{ij}$ ($2 \times 3$). The damped normal equations $(J^\top J + \lambda D)\,\boldsymbol{\delta} = -J^\top \mathbf{e}$ then have the block form

$$\begin{bmatrix} U & W \\ W^\top & V \end{bmatrix} \begin{bmatrix} \boldsymbol{\delta}_a \\ \boldsymbol{\delta}_b \end{bmatrix} = \begin{bmatrix} \boldsymbol{\epsilon}_a \\ \boldsymbol{\epsilon}_b \end{bmatrix},$$

with $U = \mathrm{diag}(U_1, \dots, U_m)$, $U_i = \sum_j A_{ij}^\top A_{ij}$, and $V = \mathrm{diag}(V_1, \dots, V_n)$, $V_j = \sum_i B_{ij}^\top B_{ij}$, each including its share of the damping. The off-diagonal block $W$ has a $c \times 3$ block $W_{ij} = A_{ij}^\top B_{ij}$ for each observation, and the right-hand sides are $\boldsymbol{\epsilon}_{a,i} = -\sum_j A_{ij}^\top \mathbf{e}_{ij}$ and $\boldsymbol{\epsilon}_{b,j} = -\sum_i B_{ij}^\top \mathbf{e}_{ij}$. The point block $V$ is block diagonal because no observation involves two points (Hartley and Zisserman, 2004, Appendix 6; Triggs et al., 2000, §5). $U$ is block diagonal only while no parameters are shared; shared intrinsics add off-diagonal blocks to it.

The Schur complement and the reduced camera system

The bottom block row gives $\boldsymbol{\delta}_b = V^{-1}(\boldsymbol{\epsilon}_b - W^\top \boldsymbol{\delta}_a)$. Substituting into the top row eliminates the points:

$$\underbrace{\left(U - W V^{-1} W^\top\right)}_{S} \boldsymbol{\delta}_a = \boldsymbol{\epsilon}_a - W V^{-1} \boldsymbol{\epsilon}_b .$$

$S$ is the Schur complement of $V$, the reduced camera system. Since $V$ is block diagonal, $V^{-1}$ needs only $n$ inversions of $3 \times 3$ matrices, and $S$ accumulates one contribution per point:

$$S_{ik} = \delta_{ik} U_i - \sum_{j} W_{ij} V_j^{-1} W_{kj}^\top ,$$

where the sum runs over points seen by both cameras $i$ and $k$. After solving for $\boldsymbol{\delta}_a$, each point is recovered independently by back-substitution:

$$\boldsymbol{\delta}_{b,j} = V_j^{-1}\Big(\boldsymbol{\epsilon}_{b,j} - \sum_i W_{ij}^\top \boldsymbol{\delta}_{a,i}\Big).$$

The block $S_{ik}$ is nonzero only if cameras $i$ and $k$ share a point (Agarwal et al., 2010). For a long video sequence $S$ is close to banded; for a photo collection of one landmark, many images overlap and $S$ fills in. Triggs et al. note that when camera parameters outnumber point parameters, as in some long image sequences, it is cheaper to eliminate the cameras instead.

Gauge freedom

The images determine a reconstruction only up to a transformation that leaves all projections unchanged. For a calibrated reconstruction this is a 3D similarity transformation, with $3 + 3 + 1 = 7$ degrees of freedom (translation, rotation, scale); for an uncalibrated, projective one it is a 3D projective transformation with 15 (Triggs et al., 2000, §9). Along these directions the cost is constant, so $J^\top J$ has a null space of dimension at least 7 or 15.

The singularity is fixed by choosing a gauge. The simplest choice, which Triggs et al. suggest for computation, is a trivial gauge: freeze enough parameters, for example the first camera’s pose and one coordinate of the second camera’s position (fixing scale). Alternatives constrain the reconstruction’s centroid, orientation, and scale, or leave the gauge free: damping still makes each linear system solvable, but the reconstruction can drift along the free directions. The choice does not change the geometry or the residuals, but it changes the parameter covariances: in a camera-centered gauge the points can look uncertain even when their shape is well determined.

Parameters

• Camera model: which intrinsics and distortion coefficients are free or shared. Freeing poorly constrained parameters lowers the cost without improving accuracy.
• Gauge: which parameters are held fixed, or which constraints are imposed.
• Robust loss and its scale: a Huber or Cauchy loss (see nonlinear least squares) limits the influence of mismatches that survived RANSAC; the problem remains a reweighted nonlinear least-squares problem (Agarwal et al., 2010).
• Damping and stopping: initial damping, an iteration limit, and tolerances on the cost, step, and gradient.
• Linear solver: dense or sparse factorization of $S$, or an iterative solver (see Complexity).

Complexity

Let $p = mc$ be the number of camera parameters and $q = 3n$ the number of point parameters. Factoring the full normal equations densely costs $O((p + q)^3)$, and since points usually far outnumber cameras, that is dominated by the points. With the Schur complement, inverting $V$ and back-substitution are linear in $n$; forming $S$ costs one $c \times c$ product for each point and each pair of cameras that observe it; and factoring $S$ densely costs $O(p^3)$ time and $O(p^2)$ memory. Hartley and Zisserman (2004, §A6.3) summarize the gain: with the structure exploited, an iteration’s cost grows linearly in the number of points, against cubically without it.

Agarwal et al. (2010) report that dense factorization of $S$ is practical up to a few hundred cameras. Beyond that, sparse Cholesky factorization with a fill-reducing ordering exploits the sparsity of $S$. For thousands of images they studied truncated Newton steps: preconditioned conjugate gradients solves the reduced system approximately, multiplying vectors by $S$ through $U$, $W$, and $V$ without forming it, with simple block diagonal preconditioners. In their experiments these solvers were up to an order of magnitude faster than factorization on problems of a few thousand images, used three to five times less memory, and solved the large Venice problem, on which the factorization solvers ran out of memory.

Implementation

The example builds a synthetic calibrated problem (5 cameras, 100 points each seen by 3 to 5 cameras, 1-pixel noise), freezes camera 0 and one translation entry of camera 1 as the gauge, and solves it with scipy.optimize.least_squares given the Jacobian’s sparsity pattern (jac_sparsity). It then checks the gauge freedom through the singular values of the unconstrained Jacobian and computes one Gauss–Newton step both directly and through the reduced camera system.

import numpy as np
from scipy.optimize import least_squares
from scipy.sparse import lil_matrix
from scipy.spatial.transform import Rotation

rng = np.random.default_rng(0)
f, sigma = 800.0, 1.0                      # focal length and pixel noise (per coordinate)
n_cams, n_pts = 5, 100

# Ground truth: points in a 4 x 4 x 2 m box, cameras on an arc about 8 m away, looking at it.
X_true = rng.uniform(-1, 1, (n_pts, 3)) * [2, 2, 1]
cams_true = []
for a in np.linspace(-0.4, 0.4, n_cams):
    R = Rotation.from_euler("y", -a).as_matrix()
    c = 8 * np.array([np.sin(a), 0.0, -np.cos(a)])          # camera centre
    cams_true.append(np.r_[Rotation.from_matrix(R).as_rotvec(), -R @ c])
cams_true = np.array(cams_true)                              # (rotation vector, t) per camera

# Each point is seen by 3 to 5 cameras.
obs = [(i, j) for j in range(n_pts)
       for i in rng.choice(n_cams, rng.integers(3, n_cams + 1), replace=False)]
cam_idx, pt_idx = np.array(obs).T

def project(cams, X):
    """Pixel coordinates (principal point at 0) of X[k] in camera cams[k]."""
    Xc = Rotation.from_rotvec(cams[:, :3]).apply(X) + cams[:, 3:]
    return f * Xc[:, :2] / Xc[:, 2:]

x_obs = project(cams_true[cam_idx], X_true[pt_idx])
x_obs += sigma * rng.standard_normal(x_obs.shape)

# Gauge: freeze camera 0 (6 parameters) and the first translation entry of camera 1 (scale).
free = np.ones(6 * n_cams + 3 * n_pts, bool)
free[:6] = False
free[9] = False

def unpack(p, full):
    full = full.copy()
    full[free] = p
    return full[:6 * n_cams].reshape(n_cams, 6), full[6 * n_cams:].reshape(n_pts, 3)

def residuals(p, full):
    cams, X = unpack(p, full)
    return (project(cams[cam_idx], X[pt_idx]) - x_obs).ravel()

def sparsity():
    """Residual k depends on 6 camera and 3 point parameters only."""
    A = lil_matrix((2 * len(obs), 6 * n_cams + 3 * n_pts), dtype=int)
    for k, (i, j) in enumerate(obs):
        A[2 * k:2 * k + 2, 6 * i:6 * i + 6] = 1
        A[2 * k:2 * k + 2, 6 * n_cams + 3 * j:6 * n_cams + 3 * j + 3] = 1
    return A.tocsc()[:, free]

# Initialization: perturbed cameras (about 2 degrees, 0.2 m) and points (0.2 m).
cams0 = cams_true + np.c_[rng.normal(0, 0.035, (n_cams, 3)), rng.normal(0, 0.2, (n_cams, 3))]
cams0[0], cams0[1, 3] = cams_true[0], cams_true[1, 3]       # gauge values are kept
X0 = X_true + rng.normal(0, 0.2, X_true.shape)
full0 = np.r_[cams0.ravel(), X0.ravel()]
p0 = full0[free]

def rms(p):
    return np.sqrt(np.mean(residuals(p, full0).reshape(-1, 2) ** 2) * 2)   # per point

A = sparsity()
N, d = len(obs), int(free.sum())
print(f"{n_cams} cameras, {n_pts} points, {N} observations, {d} free parameters")
print(f"Jacobian {A.shape[0]} x {A.shape[1]}, nonzero fraction {A.nnz / np.prod(A.shape):.4f}")
H = (A.T @ A).toarray() != 0
print(f"J^T J nonzero fraction {H.mean():.4f}")

res = least_squares(residuals, p0, jac_sparsity=A, x_scale="jac", args=(full0,))
cams, X = unpack(res.x, full0)
print(f"RMS before {rms(p0):8.3f} px, after {rms(res.x):.3f} px "
      f"(expected {sigma * np.sqrt(2 * (1 - d / (2 * N))):.3f}), {res.nfev} evaluations")
print(f"point error: before {np.sqrt(np.mean((X0 - X_true) ** 2)):.4f} m, "
      f"after {np.sqrt(np.mean((X - X_true) ** 2)):.4f} m (RMS per coordinate)")

# Gauge freedom: the Jacobian over ALL parameters (nothing frozen) is rank deficient.
full = np.r_[cams.ravel(), X.ravel()]
def r_all(q):
    c, P = q[:6 * n_cams].reshape(n_cams, 6), q[6 * n_cams:].reshape(n_pts, 3)
    return (project(c[cam_idx], P[pt_idx]) - x_obs).ravel()
def jac_all(q, h=1e-6):                                     # central differences
    return np.column_stack([(r_all(q + h * e) - r_all(q - h * e)) / (2 * h)
                            for e in np.eye(len(q))])
s = np.linalg.svd(jac_all(full), compute_uv=False)
print("smallest singular values / largest:", np.array2string(s[-9:] / s[0], precision=1))

# One Gauss-Newton step at the start, solved directly and by the Schur complement.
Jf = jac_all(full0)[:, free]
r = r_all(full0)
Hf, g = Jf.T @ Jf, -Jf.T @ r
nc = 6 * n_cams - 7                                         # free camera parameters
U, W, V = Hf[:nc, :nc], Hf[:nc, nc:], Hf[nc:, nc:]
ea, eb = g[:nc], g[nc:]
Vinv = np.zeros_like(V)
for k in range(n_pts):                                      # invert 3x3 blocks only
    b = slice(3 * k, 3 * k + 3)
    Vinv[b, b] = np.linalg.inv(V[b, b])
S = U - W @ Vinv @ W.T                                      # reduced camera system
da = np.linalg.solve(S, ea - W @ Vinv @ eb)
db = Vinv @ (eb - W.T @ da)                                 # back-substitution
d_full = np.linalg.solve(Hf, g)
print(f"reduced system {S.shape[0]} x {S.shape[0]} instead of {Hf.shape[0]} x {Hf.shape[0]}; "
      f"max difference {np.abs(np.r_[da, db] - d_full).max():.1e}")

Output:

5 cameras, 100 points, 405 observations, 323 free parameters
Jacobian 810 x 323, nonzero fraction 0.0234
J^T J nonzero fraction 0.1162
RMS before   50.052 px, after 1.107 px (expected 1.097), 6 evaluations
point error: before 0.2113 m, after 0.0112 m (RMS per coordinate)
smallest singular values / largest: [1.1e-02 8.9e-03 3.3e-11 2.0e-11 1.9e-11 1.6e-11 1.5e-11 1.3e-11 1.2e-11]
reduced system 23 x 23 instead of 323 x 323; max difference 1.7e-14

Only 2.3% of the Jacobian is nonzero, and the pattern lets SciPy estimate it by finite differences cheaply, perturbing columns that share no row together. The per-point RMS falls from 50 pixels to 1.107 pixels, close to the value $\sigma\sqrt{2}\sqrt{1 - d/2N} = 1.097$ expected at the optimum (see reprojection error); SciPy’s count of six evaluations excludes those used for the Jacobian. The Jacobian with nothing frozen has exactly seven singular values at the level of numerical noise, the similarity gauge, eight orders of magnitude below the rest. The reduced camera system has 23 unknowns instead of 323 and gives the same step to rounding error. In a real implementation, $S$ would be formed blockwise rather than from a dense $J^\top J$ as here, and least_squares does not use the Schur complement itself; its sparse trust-region solver (LSMR) works on the whole Jacobian.

Properties and Behavior

• Optimality. With independent Gaussian image noise, the minimum is the maximum-likelihood estimate of all cameras and points (Hartley and Zisserman, 2004, §18.1; see reprojection error). Robust losses trade some of this optimality for resistance to outliers.
• Local convergence. Like any Gauss–Newton method, bundle adjustment converges to a nearby local minimum; it cannot repair a wrong relative pose or many mismatched tracks. Incremental pipelines therefore adjust again as images are added (Hartley and Zisserman, 2004, §18.6).
• Second-order versus alternation. Alternating between resection (cameras from fixed points) and intersection (points from fixed cameras) minimizes the same cost, but converges slowly: Triggs et al. explain that corrections propagate through the network only about one step per iteration.

Limitations

• Initialization. Bundle adjustment is a refinement, not a reconstruction method; it needs a starting point near the solution (Hartley and Zisserman, 2004, §18.1).
• Weak geometry. Short baselines or long chains of images with little overlap leave directions that the data barely constrain; the cost is nearly flat along them, convergence slows, and a low final residual does not reveal the problem.

Variants

• Motion-only bundle adjustment fixes the points and refines cameras, which decouples into independent pose refinements, the nonlinear stage of perspective-n-point. The SBA package offers it as a separate routine, which its authors note is, strictly, not bundle adjustment.
• Structure-only bundle adjustment fixes the cameras and refines points, which decouples into independent triangulation problems, one $3 \times 3$ system per point.
• Local bundle adjustment optimizes a window of recent cameras and the points they see, holding older cameras fixed, which also fixes the gauge. Visual odometry and SLAM systems use it to bound the cost of each update. Agarwal et al. note that focusing on the recently modified part of a reconstruction works well for video and street-level sequences, whose overlap graphs are long and thin.
• Self-calibrating bundle adjustment also frees intrinsics and distortion; Triggs et al. date it to the 1960s in photogrammetry.

Among libraries, SBA (Lourakis and Argyros, 2009) is a C package applying sparse Levenberg–Marquardt with user-supplied camera and point models. Ceres Solver, an open-source C++ library used in production at Google, provides dense, sparse, and iterative Schur-based solvers, robust losses, and manifolds for rotations. g2o (Kümmerle et al., 2011) treats bundle adjustment and SLAM as optimization over a graph of variables and error terms.