用python实现线性回归

基本形式

线性回归基本形式：

其中$ \theta $代表参数， $ \theta_0 $就是通常说的偏置bias。$ x_n $代表特征，如果$ n = 1 $，则是一元线性回归；$ n > 1 $时是多元线性回归。

在实际实现中，常常采取如下的向量化:

其中$ X $代表的是训练数据，并且每个训练数据是一行。假设每条训练数据有n个特征（包含$x_0$），总共有m条训练数据，那么$ X $就是$m \times n$的矩阵。$ \theta $代表的是参数矩阵，因为训练数据有n个特征，所以$ \theta $是$n \times 1$的列向量。

模型求解

根据以上的公式，模型求解的目标就是要求出参数$\theta$。最常见的有两种方法，一种是梯度下降法，另一种是正规方程（normal equation)。

梯度下降法

介绍梯度下降算法之前需要先介绍一下线性回归模型的损失函数，常用的平方损失函数，定义如下：

一个直观的解释是，线性回归的损失函数计算的是预测点到实际点的距离的平方和。至于为什么要除以2，是因为在计算损失函数偏导的时候，会乘以2，所以刚好能够抵消，最后系数只剩下 $\frac{1}{m}$。

损失函数的向量化形式如下：

向量化的意义在于用向量的形式实现，就可以使用矩阵来进行运算，相比于普通的循环来说更加高效，特别是当数据量很大的时候。

在知道损失函数之后，我们的目标就转换为，求解参数$\theta$使得损失函数的值达到最小。梯度下降法就是对损失函数求每个参数$\theta_j$的偏导，偏导就称为梯度，然后通过将当前的$\theta_j$减去其梯度的这种方式来更新参数，直至算法收敛（算法收敛的条件可以是损失函数的值小于某个值）。公式如下：

其中$\alpha$是学习率，$\frac{1}{m} \sum\limits_{i=1}^{m} (h_\theta(x^{(i)}) - y^{(i)}) \cdot x_j^{(i)}$就是$\theta_j$对应的梯度。所以梯度下降法是一个迭代的算法，每轮迭代同时更新所有参数；另外还需要注意学习率$\alpha$的选择，$\alpha$太大的话，可能导致算法不收敛，$\alpha$太小可能导致迭代次数过多，算法收敛慢。

下面是梯度下降法更新参数的向量化形式：

正规方程

梯度下降法是个迭代的算法，需要迭代很多次，而正规方程则可以直接计算出线性回归的参数。这里先举个直观的例子，求某个二次函数的极值，例如$y = x^2 - 2x + 1$，要求它的极小值，可以对$x$求导，得到$2x - 2$，令$2x - 2 = 0$，得到$x = 1$是极小值点。正规方程也是基于这个思想，回顾前文所提到的线性回归模型损失函数的向量化形式，正规方程是根据损失函数对参数矩阵$\theta$求导（这里涉及到矩阵求导，不再介绍，具体推导过程可以参考这里），令求导式子等于0，直接得出参数矩阵$\theta$的值。正规方程如下：

正规方程相比于梯度下降法，不需要迭代，也不需要选择学习率$\alpha$，只需要对矩阵求逆。当特征数量很少时，可以采用normal equation的方式，更加直接；当特征数量很大时，矩阵求逆会耗费大量的时间，采用梯度下降法更佳。另外，使用normal equation时，可能出现矩阵$X^T X$不可逆的情况，可能的原因是存在多余的特征，特征与特征之间依赖比较严重，例如线性相关；或者特征数太多，远大于样本数量，此时应当删去不必要的特征或者采取正则化（Regularization）。

具体实现

线性回归算法的核心就是参数的求解，这里给出梯度下降算法的实现和正规方程的python实现。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59

import numpy

class LinearRegression:
    def __init__(self, n_iterations=1000, learning_rate=0.01, gradient_descent=True):
        """初始化.
        Args:
            n_iterations: 梯度下降算法的迭代次数.
            learning_rate: 梯度下降算法的学习速率.
            gradient_descent: True代表使用梯度下降算法求解，False代表使用normal equation求解.
        """
        self.n_iterationns = n_iterations
        self.learning_rate = learning_rate
        self.gradient_descent = gradient_descent

    def computeCost(self, X, y):
        """计算当前模型对于训练数据的损失(平方损失)
        Args:
            X: 训练数据中的特征
            y: 训练数据中的目标值

        Returns:
            当前模型对于训练数据的损失
        """
        n_samples = len(y)
        return (X.dot(self.theta) - y).T.dot(X.dot(self.theta) - y) / 2 / n_samples

    def fit(self, X, y):
        """用训练数据来训练模型
        Args:
            X: 训练数据中的特征
            y: 训练数据中的目标值

        Returns:
            theta: 线性回归模型的参数
            costHistory: 如果使用梯度下降算法求解，梯度下降算法过程中模型的历史损失会被返回
        """

        # 加入常量1作为偏置bias特征
        X = np.insert(X, 0, 1, axis=1)
        y = np.array(y).reshape(len(y), 1)
        n_samples, n_features = X.shape
        self.theta = np.zeros(n_features).reshape(n_features, 1)

        costHistory = np.zeros(self.n_iterationns)

        if not self.gradient_descent:
            # normal equation
            self.theta = np.linalg.pinv(X.T.dot(X)).dot(X.T).dot(y)
        else:
            for i in range(self.n_iterationns):
                tmp = X.T.dot(X.dot(self.theta) - y)
                self.theta -= tmp * self.learning_rate / n_samples
                costHistory[i] = self.computeCost(X, y)

        return self.theta, costHistory

    def predict(self, X):
        X = np.insert(X, 0, 1, axis=1)
        return X.dot(self.theta)

在实现梯度算法的过程中，我们可以计算每一步当前模型的损失，用来判断梯度下降法是否正常工作，因此在上面的fit方法实现时，返回了梯度下降算法过程中模型的历史损失，可以通过绘图的方式来观察其变化，如下图所示：

可以看到模型的损失随着迭代次数的增加而不断减少，证明梯度下降算法是正常工作的。

下面使用sklearn自带的数据diabetes来测试我们的线性回归模型，用scatterplot画出数据点的分布范围，用不同的直线画出不同迭代次数时的模型拟合情况以及使用正规方程时的模型拟合情况，如下图：

可以看到随着迭代次数不断增加，拟合的效果越来越好，越来越接近正规方程的拟合效果。图中的模型拟合时学习速率设置相同，都为0.05，读者还可以尝试不同学习速率导致的拟合效果变化，完整代码在这里。

多项式回归

这里有必要提一下多项式回归，其实多项式回归的本质也是线性回归，只不过添加了高维的特征。要实现多项式回归，上面的线性回归代码都无需改动，只需对特征$X$进行改动，添加高阶特征即可，感兴趣的读者可以自行实现并寻找一些训练数据来查看拟合效果。

总结

本文用python实现了简单的线性回归模型，实现了梯度下降算法和正规方程（normal equation）。由于梯度下降的过程中我们考虑了所有的样本点，这种梯度下降方式也称为批量梯度下降，在训练数据量较小时可以使用，但训练数据很大时会导致训练时间太长，不可取，此时可以使用随机梯度下降算法或者小批量梯度下降算法。

本文还没有涉及到模型的正则化（Regularization），用于解决模型的过拟合问题，会在以后实现其他模型的时候涉及，敬请期待！