numpy 避免在Python中进行随机矩阵乘法时的溢出

tf7tbtn2  于 2022-11-10  发布在  Python
关注(0)|答案(2)|浏览(184)

我必须将许多(大约700个)矩阵与一个随机元素(在下面,我使用的是盒分布)相乘:


# define parameters

μ=2.
σ=2.
L=700

# define random matrix

T=[None]*L

product=np.array([[1,0],[0,1]])
for i in range(L):
    m=np.random.uniform(μ-σ*3**(1/2), μ+σ*3**(1/2)) #box distribution
    T[i]=np.array([[-1,-m/2],[1,0]])
    product=product.dot(T[i])  #multiplying matrices

Det=abs(np.linalg.det(product))
print(Det)

对于这个μ和σ的选择,我获得了e^30+量级的量,但这个量应该收敛到0。我怎么知道?因为从分析上讲,它可以被证明等同于:

Y=[None]*L
product1=np.array([[1,0],[0,1]])
for i in range(L):
    m=np.random.uniform(μ-σ*(3**(1/2)), μ+σ*(3**(1/2))) #box distribution
    Y[i]=np.array([[-m/2,0],[1,0]])
    product1=product1.dot(Y[i])

l,v=np.linalg.eig(product1)

print(abs(l[1]))

这确实给出了e^-60。我认为这里存在溢出问题。我怎么才能修好它?

编辑:

预计两个印刷量相等,因为第一个是行列式的abs:

根据比奈定理(乘积的行列式是行列式的乘积):

第二个代码打印最大本征值的abs:

很容易看出,一个本征值为0,另一个等于

jyztefdp

jyztefdp1#

我将使用m而不是m/2来简化公式,但它不会改变任何事情。
前两个指标的乘积为

[-1  -m1 ] * [-1 -m2 ] = [1-m1  m2 ]
[ 1   0  ]   [ 1  0  ]   [-1   -m2 ]

如果你拿DET,那就是(-m2+m1*m2) + m2
您可以看到一种更大项(M2)的取消形式,从而产生统计上更小的m1*m2
经过两次乘法后,情况变得更糟了

[m1+m2-1  m1*m3-m3 ]
[ 1-m2       m3    ]

DET为(m1*m3+m2*m3-m3)+(m1*m2*m3-m2*m3-m1*m3+m3)
同样,这两个项的大小是m3
而结果是较小的m1*m2*m3
少数几个操作必然会导致灾难性的取消
(big+small)-big
big计算中的数值噪声大大超过了精确结果的大小。
这表明问题不能通过简单的缩放来缓解,矩阵因设计而不知何故出现了病态。
您可以尝试将随机数转换为有理数(浮点值是有理数),并用有理数(无限精度)计算乘积,但对于700项,预计会有非常大的整数和非常慢的计算速度……

s3fp2yjn

s3fp2yjn2#

这通常是一个棘手的问题。有几篇关于浮点算法和精度的好文章。这是一款著名的one
其中一个常见的技巧是使用比例变量。就像这样:

import numpy as np

# define parameters

μ=2.
σ=2.
L=700

# define random matrix

T=[None]*L

scale = 1
product=np.array([[1,0],[0,1]])
for i in range(L):
    m=np.random.uniform(μ-σ*3**(1/2), μ+σ*3**(1/2)) #box distribution
    T[i]=np.array([[-1,-m/2],[1,0]])
    product=np.matmul(product, T[i])  #multiplying matrices
    scale *= abs(product[0][0])
    product /= abs(product[0][0])

Det=abs(np.linalg.det(product*scale))
print(Det)

它让事情变得更好一些,但不幸的是,它没有帮助。
在这种情况下,您可以做的是将行列式相乘,而不是矩阵。就像这样:

import numpy as np

# define parameters

μ=2.
σ=2.
L=700

# define random matrix

T=[None]*L

scale = 1
product=1
for i in range(L):
    m=np.random.uniform(μ-σ*3**(1/2), μ+σ*3**(1/2)) #box distribution
    T[i]=np.array([[-1,-m/2],[1,0]])
    product *= np.linalg.det(T[i])  #multiplying matrices determinants

Det=abs(product)
print(Det)

输出:

1.1081329233309736e-61

因此,这就解决了问题。

相关问题