标签：梯度 long eta 下降 dx double part

本文算是对上次写的题解「洛谷P2571 [SCOI2010]传送带」中讲到的梯度下降法的整理吧。。。

非 \(O(1)\) 复杂度求解多元函数最值的方法有很多：粒子群算法、模拟退火、三分套三分、牛顿迭代法……

在此介绍梯度下降法。

梯度

了解多元微积分的各位大佬们都知道，梯度是一个向量，指向多元函数的值变化最快的方向，大小即为变化率。

「不了解的可以看看 3Blue1Brown 的作者 Grant Sanderson 讲的多元微积分课程，在b站上已有熟肉【链接】」

百度百科对梯度的定义：

设二元函数\(z=f(x,y)\) 在平面区域\(D\)上具有一阶连续偏导数，则对于每一个点\(P(x,y)\)都可定出一个向量

\[\left\{ \frac{\part f}{\part x},\frac{\part f}{\part y} \right\}=f_x(x,y)\hat{i}+f_y(x,y)\hat{j} \]

该函数「也不是不对，但称为向量更好理解」就称为函数\(z=f(x,y)\)在点\(P(x,y)\)的梯度，记作\(\text{grad}f(x,y)\)或\(\nabla f(x,y)\)。

于是有：

\[\text{grad}f(x,y)=\nabla f(x,y)=\left\{ \frac{\part f}{\part x},\frac{\part f}{\part y} \right\}=f_x(x,y)\hat{i}+f_y(x,y)\hat{j} \]

其中

\[\nabla=\frac{\part}{\part x}\hat{i}+\frac{\part}{\part y}\hat{j} \]

称为（二维的）向量微分算子或Nabla算子，且有

\[\nabla f=\frac{\part f}{\part x}\hat{i}+\frac{\part f}{\part y}\hat{j} \]

梯度下降法的主要思想

再来介绍梯度下降法。

正如我们先前对 \(f(x,y)\) 求偏导，令偏导值为 \(0\)，解出所有驻点的操作，实际上就是试图直接令梯度值为 \(0\)。

但是也正如我们遇到的问题，有时直接求解比较困难，我们不妨通过迭代，让梯度值逐渐下降，即逐渐接近极值点。

