Convergence failure using fixed_poi_fit

[APMDSLHC]

Summary

Dear pyhf developers, first of all thank you for your work, it is really helpful for the SModelS collaboration.

Second, I am trying to compute likelihoods for various values of the parameter of interest, and for some them, the computation fails below a certain value of the parameter of interest. Here for instance, when computing pyhf.infer.mle.fixed_poi_fit of ATLAS-SUSY-2019-08 with the signal strength below around -1.08, I get the error 'Positive directional derivative for linesearch'. (Please keep in mind that in SModelS we rescale the signals according to the signal strength, so we call the fixed_poi_fit method with poi_val = 1.) In the minimal test-case below, it works for mu = -1.08 but fails for mu = -1.09.

The first guess was that the yields of the main model become negative below -1.08, but I doesn't seem to be the case because there is no change of sign between the yields obtained using mu = -1.08 and mu = -1.09. (I tried 'by hand' and with model.expected_actualdata but the results didn't match. I probably misused the method so I opened another question thread: #1885.)

I tried with different backends (numpy, tensorflow, pytorch) and optimizers (scipy, minuit) but none worked.
I also tried to change the tolerance and epsilon values for the optimizer but it didn't improved anything.

The last guess would be that modifiers parameters have to go outside of their bounds to fit correctly the yields that were not supposed to be that low (even if these yields are far to 0), but I don't know if it makes sense (aren't the parameters supposed to stop at their bounds even if the minimum is outside?) and how to check that.

In the test-case below I start the fit with different initial parameters than the suggested ones because it fails using the latter for mu = -1.08, but it works if I take instead the best fit parameters obtained with mu = 0. Can the issue be that the fit didn't converge because of the initial parameters? If so, how do I choose them?

Thanks in advance for helping me understand what goes wrong.

Finally, I searched for similar questions and didn't find anything, sorry if it has already been answered.

OS / Environment

ID=fedora
VERSION_ID=35
VERSION_CODENAME=""
PLATFORM_ID="platform:f35"
PRETTY_NAME="Fedora Linux 35 (Workstation Edition)"
ANSI_COLOR="0;38;2;60;110;180"
LOGO=fedora-logo-icon
CPE_NAME="cpe:/o:fedoraproject:fedora:35"
HOME_URL="https://fedoraproject.org/"
DOCUMENTATION_URL="https://docs.fedoraproject.org/en-US/fedora/f35/system-administrators-guide/"
SUPPORT_URL="https://ask.fedoraproject.org/"
BUG_REPORT_URL="https://bugzilla.redhat.com/"
REDHAT_BUGZILLA_PRODUCT="Fedora"
REDHAT_BUGZILLA_PRODUCT_VERSION=35
REDHAT_SUPPORT_PRODUCT="Fedora"
REDHAT_SUPPORT_PRODUCT_VERSION=35
PRIVACY_POLICY_URL="https://fedoraproject.org/wiki/Legal:PrivacyPolicy"
VARIANT="Workstation Edition"
VARIANT_ID=workstation

Steps to Reproduce

import pyhf
pyhf.set_backend(b"pytorch",'scipy')
import json
import jsonpatch

# ATLAS-SUSY-2019-08 json file
with open("./BkgOnly.json", "r") as f:
    bkg = json.load(f)

# Signals produced by SModelS
nsig = [1.4349825452486582, 1.0533138982733217, 1.1158810930352483, 
        0.8276376394629201, 0.8551014982035671, 1.6913336713087228, 
        0.719039660507198, 0.9461983338046028, 1.5376094024329356]

mu=-1.09
nsig_rescaled = [sig*mu for sig in nsig]

# The patch produced by SModelS
patch = [
    {
      "op": "add",
      "path": "/channels/5/samples/0",
      "value": {
        "data": [
          nsig_rescaled[0],
          nsig_rescaled[1],
          nsig_rescaled[2]
        ],
        "modifiers": [
          {
            "data": None,
            "type": "normfactor",
            "name": "mu_SIG"
          },
          {
            "data": None,
            "type": "lumi",
            "name": "lumi"
          }
        ],
        "name": "bsm"
      }
    },
    {
      "op": "add",
      "path": "/channels/6/samples/0",
      "value": {
        "data": [
          nsig_rescaled[3],
          nsig_rescaled[4],
          nsig_rescaled[5]
        ],
        "modifiers": [
          {
            "data": None,
            "type": "normfactor",
            "name": "mu_SIG"
          },
          {
            "data": None,
            "type": "lumi",
            "name": "lumi"
          }
        ],
        "name": "bsm"
      }
    },
    {
      "op": "add",
      "path": "/channels/7/samples/0",
      "value": {
        "data": [
          nsig_rescaled[6],
          nsig_rescaled[7],
          nsig_rescaled[8]
        ],
        "modifiers": [
          {
            "data": None,
            "type": "normfactor",
            "name": "mu_SIG"
          },
          {
            "data": None,
            "type": "lumi",
            "name": "lumi"
          }
        ],
        "name": "bsm"
      }
    },
    {
      "op": "remove",
      "path": "/channels/4"
    },
    {
      "op": "remove",
      "path": "/channels/3"
    },
    {
      "op": "remove",
      "path": "/channels/2"
    },
    {
      "op": "remove",
      "path": "/channels/1"
    },
    {
      "op": "remove",
      "path": "/channels/0"
    }
]

