Re: [eigen] Issue in SelfAdjointEigenSolver for 3-by-3 matrices |

[ Thread Index | Date Index | More lists.tuxfamily.org/eigen Archives ]

*To*: eigen@xxxxxxxxxxxxxxxxxxx*Subject*: Re: [eigen] Issue in SelfAdjointEigenSolver for 3-by-3 matrices*From*: Benoit Jacob <jacob.benoit.1@xxxxxxxxx>*Date*: Thu, 15 Jul 2010 16:10:21 -0400*Dkim-signature*: v=1; a=rsa-sha256; c=relaxed/relaxed; d=gmail.com; s=gamma; h=domainkey-signature:mime-version:received:received:in-reply-to :references:date:message-id:subject:from:to:content-type :content-transfer-encoding; bh=6AKk7mF6G+iFribVgeTCcREs3hocf2jmH7M76jYfS4Q=; b=nKrVK3VImyfbQwGGUNiS9t47a9w2SkPf3Hid0ny+nMsJzNUxd5e535wsOP7u77M8k0 J20C1RCFuEofHbX3l04iZu5wCftpby0AHqE9RIxyAt73Sqg3U6U/kNFME1qNO/dAXxyV rZ4jdIFFs5GS7VsQvA58MRyIc1Xh0MW9T57aI=*Domainkey-signature*: a=rsa-sha1; c=nofws; d=gmail.com; s=gamma; h=mime-version:in-reply-to:references:date:message-id:subject:from:to :content-type:content-transfer-encoding; b=NX6Cw5NeyNzs1NESjy1/fFIykHrDOIFB+LXylfu+pyfuYV3+I8Ed8CJG0kuMx8D9PI 8VYusXpI32pNJ1rvq/Zg8UUv9zfZ7qwN0pmro//mEkdCab8vCMShgh35JbDgzf9tD85i jCtMpbShxBs4djcyNodignpmX4Hp9FAvNPMXY=

2010/7/15 Jitse Niesen <jitse@xxxxxxxxxxxxxxxxx>: > Hi all, > > There seems to be an issue with SelfAdjointEigenSolver when applied to > 3-by-3 matrices, because the unit test eigensolver_selfadjoint_1 fails > occasionally on my computer. > > I think I boiled it down to this example: > > ------------------------------------------------ > > #include <iostream> > #include <Eigen/Core> > #include <Eigen/Eigenvalues> > > using namespace Eigen; > using namespace std; > > int main() > { > Matrix3f A = Matrix3f::Identity(3,3); > A(0,2) = 0.001; > SelfAdjointEigenSolver<Matrix3f> eigOfA(A); > cout << eigOfA.eigenvalues() << endl; > } > > ------------------------------------------------ > > The output of this program is > > 1 > 1 > 1 > > This program computes the eigenvalues of the matrix > > 1 0 0.001 > 0 1 0 > 0.001 0 1 > > They are easily computed by hand; the exact eigenvalues are 0.999, 1, and > 1.001. I think this is rather too big an error. > > The issue seems to be in > > ei_tridiagonalization_inplace_selector<MatrixType,3>::run() > > in Tridiagonalization.h, where we have > > RealScalar v1norm2 = ei_abs2(mat(2,0)); > if (ei_isMuchSmallerThan(v1norm2, RealScalar(1))) > { I dont have time now, but I can already tell you that this code is out of place in a decomposition. It uses the default "dummy" precision which for float is 1e-5f, so it considers v1 to be "zero" if v1norm2 < 1e-5f, which means it considers v1 to be "zero" if its norm is less than sqrt(1e-5) which is 3.1e-3, which is big enough to account for the inaccuracy you are reporting. So this is doubly wrong: - don't use the default dummy precision when doing fuzzy compares in a decomposition. - check for formulas that amount to replacing a precision value by its sqrt, which is much bigger. Benoit > ... > > If the condition in the if statement is true, then mat(2,0) is simply > ignored. This does not seem the right condition, and I'm not even sure why > we need two branches here. Can somebody please help me here? I'd rather not > touch code I don't understand ... > > This is all old code and the same issue is present in Eigen 2, so perhaps we > want to backport the fix before the next version is released? > > Cheers, > Jitse > > >

**References**:**[eigen] Issue in SelfAdjointEigenSolver for 3-by-3 matrices***From:*Jitse Niesen

**Messages sorted by:**[ date | thread ]- Prev by Date:
**[eigen] Issue in SelfAdjointEigenSolver for 3-by-3 matrices** - Next by Date:
**Re: [eigen] Eigen2: Assert failure in QR** - Previous by thread:
**[eigen] Issue in SelfAdjointEigenSolver for 3-by-3 matrices** - Next by thread:
**Re: [eigen] Issue in SelfAdjointEigenSolver for 3-by-3 matrices**

Mail converted by MHonArc 2.6.19+ | http://listengine.tuxfamily.org/ |