使用numpy找不到矩阵的正确反值

Question

我试图使用numpy来求解矩阵逆，但无法使其工作。Python代码如下所示：

import copy
import numpy as np
from numpy.linalg import *
K=[[0.13535533905932737, -0.03535533905932737, 0.0, 0.0, -0.1],       
   [-0.03535533905932737, 0.13535533905932737, 0.0, -0.1, 0.0], 
   [0.0, 0.0, 0.13535533905932737, 0.03535533905932737, -0.03535533905932737], 
   [0.0, -0.1, 0.03535533905932737, 0.13535533905932737, -0.03535533905932737],
   [-0.1, 0.0, -0.03535533905932737, -0.03535533905932737, 0.13535533905932737]]
newK = np.array(K)
newK = inv(newK)

结果显示为：

[[ -2.74615429e+16  -2.74615429e+16   1.89368970e+00  -2.74615429e+16
   -2.74615429e+16]
 [ -2.74615429e+16  -2.74615429e+16  -1.37357926e+00  -2.74615429e+16
   -2.74615429e+16]
 [  9.23495156e-01  -1.84699031e+00   8.84484598e+00  -3.69398063e+00
    1.84699031e+00]
 [ -2.74615429e+16  -2.74615429e+16  -2.52873329e+00  -2.74615429e+16
   -2.74615429e+16]
 [ -2.74615429e+16  -2.74615429e+16   3.04884372e+00  -2.74615429e+16
   -2.74615429e+16]]

我也尝试使用MATLAB来检查这个值，但结果是：

    a =
    0.1340   -0.0350         0         0   -0.1000
   -0.0350    0.1350         0   -0.1000         0
         0         0    0.1350    0.0350   -0.0350
         0   -0.1000    0.0350    0.1350   -0.0350
   -0.1000         0   -0.0350   -0.0350    0.1350
>> inv(a)
ans =
   1.0e+03 *
   -1.0000   -1.0000    0.0000   -1.0000   -1.0000
   -1.0000   -0.9808   -0.0033   -0.9841   -0.9967
    0.0000   -0.0033    0.0089   -0.0044    0.0011
   -1.0000   -0.9841   -0.0044   -0.9785   -0.9956
   -1.0000   -0.9967    0.0011   -0.9956   -0.9911

有人能帮我做这个吗？

Answer 1

首先，两个程序的两个矩阵不相同：

numpy:

In [295]: np.linalg.det(newK)
Out[295]: -4.5051744884111939e-21

matlab：

>> det(a)
ans =
    -1.067499999999957e-07

在numpy:中使用相同的矩阵

In [296]: K2=np.array([[0.1340, -0.0350, 0, 0, -0.1000],[-0.0350, 0.1350, 0, -0.1000, 0],[0, 0, 0.1350, 0.0350, -0.0350],[0, -0.1000, 0.0350, 0.1350, -0.0350],[-0.1000, 0, -0.0350, -0.0350, 0.135]])
In [297]: np.linalg.det(K2)
Out[297]: -1.0674999999999714e-07

相反的结果是正确的：

In [298]: np.linalg.inv(K2)
Out[298]: 
array([[ -1.00000000e+03,  -1.00000000e+03,  -2.93090593e-14,
         -1.00000000e+03,  -1.00000000e+03],
       [ -1.00000000e+03,  -9.80796253e+02,  -3.27868852e+00,
         -9.84074941e+02,  -9.96721311e+02],
       [ -5.26327952e-14,  -3.27868852e+00,   8.85245902e+00,
         -4.42622951e+00,   1.14754098e+00],
       [ -1.00000000e+03,  -9.84074941e+02,  -4.42622951e+00,
         -9.78501171e+02,  -9.95573770e+02],
       [ -1.00000000e+03,  -9.96721311e+02,   1.14754098e+00,
         -9.95573770e+02,  -9.91147541e+02]])

其次，在numpy情况下，矩阵实际上是奇异的，它的行列式是1e-21。这意味着它的逆不存在，并且在结果中充其量是定义不清的（解释大小为1e16的矩阵元素）。

请注意，矩阵之间的差异无法手工消除，即"内部两个矩阵是相同的，只有matlab没有写那么多小数位数"，因为numpy 0.13535533905932737的矩阵元素在matlab中打印为0.1340，无论打印精度如何，这显然是错误的。

还要注意，你的matlab矩阵不是奇异的唯一原因是它的大多数元素是0.135，但它的第一个元素是0.134。如果你让这个元素与其他元素相等（就像在numpy的情况下一样！），你就会得到

>> b=a
>> b(1,1)=0.135
>> det(b)
ans =
     1.481453848484194e-21

你有了它。如果你的实际矩阵是numpy中的那个，那么解决方案很简单：不要反转它，因为它没有逆。

---

为了说服你，考虑一个与newK矩阵具有相同元素模式的随机矩阵：

k1=np.random.rand()*2-1 #element [0,0]
k2=np.random.rand()*2-1 #element [0,1]
k3=k1+k2 #element [0,-1]
Krand=np.array([[k1,k2,0,0,k3],[k2,k1,0,k3,0],[0,0,k1,-k2,k2],[0,k3,-k2,k1,k2],[k3,0,k2,k2,k1]])

然后你得到

In [322]: np.linalg.det(Krand)
Out[322]: 3.3121334687895703e-17

意味着任何具有这种结构的矩阵都是奇异的。

使用sympy的最终证明：

import sympy as sym
k1,k2=sym.symbols('k1,k2')
Ksym=sym.Matrix([[k1,k2,0,0,k1+k2],[k2,k1,0,k1+k2,0],[0,0,k1,-k2,k2],[0,k1+k2,-k2,k1,k2],[k1+k2,0,k2,k2,k1]])

那么符号行列式就是

In [347]: sym.det(Ksym)
Out[347]: 0