1252
2016-07-11 02:24:14 +0000
JacobiSVD gives NaN
2019-12-04 15:58:49 +0000
1
1
1
Unclassified
Eigen
SVD
unspecified
All
Linux
RESOLVED
FIXED
Normal
Crash
---
0
zh2196
eigen.nobody
chtz
gael.guennebaud
oldest_to_newest
5969
0
zh2196
2016-07-11 02:24:14 +0000
Here's the code:
#############################################################
#include <iostream>
#include "../eigen-3.2.8/Eigen/Dense"
using namespace Eigen;
using namespace std;
#define cp complex<double>
int main(){
MatrixXcd m; m.resize(6,8);
m<<cp(0,-3.97691e-58), 0, 0, 0, 0, cp(0.389209,0.000989259), cp(0.766077,-3.17912e-19), 0,
0, cp(3.95428e-17,-2.70056e-20), 0, 0, 0, cp(0,-8.30225e-91), 0, cp(-0.511509,5.44428e-19),
0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0,
0, cp(-2.3115e-92,-2.55599e-74), 0, 0, 0, cp(-1.62153e-93,8.30225e-91), 0, cp(9.82409e-18,-6.77259e-35),
0, 0, 0, 0, 0, 0, 0, cp(-5.26899e-94,8.30225e-91);
cout<<m<<endl;
JacobiSVD<MatrixXcd> svd(m, ComputeThinU | ComputeThinV);
cout<<svd.matrixU()<<endl<<endl;
cout<<svd.singularValues().transpose()<<endl<<endl;
cout<<svd.matrixV()<<endl<<endl;
return 0;
}
#############################################################
The matrixV is an 8x6 matrix of NaNs. The version of eigen was 3.2.8. The 2 big singular values agree with Matlab.
5970
1
gael.guennebaud
2016-07-11 08:01:19 +0000
Same issue as in bug 1017, I've backported the respective change to 3.2.
*** This bug has been marked as a duplicate of bug 1017 ***
5972
2
716
zh2196
2016-07-11 16:46:03 +0000
Created attachment 716
data.bin (test case)
The binary input file of the 8x8 complex matrix m. File size = 8*8*16=1024 bytes.
5973
3
zh2196
2016-07-11 16:47:52 +0000
I used the patch in bug 1017. However, I found another example for which the result is correct, but a floating-point exception was still encountered during the implementation of JacobiSVD. This was tracked using fenv.h. Here's the code:
#############################################################
#include <iostream>
#include <fstream>
#include <fenv.h>
#include "../eigen-3.2.8/Eigen/Dense"
using namespace Eigen;
using namespace std;
#define cp complex<double>
int main(){
ifstream fp("data.bin",ios::in | ios::binary); int x,y;
MatrixXcd m(8,8); feenableexcept(FE_INVALID | FE_OVERFLOW);
/* //The single-precision approximation of m:
m<<0, cp(-2.82726e-51,-8.24618e-52), 0, 0, 0, 0, cp(0.63724,0), cp(8.94661e-36,6.17364e-37),
0, 0, cp(-0.213794,0), cp(-3.00159e-36,-2.07126e-37), 0, 0, 0, cp(0.619902,-0.250179),
0, 0, 0, 0, cp(-0.213801,0), cp(-3.00169e-36,-2.07133e-37), 0, 0,
cp(0.0720597,0), cp(1.01169e-36,6.98121e-38), 0, 0, 0, cp(-0.208243,0.084121), 0, 0,
0, cp(-1.25496e-18,0), cp(-0.00163295,-0.000394592), cp(-2.25437e-38,-7.12194e-39), 0, 0, 0, cp(-0.000647887,0.000174169),
cp(3.33671e-17,-7.75707e-35), cp(-2.00658e-19,-6.94573e-20), cp(5.32267e-35,3.4142e-35), cp(1.58459e-17,1.11847e-17), 0, cp(1.17701e-17,-6.42684e-18), 0, cp(2.51903e-35,1.78191e-36),
cp(-1.15011e-17,2.71009e-31), cp(1.08395e-19,-1.56512e-19), cp(-4.23255e-35,3.77238e-35), cp(4.06081e-17,2.86604e-17), 0, cp(-4.05699e-18,2.21523e-18), 0, cp(-7.61704e-36,2.13223e-35),
cp(4.75655e-20,1.20369e-34), cp(1.85949e-17,6.41572e-19), cp(-5.33786e-33,-1.48642e-33), cp(2.56783e-20,-5.52577e-20), 0, cp(1.67786e-20,-9.16159e-21), 0, cp(-2.15546e-33,5.00435e-34),*/
for(x=0;x<8;x++)
for(y=0;y<8;y++)
fp.read((char*)&(m(x,y)),16);
fp.close(); cout<<m<<endl<<endl;
cout<<"Before SVD"<<endl;
JacobiSVD<MatrixXcd> svd(m, ComputeThinU | ComputeThinV);
cout<<"After SVD"<<endl;
cout<<svd.matrixU()<<endl<<endl;
cout<<svd.singularValues().transpose()<<endl<<endl;
cout<<svd.matrixV()<<endl; return 0;
}
#############################################################
This hidden bug can only be reproduced using the exact double-precision input from the file data.bin, which was uploaded just now as an attachment.
5974
4
gael.guennebaud
2016-07-11 20:29:39 +0000
can you reproduce with devel branch? (there are a few other changes in JacobiSVD that I have not backported yet)
5975
5
gael.guennebaud
2016-07-11 20:46:05 +0000
ok, I've tried with the devel branch, and no overflow/underflow occurs. So I'll just backport all other changes.
5976
6
gael.guennebaud
2016-07-11 20:48:40 +0000
https://bitbucket.org/eigen/eigen/commits/75abaeb43845/
Branch: 3.2
Summary: Backport numerical robustness fixes from 3.3 branch
5977
7
zh2196
2016-07-11 22:44:40 +0000
Can you tell me how I can track floating-point exceptions using devel branch? Is it possible for fenv.h to give false errors (terminating the code with "core dumped" message when there's really no problem)? I used the more robust version of JacobiSVD.h (with 1000 lines) you offered but fenv.h still gives the error, while the result is correct if fenv.h is disabled.
5978
8
zh2196
2016-07-11 22:58:16 +0000
Sorry, my fault. I forgot to change the version of Eigen to the modified one. Now it's all right. Both test cases are passed with fenv.h switched on. This 1000-line version of SVD code is now more robust. Thanks for your work!
9696
9
eigen.nobody
2019-12-04 15:58:49 +0000
-- GitLab Migration Automatic Message --
This bug has been migrated to gitlab.com's GitLab instance and has been closed from further activity.
You can subscribe and participate further through the new bug through this link to our GitLab instance: https://gitlab.com/libeigen/eigen/issues/1252.
716
2016-07-11 16:46:03 +0000
2016-07-11 16:46:03 +0000
data.bin (test case)
data.bin
application/octet-stream
1024
zh2196
AAAAAAAAAAAAAAAAAAAAALcTCIDD7HC1q8FeFeS+U7UAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA
AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAASYSwlkVk5D8AAAAAAAAAAMq4
O5a/yKc4kivv3oBCajgAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAMIoWuSdXcu/AAAA
AAAAAABjuGVCHeuPuC1wnLfOnlG4AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA
AAAAAAAAAAAAAAAAugesWj3W4z+uwFVN8ALQvwAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA
AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAnI7Ol1l3LvwAAAAAAAAAA8MXvdF/r
j7g95dRC855RuAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAaRz4PIFysj8AAAAAAAAA
AAG4ungrhHU4tciZnH7BNzgAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA
AAAAAAAAAABycPCIs6fKv9qGIIL0iLU/AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA
AAAAAAAAAAAAAAAAAAAAHaZD5l0mN7wAAAAAAAAAALR1sNcTwVq/2RPZISncOb9mDj5MU68euNRi
se1EYwO4AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAUNsj
Od46Rb9S5/9HKNQmP1vUTyUdPIM8jMi2iPvG2bgPxXUMppwNvAN/4MYKgPS7K63JOwew0TgBe3YL
9rDGOLRMReLlRHI8/xQX8kXKaTwAAAAAAAAAAAAAAAAAAAAAiIW8i94jazz/w9f1eaNdvAAAAAAA
AAAAAAAAAAAAAAA0APVN6r3AOC+RxerS8oI4P49LqRCFarzsZJD5ofyVOREnoikg/v87oOpStdcY
B7wHS15xUSHMuLgkDwFeEsk4vn4QULRohzw4op53hYWAPAAAAAAAAAAAAAAAAAAAAAC+qQG0pbVS
vKERvZ2SbkQ8AAAAAAAAAAAAAAAAAAAAAGXaUubiP6S4I7jnRZlXvDhCAKMI5BPsO/3D/dby/+M4
GF2oMT9wdTwkRsTneKsnPFVq/GY6tzu5o9dJjDDfHrncESjHyVDeOw29GrMnT/C7AAAAAAAAAAAA
AAAAAAAAALTToTwEz9M77NBbI9ihxbsAAAAAAAAAAAAAAAAAAAAA/izQWjBiJrkbKwEKickEOQ==