机器学习-梯度下降法的详细推导与代码实现

计算

对于线性回归，梯度下降法的目标就是找到一个足够好的向量$\theta$，使代价函数$J(\theta )={\sum }_{i=1}^m(\hat{y}-y_i{)}^2$取得最小值。线性回归的代价函数是关于$\theta$的多元函数。如下：
$$J(\theta )=\sum _{i=1}^m(\hat{y}-y^i{)}^2=\sum _{i=1}^m(\theta x^{(i)}-y^i{)}^2$$
$$J(\theta )=\sum _{i=1}^m(y^i-{\theta }_0-{\theta }_1{x}_1^{(i)}-{\theta }_2{x}_2^{(i)}-...-{\theta }_n{x}_n^{(i)}{)}^2$$

对代价函数$J(\theta )$的每一个$\theta$分量求偏导数，
$$\frac{\partial J(\theta )}{\partial {\theta }_i}=\left[\begin{array}{c}\frac{\partial J(\theta )}{\partial {\theta }_1}\\ \frac{\partial J(\theta )}{\partial {\theta }_2}\\ ⋮\\ \frac{\partial J(\theta )}{\partial {\theta }_n}\end{array}\right]=2\left[\begin{array}{c}{\sum }_{i=1}^m(y^{(i)}-X_b^{(i)}\theta )(-1)\\ {\sum }_{i=1}^m(y^{(i)}-X_b^{(i)}\theta )(-X_1^{(i)})\\ ⋮\\ {\sum }_{i=1}^m(y^{(i)}-X_b^{(i)}\theta )(-X_n^{(i)})\end{array}\right]=2\left[\begin{array}{c}{\sum }_{i=1}^m(X_b^{(i)}\theta -y^{(i)})\\ {\sum }_{i=1}^m(X_b^{(i)}\theta -y^{(i)})(X_1^{(i)})\\ ⋮\\ {\sum }_{i=1}^m(X_b^{(i)}\theta -y^{(i)})(X_n^{(i)})\end{array}\right]$$
这里将J()除以m可以缩小梯度的值：
$$J(\theta )=\frac{1}{m}\sum _{i=1}^m(\hat{y}-y^i{)}^2$$
此时梯度为：
$$\frac{\partial J(\theta )}{\partial {\theta }_i}=\frac{2}{m}\left[\begin{array}{c}{\sum }_{i=1}^m(X_b^{(i)}\theta -y^{(i)})\\ {\sum }_{i=1}^m(X_b^{(i)}\theta -y^{(i)})(X_1^{(i)})\\ ⋮\\ {\sum }_{i=1}^m(X_b^{(i)}\theta -y^{(i)})(X_n^{(i)})\end{array}\right]$$

• 此处也可以除以2m，因为可以方便求导数，不过，对结果没有什么影响，可以随意
• 这里的$J(\theta )$分为m个部分，因为有m个点。计算梯度的时候需要把m个点的坐标代进去关于$\theta$的表达式。这是梯度的标准求法。这个公式，是完整的梯度公式。即，在每个$\theta$方向求导数。得到的梯度向量就是下降最快的方向。
• 随机梯度下降法就是简化了这个$J(\theta )$，不需要再把m个点都代入进去，而是随机地抽取一个点，所以计算速度大为提高。
• 小批量梯度下降法就是中和了这两个算法，每次取batch个点代进去，计算梯度。BGD的梯度最后除以了m，这里的MBGD也可以除以batch。其实除和不除，都能计算出来，因为除的这个batch，在迭代的时候，会和$\alpha$相乘，只不过会对$\alpha$的取值有所影响。

以下是求梯度的函数，很简单，但是我总是记不住，这里作详细的介绍：

        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)

