二次规划笔记

最近碰到了一个二次规划问题,第一次自己写 solver 了。二次规划的一般形式是

\displaystyle \min \frac{1}{2} x^\top G x - c^\top x s.t. A^\top x = b, C^\top x \geq d

这个问题的 Lagrange 函数为

\displaystyle L(x, \alpha, \beta) = \frac{1}{2} x^\top G x - c^\top x - \alpha^\top (A^\top x - b) - \beta^\top (C^\top x - d)

这个问题我们可以先考虑只含有等式约束的情形。此时我们根据 KKT 条件需要满足以下方程组

\displaystyle \begin{pmatrix} G & -A \\ A^\top & 0 \end{pmatrix} \begin{pmatrix} x \\ \alpha\end{pmatrix} = \begin{pmatrix} c \\ b\end{pmatrix}

也就是说原问题等价于一个线性方程组求解,我们这里仅仅介绍将此问题转换为 CG 可以 handle 的形式,我们可以假定迭代更新的时候 x\gets x + p,这样变化量 p 代换之后就得到

\displaystyle \begin{pmatrix} G & A \\ A^\top & 0 \end{pmatrix} \begin{pmatrix} -p \\ \alpha\end{pmatrix} = \begin{pmatrix} Gx + c \\ A^\top x - b\end{pmatrix} \stackrel{\triangle}{=} \begin{pmatrix} g \\ h \end{pmatrix}

这是一个对称矩阵,CG 就可以用了,但是直接用也不是很合适,因为满足约束的解一定是在 A 列空间的正交补生成的子空间里面,为此我们 CG 搜索的时候只需要投影在此子空间即可。计算投影矩阵并不困难,如 P = I - A(A^\top A)^{-1} A 即为一个可以使用的结果。我们将标准的 CG 这样修改后(一般称为 PCG)问题就得到了求解。

含不等式的问题一般有两种策略,其一就是尽量使用这个 ECQP(equality-constrained quadratic programming)求解,因为不等式约束碰到边界后就变成了等式约束,这就是常见的 active set 的思路,所谓的 active set 就是指不等式约束中处于 = 状态的集合。这样我们从一个可行解出发,使用 ECQP solver 迭代的更新 x,一旦发现更新碰到边界就提前结束 PCG,这时检查看看那条边界碰到了,将对应的不等式加入到 active set 里面,这样此时的 ECQP 就出现了变化(等式约束变多),于是继续 ECQP。如果 ECQP solver 找不到更新的步长,那么有两种可能,如果对应不等式的 Lagrange 乘子全部非负,那说明这是此问题的解;否则说明某个 KKT 条件尚未满足,将对应的(active set 中)等式约束转换为不等式(从 active set 中移除)进行重新 ECQP 即可。

另一种策略就是基于 interior point method 的思路。interior point method 包括三种思路来 formulate 算法,这里仅仅讨论其一使用 primal-dual method 求解。这个方法仍然是使用 KKT 条件,但是首先对不等式约束引入 slack variable 使得我们拿到的都是等式约束,如 C^\top x - d \geq 0 可以被写做 C^\top x -d -y = 0,其中 y \geq 0,这样写出 KKT 条件后,我们需要类似 log-barrier 方式扰动 y_i \lambda = 0 这个部分,使得解落在内部(后面做详细的分析),这样我们实际求解的线性方程组为(暂时不考虑等式约束,因为这部分相当于做 projection 就行了)

\displaystyle F(x, y, \lambda; \sigma, \mu) = \begin{pmatrix} Gx - C \lambda + c \\ C^\top x - y - b \\ Y\Lambda e - \sigma\mu\end{pmatrix} = 0

其中 \mu = \frac{y^\top \lambda}{m}, \sigma\in[0, 1],且 m 是不等式约束的个数。将更新部分代入 x \gets x + \delta x, y \gets y + \delta y, \lambda \gets \lambda + \delta \lambda,获得更新满足的线性方程组

\displaystyle \begin{pmatrix} G & 0 & -C \\ C^\top & - I & 0 \\ 0 & \Lambda & Y\end{pmatrix} \begin{pmatrix} \delta x \\ \delta y \\ \delta \lambda\end{pmatrix} = \begin{pmatrix} Gx - A\lambda + c \\ A^\top x - y - b \\ -\Lambda Y e + \sigma\mu e\end{pmatrix} \stackrel{\triangle}{=} \begin{pmatrix} r_d \\ r_P \\ -\Lambda Y e + \sigma\mu e\end{pmatrix}

剩下的和前面的问题类似就是如何高效的求解这个线性系统。如果我们还希望使用 PCG 类似的技术,可以如下变换

\displaystyle \begin{pmatrix} G & 0 & C \\ 0 & -Y\Lambda & -I \\ C^\top & I & 0\end{pmatrix} \begin{pmatrix} \delta x \\ -\delta y \\ -\delta \lambda \end{pmatrix} = \begin{pmatrix} -r_d \\ \Lambda e - \sigma\mu Y^{-1}e \\ -r_P\end{pmatrix}

获得了搜索方向,剩下就是选择步长,这部分我们也留在 interior point methods 部分介绍。

为了处理形如 l \leq x \leq u 这类约束,我们通常将搜索的结果限制在这个 box 里面,比较粗暴的方式就是将越界的变量设置为边界值,甚至都不做 line search 看看先碰到哪个边界。但是这样不是我们做 line search 的好办法,比较合理的是以此为终点,而连接两者的必然是一些分段线性函数,我们可以在每一段上进行搜索函数极值,如二次规划就会得到一系列的二次函数极值问题,依次求解即可(这称为搜索 Cauchy point),这时对应的 active set 一般会被固定住,而在这一轮优化过程(的子优化问题)中保持不变(在一个子空间里面优化),此时获得的解作为这步搜索的结果。这称为 gradient projection method。

——————
And in your seed shall all the nations of the earth be blessed; because you have obeyed my voice.

Advertisements
二次规划笔记

发表评论

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / 更改 )

Twitter picture

You are commenting using your Twitter account. Log Out / 更改 )

Facebook photo

You are commenting using your Facebook account. Log Out / 更改 )

Google+ photo

You are commenting using your Google+ account. Log Out / 更改 )

Connecting to %s