32b vs 64b backend

[NiclasEich]

Hi all,

I was trying out different backends for runtime optimisation, when I noticed that my computed limits differ by a factor of 2, which goes beyond floating point precision for "normal" computations.

Below is an example, where I would not expect the "32b" backend to fail as it does.

The difference between numpy-64b backend and numpy-32b backend also occurs when I only work on np.float32 arrays, which should be save to convert to np.float64, so I suspect I either misunderstand what the backend does or there is something very sensitive to floating point precision.

There is also this GitHub-issue, aking something related

(myenv) user@machine: python3 -m IPython
Python 3.10.2 | packaged by conda-forge | (main, Feb  1 2022, 19:29:00) [GCC 9.4.0]
Type 'copyright', 'credits' or 'license' for more information
IPython 8.1.1 -- An enhanced Interactive Python. Type '?' for help.

In [1]: import numpy as np
   ...: import pyhf
   ...:
   ...: model = pyhf.simplemodels.uncorrelated_background(
   ...:     signal=[5.0, 10.0], bkg=[50.0, 60.0], bkg_uncertainty=[5.0, 12.0]
   ...: )

In [2]: CLs_obs, CLs_exp = pyhf.infer.hypotest(
   ...:     1.0,  # null hypothesis
   ...:     [53.0, 65.0] + model.config.auxdata,
   ...:     model,
   ...:     test_stat="qtilde",
   ...:     return_expected_set=True,
   ...: )
   ...: print(f"      Observed CLs: {CLs_obs:.4f}")
   ...: for expected_value, n_sigma in zip(CLs_exp, np.arange(-2, 3)):
   ...:     print(f"Expected CLs({n_sigma:2d} σ): {expected_value:.4f}")
   ...:
      Observed CLs: 0.4957
Expected CLs(-2 σ): 0.0832
Expected CLs(-1 σ): 0.1828
Expected CLs( 0 σ): 0.3705
Expected CLs( 1 σ): 0.6437
Expected CLs( 2 σ): 0.8854

In [3]: pyhf.set_backend("NumPy", precision="32B")

In [4]: CLs_obs, CLs_exp = pyhf.infer.hypotest(
   ...:     1.0,  # null hypothesis
   ...:     [53.0, 65.0] + model.config.auxdata,
   ...:     model,
   ...:     test_stat="qtilde",
   ...:     return_expected_set=True,
   ...: )
   ...: print(f"      Observed CLs: {CLs_obs:.4f}")
   ...: for expected_value, n_sigma in zip(CLs_exp, np.arange(-2, 3)):
   ...:     print(f"Expected CLs({n_sigma:2d} σ): {expected_value:.4f}")
   ...:
.../conda/envs/vhclustering/lib/python3.10/site-packages/pyhf/infer/calculators.py:369: RuntimeWarning: invalid value encountered in double_scalars
  teststat = (qmu - qmu_A) / (2 * self.sqrtqmuA_v)
      Observed CLs: nan
Expected CLs(-2 σ): 1.0000
Expected CLs(-1 σ): 1.0000
Expected CLs( 0 σ): 1.0000
Expected CLs( 1 σ): 1.0000
Expected CLs( 2 σ): 1.0000

[alexander-held]

The issue seems to appear for simple maximum likelihood fits as well, and doesn't exist with jax:

import pyhf

model = pyhf.simplemodels.uncorrelated_background(
    signal=[5.0], bkg=[50.0], bkg_uncertainty=[5.0]
)

for backend in ["numpy", "jax"]:
    for precision in ["64b", "32b"]:
        pyhf.set_backend(backend, precision=precision)
        print(
            f"{backend} {precision}: {pyhf.infer.mle.fit([53.0] + model.config.auxdata, model)}"
        )

This outputs

numpy 64b: [0.60001393 1.00000092]
numpy 32b: [1. 1.]
jax 64b: [0.60001705 1.00000086]
jax 32b: [0.600017  1.0000008]

with numpy v1.23.1 and current pyhf master.

[alexander-held]

The behavior seems minimizer-dependent too. With

pyhf.set_backend(backend, "minuit", precision=precision)

I see the following:

numpy 64b: [0.60072303 0.99995603]
numpy 32b: [0.61967224 0.99993825]  # surprisingly far off?
jax 64b: [0.60072306 0.99995579]
# jax 32b fit does not succeed

[NiclasEich]

It also appeared for my pyhf.infer.intervals.upperlimit (cannot provide a minimal example right now).

[kratsg]

