skaiuijing

前言

在编译原理一书中，有一个困扰了笔者很久的算法，就是附录B中的求解线性独立解算法。

笔者第一次翻看时就注意到了这个算法，龙书的附录A是一个编译器前端，附录B则是这个算法。这个算法困扰了笔者很久，然而笔者的数学功底实在是糟糕，初看时未得到答案，甚至到现在再看时，仍然是一知半解。

这本书，在书架中积满灰尘，一个月、两个月、一年、两年，以至于我早已遗忘了它。直到最近，笔者一瞥间，才注意到了它，拂去灰尘，重新细细翻看，在顺着笔记翻阅时，不由得又想起了这个算法。

书中也并未介绍这个算法的原理，只是将其作为工具使用，但是算法的伪代码却激发了笔者的疑惑，思来想去，还是决定写下这么一篇博客用作记录。

问题

1
2
3
4
5
6
7
8
9
附录 B 寻找线性独立解

算法 B.1：找到 Ax >= 0 的线性独立解的最大集合，并将它们表示为矩阵 B 的行。

输入：一个 m x n 矩阵 A。

输出：一个线性独立解的矩阵 B，使得 Ax >= 0。

方法：算法如下伪代码所示。注意，X[y] 表示矩阵 X 的第 y 行，X[y : z] 表示矩阵 X 的第 y 行到第 z 行，X[y : z][u : v] 表示矩阵 X 在行 y 到 z 和列 u 到 v 之间的矩形。

现在我们需要一些工具来解决它。

零空间

笔者初看时，立刻想到了线性代数中常见的Fx=0的问题，常见的解题思路是先对矩阵进行消元，这一步可以采用高斯消元法进行：

从左上角的元素（主元）开始，利用主元消去该列下面的所有元素。

对每一列重复上述过程，直到将矩阵上三角化。

第二步：回代

矩阵中有主列，也有自由列，我们可以给自由列赋值，这样就可以得到主列变量的值，也可以求得基础解。

那么方程的解空间就是两个基础解的线性组合。

但是，我们要求解的是Ax>=0，这明显是一个线性规划问题。

在运筹学中，线性规划问题可以使用单纯形法解决。

单纯形法

单纯形法就是通过设置不同的基向量，经过矩阵的线性变换，求得基可行解，并判断该解是否最优，否则继续设置另一组基向量，重复执行以上步骤，直到找到最优解。

单纯形法的一般解题步骤可归纳如下：

  1. 把线性规划问题的约束方程组表达成典范型方程组，找出基本可行解作为初始基本可行解。
  2. 若基本可行解不存在，即约束条件有矛盾，则问题无解。
  3. 若基本可行解存在，从初始基本可行解作为起点，根据最优性条件和可行性条件，引入非基变量取代某一基变量，找出目标函数值更优的另一基本可行解。
  4. 按步骤3进行迭代，直到对应检验数满足最优性条件（这时目标函数值不能再改善），即得到问题的最优解。
  5. 若迭代过程中发现问题的目标函数值无界，则终止迭代。

简单举例：

求最大利润P=3x + 5y

约束如下：

1
2
3
4
x <= 4
y <= 6
3x + 2y <= 18
x,y > 0

现在应用单纯形法：

1.从原点出发(0,0)出发，我们会发现y方向的利润增长较x方向快，沿y方向移动到(0,6)，此时利润30

2.由于约束3x + 2y成立，再沿x方向移动到(2,6)，此时利润36

3.观察其他顶点，显然，(2,6)是最优解成立。

此外，还有多种方法可以进行求解：

线性约束表示多面体

多面体可以通过一组线性不等式来定义，这些不等式可以表示为矩阵和向量的形式。

从空间角度来看，Ax=0的解是一个线性子空间，而Ax >= 0的解是一个凸锥，我们要求线性独立解的最大集合，其实就是求极射线。

对于Ax >= 0，可以表示成一个凸多面体，同理，顶点就是满足约束条件的解，使用单纯形法在多顶点之间移动，我们可以获得最优解。

但是现在还有一个问题，如何求得移动的范围？难道依靠穷举吗？由于我们往往求解的是丢番图方程，因此可以先使用Fourier-Motzkin算法求得约束范围或者简化不等式，再使用单纯性法进行求解。

Farkas引理

Farkas 引理有几种不同的形式，但其基本思想是关于线性不等式系统的可行性。最常见的一种形式如下：

考虑以下线性不等式系统： A x ≤ b

其中，A 是一个 m x n 的矩阵，x 是一个 n 维向量，b 是一个 m 维向量。Farkas 引理断言，对于这个系统：

  1. 存在一个向量 x 使得 A x ≤ b
  2. 存在一个向量 y 使得 y ≥ 0 并且 y^T A = 0 以及 y^T b < 0。

这两种情况必有一种成立，但不会同时成立。

单纯形法的最优化过程

如果 b 向量所有元素非负，则显然我们只需要令所有的变量等于 0，就可以得到一个可行解。在这种情况下，通过下述最优化过程，我们可以得到该线性规划的最优解，或者指出该线性规划的最优解为无穷大（不存在）。

  1. 任取一个非基变量 x_e，使得 (c_e > 0)。
  2. 选取一个基变量 (x_d)，使得 (A{d,e} > 0),且最小化b_d/{A_{d,e}})。
  3. 执行转轴操作 ，并转到第一步继续算法。

