Question on matrix eigendecomposition

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

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

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).

Maybe this hints to the problem?

Maarten

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

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

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

Cheers
Maarten

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

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

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