二次规划算法

二次规划定义

二次规划是指这样的问题：找到向量 x，它能在满足特定线性约束的条件下最小化二次函数：

$$\underset{x}{\min }\frac{1}{2}{x}^{T}Hx+c^{T}x$$

(1)

满足 A·x ≤ b, Aeq·x = beq, l ≤ x ≤ u。

interior-point-convex quadprog 算法

interior-point-convex 算法执行以下步骤：

• 预求解/后求解

• 生成初始点

• 预测-校正

• 停止条件

• 不可行性检测

注意

算法有两种代码路径。当黑塞矩阵 H 是由双精度值组成的普通（满）矩阵时，它采用一种代码路径；当 H 是稀疏矩阵时，它采用另一种代码路径。有关稀疏数据类型的详细信息，请参阅稀疏矩阵。通常，对于非零项相对较少的大型问题，该算法在 H 指定为 sparse 时执行更快。同样，对于较小或相对稠密的问题，该算法在 H 指定为 full 时执行更快。

预求解/后求解

该算法首先尝试通过消除冗余和简化约束来简化问题。在预求解步骤中执行的任务可以包括以下内容：

• 检查是否有任何变量具有相等的上界和下界。如果有，请检查可行性，然后固定变量的值以去除变量。

• 检查是否有任何线性不等式约束只涉及一个变量。如果是，请检查可行性，然后将线性约束更改为一个边界。

• 检查是否有任何线性等式约束只涉及一个变量。如果是，请检查可行性，然后固定该变量的值以去除该变量。

• 检查是否有任何线性约束矩阵具有全零行。如果是，请检查可行性，然后删除这些行。

• 确定边界和线性约束是否一致。

• 检查是否有任何变量仅作为线性项出现在目标函数中，而不出现在任何线性约束中。如果是这样，请检查可行性和有界性，然后将变量固定为其适当的边界值。

• 通过添加松弛变量，将任何线性不等式约束更改为线性等式约束。

如果算法检测到不可行或无界的问题，它会停止并发出相应的退出消息。

该算法可能到达一个表示解的单个可行点。

如果算法在预求解步骤中没有检测到不可行或无界的问题，并且预求解没有产生解，则算法会继续其后续步骤。到达停止条件后，算法会重新构造原始问题，这将撤消任何预求解变换。这最后一步是后求解步骤。

有关详细信息，请参阅 Gould 和 Toint 合著的 [63]。

生成初始点

该算法的初始点 x0 是：

1. 将 x0 初始化为 ones(n,1)，其中 n 是 H 中的行数。

2. 对于同时具有上界 ub 和下界 lb 的分量，如果 x0 的分量未严格位于边界之内，则分量将设置为 (ub + lb)/2。

3. 对于只有一个边界的分量，如有必要，请修改分量，使其严格位于边界内。

4. 执行预测步（请参阅预测-校正），进行微小校正以实现可行性，而不是执行完整的预测-校正步。这将使初始点更靠近中心路径，且不必花费完整预测-校正步的开销。有关中心路径的详细信息，请参阅 Nocedal 和 Wright 的论著 [7]（第 397 页）。

预测-校正

稀疏内点凸算法和满内点凸算法的主要区别在预测-校正阶段。两种算法相似，但在某些细节上有所不同。有关基本算法说明，请参阅 Mehrotra 的 [47]。

算法首先通过将 A 和 b 乘以 -1，将线性不等式 Ax <= b 转换为 Ax >= b 形式的不等式。这不影响解，但可使其形式与一些文献中的形式保持一致。

• 稀疏预测-校正算法

• 满预测-校正算法

稀疏预测-校正算法.  与 fmincon 内点算法相似，稀疏 interior-point-convex 算法尝试找到 卡鲁什-库恩-塔克 (KKT) 条件成立的点。对于二次规划定义中所述的二次规划问题，这些条件是：

$$\begin{array}{c}Hx+c-A_{eq}^{T}y-{\overline{A}}^{T}z=0\\ \overline{A}x-\overline{b}-s=0\\ A_{eq}x-b_{eq}=0\\ s_i{z}_i=0,i=1,2,\mathrm{...},m\\ s\ge 0\\ z\ge 0.\end{array}$$

此处

