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 =单位。

这是我的输出：

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

Answer 1

你的矩阵是病态的，因为

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

根据this answer的说法，您可以考虑使用Singular Value Decomposition来求这类矩阵的逆。

Answer 2

对于2个矩阵的行列式，您有

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

因此，如果det(A)很大，那么det(A^{-1})就很小。对于2个矩阵的范数(如果选择一个子乘法范数)，您可以使用have

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

其中||是对次乘法范数的合理选择。在这里，您可以直观地从数值上观察到什么:如果>=符号实际上是一个~=，那么您可以恢复对行列式完全正确的相同的观察结果。

如果你考虑这个产品，同样的推理也适用。

A * A^{-1} = 1

对于具有所有正元素的矩阵A。对于1在RHS的对角线上的元素，如果A的元素非常大，则需要来自A^{-1}的非常小的数字。

PS:但是请注意，这并不能证明这种趋势总是存在的。这只是提供了观察这种缩放的数学直觉。

编辑，回复评论：

最初的问题是“如果这在数学上是正确的，有人能给我解释一下背后的数学直觉吗？”事实上，给出一个小数字的矩阵，反矩阵就会有大数字，这在数学上是正确的，也是合理的。上面我解释了为什么会这样。

为了回答在OP的编辑中出现的另一个问题，这就是为什么inv()会导致数值错误:求矩阵的逆是一个很难的问题。这就是为什么每次我们都会避免颠倒它们的原因。例如，对于问题

A x = b

我们不计算A的倒数，而是使用其他算法(实际上，您可以在python中调用scipy.linalg.solve )。