非线性规划——无约束求解方法（三）-526互联

无约束最优化问题的解析法主要有：最速下降法、牛顿法、共轭梯度法（DFP法）和变尺度法（变度量法）。对于特殊的最小二乘问题，有最小二乘法。这些方法各有千秋，除了最小二乘法，后面的方法都针对前面方法的某个问题做了改进。这些方法的核心就是研究如何确定每一步迭代的方向和步长。

一、无约束最优化问题

最优化问题的一般形式为：

\[\begin{aligned} & \min f(x) \\ & \text { s.t.x } \in X \end{aligned} \]

其中 $x$ 为决策变量，$f(x)$ 是目标函数，$X$ 为约束集或可行域。特别地，如果 $X=R^n$ ，则最优化问题成为无约束最优化问题。
最优化方法通常采用迭代法求它的最优解，其基本思想是：给定一个初始点 $x_0$ ，按照某一迭代规则产品一个点列 $\left\{x_n\right\}$ ，使得当 $\left\{x_n\right\}$ 是有穷点列时，其最后一个点是最优化模型问题的最优解。迭代规则由迭代公式决定，迭代公式的基本表示形式如下:

\[x_{k+1}=x_k+\alpha_k d_k \]

式中， $\alpha_k$ 为步长因子， $d_k$ 为搜索方向。在最优化算法中，搜索方向 $d_k$ 是 $f$ 在 $x_k$ 点处的下降方向, 即:

\[f\left(x_k+\alpha_k d_k\right)<f\left(x_k\right) \]

最优化方法的基本结构如下:

给定初始点 $x_0$ ；

确定搜索方向 $d_k$ ，即按照一定规则，构造 $f$ 在 $x_k$ 点处的下降方向作为搜索方向；

确定步长因子 $\alpha_k$ ，使目标函数有某种意义的下降；

二、最速下降法

最速下降法(Gradient Descent Method)的基本思想是：因为连续函数沿着负梯度方向的下降速度是最快的（这一结论由梯度和方向导数的定义可以推出），所以每次迭代我们都从当前点出发，沿着负梯度方向前进一个最优步长，可以期望能较快逼近函数的极值。

最速下降法仅有三个步骤：
设置初始值。设置迭代起点$x_0\in R^n$,允许误差$\epsilon>0$和迭代变量初值$k\leftarrow0$。
检查终止条件。如果$||\nabla f(x_k)||<\epsilon$，停止迭代输出$x_k$作为近似最优解；否则转步骤3。
迭代，通过一维搜索求下一个迭代点。取搜索方向为负梯度方向$d_k=-\nabla f(x)$，求$\lambda_k$使得$$f(x_k+\lambda_k d_k)=\min_{\lambda\geq0} f(x_k+\lambda d_k)$$再令$$x_{k+1}=x_k+\lambda d_k$$转步骤２。

步骤中隐含了条件：函数$f(x)$必须可微，也就是说函数$f(x)$的梯度必须存在。这个步骤描述是最速下降发最抽象的形式。其中最关键的步骤是求解问题
\begin{equation}
\min_{\lambda\geq0}f(x_k+\lambda d_k)
\end{equation}
给出它的最优性必要条件
\begin{equation}
\frac{d f(x_k+\lambda d_k)}{d\lambda}=\nabla f(x_k+\lambda d_k)^T d_k=0
\end{equation}
当$f(x)$的形式确定，我们可以通过求解这个一元方程来获得迭代步长$\lambda$。当此方程形式复杂，解析解不存在，我们就需要使用“一维搜索”来求解$\lambda$了。一维搜索是一些数值方法，有0.618法、Fibonacci法、抛物线法等等，这里不详细解释了。

我们以一个简单的二次规划命题为例进行说，即

\[f(x)=\frac{1}{2} x^T Q x+c^T x \]

其中 $Q$ 是对称正定阵。对于目标函数 $f(x)$ ，记当前的迭代点为 $x_k ， f(x)$ 在 $x_k$ 处的一阶导数为 $\nabla f_k=Q x+c$ 。于是在 $x_k$ 处的前进方向即为 $-\nabla f_k$ 。要确定前进的步长，即要是的下面这个以 $\alpha$ 为决策变量的优化命题取得最小值:

