问题描述

请设计下述算法，求解 MNIST 数据集上的 Logistic Regression 问题：

  1. 梯度下降法；
  2. 随机梯度法。

对于每种算法，请给出每步计算结果与最优解的距离以及每步计算结果在测试集上所对应的分类精度。此外，请讨论随机梯度法中 Mini Batch 大小对计算结果的影响。可参考：Classifying MNIST digits using Logistic Regression — DeepLearning 0.1 documentation

本文地址：https://www.jeddd.com/article/use-python-to-solve-softmax-regression-on-mnist-by-sgd.html

思路

MNIST 是一个经典的手写数字数据集，其原始数据集包含 60000 个训练数据和 10000 个测试数据。每个数据由样本特征和标签组成：特征为 28×28 像素的灰度图（784 维向量），表示一个手写数字；标签为 0~9 中的一个数字。我们的目标是构建一个模型，然后在训练数据集上使用梯度下降法（GD）和随机梯度法（SGD）训练模型。模型训练好后，用测试数据集去评估模型的准确率。

实验环境为 Python 3.7 和 Jupyter Notebook，用到的外部包只有 numpy 和 matplotlib，以及用于高性能计算的 numba 库。没有使用任何机器学习框架，易于理解。

完整代码：https://github.com/Jed-Z/optimization-theory/blob/master/Opt2_MNIST_Softmax_LR.ipynb

建模

我们的目标是建立一个模型， 模型的输入是 784 维的特征向量 $x$，表示一个手写数字的灰度图像；模型的输出是一个 10 维的向量 $y$，表示该预测图像是 0~9 中某个数字的概率。输出概率最大的那个数字就是预测值，需将其与标签进行比对，从而评估模型的准确率。这是一个线性模型：

$$ y=Wx+b, $$

其中特征向量 $x\in \mathbb{R^{784}}$，预测结果 $y\in \mathbb{R^{10}}$，权重参数矩阵 $W\in \mathbb{R^{10\times 784}}$，偏置参数向量 $b\in \mathbb{R^{10}}$，预测结果 $y\in \mathbb{R^{10}}$。

在 Python 中，首先定义一些常量，用于表示各种向量和矩阵的维度：

# 定义常量
x_dim = 28 * 28  # 输入维度
y_dim = 10       # 输出维度
W_dim = (y_dim, x_dim)  # 权重矩阵参数的维度
b_dim = y_dim           # 偏置向量参数的维度

读取数据

从 http://yann.lecun.com/exdb/mnist/ 上下载四个 gz 文件放入MNIST_data目录下，然后分别将四个文件读入四个 np.array 中。

# 读取数据
import gzip
def read_idx3(filename):
    with gzip.open(filename, 'rb') as fo:
        buf = fo.read()
        index = 0
        header = np.frombuffer(buf, '>i', 4, index)
        index += header.size * header.itemsize
        data = np.frombuffer(buf, '>B', header[1] * header[2] * header[3], index).reshape(header[1], -1)
        return data
    
def read_idx1(filename):
    with gzip.open(filename, 'rb') as fo:
        buf = fo.read()
        index = 0
        header = np.frombuffer(buf, '>i', 2, index)
        index += header.size * header.itemsize
        data = np.frombuffer(buf, '>B', header[1], index)
        return data

data_path = 'MNIST_data/'

X_train = read_idx3(data_path + 'train-images-idx3-ubyte.gz')  # 训练数据集的样本特征
y_train = read_idx1(data_path + 'train-labels-idx1-ubyte.gz')  # 训练数据集的标签
X_test = read_idx3(data_path + 't10k-images-idx3-ubyte.gz')  # 测试数据集的样本特征
y_test = read_idx1(data_path + 't10k-labels-idx1-ubyte.gz')  # 测试数据集的标签

可以观察一下读入数据集的格式：

print(X_train.shape, y_train.shape)
print(X_test.shape, y_test.shape)
(60000, 784) (60000,)
(10000, 784) (10000,)

从训练数据集中随意抽取一个样本，观察它的图像和标签：

# 随机选一个数据观察
example_idx = np.random.randint(len(X_train))  # 随机选取的样本编号
# example_idx = 4221
plt.imshow(np.reshape(X_train[example_idx], (28, 28)), cmap='gray')
plt.show()
print('Index: {}; Label: {}'.format(example_idx, y_train[example_idx]))