首先上面矩阵中的每个元素其实也需要通过矩阵运算得出
$$\sum _{i=1}^m(X_b^{(i)}\theta -y^{(i)})(X_1^{(i)})=(X^{(1)}\theta -y^{(1)})X_1^{(1)}+(X^{(2)}\theta -y^{(2)})X_1^{(2)}+...+(X^{(m)}\theta -y^{(m)})X_1^{(m)}$$
上面的式子可以用用两个向量的内积表示为，此处是内积，不是矩阵乘法:
$$\left[\begin{array}{c}{\theta }_0+X_1^{(1)}{\theta }_1+X_2^{(1)}{\theta }_2+X_n^{(1)}{\theta }_n...-y^{(1)}\\ {\theta }_0+X_1^{(2)}{\theta }_1+X_2^{(2)}{\theta }_2+X_n^{(2)}{\theta }_n...-y^{(2)}\\ ⋮\\ {\theta }_0+X_1^{(m)}{\theta }_1+X_2^{(m)}{\theta }_2+X_n^{(m)}{\theta }_n...-y^{(m)}\end{array}\right]\cdot \left[\begin{array}{c}{X}_1^{(1)}\\ X_1^{(2)}\\ ⋮\\ X_1^{(m)}\end{array}\right]$$
上面的式子又等价于：
$$(\left[\begin{array}{ccccc}1& X_1^{(1)}& X_2^{(1)}& \cdots & X_n^{(1)}\\ 1& X_1^{(2)}& X_2^{(2)}& \cdots & X_n^{(2)}\\ ⋮& ⋮& \ddots & ⋮& \\ 1& X_1^{(m)}& X_2^{(m)}& \cdots & X_n^{(m)}\end{array}\right]\cdot \left[\begin{array}{c}{\theta }_0\\ {\theta }_1\\ ⋮\\ {\theta }_n\end{array}\right]-\left[\begin{array}{c}{y}^{(1)}\\ y^{(2)}\\ ⋮\\ y^{(m)}\end{array}\right])\cdot \left[\begin{array}{c}{X}_1^{(1)}\\ X_1^{(2)}\\ ⋮\\ X_1^{(m)}\end{array}\right]$$

最终可以表示为，$X[:,i]$是python中的切片的写法：
$$(X_b\theta -y)\ast X[:,i]$$

因此，除了第一行以外，每一行的偏导数都可以表示为 $(X_b\theta -y)\ast X[:,i]$，上面的代码中的循环，就是对除了第一行以外的所有行，每一行都进行一次计算。这里dot()是点乘，向量点乘之后就自动求和了。所以不需要再次求和。

            for i in range(1, len(theta)):
                res[i] = (X_b.dot(theta) - y).dot(X_b[:, i])

所以其实是两种情况:

$(X_b\theta -y)$ 第一行
$(X_b\theta -y)\ast X[:,i]$ 之后的所有行

实现

思路：先准备数据集X和y ==> 使用交叉验证选出训练样本 ==> 回归
回归的过程: 将转换为X_b(就是加上X0) ==> 设置初始化 $\theta$ ==> 进行梯度下降(计算各个方向的偏导数-> 循环计算theta = theta - eta * dj)

class LinearRegression2:
    def __init__(self):
        self.coef_ = None        # 表示参数，theta_[1:]
        self.intercept_ = None   # 表示截距 ==>theta[0]
        self._thera = None       # 表示完整的theta==> theta[:]

    def fit_gd(self, X_train, y_train, eta=0.01, n_iters=1e4):
        """使用梯度下降法寻找最小的代价函数"""
        # 格式化X和theta，加上x0 和 theta0
        X_b = np.hstack([np.ones((len(X_train), 1)), X_train])
        initial_theta = np.zeros(X_b.shape[1])

        # 调用循环的梯度下降
        self._thera = self.gradient_descent(X_b, y_train, initial_theta, eta=eta, n_iters=n_iters)
        self.intercept_ = self._thera[0]
        self.coef_ = self._thera[1:]
        return self

    def gradient_descent(self, X_b, y, initial_theta, eta, n_iters=1e4, epsilon=1e-8):
        theta = initial_theta

        i = 0
        while i < n_iters:
            i += 1
            lastTheta = theta       # 记录上一个参数向量
            dj = self.dj(theta, X_b, y)
            theta = theta - eta * dj
            if np.absolute(self.J(lastTheta, X_b, y) - self.J(theta, X_b, y)) < epsilon:
                break
        return theta

    def dj(self, theta, X_b, y):
        """计算代价函数的偏导数"""
        # res 是 长度为len的一维数组
        res = np.zeros((len(theta)))        ##### zeros(5) == zeros((5,)) != zeros((5,1))
        res[0] = np.sum(X_b.dot(theta) - y)     # 需要将向量求和。没有办法，这里只能手动分开求解。因为这里X是从列方向分割，没有X0
        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 J(self, theta, X_b, y):
        """计算代价函数"""
        return np.sum((y - X_b.dot(theta)) ** 2) / len(y)

调用效果如下：

    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_train, X_test, y_train, y_test = train_test_split(X, y, test_ratio=0.2, seed=666)

    line = LinearRegression2()
    line.fit_gd(X_train, y_train)

    ###
    4.085675667692203
[   2.97732994]

使用向量实现

博客没法显示全公式。。就用图片凑合吧.$X_b\theta -y$是列向量

即，最后的公式为：
$$J(\theta )=\frac{2}{m}{X}_b^T\cdot (X_b\theta -y)$$
使用python的numpy实现：

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

