Python中随机矩阵乘积行列式计算的溢出问题及解决

我需要在Python中计算约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}$，因此我推测存在数值溢出问题，请问该如何解决？

补充说明：
这两个输出理论等价，第一个代码输出的是乘积矩阵行列式的绝对值，根据Binet定理（乘积的行列式等于行列式的乘积），其值为各矩阵行列式的乘积；第二个代码输出的是另一乘积矩阵最大特征值的绝对值，该矩阵一个特征值为0，另一个等于前者的行列式值。

问题根源是直接计算700个矩阵的乘积会导致数值溢出——矩阵元素的绝对值会指数级增长，超出浮点数的表示范围，最终导致行列式计算失真。结合Binet定理，我们可以通过对数求和代替直接乘积来避免溢出，具体步骤如下：

1. 计算每个单个矩阵$T[i]$的行列式，取其绝对值的自然对数
2. 将所有对数结果求和，得到总对数
3. 对总对数取指数，得到最终行列式的绝对值

修改后的代码如下：

import numpy as np

# define parameters
μ=2.
σ=2.
L=700

log_det_sum = 0.0
for i in range(L):
    m = np.random.uniform(μ - σ * 3**(1/2), μ + σ * 3**(1/2))
    T_i = np.array([[-1, -m/2], [1, 0]])
    # 计算当前矩阵的行列式绝对值，取对数后累加
    det_i = abs(np.linalg.det(T_i))
    log_det_sum += np.log(det_i)

# 指数还原得到最终行列式绝对值
Det = np.exp(log_det_sum)
print(Det)

进一步优化 #

观察单个矩阵$T[i]$的行列式：$\det(T[i]) = \begin{vmatrix}-1 & -m/2 \ 1 & 0\end{vmatrix} = 0 - (-m/2)*1 = m/2$，因此可以直接计算$\log(|m/2|)$，省去调用np.linalg.det的开销：

import numpy as np

# define parameters
μ=2.
σ=2.
L=700

log_det_sum = 0.0
for i in range(L):
    m = np.random.uniform(μ - σ * 3**(1/2), μ + σ * 3**(1/2))
    # 直接利用行列式公式计算对数
    log_det_sum += np.log(abs(m / 2))

Det = np.exp(log_det_sum)
print(Det)

这样计算的结果会和第二个代码的输出一致，量级为$e^{-60}$，完全避免了数值溢出问题。

内容的提问来源于stack exchange，提问作者kBoltzmann