机器学习之梯度下降法与线性回归

摘要：表示学习率，是梯度下降法的一个超参数，其取值影响最优解的速度。因此在使用梯度下降法之前，最好进行数据归一化。同时在随机梯度下降法中学习率的取值是逐渐递减的，为了防止固定取值的学习率使得梯度下降到达最优解附近时继续跳出这个范围。

梯度下降法不是一个机器学习算法，而是一种基于搜索的最优化方法，用于最小化一个效用函数。

假设存在一个只有一个参数 $ heta$ 的损失函数 $J$，想找到最小极值处的 $ heta$，如图所示：

借助于损失函数 $J$ 在 $ heta$ 处的切线，可以直观的反映出损失函数 $J$ 在 $ heta$ 处的导数大小；导数的大小代表着 $ heta$ 变化时 $J$ 相应的变化。

同时导数也可以代表 $J$ 增大的方向，如果将导数与 $-eta$ 相乘，即 $-etafrac{dJ}{d heta}$ 代表了 $J$ 减小的方向。

$eta$ 表示学习率，是梯度下降法的一个超参数，其取值影响最优解的速度。太小会减慢收敛学习速度，太大可能导致不收敛。

如果 $J$ 中存在多处导数为零的情况，即存在一个全局最优解和多个局部最优解，此时可以多次使用梯度下降法，每次随机化一个初始点。

对于有多个 $ heta$ 的 $J$ 类似，即找出全局最优解处的这些 $ heta$ 的值。

首先，模拟一个损失曲线 $J$：

import numpy as np

plot_x = np.linspace(-1, 6, 141)
plot_y = (plot_x - 2.5) ** 2 - 1

作图表示如下：

定义函数 dJ() 用于求 $J$ 在 $ heta$ 处的导数：

def dJ(theta):
    return 2 * (theta - 2.5)

函数 J() 用于求 $J$ 在 $ heta$ 处的大小：

def J(theta):
    return (theta - 2.5) ** 2 - 1

接着使用梯度下降法，首先给 $ heta$ 赋一个初值 0 及学习率 $eta$ 为 0.1，接着在循环里进行多次梯度下降。每次循环都要求得 $J$ 在 $ heta$ 处的导数值 $gradient$，并且 $ heta$ 向导数的负方向移动，即：$ heta= heta-eta*gradient$。

由于计算机计算浮点数存在误差，对于求得的 $ heta$ 可能不能刚好等于 0，因此设定一个精度值（epsilon = 1e-8），如果新的 $ heta$ 对应的损失函数的值与上一次 $ heta$ 对应的损失函数的值的差值满足精度要求，就表示找到了要找的 $ heta$。程序如下：

theta = 0.0
eta = 0.1
epsilon = 1e-8

while True:
    gradient = dJ(theta)
    last_theta = theta
    
    theta = theta - eta * gradient
    if (abs(J(theta) - J(last_theta)) < epsilon):
        break

运行程序求得的 $ heta$ 为：2.499891109642585。

对于 $ heta$ 的取值变化，可以用图片表示，如下（红色的点）：

对于学习率 $eta$，这里取值 0.1 是没有问题的，但如果取值 1.1 程序运行就会报错：

OverflowError                             Traceback (most recent call last)
 in 
      8 
      9     theta = theta - eta * gradient
---> 10     if (abs(J(theta) - J(last_theta)) < epsilon):
     11         break
     12 

 in J(theta)
      1 def J(theta):
----> 2     return (theta - 2.5) ** 2 - 1

OverflowError: (34, "Result too large")

这是因为学习率过大会导致 J(theta) 越来越大。为了使程序不会报错，修改 J() 方法：

def J(theta):
    try:
        return (theta - 2.5) ** 2 - 1
    except:
        return float("inf")

注意，当无穷减无穷时，结果时 nan 而不是 0，此时 if (abs(J(theta) - J(last_theta)) < epsilon) 将永远无法触发而使得程序进入死循环。为了解决这个问题，增加一个新的超参数 n_iters，表示能够执行循环的最大次数。

比如使 n_iters=10，$ heta$ 取值变化如图：

最后，把梯度下降法封装到方法中：

def gradient_descent(initial_theta, eta, n_iters=1e4, epsilon=1e-8):
    theta = initial_theta
    i_ters = 0

    while i_ters < n_iters:
        gradient = dJ(theta)
        last_theta = theta

        theta = theta - eta * gradient
        if (abs(J(theta) - J(last_theta)) < epsilon):
            break
        
        i_ters += 1
    
    return theta

多元线性回归的损失函数为：

