数值优化：经典二阶确定性算法与对偶方法

我们在上一篇博客中主要介绍了单机数值优化中一些经典的一阶确定性算法，本篇文章我们将会介绍二阶确定性算法和对偶方法。

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)\)是光滑的，存在二阶导数，并且在其最优点\(x^*\)处满足\(\nabla f\left(x^{*}\right)=0, \nabla^{2} f\left(x^{*}\right) \succ 0\)且初始点\(x^0\)与最优点\(x^*\)足够近，那么牛顿法具有\(\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^*恒成立，这个关系被称为弱对偶条件；如果\(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)