无约束最优化方法

我会带着你远行 2022-08-05 01:21 163阅读 0赞

# 无约束最优化方法 #

最优化(Optimization)是工程技术、经济管理、  
科学研究、社会生活中经常遇到的问题, 如：结构设计，资源分配，生产计划，运输方案等等。  
最优化: 在一定条件下，寻求使目标最大(小)的决策  
解决优化问题的手段  
• 经验积累，主观判断  
• 作试验，比优劣  
• 建立数学模型，求解最优策略

## 优化问题的数学模型 ##

无约束优化问题标准形式：  
求n维决策变量x=\[x1,x2,…,xn\]T在可行域Ω中  
使得目标函数f(x)→min,x⊆Ω⊆Rn  
⇒

![这里写图片描述][20150529104029026]  
**可行解**只满足(2)， **最优解**满足(1)(2)  
**无约束优化**只有(1)， **约束优化**需要(1)(2)

## 无约束优化的基本方法 ##

### 求解无约束优化的基本思路：[下山法][Link 1] ###

#### 基本思想： ####

和爬山法基本思想一样，在Rn中某一点，确定一个**搜索放向**然后沿着该方向的**移动步长**，得到使目标函数下降的最新点。

#### 迭代步骤： ####

Step 1 初始化：初始点x0，终止准则等  
Step 2 迭代改进：方向dk，步长αk

xk\+1=xk\+(αd)k 且 f(xk\+1)<f(xk)

Step 3 终止检验：得到近优解或  k\+1⇒k 转2

> 选择dk,αk使得f下降更快⇒利用不同的算法

### 最速下降法（梯度法） ###

泰勒级数：

![2012091914365021.png][]  
其中  Δx 可正可负，但必须充分接近于0。  
 X 沿  D 方向移动步长  α 后，变为  X\+αD 。由泰勒展开式： ![2012072008442633.png][]  
目标函数： ![2012072008450656.png][]  
 α 确定的情况下即最小化： ![2012072008452618.png][]  
向量的内积何时最小？当然是两向量方向相反时。所以X移动的方向应该和梯度的方向相反。

> **下降方向：**∇fT(xk)dk<0  
> **最速下降方向：**dk=−∇f(xk)负梯度方向  
> **迭代改进形式：**xk\+1=xk−∇f(xk)

算法特点：在初始状态改进较快，在最优解附近改进缓慢。

求函数f(x1,x2)=x21\+4x22在极小点，以初始点X0=(1,1)T。

#include"matrix.h"//大家自己YY一个矩阵类吧
    #include<iostream>
    #include<iomanip>
    #include<cmath>
    #include<limits>
    #include<cassert>
    
    using namespace std;
    
    const int SIZE=2;
    const double ZETA=0.001;
    
    inline double func(Matrix<double> &X){
        assert(SIZE==X.getRows());
        double x1=X.get(0,0);
        double x2=X.get(1,0);
        return x1*x1+4*x2*x2;
    }
    
    inline Matrix<double> gradient(Matrix<double> &X){
        assert(SIZE==X.getRows());
        double x1=X.get(0,0);
        double x2=X.get(1,0);
        Matrix<double> rect(SIZE,1);
        rect.put(0,0,2*x1);
        rect.put(1,0,8*x2);
        return rect;
    }
    
    inline Matrix<double> Hesse(Matrix<double> &X){
        Matrix<double> rect(SIZE,SIZE);
        rect.put(0,0,2);
        rect.put(0,1,0);
        rect.put(1,0,0);
        rect.put(1,1,8);
        return rect;
    
    }
    
    int main(int agrc,char *argv[]){
        Matrix<double> X(SIZE,1);
        X.put(0,0,1);
        X.put(1,0,1);
        int iteration=0;
        double value=func(X);
        double newValue=numeric_limits<double>::max();
        while(++iteration){
            Matrix<double> G=gradient(X);
            double factor=((G.getTranspose()*G).get(0,0))/((G.getTranspose()*Hesse(X)*G).get(0,0));
            for(int i=0;i<G.getRows();++i)
                    for(int j=0;j<G.getColumns();++j)
                        G.put(i,j,G.get(i,j)*factor);
            Matrix<double> newX=X-G;
            //cout<<"X=["<<newX.get(0,0)<<","<<newX.get(1,0)<<"]"<<endl;
            newValue=func(newX);
            if(fabs(newValue-value)<ZETA)
                    break;
            else{
                    X=newX;
                    value=newValue;
            }
            //cout<<"本次迭代找到的极小值"<<value<<endl;
        }
        cout<<"迭代次数"<<iteration<<endl;
        cout<<"极小值"<<value<<endl;
        return 0;
    }