例如，对于函数 \(f(x)=x^2\)，求导得 \(f'(x)=2x\)。

不妨设 \(x_0=1\)，则从点 \((1,1)\) 开始，计算其梯度值：\(f'(1)=2>0\)。

梯度值 \(>0\)，说明在该点沿使 \(x\) 增大的方向「即 \(x\) 轴正半轴方向」，\(f(x)\) 函数值会增大，增大的速率为 \(2\)。

容易发现，在 \((3,9)\) 处，\(f'(3)=6\)，增大的速率为 \(6\)，观察图像也可得知，沿 \(x\) 轴正半轴方向，\(f(x)\) 的图像变陡，上升速率变快。

而在 \((0,0)\) 处，\(f'(0)=0\)，函数图像在这一点的切线斜率是平的，增大的速率无限接近 \(0\)，此时梯度值为 \(0\)，\((0,0)\) 是函数 \(f(x)\) 的一个驻点，同时 \(f(x)\) 也在此取到最小值。

因此，我们有如下策略：

梯度值 \(>0\)，向左移动一点；梯度值 \(<0\)，向右移动一点。

即每次让 \(x\) 朝与梯度值符号相反的方向移动，使梯度值逐渐下降，最终趋于 \(0\)。

用数学语言来描述，即：

\[x_{k+1}=x_k-\nabla f(x_k) \]

其中 \(x_k\) 为第 \(k\) 次迭代时点的横坐标，\(x_{k+1}\) 为第 \(k+1\) 次迭代移动到的点的横坐标，\(x_0\) 表示初始横坐标。

\(\nabla f(x)\) 表示在 \(x\) 处的梯度值，\(\nabla f(x)=\frac{\text df}{\text dx}\)。

\(x_k-\nabla f(x_k)\) 表示向与梯度值符号相反的方向移动，符合之前的策略

优化——「学习率」

但这样就有一个问题：

例如 \(f(x)=x^2,x_0=10\)，则 \(\nabla f(x)=2x\)：

\[\begin{aligned} x_1&=x_0-2\cdot x_0=-10 \\ x_2&=x_1-2\cdot x_1=10 \\ x_3&=x_2-2\cdot x_2=-10 \\ ... \end{aligned} \]

可以看到，虽然横坐标一直在变化，但一直在 \(10、-10\) 之间振荡，函数值始终没有降到最低点。

为此，我们需要一个参数控制移动的距离，这个参数被称作学习率，用 \(\eta\) 表示：

\[x_{k+1}=x_k-\eta\cdot\nabla f(x_k) \]

显然，在刚刚的计算中，\(\eta=1\)。

当 \(\eta>1\) 时，函数值不仅不会降到最低点，甚至会越来越大：

\(e.g.\qquad\eta=1.05\)

\[\begin{aligned} x_1&=x_0-2\cdot x_0=-11 \\ x_2&=x_1-2\cdot x_1=12.1 \\ x_3&=x_2-2\cdot x_2=-13.31 \\ ... \end{aligned} \]

而若将 \(\eta\) 调低到一个较小值，如 \(0.02\)：

\[x_1=x_0-2\cdot x_0=9.6 \\ x_2=x_1-2\cdot x_1=9.2 \\ x_3=x_2-2\cdot x_2=8.8 \\ ... \]

降低的速度太慢，容易超时。

当 \(\eta=0.2\) 时，迭代 \(10\) 次左右即可降入谷底：

\[\begin{aligned} x_1&=x_0-2\cdot x_0=6&||\nabla f(x_0)||=12 \\ x_2&=x_1-2\cdot x_1=3.6&||\nabla f(x_1)||=7.2 \\ x_3&=x_2-2\cdot x_2=2.16&||\nabla f(x_2)||=4.32 \\ ... \end{aligned} \]

可以看到，右边的梯度值在每次迭代后都会下降，故称为梯度下降法「Gradient Descent」

只要选择合适的学习率，梯度就可以下降到任意小：

\[\lim_{k\to\infty}f(x_k)=\min f(x) \]

可以用泰勒公式进行严格证明，此处不再赘述（逃

因此可以通过梯度值的大小作为终止条件「也可以直接用迭代次数控制精度」

对于二元函数，同样可以用梯度下降法求解极值：

\[f(x_{k+1},y_{k+1})=f(x_k,y_k)-\eta\cdot\nabla f(x_k,y_k) \]

\(e.g.\qquad f(x,y)=x^2+2y^2,(x_0,y_0)=(-3.5,-3.5),\eta=0.1\)，则\(\nabla f(x,y)=(2x,4y)\)

\[(x_1,y_1)=(x_0,y_0)-\eta\cdot\nabla f(x_0,y_0)=(-2.8,-2.1) \\ (x_2,y_2)=(x_1,y_1)-\eta\cdot\nabla f(x_1,y_1)=(-2.24,-1.26) \\ ... \]

代码——「梯度下降法」

code：

#include<iostream>
#include<cstdio>
#include<cmath>

using namespace std;

inline long double f(long double x)
{
    return 3*x*x*x*x-x*x*x+2*x*x-9*x+5*sqrt((x+3)*(x+3)+(5*x+6)*(5*x+6))-25;
}

inline long double numerical_diff(long double x)//数值微分法估计一阶导数
{
    long double dx=1e-6;
    return (f(x+dx)-f(x-dx))/(dx*2);
}

inline long double gradient_descent(long double x,long double eta)//梯度下降法
{
    long double eps=1e-6;
    //int cnt=0;
    while(abs(numerical_diff(x))>eps)
    {
        x-=eta*numerical_diff(x);
        //cnt++;
    }
    //cout<<cnt<<endl;
    return f(x);
}

signed main()
{
    cout<<gradient_descent(0,0.02)<<endl;
    return 0;
}

\(e.g.\qquad f(x)=3x^4-x^3+2x^2-9x+5\sqrt{(x+3)^2+(5x+6)^2}-25\)

经过\(9\)次迭代达到目标精度。

优化——动量梯度下降法「MGD」

需要指出的是，梯度值为 \(0\) 的驻点不一定是函数的极值点，如 \(f(x)=x^3\) 在 \(x=0\) 处梯度值为 \(0\)，但并不是函数的极值点：

「实际上，\((0,0)\) 是 \(f(x)\) 的拐点，此处函数的凹凸性发生改变」

同时，对于非单峰函数来说，梯度下降法的结果易受到初始值的影响，也就是陷入局部最优解。

现实生活中，一个小球从高处落下，大概率会越过比较低的坎继续下降。

我们同样可以引入惯性来优化：

用梯度模拟受力，使之不直接控制移动距离，而是给小球一个速度 \(v\)；同时可以引入阻力，也就是速度衰减率 \(\beta\) 使小球减速。

于是，整个过程就像一个有动量的小球在空间中来回滚动，故称动量梯度下降法「Gradient Descent with Momentum，简称MGD」

根据动量定理

\[F\Delta t=m\Delta v \]

移项，得：

\[\Delta v=\frac{\Delta t}{m}\cdot F \]

直接令 \(\eta=\frac{\Delta t}{m}\)，则：

\[\Delta v=\eta\cdot F \]

但在实际编程中，不必那么精细地刻画阻力，直接用一个速度衰减率 \(\beta\) 代替即可：

\[v_{k+1}=\beta\cdot v_k-\eta\cdot\nabla f(x_k) \]

code：

#include<iostream>
#include<cstdio>
#include<cmath>

using namespace std;

inline long double f(long double x)
{
    return 3*x*x*x*x-x*x*x+2*x*x-9*x+5*sqrt((x+3)*(x+3)+(5*x+6)*(5*x+6))-25;
}

inline long double numerical_diff(long double x)
{
    long double dx=1e-6;
    return (f(x+dx)-f(x-dx))/(dx*2);
}

inline long double mgd(long double x,long double eta,long double beta)//动量梯度下降法
{
    long double v=0,eps=1e-6;
    int cnt=0;
    while(abs(numerical_diff(x))>eps)
    {
        v=beta*v-eta*numerical_diff(x);
        x+=v;
        cnt++;
    }
    cout<<cnt<<endl;
    return f(x);
}

signed main()
{
    cout<<mgd(0,0.03,0.05)<<endl;
    return 0;
}

迭代次数：\(13\)。

「多峰函数的优化效果较明显」

优化——自适应梯度算法「AdaGrad」

对比一下两种算法的真实移动距离 \(\Delta x_i\)：

	梯度下降法	动量梯度下降法
\(\Delta x_i\)	\(\eta\cdot\frac{\part f}{\part x_i}\)	\(\eta\cdot v_i\)

动量梯度下降法只优化了 \(\frac{\part f}{\part x_i}\) 的部分。

很自然的想到，是否可以优化学习率 \(\eta\)？

答案是肯定的，这种算法称为自适应梯度算法「AdaGrad【Adapt+Gradient】」，基本思想是开始时快速移动接近目标，然后减速提高精度「梯度对其产生的动力效果逐渐减弱」。实现时记录每次梯度下降的梯度平方和 \(h\)，将 \(\eta\) 除以 \(\sqrt{h}\) 进行调整「为防止分母出现 \(0\)，往往加上一个较小的常量」

code：

#include<iostream>
#include<cstdio>
#include<cmath>

using namespace std;

inline long double f(long double x)
{
    return 3*x*x*x*x-x*x*x+2*x*x-9*x+5*sqrt((x+3)*(x+3)+(5*x+6)*(5*x+6))-25;
}

inline long double numerical_diff(long double x)
{
    long double dx=1e-6;
    return (f(x+dx)-f(x-dx))/(dx*2);
}

inline long double AdaGrad(long double x,long double eta,long double beta)//自适应梯度算法
{
    long double v=0,h=0,grad=numerical_diff(x),eps=1e-6;
  	//int cnt=0;
    while(abs(grad)>eps)
    {
        h+=grad*grad;
        v=beta*v-eta*grad/(sqrt(h)+eps);
        x+=v;
      	grad=numerical_diff(x);
        //cnt++;
    }
    //cout<<cnt<<endl;
    return f(x);
}

signed main()
{
    cout<<AdaGrad(0,0.5,0.05)<<endl;
    return 0;
}

迭代次数：\(12\)。

「多峰函数的优化效果较明显」

优化——自适应动量算法「Adam」

自适应梯度算法的缺陷：若初始点离极值点较远，可能还没到极值点，小球就跑不动了「相应的有 RMSProp 优化算法，通过把 \(h\) 乘以一个衰减率 \(\beta_2\) 来逐步遗忘之前的梯度，专业一点说是“指数移动平均”，与动量梯度下降法有异曲同工之妙，此处不作详细介绍」

我们也可以将动量梯度下降法「MGD」与自适应梯度算法「AdaGrad」融合在一起，得到自适应动量算法「Adaptive Momentum Estimation，简称 Adam」

在梯度下降过程中，Adam 算法用梯度模拟小球受力改变速度 \(v\)，同时会增大小球的质量而改变小球移动的难易程度，二者的作用最终影响移动距离。

code：

#include<iostream>
#include<cstdio>
#include<cmath>

using namespace std;

inline long double f(long double x)
{
    return 3*x*x*x*x-x*x*x+2*x*x-9*x+5*sqrt((x+3)*(x+3)+(5*x+6)*(5*x+6))-25;
}

inline long double numerical_diff(long double x)
{
    long double dx=1e-6;
    return (f(x+dx)-f(x-dx))/(dx*2);
}

inline long double Adam(long double x,long double eta,long double beta1,long double beta2)
{
    long double v=0,m=0,grad=numerical_diff(x),t,t1=beta1,t2=beta2,eps=1e-6;
    //int cnt=0;
    while(abs(grad)>eps)
    {
        t=eta*sqrt(1-t2)/(1-t1);
        v+=(1-beta1)*(grad-v);
        m+=(1-beta2)*(grad*grad-m);
        x-=t*v/(sqrt(m)+eps);
        t1*=beta1,t2*=beta2;
        grad=numerical_diff(x);
        //cnt++;
    }
    //cout<<cnt<<endl;
    return f(x);
}

signed main()
{
    cout<<Adam(0,0.5,0.6,0.9999)<<endl;
    return 0;
}

迭代次数：\(49\)。

综合来看，AdaGrad 的算法似乎更具优势，但毕竟各个算法特点不同，适用的函数也不同，用哪种算法还是应根据实际应用而定。

二维写法

code：

struct node{
    long double x,y;
};

inline long double f(long double x,long double y)
{
    return ...;
}

inline long double part_x(long double x,long double y)
{
    long double dx=1e-6;
    return (f(x+dx,y)-f(x-dx,y))/(dx*2);
}

inline long double part_y(long double x,long double y)
{
    long double dy=1e-6;
    return (f(x,y+dy)-f(x,y-dy))/(dy*2);
}

inline node gd(long double x,long double y,long double eta,long double beta)
{
    long double eps=1e-6;
    node res;
    int cnt=0,T=100;
    while(T--&&abs(part_x(x,y))+abs(part_y(x,y))>eps)
    {
        x-=eta*part_x(x,y);
        y-=eta*part_y(x,y);
        //cout<<x<<' '<<y<<' '<<f(x,y)<<' '<<++cnt<<endl;
    }
    res.x=x,res.y=y;
    return res;
}

inline node mgd(long double x,long double y,long double eta,long double beta)
{
    long double vx=0,vy=0,eps=1e-6;
    node res;
    int cnt=0,T=100;
    while(T--&&abs(part_x(x,y))+abs(part_y(x,y))>eps)
    {
        vx=beta*vx-eta*part_x(x,y);
        vy=beta*vy-eta*part_y(x,y);
        x+=vx,y+=vy;
        //cout<<x<<' '<<y<<' '<<f(x,y)<<' '<<++cnt<<endl;
    }
    res.x=x,res.y=y;
    return res;
}

inline node AdaGrad(long double x,long double y,long double eta,long double beta)
{
    long double vx=0,vy=0,hx=0,hy=0,gradx=part_x(x,y),grady=part_y(x,y),grad=abs(gradx)+abs(grady),eps=1e-6;
    node res;
    int cnt=0,T=100;
    while(T--&&grad>1e-5)
    {
        hx+=gradx*gradx;
        hy+=grady*grady;
        vx=beta*vx-eta*gradx/(sqrt(hx)+eps);
        vy=beta*vy-eta*grady/(sqrt(hy)+eps);
        x+=vx,y+=vy;
        gradx=part_x(x,y);
        grady=part_y(x,y);
        grad=abs(gradx)+abs(grady);
        //cout<<x<<' '<<y<<' '<<f(x,y)<<' '<<++cnt<<endl;
    }
    res.x=x,res.y=y;
    return res;
}

inline node Adam(long double x,long double y,long double eta,long double beta1,long double beta2)
{
    long double vx=0,vy=0,mx=0,my=0,gradx=part_x(x,y),grady=part_y(x,y),grad=abs(gradx)+abs(grady),t,t1=beta1,t2=beta2,eps=1e-7;
    node res;
    int cnt=0,T=100;
    while(T--&&grad>1e-5)
    {
        t=eta*sqrt(1-t2)/(1-t1);
        vx+=(1-beta1)*(gradx-vx);
        vy+=(1-beta1)*(grady-vy);
        mx+=(1-beta2)*(gradx*gradx-mx);
        my+=(1-beta2)*(grady*grady-my);
        x-=t*vx/(sqrt(mx)+eps);
        y-=t*vy/(sqrt(my)+eps);
        t1*=beta1,t2*=beta2;
        gradx=part_x(x,y);
        grady=part_y(x,y);
        grad=abs(gradx)+abs(grady);
        //cout<<x<<' '<<y<<' '<<f(x,y)<<' '<<++cnt<<endl;
    }
    res.x=x,res.y=y;
    return res;
}

例题

传送门：洛谷P2571 [SCOI2010]传送带

直接用 Adam 交了一发，AC～「可能运气较好……」

code：

#include<iostream>
#include<cstdio>
#include<cmath>
#include<climits>
#include<cstdlib>
#include<ctime>
#define Min(a,b,c,d,e) (min(a,min(b,min(c,min(d,e)))))

using namespace std;

long double P,Q,R,ans=INT_MAX;
long double c1,c2,c3,c4,c5,c6,c7,c8,c9;
long double x,y;

struct node{
    long double x,y;
}A,B,C,D,M,N,t;

inline long double dis(node a,node b)
{
    return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y));
}

inline long double f(long double x,long double y)
{
    return c1*x+c2*y+c3*sqrt((c4*x+c5*y+c6)*(c4*x+c5*y+c6)+(c7*x+c8*y+c9)*(c7*x+c8*y+c9));
}

inline long double part_x(long double x,long double y)
{
    long double dx=1e-6;
    return (f(x+dx,y)-f(x-dx,y))/(dx*2);
}

inline long double part_y(long double x,long double y)
{
    long double dy=1e-6;
    return (f(x,y+dy)-f(x,y-dy))/(dy*2);
}

inline node Adam(long double x,long double y,long double eta,long double beta1,long double beta2)
{
    long double vx=0,vy=0,mx=0,my=0,gradx=part_x(x,y),grady=part_y(x,y),grad=abs(gradx)+abs(grady),t,t1=beta1,t2=beta2,eps=1e-6;
    node res;
    int T=20;
    while(T--&&grad>1e-4)
    {
        t=eta*sqrt(1-t2)/(1-t1);
        vx+=(1-beta1)*(gradx-vx);
        vy+=(1-beta1)*(grady-vy);
        mx+=(1-beta2)*(gradx*gradx-mx);
        my+=(1-beta2)*(grady*grady-my);
        x-=t*vx/(sqrt(mx)+eps);
        y-=t*vy/(sqrt(my)+eps);
        t1*=beta1,t2*=beta2;
        gradx=part_x(x,y);
        grady=part_y(x,y);
        grad=abs(gradx)+abs(grady);
    }
    res.x=x,res.y=y;
    return res;
}

inline long double Rand()
{
    return abs(rand()*rand()*rand()%(10000000)/(long double)(10000000));
}

inline void gradient()
{
    for(int T=1;T<=200000;T++)
    {
        t=Adam(Rand(),Rand(),Rand(),Rand(),Rand());
        if(t.x>=0&&t.x<=1&&t.y>=0&&t.y<=1) ans=min(ans,f(t.x,t.y));
    }
}

signed main()
{
    srand(time(0));
    scanf("%Lf%Lf%Lf%Lf%Lf%Lf%Lf%Lf%Lf%Lf%Lf",&A.x,&A.y,&B.x,&B.y,&C.x,&C.y,&D.x,&D.y,&P,&Q,&R);
    c1=dis(A,B)/P,c2=dis(C,D)/Q,c3=1/R;
    c4=B.x-A.x,c5=D.x-C.x,c6=A.x-D.x,c7=B.y-A.y,c8=D.y-C.y,c9=A.y-D.y;
    gradient();
    ans=Min(f(0,0),f(0,1),f(1,0),f(1,1),ans);
    printf("%.2Lf\n",ans);
    return 0;
}

\[Thanks\quad for\quad reading. \]

————THE——END————

标签：梯度,long,eta,下降,dx,double,part
来源： https://www.cnblogs.com/lingyunvoid/p/16690821.html

本站声明： 1. iCode9 技术分享网（下文简称本站）提供的所有内容，仅供技术学习、探讨和分享；
2. 关于本站的所有留言、评论、转载及引用，纯属内容发起人的个人观点，与本站观点和立场无关；
3. 关于本站的所有言论和文字，纯属内容发起人的个人观点，与本站观点和立场无关；
4. 本站文章均是网友提供，不完全保证技术分享内容的完整性、准确性、时效性、风险性和版权归属；如您发现该文章侵犯了您的权益，可联系我们第一时间进行删除；
5. 本站为非盈利性的个人网站，所有内容不会用来进行牟利，也不会利用任何形式的广告来间接获益，纯粹是为了广大技术爱好者提供技术内容和技术思想的分享性交流网站。