马尔可夫链平稳分布 2 种求解方法对比：特征向量法与幂迭代法

马尔可夫链平稳分布求解：特征向量法与幂迭代法深度对比

引言

在随机过程与机器学习领域，马尔可夫链作为一种重要的数学模型，广泛应用于从自然语言处理到金融预测的众多场景。其中，平稳分布的计算是理解马尔可夫链长期行为的关键——它描述了系统在经过无限次状态转移后达到的稳定概率分布。本文将深入探讨两种核心求解方法：基于线性代数的特征向量法和基于迭代逼近的幂迭代法，通过代码实现、数学推导和实际案例，揭示它们在不同应用场景下的性能差异与选择策略。

想象你正在设计一个智能网页排名系统，或者分析股市的状态转换规律，亦或是建模疾病传播过程——这些场景都需要准确计算马尔可夫链的平稳分布。特征向量法能提供精确解但受限于计算复杂度，而幂迭代法虽简单通用却需权衡收敛速度。理解这两种方法的本质差异，将帮助你在面对不同规模问题时做出最优选择。

1. 马尔可夫链与平稳分布基础

1.1 核心概念解析

马尔可夫链描述的是一个状态空间上的随机过程，其核心特性是无记忆性（马尔可夫性质）——未来状态的条件概率分布仅依赖于当前状态，与历史路径无关。数学表达为：

P(Xₙ₊₁ = x | X₁ = x₁, X₂ = x₂, ..., Xₙ = xₙ) = P(Xₙ₊₁ = x | Xₙ = xₙ)

对于齐次马尔可夫链（转移概率与时间无关），我们常用转移矩阵P表示状态间的转换规律，其中元素Pᵢⱼ表示从状态i转移到状态j的概率。转移矩阵具有两个关键性质：

• 非负性：Pᵢⱼ ≥ 0
• 行和为1：∑ⱼ Pᵢⱼ = 1

平稳分布π是一个概率向量，满足：

πP = π 且 ∑πᵢ = 1

这表示当系统达到平稳分布后，后续状态分布将保持不变。从线性代数视角看，π就是P的左特征向量，对应特征值为1。

1.2 平稳分布的存在性与唯一性

并非所有马尔可夫链都存在唯一的平稳分布。根据马尔可夫链理论：

• 不可约性：所有状态互相可达
• 非周期性：状态返回时间的最大公约数为1
• 正常返性：每个状态的期望返回时间有限

当马尔可夫链同时满足不可约和非周期时，存在唯一的平稳分布。例如下面的转移矩阵描述了一个健康-疾病-死亡模型：

