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.
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]