BUG: Moore-Penrose Pseudoinverse numerical stability/consistency affected by recent OpenBLAS activity

@seberg

I noticed a test failure downstream where a machine learning estimator was performing a partial_fit that leverages the Moore-Penrose pseudoinverse for a closed form solution (lanl/GFDL#117). Although only 1/100 elements mismatched, the relative tolerance diff was quite substantial on that element (1.46476215) vs. previous NumPy/OpenBLAS combinations. This failure showed up starting with NumPy 2.5.0 from 3 days ago, which shipped with a newer version of OpenBLAS.

I then dissected out a plain text, but rather lengthy, reproducer on the NumPy level. I've pushed it up to GitHub at https://github.com/tylerjereddy/pinv_issue, and it prints out a single element of np.pinv output from a fixed 2D input. The input is a (25, 25) array carefully dissected out of our in house testsuite--the size of that array in plain text is the main reason I've provided it in a separate repo--it is perhaps too large to dump here and so you can confirm it is plain text, etc.

Below are the results it prints in a few different scenarios -- the main problem with consistency appears to be with the OpenBLAS that ships with the latest x86_64 NumPy Linux wheels, since building NumPy 2.5.0 from source with an older OpenBLAS passes our testing just fine. I realize that these are just "proxy indicators" for the problem rather than the full picture, especially since different values on different platforms and with different linalg backends are acceptable, but nonetheless, the outcome on NumPy 2.5.0 with OpenBLAS from wheel is not acceptable in the broader testing, so I tried to cut this down to a proxy reproducer. Just to emphasize here--it is "ok" that some of the values in the table below are different, but the 0.3609... value with NumPy 2.5.0 + the OpenBLAS that ships with that wheel, is not ok (though this may only be part of the total numerical instability story). It is likely that the input array is pathological in a variety of ways, and this is not uncommon for batched input over a design matrix. Nonetheless, most other platforms and linalg backends and older OpenBLAS versions are doing just fine here in the broader numerical testing we do, and I wanted to keep the proxy test focused on a single float value for now.

NumPy version	Result
NumPy `2.5.0` wheel on x86_64 Linux	`0.360954047299216`
NumPy `2.5.0` built from source on x86_64 Linux with OpenBLAS `0.3.20`	`0.5567409444424644` (passes our broader numerical testing)
NumPy `2.4.6` wheel on x86_64 Linux	`0.08970465208832239`
NumPy `2.3.4` wheel on x86_64 Linux	`0.08970465208832239`
NumPy `2.5.0` wheel on M series Mac	`0.1422214449163197`
NumPy `2.4.6` wheel on M series Mac	`0.1422214449163197`

I then talked to @seberg and @ngoldbaum and they suggested I bisect on OpenBLAS and open an issue here, so I bisected OpenBLAS alongside NumPy 2.5.0 to find the first point where x86_64 Linux + OpenBLAS produces the problem value above with NumPy 2.5.0 (the 0.3609... value) from wheel install (see below) when linked to NumPy. In fact, I ended up needing to use LD_PRELOAD to get the right behavior (force NumPy to link the correct source-built OpenBLAS, pkg config env vars/settings were not sufficient on their own...). I realize that this isn't a reproducer at the C level for OpenBLAS, and that these single proxy floats are just part of the story, but thought I'd get this started. I can work toward that if you really need it.

To further complicate matters, some of the commits in the bisection had to be skipped because they failed to build; presumably those commits really are broken for other reasons. So, the best I could do was this output that narrows it down to two possible problem commits

There are only 'skip'ped commits left to test.
The first bad commit could be any of:
f5f789fc529ffcc4fb7dd5f525fee1e719339a3a
d9bb8f369f5d80cb0b074f01bde156a6dda1f193
We cannot bisect more!

Those both reference Reference-LAPACK PRs. I did confirm that the latest develop commit (a36e22c) is still broken (produces 0.360954047299216 using the reproducer above on x86_64 Linux).

Sorry if this is a little tricky to follow, but it is the result of several hours of careful dissection and I did not use AI--this is a real problem affecting our estimator research at DOE, and we're happy to help out if we can, but would probably need a bit more guidance. I wonder if an expert can look at the suspicious commits with a knowledge of the ops used by pinv and get a pretty good sense for where the issue is.

I notice both of the two "possibly bad" commits belong to PR #5721 a.k.a Reference-LAPACK 1196+1019+1020 implementing the new LARF1F and LARF1L functions. It so happens that a bug in the handling of negative strides by these functions was found and fixed just this week, I have queued #5865 to copy this correction into OpenBLAS. I wonder if you could give this PR a try - otherwise I'll probably get around to merging it tomorrow and hopefully find time to run your numpy reproducer.

I didn't see an improvement at hash 562aeed from that PR, the result was still 0.360954047299216.

To check, in the same environment, when I checkout OpenBLAS v0.3.31 and clean rebuild, the LD_PRELOAD of that version gives me 0.08970465208832239, equivalent to older passing versions above.

Indeed flags d9bb8f3 as the first bad commit. I still need to check if this is a genuine bug in the upstream change or if I missed part of it. (I also still need to find out which of the changed files is on the call graph of whatever linalg.pinv does internally...)

Call graph is fairly straightforward, DGESDD calling DLARF1F via one of DGEBD2, DORM2R and DORML2 - I see no missed upstream changes there. Possibly it is a question of algorithmic stability against small numerical changes from the different BLAS codepaths taken - I notice the supposed canary value varies between 0.36 for "Haswell" and 0.10 for "Generic" OpenBLAS TARGET setting....

A numpy 2.5.0 wheel linked against Reference-LAPACK 3.11 (i.e. before the introduction of DLARF1F) returns 0.21898014163337262 - I really wish we had a better metric for what constitutes a good or bad build. Swapping in 3.12.1 gives 0.09853545542524698

There is one intervening LAPACK PR that I missed - their 1163 "add non-trivial terminating case to ?LARFT" (to be added tomorrow), but although it is on the call graph, fixing it does not change the 0.3609...

I don't think I missed any relevant patches - what the DLARF1F PR mostly did is add a dedicated DLARF variant that replaces the first matrix element with a 1 temporarily, instead of relying on the caller to do this.
It also rearranged some calls to DAXPY and DGEMV, and the latter appears to have "resonated" with the elusive inaccuracy in the DGEMV AVX microkernel for x86_64 already claimed in #4324. I can make OpenBLAS return 0.07586 by just removing the #if defined... include dgemv_t_microk_haswell-4.c at the top of kernel/x86_64/dgemv_t_4.c

If we need a better grasp of the numerical situation beyond the proxy float (which is kind of terrible as noted above), I suppose we could try to i.e., plot a heatmap of the values produced by pinv in the various scenarios to see if there's some kind of numerical trend in the data across the various pinv returns values for that input (that is preserved in the "better behaved" cases).

I guess so - but perhaps it would be sufficient for now to see if disabling the dgemv_t microkernel solves the problem in your actual code, rather than just producing a more likable value for the proxy float in the reproducer ? As said, nothing wrong has been found with it so far beyond the usual misgivings about FMA operations affecting rounding, but if there happens to be any step in your code that is particularly susceptible to small numerical changes...

So I installed your GFDL-0.2.0 wheel and ran python -m pytest --pyargs gfdl with variations of libopenblas in numpy.libs, but disabling the microkernels for transpose or both forms of gemv did not appear to help

E       AssertionError: 
E       Not equal to tolerance rtol=1e-07, atol=1e-05
E       
E       Mismatched elements: 1 / 1 (100%)
E       Max absolute difference among violations: 1.79641559e-05
E       Max relative difference among violations: 1.88590432e-05
E        ACTUAL: array(0.952531)
E        DESIRED: array(0.952549)

lib64/python3.13/site-packages/gfdl/tests/test_model.py:516: AssertionError

FAILED lib64/python3.13/site-packages/gfdl/tests/test_model.py::test_rtol_ensemble[5.0-1e-15-0.7222222222222222-0.9525486362311113] - AssertionError:
================================================== 1 failed, 6902 passed, 3 skipped, 2 xfailed, 9 warnings in 255.73s (0:04:15) ===================================================

`

did not appear to help

It was only test_partial_fit that was failing--so if your result from #5866 (comment) implies that the things you disabled allowed test_partial_fit to pass (was it failing for you before?), then it actually is helping? The test failure you cite looks less serious at first glance, the rtol violation there is a few orders of magnitude closer vs. the failure we have for test_partial_fit.

Oh right, I misread - and somehow did not anticipate additional failures. Let me retest - but indeed the unmodified develop branch generates two failures (the aforementioned one in test_rtol_ensemble and the one in test_partial_fit this issue is about) while deactivating the dgemv_t_microk_haswell.c makes test_partial_fit pass

Great -- I suppose this doesn't necessarily convince either of us that anything other than "that's floating point math" is happening here, but at least that's some traction if there is a real problem.

Yes - some code appears to be very sensitive to numerical error, but it may be too early to say whose code it is. (Of note, I can also make test_partial pass if I replace the DAXPY microkernel with a trivial C loop. This was not apparent with the pinv reproducer. On the other hand, all tests pass if I build OpenBLAS for AVX512 SKYLAKEX or newer instead of the AVX2 "Haswell", without any changes to the current develop branch - although the pinv reproducer still returns the "bad" 0.3609...)

ENH: add heatmap plot code

More detailed analysis at lanl/GFDL#117 (comment) indicates that the input arrays in the problem scenario are guaranteed to be Hermitian, and generally have extraordinarily high condition numbers.

I suspect the analysis over there narrows us down to 3-4 possible OpenBLAS routines that have changed if we pursue the non-Hermitian SVD path with extremely high condition number despite having a Hermitian matrix.

From a pragmatic standpoint, it seems we can circumvent the problem downstream by guiding the algorithm with hermitian=True, which may be a best practice in this scenario anyway. As for whether the problem/naive path should provide better stability across OpenBLAS releases, I'm still not sure on that one, but all evidence still points to a change in OpenBLAS behavior in svd with very high condition numbers.

Thank you for the careful investigation here. We encountered what appears to be a closely related downstream numerical regression on x86_64 Linux.

This investigation was AI-assisted, with the observations below manually reproduced and verified.

In our case, NumPy 2.4.6 and NumPy 2.5.x produce slightly different results for repeated decompositions of small, symmetric positive-definite 3x3 matrices. The first differences are in the last bits of both:

np.linalg.svd (DGESDD)
np.linalg.eigh (DSYEVD)

The input matrices immediately before those calls are bit-for-bit identical. The differences are very small (roughly 1e-20 in the projected matrices and 1e-16 in a subsequent square root), but they enter an iterative optimization and eventually cross a discrete rounding boundary much later in the calculation.

A few observations that may be useful:

NumPy 2.4.6 reproduces the previous result; NumPy 2.5.0 already shows the changed result, as does 2.5.2.
Running NumPy 2.5.2 while preloading the OpenBLAS library bundled with the NumPy 2.4.6 wheel restores the previous decomposition results and makes the full downstream regression pass.
Keeping NumPy 2.5.2 but routing both decompositions through SciPy with explicit gesdd and evd LAPACK drivers also restores the previous downstream result. Replacing only the SVD step was not sufficient.

I appreciate that last-bit variation across BLAS/LAPACK implementations and CPU paths is generally expected, and our iterative/discrete calculation is unusually sensitive to it. I am sharing this mainly because the NumPy-version boundary and the LD_PRELOAD result closely resemble the observations already reported here.

Can you share a small reproducer (and/or check if the OpenBLAS workaround outlined in #5866 (comment) "solves" the problem. I think instead of rebuilding OpenBLAS, it might be sufficient for you to export OPENBLAS_CORETYPE=HASWELL as a first test to exclude the AVX512 microkernel specifically) ?

I checked pretty carefully against reference LAPACK versions d134ca0070b0 (their latest master branch at time of writing), v3.12.1, and v3.9.1 linked to NumPy on x86_64 Linux (a few hours of work detailed at lanl/GFDL#117 (comment)), and found no evidence of the issue/reproduction. Let me know if I need to test another hash of reference LAPACK.

Ok, performing more analysis here today (this is mostly me, with some queries to an AI agent noted here and there below to identify likely problem routines/files at the low level). Building latest SciPy main branch (hash: 78f6a6a3e4b) against latest OpenBLAS develop hash (d3d4c6cd27689) at time of writing, and substituting scipy.linalg.lstsq in for the original reproducer reported at lanl/GFDL#117:

--- a/src/gfdl/model.py
+++ b/src/gfdl/model.py
@@ -700,7 +700,7 @@ class EnsembleGFDL(BaseEstimator):
             # = (D.T @ D)^-1 @ D.T @ y
 
             if self.reg_alpha is None:
-                coef_ = np.linalg.pinv(self.As[i], rtol=self.rtol) @ self.Bs[i]
+                coef_ = scipy.linalg.lstsq(self.As[i], self.Bs[i], cond=self.rtol)[0]

of course reproduces the failure exactly. That means that lapack_driver="gelsd" reproduces the problem, since it is the default. lapack_driver="gelsy" works "ok"/passes the test, while lapack_driver="gelss" fails identically to gelsd.

I then substituted with a roughly equivalent scipy.linalg.svd version of the algorithm,

Details

--- a/src/gfdl/model.py
+++ b/src/gfdl/model.py
@@ -700,7 +700,12 @@ class EnsembleGFDL(BaseEstimator):
             # = (D.T @ D)^-1 @ D.T @ y
 
             if self.reg_alpha is None:
-                coef_ = np.linalg.pinv(self.As[i], rtol=self.rtol) @ self.Bs[i]
+                A = self.As[i]
+                B = self.Bs[i]
+                U, s, Vt = scipy.linalg.svd(A, full_matrices=False, lapack_driver="gesdd")
+                tol = 1e-15 * s.max()
+                s_inv = np.where(s > tol, 1.0 / s, 0.0)
+                coef_ = Vt.T @ (s_inv[:, None] * (U.T @ B))
             else:
                 # scipy.linalg.solve(self.A + reg_mat, self.B)
                 # is equivalent to

which fails the test case identically with lapack_driver="gesdd". lapack_driver="gesvd" also fails identically here.

From here, I started asking an AI for "expert" views on the most likely OpenBLAS/low level routines that might be involved. The suggestion was that, based on these analyses, the "bug" is almost certainly before the bidiagonal SVD, with these candidates topping the list: dgebrd, dgebd2, dlabrd, dlarfb, dlarft, dlarf, and the dgemm/dgemv x86 kernels that get called.

I then performed a manual inspection of the bisection results I performed a while back at #5866 (comment). It narrowed down to two possible problem hashes.

The first hash, f5f789f, implements/changes tons of ?LARF1F (and 1L) stuff, which matches some of the likely problem routines flagged above. The second hash, d9bb8f3, also messes with a ton of ?LARF1F calls. These changes appear to assume first/last element constant unit (1), which may cause very subtle off-by-one bugs if there's even a single call site that isn't updated to reflect this. At least one fix related to this has already been merged, as noted at #5866 (comment), and that is gh-5865.

I then started analyzing the reference lapack sources related to *larf* source files compared to those vendored in OpenBLAS latest develop branch. It looks like OpenBLAS is missing a bug fix from reference LAPACK in this regard: Reference-LAPACK/lapack#1315. This was first identified upstream on June 19th, 5 days before I opened this issue--see: Reference-LAPACK/lapack#1314. However, pulling that upstream patch to the lapack-netlib/SRC/{c,d,s,z}larf.f OpenBLAS source had no effect whatsoever on my failing test case.

I next tried using a generic target on latest OpenBLAS develop (make TARGET=GENERIC ...)--that had no effect on the reproducing failure with scipy.linalg.lstsq. I then performed an experiment where I reverted only lapack-netlib/SRC/dgebd2.f to its state in OpenBLAS v0.3.30, which basically just reverts the DLARF -> DLARF1F update from upstream in that file. This finally allowed my test case with scipy.linalg.lstsq and the lapack_driver="gelsd" to pass.

I think I'm getting closer here now, I'll dig a little more to see if I can develop a patch.

LAPACK 1315 was merged as #5865 on June 26 so pulling that one again would be expected to have no effect.

I think we're still circling the one aspect of the DLARF -> DLARF1F transition that is more reliance on DAXPY and DGEMV - but it is curious that the fairly tame generic code would show the same "failure" I had previously attributed to (some of) the cpu-specific x86_64 assembly microkernels for GEMV

export OPENBLAS_CORETYPE=HASWELL has max rtol violation of 1.46476217, export OPENBLAS_CORETYPE=NEHALEM passes the test, export OPENBLAS_CORETYPE=SANDYBRIDGE passes the test, export OPENBLAS_CORETYPE=ZEN has a max rtol violation of 1.46476217.

This may point toward FMA business it seems, although my LLM interaction is getting rather verbose now. Perhaps this means there's not much that OpenBLAS can do?

export OPENBLAS_CORETYPE=ZEN has a max rtol violation of 1.46476217.

HASWELL and ZEN are basically the same thing under different names. ISTR that SKYLAKEX/COOPERLAKE passed the test too, haven't tried entirely different architectures like arm64 or riscv

This may point toward FMA business it seems, although my LLM interaction is getting rather verbose now. Perhaps this means there's not much that OpenBLAS can do?

Could well be, which is why I was harping on about making your code more robust against small numerical differences.
Unless, of course, there's some subtle fault in the AVX2 DGEMV microkernel that I haven't been able to find (and nobody has pointed out over the years). That code was written by a guy in his 60s who just "went away" from one day to the next, some months before the previous maintainer decided to drop out.
(And mind you I'm not maintaining OpenBLAS because I'm an assembly buff - I'm maintaining it because I was the only semi-regular contributor left at the time)

Closing here as ongoing activity in the linked GFDL ticket points to this being a algorithm stability issue with the handling of almost singular matrices, rather than a demonstrable bug in OpenBLAS

Navigation Menu

Repository navigation

BUG: Moore-Penrose Pseudoinverse numerical stability/consistency affected by recent OpenBLAS activity #5866

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

BUG: Moore-Penrose Pseudoinverse numerical stability/consistency affected by recent OpenBLAS activity #5866

Description

Activity

martin-frbg commented on Jun 24, 2026

tylerjereddy commented on Jun 24, 2026

martin-frbg commented on Jun 25, 2026

martin-frbg commented on Jun 25, 2026

martin-frbg commented on Jun 26, 2026

martin-frbg commented on Jun 26, 2026

martin-frbg commented on Jun 30, 2026

tylerjereddy commented on Jul 1, 2026

martin-frbg commented on Jul 1, 2026

martin-frbg commented on Jul 2, 2026

tylerjereddy commented on Jul 2, 2026

martin-frbg commented on Jul 2, 2026

tylerjereddy commented on Jul 2, 2026

martin-frbg commented on Jul 2, 2026

tylerjereddy commented on Jul 22, 2026

fchareyr commented on Aug 14, 2026

martin-frbg commented on Aug 14, 2026

tylerjereddy commented on Aug 23, 2026

tylerjereddy commented on Sep 6, 2026

martin-frbg commented on Sep 6, 2026

tylerjereddy commented on Sep 7, 2026

martin-frbg commented on Sep 7, 2026

martin-frbg commented on Oct 2, 2026

Metadata

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

Issue actions