线性回归--用梯度下降算法求解线性回归线性回归是回归算法里最简单的一种算法了。算法通俗介绍：举个二维平面的例子吧。初中

算法简介

线性回归是回归算法里最简单的一种算法了。算法通俗介绍：举个二维平面的例子吧。初中数学里，已知2点可以求一条直线。那么如果给定100个点呢。怎么去求一条直线？一条直线是不可能拟合100个点的。误差一定是会有的，那么线性规划的目的就是用一条直线去拟合这100个点。使得总的误差最小。回归通常用来研究目标和特征之间的线性关系。主要适用于预测分析。

算法场景

我们用代码来构造一下算法的使用场景，首先我们定义n个点。看看这n个点的分布情况。

from sklearn.datasets import load_diabetes
import matplotlib.pyplot as plt
import numpy as np

X = np.arange(0,10,0.05).reshape(-1, 1)
np.random.seed(100)
y = 0.5 * np.arange(0,10,0.05) + 3 + np.random.rand(len(X))

plt.scatter(X.reshape(1, -1), y)

结果如下图所示：

可以看到我们刻意构造的这n个点，大概是围绕 $y = \theta * x + b$ 这根直线分布的。

那么现在有一个问题来了，如果给定一个新的点，横坐标是 $x_0$ ，请预测纵坐标 $y_0$ 。

算法模型

上一步我们已经知道了给定了一批点，按照一定的线性规律分布，我们需要用一条直线去拟合这批样本。

首先假设样本是围绕着

\hat y = \theta · x + b

这条线分布的。 $\hat y$ 表示预测值。那么求解线性回归模型的问题，就转化为求解参数 $\theta$ 和 $b$ 。

上一步已经说过，线性回归的目标就是找到这样一条直线，使得所有点的误差最小。为了求解最佳参数，我们采用什么标准来表示“误差”呢？

我们现在知道某个样本的特征是 $x_i$ 。同时样本 $x_i$ 对应的目标值 $y_i$ 也是已知的。

我们用参数 $\theta$ 和 $b$ 计算得到预测值 $\hat y_i$ :

\hat y_i = \theta· x + b

然后定义单个点的损失

Cost_i =(\hat y_i - y_i)^2

也就是真实值和预测值之间的平方距离。那么n个点的误差之和如下：

L = \frac {1} {n} \sum _{i=1} ^{n} {(\hat y_i - y_i)^2}

把参数代入，得到：

L(\theta,b) = \frac {1} {n} \sum _{i=1} ^{n} (\theta ·x_i + b - y_i)^2

这就是我们要定义的损失函数。那么现在问题已经转化为

(\theta^*, b^*) = arg \min _{(\theta, b)} \sum _{i=1} ^{n} (\theta · x_i + b - y_i)^2

也就是求解 $\theta$ 和 ${b} 使得损失函数具有最小值。

模型求解

线性回归模型有2种常见的求解方式：（1）最小二乘法（2）梯度下降算法。

最小二乘法求解简单，思路是计算 $L 对 \theta 和 b$ 的到导数，使得式子为0。这时候式子具有最小值。详细细节暂不讨论。

由于梯度下降算法具有更广泛的使用场景，我们这里也采用梯度下降算法求解。

梯度下降算法求解线性回归的过程：

（1）初始化 $\theta$ 和 $b$ 为不等于0的很小的值。

（2）用当前的 $\theta$ 和 $b$ 求解预测值 $\hat y$ 。

（3）如果预测值 $\hat y$ 和实际值预测值 $y$ 的误差大于 $\epsilon$ ，跳到（4）；如果预测值 $\hat y$ 和实际值预测值 $y$ 的误差小于 $\epsilon$ ，返回当前的 $\theta$ 和 $b$ 。求解结束。（4）求解损失函数的偏导数，更新参数 $\theta$ 和 $b$ ：

\theta \leftarrow \theta - \eta \frac {\partial L} {\partial \theta}

b \leftarrow b - \eta \frac {\partial L} {\partial b}

自己实现梯度下降算法