> 这份算法使用的是利用Hesse矩阵和泰勒公式计算最优步长  
> 若f(X)具有二阶连续偏导，由泰勒展开式可得：  
> 
> 
> ![2012072010525191.png][]  
> H是f(X)的Hesse矩阵。  
> ![2012072010532298.png][]  
> 可得最优步长： ![2012072010593228.png][]  
> g是f(X)的梯度矩阵。  
> 此时： ![2012072010595473.png][]  
> 可见最速下降法中最优步长不仅与梯度有关，而且与Hesse矩阵有关。

### 牛顿迭代法 ###

将f(xk\+1)在xk点作泰勒展开至二阶项，用d替代dk，有

f(xk\+1)=f(xk\+d)=f(xk)\+∇fT(xk)d\+12dT∇2f(xk)d

求  d 使得  f(xk\+1) 极小  ⇒ 右端对  d 的导数为0  ⇒∇f(xk)\+∇2f(xk)d=0  
可得到

> **牛顿方程** ∇2f(xk)d = −∇f(xk)  
> **牛顿方向** dk = −(∇2f(xk))−1∇f(xk)  
> **迭代格式** xk\+1 = xk−(∇2f(xk))−1∇f(xk)  
> 其中∇2f(xk)为[Hessen矩阵][Hessen]

可见x的搜索方向是−(∇2f(xk))−1∇f(xk)，函数值在此方向上下降，就需要它与梯度方向相反，即−∇f(xk)(∇2f(xk))−1∇f(xk)<0，所以要求在每一个迭代点上Hessian矩阵必须是正定的。  
算法特点：  
牛顿法是二次收敛的，并且收敛阶数是2。一般目标函数在最优点附近呈现为二次函数，于是可以想像最优点附近用牛顿迭代法收敛是比较快的。而在开始搜索的几步，我们用梯度下降法收敛是比较快的。将两个方法融合起来可以达到满意的效果。  
收敛快是牛顿迭代法最大的优点，但也有致命的缺点：Hessian矩阵及其逆的求解计算量大，更何况在某个迭代点Xk处Hessian矩阵的逆可能根本就不存在（即Hessian  
矩阵奇异），这样无法求得Xk\+1。

求f(x1,x2)=(x1−2)4\+(x1−2x2)2的极小点，初始点取X=(0,3)T

#include"matrix.h" 
    #include<iostream> 
    #include<iomanip> 
    #include<cmath> 
    #include<limits> 
    #include<cassert> 
    
    using namespace std; 
    
    const int SIZE=2; 
    const double ZETA=0.001; 
    
    inline double func(Matrix<double> &X){ 
        assert(SIZE==X.getRows()); 
        double x1=X.get(0,0); 
        double x2=X.get(1,0); 
        return pow(x1-2,4)+pow(x1-2*x2,2); 
    } 
    
    inline Matrix<double> gradient(Matrix<double> &X){ 
        assert(SIZE==X.getRows()); 
        double x1=X.get(0,0); 
        double x2=X.get(1,0); 
        Matrix<double> rect(SIZE,1); 
        rect.put(0,0,4*pow(x1-2,3)+2*(x1-2*x2)); 
        rect.put(1,0,-4*(x1-2*x2)); 
        return rect; 
    } 
    
    inline Matrix<double> Hesse(Matrix<double> &X){ 
        Matrix<double> rect(SIZE,SIZE); 
        double x1=X.get(0,0); 
        double x2=X.get(1,0); 
        rect.put(0,0,12*pow(x1-2,2)+2); 
        rect.put(0,1,-4); 
        rect.put(1,0,-4); 
        rect.put(1,1,8); 
        return rect; 
    } 
    
    int main(int agrc,char *argv[]){ 
        Matrix<double> X(SIZE,1); 
        X.put(0,0,0); 
        X.put(1,0,3); 
        int iteration=0; 
        double value=func(X); 
        double newValue=numeric_limits<double>::max(); 
        while(++iteration){ 
            Matrix<double> G=gradient(X); 
            Matrix<double> H=Hesse(X); 
            Matrix<double> newX=X-H.getInverse()*G; 
            cout<<"X=["<<newX.get(0,0)<<","<<newX.get(1,0)<<"]"<<endl; 
            newValue=func(newX); 
            if(fabs(newValue-value)<ZETA) 
                break; 
            else{ 
                X=newX; 
                value=newValue; 
            } 
            cout<<"本次迭代找到的极小值"<<value<<endl; 
        } 
        cout<<"迭代次数"<<iteration<<endl; 
        cout<<"极小值"<<value<<endl; 
        return 0; 
    }

### 拟牛顿法 ###

Hessian矩阵在拟牛顿法中是不计算的，拟牛顿法是构造与Hessian矩阵相似的正定矩阵，这个构造方法，使用了目标函数的梯度（一阶导数）信息和两个点的“位移”（Xk−Xk−1）来实现。  
有人会说，是不是用Hessian矩阵的近似矩阵来代替Hessian矩阵，会导致求解效果变差呢？事实上，效果反而通常会变好。  
拟牛顿法与牛顿法的迭代过程一样，仅仅是各个Hessian矩阵的求解方法不一样。  
在远离极小值点处，Hessian矩阵一般不能保证正定，使得目标函数值不降反升。而拟牛顿法可以使目标函数值沿下降方向走下去，并且到了最后，在极小值点附近，可使构造出来的矩阵与Hessian矩阵“很像”了，这样，拟牛顿法也会具有牛顿法的二阶收敛性。

