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

Bug 1017

Summary:

Tridiagonalization often generates non-unitarian Q-matrix for rank-deficient input

Product:

Eigen

Reporter:

Christoph Hertzberg <chtz>

Component:

Householder

Assignee:

Nobody <eigen.nobody>

Status:

DECISIONNEEDED ---

Severity:

Accuracy Problem

CC:

chtz, gael.guennebaud, zh2196

Priority:

Normal

Version:

3.3 (current stable)

Hardware:

All

OS:

All

Whiteboard:

Bug Depends on:

Bug Blocks:

661, 814

Attachments:

Description	Flags
Avoid underflow problems in makeHouseholder	none

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.

https://bitbucket.org/eigen/eigen/commits/8f80a9f1e9b5/
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. ***

Comment 8 Nobody 2019-12-04 14:39:01 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/1017.