深度学习基础——线性回归（Python实现）

---

前言

本文主要对线性回归的原理以及基于Python 的代码实现进行简单的讲解。

ps: 文中主要部分来自于《动手学深度学习》(第一版) 一 书，此文为笔者学习笔记，仅作为参考，如有错误，请联系笔者更正，欢迎讨论学习。
本文参考自：《动手学深度学习》线性回归章节

环境配置：
name: gluon
dependencies:

• python=3.6
• pip:
  • mxnet==1.5.0
  • d2lzh==1.0.0
  • jupyter==1.0.0
  • matplotlib==2.2.2
  • pandas==0.23.4

其中 d2lzh 为《动手学深度学习》书籍附带的代码包。

---

一、线性回归原理

线性回归是单层神经网络，其输出值为连续值，在日常生活中比较常见，如房屋价格，气温，销售额的预测等，我们以房屋价格预测为例简单解释线性回归的原理。

1. 线性回归模型

设房屋的面积为 $x_1$，房龄为 $x_2$，售出价格为 $y$。$y$ 与 $x_1$，$x_2$ 线性相关，由此可以得到线性回归模型：
$$\hat{y}=x_1{w}_1+x_2{w}_2+b$$其中$w_1$，$w_2$ 为权重，$b$ 为偏差，$\hat{y}$ 为预测房价，一般情况下 $\hat{y}$ 与真实值 $y$ 之间有一定的误差。

当样本数为 $n$，特征数为 $d$ 时可定义广义线性回归模型：
$$\hat{\mathbf{y}}=\mathbf{X}\mathbf{w}+b$$其中模型输出 $\hat{\mathbf{y}}$$\in$${\mathbb{R}}^{n\times 1}$，样本数据特征 $\mathbf{X}$$\in$${\mathbb{R}}^{n\times d}$，权重 $\mathbf{w}$$\in$${\mathbb{R}}^{d\times 1}$，偏差 b$\in$$\mathbb{R}$。

2. 损失函数

在模型训练中，我们通常需要选取一个损失函数来监测模型的质量，以便我们训练合适的模型，本次选用平方函数作为损失函数，单个样本的损失函数 ${\ell }^{(i)}$ 形式如下：
$${\ell }^{(i)}(w_1,w_2,b)=\frac{1}{2}({\hat{y}}^{(i)}-y^{(i)}{)}^2$$其中 ${\ell }^{(i)}$ 为第 $i$ 个样本的误差，${\hat{y}}^{(i)}$ 为第 $i$ 个样本的样本预测值，$y^{(i)}$ 为第 $i$ 个样本的样本标签值,常数 $\frac{1}{2}$ 只是为了求导之后与平方的2约分为1，只是为了简便运算。
以所有样本的误差值的平均值 $\ell$ 衡量模型质量，则有：
$$\ell (w_1,w_2,b)=\frac{1}{n}\sum _{i=1}^n{\ell }^{(i)}(w_1,w_2,b)=\frac{1}{n}\sum _{i=1}^n\frac{1}{2}({\hat{y}}^{(i)}-y^{(i)}{)}^2=\frac{1}{n}\sum _{i=1}^n\frac{1}{2}(x_1^{(i)}{w}_1+x_2^{(i)}{w}_2+b-y^{(i)}{)}^2$$定义 $\theta =[w_1,w_2,\dots w_d,b{]}^T$ 为模型参数，则广义损失函数矢量形式为：
$$\ell (\theta )=\frac{1}{2n}({\hat{y}}^{(i)}-y^{(i)}{)}^T({\hat{y}}^{(i)}-y^{(i)})$$我们通过模型训练就是要找到模型平均误差最低时所对应一组参数 $${\theta }^{\ast }=[w_1^{\ast },w_2^{\ast },\dots w_d^{\ast },b^{\ast }{]}^T=\underset{\theta }{\mathrm{argmin}}\ell (\theta )$$

3. 优化算法

