Matrix Pseudo-Inverse for the HP-41

Overview
 

 1°)  The Method of Gréville ( 2nd version )
 2°)  Ben-Israel and Cohen Algorithm
 3°)  An Iterative Method of order 3
 4°)  Rank of a Matrix
 

-If A is a nxm matrix ( n rows, m columns ), the pseudoinverse A⁺ , also called Moore-Penrose inverse, of the matrix A
  A⁺ is the unique mxn matrix ( m rows, n columns ) such that:  A A⁺ A = A ,  A⁺ A A⁺  = A⁺ ,  A A⁺ and  A⁺ A  symmetric.

-The following routines only deal with real matrices.
 
 

1°)  The Method of Gréville ( 2nd Version )
 

-A program that uses Gréville's method is already listed in "Matrix Operations for the HP-41" where the formulas are given too.
-Unlike the 1st program, this variant does not call "TM*" as a subroutine.
-For the rest, the calculations are almost identical.
 

Data Registers:           R00 thru R08: temp                    ( Registers Rbb to Ree  are to be initialized before executing "GRVL" )

                                  •  Rbb thru  •  Ree  the coefficients of the matrix A, in column order   bb > 07

         >>>>   When the program stops, the pseudo-inverse A⁺ has replaced A in Rbb to Ree

                                      Registers  Ree+1 , ... are used during the calculations
Flags: /
Subroutines:    "TRN"  "M*"  "MCO"  ( cf  "Matrix Operations for the HP-41" )
                          "MCO" may be replaced by the M-Code routine "LCO"  listed in "Miscellaneous Short Routines for the HP-41"

-Line 194 is a three-byte  GTO 03
 
 

 
    ( 282 bytes / SIZE bb+3n.m+m-1 )
 
 

  where  bbb.eeenn = control number of the matrix A   bbb > 008

Example:    The same matrix as in "Matrix Operations"

                 1  1  4  2                     R10  R13  R16  R19
      A  =    0  1  2  3         =          R11  R14  R17  R20      ( respectively )      control number =  10.02103
                 3  2  6  7                     R12  R15  R18  R21
 

   •   10.02103         XEQ "GRVL"  >>>>   10.02104                  --- execution time = 71s ---

 and we get:

                  R10   R14   R18              -21/112    -85/112   43/112
                  R11   R15   R19     =         7/112      23/112   -9/112         =   A⁺   ( with some roundoff errors of course ).
                  R12   R16   R20                49/112      1/112   -15/112
                  R13   R17   R21               -35/112    29/112    13/112

-In this example,  A A⁺ = Id₃

Notes:

-SIZE 049 is enough for this matrix

-Lines 65 and 128,  E-7  is used to identify tiny real numbers.
-Another choice may be better according to the matrix elements.

-The routine stops at line 87.
-Lines 83 to 87 actually overwrite the original matrix A.
-If you want to save A, replace these lines by

    RCL 01           XEQ "MCO"       XEQ "TRN"
    RCL 06           RCL 01               RTN
    INT                 INT
 

2°)  Ben-Israel and Cohen Algorithm
 

-Ben-Israel and Cohen algorithm is an iterative method:

    A_k+1 =  2 A_k -  A_k A A_k      with  A₀ = µ A^T          where   A^T = transposed of A   and   0 < µ <= 2 / Trace(A.A^T)

-The following program uses  µ = 1 / Trace(A.A^T)
 

Data Registers:           R00 thru R05: temp                    ( Registers Rbb to Ree  are to be initialized before executing "BISC" )

                                  •  Rbb thru  •  Ree  the coefficients of the matrix A, in column order   bb > 05

         >>>>   When the program stops, the pseudo-inverse A⁺ is stored in Ree+1 to Ree+m.n

Flags: /
Subroutines:    "TRN"  "M*"  ( cf  "Matrix Operations for the HP-41" )

-Lines 31 to 38 may be replaced by   XEQ "TRACE"   ( cf  "Matrix Operations for the HP-41" )
 
 

 
    ( 127 bytes / SIZE bb+3n.m+n^2 )
 
 

  where  bbb.eeenn = control number of the matrix A   bbb > 005

Example:       The same matrix again:

                 1  1  4  2                     R10  R13  R16  R19
      A  =    0  1  2  3         =          R11  R14  R17  R20      ( respectively )      control number =  10.02103
                 3  2  6  7                     R12  R15  R18  R21
 

   •   10.02103         XEQ "BISC"  >>>>   22.03304                  --- execution time = 7m40s ---

 and we get:

                  R22   R26   R30              -21/112    -85/112   43/112
                  R23   R27   R31     =         7/112      23/112   -9/112         =   A⁺   ( with some roundoff errors ).
                  R24   R28   R32                49/112      1/112   -15/112
                  R25   R29   R33               -35/112    29/112    13/112

Notes:

-"BISC" is obviously very slow but with a good emulator...
-The advantage is that it's also much shorter than "GRVL" !

-The successive "amplitudes" of the increments are successively displayed ( line 71 )
-Though very slow in the first iterations, the convergence becomes quadratic near the end.

-Line 72,  10^(-7) is used to evaluate the difference between 2 successive approximations.
-This real may be changed according to the type of matrix.
-Anyway, after the program stops, simply execute XEQ 02 to perform another iteration if need be.

   A_k+1 =  2 A_k -  A_k A A_k  may also be used to improve the precision of the result returned by "GRVL"
 

3°)  An Iterative Method of Order 3
 

-Other methods of higher order may also be derived:
 

