深度学习之线性回归

一、概念

首先咱们先来了解一下 回归、线性等概念。
回归
回归（regression）是能为一个或多个自变量与因变量之间关系建模的一类方法。直白点，就是先一堆点里面找一根线，让线在尽可能在所有点的中间。
而回归的目的就是为了预测未来，通过找到的这根线来预测未来。
线性
而线性也比较简单，就是自变量、因变量可以用一条直线来表示。
此时应该已经明白线性回归的意思了吧。在一堆点里面寻找一条直线，尽可能让直线在所有点的中间。

二、原理

为了能够有预测某种东西的模型，我们需要收集一个真实的数据集，这个用机器学习的术语叫做 训练数据集 或 训练集。这些数据集里面包含了很多行数据，每一行数据就是一个样本/数据点/数据样本。而每一行数据里面的预测信息称为特征/协变量，需要预测的目标为 标签/目标。
比如，想预测一下房价。先去网上拿到一堆交易信息（包含面积、房龄、销售价钱），这个就是训练数据集；每一个交易信息就是一个样本；每一个样本里面的面积、房龄 被叫做特征；想要预测的房屋价钱就是标签/目的。

1、线性模型

在线性模型中，我们假设因变量可以通过自变量的线性组合来描述。
当有 d 个特征时，线性模型的一般形式可以表示为：
$$\hat{y}=w_1{x}_1+w_2{x}_2+\cdots +w_d{x}_d+b$$
其中，$\hat{y}$ 是因变量（目标变量），$x_1,x_2,\dots ,x_d$ 是自变量（特征变量） $w_1,w_2,\dots ,w_d$ 是权重（斜率），$b$ 是偏置（截距）。
将所有特征放到向量 $\mathbf{x}\in {\mathbb{R}}^d$ 中， 并将所有权重放到向量 $\mathbf{w}\in {\mathbb{R}}^d$ 中， 我们可以用点积形式来简洁地表达模型：
$$\hat{y}=\mathbf{w}\top \mathbf{x}+b$$
此时 向量表示的是一个数据样本的特征，但是数据集有n个数据样本呢。所以咱们可以通过矩阵 $X\in {\mathbb{R}}^{n\times d}$ 将所有的特征引入。
$$\hat{y}=\mathbf{w}\top \mathbf{X}+b$$
这样线性模型就做好了。但是在寻找 模型参数 $\mathbf{w}$ $b$ 和 之前，还需要两个东西： （1）一种模型质量的度量方式； （2）一种能够更新模型以提高模型预测质量的方法。
这两个是干什么的呢？下面只是说一下作用，后面会详细讲的。
第一个：模型质量的的度量方式可以帮助我们评估模型预测的准确程度。
第二个：为了提高模型预测质量，我们需要使用一种有效的方法来更新模型的参数.

2、损失函数

损失函数用来衡量模型预测结果与真实值之间的差异或误差程度。通常我们会选择非负数作为损失，且数值越小表示损失越小，完美预测时的损失为0。
损失函数
下面咱们用比较常用的平分损失函数（MSE）
$$l^{(i)}(\mathbf{w},b)=\frac{1}{n}{\left({\hat{y}}^{(i)}-y^{(i)}\right)}^2$$
其中 $l^{(i)}(\mathbf{w},b)$ 表示在 样本i 时的误差，$n$ 为样本数量，${\hat{y}}^{(i)}$ 表示在 样本 i 的预测结果， $y^{(i)}$ 为真实的标签.
为了方便计算梯度时的推导和计算，咱们采用1/2来表示。
当我们对损失函数进行优化时，使用梯度下降等方法来更新模型参数。在求解梯度时，平方的因子 2 可以抵消平方项的导数中产生的系数，简化了计算过程。
$$l^{(i)}(\mathbf{w},b)=\frac{1}{2}{\left({\hat{y}}^{(i)}-y^{(i)}\right)}^2$$
为了度量模型在整个数量集上的质量，我们需要计算训练集 n 个样本上的损失均值（等价于求和）。
$$L(\mathbf{w},b)=\frac{1}{n}\sum _{i=1}^n{l}^{(i)}(\mathbf{w},b)=\frac{1}{n}\sum _{i=1}^n\frac{1}{2}{\left({\mathbf{w}}^{\top }{\mathbf{x}}^{(i)}+b-y^{(i)}\right)}^2$$
在训练模型时，我们希望寻找一组参数 （${\mathbf{w}}^{\ast },b^{\ast }$）， 这组参数能最小化在所有训练样本上的总损失。如下式：
$${\mathbf{w}}^{\ast },b^{\ast }=\underset{\mathbf{w},b}{\mathrm{argmin}}\mathbf{L}(\mathbf{w},b)$$

