最速下降法简单分析报告

作者: 陨落的小白 | 来源:发表于2020-11-03 19:40 被阅读0次

问题分析

该问题描述如下：给定1000组数据 $(x,y)$ ，找到一条直线 $f(x)=ax+b$ 使得该直线能够最好地表达 $y$ 与 $x$ 之间的线性关系。因此，我们需要找到相应的指标对于找到的直线进行评价。常用的模型是最小化其平方损失函数，即
$(a,b)=Argmin\ \sum_{i=1}^n(y_i-f(x_i))^2=Argmin\ \sum_{i=1}^n(ax_i+b-y_i)^2$
显然，当直线越能较好地表达 $y$ 与 $x$ 之间的线性关系， $\sum_{i=1}^n(y_i-f(x_i))^2$ 的值应该越小。倘若 $\sum_{i=1}^n(y_i-\hat{f}(x_i))^2=0$ ，则说明 $y$ 与 $x$ 完全线性相关，此时的 $\hat{f}(x)$ 就是 $y$ 与 $x$ 之间的线性表达式。

对于这类线性回归问题，往往存在着相应的闭式解，即最小二乘回归法得到的解。这里不对此闭式解进行描述。我们使用最速下降，或者说梯度下降的方法求解参数，尝试使得 $\sum_{i=1}^n(y_i-f(x_i))^2$ 达到最小。

梯度下降法

梯度下降法是求解优化问题的常用方法，其步骤相对比较简单。即给定一个待求参数的初始值，对于需要最小化的目标函数，求其关于参数的梯度，之后沿着负梯度方向对参数进行移动，即可以实现使得目标函数值下降的目的。其背后的数学原理主要用到了泰勒展开式
$f(x+\Delta x)=f(x)+\nabla f(x)·\Delta x+o(|\Delta x|)$
我们希望经过 $\Delta x$ 的一段移动， $f(x+\Delta x)$ 的函数值相对于 $f(x)$ 可以有所下降。在不考虑高阶无穷小的情况下，问题则转化成找到合适的 $\Delta x$ 使得 $\nabla f(x)·\Delta x$ 最小（同时也需要使得 $\nabla f(x)·\Delta x<0$ ）。对于这样一个内积形式，我们使用余弦公式
$\nabla f(x)·\Delta x=|\nabla f(x)|·|\Delta x|cos<\nabla f(x),\Delta x>$
将其展开，可以发现在 $|\nabla f(x)|$ 和 $|\Delta x|$ 固定的情况下，如果我们使得 $cos<\nabla f(x),\Delta x>$ 的值最小，即 $\Delta x$ 与 $\nabla f(x)$ 的方向相反时， $\nabla f(x)·\Delta x$ 可以取到最小。因此，参数移动过程中，目标函数下降最快的方向
$v=-\frac{\nabla f(x)}{|\nabla f(x)|}$
即负梯度方向。

而我们的 $\Delta x$ 可以表达为 $\Delta x=\alpha ’ v=-\alpha \nabla f(x)$ ，所以最终我们得到了梯度下降的迭代格式：
$x_{t+1}=x_{t}-\alpha \nabla f(x_{t})$

凸函数

梯度下降法固然是一个很有效的方法，通过沿着负梯度方向移动参数，使得目标函数以较快的速度下降。但是当梯度下降法移动到局部最优点时，考虑到此时局部最优点处的梯度为0，参数就没有动力继续进行移动了。如图，当参数由A点移动到B点时，就无法使用梯度下降法作进一步移动了，而全局最优点位于C处。

为了避免这种情况，我们可以提出一些有效的扰动机制，使得参数可以从B点跳动到D点，再从D点通过梯度下降到达C点。但是实际问题往往比较复杂，扰动机制也不像图片中那样明显。

而如果我们的目标函数相对于参数是一个凸函数，则不存在这种收敛到局部最优的情况了。对于定义在凸集上的凸函数，其局部最优点就是全局最优点。所以，如果我们求解的目标函数 $\sum_{i=1}^n(ax_i+b-y_i)^2$ 是一个定义在凸集上的凸函数，梯度下降法必然可以下降到其全局最优点处。

