# nls.lm: Addresses NLS problems with the Levenberg-Marquardt algorithm In minpack.lm: R Interface to the Levenberg-Marquardt Nonlinear Least-Squares Algorithm Found in MINPACK, Plus Support for Bounds

## Description

The purpose of `nls.lm` is to minimize the sum square of the vector returned by the function `fn`, by a modification of the Levenberg-Marquardt algorithm. The user may also provide a function `jac` which calculates the Jacobian.

## Usage

 ```1 2``` ```nls.lm(par, lower=NULL, upper=NULL, fn, jac = NULL, control = nls.lm.control(), ...) ```

## Arguments

 `par` A list or numeric vector of starting estimates. If `par` is a list, then each element must be of length 1. `lower` A numeric vector of lower bounds on each parameter. If not given, the default lower bound for each parameter is set to `-Inf`. `upper` A numeric vector of upper bounds on each parameter. If not given, the default upper bound for each parameter is set to `Inf`. `fn` A function that returns a vector of residuals, the sum square of which is to be minimized. The first argument of `fn` must be `par`. `jac` A function to return the Jacobian for the `fn` function. `control` An optional list of control settings. See `nls.lm.control` for the names of the settable control values and their effect. `...` Further arguments to be passed to `fn` and `jac`.

## Details

Both functions `fn` and `jac` (if provided) must return numeric vectors. Length of the vector returned by `fn` must not be lower than the length of `par`. The vector returned by `jac` must have length equal to length(\code{fn}(\code{par}, …)) * length(\code{par}).

The `control` argument is a list; see `nls.lm.control` for details.

Successful completion.

The accuracy of `nls.lm` is controlled by the convergence parameters `ftol`, `ptol`, and `gtol`. These parameters are used in tests which make three types of comparisons between the approximation par and a solution par0. `nls.lm` terminates when any of the tests is satisfied. If any of the convergence parameters is less than the machine precision, then `nls.lm` only attempts to satisfy the test defined by the machine precision. Further progress is not usually possible.
The tests assume that `fn` as well as `jac` are reasonably well behaved. If this condition is not satisfied, then `nls.lm` may incorrectly indicate convergence. The validity of the answer can be checked, for example, by rerunning `nls.lm` with tighter tolerances.

First convergence test.
If |z| denotes the Euclidean norm of a vector z, then this test attempts to guarantee that

|fvec| < (1 + \code{ftol})|fvec0|,

where fvec0 denotes the result of `fn` function evaluated at par0. If this condition is satisfied with `ftol` ~ 10^(-k), then the final residual norm |fvec| has k significant decimal digits and `info` is set to 1 (or to 3 if the second test is also satisfied). Unless high precision solutions are required, the recommended value for `ftol` is the square root of the machine precision.

Second convergence test.
If D is the diagonal matrix whose entries are defined by the array `diag`, then this test attempt to guarantee that

|D*(par - par0)| < \code{ptol}|D*par0|,

If this condition is satisfied with `ptol` ~ 10^(-k), then the larger components of D*par have k significant decimal digits and `info` is set to 2 (or to 3 if the first test is also satisfied). There is a danger that the smaller components of D*par may have large relative errors, but if `diag` is internally set, then the accuracy of the components of par is usually related to their sensitivity. Unless high precision solutions are required, the recommended value for `ptol` is the square root of the machine precision.

Third convergence test.
This test is satisfied when the cosine of the angle between the result of `fn` evaluation fvec and any column of the Jacobian at par is at most `gtol` in absolute value. There is no clear relationship between this test and the accuracy of `nls.lm`, and furthermore, the test is equally well satisfied at other critical points, namely maximizers and saddle points. Therefore, termination caused by this test (`info` = 4) should be examined carefully. The recommended value for `gtol` is zero.

Unsuccessful completion.

Unsuccessful termination of `nls.lm` can be due to improper input parameters, arithmetic interrupts, an excessive number of function evaluations, or an excessive number of iterations.

Improper input parameters.
`info` is set to 0 if length(\code{par}) = 0, or length(fvec) < length(\code{par}), or `ftol` < 0, or `ptol` < 0, or `gtol` < 0, or `maxfev` <= 0, or `factor` <= 0.

Arithmetic interrupts.
If these interrupts occur in the `fn` function during an early stage of the computation, they may be caused by an unacceptable choice of par by `nls.lm`. In this case, it may be possible to remedy the situation by rerunning `nls.lm` with a smaller value of `factor`.

Excessive number of function evaluations.
A reasonable value for `maxfev` is 100*(length(\code{par}) + 1). If the number of calls to `fn` reaches `maxfev`, then this indicates that the routine is converging very slowly as measured by the progress of fvec and `info` is set to 5. In this case, it may be helpful to force `diag` to be internally set.

Excessive number of function iterations.
The allowed number of iterations defaults to 50, can be increased if desired.

The list returned by `nls.lm` has methods for the generic functions `coef`, `deviance`, `df.residual`, `print`, `residuals`, `summary`, `confint`, and `vcov`.

