This bugzilla service is closed. All entries have been migrated to https://gitlab.com/libeigen/eigen

Bug 1252

Summary: JacobiSVD gives NaN
Product: Eigen Reporter: zh2196 <zh2196>
Component: SVDAssignee: Nobody <eigen.nobody>
Status: RESOLVED FIXED    
Severity: Crash CC: chtz, gael.guennebaud
Priority: Normal    
Version: unspecified   
Hardware: All   
OS: Linux   
Whiteboard:
Attachments:
Description Flags
data.bin (test case) none

Description zh2196 2016-07-11 02:24:14 UTC
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.
Comment 1 Gael Guennebaud 2016-07-11 08:01:19 UTC
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 ***
Comment 2 zh2196 2016-07-11 16:46:03 UTC
Created attachment 716 [details]
data.bin (test case)

The binary input file of the 8x8 complex matrix m. File size = 8*8*16=1024 bytes.
Comment 3 zh2196 2016-07-11 16:47:52 UTC
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.
Comment 4 Gael Guennebaud 2016-07-11 20:29:39 UTC
can you reproduce with devel branch? (there are a few other changes in JacobiSVD that I have not backported yet)
Comment 5 Gael Guennebaud 2016-07-11 20:46:05 UTC
ok, I've tried with the devel branch, and no overflow/underflow occurs. So I'll just backport all other changes.
Comment 6 Gael Guennebaud 2016-07-11 20:48:40 UTC
https://bitbucket.org/eigen/eigen/commits/75abaeb43845/
Branch:      3.2
Summary:     Backport numerical robustness fixes from 3.3 branch
Comment 7 zh2196 2016-07-11 22:44:40 UTC
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.
Comment 8 zh2196 2016-07-11 22:58:16 UTC
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!
Comment 9 Nobody 2019-12-04 15:58:49 UTC
-- 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.