Repository navigation
BUG: Matrix multiplication in windows produces large numerical errors inconsistently in numpy >= 2 #27036
Description
Describe the issue:
We discovered a simple test failing on windows only: sunpy/sunpy#7754
I have investigated the symptoms of this in some detail but have not tried to find the cause: In short it seems like matrix multiplications with largeish numbers fails inconsistently in windows, and when it does fail it will fail a few times in a row before returning to normal. I would suspect openblas but the failure never occurs with numpy <2
Reproduce the code example:
import numpy as np
def test():
# Test whether matrix multiplication involving a large matrix always gives the same answer
# This indirectly tests whichever BLAS/LAPACK libraries that NumPy is linking to (if any)
x = np.arange(500000, dtype=np.float64)
src = np.vstack((x, -10*x)).T
matrix = np.array([[0, 1], [1, 0]])
expected = np.vstack((-10*x, x)).T # src @ matrix
mismatches = np.zeros(500, int)
for i in range(len(mismatches)):
result = src @ matrix
mismatches[i] = (~np.isclose(result, expected)).sum()
if mismatches[i] != 0:
print(f"{mismatches[i]} mismatching elements in multiplication #{i}")Error message:
This bug is absolutely wild by the way:
With numpy 2.0.1:
16 mismatching elements in multiplication #22
316 mismatching elements in multiplication #33
28 mismatching elements in multiplication #53
307 mismatching elements in multiplication #100
32 mismatching elements in multiplication #201
32 mismatching elements in multiplication #267
288 mismatching elements in multiplication #272
1177 mismatching elements in multiplication #276
1596 mismatching elements in multiplication #289
32 mismatching elements in multiplication #298
16 mismatching elements in multiplication #314
3268 mismatching elements in multiplication #407
1403 mismatching elements in multiplication #419
64 mismatching elements in multiplication #446
with numpy 2.0.0
80 mismatching elements in multiplication #0
920 mismatching elements in multiplication #183
683 mismatching elements in multiplication #186
232 mismatching elements in multiplication #187
400 mismatching elements in multiplication #190
1217 mismatching elements in multiplication #195
1515 mismatching elements in multiplication #249
593 mismatching elements in multiplication #269
108 mismatching elements in multiplication #285
120 mismatching elements in multiplication #295
204 mismatching elements in multiplication #332
78 mismatching elements in multiplication #357
1816 mismatching elements in multiplication #380
164 mismatching elements in multiplication #396
48 mismatching elements in multiplication #492
with 1.26.4 it passes.
So it is a numpy 2 issue and not a 2.0.1 issue... it is also random:
Inserting a return and a break on the first time the matmul fails:
In [52]: test()
32 mismatching elements in multiplication #27
In [53]: test()
152 mismatching elements in multiplication #44
In [54]: test()
108 mismatching elements in multiplication #75
In [55]: test()
68 mismatching elements in multiplication #44
In [56]: test()
964 mismatching elements in multiplication #1
In [57]: test()
228 mismatching elements in multiplication #89
In [58]: test()
24 mismatching elements in multiplication #0
In a given set of multiplications the errors always have the same values.
In [92]: print([i for i in x if np.any(i != 0)])
[array([-69480., 6948.]), array([-69480., 6948.]), array([-69480., 6948.]), array([-69480., 6948.]), array([-69480., 6948.]), array([-69480., 6948.]), array([-69480., 6948.]), array([-69480., 6948.]), array([-69480., 6948.]), array([-69480., 6948.]), array([-69480., 6948.]), array([-69480., 6948.]), array([-69480., 6948.]), array([-69480., 6948.]), array([-69480., 6948.]), array([-69480., 6948.]), array([-69480., 6948.]), array([-69480., 6948.]), array([-69480., 6948.]), array([-69480., 6948.]), array([-69480., 6948.]), array([-69480., 6948.]), array([-69480., 6948.]), array([-69480.,
6948.])]
In [94]: print([i for i in (res-ex) if np.any(i != 0)])
[array([-1289600., 128960.]), array([-1289600., 128960.]), array([-1289600., 128960.]), array([-1289600., 128960.]), array([-1289600., 128960.]), array([-1289600., 128960.]), array([-1289600., 128960.]), array([-1289600., 128960.])]
In [101]: print([i for i in (res-ex) if np.any(i != 0)])
[array([-208450., 20845.]), array([-208450., 20845.]), array([-208450., 20845.]), array([-208450., 20845.]), array([-208450., 20845.]), array([-208450., 20845.]), array([-208450., 20845.]), array([-208450., 20845.])]
As you can see they also have the weird symmetry where the first column is -10* the second column which I guess makes some twisted sense because thats what the calculation is.
The errors are also always in the same region i.e. with indexs > 400k or so
They also always fail multiple indexes in a row, e.g. you might have 20 or 50 or whatever failed calcualtions in a row but then it starts working again.
e.g.
[417292, 417293, 417294, 417295, 417296, 417297, 417298, 417299, 417300, 417301, 417302, 417303, 449292, 449293, 449294, 449295, 449296, 449297, 449298, 449299, 449300, 449301, 449302, 449303, 449400, 449401, 449402, 449403, 449404, 449405, 449406, 449407, 449408, 449409, 449410, 449411, 449412, 449413, 449414, 449415, 449416, 449417, 449418, 449419, 449420, 449421, 449422, 449423]
or
[487736, 487737, 487738, 487739, 487740, 487741, 487742, 487743, 487760, 487761, 487762, 487763, 487764, 487765, 487766, 487767, 490490, 490491, 490492, 490493, 490494, 490495, 490496, 490497, 496838, 496839, 496840, 496841, 496842, 496843, 496844, 496845]
You can see in one of the plots there has been the same failure for 2 separate sequences of indices.
Python and NumPy Versions:
Windows 10, cannot recreate this at all in linux.
python 3.11.8
numpy 2.0.0 or numpy 2.0.1 fail in the same manner
This also failed on our CI on a windows-latest runner.
https://github.com/sunpy/sunpy/actions/runs/10072095800/job/27843572893#step:10:5160
Runtime Environment:
[{'numpy_version': '2.0.1',
'python': '3.11.8 (tags/v3.11.8:db85d51, Feb 6 2024, 22:03:32) [MSC v.1937 '
'64 bit (AMD64)]',
'uname': uname_result(system='Windows', node='DESKTOP-0EA2O83', release='10', version='10.0.19045', machine='AMD64')},
{'simd_extensions': {'baseline': ['SSE', 'SSE2', 'SSE3'],
'found': ['SSSE3',
'SSE41',
'POPCNT',
'SSE42',
'AVX',
'F16C',
'FMA3',
'AVX2'],
'not_found': ['AVX512F',
'AVX512CD',
'AVX512_SKX',
'AVX512_CLX',
'AVX512_CNL',
'AVX512_ICL']}},
{'architecture': 'Haswell',
'filepath': 'C:\\Users\\Alasdair\\AppData\\Local\\Programs\\Python\\Python311\\Lib\\site-packages\\numpy.libs\\libscipy_openblas64_-fb1711452d4d8cee9f276fd1449ee5c7.dll',
'internal_api': 'openblas',
'num_threads': 8,
'prefix': 'libscipy_openblas',
'threading_layer': 'pthreads',
'user_api': 'blas',
'version': '0.3.27'}]
Context for the issue:
Any matrix multiplication result could go wrong at any time in windows.
Activity
mattip commentedon Jul 25, 2024
It does seem to be a problem with OpenBLAS and threading on windows. NumPy 1.26.4 uses
v0.3.23-293. Numpy 2.0 and 2.0.1 usev0.3.27.44. SettingOPENBLAS_NUM_THREADS=1makes the test pass. I didn't see any relevant issues searching for windows and threading on the OpenBLAS issue tracker, but maybe I need to search differently. @martin-frbg does this sound familiar?mattip commentedon Jul 25, 2024
I wonder how our CI passed☹️
mattip commentedon Jul 25, 2024
I don't think this is related to the kernels chosen. The test still fails when I set
OPENBLAS_CORETYPE=PRESCOTT(it was HASWELL).martin-frbg commentedon Jul 25, 2024
what does 0.3.27.44 correspond to, git-wise ? If post-0.3.27, any chance this could be related to the NAN handling changes in SCAL that I merged on June 29 ?
ayshih commentedon Jul 25, 2024
For what it's worth, sunpy has this test because 4 years ago we encountered intermittent failures with OpenBLAS and multithreading on macOS (see sunpy/sunpy#4290 (comment)). Those failures went away eventually with a subsequent release of OpenBLAS.
ayshih commentedon Jul 25, 2024
Importantly, there were and maybe still are bad interactions when a threaded program used OpenBLAS that was configured to use threads as well. Hence, that's why setting OpenBLAS to use a single thread "fixed" issues.
martin-frbg commentedon Jul 25, 2024
Windows gets a lot less attention, and there were a few PRs merged since 0.3.23 (whatever build your -293 is) and 0.3.27-whatever that (supposedly) improved multithreading on that platform. Notably OpenMathLib/OpenBLAS#4359 - not trying to blame anyone, just that this would be the most serious change (merged in December, so debuting in 0.3.26)
martin-frbg commentedon Jul 27, 2024
@mseminatore
mseminatore commentedon Jul 27, 2024
Is it possible to test this against an earlier (pre 0.3.26) OpenBLAS version to isolate whether this could be related to the Windows threading changes? @martin-frbg do we have existing validation tests that would/should show this issue?
Does this show up for any NTHREADS > 1?
Do we have a C repro for this yet? If not I will attempt to translate the Python to C and see if I can reproduce/debug it.
martin-frbg commentedon Jul 27, 2024
I have just managed to set up a virtual environment for reproducing this (with a pip-installed numpy - I did not get conda to install anything other than mkl or netlib even if it claimed to have switched to openblas). It definitely happens with as few as 4 cores (though not on every run), but I have not yet gotten it to fail with just two. (Unfortunately I haven't figured out how to downgrade the openblas version, pip tells me there are no packages that match the version requirement)
martin-frbg commentedon Jul 27, 2024
3 threads are definitely enough to reproduce it (on the fourth try), I still haven't seen a failure with 2 in about 20 runs so far
mattip commentedon Jul 27, 2024
The openblas implementation is in
<venv>/lib/site-packages/numpy.libs. You can try various scipy-openblas64 versions. Note that pip-installing the package does not replace the DLL used by numpy: the package is only used when building numpy from source. You can copy the DLL from the site-package/scipy_openblas/lib directory, overwriting the one in numpy.libs (be sure to preserve the mangled name).mattip commentedon Jul 27, 2024
Note I am close to uploading a newer scipy-openblas64 based on latest OpenBLAS HEAD. You can download a pre-release from https://anaconda.org/scientific-python-nightly-wheels/scipy-openblas64/files. You want scipy_openblas64-0.3.27.341.0-py3-none-win_amd64.whl.
37 remaining items
martin-frbg commentedon Aug 2, 2024
thanks - my bisect is four steps from completion, I am also building a current HEAD with 4359 manually reverted now "just in case".
martin-frbg commentedon Aug 2, 2024
Indeed HEAD with the changes from 4359 reverted appears to work fine - whether that means the problem is in 4359 or only uncovered by it remains to be seen. (Bisect still needs two steps but appears to be gravitating towards the same result)
mseminatore commentedon Aug 2, 2024
Thank you @mattip abd @martin-frbg for the bisect work. I am sorry that I can’t more actively participate (for personal reasons) but will continue inspecting the code to see if I can identify a cause and fix.
martin-frbg commentedon Aug 3, 2024
I think it is equally possible that the actual problem lies somewhere in driver/level3/level3_thread.c (though it has its own CriticalSection), but I lack the Windows developer experience to assess this
mseminatore commentedon Aug 3, 2024
I haven’t looked very closely at that code but I will take a look.
mattip commentedon Aug 4, 2024
Any further thoughts about #OpenMathLib/OpenBLAS#4835 and whether NumPy should build its OpenBLAS with that PR as a patch to move the NumPy2.1 release forward?
martin-frbg commentedon Aug 4, 2024
not sure if that counts as thoughts, but my current intention is to merge 4835 in order to get the 0.3.28 release out in the next couple of days (probably not today though)
mattip commentedon Aug 8, 2024
I backported the changes in OpenBLAS's 4835 in the 0.3.27.44.4 version of scipy-openblas, and used those wheels in #27140. I also added the test to prevent this from regressing again. Closing. Please reopen or open a new issue if I missed something.
charris commentedon Dec 11, 2025
There was another failure of the test for this on Python 3.11: https://github.com/numpy/numpy/actions/runs/20146759701/job/57829132899
mseminatore commentedon Dec 12, 2025
I guess we can safely rule out PR #4359 for this one.
mattip commentedon Dec 13, 2025
Wierd. That seems to be the only failure since this issue was closed, right? And it does not fail consistently.
mseminatore commentedon Dec 13, 2025
It's the only issue that I've seen since then. I felt this at the time when 4359 was in play, I think we need some stress tests to help find low probability and edge case usage issues. I think it was ~9 months after 4359 was merged before an issue was found.
Since we never found the root cause I don't know that we know for certain whether 4359 was cause or whether the increased efficiency (removing millions of kernel object create/deletes) on Windows helped to raise the probability of a latent issue.
If there is interest in expanding stress test coverage, I may be interested in helping. I would still like to see the Windows server code improved and believe that is necessary before proposing fundamental changes again.
mattip commentedon Dec 14, 2025
We could add more tests to OpenBLAS itself or to MacPython/openblas-libs (the build of OpenBLAS used by SciPy and NumPy). FWIW, I added the test at the top of this issue to NumPy in #27140. What do you suggest for testing?