• $$\overline{A}$$ 是扩展的线性不等式矩阵，其中包括以线性不等式形式编写的边界。$$\overline{b}$$ 是对应的线性不等式向量，包括边界。

• s 是将不等式约束转换为等式约束的由松弛变量组成的向量。s 的长度为 m（即线性不等式和边界的数目）。

• z 是对应于 s 的拉格朗日乘数组成的向量。

• y 是与等式约束相关联的拉格朗日乘数组成的向量。

该算法首先根据牛顿-拉夫逊公式计算一个预测步，然后再计算一个校正步。校正尝试以更好的方式实现非线性约束 s_iz_i = 0。

预测步的定义：

• r_d，对偶残差：

  $$r_d=Hx+c-A_{eq}^{T}y-{\overline{A}}^{T}z.$$

• r_eq，原始等式约束残差：

  $$r_{eq}=A_{eq}x-b_{eq}.$$

• r_ineq，原始不等式约束残差，包括边界和松弛变量：

  $$r_{ineq}=\overline{A}x-\overline{b}-s.$$

• r_sz，互补残差：

  r_sz = Sz。

  S 是松弛项组成的对角矩阵，z 是拉格朗日乘数组成的列矩阵。

• r_c，平均互补：

  $$r_c=\frac{s^{T}z}{m}.$$

在牛顿步中，x、s、y 和 z 中的变化由下式给出：

$$\left(\begin{array}{cccc}H& 0& -A_{eq}^T& -{\overline{A}}^T\\ A_{eq}& 0& 0& 0\\ \overline{A}& -I& 0& 0\\ 0& Z& 0& S\end{array}\right)\left(\begin{array}{c}\Delta x\\ \Delta s\\ \Delta y\\ \Delta z\end{array}\right)=-\left(\begin{array}{c}{r}_d\\ r_{eq}\\ r_{ineq}\\ r_{sz}\end{array}\right).$$

(2)

然而，完整的牛顿步可能不可行，因为存在对 s 和 z 的正性约束。因此，quadprog 会根据需要缩短步长以保持正性。

此外，为了在内部保持“中心化”位置，该算法不尝试求解 s_iz_i = 0，而是采用正参数 σ，并尝试求解

s_iz_i = σr_c。

quadprog 用 r_sz + ΔsΔz – σr_c1 取代牛顿步方程中的 r_sz，其中 1 是由 1 组成的向量。此外，quadprog 对牛顿方程进行重新排序，以获得对称的、数值更稳定的方程组来计算预测步。

在计算校正的牛顿步后，算法执行更多计算，以获得更长的当前步，并准备好计算更优的后续步。这些多重校正计算可以提高性能和稳健性。关详细信息，请参阅 Gondzio 的论著 [4]。

满预测-校正算法.  满预测-校正算法不会将边界合并为线性约束，因此它有另一组对应于边界的松弛变量。该算法将下界移至零。如果某个变量只有一个边界，则算法将上界不等式反向，进而将其转换为以零为下界。

$$\overline{A}$$ 是同时包含线性不等式和线性等式的扩展线性矩阵。$$\overline{b}$$ 是对应的线性等式向量。$$\overline{A}$$ 还包括用于扩展向量 x 的项，这些项使用松弛变量 s 将不等式约束转换为等式约束：

$$\overline{A}x=\left(\begin{array}{cc}{A}_{eq}& 0\\ A& I\end{array}\right)\left(\begin{array}{c}{x}_0\\ s\end{array}\right),$$

其中，x₀ 表示原始 x 向量。

KKT 条件是

$$\begin{array}{c}Hx+c-{\overline{A}}^{T}y-v+w=0\\ \overline{A}x=\overline{b}\\ x+t=u\\ v_i{x}_i=0,i=1,2,\mathrm{...},m\\ w_i{t}_i=0,i=1,2,\mathrm{...},n\\ x,v,w,t\ge 0.\end{array}$$

(3)

为了找到公式 3 的解 x、松弛变量和对偶变量，该算法基本上考虑牛顿-拉夫逊步：

