Bug in scipy minimization: init_pars and fixed_vals are not synced

[kratsg]

Description

Related: #1051 (found as part of this PR).

There is a bug in pyhf's codebase (since ~forever) where the init_par for a fixed parameter might differ from the constrained value set for that fixed parameter. See for example v0.4.4:

pyhf/src/pyhf/infer/mle.py

Lines 52 to 90 in 6ac0f62

	def fixed_poi_fit(poi_val, data, pdf, init_pars=None, par_bounds=None, **kwargs):
	"""
	Run a maximum likelihood fit with the POI value fixed.
	Example:
	>>> import pyhf
	>>> pyhf.set_backend("numpy")
	>>> model = pyhf.simplemodels.hepdata_like(
	... signal_data=[12.0, 11.0], bkg_data=[50.0, 52.0], bkg_uncerts=[3.0, 7.0]
	... )
	>>> observations = [51, 48]
	>>> data = pyhf.tensorlib.astensor(observations + model.config.auxdata)
	>>> test_poi = 1.0
	>>> pyhf.infer.mle.fixed_poi_fit(test_poi, data, model, return_fitted_val=True)
	(array([1. , 0.97224597, 0.87553894]), array([28.92218013]))
	Args:
	data: The data
	pdf (~pyhf.pdf.Model): The statistical model adhering to the schema model.json
	init_pars (`list`): Values to initialize the model parameters at for the fit
	par_bounds (`list` of `list`\s or `tuple`\s): The extrema of values the model parameters are allowed to reach in the fit
	kwargs: Keyword arguments passed through to the optimizer API
	Returns:
	See optimizer API
	"""
	_, opt = get_backend()
	init_pars = init_pars or pdf.config.suggested_init()
	par_bounds = par_bounds or pdf.config.suggested_bounds()
	return opt.minimize(
	twice_nll,
	data,
	pdf,
	init_pars,
	par_bounds,
	[(pdf.config.poi_index, poi_val)],
	**kwargs,
	)

pyhf/src/pyhf/optimize/opt_scipy.py

Lines 16 to 53 in 6ac0f62

	def minimize(
	self,
	objective,
	data,
	pdf,
	init_pars,
	par_bounds,
	fixed_vals=None,
	return_fitted_val=False,
	):
	"""
	Find Function Parameters that minimize the Objective.
	Returns:
	bestfit parameters
	"""
	fixed_vals = fixed_vals or []
	indices = [i for i, _ in fixed_vals]
	values = [v for _, v in fixed_vals]
	constraints = [{'type': 'eq', 'fun': lambda v: v[indices] - values}]
	result = minimize(
	objective,
	init_pars,
	constraints=constraints,
	method='SLSQP',
	args=(data, pdf),
	bounds=par_bounds,
	options=dict(maxiter=self.maxiter),
	)
	try:
	assert result.success
	except AssertionError:
	log.error(result)
	raise
	if return_fitted_val:
	return result.x, result.fun
	return result.x

Expected Behavior

I would've thought the underlying optimizers (particularly scipy) would be smart enough to handle or recognize when the constraint/init_pars are at odds with each other.

Actual Behavior

See description.

Steps to Reproduce

Follow along with me:

call fixed_poi_fit with:

• poi_val = 0.5
• data
• pdf
• init_pars=[1, 1, 1]

This will call opt.minimize with:

• objective=twice_nll
• data=data
• pdf=pdf
• init_pars=[1, 1, 1]
• fixed_vals=[(0, 0.5)]

at this point, opt.minimize (in the case of scipy) will do

        constraints = [{'type': 'eq', 'fun': lambda v: v[indices] - values}]

which becomes essentially

        constraints = [{'type': 'eq', 'fun': lambda v: v[0] - 0.5}]

but still calls scipy.optimize.minimize with init_pars=[1, 1, 1]. If one changes this to [0.5, 1, 1] -- one observes a slightly different result... but probably more correct.

Minuit does not have the same bug here because of

pyhf/src/pyhf/optimize/opt_minuit.py

Line 45 in 6ac0f62

initvals['p{}'.format(index)] = value

which updates the initval with the fixed value initvals['p{}'.format(index)] = value.

I was worried that the major refactor I did #951 broke this or changed this. In fact, I managed to keep the same bug in scipy, and did not introduce this bug into minuit.