牛顿法与拟牛顿法(含代码实现)

1. 牛顿法

牛顿法（英语：Newton’s method）又称为牛顿-拉弗森方法（英语：Newton-Raphson method），它是一种在实数域和复数域上近似求解方程的方法。

牛顿法的基本思想是使用函数 f(x){\displaystyle f(x)}f(x) 的泰勒级数的前面几项来寻找方程 f(x)=0{\displaystyle f(x)=0}f(x)=0 的根。

牛顿法主要应用在两个方面，1：求方程的根；2：最优化(求解最值问题)。

1.1 求方程f(x)=0{\displaystyle f(x)} = 0f(x)=0的根

选择一个接近函数f(x){\displaystyle f(x)}f(x)零点的x0\mathbf x_0x0, 牛顿法对函数进行一阶泰勒展开：
f(x)=f(x0)+f′(x0)(x−x0)f(x)=f\left(x_{0}\right)+f^{\prime}\left(x_{0}\right)\left(x-x_{0}\right) f(x)=f(x0)+f′(x0)(x−x0)

由f(x)=0{\displaystyle f(x)} = 0f(x)=0，得到迭代方式：
xn+1=xn−f(xn)/f′(xn)x_{n+1}=x_{n}-f\left(x_{n}\right) / f^{\prime}\left(x_{n}\right) xn+1=xn−f(xn)/f′(xn)
迭代后求得方程的根x∗x^*x∗：
f(x∗)=0f\left(x^{*}\right)=0 f(x∗)=0

牛顿法有一个性质，就是能够保证二次收敛到方程的根。论述如下：

假设函数fff在开区间(a,b)(a,b)(a,b)是二阶可导的,并存在函数的根x∗∈(a,b)x^*\in(a,b)x∗∈(a,b).定义牛顿迭代法：
xk+1=xk−f(xk)f′(xk),k=1,2,…x_{k+1}=x_{k}-\frac{f\left(x_{k}\right)}{f^{\prime}\left(x_{k}\right)}, \quad k=1,2, \ldots xk+1=xk−f′(xk)f(xk),k=1,2,…
假设k→∞k \rightarrow \inftyk→∞时，xkx_kxk收敛到x∗x^*x∗。若f′(x∗)≠0f^{\prime}\left(x^{*}\right) \neq 0f′(x∗)=0，对于足够大的 kkk，有：
∣xk+1−x∗∣≤M∣xk−x∗∣2if M>∣f′′(x∗)∣2∣f′(x∗)∣\left|x_{k+1}-x^{*}\right| \leq M\left|x_{k}-x^{*}\right|^{2} \quad \text { if } M>\frac{\left|f^{\prime \prime}\left(x^{*}\right)\right|}{2\left|f^{\prime}\left(x^{*}\right)\right|} ∣xk+1−x∗∣≤M∣xk−x∗∣2 if M>2∣f′(x∗)∣∣f′′(x∗)∣
于是，xkx_kxk是二次收敛到x∗x^*x∗。