$$\left(\begin{array}{ccccc}H& -{\overline{A}}^T& 0& -I& I\\ \overline{A}& 0& 0& 0& 0\\ -I& 0& -I& 0& 0\\ V& 0& 0& X& 0\\ 0& 0& W& 0& T\end{array}\right)\left(\begin{array}{c}\Delta x\\ \Delta y\\ \Delta t\\ \begin{array}{l}\Delta v\\ \Delta w\end{array}\end{array}\right)=-\left(\begin{array}{c}Hx+c-{\overline{A}}^{T}y-v+w\\ \overline{A}x-\overline{b}\\ u-x-t\\ \begin{array}{l}VX\\ WT\end{array}\end{array}\right)=-\left(\begin{array}{c}{r}_d\\ r_p\\ r_{ub}\\ \begin{array}{l}{r}_{vx}\\ r_{wt}\end{array}\end{array}\right),$$

(4)

其中，X、V、W、T 是分别对应于向量 x、v、w 和 t 的对角矩阵。方程最右边的各个残差向量分别是：

• r_d，对偶残差

• r_p，原始残差

• r_ub，上界残差

• r_vx，下界互补残差

• r_wt，上界互补残差

该算法以如下方式来求解 公式 4：先将其转换为以下对称矩阵形式

$$\left(\begin{array}{cc}-D& {\overline{A}}^T\\ \overline{A}& 0\end{array}\right)\left(\begin{array}{c}\Delta x\\ \Delta y\end{array}\right)=-\left(\begin{array}{c}R\\ r_p\end{array}\right),$$

(5)

其中

$$\begin{array}{l}D=H+X^{-1}V+T^{-1}W\\ R=-r_d-X^{-1}{r}_{vx}+T^{-1}{r}_{wt}+T^{-1}W{r}_{ub}.\end{array}$$

D 和 R 定义中的所有矩阵的逆都很容易计算，因为矩阵是对角矩阵。

要从公式 4 推导公式 5，请注意 公式 5 的第二行与 公式 4 的第二个矩阵行相同。公式 5 的第一行来自于先求解公式 4 的最后两行以得到 Δv 和 Δw，再求解以得到 Δt。

为了求解公式 5，该算法采用 Altman 和 Gondzio·的论著 [1] 中的基本思路。该算法通过 LDL 分解求解对称方程组。正如 Vanderbei 和 Carpenter 等在论著·[2] 中指出的那样，这种分解无需进行主元消去也可保持数值稳定，因此执行速度可以很快。

在计算校正的牛顿步后，算法执行更多计算，以获得更长的当前步，并准备好计算更优的后续步。这些多重校正计算可以提高性能和稳健性。关详细信息，请参阅 Gondzio 的论著 [4]。

quadprog 满预测-校正算法与 linprog 'interior-point' 算法基本相同，但前者还包括二次项。请参阅预测-校正。

参考

[1] Altman, Anna and J. Gondzio. Regularized symmetric indefinite systems in interior point methods for linear and quadratic optimization. Optimization Methods and Software, 1999. Available for download here.

[2] Vanderbei, R. J. and T. J. Carpenter. Symmetric indefinite systems for interior point methods. Mathematical Programming 58, 1993. pp. 1–32. Available for download here.

停止条件

预测-校正算法会持续迭代，直至达到一个可行（在容差范围内满足约束）且相对步长较小的点。具体来说，定义

$$\rho =\max \left(1,\Vert H\Vert ,\Vert \overline{A}\Vert ,\Vert A_{eq}\Vert ,\Vert c\Vert ,\Vert \overline{b}\Vert ,\Vert b_{eq}\Vert \right)$$

当满足所有下列条件时，算法停止：

$$\begin{array}{c}{\Vert r_p\Vert }_1+{\Vert r_{ub}\Vert }_1\le \rho \text{TolCon}\\ {\Vert r_d\Vert }_{\infty }\le \rho \text{TolFun}\\ r_c\le \text{TolFun,}\end{array}$$

其中

$$r_c=\underset{i}{\max }\left(\min \left(\left|x_i{v}_i\right|,\left|x_i\right|,\left|v_i\right|\right),\min \left(\left|t_i{w}_i\right|,\left|t_i\right|,\left|w_i\right|\right)\right).$$

r_c 本质上是在测量互补残差 xv 和 tw 的大小，它们均为由零组成的向量，位于同一个解处。

不可行性检测

quadprog 在每次迭代中计算一个评价函数 φ。评价函数是可行性的测度。如果评价函数变得过大，quadprog 将停止。在本例中，quadprog 声明该问题不可行。

