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

1 随机优化算法概述

随着大数据的出现，确定性优化算法的效率逐渐称为瓶颈。为了说明这一点，我们来看一个用梯度下降法求解线性回归的例子。

给定训练样本D={(xi,yi)}ni=1D={(xi,yi)}ni=1，线性回归的目标函数如下：

f(w)=1nn∑i=1fi(w)=1nn∑i=1(wTxi−yi)2

f(w)=1n∑i=1nfi(w)=1n∑i=1n(wTxi−yi)2

这里w∈Rd为模型参数。
梯度下降法的更新规则为：

wt+1=wt−η∇f(wt)=wt−2ηnn∑i=1xi((wt)Txi−yi)

可见，梯度下降法中每次更新模型所需要的单位计算复杂度为O(nd)。

对于更复杂的模型（比如神经网络）和更复杂的优化方法（比如二阶方法），确定性优化方法的计算量会更大。那么如何解决这个问题呢？

统计方法能给我们很大的帮助。虽然大数据的数据量和数据维度都很大，但我们可以通过对样本和维度进行随机采样来得到对更新量的有效估计或者替代。相应地，从确定性优化算法出发，我们可以开发出各种随机优化算法，如随机梯度下降法^[1]、随机坐标下降法^[2]、随机方差缩减梯度法^[3]、随机（拟）牛顿法^[4]等。注意，对于随机优化算法而言，收敛性分析与确定性算法不同，需要针对算法中的随机采样取期望。

下面就让我们先介绍经典的随机优化算法。

2 随机梯度下降法

2.1 算法描述

随机梯度下降法（SGD）^[1]对训练数据做随机采样，其更新公式如下：

wt+1=wt−ηt∇fit(wt)

其中，it是第t轮随机采样的数据标号。具体算法如下列的伪代码所示：

我们知道，机器学习问题中的经验损失函数定义为所有样本数据对应的损失函数的平均值。而我们这里用有放回随机采用获得的数据来计算梯度，是对用全部数据来计算梯度的一个无偏估计，即Eit∇itf(wt)=∇f(wt)，注意此处f(wt)=1n∑ni=1∇fi(wt))。而由于每次更新只随机采一个样本，随机梯度中的计算量大大减小。因此，随机梯度可以作为梯度的自然替代，从而大大提高学习效率。

不过正如我们前面所说，优化算法的复杂度不仅包括单位计算复杂度，还包括迭代次数复杂度（与收敛率有关）。天下没有免费的午餐，随机梯度下降单位计算复杂度降低，那么也相应地会付出迭代次数复杂度增大的代价。

考虑实际每次只采一个样本比较极端，常用的形式是随机梯度下降法的一个推广：小批量（mini-batch）随机梯度下降法。该方法可以看做是在随机优化算法和确定性优化算法之间寻找某种折中，每次采一个较小的样本集合It∈{1,2,...n}（多于单样本，少于全样本），然后执行更新公式：

wt+1=wt−ηt∇fIt(wt)=wt−ηt|It|∑i∈It∇fi(wt)

我们后面介绍的各种随机优化方法，都有相应的小批量版本。小批量采样可以有效减小方差，从而提高收敛速率，具体我们下一篇博客会讨论。

2.2 收敛性和计算复杂度分析

与梯度下降法类似，对不同性质的目标函数，随机梯度下降法具有不同的收敛速率。

L-Lipschitz连续凸函数收敛性：假设目标函数f:Rd→R是凸函数，并且L-Lipschitz连续，w∗=arg min‖w‖⩽Df(w)，当步长ηt=√D2L2t时，对于给定的迭代步数T，随机梯度下降法具有O(1√T)的次线性收敛速率：

E[1TT∑t=1f(wt)−f(w∗)]⩽LD√T

光滑强凸函数收敛性：假设目标函数f:Rd→R是α-强凸函数，并且β光滑，如果随机梯度的二阶矩有上界，即Eit‖∇fit(wt)‖2⩽G2，当步长ηt=1αt时，对于给定的迭代步数T，随机梯度下降法具有O(1T)的次线性收敛速率：

