NUMAL Section 3.1.1.2.1.4
BEGIN SECTION : 3.1.1.2.1.4 (December, 1978)
AUTHORS: A.BJOERCK AND G.H.GOLUB.
CONTRIBUTOR: J.KOOPMAN.
INSTITUTE: UNIVERSITY OF AMSTERDAM.
RECEIVED: 780701 .
BRIEF DESCRIPTION:
THIS SECTION CONTAINS TWO PROCEDURES FOR SOLVING A LINEAR LEAST
SQUARES PROBLEM SUBJECT TO LINEAR CONSTRAINTS:
LSQDECOMP , FOR THE QR-DECOMPOSITION OF A LEAST SQUARES MATRIX,
WHERE THIS MATRIX ALSO CONTAINS THE COEFFICIENTS OF
THE LINEAR CONSTRAINTS (LINEAR EQUATIONS);
LSQREFSOL , FOR THE SOLUTION OF THIS LEAST SQUARES PROBLEM, IF
THE MATRIX HAS BEEN DECOMPOSED BY LSQDECOMP.
KEYWORDS:
LEAST SQUARES,
LINEAR CONSTRAINTS,
HOUSEHOLDER TRIANGULARIZATION,
ITERATIVE REFINEMENT.
SUBSECTION: LSQDECOMP.
CALLING SEQUENCE:
THE HEADING OF THE PROCEDURE READS:
"PROCEDURE" LSQDECOMP(A, N, M, N1, AUX, AID, CI);
"VALUE"N,M,N1;"INTEGER"N,M,N1;"ARRAY" A,AUX,AID;
"INTEGER""ARRAY" CI;
"CODE" 34137;
THE MEANING OF THE FORMAL PARAMETERS IS:
A :<ARRAY IDENTIFIER>;
"ARRAY" A[1:N,1:M];
ENTRY: THE ORIGINAL LEAST SQUARES MATRIX, WHERE THE FIRST
N1 ROWS SHOULD FORM THE CONSTRAINT MATRIX (I.E. THE
FIRST N1 EQUATIONS ARE TO BE STRICTLY SATISFIED);
EXIT : IN THE UPPER TRIANGLE OF A (THE ELEMENTS A[I,J] WITH
I<J) THE SUPERDIAGONAL PART OF THE UPPER TRIANGULAR
MATRIX, PRODUCED BY HOUSEHOLDER TRANSFORMATIONS; IN
THE OTHER PART OF THE COLUMNS OF A THE SIGNIFICANT
ELEMENTS OF THE GENERATING VECTORS OF THE HOUSEHOLDER
MATRICES USED FOR THE TRIANGULARIZATION;
N :<ARITHMETIC EXPRESSION>;
NUMBER OF ROWS OF THE MATRIX;
M :<ARITHMETIC EXPRESSION>;
NUMBER OF COLUMNS OF THE MATRIX;
N1 :<ARITHMETIC EXPRESSION>;
NUMBER OF LINEAR CONSTRAINTS, I.E. THE FIRST N1 ROWS OF A
SET UP A SYSTEM OF N1 LINEAR EQUATIONS THAT MUST BE
STRICTLY SATISFIED (OF COURSE, IF THERE ARE NO CON-
STRAINTS, N1 MUST BE CHOSEN ZERO);
AUX :<ARRAY IDENTIFIER>;
"ARRAY" AUX[2:7];
ENTRY: AUX[2] CONTAINS A RELATIVE TOLERANCE FOR CALCULATING
THE DIAGONAL ELEMENTS OF THE UPPER TRIANGULAR MATRIX;
EXIT: AUX[3] CONTAINS THE NUMBER OF DIAGONAL ELEMENTS WHICH
ARE NOT NEGLIGIBLE, NORMAL EXIT AUX[3]=M;
AID :<ARRAY IDENTIFIER>;
"ARRAY" AID[1:M];
NORMAL EXIT (AUX[3]=M): THE DIAGONAL ELEMENTS OF THE UPPER
TRIANGULAR MATRIX PRODUCED BY THE HOUSEHOLDER TRANSFORMATION
CI :<ARRAY IDENTIFIER>;
"INTEGER""ARRAY" CI[1:M];
EXIT: THE PIVOTAL INDICES OF THE INTERCHANGES OF THE COLUMNS
OF THE GIVEN MATRIX;
PROCEDURES USED:
MATMAT = CP34013.
TAMMAT = CP34014.
ELMCOL = CP34023.
ICHCOL = CP34031.
RUNNING TIME:
ROUGHLY PROPORTIONAL TO N*M**2-M**3/3.
METHOD AND PERFORMANCE:
LET A DENOTE THE GIVEN MATRIX. LSQDECOMP PRODUCES AN N-TH ORDER
ORTHOGONAL MATRIX Q AND AN N*M UPPER TRIANGULAR MATRIX R SUCH THAT
R EQUALS QA WITH PERMUTED COLUMNS.
THE ORTHOGONAL MATRIX Q IS FORMED AS A PRODUCT OF AT MOST M TRANS-
FORMATIONS OF THE FORM (I-BETA U U'). THESE HOUSEHOLDER MATRICES
REDUCE A TO THE MATRIX R: AT THE K-TH STAGE THE K-TH COLUMN OF THE
(ALREADY MODIFIED) MATRIX IS INTERCHANGED WITH THE COLUMN OF
MAXIMUM EUCLIDEAN NORM.THESE INTERCHANGES ARE RECORDED IN THE ARRAY
CI. PREMATURE TERMINATION OCCURS IF AT SOME STAGE THE EUCLIDEAN
NORM OF THE PIVOTAL COLUMN IS LESS THAN SOME TOLERANCE (AUX[2])
TIMES THE MAXIMUM OF THE EUCLIDEAN NORMS OF THE COLUMNS OF THE
MATRIX.
LSQDECOMP DELIVERS THE UPPER TRIANGULAR MATRIX, WHERE THE DIAGONAL
ELEMENTS ARE STORED IN THE ARRAY AID AND THE OFF-DIAGONAL ELEMENTS
IN THE UPPER TRIANGULAR PART OF A. THE SIGNIFICANT ELEMENTS OF THE
GENERATING VECTORS OF THE HOUSEHOLDER TRANSFORMATIONS ARE STORED
IN THE COLUMNS OF A, I.E. ON AND BELOW THE LEADING DIAGONAL OF
A. IT SHOULD BE NOTED THAT FOR THE SOLUTION OF LEAST SQUARES
PROBLEMS, ONLY CALLS WITH N>=M ARE USEFUL.
REFERENCES:
A.BJOERCK AND G.H.GOLUB:
ITERATIVE REFINEMENT OF LEAST SQUARES SOLUTIONS BY HOUSEHOLDER
TRANSFORMATION, BIT 7 (1967), PP. 322-337
SUBSECTION: LSQREFSOL.
CALLING SEQUENCE:
THE HEADING OF THE PROCEDURE READS:
"PROCEDURE" LSQREFSOL(A, QR, N, M, N1, AUX,AID,CI,B,LDX,X,RES);
"VALUE"N,M,N1;"INTEGER"N,M,N1;"INTEGER""ARRAY"CI;"REAL" LDX;
"ARRAY"A,QR,AUX,AID,B,X,RES;
"CODE" 34138;
THE MEANING OF THE FORMAL PARAMETERS IS:
A :<ARRAY IDENTIFIER>;
"ARRAY" A[1:N,1:M];
THE ORIGINAL LEAST SQUARES MATRIX, WHERE THE FIRST N1
ROWS SET UP A SYSTEM OF LINEAR EQUATIONS THAT MUST BE
STRICTLY SATISFIED;
QR :<ARRAY IDENTIFIER>;
"ARRAY" QR[1:N,1:M];
THE QR-DECOMPOSITION OF THE ORIGINAL LEAST SQUARES
MATRIX AS DELIVERED BY A SUCCESSFUL CALL OF LSQDECOMP;
N :<ARITHMETIC EXPRESSION>;
NUMBER OF ROWS OF THE MATRICES A AND QR;
M :<ARITHMETIC EXPRESSION>;
NUMBER OF COLUMNS OF THE MATRICES A AND QR;
N1 :<ARITHMETIC EXPRESSION>;
NUMBER OF LINEAR CONSTRAINTS;
AUX :<ARRAY IDENTIFIER>;
"ARRAY" AUX[2:7];
ENTRY: AUX[2] CONTAINS A RELATIVE TOLERANCE AS A
CRITERION TO STOP ITERATIVE REFINING: IF THE
EUCLIDEAN NORM OF THE CORRECTION IS SMALLER
THAN AUX[2] TIMES THE CURRENT APPROXIMATION OF
THE SOLUTION, THE ITERATIVE REFINING IS STOPPED;
AUX[6]: MAXIMUM NUMBER OF ITERATIONS ALLOWED
(USUALLY AUX[6]=5 WILL BE SUFFICIENT);
EXIT : AUX[7]: THE NUMBER OF ITERATIONS PERFORMED;
AID :<ARRAY IDENTIFIER>;
"ARRAY" AID[1:M];
THE DIAGONAL ELEMENTS OF THE UPPER TRIANGULAR MATRIX
AS DELIVERED BY A SUCCESSFUL CALL OF LSQDECOMP;
CI :<ARRAY IDENTIFIER>;
"INTEGER""ARRAY" CI[1:M];
THE PIVOTAL INDICES AS PRODUCED BY LSQDECOMP;
B :<ARRAY IDENTIFIER>;
"ARRAY" B[1:N];
THE RIGHT-HAND SIDE OF THE LEAST SQUARES PROBLEM;
FIRST N1 ELEMENTS FORM THE RIGHT HAND SIDES OF THE
CONSTRAINTS;
LDX :<REAL VARIABLE>;
THE EUCLIDEAN NORM OF THE LAST CORRECTION OF THE
SOLUTION;
X :<ARRAY IDENTIFIER>;
"ARRAY" X[1:M];
EXIT: THE SOLUTION VECTOR;
RES :<ARRAY IDENTIFIER>;
"ARRAY" RES[1:N];
EXIT: THE RESIDUAL VECTOR CORRESPONDING TO THE SOLUTION;
PROCEDURES USED:
VECVEC = CP34010.
MATVEC = CP34011.
TAMVEC = CP34012.
ELMVECCOL = CP34021.
ICHCOL = CP34031.
LNG SUB = CP31106.
LNGMATVEC = CP34411.
LNGTAMVEC = CP34412.
RUNNING TIME:
ROUGHLY PROPORTIONAL TO C1*M*N-C2*M**2, WHERE C1 AND C2
ARE CONSTANTS.
METHOD AND PERFORMANCE:
LSQREFSOL SHOULD BE CALLED AFTER A SUCCESSFUL CALL OF LSQDECOMP
(I.E. AUX[3]=M). LSQREFSOL YIELDS THE LEAST SQUARES SOLUTION OF
THE OVERDETERMINED SYSTEM WITH THE DECOMPOSED COEFFICIENT MATRIX
IN THE ARRAY A AND THE RIGHT-HAND SIDE IN ARRAY B. THE ORIGINAL
LEAST SQUARES MATRIX ALSO CONTAINS THE LINEAR CONSTRAINTS (THE
FIRST N1 ROWS OF THIS MATRIX SET UP A SYSTEM THAT MUST BE STRICT-
LY SATISFIED). FIRST, THE ORTHOGONAL TRANSFORMATION WITH THE
HOUSEHOLDER MATRICES IS PERFORMED ON THE RIGHT-HAND SIDE. NEXT
THE SYSTEM IS SOLVED BY MEANS OF A BACK SUBSTITUTION. IN THIS
WAY THE FIRST APPROXIMATION OF THE SOLUTION (WITH PERMUTED
COLUMNS) IS OBTAINED.
AFTER THIS AN ITERATIVE PROCESS REFINES THE APPROXIMATION
UNTIL THE EUCLIDEAN NORM OF THE CORRECTION VECTOR IS NEGLIGIBLY
SMALL COMPARED TO THE APPROXIMATION OR UNTIL THE MAXIMUM NUMBER
OF ITERATIONS IS REACHED. FOR A MORE DETAILED DESCRIPTION OF
THE ITERATIVE IMPROVEMENT, SEE REF[1].
AFTER THE ITERATIVE PROCESS AN APPROXIMATION TO THE SOLUTION IS
FOUND. HOWEVER, THE ORDER OF THE COMPONENTS POSSIBLY IS NOT
CORRECT. THEREFORE THIS ORDER IS RESTORED AT THE END OF THE
PROCEDURE (SEE ALSO METHOD AND PERFORMANCE OF LSQDECOMP).
THE LEAST SQUARES SOLUTIONS OF SEVERAL OVERDETERMINED SYSTEMS
WITH THE SAME CONSTRAINTS AND COEFFICIENT MATRIX CAN BE SOLVED
BY SUCCESSIVE CALLS OF LSQREFSOL WITH DIFFERENT RIGHT-HAND
SIDES.
REFERENCES:
[1] A.BJOERCK AND G.H.GOLUB:
ITERATIVE REFINEMENT OF LEAST SQUARES SOLUTIONS BY HOUSEHOLDER
TRANSFORMATION, BIT 7 (1967), PP. 322-337.
EXAMPLE OF USE:
THE FOLLOWING PROGRAM SOLVES THE PROBLEM:
MINIMIZE //B - A*X//
UNDER X1 + 1000*X2 + 5*X3 = 2016 ,WHERE
A IS THE MATRIX
1 0 8
0 3 2
1 2 "-5
0 0 0
AND THE VECTOR B = (25, 12, 5.00003, 1)' , X = (X1, X2, X3)' ;
"BEGIN" "INTEGER" N,M,N1,I,J;
N := 5; M := 3; N1 := 1;
"BEGIN""INTEGER""ARRAY" CI[1:M];
"ARRAY" AUX[2:7],QR,A[1:N,1:M],B,RES[1:N],AID,X[1:M];
"REAL" LDX;
A[1,1] := 1; A[1,2] := 1000; A[1,3] := 5;
A[2,1] := 1; A[2,2] := 0; A[2,3] := 8;
A[3,1] := 0; A[3,2] := 3; A[3,3] := 2;
A[4,1] := 1; A[4,2] := 2; A[4,3] := "-5;
A[5,1]:=A[5,2]:=A[5,3]:=0;
B[1] := 2016; B[2] :=25; B[3] :=12; B[4] := 5.00003;
B[5] := 1; AUX[2] := "-14; AUX[6] := 5;
"FOR" I:=1 "STEP" 1 "UNTIL" N "DO"
"FOR" J:=1 "STEP" 1 "UNTIL" M "DO" QR[I,J] := A[I,J];
LSQDECOMP(QR,N,M,N1,AUX,AID,CI);
LSQREFSOL(A,QR,N,M,N1,AUX,AID,CI,B,LDX,X,RES);
OUTPUT(61,"(""("THE SOLUTION VECTOR: ")",//")");
"FOR" I:=1 "STEP" 1 "UNTIL" M "DO" OUTPUT(61,"("/")",X[I]);
OUTPUT(61,"("//,"("THE RESIDUAL VECTOR: ")",//")");
"FOR" J:=N1+1 "STEP" 1 "UNTIL" N "DO" OUTPUT(61,"("/")",RES[J]);
OUTPUT(61,"("///,"("NUMBER OF ITERATIONS: ")",D,//,
"("NORM LAST CORRECTION OF X: ")",N")",AUX[7],LDX)
"END"
"END"
DELIVERS:
THE SOLUTION VECTOR:
+1.0000000000000"+000
+2.0000000000000"+000
+3.0000000000000"+000
THE RESIDUAL VECTOR:
-5.2734444477081"-016
+2.1280091641666"-015
+5.3479806840033"-016
+1.0000000000000"+000
NUMBER OF ITERATIONS: 2
NORM LAST CORRECTION OF X: +2.1657844626990"-015
SOURCE TEXT(S):
"CODE" 34137;
"CODE" 34138;