4.2.2 单纯形

如果 \boldsymbol 没加载出来可以 Ctrl + F5 刷新一下页面……

实现

我们先来讨论，如果所有 \(\boldsymbol{b}\) 中的元素都 \(\ge 0\) 会怎么样。

init

inline void init(int _n, int _m) {
    n = _n, m = _m;
    for (int i = 0; i <= m; i++)
        for (int j = 0; j <= n; j++)
            A[i][j] = 0;
    for (int i = 1; i <= n; i++)
        id[i] = i;
    for (int i = 1; i <= m; i++)
        id[n + i] = 0;
    result = 0;
    for (int i = 1; i <= n; i++)
        ans[i] = 0;
}

首先我们肯定转成了松弛型，所以会新建 \(m\) 个变量 \(x_{n+1},x_{n+2},x_{n+3},\cdots, x_{n+m}\)。然后根据松弛型的矩阵，这当然可以作为一组基。

另外，因为基变量的情况系数是平凡的，我们可以直接不存，只存非基变量的系数。

特别地，我们钦定矩阵的第 \(0\) 列是限制的向量 \(\boldsymbol{b}\)，第 \(0\) 行是贡献系数向量 \(\boldsymbol{c}^T\)，然后：

\[ A_{0,0}=z=-z_0+\sum_{i=1}^{n}c_i\cdot x_i=-z_0 \]

最后一个等号是因为我们会钦定非基变量全为 \(0\)，而这个 \(z_0\) 就是我们在下面的 \(\operatorname{Pivot}\) 过程中给它添加的常数项。

需要注意的是，上面的东西全都是会变的，不会只是初始的那一个。

因为有 \(\boldsymbol{b}\) 中的元素都 \(\ge 0\)，那么初始 \(x\) 全 \(0\) 就是一组可行解。当然，答案 \(\mathrm{res}\) 也为 \(0\)。

\(\operatorname{id}(u)\) 表示最终的矩阵里第 \(u\) 个变量在初始的时候位于的是哪一列（用来构造解）。对于松弛那一组基变量，它们的值不会出现在最终构造的解里，所以把它们的 \(\operatorname{id}\) 都钦定成 \(0\) 是无所谓的。

Pivot

上面的过程基于一个很重要的东西，即 Pivot(x, y)，表示把第 \(x\) 个基变量替换成第 \(y\) 个非基变量。

由于基变量的系数非常平凡，无非就是只有一个 \(1\)，所以我们为了实现的简便，只维护非基变量的系数。

给出一个 Pivot 过程的代码

inline void pivot(int x, int y) {
    std::swap(id[n + x], id[y]);

    f128 t = A[x][y];
    A[x][y] = 1;
    for (int i = 0; i <= n; i++)
        A[x][i] /= t;

    for (int i = 0; i <= m; i++) {
        if (i != x && fabs(A[i][y]) > eps) {
            t = A[i][y], A[i][y] = 0;
            for (int j = 0; j <= n; j++)
                A[i][j] -= t * A[x][j];
        }
    }
}

第一行的内容就是交换第 \(x\) 个基变量和第 \(y\) 个非基变量对应的初始变量。

接下来的这几行无非是高斯消元，但是在被替换的第 \(y\) 列却是特殊处理的。这其实是直接在第 \(y\) 列进行了第 \(n+x\) 列变换。

或者说的更清楚一些，考虑第 \(x+n\) 列，它只有在第 \(x\) 行才为 \(1\)，其他行全部为 \(0\)。也就是只看这两列，矩阵是这样的

\[ \begin{bmatrix} \cdots &a_{0,y} &\cdots &0 &\cdots\\ \cdots &a_{1,y} &\cdots &0 &\cdots\\ \cdots &a_{2,y} &\cdots &0 &\cdots\\ \vdots &\vdots &\vdots &\vdots &\vdots\\ \cdots &a_{x,y} &\cdots &1 & \cdots\\ \vdots &\vdots &\vdots &\vdots &\vdots\\ \cdots &a_{m,y} &\cdots &0 & \cdots\\ \end{bmatrix} \]

对第 \(y\) 列消元，然后考察第 \(x+n\) 列的变化。

我们消元先在第 \(x\) 行全部同除 \(a_{x,y}\)，然后其他的每个行 \(i\) 减去 \(a_{i,y}\) 倍的第 \(x\) 行。

\[ \begin{bmatrix} \cdots &0 &\cdots &-\dfrac{a_{0,y}}{a_{x,y}} &\cdots\\ \cdots &0 &\cdots &-\dfrac{a_{1,y}}{a_{x,y}} &\cdots\\ \cdots &0 &\cdots &-\dfrac{a_{2,y}}{a_{x,y}} &\cdots\\ \vdots &\vdots &\vdots &\vdots &\vdots\\ \cdots &1 &\cdots &\dfrac{1}{a_{x,y}} & \cdots\\ \vdots &\vdots &\vdots &\vdots &\vdots\\ \cdots &0 &\cdots &-\dfrac{a_{m,y}}{a_{x,y}} & \cdots\\ \end{bmatrix} \]

