单变量线性回归（机器学习笔记二）

前言

上一章我们学习了机器学习的应用领域、定义以及算法的分类，这一章来学习最简单的一个机器学习算法——单变量线性回归算法，它是一种监督学习的算法，而且输出值是连续变化的值，因此是一种“回归”算法。

一、符号定义

我们通过吴老师课堂上给出的练习题中的例子来学习单变量线性回归算法。
　　在进行监督学习算法的时候，我们首先要有一个训练集，这里的训练集就是人们统计的50名不同年龄的儿童和他们的身高，如下图所示。
　　　　
这里我们开发机器学习算法的目的就是找出身高与年龄的大致对应关系，从而可以预测特定年龄小孩的大致身高。
　　因此，这个训练集中，我们的输入值是年龄，输出值是身高。通常来讲，我们称输入变量为特征（feature），用$x$表示，输出值定义为$y$，这里的$x$就是小孩的年龄了，因为只有这么一个输入变量，因此这个问题是一个单变量的问题，如果有多个输入变量的话可以用$x_1,x_2,x_3\cdots$表示。我们用$m$表示训练样本的数量，在本例中，$m$就是小孩的数量，也就是$m=50$。我们用符号（$x^{(i)}$，$y^{(i)}$）表示第$i$个小孩的年龄和身高，$i=1,2,3\cdots 50$。
　　接下来，为了找出x和y之间的关系，我们把x和y画成散点图，如下所示
　　　　　　　
从数据点的分布可以看出，身高和年龄之间大致呈线性关系，因此假设身高与年龄之间的函数关系是一个线性函数。这个假设函数（hypothesis）我们将其命名为$h$，即$y=h(x)$，机器学习算法的目的就是找出$h(x)$的具体表达式。在本例中，我们认为假设函数是一条直线，于是有$h(x)={\theta }_0+{\theta }_{1}x$，其中${\theta }_1$表示直线的斜率，${\theta }_0$表示直线在$x$轴上的截距。
　　这样，我们的研究目标就转化为寻找最优的${\theta }_0$和${\theta }_1$。等学习了后面的内容后，我们就会发现，监督学习算法的最终目的通常就是寻找最优的参数。
　　接下来，为了让我们的表达方法更具有普适性和简洁性，我们把变量都用向量的方式来表示。我们定义函数$h(x)$的参数为向量$\theta =$$\left[\begin{array}{c}{\theta }_0\\ {\theta }_1\end{array}\right]$，同时，我们将$x$也写成向量的形式，为了让$x$和$\theta$对应起来，我们给它增加一个常数1，变成：$\mathbf{x}=$$\left[\begin{array}{c}1\\ x\end{array}\right]$，由于只有一个输出量，因此$y$的向量形式为$\mathbf{y}=$$\left[\begin{array}{c}y\end{array}\right]$。这样，我们样本集的输入值就可以写成一个$2\times m$的矩阵：$\left[\begin{array}{ccc}1& 1& \cdots \\ x^{(1)}& x^{(2)}& \cdots \end{array}\right]$，输出值就可以写成一个$1\times m$的矩阵$\left[\begin{array}{ccc}{y}^{(1)}& y^{(2)}& \cdots \end{array}\right]$。这样一来我们的假设函数就可以用向量表示为$h={\theta }^{\mathbf{T}}\times \mathbf{x}$，这里的$\times$是矩阵的乘积，我们也可以省略这个乘号，写成$h={\theta }^{\mathbf{T}}\mathbf{x}$。这种表达方式对后面的多变量回归乃至逻辑回归都是适用的，同时也可以方便地转化为Matlab代码。因此，以后我们就用不带下标的字母$\theta$表示参数向量，用${\theta }_0,{\theta }_1\cdots$表示向量中的元素；用$\mathbf{x}$表示输入向量，$x_j$表示第$j$个输入变量，$x_j^{(i)}$表示第$i$个样本的第$j$个输入变量；同样地，$\mathbf{y}$表示输出向量，$y_j$表示第$j$个输出量，$y_j^{(i)}$表示第$i$个样本的第$j$个输出量
　　至此，我们就完成了单变量线性回归的准备工作了，下面开始考虑如何寻找最优的参数$\theta$。

二、代价函数

为了寻找最优的参数$\theta$，我们首先要确定一个标准，即什么样的参数是“优”的。直观地理解，对于每一个样本，我们的函数值$h(x^{(i)})$与实际的输出值$y^{(i)}$符合得越好，我们认为对应的参数越“好”。我们可以使用方差$[h(x^{(i)})-y^{(i)}{]}^2$来定量地描述这两个值的符合程度，方差越大，说明这两个值的差异越大，方差为零时这两个值完全相同。
　　综合考虑所有的样本，我们可以用所有样本的方差之和$J=\sum _{i=1}^m[h(x^{(i)})-y^{(i)}{]}^2$来表征函数$h$与样本集整体的符合程度，$J$越小，符合程度越高，我们称函数$J$为函数$h$的代价函数。
　　从函数$J$的定义可以看出它的大小只与参数$\theta$有关，即$J=J(\theta )$。这样一来，我们的目标就再一次转化为寻找使代价函数值$J$最小的参数$\theta$。

三、梯度下降算法