证明：
假设ek=xk−x∗e_k = x_k - x^*ek=xk−x∗, 即 $x^*= x_k - e_k $，根据泰勒展开公式，
f(xk−ek)=f(xk)−ekf′(xk)+(ek)22f′′(ξk)f\left(x_{k}-e_{k}\right)=f\left(x_{k}\right)-e_{k} f^{\prime}\left(x_{k}\right)+\frac{\left(e_{k}\right)^{2}}{2} f^{\prime \prime}\left(\xi_{k}\right) f(xk−ek)=f(xk)−ekf′(xk)+2(ek)2f′′(ξk)
其中，ξk\xi_kξk介于xkx_kxk和x∗x^*x∗之间。
由于f(x∗)=0f(x^*) = 0f(x∗)=0,所以有
0=f(xk)−(xk−x∗)f′(xk)+(ek)22f′′(ξk)0=f\left(x_{k}\right)-\left(x_{k}-x^{*}\right) f^{\prime}\left(x_{k}\right)+\frac{\left(e_{k}\right)^{2}}{2} f^{\prime \prime}\left(\xi_{k}\right) 0=f(xk)−(xk−x∗)f′(xk)+2(ek)2f′′(ξk)
由于函数fff连续可导，且f′(x∗)≠0f^{\prime}\left(x^{*}\right) \neq 0f′(x∗)=0，只要xkx_kxk和x∗x^*x∗足够接近，则有f′(xk)≠0f^{\prime}\left(x_{k}\right) \neq 0f′(xk)=0。两边除于f′(xk)f^{\prime}\left(x_{k}\right)f′(xk)得到
0=f(xk)f′(xk)−(xk−x∗)+(ek)2f′′(ξk)2f′(xk)0=\frac{f\left(x_{k}\right)}{f^{\prime}\left(x_{k}\right)}-\left(x_{k}-x^{*}\right)+\frac{\left(e_{k}\right)^{2} f^{\prime \prime}\left(\xi_{k}\right)}{2 f^{\prime}\left(x_{k}\right)} 0=f′(xk)f(xk)−(xk−x∗)+2f′(xk)(ek)2f′′(ξk)
根据牛顿迭代法的定义，可以得到下式：
xk+1−x∗=(ek)2f′′(ξk)2f′(xk)x_{k+1}-x^{*}=\frac{\left(e_{k}\right)^{2} f^{\prime \prime}\left(\xi_{k}\right)}{2 f^{\prime}\left(x_{k}\right)} xk+1−x∗=2f′(xk)(ek)2f′′(ξk)
所以有，
∣xk+1−x∗∣≤∣f′′(ξk)∣2∣f′(xk)∣∣xk−x∗∣2\left|x_{k+1}-x^{*}\right| \leq \frac{\left|f^{\prime \prime}\left(\xi_{k}\right)\right|}{2\left|f^{\prime}\left(x_{k}\right)\right|}\left|x_{k}-x^{*}\right|^{2} ∣xk+1−x∗∣≤2∣f′(xk)∣∣f′′(ξk)∣∣xk−x∗∣2

一般地，f′(xk)f^{\prime}\left(x_{k}\right)f′(xk)会收敛到f′(x∗)f^{\prime}\left(x^*\right)f′(x∗)，由于ξk\xi_kξk介于xkx_kxk和x∗x^*x∗之间，因此，ξk\xi_kξk收敛到x∗x^*x∗，f′′(ξk)f^{\prime\prime}\left(\xi_{k}\right)f′′(ξk)收敛到f′′(x∗)f^{\prime\prime}\left(x^*\right)f′′(x∗)。对应足够大的kkk，有
∣xk+1−x∗∣≤M∣xk−x∗∣2if M>∣f′′(x∗)∣2∣f′(x∗)∣\left|x_{k+1}-x^{*}\right| \leq M\left|x_{k}-x^{*}\right|^{2} \quad \text { if } M>\frac{\left|f^{\prime \prime}\left(x^{*}\right)\right|}{2\left|f^{\prime}\left(x^{*}\right)\right|} ∣xk+1−x∗∣≤M∣xk−x∗∣2 if M>2∣f′(x∗)∣∣f′′(x∗)∣

证毕。

1.2 最优化（求取极值）

解决最优化问题min⁡x∈Rnf(x)\min\limits_{x \in \mathbf{R}^{n}} f(x)x∈Rnminf(x)的结构：

给定初始点x0\mathbb x_0x0,

确定搜索方向dk\mathbb d_kdk，即依照一定规则构造fff在xk\mathbb x_kxk点处的下降方向为搜索方向;
确定步长因子αkα_kαk，使目标函数值有某种意义下降;
令xk+1=xk+αkdk\mathbb x_{k+1}=\mathbb x_k+ α_k \mathbb d_kxk+1=xk+αkdk

a) 若xk+1\mathbb x_{k+1}xk+1满足某种终止条件，则停止迭代，得到近似最优解，

b) 否则，重复以上步骤。

牛顿法解决最优化问题的基本思想是利用目标函数的二次Taylor展开，并将其极小化。

假设目标函数f(x)f(x)f(x)具有二阶连续偏导数，x∗x^{*}x∗ 为目标函数的极小点，对目标函数在第kkk次迭代值进行二阶泰勒展开：

