对偶单纯形算法

单纯形算法从一个基本可行解出发，朝着目标函数值下降的方向迭代，直到最优。从对偶的角度来看，原问题目标函数下降的方向，就是对偶问题的对偶解可行的方向，当对偶解可行时，目标函数达到最优。

本文介绍对偶单纯形法。它的思路是从对偶可行解出发，朝着原问题可行的方向迭代，直到原问题可行，于是得到最优解。

对偶可行解

考虑线性规划问题：
$$\begin{array}{rl}\min & c^{T}x\\ \text{s.t.}& Ax=b\\ & x\ge 0\end{array}$$
和它的对偶问题：
$$\begin{array}{rl}\max & b^{T}y\\ \text{s.t.}& A^{T}y\le c\end{array}$$
其中 $c,x\in {\mathbb{R}}^n$，$y\in {\mathbb{R}}^m$，$A\in {\mathbb{R}}^{m\times n}$，$b\in {\mathbb{R}}^m\ge 0$。

令 $B,N$​​​ 分别代表基矩阵和非基矩阵。对偶问题可以写成如下形式：
$$\begin{array}{rl}\max & b^{T}y\\ \text{s.t.}& B^{T}y\le c_B\\ & N^{T}y\le c_N\end{array}$$
回顾原问题的基本解 $x$​，称为 原始解，是这样的形式：
$$x=\left[\begin{array}{c}{x}_B\\ x_N\end{array}\right]=\left[\begin{array}{c}{B}^{-1}b\\ 0\end{array}\right].$$
当 $B^{-1}b\ge 0$ 时，$x\ge 0$，它就是基本可行解。

再看对偶问题，它的基本解 $y$​​，称为 对偶解，是这样的形式：
$$y=(B^T{)}^{-1}{c}_B.$$
检查对偶问题的两个约束：第一个 $B^T{y}_0\le c_B$​​​ 自然成立；第二个约束 $c_N-N^{T}y\ge 0$​​​，即
$$\begin{array}{rl}{c}_N-N^{T}y& =c_N-N^T(B^T{)}^{-1}{c}_B\\ & =c_N-N^T(B^{-1}{)}^T{c}_B\\ & =c_N-(B^{-1}N{)}^T{c}_B\\ & \ge 0.\end{array}$$
令 $\mu :=c_N-(B^{-1}N{)}^T{c}_B$，如果 $\mu \ge 0$，那么 $y$​ 就是对偶问题的基本可行解。

回顾单纯形算法，$\mu$ 实际上是原问题的 Reduced Cost，即原问题目标函数关于 $x_N$ 的导数。因此，对偶可行意味着原始“最优”，即原问题的目标函数值关于 $x_N$ 无法降低。如果原始解 $x$ 也是可行的，即 $B^{-1}b\ge 0$，那么 $x$ 和 $y$ 分别为原问题和对偶问题的最优解。

如何找到初始的对偶可行解？本文不做介绍，详情见文末的参考文献[1]。

如何迭代

对偶可行保证原始解的最优性，但可能损失可行性。迭代的思路就让原始解满足最优性的同时，也变得可行。

已知基矩阵 $B$​，我们介绍出入基的规则。

已知对偶可行解 $y=(B^T{)}^{-1}{c}_B$​。令 $\tilde{b}:=B^{-1}b$，如果 $\tilde{b}\ge 0$​，那么原始解可行，$x$ 和 $y$ 即为原问题和对偶问题的最优解。

否则存在分量 ${\tilde{b}}_i<0$ 。为了让 $x$ 变得可行，自然的想法是让 $B$ 的第 $i$ 行对应的变量 $x_{B_i}$ 出基，所以 $x_{B_i}$ 是出基变量。

那么入基变量如何计算？

