Large-Scale Linear Optimization

本文介绍大规模线性规划问题的求解思想。内容主要来源自《Introduction to Linear Optimization》Chapter 6.

Delayed column generation

考虑一个非退化的线性规划的标准问题：

$\begin{matrix} (1) & \begin{array}{cl} min & c^{T} x \\ s.t. & A x = b, x \geq 0 \end{array} \end{matrix}$

假如说矩阵 $A \in R^{m \times n}$ 的列数非常大，即 $n ≫ m$ ，这时候我们无法把所有的列都放入内存执行计算。但是，注意到问题只有 $m$ 个基变量，也就是说我们只需要找到 $n$ 列中特定的 $m$ 列，就可以完成单纯型法的迭代，找到最优解了！

我们可以先随便选取一组基变量进行计算，接着我们就要去找 entering variable，找进基变量的准则是 reduced cost ${\bar{c}}_{j} = c_{j} - c_{B}^{T} B^{- 1} A_{j} < 0$ 。当然，一般会选择 ${\bar{c}}_{j}$ 最小的指标 $j$ 进基，这就归结于问题：

$min {\bar{c}}_{j} = c_{j} - c_{B}^{T} B^{- 1} A_{j}$

对于一些具有特殊结构的线性规划，只要如上的问题可以轻松解决，即找到进基指标 $j$ 和相应的列 $A_{j}$ ，那么原问题就可以得到高效地求解！

注：矩阵 $B$ 可以借助 revised simplex 方法进行迭代，因此计算量也是小的。

下料问题是一个经典的、可以使用列生成算法来求解的问题。

Cutting-stock problem

下料问题由Kantorovich在1939年提出，是一个 NP-Hard 问题。

假设一个造纸厂生产长度为 $W$ 的纸，然而， $m$ 个顾客需要的是 $b_{i}$ 卷长度为 $w_{i} < W$ 的纸（ $i = 1, \dots, m$ ）。那么，造纸厂最少需要多少卷纸来满足顾客的需求呢？

为了解决这个问题，用一个列向量 $A_{j} = [a_{i j}]_{m \times 1}$ 表示一卷长度为 $W$ 的纸是如何切割成若干个长度为 $w_{i}$ 的纸的，其中 $a_{i j}$ 表示第 $j$ 种切法中切出长度为 $w_{i}$ 的纸的数量。这样的话，对于一个向量 $[a_{1 j}, a_{2 j}, \dots, a_{m j}]^{T}$ ，它成为一种可行的切法的充要条件是满足：

$A_{j}^{T} w = \sum_{i = 1}^{m} a_{i j} w_{i} \leq W, a_{i j} \in N$

令列分块矩阵 $A = [A_{j}]_{1 \times n} \in R^{m \times n}$ 表示所有可行的切割方法，注意到切法 $n$ 可能非常大，这会给问题带来困难。

下料问题归结为一下优化问题：

$\begin{aligned} min & \sum_{j = 1}^{n} x_{j} \\ s.t. & \sum_{j = 1}^{n} a_{i j} x_{j} \geq b_{i}, i = 1, 2, \dots, m \\ x_{j} \geq 0, x_{j} \in N, j = 1, 2, \dots, n \end{aligned}$

决策变量 $x_{j}$ 表示第 $j$ 种切割方法执行的数量。这是一个整数规划问题，把整数条件放松掉，我们求解一个线性规划问题能轻松得到原整数规划问题的一个上界。（向上取整即可），于是其线性松弛可以表示为：

$\begin{aligned} min & \sum_{j = 1}^{n} x_{j} \\ s.t. & \sum_{j = 1}^{n} a_{i j} x_{j} = b_{i}, i = 1, 2, \dots, m \\ x_{j} \geq 0, j = 1, 2, \dots, n \end{aligned}$

注意到 $n ≫ m$ ，可行的切割方案数是非常大的，所以，这里可以使用 delayed column generation 的思想。不妨初始化设置 $B = I_{m \times m}$ ，接下来解优化问题：

$\begin{aligned} min & 1 - c_{B}^{T} B^{- 1} A_{j} \\ s.t. & [w_{1}, \dots, w_{m}] A_{j} \leq W, a_{i j} \in N, i = 1, \dots, m \end{aligned}$

找到合适的进基变量和对应的列 $A_{j}$ . 这里的目标函数是变量 $x_{j}$ 在单纯型表中的系数。

注意到这里是一个整数规划问题，但是可以用动态规划算法以伪多项式的时间内求解。

接下来用 revised simplex 更新矩阵 $B$ . 迭代到上述问题的最优值大于等于0的时候，原问题达到最优，单纯型迭代停止。这样我们就解决了一个整数规划的线性松弛问题。

