文章目录

算法结构

梯度类算法有很多，本文主要学习最常见的3个算法：最速下降法、牛顿法和拟牛顿法。算法名称虽多，但是他们的算法结构都是一样的，可以描述为

（1）选定初始点 x 0 \pmb {x}_0 x0​。

（2）迭代公式为 x k + 1 = x k + α k d k \pmb{x}_{k+1}=\pmb{x}_k+\alpha_k\pmb{d}_k xk+1​=xk​+αk​dk​，其中 α k \alpha_k αk​和 d k \pmb{d}_k dk​分别为第 k k k次的迭代步长和迭代方向。

人们常说，选择大于努力。对优化算法来说，也是类似的：如果 d k \pmb{d}_k dk​设计的不好，可能无论 α k \alpha_k αk​怎么努力，都得不到最优解。所以，我们应该重点关注 d k \pmb{d}_k dk​的构造方法。

最速下降法

最朴素的想法，就是如果针对任意初始点，都能直接算出让函数下降最快的方向，就好了。下面我们看看怎么能够做到。

假设 x k = [ x k 1 , x k 2 , . . , x k n ] \pmb{x}_k=[x_{k1},x_{k2},..,x_{kn}] xk​=[xk1​,xk2​,..,xkn​]，给定方向 s \pmb{s} s上的一个增量为

Δ s = [ Δ x 1 , Δ x 2 , . . . , Δ x n ] \Delta \pmb{s}=[\Delta x_1, \Delta x_2, ..., \Delta x_n] Δs=[Δx1​,Δx2​,...,Δxn​]

定义在点 x k \pmb{x}_k xk​处沿方向 s \pmb{s} s的变化率为

∂ f ( x ) ∂ s ∣ x k = lim ⁡ Δ s → 0 f ( x k + Δ s ) − f ( x k ) ∥ Δ s ∥ \frac{\partial f(\pmb{x})}{\partial \pmb{s}} \vert _{\pmb{x}_k}= \lim\limits_{\Delta \pmb{s}\to0}\frac{f(\pmb{x}_k+\Delta \pmb{s})-f(\pmb{x}_k)}{\Vert \Delta \pmb{s} \Vert} ∂s∂f(x)​∣xk​​=Δs→0lim​∥Δs∥f(xk​+Δs)−f(xk​)​

当 Δ s \Delta \pmb{s} Δs足够小时，上式右侧的分子变为全微分形式

f ( x k + Δ s ) − f ( x k ) = ∂ f ∂ x k 1 ∣ x k Δ x 1 + ∂ f ∂ x k 2 ∣ x k Δ x 2 + ⋅ ⋅ ⋅ + ∂ f ∂ x k n ∣ x k Δ x n f(\pmb{x}_k+\Delta \pmb{s})-f(\pmb{x}_k)=\frac{\partial f}{\partial x_{k1}}\vert _{\pmb{x}_k}\Delta x_{1}+\frac{\partial f}{\partial x_{k2}}\vert _{\pmb{x}_k}\Delta x_2+···+\frac{\partial f}{\partial x_{kn}}\vert _{\pmb{x}_k}\Delta x_n f(xk​+Δs)−f(xk​)=∂xk1​∂f​∣xk​​Δx1​+∂xk2​∂f​∣xk​​Δx2​+⋅⋅⋅+∂xkn​∂f​∣xk​​Δxn​

将上式带入上上式

∂ f ( x ) ∂ s ∣ x k = lim ⁡ Δ s → 0 ( ∂ f ∂ x k 1 ∣ x k Δ x 1 Δ s + ∂ f ∂ x k 2 ∣ x k Δ x 2 Δ s + ⋅ ⋅ ⋅ + ∂ f ∂ x k n ∣ x k Δ x n Δ s ) = ∂ f ∂ x 1 ∣ x k cos ⁡ α 1 + ∂ f ∂ x 2 ∣ x k cos ⁡ α 2 + ⋅ ⋅ ⋅ + ∂ f ∂ x n ∣ x k cos ⁡ α n = [ ∇ f ( x k ) ] T e \frac{\partial f(\pmb{x})}{\partial \pmb{s}} \vert _{\pmb{x}_k} = \lim\limits_{\Delta \pmb{s}\to0}( \frac{\partial f}{\partial x_{k1}}\vert _{\pmb{x}_k}\frac{\Delta x_1}{\Delta \pmb{s}}+\frac{\partial f}{\partial x_{k2}}\vert _{\pmb{x}_k}\frac{\Delta x_2}{\Delta \pmb{s}}+···+\frac{\partial f}{\partial x_{kn}}\vert _{\pmb{x}_k}\frac{\Delta x_n}{\Delta \pmb{s}}) \\ = \frac{\partial f}{\partial x_1}\vert _{\pmb{x}_k}\cos{\alpha}_1+\frac{\partial f}{\partial x_2}\vert _{\pmb{x}_k}\cos{\alpha}_2+···+\frac{\partial f}{\partial x_n}\vert _{\pmb{x}_k}\cos{\alpha}_n \\ = [\nabla f(\pmb{x}_k)]^T \pmb{e} ∂s∂f(x)​∣xk​​=Δs→0lim​(∂xk1​∂f​∣xk​​ΔsΔx1​​+∂xk2​∂f​∣xk​​ΔsΔx2​​+⋅⋅⋅+∂xkn​∂f​∣xk​​ΔsΔxn​​)=∂x1​∂f​∣xk​​cosα1​+∂x2​∂f​∣xk​​cosα2​+⋅⋅⋅+∂xn​∂f​∣xk​​cosαn​=[∇f(xk​)]Te

