Deep Learning(二) 再论线性回归

再论线性回归

---

模型描述

在机器学习领域中的大多数任务通常都与预测有关。当我们想预测一个数值时，就会涉及到回归问题。回归(regression)是能为一个或多个自变量与因变量之间关系建模的一类方法，表示输入和输出的关系。

在机器学习日志（二）中，我们简单介绍了基于线性函数的线性回归模型，其模型假设如下

• $m$：训练集数据量
• $n$：特征数量
• $\boldsymbol{x}^{(i)}=(x_1^{(i)},x_2^{(i)},\cdots,x_n^{(i)})$：第 $i$ 个样本的特征向量
• $\boldsymbol{X}\in\mathbb{R}^{m\times n}$：样本特征矩阵，第 $i$ 行代表 $\boldsymbol{x}^{(i)}$
• $\boldsymbol{w}=(w_1,x_2,\cdots,x_n)^T$：权重向量(Weight)
• $b$：偏置标量
• $y^{(i)}$：第 $i$ 个样本的标签
• $\boldsymbol{y}=(y^{(1)},y^{(2)},\cdots,y^{(m)})^T$：样本标签矩阵
• $\hat{y}^{(i)}=\langle\boldsymbol{w},\boldsymbol{x}^{(i)}\rangle+b$：第 $i$ 个样本的估计值

线性模型可以看做一个单层神经网络，即只有一层输入层和一层输出层的神经网络模型。为了衡量该线性回归模型的预测水平，也就是考虑模型训练后，样本预测值和真实值的区别，定义平方损失函数

$$L(\boldsymbol{X},\boldsymbol{y},\boldsymbol{w},b)=\frac{1}{2m}\sum_{i=1}^m\big(y^{(i)}-\hat{y}^{(i)}\big)^2=\frac{1}{2m}\|\boldsymbol{y}-\boldsymbol{Xw}-b\|^2$$

这个损失函数应当尽可能的小，因此应当选取合适的参数

$$\boldsymbol{w}^{*},b^{*}=\arg\min_{\boldsymbol{w},b} L(\boldsymbol{X},\boldsymbol{y},\boldsymbol{w},b)$$

上式中，$\arg\min$ 表示取得最小值对应的自变量值。

正规方程解

本篇用数学方法推导正规方程法的显示解。现在只考虑关于变元 $\boldsymbol{w}$ 的损失函数，把原符号假设进行修改

$$\boldsymbol{X}\leftarrow [\boldsymbol{X},\boldsymbol{1}]\in\mathbb{R}^{m\times(n+1)},\quad\boldsymbol{w}\leftarrow\begin{bmatrix}\boldsymbol{w}\\b\end{bmatrix}\in\mathbb{R}^{n+1}$$

此时，平方损失函数就变成了只关于 $\boldsymbol{w}$ 的函数

$$L(\boldsymbol{X},\boldsymbol{y},\boldsymbol{w})=\frac{1}{2m}\|\boldsymbol{y}-\boldsymbol{Xw}\|^2$$

这个函数是一个凸函数，因此一定有一个全局最小值，满足梯度为0，即有

$$\frac{\partial L}{\partial \boldsymbol{w}}=-\frac{1}{m}\big(\boldsymbol{y}-\boldsymbol{Xw}\big)^T\boldsymbol{X}=\boldsymbol{0}^T$$

消去负号和系数，两边取转置得到

$$\boldsymbol{X}^T\big(\boldsymbol{y}-\boldsymbol{Xw}\big)=\boldsymbol{0}$$

展开括号并移项得到

$$\boldsymbol{X}^T\boldsymbol{Xw}=\boldsymbol{X}^T\boldsymbol{y}$$

当矩阵 $\boldsymbol{X}^T\boldsymbol{X}$ 可逆时，两边乘以 $(\boldsymbol{X}^T\boldsymbol{X})^{-1}$ 便得到解

$$\boldsymbol{w}^* = (\boldsymbol{X}^T \boldsymbol{X})^{-1} \boldsymbol{X}^T \boldsymbol{y}$$

梯度下降优化

机器学习篇中已经介绍，当矩阵 $\boldsymbol{X}^T\boldsymbol{X}$ 不可逆时，无法直接通过正规方程法得到解。此时采用梯度下降优化是一个较优的选择，我们重新回顾一下梯度下降的步骤，并重新规定部分符号

