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