【笔记】随机梯度下降法

随机梯度下降法

上文书我们说到批量梯度下降法

批量梯度下降法是有个问题的，如果样本量非常的大的话，计算梯度本身就是非常的耗时的，那么应该如何改进呢，基于能不能每次只对其中的一个值进行计算的思想，我们可以变幻出新的结构

对于i来说，我们每次都随机的只取一个固定的i，相应的，括号外也就不需要除以m了，同时可以进行向量化（对于X_b来说每次只取一行就行了）

那么我们设想一下，每一次随机取一个i，对于随机出的i，我们都对其进行这个式子的计算，由于这个式子也是一个向量，也可以表达一个方向，我们向着这个方向不停的搜索迭代，最后得到损失函数的最小值，那么基于这样的思想的实现的思路就是随机梯度下降法，具有不可预知性，即便是会不断地变化方向，到时最后还是可以来到最小值的附近，虽然精度可能不高，但是某些时候需要用到部分的精度来换取部分的时间，这个时候随机梯度下降法就有意义了

且使用的时候有一个很重要的技巧，即在下降的过程中，学习率的取值是很重要的，若学习率不好且为一个固定值，那么我们可能会错过最小值的位置，那么我们可以将学习率设置成一个递减的，随着循环的增加，学习率变得越来越小，这样，使用倒数的形式是可以的

但是这也是有问题的，如果学习率下降的比率差距太大也是不行的，所以一般是实现的时候，我们会在分母上添加一个常数来缓解这种情况，有的时候我们的分子一直取1也是很难达到我们想要的效果，这里我们就将分子也替换成一个常数，这两个就是对应的超参数（这里设置为a为5，b为50）

实际上这种逐渐递减的思想就是模拟的模拟退火的思想，其就是模拟打铁的时候温度是要从高到低的冷却的退火的过程，冷却的函数是和时间相关的

实现

（在notebook中）

这里一样的使用了虚拟的数据，m为十万，特征数为1，同样的随机出一个十万个元素的向量，再将其转换成m*1的矩阵，之后在生成一个y

  m = 100000
  x = np.random.normal(size=m)
  X = x.reshape(-1,1)
  y = 4. *x + 3. + np.random.normal(0,3,size=m)

同样的，虽然只有一个特征，但是后面是适用于多个特征的情况的，使用先前的拟合出来的函数，然后生成X对应的X_b，然后设置一个初始化为0的一个theta值，eta设为0.01，使用梯度下降法，寻找最好的theta值

  %%time
  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)

运行时间可以得到

输出后theta为

具体使用随机梯度下降法

注意的是，与之前不一样的是，之前传进来的是一整个矩阵，这个里面是一行的数据，y也一样只是其中的一个数值，然后对先前分析的计算方式进行实现

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

搜索过程，将常数设置为5和50，eta直接进行计算，再循环前我们还是让theta为 initial_theta，之后我们要进行循环的过程，只要判断当前的循环次数有没有达到我们之前的预设值就好的，在循环中，首先我们随机一个样本，获得这个样本的索引以后，相应的我们就可以提取出来我们所对应的这个值送去进行随机梯度的计算，之后计算theta为当前的theta减去学习率再乘上gradient，最后返回这样就完成了随机梯度下降法

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

      t0 = 5
      t1 = 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

我们测试一下随机梯度下降法，为了体现出威力，我们让循环总数为样本总数除以3，即我们只随机的检查了三分之一个样本总量，真实使用的时候是不会这样的

  %%time
  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=len(X_b)//3)

最终效果为

我们使用theta来看一下值

可以看出来，这两个的结果是非常接近的

我们对随机梯度下降法进行封装

我们设其为fit_sgd，和上面的思想相同，指定默认的循环次数为一万，设置好两个常数，其余的和上述大致相同

那么这样实现会出现什么问题呢，因为我们在使用的时候是要将所有的样本至少遍历一遍的，这里就会使得这个循环的个数很难取值，所以在这种情况下，n_iters是代表对我们的样本看几圈，所以其实这个的初始化的值可以非常的小，设置成5也没问题，这样也是有点问题的

那么我们直接做两重循环，外重循环就是对n_iters进行一次循环，那么在循环的过程中，我们都要对所有的样本看一遍，那么我们就要对原始样本进行一个乱序的排序，这样才能保证对全部的样本都有循环且还是随机的，保证了随机性也保证了我们每一个样本至少都被参照过了，那么我们对索引进行乱序，然后设置一个X_b_new来代表乱序化以后的结果，同样y也是如此，那么我们就开始对新的X和y进行遍历，这样我们就不需要再随机化索引了，直接取就行了，针对这个梯度再传theta就行了

代码如下

def fit_sgd(self, X_train, y_train, n_iters=5, t0=5, t1=50):
    """根据训练数据集使用梯度算法训练模型"""
    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, t0=5, t1=50):

        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):
                gradient = dJ_sgd(theta, X_b_new[i], y_new[i])
                theta = theta - learning_rate(cur_iter * m + i) * gradient

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

    return self

那么我们就使用自己的SGD来看一下实际的效果是怎么样的

（在notebook中）

应用在虚拟的例子中来验证算法的正确性

  m = 100000
  x = np.random.normal(size=m)
  X = x.reshape(-1,1)
  y = 4. *x + 3. + np.random.normal(0,3,size=m)

使用我们自己的算法，实例化对象，调用fit_sgd进行训练，看两次

  from LinearRegression import LinearRegression
  lin_reg = LinearRegression()
  lin_reg.fit_sgd(X,y,n_iters=2)

结果如下

对应的系数为lin_reg.coef_

对应的截距lin_reg.interception_

使用真实的数据（波士顿房价）来验证，调用以后在设置范围，最后分割数据，随机种子666

  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 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 LinearRegression import LinearRegression

  lin_reg = LinearRegression()
  %time lin_reg.fit_sgd(X_train_standard,y_train,n_iters=2)
  lin_reg.score(X_test_standard,y_test)

结果如下

很显然，没有达到最好的结果，我们使用50次

  %time lin_reg.fit_sgd(X_train_standard,y_train,n_iters=50)
  lin_reg.score(X_test_standard,y_test)

结果如下

可以发现，结果更好了

对于sklearn来说，也有这种随机梯度下降法，和我们自己写的是不一样的，其使用了很多的优化方案，也更加的复杂，需要注意的是，这个只能解决线性模型，同样进行实例化，进行fit，在训练之后看一下最终结果的准确度

  from sklearn.linear_model import SGDRegressor
  sgd_reg = SGDRegressor()
  %time sgd_reg.fit(X_train_standard,y_train)
  sgd_reg.score(X_test_standard,y_test)

结果如下

相应的我们传入一个参数，我们设为100

  sgd_reg = SGDRegressor(n_iter_no_change=100)
  %time sgd_reg.fit(X_train_standard,y_train)
  sgd_reg.score(X_test_standard,y_test)

结果如下（为啥还变低了？）