1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
//! Implementation of the [Levenberg-Marquardt](https://en.wikipedia.org/wiki/Levenberg%E2%80%93Marquardt_algorithm)
//! optimization algorithm using [nalgebra](https://nalgebra.org).
//!
//! This algorithm tries to solve the least squares optimization problem
//! ```math
//! \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
//!   ```math
//!   \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`](struct.LevenbergMarquardt.html).
//!
//! # 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](https://en.wikipedia.org/wiki/Himmelblau%27s_function)
//! for this example.
//! In this case we have `$n = 2$` and `$m = 2$` with
//!
//! ```math
//!   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.
//! ```
//!
//! ```
//! # use nalgebra::*;
//! # use nalgebra::storage::Owned;
//! # use levenberg_marquardt::{LeastSquaresProblem, LevenbergMarquardt};
//! 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`](fn.differentiate_numerically.html)
//! in a unit test to verify that your Jacobian implementation matches the residuals.
#![no_std]
#![cfg_attr(feature = "RUSTC_IS_NIGHTLY", core_intrinsics)]

extern crate alloc;

mod lm;
mod problem;
mod qr;
mod trust_region;
pub(crate) mod utils;

pub use lm::TerminationReason;
pub use problem::LeastSquaresProblem;

pub use utils::{differentiate_holomorphic_numerically, differentiate_numerically};

cfg_if::cfg_if! {
    if #[cfg(feature="minpack-compat")] {
        pub type LevenbergMarquardt = lm::LevenbergMarquardt<f64>;
        pub type MinimizationReport = lm::MinimizationReport<f64>;
    } else {
        pub use lm::{LevenbergMarquardt, MinimizationReport};
    }
}