Least squares fitting (linear/nonlinear)

Least squares method is usually mentioned in two contexts. First, its application is widely known in regression analysis. It is used to construct models on the basis of noisy experimental data. At that, besides constructing a model, it usually performed an estimate of the error with which parameters were calculated. Some other problems could be solved as well. Secondly, the least-squares method is often applied as an fitting method without any connection to statistics. Here this technique is considered as an approximation method. It should also be noted that lsfit unit solves problems of general form. There are units for polynomial, cubic spline and rational interpolation/fitting which provide specialized subroutines with similar interface.

Contents     Linear least squares fitting     Nonlinear least squares fitting: Hessian-based or Hessian-free     Nonlinear least squares fitting: reverse communication interface     Manual entries

Linear least squares fitting

Linear least squares problem can be written in matrix form as

Constrained problem can be written as

Here c denotes coefficients vector. The columns of F correspond to basis functions, rows correspond to points, F_ij  contains the value of a j-th basis function at the i-th point. y vector contains values of function being fitted. W matrix is a diagonal weight matrix. Additional matrix C defines optional constraints. In this formulation, the problem reduces to solving a set of linear equations. QR-based solver is used to solve this system. At first, A is represented by a product of orthogonal matrix Q and upper triangular matrix R. Then, a set of equations Rx = Q ^Tb is solved. If R is degenerate, algorithm switches to SVD-based solver, which is a bit slower, but more stable. Anyway, problem has O(N·M ²) complexity.

lsfit unit includes four subroutines for linear fitting: LSFitLinear (the most simple task - no constraints, identity W), LSFitLinearW (unconstrained weighted fitting), LSFitLinearC (constrained fitting without weights) and LSFitLinearWC (weighted constrained fitting).

Nonlinear least squares fitting: Hessian-based or Hessian-free

Nonlinear least squares problem is much more complex than linear one. A general form function is fitted to a data, a function, which depends on M arguments x_i  and K parameters c_j :

We need to find c that minimize difference between y_i  and f(x_i ,c):

lsfit unit uses Levenberg-Marquardt algorithm (which is implemented in the minlm unit) to solve this problem. As with minlm unit, user may choose between two optimization types: FG (using function value and its gradient) and FGH (using function value, gradient and Hessian). Points may have individual weights (W suffix) or have no individual weights. So we have four optimization subroutines: LSFitNonlinearWFG, LSFitNonlinearFG, LSFitNonlinearWFGH, LSFitNonlinearFGH.

There are two situations possible: you have "expensive" gradient whose complexity is O((M+K) ²) and you have "cheap" gradient whose complexity is significantly lower than O((M+K) ²). "Expensive" gradient is usually calculated using finite difference approximation or analytic formula. Functions with regular structure sometimes allow to calculate "cheap" gradient. If you have "cheap" gradient you may use hybrid Levenberg-Marquardt optimizer with quadratic model reuse. Hybrid algorithm allows to solve optimization problems much more faster than original Levenberg-Marquardt algorithm. It is always used to solve problems with known Hessian (FGH setting) because gradient is always cheap when compared with Hessian. It is optional in FG setting and may be turned on/off with CheapFG parameter of LSFitNonlinearWFG and LSFitNonlinearFG subroutines.

Note #1
Don't calculate the function gradient on the basis of a two-point difference formula because it is insufficiently precise. Use at least a four-point formula or analytical form of the gradient. But the most preferable is to use analytic gradient.

Nonlinear least squares fitting: reverse communication interface

The fitting algorithm shall obtain values of a function/gradient/... during its operation. This problem is solved in most program packages by transferring the pointer to the function (C++, Delphi) or delegate (C#) which is used to calculate function value/gradient/Hessian.

The ALGLIB package, differently from other libraries, makes use of reverse communication to solve this problem. When a value (or derivatives) of a function needs/need to be calculated, the algorithm state is stored within a special structure, control is returned to the calling program, which makes all calculations and recalls the computing subroutine.

Thus, the optimization algorithm is operated in accordance in the following order:

1. LSFitState data structure preparation by means of one of the algorithm initialization subrotuines (LSFitNonlinearWFG, LSFitNonlinearFG, LSFitNonlinearWFGH, LSFitNonlinearFGH).
2. LSFitNonlinearIteration subroutine call.
3. If False is returned from the subroutine, then the algorithm operation is completed, and the minimum is found (the minimum itself can be obtained by calling the LSFitNonlinearResults subprogram).
4. If True is returned from the subroutine, the latter will make a request for information on the function. The function/gradient/Hessian shall be calculated according to the field of the LSFitState structure, which is set as True (this issue is fully detailed below).
5. The LSFitNonlinearIteration subroutine needs to be called again after the requested information is loaded into the LSFitState structure.

The following fields of the LSFitState structure are used in order to exchange information with the user:

• LSFitState.X[0..M-1] – an array storing information on coordinates of the point x
• LSFitState.C[0..K-1] - an array storing function parameters
• LSFitState.F – the value of function F(x) should be stored in this field (if it is requested)
• LSFitState.G[0..K-1] – the gradient grad F(x) should be stored in this field (if it is requested)
• LSFitState.H[0..K-1,0..K-1] – the Hessian matrix H(F) should be stored in this field (if it is requested)

The LSFitNonlinearIteration subroutine can set one, and only one, of the following fields to True, subject to what exactly needs to be calculated:

• NeedF – the value of function F(x) needs to be calculated (Not to be confused with vector f_i )
• NeedFG – the value of function F(x) and gradient grad F(x) need to be calculated
• NeedFGH – the value of function F(x), gradient grad F(x) and Hessian H(F)}} need to be calculated

The table given below will show which of the fields can be set by different algorithm variants.

       NeedF   NeedFG  NeedFGH
FG     X       X       
FGH    X       X       X

Manual entries

C++		lsfit.h
C#		lsfit.cs
MPFR		lsfit.h
Delphi		lsfit.pas
FreePascal		lsfit.pas
VBA		lsfit.bas

This article is intended for personal use only.

Download ALGLIB