线性回归中的梯度下降法及其实现

梯度下降法


 

线性回归中使用梯度下降法

上式中将
的值统一表示成了$X_b^{(i)}\Theta$
把目标函数J($\theta$)当成
 

$\frac{1}{2m}$中的2会在求偏导的时候约掉。如果没有$\frac{1}{m}$,则梯度中每一行的元素都会非常大，在具体实践中就会出现问题
 

在线性回归模型中使用梯度下降法

### 在线性回归模型中使用梯度下降法

import numpy as np
import matplotlib.pyplot as plt

np.random.seed(666)
x = 2 * np.random.random(size = 100)#模拟100个样本，每个样本只有一个特征
y = x * 3. + 4. + np.random.normal(size=100)

X = x.reshape(-1,1)#X为100行1列

X.shape
输出：（100，1）

y.shape
输出：（100，）

plt.scatter(x,y)
plt.show()

输出图片：

使用梯度下降法训练

def J(theta,X_b,y):
    try:
        return np.sum((y - X_b.dot(theta)) ** 2) / len(X_b)
    except:
        return float('inf')#如果出现异常的话，返回一个float值，代表损失函数达到了最大

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[:,1])
    return res * 2 / len(X_b)

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_b = np.hstack([np.ones((len(X),1)),X])#X_b是在原来的X基础上添加一列，每一个元素都是1
initial_theta = np.zeros(X_b.shape[1])#设置一个n+1维的向量。该向量中元素的数量是特征数+1，也就是X_b矩阵中的里列数
eta = 0.01

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

theta#计算结果分别代表截距和斜率
输出结果：array([4.02145786, 3.00706277])

//所得到的截距和斜率与一开始的设定差不多，说明梯度下降法成功的训练了模型

封装我们的线性回归算法

from playML.LinearRegression import LinearRegression

lin_reg = LinearRegression()
lin_reg.fit_gd(X,y)
输出：LinearRegression()

lin_reg.coef_
输出：array([3.00706277])

lin_reg.intercept_
输出：4.021457858204859

import numpy as np
from .metrics import r2_score

class LinearRegression:

    def __init__(self):
        """初始化Linear Regression模型"""
        self.coef_ = None
        self.intercept_ = None
        self._theta = None

    def fit_normal(self, X_train, y_train):
        """根据训练数据集X_train, y_train训练Linear Regression模型"""
        assert X_train.shape[0] == y_train.shape[0], \
            "the size of X_train must be equal to the size of y_train"

        X_b = np.hstack([np.ones((len(X_train), 1)), X_train])
        self._theta = np.linalg.inv(X_b.T.dot(X_b)).dot(X_b.T).dot(y_train)

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

        return self

    def fit_gd(self, X_train, y_train, eta=0.01, n_iters=1e4):
        """根据训练数据集X_train, y_train, 使用梯度下降法训练Linear Regression模型"""
        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)
            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
            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.intercept_ = self._theta[0]
        self.coef_ = self._theta[1:]

        return self

    def fit_sgd(self, X_train, y_train, n_iters=5, t0=5, t1=50):
        """根据训练数据集X_train, y_train, 使用梯度下降法训练Linear Regression模型"""
        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.intercept_ = self._theta[0]
        self.coef_ = self._theta[1:]

        return self

    def predict(self, X_predict):
        """给定待预测数据集X_predict，返回表示X_predict的结果向量"""
        assert self.intercept_ is not None and self.coef_ is not None, \
            "must fit before predict!"
        assert X_predict.shape[1] == len(self.coef_), \
            "the feature number of X_predict must be equal to X_train"

        X_b = np.hstack([np.ones((len(X_predict), 1)), X_predict])
        return X_b.dot(self._theta)

    def score(self, X_test, y_test):
        """根据测试数据集 X_test 和 y_test 确定当前模型的准确度"""

        y_predict = self.predict(X_test)
        return r2_score(y_test, y_predict)

    def __repr__(self):
        return "LinearRegression()"

将第0项和其他项进行统一，$X_0^{(i)}$的值恒为1

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 playML.model_selection import train_test_split

X_train,X_test,y_train,y_test = train_test_split(X,y,seed=666)

from playML.LinearRegression import LinearRegression

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

输出：Wall time: 499 ms
0.8129794056212895

## 使用梯度下降法 

lin_reg2 = LinearRegression()
lin_reg2.fit_gd(X_train,y_train)
输出：
D:\Anaconda3\lib\site-packages\numpy\core\fromnumeric.py:87: RuntimeWarning: overflow encountered in reduce
  return ufunc.reduce(obj, axis, dtype, out, **passkwargs)
C:\Users\wanggege\Desktop\MachineLearning\playML\LinearRegression.py:32: RuntimeWarning: overflow encountered in square
  return np.sum((y - X_b.dot(theta)) ** 2) / len(y)
