On 3/6/2011 1:44 AM, Chris Marshall wrote:
> On 3/6/2011 12:20 AM, Dima Kogan wrote:
>>
>> It seems that the current implementation of PDL::MatrixOps::svd has a
>> bug where rank-deficient matrices can produce a non-unitary singular
>> vector matrix:
>>
>> pdl> p $x
>>
>> [
>> [1 2 3]
>> [4 5 6]
>> [4 5 6]
>> ]
>>
>> pdl> ($u,$s,$v) = svd($x)
>>
>> pdl> p $u
>>
>> [
>> [ 0.28303265 0.95911028 -4.0699249e-10]
>> [ 0.67819338 -0.2001343 -9.7522246e-10]
>> [ 0.67819338 -0.2001343 -9.7522246e-10]
>> ]
>>
>> pdl> p $s
>> [ 12.936563 0.80332812 3.9580815e-19]
>> pdl> p $v
>>
>> [
>> [ 0.44127483 -0.79913069 0.40824829]
>> [ 0.56800242 -0.10347264 -0.81649658]
>> [ 0.69473 0.59218541 0.40824829]
>> ]
>>
>>
>> We see that $u has a column of 0s, which is wrong. Here
>>
>> $x == $u x stretcher($s) x $v->transpose
>>
>> is still true despite the wrong $u because the associated singular
>> value is 0 (as it should be). PDL has 2 other implementations of the
>> SVD that do NOT have this issue: PDL::LinearAlgebra::{msvd,mdsvd}.
>> These seem to call LAPACK directly. Is this a bug worth fixing? Should
>> the non-LAPACK SVD method simply be deprecated and removed? Is there a
>> reason multiple SVD implementation exist in parallel?
Ideally, there would be common code for all PDL installations.
> I don't appear to have PDL::LinearAlgebra in my PDL
> version 2.4.7_010 so I'm, not sure to what you are
> referring as far as SVD implementations.
I was able to install PDL::LinearAlgebra. There are
a number of useful routines. The catch is that it
requires the availability of the LAPACK library to
link to.
> Bug fixes are preferred unless there is an option
> that is available for *all* PDL users. I see no
> documentation for PDL::LinearAlgebra is the latest
> PDL distribution so a working svd decomp in PDL::MatrixOps
> would be of use to me....
One possibility is to move to GSL based functionality
for the core PDL numeric components. It has the
advantage of being ANSI C based as opposed to F77
which makes things a bit easier as far as binding
to perl/PDL.
--Chris
_______________________________________________
Perldl mailing list
[email protected]
http://mailman.jach.hawaii.edu/mailman/listinfo/perldl