1. 挑选一个初始权重 $\boldsymbol{w}_0$
2. 重复执行梯度下降迭代

$$\boldsymbol{w}_t=\boldsymbol{w}_{t-1}-\eta\frac{\partial L}{\partial \boldsymbol{w}_{t-1}}$$

3. 权重收敛时停止算法

在上述过程中，$\eta$ 是学习率，决定了每次迭代的步长。我们在机器学习篇中提到：学习率不应该过小或者过大，但有时候较小值往往是一个优先考虑的选择。然而，在整个数据集上应用梯度下降往往会降低运行效率，因此我们可以考虑选取一部分样本作为梯度下降的参数数据，称作小批量随机梯度下降(Mini-Batch Stochastic Gradient Descent)。

在每次迭代中，我们首先随机抽样一个小批量 $\mathcal{B}$，其具有一个固定的批量大小 $|\mathcal{B}|$，即由 $|\mathcal{B}|$ 个训练样本组成，然后执行梯度下降的迭代步骤，即有

$$(\boldsymbol{w},b)\leftarrow(\boldsymbol{w},b)-\frac{\eta}{|\mathcal{B}|}\sum_{i\in\mathcal{B}}\partial_{(\boldsymbol{w,b})}L(\boldsymbol{x}^{(i)},y^{(i)},\boldsymbol{w},b)$$

这里，批量大小 $|\mathcal{B}|$ 是一个超参数，不应该过大或者过小。

线性回归代码

---

这一部分，我们参考沐神的课程介绍线性回归的代码实现。首先，我们需要导入必要的包。

# 导入包
%matplotlib inline
import random
import torch
from d2l import torch as d2l

我们采用人工生成数据集的形式，定义 synthetic_data 随机生成若干份样本数据。

# 定义数据集构造函数
def synthetic_data(w, b, num_examples):
    X = torch.normal(0, 1, (num_examples, len(w))) # 生成服从 N(0,1) 的样本矩阵
    y = torch.matmul(X, w) + b
    y += torch.normal(0, 0.01, y.shape) # 引入噪音 epsilon
    return X, y.reshape(-1, 1) # 自动补全为列向量

# 生成 1000 样本的数据集
true_w = torch.tensor([2, -3.4])
true_b = 4.2
features, labels = synthetic_data(true_w, true_b, 1000)

为了使用小批量随机梯度下降方法，定义 data_iter 随机抽取小批量样本。

# 定义随机读取批量函数
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):
        batch_indices = torch.tensor(indices[i: min(i + batch_size, num_examples)]) # 取最小值避免溢出
        yield features[batch_indices], labels[batch_indices] # 迭代器返回当前值

# 定义小批量随机梯度下降函数
def sgd(params, lr, batch_size):
    with torch.no_grad():
        for param in params:
            param -= lr * param.grad / batch_size
            param.grad.zero_() # 将参数的梯度清零

接着，随机初始化参数，并定义线性回归的基本模型和损失函数。

# 设置初始参数
w = torch.normal(0, 0.01, size=(2, 1), requires_grad=True) # 生成服从 N(0,0.01) 的随机权重
b = torch.zeros(1, requires_grad=True) # 初始化偏差为 0

# 定义线性计算模型
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

最后，利用 for 循环执行训练过程，并输出参数偏差以评估训练结果。

# 训练过程
lr = 0.03 # 学习率
num_epochs = 10 # 迭代次数
batch_size = 10 # 批量大小
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) # 小批量损失
        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}') # 输出当前 epoch 和损失

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

一般情况下，会得到如下形式的训练结果。

epoch 1, loss 0.036340
epoch 2, loss 0.000131
epoch 3, loss 0.000049
epoch 4, loss 0.000049
epoch 5, loss 0.000049
w的估计误差: tensor([-0.0004,  0.0003], grad_fn=<SubBackward0>)
b的估计误差: tensor([0.0002], grad_fn=<RsubBackward1>)

事实上，很多我们刚定义的函数都可以调用已有库，这里只是为了展现所有环节的原理。

Softmax 回归

---

回归与分类

在机器学习篇中，我们简单介绍了基于线性回归模型改进的Logistic回归模型。虽然其名称中包含回归二字，但实际上是一个分类模型，尤其是对于二分类问题。一般而言，我们可以把回归和分类问题看成两种不同的神经网络模型。

