处理numpy中矩阵乘法的舍入误差

Question

TL;DR: 状态转移矩阵的矩阵乘法应该保范，但是np.matmul不保范。我怎样才能解决这个问题？有更好的 python 模块吗？

我有一个正确的状态转移矩阵 A，即 s(t)A(tau)=s(t+tau)

其中 s(t) 是一个总和为 1 的列矩阵。另外，我们知道 A 的每一行加起来也是 1。

我们知道A^n也是自然数中任意n的右状态转移矩阵

找到稳态分布的一种方法是在 n 趋于无穷大时计算 A^n。以下片段计算 A^(2^n):

    def p_multiplier(A,n):
        B=copy.deepcopy(A)
        for i in range(n):
            B=numpy.matmul(B,B)
        return B[0]

这工作正常，但由于某种原因，行的总和开始不加到 1，即 numpy 开始泄漏 A[i] 的范数。

原因可能是舍入误差。一种快速解决方法可能是在每次迭代时强制保留每一行的范数：

def p_multiplier(A,n):
    B=copy.deepcopy(A)
    for i in range(n):
        B=np.matmul(B,B)
        for j in range(len(B)):
            B[j]=B[j]/sum(B[j])
    return B[0]

此修复的科学性如何？ scipy 还是 mpmath 处理得更好？

编辑：为了完全满足 MWE，您可以对 A 使用以下状态转换矩阵：

A=np.asarray([[0.76554539,0.13929202,0,0,0.09516258],[0.04877058,0.76996018,0.18126925,0,0],[0,0.09516258,0.76554539,0.13929202,0],[0,0,0.13929202,0.67943873,0.18126925],[0.09516258,0,0,0.04877058,0.85606684]])

Answer 1

问题是因为操作数值不稳定。结果，它 迅速发散 （指数）到 0 或无穷大，即使 relatively-small 的 n 值如 70。您可以使用 对角化基于特征值的方法（有关更多信息，请参见here），它在数值上更加稳定。

def stable_multiplier(A, n):
    d, P = numpy.linalg.eig(A)
    return (P @ numpy.diag(d**n) @ numpy.linalg.inv(P))[0].real

丢失的（小数）位数在 O(log10(n)) 中，因此误差与 n 呈线性关系。请注意，如果 n 很大，如 > 1000，那么您需要非常注意整个算法的数值稳定性，特别是如果它是一个迭代过程（具有 at of 迭代）。想了解更多请跟帖this and this.

处理numpy中矩阵乘法的舍入误差

Handling rounding errors in matrix multiplication in numpy

python

numpy

matrix-multiplication