Pseudo-Inverse

# 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

 01  LBL "GRVL"   02  STO 00   03  STO 01    04  FRC   05  ISG X   06  INT   07  .1   08  %   09  ST+ 01   10  +   11  STO 02    12  LASTX   13  RCL 00   14  INT   15  .1   16  %   17  +   18  +   19   E3   20  STO 08    21  1/X   22  -   23  STO 03   24  STO 04           25  RCL 00   26  FRC   27  RCL 08   28  *   29  INT   30  1   31  +   32  STO 06   33  RCL 00 34  INT   35  -   36  ST+ 01   37  STO 05    38  ST+ 06   39  RCL 02   40  INT   41  ST/ 05   42   E5   43  /   44  ST+ 04   45  RCL 05    46  RCL 06   47  +   48  STO 07   49  RCL 03   50  RCL 01    51  INT   52  XEQ "MCO"   53  ENTER^   54  ENTER^   55  CLX   56  STO 05   57  DSE 07   58  LBL 01   59  RCL IND Y   60  X^2   61  +   62  ISG Y   63  GTO 01   64  SQRT   65   E-7   66  X>Y? 67  GTO 03   68  R^   69  LASTX   70  LBL 02   71  ST/ IND Y   72  ISG Y   73  GTO 02   74  LBL 03   75  RCL 02    76  ST+ 03   77  RCL 01    78  INT   79  RCL 03   80  INT   81  X#Y?   82  GTO 03   83  RCL 01    84  RCL 00   85  INT   86  XEQ "TRN"   87  RTN   88  LBL 03   89  LASTX   90  RCL 01   91  RCL 06           92  INT   93  XEQ "M*"   94  STO 06   95  RCL 05   96  +   97  RCL 04   98  X<>Y   99  RCL 07 100  INT 101  XEQ "M*" 102  STO 07 103  RCL 08 104  * 105  INT 106  RCL 08  107  / 108  STO M 109  RCL 03  110  LBL 04 111  RCL IND X 112  RCL IND Z 113  - 114  STO IND Z 115  ISG Y 116  RDN 117  ISG Y 118  GTO 04 119  RCL M 120  0 121  LBL 05 122  RCL IND Y 123  X^2 124  + 125  ISG Y 126  GTO 05        127  SQRT 128   E-7 129  X

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

 STACK INPUT OUTPUT X bbb.eeennA bbb.eeemmA+

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+ = Id3

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:

Ak+1 =  2 Ak -  Ak A Ak      with  A0 = µ AT          where   AT = transposed of A   and   0 < µ <= 2 / Trace(A.AT)

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

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

 01  LBL "BISC"  02  ENTER^  03  STO 00   04  FRC  05   E3  06  STO 05   07  *  08  INT  09  1  10  +  11  XEQ "TRN"  12  STO 01          13  RCL 00  14  -  15  INT  16  RCL 01   17  INT 18  +  19  STO 02  20  RCL 00   21  FRC  22  ISG X  23  INT  24  X^2  25  +  26  STO 03  27  RCL 00           28  RCL 01   29  RCL 02  30  XEQ "M*"  31   E-5  32  +  33  0  34  LBL 00 35  RCL IND Y  36  +  37  ISG Y  38  GTO 00  39  RCL 01  40  RCL 05   41  *  42  INT  43  RCL 05           44  /  45  STO 04   46  X<>Y  47  LBL 01  48  ST/ IND Y  49  ISG Y  50  GTO 01  51  LBL 02 52  RCL 00  53  RCL 01  54  RCL 02  55  XEQ "M*"  56  RCL 01   57  X<>Y  58  RCL 03  59  XEQ "M*"  60  INT  61  STO M  62  RCL 04          63  STO N  64  CLX  65  LBL 03  66  RCL IND N  67  RCL IND M  68  - 69  ST+ IND N  70  ABS  71  +  72  ISG M  73  CLX  74  ISG N  75  GTO 03  76  VIEW X  77   E-7  78  X

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

 STACK INPUT OUTPUT X bbb.eeennA bbb.eeemm'A+

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.

Ak+1 =  2 Ak -  Ak A Ak  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  A0 = µ AT          where     0 < µ <= 2 / Trace(A.AT)

