37 约束优化方法

考虑约束最优化问题。 约束最优化问题比无约束优化问题复杂得多, 也没有很好的通用解决方法, 在这个方面有大量研究。 现有的方法大致上有如下三类:

  • 把约束问题转化为一系列无约束问题, 用这些无约束问题的极小值点去逼近约束极小值点, 称这样的方法为序列无约束优化方法(SUMT), 如惩罚函数法、乘子罚函数法。

  • 在迭代的每一步用一个二次函数逼近目标函数, 用线性约束近似一般约束,构造一系列二次规划来逼近原问题, 称这样的方法为序列二次规划(SQP)法。

  • 每次迭代找到可行下降方向, 保证迭代始终处于可行域内, 这样的方法称为可行方向法。

37.1 约束的化简

最简单的一些约束可以通过参数变换转化成无约束问题。 例如,对\(f(x)\)\(x \in (0, \infty)\)求最小值, 可令\(x = e^t\), \(h(u) = f(e^u), u \in (-\infty, \infty)\), 若求得\(\mathop{\text{argmin}}_{u \in \mathbb R} h(u) = u^*\), 则\(x^* \stackrel{\triangle}{=} e^{u^*}\)\(f(x)\)\((0, \infty)\)的最小值点。 类似地,对\(x\)限制在有限开区间\((a,b)\)的情形,可以作变换 \[\begin{align*} x = a + \frac{b-a}{1 + e^{-u}}, \ u \in (-\infty, \infty), \end{align*}\]\(x\)限制在有限闭区间\([a,b]\)的情形, 可以做变换 \[\begin{align*} x = a + (b-a) \sin^2 u, \ u \in (-\infty, \infty), \end{align*}\]\(\boldsymbol x = (x_1, x_2, x_3)\)限制为\(0 \leq x_1 \leq x_2 \leq x_3\)的情形, 可以做变换 \[\begin{align*} x_1 = u_1^2, \quad x_2 = u_1^2 + u_2^2, \quad x_3 = u_1^2 + u_2^2 + u_3^2, \ (u_1, u_2, u_3) \in \mathbb R^3, \end{align*}\] 对于\(\boldsymbol x = (x_1, x_2)\)限制为\(x_1 > 0, x_2 > 0, x_1 + x_2 < 1\)的情形, 可以做变换 \[\begin{align*} x_1 =& \frac{e^{t_1}}{1 + e^{t_1} + e^{t_2}} , \quad x_2 = \frac{e^{t_2}}{1 + e^{t_1} + e^{t_2}} , \ (t_1, t_2) \in \mathbb R^2. \end{align*}\] 这些变换可以把原来的约束优化问题转化为无约束优化问题, 或者减少约束个数。

对于等式约束, 如果有显式解, 可以从中解出部分自变量, 简化原来的问题, 比如下面的情况。

37.2 仅含线性等式约束的情形

如果约束仅包含线性方程组, 可以把问题(34.1)用矩阵形式写成 \[\begin{align} \begin{cases} & \mathop{\text{argmin}}\limits_{\boldsymbol x \in \mathbb R^d } f(\boldsymbol x), \ \text{s.t.} \\ & A^T \boldsymbol x = \boldsymbol b, \end{cases} \tag{37.1} \end{align}\] 其中\(A\)\(d \times p\)矩阵, 各列为\(\boldsymbol a_i, i=1,\dots,p\)\(\boldsymbol b = (b_1, \dots, b_p)^T\), 即问题仅有\(p\)个线性约束\(\boldsymbol a_i^T \boldsymbol x = b_i, i=1,2,\dots,p\)。 设\(p < d\)\(A\)列满秩, 当\(A\)规模很小时, 可以预先用消元法把\(\boldsymbol x\)\(p\)个分量用其余的\(d-p\)个线性表出, 把问题转化为\(d-p\)维的无约束优化问题。

下面讨论当\(A\)规模较大时的处理方法。

\(\mu(A)\)\(A\)的各列在\(\mathbb R^d\)中张成的线性空间, 即 \[\begin{align} \mu(A) \stackrel{\triangle}{=} \{ \boldsymbol x \in \mathbb R^d: \boldsymbol x = A \boldsymbol u, \boldsymbol u \in \mathbb R^p \} = \{ \boldsymbol x \in \mathbb R^d: \boldsymbol x = \sum_{i=1}^p c_i \boldsymbol a_i, \ c_1, \dots, c_p \in \mathbb R \} \end{align}\]\(\mathbb R^d\)可以分解为\(\mu(A)\)\(\mu(A)\)的正交补空间 \[\begin{align} \mu(A)^\perp \stackrel{\triangle}{=} \{ \boldsymbol x \in \mathbb R^d: A^T \boldsymbol x = 0 \} \end{align}\] 的直和: \[\begin{align*} \mathbb R^d = \mu(A) \oplus \mu(A)^\perp. \end{align*}\]\(A\)\(QR\)分解 \[\begin{align*} A = Q \left( \begin{array}{c} R \\ 0 \end{array} \right) = (Q_1, Q_2) \left( \begin{array}{c} R \\ 0 \end{array} \right) = Q_1 R, \end{align*}\] 其中\(Q\)\(d\times d\)正交阵, \(R\)\(p\times p\)上三角阵, \(Q_1\)\(d \times p\)矩阵, 易见\(Q_1\)各列构成\(\mu(A)\)的标准正交基, \(Q_2\)各列构成\(\mu(A)^\perp\)的标准正交基。 若\(\boldsymbol x\)满足约束\(A \boldsymbol x = b\), 设 \[\begin{align*} \boldsymbol x = Q_1 \boldsymbol u \oplus Q_2 \boldsymbol v, \ \boldsymbol u \in \mathbb R^p, \ \boldsymbol v \in \mathbb R^{d-p}, \end{align*}\]\[\begin{align*} A^T \boldsymbol x = R^T Q_1^T \boldsymbol x = R^T \boldsymbol u = \boldsymbol b \end{align*}\] 可得\(\boldsymbol u = R^{-T} \boldsymbol b\)(简记\((R^T)^{-1}\)\(R^{-T}\)), 于是 \[\begin{align} \boldsymbol x = Q_1 R^{-T} \boldsymbol b + Q_2 \boldsymbol v, \ \boldsymbol v \in \mathbb R^{d-p}, \tag{37.2} \end{align}\] 只要对目标函数\(f(Q_1 R^{-T} \boldsymbol b + Q_2 \boldsymbol v)\), \(\boldsymbol v \in \mathbb R^{d-p}\)求解无约束最优化问题即可。

