机器学习深度学习——softmax回归从零开始实现

👨‍🎓作者简介：一位即将上大四，正专攻机器学习的保研er
🌌上期文章：机器学习&&深度学习——向量求导问题
📚订阅专栏：机器学习&&深度学习
希望文章对你们有所帮助

就跟之前从零开始实现线性回归一样，softmax回归也很重要，因此也进行一次从0开始实现。之前的章节中，我们已经引入了Fashion-MNIST数据集，并设置数据迭代器的批量大小为256。

import torch
from IPython import display
from d2l import torch as d2lbatch_size = 256
train_iter, test_iter = d2l.load_data_fashion_mnist(batch_size)

初始化模型参数

和之前线性回归例子一样，每个样本都用固定长度的向量表示，则之前数据集中每个样本都是28×28的图像，将要进行展平，把他们看做是长度为784的向量。（在这里我们暂且把每个像素的位置都看作是一个特征，其实严格意义上要讨论其空间结构的，在这不做讨论）
而在softmax回归中，我们的输出和类别一样多，因为数据集由10个类别，所以网络输出维度为10。因此，权重将构成一个784×10的矩阵，偏置将构成一个1×10的行向量。与线性回归一样，我们将使用正态分布初始化我们的权重W，偏置初始化为0。

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

定义softmax操作

回顾sum运算符

按照之前的线性代数的内容，给定一个矩阵X，我们可以利用sum函数给所有元素求和（默认）。也可以对同一列（轴0）或同一行（轴1）进行求和。用例子表示：

X = torch.tensor([[1.0, 2.0, 3.0], [4.0, 5.0, 6.0]])
print(X.sum(0, keepdim=True), X.sum(1, keepdim=True))  # keepdim表示还保留着之前维度即二维

结果：

tensor([[5., 7., 9.]]) tensor([[ 6.],
[15.]])

构建softmax运算函数

回想一下softmax的三个步骤：
1、对每个项求幂（使用exp）；
2、对每一行求和（因为小批量中每一行就是一个样本），得到每个样本的规范化常数
3、将每一行除以其规范化常数，确保结果的和为1
回顾一下表达式：
$$softmax(X{)}_{ij}=\frac{exp(X_{ij})}{{\sum }_{k}exp(X_{ik})}$$

def softmax(X):X_exp = torch.exp(X)partition = X_exp.sum(1, keepdim=True)return X_exp / partition  # 广播机制

可以验证上述的代码：

X = torch.normal(0, 1, (2, 5))
X_prob = softmax(X)
print(X_prob, '\n', X_prob.sum(1))

结果：

tensor([[0.0152, 0.1212, 0.6149, 0.0877, 0.1610],
[0.1921, 0.0852, 0.1945, 0.4261, 0.1020]])
tensor([1.0000, 1.0000])

根据概率原理易得每行的和为1
注意：数学上看起来很正确，但是代码实现太草率了。矩阵中的非常大或非常小的元素可能造成数值上溢或下溢，但是这里没有采取措施来防止这一点。

定义模型

也就是直接将y=XW+b进行softmax运算得到，注意下面的X要使用reshape来将每张原始图像展平为向量（轴0放个-1让他自己根据列长度=784来进行运算，这里应为256，因为批量大小为256，每个批量（图像）都被展开成了784的向量）

def net(X):return softmax(torch.matmul(X.reshape((-1, W.shape[0])), W) + b)

定义损失函数

引入交叉熵函数，这在深度学习中很可能是最常见的损失函数了（目前分类问题数量远超回归问题数量）
回顾一下，交叉熵采用真实标签的预测概率的负对数似然。这边我们不使用for循环这种低效的方式，而是通过一个运算符选择所有函数。在这里我们先介绍下NumPy的整数数组索引。

NumPy的整数数组索引

整数数组索引，它可以选择数组中的任意一个元素，比如，选择第几行第几列的某个元素，示例如下：

import numpy as np
#创建二维数组
x = np.array([[1,  2],  [3,  4],  [5,  6]])
#[0,1,2]代表行索引;[0,1,0]代表列索引
y = x[[0,1,2],[0,1,0]] 
print (y)

结果：

[1 4 5]

对着样例做简单分析：将行、列索引组合会得到 (0,0)、(1,1) 和 (2,0) ，它们分别对应着输出结果在原数组中的索引位置。

下面，我们创建一个数据样本y_hat，其中包含2个样本在3个类别的预测概率，以及它们对应的标签y。然后使用y作为y_hat中概率的索引，我们选择第一个样本中第一个类的概率和第二个样本中第三个类的概率：

y = torch.tensor([0, 2])
y_hat = torch.tensor([[0.1, 0.3, 0.6], [0.3, 0.2, 0.5]])
print(y_hat[[0, 1], y])

输出：

tensor([0.1000, 0.5000])

交叉熵损失函数定义

那么现在只需要一行就可以实现交叉熵函数了：

def cross_entropy(y_hat, y):return -torch.log(y_hat[range(len(y_hat)), y])

注意，原来的交叉熵损失函数实际上是：
$$l(y,\hat{y})=-\sum _{j=1}^q{y}_{j}log{\hat{y}}_j$$
其中，q是独热编码的长度，那么容易知道，那个求和符号其实没啥用，因为利用独热编码的话，除了中标的那一项，其他的y中元素全是0。所以引变为代码中的：
$$l(y,\hat{y})=-log{\hat{y}}_j$$
验证：

print(cross_entropy(y_hat, y))

