Python手撸机器学习系列（五）：线性回归与岭回归（最小二乘法求解）

线性回归

一、原理

线性回归应该是最简单的机器学习算法了，我们以最小二乘法为基础来进行求解

1.1 数据输入

设数据表现形式为：

$Data=(x_1,y_1),(x_2,y_2),...,(x_N,y_N)$

其中$x_i\in R^p,y_i\in R,i=1,2,...,N$

就是说有N条数据，每条数据分为特征$x_i$和对应的标签值$y_i$，其中$x_i$为$p\ast 1$维，也就是说有$p$个特征，$x_i$表现为列向量

用$X$表示特征矩阵为
$$X=(x_1,x_2,..,x_N)=\left[\begin{array}{c}{x}_1^T\\ x_2^T\\ ⋮\\ x_N^T\end{array}\right]={\left[\begin{array}{cccc}{x}_{11}& x_{12}& \cdots & x_{1p}\\ x_{21}& x_{22}& \cdots & x_{2p}\\ ⋮& ⋮& \ddots & \cdots \\ x_{N1}& x_{N2}& \cdots & x_{Np}\end{array}\right]}_{N\times p}^T$$

用$Y$表示标签矩阵为：
$$Y={\left[\begin{array}{c}{y}_1\\ y_2\\ ⋮\\ y_N\end{array}\right]}_{N\times 1}$$
线性回归的目的在于拟合一个函数来根据特征值预测标签值，当数据为平面坐标点时（ 横坐标为x纵坐标为y），拟合的函数表现为一条坐标平面直线。

我们需要拟合函数的权重$W$，其维度为$p\times 1$，当数据为平面坐标点时为一个数（$1\times 1$的矩阵）
$$W={\left[\begin{array}{c}{W}_1\\ W_2\\ ⋮\\ W_N\end{array}\right]}_{p\times 1}$$

1.2 最小二乘法一般解法

显然，用最小二乘法求解，其损失函数为$MSE$ :
$$\begin{array}{rl}L(W)& =\sum _{i=1}^N||W^T{x}_i-y_i|{|}^2\\ & =\sum _{i=1}^N(W^T{x}_i-y_i{)}^2\\ & =\left[W^T{x}_1-y_1,W^T{x}_2-y_2,...,W^T{x}_N-y_N\right]\left[\begin{array}{c}{W}^T{x}_1-y_1\\ W^T{x}_2-y_2\\ ⋮\\ W^T{x}_N-y_B\end{array}\right]\\ & =(W^T{X}^T-Y^T)(XW-Y)\\ & =W^T{X}^{T}XW-W^T{X}^{T}Y-Y^{T}XW+Y^{T}Y\\ & =W^T{X}^{T}XW-2W^T{X}^{T}Y+Y^{T}Y\end{array}$$
使用损失函数对$W$求导：
$$\frac{\partial L(W)}{\partial W}=2X^{T}XW-2X^{T}Y=0$$
可得
$$\begin{array}{rl}{X}^{T}XW& =X^{T}Y\\ W& =(X^{T}X{)}^{-1}{X}^{T}Y\end{array}$$
即求得参数$W$，其中$(X^{T}X{)}^{-1}{X}^T$被称为$X$的伪逆矩阵

1.3 最小二乘法的几何理解

我们将$X$理解为一种特征空间，比如说，现在由$x_1、x_2$两个特征组成了一个平面，而使用$Y$向该平面投影：

投影可以理解为在向量$x_1、x_2$的线性组合，则可以找到一组$w$用来表示这段线性组合的权重，即$[x_1,x_2][w_1,w_2{]}^T$

扩展到全部的$X$和$Y$，投影就是$XW$，而根据向量的减法，垂直于平面的那一段可以描述为$Y-XW$

垂直的向量内积为0，即：
$$\begin{array}{rl}{X}^T(Y-XW)& =0\\ X^{T}Y& =X^{T}XW\\ W& =(X^{T}X{)}^{-1}{X}^{T}Y\end{array}$$
与之前使用求导算得的$W$一致

二、正则化与岭回归

对于上述问题，求得的线性回归损失函数及其$W$为：
$$\begin{array}{rl}L(W)& =\sum _{i=1}^N||W^T{x}_i-y_i|{|}^2\\ W& =(X^{T}X{)}^{-1}{X}^{T}Y\end{array}$$
而在实际应用中，$(X^{T}X)$可能是不可逆的，比如特征大于样本时，而且样本过少特征过多还容易造成过拟合，由此提出了正则化的解决方案。

正则化的框架为：
$$L(W)+\lambda P(W)$$
其中$\lambda P(w)$被称为惩罚项

一般正则化可以分为L1正则化和L2正则化：

• L1正则化：也被称作Lasso，此时$P(W)=||W|{|}_1$，即$W$的一范式
• L2正则化：也被称作Ridge和岭回归，此时$P(W)=||W|{|}_2^2=W^{T}W$

在这里我们主要介绍L2正则化，也就是岭回归

