梯度下降法解决线性回归问题

简介

梯度下降法是一种寻找目标函数最小值的优化算法

原理

假设一个函数 $f(x)$，在某点 $P_0(x_0)$ 处的值为 $f(x_0)$ ，导数为值 $\Delta{f}(x_0)$。由于要求函数最小值，导数表示$f(x)$的变化率，为正表示当 $x$ 增大 $f(x)$增大，反之减小。所以要寻找到下一个点 $x_1$，则

$$x_1 = x_0 - \eta\Delta(f(x_0))$$

$x$沿着使$f(x)$减小的方向前进，其中 $\eta$ 表示步长，或者说学习率。当 $\Delta{f(x)} = f(xn) - f(x{n-1})$ 的绝对值小于某个指定的数值，如 $0.00001$时，表示已经沿着梯度下降的很吃力了，已经逼近极小值了。

简单的代码示例

示例目标函数为 $y = (x - 2.5)^2 + 1$

 # 函数的导数 df
df = lambda x: 2 * (x - 2.5)
f = lambda x: (x - 2.5) ** 2 + 1
eta = 0.1
epsilon = 1e-8
 # 最大循环次数
max_iters=100000
 # 下降的起点
cur_x = 0.0
cur_iter = 0
while cur_iter < max_iters:
    prev_x = cur_x
    cur_x -= eta * df(cur_x)
    if abs(f(cur_x) - f(prev_x)) < epsilon:
        break
cur_iter += 1
print(cur_x)

打印结果为: 2.499891109642585，根据函数知道$x$ 在 $2.5$ 处取得最小值，结果十分接近。

线性回归中的梯度下降法

在上面的文章中谈到，线性回归的方程为

$$\hat{y}^i = \theta_0 + \theta_1x_1^i + \theta_2x_2^i + \theta_3x_3^i + \dots + \theta_nx_n^i$$

的损失函数为

$$J = \frac{1}{m}\sum_{i=1}^m(y_i - \hat{y_i})^2$$ $$J = \frac{1}{m}\sum_{i=1}^m(y_i - \theta_0 - \theta_1x_1^i - \theta_2x_2^i - \theta_3x_3^i - \dots - \theta_nx_n^i)^2$$

此时，每个点的梯度方向就是对每个维度求偏导，形成的向量

$$\Delta{J} = \begin{pmatrix} \frac{\phi(J)}{\phi(\theta_0)} \\ \frac{\phi(J)}{\phi(\theta_1)} \\ \frac{\phi(J)}{\phi(\theta_2)} \\ \cdots \\ \frac{\phi(J)}{\phi(\theta_n)} \end{pmatrix}$$ $$\Delta{J} = \frac{2}{m} \begin{pmatrix} \sum_{i=1}^m(y^i - \theta_0 - \theta_1X_1^i - \theta_2X_2^i - \theta_3X_3^i - \dots - \theta_nX_n^i)(-1) \\ \sum_{i=1}^m(y^i - \theta_0 - \theta_1X_1^i - \theta_2X_2^i - \theta_3X_3^i - \dots - \theta_nX_n^i)(-X_1^i) \\ \cdots \\ \sum_{i=1}^m(y^i - \theta_0 - \theta_1X_1^i - \theta_2X_2^i - \theta_3X_3^i - \dots - \theta_nX_n^i)(-X_n^i) \end{pmatrix}$$

将数据集添加第一列全是1，形成向量 $X_b$

$$\Delta{J} = \frac{2}{m} \begin{pmatrix} \sum_{i=1}^m(X_b^i\theta - y^i) \\ \sum_{i=1}^m(X_b^i\theta - y^i)X_1^i \\ \cdots \\ \sum_{i=1}^m(X_b^i\theta - y^i)X_n^i \\ \end{pmatrix}$$

此时已经得到$\Delta{J}$

• 我们选初始化一个随机的数据集

  x = np.linspace(-3, 10, 30)
  y = x * 2 + 3 + np.random.normal(0, 2, size=30)
  # 转换成矩阵
  X = x.reshape(-1, 1)

将x和y点散点图画出来就是



• 目标函数和梯度函数

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

• 梯度下降函数

  def gradient_descent( X_b, y, init_theta, eta, epsilon=1e-8, max_iters=1e5):
      count = 0
      theta = init_theta
      while count < max_iters:
          gradient = dJ(theta, X_b, y)
          pre = theta
          theta = theta - eta * gradient
          if abs(J(theta, X_b, y) - J(pre, X_b, y)) < epsilon:
              break
          count += 1
      return theta

