牛顿法、拟牛顿法原理

非线性方程与其最优化方法

非线性方程指方程的因变量与自变量之间的关系不是线性关系的方程，比如平方关系、对数关系、指数关系、三角函数关系等。对于此类方程，求解n元实函数f在整个n维向量空间Rn上的最优值点往往很难得到精确解，经常需要求近似解问题。
求解该类方程的最优化问题的方法大多是逐次一维搜索的迭代算法，基本思想是在一个近似点处选定一个有利于搜索方向，沿这个方向进行一维搜索，得到新的近似点。如此反复迭代，直到满足预定的精度要求为止。
常用的迭代算法有:

• 梯度法:只使用当前点一阶导数决定搜索方向，又称最速下降法，是早期的解析法，收敛速度较慢。
• 牛顿法:利用函数在当前点的一阶导数和二阶导数寻找搜索方向。从本质上看，牛顿法是二阶收敛，梯度下降是一阶收敛，所以牛顿法更快。
• 拟牛顿法:在牛顿法的迭代中，需要计算海赛矩阵的逆矩阵H−1，这一计算比较复杂，我们可以考虑用一个正定矩阵来近似代替。这就是拟牛顿法的基本思想。
• 高斯牛顿法:使用泰勒级数展开式去近似地代替非线性回归模型，然后通过多次迭代，多次修正回归系数，使回归系数不断逼近非线性回归模型的最佳回归系数，最后使原模型的残差平方和达到最小。
• 共轭梯度法: 介于最速下降法与牛顿法之间的一个方法，它的每一个搜索方向是互相共轭的，而这些搜索方向仅仅是负梯度方向与上一次迭代的搜索方向的组合，因此，存储量少，计算方便。它仅需利用一阶导数信息，但克服了最速下降法收敛慢的缺点，又避免了牛顿法需要存储和计算Hesse矩阵并求逆的缺点，速度介于最速下降法与牛顿法之间。
• 直接法:不涉及导数,只使用函数值。有交替方向法（又称坐标轮换法）、模式搜索法、旋转方向法、鲍威尔共轭方向法和单纯形加速法等。

梯度法与牛顿法的比较:
梯度下降法每次选择当前位置梯度最大的方向，然后更新；牛顿法在当前位置选择方向时，不仅考虑是否是当前位置梯度最大的方向，还考虑该方向上继续更新时梯度是否会变得更大。因此牛顿法比梯度下降法看得更远一点，能更快地走到最底部。但是因为牛顿法需要计算二阶导数，计算量更大，实际使用时不如梯度下降法计算速度快。

牛顿法原理

