不稳定性和优化问题的几何#

正如我们所看到的，每当我们用数值方法解决一个问题时，我们必须接受我们提供的输入和获得的输出可能与给定问题的理论或数学解有所不同。例如，\(0.1\)在计算机中将被表示为一个与\(0.1\)相差约\(10^{-17}\)的数字。因此，一个自然的担忧是这些微小的差异是否会导致计算解中的大差异。

这是给定问题的条件数概念背后的思想。虽然对于大多数实际优化问题来说，输入中的小扰动只会引起问题最终答案的小扰动，但在某些特殊情况下，情况并非如此。这些表现不佳的问题被称为病态或数值不稳定。

本节旨在展示，在线性优化问题的背景下，这种行为的最常见来源，以及如何完全避免这种行为。我们将首先回顾具有唯一解的线性系统的求解问题，然后进入线性优化问题的更核心问题，其几何解释，然后描述一些最常见的不良情况。我们随后提供了两个带有互动材料的思维实验，以帮助说明本节的概念。我们以关于这个主题的一些进一步思考作为结束。

请注意，尽管条件数的概念在学术界受到了广泛关注，但回顾这些文献超出了本文档的范围。如果您想开始研究这个主题，可以在维基百科的条件数页面上找到一个好的切入点。

线性系统的情况#

求解线性系统是任何MI(QC)P求解器中非常常见的子程序，因为在算法的整个执行过程中我们必须解决许多线性系统。

因此，考虑我们有一个线性系统 \(Ax=b\) 具有唯一解（即 \(A\) 是一个满秩的方阵），并且您希望评估如果我们扰动右侧 \(b\)，系统的解可能会如何变化。由于系统有唯一解，我们知道给定 \(b\)，解将是 \(A^{-1}b\)，如果我们用 \(\varepsilon\) 扰动 \(b\)，解将是 \(A^{-1}(b+\varepsilon)\)。解相对于输入的相对变化的度量将是比值

\[\eta(b,\varepsilon):=\frac{\|A^{-1}b\|}{\|A^{-1}(b+\varepsilon)\|}/\frac{\|b\|}{\|b+\varepsilon\|}.\]

请注意，上述定义与\(b\)和\(\varepsilon\)的大小无关。从那里，最差可能的比率将是以下结果

\[\kappa(A):=\max_{b,\varepsilon}\eta(b,\varepsilon).\]

这个量被称为矩阵 \(A\) 的条件数，并定义为

\[\kappa(A)={\|A\|}{\|A^{-1}\|}.\]

对\(\kappa(A)=10^k\)的一个常见解释是，在求解系统\(Ax=b\)时，你可能会在\(x\)中失去\(b\)中精度的\(k\)位数字。

LP中最佳单纯形基的条件数在KappaExact属性中捕获。非常大的\(\kappa\)值可能表明结果可能不稳定。

当这种情况确实发生时，最好的建议是缩放约束矩阵的系数，以使最终的系数范围较小。这种转换通常会减少最终基的\(\kappa\)值；请参阅容差和用户缩放部分，了解如何执行此重新缩放，以及关于缩放的一般注意事项。

线性优化问题的几何#

在展示表现出不良行为的优化模型之前，我们首先需要理解它们背后的几何。考虑一个形式的问题

\[\begin{split}\begin{array}{ll} \max & cx\\ s.t. & Ax \leq b.\\ \end{array}\end{split}\]

例如：

\[\begin{split}\begin{array}{lrrl} \max & x + y & \vec{c}=&(1,1)\\ s.t. & -x \leq 0 & A_{1\cdot}=&(-1,0)\\ & x \leq 1 & A_{2\cdot}=&(1,0)\\ & -y \leq 0 & A_{3\cdot}=&(0,-1)\\ & y \leq 1 & A_{4\cdot}=&(0,1).\\ \end{array}\end{split}\]

注意，如果我们表示 \(b^t:=(0,1,0,1)\)，那么问题可以表述为

\[\max_{x\in\mathbb{R}^2}\{ \vec{c}x:Ax\leq b\}.\]

可行区域、改进方向 \(\vec{c}\) 和最优解 \(x^*\) 可以表示为

../../_images/codedraw3.svg
../../_images/codedraw3-dark.svg

