非线性优化（二）

• 一维搜索

 Fibonacci法(分数法)

设y=f(t)是区间[a,b]上的下单峰函数(如下图),在此区间内它有唯一极小点t.若在此区间内任取两点c和d,c<d,并计算函数值f(c)和f(d),可能出现以下两种情形:

  1. f(c)<f(d),这时极小点t必在区间[a,d]内.

  2. f(c)≥f(d),这时极小点t必在区间[c,b]内.

     在区间[a,b]内取两个不同的点,算出它们的函数值加以比较,就可以把搜索区间[a,b]缩小成[a,d]或[c,b](缩小后仍包含极小点).只要在[a,d]或[c,b]中再任取一点算出其函数值,并与f(c)或f(d)加以比较,就可以继续缩小搜索区间.只要缩小后的区间仍包含极小点t,则区间缩得越小,就越接近于函数的极小点,但计算函数值的次数也就越多.

 

     在区间[a,b]内取两个不同的点,算出它们的函数值加以比较,就可以把搜索区间[a,b]缩小成[a,d]或[c,b](缩小后仍包含极小点).只要在[a,d]或[c,b]中再任取一点算出其函数值,并与f(c)或f(d)加以比较,就可以继续缩小搜索区间.只要缩小后的区间仍包含极小点t,则区间缩得越小,就越接近于函数的极小点,但计算函数值的次数也就越多.

用Fibonacci法缩短区间的步骤如下

 

1.确定试点的个数n.根据缩短率δ,即可用式    (6.30)算出Fn,然后由递推公式确定最小的n.

2.选取前两个试点的位置.由递推公式可知,第一次缩短时的两个试点位置是

      t1 =  b0+(a0- b0) Fn-1/ Fn      

      t1， = a0 +(b0- a0) Fn-1/ Fn

它们在区间内是对称的.

 3.计算函数值f(t1)和f(t1，),并比较它们的大小.

    若f(t1)<f(t1，),

则取

         a1 =a0    b1=t1，  t2，= t1

 并令t2= b1+(a1- b1) Fn-2/ Fn-1

 否则,

取   a1 = t1    b1=b0  t2= t1，

 并令t2 ， = a1+(b1 - a1) Fn-2/ Fn-1.

 4.计算f(t2)或 f(t2，),如第三步那样一步步迭代.计算的一般公式为

                     tk= bk-1+(ak-1- bk-1) Fn-k/ Fn-k+1

                     tk ， = ak-1+(bk-1 - ak-1) Fn-k/ Fn-k+1

     其中k=1,2,…,n-1.            

5.当进行至k=n-1时,

                   tn-1=tn-1，=1/2  (an-2 +bn-2)

无法比较函数值f(tn-1)和f(tn-1 ，)的大小以确定最后区间,为此 ,取

                   tn-1=1/2  (an-2 +bn-2)

                   tn-1，= an-2 +(1/2+ε) (bn-2 - an-2)

 

 

（其中ε为任意小的数,在tn-1和tn-1，中以函数值较小者为近似极小点,相应的函数值为近似极小值.）并取得最终区间[an-2 ,tn-1，]或[tn-1, bn-2].

     由上述分析可知,Fibonacci使用对称搜索的方法,逐步缩短所考察的区间,它能以尽量少的函数求值次数,达到预定的某一缩短率.

 

 

0.618法(黄金分割法)

当用Fibonacci法以n个试点来缩短某一区间时,区间长度的第一次缩短率为Fn-1/ Fn,其后各次分别为

      Fn-2/ Fn-1,Fn-3/ Fn-2 ,…,F1/F2

 把以上数列分为奇数项F2k-1/ F2k和偶数项F2k/ F2k+1,这两个数列收敛于同一个极限0.618.

   当用0.618方法时,计算n个试点的函数值可以把原区间[a0,b0]连续缩短n-1次,因为每次缩短率均为μ,故最后的区间长度为 (b0- a0)μn-1故当已知缩短的相对精度为δ时,可以用下式计算试点个数n:  μn-1≤δ

     当然,也可以不预先计算试点的数目n,而在计算过程中逐次加以判断,看是否已满足了提出的精度要求.

     0.618法是一种等速对称进行试探的方法,每次的试点均取在区间长度的0.618倍和0.328倍处.

 