3、解析解

解析解是指直接通过数学公式来计算得到的闭式解，而不是通过迭代算法进行优化。对于线性回归问题，解析解可以使用最小二乘法来求解。
首先，需要把偏置 $b$ 合并到 $\mathbf{w}$ 中，让它也成为一个特征，权重为1。合并方法是在包含所有参数的矩阵中附加一列。
这样就变成了
$$L(\mathbf{w},b)=\frac{1}{n}\sum _{i=1}^n\frac{1}{2}{\left({\mathbf{w}}^{\top }{\mathbf{x}}^{(i)}-y^{(i)}\right)}^2=\frac{1}{2n}{\Vert \mathbf{wX}-\mathbf{y}\Vert }^2$$
最后对 $\mathbf{w}$ 求偏导
$$\begin{array}{rl}\frac{\partial }{\partial \mathbf{w}}L(\mathbf{w},b)& =\frac{\partial }{\partial \mathbf{w}}\frac{1}{2n}{\Vert \mathbf{wX}-\mathbf{y}\Vert }^2\\ & \\ & =\frac{1}{2n}\frac{\partial {\Vert \mathbf{wX}-\mathbf{y}\Vert }^2}{\partial \Vert \mathbf{wX}-\mathbf{y}\Vert }\frac{\partial }{\partial \mathbf{w}}\mathbf{wX}-y\\ & \\ & =\frac{1}{n}(\mathbf{wX}-\mathbf{y}{)}^T\mathbf{X}\end{array}$$
为什么结果是这个，是因为 $\frac{\partial }{\partial \mathbf{x}}{\Vert \mathbf{x}\Vert }^2=2{\mathbf{x}}^T,\frac{\partial \mathbf{Ax}}{\partial \mathbf{x}}=\mathbf{A}$ ，这样就理解了吧。
由于是最优解，所有求偏导后的一阶导数等于0。就可以得到 ${\mathbf{w}}^{\ast }$ 了。
$$\begin{array}{r}\frac{1}{n}(\mathbf{wX}-\mathbf{y}{)}^T\mathbf{X}=0\\ \\ \frac{1}{n}({\mathbf{X}}^T{\mathbf{w}}^T-{\mathbf{y}}^T)\mathbf{X}=0\\ \\ \frac{1}{n}({\mathbf{X}}^T{\mathbf{w}}^T\mathbf{X}-{\mathbf{y}}^T\mathbf{X})=0\\ \\ \frac{1}{n}({\mathbf{X}}^T\mathbf{Xw}-{\mathbf{X}}^T\mathbf{y})=0\\ \\ {\mathbf{X}}^T\mathbf{Xw}-{\mathbf{X}}^T\mathbf{y}=0\\ \\ {\mathbf{X}}^T\mathbf{Xw}={\mathbf{X}}^T\mathbf{y}\\ \\ ({\mathbf{X}}^T\mathbf{X}{)}^{-1}({\mathbf{X}}^T\mathbf{Xw})=({\mathbf{X}}^T\mathbf{X}{)}^{-1}{\mathbf{X}}^T\mathbf{y}\\ \\ \mathbf{w}=({\mathbf{X}}^T\mathbf{X}{)}^{-1}{\mathbf{X}}^T\mathbf{y}\end{array}$$

4、随机梯度下降

但模型没有显示解（解析解）的时候，我们只能选择迭代优化算法来确定模型的参数。
我们选择的是 梯度下降（gradient descent）的方法， 这种方法几乎可以优化所有深度学习模型。 它通过不断地在损失函数递减的方向上更新参数来降低误差。
首先，我们需要确定初始化模型参数 $\mathbf{w}$ ，接下来重复迭代更新参数t=1、2、3、…、n，更新公式为
$${\mathbf{w}}_t={\mathbf{w}}_{t-1}-\eta \frac{\partial \ell }{\partial {\mathbf{w}}_{t-1}}$$