牛顿法的基本思想是:在当前位置x0对f（x）做二阶泰勒展开，进而找到下一个位置x的估计值。
$$f(x)=f\left(x_0\right)+f^{\prime }\left(x_0\right)\left(x-x_0\right)+\frac{1}{2}{f}^{\prime \prime }\left(x_0\right){\left(x-x_0\right)}^2$$
对首式左右两边求导，得
$$f^{\prime }(x)=f^{\prime }\left(x_0\right)+f^{\prime \prime }\left(x_0\right)\left(x-x_0\right)$$
由于我们要求极值，此时有
$$f^{\prime }(x)=0$$
于是上式变为
$$f^{\prime }\left(x_0\right)+f^{\prime \prime }\left(x_0\right)\left(x-x_0\right)=0$$
于是可得求下一步位置x的公式
$$x=x_0-\frac{f^{\prime }\left(x_0\right)}{f^{\prime \prime }\left(x_0\right)}$$
写成迭代公式形式为
$$x_{k+1}=x_k-\frac{f^{\prime }\left(x_k\right)}{f^{\prime \prime }\left(x_k\right)},k=0,1,\cdots$$
上面是对单变量函数的牛顿法，对于多变量函数，二阶泰勒展开公式如下:
$$f(x)=f(x_0)+\nabla f(x_0)\cdot (x-x_0)+\frac{1}{2}\cdot (x-x_0{)}^T\cdot {\nabla }^{2}f(x_0)\cdot (x-x_0)$$
其中▽f是f的梯度向量，▽2f是f的海森矩阵。即:
∇ f = [ ∂ f ∂ x 1 ∂ f ∂ x 2 ⋮ ∂ f ∂ x N ] \nabla f=\left[

$$\begin{array}{c}{\frac{\partial f}{\partial x_{1}}} \\\\ {\frac{\partial f}{\partial x_{2}}} \\\\ {\vdots} \\\\ {\frac{\partial f}{\partial x_{N}}}\end{array}$$ \right] ∇f=⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎡​∂x1​∂f​∂x2​∂f​⋮∂xN​∂f​​⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎤​
$${\nabla }^{2}f=\left[\begin{array}{cccc}\frac{{\partial }^{2}f}{\partial x_1^2}& \frac{{\partial }^{2}f}{\partial x_1\partial x_2}& \dots & \frac{{\partial }^{2}f}{\partial x_1\partial x_N}\\ & & & \\ \frac{{\partial }^{2}f}{\partial x_2\partial x_1}& \frac{{\partial }^{2}f}{\partial x_2^2}& \dots & \frac{{\partial }^{2}f}{\partial x_2\partial x_N}\\ & & & \\ \frac{{\partial }^{2}f}{\partial x_N\partial x_1}& \frac{{\partial }^{2}f}{\partial x_N\partial x_2}& \cdots & \frac{{\partial }^{2}f}{\partial x_N^2}\end{array}\right]$$
海森矩阵其实就是多变量函数的所有二阶导数组成的NXN矩阵（假设有N个变量）。若f的混合偏导数可以交换次序（这个混合偏导数满足在（x0，y0）点连续即可），则海森矩阵是一个对称矩阵。
下面将▽f记为g，▽2f记为H。
同样对上面的二阶泰勒展开式两边求导，得:
$$\nabla f(x)=g_0+H_0\cdot (x-x_0)$$
又极值点处有
$$\nabla f(x)=0$$
于是可得
$$x=x_0-H_0^{-1}\cdot g_0$$
写成迭代公式形式为
$$x_{k+1}=x_k-H_k^{-1}\cdot g_k,k=0,1,\cdots$$
这就是牛顿迭代法，其搜索方向为
$$d_k=-H_k^{-1}\cdot g_k$$
又称为牛顿方向。

拟牛顿法原理

牛顿法需要计算二阶偏导数，计算量较大。而且有时目标函数的海森矩阵无法保持正定，此时牛顿法会失效。拟牛顿法的思想就是不使用海森矩阵，而是构造一个近似海森矩阵（或其逆矩阵）的正定对称阵来代替，在“拟牛顿”的条件下优化目标函数。

拟牛顿条件

先将函数f（x）在xk+1处二阶泰勒展开:
$$f(x)\approx f(x_{k+1})+\nabla f(x_{k+1})\cdot (x-x_{k+1})+\frac{1}{2}\cdot (x-x_{k+1}{)}^T\cdot {\nabla }^{2}f(x_{k+1})\cdot (x-x_{k+1})$$
两边求一阶偏导数（这里是一个一阶偏导数矩阵），得:
$$\nabla f(x)\approx \nabla f(x_{k+1})+H_{k+1}\cdot (x-x_{k+1})$$
取x=xk，得
$$g_{k+1}-g_k\approx H_{k+1}\cdot (x_{k+1}-x_k)$$
令
$$s_k=x_{k+1}-x_k,y_k=g_{k+1}-g_k$$
则上式可写成
$$y_k\approx H_{k+1}\cdot s_k或s_k\approx H_{k+1}^{-1}\cdot y_k$$
这就是拟牛顿条件。
也可以写成
$$y_k=B_{k+1}\cdot s_k或s_k=D_{k+1}\cdot y_k$$
Bk+1和Dk+1分别代表海森矩阵和海森矩阵的逆矩阵。
在满足此条件的基础上如何构造近似海塞矩阵呢？主要有两种方法:DFP算法和BFGS算法。

DFP算法

该算法的思想是通过迭代方法，逐步近似拟合海森矩阵的逆矩阵。迭代公式如下:
$$D_{k+1}=D_k+\Delta D_k,k=0,1,2,\cdots$$
其中D0通常取单位矩阵I。
这个迭代公式的关键是每一步的校正矩阵ΔDK的构造。我们采用待定法，先将ΔDK待定为下面的形式:
$$\Delta D_k=\alpha u{u}^T+\beta v{v}^T$$
α和β为待定系数，u和v为待定向量。uuT和vvT均为对称矩阵，因此可以保证ΔDK也是对称矩阵。
将待定公式代入迭代公式，并结合上面拟牛顿条件最后的sk和yk的表达式，可得:
$$s_k=D_k{y}_k+\alpha u{u}^T{y}_k+\beta v{v}^T{y}_k$$
继续将上式化为
s k = D k y k + u ( α u T y k ) + v ( β v T y k ) = D k y k + ( α u T y k ) u + ( β v T y k ) v