本文地址：https://www.jeddd.com/article/use-python-to-solve-softmax-regression-on-mnist-by-sgd.html

Softmax 函数

狭义的逻辑回归（Logistic Regression）是解决二分类问题的，常使用 Sigmoid 函数作为激活函数。而 Softmax 回归模型是 Logistic 回归模型在多分类问题上的推广，常用的激活函数是 Softmax 函数。在 MNIST 数据集上，标签为 0~9 共 10 个类别，因此这是一个多分类问题。（注：当 Softmax 作用于二维向量时，它将退化为 Sigmoid。）

Softmax 函数的定义如下：

$$ S_i(x)=\dfrac{e^{x_i}}{\sum_{j}e^{x_j}} $$

对一个向量 $x$，Softmax 函数把它的元素映射到 $(0, 1)$ 区间内，且映射后的向量的所有元素之和为 1，这恰好符合概率的定义，因此可以将 Softmax 函数的结果看作每个元素值出现的概率。

在 Python 中定义 Softmax 函数，其输入参数为一个向量。注意该函数中首先将所有元素都减去了最大元素np.max(x)，这样做是为了防止大元素经指数变换后溢出而使结果中出现 nan。

def softmax(x):
    x_shifted = x - np.max(x)  # 防止溢出为nan
    return np.exp(x_shifted) / np.exp(x_shifted).sum()

损失函数及其梯度

理解损失函数

使用对数损失函数。上述模型在一个样本点 $(x_i, y_i)\in\mathbb{R^{784}}\times \mathbb{R}$ 处的损失函数定义如下：

$$ L_i(W, b)=-\log\left(S_{y_i}(Wx_i+b)\right) $$

对于这个损失函数可以这样理解：现在有一个样本，其特征为 $x_i$，将该特征输入我们的模型，输出的结果（预测值）为 $Wx_i+b$，这是一个 10 维的向量，现将这个 10 维向量通过 Softmax 函数映射为另一个 10 维向量，记为 $S(x)$，映射后的这个向量 $S(x)$ 的元素具有概率意义，其每个元素的值代表预测为 0~9 中某个数字的概率。预测值等于标签的概率是多少呢？我们知道标签是 $y_i$，它是 $x_i$ 对应的真实结果；而预测结果和真实结果相等的概率就是 $S(x)$ 中与 $y_i$ 对应元素的值，即 $S_{y_i}(x)$。根据最大似然，我们希望预测结果和真实结果相等的概率最大。因而损失函数就定义为上述概率的负对数。取对数的原因是：当损失函数作用于多个样本时，可以将概率的乘积转换为其对数的和；取负号的原因是可以把一个极大化问题转换为极小化问题，从而可以使用下降方法来求解。（以上仅为个人理解，如有误请指正。）

举个例子帮助理解。比如现在有一个编号为 $k$ 的样本，其特征 $x_k$ 是一个784 维的向量，假设它的标签 $y_k=3$，也就是这个手写数字表示的是“3”。现将 $x_k$ 输入我们的模型，记为 $\text{y_pred}_k=Wx_k+b$，我们将其输入 Softmax 函数，得到 $S(\text{y_pred}_k)$，假设 $S(\text{y_pred}_k)=(0.01, 0.01, 0.01,0.91, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01)^T$，发现预测值 $S(\text{y_pred}_k)$ 等于标签 3 的概率为 91%，为其它每个数字的概率都是 1%。此时我们取 91%。最后，取其负对数。

需要注意，这里的损失函数 $L(W, b)$ 是在一个样本点处计算的，如果有多个样本点 $j$，只需对每个样本点处的损失求和即可（上面说了，概率乘积的对数即为概率对数的加和），总损失为 $\sum_j L_j(W, b)$。

损失函数可以写成下面的 Python 函数。不过损失函数本身并不直接起作用，真正起作用的是损失函数的梯度。因此，损失函数的代码在我的源码中已被注释掉，仅用于帮助理解。

def loss(W, b, x, y):
    y_pred = W @ x + b
    return -np.log(softmax(y_pred)[y])  # 预测值与标签相同的概率