在公式中 ${\mathbf{w}}_{t-1}$ 为上一次的结果，$\eta$ 为学习率（步长的超参数），$\frac{\partial \ell }{\partial {\mathbf{w}}_{t-1}}$ 为当前参数的梯度。

$\eta$ 学习率不能过大、也不能过小。太小的话、步长走的有限，需要计算多次梯度，但是计算梯度比较贵；也不能过大，会越过最优点，导致一直在上下震荡。

$\frac{\partial \ell }{\partial {\mathbf{w}}_{t-1}}$ 是当前梯度，但是由于每一次更新参数前都需要变量整个数据集，比较慢。所以我们通常会在需要计算更新的时候 随机选取一小批样本，这种变体叫做小批量随机梯度下降。
在每次迭代中，我们首先随机抽样一个小批量 $\mathcal{B}$ ， 它是由固定数量的训练样本组成的。 然后，我们计算小批量的平均损失关于模型参数的导数（也可以称为梯度）。 最后，我们将梯度乘以一个预先确定的正数 $\eta$ ，并从当前参数的值中减掉。
$$(\mathbf{w},b)\leftarrow (\mathbf{w},b)-\frac{\eta }{|\mathcal{B}|}\sum _{i\in \mathcal{B}}{\partial }_{(\mathbf{w},b)}{l}^{(i)}(\mathbf{w},b)$$
对于平方损失和仿射变换，我们可以明确地写成如下形式:
$$\begin{array}{rl}\mathbf{w}& \leftarrow \mathbf{w}-\frac{\eta }{|\mathcal{B}|}\sum _{i\in \mathcal{B}}{\partial }_{\mathbf{w}}{l}^{(i)}(\mathbf{w},b)=\mathbf{w}-\frac{\eta }{|\mathcal{B}|}\sum _{i\in \mathcal{B}}{\mathbf{x}}^{(i)}\left({\mathbf{w}}^{\top }{\mathbf{x}}^{(i)}+b-y^{(i)}\right),\\ & \\ b& \leftarrow b-\frac{\eta }{|\mathcal{B}|}\sum _{i\in \mathcal{B}}{\partial }_b{l}^{(i)}(\mathbf{w},b)=b-\frac{\eta }{|\mathcal{B}|}\sum _{i\in \mathcal{B}}\left({\mathbf{w}}^{\top }{\mathbf{x}}^{(i)}+b-y^{(i)}\right)\end{array}$$
$\mathbf{w}$ 和 $\mathbf{x}$ 都是向量。
表示每个小批量中的样本数，这也称为批量大小（batch size）。$\mathcal{B}$ 表示学习率（learning rate）。 批量大小和学习率的值通常是手动预先指定，而不是通过模型训练得到的。 这些可以调整但不在训练过程中更新的参数称为超参数（hyperparameter）。 调参（hyperparameter tuning）是选择超参数的过程。 超参数通常是我们根据训练迭代结果来调整的， 而训练迭代结果是在独立的验证数据集（validation dataset）上评估得到的。
大家可以看一下下面的视频（十分形象），看看梯度下降法到底在做什么。
梯度下降法

根据这个视频，我理解的是如何走，才使得损失函数走到最小值。
已知函数沿着梯度（一个向量）走将取得最大值、变化最快、变化率最大；那么反方向走就是下降最快，走到最小值点即可！！

三、从零实现

首先，需要导入后面用到的包

!pip install d2l==0.14 #安装d2l包

%matplotlib inline # 显示Matplotlib绘图
import random
import torch
from d2l import torch as d2l #导入d2l里面的torch，别名为d2l

1、生成数据集

我们将根据带有噪声的线性模型构造一个人造数据集。 我们的任务是使用这个有限样本的数据集来恢复这个模型的参数。
下面代码生成一个包含1000个样本的数据集， 每个样本包含从标准正态分布中采样的2个特征。合成的数据集是一个矩阵 $\mathbf{X}\in {\mathbb{R}}^{1000\times 2}$
我们使用线性模型参数 $\mathbf{w}=[2,-3.4]\top 、b=4.2$ 和噪声项 $\epsilon $ 生成数据集、标签。
$$\mathbf{y}=\mathbf{Xw}+b+ϵ$$

