3. 转移概率矩阵的运算与性质

在马尔可夫链的建模中，转移概率矩阵 P 不只是把一步转移概率“排成表格”，它本身就是一个可参与运算的数学对象。理解它的乘法规则、随机性约束、正则性条件以及稀疏存储方式，是从“会写矩阵”走向“能用矩阵算”的关键一步。本节围绕四个相互关联的问题展开：状态分布如何随步数演化、C-K 方程如何把多步转移拆成一步转移的乘积、什么样的矩阵才合法且具有良好极限行为、以及在状态数很大时如何避免稠密矩阵带来的计算浪费。

一、状态分布的演化：行向量右乘转移矩阵

约定状态分布用行向量表示，即 π_t = (π_t(1), π_t(2), …, π_t(n))，其中 π_t(i) 表示第 t 步处于状态 i 的概率，各分量非负且和为 1。转移矩阵 P 的第 i 行第 j 列元素 p_ij 表示“从 i 一步到 j”的概率，因此每一行之和为 1，这类矩阵称为行随机矩阵。在这一约定下，一步演化就是右乘：

π_(t+1) = π_t · P
π_(t+k) = π_t · P^k

展开分量形式可看得更清楚：π_t+1(j) = Σ_i π_t(i)·p_ij，即“到达 j 的总概率 = 从各出发状态 i 来的概率加权和”。这正是矩阵乘法的语义：左矩阵的行与右矩阵的列做内积。若改用列向量约定，则演化为 π_t+1 = P^Tπ_t，此时 P 的列和为 1。两种约定都常见，但混用会导致转置错误，建模时必须先固定一种并在全文保持一致。本节统一采用行向量约定。

举个可手算的示例（示例参数，非实测）：设三状态链

P = [0.7  0.2  0.1
     0.3  0.5  0.2
     0.2  0.3  0.5]

初始分布 π₀ = (0.5, 0.3, 0.2)。计算 π₁：第一分量 = 0.5×0.7 + 0.3×0.3 + 0.2×0.2 = 0.35 + 0.09 + 0.04 = 0.48；第二分量 = 0.5×0.2 + 0.3×0.5 + 0.2×0.3 = 0.10 + 0.15 + 0.06 = 0.31；第三分量 = 0.5×0.1 + 0.3×0.2 + 0.2×0.5 = 0.05 + 0.06 + 0.10 = 0.21。三者之和为 1，验证了行随机矩阵右乘后仍保持概率分布。继续迭代会看到分布逐步向平稳分布收敛，这一过程正是后续稳态分析的计算基础。

二、C-K 方程：多步转移的分解

Chapman-Kolmogorov 方程（C-K 方程）指出，从 i 经 m+n 步到 j 的概率，等于“先经 m 步到某个中间状态 k，再经 n 步从 k 到 j”对所有 k 求和：

p_ij^(m+n) = Σ_k p_ik^(m) · p_kj^(n)
即  P^(m+n) = P^(m) · P^(n)

它的价值在于把“一次跨很多步”的复杂事件拆成两段可计算的转移，从而允许用矩阵乘法递推得到任意步转移矩阵：P² = P·P，P³ = P²·P，依此类推。C-K 方程成立的前提是马尔可夫性——给定当前状态后，未来与过去无关，中间状态 k 的“记忆”被完全压缩在一步转移里。若过程不满足马尔可夫性，该分解就不成立。工程上，C-K 方程也是验证实现正确性的一个手段：随机生成 P 后检查 P² 是否仍为行随机矩阵、P²·P 是否等于 P·P²（矩阵乘法满足结合律，两者应一致）。

三、随机性与正则性：矩阵的“合法性”与“收敛性”

随机性是合法性约束：P 的每个元素非负（p_ij ≥ 0），且每行之和为 1。满足这两条即为行随机矩阵。它保证任意概率分布右乘后仍是概率分布，也保证 P^k 对任意 k 仍是随机矩阵。实现中常见的错误是浮点累加导致行和偏离 1（如 0.9999999），长期迭代会放大误差，因此工程上常在每步或每隔若干步做一次行归一化，或在构造 P 时用“最后一个元素 = 1 − 其余元素之和”的方式强制闭合。

正则性是收敛性条件：若存在某个正整数 k，使 P^k 的所有元素严格大于 0，则称该链为正则链。正则链具有唯一平稳分布，且无论初始分布如何，π_t 都收敛到该平稳分布。直观理解是：经过有限步后，从任何状态出发都可能到达任何状态，链“充分混合”。注意正则性比“不可约且非周期”更容易直接检验——只需计算 P 的幂并观察是否出现全正矩阵。例如上面的三状态示例，P² 已无零元素，故为正则链。反之，若 P 是单位矩阵（每步原地不动），则 P^k 永远只有对角元为 1，不是正则链，也不存在唯一的混合极限。

四、稀疏矩阵存储与计算优化

实际建模中状态数常达上万甚至百万，而每个状态往往只与少数几个状态有直接转移，P 极度稀疏。若用 n×n 稠密数组存储，内存为 O(n²)，一次矩阵乘法为 O(n³)，完全不可行。此时应改用稀疏表示，只存非零元素及其行列索引。以坐标格式（COO）为例，三个数组分别记录行号、列号和值：

行索引 row = [0, 0, 0, 1, 1, 1, 2, 2, 2]
列索引 col = [0, 1, 2, 0, 1, 2, 0, 1, 2]
值     val = [0.7,0.2,0.1,0.3,0.5,0.2,0.2,0.3,0.5]

更常用的是按行压缩的 CSR 格式：用 indptr 标记每行在列索引与值数组中的起止位置，适合“行向量右乘矩阵”的运算模式。稀疏矩阵乘法的复杂度约为 O(nnz·平均每行非零数)，远低于稠密乘法。工程优化要点包括：

• 利用稀疏结构：若每行非零元很少，优先用 CSR 并按行遍历，避免随机访存。
• 避免显式求 P^k：求平稳分布时不必反复乘出完整矩阵，而用幂迭代 π ← π·P 直到 ‖π_t+1 − π_t‖₁ < ε（如 1e−10），每步只做一次稀疏向量乘矩阵。
• 控制数值误差：迭代中定期归一化，防止分布和漂移；用双精度浮点并设置最大迭代次数与收敛判据双重保护。
• 并行与分块：按行分块后各块独立计算再合并，适合多核或分布式场景。

以一个含 10 万个状态、平均每状态 4 个非零转移的示例（合理假设）估算：稠密存储约需 10¹⁰ 个双精度数，约 80 GB，不可接受；CSR 仅存约 4×10⁵ 个非零元，连同索引约数 MB 量级，单次幂迭代可在毫秒到秒级完成。这一对比说明，稀疏存储不是“锦上添花”，而是大规模马尔可夫链能否落地的前提。

把以上四点串起来，就形成了一条完整的计算路径：用行随机矩阵保证合法性，用 C-K 方程与矩阵乘法完成多步演化，用正则性判断收敛前景，用稀疏存储与幂迭代把计算落到可执行的规模。后续的稳态求解与参数估计，都将建立在这套矩阵运算与性质之上。