Crate levenberg_marquardt_sparse 

Sparse implementation of the Levenberg-Marquardt optimization algorithm using nalgebra.

This crate is a fork of levenberg-marquardt that removes the dense QR-based solver path and always uses a sparse normal-equations solve based on Jacobi-preconditioned Conjugate Gradient (PCG).

The port from the original dense implementation to sparse matrices was carried out via agentic coding using GitHub Copilot.

It is intended for problems where the Jacobian is sparse, such as bundle adjustment and SLAM, while still providing a dense-to-sparse bridge through SparseJacobian::from_dense for small examples and migration.

This algorithm tries to solve the least squares optimization problem

\min_{\vec{x}\in\R^n}f(\vec{x})\quad\text{where}\quad\begin{cases}\begin{aligned}
  \ f\!:\R^n &\to \R \\
 \vec{x} &\mapsto \frac{1}{2}\sum_{i=1}^m \bigl(r_i(\vec{x})\bigr)^2,
\end{aligned}\end{cases}

for differentiable residual functions $r_i\!:\R^n\to\R$.

Inputs

The problem has $n$ parameters $\vec{x}\in\R^n$ and $m$ residual functions $r_i\!:\R^n\to\R$.

You must provide an implementation of

• the residual vector $\vec{x} \mapsto (r_1(\vec{x}), \ldots, r_m(\vec{x}))^\top\in\R^m$
• and its Jacobian $\mathbf{J} \in \R^{m\times n}$, returned as a SparseJacobian in triplet (COO) format, defined as

  \mathbf{J} \coloneqq
  \def\arraystretch{1.5}
  \begin{pmatrix}
  \frac{\partial r_1}{\partial x_1} & \cdots & \frac{\partial r_1}{\partial x_n} \\
  \frac{\partial r_2}{\partial x_1} & \cdots & \frac{\partial r_2}{\partial x_n} \\
  \vdots & \ddots & \vdots \\
  \frac{\partial r_m}{\partial x_1} & \cdots & \frac{\partial r_m}{\partial x_n}
  \end{pmatrix}.

Finally, you have to provide an initial guess for $\vec{x}$. This can be a constant value, but typically the optimization result crucially depends on a good initial value.

SparseJacobian stores (row, col, value) entries. You can construct it directly with SparseJacobian::new or SparseJacobian::from_triplets, or use SparseJacobian::from_dense as a bridge when migrating existing dense code.

The algorithm also has a number of hyperparameters which are documented at LevenbergMarquardt.

This crate is #![no_std] and requires alloc.

Usage Example

We use $f(x, y) \coloneqq \frac{1}{2}[(x^2 + y - 11)^2 + (x + y^2 - 7)^2]$ as a test function for this example. In this case we have $n = 2$ and $m = 2$ with

  r_1(\vec{x}) \coloneqq x_1^2 + x_2 - 11\quad\text{and}\quad
  r_2(\vec{x}) \coloneqq x_1 + x_2^2 - 7.

struct ExampleProblem {
    // holds current value of the n parameters
    p: Vector2<f64>,
}

// We implement a trait for every problem we want to solve
impl LeastSquaresProblem<f64, U2, U2> for ExampleProblem {
    type ParameterStorage = Owned<f64, U2>;
    type ResidualStorage = Owned<f64, U2>;

    fn set_params(&mut self, p: &Vector2<f64>) {
        self.p.copy_from(p);
        // do common calculations for residuals and the Jacobian here
    }

    fn params(&self) -> Vector2<f64> { self.p }

    fn residuals(&self) -> Option<Vector2<f64>> {
        let [x, y] = [self.p.x, self.p.y];
        // vector containing residuals $r_1(\vec{x})$ and $r_2(\vec{x})$
        Some(Vector2::new(x*x + y - 11., x + y*y - 7.))
    }

    fn jacobian(&self) -> Option<SparseJacobian<f64>> {
        let [x, y] = [self.p.x, self.p.y];

        // first row of Jacobian, derivatives of first residual
        let d1_x = 2. * x; // $\frac{\partial}{\partial x_1}r_1(\vec{x}) = \frac{\partial}{\partial x} (x^2 + y - 11) = 2x$
        let d1_y = 1.;     // $\frac{\partial}{\partial x_2}r_1(\vec{x}) = \frac{\partial}{\partial y} (x^2 + y - 11) = 1$

        // second row of Jacobian, derivatives of second residual
        let d2_x = 1.;     // $\frac{\partial}{\partial x_1}r_2(\vec{x}) = \frac{\partial}{\partial x} (x + y^2 - 7) = 1$
        let d2_y = 2. * y; // $\frac{\partial}{\partial x_2}r_2(\vec{x}) = \frac{\partial}{\partial y} (x + y^2 - 7) = 2y$

        Some(SparseJacobian::from_triplets(
            2,
            2,
            vec![
                (0, 0, d1_x),
                (0, 1, d1_y),
                (1, 0, d2_x),
                (1, 1, d2_y),
            ],
        ))
    }
}

let problem = ExampleProblem {
    // the initial guess for $\vec{x}$
    p: Vector2::new(1., 1.),
};
let (_result, report) = LevenbergMarquardt::new().minimize(problem);
assert!(report.termination.was_successful());
assert!(report.objective_function.abs() < 1e-10);

Derivative checking

You should try using differentiate_numerically in a unit test to verify that your Jacobian implementation matches the residuals.

Structs

LevenbergMarquardt

Levenberg-Marquardt optimization algorithm.

MinimizationReport

Information about the minimization.

SparseJacobian

Sparse Jacobian in triplet (COO) format.

Enums

TerminationReason

Reasons for terminating the minimization.

Traits

LeastSquaresProblem

A least squares minimization problem.

Functions

differentiate_holomorphic_numerically

Compute a numerical approximation of the Jacobian for holomorphic residuals.

differentiate_numerically

Compute a numerical approximation of the Jacobian.