TMatrixDEigen in python

ROOT
1

Hi,

I’m trying to migrate my C++ fitting script into python and encountered a problem doing eigen-analysis on the covariance matrix.
The following code snippet:

fun = ROOT.TF1("func", "expo", 2.5, 9.5)                                                                                                                          
fitres = hist_data.Fit("func", "IRS")

# deal with uncertainties:                                                                                                                                       
C = fitres.GetCovarianceMatrix()
C.Print()
EigenC = ROOT.TMatrixDEigen(C)

results in:

FCN=1115.07 FROM MIGRAD STATUS=CONVERGED 54 CALLS 55 TOTAL
EDM=1.06952e-09 STRATEGY= 1 ERROR MATRIX ACCURATE
EXT PARAMETER STEP FIRST
NO. NAME VALUE ERROR SIZE DERIVATIVE
1 Constant 1.37419e+01 9.78097e-03 5.02296e-05 5.11911e-04
2 Slope -9.08180e-01 2.58631e-03 1.32823e-05 1.97096e-02

2x2 matrix is as follows

 |      0    |      1    |

0 | 9.567e-05 -2.401e-05
1 | -2.401e-05 6.689e-06

Traceback (most recent call last):
File “test.py”, line 33, in
EigenC = ROOT.TMatrixDEigen(C)
TypeError: none of the 3 overloaded methods succeeded. Full details:
TMatrixDEigen::TMatrixDEigen() =>
takes at most 0 arguments (1 given)
TMatrixDEigen::TMatrixDEigen(const TMatrixT& a) =>
could not convert argument 1
TMatrixDEigen::TMatrixDEigen(const TMatrixDEigen& another) =>
could not convert argument 1

So, I gather the TFitResult object is properly retrieved and I can produce the covariance matrix of the fit. However, TMatrixDEigen constructor does not seem to work with the argument it should accept, to my understanding. The C++ version reads:

…
TMatrixD C = fitres->GetCovarianceMatrix();
TMatrixDEigen* EigenC = new TMatrixDEigen(C);
TMatrixD E = EigenC->GetEigenValues();
TMatrixD V = EigenC->GetEigenVectors();
…

and works for me alright.
I clearly must have missed something here.

Any hints will be much appreciated.

Many thanks and best wishes,
pawel bruckman

2

Hello,

What’s the type of C in the snippet? Can you print(type(C))?

Also, what ROOT version are you using?

3

Hi,

Print(type(C) gives:

<class ‘ROOT.TMatrixTSym’>

and I’m currently using ROOT 6.16/00

Thanks!
pawel

4

Hi, I’m sorry. Perhaps I should have replied to you directly. Anyhow, I replied to your questions in a general post. I’m new to the forum :wink:

5

Hi,

First I think it would be worth trying it with a newer version of ROOT (>= 6.22).

But looking at the code, it seems that the type of the result of GetCovarianceMatrix is TMatrixTSym<double>, aka TMatrixDSym. However, the constructor of TMatrixDEigen accepts a TMatrixT.

6

Well, this is true. But I thought TMatrixDSym derived from TMatrixT. Otherwise, it wouldn’t work in C++ Root, would it? Is there a way to force an appropriate type conversion? I shall check it with a higher Root version, as you recommend.
p.

7

According to:

https://root.cern.ch/doc/v622/classTMatrixTBase.html

Both TMatrixTSym and TMatrixT are subclasses of TMatrixTBase, but there is no direct inheritance between them. Am I right that TMatrixDSym is TMatrixTSym<double>?

8

This I do not know for sure, but very likely it is TMatrixTSym. My trouble is that in C++ Root TMatrixDEigen accepts both TMatrixDSym and TMatrixD in the constructor. I just checked it explicitly.

9

I also checked that in PyRoot ROOT.TMatrixDEigen behaves the same no matter if you pass the TMatrixTSym or TMatrixT. The problem seems in interpreting the reference “&”. This is what I get when TMatrixT is passed to the constructor:
<class ‘ROOT.TMatrixT’>
Traceback (most recent call last):
File “test.py”, line 40, in
EigenC = ROOT.TMatrixDEigen(C)
TypeError: none of the 3 overloaded methods succeeded. Full details:
TMatrixDEigen::TMatrixDEigen() =>
takes at most 0 arguments (1 given)
TMatrixDEigen::TMatrixDEigen(const TMatrixT& a) =>
could not convert argument 1
TMatrixDEigen::TMatrixDEigen(const TMatrixDEigen& another) =>
could not convert argument 1

Maybe there is a simple workaround I’m not aware of?

10

I temporarily solved my problem by simply converting the TMatrixDSym object into a numpy array. Then I used scipy.linalg.eigh to solve my eigenproblem.
Still, this is not a satisfactory solution to the original issue, of course :frowning:
p.

11

The missing link

TMatrixT(const TMatrixTSym   <Element> &another);
12

Would adding the following line be helpful:

TMatrixD EigenC = ROOT.TMatrixDEigen(C)
or
EigenC = ROOT.TMatrixDEigen(TMatrixD(C))
13

Hi Eddy,
Bingo! The explicit type conversion does work indeed (just adding the ROOT namesapace):
EigenC = ROOT.TMatrixDEigen(ROOT.TMatrixD(C))
The first option throws a syntax error, though.
Thanks a lot for this very helpful suggestion :slight_smile:
Best,
pawel

14

first option should have been:

EigenC = TMatrixD(ROOT.TMatrixDEigen(C))
15

No, this version still has the same problem with ROOT.TMatrixDEigen(C) which is not accepted.
So presumably, we are left with the second option. Still fair enough! :slight_smile:
Thanks again,
p.

16

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