np.linalg.svd gives weird output

[mortonjt]

It appears that numpy.linalg.svd has different outputs between 1.15.4 and 1.16.3 -- and I don't see any hints that this should happen in the changelog.

Reproducing code example:

If you run numpy==1.15.4 with the following input, you get the following

In [20]: np.random.seed(0)

In [21]: X = np.random.randn(100, 100) * 10

In [22]: u, s, v = np.linalg.svd(X)

In [23]: s[:5]
Out[23]:
array([193.6995948 , 183.01489365, 180.94983595, 178.47342134,
       174.07146491])

In [24]: u[:5, :5]
Out[24]:
array([[-0.12777193,  0.1473861 , -0.04626332,  0.09269663,  0.12149756],
       [-0.23535239,  0.06543527, -0.02588301, -0.05401032, -0.02436388],
       [ 0.09956489, -0.07108259, -0.10083294,  0.16710062, -0.0333519 ],
       [-0.05718053, -0.11184903, -0.12280465,  0.04707958,  0.16647603],
       [ 0.25422798,  0.01464547, -0.15071626, -0.03225463, -0.12620621]])

In [25]: np.mean(X - u @ np.diag(s) @ v)
Out[25]: 1.336599234069702e-16

In [26]: np.std(X - u @ np.diag(s) @ v)
Out[26]: 2.0947185697750298e-14

However, if you run svd with numpy==1.16.3, you get the following

In [13]: np.random.seed(0)

In [14]: X = np.random.randn(100, 100) * 10

In [15]: u, s, v = np.linalg.svd(X)

In [16]: s[:5]
Out[16]:
array([193.6995948 , 183.01489365, 180.94983595, 178.47342134,
       174.07146491])

In [17]: u[:5, :5]
Out[17]:
array([[-0.07680886,  0.17809124, -0.41179481,  0.14073815,  0.14494919],
       [ 0.46612298, -0.38573887,  0.34177276, -0.11314353, -0.03170648],
       [-0.62584307,  0.10539512,  0.20991414, -0.24513665, -0.11505925],
       [-0.43937521, -0.07100245,  0.2941617 , -0.03515797,  0.04093411],
       [ 0.00359183, -0.64178218,  0.69484207, -0.09871542,  0.28272959]])

In [18]: np.mean(X - u @ np.diag(s) @ v)
Out[18]: 1.0356975839816995

In [19]: np.std(X - u @ np.diag(s) @ v)
Out[19]: 92.84521345330656

Numpy/Python version information:

numpy version 1.15

In [18]: np.__config__.show()
mkl_info:
    libraries = ['mkl_rt', 'pthread']
    library_dirs = ['/Users/jmorton/miniconda3/envs/qiime2-2019.1/lib']
    define_macros = [('SCIPY_MKL_H', None), ('HAVE_CBLAS', None)]
    include_dirs = ['/Users/jmorton/miniconda3/envs/qiime2-2019.1/include']
blas_mkl_info:
    libraries = ['mkl_rt', 'pthread']
    library_dirs = ['/Users/jmorton/miniconda3/envs/qiime2-2019.1/lib']
    define_macros = [('SCIPY_MKL_H', None), ('HAVE_CBLAS', None)]
    include_dirs = ['/Users/jmorton/miniconda3/envs/qiime2-2019.1/include']
blas_opt_info:
    libraries = ['mkl_rt', 'pthread']
    library_dirs = ['/Users/jmorton/miniconda3/envs/qiime2-2019.1/lib']
    define_macros = [('SCIPY_MKL_H', None), ('HAVE_CBLAS', None)]
    include_dirs = ['/Users/jmorton/miniconda3/envs/qiime2-2019.1/include']
lapack_mkl_info:
    libraries = ['mkl_rt', 'pthread']
    library_dirs = ['/Users/jmorton/miniconda3/envs/qiime2-2019.1/lib']
    define_macros = [('SCIPY_MKL_H', None), ('HAVE_CBLAS', None)]
    include_dirs = ['/Users/jmorton/miniconda3/envs/qiime2-2019.1/include']
