http://qs1969.pair.com?node_id=480466

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;
Does this look OK to people or not?

Replies are listed 'Best First'.
Re: Pseudoinversing matrix .. is this correct?
by etj (Chaplain) on Jun 21, 2022 at 19:56 UTC
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
    results courtesy matlab.
    ...wufnik

    -- in the world of the mules there are no rules --
      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?
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;

    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.