\(\boldsymbol x\)的分解式(37.2)中, 易见 \[\begin{align*} P \stackrel{\triangle}{=}& Q_2 Q_2^T = I_n - Q_1 Q_1^T = I_n - A R^{-1} R^{-T} A^T = I_n - A (R^T R)^{-1} A^T \\ =& I_n - A (A^T A)^{-1} A^T, \end{align*}\] 这是向正交补空间\(\mu(A)^\perp\)投影的投影矩阵, 所以(37.2)中属于正交补空间\(\mu(A)^\perp\)的部分\(Q_2 \boldsymbol v\), 在已知\(\boldsymbol x\)时可以用投影表示为\(P \boldsymbol x\)。 后文中的投影梯度法利用了这样的思想。

如果\(\boldsymbol x^{(t)}\)是一个可行点, 为了使得\(\boldsymbol x^{(t+1)} = \boldsymbol x^{(t)} + \boldsymbol d\)也是可行点, 需要\(A^T \boldsymbol d = 0\)。 注意问题(37.1)是一般约束问题(34.1)仅含等式约束的版本, 这里\(c_i(\boldsymbol x) = \boldsymbol a_i^T \boldsymbol x\), \(\boldsymbol a_i\)\(A\)的第\(i\)列, 可见\(\nabla c_i(\boldsymbol x) = \boldsymbol a_i\), 上述对\(\boldsymbol d\)的要求\(A^T \boldsymbol d=0\)可以表示为 \[\begin{align*} \boldsymbol d^T \nabla c_i(\boldsymbol x) = 0, i=1,\dots, p, \end{align*}\] 这个要求可以推广到非线性等式约束的情形, 即为了\(\boldsymbol x^{(t)} + \boldsymbol d\)继续满足等式约束, 搜索方向\(\boldsymbol d\)必须和等式约束函数的梯度都正交, 否则会使得\(c_i(\boldsymbol x)\)的值改变。

37.3 线性约束最优化方法

在约束最优化问题中, 如果所有约束都是线性函数, 这样的问题相对比较容易处理, 但是也能体现约束优化问题的困难。 如下的约束优化问题称为线性约束优化问题: \[\begin{align} \begin{cases} & \mathop{\text{argmin}}\limits_{\boldsymbol x \in \mathbb R^d } f(\boldsymbol x), \ \text{s.t.} \\ & \boldsymbol a_i^T \boldsymbol x = b_i, \ i=1,\dots, p, \\ & \boldsymbol a_i^T \boldsymbol x \leq b_i, \ i=p+1, \dots, q, \end{cases} \tag{37.3} \end{align}\] 其中\(\boldsymbol a_i \in \mathbb R^d, i=1, \dots, p+q\), \(b_i \in \mathbb R, i=1,\dots,p+q\)

前面已经讨论了仅有线性等式约束的情形, 下面讨论含有线性不等式约束的情形。

37.3.1 投影梯度法

(37.2)看出, 当约束为线性等式约束时, 搜索时应该在和各\(\boldsymbol a_i\)正交的方向上搜索。 这提示我们, 对于包含不等式约束的优化问题(37.3), 如果\(\boldsymbol x^{(t)}\)是一个可行点, 希望找到使得函数值下降的可行方向, 只需考虑所有等式约束和起作用的不等式约束。 假设\(\boldsymbol x^{(t)}\)处起作用的不等式约束下标为\(p+1, \dots, p+r\), 记 \[\begin{align} \mathcal B = \{ \boldsymbol d \in \mathbb R^d:\ \boldsymbol a_i^T \boldsymbol d = 0, \ i=1, \dots, p+r \} \tag{37.4} \end{align}\] 设负梯度\(-\nabla f(\boldsymbol x^{(t)})\) 投影到\(\mathcal B\)中得到\(\boldsymbol d\), 如果\(\boldsymbol d \neq \boldsymbol 0\), \(\boldsymbol d\)就是一个可行下降方向, 沿方向\(\boldsymbol d\)搜索使得函数值下降且保持所有\(p+q\)个约束成立。 有了可行下降方向\(\boldsymbol d\)以后, 从\(\boldsymbol x^{(t)}\)出发沿\(\boldsymbol d\)方向进行一维搜索, 得到下一个近似极小值点\(\boldsymbol x^{(t+1)}\), 如此迭代逼近。

如果投影得到的\(\boldsymbol d = \boldsymbol 0\), 必存在\(\lambda_i^{(t)}, i=1,\dots,p+r\), 满足 \[\begin{align*} \sum_{i=1}^{p+r} \lambda_i^{(t)} \boldsymbol a_i = -\nabla f(\boldsymbol x^{(t)}) . \end{align*}\] 解出\(\lambda_i^{(t)}, i=1,\dots,p+r\), 如果\(\lambda_i^{(t)} \geq 0, i=p+1, \dots, p+r\), 只要令\(\lambda_i^{(t)} = 0, i=p+r+1, \dots, p+q\), \(\boldsymbol\lambda^{(t)} = (\lambda_1^{(t)}, \dots, \lambda_{p+q}^{(t)})^T\), 则\((\boldsymbol x^{(t)}, \boldsymbol\lambda^{(t)})\)是问题(37.3)的KKT对, 算法不再继续,只要设法判断\(\boldsymbol x^{(t)}\)是否极小值点。