C:\Users\wanggege\Desktop\MachineLearning\playML\LinearRegression.py:53: RuntimeWarning: invalid value encountered in double_scalars
  if (abs(J(theta, X_b, y) - J(last_theta, X_b, y)) < epsilon):
LinearRegression()

lin_reg2.coef_
输出：array([nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan])

X_train[:10,:]
输出：
array([[1.42362e+01, 0.00000e+00, 1.81000e+01, 0.00000e+00, 6.93000e-01,
        6.34300e+00, 1.00000e+02, 1.57410e+00, 2.40000e+01, 6.66000e+02,
        2.02000e+01, 3.96900e+02, 2.03200e+01],
       [3.67822e+00, 0.00000e+00, 1.81000e+01, 0.00000e+00, 7.70000e-01,
        5.36200e+00, 9.62000e+01, 2.10360e+00, 2.40000e+01, 6.66000e+02,
        2.02000e+01, 3.80790e+02, 1.01900e+01],
       [1.04690e-01, 4.00000e+01, 6.41000e+00, 1.00000e+00, 4.47000e-01,
        7.26700e+00, 4.90000e+01, 4.78720e+00, 4.00000e+00, 2.54000e+02,
        1.76000e+01, 3.89250e+02, 6.05000e+00],
       [1.15172e+00, 0.00000e+00, 8.14000e+00, 0.00000e+00, 5.38000e-01,
        5.70100e+00, 9.50000e+01, 3.78720e+00, 4.00000e+00, 3.07000e+02,
        2.10000e+01, 3.58770e+02, 1.83500e+01],
       [6.58800e-02, 0.00000e+00, 2.46000e+00, 0.00000e+00, 4.88000e-01,
        7.76500e+00, 8.33000e+01, 2.74100e+00, 3.00000e+00, 1.93000e+02,
        1.78000e+01, 3.95560e+02, 7.56000e+00],
       [2.49800e-02, 0.00000e+00, 1.89000e+00, 0.00000e+00, 5.18000e-01,
        6.54000e+00, 5.97000e+01, 6.26690e+00, 1.00000e+00, 4.22000e+02,
        1.59000e+01, 3.89960e+02, 8.65000e+00],
       [7.75223e+00, 0.00000e+00, 1.81000e+01, 0.00000e+00, 7.13000e-01,
        6.30100e+00, 8.37000e+01, 2.78310e+00, 2.40000e+01, 6.66000e+02,
        2.02000e+01, 2.72210e+02, 1.62300e+01],
       [9.88430e-01, 0.00000e+00, 8.14000e+00, 0.00000e+00, 5.38000e-01,
        5.81300e+00, 1.00000e+02, 4.09520e+00, 4.00000e+00, 3.07000e+02,
        2.10000e+01, 3.94540e+02, 1.98800e+01],
       [1.14320e-01, 0.00000e+00, 8.56000e+00, 0.00000e+00, 5.20000e-01,
        6.78100e+00, 7.13000e+01, 2.85610e+00, 5.00000e+00, 3.84000e+02,
        2.09000e+01, 3.95580e+02, 7.67000e+00],
       [5.69175e+00, 0.00000e+00, 1.81000e+01, 0.00000e+00, 5.83000e-01,
        6.11400e+00, 7.98000e+01, 3.54590e+00, 2.40000e+01, 6.66000e+02,
        2.02000e+01, 3.92680e+02, 1.49800e+01]])

# 由上述数据可以知道：每一个特征所对应的数据的规模是不一样的，有的只用零点几，有的数据能达到几百的维度，面对这样的数据我们最终所求到的结果的维度可能也是非常大的，默认的eta（0.01）最终形成的步长还是太大，使得我们梯度下降法的过程最终是不收敛的，为验证这个假设，我们手动设置一个eta

lin_reg2.fit_gd(X_train,y_train,eta = 0.000001)
输出：LinearRegression()

lin_reg2.score(X_test,y_test)
输出：0.2758681872447726

### n_iters设置梯度下降法最多要执行几次

%time lin_reg2.fit_gd(X_train,y_train,eta=0.000001,n_iters=1e6)
输出：Wall time: 1min 7s
LinearRegression()

### 运行了1e6次最终得到的结果还是未达到理想中的最小值，如果需要更理想的值则需要执行更多次，但显然太耗时了，更好的额解决方式是进行数据归一化

lin_reg2.score(X_test,y_test)
输出：0.754293258194391

梯度下降法与数据归一化

使用梯度下降法前进行数据归一化

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

standardScaler.fit(X_train)
StandardScaler()


X_train_standard = standardScaler.transform(X_train)


lin_reg3 = LinearRegression()

%time lin_reg3.fit_gd(X_train_standard,y_train)
Wall time: 868 ms
LinearRegression()

X_test_standard = standardScaler.transform(X_test)
该结果与使用正规方程得到的结果是一致的，说明已经找到了损失函数的最小值

lin_reg3.score(X_test_standard,y_test)
0.8129873310487505