牛顿法和拟牛顿法

牛顿法与拟牛顿法是求解无约束最优化问题的常用方法，有收敛速度快的优点。

牛顿法是迭代算法，每一步需要求解目标函数的海森矩阵的逆矩阵，计算比较复杂。

拟牛顿法通过正定矩阵近似海塞矩阵的逆矩阵或海塞矩阵，简化了这一计算过程。

15.1 牛顿法

无约束最优化问题：

$$\min_{x \in R^n} f(x) \ (15.1)$$

其中$x^*$为目标函数的极小点。

假设f(x)具有二阶连续偏导数，若第k次迭代值为$x^{(k)}$，则可将$f(x)$在$x^{(k)}$附近进行二阶泰勒展开:

$$f(x) = f(x^{(k)})+g_k^{\top}(x-x^{(k)})+\frac{1}{2}(x-x^{(k)})^{\top}H(x^{(k)})(x-x^{(k)}) \ (15.2)$$

这里，$g_k=g(x^{(k)})=\nabla f(x^{(k)})$是f(x)的梯度向量在点$x^{(k)}$的值，$H(x^{(k)})$是$f(x)$的海森矩阵:

$$H(x) = \left[\frac{\partial^2 f}{\partial x_i \partial x_j}\right]_{n\times n} \ (15.3)$$

在点$x^{(k)}$的值。函数f(x)有极值的必要条件是在极值点处一阶导数为0，即梯度向量为0。特别是当$H(x^{(k)})$是正定矩阵时，函数$f(x)$的极值为极小值。

牛顿法利用极小点的必要条件:

$$\nabla f(x) = 0 \ (15.4)$$

每次迭代中，从点$x^{(k)}$开始，求目标函数的极小点，作为第k+1次迭代值$x^{(k+1)}$。具体地，假设$x^{(k+1)}$满足:

$$\nabla f(x^{(k+1)}) = 0 \ (15.5)$$

由式(15.2)有:

$$\nabla f(x) = g_k+H_k(x-x^{(k)}) \ (15.6)$$

其中$H_k = H(x^{(k)})$。这样，式(15.5)成为

$$g_k + H_k(x^{(k+1)}-x^{(k)}) = 0 \ (15.7)$$

因此，

$$x^{(k+1)}=x^{(k)}-H_k^{-1}g_k \ (15.8)$$

或者

$$x^{(k+1)} = x^{(k)} + p_k \ (15.9)$$

其中，

$$H_kp_k = -g_k \ (15.10)$$

用式(15.8)作为迭代公式的算法就是牛顿法。

算法15.1(牛顿法)

输入：目标函数f(x)，梯度$g(x)=\nabla f(x)$。海森矩阵H(x)，精度要求$\varepsilon$

输出：f(x)的极小点$x^*$

(1)取初始点$x^{(0)}$，置k=0

(2)计算$g_k=g(x^{(k)})$

(3)若$\|g_k\| < \varepsilon$，则停止计算，得近似解$x^*=x^{(k)}$

(4)计算$H_k=H(x^{(k)})$，并求$p_k$

$$H_kp_k = -g_k$$

(5)置$x^{(k+1)} = x^{(k)}+p_k$

(6)置k=k+1，转(2)

步骤(4)求$p_k$，$p_k=-H^{-1}_k g_k$，要求$H_k^{-1}$，计算比较复杂，所以有其他的改进方法。

15.2 拟牛顿法的思路

