Non-Linear Least Squares Overview

Non-linear least squares (NLLS) is the mathematical framework underlying all refinement in calibration-rs. After linear initialization provides approximate parameter estimates, NLLS minimizes the reprojection error to achieve sub-pixel accuracy.

Problem Formulation

Objective: Minimize the sum of squared residuals:

$$\underset{\mathbf{\theta }}{\min }F(\mathbf{\theta })=\frac{1}{2}\sum _{i=1}^N\Vert {\mathbf{r}}_i(\mathbf{\theta }){\Vert }^2=\frac{1}{2}\mathbf{r}(\mathbf{\theta }{)}^T\mathbf{r}(\mathbf{\theta })$$

where $\mathbf{\theta }\in {\mathbb{R}}^n$ is the parameter vector and $\mathbf{r}(\mathbf{\theta })\in {\mathbb{R}}^m$ is the stacked residual vector.

In camera calibration, a typical residual is the reprojection error: the difference between an observed pixel and the predicted projection:

$${\mathbf{r}}_i=\pi (K,\mathbf{d},T_v,{\mathbf{P}}_j)-{\mathbf{p}}_{vj}$$

where $\pi$ is the camera projection function, and the parameters $\mathbf{\theta }$ include intrinsics $K$, distortion $\mathbf{d}$, and poses $\{T_v\}$.

Gauss-Newton Method

The Gauss-Newton method exploits the least-squares structure. At the current estimate $\mathbf{\theta }$, linearize the residuals:

$$\mathbf{r}(\mathbf{\theta }+\mathbf{\delta })\approx \mathbf{r}(\mathbf{\theta })+J\mathbf{\delta }$$

where $J=\frac{\partial \mathbf{r}}{\partial \mathbf{\theta }}\in {\mathbb{R}}^{m\times n}$ is the Jacobian.

Substituting into the objective and minimizing with respect to $\mathbf{\delta }$:

$$\frac{\partial }{\partial \mathbf{\delta }}\frac{1}{2}\Vert \mathbf{r}+J\mathbf{\delta }{\Vert }^2=0$$

gives the normal equations:

$$J^{T}J\mathbf{\delta }=-J^T\mathbf{r}$$

The matrix $H_{\text{GN}}=J^{T}J$ is the Gauss-Newton approximation to the Hessian ($H_{\text{GN}}\approx {\nabla }^{2}F$ when residuals are small).

Update: $\mathbf{\theta }\leftarrow \mathbf{\theta }+\mathbf{\delta }$

Levenberg-Marquardt Method

Gauss-Newton can diverge when the linearization is poor (far from the minimum) or when $J^{T}J$ is singular. Levenberg-Marquardt (LM) adds a damping term:

$$(J^{T}J+\lambda \text{diag}(J^{T}J))\mathbf{\delta }=-J^T\mathbf{r}$$

The damping parameter $\lambda >0$ interpolates between:

• $\lambda \to 0$: Pure Gauss-Newton (fast convergence near minimum)
• $\lambda \to \infty$: Gradient descent with small step (safe far from minimum)

Trust Region Interpretation

LM is a trust region method: $\lambda$ controls the size of the region where the linear approximation is trusted.

• If the update reduces the cost: accept the step, decrease $\lambda$ (expand trust region)
• If the update increases the cost: reject the step, increase $\lambda$ (shrink trust region)

Convergence Criteria

LM terminates when any of:

• Cost threshold: $F(\mathbf{\theta })<{ϵ}_{\text{min}}$
• Absolute decrease: $|F_k-F_{k+1}|<{ϵ}_{\text{abs}}$
• Relative decrease: $|F_k-F_{k+1}|/F_k<{ϵ}_{\text{rel}}$
• Maximum iterations: $k>k_{\max }$
• Parameter change: $\Vert \mathbf{\delta }\Vert <{ϵ}_{\text{param}}$

Sparsity

In bundle adjustment problems, the Jacobian $J$ is sparse: each residual depends on only a few parameter blocks (one camera intrinsics, one distortion, one pose). $J^{T}J$ has a block-arrow structure, so calibration-rs solves the damped normal equations with a sparse Cholesky factorization: the sparsity pattern and its symbolic analysis are computed once per solve, and only the numeric factorization repeats per step. The Marquardt damping keeps the system positive definite, which is what makes Cholesky (rather than QR) sufficient. See Solver Backends.

Cost Function vs. Reprojection Error

The optimizer minimizes the cost $F=\frac{1}{2}{\sum }_i{\rho }_i(\Vert {\mathbf{r}}_i{\Vert }^2)$, which is $\frac{1}{2}\sum \Vert {\mathbf{r}}_i{\Vert }^2$ without a robust loss. The commonly reported mean reprojection error is:

$$\bar{e}=\frac{1}{N}\sum _{i=1}^N\Vert {\mathbf{r}}_i\Vert$$

These are related but not identical: the cost weights large errors quadratically, while the mean error weighs them linearly. calibration-rs reports both: the cost from the solver and the mean reprojection error computed post-optimization.

With Robust Losses

When a robust loss $\rho$ is used, the objective becomes:

$$\underset{\mathbf{\theta }}{\min }\sum _{i=1}^N\rho (\Vert {\mathbf{r}}_i\Vert )$$

The normal equations are modified to incorporate the loss function's weight:

$$J^{T}WJ\mathbf{\delta }=-J^{T}W\mathbf{r}$$

where $W=\text{diag}({\rho }^{\prime }(\Vert {\mathbf{r}}_i\Vert )/\Vert {\mathbf{r}}_i\Vert )$ is the iteratively reweighted diagonal matrix. See Robust Loss Functions for details.