\[f\left(x-\alpha \nabla f_k\right)=\frac{1}{2}\left(x-\alpha \nabla f_k\right)^T Q\left(x-\alpha \nabla f_k\right)+c^T\left(x-\alpha \nabla f_k\right) \]

令上式对 $\alpha$ 的导数为 0 即可求得对应步长为:

\[\alpha_k=\frac{\nabla f_k^T \nabla f_k}{\nabla f_k^T Q \nabla f_k} \]

因此，每一步更新迭代点的过程为:

\[x_{k+1}=x_k-\frac{\nabla f_k^T \nabla f_k}{\nabla f_k^T Q \nabla f_k} \nabla f_k \]

因此，最速下降法的迭代轨迹可以用下图描述。

在实际使用中，为了简便，也可以使用一个预定义的常数而不用一维搜索来确定步长$\lambda$。这时步长的选择往往根据经验或者通过试算来确定。步长过小则收敛慢，步长过大可能震荡而不收敛。最速下降法是最基本的迭代优化方法。它最简单，最基础，通常是收敛的。可以证明，最速下降法是一阶收敛的，往往需要多次迭代才能接近问题最优解。这是它的不足。

三、牛顿法

牛顿法
牛顿法是从函数的二阶泰勒展开式推导而来，其思想就是利用目标函数二阶近似的解去逼近目标函数的解。将函数$f$在迭代点$x_k$附近作二阶泰勒展开，有

\[f(x)\approx \phi(x)=f(x_k)+\nabla f(x_k)^T(x-x_k)+\frac{1}{2}(x-x_k)^T\nabla^2f(x_k)(x-x_k) \]

为了计算$\phi(x)$的极小值，令它的梯度为零，即

\[\nabla\phi(x)=0 \]

有

\[\nabla f(x_k)+\nabla^2f(x_k)(x-x_k)=0 \]

从而推出

\[x = x_k - [\nabla^2f(x_k)]^{-1}\nabla f(x_k) \]

我们将这里的$x$作为第$k+1$次迭代的估值，就有了迭代公式

\[x_{k+1}=x_k-[\nabla^2(fx_k)]^{-1}\nabla f(x_k) \]

牛顿法(Newton Method)利用了函数的二阶信息，即黑塞矩阵，来加速迭代收敛。具有二阶收敛速度是它的显著优势。牛顿法的步骤为：

设置初始值。给定迭代初值$x_0\in R^n$,$\epsilon>0$,令$k \leftarrow 0$。
检查终止条件。如果$||\nabla f(x)||<\epsilon$,迭代终止，$x_k$为近似最优解；否则，转步骤３。
迭代计算。取迭代方向$$d_k=-[\nabla^2f(x_k)]\nabla f(x_k)$$令$$x_{k+1}=x_k+d_k \ k \leftarrow k+1$$转步骤２。
牛顿法要求初始点在最优点附近(泰勒展开的前提就是在邻域内)，否则可能不收敛。为了使得牛顿法能够全局收敛，提出了阻尼牛顿法（Damped Newton Method）。阻尼牛顿法的改进在于每次的搜索步长不固定为１，而是通过一维搜索来确定步长。其步骤如下：

设置初始值。给定迭代初值$x_0\in R^n$,$\epsilon>0$,令$k \leftarrow 0$。
检查终止条件。如果$||\nabla f(x)||<\epsilon$,迭代终止，$x_k$为近似最优解；否则，转步骤３。
取迭代方向$$d_k=-[\nabla^2f(x_k)]\nabla f(x_k)$$
进行一维搜索确定步长$\lambda_k$使得$$f(x_k+\lambda_k d_k)=\min_{\lambda\geq0}f(x_k+\lambda d_k)$$令$$x_{k+1}=x_k+\lambda d_k\k\leftarrow k+1$$转步骤２。

可以证明，当$f(x)$具有二阶连续偏导数且$\nabla^2f(x)$正定，阻尼牛顿法是全局收敛的。

四、共轭梯度法

