机器学习算法笔记(十)：随机梯度下降法

我们上一篇文章讨论的梯度下降法是要将我们想要最优化的那个损失函数相应在某一点 θ 的梯度值准确的求出来。要想求出准确的梯度，每一项都要对所有的样本进行计算，这样的梯度下降法又被称为批量梯度下降法(Batch Gradient Descent)，每一次的计算过程都要将样本中所有的信息批量进行计算。但若我们的样本量非常大，我们的算法在计算梯度时就会非常耗时。为了解决样本量过大时的运行效率问题，我们就提出了随机梯度下降法(Stochastic Batch Gradient Descent)。

一、随机梯度下降法的基本思路

上一篇文章批量梯度下降法的思路是计算损失函数的梯度，按梯度的方向，逐步减小损失函数的变量 θ，对应的损失函数也不断减小，直到损失函数的的变化量满足精度要求。

而随机梯度下降法，是每次只对其中一个样本进行导数的计算，每一次只取出X_b中随机的一个 i （也就是随机的一行）。

如上图所示，我们随机抽取 n 个样本，在每个样本的梯度方向上逐步优化变量 θ（每随机抽取一个样本就对 θ 做一次递减优化）。批量梯度下降法的优化，是整体数据集的梯度方向逐步循环递减变量 θ；而随机梯度下降法，是数据集中的一个随机的样本的梯度方向，优化变量 θ。

与批量梯度下降法直接从整体梯度的方向搜索不同，随机梯度下降法每一次搜索的方向，不能保证是损失函数减小的方向（也就是说可能还会往损失函数增大的方向搜索），也不能保证是损失函数减小最快的方向，它的优化方向具有不可预知性。

实验结论表明，即使随机梯度下降法的优化方向具有不可预知性，通过此方法依然可以差不多来到损失函数最小值的附近。虽然不像批量梯度下降法那样，一定可以来到损失函数最小值位置，但是，如果样本数量很大时，有时可以用一定的模型精度，换取优化模型所用的时间。

二、随机梯度下降法中 η 的确定过程

在随机梯度下降法中，确定学习率 η 的取值非常重要。在优化损失函数的过程中，如果 η 一直取固定值，可能会出现已经优化到损失函数最小值位置了，但由于随机的过程不够好，η 又是固定值，导致优化时慢慢的又跳出最小值位置。一个较好的解决方案是让优化过程中让 η 逐渐递减（随着梯度下降法循环次数的增加，η 值越来越小）。

确定 η 的过程大致可以分成以下三步：

1、让 η 逐渐递减，我们很容易想到倒数的形式。如果 η = 1 / i_iters，（i_iters：当前循环次数）就有随着循环次数（i_iters）的增加，η 的变化率差别太大的问题。若循环次数从1到2，我们的 η 一下就下降了50%；若循环次数从10000到10001，η 值才下降了0.01%；
2、所以我们在分母上添加一个常数b来缓解上面的问题，即 η = 1 / (i_iters + b)（b为常量）；
3、再次变形：η = a / (i_iters + b)，（a、b为常量）分子改为 a ，增加 η 取值的灵活度。

● a和 b实际上是随机梯度下降法的两个超参数。一般在经验上我们 a 取5，b 取50。

● 学习率是随着循环次数的增加逐渐递减的，这也是搜索领域一个重要思想——模拟退火思想。在这个思想中，我们将 η 的值用 t0 和 t1 表示 为。

三、编程实现随机梯度下降法

新建一个工程，创建一个main.py文件，输入以下代码：

import numpy as np

#生成测试数据
m = 100000
x = np.random.normal(size=m)
X = x.reshape(-1,1)
y = 4.*x + 3. + np.random.normal(0, 3, size=m)

def dJ_sgd(theta, X_b_i, y_i):
    """
    :param theta: 传入θ值
    :param X_b_i: 传入随机抽取的单行向量
    :param y_i: 传入单行向量对应的标记
    """
    return 2 * X_b_i.T.dot(X_b_i.dot(theta) - y_i)

