等式约束QP命题的求解（Solving equality-constrained QP）

在有效集法（Active Set Method）中，在每次迭代中都要求解一个等式约束的 QP 命题。本文主要对这一问题的求解方法进行关注。此外，本文还会涉及到部分关于 Range Space，Null Space 等内容。主要参考的是 Nocedal 的 Numerical Optimization 一书。

等式约束 QP 命题

这里我们考虑如下的一个等式约束 QP 命题：

minxsubject toq(x)=12xTGx+cTxAx=b,(1a)(1b)\begin{align} \min_x \quad &q(x) = \frac12x^TGx+c^Tx \tag{1a} \\ \text{subject to} \quad &Ax = b, \tag{1b}\end{align}

其中，A∈Rm×nA\in \mathbf R^{m\times n} 是等式约束的系数。整个优化命题有 nn 个决策变量（优化变量），mm 条等式约束（m≤nm\le n）。写出该优化命题的 Lagrange 多项式为：

L=12xTGx+cTx−λT(Ax−b) L= \frac12x^TGx+c^Tx - \lambda^T (Ax-b)

其中，λ\lambda 被称为 Lagrange 乘子，也可以称为对偶变量。由此可以写出其 KKT 条件为：

Gx∗+c−ATλ∗=0Ax∗=b Gx^*+c - A^T\lambda^*=0 \\ Ax^*=b

将上式写成矩阵形式为：

[GA−AT0][x∗λ∗]=[−cb] \begin{bmatrix} G&-A^T \\ A & 0\end{bmatrix}\begin{bmatrix} x^* \\ \lambda^* \end{bmatrix} = \begin{bmatrix} -c \\ b \end{bmatrix}

而在实际计算中，我们往往不是直接计算最优点 x∗x^* 而是计算其由某个初始的估计点 xx 向最优点前进的一步 pp，即 x∗=x+px^*=x+p 。因此我们以 pp 而不是 xx 作为优化变量，可以将上面的方程改写为：

[GAAT0][−pλ∗]=[gh](2) \begin{bmatrix} G&A^T \\ A & 0\end{bmatrix}\begin{bmatrix} -p \\ \lambda^* \end{bmatrix} = \begin{bmatrix} g \\ h \end{bmatrix} \tag{2}

其中，

h=Ax−b,g=c+Gx,p=x∗−xh=Ax-b, \quad g = c+Gx, \quad p=x^*-x

(2) 式中的矩阵通常被称作 KKT 矩阵，在一定的假设下（矩阵 AA 行满秩，ZTGZZ^TGZ 正定，其中 ZZ 的列由矩阵 AA 的零空间的基组成），可以证明该矩阵是非奇异的，也就是说存在唯一的解 (x∗,λ∗)(x^*, \lambda^*) 满足 (2) 式。因此，我们下面讨论的问题就主要集中在如何求解 (2) 的方程组。求解这个方程组的方法主要有两大类，一类是直接进行矩阵分解来求解，另一类是通过迭代的方法来求解。这里我们只介绍直接进行矩阵分解的方法。

直接分解的方法求解 KKT 系统

首先我们记

K=[GAAT0]K= \begin{bmatrix} G&A^T \\ A & 0\end{bmatrix}

虽然一般 GG 阵都是对称正定的，但我们可以证明扩展后的 KK 阵是不定的（indefinite）。因此我们无法直接对 KK 阵进行 Cholesky 分解。直接求解 KKT 系统的方法有以下三种，Full Space / Full KKT System 方法，Range Space / Schur-Complement 方法，Null-Space 方法。关于 Range Space 和 Null Space 的介绍，可以参考 这篇文档 ，这篇文档里从几何上两个空间进行了讲解。关于本文的介绍的三种直接分解 KKT 系统进行求解的方法，出来参考 Numerical Optimization 一书之外，还可以参考 这篇文档，结合文档和书来看可能更容易理解。对三种方法在不同情况下的优劣的分析和比较超出了笔者现在的能力，可能会在进一步的学习之后补充相应的内容。

Full Space / Full KKT System 方法

虽然这里的 K 阵不具有正定性，但还有对称性这一特性可以利用，因此我们可以考虑采用对称不定分解（Symmetric Indefinite Factorization）来对 K 阵进行分解。即对于一个一般性的对称矩阵 K，有

PTKP=LBLT(3)P^TKP=LBL^T \tag{3}

其中 PP 是置换矩阵（Permutation Matrix）的引入是为了保证计算过程中的数值稳定性，LL 是单位下三角矩阵，BB 是块对角矩阵，即对角线上是 1×11\times1 或者 2×22\times2 的块矩阵。结合 (2) 式和 (3) 式，同时因为置换矩阵是正交阵，即 PPT=IPP^T=I，可以有

PTKPPT[−pλ∗]=PT[gh]P^TKPP^T\begin{bmatrix} -p \\ \lambda^* \end{bmatrix}=P^T\begin{bmatrix} g \\ h \end{bmatrix}

即

LBLTPT[−pλ∗]=PT[gh]LBL^T P^T\begin{bmatrix} -p \\ \lambda^* \end{bmatrix}=P^T\begin{bmatrix} g \\ h \end{bmatrix}