评价函数与问题的 KKT 条件相关 - 请参阅预测-校正。使用以下定义：

$$\begin{array}{c}\rho =\max \left(1,\Vert H\Vert ,\Vert \overline{A}\Vert ,\Vert A_{eq}\Vert ,\Vert c\Vert ,\Vert \overline{b}\Vert ,\Vert b_{eq}\Vert \right)\\ r_{\text{eq}}=A_{\text{eq}}x-b_{\text{eq}}\\ r_{\text{ineq}}=\overline{A}x-\overline{b}-s\\ r_{\text{d}}=Hx+c+A_{\text{eq}}^T{\lambda }_{\text{eq}}+{\overline{A}}^T{\overline{\lambda }}_{\text{ineq}}\\ g=\frac{1}{2}{x}^{T}Hx+c^{T}x-{\overline{b}}^T{\overline{\lambda }}_{\text{ineq}}-b_{\text{eq}}^T{\lambda }_{\text{eq}}.\end{array}$$

$$\overline{A}$$ 和 $$\overline{b}$$ 表示扩展的线性不等式系数，其中包含表示边界的项，可用于稀疏算法。表示法 $${\overline{\lambda }}_{\text{ineq}}$$ 与之相似，表示线性不等式约束的拉格朗日乘数，包括边界约束。这在预测-校正中被称为 z，$${\lambda }_{\text{eq}}$$ 被称为 y。

φ 评价函数是

$$\frac{1}{\rho }\left(\max \left({\Vert r_{\text{eq}}\Vert }_{\infty },{\Vert r_{\text{ineq}}\Vert }_{\infty },{\Vert r_{\text{d}}\Vert }_{\infty }\right)+g\right).$$

如果此评价函数变得太大，则 quadprog 会声明该问题不可行并停止，同时返回退出标志 -2。

trust-region-reflective quadprog 算法

Optimization Toolbox™ 求解器中使用的许多方法都基于信赖域，这是一个简单而功能强大的优化概念。

要理解信赖域优化方法，请考虑无约束最小化问题，最小化 f(x)，该函数接受向量参量并返回标量。假设您现在位于 n 维空间中的点 x 处，并且您要寻求改进，即移至函数值较低的点。基本思路是用较简单的函数 q 来逼近 f，该函数需能充分反映函数 f 在点 x 的邻域 N 中的行为。此邻域是信赖域。试探步 s 是通过在 N 上进行最小化（或近似最小化）来计算的。以下是信赖域子问题，

$$\underset{s}{\min }\left\{q(s),s\in N\right\}.$$

(6)

如果 f(x + s) < f(x)，当前点更新为 x + s；否则，当前点保持不变，信赖域 N 缩小，算法再次计算试探步。

在定义特定信赖域方法以最小化 f(x) 的过程中，关键问题是如何选择和计算逼近 q（在当前点 x 上定义）、如何选择和修改信赖域 N，以及如何准确求解信赖域子问题。本节重点讨论无约束问题。后面的章节讨论由于变量约束的存在而带来的额外复杂性。

在标准信赖域方法 ([48]) 中，二次逼近 q 由 F 在 x 处的泰勒逼近的前两项定义；邻域 N 通常是球形或椭圆形。以数学语言表述，信赖域子问题通常写作

$$\min \left\{\frac{1}{2}{s}^{T}Hs+s^{T}g\text{ such that }\Vert Ds\Vert \le \Delta \right\},$$

(7)

其中，g 是 f 在当前点 x 处的梯度，H 是黑塞矩阵（二阶导数的对称矩阵），D 是对角缩放矩阵，Δ 是正标量，‖ . ‖ 是 2-范数。存在求解公式 7 的好算法（请参阅[48]）；此类算法通常涉及计算 H 的所有特征值，并将牛顿法应用于以下久期方程

$$\frac{1}{\Delta }-\frac{1}{\Vert s\Vert }=0.$$