先把对偶问题换一种写法。令 $\tilde{y}=B^{T}y$​​​​，于是 $y=(B^T{)}^{-1}\tilde{y}$​​​​​​。代入对偶问题，得到下面的等价形式：
$$\begin{array}{rl}\max & {\tilde{b}}^T\tilde{y}\\ \text{s.t.}& \tilde{y}\le c_B\\ & (B^{-1}N{)}^T\tilde{y}\le c_N\end{array}$$
前面假设 ${\tilde{b}}_i<0$​​，那么减少 ${\tilde{y}}_i$​​ 可以增加目标函数值，同时需要保证
$$(B^{-1}N{)}^T\tilde{y}\le c_N.$$
令 $J$​ 代表非基变量的下标。$\forall j\in J$​，令 ${\tilde{a}}_j=B^{-1}{a}_j$​，其中 $a_j$​ 代表 $A$​ 的第 $j$​ 列。上面的约束条件可以写成：
$${\tilde{a}}_j^T\tilde{y}\le c_j,\forall j\in J.$$
令 ${\tilde{y}}_i$ 减小 $\delta$，其中 $\delta >0$，我们有
$${\tilde{a}}_j^T\tilde{y}-{\tilde{a}}_{ij}\delta \le c_j,\forall j\in J.$$
如果 ${\tilde{a}}_{ij}>0$​​，$\forall j\in J$，那么 $\delta$ 可以为无穷大，对偶问题的最优目标函数值为$+\infty$，那么原问题无可行解。

假设存在 $j\in J$​​​，使得 ${\tilde{a}}_{ij}<0$​​​​。为了满足约束条件，我们有​​
$$\delta :=\min \left\{\frac{c_j-{\tilde{a}}_j^T\tilde{y}}{-{\tilde{a}}_{ij}}\text{ and }{\tilde{a}}_{ij}<0,\forall j\in J\right\}.$$
注意到
$$\begin{array}{rl}{\tilde{a}}_j^T\tilde{y}& =(B^{-1}{a}_j{)}^T{B}^{T}y\\ & =a_j^T(B^{-1}{)}^T{B}^{T}y\\ & =a_j^{T}y,\end{array}$$
我们有
$$\delta =\min \left\{\frac{c_j-a_j^{T}y}{-{\tilde{a}}_{ij}}\text{ and }{\tilde{a}}_{ij}<0,\forall j\in J\right\}.$$
上式称为 Minimum Ratio Test，取得最小值对应的下标 $j$​ 即为入基变量 $x_j$​ 的下标。

算法描述

第0步：输入对偶可行的基矩阵 $B$。

第1步：判断原始解是否可行。如果 $\tilde{b}\ge 0$​​，当前是最优解，算法停止。

第2步：计算入基变量和出基变量。如果存在 ${\tilde{b}}_i<0$​，那么 $x_i$​ 是出基变量。如果存在 $j\in J$​​​​​ 使得 ${\tilde{a}}_{ij}<0$​​​​​，根据 Minimum Ratio Test， 找到入基变量 $x_j$​​​​​。

第3步：判断问题是否无界。如果 ${\tilde{a}}_{ij}\ge 0$​，$\forall i$​，则对偶问题无界，原问题无可行解，算法停止。

第4步：执行出入基操作，更新基矩阵 $B$，然后执行第1步。

算法实现

下面我们用Python来实现对偶单纯形算法。

先定义算法的输入和输出。

class DualSimplex(object):
    """
        对偶单纯形算法（基本版）。
        Note:
            1、系数矩阵满秩。
            2、未处理退化情形。
            3、输入对偶可行解（对应的列）。
    """

    def __init__(self, c, A, b, v1):
        """
        :param c: n * 1 vector
        :param A: m * n matrix
        :param b: m * 1 vector
        :param v1: dual feasible solution, list of column indices
        """
        # 输入
        self._c = np.array(c)
        self._A = np.array(A)
        self._b = np.array(b)
        self._basic_vars = v0
        self._m = len(A)
        self._n = len(c)
        self._non_basic_vars = self._init_non_basic_vars()
        # 辅助变量
        self._iter_num = 0
        self._B_inv = None
        self._N_bar = None  # N_bar = B^{-1}N
        self._y = None  # dual solution (i.e., shadow price)
        # 输出
        self._sol = None  # primal solution
        self._obj_dual = None  # dual objective
        self._status = None