采用小批量随机梯度下降法，先选取一组模型参数的初始值，如随机选取；接下来对参数进行多次迭代，使每次迭代都可能降低损失函数的值。在每次迭代中，先随机均匀采样一个由固定数目训练数据样本所组成的小批量 $\mathcal{B}$ ，然后求小批量中数据样本的平均损失有关模型参数的导数（梯度），最后用此结果与预先设定的一个正数的乘积作为模型参数在本次迭代的减小量。
在训练本次讨论的线性回归模型的过程中，模型的每个参数将作如下迭代：先迭代 $w$ ，后迭代 $b$ 。如果先迭代 $b$ ,则 $w$ 迭代项最后一项的值会变，所以顺序不能变。
$$w_1\leftarrow w_1-\frac{\eta }{|\mathcal{B}|}\sum _{i\in \mathcal{B}}\frac{\partial {\ell }^{(i)}(w_1,w_2,b)}{\partial w_1}=w_1-\frac{\eta }{|\mathcal{B}|}\sum _{i\in \mathcal{B}}{x}_1^{(i)}(x_1^{(i)}{w}_1+x_2^{(i)}{w}_2+b-y^{(i)}),$$
$$w_2\leftarrow w_2-\frac{\eta }{|\mathcal{B}|}\sum _{i\in \mathcal{B}}\frac{\partial {\ell }^{(i)}(w_{,}{w}_2,b)}{\partial w_2}=w_2-\frac{\eta }{|\mathcal{B}|}\sum _{i\in \mathcal{B}}{x}_2^{(i)}(x_1^{(i)}{w}_1+x_2^{(i)}{w}_2+b-y^{(i)}),$$
$$b\leftarrow b-\frac{\eta }{|\mathcal{B}|}\sum _{i\in \mathcal{B}}\frac{\partial {\ell }^{(i)}(w_{,}{w}_2,b)}{\partial b}=b-\frac{\eta }{|\mathcal{B}|}\sum _{i\in \mathcal{B}}(x_1^{(i)}{w}_1+x_2^{(i)}{w}_2+b-y^{(i)})。$$
在上式中，$|\mathcal{B}|$代表每个小批量中的样本个数（批量大小，batch size），$\eta$称作学习率（learning rate）并取正数。需要强调的是，这里的批量大小和学习率的值是人为设定的，并不是通过模型训练学出的，因此叫作超参数（hyperparameter）。我们通常所说的“调参”指的正是调节超参数，例如通过反复试错来找到超参数合适的值。

梯度下降的矢量形式迭代步骤写法为：
$$\theta =\theta -\frac{\eta }{|\mathcal{B}|}\sum _{i\in \mathcal{B}}{\nabla }_{\theta }{\ell }^{(i)}(\theta )$$
其中${\nabla }_{\theta }{\ell }^{(i)}(\theta )$是损失函数关于模型参数的偏导数所组成的向量。

二、线性回归的从零实现

1. 引入库

代码如下：

%matplotlib inline
from IPython import display
from matplot import pyplot as plt
from mxnet import autograd,nd
import random

2. 生成数据集

我们构造一个简单的人工训练数据集，它可以使我们能够直观比较学到的参数和真实的模型参数的区别。设训练数据集样本数为1000，输入个数（特征数）为2。给定随机生成的批量样本特征 $\mathbf{X}$$\in$${\mathbb{R}}^{1000\times 2}$，我们使用线性回归模型真实权重$w=[2,-3.4{]}^T$和偏差$b$ = 4.2，以及一个随机噪声项 $ϵ$ 来生成标签
$$\hat{\mathbf{y}}=\mathbf{X}\mathbf{w}+b$$
其中噪声项 $ϵ$ 服从均值为0、标准差为0.01的正态分布。噪声代表了数据集中无意义的干扰。
代码如下：

num_inputs = 2
num_examples = 1000
true_w = [2, -3.4]  #真实权重
true_b = 4.2  #真实偏差
features = nd.random.normal(scale=1, shape=(num_examples, num_inputs))
labels = true_w[0] * features[:, 0] + true_w[1] * features[:, 1] + true_b
labels += nd.random.normal(scale=0.01, shape=labels.shape)  #加入随机干扰

features的每一行是一个长度为2的向量，labels的每一行是一个长度为1 的标量。由于数据是随机生成，所以只要在范围内即可。

