如何计算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]]

这些结果（如果规范化）与计算它们时得到的结果很接近，但在特征值和特征向量中似乎都引入了一些错误。有没有办法用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并将其设置为显示乳胶输出，则将get a nicer display。

Answer 2

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

因为你知道矩阵是对称的，所以用八可以做得更好。 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