请注意，每当我们沿着\(\vec{c}\)的方向移动时，值\(\vec{c}x\)会增加。此外，由于我们无法从\(x^*\)移动到另一个具有更好目标值的可行点，我们可以得出结论，\(x^*\)确实是该问题的最优解。请注意，\(x^*\)是可行区域的一个角点。这不是巧合；如果可行区域有界且\(\vec{c}\)不为零，你总是会在角点找到最优解。如果目标为零，则所有可行解都是最优的；我们稍后会详细讨论零目标及其含义。

要理解输入数据的变化如何影响可行区域和最优解，考虑一个小修改： \(\tilde{b}^t=(\varepsilon,1,0,1)\), \(\tilde{\vec{c}}=(1+\varepsilon,1)\), 和 \(\tilde{A_{4\cdot}}=(\varepsilon,1)\)。那么我们的优化问题将看起来像

../../_images/codedraw4.svg
../../_images/codedraw4-dark.svg

请注意，尽管我们改变了右侧，这一改变对问题的最优解没有影响，但它确实通过扩大可行区域的底部部分改变了可行区域。

改变目标向量会倾斜图形表示中的相应向量。这当然也会改变最优目标值。扰动约束会倾斜约束的图形表示。\(A_{4\cdot}\)的变化会改变原始解本身。约束的倾斜程度取决于扰动的相对值。例如，尽管约束\(x \leq 1\)和约束\(100 x \leq 100\)诱导相同的可行区域，但扰动\(x + \varepsilon y\leq 1\)会比扰动\(100 x + \varepsilon y \leq 100\)诱导更多的倾斜。

多重最优解#

优化初学者中常见的一个误解是认为优化问题实际上只有一个解决方案。令人惊讶的是，这通常是不正确的。对于许多实际问题，目标（无论是成本还是收入或……）主要由少数几个变量决定，而大多数变量只是确保解决方案的实际操作是可能的。例如，考虑一个人员配置问题，成本通常由在给定工作日工作的人数决定，而不是由具体的人员决定。

这种情况自然会导致类似于

\[\begin{split}\begin{array}{lrrl} \max & y & \vec{c}=&(0,1)\\ s.t. & -x \leq 0 & A_{1\cdot}=&(-1,0)\\ & x \leq 1 & A_{2\cdot}=&(1,0)\\ & -y \leq 0 & A_{3\cdot}=&(0,-1)\\ & y \leq 1 & A_{4\cdot}=&(0,1).\\ \end{array}\end{split}\]

图形上可以表示为

../../_images/codedraw5.svg
../../_images/codedraw5-dark.svg

在这种情况下，很明显\(x^1\)、\(x^3\)以及位于这两点之间的所有解都是最优的。单纯形算法将返回\(x^1\)或\(x^3\)（如果更改参数，可能会切换）。障碍算法（无交叉）将返回\(x^2\)。这些解都是正确的；问题本身没有理由偏好其中一个。如果你有偏好，你需要在目标函数中明确说明。

处理Epsilon最优解#

上一节考虑了多个（真实）最优解的情况。当我们有几个\(\varepsilon\)-最优解时会发生什么？更具体地说，考虑

\[\begin{split}\begin{array}{lrrl} \max & \varepsilon x + y & \vec{c}=&(\varepsilon,1)\\ s.t. & -x \leq 0 & A_{1\cdot}=&(-1,0)\\ & x \leq 1 & A_{2\cdot}=&(1,0)\\ & -y \leq 0 & A_{3\cdot}=&(0,-1)\\ & y \leq 1 & A_{4\cdot}=&(0,1).\\ \end{array}\end{split}\]

图形上可以表示为

../../_images/codedraw6.svg
../../_images/codedraw6-dark.svg

如果 \(\varepsilon\) 为零，那么我们处于之前描述的情况。然而，请注意，目标向量的微小扰动可能会导致 \(x^1\) 或 \(x^2\) 被报告为最优解。而容差在这里可能起到很大的作用。例如，如果 \(\varepsilon\) 为负，那么 \(x^1\) 将是数学上的最优结果，但由于 最优性容差，单纯形法可能会得出 \(x^2\) 是最优的结论。更准确地说，如果 \(\varepsilon\) 小于默认的最优性容差 \(10^{-6}\)，那么单纯形法可以自由地宣布任一解为最优解（在容差范围内）。