通过生成第二个特征features[:, 1]和标签 labels 的散点图，可以更直观地观察两者间的线性关系。由于特征值是两个，图示将第一个特征值映射到该平面以展示。

def use_svg_display():
    # 用矢量图显示
    display.set_matplotlib_formats('svg')

def set_figsize(figsize=(3.5, 2.5)):
    use_svg_display()
    # 设置图的尺寸
    plt.rcParams['figure.figsize'] = figsize

set_figsize()
plt.scatter(features[:, 1].asnumpy(), labels.asnumpy(), 1);  # 加分号只显示图

3. 读取数据集

在训练模型的时候，我们需要遍历数据集并不断读取小批量数据样本。这里我们定义一个函数：它每次返回batch_size（批量大小）个随机样本的特征和标签。

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):
        j = nd.array(indices[i: min(i + batch_size, num_examples)])
        yield features.take(j), labels.take(j)  # take函数根据索引返回对应元素 #yield 函数返回迭代对象 与return有一定差别

让我们读取第一个小批量数据样本并打印。每个批量的特征形状为(10, 2)，分别对应批量大小和输入个数；标签形状为批量大小。

batch_size = 10

for X, y in data_iter(batch_size, features, labels):
    print(X, y)
    break #该break结束的是上个代码块中的for循环

类似的可以通过计数器得到第i个批次的样本数据，例如前两个批次的数据：

4. 初始化模型参数

初始化一般任意都可以，我们通常选取小一点的数，将权重初始化成均值为0、标准差为0.01的正态随机数，偏差则初始化成0。

w = nd.random.normal(scale=0.01, shape=(num_inputs, 1))
b = nd.zeros(shape=(1,))

我们需要对参数求梯度来更迭参数的值，所以我们需要对创建参数的梯度。attach_grad()函数用于申请内存。

w.attach_grad()
b.attach_grad()

5. 定义模型

定义线性回归模型
$$\hat{\mathbf{y}}=\mathbf{X}\mathbf{w}$$

def linreg(X, w, b):  
    return nd.dot(X, w) + b

6. 定义损失函数

定义平方损失函数
$$\ell (w_1,w_2,b)=\frac{1}{n}\sum _{i=1}^n{\ell }^{(i)}(w_1,w_2,b)=\frac{1}{n}\sum _{i=1}^n\frac{1}{2}({\hat{y}}^{(i)}-y^{(i)}{)}^2$$
在实现中，我们需要把真实值y变形成预测值y_hat的形状。以下函数返回的结果也将和y_hat的形状相同。故调用shape函数。

def squared_loss(y_hat, y):  # 本函数已保存在d2lzh包中方便以后使用
    return (y_hat - y.reshape(y_hat.shape)) ** 2 / 2

7. 定义优化算法

这里自动求梯度模块计算得来的梯度是一个批量样本的梯度和。我们将它除以批量大小来得到平均值。

def sgd(params, lr, batch_size): 
    for param in params:
        param[:] = param - lr * param.grad / batch_size

8. 训练模型

在训练中，我们将多次迭代模型参数。在每次迭代中，我们根据当前读取的小批量数据样本（特征 X 和标签 y ），通过调用反向函数backward计算小批量随机梯度，并调用优化算法 sgd 迭代模型参数。由于我们之前设批量大小 batch_size 为10，每个小批量的损失l的形状为(10, 1)。由于变量 l 并不是一个标量，运行 l.backward() 将对 l 中元素求和得到新的变量，再求该变量有关模型参数的梯度。

在一个迭代周期（epoch）中，我们将完整遍历一遍data_iter函数，并对训练数据集中所有样本都使用一次（假设样本数能够被批量大小整除）。这里的迭代周期个数num_epochs和学习率lr都是超参数，分别设3和0.03。在实践中，大多超参数都需要通过反复试错来不断调节。虽然迭代周期数设得越大模型可能越有效，但是训练时间可能过长。

lr = 0.03
num_epochs = 3
net = linreg
loss = squared_loss