此类算法提供公式 7 的精确解。但是，它们要耗费与 H 的几个分解成比例的时间。因此，对于大型问题，需要一种不同方法。文献（[42] 和 [50]）中提出了几种基于公式 7 的逼近和启发式方法建议。Optimization Toolbox 求解器采用的逼近方法是将信赖域子问题限制在二维子空间 S 内（[39] 和 [42]）。一旦计算出子空间 S，即使需要完整的特征值/特征向量信息，求解公式 7 的工作量也不大（因为在子空间中，问题只是二维的）。现在的主要工作已转移到子空间的确定上。

二维子空间 S 是借助下述预条件共轭梯度法确定的。求解器将 S 定义为由 s₁ 和 s₂ 确定的线性空间，其中 s₁ 是梯度 g 的方向，s₂ 是近似牛顿方向，即下式的解

$$H\cdot s_2=-g,$$

(8)

或是负曲率的方向，

$$s_2^T\cdot H\cdot s_2<0.$$

(9)

以此种方式选择 S 背后的理念是强制全局收敛（通过最陡下降方向或负曲率方向）并实现快速局部收敛（通过牛顿步，如果它存在）。

现在，我们可以很容易地给出基于信赖域的无约束最小化的大致框架：

1. 构造二维信赖域子问题。

2. 求解公式 7 以确定试探步 s。

3. 如果 f(x + s) < f(x)，则 x = x + s。

4. 调整 Δ。

重复这四个步骤，直到收敛。信赖域维度 Δ 根据标准规则进行调整。具体来说，它会在试探步不被接受（即 f(x + s) ≥ f(x)）时减小。有关这方面的讨论，请参阅[46] 和 [49]。

Optimization Toolbox 求解器用专用函数处理 f 的一些重要特例：非线性最小二乘、二次函数和线性最小二乘。然而，其底层算法思路与一般情况相同。这些特例将在后面的章节中讨论。

此处使用子空间信赖域方法来确定搜索方向。然而，它不是像在非线性求最小值情形中那样可能将步长限制为一个反射步，而是在每次迭代中进行分段反射线搜索。请参阅[45] 了解线搜索的详细信息。

预条件共轭梯度法

求解大型对称正定线性方程组 Hp = –g 的一种常用方法是预条件共轭梯度法 (PCG)。这种迭代方法需要能够计算 H·v 形式的矩阵-向量积，其中 v 是任意向量。对称正定矩阵 M 是 H 的预条件子。也就是说，M = C²，其中 C^–1HC^–1 是良态矩阵或具有聚类特征值的矩阵。

在最小化上下文中，您可以假设黑塞矩阵 H 是对称矩阵。然而，只有在强最小解的邻域内才能保证 H 是正定矩阵。算法 PCG 在遇到负（或零）曲率方向（即 d^THd ≤ 0）时退出。PCG 输出方向 p 要么是负曲率方向，要么是牛顿方程组 Hp = –g 的近似解。无论哪种情况，p 都有助于定义在非线性最小化信赖域方法中讨论的信赖域方法中使用的二维子空间。

线性等式约束

线性约束使无约束最小化中所述的情形复杂化了。不过，我们仍可以简洁高效地沿用前述的基本思路。Optimization Toolbox 求解器中的信赖域方法生成严格可行的迭代。

一般的线性等式约束最小化问题可以写作

$$\min \left\{f(x)\text{ such that }Ax=b\right\},$$

(10)

其中 A 是 m×n 矩阵 (m ≤ n)。一些 Optimization Toolbox 求解器对 A 进行预处理，以使用基于 A^T 的 LU 分解的方法来消除严格的线性相关性 [46]。此处假定 A 的秩为 m。

用于求解公式 10 的方法在两个重要方面不同于无约束方法。首先，它使用稀疏最小二乘步计算初始可行点 x₀，满足 Ax₀ = b。然后，它用简化的预条件共轭梯度 (RPCG) 代替算法 PCG（请参阅[46]），以便计算近似的简化牛顿步（或 A 的零空间中的负曲率方向）。关键的线性代数步涉及求解以下形式的方程组

$$\left[\begin{array}{cc}C& {\tilde{A}}^T\\ \tilde{A}& 0\end{array}\right]\left[\begin{array}{c}s\\ t\end{array}\right]=\left[\begin{array}{c}r\\ 0\end{array}\right],$$

(11)

其中，$$\tilde{A}$$ 逼近 A（在不丢失秩的前提下，A 的小非零值会设置为零）且 C 是 H 的稀疏对称正定逼近，即 C = H。有关详细信息，请参阅[46]。

