Bundle Adjustment (BA) Overview

• Bundle adjustment (BA) is a nonlinear optimization technique that simultaneously refines camera parameters and 3D structure points, minimizing reprojection errors in visual data, especially useful in structure from motion (SfM) and SLAM.
• BA leverages sparse camera–landmark incidence graphs, enabling efficient solving using techniques like Gauss–Newton, Levenberg–Marquardt, and distributed optimization, tailored to handle large-scale datasets in real-time across various applications such as photogrammetry and sensor fusion.
• Successful BA addresses key challenges like gauge freedoms, linearity of residuals, and robust estimation to enhance accuracy, usage examples include satellite imagery registration, robust to scenarios with outliers

Bundle adjustment (BA) is a nonlinear optimization problem that jointly refines camera or sensor parameters and scene structure so that measured observations agree with projections or geometric predictions. In classical visual BA, the unknowns are camera parameters and 3D points, and the residuals are reprojection errors over the observed camera–point pairs. BA is central to structure from motion (SfM), visual simultaneous localization and mapping (SLAM), photogrammetry, camera calibration, robotic localization, 3D reconstruction, satellite-image registration, and sensor-fusion systems. Its characteristic computational structure is a sparse camera–landmark incidence graph: each observation depends on one camera and one scene element. This structure enables sparse Levenberg–Marquardt (LM), Gauss–Newton (GN), Schur-complement elimination, iterative linear solvers, distributed optimization, specialized hardware, and reduced geometric formulations.

1. Mathematical formulation and geometric structure

Let $m$ images or sensor poses observe $n$ scene points. Denote the parameters of camera $j$ by $x_j$, the world coordinates of point $i$ by $y_i\in\mathbb{R}^3$, and the measured image position by $z_{ij}\in\mathbb{R}^2$. If $h(x_j,y_i)$ is the projection function, the observation model is

$z_{ij}=h(x_j,y_i)+\epsilon_{ij},$

with reprojection residual

$\epsilon_{ij}=z_{ij}-h(x_j,y_i).$

The set of valid observations is commonly represented by $n$0, where $n$1 means that point $n$2 is observed in image $n$3. A conventional Gaussian model assumes

$n$4

With optional Gaussian priors on cameras and points, the maximum a posteriori estimate minimizes a weighted nonlinear least-squares objective,

$n$5

If no prior information is available, the corresponding precision matrix is set to zero. Priors or explicit constraints can also resolve gauge freedoms: without gauge fixing, translations, rotations, and—in monocular reconstruction—scale transformations can leave the reprojection geometry unchanged.

For a calibrated pinhole camera, a world point is transformed by

$n$6

and projected after division by depth. With intrinsic matrix $n$7, the camera model may include focal lengths, principal point, skew, and lens-distortion parameters. In monocular visual SLAM, the camera pose consists of $n$8 and $n$9, while the reconstructed points are $j$0. The standard objective is

$j$1

where the sum is restricted to valid feature-image observations (Parra et al., 2019).

The Jacobian has a bipartite block structure. If

$j$2

each observation contributes only one camera block and one point block. Derivatives with respect to unrelated cameras or points are zero. This sparsity is the defining algebraic property of BA and underlies both classical and specialized solvers.

2. Linearization, LM optimization, and gauge handling

At an iterate $j$3, where $j$4 stacks all camera and structure variables, nonlinear residuals are linearized as

$j$5

GN uses the normal equations

$j$6

LM adds damping,

$j$7

where $j$8 may be the identity or a diagonal scaling derived from $j$9. Large damping produces gradient-descent-like behavior, while small damping approaches GN. Trial steps are accepted or rejected according to the reduction in the nonlinear objective, and the damping parameter is adjusted accordingly.

The damped normal matrix may be partitioned into camera and point variables:

$x_j$0

The point block $x_j$1 is generally block diagonal because each point contributes independently before camera coupling is induced through shared observations. Eliminating point increments gives the camera Schur system

