From a32e4929b73681acd9e7b7678c9f9328749e488a Mon Sep 17 00:00:00 2001 From: Jonas Rembser Date: Sun, 6 Sep 2026 21:57:38 +0000 Subject: [PATCH 1/2] [Minuit2] Skip identically-vanishing second derivatives in MnHesse Objective functions can now advertise pairs of parameters whose mixed second derivative is identically zero, via the new virtual function FCNBase::SecondDerivativeAlwaysVanishes(). The numerical Hessian computation in MnHesse skips the finite-difference evaluations for such parameter pairs, which can speed up Hesse significantly for likelihoods with many mutually independent parameters. The parameter indices in this interface refer to the FCN's own full (external) parameter space, and the advertised information must hold for all parameter values; see the FCNBase documentation for the complete contract. Since MnHesse loops over Minuit-internal indices that exclude fixed parameters, it translates them via MnUserTransformation::ExtOfInt() before querying the predicate, so results stay correct when parameters are fixed in the minimizer. Clients that drive Minuit2 through the ROOT::Math interfaces can inject the predicate with the new Minuit2Minimizer::SetSecondDerivativeAlwaysVanishesFunc(), which follows the same pattern as SetHessianFunction(). This keeps the feature out of the general ROOT::Math function interfaces on purpose: it is only consumed by Minuit2. --- README/ReleaseNotes/v642/index.md | 6 + math/minuit2/inc/Minuit2/FCNAdapter.h | 23 ++++ math/minuit2/inc/Minuit2/FCNBase.h | 20 +++ math/minuit2/inc/Minuit2/Minuit2Minimizer.h | 9 ++ math/minuit2/src/Minuit2Minimizer.cxx | 9 ++ math/minuit2/src/MnHesse.cxx | 14 +- math/minuit2/test/testMinuit2.cxx | 134 ++++++++++++++++++++ 7 files changed, 211 insertions(+), 4 deletions(-) diff --git a/README/ReleaseNotes/v642/index.md b/README/ReleaseNotes/v642/index.md index bcfa03a21d92f..a02d965e609c7 100644 --- a/README/ReleaseNotes/v642/index.md +++ b/README/ReleaseNotes/v642/index.md @@ -122,6 +122,12 @@ maps) will now obtain different, mathematically consistent values. ## Math +### Skipping identically-vanishing second derivatives in numerical Hessian evaluation + +Minuit2 objective functions can now advertise pairs of parameters whose mixed second derivative is identically zero, via the new virtual function `ROOT::Minuit2::FCNBase::SecondDerivativeAlwaysVanishes()` (with corresponding setters on `ROOT::Minuit2::FCNAdapter` and `ROOT::Minuit2::Minuit2Minimizer`). +The numerical Hessian computation in Minuit2 (`MnHesse`) skips the finite-difference evaluations for such parameter pairs, which can speed up Hesse significantly for likelihoods with many mutually independent parameters. +The parameter indices in this interface always refer to the function's own full (external) parameter space, including parameters that are fixed in the minimizer, and the advertised information must hold for all parameter values; see the `FCNBase` documentation for the full contract. + ## RDataFrame * Added `RedefinePerSample` transformation. Works similarly to `DefinePerSample`, but allows to redefine existing values diff --git a/math/minuit2/inc/Minuit2/FCNAdapter.h b/math/minuit2/inc/Minuit2/FCNAdapter.h index d1725a6d0a30b..1b2c87bead441 100644 --- a/math/minuit2/inc/Minuit2/FCNAdapter.h +++ b/math/minuit2/inc/Minuit2/FCNAdapter.h @@ -160,6 +160,27 @@ class FCNAdapter : public FCNBase { /// @param up New error definition value. void SetErrorDef(double up) override { fUp = up; } + /// Indicate whether the mixed second order derivative with respect to + /// parameters i and j is identically zero, forwarding to the user-provided + /// predicate if set. + bool SecondDerivativeAlwaysVanishes(unsigned int i, unsigned int j) const override + { + return fSecondDerivAlwaysVanishesFunc ? fSecondDerivAlwaysVanishesFunc(i, j) : false; + } + + /// Set the predicate advertising which mixed second derivatives are + /// identically zero. + /// + /// @param f Function taking two parameter indices (in the FCN's full + /// external parameter space) and returning `true` if the + /// corresponding second derivative is zero for all parameter + /// values. It must fulfill the contract documented for + /// FCNBase::SecondDerivativeAlwaysVanishes(). + void SetSecondDerivativeAlwaysVanishesFunc(std::function f) + { + fSecondDerivAlwaysVanishesFunc = std::move(f); + } + private: using Function = std::function; using GradFunction = std::function; @@ -173,6 +194,8 @@ class FCNAdapter : public FCNBase { GradFunction fGradFunc; ///< Optional gradient function. G2Function fG2Func; ///< Optional diagonal second-derivative function. mutable HessianFunction fHessianFunc; ///< Optional Hessian function. + /// Optional predicate for identically-vanishing second derivatives. + std::function fSecondDerivAlwaysVanishesFunc; }; } // namespace ROOT::Minuit2 diff --git a/math/minuit2/inc/Minuit2/FCNBase.h b/math/minuit2/inc/Minuit2/FCNBase.h index 245c02ac05a6c..f9798404d1152 100644 --- a/math/minuit2/inc/Minuit2/FCNBase.h +++ b/math/minuit2/inc/Minuit2/FCNBase.h @@ -127,6 +127,26 @@ class FCNBase { virtual bool HasHessian() const { return false; } virtual bool HasG2() const { return false; } + + /// Indicate whether the mixed second order derivative with respect to + /// parameters i and j is identically zero, i.e. zero for *all* parameter + /// values. This can help to avoid expensive function calls in numerical + /// Hessian evaluations (see MnHesse). + /// + /// The contract for implementations: + /// + /// * The indices refer to the FCN's own full (external) parameter + /// space, including any parameters that are fixed in the minimizer. + /// Callers that work with internal indices, like MnHesse, must + /// translate them via MnUserTransformation::ExtOfInt() before calling + /// this function. + /// * The result must be symmetric in i and j. + /// * Only return `true` if the second derivative vanishes for *any* + /// value of the parameters: the result may be cached and is used at + /// arbitrary points in parameter space. + /// * The diagonal (i == j) is never queried by Minuit2, so its return + /// value has no effect. + virtual bool SecondDerivativeAlwaysVanishes(unsigned int /*i*/, unsigned int /*j*/) const { return false; } }; } // namespace ROOT::Minuit2 diff --git a/math/minuit2/inc/Minuit2/Minuit2Minimizer.h b/math/minuit2/inc/Minuit2/Minuit2Minimizer.h index dd1669a4f1fd4..c188f6038df53 100644 --- a/math/minuit2/inc/Minuit2/Minuit2Minimizer.h +++ b/math/minuit2/inc/Minuit2/Minuit2Minimizer.h @@ -82,6 +82,15 @@ class Minuit2Minimizer : public ROOT::Math::Minimizer { /// set the function implementing Hessian computation void SetHessianFunction(std::function, double *)> hfunc) override; + /// Set a predicate advertising which mixed second derivatives of the + /// minimized function are identically zero, so that MnHesse can skip the + /// corresponding finite-difference evaluations. The predicate must + /// fulfill the contract documented for + /// FCNBase::SecondDerivativeAlwaysVanishes(); in particular, the indices + /// refer to the full (external) parameter space. Like SetHessianFunction, + /// this must be called after SetFunction. + void SetSecondDerivativeAlwaysVanishesFunc(std::function func); + /// set free variable bool SetVariable(unsigned int ivar, const std::string &name, double val, double step) override; diff --git a/math/minuit2/src/Minuit2Minimizer.cxx b/math/minuit2/src/Minuit2Minimizer.cxx index 9788394f02d09..54b87884ae7b1 100644 --- a/math/minuit2/src/Minuit2Minimizer.cxx +++ b/math/minuit2/src/Minuit2Minimizer.cxx @@ -416,6 +416,15 @@ void Minuit2Minimizer::SetHessianFunction(std::functionSetHessianFunction(hfunc); } +void Minuit2Minimizer::SetSecondDerivativeAlwaysVanishesFunc(std::function func) +{ + // not supported for Fumili, whose FCN is not an FCNAdapter + if (fUseFumili) return; + auto fcn = static_cast(fMinuitFCN.get()); + if (!fcn) return; + fcn->SetSecondDerivativeAlwaysVanishesFunc(std::move(func)); +} + namespace { ROOT::Minuit2::MnStrategy customizedStrategy(unsigned int strategyLevel, ROOT::Math::MinimizerOptions const &options) diff --git a/math/minuit2/src/MnHesse.cxx b/math/minuit2/src/MnHesse.cxx index a5fa41edc3c61..2f9591c15fc77 100644 --- a/math/minuit2/src/MnHesse.cxx +++ b/math/minuit2/src/MnHesse.cxx @@ -342,14 +342,20 @@ 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) { + // The FCN advertises vanishing second derivatives in its own + // (external) parameter space, but i and j are internal indices that + // exclude fixed parameters, so translate before asking. + if(mfcn.Fcn().SecondDerivativeAlwaysVanishes(trafo.ExtOfInt(i), trafo.ExtOfInt(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); diff --git a/math/minuit2/test/testMinuit2.cxx b/math/minuit2/test/testMinuit2.cxx index 7f3624f15e74c..f33bb0e7494ad 100644 --- a/math/minuit2/test/testMinuit2.cxx +++ b/math/minuit2/test/testMinuit2.cxx @@ -3,6 +3,7 @@ // Author: Jonas Rembser, CERN 01/2026 #include +#include #include #include #include @@ -411,3 +412,136 @@ TEST(Minuit2, AnalyticalHessianLimitTransformation) } } } + +// ---------------------------------------------------------------------- +// Tests for FCNBase::SecondDerivativeAlwaysVanishes() +// ---------------------------------------------------------------------- + +// Sum of two additive terms: the first couples parameters (0, 2) and the +// second couples (1, 2), so only the mixed second derivative with respect to +// the pair (0, 1) is identically zero. The full Hessian is: +// +// [ 4 0 -2] +// [ 0 4 -2] +// [-2 -2 6] +class SparseQuadraticFCN : public ROOT::Minuit2::FCNBase { +public: + double operator()(const std::vector &p) const override + { + ++fNCalls; + const double a = p[0]; + const double b = p[1]; + const double c = p[2]; + // clang-format off + return (a-1)*(a-1) + (a-c)*(a-c) + (c-1)*(c-1) + (b-c)*(b-c) + (b-1)*(b-1); + // clang-format on + } + + double Up() const override { return 1.0; } + + mutable int fNCalls = 0; +}; + +// Same function, but advertising the identically-vanishing mixed second +// derivative with respect to the (external) parameter pair (0, 1). +class SparseQuadraticFCNWithInfo : public SparseQuadraticFCN { +public: + bool SecondDerivativeAlwaysVanishes(unsigned int i, unsigned int j) const override + { + return (i == 0 && j == 1) || (i == 1 && j == 0); + } +}; + +namespace { + +/// Minimize, run MnHesse explicitly, and report the number of FCN calls that +/// the MnHesse run needed in the output parameter. +template +ROOT::Minuit2::FunctionMinimum RunMigradAndHesse(const FCN &fcn, int &nHesseCalls) +{ + using namespace ROOT::Minuit2; + + MnUserParameters upar; + upar.Add("a", 0.5, 0.1); + upar.Add("b", 0.5, 0.1); + upar.Add("c", 0.5, 0.1); + + MnMigrad migrad(fcn, upar); + FunctionMinimum min = migrad(); + + fcn.fNCalls = 0; + MnHesse hesse; + hesse(fcn, min); + nHesseCalls = fcn.fNCalls; + return min; +} + +} // namespace + +// The Hessian from MnHesse must be correct whether or not the FCN advertises +// its identically-vanishing second derivative, and advertising it must save +// function calls. +TEST(Minuit2, SecondDerivativeAlwaysVanishes) +{ + using ROOT::Minuit2::FunctionMinimum; + + const double tol = 1e-5; + + SparseQuadraticFCN fcnDense; + int nCallsDense = 0; + const FunctionMinimum minDense = RunMigradAndHesse(fcnDense, nCallsDense); + + SparseQuadraticFCNWithInfo fcnSparse; + int nCallsSparse = 0; + const FunctionMinimum minSparse = RunMigradAndHesse(fcnSparse, nCallsSparse); + + for (FunctionMinimum const *min : {&minDense, &minSparse}) { + const auto &hessian = min->Error().Hessian(); + ASSERT_EQ(hessian.Nrow(), 3); + EXPECT_NEAR(hessian(0, 0), 4.0, tol); + EXPECT_NEAR(hessian(1, 1), 4.0, tol); + EXPECT_NEAR(hessian(2, 2), 6.0, tol); + EXPECT_NEAR(hessian(0, 1), 0.0, tol); + EXPECT_NEAR(hessian(0, 2), -2.0, tol); + EXPECT_NEAR(hessian(1, 2), -2.0, tol); + } + + // The advertised element must be skipped, not computed: zero up to the + // round-off of the covariance inversion round-trip (much smaller than + // finite-difference noise), and with fewer function calls. + EXPECT_NEAR(minSparse.Error().Hessian()(0, 1), 0.0, 1e-12); + EXPECT_LT(nCallsSparse, nCallsDense); +} + +// With parameter "a" fixed, Minuit's internal parameter indices no longer +// coincide with the external ones that the FCN advertises the vanishing +// second derivative in. If MnHesse did not translate the indices via +// MnUserTransformation::ExtOfInt(), the internal pair (0, 1) = (b, c) would +// be looked up as the advertised external pair (0, 1) = (a, b), wrongly +// zeroing the genuinely non-vanishing (b, c) Hessian element. +TEST(Minuit2, SecondDerivativeAlwaysVanishesFixedParameter) +{ + using namespace ROOT::Minuit2; + + const double tol = 1e-5; + + SparseQuadraticFCNWithInfo fcn; + + MnUserParameters upar; + upar.Add("a", 1.0, 0.1); + upar.Add("b", 0.5, 0.1); + upar.Add("c", 0.5, 0.1); + upar.Fix("a"); + + MnMigrad migrad(fcn, upar); + FunctionMinimum min = migrad(); + MnHesse hesse; + hesse(fcn, min); + + // Reduced Hessian in the internal space (b, c) + const auto &hessian = min.Error().Hessian(); + ASSERT_EQ(hessian.Nrow(), 2); + EXPECT_NEAR(hessian(0, 0), 4.0, tol); + EXPECT_NEAR(hessian(0, 1), -2.0, tol); + EXPECT_NEAR(hessian(1, 1), 6.0, tol); +} From 418e5251fc04e11122fca6b8b4f9b71ffea31f69 Mon Sep 17 00:00:00 2001 From: Jonas Rembser Date: Sun, 6 Sep 2026 21:57:57 +0000 Subject: [PATCH 2/2] [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. RooFit analyzes the computation graph of the minimized function to find pairs of parameters that never appear in the same additive term of the likelihood, meaning their mixed second derivative is identically zero. The analysis is implemented privately in RooMinimizerFcn, on purpose without adding any new public interfaces, so that the design can still evolve: only strictly additive nodes (RooAddition, RooConstraintSum, and the RooEvaluatorWrapper forwarding to its wrapped function) are recursed into; everything else conservatively contributes all of its graph leaves as a single term, which advertises no independence but is always correct. The resulting pairwise mask is stored as a packed bitvector and built lazily on the first query, so that migrad-only fits and toy loops do not pay for it. It is handed to Minuit2 in RooMinimizerFcn::initMinimizer() via Minuit2Minimizer::SetSecondDerivativeAlwaysVanishesFunc(), which is why RooFitCore now has a private link dependency on Minuit2. The mask indices follow Minuit's external parameter convention, so the result is also correct when parameters are fixed after the RooMinimizer was constructed, which is covered by a unit test. 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). --- README/ReleaseNotes/v642/index.md | 7 + roofit/roofitcore/CMakeLists.txt | 4 + roofit/roofitcore/inc/RooConstraintSum.h | 2 +- roofit/roofitcore/inc/RooEvaluatorWrapper.h | 3 + roofit/roofitcore/src/RooMinimizer.cxx | 2 +- roofit/roofitcore/src/RooMinimizerFcn.cxx | 147 +++++++++++++++++- roofit/roofitcore/src/RooMinimizerFcn.h | 7 + .../src/TestStatistics/MinuitFcnGrad.cxx | 3 + roofit/roofitcore/test/testRooMinimizer.cxx | 77 +++++++++ 9 files changed, 247 insertions(+), 5 deletions(-) diff --git a/README/ReleaseNotes/v642/index.md b/README/ReleaseNotes/v642/index.md index a02d965e609c7..293187bde9380 100644 --- a/README/ReleaseNotes/v642/index.md +++ b/README/ReleaseNotes/v642/index.md @@ -151,6 +151,13 @@ the cut instead of being selected based on `sqrt(abs(x))`. * The `RooMinimizer::Strategy` enum has been removed. It named the Minuit strategies that are usually referred to just by integers, but caused confusion because it didn't include the unnamed "Strategy 3". Since people usually set the strategy with integer values anyway, it was decided that the simplest solution to avoid the confusion was simply to remove the `RooMinimizer::Strategy` enum +### Faster Hesse for likelihoods with many independent parameters + +RooFit now analyzes the computation graph of the minimized function to find pairs of parameters that never appear in the same additive term of the likelihood, meaning their mixed second derivative is identically zero. +This information is forwarded to Minuit2, which skips the corresponding finite-difference evaluations in the numerical Hessian computation (see the Math section above). +For likelihoods with many mutually independent parameters, such as the per-channel nuisance parameters of large combined HistFactory models, this can speed up `RooMinimizer::hesse()` by 30 % or more, with results identical up to floating-point noise. +This optimization is automatic and requires no user action. + ### Deprecation of the legacy evaluation backend The `legacy` evaluation backend for likelihood and chi-square fits is deprecated and will be removed in ROOT 6.44. diff --git a/roofit/roofitcore/CMakeLists.txt b/roofit/roofitcore/CMakeLists.txt index d8143917c8073..914accf2d899e 100644 --- a/roofit/roofitcore/CMakeLists.txt +++ b/roofit/roofitcore/CMakeLists.txt @@ -486,6 +486,10 @@ ROOT_STANDARD_LIBRARY_PACKAGE(RooFitCore ${EXTRA_DICT_OPTS} ) +# Minuit2 is needed in the implementation of RooMinimizerFcn, to advertise +# vanishing second derivatives to Minuit2Minimizer. +target_link_libraries(RooFitCore PRIVATE Minuit2) + # The following definitions are PUBLIC so they can also be used in ROOT-internal tests if(roofit_legacy_eval_backend) diff --git a/roofit/roofitcore/inc/RooConstraintSum.h b/roofit/roofitcore/inc/RooConstraintSum.h index f66d16f197524..b03ca80063bff 100644 --- a/roofit/roofitcore/inc/RooConstraintSum.h +++ b/roofit/roofitcore/inc/RooConstraintSum.h @@ -33,7 +33,7 @@ class RooConstraintSum : public RooAbsReal { RooConstraintSum(const RooConstraintSum& other, const char* name = nullptr); TObject* clone(const char* newname=nullptr) const override { return new RooConstraintSum(*this, newname); } - const RooArgList& list() { return _set1 ; } + const RooArgList& list() const { return _set1 ; } bool setData(RooAbsData const& data, bool cloneData=true); /// \copydoc setData(RooAbsData const&, bool) diff --git a/roofit/roofitcore/inc/RooEvaluatorWrapper.h b/roofit/roofitcore/inc/RooEvaluatorWrapper.h index 7b5fe148349ce..cd0ec26b7c769 100644 --- a/roofit/roofitcore/inc/RooEvaluatorWrapper.h +++ b/roofit/roofitcore/inc/RooEvaluatorWrapper.h @@ -70,6 +70,9 @@ class RooEvaluatorWrapper final : public RooAbsReal { RooFit::Evaluator &evaluator() const { return *_evaluator; } + /// The RooFit object that this wrapper evaluates. + RooAbsReal const &topNode() const { return *_topNode; } + protected: double evaluate() const override; diff --git a/roofit/roofitcore/src/RooMinimizer.cxx b/roofit/roofitcore/src/RooMinimizer.cxx index 6a62be7ef4725..2246a8f3d3529 100644 --- a/roofit/roofitcore/src/RooMinimizer.cxx +++ b/roofit/roofitcore/src/RooMinimizer.cxx @@ -205,7 +205,7 @@ RooMinimizer::RooMinimizer(RooAbsReal &function, Config const &cfg) : _function{ _fcn = std::make_unique(&function, this); } initMinimizerFcnDependentPart(function.defaultErrorLevel()); -}; +} /// Initialize the part of the minimizer that is independent of the function to be minimized void RooMinimizer::initMinimizerFirstPart() diff --git a/roofit/roofitcore/src/RooMinimizerFcn.cxx b/roofit/roofitcore/src/RooMinimizerFcn.cxx index 71fedcf48bcf8..6a3235950f2b3 100644 --- a/roofit/roofitcore/src/RooMinimizerFcn.cxx +++ b/roofit/roofitcore/src/RooMinimizerFcn.cxx @@ -24,23 +24,49 @@ #include "RooAbsArg.h" #include "RooAbsPdf.h" +#include "RooAddition.h" #include "RooArgSet.h" -#include "RooRealVar.h" -#include "RooMsgService.h" +#include "RooConstraintSum.h" +#include "RooEvaluatorWrapper.h" #include "RooMinimizer.h" +#include "RooMsgService.h" #include "RooNaNPacker.h" #include "RooCategory.h" +#include "RooRealVar.h" #include "Math/Functor.h" +#include "Minuit2/Minuit2Minimizer.h" #include "TMatrixDSym.h" #include #include +#include +#include using std::setprecision; namespace { +// Check whether two sorted ranges have at least one element in common. +// Like std::set_intersection, both input ranges must be sorted; the early +// return on the first match keeps the common case cheap. +template +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. In the case of figuring out all parameters for the minimizer, @@ -53,6 +79,70 @@ RooArgSet getAllParameters(RooAbsReal const &funct) return out; } +// Groups of computation-graph leaves by the additive term of the minimized +// function they appear in. Two parameters that share no term index have a +// mixed second derivative that is identically zero. +// +// Note that distinct objects with the same name share their pointer from the +// name registry, so they end up merged in the same map entry. This is +// conservative: it can only add co-occurrences, never remove any. +struct VariableGroups { + /// For each leaf (keyed by its unique name pointer), the sorted list of + /// indices of the additive terms it appears in. + std::unordered_map> groups; + + /// Register one additive term: record for every leaf in the collection + /// that it appears in this term. + void registerTerm(RooAbsCollection const &leaves) + { + for (RooAbsArg const *arg : leaves) { + groups[arg->namePtr()].push_back(_nextIndex); + } + ++_nextIndex; + } + +private: + int _nextIndex = 0; +}; + +// Fill the map from computation-graph leaves to the additive terms of the +// minimized function they appear in. +// +// Recursing into the components of a node is only correct if the value of +// the node is a *strictly additive* combination of them: recursing into +// anything else would wrongly advertise vanishing second derivatives and +// silently corrupt Hessian results. Any other node contributes all of its +// leaves as one single term, which advertises no independence but is always +// correct. +void fillVariableGroups(RooAbsArg const &arg, VariableGroups &out) +{ + if (auto addition = dynamic_cast(&arg)) { + for (RooAbsArg *component : addition->list()) { + fillVariableGroups(*component, out); + } + return; + } + if (auto constraintSum = dynamic_cast(&arg)) { + for (RooAbsArg *component : constraintSum->list()) { + fillVariableGroups(*component, out); + } + return; + } + if (auto wrapper = dynamic_cast(&arg)) { + fillVariableGroups(wrapper->topNode(), out); + return; + } + + // Get the set of leaves in the computation graph. Do the detour via + // RooArgList to avoid deduplication done after adding each element. + RooArgSet leafSet; + RooArgList leafList; + arg.treeNodeServerList(&leafList, nullptr, /*branches*/ false, /*leaves*/ true, /*valueOnly*/ false, + /*recurseFundamental*/ true); + leafSet.add(leafList.begin(), leafList.end()); + out.registerTerm(leafSet); +} + } // namespace // use reference wrapper for the Functor, such that the functor points to this RooMinimizerFcn by reference. @@ -66,7 +156,7 @@ RooMinimizerFcn::RooMinimizerFcn(RooAbsReal *funct, RooMinimizer *context) _multiGenFcn = std::make_unique(this, &RooMinimizerFcn::operator(), &RooMinimizerFcn::evaluateGradient, nDim); } else { - _multiGenFcn = std::make_unique(std::cref(*this), getNDim()); + _multiGenFcn = std::make_unique(std::cref(*this), nDim); } if (context->_cfg.useHessian) { _hessianOutput.resize(_allParams.size() * _allParams.size()); @@ -228,6 +318,57 @@ void RooMinimizerFcn::initMinimizer(ROOT::Math::Minimizer &minim, RooMinimizer * minim.SetHessianFunction( std::bind(&RooMinimizerFcn::evaluateHessian, this, std::placeholders::_1, std::placeholders::_2)); } + // The independence information for skipping vanishing second derivatives + // in numerical Hessian computations is a Minuit2-only feature, so it is + // wired up directly with the concrete minimizer type instead of going + // through the ROOT::Math::Minimizer interface. + if (auto *minuit2 = dynamic_cast(&minim)) { + minuit2->SetSecondDerivativeAlwaysVanishesFunc( + [this](unsigned int i, unsigned int j) { return secondDerivativeAlwaysVanishes(i, j); }); + } +} + +//////////////////////////////////////////////////////////////////////////////// +/// Fill the bitvector that flags for each pair of floatable parameters +/// whether they appear together in at least one additive term of the +/// minimized function, i.e. whether their mixed second derivative can be +/// non-vanishing. Built lazily because it is only needed for Hessian +/// evaluations, and building it for models with many parameters is not free. +void RooMinimizerFcn::buildSecondDerivMask() const +{ + VariableGroups groups; + fillVariableGroups(*_funct, groups); + + std::size_t nParams = getNDim(); + + // Packed bitvector: bit set means the parameter pair shares an additive + // term, so the mixed second derivative can be non-zero. + _secondDerivMask.assign(nParams * nParams, false); + for (std::size_t i = 0; i < nParams; ++i) { + _secondDerivMask[nParams * i + i] = true; + auto found1 = groups.groups.find(floatableParam(i).namePtr()); + for (std::size_t j = 0; j < i; ++j) { + auto found2 = groups.groups.find(floatableParam(j).namePtr()); + // A parameter that was not seen in the computation graph traversal + // is conservatively treated as intersecting with everything. + bool canBeNonZero = found1 == groups.groups.end() || found2 == groups.groups.end() || + intersect(found1->second.begin(), found1->second.end(), found2->second.begin(), + found2->second.end()); + _secondDerivMask[nParams * i + j] = canBeNonZero; + _secondDerivMask[nParams * j + i] = canBeNonZero; + } + } +} + +//////////////////////////////////////////////////////////////////////////////// +/// Report whether the second derivative with respect to parameters i and j +/// (indices in the space of all floatable parameters, matching Minuit's +/// external parameter indices) is identically zero because the parameters +/// share no additive term of the minimized function. +bool RooMinimizerFcn::secondDerivativeAlwaysVanishes(unsigned int i, unsigned int j) const +{ + std::call_once(_secondDerivMaskOnce, &RooMinimizerFcn::buildSecondDerivMask, this); + return !_secondDerivMask[getNDim() * i + j]; } /// \endcond diff --git a/roofit/roofitcore/src/RooMinimizerFcn.h b/roofit/roofitcore/src/RooMinimizerFcn.h index 48c3cad015046..de7052c161198 100644 --- a/roofit/roofitcore/src/RooMinimizerFcn.h +++ b/roofit/roofitcore/src/RooMinimizerFcn.h @@ -23,6 +23,7 @@ #include "RooArgList.h" #include +#include #include #include "RooAbsMinimizerFcn.h" @@ -48,11 +49,17 @@ class RooMinimizerFcn : public RooAbsMinimizerFcn { RooArgSet freezeDisconnectedParameters() const override; + bool secondDerivativeAlwaysVanishes(unsigned int i, unsigned int j) const; + private: + void buildSecondDerivMask() const; + RooAbsReal *_funct = nullptr; std::unique_ptr _multiGenFcn; mutable std::vector _gradientOutput; mutable std::vector _hessianOutput; + mutable std::vector _secondDerivMask; ///< Lazily built, see buildSecondDerivMask(). + mutable std::once_flag _secondDerivMaskOnce; }; #endif diff --git a/roofit/roofitcore/src/TestStatistics/MinuitFcnGrad.cxx b/roofit/roofitcore/src/TestStatistics/MinuitFcnGrad.cxx index 2ee0be304f4ad..7eac1268673c8 100644 --- a/roofit/roofitcore/src/TestStatistics/MinuitFcnGrad.cxx +++ b/roofit/roofitcore/src/TestStatistics/MinuitFcnGrad.cxx @@ -54,6 +54,9 @@ class MinuitGradFunctor : public ROOT::Minuit2::FCNBase { : ROOT::Minuit2::GradientParameterSpace::External; } + // TODO: Implement this + bool SecondDerivativeAlwaysVanishes(unsigned int /*i*/, unsigned int /*j*/) const override { return false; } + private: MinuitFcnGrad const &_fcn; double _up; diff --git a/roofit/roofitcore/test/testRooMinimizer.cxx b/roofit/roofitcore/test/testRooMinimizer.cxx index f565816adb474..05004ad171434 100644 --- a/roofit/roofitcore/test/testRooMinimizer.cxx +++ b/roofit/roofitcore/test/testRooMinimizer.cxx @@ -1,15 +1,19 @@ // Tests for the RooProdPdf // Author: Jonas Rembser, CERN, October 2024 +#include #include #include #include +#include #include #include #include #include #include +#include + #include "gtest_wrapper.h" class EvalBackendParametrizedTest : public testing::TestWithParam> { @@ -170,3 +174,76 @@ INSTANTIATE_TEST_SUITE_P(RooMinimizer, EvalBackendParametrizedTest, testing::Val ss << "EvalBackend" << std::get<0>(paramInfo.param).name(); return ss.str(); }); + +// Check the vanishing-second-derivative optimization in MnHesse, which is +// driven by the parameter independence information that RooMinimizerFcn +// derives from the computation graph. The covariance matrix from hesse() must +// agree with the analytical one, both when all parameters float and when a +// parameter is fixed after the RooMinimizer was constructed. The latter case +// is a regression test for the translation between Minuit-internal parameter +// indices (which exclude fixed parameters) and the external indices that +// RooMinimizerFcn::secondDerivativeAlwaysVanishes() is defined in. +TEST(RooMinimizer, SecondDerivativeAlwaysVanishesHesse) +{ + RooHelpers::LocalChangeMsgLevel chmsglvl{RooFit::WARNING}; + + RooRealVar a("a", "a", 1, -10, 10); + RooRealVar b("b", "b", 1, -10, 10); + RooRealVar c("c", "c", 1, -10, 10); + + // Two additive terms: t1 couples (a,c) and t2 couples (b,c), so only the + // mixed second derivative w.r.t. (a,b) vanishes identically and can be + // skipped in the numerical Hessian. + RooFormulaVar t1("t1", "t1", "(a-1)^2 + (a-c)^2 + (c-1)^2", {a, c}); + RooFormulaVar t2("t2", "t2", "(b-c)^2 + (b-1)^2", {b, c}); + RooAddition f("f", "f", {t1, t2}); + + auto resetParams = [&]() { + for (RooRealVar *var : {&a, &b, &c}) { + var->setVal(0.5); + var->setError(0.0); + var->setConstant(false); + } + }; + + // The function is quadratic, so the finite-difference Hessian is exact up + // to numerical noise. + const double tol = 1e-3; + + // All parameters floating: with errorDef = 1, the covariance is 2 * H^-1 + // with H = [[4, 0, -2], [0, 4, -2], [-2, -2, 6]] in the order (a, b, c). + { + resetParams(); + RooMinimizer m(f); + m.setPrintLevel(-1); + m.migrad(); + m.hesse(); + EXPECT_NEAR(a.getError(), std::sqrt(0.625), tol); + EXPECT_NEAR(b.getError(), std::sqrt(0.625), tol); + EXPECT_NEAR(c.getError(), std::sqrt(0.5), tol); + } + + // With "a" fixed, the reduced Hessian in (b, c) is [[4, -2], [-2, 6]], + // so the covariance is 2 * H^-1 = [[0.6, 0.2], [0.2, 0.4]]. This must not + // depend on whether "a" was fixed before or after constructing the + // minimizer: in the latter case, Minuit's internal parameter indices no + // longer coincide with the external ones that the vanishing-second- + // derivative mask is defined in, so wrong bookkeeping would zero out the + // genuinely non-vanishing (b, c) element here. + for (bool fixAfterConstruction : {false, true}) { + resetParams(); + a.setVal(1.0); + if (!fixAfterConstruction) { + a.setConstant(true); + } + RooMinimizer m(f); + m.setPrintLevel(-1); + if (fixAfterConstruction) { + a.setConstant(true); + } + m.migrad(); + m.hesse(); + EXPECT_NEAR(b.getError(), std::sqrt(0.6), tol) << "fixAfterConstruction = " << fixAfterConstruction; + EXPECT_NEAR(c.getError(), std::sqrt(0.4), tol) << "fixAfterConstruction = " << fixAfterConstruction; + } +}