SciPy v1.16.0 gives different numerical results and Singular matrix E in LSQ subproblem

[matthewfeickert]

Summary

As found in PR #2592, scipy v1.16.0 gives different numerical results than scipy v1.15.3 for tests/test_calculator.py and for tests/test_optim.py it fails with

E           pyhf.exceptions.FailedMinimization: Singular matrix E in LSQ subproblem

OS / Environment

PRETTY_NAME="Ubuntu 22.04.5 LTS"
NAME="Ubuntu"
VERSION_ID="22.04"
VERSION="22.04.5 LTS (Jammy Jellyfish)"
VERSION_CODENAME=jammy
ID=ubuntu
ID_LIKE=debian
HOME_URL="https://www.ubuntu.com/"
SUPPORT_URL="https://help.ubuntu.com/"
BUG_REPORT_URL="https://bugs.launchpad.net/ubuntu/"
PRIVACY_POLICY_URL="https://www.ubuntu.com/legal/terms-and-policies/privacy-policy"
UBUNTU_CODENAME=jammy

Steps to Reproduce

$ uv pip install --upgrade -e '.[develop]'
$ uv pip install --upgrade scipy
$ pytest --ignore tests/benchmarks --ignore tests/test_notebooks.py tests/

=========================================================================================== test session starts ============================================================================================
platform linux -- Python 3.12.11, pytest-8.4.1, pluggy-1.6.0
Matplotlib: 3.10.3
Freetype: 2.6.1
benchmark: 5.1.0 (defaults: timer=time.perf_counter disable_gc=False min_rounds=5 min_time=0.000005 max_time=1.0 calibration_precision=10 warmup=False warmup_iterations=100000)
rootdir: /home/feickert/Code/GitHub/scikit-hep/pyhf
configfile: pyproject.toml
plugins: mock-3.14.1, requests-mock-1.12.1, mpl-0.17.0, console-scripts-1.4.1, socket-0.7.0, anyio-4.9.0, benchmark-5.1.0
collected 1518 items                                                                                                                                                                                       

tests/contrib/test_contrib_utils.py .....                                                                                                                                                            [  0%]
tests/contrib/test_viz.py .........                                                                                                                                                                  [  0%]
tests/test_backend_consistency.py ....                                                                                                                                                               [  1%]
tests/test_backends.py ............                                                                                                                                                                  [  1%]
tests/test_calculator.py ...FF.                                                                                                                                                                      [  2%]
tests/test_cli.py ..                                                                                                                                                                                 [  2%]
tests/test_combined_modifiers.py ................................................                                                                                                                    [  5%]
tests/test_compat.py ...                                                                                                                                                                             [  5%]
tests/test_constraints.py ............                                                                                                                                                               [  6%]
tests/test_custom_mods.py ..                                                                                                                                                                         [  6%]
tests/test_events.py .......                                                                                                                                                                         [  7%]
tests/test_examples.py .                                                                                                                                                                             [  7%]
tests/test_export.py ...........................                                                                                                                                                     [  9%]
tests/test_import.py ......................                                                                                                                                                          [ 10%]
tests/test_infer.py ...........................................................................................                                                                                      [ 16%]
tests/test_init.py ......                                                                                                                                                                            [ 16%]
tests/test_interpolate.py .........................................................................................................                                                                  [ 23%]
tests/test_jit.py ................................                                                                                                                                                   [ 25%]
tests/test_mixins.py ...                                                                                                                                                                             [ 26%]
tests/test_modifiers.py ............                                                                                                                                                                 [ 26%]
tests/test_optim.py .sss.s.......................................................................................FFFFF.FFFFF.sssss.sssss...................                                          [ 35%]
tests/test_paramsets.py .........                                                                                                                                                                    [ 36%]
tests/test_paramviewer.py .............                                                                                                                                                              [ 37%]
tests/test_patchset.py ...........................                                                                                                                                                   [ 39%]
tests/test_pdf.py ............................sssss.........................sssss.......sssss............FFFFF.....                                                                                  [ 45%]
tests/test_probability.py .........................                                                                                                                                                  [ 47%]
tests/test_public_api.py ...................................................................                                                                                                         [ 51%]
tests/test_public_api_repr.py ...............................                                                                                                                                        [ 53%]
tests/test_regression.py FFF.F.                                                                                                                                                                      [ 53%]
tests/test_schema.py .........................................................................                                                                                                       [ 58%]
tests/test_scripts.py ..............................................................                                                                                                                 [ 62%]
tests/test_simplemodels.py .............xxxx..xxxx.                                                                                                                                                  [ 64%]
tests/test_tensor.py .................................................................................................xx...s..sss.............xx.............................................sss.ss. [ 75%]
...................................                                                                                                                                                                  [ 78%]
tests/test_tensorviewer.py ......                                                                                                                                                                    [ 78%]
tests/test_teststats.py ..F.F..............                                                                                                                                                          [ 79%]
tests/test_toys.py ......                                                                                                                                                                            [ 80%]
tests/test_utils.py ..................                                                                                                                                                               [ 81%]
tests/test_validation.py .....F..F..F................                                                                                                                                                [ 83%]
tests/test_workspace.py ............................................................................................................................................................................ [ 94%]
.................................................................................                                                                                                                    [100%]