def synthetic_data(w, b, num_examples):  #@save
    """生成y=Xw+b+噪声"""
    X = torch.normal(0, 1, (num_examples, len(w))) #生成均值为0、标准差为1、num_example * len(w)的矩阵
    y = torch.matmul(X, w) + b
    y += torch.normal(0, 0.01, y.shape) 添加噪声
    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)

torch.normal 用于生成服从正态分布（高斯分布）的随机张量。函数语法如下

torch.normal(mean, std, size, *, out=None, dtype=None, layout=torch.strided, device=None, requires_grad=False) -> Tensor
# mean：表示正态分布的均值。
# std：表示正态分布的标准差。
# size：表示要生成的随机张量的形状。

torch.matmul 用于执行矩阵相乘或张量相乘的操作。它支持多种维度的输入。

torch.matmul(input, other, *, out=None) -> Tensor
# input：要进行相乘操作的第一个输入张量或矩阵。
# other：要进行相乘操作的第二个输入张量或矩阵。
# 返回值是两个输入张量或矩阵相乘的结果。

注意：features是一个矩阵，每一行有两个特征，共有1000行；y是一个列向量，有1000行，每一行是一个标量（标签值）。
接下来咱们可以看一下散点图

d2l.set_figsize()
d2l.plt.scatter(features[:, 1].detach().numpy(), labels.detach().numpy(), 1);

2、读取数据集

还记得咱们更新参数的时候，需要选取一小部分样本来代替这个数据集进行更新模型。所以我们最好封装一个函数，来进行随机获取小批量样本。

def data_iter(batch_size, features, labels): # batch_size需要样本的长度
    num_examples = len(features) 
    indices = list(range(num_examples)) #创建一个下标索引数组
    # 这些样本是随机读取的，没有特定的顺序
    random.shuffle(indices) #把下标索引数组里面的顺序打乱
    for i in range(0, num_examples, batch_size): # 遍历下标索引数组，跨度为 batch_size
        batch_indices = torch.tensor(
            indices[i: min(i + batch_size, num_examples)]) #创建一个从 i ~ i + batch_size，min是为了防止越界
        yield features[batch_indices], labels[batch_indices]

yield 是一个Python关键字，用于定义生成器函数。当在函数中使用yield语句时，函数将成为一个生成器函数。
它的作用是使函数能够暂停执行并返回一个值，然后在下一次迭代时从离开的地方继续执行。每次调用生成器函数时，它会返回一个生成器对象，该对象可以用于迭代获取生成器函数产生的值。
可以读取一组试试

batch_size = 10

for X, y in data_iter(batch_size, features, labels):
    print(X, '\n', y)
    break

3、初始化模型参数

接下来我们利用 normal函数随机初始化 w，并把偏差设置为0

w = torch.normal(0, 0.01, size=(2,1), requires_grad=True) 
b = torch.zeros(1, requires_grad=True) 

接下来就是通过随机梯度下降来不断更新这些参数了。

4、定义模型

接下来就是定义一下模型，根据咱们上面的公式来用py写

def linreg(X, w, b): 
    """线性回归模型"""
    return torch.matmul(X, w) + b

5、定义损失函数

这里和咱们之前讲的一样，用平分损失函数。
注意：需要把两个张量形状转换为一样的。

def squared_loss(y_hat, y):  
    """均方损失"""
    return (y_hat - y.reshape(y_hat.shape)) ** 2 / 2 #没有除以均值，在计算梯度的时候再除

6、定义优化算法

这里咱们就不使用解析解了，直接用迭代优化算法来确定模型参数的值。这里采用随机梯度下降法。

def sgd(params, lr, batch_size):  #params是所有参数w、b的列表；lr为学习率、baatch_size为样本数量
    """小批量随机梯度下降"""
    with torch.no_grad(): #上下文，关闭梯度计算，确保在更新参数时不会计算梯度。
        for param in params:
            param -= lr * param.grad / batch_size #更新
            param.grad.zero_() #清空梯度，因为梯度默认累积

7、训练

接下来设置超参数后，就可以进行训练了。
在每次迭代中，我们读取一小批量训练样本，并通过我们的模型来获得一组预测。 计算完损失后，我们开始反向传播，存储每个参数的梯度。 最后，我们调用优化算法sgd来更新模型参数。
在每个迭代周期（epoch）中，我们使用data_iter函数遍历整个数据集， 并将训练数据集中所有样本都使用一次（假设样本数能够被批量大小整除）。 这里的迭代周期个数num_epochs和学习率lr都是超参数，分别设为3和0.03。

