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 :

 

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 non-iteratively such that A simply gets into upper triangular shape.

We need that rotation later so we append a 6x6-unit-matrix unit(6) to the right of A. In MatMate we have only column-rotations so we need also to transpose that aggregate such that we have a 10 x 6 – matrix, which gets column-rotated and then backtransposed to 6 x 10 :

 [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 t which contains the applied rotation.

      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 V:

[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