线性回归的代码实现

1. 线性回归推导

​ 对于给定了$m$个样本点($x^{(1)}$,$y^{(1)}$)，($x^{(2)}$,$y^{(2)}$)，$\dots$，($x^{(m)}$,$y^{(m)}$)，其中$x^{(i)}\in\R^{n}$。定义假设函数为$h_\theta(x)$，即$h_\theta(x)$为最终的拟合函数，$\theta$为待拟合参数也称作权重。
$$
\begin{aligned}
h_\theta(x)&=\theta_0+\theta_1x_1+\cdots+\theta_nx_n \
&=\begin{bmatrix}
\theta_0, & \theta_1,& \cdots ,& \theta_n
\end{bmatrix}
\begin{bmatrix} 1\x_1 \ \vdots \ x_n\end{bmatrix}\
&=\theta^Tx

\end{aligned}
$$

1.1 从样本数据出发推导损失函数

​ 在样本数据中$y^{(i)}$是实际存在值而$h_\theta(x^{(i)})$对应的是模型预测值，显然如果想要模型预测的效果好，那么对应的误差就要小，假设函数在任意样本点的误差为$|h_\theta(x^{(i)})-y^{(i)}|$则$m$个样本点的误差和为$\ \sum\limits_{i=1}^m|h_\theta(x^{(i)})-y^{(i)}|$，因此问题就转化为求解$\ \begin{aligned}\mathop{\arg\min}{\theta}\ \sum\limits{i=1}^m|h_\theta(x^{(i)}-y^{(i)}|\end{aligned}$为了后续求解最优值(绝对值函数不好求导)，所以损失函数采用了误差平方和的形式$\ \begin{aligned}\mathop{\arg\min}{\theta} \sum\limits{i=1}^m(h_\theta(x^{(i)}-y^{(i)})^2\end{aligned}$。

1.2 从统计学理论出发推导损失函数

​ 为什么线性回归问题的损失函数会是误差平方和的形式？这里从统计理论上进行解释，对于给定的$y^{(i)}$总能找到$\varepsilon^{(i)}$使得这个等式成立$y^{(i)}=h_\theta(x^{(i)})+\varepsilon^{(i))}$，$\ \varepsilon^{(i)}$代表真实值和预测值之间的误差且$\varepsilon^{(i)} \sim \mathcal{N}(0,\sigma^2)$,简单解释下为什么误差$\varepsilon^{(i)}$会服从均值为零的正态分布，误差的产生有很多种因素的影响，误差可以看作是这些因素(随机变量)之和共同作用而产生的，由中心极限定理可知随机变量和的分布近似的服从正态分布；更通俗易懂的解释是，当你在选择$h_\theta(x)$时主观的会认定这个$h_\theta(x)$是比较符合样本数据的，比如对一些样本数据可视化后，发现样本数据明显是趋近于一条直线，而你在对$h_\theta(x)$的选择上肯定会选择直线方程作为$h_\theta(x)$而不会选择多项式函数作为$h_\theta(x)$。而这种$h_\theta(x)$一旦选定，可以认为大部分数据都在$h_\theta(x)$的附近，因此误差大部分集中在零值附近所以$\mathcal{N}(0,\sigma^2)$作为$\varepsilon^{(i)}$的先验分布是比较合理的。

​ 随机变量$\varepsilon^{(i)}$的概率密度函数为：
$$
p(\varepsilon^{(i)}) = \frac{1}{\sqrt{2\pi}\sigma} e^{-\frac{(\varepsilon^{(i)})^2}{2\sigma^2}}
$$
​ 代入$h_\theta(x^{(i)}), \ y^{(i)}$则
$$
p(y^{(i)}|x^{(i)};\theta) = \frac{1}{\sqrt{2\pi}\sigma} e^{-\frac{(\theta^Tx^{(i)}-y^{(i)})^2}{2\sigma^2}}
$$
这里的$p(y^{(i)}|x^{(i)};\theta)$并不代表条件概率，只是一个记号它表示给定$x^{(i)},\ y^{(i)}\ $和一组$\theta$后的概率密度函数。由最大似然估计可知：
$$
\begin{aligned}
L(\theta)&=\prod\limits_{i=1}^mp(y^{(i)}|x^{(i)};\theta)\
&=\prod\limits_{i=1}^m\frac{1}{\sqrt{2\pi}\sigma} e^{-\frac{(\theta^Tx^{(i)}-y^{(i)})^2}{2\sigma^2}}
\end{aligned}
$$
对$L(\theta)$取对数从而得到对数化最大似然估计函数
$$
\begin{aligned}
\mathcal L(\theta) &= \mathcal log(L(\theta)) \
&=\mathcal \log \prod\limits_{i=1}^m\frac{1}{\sqrt{2\pi}\sigma} e^{-\frac{(\theta^Tx^{(i)}-y^{(i)})^2}{2\sigma^2}}\
&=\sum\limits_{i=1}^m \mathcal \log\frac{1}{\sqrt{2\pi}\sigma}e^{-\frac{(\theta^Tx^{(i)}-y^{(i)})^2}{2\sigma^2}}\
&=m \mathcal \log \frac{1}{\sqrt{2\pi}\sigma} - \frac{1}{2\sigma^2}\sum\limits_{i=1}^m(\theta^Tx^{(i)}-y^{(i)})^2
\end{aligned}
$$
模型的损失函数$\ J(\theta)=\dfrac{1}{2}\sum\limits_{i=1}^m(h_\theta(x^{(i)})-y^{(i)})^2$，其实只是在先前的最优化模型前面乘了一个常系数$\dfrac{1}{2}$对最优化的结果并不会产生影响。

2. 梯度下降算法求解凸优化问题

​ 在机器学习算法中，求解最小化损失函数时，可以通过梯度下降法来一步步进行迭代求解，从而得到最小化损失函数的模型参数值，梯度下降算法不一定能够找到全局的最优解，有可能是一个局部最优解。然而，如果损失函数是凸函数，那么梯度下降法得到的解就一定是全局最优解。 在这里不加证明的给出梯度下降法的迭代格式:
$$
\theta_j:=\theta_j-\alpha\nabla_{\theta_j} {J(\theta)}
$$
其中$\theta_j$为假设函数$h_\theta(x)$的参数值，$\alpha \in (0,1]$且为常数代表模型学习速率，$\nabla _{\theta_j} J(\theta)$为损失函数$J(\theta)$对参数$\theta_j$的梯度。

​ 下面使用梯度下降算法求解上面推导出的线性回归模型的损失函数$\ J(\theta)=\dfrac{1}{2}\sum\limits_{i=1}^m(h_\theta(x^{(i)})-y^{(i)})^2$，为了实现该算法首先要求出损失函数$J(\theta)$对参数$\theta_j$的梯度：
$$
\begin{aligned}
\nabla_{\theta_j} J(\theta) &=\nabla_{\theta_j} \dfrac{1}{2} \sum\limits_{i=1}^m(h_\theta(x^{(i)})-y^{(i)})^2\
&=2\cdot\dfrac{1}{2}\sum\limits_{i=1}^m(h_\theta(x^{(i)})-y^{(i)}) \nabla_{\theta_j}(h_\theta(x^{(i)})-y^{(i)})\
&=\sum\limits_{i=1}^m(h_\theta(x^{(i)})-y^{(i)}) \nabla_{\theta_j}\theta^Tx^{(i)}\
&=\sum\limits_{i=1}^m(h_\theta(x^{(i)})-y^{(i)})x_j^{(i)}
\end{aligned}
$$
因此在线性回归模型中利用所有的样本数据，训练梯度下降算法的完整迭代格式为
$$
\theta_j:=\theta_j-\alpha \sum\limits_{i=1}^m(h_\theta(x^{(i)})-y^{(i)})x_j^{(i)} \ \ \ ( j\ \ for \ \ 0 \ \sim \ n)
$$
上述迭代过程每次迭代都会使用所有的样本数据，数学上已经证明线性回归模型的损失函数通过梯度下降算法求解一定会全局收敛，所以如果要编程实现该算法只需要控制迭代次数即可，不过对于线性回归模型的求解一般不用梯度下降算法，还有更容易实现且更快捷的形式—正规方程。

3. 正规方程

多元线性回归的损失函数为：
$$
J=\sum_{i=1}^{m}(y^{(i)} - \hat{y}^{(i)})^2
$$
其中：$\hat{y}^{(i)} = \theta_{0} + \theta_{1}X_{1}^{(i)} + \theta_{2}X_{2}^{(i)} + … + \theta_{n}X_{n}^{(i)}$ 。

对 $J$ 求导为：
$$
\nabla J=(\frac{\partial J}{\partial \theta_0},\frac{\partial J}{\partial \theta_1},…,\frac{\partial J}{\partial \theta_n})
$$
其中：$\frac{\partial J}{\partial \theta_i}$ 为偏导数，与导数的求法一样。

对 $\nabla J$ 进一步计算：
$$
\nabla J(\theta) = \begin{pmatrix} \frac{\partial J}{\partial \theta_0} \\ \frac{\partial J}{\partial \theta_1} \\ \frac{\partial J}{\partial \theta_2} \\ \cdots \\ \frac{\partial J}{\partial \theta_n} \end{pmatrix} = \begin{pmatrix} \sum_{i=1}^{m}2(y^{(i)} - X_b^{(i)}\theta)·(-1) \\ \sum_{i=1}^{m}2(y^{(i)} - X_b^{(i)}\theta)·(-X_1^{(i)}) \\ \sum_{i=1}^{m}2(y^{(i)} - X_b^{(i)}\theta)·(-X_2^{(i)}) \\ \cdots \\ \sum_{i=1}^{m}2(y^{(i)} - X_b^{(i)}\theta)·(-X_n^{(i)}) \end{pmatrix} = 2·\begin{pmatrix} \sum_{i=1}^{m}(X_b^{(i)}\theta - y^{(i)}) \\ \sum_{i=1}^{m}(X_b^{(i)}\theta - y^{(i)})·X_1^{(i)} \\ \sum_{i=1}^{m}(X_b^{(i)}\theta - y^{(i)})·X_2^{(i)} \\ \cdots \\ \sum_{i=1}^{m}(X_b^{(i)}\theta - y^{(i)})·X_n^{(i)} \end{pmatrix}
$$

其中：$X_b = \begin{pmatrix}
1 & X_1^{(1)} & X_2^{(1)} & \cdots & X_n^{(1)} \
1 & X_1^{(2)} & X_2^{(2)} & \cdots & X_n^{(2)} \
\vdots & \vdots & \vdots & \cdots & \vdots \
1 & X_1^{(m)} & X_2^{(m)} & \cdots & X_n^{(m)}
\end{pmatrix}$

​ 相应的对$J$上对$θ$这个向量去求梯度值，也就是损失函数$J$对$θ$每一个维度的未知量去求导。此时需要注意，求导过程中，$θ$是未知数，相应的$X$和$y$都是已知的，都是在监督学习中获得的样本信息。对于最右边式子的每一项都是m项的求和，显然梯度的大小和样本数量有关，样本数量越大，求出来的梯度中，每一个元素相应的也就越大，这个其实是不合理的，求出来的梯度中每一个元素的值应该和$m$样本数量是无关的，为此将整个梯度值再除上一个m，相应的目标函数的式子变成了下面的式子即：
$$
\nabla J(\theta) = \frac{2}{m}\begin{pmatrix} \sum_{i=1}^{m}(X_b^{(i)}\theta - y^{(i)}) \\ \sum_{i=1}^{m}(X_b^{(i)}\theta - y^{(i)})·X_1^{(i)} \\ \sum_{i=1}^{m}(X_b^{(i)}\theta - y^{(i)})·X_2^{(i)} \\ \cdots \\ \sum_{i=1}^{m}(X_b^{(i)}\theta - y^{(i)})·X_n^{(i)} \end{pmatrix}
$$
​ 此时，目标函数就成了使 $\frac{1}{m}\sum_{i=1}^{m}(y^{(i)} - \hat{y}^{(i)})^2$ 尽可能小，即均方误差尽可能小：
$$
J(\theta) = MSE(y, \hat{y})
$$
​

注1. 有时候目标函数也去 $\frac{1}{\boldsymbol{2}m}\sum_{i=1}^{m}(y^{(i)} - \hat{y}^{(i)})^2$ ,其中的2是为了抵消梯度中的2，实际的效果有限。

注2. 将梯度除以m相当于目标函数本身变成了MSE，也就是对原来的目标函数除上m。如果没有1/m的话，梯度中的元素就会特别的大，在具体编程实践中就会出现一些问题。当我们在使用梯度下降法来求函数的最小值的时候，有时候需要对目标函数进行一些特殊的设计，不见得所有的目标函数都非常的合适，虽然理论上梯度中每一个元素都非常大的话，我们依然可以通过调节learning_rate(学习率)来得到我们想要的结果，但是那样的话可能会影响效率。

2.实现代码

2.1 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
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
# -*- coding: utf-8 -*-
"""
Created on Wed Jan  8 17:25:39 2020

@author: lixin
"""

import numpy as np

class LinearRegression():
    def __init__(self,lr=.01, num_iters=10000,tolerance=1e-8):
        self.lr = lr
        self.num_iters = num_iters
        self.w = None
        self.toterance = tolerance

    def fit(self,x_train, y_train):
        n_samples = len(x_train)
        x_train = np.column_stack((np.ones(n_samples), x_train)) #也可以写成np.c_
        
        self.w = .01 * np.ones(x_train.shape[1])
        self.loss_ = [0]
        
        # w_{i} = w_{i} - lr * (h(w_{i}) - y)*x_{i} 迭代公式
        self.count = 0
        for iteration in range(self.num_iters):
            self.count += 1
            raw_output = np.matmul(x_train, self.w)
            errors =  raw_output - y_train
            loss = 1/(2 * n_samples) * errors.dot(errors)
            delta_loss = loss - self.loss_[-1]

            self.loss_.append(loss)
            if np.abs(delta_loss) < self.toterance:
                break
            else:
                grad = (1.0 /n_samples) *np.matmul(x_train.T, np.array(errors))
                self.w -= self.lr * grad

    def predict(self,x_test):
        x_test = np.column_stack((np.ones(len(x_test)), x_test))
        
        output = np.matmul(x_test, self.w)
        return output
    
if __name__ == '__main__':
    
    num_inputs = 2
    num_examples = 10000
    true_w = [6.4, -3.2]
    true_b = 2.3
    features = np.random.random((num_examples, num_inputs))
    labels = true_w[0] * features[:, 0] + true_w[1] * features[:, 1] + true_b
    labels += np.random.normal(0, 0.1,size = len(labels))
    
    import sklearn.model_selection
    x_train, x_test, y_train, y_test = sklearn.model_selection.train_test_split(features, labels, test_size = .20, random_state=42)

    lr = LinearRegression()
    lr.fit(x_train,y_train)
    
    print("回归系数:",lr.w)
    print("迭代次数:",lr.count)
    
    y_pred = lr.predict(x_test)
    from sklearn import metrics
    mse = metrics.mean_squared_error(y_test, y_pred)
    print("MSE: %.4f" % mse)

    mae = metrics.mean_absolute_error(y_test, y_pred)
    print("MAE: %.4f" % mae)

    R2 = metrics.r2_score(y_test,y_pred)
    print("R2: %.4f" % R2)
2.2 Scala版本实现
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 breeze.linalg.{DenseMatrix => BDM, DenseVector => BDV}
import scala.collection.mutable.ArrayBuffer

/**
 * Scala 版本的实现
 */

object LinearRegression {
  def main(args: Array[String]): Unit = {
    val num_inputs = 2
    val num_examples = 1000
    val x_train: BDM[Double] = BDM.rand(num_examples, num_inputs)
    val ones = BDM.ones[Double](num_examples, 1)
    val x_cat = BDM.horzcat(ones, x_train)
    val y_train = x_cat * BDV(2.3, 6.4, -3.2)

    val model = new LinearRegression(num_iters = 10000)
    val weights = model.fit(x_train, y_train)
    val predictions = model.predict(weights, x_train)
    println("梯度下降求解的权重为：" + weights)
    println(predictions)
  }
}

class LinearRegression(var lr: Double = 0.01, var tolerance: Double = 1e-6, var num_iters: Int = 1000) {

  def fit(x: BDM[Double], y_train: BDV[Double]): BDV[Double] = {
    val ones = BDM.ones[Double](x.rows, 1)
    val x_train = BDM.horzcat(ones, x)
    val n_samples = x_train.rows
    val n_features = x_train.cols
    var weights = BDV.ones[Double](n_features) :* .01 // 注意是:*

    val loss_lst: ArrayBuffer[Double] = new ArrayBuffer[Double]()
    loss_lst.append(0.0)

    var flag = true
    for (i <- 0 to num_iters if flag) {
      val raw_output = x_train * weights
      val error = raw_output - y_train
      val loss: Double = error.t * error
      val delta_loss = loss - loss_lst.apply(loss_lst.size - 1)
      loss_lst.append(loss)
      if (scala.math.abs(delta_loss) < tolerance) {
        flag = false
      } else {
        val gradient = (error.t * x_train) :/ n_samples.toDouble
        weights = weights - (gradient :* lr).t
      }
    }
    weights
  }

  def predict(weights: BDV[Double], x: BDM[Double]): BDV[Double] = {
    val x_test = BDM.horzcat(BDM.ones[Double](x.rows, 1), x)
    val output = x_test * weights
    output
  }
}

参考文献：1. 机器学习之梯度下降法与线性回归