numpy.linalg.inv()是否给出了正确的矩阵逆？编辑：为什么inv()会出现数值错误？

Question

我有一个形状为（4000，4000）的矩阵，我想取它的逆。（对于如此大的矩阵，我的矩阵求逆的直觉会崩溃。）起始矩阵具有数量级为e-10的值，以下是一些值：print matrix给出了输出。

[[  2.19885119e-10   2.16462810e-10   2.13062782e-10 ...,  -2.16462810e-10
   -2.19885119e-10  -2.16462810e-10]
 [  2.16462810e-10   2.19885119e-10   2.16462810e-10 ...,  -2.13062782e-10
   -2.16462810e-10  -2.19885119e-10]
 [  2.13062782e-10   2.16462810e-10   2.19885119e-10 ...,  -2.16462810e-10
   -2.13062782e-10  -2.16462810e-10]
 ..., 
 [ -2.16462810e-10  -2.13062782e-10  -2.16462810e-10 ...,   2.19885119e-10
    2.16462810e-10   2.13062782e-10]
 [ -2.19885119e-10  -2.16462810e-10  -2.13062782e-10 ...,   2.16462810e-10
    2.19885119e-10   2.16462810e-10]
 [ -2.16462810e-10  -2.19885119e-10  -2.16462810e-10 ...,   2.13062782e-10
    2.16462810e-10   2.19885119e-10]]

我使用NumPy的numpy.linalg.inv()函数来求矩阵的逆。

import numpy as np
new_matrix = np.linalg.inv(matrix)
print new_matrix

以下是我得到的输出结果：

这是我得到的返回结果：

[[  1.95176541e+25   9.66643852e+23  -1.22660930e+25 ...,  -1.96621184e+25
   -9.41413909e+24   1.33500310e+25]
 [  2.01500967e+25   1.08946558e+24  -1.25813014e+25 ...,  -2.07717912e+25
   -9.86804459e+24   1.42950556e+25]
 [  3.55575106e+25   2.11333704e+24  -2.25333936e+25 ...,  -3.68616202e+25
   -1.72651875e+25   2.51239524e+25]
 ..., 
 [  3.07255588e+25   1.61759838e+24  -1.95678425e+25 ...,  -3.15440712e+25
   -1.47472306e+25   2.13570651e+25]
 [ -7.24380790e+24  -8.63730581e+23   4.90519245e+24 ...,   8.30663797e+24
    3.70858694e+24  -5.32291734e+24]
 [ -1.95760004e+25  -1.12341031e+24   1.23820305e+25 ...,   2.01608416e+25
    9.40221886e+24  -1.37605863e+25]]

这差别也太大了吧！怎么会这样？一个幅度为e-10的矩阵被反转成了一个幅度为e+25的矩阵？

这是否数学上正确，还是IEEE浮点值出现问题了？

如果这是数学上正确的，有人能解释一下背后的数学原理吗？

编辑：

根据下面的评论，我决定进行测试。

np.dot(matrix, new_matrix)应该给出单位矩阵，A * A^T = Identity。

这是我的输出结果：

[[  0.   -3.  -16.  ...,  16.    8.   12. ]
 [-24.   -1.5  -8.  ...,  32.   -4.   36. ]
 [ 40.    1.  -64.  ...,  24.   20.   24. ]
 ..., 
 [ 32.   -0.5  48.  ..., -16.  -20.   16. ]
 [ 40.    7.   16.  ..., -48.  -36.  -28. ]
 [ 16.    3.   12.  ..., -80.   16.    0. ]]

为什么numpy.linalg.inv()会产生数值误差？

np.allclose( np.dot(matrix, new_matrix), np.identity(4000) )

返回 False。

- ShanZhengYang

---

你的矩阵的行列式(numpy.linalg.det(m))和条件数(numpy.linalg.cond(m))是多少？ - DSM

