Skip to content

[RF] Test RooFit Hessians with Clad - #21622

Draft
guitargeek wants to merge 9 commits into
root-project:masterfrom
guitargeek:roofit_clad_hessians
Draft

guitargeek wants to merge 9 commits into
root-project:masterfrom
guitargeek:roofit_clad_hessians

Conversation

@guitargeek

Copy link
Copy Markdown
Contributor

No description provided.

@guitargeek guitargeek self-assigned this Mar 16, 2026
@github-actions

github-actions Bot commented Mar 17, 2026 •

Copy link
Copy Markdown

Test Results

    24 files      24 suites   4d 0h 38m 2s ⏱️
 3 886 tests  3 883 ✅ 0 💤 3 ❌
83 181 runs  83 178 ✅ 0 💤 3 ❌

For more details on these failures, see this check.

Results for commit 09938e1.

♻️ This comment has been updated with latest results.

@guitargeek
guitargeek force-pushed the roofit_clad_hessians branch from e54bf52 to 2a2b0cc Compare March 17, 2026 09:52
@guitargeek guitargeek added the clean build Ask CI to do non-incremental build on PR label Mar 29, 2026
@guitargeek
guitargeek force-pushed the roofit_clad_hessians branch 2 times, most recently from 6368d0c to 11eab4c Compare March 29, 2026 18:19
@guitargeek guitargeek removed the clean build Ask CI to do non-incremental build on PR label Apr 20, 2026
@guitargeek
guitargeek force-pushed the roofit_clad_hessians branch 2 times, most recently from 92c3589 to a79ccd5 Compare April 26, 2026 13:45
@guitargeek
guitargeek force-pushed the roofit_clad_hessians branch from a79ccd5 to db537ff Compare June 17, 2026 06:56
@guitargeek
guitargeek force-pushed the roofit_clad_hessians branch from db537ff to 35d56ad Compare July 5, 2026 18:38
@guitargeek guitargeek added the clean build Ask CI to do non-incremental build on PR label Jul 27, 2026
@guitargeek
guitargeek force-pushed the roofit_clad_hessians branch 7 times, most recently from c63fa60 to 5e627d1 Compare August 13, 2026 20:48
@guitargeek guitargeek closed this Aug 13, 2026
@guitargeek guitargeek reopened this Aug 13, 2026
@guitargeek
guitargeek force-pushed the roofit_clad_hessians branch from 5e627d1 to 56addac Compare August 14, 2026 16:01
@guitargeek guitargeek closed this Aug 14, 2026
@guitargeek guitargeek reopened this Aug 14, 2026
@guitargeek guitargeek closed this Aug 16, 2026
@guitargeek guitargeek reopened this Aug 16, 2026
@guitargeek
guitargeek force-pushed the roofit_clad_hessians branch 2 times, most recently from d4726b8 to 29e1d3e Compare August 17, 2026 16:41
@guitargeek guitargeek closed this Aug 17, 2026
@guitargeek
guitargeek force-pushed the roofit_clad_hessians branch 2 times, most recently from 7449980 to 1fcc494 Compare August 30, 2026 09:20
@guitargeek
guitargeek force-pushed the roofit_clad_hessians branch 3 times, most recently from 994beed to 9201d97 Compare September 7, 2026 08:18
vgvassilev pushed a commit to guitargeek/clad that referenced this pull request Sep 7, 2026
Two instantiations of one template ask for the same derivative name, so
the lookup for the second finds the derivative built for the first:

  template <typename T> double scale(double x, T n) { return x * x * n; }
  double sumOfScales(double x) { return scale(x, 2.) + scale(x, 3); }

clad::hessian builds scale_pushforward_pullback for T = double and injects
it into the enclosing namespace. The T = int call finds only that overload.
With one candidate Sema builds a plain DeclRefExpr instead of an overload
set, and noOverloadExists compared only the argument count, so the call was
built anyway:

  error: cannot initialize a parameter of type 'double *' with an rvalue
         of type 'int *'
  note: passing argument to parameter '_d_n' here

