最优化 | 无约束优化方法 | C++实现

参考资料

• https://blog.csdn.net/tangshishe/article/details/116670314
• 无约束优化方法
• 模拟退火算法
• 基于matlab模拟退火算法求解函数极值问题

1. 前言

无约束问题是指只有优化目标，而不存在约束条件的问题，用数学模型可表示为
$$\underset{x\in R^n}{\min }f(x)$$

通常假设$f(x)$具有二阶连续偏导数。

无约束优化方法是通常是由给定的初始点出发，按一定的规则不断向最优解趋近。在前文最优化 | 一维搜索与方程求根中提到的一维搜索方法、梯度下降法、牛顿法等都可以用来解决无约束优化问题。

无约束最优化方法的核心问题是选择搜索方向。

下面主要介绍以及实现几种常见的无约束优化方法。

所涉及c++代码均存于github仓库。
部分python版本代码见于github仓库。

2. 梯度下降法

2.1 原理

梯度下降法的优化思想是用当前位置负梯度方向作为搜索方向，因为该方向为当前位置的最快下降方向，所以也被称为是最速下降法。

基本的算法为:

在算法中，迭代收敛可以用梯度向量长度小于某个预定值， 或者两次迭代间函数值变化小于某个预定值来判断。迭代中需要用某种一维搜索方法求步长，不过实际应用时，绝大部分是使用预先确定的步长，所以步长也叫做学习率。

在适当正则性条件下， 最速下降法可以收敛到局部极小值点并且具有线性收敛速度。 但是， 最速下降法连续两次的下降方向$-\nabla f(x^{(t)})$和$-\nabla f(x^{(t+1)})$是正交的， 这是因为在第$t$步沿$-\nabla f(x^{(t)})$搜索时找到点时已经使得函数值在该方向上不再变化， 下一个搜索方向如果与该方向不正交就不是下降最快的方向。 这样，最速下降法收敛速度比较慢。 所以， 在每次迭代的一维搜索时， 可以不要求找到一维搜索的极小值点， 而是函数值足够降低就结束一维搜索。

2.2 c++实现

• 求函数$f(x)=x^2-4x$的极值点。

  //
  // Created by CHH3213 on 2022/9/12.
  //
  #include<iostream>
  #include <valarray>
  /*
   * 使用梯度下降法求解函数极值
   */
  using namespace std;
  
  #define DELTA 1e-4
  #define EPS 1e-6
  
  double func(double x){
      /**
       * 函数
       */
      return x*x-4*x;
  }
  
  
  double func_prime(double x){
      /**
       * 函数导数
       */
      return (func(x+DELTA)-func(x-DELTA))/(2*DELTA);
  }
  
  
  
  double GradientDescent(double x0,double learning_rate ){
      /**
       * 梯度下降算法
       * x0为迭代初值
       * learning_rate为学习率
       */
      double x = x0;
      double x_last = 0.5;
      while(abs(x-x_last)>EPS){
          x_last=x;
          x = x_last-learning_rate* func_prime(x_last);
      }
      return x;
  }
  int main(){
  
      double x1= GradientDescent(0,0.01);
      cout<<x1<<endl;
  }
  
  

• 求解$f(x)=4(x_1-1{)}^2+(x_2-2{)}^4$的极值点，假设从初始点$[0,0]$出发。（显然，函数具有全局极小值点(1,2)。）
  非常容易求得
  $$\nabla f(x)=[8(x_1-1),4(x_2-2{)}^3{]}^T$$
  然后使用梯度更新公式即可，和一元函数差别不大（后面举例均已一元函数为主）。

  c++代码实现如下：

  //
  // Created by CHH3213 on 2022/9/12.
  //
  #include<iostream>
  #include <valarray>
  #include <vector>
  /*
   * 使用梯度下降法求解多元函数极值：f(x)=4(x1-1)^2+(x2-2)^4
   */
  using namespace std;
  
  #define DELTA 1e-4
  #define EPS 1e-9
  
  double func(double x1, double x2){
      /**
       * 函数
       */
      return 4*pow(x1-1,2)+pow(x2-2,4);
  }
  
  
  vector<double> func_prime(double x1,double x2){
      /**
       * 函数导数
       */
       double nabla_x1 = (func(x1+DELTA,x2)-func(x1-DELTA,x2))/(2*DELTA);//对x1的偏导数
       double nabla_x2 = (func(x1,x2+DELTA)-func(x1,x2-DELTA))/(2*DELTA);//对x2的偏导数
      return {nabla_x1,nabla_x2};
  }
  
  
  
  vector<double> GradientDescent(vector<double> x0,double learning_rate ){
      /**
       * 梯度下降算法
       * x0为迭代初值
       * learning_rate为学习率
       */
      vector<double> x = x0;
      vector<double> x_last = {0.5,0.5};
      while(abs(x[0]-x_last[0])>EPS||abs(x[1]-x_last[1])>EPS){
          x_last=x;
          vector<double>prime = func_prime(x_last[0],x_last[1]);
          x[0] = x_last[0]-learning_rate*prime[0] ;
          x[1] = x_last[1]-learning_rate*prime[1] ;
  //        cout<<x[0]<<' '<<x[1]<<endl;
      }
      return x;
  }
  int main(){
  
      vector<double> x= GradientDescent({0,0},0.01);
      cout<<x[0]<<' '<<x[1]<<endl;
  }
  
  

