QP问题

发表于   |   分类于   |  

QP问题，即带线性约束的二次方程最小化问题。QP问题的一般形式为：
$$\label{origin} \min_x\ q(x) = \frac{1}{2}x^TGx+x^Tc\\ s.t.\ a^T_ix = b_i,\quad i\in\mathcal{E}\\ a^T_ix \ge b_i,\quad i\in\mathcal{I}$$
对于 Hessian 矩阵 $G$ 来说，如果是半正定的，则QP称为凸QP；否则称为非凸QP。显然，非凸QP问题更麻烦一些。不过我们这里只关注凸QP问题。

只有等式约束的QP问题

先从等式约束开始说起。等式约束的QP问题可以表示为：
$$\min_x\ q(x) = \frac{1}{2}x^TGx+x^Tc\\ s.t.\ Ax = b$$
先看A满秩的情况。根据 KKT 条件，方程的解需要满足下面的条件：
$\begin{bmatrix}G & -A^T\\A&0\end{bmatrix}\begin{bmatrix}x^* \\\lambda^* \end{bmatrix} = \begin{bmatrix}-c \\b \end{bmatrix}$
我们令 $x^* = x+p$，对上面的等式稍加变换，就可以得到
$$\label{KKT} \begin{bmatrix}G & A^T\\A&0\end{bmatrix}\begin{bmatrix}-p \\\lambda^* \end{bmatrix} = \begin{bmatrix}g \\h \end{bmatrix}\\ h = Ax-b\\ g = c+Gx\\ p = x^* - x$$
上面的矩阵被称为 KKT 矩阵。令 Z 为 A 的 null space(即 AZ = 0)，维数为 $n\times(n-m)$ 。如果 $Z^TGZ$ 正定，则 KKT 矩阵必然非奇异，也非正定。而满足 KKT 矩阵的解必然为 QP 的解。

KKT 系统求解

所以现在的问题就是如何求解上面的 KKT 系统。我们先来看直接求解的方法。为了方便说明，将 KKT 矩阵记为 $K$

直接求解

schur-complement

直接用 schur-complement 暴力求解。利用 schur-complement 可以直接拿到 $K$ 的逆：
$$\begin{aligned} K^{-1} &= \begin{bmatrix}C & E\\E^T&F\end{bmatrix}\\ C &= G^{-1}-G^{-1}A^T(AG^{-1}A^T)^{-1}AG^{-1}\\ E &= G^{-1}A^T(AG^{-1}A^T)^{-1}\\ F &= -(AG^{-1}A^T)^{-1} \end{aligned}$$
当然了其实你用高斯消元是可以直接解出来的，和这么算的结果也太不太多。不过这么暴力真的好么= =

矩阵分解

由于 KKT 矩阵非正定，无法使用 cholesky 分解，所以用 symmetric indefinite factorization。对于一个对称矩阵 $K$， symmetric indefinite factorization 可以表示为
$P^TKP=LBL^T$
其中 $P$ 为 permutation 矩阵， $L$ 为下三角矩阵，$B$ 为 block-diagonal 矩阵， block 的大小一般为 1*1 或者 2*2.这种情况下，可以将 \ref{KKT} 的解表示为：
$$Lz = P^T\begin{bmatrix}g \\h \end{bmatrix}\\ B\hat{z} = z\\ L^T \bar{z} = \hat{z}\\ \begin{bmatrix}-p \\\lambda^* \end{bmatrix} = P\bar{z}\\$$

null space

\ref{KKT} 中的向量 $p$ 其实可以被分解为2个部分：
$p = Yp_Y + Zp_Z$
其中 $Z$ 为 n(n-m) 的矩阵，是 $A$ 的 null space。 Y 为 nm 的矩阵，保证 [Y|Z] 非奇异。由于 $AZ=0$，所以
$$(AY)p_Y = -h\\ -GYp_Y-GZp_Z+A^T\lambda^* = g$$
将后一个式子乘以 $Z^T$ 可以得到
$(Z^TGZ)p_z = -Z^TGYp_Y-Z^Tg$
$Z^TGZ$ 是对称正定的，所以可以用 cholesky 分解。所以我们可以通过解上面的等式得到 \ref{KKT} 的解。最后
$$p = Yp_Y+Zp_Z\\ (AY)^T\lambda^* = Y^T(g+Gp)$$
当然 Z 的维度越小，算法跑的越快。但是如果 $AG^{-1}A^T$ 很好算的话，还是 schur-complement 暴力硬上更快一些。

