numpy 1 e-20小元素矩阵求逆的正确方法

xn1cxnb4  于 2023-10-19  发布在  其他
关注(0)|答案(1)|浏览(92)

我生成的矩阵是对称的,它的元素的大小取决于我处理的问题的长度。通常,矩阵是9 x9,但也可以使用6x6版本。
矩阵的一个例子是这样的

// A
2.8538093853e-21 7.1747731778e-22 7.1747731778e-22 7.9791996325e-22 7.9791996325e-22 3.0550862084e-22
7.1747731778e-22 2.8538093853e-21 7.1747731778e-22 7.9791996325e-22 3.0550862084e-22 7.9791996325e-22
7.1747731778e-22 7.1747731778e-22 2.8538093853e-21 3.0550862084e-22 7.9791996325e-22 7.9791996325e-22
7.9791996325e-22 7.9791996325e-22 3.0550862084e-22 7.1747731778e-22 3.0550862084e-22 3.0550862084e-22
7.9791996325e-22 3.0550862084e-22 7.9791996325e-22 3.0550862084e-22 7.1747731778e-22 3.0550862084e-22
3.0550862084e-22 7.9791996325e-22 7.9791996325e-22 3.0550862084e-22 3.0550862084e-22 7.1747731778e-22

它包含相当小的值(1 e-22的数量级),基于Numpy,它们的行列式和条件数如下所示

cond 18.552692061099496
det  8.466833864265427e-127

所以条件数不是很糟糕,但是行列式非常非常小。因此,我预计Numpy将很难将其反转,并根据其内部公差为矩阵求逆的数值方案提供一些随机解决方案。
作为一种非常简单的方法来改善这种情况,我只是将这个矩阵缩放了一个特定的值(这里,我取1.3972450422e+21)。

// A2
3.9874710151e+00 1.0024916252e+00 1.0024916252e+00 1.1148897128e+00 1.1148897128e+00 4.2687040583e-01
1.0024916252e+00 3.9874710151e+00 1.0024916252e+00 1.1148897128e+00 4.2687040583e-01 1.1148897128e+00
1.0024916252e+00 1.0024916252e+00 3.9874710151e+00 4.2687040583e-01 1.1148897128e+00 1.1148897128e+00
1.1148897128e+00 1.1148897128e+00 4.2687040583e-01 1.0024916252e+00 4.2687040583e-01 4.2687040583e-01
1.1148897128e+00 4.2687040583e-01 1.1148897128e+00 4.2687040583e-01 1.0024916252e+00 4.2687040583e-01
4.2687040583e-01 1.1148897128e+00 1.1148897128e+00 4.2687040583e-01 4.2687040583e-01 1.0024916252e+00

它的条件数和行列式是

cond 18.552692061099485
det  6.300231416456931

很明显,条件数和行列式看起来很正常。
但我现在不能得到的是,这两个案件的结果在质量上并没有多大区别。逆矩阵值非常相似,并且A*invA显示出与单位矩阵的类似偏差(~ 1 e-16)。

inv(A)
6.4337397752e+20 -2.7166291469e+18 -2.7166291469e+18 -5.6175766999e+20 -5.6175766999e+20 2.1049115695e+20
-2.7166291469e+18 6.4337397752e+20 -2.7166291469e+18 -5.6175766999e+20 2.1049115695e+20 -5.6175766999e+20
-2.7166291469e+18 -2.7166291469e+18 6.4337397752e+20 2.1049115695e+20 -5.6175766999e+20 -5.6175766999e+20
-5.6175766999e+20 -5.6175766999e+20 2.1049115695e+20 2.9200923433e+21 -4.3031775290e+20 -4.3031775290e+20
-5.6175766999e+20 2.1049115695e+20 -5.6175766999e+20 -4.3031775290e+20 2.9200923433e+21 -4.3031775290e+20
2.1049115695e+20 -5.6175766999e+20 -5.6175766999e+20 -4.3031775290e+20 -4.3031775290e+20 2.9200923433e+21

inv(A2)*fnorm
6.4337397752e+20 -2.7166291469e+18 -2.7166291469e+18 -5.6175766999e+20 -5.6175766999e+20 2.1049115695e+20
-2.7166291469e+18 6.4337397752e+20 -2.7166291469e+18 -5.6175766999e+20 2.1049115695e+20 -5.6175766999e+20
-2.7166291469e+18 -2.7166291469e+18 6.4337397752e+20 2.1049115695e+20 -5.6175766999e+20 -5.6175766999e+20
-5.6175766999e+20 -5.6175766999e+20 2.1049115695e+20 2.9200923433e+21 -4.3031775290e+20 -4.3031775290e+20
-5.6175766999e+20 2.1049115695e+20 -5.6175766999e+20 -4.3031775290e+20 2.9200923433e+21 -4.3031775290e+20
2.1049115695e+20 -5.6175766999e+20 -5.6175766999e+20 -4.3031775290e+20 -4.3031775290e+20 2.9200923433e+21
check A*invA
1.0000000000e+00 -1.9428902931e-16 -2.7755575616e-17 -5.5511151231e-17 5.5511151231e-17 -5.5511151231e-17
-2.7755575616e-17 1.0000000000e+00 -5.5511151231e-17 -2.7755575616e-17 -5.5511151231e-17 0.0000000000e+00
1.6653345369e-16 0.0000000000e+00 1.0000000000e+00 2.7755575616e-17 5.5511151231e-17 0.0000000000e+00 
1.1102230246e-16 3.3306690739e-16 0.0000000000e+00 1.0000000000e+00 1.6653345369e-16 5.5511151231e-17 
-2.2204460493e-16 0.0000000000e+00 2.2204460493e-16 1.1102230246e-16 1.0000000000e+00 1.1102230246e-16 
0.0000000000e+00 3.3306690739e-16 0.0000000000e+00 -2.7755575616e-17 0.0000000000e+00 1.0000000000e+00 

