> On Wed, 09 Mar 2011 00:28:01 -0500 > Chris Marshall <[email protected]> wrote: > > 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
If the only thing preventing PDL::LinearAlgebra from being included by default is a dependency on LAPACK, I really think it should be included. LAPACK is VERY standard and has plenty of mature and efficient implementations. Rewriting LAPACK in parallel to avoid this dependency really does not seem like the right thing to do. Is it difficult to package? Are there licensing problems? dima _______________________________________________ Perldl mailing list [email protected] http://mailman.jach.hawaii.edu/mailman/listinfo/perldl
