JacobiSVD fails on small rank-deficient float matrices (U is not unitary)

Submitted by Charles Karney

Assigned to Nobody

Link to original bugzilla bug (#790)
Version: 3.2
Operating system: Linux

Description

(Found during testing of the new alignPoints function. This bug will
also cause umeyama to fail.)

JacobiSVD can fail spectacularly with small rank-deficient float
matrices. Here's an example (Eigen 3.2.1, Linux 64-bit, g++ (GCC)
4.7.2 20121109, optimize or debug):

Matrix<float,4,4> M;  
M <<   
  0.030851028,0.275531113,-0.451450854,-0.126374155,  
  -0.015783134,-0.140959471,0.230958566,0.064651995,  
  0.051904849,0.463563203,-0.759536683,-0.212616309,  
  0.037682864,0.336546361,-0.551422834,-0.154359206;  
JacobiSVD< Matrix<float,4,4> > svd(M, ComputeFullU | ComputeFullV);  
std::cout << "U\n" << svd.matrixU() << "\n\n";  
std::cout << "U.determinant() " << svd.matrixU().determinant() << "\n\n";  

This prints out

U  
   -0.423173    -0.740408 -1.07744e-11            0  
    0.216492     0.378787  5.51207e-12            0  
   -0.711961     0.555262  8.08013e-12            0  
   -0.516883            0            0            0  

U.determinant() 0  

JacobiSVD does OK on M.transpose() and if float is promoted to double.
I realize that there are limits to what I can expect from floats. But
it would be nice if the basic properties of the SVD were more accurately
met (specifically U should be nearly unitary).

Blocking

#558 (closed)