NUMAL Section 3.1.1.2.1.3

BEGIN SECTION : 3.1.1.2.1.3 (October, 1974)

CONTRIBUTOR :    J. KOK.

INSTITUTE   :    MATHEMATICAL CENTRE.

RECEIVED    :    740617.

BRIEF DESCRIPTION    :

    THIS SECTION CONTAINS ONE PROCEDURE,
    LSQINV, FOR THE CALCULATION OF THE INVERSE OF THE MATRIX  S'S,
    WHERE  S  IS THE COEFFICIENT MATRIX OF A LINEAR LEAST SQUARES
    PROBLEM.

KEYWORDS   :

    INVERSE MATRIX,
    LINEAR LEAST SQUARES PROBLEM.

CALLING SEQUENCE    :

    THE HEADING OF THE PROCEDURE READS :
    "PROCEDURE" LSQINV(A, M, AID, C); "VALUE" M; "INTEGER" M;
    "ARRAY" A, AID; "INTEGER""ARRAY" C;
    "CODE" 34136;

    THE MEANING OF THE FORMAL PARAMETERS IS   :

        A   :   <ARRAY IDENTIFIER>;
                "ARRAY" A[1 : M, 1 : M];
                ENTRY :  IN THE UPPER TRIANGLE OF A (THE ELEMENTS
                A[I,J] WITH  1 <= I < J <= M) THE SUPERDIAGONAL
                ELEMENTS SHOULD BE GIVEN OF THE UPPERTRIANGULAR MATRIX
                THAT IS PRODUCED BY THE HOUSEHOLDER TRIANGULARIZATION
                IN A CALL OF THE PROCEDURE  LSQORTDEC  (SECTION
                3.1.1.2.1.1.) WITH A NORMAL EXIT (AUX[3] = M).
                SEE ALSO THE MEANING OF THE PARAMETER  AID;
                EXIT :  THE UPPER TRIANGLE OF THE (SYMMETRIC) INVERSE
                MATRIX IS DELIVERED IN THE UPPERTRIANGULAR ELEMENTS OF
                THE ARRAY A (A[I,J] FOR  1 <= I <= J <= M);
        M   :   <ARITHMETIC EXPRESSION>;
                NUMBER OF COLUMNS OF THE MATRIX OF THE LINEAR LEAST
                SQUARES PROBLEM;
        AID :   <ARRAY IDENTIFIER>;
                "ARRAY" AID[1 : M];
                ENTRY :  AID CONTAINS THE DIAGONAL ELEMENTS OF THE
                UPPERTRIANGULAR MATRIX THAT IS PRODUCED BY LSQORTDEC;
        C   :   <ARRAY IDENTIFIER>;
                "INTEGER""ARRAY" C[1 : M];
                ENTRY :  C CONTAINS THE PIVOTAL INDICES AS PRODUCED BY
                A CALL OF LSQORTDEC.

PROCEDURES USED    :

    CHLINV2   = CP34400,
    ICHCOL    = CP34031,
    ICHROW    = CP34032,
    ICHROWCOL = CP34033.

REQUIRED CENTRAL MEMORY    :

    EXECUTION FIELD LENGTH :  A REAL ARRAY OF M ELEMENTS IS
    DECLARED (IN THE CALL OF CHLINV2).

RUNNING TIME    :  PROPORTIONAL TO M CUBED.

LANGUAGE    :  ALGOL 60.

METHOD AND PERFORMANCE  :

    LSQINV  SHOULD  BE CALLED  AFTER  A SUCCESSFUL  CALL  OF  LSQORTDEC
    (SECTION  3.1.1.2.1.1.).  LSQINV CAN BE USED FOR THE CALCULATION OF
    THE COVARIANCE MATRIX OF A LINEAR LEAST SQUARES PROBLEM.
    LET  S  BE THE MATRIX  OF THE LEAST SQUARES  SYSTEM  WITH  PERMUTED
    COLUMNS  AND  Q * R  THE  HOUSEHOLDER DECOMPOSITION OF  S. THEN THE
    INVERSE  OF   S'S  ALSO IS  THE INVERSE  OF  R'R.  SINCE  R  IS THE
    CHOLESKY MATRIX OF  R'R,  THE INVERSE MATRIX IS COMPUTED  IN A CALL
    OF  CHLINV2  (SECTION 3.1.1.1.1.2.4.).  AFTERWARDS  THE  COVARIANCE
    MATRIX IS OBTAINED BY INTERCHANGES  OF THE COLUMNS AND ROWS  OF THE
    INVERSE MATRIX.

EXAMPLE OF USE  :

    THE FOLLOWING PROGRAM COMPUTES THE INVERSE  T  OF  S'S, WHERE  S IS
    THE COEFFICIENT  MATRIX  OF THE SYSTEM  IN THE EXAMPLE  OF  USE  OF
    LSQORTDEC AND  LSQDGLINV (SECTION 3.1.1.2.1.1.). THE DIAGONAL OF  S
    CAN BE COMPARED WITH THE RESULT OF  LSQDGLINV. TO CHECK THE ANSWERS
    S' * (S * T)  IS PRINTED.

    "BEGIN""COMMENT" JKOK, 740530, EXAMPLE OF USE OF  LSQORTDEC AND
        LSQINV;

        "ARRAY" A, C[1 : 5, 1 : 2], AID[1 : 2], T[1 : 2, 1 : 2],
        AUX[2 : 5];
        "INTEGER""ARRAY" PIV[1 : 2];
        "INTEGER" I, J; "REAL" H;

        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";

        LSQORTDEC(A, 5, 2, AUX, AID, PIV); "IF" AUX[3] = 2 "THEN"
        "BEGIN" LSQINV(A, 2, AID, PIV);
            T[1,1]:= A[1,1]; T[2,2]:= A[2,2]; T[2,1]:= T[1,2]:= A[1,2];
            "FOR" J:= 1, 2 "DO""FOR" I:= 1 "STEP" 1 "UNTIL" 5 "DO"
            A[I,J]:= MATMAT(1, 2, I, J, C, T);
            OUTPUT(61, "("/4B, "(" AUX[2, 3, 5] = ")",
                -.4D"+DD5B, 3ZD5B, +.4D"+DD/,
                2(/4B, 30S, /, 2(/4B, 2(2B+.8D"+DD)), /) ")",
                AUX[2], AUX[3], AUX[5],
                "(" INVERSE :")", T,
                "(" CHECK : S' * (S * T) :")",
                 TAMMAT(1, 5, 1, 1, C, A),
                 TAMMAT(1, 5, 1, 2, C, A),
                 TAMMAT(1, 5, 2, 1, C, A),
                 TAMMAT(1, 5, 2, 2, C, A) )
        "END"
    "END"

    OUTPUT  :

     AUX[2, 3, 5] =  .1000"-11        2     +.3317"+01

     INVERSE :

      +.95238095"-01  -.23809524"-01
      -.23809524"-01  +.13095238"+00

     CHECK : S' * (S * T) :

      +.10000000"+01  +.17763568"-14
      +.00000000"+00  +.10000000"+01

SOURCE TEXT(S):

"CODE" 34136;