f(x)=f(xk)+gkT(x−xk)+12(x−xk)TH(xk)(x−xk)f(x)=f\left(x_{k}\right)+g_{k}^{T}\left(x-x_{k}\right)+\frac{1}{2}\left(x-x_{k}\right)^{T} H\left(x_{k}\right)\left(x-x_{k}\right) f(x)=f(xk)+gkT(x−xk)+21(x−xk)TH(xk)(x−xk)
其中，gk=g(xk)=∇f(xk)g_{k}=g\left(x_{k}\right)=\nabla f\left(x_{k}\right)gk=g(xk)=∇f(xk)是f(x)f(x)f(x)在xkx_{k}xk的一阶导数值，H(xk)H\left(x_{k}\right)H(xk)是f(x)f(x)f(x)的海森矩阵：
H(x)=[∂2f∂xi∂xj]n×nH(x)=\left[\frac{\partial^{2} f}{\partial x_{i} \partial x_{j}}\right]_{n \times n} H(x)=[∂xi∂xj∂2f]n×n

函数 f(x)f(x)f(x)有极值的必要条件是在极值点处一阶导数为0。特别的当 H(xk)H(x_k)H(xk)是正定矩阵时，函数f(x)f(x)f(x)的极值为极小值。

对方程∇f(x)=0\nabla f(x)=0∇f(x)=0,根据上述牛顿迭代法可以求解，也可以二阶泰勒展开公式再进行求导：
∇f(x)=gk+Hk(x−xk)\nabla f(x)=g_{k}+H_{k}\left(x-x_{k}\right) ∇f(x)=gk+Hk(x−xk)
其中，记Hk=H(xk)H_{k}=H\left(x_{k}\right)Hk=H(xk),则有
gk+Hk(xk+1−xk)=0\begin{array}{l} g_{k}+H_{k}\left(x_{k+1}-x_{k}\right)=0 \end{array} gk+Hk(xk+1−xk)=0
迭代公式：
xk+1=xk−Hk−1gkx_{k+1}=x_{k}-H_{k}^{-1} g_{k}xk+1=xk−Hk−1gk

对于一元函数，上述迭代公式也可以写成：
xk+1=xk−f′(xk)f′′(xk)x_{k+1}=x_{k}-\frac{f^{\prime}\left(x_{k}\right)}{f^{\prime \prime}\left(x_{k}\right)} xk+1=xk−f′′(xk)f′(xk)

1.3 牛顿法最优化的示例

import numpy as np
from sklearn import datasets
from sklearn.linear_model import LinearRegression
class Newton(object):def __init__(self,epochs=50):self.W = Noneself.epochs = epochsdef get_loss(self, X, y, W,b):"""计算损失 0.5*sum(y_pred-y)^2input: X(2 dim np.array):特征y(1 dim np.array):标签W(2 dim np.array):线性回归模型权重矩阵output：损失函数值"""#print(np.dot(X,W))loss = 0.5*np.sum((y - np.dot(X,W)-b)**2)return lossdef first_derivative(self,X,y):"""计算一阶导数g = (y_pred - y)*xinput: X(2 dim np.array):特征y(1 dim np.array):标签W(2 dim np.array):线性回归模型权重矩阵output：损失函数值"""y_pred = np.dot(X,self.W) + self.bg = np.dot(X.T, np.array(y_pred - y))g_b = np.mean(y_pred-y)return g,g_bdef second_derivative(self,X,y):"""计算二阶导数 Hij = sum(X.T[i]*X.T[j])input: X(2 dim np.array):特征y(1 dim np.array):标签output：损失函数值"""H = np.zeros(shape=(X.shape[1],X.shape[1]))H = np.dot(X.T, X)H_b = 1return H, H_bdef fit(self, X, y):"""线性回归 y = WX + b拟合，牛顿法求解input: X(2 dim np.array):特征y(1 dim np.array):标签output：拟合的线性回归"""np.random.seed(10)self.W = np.random.normal(size=(X.shape[1]))self.b = 0for epoch in range(self.epochs):g,g_b = self.first_derivative(X,y)  # 一阶导数H,H_b = self.second_derivative(X,y)  # 二阶导数self.W = self.W - np.dot(np.linalg.pinv(H),g)self.b = self.b - 1/H_b*g_bprint("itration:{} ".format(epoch), "loss:{:.4f}".format(self.get_loss(X, y , self.W,self.b)))def normalize(x):return (x - np.min(x))/(np.max(x) - np.min(x))
if __name__ == "__main__":np.random.seed(2)X = np.random.rand(100,5)y = np.sum(X**3 + X**2,axis=1)print(X.shape, y.shape)# 归一化X_norm = normalize(X)X_train = X_norm[:int(len(X_norm)*0.8)]X_test = X_norm[int(len(X_norm)*0.8):]y_train = y[:int(len(X_norm)*0.8)]y_test = y[int(len(X_norm)*0.8):]# 牛顿法求解回归问题newton=Newton()newton.fit(X_train, y_train)y_pred = newton.predict(X_test,y_test)print(0.5*np.sum((y_test - y_pred)**2))reg = LinearRegression().fit(X_train, y_train)y_pred = reg.predict(X_test)print(0.5*np.sum((y_test - y_pred)**2))