那么，如何寻找$J$最小时的参数$\theta$呢？这就要用到梯度下降算法了。
　　梯度下降算法的思想很简单，$J({\theta }_0,{\theta }_1)$实际上是一个关于${\theta }_0和{\theta }_1$的二维曲面，每一对$({\theta }_0,{\theta }_1)$都对应于曲面上的一个点，如下图所示：
　　　　
我们先任意初始化一个$\theta$值，然后找出该点处$J$下降最快的方向，将$\theta$向这个方向移动一小步，然后再重新寻找$J$下降最快的方向，再向新的方向移动一小步，如此循环，直到$J$达到最小值，这个过程中$\theta$的移动路径如上图中黑色四芒星所示。
　　而函数$J$下降最快的方向其实就是函数$J({\theta }_0,{\theta }_1)$的梯度的反方向。根据梯度的定义，函数$J$的梯度为$\nabla J=\frac{\partial J}{\partial {\theta }_0}\overrightarrow{{\theta }_0}+\frac{\partial J}{\partial {\theta }_1}\overrightarrow{{\theta }_1}$，其中，$\overrightarrow{{\theta }_0}$和$\overrightarrow{{\theta }_1}$分别代表沿${\theta }_0$和${\theta }_1$正方向的单位向量，为了使$\theta$往$J$下降最快的方向移动，应该使${\theta }_0$和${\theta }_1$分别减去$\alpha \ast \frac{\partial J}{\partial {\theta }_0}$和$\alpha \ast \frac{\partial J}{\partial {\theta }_1}$，其中$\alpha$(>0)代表移动的步长。于是，梯度下降算法的计算过程如下所示：
　　Step1、初始化$\theta$：例如令${\theta }_0={\theta }_1=0$;
　　Step2、计算当前$\theta$值对应的$J$的梯度值：$\frac{\partial J}{\partial {\theta }_0}$,$\frac{\partial J}{\partial {\theta }_1}$；
　　Step3、更新$\theta$：${\theta }_0:={\theta }_0-\alpha \ast \frac{\partial J}{\partial {\theta }_0}$；${\theta }_1:={\theta }_1-\alpha \ast \frac{\partial J}{\partial {\theta }_1}$;
　　Step4、判断是否结束循环，若为否则回到Step2。
　　从以上过程可以看出，每一步的移动步长是由该点处$J$的梯度大小和$\alpha$共同决定的。在$\alpha$大小不变的情况下，当$J$离最小值较远时，对应的梯度值较大，$\theta$移动的较快；$J$越接近最小值，$\theta$移动步长越小，最终逐渐收敛到$J$最小处。
　　而第二步中偏导数$\frac{\partial J}{\partial {\theta }_0}$,$\frac{\partial J}{\partial {\theta }_1}$的值可以将$h(x)={\theta }_0+{\theta }_{1}x$代入$J=\sum _{i=1}^m[h(x^{(i)})-y^{(i)}{]}^2$中计算得出。根据偏导数的运算规则，可以得出：
$\frac{\partial J}{\partial {\theta }_0}=2\sum _{i=1}^m[h(x^{(i)})-y^{(i)}]$；$\frac{\partial J}{\partial {\theta }_1}=2\sum _{i=1}^m{x}^{(i)}[h(x^{(i)})-y^{(i)}]$，实际执行时，会在前面乘上一个系数$\frac{1}{2m}$,并不会影响最终的优化结果。

四、Matlab算法实现

接下来，用Matlab来实现以上过程，代码如下：

%导入训练集；
x=load('ex2x.dat');
y=load('ex2y.dat');

%画出训练集的散点图
figure,
subplot(121),
plot(x,y,'o');
ylabel('Height in meters');
xlabel('Age in years');
hold on;

%step1、初始化相关参数，注意：这里我们将x和y写成了m*2和m*1的向量形
%式，这样一来h=x*theta，与文中所述有所不同。
m=length(y);
x=[ones(m,1),x];
n=2;
theta=zeros(n,1);
alpha=0.07;

%开始单变量梯度下降循环，并实时显示J的变化
subplot(122);
for step=1:800
    %先计算h-y的值，h-y是一个m×1的向量
    h_y=x*theta-y;
    %step2、3 计算梯度并更新θ的值，通过向量乘法实现样本的求和。
    theta(1)=theta(1)-alpha/m*(h_y'*x(:,1));
    theta(2)=theta(2)-alpha/m*(h_y'*x(:,2)); 
    %计算并显示代价函数值
    J=sum(h_y.^2)/m/2;
    if mod(step,10)==9
        plot(step,log10(J),'o');
        axis([0,800,-inf,inf]);
        ylabel('log(J)');
        xlabel('Step');
        hold on;
        drawnow;
    end   
end

%显示最终的优化结果 
subplot(121),
tempx=linspace(2,8,100);
tempy=theta(1)+tempx*theta(2);
plot(tempx,tempy);

%准备数据，画出代价函数曲面
Jval=zeros(100,100);
theta0=linspace(-3,3,100);
theta1=linspace(-1,1,100);
for xs=1:length(theta0)
    for ys=1:length(theta1)
        theta=[theta0(xs);theta1(ys)];
        h_y=x*theta-y;
        Jval(xs,ys)=sum(h_y.^2)/m/2;
    end
end
figure,
surf(theta0,theta1,Jval');
axis([-3,3,-1,1]);
xlabel('\theta_0'); ylabel('\theta_1')

程序的运行结果如下图所示，红色的直线为最终的优化结果。

五、结语

以上Matlab代码中用到的数据和代码文件可以在这里下载，提取码m806
　　这一章主要介绍了单变量线性回归算法的原理和实现过程，定义了相关的符号和函数。随着后续的学习我们会发现，本章所学习的算法过程也是其它监督学习算法的大致实现过程，即定义特征量和假设函数、推导代价函数和偏导数、应用梯度下降算法寻找最优参数。
　　下一章介绍多变量的线性回归算法的原理和实现方法。