区别	回归模型	分类模型
输出个数	单个数值	多个数值
输出区间	自然区间 $\mathbb{R}$	置信度 $[0,1]$
损失	与真实值的区别	与真实类别是否一致

[{"url":"/img/Pytorch/p2f1.png","alt":"回归与分类","title":""}]

若一个分类问题中包含 $K$ 个类别，可以采用独热编码(One-Hot Encoding)的形式。独热编码是一个二进制向量，它的分量和类别一样多，对于任意一个表示分类结果的标签 $y=1,2,\cdots,K$，可以将其看做

$$y\rightarrow \boldsymbol{y}=(y_1,y_2,\cdots,y_K)^T,\quad y_i=\begin{cases}1,&i=y\\0,&else\end{cases}$$

分类模型也可以视作一个单层神经网络模型，输出层分别为 $o_1,o_2,\cdots,o_K$，如果使用线性模型来表示这个神经网络的参数传递，那么就有

$$o_i=\sum_{j=1}^nx_jw_{ij}+b_i$$

为了更简洁的表示这个模型，写成向量的形式有

$$\boldsymbol{o}=\boldsymbol{Wx}+\boldsymbol{b}$$

Softmax 函数

由于输出层的值 $o_i$ 表示了 $y=i$ 的置信度，因此预测结果可以表示为

$$\hat{y}=\arg\max_i o_i$$

不仅如此，我们还希望模型对正确类别的置信度大于错误类别的置信度。举例而言，如果我们要实现一个从动物种类中识别猫咪的分类模型，我们希望：对于一个猫咪，模型预测为猫的置信度远大于预测为狗的置信度、或者远大于预测为鸡的置信度。用数学语言表示为

$$o_y-o_i\geqslant\Delta(y,i),\quad\Delta\text{ is a threshold}$$

然而，线性模型容易带来一个缺陷，即所有输出之和不一定为1，且部分输出可能为负。如果真的是这样，这显然违背了概率论的原则要求，因为预测结果的事件构成了一个互斥完备事件组。为解决这一问题，人们定义了Softmax函数，即有

$$\hat{\boldsymbol{y}}=\mathrm{softmax}(\boldsymbol{o}),\quad \hat{y}_i=\frac{\exp(o_i)}{\displaystyle\sum_{j=1}^K \exp(o_j)}$$

注意到 softmax 函数有如下性质

• 对于每一个处理后的输出 $\hat{y}_i$ 都有 $0\leqslant \hat{y}_i\leqslant 1$
• 处理后的输出 $\hat{\boldsymbol{y}}$ 未改变 $\boldsymbol{o}$ 中元素的大小次序，即

$$\arg\max_i \hat{y}_i=\arg\max_io_i$$

交叉熵损失

在 softmax 回归或者是之前提到的 Logistic 回归中，运用平方损失函数（即L2损失）往往会导致优化结果收敛于局部最优解，因此需要改变损失函数的形式。softmax 函数给出了一个处理后的向量 $\hat{\boldsymbol{y}}$，根据置信度的概念有

$$\hat{y}_i=P(y=i|\boldsymbol{x})$$

假设整个数据集 $\{\boldsymbol{X},\boldsymbol{Y}\}$ 有 $m$ 个样本，其中第 $i$ 个样本由特征向量 $\boldsymbol{x}^{(i)}\in\mathbb{R}^n$ 和独热标签向量 $\boldsymbol{y}^{(i)}\in\mathbb{R}^K$ 组成。根据数理统计的知识，可以构造似然函数来比较预测值和真实值的区别

$$P(\boldsymbol{Y}|\boldsymbol{X})=\prod_{i=1}^mP(\boldsymbol{y}^{(i)}|\boldsymbol{x}^{(i)})$$

想要最大化似然估计，即最小化负对数似然估计

$$-\log P(\boldsymbol{Y}|\boldsymbol{X})=-\sum_{i=1}^m\log P(\boldsymbol{y}^{(i)}|\boldsymbol{x}^{(i)})$$

由于 $P(\boldsymbol{y}^{(i)}|\boldsymbol{x}^{(i)})=\hat{y}^{(i)}_y$，其中下标 $y\in\{1,2,\cdots,K\}$ 是该样本的真实类别。定义交叉熵损失(Cross-Entropy Loss)函数

$$L(\boldsymbol{y},\hat{\boldsymbol{y}})=-\sum_{i=1}^K y_i\log \hat{y}_i=-\log \hat{y}_y$$

因此，代入原式有