框约束

框约束问题具有以下形式

$$\min \left\{f(x)\text{ such that }l\le x\le u\right\},$$

(12)

其中，l 是由下界组成的向量，u 是由上界组成的向量。l 的部分（或全部）分量可以等于 –∞，u 的部分（或全部）分量可以等于 ∞。该方法生成一系列严格可行点。此处用到了两种技巧，在保持可行性的同时实现稳健的收敛行为。其一，用缩放的修正牛顿步取代无约束牛顿步（以定义二维子空间 S）。其次，反射用于增大步长。

缩放的修正牛顿步源于对公式 12 的库恩-塔克必要条件的检验，

$${\left(D(x)\right)}^{-2}g=0,$$

(13)

其中

$$D(x)=\text{diag}\left({\left|v_k\right|}^{-1/2}\right),$$

对于每个 1 ≤ i ≤ n，向量 v(x) 定义如下：

• 如果 g_i < 0 且 u_i < ∞，则 v_i = x_i – u_i

• 如果 g_i ≥ 0 且 l_i > –∞，则 v_i = x_i – l_i

• 如果 g_i < 0 且 u_i = ∞，则 v_i = –1

• 如果 g_i ≥ 0 且 l_i = –∞，则 v_i = 1

非线性系统公式 13 不是处处可微的。不可微分性发生在 v_i = 0 时。您可以通过保持严格的可行性（例如限制 l < x < u）来避免这些不可微分的点。

公式 13 给出的非线性方程组的缩放的修正牛顿步 s_k 定义为以下线性方程组的解

$$\widehat{M}D{s}^N=-\widehat{g}$$

(14)

在第 k 次迭代中的解，其中

$$\widehat{g}=D^{-1}g=\text{diag}\left({\left|v\right|}^{1/2}\right)g,$$

(15)

并且

$$\widehat{M}=D^{-1}H{D}^{-1}+\text{diag}(g)J^v.$$

(16)

在此处，J^v 表示 |v| 的雅可比。对角矩阵 J^v 的每个对角分量等于 0、–1 或 1。如果 l 和 u 的所有分量均为有限的，则 J^v = diag(sign(g))。在满足 g_i = 0 的点上，v_i 可能不可微分。$$J_{ii}^v=0$$ 在此类点上定义。这种类型的不可微分性不值得关注，因为对于这种分量，v_i 采用哪个值并不重要。此外，|v_i| 在此点上仍是不连续的，但函数 |v_i|·g_i 是连续的。

其次，反射用于增大步长。单个反射步定义如下。给定与边界约束相交的步 p，假设与 p 相交的第一个边界约束是第 i 个边界约束（可以是第 i 个上界，也可以是第 i 个下界），则反射步 p^R = p，但在第 i 个分量中除外，其中 p^R_i = –p_i。

active-set quadprog 算法

完成预求解步骤后，active-set 算法分两个阶段进行。

• 阶段 1 - 获得一个在所有约束下可行的点。

• 阶段 2 - 以迭代方式降低目标函数，每次迭代均保持一组活动约束并保持可行性。

active-set 策略（也称为投影方法）类似于吉尔等人在 [18] 和 [17] 中所述的策略。

预求解步骤

该算法首先尝试通过消除冗余和简化约束来简化问题。在预求解步骤中执行的任务可以包括以下内容：

• 检查是否有任何变量具有相等的上界和下界。如果有，请检查可行性，然后固定变量的值以去除变量。

• 检查是否有任何线性不等式约束只涉及一个变量。如果是，请检查可行性，然后将线性约束更改为一个边界。

• 检查是否有任何线性等式约束只涉及一个变量。如果是，请检查可行性，然后固定该变量的值以去除该变量。

• 检查是否有任何线性约束矩阵具有全零行。如果是，请检查可行性，然后删除这些行。

• 确定边界和线性约束是否一致。

• 检查是否有任何变量仅作为线性项出现在目标函数中，而不出现在任何线性约束中。如果是这样，请检查可行性和有界性，然后将变量固定为其适当的边界值。

• 通过添加松弛变量，将任何线性不等式约束更改为线性等式约束。

如果算法检测到不可行或无界的问题，它会停止并发出相应的退出消息。

