I am trying to experiment with JacobiSVD of Eigen. In particular I am trying to reconstruct the input matrix from its singular value decomposition. http://eigen.tuxfamily.org/dox/classEigen_1_1JacobiSVD.html.
Eigen::MatrixXf m = Eigen::MatrixXf::Random(3,3);
Eigen::JacobiSVD<Eigen::MatrixXf, Eigen::NoQRPreconditioner> svd(m, Eigen::ComputeFullU | Eigen:: ComputeFullV);
Eigen::VectorXf SVec = svd.singularValues();
Eigen::MatrixXf S = Eigen::MatrixXf::Identity(3,3);
S(0,0) = SVec(0);
S(1,1) = SVec(1);
S(2,2) = SVec(2);
Eigen::MatrixXf recon = svd.matrixU() * S * svd.matrixV().transpose();
cout<< "diff : \n"<< m - recon << endl;
I know that internally the SVD is computed by an iterative method and can never get a perfect reconstruction. The errors are in order of 10^-7. With the above code the output is --
diff :
9.53674e-07 1.2517e-06 -2.98023e-07
-4.47035e-08 1.3113e-06 8.9407e-07
5.96046e-07 -9.53674e-07 -7.7486e-07
For my application this error is too high, I am aiming for an error in the range 10^-10 - 10^-12. My question is how to set the threshold for the decomposition.
NOTE : In the docs I have noted that there is a method setThreshold() but it clearly states that this does not set a threshold for the decomposition but for singular values comparison with zero.
NOTE : As far as possible I do not wish to go for double. Is it even possible with float?
A single precision floating point (a 32 bit float) has between six to nine significant decimal figures, so your requirement of 10^{-10} is impossible (assuming the values are around 0.5f). A double precision floating point (a 64 bit double) has 15-17 significant decimal figures, so should work as long as the values aren't 10^6.