$$-\log P(\boldsymbol{Y}|\boldsymbol{X})=\sum_{i=1}^mL(\boldsymbol{y}^{(i)},\hat{\boldsymbol{y}}^{(i)})$$

可以证明，使用交叉熵损失函数能够把参数训练过程变为一个凸优化问题。

有关凸优化的问题，我会后续单独开一篇日志。

梯度计算

为了使用梯度下降法，需要计算损失函数的梯度。根据上一小节有

$$\begin{aligned}L(\boldsymbol{y},\hat{\boldsymbol{y}})&=-\sum_{i=1}^Ky_i\log\frac{\exp(o_i)}{\displaystyle\sum_{j=1}^K \exp(o_j)}\\&=\sum_{i=1}^Ky_i\log\sum_{j=1}^K \exp(o_j)-\sum_{i=1}^K y_io_i\\&=\log\sum_{j=1}^K\exp(o_j)-\sum_{i=1}^Ky_io_i\\ \end{aligned}$$

因此，对任何未规范化的预测值 $o_i$，其导数为

$$\frac{\partial L}{\partial o_i}=\frac{\exp(o_i)}{\displaystyle\sum_{j=1}^K \exp(o_j)}-y_i=\mathrm{softmax}(\boldsymbol{o})_i-y_i$$

可以发现，这个导数结果是 softmax 模型分配的概率与实际发生的情况（由独热标签向量表示）之间的差异。从这个意义上讲，这与我们在回归中看到的非常相似，其中梯度是观测值和估计值之间的差异。

Softmax 回归代码

---

下载数据集

我们使用 Fashion-MNIST 数据集作为训练集，点此前往下载数据集，把四个 .gz 文件下载到本地。

FashionMNIST 是一个替代 MNIST 手写数字集的图像数据集。 它是由 Zalando 公司旗下的研究部门提供。其涵盖了来自 10 种类别的共 7 万个不同商品的正面图片，图片数据和 MNIST 数据集一样均为 $28\times 28$ 像素尺寸，是目前分类模型实战项目的优秀选择。

需要把数据集文件按照正确地址存放，data\ 及下属目录命名必须一致，其他无要求。

Root\
  ├── code\                   
  │    └── softmax.ipynb
  └── data\
       └── FashionMNIST\
            └── raw\
                 ├── t10k-images-idx3-ubyte.gz
                 ├── t10k-labels-idx1-ubyte.gz
                 ├── train-images-idx3-ubyte.gz
                 └── train-labels-idx1-ubyte.gz

代码

Softmax 回归代码和线性回归代码逻辑基本一致，这里不再分块讲解，可参考沐神原视频。

import torch
from IPython import display
from d2l import torch as d2l

# 检测当前设备
device = torch.device('cuda' if torch.cuda.is_available() else 'cpu')
print(f"当前使用的设备: {device}") 

# 调用 Fashion-MNIST 数据集
batch_size = 256
train_iter, test_iter =  d2l.load_data_fashion_mnist(batch_size)

# 原始数据为 28*28 的图像，展平为 784 个特征，输出类别为 10 个
num_inputs = 784
num_outputs = 10

W = torch.normal(0, 0.01, size=(num_inputs, num_outputs), requires_grad=True, device=device)
b = torch.zeros(num_outputs, requires_grad=True, device=device)

# 定义 softmax 函数
def softmax(X):
    X_exp = torch.exp(X)
    partition = X_exp.sum(1, keepdim=True)
    return X_exp / partition  # 这里应用了广播机制

# 定义 softmax 回归模型
def net(X):
    return softmax(torch.matmul(X.reshape((-1, W.shape[0])), W) + b) # reshape 展平为向量

# 定义交叉熵损失函数
def cross_entropy(y_hat, y):
    return - torch.log(y_hat[range(len(y_hat)), y])

# 定义分类准确率函数
def accuracy(y_hat, y):
    if len(y_hat.shape) > 1 and y_hat.shape[1] > 1:
        y_hat = y_hat.argmax(axis=1)
    cmp = y_hat.type(y.dtype) == y
    return float(cmp.type(y.dtype).sum())

# 定义一个累加器的类 Accumulator
class Accumulator:
    def __init__(self, n):
        self.data = [0.0] * n

    def add(self, *args):
        self.data = [a + float(b) for a, b in zip(self.data, args)]

    def reset(self):
        self.data = [0.0] * len(self.data)

    def __getitem__(self, idx):
        return self.data[idx]

