numpy计算的特征值不正确

Question

我正在尝试使用numpy在Python中获得正半定矩阵的特征值。

所有的特征值都应该是非负的和实的。最小特征值应为零。有时，numpy.eig返回复数值作为特征值，即使它们被认为是实数(尽管复数在虚部中有0j(。对于相关SO问题，接受的答案建议使用numpy.eigh而不是numpy.eig。我尝试了numpy.eigh，但它产生了不正确的结果：在下面的MEW的第二种情况下，最小特征值不为零。

我知道numpy.eigh是针对复埃尔米特矩阵或实对称矩阵的。下面的第二种情况不是这些。但似乎numpy.eig已知会返回复数特征值，即使在不应该返回的情况下也是如此(见相关问题(！一个人应该怎么做？

MWE(见第二种情况的输出(：

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

这是输出：

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

有人能帮我做这个吗？

---

EDIT：以前，MWE中的图形是定向的。我用无向图编辑了这个问题。

Answer 1

eig在此处执行正确。

您看到的是浮点精度的局限性。

在使用64位ieee753浮点数的计算机上，当涉及到数值计算结果时，+-1e-16为0。您只能获得15到16位小数精度。

因此，如果你有外部知识，你的矩阵是半正定的，因此应该有正的实本征值，你可以四舍五入这些值。

Numpy对复杂值有一个方便的函数，对负值也很容易做到。

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)

输出：

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