#随机梯度下降法函数实现
def sgd(X_b, y, initial_theta, n_iters):

    t0, t1 = 5, 50
    def learning_rate(t): #学习率的计算
        return t0 / (t + t1)

    theta = initial_theta
    for cur_iter in range(n_iters):
        rand_i = np.random.randint(len(X_b))
        gradient = dJ_sgd(theta, X_b[rand_i], y[rand_i])
        theta = theta - learning_rate(cur_iter) * gradient

    return theta

X_b = np.hstack([np.ones((len(X), 1)), X])
initial_theta = np.zeros(X_b.shape[1])
theta = sgd(X_b, y, initial_theta, n_iters=m//3) #我们这里随机抽取样本总量的1/3来训练

print(theta)

输出结果如下：

我们可以看到在仅仅使用 1 / 3 的样本时，随机梯度下降的效果也非常好，大致也拟合出了我们示例数据的斜率和截距，可见随机梯度下降法的强大之处。

由于梯度改变方向是随机的，损失函数不能保证一直减小（可能损失函数是跳跃的）。正因如此，abs(J(theta, X_b, y) - J(last_theta, X_b, y)) < epsilon 这个条件用于随机梯度下降来说不是很好：这一次搜索和上一次搜索的差距特别小，并不代表离损失函数的中心更近。所以我们在代码中删去了这个退出循环的条件。

四、随机梯度下降法的封装与sklearn中提供的接口

在MyML的LinearRegression.py的LinearRegression类中，实现一个fit_sgd函数，具体代码如下：

    def fit_sgd(self, X_train, y_train, n_iters=50, t0=5, t1=50):
        #根据训练数据集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"
        assert n_iters >= 1

        def dJ_sgd(theta, X_b_i, y_i):
            return X_b_i * (X_b_i.dot(theta) - y_i) * 2.

        def sgd(X_b, y, initial_theta, n_iters=5, t0=5, t1=50):
            '''
            :param X_b: 测试集矩阵
            :param y: 测试集的标记矩阵
            :param initial_theta: θ的起始值
            :param n_iters: 对所有的样本看几遍
            :param t0: 学习率 η 的参数
            :param t1: 学习率 η 的参数
            '''
            def learning_rate(t):
                return t0 / (t + t1)

            theta = initial_theta
            m = len(X_b)
            for i_iter in range(n_iters): #进行n_iters次搜索
                indexes = np.random.permutation(m) #将样本乱序以便下一层for循环能够将训练集中每一个数据都以随机顺序遍历一遍
                X_b_new = X_b[indexes, :]
                y_new = y[indexes]
                for i in range(m): #每轮搜索进行m次
                    gradient = dJ_sgd(theta, X_b_new[i], y_new[i])
                    theta = theta - learning_rate(i_iter * m + i) * gradient #这是第i_iter * m + i 次搜索

            return theta

        X_b = np.hstack([np.ones((len(X_train), 1)), X_train])
        initial_theta = np.random.randn(X_b.shape[1])
        self._theta = sgd(X_b, y_train, initial_theta, n_iters, t0, t1)

        self.intercept_ = self._theta[0]
        self.coef_ = self._theta[1:]

        return self

在前面的编程中我们把循环次数(n_iters)设为总数除以3是不科学的，我们应该把所有的样本至少都看一遍，这样才能保证我们所有样本的信息都考虑进来了。基于这样的想法，上面代码的18至30行我们做了一些必要的调整：n_iters改为对所有的样本看几圈；第25行打乱顺序是为了保证在每一次大循环中，内循环对每一个数据都遍历了一次并且顺序随机。

在main.py中实现以下代码，以验证我们封装的正确性：

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)
X_test_standard = standardScaler.transform(X_test)

from MyML.LinearRegression import LinearRegression

lin_reg = LinearRegression()
lin_reg.fit_sgd(X_train_standard, y_train, n_iters=100)
print(lin_reg.score(X_test_standard, y_test))

运行结果如下：

结果已经非常接近于使用批量梯度下降法得出的准确度。

在sklearn中同样为我们封装了随机梯度下降的函数，调用如下：

from sklearn.linear_model import SGDRegressor

#sklearn中n_iters被重命名为max_iter
sgd_reg = SGDRegressor(max_iter=50, tol=1e-3)
sgd_reg.fit(X_train_standard, y_train)
print(sgd_reg.score(X_test_standard, y_test))

运行结果如下：