高斯消元法在复杂动态规划中的应用与实现

核心逻辑与数学原理

高斯消元（Gaussian Elimination）是线性代数中求解线性方程组的基础算法。在算法竞赛中，当概率/期望动态规划（DP）中存在复杂的环路依赖（例如图上的随机游走），导致无法通过拓扑排序或记忆化搜索顺序填表时，高斯消元便成为了解决问题的有效工具。

1. 期望矩阵方程组的构建

设图中有 $n$ 个状态，我们要维护每个状态到终点的期望代价 $f[i]$。根据期望的单步转移规律（参见 21.3 章），每个状态都可以写成关于其他状态的一次代数式：

$$f[i] = \sum_{j=1}^n p_{i \to j} \cdot f[j] + w_i$$

其中 $p_{i \to j}$ 是从状态 $i$ 转移到状态 $j$ 的概率，$w_i$ 是单步付出的期望代价。

我们将未知数 $f[j]$ 统一移项到等号左侧，常数项留在右侧，即可变换为标准线性方程组形式：

$$(1 - p_{i \to i})f[i] - \sum_{j \neq i} p_{i \to j} f[j] = w_i$$

这就构成了包含 $n$ 个未知数、 @@@MATH_BLOCK10@@@ 个方程的增广矩阵（Augmented Matrix）： $$A = \begin{pmatrix} a{1,1} & a{1,2} & \dots & a{1,n} & | & b1 \ a{2,1} & a{2,2} & \dots & a{2,n} & | & b2 \ \vdots & \vdots & \ddots & \vdots & | & \vdots \ a{n,1} & a{n,2} & \dots & a{n,n} & | & b_n \end{pmatrix}$$

2. 消元算法核心步骤

高斯消元的核心是通过矩阵的初等行变换，将增广矩阵转化为上三角矩阵（Upper Triangular Matrix），进而通过回代（Back-substitution）求出所有未知数。

1. 主元选择（Pivoting）：在处理第 $r$ 列时，为了减小浮点数除法带来的数值算术误差，并在中途拦截系数为 $0$ 的情况，需要从第 $r$ 行到第 $n$ 行中，找出第 $r$ 列绝对值最大的那一行，将其交换到第 $r$ 行。
2. 消元发射：利用当前第 $r$ 行的方程，将下方所有行（第 $r+1$ 到 $n$ 行）的第 $r$ 列系数全部强行消减为 $0$。
3. 逆向回代：完成上三角化后，最后一行的方程已经变成了单变量方程 $a_{n,n} \cdot f[n] = b_n$，可直接算出 $f[n]$。然后从下往上层层带入前驱方程，即可在 $O(n^3)$ 时间内解除整个期望网络。

算法推导与状态设计

概率游走模型的矩阵映射

假设在一个包含自环和双向边的图上，我们要求从起点 $1$ 走到终点 $n$ 的期望总步数。

• 状态设计：设 $f[i]$ 表示从节点 $i$ 出发走到终点 $n$ 的期望步数。
• 边界状态方程：终点 $n$ 是明确的，它到自己的期望步数为 $0$。因此第 $n$ 行的方程非常特殊，不需要看任何连边，直接写死为：

$$1 \cdot f[n] = 0$$

• 普通状态方程：对于任意非终点节点 $i$，设其出度为 $deg[i]$。它会等概率地走向任何一个邻接点 $v \in to(i)$。单步代价为 $1$。方程为：

$$f[i] = \sum_{v \in to(i)} \frac{1}{deg[i]} \cdot f[v] + 1 \implies f[i] - \sum_{v \in to(i)} \frac{1}{deg[i]} f[v] = 1$$

将这些系数依次填入矩阵的 A[i][v] 中，A[i][i] 初始化为 $1$，常数项 A[i][n+1] 初始化为 $1$。

C++ 标准源码

#include <iostream>
#include <vector>
#include <cmath>
#include <iomanip>
#include <algorithm>

const int MAXN = 200;
const double EPS = 1e-9;

double A[MAXN + 5][MAXN + 6]; // 增广矩阵，大小为 N * (N + 1)
double ans[MAXN + 5];         // 存储最终每个状态的期望值

// 高斯消元核心算法，返回 false 表示无解或有无数解
bool gaussian_elimination(int n) {
    for (int r = 1; r <= n; ++r) {
        // 1. 寻找主元
        int pivot = r;
        for (int i = r + 1; i <= n; ++i) {
            if (std::abs(A[i][r]) > std::abs(A[pivot][r])) {
                pivot = i;
            }
        }

        // 交换主元行到当前行
        if (pivot != r) {
            for (int j = 1; j <= n + 1; ++j) {
                std::swap(A[r][j], A[pivot][j]);
            }
        }

        // 主元为 0，说明该列无法消元（矩阵奇异）
        if (std::abs(A[r][r]) < EPS) {
            return false;
        }

        // 2. 将当前行的主元系数化为 1
        double divisor = A[r][r];
        for (int j = r; j <= n + 1; ++j) {
            A[r][j] /= divisor;
        }

        // 3. 向下消元
        for (int i = r + 1; i <= n; ++i) {
            if (std::abs(A[i][r]) > EPS) {
                double factor = A[i][r];
                for (int j = r; j <= n + 1; ++j) {
                    A[i][j] -= factor * A[r][j];
                }
            }
        }
    }

    // 4. 逆向回代
    for (int i = n; i >= 1; --i) {
        ans[i] = A[i][n + 1];
        for (int j = i + 1; j <= n; ++j) {
            ans[i] -= A[i][j] * ans[j];
        }
    }
    return true;
}