如果\(\boldsymbol d = \boldsymbol 0\)但是解出的\(\lambda_i^{(t)}, i=p+1,\dots,p+r\)中有负值, 则设\(\lambda_{i_0}^{(t)} = \min\{\lambda_i^{(t)}, i=p+1,\dots,p+r \}\), 令 \[\begin{align*} \mathcal B' = \{ \boldsymbol d \in \mathbb R^d:\ \boldsymbol a_i^T \boldsymbol d = 0, \ i=1, \dots, p+r, i \neq i_0 \} \end{align*}\]\(\tilde{\boldsymbol d}\)\(-\nabla f(\boldsymbol x^{(t)})\)\(\mathcal B'\)的投影, 则\(\tilde{\boldsymbol d}\)一定是可行下降方向。

这种方法称为投影梯度法, 是一种比较早期的可行方向法, 具体算法略。 投影梯度法以及与之类似的简约梯度法优点是比较简单, 缺点是仅利用梯度信息而没有试图采用高阶逼近, 收敛速度慢。

37.3.2 简约梯度法

在一般线性约束优化问题(37.3)中, 一般无约束的分量\(x_i\)可以写成 \(x_j = x_j^+ - x_j^-\), 其中\(x_j^+ = \max(x_j, 0)\), \(x_j^- = -\min(x_j, 0)\), 这样,可以设所有的分量都满足\(x_j \geq 0\)。 对不等式约束\(\boldsymbol a_i^T \boldsymbol x \geq b_i\), 可以引入新自变量\(z_i \geq 0\)使得\(\boldsymbol a_i^T \boldsymbol x - z_i = b_i\), 这样,可以设所有自变量非负, 所有的其它约束都是等式约束, 于是问题(37.3)可以写成 \[\begin{align} \begin{cases} & \mathop{\text{argmin}}\limits_{\boldsymbol x \in \mathbb R^d } f(\boldsymbol x), \ \text{s.t.} \\ & C \boldsymbol x = \boldsymbol b, \\ & \boldsymbol x \geq \boldsymbol 0 . \end{cases} \tag{37.5} \end{align}\]\(C\)\(p \times d\)矩阵, 并设\(C\)行满秩, 不妨设\(C = (C_B, C_N)\), 其中\(C_B\)\(p \times p\)满秩方阵, 把\(\boldsymbol x\)拆分为\((\boldsymbol x_B, \boldsymbol x_N)\), 由\(C \boldsymbol x = \boldsymbol b\)可解出 \[\begin{align*} \boldsymbol x_B = \boldsymbol x_B(\boldsymbol x_N) = C_B^{-1} \boldsymbol b - C_B^{-1} C_N \boldsymbol x_N, \ \boldsymbol x_N \in \mathbb R^{d-p}, \end{align*}\] 问题可简化为 \[\begin{align} \begin{cases} & \mathop{\text{argmin}}\limits_{\boldsymbol x_N \in \mathbb R^{d-p} } h(\boldsymbol x_N) = f(\boldsymbol x_B(\boldsymbol x_N), \boldsymbol x_N), \ \text{s.t.} \\ & \boldsymbol x_N \geq \boldsymbol 0 , \boldsymbol x_B(\boldsymbol x_N) \geq \boldsymbol 0. \end{cases} \end{align}\]\(f(\boldsymbol x)\)的梯度前\(p\)个为\(G_B(\boldsymbol x)\), 后\(d-p\)个为\(G_N(\boldsymbol x)\), 则目标函数\(h(\boldsymbol x_N)\)的梯度为 \[\begin{align*} \nabla h(\boldsymbol x_N) = G_N(\boldsymbol x) - (C_B^{-1} C_N)^T G_B(\boldsymbol x), \end{align*}\] 称为\(f(\boldsymbol x)\)的简约梯度, 其中\(\boldsymbol x = (\boldsymbol x_B(\boldsymbol x_N), \boldsymbol x_N)\)。 记\(\tilde {\boldsymbol d}_N = -\nabla h(\boldsymbol x_N)\), 这是目标函数\(h(\boldsymbol x_N)\)的下降方向。 在迭代逼近约束最小值点的过程中, 如果负简约梯度方向\(\tilde {\boldsymbol d}_N\)所有分量都大于零或等于零, 则它指向\(\boldsymbol x_N\)的分量都增加或不变的方向, 关于\(\boldsymbol x_N\)部分是可行下降方向, \(\boldsymbol x_N\)部分只要延\(\tilde {\boldsymbol d}\)方向搜索即可; 如果\(\tilde {\boldsymbol d}_N\)的某个分量(设为第\(j\)个)小于零, 因为则该分量指向\(x_j\)减小的方向, 有可能使得搜索跳出可行域, 所以要把\(\tilde {\boldsymbol d}_N\)的这样的分量修改为\(x_j \tilde d_j\), 这样,当\(x_j>0\)时, 搜索方向是使得\(x_j\)减少的方向, 当\(x_j = 0\)时,搜索方向使得\(x_j\)不变从而不会离开可行域。 设修改后的\(\tilde{\boldsymbol d}_N\)\(\boldsymbol d_N\), 对于\(\boldsymbol x_B\)部分, 令相应的搜索方向为\(\boldsymbol d_B = - C_B^{-1} C_N \boldsymbol d_N\), 令\(\boldsymbol d = (\boldsymbol d_B^T, \boldsymbol d_N^T)^T\)为搜索方向, 进行一维搜索,搜索要限制在可行域内。 在每次迭代中,都重新把\(\boldsymbol x\)拆分为\(\boldsymbol x_B\)\(\boldsymbol x_N\)\(\boldsymbol x_B\)部分要求分量都取正值。 利用这样的方法对原目标函数\(f(\boldsymbol x)\)进行迭代搜索, 可以保证每次迭代延可行下降方向搜索, 当满足KKT条件时停止。 这样的方法叫做简约梯度法。 详见 徐成贤, 陈志平, and 李乃成 (2002) §6.1.2。