lr = 0.03 # 学习率
num_epochs = 3 #迭代周期
# 写不写下面两行都可以，写是因为下面这两个已经成为一个约定了，net默认是模型，loss默认为损失率
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中的所有元素被加到一起，backward的对象只能是标量
        # 并以此计算关于[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、生成数据集

还是和之前一样，手动生成。

def synthetic_data(w, b, num_examples): 
    """生成y=Xw+b+噪声"""
    X = torch.normal(0, 1, (num_examples, len(w))) #生成均值为0、标准差为1、num_example * len(w)的矩阵
    y = torch.matmul(X, w) + b
    y += torch.normal(0, 0.01, y.shape) 添加噪声
    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)

2、读取数据集

可以调用API来随机读取小批量数据，不需要手动创建下标数组，打乱随机取。

def load_array(data_arrays, batch_size, is_train=True):
    """构造一个PyTorch数据迭代器"""
    dataset = data.TensorDataset(*data_arrays) #创建一下数据集对象
    return data.DataLoader(dataset, batch_size, shuffle=is_train) #每一次随机返回 batch_size个样本

batch_size = 10
data_iter = load_array((features, labels), batch_size)

关于 data.TensorDataset、data.DataLoader这两个，下面来浅显解释一下。
data.TensorDataset就是把特征和标签打包在一起，使之形成一一对应的关系。
*data_arrays是解包，把元祖里面的内容解出来。
data.DataLoader 就是随机返回 batch_size个样本。
关于更详细的，请看下面这篇文章。
TensorDataSet 和 DataLoader
接下来可以看看是否可以正确读取

next(iter(data_iter)) #用于从数据迭代器 data_iter 中获取下一个批量的数据。

3、定义模型

在前面，我们是通过数学公式来编写代码定义模型的；但是当公式复杂的时候，这时候就需要简化这个过程了。
对于标准深度学习模型，我们可以使用框架的预定义好的层。这使我们只需关注使用哪些层来构造模型，而不必关注层的实现细节。

# nn是神经网络的缩写
from torch import nn

net = nn.Sequential(nn.Linear(2, 1)) # 输入是2维、输出为1维的模型

nn.Sequential 是 torch.nn 模块中的一个类，它可以按照顺序将多个层组合起来形成一个神经网络模型。
在这里，nn.Linear(2, 1) 创建了一个线性层（全连接层），输入特征的维度为 2，输出特征的维度为 1。这个线性层用于将输入特征的维度从 2 维降低到 1 维。
然后，通过 nn.Sequential 将这个线性层添加到神经网络模型中。nn.Sequential 接受一个或多个层作为参数，并按照传入的顺序将这些层组合起来形成一个序列化的神经网络模型。
最后，将创建的神经网络模型赋值给变量 net。

4、初始化模型参数

深度学习框架通常有预定义的方法来初始化参数。 在这里，我们指定每个权重参数应该从均值为0、标准差为0.01的正态分布中随机采样， 偏置参数将初始化为零。

net[0].weight.data.normal_(0, 0.01)
net[0].bias.data.fill_(0)

5、定义损失函数

计算均方误差使用的是MSELoss类，也称为平方 $L_2$ 范数。 默认情况下，它返回所有样本损失的平均值。

loss = nn.MSELoss()

6、定义优化算法

小批量随机梯度下降算法是一种优化神经网络的标准工具， PyTorch在optim模块中实现了该算法的许多变种。
只需要传入优化的参数、学习率即可

trainer = torch.optim.SGD(net.parameters(), lr=0.03)

7、训练

在每个迭代周期里，我们将完整遍历一次数据集（train_data）， 不停地从中获取一个小批量的输入和相应的标签。 对于每一个小批量，我们会进行以下步骤:

• 通过调用net(X)生成预测并计算损失l（前向传播）
• 通过进行反向传播来计算梯度
• 通过调用优化器来更新模型参数

num_epochs = 3
for epoch in range(num_epochs):
    for X, y in data_iter:
        l = loss(net(X) ,y) #计算损失
        trainer.zero_grad() #清空梯度
        l.backward() #pytorch已经帮我们求和降维了，求导
        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)

参考文章
动手学深度学习