Tk = I - Ak A
Mk = I + Tk + ..... + Tkq-1
Ak+1 = Mk Ak                             the sequence  ( Ak )  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" )

 01  LBL "INV3"  02  ENTER^  03  STO 00          04  FRC  05   E3  06  STO 06   07  *  08  INT  09  1  10  +  11  XEQ "TRN"  12  STO 01  13  RCL 00   14  -  15  INT  16  STO 04  17  RCL 01  18  INT  19  +  20  STO 02 21  RCL 00          22  FRC  23  ISG X  24  INT  25  X^2  26  +  27  STO 03   28  ST+ 04  29  RCL 00  30  RCL 01  31  RCL 02   32  XEQ "M*"  33   E-5  34  +  35  0  36  LBL 00  37  RCL IND Y  38  +  39  ISG Y  40  GTO 00 41  RCL 01  42  RCL 06  43  *  44  INT  45  RCL 06          46  /  47  STO 05   48  X<>Y  49  LBL 01  50  ST/ IND Y  51  ISG Y  52  GTO 01  53  LBL 02  54  RCL 00   55  RCL 01  56  RCL 02  57  INT  58  XEQ "M*"  59  STO 02  60  RCL 01 61  X<>Y  62  RCL 03  63  XEQ "M*"  64  RCL 02   65  RCL 04          66  XEQ "M*"  67  INT  68  STO M  69  RCL 05   70  STO N  71  RCL 03   72  STO O  73  CLX  74  LBL 03  75  RCL IND N  76  ST+ X  77  RCL IND M  78  +  79  RCL IND O  80  3 81  *  82  -  83  ST+ IND N  84  ABS  85  +  86  ISG M  87  CLX  88  ISG O          89  CLX  90  ISG N  91  GTO 03          92  VIEW X  93   E-7  94  X

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

 STACK INPUT OUTPUT X bbb.eeennA bbb.eeemm'A+

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.

 01  LBL "RANK"  02  XEQ "GRVL"  03  RCL 00  04  X<>Y  05  RCL 06  06  INT  07  XEQ "M*"  08  XEQ "TRACE"  09  FIX 0  10  RND  11  FIX 4  12  END 01  LBL "RANK"  02  XEQ "BISC"  03  RCL 05  04  XEQ "TRACE"  05  FIX 0  06  RND  07  FIX 4  08  END 01  LBL "RANK"  02  XEQ "INV3"  03  RCL 02  04  XEQ "TRACE"  05  FIX 0  06  RND  07  FIX 4  08  END

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

 STACK INPUT OUTPUT X bbb.eeennA r(A)

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:

B0   =  ( A AT )  /   Trace ( A AT )
Bk+1 = 2 Bk - Bk2

>>>     r(A) = Limk->infinity Tr ( Bk )

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"

 01  LBL "RANK"  02  ENTER^  03  STO 00          04  FRC  05   E3  06  STO 04   07  *  08  INT  09  1  10  +  11  XEQ "TRN"  12  STO 01   13  ST+ X  14  RCL 00  15  -  16  INT  17   E-5  18  STO 03  19  RCL 00 20  RCL 01  21  R^  22  XEQ "M*"  23  RCL 00          24  INT  25  XEQ "MCO"  26  STO 00  27  FRC  28  RCL 04   29  *  30  INT  31  STO 01   32  RCL 00  33  RCL 04  34  *  35  INT  36  RCL 04  37  /  38  STO 02 39  RCL 00          40  XEQ 00  41  RCL 02  42  X<>Y  43  LBL 01  44  ST/ IND Y  45  ISG Y  46  GTO 01  47  CHS  48  STO 04   49  ISG 01  50  LBL 10  51  RCL 00   52  RCL 00  53  RCL 01  54  XEQ "M*"  55  RCL 01  56  RCL 02  57  LBL 02 58  RCL IND X  59  RCL IND Z  60  -  61  ST+ IND Y  62  RDN  63  ISG Y  64  CLX  65  ISG X  66  GTO 02  67  RCL 00          68  XEQ 00  69  VIEW X  70  ENTER^  71  X<> 04  72  -   73  RND  74  X#0?  75  GTO 10  76  RCL 04 77  .5  78  +  79  INT  80  RTN  81  LBL 00  82  RCL 03          83  +   84  ENTER^  85  CLX  86  LBL 03  87  RCL IND Y  88  +  89  ISG Y  90  GTO 03  91  RTN  92  END

( 140 bytes )

 STACK INPUT OUTPUT X bbb.eeennA r(A)

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 AT 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 ( Bk ) 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" -