投影梯度法和简约梯度法优点是比较简单, 缺点是仅利用梯度信息而没有试图采用高阶逼近, 收敛速度慢。

37.3.3 有效集方法

有效集方法可以看成是投影梯度法的推广。 对问题(37.3), 假设已有约束最小值点的近似值\(\boldsymbol x^{(t)}\)是可行点(满足约束条件), 设\(\boldsymbol x^{(t)}\)处不等式约束中\(c_i, i=p+1, \dots, p+r\)是起作用约束, \(c_i, i=p+r+1, \dots, p+q\)不起作用。 如前所述,搜索方向必须与\(\boldsymbol a_i, i=1, \dots, p+r\)正交, 仍按(37.4)定义\(\mathcal B\)为这些方向的集合。 投影梯度法是把负梯度向量投影到线性子空间\(\mathcal B\)上作为搜索方向, 有效集方法推广了这种方法,仍要求搜索方向在\(\mathcal B\)中, 但不限于负梯度投影方向,即在第\(t+1\)步求解如下的仅含线性约束的子问题 \[\begin{align} \begin{cases} & \mathop{\text{argmin}}\limits_{\boldsymbol d \in \mathbb R^d } f(\boldsymbol x^{(t)} + \boldsymbol d), \ \text{s.t.} \\ & \boldsymbol a_i^T \boldsymbol d = 0, \ i=1,\dots, p+r \end{cases} \tag{37.6} \end{align}\] 这个子问题可以比较容易地求解。

设解(37.6)得到\(\boldsymbol d^{(t)}\)。 如果\(\boldsymbol d^{(t)} = \boldsymbol 0\), 这说明仅含等式约束的子问题 \[\begin{align} \begin{cases} & \mathop{\text{argmin}}\limits_{\boldsymbol x \in \mathbb R^d } f(\boldsymbol x), \ \text{s.t.} \\ & \boldsymbol a_i^T \boldsymbol x = b_i, \ i=1,\dots, p+r \end{cases} \tag{37.7} \end{align}\]\(\boldsymbol x^{(t)}\)为最小值点, 由定理, 可以解出拉格朗日乘子\(\lambda_i^{(t)}, i=1,\dots,p+r\)使得 \[\begin{align*} \sum_{i=1}^{p+r} \lambda_i^{(t)} \boldsymbol a_i = -\nabla f(\boldsymbol x) . \end{align*}\] 如果\(\lambda_i^{(t)} \geq 0, i=p+1,\dots,p+r\), 只要令\(\lambda_i^{(t)} = 0, i=p+r+1, \dots, p+q\), \(\boldsymbol\lambda^{(t)} = (\lambda_1^{(t)}, \dots, \lambda_{p+q}^{(t)})^T\), 则由定理可知 \((\boldsymbol x^{(t)}, \boldsymbol\lambda^{(t)})\)是问题(37.3)的KKT对, 只要判断\(\boldsymbol x^{(t)}\)是不是极小值点。

如果\(\boldsymbol d^{(t)} = \boldsymbol 0\)但是解出的\(\lambda_i^{(t)}, i=p+1,\dots,p+r\)中有负值, 设\(\lambda_{i_0}^{(t)} = \min\{\lambda_i^{(t)}, i=p+1,\dots,p+r \}\), 只要从子问题(37.6)的约束中删去第\(i_0\)个重新求解。

如果子问题(37.6)的解\(\boldsymbol d^{(t)} \neq \boldsymbol 0\), 这时如果\(\boldsymbol x^{(t)} + \boldsymbol d^{(t)}\)是可行点, 满足所有的等式约束和不等式约束(包括不起作用的不等式约束), 则令\(\boldsymbol x^{(t+1)} = \boldsymbol x^{(t)} + \boldsymbol d^{(t)}\), 否则就把\(\boldsymbol d\)当作一个可行下降方向做一维搜索, 一维搜索时要注意不起作用约束也要满足。

具体算法从略,详见 徐成贤, 陈志平, and 李乃成 (2002) §6.3.1。

37.4 二次规划问题

在仅含线性约束的非线性规划问题中, 如果\(f(\boldsymbol x)\)是二次多项式函数, 称这样的问题为二次规划问题, 这是最简单的非线性规划问题, 这种问题有很多实际应用, 另外, 非线性约束优化问题的求解也往往需要通过求解一系列二次规划子问题实现。

37.4.1 仅含等式约束的二次规划问题

考虑如下的二次规划问题: \[\begin{align} \begin{cases} \mathop{\text{argmin}}\limits_{\boldsymbol x \in \mathbb R^d} q(\boldsymbol x) \stackrel{\triangle}{=} \frac12 \boldsymbol x^T H \boldsymbol x + \boldsymbol g^T \boldsymbol x, \quad \text{s.t.}\\ A^T \boldsymbol x = \boldsymbol b, \end{cases} \tag{37.8} \end{align}\] 其中\(H\)\(d\)阶实对称矩阵, \(\boldsymbol g \in \mathbb R^d\), \(A\)\(n\times p\)列满秩矩阵,各列为\(\boldsymbol a_i, i=1,\dots, p\)\(\boldsymbol b \in \mathbb R^p\)。 此问题可以用§37.2的方法求解。

