线性回归模型的构建与训练

1.基本的导入与配置

# To support both python 2 and python 3
from __future__ import division, print_function, unicode_literals# Common imports
import numpy as np
import pandas as pd
import os# to make this notebook's output stable across runs
np.random.seed(42)# To plot pretty figures
%matplotlib inline
import matplotlib as mpl
import matplotlib.pyplot as plt
mpl.rc('axes', labelsize=14)
mpl.rc('xtick', labelsize=12)
mpl.rc('ytick', labelsize=12)# Where to save the figures
PROJECT_ROOT_DIR = "."
CHAPTER_ID = "training_linear_models"# Ignore useless warnings (see SciPy issue #5998)
import warnings
warnings.filterwarnings(action="ignore", message="^internal gelsd")

2.正规方程处理线性回归

用y=4+3x+高斯噪声来生成100个数据点（x0为1）

# 首先构造X，维度X中所有点的范围在0-2之间，可利用np.random.rand生成符合均匀分布的0-1之间的点。
X = 2 * np.random.rand(100, 1) 
# 构造y, 可用np.random.randn构造高斯噪声
y = 4 + 3 * X + np.random.randn(100, 1)

用matplotlib查看生成的点

#2
plt.plot(X, y, "b.")
plt.xlabel("$x_1$", fontsize=18)
plt.ylabel("$y$", rotation=0, fontsize=18)
plt.axis([0, 2, 0, 15])
# save_fig("generated_data_plot")
plt.show()

通过求解MSE（均方误差）的最小值^θ=(XT⋅X)−1⋅XT⋅y 求解θ值（这就是一个模型的训练过程）

# add x0 = 1 to each instance
X_b = np.c_[np.ones((100, 1)), X]  
# 通过正规方程求解最佳theta值
theta_best = np.linalg.inv(X_b.T.dot(X_b)).dot(X_b.T).dot(y)
theta_best

说明：

X_b = np.c_[np.ones((100, 1)), X]这行代码创建了一个新矩阵X_b，它是原始特征矩阵X的扩展。np.ones((100, 1))生成了一个100行1列的全1矩阵，这在线性代数中常用来表示截距项（也就是当所有特征值为0时的偏置项）。np.c_是NumPy中的一个函数，用于沿着第二轴（列）连接两个数组。所以，X_b实际上是将全1列向量添加到X的每一行前面，以便在模型中包含截距项。

theta_best：

用matplotlib查看拟合的回归函数（的预测结果）

plt.plot(X_new, y_predict, "r-", linewidth=2, label="Predictions")
plt.plot(X, y, "b.")
plt.xlabel("$x_1$", fontsize=18)
plt.ylabel("$y$", rotation=0, fontsize=18)
plt.legend(loc="upper left", fontsize=14)
plt.axis([0, 2, 0, 15])
plt.show()

LinearRegression是基于最小二乘方法的线性模型（训练过程使用正规方程求解）

因此使用LinearRegression求解的theta和y值应该与前面正规方程求出的结果一致。

LinearRegression(

*, # 后面参数都是可选参数

fit_intercept=True, # bool值，指定是否应该计算截距项（即线性方程中的常数项）

copy_X=True, # 指定是否应该复制 X

n_jobs=None, # 指定用于计算的作业数，-1表示使用所有可用的CPU核心，1表示不使用并行计算，None，则作业数由底层实现决定（通常是一个合理的默认值，取决于是否安装了支持并行计算的库，如 joblib

positive=False, # 仅在 fit_intercept=False 时有效，False，则不对系数的符号进行限制，True，则强制系数为正数

)

from sklearn.linear_model import LinearRegression
# Step1：构建LinearRegression()估计器的实例对象lin_reg
lin_reg = LinearRegression()
# Step2:调用实例对象lin_reg的fit方法进行训练
lin_reg.fit(X, y)
# Step2:通过intercept_和coef_查看训练得到的参数值，即线性方程的截断和系数
lin_reg.intercept_, lin_reg.coef_

intercept_属性存储了模型的截距项，coef_属性存储了模型的系数。

说明：

array([4.21509616])这是一个只包含一个元素的数组，表示线性回归模型的截距（intercept）。截距是线性方程中的常数项，也就是当所有特征变量（x值）都为0时，因变量（y值）的预测值。在这个例子中，截距是4.21509616。

