SandMatrix Routine of the week: Pseudo-Inverse

[Ángel Martin]

Even if the convergence is finicky, using the BenIsrael-Cohen algorithm it takes a relatively short program to calculate the Moore-Penrose matrix pseudo-inverse. Once the initial estimation is calculated it's a matter of a few iterations until the final result is obtained - although each iteration also takes its time.

For additional references see:
http://benisrael.net/COHEN-BI-ITER-GI.pdf

and JM Baillard's page:
http://hp41programs.yolasite.com/pseudoinverse.php

The code is listed below; lots of Alpha stuff as they are the input control strings that drive the matrix operations. Also note three auxiliary matrices are needed - mostly due to the M*M instruction that requires a result matrix different from the operands.

Code:


01    LBL "MPINV"
02    ASTO O1
03    "|-,#"
04    MAT=
05    DIM?
06    INT
07    I<>J
08    LASTX
09    +
10    "P"
11    MATDIM
12    "#"
13    TRNPS
14    XEQ 00
15    "P"
16    MTRACE
17    "#,X,#"
18    MAT/
19    FNRM
20    STO 00
21    DIM?
22    "Q"
23    MATDIM
24    LBL 01
25    XEQ 00
26    "#,P,Q"
27    M*M
28    "#"
29    MAT+
30    "#,Q,#"
31    MAT-
32    FNRM
33    X<> 00
34    RCL- (00)
35    ABS
36    VIEW X
37    E-7
38    X<=Y?
39    GTO 01
40    "P"
41    PURFL
42    "Q"
43    PURFL
44    MNAME?
45    RTN
46    LBL 00
47    CLA
48    ARCL 01
49    "|-,#,P"
50    M*M
51    END

You'll need the SandMatrix for functions I<>J, MTRACE and MAT=

Auxiliary matrices "P" and "Q" are purged upon completion.
The result matrix "#" will not overwrite the initial one.

Cheers,
'AM