================================================================================================= FAILURES =================================================================================================
____________________________________________________________________________ test_asymptotic_calculator_has_fitted_pars[qtilde] ____________________________________________________________________________

test_stat = 'qtilde'

    @pytest.mark.parametrize('test_stat', ['qtilde', 'q', 'q0'])
    def test_asymptotic_calculator_has_fitted_pars(test_stat):
        model = pyhf.simplemodels.uncorrelated_background([1], [1], [1])
        data = [2, 1]  # [main, aux]
    
        calc = pyhf.infer.calculators.AsymptoticCalculator(data, model, test_stat=test_stat)
        calc.teststatistic(0 if test_stat == 'q0' else 1)
    
        assert hasattr(calc, 'fitted_pars')
        fitted_pars = calc.fitted_pars
        assert hasattr(fitted_pars, 'asimov_pars')
        assert hasattr(fitted_pars, 'fixed_poi_fit_to_data')
        assert hasattr(fitted_pars, 'fixed_poi_fit_to_asimov')
        assert hasattr(fitted_pars, 'free_fit_to_data')
        assert hasattr(fitted_pars, 'free_fit_to_asimov')
    
        rtol = 1e-5
        if test_stat == 'q0':
            assert pytest.approx([1.0, 1.0], rel=rtol) == pyhf.tensorlib.tolist(
                fitted_pars.asimov_pars
            )
            assert pytest.approx([0.0, 1.5], rel=rtol) == pyhf.tensorlib.tolist(
                fitted_pars.fixed_poi_fit_to_data
            )
            assert pytest.approx([0.0, 1.5], rel=rtol) == pyhf.tensorlib.tolist(
                fitted_pars.fixed_poi_fit_to_asimov
            )
            assert pytest.approx([1.0, 1.0], rel=rtol) == pyhf.tensorlib.tolist(
                fitted_pars.free_fit_to_data
            )
            assert pytest.approx([1.0, 1.0], rel=rtol) == pyhf.tensorlib.tolist(
                fitted_pars.free_fit_to_asimov
            )
        else:
            assert pytest.approx([0.0, 1.5], rel=rtol) == pyhf.tensorlib.tolist(
                fitted_pars.asimov_pars
            )
            assert pytest.approx([1.0, 1.0], rel=rtol) == pyhf.tensorlib.tolist(
                fitted_pars.fixed_poi_fit_to_data
            )
            assert pytest.approx([1.0, 1.1513553], rel=rtol) == pyhf.tensorlib.tolist(
                fitted_pars.fixed_poi_fit_to_asimov
            )
            assert pytest.approx([1.0, 1.0], rel=rtol) == pyhf.tensorlib.tolist(
                fitted_pars.free_fit_to_data
            )
            # lower tolerance for amd64 and arm64 to agree
            # FIXME: SciPy v1.16.0 gives a different result from SciPy v1.15.3
>           assert pytest.approx(
                [7.6470499e-05, 1.4997178], rel=1e-3
            ) == pyhf.tensorlib.tolist(fitted_pars.free_fit_to_asimov)