由于目标函数 $g(a,b)=\sum_{i=1}^n(ax_i+b-y_i)^2$ 是一个关于 $(a,b)$ 的二元函数，所以令
$A=g''_{aa}(a,b),B=g''_{ab}(a,b),C=g''_{bb}(a,b)$
，则对于 $\forall (a,b)\in R^2$ ，若有
$A>0,AC-B^2 \geq 0$
则 $g(a,b)$ 是一个凸函数（来源于微积分课本）。

经过计算，
$A=2\sum_{i=1}^nx_i^2,B=2\sum_{i=1}^nx_i,C=2n$
则
$AC-B^2=4n\sum_{i=1}^n(x_i^2-2x_i\bar{x}+\bar{x}^2)\geq 0,A>0$
所以，目标函数 $g(a,b)=\sum_{i=1}^n(ax_i+b-y_i)^2$ 是一个定义在凸集上的凸函数，梯度下降法可以得到其全局最优解。

代码实现

原理已经完全清晰，并且证明了梯度下降法的有效性，下一步就是代码实现。代码实现相对比较简单，首先给定参数的初始值 $(a_0,b_0)$ 和学习率 $\alpha$ 。然后计算出相应的梯度分量
$\frac{\partial g}{\partial a}=2\sum_{i=1}^nx_i(ax_i+b-y_i)$

$\frac{\partial g}{\partial b}=2\sum_{i=1}^n(ax_i+b-y_i)$

之后便可以在终止条件达到之前，不断进行迭代，即
$(a_{t+1},b_{t+1})=(a_t,b_t)-\alpha (\frac{\partial g}{\partial a},\frac{\partial g}{\partial b})$
至于终止准则，考虑到对于该目标函数，梯度下降可以收敛到全局最优点，所以我们设定终止准则为 $||(a_t-a_{t+1},b_t-b_{t+1})||_{2}<eps$ ，即迭代前后参数向量差的二范数小于某个较小的值时，迭代终止。

代码见附件。

结果及参数分析

结果分析

我们取初值为 $(a_0,b_0)=(0,0)$ ，学习率 $\alpha=10^{-5}$ ， $eps=10^{-8}$ ，最终结果收敛于 $(0.7009,1.4188)$ 。

image.png

我们使用matlab的curve Fitting工具箱中的线性拟合函数，直接对给定数据进行拟合，得到的结果是 $(0.7009,1.419)$ 。说明我们使用的梯度下降法基本收敛到最优解。

image.png

参数分析

由于我采用的是定步长梯度下降，所以算法中主要有初值、学习率以及 $eps$ 这三个参数。接下来我们分析一下参数的变化对于结果的影响。为了直观的展示参数的不同取值对于结果的影响，这里使用“均方误差的下降曲线”来呈现梯度下降法的下降效果。

初值的变化

考虑到最终的结果是 $(0.7009,1.4188)$ ，位于 $(0,0)$ 附近，所以可以以 $(0,0)$ 为参考点，选择不同的初值观察其迭代情况。这里 $\alpha=10^{-5},eps=10^{-8}$ .

image.png

首先选择 $(0,0)$ 周围八个点作为初始值，观察其迭代情况。

image.png

从最终的收敛性来看，八个不同方位的初值经过不同次数的梯度下降迭代后，都可以收敛到 $(0.7009,1.4188)$ 处。说明对于该问题给出的凸目标函数，从任何一个方位进行出发，选择恰当的学习率后，都可以很好地迭代到最优解的位置。

这八个初始值虽然最终都可以收敛到最优解，但其迭代次数具有一定的区别，其中有三个位置迭代了948次，有三个位置迭代了947次，还有两个位置分别迭代了734和740次。迭代次数的差异，我认为主要与数据集本身的特征有关，因为八个初始位置距离最优解的距离是差不太多的，而采取的学习率与终止准则又是相同的，所以应该是数据集本身的特征带来了迭代次数的差异。