$x_j$2

After solving for $x_j$3, the point increment is recovered by back-substitution. The reduced system is governed by the camera covisibility graph: two cameras are coupled when they observe at least one common point.

Gauge freedoms can make the normal matrix singular or nearly singular. Common remedies include fixing one camera pose, fixing a second camera orientation or scale, imposing priors, fixing reference points, or projecting out null-space directions. In monocular BA, the similarity gauge consists of three global rotations, three global translations, and one global scale. The multigrid formulation of BA explicitly models these seven modes as near-nullspace directions (Konolige et al., 2020).

BA is nonconvex because of camera rotations, perspective division, and the joint camera–structure parameterization. GN and LM therefore provide local optimization methods rather than global-optimality guarantees. Initialization from feature matching, pose estimation, triangulation, tracking, and map construction is consequently important. In monocular SLAM, poor local BA solutions can affect subsequent feature filtering and map maintenance (Parra et al., 2019).

3. Robust estimation and alternatives to the Gaussian objective

The quadratic Gaussian penalty grows as the square of the residual. A single mismatched feature with a large reprojection error can therefore dominate many correct observations. Because camera and point variables are coupled, the optimizer may distort poses and structure in an attempt to reduce the residual of an incorrect correspondence.

A Student’s-$x_j$4 robust formulation replaces the Gaussian model with a heavy-tailed distribution. For Mahalanobis distance

$x_j$5

the negative log-likelihood, up to parameter-independent constants, is

$x_j$6

where $x_j$7 is the degrees-of-freedom parameter and $x_j$8 is the residual dimension. Large $x_j$9 approaches Gaussian behavior; smaller $i$0 produces heavier tails and reduces the influence of large residuals. The resulting robust objective uses logarithmic rather than quadratic growth for very large residuals. Iteratively reweighted least-squares subproblems assign small weights to large residuals while retaining the observations in the optimization.

The Student’s-$i$1 robust bundle adjustment algorithm, RST-BA, applies this model to reprojection errors and also to camera and point priors. It preserves the camera–point sparsity pattern, uses an implicit trust-region method related to LM, has the same asymptotic time complexity as conventional $i$2-BA, and requires only limited modifications to a least-squares implementation (Aravkin et al., 2011). In the reported synthetic experiments, RST-BA was substantially more accurate than $i$3-BA and $i$4 editing under contaminated Gaussian and Student’s-$i$5 errors. With $i$6 contamination from a high-variance $i$7 component, relative world-point MSEs were approximately $i$8, $i$9, and $y_i\in\mathbb{R}^3$0 for $y_i\in\mathbb{R}^3$1-BA, $y_i\in\mathbb{R}^3$2-BA, and RST-BA, respectively; corresponding camera-coordinate errors were approximately $y_i\in\mathbb{R}^3$3, $y_i\in\mathbb{R}^3$4, and $y_i\in\mathbb{R}^3$5. On unprocessed Apollo 15 imagery containing nearly $y_i\in\mathbb{R}^3$6 outliers, RST-BA produced a lunar DEM without ground-control points, whereas conventional $y_i\in\mathbb{R}^3$7-BA failed.

Other robust losses are used in distributed BA. A distributed formulation based on alternating direction method of multipliers (ADMM) compares squared $y_i\in\mathbb{R}^3$8 and Huber reprojection losses. Huber penalties are useful not only for incorrect correspondences but also because each local subproblem may contain too few observations to estimate camera and point variables stably. Extreme local estimates can act as algorithmic outliers during consensus even when measurements are Gaussian (Ramamurthy et al., 2017).

Robust losses do not resolve every failure mode. Systematic calibration errors, incorrect camera models, rolling-shutter effects, structured mismatches, gauge ambiguity, and poor initialization may not be adequately handled by residual downweighting alone. Robustness is therefore distinct from correspondence generation, geometric validity, and observability.

4. Sparse, distributed, and reduced formulations

