NUMAL Section 5.1.3.1.3
BEGIN SECTION : 5.1.3.1.3 (December, 1975)
AUTHORS: J.C.P. BUS, B. VAN DOMSELAAR, J. KOK.
CONTRIBUTORS: B. VAN DOMSELAAR, J. KOK.
INSTITUTE: MATHEMATICAL CENTRE.
RECEIVED: 741101.
BRIEF DESCRIPTION:
THIS SECTION CONTAINS TWO PROCEDURES;
MARQUARDT CALCULATES THE LEAST SQUARES SOLUTION OF AN
OVERDETERMINED SYSTEM OF NONLINEAR EQUATIONS WITH MARQUARDT'S
METHOD;
GSSNEWTON CALCULATES THE LEAST SQUARES SOLUTION OF AN
OVERDETERMINED SYSTEM OF NONLINEAR EQUATIONS WITH THE GAUSS-NEWTON
METHOD;
THE USER HAS TO PROGRAM THE EVALUATION OF THE RESIDUAL VECTOR OF
THE SYSTEM AND THE JACOBIAN MATRIX OF ITS PARTIAL DERIVATIVES.
KEYWORDS:
NONLINEAR EQUATIONS,
LEAST SQUARES SOLUTION,
OVERDETERMINED SYSTEM,
MARQUARDT'S METHOD,
GAUSS-NEWTON METHOD,
CURVE FITTING.
SUBSECTION: MARQUARDT.
CALLING SEQUENCE:
THE HEADING OF THIS PROCEDURE IS:
"PROCEDURE" MARQUARDT(M, N, PAR, RV, JJINV, FUNCT, JACOBIAN, IN,
OUT); "VALUE" M, N; "INTEGER" M, N;
"ARRAY" PAR, RV, JJINV, IN, OUT; "BOOLEAN" "PROCEDURE" FUNCT;
"PROCEDURE" JACOBIAN;
"CODE" 34440;
THE MEANING OF THE FORMAL PARAMETERS IS:
M: <ARITHMETIC EXPRESSION>;
THE NUMBER OF EQUATIONS;
N: <ARITHMETIC EXPRESSION>;
THE NUMBER OF UNKNOWN VARIABLES; N SHOULD SATISFY N<=M;
PAR: <ARRAY IDENTIFIER>;
"ARRAY" PAR[1 : N];
THE UNKNOWN VARIABLES OF THE SYSTEM;
ENTRY: AN APPROXIMATION TO A LEAST SQUARES SOLUTION
OF THE SYSTEM;
EXIT: THE CALCULATED LEAST SQUARES SOLUTION;
RV: <ARRAY IDENTIFIER>;
"ARRAY" RV[1 : M];
EXIT: THE RESIDUAL VECTOR AT THE CALCULATED SOLUTION;
JJINV: <ARRAY IDENTIFIER>;
"ARRAY" JJINV[1 : N, 1 : N];
EXIT: THE INVERSE OF THE MATRIX J' * J WHERE J DENOTES
THE MATRIX OF PARTIAL DERIVATIVES DRV[I] / DPAR[J]
(I=1,...,M; J=1,...,N) AND J' DENOTES THE
TRANSPOSE OF J.
FUNCT: <PROCEDURE IDENTIFIER>;
THE HEADING OF THIS PROCEDURE SHOULD BE:
"BOOLEAN" "PROCEDURE" FUNCT(M, N, PAR, RV); "VALUE" M, N;
"INTEGER" M, N; "ARRAY" PAR, RV;
ENTRY: M, N, PAR;
M, N HAVE THE SAME MEANING AS IN THE PROCEDURE
MARQUARDT;
"ARRAY" PAR[1:N] CONTAINS THE CURRENT VALUES OF
THE UNKNOWNS AND SHOULD NOT BE ALTERED;
EXIT: "ARRAY" RV[1 : M];
UPON COMPLETION OF A CALL OF FUNCT, THIS ARRAY RV
SHOULD CONTAIN THE RESIDUAL VECTOR, OBTAINED WITH
THE CURRENT VALUES OF THE UNKNOWNS;
E.G. IN CURVE FITTING PROBLEMS:
RV[I] := THEORETICAL VALUE F(X[I], PAR) -
OBSERVED VALUE Y[I];
AFTER A SUCCESSFUL CALL OF FUNCT, THE BOOLEAN PROCEDURE
SHOULD DELIVER THE VALUE TRUE;
HOWEVER, IF FUNCT DELIVERS THE VALUE FALSE, THEN IT IS
ASSUMED THAT THE CURRENT ESTIMATES OF THE UNKNOWNS LIE
OUTSIDE A FEASIBLE REGION AND THE PROCESS IS TERMINATED
(SEE OUT[1]);
HENCE, PROPER PROGRAMMING OF FUNCT MAKES IT POSSIBLE TO
AVOID CALCULATION OF A RESIDUAL VECTOR WITH VALUES OF THE
UNKNOWN VARIABLES WHICH MAKE NO SENSE OR WHICH EVEN MAY
CAUSE OVERFLOW IN THE COMPUTATION;
JACOBIAN: <PROCEDURE IDENTIFIER>;
THE HEADING OF THIS PROCEDURE SHOULD BE:
"PROCEDURE" JACOBIAN(M, N, PAR, RV, JAC, LOCFUNCT);
"VALUE" M, N; "INTEGER" M, N; "ARRAY" PAR, RV, JAC;
"PROCEDURE" LOCFUNCT;
ENTRY: M, N, PAR, RV, LOCFUNCT;
FOR M,N,PAR SEE: FUNCT;
RV CONTAINS THE RESIDUAL VECTOR OBTAINED WITH THE
CURRENT VALUES OF THE UNKNOWNS AND SHOULD NOT BE
ALTERED;
A CALL OF LOCFUNCT(M,N,PAR,RV) IS EQUIVALENT WITH
A CALL OF THE USER-DEFINED PROCEDURE
FUNCT(M,N,PAR,RV), BUT, IN ADDITION, THIS CALL IS
COUNTED TO THE TOTAL NUMBER OF CALLS OF FUNCT
(SEE OUT[4]) AND, MOREOVER, IF FUNCT DELIVERS THE
VALUE FALSE THEN THE PROCESS IS TERMINATED;
EXIT: "ARRAY" JAC[1 : M, 1 : N];
UPON COMPLETION OF A CALL OF JACOBIAN, JAC SHOULD
CONTAIN THE PARTIAL DERIVATIVES DRV[I] / DPAR[J],
OBTAINED WITH THE CURRENT VALUES OF THE UNKNOWN
VARIABLES GIVEN IN PAR[1:N];
IT IS A PREREQUISITE FOR THE PROPER OPERATION OF THE
PROCEDURE MARQUARDT THAT THE PRECISION OF THE ELEMENTS OF
THE MATRIX JAC IS AT LEAST THE PRECISION DEFINED BY
IN[3] AND IN[4];
IN: <ARRAY IDENTIFIER>;
"ARRAY" IN[0 : 6];
ENTRY: IN THIS ARRAY THE USER SHOULD GIVE SOME DATA TO
CONTROL THE PROCESS;
IN[0]: THE MACHINE PRECISION;
FOR THE CYBER 73 A SUITABLE VALUE IS "-14;
IN[1], IN[2] ARE NOT USED BY THE PROCEDURE MARQUARDT;
IN[3], IN[4]:
THE RELATIVE AND ABSOLUTE TOLERANCE FOR THE
DIFFERENCE BETWEEN THE EUCLIDEAN NORM OF THE
ULTIMATE AND PENULTIMATE RESIDUAL VECTOR;
THE PROCESS IS TERMINATED IF THE IMPROVEMENT OF
THE SUM OF SQUARES IS LESS THAN
IN[3] * (SUM OF SQUARES) + IN[4] * IN[4];
THESE TOLERANCES SHOULD BE CHOSEN GREATER THAN
THE CORRESPONDING ERRORS OF THE CALCULATED
RESIDUAL VECTOR;
NOTE THAT THE EUCLIDEAN NORM OF THE RESIDUAL
VECTOR IS DEFINED AS THE SQUARE ROOT OF THE SUM
OF SQUARES;
IN[5]: THE MAXIMUM NUMBER OF CALLS OF FUNCT ALLOWED;
IN[6]: A STARTING VALUE USED FOR THE RELATION BETWEEN
THE GRADIENT AND THE GAUSS-NEWTON DIRECTION (SEE
[2]); IF THE PROBLEM IS WELL CONDITIONED THEN A
SUITABLE VALUE FOR IN[6] WILL BE 0.01; IF THE
PROBLEM IS ILL CONDITIONED THEN IN[6] SHOULD BE
GREATER, BUT THE VALUE OF IN[6] SHOULD SATISFY:
IN[0] < IN[6] <= 1/IN[0];
OUT: <ARRAY IDENTIFIER>; "ARRAY" OUT[1 : 7];
EXIT : IN ARRAY OUT SOME BY-PRODUCTS ARE DELIVERED;
OUT[1]: THIS VALUE GIVES INFORMATION ABOUT THE
TERMINATION OF THE PROCESS;
OUT[1]=0: NORMAL TERMINATION;
OUT[1]=1: THE PROCESS HAS BEEN BROKEN OFF,
BECAUSE THE NUMBER OF CALLS OF FUNCT
EXCEEDED THE NUMBER GIVEN IN IN[5];
OUT[1]=2: THE PROCESS HAS BEEN BROKEN OFF,
BECAUSE A CALL OF FUNCT DELIVERED THE
VALUE FALSE;
OUT[1]=3: FUNCT BECAME FALSE WHEN CALLED WITH
THE INITIAL ESTIMATES OF PAR[1:N];
THE ITERATION PROCESS WAS NOT STARTED
AND SO JJINV[1:N,1:N] CAN NOT BE USED;
OUT[1]=4: THE PROCESS HAS BEEN BROKEN OFF,
BECAUSE THE PRECISION ASKED FOR CAN
NOT BE ATTAINED; THIS PRECISION IS
POSSIBLY CHOSEN TOO HIGH, RELATIVE TO
THE PRECISION IN WHICH THE RESIDUAL
VECTOR IS CALCULATED (SEE IN[3]);
OUT[2]: THE EUCLIDEAN NORM OF THE RESIDUAL VECTOR
CALCULATED WITH VALUES OF THE UNKNOWNS DELIVERED;
OUT[3]: THE EUCLIDEAN NORM OF THE RESIDUAL VECTOR
CALCULATED WITH THE INITIAL VALUES OF THE
UNKNOWN VARIABLES;
OUT[4]: THE NUMBER OF CALLS OF FUNCT NECESSARY TO OBTAIN
THE CALCULATED RESULT;
OUT[5]: THE TOTAL NUMBER OF ITERATIONS PERFORMED; NOTE
THAT IN EACH ITERATION ONE EVALUATION OF THE
JACOBIAN MATRIX HAD TO BE MADE;
OUT[6]: THE IMPROVEMENT OF THE EUCLIDEAN NORM OF THE
RESIDUAL VECTOR IN THE LAST ITERATION STEP;
OUT[7]: THE CONDITION NUMBER OF J' * J , I.E. THE RATIO
OF ITS LARGEST TO SMALLEST EIGENVALUES;
DATA AND RESULTS:
IF THIS PROCEDURE IS USED FOR CURVE FITTING THEN THE RELATIVE
ACCURACY IN THE CALCULATION OF THE RESIDUAL VECTOR DEPENDS STRONGLY
ON THE ERRORS IN THE EXPERIMENTAL DATA AND THIS SHOULD BE REFLECTED
IN THE PARAMETERS IN[3] AND IN[4];
THE MATRIX JJINV CAN BE USED IF SOME STATISTICAL INFORMATION
ABOUT THE FITTED PARAMETERS IS REQUIRED; THE STANDARD DEVIATION,
COVARIANCE MATRIX AND CORRELATION MATRIX MAY BE CALCULATED EASILY
FROM JJINV ;
PROCEDURES USED:
MULCOL = CP31022,
DUPVEC = CP31030,
VECVEC = CP34010,
MATVEC = CP34011,
TAMVEC = CP34012,
MATTAM = CP34015,
QRISNGVALDEC = CP34273.
REQUIRED CENTRAL MEMORY:
EXECUTION FIELD LENGTH: ONE ARRAY OF LENGTH M * N AND FOUR
ARRAYS OF LENGTH N ARE DECLARED;
RUNNING TIME:
THE NUMBER OF ITERATION STEPS DEPENDS STRONGLY ON THE PROBLEM TO BE
SOLVED; HOWEVER, THE WORK DONE EACH ITERATION STEP IS PROPORTIONAL
TO N CUBED;
LANGUAGE: ALGOL 60.
METHOD AND PERFORMANCE:
MARQUARDT USES MARQUARDT'S METHOD; FOR DETAILS SEE [2];
REFERENCES:
[1] D. W. MARQUARDT,
AN ALGORITHM FOR LEAST-SQUARES ESTIMATION OF NONLINEAR
PARAMETERS,
J. SIAM 11 (1963), PP.431-441.
[2] J. C. P. BUS, B. VAN DOMSELAAR, J. KOK,
NONLINEAR LEAST SQUARES ESTIMATION,
MATHEMATICAL CENTRE, AMSTERDAM. (TO APPEAR)
EXAMPLE OF USE:
THE PARAMETERS PAR[1 : 3] IN THE CURVE FITTING PROBLEM:
RV[I] = PAR[1] + PAR[2] * EXP(PAR[3] * X[I]) - Y[I]
WITH (X[I], Y[I]), I=1,...,6 AS THE EXPERIMENTAL DATA, MAY BE
OBTAINED BY THE FOLLOWING PROGRAM:
"BEGIN"
"INTEGER" I;
"ARRAY" IN[0:6],OUT[0:7],X,Y,RV[1:6],JJINV[1:3,1:3],PAR[1:3];
"BOOLEAN" "PROCEDURE" EXPFUNCT(M,N,PAR,RV);
"VALUE" M,N; "INTEGER" M,N; "ARRAY" PAR,RV;
"BEGIN" "INTEGER" I;
"FOR" I:=1 "STEP" 1 "UNTIL" M "DO"
"BEGIN" "IF" PAR[3]*X[I]>680 "THEN"
"BEGIN" EXPFUNCT:="FALSE";
"GOTO" STOP
"END";
RV[I]:=PAR[1]+PAR[2]*EXP(PAR[3]*X[I])-Y[I]
"END"; EXPFUNCT:="TRUE";
STOP:
"END" EXPFUNCT;
"PROCEDURE" JACOBIAN(M,N,PAR,RV,JAC,LOCFUNCT);
"VALUE" M,N; "INTEGER" M,N; "ARRAY" PAR,RV,JAC;
"PROCEDURE" LOCFUNCT;
"BEGIN" "INTEGER" I; "REAL" EX;
"FOR" I:=1 "STEP" 1 "UNTIL" M "DO"
"BEGIN" JAC[I,1]:=1;
JAC[I,2]:=EX:=EXP(PAR[3]*X[I]);
JAC[I,3]:=X[I]*PAR[2]*EX
"END"
"END" JACOBIAN;
IN[0]:="-14; IN[3]:="-4; IN[4]:="-1; IN[5]:=75; IN[6]:="-2;
"FOR" I:=1 "STEP" 1 "UNTIL" 6 "DO"
"BEGIN" INREAL(70,X[I]); INREAL(70,Y[I]) "END";
PAR[1]:=580; PAR[2]:=-180; PAR[3]:=-0.160;
MARQUARDT(6,3,PAR,RV,JJINV,EXPFUNCT,JACOBIAN,IN,OUT);
OUTPUT(61,"("3/,"("X[I] Y[I]")",/,6(B,+D,5B,3D.D,/),2/,
"("PARAMETERS")",/,3(+D.3D"+ZD,/),2/,"("OUT:")",/,7(5B+D.6D"+ZD,
/),2/,"("LAST RESIDUAL VECTOR")",/,6(6B+3Z.D,/)")",X[1],Y[1],
X[2],Y[2],X[3],Y[3],X[4],Y[4],X[5],Y[5],X[6],Y[6],PAR[1],PAR[2],
PAR[3],OUT[7],OUT[2],OUT[6],OUT[3],OUT[4],OUT[5],OUT[1],RV[1],
RV[2],RV[3],RV[4],RV[5],RV[6])
"END"
WITH THE DATA GIVEN IN THE X - Y - TABLE THIS PROGRAM DELIVERS:
X[I] Y[I]
-5 127.0
-3 151.0
-1 379.0
+1 421.0
+3 460.0
+5 426.0
PARAMETERS
+5.232" +2
-1.568" +2
-1.998" -1
OUT:
+7.220828" +7
+1.157156" +2
+1.728008" -3
+1.654588" +2
+2.300000" +1
+2.200000" +1
+0.000000" +0
LAST RESIDUAL VECTOR
-29.6
+86.6
-47.3
-26.2
-22.9
+39.5
SUBSECTION : GSSNEWTON.
CALLING SEQUENCE :
THE HEADING OF THE PROCEDURE READS :
"PROCEDURE" GSSNEWTON(M, N, PAR, RV, JJINV, FUNCT, JACOBIAN, IN,
OUT);
"VALUE" M, N; "INTEGER" M, N; "ARRAY" PAR, RV, JJINV, IN, OUT;
"BOOLEAN""PROCEDURE" FUNCT; "PROCEDURE" JACOBIAN;
"CODE" 34441;
THE MEANING OF THE FORMAL PARAMETERS IS :
M : <ARITHMETIC EXPRESSION>;
THE NUMBER OF EQUATIONS;
N : <ARITHMETIC EXPRESSION>;
THE NUMBER OF UNKNOWNS IN THE M EQUATIONS (N <= M);
PAR : <ARRAY IDENTIFIER>; "ARRAY" PAR[1 : N];
THE UNKNOWNS OF THE EQUATIONS.
ENTRY : AN APPROXIMATION TO A LEAST SQUARES SOLUTION
OF THE SYSTEM.
EXIT : THE CALULATED LEAST SQUARES SOLUTION;
RV : <ARRAY IDENTIFIER>; "ARRAY" RV[1 : M];
EXIT : THE RESIDUAL VECTOR OF THE SYSTEM AT THE
CALCULATED SOLUTION;
JJINV : <ARRAY IDENTIFIER>; "ARRAY" JJINV[1 : N,1 : N];
EXIT : THE INVERSE OF THE MATRIX J' * J, WHERE J
IS THE JACOBIAN MATRIX AT THE SOLUTION AND J' IS
J TRANSPOSED;
FUNCT : <PROCEDURE IDENTIFIER>;
THE HEADING OF THIS PROCEDURE SHOULD BE :
"BOOLEAN""PROCEDURE" FUNCT(M, N, PAR, RV); "VALUE" M,
N; "INTEGER" M, N; "ARRAY" PAR, RV;
ENTRY: M, N, PAR;
M, N HAVE THE SAME MEANING AS IN THE PROCEDURE
GSSNEWTON;
"ARRAY" PAR[1:N] CONTAINS THE CURRENT VALUES OF
THE UNKNOWNS AND SHOULD NOT BE ALTERED.
EXIT: "ARRAY" RV[1 : M];
UPON COMPLETION OF A CALL OF FUNCT, THIS ARRAY RV
SHOULD CONTAIN THE RESIDUAL VECTOR, OBTAINED WITH
THE CURRENT VALUES OF THE UNKNOWNS.
THE PROGRAMMER OF FUNCT MAY DECIDE THAT SOME CURRENT
ESTIMATES OF THE UNKNOWNS LIE OUTSIDE A FEASIBLE
REGION; IN THIS CASE FUNCT SHOULD DELIVER THE VALUE
FALSE AND THE PROCESS IS TERMINATED (SEE OUT[1]).
OTHERWISE FUNCT SHOULD DELIVER THE VALUE TRUE;
JACOBIAN : <PROCEDURE IDENTIFIER>;
THE HEADING OF THIS PROCEDURE SHOULD BE :
"PROCEDURE" JACOBIAN(M, N, PAR, RV, JAC, LOCFUNCT);
"VALUE" M, N; "INTEGER" M, N; "ARRAY" PAR, RV, JAC;
"PROCEDURE" LOCFUNCT;
THE MEANING OF THE PARAMETERS OF JACOBIAN IS :
M, N : SEE GSSNEWTON.
PAR : <ARRAY IDENTIFIER>; "ARRAY" PAR[1 : N];
ENTRY : CURRENT ESTIMATES OF THE UNKNOWNS.
THESE VALUES SHOULD NOT BE CHANGED.
RV : <ARRAY IDENTIFIER>; "ARRAY" RV[1 : M];
ENTRY : THE RESIDUAL VECTOR OF THE SYSTEM OF
EQUATIONS CORRESPONDING TO THE VECTOR OF UNKNOWNS
AS GIVEN IN PAR.
EXIT : THE ENTRY VALUES.
JAC : <ARRAY IDENTIFIER>; "ARRAY" JAC[1 : M, 1 : N];
EXIT : THE JACOBIAN MATRIX AT THE CURRENT
ESTIMATES GIVEN IN PAR, I.E. THE MATRIX OF PARTIAL
DERIVATIVES
D(RV)[I] / DPAR[J], I = 1(1)M, J = 1(1)N.
LOCFUNCT : <PROCEDURE IDENTIFIER>; THE HEADING OF THIS
PROCEDURE IS THE SAME AS THE HEADING OF FUNCT.
A CALL OF THE PROCEDURE JACOBIAN SHOULD DELIVER THE
JACOBIAN MATRIX EVALUATED WITH THE CURRENT ESTIMATES
OF THE UNKNOWN VARIABLES GIVEN IN PAR
IN SUCH A WAY, THAT THE PARTIAL DERIVATIVE
D(RV)[I] / DPAR[J] IS DELIVERED IN JAC[I,J], I = 1(1)M,
J = 1(1)N.
FOR THE CALCULATION OF THE DERIVATIVES ONE CAN USE THE
VALUES OF THE CURRENT ESTIMATES OF THE
UNKNOWNS AS GIVEN IN PAR AND THE RESIDUAL VECTOR AS
GIVEN IN RV.
ONE CAN ALSO USE THE PROCEDURE FUNCT
(PARAMETER OF GSSNEWTON) THROUGH CALLS OF THE PROCEDURE
LOCFUNCT (PARAMETER OF JACOBIAN). THIS PARAMETER OF
JACOBIAN MAY BE USED WHEN THE JACOBIAN MATRIX IS
APPROXIMATED USING (FORWARD) DIFFERENCES.
AN APPROPRIATE PROCEDURE TO THIS PURPOSE IS JACOBNMF
(SECTION 4.3.2.1). SUCH A PROCEDURE MAY BE USED ONLY IF
THE MATRIX ELEMENTS ARE COMPUTED SUFFICIENTLY ACCURATE;
IN : <ARRAY IDENTIFIER>; "ARRAY" IN[0 : 7];
IN THIS ARRAY TOLERANCES AND CONTROL PARAMETERS SHOULD
BE GIVEN.
ENTRY :
IN[0] : THE MACHINE PRECISION. FOR CALCULATION ON THE
CYBER 73 A SUITABLE VALUE IS "-14.
IN[1], IN[2] :
RELATIVE AND ABSOLUTE TOLERANCE FOR THE STEP VECTOR
(RELATIVE TO THE VECTOR OF CURRENT ESTIMATES IN
PAR).
THE PROCESS IS TERMINATED IF IN SOME ITERATION (BUT
NOT THE FIRST) THE EUCLIDEAN NORM OF THE CALCULATED
NEWTON STEP IS LESS THAN IN[1] * NORM(PAR) + IN[2].
IN[1] SHOULD NOT BE CHOSEN SMALLER THAN IN[0].
IN[3] IS NOT USED BY THE PROCEDURE GSSNEWTON;
IN[4] : ABSOLUTE TOLERANCE FOR THE EUCLIDEAN NORM OF
THE RESIDUAL VECTOR. THE PROCESS IS TERMINATED WHEN
THIS NORM IS LESS THAN IN[4].
IN[5] : THE MAXIMUM ALLOWED NUMBER OF FUNCTION
EVALUATIONS (I.E. CALLS OF FUNCT).
IN[6] : THE MAXIMUM ALLOWED NUMBER OF HALVINGS OF A
CALCULATED NEWTON STEP VECTOR (SEE METHOD AND
PERFORMANCE). A SUITABLE VALUE IS 15.
IN[7] : THE MAXIMUM ALLOWED NUMBER OF SUCCESSIVE IN[6]
TIMES HALVED STEP VECTORS. SUITABLE VALUES ARE 1
AND 2;
OUT : <ARRAY IDENTIFIER>; "ARRAY" OUT[1 : 9];
IN ARRAY OUT INFORMATION ABOUT THE TERMINATION OF THE
PROCESS IS DELIVERED.
EXIT :
OUT[1] :
THE PROCESS WAS TERMINATED BECAUSE (OUT[1] = )
1.THE NORM OF THE RESIDUAL VECTOR IS SMALL WITH
RESPECT TO IN[4],
2.THE CALCULATED NEWTON STEP IS SUFFICIENTLY SMALL
(SEE IN[1], IN[2]),
3.THE CALCULATED STEP WAS COMPLETELY DAMPED (HALVED)
IN IN[7] SUCCESSIVE ITERATIONS,
4.OUT[4] EXCEEDS IN[5], THE MAXIMUM ALLOWED NUMBER OF
CALLS OF FUNCT,
5.THE JACOBIAN WAS NOT FULL-RANK (SEE OUT[8]),
6.FUNCT DELIVERED FALSE AT A NEW VECTOR OF
ESTIMATES OF THE UNKNOWNS,
7.FUNCT DELIVERED FALSE IN A CALL FROM JACOBIAN.
OUT[2] : THE EUCLIDEAN NORM OF THE LAST RESIDUAL
VECTOR.
OUT[3] : THE EUCLIDEAN NORM OF THE INITIAL RESIDUAL
VECTOR.
OUT[4] : THE TOTAL NUMBER OF CALLS OF FUNCT.
OUT[4] WILL BE LESS THAN IN[5] + IN[6].
OUT[5] : THE TOTAL NUMBER OF ITERATIONS.
OUT[6] : THE EUCLIDEAN NORM OF THE LAST STEP VECTOR.
OUT[7] : ITERATION NUMBER OF THE LAST ITERATION IN
WHICH THE NEWTON STEP WAS HALVED.
OUT[8], OUT[9] :
RANK AND MAXIMUM COLUMN NORM OF THE JACOBIAN MATRIX
IN THE LAST ITERATION, AS DELIVERED BY LSQORTDEC
(SECTION 3.1.1.2.1.1) IN AUX[3] AND AUX[5].
DATA AND RESULTS :
THE PROCEDURE GSSNEWTON CAN BE USED FOR APPROXIMATING AN EXACT OR A
LEAST SQUARES SOLUTION OF A SYSTEM OF NONLINEAR EQUATIONS. WHEN AN
EXACT SOLUTION IS REQUIRED, THE PROCEDURE MAY TERMINATE ONLY WITH
OUT[1] = 1, AND VERY SMALL VALUES SHOULD BE ASSIGNED TO IN[1] AND
IN[2]. WHEN A LEAST SQUARES SOLUTION IS REQUIRED, POSITIVE RESULTS
OF THE PROCEDURE ARE SIGNALED BY OUT[1] = 1 OR 2. WHENEVER THE
PROCEDURE TERMINATES WITH OUT[1] < 5, THEN THE INVERSE OF J' * J
(SEE MEANING OF THE PARAMETER JJINV) IS DELIVERED IN JJINV. IN
THAT CASE THE COVARIANCE MATRIX AND THE STANDARD DEVIATIONS OF THE
SOLUTION CAN BE CALCULATED.
FOR A CURVE FITTING PROBLEM, SAY :
"ESTIMATE PARAMETERS PAR[1], ... , PAR[N] OF A FUNCTION
"Y = F(X; PAR[1], ... , PAR[N]), WHEN A SET OF DATA (X[I],Y[I]),
"I = 1(1)M, HAS TO BE FITTED,"
THE FOLLOWING SYSTEM OF M EQUATIONS IN THE N UNKNOWN PARAMETERS
PAR[1], ... , PAR[N] CAN BE DERIVED :
F(X[I]; PAR[1], ... , PAR[N]) - Y[I] = 0, I = 1(1)M.
PROCEDURES USED :
VECVEC = CP34010,
DUPVEC = CP31030,
ELMVEC = CP34020,
LSQORTDEC = CP34134,
LSQSOL = CP34131,
LSQINV = CP34136.
REQUIRED CENTRAL MEMORY :
EXECUTION FIELD LENGTH : AN ARRAY OF (M + 1) * N ELEMENTS, FOUR
ARRAYS OF N ELEMENTS AND ONE ARRAY OF M ELEMENTS ARE DECLARED.
RUNNING TIME :
THE RUNNING TIME IS PROPORTIONAL TO THE TOTAL NUMBER OF CALLS OF
FUNCT. BESIDES, IN EACH ITERATION A LINEAR LEAST SQUARES SYSTEM
IS SOLVED (NUMBER OF OPERATIONS PROPORTIONAL TO M * (N SQUARED)).
LANGUAGE : ALGOL 60.
METHOD AND PERFORMANCE :
THE PROCEDURE GSSNEWTON IS BASED UPON THE GAUSS-NEWTON METHOD FOR
CALCULATING A LEAST SQUARES SOLUTION OF A SYSTEM OF NONLINEAR
EQUATIONS (SEE E.G. [1], [2]). IN SEVERAL ITERATIONS A STEP VECTOR
(FOR THE VECTOR OF UNKNOWNS) IS OBTAINED BY SOLVING A LINEARIZED
SYSTEM OF EQUATIONS. THIS STEP IS PERFORMED ONLY IF THE NORM OF THE
RESIDUAL VECTOR DECREASES. OTHERWISE THE NEWTON STEP VECTOR IS
HALVED A NUMBER OF TIMES UNTIL THE NORM OF THE RESIDUAL VECTOR IS
SMALLER THAN THE NORMS OF THE LAST AND SUBSEQUENT RESIDUAL VECTORS
(THIS PROCESS IS CALLED STEP SIZE CONTROL).
FOR FURTHER DETAILS SEE [3] (SEE ALSO [2]).
REFERENCES :
[1] H.O. HARTLEY :
THE MODIFIED GAUSS-NEWTON METHOD.
TECHNOMETRICS, V.3 (1961), PP. 269 - 280.
[2] H. SPAETH :
THE DAMPED TAYLOR'S SERIES METHOD FOR MINIMIZING A SUM OF
SQUARES AND FOR SOLVING SYSTEMS OF NONLINEAR EQUATIONS.
ALGORITHM 315, COLLECTED ALGORITHMS FROM CACM,
COMMUNICATIONS OF THE ACM, VOL. 10 (NOV. 1967), PP. 726 - 728.
[3] J.C.P. BUS, B. VAN DOMSELAAR, J. KOK :
NONLINEAR LEAST SQUARES ESTIMATION.
MATHEMATICAL CENTRE (TO APPEAR).
EXAMPLE OF USE :
THE PARAMETERS PAR[1 : 3] IN THE CURVE FITTING PROBLEM:
G[I] = PAR[1] + PAR[2] * EXP(PAR[3] * X[I]) - Y[I]
WITH (X[I], Y[I]), I=1,...,6 AS THE EXPERIMENTAL DATA, MAY BE
OBTAINED BY THE FOLLOWING PROGRAM:
"BEGIN"
"INTEGER" I;
"ARRAY" IN[0:7], OUT[1:9], X, Y, G[1:6], V[1:3, 1:3], PAR[1:3];
"BOOLEAN""PROCEDURE" EXPFUNCT(M, N, PAR, G);
"VALUE" M, N; "INTEGER" M, N; "ARRAY" PAR, G;
"BEGIN""INTEGER" I;
"FOR" I:= 1 "STEP" 1 "UNTIL" M "DO"
"BEGIN""IF" PAR[3] * X[I] > 680 "THEN"
"BEGIN" EXPFUNCT:= "FALSE"; "GO TO" STOP "END";
G[I]:= PAR[1] + PAR[2] * EXP(PAR[3] * X[I]) - Y[I]
"END"; EXPFUNCT:="TRUE";
STOP:
"END" EXPFUNCT;
"PROCEDURE" JACOBIAN(M, N, PAR, G, JAC, LOCFUNCT);
"VALUE" M, N; "INTEGER" M, N; "ARRAY" PAR, G, JAC;
"PROCEDURE" LOCFUNCT;
"BEGIN" "INTEGER" I; "REAL" EX;
"FOR" I:= 1 "STEP" 1 "UNTIL" M "DO"
"BEGIN" JAC[I,1]:=1; JAC[I,2]:= EX:= EXP(PAR[3] * X[I]);
JAC[I,3]:= X[I] * PAR[2] * EX
"END"
"END" JACOBIAN;
IN[0]:= "-14; IN[1]:= IN[2]:= "-6; IN[5]:= 75; IN[4]:="-10;
IN[6]:= 14; IN[7]:= 1;
"FOR" I:= 1 "STEP" 1 "UNTIL" 6 "DO"
"BEGIN" INREAL(70,X[I]); INREAL(70,Y[I]) "END";
PAR[1]:= 580; PAR[2]:= - 180; PAR[3]:= - 0.160;
GSSNEWTON(6, 3, PAR, G, V, EXPFUNCT, JACOBIAN, IN, OUT);
OUTPUT(61,"("3/4B,"("X[I] Y[I]")",/, 6(5B+D, 5B3D.D, /), 2/,
4B"("PARAMETERS")",/,3(4B+D.3D"+ZD,/),2/,4B"("OUT:")",/,
3(9B+D.6D"+ZD,/), 5(14B3ZD,/), 9B+D.6D"+ZD,2/4B,
"("LAST RESIDUAL VECTOR")",/,6(10B+3Z.D,/)")", X[1], Y[1],
X[2],Y[2],X[3],Y[3],X[4],Y[4],X[5],Y[5],X[6],Y[6],PAR[1],PAR[2],
PAR[3],OUT[6],OUT[2],OUT[3],OUT[4],OUT[5],OUT[1],OUT[7],OUT[8],
OUT[9], G[1],G[2],G[3],G[4],G[5],G[6])
"END"
WITH THE DATA GIVEN IN THE X - Y - TABLE THIS PROGRAM DELIVERS:
X[I] Y[I]
-5 127.0
-3 151.0
-1 379.0
+1 421.0
+3 460.0
+5 426.0
PARAMETERS
+5.233" +2
-1.569" +2
-1.997" -1
OUT:
+5.260478" -4
+1.157156" +2
+1.654588" +2
16
16
2
0
3
+2.339529" +3
LAST RESIDUAL VECTOR
-29.6
+86.6
-47.3
-26.2
-22.9
+39.5
SOURCE TEXTS :
"CODE" 34440;
"CODE" 34441;