NumPy矩阵连乘特征值预期收敛但实际发散问题排查

这个问题是双精度浮点数精度限制+连乘运算误差累积导致的数值不稳定，和数学推导无关：

• 你对T的特征值模长的判断有遗漏：4阶矩阵T共有4个特征值，其中3个模长小于1，剩余1个模长大于1，是连乘过程的主导项。连乘n次后，主导特征值对应的分量会以指数速度增长，当n超过40时，这个分量的量级已经远大于其余3个衰减分量：双精度浮点数仅能保留15~17位有效十进制数字，当主导分量达到1e20量级时，所有理论值小于1e5的分量都会被浮点舍入误差完全覆盖，本来应该衰减到极小值的3个特征值分量，计算结果完全是随机舍入噪声，看起来就像出现了发散。
• 当σ=0时所有T完全相同，理论上$P=T^n$的3个小特征值确实会趋近于0，但朴素矩阵连乘过程中，每一步运算的舍入误差都会被主导特征值的指数增长同步放大，当n足够大时，舍入误差的量级会超过小特征值分量的理论值，最终计算得到的特征值自然不符合预期。
• 额外的代码逻辑bug：你初始化P时已经生成了第一个T矩阵，后续循环又乘了n次T，最终得到的是n+1个矩阵的乘积，和你预期的n次连乘不符，会进一步放大数值误差。

根据使用场景选择对应方法即可，核心思路是避免直接做高次朴素矩阵连乘，从根源上控制数值量级差：

• 固定矩阵场景（σ=0）：不需要做循环连乘，直接对单个T做特征分解，$T^n$的特征值就是T的特征值的n次幂，特征向量和T完全一致，完全规避连乘带来的误差累积。
• 随机矩阵场景（σ>0）：如果需要计算多组随机T的连乘结果，每完成一次矩阵乘法就对P做归一化，将P的整体量级控制在1附近，避免元素量级差过大吞掉小分量的精度。归一化不会改变特征值的相对比例，最后计算特征值时，把每一步累计的归一化系数乘回结果即可。
• 修正初始化逻辑：将P的初始值设为4阶单位矩阵np.eye(4)，再进入循环完成n次矩阵乘法，保证最终结果是n个T的连乘。

修正后参考代码 #

固定矩阵场景（σ=0） #

import numpy as np

μ = 1.5
σ = 0
m = μ
T = np.array([[-1, m, 0.1, 0.1], [1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 1, 0]])
n = 200

# 直接通过特征分解计算T^n的特征值
l, v = np.linalg.eig(T)
lambda_P = l ** n
print("单个矩阵T的特征值：", l)
print(f"{n}个T连乘的特征值：", lambda_P)

随机矩阵场景（σ>0） #

import numpy as np

μ = 1.5
σ = 0.5
n = 200
P = np.eye(4)  # 初始化为单位矩阵
total_norm = 1  # 记录累计归一化系数

for i in range(n):
    m = np.random.normal(μ, σ)
    T = np.array([[-1, m, 0.1, 0.1], [1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 1, 0]])
    P = T @ P
    # 用谱范数做归一化，控制矩阵量级
    step_norm = np.linalg.norm(P, ord=2)
    P = P / step_norm
    total_norm *= step_norm

# 还原归一化对特征值的缩放
λ, w = np.linalg.eig(P)
λ = λ * total_norm
print(f"{n}个随机T连乘的特征值：", λ)

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