check A2*invA2
1.0000000000e+00 -2.7755575616e-17 -1.1102230246e-16 4.1633363423e-17 6.9388939039e-17 2.7755575616e-17 
-2.7755575616e-17 1.0000000000e+00 5.5511151231e-17 0.0000000000e+00 -2.7755575616e-17 0.0000000000e+00 
5.5511151231e-17 0.0000000000e+00 1.0000000000e+00 0.0000000000e+00 -1.1102230246e-16 -5.5511151231e-17 
-8.3266726847e-17 -1.1102230246e-16 0.0000000000e+00 1.0000000000e+00 2.7755575616e-17 5.5511151231e-17 
-1.1102230246e-16 0.0000000000e+00 2.2204460493e-16 -5.5511151231e-17 1.0000000000e+00 0.0000000000e+00 
-2.2204460493e-16 -7.2164496601e-16 -2.2204460493e-16 8.3266726847e-17 1.1102230246e-16 1.0000000000e+00

虽然1 e-16阶的差异可以被认为是微不足道的,但我很担心,因为他们有迹象。如果这种错误导致了逆矩阵元素符号的改变,那么在我的代码中,我将得到一个反向力。在这个矩阵求逆过程中,我能做些什么呢?
也许我的问题首先是不适定的,所以我基本上无能为力?
下面是我在这里使用的完整Python代码

import numpy as np

def printA(A, n):
    for i in range(n):
        for j in range(n):
            print("{:.10e}".format(A[i,j]), end=" ")

        print("")

A = np.array([[2.8538093853e-21, 7.1747731778e-22, 7.1747731778e-22, 7.9791996325e-22, 7.9791996325e-22, 3.0550862084e-22],
              [7.1747731778e-22, 2.8538093853e-21, 7.1747731778e-22, 7.9791996325e-22, 3.0550862084e-22, 7.9791996325e-22],
              [7.1747731778e-22, 7.1747731778e-22, 2.8538093853e-21, 3.0550862084e-22, 7.9791996325e-22, 7.9791996325e-22],
              [7.9791996325e-22, 7.9791996325e-22, 3.0550862084e-22, 7.1747731778e-22, 3.0550862084e-22, 3.0550862084e-22],
              [7.9791996325e-22, 3.0550862084e-22, 7.9791996325e-22, 3.0550862084e-22, 7.1747731778e-22, 3.0550862084e-22],
              [3.0550862084e-22, 7.9791996325e-22, 7.9791996325e-22, 3.0550862084e-22, 3.0550862084e-22, 7.1747731778e-22]])

print("A - cond / det")
printA(A,6)
print("cond {}".format(np.linalg.cond(A)))
print("det  {}".format(np.linalg.det(A)))
#fnorm = 2.7944900844667493e+20
fnorm = 2.7944900844667493e+21 / 2.
A2 = A*fnorm
print("\nA2 = fnorm * A - cond / det")
print("fnorm {:.10e}".format(fnorm))
printA(A2,6)
print("cond {}".format(np.linalg.cond(A2)))
print("det  {}".format(np.linalg.det(A2)))

print("\ninv(A)")
printA(np.linalg.inv(A),6)

print("\ninv(A2)*fnorm")
printA(np.linalg.inv(A2)*fnorm,6)

print("\ncheck A*invA")
printA(np.matmul(np.linalg.inv(A),A),6)

print("\ncheck A2*invA2")
printA(np.matmul(np.linalg.inv(A2),A2),6)
gkl3eglg

gkl3eglg1#

你的代码中有个错误。我无法重现这个问题:

>>> A = np.array([[2.8538093853e-21, 7.1747731778e-22, 7.1747731778e-22, 7.9791996325e-22, 7.9791996325e-22, 3.0550862084e-22],              [7.1747731778e-22, 2.8538093853e-21, 7.1747731778e-22, 7.9791996325e-22, 3.0550862084e-22, 7.9791996325e-22],
              [7.1747731778e-22, 7.1747731778e-22, 2.8538093853e-21, 3.0550862084e-22, 7.9791996325e-22, 7.9791996325e-22],              [7.9791996325e-22, 7.9791996325e-22, 3.0550862084e-22, 7.1747731778e-22, 3.0550862084e-22, 3.0550862084e-22],
              [7.9791996325e-22, 3.0550862084e-22, 7.9791996325e-22, 3.0550862084e-22, 7.1747731778e-22, 3.0550862084e-22],              [3.0550862084e-22, 7.9791996325e-22, 7.9791996325e-22, 3.0550862084e-22, 3.0550862084e-22, 7.1747731778e-22]])

>>> np.linalg.inv(A*1e22).round(4)

array([[ 0.0643, -0.0003, -0.0003, -0.0562, -0.0562,  0.021 ],
       [-0.0003,  0.0643, -0.0003, -0.0562,  0.021 , -0.0562],
       [-0.0003, -0.0003,  0.0643,  0.021 , -0.0562, -0.0562],
       [-0.0562, -0.0562,  0.021 ,  0.292 , -0.043 , -0.043 ],
       [-0.0562,  0.021 , -0.0562, -0.043 ,  0.292 , -0.043 ],
       [ 0.021 , -0.0562, -0.0562, -0.043 , -0.043 ,  0.292 ]])

相关问题