數值優化：經典二階確定性演演算法與對偶方法

我們在上一篇部落格《數值優化：經典一階確定性演演算法及其收斂性分析》中主要介紹了單機數值優化中一些經典的一階確定性演演算法，本篇文章我們將會介紹二階確定性演演算法和對偶方法。

1 牛頓法

1.1 演演算法描述

牛頓法^[1]的基本思想是將目標函數在當前迭代點處進行二階泰勒展開，然後最小化這個近似目標函數，即

\[\underset{w\in \mathcal{W}}{\text{min}} f(w) \approx \underset{w \in W}{\text{min}} f(w^t) + \nabla f(w^t)^T(w - w^t) + \frac{1}{2}(w - w^t)^T\nabla^2f(w^t)(w-w^t) \]

此處\(\nabla^2f(w^t)\)是目標函數在當前迭代點\(w^t\)處的海森矩陣。如果該海森矩陣是正定的，則上述問題的最優值 \(w^t -[\nabla^2f(w^t)]^{-1}\nabla f(w^t)\)處取到，牛頓法將其做為下一時刻的狀態，即

\[w^{t+1} = w^t -[\nabla^2f(w^t)]^{-1}\nabla f(w^t) \]

下圖給了一個簡單範例：

牛頓法的虛擬碼如下所示：

1.2 收斂性分析

相比於一階方法中的梯度下降法，二階的牛頓法提供了更為精細的步長調節，即利用當前狀態海森矩陣的逆矩陣。因為步長更為精細，牛頓法的收斂速率比梯度下降法的收斂速率顯著加快，具有二次收斂速率。

牛頓法收斂性 假設目標函數\(f: \mathbb{R}^d \rightarrow \mathbb{R}\)的導數\(\nabla f(w)\)是光滑的，存在二階導數，並且在其最優點處的導數為零，那麼牛頓法具有\(\mathcal{O}(e^{-2^t})\)的超線性（二階）收斂速率：

\[\lVert w^t - w^*\rVert \leqslant \mathcal{O}(e^{-2^t}) \]

然而，事物總有兩面性，相比一階方法，雖然牛頓法的收斂速度更快，但是存在下面兩個問題：

在每個時刻需要計算當前狀態的海森逆矩陣，計算量和儲存量都顯著增大；
海森矩陣不一定是正定的。
為了解決這個問題，人們提出了擬牛頓法。

2 擬牛頓法

2.1 演演算法描述

既然海森矩陣不一定正定，那就構造一個與海森矩陣相差不太遠的正定矩陣作為其替代品，這正是擬牛頓法^[5]的主要思想。此外，逆牛頓法可以迭代更新海森逆矩陣，而不是每次迭代都需要重新進行逆矩陣的計算。

記\(H_t = \nabla^2 f(w^t)\)，\(M^t = [\nabla^2 f(w^t)]^{-1}\)。在第\(t+1\)輪迭代，對第\(t\)輪迭代的海森矩陣\(H_t\)加上一個或者兩個秩為1的矩陣作為對\(t+1\)時刻的海森矩陣\(H^{t+1}\)的估計。例如：

\[H^{t+1} = H^t + aa^T + bb^T \]

當然，海森矩陣的更新需要滿足一定的條件，比如下面推導的擬牛頓條件。
對\(f(w)\)在\(w^t\)處的二階泰勒展開的左右兩邊計算梯度，可得

\[\nabla f(w) \approx \nabla f(w^t) + H^t(w-w^t) \]

令\(\delta^t_1 = w^{t+1}-w^t\)為模型的更新量，\(\delta^t_2=\nabla f(w^{t+1}) - \nabla f(w^t)\)為目標函數導數的更新量，則我們有：

\[(H^t)^{-1}\delta^t_2 \approx \delta^t_1 \]

在擬牛頓條件下，更新規則為

\[\begin{aligned} &H^0 = I\\ &H^{t+1}= H^t + \frac{\delta^t_2 (\delta^t_2)^{T}}{(\delta^t_2)^T\delta^t_1} - \frac{H^t\delta^t_1(H^t\delta^t_1)^T }{(\delta^t_1)^TH^t\delta^t_1} \\ &M^{t+1} = \left(I - \frac{\delta^t_1(\delta^t_2)^T}{(\delta^t_2)^T\delta^t_1}\right)M^t(I - \frac{\delta^t_2 (\delta^t_1)^T}{(\delta^t_2)^T\delta^t_1}) + \frac{\delta_1^t (\delta_1)^T}{(\delta^t_2)^T\delta^t_1} \end{aligned} \]