Formula of order q:

    with  A₀ = µ A^T          where     0 < µ <= 2 / Trace(A.A^T)

    T_k = I - A_k A
    M_k = I + T_k + ..... + T_k^q-1
   A_k+1 = M_k A_k                             the sequence  ( A_k )  tends to  A⁺  when k tends to infinity
 

-The routine hereunder uses  q = 3
 
 

Data Registers:           R00 thru R06: temp                     ( Registers Rbb to Ree  are to be initialized before executing "INV3" )

                                  •  Rbb thru  •  Ree  the coefficients of the matrix A, in column order   bb > 06

         >>>>   When the program stops, the pseudo-inverse A⁺ is stored in Ree+1 to Ree+m.n

Flags: /
Subroutines:    "TRN"  "M*"  ( cf  "Matrix Operations for the HP-41" )

-Lines 33 to 40 may be replaced by   XEQ "TRACE"   ( cf  "Matrix Operations for the HP-41" )
 
 

 
    ( 151 bytes / SIZE bb+4n.m+n^2 )
 
 

  where  bbb.eeenn = control number of the matrix A   bbb > 006

Example:    Still the same matrix

                 1  1  4  2                     R10  R13  R16  R19
      A  =    0  1  2  3         =          R11  R14  R17  R20      ( respectively )      control number =  10.02103
                 3  2  6  7                     R12  R15  R18  R21
 

   •   10.02103         XEQ "INV3"  >>>>   22.03304                  --- execution time = 7m40s ---

 and we get:

                  R22   R26   R30              -21/112    -85/112   43/112
                  R23   R27   R31     =         7/112      23/112   -9/112         =   A⁺   ( with some roundoff errors ).
                  R24   R28   R32                49/112      1/112   -15/112
                  R25   R29   R33               -35/112    29/112    13/112

Notes:

-In this example, the execution is the same !
-Convergence is faster, so less iterations are needed, but each iteration is more time consuming.

-If you want to execute another loop, simply press  XEQ 02
 

4°)  Rank of a Matrix
 

-The rank r(A) of an nxm matrix A may be computed by    r(A) = Trace ( A A⁺ )
-So we can use one of the 3 programs above to calculate  r(A)

-"TRACE" is listed in "Matrix Operations for the HP-41"

-If you want to use the 1st version, on the left, modify "GRVL" to save the matrix A, as suggested in the note, paragraph 1.
-If you prefer the 2nd version, add  STO 05  after line 55 in the "BISC" listing.
-In the 3rd version on the right, "INV3" doesn't need any change.
 
 
 

 
        ( 37 bytes )                 ( 30 bytes )                ( 30 bytes )
 
 

  where  bbb.eeenn = control number of the matrix A

Example1:

                 1  1  4  2                     R10  R13  R16  R19
      A  =    0  1  2  3         =          R11  R14  R17  R20           control number =  10.02103
                 3  2  6  7                     R12  R15  R18  R21
 

   •   10.02103         XEQ "RANK"  >>>>   3.0000

Example2:

                 1  1  4  2                     R10  R13  R16  R19
      A  =    0  1  2  3         =          R11  R14  R17  R20           control number =  10.02103
                 1  3  8  8                     R12  R15  R18  R21
 

   •   10.02103         XEQ "RANK"  >>>>   2.0000
 

Notes:

-Due to roundoff-errors, "TRACE" doesn't always return an integer.
-That's why  FIX 0  RND  FIX 4  are added before the  END
-Alternatively, these 3 instructions may be replaced by  .5  +   INT

-The 1st version  - that calls "GRVL" - is of course much faster.

-But if you like the iterative method of Ben-Israel & Cohen, the following program can be used to get the rank only.
 

Formula:

    B₀   =  ( A A^T )  /   Trace ( A A^T )
   B_k+1 = 2 B_k - B_k²

  >>>     r(A) = Lim_k->infinity Tr ( B_k )
 
 

Data Registers:           R00 thru R04: temp                    ( Registers Rbb to Ree  are to be initialized before executing "RANK" )

                                  •  Rbb thru  •  Ree  the coefficients of the matrix A, in column order   bb > 04

         >>>>   When the program stops, the rank  r(A) is in X ( the non-rounded value is in R04 )

Flags: /
Subroutines:    "TRN"  "M*"  "MCO"
 
 

 
      ( 140 bytes )
 
 

  where  bbb.eeenn = control number of the matrix A

Example:

                 1  1  4  2                     R10  R13  R16  R19
      A  =    0  1  2  3         =          R11  R14  R17  R20           control number =  10.02103
                 1  3  8  8                     R12  R15  R18  R21
 

   •   FIX 7    10.02103         XEQ "RANK"  >>>>   r(A) = 2        ---Execution time = 3m22s---
 

Notes:

-The termination criterion depends on the display setting.
-The program stops when the rounded difference between 2 successive values = 0

-Lines 40 & 68 may be replaced by  XEQ "TRACE"
-In this case, delete lines 80 to 91.

-If there are more rows than columns, store A^T instead of A.

Warning:

-In the example above, if you execute "RANK" in FIX 9, it first seems to converge to 2
  but then, Tr ( B_k ) starts to decrease and finally tends to minus infinity !
-So, the formula becomes unstable, especially if the matrix is badly ill-conditionned.
 

References:

[1]  Adi Ben-Israel & Dan Cohen - "On Iterative Computation of Generalized Inverses and Associated Projections" -
          J. Siam Numer. Anal. Vol 3, n°3 , 1966
[2]  R. Penrose - "A Generalized Inverse for Matrices"
[3]  Ben-Israel & Greville - "Generalized Inverses: Theory and Applications" -