Do numpy and scipy use different LAPACK drivers for eigvalsh ? 3 Jul 2020

0-eigvals-numpy-scipy.md

Do numpy and scipy use different LAPACK drivers for eigvalsh on macos ?

eigvals-numpy-scipy.py below runs numpy vs scipy eigvalsh on a dense random matrix:

np.linalg.eigvalsh: 201 sec
scipy.linalg.eigvalsh: 83 sec  driver ev
scipy.linalg.eigvalsh: 84 sec  driver evd

Maybe this is only on macos ? (Gnutime: 381% of 4 cores, so it's not that.)

A suggested doc change, if it's correct --

the numpy.linalg functions all work on 64-bit floating vectors and arrays; 32-bit input arrays (dtype np.float32) are silently converted to np.float64. The corresponding scipy.linalg functions work on either, so e.g. scipy.linalg.eigvalsh( A.astype( np.float32 )) may run twice as fast as ( A ).

cheers
-- denis

16jul-eigh-numpy-vs-scipy-randomsparse-n10000.log

# from: numpy-vs-scipy-linalg.py run="eig.*h" n=10000 test="randomsparse" incdiag=0 dtype=np.float64
# 16 Jul 2020 16:13z  ~bz/py/np/sparse/eig/issues/dense  Denis-iMac 10.10.5

▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄
versions: numpy 1.19.0  scipy 1.5.1  python 3.7.6  plaform Darwin-14.5.0-x86_64-i386-64bit 
scipy lapack ilaver: (3, 9, 0)
randomsparse A (10000, 10000) density 0.001 float64 order C  incdiag 0
800 mb dense
np_eigvalsh: 201 sec  [-4.02 -4.01 -3.97 ... 4.01 4.03 5.75]
sp_eigvalsh_ev: 83 sec  [-4.02 -4.01 -3.97 ... 4.01 4.03 5.75]
sp_eigvalsh_evd: 84 sec  [-4.02 -4.01 -3.97 ... 4.01 4.03 5.75]
sp_eigvalsh_evr: 85 sec  [-4.02 -4.01 -3.97 ... 4.01 4.03 5.75]
np_eigh: 148 sec  [-4.02 -4.01 -3.97 ... 4.01 4.03 5.75]
sp_eigh_evd: 149 sec  [-4.02 -4.01 -3.97 ... 4.01 4.03 5.75]
sp_eigh_evr: 162 sec  [-4.02 -4.01 -3.97 ... 4.01 4.03 5.75]

919.04 real time  3356.06 user  65.30 sys  372% of 4 cores  imac 2.7 GHz

16jul-numpy-vs-scipy-linalg.py

#!/usr/bin/env python
""" time numpy / scipy linalg eigvalsh solve ... on sparse.random .A, dense """
# https://github.com/numpy/numpy/issues/16741
# scipy lapack-3.9.0.tar.gz 2019, numpy lapack-lite-3.1.1.tgz 2007
# http://www.netlib.org/lapack  https://github.com/Reference-LAPACK/lapack/issues

from collections import OrderedDict
from functools import partial
import platform
import re
import sys
from time import time
import numpy as np
import scipy
import scipy.linalg, scipy.linalg.lapack, scipy.sparse

try:
    from testmatrix import testmatrix  # powerlaw, suitesparse ...
except ImportError:
    testmatrix = None

__version__ = "2020-07-16 July  denis-bz-py t-online.de"

np.set_printoptions( threshold=10, edgeitems=3, linewidth=120,
        formatter = dict( float = lambda x: "%.3g" % x,  # float arrays %.3g
                        complexfloat = " {0.real:.3g} {0.imag:+.3g}j " .format ))

print( 80 * "▄" )
print( "versions: numpy %s  scipy %s  python %s  plaform %s " % (
    np.__version__, scipy.__version__ , sys.version.split()[0], platform.platform() ))
print( "scipy lapack ilaver:", scipy.linalg.lapack.ilaver() )  # (3, 9, 0)

#...........................................................................
run = "."  # "eigvals" ... match np_sp_linalg keys below
test = "randomsparse"
# test = "Norris-fv1"  # suitesparse https://sparse.tamu.edu/Norris
n = 10000
density = .001
seed = 0
incdiag = .01
dtype = np.float64
order = 'C'  # 'F' same
rcond = 1e-6

# to change these params, run this.py  a=1  b=None  'c = expr' ... in sh or ipython --
for arg in sys.argv[1:]:
    exec( arg )

#...........................................................................
if test == "randomsparse" or testmatrix is None:
    Asparse = scipy.sparse.random( n, n, density=density, random_state=seed )  # uniform 0 - 1
else:
    Asparse = testmatrix( test, n=n )
Asparse += incdiag * scipy.sparse.eye( n )  # shift away from singular, wikipedia Circular_law
A = Asparse.A.astype( dtype=dtype, order=order )  # dense
print( "%s A %s density %.3g %s order %s  incdiag %g" % (
        test, A.shape, density, A.dtype, order, incdiag ))
print( "%.0f mb dense" % (A.nbytes / 1e6 ))
b = np.ones( n, dtype=dtype )

#...........................................................................
    # run some or all of these, e.g. run="solve"
np_sp_linalg = OrderedDict(
    np_eigvalsh = np.linalg.eigvalsh,  # _syevd
    sp_eigvalsh_ev = partial( scipy.linalg.eigvalsh, driver="ev" ),  # ev evd evr evx
    sp_eigvalsh_evd = partial( scipy.linalg.eigvalsh, driver="evd" ),
    sp_eigvalsh_evr = partial( scipy.linalg.eigvalsh, driver="evr" ),

        # evecs too --
    np_eigh = np.linalg.eigh,
    sp_eigh_evd = partial( scipy.linalg.eigh, driver="evd" ),
    sp_eigh_evr = partial( scipy.linalg.eigh, driver="evr" ),
        # ev evd evr evx / gv gvd gvx generalized

        # complex evals --
    np_eigvals = np.linalg.eigvals,  # _geev
    sp_eigvals = scipy.linalg.eigvals,

    np_lstsq = partial( np.linalg.lstsq, b=b, rcond=rcond ),
    sp_lstsq = partial( scipy.linalg.lstsq, b=b, cond=rcond ),

    np_solve = partial( np.linalg.solve, b=b ),
    sp_solve = partial( scipy.linalg.solve, b=b ),
    np_svd = partial( np.linalg.svd, compute_uv=False ),  # gesdd
    sp_svd = partial( scipy.linalg.svd, compute_uv=False ),  # lapack_driver : {'gesdd', 'gesvd'}
)

for name, func in np_sp_linalg.items():
    if not re.search( run, name ):
        continue
    t0 = time()

    x = func( A )
    if isinstance( x, tuple ):
        xtup = x
        x = x[0]
    t = time() - t0
    x = x.real
    print( "%s: %.0f sec  %s" % (name, t, np.sort( x )))

3jul-eigvals-numpy-scipy-randomsparse-n10000.log

# from: eigvals-numpy-scipy.py n=10000 test="randomsparse" dtype=np.float64
#  3 Jul 2020 10:16z  ~bz/py/np/sparse/eig/dense  Denis-iMac 10.10.5

▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄
versions: numpy 1.19.0  scipy 1.5.0  python 3.7.6  plaform Darwin-14.5.0-x86_64-i386-64bit 
randomsparse A (10000, 10000) density 0.001 float64: dense 800 mb
np.linalg.eigvalsh: 201 sec  [-4.02 -4.01 -3.97 -3.96 -3.95 ... 3.96 3.97 4.01 4.03 5.75]
scipy.linalg.eigvalsh: 83 sec  driver ev  [-4.02 -4.01 -3.97 -3.96 -3.95 ... 3.96 3.97 4.01 4.03 5.75]
scipy.linalg.eigvalsh: 84 sec  driver evd  [-4.02 -4.01 -3.97 -3.96 -3.95 ... 3.96 3.97 4.01 4.03 5.75]
max |npevals - scipyevals|: 8.5e-14

375.02 real time  1412.24 user  18.74 sys  381% of 4 cores  imac 2.7 GHz

3jul-eigvals-numpy-scipy.py

#!/usr/bin/env python
""" time eigvalsh numpy / scipy of sparse.random() .A, dense """
# wibni, wouldn't it be nice if numpy and scipy eig* used the same default LAPACK drivers

from __future__ import division, print_function
import platform
import sys
from time import time
import numpy as np
from numpy.linalg import norm
import scipy
import scipy.linalg, scipy.sparse

try:
    from testmatrix import testmatrix  # powerlaw, suitesparse ...
except ImportError:
    testmatrix = None

__version__ = "2020-07-03 July  denis-bz-py t-online.de"
np.set_printoptions( threshold=10, edgeitems=5, linewidth=120,
        formatter = dict( float = lambda x: "%.3g" % x,  # float arrays %.3g
                        complexfloat = " {0.real:.3g} {0.imag:+.3g}j " .format ))
print( 80 * "▄" )
print( "versions: numpy %s  scipy %s  python %s  plaform %s " % (
    np.__version__, scipy.__version__ , sys.version.split()[0], platform.platform() ))

#...........................................................................
test = "randomsparse"
# test = "Norris-fv1"  # suitesparse https://sparse.tamu.edu/Norris
n = 5000
density = .001
seed = 0
dtype = np.float64
drivers = "np ev evd"  # evr evx"

# to change these params, run this.py  a=1  b=None  'c = expr' ... in sh or ipython --
for arg in sys.argv[1:]:
    exec( arg )

#...........................................................................
if test == "randomsparse" or testmatrix is None:
    Asparse = scipy.sparse.random( n, n, density=density, random_state=seed )  # uniform 0 - 1
else:
    Asparse = testmatrix( test, n=n )
Asparse = Asparse.astype( dtype )
# A += A.T  -- only lower tri used
A = Asparse.A
print( "%s A %s density %.3g %s: dense %.0f mb" % (
        test, A.shape, density, A.dtype, A.nbytes / 1e6 ))

for driver in drivers.split():
    t0 = time()

    if driver == "np":
        npevals = np.linalg.eigvalsh( A )
        npevals = np.sort( npevals )
        print( "np.linalg.eigvalsh: %.0f sec  %s" % (
                time() - t0, npevals ))
    else:
        scipyevals = scipy.linalg.eigvalsh( A, driver=driver )  # driver= ev evd evr evx
        npevals = np.sort( npevals )
        print( "scipy.linalg.eigvalsh: %.0f sec  driver %s  %s" % (
                time() - t0, driver, scipyevals ))

diff = npevals - scipyevals
print( "max |npevals - scipyevals|: %.2g" % norm( diff, ord=np.inf ))