该算法可能到达一个表示解的单个可行点。

如果算法在预求解步骤中没有检测到不可行或无界的问题，并且预求解没有产生解，则算法会继续其后续步骤。到达停止条件后，算法会重新构造原始问题，这将撤消任何预求解变换。这最后一步是后求解步骤。

有关详细信息，请参阅 Gould 和 Toint 合著的 [63]。

阶段 1 算法

在阶段 1 中，算法尝试找到一个点 x，它满足所有约束而不考虑目标函数。仅当提供的初始点 x0 不可行时，quadprog 才运行阶段 1 算法。

首先，算法尝试找到一个在所有等式约束下都可行的点，例如 X = Aeq\beq。如果方程 Aeq*x = beq 没有解 x，则算法停止。如果有解 X，下一步是满足边界和线性不等式。在没有等式约束的情况下，设置 X = x0，即初始点。

从 X 开始，算法计算 M = max(A*X – b, X - ub, lb – X)。如果 M <= options.ConstraintTolerance，则点 X 是可行的，并且阶段 1 算法停止。

如果 M > options.ConstraintTolerance，算法为辅助线性规划问题引入非负松弛变量 γ

$$\underset{x,\gamma }{\min }\gamma$$

使得

$$\begin{array}{c}Ax-\gamma \le b\\ \text{Aeq}x=\text{beq}\\ \text{lb}-\gamma \le x\le \text{ub}+\gamma \\ \gamma \ge -\rho .\end{array}$$

此处，ρ 是 ConstraintTolerance 选项乘以 A 和 Aeq 中最大元素的绝对值。如果算法找到 γ = 0，并且一个点 x 满足方程和不等式，则 x 是一个可行的阶段 1 点。如果在 γ = 0 的情况下，辅助线性规划问题没有解 x，则阶段 1 问题不可行。

为了求解辅助线性规划问题，算法设置 γ₀ = M + 1，设置 x₀ = X，并将活动集初始化为固定变量（如果有）和所有等式约束。算法将线性规划变量 p 重新表示为 x 相对于当前点 x₀ 的偏移量，即 x = x₀ + p。该算法使用阶段 2 中求解二次规划问题相同的迭代来求解线性规划问题，并使用适当的修正黑塞矩阵。

阶段 2 算法

使用变量 d，问题可表示为

$$\begin{array}{c}\underset{d\in {\Re }^n}{\min }q(d)=\frac{1}{2}{d}^{T}Hd+c^{T}d,\\ A_{i}d=b_i,i=1,\mathrm{...},m_e\\ A_{i}d\le b_i,i=m_e+1,\mathrm{...},m.\end{array}$$

(17)

此处，A_i 指 m×n 矩阵 A 的第 i 行。

在阶段 2 中，活动集 $${\overline{A}}_k$$，它是在解点处的活动约束（约束边界上的约束）的估计值。

算法在每一迭代 k 中更新 $${\overline{A}}_k$$，以构成搜索方向 d_k 的基础。等式约束始终维持在活动集 $${\overline{A}}_k$$ 中。计算搜索方向 d_k 以最小化目标函数，同时确保该方向位于活动约束边界上。算法基于 Z_k 形成 d_k 的可行子空间，Z_k 的列与活动集 $${\overline{A}}_k$$ 的估计值正交（即 $${\overline{A}}_k{Z}_k=0$$）。这就确保了搜索方向（由 Z_k 的列的任意组合经线性求和形成）始终位于活动约束的边界上。

算法基于矩阵 $${\overline{A}}_k^T$$ 的 QR 分解的最后 n – l 个列形成矩阵 Z_k，其中 l 是活动约束的数目且 l < n。也就是说，Z_k 由下式给出

$$Z_k=Q\left[:,l+1:n\right],$$

(18)

其中

$$Q^T{\overline{A}}_k^T=\left[\begin{array}{c}R\\ 0\end{array}\right].$$

在找到 Z_k 后，算法寻找新搜索方向 d_k 以最小化 q(d)，其中 d_k 在活动约束的零空间中。也就是说，d_k 是 Z_k 的列的线性组合：对于某个向量 p，$${\widehat{d}}_k=Z_{k}p$$。

通过代换 d_k，将二次函数视为 p 的函数，可得

