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

1. 牛顿法

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

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

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

1.1 求方程$f(x)=0$的根

选择一个接近函数$f(x)$零点的${\mathbf{x}}_0$, 牛顿法对函数进行一阶泰勒展开：
$$f(x)=f\left(x_0\right)+f^{\prime }\left(x_0\right)\left(x-x_0\right)$$

由$f(x)=0$，得到迭代方式：
$$x_{n+1}=x_n-f\left(x_n\right)/f^{\prime }\left(x_n\right)$$
迭代后求得方程的根$x^{\ast }$：
$$f\left(x^{\ast }\right)=0$$

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

假设函数$f$在开区间$(a,b)$是二阶可导的,并存在函数的根$x^{\ast }\in (a,b)$.定义牛顿迭代法：
$$x_{k+1}=x_k-\frac{f\left(x_k\right)}{f^{\prime }\left(x_k\right)},k=1,2,\dots$$
假设$k\to \infty$时，$x_k$收敛到$x^{\ast }$。若$f^{\prime }\left(x^{\ast }\right)\ne 0$，对于足够大的 $k$，有：
$$\left|x_{k+1}-x^{\ast }\right|\le M{\left|x_k-x^{\ast }\right|}^2\text{ if }M>\frac{\left|f^{\prime \prime }\left(x^{\ast }\right)\right|}{2\left|f^{\prime }\left(x^{\ast }\right)\right|}$$
于是，$x_k$是二次收敛到$x^{\ast }$。

证明：
假设$e_k=x_k-x^{\ast }$, 即 $x^{\ast }=x_k-e_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)$$
其中，${\xi }_k$介于$x_k$和$x^{\ast }$之间。
由于$f(x^{\ast })=0$,所以有
$$0=f\left(x_k\right)-\left(x_k-x^{\ast }\right)f^{\prime }\left(x_k\right)+\frac{{\left(e_k\right)}^2}{2}{f}^{\prime \prime }\left({\xi }_k\right)$$
由于函数$f$连续可导，且$f^{\prime }\left(x^{\ast }\right)\ne 0$，只要$x_k$和$x^{\ast }$足够接近，则有$f^{\prime }\left(x_k\right)\ne 0$。两边除于$f^{\prime }\left(x_k\right)$得到
$$0=\frac{f\left(x_k\right)}{f^{\prime }\left(x_k\right)}-\left(x_k-x^{\ast }\right)+\frac{{\left(e_k\right)}^2{f}^{\prime \prime }\left({\xi }_k\right)}{2f^{\prime }\left(x_k\right)}$$
根据牛顿迭代法的定义，可以得到下式：
$$x_{k+1}-x^{\ast }=\frac{{\left(e_k\right)}^2{f}^{\prime \prime }\left({\xi }_k\right)}{2f^{\prime }\left(x_k\right)}$$
所以有，
$$\left|x_{k+1}-x^{\ast }\right|\le \frac{\left|f^{\prime \prime }\left({\xi }_k\right)\right|}{2\left|f^{\prime }\left(x_k\right)\right|}{\left|x_k-x^{\ast }\right|}^2$$

一般地，$f^{\prime }\left(x_k\right)$会收敛到$f^{\prime }\left(x^{\ast }\right)$，由于${\xi }_k$介于$x_k$和$x^{\ast }$之间，因此，${\xi }_k$收敛到$x^{\ast }$，$f^{\prime \prime }\left({\xi }_k\right)$收敛到$f^{\prime \prime }\left(x^{\ast }\right)$。对应足够大的$k$，有
$$\left|x_{k+1}-x^{\ast }\right|\le M{\left|x_k-x^{\ast }\right|}^2\text{ if }M>\frac{\left|f^{\prime \prime }\left(x^{\ast }\right)\right|}{2\left|f^{\prime }\left(x^{\ast }\right)\right|}$$

证毕。

1.2 最优化（求取极值）

解决最优化问题$\underset{x\in {\mathbf{R}}^n}{\min }f(x)$的结构：

给定初始点$x_0$,

1. 确定搜索方向$d_k$，即依照一定规则构造$f$在$x_k$点处的下降方向为搜索方向;
2. 确定步长因子${\alpha }_k$，使目标函数值有某种意义下降;
3. 令$x_{k+1}=x_k+{\alpha }_k{d}_k$

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

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

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

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

$$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)$$
其中，$g_k=g\left(x_k\right)=\nabla f\left(x_k\right)$是$f(x)$在$x_k$的一阶导数值，$H\left(x_k\right)$是$f(x)$的海森矩阵：
$$H(x)={\left[\frac{{\partial }^{2}f}{\partial x_i\partial x_j}\right]}_{n\times n}$$

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

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

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

1.3 牛顿法最优化的示例

已知(X, y) 样本，X 为特征值，y为label。

建立线性回归模型对样本进行拟合：y_pred = Wx + b.

则，用平方差的一半作为损失（乘于0.5后，W一阶导的系数为1）， loss = $0.5\ast (y-y_{pred}{)}^2=0.5\ast (y-Wx-b{)}^2$

最小化损失，即 $Mi{n}_{W,b}(y-Wx-b{)}^2$.

根据上述最优化理论，求取，W,b 一阶导，二阶导，进行迭代更新即可.

注意：sklearn Python包安装为：
pip install scikit-learn

