线性回归 + 基础优化算法

线性回归

线性回归是机器学习最基础的模型，也是理解后续所有深度学习的基础。


线性模型可以看做是单层神经网络。

上述有个0.5是在求导的时候可以很方便的将2消去。
实际上，这里的数据样本受限很大，比如地球上房子就那么多，肯定不会上百亿的，那么这样具体如何处理，后续会有对应的算法对这部分进行处理，此处不做深入研究，因为笔者还不会（绷）。

因为上述是一个线性模型，因此是有显式解的。

这是一个唯一有最优解的模型，后续较为复杂的模型都不会有最优解了。

下面总结一下：

基础优化算法

优化方法里面最常见的算法叫做梯度下降，即当我的模型没有显示解的时候，需要挑选一个参数的随机初始值，随便选择，记为$w_0$，在接下来的时刻，不断的更新$w_0$，使其接近我们的最优解。具体解释如下：
上述是一个很简单的二次函数的等高图，最优点在最小的圈里，假设$w_0$随机取在图中的位置，去梯度的反方向，和学习率相乘，然后相加，就可以得到$w_1$的位置，同样，在$w_1$处，继续计算其梯度，如此反复迭代，知道到达最优点附近。

上述的学习率是一个超参数，是需要人为指定的值，对于该值，不能太小，要是太小的话，每次走的步长很有限，到达一个点需要走很多很多步，这样需要计算很多梯度，浪费时间和资源，同时也不能太大，要是迈过了，会使得结果一直在震荡，没有真正的在下降。

从 0 开始实现

虽然现在的深度学习框架可以自动化的进行上述线性回归的实现，但是我们还是需要明白线性回归到底是回归啥，因此我们这里选择从零开始实现。
首先导入必要的环境：

%matplotlib inline
import random
import torch
from d2l import torch as d2l

上述的 %matplotlib inline 是一个魔法命令（magic command），主要用于IPython环境和Jupyter Notebook中。它的主要作用是将matplotlib生成的图表直接嵌入到Notebook页面中，而不是在一个弹出的新窗口中显示这些图表。这对于数据分析和可视化工作特别方便，因为它允许用户直接在Notebook中查看图表结果，而无需切换窗口或上下文。

下面我们将根据带有噪声的线性模型构造一个人造数据集（好处是可以知道真实的$w$和$b$.）：

def synthetic_data(w, b, num_examples):  #@save"""生成y=Xw+b+噪声"""X = torch.normal(0, 1, (num_examples, len(w)))# 生成一个均值为0，方差为1的随机数，后面的两个参数代表的是样本数量和样本长度y = torch.matmul(X, w) + b # torch.matmul()用于计算矩阵乘法。# 给y加上一些均值为0，标准差为0.01的噪声，模拟真实世界数据中的误差。# torch.normal()生成与y相同形状的噪声张量，并将其加到y上。y += torch.normal(0, 0.01, y.shape)# 最后，将y转换成列向量（即形状为(num_examples, 1)的张量）。return X, y.reshape((-1, 1))# 然后就是将真实的数据传入到对应函数的位置即可
true_w = torch.tensor([2, -3.4])
true_b = 4.2
features, labels = synthetic_data(true_w, true_b, 1000)

执行上述代码后，features中的每一行都包含一个二维数据样本，labels中的每一行都包含一维标签值（一个标量）。

通过生成第二个特征features[:, 1]和labels的散点图，可以直观观察到两者之间的线性关系。

d2l.set_figsize()
d2l.plt.scatter(features[:, (1)].detach().numpy(), labels.detach().numpy(), 1);
# detach用于从计算图中分离出来一个张量，使得该张量不会被追踪计算历史

下面，我们定义一个data_iter函数，该函数接收批量大小、特征矩阵和标签向量作为输入，生成大小为batch_size的小批量。每个小批量包含一组特征和标签。

def data_iter(batch_size, features, labels):num_examples = len(features)indices = list(range(num_examples))# 这些样本是随机读取的，没有特定的顺序，shuffle用于打乱数据顺序random.shuffle(indices)# 从0开始到num_examples每一次跳跃batch_size 的大小for i in range(0, num_examples, batch_size):batch_indices = torch.tensor(indices[i: min(i + batch_size, num_examples)])# 后面使用min是为了防止超出范围yield features[batch_indices], labels[batch_indices]

