TST, CI: `test_partial_fit` is failing in Linux CI (with NumPy `2.5.0`)

[tylerjereddy]

Here is a sample log of the test failure: https://github.com/lanl/GFDL/actions/runs/28083278183/job/83143018162?pr=116. The PR where that occurs is clearly unrelated, since it is just a dependabot bump of one of our actions versions (#116).

It is reproducible for me on x86_64 Linux locally via this incantation: python -m pytest "src/gfdl/tests/test_model.py::test_partial_fit[EnsembleGFDLClassifier-None-coeffs_-None-range-log_softmax-2-hidden_layer_sizes1]"

The error is only reproducible with NumPy 2.5.0, which was released 3 days ago (https://github.com/numpy/numpy/releases/tag/v2.5.0). Downgrading to NumPy 2.4.6 allows the test case to pass again.

A sample traceback from running locally on x86_64 Linux is shown below the fold. The relative tolerance violation is actually rather large at 1.46476215 -- I would not be comfortable adjusting the rtol to be that loose.

Details

>       assert_allclose(
            getattr(pf_model, attr), getattr(ff_model, attr), rtol=1e-5, atol=1e-3
            )
E       AssertionError: 
E       Not equal to tolerance rtol=1e-05, atol=0.001
E       
E       Mismatched elements: 1 / 100 (1%)
E       Mismatch at index:
E        [1, 3, 1]: 0.0003323436399264836 (ACTUAL), -0.0007150832738717083 (DESIRED)
E       Max absolute difference among violations: 0.00104743
E       Max relative difference among violations: 1.46476215
E        ACTUAL: array([[[-1.215602e-01, -1.612826e-01],
E               [-6.078257e-02, -8.063879e-02],
E               [-4.978902e-06,  4.978902e-06],...
E        DESIRED: array([[[-1.215602e-01, -1.612826e-01],
E               [-6.078257e-02, -8.063879e-02],
E               [-4.978902e-06,  4.978902e-06],...

src/gfdl/tests/test_model.py:581: AssertionError

Some things that might be helpful to do here:

• one could git bisect against NumPy to determine the root cause of the introduction of this test failure--is it an upstream bug?
• one could check if i.e., it is the OpenBLAS version associated with NumPy that is the issue, or NumPy source proper
• once the root cause is determined, propose a course of action--is it a bug on our side or upstream? is it an interaction with our library directly or with sklearn, etc.

[tylerjereddy]

The test error is not reproducible when NumPy 2.5.0 is built from source locally on x86_64 Linux, but it is reproducible from the pip installed wheel for 2.5.0, which may point a bit of suspicion toward the bundled OpenBLAS.

[tylerjereddy]

For the failing test case, the partial_fit coefficient value at the problem index reported ([1, 3, 1]) with NumPy 2.4.6 is 3.3265039510005064e-05 (vs. 3.323436E-04 with 2.5.0), while the expected full fit value looks pretty stable on 2.4.6 and 2.5.0, so the issue seems to be on the partial_fit side.

I spent some more time debugging this morning and discussed with the NumPy core team today, and they'd prefer if I bisect OpenBLAS and open an issue with them.

[tylerjereddy]

I reported this upstream at: OpenMathLib/OpenBLAS#5866.

This is the second time I've spent quite a bit of time dissecting out reproducers for upstream projects that affect us--would be great to have folks stepping up to help on these kinds of things. It is also a good learning experience.

[tylerjereddy]

@eiviani-lanl @nray assistance in (further) root causing the upstream matter and resolving it would be helpful here. This one is likely to be extremely challenging, but relevant to our long term sustainability and contributing back to the ecosystem we're built on is productive.

[tylerjereddy]

@eiviani-lanl @nray any progress here? Would be good to have some assistance beyond the project leader driving complex debugging tasks that involve upstream.

[nray]

@eiviani-lanl @nray any progress here? Would be good to have some assistance beyond the project leader driving complex debugging tasks that involve upstream.

I am working on this issue this week.

[tylerjereddy]

One thing that came up in the cross-linked sklearn PR is whether this issue should block that one, since the failing case here is for the ensemble estimator that isn't ported over there.

I'm a bit tired from travel, but poking around it certainly is possible to change the test parameters and cause platform-specific differences/failures on the "regular" estimators (probably because of Apple Accelerate vs. OpenBLAS backends), mostly via massaging alpha to lower values it seems. However, I'm not sure that's really the same problem, especially since that shouldn't transit through pinv anymore. Perhaps we could alleviate concerns to some degree by expanding the testing temporarily/a bit more exhaustively on the pinv codepath for regular/non-ensemble estimators (could use hypothesis library maybe) to confirm that we don't reproduce the problem for the non-ensemble estimator(s) on x86_64 Linux with NumPy 2.5.1 wheels.

That's still a bit of a pain to do convincingly, since we don't really know how rare the numerical issue is in practice, and if pinv is temporarily busted on Linux because of OpenBLAS one might presume that we can't really ever fully guarantee safety when we transit through.

[martin-frbg]

Are we even sure that it is pinv that is causing the problem ? I got the impression that the magic value from the "minimal reproducer" wasn't that helpful at all for predicting whether the full calculation would pass or fail, and I'd rather not have to sift through all BLAS calls that occur in the course of your CI job to see if there is significant divergence between BLAS kernels for "working" and "non-working" cpus, or if it is just creeping accumulation of FMA rounding errors driving an unstable algorithm over the edge

[tylerjereddy]

accumulation of FMA rounding errors driving an unstable algorithm over the edge

Yeah, agreed that this should be dissected carefully on our side. I asked my team to help me out a bit. Failing that I'll do another deep dive soon-ish.

[tylerjereddy]

For the failing test case on Linux the full fit array condition numbers before pinv vary from ...e+16 to ...e+17, and for partial fit it is even worse, hitting ...e+22.

The situation is pretty much the same on ARM Mac where the test case passes. Doesn't this mean that these arrays are basically singular? The small diff to verify that is available on my fork in branch treddy_more_debug_pinv.

Is it reasonable to expect that reconstructions from partial_fit and fit can match in these cases? Why/why not? In practice, pinv seems to mostly allow reconstruction in these cases, but apparently not always. I suppose someone could plot the trends in condition numbers and other related matrix/array properties to look for degrading trends with one backing OpenBLAS version vs. i.e., Accelerate, etc.