最后，对线性松弛最优解向上取证，可得： IP-optimal cost $\leq$ LP-optimal cost + $m$ . 在 $m$ 较小的时候，这是一个良好的近似解。

Cutting-plane method

列生成算法针对的是线性规划中列数特别多的情况，而切平面方法针对的是约束条件特别多的情况。这两种方法的联系可以理解为，列生成针对的是原问题，而切平面解决的是对偶问题：

$\begin{matrix} (2) & \begin{array}{cl} max & p^{T} b \\ s.t. & p^{T} A \leq c^{T} \end{array} \end{matrix}$

注意到，如果 $A_{m \times n}$ 的列数 $n$ 非常大时，上述问题的约束条件非常之多，使得单纯型法的基变量个数也非常多。

类似地，我们没有必要同时考虑所有的约束条件，如果我们考虑一个子集上的：

$\begin{array}{cl} max & p^{T} b \\ s.t. & p^{T} A_{i} \leq c_{i}, i \in I \end{array}$

计算出的最优值 $\bar{p}$ 是问题(2)的可行解，那么 $\bar{p}$ 就一定是(2)的最优解！

如果不可行，那么，求解

$min c_{i} - p^{T} A_{i}$

找到一组使 $\bar{p}$ 不可行的约束条件和对应的 $A_{i}$ ，继续迭代就可以了。因此，cutting-plane method 能否应用归结于如上问题是否能高效求解！

加入新的约束条件时，用对偶单纯型法继续迭代。

对原问题执行列生成等价于对对偶问题执行切平面！

Dantzig-Wolfe decomposition

考虑如下的线性规划：

$\begin{matrix} (3) & \begin{aligned} min & c_{1}^{'} x_{1} + c_{2}^{'} x_{2} \\ s.t. & D_{1} x_{1} + D_{2} x_{2} = b_{0} \\ F_{1} x_{1} = b_{1} \\ F_{2} x_{2} = b_{2} \\ x_{1}, x_{2} \geq 0 \end{aligned} \end{matrix}$

假设 $b_{0}, b_{1}, b_{2}$ 的维度分别是 $m_{0}, m_{1}, m_{2}$ ，对于 $m_{1}, m_{2} ≫ m_{0}$ 的情况，可以设计恰当的分解算法来高效求解。

定义多面体

$P_{i} = {x ∣ F_{i} x = b_{i}}$

于是原问题可以写成：

$\begin{matrix} (4) & \begin{aligned} min & c_{1}^{'} x_{1} + c_{2}^{'} x_{2} \\ s.t. & D_{1} x_{1} + D_{2} x_{2} = b_{0} \\ x_{1} \in P_{1}, x_{2} \in P_{2} \end{aligned} \end{matrix}$

根据 Minkowski-Weyl 定理，一个 polyhedron 可以由若干个极点和极线构成，于是可以把 $x_{1}, x_{2}$ 改写成：

$x_{i} = \sum_{j \in J_{i}} λ_{i}^{j} x_{i}^{j} + \sum_{k \in K_{i}} θ_{i}^{k} w_{i}^{k}, \sum_{j \in J_{i}} λ_{i}^{j} = 1, λ_{i}^{j}, θ_{i}^{k} \geq 0 i = 1, 2$

其中 $J_{i}, K_{i}$ 分别表示 $P_{1}, P_{2}$ 的极点集，极线集。代入原问题，可得：

$\begin{matrix} (DW-MP) & \begin{aligned} min & \sum_{j \in J_{1}} λ_{1}^{j} c_{1}^{'} x_{1}^{j} + \sum_{k \in K_{1}} θ_{1}^{k} c_{1}^{'} w_{1}^{k} + \sum_{j \in J_{2}} λ_{2}^{j} c_{2}^{'} x_{2}^{j} + \sum_{k \in K_{2}} θ_{2}^{k} c_{2}^{'} w_{2}^{k} \\ s.t. & \sum_{j \in J_{1}} λ_{1}^{j} D_{1} x_{1}^{j} + \sum_{k \in K_{1}} θ_{1}^{k} D_{1} w_{1}^{k} + \sum_{j \in J_{2}} λ_{2}^{j} D_{2} x_{2}^{j} + \sum_{k \in K_{2}} θ_{2}^{k} D_{2} w_{2}^{k} = b_{0} \\ \sum_{j \in J_{1}} λ_{1}^{j} = 1 \\ \sum_{j \in J_{2}} λ_{2}^{j} = 1 \\ λ_{i}^{j} \geq 0, θ_{i}^{k} \geq 0, \forall i, j, k . \end{aligned} \end{matrix}$