2.3 共轭梯度法

最速下降法收敛速度较慢，而且越靠近最小值点步长越短，速度越慢。对于一个函数曲面为狭长的山谷形 状的目标函数，为了达到山谷的最低点，会沿相互正交的梯度方向反复跳跃。

一种解决的方法是前面所说的不使用最优一维搜索而是使用预先确定的步长或步长序列，还有一种方法是将两步之间的梯度正交， 推广到两步之间为“共轭”关系，即对某个矩阵 $A_t$ 使得
$${\left[{\mathbf{g}}^{(t)}\right]}^T{A}_t{\mathbf{g}}^{(t+1)}=0.$$
这种方法的来源是对目标函数 $f(\mathbf{x})$ 用二次多项式近似（泰勒展开，忽略高次项）：
$$f(\mathbf{x})\approx \frac{1}{2}{\mathbf{x}}^{T}A\mathbf{x}+{\mathbf{b}}^T\mathbf{x}+c,$$
而为了将一个 $k$元的凸二次多项式达到最小值点，只要从负梯度方向出发，沿着 $k$ 个相互都关于 $A$ 共轭的方向进行一维搜索就保证达到最小值点。连续可微函数在最小值点附近一般都很接近于一个凸二次多项式函数， 所以共轭梯度法能够克服最速下降法在靠近最小值点时收敛变慢的问题。

$A_t$ 矩阵选取最好是Hessian阵 ${\nabla }^{2}f\left({\mathbf{x}}^{(t)}\right)$ ， 但Hessian阵一般比梯度更难获得。共轭梯度法实际的算法并不去求 $A_t$ ，而是从 ${\mathbf{x}}^{(1)}$ 的负梯度出发， 然后每次用当前 ${\mathbf{x}}^{(t)}$ 处的负梯度和 上一步的共轭方向的线性组合作为 $t$ 到 $t+1$ 步的搜索方向。即：
$$\begin{array}{rl}& {\mathbf{d}}^{(1)}=-{\mathbf{g}}^{(1)};\\ & {\mathbf{d}}^{(t)}=-{\mathbf{g}}^{(t)}+{\beta }^{(t)}{\mathbf{d}}^{(t-1)},t=2,3,\dots .\end{array}$$
对 $k$元凸二次多项式，其中的系数 ${\beta }^{(t)}$ 可以用 $A$ 解出精确公式， 使得各个 ${\mathbf{d}}^{(t)}$ 彼此共轭。但是对一般标函 数 $f(\mathbf{x})$ ，只能使用一些近似的 ${\beta }^{(t)}$ ，方法如：

• Fletcher-Reeves: ${\beta }^{(t)}=\max \left(0,\frac{{\left[{\mathbf{g}}^{(t)}\right]}^T{\mathbf{g}}^{(t)}}{{\left[{\mathbf{g}}^{(t-1)}\right]}^T{\mathbf{g}}^{(t-1)}}\right)$.
• Polak-Ribière: ${\beta }^{(t)}=\frac{{\left[{\mathbf{g}}^{(t)}\right]}^T\left[{\mathbf{g}}^{(t)}-{\mathbf{g}}^{(t-1)}\right]}{{\left[{\mathbf{g}}^{(t-1)}\right]}^T{\mathbf{g}}^{(t-1)}}$.

Fletcher-Reeves共轭梯度法在极小值点附近可以确保收敛。

3. 牛顿法