损失函数的梯度

损失函数梯度的计算需要一些数理技巧，详细推导见这篇博客，这里直接使用其结果。注意下面的 $i$ 和 $j$ 仅代表矩阵或向量的下标，$x$ 和 $y$ 分别代表一个样本的特征向量及其标签，$S_i=S_i(Wx+b)$ 为预测值经 Softmax 映射后的概率。

$$ \begin{equation} \nabla W_{ij}= \begin{cases} (S_i-1)x_j,& y=i \\ S_i x_j,& y\ne i \end{cases} \end{equation} $$

$$ \begin{equation} \nabla b_i= \begin{cases} (S_i-1),& y=i \\ S_i,& y\ne i \end{cases} \end{equation} $$

在 Python 中实现：

@numba.jit
def loss_gradient(W, b, x, y):
    """
    W和b为当前的参数，x为一个样本，y为x对应的标签。
    """
    W_gradient = np.zeros(W.shape)
    b_gradient = np.zeros(b.shape)
    sm = softmax(W @ x + b)
    for i in range(W.shape[0]):
        for j in range(W.shape[1]):
            if y == i:
                W_gradient[i][j] = (sm[i] - 1) * x[j]
            else:
                W_gradient[i][j] = sm[i] * x[j]
                
    for i in range(b.shape[0]):
        if y == i:
            b_gradient[i] = sm[i] - 1
        else:
            b_gradient[i] = sm[i]
    return W_gradient, b_gradient

装饰器@numba.jit是用于加速计算的，请参考Numba: A High Performance Python Compiler。

梯度下降与随机梯度下降

随机梯度与批量梯度

梯度下降法（Gradient descent, GD），或称“全批量梯度下降法”，其每次更新参数时所用的梯度是整个样本的梯度，也就是所有样本点处梯度的平均值。在 MNIST 数据集上，一次全批量梯度下降需要计算 60000 次梯度。

随机梯度下降（Stochastic gradient descent, SGD）在每次更新时只使用 1 个样本。这样，尽管不是每次迭代得到的损失函数都向着全局最优方向， 但是大的整体的方向是向全局最优解的，最终的结果往往是在全局最优解附近。

介于二者之间存在一种称为“小批量（mini-batch）梯度下降”（MSGD）的方法，它将整个样本分为多个 batch，每个 batch 中有batch_size个样本。这种方法中，每次更新使用的是一个 batch 内的所有样本的平均梯度。

当batch_size等于样本数量时，即为全批量梯度下降 GD；当batch_size等于 1 时，即为狭义的随机梯度下降 SGD。因此可以认为 MSGD 是 GD 和 SGD 的一般情况。

概念辨析：Epoch、Batch、Iteration

  • Epoch：使用整个训练集的全部数据对模型进行了一次完整的训练，称为“一代训练”；
  • Batch：训练集被分成若干批，每一批数据中包含 batch_size 个样本，称为“一批数据”；
  • Iteration：使用一个 batch 对模型进行一次参数更新，称为“一次训练”。

重点要区分 epoch 和 iteration 的区别，这些概念是后面作图的基础。

小批量梯度下降

通过设置不同的batch_size，MSGD 可以转化为 GD 或 SGD，因此我们只需要实现一个 MSGD 的函数即可。为了方便我们测试不同的步长、batch_size、epochs_num 等参数对模型效果的影响，考虑将 MSGD 写为一个函数msgd。

在msgd中，首先要将全局变量W和b清零，这样是为了避免多次运行msgd时上次运行的结果对本次运行造成影响。为了实现小批量，我们使用列表解析和列表分片的语法将训练数据集X_train和y_train分割成 batch，即下面代码中的batches_x和batches_y。下面有三层嵌套的 for 循环，最外层的循环变量为epoch，表示一代训练；中间一层的循环变量为i，表示一次 iteration，即一次参数更新；最内层的j是在 batch 内求平均梯度用的。每次 iteration 结束后，立即用score测试当前的 $W$ 和 $b$ （score函数的定义见下方“评估模型”），并将结果放入列表plot_scores中；此外还要把W和b分别放入列表plot_iter_W和plot_iter_b中。注意加入参数时可能需要使用copy()。