$$\begin{aligned} s_{k} &=D_{k} y_{k}+u(\alpha u^{T} y_{k})+v(\beta v^{T} y_{k}) \\\\ &=D_{k} y_{k}+(\alpha u^{T} y_{k})u+(\beta v^{T} y_{k})v \end{aligned}$$ sk​​=Dk​yk​+u(αuTyk​)+v(βvTyk​)=Dk​yk​+(αuTyk​)u+(βvTyk​)v​
上面之所以可以这样交换，是因为
$$\alpha u^T{y}_k，\beta v^T{y}_k$$
是两个数，我们不妨将其设为
$$\alpha u^T{y}_k=1,\beta v^T{y}_k=-1$$
即
$$\alpha =\frac{1}{u^T{y}_k},\beta =-\frac{1}{v^T{y}_k}$$
将上式代入sk转化后的公式，得
$$u-v=s_k-D_k{y}_k$$
我们不妨直接取
$$u=s_k,v=D_k{y}_k$$
代回上面α和β的表达式，得:
$$\alpha =\frac{1}{s_k^T{y}_k},\beta =-\frac{1}{(D_k{y}_k{)}^T{y}_k}=-\frac{1}{y_k^T{D}_k{y}_k}$$
将上两式代回ΔDK的表达式，得
$$\Delta D_k=\frac{s_k{s}_k^T}{s_k^T{y}_k}-\frac{D_k{y}_k{y}_k^T{D}_k}{y_k^T{D}_k{y}_k}$$
这样我们就可以使用迭代公式来进行拟牛顿法的计算了。
DFP算法步骤:

• 给定初值x0和阈值e，令B0=I，k=0；
• 每一次迭代时我们计算出搜索方向dk:
  $$d_k=-D_k{g}_k$$
• 由下式得到步长，然后计算出sk和xk+1；
  $${\lambda }_k=\arg \underset{\lambda \in R}{\min }f(x_k+\lambda d_k)$$
  $$s_k={\lambda }_k{d}_k,x_{k+1}=x_k+s_k$$
• 然后计算gk+1，若其值小于设定阈值，则停止计算；
  $$g_{k+1}=g(x_{k+1})$$
• 否则计算yk；
  $$y_k=g_{k+1}-g_k$$
• 然后计算Dk+1；
  $$D_{k+1}=D_k+\frac{s_k{s}_k^T}{s_k^T{y}_k}-\frac{D_k{y}_k{y}_k^T{D}_k}{y_k^T{D}_k{y}_k}$$
  令k=k+1，重新开始计算搜索方向dk。

BFGS算法

BFGS算法是构造近似海森矩阵（而不是构造海森矩阵的逆矩阵），同样，我们使用迭代的方式来逐步逼近。我们使用B来表示海森矩阵的近似矩阵，而在DFP算法中我们是直接使用D来构造近似海森矩阵的逆矩阵。
BFGS算法的迭代公式:
$$B_{k+1}=B_k+\Delta B_k,k=0,1,2,\cdots$$
与DFP一样，我们设
$$\Delta B_k=\alpha u{u}^T+\beta v{v}^T$$
于是有
$$y_k=B_k{s}_k+(\alpha u^T{s}_k)u+(\beta v^T{s}_k)v$$
也令
$$\alpha u^T{s}_k=1,\beta v^T{s}_k=-1$$
$$u=y_k,v=B_k{s}_k$$
可得
$$\alpha =\frac{1}{y_k^T{s}_k},\beta =-\frac{1}{s_k^T{B}_k{s}_k}$$
于是ΔBk为
$$\Delta B_k=\frac{y_k{y}_k^T}{y_k^T{s}_k}-\frac{B_k{s}_k{s}_k^T{B}_k}{s_k^T{B}_k{s}_k}$$
该算法的流程和DFP算法完全一样，只是将迭代公式换为BFGS算法的迭代公式，搜索方向换为:
$$d_k=-B_k^{-1}\cdot g_k$$
计算Bk+1公式为:
$$B_{k+1}=B_k+\frac{y_k{y}_k^T}{y_k^T{s}_k}-\frac{B_k{s}_k{s}_k^T{B}_k}{s_k^T{B}_k{s}_k}$$
其他步骤与DFP算法完全一样，具体过程可以参考上面的DFP算法。