數值優化：經典隨機優化演演算法及其收斂性與複雜度分析

1 隨機優化演演算法概述

隨著巨量資料的出現，確定性優化演演算法的效率逐漸稱為瓶頸。為了說明這一點，我們來看一個用梯度下降法求解線性迴歸的例子。

給定訓練樣本\(D = \{(x_i, y_i)\}_{i=1}^n\)，線性迴歸的目標函數如下：

\[f(w) = \frac{1}{n}\sum_{i=1}^nf_i(w)= \frac{1}{n}\sum_{i=1}^n(w^T x_i - y_i)^2 \]

這裡\(w\in \mathbb{R}^d\)為模型引數。
梯度下降法的更新規則為：

\[w^{t+1} = w^t - \eta \nabla f(w^t) = w^t - \frac{2\eta}{n}\sum_{i=1}^nx_i\left((w^t)^Tx_i - y_i\right) \]

可見，梯度下降法中每次更新模型所需要的單位計算複雜度為\(\mathcal{O}(nd)\)。

對於更復雜的模型（比如神經網路）和更復雜的優化方法（比如二階方法），確定性優化方法的計算量會更大。那麼如何解決這個問題呢？

統計方法能給我們很大的幫助。雖然巨量資料的資料量和資料維度都很大，但我們可以通過對樣本和維度進行隨機取樣來得到對更新量的有效估計或者替代。相應地，從確定性優化演演算法出發，我們可以開發出各種隨機優化演演算法，如隨機梯度下降法^[1]、隨機座標下降法^[2]、隨機方差縮減梯度法^[3]、隨機（擬）牛頓法^[4]等。注意，對於隨機優化演演算法而言，收斂性分析與確定性演演算法不同，需要針對演演算法中的隨機取樣取期望。

下面就讓我們先介紹經典的隨機優化演演算法。

2 隨機梯度下降法

2.1 演演算法描述

隨機梯度下降法（SGD）^[1]對訓練資料做隨機取樣，其更新公式如下：

\[w^{t+1} = w^t - \eta^t \nabla f_{i^t}(w^t) \]

其中，\(i^t\)是第\(t\)輪隨機取樣的資料標號。具體演演算法如下列的虛擬碼所示：

我們知道，機器學習問題中的經驗損失函數定義為所有樣本資料對應的損失函數的平均值。而我們這裡用有放回隨機採用獲得的資料來計算梯度，是對用全部資料來計算梯度的一個無偏估計，即\(\mathbb{E}_{i^t} {\nabla_{i^t}f(w^t)} = \nabla f(w^t)\)，注意此處\(f(w^t)=\frac{1}{n}\sum_{i=1}^n\nabla f_i(w^t)\))。而由於每次更新只隨機採一個樣本，隨機梯度中的計算量大大減小。因此，隨機梯度可以作為梯度的自然替代，從而大大提高學習效率。

不過正如我們前面所說，優化演演算法的複雜度不僅包括單位計算複雜度，還包括迭代次數複雜度（與收斂率有關）。天下沒有免費的午餐，隨機梯度下降單位計算複雜度降低，那麼也相應地會付出迭代次數複雜度增大的代價。

考慮實際每次只採一個樣本比較極端，常用的形式是隨機梯度下降法的一個推廣：小批次（mini-batch）隨機梯度下降法。該方法可以看做是在隨機優化演演算法和確定性優化演演算法之間尋找某種折中，每次採一個較小的樣本集合\(\mathcal{I}^t\in \{1,2,...n\}\)（多於單樣本，少於全樣本），然後執行更新公式：

\[w^{t+1} = w^t - \eta^t\nabla f_{\mathcal{I^t}}(w^t) = w^t-\frac{\eta^t}{|\mathcal{I}^t|} \sum_{i\in\mathcal{I}^t} \nabla f_i (w^t) \]

我們後面介紹的各種隨機優化方法，都有相應的小批次版本。小批次取樣可以有效減小方差，從而提高收斂速率，具體我們下一篇部落格會討論。

2.2 收斂性和計算複雜度分析

與梯度下降法類似，對不同性質的目標函數，隨機梯度下降法具有不同的收斂速率。

\(L\)-Lipschitz連續凸函數收斂性：假設目標函數\(f: \mathbb{R}^d\rightarrow \mathbb{R}\)是凸函數，並且\(L\)-Lipschitz連續，\(w^* = \underset{\lVert w\rVert\leqslant D}{\text{arg min}}f(w)\)，當步長\(\eta^t = \sqrt{\frac{D^2}{L^2t}}\)時，對於給定的迭代步數\(T\)，隨機梯度下降法具有\(\mathcal{O}(\frac{1}{\sqrt{T}})\)的次線性收斂速率：

