ML-梯度下降法的详细推导与代码实现

计算

对于线性回归，梯度下降法的目标就是找到一个足够好的向量( heta)，使代价函数(J( heta) = sum_{i=1}^{m}(hat{y}-y_{i})^{2})取得最小值。线性回归的代价函数是关于( heta)的多元函数。如下：

[J( heta) = sum_{i=1}^{m}(hat{y}-y^{i})^{2} = sum_{i=1}^{m}( heta x^{(i)}-y^{i})^{2} ]

[J( heta) = sum_{i=1}^{m}(y^{i} - heta_{0} - heta_{1}x^{(i)}_1 - heta_{2}x^{(i)}_2 - ... - heta_{n}x^{(i)}_n )^{2} ]

对代价函数(J( heta))的每一个( heta)分量求偏导数，

[frac{partial J( heta)}{partial heta_i}= egin{bmatrix} frac{partial J( heta)}{partial heta_1} \ frac{partial J( heta)}{partial heta_2} \ vdots \ frac{partial J( heta)}{partial heta_n} end{bmatrix} = 2 egin{bmatrix} sum_{i=1}^{m}(y^{(i)} - X^{(i)}_b heta) (-1) \ sum_{i=1}^{m}(y^{(i)} - X^{(i)}_b heta) (-X^{(i)}_1) \ vdots \ sum_{i=1}^{m}(y^{(i)} - X^{(i)}_b heta) (-X^{(i)}_n) end{bmatrix} = 2 egin{bmatrix} sum_{i=1}^{m}(X^{(i)}_b heta -y^{(i)}) \ sum_{i=1}^{m}(X^{(i)}_b heta -y^{(i)}) (X^{(i)}_1) \ vdots \ sum_{i=1}^{m}(X^{(i)}_b heta - y^{(i)}) (X^{(i)}_n) end{bmatrix} ]

这里将J()除以m可以缩小梯度的值：

[J( heta) = frac{1}{m}sum_{i=1}^{m}(hat{y}-y^{i})^{2} ]

此时梯度为：

[frac{partial J( heta)}{partial heta_i}= frac{2}{m} egin{bmatrix} sum_{i=1}^{m}(X^{(i)}_b heta -y^{(i)}) \ sum_{i=1}^{m}(X^{(i)}_b heta -y^{(i)}) (X^{(i)}_1) \ vdots \ sum_{i=1}^{m}(X^{(i)}_b heta - y^{(i)}) (X^{(i)}_n) end{bmatrix} ]

• 此处也可以除以2m，因为可以方便求导数，不过，对结果没有什么影响，可以随意
• 这里的(J( heta))分为m个部分，因为有m个点。计算梯度的时候需要把m个点的坐标代进去关于( heta)的表达式。这是梯度的标准求法。这个公式，是完整的梯度公式。即，在每个( heta)方向求导数。得到的梯度向量就是下降最快的方向。
• 随机梯度下降法就是简化了这个(J( heta))，不需要再把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^{(i)}_b heta -y^{(i)}) (X^{(i)}_1) = (X^{(1)} heta - y^{(1)})X^{(1)}_1 + (X^{(2)} heta - y^{(2)})X^{(2)}_1 + ... + (X^{(m)} heta - y^{(m)})X^{(m)}_1 ]

上面的式子可以用用两个向量的内积表示为，此处是内积，不是矩阵乘法:

[egin{bmatrix} heta_0 + X^{(1)}_1 heta_1 + X^{(1)}_2 heta_2 + X^{(1)}_n heta_n ... - y^{(1)} \ heta_0 + X^{(2)}_1 heta_1 + X^{(2)}_2 heta_2 + X^{(2)}_n heta_n ... - y^{(2)} \ vdots \ heta_0 + X^{(m)}_1 heta_1 + X^{(m)}_2 heta_2 + X^{(m)}_n heta_n ... - y^{(m)} end{bmatrix} cdot egin{bmatrix} X^{(1)}_1 \ X^{(2)}_1 \ vdots \ X^{(m)}_1 end{bmatrix} ]

上面的式子又等价于：

[(egin{bmatrix} 1 & X^{(1)}_{1} & X^{(1)}_{2} & cdots & X^{(1)}_{n} \ 1 & X^{(2)}_{1} & X^{(2)}_{2} & cdots & X^{(2)}_{n} \ vdots & vdots & ddots & vdots \ 1 & X^{(m)}_{1} & X^{(m)}_{2} & cdots & X^{(m)}_{n} end{bmatrix} cdot egin{bmatrix} heta_0 \ heta_1 \ vdots \ heta_n end{bmatrix} - egin{bmatrix} y^{(1)} \ y^{(2)} \ vdots \ y^{(m)} end{bmatrix})cdot egin{bmatrix} X^{(1)}_1 \ X^{(2)}_1 \ vdots \ X^{(m)}_1 end{bmatrix} ]

最终可以表示为，(X[:,i])是python中的切片的写法：

[(X_b heta - y)* X[:,i] ]

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

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

所以其实是两种情况:

((X_b heta - y)) 第一行
((X_b heta - y)* X[:,i]) 之后的所有行

实现

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

即，最后的公式为：

[J( heta) = frac{2}{m}X^T_b cdot (X_b heta-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)，并初始化( heta)。最后使调用梯度下降方法gradient_descent()来找到最优的( heta)
  2. gradient_descent() 循环调用dj()计算对( heta)的偏导数,并每一次都对( heta)的值进行更新。
  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( heta)}{partial heta_i}= frac{2}{m} egin{bmatrix} sum_{i=1}^{m}(X^{(i)}_b heta -y^{(i)}) \ sum_{i=1}^{m}(X^{(i)}_b heta -y^{(i)}) (X^{(i)}_1) \ vdots \ sum_{i=1}^{m}(X^{(i)}_b heta - y^{(i)}) (X^{(i)}_n) end{bmatrix} ]

随机梯度下降法：

[ ext { 2. }left(egin{array}{c}{left(X_{b}^{(i)} heta-y^{(i)} ight) cdot X_{0}^{(i)}} \ {left(X_{b}^{(i)} heta-y^{(i)} ight) cdot X_{1}^{(i)}} \ {left(X_{b}^{(i)} heta-y^{(i)} ight) cdot X_{2}^{(i)}} \ {cdots} \ {left(X_{b}^{(i)} heta-y^{(i)} ight) cdot X_{n}^{(i)}}end{array} ight)=2 cdotleft(X_{b}^{(i)} ight)^{T} cdotleft(X_{b}^{(i)} heta-y^{(i)} ight) ]

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

[eta=frac{t_{0}}{i_{-} i t e r s+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)

调试梯度下降法

对每个( heta)分量，都要微积分，最后，用算出来的微积分向量，作为梯度。

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

[frac{d J}{d heta}=frac{J( heta+varepsilon)-J( heta-varepsilon)}{2 varepsilon} ]

关键代码在这：

    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

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