## Value

A list with components:

 `par` The best set of parameters found. `hessian` A symmetric matrix giving an estimate of the Hessian at the solution found. `fvec` The result of the last `fn` evaluation; that is, the residuals. `info` `info` is an integer code indicating the reason for termination. 0Improper input parameters. 1Both actual and predicted relative reductions in the sum of squares are at most `ftol`. 2Relative error between two consecutive iterates is at most `ptol`. 3Conditions for `info` = 1 and `info` = 2 both hold. 4The cosine of the angle between `fvec` and any column of the Jacobian is at most `gtol` in absolute value. 5Number of calls to `fn` has reached `maxfev`. 6`ftol` is too small. No further reduction in the sum of squares is possible. 7`ptol` is too small. No further improvement in the approximate solution `par` is possible. 8`gtol` is too small. `fvec` is orthogonal to the columns of the Jacobian to machine precision. 9The number of iterations has reached `maxiter`. `message` character string indicating reason for termination

.

 `diag` The result list of `diag`. See Details. `niter` The number of iterations completed before termination. `rsstrace` The residual sum of squares at each iteration. Can be used to check the progress each iteration. `deviance` The sum of the squared residual vector.

## Note

The public domain FORTRAN sources of MINPACK package by J.J. Moré, implementing the Levenberg-Marquardt algorithm were downloaded from http://netlib.org/minpack/, and left unchanged. The contents of this manual page are largely extracted from the comments of MINPACK sources.

## References

J.J. Moré, "The Levenberg-Marquardt algorithm: implementation and theory," in Lecture Notes in Mathematics 630: Numerical Analysis, G.A. Watson (Ed.), Springer-Verlag: Berlin, 1978, pp. 105-116.

`optim`, `nls`, `nls.lm.control`
 ``` 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``` ```###### example 1 ## values over which to simulate data x <- seq(0,5,length=100) ## model based on a list of parameters getPred <- function(parS, xx) parS\$a * exp(xx * parS\$b) + parS\$c ## parameter values used to simulate data pp <- list(a=9,b=-1, c=6) ## simulated data, with noise simDNoisy <- getPred(pp,x) + rnorm(length(x),sd=.1) ## plot data plot(x,simDNoisy, main="data") ## residual function residFun <- function(p, observed, xx) observed - getPred(p,xx) ## starting values for parameters parStart <- list(a=3,b=-.001, c=1) ## perform fit nls.out <- nls.lm(par=parStart, fn = residFun, observed = simDNoisy, xx = x, control = nls.lm.control(nprint=1)) ## plot model evaluated at final parameter estimates lines(x,getPred(as.list(coef(nls.out)), x), col=2, lwd=2) ## summary information on parameter estimates summary(nls.out) ###### example 2 ## function to simulate data f <- function(TT, tau, N0, a, f0) { expr <- expression(N0*exp(-TT/tau)*(1 + a*cos(f0*TT))) eval(expr) } ## helper function for an analytical gradient j <- function(TT, tau, N0, a, f0) { expr <- expression(N0*exp(-TT/tau)*(1 + a*cos(f0*TT))) c(eval(D(expr, "tau")), eval(D(expr, "N0" )), eval(D(expr, "a" )), eval(D(expr, "f0" ))) } ## values over which to simulate data TT <- seq(0, 8, length=501) ## parameter values underlying simulated data p <- c(tau = 2.2, N0 = 1000, a = 0.25, f0 = 8) ## get data Ndet <- do.call("f", c(list(TT = TT), as.list(p))) ## with noise N <- Ndet + rnorm(length(Ndet), mean=Ndet, sd=.01*max(Ndet)) ## plot the data to fit par(mfrow=c(2,1), mar = c(3,5,2,1)) plot(TT, N, bg = "black", cex = 0.5, main="data") ## define a residual function fcn <- function(p, TT, N, fcall, jcall) (N - do.call("fcall", c(list(TT = TT), as.list(p)))) ## define analytical expression for the gradient fcn.jac <- function(p, TT, N, fcall, jcall) -do.call("jcall", c(list(TT = TT), as.list(p))) ## starting values guess <- c(tau = 2.2, N0 = 1500, a = 0.25, f0 = 10) ## to use an analytical expression for the gradient found in fcn.jac ## uncomment jac = fcn.jac out <- nls.lm(par = guess, fn = fcn, jac = fcn.jac, fcall = f, jcall = j, TT = TT, N = N, control = nls.lm.control(nprint=1)) ## get the fitted values N1 <- do.call("f", c(list(TT = TT), out\$par)) ## add a blue line representing the fitting values to the plot of data lines(TT, N1, col="blue", lwd=2) ## add a plot of the log residual sum of squares as it is made to ## decrease each iteration; note that the RSS at the starting parameter ## values is also stored plot(1:(out\$niter+1), log(out\$rsstrace), type="b", main="log residual sum of squares vs. iteration number", xlab="iteration", ylab="log residual sum of squares", pch=21,bg=2) ## get information regarding standard errors summary(out) ```