对比了牛顿法求解线性回归和sklearn中的线性回归（最小二乘法求解）的结果：

itration:45  loss:9.584074
itration:46  loss:9.198050
itration:47  loss:8.844083
itration:48  loss:8.519513
itration:49  loss:8.221897
predict_loss 3.770617848429981 # 牛顿法 LR
3.2062598186600217#sklearn LR

如果迭代次数epoch调大，两者的效果一样：

itration:495  loss:4.935776
itration:496  loss:4.935776
itration:497  loss:4.935776
itration:498  loss:4.935776
itration:499  loss:4.935776
predict_loss 3.20625981684797
3.2062598186600217

2. 拟牛顿法

由于Hesse矩阵的计算工作量大，有时目标函数的Hesse阵很难计算。

拟牛顿法利用目标函数和一阶导数，来构造目标函数的曲率近似，而不需要明显形成Hesse阵，同时具有收敛速度快的优点。

2.1 一般拟牛顿法

由
g(x)=∇f(x)≈gk+Hk(x−xk)g(x) = \nabla f(x) \approx g_{k}+H_{k}\left(x-x_{k}\right) g(x)=∇f(x)≈gk+Hk(x−xk)
令x=xkx = x_kx=xk,得：
gk−1−gk≈Hk(xk−1−xk)g_{k-1}-g_{k} \approx H_{k}\left(x_{k-1}-x_{k}\right) gk−1−gk≈Hk(xk−1−xk)
令sk−1=xk−xk−1,yk−1=gk−gk−1\mathbf{s}_{k-1}=\mathbf{x}_{k}-\mathbf{x}_{k-1}, \mathbf{y}_{k-1}=\mathbf{g}_{k}-\mathbf{g}_{k-1}sk−1=xk−xk−1,yk−1=gk−gk−1,得
Hk−1yk−1≈sk−1H_{k}^{-1} \mathbf{y}_{k-1} \approx \mathbf{s}_{k-1} Hk−1yk−1≈sk−1
或者记为：
Hk+1−1yk≈skH_{k+1}^{-1} \mathbf{y}_{k} \approx \mathbf{s}_{k} Hk+1−1yk≈sk
对于二次函数fff,上述关系式精确成立。

要求在拟牛顿法中构造出Hesse逆近似 Gk+1G_{k+1}Gk+1（如DFP算法）,满足
Gk+1yk=skG_{k+1} \mathbf{y}_{k} = \mathbf{s}_{k}Gk+1yk=sk

或者构建Hesse近似（如BFGS算法）
Bk+1sk≈ykB_{k+1} \mathbf{s}_{k} \approx \mathbf{y}_{k} Bk+1sk≈yk

这称为拟牛顿法条件 。

一般拟牛顿法：

给定初始点x0∈Rn,H0∈Rn×n,0≤ε≤1,k=0\mathbf{x}_{0} \in R^{n}, H_{0} \in R^{n \times n}, 0 \leq \varepsilon \leq 1, k=0x0∈Rn,H0∈Rn×n,0≤ε≤1,k=0;
若∥gk∥≤ε\left\|\mathbf{g}_{k}\right\| \leq \varepsilon∥gk∥≤ε，则停止；否则，计算dk=−Hkgk\mathbf{d}_{k}=-H_{k} \mathbf{g}_{k}dk=−Hkgk；
沿方向dk\mathbb d_kdk线性搜索求步长因子αkα_kαk，令xk+1=xk+αkdk\mathbb x_{k+1}=\mathbb x_k+ α_k \mathbb d_kxk+1=xk+αkdk；
校正HkH_kHk产生Hk+1H_{k+1}Hk+1,使得拟牛顿法条件满足。
k=k+1k=k+1k=k+1，转2.

