算法细节系列（3）：梯度下降法，牛顿法，拟牛顿法

2019-05-26 19:43:11 浏览数 (1)

算法细节系列（3）：梯度下降法，牛顿法，拟牛顿法

迭代算法原型

话不多说，直接进入主题。在我看来，不管是梯度下降法还是牛顿法，它们都可以归结为一个式子，即

x=ϕ(x)

x = phi(x) 也就是我们的不动点迭代法（fixed pointed iteration）最核心的迭代公式。神奇的式子，它该如何操作呢？用来干什么呢？不动点迭代法主要用于求解函数的零点。如求以下函数的零点，

f(x)=x3−x−1

f(x) = x^3 - x -1 该怎么做？中学的方法很简单，多项式方程吗，令f(x)=0f(x) = 0，利用公式解出来即可。那求f(x)=x3−logx−1f(x) = x^3 - log x -1的零点呢？貌似就难求了，没关系，不动点迭代法就是用来求解这些超越方程的，或者说可以用计算的方法，不断迭代逼近正确值。如下，把f(x)f(x)写成如上形式x=ϕ(x)x = phi(x)，即：

x=x 1−−−−−√3

x = sqrt[3]{x 1} 为了方便迭代我们重新改写下公式，即

xk 1=xk 1−−−−−√3

x_{k 1} = sqrt[3]{x_k 1} 好了，有了这个最基础的公式我们就可以迭代了，如下随意取一个值x0=1.5x_0 = 1.5，代入公式的右侧得x1=x0 1−−−−−√3=1.35721x_1 = sqrt[3]{x_0 1} = 1.35721，在把求得的x1x_1代入右式，如此往复，如下表所示：

k

xkx_k

0

1.5

1

1.35721

2

1.33085

3

1.32588

4

1.32494

5

1.32475

6

1.32473

7

1.32472

8

1.32742

经过几次迭代后，你会发现一个神奇的现象，即在指定误差内，xkx_k将收敛于一个稳定值。而这个稳定值从直观上来理解，就是待求函数f(x)f(x)的零点。因为左边和右边均等于0，为该方程的唯一解。该迭代过程可以用几何图形直观的感觉一下，如下图：

你可能要问了，初始值的取法是随意的吗？该过程最终一定会收敛么？可以明确的告诉你，不动点迭代法对其初始值是敏感的，并不是任何值都能让该过程收敛。所以从数学的角度，是需要研究收敛性的充分必要条件是什么，方便更好得使用该方法。

这里直接给出在某区域内满足不动点迭代收敛性的充分必要条件，证明可以参考《数值分析》李庆扬版。

设迭代函数ϕ(x)在[a,b]phi(x)在[a,b]上具有连续的一阶导数，且 (1) 当x∈[a,b]xin[a,b]时，a≤ϕ(x)≤b;a le phi(x) le b; (2) 存在正数L<1L lt 1，使对任意x∈[a,b]xin [a,b]，有

|ϕ′(x)|≤L<1成立

vert phi'(x)vert le L lt 1 成立

书上的证明还是相对比较繁琐的，简单说一下这两个条件在几何中的意义。第一个是为了说明x=ϕ(x)x = phi(x)在区间[a,b][a,b]中存在唯一解，第二个条件则说明了在区间[a,b][a,b]之内的所有点的斜率小于1，让迭代过程收敛。在图上画几条线比划比划下就能明了，这里就不再过多的阐述了。

细心的同学读到这里有没有发现，其实x=ϕ(x)x = phi(x)这种表达式可以有多种，比如同样的一个问题，我们可以以

x=x3−1

x = x^3 -1 进行迭代，如果符合上述两个收敛条件，那么它就是一个有效的迭代方程。所以这里就引出了另外一点，ϕ(x)phi(x)的选取是有讲究的，试想一下，不同的ϕ(x)phi(x)是否能够影响收敛速度？

梯度下降法

刚接触梯度下降法时，是在上吴恩达Stanford机器学习公开课时了解的，梯度下降法的公式很简单，为

x:=x−αf′(x)

