深度学习基础之线性回归学习

线性回归

理论

线性回归最简单和最经典的机器学习模型之一。

任何一个机器学习模型都会有如下4个要素：

1. 训练数据
2. 数学模型
3. 损失函数
4. 计算方法

训练数据

下面先阐述一下训练数据的数学定义，假设我们训练数据中有 m 个样本，每个输入样本中有 n 个特征输入和一个标注的输出：

第 i 个样本，记作：

一个模型的表现好坏，训练数据的构建往往是最重要的一步，我们在做模型的时候往往会所这样一句话：

Garbage In, Garbage Out

不过如何构造和优化训练数据往往跟要解决的业务问题领域以及人工标注的质量有关，且跟采用哪种机器学习模型无关，所以就不在本文中详细描述。

而其他三个 数学模型  损失函数 计算方法 都是定义了某一种机器学习模型的关键，所以本文就详细介绍这几个点。

数学模型

线性回归的数学模型定义很简单，就是一维的线性方程：

而训练数据中的第 i 个样本，我们记作：

从上可以看出线性回归的数学假设是，训练数据中所有输出和所有输入的特征之间是一个线性表达式的关系。

您很有可能会问，那如果我输入的训练数据之间他本身不是一个线性关系的话，那这个线性回归不就完全失效了么？
确实是这样的，很多情况下我们拿到的数据并不一定会是线性关系，针对这种情况，我们一般会采取两种办法来解决：

1. 调整数学模型，不用线性方程来表达，用二维，三维甚至更高维度的方程式来表达这个数学关系
2. 调整输入特征，将输入特征做各种函数变化，比如可以将 当成一个特征的话，我们就可以将一个二维模型才能表达的公式转换成一维方程就能表达了

在实际生产环境中，我们一般不会升高数学模型的维度，但是可能会做变化，比如我们之后会讲到的深度学习，就是将多个一维线性模型串联起来作为一个大的模型。

但是不管采用什么样的模型，都是基于假设的，所以机器学习的领域不管怎么优化，都是会有一个准确率和召回率的指标的，因为有可能这个假设就是不对的，当然学习出来的模型也是由偏差的。

说远了，让我们回到我们的线性回归的模型上来，在本文中，为了描述简单，我们假设 n=2 来演示计算过程，比如我们要做一个房价预测系统，特征有两个，一个是房屋面积，一个是房屋年龄，那么我们的计算公式如下：

那么我们的任务就变成了**如何利用训练数据求解 ****[w1, w2, b]** 这几个参数。

在数据量比较少的情况下，我们可以通过初中学过的解方程的思想来求解这个参数，一般来说，只要有3个(即求解参数个数)样本，我们就能求解这个参数了，这个解，我们叫 解析解 。

但是在实际情况中，我们的数据量远远比3个要多，选不同的3个样本的时候求出来的解是不一样的，而最终又无法通过一个 [w1, w2, b] 的参数值来满足所有训练数据集合，即通过解方程得到的结论是 无解 ，那么我们的任务就变成**在所有解当中找到一个最优的近似解，我们称这个解为 数值解 **。

下面两章我们就来介绍，我们如何寻找到最优的近似解。

损失函数

要找到最优解，我们就需要来衡量一个解相比于另一个解到底哪个解更好。怎么衡量呢，就是我们本节要讲的 损失函数 了。

以如上的房屋面积的例子来距离，假设我么有两组解：

然后我们用训练数据，分别用 w 和 w' 来计算出一个房屋价格的数组，分别为：

而输入的训练样本的真实的房屋面积的数组为：

那么我们就是比较  和 哪个离我们的真实数据   更相近，哪个就更好，那么怎么定义呢，一般我们采用每个元素之差的平方和来定义，如下：

这个函数，我们也就称之为损失函数。如上两组参数哪个好，只需要比较  和 对应的损失函数取得的值，哪个损失函数的值比较小，我们就选择哪一组参数。

当然，损失函数也不一定要采用平方，用绝对值，用四次方都可以，用什么函数，就跟下一步求解过程有关了，针对线性回归问题，我们一般都采用的是平方和损失函数。