Schur-complement and iterative solvers

Since the number of 3D points is often much larger than the number of cameras, point elimination is usually the main reduction. Simplified complexity descriptions characterize dense formulations as approximately $y_i\in\mathbb{R}^3$9 and Schur-reduced systems as approximately $z_{ij}\in\mathbb{R}^2$0, although actual cost depends on visibility sparsity, fill-in, block dimensions, and solver choice (Chen et al., 2019).

For very large camera graphs, conjugate gradient (CG) and preconditioned CG (PCG) avoid explicit factorization of the reduced system. Preconditioners include block-Jacobi, visibility-based block-Jacobi, visibility-based block-tridiagonal, SSOR-style, and QR-based constructions. Multigrid methods address the global low-frequency modes that local preconditioners handle poorly. An unsmoothed aggregation multigrid preconditioner uses camera visibility similarity to form aggregates and explicitly represents the seven similarity-gauge modes together with nine camera-parameter constant modes. On challenging synthetic city-scale problems, it achieved up to $z_{ij}\in\mathbb{R}^2$1 speedup over point-block Jacobi and up to $z_{ij}\in\mathbb{R}^2$2 over visibility-based preconditioners (Konolige et al., 2020).

Distributed BA

Distributed BA replicates camera and point variables locally for each observation or observation block and imposes consensus constraints. In the ADMM formulation, $z_{ij}\in\mathbb{R}^2$3 is a local copy of point $z_{ij}\in\mathbb{R}^2$4 for camera $z_{ij}\in\mathbb{R}^2$5, and $z_{ij}\in\mathbb{R}^2$6 is a local copy of camera $z_{ij}\in\mathbb{R}^2$7 for point $z_{ij}\in\mathbb{R}^2$8. Consensus constraints require

$z_{ij}\in\mathbb{R}^2$9

The local observation problems are independent, while global point and camera values are obtained by averaging shifted local estimates. Dual variables accumulate consensus violations. The local-update stage has $h(x_j,y_i)$0 work per iteration, where

$h(x_j,y_i)$1

is the number of observations; consensus updates require $h(x_j,y_i)$2 work, and dual updates require $h(x_j,y_i)$3 work. Communication remains a limiting factor because replicated variables must be exchanged across processors (Ramamurthy et al., 2017).

The convergence analysis is for a nonconvex and nonsmooth ADMM framework. Under bounded projection depth, a convex smooth measurement loss, squared Euclidean consensus penalties, and sufficiently large penalty parameters, the method converges to a stationary point of the augmented Lagrangian. It does not guarantee a global optimum.

A related distributed treatment uses partial Douglas–Rachford splitting. Cameras are partitioned across workers, each worker maintains local copies of shared 3D points, and consensus averages those copies. This approach preserves local BA structure but remains subject to nonconvexity, communication cost, gauge handling, and initialization dependence (Chen et al., 2019).

Structure compression and pointless BA

Several formulations reduce or eliminate explicit scene variables. Pointless Global BA computes local pair or triplet adjustments, eliminates local points with a Schur complement, and retains the camera-space reduced Hessian. If $h(x_j,y_i)$4 is the relative motion induced by global camera poses, $h(x_j,y_i)$5 is the locally estimated motion, and $h(x_j,y_i)$6 is the local reduced Hessian, the global objective is

$h(x_j,y_i)$7

Unlike unweighted motion averaging, this formulation retains anisotropic uncertainty, camera-parameter correlations, feature geometry, and depth observability through $h(x_j,y_i)$8. On a $h(x_j,y_i)$9-image photogrammetric dataset, it achieved a final reprojection error of approximately $z_{ij}=h(x_j,y_i)+\epsilon_{ij},$0 compared with $z_{ij}=h(x_j,y_i)+\epsilon_{ij},$1 for MicMac BA, while using approximately $z_{ij}=h(x_j,y_i)+\epsilon_{ij},$2k rather than $z_{ij}=h(x_j,y_i)+\epsilon_{ij},$3k parameters (Rupnik et al., 2023).

