最优化算法

最新推荐文章于 2026-07-09 13:49:15 发布

原创最新推荐文章于 2026-07-09 13:49:15 发布 · 2.1k 阅读

9 ·

本内容遵循CC 4.0 BY-SA版权协议

最优化算法

1. 一元函数最优化算法

1.1 三分法

普通的三分法是比较常用的求一元函数最值的方法

以求函数 $f (x)$ 的最小值为例，假设函数 $f (x)$ 在区间 $[a, b]$ 上先严格单调递减，后严格单调递增，那么可以通过迭代缩小最值点可能在的区间 $[l, r]$

令 $l = a, r = b$
$m_l=\frac{2l+r}{3},m_r=\frac{l+2r}{3}$ ，即取 $[l, r]$ 的三等分点
如果 $f(m_l)>f(m_r)$ ，那么 $l=m_l$ ，否则 $r=m_r$
如果达到最大迭代次数，那么停止，否则转2

每经过一次迭代，区间 $[l, r]$ 的长度变为原来的 $\frac{2}{3}$ ，因此其精度是一阶收敛的

/*三分法

    f:目标函数
    l:区间下界
    r:区间上界
    max_step:最大迭代次数

*///2019.2.23
template <class fun>
double oneThirdIterate(fun f,double l,double r,int max_step){
    double ml,mr;
    rep(it,1,max_step){
        ml=(2*l+r)/3,mr=(l+2*r)/3;
        if(f(ml)>f(mr))l=ml;
        else r=mr;
    }
    return l;
}

1.2 黄金分割三分法

黄金分割三分法是在三分法的基础上改进得到的算法，在三分法中，我们可以把 $m_l,m_r$ 的位置选得更接近 $\frac{l+r}{2}$ ，这样收敛速度会更快，在此基础上，还可以减少目标函数的调用来提高效率

取 $m_l=gl+(1-g)r,m_r=(1-g)l+gr$ ，其中 $g=\frac{\sqrt 5-1}{2}$

在这里插入图片描述

如图所示，这样选取可以沿用前一次迭代的 $f(m_l)$ 值作为新的 $f(m_r)$ 值（或沿用前一次迭代的 $f(m_r)$ 值作为新的 $f(m_l)$ 值），避免了目标函数的重复调用

这样可以提高算法效率，但会带来较大的舍入误差

/*黄金分割迭代

    f:目标函数
    l:区间下界
    r:区间上界
    max_step:最大迭代次数

*///2019.2.20
template <class fun>
double GoldenIterate(fun f,double l,double r,int max_step){
    double g=(sqrt(5.0)-1)/2,ml=g*l+(1-g)*r,mr=(1-g)*l+g*r,lval=f(ml),rval=f(mr);
    rep(it,1,max_step){
        if(lval<rval)r=mr,mr=ml,ml=g*l+(1-g)*r,rval=lval,lval=f(ml);
        else l=ml,ml=mr,mr=(1-g)*l+g*r,lval=rval,rval=f(mr);
    }
    return l;
}

2. 多元函数最优化算法

2.1 模拟退火（SAA）

模拟退火算法是来自冶金的算法，当温度过高时，粒子会剧烈运动，随着时间延长，温度逐渐降低，粒子运动速度下降。在算法中，随机确定下一步的位置，如果下一步更优，那就移动到该位置，否则根据当前温度计算粒子运动的概率，以此概率随机确定是否移动到该位置