\(A\)有如下QR分解 \[\begin{align*} A = Q \left( \begin{array}{c} R \\ 0 \end{array} \right) = (Q_1, Q_2) \left( \begin{array}{c} R \\ 0 \end{array} \right) = Q_1 R, \end{align*}\] 其中\(Q\)\(d\times d\)正交阵, \(R\)\(p\times p\)满秩上三角阵, \(Q_1\)\(d \times p\)矩阵, \(Q_1^T Q_1 = I_p\), \(Q_1^T Q_2 = \boldsymbol 0\)。 满足约束的可行点\(\boldsymbol x\)的通解为 \[\begin{align} \boldsymbol x = Q_1 R^{-T} \boldsymbol b + Q_2 \boldsymbol v, \quad \boldsymbol v \in \mathbb R^{d-p}, \tag{37.9} \end{align}\]\[\begin{align} \tilde q(\boldsymbol v) = q(Q_1 R^{-T} \boldsymbol b + Q_2 \boldsymbol v) \stackrel{\triangle}{=} \frac12 \boldsymbol v^T \tilde H \boldsymbol v + \tilde{\boldsymbol g}^T \boldsymbol v + \tilde c, \quad \boldsymbol v \in \mathbb R^{d-p}, \tag{37.10} \end{align}\] 其中 \[\begin{align} \tilde H =& Q_2^T H Q_2, \quad \tilde{\boldsymbol g} = Q_2^T \boldsymbol g + Q_2^T H Q_1 R^{-T} \boldsymbol b, \end{align}\] \(\tilde c\)为与\(\boldsymbol v\)无关的常数项。 只要求解无约束优化问题\(\mathop{\text{argmin}} \tilde q(\boldsymbol v)\)得到\(\boldsymbol v^*\), 再用(37.9)就可以得到问题(37.8)的最小值点\(\boldsymbol x^*\)

(37.10)的目标函数\(\tilde q(\boldsymbol v)\)中, 海色阵\(\tilde H\)如果有负特征值, 则\(\tilde q(\boldsymbol v)\)的最小值是\(-\infty\), 问题无解。 事实上,设有\(\lambda<0\)以及\(\boldsymbol d \neq \boldsymbol 0\)使得\(\tilde H \boldsymbol d = \lambda \boldsymbol d\), 取\(\boldsymbol v^{(k)} = k \boldsymbol d\), 则 \[\begin{align} \tilde q(\boldsymbol v^{(k)}) = \frac12 \lambda \| \boldsymbol d \|^2 \cdot k^2 + \tilde{\boldsymbol g}^T \boldsymbol d \cdot k \to -\infty, \quad \text{当}k \to \infty. \end{align}\] 所以,不妨设\(\tilde H\)为非负定阵。

如果\(\tilde H\)是正定阵, 则\(\tilde q(\boldsymbol v)\)是严格凸函数, 有全局严格最小值点\(\boldsymbol v^* = - \tilde H^{-1} \tilde{\boldsymbol g}\), 代入(37.9)可得问题(37.8)的最小值点\(\boldsymbol x^*\)。 这时,再考虑拉格朗日函数 \[\begin{align} L(\boldsymbol x, \boldsymbol\lambda) = q(\boldsymbol x) - \boldsymbol\lambda^T (A^T \boldsymbol x - \boldsymbol b), \end{align}\] 来求拉格朗日函数的稳定点, 其中\(\boldsymbol x^*\)已知, 由定理知拉格朗日乘子\(\boldsymbol\lambda^*\)满足 \[\begin{align*} \nabla q(\boldsymbol x^*) - A \boldsymbol\lambda^* =& 0, \\ A \boldsymbol\lambda^* =& H \boldsymbol x^* + \boldsymbol g, \end{align*}\]\(A = Q_1 R\)\[\begin{align} R \boldsymbol\lambda^* = Q_1^T (H \boldsymbol x^* + \boldsymbol g), \end{align}\] 用回代法解此方程可得拉格朗日乘子\(\boldsymbol\lambda^*\)

如果\(\tilde H\)仅为非负定阵而不是正定阵, \(\tilde q(\boldsymbol v)\)仍然是凸函数, 最小值点与稳定点等价, 稳定点存在当且仅当\(\tilde H \boldsymbol v = -\tilde{\boldsymbol g}\)有解, 这等价于\(\tilde{\boldsymbol g}\)属于\(\mu(\tilde H)\)\(\mu(\tilde H)\)\(\tilde H\)的各列张成的线性子空间), 又等价于 \[\begin{align} (I - \tilde H \tilde H^+) \tilde{\boldsymbol g} = 0, \tag{37.11} \end{align}\] 其中\(\tilde H^+\)\(\tilde H\)的加号逆。 由定理, 当(37.11)成立时\(\tilde q(\boldsymbol v)\)有无穷多个最小值点, 可以表达为 \[\begin{align} \boldsymbol v^* = - \tilde H^+ \tilde g + (I - \tilde H^+ \tilde H) \boldsymbol u, \quad \forall \boldsymbol u \in \mathbb R^{n-p} , \end{align}\] 再代入(37.9)可得问题(37.8)的无穷多个最小值点\(\boldsymbol x^*\)

在问题(37.8)规模很大的时候, 求解线性方程组计算量很大, 可以用共轭梯度法迭代地求解, 由于二次目标函数\(q(\boldsymbol x)\)的特殊性, 迭代最多只需要\(n-p\)次就可以收敛到最小值点。 参见 徐成贤, 陈志平, and 李乃成 (2002) §6.4.1。

37.4.2 含有不等式约束的二次规划问题

