From e1137acadc8e03d45775ae16e2bc79c9fb652635 Mon Sep 17 00:00:00 2001 From: Jonas Rembser Date: Wed, 9 Sep 2026 08:15:27 +0000 Subject: [PATCH 1/4] [RF] Vectorize the batched reduction in RooUnbinnedL MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The evaluator-backed RooUnbinnedL computed the negative log-likelihood from the batch of probabilities with a scalar std::log loop, costing several times more than one vectorized likelihood pass now that the rest of the evaluation is SIMD. Use the same RooBatchCompute::reduceNLL() reduction as RooNLLVarNew, which also reproduces the RooNaNPacker-based error propagation of the loop (badness packed into the returned NaN). 🤖 Done with the help of AI --- .../inc/RooFit/TestStatistics/RooUnbinnedL.h | 1 + .../src/TestStatistics/RooUnbinnedL.cxx | 57 +++++++++---------- 2 files changed, 28 insertions(+), 30 deletions(-) diff --git a/roofit/roofitcore/inc/RooFit/TestStatistics/RooUnbinnedL.h b/roofit/roofitcore/inc/RooFit/TestStatistics/RooUnbinnedL.h index 4888f883aecd0..df5e7c2d9ceb1 100644 --- a/roofit/roofitcore/inc/RooFit/TestStatistics/RooUnbinnedL.h +++ b/roofit/roofitcore/inc/RooFit/TestStatistics/RooUnbinnedL.h @@ -52,6 +52,7 @@ class RooUnbinnedL : public RooAbsL { mutable ROOT::Math::KahanSum cachedResult_{0.}; std::shared_ptr evaluator_; ///> _vectorBuffers; // used for preserving resources in batched evaluation + std::vector _unitWeights; /// #include #include +#include #include #include #include @@ -148,42 +149,38 @@ ComputeResult computeScalarFunc(const RooAbsPdf *pdfClone, RooAbsData *dataClone // Similar to computeScalarFunc, but the probabilities were already evaluated // as a batch, and the weights are also retrieved as batches instead of looping // over RooAbsData::get(i), which loads every column of the dataset only to -// then read a single weight. +// then read a single weight. The reduction is done with the same vectorized +// RooBatchCompute::reduceNLL() that RooNLLVarNew uses in the standard +// evaluation backend, including its RooNaNPacker-based error propagation. ComputeResult computeBatchFunc(std::span probas, RooAbsData *dataClone, bool weightSq, - std::size_t stepSize, std::size_t firstEvent, std::size_t lastEvent) + std::size_t firstEvent, std::size_t lastEvent, std::vector &unitWeights, + RooBatchCompute::Config const &cfg) { - ROOT::Math::KahanSum kahanWeight; - ROOT::Math::KahanSum kahanProb; - RooNaNPacker packedNaN(0.f); - const std::size_t nEvents = lastEvent - firstEvent; // Empty spans mean the dataset is unweighted, i.e. all weights are one. - std::span weights = dataClone->getWeightBatch(firstEvent, nEvents, /*sumW2=*/false); - std::span weightsSumW2 = - weightSq ? dataClone->getWeightBatch(firstEvent, nEvents, /*sumW2=*/true) : std::span{}; - - for (auto i = firstEvent; i < lastEvent; i += stepSize) { - double weight = weights.empty() ? 1.0 : weights[i - firstEvent]; - - if (0. == weight * weight) - continue; - if (weightSq) - weight = weightsSumW2.empty() ? 1.0 : weightsSumW2[i - firstEvent]; - - double logProba = std::log(probas[i]); - const double term = -weight * logProba; + std::span dataWeights = dataClone->getWeightBatch(firstEvent, nEvents, /*sumW2=*/weightSq); - kahanWeight.Add(weight); - kahanProb.Add(term); - packedNaN.accumulate(term); + double sumWeight; + const double *weightData = nullptr; + std::size_t nWeights = 0; + if (dataWeights.empty()) { + if (unitWeights.size() < nEvents) { + unitWeights.assign(nEvents, 1.0); + } + weightData = unitWeights.data(); + nWeights = nEvents; + sumWeight = nEvents; + } else { + weightData = dataWeights.data(); + nWeights = dataWeights.size(); + sumWeight = RooBatchCompute::reduceSum(cfg, weightData, nWeights); } + std::span weights{weightData, nWeights}; - if (packedNaN.getPayload() != 0.) { - // Some events with evaluation errors. Return "badness" of errors. - return {ROOT::Math::KahanSum{packedNaN.getNaNWithPayload()}, kahanWeight.Sum()}; - } + std::span probasInRange{probas.data() + firstEvent, nEvents}; - return {kahanProb, kahanWeight.Sum()}; + auto out = RooBatchCompute::reduceNLL(cfg, probasInRange, weights, {}); + return {ROOT::Math::KahanSum{out.nllSum, out.nllSumCarry}, sumWeight}; } } // namespace @@ -214,8 +211,8 @@ RooUnbinnedL::evaluatePartition(Section events, std::size_t /*components_begin*/ // Here, we have a memory allocation that should be avoided when this // code needs to be optimized. std::span probas = evaluator_->run(); - std::tie(result, sumWeight) = - computeBatchFunc(probas, data_.get(), apply_weight_squared, 1, events.begin(N_events_), events.end(N_events_)); + std::tie(result, sumWeight) = computeBatchFunc(probas, data_.get(), apply_weight_squared, events.begin(N_events_), + events.end(N_events_), _unitWeights, RooBatchCompute::Config{}); } else { std::tie(result, sumWeight) = computeScalarFunc(pdf_.get(), data_.get(), normSet_.get(), apply_weight_squared, 1, events.begin(N_events_), events.end(N_events_)); From 4e370fa6b4a996269f4cc98ee780aa7520b1a77a Mon Sep 17 00:00:00 2001 From: Jonas Rembser Date: Tue, 8 Sep 2026 23:13:35 +0000 Subject: [PATCH 2/4] [RF] Forward the EvalBackend option to modular likelihoods in createNLL MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The ModularL branch of createNLL never passed the parsed EvalBackend on to the NLLFactory, whose default is the legacy backend. As a result, every parallel fit (fitTo with Parallelize(), which implies ModularL) silently evaluated the deprecated legacy scalar likelihood on the workers, costing about 7x per likelihood evaluation compared to the vectorized CPU backend and making the parallel gradient lose against any serial fit (benchmark: 8-channel unbinned simultaneous fit with 64 correlated constrained systematics, 96 free parameters, 200k events: serial fitTo 9.9s, Parallelize(8) 23.7s before, 4.7s after this and the accompanying gradient-job commits). The bitwise legacy-vs-modular comparisons in testLikelihoodGradientJob now request EvalBackend(Legacy) explicitly on the modular side: they compare against a legacy reference fit, so both likelihoods must use the same arithmetic (previously that happened by virtue of this bug). 🤖 Done with the help of AI --- roofit/roofitcore/src/FitHelpers.cxx | 3 ++- .../src/TestStatistics/buildLikelihood.cxx | 15 ++++++++++++++- .../TestStatistics/testLikelihoodGradientJob.cxx | 12 +++++++++--- 3 files changed, 25 insertions(+), 5 deletions(-) diff --git a/roofit/roofitcore/src/FitHelpers.cxx b/roofit/roofitcore/src/FitHelpers.cxx index fc5759c8c2553..99f87a3f33592 100644 --- a/roofit/roofitcore/src/FitHelpers.cxx +++ b/roofit/roofitcore/src/FitHelpers.cxx @@ -825,7 +825,8 @@ std::unique_ptr createNLL(RooAbsPdf &pdf, RooAbsData &data, const Ro .ConstrainedParameters(cParsSet) .ExternalConstraints(extConsSet) .GlobalObservables(glObsSet) - .GlobalObservablesTag(rangeName.c_str()); + .GlobalObservablesTag(rangeName.c_str()) + .EvalBackend(RooFit::EvalBackend(static_cast(pc.getInt("EvalBackend")))); return std::make_unique("likelihood", "", builder.build()); } diff --git a/roofit/roofitcore/src/TestStatistics/buildLikelihood.cxx b/roofit/roofitcore/src/TestStatistics/buildLikelihood.cxx index bac8a76ab0bb3..8db7e3dc49d3e 100644 --- a/roofit/roofitcore/src/TestStatistics/buildLikelihood.cxx +++ b/roofit/roofitcore/src/TestStatistics/buildLikelihood.cxx @@ -294,7 +294,20 @@ std::vector> NLLFactory::getSimultaneousComponents() RooArgSet selTargetParams; params.selectCommon(*actualParams, selTargetParams); - assert(selTargetParams.equals(*components.back()->getParameters())); + // Sanity check that the parameters of the component likelihood are among the parameters that the component + // pdf has in common with the full model. With a non-legacy evaluation backend, the component pdf is + // additionally "compiled" for the dataset, which can move constant constraint factors out of the pdf (see + // RooProdPdf::compileForNormSet and RooFixedProdPdf). The component likelihood can then legitimately report + // fewer parameters than the original component pdf, so only containment can be asserted, not set equality. + assert([&] { + std::unique_ptr componentParams{components.back()->getParameters()}; + for (auto *param : *componentParams) { + if (!selTargetParams.find(param->GetName())) { + return false; + } + } + return true; + }()); ++n; } else { diff --git a/roofit/roofitcore/test/TestStatistics/testLikelihoodGradientJob.cxx b/roofit/roofitcore/test/TestStatistics/testLikelihoodGradientJob.cxx index a2e43b78e7e65..db0641eba7977 100644 --- a/roofit/roofitcore/test/TestStatistics/testLikelihoodGradientJob.cxx +++ b/roofit/roofitcore/test/TestStatistics/testLikelihoodGradientJob.cxx @@ -651,7 +651,11 @@ TEST_P(LikelihoodGradientJobErrorTest, ErrorHandling) values.assign(savedValues); - std::unique_ptr likelihoodAbsReal{pdf->createNLL(*data, RooFit::ModularL(true))}; + // Explicitly request the legacy backend also here: this test compares + // bitwise against the legacy nominal fit above, so both likelihoods must + // use the same arithmetic. + std::unique_ptr likelihoodAbsReal{ + pdf->createNLL(*data, RooFit::ModularL(true), RooFit::EvalBackend(RooFit::EvalBackend::Value::Legacy))}; RooMinimizer::Config cfg; cfg.parallelize = NWorkers; @@ -709,8 +713,10 @@ TEST_P(LikelihoodGradientJobErrorTest, FitSimpleLinear) std::unique_ptr fitResult{minim.save()}; auto a1Result = a1.getVal(); - // now with multiprocess - std::unique_ptr nll_mp(pdf.createNLL(*data, RooFit::ModularL(true))); + // now with multiprocess; explicitly request the legacy backend to compare + // bitwise against the legacy nominal fit above + std::unique_ptr nll_mp( + pdf.createNLL(*data, RooFit::ModularL(true), RooFit::EvalBackend(RooFit::EvalBackend::Value::Legacy))); a1.setVal(-5.); a1.removeError(); From 97acd18bb768359de22c0046b5a5738d44f63809 Mon Sep 17 00:00:00 2001 From: Jonas Rembser Date: Tue, 8 Sep 2026 23:14:32 +0000 Subject: [PATCH 3/4] [RF] Skip the central-point evaluation in parallel gradient workers MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Every gradient calculation in LikelihoodGradientJob started with each worker evaluating the full likelihood once at the central point, only to obtain the scalar function value that NumericalDerivator needs for its step-size tolerances (SetupDifferentiate). Minuit already knows this exact value: the line search that precedes each gradient request stores it in MinimumParameters::Fval(). Hand that value through a new fifth argument of FCNBase::GradientWithPrevResult() (a backward-compatible overload that falls back to the old virtual), broadcast it to the workers along with the rest of the minimizer state, and pre-seed the derivator's existing central-value cache (fVxFValCache) with it, so SetupDifferentiate skips its function call. When the likelihood offsets changed in the same state update, NaN is broadcast instead and the workers evaluate as before, since the known value corresponds to the previous offsets. This is bitwise-transparent: the exact-equality comparisons in testLikelihoodGradientJob pass unchanged with the pre-seeding active. 🤖 Done with the help of AI --- math/minuit2/inc/Minuit2/FCNBase.h | 14 ++++++++++ math/minuit2/inc/Minuit2/NumericalDerivator.h | 12 +++++++++ .../ExternalInternalGradientCalculator.cxx | 3 ++- .../LikelihoodGradientWrapper.h | 7 +++-- .../TestStatistics/LikelihoodGradientJob.cxx | 26 ++++++++++++++++--- .../TestStatistics/LikelihoodGradientJob.h | 8 ++++-- .../src/TestStatistics/MinuitFcnGrad.cxx | 8 +++--- .../src/TestStatistics/MinuitFcnGrad.h | 2 +- 8 files changed, 67 insertions(+), 13 deletions(-) diff --git a/math/minuit2/inc/Minuit2/FCNBase.h b/math/minuit2/inc/Minuit2/FCNBase.h index f9798404d1152..c18c3785e0783 100644 --- a/math/minuit2/inc/Minuit2/FCNBase.h +++ b/math/minuit2/inc/Minuit2/FCNBase.h @@ -102,6 +102,20 @@ class FCNBase { return Gradient(parameters); }; + /// Variant of GradientWithPrevResult() that additionally receives the + /// already-known function value at \p parameters (e.g. from the line search + /// that preceded the gradient request), so implementations can avoid + /// re-evaluating the function at the central point. + /// + /// \warning Not meant to be overridden! This is a requirement for an + /// internal optimization in RooFit that might go away with any refactoring. + virtual std::vector GradientWithPrevResult(std::vector const ¶meters, double *previous_grad, + double *previous_g2, double *previous_gstep, + double /*fValAtParameters*/) const + { + return GradientWithPrevResult(parameters, previous_grad, previous_g2, previous_gstep); + }; + /// \warning Not meant to be overridden! This is a requirement for an /// internal optimization in RooFit that might go away with any refactoring. virtual GradientParameterSpace gradParameterSpace() const { return GradientParameterSpace::External; }; diff --git a/math/minuit2/inc/Minuit2/NumericalDerivator.h b/math/minuit2/inc/Minuit2/NumericalDerivator.h index c4ad937fe594f..f8e42fc35fd7a 100644 --- a/math/minuit2/inc/Minuit2/NumericalDerivator.h +++ b/math/minuit2/inc/Minuit2/NumericalDerivator.h @@ -46,6 +46,18 @@ class NumericalDerivator { void SetupDifferentiate(unsigned int nDim, const FCNBase *function, const double *cx, std::span parameters); + + /// Pre-seed the cache of the function value at the central point, so that a + /// subsequent SetupDifferentiate() at the same point \p cx (in Minuit-internal + /// coordinates) can skip its function evaluation. Callers that already know + /// the function value at the gradient point (e.g. from the line search that + /// Minuit just completed) can use this to avoid one full function call. + void PreseedFVal(double fval, std::span cx) + { + fVxFValCache.assign(cx.begin(), cx.end()); + fVal = fval; + } + std::vector Differentiate(unsigned int nDim, const FCNBase *function, const double *x, std::span parameters, std::span previous_gradient); diff --git a/math/minuit2/src/ExternalInternalGradientCalculator.cxx b/math/minuit2/src/ExternalInternalGradientCalculator.cxx index c5de47eab1bff..5e11abf2d9025 100644 --- a/math/minuit2/src/ExternalInternalGradientCalculator.cxx +++ b/math/minuit2/src/ExternalInternalGradientCalculator.cxx @@ -54,7 +54,8 @@ ExternalInternalGradientCalculator::operator()(const MinimumParameters &par, con std::vector previous_g2(functionGradient.G2().Data(), functionGradient.G2().Data() + functionGradient.G2().size()); std::vector previous_gstep(functionGradient.Gstep().Data(), functionGradient.Gstep().Data() + functionGradient.Gstep().size()); - std::vector grad = fGradFunc.GradientWithPrevResult(par_vec, previous_grad.data(), previous_g2.data(), previous_gstep.data()); + std::vector grad = fGradFunc.GradientWithPrevResult(par_vec, previous_grad.data(), previous_g2.data(), + previous_gstep.data(), par.Fval()); assert(grad.size() == fTransformation.Parameters().size()); MnAlgebraicVector v(par.Vec().size()); diff --git a/roofit/roofitcore/inc/RooFit/TestStatistics/LikelihoodGradientWrapper.h b/roofit/roofitcore/inc/RooFit/TestStatistics/LikelihoodGradientWrapper.h index 825b01576c10d..7a4a480b037e8 100644 --- a/roofit/roofitcore/inc/RooFit/TestStatistics/LikelihoodGradientWrapper.h +++ b/roofit/roofitcore/inc/RooFit/TestStatistics/LikelihoodGradientWrapper.h @@ -52,8 +52,11 @@ class LikelihoodGradientWrapper { SharedOffset offset); virtual void fillGradient(double *grad) = 0; - virtual void - fillGradientWithPrevResult(double *grad, double *previous_grad, double *previous_g2, double *previous_gstep) = 0; + /// \param[in] fValAtX The function value at the current parameter point, if known by the + /// caller (e.g. from the preceding line search); NaN when unknown. Implementations + /// can use it to avoid re-evaluating the function at the central point. + virtual void fillGradientWithPrevResult(double *grad, double *previous_grad, double *previous_g2, + double *previous_gstep, double fValAtX) = 0; /// Synchronize minimizer settings with calculators in child classes. virtual void synchronizeWithMinimizer(const ROOT::Math::MinimizerOptions &options); diff --git a/roofit/roofitcore/src/TestStatistics/LikelihoodGradientJob.cxx b/roofit/roofitcore/src/TestStatistics/LikelihoodGradientJob.cxx index 58ae110730382..223e75824fb83 100644 --- a/roofit/roofitcore/src/TestStatistics/LikelihoodGradientJob.cxx +++ b/roofit/roofitcore/src/TestStatistics/LikelihoodGradientJob.cxx @@ -24,6 +24,8 @@ #include "Minuit2/Minuit2Minimizer.h" #include "Minuit2/MnStrategy.h" +#include + namespace RooFit { namespace TestStatistics { @@ -122,14 +124,18 @@ void LikelihoodGradientJob::update_workers_state() ++state_id_; if (shared_offset_.offsets() != offsets_previous_) { + // The function value known by the master corresponds to the previous + // offsets, so it cannot be used to skip the central-point evaluation on + // the workers in this update. + double fValAtX = std::numeric_limits::quiet_NaN(); zmq::message_t offsets_message(shared_offset_.offsets().begin(), shared_offset_.offsets().end()); get_manager()->messenger().publish_from_master_to_workers( - id_, state_id_, isCalculating_, maxFCN, fcnOffset, std::move(gradient_message), + id_, state_id_, isCalculating_, maxFCN, fcnOffset, fValAtX, std::move(gradient_message), std::move(minuit_internal_x_message), std::move(offsets_message)); offsets_previous_ = shared_offset_.offsets(); } else { get_manager()->messenger().publish_from_master_to_workers(id_, state_id_, isCalculating_, maxFCN, fcnOffset, - std::move(gradient_message), + fval_at_x_, std::move(gradient_message), std::move(minuit_internal_x_message)); } } @@ -157,6 +163,9 @@ void LikelihoodGradientJob::update_state() minimizer_->fcnOffset() = fcnOffset; assert(more); + auto fValAtX = get_manager()->messenger().receive_from_master_on_worker(&more); + assert(more); + auto gradient_message = get_manager()->messenger().receive_from_master_on_worker(&more); assert(more); auto gradient_message_begin = gradient_message.data(); @@ -184,6 +193,14 @@ void LikelihoodGradientJob::update_state() // Since the gradient parallelization only support Minuit 2, we can do this cast auto &minim = static_cast(*minimizer_->_minimizer); + // The master already knows the function value at the current point from + // the line search that Minuit just completed; pre-seeding the derivator's + // cache with it makes the SetupDifferentiate call below skip its full + // likelihood evaluation at the central point. + if (!std::isnan(fValAtX)) { + gradf_.PreseedFVal(fValAtX, {minuit_internal_x_.data(), static_cast(minimizer_->getNPar())}); + } + // note: the next call must stay after the (possible) update of the offset, because it // calls the likelihood function, so the offset must be correct at this point gradf_.SetupDifferentiate(minimizer_->getNPar(), minim.GetFCN(), minuit_internal_x_.data(), @@ -231,6 +248,7 @@ void LikelihoodGradientJob::fillGradient(double *grad) { if (get_manager()->process_manager().is_master()) { if (!calculation_is_clean_->gradient) { + fval_at_x_ = std::numeric_limits::quiet_NaN(); calculate_all(); } @@ -242,13 +260,15 @@ void LikelihoodGradientJob::fillGradient(double *grad) } void LikelihoodGradientJob::fillGradientWithPrevResult(double *grad, double *previous_grad, double *previous_g2, - double *previous_gstep) + double *previous_gstep, double fValAtX) { if (get_manager()->process_manager().is_master()) { for (std::size_t i_component = 0; i_component < N_tasks_; ++i_component) { grad_[i_component] = {previous_grad[i_component], previous_g2[i_component], previous_gstep[i_component]}; } + fval_at_x_ = fValAtX; + if (!calculation_is_clean_->gradient) { if (RooFit::MultiProcess::Config::getTimingAnalysis()) { RooFit::MultiProcess::ProcessTimer::start_timer("master:gradient"); diff --git a/roofit/roofitcore/src/TestStatistics/LikelihoodGradientJob.h b/roofit/roofitcore/src/TestStatistics/LikelihoodGradientJob.h index b54f61606cdeb..ea9164e28685d 100644 --- a/roofit/roofitcore/src/TestStatistics/LikelihoodGradientJob.h +++ b/roofit/roofitcore/src/TestStatistics/LikelihoodGradientJob.h @@ -20,6 +20,7 @@ #include "Minuit2/NumericalDerivator.h" #include "Minuit2/MnMatrix.h" +#include #include namespace RooFit { @@ -32,8 +33,8 @@ class LikelihoodGradientJob : public MultiProcess::Job, public LikelihoodGradien RooMinimizer *minimizer, SharedOffset offset); void fillGradient(double *grad) override; - void fillGradientWithPrevResult(double *grad, double *previous_grad, double *previous_g2, - double *previous_gstep) override; + void fillGradientWithPrevResult(double *grad, double *previous_grad, double *previous_g2, double *previous_gstep, + double fValAtX) override; void update_state() override; @@ -80,6 +81,9 @@ class LikelihoodGradientJob : public MultiProcess::Job, public LikelihoodGradien std::size_t N_tasks_ = 0; std::size_t N_tasks_at_workers_ = 0; std::vector minuit_internal_x_; + /// Function value at minuit_internal_x_ as known by the master (NaN when unknown); broadcast + /// to workers so their NumericalDerivator setup can skip the central-point evaluation. + double fval_at_x_ = std::numeric_limits::quiet_NaN(); mutable bool isCalculating_ = false; diff --git a/roofit/roofitcore/src/TestStatistics/MinuitFcnGrad.cxx b/roofit/roofitcore/src/TestStatistics/MinuitFcnGrad.cxx index 7eac1268673c8..4fcb3b76a280b 100644 --- a/roofit/roofitcore/src/TestStatistics/MinuitFcnGrad.cxx +++ b/roofit/roofitcore/src/TestStatistics/MinuitFcnGrad.cxx @@ -42,10 +42,10 @@ class MinuitGradFunctor : public ROOT::Minuit2::FCNBase { return grad; } std::vector GradientWithPrevResult(std::vector const &v, double *previous_grad, double *previous_g2, - double *previous_gstep) const override + double *previous_gstep, double fValAtV) const override { std::vector output(v.size()); - _fcn.GradientWithPrevResult(v.data(), output.data(), previous_grad, previous_g2, previous_gstep); + _fcn.GradientWithPrevResult(v.data(), output.data(), previous_grad, previous_g2, previous_gstep, fValAtV); return output; } ROOT::Minuit2::GradientParameterSpace gradParameterSpace() const override @@ -259,12 +259,12 @@ void MinuitFcnGrad::Gradient(const double *x, double *grad) const } void MinuitFcnGrad::GradientWithPrevResult(const double *x, double *grad, double *previous_grad, double *previous_g2, - double *previous_gstep) const + double *previous_gstep, double fValAtX) const { _calculatingGradient = true; syncParameterValuesFromMinuitCalls(x, returnsInMinuit2ParameterSpace()); syncOffsets(); - _gradient->fillGradientWithPrevResult(grad, previous_grad, previous_g2, previous_gstep); + _gradient->fillGradientWithPrevResult(grad, previous_grad, previous_g2, previous_gstep, fValAtX); _calculatingGradient = false; } diff --git a/roofit/roofitcore/src/TestStatistics/MinuitFcnGrad.h b/roofit/roofitcore/src/TestStatistics/MinuitFcnGrad.h index 902381c9c98aa..f345fdf67e188 100644 --- a/roofit/roofitcore/src/TestStatistics/MinuitFcnGrad.h +++ b/roofit/roofitcore/src/TestStatistics/MinuitFcnGrad.h @@ -46,7 +46,7 @@ class MinuitFcnGrad : public RooAbsMinimizerFcn { /// IMultiGradFunction overrides necessary for Minuit void Gradient(const double *x, double *grad) const; void GradientWithPrevResult(const double *x, double *grad, double *previous_grad, double *previous_g2, - double *previous_gstep) const; + double *previous_gstep, double fValAtX) const; inline std::string getFunctionName() const override { return _likelihood->GetName(); } From 1e28d69bf719c6f3bc6160d3000c2ab667ed0ca4 Mon Sep 17 00:00:00 2001 From: Jonas Rembser Date: Tue, 8 Sep 2026 23:14:52 +0000 Subject: [PATCH 4/4] [RF] Chunk the parallel gradient into multi-parameter tasks MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit One gradient task used to be one partial derivative, so a gradient call in LikelihoodGradientJob costed one dequeue round trip through the queue process plus one result message per parameter. With the vectorized CPU backend a partial derivative is only a few milliseconds of work, and for fits with around a hundred parameters this scheduling overhead dominated the parallel gradient and capped its speedup. Each task now covers a strided block of parameters (indices congruent to the task id), so parameters with expensive derivatives - typically adjacent blocks of correlated systematics - spread evenly over the tasks, and the per-task results travel in one message. The number of tasks is configurable through the new Config::LikelihoodGradientJob::defaultNParamTasks; the automatic default of four tasks per worker keeps enough granularity for the queue to balance load between workers. Benchmark (96-parameter unbinned simultaneous fit, 8 channels x 25k events, medians of 3 interleaved repeats on 12 cores): Parallelize(8) 9.5s with per-parameter tasks, 4.7s with the automatic chunking, against 9.8s for the serial fit. 🤖 Done with the help of AI --- .../inc/RooFit/MultiProcess/Config.h | 12 +++++ roofit/multiprocess/src/Config.cxx | 1 + .../TestStatistics/LikelihoodGradientJob.cxx | 45 ++++++++++++++++--- .../TestStatistics/LikelihoodGradientJob.h | 4 +- 4 files changed, 54 insertions(+), 8 deletions(-) diff --git a/roofit/multiprocess/inc/RooFit/MultiProcess/Config.h b/roofit/multiprocess/inc/RooFit/MultiProcess/Config.h index 42508cc2207ad..bf1ae18a44f83 100644 --- a/roofit/multiprocess/inc/RooFit/MultiProcess/Config.h +++ b/roofit/multiprocess/inc/RooFit/MultiProcess/Config.h @@ -38,6 +38,18 @@ class Config { static std::size_t defaultNComponentTasks; }; + struct LikelihoodGradientJob { + // magic value to indicate that the number of tasks will be set automatically + constexpr static std::size_t automaticNParamTasks = 0; + + /// Number of tasks to split a gradient calculation into. Each task + /// covers a contiguous block of parameters, so fewer tasks mean less + /// scheduling overhead per gradient, but also less opportunity for load + /// balancing between the workers. The automatic default uses a small + /// multiple of the number of workers. + static std::size_t defaultNParamTasks; + }; + struct Queue { enum class QueueType {FIFO, Priority}; static bool setQueueType(QueueType queueType); diff --git a/roofit/multiprocess/src/Config.cxx b/roofit/multiprocess/src/Config.cxx index 1dfe3adb40687..23d6bce1697d9 100644 --- a/roofit/multiprocess/src/Config.cxx +++ b/roofit/multiprocess/src/Config.cxx @@ -146,6 +146,7 @@ void Config::Queue::suggestTaskOrder(std::size_t job_id, const std::vector unsigned int Config::defaultNWorkers_ = std::thread::hardware_concurrency(); std::size_t Config::LikelihoodJob::defaultNEventTasks = Config::LikelihoodJob::automaticNEventTasks; std::size_t Config::LikelihoodJob::defaultNComponentTasks = Config::LikelihoodJob::automaticNComponentTasks; +std::size_t Config::LikelihoodGradientJob::defaultNParamTasks = Config::LikelihoodGradientJob::automaticNParamTasks; Config::Queue::QueueType Config::Queue::queueType_ = Config::Queue::QueueType::FIFO; bool Config::timingAnalysis_ = false; diff --git a/roofit/roofitcore/src/TestStatistics/LikelihoodGradientJob.cxx b/roofit/roofitcore/src/TestStatistics/LikelihoodGradientJob.cxx index 223e75824fb83..cffdee74e2e2a 100644 --- a/roofit/roofitcore/src/TestStatistics/LikelihoodGradientJob.cxx +++ b/roofit/roofitcore/src/TestStatistics/LikelihoodGradientJob.cxx @@ -24,6 +24,7 @@ #include "Minuit2/Minuit2Minimizer.h" #include "Minuit2/MnStrategy.h" +#include #include namespace RooFit { @@ -34,13 +35,33 @@ LikelihoodGradientJob::LikelihoodGradientJob(std::shared_ptr likelihood std::size_t N_dim, RooMinimizer *minimizer, SharedOffset offset) : LikelihoodGradientWrapper(std::move(likelihood), std::move(calculation_is_clean), N_dim, minimizer, std::move(offset)), - grad_(N_dim), - N_tasks_(N_dim) + grad_(N_dim) { + // Each task covers multiple parameters to reduce the scheduling overhead + // per gradient: every task costs a dequeue round trip through the queue + // process plus a result message to the master, which for cheap partial + // derivatives can otherwise dominate over the actual calculation. A few + // tasks per worker still leave the queue room for load balancing. Parameters + // are assigned to tasks in strides, so parameters with expensive partial + // derivatives (that often sit next to each other in the parameter order, + // e.g. a block of correlated systematics) spread evenly over the tasks. + N_tasks_ = MultiProcess::Config::LikelihoodGradientJob::defaultNParamTasks; + if (N_tasks_ == MultiProcess::Config::LikelihoodGradientJob::automaticNParamTasks) { + N_tasks_ = 4 * MultiProcess::Config::getDefaultNWorkers(); + } + N_tasks_ = std::min(N_tasks_, N_dim); + minuit_internal_x_.reserve(N_dim); offsets_previous_ = shared_offset_.offsets(); } +/// Number of parameters assigned to task \p task (parameter indices congruent +/// to \p task modulo N_tasks_). +std::size_t LikelihoodGradientJob::taskSize(std::size_t task) const +{ + return grad_.size() / N_tasks_ + (task < grad_.size() % N_tasks_ ? 1 : 0); +} + void LikelihoodGradientJob::synchronizeParameterSettingsImpl( const std::vector ¶meter_settings) { @@ -88,23 +109,33 @@ void LikelihoodGradientJob::setErrorLevel(double error_level) const void LikelihoodGradientJob::evaluate_task(std::size_t task) { - run_derivator(task); + for (std::size_t ix = task; ix < grad_.size(); ix += N_tasks_) { + run_derivator(ix); + } } // SYNCHRONIZATION FROM WORKERS TO MASTER void LikelihoodGradientJob::send_back_task_result_from_worker(std::size_t task) { - task_result_t task_result{id_, task, grad_[task]}; - zmq::message_t message(sizeof(task_result_t)); + task_result_t task_result{id_, task}; + zmq::message_t message(sizeof(task_result_t) + taskSize(task) * sizeof(ROOT::Minuit2::DerivatorElement)); memcpy(message.data(), &task_result, sizeof(task_result_t)); + auto elements = reinterpret_cast(message.data() + sizeof(task_result_t)); + for (std::size_t ix = task; ix < grad_.size(); ix += N_tasks_) { + *elements++ = grad_[ix]; + } get_manager()->messenger().send_from_worker_to_master(std::move(message)); } bool LikelihoodGradientJob::receive_task_result_on_master(const zmq::message_t &message) { auto result = message.data(); - grad_[result->task_id] = result->grad; + auto elements = + reinterpret_cast(message.data() + sizeof(task_result_t)); + for (std::size_t ix = result->task_id; ix < grad_.size(); ix += N_tasks_) { + grad_[ix] = *elements++; + } --N_tasks_at_workers_; bool job_completed = (N_tasks_at_workers_ == 0); return job_completed; @@ -263,7 +294,7 @@ void LikelihoodGradientJob::fillGradientWithPrevResult(double *grad, double *pre double *previous_gstep, double fValAtX) { if (get_manager()->process_manager().is_master()) { - for (std::size_t i_component = 0; i_component < N_tasks_; ++i_component) { + for (std::size_t i_component = 0; i_component < grad_.size(); ++i_component) { grad_[i_component] = {previous_grad[i_component], previous_g2[i_component], previous_gstep[i_component]}; } diff --git a/roofit/roofitcore/src/TestStatistics/LikelihoodGradientJob.h b/roofit/roofitcore/src/TestStatistics/LikelihoodGradientJob.h index ea9164e28685d..8c5672726da62 100644 --- a/roofit/roofitcore/src/TestStatistics/LikelihoodGradientJob.h +++ b/roofit/roofitcore/src/TestStatistics/LikelihoodGradientJob.h @@ -60,11 +60,13 @@ class LikelihoodGradientJob : public MultiProcess::Job, public LikelihoodGradien // Job overrides: void evaluate_task(std::size_t task) override; + /// Message header for a task result; followed in the same message by the + /// DerivatorElement results for the parameters of that task. struct task_result_t { std::size_t job_id; std::size_t task_id; - ROOT::Minuit2::DerivatorElement grad; }; + std::size_t taskSize(std::size_t task) const; void send_back_task_result_from_worker(std::size_t task) override; bool receive_task_result_on_master(const zmq::message_t &message) override;