Can ROOT compute a covariance matrix?

Dear ROOTers,

I have 1000 events, and for each event the values of 5 variables in an ntuples.

I want to get the covariance matrix by using the distributions of these variables. Is there already a class that does that in ROOT ?

Thanks in advance,

Xavier

You can look at the class TPrincipal and example at $ROOTSYS/tutorials/math/principal.C
root.cern.ch/root/html/TPrincipal.html
root.cern.ch/root/html/tutorials … pal.C.html

Rene

There is also the class TRobustEstimator for a robust estimation

See root.cern.ch/root/htmldoc/TRobustEstimator.html

Lorenzo

Thanks for your answers,

I did a simple test with TPrincipal:

First filling an ntuple with 3 variables 100% correlated

void fill()
{
  TFile *f = new TFile("test.root", "RECREATE");

  TNtuple* ntuple = new TNtuple("ntuple","ntuple","x:y:z");
  for( int i=0; i<100000; i++){
    float x = i;
    float y = i;
    float z = i;
    ntuple->Fill(x,y,z); 
  }
    ntuple->Write();
  f->Close();
   return;
}

Then computing the covariance matrix and printing it:

 void read()
{
  TFile* f = new TFile("test.root");
  TNtuple* ntuple = (TNtuple*) f->Get("ntuple");
  Float_t x;
  Float_t y;
  Float_t z;
  ntuple->SetBranchAddress("x",&x);
  ntuple->SetBranchAddress("y",&y);
  ntuple->SetBranchAddress("z",&z);

  TPrincipal p( 3, "ND");
  Double_t data[3];

  for ( int ient = 0; ient < ntuple->GetEntries(); ient++ ) {
    ntuple->GetEntry(ient);
    data[0] = x;
    data[1] = y;
    data[2] = z;
    p.AddRow( data );
  }

  TMatrixD* m = p.GetCovarianceMatrix();
  m->Print();
  return;
}

The result is:

3x3 matrix is as follows

 |      0    |      1    |      2    |

0 | 8.333e+08 0 0
1 | 8.333e+08 8.333e+08 0
2 | 8.333e+08 8.333e+08 8.333e+08

I have then 2 questions:
1/ the matrix elements are not 1, obviously a normalization concern, but should’nt the option “N” normalize the covariance matrix ?

2/ the matrix is not symmetric, is that made to save computing ressources ? is there an automated way to symmetrize it ?

Xavier

The doc from TPrincipal claims that the matrix released by GetCovarianceMatrix() is symmetric, but checking that with

  TPrincipal* p = new TPrincipal( n, "ND");
  Double_t data[4];
  p->AddRow( data );
  TMatrixD* p->GetCovarianceMatrix();
  p->Print();
  cout << p->IsSymmetric()  << endl;

returns 0.

Hi Xavier,

For reasons unknown, the author stored the only the lower left
part of the matrix in a TMatrixD.

However, when calculating the eigen vectors, symmetric algorithms
are used and the matrix is converted to a TMatrixDSym:

TMatrixD* m = p.GetCovarianceMatrix();
TMatrixDSym sym; sym.Use(m->GetNrows(),m->GetMatrixArray());

(Notice the “Use” which makes sure that we do not copy data back and
forth)

You can do this too in your code and of course this should be updated
in the TPrincipal class.

Concerning your normalization question:

Looking at the code the “N” option makes sure that the covariance
matrix is normalized before calculating the eigen vectors and values.
HOWEVER, the matrix stored in the class is not touched in other
words GetCovarianceMatrix will always give you the same matrix.

You can create yourself this normalized matrix by doing the
following operation on your matrix m:

const TVectorD diag = TMatrixDDiag_const(*m);
m->NormbyDiag(diag);

Eddy

Hi !

I have found this thread about covariance matrix and it helps. But I have some more questions.
When I’m doing

 const TMatrixD *ms=ps.GetCovarianceMatrix();
 TMatrixDSym sym; 
 sym.Use(ms->GetNrows(),ms->GetMatrixArray());
 ms->Print(); 
 sym.Print();

I get :

4x4 matrix is as follows

     |      0    |      1    |      2    |      3    |
---------------------------------------------------------
   0 |       1.59           0           0           0 
   1 |   -0.03833       1.209           0           0 
   2 |   -0.02875      0.2255      0.3314           0 
   3 |   0.007018      0.8823      0.2156       3.461 


4x4 matrix is as follows

     |      0    |      1    |      2    |      3    |
---------------------------------------------------------
   0 |       1.59           0           0           0 
   1 |   -0.03833       1.209           0           0 
   2 |   -0.02875      0.2255      0.3314           0 
   3 |   0.007018      0.8823      0.2156       3.461 

The supposed symetric matrix is not.
Am I doing things wrong ?

Hi Maddalena,

These operations will not force a matrix to be symmetric. For that do something like:

const TMatrixD* m = p.GetCovarianceMatrix();
TMatrixD mt = *m; mt.T();
TMatrixDSym sym = *m+mt;

Eddy

Hi

Is there something in root that allows to do the same thing as here with a multivariate_normal_cdf :
roguewave.com/portals/0/prod … malcdf.htm

Thank you !

Hi,

I know this post is from 2012 but:

  1. M + M^T does not symmetrize the matrix because you are adding the diagonal twice.
  2. Maddalena is right, the covariance matrix IS symmetric, what you have here is what we call a bug. The function should return a TMatrixDSym. It’s been years and no one has done anything about this.

I think the safest way would be to just loop over the matrix to make a symmetric version.

Cheers.

rooter_03,

  1. the symmetric matrix has the wrong diagonal, just make this one line modification :
const TMatrixD* m = p.GetCovarianceMatrix();
TMatrixD mt = *m; mt.T();
TMatrixDDiag d(mt); d = 0;
TMatrixDSym sym = *m+mt;
  1. either GetCovarianceMatrix() should be documented properly or (better) the
    code should be updated; any volunteers ?

-Eddy

Hi Eddy,

I will do it. As I understand it I should download the ROOT sources, find the place in the code where this is done, modify it, recompile it and test it. When I be sure that it is fixed I should tell you what the changes should be, so that you update the source code. Is this correct?

Cheers.

Hi,

My contributions used to flow through Rene Brun. The proper way to contribute
now is described here:

root.cern.ch/collaborate-with-us

-Eddy

Hi,
Is there some news about this issue? I am trying to use TPrincipal::GetCovarianceMatrix right now but I have the same problem. The matrix is not symetric and show me just its lower left. I need this to be symetric in order to do some math. Could you guys please suggest another method?

Thank you in advance

Dear Gabriel,

Any reason not to use my recipe proposed in Jan 2017 :slight_smile:

const TMatrixD* m = p.GetCovarianceMatrix();
TMatrixD mt = *m; mt.T();
TMatrixDDiag d(mt); d = 0;
TMatrixDSym sym = *m+mt;  

As stated in the TPrincipal description, it is an implementation of the LINTRA package which was developed at a time when memory was a KB and not a GB issue …

Changing GetCovariance could be a backward compatibility issue, unless it is done by adding a new option variable.
Probably better to update the documentation.

Any volunteers ?

rooter_03 …:

  • saw that his analysis is not Nobel-prize suspicious
  • dropped this “physics - thing” and created a startup with a planned IPO in 2020
  • rewrote his code in Python and is still waiting for an answer …

-Eddy