Bugzilla – Bug 1017
Tridiagonalization often generates non-unitarian Q-matrix for rank-deficient input
Last modified: 2016-07-11 08:01:19 UTC
This seems to be the cause of eigensolver_selfadjoint_3 failing most of the time.
As far as I see what happens is that (in src/Householder/Householder.h) makeHouseholder is not robust to underflows, i.e. if c0!=0 && c0*c0==0 then
beta = sqrt(numext::abs2(c0) + tailSqNorm);
is slightly too small (especially if tailSqNorm is quite small as well).
Created attachment 582 [details]
Avoid underflow problems in makeHouseholder
The attached patch resolves the issue (at least all eigensolver_selfadjoint pass with that). However, I'd prefer adding numeric_limits::min (which is the smallest normalized positive number) to our own NumTraits struct.
Perhaps, an alternative fix would be:
beta = numext::hypot(c0,sqrt(tailSqNorm));
hm, after all your suggestion is probably good enough as comparing tailSqNorm to tol would still be needed to avoid denormalized. Using hypot could save squaring imag(c0), but I'm not 100% confident with hypot on complexes (abs and div on complexes are not safe on some compiler) and if we use hypot here, then should also use the costly stableNorm in place of squaredNorm as in bug 661.
I think it is better to make sure that makeHouseholder is only called after the matrix has been normalized (i.e. divided by maxCoeff), which is what the SVD and Eigensolver modules already do.
Then only underflows may occur, and only if the columns are linearly dependent (up to precision). For these cases, using the identity reflection should be good enough. Actually, that's not a reflection, so probably tau=2 and beta=-real(c0) would be better (in case the number of reflections is important).
yes, I agree with your patch. Feel free to push it.
I don't see use cases for which having reflections is needed, so I applied it as is.
Summary: Bug 1017: apply Christoph's patch preventing underflows in makeHouseholder
We could still think about implementing scaling in *HouseholderQR, but have to find an elegant way to avoid applying scaling twice, for instance when QR is called from JacobiSVD.
*** Bug 1252 has been marked as a duplicate of this bug. ***