[Math][Minuit2] New Minuit2 strategy for improved Hessian calculation

New Minuit2 strategy for improved Hessian calculation and return without
making positive definite.

This proposed new Strategy in minuit2 is the same migrad behaviour as
Strategy=2 but with the following changes to the Hesse calculation:

  * The step and g2 tolerances have been zeroed so that the 7 cycles are
    forced in the diagonal-term calculation. This was found to be
    necessary in cases where asimov datasets were used for the
    minimization and there were very few iterations for the approximate
    covariance to be determined from.

  * Central finite difference is used for mixed partial derivatives.
    This requires 3 extra function evaluations per derivative, but is
    necessary in the case of minima where there is high curvature (in
    the case of high stats) and the forward finite difference (default)
    behaviour leads incorrectly to a non-positive-definite covariance
    matrix

  * Return the uncorrected covariance matrix, even if it is not positive
    definite. This useful for checking just how far from
    positive-definiteness the matrix is by being able to examine the
    eigenvalues.

Additionally, a lower bound on the precision allowed for the spread of
eigenvalues of the "hessian" correlation matrix (computing a correlation
matrix with the hessian as if it was a covariance matrix) was reduced
from 1e-6 to 1e-12 (see MnHesse.cxx) ... it is not clear why 1e-6 was
the lower bound previously, but current machine precision can beat that
(I get locally 1e-8). I left a comment about whether this lower bound
should be made configurable or not...

This new strategy was tested with a model with high statistics (almost
50 million events) where the migrad minimization was successful but the
hessian was being forced positive definite. With this new Strategy 3 the
hessian is accurate and positive definite in all cases tested.

Author: will-cern

math/minuit2/inc/Minuit2/MinimumError.h