为了进一步说明梯度下降在该问题中的强大的效果，我们采取几个较为“遥远”的位置作为初始的参数，分别为 $(-23749874,43874),(32479,78),(97427,2987),(4278947,-378)$ 并观察他们的下降情况。

image.png

实验后可知，这四个初始位置最终都顺利收敛到最优解，其迭代次数分别为 $1267,912,1114,1130$ 次。由此我们可以认为，对于不加约束的二元凸优化问题，初始位置的选取对于是否收敛并没有太大的影响。

学习率的变化

学习率 $\alpha$ 是一个非常重要的参数，决定了梯度下降能否收敛到最优解以及其收敛速度。由于该问题数据集较大，一般点处的梯度分量较大，所以为了使得解的每一步变化在一个合适的范围内，学习率应该取一个较小的值，最终我取了 $10^{-5}$ 。我们取 $\alpha\in \{10^{-20},10^{-10},10^{-8},10^{-7},10^{-6},10^{-5},10^{-4},10^{-3}\}$ ，初值取 $(0,0)$ ，终止准则取 $eps=10^{-8}$ 。同时，考虑到发散的情况，我们设置一个最大迭代次数为2千万次。

由于最终的结果较难用均方误差下降曲线表示，我们使用表格来表示这一结果

$\alpha$ 取值	迭代次数	是否收敛到最优解
$10^{-20}$	1	否
$10^{-10}$	16725705	接近最优解
$10^{-8}$	397517	是
$10^{-7}$	51261	是
$10^{-6}$	6273	是
$10^{-5}$	737	是
$10^{-4}$	20000000	否
$10^{-3}$	20000000	否

由表格可以发现，学习率 $\alpha$ 对于能够收敛以及收敛速度都具有非常大的影响。虽然还未考虑终止准则，但是理论分析来看，终止准则主要影响收敛时的解的精度，初始位置主要影响迭代次数，所以学习率对于梯度下降法的影响是最大的。

对于该问题，当 $\alpha=10^{-20}$ ，由于一次迭代过程中的变动太小，很容易就小于 $eps$ ，所以仅迭代一次就收敛了。倘若改变 $eps$ ，也有可能在很久的迭代之后收敛。当 $\alpha =10^{-10}$ 次方时，同样的，由于每次迭代过程中的变化幅度太小，所以需要迭代一千多万次才能收敛到一个近似最优解的地方。如果我们将 $eps$ 设置的再小一些，也是有可能收敛到最优解的。当 $10^{-10}\leq \alpha \leq 10^{-5}$ 时，我们可以发现，经过足够多的迭代次数，最终总能收敛到最优解的位置。但是 $\alpha$ 越大时，每一步的变化幅度越大，也就能以更快地速度收敛到最优解，相应的迭代次数也就越少。当 $\alpha\geq 10^{-4}$ 时，由于这时候一步的跨度太大，以至于其偏离最优解越来越远，最终导致结果发散无法收敛。

image.png

所以， $\alpha$ 的选取必须要合适。太小会导致一次就收敛，较小会导致收敛较慢，较大又会导致结果发散。我认为 $\alpha$ 的选取主要与梯度分量的大小有关。我们可以观察 $\alpha$ 取不同的值的时候，梯度分量，例如 $\frac{\partial g}{\partial a}$ 的变化。由图可知，一个恰当的 $\alpha$ 应该使得梯度分量在迭代的过程中慢慢趋于0，也就是 $\alpha$ 需要与梯度的大小相匹配。

image.png

如果 $\alpha$ 与梯度不匹配的话，例如 $\alpha=10^{-4}$ 时，我们发现即使迭代次数在十次以内，其梯度分量就开始呈震荡趋势，且其绝对值越来越大。因此，一个合适的 $\alpha$ 必须与其梯度的规模相匹配，才能以较好的步伐收敛到最优解。

image.png

$eps$ 的变化

$eps$	迭代次数	收敛结果
$10^{-2}$	53	$(0.70134,0.93267)$
$10^{-4}$	281	$(0.70093,1.4193)$
$10^{-6}$	509	$(0.70093,1.4187)$
$10^{-8}$	737	$(0.70093,1.4188)$
$10^{-10}$	965	$(0.70093,1.4188)$