$$ J=sum_{i=1}^{m}(y^{(i)} - hat{y}^{(i)})^2 $$

其中：$hat{y}^{(i)} = heta_{0} + heta_{1}X_{1}^{(i)} + heta_{2}X_{2}^{(i)} + ... + heta_{n}X_{n}^{(i)}$ 。

对 $J$ 求导为：

$$ abla J=(frac{partial J}{partial heta_0},frac{partial J}{partial heta_1},...,frac{partial J}{partial heta_n}) $$

其中：$frac{partial J}{partial heta_i}$ 为偏导数，与导数的求法一样。

对 $ abla J$ 进一步计算：

$$ abla J( heta) = egin{pmatrix} frac{partial J}{partial heta_0} frac{partial J}{partial heta_1} frac{partial J}{partial heta_2} cdots frac{partial J}{partial heta_n} end{pmatrix} = egin{pmatrix} sum_{i=1}^{m}2(y^{(i)} - X_b^{(i)} heta)·(-1) sum_{i=1}^{m}2(y^{(i)} - X_b^{(i)} heta)·(-X_1^{(i)}) sum_{i=1}^{m}2(y^{(i)} - X_b^{(i)} heta)·(-X_2^{(i)}) cdots sum_{i=1}^{m}2(y^{(i)} - X_b^{(i)} heta)·(-X_n^{(i)}) end{pmatrix} = 2·egin{pmatrix} sum_{i=1}^{m}(X_b^{(i)} heta - y^{(i)}) sum_{i=1}^{m}(X_b^{(i)} heta - y^{(i)})·X_1^{(i)} sum_{i=1}^{m}(X_b^{(i)} heta - y^{(i)})·X_2^{(i)} cdots sum_{i=1}^{m}(X_b^{(i)} heta - y^{(i)})·X_n^{(i)} end{pmatrix} $$

其中：$X_b = egin{pmatrix} 1 & X_1^{(1)} & X_2^{(1)} & cdots & X_n^{(1)} 1 & X_1^{(2)} & X_2^{(2)} & cdots & X_n^{(2)} cdots & & & & cdots 1 & X_1^{(m)} & X_2^{(m)} & cdots & X_n^{(m)} end{pmatrix}$

这个结果是与样本数量 m 相关的，为了使结果与 m 无关，对这个梯度除以 m，即：

$$ abla J( heta) = frac{2}{m}·egin{pmatrix} sum_{i=1}^{m}(X_b^{(i)} heta - y^{(i)}) sum_{i=1}^{m}(X_b^{(i)} heta - y^{(i)})·X_1^{(i)} sum_{i=1}^{m}(X_b^{(i)} heta - y^{(i)})·X_2^{(i)} cdots sum_{i=1}^{m}(X_b^{(i)} heta - y^{(i)})·X_n^{(i)} end{pmatrix} $$

此时，目标函数就成了使 $frac{1}{m}sum_{i=1}^{m}(y^{(i)} - hat{y}^{(i)})^2$ 尽可能小，即均方误差尽可能小：

$$ J( heta) = MSE(y, hat{y}) $$

首先模拟训练数据集：

import numpy as np

x = np.random.random(size=100)
y = x * 3.0 + 4.0 + np.random.normal(size=100)

X = x.reshape(-1, 1)

定义函数 J() 计算损失函数的值：

def J(theta, X_b, y):
    try:
        return np.sum((y - X_b.dot(theta)) ** 2) / len(X_b)
    except:
        return float("inf")

函数 dJ() 对 $ heta$ 求导数：

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)

注意：对 $J$ 求导更好的方式是进行向量化处理，即 $ abla J( heta) = frac{2}{m}·X_b^T·(X_b heta-y)$，dJ() 改写为：

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

梯度下降的过程为：

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

    while i_ters < 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
        
        i_ters += 1
        
    return theta

执行程序找出最优的 $ heta$ 如下：

X_b = np.hstack([np.ones((len(X), 1)), X])
initial_theta = np.zeros(X_b.shape[1])
eta = 0.01

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

$ heta$ 结果为：

array([4.0269033, 3.0043078])

将梯度下降法封装到线性回归算法的模型训练方法 fit_gd() 中：

class LinearRegression:
    # other codes here
    
    def fit_gd(self, X_train, y_train, eta=0.01, n_iters=1e4):
        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):
            return X_b.T.dot(X_b.dot(theta) - y) * 2 /len(X_b)

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

            while i_ters < 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
                
                i_ters += 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)
        self.interception_ = self._theta[0]
        self.coef_ = self._theta[1:]

        return self