那么到目前为止，我们的任务已经变成了，在所有可能的 [w1, w2, b] 的组合中，找寻到一组解，使得 loss(y) 损失函数的取值最小。

计算过程

如上一节损失函数中描述，我们已经将我们的问题求解过程转换成了一个数学求解的任务，找到损失函数最小的参数解。

这个求解过程，最广泛被使用的方法就是随机梯度下降算法。该算法的指导思想也很简单，就是随机选择一组参数值 [w]  作为初始值，然后执行N次循环，每次循环都根据优化算法改变一点点 w  的值，使得损失函数能下降一点点，经过有限次循环之后，如果发现损失函数基本上不怎么变化了，就可以停止循环，将当前计算出来的 w 当成是近似最优解返回。

那么您可能也看出来了，如上的描述中提到了随机，也提到了下降，但是名字中的梯度没提到，同时可能您也看出来了，上面描述中还有一个黑盒一般的描述：**每次根据优化算法改变一点点 ****w** 的值。那其实这里的优化算法就是对损失函数求梯度，梯度是微积分上的一个数学定义，定义参考我们在数学基础一章中的描述。这里简单说一下，就是针对每个 w 中的参数，我们针对损失函数求该参数的偏导数，然后将该参数沿着偏导数的方向往下走一小段，这样大概率下一次我们的损失函数是会下降的。

以房价为例子，我们的优化算法的数学表达式如下：

那么我们根据如上定义的平方损失函数，针对w1计算偏导数如下：

注意我们在损失函数前面增加了一个 1/2 ，主要为了求导更方便。

然后我们将偏导数代入上面的公式中，可以得到如下的公式：

如上公式中的  的定义为学习率，即每次都将该参数的偏导数乘以一个系数再下降，如果学习率选得大，那么迭代次数会更少，学习得会更快，但是有可能会导致跳过最优解，甚至无法收敛，而最终有可能训练时间比更小的学习率的时间更长；如果学习率选得小，那么会增加循环迭代次数，增加模型训练的时间。而如何选择这个学习率，并不是自动的，而是靠经验主义，在训练模型之初人工设定的，像这种靠人工一开始输入而不是学习出来的参数，就叫超参数。

在如上的更新参数的过程中，我们每次都全量计算了所有的训练数据集合，在实际应用场景中，为了提高模型训练效率，我们一般不会这样做，而是每次迭代只选择一小部分数据进行计算，那么怎么选取这一小部分呢，我们一般会将训练数据随机打乱，按照窗口大小从前往后遍历，每次都取一小部分数据进行计算，比如假设每次取100条数据，那么第一次我们会取1_{100条，第二次回取101}200条，以此类推，直到将训练数据都遍历完。当遍历完了之后很有可能模型还是没有训练好的，那么就再次重新打乱所有训练数据，按照窗口从前往后取数据再进行第二轮训练。所以综上，我们的求最佳参数的过程都是 **按轮 按次 **重复进行训练的。

在增加了按轮按次进行了之后，我们的迭代计算公式就需要调整了，如下：

如上公式中的 B ，就是对应轮次的选取的样本集合，并且我们将学习率除以了窗口大小，来保证学习率的经验数字不需要随着窗口选择大小而改变。

最后我们来说一下如何判断模型训练是否结束，一般有两种方法：

1. 超参数指定训练的轮数和次数，到了次数就结束
2. 每轮次根据计算出来的损失函数的变化率来判断是否应该结束

向量表示

我们计算机在执行运算的时候，对矩阵预算是比较友好的，尤其是GPU机器，所以我们在实现模型训练的代码的时候，尽量要用向量(矩阵)计算的方式，而不是用 for 循环逐个计算。

from mxnet import nd
from time import time

# 初始化几个数组
a = nd.ones(1000)
b = nd.ones(1000)
c = nd.zeros(1000)

# 普通操作，通过循环赋值的方式操作
start1 = time()
for i in range(1000):
    c[i] = a[i] + b[i]
end1 = time()

print("normal time: ", end1-start1)

# 向量操作，利用API一次操作
start2 = time()
d = a + b
end2 = time()