for epoch in range(num_epochs):  # 训练模型一共需要num_epochs个迭代周期
    # 在每一个迭代周期中，会使用训练数据集中所有样本一次（假设样本数能够被批量大小整除）。X
    # 和y分别是小批量样本的特征和标签
    for X, y in data_iter(batch_size, features, labels):
        with autograd.record():
            l = loss(net(X, w, b), y)  # l是有关小批量X和y的损失
        l.backward()  # 小批量的损失对模型参数求梯度
        sgd([w, b], lr, batch_size)  # 使用小批量随机梯度下降迭代模型参数
        #这里可以记录每批次迭代的结果
    #这里可以记录每周期迭代的结果
    train_l = loss(net(features, w, b), labels)
    print('epoch %d, loss %f' % (epoch + 1, train_l.mean().asnumpy()))

训练完成后，我们可以比较学到的参数和用来生成训练集的真实参数。它们应该很接近。

完整代码如下：

#完整代码
%matplotlib inline
from IPython import display
from matplotlib import pyplot as plt
from mxnet import autograd,nd
import random
num_inputs=2
num_examples=1000
true_w=nd.array([2, -3.4])
true_b=4.2
features = nd.random.normal(scale=1,shape=(num_examples, num_inputs))
labels = nd.dot(features, true_w.T) + true_b
labels += nd.random.normal(scale=0.01 , shape = labels.shape)
def use_svg_display():
    # 用矢量图显示
    display.set_matplotlib_formats('svg')

def set_figsize(figsize=(3.5, 2.5)):
    use_svg_display()
    # 设置图的尺寸
    plt.rcParams['figure.figsize'] = figsize

set_figsize()
plt.scatter(features[:, 1].asnumpy(), labels.asnumpy(), 1);  # 加分号只显示图
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):
        j = nd.array(indices[i: min(i + batch_size, num_examples)])
        yield features.take(j), labels.take(j)  # take函数根据索引返回对应元素
w = nd.random.normal(scale=0.01, shape=(num_inputs, 1))
b = nd.zeros(shape=(1,))
w.attach_grad()
b.attach_grad()
def linreg(X, w, b): 
    return nd.dot(X, w) + b
def squared_loss(y_hat, y):  
    return (y_hat - y.reshape(y_hat.shape)) ** 2 / 2
def sgd(params, lr, batch_size):  
    for param in params:
        param[:] = param - lr * param.grad / batch_size
lr = 0.03
num_epochs = 3
net = linreg
loss = squared_loss

for epoch in range(num_epochs):  # 训练模型一共需要num_epochs个迭代周期
    # 在每一个迭代周期中，会使用训练数据集中所有样本一次（假设样本数能够被批量大小整除）。X
    # X和y分别是小批量样本的特征和标签
    for X, y in data_iter(batch_size, features, labels):
        with autograd.record():
            l = loss(net(X, w, b), y)  # l是有关小批量X和y的损失
        l.backward()  # 小批量的损失对模型参数求梯度
        sgd([w, b], lr, batch_size)  # 使用小批量随机梯度下降迭代模型参数
    train_l = loss(net(features, w, b), labels)
    print('epoch %d, loss %f' % (epoch + 1, train_l.mean().asnumpy()))

三. 线性回归的简洁实现

该节借用MXNet提供的Gluon接口更快捷的实现线性回归的训练。环境配置请参考前言部分超链接。关于mxnet的详细函数说明请参见：mxnet文档 。

1. 生成数据集

本节依然与上节生成数据集的方法相同。

from mxnet import autograd, nd

num_inputs = 2
num_examples = 1000
true_w = [2, -3.4]
true_b = 4.2
features = nd.random.normal(scale=1, shape=(num_examples, num_inputs))
labels = true_w[0] * features[:, 0] + true_w[1] * features[:, 1] + true_b
labels += nd.random.normal(scale=0.01, shape=labels.shape)

2. 读取数据集

Gluon提供了 data 包来读取数据。由于 data 常用作变量名，我们将导入的data模块用添加了 Gluon 首字母的假名 gdata 代替。在每一次迭代中，我们将随机读取包含10个数据样本的小批量。

from mxnet.gluon import data as gdata

batch_size = 10
# 将训练数据的特征和标签组合
dataset = gdata.ArrayDataset(features, labels)
# 随机读取小批量
data_iter = gdata.DataLoader(dataset, batch_size, shuffle=True)

