线性规划：单纯形算法之处理退化

上一篇文章介绍了单纯形算法（《线性规划：单纯形算法》），但是还有一些遗留问题没有解决，比如退化情形（Degeneracy）。

本文介绍如何处理退化情形。

退化情形

从几何上看上看，造成退化情形的原因是顶点重合。下图是三维空间中非退化情形下的多面体，每个顶点由三个面交汇。

当我们把线段 $ab$缩短至一点，于是得到退化情形。

如上图所示，$a$和$b$重合时四个面交于一点，于是迭代过程可能在 $a$ 和 $b$ 之间转圈（Cycling）。

处理方法

回顾单纯形算法， 根据 Minimum Ratio Test 计算入基变量 $x_j$ 的值。

$$x_j:=\min \left\{\frac{{\tilde{b}}_i}{{\tilde{a}}_{j_i}}\text{ and }{\tilde{a}}_{j_i}>0,i=1,2,\dots ,m\right\}.$$
其中 $\tilde{b}=B^{-1}b$，${\tilde{a}}_j=B^{-1}{a}_j$，$a_j$ 代表矩阵 $A$ 的列向量。

令 $I_0$ 代表上式达到最小值的下标集合。
$$I_0=\left\{r:\frac{{\tilde{b}}_r}{{\tilde{a}}_{j_r}}=\min \left[\frac{{\tilde{b}}_i}{{\tilde{a}}_{j_i}}\text{ and }{\tilde{a}}_{j_i}>0,i=1,2,\dots ,m\right]\right\}.$$

如果 $I_0$ 只包含一个元素 $r$，那么出基变量为 $x_{B_r}$，不会出现退化情形。否则可能出现退化情形，即变量 $x_{B_r}$ 先出基，迭代数次又入基，这样一直循环下去导致算法死循环。

解决思路是防止同一个变量反复出基和入基。当 Minimum Ratio Test 发现 $I_0$ 包含多个出基变量时，用一些规则防止它们反复出基和入基。

下面介绍 Lexicographic Minimum Ratio Test。

在退化的情形下，$I_0$ 包含了多个元素，这此时我们定义新的下标集合 $I_1$：
$$I_1=\left\{r:\frac{{\tilde{a}}_{1_r}}{{\tilde{a}}_{j_r}}=\min \left[\frac{{\tilde{a}}_{1_i}}{{\tilde{a}}_{j_i}}:i\in I_0\right]\right\}.$$
注意：上述分式中的分子 ${\tilde{a}}_1=B^{-1}{a}_1$，而 ${\tilde{a}}_{1_i}$ 代表 ${\tilde{a}}_1$ 中下标为 $i$ 的分量。

如果 $I_1$ 只包含一个元素 $r$，那么交换 $x_j$ 与 $x_{B_r}$ 从而得到新的基。否则，我们定义类似的集合 $I_2,I_3$，直到 $I_k$，使得 $I_k$ 只包含一个元素。
$$I_k=\left\{r:\frac{{\tilde{a}}_{k_r}}{{\tilde{a}}_{j_r}}=\min \left[\frac{{\tilde{a}}_{k_i}}{{\tilde{a}}_{j_i}}:i\in I_{k-1}\right]\right\}.$$
可以证明，这样的 $I_k$ 始终存在。

设 $r$ 是 $I_k$ 中唯一的元素，那么交换 $x_j$ 与 $x_{B_r}$ 从而得到新的基。

通过这种方式，单纯形算法可以在有限步内返回最优解（证明略）。

算法实现

下面用Python实现 Lexicographic Minimum Ratio Test。

单纯形算法的实现这里不重复展示，可以参考上一篇文章。只需要继承之前的类SimplexA，然后重写方法 _minimum_ratio_test(self, j) 即可。

class SimplexAD(SimplexA):
    """
    单纯形算法：处理退化情形。
    """

    def _minimum_ratio_test(self, j):
        """ Lexicographically minimum ratio test.
        :param j: 入基变量 x_j 的下标 j
        :return: 出基变量 x_i 的下标 i
        """
        a_in = np.dot(self._B_inv, self._A[:, j])
        d0 = self._get_d0(a_in)  # 计算 I_0 的数组下标
        if d0 is None:
            return None
        # The index of "basic_vars".
        # The associated basic variable will be moved out.
        i_ind = self._get_out_ind(d0, 0, a_in)
        return self._basic_vars[i_ind]

    def _get_d0(self, a_in):
        """ 根据 Minimum Ratio Test 计算 I_0.
        实际上计算 d0 即可（定义如下）：
        1、计算 ratios = [r_1, r_2, ..., r_m].
        2、计算 d0 = arg min(ratios)，即 ratios 中最小值的下标
        3、I_0 = [self._basic_vars[i] for i in d0]
        """
        b_bar = np.dot(self._B_inv, self._b)
        ratios = list(map(lambda b, a: b / a if a > 1e-6 else np.infty, b_bar, a_in))
        d0 = arg_min(ratios)
        if ratios[d0[0]] == np.infty:
            # a_in 的分量 <= 0
            # 最优目标函数值无界
            return None
        return d0

    def _get_out_ind(self, d0, it, a_in):
        """ 用递归的方法计算出基变量（在self._basic_vars[]中的index）.
        :param d0: I_0 （在self._basic_vars[]的indices）.
        :param it: iteration number, equals 0 at the beginning.
        :param a_in: B_inv * A[:, j]
        :return: the index to "basic_vars" (to be moved out)
        """
        if len(d0) == 1:
            return d0[0]
        a_it = np.dot(self._B_inv, self._A[:, it])
        ratios = [a_it[k] / a_in[k] for k in d0]
        indices = arg_min(ratios)
        d1 = [d0[k] for k in indices]
        return self._get_out_ind(d1, it + 1, a_in)

完整代码

参考文献

[1] Christopher Griffin. Linear Programming: Penn State Math 484 Lecture Notes. Version 1.8.3. Chapter 7.（点此下载）