3.1 原理

牛顿法利用Hessian阵和梯度向量求下一步，不需要一维搜索，对于二次多项式函数可以一步得到最小值点。 许多光滑的目标函数在靠近极小值点的邻域内可以很好地用凸二阶多项式逼近。对一个存在二阶连续偏导的多元函数 $f(\mathbf{x})$, 有如下的泰勒近似
$$f(\mathbf{x})\approx f\left({\mathbf{x}}^{(0)}\right)+\nabla f{\left({\mathbf{x}}_0\right)}^T\left(\mathbf{x}-{\mathbf{x}}^{(0)}\right)+\frac{1}{2}{\left(\mathbf{x}-{\mathbf{x}}^{(0)}\right)}^T{\nabla }^{2}f\left({\mathbf{x}}^{(0)}\right)\left(\mathbf{x}-{\mathbf{x}}^{(0)}\right),\text{\hspace{1em}(1)}$$
若Hessian阵 ${\nabla }^{2}f\left({\mathbf{x}}^{(0)}\right)$ 为正定矩阵，上式右侧的二次多项式函数的最小值点在
$${\mathbf{x}}^{\ast }={\mathbf{x}}^{(0)}-{\left[{\nabla }^{2}f\left({\mathbf{x}}^{(0)}\right)\right]}^{-1}\nabla f\left({\mathbf{x}}^{(0)}\right)\text{\hspace{1em}(2)}$$
处达到。所以牛顿法取初值 ${\mathbf{x}}^{(0)}$ 后，用公式
$${\mathbf{x}}^{(t+1)}={\mathbf{x}}^{(t)}-{\left[{\nabla }^{2}f\left({\mathbf{x}}^{(t)}\right)\right]}^{-1}\nabla f\left({\mathbf{x}}^{(t)}\right),t=0,1,2,\dots \text{\hspace{1em}(3)}$$
进行迭代，直至收敛。收敛的判断准则可以取为 $\Vert \nabla \mathbf{f}(x)\Vert$ 足够小，或者取为 $\left|f\left({\mathbf{x}}^{(t+1)}\right)-f\left({\mathbf{x}}^{(t)}\right)\right|$ 足够 小。

公式(3)在算法实现时，应将逆矩阵表示改为如下的线性方程组求解:
$$\begin{gathered} {\nabla }^{2}f\left({\mathbf{x}}^{(t)}\right){\mathbf{d}}^{(t)}=-\nabla f\left({\mathbf{x}}^{(t)}\right), \\ {\mathbf{x}}^{(t+1)}={\mathbf{x}}^{(t)}-{\mathbf{d}}^{(t)}\text{\hspace{1em}(4)} \end{gathered}$$
如果初值接近于最小值点并且目标函数满足一定正则性条件，牛顿法有二阶收敛速度。

牛顿法的性能对Hessian阵 $H(\mathbf{x})$ 十分依赖。

• 如果 $H\left({\mathbf{x}}^{(t)}\right)$ 不满秩，则公式(4)的解不唯一，算法无法进行；
• 如果 $H\left({\mathbf{x}}^{(t)}\right)$ 很接近于不满秩，则 $H{\left({\mathbf{x}}^{(t)}\right)}^{-1}\nabla f\left({\mathbf{x}}^{(t)}\right)$ 的计算很不稳定，而且会将下一个近似点 ${\mathbf{x}}^{(t+1)}$ 远离 ${\mathbf{x}}^{(t)}$ ，这很可能会使得逼近失败。

另外，从公式(1)看出，Hessian阵也是函数的局部性质的重要特征。

• 如果 ${\mathbf{x}}_{\ast }$ 是一个稳定点， 即 $\nabla f\left({\mathbf{x}}_{\ast }\right)=0$ ， Hessian阵与此处目标函数的曲面形状密切相关。
• 如果 $H\left({\mathbf{x}}_{\ast }\right)$ 是正定的， 则 ${\mathbf{x}}_{\ast }$ 是局部极小值点。
• 如果 $H\left({\mathbf{x}}_{\ast }\right)$ 是负定的，则 ${\mathbf{x}}_{\ast }$ 是局部极大值点。如果 $H\left({\mathbf{x}}_{\ast }\right)$ 既有正特征值又有负特征值，则 ${\mathbf{x}}_{\ast }$ 是目标函数的鞍点，不是局部极值点。

