@@ -11,77 +11,94 @@ Nonlinear function fitting made simple. This library provides robust and fast
1111least-squares fitting of a wide class of model functions to data.
1212It uses the VarPro algorithm to achieve this, hence the name.
1313
14- ## Brief Introduction
14+ ## Introduction
1515
16- This crate implements a powerful algorithm
17- to fit model functions to data, but it is restricted to so called _ separable_
18- models. See the next section for an explanation. The lack of formulas on this
19- site makes it hard to get into the depth of the what and how of this crate at this point.
20- [ Refer to the documentation ] ( https://docs.rs/varpro/ ) for all the meaty details including the math.
16+ This crate implements a powerful algorithm to fit model functions to data,
17+ but it is restricted to so called _ separable_ models. The lack of formulas on
18+ this site makes it hard to go into detail, but a brief overview is provided in
19+ the next sections. [ Refer to the documentation ] ( https://docs.rs/varpro/ ) for all
20+ the meaty details including the math.
2121
2222### What are Separable Models?
2323
24- Put simply, separable models are nonlinear functions that can be
24+ Put simply, separable models are nonlinear functions which can be
2525written as a * linear combination* of some * nonlinear* basis functions.
2626A common use case for VarPro is e.g. fitting sums of exponentials,
2727which is a notoriously ill-conditioned problem.
2828
2929### What is VarPro?
3030
31- Variable Projection (VarPro) is an algorithm that takes advantage of the fact
32- that the given fitting problem can be separated into linear and truly nonlinear parameters.
33- The linear parameters are eliminated using some clever linear algebra
34- and the fitting problem is cast into a problem that only depends on the nonlinear parameters.
35- This reduced problem is then solved by using a common nonlinear fitting algorithm,
31+ Variable Projection (VarPro) is an algorithm that exploits that its fitting
32+ problem can be separated into linear and nonlinear parameters.
33+ First, the linear parameters are eliminated using some clever linear algebra. Then,
34+ the fitting problem is rewritten so that it depends only on the nonlinear parameters.
35+ Finally, this reduced problem is solved by using a general purpose nonlinear minimization algorithm,
3636such as Levenberg-Marquardt (LM).
3737
3838### When Should You Give it a Try?
3939
4040VarPro can dramatically increase the robustness and speed of the fitting process
41- compared to using a "normal" nonlinear least squares fitting algorithm. When
41+ compared to using a general purpose nonlinear least squares fitting algorithm. When
4242
43- * the model function you want to fit is a linear combination of nonlinear functions
43+ * the model function you want to fit is a linear combination of nonlinear functions,
4444* _ and_ you know the analytical derivatives of all those functions
4545
46- _ then_ you should give it a whirl.
46+ _ then_ you should give it a whirl. Also consider the section on global fitting below,
47+ which provides another great use case for this crate.
4748
4849## Example Usage
4950
50- The following example shows how to use varpro to fit a double exponential decay
51- with constant offset to a data vector ` y ` obtained at grid points ` x ` .
51+ The following example shows, how to use this crate to fit a double exponential decay
52+ with constant offset to a data vector ` y ` obtained at time points ` t ` .
5253[ Refer to the documentation] ( https://docs.rs/varpro/ ) for a more in-depth guide.
5354
54- The exponential decay and it's derivative are given as:
55-
5655``` rust
57- use nalgebra :: DVector ;
58- fn exp_decay (x : & DVector <f64 >, tau : f64 ) -> DVector <f64 > {
59- x . map (| x | (- x / tau ). exp ())
56+ use varpro :: prelude :: * ;
57+ use varpro :: solvers :: levmar :: {LevMarProblemBuilder , LevMarSolver };
58+ use nalgebra :: {dvector,DVector };
59+
60+ // Define the exponential decay e^(-t/tau).
61+ // Both of the nonlinear basis functions in this example
62+ // are exponential decays.
63+ fn exp_decay (t : & DVector <f64 >, tau : f64 )
64+ -> DVector <f64 > {
65+ t . map (| t | (- t / tau ). exp ())
6066}
6167
62- fn exp_decay_dtau (tvec : & DVector <f64 >,tau : f64 ) -> DVector <f64 > {
63- tvec . map (| t | (- t / tau ). exp () * t / tau . powi (2 ))
68+ // the partial derivative of the exponential
69+ // decay with respect to the nonlinear parameter tau.
70+ // d/dtau e^(-t/tau) = e^(-t/tau)*t/tau^2
71+ fn exp_decay_dtau (t : & DVector <f64 >,tau : f64 )
72+ -> DVector <f64 > {
73+ t . map (| t | (- t / tau )
74+ . exp () * t / tau . powi (2 ))
6475}
65- ```
66-
67- The steps to perform the fitting are:
6876
69- ``` rust
70- use varpro :: prelude :: * ;
71- use varpro :: solvers :: levmar :: {LevMarProblemBuilder , LevMarSolver };
72-
73- let x = /*time or spatial coordinates of the observations*/ ;
74- let y = /*the observed data we want to fit*/ ;
77+ // temporal (or spatial) coordintates of the observations
78+ let t = dvector! [0 . ,1 . ,2 . ,3 . ,4 . ,5 . ,6 . ,7 . ,8 . ,9 . ,10 . ];
79+ // the observations we want to fit
80+ let y = dvector! [6.0 ,4.8 ,4.0 ,3.3 ,2.8 ,2.5 ,2.2 ,1.9 ,1.7 ,1.6 ,1.5 ];
7581
7682// 1. create the model by giving only the nonlinear parameter names it depends on
7783let model = SeparableModelBuilder :: <f64 >:: new (& [" tau1" " tau2"
84+ // provide the nonlinear basis functions and their derivatives.
85+ // In general, base functions can depend on more than just one parameter.
86+ // first function:
7887 . function (& [" tau1" exp_decay )
7988 . partial_deriv (" tau1" exp_decay_dtau )
89+ // second function and derivatives with respect to all parameters
90+ // that it depends on (just one in this case)
8091 . function (& [" tau2" exp_decay )
8192 . partial_deriv (" tau2" exp_decay_dtau )
82- . invariant_function (| x | DVector :: from_element (x . len (),1 . ))
83- . independent_variable (x )
84- . initial_parameters (initial_params )
93+ // a constant offset is added as an invariant basefunction
94+ // as a vector of ones. It is multiplied with its own linear coefficient,
95+ // creating a fittable offset
96+ . invariant_function (| v | DVector :: from_element (v . len (),1 . ))
97+ // give the coordinates of the problem
98+ . independent_variable (t )
99+ // provide guesses only for the nonlinear parameters in the
100+ // order that they were given on construction.
101+ . initial_parameters (vec! [2.5 ,5.5 ])
85102 . build ()
86103 . unwrap ();
87104// 2. Cast the fitting problem as a nonlinear least squares minimization problem
@@ -90,7 +107,7 @@ let problem = LevMarProblemBuilder::new(model)
90107 . build ()
91108 . unwrap ();
92109// 3. Solve the fitting problem
93- let fit_result = LevMarSolver :: new ()
110+ let fit_result = LevMarSolver :: default ()
94111 . fit (problem )
95112 . expect (" fit must exit successfully"
96113// 4. obtain the nonlinear parameters after fitting
@@ -99,37 +116,42 @@ let alpha = fit_result.nonlinear_parameters();
99116let c = fit_result . linear_coefficients (). unwrap ();
100117```
101118
102- For more examples please refer to the crate documentation.
119+ For more in depth examples, please refer to the crate documentation.
103120
104121### Fit Statistics
105122
106- Additionally to the ` fit ` member function, the ` LevMarSolver ` also provides a
107- ` fit_with_statistics ` function that calculates some additional statistical
108- information after the fit has finished .
123+ Additionally to the ` fit ` member function, the ` LevMarSolver ` provides a
124+ ` fit_with_statistics ` function that calculates quite a bit of useful additional statistical
125+ information.
109126
110127### Global Fitting of Multiple Right Hand Sides
111128
112- Before, we have passed a single column vector as the observations. It is also
113- possible to pass a matrix, whose columns represent individual observations. We
114- are now fitting a problem with multiple right hand sides. ` vapro ` will performa a _ global fit_
115- in which the nonlinear parameters are optimized across all right hand sides, but
116- linear parameters of the fit are optimized for each right hand side individually.
129+ In the example above, we have passed a single column vector as the observations.
130+ The library also allows fitting multiple right hand sides, by constructing a
131+ problem via ` LevMarProblemBuilder::mrhs ` . When fitting multiple right hand sides,
132+ ` vapro ` will performa a _ global fit_ , in which the nonlinear parameters are optimized
133+ across all right hand sides, but linear coefficients of the fit are optimized for
134+ each right hand side individually.
117135
118- This is an application where varpro really shines because it can take advantage
119- of the separable nature of the problem. It allows us to perform a global fit over thousands
136+ This is another application where varpro really shines, since it can take advantage
137+ of the separable nature of the problem. It allows us to perform a global fit over thousands,
120138or even tens of thousands of right hand sides in reasonable time (fractions of seconds to minutes),
121139where conventional purely nonlinear solvers must perform much more work.
122140
123- ### Maximum Performance
141+ ### Maximum Performance and Advanced Use Cases
124142
125143The example code above will already run many times faster
126144than just using a nonlinear solver without the magic of varpro.
127- But this crate offers an additional way to eek out the last bits of performance.
145+ But this crate offers an additional way to eek out the last bits of performance.
146+
147+ The ` SeparableNonlinearModel ` trait can be manually implemented to describe a
148+ model function. This often allows us to shave of the last hundreds of microseconds
149+ from the computation, e.g. by caching intermediate calculations. The crate documentation
150+ contains detailed examples.
128151
129- The ` SeparableNonlinearModel ` trait can be manually
130- implemented to describe a model function. This often allows us to shave of the
131- last hundreds of microseconds from the computation e.g. by caching intermediate
132- calculations. The crate documentation contains detailed examples.
152+ This is not only useful for performance, but also for use cases that are difficult
153+ or impossible to accomodate using only the ` SeparableModelBuilder ` . The builder
154+ was created for ease of use _ and_ performance, but it has some limitations by design.
133155
134156## Acknowledgements
135157
0 commit comments