此時，所對應的演演算法被稱作BFGS演演算法，如下面虛擬碼所示：

在實際應用中基於\(M^t\)的擬牛頓法更實用，應為根據\(M^t\)計算下降方向的方法不需要對\(H^t\)矩陣求逆（解線性方程，其複雜度為\(\mathcal{O}(d^3)\)，非常耗時）。但基於\(H^t\)的擬牛頓法有比較好的理論性質，產生的迭代序列比較穩定，但如果有辦法快速求解線性方程組，我們也可以採用基於\(H^t\)的擬牛頓法。此外在某些場景下，比如有些帶約束的優化問題的演演算法設計，由於需要用到海森矩陣的近似，\(H^t\)的使用也很常見。

當使用其它的海森矩陣近似方式和約束條件時，我們還可以設計出其他擬牛頓法，比如DFP^[3][4]、Broyden^[5]、SR1^[6]等。

2.2 收斂性分析

可以證明，當初始點最優點足夠近時，擬牛頓法和牛頓法同樣具有二階收斂速率。

3 對偶方法

3.1 演演算法描述

我們到目前為止提到的各種優化演演算法都是直接求解原始優化問題。某些時候，如果把原始優化問題轉化成對偶優化問題，會更容易求解。比如，當原始問題的變數維度很高，但是約束條件個數不太多時，對偶問題的複雜度（對偶變數的維度對應於約束條件的個數）會遠小於原始問題的複雜度，因而更容易求解（這也是支援向量機高效的主要原因）。本節我們將介紹與此相關的對偶理論以及常見的對偶優化演演算法^[7]。

考慮下述原始優化問題\(P_0\)：

\[\underset{w}{\text{min}} f(w) \\ \begin{aligned} \quad\quad\quad\quad\quad\text{s.t.} \quad\quad\quad &g_i(w) \leqslant 0, \quad i=1,\cdots,m_1 \\ &h_j(w) = 0, \quad j=1,\cdots,m_2 \end{aligned} \]

包含\(m_1\)個不等式約束和\(m_2\)個等式約束。假設問題\(P_0\)的最優解為\(p^*\)。

首先，我們將有約束問題轉化為無約束優化問題，即構造\(P_0\)的拉格朗日函數，將約束條件轉化到目標函數中。考慮原始問題\(P_0\)中帶有的約束條件，例如，當\(g_1(w^0)\geqslant0\)時，\(w^0\)就不是一個可行解，這時我們需要把對應的目標函數之補充定義為正無窮。於是可將原始優化問題轉換為如下的無約束優化形式：

\[\underset{w}{\text{min}} f(w) + \infty \sum_{i=1}^{m_1} I_{g_i(w)>0} + \infty \sum_{j=1}^{m_2}I_{h_j(w)\neq 0} \]

進一步地，我們用線性函數\(\lambda x (\lambda > 0)\)和\(\nu x\)來近似指示函數\(I_{[x>0]}\)和\(I_{[x\neq 0]}\)，得到原始問題\(P_0\)所對應的拉格朗日函數：

\[L(w, \lambda, \nu) \triangleq f(w) + \sum_{i=1}^{m_1}\lambda_ig_i(w) + \sum_{j=1}^{m_2}\nu_jh_j(w) \]

其中，\(\lambda_i\geqslant 0\)，\(\nu_j\in \mathbb{R}\)稱為拉格朗日乘子。

接下來，我們討論拉格朗日函數和原目標函數取值的大小關係，引出拉格朗日對偶函數的定義。如果假設\(w^0\)滿足所有約束條件，那額拉格朗日函數中的第二項小於等於零，第三項等於零，即

\[L(w, \lambda, \nu) \leqslant f(w),\quad \text{對任意可行解}w \]

在可行域\(\mathcal{W}\)上對上述不等式取最小值，得到

\[\underset{w\in \mathcal{W}}{\text{inf}}L(w, \lambda, v) \leqslant \underset{w\in\mathcal{W}}{\text{min}} f(w) = p^* \]

我們將上式左端稱為拉格朗日對偶函數，簡記為

\[h(\lambda, \nu) \triangleq \underset{w\in \mathcal{W}}{\text{inf}}L(w, \lambda, v) \leqslant p^*, \quad 其中\lambda_i\geqslant0, \nu_j \in \mathbb{R} \]

由於拉格朗日對偶函數\(h(\lambda, \nu)\)是一族關於\((\lambda, \nu)\)仿射函數的逐點下确界，所以即使原始問題\(P_0\)不是凸的，拉格朗日對偶函數\(h(\lambda, \nu)\)也是凹的。

