如何计算 numpy 特征值和特征向量

Question

我正在使用 numpy 计算对称方形数组的特征值和特征向量。我的数组是：

L = [[ 2. -1. -1.  0.  0.  0.]
     [-1.  3.  0. -1.  0. -1.]
     [-1.  0.  2. -1.  0.  0.]
     [ 0. -1. -1.  3. -1.  0.]
     [ 0.  0.  0. -1.  2. -1.]
     [ 0. -1.  0.  0. -1.  2.]]

执行numpy.linalg.eig(L)时的结果如下所示

特征值：

[ 5.00000000e+00,   
  3.96872205e-16,   
  1.00000000e+00,
  2.00000000e+00,   
  3.00000000e+00,   
  3.00000000e+00 ]

特征向量：

[[ -2.88675135e-01   4.08248290e-01  -5.00000000e-01  4.08248290e-01   -4.36632863e-01   4.44614891e-01]
 [  5.77350269e-01   4.08248290e-01  -3.34129212e-16  4.08248290e-01   -1.08813217e-01  -5.41271705e-01]
 [  2.88675135e-01   4.08248290e-01  -5.00000000e-01  4.08248290e-01    5.45446080e-01   9.66568140e-02]
 [ -5.77350269e-01   4.08248290e-01   1.06732810e-16  4.08248290e-01   -1.08813217e-01  -5.41271705e-01]
 [  2.88675135e-01   4.08248290e-01   5.00000000e-01  4.08248290e-01   -4.36632863e-01   4.44614891e-01]
 [ -2.88675135e-01   4.08248290e-01   5.00000000e-01 -4.08248290e-01    5.45446080e-01   9.66568140e-02]]

结果与您在analytically 计算时得到的结果接近（如果标准化），但似乎在特征值和特征向量中都引入了一些错误。有没有办法使用 numpy 绕过这些错误？

这些错误来自哪里？ numpy 使用什么算法？

Answer 1

如果您想要解析推导的精度，您将需要使用 symbolic computation，这是 Wolfram Alpha、Mathematica 和相关系统使用的。例如，在 Python 中，您可能需要查看 SymPy。

嵌入到您正在使用的 NumPy 包中的 numerical computation 本质上会受到 floating point numerical representations 的小错误和变迁的影响。这种误差和近似值在数值计算中是不可避免的。

这是一个例子：

from sympy import Matrix, pretty

L = Matrix([[ 2, -1, -1,  0,  0,  0,],
     [-1,  3,  0, -1,  0, -1,],
     [-1,  0,  2, -1,  0,  0,],
     [ 0, -1, -1,  3, -1,  0,],
     [ 0,  0,  0, -1,  2, -1,],
     [ 0, -1,  0,  0, -1,  2,]])

print "eigenvalues:"
print pretty(L.eigenvals())
print
print "eigenvectors:"
print pretty(L.eigenvects(), num_columns=132)

产量：

eigenvalues:
{0: 1, 1: 1, 2: 1, 3: 2, 5: 1}

eigenvectors:
⎡⎛0, 1, ⎡⎡1⎤⎤⎞, ⎛1, 1, ⎡⎡-1⎤⎤⎞, ⎛2, 1, ⎡⎡1 ⎤⎤⎞, ⎛3, 2, ⎡⎡1 ⎤, ⎡0 ⎤⎤⎞, ⎛5, 1, ⎡⎡1 ⎤⎤⎞⎤
⎢⎜      ⎢⎢ ⎥⎥⎟  ⎜      ⎢⎢  ⎥⎥⎟  ⎜      ⎢⎢  ⎥⎥⎟  ⎜      ⎢⎢  ⎥  ⎢  ⎥⎥⎟  ⎜      ⎢⎢  ⎥⎥⎟⎥
⎢⎜      ⎢⎢1⎥⎥⎟  ⎜      ⎢⎢0 ⎥⎥⎟  ⎜      ⎢⎢1 ⎥⎥⎟  ⎜      ⎢⎢-1⎥  ⎢-1⎥⎥⎟  ⎜      ⎢⎢-2⎥⎥⎟⎥
⎢⎜      ⎢⎢ ⎥⎥⎟  ⎜      ⎢⎢  ⎥⎥⎟  ⎜      ⎢⎢  ⎥⎥⎟  ⎜      ⎢⎢  ⎥  ⎢  ⎥⎥⎟  ⎜      ⎢⎢  ⎥⎥⎟⎥
⎢⎜      ⎢⎢1⎥⎥⎟  ⎜      ⎢⎢-1⎥⎥⎟  ⎜      ⎢⎢-1⎥⎥⎟  ⎜      ⎢⎢0 ⎥  ⎢1 ⎥⎥⎟  ⎜      ⎢⎢-1⎥⎥⎟⎥
⎢⎜      ⎢⎢ ⎥⎥⎟  ⎜      ⎢⎢  ⎥⎥⎟  ⎜      ⎢⎢  ⎥⎥⎟  ⎜      ⎢⎢  ⎥  ⎢  ⎥⎥⎟  ⎜      ⎢⎢  ⎥⎥⎟⎥
⎢⎜      ⎢⎢1⎥⎥⎟  ⎜      ⎢⎢0 ⎥⎥⎟  ⎜      ⎢⎢-1⎥⎥⎟  ⎜      ⎢⎢-1⎥  ⎢-1⎥⎥⎟  ⎜      ⎢⎢2 ⎥⎥⎟⎥
⎢⎜      ⎢⎢ ⎥⎥⎟  ⎜      ⎢⎢  ⎥⎥⎟  ⎜      ⎢⎢  ⎥⎥⎟  ⎜      ⎢⎢  ⎥  ⎢  ⎥⎥⎟  ⎜      ⎢⎢  ⎥⎥⎟⎥
⎢⎜      ⎢⎢1⎥⎥⎟  ⎜      ⎢⎢1 ⎥⎥⎟  ⎜      ⎢⎢-1⎥⎥⎟  ⎜      ⎢⎢1 ⎥  ⎢0 ⎥⎥⎟  ⎜      ⎢⎢-1⎥⎥⎟⎥
⎢⎜      ⎢⎢ ⎥⎥⎟  ⎜      ⎢⎢  ⎥⎥⎟  ⎜      ⎢⎢  ⎥⎥⎟  ⎜      ⎢⎢  ⎥  ⎢  ⎥⎥⎟  ⎜      ⎢⎢  ⎥⎥⎟⎥
⎣⎝      ⎣⎣1⎦⎦⎠  ⎝      ⎣⎣1 ⎦⎦⎠  ⎝      ⎣⎣1 ⎦⎦⎠  ⎝      ⎣⎣0 ⎦  ⎣1 ⎦⎦⎠  ⎝      ⎣⎣1 ⎦⎦⎠⎦

虽然 ASCII 漂亮的打印机，嗯，努力工作以提供看起来很漂亮的输出，但您可以看到您正在获得符号计算的精确输出。如果你使用 IPython 并设置它显示 LaTeX 输出，你会get a nicer display。

Answer 2

看起来它正在使用 LAPACK 的迭代方法。它收敛到一个解决方案。如果不收敛，则抛出异常。

既然你知道矩阵是对称的，你可能会用 eigh 做得更好。 http://docs.scipy.org/doc/numpy/reference/generated/numpy.linalg.eigh.html

文档页面： http://docs.scipy.org/doc/numpy/reference/generated/numpy.linalg.eig.html

源代码：https://github.com/numpy/numpy/blob/v1.9.1/numpy/linalg/linalg.py#L982