每当\(x^1\)和\(x^2\)之间的距离不太大时，上述陈述是正确的。为了理解这一点，考虑当我们将\(A_{4\cdot}\)的右侧从1更改为\(10^6\)时会发生什么。那么可行区域将是一个非常长的矩形盒子，顶点为\((0,0)\)、\((0,1)\)、\((10^6,1)\)和\((10^6,0)\)。可能有些令人惊讶的是，如果\(\varepsilon\)低于对偶容差，单纯形法可能会认为\((10^6,1)\)是最优的，即使其目标值为\(1-10^6\varepsilon\)，这在最终目标值方面可能非常相关。

注意，这两种情况都有一个共同点：目标函数几乎与可行区域的一边平行。在第一种情况下，这一边相对较短，因此从\(x^2\)跳到\(x^1\)会导致目标值的微小变化。在第二种情况下，几乎与目标函数平行的一边非常长，现在从\(x^2\)跳到\(x^1\)可能会对最终的目标函数产生重大影响。

如果你去掉这两个因素中的任何一个，即目标向量几乎与约束平行，或者由这个几乎平行的约束引起的边缘非常长，那么这个问题就不会出现。由于本节开头讨论的原因，目标函数通常接近与一个或多个约束平行。因此，避免这种情况的最佳方法是避免第二个条件。最简单的方法是确保变量的范围不要太大。请参考容差和用户缩放部分以获取相关指导。

薄可行区域#

我们现在考虑另一种可能导致意外结果的极端情况。考虑定义为

\[\begin{split}\begin{array}{lrrl} \max & y & \vec{c}=&(0,1)\\ s.t. & -x + \varepsilon y \leq 1 & A_{1\cdot}=&(-1,\varepsilon)\\ & x + \varepsilon y \leq 1 & A_{2\cdot}=&(1,\varepsilon)\\ & -y \leq 0 & A_{3\cdot}=&(0,-1)\\ \end{array}\end{split}\]

及其图形表示

../../_images/codedraw7.svg
../../_images/codedraw7-dark.svg

从图形表示中可以清楚地看出，问题的最优解将在约束条件\(A_{1\cdot}\)与\(A_{2\cdot}\)的交点处；如果我们进行代数运算，我们将得到\(x^*=(0,\frac{1}{\varepsilon})\)。还要注意，随着\(\varepsilon\)的减小，可行区域向上延伸，而其底部保持不变。我们将考虑\(\varepsilon\)是一个非常小的正数（在\(10^{-9}\)到\(10^{-6}\)之间）的情况。

如果我们扰动右侧向量 \(b\) 从 \((1,1)\) 到 \((1+\delta,1)\)，新的解将是 \(\tilde{x}^*=(-\frac{\delta}{2},\frac{2+\delta}{2\varepsilon})\)。 为了评估这种扰动的影响，我们计算修改后的解与之前解之间的 \(L_1\) 距离，其由以下公式给出：

\[\|x^*-\tilde{x}^*\|_1 = \frac{|\delta|}{2}+\frac{|\delta|}{\varepsilon}\]

这个量可以很小也可以非常大，取决于\(\delta\)和\(\varepsilon\)之间的相对大小。如果\(\delta\)远小于\(\varepsilon\)，那么这个量将会很小。然而，如果\(\delta\)大于或甚至与\(\varepsilon\)同数量级，情况将会相反。输入数据中的非常小的扰动可能会导致最优解中的大变化。

如果我们扰动\(A_{1\cdot}\)到\((-1,\delta)\)，也会出现类似的问题；新的最优解变为\(\tilde{x}^*=(1-\frac{2\varepsilon}{\varepsilon+\delta},\frac{2}{\varepsilon+\delta})\)。但现在，如果\(\delta=\varepsilon/2\)，那么\(y\)的新解将从\(\frac{1}{\varepsilon}\)变为\(\frac{4}{3\varepsilon}\)（相对差异为33%）。同样，输入的微小变化可能会导致最优解的巨大变化。

是什么导致了这种不良行为？问题在于最优解是由两个几乎平行的约束条件定义的。当\(\varepsilon\)越小，这两个约束条件就越接近平行。当约束条件如此接近平行时，斜率的微小变化可能导致它们交点的巨大移动。从数学上讲：

