import numpy as np
def project_point(R, t, K, X):
"""Project one 3D point into one calibrated camera."""
Y = R @ X + t
u = K[0, 0] * Y[0] / Y[2] + K[0, 2]
v = K[1, 1] * Y[1] / Y[2] + K[1, 2]
return np.array([u, v])
K = np.array([[800.0, 0.0, 320.0],
[0.0, 800.0, 240.0],
[0.0, 0.0, 1.0]])
R = np.eye(3)
t = np.array([0.2, -0.1, 0.0])
X = np.array([0.5, 0.3, 4.0])
x_true = project_point(R, t, K, X)
x_obs = x_true + np.array([0.7, -0.4])
residual = x_obs - project_point(R, t, K, X)
print("predicted pixel:", x_true)
print("observed pixel:", x_obs)
print("residual:", residual)
print("reprojection error:", np.linalg.norm(residual))18 Bundle Adjustment
Bundle adjustment is the optimization heart of geometric computer vision. After image matching, epipolar geometry, triangulation, and Structure from Motion have produced an initial reconstruction, bundle adjustment refines all camera parameters and all 3D points so that the predicted image measurements agree with the observed image measurements. This chapter develops bundle adjustment as a sparse nonlinear least-squares problem with a clear statistical interpretation. The emphasis is on the mathematics: reprojection residuals, Jacobians, Gauss–Newton, Levenberg–Marquardt, block sparsity, Schur complements, gauge freedom, robust losses, uncertainty, and modern differentiable versions used in AI systems.
18.1 Learning goals
After reading this chapter, students should be able to:
- formulate bundle adjustment as a nonlinear least-squares problem over camera and point variables;
- derive reprojection residuals for calibrated and uncalibrated camera models;
- compute the structure of the Jacobian and explain why it is sparse;
- derive the Gauss–Newton normal equations from a first-order Taylor approximation;
- explain Levenberg–Marquardt damping as a trust-region or regularized least-squares method;
- use block matrix algebra to derive the Schur complement that eliminates point variables;
- explain gauge freedom and why the Hessian is singular unless a coordinate gauge is fixed;
- interpret bundle adjustment statistically as maximum likelihood or MAP estimation;
- compare squared loss, Huber loss, Cauchy loss, and other robust M-estimators;
- connect classical bundle adjustment to differentiable optimization layers, neural rendering, learned correspondence uncertainty, and modern 3D foundation models.
18.2 17.1 Why bundle adjustment is needed
Chapters 14–16 developed the main geometric ingredients for 3D reconstruction:
- epipolar geometry estimates pairwise relative camera motion;
- stereo and triangulation estimate depth or sparse 3D points;
- Structure from Motion organizes many images into a global camera-and-point reconstruction.
Each step is approximate. Image features are noisy. Matches contain outliers. Essential matrices are estimated from finite samples. Triangulated points may be inconsistent across views. Incremental SfM accumulates drift. Bundle adjustment is the global refinement step that asks all cameras and all points to agree at once.
Suppose there are \(m\) cameras and \(n\) 3D points. Camera \(i\) has parameter vector
\[ \theta_i \in \mathbb{R}^{p_i}, \]
and point \(j\) has coordinates
\[ X_j \in \mathbb{R}^3. \]
The set of observed camera–point pairs is
\[ \mathcal{O}\subseteq \{1,\dots,m\}\times \{1,\dots,n\}. \]
For each observed pair \((i,j)\in\mathcal{O}\), we have an image measurement
\[ x_{ij}\in\mathbb{R}^2. \]
The camera model predicts the image point
\[ \widehat x_{ij}=h(\theta_i,X_j), \]
where \(h\) includes projection, intrinsics, distortion, and possibly other calibration parameters. The reprojection residual is
\[ r_{ij}(\theta_i,X_j)=x_{ij}-h(\theta_i,X_j)\in\mathbb{R}^2. \]
The basic bundle-adjustment objective is
\[ \min_{\theta_1,\dots,\theta_m,\,X_1,\dots,X_n} \frac{1}{2}\sum_{(i,j)\in\mathcal{O}} \|r_{ij}(\theta_i,X_j)\|_2^2. \]
This is called bundle adjustment because each 3D point generates a bundle of rays from multiple cameras, and the optimization adjusts the ray bundles so that they intersect the scene consistently.
Bundle adjustment is not a new camera model. It is the global nonlinear optimization problem that refines camera and scene variables by minimizing reprojection error.
Interactive 17.1. Reprojection error landscape. The objective is nonlinear. Near a good initialization, Gauss–Newton uses a local quadratic approximation; far away, the landscape may have valleys, flat directions, and local minima.
18.3 17.2 Reprojection residuals
The simplest calibrated perspective camera maps a world point \(X\in\mathbb{R}^3\) into camera coordinates by
\[ Y=RX+t, \]
where \(R\in SO(3)\) and \(t\in\mathbb{R}^3\). Write
\[ Y=(Y_1,Y_2,Y_3)^T. \]
The normalized image coordinates are
\[ \pi(Y)= \begin{bmatrix} Y_1/Y_3\\ Y_2/Y_3 \end{bmatrix}. \]
If the intrinsic matrix is
\[ K= \begin{bmatrix} f_x & s & c_x\\ 0 & f_y & c_y\\ 0 & 0 & 1 \end{bmatrix}, \]
then the predicted pixel coordinates are
\[ h(\theta,X)= \begin{bmatrix} f_xY_1/Y_3+sY_2/Y_3+c_x\\ f_yY_2/Y_3+c_y \end{bmatrix}. \]
The residual for an observed pixel \(x=(u,v)^T\) is
\[ r(\theta,X)= \begin{bmatrix} u\\v \end{bmatrix} - \begin{bmatrix} f_xY_1/Y_3+sY_2/Y_3+c_x\\ f_yY_2/Y_3+c_y \end{bmatrix}. \]
This residual is measured in pixels. A residual of length \(1\) means a one-pixel reprojection error. In practice, measurement noise may be anisotropic or may depend on the feature detector. If the covariance of \(x_{ij}\) is \(\Sigma_{ij}\), a statistically appropriate residual is the whitened residual
\[ \widetilde r_{ij}=\Sigma_{ij}^{-1/2}r_{ij}. \]
The corresponding objective is
\[ \min \frac{1}{2}\sum_{(i,j)\in\mathcal{O}} r_{ij}^T\Sigma_{ij}^{-1}r_{ij}. \]
This is weighted nonlinear least squares.
18.3.1 Residual dimension
Each observation contributes two scalar residuals. If there are \(N_{obs}=|\mathcal{O}|\) observations, the full residual vector has dimension
\[ r(z)\in\mathbb{R}^{2N_{obs}}, \]
where
\[ z=(\theta_1,\dots,\theta_m,X_1,\dots,X_n) \]
is the stacked optimization variable.
A typical reconstruction can have millions of residual components. The problem is large, but it is highly sparse.
18.4 17.3 Camera and pose parameterization
Bundle adjustment needs a finite-dimensional parameterization of camera geometry. The choice matters because rotations live on a curved manifold, not a Euclidean vector space.
18.4.1 Camera extrinsics
A camera pose can be written as a rigid transformation
\[ T=(R,t)\in SE(3), \]
where
\[ SE(3)=\left\{ \begin{bmatrix} R&t\\ 0&1 \end{bmatrix}:R\in SO(3),\;t\in\mathbb{R}^3 \right\}. \]
The rotation group is
\[ SO(3)=\{R\in\mathbb{R}^{3\times 3}:R^TR=I,\;\det R=1\}. \]
A small rotation update can be represented by a vector \(\omega\in\mathbb{R}^3\) through the exponential map
\[ R \leftarrow \exp([\omega]_\times)R, \]
where
\[ [\omega]_\times= \begin{bmatrix} 0&-\omega_3&\omega_2\\ \omega_3&0&-\omega_1\\ -\omega_2&\omega_1&0 \end{bmatrix}. \]
For a full rigid-motion update, write
\[ \xi=(\omega,v)\in\mathbb{R}^6, \]
and update the camera pose by
\[ T\leftarrow \exp(\xi^\wedge)T. \]
This local Euclidean update respects the geometry of \(SE(3)\) while allowing Gauss–Newton or Levenberg–Marquardt to solve for an increment \(\delta\xi\in\mathbb{R}^6\).
18.4.2 Camera intrinsics
A complete camera parameter vector may include
\[ \theta_i=(R_i,t_i,f_{x,i},f_{y,i},c_{x,i},c_{y,i},s_i,k_{1,i},k_{2,i},p_{1,i},p_{2,i},\dots). \]
In many SfM pipelines, intrinsics are fixed or partially fixed. For example, one may optimize focal length but keep principal point fixed. The more parameters are optimized, the more difficult the problem becomes and the more important calibration constraints become.
18.4.3 Point parameterization
The simplest point parameterization is
\[ X_j=(X_j,Y_j,Z_j)^T\in\mathbb{R}^3. \]
For some problems, especially monocular depth, it is better to use inverse depth
\[ \rho_j=\frac{1}{Z_j}. \]
Inverse depth is often numerically better when points are far away, because small changes in disparity correspond more linearly to changes in inverse depth than to changes in depth.
18.5 17.4 Nonlinear least squares
Let \(z\in\mathbb{R}^q\) denote all unknown parameters. The bundle-adjustment objective has the form
\[ F(z)=\frac{1}{2}\|r(z)\|_2^2 = \frac{1}{2}\sum_{k=1}^{2N_{obs}}r_k(z)^2. \]
Because projection contains divisions by depth and rotations are nonlinear, \(F\) is nonlinear.
At a current estimate \(z\), approximate the residual by first-order Taylor expansion:
\[ r(z+\delta)\approx r(z)+J(z)\delta, \]
where
\[ J(z)=\frac{\partial r}{\partial z} \in\mathbb{R}^{2N_{obs}\times q} \]
is the Jacobian matrix. The Gauss–Newton step solves
\[ \min_\delta \frac{1}{2}\|r+J\delta\|_2^2. \]
The normal equations are
\[ J^TJ\delta=-J^Tr. \]
Define
\[ H=J^TJ, \qquad b=J^Tr. \]
Then the step satisfies
\[ H\delta=-b. \]
The matrix \(H\) is the Gauss–Newton approximation to the Hessian of \(F\).
18.5.1 Why Gauss–Newton works well for reprojection problems
The exact Hessian of \(F\) is
\[ \nabla^2 F(z)=J^TJ+\sum_{k=1}^{2N_{obs}}r_k(z)\nabla^2 r_k(z). \]
Gauss–Newton drops the second term. This is a good approximation when either:
- the residuals are already small; or
- the residual functions are nearly linear in the current neighborhood.
Bundle adjustment usually starts from an SfM initialization, so residuals are often moderately small. This makes Gauss–Newton and Levenberg–Marquardt effective.
18.5.2 Levenberg–Marquardt damping
Pure Gauss–Newton can fail when the initialization is poor or when \(J^TJ\) is ill-conditioned. Levenberg–Marquardt solves the damped system
\[ (J^TJ+\lambda D)\delta=-J^Tr, \]
where \(\lambda>0\) is a damping parameter and \(D\) is often \(I\) or \(\operatorname{diag}(J^TJ)\). Large \(\lambda\) makes the step more gradient-descent-like; small \(\lambda\) makes the step more Gauss–Newton-like.
A useful equivalent formulation is
\[ \min_\delta \frac{1}{2}\|r+J\delta\|_2^2 + \frac{\lambda}{2}\|D^{1/2}\delta\|_2^2. \]
Thus damping is regularization of the local linearized step.
Interactive 17.2. Gauss–Newton versus damping. Damping changes the local quadratic step. With strong damping, the update moves cautiously along a stable direction; with weak damping, the update trusts the curvature model more strongly.
18.6 17.5 Projection Jacobians
The Jacobian is the mathematical object that connects geometry to optimization. For a calibrated camera, let
\[ Y=RX+t, \qquad Y=(Y_1,Y_2,Y_3)^T. \]
The normalized projection is
\[ \pi(Y)= \begin{bmatrix} Y_1/Y_3\\ Y_2/Y_3 \end{bmatrix}. \]
The derivative of \(\pi\) with respect to \(Y\) is
\[ \frac{\partial \pi}{\partial Y} = \begin{bmatrix} 1/Y_3 & 0 & -Y_1/Y_3^2\\ 0 & 1/Y_3 & -Y_2/Y_3^2 \end{bmatrix}. \]
If pixel intrinsics are included with no skew, then
\[ \frac{\partial h}{\partial Y} = \begin{bmatrix} f_x/Y_3 & 0 & -f_xY_1/Y_3^2\\ 0 & f_y/Y_3 & -f_yY_2/Y_3^2 \end{bmatrix}. \]
Because
\[ Y=RX+t, \]
the point derivative is
\[ \frac{\partial h}{\partial X} = \frac{\partial h}{\partial Y}R. \]
For pose derivatives, a small camera-frame perturbation gives approximately
\[ \delta Y \approx -[Y]_\times\delta\omega + \delta t. \]
Therefore the derivative with respect to a six-dimensional pose increment \(\delta\xi=(\delta\omega,\delta t)\) has the block form
\[ \frac{\partial h}{\partial \xi} = \frac{\partial h}{\partial Y} \begin{bmatrix} -[Y]_\times & I_3 \end{bmatrix}. \]
Since the residual is \(r=x-h\), the residual Jacobian is the negative of the prediction Jacobian:
\[ \frac{\partial r}{\partial z}=-\frac{\partial h}{\partial z}. \]
This sign does not change the least-squares objective, but it matters when writing the normal equations.
18.7 17.6 Sparsity of the Jacobian
A single residual \(r_{ij}\) depends only on camera \(i\) and point \(j\):
\[ r_{ij}=r_{ij}(\theta_i,X_j). \]
It does not depend on any other camera \(\theta_k\) with \(k\ne i\), and it does not depend on any other point \(X_\ell\) with \(\ell\ne j\). Therefore, each residual block has a Jacobian row block of the form
\[ \begin{bmatrix} 0&\cdots&A_{ij}&\cdots&0&0&\cdots&B_{ij}&\cdots&0 \end{bmatrix}, \]
where
\[ A_{ij}=\frac{\partial r_{ij}}{\partial \theta_i}, \qquad B_{ij}=\frac{\partial r_{ij}}{\partial X_j}. \]
This structure is the key to large-scale bundle adjustment. If camera \(i\) has \(p\) parameters, then
\[ A_{ij}\in\mathbb{R}^{2\times p}, \qquad B_{ij}\in\mathbb{R}^{2\times 3}. \]
The full Jacobian has many zeros. The normal matrix
\[ H=J^TJ \]
also has block sparsity. Camera variables become connected when they observe common points. Point variables are independent of each other except through cameras.
Interactive 17.3. Jacobian and normal-equation sparsity. Each observation touches exactly one camera block and one point block. This creates a sparse bipartite pattern in \(J\) and a block arrow pattern in \(J^TJ\).
18.8 17.7 Block normal equations
Separate camera increments and point increments:
\[ \delta= \begin{bmatrix} \delta_c\\ \delta_p \end{bmatrix}. \]
The normal equations have the block form
\[ \begin{bmatrix} U&W\\ W^T&V \end{bmatrix} \begin{bmatrix} \delta_c\\ \delta_p \end{bmatrix} = - \begin{bmatrix} b_c\\ b_p \end{bmatrix}. \]
Here:
- \(U\) is the camera–camera block;
- \(V\) is the point–point block;
- \(W\) couples cameras to points;
- \(b_c\) and \(b_p\) are gradient blocks.
The crucial property is that \(V\) is block diagonal. Each 3D point has its own small \(3\times 3\) block because different point variables do not directly appear in the same residual.
Therefore \(V^{-1}\) is cheap to apply: invert many small \(3\times 3\) matrices independently.
18.9 17.8 Schur complement elimination
From the block equations,
\[ U\delta_c+W\delta_p=-b_c, \]
and
\[ W^T\delta_c+V\delta_p=-b_p. \]
Assuming \(V\) is invertible, solve the second equation for \(\delta_p\):
\[ \delta_p=-V^{-1}b_p-V^{-1}W^T\delta_c. \]
Substitute into the first equation:
\[ U\delta_c+W(-V^{-1}b_p-V^{-1}W^T\delta_c)=-b_c. \]
Collect terms:
\[ (U-WV^{-1}W^T)\delta_c=-(b_c-WV^{-1}b_p). \]
The matrix
\[ S=U-WV^{-1}W^T \]
is the Schur complement. The reduced camera system is
\[ S\delta_c=-(b_c-WV^{-1}b_p). \]
After solving for \(\delta_c\), recover point increments by back-substitution:
\[ \delta_p=-V^{-1}(b_p+W^T\delta_c). \]
This is the central computational trick of bundle adjustment.
The Schur complement converts one huge camera-and-point system into a smaller camera-only system, using the fact that point blocks are conditionally independent given cameras.
Interactive 17.4. Schur complement intuition. Eliminating points creates direct camera–camera couplings whenever two cameras observe the same 3D point.
18.10 17.9 Gauge freedom and singular normal equations
Bundle adjustment has built-in coordinate ambiguity. If a reconstructed scene and camera trajectory explain the images, then a global similarity transformation also explains the images.
For calibrated SfM, the ambiguity is
\[ X_j' = sQX_j+a, \]
where \(s>0\), \(Q\in SO(3)\), and \(a\in\mathbb{R}^3\). Camera poses transform consistently. Image projections are unchanged.
This creates a seven-dimensional gauge freedom:
- 3 dimensions for global rotation;
- 3 dimensions for global translation;
- 1 dimension for global scale.
Therefore, the normal matrix has a null space unless the gauge is fixed. In practice one may:
- fix one camera pose and one scale variable;
- add weak priors on selected camera poses or points;
- impose a coordinate convention after each iteration;
- use damping and solve the system in a way that handles rank deficiency.
Gauge freedom is not a numerical bug. It is a mathematical property of the inverse problem.
18.10.1 Gauge and statistical uncertainty
Because of gauge freedom, absolute covariance is not identifiable unless a coordinate frame is chosen. Meaningful uncertainty statements should be gauge-invariant or should explicitly state the gauge convention.
For example, relative distance between two points can be meaningful after scale is fixed, but the absolute coordinate of one point depends on the chosen origin and orientation.
18.11 17.10 Robust bundle adjustment
Squared reprojection error is sensitive to outliers. A single wrong feature match may produce a huge residual and dominate the objective. Robust bundle adjustment replaces the squared loss with a robust penalty:
\[ \min_z \sum_{(i,j)\in\mathcal{O}}\rho\left(\|r_{ij}(z)\|_2^2\right). \]
Common choices include:
\[ \rho_{\text{quadratic}}(s)=s, \]
\[ \rho_{\text{Huber}}(s)= \begin{cases} s, & s\le c^2,\\ 2c\sqrt{s}-c^2, & s>c^2, \end{cases} \]
and
\[ \rho_{\text{Cauchy}}(s)=c^2\log\left(1+\frac{s}{c^2}\right). \]
A robust loss grows more slowly than \(s\) for large residuals. This reduces the influence of outliers.
18.11.1 Iteratively reweighted least squares
Many robust objectives can be solved by iteratively reweighted least squares. If \(s_k=\|r_k\|^2\), define a weight
\[ w_k=\rho'(s_k). \]
The local weighted least-squares problem is
\[ \min_\delta \frac{1}{2}\sum_k w_k\|r_k+J_k\delta\|^2. \]
Large residuals receive smaller weights. For the Huber loss, residuals smaller than \(c\) behave like ordinary least squares, while residuals larger than \(c\) receive reduced influence.
Interactive 17.5. Robust losses. Quadratic loss grows rapidly. Huber and Cauchy losses keep small residuals approximately quadratic but reduce the influence of large residuals.
18.12 17.11 Statistical interpretation
Assume each observation satisfies
\[ x_{ij}=h(\theta_i,X_j)+\varepsilon_{ij}, \]
where
\[ \varepsilon_{ij}\sim N(0,\sigma^2I_2) \]
independently. Then the likelihood is
\[ p(\{x_{ij}\}\mid z) \propto \exp\left( -\frac{1}{2\sigma^2}\sum_{(i,j)\in\mathcal{O}} \|x_{ij}-h(\theta_i,X_j)\|^2 \right). \]
Maximizing this likelihood is equivalent to minimizing the sum of squared reprojection errors.
If observation covariance is \(\Sigma_{ij}\), then the negative log-likelihood becomes
\[ \frac{1}{2}\sum_{(i,j)\in\mathcal{O}} r_{ij}^T\Sigma_{ij}^{-1}r_{ij} + ext{constant}. \]
Thus weighted bundle adjustment is maximum likelihood under Gaussian measurement noise.
18.12.1 MAP estimation
If we include priors, we obtain maximum a posteriori estimation. For example, a prior on camera pose increments might be
\[ \theta_i\sim N(\theta_i^0,\Lambda_i^{-1}). \]
This adds a quadratic penalty
\[ \frac{1}{2}(\theta_i-\theta_i^0)^T\Lambda_i(\theta_i-\theta_i^0). \]
Priors are useful when external sensors are available, such as GPS, IMU, wheel odometry, or known camera rig constraints.
18.12.2 Robust losses as heavy-tailed noise models
A robust loss corresponds to a non-Gaussian likelihood. For example, the Cauchy loss is related to a heavy-tailed distribution in which large residuals are more probable than under a Gaussian model. This explains why robust bundle adjustment is statistically reasonable when matches contain outliers.
18.13 17.12 Uncertainty and covariance
Near a local optimum \(\widehat z\), approximate the residual model by
\[ r(z)\approx r(\widehat z)+J(\widehat z)(z-\widehat z). \]
Under Gaussian noise with variance \(\sigma^2\), the approximate covariance of the estimator is
\[ \operatorname{Cov}(\widehat z) \approx \sigma^2(J^TJ)^{-1}. \]
More generally, with weights,
\[ \operatorname{Cov}(\widehat z) \approx (J^TWJ)^{-1}. \]
This formula must be interpreted carefully because of gauge freedom. After gauge fixing, covariance can describe uncertainty in camera poses, 3D points, or derived quantities.
18.13.1 Point uncertainty
A point seen by many cameras from diverse viewpoints usually has small uncertainty. A point seen by only two nearby cameras may have large depth uncertainty. This matches the stereo formula from Chapter 15: small baseline leads to unstable depth.
18.13.2 Camera uncertainty
A camera observing many well-distributed points has better pose conditioning. If all points lie on a plane, or if the camera motion is close to pure rotation, some parameters become weakly determined.
Interactive 17.6. Baseline and point uncertainty. Larger triangulation angle usually reduces depth uncertainty, while small baseline creates a long uncertainty direction along the viewing ray.
18.14 17.13 Conditioning and scaling
Bundle adjustment often contains variables with different units and scales:
- pixel residuals are measured in pixels;
- translations may be measured in meters or arbitrary SfM units;
- rotations are measured in radians;
- focal lengths may be hundreds or thousands of pixels;
- radial distortion parameters may be small numbers.
Poor scaling can make \(J^TJ\) ill-conditioned. Common remedies include:
- normalize image coordinates before initialization;
- scale 3D points to have reasonable magnitude;
- use diagonal damping \(D=\operatorname{diag}(J^TJ)\);
- choose minimal and well-scaled parameterizations;
- avoid optimizing poorly identifiable parameters unless data support them.
18.14.1 Degenerate configurations
Bundle adjustment cannot recover information that is absent from the data. Common degeneracies include:
- pure camera rotation with no translation;
- all points on a single plane;
- very small baseline;
- too few feature tracks;
- repeated texture causing ambiguous correspondences;
- incorrect intrinsics absorbed by structure deformation;
- rolling-shutter effects modeled as global-shutter cameras.
A clean mathematical formulation helps diagnose these failures: degeneracy appears as rank deficiency, high condition number, large covariance, or unstable updates.
18.15 17.14 Algorithmic structure
A typical Levenberg–Marquardt bundle-adjustment iteration is:
- project all visible points into all observing cameras;
- compute residuals \(r_{ij}\);
- compute Jacobian blocks \(A_{ij}\) and \(B_{ij}\);
- apply robust weights if needed;
- assemble block normal equations;
- add damping;
- eliminate point variables by Schur complement;
- solve the reduced camera system;
- recover point updates by back-substitution;
- update cameras and points;
- accept or reject the step based on objective decrease.
In pseudocode:
Given initial cameras theta_i and points X_j
repeat
for each observation (i,j):
predict xhat_ij = h(theta_i, X_j)
compute residual r_ij = x_ij - xhat_ij
compute Jacobian blocks A_ij, B_ij
compute robust weight w_ij
assemble U, V, W, b_c, b_p
add Levenberg--Marquardt damping
form S = U - W V^{-1} W^T
solve S delta_c = -(b_c - W V^{-1} b_p)
recover delta_p = -V^{-1}(b_p + W^T delta_c)
update cameras and points
until convergenceThe expensive operation is solving the reduced camera system. Large systems use sparse Cholesky, preconditioned conjugate gradients, or hybrid direct-iterative methods.
18.16 17.15 Python example: synthetic reprojection residuals
The following example constructs a simple camera projection model and computes reprojection residuals. It is intentionally small so that the mathematical objects are visible.
18.17 17.16 Python example: finite-difference Jacobian
For teaching purposes, finite differences are a useful way to check analytic derivatives.
import numpy as np
def residual_from_point(X, x_obs, R, t, K):
return x_obs - project_point(R, t, K, X)
def finite_difference_jacobian(fun, x, eps=1e-6):
f0 = fun(x)
J = np.zeros((len(f0), len(x)))
for k in range(len(x)):
xp = x.copy()
xm = x.copy()
xp[k] += eps
xm[k] -= eps
J[:, k] = (fun(xp) - fun(xm)) / (2 * eps)
return J
x_obs = np.array([460.7, 279.6])
X0 = np.array([0.5, 0.3, 4.0])
J_fd = finite_difference_jacobian(
lambda X: residual_from_point(X, x_obs, R, t, K),
X0
)
print(J_fd)18.18 17.17 Python example: a small least-squares problem
SciPy can solve small bundle-adjustment-like problems. The example below estimates 3D points when the cameras are known. Full bundle adjustment additionally optimizes camera variables.
import numpy as np
from scipy.optimize import least_squares
rng = np.random.default_rng(7243)
K = np.array([[700.0, 0.0, 320.0],
[0.0, 700.0, 240.0],
[0.0, 0.0, 1.0]])
camera_ts = [
np.array([-0.4, 0.0, 0.0]),
np.array([0.0, 0.0, 0.0]),
np.array([0.4, 0.0, 0.0]),
]
R = np.eye(3)
points_true = rng.normal(size=(20, 3))
points_true[:, 2] += 5.0
observations = []
for i, t in enumerate(camera_ts):
for j, X in enumerate(points_true):
x = project_point(R, t, K, X)
x_noisy = x + rng.normal(scale=0.5, size=2)
observations.append((i, j, x_noisy))
points_init = points_true + rng.normal(scale=0.15, size=points_true.shape)
def residual_vector(flat_points):
P = flat_points.reshape((-1, 3))
residuals = []
for i, j, x_obs in observations:
x_hat = project_point(R, camera_ts[i], K, P[j])
residuals.extend(x_obs - x_hat)
return np.array(residuals)
result = least_squares(residual_vector, points_init.ravel(), method="trf")
points_est = result.x.reshape((-1, 3))
print("initial RMSE:", np.sqrt(np.mean(residual_vector(points_init.ravel())**2)))
print("final RMSE:", np.sqrt(np.mean(residual_vector(result.x)**2)))
print("mean 3D error:", np.mean(np.linalg.norm(points_est - points_true, axis=1)))18.19 17.18 Python example: sparse Jacobian pattern
For large problems, sparse structure is as important as the formulas.
import numpy as np
from scipy.sparse import lil_matrix
import matplotlib.pyplot as plt
n_cameras = 5
n_points = 12
camera_dim = 6
point_dim = 3
rng = np.random.default_rng(17)
observations = []
for i in range(n_cameras):
visible = rng.choice(n_points, size=7, replace=False)
for j in visible:
observations.append((i, j))
m = 2 * len(observations)
n = camera_dim * n_cameras + point_dim * n_points
J_pattern = lil_matrix((m, n), dtype=int)
for obs_index, (i, j) in enumerate(observations):
rows = slice(2 * obs_index, 2 * obs_index + 2)
cam_cols = slice(camera_dim * i, camera_dim * (i + 1))
pt_offset = camera_dim * n_cameras
pt_cols = slice(pt_offset + point_dim * j, pt_offset + point_dim * (j + 1))
J_pattern[rows, cam_cols] = 1
J_pattern[rows, pt_cols] = 1
plt.figure(figsize=(8, 4))
plt.spy(J_pattern, markersize=3)
plt.title("Bundle-adjustment Jacobian sparsity pattern")
plt.xlabel("parameters")
plt.ylabel("residual components")
plt.show()18.20 17.19 TensorFlow example: differentiable reprojection loss
Modern AI systems often include geometry losses inside neural training. The following sketch shows a differentiable reprojection loss. The network might predict depth, pose, uncertainty weights, or feature correspondences.
import tensorflow as tf
def tf_project(points, R, t, K):
"""points: [N, 3], R: [3, 3], t: [3], K: [3, 3]."""
Y = tf.linalg.matmul(points, R, transpose_b=True) + t[None, :]
x = Y[:, :2] / Y[:, 2:3]
fx = K[0, 0]
fy = K[1, 1]
cx = K[0, 2]
cy = K[1, 2]
u = fx * x[:, 0] + cx
v = fy * x[:, 1] + cy
return tf.stack([u, v], axis=1)
points = tf.Variable([[0.4, 0.2, 4.0],
[-0.3, 0.1, 5.0]], dtype=tf.float32)
R = tf.eye(3, dtype=tf.float32)
t = tf.constant([0.1, 0.0, 0.0], dtype=tf.float32)
K = tf.constant([[700.0, 0.0, 320.0],
[0.0, 700.0, 240.0],
[0.0, 0.0, 1.0]], dtype=tf.float32)
x_obs = tf.constant([[407.0, 275.0],
[292.0, 254.0]], dtype=tf.float32)
with tf.GradientTape() as tape:
x_hat = tf_project(points, R, t, K)
residuals = x_obs - x_hat
loss = tf.reduce_mean(tf.reduce_sum(residuals**2, axis=1))
grad_points = tape.gradient(loss, points)
print("loss:", loss.numpy())
print("gradient wrt points:\n", grad_points.numpy())This is not a replacement for sparse second-order bundle adjustment. Instead, it shows how reprojection error becomes a differentiable learning signal.
18.21 17.20 AI components and modern connections
Bundle adjustment appears in modern AI systems in several forms.
18.21.1 Learned correspondences with classical optimization
Modern feature matchers may produce:
\[ (x_{ij},x_{kj},s_j), \]
where \(s_j\) is a confidence score. Bundle adjustment can convert confidence into weights:
\[ w_j=g(s_j), \]
and minimize
\[ \sum_j w_j\|r_j\|^2. \]
The learned model handles appearance; the geometric optimizer enforces multi-view consistency.
18.21.2 Differentiable bundle adjustment
A differentiable BA layer computes
\[ z^*(\alpha)=\arg\min_z F(z,\alpha), \]
where \(\alpha\) are network outputs such as correspondences, depths, or weights. Gradients can be obtained by unrolling optimization iterations or by implicit differentiation.
At an optimum,
\[ \nabla_zF(z^*(\alpha),\alpha)=0. \]
Implicit differentiation gives
\[ \frac{\partial z^*}{\partial \alpha} = - \left(\nabla_{zz}^2F\right)^{-1} \nabla_{z\alpha}^2F. \]
This equation shows why differentiable optimization is mathematically connected to Hessians, Schur complements, and linear solves.
18.21.3 Neural rendering and camera refinement
Neural rendering methods often optimize camera poses together with scene parameters. The residual may be photometric rather than sparse geometric:
\[ r_{ip}=I_i(p)-\widehat I_i(p;\theta,\phi), \]
where \(\phi\) denotes a neural scene representation. This is a dense cousin of bundle adjustment. The same ideas appear: reprojection, nonlinear least squares, gauge, conditioning, robust losses, and uncertainty.
18.21.4 Learned uncertainty
A network may predict an observation variance \(\sigma_{ij}^2\). Then the loss becomes
\[ \sum_{(i,j)} \left( \frac{\|r_{ij}\|^2}{2\sigma_{ij}^2} + \log \sigma_{ij} \right). \]
The first term downweights uncertain observations; the second prevents the network from setting all uncertainties to infinity.
18.22 17.21 Common failure modes
Bundle adjustment is powerful but not magical. It refines a model; it cannot fix a completely wrong model.
18.22.1 Poor initialization
Nonlinear least squares is local. If the initial cameras and points are far from the correct basin, the optimizer can converge to a wrong local minimum or diverge.
18.22.2 Outlier correspondences
Wrong matches can create large residuals. Robust losses help, but too many outliers can still corrupt the reconstruction.
18.22.3 Incorrect camera model
If the real camera has radial distortion, rolling shutter, or changing intrinsics, but the optimizer assumes a simple pinhole camera, reprojection residuals may be systematically biased.
18.22.4 Weak geometry
Small baseline, planar scenes, and poor feature distribution create ill-conditioned systems. The optimizer may reduce residuals while producing unstable 3D structure.
18.22.5 Overparameterization
Optimizing every intrinsic and distortion parameter without enough data can lead to parameter coupling. For example, focal length, depth scale, and radial distortion can partially compensate for each other.
18.23 17.22 Summary
Bundle adjustment is the mathematical center of geometric reconstruction. Its core ingredients are:
- residuals from image measurements to projected 3D points;
- nonlinear least squares over camera and point variables;
- Jacobian blocks from projection derivatives;
- sparse normal equations induced by the camera–point observation graph;
- Schur complement elimination of point variables;
- gauge fixing for coordinate ambiguity;
- robust losses for outlier resistance;
- statistical interpretation as MLE or MAP estimation;
- uncertainty estimates from local curvature;
- differentiable optimization connections to modern AI.
A student who understands bundle adjustment understands how linear algebra, geometry, statistics, and optimization meet in modern computer vision.
18.24 Exercises
18.24.1 Conceptual exercises
- Explain why minimizing 3D point error is usually less appropriate than minimizing reprojection error.
- What does it mean for a residual to depend only on one camera and one point? Why does this produce sparsity?
- Explain gauge freedom in calibrated SfM. Why does it have seven dimensions?
- Why can bundle adjustment reduce drift in incremental Structure from Motion?
- Give two examples where low reprojection error does not guarantee a correct reconstruction.
- Explain the difference between a camera parameter update in Euclidean coordinates and an update on \(SO(3)\) or \(SE(3)\).
- Why does robust loss improve stability when some correspondences are wrong?
- Explain why planar scenes can be degenerate for camera reconstruction.
18.24.2 Mathematical exercises
- Derive the derivative of \(u=fX/Z+c_x\) and \(v=fY/Z+c_y\) with respect to \((X,Y,Z)\).
- Starting from \(\min_\delta \|r+J\delta\|^2/2\), derive the Gauss–Newton normal equations.
- Derive the Levenberg–Marquardt equation from the regularized local problem \[ \min_\delta \frac{1}{2}\|r+J\delta\|^2+ rac{\lambda}{2}\|\delta\|^2. \]
- For the block system \[ \begin{bmatrix}U&W\\W^T&V\end{bmatrix} \begin{bmatrix}\delta_c\\\delta_p\end{bmatrix} =- \begin{bmatrix}b_c\\b_p\end{bmatrix}, \] derive the Schur complement camera system.
- Show that if two cameras observe no common points, their off-diagonal Schur complement block is zero.
- Suppose \(\varepsilon\sim N(0,\sigma^2I)\). Show that maximum likelihood gives ordinary least squares.
- Derive the Huber weight used in iteratively reweighted least squares.
- Explain why \(J^TJ\) is singular when the global coordinate system is not fixed.
- Derive the approximate covariance formula \(\sigma^2(J^TJ)^{-1}\) from the linearized residual model.
- For a point seen by two cameras, explain mathematically why small triangulation angle gives large depth uncertainty.
18.24.3 Python exercises
- Implement projection and reprojection residuals for three cameras and twenty points.
- Use finite differences to verify the analytic point Jacobian for the pinhole camera.
- Use
scipy.optimize.least_squaresto estimate 3D points with fixed cameras. - Extend the previous exercise to estimate translations while keeping rotations fixed.
- Construct the sparse Jacobian pattern for a synthetic camera–point observation graph.
- Implement Huber weights and compare least squares with robust least squares when some observations are outliers.
- Simulate small baseline and large baseline stereo. Compare point uncertainty numerically.
- Use the Schur complement to solve a small block system and compare with solving the full system directly.
- Plot objective value versus iteration for Gauss–Newton and damped Gauss–Newton on a toy problem.
- Add Gaussian priors on camera positions and show how they change the normal equations.
18.24.4 AI/modeling exercises
- Ask an AI assistant to explain bundle adjustment in one paragraph. Then rewrite the explanation using residuals, Jacobians, and Schur complements.
- Ask an AI assistant to generate pseudocode for Levenberg–Marquardt bundle adjustment. Check whether it correctly handles damping and step acceptance.
- Ask an AI assistant to compare bundle adjustment and neural rendering optimization. Identify the shared mathematical structure.
- Design a learned feature-matching pipeline that outputs both correspondences and uncertainty weights for robust BA.
- Ask an AI assistant for possible degeneracies in bundle adjustment. For each answer, write the corresponding rank, conditioning, or gauge issue.
- Propose a differentiable BA layer inside a neural network. State what is differentiated, what is optimized internally, and what loss supervises the outer network.
18.25 References and further reading
Bundle adjustment is surveyed in Triggs et al. (2000) and is developed as part of multi-view geometry in Hartley and Zisserman (2004), Faugeras (1993), Ma et al. (2004), and Szeliski (2022). Large-scale SfM systems and practical sparse reconstruction pipelines are represented by Snavely et al. (2006), Agarwal et al. (2011), and Sch"onberger and Frahm (2016). The numerical optimization background includes Gauss–Newton, Levenberg–Marquardt, sparse least squares, and trust-region methods; see Levenberg (1944), Marquardt (1963), and Nocedal and Wright (2006). Robust estimation and M-estimators are represented by Huber (1964) and Black and Anandan (1996). Practical implementations of large sparse nonlinear least squares include Agarwal et al. (2012).
Agarwal, Sameer, Yasutaka Furukawa, Noah Snavely, et al. 2011. “Building Rome in a Day.” Communications of the ACM 54 (10): 105–12.
Agarwal, Sameer, Keir Mierle, et al. 2012. “Ceres Solver: Tutorial and Reference.” Google Research.
Black, Michael J., and P. Anandan. 1996. “The Robust Estimation of Multiple Motions: Parametric and Piecewise-Smooth Flow Fields.” Computer Vision and Image Understanding 63 (1): 75–104.
Faugeras, Olivier. 1993. Three-Dimensional Computer Vision: A Geometric Viewpoint. MIT Press.
Hartley, Richard, and Andrew Zisserman. 2004. Multiple View Geometry in Computer Vision. 2nd ed. Cambridge University Press. https://www.robots.ox.ac.uk/~vgg/hzbook/.
Huber, Peter J. 1964. “Robust Estimation of a Location Parameter.” The Annals of Mathematical Statistics 35 (1): 73–101.
Levenberg, Kenneth. 1944. “A Method for the Solution of Certain Non-Linear Problems in Least Squares.” Quarterly of Applied Mathematics 2 (2): 164–68.
Ma, Yi, Stefano Soatto, Jana Kosecka, and Shankar S. Sastry. 2004. An Invitation to 3-d Vision: From Images to Geometric Models. Springer.
Marquardt, Donald W. 1963. “An Algorithm for Least-Squares Estimation of Nonlinear Parameters.” Journal of the Society for Industrial and Applied Mathematics 11 (2): 431–41.
Nocedal, Jorge, and Stephen J. Wright. 2006. Numerical Optimization. 2nd ed. Springer.
Sch"onberger, Johannes L., and Jan-Michael Frahm. 2016. “Structure-from-Motion Revisited.” IEEE Conference on Computer Vision and Pattern Recognition, 4104–13.
Szeliski, Richard. 2022. Computer Vision: Algorithms and Applications. 2nd ed. Springer. https://szeliski.org/Book/.
Triggs, Bill, Philip F. McLauchlan, Richard I. Hartley, and Andrew W. Fitzgibbon. 2000. “Bundle Adjustment – a Modern Synthesis.” Vision Algorithms: Theory and Practice, 298–372.