大规模优化

为啥需要大规模优化¶

有些线性规划模型的约束条件或者变量的数量, 会随着问题规模的扩大而呈现指数级增长, 导致模型异常庞大, 难以直接求解. 本节课将学习一些专门用于处理这类大规模问题的新技术, 主要分为两种情况: 处理约束条件过多的情况和处理变量过多的情况. 除了上述的技术, 还会介绍一种名为椭球算法的, 可用于求解线性规划的替代性算法.

大规模问题¶

大规模约束条件¶

旅行商问题背景¶

首先, 我们来交代一下旅行商问题的背景. 目标函数是找到一组边, 让总权重最小, 表示为\(\sum_{e\in E}w_ex_e\), \(w_e\)是边\(e\)的权重, \(x_e\)表示\(e\)是否在路径中.

第一个约束条件是每个节点的度数为2, 也就是\(\sum_{e\in \delta(u)} x_e = 2, \forall u\in V\), 这类约束的数量等于顶点的数量\(|V|\), 是可以控制的;

第二个约束条件是子图消除约束(如果只有第一个约束, 可能会形成几个互不相关的小圈子, 例如小圈子A: 城市1->2->3->1, 小圈子B: 城市4->5->6->4, 这个满足第一个约束条件, 但是旅行商没有访问所有6个城市, 它被困在了其中的一个小圈子里面), 把任意的一部分城市圈起来(子集\(S\)), 连接这一部分城市和和其他城市的边的数量至少有两条, 也就是\(\sum_{e\in \delta(S)} x_e \geq 2, \forall \emptyset \subseteq S \subseteq V\). 这是一个复杂约束, 子集\(S\)的数量是指数级的(\(2^{|V|}\)). 一个包含60个城市的图, 其子集的数量比宇宙中的原子还多, 因此, 我们无法把这些所有的约束条件都写出来并让计算机求解.

延迟约束生成以及分离预言机¶

那么, 我们如何解决约束过多的问题呢? 这里介绍一种名为延迟约束生成(Delayed Constraint Generation)的方法. 该方法的核心思路是: 先忽略掉大部分的约束条件, 只保留一小部分, 求解这个简化后的模型, 然后检查当前解是否违反了被忽略掉的约束条件, 如果有违反的约束条件, 就把它们加回模型中, 重新求解. 重复这个过程, 直到没有违反的约束条件为止. 这个过程就像不断地用平面去切割掉不满足条件的解空间, 一步步逼近最终的正确解, 从而避免了一开始就处理海量约束的难题, 故又叫做切割平面法(Cuttinng-Plane Method).

这个策略成功的关键在于, 我们需要一个高效的方法来检查当前解是否满足约束, 如果不满足就找出一个来, 这个高效的检查工具就叫做分离预言机(Separation Oracle).

那么分离预言机是如何工作的呢? 首先, 我们要看一下, 如何高效地检查解是否有效, 一个解是否完美(即是否在解空间\(P\)中), 需要满足两种约束, 首先是度约束, 检查很简单, 因为数量不多, 一个一个核对就可以; 连通性约束, 也叫做子回路消除约束, 它们的数量是天文数字, 不可能一个一个去检查. 那么如何解决连通性约束检查难得问题呢? 这里我们会进行一个绝妙的转换: 不要去检查所有的连通性约束, 而是把它变为一个最小化的问题: 我们寻找一个子集\(S^*\), 使得连接它内外的所有边的权重之和\(\sum_{e\in \delta(S^*)} x_e\)达到最小, 这个问题可以被等价地看作是一个经典的图论问题, 全局最小割, 因为最小割问题已经被证明是可以在多项式时间内解决的, 所以我们有了一个高效的方法, 如果我们找到的这个最小割的值小于约束要求的2, 我们就成功地, 高效地找到了一个被违反地连通性约束, 这就证明了高效的分离预言是存在的.

总结一下, 在延迟约束生成这种方法中, 我们每添加一个新的约束, 之前的约束就作废了, 我们必须从头开始新的优化, 这可能导致效率低下. 一个更好的替代方案是去解它的对偶问题, 可以避免从头计算.

大规模变量¶

图的边界着色问题¶

想象一个由点和线组成的网络图, 边着色问题的目标是: 给图的每一条边染色, 要求任何两条连接到同一个点的边的颜色必须不同. 目标是使用的颜色的种类越少越好.

在这个问题中, 一种颜色对应的就是一个匹配, 匹配(Matching)是图中的一组边的集合, 在这个集合中, 没有任何两条边共享同一个端点. 所以, 用最少的颜色给边颜色的问题, 就等价于"把图中所有的边, 划分为数量最少的几个匹配".

