【深度学习】线性回归的从零开始实现与简洁实现

前言

我原本后面打算用李沐老师那本《动手学深度学习》继续“抄书”，他们团队也免费提供了电子版(https://zh-v2.d2l.ai/d2l-zh-pytorch.pdf)。但书里涉及到代码，一方面展示起来不太方便，另一方面我自己也有很多地方看不太懂。

这让我开始思考起我“抄书”的意义了。如果都是文字，我感觉抄起来可以加深印象，在抄的同时理解并思考。

但是如果涉及到代码，我没办法在抄的时候就能理解，就能很顺畅地把代码解释出来。需要自己先去调试一遍，才能心安理得地介绍起来。

于是，我开始停止自己看书，打算先去看看视频。土堆的视频先从加载数据讲起，然后是 TensorBoard，我感觉和前面学到的训练/测试什么的，根本没关系，也不知道为什么要讲这些，就听不下去了，想着先去看看实操。

然后我到看刘二大人的教程，他先讲线性模型和梯度下降等等的原理，再上代码实操，看起来就比较舒服。我把他的代码自己敲了一遍，虽说能跑起来，但还是存在很多疑惑。

控制台不停地打印着训练轮数和损失，看起来有点样子，但是不直观。我想着要是能把训练过程中的 loss 作图出来就好。虽说可以利用 python 的作图包绘制，但我感觉应该会有专门的工具来实现。

我上网搜了下，发现 pytorch 就自带了 tensorboard 这个工具可以可视化训练过程。哎，这不就是土堆视频里讲的那个么。我就又老老实实去看土堆的视频了。

经过这些过程，我早已经变得很挫败，甚至感觉自己不是学这个的料。东看一点，西看一点，最后还是回到了李沐老师的书上。

我感觉自已经走了很多弯路，《动手学深度学习》这本书里最前面一章就是介绍预备知识，我嫌繁琐，跳过没看，直接看后面的内容，但我又看不懂，又要倒回去重新看预备知识，这样就会产生很多挫败情绪。

虽说经过这几天折腾，隐隐约约学到了一些东西，但显然是划不来的。我以后就老老实实跟着李沐老师，先认认真真搞懂基础的几个模型，不去担心可能学了以后用不到之类的问题，之后再重点看看自己导师研究的方向。

下面我们从线性回归开始，先介绍它的关键思想，然后手动实现它，体会这一过程，再借助 pytorch 来实现。

通过对比这两种方式，可以更好地理解 pytorch 的用法。

一、线性回归关键思想

回归(regression)是能为一个或多个自变量与因变量之间关系建模的一类方法。在自然科学和社会科学领域，回归经常用来表示输入和输出之间的关系。

线性回归(linear regression)指输入 $\pmb{x}$ 和输出 $y$ 之间的关系是线性的，即 $y$ 可以表示为 $\pmb{x}$ 中元素的加权和，如下式： $y=w_1\cdot x_{1}+w_{2}\cdot x_{2}+b.$ 其中， $w_{1}$ 和 $w_{2}$ 称为权重，其决定了每个 $x$ 对 $y$ 的影响； $b$ 称为偏置(bias)。

偏置是指当所有 $x$ 都取 0 时，输出值应该为多少。如果没有偏置项，模型的表达能力将受到限制。

给定一个数据集，我们的目标是寻找模型的权重 $\pmb{w}$ 和偏置 $b$ ，使得根据模型做出的预测大体符合真实情况。

在开始寻找最好的模型参数(model parameters)前，我们还需要两个东西：

一种模型质量的度量方式；（怎么判断这个模型好不好？）
一种能够更新模型以提高模型质量的方法。（如果模型质量不太行，怎么提高？）

损失函数（模型质量度量方式）

损失函数(loss function)能够量化目标的实际值与预测值之间的差距。通常我们会选择非负数作为损失，且数值越小表示损失越小。

回归问题中最常用的损失函数是平方误差函数。

可以想想我们以前学直线拟合时的最小二乘法。

随机梯度下降（提高模型质量方法）

梯度下降(gradient descent)通过不断地在损失函数递减的方向上更新参数来降低误差。

其最简单的方法是计算损失函数（所有样本的损失均值）关于模型参数的导数。但实际中的执行可能会非常慢：

因为在每一次更新参数前，我们必须遍历整个数据集。因此，我们通常会在每次需要计算更新的时候随机抽取一小批样本，这种变体称为小批量随机梯度下降(minibatch stochastic gradient descent)。

在每次迭代中，我们首先随机抽样一个小批量 $B$ ，它是又固定数量的训练样本组成。然后，我们计算小批量的平均损失关于模型参数的导数。最后，我们将梯度乘以一个预先确定的正数 $\eta$ ，并从当前参数的值中减掉。

总结起来，算法的步骤为：

初始化模型参数；
从数据集中随机抽取小批量样本且在负梯度的方向上更新参数，并不断迭代这一步骤。

正态分布与平方损失

正态分布和线性回归之间的关系很密切。正态分布(normal distribution)，也称为高斯分布(Gaussian distribution)，最早由德国数学家高斯(Gauss)应用于天文学研究。

简单地说，若随机变量 $X$ 具有均值 $\mu$ 和方差 $\sigma^2$ ，其正态概率密度函数为： $p(x)=\frac{1}{\sqrt{2\pi\sigma^2}}\exp\bigg(-\frac{1}{2\sigma^2}(x-\mu)^2\bigg).$ 正态分布的可视化图像如下图所示。改变均值会产生沿 $x$ 轴的偏移，增加方差将会分散分布、降低其峰值。

不同参数的正态分布图像

均方误差损失函数可以用于线性回归的一个原因是：我们假设了观测中包含噪声，其中噪声服从正态分布。

在高斯噪声的假设下，最小化均方误差等价于对线性模型的极大似然估计。

二、线性回归的从零开始实现

导入的包有：

import random
import torch
import matplotlib.pyplot as plt

生成数据集

我们可以自己用带有噪声的线性模型生成一系列数据，通过优化算法不断减小损失，最终尽可能恢复用来生成数据的参数。

例如，用 $y=2x_1-3.4x_2+4.2+\epsilon$ 来生成一系列数据点，其中 $\epsilon$ 服从均值为 0，标准差为 0.01 的正态分布。

为提高计算性能，我们可以用矩阵进行表示，即 $\begin{bmatrix} y_1 \\\ y_2 \\\ \vdots \\\ y_n \end{bmatrix} = \begin{bmatrix} x_1 & x_2 \\\ x_1 & x_2 \\\ \cdots & \cdots \\\ x_1 & x_2 \end{bmatrix}_{n\times 2} \begin{bmatrix} 2 \\\ -3.4 \end{bmatrix}+4.2+\epsilon$ 后面两项常数项可以看成 $n\times1$ 的矩阵，由于张量具有广播机制，可以直接相加。

生成数据并可视化 $y$ 与 $x_2$ 关系的代码如下：

import torch
import matplotlib.pyplot as plt

num_samples = 1000  # 生成的数据点个数
# 原模型的两个参数
w_true = torch.tensor([2, -3.4])
b_true = 4.2
# 设置随机种子，方便验证和复现
random_seed = 326
torch.manual_seed(random_seed)
# 随机生成自变量
X = torch.randn((num_samples, len(w_true)))
# 生成不含噪声的 y
y = torch.matmul(X, w_true) + b_true
# 生成带有噪声的y
y += torch.normal(0, 0.01, y.shape)
# 绘制数据点，看看长什么样
# 自变量有两个，只绘制一个自变量与y的关系
plt.scatter(X[:, 1], y)
plt.show()

y与第二个自变量的散点图

为提高代码的复用性，可以把生成数据的部分封装成一个函数，这样如果要生成其他参数的数据时就更加方便。后续很多代码都会尽量用函数的形式表达。

# 生成数据集
def generate_data(num_samples, w, b):
    # 随机生成自变量
    x = torch.randn((num_samples, len(w)))
    # 生成不含噪声的 y
    y = torch.matmul(x, w) + b
    # 生成带有噪声的y
    y += torch.normal(0, 0.01, y.shape)
    print(y.shape)
    return x, y.reshape((-1, 1))

读取数据

如果我们采用小批量随机梯度下降作为优化算法，需要读取数据集并每次抽取一小批。

我们采用随机抽取的方式，即需要将数据集的样本打乱。编写的读取数据函数如下：

def load_data(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):
        batch_indices = torch.tensor(
            indices[i: min(i + batch_size, num_examples)])
        yield features[batch_indices], labels[batch_indices]

初始化参数

我们得先有一个最初的参数放进去，这样才可以去改进。如果这个初始参数选的好，可能很快就接近真实参数了。

初始权重 w 我们采用均值为 0，标准差为 0.01 的正态分布随机数，偏置 b 初始设为 0。

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

requires_grad=True表示需要对这些参数进行梯度的计算。

模型的定义

想想之前学过的拟合，可以用线性拟合，也可以用二次曲线拟合，都是为了将输入与输出关联起来。

可以把模型理解为一个函数，把输入放进去，可以得到一个输出。我们通过对原始数据的可视化分析，选取线性模型。

# 定义模型
def my_linear(x, w, b):
    # 线性模型
    return torch.matmul(x, w) + b

损失函数

采用平方损失函数来衡量模型好坏。注意需要对 y 使用 reshape，存在行向量与列向量的区别。

# 定义损失函数
def my_loss(y_hat, y):
    # 平方损失函数，y_hat为模型预测值
    # 这里reshape是因为读取出来的y是行向量，而预测值为列向量
    return (y_hat - y)**2 / 2

优化算法

采用小批量随机梯度下降算法进行优化。对参数进行更新，并清零梯度。

# 定义优化器
def my_optimizer(params, lr, batch_size):
    """小批量随机梯度下降"""
    with torch.no_grad(): # 以下操作不计算梯度
        for param in params:
            param -= lr * param.grad / batch_size
            param.grad.zero_()

训练模型

训练是一个反复的过程，从初始的模型参数出发，计算其预测值与损失，随后进行反向传播，得到损失函数对各个参数的梯度，最后利用优化算法对参数进行更新，转入下一轮训练。

学习率lr和训练轮数num_epochs是两个非常重要的超参数，需要反复试验不断调整，以获得较好的训练效果。

# 小批量的大小
batch_size = 10
# 学习率和训练轮数
lr = 0.03
num_epochs = 3

# 开始训练
for epoch in range(num_epochs):
    for x, y in load_data(batch_size, features, labels):
        l = my_loss(my_linear(x, w, b), y)  # X和y的小批量损失
        # 因为l形状是(batch_size,1)，而不是一个标量。
        # l中的所有元素被加到一起，并以此计算关于[w,b]的梯度
        l.sum().backward()
        my_optimizer([w, b], lr, batch_size)  # 使用参数的梯度更新参数
    with torch.no_grad():
        train_l = my_loss(my_linear(features, w, b), labels)
        # 打印出每轮训练的损失
        print(f'epoch {epoch + 1}, loss {float(train_l.mean()):f}')

# 评估训练效果
print(f'w的估计误差: {w_true - w.reshape(w_true.shape)}')
print(f'b的估计误差: {b_true - b}')

从打印出的结果可以发现，第一轮训练后的损失已经比较小了，且每轮训练后损失都在减小。最终得到的两个参数非常接近于我们用于生成数据的参数了。

epoch 1, loss 0.041334
epoch 2, loss 0.000153
epoch 3, loss 0.000048
w的估计误差: tensor([4.0889e-05, 1.1015e-04], grad_fn=<SubBackward0>)
b的估计误差: tensor([0.0004], grad_fn=<RsubBackward1>)

需要注意的是，我们不应该认为我们能完美地恢复参数，我们更关心如何高度准确预测参数。

三、线性回归的简洁实现

下面我们利用 pytorch 来实现线性回归模型，它可以帮助避免一些重复性的工作。需要导入的包有：

import torch
from torch.utils import data
from torch import nn

生成数据集

与从零实现同样，我们需要先有一个数据集。由于线性模型所需要的数据较为简单，可以自己生成。

我们采用前面写好的生成数据函数，利用随机种子生成一个相同的数据集。

读取数据集

pytorch 提供了 DataLoader 供我们进行数据读取，我们只需要传入相应格式的数据，传入参数batch_size和shuffle的值，就可以得到一个随机小批量数据迭代器。

def load_data(data_arrays, batch_size, is_train=True):     # 读取数据
    """构造一个PyTorch数据迭代器"""
    dataset = data.TensorDataset(*data_arrays)
    return data.DataLoader(dataset, batch_size, shuffle=is_train)

定义模型

在从零实现中，我们明确定义了线性模型的参数，并编写了内部进行计算的代码。如果模型的计算变得复杂，并且每天需要经常使用时，我们就会考虑对这一过程进行简化。

我们可以使用 pytorch 预先定义好的“层”，这样我们就只需要关注使用哪些层来构建我们的模型，而不用去关注怎么实现这个层。

我们首先定义一个模型变量 net，它是一个 Sequential 类的实例。Sequential 类可以将很多个层串联在一起。

当给定输入数据时，net 将数据传入第一层，然后将第一层的输出作为第二层的输入，以此类推。

线性回归可以看作一层的神经网络，实际上是不需要使用 Sequential 的。但由于以后使用的模型都会是很多层的，在这里使用 Sequtntial 可以帮助我们理解。

单层神经网络（线性回归）

上图中，它的每一个输入都通过矩阵-向量乘法得到它的输出，被称为全连接层(fully-connected layer)。

在 pytorch 中，全连接层在 nn.Linear 类中定义，它的第一参数是输入的形状，我们有两个 $x$ ，所以是 2。它的第二个参数是输出的形状，线性模型的 $y$ 是标量，所以是 1。

net = nn.Sequential(nn.Linear(2, 1))    # 定义模型

初始化模型参数

可以把 net 看成一个列表，用 net[0] 访问模型中的第一层，用 weight.data 和 bias.data 访问里面的参数。

我们还可以用 noraml_ 和 fill_ 来重写参数值，完成初始化。

同样，初始权重 w 我们采用均值为 0，标准差为 0.01 的正态分布随机数，偏置 b 初始设为 0。

# 初始化模型参数
net[0].weight.data.normal_(0, 0.01)
net[0].bias.data.fill_(0)

定义损失函数

同定义模型类似，我们可以直接调用现成的类，平方误差使用的是 nn.MSELoss 类。

loss = nn.MSELoss()   # 损失函数

定义优化器

同样，优化算法直接调用现成的 optim.SGD，利用 net.paremeters() 可以直接获取模型中的参数传入，再传入 lr 值即可定义一个小批量随机梯度下降优化器。

lr = 0.03    # 学习率
optimizer = torch.optim.SGD(net.parameters(), lr=lr)   # 优化器

训练模型

和前面的训练逻辑类似，对于每一个小批量，我们会经历以下步骤：

出初始的参数出发，将数据代入模型生成输出并计算损失（前向传播）；
通过反向传播计算梯度；
通过调用优化器来更新参数。

num_epochs = 3    # 训练轮数
for epoch in range(num_epochs):      # 开始训练
    for X, y in data_iter:
        l = loss(net(X) ,y)      # 前向传播
        optimizer.zero_grad()    # 避免上一轮影响，清零梯度
        l.backward()             # 反向传播
        optimizer.step()         # 更新参数
    l = loss(net(features), labels)
    print(f'epoch {epoch + 1}, loss {l:f}')
w = net[0].weight.data
print('w的估计误差：', w_true - w.reshape(w_true.shape))
b = net[0].bias.data
print('b的估计误差：', b_true - b)

同样我们可以在每轮训练后输出当前的损失，以及打印训练完成时的参数和真实参数的偏差。

epoch 1, loss 0.000273
epoch 2, loss 0.000096
epoch 3, loss 0.000096
w的估计误差： tensor([0.0004, 0.0002])
b的估计误差： tensor([-0.0002])

和从零开始实现一样，最后估计得到的参数非常接近真实参数。

两种方式实现线性神经网络的完整 py 文件见附件。