【笔记】线性回归中的梯度下降法（实现以及向量化并进行数据归一化）

线性回归中的梯度下降法（实现以及向量化并进行数据归一化）

多元线性回归中的梯度下降法

我们试一下应用在多元线性回归中，对于线性回归的问题，前面的基本都是使每一次模型预测出的结果和数据所对应的真值的差的平方的和为损失函数，对于参数来说，有n+1个元素，这种情况下，我们就需要变换式子

这实际上就是求对应的梯度值，梯度本身也是一个向量，含有n+1个元素，即对每一个参数进行一次偏导数，这样一来，梯度就代表了方向，对应J增大最快的方向，区别在于这里变成了向量操作

原理终了，这样我们的目标就变成了使下面的式子尽可能的小

梯度就是J对每一个维度上的参数进行求导，导数怎么求呢，我们除了对截距求导的时候是要在其基础上乘上-1以外

其余的n个参数都是按照链式求导来计算

我们将其排成一个列向量，最后整理就可以得到

很明显的是，我们不难发现其中每一项都是若干数量的样本进行计算以后的求和结果，这就说明了样本数量是影响着梯度的的，也就是说样本数量越大，求出来的梯度中，元素相应的也就越大，这是不太合理的，应该和m无关才比较好，这样我们再为其除以m，这样我们整体的式子就变成了

这很明显，J实际上就是MSE的值
有的时候会对J取不一样的公式

实际上出来的式子的结果差别不大，但是如果我们没有1/m的话，我们数据过大就会导致一些问题，这就说明了，当我们使用梯度下降法来进行求函数最小值，有时候我们需要对目标函数进行特殊的设计，虽然理论上没问题，但是实际上数据过大可能影响效率，这种一般被称为批量梯度下降法

实现：

（在notebook中）

首先使用我们自己设计的虚假的例子来实现梯度下降法的方式

  import numpy as np
  import matplotlib.pyplot as plt

生成一组一维数组x，size为100，对于y，使用线性的方式生成，在系数截距的上加上一个噪音

  x = 2 * np.random.random(size=100)
  y = x * 3. + 4. + np.random.normal(size=100)

为了代码容易的扩展到多维，用X来代替x，设置成一百行一列这个结构，然后看一下shape

  X = x.reshape(-1,1)
  X.shape
  y.shape

输出如下

然后简单的绘制出图像

  plt.scatter(x,y)

图像如下

为了保证随机性的可重复性，可以在随机生成前用np.random.seed（）的方式设置一个种子，我这里就没有设置出来

下面使用梯度下降法训练

首先我们先设置一个计算损失函数J的值，我们要存进去theta，X_b这个矩阵以及y，相应的直接返回真值与预测值的差值的平方再除以样本数这个结果，设置一个异常检测来防止数据过大导致出现问题

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

然后我们来设置dJ来作为计算其导数的函数，参数与J一样，求导方面，设置最后导数为res，其空间为theta的量，然后我们对第一个元素进行特殊处理，对X_b.dot(theta) - y在进行求和，剩下的n项使用循环就可以得出，对X_b.dot(theta) - y和X_b[:,i]每一个样本取出第i个特征的这两个向量进行点乘，然后返回乘上2/m的结果

  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)

尽管例子是简单的线性回归，但是这个方法对于多元来说也是通用的，然后我们就可以开始梯度下降的过程，我们需要传进X_b和y，相应的计算也是要更改一下，其他的和之前一样不变，相应的，将之前设置用来跟踪theta的删除掉即可，最后返回theta值

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

      theta = initial_theta
      i_iter = 0

      while i_iter < 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_iter += 1

      return theta

下面我们就可以开始调用了，我们先为原先的X添加一列，每一个元素都为1，相应的y值不变，我们需要一个初始化theta，将initial_theta设置成元素数量为特征数的数量加1的向量，元素数量即为X_b的列数，对eta取0.01

  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值是多少
结果如下，其对应的是截距和斜率，虽然有些差距，但是可以发现还是成功地训练了

我们在LinearRegression中新加一个函数fit_gd，将上面的代码封装进取，包括了计算损失函数的函数，计算梯度的函数以及梯度下降法的函数

      def fit_gd(self, X_train, y_train, eta=0.01, n_iters=1e4):
          """根据训练数据集使用梯度算法训练模型"""
          assert X_train.shape[0] == y_train.shape[0], 
              "the size of X_train must be equal to the size of y_train"

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

          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)


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

              theta = initial_theta
              cur_iter = 0

              while cur_iter < 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
  
                  cur_iter += 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, n_iters)

          self.interception_ = self._theta[0]
          self.coef_ = self._theta[1:]

          return self

然后我们试一下使用我们自己的代码

  from LinearRegression import LinearRegression

实例化对象并调用一下梯度下降法进行训练

  lin_reg = LinearRegression()
  lin_reg.fit_gd(X,y)

使用lin_reg.coef_可以看到此时系数结果

通过lin_reg.intercept_可以看到此时截距结果

可以发现和之前是一样的

梯度下降法的向量化和数据标准化

向量化的处理主要集中在梯度的过程中，我们之前是一项一项的将对应的元素给求出来，那么我们尝试将其转换成矩阵的运算，这样算起来就会简单很多，我们对第一项乘上一个恒等于1的X0，这样J就会变成

那么我们开始进行向量化，我们将这看作为两个向量进行点乘的形式，我们可以对梯度里面的式子理解成一个矩阵的运算，我们将前一个向量给单独提出来并将m项给拆开

