数值优化：经典随机优化算法及其收敛性与复杂度分析

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)