\[\mathbb{E}[ \frac{1}{T}\sum_{t=1}^T f(w^t) - f(w^*) ] \leqslant \frac{LD}{\sqrt{T}} \]

光滑強凸函數收斂性：假設目標函數\(f: \mathbb{R}^d\rightarrow \mathbb{R}\)是\(\alpha\)-強凸函數，並且\(\beta\)光滑，如果隨機梯度的二階矩有上界，即\(\mathbb{E}_{i^t}{\lVert\nabla f_{i^t}(w^t) \rVert^2\leqslant G^2}\)，當步長\(\eta^t = \frac{1}{\alpha t}\)時，對於給定的迭代步數\(T\)，隨機梯度下降法具有\(\mathcal{O}(\frac{1}{T})\)的次線性收斂速率：

\[\mathbb{E}[ f(w^T) - f(w^*) ] \leqslant \frac{2\beta G^2}{\alpha^2T} \]

通過與梯度下降法的收斂速率進行對比，我們可以發現隨機梯度下降法的收斂速率更慢。這主要是由於雖然隨機梯度是是全梯度的無偏估計，但這種估計存在一定的方差，會引入不確定性，導致最終演演算法的收斂速率下降。

雖然隨機梯度下降法的收斂速率慢於梯度下降法，但因為其每一輪的單位計算複雜度為\(\mathcal{O}(d)\)，而梯度下降法每一輪的單位計算複雜度為\(\mathcal{O}(nd)\)，所以當樣本量很大時，隨機梯度下降法的總計算複雜度\(\mathcal{O}(d(\frac{1}{\varepsilon}))\)比梯度下降法的總計算複雜度\(\mathcal{O}\left(ndQ\text{log}(\frac{1}{\epsilon})\right)\)要低。

3 隨機座標下降法

3.1 演演算法描述

除了對樣本進行隨機取樣外，還可以對模型的維度進行取樣，相應的演演算法稱為隨機座標下降法^[2]，其更新公式如下：

\[w^{t+1}_{j^t} = w^t_{j^t} - \eta^t \nabla_{j^t}f(w^t) \]

其中\(j^t\)表示第\(t\)輪迭代中隨機採的維度標號。\(\nabla_{j^t}f(w^t)\)是損失函數對於模型\(w^t\)中的第\(j^t\)個維度的偏導數。

隨機座標下降法虛擬碼如下：

一方面，如果取樣方法是又放回取樣，那麼可以得到\(\mathbb{E}_{j^t}\nabla_{j^t}f(w^t) = \frac{1}{d} \nabla f(w^t)\)（為了不引入新的記號，設\(\nabla_{j^t}f(w^t)\)是\(d\)維向量，其第\(j^t\)個維度是\(\nabla_{j^t}f(w^t)\)，其它維度是0）。也就是說，在期望意義上，隨機座標梯度與梯度方向是一致的。另一方面，對於線性模型，計算一個維度的梯度所需要的計算量只有整個梯度向量的\(\frac{1}{d}\)。因此，隨機座標梯度可以作為原始梯度的高效替代品（尤其是在引數維度較高時）。

隨機座標下降法的一個推廣版本是隨機塊座標下降法，也就是將引數維度分為\(J\)個塊，每次採一個塊\(J^t\in \{1,2,...J\}\)，然後執行更新公式：

\[w^{t+1}_j = w^t_j - \eta^t\nabla_{J^t} f(w^t) \]

我們後面介紹的各種隨機優化方法，都有相應的小批次版本。小批次取樣可以有效減小方差，從而提高收斂速率，具體我們下一篇部落格會討論。

3.2 收斂性和計算複雜度分析

對於不同性質的目標函數，隨機座標下降法也有不同的收斂速率。

由於模型每次只更新一個維度，為了對隨機座標下降法進行分析，我們先來刻畫偏導數的連續性質。

偏導數的連續性質 如果對任意模型\(w\in \mathbb{R}^d\)，對於維度\(j\)存在常數\(\beta_j\)，使得\(\forall \delta \in \mathbb{R}\)有下面不等式成立：

\[|\nabla_j f(w + \delta e^j) - \nabla_j f(w)| \leqslant \beta_j |\delta| \]

則稱目標函數\(f\)對於維度\(j\)具有\(\beta_j\)-Lipschitz連續的偏導數。

