算法细节系列（3）：梯度下降法，牛顿法，拟牛顿法
算法细节系列（3）：梯度下降法，牛顿法，拟牛顿法
算法细节系列（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]x\in[a,b]时，a≤ϕ(x)≤b;a \le \phi(x) \le b; (2) 存在正数L<1L \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) 公式所表达的含义很简单，首先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}}]_{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 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()输出结果为：
初始点为:
[-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是如何更新和选取的呢？不解，等学习到具体的拟牛顿方法再来完善吧。