随机确定起点 $x=x_0$ ，确定初始温度 $T=T_0$
温度 $T=T\eta(0<\eta<1)$
随机选取 $x$ 附近的点 $x^{'}$ （ $d (x, x^{'})$ 需要收敛到0）
若 $f (x^{'}) < f (x)$ ，那么 $x = x^{'}$ ，否则若 $P (f (x), f (x^{'}), T) > r a n d o m (0, 1)$ ，那么 $x = x^{'}$
如果达到最大迭代次数，那么停止，否则转2

/*SAA模拟退火算法

    f:目标函数
    x:初始点
    d:初始步长
    mi:变量下界
    ma:变量上界
    max_step:最大迭代次数
    T:初始温度

*///2019.2.15
double uniformRand(double l,double r){return rand()/32767.0*(r-l)+l;}
ll uniformRandInt(int l,int r){return rand()%(r-l+1)+l;}
template <class fun>
tensor SAA(fun f,tensor x,tensor d,tensor mi,tensor ma,int max_step,double T){
    int n=x.size();
    double yita=0.97,fval=f(x),next_fval;
    tensor next;
    rep(it,1,max_step){
        next=x;
        rep(i,0,n-1)next[i]=(rand()&1)?min(next[i]+d[i],ma[i]):max(next[i]-d[i],mi[i]);
        next_fval=f(next);
        if(next_fval<fval || uniformRand(0,1)<exp((fval-next_fval)/T))x=next,fval=next_fval;
        rep(i,0,n-1)d[i]*=yita;
        T*=yita;
    }
    return x;
}

2.2 粒子群算法（PSO）

粒子群算法是由鸟群的运动启发得到的算法，想象一群鸟在寻找温暖的地方，每只鸟有自己的位置 $P_i$ ，有自己的速度 $v_i$ ，同时记录了自己飞过的最温暖的位置 $pbest_i$ ，鸟群间有信息交流，共享所有鸟飞过的位置中最温暖的位置 $g b e s t$ 。每只鸟在飞行时，有沿自身速度方向飞的趋势，有向 $pbest_i$ 飞的趋势，有向 $gbest_i$ 飞的趋势，因此，其位置由以上三个量确定

随机确定 $m$ 个粒子的位置 $P_i$
更新每个粒子的 $pbest_i$ ，更新全局最优解 $g b e s t$
对于每个粒子，其下一步位置 $P_i'=P_i+wv_i+c_1r_1(pbest_i-P_i)+c_2r_2(gbest-P_i)$ ，其中 $w$ 是惯性系数， $c_1$ 是认知加速常数， $c_2$ 是社会加速常数， $r_1$ 和 $r_2$ 是 $[0, 1]$ 范围内的随机数
如果达到最大迭代次数，那么停止，否则转2

/*PSO粒子群算法

    f:目标函数
    vmax:最大飞行速度（参考值：变量变化范围的10%）
    mi:变量下界
    ma:变量上界
    max_step:最大迭代次数
    m:粒子群规模（参考值:50）
    w:惯性系数（参考值:0.6）
    c1:"认知"加速常数（参考值：1.7）
    c2:"社会"加速常数（参考值：1.7）

*///2019.2.17
template <class fun>
tensor PSO(fun f,tensor vmax,tensor mi,tensor ma,int max_step,int m,double w,double c1,double c2){
    int n=vmax.size();
    tensor p[m],g,x[m],v[m];
    double pbest[m],gbest,now;
    rep(i,0,m-1){
        x[i].resize(n),v[i].resize(n),p[i].resize(n);
        rep(j,0,n-1){
            x[i][j]=uniformRand(mi[j],ma[j]);
            v[i][j]=uniformRand(-vmax[j],vmax[j]);
        }
        p[i]=x[i],pbest[i]=f(x[i]);
        if(i==0 || pbest[i]<gbest)gbest=pbest[i],g=x[i];
    }
    rep(it,1,max_step){
        rep(i,0,m-1){
            v[i]=w*v[i]+c1*uniformRand(0,1)*(p[i]-x[i])+c2*(g-x[i]);
            rep(j,0,n-1)v[i][j]=max(min(v[i][j],vmax[j]),-vmax[j]);
            x[i]+=v[i];
            now=f(x[i]);
            if(now<pbest[i])pbest[i]=now,p[i]=x[i];
            if(now<gbest)gbest=now,g=x[i];
        }
    }
    return g;
}

2.3 梯度下降法（GD）

梯度下降法是机器学习中常用的算法，但是需要依赖函数的导数。函数的梯度反方向是函数值下降最快的方向，因此可以以公式
$x_{n+1}=x_n-\gamma\ \text{gradient}\left(f(x)\right)$
迭代若干次，逼近最优解

（函数的梯度即为函数对每个自变量的偏导数组成的向量）

/*GD梯度下降法

    f_:目标函数的梯度函数
    x:起始点
    gamma:学习率（参考值：0.4）
    max_step:最大迭代次数

*///2019.2.21
template <class fun>
tensor GD(fun f_,tensor x,double gamma,int max_step){
    rep(it,1,max_step)x-=gamma*f_(x);
    return x;
}

2.4 坐标下降法（CD）

坐标下降法是依赖一元函数最优化算法的多元函数最优化算法，设目标函数的维数为 $n$ ，在每一次迭代中，固定自变量中的 $n - 1$ 维，调整剩下的一维达到局部最优，重复上述过程即可

注意，某些情况下，坐标下降法的收敛速度可能较慢

/*CD坐标下降法

    f:目标函数
    x:初始点
    mi:变量下界
    ma:变量上界
    max_step:最大迭代次数
    max_step2:一维搜索的最大迭代次数

*///2019.2.23
template <class fun>
tensor CD(fun f,tensor x,tensor mi,tensor ma,int max_step,int max_step2){
    ll d=mi.size(),xl,xr;
    static fun _f=f;
    class Fun{
    public:
        int idx;
        tensor x;
        double operator () (double y){
            x[idx]=y;
            return _f(x);
        }
    }F;
    F.x=x;
    rep(it,1,max_step){
        rep(i,0,d-1){
            F.idx=i;
            double l=GoldenIterate(F,mi[i],x[i],max_step2);
            double r=GoldenIterate(F,x[i],ma[i],max_step2);
            x[i]=F(l)<F(r)?l:r;
            F.x[i]=x[i];
        }
    }
    return x;
}

2.5 Nelder–Mead方法（NM）

Nelder–Mead方法又名下山单纯形法，维护高维空间中的多胞体包裹住最优解，通过迭代减小多胞体的测度，确定最优解

随机确定 $x_1,x_2,\dots,x_{n+1}$
根据 $f(x_i)$ 对其从小到大进行排序
计算 $x_1,x_2,\dots,x_n$ 的中心 $x_o$
Reflection 计算Reflection点 $x_r=x_o+\alpha(x_o-x_{n+1}),\alpha>0$ ，如果 $f(x_1)\le f(x_r)<f(x_n)$ ，那么将 $x_{n+1}$ 替换为 $x_r$ 并转2
Expansion 如果 $f(x_r)<f(x_1)$ ，那么计算 $x_e=x_o+\gamma(x_r-x_o),\gamma>1$ ，如果 $f(x_e)<f(x_r)$ ，那么将 $x_{n+1}$ 替换为 $x_e$ 并转2 ，否则将 $x_{n+1}$ 替换为 $x_r$ 并转2
Contraction 计算contracted点 $x_c=x_o+\rho(x_{n+1}-x_o),0<\rho \le0.5$ ，如果 $f(x_c)<f(x_{n+1})$ ，那么将 $x_{n+1}$ 替换为 $x_c$ 并转2
Shrink 将 $x_i$ 替换为 $x_1+\gamma(x_i-x_1),i=2,3,\dots,n+1$
如果达到最大迭代次数，那么停止，否则转2

/*Nelder–Mead method

    f:目标函数
    mi:变量下界
    ma:变量上界
    max_step:最大迭代次数
    n:多胞体顶点数
    alpha（参考值：1）
    gamma（参考值：2）
    rho（参考值：0.5）
    sigma（参考值：0.5）

*///2019.2.22
template <class fun>
tensor NM(fun f,tensor mi,tensor ma,int max_step,int n,double alpha,double gamma,double rho,double sigma){
    tensor xr,xe,xc,x[n+2];
    /*init*/
    ll d=mi.size();
    rep(i,1,n+1){
        x[i].clear();
        rep(j,0,d-1)x[i].push_back(uniformRand(mi[j],ma[j]));
    }
    rep(it,1,max_step){
        /*sort*/
        if(it==1){
            rep(i,1,n+1)rep(j,1,n)if(f(x[j])>f(x[j+1]))swap(x[j],x[j+1]);
        }
        else{
            int pos=-1;
            double fx=f(x[n+1]);
            rep(i,1,n)if(fx<f(x[i])){
                pos=i;
                break;
            }
            if(pos!=-1){
                tensor temp=x[n+1];
                for(ll i=n+1;i>pos;i--)x[i]=x[i-1];
                x[pos]=temp;
            }
        }
        /*calculate centroid*/
        x[0]=tensor(d,0);
        rep(i,1,n)x[0]+=x[i]/n;
        /*Reflection*/
        xr=x[0]+alpha*(x[0]-x[n+1]);
        if(f(x[1])<=f(xr) && f(xr)<f(x[n])){
            x[n+1]=xr;
            continue;
        }
        /*Expansion*/
        if(f(xr)<f(x[1])){
            xe=x[0]+gamma*(xr-x[0]);
            x[n+1]=f(xe)<f(xr)?xe:xr;
            continue;
        }
        /*Contraction*/
        xc=x[0]+rho*(x[n+1]-x[0]);
        if(f(xc)<f(x[n+1])){
            x[n+1]=xc;
            continue;
        }
        /*Shrink*/
        rep(i,2,n+1)x[i]=x[1]+sigma*(x[i]-x[1]);
    }
    return x[1];
}