如果\(f\)對於每個維度的偏導數都是Lipchitz連續的，我們記\(\beta_{\text{max}} = \underset{j=1,2,\cdots, d}{\text{max}}\beta_j\)。可以驗證，如果目標函數是\(\beta\)光滑的，那麼\(\beta_{\text{max}}\leqslant \beta \leqslant \sqrt{d}\beta_j, \forall j=1,2,\cdots, d\)。

偏導數滿足\(\beta_j\)-Lipschitz連續的凸函數收斂性分析: 假設目標函數\(\mathbb{R}^d\rightarrow \mathbb{R}\)是凸函數，並且具有具有\(\beta_j\)-Lipschitz連續的偏導數，記\(w^* = \underset{\lVert w\rVert \leqslant D}{\text{arg min}} f(w)\)，當步長\(\eta = \frac{1}{\beta_{\text{max}}}\)時，對於給定的迭代步數\(T\)，隨機座標下降法具有\(\mathcal{O}(\frac{d\beta_{\text{max}}}{T})\)的次線性收斂速率：

\[\mathbb{E}[f(w^T) - f(w^*)] \leqslant \frac{2d\beta_{\text{max}}D^2}{T} \]

偏導數滿足\(\beta_j\)-Lipschitz連續的強凸函數收斂性分析: 假設目標函數\(\mathbb{R}^d\rightarrow \mathbb{R}\)是強凸函數，並且具有具有\(\beta_j\)-Lipschitz連續的偏導數，當步長\(\eta = \frac{1}{\beta_{\text{max}}}\)時，對於給定的迭代步數\(T\)，隨機座標下降法具有如下的線性收斂速率：

\[\mathbb{E}[f(w^T) - f(w^*)] \leqslant （1-\frac{\alpha}{d\beta_{\text{max}}}）^T (f(w^0) - f(w^*)) \]

對比梯度下降法的收斂速率，我們發下隨機梯度下降法的收斂速率關於迭代次數\(T\)的階數與梯度下降法的是一致的。從這個意義上講，隨機座標下降法的理論性質優於隨機梯度下降法。

在計算複雜度上，雖然隨機座標下降法線上性迴歸問題中的更新公式為

\[\begin{aligned} w^{t+1}_j & = w^t_j - \eta^t\nabla_j f(w^t)\\ & = w^t_j - \frac{2\eta_t}{n}\sum_{i=1}^nx_{i,j}((w^t)^Tx_i - y_i) \end{aligned} \]

雖然隨機座標下降法每輪迭代也需要像梯度下降法一樣計算\(n\)個內積\(\{(w^t)^Tx_i\}_{i=1}^n\)，但每個內積的計算量不再是\(\mathcal{O}(d)\) 。因為\(w^t\)每次只更新一維，可通過引入輔助變數的形式將\(\mathcal{O}(d)\)的計算量降為\(\mathcal{O}(1)\)，最終得到\(\mathcal{O}(n)\)的單位計算複雜度。這種偏導數的計算量小於梯度的計算量的情形一般被稱為「可分離」情形。可以證明，對於線性模型，常用損失函數都是可分離的。

至於隨機塊座標下降法的收斂性則與隨機座標下降法基本相同，只是其中的維度數目\(d\)會被塊的個數\(J\)所取代。

4 隨機擬牛頓法

4.1 演演算法描述

隨機擬牛頓法^[4]的思想與一階隨機演演算法類似，隨機採一個樣本或者一個小批次樣本來計算梯度，然後更新海森逆矩陣。小批次隨機擬牛頓法的更新公式如下：

\[w^{t+1} = w^t - \eta^t M^t \left(\frac{1}{|\mathcal{I}^t|} \sum_{i\in \mathcal{I}^t}\nabla f_i(w^t) \right) \]

其中\(\mathcal{I}^t\)是所採的小批次資料子集，\(M^t\)為\(\mathcal{I}\)上目標函數的海森逆矩陣。

類似於上一章介紹的擬牛頓法，雖然直接計算海森逆矩陣\(M^t\)的複雜度很高，但是利用歷史資訊迭代更新上一輪海森逆矩陣\(M^{t-1}\)可以得到對\(M^t\)的良好逼近，其計算複雜度要低很多。具體而言，首先在演演算法執行過程中記錄下表徵模型變化量和梯度變化量的修正組合\((\delta^s_1, \delta^s_2)\)，也就是

\[\begin{aligned} &\delta^s_1 = w^s - w^{s-1} \\ &\delta^s_2 = \frac{1}{|\mathcal{I^t}|}\sum_{i\in \mathcal{I^t}}\nabla^2f_i(w^s)(w^s - w^{s-1}) \end{aligned} \]

