机器学习算法之线性回归

3 min read

线性回归是统计学总最常用的算法之一。从根本上来说,当你想表示两个变量间数学关系时,就可以使用线性回归。当你使用它时,你首先假设输出变量(有时称为响应变量、因变量或标签)和预测变量(有时称为自变量、解释变量或特征)之间存在线性关系。当然这种线性关系也可能存在于一个输出变量和数个预测变量之间。输出变量于预测变量之间存在线性关系是一个大胆的假设,同时也是一个最简单的假设。从数学表示形式来看,线性函数比非线性函数更加简单。线性模型作为最简单的参数化方法,始终值得关注。这是因为很多问题,甚至本质是非线性的问题,也可以采用线性模型解决。

在线性回归中,数据采用线性函数的方式进行数据建模,对于模型中的未知参数也采用数据进行估计。对于一个多变量的线性回归模型可以表示为如下公式:

    \[Y_i = \beta _0 + \beta _1X_{i1}+\beta _2X_{i2}+...+\beta _pX_{ip}+\varepsilon _i, i = 1,2,...,n\]

其中,Y是X的线性函数,\varepsilon _i是误差项。线性则是Y_i的条件均值,在参数\beta里是线性的。有些函数模型看起来并不一定是线性回归 但是通过代数转换可以转换为线性回归模型。对于一个简单的线性回归,可以表示为如下:

    \[Y_i = \alpha +\beta x_i\]

对于单变量的线性回归,是目前最容易理解的线性模型,其回归式如下:

    \[Y = \alpha +\beta x + \varepsilon\]

线性回归是回归分析中最为广泛使用的模型, 在结果预测及函数关系的应用中较为频繁。

对于线性回归算法,我们希望从训练数据中学习到线性回归方程,即:

    \[y = b + \sum_{i=1}^{n}w_ix_i\]

其中,b称为偏置,w_i为回归系数。对于上式,令x_0=1则上式可以表示为:

    \[y = \sum_{i=0}^{n}w_ix_i\]

在线性回归模型中,其目的是 求出线性模型回归方程,即求出线性回归方程中的回归系数w_i。线性回归的评价是指如何度量预测值(Prediction)于标签(Label)之间的接近程序,线性回归模型的损失函数可以是绝对损失(Absolute Loss)或者平方损失(Squared Loss)。其中,绝对损失函数为:

    \[l = |y-\hat{y}|\]

其中,\hat{y}为预测值,且:\hat{y} = \sum_{i=0}^{n}w_ix_i

平方损失函数为:

    \[l =(y- \hat{y})^2\]

由于平方损失处处可导,通常使用平方误差作为线性回归模型的损失函数。假设有m个训练样本,每个样本中有n-1个特征,则平方误差可以表示为:

    \[l = \frac{1}{2}\sum_{i=1}^{m}(y^i-\sum_{j=0}^{n-1}w_jx_j^i)^2\]

对于如上的损失函数,线性回归的求解是希望求得平方误差最小值。

最小二乘法

最小二乘法也称作最小平方法, 最常用的是普通最小二乘法(Ordinary Least Square) ,它是一种数学中的优化方法,试图找到一个或一组估计值,使得实际值与估计值尽可能相似, 距离最小,目的是通过已有的数据来预测未知数据。一般通过一条多元一次的直线方程来计算,在二维坐标中即二元一次方程。 例如,在二维坐标中, 非常多的点分散在其中,试图绘制一条直线,使得这些分散的点到直线上的距离最小。这里的距离最小并非点到直线的垂直距离最短,而是点到直接的y轴距离最短,即通过该点并与y轴平行的直线,点到该y轴平行线与直线交点的距离最短。