结果：

tensor([2.3026, 0.6931])

分类精度

给定预测概率分布y_hat，我们要给出硬预测时，通常选择预测概率最高的类。
当预测和标签分类y一致时，就是正确的。分类精度即正确预测数量与总预测数量之比。虽然直接优化精度可能很难（精度计算不可导），但我们总是要关注他。
我们可以进行下面的操作：
若y_hat是矩阵，假定第二维度存储每个类的预测分数，我们就可以使用argmax来获得每行的最大元素的索引，用来获得预测的类别。然后和真实的y比较。（注意，由于等式运算符号"=="对数据类型很敏感，因此我们需要将数据类型转换为一致的。）结果会是一个包含0和1的张量，求和就可以得到正确预测的数量了。

def accuracy(y_hat, y):  #@save"""计算预测正确的数量"""if len(y_hat.shape) > 1 and y_hat.shape[1] > 1:  # 判断是矩阵y_hat = y_hat.argmax(axis=1)cmp = y_hat.astype(y.dtype) == yreturn float(cmp.astype(y.dtype).sum())

我们将继续使用之前定义的变量y_hat和y分别作为预测的概率分布和标签。 可以看到，第一个样本的预测类别是2（该行的最大元素为0.6，索引为2），这与实际标签0不一致。 第二个样本的预测类别是2（该行的最大元素为0.5，索引为2），这与实际标签2一致。 因此，这两个样本的分类精度率为0.5。

print(accuracy(y_hat, y) / len(y))

结果：

0.5

同样，对于任意数据迭代器data_iter可访问的数据集，我们可以评估在任意模型net的精度。
我们先定义一个实用程序类Accumulator用于对多个变量进行累加：

class Accumulator:  #@save"""在n个变量上累加"""def __init__(self, n):self.data = [0.0] * ndef add(self, *args):  # *号会拆解为元组self.data = [a + float(b) for a, b in zip(self.data, args)]  # zip就是把两元组组合起来def reset(self):self.data = [0.0] * len(self.data)def __getitem__(self, idx):return self.data[idx]

接着我们定义evaluate_accuracy函数用于计算在指定数据集上模型的精度：

def evaluate_accuracy(net, data_iter):  #@save"""计算在指定数据集上模型的精度"""if isinstance(net, torch.nn.Module):net.eval()  # 将模型设置为评估模式metric = Accumulator(2)  # 正确预测数、预测总数with torch.no_grad():for X, y in data_iter:metric.add(accuracy(net(X), y), y.numel())return metric[0] / metric[1]

在上面的evaluate_accuracy函数中，我们在Accumulator实例中创建了2个变量，分别用于存储正确预测的数量和预测的总数量。当我们遍历数据集时，两者都将随着时间的推移而累加。

训练

首先，我们定义一个函数来训练一个迭代周期。（注意：updater是更新模型参数的常用函数，它接受批量大小作为参数。它可以是d2l.sgd函数，也可以是框架的内置优化函数。）

def train_epoch_ch3(net, train_iter, loss, updater):  #@save"""训练模型一个迭代周期"""# 将模型设置为训练模式if isinstance(net, torch.nn.Module):net.train()# 训练损失总和、训练准确度总和、样本数metric = Accumulator(3)for X, y in train_iter:# 计算梯度并更新参数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]

在展示训练函数实现前，定义一个在动画中绘制数据的应用程序类Animator（会用就行）：

class Animator:  #@save"""在动画中绘制数据"""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, fmtsdef add(self, x, y):# 向图表中添加多个数据点if not hasattr(y, "__len__"):y = [y]n = len(y)if not hasattr(x, "__len__"):x = [x] * nif 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)d2l.plt.draw()d2l.plt.pause(0.001)display.clear_output(wait=True)

接下来，实现一个训练函数，它会在train_iter访问到的训练数据集上训练一个模型net。该训练函数会运行多个迭代周期。在每个迭代周期结束时，利用test_iter访问到的测试数据集对模型进行评估。我们利用Animator类来可视化训练进度。

def train_ch3(net, train_iter, test_iter, loss, num_epochs, updater):  #@save"""训练模型"""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_ch3(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语句表示断言，表达式为False时会触发AssertionError异常assert train_loss < 0.5, train_lossassert train_acc <= 1 and train_acc > 0.7, train_accassert test_acc <= 1 and test_acc > 0.7, test_acc

我们使用之前定义的小批量随机梯度下降来优化模型的损失函数，设学习率为0.1：

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

现在，训练10个迭代周期：

num_epochs = 10
train_ch3(net, train_iter, test_iter, cross_entropy, num_epochs, updater)

这边是可以跑出动图的，如果跑不出来动态的效果，解决方案：
File ——> Settings ——> Tools ——> Python Scientific ——> 取消勾选 Show plots in toolwindow
（电脑快跑炸了）

预测

训练已经完成，我们的模型可以进行分类预测了，给定一系列图像，我们将比较它们的实际标签（文本输出的第一行）和模型预测（文本输出的第二行）。

def predict_ch3(net, test_iter, n=6):  #@save"""预测标签（定义见第3章）"""for X, y in test_iter:breaktrues = d2l.get_fashion_mnist_labels(y)preds = d2l.get_fashion_mnist_labels(net(X).argmax(axis=1))titles = [true +'\n' + pred for true, pred in zip(trues, preds)]d2l.show_images(X[0:n].reshape((n, 28, 28)), 2, n, titles=titles[0:n])predict_ch3(net, test_iter)