最速下降法收敛速度慢，牛顿法虽然收敛速度快但是需要计算黑塞矩阵的逆矩阵，计算复杂度比较高。有的时候我们希望有一种收敛快且不用计算二阶导数的方法，共轭梯度法（Conjugate Gradient Method）应运而生了。当然，计算量减小的代价是理论推导更加繁琐。

要了解共轭梯度法，首先要了解共轭的概念。设有两个$n$维向量$d_1$、$d_2$和一个$n$阶正定矩阵$Q$，如果它们满足

\[d_1^T Q d_2=0 \]

则称向量$d_1$、$d_2$关于$Q$共轭。当$Q$为单位阵时，共轭就成了正交。可将共轭看成是正交概念的推广。共轭的概念还可以推广到含有m个向量的向量组，即

\[d_i^TQd_j=0,\forall i,j=1,...m,i\neq j \]

我们用正定二次函数$$f(x)=\frac{1}{2}x^TQx+bTx+c$$来推导，最后推广到一般问题。注意这里的$Q\in R^{n\times n}$为正定矩阵，$b\in R^n$，$c\in R$。

我们从$x_0$出发，沿着某个下降方向$d_0$作一维搜索得到下一个迭代点$x_1$，那么有

\[f(x_1)=f(x_0+\lambda_0 d_0)=\min_{\lambda\geq0}f(x_0+\lambda d_0) \]

它的最优性条件是对$\lambda$求导，于是有

\[\nabla f(x_0+\lambda_0d_0)^Td_0=0 \]

即

\[\nabla f(x_1)^Td_0=0 \]

新的迭代方向$d_1$如果取负梯度方向，就成了最速下降法。负梯度方向虽然是当前函数值下降最快的方向，却未必是直指函数最优点的方向。设函数最优点为$\bar{x}$，我们希望搜索方向$d_1$能直接指向最优值，即

\[\bar{x}=x_1+\lambda_1 d_1 \]

而作为最优点，$\bar{x}$应满足

\[\begin{align} 0=\nabla f(\bar{x})&=Q\bar{x}+b \\ &=Q(x_1+\lambda_1 d_1)+b \\ &=(Qx_1+b)+\lambda_1 Q d_1 \\ &=\nabla f(x_1)+\lambda_1 Q d_1 \end{align} \]

两边左乘$d_0^T$并注意到$d_0^T \nabla f(x_1)=0$就有

\[d_0^TQd_1=0 \]

上式表明，搜索方向$d_1$与$d_0$是关于$Q$共轭的。

可以证明，如果能构造出与$Q$共轭的向量组，至多通过$n$次迭代就可以求得正定二次函数最优化问题的最优解。由此产生一类称为共轭方向法的方法，当我们利用梯度信息来构造与$Q$共轭的向量组，就产生了共轭梯度法。

下面推导如何利用梯度构造与$Q$共轭的向量组。首先初始搜索方向取为负梯度方向，即

\[d_0=-\nabla f(x_0) \]

新的搜索方向由负梯度方向和上一次搜索方向的线型组合产生，即

\[d_{k+1}=-\nabla f(x_{k+1})+\alpha_k d_k, k=0,1,...,n-2 \]

其中$\alpha_k$待定。我们利用$d_{k+1}$与$d_k$关于$Q$共轭作为约束条件，有

\[\begin{align} 0&=d_{k+1}^TQd_k\\ &=(-\nabla f(x_{k+1})^T+\alpha_k d_k^T)Qd_k\\ &=-\nabla f(x_{k+1})^TQd_k+\alpha_k d_k^TQd_k \end{align} \]

从而解得

\[\alpha_k=\frac{\nabla f(x_{k+1})^TQd_k}{d_k^TQd_k} \]

总结一下，我们得到了n个搜索方向的递推公式如下：

\[\begin{align} d_0 & =-\nabla f(x_0) \\ d_{k+1}&=-\nabla f(x_{k+1})+\alpha_k d_k \\ \alpha_k &=-\frac{\nabla f(x_{k+1})^TQd_k}{d_k^TQd_k} \end{align} \]

可以证明，上式得到的${d_k}$是一组关于$Ｑ$共轭的方向向量。