• 无约束极值问题的解法

无约束极值问题可表述为    

            (MP)   minf(X) , X∈ En               (3.1)

   解上述问题常用迭代法,迭代法大体分为两类:

一类要用到函数值的一阶导数(或)二阶导数,由于用到了函数的解析性质,故称为解析法;

另一类在迭代过程中仅用到函数值,而不要求函数的解析性质,这类方法称为直接法.本节介绍几种常用的基本方法,其中前三种属解析法,后面一种属直接法.

 

1.梯度法(最速下降法)

假定无约束极值问题 minf(X)，X∈ En 中f(X)有一阶连续偏导数,具有极小点X*.以X(k)表示极小点的第k次近似,为了求其第k+1次近似点X(k+1),在X(k)点沿方向P(k)作射线

            X=X(k)+λP(k)     (λ≥0)

将f(X)在X(k)点处展成泰勒级数f(X)=f(X(k)+λP(k))=f(X(k))+λ▽f(X(k))T P(k)+o(λ)

 对于充分小的λ,只要

▽f(X(k))T P(k)<0            (3.2)

 即可保证 f(X(k)+λP(k))< f(X(k)).

 这时若取 X(k+1)=X(k)+λP(k),就能使目标值下降.    

 

下降方向的选取

  设‖P(k)‖≠0 , ▽f(X(k))≠0.由

  ▽f(X(k))TP(k)=‖▽f(X(k))‖·‖P(k)‖cosθ

式中θ为向量▽f(X(k))与P(k)的夹角,当▽f(X(k))与

P(k)反向时,θ=180。，cosθ=-1.这时式(3.2)成立,而且其左端取最小值,称方向 P(k)=- ▽f(X(k)) 为负梯度方向.它是使函数值(在X(k)的附近)下降最快的方向.

 

步长的选取

     若选取负梯度方向,则存在λ满足不等式

       f(X(k)-λ▽f(X(k)))< f(X(k))

   通过负梯度方向的一维搜索,来确定使f(X)最小的λk,这种梯度法就是最速下降法.近似最佳步长:                                                         

                                                                    (3.3)

 

 

最速下降算法步骤:

     (1) 给定初始近似点X(0)及精度ε>0 及k=1.

     (2) 求搜索方向 P(k)=-▽f(X(k)).

     (3)若‖▽f(X(k))‖<ε,则X(k)为近似极小点.

否则求λk ，使

          f(X(k)+λkP(k))=minλ≥0 f(X(k)+λP(k))

1. 令X(k+1):=X(k)+λkP(k), k:=k+1,转(2).

 

有时将搜索方向P(k)的模规范化为1，即

      P(k)=-▽f(X(k))/‖▽f(X(k))‖

 

 

 

 

 

 

 

 

设凸函数f(X)有一阶连续偏导数,若由最速下降法所得的点列X(k)有界,则有

(1). f(X(k))单调下降;

(2).X(k)的极限点X* 满足▽f(X*)=0;

(3). X*为全局极小点.

 

2.共轭梯度法

 共轭方向

      定义1  设An*n为对称正定矩阵,若对向量0≠P(1),P(2)…,P(n)∈En,满足

          (P(i))TAP(j)=0,i≠j,i=1,2,…,n

  则称该向量组为A共轭.或称为A的n个共轭方向.若A=In(单位阵),则A共轭就为正交.

说明：

设A为阶实对称正定矩阵，如果有两个n维向量S1和满足S2

S1AS2=0 （1）

则称向量S1与S2对于矩阵A共轭。如果A为单位矩阵，则式(1)即成为S1S2，这样两个向量的点积（或称内积）为零，此二向量在几何上是正交的，它是共轭的一种特例。

设A为对称正定矩阵，若一组非零向量S1，S2，…Sn满足

SiASj=0 （i≠j） （2）

则称向量系Si（i=1,2，…n）为关于矩阵A共轭。

共轭向量的方向称为共轭方向

 

定理8 设A为n×n对称正定阵,P(1),P(2),…,P(n)为A共轭的非零向量,则这一组向量线性无关.

 

正定二次函数极小化问题

       minf(X)=(1/2)XTAX+BTX+c              (3.4)