考虑如下的二次规划问题: \[\begin{align} \begin{cases} \mathop{\text{argmin}}\limits_{\boldsymbol x \in \mathbb R^d} q(\boldsymbol x) \stackrel{\triangle}{=} \frac12 \boldsymbol x^T H \boldsymbol x + \boldsymbol g^T \boldsymbol x, \quad \text{s.t.}\\ \boldsymbol a_i^T \boldsymbol x = b_i, \ i=1,\dots,p \\ \boldsymbol a_i^T \boldsymbol x \leq b_i, \ i=p+1,\dots,p+q \\ \end{cases} \tag{37.12} \end{align}\] 其中\(H\)\(d\)阶实对称矩阵, \(\boldsymbol g \in \mathbb R^d\), \(\boldsymbol a_i \in \mathbb R^d, i=1, \dots, p+q\), \(b_i \in \mathbb R, i=1,\dots,p+q\)

用前面所述关于一般线性约束的有效集方法可以把含有不等式约束的二次规划问题转化为一系列仅含等式约束的二次规划问题。 假设经过\(t\)步迭代得到问题(37.12)的一个可行点\(\boldsymbol x^{(t)}\)\(\boldsymbol x^{(t)}\)处不等式约束中\(c_i, i=p+1, \dots, p+r\)是起作用约束, \(c_i, i=p+r+1, \dots, p+q\)不起作用。 考虑如下只有等式约束的子问题 \[\begin{align} \begin{cases} \mathop{\text{argmin}}\limits_{\boldsymbol x \in \mathbb R^d} \frac12 \boldsymbol x^T H \boldsymbol x + \boldsymbol g^T \boldsymbol x, \quad \text{s.t.}\\ \boldsymbol a_i^T \boldsymbol x = b_i, \ i=1,\dots,p+r \\ \end{cases} \tag{37.13} \end{align}\]\(\boldsymbol x = \boldsymbol x^{(t)} + \boldsymbol d\), 问题(37.13)可以等价地写成 \[\begin{align} \begin{cases} \mathop{\text{argmin}}\limits_{\boldsymbol d \in \mathbb R^d} \frac12 \boldsymbol d^T H \boldsymbol d + (\boldsymbol g + H \boldsymbol x^{(t)})^T \boldsymbol d, \quad \text{s.t.}\\ \boldsymbol a_i^T \boldsymbol d = 0, \ i=1,\dots,p+r \\ \end{cases} \tag{37.14} \end{align}\]\(\mathcal B = \{ \boldsymbol d \in \mathbb R^d: \boldsymbol a_i^T \boldsymbol d = 0, i=1,\dots, p+r \}\), 如果以\(\boldsymbol a_i, i=1, \dots, p+r\)为列向量组成矩阵\(A\), 则\(\mathcal B = \mu(A)^\perp\)

通过前面对仅含等式约束的二次规划问题的讨论可以知道, 如果 \[\begin{align} \boldsymbol d^T H \boldsymbol d >0, \ \forall \boldsymbol d \in \mathcal B \backslash \{\boldsymbol 0 \}, \tag{37.15} \end{align}\] 则子问题(37.14)存在唯一的全局严格最小值点, 设其为\(\boldsymbol d^{(t)}\)。 如果\(\boldsymbol d^{(t)}= \boldsymbol 0\), 就说明\(\boldsymbol x^{(t)}\)已经是子问题(37.13)的约束最小值点, 可以解出相应的拉格朗日乘子\(\lambda_i^{(t)}, i=1,\dots, p+r\)。 如果解出的\(\lambda_i^{(t)} \geq 0, i=p+1,\dots, p+r\), 则令\(\lambda_i^{(t)} = 0, i=p+r+1, \dots, p+q\), \(\boldsymbol\lambda^{(t)} = (\lambda_1^{(t)}, \dots, \lambda_{p+q}^{(t)})^T\), 这时\((\boldsymbol x^{(t)}, \boldsymbol\lambda^{(t)})\)为原始问题(37.12)的KKT对, 由于二次目标函数的特殊性以及(37.15)条件可知\(\boldsymbol x^{(t)}\) 是问题(37.12)的约束全局严格最小值点。

如果\(\boldsymbol d^{(t)} = \boldsymbol 0\)但是\(\lambda_{i_0}^{(t)} = \min \{\lambda_i^{(t)}: i=p+1, \dots, p+r \} < 0\), 把对应于\(i_0\)的约束从子问题(37.14)中删除, 可以证明新的子问题有非零解\(\tilde{\boldsymbol d}^{(t)}\), 且\(\tilde{\boldsymbol d}^{(t)}\)是原问题(37.12)的可行下降方向。

如果\(\boldsymbol d^{(t)} \neq \boldsymbol 0\), 则\(\boldsymbol d^{(t)}\)是原问题(37.12)的可行下降方向。

\(\boldsymbol d^{(t)}\)是原问题(37.12)的可行下降方向。 若\(\boldsymbol x^{(t)} + \boldsymbol d^{(t)}\)是原问题(37.12)的可行点, 直接取为\(\boldsymbol x^{(t+1)}\)继续迭代即可。 如果\(\boldsymbol x^{(t)} + \boldsymbol d^{(t)}\)超出了(37.12)的可行域, 这一定是不起作用的不等式约束中的某些被突破了, 可以从\(\boldsymbol x^{(t)}\)出发沿\(\boldsymbol d^{(t)}\)方向进行线性搜索, 由二次目标函数的特点知道线性搜索的最小值点在可行域边界处达到, 直接取步长为 \[\begin{align} \alpha_t = \min\left\{ \frac{b_i - \boldsymbol a_i^T \boldsymbol x^{(t)}}{\boldsymbol a_i^T \boldsymbol d^{(t)}} : \ p+r+1 \leq i \leq p+q \text{且} \boldsymbol a_i^T \boldsymbol d^{(t)} > 0 \right\} \tag{37.16} \end{align}\] 注意到\(\boldsymbol d^{(t)}\)已经是子问题(37.14)的约束最小值点, 所以如果得到的\(\alpha_t > 1\)只要取\(\alpha_t = 1\)。 如果\(\alpha_t<1\), 设(37.16)的最小值是在第\(i_1\)个约束处达到的, 则\(\boldsymbol x^{(t+1)} = \boldsymbol x^{(t)} + \alpha_t \boldsymbol d^{(t)}\)\(\boldsymbol x^{(t)}\)增加了一个起作用的不等式约束\(i_1\), 下一步迭代时子问题(37.14)的约束中需要增加约束\(\boldsymbol a_{i_1}^T \boldsymbol x = b_{i_1}\)

