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