lapack_opt_info:
    libraries = ['mkl_rt', 'pthread']
    library_dirs = ['/Users/jmorton/miniconda3/envs/qiime2-2019.1/lib']
    define_macros = [('SCIPY_MKL_H', None), ('HAVE_CBLAS', None)]
    include_dirs = ['/Users/jmorton/miniconda3/envs/qiime2-2019.1/include']
In [27]:  import sys, numpy; print(numpy.__version__, sys.version)
1.15.4 3.6.5 | packaged by conda-forge | (default, Apr  6 2018, 13:44:09)
[GCC 4.2.1 Compatible Apple LLVM 6.1.0 (clang-602.0.53)]

numpy version 1.16

In [3]: np.__config__.show()
blas_mkl_info:
  NOT AVAILABLE
blis_info:
  NOT AVAILABLE
openblas_info:
    libraries = ['openblas', 'openblas']
    library_dirs = ['/usr/local/lib']
    language = c
    define_macros = [('HAVE_CBLAS', None)]
blas_opt_info:
    libraries = ['openblas', 'openblas']
    library_dirs = ['/usr/local/lib']
    language = c
    define_macros = [('HAVE_CBLAS', None)]
lapack_mkl_info:
  NOT AVAILABLE
openblas_lapack_info:
    libraries = ['openblas', 'openblas']
    library_dirs = ['/usr/local/lib']
    language = c
    define_macros = [('HAVE_CBLAS', None)]
lapack_opt_info:
    libraries = ['openblas', 'openblas']
    library_dirs = ['/usr/local/lib']
    language = c
    define_macros = [('HAVE_CBLAS', None)]

In [20]:  import sys, numpy; print(numpy.__version__, sys.version)
1.16.3 3.6.7 | packaged by conda-forge | (default, Feb 28 2019, 02:16:08)
[GCC 4.2.1 Compatible Clang 4.0.1 (tags/RELEASE_401/final)]

As you can see, the output from 1.16.3 is broken, and the svd is no longer a true reconstruction.

[charris]

This is probably a conda issue related to MKL. Can you check if you are linked to MKL? You might be able to see that with numpy.__config__.show(). BTW, it helps for reproducting these examples if you explicitly set the random generator seed.

[mortonjt]

@charris that would make a lot of sense. Except the 1.16.3 version was pip installed (within a conda environment) -- config is attached. The examples have been updated with random seeds (sorry about that ...).

[mortonjt]

It also looks like scipy is having similar issues

In [12]: from scipy.linalg import svd

In [16]: np.random.seed(0)

In [17]: X = np.random.randn(100, 100) * 10

In [18]: u, s, v = svd(X)

In [19]: np.mean(X - u @ np.diag(s) @ v)
Out[19]: 1.0356975839819398

In [20]: np.std(X - u @ np.diag(s) @ v)
Out[20]: 92.84521345330933

In [22]: import scipy

In [23]: scipy.__version__
Out[23]: '1.2.1'

[charris]

I can't duplicate the problem on linux with pip installed 1.16.3. @eric-wieser Do you see this on Windows? @mortonjt Are you running 32 or 64 bits python.

[charris]

Note that pip and conda do not play well together. It would be helpful if you could duplicate the results in a non-conda environment.

[mattip]

The 1.15.4 results are using mkl, the 1.16.4 are without any acceleration. They also use different clang and python

1.15.4 3.6.5 | packaged by conda-forge | (default, Apr  6 2018, 13:44:09)
[GCC 4.2.1 Compatible Apple LLVM 6.1.0 (clang-602.0.53)]

The 1.16.4 results seem to use a different clang

1.16.3 3.6.7 | packaged by conda-forge | (default, Feb 28 2019, 02:16:08)
[GCC 4.2.1 Compatible Clang 4.0.1 (tags/RELEASE_401/final)]

[mortonjt]

@charris , I was able to reproduce the results with just pip outside of my conda environment