迭代式解法

迭代式解法主要用的是 CG。 根据前面的说明，最优解 $x^* = Yx_Y +Zx_Z$
这样的话，根据 $AYx_Y = b$，我们可以直接解出 $x_Y$。不过这涉及到 Y 的选择问题。

对于Y的选择，方法有很多。这里只简单列出两种：

  • A为 m*n 维矩阵，如果 A 满秩，就会有 m 个线性无关组(假设 m < n)。那么就可以找到一个 permutation 矩阵 $P$ 使得 $AP = [B|N]$ 其中 $B$ 为线性无关组所在的矩阵。我们可以令 $Y = \begin{bmatrix}B^{-1}\\0\end{bmatrix}$
  • 对 $A$ 做 QR 分解： $A^T\Pi = \begin{bmatrix}Q_1&Q_2\end{bmatrix} \begin{bmatrix}R\\0\end{bmatrix}$，令 $Y=Q_1$

这种情况下，去掉只和 $x_Y$ 有关的部分，可以将原问题化简为
$$\min_{x_Z} \frac{1}{2}x_Z^TZ^TGZx_Z+x_Z^Tc_Z\\ c_Z = Z^TGYx_y+Z^Tc$$
上面的解很简单：$Z^TGZx_Z=-c_Z$ 所以我们只要用 CG 解这个方程就行了。CG的算法如下：

  1. choose initial point $x_Z$
  2. compute $r_Z=Z^TGZx_Z,\ g_Z=W_{ZZ}^{-1}r_Z,\ d_Z=-g_Z$
  3. repeat $$\begin{aligned} \alpha &= r_Z^Tg_Z/d^T_ZZ^TGZd_Z\\ x_Z &= x_Z+\alpha d_Z\\ r_Z^+ &= r_Z+\alpha Z^TGZd_Z\\ g_Z^+ &= W_{ZZ}^{-1}r_Z^+\\ \beta &= (r_Z^+)^Tg_Z^+/r_Z^Tg_Z\\ d_Z &= -g_Z^+ + \beta d_Z\\ g_Z &= g_Z^+\\ r_Z &= r_Z^+ \end{aligned}$$ until convergence

和CG算法基本一致，不过进行了 precondition。 终止条件可以是 $r_Z^TW_{ZZ}^{-1}r_Z$ 足够小或者其他表示残差小的东西。至于 preconditioner $W_{ZZ}$ ，它本身的存在就是为了防止 $Z^TGZ$ 的条件数过高，所以取值方面自然是要降低条件数。最极限的做法自然是令 $W_{ZZ} = Z^TGZ$， 一下给降到1。当然也可以用缓和点的取法，$W_{ZZ} = Z^THZ$ 其中 $H$ 是对称矩阵保证$W_{ZZ}$ 正定。 $H$ 可以取值为 $diag(|G_{ii}|)$ 或者 $I$ 。有时也可以取值为 $G$ 的 block diagonal submatrix.

算法中其实没有对 $Z$ 单独做操作，所以只要能求出来 $Z^TGZ$ 或者它和其他向量的乘法计算就可以。

但是算法上避免 $Z$ 的计算仍然是可能的。不过我们最后求出来的直接是 $x$ 不是 $x_Z$ 。不过介绍算法之前，需要引入几个记号：
$$\begin{aligned} Z^Tr&=r_Z\\ g&=Zg_Z\\ d&=Zd_Z\\ P&=Z(Z^THZ)^{-1}Z^T \end{aligned}$$
算法叫做 Projected CG：

  1. choose init point $x$ satisfied $Ax=b$
  2. compute $r=Gx+c, g=Pr, d=-g$
  3. repeat $$\begin{aligned} \alpha &= r^Tg/d^TGd\\ x &= x+\alpha d\\ r^+ &= r + \alpha Gd\\ g^+ &= Pr^+\\ \beta &= (r^+)^Tg^+/r^Tg\\ d &= -g^+ + \beta d\\ g &= g^+\\ r &= r^+ \end{aligned}$$ until convergence