print("vector oper time: ", end2-start2)

该段代码的输出为：

normal time:  0.16321206092834473
vector oper time:  0.0009930133819580078

可以看出，普通操作是通过向量操作方式的 200倍！

所以，我们尽可能的将所有数学操作都通过向量表达式来操作，来提高模型训练效率。

在线性回归的计算式当中，主要是有两个需要计算的数学公式，一个是计算损失函数，一个是计算损失函数的梯度。

首先用向量来定义训练样本，定义如下：

定义要学习的参数向量为：

那么给定当前参数，计算出来的输出数组的向量计算表达式看起来就很简单了，如下：

然后我们再来看如上线性回归定义的损失函数，原始定义如下：

调整为向量表达式之后变成如下了：

然后参数的递归梯度下降的公式也可以用向量公式来替换，先回忆一下原始的定义：

然后再来看一下新的向量计算方式：

这样一来，可以看出，原本来起来还比较复杂的公式一下子看起来就很简单清晰了。

下面我们就摩拳擦掌，开始来按照向量表示的方式自己从零实现一个线性回归的模型训练吧。

普通实现

构造训练数据

为了理解简单，我们在实现过程中，演示一个只有两个参数的例子，下面我们首先来构造一下训练数据，具体请参考注释：

from mxnet import autograd, nd
import random

true_w = nd.array([[2,-3.4]])
true_b = 4.2
num_inputs = 2
num_examples = 1000

# 生成随机的特征数组
features = nd.random.normal(scale=1, shape=((num_examples, num_inputs)))
# 根据特征数组生成label数组
labels = nd.dot(features, true_w.T) + true_b
# 再给label数组加上一个随机的噪音
labels += nd.random.normal(scale=0.01, shape=labels.shape)

定义遍历函数和超参数

# 我们定义一个按轮按次轮询数据的函数                                                                          
def data_iter(batch_size, X, Y):                                                             
    # 定义一个下标数组并打印出来，并将下标做随机打乱                                                                
    num_inputs = len(X)                                                                      
    indices = list(range(num_inputs))                                                        
    random.shuffle(indices)                                                                  
    # 每次都从下标数组中取一个 batch_size 窗口大小的数据出来                                                      
    for i in range(0, num_inputs, batch_size):                                               
        cur_indices = nd.array(indices[i:min(i+batch_size, num_inputs)])                     
        yield X.take(cur_indices), Y.take(cur_indices)                                       
                                                                                             
# 定义超参数轮和每轮需要遍历的次                                                                            
# 对训练数据总共遍历多少轮                                                                               
loop_num = 3                                                                                 
# 每轮的小批量的窗口大小                                                                                
batch_size = 10                                                                              
# 梯度下降的学习率                                                                                   
learn_rate = 0.5                                                                             
                                                                                             
# 定义我们要学习的参数并做随机初始化, shape (3, 1)                                                            
W = nd.random.normal(scale=1, shape=(num_features+1, 1))                                     

开始训练

start_time = time()
for epoch_num in range(loop_num):
    for batchX, batchY in data_iter(batch_size, X, Y):
      	# 按照上面向量表示的梯度下降参数执行
        W[:] = W - (learn_rate / batch_size) * nd.dot(batchX.T, nd.dot(batchX, W) - batchY)
end_time = time()

# 打印最初的初始值和训练出来的值
print("train time: ", end_time-start_time)
print(true_W, W)

到此为止，我们的线性回归实现就完成了，该脚本的实现输出为：

train time:  0.09574174880981445

[[2.1]
 [3.7]
 [1.2]]
<NDArray 3x1 @cpu(0)>
[[2.1006684]
 [3.6992905]
 [1.1961068]]
<NDArray 3x1 @cpu(0)>

可以看到针对一个1000条数据的训练过程还是很快的，只需要不到 0.1 秒就训练完了，同时训练出来的结果跟一开始构造数据的参数是几乎相同的。

下面，我们引入mxnet的自动求梯度的机制，利用这个机制，我们甚至可以不用了解损失函数求导的过程，只需要通过 mxnet 提供的API接口计算出损失函数，然后 mxnet 就能自动推导出要计算的参数的梯度。