>>> import numpy as np
>>> np.__version__
'1.16.3'
>>> np.random.seed(0)
>>> X = np.random.randn(100, 100) * 10
>>> u, s, v = np.linalg.svd(X)
>>> np.mean(X - u @ np.diag(s) @ v)
1.0356975839816995
>>> np.std(X - u @ np.diag(s) @ v)
92.84521345330656
>>> np.__config__.show()
blas_mkl_info:
  NOT AVAILABLE
blis_info:
  NOT AVAILABLE
openblas_info:
    libraries = ['openblas', 'openblas']
    library_dirs = ['/usr/local/lib']
    language = c
    define_macros = [('HAVE_CBLAS', None)]
blas_opt_info:
    libraries = ['openblas', 'openblas']
    library_dirs = ['/usr/local/lib']
    language = c
    define_macros = [('HAVE_CBLAS', None)]
lapack_mkl_info:
  NOT AVAILABLE
openblas_lapack_info:
    libraries = ['openblas', 'openblas']
    library_dirs = ['/usr/local/lib']
    language = c
    define_macros = [('HAVE_CBLAS', None)]
lapack_opt_info:
    libraries = ['openblas', 'openblas']
    library_dirs = ['/usr/local/lib']
    language = c
    define_macros = [('HAVE_CBLAS', None)]
>>> import sys, numpy; print(numpy.__version__, sys.version)
1.16.3 3.7.3 (default, Mar 27 2019, 09:23:39)
[Clang 10.0.0 (clang-1000.11.45.5)]

@mattip - I don't think it is compiler specific, because in that same environment with numpy==1.15.2, I get correct results

>>> import numpy as np
>>> np.random.seed(0)
>>> X = np.random.randn(100, 100) * 10
>>> u, s, v = np.linalg.svd(X)
>>> np.mean(X - u @ np.diag(s) @ v)
1.0218109995283053e-16
>>> np.std(X - u @ np.diag(s) @ v)
2.1302064719026226e-14
>>> np.__config__.show()
blas_mkl_info:
  NOT AVAILABLE
blis_info:
  NOT AVAILABLE
openblas_info:
    libraries = ['openblas', 'openblas']
    library_dirs = ['/usr/local/lib']
    language = c
    define_macros = [('HAVE_CBLAS', None)]
blas_opt_info:
    libraries = ['openblas', 'openblas']
    library_dirs = ['/usr/local/lib']
    language = c
    define_macros = [('HAVE_CBLAS', None)]
lapack_mkl_info:
  NOT AVAILABLE
openblas_lapack_info:
    libraries = ['openblas', 'openblas']
    library_dirs = ['/usr/local/lib']
    language = c
    define_macros = [('HAVE_CBLAS', None)]
lapack_opt_info:
    libraries = ['openblas', 'openblas']
    library_dirs = ['/usr/local/lib']
    language = c
    define_macros = [('HAVE_CBLAS', None)]
>>> import sys, numpy; print(numpy.__version__, sys.version)
1.15.4 3.7.3 (default, Mar 27 2019, 09:23:39)
[Clang 10.0.0 (clang-1000.11.45.5)]

[charris]

Looks like it is probably an openblas issue on windows. We have a different openblas version for the forthcoming 1.17 and there are wheels for that at https://7933911d6844c6c53a7d-47bd50c35cd79bd838daf386af554a83.ssl.cf2.rackcdn.com/, although the most recent windows version seems to be from April 2. If you could try that as well it would be helpful.

[mortonjt]

@charris the errors above were run on a mac.

But I can confirm that this error will not arise on linux - see below for output.

In [1]: import numpy as np

In [2]: np.random.seed(0)

In [3]: X = np.random.randn(100, 100) * 10

In [4]: u, s, v = np.linalg.svd(X)

In [5]: np.mean(X - u @ np.diag(s) @ v)
Out[5]: 9.420362710385599e-17

In [6]: np.std(X - u @ np.diag(s) @ v)
Out[6]: 2.086820837986988e-14

In [7]: np.__config__.show()
blas_mkl_info:
  NOT AVAILABLE