我们可以将后面的向量写成这个样式，其就是X_b这个矩阵

这就相当于是我们对这个矩阵中每一个对应的列进行点乘，其实际上就是原来式子中的计算方式，即对第0列进行点乘，得到的就是以前的式子的第0项，以此类推，知道点乘到第n列，这样我们就有了n+1项，这样写其实就是将前一个向量写作成一个1m的向量，矩阵则是mn+1的矩阵

将求之前的梯度过程转换成了两个矩阵的乘法运算，最终就可以得到一个1*n+1的向量，而我们原来要求的梯度中也有n+1个元素，这两个中的元素是相等一一对应的，这样我们就得到了新的式子（由于之后使用的时候都是列向量，但是我们上面的分析中其是一个行向量，所以要进行一下转置）

虽然在numpy中的计算是不区分行列向量的，但是为了直观以及好看，我们将其也进行一下转换，对最后的结果整体进行一下转置，很好理解，不再赘述，新的式子即为

这个列向量就是我们要求的梯度

这里我们只要将对梯度的求解部分dJ进行一下改变，将其变成

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

那么我们使用真实的数据进行训练（老波士顿房价了）

（notebook中）

老样子引用加限制范围再进行分割，分割为测试数据集以及训练数据集，随机种子为666

  import numpy as np
  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 LinearRegression import LinearRegression

  lin_reg1 = LinearRegression()
  %time lin_reg1.fit_normal(X_train,y_train)
  lin_reg1.score(X_test,y_test)

结果如下

然后我们使用梯度下降法，首先先进行实例化，调用fit_gd并将训练集传进来

  lin_reg2 = LinearRegression()
  lin_reg2.fit_gd(X_train,y_train)

可以发现，又出现了一些报错

我们使用lin_reg2.coef_来看一下系数，可以发现结果都是nan

原因是我们将其放在了一个真实的数据集中，并没有观察这个数据的情况就进行了操作

我们可以使用X_train[:10,:]来看看前十行的数据

我们可以发现，在其中每一个特征的数据规模是不一样的，面对这样的一个数据，我们最后求出来的结果很有可能也是非常大的，我们使用的默认的eta为0.01，其出来的结果还是发散的，因此我们降低eta，设置eta为0.000001

  lin_reg2.fit_gd(X_train,y_train,eta=0.000001)

输出结果可以看出来，此时是没有报错的

我们可以看看这个情况下的值是多少

  lin_reg2.score(X_test,y_test)

结果为

通过这个值我们可以发现，我们还没有达到这个的最小值，很有可能是因为现在的eta太小了，导致行进的步长很小，那么我们使用更多的循环次数可能就可以了
我们设置成循环次数为一百万次

%time lin_reg2.fit_gd(X_train,y_train,eta=0.000001,n_iters=1e6)

结果为

此时我们用lin_reg2.score(X_test,y_test)来看一下这个最后的值是多少

很显然还是没有达到对应的最小值，虽然我们可以使用更多的循环来使其达到最小值，但是很显然，这太耗时了，因此我们在使用梯度下降前，使用数据归一化就可以很好地解决这个情况

数据归一化

先加载上StandardScaler这个类，然后就是实例化加fit操作就可以了

  from sklearn.preprocessing import StandardScaler
  standardScaler = StandardScaler()
  standardScaler.fit(X_train)

可以得到

然后我们就可以得到标准化以后的结果X_train_standard

  X_train_standard = standardScaler.transform(X_train)

这样我们就对原始数据进行了归一化处理，同样的进行fit，使用默认值并进行计时

  lin_reg3 = LinearRegression()
  %time lin_reg3.fit_gd(X_train_standard,y_train)

得出的输出结果

可以看出来这是非常快就训练完成了，由于我们之前已经进行了数据归一化，同样的也要对test进行数据归一化

  X_test_standard = standardScaler.transform(X_test)
  lin_reg3.score(X_test_standard,y_test)

得到的输出结果为

这是我们已经找到了这个损失函数的最小值，同时这个速度是非常快的

我们通过设计一个虚拟的测试用例来说明梯度下降法的优势在哪里

我们设置样本数为1000个，相应的特征值为5000，为什么样本数会少于特征值呢，这是因为我们在计算梯度的时候，其实是需要所有的样本都来参与计算的，当样本量比较大的时候，我们计算也是会比较慢的，后面可以进行随机梯度来改进，不过这里就不变了

我们设置一个使用随机化的正态分布的基于这个规模的量，这里由于之前归一化过，这里就不需要归一化了，在设置一个初始化的随机生成的theta，然后我们就可以生成一个y了，将矩阵点乘上系数再加上true_theta第0位元素再加上一个噪音即可得到

  m =1000
  n = 5000

  big_X = np.random.normal(size=(m,n))
  true_theta = np.random.uniform(0.0,100.0,size=n+1)
  big_y = big_X.dot(true_theta[1:]) + true_theta[0] + np.random.normal(0.,10.,size=m)

基于上述的数据，我们先实例化后在使用正规方程来解析

  big_reg1 = LinearRegression()
  %time big_reg1.fit_normal(big_X,big_y)

结果如下

然后我们再实例化一个对象，使用梯度下降法来解析

  big_reg2 = LinearRegression()
  %time big_reg2.fit_gd(big_X,big_y)

结果为

可以看出来，我们使用梯度下降法的方式来训练是比正规方程的耗时要短的，加大数据集的话，优势会更加的明显