• 使用测试的数据进行操作

  # 将原矩阵添加一列1，使得与theta的维度相同
  X_b = np.hstack([np.ones((len(X), 1)), X])
  init_theta = np.zeros(X_b.shape[1])
  eta = 0.01
  theta = gradient_descent(X_b, y, init_theta, eta)

求出点theta值为 $[2.26867578, 2.20822388]$，将原数据点和 $[2.26867578, 2.20822388]$构成点直线作图，获得

• 向量化的处理

  $$\Delta{J} = \frac{2}{m} \begin{pmatrix} \sum_{i=1}^m(X_b^i\theta - y^i) \\ \sum_{i=1}^m(X_b^i\theta - y^i)X_1^i \\ \cdots \\ \sum_{i=1}^m(X_b^i\theta - y^i)X_n^i \\ \end{pmatrix} = \frac{2}{m} \begin{pmatrix} \sum_{i=1}^m(X_b^i\theta - y^i)X_0^i \\ \sum_{i=1}^m(X_b^i\theta - y^i)X_1^i \\ \cdots \\ \sum_{i=1}^m(X_b^i\theta - y^i)X_n^i \\ \end{pmatrix}$$ $$= \frac{2}{m} \cdot X_b^T \cdot (X_b\cdot \theta - y)$$

  求梯度函数修改为

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

• StandardScaler 预处理数据

  由于boston房价数据各个特征的单位不同，数值大小差异较大，故必须进行归一化处理，使用 StandardScaler 处理，在调用fit方法前

  from sklearn.preprocessing import StandardScaler
  scaler = StandardScaler()
  scaler.fit(X_train)
  X_train = scaler.transform(X_train)
  X_test = scaler.transform(X_test)

  • 随机梯度下降法 Stochastic Gradient Descent

  每次计算梯度都要用到数据集中所有的样本，如果样本量过大，会导致速度受到很大的影响，假设将方程

  $$\Delta{J} = \frac{2}{m} \begin{pmatrix} \sum_{i=1}^m(X_b^i\theta - y^i)X_0^i \\ \sum_{i=1}^m(X_b^i\theta - y^i)X_1^i \\ \cdots \\ \sum_{i=1}^m(X_b^i\theta - y^i)X_n^i \\ \end{pmatrix}$$

  中的 $X_b^i$ 随机为某个样本，那么m个相同的值相加再除以m，可以得到

  $$\Delta{J} = 2 \begin{pmatrix} (X_b^i\theta - y^i)X_0^i \\ (X_b^i\theta - y^i)X_1^i \\ \cdots \\ (X_b^i\theta - y^i)X_n^i \\ \end{pmatrix}$$ $$= 2 \cdot (X_b^i)^T \cdot (X_b^i\cdot \theta - y^i)$$

  随机梯度下降的趋势如下，下降的方向并不是梯度的方向，但是总体的趋势是向最小值的方向前进

对于 $\eta$ 的取值，在刚开始时候较大，随着越来越靠近最小区域，$\eta$ 的值应该慢慢减小，故设定 $\eta$ 的获取函数

$$\eta = \frac{a}{iters + b}$$

其中 $a, b$ 两个是算法的超参数，如果不去调节该值，可以使用经验值 $a = 5, b = 50$, $iters$ 是循环的次数，随着次数的增多，$\eta$ 越来越小

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

def gradient_descent_sgd(X_b, y, initial_theta, n_iters=1e5, epsilon=1e-8):
    t0 = 5
    t1 = 50
    theta = initial_theta
    for n in range(n_iters):
        eta = t0 / (n + t1)
        pre_theta = theta
        # 随机查询一个样本
        ind = np.random.randint(len(X_b))
        gradient = dJ_sgd(theta, X_b[ind], y[ind])
        theta -= eta * gradient
    return theta

上式无法让所有的样本都对算法产生影响，修改梯度下降的函数：

# n_iters 表示要将所有的数据循环几次
def gradient_descent_sgd(X_b, y, initial_theta, n_iters=5, epsilon=1e-8, t0=5, t1=50):
    theta = initial_theta
    length = len(X_b)
    def learning_rate(n):
        return t0 / (n + t1)
    for cur_iter in range(n_iters):
        random_indexes = np.random.permutation(length)
        for i in range(length):
            eta = learning_rate(cur_iter * length + i)
            pre_theta = theta
            gradient = dJ_sgd(theta, X_b[random_indexes[i]], y[random_indexes[i]])
            theta -= eta * gradient
    return theta