首页 > 其他分享 >非线性规划——无约束求解方法(三)

非线性规划——无约束求解方法(三)

时间:2023-06-08 23:45:21浏览次数:86  
标签:非线性 迭代 求解 nabla 无约束 Delta alpha 共轭 lambda

无约束最优化问题的解析法主要有:最速下降法、牛顿法、共轭梯度法(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$$转步骤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\)附近作二阶泰勒展开,有

\[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\)为近似最优解;否则,转步骤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,而是通过一维搜索来确定步长。其步骤如下:

  1. 设置初始值。给定迭代初值\(x_0\in R^n\),\(\epsilon>0\),令\(k \leftarrow 0\)。
  2. 检查终止条件。如果\(||\nabla f(x)||<\epsilon\),迭代终止,\(x_k\)为近似最优解;否则,转步骤3。
  3. 取迭代方向$$d_k=-[\nabla2f(x_k)]\nabla f(x_k)$$
  4. 进行一维搜索确定步长\(\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)

\[\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\)次迭代就重新初始化可以减少计算误差的积累,有利于快速收敛。
共轭梯度法与牛顿法类似,具有二阶收敛速度。

参考文献

  1. 机器学习之数学
  2. 无约束最优化方法学习笔记
  3. 最优化方法总结

标签:非线性,迭代,求解,nabla,无约束,Delta,alpha,共轭,lambda
From: https://www.cnblogs.com/haohai9309/p/17465178.html

相关文章

  • 非线性规划凸优化——凸函数、凸规划(二)
    凸规划是指若最优化问题的目标函数为凸函数,不等式约束函数也为凸函数,等式约束函数是仿射的。凸规划的可行域为凸集,因而凸规划的局部最优解就是它的全局最优解。当凸规划的目标函数为严格凸函数时,若存在最优解,则这个最优解一定是唯一的最优解。一、凸集凸集:设\(C\)为\(n\)维欧式......
  • 非线性规划——非线性规划的标准型(一)
    非线性规划是一种求解目标函数或约束条件中有一个或几个非线性函数的最优化问题的方法。运筹学的一个重要分支。20世纪50年代初,库哈(H.W.Kuhn)和托克(A.W.Tucker)提出了非线性规划的基本定理,为非线性规划奠定了理论基础。20世纪80年代以来,随着计算机技术的快速发展,非线性规划方......
  • 2.数值计算(1) --求解连续微分系统和混沌系统
    ✅作者简介:热爱科研的算法开发者,Python、Matlab项目可交流、沟通、学习。......
  • SCI1022 matlab求解
    MATLABforSCI1022assignment4:SolvingSudokuThisassignmentcontributes20%ofyourfinalmarkfortheMATLABmodule.Tocompletetheassignment,writeaMATLABlive-scriptornormalscript,andprovideapdf-filewithyourresultsincludingthefigures.......
  • 信号与系统基础复习:系统分析、求解方程、电路基础
    信息与系统总论信息是人类社会和自然界中需要传送、交换、存储和提取的抽象内容。信息存在于一切事物之中,事物的一切变化和运动都伴随着信息的交换和传送。各种各样的社会活动、无线电波的传播、计算机的运算等都是信息交换和传输的过程。信息是抽象的内容,为了传送和......
  • matlab中通过ode函数求解常微分方程附加简单的钟摆模型
    ✅作者简介:热爱科研的算法开发者,Python、Matlab项目可交流、沟通、学习。......
  • 问题求解
     为什么pppyyy2提示没有被使用【实际是调的到的】......
  • python spark 求解最大 最小 平均
    rdd=sc.parallelizeDoubles(testData);Nowwe’llcalculatethemeanofourdataset. 1LOGGER.info("Mean:"+rdd.mean());Therearesimilarmethodsforotherstatisticsoperationsuchasmax,standarddeviation,…etc.Everytimeoneofthismethodisin......
  • 最大信息系数——检测变量之间非线性相关性
    最后的效果就是这样的。很明显可以看到,左下角那个有点像三角函数的关系,Pearson系数(就是线性相关系数)为0,而MIC则有0.8。 摘自:http://tech.ifeng.com/a/20180323/44917506_0.shtml最大信息系数最大信息系数(MIC)于2011年提出,它是用于检测变量之间非线性相关性的最新方法。用于进行......
  • 迭代加深搜索与埃及分数求解
    迭代加深搜索,实质上是限定下界的深度优先搜索。即首先允许深度优先搜索K层,若没有发现可行解,再将K+1后重复以上步骤搜索,直到搜索到可行解。在迭代加深搜索的算法中,连续的深度优先搜索被引入,每一个深度约束逐次加1,直到搜索到目标为止。这样可以看出重复搜索了好多。但是它的好处在于:......