然後依據第\(t\)輪迭代之前的多個修正組合按照如下公式迭代更新海森逆矩陣：

\[M \leftarrow(I - \rho^s\delta^s_1(\delta^s_2)^T)M(I - \rho^s\delta^s_2(\delta^s_1)^T) + \rho^s \delta^s_1(\delta^s_1)^T \]

其中\(\rho^s = \frac{1}{(\delta^s_2)^T\delta^s_1}\)。

隨機擬牛頓法及海森逆矩陣更新的具體演演算法的虛擬碼如下：

4.2 收斂性和計算複雜度分析

在以下四條假設之下，我們可以證明隨機擬牛頓法具有與隨機梯度下降法相同的收斂速率：

目標函數\(f(w)\)是二階連續可微的；
存在\(\lambda_2 > \lambda_1 > 0\)，使得對於任意\(w\in \mathbb{R}^d\)，\(\lambda_1 I \prec \nabla^2 f(w) \prec \lambda_2I\)；
存在\(\mu_2 > \mu_1 > 0\)，使得對於任意\(w\in \mathbb{R}^d\)，\(\mu_1 I \prec (\nabla^2 f(w))^{-1} \prec \mu_2I\)；
存在\(G>0\)，使得對於任意\(w\in \mathbb{R}^d\)，\(\mathbb{E}[\lVert \nabla f_{i}(w) \rVert^2] \leqslant G^2\)。

假設上述條件成立，當步長\(\eta^t = \frac{a}{t}\)並且\(a > \frac{1}{2\mu_1\lambda_1}\)時，對於給定的迭代步數\(T\)，隨機擬牛頓法具有\(\mathcal{O}(\frac{1}{T})\)的次線性收斂速率：

\[\mathbb{E}[f(w^T)-f(w^*)] \leqslant \frac{Q(a)}{T} \]

其中\(Q(a) = \max \{ \frac{\lambda_2\mu_2^2a^2G^2}{2(2\mu_1\lambda_1a-1)}, f(w^1) - f(w^*) \}\)。

下面我們以邏輯迴歸為例對比一下隨機擬牛頓法和隨機梯度下降法的複雜度。隨機擬牛頓法每一輪的計算量（浮點運算次數）為\(2bd+4Sd+\frac{3b_Hd}{L}\)，其中\(b\)是隨機梯度下降小批次資料子集的大小，\(b_h\)是計算修正項的小批次資料子集的大小，\(S\)是記憶體引數，\(L\)是修正步驟的輪數。由此，可以得到隨機擬牛頓法和隨機梯度下降法每一輪計算的複雜度之比為

\[1 + \frac{2M}{b} + \frac{3b_H}{2bL} \]

於是，我們可以通過設定合適的引數使得隨機擬牛頓法的複雜度與隨機梯度下降法同階。
依據上面的討論，二階隨機演演算法在收斂速率和複雜度上都與一階隨機演演算法差不多，不像確定性演演算法那樣收斂速率有顯著的提高。原因是，對於更精細的二階演演算法，隨機取樣的方差會影響收斂精度。如何提高二階隨機優化演演算法的效率，仍然是未解決的問題。

5 隨機對偶座標上升法

5.1 演演算法描述

考慮線性模型，其對應的正則化風險經驗最小化（R-ERM）過程如下：

\[\underset{w\in\mathbb{R}^d}{\text{min}} f(w) = \frac{1}{n}\sum_{i=1}^n\phi_i(w^Tx_i) + \frac{\lambda}{2}\lVert w\rVert^2 \]

其中\(\phi_i(w^Tx_i)=\mathcal{l}(w; x_i, y_i)\)為線性模型\(w\in\mathbb{R}^d\)在樣本\((x_i, y_i)\)上的損失函數。

上述原始問題的對偶問題為：

\[\underset{\alpha \in\mathbb{R}^n}{\text{max}} D(\alpha) = \frac{1}{n}\sum_{i=1}^n - \phi_i^*(-\alpha_i) - \frac{\lambda n}{2}\left\lVert \frac{1}{\lambda n}\sum_{i=1}^n\alpha_ix_i \right\rVert^2 \]

其中\(\phi_i^*(u)=\underset{z}{\text{max}}(zu - \phi_i(z))\)。