<cmath> reaches the same state, where a floating-point and an integral
exponent are two instantiations of pow_pushforward. That broke the Hessian
of every RooFit likelihood containing a RooBernstein, whose integral mixes
std::pow(x, j + 1.) with std::pow(-1., j - i) (root-project/root#21622).

Run the single candidate through overload resolution and report no overload
when the arguments are not viable, so the caller derives the one that fits.
guitargeek added a commit to guitargeek/clad that referenced this pull request Sep 7, 2026
Two instantiations of one template ask for the same derivative name, so
the lookup for the second finds the derivative built for the first:

  template <typename T> double scale(double x, T n) { return x * x * n; }
  double sumOfScales(double x) { return scale(x, 2.) + scale(x, 3); }

clad::hessian builds scale_pushforward_pullback for T = double and injects
it into the enclosing namespace. The T = int call finds only that overload.
With one candidate Sema builds a plain DeclRefExpr instead of an overload
set, and noOverloadExists compared only the argument count, so the call was
built anyway:

  error: cannot initialize a parameter of type 'double *' with an rvalue
         of type 'int *'
  note: passing argument to parameter '_d_n' here

<cmath> reaches the same state, where a floating-point and an integral
exponent are two instantiations of pow_pushforward. That broke the Hessian
of every RooFit likelihood containing a RooBernstein, whose integral mixes
std::pow(x, j + 1.) with std::pow(-1., j - i) (root-project/root#21622).

Run the single candidate through overload resolution and report no overload
when the arguments are not viable, so the caller derives the one that fits.
guitargeek added a commit to vgvassilev/clad that referenced this pull request Sep 7, 2026
Two instantiations of one template ask for the same derivative name, so
the lookup for the second finds the derivative built for the first:

  template <typename T> double scale(double x, T n) { return x * x * n; }
  double sumOfScales(double x) { return scale(x, 2.) + scale(x, 3); }

clad::hessian builds scale_pushforward_pullback for T = double and injects
it into the enclosing namespace. The T = int call finds only that overload.
With one candidate Sema builds a plain DeclRefExpr instead of an overload
set, and noOverloadExists compared only the argument count, so the call was
built anyway:

  error: cannot initialize a parameter of type 'double *' with an rvalue
         of type 'int *'
  note: passing argument to parameter '_d_n' here

<cmath> reaches the same state, where a floating-point and an integral
exponent are two instantiations of pow_pushforward. That broke the Hessian
of every RooFit likelihood containing a RooBernstein, whose integral mixes
std::pow(x, j + 1.) with std::pow(-1., j - i) (root-project/root#21622).

Run the single candidate through overload resolution and report no overload
when the arguments are not viable, so the caller derives the one that fits.
@guitargeek
guitargeek force-pushed the roofit_clad_hessians branch 2 times, most recently from 14475ea to 08428cb Compare September 10, 2026 10:38
@guitargeek
guitargeek force-pushed the roofit_clad_hessians branch 2 times, most recently from 0200129 to 0a201b6 Compare September 25, 2026 09:19
@guitargeek guitargeek removed the clean build Ask CI to do non-incremental build on PR label Sep 28, 2026
@guitargeek
guitargeek force-pushed the roofit_clad_hessians branch 2 times, most recently from 8c5dc9a to 436d012 Compare September 29, 2026 14:27
@guitargeek
guitargeek force-pushed the roofit_clad_hessians branch from 436d012 to 872b2bf Compare October 9, 2026 07:44
Add second-derivative support to the function interfaces, following
the existing HasGradient()/Gradient() pattern:

- IBaseFunctionMultiDimTempl gets HasHessian() (default false) and
  bool Hessian(x, hess), filling a row-major NDim() x NDim() array
  and returning false when not implemented. The full-matrix layout
  and bool return match the Minuit2 FCNBase convention. Note that
  FitMethodFunction already had a bool Hessian() override filling a
  packed lower triangle; that pre-existing layout now deviates from
  the base-class contract and is left untouched here.
- IBaseFunctionOneDim gets HasHessian(), SecondDerivative(), and a
  multi-dim-compatible Hessian(), implemented via a new private
  DoSecondDerivative() that throws by default.
- GradFunctor takes an optional Hessian std::function in its
  (f, dim, gradient) constructor.
- GradFunctor1D takes an optional second-derivative std::function in
  its two-function constructor. Its member-pointer constructor
  template needed an is_member_pointer constraint so that three
  plain function pointers select the new std::function overload
  instead.

This lets minimizers and, in particular, the RooFit codegen backend
query externally provided second derivatives, which is needed to
support Clad Hessians through opaque functor calls.

🤖 Done with the help of AI
The codegen backend already emitted a custom pullback for
RooFunctorBinding-style pdfs, forwarding to the wrapped functor's
Gradient(). That covers gradients, but Clad Hessians run in
reverse-over-forward mode: the forward pass needs a
`<name>_pushforward`, and the reverse pass over that opaque call
needs a `<name>_pushforward_pullback`. Without them, Hessian
generation failed for any model containing a bound functor.

When the wrapped function reports HasHessian(), additionally emit

- roo_functor_<addr>_pushforward: value plus grad . dx, from the
  functor's operator() and Gradient(), and
- roo_functor_<addr>_pushforward_pullback: the exact adjoint of the
  pushforward, using the functor's Hessian(),
      d_x[i]  += d_y.value * grad[i] + d_y.pushforward * (H dx)[i]
      d_dx[i] += d_y.pushforward * grad[i],

following the signature convention Clad uses for custom pushforward
pullbacks (see clad's test/Hessian/NestedArrays.C). The declared code
includes Clad's BuiltinDerivatives.h itself, since the functor
declaration is JIT-ed before RooFuncWrapper includes CladDerivator.h.

Functors that do not implement Hessian() keep the previous
gradient-only behavior.

🤖 Done with the help of AI
Enable the Hessian cross-check (Minuit's numeric Hesse on the
reference fit vs. clad::hessian on the generated code) for every test
in the suite. All Clad Hessians were validated against finite
differences of the exact Clad gradient using the writeDebugMacro()
output; where fit errors still differ, the numeric Hesse of the
reference fit is the imprecise side.

- Add a per-test hesseTolerance for the parameter-error comparison,
  since Minuit's numeric Hesse and the analytic Clad Hessian agree
  less tightly (1e-3 .. 5e-2 relative) than the fitted values do.
- Implement analytic Hessians for the RooFunctor test functors, so
  that test exercises the new externally-bound-functor support
  end-to-end.
- Fix degenerate test models whose Hessians were singular, making the
  error comparison meaningless: the Gaussian/RooFormulaVar models
  depended only on mu + shift (shift is now constant), and the
  Bernstein pdf was invariant under a common coefficient rescaling
  (c0 is now constant).
- Improve the conditioning of ill-defined fits: RooLandau3 (sl = 10,
  comparable to the observable window) generates 10k events via a new
  nEvents argument, and the RooFunctor Gaussian starts at mu = 2,
  sigma = 1.5 so the data constrain all parameters and the fitted
  values are away from zero.
- Skip constant parameters when randomizing initial values; factory
  constants have an infinite range, so randomization pushed them to
  +-inf.

🤖 Done with the help of AI
Clad computes Hessians in reverse-over-forward mode: the forward pass
needs a `<name>_pushforward` for the opaque roo_outer_wrapper call, and
the reverse pass over that call needs the matching
`<name>_pushforward_pullback`. Without them, Hessian generation failed
for any model containing a RooONNXFunc.

Emit both from RooONNXFunc::initialize(), into the same
clad::custom_derivatives namespace as the existing custom pullback.
Both are plain C++ built on the already-emitted exact gradient
pullback, so no additional Clad differentiation happens:

- roo_outer_wrapper_pushforward returns the function value plus the
  directional derivative grad . d_input, both exact.
- roo_outer_wrapper_pushforward_pullback needs second derivatives only
  as the Hessian-vector product H . d_input. SOFIE emits Clad pullbacks
  for its operators but no pushforwards, so differentiating the model
  code again would silently fall back to numerical differentiation of
  the value. Instead, evaluate the product as a central finite
  difference of the exact generated gradient along the normalized
  tangent direction. The step size (~cbrt of the float machine epsilon,
  scaled to the input magnitude) balances the float-precision noise of
  the SOFIE gradient against the truncation error; the result agrees
  with a float64 PyTorch reference Hessian to ~2e-4 relative accuracy.

The existing test models cannot validate this: a ReLU MLP is
piecewise-linear in its inputs, so its input Hessian vanishes almost
everywhere. Make the activation function a model-class argument in
create_onnx_model.py and add tanh variants of both models, together
with a torch.autograd.functional.hessian reference saved in the RooFit
parameter ordering (SOFIE generates Tanh as a plain std::tanh loop,
which Clad reverse-differentiates natively).

New tests: the single-tensor 10x10 Hessian, the two-tensor 15x15
Hessian including the cross-tensor blocks, and the square of a
RooONNXFunc via RooProduct, whose H(f^2) = 2 (f H + grad grad^T)
exercises the primal-value adjoint term that stays dormant when the
ONNX function is the top-level function.

🤖 Done with the help of AI
Now that the SOFIE-generated inference code supports Clad forward-mode
differentiation, replace the finite-difference-of-gradient
approximation for RooONNXFunc second derivatives with exact
Hessian-vector products.

RooONNXFunc::initialize() now additionally requests a forward-mode
derivative of the inner wrapper, keeps a persistent zero-weight
tangent Session (via the generated SetWeightsToZero()), declares a
directional-derivative wrapper around the generated pushforward, and
reverse-differentiates it with respect to both the inputs and the
tangent direction. The resulting pullback yields the exact
Hessian-vector product and the gradient in a single reverse pass, so
the emitted pushforward_pullback is both exact and cheaper than the
previous three-gradient finite-difference stencil.

The emitted second-order pullback constructs fresh adjoint Session
objects on every call, working around a Clad bug where generated
second-order pullbacks do not restore intermediate adjoint state.
Differentiating with respect to the tangent direction as well works
around a Clad custom-derivative lookup limitation with non-varied
pointer arguments (and its adjoint is the gradient, which is needed
anyway).

Measured accuracy against the float64 PyTorch reference Hessian
improves from 2.3e-6 (finite differences) to 1.9e-9, so the Hessian
tolerances in testRooONNXFunc are tightened from 1e-3 to 1e-5,
matching the gradient checks.

🤖 Done with the help of AI
RooMinimizer used to reach Minuit2 through the ROOT::Math::Minimizer plugin
interface, while at the same time using Minuit2 directly as a library for
things the generic interface can't express (vanishing second derivatives,
the multiprocess gradient FCN). This commit removes that split: when the
minimizer type is "Minuit2" (the default), RooMinimizer now uses the Minuit2
library interface directly (ModularFunctionMinimizer, MnHesse, MnMinos,
MnContours, MnUserParameterState, FunctionMinimum). Other minimizer types
still go through ROOT::Math::Minimizer.

To make that possible, RooAbsMinimizerFcn now implements
ROOT::Minuit2::FCNBase. RooMinimizerFcn provides gradient, Hessian and the
vanishing-second-derivative mask through that interface, and the FCNBase
adapter that MinuitFcnGrad carried for the multiprocess gradient is folded
into MinuitFcnGrad itself. LikelihoodGradientJob gets the FCN from the
RooMinimizer instead of from the Minuit2Minimizer wrapper.

The internal FitResult is now the single source of truth for fit results on
both paths; it additionally stores the global correlation coefficients, and
exposes the covariance matrix status, which xRooFit uses instead of asking
the (no longer existing) ROOT::Math::Minimizer. FitterInterface::GetMinimizer
returns nullptr when Minuit2 is used directly.

Behaviour is kept as in Minuit2Minimizer, including its status codes, print
levels, strategy customization from the Minuit2 extra options, and the
re-minimization when Minos finds a new minimum. Differences: Hesse now
applies the current error definition, the "Minuit2Minimizer:" std::cout
lines are now prefixed "RooMinimizer: Minuit2", and the Fumili algorithm,
which never worked with RooFit FCNs, is rejected with an explicit error.

🤖 Done with the help of AI
Now that RooMinimizer drives Minuit2 directly, the entry points that were
added to the Minuit2Minimizer wrapper only for RooFit have no users left:

  - Minuit2Minimizer::SetFCN(unsigned int, std::unique_ptr<FCNBase>) and
    both GetFCN() overloads, which let RooFit's multiprocess gradient FCN
    bypass the ROOT::Math function interface.
  - Minuit2Minimizer::SetSecondDerivativeAlwaysVanishesFunc(), which just
    forwarded to the FCNAdapter setter. RooMinimizerFcn now overrides
    FCNBase::SecondDerivativeAlwaysVanishes() itself. The FCNAdapter setter
    stays, since it is useful independently of RooFit.

The NumericalDerivator class is a RooFit component that happened to live in
Minuit2: nothing in Minuit2 uses it, its only user is LikelihoodGradientJob.
It moves next to that user as RooFit::TestStatistics::NumericalDerivator,
with the original authorship notice preserved.

🤖 Done with the help of AI
A collection of small simplifications that fall out of RooMinimizer using
Minuit2 directly and the FCN classes implementing ROOT::Minuit2::FCNBase:

  - RooMinimizerFcn::initMinimizer() no longer passes the Hessian through
    ROOT::Math::Minimizer::SetHessianFunction(). Only the Minuit2 wrapper
    implements that method, and Minuit2 doesn't go through this path
    anymore, so it was a no-op. The Hessian is now computed directly in
    the FCNBase::Hessian() override.
  - initMinimizer() is a plain member of RooMinimizerFcn instead of a
    virtual on RooAbsMinimizerFcn, since the multiprocess MinuitFcnGrad
    could only throw there. RooMinimizer checks the function type when it
    creates a ROOT::Math::Minimizer.
  - The raw-pointer evaluation and gradient methods of MinuitFcnGrad were
    only wrappers called by the FCNBase overrides; the implementation
    moves into the overrides.
  - The RooMinimizer::FitResult constructor from a FitConfig and the
    fallback that used it in calculateHessErrors() were unreachable:
    fitFCN() rejects functions without parameters before a result is
    created, and hesse() refuses to run without a result.
  - The FCN classes and LikelihoodGradientJob access the RooMinimizer
    members directly instead of going through the fitter() compatibility
    shim, which allocated a FitterInterface on every call.

🤖 Done with the help of AI
@guitargeek
guitargeek force-pushed the roofit_clad_hessians branch from e6d809c to 09938e1 Compare October 13, 2026 11:21

This branch has not been deployed

No deployments
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Projects

None yet

Development

Successfully merging this pull request may close these issues.

1 participant