Matrix eigenvalues

Roofit and RooStats
1

Hi,

[using ROOT version 5.22/00 on linux and windows]

I want to use ROOT to calculate the eigenvalues and eigenvectors of a symmetric matrix. In my simple example I do not get the expected values. I’m not sure if I did something wrong in my code or if there is a problem in the internal calculation of the eigenvalues.

The output of the attached script gives

2x2 matrix is as follows

 |      0    |      1    |

0 | 1 2
1 | 3 4

Vector (2) is as follows

 |        1  |

0 |5.8541
1 |-0.854102

2x2 matrix is as follows

 |      0    |      1    |

0 | 0.5257 -0.8507
1 | 0.8507 0.5257

But the true eigenvalues are 5.3723 and -0.3723.
The true eigenvectors are (-0.8246, 0.5658) and (-0.4160, -0.9094).

Thanks for any help,
Jochen
matrix2.C (360 Bytes)

2

I just realized that the matrix I used in the example is not good, because it is not symmetric. Using a symmetric matrix I get reasonable results.

However, when I test that the inverted matrix of eigenvectors times the original matrix times the matrix of eigenvectors indeed gives a diagonal matrix with the eigenvalues on the diagonal, then my test fails.

Here is the modified example:

#include <TMatrixDSym.h>
#include <TMatrixDSymEigen.h>
#include <TVectorD.h>

void matrix2(){

  const int N = 3;
  
  double e[N*N] = {
    1, 4, 5,
    4, 2, 6,
    5, 6, 3
  };

  TMatrixDSym m(N, e);
  TMatrixDSymEigen me(m);

  TVectorD eigenval = me.GetEigenValues();
  TMatrixD eigenvec = me.GetEigenVectors();

  m.Print();
  eigenval.Print();
  eigenvec.Print();

  eigenvec.Invert().Print();

  TMatrixD MDiag = eigenvec.Invert() * (m * eigenvec);

  MDiag.Print();

}

The resulting matrix MDIag gives:


     |      0    |      1    |      2    |
--------------------------------------------
   0 |      11.16      0.2471      -1.419 
   1 |      3.988      -1.854       2.157 
   2 |     -2.777      -1.669      -2.606

Thanks again,
Jochen

3

Hi Jochen,

Just noticed your problem in the Stat section.

Your operation eigenvec.Invert() replaces the matrix with its invert !
Replace this expression by TMatrixD(eigenvec,TMatrixD::kInverted).

Eddy