\[\lim_{\varepsilon\rightarrow0^+}\|x^*\|=\infty\]

然而请注意，如果原始问题有一个额外的变量约束形式为 \(y \leq 10^4\)，那么这些不良行为都不会发生。对于任何小于 \(10^{-4}\) 的 \(\varepsilon\) 值，最优解将由新约束和其中一个约束 \(A_{2\cdot}\) 或 \(A_{1\cdot}\) 定义，这将再次导致行为良好（即稳定）的解。总之，这种问题只有在可行区域无界或非常大时才会出现。有关限制可行区域的进一步指导，请参阅 容差和用户缩放 部分。

优化圆上的问题#

现在我们提供我们的第一个思维实验：考虑在由约束定义的可行区域上优化线性函数的问题

\[\sin(2\pi\frac{i}{10^6}) x + \cos(2\pi\frac{i}{10^6}) y \leq 1,\,\forall i\in\{1,\ldots,10^6\},\]

即可行域本质上是\(\mathbb{R}^2\)中的一个单位圆。请注意，对于所有目标函数，相应的最优点将由两个非常接近平行的线性约束定义。问题的数值解会发生什么？你能猜到吗？这种情况如下图所示：

../../_images/codedraw1.svg
../../_images/codedraw1-dark.svg

为了进行实验，我们执行代码 circleOpt.py，在其中我们随机选择一个目标 向量，找到结果优化问题的最优解， 并计算几个相关的量：

  • 在所有之前的运行中，报告的原解与实际优化完美圆问题的理论解之间的最差距离。

  • Gurobi 在所有先前运行中报告的最严重的边界违规。

  • Gurobi 在所有先前运行中报告的最严重的约束违反情况。

  • Gurobi 在所有先前运行中报告的最差对偶违反。

  • 之前的实验数量。

  • 累计的单纯形迭代次数。

  • 当前最优基的 \(\kappa\) (KappaExact 属性) 值。

示例输出如下所示：

Added 2 Vars and 1048576 constraints in 19.19 seconds
Errors: 8.65535e-08 0 2.94137e-07 2.77556e-17 Iter 0 10 Kappa 3150.06
Errors: 4.81978e-07 0 3.22359e-07 2.77556e-17 Iter 1 21 Kappa 3009.12
Errors: 4.81978e-07 0 3.4936e-07 1.11022e-16 Iter 2 33 Kappa 2890.58
Errors: 1.53201e-06 0 9.78818e-07 1.11022e-16 Iter 6 79 Kappa 1727.89
Errors: 1.61065e-06 0 8.26005e-07 1.11022e-16 Iter 46 536 Kappa 1880.73
Errors: 1.61065e-06 0 8.84782e-07 1.11022e-16 Iter 52 602 Kappa 1817.27
Errors: 1.61065e-06 0 9.4557e-07 1.11022e-16 Iter 54 625 Kappa 1757.96
Errors: 1.69167e-06 0 9.78818e-07 1.11022e-16 Iter 64 742 Kappa 1727.89
Errors: 1.69167e-06 0 3.8268e-07 1.66533e-16 Iter 88 1022 Kappa 2761.99
Errors: 1.69167e-06 0 9.04817e-07 1.66533e-16 Iter 92 1067 Kappa 1797.06
Errors: 1.69167e-06 0 2.94137e-07 2.22045e-16 Iter 94 1089 Kappa 3150.06
Errors: 1.69167e-06 0 3.29612e-07 2.22045e-16 Iter 95 1101 Kappa 2975.84
Errors: 1.69167e-06 0 3.4936e-07 2.22045e-16 Iter 98 1137 Kappa 2890.58
Errors: 1.69167e-06 0 9.25086e-07 2.22045e-16 Iter 99 1147 Kappa 1777.3
Errors: 1.69167e-06 0 9.78818e-07 2.22045e-16 Iter 107 1237 Kappa 1727.89
Errors: 1.69167e-06 0 9.99895e-07 2.22045e-16 Iter 112 1293 Kappa 1709.61
Errors: 1.84851e-06 0 9.78818e-07 2.22045e-16 Iter 132 1523 Kappa 1727.89
Errors: 1.96603e-06 0 9.99895e-07 2.22045e-16 Iter 134 1545 Kappa 1709.61