import numpy as np
from sklearn.linear_model import LinearRegression
class Newton(object):
	def __init__(self,epochs=50):
		self.W = None
		self.b = 0
		self.epochs = epochs
    
	def get_loss(self, X, y, W,b):
		"""
		计算损失 0.5*sum(y_pred-y)^2
		input: 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 loss
		
	def first_derivative(self,X,y):
		"""
		计算一阶导数g = (y_pred - y)*x
		input: X(2 dim np.array):特征
				y(1 dim np.array):标签
				W(2 dim np.array):线性回归模型权重矩阵
		output：损失函数值
		"""
		y_pred = np.dot(X,self.W) + self.b
		g = np.dot(X.T, np.array(y_pred - y))
		g_b = np.mean(y_pred-y)
		return g,g_b
		 
	def 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 = 1
		return H, H_b
		
	def 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 = 0
		for 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_b
			print("itration:{} ".format(epoch), "loss:{:.4f}".format(
			self.get_loss(X, y , self.W,self.b)))
		
    def predict(self, X):
    	"""
    	预估未知的样本 X, y = W * X + b.
    	"""
        return np.dot(X, 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)
	print("newton predict mse: ", 1 / len(y_test) * np.sum((y_test - y_pred)**2))
	reg = LinearRegression().fit(X_train, y_train)
	y_pred = reg.predict(X_test)
	print("lr predict mse: ", 1 / len(y_test) * np.sum((y_test - y_pred)**2))

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

itration:48  loss:8.5195
itration:49  loss:8.2219
newton predict mse:  0.3770617848430015
lr predict mse:  0.32062598186600233

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

itration:498  loss:4.9358
itration:499  loss:4.9358
newton predict mse:  0.3206259816847965
lr predict mse:  0.32062598186600233

2. 拟牛顿法

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

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

2.1 一般拟牛顿法

由
$$g(x)=\nabla f(x)\approx g_k+H_k\left(x-x_k\right)$$
令$x=x_k$,得：
$$g_{k-1}-g_k\approx H_k\left(x_{k-1}-x_k\right)$$
令${\mathbf{s}}_{k-1}={\mathbf{x}}_k-{\mathbf{x}}_{k-1},{\mathbf{y}}_{k-1}={\mathbf{g}}_k-{\mathbf{g}}_{k-1}$,得
$$H_k^{-1}{\mathbf{y}}_{k-1}\approx {\mathbf{s}}_{k-1}$$
或者记为：
$$H_{k+1}^{-1}{\mathbf{y}}_k\approx {\mathbf{s}}_k$$
对于二次函数$f$,上述关系式精确成立。

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

或者构建Hesse近似（如BFGS算法）
$$B_{k+1}{\mathbf{s}}_k\approx {\mathbf{y}}_k$$

这称为拟牛顿法条件 。

一般拟牛顿法：

1. 给定初始点${\mathbf{x}}_0\in R^n,H_0\in R^{n\times n},0\le \epsilon \le 1,k=0$;
2. 若$\Vert {\mathbf{g}}_k\Vert \le \epsilon$，则停止；否则，计算${\mathbf{d}}_k=-H_k{\mathbf{g}}_k$；
3. 沿方向$d_k$线性搜索求步长因子${\alpha }_k$，令$x_{k+1}=x_k+{\alpha }_k{d}_k$；
4. 校正$H_k$产生$H_{k+1}$,使得拟牛顿法条件满足。
5. $k=k+1$， 转2.

优点：

（1）只需要一阶导数；

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

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

2.2 DFP算法

设秩二校正为：
$$G_{k+1}=G_k+au{u}^T+bv{v}^T$$
若要拟牛顿法条件$G_{k+1}{y}_k=(G_k+au{u}^T+bv{v}^T)y_k=s_k$成立，对于$u,v$一个取法为：
$$u=s_k,v=G_k{y}_k$$
$$a{u}^T{y}_k=1,b{v}^T{y}_k=-1$$
可以得到，
$$a=\frac{1}{s_k^T{y}_k},b=-\frac{1}{y_k^T{G}_k{y}_k}$$
那么，
$$G_{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}$$

2.2 L-BFGS算法

BFGS算法推导和DFP算法差不多，直接给出公式：
$$B_{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}$$
Limited memory Broyden–Fletcher–Goldfarb–Shanno (L-BFGS) 是对BFGS的一种优化，BFGS需要存储$n\ast n$的方阵$B_k$来近似Hessian矩阵的逆矩阵，而L-BFGS只存储最近m（m约为10）个（$y_k,s_k$）用于近似$B_k$. 因此，L-BFGS的空间复杂度是$O(mn)$.

L-BFGS算法：
$$\begin{array}{l}{y}_k=g_{k+1}-g_k{\rho }_k=\frac{1}{y_k^{\mathrm{T}}{s}_k}\\ q=g_k\\ \text{ For }i=k-1,k-2,\dots ,k-m\\ {\alpha }_i={\rho }_i{s}_i^{\top }q\\ 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,\dots ,k-1\\ {\beta }_i={\rho }_i{y}_i^{\top }z\\ z=z+s_i\left({\alpha }_i-{\beta }_i\right)\\ Stopwith{H}_k{g}_k=z\end{array}$$

3. 总结

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

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

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

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

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

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

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

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

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

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

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

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

---

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

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

---

参考：

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