上述原始問題和對偶問題相比我們在部落格《數值優化：經典二階確定性演演算法與對偶方法》中介紹的問題更加特殊，利用Fenchel對偶定理，如果定義 \(w(\alpha)= \frac{1}{\lambda n}\sum_{i=1}^n\alpha_ix_i\)，並且原始目標函數及其對偶函數分別存在最優解\(w^*\)和\(\alpha^*\)。

\[w(\alpha^*) = w^*, f(w^*) = D(\alpha^*) \]

優於對偶問題是線性可分的，隨機座標上升法比確定性的梯度上升法更能有效地對其進行優化，我們稱之為隨機對偶座標上升法（SDCA）^[5]。其主要計算步驟為：

隨機採一個樣本\(i\)，計算

\[\Delta \alpha_i = \underset{z}{\text{argmax}} \left\{ -\phi_i^*(-\alpha_i^t + z) - \frac{\lambda n}{2} \lVert w^t + \frac{1}{\lambda n} zx_i \rVert^2 \right\} \]

更新對偶變數，\(\alpha^{t+1} = \alpha^t + \Delta \alpha_ie_i\)。
更新原始變數，\(w^{t+1} = w^t + \frac{1}{\lambda n} \Delta \alpha_i x_i\)。

隨機對偶座標上升法虛擬碼如下：

請注意，對於隨機座標上升法，損失函數的光滑性質對收斂速率有顯著的影響，因為如果損失函數\(\phi_i(\alpha)\)是\(\beta\)-光滑的，那麼其對偶函數\(\phi^*_i(u)\)是\(\frac{1}{\beta}\)-強凹的，於是對偶問題的凹性得到了加強。

5.2 收斂性和計算複雜度分析

如果損失函數是凸函數且是\(L\)-Lipschitz連續的，隨機對偶座標上升法具有次線性收斂速率\(\mathcal{O}(\frac{L^2}{\lambda n }+\frac{L^2}{\lambda t})\)；如果損失函數進一步是\(\beta\)-光滑的，隨機對偶座標上升法具有線性收斂速率\(\mathcal{O}(e^{-\frac{\lambda t}{\beta + \lambda n}})\)。

隨機對偶梯度上升法的確定性版本的收斂速率與上面所述相同，但是由於線性模型的正則化損失函數是線性可分的，隨機對偶座標上升法中每次迭代的計算量從\(\mathcal{O}(d)\)減小為\(\mathcal{O}(1)\)，從而提高了演演算法效率。

6 隨機優化演演算法總結

前面已介紹的幾種隨機優化演演算法的收斂速率、單位計算複雜度和總計算複雜度總結如下表所示：

上表中，\(T\)為迭代次數，\(\beta\)和\(\beta_{\text{max}}\)為光滑係數和各個維度對應的最大光滑係數，\(L\)和\(L_{\text{max}}\)為Lipschitz係數和各個維度對應的最大Lipschitz係數，\(Q\)為條件數，\(n\)為資料量，\(d\)為資料維度，\(b\)和\(b_H\)為隨機演演算法中求海森逆矩陣的小批次資料集的大小，\(\lambda\)為拉格朗日系數。

從表中可以得出以下幾點結論：

當資料量較大時，隨機梯度下降法比梯度下降法更高效。
如果目標函數是可分離的，隨機座標下降法比梯度下降法更高效。
如果目標函數是可分離的，並且資料維度較高，隨機座標下降法比隨機梯度下降法更高效。
隨機擬牛頓法的效率與隨機梯度下降法的效率相同。

參考

[1] Robbins H, Monro S. A stochastic approximation method[J]. The annals of mathematical statistics, 1951: 400-407.
[2] Nesterov Y. Efficiency of coordinate descent methods on huge-scale optimization problems[J]. SIAM Journal on Optimization, 2012, 22(2): 341-362.
[3] Johnson R, Zhang T. Accelerating stochastic gradient descent using predictive variance reduction[J]. Advances in neural information processing systems, 2013, 26.
[4] Byrd R H, Hansen S L, Nocedal J, et al. A stochastic quasi-Newton method for large-scale optimization[J]. SIAM Journal on Optimization, 2016, 26(2): 1008-1031.
[5] Shalev-Shwartz S, Zhang T. Stochastic dual coordinate ascent methods for regularized loss minimization[J]. Journal of Machine Learning Research, 2013, 14(Feb):567-599.
[6] 劉浩洋，戶將等. 最佳化：建模、演演算法與理論[M]. 高教出版社, 2020.
[7] 劉鐵巖，陳薇等. 分散式機器學習：演演算法、理論與實踐[M]. 機械工業出版社, 2018.
[8] Stanford CME 323: Distributed Algorithms and Optimization (Lecture 7)