Fit succeeds with scipy, but fails with minuit

[gollumben]

Dear pyhf developers,

for the attached workspace the fit with scipy succeeds, but fails with minuit. The workspace was imported from HistFitter, where the fit converges without issues. I have tried a variety of changes, including the tolerance between 1e-3 and 200 and the strategy.

For the minimum working example below, the code exists fine for useMinuit = 0, but gives an error pyhf.exceptions.FailedMinimization: Optimization failed. Estimated distance to minimum too large for useMinuit = 1

Do you have any tip what could be going wrong? Some of the NPs have different types, could that be an issue?
Cheers,
Ben

Versions: pyhf - master; iminuit - 2.8.2

Minimum usage examples:

import pyhf
import json

jsonObj = None
with open("NormalMeasurement_combined.txt") as iF:
                jsonObj = json.load(iF)
workspace = pyhf.Workspace(jsonObj)

useMinuit = 1
if useMinuit:
        pyhf.set_backend("numpy", "minuit") ## get uncertainties
        pyhf.set_backend(pyhf.tensorlib, pyhf.optimize.minuit_optimizer(tolerance=1,verbose=2,strategy=1))

model = workspace.model(measurement_name="NormalMeasurement")
data = workspace.data(model)

result = pyhf.infer.mle.fixed_poi_fit(1.0, data, model, return_uncertainties=True)

NormalMeasurement_combined.txt

[kratsg]

ping @alexander-held . Not sure if you've tried changing the tolerance down to 1e-4 or something like that. Looking at your workspace, you seem to have most systematics with very little impact.

[matthewfeickert]

There have been some other documented instances of minuit failing where it shouldn't (c.f. Issue #1591). This is something that we need to actively work on.

@gollumben If you just want something that will get you the minuit uncertainties so you can keep moving though for the time being, and you don't have any hesitancy about the fit parameter values from the other backends optimizer results then you could use those fit results as the initial parameters for minuit like this:

import json

import numpy as np
import pyhf

if __name__ == "__main__":
    backend_name = "numpy"
    pyhf.set_backend(backend_name)

    with open("NormalMeasurement_combined.json") as input_file:
        workspace = pyhf.Workspace(json.load(input_file))

    model = workspace.model(measurement_name="NormalMeasurement")
    data = workspace.data(model)
    test_poi = 1.0

    fit_pars = pyhf.infer.mle.fixed_poi_fit(
        test_poi, data, model, return_uncertainties=True
    )
    print(f"{fit_pars=}")
    np.savetxt(f"fit_pars_{backend_name}.txt", fit_pars, delimiter=", ")

    pyhf.set_backend(
        backend_name,
        pyhf.optimize.minuit_optimizer(tolerance=1, verbose=2, strategy=1),
    )

    fit_pars_minuit = pyhf.infer.mle.fixed_poi_fit(
        test_poi,
        data,
        model,
        init_pars=fit_pars.tolist(),
        return_uncertainties=True,
    )
    print(f"{fit_pars_minuit=}")
    np.savetxt(f"fit_pars_minuit_{backend_name}.txt", fit_pars_minuit, delimiter=", ")

    pars_difference = fit_pars - fit_pars_minuit[:, 0]
    np.savetxt(f"pars_difference_{backend_name}.txt", pars_difference, delimiter=", ")

If I use JAX

backend_name = "jax"

I can complete both fits in under a minute on my laptop on CPU (and in 30 seconds on GPU).

[gollumben]

Hi @kratsg and @matthewfeickert,
again, thank you very much for your input! I'll use some bullet points to structure my answer.

0. Unfortunately reducing the tolerance=1e-4 did not lead to any conversion :/.

1. Using the output of scipy as initial values of minuit did the trick: all of my fits converge now and I have uncertainties. Thank you so much! :) However, the fit results are slightly different w.r.t. using HistFitter standalone. I suppose this is expected and all values agree within uncertainies, but I thought this might be worth mentioning. If I can help you somehow to debug the issues you seem to be having with minuit, please let me know.

2. Also, thank you for pointing me to the jax acceleration! I suppose I'd turn on GPU usage by doing following?

import jax
jax.devices()

3. If I may draw your attention to one further "issue" I encountered: it's about the documentation of the combination strategy of workspaces. The options for join are documented, yet, it is not documented what exactly none, outer, left outer, right outer stand for. While I quickly managed to find this PR feat: Combine is taught about formal data (outer) join methods #752, which explains that outer etc. follow the SQL syntax, I feel like a reference to the SQL syntax or the PR or a short explanation would be useful.

Cheers,
Ben

[alexander-held]

Hi, I have two small follow-up points regarding what you mention.

0. The fit fails presumably since it hits regions of parameter space where the NLL is NaN and then does not recover afterwards. The tolerance sets the cut-off condition for minimization, so in this context unfortunately that cannot help much since the setup is far away from convergence regardless. I think it would be interesting to determine what causes the NaN region of parameter space. I have also observed such behavior in ROOT previously. It may very well be that the ROOT model would also be NaN for the same parameter values, but that the fit takes a slightly different route to the minimum and misses that region. From a brief look at the workspace I did not manage to isolate any obvious element responsible for the NaN result, but I will try to have a more detailed look. If you manage to create a simplified example (less channels / samples / systematics) that still fails consistently, it would be useful to try and track down potentially problematic parts of the setup.
1. Differences in best-fit results are generally expected, the important part is that they are not very large. Especially with a high tolerance, repeated fits within the same framework from different initial conditions can give slight variations in the parameter results. If the differences are larger than ~1% or so, and do not decrease with a lower tolerance, then this is something we should investigate in detail.

