5.3 线性回归中的梯度下降学习笔记

1、代码演示

1.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

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

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

import numpy as np
from sklearn import datasets

boston = datasets.load_boston()
X = boston.data
y = boston.target

X = X[y < 50.0]
y = y[y < 50.0]

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

from LinearRegression import LinearRegression
from model_selection import train_test_split

X_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 ms
Wall time: 116 ms
0.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。并且得到的$\theta$全都是空。

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

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

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
"""

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

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

1.3 数据归一化

from sklearn.preprocessing import StandardScaler

standardScaler = 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一样，这就说明我们找到了损失函数的最小值。

2.总结

在本篇文章中，对多元线性回归的损失函数进行求导，其结果是$\frac{2}{m}(X_b\theta -y{)}^T{X}_b)$ 。

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