虽然用在线性回归上是有点大材小用了，但是对一些比较复杂的数学模型（比如深度学习），这个机制就会显得非常有用了，大大提高了我们代码的编写效率。
 

利用mxnet的自动求梯度

# 我们按照mxnet自动求梯度的接口先调用attach_grad()函数对要求梯度的矩阵/向量进行标记
W.attach_grad()
start_time = time()
for epoch_num in range(loop_num):
    for batchX, batchY in data_iter(batch_size, X, Y):
        # 在这个标记片段内所有计算操作都会被记录下来并针对标记的参数求梯度
        with autograd.record():
            # 计算损失函数
            batchY_hat = nd.dot(batchX, W)
            loss = (batchY_hat - batchY)**2/2
        # 针对损失函数求梯度, 经过这个函数之后，后面就可以一共W.grad来获取此次计算针对W的梯度了
        loss.backward()
        # 然后将参数按照梯度进行下降
        W[:] = W - (learn_rate / batch_size) * W.grad

end_time = time()

# 打印最初的初始值和训练出来的值
print("train time: ", end_time-start_time)
print(true_W, W)

这段代码的输出如下：

train time:  0.19647645950317383

[[2.1]
 [3.7]
 [1.2]]
<NDArray 3x1 @cpu(0)>
[[2.102761 ]
 [3.701326 ]
 [1.1980474]]
<NDArray 3x1 @cpu(0)>

看起来代码要更复杂一些了，训练效率也要低一些了，下降了差不多一倍，训练参数的结果是差不多的，主要的好处是对将来的更复杂的模型打下了代码实现的基础。

加上中间输出

我们在如下的代码中，加上了一些中间输出结果，来方便我们下一步来调试不同的超参数对模型训练过程的影响：

#  求线性回归的损失函数
def get_loss(X, Y, W):
    YHat = nd.dot(X, W)
    return (YHat - Y)**2/2
  