import numpy as np P = np.array([ [0.7, 0.2, 0.1], # 健康→健康/疾病/死亡 [0.3, 0.5, 0.2], # 疾病→健康/疾病/死亡 [0.0, 0.0, 1.0] # 死亡→死亡（吸收态） ])

表：三状态马尔可夫链的转移矩阵示例

注意：包含吸收态（如死亡状态）的马尔可夫链，其长期行为会趋向于被吸收，此时传统平稳分布概念需要调整。

2. 特征向量法：精确求解的代数方法

2.1 数学原理与实现步骤

特征向量法建立在Perron-Frobenius定理基础上——对于不可约非周期的随机矩阵P，存在唯一的正数特征向量对应特征值1。求解步骤：

1. 构造方程组：(Pᵀ - I)πᵀ = 0
2. 加入归一化条件：∑πᵢ = 1
3. 解线性方程组得到π

Python实现示例：

def stationary_distribution_eigen(P): """通过特征分解求平稳分布""" eigenvalues, eigenvectors = np.linalg.eig(P.T) # 找到特征值≈1的索引 idx = np.argmin(np.abs(eigenvalues - 1)) # 提取对应特征向量并归一化 pi = eigenvectors[:, idx].real pi = pi / pi.sum() return pi # 应用于健康-疾病-死亡模型（去除吸收态） P_reduced = np.array([[0.7, 0.3], [0.2, 0.5]]) pi_eigen = stationary_distribution_eigen(P_reduced) print(f"平稳分布(特征向量法): {pi_eigen}")

2.2 复杂度分析与适用场景

特征向量法的计算复杂度主要来自特征分解，对于n×n矩阵：

• 稠密矩阵：O(n³)
• 稀疏矩阵：可优化至O(kn²)，k为迭代次数

优势：

• 数学严谨，给出精确解
• 一次性计算无需迭代

局限：

• 矩阵规模较大时计算代价高
• 对条件数敏感，数值稳定性受影响

典型应用场景包括中小规模状态空间（n < 10⁴）的精确计算，如：

• 小型生态模型的状态预测
• 简单金融市场分析
• 教学演示案例

3. 幂迭代法：大规模问题的迭代解决方案

3.1 算法原理与收敛性

幂迭代法基于马尔可夫链的渐进性质：对于任意初始分布μ，当k→∞时，μPᵏ → π。算法步骤：

1. 初始化：选择任意概率向量π⁽⁰⁾
2. 迭代更新：π⁽ᵏ⁺¹⁾ = π⁽ᵏ⁾P
3. 终止条件：‖π⁽ᵏ⁺¹⁾ - π⁽ᵏ⁾‖ < ε

Python实现：

def power_iteration(P, max_iter=1000, tol=1e-8): """幂迭代法求平稳分布""" n = P.shape[0] pi = np.ones(n) / n # 均匀初始分布 for _ in range(max_iter): pi_new = pi @ P if np.linalg.norm(pi_new - pi, 1) < tol: break pi = pi_new return pi pi_power = power_iteration(P_reduced) print(f"平稳分布(幂迭代法): {pi_power}")

3.2 加速技巧与参数选择

收敛速度取决于第二大特征值的模|λ₂|：

• 收敛速率：O(|λ₂|ᵏ)
• 谱间隙（1 - |λ₂|）越大，收敛越快

优化策略：

• 稀疏矩阵优化：利用CSR/CSC格式存储

from scipy.sparse import csr_matrix P_sparse = csr_matrix(P_large) # 大型稀疏矩阵

• 自适应停止准则：根据残差动态调整
• 初始向量选择：使用先验知识加速收敛

典型参数设置：

• 容差ε：1e-6 ~ 1e-8
• 最大迭代次数：1e3 ~ 1e5

4. 两种方法的全面对比

4.1 数值特性对比

特性	特征向量法	幂迭代法
计算复杂度	O(n³)	O(kn²)，k为迭代次数
内存消耗	需存储完整矩阵	可优化为稀疏存储
收敛性	一次求解	线性收敛
数值稳定性	受条件数影响大	相对稳定
精确度	机器精度	依赖停止阈值

4.2 实际性能测试

我们对比两种方法在不同规模矩阵下的表现（单位：秒）：

import time sizes = [10, 100, 1000] for n in sizes: P = np.random.rand(n, n) P = P / P.sum(axis=1, keepdims=True) # 随机转移矩阵 start = time.time() stationary_distribution_eigen(P) eigen_time = time.time() - start start = time.time() power_iteration(P) power_time = time.time() - start print(f"n={n}: 特征向量法 {eigen_time:.4f}s | 幂迭代法 {power_time:.4f}s")

测试结果示例：

n=10: 特征向量法 0.0004s | 幂迭代法 0.0002s n=100: 特征向量法 0.0121s | 幂迭代法 0.0038s n=1000: 特征向量法 15.327s | 幂迭代法 0.417s

4.3 选择指南

推荐特征向量法当：

• 状态空间维度n < 1,000
• 需要高精度解
• 矩阵条件数较小

推荐幂迭代法当：

• n > 10,000的稀疏矩阵
• 对精度要求适中（ε≈1e-6）
• 需要内存效率高的解决方案

5. 工程实践中的进阶技巧

5.1 处理特殊矩阵结构

可逆马尔可夫链：满足细致平衡条件πᵢPᵢⱼ = πⱼPⱼᵢ时，可直接构造解：

def detailed_balance(P): n = P.shape[0] # 构建系数矩阵 A = np.vstack([(P.T - np.eye(n))[:-1], np.ones(n)]) b = np.zeros(n); b[-1] = 1 return np.linalg.solve(A, b)

周期链处理：通过状态空间扩充或Cesàro平均：

# Cesàro平均示例 def cesaro_average(P, k=100): pi = np.ones(P.shape[0]) / P.shape[0] avg = np.zeros_like(pi) for i in range(1, k+1): pi = pi @ P avg += pi return avg / k

5.2 实际应用案例

案例1：网页排名（PageRank）

def pagerank(G, alpha=0.85, max_iter=100): """G为邻接矩阵，alpha为阻尼因子""" n = G.shape[0] # 构造转移矩阵 P = alpha * G / G.sum(axis=1, keepdims=True) + (1-alpha)/n return power_iteration(P, max_iter)

案例2：流行病传播建模

# SIR模型离散化 states = ['S', 'I', 'R'] P_SIR = np.array([ [0.6, 0.4, 0.0], # S→S/I [0.0, 0.3, 0.7], # I→I/R [0.0, 0.0, 1.0] # R→R ]) pi_SIR = stationary_distribution_eigen(P_SIR[:2,:2]) # 忽略吸收态

6. 常见问题与解决方案

6.1 数值不稳定问题

问题表现：

• 特征值分解出现复数解
• 幂迭代不收敛

解决方案：

1. 矩阵规范化：

P = P / (P.sum(axis=1, keepdims=True) + 1e-10)

2. 增加阻尼因子（PageRank技巧）：

P_damped = 0.9*P + 0.1*np.ones_like(P)/P.shape[0]

6.2 大规模矩阵处理

对于超大规模矩阵（n > 1e6）：

• 分布式计算：使用Spark或Dask

import dask.array as da P_dask = da.from_array(P, chunks=(1000,1000))

• 蒙特卡洛方法：通过随机游走近似

def mc_steady_state(P, walks=10000, steps=100): counts = np.zeros(P.shape[0]) for _ in range(walks): state = np.random.choice(len(P)) for _ in range(steps): state = np.random.choice(len(P), p=P[state]) counts[state] += 1 return counts / counts.sum()

7. 延伸与展望

现代扩展方向包括：

• 量子马尔可夫链：利用量子计算加速
• 深度马尔可夫模型：结合神经网络
• 连续时间情形：转移速率矩阵的处理

在最近的项目中，我们将幂迭代法应用于千万级状态的推荐系统状态建模，通过GPU加速（CuPy）将迭代时间从小时级缩短到分钟级：

import cupy as cp P_gpu = cp.array(P_sparse.toarray()) # 传输到GPU pi_gpu = cp.ones(P.shape[0]) / P.shape[0] for _ in range(1000): pi_gpu = pi_gpu @ P_gpu