[−][src]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.