此处， cos ⁡ α i = lim ⁡ Δ s → 0 Δ x i Δ s ( i = 1 , 2 , ⋅ ⋅ ⋅ , n ) \cos{\alpha}_i=\lim\limits_{\Delta \pmb{s}\to0}\frac{\Delta x_i}{\Delta \pmb{s}}(i=1,2,···,n) cosαi​=Δs→0lim​ΔsΔxi​​(i=1,2,⋅⋅⋅,n)，表示方向 s \pmb{s} s与坐标轴 x i 夹 \pmb{x}_i夹 xi​夹角的余弦值。 [ ∇ f ( x k ) ] [\nabla f(\pmb{x}_k)] [∇f(xk​)]是函数 f f f在 x k \pmb{x}_k xk​处的梯度向量，表达式为

[ ∇ f ( x k ) ] T = [ ∂ f ∂ x k 1 ∣ x k , ∂ f ∂ x k 2 ∣ x k , ⋅ ⋅ ⋅ , ∂ f ∂ x k n ∣ x k ] [\nabla f(\pmb{x}_k)]^T=[\frac{\partial f}{\partial x_{k1}}\vert _{\pmb{x}_k},\frac{\partial f}{\partial x_{k2}}\vert _{\pmb{x}_k},···,\frac{\partial f}{\partial x_{kn}}\vert _{\pmb{x}_k}] [∇f(xk​)]T=[∂xk1​∂f​∣xk​​,∂xk2​∂f​∣xk​​,⋅⋅⋅,∂xkn​∂f​∣xk​​]

e = [ cos ⁡ α 1 , cos ⁡ α 2 , ⋅ ⋅ ⋅ , cos ⁡ α n ] \pmb{e}=[\cos{\alpha _1},\cos{\alpha _2},···,\cos{\alpha _n}] e=[cosα1​,cosα2​,⋅⋅⋅,cosαn​]，且 ∥ e ∥ = 1 \Vert \pmb{e} \Vert=1 ∥e∥=1，所以 e \pmb{e} e是 s \pmb{s} s方向上的单位向量。

在 x k \pmb{x_k} xk​处，由于梯度向量是固定的，而 e \pmb{e} e的模值为1，所以变化率的值仅随梯度向量和 e \pmb{e} e之间的相对空间关系变化。当两者为同一方向时，变化率为最大正值；当两者为相反方向时，变化率为最小负值。

在最速下降法中，就设定

d k = − [ ∇ f ( x k ) ] \pmb{d}_k=-[\nabla f(\pmb{x}_k)] dk​=−[∇f(xk​)]

牛顿法

看起来，最速下降法已经非常优秀了，毕竟都已经沿着让函数下降最快的方向了，还能更好的设计方案嘛？

答案是有的。单从下一步来看，负梯度方向确实是最好的；但是如果看两步呢？两次负梯度方向的矢量和能否通过一次计算得到？甚至是，出现比两次负梯度方向矢量和更优的方向？

下面看一下牛顿法是如何解决以上问题的。

要看两步，就需要用到 x \pmb{x} x的二次项。针对任意函数 f ( x ) f(\pmb{x}) f(x)，使用泰勒公式展开，并保留二次项

f ( x ) = f ( x k ) + [ ∇ f ( x k ) ] T [ x − x k ] + 1 2 [ x − x k ] T H ( x k ) [ x − x k ] f(\pmb{x})=f(\pmb{x_k})+[\nabla f(\pmb{x}_k)]^T[\pmb{x}-\pmb{x_k}]+\frac{1}{2}[\pmb{x}-\pmb{x_k}]^T\pmb{H}(\pmb{x}_k)[\pmb{x}-\pmb{x_k}] f(x)=f(xk​)+[∇f(xk​)]T[x−xk​]+21​[x−xk​]TH(xk​)[x−xk​]

其中， H ( x k ) \pmb{H}(\pmb{x}_k) H(xk​)为海森矩阵。