So this is expected. We don't promise good support for 32b backends (part of this is just due to the fact that it's very hard to really keep everything perfectly consistent, but we've done our best for now). You can see this in the test we have here (

pyhf/tests/test_optim.py

Lines 48 to 94 in 8e242c9

	@pytest.mark.parametrize('do_grad', [False, True], ids=['no_grad', 'do_grad'])
	def test_minimize(tensorlib, optimizer, do_grad, do_stitch):
	pyhf.set_backend(tensorlib(precision="64b"), optimizer())
	m = pyhf.simplemodels.uncorrelated_background([50.0], [100.0], [10.0])
	data = pyhf.tensorlib.astensor([125.0] + m.config.auxdata)
	# numpy does not support grad
	if pyhf.tensorlib.name == 'numpy' and do_grad:
	with pytest.raises(pyhf.exceptions.Unsupported):
	pyhf.infer.mle.fit(data, m, do_grad=do_grad)
	else:
	identifier = f'{"do_grad" if do_grad else "no_grad"}-{pyhf.optimizer.name}-{pyhf.tensorlib.name}'
	expected = {
	# numpy does not do grad
	'do_grad-scipy-numpy': None,
	'do_grad-minuit-numpy': None,
	# no grad, scipy, 64b
	'no_grad-scipy-numpy': [0.49998815367220306, 0.9999696999038924],
	'no_grad-scipy-pytorch': [0.49998815367220306, 0.9999696999038924],
	'no_grad-scipy-tensorflow': [0.49998865164653106, 0.9999696533705097],
	'no_grad-scipy-jax': [0.4999880886490433, 0.9999696971774877],
	# do grad, scipy, 64b
	'do_grad-scipy-pytorch': [0.49998837853531425, 0.9999696648069287],
	'do_grad-scipy-tensorflow': [0.4999883785353142, 0.9999696648069278],
	'do_grad-scipy-jax': [0.49998837853531414, 0.9999696648069285],
	# no grad, minuit, 64b - quite consistent
	'no_grad-minuit-numpy': [0.5000493563629738, 1.0000043833598724],
	'no_grad-minuit-pytorch': [0.5000493563758468, 1.0000043833508256],
	'no_grad-minuit-tensorflow': [0.5000493563645547, 1.0000043833598657],
	'no_grad-minuit-jax': [0.5000493563528641, 1.0000043833614634],
	# do grad, minuit, 64b
	'do_grad-minuit-pytorch': [0.500049321728735, 1.00000441739846],
	'do_grad-minuit-tensorflow': [0.5000492930412292, 1.0000044107437134],
	'do_grad-minuit-jax': [0.500049321731032, 1.0000044174002167],
	}[identifier]
	result = pyhf.infer.mle.fit(data, m, do_grad=do_grad, do_stitch=do_stitch)
	rel_tol = 1e-5 if "no_grad" in identifier else 1e-6
	# Fluctuations beyond precision shouldn't matter
	abs_tol = 1e-8
	# check fitted parameters
	assert pytest.approx(
	expected, rel=rel_tol, abs=abs_tol
	) == pyhf.tensorlib.tolist(
	result
	), f"{identifier} = {pyhf.tensorlib.tolist(result)}"

) to check consistency across the minimization. In particular, we dropped 32b in #1503 / #1504. Even worse was that the 32b tests were significantly different across architectures... even between my mac and the Mac OS runners GitHub provided... and similarly for @matthewfeickert's linux machines vs GitHub's linux runners... It was just a headache and a mess and was pretty clear to us that there were a lot more things happening under the hood.

[alexander-held]

Would it be worth to show a warning when 32b is used? From the few very simple examples here this looks very dangerous to use in practice. Perhaps a mention in the documentation could also be useful.

Related to documentation: https://pyhf.readthedocs.io/en/v0.6.3/_generated/pyhf.set_backend.html says the precision is backend-dependent, as far as I can tell it's all 64b at the moment (but technically true that the default is set on a per-backend level). Were some of the backends 32b by default in the past?

[kratsg]

Yes, some of the backends were 32b by default. Can you open an issue with the stale documentation link so we can fix that line up? A warning would be good to add as well and we can certainly do that.

[alexander-held]

done: #1923

[NiclasEich]

This is of course not to me to decide, but I would rather disable 32b than having wrong results that are different by factors than 64b. The documentation is not clear about it and I can imagine that there are cases where people would not realise that something is going wrong.

[kratsg]

This is one of the reasons we switched to 64b by default across the board. We can add a warning in for 32b, but we would prefer not to disable because this is something the user should decide for themselves (and we want to be able to explore this more in the future).