W = np.zeros(W_dim)  # W的初值
b = np.zeros(b_dim)  # b的初值

# mini-batch随机梯度下降
def msgd(alpha, batch_size, epochsnum):
    global W, b
    W = np.zeros(W_dim)  # 清零，避免上次运行结果的影响
    b = np.zeros(b_dim)  # 清零
    
    # batches_x和batches_y分别是将原样本和标签按batch分批后的结果
    batches_x = np.array([X_train[i:i+batch_size] for i in range(0, X_train.shape[0], batch_size)])
    batches_y = np.array([y_train[i:i+batch_size] for i in range(0, y_train.shape[0], batch_size)])

    batches_num = batches_x.shape[0]  # batch的总数量
    plot_scores = []  # 记录每次迭代的分类精度
    plot_iter_W = []  # 记录每次迭代的参数W
    plot_iter_b = []  # 记录每次迭代的参数b

    print('Start training...')

    time_start = time.time()  # 开始计时

    for epoch in range(epochs_num):
        for i in range(batches_num):
            batch_x = batches_x[i]
            batch_y = batches_y[i]
            W_gradient = np.zeros(W_dim)
            b_gradient = np.zeros(b_dim)
            for j in range(batch_size):
                gradients = loss_gradient(W, b, batch_x[j], batch_y[j])
                W_gradient += gradients[0]
                b_gradient += gradients[1]
            W_gradient /= batch_size  # 该batch的平均梯度
            b_gradient /= batch_size  # 该batch的平均梯度

            W -= alpha * W_gradient
            b -= alpha * b_gradient
            
            current_score = score(W, b, X_test, y_test)
            plot_scores.append(current_score)  # 记录每次迭代的分类精度
            plot_iter_W.append(W.copy())
            plot_iter_b.append(b.copy())

    time_end = time.time()  # 结束计时
    return (time_end - time_start), plot_scores, plot_iter_W, plot_iter_b

本文地址：https://www.jeddd.com/article/use-python-to-solve-softmax-regression-on-mnist-by-sgd.html

评估模型

检验模型在测试数据集上的预测分类精度（准确率）。定义一个测试函数score，该函数在测试集X_test和y_test上运行由 $W$ 和 $b$ 确定的模型，返回预测分类精度，也就是预测准确率。

def score(W, b, X_test, y_test):
    """
    W和b为待评估的模型参数，x_test和y_test为测试数据集的特征和标签
    """
    assert len(X_test) == len(y_test)
    
    tests_num = len(X_test)  # 测试集的大小
    results = []             # 保存每个测试的结果，0为错误预测；1为正确预测
    for i in range(tests_num):
        y_pred = W @ X_test[i] + b  # 预测结果
        sm = softmax(y_pred)
        results.append(sm.argmax() == y_test[i])
    return np.mean(results)  # 返回分类精度

首先设定参数（调参)，然后调用msgd，训练模型并绘图：

alpha = 0.001    # 步长
batch_size = 60000  # mini-batch的大小
epochs_num = 10000  # 遍历完整样本的次数

result = msgd(alpha, batch_size, epochs_num)
time_cost = result[0]
plot_scores = result[1]
plot_iter_W = result[2]
plot_iter_b = result[3]
print('Alpha:{}; Batch_size:{}; Epochs_num:{}; Time cost:{:.2f} s.'.format(alpha, batch_size, epochs_num, time_cost))

for i in range(len(plot_iter_W)):
    plot_iter_W[i] = np.linalg.norm(plot_iter_W[i] - W)
    plot_iter_b[i] = np.linalg.norm(plot_iter_W[i] - b)

# 作图
plt.title('Classification accuracy with the number of iterations')
plt.xlabel('Iterations')
plt.ylabel('Accuracy')
plt.plot(plot_scores)
plt.grid()
# plt.savefig('1.svg')
plt.show()

plt.title('Distance')
plt.xlabel('Iterations')
plt.ylabel('Distance')
plt.plot(plot_iter_W)
plt.plot(plot_iter_b)
plt.grid()
# plt.savefig('2.svg')
plt.show()

结果分析

下面用不同的参数组合 $(\text{batch_size}, \alpha, \text{epochs_num})$ 训练模型。

SGD：batch_size = 1

