[npnet] 逻辑回归
逻辑回归 (logistic regression) 是广义线性模型 (generalized linear model) 的一种，对应的是伯努利分布 (Bernoulli distribution)。如果我们在广义线性模型中使用的是多项分布 (multinomial distribution)，那么我们将得到逻辑回归的推广——多项逻辑回归 (multinomial logistic regression)。我们在这里用逻辑回归来代表逻辑回归和多项逻辑回归。
import numpy as np
import matplotlib.pyplot as plt
%run import_npnet.py
npnet = import_npnet(2)
我之前已经写过三篇文章介绍广义线性模型以及多项逻辑回归了，所以我并不打算重复介绍理论部分。关于广义线性模型的理论介绍可以在我的一篇英文博文 Generalized Linear Model 找到。在另一篇英文博文 Generalized Linear Model (Examples)，我给出了两个广义线性模型的例子：线性模型以及逻辑回归。我在一篇中文博客 Digit recognition, Softmax 使用多项逻辑回归进行手写数字识别。
逻辑回归解决的是分类问题。比如手写数字识别的任务是识别图片中的数字，或句话说就是将图片分成0到9共10类。抛开理论来说，逻辑回归是一个使用特殊损失函数的线性模型。还是以手写数字识别任务为例。每一个输入是一个 \(28 \times 28\) 的图片，也就是一个长度为784的向量。因为有10个类别（0-9），所以我们可以使用一个输入规模784输出规模为10的线性模型。当我们给这个线性模型一个输入时，我们会得到一个长度为10的向量。我们规定向量的位置是从0开始数起的，那么这个长度为10的向量中数字最大的位置便作为此向量代表的类别。比如向量 \(x = (-3.4, 4, 2.9, -23, 3.9, 12, -9.4, 9.9, 8.3, -5.3)\) 最大的数字是12，在位置5上，所以这个向量代表数字5。
我们还能简单地使用均方误差 (mean square error) 吗？假设上述的向量 \(x\) 事实上对应的是1，那么我们能用
来衡量这个局部误差吗？能，但是这是一个很差的衡量标准。问题出现在找向量最大数字所在的位置不是一个可导函数，所以我们不能对这种标准使用梯度下降。如果我们坚持用均方误差怎么办。我们先来回顾一下 softmax 函数 \(\sigma\)。假设 \(x = (x_1, \dots, x_n)\), 那么
显然 \(\sigma(x)\) 的所有坐标都是大于0，而且所有坐标之和为1。所以我们可以将 \(\sigma(x)\) 的每一个坐标看作一个概率，而 \(\sigma(x)\) 是一个概率分布。一个坚持使用均方误差的可行解决方案是首先将上面那个向量\(x\)用 softmax 函数 \(\sigma\) 变成一个概率分布，然后用
来衡量这个局部误差。这里 \((0, 1, 0, 0, 0, 0, 0, 0, 0, 0)\) 是对数字1的独热编码 (one-hot encoding)。
我们看到 \(\sigma(x)\) 是一个概率分布，而独热编码也是一个概率分布（所有坐标非负，而且和是1），我们可以用交叉熵 (cross entropy)。交叉熵是专门用来比较两个概率分布的。假设我们要比较（多项）概率分布 \(p\) 和 \(q\)，交叉熵的公式是
这里我将 \(p\) 看作预测，而 \(q\) 是观测值。在我们的应用中，\(p = \sigma(x), q = e_k\)，其中 \(e_k\) 代表类别 \(k\) 的独热编码。假设 \(x = (x_1, \dots, x_n)\)，那么
这里给出了一个局部误差公式 \(CE(x, e_k)\)。我们用符号 \(s = (s_1, \dots, s_n)\) 代表 \(\sigma(x)\) 的计算结果，也就是说
那么上述局部误差公式可以简化成 \(CE(x, e_k) = -\log s_k\)。
接下来我们要求出交叉熵对于输入 \(x\) 的偏导数。假设 \(j \ne k\)，
另一方面，
从上面的公式可以看到，无论是计算 \(CE(x, e_k)\) 或是偏导数，我们只需要用到 \(\sigma(x)\) 的计算结果 \(s\)。因此，在设计交叉熵标准时，我们要保存 \(\sigma(x)\) 的计算结果。与往常一样，CrossEntroy.forward 接受两个输入 output 和 target 。其中 output 是多个输入组合而成的矩阵，每一行代表一个输入 \(x\)。而 target 的每个元素代表 output 相对应行的分类索引（从0开始数起）。
# import
class CrossEntropy(npnet.Criterion):
def forward(self, output, target):
# 先计算 σ(x)
s = np.exp(output - np.max(output, axis=1, keepdims=True))
s = s / np.sum(s, axis=1, keepdims=True)
self.cached = s
self.target = target.ravel()
# 计算交叉熵 -log s_k
s = s[np.arange(self.target.size), self.target]
s = np.log(s)
return -np.sum(s) / output.shape[0]
def backward(self, grad=1.0):
# 计算偏导数
target, cached = self.target, self.cached
cached[np.arange(target.size), target] -= 1
cached /= cached.shape[0]
cached *= grad
return cached
Fashion MNIST
手写数字识别是机器学习经典的例子，关于 MNIST 数据集的教程已经烂大街了。我自己就在这个博客中写了两篇相关的文章：Digit recognition, Softmax 和 Digit recognition, CNN。最近看到一个叫做 Fashion-MNIST 的数据集，它跟经典的 MNIST 数据集很像，但是更有趣，我也没有玩过。作为例子用在这个系列笔记我感觉是个不错的选择。
import os
from urllib.request import urlopen
def download():
url = 'http://fashion-mnist.s3-website.eu-central-1.amazonaws.com/'
folder = 'data'
files = [
'train-images-idx3-ubyte.gz',
'train-labels-idx1-ubyte.gz',
't10k-images-idx3-ubyte.gz',
't10k-labels-idx1-ubyte.gz'
]
if not os.path.exists(folder):
os.makedirs(folder)
for file in files:
path = os.path.join(folder, file)
if os.path.exists(path):
continue
file_url = url + file
with open(path, 'wb') as f:
f.write(urlopen(file_url).read())
download()
数据集已经下载下来了，我们要将数据集转化成 numpy 格式。
# import
# credit: https://github.com/zalandoresearch/fashion-mnist/blob/master/utils/mnist_reader.py
# kind = 'train' or 't10k'
def load_fashion_mnist(kind='train'):
import os
import gzip
import numpy as np
labels_npy = f'data/{kind}-labels-idx1-ubyte.npy'
images_npy = f'data/{kind}-images-idx3-ubyte.npy'
try:
images = np.load(images_npy)
labels = np.load(labels_npy)
except Exception:
labels_path = f'data/{kind}-labels-idx1-ubyte.gz'
images_path = f'data/{kind}-images-idx3-ubyte.gz'
with gzip.open(labels_path, 'rb') as lbpath:
labels = np.frombuffer(lbpath.read(), dtype=np.uint8, offset=8)
with gzip.open(images_path, 'rb') as imgpath:
images = np.frombuffer(imgpath.read(), dtype=np.uint8,
offset=16).reshape(len(labels), 784)
np.save(images_npy, images)
np.save(labels_npy, labels)
return images, labels
fashion_mnist_labels = [
'T-shirt/top','Trouser', 'Pullover', 'Dress', 'Coat',
'Sandal', 'Shirt', 'Sneaker', 'Bag', 'Ankle boot'
]
好了，让我们看看前25个训练数据。
images, labels = load_fashion_mnist(kind='train')
r = c = 5
fig, ax = plt.subplots(r, c)
for i in range(r):
for j in range(c):
k = i * c + j
img = images[k]
img = img.reshape((28, 28))
lbl = labels[k]
lbl = fashion_mnist_labels[lbl]
ax[i, j].imshow(img, cmap='gray_r')
ax[i, j].set_title(lbl)
ax[i, j].axis('off')
plt.subplots_adjust(hspace=0.7)
plt.show()
我们现在要用逻辑回归去解决 Fashion-MNIST 分类。我们将要建一个输入规模784输出规模10的线性模型，使用 CrossEntropy 作为标准，并用梯度下降法 SGD 来训练模型。训练集一共有60000条数据，我们训练的时候不会一次输入60000条数据，因为这样的话花费的时间和内存都特别高。我们会将数据集分成多个批次，每个批次的数据量都为100。这种训练方式叫做批梯度下降 (mini-batch gradient descent)。这是两种极端的中和：每次使用所有数据的梯度下降 (gradient descent)，以及每次使用一个数据的随机梯度下降 (stochastic gradient descent)。
我们将测试的数据集也加载进来：
test_images, test_labels = load_fashion_mnist(kind='t10k')
在训练过程中，一个 epoch （实在找不到好的翻译）是指将所有数据都完整地使用一次；而一次迭代 (iteration) 是指一次更新参数的过程。在下面的训练代码中，我在每个 epoch 之前都先随机打乱数据集的顺序，然后将数据分成若干份进行迭代。在每个 epoch 结束之后我们都会保存一下训练与测试各自的损失以及准确率。
因为训练模型的代码具有复用性，我们就将训练的过程写成了一个函数，这样以后也不用重复这段代码了。
# import
def mnist_loss_and_accuracy(model, criterion,
images, labels, batch_size=None):
if batch_size is None:
batch_size = len(images)
batch_num = len(images) // batch_size
if len(images) % batch_size != 0:
batch_num += 1
loss = accuracy = 0
for i in range(batch_num):
input_batch = images[i * batch_size: (i + 1) * batch_size]
target_batch = labels[i * batch_size: (i + 1) * batch_size]
output = model.forward(input_batch)
loss += criterion.forward(output, target_batch)
output = np.argmax(output, axis=1)
accuracy += np.sum(output == target_batch)
loss /= batch_num
accuracy /= len(images)
return loss, accuracy
def train_fashion_mnist(model, criterion, optim,
train_images, train_labels,
batch_size=100, epoch=5,
profile=False,
test_images=None, test_labels=None):
batch_num = len(train_images) // batch_size
if len(train_images) % batch_size != 0:
batch_num += 1
# 存储训练损失以及测试损失、测试准确率
train_loss = []
train_accuracy = []
test_loss = []
test_accuracy = []
# 开始训练
for ep in range(epoch):
# 随机化训练数据的顺序
idx = np.arange(len(train_images))
np.random.shuffle(idx)
images = train_images[idx]
labels = train_labels[idx]
for i in range(batch_num):
# 一次迭代
batch = images[i * batch_size: (i + 1) * batch_size]
target = labels[i * batch_size: (i + 1) * batch_size]
output = model.forward(batch)
criterion.forward(output, target)
grad = criterion.backward()
optim.zero_grad()
model.backward(grad)
optim.step()
# 是否计算损失与准确率
if not profile:
continue
# 保存训练以及测试各自的损失、准确率
loss, accuracy = mnist_loss_and_accuracy(
model, criterion, train_images, train_labels, batch_size
)
train_loss.append(loss)
train_accuracy.append(accuracy)
if type(test_images) is not np.ndarray:
continue
if type(test_labels) is not np.ndarray:
continue
loss, accuracy = mnist_loss_and_accuracy(
model, criterion, test_images, test_labels, batch_size
)
test_loss.append(loss)
test_accuracy.append(accuracy)
return train_loss, train_accuracy, test_loss, test_accuracy
我们先建立模型，然后使用 train_fashion_mnist 训练模型。
# 超参数 (hyperparameters)
lr = 1e-5
batch_size = 100
epoch = 50
# 模型、标准、优化器
model = npnet.Linear(784, 10)
criterion = CrossEntropy()
optim = npnet.SGD(model, lr=lr)
train_loss, train_accuracy, test_loss, test_accuracy = train_fashion_mnist(
model=model, criterion=criterion, optim=optim,
train_images=images, train_labels=labels,
batch_size=batch_size, epoch=epoch,
test_images=test_images, test_labels=test_labels,
profile=True
)
接下来看看训练过程中训练损失以及测试损失是如何随着 epoch 数量而变化的。训练损失与测试损失很相近，这是一件好事，说明我们的模型没有过拟合 (overfitting)。
x = np.arange(len(train_loss))
plt.plot(x, train_loss, 'b', label='train')
plt.plot(x, test_loss, 'g', label='test')
plt.legend()
plt.title('Model Loss')
plt.xlabel('epoch')
plt.ylabel('loss')
plt.show()
我们还可以看测试准确率的变化：
plt.plot(x, train_accuracy, 'b', label='train')
plt.plot(x, test_accuracy, 'g', label='test')
plt.legend(loc='lower right')
plt.title('Model Accuracy')
plt.xlabel('epoch')
plt.ylabel('accuracy')
plt.show()
我们看到模型的准确率大概在75.0%到80.0%之间。训练损失一开始的时候下降得很快，但是后面下降速度就减缓了。这是正常的，因为一开始的时候模型很差，所以一个 epoch 能对模型有比较大的改进。但是后来模型逐渐收敛之后，一个 epoch 起到的改进作用就越来越小了。下一篇文章我们将会介绍梯度下降的变种。这些变种一般会有更快的收敛速度。
最后我们来看一下不同学习率 (learning rate) 的训练效果。
batch_size = 100
epoch = 50
lrs = [1e-4, 1e-5, 1e-6]
for lr in lrs:
model = npnet.Linear(784, 10)
criterion = CrossEntropy()
optim = npnet.SGD(model, lr=lr)
loss, *_ = train_fashion_mnist(
model=model, criterion=criterion, optim=optim,
train_images=images, train_labels=labels,
profile=True, batch_size=batch_size, epoch=epoch
)
plt.plot(np.arange(len(loss)), loss, label=f'lr={lr:.0e}')
plt.title('Training loss with respect to learning rate')
plt.xlabel('epoch')
plt.ylabel('loss')
plt.legend()
plt.show()
当学习率过大的时候训练误差是随机上下浮动没有什么规律可言的，当学习率太小了训练误差收敛的速度会很慢。因此选择一个合适的学习率在机器学习中非常重要。事实上，在机器学习中如学习率等的超参数 (hyperparameter) 的选取是一门高深的学问，而且模型的训练很大一部分时间就花在超参数的选取中。所以在业界，从事机器学习的人员被称为调参师（手动滑稽）。
