Is BDCSVD numerical stable?
Page 1 of 1 (6 posts)
Tags:
None
Registered Member 
|
I have been applying BDCSVD recently, and I've found that for some matrices, it will return a ZERO singular value, which should not be possible.
For example, I start with a 128*128 tridiagonal matrix with all the non-zero elements being k. (Denote it as T[k]). Then I compute the matrix exponential of T[k]: exp(T[k]), which is the easiest way to generate some very ill matrices. After that, I compute the SVD of exp(T[k]). I found that when k=-9 the smallest singular value will be zero, which I think will not be possible. For k=-8, however, the smallest singular value us 3.75658e-11. This is not a very small value. Some clarifications: 1. One may wonder that this may because the matrix exponential is not numerical stable. But I can start with exp(T[-1]), which is fairly well-behaved. Then I use a numerical stable way to compute exp(T[-1])exp(T[-1])exp(T[-1])...exp(T[-1]), (say, do a SVD every matrix multiplication and always leave the big or small numbers in the singular values. ) However, this gives essentially the same result. 2. I am using MatrixXd. So I don't think 3.75658e-11 is some small value that may cause a overflow. I write exactly the same code in MATLAB using their expm() to compute matrix exponential and svd() to compute SVD. For k=-8, the smallest singular value is 3.7950e-11, and for k=-9 it gives 1.7450e-12. Questions: I don't know what is going wrong with my test. The only thing I could guess is that BDCSVD is not numerical stable? Codes I used to test the BDCSVD:
|
Moderator
|
You must always compare to the largest singular values. 1e-11 can extremely huge or extremely small depending on the context.
So let s_0 be the largest one, then any singular values s_i < epsilon * s_0 can be considered as numerically zero. With double precision, epsilon is about 2.22e-16. Since for k=-9, s_0=8059.95, and if we take MAtlab's s_127=1.7450e-12, then we have: s_127/s_0=2.16503e-16 that is indeed numerically zero. BDCSVD has to deflate such numerically small singular-value to be numerically stable. If you want to preserve them then you need another algorithm, for instance here Eigen::JacobiSVD returns for k=-9 , s_127=4.64366e-12. |
|
Registered Member
|
I'm applying the SVD implemented in Eigen 3.3.4 to 66x66 square matrices with no particular features except they are not singular. If I use the Jacobi SVD I have no problems, but I wanna adopt the BDCSVD decomposition to get a faster result. Unfortunately I'm obtaining this floating point error:
Eigen::BDCSVD<Eigen::Matrix<double, -1, -1, 0, -1, -1> >::secularEq(double, Eigen::Ref<Eigen::Array<double, -1, 1, 0, -1, 1>, 0, Eigen::InnerStride<1> > const&, Eigen::Ref<Eigen::Array<double, -1, 1, 0, -1, 1>, 0, Eigen::InnerStride<1> > const&, Eigen::Ref<Eigen::Array<long, 1, -1, 1, 1, -1>, 0, Eigen::InnerStride<1> > const&, Eigen::Ref<Eigen::Array<double, -1, 1, 0, -1, 1>, 0, Eigen::InnerStride<1> > const&, double) Floating point exception Are there some matrices for which the BDCSVD algorithm becomes numerical unstable? |
|
Moderator
|
I've recently fixed a similar issue, so check with the head of the 3.3 branch (will be 3.3.5 soon): http://bitbucket.org/eigen/eigen/get/3.3.tar.gz
If that does not work, then please share your input matrix with full accuracy. |
|
Registered Member
|
I have followed your suggestion, but I am still experiencing the following floating point error:
Eigen::BDCSVD<Eigen::Matrix<double, -1, -1, 0, -1, -1> >::computeSingVals(Eigen::Ref<Eigen::Array<double, -1, 1, 0, -1, 1>, 0, Eigen::InnerStride<1> > const&, Eigen::Ref<Eigen::Array<double, -1, 1, 0, -1, 1>, 0, Eigen::InnerStride<1> > const&, Eigen::Ref<Eigen::Array<long, 1, -1, 1, 1, -1>, 0, Eigen::InnerStride<1> > const&, Eigen::Matrix<double, -1, 1, 0, -1, 1>&, Eigen::Ref<Eigen::Array<double, -1, 1, 0, -1, 1>, 0, Eigen::InnerStride<1> >, Eigen::Ref<Eigen::Array<double, -1, 1, 0, -1, 1>, 0, Eigen::InnerStride<1> >) at ??:? Floating point exception I have attached my input matrix as requested by you.
|
|
Moderator
|
works for me:
|
Page 1 of 1 (6 posts)
Bookmarks
Who is online
Registered users: Bing [Bot], daret, Google [Bot], sandyvee, Sogou [Bot]