注意：
在真实的数据集 (X, y) 中，X 整体不在一个规模上会影响梯度的结果，而梯度的结果再乘以 $eta$ 得到步长就太大或太小，从而导致训练出的模型可能很差。因此在使用梯度下降法之前，最好进行数据归一化。

$ abla J( heta) = frac{2}{m}·egin{pmatrix} sum_{i=1}^{m}(X_b^{(i)} heta - y^{(i)}) sum_{i=1}^{m}(X_b^{(i)} heta - y^{(i)})·X_1^{(i)} sum_{i=1}^{m}(X_b^{(i)} heta - y^{(i)})·X_2^{(i)} cdots sum_{i=1}^{m}(X_b^{(i)} heta - y^{(i)})·X_n^{(i)} end{pmatrix}$ 中每一项都要对所有样本进行计算，因此这种梯度下降法称为批量梯度下降法（Batch Gradient Descent）。如果 m 非常大，使用批量梯度下降法计算梯度的计算量就会非常大。

改进的方法是对每一项计算时不再计算所有的样本而是选取其中一个样本进行计算，即：

$$ 2·egin{pmatrix} (X_b^{(i)} heta - y^{(i)})·X_0^{(i)} (X_b^{(i)} heta - y^{(i)})·X_1^{(i)} (X_b^{(i)} heta - y^{(i)})·X_2^{(i)} cdots (X_b^{(i)} heta - y^{(i)})·X_n^{(i)} end{pmatrix} = 2·(X_b^{(i)})^T·(X_b^{(i)} heta-y^{(i)}) $$

这样的方式就是随机梯度下降法（Stochastic Gradient Descent），此时搜索路径如图所示：

随机梯度下降法不能保证梯度下降的方向就是损失函数减小最快的方向（有时会是增大的方向），但整体上是会到达最小值附近的（满足一定的精度）。

同时在随机梯度下降法中学习率 $eta$ 的取值是逐渐递减的，为了防止固定取值的学习率使得梯度下降到达最优解附近时继续跳出这个范围。一个合理的取值方式为：

$$ eta = frac{t0}{i\_iter + t1} $$

其中：$t0、t1$ 为超参数。

定义函数 dJ_sgd() 对应批量梯度下降法中对损失函数求导的过程，此时传入函数的就不再是所有样本 $(X_b, y)$ 了，而是其中一个样本 $(X_b^{(i)}, y^{(i)})$：

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

函数 sgd() 中定义了一个 learning_rate() 方法用来计算学习率，传入参数为当前迭代次数。

因为要随机的选取样本中的一个，但又要将所有样本看一遍，所以我们将这个样本集打乱形成一个新的样本集 $(X_{b,new}^{(i)}, y_{new}^{(i)})$，同时指定参数 n_iters 表示将这个样本集看几遍：

def sgd(X_b, y, initial_theta, n_iters, t0, t1):

    def learning_rate(t):
        return t0 / (t + t1)

    theta = initial_theta
    m = len(X_b)

    for cur_iter in range(n_iters):
        indexes = np.random.permutation(m)
        X_b_new = X_b[indexes]
        y_new = y[indexes]

        for i in range(m):
            grandient = dJ_sgd(theta, X_b_new[i], y_new[i])
            theta = theta - learning_rate(cur_iter * m + i) * grandient

            return theta

将随机梯度下降法封装到线性回归算法的模型训练方法 fit_sgd() 中：

class LinearRegression:
    # other codes here
    
    def fit_sgd(self, X_train, y_train, n_iters=5, t0=5, t1=50):
        def dJ_sgd(theta, X_b_i, y_i):
            return X_b_i.T.dot(X_b_i.dot(theta) - y_i) * 2.

        def sgd(X_b, y, initial_theta, n_iters, t0, t1):
            
            def learning_rate(t):
                return t0 / (t + t1)
            
            theta = initial_theta
            m = len(X_b)

            for cur_iter in range(n_iters):
                indexes = np.random.permutation(m)
                X_b_new = X_b[indexes]
                y_new = y[indexes]

                for i in range(m):
                    grandient = dJ_sgd(theta, X_b_new[i], y_new[i])
                    theta = theta - learning_rate(cur_iter * m + i) * grandient
                
            return theta

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

        self._theta = sgd(X_b, y_train, initial_theta, n_iters, t0, t1)
        self.interception_ = self._theta[0]
        self.coef_ = self._theta[1:]

        return self

在 Scikit Learn 的 linear_model 模块中提供了一个使用随机梯度下降法的回归算法 SGDRegressor：

from sklearn.linear_model import SGDRegressor

Github | ML-Algorithms-Action