Crate levenberg_marquardt

Implementation of the Levenberg-Marquardt optimization algorithm using nalgebra.

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}$, 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.

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

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>;
    type JacobianStorage = Owned<f64, U2, U2>;
     
    fn set_params(&mut self, p: &VectorN<f64, U2>) {
        self.p.copy_from(p);
        // do common calculations for residuals and the Jacobian here
    }
     
    fn params(&self) -> VectorN<f64, U2> { 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<Matrix2<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(Matrix2::new(
            d1_x, d1_y,
            d2_x, 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.

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.