Negative eigenvalues in call to dsyevr

Coming from the as a reference: How to include LAPACK functions with complex inputs · Issue #5301 · numba/numba · GitHub

I have the following attempt at calling dsyevr from scipy using numba:

import ctypes
import numpy as np
import numba as nb
from numba import vectorize, njit, jit
from numba.extending import get_cython_function_address

addr = get_cython_function_address('scipy.linalg.cython_lapack','dsyevr')
_PTR, _dble, _char, _int  = ctypes.POINTER, ctypes.c_double, ctypes.c_char, ctypes.c_int
_ptr_select = ctypes.c_voidp
_ptr_dble = _PTR(_dble)
_ptr_char = _PTR(_char)
_ptr_int  = _PTR(_int)

## Retrieve function address
functype = ctypes.CFUNCTYPE(None,
	_ptr_int, # JOBVS
	_ptr_int, # RANGE
	_ptr_int, # UPLO
	_ptr_int, # N
	_ptr_dble, # A
	_ptr_int, # LDA
	_ptr_dble, # VL
	_ptr_dble, # VU
	_ptr_int, # IL
	_ptr_int, # IU
	_ptr_dble, # ABSTOL
	_ptr_int, # M
	_ptr_dble, # W
	_ptr_dble, # Z
	_ptr_int, # LDZ
	_ptr_int, # ISUPPZ 
	_ptr_dble, # WORK
	_ptr_int, # LWORK
	_ptr_int, # IWORK 
	_ptr_int, # LIWORK 
	_ptr_int, # INFO
)
dsyevr_fn = functype(addr)

_ORD_JOBVS = ord('N')
_ORD_RNG = ord('I')
_ORD_UPLO = ord('U')

@nb.jit(nb.types.Tuple((nb.float64[:], nb.float64[:,:], nb.int32))(nb.float64[:,:], nb.int32, nb.int32), nopython=True)
def numba_dsyevr(x, I1, I2):
	assert (x.shape[0] == x.shape[1]) and x[0,0] == 0.0
	assert I1 <= I2 and I1 >= 1 and I2 <= x.shape[0]
	
	## Setup arguments
	i1, i2 = I1, I2
	n, m = x.shape[0], abs(i2-i1)+1
	JOBVS = np.array([_ORD_JOBVS], np.int32)
	RNG = np.array([_ORD_RNG], np.int32)
	UPLO = np.array([_ORD_UPLO], np.int32)
	N = np.array([n], np.int32)
	A = x.copy()     # in & out
	LDA = np.array([n], np.int32)
	VL, VU, IL, IU =  np.array(0, np.float64), np.array(0, np.float64), np.array(i1, np.int32), np.array(i2, np.int32)
	ABS_TOL = np.array([0.0], np.float64)
	M = np.array([m], np.int32)
	W, Z, LDZ, ISUPPZ, WORK, LWORK = np.zeros(n, np.float64), np.zeros((n, m), np.float64).T, np.array([n], np.int32), np.zeros(2*m, np.int32), np.zeros(26*n, np.float64), np.array([26*n], np.int32) 
	IWORK, LIWORK = np.zeros(10*n, np.int32), np.array([10*n], np.int32) 
	INFO = np.empty(1, dtype=np.int32)
	
	## Call the function
	# jobvs = ctypes.POINTER(ctypes.c_int)()
	dsyevr_fn(
		JOBVS.ctypes,
		RNG.ctypes,
		UPLO.ctypes,
		N.ctypes,
		A.view(np.float64).ctypes,
		LDA.ctypes,
		VL.ctypes,
		VU.ctypes,
		IL.ctypes,
		IU.ctypes,
		ABS_TOL.ctypes,
		M.ctypes,
		W.view(np.float64).ctypes,
		Z.view(np.float64).ctypes,
		LDZ.ctypes,
		ISUPPZ.view(np.int32).ctypes,
		WORK.view(np.float64).ctypes,
		LWORK.ctypes,
		IWORK.view(np.int32).ctypes,
		LIWORK.ctypes,
		INFO.ctypes
	)
	return((W, Z.T, INFO[0]))

