New user self-registration is currently disabled. Please email eigen-core-team @ lists.tuxfamily.org if you need an account.
Bug 1014 - Eigenvalues 3x3 matrix
Eigenvalues 3x3 matrix
Status: CONFIRMED
Product: Eigen
Classification: Unclassified
Component: Eigenvalues
3.2
All All
: Normal Feature Request
Assigned To: Nobody
:
Depends on:
Blocks: 3.4
  Show dependency treegraph
 
Reported: 2015-05-07 12:45 UTC by thomas.rueberg
Modified: 2015-09-03 10:09 UTC (History)
4 users (show)



Attachments
Test code (780 bytes, text/x-c++src)
2015-05-07 12:45 UTC, thomas.rueberg
no flags Details
Patch (447 bytes, patch)
2015-05-07 12:46 UTC, thomas.rueberg
no flags Details | Diff
Make direct eigenvalue computation of 3x3 matrices more stable by shifting the eigenvalues (922 bytes, patch)
2015-05-11 15:28 UTC, Christoph Hertzberg
no flags Details | Diff
Updated root computation (5.02 KB, patch)
2015-05-12 15:10 UTC, Christoph Hertzberg
no flags Details | Diff
Compute 3x3 eigenvectors using nullspace from specialized LU (8.28 KB, patch)
2015-05-12 16:56 UTC, Gael Guennebaud
no flags Details | Diff
Stablely working version (9.15 KB, patch)
2015-05-16 10:42 UTC, Christoph Hertzberg
no flags Details | Diff

Description thomas.rueberg 2015-05-07 12:45:35 UTC
Created attachment 575 [details]
Test code

There is floating point comparison issue with the routine in SelfAdjointEigenSolver.h (Version 3.2.4) for fixed-size 3x3 matrices.
A tolerance is created as follows:
Scalar safeNorm2 = Eigen::NumTraits<Scalar>::epsilon();

In an older version (3.2.2), if followed the line
safeNorm2 *= safeNorm2;
which has been omitted. For matrices which are close to a multiple of the identity matrix, this change becomes crucial in some comparisons.
Comment 1 thomas.rueberg 2015-05-07 12:46:22 UTC
Created attachment 576 [details]
Patch
Comment 2 Gael Guennebaud 2015-05-11 14:26:09 UTC
this seems to conflict with bug 854
Comment 3 Christoph Hertzberg 2015-05-11 15:28:09 UTC
Created attachment 577 [details]
Make direct eigenvalue computation of 3x3 matrices more stable by shifting the eigenvalues

We can make the computation significantly more stable by shifting the matrix by the mean of the eigenvalues (i.e. mat.trace()/3) -- note that (in exact math) this shifts the eigenvalues but does not influence the eigenvectors. We could then also safe the calculation of c2 in computeRoots.

