利用基本梯度下降法和牛顿法对样本进行判别

机器学习实验五，详情请参考《模式分类》第二版第五章课后上机练习5.4节

实验环境：

Matlab 2016a

基本梯度下降法和牛顿法：

我们在寻找能将两类类别分开的权向量时采用的方法是：定义一个准则函数J（a），当a是解向量时，J（a）最小。这样就将问题简化为一个标量函数的极小化问题——通常可以用梯度下降法来解决。梯度下降法的原理非常简单，首先从一个任意选择的权向量a（1）开始，计算其梯度向量，下一个值a（2）由自a（1）向下降最陡的方向移一段距离得到，即沿梯度的负方向。通常a（k+1）由下等式确定：

其中 是正的比例因子，或者说是用于设定学习步长的“学习率”。我们希望这样得到的一个权向量序列：最终收敛到一个使J（a）极小化的解上。具体步骤是初始化权向量a、阈值和，不断迭代更新a，直到<，使得准则函数达到一个极小值，a收敛。

而牛顿法权向量的更新公式为：

其中，H为准则函数的赫森矩阵。因为牛顿法使用了准则函数的二次偏导，因此牛顿算法比梯度下降算法每一步都给出了更好的步长，也就更快收敛。

实验内容：

给定以下两个类别w1、w3数据，应用基本梯度法和牛顿法：

用这两种算法对二维数据给出w1、w3的判别。对梯度下降法取=0.1。画出准则函数-迭代次数的曲线。

实验分析：

实验的关键就是准则函数的选取，首先我们知道当错分点减少时，准则函数的值应该减小，并且准则函数至少是关于权向量二次的，因为牛顿算法要二次求导。最终使用书中的如下函数作为准则函数：

这里的是满足的样本集。如果为空的话，我们定义J（a）为0。这样J（a）就不是负的。J（a）的梯度为：

J（a）的二次偏导为：

将其代入各自的更新公式即可。

实验步骤：

确定了准则函数后，剩下的就是编程的事。首先，我们知道这里的权向量和样本向量都是增广的，而且样本集需要先规范化：

load sample_ex5.mat;%规范化矩阵Y[rows,cols]=size(w1);one=ones(rows,1);w1=[one,w1];w3=[one,w3];Y=[w1;(-1)*w3];%初始化权向量a，阈值u，初始学习率sa=[1;1;1];a1=a;u=0.01;s=0.1;

现在开始迭代，首先需要选出的那些样本集，代入式子（1）叠加求和，因为我们还需要计算准则函数的值，因此用value1数组存储准则函数的值，用delta1变量存储每次更新a的增量：

for i=1:2*rows        if(a1'*Y(i,:)'<=0)        value1(1,count1)=value1(1,count1)+(a1'*Y(i,:)')^2/sum(Y(i,:).*Y(i,:));        delta1=delta1+s*(a1'*Y(i,:)'/sum(Y(i,:).*Y(i,:))*Y(i,:)');        end    endvalue1(1,count1)=value1(1,count1)/2;  %计算准则函数的值    a1=a1-delta1;    %更新a

然后我们判断是否需要停止迭代：

if(norm(delta1)<u)        for i=1:2*rows            if(a1'*Y(i,:)'<=0)                value1(1,count1+1)=value1(1,count1+1)+(a1'*Y(i,:)')^2/sum(Y(i,:).*Y(i,:));%收敛时计算最后一次准则函数的值            end        end        value1(1,count1+1)=value1(1,count1+1)/2;        break;    end

如果不满足条件，继续迭代直到收敛，此时的a即为最终的权向量。

对于牛顿法，基本思路差不多，先求出，再求出，用的逆乘，算出增量更新a。
同样，准则函数的值存于value2数组中：

 for i=1:2*rows        if(a2'*Y(i,:)'<=0)            value2(1,count2)=value2(1,count2)+(a2'*Y(i,:)')^2/sum(Y(i,:).*Y(i,:));            delta_J=delta_J+(a2'*Y(i,:)'/sum(Y(i,:).*Y(i,:)))*Y(i,:)';            H=H+(Y(i,:)'*Y(i,:))/sum(Y(i,:).*Y(i,:));        end    end

更新a：

delta2=H\delta_J;  %inv(H)*delta_J      a2=a2-delta2;  %更新a

更新完后判断是否需要继续迭代，如果收敛就跳出，并求出最后一次准则函数的值：

if(norm(delta2)<u)        for i=1:2*rows            if(a2'*Y(i,:)'<=0)                value2(1,count2+1)=value2(1,count2+1)+(a2'*Y(i,:)')^2/sum(Y(i,:).*Y(i,:));%收敛时计算最后一次准则函数的值            end        end        value2(1,count2+1)=value2(1,count2+1)/2;        break;    end

实验结果：

最后画出准则函数-迭代次数的曲线：

从图中明显看出，基本梯度下降法迭代20次才收敛，而牛顿法迭代2次就收敛了，并且我们注意到，最后准则函数的值并不为0，这是因为原始样本不是线性可分的，见下图，不存在一个权向量使得所有样本都分类正确。

虽然牛顿算法在每一步都给出更好的步长，但是每次递归都要计算赫森矩阵，时间复杂度为O（d^3），运算量过大。实际上，将学习率设置为一个较小的常数，虽然比每一步使用最优的将需要更多的步骤来校正，但通常总的时间花销更少。

现在我们来求收敛时间-学习率曲线。收敛时间可用迭代次数来代替，即求迭代次数-学习率曲线，我们可以大概的猜想一下，在一定范围内，学习率增大，迭代次数应该减小，当学习率过大时，可能会发生过冲甚至发散，导致不收敛。我们在第一问的基础上，设置学习率从0.01增大到1，每次增加0.01（0.01：0.01：1），求出每次的迭代次数：

for s=0.01:0.01:1    count1=0;    while(1)        count1=count1+1;        delta1=0;        for i=1:2*rows            if(a1'*Y(i,:)'<=0)                delta1=delta1+s*(a1'*Y(i,:)'/sum(Y(i,:).*Y(i,:))*(Y(i,:)'));            end        end        a1=a1-delta1;        if(norm(delta1)<u)            break;        end    end     count(1,uint8(s*100))=count1;%     if(count1>count(1,1))%         count(1,uint8(s*100))=inf;%     endEnd

由于当学习率增大到一定程度，a不收敛即迭代次数无穷大，判断当迭代次数过大时直接设置为inf：

由于matlab中不画出值为inf的点，因此我们可以看到后面的点都没有，因此不收敛的最小学习率为0.63。

PS: 个人实验结果可能有误,望指正。

总结：

基本梯度下降法和牛顿法都能求得最终的权向量，使得准则函数取得一个极小值。基本梯度下降法利用一阶导，而牛顿法利用二阶导，牛顿法每次都是使用基本梯度的最优情况，因此牛顿法下降步长大，较快收敛。但是当赫森矩阵为奇异矩阵就不能使用牛顿法了。而且，即使赫森矩阵非奇异，每次递归计算H逆矩阵所需的O（d^3）时间可轻易地将牛顿算法带来的好处给抵消了。实际上，设置学习率为一个较小的值，虽然需要更多步骤来达到收敛，但通常总的时间花销更小。
在使用这些方法的时候还有一个问题：如何选择学习率，如果学习率太小，收敛非常慢，如果太大，可能会发生过冲甚至发散。

附（矩阵函数对矩阵求导公式）：

A, B, C 是不依赖于 X 的矩阵，a,b 是不依赖于x 的向量。