N2fb N2fb_r N2gb N2gb_r N2pb N2pb_r

N2pb_r() Sub N2pb_r ( M As  Long, Md As  Long, N As  Long, X() As  Double, B() As  Double, Info As  Long, M1 As  Long, M2 As  Long, YY() As  Double, YYp() As  Double, IRev As  Long, Optional Iout As  Long = 0, Optional Info2 As  Long, Optional NFcall As  Long, Optional NFjcall As  Long, Optional Niter As  Long, Optional S As  Double, Optional NFGcal As  Long, Optional Rtol As  Double = -1, Optional Atol As  Double = -1, Optional MaxFcall As  Long = -1, Optional MaxIter As  Long = -1, Optional Dtype As  Long = -1, Optional Dfac As  Double = -1, Optional Dtol As  Double = -1, Optional D0 As  Double = -1, Optional Tuner1 As  Double = -1, Optional Xctol As  Double = -1, Optional Xftol As  Double = -1, Optional Lmax0 As  Double = -1, Optional Lmaxs As  Double = -1, Optional Sctol As  Double = -1  ) Nonlinear least squares approximation (adaptive algorithm) (limited storage version) (reverse communication version) Purpose This routine minimizes the sum of the squares of M nonlinear functions in N variables by the adaptive algorithm which combines and augments a Gauss-Newton, Levenberg-Marquardt and other techniques for better convergence. minimize the sum of fi(x1, x2, ..., xn)^2 (sum for i = 1 to M) The user must provide the function and Jacobian values in accordance with IRev. The function and Jacobian can be calculated in chunks rather than all at once. Parameters [in] M Number of data. (M > 0) [in] Md Maximum number of residual components provided by a single evaluation of function values with IRev = 1. (0 < Md <= M) [in] N Number of parameters. (0 < N <= M) [in,out] X() Array X(LX - 1) (LX >= N) [in] X must contain an initial estimate of the solution vector. [out] IRev = 0: Solution vector, i.e. best point so far found.   IRev = 1, 2: The abscissa where the function value or derivatives shoule be evaluated.   IRev = 3: Recent approximation of the solution vector. [in] B() Array B(LB1 - 1, LB2 - 1) (LB1 = 2, LB2 >= N) Lower and upper bounds on solution vector X.   B(0, I) <= X(I) <= B(1, I) (I = 0 to N-1) [out] Info = 0: Successful exit. (Sub-code is set to Info2) = -1: The argument M had an illegal value. (M < N) = -2: The argument Md had an illegal value. (Md < 1 or Md > M) = -3: The argument N had an illegal value. (N < 1) = -4: The argument X() is invalid. = -5: The argument B() is invalid. = -9: The argument YY() is invalid. = -10: The argument YYp() is invalid. = 7: Singular convergence. (Hessian near the current iterate appears to be singular) = 8: False convergence. (Iterate appears to be converging to a noncritical point. Tolerances may be too small) = 9: Function evaluation limit reached. = 10: Iteration limit reached. = 63: F(X) cannot be computed at the initial X. = 65: The gradient could not be computed at X. [out] M1 Index of the first function value or the first row of Jacobian matrix to be provided when returned with IRev = 1 or 2. (1 <= M1 <= M) [out] M2 Index of the last function value or the last row of Jacobian matrix to be provided when returned with IRev = 1 or 2. (M2 = min(M, M1+Md-1)) [in] YY() Array YYp(LYY - 1) (LYY >= Nd) When returned with IRev = 1, the function values with the indices M1+i-1 (i = 1 to M2-M1+1) at X() should be provided in YY(i-1) in the next call. [in] YYp() Array YYp(LYYp1 - 1, LYYp2 - 1) (LYYp1 >= Nd, LYYp2 >= N) When returned with IRev = 2, (M1+i-1)-th rows of Jacobian matrix at X() should be provided in YYp(i-1,j-1) in the next call (i = 1 to M2-M1+1, j = 1 to N). [in,out] IRev Control variable for reverse communication. [in] Before first call, IRev should be initialized to zero. On succeeding calls, IRev should not be altered. [out] If IRev is not zero, set the required values to the specified variables or display the intermediate result as follows and call the routine again. = 0: Computation finished. See return code in Info. = 1: User should provide the function values with the indices M1+i-1 (i = 1 to M2-M1+1) at X() in YY(i-1). Do not alter any variables other than YY(). = 2: User should provide the (M1+i-1)-th rows of Jacobian matrix at X() in YYp(i-1,j-1) (i = 1 to M2-M1+1, j = 1 to N). Do not alter any variables other than YYp(). = 3: To be returned with IRev = 3 on each iteration if Iout = 1. Output the intermediate result (X(), NFcall, NFjcall, Niter, S, etc.). Do not alter any variables. [in] Iout (Optional) Whether the intermediate result output is necessary. (default = 0) = 0: Output is not necessary (never return with IRev = 3). = 1: Return on each iteration with IRev = 3 to output the intermediate results. (If other value is specified, Iout = 0 will be assumed) [out] Info2 (Optional) Sub-code for Info = 0. = 1: X convergence. = 2: Relative function convergence. = 3: Both X and relative function convergence. = 4: Absolute function convergence. [out] NFcall (Optional) Number of function evaluations with IRev = 1 including those for computing the covariance matrix. [out] NFjcall (Optional) Number of Jacobian evaluations with IRev = 2 including those for computing the covariance matrix. [out] Niter (Optional) Number of iterations. [out] S (Optional) Residual sum of squares at the obtained solution vector X(). [out] NFGcal (Optional) Invocation counter Nf which is used for subroutines F and Fj in N2g. [in] Rtol (Optional) Relative function convergence tolerance. (Eps <= Rtol <= 0.1) (default = 1e-10) (Eps shows the machine epsilon hereafter) (If Rtol < Eps or Rtol > 0.1, the default value will be used) [in] Atol (Optional) Absolute function convergence tolerance. (default = 1e-20) (If Atol < 0, the default value will be used) [in] MaxFcall (Optional) Maximum number of function evaluations of F. (default = 200) (If MaxFcall <= 0, the default value will be used) [in] MaxIter (Optional) Maximum number of iterations. (default = 150) (If MaxIter <= 0, the default value will be used) [in] Dtype (Optional) Choice of adaptive scaling. (Dtype = 0, 1 or 2) (default = 1) = 0: Disable adaptive scaling. (scale factor = 1) = 1: Enable adaptive scaling during all iterations. = 2: Enable adaptive scaling during the first iteration and scale factor is left unchanged thereafter. (For other values, the default value will be used) [in] Dfac (Optional) Factor for adaptive scaling (0 <= Dfac <= 1) (default = 0.6) A scale factor D(i) is chosen by adaptive scaling so that D(i)*X(i) has about the same magnitude for all i. Let D1(i) = max(||Ji||, Dfac*D(i)) where ||Ji|| is the 2-norm of the i-th column of Jacobian matrix, then D(i) is chosen as follows.   if D1(i) >= Dtol: D(i) = D1(i)   if D1(i) < Dtol: D(i) = D0 (If Dfac < 0 or Dfac > 1, the default value will be used) [in] Dtol (Optional) Tolerance for adaptive scaling. (Dtol > 0) (default = 1.0e-6) (If Dtol <= 0, the default value will be used) [in] D0 (Optional) Initial value for adaptive scaling. (D0 > 0) (default = 1) (If D0 <= 0, the default value will be used) [in] Tuner1 (Optional) Parameter to check for false convergence. (0 <= Tuner1 <= 0.5) (default = 0.1) (If Tuner1 < 0 or Tuner1 > 0.5, the default value will be used) [in] Xctol (Optional) X convergence tolerance. (0 <= Xctol <= 1) (default = Eps^(1/2)) (If Xctol < 0 or Xctol > 1, the default value will be used) [in] Xftol (Optional) False convergence tolerance. (0 <= Xftol <= 1) (default = 100*Eps) (If Xftol < 0 or Xftol > 1, the default value will be used) [in] Lmax0 (Optional) Maximum 2-norm allowed for scaled very first step. (Lmax0 > 0) (default = 1) (If Lmax0 <= 0, the default value will be used) [in] Lmaxs (Optional) [in] Sctol (Optional) Lmaxs and Sctol are the singular convergence test parameters. (Lmaxs > 0, 0 <= Sctol <= 1) (default: Lmaxs = 1, Sctol = 1e-10) If the function reduction predicted for a step of length bounded by Lmaxs is less than Sctol*abs(f), returns with Info = 7 (f is the function value at the start of the current iteration). (If Lmaxs <= 0, the default value will be used) (If Sctol < 0 or Sctol > 1, the default value will be used) Reference netlib/port Example Program Approximate the following data with model function f(x) = c1*(1 - exp(-c2*x)). Determine two parameters c1 and c2 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. Lower and upper bounds on solution are 100 <= c1 <= 500 and 0 <= c2 <= 0.1. Sub FN2p(M As Long, Md1 As Long, M1 As Long, M2 As Long, N As Long, X() As Double, Nf As Long, Fvec() As Double) 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 For I = 0 To M2 - M1 Fvec(I) = Ydata(M1 + I - 1) - X(0) * (1 - Exp(-Xdata(M1 + I - 1) * X(1))) Next End Sub Sub JN2p(M As Long, Md1 As Long, M1 As Long, M2 As Long, N As Long, X() As Double, Nf As Long, Fjac() As Double) 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 For I = 0 To M2 - M1 Fjac(I, 0) = Exp(-Xdata(M1 + I - 1) * X(1)) - 1 Fjac(I, 1) = -Xdata(M1 + I - 1) * X(0) * Exp(-X(1) * Xdata(M1 + I - 1)) Next End Sub Sub Ex_N2pb_r() Const M = 4, Md = 2, N = 2 Dim X(N - 1) As Double, B(1, N - 1) As Double, Info As Long, M1 As Long, M2 As Long Dim YY(M - 1) As Double, YYp(M - 1, N - 1) As Double, IRev As Long X(0) = 500: X(1) = 0.0001 B(0, 0) = 200: B(1, 0) = 500 B(0, 1) = 0: B(1, 1) = 0.1 IRev = 0 Do Call N2pb_r(M, Md, N, X(), B(), Info, M1, M2, YY(), YYp(), IRev) If IRev = 1 Then Call FN2p(M, Md, M1, M2, N, X(), 0, YY()) ElseIf IRev = 2 Then Call JN2p(M, Md, M1, M2, N, X(), 0, YYp()) End If Loop While IRev <> 0 Debug.Print "C1, C2 =", X(0), X(1) Debug.Print "Info =", Info End Sub Example Results C1, C2 = 241.084896094489 5.44942234107799E-04 Info = 0