blis_info:
  NOT AVAILABLE
openblas_info:
    libraries = ['openblas', 'openblas']
    library_dirs = ['/usr/local/lib']
    language = c
    define_macros = [('HAVE_CBLAS', None)]
blas_opt_info:
    libraries = ['openblas', 'openblas']
    library_dirs = ['/usr/local/lib']
    language = c
    define_macros = [('HAVE_CBLAS', None)]
lapack_mkl_info:
  NOT AVAILABLE
openblas_lapack_info:
    libraries = ['openblas', 'openblas']
    library_dirs = ['/usr/local/lib']
    language = c
    define_macros = [('HAVE_CBLAS', None)]
lapack_opt_info:
    libraries = ['openblas', 'openblas']
    library_dirs = ['/usr/local/lib']
    language = c
    define_macros = [('HAVE_CBLAS', None)]

In [8]: import sys, numpy; print(numpy.__version__, sys.version)
1.16.3 3.6.2 (default, Aug  2 2017, 14:29:03)
[GCC 5.4.0]

[charris]

Hmm. @tylerjereddy @mathew-brett Is there anything special about how OpenBLAS is compiled on the Mac? Compiler used, that sort of thing.

[matthew-brett]

I believe it uses clang - flags etc here: https://github.com/MacPython/openblas-libs/blob/master/travis-ci/build_steps.sh#L81

I wonder if it's an interaction with the OpenBLAS kernel. @mortonjt - can you try some different values for the OpenBLAS kernel, with:

export OPENBLAS_VERBOSE=2
export OPENBLAS_CORETYPE=haswell

in the shell, before starting Python?

My Mac laptop and Mac desktop give correct 1e-15ish results, choosing default values of 'haswell' and 'sandybridge' for the OpenBLAS kernel, respectively.

Examples of values for OPENBLAS_CORETYPE in https://github.com/matthew-brett/manylinux-testing/blob/master/.travis.yml

[tylerjereddy]

For diagnostic purposes, with 1.16.3, you should also be able to directly probe the OpenBLAS config data as well. From Matti's code block there:

import numpy, ctypes, os
dll = ctypes.CDLL(numpy.core._multiarray_umath.__file__)
get_config = dll.openblas_get_config
get_config.restype=ctypes.c_char_p
res = get_config()
res

On Linux I get something like this: b'OpenBLAS 0.3.5.dev DYNAMIC_ARCH NO_AFFINITY Haswell MAX_THREADS=64'
I get the same config result on Mac & I can't reproduce your original issue. Maybe providing us with that ctypes output would be helpful too.

[matthew-brett]

I can't get a failure on any kernel on my laptop or desktop:

#!/bin/bash
cat << EOF > np_check.py
import numpy as np
print(np.__version__)
np.random.seed(0)
X = np.random.randn(100, 100) * 10
u, s, v = np.linalg.svd(X)
print(np.mean(X - u @ np.diag(s) @ v))
print(np.std(X - u @ np.diag(s) @ v))
EOF
export OPENBLAS_CORETYPE
export OPENBLAS_VERBOSE=2
for kernel in penryn sandybridge prescott nehalem atom haswell; do
    echo Testing with kernel $kernel
    export OPENBLAS_CORETYPE=$kernel
    python np_check.py
done

With Tyler's script I get:

b'OpenBLAS 0.3.5.dev DYNAMIC_ARCH NO_AFFINITY Core2 MAX_THREADS=64'

(Laptop)

b'OpenBLAS 0.3.5.dev DYNAMIC_ARCH NO_AFFINITY Sandybridge MAX_THREADS=64'

(Desktop)

[matthew-brett]

@mortonjt - what do you get from my script above?

[charris]

I wonder if there is another OpenBLAS library floating around? @mortonjt Could you do a file search for '*openblas*'.

[mortonjt]

@matthew-brett , I think you are on to something. When setting the flags,

export OPENBLAS_VERBOSE=2
export OPENBLAS_CORETYPE=haswell

I get the right answer. And when I run the script provided, it seems to nail the right answer every time.

