线性回归梯度下降法python_线性回归中的梯度下降

0x00 前言

在上一篇文章《手动实现梯度下降(可视化)》的多元线性回归中，使用梯度下降法求损失函数的最小值。在多元线性回归中，我们已经使用“标准方程解”求了损失函数最小值，那么如何使用梯度下降来计算最小值呢？在本篇文章中，对多元线性回归的损失函数进行求导，并使用向量化的方式编写代码。最后证明了，在梯度下降之前需要使用归一化。

0x01 理论推导

1.1 目标

在多元线性回归中，其预测值为： 

已知训练数据样本、 ，找到，使得损失函数尽可能小.：

下面就要使用梯度下降法求多元线性回归中损失函数的最小值。

1.2 整理损失函数

使用梯度下降法来求损失函数最小值时，要对损失函数进行一定的设计。

下面我们要进行一下化简。

将预测值  代入到损失函数中，可以得到：

其中  是列向量列向量，而且我们注意到，可以虚构第0个特征，另其恒等于1，推导时结构更整齐，也更加方便：

这样我们就可以改写成向量点乘的形式：

将带入损失函数中，可以将其改写为：

对进行求导，可以得到:

在上式中梯度大小实际上是和样本数量m相关，m越大，累加之和越大，这是不合理的；为了使与m无关，除以m。

于是演变成：使损失函数尽可能小。同时注意到，损失函数变成了均方误差MSE：。

对其求导，有：

我们将这个复杂的表达式再一次进行分解，可以写成向量  再去点乘矩阵：

即：将求梯度的过程转换为两个矩阵之间的乘法运算：1*m的向量 乘以 m*(n+1)的矩阵，得到1*(n+1)的向量。其结果是 。

然后这里还需要一个转秩。为什么转秩？本来应该是一个列向量，但是这里需要转秩成行向量。因为得到的是1*(n+1)的行向量，梯度是(n+1)*1的列向量，因此需要再次进行转秩(分别转秩，再交换位置)。

1.3 最终得到的梯度

最终得到梯度结果为： 

0x02 代码演示

2.1 梯度下降代码

下面就可以整合梯度下降的代码，将上面推导出来的梯度结果的表达式带入梯度下降的流程中。

def fit_gd(self, X_train, y_train, eta=0.01, n_iters=1e4):    """根据训练数据集X_train, y_train, 使用梯度下降法训练Linear Regression模型"""    assert X_train.shape[0] == y_train.shape[0], \        "the size of X_train must be equal to the size of y_train"    def J(theta, X_b, y):        try:            return np.sum((y - X_b.dot(theta)) ** 2) / len(y)        except:            return float('inf')            def dJ(theta, X_b, y):        return X_b.T.dot(X_b.dot(theta) - y) * 2. / len(y)    def gradient_descent(X_b, y, initial_theta, eta, n_iters=1e4, epsilon=1e-8):        theta = initial_theta        cur_iter = 0        while cur_iter < n_iters:            gradient = dJ(theta, X_b, y)            last_theta = theta            theta = theta - eta * gradient            if (abs(J(theta, X_b, y) - J(last_theta, X_b, y)) < epsilon):                break            cur_iter += 1        return theta    X_b = np.hstack([np.ones((len(X_train), 1)), X_train])    initial_theta = np.zeros(X_b.shape[1])    self._theta = gradient_descent(X_b, y_train, initial_theta, eta, n_iters)    self.intercept_ = self._theta[0]    self.coef_ = self._theta[1:]    return self

2.2 使用真实数据出现的问题

下面准备波士顿房产数据，准备进行验证

import numpy as npfrom sklearn import datasetsboston = datasets.load_boston()X = boston.datay = boston.targetX = X[y < 50.0]y = y[y < 50.0]

首先使用线性回归中的“标准方程解”进行计算，得到相应系数与截距的最优值，然后计算回归评分(R方)：

from myAlgorithm.LinearRegression import LinearRegressionfrom myAlgorithm.model_selection import train_test_splitX_train, X_test, y_train, y_test = train_test_split(X, y, seed=666)lin_reg1 = LinearRegression()%time lin_reg1.fit_normal(X_train, y_train)lin_reg1.score(X_test, y_test)"""输出：CPU times: user 25.7 ms, sys: 25.5 ms, total: 51.2 msWall time: 116 ms0.8129802602658466"""

然后创建一个lin_reg2，使用梯度下降法得到最优参数：

lin_reg2 = LinearRegression()lin_reg2.fit_gd(X_train, y_train, )lin_reg2.coef_"""输出：array([nan,nan,nan,nan,nan,nan,nan,nan,nan,nan,nan,nan,nan])"""

发现在执行的过程中有一个RuntimeWarning的警告overflow。并且得到的全都是空。

这是因为，在一个真实的数据中，查看前3行数据，就可以观察到，每一个特征所对应的规模是不一样的，有一些很小，有一些很大，使用默认的学习率可能过大，导致不收敛。

下面我们传递一个很小的学习率来看一下结果：

lin_reg2.fit_gd(X_train, y_train, eta=0.000001)lin_reg2.score(X_test, y_test)"""输出：0.27556634853389195"""

我们发现得到的R方评分很差，这可能是因为学习率设置的过小，步长太小没有到最优值。为了验证这个假设，可以再进行一次验证：

lin_reg2.fit_gd(X_train, y_train, eta=0.000001, n_iters=1e6)lin_reg2.score(X_test, y_test)"""输出：0.75418523539807636"""

我们发现，即使计算了这么久，其结果还是没有“标准方程解”的评分高。

对于这种情况应该怎么办呢？之前已经分析出来了，之所以出现这种现象，是因为真实数据整体不在一个规模上，解决的方式，就是在梯度下降之前进行数据归一化。

2.3 数据归一化

from sklearn.preprocessing import StandardScalerstandardScaler = StandardScaler()standardScaler.fit(X_train)X_train_std = standardScaler.transform(X_train)lin_reg3 = LinearRegression()lin_reg3.fit_gd(X_train_std, y_train)X_test_std = standardScaler.transform(X_test)lin_reg2.score(X_test, y_test)"""输出：0.8129802602658466"""

这和“正规方程解”得到的score一样，这就说明我们找到了损失函数的最小值。

0xFF 总结

在本篇文章中，对多元线性回归的损失函数进行求导，其结果是 。

然后使用向量化的方式编写代码，但是发现在真实数据中效果比较差，这是因为数据的规模不一样，因此在梯度下降之前需要使用归一化。