E           assert approx([7.647...± 0.00149972]) == [7.6321607829...7178921881127]
E             
E             comparison failed. Mismatched elements: 1 / 2:
E             Max absolute difference: 1.488911703012493e-07
E             Max relative difference: 0.0019508390157801656
E             Index | Obtained              | Expected               
E             0     | 7.632160782969875e-05 | 7.6470499e-05 ± 7.6e-08

calc       = <pyhf.infer.calculators.AsymptoticCalculator object at 0x7c782ca5e720>
data       = [2, 1]
fitted_pars = HypoTestFitResults(asimov_pars=array([0.        , 1.50000576]), free_fit_to_data=array([1., 1.]), free_fit_to_asimov=a...e-05, 1.49971789e+00]), fixed_poi_fit_to_data=array([1., 1.]), fixed_poi_fit_to_asimov=array([1.        , 1.15135536]))
model      = <pyhf.pdf.Model object at 0x7c76fc3c6900>
rtol       = 1e-05
test_stat  = 'qtilde'

tests/test_calculator.py:89: AssertionError
-------------------------------------------------------------------------------------------- Captured log call ---------------------------------------------------------------------------------------------
INFO     pyhf.pdf:pdf.py:782 Validating spec against schema: model.json
INFO     pyhf.pdf:pdf.py:478 adding modifier mu (1 new nuisance parameters)
INFO     pyhf.pdf:pdf.py:478 adding modifier uncorr_bkguncrt (1 new nuisance parameters)

File Upload (optional)

No response

Expected Results

For the tests to run with the same numerical results within accepted floating point precision.

Actual Results

>           raise exceptions.FailedMinimization(result)
E           pyhf.exceptions.FailedMinimization: Singular matrix E in LSQ subproblem

pyhf Version

pyhf, version 0.7.1.dev308

Code of Conduct

• I agree to follow the Code of Conduct

[matthewfeickert]

@MoAly98 I know that you're going to be busy over the next week, but would you have time to look into this and see if you can understand what differences exist between the two SciPy versions and why?

[matthewfeickert]

Probably relevant is the SciPy v1.16.0 release notes on scipy.optimize improvements.

If I search SciPy for "Singular matrix E in LSQ subproblem" it matches in scipy/optimize/_slsqp_py.py, and it mentions in the release notes that

The SLSQP Fortran 77 code was ported to C. When this method is used now the constraint multipliers are exposed to the user through the multiplier keyword of the returned OptimizeResult object.

and for pyhf/optimize/opt_scipy.py

pyhf/src/pyhf/optimize/opt_scipy.py

Line 62 in c30e337

* method (:obj:`str`): Minimization routine. Default is ``'SLSQP'``.

[matthewfeickert]

Asked the SciPy team on the Scientific Python forums: Minimization result success change after SLSQP Fortran to C rewrite

[matthewfeickert]

Full output (test_failure.txt) of

pytest --ignore tests/benchmarks --ignore tests/test_notebooks.py tests/ &> test_failure.txt

working from commit c30e337.

[ilayn]

That's me who did the translation so apologies if it caused any trouble. Could you direct me to the problem that I can test things out why the difference might be happening?

The main difference between the old Fortran version and the new C code is the inclusion of BLAS/LAPACK calls to the underlying library instead of using the very old hand-baked code that came with SLSQP. So that would be the first point of contact to check, otherwise it might be my error in translation somewhere that we missed.

However we are seeing a very strange error in Ubuntu platform as reported in scipy/scipy#23161 and we are not able to reproduce it anywhere else. So hopefully it is not that.

[MoAly98]

I've been trying to get a pure scipy reproducer but not there yet. Here is a standalone example that fails in 1.16 but not 1.15.3

import pyhf
import scipy

print(f"pyhf version: {pyhf.__version__}")
print(f"scipy version: {scipy.__version__}")

pyhf.set_backend("numpy", "scipy")

# Test data from test_optim fixture
source = {
    "bindata": {
        "data": [120.0, 180.0],
        "bkg": [100.0, 150.0],
        "bkgsys_up": [102, 190],
        "bkgsys_dn": [98, 100],
        "sig": [30.0, 95.0],
    },
}

