How to print each matrix element of the covariance matrix saved in workspace

nsahoo · October 20, 2016, 9:02am

Hi all,

I have saved the co-variance matrix (21x21) information into workspace. In the next step, I want to print the matrix elements of the above co-variance matrix but how to do that ??

Here are the steps followed to print the whole matrix, but I want to print them one by one that means 21 columns and 21 rows …

[nsahoo@lxplus0028 plugins]$ root -l wspace_sigA_bin10.root
root [0]
Attaching file wspace_sigA_bin10.root as _file0…

RooFit v3.60 – Developed by Wouter Verkerke and David Kirkby
Copyright © 2000-2013 NIKHEF, University of California & Stanford University
All rights reserved, please read roofit.sourceforge.net/license.txt

root [1] wspace->Print();

RooWorkspace(wspace) wspace contents

variables

(CosThetaK,CosThetaL,afb,as,effNorm,fl,fs,recK0L0,recK0L2,recK0L4,recK0L6,recK1L0,recK1L2,recK1L4,recK1L6,recK2L0,recK2L2,recK2L4,recK2L6,recK3L0,recK3L2,recK3L4,recK3L6)

p.d.f.s

RooGenericPdf::f_sigA[ actualVars=(CosThetaL,CosThetaK,fl,afb,fs,as,recK0L0,recK1L0,recK2L0,recK3L0,recK0L2,recK1L2,recK2L2,recK3L2,recK0L4,recK1L4,recK2L4,recK3L4,recK0L6,recK1L6,recK2L6,recK3L6) formula=“1.733885e+04*((3.359553e-05exp(-0.5((CosThetaL-(-5.161795e-01))/1.967249e-01)2)+7.791949e-05exp(-0.5((CosThetaL-(5.960371e-04))/3.270870e-01)2)+4.356574e-05exp(-0.5((CosThetaL-(5.628101e-01))/2.302873e-01)2))(4.908191e-05-1.457834e-05CosThetaK+6.357759e-05*CosThetaK2+5.700680e-06*CosThetaK3-9.847393e-05*CosThetaK4-9.202929e-06CosThetaK**5+4.441154e-05CosThetaK6))(1+(recK0L0+recK1L0CosThetaK+recK2L0*(3*CosThetaK2-1)/2+recK3L0*(5CosThetaK**3-3CosThetaK)/2)+(recK0L2+recK1L2CosThetaK+recK2L2(3CosThetaK**2-1)/2+recK3L2(5CosThetaK**3-3CosThetaK)/2)CosThetaL**2+(recK0L4+recK1L4CosThetaK+recK2L4*(3CosThetaK**2-1)/2+recK3L4(5CosThetaK**3-3CosThetaK)/2)CosThetaL**4+(recK0L6+recK1L6CosThetaK+recK2L6*(3CosThetaK**2-1)/2+recK3L6(5CosThetaK**3-3CosThetaK)/2)CosThetaL**6)9/16((2/3fs+4/3as2sqrt(3fs*(1-fs)(0.5+TMath::ATan(fl)/TMath::Pi()))CosThetaK)(1-CosThetaL**2)+(1-fs)(2*(0.5+TMath::ATan(fl)/TMath::Pi())CosThetaK**2(1-CosThetaL2)+0.5*(0.5-TMath::ATan(fl)/TMath::Pi())*(1-CosThetaK2)(1+CosThetaL**2)+4/3(3/2*(1/2-TMath::ATan(fl)/TMath::Pi())TMath::ATan(afb)/TMath::Pi())(1-CosThetaK2)CosThetaL))" ] = 7.03057e-06
RooGenericPdf::f_sigA_MK[ actualVars=(CosThetaK,fl,fs,as) formula="(4.908191e-05-1.457834e-05CosThetaK+6.357759e-05*CosThetaK2+5.700680e-06CosThetaK**3-9.847393e-05CosThetaK4-9.202929e-06*CosThetaK5+4.441154e-05CosThetaK**6)3/4((2/3fs+4/3asCosThetaK)+(1-fs)(2(1/2+TMath::ATan(fl)/TMath::Pi())CosThetaK**2+(1/2-TMath::ATan(fl)/TMath::Pi())(1-CosThetaK**2)))” ] = 1.14681e-05

generic objects

TMatrixTSym::errMtx

root [2] errMtx = (TMatrixDSym*)wspace->obj(“errMtx”);
root [3] errMtx->Print();

21x21 matrix is as follows

 |       0    |       1    |       2    |       3    |       4    |

0 | 6.967e-05 7.145e-06 -6.418e-05 1.015e-06 -0.0003515
1 | 7.145e-06 0.0001311 7.502e-06 -8.256e-05 -4.814e-05
2 | -6.418e-05 7.502e-06 0.0002315 1.46e-05 0.0002326
3 | 1.015e-06 -8.256e-05 1.46e-05 0.0003831 1.946e-06
4 | -0.0003515 -4.814e-05 0.0002326 1.946e-06 0.003202
5 | -6.149e-05 -0.0009499 -2.36e-05 0.00035 0.000619
6 | 0.000204 -2.48e-05 -0.0007389 -5.195e-05 -0.001408
7 | -2.286e-06 0.000265 -5.287e-05 -0.001207 3.186e-05
8 | 1e-30 1e-30 1e-30 1e-30 1e-30
9 | -4.04e-05 4.219e-06 0.000152 1.66e-05 0.0004563
10 | -6.859e-07 -5.195e-05 2.011e-05 0.0002564 -3.9e-05
11 | 7.778e-05 0.001157 1.132e-05 -0.0002759 -0.0008665
12 | 1e-30 1e-30 1e-30 1e-30 1e-30
13 | 1.038e-322 6.917e-323 6.953e-310 1.363e+69 4.52e-153
14 | 1e-30 1e-30 1e-30 1e-30 1e-30
15 | 1e-30 1e-30 1e-30 1e-30 1e-30
16 | 2.356e-310 6.538e-315 0.1185 1.363e+69 4.52e-153
17 | 1e-30 1e-30 1e-30 1e-30 1e-30
18 | 1e-30 1e-30 1e-30 1e-30 1e-30
19 | 1e-30 1e-30 1e-30 1e-30 1e-30
20 | 1.358e-312 6.953e-310 6.953e-310 8.635e-320 6.953e-310

 |       5    |       6    |       7    |       8    |       9    |

0 | -6.149e-05 0.000204 -2.286e-06 1e-30 -4.04e-05
1 | -0.0009499 -2.48e-05 0.000265 1e-30 4.219e-06
2 | -2.36e-05 -0.0007389 -5.287e-05 1e-30 0.000152
3 | 0.00035 -5.195e-05 -0.001207 1e-30 1.66e-05
4 | 0.000619 -0.001408 3.186e-05 1e-30 0.0004563
5 | 0.0131 0.0002252 -0.002215 1e-30 -7.627e-05
6 | 0.0002252 0.004829 0.0003727 1e-30 -0.001719
7 | -0.002215 0.0003727 0.007835 1e-30 -0.0001734
8 | 1e-30 1e-30 1e-30 1e-30 1e-30
9 | -7.627e-05 -0.001719 -0.0001734 1e-30 0.004708
10 | 0.0006632 -0.0001854 -0.002946 1e-30 0.0004721
11 | -0.01982 -0.0002042 0.002125 1e-30 6.776e-05
12 | 1e-30 1e-30 1e-30 1e-30 1e-30
13 | 2.356e-310 1.728e-314 6.953e-310 1e-30 0
14 | 1e-30 1e-30 1e-30 1e-30 1e-30
15 | 1e-30 1e-30 1e-30 1e-30 1e-30
16 | 6.538e-315 4.244e-314 6.953e-310 1e-30 0
17 | 1e-30 1e-30 1e-30 1e-30 1e-30
18 | 1e-30 1e-30 1e-30 1e-30 1e-30
19 | 1e-30 1e-30 1e-30 1e-30 1e-30
20 | 2.356e-310 6.538e-315 6.953e-310 1e-30 3.265e-314

 |      10    |      11    |      12    |      13    |      14    |

 |      15    |      16    |      17    |      18    |      19    |

 |      20    |

0 | 1.358e-312
1 | 6.953e-310
2 | 6.953e-310
3 | 8.635e-320
4 | 6.953e-310
5 | 2.356e-310
6 | 6.538e-315
7 | 6.953e-310
8 | 1e-30
9 | 3.265e-314
10 | 0
11 | 6.933e-315
12 | 1e-30
13 | 9.881e-324
14 | 1e-30
15 | 1e-30
16 | 9.881e-324
17 | 1e-30
18 | 1e-30
19 | 1e-30
20 | 6.953e-310

root [4]

cheers - Niladri

Eddy_Offermann · October 24, 2016, 6:40pm

Hi Niladri,

First of all, please have a look at your diagonal entries; you have a dynamic
range of 10^300 !!!
You should scale your data so that this range becomes on the order of 10^3.

Now coming to your question about accessing individual matrix elements.
Element i,j of wspace can be accessed as wspace(i,j) or wspace[i][j].
If you want to get row i : const TVectorD v = TMatrixDRow_const(wspace,i)

-Eddy

nsahoo · October 25, 2016, 9:45am

Hi Eddy,

Many thanks for your reply. I used another way GetMatrixArray() to access each matrix elements
root.cern.ch/root/html528/TMatr … uble_.html

You told me in previous reply to scale the data so that the range becomes order of 10^3. Could you please tell me how to do ?

cheers - Niladri

Eddy_Offermann · October 25, 2016, 8:38pm

Hi Niladri,

Assuming that you want to invert the covariance as you seem to point out
in another forum entry

[url]Error in <TDecompChol::Decompose()>: matrix not positive definite

the idea is the following. Suppose you have positive definite matrix M,
define a new matrix M_sc, where

M_sc[i][j] = M[i][j]/sqrt(M[i][i]*M[j][j])

The matrix M_sc has now 1 on every diagonal position and off-diagonal <1and > -1,
so numerically a well-behaved matrix

Now invert M_sc => M_sc^-1 and repeat the scaling:

(M_sc^-1)_sc = M_sc^-1[i][j]/sqrt(M[i][i]*M[j][j])

Now the claim is M^-1 = (M_sc^-1)_sc

Proof:

M is positive definite => M = L^T L
define a diagonal matrix S with entry S[i][i] = 1/sqrt(M[i][i])

M_sc = S L^T L S
M_sc^-1 = ( S L^T L S)^-1 = S^-1 L^-1 (L^-1)^T (S^-1)

S M_sc^-1 S = S (S^-1 L^-1 (L^-1)^T (S^-1)) S = L^-1 (L^-1)^T = (L^T L)^-1 = M^-1

Implementation in ROOT:

const TVectorD diag = TMatrixDDiag(m);
const TMatrixDSym m_sc = m;
m.NormByDiag(diag,“D”);
TDecompChol chol(m_sc);
TMatrixDSym m_sc_inv = chol.Invert();
m_inv = m_sc_inv;
m_inv.NormByDiag(diag,“M”);