使用牛顿法求最小值时，最好的情况是Hessian阵始终为正定阵， 这时函数是严格凸函数，极小值点也是全局最小值点。在线性最小二乘问题中，目标函数的Hessian阵是 $X^{T}X$ 矩阵。

如果目标函数非负定但是接近于不满秩，可以考虑给 $H\left({\mathbf{x}}^{(t)}\right)$ 加一个对角的正定阵，例如，岭回归就是将 $X^{T}X$ 替换成 $X^{T}X+\lambda I$ 。

3.2 c++实现

求函数$f(x)=x^2-4x$的极值点。

//
// Created by CHH3213 on 2022/9/12.
//
#include<iostream>
#include <valarray>

using namespace std;

#define DELTA 1e-4
#define EPS 1e-6
/*
 * 使用牛顿法求函数极值点
 */
double func(double x){
    /**
     * 构造函数
     * x:自变量
     * return：函数值
     */
    return x*x-4*x;
}

double func_prime(double x){
    /**
     * 计算函数导数,以导数定义的方式求解
     */
    return (func(x+DELTA)-func(x-DELTA))/(2*DELTA);
}
double func_prime_prime(double x){
    /**
     * 计算函数二阶导,以导数定义的方式求解
     */
    return (func_prime(x+DELTA)-func_prime(x-DELTA))/(2*DELTA);
}


double Newton(double x0){
    double x = x0;
    double x_last=0.1;
    while(abs(x-x_last)>EPS){
        x_last=x;
        x=x_last-func_prime(x_last)/ func_prime_prime(x_last);
    }
    return x;
}

int main(){
    cout<<"extreme value x: "<<Newton(1)<<endl;
    return 0;
}

4. 模拟退火算法

4.1 原理

退火是一种机械加工技术，加热金属零件到一定温度后进行加工，然后控制温度缓慢地下降到常温。在高 温时，分子运动速度较快， 在缓慢降温时可以调整机械应力、减少缺陷。

模拟退火算法的名称是利用了“缓慢降温”的做法和“温度”参数，这种算法利用随机性进行全局优化，在目标函数定义域进行随机游动，试图找到全局最小值点。