int main() {
    std::ios_base::sync_with_stdio(false);
    std::cin.tie(NULL);

    int n, m;
    if (std::cin >> n >> m) {
        std::vector<int> deg(n + 1, 0);
        std::vector<std::pair<int, int>> edges;

        for (int i = 0; i < m; ++i) {
            int u, v;
            std::cin >> u >> v;
            edges.push_back({u, v});
            deg[u]++; // 统计出度
        }

        // 初始化增广矩阵
        // 为 1 到 n-1 号节点构建期望转移方程
        for (int i = 1; i < n; ++i) {
            A[i][i] = 1.0;
            A[i][n + 1] = 1.0; // 单步代价 1
        }

        // 填入转移概率系数
        for (const auto &edge : edges) {
            int u = edge.first;
            int v = edge.second;
            if (u < n) { // 终点不需要转移
                A[u][v] -= 1.0 / deg[u];
            }
        }

        // 特殊边界：终点 n 的方程写死为 1 * f[n] = 0
        A[n][n] = 1.0;
        A[n][n + 1] = 0.0;

        if (gaussian_elimination(n)) {
            // 输出从起点 1 出发走到终点 n 的期望步数
            std::cout << std::fixed << std::setprecision(6) << ans[1] << "\n";
        } else {
            std::cout << "INF\n";
        }
    }
    return 0;
}

NOIP 实战避坑指南

1. 错把终点的行方程和其他普通行混合消元 终点 $n$ 到终点的期望代价是明确的 $0$。在构建矩阵时，绝对不能像普通节点那样对终点也执行 A[n][v] -= 1.0/deg[n] 这种转移更新。如果允许终点向后继更新，相当于允许流程在到达终点后“死而复生”继续游走，这会彻底破坏期望的定义，直接导致消元出来的结果演变成毫无意义的垃圾随机数。
2. 没有选择主元引发的浮点数大数吞小数（精度溃败） 部分选手为了省事，在消元时省略了“寻找绝对值最大行并交换”的步骤（即不进行 Pivoting ）。这在基准测试中可能混过去，但在 NOIP 的刁钻数据下，若恰好某行的 A[r][r] 接近 $10^{-14}$，直接拿它作除数去消其他行，会把微小的浮点误差放大到 $10^{14}$ 级别，导致整个矩阵的精度在瞬间全面崩溃。

经典 NOIP/洛谷 真题

1. 洛谷 P3232 [HNOI2013] 游走

• 题意描述：给定一个无向连通图，从 $1$ 号节点出发随机游走，走到 $n$ 号节点时停止。每一步都会等概率选择当前节点的一条边走过去。现在你可以对图中的每条边进行 $1 \sim m$ 的编号（即重新给每条边分配边权），要求使得游走过程中经过的边权总和的数学期望最小。
• 问题本质：节点经过次数期望的矩阵求解与贪心赋权。
• 核心解题思路：根据期望的线性性质，总边权期望等于 $\sum (\text{边的期望经过次数} \times \text{边权})$。因此核心是求出每条边的期望经过次数。而一条边 $(u,v)$ 的期望经过次数，本质上等于从 $u$ 出发的期望次数 $\times \frac{1}{deg[u]}$ + 从 $v$ 出发的期望次数 $\times \frac{1}{deg[v]}$。 于是问题转化为求每个节点的期望经过次数 $f[i]$。通过邻接关系可以建立方程：

$$f[i] = \sum_{(j,i) \in E, j \neq n} \frac{f[j]}{deg[j]} + [i == 1]$$

利用高斯消元求出所有 $f[i]$ 后，反向推导出每条边的期望经过次数。最后根据贪心策略：给期望经过次数越大的边，分配越小的编号，乘积求和即为最小期望。

2. 洛谷 P2973 [USACO10HOL] driving 出行

• 题意描述：一个炸弹在 $1$ 号节点被引爆，它有 $P$ 的概率原地爆炸，有 $1-P$ 的概率随机等概率选择一条边移动到相邻节点。求炸弹最终在每个节点爆炸的绝对概率。
• 问题本质：无终点状态收敛的概率矩阵高斯消元。
• 核心解题思路：本题不属于期望倒推，而是纯粹的概率前向收敛。设 $f[i]$ 表示炸弹到达节点 $i$ 的期望次数。由于每到一个节点，都有 $P$ 的概率在当前位置爆炸，因此第 $i$ 个节点爆炸的真实概率就是 $f[i] \cdot P$。 利用前驱节点的贡献建立方程：

$$f[i] = (1 - P) \sum_{(j,i) \in E} \frac{f[j]}{deg[j]} + [i == 1]$$

构建出 $n$ 个方程的增广矩阵后，直接跑高斯消元，算出的 $f[i] \times P$ 即为每个节点的最终爆炸概率。