其中，$X_b$是在原始数据的每一列前面加了一列后的矩阵，y是训练数据中给定的$y\_train$。

小结

• 不管是用向量计算还是用手动计算。都需要先把训练数据$X$转换为$X_b$
• 一般来说，$X$都是ndarray类型的二维数组，使用hstack([ist1, list2])函数可以将两个数组摊在一起。将$X$转换为$X_b$的代码为：

  X_b = np.hstack([np.ones((len(X_train), 1)), X_train])

• 大体思路就是:
  1. 调用fit()方法拟合，该方法中产生$X_b$，并初始化$\theta$。最后使调用梯度下降方法gradient_descent()来找到最优的$\theta$
  2. gradient_descent() 循环调用dj()计算对$\theta$的偏导数,并每一次都对$\theta$的值进行更新。
  3. gradient_descent() 中会调用J()来判断梯度的增量是否已经足够小。

数据标准化

在这里使用波士顿房价数据测试，发现一些很有意思的东西，如下。为什么会这样？，因为数据没有归一化，在数据集中存在一下特别大和特别小的数字，导致在计算梯度时，可能会导致梯度的跨度太大，而无法收敛。也有可能时计算式出现了过大的数inf。

    linear.fit_gd(X, y, n_iters=1e4, eta=0.001)
    print(linear.coef_)

    ###
    [nan nan nan nan nan nan nan nan nan nan nan nan nan]

只要将学习系数设定的足够小，就不会报错。

linear.fit_gd(X, y, n_iters=1e4, eta=0.0000001)

使用sklearn的StandardScaler归一化数据。注意归一化需要将所有的X都归一化，包括用来训练的x_train和用来测试的x_test

    from sklearn.preprocessing import StandardScaler

    standardScaler = StandardScaler()
    standardScaler.fit(X)                       # 导入数据
    X_standard = standardScaler.transform(X)    # 转换数据

所以，到目前为止，要预测一个波士顿房价数据的完整过程是这样的：

    from sklearn import datasets

    # 加载数据集
    X = datasets.load_boston().data
    y = datasets.load_boston().target

    # 剔除噪音
    X = X[y < 50]
    y = y[y < 50]

    # 数据归一化处理
    from sklearn.preprocessing import StandardScaler

    standardScaler = StandardScaler()
    standardScaler.fit(X)
    X_standard = standardScaler.transform(X)

    # 交叉验证
    X_train_standard, X_test_standard, y_train, y_test = train_test_split(X_standard, y, test_ratio=0.2, seed=666)

    # 回归，拟合
    linear = LinearRegression2()
    linear.fit_gd(X_train_standard, y_train, n_iters=3e5)

    # 使用score计算拟合的效果
    print(linear.score(X_test_standard, y_test))

随机梯度下降法

原本的梯度下降法：
$$\frac{\partial J(\theta )}{\partial {\theta }_i}=\frac{2}{m}\left[\begin{array}{c}{\sum }_{i=1}^m(X_b^{(i)}\theta -y^{(i)})\\ {\sum }_{i=1}^m(X_b^{(i)}\theta -y^{(i)})(X_1^{(i)})\\ ⋮\\ {\sum }_{i=1}^m(X_b^{(i)}\theta -y^{(i)})(X_n^{(i)})\end{array}\right]$$

随机梯度下降法：
$$\text{ 2. }\left(\begin{array}{c}\left(X_b^{(i)}\theta -y^{(i)}\right)\cdot X_0^{(i)}\\ \left(X_b^{(i)}\theta -y^{(i)}\right)\cdot X_1^{(i)}\\ \left(X_b^{(i)}\theta -y^{(i)}\right)\cdot X_2^{(i)}\\ \cdots \\ \left(X_b^{(i)}\theta -y^{(i)}\right)\cdot X_n^{(i)}\end{array}\right)=2\cdot {\left(X_b^{(i)}\right)}^T\cdot \left(X_b^{(i)}\theta -y^{(i)}\right)$$

原本的计算偏导数的函数用这个替代。里只有一个样本哦，所有的误差和X都是一个同一个样本的。而原来的函数，每计算一次偏导数都要对n个特征，每个循环m次，共mn次计算。计算量确实是大大提高，但是，不能保证每次都找到最好的梯度。使用随机产生的梯度。相比较原来的函数，去掉了步长。
学习率：
$$\eta =\frac{t_0}{i_{-}iters+t_1}$$

• 去掉了参数中的步长
• n_iters 表示遍历所有样本的轮数
• 学习率需要不断缩小，因为，算出来的梯度不能保证是越来越小的（因为随机）

