机器学习算法之线性回归

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

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

 

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

 

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

 

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

对于线性回归算法，我们下网从训练数据中学习到线性回归方程，即：

 

其中，b称为偏置，为回归系数。对于上式，令则上式可以表示为：

 

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

 

其中，为预测值，且：

平方损失函数为：

 

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

 

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

最小二乘法

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

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

 

 

 

 

 

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

 

因此求最小值即可通过对 求偏导数获得，并得一阶倒数得值为0，因此：

 

 

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

 

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

最小二乘解法

对于线性回归模型，假设训练集中有m个训练样本，每个训练样本中有n-1个特征，可以使用矩阵的表示方法，预测函数可以表示为：，其损失函数表示为： 。

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

 

令其为0，得到：

 

Python实现：

def least_square(feature, label):
    '''最小二乘法
    input:  feature(mat):特征
            label(mat):标签
    output: w(mat):回归系数
    '''
    w = (feature.T * feature).I * feature.T * label
    return w

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

梯度下降法

什么是梯度下降法

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

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

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

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

1. 当，自己设置初始点 ，设置步长（也叫做学习率），设置迭代终止的误差忍耐度tol。
2. 计算目标函数在点上的梯度 。
3. 计算 ，公式如下：
4. 计算梯度 。如果 则迭代停止，最优解的取值为 ；如果梯度的二范数大于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）的数据集。

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

 

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

 

其中，m是训练样本（training examples）的数量。

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

优点：全局最优解；易于并行实现；

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

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

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

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

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

import numpy as np
 
 
def gradient_descent(alpha, x, y, ep=0.0001, max_iter=10000):
    converged = False
    iter = 0
    m = x.shape[0]  # number of samples  
 
    # initial theta  
    t0 = np.random.random(x.shape[1])
    t1 = np.random.random(x.shape[1])
 
    # total error, J(theta)  
    J = sum([(t0 + t1 * x[i] - y[i]) ** 2 for i in range(m)])
 
    # Iterate Loop  
    while not converged:
        # for each training sample, compute the gradient (d/d_theta j(theta))  
        grad0 = 1.0 / m * sum([(t0 + t1 * x[i] - y[i]) for i in range(m)])
        grad1 = 1.0 / m * sum([(t0 + t1 * x[i] - y[i]) * x[i] for i in range(m)])
        # update the theta_temp  
        temp0 = t0 - alpha * grad0
        temp1 = t1 - alpha * grad1
 
        # update theta  
        t0 = temp0
        t1 = temp1
 
        # mean squared error  
        e = sum([(t0 + t1 * x[i] - y[i]) ** 2 for i in range(m)])
 
        if abs(J - e) <= ep:
            print('Converged, iterations: ', iter, '!!!')
            converged = True
 
        J = e  # update error   
        iter += 1  # update iter  
 
        if iter == max_iter:
            print('Max interactions exceeded!')
            converged = True
 
    return t0, t1

小批量梯度下降法MBGD

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

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

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

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

随机梯度下降法SGD

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

 

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

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

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

优点：训练速度快；

缺点：准确度下降，并不是全局最优；不易于并行实现。

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

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

牛顿法

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

基本牛顿法

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

 

对上式求导并令其为0，则为：

 

即得到：

 

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

基本牛顿法的流程

1. 给定终止误差值，初始点，令k=0；
2. 计算 ，若 ，则停止，输出 ；
3. 计算 ，并求解线性方程组 得解 ；
4. 令 ，并转2。

全局牛顿法

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

1. 给定终止误差值， ， ，初始点，令k=0；
2. 计算 ，若 ，则停止，输出 ；
3. 计算 ，并求解线性方程组 得解 ；
4. 设是不满足下列不等式的最小非负整数m：
5. 令 ，并转2。

全局牛顿法的具体实现：

def newton(feature, label, iterMax, sigma, delta):
    '''牛顿法
    input:  feature(mat):特征
            label(mat):标签
            iterMax(int):最大迭代次数
            sigma(float), delta(float):牛顿法中的参数
    output: w(mat):回归系数
    '''
    n = np.shape(feature)[1]
    w = np.mat(np.zeros((n, 1)))
    it = 0
    while it <= iterMax:
        # print it
        g = first_derivativ(feature, label, w)  # 一阶导数
        G = second_derivative(feature)  # 二阶导数
        d = -G.I * g
        m = get_min_m(feature, label, sigma, delta, d, w, g)  # 得到最小的m
        w = w + pow(sigma, m) * d
        if it % 10 == 0:
            print "\t---- itration: ", it, " , error: ", get_error(feature, label , w)[0, 0]
        it += 1       
    return w

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