然后，我们交换两列。但是注意到交换后的第 \(n+x\) 列和原来一模一样，只有第 \(y\) 列变了。那我们就在第 \(y\) 列原地修改就好了。从而有了上面的代码。

Simplex

然后，是整个程序的主体过程

inline bool simplex() {
    while (true) {
        int u = 0, v = 0;
        f128 minval = 1e18;
        for (int i = 1; i <= n; i++)
            if (A[0][i] > eps) {
                v = i;
                break;
            }

        if (!v)
            break;

        for (int i = 1; i <= m; i++)
            if (A[i][v] > eps && (!u || A[i][0] / A[i][v] < minval))
                minval = A[i][0] / A[i][v], u = i;

        if (!u) {
            std::cout << "Unbounded" << '\n';
            return false;
        } else
            pivot(u, v);
    }
    return true;
}

前几行就是随便找一个贡献系数 \(>0\) 的非基变量，设其编号为 \(v\)，然后接下来是要找到某一个基变量，设其编号为 \(u\)，从而使得 \(\operatorname{Pivot}(u,v)\) 后 \(v\) 变成的基变量的值最小。

假设约束向量是 \(\boldsymbol{b}\)，回忆上面 \(\operatorname{Pivot}(x,y)\) 的消元后的矩阵，是：

\[ \begin{bmatrix} b_1-\dfrac{b_x\cdot a_{0,y}}{a_{x,y}}&\cdots &-\dfrac{a_{0,y}}{a_{x,y}} &\cdots &0 &\cdots\\ b_2-\dfrac{b_x\cdot a_{1,y}}{a_{x,y}}&\cdots &-\dfrac{a_{1,y}}{a_{x,y}} &\cdots &0 &\cdots\\ b_3-\dfrac{b_x\cdot a_{2,y}}{a_{x,y}}&\cdots &-\dfrac{a_{2,y}}{a_{x,y}} &\cdots &0 &\cdots\\ \vdots &\vdots &\vdots &\vdots &\vdots &\vdots\\ \dfrac{b_{x}}{a_{x,y}}&\cdots &\dfrac{1}{a_{x,y}} &\cdots &1 & \cdots\\ \vdots&\vdots &\vdots &\vdots &\vdots &\vdots\\ b_{m}-\dfrac{b_x\cdot a_{m,y}}{a_{x,y}}&\cdots &-\dfrac{a_{m,y}}{a_{x,y}} &\cdots &0 & \cdots\\ \end{bmatrix} \]

着眼于第 \(x\) 行，因为非基变量全都为 \(0\)，而且所有基变量（除了第 \(x\) 个自身）的第 \(x\) 行系数都为 \(0\)，不难发现新的基变量 \(x\) 的值就是第 \(0\) 列的 \(\dfrac{b_x}{a_{x,y}}\)。这也是符合我们上面的推理的（每个基变量都能被表达成一组非基变量的和，并且能够钦定非基变量都为 \(0\)）

如果找不到这样的 \(u\)，那就说明了我们有一个全 \(0\) 行，它完全不受所有等式的约束，可以取到任意大的值，所以我们的线性规划无界，只要让这个变量取到 \(+\infty\) 就好了。

check

我们上面的讨论都建立在给定一个初始解的情况下，那对于一个任意的线性规划，我们不一定有 \(\boldsymbol{b}\ge 0\)，那么怎么构造出这个初始解？

inline bool check() {
    while (true) {
        int u = 0, v = 0;
        for (int i = 1; i <= m; i++)
            if (A[i][0] < -eps && (!u || rng() & 1))
                u = i;
        if (!u)
            break;
        for (int i = 1; i <= n; i++)
            if (A[u][i] < -eps && (!v || rng() & 1))
                v = i;
        if (!v) {
            std::cout << "Infeasible" << '\n';
            return false;
        }
        pivot(u, v);
    }
    return true;
}

这个事情其实非常简单，我们随便找到一个 \(\boldsymbol{b}_i<0\) 的 \(i\)，如果找不到就说明已经构造出一组初始解了。

然后再找有没有某一个非基变量 \(v\) 在替换成基变量的时候能够让这一项 \(\boldsymbol{b_i}\) 变成非负数的。这里的原理可以参照上一节。

如果找不到这样的 \(v\)，那肯定就无解了。否则 \(\operatorname{Pivot}(u,v)\)。

使用一点点随机化，可以让程序跑得更快（虽然肯定还是过不去 \(\text{UOJ}\) 的 \(\text{hack}\)）。