E[f(wT)−f(w∗)]⩽2βG2α2T

通过与梯度下降法的收敛速率进行对比，我们可以发现随机梯度下降法的收敛速率更慢。这主要是由于虽然随机梯度是是全梯度的无偏估计，但这种估计存在一定的方差，会引入不确定性，导致最终算法的收敛速率下降。

虽然随机梯度下降法的收敛速率慢于梯度下降法，但因为其每一轮的单位计算复杂度为O(d)，而梯度下降法每一轮的单位计算复杂度为O(nd)，所以当样本量很大时，随机梯度下降法的总计算复杂度O(d(1ε))比梯度下降法的总计算复杂度O(ndQlog(1ϵ))要低。

3 随机坐标下降法

3.1 算法描述

除了对样本进行随机采样外，还可以对模型的维度进行采样，相应的算法称为随机坐标下降法^[2]，其更新公式如下：

wt+1jt=wtjt−ηt∇jtf(wt)

其中jt表示第t轮迭代中随机采的维度标号。∇jtf(wt)是损失函数对于模型wt中的第jt个维度的偏导数。

随机坐标下降法伪代码如下：

一方面，如果采样方法是又放回采样，那么可以得到Ejt∇jtf(wt)=1d∇f(wt)（为了不引入新的记号，设∇jtf(wt)是d维向量，其第jt个维度是∇jtf(wt)，其它维度是0）。也就是说，在期望意义上，随机坐标梯度与梯度方向是一致的。另一方面，对于线性模型，计算一个维度的梯度所需要的计算量只有整个梯度向量的1d。因此，随机坐标梯度可以作为原始梯度的高效替代品（尤其是在参数维度较高时）。

随机坐标下降法的一个推广版本是随机块坐标下降法，也就是将参数维度分为J个块，每次采一个块Jt∈{1,2,...J}，然后执行更新公式：

wt+1j=wtj−ηt∇Jtf(wt)

我们后面介绍的各种随机优化方法，都有相应的小批量版本。小批量采样可以有效减小方差，从而提高收敛速率，具体我们下一篇博客会讨论。

3.2 收敛性和计算复杂度分析

对于不同性质的目标函数，随机坐标下降法也有不同的收敛速率。

由于模型每次只更新一个维度，为了对随机坐标下降法进行分析，我们先来刻画偏导数的连续性质。

偏导数的连续性质 如果对任意模型w∈Rd，对于维度j存在常数βj，使得∀δ∈R有下面不等式成立：

|∇jf(w+δej)−∇jf(w)|⩽βj|δ|

则称目标函数f对于维度j具有βj-Lipschitz连续的偏导数。

如果f对于每个维度的偏导数都是Lipchitz连续的，我们记βmax=maxj=1,2,⋯,dβj。可以验证，如果目标函数是β光滑的，那么βmax⩽β⩽√dβj,∀j=1,2,⋯,d。

偏导数满足βj-Lipschitz连续的凸函数收敛性分析: 假设目标函数Rd→R是凸函数，并且具有具有βj-Lipschitz连续的偏导数，记w∗=arg min‖w‖⩽Df(w)，当步长η=1βmax时，对于给定的迭代步数T，随机坐标下降法具有O(dβmaxT)的次线性收敛速率：

E[f(wT)−f(w∗)]⩽2dβmaxD2T

偏导数满足βj-Lipschitz连续的强凸函数收敛性分析: 假设目标函数Rd→R是强凸函数，并且具有具有βj-Lipschitz连续的偏导数，当步长η=1βmax时，对于给定的迭代步数T，随机坐标下降法具有如下的线性收敛速率：

E[f(wT)−f(w∗)]⩽（1−αdβmax）T(f(w0)−f(w∗))