优点：

（1）只需要一阶导数；

（2）HkH_kHk保持正定，具有下降性；

（3）迭代每次需要O(n2)O(n^2)O(n2)次乘法；牛顿法是O(n3)O(n^3)O(n3)次（因为牛顿法需要求逆）。

2.2 DFP算法

设秩二校正为：
Gk+1=Gk+auuT+bvvTG_{k+1}=G_{k}+a u u^{T}+b v v^{T} Gk+1=Gk+auuT+bvvT
若要拟牛顿法条件Gk+1yk=(Gk+auuT+bvvT)yk=skG_{k+1} y_k =(G_{k}+a u u^{T}+b v v^{T}) y_k = s_kGk+1yk=(Gk+auuT+bvvT)yk=sk成立，对于u,vu,vu,v一个取法为：
u=sk,v=Gkyku=s_{k}, \quad v=G_{k} y_{k} u=sk,v=Gkyk
auTyk=1,bvTyk=−1a u^{T} y_{k}=1, b v^{T} y_{k}=-1 auTyk=1,bvTyk=−1
可以得到，
a=1skTyk,b=−1ykTGkyka=\frac{1}{s_{k}^{T} y_{k}}, b=-\frac{1}{y_{k}^{T} G_{k} y_{k}} a=skTyk1,b=−ykTGkyk1
那么，
Gk+1=Gk+skskTskTyk−GkykykTGkTykTGkykG_{k+1}=G_{k}+\frac{s_{k} s_{k}^{T}}{s_{k}^{T} y_{k}}-\frac{G_{k} y_{k} y_{k}^{T} G_{k}^{T}}{y_{k}^{T} G_{k} y_{k}} Gk+1=Gk+skTykskskT−ykTGkykGkykykTGkT

2.2 L-BFGS算法

BFGS算法推导和DFP算法差不多，直接给出公式：
Bk+1=Bk+ykykTykTsk−BkskskTBkTskTBkskB_{k+1}=B_{k}+\frac{\mathbf{y}_{k} \mathbf{y}_{k}^{\mathrm{T}}}{\mathbf{y}_{k}^{\mathrm{T}} \mathbf{s}_{k}}-\frac{B_{k} \mathbf{s}_{k} \mathbf{s}_{k}^{\mathrm{T}} B_{k}^{\mathrm{T}}}{\mathbf{s}_{k}^{\mathrm{T}} B_{k} \mathbf{s}_{k}} Bk+1=Bk+ykTskykykT−skTBkskBkskskTBkT
Limited memory Broyden–Fletcher–Goldfarb–Shanno (L-BFGS) 是对BFGS的一种优化，BFGS需要存储n∗nn* nn∗n的方阵BkB_kBk来近似Hessian矩阵的逆矩阵，而L-BFGS只存储最近m（m约为10）个（yk,sky_k, s_kyk,sk）用于近似BkB_kBk. 因此，L-BFGS的空间复杂度是O(mn)O(mn)O(mn).