I did not check yet if we could also save some cases of the eigenvector computation, i.e., is it even possible that the 3rd cross product is small? If the original matrix had NaNs or Infs, computeRoot should have failed already.
Comment 4 Gael Guennebaud 2015-05-11 15:43:57 UTC
This shift definitely helps a lot in computing more accurate eigenvalues. I'd only suggest doing ".diagonal().array() -= shift" instead of using Identity.
Comment 5 Gael Guennebaud 2015-05-11 16:06:58 UTC
(In reply to Christoph Hertzberg from comment #3)
> I did not check yet if we could also save some cases of the eigenvector
> computation, i.e., is it even possible that the 3rd cross product is small?
> If the original matrix had NaNs or Infs, computeRoot should have failed
> already.

I've checked for several matrices, and any one or twos of the 3 cases can be arbitrarily small.
Comment 6 Gael Guennebaud 2015-05-12 12:14:23 UTC
After more tests, the current cross-product based approach is not very good. It is really reliable only if we compute all three possibilities, and then pick the best, but such a approach kills the performance (recall that high performance is the main motivation for this "direct" path).

So here is an alternative to the cross product approach based on the direct extraction of the nullspace:

{
  // Compute the kernel of 'tmp' using a LU factorization of tmp with symmetric pivoting.
  // Note that we only need to compute the factor U, which is of the form:
  //   x x x
  //   0 x x
  //   0 0 0
  Index i0, i1, i2;
  Scalar p0, p1, p2, f;
  // Find the first pivot of index i0:
  tmp.diagonal().cwiseAbs().maxCoeff(&i0);
  p0 = Scalar(1)/tmp.coeff(i0,i0);
  // The other two rows/columns are:
  i1 = (i0+1)%3;
  i2 = (i0+2)%3;
  // Update the first row of U:
  tmp.coeffRef(i0,i1) *= p0;
  tmp.coeffRef(i0,i2) *= p0;
  // Thanks to symmetry, the remaining 2x2 block is of the form:
  //     p1 f
  //     f  p2
  // with:
  p1 = tmp.coeffRef(i1,i1) - tmp.coeff(i1,i0) * tmp.coeff(i0,i1);
  p2 = tmp.coeffRef(i2,i2) - tmp.coeff(i2,i0) * tmp.coeff(i0,i2);
  f  = tmp.coeffRef(i1,i2) - tmp.coeff(i1,i0) * tmp.coeff(i0,i2);
  // Second pivoting:
  if(std::abs(p1) < std::abs(p2))
  {
    std::swap(i1,i2);
    p1 = p2;
  }
  // Back-subtitutions to form the 1D kernel of 'tmp':
  eivecs.col(k).coeffRef(i2) = 1;
  eivecs.col(k).coeffRef(i1) = - f / p1;
  eivecs.col(k).coeffRef(i0) = - (tmp.coeff(i0,i2) + tmp.coeff(i0,i1) * eivecs.col(k).coeff(i1) );
  eivecs.col(k).normalize();
}


It is clearly more costly than doing a single cross product, but it should be faster than doing two.

With both patches, I still haven't find a failing case after running the unit tests for hours...
Comment 7 Christoph Hertzberg 2015-05-12 15:03:32 UTC
(In reply to Gael Guennebaud from comment #6)
> So here is an alternative to the cross product approach based on the direct
> extraction of the nullspace: [...]

Sounds good to me. Better a reliable, reasonably fast implementation than a slightly faster implementation which occasionally fails. We will still provide a determinable maximal running time.
A possible flaw with the proposed implementation is that diagonal().maxCoeff() might be zero (though probably that's caught outside the fragment you posted).

> With both patches, I still haven't find a failing case after running the
> unit tests for hours...

We probably also should extend our unit tests by some 'trivial' border cases (e.g., diagonal matrices with very close eigenvalues, as in the original report -- probably with small off-diagonal noise).


I'll provide an updated patch for root computation in the next comment.
Comment 8 Christoph Hertzberg 2015-05-12 15:10:01 UTC
Created attachment 578 [details]
Updated root computation

For the computation of roots exploit that m.trace()==0 after shifting the matrix.

Actually, the computation should work for complex (self-adjoint) matrices as well, after replacing every other m(i,j) by conj(m(i,j)) [which of course can be done implicitly], and knowing that imag(m(i,i))==0.
Comment 9 Gael Guennebaud 2015-05-12 16:56:37 UTC
Created attachment 579 [details]
Compute 3x3 eigenvectors using nullspace from specialized LU

Here is a more complete patch for the new eigenvector extraction method.

In the two calls to extract_kernel we are guaranteed that there is at least one meaningful non zero on the diagonal (because we are shifting wrt the minimum and maximum eigenvalues and the case were all eigenvalues are the same is already handled).

However, I'm 100% sure that a symmetric pivoting is guaranteed to succeed (compared to a full pivoting), that is, that we are guaranteed to find a second meaningful pivot along the diagonal...
Comment 10 Gael Guennebaud 2015-05-12 17:00:53 UTC
Regarding the unit tests:

https://bitbucket.org/eigen/eigen/commits/01a9e56db770/
Summary:     Make test matrices for eigensolver/selfadjoint even more tricky

https://bitbucket.org/eigen/eigen/commits/0ba2cd3d31c8/
Summary:     Add regression test for bugs 854 and 1014, and check that the eigenvector matrix is unitary.
Comment 11 Gael Guennebaud 2015-05-12 17:01:44 UTC
I haven't tested it yet, but your updated computeRoots looks good to me.
Comment 12 Christoph Hertzberg 2015-05-12 21:55:06 UTC
I found a problem (with both patches applied) when running:
$ ./test/eigensolver_selfadjoint_13 s1000

In that test a matrix is generated with a difference in the diagonal elements of just 4 ULPs. This makes the shift inexact and results in a scaledMatrix with a trace of about 0.66, thus returning bogus eigenvalues (with the patched computeRoots). This barely matters for the final eigenvalues but causes problems computing the eigenvectors.

A pragmatic solution would be to short-circuit the eigenvector computation if scale < abs(shift)*epsilon*8 (maybe 16) and return an identity matrix for the eigenvectors with the shift value as eigenvalues.
Or we just keep the current computeRoots function (it's not too many additional computations ...).
Comment 13 Gael Guennebaud 2015-05-13 10:09:08 UTC
Right, I also observed this kind of problem. I've also locally added a comparison between the "direct" eigenvalues and the ones coming from the iterative algorithm and assuming that trace==0 severely affects the accuracy of the eigenvalues, not only the one of the eigenvector computation.

I've also run some benchmark, and the gain of assuming that trace==0 is negligible: computeRoots is dominated by the trigonometric function calls.

I've also observed that my new eigenvector computation is really slower, so I came up with the following version where the idea is to pick the column having the biggest diagonal element, and then pick the best of the two cross-product possibilities: (seems to be 15% faster even though this imply more muls and adds):

Index i0;
tmp.diagonal().cwiseAbs().maxCoeff(&i0);
Scalar n0, n1;
VectorType c0, c1;
eivecs.col(l) = tmp.col(i0);
n0 = (c0 = tmp.col(i0).cross(tmp.col((i0+1)%3))).squaredNorm();
n1 = (c1 = tmp.col(i0).cross(tmp.col((i0+2)%3))).squaredNorm();
if(n0>n1) eivecs.col(k) = c0/std::sqrt(n0);
else      eivecs.col(k) = c1/std::sqrt(n1);

We can also save some redundant computation between the two cross products (at the cost of readability, so I'll check whether the compiler already managed to get rid of them).
Comment 14 Christoph Hertzberg 2015-05-13 16:01:27 UTC
I guess to determine i0 we rather need to select the column with the largest element (or the largest norm), e.g. consider this matrix (with e=epsilon):
 2e  0  0
  0 -e  1
  0  1 -e

But then I'm wondering if simply computing all 3 cross-products may not be faster? (exploiting the symmetry)
Comment 15 Christoph Hertzberg 2015-05-16 10:42:36 UTC
Created attachment 580 [details]
Stablely working version

It seems my last concerns are not justified. A near-zero diagonal can only occur for the "middle" eigenvector, which is calculated differently anyways.

I attached a patch combining Gael's last proposal with only minor improvements in computeRoots.
extractKernel now always returns true, so the if-branches should actually be removed. Also, we can save some computations and copies by exploiting the symmetry better.

Furthermore, all computations could be ported to complex self-adjoint matrices.
Comment 16 Christoph Hertzberg 2015-05-17 20:01:21 UTC
I pushed the last patch after some clean-up here:
https://bitbucket.org/eigen/eigen/commits/ee4084b1f42c94b3ed820

Still must be backported to 3.2 and probably optimized before that.
Comment 17 Gael Guennebaud 2015-05-18 20:15:06 UTC
(In reply to Christoph Hertzberg from comment #15)
> It seems my last concerns are not justified. A near-zero diagonal can only
> occur for the "middle" eigenvector, which is calculated differently anyways.

Yes, that's exactly why I reversed the computation of the closest two eigenvectors.

Thanks for pushing the changes.
Comment 18 Gael Guennebaud 2015-06-15 20:17:43 UTC
I backported it to 3.2, the optimization opportunities are very marginal, so good enough for 3.2.5, and not a priority for 3.3 either.

https://bitbucket.org/eigen/eigen/commits/5475a5ded9b0/
Comment 19 Gael Guennebaud 2015-09-03 10:09:37 UTC
To summarize, it remains to add support for complexes, and check for some speed improvement opportunities.

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