机器学习笔记(十四)——线性回归及其两种常用的优化方法

何为回归

回归的目的是预测数值型的目标值，最直接的办法是依据输入写出一个目标值的计算公式，比如要计算一个男生可以找到女朋友的概率：

$$P=0.4\ast 财产+0.3\ast 长相+0.3\ast 身高$$

这意味着要综合财产、长相、身高三个因素来判断概率，其中财产也是最重要的因素。

这个式子就可以被称作回归方程，其中0.4和0.3也被称作回归系数，一般来说回归系数都是未知的，我们通过输入数据求得回归系数的过程就是回归，得到回归系数之后，就可以通过公式得到最后的预测值。

这里给出的例子属于线性回归，而回归还有另一种较为复杂的形式——非线性回归，本文只介绍线性回归的相关知识。

线性回归

回归系数推导

线性回归(LR)可分为简单一元线性回归和多元线性回归，也就是我们平时接触的一次线性方程和多次线性方程，二者的主要区别也就是未知项的个数，为了便于理解，这里主要利用一元线性回归为例。

一元线性方程的公式应该是非常熟悉的：
$$y=wx+b$$

如果将输入数据都存放在矩阵X中，而回归系数都存放在向量$\omega$中，这样就可以得到矩阵形式的表达式：
$$y=X^T\omega$$

现在的问题是如何找到$\omega$，我们已经知道了如何度量一个分类器的性能，而回归模型的性能通常使用度量方法是“均方误差”，我们可以利用这个公式找到误差最小时的$\omega$，这里的误差是指预测值与真实值之间的差值。

均方误差的表示形式如下：
$$E=\sum _{i=0}^m(y_i-x_i^T\omega {)}^2$$

可以将其转化为矩阵形式：
$$E=(y-X\omega {)}^T(y-X\omega )$$

对$\omega$求导得出下面式子，并令其等于0：
$$2X^T(X\omega -y)=0$$

最后解出$\omega$如下：
$${\omega }^{\ast }=(X^{T}X{)}^{-1}{X}^{T}y$$

可以看到这里涉及到了对矩阵求逆，所以这个公式只有在可逆矩阵中才适用，就是说只有$X^{T}X$为满秩矩阵时，上述等式才成立。矩阵中非零行的个数定义为这个矩阵的秩, 记为R(A),对于矩阵$A_{n\times n}$，若R(A)=n，则称A为满秩矩阵。

线性拟合

现在有上图这样一个数据集，按照上文所提及的方法计算出最佳拟合直线的回归系数，就可以获得这个数据集的回归曲线。
这部分代码如下：

def standRegres(xMat,yMat):
    # 根据公式计算回归系数
    xTx = xMat.T * xMat
    if np.linalg.det(xTx) == 0.0:
        print("矩阵为奇异矩阵,不能求逆")
        return
    ws = xTx.I * (xMat.T*yMat)
    return ws
def plotRegression():
    xMat, yMat = loadDataSet('regression_data.txt')
    ws = standRegres(xMat, yMat) #计算回归系数
    xCopy = xMat.copy() #拷贝一个xMat矩阵
    xCopy.sort(0) # 对xCopy排序，方便绘图
    yHat = xCopy * ws # 计算对应的y值
    fig = plt.figure()
    ax = fig.add_subplot(111)
    ax.plot(xCopy[:, 1], yHat, c = 'red') # 绘制回归曲线
    plt.scatter(xMat.A[:,1],yMat.A,c = 'b',s = 10) # 绘制数据集的样本
    plt.show()

这里省略了加载数据集函数，第一个函数standRegres用来计算回归系数，先将x和y以矩阵形式传入，然后计算$X^{T}X$，接下来这个操作可能会比较陌生，因为上文说过了只有可逆矩阵才能求逆，这里linalg.det()方法可以计算矩阵的行列式，在行列式非零情况下，计算回归系数ws并返回。

第二个函数是绘制函数，在第一个函数计算出的回归系数基础上绘制回归曲线，最后绘制图像如下：

几乎任一数据集都可以用上述方法建立一个模型，那么这些模型的好坏程度如何评断呢？“相关系数”就可以计算预测值序列和真实值序列的匹配程度，Numpy中corrcoef方法就刚好可以计算出两个序列的相关性。

图中主对角线上的1代表自己与自己完全匹配，而次对角线则代表预测值和真实值间的相关性。

局部加权线性回归

上面拟合的回归曲线难免有些粗糙，而且会有一些欠拟合的现象，比如波浪处到直线的距离还是有点远，这很难取得最好的预测效果，而局部加权线性回归(LWLR)通过在估计值中引入一些偏差，从而降低预测的均方误差。

这种方法的基本思想就是给待预测点附近的每个点赋予一定的权重，将这些权重用一个新的矩阵W存储，表现形式如下：
$${\omega }^{\ast }=(X^{T}WX{)}^{-1}{X}^{T}Wy$$

加权模型也会认为样本点之间距离越近，越有可能符合同一个线性模型，所以对于当前预测样本点$x^{(i)}$来说，离它越近的样本点$x$将会被赋予更大的权重，利用高斯核就可实现这种机制：
$$\omega (i,i)=exp\left(\frac{x^{(i)}-x}{-2k^2}\right)$$