5. Extensions beyond point-reprojection BA

Photometric, semantic, and event-based BA

Photometric BA minimizes intensity consistency across warped images rather than only feature reprojection errors. Semantic photometric BA incorporates a learned 3D object prior into photometric reconstruction, addressing limitations associated with incomplete visibility and weak surface texture; the supplied source does not provide sufficient technical detail to verify its exact objective or optimization procedure (Zhu et al., 2017).

Photometric LiDAR and RGB-D BA represents observations as intensity, depth or range, and surface-normal images. It optimizes sensor poses using differentiable warping and Huber residuals without explicit landmarks or surfels. RGB-D uses pinhole projection and optical-axis depth, whereas LiDAR uses spherical projection and Euclidean range. The system supports coupled or consecutive LiDAR–RGB-D refinement and uses coarse-to-fine LM optimization (Giammarino et al., 2023).

Event-based Photometric BA directly uses the event-generation model. For event $z_{ij}=h(x_j,y_i)+\epsilon_{ij},$4, the log-brightness relation is

$z_{ij}=h(x_j,y_i)+\epsilon_{ij},$5

Under pure rotational motion, the method jointly optimizes camera orientations and a semi-dense panoramic brightness map. Events are warped into the map using camera calibration and orientation interpolation, without conversion into frames, time surfaces, voxel grids, or other image-like representations. The reported experiments include synthetic and real event-camera data and reductions in photometric error of up to $z_{ij}=h(x_j,y_i)+\epsilon_{ij},$6 (Guo et al., 2024).

Planar and LiDAR BA

Planar Bundle Adjustment jointly optimizes depth-sensor poses and plane parameters using point-to-plane residuals. A plane observation may contain hundreds or thousands of points, making direct BA inefficient. Because the point coordinates enter only through $z_{ij}=h(x_j,y_i)+\epsilon_{ij},$7, $z_{ij}=h(x_j,y_i)+\epsilon_{ij},$8, $z_{ij}=h(x_j,y_i)+\epsilon_{ij},$9, and $\epsilon_{ij}=z_{ij}-h(x_j,y_i).$0, a QR factorization compresses each pose–plane observation to four equivalent residual rows. The reduced Jacobian and residual preserve exactly $\epsilon_{ij}=z_{ij}-h(x_j,y_i).$1 and $\epsilon_{ij}=z_{ij}-h(x_j,y_i).$2, and therefore produce the same LM step as the original point-by-point formulation (Zhou et al., 2020).

BALM eliminates plane and edge parameters analytically. For a transformed feature point set with covariance matrix $\epsilon_{ij}=z_{ij}-h(x_j,y_i).$3, plane cost is $\epsilon_{ij}=z_{ij}-h(x_j,y_i).$4 and edge cost is $\epsilon_{ij}=z_{ij}-h(x_j,y_i).$5. The optimization therefore contains only scan poses. Adaptive voxelization creates local geometric features, and closed-form first- and second-order derivatives support LM. The method operates as a sliding-window back end for LOAM and reports near-real-time optimization of 20 scans (Liu et al., 2020).

A later point-cluster formulation stores each feature’s local point statistics in a $\epsilon_{ij}=z_{ij}-h(x_j,y_i).$6 matrix containing second-order moments, first-order moments, and point count. Rigid transformations and cluster merging become matrix operations, so optimization, derivatives, and uncertainty estimation avoid enumerating raw points. On 19 real sequences, the reported average absolute trajectory RMSE was $\epsilon_{ij}=z_{ij}-h(x_j,y_i).$7 m for plane-plus-edge features and $\epsilon_{ij}=z_{ij}-h(x_j,y_i).$8 m for plane-only features, with average optimization times of $\epsilon_{ij}=z_{ij}-h(x_j,y_i).$9 s and $n$00 s, respectively (Liu et al., 2022).