参数：batch_size = 1，alpha = 1e-5，epochs_num = 1。将会迭代（参数更新） 60000 次。最终在测试集上的分类精度为 85.87%。

分类精度随迭代次数的变化曲线：

参数 $W$ 和参数 $b$ 距最优解的距离的变化曲线：

小批量：batch_size = 10

参数：batch_size = 10，alpha = 1e-5，epochs_num = 1。将会迭代（参数更新） 6000 次。最终在测试集上的分类精度为 85.92%。

分类精度随迭代次数的变化曲线：

参数 $W$ 和参数 $b$ 距最优解的距离的变化曲线：

小批量：batch_size = 100

参数：batch_size = 100，alpha = 1e-5，epochs_num = 2。将会迭代（参数更新） 1200 次。最终在测试集上的分类精度为 91.34%。

分类精度随迭代次数的变化曲线：

参数 $W$ 和参数 $b$ 距最优解的距离的变化曲线：

小批量：batch_size = 1000

参数：batch_size = 1000，alpha = 1e-5，epochs_num = 10。将会迭代（参数更新） 600 次。最终在测试集上的分类精度为 91.95%。

分类精度随迭代次数的变化曲线：

参数 $W$ 和参数 $b$ 距最优解的距离的变化曲线：

全批量：batch_size = 60000

参数：batch_size = 60000，alpha = 1e-5，epochs_num = 500。将会迭代（参数更新） 500 次。最终在测试集上的分类精度为 91.82%。

分类精度随迭代次数的变化曲线：

参数 $W$ 和参数 $b$ 距最优解的距离的变化曲线：

分析与总结

将以上结果总结在一张表格内：

mini-batch步长 $\alpha$epochs 数iterations 数精度训练时间
1$1\times 10^{-5}$16000085.87%12890.55 秒*
10$1\times 10^{-5}$1600085.92%1084.3 秒
100$1\times 10^{-5}$2600091.34%212.22 秒
1000$1\times 10^{-5}$1060091.95%132.14 秒
60000$1\times 10^{-5}$50050091.82%1289.94 秒

*：mini-batch 大小为 1 的模型是在 Kaggle kernel 上训练的，其余均是在本地的 i7-7700HQ 上运行的。

事实上，我尝试了更多参数组合，上表列出的这些参数组合只是经过多次实验得出选出来的效果比较好的，或是比较具有代表性的参数组合。调参的技巧性非常强，需要丰富的经验作为支撑。

当 batch_size 为 60000 时，MSGD 即为全批量梯度下降法，我发现取步长 $\alpha=1\times 10^{-5}$ 时效果较好，经过 100 次迭代（全批量中 iterations 数与 epochs 数相等），准确率可达 90.61%。通过改变 mini-batch 的大小测试模型的性能和效果，可以发现，全批量梯度下降只需迭代一代（一个 epoch）就能得到一个令人满意的效果——随着迭代次数（iterations）的增加，$W$ 和 $b$ 都在逐渐向最优解靠近，分类精度（准确率）也稳步提升至 90% 以上，且变化几乎是单调的。而当 mini-batch 较小时，由于 SGD 易受噪声影响，可以看出参数曲线和精度曲线的变化趋势不是单调的，而是有一定程度的振荡，但总体来说精度呈逐渐提高的趋势。理论上，当 mini-batch 比较小时，步长也应该取得比较小，否则会噪声的影响会被放大，从而导致严重的振荡。我上面的实验取 mini-batch 为 1 和 10 时，从图中可以看出精度曲线振荡得非常激烈，我认为这可能是因为步长取得过大了。另外，当 mini-batch 较小时，精度提升的速度明显比全批量更快，经过很少的迭代就可以达到令人满意的精度了。

通常情况下，当 mini-batch 的大小较为适中时，模型能够训练得又快又好，这是我们所希望的。比如上面 mini-batch 为 1000 时，只需迭代 10 代（600 次 iterations），花费约 2 分钟就能得到一个精度约为 92% 的模型了。

References

  1. The Softmax function and its derivative
  2. Getting Started — DeepLearning 0.1 documentation

本文地址：https://www.jeddd.com/article/use-python-to-solve-softmax-regression-on-mnist-by-sgd.html