用数学语言表示, \(\mathcal{M}\)表示这个图中所有可能的匹配的集合, \(x_M\)是一个indicator, 表示你要不要使用匹配\(M\). 目标是最小化所有被使用的匹配\(x_M\)的总和, 简单来说, 就是使用最少的匹配. 约束条件为\(\sum_{M\in \mathcal{M}:e\in M} x_M = 1\forall e\in E\), 就是说, 对于图中的任意一条边\(e\), 任意包含了这条边的匹配\(M\), 它们的权重加起来必须正好等于\(1\), 这保证了每一条边都恰好被覆盖一次, 不重复不漏下. 另一个约束是一个匹配的使用权重不能是负数.

上述线性规划问题的最大挑战在于, 一个图中所有可能的匹配的集合的数量可以是天文数字, 无法直接求解. 对于这个问题, 我们的任务是从这个包含所有可能性的集合\(\mathcal{M}\)中, 找出一个最优的子集.

列生成算法(最小化简成本预言机)¶

为了解决变量过多的问题, 我们可以使用一种名为列生成(Column Generation)的方法. 该方法的核心思路是: 先从所有可能的变量中选取一小部分, 构建一个简化的线性规划模型, 求解这个简化模型, 然后检查是否有其他被忽略的变量可以改善当前解, 如果有, 就把这些变量加回模型中, 重新求解. 重复这个过程, 直到没有更多的变量可以改善当前解为止.

首先, 我们的算法会找到一个初始可行基, 这通常很容易, 例如, 可以使用一系列只包含一条边的匹配构成的集合作为初始基. 在每个步骤中, 算法会决定引入一个新的变量, 即一个新的匹配来替换现有的基中的一个变量, 以改进当前的解.

我们如何选择那个变量呢? 取决于一个叫做判别数或者reduced cost的值. 算法会自动寻找一个新的匹配, 使得判别数为负数, 这个reduced cost公式是在前几节中就出现过的, \(c_M-c_B^TA_B^{-1}A_M\), \(c_M\)是目标函数中\(x_M\)的系数. \(A_M\)是\(x_M\)对应的列, \(A_M\)是一个由0和1组成的列下能量, 如果图中的某条边\(e\)属于我们正在考虑的这个匹配\(M\), 那么\(A_M\)在\(e\)对应的位置就是\(1\), 否则就是\(0\), 它起到了一个选择器的作用. 为了简化公式, 我们做了一个标准和聪明的代换, 定义了一个新的变量\(y\), \(y^T = c_B^TA_B^{-1}\), \(y\)是对偶变量, 把\(y^T\)代入公式, 并且\(c_M=1\), 判别数就变成了\(1-y^TA_M\). 现在我们来分析后面这一项\(y^TA_M\), \(y\)是一个于图中所有边相关的向量, 每一条边\(e\)都由一个对应的\(y_e\)值, \(A_M\)是一个只在匹配\(M\)包含的边\(e\)的位置上取值为\(1\)的向量. 所以, 这两者的乘积\(y^TA_M\), 实际上就是在计算, 把匹配\(M\)中所有边\(y_e\)值加起来的总和, 因此, \(y^TA_M=\sum_{e\in M} y_e\), 所以, reduced cost就变成了\(1-\sum_{e\in M} y_e\), 在所有判别数为负数的匹配中, 我们想要的是那个使得判别数最小(负得最厉害)的匹配, 因为它能够最大化程度的优化我们的解, 要让\(1-\sum_{e\in M} y_e\)尽可能小, 就必须让\(\sum_{e\in M} y_e\)尽可能大. 这就把我们的问题转换成了一个经典的图论问题, 最大权匹配(Maximum Weight Matching), 因为我们想要找到一个匹配\(M\), 使得它包含的边的\(y_e\)值之和最大. 最大权匹配问题已经被证明可以在多项式时间内解决, 所以我们有了一个高效的方法来选择新的变量. 这个算法的输出 -- 那个具体的匹配\(M\), 是被创造或者寻找出来的, 而不是从一个现成的清单里面选出来的, 这是和普通的单纯性法的不同之处.

椭球算法¶

经过上述的小节, 我们知道, 在解决变量或约束几多的线性规划问题的时候, 我们不必一开始就写出所有的变量或约束. 我我们可以借助预言机来运行单纯形法. 分为分离预言机: 给出一个点, 它能判断这个点是否在可行解区域内, 如果不在, 它会给出一个将该点和可行域分开的超平面(即一个新的约束); 最小化简成本预言机: 帮助找到能最有效改进当前解的变量.

虽然预言机方法很巧妙, 但是在最坏的情况下, 它可能会不断生成大量的约束或者变量, 导致算法效率低下, 依然无法保证高效求解. 为了解决这个问题, 我们接下来要引出一种算法: 椭球算法.

椭球的定义¶

首先, 我们需要用数学语言给椭球下一个精确的定义. 简单的来说, 它的定义分为两步:

基础形状: 单位球. \(B^n=\{x\in \mathbb{R}^n: x^Tx\leq 1\}\). 这代表在\(n\)维空间中, 所有到源点距离小于等于\(1\)的点的集合, 想象一个以原点为中心, 半径为\(1\)的完美球体.
仿射变换: 公式为\(f(x)=Mx+s\), 这个变换做了两件事情, \(Mx\)(拉伸/旋转), \(M\)是一个矩阵, 负责把单位球沿着不同方向进行拉伸, 挤压或者旋转, 让它从一个完美的球体变为一个椭的形状. \(+s\)表示平移, \(s\)是一个向量, 负责把这个变形后的椭球中心挪到空间中的另一个位置\(s\).

算法¶

椭球算法的最基本功能不是直接求解线性规划的最优解, 而是判断一个由不等式组\(Ax\leq b\)定义的多面体区域\(P\)是否可行. 也就是说, 它只回答一个问题, 这个区域\(P\)是空的, 还是至少包含一个点.

算法的思路非常直观, 从一个巨大的, 已知的能完全包住目标区域\(P\)的椭球开始, 不断找到一个更小的椭球\(E'\), 它能够同样包住目标区域, 重复这个过程, 直到椭球缩小到足够小, 我们就能判断出\(P\)是否为空. 算法得以成立的三个前提:

前提1: 我们可以找到一个足够大的初始球, 将这个目标区域\(P\)完全包裹起球来, 这是算法的起点
前提2: 如果目标区域\(P\)不是空的, 那么它就有一定的体积, 我们至少可以在里面放一个半径为\(\epsilon\)的球, 这保证了如果解存在, 我们最终不会错过它
前提3: 当我们一个椭球\(E\)的时候, 我们可以检查它的中心点\(s\)是否在目标区域\(P\)内, 如果中心点\(s\)不在\(P\)里面, 我们就能找到一个超平面\(H\)将\(s\)和\(P\)分开. 最关键的是, 我们能够高效地计算出一个新的, 更小的椭球\(E'\), 这个新椭球只包裹旧椭球被切掉一般之后的有效部分\(E\cap H\), 这个新椭球\(E'\)的体积会比旧椭球\(E\)的体积显著减小. 最终要么找到一个可行点, 要么证明这个区域是空的(通过判定是否超时). 这个体积缩减的速度满足\(\text{volume}(E')<\frac{\text{volume}(E)}{e^{\frac{1}{2n+2}}}\).

正确性证明¶

Theorem 8.1说明了椭球算法的正确性. 这是用反证法来证明的. 我们线假设一种情况: 算法运行了足够长的时间后, 报告说\(P\)是空的, 但实际上\(P\)不是空的, 然后我们从这个错误的假设出发进行推导, 最终得到一个荒谬的, 自相矛盾的结论, 这个矛盾证明了我们最初的假设是错误的. 因此, \(P\)不可能不是空的, 它只能是空的.

假设算法运行了\(k=n(2n+2)\ln (R/\epsilon)\)次迭代后超时, 并返回\(P\)是空的. 同时, 我们假设\(P\)实际上不是空的.

事实1: 算法从一个包含\(P\)的大椭球\(E_0\)开始, 每一次迭代, 我们都只切掉不包含\(P\)的部分, 所以\(P\)始终被包含在每一个后续的椭球\(E_i\)中. 因此, 当算法结束的时候, 非空的\(P\)必定完整地存在于最后一个, 也是最小地椭球\(E_k\)地内部, 即\(P\subseteq E_k\), 既然\(P\)在\(E_k\)内部, 那么\(P\)的体积不可能大于\(E_k\)的体积, 有\(\text{volume}(P)\leq \text{volume}(E_k)\).

事实2: 根据第二个前提, 如果\(P\)不是空的, 那么它的体积至少有\(\text{volume}(B^n(\epsilon))\)这个大. 而根据椭球算法的收缩保证, 经过\(k\)次迭代后, 椭球\(E_k\)的体积已经变得非常小了. 中间的数学推到表明: \(\text{volume}(E_k)<\text{volume}(B^n)(\epsilon)\)(见note推导).

根据事实1和事实2. \(\text{volume}(P)<\text{volume}(P)\)是不成立的, 所以是一个矛盾. 因此, 我们的假设是错误的, 也就是说, 如果算法报告\(P\)是空的, 那么\(P\)确实是空的.

前面的算法有一些简化的前提, 是故意为了让证明过程变得简单移动, 在实际应用中要移除这些理想化的前提, 需要非常复杂和繁琐的数学技巧, 这些技巧已经超出了这门课的范围, 所以不进行深入探讨. 最重要的结论是, 尽管技术细节复杂, 但是最终的结论非常重要: 椭球算法是一个高效的算法, 它的运行时间是多项式时间, 这个时间复杂度的增加取决于三个维度, 变量的数量, 约束条件的数量和所有约束条件中的数字所需的总数据量.