Limit determination fails to converge

[nhartman94]

Hello @matthewfeickert @alexander-held !

Daarii is doing some fits, and is getting an intermittent error depending on what variable we fit on. Could you guide us to how to interface w/ scipy.optimize via cabinetry to converge?

DEBUG - cabinetry.fit - Signal_norm = 0.1000, expected -2 sigma CLs = 0.0500 (cached)
DEBUG - cabinetry.fit - Signal_norm = 0.1000, expected -2 sigma CLs = 0.0500 (cached)
ERROR - cabinetry.fit - CLs values at 0.1000 and 0.1000 do not bracket CLs=0.0500, try a different starting bracket
Traceback (most recent call last):
  File "/lustre/fs22/group/atlas/dbattulga/MyAnalysis/cabinetry-tutorials/example.py", line 306, in <module>
    limit_results = cabinetry.fit.limit(model, data, confidence_level=0.95, poi_tolerance=0.01)
  File "/cvmfs/sft.cern.ch/lcg/views/LCG_104/x86_64-centos7-gcc12-opt/lib/python3.9/site-packages/cabinetry/fit/__init__.py", line 959, in limit
    res = scipy.optimize.root_scalar(
  File "/cvmfs/sft.cern.ch/lcg/views/LCG_104/x86_64-centos7-gcc12-opt/lib/python3.9/site-packages/scipy/optimize/_root_scalar.py", line 275, in root_scalar
    r, sol = methodc(f, a, b, args=args, **kwargs)
  File "/cvmfs/sft.cern.ch/lcg/views/LCG_104/x86_64-centos7-gcc12-opt/lib/python3.9/site-packages/scipy/optimize/_zeros_py.py", line 784, in brentq
    r = _zeros._brentq(f, a, b, xtol, rtol, maxiter, args, full_output, disp)
ValueError: f(a) and f(b) must have different signs

Thank you so much!!

[alexander-held]

Hi @nhartman94! What happens here is that cabinetry sets a default lower bound of 0.1 for the interval it probes to find all the CLs=0.05 crossings. The default upper bound is the upper POI bound set in the workspace (which defaults to 10 if not set). That means by default the interval being searched is [0.1, 10]. The algorithm being used to identify the crossings, scipy.optimize.root_scalar, does not handle cases where the crossing (which is set up as a root in this case) lies outside the interval provided¹.

This causes the crash here: to get the crossing for -2 sigma, the starting point (lower bound) of 0.1 already results in a CLs of 0.0500²:

DEBUG - cabinetry.fit - Signal_norm = 0.1000, expected -2 sigma CLs = 0.0500 (cached)

My guess is that this value is very slightly lower than 0.05 (just rounded in the log), so the algorithm would want to next probe a POI value below 0.1 (and find a CLs larger than 0.05 to confirm the root), but cannot because it is already at the lower bound.

Setting a different starting bracket via the bracket keyword argument should help solve this issue. You can use e.g. fit.limit(model, data, bracket=(0.01, 10))³. This is what the last ERRROR message from cabinetry is meant to communicate. When reading this in this context again myself, it might be confusing: the interval [0.1, 0.1] looks very unintuitive² and the message does not say that these values need to bracket the target POI value for the algorithm to succeed. We can try to make the message a bit more detailed. Suggestions are welcome!

Footnotes

1. We previously used another algorithm that could also handle such cases, but anecdotally it did not work very well: it took longer to converge when the root was within the interval, and frequently missed the root when moving too far outside of the interval (then running into numeric issues when the CLs becomes really small) and failed to converge. ↩

2. It probes this point where the POI is 0.1 twice, as it starts out by evaluating lower and upper bound as set from previous information (the observed limit was obtained, so some fits are already done and that info is used). In this case they coincide: the upper bound to start with is the lowest POI value for which CLs<0.05, the lower bound is the highest POI value for which CLs>0.05. If no such POI values exist for the lower bound (as in this case) it defaults to the default bracket setting, here 0.1. ↩ ↩²

3. The default lower bound of 0.1 might seem quite large, and it is: with typical POIs of O(1) it is not too unlikely that limits end up requiring using lower POI values. In practice I have found that setting a lower default bound of e.g. 0.01 can cause convergence issues when the limits are O(10). This is something that should be solvable: more clever algorithms can probably handle all this. If anyone reading this is interested in developing better approaches, please get in touch! ↩

[DaariiB]

It is working now. Thank you very much!

[nhartman94]

Thank you so much Alex for the quick response, with the simultaneously very detailed answer.
It's so nice this was such a quick fix :)
If there was any way for the error message to suggest the next bound, e.g, like you said to us just now:

" Try setting a different starting bracket via thefit.limit(model, data, bracket=(0.01, 10)) "

I think would have been helpful for me to understand how to solve the issue in case that's easy to add to the message.

Thank you again so much for your support!!!