终止条件是 $r^Tg - r^TPr$ 很小。当然了，这么搞其实和上面一样，因为 $P$ 的计算需要依靠 $Z$ 。但是，$P$ 其实可以用其他方法计算：
$P = H^{-1}(I-A^T(AH^{-1}A^T)^{-1}AH^{-1})$
这种情况下 $g^+ = Pr^+$ 可以表示为：
$\begin{bmatrix}H & A^T \\ A & 0\end{bmatrix}\begin{bmatrix}g^+\\ v^+\end{bmatrix} = \begin{bmatrix}r^+\\ 0\end{bmatrix}$
对于这个等式，可以直接用前面的矩阵分解来求。

不等式约束的情况

讨论完了只有等式的部分，现在开始处理不等式的部分。我们令 $\mathcal{A}(x^* )=\{i\in \mathcal{E}\cup\mathcal{I}|a_i^Tx^* = b_i\}$。 那么对于 \ref{origin} 来说，利用 KKT 条件可以得到最优解满足的条件：
$$\begin{aligned} Gx^* + c - \sum_{i \in \mathcal{A}(x^* )}\lambda_i^* a_i &= 0\\ a_i^Tx^* &= b_i \quad \text{for all } i \in \mathcal{A}(x^* )\\ a_i^Tx^* &\ge b_i\quad \text{for all } i \in \mathcal{I}\backslash \mathcal{A}(x^* )\\ \lambda_i^* &\ge 0\quad \text{for all } i \in \mathcal{I}\cap \mathcal{A}(x^* ) \end{aligned}$$

当 $G$ 正定时，目标函数是凸函数，没有任何问题；但是如果不满足条件，就比较麻烦了，这时候会出现多个局部最优解。

但是麻烦远不止这一个。有时候会出现 degeracy 的情况。这种情况一般分为两种：

  1. $a_i, i\in \mathcal{A}(x^* )$ 线性相关
  2. 一些 $\lambda_i^* = 0, i\in \mathcal{A}(x^* )$

第一点的问题主要是不正定所带来的数值问题；而第二点的问题主要是无法判断哪些条件属于 active set(即 $\mathcal{A}(x)$)

下面给出3种解法。

active-set method

active-set method 的思想和线性规划的 simplex method 很像，都需要找到最后的 active set。

假设我们在第k步得到了 $x_k$ 和当前的 active set $\mathcal{W}_k$ 。首先要检查 $x_k$ 是否是当前条件下的最小值，如果不是，则需要计算 $x$ 的更新量 $p$：
$$\label{p_min} \min_p \frac{1}{2}p^TGp + g_k^Tp\\ s.t.\ a_i^Tp = 0\quad i\in\mathcal{W}_k$$
如果更新后的 $x_{k+1}$ 满足所有的不等式约束，就可以直接更新，否则我们需要缩小更新量直到所有的约束都满足(所以是有可能一次更新的更新量为0)。对于每个不等式约束来说，都需要满足：
$a_i^T(x_k+\alpha_kp_k)\ge b_i$
因此
$\alpha_k \le \dfrac{b_i-a_i^Tx_k}{a_i^Tp_k}$
对于所有的不等式约束都算一遍，就能得到最后的更新率：
$\alpha_k = \min(1, \min_{i\notin\mathcal{W}_k,a_i^Tp_k<0}\dfrac{b_i-a_i^Tx_k}{a_i^Tp_k})$
所以很显然，当 $\alpha_k$ 小于1的时候，肯定会有不等书约束被 active 。这时我们需要更新 $\mathcal{W}_k$, 加上刚刚 active 的约束。显然这样一直往下加，必然会结束。要不就是约束全部 active， 要不就是找到了最小值。这种情况下的 $\hat{x}$ 和 $\hat{\lambda}$ 是满足 KKT 条件的，所以我们可以计算出部分的 $\hat{\lambda}$。如果计算的结果全部非负，那么找到最优解；如果有负值，说明负值对应的约束不在实际的 active set 中，需要去掉再进行计算。