接下来要实现对偶单纯形算法 DualSimplex.solve()，思路如下。

  class DualSimplex(object):
      
      # ...
      # 其它函数省略……
      
      def solve(self):
          self._iter_num = 0  # 记录迭代次数
          self._check_init_solution()  # 检查初始的对偶解是否可行
          self._update_solutions()
          self._update_obj()
          self._print_info()
          while not self._is_optimal():  # 判断是否最优或者不可行
              if self._status == "INFEASIBLE":
                  break
              self._pivot()  # 迭代（选主元入基，执行Minimum Ratio Test，然后出基）
              self._update_solutions()
              self._update_obj()
              self._iter_num += 1
          if self._status != 'INFEASIBLE':
              self._status = 'OPTIMAL'

上面的关键步骤是实现每次迭代的出入基操作，即 SimplexA._pivot()。

  class DualSimplex(object):
      
      # ...
      # 其它函数省略……
      
      def _pivot(self):
            # 出基变量 x_i
            i = int(np.argmin(self._sol))
            # 出基变量 x_i 在 self._basic_vars 中的index
            i_ind = None
            for k in range(self._m):
                if self._basic_vars[k] == i:
                    i_ind = k
                    break
            # 判断原问题是否可行
            # 对偶问题无界，则原问题不可行
            if not self._check_feasibility(i_ind):
                self._status = 'INFEASIBLE'
                return
            # 入基变量 x_j 在 self._basic_vars 中的index
            j_ind = self._minimum_ratio_test(i_ind)
            # 入基变量 x_j
            j = self._non_basic_vars[j_ind]
            # update basic vars
            for k in range(self._m):
                if self._basic_vars[k] == i:
                    self._basic_vars[k] = j
                    break
            # update non basic vars
            self._non_basic_vars[j_ind] = i

        def _check_feasibility(self, row_ind):
            N = np.array([self._A[:, j] for j in self._non_basic_vars]).transpose()
            self._N_bar = self._B_inv @ N
            for x in self._N_bar[row_ind]:
                if x < 0:
                    return True
            return False

        def _minimum_ratio_test(self, ind_out):
            N = np.array([self._A[:, j] for j in self._non_basic_vars]).transpose()
            c_N = np.array([self._c[j] for j in self._non_basic_vars])
            c_bar = c_N - self._y @ N
            self._N_bar = self._B_inv @ N
            a_bar = self._N_bar[ind_out] * -1
            ratios = list(map(lambda c, a: c/a if a > 0 else np.infty, c_bar, a_bar))
            return int(np.argmin(ratios))

完整代码

结语

本文介绍了对偶单纯形法，但是还有两个问题没有解决：一个是对偶可行解的初始化；另一个处理退化情形。解决思路与单纯形法的处理思路类似，本文不做介绍，详情可以看下面的参考文献[1]。

我们已经知道了单纯形算法，为什么还需要对偶单纯形算法？它有什么样的应用场景？

具体来说，有如下两点：

1、对偶单纯形法在实际中表现不错。不仅如此，原问题如果退化，其对偶问题可能是正常的。

2、对偶单纯形法可应用于整数规划的求解。分支定界是求解整数规划问题的经典方法，每个分支需要求解一个线性规划子问题。子问题相比原来问题，增加了新的约束，一般导致原始解不可行，却不会违背对偶可行。于是可以用之前的对偶解作为子问题的初始可行解。

参考文献

[1] Mihai Banciu. Dual Simplex (lecture notes). Bucknell University, 2010.
[2] David P. Williamson. ORIE 6300 Mathematical Programming I (lecture notes). 2014.