通过这种方式构建了一个只含对角元素的权重矩阵W，上述公式中只包含了一个需要调节的参数k，它决定了对附近的点赋予多大权重。

实际上当k值越大，参与训练的样本点也就越多；反之k值越小，只会有很少的局部点会参与训练。相应的k值过大可能就会出现欠拟合现象，k值过小可能就会出现过拟合现象，所以，k值的选择决定了模型的好坏。

具体代码如下：

def LWLR(testMat,xMat, yMat,k = 1.0):
    m = xMat.shape[0]
    n = testMat.shape[0]
    weights = np.mat(np.eye(m)) # 创建一个单位对角矩阵
    yHat = np.zeros(n) # 创建一个用来存预测值的矩阵
    for i in range(n):
        for j in range(m): # 遍历，根据公式求权重
            diffMat = testMat[i] - xMat[j]
            weights[j, j] = np.exp(diffMat * diffMat.T / (-2.0 * k ** 2))
        xTx = xMat.T * (weights * xMat)
        if np.linalg.det(xTx) == 0.0:
            print("矩阵为奇异矩阵,不能求逆")
            return
        ws = xTx.I * (xMat.T * (weights * yMat))
        yHat[i]=testMat[i]*ws # 求出预测值
    return yHat

这里先初始化了一个权重矩阵weights，形式是主对角线上都为1，其他位置都为0。接着遍历数据集，计算每个样本点对应的权重值，当样本点与待预测点距离越来越远时，权重将会衰减，而k控制衰减的速度。

这里需要注意的是testMat正是xMat本身，因为实际上每个样本的权重是数据集内的样本点之间相互计算得出的，用testMat命名只是方便区分而已，最后可绘制出不同k值对应的回归曲线如下：

可以看到当k=1.0时和普通的回归曲线没有什么差别；当k=0.01时回归直线拟合的就比较不错了；当k=0.002时回归曲线开始出现棱角，证明曲线的部分受其附近样本点影响很大，导致了过拟合的现象。

随着拟合变准确的同时，该模型也付出了相应的代价，即增加了计算量，因为在对每个点预测的同时都要使用整个数据集，下面将会介绍可以解决该问题的方法。

岭回归

现实生活中$X^{T}X$往往不是满秩矩阵，例如在数据集中可能遇到非常多的特征，其数目甚至超过了样本个数，导致X的列数多于行数，此时$X^{T}X$显然不满秩，为了解决这个问题，就引入了岭回归(ridge regression)的概念。

岭回归的思想非常简单，就是通过引入一个矩阵$\lambda I$，并且将这个单位矩阵和$X^{T}X$相加，从而将$X^{T}X$转化成一个可逆矩阵，进而可以对$X^{T}X+\lambda I$求逆，在这种情况下回归系数表达式就可以写成：$${\omega }^{\ast }=(X^{T}X+\lambda I{)}^{-1}{X}^{T}y$$

这里$I$是一个单位矩阵，即主对角线上元素为1，其余元素都为0。$\lambda$称为惩罚项，它的作用就是减少不重要的参数，这个技术被称作缩减，很明显缩减后能取得更好的预测结果。

在使用岭回归和缩减技术之前，需要对特征做标准化处理，这部分我们曾在KNN中也使用过，它的目的就是让每个特征具有相同的重要性，这部分代码如下：

yMean = np.mean(yMat, axis = 0) # 求真实y均值
yMat = yMat - yMean # 真实y值减去均值
xMeans = np.mean(xMat, axis = 0)# 求x均值
xVar = np.var(xMat, axis = 0) #求方差
xMat = (xMat - xMeans) / xVar # 实现标准化

代码的剩余部分和局部加权回归相似，所以就不再展示了，为了方便我们所用数据集还是上文提及的，这里最好使用特征数比样本数多的数据集，但我们只为展示得出的结论，如下：

这张图绘制了回归系数与$log(\lambda )$的关系，当$\lambda$非常小时，得到的系数是和最初线性回归一致的；当$\lambda$达到一定值时，系数全部缩减成0；所以在中间部分的某值将会取得最好的预测结果。

这里有一条恒为0的蓝色线是因为我们开始时设定了一个特征全为1，所以$\lambda$对这个特征是没有影响的。若想找到最佳的参数值还需要进行交叉验证，即多次调整参数寻求最优，并且根据每个特征对应系数的大小，也能判断出这个特征对预测结果影响的大小。

缩减方法除了岭回归之外还有lasso、LAR、PCA回归以及子集选择等，这些方法都可以提高预测精准率，感兴趣的伙伴可以自己了解一下。

总结

普通线性回归虽然可以拟合回归曲线，但是过于粗糙，有欠拟合的现象，通过局部加权回归可以得到更好的预测结果，但调参k是关键，k过小可能会出现过拟合的现象。针对现实任务中总出现的不可逆矩阵，缩减技术中的岭回归可以很好的解决该问题，它的主要思想就是通过消除多余的特征降低预测误差。

公众号【喵说Python】回台回复“线性回归”可获取源码供参考，感谢阅读