# Trying to reproduce the issue
llhdSpec = jsonpatch.apply_patch(bkg, patch)
msettings = {'normsys': {'interpcode': 'code4'}, 'histosys': {'interpcode': 'code4p'}}
workspace = pyhf.Workspace(llhdSpec)
model = workspace.model(modifier_settings=msettings)

init, n_ = pyhf.infer.mle.fixed_poi_fit(0.0, workspace.data(model), model, return_fitted_val=True, maxiter=200)
init[model.config.poi_index] = 1.

_, nllh = pyhf.infer.mle.fixed_poi_fit(1., workspace.data(model), model, return_fitted_val=True, init_pars=init.tolist(), maxiter=200)
print(nllh)

File Upload (optional)

The ATLAS-SUSY-2019-08 background only json file: BkgOnly.zip

Expected Results

I expect the pyhf.infer.mle.fixed_poi_fit method to converge for mu = -1.09 since it did for mu = -1.08 and the main model yields didn't change their signs.

Actual Results

ERROR in mixins._internal_minimize() in 63:      fun: nan
     jac: array([-1.60974331e+03, -1.26528416e+03, -3.62095746e+00, -2.96069603e+00,
       -9.31107394e+00, -5.92875301e+00, -3.95191172e+00, -1.09481013e+00,
        1.98526130e+00,  0.00000000e+00,  0.00000000e+00, -2.61623611e-01,
       -2.91744713e-01,  4.78787438e-02, -1.14895247e-01, -2.19530725e+00,
       -3.08090779e-01, -2.01610681e-02,  5.56939377e-02, -1.23923733e+01,
        1.82827884e+01, -2.31160044e+00, -7.01946576e+00,  8.12185110e+01,
       -7.43248366e+01, -2.18593652e+00,  4.20367681e-07, -3.01232848e+00,
        1.34875544e+01,  7.74696699e+00, -1.03325355e+01, -1.97692861e+00,
       -5.14075418e+01, -5.01545550e+01,  1.20103830e+01, -3.57878791e+00,
       -2.30993275e+01,  4.44221067e+00, -7.52509470e+00, -3.66384988e+00,
       -1.03347115e+01,  1.11810778e-05,  3.51899237e-05, -2.34341867e-05,
        1.19292659e-05, -3.13825203e-08,  3.03152521e+00, -1.28334787e+00,
       -4.32407676e+00, -5.09903434e+00,  2.04141428e-01, -1.89871618e+01,
        1.04031278e+00,  1.73120504e+00,  4.60851018e+00,  5.19418345e-01,
        1.14268112e+01, -3.57270393e+00, -5.41869452e+00, -3.87590570e+00,
        9.41609612e-01,  4.02259363e+00, -6.61617101e-01,  2.74473331e-01,
        5.59316983e+00, -3.15191816e+00, -8.08768057e+00, -2.10667595e+01,
       -1.95631376e+00, -6.06152697e+00,  6.26285962e-02,  5.39568481e+00,
       -4.82324499e-01,  2.30374024e-05,  1.47277931e+01,  2.33850727e+00,
        6.49301932e+00, -4.19677148e+00,  3.98913374e+00,  9.73124241e-01,
        2.47022211e-07, -3.32185078e+00, -4.90462987e+01,  1.11342365e+02,
        1.03764876e+02,  5.18201994e+01,  3.17514693e+01,  2.58190318e+02,
        3.65515541e+02,  9.87664763e+01,  3.80348579e+02,  6.22617331e+01,
       -8.53710176e+01,  2.60184869e+02,  7.60798866e+00,  1.06295164e+01,
       -9.39105087e+00,  1.85471213e+01,  4.61838002e+02,  1.26344259e+01,
        7.19938698e+00,  3.08362251e+00,  5.67613404e+00,  2.22176951e+01,
       -1.78534368e+00, -2.21602320e+00,  1.59233794e+03, -1.15381966e-04,
       -6.85014709e-04, -1.75581079e+00, -2.23552249e+00, -4.23820705e+00,
        1.84184018e-03, -5.41821918e-04, -1.22347766e-03,  1.07853930e-03,
       -2.26055107e-03])
 message: 'Positive directional derivative for linesearch'
    nfev: 408
     nit: 42
    njev: 38
  status: 8
 success: False
       x: array([ 1.00000000e+00,  1.00000000e+00,  9.58886570e-01,  9.97696558e-01,
        1.08981802e+00, -3.60593775e-03, -1.14103890e-02, -4.96081000e-04,
        8.75706148e-03,  0.00000000e+00,  0.00000000e+00,  3.97060272e-05,
        1.86660588e-03,  3.71200829e-04,  3.34510586e-06, -7.00390783e-04,
       -3.00302912e-05,  1.41270604e-04,  1.67682598e-04,  5.26169671e-03,
        4.08688609e-03,  5.70512574e-03, -2.26821170e-04,  1.74170870e-02,
       -8.90834102e-03, -1.31640865e-02,  1.48830609e-08, -1.19241950e-02,
        1.00610484e-01,  4.34449496e-01, -2.38917682e-02, -2.63947287e-01,
       -2.75862486e-02, -2.29484271e-01, -7.76770744e-02,  6.39863578e-03,
        1.57811398e-01, -1.21716944e-02,  1.30392004e-01, -2.37317059e-02,
        1.54718811e-02, -5.50983002e-09, -2.70633098e-08,  1.18182744e-08,
        2.38744162e-08, -2.12351237e-08, -4.26668255e-02, -7.09514843e-04,
        6.23132923e-02, -3.17626867e-02, -2.91056610e-03,  2.43762272e-01,
       -7.17585624e-02,  2.37442108e-02, -4.70766739e-02, -4.23360101e-02,
       -7.14426477e-02,  4.25623821e-02,  5.12821275e-02,  4.77077237e-02,
        8.57077165e-03, -1.76747770e-02,  1.93550820e-06,  1.45023060e-04,
        3.78501837e-04, -1.57105215e-02,  3.42261725e-01, -3.34979529e-01,
        1.65563110e-02,  1.66380135e-01,  1.10685827e-01, -8.53969855e-02,
        9.02734668e-05,  2.85312971e-08, -4.63458866e-02,  3.48438356e-02,
       -1.66824223e-01, -3.59695286e-02,  8.09042950e-02,  1.27854537e-02,
       -5.30580051e-08,  5.00428015e-02,  2.34490607e-02,  5.49906272e-02,
       -9.12593032e-02,  7.49021757e-02,  1.18206572e-03,  8.92087399e-01,
        2.82968811e-02,  1.17650301e-01,  4.56959847e-01,  2.55444945e-01,
        2.17030067e-01,  2.65345683e+00,  8.24716469e-03,  6.86620028e-03,
       -3.55064521e-03, -1.47798591e-02,  6.97768957e-01,  4.51961102e-03,
        4.05470915e-03,  7.47849933e-04, -6.90514422e-03, -3.42759564e-04,
        9.66281806e-01,  1.01714164e+00,  1.01506935e+00, -5.29982261e-05,
       -2.17083661e-04,  1.01093401e+00,  9.97608185e-01,  9.97469708e-01,
       -8.15874543e-05,  6.83610316e-05, -2.89910275e-05,  1.56607461e-04,
        2.74800236e-05])
