基于matlab的梯度下降法实现线性回归

基于matlab的梯度下降法实现线性回归

1 绪论

1.1线性回归的定义

1.2单变量线性回归

1.3多变量线性回归

2 梯度下降

2.1 cost function

2.2 梯度下降：解决线性回归的方法之一

2.3 feature scaling：加快梯度下降执行速度的方法

2.4 梯度下降的实现

3 梯度下降的扩展和比较

摘要:本文首先说明线性回归的定义，以及在统计机器学习中的意义和现实生活中的影响，引出梯度下降法的需求和好处。然后分析了梯度下降的基本原理，给出基于matlab的具体实现方法。总结梯度下降算法的特点，并分析相关参数对算法运行结果的影响。最后，扩展并比较了其他几种算法(批梯度下降算法、增量梯度下降算法、最小二乘法)以及相关实现。

关键词: 线性回归；梯度下降；机器学习

1 绪论

1.1线性回归的定义

线性回归，是利用数理统计中回归分析，来确定两种或两种以上变量间相互依赖的定量关系的一种统计分析方法，运用十分广泛。其表达形式为y = w’x+e，e为误差服从均值为0的正态分布。

具体来说，如果输入 x是列向量，目标 y 是连续值（实数或连续整数），预测函数 f(x)的输出也是连续值。这种机器学习问题是回归问题。如果我们定义 f(x)是线性函数，f(x) = wTx + b, 这就是线性回归问题（Linear Regression）。

1.2单变量线性回归

单变量线性回归就是从一个输入值预测一个输出值。输入/输出的对应关系就是一个线性函数。

1.3多变量线性回归

多变量线性回归就是说输出值受多个输入变量x影响，每个变量的影响大小用w(weight)来刻画，w与x呈线性关系，然后把多个wx线性组合(w1x1,w2x2,,,)相加，再和输出y的对应关系就是一个多变量线性回归问题。注意：多变量线性回归不是指y与X是线性的，而是指各自的w与x呈线性关系。所以最终的结果，在几何上甚至可能呈现出曲线，而不是必须是直线。

2 梯度下降

2.1 cost function

线性回归属于监督学习，因此方法和监督学习应该是一样的，先给定一个训练集，根据这个训练集学习出一个线性函数，然后测试这个函数训练的好不好（即此函数是否足够拟合训练集数据），挑选出最好的函数（cost function最小）即可。

1)下面给出单变量线性回归的模型：

从上面“方法”中，我们肯定有一个疑问，怎么样能够看出线性函数拟合的好不好呢？

我们需要使用到Cost Function（代价函数），代价函数越小，说明线性回归地越好（和训练集拟合地越好），当然最小就是0，即完全拟合。

2)举个实际的例子：

我们想要根据房子的大小，预测房子的价格，给定如下数据集：



根据以上的数据集画在图上，如下图所示：



我们需要根据这些点拟合出一条直线，使得cost Function最小。虽然我们现在还不知道Cost Function内部到底是什么样的，但是我们的目标是：给定输入向量x，输出向量y，theta向量，输出Cost值。

3)Cost Function的用途：对假设的函数进行评价，cost function越小的函数，说明拟合训练数据拟合的越好。

下图详细说明了当cost function为黑盒的时候，cost function 的作用。



4)但是我们肯定想知道cost Function的内部构造是什么？因此我们下面给出公式：



其中：


表示向量x中的第i个元素；

表示向量y中的第i个元素；


表示已知的假设函数； m 为训练集的数量。

比如给定数据集(1,1)、(2,2)、(3,3)

则x = [1;2;3]，y = [1;2;3]（此处的语法为Octave语言的语法，表示3*1的矩阵）

如果我们预测theta0 = 0，theta1 = 1，则h(x) = x，则cost function：

J(0,1) = 1/(2*3) * [(h(1)-1)^2+(h(2)-2)^2+(h(3)-3)^2] = 0；

如果我们预测theta0 = 0，theta1 = 0.5，则h(x) = 0.5x，则cost function：

J(0,0.5) = 1/(2*3) * [(h(1)-1)^2+(h(2)-2)^2+(h(3)-3)^2] = 0.58；

如果theta0 一直为 0， 则theta1与J的函数为：



如果有theta0与theta1都不固定(即都为变量)，则theta0、theta1、J 的函数为：



当然我们也能够用二维的图来表示，即等高线图。



注意：如果是线性回归，则costfunctionJ与theta0、theta1的函数一定是碗状的，即只有一个最小点。

2.2 梯度下降：解决线性回归的方法之一

但是又一个问题引出了，虽然给定一个函数，我们能够根据cost function知道这个函数拟合的好不好，但是毕竟函数有这么多，总不可能一个一个试吧？因此我们引出了梯度下降：能够找出cost function函数的最小值；梯度下降原理：将函数比作一座山，我们站在某个山坡上，往四周看，从哪个方向向下走一小步，能够下降的最快。当然解决问题的方法有很多，梯度下降只是其中一个，还有一种方法叫Normal Equation。

方法：

(1)先确定向下一步的步伐大小，我们称为Learning rate；

(2)任意给定一个初始值：



；

(3)确定一个向下的方向，并向下走预先规定的步伐，并更新



；

(4)当下降的高度小于某个定义的值，则停止下降。

算法：



特点：

(1)初始点不同，获得的最小值也不同，因此梯度下降求得的只是局部最小值；