上式的一阶导数为

∇ f ( x ) = [ ∇ f ( x k ) ] + H ( x k ) [ x − x k ] \nabla f(\pmb{x})=[\nabla f(\pmb{x}_k)]+\pmb{H}(\pmb{x}_k)[\pmb{x}-\pmb{x_k}] ∇f(x)=[∇f(xk​)]+H(xk​)[x−xk​]

令导数值为0，并写成迭代形式，得到

x k + 1 = x k − [ H ( x k ) ] − 1 [ ∇ f ( x k ) ] \pmb{x}_{k+1}=\pmb{x}_k-[\pmb{H}(\pmb{x}_k)]^{-1}[\nabla f(\pmb{x}_k)] xk+1​=xk​−[H(xk​)]−1[∇f(xk​)]

以上即为牛顿法。相比最速下降法，牛顿法在 [ ∇ f ( x k ) ] [\nabla f(\pmb{x}_k)] [∇f(xk​)]的前面，多了一项 [ H ( x k ) ] − 1 [\pmb{H}(\pmb{x}_k)]^{-1} [H(xk​)]−1。

牛顿法因为利用了泰勒展开式，所以当 x k \pmb{x}_k xk​接近（局部）最优解时，收敛速度非常快。

拟牛顿法

既然牛顿法相比最速下降法，可以理解为多看了一步，那是不是照葫芦画瓢，再多看几步，然后不断更新迭代方向呢？原则上是可以的，但是牛顿法已经有自己的问题了：（1）函数 f ( x ) f(\pmb x) f(x)必须二阶可导；（2） H \pmb{H} H和 H − 1 \pmb{H}^{-1} H−1的计算过程较为复杂。再多看几步，对函数 f ( x ) f(\pmb x) f(x)的要求只会更高，所以更靠谱的优化方案是去降低 H − 1 \pmb{H}^{-1} H−1的计算复杂度。

接下来要介绍的拟牛顿法，其主要思路是通过构造复杂度低的函数来替代 H − 1 \pmb{H}^{-1} H−1，从而降低复杂度。

把牛顿法推导过程中用到的一阶导数公式再抄一遍，并写成迭代形式

∇ f ( x k + 1 ) = [ ∇ f ( x k ) ] + H ( x k ) [ x k + 1 − x k ] \nabla f(\pmb{x}_{k+1})=[\nabla f(\pmb{x}_k)]+\pmb{H}(\pmb{x}_k)[\pmb{x}_{k+1}-\pmb{x_k}] ∇f(xk+1​)=[∇f(xk​)]+H(xk​)[xk+1​−xk​]

令 y k = [ ∇ f ( x k + 1 ) ] − [ ∇ f ( x k ) ] \pmb{y}_k=[\nabla f(\pmb{x}_{k+1})]-[\nabla f(\pmb{x}_k)] yk​=[∇f(xk+1​)]−[∇f(xk​)]， s k = x k + 1 − x k \pmb{s}_k=\pmb{x}_{k+1}-\pmb{x_k} sk​=xk+1​−xk​，上式可以化简为

s k = [ H ( x k ) ] − 1 y k \pmb{s}_k=[\pmb{H}(\pmb{x}_k)]^{-1}\pmb{y}_k sk​=[H(xk​)]−1yk​

令 [ H ( x k ) ] − 1 = G ( x k ) [\pmb{H}(\pmb{x}_k)]^{-1}=\pmb{G}(\pmb{x}_k) [H(xk​)]−1=G(xk​)，并构造迭代公式

G k + 1 = G k + E k \pmb{G}_{k+1}=\pmb{G}_{k}+\pmb{E}_{k} Gk+1​=Gk​+Ek​

如果设定 G 0 \pmb{G}_0 G0​为单位阵，只需要确定 E k \pmb{E}_{k} Ek​，便可以替代海森矩阵了。

令 E k = α μ k μ k T + β ν k ν k T \pmb{E}_{k}=\alpha \pmb{\mu}_k \pmb{\mu}_k^T+\beta \pmb{\nu}_k \pmb{\nu}_k^T Ek​=αμk​μkT​+βνk​νkT​，其中 μ k \pmb{\mu}_k μk​和 ν k \pmb{\nu}_k νk​均为 n × 1 n\times1 n×1向量(构造成两项的原因会在后面描述)。

s k = ( G k + α μ k μ k T + β ν k ν k T ) y k \pmb{s}_k=(\pmb{G}_{k}+\alpha \pmb{\mu}_k \pmb{\mu}_k^T+\beta \pmb{\nu}_k \pmb{\nu}_k^T)\pmb{y}_k sk​=(Gk​+αμk​μkT​+βνk​νkT​)yk​

上式做一下变换

