无约束最优化问题的解析法主要有:最速下降法、牛顿法、共轭梯度法(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\}\) 是有穷点列时,其最后一个点是最优化模型问题的最优解。迭代规则由迭代公式决定,迭代公式的基本表示形式如下:
式中, \(\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$$转步骤2。
步骤中隐含了条件:函数\(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\)附近作二阶泰勒展开,有
为了计算\(\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\)为近似最优解;否则,转步骤3。
迭代计算。取迭代方向$$d_k=-[\nabla2f(x_k)]\nabla f(x_k)$$令$$x_{k+1}=x_k+d_k \ k \leftarrow k+1$$转步骤2。
牛顿法要求初始点在最优点附近(泰勒展开的前提就是在邻域内),否则可能不收敛。为了使得牛顿法能够全局收敛,提出了阻尼牛顿法(Damped Newton Method)。阻尼牛顿法的改进在于每次的搜索步长不固定为1,而是通过一维搜索来确定步长。其步骤如下:
- 设置初始值。给定迭代初值\(x_0\in R^n\),\(\epsilon>0\),令\(k \leftarrow 0\)。
- 检查终止条件。如果\(||\nabla f(x)||<\epsilon\),迭代终止,\(x_k\)为近似最优解;否则,转步骤3。
- 取迭代方向$$d_k=-[\nabla2f(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$$转步骤2。
可以证明,当\(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}xTQx+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\)。这里直接给出三个消去\(Q\)后\(\alpha_k\)的表达式:
(1)FR公式(Fletcher-Reeves)
(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\)次迭代就重新初始化可以减少计算误差的积累,有利于快速收敛。
共轭梯度法与牛顿法类似,具有二阶收敛速度。