def get_min_m(feature, label, sigma, delta, d, w, g):
    '''计算步长中最小的值m
    input:  feature(mat):特征
            label(mat):标签
            sigma(float),delta(float):全局牛顿法的参数
            d(mat):负的一阶导数除以二阶导数值
            g(mat):一阶导数值
    output: m(int):最小m值
    '''
    m = 0
    while True:
        w_new = w + pow(sigma, m) * d
        left = get_error(feature, label , w_new)
        right = get_error(feature, label , w) + delta * pow(sigma, m) * g.T * d
        if left <= right:
            break
        else:
            m += 1
    return m

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

def get_error(feature, label, w):
    '''计算误差
    input:  feature(mat):特征
            label(mat):标签
            w(mat):线性回归模型的参数
    output: 损失函数值
    '''
    return (label - feature * w).T * (label - feature * w) / 2

Armijo搜索

给定 ， ，令步长因子 ，其中是满足下列不等式的最小非负整数：

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

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

 

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

 

损失函数的二阶导数为：

 

相关代码：

def first_derivativ(feature, label, w):
    '''计算一阶导函数的值
    input:  feature(mat):特征
            label(mat):标签
    output: g(mat):一阶导数值
    '''
    m, n = np.shape(feature)
    g = np.mat(np.zeros((n, 1)))
    for i in xrange(m):
        err = label[i, 0] - feature[i, ] * w
        for j in xrange(n):
            g[j, ] -= err * feature[i, j]
    return g     
 
def second_derivative(feature):
    '''计算二阶导函数的值
    input:  feature(mat):特征
    output: G(mat):二阶导数值
    '''
    m, n = np.shape(feature)
    G = np.mat(np.zeros((n, n)))
    for i in xrange(m):
        x_left = feature[i, ].T
        x_right = feature[i, ]
        G += x_left * x_right
    return G

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

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

局部加权线性回归

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

 

M为给每个点的权重。

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

 

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

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

import numpy as np
 
def lwlr(feature, label, k):
    '''局部加权线性回归
    input:  feature(mat):特征
            label(mat):标签
            k(int):核函数的系数
    output: predict(mat):最终的结果
    '''
    m = np.shape(feature)[0]
    predict = np.zeros(m)
    weights = np.mat(np.eye(m))
    for i in xrange(m):
        for j in xrange(m):
            diff = feature[i, ] - feature[j, ]
            weights[j,j] = np.exp(diff * diff.T / (-2.0 * k ** 2))
        xTx = feature.T * (weights * feature)
        ws = xTx.I * (feature.T * (weights * label))
        predict[i] = feature[i, ] * ws
    return predict

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

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

from sklearn.datasets import load_boston
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LinearRegression
 
boston = load_boston()
X = boston.data
y = boston.target
 
print(X.shape)  # 样本个数及特征数
print(boston.feature_names)  # 特征标签名字
 
# 把数据分为训练数据集和测试数据集(20%数据作为测试数据集）
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=3)
 
model = LinearRegression()
model.fit(X_train, y_train)
 
train_score = model.score(X_train, y_train)  # 模型对训练样本得准确性
test_score = model.score(X_test, y_test)  # 模型对测试集的准确性
print(train_score)
print(test_score)

模型优化

多项式于线性回归

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

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

数据归一化

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

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

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

优化后代码：

from sklearn.datasets import load_boston
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LinearRegression
from sklearn.preprocessing import PolynomialFeatures
from sklearn.pipeline import Pipeline
 
boston = load_boston()
X = boston.data
y = boston.target
 
print(X.shape)  # 样本个数及特征数
print(boston.feature_names)  # 特征标签名字
 
# 把数据分为训练数据集和测试数据集(20%数据作为测试数据集）
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=3)
 
 
def polynomial_model(degree=1):
    polynomial_features = PolynomialFeatures(degree=degree, include_bias=False)
    linear_regression = LinearRegression(normalize=True)
    pipeline = Pipeline([("polynomial_features", polynomial_features), ("linear_regression", linear_regression)])
    return pipeline
 
 
model = polynomial_model(degree=2)
model.fit(X_train, y_train)
 
train_score = model.score(X_train, y_train)  # 模型对训练样本得准确性
test_score = model.score(X_test, y_test)  # 模型对测试集的准确性
print(train_score)
print(test_score)

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