data_iter 的使用与上一节中的一样。读取并打印第一个小批量数据样本。

3. 定义模型

首先，导入nn模块。实际上，“nn” 是neural networks（神经网络）的缩写。顾名思义，该模块定义了大量神经网络的层。我们先定义一个模型变量net，它是一个Sequential实例。在Gluon中，Sequential实例可以看作是一个串联各个层的容器。在构造模型时，我们在该容器中依次添加层。当给定输入数据时，容器中的每一层将依次计算并将输出作为下一层的输入。

from mxnet.gluon import nn

net = nn.Sequential()

作为一个单层神经网络，线性回归输出层中的神经元和输入层中各个输入完全连接。因此，线性回归的输出层又叫全连接层。在Gluon中，全连接层是一个Dense实例。我们定义该层输出个数为1。

net.add(nn.Dense(1))

在Gluon中我们无须指定每一层输入的形状，例如线性回归的输入个数。当模型得到数据时，例如后面执行net(X)时，模型将自动推断出每一层的输入个数。

4. 初始化模型参数

在Gluon中，loss模块定义了各种损失函数。我们用假名gloss代替导入的loss模块，并直接使用它提供的平方损失作为模型的损失函数。

from mxnet.gluon import loss as gloss

loss = gloss.L2Loss()  # 平方损失又称L2范数损失

5. 定义优化算法

导入 Gluon 后，我们创建一个Trainer实例，并指定学习率为 0.03 的小批量随机梯度下降（sgd）为优化算法。该优化算法将用来迭代 net 实例所有通过add函数嵌套的层所包含的全部参数。这些参数可以通过collect_params 函数获取。

from mxnet import gluon

trainer = gluon.Trainer(net.collect_params(), 'sgd', {'learning_rate': 0.03})

6. 训练模型

在使用 Gluon 训练模型时，我们通过调用 Trainer 实例的 step 函数来迭代模型参数。上一节中我们提到，由于变量 l 是长度为 batch_size 的一维 NDArray，执行 l.backward() 等价于执行 l.sum().backward() 。按照小批量随机梯度下降的定义，我们在 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)
        l.backward()
        trainer.step(batch_size)#参数迭代，梯度会以1/batch.size标准化，如果loss函数已经标准化则设置为1.
    l = loss(net(features), labels)
    print('epoch %d, loss: %f' % (epoch, l.mean().asnumpy()))

我们分别比较学习到的模型参数和真实的模型参数。我们从net获得需要的层，并访问其权重（weight）和偏差（bias）。学到的参数和真实的参数很接近。

完整代码如下：

#完整代码
from mxnet import autograd, nd
from mxnet.gluon import data as gdata
from mxnet.gluon import nn
from mxnet import init
from mxnet.gluon import loss as gloss
from mxnet import gluon

#生成数据集
num_inputs = 2
num_examples = 1000
true_w = [2, -3.4]
true_b = 4.2
features = nd.random.normal(scale=1, shape=(num_examples, num_inputs))
labels = true_w[0] * features[:, 0] + true_w[1] * features[:, 1] + true_b
labels += nd.random.normal(scale=0.01, shape=labels.shape)

#读取数据集
batch_size = 10
# 将训练数据的特征和标签组合
dataset = gdata.ArrayDataset(features, labels)
# 随机读取小批量
data_iter = gdata.DataLoader(dataset, batch_size, shuffle=True)

#建立模型
net = nn.Sequential()
net.add(nn.Dense(1))

#初始化模型参数
net.initialize(init.Normal(sigma=0.01))

#定义损失函数
loss = gloss.L2Loss()  # 平方损失又称L2范数损失

#定义优化算法
trainer = gluon.Trainer(net.collect_params(), 'sgd', {'learning_rate': 0.03})

#模型训练
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)
        l.backward()
        trainer.step(batch_size)
    l = loss(net(features), labels)
    print('epoch %d, loss: %f' % (epoch, l.mean().asnumpy()))

---

总结

本文分别介绍了两种线性回归的实现方法，显然后者相比于前者实现训练较为简单，重点关注一下求梯度的步骤实现。