显然，一般情况下， $eps$ 越小，迭代次数越多，相应的收敛结果也越精确。

总结

首先，我们可以作出一张 $g(a,b)=\sum_{i=1}^n(ax_i+b-y_i)^2$ 的函数图像。

image.png

可以发现，其确实明显存在一个下凸的形状，这也直观解释了为什么初始位置的选择不太影响其收敛的结果。

对于定步长的梯度下降法，我们发现它可以很好地解决一元函数线性拟合的问题，且往往可以快速收敛到最优解。在所有的参数中，初始位置的重要性不大，因为该目标函数是一个典型的凸函数，且其定义域是 $R^n$ ，也不存在迭代过程中离开定义域的可能。终止准则的设定对于结果的精度有较大的影响，但是对于结果的收敛性影响不太大，一般而言，精度越大，需要的迭代次数越多。学习率 $\alpha$ 对于是否收敛以及收敛的速度具有十分重大的影响。一般而言， $\alpha$ 太小或太大时，往往导致梯度下降法迭代效率极低，亦或者无法收敛。 $\alpha$ 的选取与梯度有很大的关系，但对于更加复杂的函数，其相关因素应该会更多。所以 $\alpha$ 的值需要多多尝试。

对于变步长梯度下降法，我也进行了实验，例如每次取 $\alpha=0.99*\alpha$ ，但是实际效果甚至没有定步长好，所以在此不做过多阐述。事实上对于这一题，定步长已经可以在很快的速度内得到最优解了。

以上就是本次分析报告的全部内容。

下面附上代码（写得有点儿烂，请勿介意）

%主函数
clear;clc;
load a1data.mat
a=[0;0];% 初始值，a(1)=a;   a(2)=b; y=ax+b
alpha=1e-5; %学习率 alpha
eps=1e-8; %终止准则 eps
[q,s,num]=MSE(a,x,y,alpha,eps); %q是最终的解，s的均方误差的变化向量，num是迭代次数
figure(1);
plot(1:length(s),s);
disp("迭代了"+num+"次"+"   收敛结果是"+q)
title("均方误差下降曲线")
xlabel("迭代次数")
ylabel("均方误差")
y1=q(1)*x+q(2);
figure(2);
plot(x,y,'or',x,y1,'b');

%作出SSE关于参数的三维图像
figure(3)
[n,m]=meshgrid(-10:0.01:10);
v=n.^2*sum(x.^2)+2*n.*m*sum(x)+1000*m.^2+sum(y.^2)-2*sum(x.*y)*n-2*m*sum(y);
mesh(n,m,v);

%迭代函数
function [q,s,num,daa]=MSE(a,x,y,alpha,eps)
i=1; %迭代序号
s=1:10;  %MSE向量
num=0; %迭代次数
q=[inf;inf]; %返回结果的初值
daa=1:10; % 梯度对a的分量向量
while 1
   % disp(i);
    da=2*(a(1)*sum(x.^2)+a(2)*sum(x)-sum(x.*y)); %梯度对a的分量
    db=2*(a(1)*sum(x)+length(x)*a(2)-sum(y));%梯度对b的分量
    a1=a-alpha*[da;db]; %迭代格式
    %disp(a1);
    num=i; 
    if norm(a1-a)<eps  %终止准则
        q=a1;
        break;
    end
    if i>20000
        disp("迭代了20000次也没有收敛，停下来改一下最大迭代次数试试")
        break;
    end
    a=a1;
    s(i)=sum((a1(1)*x+a1(2)-y).^2)/length(x); %计算MSE
    daa(i)=da;
    i=i+1;
end
end

原始数据

链接：https://pan.baidu.com/s/1ZdnBWuG8ToCYD4d1VbJc8g

提取码：cybr

最速下降法简单分析报告

问题分析

梯度下降法

凸函数

代码实现

结果及参数分析

结果分析

参数分析

初值的变化

学习率的变化

$eps$ 的变化

总结

相关文章

网友评论

延伸阅读

深度阅读

栏目导航

热点阅读