详细算法从略。

37.5 非线性约束优化问题

当约束中含有非线性函数时, 问题比无约束和线性约束问题都要复杂得多。 这时, 很难求得可行的下降方向, 而且由于数值误差的因素, 对于不等式约束很难判断是否起作用。 所以一般的非线性约束问题没有通用的解决方法, 更没有可靠的软件。 来自于线性约束优化的一些做法仍可以采用, 但是算法会很复杂。 比较容易实现的方法是各类罚函数法, 定义一系列增加了惩罚项的目标函数, 用一系列无约束优化问题的解逼近约束问题的解, 这样的方法又称为序列无约束最小化方法 (Sequential Unconstrained Minimization Teckniques, SUMT), 特点是实现简单, 但是收敛慢而且惩罚很重时问题适定性很差, 所以难以求得比较精确的结果。 另一类方法是在每步迭代时构造一个二次规划子问题作为近似, 通过求解一系列二次规划问题逼近一般非线性约束问题的解, 这种方法称为序列二次规划法(Sequential Quadratic Programming, SQP), 是现在较好的方法。

37.5.1 外点罚函数法

罚函数法引入包括原来的目标函数和对偏离约束的惩罚两部分的新目标函数, 对新目标函数求无约束最小值点。

对一般的约束优化问题(34.1), 等式约束可以用\(| c_i(\boldsymbol x) |\)惩罚, 不等式约束可以用\(\max(0, c_i(\boldsymbol x))\)惩罚, 取惩罚函数为 \[\begin{align} \bar p(\boldsymbol x) = \sum_{i=1}^p c_i^2(\boldsymbol x) + \sum_{i=p+1}^{p+q} \left[ \max(0, c_i(\boldsymbol x)) \right]^2, \end{align}\] 当且仅当\(\boldsymbol x\)是可行点时\(\bar p(\boldsymbol x)=0\)。 定义新的目标函数 \[\begin{align} \bar f(\boldsymbol x) = f(\boldsymbol x) + \sigma \bar p(\boldsymbol x), \end{align}\] 其中\(\sigma>0\)称为罚因子。 若\(\bar{\boldsymbol x}\)\(\bar f(\boldsymbol x)\) 的一个无约束极值点且\(\bar{\boldsymbol x}\)(34.1)的可行点, 则对(34.1)的任意可行点\(\tilde{\boldsymbol x}\), 都有 \[\begin{align} \bar f(\bar{\boldsymbol x}) =& f(\bar{\boldsymbol x}) + \sigma \bar p(\bar{\boldsymbol x}) = \min \big[ f(\boldsymbol x) + \sigma \bar p(\boldsymbol x) ] \leq f(\tilde{\boldsymbol x}) + \sigma \bar p(\tilde{\boldsymbol x}) = f(\tilde{\boldsymbol x}) \end{align}\] 即只要无约束问题\(\mathop{\text{argmin}} \bar f(\boldsymbol x)\)的解是约束问题的可行点就是约束极小值点。 \(\sigma\)越大,对无约束极小值点不在可行域内的惩罚越大, 所以\(\mathop{\text{argmin}} \bar f(\boldsymbol x)\)也越可能落在可行域内。 算法迭代进行,只要找到的无约束解不可行就增大罚因子\(\sigma\)然后继续求解无约束极值点。 算法如下:

算法实现时可以取\(\sigma_0=1\)\(c=10\)。 求解无约束问题时迭代停止的条件可以取为梯度向量长度\(\| \nabla \bar f(\boldsymbol x) \|\)小于预定精度值。

运用这样的算法, 可以任意选取初始点, 不需要初始点在可行域内, 每次迭代得到的近似解都在可行域外, 只有最后一个解近似地落入可行域而且是近似约束最小值点, 所以这种方法称为外点罚函数法。 这种方法要求目标函数在\(\mathbb R^d\)上都有定义, 如果目标函数仅在与可行域有关的子集上有定义则无法使用; 最后的解不一定严格满足可行性条件, 如果要求解严格可行此种方法也无法使用。 另外,随着\(\sigma\)的增大, 在新目标函数\(\bar f(\boldsymbol x)\)中目标函数占的比例越来越小, 约束惩罚占的比例越来越大, 使得优化问题的适定性变得越来越差, 求解无约束优化问题会变得很困难, 计算误差变得很大。 所以,罚函数法虽然看起来很简单, 但是不一定能得到好的结果。

37.5.2 内点罚函数法