Traceback (most recent call last):
  File "/home/pascal/.local/lib/python3.9/site-packages/pyhf/optimize/mixins.py", line 61, in _internal_minimize
    assert result.success
AssertionError

---------------------------------------------------------------------------
AssertionError                            Traceback (most recent call last)
File ~/.local/lib/python3.9/site-packages/pyhf/optimize/mixins.py:61, in OptimizerMixin._internal_minimize(self, func, x0, do_grad, bounds, fixed_vals, options, par_names)
     60 try:
---> 61     assert result.success
     62 except AssertionError:

AssertionError: 

During handling of the above exception, another exception occurred:

FailedMinimization                        Traceback (most recent call last)
Input In [23], in <cell line: 3>()
      1 bounds = model.config.suggested_bounds()
----> 3 best_fit, nllh = pyhf.infer.mle.fixed_poi_fit( 1., workspace.data(model), model, return_fitted_val=True, 
      4                                        init_pars = initpars, par_bounds = bounds, maxiter=200 )
      5 print(nllh)

File ~/.local/lib/python3.9/site-packages/pyhf/infer/mle.py:202, in fixed_poi_fit(poi_val, data, pdf, init_pars, par_bounds, fixed_params, **kwargs)
    199 init_pars[pdf.config.poi_index] = poi_val
    200 fixed_params[pdf.config.poi_index] = True
--> 202 return fit(data, pdf, init_pars, par_bounds, fixed_params, **kwargs)

File ~/.local/lib/python3.9/site-packages/pyhf/infer/mle.py:131, in fit(data, pdf, init_pars, par_bounds, fixed_params, **kwargs)
    124 # get fixed vals from the model
    125 fixed_vals = [
    126     (index, init)
    127     for index, (init, is_fixed) in enumerate(zip(init_pars, fixed_params))
    128     if is_fixed
    129 ]