在牛顿法的迭代中，需要计算海塞矩阵的逆矩阵$H^{-1}$，这一计算比较复杂，考虑用一个n阶矩阵$G_k=G(x^{(k)}$来近似代替$H_k^{(-1)}=H^{(-1)}(x^{(k)}$。这就是拟牛顿法的基本想法。

在牛顿法迭代中，海森矩阵$H_k$需要满足以下条件：

在(15.6)中取$x=x^{(k+1)}$，即得

$$g_{k+1}-g_k = H_k(x^{(k+1)}-x^{(k)}) \ (15.11)$$

记$y_k=g_{k+1}-g_k$，$\delta_k = x^{(k+1)}-x^{(k)}$，则

$$y_k = H_k \delta_k \ (15.12)$$

或

$$H^{-1}_k y_k = \delta_k \ (15.13)$$

式(15.12)或式(15.13)称为拟牛顿条件。

如果$H_k$是正定的($H_k^{-1}$也是正定的)，那么可以保证牛顿法搜索方向$p_k$是下降方向。这是因为搜索方向是$p_k=-H_k^{-1}g_k$，由式(15.8)有：

$$x = x^{(k)}+\lambda p_k = x^{(k)}-\lambda H_k^{-1}g_k \ (15.14)$$

所以$f(x)$在$x^{(k)}$的泰勒展开式(15.2)可以近似写成：

$$f(x) = f(x^{(k)})-\lambda g_k^{\top}H_k^{-1}g_k \ (15.15)$$

因$H_k^{-1}$正定，故有$g_k^{\top}H_k^{-1}g_k>0$。当$\lambda$为一个充分小的正数时，总有$f(x)<f(x^{(k)})$，也就是说$p_k$是下降方向。

将$G_k$作为$H_k^{-1}$的近似或选择$B_k$作为$H_k$的近似的算法称为拟牛顿法。

首先，每次迭代矩阵$G_k$是正定的。同时，$G_k$满足下面的拟牛顿条件:

$$G_{k+1}y_k = \delta_k \ (15.16)$$

按照拟牛顿条件选择$G_k$作为$H_k^{-1}$的近似或选择$B_k$作为$H_k$的近似的算法称为拟牛顿法。

按照拟牛顿条件，在每次迭代中可以选择更新矩阵$G_{k+1}$:

$$G_{k+1} = G_k+\Delta G_k \ (15.17)$$

关于如何更新矩阵$G_k$，有多种实现方法：

15.2.1 DEP(Davidon-Flecher-Powell)算法

DFP算法选择$G_{k+1}$的方法是，假设每一步迭代中矩阵$G_{k+1}$是由$G_k$加上两个附加项构成的，即

$$G_{k+1}=G_k+P_k+Q_k \ (15.18)$$

其中$P_k$，$Q_k$是待定矩阵。这时，

$$G_{k+1}y_k = G_ky_k+P_ky_k+Q_ky_k \ (15.19)$$

为使$G_{k+1}$满足拟牛顿条件，可使$P_k$和$Q_k$满足：

$$P_ky_k = \delta_k \ (15.20)$$ $$Q_ky_k = -G_k y_k \ (15.21)$$

事实上，不难找出这样的$P_k$和$Q_k$，比如：

$$P_k = \frac{\delta_k \delta_k^{\top}}{\delta_k^{\top}y_k} \ (15.22)$$ $$Q_k = - \frac{G_ky_ky_k^{\top}G_k}{y_k^{\top}G_k y_k} \ (15.23)$$

这样既可得到矩阵$G_{k+1}$的迭代公式:

$$G_{k+1} = G_k+\frac{\delta_k \delta_k^{\top}}{\delta_k^{\top}y_k} - \frac{G_ky_ky_k^{\top}G_k}{y_k^{\top}G_k y_k} \ (15.24)$$

称为DFP算法。

可以证明，如果初始矩阵$G_0$是正定的，则迭代过程中的每个矩阵$G_k$都是正定的。

算法15.2(DFP算法)

输入：目标函数$f(x)$，梯度$g(x)=\nabla f(x)$，精度要求$\varepsilon$；

输出：$f(x)$的极小点$x^*$。

(1)选定初始点$x^{(0)}$，取$G_0$为正定对称矩阵，置$k=0$

(2)计算$g_k = g(x^{(k)})$。若$\|g_k\|<\varepsilon$，则停止计算，得近似解$x^* = x^{(k)}$；否则转(3)

(3)置$p_k = -G_kg_k$

(4)一维搜索:求$\lambda_k$使得

$$f(x^{(k)}+\lambda_k p_k) = \min_{\lambda \ge 0} f(x^{(k)}+\lambda p_k)$$

(5)置$x^{(k+1)}=x^{(k)}+\lambda_kp_k$

(6)计算$g_{k+1}=g(x^{(k+1)})$，若$\|g_{k+1}\|<\varepsilon$，则停止计算，得近似解$x^* = x^{(k+1)}$；否则，按式(15.24)算出$G_{k+1}$。

(7)置k=k+1，转(3)

15.2.2 BFGS(Broyden-Fletcher-Goldfarb-Shanno)算法

用$B_k$逼近海塞矩阵H。这时，相应的拟牛顿条件是:

$$B_{k+1}\delta_k = y_k \ (15.25)$$

可以用同样的方法得到另一迭代公式。

首先令：

$$B_{k+1} = B_k+P_k+Q_k \ (15.26)$$ $$B_{k+1}\delta_k = B_k \delta_k+P_k\delta_k+Q_k\delta_k \ (15.27)$$

考虑使$P_k$和$Q_k$满足:

$$P_k \delta_k = y_k \ (15.28)$$ $$Q_k \delta_k = -B_k \delta_k \ (15.29)$$

找出适合条件的$P_k$和$Q_k$，得到BFGS算法矩阵$B_{k+1}$的迭代公式:

$$B_{k+1} = B_k + \frac{y_k y_k^{\top}}{y_k^{\top}\delta_k}-\frac{B_k\delta_k \delta_k^{\top}B_k}{\delta_k^{\top}B_k \delta_k} \ (15.30)$$

可以证明，如果初始矩阵$B_0$是正定的，则迭代过程中的每个矩阵$B_k$都是正定的。

下面写出BFGS拟牛顿法

算法15.3(BFGS算法)

输入：目标函数$f(x)$，$g(x)=\nabla f(x)$，精度要求$\varepsilon$；

输出：$f(x)$的极小点$x^*$.

(1)选定初始点$x^{(0)}$，取$B_0$为正定对称矩阵，置k=0

(2)计算$g_k=g(x^{(k)})$。若$\|g_k\|<\varepsilon$，则停止计算，得近似解$x^* = x^{(k)}$；否则转(3)

(3)由$B_kp_k = -g_k$求得$p_k$

(4)一维搜索:求$\lambda_k$使得

$$f(x^{(k)}+\lambda_k p_k) = \min_{\lambda \ge 0} f(x^{(k)}+\lambda p_k)$$

(5)置$x^{(k+1)}=x^{(k)}+\lambda_k p_k$

(6)计算$g_{k+1}=g(x^{(k+1)})$，若$\|g_{k+1}\|<\varepsilon$，则停止计算，得近似解$x^* = x^{(k+1)}$；否则，按式(15.30)算出$B_{k+1}$

(7)置k=k+1，转(3)

15.2.3 Broyden类算法

从BFGS算法矩阵$B_k$的迭代公式(15.30)得到BFGS算法关于$G_k$的迭代公式。

Sherman-Morrison公式：假设A是n阶可逆矩阵，u,v是n维向量，且$A+uv^{\top}$也是可逆矩阵，则：

$$(A+uv^{\top})^{-1} = A^{-1}-\frac{A^{-1}uv^{\top}A^{-1}}{1+v^{\top}A^{-1}u}$$

事实上，若记$G_k = B_k^{-1}$，$G_{k+1}=B_{k+1}^{-1}$，那么对式(15.30)应用Sherman-Morrison公式即得：

$$G_{k+1} = \left( I-\frac{\delta_k y_k^{\top}}{\delta_k^{\top}y_k} \right)G_k \left(I-\frac{\delta_k y_k^{\top}}{\delta_k^{\top}y_k}\right)^{\top}+\frac{\delta_k y_k^{\top}}{\delta_k^{\top}y_k} \ (15.31)$$

称为BFGS算法关于$G_k$的迭代公式。

由DFP算法$G_k$的迭代公式(15.23)得到的$G_{k+1}$记作$G^{\mathrm{DFP}}$，由BFGS算法$G_k$的迭代公式(15.31)得到的$G_{k+1}$记作$G^{\mathrm{BFGS}}$，他们都满足都满足方程拟牛顿条件公式，所以他们的线性组合：

$$G_{k+1}=\alpha G^{\mathrm{DFP}}+(1-\alpha)G^{\mathrm{BFGS}} \ (15.32)$$

也满足拟牛顿条件式，而且是正定的。

其中$0\le \alpha \le 1$。这样就得到了一类拟牛顿法，称为Broyden类算法。