[alexander-held]

Another strategy that may help with convergence is randomizing initial parameter values at the start of the fit. This will cause the minimization to take a different route, and potentially succeed. For that you may want to use a fixed seed so your fits fail / succeed consistently.

[alexander-held]

Small observation: using the script provided at the top, I see a difference between the pyhf master (reproducing the failure) and pyhf 0.6.3, which hits a problematic region of parameter space briefly early on but then recovers. This is with iminuit 2.8.4.

This may be related to the modifier changes introduced recently, the reordering might have introduced some tiny floating point changes which are enough to alter the path to the minimum such that 0.6.3 succeeds.

[alexander-held]

I managed to isolate a problem in this workspace, and expect that it's likely the full explanation for what you see. It seems to come down to the use of histosys modifiers that have normalization effects (not just shape), but due to the way the HistFactory extrapolation works, this means that strong pulls can result in negative yield predictions. The exponential extrapolation for normsys modifiers protects against this, and I would recommend using it especially for single-bin cases, like the one causing an issue here.

To get started, it is useful to have parameters where the NaN logpdf value happens. One way to do this is to add a breakpoint to Model.logpdf if np.isnan(result). From there, the problematic parameters (and data) can be extracted. To simplify things, it's then helpful to prune out as much as possible from the model. You can find a script in this gist that does exactly that, leaving you with a significantly simplified model in the end: single channel / sample and 7 parameters:

channels: ['bin3_cuts']
samples: ['zjetsBoosted']
parameter names and values:
Sherpa_muR_zjets = -1.81622979
WTag_BGSF_Gammajet_Stat = -3.84084271
WTag_BGSF_Propagated_AllOthers = -3.39250537
WTag_JetTagSF_Dijet_Modelling = -3.6207387
WTag_JetTagSF_Radiation = -4.60027512
WTag_SigSF_BinVariation = -3.41855226
staterror_bin3_cuts[0] = 0.934653809

These are multiple histosys modifiers (and a staterror), and jointly they are pulling the model prediction to a negative value.

The resulting model can be simplified even further by removing measurement config pieces that have no impact here. Here is what remains:

import pyhf

spec = {
    "channels": [
        {
            "name": "bin3_cuts",
            "samples": [
                {
                    "data": [3.490729331970215],
                    "modifiers": [
                        {
                            "data": [0.277386256344252],
                            "name": "staterror_bin3_cuts",
                            "type": "staterror",
                        },
                        {
                            "data": {
                                "hi_data": [3.8195865154266357],
                                "lo_data": [3.2028777599334717],
                            },
                            "name": "Sherpa_muR_zjets",
                            "type": "histosys",
                        },
                        {
                            "data": {
                                "hi_data": [3.546879291534424],
                                "lo_data": [3.3108692169189453],
                            },
                            "name": "WTag_BGSF_Gammajet_Stat",
                            "type": "histosys",
                        },
                        {
                            "data": {
                                "hi_data": [3.530710220336914],
                                "lo_data": [3.3270678520202637],
                            },
                            "name": "WTag_BGSF_Propagated_AllOthers",
                            "type": "histosys",
                        },
                        {
                            "data": {
                                "hi_data": [3.5470385551452637],
                                "lo_data": [3.3109240531921387],
                            },
                            "name": "WTag_JetTagSF_Dijet_Modelling",
                            "type": "histosys",
                        },
                        {
                            "data": {
                                "hi_data": [3.514037847518921],
                                "lo_data": [3.343658924102783],
                            },
                            "name": "WTag_JetTagSF_Radiation",
                            "type": "histosys",
                        },
                        {
                            "data": {
                                "hi_data": [3.5140371322631836],
                                "lo_data": [3.3436598777770996],
                            },
                            "name": "WTag_SigSF_BinVariation",
                            "type": "histosys",
                        },
                    ],
                    "name": "zjetsBoosted",
                }
            ],
        }
    ],
    "measurements": [
        {"config": {"parameters": [], "poi": ""}, "name": "NormalMeasurement"}
    ],
    "observations": [{"data": [3.0], "name": "bin3_cuts"}],
    "version": "1.0.0",
}

ws = pyhf.Workspace(spec)
model = ws.model(poi_name=None)
data = ws.data(model)

# parameter values from other script
pars = [
    -1.81622979,
    -3.84084271,
    -3.39250537,
    -3.6207387,
    -4.60027512,
    -3.41855226,
    0.934653809,
]
print(model.expected_data(pars, include_auxdata=False))

resulting in

[-0.10138098]

which is a negative (unphysical) yield prediction.

The setup in the full workspace causing a NaN may be slightly different from this, but I expect it comes down to the same. In principle other samples can cancel the negative prediction from this sample, but I expect jointly it is still something that plays out very similarly.

Since the final fit results are unlikely to include such strong pulls, this is likely just a minimization issue that arises while the minimizer explores the parameter space. Different initialization can avoid this part of the parameter space and cause the result to be ok regardless. For stability reasons I would recommend the use of normsys modifiers though to protect against this issue.

[kratsg]

As @alexander-held mentioned, this workspace seems to be impacted by the issue of negative yields. This is fixed in #1845 and will be in v0.7.0. There is a release candidate 0.7.0rc1 out now if you would like to try. You can clip your model's expected data yields during minimization by building it like so:

model_clipped = workspace.model(clip_bin_data=0.0)

and that should allow minuit to minimize properly (although note that by doing this: you will introduce discontinuities in your likelihood). If this works for you, accept either Alex's answer or my reply so we can mark this as answered/resolved.