这个问题叫做 Dantzig-Wolfe master problem。注意到这个等价的问题只有 $m_{0} + 2$ 个约束条件，当 $m_{1}, m_{2}$ 比较大的时候，它很好地降低了单纯型法的存储规模！但是呢，它的列数非常大（变量个数很多），一个多面体的极点、极线的个数是阶乘级别的；幸运的是，我们能够使用列生成的思想去解决它。

首先我们没有必要同时考虑那么多个极点和极限，先只取一个很小的子集，构成一个 restricted master problem（RMP）。

设 DW-RMP 的对偶最优解是 $p = c_{B}^{T} B^{- 1} = [\begin{matrix} q & r_{1} & r_{2} \end{matrix}]$ ，变量 $λ_{1}^{j}$ 的 reduced cost 是 $(c_{1}^{'} - q^{'} D_{1}) x_{1}^{j} - r_{1}$ ，变量 $θ_{1}^{k}$ 的 reduced cost 是 $(c_{1}^{'} - q^{'} D_{1}) w_{1}^{k}$ . 接下来我们要去检验这些变量的检验数是否小于0，这归结于一个更小规模的线性规划问题：

$\begin{matrix} (DW-SP) & \begin{aligned} min & (c_{1}^{'} - q^{'} D_{1}) x_{1} \\ s.t. & x_{1} \in P_{1} \end{aligned} \end{matrix}$

这个问题叫做 Dantzig-Wolfe subproblem （DW-SP），子问题用来检验最优性！

如果最优值是 $- \infty$ ，我们能找到一条极线使得 $(c_{1}^{'} - q^{'} D_{1}) w_{i}^{k} < 0$ ，这说明 $θ_{i}^{k}$ 应该进基。

又或者最优值小于 $r_{1}$ ，也就是能找到一个极点 $x_{i}^{j}$ 使得 $(c_{1}^{'} - q^{'} D_{1}) x_{1}^{j} - r_{1} < 0$ ，这说明 $λ_{i}^{j}$ 应该进基。如果算出的最优值不小于 $r_{1}$ ，说明达到最优了~~！

对另一组问题也是类似的：

$\begin{matrix} (DW-SP) & \begin{aligned} min & (c_{2}^{'} - q^{'} D_{2}) x_{2} \\ s.t. & x_{2} \in P_{2}, \end{aligned} \end{matrix}$

列生成在 DW 分解中扮演着非常重要的作用，把求解一个大问题分解成若干个小问题，使计算得到了简化。

这个方法可以推广到多个可分离变量上：

