scipy.linalg.eigh and numpy.linalg.eigh calculates different eigenvalues for a symmetric matrix !

[Biswas__Sounak]

Hi,

  I am using Intel Python 2.7.14. This is what I get on my loading the python prompt :

Python 2.7.14 |Intel Corporation| (default, May  4 2018, 04:27:35).

I have come across a surprising case, where the eigenvalues of a symmetric 500 X 500 matrix calculated using scipy.linalg.eigh differs from the ones calculated using numpy.linalg.eigh. Further, the eigenvalues calculated by the scipy.linalg.eigh routine seem to be wrong, and two eigenvectors (v[:,449] and v[:,451] have NaN entries. The eigenvalues calculated using the numpy.linalg.eigh routine matches the results of the the general scipy.linalg.eig routine as well. I would be grateful for any suggestions as to what might be going on.

I have attached the code and the array(badmat.npy) in a zip file.
Here is the code:
########################################

import numpy as np
import scipy.linalg as LA
import numpy.linalg as nLA

Hup=np.load("badmat.npy")
#Verify that Hup is symmetric,returns zero
print LA.norm(Hup-np.transpose(Hup))
print np.shape(Hup)

e,v=nLA.eigh(Hup)
e=np.real(e)
idx=np.argsort(e)
e=e[idx]
e2,v2=LA.eigh(Hup)
#this does not return zero
print LA.norm(e2-e)

[Oleksandr_P_Intel]

Thank you for providing the script and the dataset. Please provide output of conda list --explicit, as well as your processor type. If you are on Linux, these could be obtained with lscpu. Thank you, Oleksandr

[Biswas__Sounak]

Oleksandr P. (Intel) wrote:

Thank you for providing the script and the dataset.

Please provide output of conda list --explicit, as well as your processor type. If you are on Linux, these could be obtained with lscpu.

Thank you,
Oleksandr

Hi, 
  I do not have conda. My distribution is installed from the tarball l_python2_pu3_2018.3.039.tgz (downloaded from intel), if that helps.

Here is my lscpu output :

Architecture:        x86_64
CPU op-mode(s):      32-bit, 64-bit
Byte Order:          Little Endian
CPU(s):              8
On-line CPU(s) list: 0-7
Thread(s) per core:  1
Core(s) per socket:  4
Socket(s):           2
NUMA node(s):        2
Vendor ID:           GenuineIntel
CPU family:          6
Model:               45
Model name:          Intel(R) Xeon(R) CPU E5-2609 0 @ 2.40GHz
Stepping:            7
CPU MHz:             1196.826
CPU max MHz:         2400.0000
CPU min MHz:         1200.0000
BogoMIPS:            4788.16
Virtualization:      VT-x
L1d cache:           32K
L1i cache:           32K
L2 cache:            256K
L3 cache:            10240K
NUMA node0 CPU(s):   0-3
NUMA node1 CPU(s):   4-7
Flags:               fpu vme de pse tsc msr pae mce cx8 apic sep mtrr pge mca cmov pat pse36 clflush dts acpi mmx fxsr sse sse2 ss ht tm pbe syscall nx pdpe1gb rdtscp lm constant_tsc arch_perfmon pebs bts rep_good nopl xtopology nonstop_tsc cpuid aperfmperf pni pclmulqdq dtes64 monitor ds_cpl vmx smx est tm2 ssse3 cx16 xtpr pdcm pcid dca sse4_1 sse4_2 x2apic popcnt tsc_deadline_timer aes xsave avx lahf_lm epb pti ibrs ibpb stibp tpr_shadow vnmi flexpriority ept vpid xsaveopt dtherm arat pln pts

Thanks

[Oleksandr_P_Intel]

Thank your for providing that information. I was able to locate the installer, and recreate your python installation. 

Intel Distribution for Python is built on conda, so you can do `source ./intelpython2/bin/activate` to activate conda environment. 

After this `conda list` allows to query environment for version of installed packages and libraries, which was MKL 2018.0.3, running NumPy 1.14.3, and SciPy 1.0.1

Running your script on `Intel(R) Xeon(R) CPU E5-2630 v3 @ 2.40GHz` I see 
 

[15:17:28 vmlin03 eig_vs_eigh]$ ./intelpython2/bin/python2 test.py
0.0
(500, 500)
1.16100171441e-14

I used flags to try to reproduce your hardware behavior as well. Your processor supports sse4_2, and avx vector instructions.

Using SSE4.2 produces expected result:
 

[15:20:50 vmlin03 eig_vs_eigh]$ MKL_ENABLE_INSTRUCTIONS=SSE4_2 MKL_CBWR=SSE4_2 MKL_NUM_THREADS=1 ./intelpython2/bin/python2 test.py
0.0
(500, 500)
6.01311978086e-15

But using avx, there seem to be a fault:
 

[15:21:07 vmlin03 eig_vs_eigh]$ MKL_ENABLE_INSTRUCTIONS=AVX MKL_CBWR=AVX MKL_NUM_THREADS=1 ./intelpython2/bin/python2 test.py
0.0
(500, 500)
0.237262552795

The problem with AVX path is still reproducible in the latest version of MKL. I will file a bug on your behalf.

As a work-around, try using MKL_ENABLE_INSTRUCTIONS=SSE4_2

[Biswas__Sounak]

Dear Oleksandr,

   Thanks  a lot for taking a look at this. I will use SSE4.2 in all cases from now on.
-Sounak.

[Oleksandr_P_Intel]

It turns out that if eigenvalues produced by call to scipy.linalg.eigh are also sorted

 

e2,v2=LA.eigh(Hup)

# idx2=np.argsort(e2)
# e2 = e[idx2]

Then the discrepancy goes away even when using AVX instruction set.

[Pamela_H_Intel]

Sounak,

Thank you for bringing this issue to our attention. Oleksandr was able to reproduce it. The issue is quite specific and only appears for CBWR mode + AVX codepath + sequential mode, and only for very specific input matrices. This is fixed in MKL 2020 update 1 which releases in a few weeks.

Pamela