下面设置 batch_size = 10 ， 读取一下数据集，并展示出来：

下面初始化模型的权重与偏置：

w = torch.normal(0, 0.01, size=(2,1), requires_grad=True)
# 创建一个二维张量w，形状为2×1，并启用了自动求导的功能
b = torch.zeros(1, requires_grad=True)

下面定义模型，将模型的输入和参数同模型的输出关联起来：

def linreg(X, w, b):  #@save"""线性回归模型"""return torch.matmul(X, w) + b # 触发广播机制

下面定义损失函数：

def squared_loss(y_hat, y):  #@save"""均方损失"""return (y_hat - y.reshape(y_hat.shape)) ** 2 / 2

下面定义优化算法：
尽管线性回归有显式解，但是除了这个模型，其他的模型没有，因此此处选择使用小批量随机梯度下降。

def sgd(params, lr, batch_size):  #@save"""小批量随机梯度下降"""with torch.no_grad():#更新的时候不需要参与梯度计算for param in params:param -= lr * param.grad / batch_sizeparam.grad.zero_()# 手动将梯度设置成0，这样下一次计算梯度的时候就不会和上一次相关了

训练

现在我们已经准备好了模型训练所有需要的要素，可以实现主要的训练过程部分了。
理解这段代码至关重要，因为从事深度学习后，相同的训练过程几乎一遍又一遍地出现。在每次迭代中，我们读取一小批量训练样本，并通过我们的模型来获得一组预测。计算完损失后，我们开始反向传播，存储每个参数的梯度。最后，我们调用优化算法 sgd 来更新模型参数。
概括一下，我们将执行以下循环：

• 初始化参数
• 重复以下训练，直到完成
  • 计算梯度$\mathbf{g}\leftarrow {\partial }_{(\mathbf{w},b)}\frac{1}{|\mathcal{B}|}{\sum }_{i\in \mathcal{B}}l({\mathbf{x}}^{(i)},y^{(i)},\mathbf{w},b)$
  • 更新参数$(\mathbf{w},b)\leftarrow (\mathbf{w},b)-\eta \mathbf{g}$

在每个 迭代周期（epoch）中，我们使用data_iter函数遍历整个数据集，并将训练数据集中所有样本都使用一次（假设样本数能够被批量大小整除）。这里的迭代周期个数num_epochs和学习率lr都是超参数，分别设为3和0.03。设置超参数很棘手，需要通过反复试验进行调整。

lr = 0.03
num_epochs = 3 # 把整个数据扫三遍
net = linreg
loss = squared_loss  # 均方损失for epoch in range(num_epochs):for X, y in data_iter(batch_size, features, labels):l = loss(net(X, w, b), y)  # X和y的小批量损失# 因为l形状是(batch_size,1)，而不是一个标量。l中的所有元素被加到一起，# 并以此计算关于[w,b]的梯度l.sum().backward()sgd([w, b], lr, batch_size)  # 使用参数的梯度更新参数with torch.no_grad():train_l = loss(net(features, w, b), labels)print(f'epoch {epoch + 1}, loss {float(train_l.mean()):f}')

因为我们使用的是自己合成的数据集，所以我们知道真正的参数是什么。因此，我们可以通过比较真实参数和通过训练学到的参数来评估训练的成功程度。事实上，真实参数和通过训练学到的参数确实非常接近。

注意，我们不应该想当然地认为我们能够完美地求解参数。在机器学习中，我们通常不太关心恢复真正的参数，而更关心如何高度准确预测参数。幸运的是，即使是在复杂的优化问题上，随机梯度下降通常也能找到非常好的解。其中一个原因是，在深度网络中存在许多参数组合能够实现高度精确的预测。

简洁实现

在上述的实现中，我们只运用了：

• （1）通过张量来进行数据存储和线性代数；
• （2）通过自动微分来计算梯度。

实际上，由于数据迭代器、损失函数、优化器和神经网络层很常用，现代深度学习库也为我们实现了这些组件。