对比梯度下降法的收敛速率，我们发下随机梯度下降法的收敛速率关于迭代次数T的阶数与梯度下降法的是一致的。从这个意义上讲，随机坐标下降法的理论性质优于随机梯度下降法。

在计算复杂度上，虽然随机坐标下降法在线性回归问题中的更新公式为

wt+1j=wtj−ηt∇jf(wt)=wtj−2ηtnn∑i=1xi,j((wt)Txi−yi)

虽然随机坐标下降法每轮迭代也需要像梯度下降法一样计算n个内积{(wt)Txi}ni=1，但每个内积的计算量不再是O(d) 。因为wt每次只更新一维，可通过引入辅助变量的形式将O(d)的计算量降为O(1)，最终得到O(n)的单位计算复杂度。这种偏导数的计算量小于梯度的计算量的情形一般被称为“可分离”情形。可以证明，对于线性模型，常用损失函数都是可分离的。

至于随机块坐标下降法的收敛性则与随机坐标下降法基本相同，只是其中的维度数目d会被块的个数J所取代。

4 随机拟牛顿法

4.1 算法描述

随机拟牛顿法^[4]的思想与一阶随机算法类似，随机采一个样本或者一个小批量样本来计算梯度，然后更新海森逆矩阵。小批量随机拟牛顿法的更新公式如下：

wt+1=wt−ηtMt(1|It|∑i∈It∇fi(wt))

其中It是所采的小批量数据子集，Mt为I上目标函数的海森逆矩阵。

类似于上一章介绍的拟牛顿法，虽然直接计算海森逆矩阵Mt的复杂度很高，但是利用历史信息迭代更新上一轮海森逆矩阵Mt−1可以得到对Mt的良好逼近，其计算复杂度要低很多。具体而言，首先在算法运行过程中记录下表征模型变化量和梯度变化量的修正组合(δs1,δs2)，也就是

δs1=ws−ws−1δs2=1|It|∑i∈It∇2fi(ws)(ws−ws−1)

然后依据第t轮迭代之前的多个修正组合按照如下公式迭代更新海森逆矩阵：

M←(I−ρsδs1(δs2)T)M(I−ρsδs2(δs1)T)+ρsδs1(δs1)T

其中ρs=1(δs2)Tδs1。

随机拟牛顿法及海森逆矩阵更新的具体算法的伪代码如下：

4.2 收敛性和计算复杂度分析

在以下四条假设之下，我们可以证明随机拟牛顿法具有与随机梯度下降法相同的收敛速率：

• 目标函数f(w)是二阶连续可微的；
• 存在λ2>λ1>0，使得对于任意w∈Rd，λ1I≺∇2f(w)≺λ2I；
• 存在μ2>μ1>0，使得对于任意w∈Rd，μ1I≺(∇2f(w))−1≺μ2I；
• 存在G>0，使得对于任意w∈Rd，E[‖∇fi(w)‖2]⩽G2。

假设上述条件成立，当步长ηt=at并且a>12μ1λ1时，对于给定的迭代步数T，随机拟牛顿法具有O(1T)的次线性收敛速率：

E[f(wT)−f(w∗)]⩽Q(a)T

其中Q(a)=max{λ2μ22a2G22(2μ1λ1a−1),f(w1)−f(w∗)}。

下面我们以逻辑回归为例对比一下随机拟牛顿法和随机梯度下降法的复杂度。随机拟牛顿法每一轮的计算量（浮点运算次数）为2bd+4Sd+3bHdL，其中b是随机梯度下降小批量数据子集的大小，bh是计算修正项的小批量数据子集的大小，S是内存参数，L是修正步骤的轮数。由此，可以得到随机拟牛顿法和随机梯度下降法每一轮计算的复杂度之比为

1+2Mb+3bH2bL

于是，我们可以通过设置合适的参数使得随机拟牛顿法的复杂度与随机梯度下降法同阶。
依据上面的讨论，二阶随机算法在收敛速率和复杂度上都与一阶随机算法差不多，不像确定性算法那样收敛速率有显著的提高。原因是，对于更精细的二阶算法，随机采样的方差会影响收敛精度。如何提高二阶随机优化算法的效率，仍然是未解决的问题。

