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<Y? 130 GTO 03 131 RCL 01 132 RCL 06 |
133 RCL 05
134 + 135 RCL 07 136 INT 137 XEQ "M*" 138 STO 07 139 RCL 08 140 * 141 INT 142 RCL 08 143 / 144 STO M 145 RCL 06 146 1 147 GTO 05 148 LBL 03 149 RCL M 150 LASTX 151 LBL 06 152 ST/ IND Y 153 ISG Y 154 GTO 06 155 RCL M 156 RCL 01 157 FRC 158 RCL 08 159 * 160 INT 161 1 162 + 163 XEQ "MCO" 164 RCL 07 165 RCL 06 |
166 RCL 07
167 FRC 168 RCL 08 169 * 170 INT 171 1 172 + 173 XEQ "M*" 174 RCL 08 175 * 176 INT 177 RCL 08 178 / 179 RCL 01 180 LBL 07 181 RCL IND Y 182 ST- IND Y 183 CLX 184 SIGN 185 + 186 ISG Y 187 GTO 07 188 RCL 02 189 FRC 190 ST+ 01 191 ST+ 04 192 E-5 193 ST+ 05 194 GTO 03 195 END |
( 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<Y? 79 GTO 02 80 CLA 81 CLD 82 RCL 01 83 END |
( 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<Y? 95 GTO 02 96 CLA 97 CLD 98 RCL 01 99 END |
( 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" -