L-BFGS算法：
yk=gk+1−gkρk=1ykTskq=gkFor i=k−1,k−2,…,k−mαi=ρisi⊤qq=q−αiyiγk=sk−1yk−1yk−1⊤yk−1Hk0=γkIz=Hk0qFor i=k−m,k−m+1,…,k−1βi=ρiyi⊤zz=z+si(αi−βi)StopwithHkgk=z\begin{array}{l} y_{k}=g_{k+1}-g_{k} \quad \rho_{k}=\frac{1}{y_{k}^{\mathrm{T}} s_{k}}\\ q=g_{k} \\ \text { For } i=k-1, k-2, \ldots, k-m \\ \quad \alpha_{i}=\rho_{i} s_{i}^{\top} q \\ \quad q=q-\alpha_{i} y_{i} \\ \gamma_{k}=\frac{s_{k-1} y_{k-1}}{y_{k-1}^{\top} y_{k-1}} \\ H_{k}^{0}=\gamma_{k} I \\ z=H_{k}^{0} q \\ \text { For } i=k-m, k-m+1, \ldots, k-1 \\ \quad \beta_{i}=\rho_{i} y_{i}^{\top} z \\ \quad z=z+s_{i}\left(\alpha_{i}-\beta_{i}\right)\\ Stop \;with \; H_k g_k = z \end{array} yk=gk+1−gkρk=ykTsk1q=gk For i=k−1,k−2,…,k−mαi=ρisi⊤qq=q−αiyiγk=yk−1⊤yk−1sk−1yk−1Hk0=γkIz=Hk0q For i=k−m,k−m+1,…,k−1βi=ρiyi⊤zz=z+si(αi−βi)StopwithHkgk=z

3. 总结

（1）牛顿法是通过求解一阶导数为0，从而引入Hessian矩阵指导参数优化；拟牛顿法是对牛顿法的性能优化，思路是近似求解Hessian矩阵。

（2）牛顿法在什么时候只需要迭代一次就能求解，什么时候牛顿法不能适用？

对于正定二次函数，一步即可得最优解。

当初始点远离最优解时，HkH_kHk不一定是正定的，则牛顿方向不一定为下降方向，其收敛性不能保证。这说明恒取步长因子为1是不合适的，应该采用一维搜索（仅当步长因子αk{αk}αk收敛1时，牛顿法才是二阶收敛的），此时迭代公式是：
dk=−Hk−1gk,xk+1=xk+αkdk\mathbf{d}_{k}=-H_{k}^{-1} \mathbf{g}_{k}, \mathbf{x}_{k+1}=\mathbf{x}_{k}+\alpha_{k} \mathbf{d}_{k} dk=−Hk−1gk,xk+1=xk+αkdk
带步长因子的牛顿法是总体收敛的。

（3）牛顿法和最小二乘法的区别？

牛顿法和最小二乘法相比，牛顿法需要求解Hessian矩阵及确定步长，而最小二乘法不需要。牛顿法是迭代求解，最小二乘法是计算解析解。如果样本量不算很大，且存在解析解，最小二乘法比起牛顿法要有优势，计算速度很快。但是如果样本量很大，用最小二乘法由于需要求一个超级大的逆矩阵，这时就很难或者很慢才能求解解析解了，使用迭代的牛顿法比较有优势。

（4）牛顿法和梯度下降法的区别，以及为什么深度学习不使用牛顿法进行参数优化？

两者都是迭代求解，不过梯度下降法是梯度求解，而牛顿法/拟牛顿法是用二阶的海森矩阵的逆矩阵或伪逆矩阵求解。相对而言，使用牛顿法/拟牛顿法收敛更快。但是每次迭代的时间比梯度下降法长。

深度学习一般不使用牛顿法的原因:

原因一：牛顿法需要用到Hessian矩阵，这难以求解，因为很难写出深度神经网络拟合函数的表达式。

原因二：即使可以得到梯度和Hessian矩阵，当输入向量的维度nnn较大时，Hessian矩阵的大小是n∗nn*nn∗n，所需要的内存非常大(L-BFGS会有所改善）。

原因三：在高维非凸优化问题中，鞍点相对于局部最小值的数量非常多，而且鞍点处的损失值相对于局部最小值处也比较大。而二阶优化算法是寻找梯度为0的点，所以很容易陷入鞍点。

最近开通了个公众号，主要分享推荐系统，风控等算法相关的内容，感兴趣的伙伴可以关注下。

公众号相关的学习资料会上传到QQ群596506387，欢迎关注。

参考：

Eureka 梯度下降法、牛顿法和拟牛顿法；
TangowL 数学优化入门：梯度下降法、牛顿法、共轭梯度法;
AM 221：Advanced Optimization ;
最优化方法（III）(推荐阅读）；
Quadratic Convergence of Newton’s Method Michael Overton, Numerical Computing, Spring 2017;
ysh329 Newton method;
wiki Limited-memory_BFGS.