spec = {
    "channels": [
        {
            "name": "singlechannel",
            "samples": [
                {
                    "name": "signal",
                    "data": source["bindata"]["sig"],
                    "modifiers": [{"name": "mu", "type": "normfactor", "data": None}],
                },
                {
                    "name": "background",
                    "data": source["bindata"]["bkg"],
                    "modifiers": [
                        {
                            "name": "bkg_norm",
                            "type": "histosys",
                            "data": {
                                "lo_data": source["bindata"]["bkgsys_dn"],
                                "hi_data": source["bindata"]["bkgsys_up"],
                            },
                        }
                    ],
                },
            ],
        }
    ]
}

print("Creating model...")
pdf = pyhf.Model(spec, poi_name="mu")
data = source["bindata"]["data"] + pdf.config.auxdata

init_pars = pdf.config.suggested_init()
par_bounds = pdf.config.suggested_bounds()

print(f"Init pars: {init_pars}")
print(f"Par bounds: {par_bounds}")
print(f"Data: {data}")
print(f"Auxdata: {pdf.config.auxdata}")
print(f"POI index: {pdf.config.poi_index}")
print(f"Par names: {pdf.config.par_names}")

try:
    result = pyhf.optimizer.minimize(
        pyhf.infer.mle.twice_nll, data, pdf, init_pars, par_bounds
    )
    print(f"SUCCESS: {pyhf.tensorlib.tolist(result)}")
except pyhf.exceptions.FailedMinimization as e:
    print(f"FAILED: {e}")
    print(f"Result: {e.result}")

[matthewfeickert]

That's me who did the translation so apologies if it caused any trouble.

Thanks very much for the work that you did on this port, @ilayn! No small feat!

Could you direct me to the problem that I can test things out why the difference might be happening?

@MoAly98 has a standalone reproducer in #2593 (comment), but as that is using our pyhf models that is probably still pretty opaque to anyone expecting just NumPy and SciPy. We can work with you to try to get something that is even more simple, but that might take a bit to isolate what the issue is, given that there aren't failures for every use of SLSQP given that the following works for both SciPy versions

import pyhf


pyhf.set_backend("numpy", "scipy")

model = pyhf.simplemodels.uncorrelated_background(
    signal=[12.0, 11.0], bkg=[50.0, 52.0], bkg_uncertainty=[3.0, 7.0]
)
data = [51, 48] + model.config.auxdata

init_pars = model.config.suggested_init()
par_bounds = model.config.suggested_bounds()

result = pyhf.optimizer.minimize(
    pyhf.infer.mle.twice_nll, data, model, init_pars, par_bounds
)
print(result)

However we are seeing a very strange error in Ubuntu platform as reported in scipy/scipy#23161 and we are not able to reproduce it anywhere else. So hopefully it is not that.

We're seeing this happen across multiple OS, so there's nothing to indicate yet that this is related to the Ubuntu specific issue in that problem.

[ilayn]

OK then at least we know what the problem would be since the error is coming from a well-defined place (famous last words)

https://github.com/scipy/scipy/blob/51dfbcc7b38b251b274cdeb1d6e028a83f3f0cb8/scipy/optimize/__slsqp.c#L815-L823

Here I can add some print statements and see what the difference is. Probably some machine specific detail I might have missed.

[ilayn]

A preliminary run shows that

pyhf version: 0.7.6
scipy version: 1.17.0.dev0+git20251013.51dfbcc
Creating model...
Init pars: [0.0, 1.0]
Par bounds: [(-5.0, 5.0), (0, 10)]
Data: [120.0, 180.0, 0.0]
Auxdata: [0.0]
POI index: 1
Par names: ['bkg_norm', 'mu']
[...]
R diag 0 : 1.0000000000000000e+00
R diag 1 : 1.0000000000000000e+00
R diag 0 :                   -nan

Something is NaN contaminated which can happen since Fortran and C treat them a bit differently. This is treated as a singular matrix and then you get the error.

[MoAly98]

I can build a scipy example which runs with 1.16 but not 1.15. There seems to be a point in the minimization at which one of the parameters makes a large jump (4th evaluation in the example below) to parameter values with a nan likelihood. In 1.15.3 it seems like the optimizer recovers, but in 1.16.0 it just gets stuck.

import numpy as np
from scipy.optimize import minimize
import scipy