接下來我們定義對偶優化問題。既然拉格朗日對偶函數的取值是原始物體最優值的下界，在拉格朗日乘子\(\lambda\)，\(\nu\)的取值空間對函數\(h(\lambda, \nu)\)取最大值，將會得到原始問題最優解的最大下界。這個問題通常被定義為原始問題\(P_0\)的對偶問題，記為\(D_0\)，具體如下：

\[\underset{\lambda, \nu}{\text{max}} \space h(\lambda, \nu) \\ \begin{aligned} \quad\quad\quad\text{s.t.} \quad\quad\quad &\lambda_i \geqslant 0, \quad i=1,\cdots,m \\ \end{aligned} \]

對偶問題是一個對凹函數在可行域（凸集）內求解最大值的問題，記最優值為\(d^*\)。通過上面的討論，\(d^*<p^*\)恆成立，這個關係被稱為弱對偶條件；如果\(d^* = p^*\)，則稱強對偶條件成立。很多研究工作討論了強對偶條件成立的前提，比如Slater條件是強對偶的充分條件，KKT條件是強對偶的必要條件。

最後，可通過求解對偶問題來得到原始問題的解。在求解對偶問題的過程中，我們可以使用前面的各種一階或二階優化方法。假設對偶問題求得的解為\((\lambda^*, \nu^*)\)，將其代入拉格朗日函數中，對原始變數求拉格朗日函數的最小值，即:

\[\underset{w}{\text{min}}\space L(w, \lambda^*, \nu^*) \]

如果所求得的解是原始問題的一個可行解，那麼它就是原始問題的最優解。

如果一個優化問題需要通過求解其對偶問題來解決，我們首先推匯出它的對偶問題，並用各種一階或者二階方法來最大化對偶目標函數。比如，對偶座標上升法使用梯度上升法最大化對偶目標函數，如下面虛擬碼所示：

3.2 收斂性分析

對於帶有線性約束的凸優化問題，對偶座標上升法被證明至少具有線性收斂速率^[6]。

4 對確定性演演算法方法的小結

到目前為止，我們介紹了求解數值優化問題的一階、二階確定性演演算法及對偶方法。作為總結，我們來對比一下這些演演算法的效能。

給定一個演演算法，對於不同凸性和光滑性質的目標函數，其收斂速率有所不同。通常來講，更好的凸性和光滑性質會加速演演算法的收斂。同樣對於強凸且光滑的目標函數，不同演演算法的性質也不同，我們有如下結論：

一階方法具有線性收斂速率，二階方法具有超線性(二階)收斂速率。
Nesterov加速法將梯度下降法速率中關於條件數的階數進一步改進（但仍然是線性收斂速率）。
投影次梯度法和Frank-Wolfe方法都可以用於解決帶有約束的優化問題，兩種方法的收斂速率都為次線性，具體可以依據投影操作的難易程度來選擇。
一些擬牛頓法（比如BFGS演演算法）也可以和牛頓法一樣達到二階收斂速率。

確定性演演算法是數值優化的基石。機器學習實踐中更為常見也更實用的隨機優化演演算法就是建立在確定性演演算法之上的，我們將在後面的部落格中進行詳細介紹。

參考

[1] Polyak B T. Newton’s method and its use in optimization[J]. European Journal of Operational Research, 2007, 181(3): 1086-1096.
[2] Dennis, Jr J E, Moré J J. Quasi-Newton methods, motivation and theory[J]. SIAM review, 1977, 19(1): 46-89.
[3] Davidon W C. Variable metric method for minimization[J]. SIAM Journal on Optimization, 1991, 1(1): 1-17.
[4] Fletcher R, Powell M J D. A rapidly convergent descent method for minimization[J]. The computer journal, 1963, 6(2): 163-168.
[5] Broyden C G. A class of methods for solving nonlinear simultaneous equations[J]. Mathematics of computation, 1965, 19(92): 577-593.
[6] Luo Z Q, Tseng P. Error bounds and convergence analysis of feasible descent methods: a general approach[J]. Annals of Operations Research, 1993, 46(1): 157-178.
[7] Numerical optimization[M]. New York, NY: Springer New York, 1999.
[8] 劉浩洋，戶將等. 最佳化：建模、演演算法與理論[M]. 高教出版社, 2020.
[9] 劉鐵巖，陳薇等. 分散式機器學習：演演算法、理論與實踐[M]. 機械工業出版社, 2018.
[10] Stanford CME 323: Distributed Algorithms and Optimization (Lecture 7)