--> 131 return opt.minimize(
    132     twice_nll, data, pdf, init_pars, par_bounds, fixed_vals, **kwargs
    133 )

File ~/.local/lib/python3.9/site-packages/pyhf/optimize/mixins.py:183, in OptimizerMixin.minimize(self, objective, data, pdf, init_pars, par_bounds, fixed_vals, return_fitted_val, return_result_obj, return_uncertainties, return_correlations, do_grad, do_stitch, **kwargs)
    180         par_names[index] = None
    181     par_names = [name for name in par_names if name]
--> 183 result = self._internal_minimize(
    184     **minimizer_kwargs, options=kwargs, par_names=par_names
    185 )
    186 result = self._internal_postprocess(
    187     result, stitch_pars, return_uncertainties=return_uncertainties
    188 )
    190 _returns = [result.x]

File ~/.local/lib/python3.9/site-packages/pyhf/optimize/mixins.py:64, in OptimizerMixin._internal_minimize(self, func, x0, do_grad, bounds, fixed_vals, options, par_names)
     62 except AssertionError:
     63     log.error(result, exc_info=True)
---> 64     raise exceptions.FailedMinimization(result)
     65 return result

FailedMinimization: Positive directional derivative for linesearch

pyhf Version

pyhf verion: 0.6.3
scipy version: 1.8.0
minuit version: 2.11.2
pythorch version: 1.11.0+cu102
tensorflow version: 2.9.0
numpy version: 1.21.4

Code of Conduct

• I agree to follow the Code of Conduct

[alexander-held]

I think the difference you see between -1.08 and -1.09 may not be particularly meaningful. When using your script, the second fit fails for me with -1.08 as well.

With mu_SIG=1 fixed, your model predicts negative yields in two bins (for both the -1.08 and -1.09 scaling). I have not investigated the model here more closely, but I assume what happens is that the fit starts out in a region of parameter space where your likelihood is NaN and does not manage to leave it. Appropriate initialization could help, but it is not clear to me how to generalize this. I think the problem here is that the model is ill-defined to begin with.

#1845 is somewhat related here. clip_sample_data=0.0 allows the second fit to succeed, but presumably does not do what you actually want (because it would set signal in those bins back to 0 where it is below 0). You can also use clip_bin_data=1e-6 (note: zero won't work here). I would not generally recommend this as a "solution" for what you see, it's a very ad-hoc method to work around what in my opinion should be fixed on a per-model level when it arises to ensure that the models make sense. What is the reason for testing large negative values of the signal strength where this comes up?

[APMDSLHC]

Thank you for the answer.

I think the difference you see between -1.08 and -1.09 may not be particularly meaningful. When using your script, the second fit fails for me with -1.08 as well.

It is odd, I tried again and -1.08 works for me.

The reason for testing large negative values is that I wanted to visualize the product of likelihoods coming from different analyses. For some analyses, the signal strength that maximises the likelihood is for such negative values. So it can be an issue when one of them cannot be computed while another is at its maximum.
Changing the initial parameters seems to works, so a solution to my problem might be to set them equal to the best fit of the nearest mu that converged.

Of course, I don't want to compute likelihoods that have negative yields (that's why I opened the discussion about model.expected_actualdata), but the main model yields can be negative when computed with the suggested initial parameters while they can be positive when computed with the best fit parameters. Anyway, that's not an issue, it just makes it a bit harder to perform a 'security check' to ensure that the yields are positive using the best fit if the fit cannot be done because the yields are negative with the suggested initial parameters.

[alexander-held]

It is odd, I tried again and -1.08 works for me.

I used numpy instead of pytorch, which might have made the difference. The path taken by the minimizer will generally vary between different machines, so it's also possible that for edge cases the fit will fail on one machine and succeed on the next with the same minimizer configuration.

Changing the initial parameters seems to works, so a solution to my problem might be to set them equal to the best fit of the nearest mu that converged.

That sounds like a good idea to me. Another thing that could help (assuming that the signal normalization is typically responsible for the initial issues) is initializing the signal normalization at 0.

[APMDSLHC]

That sounds like a good idea to me. Another thing that could help (assuming that the signal normalization is typically responsible for the initial issues) is initializing the signal normalization at 0.

I'm not sure to understand. I am using fixed_poi_fit with poi_val = 1, so I can't set the signal normalization ('mu_SIG', right?) to 0, it would create a conflict with the fixed poi_val, wouldn't it?

[alexander-held]

You're completely correct, I forgot that the context here was a fit with POI = 1 fixed. My suggestion doesn't do anything here then. Sorry for the confusion!

[APMDSLHC]

No problem, thanks for your help!