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;