Testing with kernel penryn
Core: Core2
1.16.3
3.666055848651484e-16
2.1619569755102583e-14
Testing with kernel sandybridge
Core: Sandybridge
1.16.3
4.3135485417941654e-16
2.2037472763482323e-14
Testing with kernel prescott
Core: Prescott
1.16.3
3.654898107254001e-16
2.162364554532177e-14
Testing with kernel nehalem
Core: Atom
1.16.3
1.424538330178915e-16
2.1591982659274054e-14
Testing with kernel atom
Core: Atom
1.16.3
1.424538330178915e-16
2.1591982659274054e-14
Testing with kernel haswell
Core: Haswell
1.16.3
1.0218109995283053e-16
2.1302064719026226e-14

So, I'm guessing I had a faulty openblas library floating around somewhere. Still not sure what is going on there, but I'm going to close this issue, since it looks like I have a working solution now.

Thank you @matthew-brett @charris @tylerjereddy for the detective work!

[matthew-brett]

@mortonjt - just to check - what do you get after:

export OPENBLAS_VERBOSE=2
unset OPENBLAS_CORETYPE
python np_check.py

What do you get for Tyler's script? It would be good to be sure what the original problem was.

[mortonjt]

Below are the results for the np_check.py

python3 np_check.py
Core: SkylakeX
1.16.3
1.0356975839816995
92.84521345330656

And here are the results for Tyler's script

>>> import numpy, ctypes, os
Core: SkylakeX
>>> dll = ctypes.CDLL(numpy.core._multiarray_umath.__file__)
>>> get_config = dll.openblas_get_config
>>> get_config.restype=ctypes.c_char_p
>>> res = get_config()
>>> res
b'OpenBLAS 0.3.5.dev DYNAMIC_ARCH NO_AFFINITY SkylakeX MAX_THREADS=64'

[matthew-brett]

OK - so I suppose this must be an error in the OpenBLAS SkylakeX kernel - that I can't test - I don't have the hardware.

Are you running Linux on the same machine? Can you test on Linux with:

export OPENBLAS_VERBOSE=2
export OPENBLAS_CORETYPE=skylakex
python np_check.py

?

[mortonjt]

I'm running linux on another machine - and it doesn't seem to have the openblas skylakeX kernel

[charris]

There is a reported skylake openblas bug that should be fixed in 3.5. So maybe this is related.

[matthew-brett]

Chuck - which bug were you thinking of? I was just looking at OpenMathLib/OpenBLAS#1955 . I think the relevant kernel is being disabled in OpenBLAS 0.3.6 - OpenMathLib/OpenBLAS#2061 .

[charris]

@matthew-brett I didn't have a particular one in mind, but I note that one of the referenced issues in your link involves pinv (Julia), so that is probably svd related also.

[charris]

Sounds like we want OpenBLAS 3.6 until the skylake bugs are ironed out, probably for NumPy 1.16.4 also.

[tylerjereddy]

@stefanv was asking about unit testing for this issue. I suspect we do have tests that would fail here (perhaps many things in linalg?), but they don't fail in CI because of a specific CPU architecture requirement to induce the failure? Any creative way we could detect the issue in CI with tests, and is that worth the effort?

@rgommers It looks like this means that SciPy 1.2.1 / 1.3.x are hosed for some linear algebra calculations on Skylake, since the user report above is for 1.2.1 & I don't think I've bumped OpenBLAS for 1.3.x series (yet anyway). Proof / regression guards with tests would also be better for SciPy of course.

But do we really want to invest effort in generating creative tests that somehow guard against regressions specific to N architectures upstream?

[matthew-brett]

We can try to run specific blas tests with different archictures, as set by the OPENBLAS_CORETYPE variable, and that seems like a good idea to me - but we can only test the architectures supported by the CI virtual machines, which might not include SkylakeX, for example. We could also ask people with specific hardware to test for us, before release.

[rgommers]

they don't fail in CI because of a specific CPU architecture requirement to induce the failure?

looks like it

