◆ Lmstr1()

Sub Lmstr1	(	F As	LongPtr,
		M As	Long,
		N As	Long,
		X() As	Double,
		Fvec() As	Double,
		Fjac() As	Double,
		Tol As	Double,
		Ipvt() As	Long,
		Info As	Long,
		Optional Info2 As	Long
	)

Nonlinear least squares approximation by Levenberg-Marquardt method (limited storage) (simple driver)

Purpose: This routine minimizes the sum of the squares of M nonlinear functions in N variables by a modification of the Levenberg-Marquardt algorithm.
minimize the sum of fi(x1, x2, ..., xn)^2 (sum for i = 1 to M)

The user must provide a subroutine which calculates the function values and the Jacobian (row-by-row).

Lmstr1 is equivalent to using Lmstr with setting FTol = Tol, XTol = Tol, GTol = 0, Maxfev = 100*(n+1), Mode = 1, Factor = 100 and Nprint = 0.

Parameters

[in]	F	User supplied subroutine which calculates the functions fi(x) defined as follows. Sub F(M As Long, N As Long, X() As Double, Fvec() As Double, Fjrow() As Double, IFlag As Long) If IFlag = 1: Calculate the function values fi(x) at X() and return in Fvec(i-1) (i = 1 to M). Other variables should not be changed. If IFlag = i (2 <= i <= M+1): Calculate the (i-1)th row of the Jacobian (df(i-1)/dxj) and return in Fjrow(j-1) (j = 1 to N). Other variables should not be changed. End Sub The value of IFlag should not be changed unless the user wants to terminate the execution. In this case, set IFlag to a negative integer.
[in]	M	Number of functions. (M > 0)
[in]	N	Number of variables. (0 < N <= M)
[in,out]	X()	Array X(LX - 1) (LX >= N) [in] An initial estimate of the solution vector. [out] The obtained solution vector.
[out]	Fvec()	Array Fvec(LFvec - 1) (LFvec >= M) The function values evaluated at the solution vector X().
[out]	Fjac()	Array Fjac(LFjac1 - 1, LFjac2 - 1) (LFjac1 >= M, LFjac2 >= N) The upper N x N submatrix contains an upper triangular matrix R with diagonal elements of nonincreasing magnitude such that P^T * (J^T * J)P = R^T R where P is a permutation matrix and J is the final calculated Jacobian.
[in]	Tol	Relative error desired in the sum of squares and the approximate solution. (Tol >= 0)
[out]	Ipvt()	Array Ipvt(LIpvt - 1) (LIpvt >= N) The pivot indices that define the permutation matrix P.
[out]	Info	= 0: Successful exit. (Sub-code is set to Info2) = -2: The argument M had an illegal value. (M < N) = -3: The argument N had an illegal value. (N <= 0) = -4: The argument X() is invalid. (Array X() is not big enough) = -5: The argument Fvec() is invalid. (Array Fvec() is not big enough) = -6: The argument Fjac() is invalid. (Array Fjac() is not big enough) = -7: The argument Tol had an illegal value. (Tol < 0) = -8: The argument Ipiv() is invalid. (Array Ipiv() is not big enough) = 1: Number of calls to F with IFlag = 1 or 2 has reached Maxfev. = 2: FTol is too small. No further reduction in the sum of squares is possible. = 3: XTol is too small. No further improvement in the approximate solution X is possible. = 4: GTol is too small. Fvec is orthogonal to the columns of the Jacobian to machine precision. = 5: User imposed termination (returned from F with IFlag < 0).
[out]	Info2	(Optional) Sub-code on return with Info = 0. = 1: Both actual and predicted relative reductions in the sum of squares are at most FTol. = 2: Relative error between two consecutive iterates is at most XTol. = 3: Both of above are satisfied.

Reference: netlib/minpack

Example Program: Approximate the following data with model function f(x) = c1*(1 - exp(-c2*x)). Two parameters c1 and c2 are determined by the nonlinear least squares method.
f(x) x

10.07 77.6

29.61 239.9

50.76 434.8

81.78 760.0

The initial estimates of the solution are c1 = 500 and c2 = 0.0001.
Sub FLmstr(M As Long, N As Long, X() As Double, Fvec() As Double, Fjrow() As Double, IFlag As Long)

Dim Xdata(3) As Double, Ydata(3) As Double, I As Long

Ydata(0) = 10.07: Xdata(0) = 77.6

Ydata(1) = 29.61: Xdata(1) = 239.9

Ydata(2) = 50.76: Xdata(2) = 434.8

Ydata(3) = 81.78: Xdata(3) = 760

If IFlag = 1 Then

For I = 0 To M - 1

Fvec(I) = Ydata(I) - X(0) * (1 - Exp(-Xdata(I) * X(1)))

Next

ElseIf IFlag >= 2 And IFlag <= M + 1 Then

Fjrow(0) = Exp(-Xdata(IFlag - 2) * X(1)) - 1

Fjrow(1) = -Xdata(IFlag - 2) * X(0) * Exp(-X(1) * Xdata(IFlag - 2))

End If

End Sub

Sub Ex_Lmstr1()

Const M = 4, N = 2

Dim X(N - 1) As Double, Fvec(M - 1) As Double, Fjac(M - 1, N - 1) As Double

Dim Tol As Double, Ipvt(N - 1) As Long, Info As Long

Tol = 0.00000001 '1.0e-8

X(0) = 500: X(1) = 0.0001

Call Lmstr1(AddressOf FLmstr, M, N, X(), Fvec(), Fjac(), Tol, Ipvt(), Info)

Debug.Print "C1, C2 =", X(0), X(1)

Debug.Print "Info =", Info

End Sub

Example Results: C1, C2 = 241.084896089973 5.44942234119956E-04

Info = 0