New user self-registration is currently disabled. Please email eigen-core-team @ if you need an account.
Bug 1017 - Tridiagonalization often generates non-unitarian Q-matrix for rank-deficient input
Tridiagonalization often generates non-unitarian Q-matrix for rank-deficient ...
Product: Eigen
Classification: Unclassified
Component: Householder
3.3 (current stable)
All All
: Normal Accuracy Problem
Assigned To: Nobody
Depends on:
Blocks: 661 3.4
  Show dependency treegraph
Reported: 2015-05-18 18:02 UTC by Christoph Hertzberg
Modified: 2016-07-11 08:01 UTC (History)
3 users (show)

Avoid underflow problems in makeHouseholder (558 bytes, patch)
2015-05-19 21:37 UTC, Christoph Hertzberg
no flags Details | Diff

Description Christoph Hertzberg 2015-05-18 18:02:35 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).
Comment 1 Christoph Hertzberg 2015-05-19 21:37:11 UTC
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.
Comment 2 Gael Guennebaud 2015-06-05 13:07:50 UTC
Perhaps, an alternative fix would be:

beta = numext::hypot(c0,sqrt(tailSqNorm));
Comment 3 Gael Guennebaud 2015-06-05 13:40:32 UTC
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.
Comment 4 Christoph Hertzberg 2015-06-06 15:28:45 UTC
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).
Comment 5 Gael Guennebaud 2015-06-16 21:08:44 UTC
yes, I agree with your patch. Feel free to push it.
Comment 6 Gael Guennebaud 2015-06-22 14:54:06 UTC
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.
Comment 7 Gael Guennebaud 2016-07-11 08:01:19 UTC
*** Bug 1252 has been marked as a duplicate of this bug. ***

Note You need to log in before you can comment on or make changes to this bug.