A proposal for a cheap implementation for the Penrose-pseudoinverse

This is an attempt to implement a method along
the idea of Adam W. in *MSE, 27.2.2013*, using the notation in
the matrix-language/software MatMate.

*Gottfried Helms/Univ Kassel* 27.2.2013
.

We want to find the left-inverse for a *n x m* -
matrix, where *m<n*
. As an example we use a *6x4* random-matrix *A*

** A : **

–6.7475 |
–0.3647 |
–0.1494 |
0.0162 |

–3.4101 |
–7.0163 |
7.7587 |
–9.5665 |

–0.6796 |
7.4867 |
6.5455 |
1.8587 |

–5.0669 |
–4.2541 |
–9.5941 |
–9.8070 |

6.5135 |
5.4551 |
–7.1789 |
5.4896 |

–4.4194 |
9.5292 |
–7.1299 |
3.0132 |

If we would do this along the *SVD-decomposition*,
this would be costful because we needed iterative methods to find the left and
right diagonalizing rotation-matrices, which requires either two
eigen-decompositions or two iterative Jacobi-rotations towards orthogonality in
rows and in columns.

Conversely in this proposal it suffices to ** rotate
the rows only**, and also it suffices to rotate them

We need that rotation later so we append a *6x6*-unit-matrix
** unit(6)** to the right of

[19] tmp = rotate( (A || unit(6)) ' , "triangular") '

Note, that the applied rotation is saved by
this in the rightmost *6x6 *submatrix now.

Then we need only the *4* rows of it, call the resulting submatrix
** C** :

[20] C = tmp[1..4,*]

we see, that the rotated submatrix ** A** has become upper triangular in the left,
white area. The yellow area is the matrix

** C : **

12.0520 |
3.0097 |
0.2875 |
8.5780 |
–0.5599 |
–0.2829 |
–0.0564 |
–0.4204 |
0.5405 |
–0.3667 |

0.0 |
15.3302 |
–4.7314 |
10.1494 |
0.0861 |
–0.4021 |
0.4994 |
–0.1950 |
0.2497 |
0.6936 |

0.0 |
0.0 |
16.5835 |
1.0065 |
0.0253 |
0.3580 |
0.5382 |
–0.6269 |
–0.3710 |
–0.2257 |

0.0 |
0.0 |
0.0 |
7.2635 |
0.5396 |
–0.4706 |
–0.4500 |
–0.4944 |
–0.1800 |
–0.0900 |

After that we need only three
Gaussian-eliminations, beginning with the pivot *(4,4)* and involving the above rows:

[22] C = gausselim(C , {4,4},1..3)

**
****C : **

12.0520 |
3.0097 |
0.2875 |
0.0 |
–1.1971 |
0.2728 |
0.4750 |
0.1634 |
0.7531 |
–0.2604 |

0.0 |
15.3302 |
–4.7314 |
0.0 |
–0.6678 |
0.2555 |
1.1282 |
0.4959 |
0.5013 |
0.8193 |

0.0 |
0.0 |
16.5835 |
0.0 |
–0.0495 |
0.4232 |
0.6005 |
–0.5584 |
–0.3461 |
–0.2132 |

0.0 |
0.0 |
0.0 |
7.2635 |
0.5396 |
–0.4706 |
–0.4500 |
–0.4944 |
–0.1800 |
–0.0900 |

… now with pivot *(3,3)* and involving the above two
rows

[23] C = gausselim( C , {3,3},1..2)

**
****C : **

12.0520 |
3.0097 |
0.0 |
0.0 |
–1.1962 |
0.2655 |
0.4646 |
0.1731 |
0.7591 |
–0.2567 |

0.0 |
15.3302 |
0.0 |
0.0 |
–0.6819 |
0.3762 |
1.2995 |
0.3365 |
0.4026 |
0.7585 |

0.0 |
0.0 |
16.5835 |
0.0 |
–0.0495 |
0.4232 |
0.6005 |
–0.5584 |
–0.3461 |
–0.2132 |

0.0 |
0.0 |
0.0 |
7.2635 |
0.5396 |
–0.4706 |
–0.4500 |
–0.4944 |
–0.1800 |
–0.0900 |

Finally with pivot *(2,2)* and involving the remaining
upper row:

[24] C = gausselim( C , {2,2},1..1)

**
****C : **

12.0520 |
0.0 |
0.0 |
0.0 |
–1.0623 |
0.1916 |
0.2095 |
0.1070 |
0.6800 |
–0.4056 |

0.0 |
15.3302 |
0.0 |
0.0 |
–0.6819 |
0.3762 |
1.2995 |
0.3365 |
0.4026 |
0.7585 |

0.0 |
0.0 |
16.5835 |
0.0 |
–0.0495 |
0.4232 |
0.6005 |
–0.5584 |
–0.3461 |
–0.2132 |

0.0 |
0.0 |
0.0 |
7.2635 |
0.5396 |
–0.4706 |
–0.4500 |
–0.4944 |
–0.1800 |
–0.0900 |

We see, that the left white segement is diagonal now; to make it a unit-matrix we scale the rows by the inverses of the diagonal:

[25] D = 1 /# diag(C2) *# C // the #-sign at the mathematical operators indicates elementwise operation

** D : **

1.0000 |
0.0 |
0.0 |
0.0 |
–0.0881 |
0.0159 |
0.0174 |
0.0089 |
0.0564 |
–0.0337 |

0.0 |
1.0000 |
0.0 |
0.0 |
–0.0445 |
0.0245 |
0.0848 |
0.0220 |
0.0263 |
0.0495 |

0.0 |
0.0 |
1.0000 |
0.0 |
–0.0030 |
0.0255 |
0.0362 |
–0.0337 |
–0.0209 |
–0.0129 |

0.0 |
0.0 |
0.0 |
1.0000 |
0.0743 |
–0.0648 |
–0.0619 |
–0.0681 |
–0.0248 |
–0.0124 |

… and in the yellow area we have the
pseudoinverse of ** A**, let's call it

[26] V = D [*,5..10]

We check it by multiplication ** V*A** should be the unit matrix

[27] chk = V*A

** chk : **

1.0000 |
0.0 |
0.0 |
0.0 |

0.0 |
1.0000 |
0.0 |
0.0 |

0.0 |
0.0 |
1.0000 |
0.0 |

0.0 |
0.0 |
0.0 |
1.0000 |

email to the author: helms(at)uni-kassel.de