令人惊讶的是，报告的错误相当小。为什么会这样呢？至少有两个因素在起作用：模型有一个有界的可行区域（在这种情况下，两个变量的范围都是\([-1,1]\)）。此外，从一个极值点（两个相邻约束的交点）到其邻居的距离也相对较小，因此所有\(\varepsilon\)-最优解都彼此接近。

我们鼓励您尝试这段代码，扰动一些输入数据，并分析结果。您将看到理论与数值最优解之间的差异将与扰动的大小相当。

优化薄区域#

现在我们转向第二个思想实验：考虑一个由三角形组成的可行区域，该三角形在\(\mathbb{R}^2\)中具有非常宽的底边和非常短的高度，如下图所示：

../../_images/codedraw2.svg
../../_images/codedraw2-dark.svg

考虑底与高之比约为\(10^5\)的情况，并且我们考虑一个如图所示的名义目标函数\(\vec{c}_1\)。

理论上，最优解应该是三角形的顶点，但假设我们随机扰动右侧和目标函数，扰动量在\(10^{-6}\)的数量级。数值解会发生什么变化？

为了进行实验，我们执行代码 thinOpt.py，在其中我们按照上述描述进行一系列不同扰动的重新优化。更准确地说，每当新计算的解比之前的试验中离数学解更远时，我们打印：

  • 解决方案之间的新最大距离。

  • 当前迭代。

  • 当前最优基的 \(\kappa\) (KappaExact 属性) 值。

  • Gurobi报告的当前解的边界违规情况。

  • Gurobi 报告的当前解的约束违规情况。

  • Gurobi 报告的当前解的对偶违反。

示例输出如下所示：

New maxdiff 4e+16 Iter 0 Kappa 3.31072 Violations: 0 0 0
New maxdiff 4e+16 Iter 1 Kappa 3.31072 Violations: 0 0 0
New maxdiff 4e+16 Iter 2 Kappa 3.31072 Violations: 0 0 0
New maxdiff 4e+16 Iter 7 Kappa 3.31072 Violations: 0 0 0
New maxdiff 4e+16 Iter 83 Kappa 3.31072 Violations: 0 0 2.64698e-23
New maxdiff 4e+16 Iter 194 Kappa 3.31072 Violations: 0 0 0
New maxdiff 4e+16 Iter 1073 Kappa 3.31072 Violations: 0 1.13687e-13 0
New maxdiff 4e+16 Iter 4981 Kappa 3.31072 Violations: 0 0 0
New maxdiff 4e+16 Iter 19514 Kappa 3.31072 Violations: 0 0 0
New maxdiff 4e+16 Iter 47117 Kappa 3.31072 Violations: 0 0 0
New maxdiff 4e+16 Iter 429955 Kappa 3.31072 Violations: 0 0 0
New maxdiff 4e+16 Iter 852480 Kappa 3.31072 Violations: 0 0 0

结果与我们第一次测试中看到的非常不同。未扰动模型的解与扰动模型的解之间的距离非常大，即使从第一次迭代开始也是如此。此外，\(\kappa\) 值相对较小，报告的原问题、对偶问题和边界违规几乎为零。那么，发生了什么？请注意，当我们选择 \(\vec{c}_1=(0,1)\) 时，我们选择了一个最优解点，目标函数的一个小倾斜可能会将我们移动到另一个非常远的极值点，因此范数很大。这是可能的，因为区域非常大，原则上没有任何边界，即这与 \(\varepsilon\)-最优解和非常长的边的情况有关。

再次，我们鼓励您尝试这个例子。例如，如果名义目标函数是 \(\vec{c}_2=(1,0)\)，会发生什么？

稳定性和收敛性#

用于解决线性规划问题的算法都不得不做出一个假设：系统的小变化（例如，在障碍中迈出一小步）会导致解决方案的小变化。如果这不是真的（由于病态条件），那么算法可能会在解空间中跳跃，并且难以收敛。

最后，改进问题几何形状的一种方法是适当地缩放变量和约束，如容差和用户缩放部分所述，并使用所有变量具有合理范围的有界可行集。