# Configuration constants
OBSERVED_DATA_MAIN = np.array([120.0, 180.0])
OBSERVED_DATA_AUX = np.array([0.0])
BACKGROUND_NOMINAL = np.array([100.0, 150.0])
SIGNAL_NOMINAL = np.array([30.0, 95.0])

# Histosys modifier bounds for bkg_norm
HISTOSYS_LO = np.array([98.0, 100.0])
HISTOSYS_HI = np.array([102.0, 190.0])

# Optimization settings
INITIAL_PARAMETERS = np.array([0.0, 1.0])  # [bkg_norm, mu]
PARAMETER_BOUNDS = [(-5.0, 5.0), (0.0, 10.0)]  # [(bkg_norm_min, bkg_norm_max), (mu_min, mu_max)]

# Test grid settings
TEST_GRID_STEP = 0.5

print(f"scipy version: {scipy.__version__}\n")

# =============================================================================
# Core NLL Functions
# =============================================================================

# The expected counts as a function of [bkg_norm, mu]:
def expected_main(bkg_norm, mu):
    """Expected counts in main bins"""
    # histosys interpolation
    if bkg_norm >= 0:
        # interpolate with hi
        bkg_mod = BACKGROUND_NOMINAL + bkg_norm * (HISTOSYS_HI - BACKGROUND_NOMINAL)
    else:
        # interpolate with lo
        bkg_mod = BACKGROUND_NOMINAL + bkg_norm * (BACKGROUND_NOMINAL - HISTOSYS_LO)

    return bkg_mod + mu * SIGNAL_NOMINAL

def expected_aux(bkg_norm):
    """Expected auxdata (constraint term)"""
    # Gaussian constraint on bkg_norm centered at 0 with width 1
    return OBSERVED_DATA_AUX

def nll(pars):
    """Negative log-likelihood"""
    bkg_norm, mu = pars

    exp_main = expected_main(bkg_norm, mu)
    exp_aux = expected_aux(bkg_norm)

    # Poisson likelihood for main bins
    # -log L = sum(expected - observed * log(expected))
    # Need to handle negative expected counts
    if np.any(exp_main <= 0):
        return np.nan

    nll_main = np.sum(exp_main - OBSERVED_DATA_MAIN * np.log(exp_main))

    # Gaussian constraint for auxdata
    # -log L = 0.5 * (obs - exp)^2 / sigma^2, sigma=1 (?)
    nll_aux = 0.5 * np.sum((OBSERVED_DATA_AUX - exp_aux)**2)

    return nll_main + nll_aux

# =============================================================================
# Test SLSQP
# =============================================================================
print("\n" + "=" * 70)
print("SLSQP optimization:")
print("=" * 70)

eval_count = [0]
def nll_traced(pars):
    eval_count[0] += 1
    result = nll(pars)
    exp_main = expected_main(pars[0], pars[1])
    print(f"  Eval {eval_count[0]:2d}: pars={pars}, exp_main={exp_main}, nll={result}")
    return result

result = minimize(
    nll_traced,
    INITIAL_PARAMETERS,
    method='SLSQP',
    bounds=PARAMETER_BOUNDS,
    constraints=None,
    options={'disp': True, 'maxiter': 100},
)

print(f"\nResult:")
print(f"  Success: {result.success}")
print(f"  Message: {result.message}")
print(f"  Final x: {result.x}")
print(f"  Final nll: {result.fun}")

# =============================================================================
# Test with another minimizer that works in both versions
# =============================================================================
print("\n" + "=" * 70)
print("L-BFGS-B optimization:")
print("=" * 70)

eval_count = [0]
def nll_traced(pars):
    eval_count[0] += 1
    result = nll(pars)
    exp_main = expected_main(pars[0], pars[1])
    print(f"  Eval {eval_count[0]:2d}: pars={pars}, exp_main={exp_main}, nll={result}")
    return result

result = minimize(
    nll_traced,
    INITIAL_PARAMETERS,
    method='trust-constr',
    bounds=PARAMETER_BOUNDS,
    constraints=None,
    options={'disp': True, 'maxiter': 100},
)

print(f"\nResult:")
print(f"  Success: {result.success}")
print(f"  Message: {result.message}")
print(f"  Final x: {result.x}")
print(f"  Final nll: {result.fun}")