Pseudoinversing matrix .. is this correct?

[Angharad]

Hi Monks!
Further to a post I made a day or so ago regarding the pseudoinverse of a mxn matrix, I proceeded to create the following script .. attempting to implement the equation pseudoinvB = (B^tB)^-1B^t where B^t = transpose of non square matrix B.

#!/usr/bin/perl 

use PDL;

my $matrix = pdl [
           [4, 2, 9],
           [2, 5, 5],
           [6, 3, 5],
           [8, 4, 6],
           ];

($a, $b, $c) = svd($matrix);

$bt = transpose($b);

$B1 = $bt*$b;

$B2 = $B1^-1;

$pseudoB = $B2*$bt;

Does this look OK to people or not?

[wufnik]

methodology aok, for the pseudoinverse of b.

check against pseudoinverse(b) ==



 
0.0565 0.0000 0.0000 0.0000
0.0000 0.2405 0.0000 0.0000
0.0000 0.0000 0.3122 0.0000

results courtesy matlab.

[Angharad]

Thanks a lot. Much appreciated. I wonder though why I dont get this answer. lol.
I am right in thinking though that I can calculate the pseudoinverse using SVD too arent I?

[tall_man]

I see a couple of problems here. First, there seems to be confusion between the mxn matrix B (which you assigned to the variable $matrix) and the diagonal matrix in the decomposition (which you assigned to the variable $b).

Secondly, the diagonal result in $b is returned in the form of a single row. It must be converted to a diagonal matrix before you can work on it.

I would do something like this instead:

my $mt = transpose($matrix);
my $m2 = $mt x $matrix;
my $m2i = $m2->inv;
my $psuedo = $m2i x $mt;

You will only need svd if the m2 matrix is singular or ill-conditioned. It will allow you to "fix" any near-zeros on the diagonal.

[etj]

To do this the very easy way, use https://metacpan.org/pod/PDL::LinearAlgebra#mpinv which wraps fast LAPACK functions for PDL, in this case for the pseudo-inverse.