x := x - alpha f'(x) 公式所表达的含义很简单，首先xx是一个迭代的过程，即自己再不断地代入，直到求得的解收敛。当初在证明解收敛时是非常感性的认识，我们是要求解f(x)f(x)的极值，不管是极小值还是极大值，通过上述公式，在初值x0x_0的情况下，只要αalpha满足一定的条件，最终xx将收敛，且是f(x)f(x)的局部极值。

现在利用不动点的迭代公式似乎就显得更加有理有据了。首先f(x)f(x)的极值问题可以表示为求解f′(x)f'(x)的零点问题。那么我们只要凑出

x=ϕ(x)

x =phi(x) 即可，而梯度下降算法，凑的公式为：

x=x−αf′(x)→f′(x)=0

x = x - alpha f'(x) rightarrow f'(x) = 0 所以只要有了f′(x)f'(x)的表达式，我们既而能求出αalpha的取值范围。详细的可以参看知乎上的一篇回答，还是比较不错的。最优化问题中，牛顿法为什么比梯度下降法求解需要的迭代次数更少？

牛顿法

牛顿迭代法是求解非线性方程f(x)=0f(x) = 0的一种重要和常用的迭代法，它的基本思想是将非线性函数f(x)f(x)逐步线性化，从而将非线性方程f(x)=0f(x) = 0近似地转化为线性方程求解。它对应的方程为：

x=x−f(x)f′(x)

x = x - frac{f(x)}{f'(x)} 牛顿法有明显的几何意义。方程f(x)=0f(x) = 0的根x∗x^*在几何上表示曲线y=f(x)y = f(x)与xx轴的交点。当我们求得x∗x^*的近似值xkx_k以后，过曲线y=f(x)y = f(x)上对应点(xk,f(xk))(x_k,f(x_k))作曲线f(x)f(x)的切线，其切线方程为：

y−f(xk)=f′(xk)(x−xk)

y - f(x_k) = f'(x_k)(x- x_k) 则该切线与xx轴交点横坐标正好是xk−f(xk)f′(xk)x_k - frac{f(x_k)}{f'(x_k)}，这就是牛顿迭代公式的计算结果。

牛顿迭代法在几何图形上的意义也是显而易见的。它的收敛速度比梯度下降算法要快得多，这里我们也不去证明了，书中主要应用了一个新的定义来论证两者的收敛速度，叫收敛阶，有兴趣的可以继续研究。

其实牛顿迭代法也是特殊的不动点迭代过程。即迭代方程依旧是x=ϕ(x)x = phi (x)，可见，不管是牛顿迭代还是梯度下降，迭代的核心思想并没有发生变化，变得只是ϕ(x)phi(x)而已。如果要求解f(x)f(x)的极值呢？很简单，牛顿迭代方程稍微做下改动即可：

x=x−f′(x)f′′(x)