# 我们按照mxnet自动求梯度的接口先调用attach_grad()函数对要求梯度的矩阵/向量进行标记
W.attach_grad()
start_time = time()
for epoch_num in range(loop_num):
    cur_batch_num = 0
    for batchX, batchY in data_iter(batch_size, X, Y):
        # 在这个标记片段内所有计算操作都会被记录下来并针对标记的参数求梯度
        with autograd.record():
            # 计算损失函数
            loss = get_loss(batchX, batchY, W)
        # 针对损失函数求梯度, 经过这个函数之后，后面就可以一共W.grad来获取此次计算针对W的梯度了
        loss.backward()
        # 然后将参数按照梯度进行下降
        W[:] = W - (learn_rate / batch_size) * W.grad

        cur_batch_num+=1
        loss = get_loss(batchX, batchY, W)
        if(cur_batch_num%50 == 0):
            print("epoch_num: %d, batch_num: %d, loss: %f" %
                    (epoch_num, cur_batch_num, loss.mean().asnumpy()))
    # 打印一些训练过程中的信息
    cur_loss = get_loss(X, Y, W)
    print("## epoch_num: %d, loss: %f
" % (epoch_num, cur_loss.mean().asnumpy()))
end_time = time()

# 打印最初的初始值和训练出来的值
print("train time: ", end_time-start_time)
print(true_W, W)

我们再来看看这个代码的输出：

epoch_num: 0, batch_num: 50, loss: 0.000046
epoch_num: 0, batch_num: 100, loss: 0.000050
## epoch_num: 0, loss: 0.000052

epoch_num: 1, batch_num: 50, loss: 0.000061
epoch_num: 1, batch_num: 100, loss: 0.000029
## epoch_num: 1, loss: 0.000062

epoch_num: 2, batch_num: 50, loss: 0.000052
epoch_num: 2, batch_num: 100, loss: 0.000022
## epoch_num: 2, loss: 0.000053

train time:  0.2802696228027344

[[2.1]
 [3.7]
 [1.2]]
<NDArray 3x1 @cpu(0)>
[[2.1014354]
 [3.6990693]
 [1.1987416]]
<NDArray 3x1 @cpu(0)>

不同超参数对训练过程的影响

学习率
固定 轮次：3， 小批量窗口 10
我们来看不同的学习率对训练过程的影响，如下表：

超参数	第一轮损失值	第二轮损失值	第三轮损失值
0.01	1.89	0.26	0.03
0.02	0.156179	0.003019	0.000111
0.03	0.027368	0.000096	0.000049
0.05	0.000423	0.000045	0.000045
0.1	0.000051	0.000050	0.000051
0.5	0.000054	0.000065	0.000055
1	0.000073	0.000081	0.000107
1.5	0.000081	0.000101	0.000296
1.8	0.001779	0.001116	0.001254
1.9	22.656252	8493858816	4272635665084055552
2.0	35379048448	nan	nan

从如上表我们可以得出如下的结论：

1. 在学习率数值较小的一段范围内，学习率越大，模型收敛得就越快，理论上最终模型的模拟效果会越好
2. 在学习率数值适中的一段范围内，最终学习效果都差不多
3. 当学习率数值较大的范围内，效果反而更不可控了，甚至到最后损失函数的值都超过了整数范围值

所以学习率千万不能设置得太大，设置太大之后不收敛了；设置太小模型还是能收敛只是训练时间要更长，最好是要选择得适中，比如在我们的实现例子中， 0.05 就是一个不错的case。

小批量窗口
固定学习率0.05，选择不同的小批量窗口对应的结果如下：

批量窗口	第一轮损失值	第二轮损失值	第三轮损失值
10	0.000288	0.000047	0.000047
20	0.076830	0.000432	0.000048
50	0.759479	0.098075	0.012729

可以看到，在合适范围内，窗口值越小就收敛得越快。

通用实现

虽然我们在如上的代码中实现线性回归的逻辑已经很简单了，但是我们还是自己写了不少的一些函数，比如我们自己实现了数据遍历的函数，自己实现了损失计算函数，自己实现了梯度下降等逻辑。

而实际上一些通用的深度学习基础库，比如 mxnet 已经为我们预先实现了很多常用的函数库了，我们可以直接调用他们的服务来重新实现一个线性回归训练脚本。

并且在这个实现中，我们会按照典型的训练模型的代码框架来实现线性回归的训练，按照这个代码框架，即使后面我们实现更复杂的深度学习的多层神经网络模型的时候，代码看起来都是类似的。

from mxnet import autograd, nd
import random

# 构造数据
true_w = nd.array([[2,-3.4]])
true_b = 4.2
num_inputs = 2
num_examples = 1000

# 生成随机的特征数组
features = nd.random.normal(scale=1, shape=((num_examples, num_inputs)))
# 根据特征数组生成label数组
labels = nd.dot(features, true_w.T) + true_b
# 再给label数组加上一个随机的噪音
labels += nd.random.normal(scale=0.01, shape=labels.shape)

# 引入gluon的管理训练集的库
from mxnet.gluon import data as gdata

batch_size = 10
# 将样本的特征值和标签值都传入
dataset = gdata.ArrayDataset(features, labels)
# 利用gdata提供的数据抽取小批量遍历函数指针
data_iter = gdata.DataLoader(dataset, batch_size, shuffle=True)

# 接下来将是最主要的区别，在手动实现中，我们自己实现了线性函数和损失函数
# 而gluon当中提供了一系列默认的函数，可以快速的模型训练
# 第一步，搭建一个神经网络, 虽然我们现在的神经网络还很简单，只有一层
from mxnet.gluon import nn
net = nn.Sequential()
net.add(nn.Dense(1))

# 第二步，将参数矩阵给初始化好
from mxnet import init
net.initialize(init.Normal(sigma=0.01))

# 定义损失函数
from mxnet.gluon import loss as gloss
loss = gloss.L2Loss()

# 定义参数梯度下降的函数
from mxnet import gluon
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()))

# 对比训练出来的值的效果
dense = net[0]
print(true_w, dense.weight.data())
print(true_b, dense.bias.data())