式中A为n×n对称正定阵;X,B∈En;c为常数.

 定理9 设向量P(i) ,i=1,2,…,n-1,为A共轭,则从任一点X(0)出发,相继以P(0),P(1),…,P(n-1)为搜索方向的下述算法:

             Minλ f(X(k)+λP(k))=f (X(k)+λkP(k))

               X(k+1) =X(k)+λkP(k)

经n次一维搜索收敛于问题(3.4) 的极小点X *.

 

 对于问题(3.4)来说,由于A为对称正定阵,故存在唯一极小点X* ,它满足方程组

           ▽f(X)=AX+B=0

且具有形式   X* =-A-1B

    如果已知某共轭向量组P(0),P(1),…,P(n-1),由定理9可知,问题(3.4)的极小点X*可通过下列算法得到:

         X(k+1) =X(k)+λkP(k) , k=0,1,2,…,n-1     

     λk:Minλ f(X(k)+λP(k))                   (3.8)

         X(n) =X*

     这种从任一点X(0)∈En出发,依次沿某组共轭方向进行最优一维搜索,求解无约束极值的方法,称为共轭方向法.

构造正定二次函数的共轭梯度法.

      任取初始近似点X(0),并取初始搜索方向为此点的负梯度方向,即  P(0) = -▽f(X(0))沿射线X(0) +λP(0)进行一维搜索,得

         X(1) =X(0)+λ0P(0)

     λ0:Minλ f(X(0)+λP(0))

   一般地,从点X(k)出发,沿方向P(k)进行最优一维搜索,有  ▽f(X(k+1))TP(k) =0 （结合3.5）.从而  

     λk=-▽f(X(k))TP(k)/(P(k))TAP(k)   (3.9)

   按 P(k+1)=-▽f(X(k+1))+βk P(k)   (3.10)

来产生搜索方向.欲选择βk使P(k+1)和P(k)为A共轭

则应有  0=(P(k+1))TAP(k)

         =-▽f(X(k+1))TAP(k)+βk( P(k))TAP(k)

从而对 K=0,1,…,n-2,取

  βk=▽f(X(k+1))TAP(k)/[(P(k))TAP(k)] (3.11)

   下面简化βk的表示式:记gk= ▽f(X(k)),由(3.5)知,

 AP(k) =(gk+1-gk)/λk,代入(3.11)得

 βk=gk+1T(gk+1-gk)/[(P(k))T(gk+1-gk)]  (3.12)

∵0=gk+2T P(k+1)=-gk+2Tgk+1+βkgk+2TP(k)=-gk+2Tgk+1

即 gk+2Tgk+1=0 , K=0,1,…,n-2.      (3.13)

共轭梯度法的计算步骤:

   (1) 选取初始近似X(0),给出允许误差ε>0.

   (2)计算 gk=▽f(X(k)),若‖gk+1‖=0,则stop,

得 X*=X(k).否则转下一步.

 (3)按(3.16)—（3.19)计算X(k+1) 和P(k+1).

 (4)若k=n,则stop,得 X*=X(k).否则goto(2).

附注:

  对于二次函数的情形,从理论上说,进行n次迭代即可达到极小点.但是,在实际计算中,由

于数据的舍入以及计算误差的积累,往往做不到这点此外,由于n维问题的共轭方向最多只

有n个,在n步以后继续如上进行是没有意义的.因此,在实际应用时,如迭代到n步还不收敛,

就将X(n)作为新的初始近似,重新开始迭代.

 

 

3.变尺度法

假定无约束极值问题的目标函数f(X)∈C(2) ,X(k)为其极小点X*的某一近似.由点X(k)附近取

f(X)的二阶泰勒多项式逼近

       f(X)≈f(X(k))+▽f(X(k))T△X

               +(1/2)△XTH(X(k))△X                                   

则其梯度为 ▽f(X)≈▽f(X(k))+H(X(k))△X    

若H(X(k))正定,则这个近似函数的极小点满足

          ▽f(X(k))+H(X(k))△X=0

从而   X= X(k)-H(X(k))-1▽f(X(k))

 

