[RF] Faster Hesse in RooFit by advertising which params are independent

This reduces the time to run Hesse in the ATLAS Higgs benchmark from
123 s to 92 seconds.

Given that some models take hours for this, this is a significant
improvement for the user experience.

Further improvement is possible by analyzing the computation graph a bit
more to find more independent parameters (e.g., the different gammas for
stat uncertainties from different bins).

Author: guitargeek

math/mathcore/inc/Math/Functor.h

@@ -67,6 +67,10 @@ class Functor : public IBaseFunctionMultiDim  {
    // for multi-dimensional functions
    unsigned int NDim() const override { return fDim; }
 
+   bool VanishingSecondDerivative(int i, int j) const override { return fVanishSecondDerivFunc ? fVanishSecondDerivFunc(i, j) : false; }
+
+   void SetVanishingSecondDerivativeFunc(std::function<bool(int, int)> func) { fVanishSecondDerivFunc = std::move(func); }
+
 private :
 
    inline double DoEval (const double * x) const override {
@@ -75,6 +79,7 @@ private :
 
    unsigned int fDim;
    std::function<double(double const *)> fFunc;
+   std::function<bool(int, int)> fVanishSecondDerivFunc;
 };
 
 /**
@@ -222,6 +227,10 @@ class GradFunctor : public IGradientFunctionMultiDim  {
       fGradFunc(x, g);
    }
 
+   bool VanishingSecondDerivative(int i, int j) const override { return fVanishSecondDerivFunc ? fVanishSecondDerivFunc(i, j) : false; }
+
+   void SetVanishingSecondDerivativeFunc(std::function<bool(int, int)> func) { fVanishSecondDerivFunc = std::move(func); }
+
 private :
 
    inline double DoEval (const double * x) const override {
@@ -244,6 +253,7 @@ private :
    std::function<double(const double *)> fFunc;
    std::function<double(double const *, unsigned int)> fDerivFunc;
    std::function<void(const double *, double*)> fGradFunc;
+   std::function<bool(int, int)> fVanishSecondDerivFunc;
 };
 
 

math/mathcore/inc/Math/IFunction.h

@@ -91,6 +91,11 @@ namespace ROOT {
          // if it inherits from ROOT::Math::IGradientFunctionMultiDim.
          virtual bool HasGradient() const { return false; }
 
+         /// Indicate whether a given second order derivative with respect to
+         /// parameters i and j is always zero. This can help to avoid
+         /// expensive function calls in Hessian evaluations.
+         virtual bool VanishingSecondDerivative(int /*i*/, int /*j*/) const { return false; }
+
       private:
 
          /// Implementation of the evaluation function. Must be implemented by derived classes.

math/minuit2/inc/Minuit2/FCNAdapter.h

@@ -43,6 +43,8 @@ class FCNAdapter : public FCNBase {
 
    void SetErrorDef(double up) override { fUp = up; }
 
+   bool VanishingSecondDerivative(int i, int j) const override { return fFunc.VanishingSecondDerivative(i, j); }
+
 private:
    const Function &fFunc;
    double fUp;

math/minuit2/inc/Minuit2/FCNBase.h

@@ -132,6 +132,11 @@ class FCNBase : public GenericFunction {
    virtual bool HasHessian() const { return false; }
 
    virtual bool HasG2() const { return false; }
+
+   /// Indicate whether a given second order derivative with respect to
+   /// parameters i and j is always zero. This can help to avoid
+   /// expensive function calls in Hessian evaluations.
+   virtual bool VanishingSecondDerivative(int /*i*/, int /*j*/) const { return false; }
 };
 
 } // namespace Minuit2

math/minuit2/inc/Minuit2/FCNGradAdapter.h