算法终止有两种情况：

  1. 对于所有的非基变量，c 均非正。
  2. 对于某一个 e，所有的 (A_{d,e}) 均非正。

可以证明，对于第一种情况，我们已经得到了该线性规划的最优解。当前的 v 即为答案。严格证明比较复杂，但是直观上是很容易理解的。因为所有的非基变量都是非负的，而所有的 c 都是非正的，因此只要某个非基变量不为 0，就会使得目标函数更小。

对于第二种情况来说，很容易证明此时线性规划的最优解是无穷大。只要让其他所有变量均为 0，变量 (x_e) 为正无穷。由于所有的 (A_{d,e}) 都非正，因此非基变量的非负性得到保证。同时由于 (c_e > 0)，目标函数值为正无穷。

实现

由于我们要将极射线表示为矩阵B的行，因此先将A转化为转置矩阵AT,

实现思路如下：

消元

使用高斯消元法，转化为上三角对角矩阵。

寻找解

我们可以使用上面提到过的单纯形法和Fourier-Motzkin算法寻找满足kro, …, kr’ - 1 ≥ 0的组合。

理论上算法到这一步就算结束了。但是，需要考虑解是否合理。

考虑解的合理性

考虑到我们的算法将会用来解决实际优化问题，因此，负数是不合理的，我们需要进行非负性强制调整。可以通过行操作来消除负元素，确保子矩阵的非负性。

这一步类似于单纯形法的转轴操作，调整基变量以保证解的可行性。

终止是否成立

当无法生成新解或处理完毕时，返回矩阵，否则继续迭代。

伪代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
M = AT
ro = 1
co = 1
B = Inxn  // n x n 单位矩阵

while (true) {
    // 1. 将 M[ro : r' - 1][co : c' - 1] 转换为对角矩阵
    // 并且 M[r' : n][co : m] = 0，M[r' : n] 是解
    r' = ro
    c' = co
    while (存在 M[r][c] ≠ 0 且 r - r' 和 c - c' 均 ≥ 0) {
        // 将主元 M[r][c] 移动到 M[r'][c']，并将行和列交换
        交换行 r 和 r' 在 B 中
        if (M[r'][c'] < 0) {
            M[r'] = -1 * M[r']
            B[r'] = -1 * B[r']
        }
        for (row = ro to n) {
            if (row ≠ r' and M[row][c'] ≠ 0) {
                u = -(M[row][c'] / M[r'][c'])
                M[row] = M[row] + u * M[r']
                B[row] = B[row] + u * B[r']
            }
        }
        r' = r' + 1
        c' = c' + 1
    }

    // 2. 找到 M[r' : n] 之外的解。它必须是 M[ro : r' - 1][co : m] 的非负组合
    
    寻找 kro, ..., kr' - 1 ≥ 0 使得
    kroM[ro][c' : m] + ... + kr' - 1M[r' - 1][c' : m] ≥ 0
    
    if (存在非平凡解, 假设 kr > 0) {
        M[r] = kroM[ro] + ... + kr' - 1M[r' - 1]
        NoMoreSoln = false
    } else { // M[r' : n] 为唯一解
        NoMoreSoln = true
    }

    // 3. 令M[ro : rn - 1][co : m] ≥ 0
    if (NoMoreSoln) {
        // 将解 M[r' : n] 移动到 M[ro : rn - 1]
        for (r = r' to n) {
            交换行 r 和 ro + r - r' 在 M 和 B 中
        }
        rn = ro + n - r' + 1
    } else {
        // 使用行相加找到更多解
        rn = n + 1
        for (col = c' to m) {
            if (存在 M[row][col] < 0 且 row ≥ ro) {
                if (存在 M[r][col] > 0 且 r ≥ ro) {
                    for (row = ro to rn - 1) {
                        if (M[row][col] < 0) {
                            u = ceil(-M[row][col] / M[r][col])
                            M[row] = M[row] + u * M[r]
                            B[row] = B[row] + u * B[r]
                        }
                    }
                }
                for (row = rn - 1 to ro step -1) {
                    if (M[row][col] < 0) {
                        rn = rn - 1
                        交换 M[row] 和 M[rn]
                        交换 B[row] 和 B[rn]
                    }
                }
            }
        }
    }

    // 4. 令M[ro : rn - 1][1 : co - 1] ≥ 0
    for (row = ro to rn - 1) {
        for (col = 1 to co - 1) {
            if (M[row][col] < 0) {
                选择 r 使得 M[r][col] > 0 且 r < ro
                u = ceil(-M[row][col] / M[r][col])
                M[row] = M[row] + u * M[r]
                B[row] = B[row] + u * B[r]
            }
        }
    }

    // 5. 如有必要，重复处理 M[rn : n] 的行
    if (NoMoreSoln or rn > n or rn == ro) {
        删除 B 中 rn 到 n 的行
        返回 B
    } else {
        Cn = m + 1
        for (col = m to 1 step -1) {
            if (不存在 M[r][col] > 0 且 r < rn) {
                Cn = Cn - 1
                交换列 col 和 Cn 在 M 中
            }
        }
        ro = rn
        co = Cn
    }
}