Bundle adjustment

Bundle adjustment is the estimator at the end of every geometric vision pipeline: refine all camera poses and all 3-D points simultaneously so that every observation is explained as well as possible. Stated plainly it is a nonlinear least-squares problem with hundreds of thousands of unknowns, which should be hopeless. It is not, because the problem has a very particular sparsity — no residual couples two cameras or two points — and exploiting that structure is what this page is really about.

Symbols used on this page are defined in the Nomenclature → optimisation & robust fitting.

The reconstruction that BA refines is built on the SfM page; the pose-only special case is on the pose-graph page.

1 · The problem

1.1 The objective

With $m$ cameras of parameters $\mathbf c_i$ and $n$ points $\mathbf X_j$, and a visibility indicator $v_{ij} \in \{0,1\}$ recording whether camera $i$ observes point $j$:

$$ \min_{\{\mathbf c_i\},\,\{\mathbf X_j\}} \sum_{i=1}^{m}\sum_{j=1}^{n} v_{ij}\, \rho\!\left(\left\| \mathbf x_{ij} - \pi\!\left(\mathbf c_i, \mathbf X_j\right) \right\|^2_{\boldsymbol\Sigma_{ij}}\right) . $$

With $\rho$ the identity and Gaussian image noise this is the maximum-likelihood estimate of the whole reconstruction. The Mahalanobis norm weights each observation by its own uncertainty, which matters because features at different scales have genuinely different localisation accuracy.

1.2 How big it gets

A dense Cholesky factorisation costs $O(N^3/3)$. For the building case that is $\sim 10^{18}$ flops — years of compute. The same problem solved with the structure exploited takes seconds. The entire practical subject is that gap.

2 · Nonlinear least squares

2.1 Gauss-Newton and Levenberg–Marquardt

Linearise the residual vector $\mathbf r(\mathbf p)$ about the current estimate, $\mathbf r(\mathbf p + \delta) \approx \mathbf r + \mathbf J\delta$, and minimise the resulting quadratic:

$$ \left(\mathbf J^\top\mathbf J\right)\delta = -\mathbf J^\top\mathbf r \qquad\text{(Gauss-Newton)} $$

Gauss-Newton converges quadratically near the solution but diverges when the linearisation is poor. Levenberg–Marquardt adds damping to interpolate toward gradient descent:

$$ \left(\mathbf J^\top\mathbf J + \lambda\,\mathbf D\right)\delta = -\mathbf J^\top\mathbf r , \qquad \mathbf D = \mathbf I \ \text{or}\ \operatorname{diag}\!\left(\mathbf J^\top\mathbf J\right). $$

Large $\lambda$ gives small, safe, gradient-like steps; small $\lambda$ recovers fast Gauss-Newton convergence. The standard adaptive rule increases $\lambda$ when a step fails to reduce the cost and decreases it when it succeeds. Using $\mathbf D = \operatorname{diag}(\mathbf J^\top\mathbf J)$ rather than $\mathbf I$ makes the damping scale-invariant, which matters here because rotation, translation and point coordinates have wildly different units and magnitudes.

Crucially, damping also guarantees invertibility. BA's normal-equation matrix is often singular or near-singular — from gauge freedom (§5), poorly triangulated points, or cameras with too few observations — and the $\lambda\mathbf D$ term is what keeps the solve well posed. LM is not merely a convergence safeguard here; it is load-bearing.

2.2 Rotations on the manifold

Rotations do not live in a vector space, so the update must respect $SO(3)$. Parameterise the increment in the Lie algebra and apply it multiplicatively:

$$ \mathbf R \leftarrow \exp\!\left(\left[\delta\boldsymbol\omega\right]_\times\right)\mathbf R . $$

This keeps the rotation exactly orthonormal at every iteration, has no singularities (unlike Euler angles), and needs no normalisation constraint (unlike a 4-parameter quaternion in the state). The attitude page derives the machinery; the same multiplicative-increment discipline appears in the error-state filters on the IMU page.

3 · The sparsity structure

3.1 Where it comes from

The Jacobian of one observation $\mathbf x_{ij}$ is nonzero in exactly two places — the block for camera $i$ and the block for point $j$:

$$ \frac{\partial \mathbf r_{ij}}{\partial \mathbf c_k} = \mathbf 0 \ \ \forall k \ne i, \qquad \frac{\partial \mathbf r_{ij}}{\partial \mathbf X_l} = \mathbf 0 \ \ \forall l \ne j . $$