x = x - frac{f'(x)}{f''(x)} 同样，我们能得到收敛的唯一解。上述所讲的迭代方法都是一维简单情形下的，接下来就推广牛顿迭代法在多维的情况。

海赛矩阵（hesse matrix）

对于高维函数，用牛顿法求极值也是这个形式，只不过这里的f′′(x)f''(x)和f′(x)f'(x)都变成了矩阵和向量。而且你也可以想到，高维函数二阶导有多个，写成矩阵的形式就是这样

H(x)=[∂2f∂xi∂xj]n×n=⎡⎣⎢⎢f11⋮fn1⋯⋱⋯f1n⋮fnn⎤⎦⎥⎥

H(x) = [frac{partial^2f}{partial{x_i}partial{x_j}}]_{ntimes n} = begin{bmatrix} f_{11} & cdots & f_{1n} \ vdots & ddots & vdots \ f_{n1} & cdots & f_{nn} \ end{bmatrix} 所以有了海赛矩阵，牛顿迭代法就有如下形式：

X:=X−H−1(X)∇f(X),X=x1,x2,...,xn

X := X - H^{-1}(X)nabla f(X), X = x_1,x_2,...,x_n 对于H(x)H(x)有如下结论：

  • 如果H(x)是正定矩阵，则临界点M处是一个局部的极小值。
  • 如果H(x)是负定矩阵，则临界点M处是一个局部的极大值。
  • 如果H(x)是不定矩阵，则临界点M处不是极值。

Code Time

我们以一个多元函数为例，用代码来求解该函数的极值。如函数：

f(x)=100(x2−x21)2 (1−x1)2

f(x) = 100(x_2-x_1^2)^2 (1- x_1)^2 于是可得函数的梯度：

∇f(x)=(−400(x2−x21)x1−2(1−x1),200(x2−x21))T

nabla f(x) = (-400(x_2-x_1^2)x_1-2(1-x_1),200(x_2-x_1^2))^T 函数f(x)f(x)的海赛矩阵为

H(x)=[−400(x2−3x21) 2−400x1−400x1200]

H(x) = begin{bmatrix} -400(x_2-3x_1^2) 2 & -400x_1 \ -400x_1 & 200 \ end{bmatrix}

编写Python代码如下：

代码语言：javascript复制
"""
Newton法
Rosenbrock函数
函数 f(x)=100*(x(2)-x(1).^2).^2 (1-x(1)).^2
梯度 g(x)=(-400*(x(2)-x(1)^2)*x(1)-2*(1-x(1)),200*(x(2)-x(1)^2))^(T)
"""

import numpy as np
import matplotlib.pyplot as plt

def jacobian(x):
    return np.array([-400*x[0]*(x[1]-x[0]**2)-2*(1-x[0]),200*(x[1]-x[0]**2)])

def hessian(x):
    return np.array([[-400*(x[1]-3*x[0]**2) 2,-400*x[0]],[-400*x[0],200]])

X1=np.arange(-1.5,1.5 0.05,0.05)
X2=np.arange(-3.5,2 0.05,0.05)
[x1,x2]=np.meshgrid(X1,X2)
f=100*(x2-x1**2)**2 (1-x1)**2; # 给定的函数
plt.contour(x1,x2,f,20) # 画出函数的20条轮廓线


def newton(x0):

    print('初始点为:')
    print(x0,'n')
    W=np.zeros((2,10**3))
    i = 1
    imax = 1000
    W[:,0] = x0 
    x = x0
    delta = 1
    alpha = 1

    while i<imax and delta>10**(-5):
        p = -np.dot(np.linalg.inv(hessian(x)),jacobian(x))
        x0 = x
        x = x   alpha*p
        W[:,i] = x
        delta = sum((x-x0)**2)
        print('第',i,'次迭代结果:')
        print(x,'n')
        i=i 1
    W=W[:,0:i]  # 记录迭代点
    return W

x0 = np.array([-1.2,1])
W=newton(x0)

plt.plot(W[0,:],W[1,:],'g*',W[0,:],W[1,:]) # 画出迭代点收敛的轨迹
plt.show()

输出结果为：

代码语言：javascript复制
初始点为:
[-1.2  1. ] 

第 1 次迭代结果:
[-1.1752809   1.38067416] 

第 2 次迭代结果:
[ 0.76311487 -3.17503385] 

第 3 次迭代结果:
[ 0.76342968  0.58282478] 

第 4 次迭代结果:
[ 0.99999531  0.94402732] 

第 5 次迭代结果:
[ 0.9999957   0.99999139] 

第 6 次迭代结果:
[ 1.  1.]

即迭代了6次得到了最优解，画出的迭代点的轨迹如下：

上述内容摘自博文用Python实现牛顿法求极值。

拟牛顿法

摘自博文牛顿法与拟牛顿法学习笔记（二）拟牛顿条件

神奇，根据gk 1−gk=Hk(xk 1−xk)g_{k 1}-g_{k} = H_k(x_{k 1}-x_k)，就得到了拟牛顿的条件了， 为何选择它们之间的差作为约束条件来选取近似的海赛矩阵呢？其次，按照拟牛顿条件D是如何更新和选取的呢？不解，等学习到具体的拟牛顿方法再来完善吧。

参考文献

  1. 最优化问题中，牛顿法为什么比梯度下降法求解需要的迭代次数更少？
  2. 用Python实现牛顿法求极值。
  3. 牛顿法与拟牛顿法学习笔记（二）拟牛顿条件
python alpha partial

0 人点赞