$$q(p)=\frac{1}{2}{p}^T{Z}_k^{T}H{Z}_{k}p+c^T{Z}_{k}p.$$

(19)

求此方程关于 p 的微分，可得

$$\nabla q(p)=Z_k^{T}H{Z}_{k}p+Z_k^{T}c.$$

(20)

∇q(p) 称为二次函数的投影梯度，因为它是投影到 Z_k 定义的子空间中的梯度。项 $$Z_k^{T}H{Z}_k$$ 称为投影黑塞矩阵。假设投影黑塞矩阵 $$Z_k^{T}H{Z}_k$$ 是半正定矩阵，在由 Z_k 定义的子空间中函数 q(p) 的最小值出现在 ∇q(p) = 0 时，它是以下线性方程组的解

$$Z_k^{T}H{Z}_{k}p=-Z_k^{T}c.$$

(21)

然后，算法执行以下求解步

$$x_{k+1}=x_k+\alpha d_k,$$

其中

$$d_k=Z_{k}p.$$

由于目标函数的二次性质，在每次迭代中，步长 α 只有两种选择。沿 d_k 的单位步是 $${\overline{A}}_k$$ 的零空间内到函数最小值的精确步。如果算法可以采取这样的步而不违反约束，则此步就是二次规划的解（公式 18）。否则，沿 d_k 到最近约束的步小于单位步，并且下一次迭代时，算法将在活动集中加入一个新的约束。沿任一方向 d_k 到约束边界的距离由下式给出

$$\alpha =\underset{i\in \left\{1,\mathrm{...},m\right\}}{\min }\left\{\frac{-\left(A_i{x}_k-b_i\right)}{A_i{d}_k}\right\},$$

该定义针对不在活动集中的约束，其中的方向 d_k 朝向约束边界，即 $$A_i{d}_k>0,i=1,\mathrm{...},m$$。

如果活动集包含 n 个独立约束，且最小值的位置未知，则算法会计算满足以下非奇异线性方程组的拉格朗日乘数 λ_k

$${\overline{A}}_k^T{\lambda }_k=c+H{x}_k.$$

(22)

如果 λ_k 的所有元素均为非负，则 x_k 是二次规划问题 公式 1 的最优解。但是，如果 λ_k 的任一分量为负，并且该分量不对应于等式约束，则最小化未完成。算法会删除对应于最小的负乘数的元素，并开始一次新迭代。

有时，即使求解器用尽了所有非负拉格朗日乘数，一阶最优性测度仍高于容差，或仍不满足约束容差。在这些情况下，求解器会尝试按照 [1] 中所述的重启过程来获得更好的解。在此过程中，求解器会放弃当前活动约束集，然后稍微放松约束，并构造新的活动约束集，同时尝试以避免循环（重复返回到相同状态）的方式求解问题。如有必要，求解器可以多次执行重启过程。

注意

不要尝试通过为 ConstraintTolerance 和 OptimalityTolerance 选项设置较大的值来提前停止算法。通常，求解器在进行迭代时不检查这些值，直到到达一个潜在的停止点时才检查是否满足容差。

有时，active-set 算法可能难以检测问题何时是无界的。如果无界性 v 的方向是二次项 v'Hv = 0 的方向，则可能出现此问题。为了数值稳定性和实现乔列斯基分解，active-set 算法对二次目标加上一个值很小的严格凸项。这个值很小的项使得目标函数的边界不会为 –Inf。在这种情况下，active-set 算法会达到一个迭代限制，而不是报告解是无界的。换句话说，该算法停止时会返回退出标志 0（而不是 –3）。

参考

[1] Gill, P. E., W. Murray, M. A. Saunders, and M. H. Wright. A practical anti-cycling procedure for linearly constrained optimization. Math. Programming 45 (1), August 1989, pp. 437–474.

热启动

当您以热启动对象作为起点运行 quadprog 或 lsqlin 'active-set' 算法时，求解器会尝试跳过第 1 阶段和第 2 阶段的许多步骤。热启动对象包含处于活动状态的约束集，此活动约束集对于新问题可能是正确的或接近正确。因此，求解器可以避免运行一些迭代而直接向活动约束集添加约束。此外，初始点可能接近新问题的解。有关详细信息，请参阅 optimwarmstart。