因此，我们可以通过下面的一系列步骤来求解

solve Lzsolve Bz^solve LTz¯set [−pλ∗]=PT[gh]=z=z^=Pz¯to obtain z;to obtain z^to obtain z¯\begin{align}\text{solve } Lz&=P^T\begin{bmatrix} g \\ h \end{bmatrix} &\text{to obtain } z; \\ \text{solve } B\hat z&=z &\text{to obtain }\hat z \\ \text{solve } L^T\bar z &= \hat z &\text{to obtain }\bar z \\ \text{set } \begin{bmatrix} -p \\ \lambda^* \end{bmatrix} &=P\bar z \end{align}

可以看出，在上面的方程组中，因为 LL 和 LTL^T 都是三角阵，因此计算量都较小；BB 是块对角矩阵，也只需要求解一些小型的 1×11\times1 或者 2×22\times2 的方程组，因此计算量也不大。

Range Space / Schur-Complement 方法

采用这种方法需要保证 K 阵中的 G 阵是正定的，这样我们可以 (2) 式中的第一个方程得到 pp 与 λ\lambda 之间的关系，然后代入第二个方程中进行求解。即向第一个方程的等号左右两边同乘 AG−1AG^{-1} 可以得到

AG−1ATλ∗=AG−1g+ApAG^{-1}A^T\lambda^*=AG^{-1}g+Ap

将第二个方程代入即可得：

AG−1ATλ∗=AG−1g−hAG^{-1}A^T\lambda^*=AG^{-1}g-h

其中，AG−1ATAG^{-1}A^T 是对称正定阵，因此可以采用 Cholesky 等较为高效的分解算法来进行分解，提高方程组求解的效率。当通过上式求解得到 λ∗\lambda^* 之后，可以通过求解下面的方程来还原 pp 。

Gp=ATλ∗−gGp=A^T\lambda^*-g

在该方法中，因为需要对 GG 阵进行求逆，同时需要对 AG−1ATAG^{-1}A^T 进行分解，因此该方法更适合以下几种情况：

GG 阵条件数较小且易于求逆，比如 GG 是对角阵或者分块对角阵，或者

GG 阵的逆可以用其他方法显式地表达出来；

原优化命题的约束数目 mm 较小，这样需要分解的 AG−1ATAG^{-1}A^T 阵的维数就较小。

另外，Schur-Complement 这个名称来源于线性代数中的一个术语，因为如果我们对原矩阵 KK 以 GG 阵为轴（pivot）进行分块的高斯消去，就可以得到下式

[G0AT−AG−1AT]\begin{bmatrix} G & A^T \\ 0 & -AG^{-1}A^T \end{bmatrix}

其中，AG−1ATAG^{-1}A^T 就是矩阵 KK 中 GG 阵的 Schur-Complement。

Null-Space 方法

这里我们不加深入讨论地引入关于 Null Space 的一些理论和方法，关于这些方法的进一步讨论，可以参看原书，我也会在后续的另一篇博文中进行更深入地讨论。

Null Space 方法不需要 GG 正定的假设，因此有更广的适用面，同时也不需要对 GG 进行求逆，但是需要了解 AA 阵的零空间（Null Space）的基矩阵 ZZ 的信息。

假如我们把向量 pp 划分成如下的两部分，即

p=YpY+ZpZ(4)p=Yp_Y+Zp_Z \tag4

其中，Z∈Rn×(n−m)Z\in\mathbf R^{n\times (n-m)} 是零空间矩阵，Y∈Rn×mY\in\mathbf R^{n\times m} 是任意可以使 [Y|Z][Y|Z] 非奇异的矩阵。pYp_Y 是长度为 mm 的向量，pZp_Z 是长度为 n−mn-m 的向量。将 (4) 式代入到 (2) 式的第二条方程中，同时因为有 AZ=0AZ=0，我们可以得到

(AY)pY=−h(AY)p_Y=-h

可以证明 AYAY 是 m×mm\times m 的非奇异矩阵，因此可以哦通过求解上面的方程组得到 pYp_Y。此外，我们将 (4) 式代入到 (2) 式的第一条方程中去可以得到

−GYpY−GZpZ+ATλ∗=g-GYp_Y-GZp_Z+A^T\lambda^*=g

等式两边同乘 ZTZ^T 可以得到

(ZTGZ)pZ=−ZTGYpY−ZTg(Z^TGZ)p_Z=-Z^TGYp_Y-Z^Tg

上面的方程组可以通过对 ZTGZZ^TGZ 进行 Cholesky 分解来求解。然后我们可以用下面的方程来求解对偶变量 λ∗\lambda^* :

(AY)Tλ∗=YT(g+Gp)(AY)^T\lambda^*=Y^T(g+Gp)

当 n−mn-m 较小时，Null Space 方法会非常有效。这种方法的缺陷在于我们需要计算零空间矩阵 ZZ 而这对于一些规模较大的问题可能会比较麻烦。计算 ZZ 的一种方法是对 AA 阵进行 QR 分解。相关内容可能在后续的博文中介绍。