(2)越接近最小值时，下降速度越慢。

问题a：如果初始值就在local minimum的位置，则会如何变化？

答：因为已经在local minimum位置，所以derivative肯定是0，因此不会变化。

如果取到一个正确的值，则cost function应该越来越小。

问题b：怎么取值？

答：随时观察值，如果cost function变小了，则ok，反之，则再取一个更小的值。

下图就详细的说明了梯度下降的过程：



从上面的图可以看出：初始点不同，获得的最小值也不同，因此梯度下降求得的只是局部最小值。

注意：下降的步伐大小非常重要，因为如果太小，则找到函数最小值的速度就很慢，如果太大，则可能会出现overshoot the minimum的现象；下图就是overshoot minimum现象：



如果Learning rate取值后发现J function 增长了，则需要减小Learning rate的值。

因此我们能够对cost function运用梯度下降，即将梯度下降和线性回归进行整合，如下图所示：



梯度下降是通过不停的迭代，而我们比较关注迭代的次数，因为这关系到梯度下降的执行速度，为了减少迭代次数，因此引入了Feature Scaling。

2.3 feature scaling：加快梯度下降执行速度的方法

思想：将各个feature的值标准化，使得取值范围大致都在-1<=x<=1之间。常用的方法是Mean Normalization，即



或者：[X-mean(X)]/std(X)。

举个实际的例子，有两个Feature：

(1)size，取值范围0~2000；

(2)#bedroom，取值范围0~5；

则通过feature scaling后，





2.4 梯度下降的实现

1)代码：

%绘制出来最原始的方程的图像，方程为f=x^2+y^2；定义域分别设置成为[-10, 10]，然后随便的设置一个初始点以及迭代的值和步长，来查看使用梯度下降算法的迭代过程%
[x, y] = meshgrid(-10:0.5:10, -10:0.5:10);
z = x.^2 + y.^2;
mesh(x, y, z);
%分别求x和y的偏导数%
%dx = 2*x;
%dy = 2*y;
%设置x，y的初始值，并且每步骤更新的值都加入到这个记录中来，第一个是初始值，为了迭代方便设置下第二个元素的值
record_values = [-9;-9];
%设置每次更新的步长%
step = 0.1;
%设置迭代次数%
count = 20;
%开始迭代%
for i=1:count
current_x = record_values(1,i);
current_y = record_values(2,i);
temp = [current_x - step*2*current_x, current_y - step*2*current_y];
%把temp附加到record_values后面%
record_values(:, i+1) = temp;
end
%计算出来点绘制出来这些点按照顺序连线%
hold on;
x_values = record_values(1,:);
y_values = record_values(2,:);
z_values = x_values.*x_values + y_values.*y_values;
plot3(x_values, y_values,z_values, 'k--o');
%我们知道当x,y等于0的时候是最小点，标示出来%
plot3(0, 0, 0, 'rp');

2)截图：





3 梯度下降的扩展和比较

1)代码：

clc

close all

n=1;%特征向量个数
m=20;%样本个数
x=ones(n+1,m);
x(2,:)=linspace(0,10,m);%采样点
y=13*ones(1,m)-3.36*x(2,:)+4*rand(1,m)-2;%样本值
scatter(x(2,:),y)
hold on
theta=[0,1]';%设置参数
xx=ones(n+1,m);
xx(2,:)=linspace(0,10,m);
yy=theta'*xx;%绘制算法运行生成的曲线
scatter(x(2,:),y)
title('批梯度下降算法结果')
hold on
plot(xx(2,:),yy,'k')
hold off
row=length(xx);
times=0;
%%
tic;
while(1)
times=times+1;
step=1/(times*row);
temp0=theta(1);
temp1=theta(2);
H=theta'*x;
theta=theta-step*((H-y)*x')';%批梯度下降算法
distant=abs(theta(1)-temp0)+abs(theta(2)-temp1);
if(distant<0.001)%判断是否已经收敛
break;
end
%     if(times>500)
%         break;
%     end

scatter(x(2,:),y)
hold on
yy=theta'*xx;
plot(xx(2,:),yy)
hold off
end
toc;
title('批梯度下降算法结果')

display('批梯度下降')
theta(1)
theta(2)
times
%%
theta=[0,1]';%设置参数
times=0;
figure(2)
tic;
while(1)
times=times+1;
step=1/(times*row);
temp0=theta(1);
temp1=theta(2);
H=theta'*x;
for j=1:m
theta=theta-step*((H(j)-y(j))*x(:,j));%增量梯度下降算法
end
distant=abs(theta(1)-temp0)+abs(theta(2)-temp1);
if(distant<0.001)%判断是否已经收敛
break;
end
%     if(times>500)
%         break;
%     end

scatter(x(2,:),y)
hold on
yy=theta'*xx;
plot(xx(2,:),yy)
hold off
end
toc;
title('增量梯度下降算法结果')
display('增量梯度下降')
theta(1)
theta(2)
times

%%
%用广义逆的最小二乘
%theta=(X' *X)^(-1)*y;
tic;
theta=inv(x*x')*x*y';
toc;
figure(3)
scatter(x(2,:),y)
hold on
yy=theta'*xx;
plot(xx(2,:),yy)
hold off
title('最小二乘解析式结果')
display('最小二乘解析式梯度下降')
theta(1)
theta(2)

2)截图：









附录本篇论文下载链接