梯度下降法详解

在学习、生活、科研以及工程应用中，人们经常要求解一个问题的最优解，通常做法是对该问题进行数学建模，转换成一个目标函数，然后通过一定的算法寻求该函数的最小值，最终寻求到最小值时的输入参数就是问题的最优解。

一个简单的例子，就是求解y=x²的最优解，也就是求当y取得最小值时x的取值。这个问题初中生都会解，谁都知道，直接对函数求导得到导数y=2x，令y=2x=0解得x=0，这就是最优解。然而，很多实际问题的目标函数是很复杂的，很难求出其导数表达式，这种情况下如果还想通过求导数表达式并令其等于0来求最优解就很困难了。不过虽然导数表达式求解困难，但某一确定点的近似导数值还是比较容易求出来的，所以人们通常利用某一确定点的近似导数值来逼近最优解。梯度下降法就是这样的一种算法，在一步一步逼近的过程中，目标函数值呈下降趋势，输入参数的值也一步一步逼近最优解，其整个迭代逼近过程如下图所示：

下面将从数学的角度解释其原理。假设函数f具有n个输入参数：x₁，x₂，…，x_n，第k次逼近的目标函数值可以表达如下：

对第k+1次逼近的目标函数进行泰勒展开，并忽略余项，得到下式：

其中▽f为梯度向量，也即在该点处所有输入参数的偏导数组成的向量，△x为从第k次到第k+1次逼近时输入参数的变化向量。

我们知道，多维函数的偏导数的定义为：

所以可以取一个很小的△x_i值（比如0.0001）来近似计算x_i的偏导数：

向量的点乘可以转换为向量模与向量夹角的相乘：

于是有：

由上式可知，当向量夹角cosθ=-1时向量的点乘结果最小，也即函数值下降最多，此时向量夹角为180度，说明梯度向量与输入参数的变化向量方向相反。所以要使下一次逼近时目标函数值尽可能地相对于当前次降低，也就是使目标函数值下降最快，就应该沿着梯度（偏导数）的负方向寻找下一个逼近点。即：

上式中，α为沿梯度负方向前进的步长，α取值过小会使逼近最优解的速度很慢，从而迭代次数增加，反之如果取值过大，则容易跳过最优解。所以α必须取合适的值才能很好的逼近最优解。通常α有两种取值方法：

1. 刚开始α取一个很大的值，在每一次迭代逼近过程中，如果当前次计算得到的目标函数值大于上次的目标函数值，则减小α步长并重新计算目标函数值，直到当前次目标函数值小于上次目标函数值则停止减小α，并将其作为当前次下降的步长。每一轮迭代逼近的流程图如下图所示：

2. 回溯线性搜索法。首先给α取一个较小的值，先按照一定比例（通常是2倍）增大α，找到第一个使当前次迭代的目标函数值大于上次迭代的目标函数值的α。其次，按照一定比例（通常是0.5倍）减小α，使当前次迭代的目标函数值满足Armijo条件，找到第一个使当前次迭代的目标函数值满足Armijo条件的α，作为本次迭代的步长α。

Armijo条件如下，其中c为0~1之间的固定值，通常取0.3。

对于一维函数：

对于多维函数的某一个输入参数：

对于一维函数，回溯线性搜索法的流程如下图所示（多维函数同理）：

下面我们举了例子，用C++实现梯度下降法来解一个三维函数的最优解。三维函数为：

很明显，其最优解为(2.5, -55.8, 300.25)。不过计算机比较笨，它是不会一眼看出来这个最优解的，但是它的优势在于计算能力强，不管是三维函数还是更复杂的多维函数，它可以通过一定的算法来一步一步逼近这个最优解：

1. 目标函数的实现代码：

//目标函数
float F_fun(float x, float y, float z)
{
  float f = (x-2.5)*(x-2.5) + 5*(y+55.8)*(y+55.8) + 100*(z-300.25)*(z-300.25);
  return f;
}

2. 求输入参数的近似偏导数的代码：

//求x，y，z的偏导数
void G_fun(float x, float y, float z, float &gx, float &gy, float &gz)
{
  float EPS = 1e-4f;   //0.0001


  float f = F_fun(x, y, z);
  float fx = F_fun(x+EPS, y, z);
  float fy = F_fun(x, y+EPS, z);
  float fz = F_fun(x, y, z+EPS);


  gx = (fx - f)/EPS;
  gy = (fy - f)/EPS;
  gz = (fz - f)/EPS;
}

3. 逐渐减小α步长的梯度下降法实现：