import numpy as np
import torch
from torch.utils import data
from d2l import torch as d2l
# 生成数据集
true_w = torch.tensor([2, -3.4])
true_b = 4.2
features, labels = d2l.synthetic_data(true_w, true_b, 1000)
# 读取数据集
def load_array(data_arrays, batch_size, is_train=True):  #@save"""构造一个PyTorch数据迭代器"""dataset = data.TensorDataset(*data_arrays)# 将多个数据数组打包成一个TensorDataset# *data_arrays 会将列表中的元素解包，作为TensorDataset的参数传入return data.DataLoader(dataset, batch_size, shuffle=is_train)# 创建一个DataLoader实例，用于加载数据# 如果is_train为True，则在每个epoch开始时会对数据进行随机洗牌batch_size = 10
data_iter = load_array((features, labels), batch_size)next(iter(data_iter))# 定义模型
from torch import nnnet = nn.Sequential(nn.Linear(2, 1))net[0].weight.data.normal_(0, 0.01) #权重还是随机采样
net[0].bias.data.fill_(0)loss = nn.MSELoss()trainer = torch.optim.SGD(net.parameters(), lr=0.03)num_epochs = 3# 设置训练次数
for epoch in range(num_epochs):for X, y in data_iter:l = loss(net(X) ,y) # 前向传播计算损失trainer.zero_grad()# 在每次计算梯度之前，需要先将之前的梯度清零，因为默认情况下，.backward() 方法会累积梯度。l.backward()trainer.step()# 根据当前的梯度和优化器的设置来更新网络的参数。l = loss(net(features), labels) # 评估当前轮次的损失print(f'epoch {epoch + 1}, loss {l:f}')w = net[0].weight.data
print('w的估计误差：', true_w - w.reshape(true_w.shape)) # 需要改变形状
b = net[0].bias.data
print('b的估计误差：', true_b - b) # b就是一个标量，无需改变形状

思考QA

Q1 ：为啥使用平方损失而不是绝对差值？
A1：实际上没有啥区别，最开始选择使用平方差值的原因是因为绝对差值在原点处不可导，导数可能处理起来有点困难，但是实际上两者都是可以的。

Q2：损失为啥需要求平均？
A2：实际上求与不求都是没有关系的，除以 N 的好处是不管我的样本有多大，求出来的梯度的值都是差不多的，方便于我调整学习率。要是没有除以 N 的话，在学习率上除也是一样的。除以 N 使得 Loss 与样本数无关。

Q3：batchsize是否会最终影响模型结果？batchsize过小是否可能导致最终累计的梯度计算不准确？
A3：实际上在同样的 epoch 下 batch_size 越小的话其实是对收敛越好的。一定的噪音使得我的神经网络不会走偏，可以减少过拟合，因为要是 batch_size 越小的话，反应出来的噪音是越大的。

Q4：随机梯度下降的“随机”是啥随机啊？
A4：指的是样本随机取，也就是代码中的 shuffle = True

Q5：为什么机器学习优化算法都采取梯度下降（一阶导算法），而不采用牛顿法（二阶导算湿），收敛速度更快，一般能算出一阶导，二阶导也应该能算。
A5：求解复杂，收敛不一定快，就算快，结果不一定好。

Q6：这样的data-iter写法，每次都把所有输入load进去，如果数据多的话，最后内存会爆掉吧？有什么好的办法吗？
A6：确实有这个问题，但是这里只是展示，实际上数据是放在硬盘上的，每次读取部分数据进去，而不是全部读入内存。

Q7：如果样本大小不是批量数的整数倍，那需要随机删除多余的样本数吗？
A7：有三种做法：

• 最常见的做法就是拿到一个小一点的样本，也就是最后剩下多少样本都拿来
• 丢弃最后不完整的样本数
• 从下一个 epoch 中补其所需要的数量使其满足批量数

Q8：如何进行收敛的判断？
A8：有以下几种方法：

• 比较两个 epoch 之间，要是目标函数变化不大的时候，就可以认为是收敛的（比如相对变化是0.01就停止迭代）
• 拿一个验证数据集，要是验证数据集上的精度没增加的话，就可以认为是收敛的。