$\begin{aligned} min & c_{1}^{'} x_{1} + c_{2}^{'} x_{2} + \dots + c_{t}^{'} x_{t} \\ s.t. & D_{1} x_{1} + D_{2} x_{2} + \dots + D_{t} x_{t} = b_{0} \\ F_{i} x_{i} = b_{i}, i = 1, 2, \dots, t, \\ x_{1}, x_{2}, \dots, x_{t} \geq 0 . \end{aligned}$

甚至当 $t = 1$ 的时候：

$\begin{aligned} min & c^{'} x \\ s.t. & D x = b_{0} \\ F x = b \\ x \geq 0 \end{aligned}$

依然可以定义 $P = {x ∣ F x = b}$ ，再利用 $P$ 的极点、极线来简化问题。当然，DW 分解并不拘泥于等式约束，只要约束条件的“较难”的一部分能够表达成 polyhedron 的形式即可。

经验上，DW 分解在迭代的一开始效果比较好，但是最优值的提升可能逐渐变慢，所以经常提前终止并输出一个次优解。

此外，在迭代的过程中，我们还能逐步更新问题的上界和下界。

DW-RMP 输出原问题的一个上界。

Stochastic programming and Benders decomposition

Benders 分解借用的是 cutting-plane 的思想。

对如下的优化问题：

$\begin{matrix} (5) & \begin{aligned} min & c^{T} x + f^{T} y \\ s.t. & A x + B y = b \\ y \geq 0, x \in X \subseteq R^{n} \end{aligned} \end{matrix}$

其中 $X$ 是一个 polyhedron。鉴于 $x, y$ 是关联（coupling）着的，问题的难度会变大。如果确定了 $x$ 能方便地计算出 $y$ ，那么就可以用 Benders 分解来使问题得到简化。

这里的 $x$ 还可以是整数变量， $y$ 是连续取值的。 $x$ 确定了，求 $y$ 就是一个简单的线性规划了。

首先，改写为关于 $x$ 的优化问题：

$\begin{matrix} (6) & \begin{aligned} min & c^{T} x + z \\ s.t. & x \in X \subseteq R^{n} \end{aligned} \end{matrix}$

式(6)称为 BD-Master Problem，其中：

$z = z (x) = [\begin{aligned} min & f^{T} y \\ s.t. & A x + B y = b, y \geq 0 \end{aligned}] = [\begin{aligned} max & α^{T} (b - A x) \\ s.t. & B^{T} α \leq f \end{aligned}]$

根据假定（原问题关于 $x$ 困难而关于 $y$ 容易），给出 $x$ ，计算 $g (x)$ 是一件轻松的事情。但是，我们还得从对偶问题来计算 $g (x)$ ，注意到右侧的对偶问题的可行域是与 $x$ 无关的。记右侧这个对偶问题为 BD-Subproblem。

记多面体 $P = {α ∣ B^{T} α \leq f}$ ，如果 $P$ 是空集，那么，由线性规划的弱对偶性，原问题无解（不可行）或者最优值无界。以下假设 $P$ 是非空的。

式(6)可以改写为：

$\begin{matrix} (7) & \begin{aligned} min & c^{T} x + z \\ s.t. & x \in X \subseteq R^{n} \\ α_{i}^{T} (b - A x) \leq z, \forall α_{i} \in J_{P} \\ α_{j}^{T} (b - A x) \leq 0, \forall α_{j} \in K_{P} \end{aligned} \end{matrix}$

$J_{P}, K_{P}$ 分别表示 $P$ 的极点、极线组成的集合。注意到 Master Problem 的约束条件非常多。式(7)的最优值和式(5)是一样的。

令 $J_{P}^{'} = K_{P}^{'} = \emptyset$ ，取一个 $x_{0}$ 开始迭代，如果对偶子问题无上界，那么这意味着原问题无解，解对偶子问题能得到一条极线 $w_{0}$ ，将 $w_{0}$ 加入到 $K_{P}^{'}$ ，这样的 $x_{0}$ 对于原始的问题(5)是不可行的！

$w_{0}$ 导致了一个 Benders feasibility cut.

现在假定对偶子问题是有界的，那么说明 $x_{0}$ 是式(5)的可行解，于是我们能得到原始问题的一个上界 $U B = c^{T} x_{0} + z (x_{0})$ ；同时，解 relaxed master problem：

$\begin{aligned} min & c^{T} x + z \\ s.t. & x \in X \subseteq R^{n} \\ α_{i}^{T} (b - A x) \leq z, \forall α_{i} \in J_{P}^{'} \\ α_{j}^{T} (b - A x) \leq 0, \forall α_{j} \in K_{P}^{'} \end{aligned}$

可以得到一个下界 $L B$ ，如果 $L B = U B$ ，那么迭代停止，输出最优解。如果 $U B \geq L B$ ，就要向 BD-RMP 添加 Benders optimality cut，在计算用对偶子问题计算 $z (x_{0})$ 的时候，得到的极点 $v_{0}$ 加入到 $J_{P}^{'}$ ，继续迭代即可。

Benders 分解解决的是一类具有特殊形式的规划问题，大规模随机规划刚好具有这样的形式。 $\begin{aligned} min & c^{T} x & + α_{1} f^{T} y_{1} & + α_{2} f^{T} y_{2} + \dots & + α_{K} f^{T} y_{K} \\ s.t. & A x & = b \\ B_{1} x & + D y_{1} & = d_{1} \\ B_{2} x & + D y_{2} & = d_{2} \\ ⋮ & ⋮ \\ B_{K} x & + D y_{K} & = d_{K} \\ x, y_{1}, \dots, y_{k} \geq 0 \end{aligned}$

Scenario 的个数 $K$ 一般很大，而第一阶段的决策 $x$ 确定下来之后，每个 $y_{k}$ 的计算是若干个小的子问题。

令：

$\begin{aligned} z_{ω} (x) = [\begin{aligned} min & f^{⊤} y_{ω} \\ s.t. & B_{ω} x + D y_{ω} = d_{ω} \\ y_{ω} \geq 0 . \end{aligned}] \end{aligned} = [\begin{aligned} max & p_{ω}^{⊤} (d_{ω} - B_{ω} x) \\ s.t. & p_{ω}^{⊤} D \leq f^{⊤} \end{aligned}]$

问题化为关于 $x$ 的：

$\begin{array}{ll} min & c^{⊤} x + \sum_{w = 1}^{K} α_{ω} z_{ω} (x) \\ s.t. & A x = b \\ x \geq 0 . \end{array}$

Summary

本文总结了求解大规模线性规划问题的算法思路。