Any creative way we could detect the issue in CI with tests, and is that worth the effort?

not easy, and no I don't think that that is worth the effort to investigate. there's no way we can cover all common architectures, and we anyway have no control over what hardware CI runs on.

@mortonjt could you check if numpy.linalg.test() and scipy.linalg.test() fail for you? If not, we should add these as regression tests.

We can try to run specific blas tests with different archictures, as set by the OPENBLAS_CORETYPE variable, and that seems like a good idea to me - but we can only test the architectures supported by the CI virtual machines, which might not include SkylakeX, for example. We could also ask people with specific hardware to test for us, before release.

Really this is an OpenBLAS issue though. It is with pretty much every release. Rather than fix those upstream of OpenBLAS, it would make more sense to help OpenBLAS improve its CI and release QA. Then these things are a lot less likely to hit us in the first place, and it also helps other libraries and communities (R for example).

We can also try to downgrade back to the same OpenBLAS version we used for 1.15.3, but I guess we upgraded to avoid some other bug........

[matthew-brett]

Yes, but we use OpenBLAS for all our binary releases, so it becomes our problem too. It's surely easier for us to make sure that OpenBLAS is working for Numpy, than for us to make sure that they are doing it.

Actually, they have been running some tests here:

http://build.openblas.net/builders

But it doesn't look as if they are testing Skylake.

[rgommers]

It's surely easier for us to make sure that OpenBLAS is working for Numpy, than for us to make sure that they are doing it.

I don't think it is. At least not in CI. We can push people to test if it works for them with a release candidate, but that's apparently not working and if we find a problem, it then becomes a release blocker straight away and there may not be a working OpenBLAS version at all that we can switch to.

detecting problems where they're created is almost always the best idea ....

[matthew-brett]

I'm sure that's better in the abstract, but the issue is agency. We've got lots of agency testing Numpy, but much less testing OpenBLAS. It's difficult for me to imagine that any of us are going to put lots of effort into - for example - trying to keep the OpenBLAS buildbots maintained. So I suspect that won't happen. Meaning that we're at risk of releasing wheels that give wrong results, unless we check. And it's pretty easy for us to check - maybe something like:

export OPENBLAS_CORETYPE
export OPENBLAS_VERBOSE=2
for kernel in penryn sandybridge prescott nehalem atom haswell skylake skylakex; do
    echo Testing with kernel $kernel
    export OPENBLAS_CORETYPE=$kernel
    python -c 'import numpy.linalg as npl; npl.test()'
done

It could go, for example, in nightly wheel tests.

[rgommers]

That very likely won't work because the hardware doesn't support all that (no control on what you get with TravisCI & co). If it was that easy, OpenBLAS wouldn't need buildbots that run outside regular CI right?

[matthew-brett]

Try it on your own machine. I think you'll find that all the kernels work up until Skylake, where they might segfault, trying to use assembly instructions your CPU does not have. So the trick is only working out what kernels the CI machine will run, and then running those. From time to time, we can prod the machines to see if we can get away with another kernel. I imagine the chances of losing support for a kernel in a VM update is pretty small.

For example, I used this trick a lot for the initial manylinux / OpenBLAS testing.

So no, you don't need special hardware for most of the kernels, but you will for the most cutting edge ones, like SkylakeX.

[rgommers]

Try it on your own machine.

I have a Ryzen CPU here ....

So no, you don't need special hardware for most of the kernels, but you will for the most cutting edge ones, like SkylakeX.

The most cutting edge ones are the ones most likely to be a problem though probably. At least it is here. Anyway, if it is that easy then we should do that. Just doesn't solve the real issue imho, so we should also do something more systematic.

[matthew-brett]

It would be good to try and gather interest for supporting wider testing by OpenBLAS, as well. I guess the Julia people will be interested, they found the SkylakeX bug first. But I suppose we might be OpenBLAS's biggest customers, in terms of users.

[rgommers]

Yep, we have the biggest user base, but Julia may have as much or more developer time available. R uses OpenBLAS too IIRC. There really aren't too many flavors to choose from ....