这样我们就得到了最后的 active-set method：

  1. get init point $x_0$ and get active set $\mathcal{W}_0$
  2. repeat
    • solve \ref{p_min} to get $p_k$
    • if $p_k$ = 0
      • compute $\hat{\lambda}$
      • if $\hat{\lambda}_i>0$ for all $i \in \mathcal{I}\cap \mathcal{W}_k$, stop
      • else$$j = \arg\min_{j\in\mathcal{I}\cap \mathcal{W}_k}\hat{\lambda}_j\\ x_{k+1} = x_k\\ \mathcal{W}_{k+1}=\mathcal{W}_{k}\backslash\{j\}$$
    • else
      • compute $\alpha_k$
      • $x_{k+1} = x_k+\alpha_kp_k$
      • if $\alpha_k < 1$, update $\mathcal{W}_k$
      • $\mathcal{W}_{k+1}$ = $\mathcal{W}_k$

和 simplex 不同的是，该算法不会出现 cycling 的情况，所以必然会在固定步数内结束。

现在我们还剩下最后一步，初值。对于初值的获得，和 simplex method 一样，我们需要一个 Phase I。下面说一些 Phase I 的算法。

Phase I

和 simplex method 一样， phase I 同样是最小化的过程。比如可以用
$$\min_{(x,z)} e^Tz\\ s.t.\ a_i^Tx + \gamma_iz_i = b_i\quad i\in \mathcal{E}\\ a_i^Tx + \gamma_iz_i \ge b_i\quad i\in \mathcal{I}\\ z \ge 0$$
其中 $e=(1,1,...,1)^T,\gamma_i=\begin{cases}-\text{sign}(a^T_ix-b) & i\in\mathcal{E}\\ 1 & i\in\mathcal{I}\end{cases}$

或者可以用 penalty(or bigM) method:
$$\min_{(x, \eta)}\ \frac{1}{2}x^TGx+x^Tc+M\eta\\ s.t. \ a_i^Tx-b_i\le\eta\quad i \in \mathcal{E}\\ -(a_i^Tx-b_i)\le\eta\quad i \in \mathcal{E}\\ b_i-a_i^Tx\le\eta\quad i \in \mathcal{I}\\ 0\le\eta\\$$
当然这种情况下我们又面临 $M$ 的选择问题。当然这个就简单很多。我们的最终目的是让求解的 $\eta = 0$ 所以如果条件不满足，加大 $M$ 就可以了。

与上面的方法类似，可以用 $l_1$ 作为惩罚项：
$$\min_{(x,s,t,v)} \ \frac{1}{2}x^TGx + x^Tc+Me^T_{\mathcal{E}}(s+t) + Me^T_{\mathcal{I}}v\\ s.t.\ a_i^Tx-b_i+s_i-t_i = 0\quad i\in\mathcal{E}\\ a_i^Tx-b_i+v_i \ge 0\quad i\in\mathcal{I}\\ s\ge 0, t\ge 0, v\ge 0$$
$e_{\mathcal{E}},e_{\mathcal{I}}$ 都是全1向量

update 的问题

最后还有一个问题。在算法中，每次迭代我们都需要计算 $p_k$ ，但是其实每次更新的 $\mathcal{W}_k$ 都只是单个条件的加减而已。所以 $p_k$ 的计算其实是可以复用的。由于 $p_k$ 的计算也是带有等式约束的QP问题，所以可以用前面提到的方法去解。这里我们对 null space method 进行讨论。