class LinearRegression:
    def __init__(self, eta=0.02, n_iters = 1e4, epsilon = 1e-3):
        
        self._theta = None
        self._eta = eta
        self._n_iters = n_iters
        self._epsilon = epsilon
        
    def Cost(self, theta, X_b, y):
        try:
            return np.sum(np.square(y - X_b.dot(theta))) / len(X_b)
        except:
            return float('inf')

    def dCost(self, theta, X_b, y):
        return (X_b.dot(theta) - y).dot(X_b) * 2 / len(X_b)

    def gradient_descent(self, X_b, y, initial_theta, eta, n_iters, epsilon):
        theta = initial_theta
        eta = self._eta
        n_iters = self._n_iters
        epsilon = self._epsilon
        
        i_iter = 0
        try:
            while i_iter < n_iters:
                gradient = self.dCost(theta, X_b, y)
                last_theta = theta
                theta = theta - eta * gradient

                diff = abs(self.Cost(theta, X_b, y) - self.Cost(last_theta, X_b, y))
                if (diff > 1e100):
                    print('eta is too large!')
                    break
                if (diff  < epsilon):
                    break
                i_iter += 1

            self._theta = theta
        except:
            self._theta = float('inf')
        return self

    def fit(self, X_train, y_train):
        assert X_train.shape[0] == y_train.shape[0], 'train data has wrong dimenssion'
        self._X_train = X_train
        self._y_train = y_train
        
        X_b = np.hstack([np.ones((len(self._X_train), 1)), self._X_train])
        initial_theta = np.zeros(X_b.shape[1])
        eta = 0.01
        
        self.gradient_descent(X_b, y_train, initial_theta, self._eta, self._n_iters, self._epsilon)
        
        return self
    
    def predict(self, X_test):
        X_b = np.hstack([np.ones((len(X_test), 1)), X_test])
        return np.dot(X_b, self._theta)

上述代码保存为独立文件，之后就可以引用并且调用了。注意我们模型中的

X\_train = \left [ \begin{matrix} X\_train[0][0] & X\_train[0][1] & X\_train[0][2] & ... & X\_train[0][m] \\ X\_train[1][0] & X\_train[1][1] & X\_train[1][2] & ... & X\_train[1][m] \\ ... \\ X\_train[n][0] & X\_train[0][1] & X\_train[0][2] & ... & X\_train[n][m] \\ \end{matrix} \right]

对应代码中的

X\_b = \left [ \begin{matrix} 1 & X\_train[0][0] & X\_train[0][1] & X\_train[0][2] & ... & X\_train[0][m] \\ 1 & X\_train[1][0] & X\_train[1][1] & X\_train[1][2] & ... & X\_train[1][m] \\ ... \\ 1 & X\_train[n][0] & X\_train[0][1] & X\_train[0][2] & ... & X\_train[n][m] \\ \end{matrix} \right]

在我们的单特征的模型里,其实就是：

X\_b = \left [\begin{matrix} 1 & X\_train[0]\\ 1 & X\_train[1]\\ ... \\ 1 & X\_train[n]\\ \end{matrix} \right]

\_theta = \left [\begin{matrix} b & \theta\\ \end{matrix} \right]

经过矩阵运算后，可以发现

X\_train ·\theta + b == X\_b ·\_theta

等式左边是模型中的表示，等式右边是代码中的表示，这样做是为了方便矩阵运算。

调用自己写的线性规划模型

linear_regression = LinearRegression(
    epsilon = 1e-6,
    eta = 0.02
)

linear_regression.fit(X_train, y_train)
y_predict = linear_regression.predict(X_test)
params = linear_regression._theta

结果如下：

array([3.45402168, 0.51257892])

也就是说,求解得到：

\theta = 0.51257892

b = 3.45402168

其实和当初构造这一批样本时的分布曲线是很接近的。我们把点和直线画在同一个图里，看看效果吧：

plt.scatter(X.reshape(1, -1), y)

theta = params[1]
b = params[0]

X = np.arange(0,10,0.05).reshape(-1, 1)
y_predict = theta * X + b

plt.scatter(X, y_predict)

结果如下图：

其实已经很好的拟合了这一批样本了。

注意：

eta的取值一定要适当，值太小，结果就收敛太慢了。在给定迭代次数运行结束后，还没有收敛到我们想要的值；如果太大了，不但不会收敛，还会发生爆炸式的发散。

epsilon是我们给定的误差范围，当误差epsilon小于一定值得时候，我们认为直线已经满足拟合的准确度了。

调用sklearn中的线性规划模型

接下来我们再调用sklearn中的模型看看效果：

from sklearn.linear_model import LinearRegression
linear_regression = LinearRegression()

linear_regression.fit(X_train, y_train)
y_predict = linear_regression.predict(X_test)
linear_regression.coef_, linear_regression.intercept_

结果如下：

(array([0.51042252]), 3.4682654875517605)

也就是说,通过sklearn中的线性模型，求解得到：

\theta = 0.51042252

b = 3.4682654875517605

看个拟合图：

发现结果是肉眼不可见的差别。

调用sklearn中的线性规划模型时，我们并没有设定具体参数，使用了默认参数。发现通过我们自己模型得到的结果，和sklearn中的模型差别还是很小的，有一定误差，是因为两个模型的可选参数取值不一致。如果各种参数取值一样，相信结果也会更加接近。

总结

线性规划是比较容易理解的算法。本文只介绍了梯度下降算法，对最小二乘法感兴趣的同学，可以自行研究。