将此公式向一般无约束最优化问题推广，必须消去$Q$。这里直接给出三个消去$Q$后$\alpha_k$的表达式：
(1)FR公式（Fletcher-Reeves）

\[\alpha_k=\frac{||\nabla f(x_{k+1})||^2}{||\nabla f(x_k)||^2} \]

(2)DM公式（Dixon-Myers）

\[\alpha_k=-\frac{||\nabla f(x_{k+1})||^2}{d_k^T\nabla f(x_k)} \]

(3)PRP公式（Polak-Ribiere-Polyak）

\[\begin{align} \alpha_k&=\frac{\nabla f(x_{k+1})^T[\nabla f(x_{k+1})-\nabla f(x_k)]}{||\nabla f(x_k)||^2} \\ &=\frac{||\nabla f(x_{k+1})||^2 - \nabla f(x_{k+1})^T \nabla f(x_k)}{||\nabla f(x_k)||^2} \end{align} \]

当问题为正定二次函数优化问题时，这三个公式是等价的。对非二次函数最优化问题，它们将产生不同的搜索方向。

共轭梯度法的步骤总结如下：

设置初始值。迭代初始点$x_0$和允许误差$\epsilon$。
检查终止条件。如果$||\nabla f(x_0)||<\epsilon$，迭代终止并输出$x_0$；否则转步骤3。
构造初始搜索方向。$$ d_0=-\nabla f(x_0)$$　并令$k=0$
进行一维搜索。求$\lambda_k$使得 $$ f(x_k+\lambda_k d_k)=\min_{\lambda\geq0}f(x_k+\lambda d_k)$$并令$$ x_{k+1}=x_k+\lambda_k d_k$$
检查终止条件。如果$||\nabla f(x_{k+1})||<\epsilon$，迭代终止并输出$x_{k+1}$；否则转步骤6。
检查迭代次数。如果$k+1=n$，令$$ x_0=x_n$$转步骤3；否则转步骤7。
构造共轭方向。令$$ \begin{align} d_{k+1}&=-\nabla f(x_{k+1})+\alpha_k d_k \
\alpha_k &=\frac{||\nabla f(x_{k+1})||^2}{||\nabla f(x_k)||^2} \end{align}$$
这里$\alpha_k$的计算取了FR公式，也可以选取另外两个公式。再令$k=k+1$，转步骤4。
共轭梯度法在$n$次迭代后将无法产生新的共轭方向（因为$R^n$中的共轭方向至多只能有$n$个），故将第$n$次迭代点$x_n$作为新的起点，重新进行一轮共轭梯度迭代。共轭梯度法比最速下降法的收敛条件更弱。当$f$具有一阶连续偏导数且有极值点时，共轭梯度法就是收敛的。对于二次函数最优化，理论上共轭梯度法$n$次迭代即可收敛。在一定条件下，共轭梯度法具有二次收敛速度。

共轭梯度法的Matlab实现请参考我的另一篇文章：Matlab实现FR共轭梯度法。

五、变尺度法

最速下降法和阻尼牛顿法的迭代公式可以写成

\[\begin{aligned} x_{k+1}=x_k+\lambda_k d_k \\ d_k=-G_k \nabla f(x_k) \end{aligned} \]

当$G_k$取单位阵时是最速下降法，当$G_k=[\nabla^2 f(x_k)]^{-1}$时是阻尼牛顿法。事实上，我们把上式确定的迭代法统称为拟牛顿法(Quasi-Newton Method)。因此改进的思路可以让$G_k$近似$[\nabla^2 f(x_k)]^{-1}$，从而避免计算二阶导数和矩阵求逆，既简化计算又保持快速收敛性。

我们同样利用二阶泰勒展开来考虑如何近似。将$f$在点$x_{k+1}$处作二阶泰勒展开，有

\[f(x)\approx f(x_{k+1})+\nabla f(x_{k+1})^T(x-x_{k+1})+\frac{1}{2}(x-x_{k+1})^T\nabla^2 f(x_{k+1})(x-x_{k+1}) \]

对上式两边求梯度，有

\[\nabla f(x) \approx \nabla f(x_{k+1}) + \nabla^2f(x_{k+1})(x-x_{k+1}) \]

