### 10 releases (breaking)

0.7.1 | Oct 16, 2023 |
---|---|

0.6.0 | Apr 30, 2023 |

0.5.0 | Jan 28, 2023 |

0.4.0 | Sep 25, 2022 |

0.2.0 | Mar 16, 2021 |

#**61** in Math

**681** downloads per month

**MIT**license

290KB

4.5K
SLoC

# varpro

Nonlinear function fitting made simple. This library provides robust and fast least-squares fitting of a wide class of model functions to data. It uses the VarPro algorithm to achieve this, hence the name.

## Brief Introduction

This crate implements a powerful algorithm
to fit model functions to data, but it is restricted to so called *separable*
models. See the next section for an explanation. The lack of formulas on this
site makes it hard to get into the depth of the what and how of this crate at this point.
Refer to the documentation for all the meaty details including the math.

### What are Separable Models?

Put simply, separable models are nonlinear functions that can be
written as a *linear combination* of some *nonlinear* basis functions.
A common use case for VarPro is e.g. fitting sums of exponentials,
which is a notoriously ill-conditioned problem.

### What is VarPro?

Variable Projection (VarPro) is an algorithm that takes advantage of the fact that the given fitting problem can be separated into linear and truly nonlinear parameters. The linear parameters are eliminated using some clever linear algebra and the fitting problem is cast into a problem that only depends on the nonlinear parameters. This reduced problem is then solved by using a common nonlinear fitting algorithm, such as Levenberg-Marquardt (LM).

### When Should You Give it a Try?

VarPro can dramatically increase the robustness and speed of the fitting process compared to using a "normal" nonlinear least squares fitting algorithm. When

- the model function you want to fit is a linear combination of nonlinear functions
*and*you know the analytical derivatives of all those functions

*then* you should give it a whirl.

## Example Usage

The following example shows how to use varpro to fit a double exponential decay
with constant offset to a data vector

obtained at grid points `y`

.
Refer to the documentation for a more in-depth guide.`x`

The exponential decay and it's derivative are given as:

`use` `nalgebra``::`DVector`;`
`fn` `exp_decay``(``x` `:``&``DVector``<``f64``>`, `tau` `:` `f64``)`` ``->` `DVector``<``f64``>` `{`
x`.``map``(``|``x``|``(``-`x`/`tau`)``.``exp``(``)``)`
`}`
`fn` `exp_decay_dtau``(``tvec``:` `&``DVector``<``f64``>`,`tau``:` `f64``)`` ``->` `DVector``<``f64``>` `{`
tvec`.``map``(``|``t``|` `(``-`t `/` tau`)``.``exp``(``)` `*` t `/` tau`.``powi``(``2``)``)`
`}`

The steps to perform the fitting are:

`use` `varpro``::``prelude``::``*``;`
`use` `varpro``::``solvers``::``levmar``::``{`LevMarProblemBuilder`,` LevMarSolver`}``;`
`let` x `=` `/*`time or spatial coordinates of the observations`*/``;`
`let` y `=` `/*`the observed data we want to fit`*/``;`
`//` 1. create the model by giving only the nonlinear parameter names it depends on
`let` model `=` `SeparableModelBuilder``::``<``f64``>``::`new`(``&``[``"`tau1`"``,` `"`tau2`"``]``)`
`.``function``(``&``[``"`tau1`"``]``,` exp_decay`)`
`.``partial_deriv``(``"`tau1`"``,` exp_decay_dtau`)`
`.``function``(``&``[``"`tau2`"``]``,` exp_decay`)`
`.``partial_deriv``(``"`tau2`"``,` exp_decay_dtau`)`
`.``invariant_function``(``|``x``|``DVector``::`from_element`(`x`.``len``(``)``,``1.``)``)`
`.``independent_variable``(`x`)`
`.``initial_parameters``(`initial_params`)`
`.``build``(``)`
`.``unwrap``(``)``;`
`//` 2. Cast the fitting problem as a nonlinear least squares minimization problem
`let` problem `=` `LevMarProblemBuilder``::`new`(`model`)`
`.``observations``(`y`)`
`.``build``(``)`
`.``unwrap``(``)``;`
`//` 3. Solve the fitting problem
`let` fit_result `=` `LevMarSolver``::`new`(``)`
`.``fit``(`problem`)`
`.``expect``(``"`fit must exit successfully`"``)``;`
`//` 4. obtain the nonlinear parameters after fitting
`let` alpha `=` fit_result`.``nonlinear_parameters``(``)``;`
`//` 5. obtain the linear parameters
`let` c `=` fit_result`.``linear_coefficients``(``)``.``unwrap``(``)``;`

For more examples please refer to the crate documentation.

### Fit Statistics

Additionally to the

member function, the `fit`

also provides a
`LevMarSolver`

function that calculates some additional statistical
information after the fit has finished.`fit_with_statistics`

### Maximum Performance

The example code above will already run many times faster than just using a nonlinear solver without the magic of varpro. But this crate offers an additional way to eek out the last bits of performance.

The

trait can be manually
implemented to describe a model function. This often allows us to shave of the
last hundreds of microseconds from the computation e.g. by caching intermediate
calculations. The crate documentation contains detailed examples.`SeparableNonlinearModel`

## Acknowledgements

I am grateful to Professor Dianne P. O'Leary and Bert W. Rust ✝ who published the paper that enabled me to understand varpro and come up with this implementation. Professor O'Leary also graciously answered my questions on her paper and some implementation details.

## References and Further Reading

(O'Leary2013) O’Leary, D.P., Rust, B.W. Variable projection for nonlinear least squares problems.
*Comput Optim Appl* **54**, 579–593 (2013). DOI: 10.1007/s10589-012-9492-9

**attention**: the O'Leary paper contains errors that are fixed (so I hope)
in this blog article of mine.

(Golub2003) Golub, G. , Pereyra, V Separable nonlinear least squares:
the variable projection method and its applications. Inverse Problems **19** R1 (2003)
https://iopscience.iop.org/article/10.1088/0266-5611/19/2/201

#### Dependencies

~4MB

~88K SLoC