@@ -76,11 +76,12 @@ class MinimumError {
    double Dcovar() const { return fPtr->fDCovar; }
    Status GetStatus() const { return fPtr->fStatus; }
 
-   bool IsValid() const { return IsAvailable() && (IsPosDef() || IsMadePosDef()); }
+   bool IsValid() const { return IsAvailable() && (IsPosDef() || IsMadePosDef() || IsNotPosDef()); }
    bool IsAccurate() const { return IsPosDef() && Dcovar() < 0.1; }
 
    bool IsPosDef() const { return GetStatus() == MnPosDef; }
    bool IsMadePosDef() const { return GetStatus() == MnMadePosDef; }
+   bool IsNotPosDef() const { return GetStatus() == MnNotPosDef; }
    bool HesseFailed() const { return GetStatus() == MnHesseFailed; }
    bool InvertFailed() const { return GetStatus() == MnInvertFailed; }
    bool HasReachedCallLimit() const { return GetStatus() == MnReachedCallLimit; }

math/minuit2/inc/Minuit2/MnStrategy.h

@@ -16,8 +16,8 @@ namespace Minuit2 {
 
 //_________________________________________________________________________
 /**
-    API class for defining three levels of strategies: low (0), medium (1),
-    high (>=2);
+    API class for defining four levels of strategies: low (0), medium (1),
+    high (2), very high (>=3);
     acts on: Migrad (behavioural),
              Minos (lowers strategy by 1 for Minos-own minimization),
              Hesse (iterations),
@@ -30,7 +30,7 @@ class MnStrategy {
    // default strategy
    MnStrategy();
 
-   // user defined strategy (0, 1, >=2)
+   // user defined strategy (0, 1, 2, >=3)
    explicit MnStrategy(unsigned int);
 
    ~MnStrategy() {}
@@ -45,16 +45,20 @@ class MnStrategy {
    double HessianStepTolerance() const { return fHessTlrStp; }
    double HessianG2Tolerance() const { return fHessTlrG2; }
    unsigned int HessianGradientNCycles() const { return fHessGradNCyc; }
+   unsigned int HessianCentralFDMixedDerivatives() const { return fHessCFDG2; }
+   unsigned int HessianForcePosDef() const { return fHessForcePosDef; }
 
    int StorageLevel() const { return fStoreLevel; }
 
    bool IsLow() const { return fStrategy == 0; }
    bool IsMedium() const { return fStrategy == 1; }
-   bool IsHigh() const { return fStrategy >= 2; }
+   bool IsHigh() const { return fStrategy == 2; }
+   bool IsVeryHigh() const { return fStrategy >= 3; }
 
    void SetLowStrategy();
    void SetMediumStrategy();
    void SetHighStrategy();
+   void SetVeryHighStrategy();
 
    void SetGradientNCycles(unsigned int n) { fGradNCyc = n; }
    void SetGradientStepTolerance(double stp) { fGradTlrStp = stp; }
@@ -65,6 +69,14 @@ class MnStrategy {
    void SetHessianG2Tolerance(double toler) { fHessTlrG2 = toler; }
    void SetHessianGradientNCycles(unsigned int n) { fHessGradNCyc = n; }
 
+   // 1 = calculate central finite difference mixed derivatives (involves 3 extra evaluations per derivative)
+   // 0 = use forward finite difference (default)
+   void SetHessianCentralFDMixedDerivatives(unsigned int flag) { fHessCFDG2 = flag; }
+
+   // 1 = returned matrix from Hesse should be forced positive definite (default)
+   // 0 = do not force matrix positive definite
+   void SetHessianForcePosDef(unsigned int flag) { fHessForcePosDef = flag; }
+
    // set storage level of iteration quantities
    // 0 = store only last iterations 1 = full storage (default)
    void SetStorageLevel(unsigned int level) { fStoreLevel = level; }
@@ -79,6 +91,8 @@ class MnStrategy {
    double fHessTlrStp;
    double fHessTlrG2;
    unsigned int fHessGradNCyc;
+   int fHessCFDG2;
+   int fHessForcePosDef;
    int fStoreLevel;
 };
 

math/minuit2/src/Minuit2Minimizer.cxx

@@ -467,7 +467,10 @@ bool Minuit2Minimizer::Minimize()
 
    // set strategy and add extra options if needed
    ROOT::Minuit2::MnStrategy strategy(strategyLevel);
-   ROOT::Math::IOptions *minuit2Opt = ROOT::Math::MinimizerOptions::FindDefault("Minuit2");
+   const ROOT::Math::IOptions *minuit2Opt = fOptions.ExtraOptions();
+   if (!minuit2Opt) {
+      minuit2Opt = ROOT::Math::MinimizerOptions::FindDefault("Minuit2");
+   }
    if (minuit2Opt) {
       // set extra  options
       int nGradCycles = strategy.GradientNCycles();
@@ -495,7 +498,7 @@ bool Minuit2Minimizer::Minimize()
       strategy.SetGradientTolerance(gradTol);
       strategy.SetGradientStepTolerance(gradStepTol);
       strategy.SetHessianStepTolerance(hessStepTol);
-      strategy.SetHessianG2Tolerance(hessStepTol);
+      strategy.SetHessianG2Tolerance(hessG2Tol);
 
       int storageLevel = 1;
       bool ret = minuit2Opt->GetValue("StorageLevel", storageLevel);
@@ -551,7 +554,7 @@ bool Minuit2Minimizer::Minimize()
       fMinimum = new ROOT::Minuit2::FunctionMinimum(min);
    }
 
-   // check if Hesse needs to be run. We do it when is requested (IsValidError() == true)
+   // check if Hesse needs to be run. We do it when is requested (IsValidError() == true , set by SetParabError(true) in fitConfig)
    // (IsValidError() means the flag to get correct error from the Minimizer is set (Minimizer::SetValidError())
    // AND when we have a valid minimum,
    // AND  when the the current covariance matrix is estimated using the iterative approximation (Dcovar != 0 , i.e. Hesse has not computed  before)
@@ -1202,7 +1205,43 @@ bool Minuit2Minimizer::Hesse()
       return false;
    }
 
-   const int strategy = Strategy();
+
+   // set strategy and add extra options if needed
+   ROOT::Minuit2::MnStrategy strategy(Strategy());
+   const ROOT::Math::IOptions *minuit2Opt = fOptions.ExtraOptions();
+   if (!minuit2Opt) {
+      minuit2Opt = ROOT::Math::MinimizerOptions::FindDefault("Minuit2");
+   }
+   if (minuit2Opt) {
+      // set extra  options
+      int nGradCycles = strategy.GradientNCycles();
+      int nHessCycles = strategy.HessianNCycles();
+      int nHessGradCycles = strategy.HessianGradientNCycles();
+
+      double gradTol = strategy.GradientTolerance();
+      double gradStepTol = strategy.GradientStepTolerance();
+      double hessStepTol = strategy.HessianStepTolerance();
+      double hessG2Tol = strategy.HessianG2Tolerance();
+
+      minuit2Opt->GetValue("GradientNCycles", nGradCycles);
+      minuit2Opt->GetValue("HessianNCycles", nHessCycles);
+      minuit2Opt->GetValue("HessianGradientNCycles", nHessGradCycles);
+
+      minuit2Opt->GetValue("GradientTolerance", gradTol);
+      minuit2Opt->GetValue("GradientStepTolerance", gradStepTol);
+      minuit2Opt->GetValue("HessianStepTolerance", hessStepTol);
+      minuit2Opt->GetValue("HessianG2Tolerance", hessG2Tol);
+
+      strategy.SetGradientNCycles(nGradCycles);
+      strategy.SetHessianNCycles(nHessCycles);
+      strategy.SetHessianGradientNCycles(nHessGradCycles);
+
+      strategy.SetGradientTolerance(gradTol);
+      strategy.SetGradientStepTolerance(gradStepTol);
+      strategy.SetHessianStepTolerance(hessStepTol);
+      strategy.SetHessianG2Tolerance(hessG2Tol);
+   }
+
    const int maxfcn = MaxFunctionCalls();
    print.Info("Using max-calls", maxfcn);
 
@@ -1254,6 +1293,8 @@ bool Minuit2Minimizer::Hesse()
       covStatusType = "full but made positive defined";
    if (covStatus == 3)
       covStatusType = "accurate";
+   if (covStatus == 0)
+      covStatusType = "full but not positive defined";
 
    if (!fState.HasCovariance()) {
       // if false means error is not valid and this is due to a failure in Hesse

math/minuit2/src/MnHesse.cxx

@@ -334,6 +334,7 @@ MinimumState MnHesse::ComputeNumerical(const MnFcn &mfcn, const MinimumState &st
 
    // off-diagonal Elements
    // initial starting values
+   bool doCentralFD = fStrategy.HessianCentralFDMixedDerivatives();
    if (n > 0) {
       MPIProcess mpiprocOffDiagonal(n * (n - 1) / 2, 0);
       unsigned int startParIndexOffDiagonal = mpiprocOffDiagonal.StartElementIndex();
@@ -357,10 +358,19 @@ MinimumState MnHesse::ComputeNumerical(const MnFcn &mfcn, const MinimumState &st
          x(j) += dirin(j);
 
          double fs1 = mfcn(x);
-         double elem = (fs1 + amin - yy(i) - yy(j)) / (dirin(i) * dirin(j));
-         vhmat(i, j) = elem;
-
-         x(j) -= dirin(j);
+         if(!doCentralFD) {
+            double elem = (fs1 + amin - yy(i) - yy(j)) / (dirin(i) * dirin(j));
+            vhmat(i, j) = elem;
+            x(j) -= dirin(j);
+         } else {
+            // three more function evaluations required for central fd
+            x(i) -= dirin(i); x(i) -= dirin(i);double fs3 = mfcn(x);
+            x(j) -= dirin(j); x(j) -= dirin(j);double fs4 = mfcn(x);
+            x(i) += dirin(i); x(i) += dirin(i);double fs2 = mfcn(x);
+            x(j) += dirin(j);
+            double elem = (fs1 - fs2 - fs3 + fs4)/(4.*dirin(i)*dirin(j));
+            vhmat(i, j) = elem;
+         }
 
          if (j % (n - 1) == 0 || in == endParIndexOffDiagonal - 1)
             x(i) -= dirin(i);
@@ -370,11 +380,17 @@ MinimumState MnHesse::ComputeNumerical(const MnFcn &mfcn, const MinimumState &st
    }
 
    // verify if matrix pos-def (still 2nd derivative)
+   // Note that for cases of extreme spread of eigenvalues, numerical precision
+   // can mean the hessian is computed as being not pos-def
+   // but the inverse of it is.
 
    print.Debug("Original error matrix", vhmat);
 
-   MinimumError tmpErr = MnPosDef()(MinimumError(vhmat, 1.), prec);
-   vhmat = tmpErr.InvHessian();
+   MinimumError tmpErr = MnPosDef()(MinimumError(vhmat, 1.), prec); // pos-def version of hessian
+
+   if(fStrategy.HessianForcePosDef()) {
+      vhmat = tmpErr.InvHessian();
+   }
 
    print.Debug("PosDef error matrix", vhmat);
 
@@ -398,7 +414,7 @@ MinimumState MnHesse::ComputeNumerical(const MnFcn &mfcn, const MinimumState &st
 
    // if matrix is made pos def returns anyway edm
    if (tmpErr.IsMadePosDef()) {
-      MinimumError err(vhmat, MinimumError::MnMadePosDef);
+      MinimumError err(vhmat, fStrategy.HessianForcePosDef() ? MinimumError::MnMadePosDef : MinimumError::MnNotPosDef);
       double edm = estim.Estimate(gr, err);
       return MinimumState(st.Parameters(), err, gr, edm, mfcn.NumOfCalls());
    }

math/minuit2/src/MnPosDef.cxx

@@ -43,7 +43,7 @@ MinimumError MnPosDef::operator()(const MinimumError &e, const MnMachinePrecisio
    }
    //   std::cout<<"MnPosDef init matrix= "<<err<<std::endl;
 
-   double epspdf = std::max(1.e-6, prec.Eps2());
+   double epspdf = std::max(1.e-12, prec.Eps2()); // should this lower bound be configurable?
    double dgmin = err(0, 0);
 
    for (unsigned int i = 0; i < err.Nrow(); i++) {

math/minuit2/src/MnStrategy.cxx

@@ -13,21 +13,23 @@ namespace ROOT {
 
 namespace Minuit2 {
 
-MnStrategy::MnStrategy() : fStoreLevel(1)
+MnStrategy::MnStrategy() : fHessCFDG2(0), fHessForcePosDef(1), fStoreLevel(1)
 {
    // default strategy
    SetMediumStrategy();
 }
 
-MnStrategy::MnStrategy(unsigned int stra) : fStoreLevel(1)
+MnStrategy::MnStrategy(unsigned int stra) : fHessCFDG2(0), fHessForcePosDef(1), fStoreLevel(1)
 {
-   // user defined strategy (0, 1, >=2)
+   // user defined strategy (0, 1, 2, >=3)
    if (stra == 0)
       SetLowStrategy();
    else if (stra == 1)
       SetMediumStrategy();
-   else
+   else if (stra == 2)
       SetHighStrategy();
+   else
+      SetVeryHighStrategy();
 }
 
 void MnStrategy::SetLowStrategy()
@@ -41,6 +43,7 @@ void MnStrategy::SetLowStrategy()
    SetHessianStepTolerance(0.5);
    SetHessianG2Tolerance(0.1);
    SetHessianGradientNCycles(1);
+   SetHessianCentralFDMixedDerivatives(0);
 }
 
 void MnStrategy::SetMediumStrategy()
@@ -54,6 +57,7 @@ void MnStrategy::SetMediumStrategy()
    SetHessianStepTolerance(0.3);
    SetHessianG2Tolerance(0.05);
    SetHessianGradientNCycles(2);
+   SetHessianCentralFDMixedDerivatives(0);
 }
 
 void MnStrategy::SetHighStrategy()
@@ -67,6 +71,22 @@ void MnStrategy::SetHighStrategy()
    SetHessianStepTolerance(0.1);
    SetHessianG2Tolerance(0.02);
    SetHessianGradientNCycles(6);
+   SetHessianCentralFDMixedDerivatives(0);
+}
+
+void MnStrategy::SetVeryHighStrategy()
+{
+    // set very high strategy (3)
+    fStrategy = 3;
+    SetGradientNCycles(5);
+    SetGradientStepTolerance(0.1);
+    SetGradientTolerance(0.02);
+    SetHessianNCycles(7);
+    SetHessianStepTolerance(0.);
+    SetHessianG2Tolerance(0.);
+    SetHessianGradientNCycles(6);
+    SetHessianCentralFDMixedDerivatives(1);
+    SetHessianForcePosDef(0);
 }
 
 } // namespace Minuit2

math/minuit2/src/MnUserParameterState.cxx

@@ -176,12 +176,14 @@ MnUserParameterState::MnUserParameterState(const MinimumState &st, double up, co
 
       assert(fCovariance.Nrow() == VariableParameters());
 
-      fCovStatus = 1; // when is valid
+      if (!st.Error().IsNotPosDef())
+         fCovStatus = 1; // when is valid and not NotPosDef
    }
    if (st.Error().IsMadePosDef())
       fCovStatus = 2;
    if (st.Error().IsAccurate())
       fCovStatus = 3;
+
 }
 
 MnUserCovariance MnUserParameterState::Hessian() const