对目标函数f(X)做二阶泰勒展开：

f(x)=f(xk\+1)\+(x−xk\+1)T∇f(xk\+1)\+12(x−xk\+1)TG(xk\+1)(x−xk\+1)

其中  ∇2f(xk) ⇔G(xk\+1)  
上式两边对  x 求导：

∇f(x)=∇f(xk\+1)\+G(xk\+1)(x−xk\+1)

当  x=xk 的时候，有

G−1(xk\+1)\[∇f(xk\+1)−∇f(xk)\]=xk\+1−xk

这里我们用  Hk 表示在点  xk 处的Hessian矩阵的 **逆矩阵**

Hk\+1\[∇fd(xk\+1)−∇f(xk)\]=xk\+1−xk ......(5)

(5)式就是 **拟牛顿方程**  
下面提供两种拟牛顿法

![这里写图片描述][20150529120043981]

### 非线性最小二乘(Least Square)拟合 ###

问题：  
给定(ti,yi),i=1,…n, 拟合一个函数y=f(t,x)，  
其中x为待定的参数向量, f 对x非线性。  
记误差ri(x)=yi−f(ti,x)     r(x)=(r1(x),…,rn(x))T  
优化模型：

minxR(x)=rT(x)r(x)

记r(x)的雅克比行列式为J(x)=(∂ri∂xj)n×m  
∇R=2J(x)Tr(x)∇2R=2J(x)TJ(x)\+2SS=∑ni=1ri(x)∇2ri(x)∇2ri(x)=(∂2ri∂xk∂xl)n×m

牛顿法要计算Hessian矩阵，其中S计算量大  
若f 对x线性, 则化为线性最小二乘拟合, 此时  
S=0  
特定算法考虑如何忽略或近似矩阵S

#### Gauss-Newton算法(GN)：忽略矩阵S  ∇R=2J(x)Tr(x)∇2R=2J(x)TJ(x)   牛顿方程：∇2f(xk)d=−∇f(xk)   f用R代替，下降方向dk满足     J(x)TJ(x)dk=−J(x)Tr(x)     GN算法收敛性特点:    S=∑ni=1ri(x) 收敛性依赖  f对x的线性程度，即偏差r的大小（”小残差问题“） ####

#### LM算法: Levenbery (1944), Marquardt (1963) ####

G-N算法修正

J(x)TJ(x)dk=−J(x)Tr(x)

⇓

(J(x)TJ(x)\+αkI)dk=−J(x)Tr(x)

其中  αk>0 为修正参数，防止  JTJ 出现病态。  
LM算法特点  
 dk 位于GN方向（  αk 很小）和负梯度方向(  αk 很大)之间

LM算法的另一种观点  
min∥r(xk)\+J(xk)(x−xk)∥2  
s.t.     ∥x−xk∥2≤hk    (∗)置信域  
其解为

(J(x)TJ(x)\+αkI)dk=−J(x)Tr(x)

若GN方向满足约束(\*)，  αk=0 ,否则  ak>0

置信域算法特点：  
先确定最长步长hk，再确定方向dk

[20150529104029026]: /images/20220805/2167890a38fa40e1be90ab71c707c2e1.png
[Link 1]: http://blog.csdn.net/u013007900/article/details/43525733
[2012091914365021.png]: /images/20220805/9d0a1489393a4bc394671cd692f910c9.png
[2012072008442633.png]: /images/20220805/43473f756e33485697be555b34e306d9.png
[2012072008450656.png]: /images/20220805/e1ef4ce9d4a449a6b2ea425fdb4c1c38.png
[2012072008452618.png]: /images/20220805/d1cf31918b344d5f98217671d718e1f6.png
[2012072010525191.png]: /images/20220805/102cfb7878514474bc7ceb3dedaa5f2a.png
[2012072010532298.png]: /images/20220805/d5721dcb2e11451eb194c781bb43cb02.png
[2012072010593228.png]: /images/20220805/05c2b14c58b6407d89368c40d68b6134.png
[2012072010595473.png]: /images/20220805/9dad070f59fd4bd59b3b182837ce0dad.png
[Hessen]: http://baike.baidu.com/link?url=E2wdH1BhTKFos1JYd_NhHvkmGAyABsHJDIgSyygvzTuoAs41tvb4OGzzWThbUMwobdBAEPpS7nzia-GQUHHf7a
[20150529120043981]: /images/20220805/0cd4338dc2524f34943b6da1e27d65e4.png