Q9：每个batch计算的时候，为什么需要把剃度先清零呀?
A9：因为 PyTorch 不会自动清零，要是不清零的话，后续求出来的梯度会在原来已有梯度的基础上进行累加，这样会导致出现意想不到的错误。

本质上我们为什么要用SGD，是因为大部分的实际Ioss太复杂，推导不出导数为0的解？只能逐个batch去逼近。

后记

考虑到我的Pycharm上没有d2l，因此我根据上述代码编写了如下可以运行的代码，也是从0开始实现的线性回归：

import random
import torch
import matplotlib.pyplot as plt# 生成数据集(weight, bias, data_num)
def synthetic_data(w, b, num_examples):# X 的形状为 (num_examples, len(w))   w 的形状为 (len(w), 1)X = torch.normal(0, 1, (num_examples, len(w)))y = torch.matmul(X, w) + by += torch.normal(0, 0.01, y.shape)  # 加噪return X, y.reshape((-1, 1))  # 将y转换为列向量# 读取数据集
def data_iter(batch_size, features, labels):num_examples = len(features)  # 返回特征矩阵中样本的数量（即第一维的大小）indices = list(range(num_examples))  # 创建所有样本的索引random.shuffle(indices)  # 打乱顺序for i in range(0, num_examples, batch_size):# 这里的 indices 不一定需要转化为 tensor ，直接使用列表也是可以的batch_indices = torch.tensor(indices[i:min(i + batch_size, num_examples)]  # 确保不会超出样本总数)yield features[batch_indices], labels[batch_indices]# 下一次调用时，函数会从上次暂停的地方继续执行，直到所有批次都被处理完# 初始化模型参数
def initialize_params(input_dim):# w 为一个二维张量， 形状为 input_num × 1w = torch.normal(0, 0.01, size=(input_dim, 1), requires_grad=True)b = torch.zeros(1, requires_grad=True)return w, b# 定义模型
def linreg(X, w, b):return torch.matmul(X, w) + b# 定义损失函数
def squared_loss(y_hat, y):return (y_hat - y.reshape(y_hat.shape)) ** 2 / 2# 定义优化算法
# params: 一个包含所有需要优化的模型参数的列表,每个参数是一个张量，且能自动计算梯度
def sgd(params, lr, batch_size):with torch.no_grad():for param in params:param -= lr * param.grad / batch_size# 由于梯度通常是对一个小批次数据计算的，因此需要对梯度进行归一化param.grad.zero_()  # 梯度清零# 绘制散点图
def plot_scatter(features, labels):# features[:,1} 提取特征列 ，s 设置散点的大小plt.scatter(features[:, 1].detach().numpy(), labels.detach().numpy(), s=1)plt.xlabel("Feature 1")plt.ylabel("Labels")plt.title("Scatter Plot of Feature 1 vs Labels")plt.show()# 训练过程
def train(features, labels, batch_size, lr, num_epochs, w, b):net = linregloss = squared_lossfor epoch in range(num_epochs):for X, y in data_iter(batch_size, features, labels):l = loss(net(X, w, b), y)  # 计算小批量损失l.sum().backward()  # 将所有样本的损失相加，并反向传播sgd([w, b], lr, batch_size)  # 更新参数 w,bwith torch.no_grad():train_l = loss(net(features, w, b), labels)  # 更新后的参数和正确参数进行比较print(f'epoch {epoch + 1}, loss {float(train_l.mean()):f}')  # 打印每次迭代后的训练损失值if __name__ == '__main__':# 超参数true_w = torch.tensor([2, -3.4], dtype=torch.float32)true_b = 4.2num_examples = 1000batch_size = 10lr = 0.03num_epochs = 3# 生成数据集features, labels = synthetic_data(true_w, true_b, num_examples)# 打印第一个样本print('features:{}, labels:{}'.format(features[0], labels[0]))# 绘制散点图plot_scatter(features, labels)# 初始化参数w, b = initialize_params(len(true_w))# 训练模型train(features, labels, batch_size, lr, num_epochs, w, b)# 输出误差print(f'w的估计误差: {true_w - w.reshape(true_w.shape)}')print(f'b的估计误差: {true_b - b}')