I’m trying to do this to improve an MDS calculation. The function compiles, returns INFO=0, and I have the tolerance down to 0.0 (which the docs say you can do). I’ve tripled checked the types, the sizes of the inputs, etc… nonetheless, I’m getting negative eigenvalues w/ the following:

import numpy as np 
from scipy.spatial.distance import squareform, pdist, cdist

a = np.random.uniform(size=(15,8))
x = squareform(pdist(a))**2
z = numba_dsyevr(x, 1, x.shape[0]-1)

Does anyone have any idea why this is happening or possible things I can do to debug this? With INFO returning 0 I’m not sure how to approach

Hi @peekxc

I’ve not looked at whether the binding is valid but have run the following against the above:

import numpy as np
from scipy.spatial.distance import squareform, pdist, cdist

np.random.seed(0) # for reproducibility
a = np.random.uniform(size=(15, 8))
x = squareform(pdist(a))**2

# for stable comparisons
def sort(a):
    return np.sort(np.real(np.abs(a)))[::-1]

numba_result = sort(numba_dsyevr(x.copy(), 1, x.shape[0])[0])
numpy_result = sort(np.linalg.eigvals(x))
u, s, vt = np.linalg.svd(x)

print("NumPy", numpy_result)
print("Numba", numba_result)
print("diff", np.abs(numpy_result - numba_result))
print("condition", np.linalg.cond(x))
print("singular values", sort(s))

which gives me:

NumPy [1.94882095e+01 5.98132395e+00 4.86940486e+00 3.11813389e+00
 2.41255451e+00 1.08677471e+00 7.81473597e-01 6.81483528e-01
 3.42723564e-01 2.14336884e-01 6.57949205e-16 5.17434069e-16
 2.08111161e-16 2.08111161e-16 9.14223340e-17]
Numba [1.94882095e+01 5.98132395e+00 4.86940486e+00 3.11813389e+00
 2.41255451e+00 1.08677471e+00 7.81473597e-01 6.81483528e-01
 3.42723564e-01 2.14336884e-01 7.47479613e-15 9.02141586e-16
 3.82164695e-16 1.24155953e-16 2.74129994e-17]
diff [3.19744231e-14 4.44089210e-15 8.88178420e-16 3.10862447e-15
 1.33226763e-15 1.77635684e-15 1.11022302e-16 0.00000000e+00
 6.10622664e-16 7.49400542e-16 6.81684692e-15 3.84707517e-16
 1.74053534e-16 8.39552081e-17 6.40093346e-17]
condition 1.4356530015086147e+18
singular values [1.94882095e+01 5.98132395e+00 4.86940486e+00 3.11813389e+00
 2.41255451e+00 1.08677471e+00 7.81473597e-01 6.81483528e-01
 3.42723564e-01 2.14336884e-01 1.30265001e-15 5.22229917e-16
 2.47226779e-16 1.96053786e-16 1.35744567e-17]

i.e. the result matches NumPy reasonably given the system is fairly poorly conditioned. Singular values also match.

Also, I think a real symmetric matrix has real eigenvalues, they don’t have to be positive? For example:

| 0 2|
| 2 0|

has eigenvalues (-2, 2).

Ahh, I misremembered the properties of a real symmetric matrix. Indeed, I had a distance matrix which was not being re-centered (thus, it had all zeros along the diagonal). I apologize for thinking it was a software fault.

In case anybody comes across this wondering how to bind Fortran LAPACK routines w/ numba, here is my implementation which wraps the dsyevr for calculation of a pre-specified range of eigenvalues/vectors.

Note that interactive sessions require the *.ctypes parameters to be converted to their appropriate data instances if one is not JIT/AOT compiling in a file.

@peekxc Great, glad you got it working as expected.