Satellite RPC BA

For high-resolution satellite imagery, BA often refines image-space RPC biases and object-space tie-point coordinates. A large-scale RPC formulation estimates two image biases per image and eliminates the three-dimensional object-point corrections analytically, leaving a reduced system whose dominant dimension depends on the number of images rather than the number of tie points. Plane rectification, epipolar-constrained matching, multi-block Census descriptors, multi-view tracking, and reprojection filtering support heterogeneous images acquired by different sensors, dates, viewing angles, and resolutions (Qin et al., 2019).

A unified satellite framework couples RPC bias refinement with photometric least-squares matching. Its global energy combines geometric reprojection and local intensity consistency. Direct joint optimization is degenerate because global RPC biases and local feature translations can explain the same displacement. The method alternates between a BA-like RPC-bias update and a geometry-guided matching subproblem, using reprojection-dependent weights and virtual ground-control-point constraints (Ling et al., 2021).

6. Scaling, implementation, and emerging optimization strategies

BA performance is constrained by both nonlinear convergence and linear-system cost. Hardware and software systems therefore exploit different parts of the same structure.

The $n$01-BA accelerator targets Schur elimination on an embedded FPGA–SoC. Its co-observation optimization records the camera set observing each point and computes only the camera-pair contributions induced by that set. If $n$02 is the number of observing cameras for point $n$03, the work is characterized by

$n$04

rather than the worst-case $n$05. A specialized pipeline accelerates point-block inversion, Jacobian products, and Schur updates. On five BAL datasets, the Schur-specific design was approximately $n$06 faster than its single-processing-element version and approximately $n$07 faster than the ARM implementation; end-to-end $n$08-BA averaged $n$09 s per iteration at $n$10 W, compared with $n$11 s at $n$12 W for ARM-only execution (Qin et al., 2019).

Quantum LM replaces the classical linear solve inside each LM iteration with an HHL-based quantum linear-system solver. The method remains a hybrid algorithm: residual and Jacobian construction, damping, and nonlinear updates remain classical, while HHL solves the damped normal equations. The reported experiments involve nine toy problems with ten points and two cameras. Theoretical complexity claims depend on matrix access, condition number, state preparation, Hamiltonian simulation, precision, and measurement, while the demonstrated circuits are small, padded, and simulated. The work therefore establishes a quantum formulation rather than a practical large-scale speedup (Bernecker et al., 2022).

Learned convergence methods modify the LM damping schedule rather than the BA objective. A Soft Actor-Critic policy observes the five most recent reprojection errors and directly selects the damping factor $n$13. On the reported KITTI and BAL experiments, the learned policy reduced BA iterations by approximately $n$14–$n$15 and produced runtime reductions while preserving final accuracy. Its benefit is scale-dependent: for small problems, neural-policy overhead can exceed the saved optimization time (Belder et al., 2023).

An eager-mode GPU implementation integrates BA with PyTorch and PyPose. It constructs sparse block-Jacobian matrices through batched autodifferentiation, stores them in sparse BSR format, and solves damped normal equations using sparse Cholesky or PCG. On BAL and 1DSfM, it reports comparable reprojection MSE to GTSAM, g2o, and Ceres with substantially lower runtimes on an RTX 4090. The reported comparisons are asymmetric because the proposed implementation uses GPU execution while baselines use CPU systems; the method also does not present an explicit Schur-complement derivation (Zhan et al., 2024).

Across these approaches, the dominant design principles are consistent: preserve camera–structure sparsity, reduce repeated processing of observations, represent global modes or uncertainty when iterative solvers are used, and match the optimization representation to the sensor geometry. Conventional BA remains the most direct formulation when explicit structure refinement and raw observation access are required. Reduced, distributed, photometric, geometric, learned, and hardware-specialized variants exchange some of that directness for lower memory use, lower computational cost, or compatibility with particular sensing and deployment regimes.