Hi
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?
dima
_______________________________________________
Perldl mailing list
[email protected]
http://mailman.jach.hawaii.edu/mailman/listinfo/perldl
