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?