机器学习算法笔记(九)：多元线性回归中的梯度下降法

上一篇文章初步讨论了梯度下降法，本篇文章将探讨梯度下降法在线性回归中的应用。

一、多元线性回归中使用梯度下降法的基本思路

上一篇文章中我们提到，梯度下降法就是在损失函数中取一个合适的点作为起始点，不断的递进改变参数值，得到目标函数的极值的过程。多元线性回归也是类似的，但因为多元线性回归中有多个参数，我们不能再以“导数”来衡量函数变化的方向与快慢，而是要以“梯度”为衡量指标。

我们下面来求损失函数 J 的梯度。线性回归的目标是使(y(i) –(i))2尽可能小。注意到，那么目标就变成了使尽可能小。其实就是我们的损失函数 J 。对其各个变量求偏导数并合并成梯度向量，就有：

其中Xb就是较原数据集X左边新增了一列“1”的矩阵

我们可以看到最右边求导的结果，发现求导结果和m有关，m越大损失函数 J 的梯度也就越大，这其实是不合理的。我们希望我们求出的梯度应该是与m无关，为此，我们将整个梯度值除以m，这样求导结果就与m无关了。这时我们的目标函数恰好就是y(i)和(i)的MSE。

有时取是因为分母上的2正好可以与分子上提出的2约去，在实际使用是差别不大，可酌情选用。

二、初步实现多元线性回归中的梯度下降法

新建一个工程并新建一个main.py文件，实现如下代码：

import numpy as np

#生成测试数据
np.random.seed(666)
x = 2 * np.random.random(size=100)
y = x * 3. + 4. + np.random.normal(size=100)

X = x.reshape(-1, 1) #将测试数据的x转化成矩阵

#实现函数J和dJ
def J(theta, X_b, y):
    try:
        return np.sum((y - X_b.dot(theta))**2) / len(X_b)
    except:
        return float('inf')

def dJ(theta, X_b, y):
    res = np.empty(len(theta))
    res[0] = np.sum(X_b.dot(theta) - y)
    for i in range(1, len(theta)):
        res[i] = (X_b.dot(theta) - y).dot(X_b[:,i])
    return res * 2 / len(X_b)

#实现梯度下降法
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), 1)), x.reshape(-1,1)])
initial_theta = np.zeros(X_b.shape[1])
eta = 0.01

theta = gradient_descent(X_b, y, initial_theta, eta)
print(theta)

运行结果如下：

我们当然还是可以将上面的代码封装起来。在上面的MyML包的LinearRegression.py的LinearRegression类中新增一个fit_gd函数，实现如下：

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):
            res = np.empty(len(theta))
            res[0] = np.sum(X_b.dot(theta) - y)
            for i in range(1, len(theta)):
                res[i] = (X_b.dot(theta) - y).dot(X_b[:, i])
            return res * 2 / len(X_b)

        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

修改main.py文件，使用封装好的代码：

import numpy as np
from MyML.LinearRegression import LinearRegression

#生成测试数据
np.random.seed(666)
x = 2 * np.random.random(size=100)
y = x * 3. + 4. + np.random.normal(size=100)

X = x.reshape(-1, 1) #将测试数据的x转化成矩阵

lin_reg = LinearRegression()
lin_reg.fit_gd(X, y)

print(lin_reg.coef_)
print(lin_reg.intercept_)

运行结果应该与上面的代码完全一致。

三、求梯度过程的向量化处理

注意到梯度的结果向量，能否将里面的分量进行向量化的处理，以提高我们程序的执行效率呢？

我们不难发现，梯度的结果中均是向量点乘的形式，并且前半部分相同，只有后半部分不同。联想到矩阵的乘法法则，我们很容易就可以把这个向量拆分成两个矩阵的乘积。

向量化可以大大提升算法运行效率

分解后右边的矩阵（图示其实在下面）就是我们前面提到的Xb，整个结果可以化简成。这个结果是一个行向量，我们按照规范把它化作列向量，根据矩阵转置的运算法则，最终结果就是如下的形式：

虽然在numpy中不区分行向量和列向量，按照这个结果编写程序也完全不会报错，但是为了规范，我们还是化成了列向量。

有了向量化后的结果，我们把它体现在我们的代码上。修改LinearRegression类的dJ函数如下：

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

经过向量化后，我们就可以更快的运行我们的算法了。

四、运用梯度下降法优化后的线性回归预测波士顿房产数据

我们现在开始使用真实的波士顿房产数据来验证我们改进后的线性回归法。新建一个工程并创建一个main.py文件，实现如下：

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]

from MyML.model_selection import train_test_split

X_train, X_test, y_train, y_test = train_test_split(X, y, seed=666)

from sklearn.preprocessing import StandardScaler

#使用梯度下降法前进行数据归一化，否则会因为数据不在同一尺度导致不收敛
standardScaler = StandardScaler()
standardScaler.fit(X_train)
X_train_standard = standardScaler.transform(X_train)

from MyML.LinearRegression import LinearRegression
lin_reg = LinearRegression()


lin_reg.fit_gd(X_train_standard, y_train)

X_test_standard = standardScaler.transform(X_test)
print(lin_reg.score(X_test_standard, y_test))

输出结果如下：

我们发现这个预测准确度与我们先前运用正规方程解求出的准确度完全一致，说明我们成功通过梯度下降找到了损失函数的极小值。同时，相比运用正规方程解，梯度下降法的效率会更高一些。

由于有 η 这个变量，如果最终数据集的数据不在一个维度上（大小相差很大）将会影响梯度的结果。梯度的结果乘以 η 是损失函数真正每次变化的步长，则步长有可能或者太大或者太小。如果步长太大会导致损结果不收敛，如果步长太小会导致搜索过程太慢或者没有找到极小值就因为循环上限而退出循环（欠拟合）。所以在运用梯度下降法前，对数据进行归一化处理显得尤为重要。

打赏

发表回复

您的邮箱地址不会被公开。 必填项已用 * 标注