内点罚函数法要求每次迭代严格地在可行域内进行。 定义新的无约束目标函数, 当近似极小值点接近可行域边界时新目标函数变得很大(通常在可行域边界处趋于无穷大), 相当于在可行域边界处设立了无限高的墙壁不允许搜索越过可行域边界, 所以这种方法又称为障碍罚函数法。 内点罚函数法不能处理含有等式约束的问题, 而且要求可行域含有内点。 对于如下仅含不等式约束的优化问题 \[\begin{align} \begin{cases} & \mathop{\text{argmin}}\limits_{\boldsymbol x \in \mathbb R^d } f(\boldsymbol x), \ \text{s.t.} \\ & c_i(\boldsymbol x) \leq 0, \ i=1, \dots, q \end{cases} \tag{37.17} \end{align}\] 定义罚函数为\(\bar p(\boldsymbol x) = -\sum_{i=1}^q \log |c_i(\boldsymbol x)|\)\(\bar p(\boldsymbol x) = \sum_{i=1}^q \frac{1}{|c_i(\boldsymbol x)|}\)\(\boldsymbol x\)只能是严格可行点, 即所有\(c_i(\boldsymbol x)<0\), 当某个\(c_i(\boldsymbol x)=0\)\(\bar p(\boldsymbol x) = +\infty\)。 取新的目标函数 \[\begin{align} \bar f(\boldsymbol x) = f(\boldsymbol x) + \sigma \bar p(\boldsymbol x), \end{align}\] 因为真正的约束最小值点允许取边界处的值, 需要取\(\{ \sigma_t \}\)使得\(\sigma_t \to 0\), 如此减轻对接近可行域边界的惩罚。 内点罚函数法和外点罚函数法类似, 先取较大的\(\sigma\), 求\(\bar f(\boldsymbol x)\)的无约束最小值点, 每次迭代减小\(\sigma\)的值后求无约束极值点, 直到近似约束最小值点达到足够精度。

由于可行域边界处的高墙的存在, 内点罚函数法在近似点靠近可行域边界时加罚的目标函数\(\bar f(\boldsymbol x)\)也会变得病态, 难以求得最小值点, 并且求\(\bar f(\boldsymbol x)\)的最小值点时要注意搜索不能跳出可行域。 内点罚函数法要求初值在可行域内, 可以设法先找可行点, 另一种办法是把内点罚函数和外点罚函数结合使用, 对等式约束和初始点不满足的不等式约束, 可以采用外点罚函数, 其它不等式约束采用内点罚函数, 这样的方法称为混合罚函数法

乘子罚函数法针对普通罚函数方法的病态问题做了改进, 参见 高立 (2014) §7.3、§7.4。

37.5.3 序列二次规划法(SQP)

序列二次规划法(SQP方法)每一步迭代解决一个二次规划子问题, 这种方法在中小规模的非线性约束规划问题中是比较好的方法。 SQP的思想类似于无约束规划中的牛顿法, 在局部对目标函数用二次函数逼近, 对约束用线性函数逼近, 得到二次规划子问题。

对一般约束优化问题(34.1), 设\(f\)有二阶连续偏导数, 设各\(c_i\)有一阶连续偏导数。 设\(\boldsymbol x^{(t)}\)为迭代\(t\)步后得到的近似约束最小值点, 在\(\boldsymbol x^{(t)}\)处对\(f(\boldsymbol x)\)作如下的二阶泰勒展开近似: \[\begin{align*} f(\boldsymbol x^{(t)} + \boldsymbol d) \approx& f(\boldsymbol x^{(t)}) + [ \nabla f(\boldsymbol x^{(t)})]^T \boldsymbol d + \frac12 \boldsymbol d^T [ \nabla^2 f(\boldsymbol x^{(t)})]^T \boldsymbol d, \ \boldsymbol d \in \mathbb R^d, \end{align*}\]\(c_i(\boldsymbol x)\)作如下的一阶泰勒展开近似: \[\begin{align*} c_i(\boldsymbol x^{(t)} + \boldsymbol d) \approx& c_i(\boldsymbol x^{(t)}) + [ \nabla c_i(\boldsymbol x^{(t)})]^T \boldsymbol d, \ \boldsymbol d \in \mathbb R^d. \end{align*}\] 进行这样的近似后, 考虑如下的二次规划子问题 \[\begin{align} \begin{cases} \mathop{\text{argmin}}\limits_{\boldsymbol d \in \mathbb R^d} \ \frac12 \boldsymbol d^T [ \nabla^2 f(\boldsymbol x^{(t)})]^T \boldsymbol d + [ \nabla f(\boldsymbol x^{(t)})]^T \boldsymbol d, \quad \text{s.t.} \\ c_i(\boldsymbol x^{(t)}) + [ \nabla c_i(\boldsymbol x^{(t)})]^T \boldsymbol d = 0, i=1,\dots, p, \\ c_i(\boldsymbol x^{(t)}) + [ \nabla c_i(\boldsymbol x^{(t)})]^T \boldsymbol d \leq 0, i=p+1,\dots, p+q \\ \end{cases} \tag{37.18} \end{align}\] 通过求解这样的二次规划子问题得到下一个近似点\(\boldsymbol x^{(t+1)}\)

SQP算法需要考虑的问题比较多, 这里略去详细的讨论, 感兴趣的读者可以参考 高立 (2014) 第九章。

习题

习题1

用外点罚函数法求解如下约束优化问题: \[\begin{align*} \begin{cases} \mathop{\text{argmin}} x_1^2 + x_2^2, \quad \text{s.t.}\\ x_1 - x_2 + 1 = 0 . \end{cases} \end{align*}\]

习题2

用外点罚函数法求解如下约束优化问题: \[\begin{align*} \begin{cases} \mathop{\text{argmin}} x_1^2 + x_2^2, \quad \text{s.t.}\\ -x_1 + x_2 - 1 \geq 0 . \end{cases} \end{align*}\]

习题3

用内点罚函数法求解上一题目。

References

徐成贤, 陈志平, and 李乃成. 2002. 近代优化方法. 科学出版社.

高立. 2014. 数值最优化方法. 北京大学出版社.