numpy computed eigenvalues incorrect?

Question

I am trying to get the eigenvalues of a positive semi-definite matrix in Python using numpy.

All eigenvalues should be non-negative and real. The minimum eigenvalue should be zero. Sometimes, numpy.eig returns complex values as eigenvalues even when they're supposed to be real (although the complex numbers have 0j in the imaginary part). To a related SO question, the accepted answer suggests using numpy.eigh instead of numpy.eig. I tried numpy.eigh but it produces incorrect results: in the second case in the MEW below, the smallest eigenvalue is not zero.

I know that numpy.eigh is for a complex Hermitian or a real symmetric matrix. The second case below is none of these. But it seems numpy.eig is known to return complex eigenvalues even in cases where it shouldn't (see linked question)! What should one do?

MWE (see output for second case):

import numpy as np
from numpy.random import RandomState

rng = RandomState(123456)

num_v = 4

weights = rng.uniform(0, 1, num_v * num_v).reshape((num_v, num_v))
weights = np.tril(weights) + np.triu(weights.T, 1)
np.fill_diagonal(weights, 0)

degrees = weights.sum(axis=0)

degrees_rw = 1 / degrees
D_rw = np.diag(degrees_rw)

degrees = 1 / np.sqrt(degrees)
D = np.diag(degrees)

I = np.eye(num_v)

laplacian = I - np.matmul(np.dot(D, weights), D)
laplacian_rw = I - np.matmul(D_rw, weights)

eigs1, _ = np.linalg.eig(laplacian)
eigs2, _ = np.linalg.eigh(laplacian)

print(np.sort(eigs1)[:10], "\n", np.sort(eigs2)[:10])

eigs1, _ = np.linalg.eig(laplacian_rw)
eigs2, _ = np.linalg.eigh(laplacian_rw)

print("\n", np.sort(eigs1)[:10], "\n", np.sort(eigs2)[:10])

Here's the output:

[0.         1.01712657 1.40644541 1.57642803] 
 [-1.08192365e-16  1.01712657e+00  1.40644541e+00  1.57642803e+00]

 [-2.22044605e-16  1.01712657e+00  1.40644541e+00  1.57642803e+00] 
 [0.08710134 1.01779168 1.38159284 1.51351414]

Could someone please help me with this?

EDIT: previously, the graph in the MWE was directed. I edited the question with an undirected graph.

Did you intend `laplacian` to be a symmetric matrix? The array you give in your MWE is not symmetric. (You _can_ define "positive semi-definite" for non-symmetric matrices, but it's a bit of an odd thing to do.) — Mark Dickinson, Jul 29 '22 at 07:12
@MarkDickinson you're right, thanks for checking! I created a directed graph by mistake. When I remove directionality, the laplacian is symmetric, and the eigenvalues are also as expected. However, for a similar matrix (which is positive semi-definite but not symmetric), I still get different values with these two functions. I edited the question with these changes. — vpk, Jul 29 '22 at 07:30
Hmm. I don't think you can expect meaningful results from `numpy.linalg.eigh` if you give a non-symmetric (or non-Hermitian in the complex case) input, so I wouldn't expect to see the same results from `eig` and `eigh` in this case. I'm still a bit confused: are you wanting this to work for a non-symmetric input? Why would you expect the eigenvalues to all be real in that case? — Mark Dickinson, Jul 29 '22 at 07:37
@MarkDickinson please see edited question. Yes I need this for a non-symmetric input. I expect eigenvalues to be real and non-negative because the input is positive semi-definite. I understand that `numpy.eigh` is for a complex Hermitian or a real symmetric matrix, but it seems `numpy.eig` is known to return complex eigenvalues even in cases where it shouldn't (see linked question)! What should one do? — vpk, Jul 29 '22 at 07:47
"because the input is positive semi-definite" <- The theorem about nonnegative real eigenvalues depends on symmetry: you can't expect the eigenvalues to be real in general. For example, consider the 2-by-2 matrix `A = [[1, 1], [-1, 1]]`. This is positive semi-definite in the sense that `x^T A x >= 0` for any 2d vector `x`, but its eigenvalues have nonzero imaginary part. — Mark Dickinson, Jul 29 '22 at 08:48
@MarkDickinson In this case, I can prove that the matrix I am dealing with has non-negatnve real eigenvalues. — vpk, Jul 29 '22 at 17:00

MaxNoe · Accepted Answer · 2022-07-29T08:14:11.473

eig is performing correctly here.

What you see are the limitations of floating point precision.

For all intents and purposes on a computer using 64 bit ieee753 floating point numbers, +-1e-16 is 0 when it comes to the results of numeric computations. You only get 15 to 16 decimal digits precision.

So if you have the external knowledge, that your matrices are positive semidefinite and thus should have positive real eigen values, you can round those values away.

Numpy has a convenience function for the complex values and its easy to do for the negative ones.

import numpy as np

# this is the new, recommended numpy random API
rng = np.random.default_rng(123456)

num_v = 4

weights = rng.uniform(0, 1, num_v * num_v).reshape((num_v, num_v))
weights = np.tril(weights) + np.triu(weights.T, 1)
np.fill_diagonal(weights, 0)

degrees = weights.sum(axis=0)
degrees_rw = 1 / degrees
D_rw = np.diag(degrees_rw)

I = np.eye(num_v)

laplacian_rw = I - np.matmul(D_rw, weights)

eigvals_rw, eigvecs_rw = np.linalg.eig(laplacian_rw)
print(eigvals_rw)


# cast to real values if the imaginary parts are in floating precision to 0
eigvals_rw = np.real_if_close(eigvals_rw)

# round values that are at the floating point precision to exactly 0
# might need more than 2 epsilon, but in this example, it was enough
small = np.abs(eigvals_rw) < 2 * np.finfo(eigvals_rw.dtype).eps
eigvals_rw[small] = 0

print(eigvals_rw)

Output:

[-2.22044605e-16  9.98074713e-01  1.62494665e+00  1.37697864e+00]
[0.         0.99807471 1.62494665 1.37697864]

That's what I first thought too when reading this question. However, I think the OP is referring to the difference between `-2.22044605e-16` and `0.08710134` (two last lines of the output). — MB-F, Jul 29 '22 at 08:05
@MB-F The different outputs are expected, IMO because the second matrix is not symmetric. I was asking how do I get the correct answer using `numpy.eig` itself, since `numpy.eigh` is not always suitable (as in the second case). — vpk, Jul 29 '22 at 08:07
I added example code, using both `np.real_if_close` (which is not needed in this specific example) and the rounding of small values to 0. — MaxNoe, Jul 29 '22 at 08:15

numpy computed eigenvalues incorrect?

1 Answers1