如果矩阵的行列式为零，则它没有逆矩阵。您真的确定 A 的零行列式是不正确的吗？我尝试了一些随机数矩阵大小为 400x400 的测试，它可以工作，只因浮点运算而有轻微误差。 - user3371637

@Sauruxum numpy.linalg.det() 经常给出错误的答案。该操作首先执行 slogdet()，然后取指数。请参见代码。 此外，矩阵确实具有逆矩阵。没有任何理由不具备。 - ShanZhengYang

@RobFoley 不是的，你怎么会这么想呢？Python和numpy都使用标准的IEEE浮点数。 - ali_m

你真的需要求逆矩阵吗？如果你想解一个线性方程组，使用np.linalg.solve通过分解计算解比直接求逆矩阵更快更准确。参考链接：http://www.johndcook.com/blog/2010/01/19/dont-invert-that-matrix/ - ali_m

Answer 1

您的矩阵病态，因为

np.linalg.cond(matrix) > np.finfo(matrix.dtype).eps

根据这个答案，您可以考虑使用奇异值分解来求解这样的矩阵的逆。

- rth

---

应该根据数组的数据类型（即np.finfo(matrix).eps）来确定epsilon，这可能与sys.float_info.epsilon不同。 - ali_m

抱歉，应该是 np.finfo(matrix.dtype).eps。 - ali_m

Answer 2

对于这2个矩阵的行列式，你有以下结果

det(A) * det(A^{-1}) = 1

如果det(A)很大，那么det(A^{-1})就很小。对于这2个矩阵的范数（如果您选择了次乘性范数），您有：

1  =  |A*A^{-1}| >= |A| |A^-1|

其中||是一个次乘性的合理范数选择。在这里，你可以直观地看到你所观察到的数字：如果>=符号实际上是~=，那么你会得到相同的观察结果，这对于行列式来说是严格正确的。

如果考虑乘积，则同样适用相同的推理。

A * A^{-1} = 1

对于所有元素为正数的矩阵A。对于右手边的对角线上的元素，如果A的元素很大，则需要从A^{-1}中取非常小的数字。

PS：注意，这并不证明这种趋势一直存在。这只是提供了数学直觉，解释为什么会观察到这种缩放。

编辑，回复评论:

最初的问题是“如果这在数学上是正确的，能否有人向我解释一下其数学直觉？”。事实上，对于一个具有小数字的矩阵，逆将具有大数字是数学上正确和可靠的。上面我解释了为什么会出现这种情况。

回答OP编辑时提出的另一个问题，即为什么inv()会导致数值误差：求逆矩阵是一个困难的问题。这就是为什么我们每次尽可能地避免求逆矩阵的原因。例如，对于以下问题：

A x = b

我们不会计算 A 的逆，而是使用其他算法（在实践中，例如在 Python 中调用 scipy.linalg.solve）。

- gg349

---

如果这是正确的，我的测试应该是取乘积，对吗？ 原始矩阵是 matrix。逆矩阵是 new_matrix。因此，np.dot(matrix, new_matrix) 应该给我们单位矩阵，对吗？ - ShanZhengYang

你不应该对这个大小的矩阵的数值误差如此确定：使用通过 A=rand(400,400) 生成的随机矩阵计算 res = dot(A,inv(A))，其非对角线元素绝对值小于 1e-11（当然这取决于矩阵本身，而不仅仅是大小）。 - gg349

请看这里的预处理器是什么。根据您要解决的问题，我建议不要求矩阵的逆。 - gg349

@rth "你不能使用scipy.linalg.solve来求矩阵的逆" - 当然可以！事实上，这就是np.linalg.inv的工作原理（它执行的是np.linalg.solve(A, np.eye(A.shape[0]))的等价操作）。 - ali_m

至少在当前的numpy开发版本中，“inv”也使用“gesv”。请参见此处 - “inv”中使用的“lapack_func”是“gesv” - 它只是将一个单位矩阵作为“B”参数传递。 - ali_m