NUMAL Section 3.1.1.2.1.2
BEGIN SECTION : 3.1.1.2.1.2 (May, 1974)
AUTHOR : T.J. DEKKER.
CONTRIBUTOR : J. KOK.
INSTITUTE : MATHEMATICAL CENTRE.
RECEIVED : 731015.
BRIEF DESCRIPTION :
THIS SECTION CONTAINS TWO PROCEDURES :
A) LSQSOL, FOR THE SOLUTION OF A LINEAR LEAST SQUARES PROBLEM IF
THE COEFFICIENT MATRIX HAS BEEN DECOMPOSED BY LSQORTDEC
(SECTION 3.1.1.2.1.1.);
B) LSQORTDECSOL, FOR THE SOLUTION OF A LINEAR LEAST SQUARES PROBLEM
BY HOUSEHOLDER TRIANGULARIZATION WITH COLUMN INTERCHANGES AND FOR
THE CALCULATION OF THE DIAGONAL OF THE INVERSE OF M'M, WHERE
M IS THE COEFFICIENT MATRIX.
KEY WORDS :
LINEAR LEAST SQUARES PROBLEM,
HOUSEHOLDER TRIANGULARIZATION.
SUBSECTION : LSQSOL.
CALLING SEQUENCE:
THE HEADING OF THE PROCEDURE IS :
"PROCEDURE" LSQSOL(A, N, M, AID, CI, B); "VALUE" N, M;
"INTEGER" N, M; "INTEGER""ARRAY" CI; "ARRAY" A, AID, B;
"CODE" 34131;
THE MEANING OF THE FORMAL PARAMETERS IS :
A, N, M, AID, CI : SEE 'CALLING SEQUENCE' OF LSQORTDEC
(SECTION 3.1.1.2.1.1.); THE CONTENTS OF THE ARRAYS
A, AID AND CI SHOULD BE PRODUCED BY A SUCCESSFUL CALL
OF LSQORTDEC, I.E. IF AUX[3] = M;
B : <ARRAY IDENTIFIER>;
"ARRAY" B[1 : N];
ENTRY : B CONTAINS THE RIGHT HAND SIDE OF A LINEAR
LEAST SQUARES PROBLEM;
EXIT : B[1 : M] CONTAINS THE SOLUTION OF THE PROBLEM;
B[M + 1 : N] CONTAINS A VECTOR WITH EUCLIDEAN
LENGTH EQUAL TO THE EUCLIDEAN LENGTH OF THE RESIDUAL
VECTOR.
PROCEDURES USED :
MATVEC = CP34011,
TAMVEC = CP34012,
ELMVECCOL = CP34021.
RUNNING TIME:
(C5 * M + C6) * N;
THE CONSTANTS C5 AND C6 DEPEND UPON THE ELEMENTARY
ARITHMETIC OF THE COMPUTER.
LANGUAGE : ALGOL 60.
METHOD AND PERFORMANCE :
LSQSOL SHOULD BE CALLED AFTER A SUCCESSFUL CALL OF LSQORTDEC
(SECTION 3.1.1.2.1.1.), I.E. IF AUX[3] = M. LSQSOL YIELDS
THE LEAST SQUARES SOLUTION OF THE OVERDETERMINED SYSTEM WITH THE
DECOMPOSED COEFFICIENT MATRIX IN ARRAY A AND THE RIGHT HAND SIDE IN
ARRAY B.
FIRST THE ORTHOGONAL TRANSFORMATION WITH THE HOUSEHOLDER MATRICES
IS PERFORMED ON THE RIGHT HAND SIDE. NEXT THE SYSTEM OF THE FIRST M
EQUATIONS AND WITH AN UPPER-TRIANGULAR COEFFICIENT MATRIX IS SOLVED
BY BACK SUBSTITUTION, YIELDING A SOLUTION WITH M PERMUTED
COMPONENTS DUE TO THE COLUMN INTERCHANGES OF THE TRIANGULARIZATION.
FINALLY THE ORDER OF THE M COMPONENTS IS RESTORED. SEE ALSO METHOD
AND PERFORMANCE OF LSQORTDEC (SECTION 3.1.1.2.1.1.).
THE LEAST SQUARES SOLUTIONS OF SEVERAL OVERDETERMINED SYSTEMS WITH
THE SAME COEFFICIENT MATRIX CAN BE OBTAINED BY SUCCESSIVE CALLS OF
LSQSOL WITH DIFFERENT RIGHT HAND SIDES.
EXAMPLE OF USE :
THE NEXT PROGRAM SOLVES THE SYSTEM
- 2 * X1 + X2 = 0
- X1 + X2 = 1
X1 + X2 = 2
2 * X1 + X2 = 2
X1 + 2 * X2 = 3
"BEGIN""COMMENT" 730912, TEST LSQORTDEC, LSQSOL, LSQDGLINV;
"ARRAY" A, C[1 : 5,1 : 2], B, X[1 : 5], DIAG, AID[1 : 2],
AUX[2 : 5];
"INTEGER""ARRAY" PIV[1 : 2];
"INTEGER" I, J;
"REAL" H;
"REAL""PROCEDURE" SUM(I, A, B, X); "VALUE" A, B;
"INTEGER" I, A, B; "REAL" X;
"BEGIN""REAL" S; S:= 0; "FOR" I:= A "STEP" 1 "UNTIL" B "DO"
S:= S + X; SUM:= S
"END" SUM;
AUX[2]:= "-12; I:= J:= 1;
"FOR" H:= - 2, - 1, 1, 2, 1, 1, 1, 1, 1, 2 "DO"
"BEGIN" A[I,J]:= C[I,J]:= H; "IF" I < 5 "THEN" I:= I + 1 "ELSE"
"BEGIN" I:= 1; J:= J + 1 "END"
"END";
"FOR" H:= 0, 1, 2, 2, 3 "DO"
"BEGIN" B[I]:= X[I]:= H; I:= I + 1 "END";
LSQORTDEC(A, 5, 2, AUX, AID, PIV);
"IF" AUX[3] = 2 "THEN"
"BEGIN" LSQSOL(A, 5, 2, AID, PIV, X);
LSQDGLINV(A, 2, AID, PIV, DIAG);
OUTPUT(61, "("/, "("AUX[2, 3, 5] = ")" +.4D"+DD5B, 3ZD5B,
+.4D"+DD/, "("LSQ SOLUTION :")", 2(2B+.8D"+DD), /
"("RESIDUE (DELIVERED) :")" +.8D"+DD/,
"("RESIDUE (CHECKED) :")" +.8D"+DD/,
"("DIAGONAL OF INVERSE M'M :")", 2(2B+.8D"+DD)")",
AUX[2], AUX[3], AUX[5], X[1], X[2],
SQRT(VECVEC(3, 5, 0, X, X)),
SQRT(SUM(I, 1, 5, (B[I] - C[I,1] * X[1] - C[I,2] * X[2])
** 2)), DIAG[1], DIAG[2])
"END"
"END"
DELIVERS :
AUX[2, 3, 5] = +.1000"-11 2 +.3317"+01
LSQ SOLUTION : +.50000000"+00 +.12500000"+01
RESIDUE (DELIVERED) :+.50000000"+00
RESIDUE (CHECKED) :+.50000000"+00
DIAGONAL OF INVERSE M'M : +.95238095"-01 +.13095238"+00
SUBSECTION : LSQORTDECSOL.
CALLING SEQUENCE :
THE HEADING OF THE PROCEDURE IS :
"PROCEDURE" LSQORTDECSOL(A, N, M, AUX, DIAG, B); "VALUE" N, M;
"INTEGER" N, M; "ARRAY" A, AUX, DIAG, B;
"CODE" 34135;
THE MEANING OF THE FORMAL PARAMETERS IS :
A : <ARRAY IDENTIFIER>;
"ARRAY" A[1 : N,1 : M];
ENTRY : A CONTAINS THE COEFFICIENT MATRIX OF THE
LINEAR LEAST SQUARES PROBLEM;
EXIT : IN THE UPPER TRIANGLE OF A (THE ELEMENTS
A[I,J] WITH I < J) THE SUPERDIAGONAL
ELEMENTS OF THE UPPER-TRIANGULAR MATRIX, PRODUCED BY
THE HOUSEHOLDER TRANSFORMATION; IN THE OTHER PART OF
THE COLUMNS OF A THE SIGNIFICANT ELEMENTS OF THE
GENERATING VECTORS OF THE HOUSEHOLDER MATRICES USED
FOR THE HOUSEHOLDER TRIANGULARIZATION;
N : <ARITHMETIC EXPRESSION>;
NUMBER OF ROWS OF THE MATRIX;
M : <ARITHMETIC EXPRESSION>;
NUMBER OF COLUMNS OF THE MATRIX (N >= M);
AUX : <ARRAY IDENTIFIER>;
"ARRAY" AUX[2 : 5];
ENTRY : AUX[2] CONTAINS A RELATIVE TOLERANCE USED FOR
CALCULATING THE DIAGONAL ELEMENTS OF THE
UPPER-TRIANGULAR MATRIX;
EXIT :
AUX[3] DELIVERS THE NUMBER OF THE DIAGONAL ELEMENTS OF
THE UPPER-TRIANGULAR MATRIX WHICH ARE FOUND NOT
NEGLIGIBLE; NORMAL EXIT AUX[3] = M;
AUX[5] := THE MAXIMUM OF THE EUCLIDEAN NORMS OF THE
COLUMNS OF THE GIVEN MATRIX;
DIAG : <ARRAY IDENTIFIER>;
"ARRAY" DIAG[1 : M];
EXIT : THE DIAGONAL ELEMENTS OF THE INVERSE OF M'M
WHERE M IS THE MATRIX OF THE LINEAR LEAST SQUARES
PROBLEM;
B : <ARRAY IDENTIFIER>;
"ARRAY" B[1 : N];
ENTRY : B CONTAINS THE RIGHT HAND SIDE OF A LINEAR
LEAST SQUARES PROBLEM;
EXIT : B[1 : M] CONTAINS THE SOLUTION OF THE PROBLEM;
B[M + 1 : N] CONTAINS A VECTOR WITH EUCLIDEAN LENGTH
EQUAL TO THE EUCLIDEAN LENGTH OF THE RESIDUAL VECTOR.
PROCEDURES USED :
LSQORTDEC = CP34134,
LSQDGLINV = CP34132,
LSQSOL = CP34131.
REQUIRED CENTRAL MEMORY :
EXECUTION FIELD LENGTH : AN INTEGER ARRAY AND A REAL ARRAY,
BOTH OF M ELEMENTS, ARE DECLARED.
RUNNING TIME : ROUGHLY PROPORTIONAL TO N * (M ** 2).
LANGUAGE : ALGOL 60.
METHOD AND PERFORMANCE :
LSQORTDECSOL SOLVES AN OVERDETERMINED SYSTEM OF N LINEAR EQUATIONS
IN M UNKNOWNS BY CALLING LSQORTDEC AND, IF THIS CALL WAS SUCCESSFUL
LSQDGLINV AND LSQSOL. LSQORTDECSOL DELIVERS THE LEAST SQUARES
SOLUTION AND THE DIAGONAL OF THE INVERSE OF M'M, WHERE M IS
THE COEFFICIENT MATRIX OF THE SYSTEM. SEE SECTION 3.1.1.2.1.1.,
AND LSQSOL (THIS SECTION).
EXAMPLE OF USE :
THE PROGRAM
"BEGIN""COMMENT" 730914, TEST LSQORTDECSOL;
"ARRAY" A, C[1 : 5,1 : 2], B, X[1 : 5], DIAG[1 : 2],
AUX[2 : 5];
"INTEGER" I, J;
"REAL" H;
"REAL""PROCEDURE" SUM(I, A, B, X); "VALUE" A, B;
"INTEGER" I, A, B; "REAL" X;
"BEGIN""REAL" S; S:= 0; "FOR" I:= A "STEP" 1 "UNTIL" B "DO"
S:= S + X; SUM:= S
"END" SUM;
AUX[2]:= "-12; I:= J:= 1;
"FOR" H:= - 2, - 1, 1, 2, 1, 1, 1, 1, 1, 2 "DO"
"BEGIN" A[I,J]:= C[I,J]:= H; "IF" I < 5 "THEN" I:= I + 1 "ELSE"
"BEGIN" I:= 1; J:= J + 1 "END"
"END";
"FOR" H:= 0, 1, 2, 2, 3 "DO"
"BEGIN" B[I]:= X[I]:= H; I:= I + 1 "END";
LSQORTDECSOL(A, 5, 2, AUX, DIAG, X);
"IF" AUX[3] = 2 "THEN"
OUTPUT(61, "("/, "("AUX[2, 3, 5] = ")" +.4D"+DD5B, 3ZD5B,
+.4D"+DD/, "("LSQ SOLUTION :")", 2(2B+.8D"+DD), /
"("RESIDUE (DELIVERED) :")" +.8D"+DD/,
"("RESIDUE (CHECKED) :")" +.8D"+DD/,
"("DIAGONAL OF INVERSE M'M :")", 2(2B+.8D"+DD)")",
AUX[2], AUX[3], AUX[5], X[1], X[2],
SQRT(VECVEC(3, 5, 0, X, X)),
SQRT(SUM(I, 1, 5, (B[I] - C[I,1] * X[1] - C[I,2] * X[2])
** 2)), DIAG[1], DIAG[2])
"END"
WHICH SOLVES THE PROBLEM OF THE EXAMPLE OF USE OF LSQSOL,
DELIVERS :
AUX[2, 3, 5] = +.1000"-11 2 +.3317"+01
LSQ SOLUTION : +.50000000"+00 +.12500000"+01
RESIDUE (DELIVERED) :+.50000000"+00
RESIDUE (CHECKED) :+.50000000"+00
DIAGONAL OF INVERSE M'M : +.95238095"-01 +.13095238"+00
SOURCE TEXT(S) :
"CODE" 34131;
"CODE" 34135;