α ( μ k T y k ) μ k + β ( ν k T y k ) ν k = s k − G k y k \alpha (\pmb{\mu}_k^T \pmb{y}_k)\pmb{\mu}_k+\beta (\pmb{\nu}_k^T \pmb{y}_k)\pmb{\nu}_k = \pmb{s}_k-\pmb{G}_{k}\pmb{y}_k α(μkT​yk​)μk​+β(νkT​yk​)νk​=sk​−Gk​yk​

其中， μ k T y k \pmb{\mu}_k^T \pmb{y}_k μkT​yk​和 ν k T y k \pmb{\nu}_k^T \pmb{y}_k νkT​yk​为实数， s k − G k \pmb{s}_k-\pmb{G}_{k} sk​−Gk​为 n × 1 n\times1 n×1向量， α \alpha α和 β \beta β可以任意选取，我们取特殊的一组： μ k = r G k y k \pmb{\mu}_k=r\pmb{G}_{k}\pmb{y}_k μk​=rGk​yk​， ν k = θ s k \pmb{\nu}_k=\theta\pmb{s}_k νk​=θsk​，此时 E k \pmb{E}_k Ek​变为

E k = α r 2 G k y k y k T G k T + β θ 2 s k s k T \pmb{E}_{k}=\alpha r^2\pmb{G}_{k}\pmb{y}_k\pmb{y}_k^T\pmb{G}_{k}^T+\beta \theta^2\pmb{s}_k\pmb{s}_k^T Ek​=αr2Gk​yk​ykT​GkT​+βθ2sk​skT​

将 μ k \pmb{\mu}_k μk​和 ν k \pmb{\nu}_k νk​带入上上式

α r 2 ( y k T G k T y k ) G k y k + β θ 2 ( s k T y k ) s k = s k − G k y k \alpha r^2(\pmb{y}_k^T\pmb{G}_{k}^T \pmb{y}_k)\pmb{G}_{k}\pmb{y}_k+\beta \theta^2(\pmb{s}_k^T \pmb{y}_k)\pmb{s}_k = \pmb{s}_k-\pmb{G}_{k}\pmb{y}_k αr2(ykT​GkT​yk​)Gk​yk​+βθ2(skT​yk​)sk​=sk​−Gk​yk​

化简一下

[ α r 2 ( y k T G k T y k ) G k + 1 ] G k y k + [ β θ 2 s k T y k − 1 ] s k = 0 [\alpha r^2(\pmb{y}_k^T\pmb{G}_{k}^T \pmb{y}_k)\pmb{G}_{k}+1]\pmb{G}_{k}\pmb{y}_k+[\beta \theta^2\pmb{s}_k^T \pmb{y}_k-1]\pmb{s}_k=0 [αr2(ykT​GkT​yk​)Gk​+1]Gk​yk​+[βθ2skT​yk​−1]sk​=0

由于 G k y k \pmb{G}_{k}\pmb{y}_k Gk​yk​和 s k \pmb{s}_k sk​是任意的，所以需要 α r 2 ( y k T G k T y k ) G k + 1 = 0 \alpha r^2(\pmb{y}_k^T\pmb{G}_{k}^T \pmb{y}_k)\pmb{G}_{k}+1=0 αr2(ykT​GkT​yk​)Gk​+1=0和 β θ 2 s k T y k − 1 = 0 \beta \theta^2\pmb{s}_k^T \pmb{y}_k-1=0 βθ2skT​yk​−1=0，得到

α r 2 = − 1 y k T G k T y k , β θ 2 = 1 s k T y k \alpha r^2=-\frac{1}{\pmb{y}_k^T\pmb{G}_{k}^T \pmb{y}_k},\beta \theta^2=\frac{1}{\pmb{s}_k^T \pmb{y}_k} αr2=−ykT​GkT​yk​1​,βθ2=skT​yk​1​

从上上式可以看出，如果 E k \pmb{E}_k Ek​在构造时只有一项，是无法保证恒等式的。

所以，最终的拟牛顿法迭代公式为

G k + 1 = G k − G k y k y k T G k T y k T G k T y k + s k s k T s k T y k \pmb{G}_{k+1}=\pmb{G}_{k}-\frac{\pmb{G}_{k}\pmb{y}_k\pmb{y}_k^T\pmb{G}_{k}^T}{\pmb{y}_k^T\pmb{G}_{k}^T \pmb{y}_k}+\frac{\pmb{s}_k\pmb{s}_k^T}{\pmb{s}_k^T \pmb{y}_k} Gk+1​=Gk​−ykT​GkT​yk​Gk​yk​ykT​GkT​​+skT​yk​sk​skT​​

相比 H − 1 \pmb{H}^{-1} H−1，上式的计算不仅复杂度低，而且不需要 f ( x ) f(\pmb x) f(x)为二阶可导。