求解100×100概率矩阵的实左特征向量与特征值遇性能瓶颈

首先明确：你的矩阵是行随机矩阵（每行和为1），属于马尔可夫链转移矩阵范畴，这类矩阵有明确的数值求解方法，完全不需要用SymPy做符号计算（100阶符号计算的复杂度是阶乘级，根本无法完成）。针对你遇到的问题，给出以下可行方案：

核心思路 #

1. 行随机矩阵的左特征值λ=1必然存在，对应的左特征向量就是马尔可夫链的稳态分布，这是你最可能关心的结果。
2. NumPy/SciPy返回复数是浮点计算的精度误差，实特征值会带极小虚部，只需筛选处理即可。
3. 大矩阵优先用数值迭代法或稀疏矩阵求解器，避免直接计算高阶行列式或特征多项式。

方案1：幂迭代法求稳态分布（最高效） #

幂迭代专门用于求矩阵的主特征值（模最大的特征值，这里就是λ=1）和对应特征向量，计算量极小，100阶矩阵迭代几百次就能收敛：

import numpy as np

# 加载矩阵
T = np.loadtxt('rep10_T_ij.dat', delimiter=' ')
n = T.shape[0]

# 初始化迭代向量（均匀分布）
pi = np.ones(n) / n

# 迭代求解稳态分布（左特征向量λ=1）
for _ in range(1000):
    pi_new = pi @ T
    # 归一化保证是概率分布
    pi_new /= np.sum(pi_new)
    # 收敛判断（误差小于1e-10时停止）
    if np.linalg.norm(pi_new - pi) < 1e-10:
        pi = pi_new
        break
    pi = pi_new

# 输出结果
print("稳态分布（左特征向量λ=1）：")
print(pi)
print("\n验证πT=π：", np.allclose(pi @ T, pi))

方案2：处理SciPy的复数结果，筛选实特征值 #

如果需要所有实特征值和对应左特征向量，可以对SciPy的计算结果做虚部过滤：

import numpy as np
from scipy.linalg import eig

# 加载矩阵
T = np.loadtxt('rep10_T_ij.dat', delimiter=' ')

# 求解左特征值与左特征向量
eigenvalues, left_eigenvectors = eig(T, left=True, right=False)

# 筛选虚部极小的"近似实特征值"（浮点误差导致的虚部）
real_threshold = 1e-10
real_mask = np.abs(np.imag(eigenvalues)) < real_threshold
real_eigs = np.real(eigenvalues[real_mask])
real_left_vecs = np.real(left_eigenvectors[:, real_mask])

# 可选：筛选马尔可夫链中模为1的实特征值（这类特征值对应长期行为）
unit_threshold = 1e-10
unit_mask = np.abs(np.abs(real_eigs) - 1) < unit_threshold
unit_eigs = real_eigs[unit_mask]
unit_left_vecs = real_left_vecs[:, unit_mask]

# 归一化特征向量（概率分布求和为1）
for i in range(unit_left_vecs.shape[1]):
    vec = unit_left_vecs[:, i]
    vec /= np.sum(vec)
    unit_left_vecs[:, i] = vec

# 输出结果
print("模为1的实特征值：", unit_eigs)
print("对应的归一化左特征向量：")
print(unit_left_vecs)

方案3：稀疏矩阵优化（如果矩阵存在大量0元素） #

如果你的概率矩阵是稀疏的（很多元素接近0），转成稀疏矩阵用scipy.sparse.linalg.eigs求解会更快：

import numpy as np
from scipy.sparse import csr_matrix
from scipy.sparse.linalg import eigs

# 加载矩阵并转为稀疏格式
T = np.loadtxt('rep10_T_ij.dat', delimiter=' ')
T_sparse = csr_matrix(T)

# 求解T^T的最大模特征值（对应T的左特征向量）
eigenvalues, left_eigenvectors = eigs(T_sparse.T, k=1, which='LM')

# 处理结果，取实部并归一化
lambda_1 = np.real(eigenvalues[0])
pi = np.real(left_eigenvectors[:, 0])
pi /= np.sum(pi)

# 输出结果
print("主特征值λ=1：", lambda_1)
print("稳态分布：")
print(pi)
print("\n验证πT=π：", np.allclose(pi @ T, pi))

为什么SymPy行不通？ #

SymPy的符号计算基于精确代数运算，100阶矩阵的行列式展开涉及100!项，这是天文数字级的计算量，无论多少硬件资源都无法在合理时间内完成。大矩阵的特征值求解必须用数值方法，这是行业共识。

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