最小二乘法的核心思想是通过最小化误差的平方和,试图找到最可能的函数方程。 例如,在二维坐标系中存在五个数据点(10,20), (11,23), (12,25), (13,27), (14,26), 希望找出一条该五个点距离最短的直线,根据二元一次方程:y = ax + b。将五个点分别带入该二元方程得到如下:

    \[20=10a+b\]

    \[23=11a+b\]

    \[25=12a+b\]

    \[27=13a+b\]

    \[26=14a+b\]

由于最小二乘法是尽可能使得等号两边的方差值最小,因此:

    \[S(a,b) = [20-(10a+b)]^2 + [23-(11a+b)]^2+[25-(12a+b)]62+[27-(13a+b)]^2+[26-(14a+b)]^2\]

因此求最小值即可通过对S(a,b)求偏导数获得,并得一阶倒数得值为0,因此:

    \[\frac{\partial S }{\partial a}= 1460a+120b-2936\]

    \[\frac{\partial S }{\partial b}= 12a + 10b -242\]

即得到关于求解未知变量a、b的二元一次方程:

    \[\left\{\begin{matrix}1460a+120b = 2936\\ 12a+10b = 242\end{matrix}\right.\]

通过计算上述二元一次方法即得到a=0.0243, b=24.1708。 因此,在上述五个点中,通过最小二乘法得到直线方程:y=24.l708+0.0243x是使得五个点到该直线距离最小的直线。

最小二乘解法

对于线性回归模型,假设训练集中有m个训练样本,每个训练样本中有n-1个特征,可以使用矩阵的表示方法,预测函数可以表示为:Y=XW,其损失函数表示为:(Y-XW)^T(Y-XW)

其中,标签Y为m×1的矩阵,训练特征X为m×n的矩阵,回归系数W 为n×1的矩阵。在最小二乘法中,对W 求导,即:

    \[\frac{d}{dW}(Y-XW)^T(Y-XW) = X^T(Y-XW)\]

令其为0,得到:

    \[\hat{W} = (X^TX)^{-1}X^TY\]

Python实现:

最小二乘法虽然看似是一个直线方程的问题,但是在实际应用中却应用非常广泛,因为它得到的方程可以视为一个函数模型,该函数模型可以对后续的工作带来极大的便利。上述过程均是通过线性问题的求解方式进行阐述,但是在更多时候,需要解决的问题不是一个线性问题,它需要通过多项式拟合的方式进行处理,但是原理和求解方式均一致。当样本数据不断增加后,计算量会明显增加,在阶数更高时计算扯则更为复杂,为解决更多问题,基千最小二乘法衍生出了移动最小二乘法、 加权最小二乘法及偏最小二乘法等。除最小二乘法外, 梯度下降也是线性回归中的方法之一。

梯度下降法

什么是梯度下降法

梯度下降法的基本思想可以类比为一个下山的过程。假设这样一个场景:一个人被困在山上,需要从山上下来(i.e. 找到山的最低点,也就是山谷)。但此时山上的浓雾很大,导致可视度很低。因此,下山的路径就无法确定,他必须利用自己周围的信息去找到下山的路径。这个时候,他就可以利用梯度下降算法来帮助自己下山。具体来说就是,以他当前的所处的位置为基准,寻找这个位置最陡峭的地方,然后朝着山的高度下降的地方走,同理,如果我们的目标是上山,也就是爬到山顶,那么此时应该是朝着最陡峭的方向往上走。然后每走一段距离,都反复采用同一个方法,最后就能成功的抵达山谷。

在微积分里面,对多元函数的参数求\partial偏导数,把求得的各个参数的偏导数以向量的形式写出来,就是梯度。比如函数f(x,y),分别对x,y求偏导数,求得的梯度向量就是(\frac{\partial f }{\partial x},\frac{\partial f }{\partial y})^T,简称grad f(x,y)或者\bigtriangledown f(x,y)。对于在点(x_0,y_0)的具体梯度向量就是(\frac{\partial f }{\partial x_0},\frac{\partial f }{\partial y_0})^T,或者\bigtriangledown f(x_0,y_0)。如果是3个参数的向量梯度,就是(\frac{\partial f }{\partial x},\frac{\partial f }{\partial y},\frac{\partial f }{\partial z})^T,以此类推。

那么这个梯度向量求出来有什么意义呢?他的意义从几何意义上讲,就是函数变化增加最快的地方。具体来说,对于函数f(x,y),在点(x_0,y_0),沿着梯度向量的方向就是(\frac{\partial f }{\partial x_0},\frac{\partial f }{\partial y_0})^T的方向是f(x,y)增加最快的地方。或者说,沿着梯度向量的方向,更加容易找到函数的最大值。反过来说,沿着梯度向量相反的方向,也就是-(\frac{\partial f }{\partial x_0},\frac{\partial f }{\partial y_0})^T的方向,梯度减少最快,也就是更加容易找到函数的最小值。

对于凸优化问题来说,导数为0(梯度为0向量)的点,就是优化问题的解。为了找到这个解,我们沿着梯度的反方向进行线性搜索,每次搜索的步长为某个特定的数值\alpha,直到梯度与0向量非常接近为止。上面描述的这个方法就是梯度下降法。迭代算法的步骤如下:

  1. i=0,自己设置初始点\textbf{x}^0=(x_1^0,x_2^0,\cdots, x_n^0),设置步长(也叫做学习率)\alpha,设置迭代终止的误差忍耐度tol。
  2. 计算目标函数f在点x_i上的梯度\nabla f_{\textbf{x}^i}
  3. 计算\textbf{x}^{i+1},公式如下:\textbf{x}^{i+1} = \textbf{x}^{i} - \alpha \nabla f_{\textbf{x}^i}
  4. 计算梯度\nabla f_{\textbf{x}^{i+1}}。如果\|\nabla f_{\textbf{x}^{i+1}}\|_2\leq tol则迭代停止,最优解的取值为\textbf{x}^{i+1};如果梯度的二范数大于tol,那么i=i+1,并返回第(3)步。

梯度下降法的类型

基于如何使用数据计算代价函数的导数,梯度下降法可以被定义为不同的形式(various variants)。确切地说,根据使用数据量的大小(the amount of data),时间复杂度(time complexity)和算法的准确率(accuracy of the algorithm),梯度下降法可分为:

  • 批量梯度下降法(Batch Gradient Descent, BGD);
  • 随机梯度下降法(Stochastic Gradient Descent, SGD);
  • 小批量梯度下降法(Mini-Batch Gradient Descent, MBGD)。

批量梯度下降法BGD

这是梯度下降法的基本类型,这种方法使用整个数据集(the complete dataset)去计算代价函数的梯度。每次使用全部数据计算梯度去更新参数,批量梯度下降法会很慢,并且很难处理不能载入内存(don’t fit in memory)的数据集。

批量梯度下降法的损失函数为:

    \[\frac{\partial J(\theta) }{\partial {\theta}_j} = \sum_{i=1}^{m}(h_\theta(x^i)-y^i )x_j^i\]

进一步得到批量梯度下降的迭代式为:

    \[{\theta}'_j = \theta _j - \alpha \frac{\partial J(\theta )}{\partial \theta _j} = \theta _j = \alpha \sum_{i=1}^{m}(h_\theta (x^i)-y^i)x_j^i\]

其中,m是训练样本(training examples)的数量。

每迭代一步,都要用到训练集所有的数据,如果样本数目很大,那么可想而知这种方法的迭代速度!

优点:全局最优解;易于并行实现;

缺点:当样本数目很多时,训练过程会很慢。

  • 如果训练集有3亿条数据,你需要从硬盘读取全部数据到内存中;
  • 每次一次计算完求和后,就进行参数更新;
  • 然后重复上面每一步;
  • 这意味着需要较长的时间才能收敛;
  • 特别是因为磁盘输入/输出(disk I/O)是系统典型瓶颈,所以这种方法会不可避免地需要大量的读取。

从迭代的次数上来看,BGD迭代的次数相对较少。其迭代的收敛曲线示意图可以表示如下:

上图是每次迭代后的等高线图,每个不同颜色的线表示代价函数不同的值。运用梯度下降会快速收敛到圆心,即唯一的一个全局最小值。

下面的Python代码实现了批量梯度下降法:

小批量梯度下降法MBGD

小批量梯度下降法是最广泛使用的一种算法,该算法每次使用m个训练样本(称之为一批)进行训练,能够更快得出准确的答案。小批量梯度下降法不是使用完整数据集,在每次迭代中仅使用m个训练样本去计算代价函数的梯度。一般小批量梯度下降法所选取的样本数量在50到256个之间,视具体应用而定。

  • 减少了参数更新时的变化,能够更加稳定地收敛。
  • 能利用高度优化的矩阵,进行高效的梯度计算。

随机初始化参数后,按如下伪码计算代价函数的梯度:

这里b表示一批训练样本的个数,m是训练样本的总数。

随机梯度下降法SGD

批量梯度下降法被证明是一个较慢的算法,所以,我们可以选择随机梯度下降法达到更快的计算。随机梯度下降法的第一步是随机化整个数据集。在每次迭代仅选择一个训练样本去计算代价函数的梯度,然后更新参数。即使是大规模数据集,随机梯度下降法也会很快收敛。随机梯度下降法得到结果的准确性可能不会是最好的,但是计算结果的速度很快。在随机化初始参数之后,使用如下方法计算代价函数的梯度:

    \[{\theta }'_j = \theta _j + \frac{1}{m}(y^i-h_\theta (x^i))x_j^i\]

批量梯度下降会最小化所有训练样本的损失函数,使得最终求解的是全局的最优解,即求解的参数是使得风险函数最小。随机梯度下降会最小化每条样本的损失函数, 虽然不是每次迭代得到的损失函数都向着全局最优方向, 但是大的整体的方向是向全局最优解的,最终的结果往往是在全局最优解附近。

随机梯度下降法的优化过程为:

随机梯度下降是通过每个样本来迭代更新一次,如果样本量很大的情况(例如几十万),那么可能只用其中几万条或者几千条的样本,就已经将theta迭代到最优解了,对比上面的批量梯度下降,迭代一次需要用到十几万训练样本,一次迭代不可能最优,如果迭代10次的话就需要遍历训练样本10次。但是,SGD伴随的一个问题是噪音较BGD要多,使得SGD并不是每次迭代都向着整体最优化方向。

优点:训练速度快;

缺点:准确度下降,并不是全局最优;不易于并行实现。

从迭代的次数上来看,SGD迭代的次数较多,在解空间的搜索过程看起来很盲目。其迭代的收敛曲线示意图可以表示如下:

参考链接:http://sofasofa.io/tutorials/python_gradient_descent/index.php

牛顿法

除了前面说的梯度下降法,牛顿法也是机器学习中用的比较多的一种优化算法。牛顿法的基本思想是利用迭代点x_k处的一阶导数(梯度)和二阶导数(Hessen矩阵)对目标函数进行二次函数近似,然后把二次模型的极小点作为新的迭代点,并不断重复这一过程,直至求得满足精度的近似极小值。牛顿法下降的速度比梯度下降的快,而且能高度逼近最优值。牛顿法分为基本的牛顿法和全局牛顿法。

基本牛顿法

基本牛顿法是一种基于导数的算法,它每一步的迭代方向都是沿着当前点函数值下降的方向。对于一维的情形,对于一个需要求解的优化函数f(x),求函数的极值的问题可以转化为求导函数{f}'(x)=0。对函数f(x)进行泰勒展开到二阶,得到:

    \[f(x)=f(x_k) + {f}'(x_k)(x-x_k) + \frac{1}{2}f^n(x_k)(x-x_k)^2\]

对上式求导并令其为0,则为:

    \[{f}'(x) + {f}''(x_k)(x-x_k) = 0\]

即得到:

    \[x = x_k - \frac{{f}'(x_k)}{{f}''(x_k)}\]

这就是牛顿法的更新公式。

基本牛顿法的流程

  1. 给定终止误差值0 \leq \varepsilon \leq 1,初始点x_0 \in \mathbb{R}^n ,令k=0;
  2. 计算g_k = \triangledown f(x_k),若\left \| g_k \right \| \leq \varepsilon,则停止,输出x^* \approx x_k
  3. 计算G_k = \triangledown ^2 f(x_k),并求解线性方程组G_kd=-g_k得解d_k
  4. x_{k+1} = x_k + d_k, k=k+1,并转2。

全局牛顿法

牛顿法最突出的优点是收敛速度快,具有局部二阶收敛性,但是,基本牛顿法初始点需要足够“靠近”极小点,否则,有可能导致算法不收敛,此时就引入了全局牛顿法。全局牛顿法的流程为:

  1. 给定终止误差值0 \leq \varepsilon \leq 1\delta \in (0,1)\sigma \in (0,0.5),初始点x_0 \in \mathbb{R}^n ,令k=0;
  2. 计算g_k = \triangledown f(x_k),若\left \| g_k \right \| \leq \varepsilon,则停止,输出x^* \approx x_k
  3. 计算G_k = \triangledown ^2 f(x_k),并求解线性方程组G_kd=-g_k得解d_k
  4. m_k是不满足下列不等式的最小非负整数m:f(x_k+\delta ^md_k)\leq f(x_k)+\sigma \delta ^mg_k^Td_k
  5. \alpha _k = \delta ^{m_k},x_{k+1}=x_k+\alpha _kd_k,k=k+1,并转2。

全局牛顿法的具体实现:

在程序中,函数newton利用全局牛顿法对线性回归模型中的参数进行学习,函数newton的输入为训练特征feature、训练的目标值label、全局牛顿法的最大迭代次数iterMax以及全局牛顿法的两个参数sigma和delta。函数newton的输出是线性回归模型的参数 w。在函数 newton 中需要计算损失函数的一阶导数,计算损失函数的二阶导数,,同时需要计算最小的m值,最终根据上述的值更新权重。

程序实现了全局牛顿法中最小m值的确定,在函数get_min_m中,其输入为训练数据的特征 feature,训练数据的目标值 label,全局牛顿法的参数 sigma、delta、d以及损失函数的一阶导数值g。其输出是最小的m值m。在计算的过程中,计算损失函数值时使用到了get_error函数,其具体实现:

Armijo搜索

给定\delta \in (0,1)\sigma  \in (0,0.5),令步长因子\alpha _k = \delta ^{m_k} ,其中m_k是满足下列不等式的最小非负整数:f(x_k+\delta ^md_k)\leq f(x_k)+\sigma \delta ^mg_k^Td_k

利用全局牛顿法求解线性回归模型

假设有m个训练样本,其中,每个样本有n-1个特征,则线性回归模型的损失函数为:

    \[l = \frac{1}{2}\sum_{i=1}^{m}(y^i-\sum_{j=0}^{n-1}w_jx_j^i)^2\]

若是利用全局牛顿法求解线性回归模型,需要计算线性回归模型损失函数的一阶导数和二阶导数,其一阶导数为:

    \[\frac{\partial l}{\partial w_j} = -\sum_{i=1}^{m}[(y^i-\sum_{j=0}^{n-1}w_j\cdot x_j^i)*x_j^i]\]

损失函数的二阶导数为:

    \[\frac{\partial l}{\partial w_j\partial w_k} = \sum_{i=1}^{m}(x_j^i\cdot x_k^i)\]

相关代码:

first_derivativ 实现了损失函数一阶导数值的求解。在函数first_derivativ 中,其输入为训练数据的特征feature 和训练数据的目标值 label,其输出为损失函数的一阶导数g,其中g是一个n×1的向量。

second_derivative函数实现了损失函数二阶导数值的计算。在函数second_derivative中,其输入为训练数据的特征feature,输出为损失函数的二阶导数G,其中G是一个n×n的矩阵。

局部加权线性回归

在线性回归中会出现欠拟合的情况,有些方法可以用来解决这样的问题。局部加权线性回归(LWLR)就是这样的一种方法。局部加权线性回归采用的是给预测点附近的每个点赋予一定的权重,此时的回归系数可以表示为:

    \[\hat{W}=(X^TMX)^{-1}X^TMY\]

M为给每个点的权重。

LWLR使用核函数来对附近的点赋予更高的权重,常用的有高斯核,对应的权重为:

    \[w(i,i)=\exp(\frac{\left | X^{(i)-x} \right |}{-2k^2})\]

这样的权重矩阵只含对角元素。

局部加权线性回归的具体实现:

使用Scikit-Learn进行线性回归分析

这里使用SkLearn自带得数据集进行测试:

模型优化

多项式与线性回归

当线性回归模型太简单导致欠拟合时,我们可以增加特征多项式来让线性回归模型更好的拟合数据。比如两个特征x_1,x_2,可以增加两特征的乘积x_1 \cdot x_2作为新特征x_3。我们还可以增加x_1^2作为另外一个新特征x_4

在scikit-learn里,线性回归是由类sklearn.linear_model.LinearRegression实现,多项式由类sklearn.preprocessing.PolynomialFeatures实现。那么要怎么添加多项式特征呢?我们需要用一个管道把两个类串起来,即用sklearn.pipeline.Pipeline把两个模型串起来。一个Pipeline可以包含多个处理节点,在scikit-learn里,除了最后一个节点外,其他的节点都必须实现’fit()’和’transform()’方法,最后一个节点只需要实现fit()方法即可。当训练样本数据送进Pipeline里进行处理时,它会逐个调用节点的fit()和transform()方法,然后调用最后一个节点的fit()方法来拟合数据。

数据归一化

当线性回归模型有多个输入特征时,特别是使用多项式添加特征时,需要对数据进行归一化处理。比如,特征x_1的范围在[1,4]之间,特征x_2的范围在[1,2000]之间,这种情况下,可以让x_1除以4来作为新特征x_1,同时让x_2除以2000来作为新特征x_2,该过程称为特征缩放(feature scaling)。可以使用特征缩放来对训练样本进行归一化处理,处理后的特征值范围在[0,1]之间。

归一化处理的目的是让算法收敛更快,提升模型拟合过程中的计算效率。进行归一化处理后,当有个新的样本需要计算预测值时,也需要先进行归一化处理,再通过模型来计算预测值,计算出来的预测值要再乘以归一化处理的系数,这样得到的数据才是真实的预测数据。

在scikit-learn里,使用LinearRegression进行线性回归时,可以指定normalize=True来对数据进行归一化处理。具体可查阅scikit-learn文档。

优化后代码:

使用二阶多项式训练后样本分数和测试分数都提高了,看来模型得到了优化,改为三阶后,针对训练样本的分数达到了1,而针对测试样本的得分数却是负数。说明这个模型已经过拟合。

打赏作者
微信支付标点符 wechat qrcode
支付宝标点符 alipay qrcode

使用GridSearchCV进行网格搜索

GridSearchCV简介 在机器学习模型中,需要人工选择的参数称为超参数。比如随机森林中决策树的个数,人工
49 sec read

PageRank算法学习与研究

什么是PageRank PageRank,简称PR,是Google排名运算法则(排名公式)的一部分,是Goog
2 min read

多经纬度坐标的中心点计算方法

在实际的应用场景,通常会遇到计算多个经纬度中心的需求。而在计算经纬度中心点通常有三种方式,每种方式对应不同的需
1 min read

发表评论

电子邮件地址不会被公开。 必填项已用*标注