This is the whole reason BA is tractable. Ordering the unknowns as cameras first, then points, gives the famous block structure

$$ \mathbf J^\top\mathbf J = \begin{bmatrix} \mathbf U & \mathbf W \\ \mathbf W^\top & \mathbf V \end{bmatrix}, \qquad \begin{bmatrix} \mathbf U & \mathbf W \\ \mathbf W^\top & \mathbf V \end{bmatrix} \begin{bmatrix} \delta_{\mathbf c} \\ \delta_{\mathbf X} \end{bmatrix} = \begin{bmatrix} \mathbf g_{\mathbf c} \\ \mathbf g_{\mathbf X} \end{bmatrix} . $$

with three properties that get exploited:

• $\mathbf U$ is block diagonal with $6\times6$ blocks — no camera couples to another camera directly.
• $\mathbf V$ is block diagonal with $3\times3$ blocks — no point couples to another point directly. This is the important one.
• $\mathbf W$ is sparse, with one $6\times3$ block per actual observation.

3.2 The Schur complement

Because $\mathbf V$ is block diagonal, $\mathbf V^{-1}$ is just $n$ independent $3\times3$ inversions — trivial and parallel. So eliminate the points. From the second block row, $\delta_{\mathbf X} = \mathbf V^{-1}\left(\mathbf g_{\mathbf X} - \mathbf W^\top\delta_{\mathbf c}\right)$, and substituting into the first gives the reduced camera system:

$$ \boxed{\;\underbrace{\left(\mathbf U - \mathbf W\mathbf V^{-1}\mathbf W^\top\right)}_{\mathbf S} \delta_{\mathbf c} = \mathbf g_{\mathbf c} - \mathbf W \mathbf V^{-1}\mathbf g_{\mathbf X}\;} $$

Solve this $6m \times 6m$ system for the cameras, then back-substitute for the points. The size now depends only on the number of cameras, and in a typical problem $n \gg m$ by two or three orders of magnitude.

The Schur complement is the single most important algorithmic idea in bundle adjustment, and it is exactly the same operation as marginalising the points out of a Gaussian — see the pose-graph page, where marginalisation is used deliberately to keep a sliding window bounded.

3.3 Fill-in: the catch

$\mathbf S$ is not as sparse as $\mathbf U$. The term $\mathbf W\mathbf V^{-1}\mathbf W^\top$ creates an entry between cameras $i$ and $k$ whenever they share a point. So the reduced system's sparsity is the camera co-visibility graph. Two regimes follow, and they call for different solvers:

• Sequential capture (video, driving, a SLAM trajectory): a camera shares points only with its temporal neighbours, so $\mathbf S$ is banded and a sparse Cholesky is near-linear. This is the happy case.
• Wide co-visibility (photos of one landmark from everywhere, a small object scanned all round): most cameras see most points, $\mathbf S$ is nearly dense, and direct factorisation returns to $O(m^3)$.

The simulation makes this visible: the fill-in pattern in the reduced system is drawn explicitly, and driving the visibility fraction up turns $\mathbf S$ dense while driving it down keeps it sparse.

4 · Solving the reduced system

A useful practical note on the iterative route: an inexact step is usually fine. Inside LM there is no need to solve the linear system to machine precision, because the linearisation itself is approximate — a loose CG tolerance, tightened as the outer iteration converges, saves a great deal of work. The standard preconditioners are block-Jacobi on $\mathbf S$ (cheap) and the "visibility-based" preconditioners that approximate the co-visibility graph with a cheaper structure.

5 · Gauge freedom

5.1 What is unobservable

Bundle adjustment from images alone cannot determine the absolute placement of the reconstruction. Rotating, translating and (monocular) scaling the entire scene while moving all cameras identically leaves every projection unchanged. So the cost function has a flat valley, and $\mathbf J^\top\mathbf J$ is rank deficient by:

5.2 Fixing it

A singular system must be regularised or the solve fails — or worse, drifts along the null space accumulating meaningless numerical motion. The options:

• Fix the first camera entirely and one distance (monocular). Simple and standard. The drawback is that it makes the reported covariance conditional on that arbitrary choice.
• Trivial gauge / anchoring. Add a weak prior pulling the solution toward the current estimate. Cheap, and keeps the system positive definite.
• Free-network adjustment. Constrain the sum of corrections to zero, distributing the gauge symmetrically rather than pinning one camera. This is the photogrammetry tradition and gives more meaningful covariances.
• External measurements. GNSS positions, a known baseline, a scale bar. These remove the gauge physically rather than by convention — the right answer when available.
• Rely on LM damping. The $\lambda\mathbf D$ term makes the system invertible, so BA "works" without any explicit gauge fix. This is what most implementations quietly do, and it is acceptable — but the resulting covariance is then meaningless in the gauge directions.

