家乐的深度学习笔记「3」 - 线性回归
线性回归
线性回归的输出是一个连续值，适用于回归问题，如房价预测、气温、销售额等连续值问题。
基本要素
以最基础的房价预测为例，假设价格只取决于面积（平方米）和房龄（年），下面来介绍大多数深度学习模型的基本要素和表示方法。
模型
设房屋面积为x，房龄为x，售出价格为y。我们需要建立基于输入x和x来计算输出y的表达式，也就是模型。
y = xw+ xw + b
其中w1和w2是权重，b是偏差，且均为标量。它们是线性回归模型的参数。y是线性回归对真实价格y的预测或估计。通常允许它们之间有一定误差。
模型训练
构建好模型后，需要通过数据来寻找特定的模型参数值，使模型在数据上的误差尽可能小。这个过程称为模型训练。
训练数据
收集到的一系列真实数据被称作训练数据集或训练集。以房价数据集为例，一栋房屋被称为一个样本，其真实出售价格叫作标签，用来预测标签的两个因素叫做特征，特征用来表征样本的特点。
假设采集的样本数为n，索引为i的样本的特征为x和x，标签为y。对于索引为i的房屋，线性回归模型的房屋价格预测表达式为
y(i) = x(i)w+ x(i)w + b
损失函数
平方函数通常用来衡量价格预测值与真实值之间的误差，其在评估索引为i的样本误差的表达式为
给定训练数据集，这个误差只与模型参数相关，因此将其记为以模型参数为参数的函数。
在机器学习里，将衡量误差的函数称为损失函数。这里使用的平方误差函数也称为平方损失。
通常使用训练数据集中所有样本误差的平均来衡量模型预测的质量，即
在模型训练中，我们希望找出一组模型参数，记为$$w_1^, w_2^, b^*$$，来使训练样本平均损失最小：
2020-03-29：使用训练数据集中所有样本的误差的平均来衡量模型预测的质量，表示刚想通这个概念，刚开始把小批量的概念混在一起了，但是，使用小批量的误差来反向传播更新模型的梯度，而不是整体所有样本。
优化算法
当模型和损失函数形式较为简单时，上面的误差最小化问题的解可以直接用公式表达出来。这类解叫作解析解(analytical solution)。
然而大多数深度学习模型并没有解析解，只能通过优化算法有限次迭代模型参数来尽可能降低损失函数的值。这类解叫作数值解(numerical solution)。
以下为抄书内容：
在求解数值解的优化算法中，小批量随机梯度下降(mini-batch stochastic gradient descent)在深度学习中被广泛使用。其先选取一组模型参数的初始值，如随机选取；接下来对参数进行多次迭代，使每次迭代都可能降低损失函数的值。在每次迭代中，先随机均匀采样一个由固定数目训练数据样本所组成的小批量（mini-batch）\(\mathcal{B}\)，然后求小批量中数据样本的平均损失有关模型参数的导数（梯度），最后用此结果与预先设定的一个正数的乘积作为模型参数在本次迭代的减小量。
在训练本节线性回归模型过程中，模型的每个参数将作如下迭代：
在上式中，\(|\mathcal{B}|\)代表每个小批量中的样本个数（批量大小，batch size），\(\eta\)称作学习率（learning rate）并取正数。需要强调的是，这里的批量大小和学习率的值是人为设定的，并不是通过模型训练学出的，因此叫作超参数（hyperparameter）。我们通常所说的“调参”指的正是调节超参数，例如通过反复试错来找到超参数合适的值。
在少数情况下，超参数也可以通过模型训练学出。
模型预测
模型训练完成后，将模型参数\(w_1, w_2, b\)在优化算法停止时的值记为\(\hat{w}_1, \hat{w}_2, \hat{b}\)，这里得到的不一定是最小化损失函数的最优解，而是对最优解的一个近似。然后，就可以用学出的线性回归模型来估算训练数据集以外任意一个样本的标签了。
这里的估算也叫作模型预测、模型推断或模型测试。
表示方法
下面解释线性回归与神经网络的联系，以及线性回归的矢量计算表达式。
神经网络图
在深度学习中，可以使用神经网络图直观地表现模型结构。
神经网络图隐去了模型参数权重和偏差。
上述神经网络输入分别为x和x，因此输入层的输入个数为2。输入个数也叫特征数或特征向量维度。
上述神经网络中，直接使用神经网络的输出o作为线性回归的输出，即$$\hat{y} = o$$。
由于输入层并不涉及计算，按照惯例，上述神经网络的层数为1。所以，线性回归是一个单层神经网络，输出层中负责计算o的单元又叫神经元。在线性回归中，o的计算依赖于x和x。也就是说，输出层中的神经元和输入层中各个输入完全连接。因此，这里的输出层又叫全连接层(fully-connected layer)，或稠密层(dense layer)。
矢量计算表达式
在模型训练或预测时，常常会同时处理多个数据样本并用到矢量计算。
矢量加法快于标量加法，应该尽可能使用矢量运算以提升效率。
广义上讲，当数据样本数为\(n\)，特征数为\(d\)时，线性回归的矢量计算表达式为
\(\boldsymbol{\hat{y}} = \boldsymbol{X} \boldsymbol{w} + b\)
其中模型输出\(\boldsymbol{\hat{y}} \in \mathbb{R}^{n \times 1}\)， 批量数据样本特征\(\boldsymbol{X} \in \mathbb{R}^{n \times d}\)，权重\(\boldsymbol{w} \in \mathbb{R}^{d \times 1}\)， 偏差\(b \in \mathbb{R}\)。相应地，批量数据样本标签\(\boldsymbol{y} \in \mathbb{R}^{n \times 1}\)。设模型参数\(\boldsymbol{\theta} = [w_1, w_2, b]^\top\)，我们可以重写损失函数为
小批量随机梯度下降的迭代步骤将相应地改写为
其中梯度是损失有关3个为标量的模型参数的偏导数组成的向量：
从零实现
这里不利用高级框架，只用NDArray和autograd实现线性回归。
from mxnet import gpu, nd, autograd
import matplotlib.pyplot as plt
import random
这里可以顺带固定下随机数种子。
random.seed(233)
生成数据集
设训练数据集样本数为1000，输入个数（特征数）为2。给定随机生成的批量样本特征\(\boldsymbol{X} \in \mathbb{R}^{1000 \times 2}\)，使用线性回归模型真实权重\(\boldsymbol{w} = [2, -3.4]^\top\)和偏差\(b = 4.2\)，以及一个随机噪声项\(\epsilon\)来生成标签，其中噪声项服从均值为0、标准差为0.01的正态分布。噪声代表了数据集中无意义的干扰。
因为使用了一台四卡2080Ti的服务器学习，为了尽量不对别人带来困扰，这里使用第四张卡，MXNet不会像tf一样默认占用整张卡的显存，而是一点点占用，这点值得表扬。
ctx = gpu(3)
num_examples = 1000
num_inputs = 2
true_w = [2, -3.4]
true_b = 4.2
features = nd.random.normal(scale=1, shape=(num_examples, num_inputs), ctx=ctx)
labels = true_w[0] * features[:, 0] + true_w[1] * features[:, 1] + true_b
labels += nd.random.normal(scale=0.01, shape=labels.shape, ctx=ctx)
看一眼生成的训练数据集部分特征与标签的样子：
features[:5], labels[:5]
(
[[ 0.34113222 0.76659095]
[ 0.6508758 -1.5585263 ]
[ 1.9786316 -0.8465744 ]
[-0.05428567 0.64376295]
[-1.7093595 0.12594718]]
<NDArray 5x2 @gpu(3)>,
[ 2.2640598 10.818717 11.050849 1.8893591 0.35090753]
<NDArray 5 @gpu(3)>)
生成一个散点图可以直观的观察两者之间的线性关系。
plt.scatter(features[:, 1].asnumpy(), labels.asnumpy(), 1);
改一下生成格式为svg，即矢量图，这样可以省下图片的空间。
如果不保存的话，画图语句末尾的分号可以省去打印该图的类型信息，只显示图。
from IPython import display
display.set_matplotlib_formats('svg')
point_size = 1
plt.scatter(features[:, 1].asnumpy(), labels.asnumpy(), point_size);
plt.savefig('散点图.svg')
书里改变了图的大小，我认为没什么必要，这样的大小在屏幕正好观察。
# plt.rcParams['figure.figsize'] = (3.5, 2.5)
# plt.scatter(features[:, 1].asnumpy(), labels.asnumpy(), 1);
读取数据集
这里使用到了python的生成器知识点，详见我的博客Python高级特性介绍 - 迭代的99种姿势与协程里有详细介绍。
按小批读取训练数据集，这里每次都判断了min，感觉有点多余但是自己又不知道怎么更好的修改。
# count = []
def data_iter(batch_size, features, labels, ctx):
num_examples = len(features)
indices = list(range(num_examples))
random.shuffle(indices)
for i in range(0, num_examples, batch_size):
j = nd.array(indices[i:min(i + batch_size, num_examples)], ctx=ctx)
# count.append(j)
yield features.take(j), labels.take(j)
# print(count[-1])
这个函数建好之后让我们试试效果：
batch_size = 10
for X, y in data_iter(batch_size, features, labels, ctx):
print(X, y)
break
[[ 0.1555565 -0.89841115]
[ 0.19113457 -0.23471509]
[-0.1649106 -1.0402957 ]
[ 0.59347206 -0.04631801]
[-1.473519 0.7591845 ]
[ 1.1765232 -0.4977442 ]
[ 0.35199913 -0.9372813 ]
[ 0.9608316 0.19639738]
[-1.253151 -0.55073833]
[-0.09366591 -0.84877795]]
<NDArray 10x2 @gpu(3)>
[ 7.580011 5.38007 7.406578 5.5393977 -1.3400362 8.245509
8.093137 5.4496117 3.5762053 6.8978252]
<NDArray 10 @gpu(3)>
初始化模型参数
将模型的权重初始化成均值为0、标准差为0.01的正态随机数，偏差则初始化成0。
w = nd.random.normal(scale=0.01, shape=(num_inputs, 1), ctx=ctx)
b = nd.zeros(shape=(1, ), ctx=ctx)
在模型训练中，需要对这些参数求梯度来迭代参数，所以在这里需要创建他们的梯度。
w.attach_grad()
b.attach_grad()
定义模型
模型就是我们期望拟合的出来的函数，在这即线性回归的矢量计算表达式的实现，下面使用dot做矩阵乘法。
def linreg(X, w, b):
return nd.dot(X, w) + b
定义损失函数
使用前述平方损失来定义线性回归的损失函数，在这里真实值y要reshape成预测值y_hat的形状，实测预测值为二维数组[batch_size, 1]，而真实值为一维数组[batch_size]，这是对书后练习一的解答。
def squared_loss(y_hat, y):
# print(y_hat, y, y.reshape(y_hat.shape))
return (y_hat - y.reshape(y_hat.shape)) ** 2 / 2
定义优化算法
以下sgd函数实现了前述小批量随机梯度下降算法。它通过不断迭代模型参数来优化损失函数。这里自动求梯度模块计算得来的梯度是一个批量样本的梯度和，在笔记「2」中有所介绍，将其除以批量大小即得出梯度的平均值用于迭代。
def sgd(params, lr, batch_size):
for param in params:
param -= lr * param.grad / batch_size
训练模型
训练需要确定（学习速率、批量大小、迭代轮数、模型、损失函数、优化算法），这样就可以开始我们的训练，俗称跑模型。
前述笔记「2」中提到，由于\(l\)并不是一个标量，运行\(l.backward()\)将对\(l\)中元素求和得到新的变量，再求该变量有关模型参数的梯度。
在一个迭代周期(epoch)中，将完整遍历一遍data_iter函数，即对训练数据集中所有样本都使用一次，如果样本数不能被批量大小整除，最后一个批次的样本数目会减少，这是对书后练习三的回答。
这里的迭代周期个数num_epochs和学习速率lr都是超参数，在实践中，大多数超参数都需要通过反复试错来不断调节，俗称调参。
lr = 0.03
num_epochs = 5
net = linreg
loss = squared_loss
for epoch in range(num_epochs):
for X, y in data_iter(batch_size, features, labels, ctx):
with autograd.record():
l = loss(net(X, w, b), y)
l.backward()
sgd([w, b], lr, batch_size)
train_l = loss(net(features, w, b), labels)
print('epoch %d, loss %f' % (epoch, train_l.mean().asnumpy()))
epoch 0, loss 0.027449
epoch 1, loss 0.000099
epoch 2, loss 0.000052
epoch 3, loss 0.000051
epoch 4, loss 0.000051
训练完后就能看到我们训练的数据与真实数据的距离了。
print(true_w, w)
print(true_b, b)
[2, -3.4]
[[ 2.0001543]
[-3.4000373]]
<NDArray 2x1 @gpu(3)>
4.2
[4.1994386]
<NDArray 1 @gpu(3)>
全部代码
这里给出我的notebook的全部代码。
线性回归从零实现.html
简洁实现
上面手写的适合理解原理，但是人和动物的区别就是人会使用工具，可以在前人的基础上进行自己的工作，MXNet就提供了Gluon接口更方便的实现线性回归的训练。
from mxnet.gluon import nn, loss as gloss, data as gdata
from mxnet import gpu, nd, autograd, init, gluon
import matplotlib.pyplot as plt
import random
设一下随机数种子，方便复现。
random.seed(233)
生成数据集
和上面一样，就不多解释了。
ctx = gpu(3)
num_examples = 1000
num_inputs = 2
true_w = [2, -3.4]
true_b = 4.2
features = nd.random.normal(scale=1, shape=(num_examples, num_inputs), ctx=ctx)
labels = true_w[0] * features[:, 0] + true_w[1] * features[:, 1] + true_b
labels += nd.random.normal(scale=0.01, shape=labels.shape, ctx=ctx)
查看一下生成的训练数据集的样子。
features[:5], labels[:5]
(
[[ 0.34113222 0.76659095]
[ 0.6508758 -1.5585263 ]
[ 1.9786316 -0.8465744 ]
[-0.05428567 0.64376295]
[-1.7093595 0.12594718]]
<NDArray 5x2 @gpu(3)>,
[ 2.2640598 10.818717 11.050849 1.8893591 0.35090753]
<NDArray 5 @gpu(3)>)
图像方式查看：
from IPython import display
display.set_matplotlib_formats('svg')
point_size = 1
plt.scatter(features[:, 1].asnumpy(), labels.asnumpy(), point_size);
plt.savefig('散点图.svg')
读取数据集
这里就可以使用Gluon提供的data包来读取数据了。
batch_size = 10
dataset = gdata.ArrayDataset(features, labels)
data_iter = gdata.DataLoader(dataset, batch_size, shuffle=True)
需要注意的是，这里X是来自我们指定的GPU，而y来自了CPU，虽然不明白这是feature还是bug，但是官方论坛上有相应讨论：论坛链接，比较可惜的是，在发文时间2020-03-26，官网论坛炸了，很奇怪为什么亚马逊背书的MXNet的aws服务器会挂掉，从谷歌快照上得到的解决方案还是把数据重新拷贝到GPU上，well。
for X, y in data_iter:
print(X, y)
break
[[ 0.5275261 -1.0045179 ]
[ 0.44837257 -0.43323204]
[-0.01466356 -0.48579317]
[ 1.1003594 0.7395018 ]
[-0.22922567 -1.9637256 ]
[ 0.4390832 -0.60523504]
[-0.75654405 -1.4318156 ]
[-1.1557755 0.4982119 ]
[ 0.5148314 0.9012168 ]
[ 0.42361584 -1.2854003 ]]
<NDArray 10x2 @gpu(3)>
[ 8.658156 6.5668592 5.829563 3.9018939 10.418489 7.1303806
7.569174 0.18902431 2.1697466 9.418344 ]
<NDArray 10 @cpu(0)>
定义模型
Gluon提供的大量预定义层让我们可以只关注使用哪些层来构造模型，是不是和Keras看起来很像（
有一种在搭积木的感觉。
net = nn.Sequential()
线性回归作为单层神经网络，输出层的神经元和输入层中各个输入完全连接。因此，线性回归的输出层又叫全连接层，这里添加一层全连接层，并定义该层的输出个数为1。
net.add(nn.Dense(1))
初始化模型参数
init模块是MXNet提供的，通过其指定权重参数每个元素的初始化值。
net.initialize(init.Normal(sigma=0.01), ctx=ctx)
定义损失函数
Gluon里提供了loss模块，在这里我们直接使用L2范数损失，其实L2就是前面用的平方损失，范数的定义就是这样，比如L3范数就是真实值与预测值差的立方的立方根。
loss = gloss.L2Loss()
定义优化算法
优化算法依然选择sgd，在这里直接使用字符串指定即可。该算法用来迭代net实例所有通过add函数嵌套的层所包含的全部参数。这些参数使用collect_params()函数获取。
trainer = gluon.Trainer(net.collect_params(), 'sgd', {'learning_rate': 0.03})
训练模型
使用Gluon训练时，使用Trainer实例的step函数来迭代训练模型。
num_epochs = 3
for epoch in range(1, num_epochs + 1):
for X, y in data_iter:
with autograd.record():
l = loss(net(X), y.as_in_context(ctx))
l.backward()
trainer.step(batch_size)
l = loss(net(features), labels)
print('epoch %d, loss %f' % (epoch, l.mean().asnumpy()))
epoch 1, loss 0.027621
epoch 2, loss 0.000095
epoch 3, loss 0.000051
从net中获得需要的层后即可访问其权重和偏差，另外也可以获取其梯度，这是对练习三的回答。
dense = net[0]
print(true_w, dense.weight.data())
print(true_b, dense.bias.data())
dense.weight.grad()
[2, -3.4]
[[ 2.000109 -3.3995216]]
<NDArray 1x2 @gpu(3)>
4.2
[4.1992593]
<NDArray 1 @gpu(3)>
Out[14]:
[[0.02693137 0.00585834]]
<NDArray 1x2 @gpu(3)>
练习
练习一说，如果\(l\)替换为mean()，trainer.step(batch_size)相应的需要调整为trainer.step(1)，这是因为迭代求得的梯度会被\(1/batch\_size\)正则化（Gradient will be normalized by 1 / batch_size）。
并且如果需要获取其中的梯度值，也可以自己调用allreduce_grads() + update()实现，step()函数本身就相当于这两个函数的组合，并且需要在autograde.backward()之后、record()作用域之外调用。
net.initialize(init.Normal(sigma=0.01), ctx=ctx, force_reinit=True)
num_epochs = 3
for epoch in range(1, num_epochs + 1):
for X, y in data_iter:
with autograd.record():
l = loss(net(X), y.as_in_context(ctx)).mean()
l.backward()
trainer.step(1)
l = loss(net(features), labels)
print('epoch %d, loss %f' % (epoch, l.mean().asnumpy()))
epoch 1, loss 0.027576
epoch 2, loss 0.000097
epoch 3, loss 0.000051
这里查看的结果和刚才是几乎相同的。
dense = net[0]
print(true_w, dense.weight.data())
print(true_b, dense.bias.data())
[2, -3.4]
[[ 2.0000014 -3.3993788]]
<NDArray 1x2 @gpu(3)>
4.2
[4.1991076]
<NDArray 1 @gpu(3)>
另外关于练习二，官方文档的链接在这gluon.loss提供的损失函数、init提供的初始化方法，笔记用到了再另外介绍。
全部代码
这里给出我的notebook的全部代码。
线性回归简洁实现.html
