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

迭代算法原型

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

x=ϕ(x)

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

f(x)=x3−x−1

f(x) = x^3 - x -1
该怎么做？中学的方法很简单，多项式方程吗，令 f(x)=0 f(x) = 0，利用公式解出来即可。那求 f(x)=x3−logx−1 f(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.5 x_0 = 1.5，代入公式的右侧得 x1=x0+1−−−−−√3=1.35721 x_1 = \sqrt[3]{x_0 + 1} = 1.35721，在把求得的 x1 x_1代入右式，如此往复，如下表所示：

k	xk x_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

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

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

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

设迭代函数 ϕ(x)在[a,b] \phi(x)在[a,b]上具有连续的一阶导数，且
(1) 当 x∈[a,b] x\in[a,b]时， a≤ϕ(x)≤b; a \le \phi(x) \le b;
(2) 存在正数 L<1 L \lt 1，使对任意 x∈[a,b] x\in [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)
公式所表达的含义很简单，首先 x x是一个迭代的过程，即自己再不断地代入，直到求得的解收敛。当初在证明解收敛时是非常感性的认识，我们是要求解f(x)f(x)的极值，不管是极小值还是极大值，通过上述公式，在初值 x0 x_0的情况下，只要 α \alpha满足一定的条件，最终 x x将收敛，且是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)=0 f(x) = 0的一种重要和常用的迭代法，它的基本思想是将非线性函数 f(x) f(x)逐步线性化，从而将非线性方程 f(x)=0 f(x) = 0近似地转化为线性方程求解。它对应的方程为：

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

x = x - \frac{f(x)}{f'(x)}
牛顿法有明显的几何意义。方程 f(x)=0 f(x) = 0的根 x∗ x^*在几何上表示曲线 y=f(x) y = f(x)与 x x轴的交点。当我们求得x∗x^*的近似值 xk x_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)
则该切线与 x x轴交点横坐标正好是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}}]_{n\times 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代码如下：

"""
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 pltdef 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 = 1imax = 1000W[:,0] = x0 x = x0delta = 1alpha = 1while i<imax and delta>10**(-5):p = -np.dot(np.linalg.inv(hessian(x)),jacobian(x))x0 = xx = x + alpha*pW[:,i] = xdelta = sum((x-x0)**2)print('第',i,'次迭代结果:')print(x,'\n')i=i+1W=W[:,0:i]  # 记录迭代点return Wx0 = np.array([-1.2,1])
W=newton(x0)plt.plot(W[0,:],W[1,:],'g*',W[0,:],W[1,:]) # 画出迭代点收敛的轨迹
plt.show()

输出结果为：

初始点为:
[-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是如何更新和选取的呢？不解，等学习到具体的拟牛顿方法再来完善吧。

参考文献

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