在 null space method 中。一旦 $Y$ 和 $Z$ 确定，$p$ 就可以确定。一种 $Y$ 和 $Z$ 的取法是对 $A$ 做QR分解：
$A^T\Pi = Q\begin{bmatrix}R\\0\end{bmatrix} = \begin{bmatrix}Q_1&Q_2\end{bmatrix}\begin{bmatrix}R\\0\end{bmatrix}$
其中 $\Pi$ 为 permutation 矩阵，$Q$ 为正交矩阵，$R$ 为非奇异的上三角矩阵。令 $Y=Q_1,\ Z=Q_2$ 我们的更新也从这里入手。
对于增加约束的情况，将约束增加后的矩阵记为 $\bar{A}$, 则 $\bar{A}^T = \begin{bmatrix}A^T&a\end{bmatrix}$ $a$ 为行向量，使 $\bar{A}$ 保持行满秩。则
$\bar{A}^T\begin{bmatrix}\Pi&0\\0&1\end{bmatrix} = \begin{bmatrix}A^T\Pi&a\end{bmatrix} = Q\begin{bmatrix}R&Q_1^Ta\\0&Q_2^Ta\end{bmatrix}$
现令正交矩阵 $\hat{Q}$ 满足：
$\hat{Q}(Q^T_2a) = \begin{bmatrix}\gamma\\0\end{bmatrix}$
则
$\bar{A}^T\begin{bmatrix}\Pi&0\\0&1\end{bmatrix} = Q\begin{bmatrix}R&0\\0&\hat{Q}^T\begin{bmatrix}\gamma\\0\end{bmatrix}\end{bmatrix} = Q\begin{bmatrix}I&0\\0&\hat{Q}^T\end{bmatrix}\begin{bmatrix}R&Q_1^Ta\\0&\gamma\\0&0\end{bmatrix}$
所以
$$\bar{\Pi} = \begin{bmatrix}\Pi&0\\0&1\end{bmatrix}\\ \bar{Q} = Q\begin{bmatrix}I&0\\0&\hat{Q}^T\end{bmatrix} = \begin{bmatrix}Q_1&Q_2\hat{Q}^T\end{bmatrix}\\ \bar{R} = \begin{bmatrix}R&Q_1^Ta\\0&\gamma\end{bmatrix}$$
对于约束的删除来说，则是 $R$ 中的某一行被删除。这样的话 $R$ 直接就不是上三角矩阵了，所以我们需要一个 permutation 矩阵进行行变换以恢复 $R$ 的上三角特征。

gradient projection method

active-set method 在大规模的约束情况下会收敛的很慢。毕竟每次只能增删一个约束，要是碰到个2千约束的，至少都要2000次迭代。gradient projection 就不一样。

我们将 gradient projection 处理的问题描述为：
$$\min_x\ q(x) = \frac{1}{2}x^TGx+x^Tc\\ s.t.\ l\le x \le u$$
其中 $G$ 是对称的，但是不需要正定。

gradient projection 的思路类似于 trust-region。算法被分为2个阶段：首先沿着导数方向进行更新，即 $g=Gx+c$ 。当 $x$ 到达边界时，我们将导数在约束空间内做投影。之后利用得到的约束进行再次计算。

我们将 $x$ 在约束空间内的投影表示为
$P(x,l,u)=\begin{cases}l_i & \text{if } x_i<l_i\\x_i& \text{if } x_i\in[l_i,u_i]\\u_i& \text{if } x_i>u_i\end{cases}$
这种情况下我们得到的导数 $x(t) = P(x-tg,l,u)$ 对于每个约束，我们可以拿到的 $t$ 值为
$\bar{t}_i=\begin{cases}(x_i-u_i)/g_i&\text{if } g_i<0,u_i<\infty \\(x_i-l_i)/g_i&\text{if } g_i>0,l_i>-\infty\\\infty&\text{otherwise }\end{cases}$
为了确定t的值，我们首先要对 $\bar{t}$ 顺序排列，在每个间隔内进行搜索：
$$x(t) = x(t_{j-1}) + (\Delta t)p^{j-1}\\ p^{j-1}_i=\begin{cases}-g_i&\text{if } t_{j-1}<\bar{t}_i\\0&\text{else}\end{cases}$$
搜索目标为最小化
$q(x(t)) = c^T(x(t))+\frac{1}{2}(x(t))^TG(x(t))$
这样的情况下我们就可以得到一组 active-set。更新后的目标为：
$$\min_x \frac{1}{2}x^TGx + x^Tc\\ s.t.\ x_i=x_i^c\quad i\in\mathcal{A}(x^c)\\ l_i\le x_i\le u_i\quad i\notin\mathcal{A}(x^c)$$
不过得到这个问题的解的难度和原问题基本一样= =。当然，大数据量的情况下也不会要求精确解，所以粗略的用CG什么的解下就好了