令$x=x_k$，有

\[\nabla f(x_k) \approx \nabla f(x_{k+1}) + \nabla^2f(x_{k+1})(x_k-x_{k+1}) \]

可得如下近似关系

\[[\nabla^2 f(x_{k+1})]^{-1}[\nabla f(x_{k+1}) - \nabla f(x_k)]=x_{k+1}-x_k \]

我们可以令$G_{k+1}$满足此关系式，即

\[G_{k+1}[\nabla f(x_{k+1}) - \nabla f(x_k)]=x_{k+1}-x_k \]

从而近似阻尼牛顿法。为了简化后续推导，记

\[\begin{aligned} \Delta x_k &=x_{k+1}-x_k \\ \Delta g_k &= \nabla f(x_{k+1}) - \nabla f(x_k) \\ \Delta G_k &= G_{k+1}-G_k \end{aligned} \]

则上式可以改写为

\[\Delta G_k \Delta g_k = \Delta x_k - G_k \Delta g_k \]

为了求出$\Delta G_k$的表达式，我们采用待定系数法，设

\[\Delta G_k = \Delta x_k p_k^T - G_k \Delta g_k q_k^T \]

其中$p_k, q_k \in R^n$为待定向量。两边右乘$\Delta g_k$有

\[\Delta G_k \Delta g_k = (p_k^T \Delta g_k) \Delta x_k + (q_k^T \Delta g_k) G_k \Delta g_k \]

与前面推导的近似式子对比，则

\[p_k^T \Delta g_k=1 \\ q_k^T \Delta g_k = 1 \]

考虑到$G_k$近似于黑塞矩阵的逆（为对称矩阵），也应该为对阵矩阵，可构造性地设

\[\begin{aligned} p_k&=\alpha_k \Delta x_k \ q_k&=\beta_kG_k \Delta g_k \end{aligned} \]

解得

\[\begin{aligned} \alpha_k&=\frac{1}{\Delta x_k^T\Delta g_k} \\ \beta_k&=\frac{1}{\Delta g_k G_k^T \Delta g_k} \end{aligned} \]

回代可得

\[\Delta G_k=\frac{\Delta x_k \Delta x_k^T}{\Delta x_k^T \Delta g_k}-\frac{G_k\Delta g_k \Delta g_k^T G_k}{\Delta g_k^T G_k \Delta g_k} \]

上述公式由Davidon提出，被Fletcher和Powell改进得到，称为DFP公式。应用DFP公式的拟牛顿法称为DFP法或者变尺度法（变度量法）。

最后总结一下变尺度法的步骤：

选取初始值。初始点$x_0 \in R^n$，允许误差 $\epsilon>0$和$G_0=I_n$为单位阵。
检查终止条件。如果$||\nabla f(x_0)||<\epsilon$，终止迭代并输出$x_0$；否则转步骤6。
构造初始DFP方向。取$d_0=-\nabla f(x_0)$，令$k=0$
一维搜索确定步长$\lambda_k$，使得$$ f(x_k+\lambda_k d_k)=\min_{\lambda\geq0}f(x_k+\lambda d_k)$$并令$$x_{k+1}=x_k+\lambda_k d_k$$
检查终止条件。如果$||\nabla f(x_{k+1})||<\epsilon$，终止迭代并输出$x_{k+1}$；否则转步骤6。
检查迭代次数。如果$k+1=n$，令$x_0 \leftarrow x_n$并转步骤(3)；否则转步骤7。
构造DFP方向。令$$\begin{aligned} G_{k+1}&=G_k+\frac{\Delta x_k \Delta x_k^T}{\Delta x_k^T \Delta g_k}-\frac{G_k\Delta g_k \Delta g_k^T G_k}{\Delta g_k^T G_k \Delta g_k} \ d_{k+1}&=-G_{k+1} \nabla f(x_k) \end{aligned}$$令$k \leftarrow k+1$，转步骤4。
这里每$n$次迭代就重新初始化可以减少计算误差的积累，有利于快速收敛。
共轭梯度法与牛顿法类似，具有二阶收敛速度。

参考文献