Pseudoinversing matrix .. is this correct?

Angharad has asked for the wisdom of the Perl Monks concerning the following question:

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;
[download]

Does this look OK to people or not?

Comment on Pseudoinversing matrix .. is this correct? Download Code

Replies are listed 'Best First'.
Re: Pseudoinversing matrix .. is this correct? by etj (Priest) on Jun 21, 2022 at 19:56 UTC
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.	[reply]
Re: Pseudoinversing matrix .. is this correct? by wufnik (Friar) on Aug 03, 2005 at 14:52 UTC
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` [download] results courtesy matlab. ...wufnik -- in the world of the mules there are no rules --	[reply] [d/l]
Re^2: Pseudoinversing matrix .. is this correct? by Angharad (Pilgrim) on Aug 03, 2005 at 15:02 UTC
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?	[reply]
Re: Pseudoinversing matrix .. is this correct? by tall_man (Parson) on Aug 11, 2005 at 23:18 UTC
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;` [download] 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.	[reply] [d/l] [select]