算法从当前位置 ${\mathbf{x}}^{(t)}$ 按照某种抽样分布生成随机数 ${\mathbf{ϵ}}^{(t)}$ ，令 ${\mathbf{x}}^{\prime }={\mathbf{x}}^{(t)}+{\mathbf{ϵ}}^{(t)}$ ， 记 $\Delta y=f\left({\mathbf{x}}^{\prime }\right)-f\left({\mathbf{x}}^{(t)}\right)$ ， 是否接受${\mathbf{x}}^{\prime }$作为下一步迭代的解取决于下面计算的概率 :
$$p^{(t+1)}=\{\begin{array}{ll}1,& \text{ 若 }\Delta y\le 0,\\ \exp \left(-\Delta y/T^{(t+1)}\right),& \text{ 若 }\Delta y>0.\end{array}$$

换句话说，如果生成的解 ${\mathbf{x}}^{\prime }$的函数值比前一个解的函数值更小，则接受${\mathbf{x}}^{(t+1)}={\mathbf{x}}^{\prime }$作为一个新解。否则以概率$\exp \left(-\Delta y/T^{(t+1)}\right)$接受${\mathbf{x}}^{\prime }$作为一个新解。否则令 ${\mathbf{x}}^{(t+1)}={\mathbf{x}}^{(t)}$ ，停留不动，等待下一次随机转移。

这里 $T^{(t+1)}$ 是转移概率的参数， 类似于退火时 的温度， $T^{(t+1)}$ 越大，则 $p^{(t+1)}$ 越大，逆向的随机跳转越频繁，可以用来从局部极小值点跳开；

随着迭代的进行， $T^{(t+1)}$ 应该逐步变小，使得逆向随机跳转变稀少， 这样迭代能最终收敛。

$T^{(t+1)}$ 应选取为趋于 0 的数列， 在一定条件下如果取
$$T^{(t)}=T^{(0)}\frac{\ln 2}{\ln (t+1)},$$
可以确保迭代收收敛到全局最小值点， 但显然这个速度太慢了。实际应用中往往用
$$T^{(t+1)}=\gamma T^{(t)},$$
其中 $0<\gamma <1$ ， 或取
$$T^{(t)}=T^{(0)}/t$$

4.2 c++实现

求函数$f(x)=x^2-4x$的极值点。

//
// Created by CHH3213 on 2022/9/12.
//
#include<iostream>
#include <valarray>

using namespace std;

#define DELTA 1e-4
#define EPS 1e-6
/*
 * 使用模拟退火求函数极值点
 */
double func(double x){
    /**
     * 构造函数
     * x:自变量
     * return：函数值
     */
    return x*x-4*x;
}


double rnd(double low,double high){
    /**
     * 返回指定范围内的随机浮点数
     */
    double rd=rand()/((double)RAND_MAX+1.0);
    return low+rd*(high-low);
}

double Annealing(double T0, double T_min, double x0,double gamma){
    /**
     * T0:初始温度
     * T_min:温度的下限，若温度T达到T_min，则停止搜索
     * x0:初始迭代值
     * gamma:T_{k+1} = gamma * T_k, 0<gamma<1
     * return:极值点
     */
     double T=T0;
     double x = x0;
     double x_ = 2;
     double delta_y=1.0;
     while(abs(delta_y)>EPS){//如果函数值满足精度，停止搜索
//     while(T>T_min){//温度T达到T_min，则停止搜索

         double x_ = x+rnd(-1,1);
         cout<<x_<<endl;
         delta_y = func(x_)-func(x);
         if(delta_y<=0){
             x = x_;
         }else{
             if(exp(-delta_y/T)>rnd(0,1)){
                 x=x_;
             }
         }
         T=gamma*T;
         cout<<x<<endl;
     }
    return x;
}

int main(){
    double T0=200,T_min = 1,x0=2,gamma=0.9;
    cout<<"extreme  value: "<<Annealing(T0,T_min,x0,gamma)<<endl;
    return 0;
}

5. 遗传算法

模拟退火算法是利用随机移动来寻找全局最优点。遗传算法是另一类算法的典型代表，它产生多个搜索点， 利用这些搜索点以及搜索点之间的交互去接近全局最优点。

遗传算法利用了生物学中择优配对、交叉遗传、遗传突变的思想使得整个种群（搜索点集合）不断进化 (接 近全局最优点）。设要求 ${\min }_{\mathbf{x}\in {\mathbb{R}}^d}f(\mathbf{x})$ ， 算法的核心包括：

• 初始化。可以预先获得包含最优点的范围，比如一个超长方体，然后在其中均匀产生 $m$ 个初始的搜索 点 ${\mathbf{x}}^{(i)},i=1,2,\dots ,m$ 。也可以用多元正态分布抽样，方差阵用对角阵。也可以用重要抽样思想， 使用对于可能包含最优点的区域的先验知识构造初始搜索点的抽样分布。每个搜索点称为一个“染色 体”。
• 每一步先选出 $m$ 对父本，从父本用交叉遗传产生下一代的搜索点。在选父本时，要体现择优配对的思 想，比如，仅从 $m$ 个现有搜索点中的 $k$ 个表现最优的点（目标函数值最小）中抽取父本; 或者，令被抽 取为父本的概率，与目标函数值反向相关，等等。
• 对选出的 $m$ 对父本中的每一对，如 ${\mathbf{x}}^{(i)},{\mathbf{x}}^{(j)}$ ，进行染色体交叉，获得一个子代的搜索点 (染色体)。 交叉的方法有多种，比如 $d$ 维的每一分量随机地从两个父本对应分量中任选一个。
• 经过交叉遗传获得下一代染色体后，如果没有别的步骤，这个种群（搜索点集合）就会仅在有限个点 上活动，不可能搜索全定义域。所以，还要有一个“突变”步骤，使得搜索点能够转移到附近的位置。 一般是选择对每个搜索点加多元标准正态噪声乘以一个小的步长。

这种利用多个搜索点的方法还有许多， 例如“粒子蜂拥”方法，记录每一个 搜索点的位置、速度、历史最优点，更新速度为倾向于自己的历史最优点和群体的历史最优点。“茧火虫”方 法，是每个搜索点都逐次被每一个优于自己的点吸引。

这些利用群体随机变化寻找全局最优的方法有利用查找比较广泛的定义域，但是收敛到最小值点比较慢， 所以可以配合一些局部搜索算法；可以对每一个搜索点都作局部搜索，并更新搜索点和目标函数值；也可以对每一个搜索点都作局部搜索后，仅更新其搜索得到的局部更优的目标函数值，但并不更新搜索点本 身，这种做法更能避开局部最小点的吸引。