int cal_gradient_down(float &x0, float &y0, float &z0)
{
  float x, y, z;
  int max_iter = 50000;   //最多迭代次数
  float g1, g2, g3, s1, s2, s3;
  //输入参数的初始值
  x = x0;
  y = y0;
  z = z0;


  float e = 0.000001;//定义迭代精度
  float ret1 = 0.0;   //保存上一次迭代的目标函数值
  float ret2 = 0.0;    //保存当前次迭代的目标函数值
  int cnt = 0;    //迭代次数计数器
  float alpha = 50000;   //初始α值
  G_fun(x, y, z, g1, g2, g3);    //求偏导数
  float tmp_x, tmp_y, tmp_z;
  
  while (cnt < max_iter)
  {
    tmp_x = x;   //保存上一次的x,y,x坐标
    tmp_y = y;
    tmp_z = z;
    x = x - alpha*g1;   //根据偏导数更新输入参数
    y = y - alpha*g2;
    z = z - alpha*g3;


    ret1 = F_fun(tmp_x, tmp_y, tmp_z);  //计算上一次迭代的目标函数值
    ret2 = F_fun(x, y, z);    //计算当前次迭代的目标函数值


    if (ret2 > ret1)   //如果当前次目标函数值大于上一次目标函数值，则减小α并重新计算当前次目标函数值
    {
      alpha *= 0.9;
      x = tmp_x;    //将上一次的x,y,x坐标赋值给当前的x,y,x坐标
      y = tmp_y;
      z = tmp_z;
      continue;
    }


    printf("f = %f\n", ret2);


    if(abs(ret2 - ret1) < e)   //如果2次迭代的结果变化很小，则认为获取到最优解，并结束迭代
    {
      x0 = x;    
      y0 = y;
      z0 = z;
      return 0;
    }


    G_fun(x, y, z, g1, g2, g3);    //更新偏导数
    cnt++;    //计数器加1
  }


  return -1;
}

运行上述代码，求得最优解(2.506215, -55.800034,  300.249939)，可以看到其与真实解还是非常接近的。

4. 回溯线性搜索法确定α步长的梯度下降法实现：

首先是回溯线性搜索法确定α步长代码：

void GetA_Armijo(float x, float y, float z, float g1, float g2, float g3, float &alpha1, float &alpha2, float &alpha3)
{
  float c1 = 0.3;
  float now = F_fun(x, y, z);
  float next1 = F_fun(x - alpha1*g1, y, z);
  float next2 = F_fun(x, y - alpha2*g2, z);
  float next3 = F_fun(x, y, z - alpha3*g3);


  //求输入参数x的步长alpha1
  int count =30;
  while (next1 < now)
  {
    alpha1 *= 2.0;
    next1 = F_fun(x - alpha1*g1, y, z);
    count--;
    if (count == 0)
      break;
  }


  count = 30;
  while (next1 > now - alpha1*c1*g1*g1)
  {
    alpha1 *= 0.5;
    next1 = F_fun(x - alpha1*g1, y, z);
    count--;
    if (count == 0)
      break;
  }


  //求输入参数y的步长alpha2
  count = 30;
  while (next2 < now)
  {
    alpha2 *= 2.0;
    next2 = F_fun(x, y - alpha2*g2, z);
    count--;
    if (count == 0)
      break;
  }


  count = 30;
  while (next2 > now - alpha2*c1*g2*g2)
  {
    alpha2 *= 0.5;
    next2 = F_fun(x, y - alpha2*g2, z);
    count--;
    if (count == 0)
      break;
  }


  //求输入参数z的步长alpha3
  count = 30;
  while (next3 < now)
  {
    alpha3 *= 2.0;
    next3 = F_fun(x, y, z - alpha3*g3);
    count--;
    if (count == 0)
      break;
  }


  count = 30;
  while (next3 > now - alpha3*c1*g3*g3)
  {
    alpha3 *= 0.5;
    next3 = F_fun(x, y, z - alpha3*g3);
    count--;
    if (count == 0)
      break;
  }
}

上述求得步长之后，则可以直接进行梯度下降：

int cal_gradient_down(float &x0, float &y0, float &z0)
{
  float x, y, z;
  int max_iter = 50000;   //最多迭代次数
  float g1, g2, g3, s1, s2, s3;
  x = x0;
  y = y0;
  z = z0;


  float e = 0.000001;//定义迭代精度
  float ret_1 = 0.0;
  float ret0 = 11111.0;
  float ret1 = 666775.0;
  float ret2 = 12.0;
  int cnt = 0;


  float alpha1 = 200000;
  float alpha2 = 200000;
  float alpha3 = 200000;


  float tmp_x, tmp_y, tmp_z;
 
  while (cnt < max_iter)
  {
    G_fun(x, y, z, g1, g2, g3);    //求偏导数
    //回溯线性法求步长
    GetA_Armijo(x, y, z, g1, g2, g3, alpha1, alpha2, alpha3);
    //更新输入参数
    x = x - alpha1*g1;
    y = y - alpha2*g2;
    z = z - alpha3*g3;




    ret_1 = ret0;
    ret0 = ret1;
    ret1 = ret2;
    ret2 = F_fun(x, y, z);


    if (abs(ret0 - ret_1) < e && abs(ret1 - ret0) < e && abs(ret2 - ret1) < e)//如果2次迭代的结果变化很小，结束迭代
    {
      x0 = x;    
      y0 = y;
      z0 = z;
      return 0;
    }


    printf("f = %f\n", ret2);


    cnt++;
  }


  return -1;
}

运行上述代码得到最优解(2.500022, -55.799988, 300.250000)。

对比上述结果可以看到，减少步长的方法需要469次迭代，而回溯线性法需要的迭代次数相比减少步长的方法少得多，只需要10次迭代就能找到最优解。所以后者得收敛速度要比前者快得多。