UNLSF
Solves a nonlinear least-squares problem using a modified Levenberg-Marquardt algorithm and a finite-difference Jacobian.
Required Arguments
FCN — User-supplied subroutine to evaluate the function that defines the least-squares problem. The usage is
CALL FCN (M, N, X, F), where
M – Length of F. (Input)
N – Length of X. (Input)
X – Vector of length N at which point the function is evaluated. (Input)
X should not be changed by FCN.
F – Vector of length M containing the function values at X. (Output)
FCN must be declared EXTERNAL in the calling program.
M — Number of functions. (Input)
X — Vector of length N containing the approximate solution. (Output)
Optional Arguments
N — Number of variables. N must be less than or equal to M. (Input)
Default: N = SIZE (X,1).
XGUESS — Vector of length N containing the initial guess. (Input)
Default: XGUESS = 0.0.
XSCALE — Vector of length N containing the diagonal scaling matrix for the variables. (Input)
XSCALE is used mainly in scaling the gradient and the distance between two points. By default, the values for XSCALE are set internally. See IPARAM(6) in Comment 4.
Default: XSCALE = 1.0.
FSCALE — Vector of length M containing the diagonal scaling matrix for the functions. (Input)
FSCALE is used mainly in scaling the gradient. In the absence of other information, set all entries to 1.0.
Default: FSCALE = 1.0.
IPARAM — Parameter vector of length 6. (Input/Output)
Set IPARAM(1) to zero for default values of IPARAM and RPARAM. See Comment 4.
Default: IPARAM = 0.
RPARAM — Parameter vector of length 7. (Input/Output)
See Comment 4.
FVEC — Vector of length M containing the residuals at the approximate solution. (Output)
FJAC — M by N matrix containing a finite difference approximate Jacobian at the approximate solution. (Output)
LDFJAC — Leading dimension of FJAC exactly as specified in the dimension statement of the calling program. (Input)
Default: LDFJAC = SIZE (FJAC,1).
FORTRAN 90 Interface
Generic: CALL UNLSF (FCN, M, X [, …])
Specific: The specific interface names are S_UNLSF and D_UNLSF.
FORTRAN 77 Interface
Single: CALL UNLSF (FCN, M, N, XGUESS, XSCALE, FSCALE, IPARAM, RPARAM, X, FVEC, FJAC, LDFJAC)
Double: The double precision name is DUNLSF.
Description
The routine UNLSF is based on the MINPACK routine LMDIF by Moré et al. (1980). It uses a modified Levenberg‑Marquardt method to solve nonlinear least squares problems. The problem is stated as follows:
where m ≥ n, F : Rn→Rm, and fi(x) is the i-th component function of F(x). From a current point, the algorithm uses the trust region approach:
subject to ∥xn ‑ xc∥2 ≤ δc
to get a new point xn, which is computed as
where μc = 0 if δc ≥ ∥(J(xc)T J(xc))-1 J(xc)T F(xc)∥2 and μc > 0 otherwise. F(xc) and J(xc) are the function values and the Jacobian evaluated at the current point xc. This procedure is repeated until the stopping criteria are satisfied. For more details, see Levenberg (1944), Marquardt (1963), or Dennis and Schnabel (1983, Chapter 10).
Since a finite-difference method is used to estimate the Jacobian for some single precision calculations, an inaccurate estimate of the Jacobian may cause the algorithm to terminate at a noncritical point. In such cases, high precision arithmetic is recommended. Also, whenever the exact Jacobian can be easily provided, routine UNLSJ should be used instead.
Comments
1. Workspace may be explicitly provided, if desired, by use of U2LSF/DU2LSF. The reference is:
CALL U2LSF (FCN, M, N, XGUESS, XSCALE, FSCALE, IPARAM, RPARAM, X, FVEC, FJAC, LDFJAC, WK, IWK)
The additional arguments are as follows:
WK — Real work vector of length 9 * N + 3 * M ‑ 1. WK contains the following information on output: The second N locations contain the last step taken. The third N locations contain the last Gauss-Newton step. The fourth N locations contain an estimate of the gradient at the solution.
IWK — Integer work vector of length N containing the permutations used in the QR factorization of the Jacobian at the solution.
2. Informational errors
|
Type |
Code |
Description |
|
3 |
1 |
Both the actual and predicted relative reductions in the function are less than or equal to the relative function convergence tolerance. |
|
3 |
2 |
The iterates appear to be converging to a noncritical point. |
|
4 |
3 |
Maximum number of iterations exceeded. |
|
4 |
4 |
Maximum number of function evaluations exceeded. |
|
3 |
6 |
Five consecutive steps have been taken with the maximum step length. |
|
2 |
7 |
Scaled step tolerance satisfied; the current point may be an approximate local solution, or the algorithm is making very slow progress and is not near a solution, or STEPTL is too big. |
3. The first stopping criterion for UNLSF occurs when the norm of the function is less than the absolute function tolerance (RPARAM(4)). The second stopping criterion occurs when the norm of the scaled gradient is less than the given gradient tolerance (RPARAM(1)). The third stopping criterion for UNLSF occurs when the scaled distance between the last two steps is less than the step tolerance (RPARAM(2)).
4. If the default parameters are desired for UNLSF, then set IPARAM(1) to zero and call the routine UNLSF. Otherwise, if any nondefault parameters are desired for IPARAM or RPARAM, then the following steps should be taken before calling UNLSF:
CALL U4LSF (IPARAM, RPARAM)
Set nondefault values for desired IPARAM, RPARAM elements.
Note: The call to U4INF will set IPARAM and RPARAM to their default values so only nondefault values need to be set above.
The following is a list of the parameters and the default values:
IPARAM — Integer vector of length 6.
IPARAM(1) = Initialization flag.
IPARAM(2) = Number of good digits in the function.
Default: Machine dependent.
IPARAM(3) = Maximum number of iterations.
Default: 100.
IPARAM(4) = Maximum number of function evaluations.
Default: 400.
IPARAM(5) = Maximum number of Jacobian evaluations.
Default: Not used in UNLSF.
IPARAM(6) = Internal variable scaling flag.
If IPARAM(6) = 1, then the values for XSCALE are set internally.
Default: 1.
RPARAM — Real vector of length 7.
RPARAM(1) = Scaled gradient tolerance.
The i-th component of the scaled gradient at x is calculated as
where
J(x) is the Jacobian, s = XSCALE, and fs = FSCALE.
Default:
in double where ɛ is the machine precision.
RPARAM(2) = Scaled step tolerance. (STEPTL)
The i-th component of the scaled step between two points x and y is computed as
where s = XSCALE.
Default: ɛ2/3 where ɛ is the machine precision.
RPARAM(3) = Relative function tolerance.
Default: max(10-10, ɛ2/3), max (10-20, ɛ2/3) in double where ɛ is the machine precision.
RPARAM(4) = Absolute function tolerance.
Default: max (10-20, ɛ2), max(10-40, ɛ2) in double where ɛ is the machine precision.
RPARAM(5) = False convergence tolerance.
Default: 100ɛ where ɛ is the machine precision.
RPARAM(6) = Maximum allowable step size.
Default: 1000 max(ɛ1, ɛ2) where
ɛ2 = ∥s∥2, s = XSCALE, and t = XGUESS.
RPARAM(7) = Size of initial trust region radius.
Default: based on the initial scaled Cauchy step.
If double precision is desired, then DU4LSF is called and RPARAM is declared double precision.
5. Users wishing to override the default print/stop attributes associated with error messages issued by this routine are referred to “Error Handling” in the Introduction.
Example
The nonlinear least squares problem
where
is solved. RPARAM(4) is changed to a non-default value.
USE UMACH_INT
USE U4LSF_INT
USE UNLSF_INT
IMPLICIT NONE
! Declaration of variables
INTEGER LDFJAC, M, N
PARAMETER (LDFJAC=2, M=2, N=2)
!
INTEGER IPARAM(6), NOUT
REAL FVEC(M), RPARAM(7),X(N), XGUESS(N)
EXTERNAL ROSBCK
! Compute the least squares for the
! Rosenbrock function.
DATA XGUESS/-1.2E0, 1.0E0/
!
! Relax the first stopping criterion by
! calling U4LSF and scaling the
! absolute function tolerance by 10.
CALL U4LSF (IPARAM, RPARAM)
RPARAM(4) = 10.0E0*RPARAM(4)
!
CALL UNLSF (ROSBCK, M, X,xguess=xguess, iparam=iparam, rparam=rparam,&
fvec=fvec)
! Print results
CALL UMACH (2, NOUT)
WRITE (NOUT,99999) X, FVEC, IPARAM(3), IPARAM(4)
!
99999 FORMAT (' The solution is ', 2F9.4, //, ' The function ', &
'evaluated at the solution is ', /, 18X, 2F9.4, //, &
' The number of iterations is ', 10X, I3, /, ' The ', &
'number of function evaluations is ', I3, /)
END
!
SUBROUTINE ROSBCK (M, N, X, F)
INTEGER M, N
REAL X(N), F(M)
!
F(1) = 10.0E0*(X(2)-X(1)*X(1))
F(2) = 1.0E0 - X(1)
RETURN
END
The solution is 1.0000 1.0000
The function evaluated at the solution is
0.0000 0.0000
The number of iterations is 24
The number of function evaluations is 33