5 随机对偶坐标上升法

5.1 算法描述

考虑线性模型，其对应的正则化风险经验最小化（R-ERM）过程如下：

minw∈Rdf(w)=1nn∑i=1ϕi(wTxi)+λ2‖w‖2

其中ϕi(wTxi)=l(w;xi,yi)为线性模型w∈Rd在样本(xi,yi)上的损失函数。

上述原始问题的对偶问题为：

maxα∈RnD(α)=1nn∑i=1−ϕ∗i(−αi)−λn2‖1λnn∑i=1αixi‖2

其中ϕ∗i(u)=maxz(zu−ϕi(z))。

上述原始问题和对偶问题相比我们在博客《数值优化：经典二阶确定性算法与对偶方法》中介绍的问题更加特殊，利用Fenchel对偶定理，如果定义 w(α)=1λn∑ni=1αixi，并且原始目标函数及其对偶函数分别存在最优解w∗和α∗。

w(α∗)=w∗,f(w∗)=D(α∗)

优于对偶问题是线性可分的，随机坐标上升法比确定性的梯度上升法更能有效地对其进行优化，我们称之为随机对偶坐标上升法（SDCA）^[5]。其主要计算步骤为：

• 随机采一个样本i，计算

Δαi=argmaxz{−ϕ∗i(−αti+z)−λn2‖wt+1λnzxi‖2}

• 更新对偶变量，αt+1=αt+Δαiei。
• 更新原始变量，wt+1=wt+1λnΔαixi。

随机对偶坐标上升法伪代码如下：

请注意，对于随机坐标上升法，损失函数的光滑性质对收敛速率有显著的影响，因为如果损失函数ϕi(α)是β-光滑的，那么其对偶函数ϕ∗i(u)是1β-强凹的，于是对偶问题的凹性得到了加强。

5.2 收敛性和计算复杂度分析

如果损失函数是凸函数且是L-Lipschitz连续的，随机对偶坐标上升法具有次线性收敛速率O(L2λn+L2λt)；如果损失函数进一步是β-光滑的，随机对偶坐标上升法具有线性收敛速率O(e−λtβ+λn)。

随机对偶梯度上升法的确定性版本的收敛速率与上面所述相同，但是由于线性模型的正则化损失函数是线性可分的，随机对偶坐标上升法中每次迭代的计算量从O(d)减小为O(1)，从而提高了算法效率。

6 随机优化算法总结

前面已介绍的几种随机优化算法的收敛速率、单位计算复杂度和总计算复杂度总结如下表所示：

上表中，T为迭代次数，β和βmax为光滑系数和各个维度对应的最大光滑系数，L和Lmax为Lipschitz系数和各个维度对应的最大Lipschitz系数，Q为条件数，n为数据量，d为数据维度，b和bH为随机算法中求海森逆矩阵的小批量数据集的大小，λ为拉格朗日系数。

从表中可以得出以下几点结论：

• 当数据量较大时，随机梯度下降法比梯度下降法更高效。
• 如果目标函数是可分离的，随机坐标下降法比梯度下降法更高效。
• 如果目标函数是可分离的，并且数据维度较高，随机坐标下降法比随机梯度下降法更高效。
• 随机拟牛顿法的效率与随机梯度下降法的效率相同。

参考

• [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)

• 1 随机优化算法概述
• 2 随机梯度下降法
•     2.1 算法描述
•     2.2 收敛性和计算复杂度分析
• 3 随机坐标下降法
•     3.1 算法描述
•     3.2 收敛性和计算复杂度分析
• 4 随机拟牛顿法
•     4.1 算法描述
•     4.2 收敛性和计算复杂度分析
• 5 随机对偶坐标上升法
•     5.1 算法描述
•     5.2 收敛性和计算复杂度分析
• 6 随机优化算法总结
• 参考

__EOF__