I’m trying to use CutsD and LikelihoodD in the TMVAClassification Macro. When I run root TMVAClassification.C I get this warning many times.
— Tools : error in matrix diagonalization; printed S and B
In fact I get it 216 times, while I only have 206 events in my signal training sample, and 133256 in background. Now normally I wouldn’t be perturbed by a warning, but all my Cut* optimizations are failing (after training eff_b = eff_s for all points). So I’m trying to take apart my log file one warning at a time. Thanks so much for any thoughts!
I have exaclty the same problem, if I try to run my TMVAClassification.C with MLP and MLPBFGS method I get this error:
— Tools : error in matrix diagonalization; printed S and B
— Tools : error in matrix diagonalization; printed S and B
— Tools : error in matrix diagonalization; printed S and B
— Tools : error in matrix diagonalization; printed S and B
a lot of time…
what is the meaning of it? And how it can be solve? Is the issue connected with the number of entries in the tree?
I would like to ask the same thing. could someone please shed some light on this cryptic message? Is there a way to get more info on it and on what we need to fix?
The code generating this error is in tmva/tmva/src/Tools.cxx:
278 ////////////////////////////////////////////////////////////////////////////////
279 /// square-root of symmetric matrix
280 /// of course the resulting sqrtMat is also symmetric, but it's easier to
281 /// treat it as a general matrix
282
283 TMatrixD* TMVA::Tools::GetSQRootMatrix( TMatrixDSym* symMat )
284 {
285 Int_t n = symMat->GetNrows();
286
287 // compute eigenvectors
288 TMatrixDSymEigen* eigen = new TMatrixDSymEigen( *symMat );
289
290 // D = ST C S
291 TMatrixD* si = new TMatrixD( eigen->GetEigenVectors() );
292 TMatrixD* s = new TMatrixD( *si ); // copy
293 si->Transpose( *si ); // invert (= transpose)
294
295 // diagonal matrices
296 TMatrixD* d = new TMatrixD( n, n);
297 d->Mult( (*si), (*symMat) ); (*d) *= (*s);
298
299 // sanity check: matrix must be diagonal and positive definit
300 Int_t i, j;
301 Double_t epsilon = 1.0e-8;
302 for (i=0; i<n; i++) {
303 for (j=0; j<n; j++) {
304 if ((i != j && TMath::Abs((*d)(i,j))/TMath::Sqrt((*d)(i,i)*(*d)(j,j)) > epsilon) ||
305 (i == j && (*d)(i,i) < 0)) {
306 //d->Print();
307 Log() << kWARNING << "<GetSQRootMatrix> error in matrix diagonalization; printed S and B" << Endl;
308 }
309 }
310 }
311
312 // make exactly diagonal
313 for (i=0; i<n; i++) for (j=0; j<n; j++) if (j != i) (*d)(i,j) = 0;
314
315 // compute the square-root C' of covariance matrix: C = C'*C'
316 for (i=0; i<n; i++) (*d)(i,i) = TMath::Sqrt((*d)(i,i));
317
318 TMatrixD* sqrtMat = new TMatrixD( n, n );
319 sqrtMat->Mult( (*s), (*d) );
320 (*sqrtMat) *= (*si);
321
322 // invert square-root matrices
323 sqrtMat->Invert();
324
325 delete eigen;
326 delete s;
327 delete si;
328 delete d;
329
330 return sqrtMat;
331 }
This piece of code is apparently calculating a transformed positive definite matrix whose eigenvalues are the square root of the original matrix. Unfortunately, this is done in a very round-about way with numerical issues.
The eigen values/vectors are calculated of this symmetric matrix and then matrix multiplications are performed to end up with a diagonal matrix. It is checked that indeed the final matrix is diagonal and its elements positive.
What is actually tested are the numerical limitations of the code that performs the eigen decomposition (TMatrixDSymEigen).
A much more direct way would be to make a Cholesky decomposition with the class TDecompChol. Along the way it also checks that the matrix is positive definite.
@axel , @moneta : The above code should be replaced by these more streamlined lines:
TMatrixD *GetSQRootMatrix2(TMatrixDSym* symMat)
{
TDecompChol chol(*symMat);
if (!chol.Decompose())
{
Log() << kWARNING << "<GetSQRootMatrix> matrix not positive definite; printed S and B" << Endl;
}
TMatrixDSymEigen eigen = TMatrixDSymEigen(*symMat);
TMatrixD s = eigen.GetEigenVectors();
TMatrixDSym d_sqrt(symMat->GetNrows());
TMatrixDDiag d_diag(d_sqrt);
d_diag = eigen.GetEigenValues();
d_sqrt.Sqrt();
// s d s^T
d_sqrt.Similarity(s);
// invert square-root matrices
d_sqrt.Invert();
return new TMatrixD(d_sqrt);
}