array([[2.77011339])这是一个二维数组，包含一个元素的一维数组，表示线性回归模型的系数（coefficients）。系数是线性方程中特征变量的系数，表示每个特征变量对因变量的影响程度。在这个例子中，只有一个特征变量，其系数是2.77011339。

所以，这两个参数定义了一个简单的线性回归方程：

y=4.21509616+2.77011339⋅x

其中，x 是特征变量，y 是因变量。这个方程可以用来预测给定特征变量 x 值的因变量 y 值。

# Step3:调用估计器的predict方法对X_new进行预测
lin_reg.predict(X_new)

第二值是预测的结果y

3.梯度下降训练线性回归

# 学习率
eta = 0.1
n_iterations = 1000# 有100个数据点 
m = 100# 随机初始化参数theta，模型有两个参数（包括截距项），因此theta是一个2x1的向量 
theta = np.random.randn(2,1)for iteration in range(n_iterations):# 计算梯度，这里使用了批量梯度下降（BGD）方法  # 2/m是归一化因子，确保梯度是平均梯度  # X_b.T是扩展特征矩阵的转置，它包含了截距项  # X_b.dot(theta) - y计算了模型预测与真实标签之间的误差  # 整个表达式计算了损失函数关于theta的梯度  
    gradients = 2/m * X_b.T.dot(X_b.dot(theta) - y)# 使用梯度下降法更新参数theta  # eta * gradients计算了参数应该更新的量（方向和大小）  # theta - eta * gradients则是更新后的参数值  
    theta = theta - eta * gradients

解释：

1.初始化参数：

eta = 0.1：这是学习率，它控制着在每次迭代中参数更新的步长。较小的学习率意味着更新步长较小，可能需要更多的迭代次数来收敛；较大的学习率可能导致更新步长过大，从而越过最优解或导致训练不稳定。

n_iterations = 1000：这是迭代次数，即梯度下降算法将执行1000次迭代。

m = 100：这是数据点的数量，表示训练集中有100个样本。

theta = np.random.randn(2,1)：这是参数theta的初始化，theta是一个2x1的向量，包含两个参数（一个是截距项，另一个是特征的系数）。这里使用正态分布随机初始化参数。

2.迭代过程：

for iteration in range(n_iterations)：开始1000次迭代的循环。

gradients = 2/m * X_b.T.dot(X_b.dot(theta) - y)：计算梯度。

X_b.dot(theta)：计算模型预测值，即线性组合X_b（扩展特征矩阵）和theta（参数向量）的点积。

X_b.dot(theta) - y：计算预测值和真实值之间的误差。

X_b.T.dot(...)：计算误差与特征矩阵X_b的点积，得到梯度向量。

2/m * ...：将梯度向量除以样本数量m，得到平均梯度，这是批量梯度下降的关键步骤。

theta = theta - eta * gradients：更新参数theta。

eta * gradients：计算参数更新的量，即学习率eta乘以梯度向量gradients。

theta - ...：根据梯度下降的方向更新参数theta，减去计算出的更新量。

3.结果：

经过1000次迭代后，theta将收敛到最小化损失函数的值，即找到了最佳拟合参数。

# 查看梯度下降训练之后的theta
theta

X_new_b.dot(theta)

绘制不同学习率下的学习过程

# 设置随机种子以保证结果的可重复性  
np.random.seed(42)  
# 随机初始化参数theta（模型有两个参数，包括截距项，所以维度为（2，1））  
theta = np.random.randn(2,1)  # 设置绘图的大小  
plt.figure(figsize=(10,4))  
# 绘制三个子图，分别对应不同的学习率eta  
plt.subplot(131); plot_gradient_descent(theta, eta=0.02)  # 学习率0.02  
plt.ylabel("$y$", rotation=0, fontsize=18)  # 设置y轴标签  
plt.subplot(132); plot_gradient_descent(theta, eta=0.1, theta_path=theta_path_bgd)  # 学习率0.1，并记录theta路径  
plt.subplot(133); plot_gradient_descent(theta, eta=0.5)  # 学习率0.5  # 显示图表  
plt.show()

这里面绘制了前10次迭代的情况，可以看到第一幅图，由于学习率很小，梯度下降进展很慢。第二幅图，可以看到当学习率适中时，梯度下降可以快速开始寻找最优解。第三幅图，当学习率比较大时，梯度下降很可能永远找不到最优解。

所以我们就记录了效果比较好的学习率0.1的theta变化过程，用来与其他方法的theta变化过程进行比较

4.随机梯度下降

# 初始化一个空列表，用于存储随机梯度下降过程中theta的路径
theta_path_sgd = []
# 获取数据点的数量
m = len(X_b)
np.random.seed(42)# 设置训练的轮数（epoch） 
n_epochs = 50
# 定义学习率调度（learning schedule）的超参数 ，这些参数用于调整学习率随时间下降的速度
t0, t1 = 5, 50  # 定义一个函数，用于根据当前的迭代次数计算学习率
def learning_schedule(t):return t0 / (t + t1)# 随机初始化参数theta
theta = np.random.randn(2,1)  # random initializationfor epoch in range(n_epochs):for i in range(m):# 以下4行代码仅用于在第一个epoch的前20次迭代中绘制预测线if epoch == 0 and i < 20:                    # not shown in the book
            y_predict = X_new_b.dot(theta)           # not shown
            style = "b-" if i > 0 else "r--"         # not shown
            plt.plot(X_new, y_predict, style)        # not shown# 随机选择一个数据点的索引 
        random_index = np.random.randint(m)# 获取随机选择的数据点xi及其对应的标签yi
        xi = X_b[random_index:random_index+1]
        yi = y[random_index:random_index+1]# 计算梯度（注意这里由于只考虑了一个数据点，所以梯度是基于单个样本的）
        gradients = 2 * xi.T.dot(xi.dot(theta) - yi)# 根据当前的迭代次数（epoch * m + i）计算学习率 
        eta = learning_schedule(epoch * m + i)# 更新参数theta 
        theta = theta - eta * gradients# 将更新后的theta添加到路径列表中
        theta_path_sgd.append(theta)                 plt.plot(X, y, "b.")                                
plt.xlabel("$x_1$", fontsize=18)                     
plt.ylabel("$y$", rotation=0, fontsize=18)           plt.axis([0, 2, 0, 15])                              
plt.show()  

说明：

初始化和设置：

theta_path_sgd = []：初始化一个空列表，用于存储随机梯度下降过程中theta的路径，即每次更新后的参数值。

m = len(X_b)：获取数据点的数量。

np.random.seed(42)：设置随机数种子，确保结果的可重复性。

n_epochs = 50：设置训练的轮数（epoch），即遍历整个数据集的次数。

t0, t1 = 5, 50：定义学习率调度的超参数，用于调整学习率随时间下降的速度。

学习率调度函数：

def learning_schedule(t):：定义一个函数，根据当前的迭代次数t计算学习率。这里的学习率随着迭代次数的增加而减小，这是一种常见的学习率衰减策略。

参数初始化：

theta = np.random.randn(2,1)：随机初始化参数theta，这是一个2x1的向量，包含两个参数（一个是截距项，另一个是特征的系数）。

训练过程：

for epoch in range(n_epochs):：开始50轮的训练过程。

for i in range(m):：在每个epoch中，遍历每个数据点。

random_index = np.random.randint(m)：随机选择一个数据点的索引，这是SGD中“随机”的由来。

xi = X_b[random_index:random_index+1]：获取随机选择的数据点xi。

yi = y[random_index:random_index+1]：获取对应的标签yi。

gradients = 2 * xi.T.dot(xi.dot(theta) - yi)：计算梯度，这里只针对一个数据点计算。

eta = learning_schedule(epoch * m + i)：根据当前的迭代次数计算学习率。

theta = theta - eta * gradients：更新参数theta。

theta_path_sgd.append(theta)：将更新后的theta添加到路径列表中。

绘图：

plt.plot(X, y, "b.")：绘制原始数据点。

plt.xlabel("$x_1$", fontsize=18)：设置x轴标签。

plt.ylabel("$y$", rotation=0, fontsize=18)：设置y轴标签。

plt.axis([0, 2, 0, 15])：设置坐标轴的范围。

plt.show()：显示图表。

这段代码通过随机梯度下降算法训练线性回归模型，并在训练过程中记录参数的更新路径。学习率随着迭代次数的增加而减小，这有助于模型在训练初期快速收敛，在训练后期进行微调。最终，代码绘制了原始数据点和训练后的预测线。

# 观察经过随机梯度下降之后最后的theta
theta

sklearn下可以使用函数SGDRegressor实现随机梯度下降

SGDRegressor(

loss='squared_error', # 损失函数类型。可选的有 'squared_error', 'huber', 'epsilon_insensitive', 或 'squared_epsilon_insensitive'

*,

penalty='l2', # 惩罚（正则化）类型。可选的有 'none', 'l2', 'l1', 'elasticnet'。有助于防止模型过拟合。

alpha=0.0001, # 正则化强度的倒数（越大表示正则化强度越低）。仅在 penalty 不为 'none' 时使用。

l1_ratio=0.15, # 当 penalty 为 'elasticnet' 时，L1 和 L2 惩罚项之间的比例。仅在 penalty='elasticnet' 时使用

fit_intercept=True, # 是否计算截距项。如果为 False，则假定数据已经中心化

max_iter=1000, # 最大迭代次数。迭代将在达到最大迭代次数或满足其他停止条件时停止

tol=0.001, # 优化的容差。当损失或分数在 tol 内没有改进时，停止训练

shuffle=True, # 是否在每个迭代中随机打乱训练样本，有助于减少收敛到局部最小值的风险。

verbose=0, # 是否打印进度消息到 stdout。0 表示不打印，1 表示偶尔打印，>=2 表示对每个迭代都打印。

epsilon=0.1, # epsilon_insensitive 和 squared_epsilon_insensitive 损失函数的参数。仅当 loss 为这些类型之一时使用。

random_state=None, # 控制伪随机数生成器的种子

learning_rate='invscaling', # 学习率策略。可选的有 'constant', 'optimal', 'invscaling', 'adaptive'。

eta0=0.01, # 初始学习率。仅在 learning_rate='constant' 或 'invscaling' 时使用。

power_t=0.25, # 逆尺度学习率的指数。仅在 learning_rate='invscaling' 时使用。

early_stopping=False, # 是否使用早停来终止训练。如果为 True，则使用验证集上的性能来提前停止训练

validation_fraction=0.1, # 保留作为验证集的数据比例。仅在 early_stopping=True 时使用。

n_iter_no_change=5, # 在早停之前，验证分数在没有改进的情况下需要保持的迭代次数。仅在 early_stopping=True 时使用。

warm_start=False, # 当设置为 True 时，重用上一次调用的解决方案作为初始化，否则，重新分配并拟合一个新的模型

average=False, # 当设置为 True 时，计算所有迭代的平均系数。如果是一个整数，则计算最后 average 次迭代的平均系数。

)

from sklearn.linear_model import SGDRegressor
# Step1：构建SGDRegressor()估计器的实例对象sgd_reg
sgd_reg = SGDRegressor(max_iter=50, tol=None, penalty=None, eta0=0.1, random_state=42)
# Step2:调用实例对象sgd_reg的fit方法进行训练
sgd_reg.fit(X, y.ravel())

# Step2:通过intercept_和coef_查看训练得到的参数值，即线性方程的截断和系数
sgd_reg.intercept_, sgd_reg.coef_

第一个参数是截距b，第二个参数是系数k

5.小批量梯度下降

在迭代的每一步，批量梯度使用整个训练集，随机梯度使用仅仅一个实例，在小批量梯度下降中，它则使用一个随机的小型实例集。它比随机梯度的主要优点在于你可以通过矩阵运算的硬件优化得到一个较好的训练表现，尤其当你使用 GPU 进行运算的时候

# 初始化一个空列表，用于存储小批量梯度下降过程中theta的路径 
theta_path_mgd = []# 设置迭代次数以及每个小批量（minibatch）的大小
n_iterations = 50
minibatch_size = 20np.random.seed(42)
theta = np.random.randn(2,1)  # random initializationt0, t1 = 200, 1000
def learning_schedule(t):return t0 / (t + t1)# 初始化时间步t（用于学习率调度）  
t = 0
for epoch in range(n_iterations):# 在每一个迭代开始时，随机打乱数据点的索引 
    shuffled_indices = np.random.permutation(m)# 根据打乱后的索引，重新排列数据X_b和标签y 
    X_b_shuffled = X_b[shuffled_indices]
    y_shuffled = y[shuffled_indices]# 使用小批量梯度下降，遍历数据（以minibatch_size为步长）  for i in range(0, m, minibatch_size):# 更新时间步t
        t += 1# 获取当前小批量的数据点和对应的标签
        xi = X_b_shuffled[i:i+minibatch_size]
        yi = y_shuffled[i:i+minibatch_size]# 计算当前小批量的梯度（注意这里要除以minibatch_size进行归一化）
        gradients = 2/minibatch_size * xi.T.dot(xi.dot(theta) - yi)# 根据当前的时间步计算学习率 
        eta = learning_schedule(t)# 更新参数theta
        theta = theta - eta * gradients# 将更新后的theta添加到路径列表中  
        theta_path_mgd.append(theta)

# 观察经过小批量梯度下降之后最后的theta
theta

theta在此处是线性方程的两个参数b和k

6.三个方法的训练过程数据

# 转化为np.array数据类型
theta_path_bgd = np.array(theta_path_bgd)
theta_path_sgd = np.array(theta_path_sgd)
theta_path_mgd = np.array(theta_path_mgd)

plt.figure(figsize=(7,4))
plt.plot(theta_path_sgd[:, 0], theta_path_sgd[:, 1], "r-s", linewidth=1, label="Stochastic")
plt.plot(theta_path_mgd[:, 0], theta_path_mgd[:, 1], "g-+", linewidth=2, label="Mini-batch")
plt.plot(theta_path_bgd[:, 0], theta_path_bgd[:, 1], "b-o", linewidth=3, label="Batch")
plt.legend(loc="upper left", fontsize=16)
plt.xlabel(r"$\theta_0$", fontsize=20)
plt.ylabel(r"$\theta_1$   ", fontsize=20, rotation=0)
plt.axis([2.5, 4.5, 2.3, 3.9])
plt.show()

小结：全量梯度下降法的收敛效果是最好的，因为它将所有的情况都进行了考虑；随机梯度下降法的效果从图中也可以看出‘随机’的特点，收敛不稳定，可能并不能得到比较好的解；小批量梯度下降法则是两者的折中，在减少计算资源消耗的同时也在一定程度上确保了模型的收敛。