@@ -119,6 +119,8 @@ class FCNGradAdapter : public FCNBase {
 
    void SetErrorDef(double up) override { fUp = up; }
 
+   bool VanishingSecondDerivative(int i, int j) const { return fFunc.VanishingSecondDerivative(i, j); }
+
 private:
    const Function &fFunc;
    double fUp;

math/minuit2/src/MnHesse.cxx

@@ -383,14 +383,17 @@ MinimumState ComputeNumerical(const MnFcn &mfcn, const MinimumState &st, const M
          if ((i + 1) == j || in == startParIndexOffDiagonal)
             x(i) += dirin(i);
 
-         x(j) += dirin(j);
-
-         double fs1 = mfcnCaller(x);
-         if(!doCentralFD) {
+         if(mfcn.Fcn().VanishingSecondDerivative(i, j)) {
+            vhmat(i, j) = 0.;
+         } else if(!doCentralFD) {
+            x(j) += dirin(j);
+            double fs1 = mfcnCaller(x);
             double elem = (fs1 + amin - yy(i) - yy(j)) / (dirin(i) * dirin(j));
             vhmat(i, j) = elem;
             x(j) -= dirin(j);
          } else {
+            x(j) += dirin(j);
+            double fs1 = mfcnCaller(x);
             // three more function evaluations required for central fd
             x(i) -= dirin(i);
             x(i) -= dirin(i);

roofit/roofitcore/CMakeLists.txt

@@ -136,6 +136,7 @@ ROOT_STANDARD_LIBRARY_PACKAGE(RooFitCore
     RooFit/TestStatistics/RooUnbinnedL.h
     RooFit/TestStatistics/SharedOffset.h
     RooFit/TestStatistics/buildLikelihood.h
+    RooFit/VariableGroups.h
     RooFitLegacy/RooCatTypeLegacy.h
     RooFitLegacy/RooCategorySharedProperties.h
     RooFitLegacy/RooTreeData.h

roofit/roofitcore/inc/RooAbsArg.h

@@ -57,7 +57,8 @@ namespace RooFit {
 namespace Experimental {
 class CodegenContext;
 }
-}
+struct VariableGroups;
+} // namespace RooFit
 
 class RooRefArray : public TObjArray {
  public:
@@ -248,7 +249,7 @@ class RooAbsArg : public TNamed, public RooPrintable {
   bool recursiveCheckObservables(const RooArgSet* nset) const ;
   RooFit::OwningPtr<RooArgSet> getComponents() const ;
 
-
+  virtual void fillVariableGroups(RooFit::VariableGroups &out) const;
 
   void attachArgs(const RooAbsCollection &set);
   void attachDataSet(const RooAbsData &set);

roofit/roofitcore/inc/RooAddition.h

@@ -56,6 +56,8 @@ class RooAddition : public RooAbsReal {
 
   void doEval(RooFit::EvalContext &) const override;
 
+  void fillVariableGroups(RooFit::VariableGroups &out) const override;
+
 protected:
 
   RooArgList   _ownedList ;      ///< List of owned components

roofit/roofitcore/inc/RooConstraintSum.h

@@ -45,6 +45,8 @@ class RooConstraintSum : public RooAbsReal {
 
   std::unique_ptr<RooAbsArg> compileForNormSet(RooArgSet const &normSet, RooFit::Detail::CompileContext & ctx) const override;
 
+  void fillVariableGroups(RooFit::VariableGroups &out) const override;
+
 protected:
 
   RooListProxy _set1 ;    ///< Set of constraint terms

roofit/roofitcore/inc/RooFit/VariableGroups.h

@@ -0,0 +1,33 @@
+#ifndef RooFit_VariableGroups_h
+#define RooFit_VariableGroups_h
+
+#include <TNamed.h>
+
+#include <iostream>
+#include <unordered_map>
+#include <vector>
+
+class TNamed;
+
+namespace RooFit {
+
+struct VariableGroups {
+
+   std::unordered_map<TNamed const*, std::vector<int>> groups;
+
+   inline void print() {
+      for (auto const& item : groups) {
+          std::cout << item.first->GetName() << " :";
+          for (int n : item.second) {
+             std::cout << " " << n;
+          }
+          std::cout << std::endl;
+      }
+   }
+
+   int currentIndex = 0;
+};
+
+}
+
+#endif

roofit/roofitcore/src/RooAbsArg.cxx

@@ -75,8 +75,8 @@ for single nodes.
 #include <RooArgSet.h>
 #include <RooConstVar.h>
 #include <RooExpensiveObjectCache.h>
+#include <RooFit/VariableGroups.h>
 #include <RooHelpers.h>
-#include "RooFitImplHelpers.h"
 #include <RooListProxy.h>
 #include <RooMsgService.h>
 #include <RooRealIntegral.h>
@@ -87,6 +87,8 @@ for single nodes.
 #include <RooVectorDataStore.h>
 #include <RooWorkspace.h>
 
+#include "RooFitImplHelpers.h"
+
 #include <TBuffer.h>
 #include <TClass.h>
 #include <TVirtualStreamerInfo.h>
@@ -2529,3 +2531,19 @@ void RooAbsArg::setDataToken(std::size_t index)
    }
    _dataToken = index;
 }
+
+void RooAbsArg::fillVariableGroups(RooFit::VariableGroups &out) const
+{
+   // Get the set of nodes in the computation graph. Do the detour via
+   // RooArgList to avoid deduplication done after adding each element.
+   RooArgSet serverSet;
+   RooArgList serverList;
+   treeNodeServerList(&serverList, nullptr, /*branches*/ false, /*leaves*/ true, /*valueOnly*/ false,
+                           /*recurseFundamental*/ true);
+   serverSet.add(serverList.begin(), serverList.end());
+
+   for (RooAbsArg const *arg : serverSet) {
+      out.groups[arg->namePtr()].push_back(out.currentIndex);
+   }
+   out.currentIndex++;
+}

roofit/roofitcore/src/RooAddition.cxx

@@ -297,3 +297,10 @@ std::list<double>* RooAddition::plotSamplingHint(RooAbsRealLValue& obs, double x
 {
   return RooRealSumPdf::plotSamplingHint(_set, obs, xlo, xhi);
 }
+
+void RooAddition::fillVariableGroups(RooFit::VariableGroups &out) const
+{
+   for (RooAbsArg *arg : _set) {
+      arg->fillVariableGroups(out);
+   }
+}

roofit/roofitcore/src/RooConstraintSum.cxx

@@ -114,3 +114,10 @@ bool RooConstraintSum::setData(RooAbsData const& data, bool /*cloneData=true*/)
   }
   return true;
 }
+
+void RooConstraintSum::fillVariableGroups(RooFit::VariableGroups &out) const
+{
+   for (RooAbsArg *arg : _set1) {
+      arg->fillVariableGroups(out);
+   }
+}

roofit/roofitcore/src/RooEvaluatorWrapper.cxx

@@ -136,4 +136,9 @@ bool RooEvaluatorWrapper::setData(RooAbsData &data, bool /*cloneData*/)
    return true;
 }
 
+void RooEvaluatorWrapper::fillVariableGroups(RooFit::VariableGroups &out) const
+{
+   _topNode->fillVariableGroups(out);
+}
+
 /// \endcond

roofit/roofitcore/src/RooEvaluatorWrapper.h

@@ -61,6 +61,8 @@ class RooEvaluatorWrapper final : public RooAbsReal {
    /// The RooFit::Evaluator is dealing with constant terms itself.
    void constOptimizeTestStatistic(ConstOpCode /*opcode*/, bool /*doAlsoTrackingOpt*/) override {}
 
+   void fillVariableGroups(RooFit::VariableGroups &out) const override;
+
 protected:
    double evaluate() const override;
 

roofit/roofitcore/src/RooMinimizer.cxx

@@ -134,7 +134,7 @@ RooMinimizer::RooMinimizer(RooAbsReal &function, Config const &cfg) : _cfg(cfg)
       _fcn = std::make_unique<RooMinimizerFcn>(&function, this);
    }
    initMinimizerFcnDependentPart(function.defaultErrorLevel());
-};
+}
 
 /// Initialize the part of the minimizer that is independent of the function to be minimized
 void RooMinimizer::initMinimizerFirstPart()

roofit/roofitcore/src/RooMinimizerFcn.cxx

@@ -25,10 +25,11 @@
 #include "RooAbsArg.h"
 #include "RooAbsPdf.h"
 #include "RooArgSet.h"
-#include "RooRealVar.h"
-#include "RooMsgService.h"
+#include "RooFit/VariableGroups.h"
 #include "RooMinimizer.h"
+#include "RooMsgService.h"
 #include "RooNaNPacker.h"
+#include "RooRealVar.h"
 
 #include "Math/Functor.h"
 #include "TMatrixDSym.h"
@@ -40,6 +41,23 @@ using std::setprecision;
 
 namespace {
 
+template <class InputIt1, class InputIt2>
+bool intersect(InputIt1 first1, InputIt1 last1, InputIt2 first2, InputIt2 last2)
+{
+   while (first1 != last1 && first2 != last2) {
+      if (*first1 < *first2) {
+         ++first1;
+         continue;
+      }
+      if (*first2 < *first1) {
+         ++first2;
+         continue;
+      }
+      return true;
+   }
+   return false;
+}
+
 // Helper function that wraps RooAbsArg::getParameters and directly returns the
 // output RooArgSet. To be used in the initializer list of the RooMinimizerFcn
 // constructor.
@@ -56,14 +74,38 @@ RooArgSet getParameters(RooAbsReal const &funct)
 RooMinimizerFcn::RooMinimizerFcn(RooAbsReal *funct, RooMinimizer *context)
    : RooAbsMinimizerFcn(getParameters(*funct), context), _funct(funct)
 {
-   unsigned int nDim = getNDim();
+   RooFit::VariableGroups groups;
+   funct->fillVariableGroups(groups);
+
+   RooArgList parameters;
+   for (std::size_t i = 0; i < getNDim(); ++i) {
+      parameters.add(floatableParam(i));
+   }
+
+   std::size_t nParams = parameters.size();
+
+   _secondDerivMask.resize(nParams * nParams);
+   for (std::size_t i = 0; i < nParams; ++i) {
+      _secondDerivMask[nParams * i + i] = 1;
+      for (std::size_t j = 0; j < i; ++j) {
+         // std::cout << parameters[i].GetName() << "   " << parameters[j].GetName() << std::endl;
+         auto const &gr1 = groups.groups.at(parameters[i].namePtr());
+         auto const &gr2 = groups.groups.at(parameters[j].namePtr());
+         _secondDerivMask[nParams * i + j] = intersect(gr1.begin(), gr1.end(), gr2.begin(), gr2.end());
+         _secondDerivMask[nParams * j + i] = _secondDerivMask[nParams * i + j];
+      }
+   }
 
    if (context->_cfg.useGradient && funct->hasGradient()) {
       _gradientOutput.resize(_allParams.size());
-      _multiGenFcn = std::make_unique<ROOT::Math::GradFunctor>(this, &RooMinimizerFcn::operator(),
-                                                               &RooMinimizerFcn::evaluateGradient, nDim);
+      auto functor = std::make_unique<ROOT::Math::GradFunctor>(this, &RooMinimizerFcn::operator(),
+                                                               &RooMinimizerFcn::evaluateGradient, getNDim());
+      functor->SetVanishingSecondDerivativeFunc([this](int i, int j) { return this->vanishingSecondDerivative(i, j); });
+      _multiGenFcn = std::move(functor);
    } else {
-      _multiGenFcn = std::make_unique<ROOT::Math::Functor>(std::cref(*this), getNDim());
+      auto functor = std::make_unique<ROOT::Math::Functor>(std::cref(*this), getNDim());
+      functor->SetVanishingSecondDerivativeFunc([this](int i, int j) { return this->vanishingSecondDerivative(i, j); });
+      _multiGenFcn = std::move(functor);
    }
 }
 
@@ -148,4 +190,9 @@ void RooMinimizerFcn::setOffsetting(bool flag)
    _funct->enableOffsetting(flag);
 }
 
+bool RooMinimizerFcn::vanishingSecondDerivative(int i, int j) const
+{
+   return _secondDerivMask[getNDim() * i + j] == 0;
+}
+
 /// \endcond

roofit/roofitcore/src/RooMinimizerFcn.h

@@ -46,10 +46,13 @@ class RooMinimizerFcn : public RooAbsMinimizerFcn {
    double operator()(const double *x) const;
    void evaluateGradient(const double *x, double *out) const;
 
+   bool vanishingSecondDerivative(int i, int j) const;
+
 private:
    RooAbsReal *_funct = nullptr;
    std::unique_ptr<ROOT::Math::IBaseFunctionMultiDim> _multiGenFcn;
    mutable std::vector<double> _gradientOutput;
+   std::vector<int> _secondDerivMask;
 };
 
 #endif

roofit/roofitcore/src/TestStatistics/MinuitFcnGrad.cxx

@@ -43,6 +43,9 @@ class MinuitGradFunctor : public ROOT::Math::IMultiGradFunction {
 
    bool returnsInMinuit2ParameterSpace() const override { return _fcn.returnsInMinuit2ParameterSpace(); }
 
+   // TODO: Implement this
+   bool VanishingSecondDerivative(int /*i*/, int /*j*/) const override { return false; }
+
 private:
    double DoEval(const double *x) const override { return _fcn(x); }