class StochasticDescent:
    """随机梯度下降法"""

    def __init__(self):
        self.coef_ = None       # 表示参数，theta_[1:]
        self.intercept_ = None  # 表示截距 ==>theta[0]
        self._thera = None      # 表示完整的theta==> theta[:]

    def fit_SGD(self, X_train, y_train, eta=0.01, n_iters=1):
        """使用梯度下降法寻找最小的代价函数"""
        # 格式化X和theta，加上x0 和 theta0
        X_b = np.hstack([np.ones((len(X_train), 1)), X_train])
        initial_theta = np.zeros(X_b.shape[1])

        # 调用循环的梯度下降
        self._thera = self.sgd(X_b, y_train, initial_theta, n_iters, 5, 50)
        self.intercept_ = self._thera[0]
        self.coef_ = self._thera[1:]
        return self

    def sgd(self, X_b, y_train, initial_theta, n_iters, t0=5, t1=50):
        """随机梯度下降, niters是轮数"""

        def get_study_rate(i_iters):
            """把学习率和迭代次数联系起来"""
            return t0 / (t1 + i_iters)

        def dJ_sgd(X_b_i, theta, y):
            """计算随机一个元素的梯度"""
            return 2 * X_b_i.T.dot(X_b_i.dot(theta) - y)

        m = len(X_b)
        theta = initial_theta
        for n in range(int(n_iters)):
            indexes = np.random.permutation(m)
            X_b_new = X_b[indexes, :]
            y_new = y[indexes]
            for i in range(m):
                sgd = dJ_sgd(X_b_new[i], theta, y_new[i])
                theta = theta - get_study_rate(m * n + i) * sgd
        return theta

sklearn 中的随机梯度下降

SGDRegressor位于线性模型下，SGDRegressor(max_iter=50)可以传入参数 max_iter指定迭代的次数。

    from sklearn.linear_model import SGDRegressor

    sgd_reg = SGDRegressor(max_iter=50)
    sgd_reg.fit(X, y)

调试梯度下降法

对每个$\theta$分量，都要微积分，最后，用算出来的微积分向量，作为梯度。

• 计算量巨大，调试成功后，要关闭
• 分子的两个顺序不能错，我把顺序搞错了，结果导数算反了。。。

$$\frac{dJ}{d\theta }=\frac{J(\theta +\epsilon )-J(\theta -\epsilon )}{2\epsilon }$$

关键代码在这：

    def dj_debug(self, X_b, theta, y, epsilon=0.001):
        res = np.zeros((len(theta), ))
        for i in range(len(theta)):
            theta_1 = theta.copy()
            theta_2 = theta.copy()
            theta_1[i] += epsilon
            theta_2[i] -= epsilon
            res[i] = (self.J(theta_1, X_b, y) - self.J(theta_2, X_b, y)) / (2 * epsilon)
        return res

完整代码：

class aDebugGradient:
    """随机梯度下降法"""

    def __init__(self):
        self.coef_ = None       # 表示参数，theta_[1:]
        self.intercept_ = None  # 表示截距 ==>theta[0]
        self._thera = None      # 表示完整的theta==> theta[:]

    def fit_debug(self, X_train, y_train, eta=0.01, n_iters=1e4):
        """使用梯度下降法寻找最小的代价函数"""
        # 格式化X和theta，加上x0 和 theta0
        X_b = np.hstack([np.ones((len(X_train), 1)), X_train])
        initial_theta = np.zeros(X_b.shape[1])

        # 调用循环的梯度下降
        self._thera = self.gradient_Descent(self.dj_debug, X_b, y_train, initial_theta, eta, n_iters)
        self.intercept_ = self._thera[0]
        self.coef_ = self._thera[1:]
        return self

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

    def dj_debug(self, X_b, theta, y, epsilon=0.001):
        res = np.zeros((len(theta), ))
        for i in range(len(theta)):
            theta_1 = theta.copy()
            theta_2 = theta.copy()
            theta_1[i] += epsilon
            theta_2[i] -= epsilon
            res[i] = (self.J(theta_1, X_b, y) - self.J(theta_2, X_b, y)) / (2 * epsilon)
        return res

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

    def gradient_Descent(self, dj, X_b, y_train, initial_theta, eta, n_iters,  epsilon=1e-8):
        """梯度下降"""
        theta = initial_theta
        for i in range(int(n_iters)):
            gradient = dj(X_b, theta, y_train)      # 计算得到梯度
            last_theta = theta
            theta = theta - eta * gradient
            if np.absolute(self.J(theta, X_b, y_train) - self.J(last_theta, X_b, y_train)) < epsilon:
                break
        return theta

最后，，，这垃圾公式真是烦死人！！！