# 定义给定模型的准确率函数
def evaluate_accuracy(net, data_iter):
    if isinstance(net, torch.nn.Module):
        net.eval()  # 将模型设置为评估模式
    metric = Accumulator(2)  # 正确预测数、预测总数
    with torch.no_grad():
        for X, y in data_iter:
            X = X.to(device)
            y = y.to(device)
            metric.add(accuracy(net(X), y), y.numel())
    return metric[0] / metric[1]

# 定义一个方便绘图的类 Animator
class Animator:
    def __init__(self, xlabel=None, ylabel=None, legend=None, xlim=None,
                 ylim=None, xscale='linear', yscale='linear',
                 fmts=('-', 'm--', 'g-.', 'r:'), nrows=1, ncols=1,
                 figsize=(3.5, 2.5)):
        # 增量地绘制多条线
        if legend is None:
            legend = []
        d2l.use_svg_display()
        self.fig, self.axes = d2l.plt.subplots(nrows, ncols, figsize=figsize)
        if nrows * ncols == 1:
            self.axes = [self.axes, ]
        # 使用lambda函数捕获参数
        self.config_axes = lambda: d2l.set_axes(
            self.axes[0], xlabel, ylabel, xlim, ylim, xscale, yscale, legend)
        self.X, self.Y, self.fmts = None, None, fmts

    def add(self, x, y):
        # 向图表中添加多个数据点
        if not hasattr(y, "__len__"):
            y = [y]
        n = len(y)
        if not hasattr(x, "__len__"):
            x = [x] * n
        if not self.X:
            self.X = [[] for _ in range(n)]
        if not self.Y:
            self.Y = [[] for _ in range(n)]
        for i, (a, b) in enumerate(zip(x, y)):
            if a is not None and b is not None:
                self.X[i].append(a)
                self.Y[i].append(b)
        self.axes[0].cla()
        for x, y, fmt in zip(self.X, self.Y, self.fmts):
            self.axes[0].plot(x, y, fmt)
        self.config_axes()
        display.display(self.fig)
        display.clear_output(wait=True)

# 定义单 epoch 训练函数
def train_epoch(net, train_iter, loss, updater):
    # 将模型设置为训练模式
    if isinstance(net, torch.nn.Module):
        net.train()
    # 训练损失总和、训练准确度总和、样本数
    metric = Accumulator(3)
    for X, y in train_iter:
        X = X.to(device)
        y = y.to(device)

        # 计算梯度并更新参数
        y_hat = net(X)
        l = loss(y_hat, y)
        if isinstance(updater, torch.optim.Optimizer):
            # 使用PyTorch内置的优化器和损失函数
            updater.zero_grad()
            l.mean().backward()
            updater.step()
        else:
            # 使用定制的优化器和损失函数
            l.sum().backward()
            updater(X.shape[0])
        metric.add(float(l.sum()), accuracy(y_hat, y), y.numel())
    # 返回训练损失和训练精度
    return metric[0] / metric[2], metric[1] / metric[2]

# 定义训练和可视化函数
def train(net, train_iter, test_iter, loss, num_epochs, updater):
    animator = Animator(xlabel='epoch', xlim=[1, num_epochs], ylim=[0.3, 0.9],
                        legend=['train loss', 'train acc', 'test acc'])
    for epoch in range(num_epochs):
        train_metrics = train_epoch(net, train_iter, loss, updater)
        test_acc = evaluate_accuracy(net, test_iter)
        animator.add(epoch + 1, train_metrics + (test_acc,))
    train_loss, train_acc = train_metrics
    assert train_loss < 0.5, train_loss
    assert train_acc <= 1 and train_acc > 0.7, train_acc
    assert test_acc <= 1 and test_acc > 0.7, test_acc

# 设置超参数
lr = 0.1
num_epochs = 10

def updater(batch_size):
    return d2l.sgd([W, b], lr, batch_size)

train(net, train_iter, test_iter, cross_entropy, num_epochs, updater)

运行代码，应当得到和下图类似的训练结果。

[{"url":"/img/Pytorch/p2f2.png","alt":"训练结果示意图","title":""}]

若出现运行缓慢等情况是正常的，耐心等待即可，一般是由于CPU瓶颈导致的。

参考网站

课程主页

沐神伟大

课程教材

无需多言

PyTorch安装

PyTorch安装教程