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