06-等式约束优化算法
凸优化从入门到放弃完整教程地址：https://www.cnblogs.com/nickchen121/p/14900036.html
一、简介
注：这里通过 KKT 条件对等式约束问题的一种引入
前面的系列我们介绍了Gradient Descent和Newton Method求解无约束的凸优化问题，并且给出了带约束问题的可行解的存在条件——KKT条件。在这篇文章里，我们考虑带等式约束的凸优化问题：
(egin{align} min quad &f(x) \ mathrm{s.t.} quad & Ax=b end{align} ag{1})
其中 (f) 二阶可导， (A in R^{p imes n}, mathrm{rank}(A) = p < n) 。根据KKT条件，我们知道 (x^* in dom(f)) 是优化问题(1)的最优解的充要条件是，存在 ( u^* in mathbb{R}^n) 满足
(Ax^* = b quad abla f(x^*) + A^T u^* = 0 ag{2})
所以，求解优化问题(1)等价于求解 (n+p) 个变量组成的方程(2)。一般而言 $ abla f(x^) + A^T u^ = 0$ 是非线性的，很难求出其解析解。但在下面的部分，我们将考虑 (f) 的二阶近似，从而使得 ( abla f(x^*)) 线性化。
二、等式约束凸二次规划
注：此处的引入非常重要，因为上面说到了我们将考虑 (f) 的二阶近似 (hat f)，而 (hat f) 就是一个二次函数
考虑 (f) 为二次函数的情况
(egin{aligned} min quad & f(x) = frac{1}{2}x^TPx + q^Tx + r \ mathrm{s.t.} quad & Ax=b end{aligned} ag{3})
根据凸性， (P in S_+^n, A in mathbb{R}^{p imes n}) 。注：该问题是扩展Newton Method处理等式约束问题的基础。此时最优条件可以写为
(Ax^* = b, quad Px^* + q + A^Tv^* = 0 \)
用矩阵表示为
(Big[egin{array}{cc} P & A^T \ A & 0 end{array} Big] Big[ egin{array}{c} x^* \ u^* end{array}Big] = Big[ egin{array}{c} -q \ b end{array}Big] ag{4})
这 (n+p) 个变量组成的线性方程组称为KKT矩阵。如果KKT矩阵非奇异，存在唯一最优的原对偶对 ((x^*, u^*)) ；如果KKT矩阵奇异，但是有解，任何解都构成最优对偶对 ((x^*, u^*)) ；如果KKT矩阵无解，那么二次优化问题或者无解或者无下界。
注：这里用到了矩阵乘法 (AX = B) 的知识点，非奇异矩阵可以理解为可逆矩阵
下面指出KKT矩阵非奇异的一个充分条件：(P) 正定。
三、等式约束的Newton方法
前面指出，一般情况下(2)是非线性方程，很难有解析解。为了方便求解，我们对目标函数 (f) 在 (x) 处做二阶近似。
(egin{align} min quad &f(x+v) = f(x) + abla f(x)^Tv + frac{1}{2} v^T abla^2 f(x) v \ mathrm{s.t.} quad & A(x+v)=b quad mathrm{or} quad Av=0 end{align} ag{5})
该问题的变量为 (v) ，且是关于 (v) 的等式约束凸二次规划问题。我们假定KKT矩阵非奇异，在此基础上定义 (x) 处的Newton下降方向 (Delta x_{nt}) 为凸二次问题(5)的解。根据(4)，Newton下降方向 (Delta x_{nt}) 由以下KKT方程决定。注：KKT 方程写成这样主要就是参考公式 (4) 后，对公式 (4) 内的变量进行了一个替换
(Big[egin{array}{cc} abla^2f(x) & A^T \ A & 0 end{array} Big] Big[ egin{array}{c} Delta x_{nt} \ w end{array}Big] = Big[ egin{array}{c} - abla f(x) \ 0 end{array}Big] ag{6})
其中 (w) 是该二次问题的最优对偶变量(防止符号滥用，我们这里没有再使用 ( u) )。所以求解问题(5)等于求解方程组(6)。
首先我们说明这样的一个下降方向 (Delta x_{nt}) 是满足迭代算法要求的。首先，可以看出 (A Delta x_{nt}=0) 满足等式约束；此外， (f) 沿 (Delta x_{nt}) 处的方向导数是小于0的，说明 (Delta x_{nt}) 是一个下降方向。
(frac{d}{dt} f(x+tDelta x_{nt}) Big|_{t=0} = abla f(x)^T Delta x_{nt} = -Delta x_{nt}^T abla^2 f(x) Delta x_{nt} < 0 \)
下面，我们给出等式约束的Newton Method算法的框架，可以看出其和无约束情况下完全一样。
重复进行：
- 计算Newton方向(Delta x_{nt}) ，即求解KKT方程(6)
- 计算Newton减量 (lambda(x) =(Delta x_{nt}^T abla^2 f(x) Delta x_{nt})^{1/2})
- 停止准则：如果 (lambda^2/2 leq epsilon) ，则退出
- 直线搜索：通过回溯直线法确定步长 (t)
- 改进： (x:= x+ t Delta x_{nt})
同样地，Newton Method收敛也存在阻尼Newton阶段和纯Newton阶段。阻尼Newton阶段收敛较慢，但有界；纯Newton阶段收敛十分迅速。下面说明如何求解KKT系统得到 (Delta x_{nt}) 。
四、求解KKT系统
注：有兴趣的可以看看书上的，这里写的有点简单
这里专门介绍求解线性方程组(6)是因为我们可以利用 ( abla^2 f(x)) 的正定性，加速计算。
(Big[egin{array}{cc} abla^2f(x) & A^T \ A & 0 end{array} Big] Big[ egin{array}{c} Delta x_{nt} \ w end{array}Big] = Big[ egin{array}{c} - abla f(x) \ 0 end{array}Big] \)
不失一般性，我们可以直接求解这个线性方程组。不利用矩阵的结构，计算量是 ((1/3)(n+p)^3) 次浮点运算。当 (n) 和 (p) 都不大时，这是一个合理的处理方法。
此外，我们可以采用消元法。假设 ( abla^2 f(x)) 正定，利用KKT方程组第一个方程 ( abla^2 f(x) Delta x_{nt} + A^T w = - abla f(x) \)
可以解出 (Delta x_{nt})
(Delta x_{nt} = - { abla^2 f(x)}^{-1} ( abla f(x) + A^T w) \)
然后将其带入KKT方程组的第二个方程，可以解出 (w)
(w = (A { abla^2 f(x)}^{-1}A^T)^{-1} (-A{ abla^2 f(x)}^{-1} abla f(x)) \)
下面给出基于消元法的求解过程：
- 计算 ({ abla^2 f(x)}^{-1} A^T) 和 ({ abla^2 f(x)}^{-1} abla f(x))
- 计算Schur补 (S = -A { abla^2 f(x)}^{-1} A^T)
- 求解 $Sw = A{ abla^2 f(x)}^{-1} abla f(x) $ 确定 (w)
- 求解 ( abla^2 f(x) Delta x_{nt} = - A^Tw - abla f(x)) 确定 (Delta x_{nt})
采用合适的矩阵分解方法(如Cholesky因式分解)，整体的浮点计算次数大概为：
(f +ps + p^2n + (1/3)p^3)
其中 (f) 是对 ({ abla^2 f(x)}) 进行Cholesky因式分解所需要的计算量。
五、总结
对于带等式约束的凸优化问题，我们将目标函数进行了二次近似，根据KKT条件，确定了最优解的存在条件——KKT方程。然后通过求解KKT方程确定Newton Method需要的下降方向 (Delta x_{nt}) ，并且对快速求解KKT方程做了一定的分析。
参考文献：Stephen Boyd, Lieven Vandenberghe: Convex Optimization