NUMAL Section 3.1.1.3.1.4
BEGIN SECTION : 3.1.1.3.1.4 (December, 1979)
AUTHOR : D.T.WINTER
INSTITUTE : MATHEMATICAL CENTRE
RECEIVED : 731217
BRIEF DESCRIPTION :
THIS SECTION CONTAINS TWO PROCEDURES FOR THE
CALCULATION OF THE PSEUDO-INVERSE OF A MATRIX:
PSDINVSVD ASSUMES THAT THE MATRIX IS GIVEN AS SINGULAR VALUES
DECOMPOSITION.
PSDINV FIRST CALCULATES THIS DECOMPOSITION.
KEYWORDS :
PSEUDO-INVERSE
SINGULAR VALUES
SUBSECTION : PSDINVSVD
CALLING SEQUENCE :
THE HEADING OF THE PROCEDURE IS :
"PROCEDURE" PSDINVSVD(U, VAL, V, M, N, EM);
"VALUE" M, N; "INTEGER" M, N; "ARRAY" U, VAL, V, EM; "CODE" 34286;
THE MEANING OF THE FORMAL PARAMETERS IS :
U: <ARRAY IDENTIFIER>;
"ARRAY" U[1:M,1:N];
ENTRY: THE MATRIX U IN THE SINGULAR VALUES DECOMPOSITION
U * S * V';
EXIT: THE TRANSPOSE OF THE PSEUDO-INVERSE.
VAL: <ARRAY IDENTIFIER>;
"ARRAY" VAL[1:N];
THE SINGULAR VALUES.
V: <ARRAY IDENTIFIER>;
"ARRAY" V[1:N,1:N];
THE MATRIX V IN THE SINGULAR VALUES DECOMPOSITION U * S * V'.
M: <ARITHMETIC EXPRESSION>;
THE NUMBER OF ROWS OF U.
N: <ARITHMETIC EXPRESSION>;
THE NUMBER OF COLUMNS OF V.
EM: <ARRAY IDENTIFIER>;
"ARRAY" EM[6:6];
ENTRY: EM[6]:THE MINIMAL NON-NEGLECTABLE SINGULAR VALUE.
PROCEDURES USED :
MATVEC = CP34011
REQUIRED CENTRAL MEMORY : AN AUXILIARY ARRAY OF N REALS IS DECLARED
RUNNING TIME : ROUGHLY PROPORTIONAL TO M * N * N
METHOD AND PERFORMANCE :
THE PSEUDO-INVERSE IS CALCULATED IN TWO STEPS :
1. THE MATRIX X = VAL+ * U' IS CALCULATED, WHERE VAL+ DENOTES THE
DIAGONAL MATRIX OBTAINED FROM VAL BY PUTTING
VAL+[I,I] = 1/VAL[I] IF VAL[I] GREATER THAN OR EQUAL TO EM[6],
AND VAL+[I,I] = 0 OTHERWISE.
2. THE PSEUDO INVERSE (V * X) IS CALCULATED.
SUBSECTION : PSDINV
CALLING SEQUENCE :
THE HEADING OF THE PROCEDURE IS :
"INTEGER" "PROCEDURE" PSDINV(A, M, N, EM);
"VALUE" M, N; "INTEGER" M, N; "ARRAY" A, EM; "CODE" 34287;
PSDINV:= THE NUMBER OF SINGULAR VALUES NOT FOUND, I.E. ZERO IF ALL
SINGULAR VALUES ARE CALCULATED.
THE MEANING OF THE FORMAL PARAMETERS IS :
A <ARRAY IDENTIFIER>;
"ARRAY" A[1:M,1:N];
ENTRY : THE GIVEN MATRIX;
EXIT : THE TRANSPOSE OF THE PSEUDO-INVERSE;
M: <ARITHMETIC EXPRESSION>;
THE NUMBER OF ROWS OF A;
N: <ARITHMETIC EXPRESSION>;
THE NUMBER OF COLUMNS OF A, N<= M;
EM: <ARRAY IDENTIFIER>;
"ARRAY" EM[0:7];
ENTRY: EM[0]: THE MACHINE PRECISION;
EM[2]: THE RELATIVE PRECISION FOR THE SINGULAR VALUES;
EM[4]: THE MAXIMAL NUMBER OF ITERATIONS TO BE PERFORMED;
EM[6]: THE MINIMAL NON-NEGLECTABLE SINGULAR VALUE;
EXIT: EM[1]: THE INFINITY NORM OF THE MATRIX;
EM[3]: THE MAXIMAL NEGLECTED SUPERDIAGONAL ELEMENT;
EM[5]: THE NUMBER OF ITERATIONS PERFORMED IN THE SINGULAR
VALUES DECOMPOSITION;
EM[7]: THE NUMERICAL RANK OF THE MATRIX, I.E. THE NUMBER OF
SINGULAR VALUES GREATER THAN OR EQUAL TO EM[6];
PROCEDURES USED :
QRISNGVALDEC = CP34273
PSDINVSVD = CP34286
REQUIRED CENTRAL MEMORY :
AUXILIARY ARRAYS ARE DECLARED TO A TOTAL OF (N + 1) * N REALS
RUNNING TIME : ROUGHLY PROPORTIONAL TO (2M + N) * N * N
METHOD AND PERFORMANCE :
FIRST THE SINGULAR VALUES DECOMPOSITION IS CALCULATED, AND THEN THE
PSEUDO-INVERSE IS CALCULATED BY PSDINVSVD.
REFERENCES :
WILKINSON, J.H. AND C.REINSCH
HANDBOOK OF AUTOMATIC COMPUTATION, VOL.2
LINEAR ALGEBRA
HEIDELBERG (1971)
EXAMPLE OF USE :
FIRST WE GIVE A PROGRAM, AND THEN THE RESULTS OF THIS PROGRAM :
"BEGIN" "ARRAY" A[1:8,1:5], EM[0:7];
"INTEGER" I, J;
A[1,1]:=22; A[1,2]:= A[2,3]:=10; A[1,3]:= A[7,1]:= A[8,5]:=2;
A[1,4]:= A[3,5]:=3; A[1,5]:= A[2,2]:=7; A[2,1]:=14; A[2,5]:=8;
A[2,4]:= A[8,3]:=0; A[3,1]:= A[3,3]:= A[6,5]:=-1; A[3,2]:=13;
A[3,4]:=-11; A[4,1]:=-3; A[4,2]:= A[4,4]:= A[5,4]:= A[8,4]:=-2;
A[4,3]:=13; A[4,5]:= A[5,5]:= A[8,1]:=4; A[5,1]:= A[6,1]:=9;
A[5,2]:=8; A[5,3]:= A[6,2]:= A[7,5]:=1; A[6,3]:=-7;
A[6,4]:= A[7,4]:= A[8,2]:=5; A[7,2]:=-6; A[7,3]:=6;
EM[0]:="-14; EM[2]:="-12; EM[4]:=80; EM[6]:="-10;
I:= PSDINV(A, 8, 5, EM);
OUTPUT(61, "("4B, "("NUMBER SINGULAR VALUES NOT FOUND : ")",
3ZD,/, 4B, "("NORM : ")", N,/, 4B, "("MAX NEGL SUBD ELEM : ")",
N,/, 4B, "("NUMBER ITERATIONS : ")", 3ZD, /, 4B, "("RANK : ")",
3ZD, /")", I, EM[1], EM[3], EM[5], EM[7]);
OUTPUT(61, "("/, 4B, "("TRANSPOSE OF PSEUDO-INVERSE")", /,
4B, "("FIRST THREE COLUMNS")", /, /, 8(4B, 3(N), /), /, /,
4B, "("LAST TWO COLUMNS")", /, /, 8(15B, 2(N), /)")",
A[1,1], A[1,2], A[1,3], A[2,1], A[2,2], A[2,3], A[3,1], A[3,2],
A[3,3], A[4,1], A[4,2], A[4,3], A[5,1], A[5,2], A[5,3], A[6,1],
A[6,2], A[6,3], A[7,1], A[7,2], A[7,3], A[8,1], A[8,2], A[8,3],
A[1,4], A[1,5], A[2,4], A[2,5], A[3,4], A[3,5], A[4,4], A[4,5],
A[5,4], A[5,5], A[6,4], A[6,5], A[7,4], A[7,5], A[8,4], A[8,5])
"END"
NUMBER SINGULAR VALUES NOT FOUND : 0
NORM : +4.4000000000000"+001
MAX NEGL SUBD ELEM : +4.3977072741076"-014
NUMBER ITERATIONS : 6
RANK : 3
TRANSPOSE OF PSEUDO-INVERSE
FIRST THREE COLUMNS
+2.1129807692308"-002 +4.6153846153850"-003 -2.1073717948727"-003
+9.3108974358974"-003 +2.2115384615376"-003 +2.0528846153848"-002
-1.1097756410256"-002 +2.7403846153848"-002 -3.8862179487199"-003
-7.9166666666667"-003 -5.0000000000007"-003 +3.3750000000001"-002
+5.5128205128205"-003 +9.8076923076935"-003 -8.9743589743826"-004
+1.4318910256410"-002 -2.5961538461548"-003 -2.0136217948716"-002
+4.8958333333335"-003 -1.4999999999998"-002 +1.5312499999996"-002
+1.5064102564102"-003 +7.4038461538447"-003 -1.6987179487147"-003
LAST TWO COLUMNS
+7.6041666666662"-003 +3.8060897435894"-003
-2.0833333333295"-004 +1.0016025641026"-002
-2.7604166666667"-002 +4.2067307692303"-003
-5.4166666666662"-003 +1.0416666666667"-002
-5.0000000000005"-003 +3.2051282051275"-003
+1.2812500000000"-002 -6.2099358974354"-003
+1.2395833333332"-002 +2.6041666666656"-003
-4.9999999999993"-003 +1.6025641025649"-003
"EOP"
SOURCE TEXT(S):
"CODE" 34286;
"CODE" 34287;