在岭回归的情况下，损失函数就变为：
$$\begin{array}{rl}J(W)& =\sum _{i=1}^N||W^T{x}_i-y_i|{|}^2+\lambda W^{T}W\\ & =(W^T{X}^T-Y^T)(XW-Y)+\lambda W^{T}W\\ & =W^T{X}^{T}XW-2W^T{X}^{T}Y+Y^{T}Y+\lambda W^{T}W\\ & =W^T(X^{T}X+\lambda I)W-2W^T{X}^{T}Y+Y^{T}Y\end{array}$$
对$W$进行求导：
$$\frac{\partial J(W)}{\partial W}=2(X^{T}X+\lambda I)W-2X^{T}Y=0$$
求得
$$W=(X^{T}X+\lambda I{)}^{-1}{X}^{T}Y$$
此时$(X^{T}X+\lambda I)$一定可逆，理由如下：

假设$X^{T}X$不可逆，即不满秩，其可以变换为某一行或某一列为全0的矩阵：
$$X^{T}X=\left[\begin{array}{ccccc}{a}_{11}& a_{12}& a_{13}& \cdots & a_{1n}\\ 0& a_{22}& a_{23}& \cdots & a_{2n}\\ 0& 0& a_{33}& \cdots & a_{3n}\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ 0& 0& 0& \cdots & 0\end{array}\right]$$

$$X^{T}X+\lambda I=\left[\begin{array}{ccccc}{a}_{11}+\lambda & a_{12}& a_{13}& \cdots & a_{1n}\\ 0& a_{22}+\lambda & a_{23}& \cdots & a_{2n}\\ 0& 0& a_{33}+\lambda & \cdots & a_{3n}\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ 0& 0& 0& \cdots & \lambda \end{array}\right]$$

此时不存在全为0的行或列，矩阵满秩可逆，当然，还是有特殊清况，比如：

• $\lambda =0$
• 原本矩阵中存在对象线上元素为$-\lambda$，其他元素都为0的行或列

当然，这些情况可以通过调节$\lambda$来解决

三、代码实现

梯度下降法，使用之前公式中损失函数$L$对$W$的求导作为梯度，每次乘上学习率进行衰减，即：

$d_W=2X^{T}XW-2X^{T}Y$

$W=W-lr\ast d_W$

公式法则是直接用最后的公式计算$X$的伪逆矩阵乘上$Y$

岭回归只需要在原始基础上加上惩罚项

其实两种方法实现起来几乎是一模一样的：

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
from sklearn.datasets import load_boston
from sklearn.linear_model import  LinearRegression

def get_points():
    #获取坐标点
    point = [[3,3],[4,3],[2,5],[6,7],[8,10]]
    X = np.array([x[0] for x in point])
    Y = np.array([x[1] for x in point])
    X = X.reshape(X.shape[0],1)
    Y = Y.reshape(Y.shape[0],1)

    #另一组点，此时x为二维，画图函数无效，记得注释掉
    # X = np.array([[1, 2], [3, 4], [3, 6]])
    # Y = np.array([3, 6, 9])
    # X = X.reshape(X.shape[0],2)
    # Y = Y.reshape(Y.shape[0],1)
    # return X,Y.reshape(Y.shape[0],1)
    return X,Y

def get_house():
    #波士顿房价数据集
    data,target = load_boston(return_X_y=True)
    return data,target.reshape(target.shape[0],1)

def cal_w(X,Y):
    #公式直接计算W
    w = np.linalg.inv(X.T.dot(X)).dot(X.T).dot(Y)
    return w

def Ridge_w(X,Y):
    #岭回归，可以通过调节lambda，此处为a来改变直线
    a = 5
    w = np.linalg.inv(X.T.dot(X)+a*np.eye(X.shape[1])).dot(X.T).dot(Y)
    return w

def gradient(X,Y):
    #梯度下降法
    lr = 0.01
    epochs = 2
    w = np.zeros((X.shape[1],1))
    for epoch in range(epochs):
        dw = X.T.dot(X).dot(w)-X.T.dot(Y)
        w -= lr*dw
    return w

if __name__ == '__main__':
    X,Y = get_points()
    w1 = gradient(X,Y)
    w2 = cal_w(X,Y)
    w3 = Ridge_w(X,Y)
    plt.scatter(X, Y, c='b',label = 'points')
    plt.plot(np.arange(1, 9), w1[0][0] * np.arange(1, 9), c='r',label = 'line_gradient')
    plt.plot(np.arange(1, 9), w2[0][0] * np.arange(1, 9), c='g',label = 'line_cal_w')
    plt.plot(np.arange(1, 9), w3[0][0] * np.arange(1, 9), c='y',label = 'Ridge_cal_w')
    plt.legend()
    plt.show()

结果如下图所示：

注意，调节梯度下降的epochs，lr 以及岭回归的惩罚系数可以使三条直线完全吻合，但为了方便展示这里做了区分。

在实验中发现梯度下降法当学习率过大或者过小时会发生梯度消失或是梯度爆炸，原因可能是数据量太小。

同时，线性回归也可以尝试其他更复杂的点以及波士顿房价预测数据集，代码中已给出数据集，请自行尝试。

四、结语

参考内容:

• B站视频：https://www.bilibili.com/video/BV1aE411o7qd?from=search&seid=1227152557997491474&spm_id_from=333.337.0.0
• 博客：https://www.cnblogs.com/oten/p/13953220.html

如有问题欢迎评论区指出，有问题也可以联系我：
1759412770@qq.com