代码实现线性回归中的梯度下降法

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

定义数据集

import numpy as np
import matplotlib.pyplot as plt

np.random.seed(666)  # 添加随机种子

x = 2 * np.random.random(size=100)   # 自己定义数据集
y = x * 3. + 4. + np.random.normal(size=100)

X = x.reshape(-1,1)   # 将 x 转换为 100 行 1 列的形式

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

绘制出如图所示的样本点

使用梯度下降法训练

结合上一节

根据公式定义函数

def J(theta,X_b,y):   # 由目标式写出，其中包括 y,theta,X_b
    try:
        return np.sum((y - X_b.dot(theta)) ** 2) / len(X_b)
    except:
        return float('inf')

def dJ(theta,X_b,y):   # 有求导式写出
    res = np.empty(len(theta))   # 对每个 theta 求导，所以开辟 len(theta) 个空间
    res[0] = np.sum(X_b.dot(theta) - y)   #矩阵第 0 行
    for i in range(1,len(theta)):
       res[i] = np.sum((X_b.dot(theta) - y).dot(X_b[:,i]))  #矩阵第 1-最后 行
             # 以 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_iters = 0    # 初值设为0
    
    while i_iters < 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_iters += 1
        
    return theta

定义如下参数并调用 gradient_descent：

X_b = np.hstack([np.ones((len(x),1)),X])   # 构造 X_b 矩阵
initial_theta = np.zeros(X_b.shape[1])   # 定义初始 theta 的个数
eta = 0.01
theta = gradient_descent(X_b,y,initial_theta,eta)

可得

theta
array([4.02145786, 3.00706277])

接下来将上面的方法封装在 LinearRegression 中：

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))   # 对每个 theta 求导，所以开辟 len(theta) 个空间
        res[0] = np.sum(X_b.dot(theta) - y)   #矩阵第 0 行
        for i in range(1,len(theta)):
            res[i] = (X_b.dot(theta) - y).dot(X_b[:,1])  #矩阵第 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
        cur_iter = 0    # 初值设为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])   # 构造 X_b 矩阵
    initial_theta = np.zeros(X_b.shape[1])
    self._theta = gradient_descent(X_b,y_train,initial_theta,eta,n_iters = 1e4,epsilon=1e-8)
    self.intercept_ = self._theta[0]
    self.coef_ = self._theta[1:]

    return self

测试自己封装的线性回归算法

from machine_learning.playML.LinearRegression import LinearRegression

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

可得结果

lin_reg.coef_
array([3.00706277])
lin_reg.intercept_
4.021457858204859

即截距和斜率与我们刚设置的 y 式

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

结果一致。