6 · Robustness

6.1 Robust kernels

Even after RANSAC, surviving outliers — a mismatch that happened to be geometrically consistent, a badly triangulated point, a moving object — will drag a least-squares solution. The fix is a robust loss $\rho$, most commonly Huber:

$$ \rho_\delta(e) = \begin{cases} \tfrac12 e^2, & |e| \le \delta \\ \delta\left(|e| - \tfrac12\delta\right), & |e| \gt \delta \end{cases} $$

Quadratic near zero so that good observations behave normally, linear beyond $\delta$ so the influence of a bad one is bounded. It is implemented as iteratively reweighted least squares: compute a weight $w_i = \rho'(e_i)/e_i$ per residual, solve the weighted normal equations, recompute. This changes nothing structurally — the sparsity and the Schur complement are untouched — so robustness is essentially free.

Set $\delta$ from the expected noise, typically at the $\chi^2$ quantile for a couple of pixels, exactly as in RANSAC §4.1. Cauchy and Tukey losses reject harder (their influence decays to zero) but are non-convex, so they can create local minima and should be introduced only after a decent solution exists — a graduated schedule that starts near-quadratic and sharpens is the robust practice.

6.2 Outlier removal between iterations

Robust kernels down-weight; they do not remove. It is standard to alternate: run BA to convergence, discard observations whose residual exceeds a threshold, re-run. Also worth removing are points with too few observations (under-constrained) and points with insufficient parallax (see SfM §3.2), both of which add unknowns without adding real information and can make $\mathbf V$ ill-conditioned.

7 · Variants in practice

8 · Practical checklist

• Exploit the sparsity. A dense solver is not a simpler version of BA, it is an unusable one. Use an established library — Ceres, g2o, GTSAM — rather than writing the linear algebra.
• Use LM, not plain Gauss-Newton, with scale-invariant damping.
• Update rotations multiplicatively on the manifold.
• Fix the gauge explicitly, and know which convention you used before quoting any covariance.
• Use a robust kernel — Huber by default — with $\delta$ derived from the pixel noise.
• Prune before solving: drop low-parallax points, points with too few views, and observations with huge residuals.
• Weight observations by feature-scale uncertainty; a coarse-scale detection is not as accurate as a fine one.
• Initialise well. BA is a refiner, not a solver — it converges to a local minimum near where it started, so a bad initialisation gives a confidently wrong reconstruction.
• Choose the ordering for sparse factorisation (AMD/nested dissection); the wrong ordering can cost orders of magnitude.
• Monitor the cost per observation, not the total. Total cost grows with problem size and tells you nothing about fit quality.

9 · Further reading

• Triggs, McLauchlan, Hartley & Fitzgibbon (1999), "Bundle Adjustment — A Modern Synthesis," Vision Algorithms workshop. The definitive reference. Long, dense, and worth the effort — especially the sections on gauge freedom and covariance interpretation.
• Hartley & Zisserman, Multiple View Geometry, 2nd ed., appendix 6. The sparse solution and the Schur complement worked through concretely.
• Agarwal, Snavely, Seitz & Szeliski (2010), "Bundle Adjustment in the Large," ECCV. The inexact-Schur and preconditioned-CG approach that made city-scale BA feasible.
• Lourakis & Argyros (2009), "SBA: A Software Package for Generic Sparse Bundle Adjustment," ACM TOMS. A clear description of a real implementation.
• Kümmerle et al. (2011), "g2o: A General Framework for Graph Optimization," ICRA.
• Agarwal, Mierle et al., Ceres Solver documentation. The best practical guide to writing a real BA problem, including ordering, preconditioners and loss functions.
• Dellaert & Kaess (2006), "Square Root SAM," IJRR. The factor-graph/QR view of the same problem, and the bridge to GraphSLAM.
• Forster et al. (2017), "On-Manifold Preintegration for Real-Time Visual-Inertial Odometry," IEEE T-RO. How IMU factors enter a BA problem properly.

Related pages: Structure from motion · PnP · GraphSLAM & pose graphs · Attitude & rotations · Extended Kalman filter