Question on matrix eigendecomposition

ROOT
1

Hello experts,

I try to eigen-decompose a (seemingly simple) symmetric 3x3 matrix M, defined like this:

        1        0.01        0.01 
     0.01           1        0.01 
     0.01        0.01           1

But TMatrixDEigen seems to fail. Fetching the eigenvalues D and eigenvectors Q, I calculate Q D Q^T and get:

    1.055     -0.1825      0.1475 
  -0.1825       1.179     0.02375 
   0.1475     0.02375      0.8487

which is not equal to M, so the decomposition failed.
I tried this on several platforms, with identical results.
I attach a macro to reproduce the issue.

Strangely enough the decomposition most often works normally (for example, changing the off-diagonal elements to 0.1 or 0.0001), but with some specific cases like the one above it seems to break. I was not able to determine more precise criteria.

Can you help?
Thanks,
Maarten

test.C (672 Bytes)
ROOT Version: 6.22
Platform: EL7
Compiler: Not Provided

2

Hi,

What you are doing to verify the decomposition is not correct.

M = Q D Q^-1 and not Q^T

see for example Eigendecomposition of a matrix - Wikipedia

If I change your test program using the inverse of the eigenvectors matrix Q, it works fine

Cheers

Lorenzo

3

Hi!

M is a square symmetric matrix, in which case the relation should hold.
I should precise that using TMatrixDSymEigen I get the same problem (if this is relevant).

(Eigendecomposition of a matrix - Wikipedia)

Maybe this hints to the problem?

Maarten

4

I wanted to add one more observation. In the example above, Q (the matrix whose columns are the eigenvectors of M) comes out as

  0.5774        -0.8165       0.2357 
  0.5774         0.4082       -0.825  
  0.5774         0.4082       0.5893

In python, when I do this :

import numpy as np
B = np.array([[1,.01,0.01], [.01,1,.01], [0.01,.01,1]])
np.linalg.eig(B)

The matrix of eigenvectors comes out as

   array([[-0.81649658,  0.57735027,  0.40792387],
          [ 0.40824829,  0.57735027, -0.8164965 ],
          [ 0.40824829,  0.57735027,  0.40857263]]))

The order of the columns is arbitrary and depend on how the eigenvalues are ordered.
The 1st and 2nd columns are switched but match; the 3rd columns however differ, and not by just an overall scale. So the two outcomes are not equivalent.

Cheers,
Maarten

5

Hi,
It is true that M i symmetric, however I am not sure this is always true, that you get orthonormal eigenvectors, especially when some of the eigenvalues are the same as in this case.
For example computing eigenvalues and eigenvector in Mathematica, see

I don’t have an orthogonal matrix for Q (called S ) in this case

Lorenzo

6

OK, thanks for the clarification. It might then be the specific purpose of TMatrixDSym / TMatrixDSymEigen to return orthonormal eigenvectors?

I will change the title of this post :slight_smile:

Cheers
Maarten

7

Hi,
This is correct, when you use TMatrixDSymEigen, the return eigenvectors are orthonormal and your test programs returns correctly M = Q D Q^-T.

I attached here your modified code for reference
test_TMatrixDSymEigen.C (790 Bytes)

Best regards

Lorenzo

8

Hi Maarten,

Of course TMatrixD / TMatrixDEigen should have returned the correct orthonormal eigenvectors. However, it seems that the (Wilkinson/EISPACK) algorithm fails and the third eigenvector for the degenerate eigenvalue 0.99 is wrong. I have no idea why but in general for accuracy and speed it is wise to use algorithms that are aware of certain properties like symmetry.
-Eddy

1 Like
9

This topic was automatically closed 14 days after the last reply. New replies are no longer allowed.