算法 X(k+1)= X(k)-H(X(k))-1▽f(X(k))称为Newton法.Newton法的计算步骤：

   1 给定初始点X(0),梯度允许误差ε>0,k:=0   

   2 若‖▽f(X(0))‖≤ε,则X(0)即为近似极小

点,停止迭代.否则,转下一步.

   3 取P(k)=-[▽2f(X(k))]-1 ▽f(X(k)),

    令X(k+1):=X(k)+λkP(k), k:=k+1,转(2)

当X(0)靠近X*时Newton法二次收敛， X(0)远离X*时Newton法可能发散。

算法  P(k)=-H(X(k))-1gk         (3.25)

 X(k+1)= X(k)+ λk P (k)          (3.26)

  λk:minf(X(k)+ λ P (k))    (3.27)

称为广义牛顿法.式(3.25)确定的搜索方向称为f(X)在X(k)的Newton方向.

     在实际问题中为了避免计算二阶导数H(X(k))和H(X(k))-1,以下构造另一个矩阵H(k),用它来直接逼近二阶导数矩阵的逆阵H(X(k))-1.

当f(X)是二次函数时,其Hesse矩阵为常数阵,可知 ▽f(X(k+1))-▽f(X(k))=A(X(k+1)-X(k)),或

 X(k+1)-X(k)=A-1[▽f(X(k+1))-▽f(X(k))]  (3.28)

对于非二次函数,要求其Hesse矩阵的逆阵的第k+1次近似矩阵H(k+1)满足关系式:

X(k+1)-X(k)=H(k+1)[▽f(X(k+1))-▽f(X(k))] (3.29)

(3.29)称为拟牛顿条件.

    若令   △G(k)=▽f(X(k+1))-▽f(X(k))

         △X(k)= X(k+1)-X(k)

则式(3.29)变为 △X(k)= H(k+1)△G(k)

设H(k)已知,并用下式求H(k+1)(假定H(k)和H(k+1)都为对称矩阵):

        H(k+1)= H(k)+△H(k)          (3.30)

上式中△H(k)为第k次校正矩阵.H(k+1)应满足拟牛顿条件,即要求

    △X(k) =(H(k)+△H(k))△G(k)或  △H(k)△G(k)=△X(k)-H(k)△G(k)   (3.31)

设△H(k)=△X(k)(Q(k))T- H(k)△G(k)(W(k))T (3.32)

其中Q(k)和W(k)为两个待定向量.

    将式(3.32)代入(3.31)得

 

  于是得到校正矩阵

   △H(k)=[△X(k)(△X(k))T ]/[(△G(k))T△X(k)]

-[H(k)△G(k)(△G(k))TH(k)]/[(△G(k))TH(k)△G(k)]

   从而得到

  H(k+1)=H(k)+[△X(k)(△X(k))T]/[(△G(k))T△X(k)]

-[H(k)△G(k)(△G(k))TH(k)]/[(△G(k))TH(k)△G(k)]

                         -------(3.36)

矩阵H(k+1)称为尺度矩阵.

变尺度法的计算步骤

   1 给定初始点X(0),梯度允许误差ε>0.   

   2 若‖▽f(X(0))‖≤ε,则X(0)即为近似极小

点,停止迭代.否则,转下一步.

   3 令H(0)=I(单位矩阵),P(0)=-H(0)▽f(X(0)),

在P(0)方向进行一维搜索,确定最佳步长λ0 :

           minλ f(X(0)+λP(0))=f(X(0)+λ0P(0)).

    4  一般地,设已得到近似点X(k),算出▽f(X(k)),

若  ‖▽f(X(k))‖≤ε,则X(k)即为近似值,停止

迭代；否则,按式(3.36)计算H(k),并令

         P(k)= -H(k)▽f(X(k))

在P(k)方向进行一维搜索,确定最佳步长λk:

             minλf(X(k)+λP(k)) =f(X(k)+λkP(k))

于是下一个近似点为

             X(k+1)= X(k)+λkP(k)

   5 若X(k+1)满足精度要求,则X(k+1)为所求近

似解.否则,转回第4步,直至某点满足精度为止.

       与共轭梯度法类似,若迭代n次未收敛,则取X(0):=X(n),以新的X(0)为初点,开始新一轮迭代.上述方法亦称为DFP法或 DFP变尺度法.