神经网络基本概念

• 神经元、多层连接、前向计算、计算图

模型结构三要素

• 模型假设、评价函数和优化算法

一、问题分析

对于预测问题，可以根据预测输出的类型是连续的实数值，还是离散的标签，区分为回归任务和分类任务。因为房价是一个连续值，所以房价预测显然是一个回归任务。下面我们尝试用最简单的多元线性回归模型解决这个问题，并用神经网络来实现这个模型。

1.1 线性回归模型

假设房价和各影响因素之间能够用线性关系来描述：
$$y=\sum _{j=1}^M{x}_j{w}_j+b(公式1)$$
模型的求解即是通过数据拟合出每个$w_j$和$b$。其中，$w_j$和$b$分别表示该线性模型的权重和偏置。一维情况下，$w_j$ 和 $b$ 是直线的斜率和截距。

线性回归模型使用均方误差MSE作为损失函数（Loss），用以衡量预测房价和真实房价的差异，公式如下：

$$MSE=\frac{1}{n}\sum _{i=1}^n(\widehat{Y_i}-Y_i{)}^2(公式2)$$
**神经网络的标准结构中每个神经元由加权和与非线性变换构成，然后将多个神经元分层的摆放并连接形成神经网络。**线性回归模型可以认为是神经网络模型的一种极简特例，是一个只有加权和、没有非线性变换的神经元（无需形成网络），如 图2 所示。

图2：线性回归模型的神经网络结构

1.2 五步法

不同场景的深度学习模型具具备一定的通用性，五个步骤即可完成模型的构建和训练。在构建不同的模型时，只有模型三要素不同，其它步骤基本一致。

• 数据处理：从本地或URL读取数据，并完成预处理操作(如数据校验、格式转化等)，保证模型可读取。
• 模型设计: (模型要素1 )网络结构设计，相当于模型的假设空间，即模型能够表达的关系集合。
• 训练配置：(模型要素2 )设定模型采用的寻解算法，即优化器，并指定计算资源。
• 训练过程：(模型要素3 )循环调用训练过程，每轮都包括前向计算、损失函数(优化目标)和后向传播三个步骤。
• 模型保存：将训练好的模型保存，模型预测时调用。

图3：构建神经网络/深度学习模型的基本步骤

二、数据处理

2.1 数据导入

通过如下代码读入数据，了解下波士顿房价的数据集结构，数据存放在本地目录下housing.data文件中。

# 导入需要用到的package
import numpy as np
import json
import math 
# 读入训练数据
datafile = './work/housing.data'
data = np.fromfile(datafile, sep=' ')
data

2.2 数据形状变换

由于读入的原始数据是1维的，所有数据都连在一起。因此需要我们将数据的形状进行变换，形成一个2维的矩阵，每行为一个数据样本（14个值），每个数据样本包含13个XXX（影响房价的特征）和一个YYY（该类型房屋的均价）。

# 读入之后的数据被转化成1维array，其中array的第0-13项是第一条数据，第14-27项是第二条数据，以此类推.... 
# 这里对原始数据做reshape，变成N x 14的形式
feature_names = [ 'CRIM', 'ZN', 'INDUS', 'CHAS', 'NOX', 'RM', 'AGE','DIS', 
                 'RAD', 'TAX', 'PTRATIO', 'B', 'LSTAT', 'MEDV' ]
feature_num = len(feature_names)
data = data.reshape([data.shape[0] // feature_num, feature_num])

# 查看数据
x = data[0]
print(x.shape)
print(x)

2.3 数据集划分

**将数据集划分成训练集和测试集，其中训练集用于确定模型的参数，测试集用于评判模型的效果。**一般我们将80%的数据用作训练集，20%用作测试集，实现代码如下。通过打印训练集的形状，可以发现共有404个样本，每个样本含有13个特征和1个预测值。

ratio = 0.8
offset = int(data.shape[0] * ratio)
training_data = data[:offset]
training_data.shape

2.4 数据归一化处理

**对每个特征进行归一化处理，使得每个特征的取值缩放到0~1之间。**这样做有两个好处：一是模型训练更高效，二是特征前的权重大小可以代表该变量对预测结果的贡献度（因为每个特征值本身的范围相同）。

# 计算train数据集的最大值，最小值，平均值
maximums, minimums, avgs = \
                     training_data.max(axis=0), \
                     training_data.min(axis=0), \
     training_data.sum(axis=0) / training_data.shape[0]
# 对数据进行归一化处理
for i in range(feature_num):
    #print(maximums[i], minimums[i], avgs[i])
    data[:, i] = (data[:, i] - minimums[i]) / (maximums[i] - minimums[i])

2.5 封装成load data函数

将上述几个数据处理操作封装成load data函数，以便下一步模型的调用，实现方法如下。

def load_data():
    # 从文件导入数据
    datafile = './work/housing.data'
    data = np.fromfile(datafile, sep=' ')

    # 每条数据包括14项，其中前面13项是影响因素，第14项是相应的房屋价格中位数
    feature_names = [ 'CRIM', 'ZN', 'INDUS', 'CHAS', 'NOX', 'RM', 'AGE', \
                      'DIS', 'RAD', 'TAX', 'PTRATIO', 'B', 'LSTAT', 'MEDV' ]
    feature_num = len(feature_names)

    # 将原始数据进行Reshape，变成[N, 14]这样的形状
    data = data.reshape([data.shape[0] // feature_num, feature_num])

    # 将原数据集拆分成训练集和测试集
    # 这里使用80%的数据做训练，20%的数据做测试
    # 测试集和训练集必须是没有交集的
    ratio = 0.8
    offset = int(data.shape[0] * ratio)
    training_data = data[:offset]

    # 计算训练集的最大值，最小值，平均值
    maximums, minimums, avgs = training_data.max(axis=0), training_data.min(axis=0), \
                                 training_data.sum(axis=0) / training_data.shape[0]

    # 对数据进行归一化处理
    for i in range(feature_num):
        #print(maximums[i], minimums[i], avgs[i])
        data[:, i] = (data[:, i] - minimums[i]) / (maximums[i] - minimums[i])

    # 训练集和测试集的划分比例
    training_data = data[:offset]
    test_data = data[offset:]
    return training_data, test_data

2.6 获取归一化后的训练集和测试集

# 获取数据
training_data, test_data = load_data()
x = training_data[:, :-1]
y = training_data[:, -1:]

# 查看数据
print(x[0])
print(y[0])

三、模型设计

模型设计是深度学习模型关键要素之一，也称为网络结构设计，相当于模型的假设空间，即实现模型“前向计算”（从输入到输出）的过程。

3.1 前向计算

• 如果将输入特征和输出预测值均以向量表示，输入特征$x$有13个分量，$y$有1个分量，那么参数权重的形状（shape）是$13\times 1$。假设我们以如下任意数字赋值参数做初始化：

$$w=[0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,-0.1,-0.2,-0.3,-0.4,0.0]$$

w = [0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, -0.1, -0.2, -0.3, -0.4, 0.0]
w = np.array(w).reshape([13, 1])

• 取出第1条样本数据，观察样本的特征向量与参数向量相乘的结果。

x1=x[0]
t = np.dot(x1, w)
print(t)

• 完整的线性回归公式，还需要初始化偏移量$b$，同样随意赋初值-0.2。那么，线性回归模型的完整输出是$z=t+b$，这个从特征和参数计算输出值的过程称为“前向计算”。

b = -0.2
z = t + b
print(z)

3.2 以类的方式来实现网络结构

将上述计算预测输出的过程以“类和对象”的方式来描述，类成员变量有参数$w$和$b$。通过写一个forward函数（代表“前向计算”）完成上述从特征和参数到输出预测值的计算过程，代码如下所示。

class Network(object):
    def __init__(self, num_of_weights):
        # 随机产生w的初始值
        # 为了保持程序每次运行结果的一致性，
        # 此处设置固定的随机数种子
        np.random.seed(0)
        self.w = np.random.randn(num_of_weights, 1)
        self.b = 0.
        
    def forward(self, x):
        z = np.dot(x, self.w) + self.b
        return z

# 生成类的实例，调用其方法来完成前向计算的过程
net = Network(13)
x1 = x[0]
y1 = y[0]
z = net.forward(x1)
print(z)

从上述前向计算的过程可见，**线性回归也可以表示成一种简单的神经网络（只有一个神经元，且激活函数为恒等式）。**这也是机器学习模型普遍为深度学习模型替代的原因：由于深度学习网络强大的表示能力，很多传统机器学习模型的学习能力等同于相对简单的深度学习模型。

四、训练配置

模型设计完成后，需要通过训练配置寻找模型的最优值，即通过损失函数来衡量模型的好坏。训练配置也是深度学习模型关键要素之一。

4.1 损失函数——均方误差

通过模型计算$x_1$表示的影响因素所对应的房价应该是$z$, 但实际数据告诉我们房价是$y$。这时我们需要有某种指标来衡量预测值$z$跟真实值$y$之间的差距。对于回归问题，最常采用的衡量方法是使用均方误差作为评价模型好坏的指标，具体定义如下：

$$Loss=(y-z{)}^2(公式3)$$
**上式中的$Loss$（简记为: $L$）通常也被称作损失函数，它是衡量模型好坏的指标。**在回归问题中，均方误差是一种比较常见的形式，分类问题中通常会采用交叉熵作为损失函数，在后续的章节中会更详细的介绍。对一个样本计算损失函数值的实现如下。

Loss = (y1 - z)*(y1 - z)
print(Loss)

4.2 同时计算多个样本的损失函数

因为计算损失函数时需要把每个样本的损失函数值都考虑到，所以我们需要对单个样本的损失函数进行求和，并除以样本总数$N$。
$$L=\frac{1}{N}\sum _{i=1}^N(y_i-z_i{)}^2(公式4)$$
在Network类下面添加损失函数的计算过程如下。

class Network(object):
    def __init__(self, num_of_weights):
        # 随机产生w的初始值
        # 为了保持程序每次运行结果的一致性，此处设置固定的随机数种子
        np.random.seed(0)
        self.w = np.random.randn(num_of_weights, 1)
        self.b = 0.
        
    def forward(self, x):
        z = np.dot(x, self.w) + self.b
        return z
    
    # 求均方误差
    def loss(self, z, y):
        cost = np.mean((z - y)**2)
        return cost

使用定义的Network类，可以方便的计算预测值和损失函数。需要注意的是，类中的变量$x$, $w$，$b$, $z$, $error$等均是向量。以变量$x$为例，共有两个维度，一个代表特征数量（值为13），一个代表样本数量，代码如下所示。

net = Network(13)
# 此处可以一次性计算多个样本的预测值和损失函数
x1 = x[0:3]
y1 = y[0:3]
z = net.forward(x1)
print('predict: ', z)
loss = net.loss(z, y1)
print('loss:', loss)

五、训练过程

上述计算过程描述了如何构建神经网络，通过神经网络完成预测值和损失函数的计算。接下来介绍如何**求解参数$w$和$b$的数值，这个过程也称为模型训练过程。**训练过程是深度学习模型的关键要素之一，其目标是让定义的损失函数$Loss$尽可能的小，也就是说找到一个参数解$w$和$b$，使得损失函数取得极小值。

其中$L$表示的是损失函数的值，$\mathbf{w}$为模型权重，$b$为偏置项。$\mathbf{w}$和$b$均为要学习的模型参数。
把损失函数表示成矩阵的形式为
$$L=\frac{1}{N}||\mathbf{y}-(\mathbf{X}\mathbf{w}+\mathbf{b})|{|}^2,(公式7)$$

其中$\mathbf{y}$为$N$个样本的标签值构成的向量，形状为$N\times 1$；$\mathbf{X}$为$N$个样本特征向量构成的矩阵，形状为$N\times D$，$D$为数据特征长度；$\mathbf{w}$为权重向量，形状为$D\times 1$；$\mathbf{b}$为所有元素都为$b$的向量，形状为$N\times 1$。
计算公式7对参数$b$的偏导数
$$\frac{\partial L}{\partial b}=1^T(\mathbf{y}-(\mathbf{X}\mathbf{w}+\mathbf{b})),(公式8)$$
请注意，上述公式忽略了系数$\frac{2}{N}$，并不影响最后结果。其中$1$为$N$维的全1向量。
令公式8等于0，得到
$$b^{\ast }={\bar{\mathbf{x}}}^T\mathbf{w}-\bar{y}(公式9)$$
**其中$\bar{y}=\frac{1}{N}{1}^T\mathbf{y}$为所有标签的平均值，$\bar{\mathbf{x}}=\frac{1}{N}(1^T\mathbf{X}{)}^T$为所有特征向量的平均值。**将$b^{\ast }$带入公式7中并对参数$\mathbf{w}$求偏导得到
$$\frac{\partial L}{\partial \mathbf{w}}=(\mathbf{X}-{\bar{\mathbf{x}}}^T{)}^T((\mathbf{y}-\bar{y})-(\mathbf{X}-{\bar{\mathbf{x}}}^T)\mathbf{w})(公式10)$$

令公式10等于0，得到最优参数

$$\begin{gathered} {\mathbf{w}}^{\ast }=((\mathbf{X}-{\bar{\mathbf{x}}}^T{)}^T(\mathbf{X}-{\bar{\mathbf{x}}}^T){)}^{-1}(\mathbf{X}-{\bar{\mathbf{x}}}^T{)}^T(\mathbf{y}-\bar{y})(公式11) \\ {b}^{\ast }={\bar{\mathbf{x}}}^T{\mathbf{w}}^{\ast }-\bar{y}(公式12) \end{gathered}$$

将样本数据$(x,y)$带入上面的公式11和公式12中即可求解出$w$和$b$的值，但是这种方法只对线性回归这样简单的任务有效。如果模型中含有非线性变换，或者损失函数不是均方差这种简单的形式，则很难通过上式求解。为了解决这个问题，下面我们将引入更加普适的数值求解方法：梯度下降法。

5.1 梯度下降法

**在现实中存在大量的函数正向求解容易，但反向求解较难，被称为单向函数，这种函数在密码学中有大量的应用。**密码锁的特点是可以迅速判断一个密钥是否是正确的(已知$x$，求$y$很容易)，但是即使获取到密码锁系统，也无法破解出正确得密钥（已知$y$，求$x$很难）。

这种情况特别类似于一位想从山峰走到坡谷的盲人，他看不见坡谷在哪（无法逆向求解出$Loss$导数为0时的参数值），但可以伸脚探索身边的坡度（当前点的导数值，也称为梯度）。那么，求解Loss函数最小值可以这样实现：从当前的参数取值，一步步的按照下坡的方向下降，直到走到最低点。这种方法笔者称它为“盲人下坡法”。哦不，有个更正式的说法“梯度下降法”。

训练的关键是找到一组$(w,b)$，使得损失函数$L$取极小值。我们先看一下损失函数$L$只随两个参数$w_5$、$w_9$变化时的简单情形，启发下寻解的思路。
$$L=L(w_5,w_9)(公式13)$$
这里将$w_0,w_1,...,w_{12}$中除$w_5,w_9$之外的参数和$b$都固定下来，可以用图画出$L(w_5,w_9)$的形式，并在三维空间中画出损失函数随参数变化的曲面图。

net = Network(13)
losses = []
#只画出参数w5和w9在区间[-160, 160]的曲线部分，以及包含损失函数的极值
w5 = np.arange(-160.0, 160.0, 1.0)	#[-160,-159 , ... , 159],共320个数
w9 = np.arange(-160.0, 160.0, 1.0)
losses = np.zeros([len(w5), len(w9)])

#计算设定区域内每个参数取值所对应的Loss
for i in range(len(w5)):
    for j in range(len(w9)):
        net.w[5] = w5[i]
        net.w[9] = w9[j]
        z = net.forward(x)
        loss = net.loss(z, y)
        losses[i, j] = loss

#使用matplotlib将两个变量和对应的Loss作3D图
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
fig = plt.figure()
ax = Axes3D(fig)

# 生成网格点坐标矩阵
w5, w9 = np.meshgrid(w5, w9)

# 绘制3D图，rstride（row）指定行的跨度，设置颜色映射
ax.plot_surface(w5, w9, losses, rstride=1, cstride=1, cmap='rainbow')
plt.show()

观察上述曲线呈现出“圆滑”的坡度，这正是我们选择以均方误差作为损失函数的原因之一。图6 呈现了只有一个参数维度时，均方误差和绝对值误差（只将每个样本的误差累加，不做平方处理）的损失函数曲线图。

图6：均方误差和绝对值误差损失函数曲线图

由此可见，均方误差表现的“圆滑”的坡度有两个好处：

• 曲线的最低点是可导的。
• 越接近最低点，曲线的坡度逐渐放缓，有助于通过当前的梯度来判断接近最低点的程度（是否逐渐减少步长，以免错过最低点）。

现在我们要找出一组$[w_5,w_9]$的值，使得损失函数最小，实现梯度下降法的方案如下：

• 步骤1：随机的选一组初始值，例如：$[w_5,w_9]=[-100.0,-100.0]$
• 步骤2：选取下一个点$[w_5^{{}^{\prime }},w_9^{{}^{\prime }}]$，使得$L(w_5^{{}^{\prime }},w_9^{{}^{\prime }})<L(w_5,w_9)$
• 步骤3：重复步骤2，直到损失函数几乎不再下降。

如何选择$[w_5^{{}^{\prime }},w_9^{{}^{\prime }}]$是至关重要的，第一要保证$L$是下降的，第二要使得下降的趋势尽可能的快。微积分的基础知识告诉我们：**沿着梯度的反方向，是函数值下降最快的方向。**简单理解，函数在某一个点的梯度方向是曲线斜率最大的方向，但梯度方向是向上的，所以下降最快的是梯度的反方向。

5.2 梯度计算

上文已经介绍了损失函数的计算方法，这里稍微改写。为了使梯度计算更加简洁，引入因子$\frac{1}{2}$，定义损失函数如下：

$$L=\frac{1}{2N}\sum _{i=1}^N(y_i-z_i{)}^2(公式14)$$
其中$z_i$是网络对第$i$个样本的预测值：

$$z_i=\sum _{j=0}^{12}{x}_i^j\cdot w_j+b(公式15)$$
梯度的定义：

$$gradient=(\frac{\partial L}{\partial w_0},\frac{\partial L}{\partial w_1},...,\frac{\partial L}{\partial w_{12}},\frac{\partial L}{\partial b})(公式16)$$
可以

计算出$L$对$w$和$b$的偏导数：
$$\frac{\partial L}{\partial w_j}=\frac{1}{N}\sum _{i=1}^N(z_i-y_i)\frac{\partial z_i}{\partial w_j}=\frac{1}{N}\sum _{i=1}^N(z_i-y_i)x_i^j(公式17)$$

$$\frac{\partial L}{\partial b}=\frac{1}{N}\sum _{i=1}^N(z_i-y_i)\frac{\partial z_i}{\partial b}=\frac{1}{N}\sum _{i=1}^N(z_i-y_i)(公式18)$$

从导数的计算过程可以看出，因子$\frac{1}{2}$被消掉了，这是因为二次函数求导的时候会产生因子$2$，这也是我们将损失函数改写的原因。

下面我们考虑只有一个样本的情况下，计算梯度：

$$L=\frac{1}{2}(y_i-z_i{)}^2(公式19)$$

$$z_1=x_1^0\cdot w_0+x_1^1\cdot w_1+...+x_1^{12}\cdot w_{12}+b(公式20)$$

可以计算出：

$$L=\frac{1}{2}(x_1^0\cdot w_0+x_1^1\cdot w_1+...+x_1^{12}\cdot w_{12}+b-y_1{)}^2(公式21)$$
可以计算出$L$对$w$和$b$的偏导数：

$$\frac{\partial L}{\partial w_0}=(x_1^0\cdot w_0+x_1^1\cdot w_1+...+x_1^{12}\cdot w_{1}2+b-y_1)\cdot x_1^0=(z_1-y_1)\cdot x_1^0(公式22)$$

$$\frac{\partial L}{\partial b}=(x_1^0\cdot w_0+x_1^1\cdot w_1+...+x_1^{12}\cdot w_{12}+b-y_1)\cdot 1=(z_1-y_1)(公式23)$$

x1 = x[0]
y1 = y[0]
z1 = net.forward(x1)
print('x1 {}, shape {}'.format(x1, x1.shape))
print('y1 {}, shape {}'.format(y1, y1.shape))
print('z1 {}, shape {}'.format(z1, z1.shape))

# 按上面的公式，当只有一个样本时，可以计算某个w_j，比如w_0的梯度。
gradient_w0 = (z1 - y1) * x1[0]
print('gradient_w0 {}'.format(gradient_w0))

# 同样我们可以计算w1的梯度
gradient_w1 = (z1 - y1) * x1[1]
print('gradient_w1 {}'.format(gradient_w1))

# 依次计算w_2的梯度
gradient_w2= (z1 - y1) * x1[2]
print('gradient_w1 {}'.format(gradient_w2))

# 因此可以写一个for循环即可计算从 w0到w_12 的所有权重的梯度

5.3 使用NumPy进行梯度计算

基于NumPy广播机制（对向量和矩阵计算如同对1个单一变量计算一样），可以更快速的实现梯度计算。计算梯度的代码中直接用$(z_1-y_1)\cdot x_1$，得到的是一个13维的向量，每个分量分别代表该维度的梯度。

gradient_w = (z1 - y1) * x1
print('gradient_w_by_sample1 {}, gradient.shape {}'.format(gradient_w, gradient_w.shape))

输入数据中有多个样本，每个样本都对梯度有贡献。如上代码计算了只有样本1时的梯度值，同样的计算方法也可以计算样本2和样本3对梯度的贡献。

x2 = x[1]
y2 = y[1]
z2 = net.forward(x2)
gradient_w = (z2 - y2) * x2
print('gradient_w_by_sample2 {}, gradient.shape {}'.format(gradient_w, gradient_w.shape))

可能有的读者再次想到可以使用for循环把每个样本对梯度的贡献都计算出来，然后再作平均。但是我们不需要这么做，仍然可以使用NumPy的矩阵操作来简化运算，如3个样本的情况。

# 注意这里是一次取出3个样本的数据，不是取出第3个样本
x3samples = x[0:3]
y3samples = y[0:3]
z3samples = net.forward(x3samples)

print('x {}, shape {}'.format(x3samples, x3samples.shape))
print('y {}, shape {}'.format(y3samples, y3samples.shape))
print('z {}, shape {}'.format(z3samples, z3samples.shape))

# 下面计算这3个样本对梯度的贡献。
gradient_w = (z3samples - y3samples) * x3samples
print('gradient_w {}, gradient.shape {}'.format(gradient_w, gradient_w.shape))

# 计算这N个样本对梯度的贡献。
z = net.forward(x)
gradient_w = (z - y) * x
print('gradient_w shape {}'.format(gradient_w.shape))
print(gradient_w)

**此处可见，计算梯度gradient_w的维度是$3\times 13$，**并且其第1行与上面第1个样本计算的梯度gradient_w_by_sample1一致，第2行与上面第2个样本计算的梯度gradient_w_by_sample2一致，第3行与上面第3个样本计算的梯度gradient_w_by_sample3一致。这里使用矩阵操作，可以更加方便的对3个样本分别计算各自对梯度的贡献。

那么对于有N个样本的情形，我们可以直接使用如下方式计算出所有样本对梯度的贡献，这就是使用NumPy库广播功能带来的便捷。
小结一下这里使用NumPy库的广播功能：

• 一方面可以扩展参数的维度，代替for循环来计算1个样本对从$w_0$到$w_{12}$的所有参数的梯度。
• 另一方面可以扩展样本的维度，代替for循环来计算样本0到样本403对参数的梯度。

上面gradient_w的每一行代表了一个样本对梯度的贡献。根据梯度的计算公式，总梯度是对每个样本对梯度贡献的平均值。
$$\frac{\partial L}{\partial w_j}=\frac{1}{N}\sum _{i=1}^N(z_i-y_i)\frac{\partial z_i}{\partial w_j}=\frac{1}{N}\sum _{i=1}^N(z_i-y_i)x_i^j(公式24)$$
可以使用NumPy的均值函数来完成此过程，代码实现如下。

# axis = 0 表示把每一行做相加然后再除以总的行数
gradient_w = np.mean(gradient_w, axis=0)
print('gradient_w ', gradient_w.shape)
print('w ', net.w.shape)
print(gradient_w)
print(net.w)

使用NumPy的矩阵操作方便地完成了gradient的计算，但引入了一个问题，gradient_w的形状是(13,)，而$w$的维度是(13, 1)。导致该问题的原因是使用np.mean函数时消除了第0维。为了加减乘除等计算方便，gradient_w和$w$必须保持一致的形状。因此我们将gradient_w的维度也设置为(13,1)，代码如下：

gradient_w = gradient_w[:, np.newaxis]
print('gradient_w shape', gradient_w.shape)

综合上面的剖析，计算梯度的代码如下所示。

z = net.forward(x)
gradient_w = (z - y) * x
gradient_w = np.mean(gradient_w, axis=0)
gradient_w = gradient_w[:, np.newaxis]
gradient_w

上述代码非常简洁地完成了$w$的梯度计算。同样，计算$b$的梯度的代码也是类似的原理。

gradient_b = (z - y)
gradient_b = np.mean(gradient_b)
# 此处b是一个数值，所以可以直接用np.mean得到一个标量
gradient_b

将上面计算$w$和$b$的梯度的过程，写成Network类的gradient函数，实现方法如下所示。

class Network(object):
    def __init__(self, num_of_weights):
        # 随机产生w的初始值
        # 为了保持程序每次运行结果的一致性，此处设置固定的随机数种子
        np.random.seed(0)
        self.w = np.random.randn(num_of_weights, 1)
        self.b = 0.
        
    def forward(self, x):
        z = np.dot(x, self.w) + self.b
        return z
    
    def loss(self, z, y):
        error = z - y
        num_samples = error.shape[0]
        cost = error * error
        cost = np.sum(cost) / num_samples
        return cost
    
    def gradient(self, x, y):
        z = self.forward(x)
        gradient_w = (z-y)*x
        gradient_w = np.mean(gradient_w, axis=0)
        gradient_w = gradient_w[:, np.newaxis]
        gradient_b = (z - y)
        gradient_b = np.mean(gradient_b)
        
        return gradient_w, gradient_b

# 调用上面定义的gradient函数，计算梯度
# 初始化网络
net = Network(13)
# 设置[w5, w9] = [-100., -100.]
net.w[5] = -100.0
net.w[9] = -100.0

z = net.forward(x)
loss = net.loss(z, y)
gradient_w, gradient_b = net.gradient(x, y)
gradient_w5 = gradient_w[5][0]
gradient_w9 = gradient_w[9][0]
print('point {}, loss {}'.format([net.w[5][0], net.w[9][0]], loss))
print('gradient {}'.format([gradient_w5, gradient_w9]))

5.4 梯度更新

下面研究更新梯度的方法，确定损失函数更小的点。首先沿着梯度的反方向移动一小步，找到下一个点P1，观察损失函数的变化。

# 在[w5, w9]平面上，沿着梯度的反方向移动到下一个点P1
# 定义移动步长 eta
eta = 0.1
# 更新参数w5和w9
net.w[5] = net.w[5] - eta * gradient_w5
net.w[9] = net.w[9] - eta * gradient_w9
# 重新计算z和loss
z = net.forward(x)
loss = net.loss(z, y)
gradient_w, gradient_b = net.gradient(x, y)
gradient_w5 = gradient_w[5][0]
gradient_w9 = gradient_w[9][0]
print('point {}, loss {}'.format([net.w[5][0], net.w[9][0]], loss))
print('gradient {}'.format([gradient_w5, gradient_w9]))

运行上面的代码，可以发现沿着梯度反方向走一小步，下一个点的损失函数的确减少了。感兴趣的话，大家可以尝试不停的点击上面的代码块，观察损失函数是否一直在变小。

在上述代码中，每次更新参数使用的语句：
net.w[5] = net.w[5] - eta * gradient_w5

• 相减：参数需要向梯度的反方向移动。
• eta：控制每次参数值沿着梯度反方向变动的大小，即每次移动的步长，又称为学习率。

大家可以思考下，为什么之前我们要做输入特征的归一化，保持尺度一致？这是为了让统一的步长更加合适，使训练更加高效。

如 图8 所示，特征输入归一化后，不同参数输出的Loss是一个比较规整的曲线，学习率可以设置成统一的值 ；特征输入未归一化时，不同特征对应的参数所需的步长不一致，尺度较大的参数需要大步长，尺寸较小的参数需要小步长，导致无法设置统一的学习率。

图8：未归一化的特征，会导致不同特征维度的理想步长不同

5.5 封装Train函数

将上面的循环计算过程封装在train和update函数中，实现方法如下所示。

class Network(object):
    def __init__(self, num_of_weights):
        # 随机产生w的初始值
        # 为了保持程序每次运行结果的一致性，此处设置固定的随机数种子
        np.random.seed(0)
        self.w = np.random.randn(num_of_weights,1)
        self.w[5] = -100.
        self.w[9] = -100.
        self.b = 0.
        
    def forward(self, x):
        z = np.dot(x, self.w) + self.b
        return z
    
    def loss(self, z, y):
        error = z - y
        num_samples = error.shape[0]
        cost = error * error
        cost = np.sum(cost) / num_samples
        return cost
    
    def gradient(self, x, y):
        z = self.forward(x)
        gradient_w = (z-y)*x
        gradient_w = np.mean(gradient_w, axis=0)
        gradient_w = gradient_w[:, np.newaxis]
        gradient_b = (z - y)
        gradient_b = np.mean(gradient_b)        
        return gradient_w, gradient_b
    
    def update(self, gradient_w5, gradient_w9, eta=0.01):
        net.w[5] = net.w[5] - eta * gradient_w5
        net.w[9] = net.w[9] - eta * gradient_w9
        
    def train(self, x, y, iterations=100, eta=0.01):
        points = []
        losses = []
        for i in range(iterations):
            points.append([net.w[5][0], net.w[9][0]])
            z = self.forward(x)
            L = self.loss(z, y)
            gradient_w, gradient_b = self.gradient(x, y)
            gradient_w5 = gradient_w[5][0]
            gradient_w9 = gradient_w[9][0]
            self.update(gradient_w5, gradient_w9, eta)
            losses.append(L)
            if i % 50 == 0:
                print('iter {}, point {}, loss {}'.format(i, [net.w[5][0], net.w[9][0]], L))
        return points, losses

# 获取数据
train_data, test_data = load_data()
x = train_data[:, :-1]
y = train_data[:, -1:]
# 创建网络
net = Network(13)
num_iterations=2000
# 启动训练
points, losses = net.train(x, y, iterations=num_iterations, eta=0.01)

# 画出损失函数的变化趋势
plot_x = np.arange(num_iterations)
plot_y = np.array(losses)
plt.plot(plot_x, plot_y)
plt.show()

5.6 训练过程扩展到全部参数

为了能给读者直观的感受，上文演示的梯度下降的过程仅包含$w_5$和$w_9$两个参数。但房价预测的模型必须要对所有参数$w$和$b$进行求解，这需要将Network中的update和train函数进行修改。由于不再限定参与计算的参数（所有参数均参与计算），修改之后的代码反而更加简洁。

实现逻辑：“前向计算输出、根据输出和真实值计算Loss、基于Loss和输入计算梯度、根据梯度更新参数值”四个部分反复执行，直到到损失函数最小。

class Network(object):
    def __init__(self, num_of_weights):
        # 随机产生w的初始值
        # 为了保持程序每次运行结果的一致性，此处设置固定的随机数种子
        np.random.seed(0)
        self.w = np.random.randn(num_of_weights, 1)
        self.b = 0.
        
    def forward(self, x):
        z = np.dot(x, self.w) + self.b
        return z
    
    def loss(self, z, y):
        error = z - y
        num_samples = error.shape[0]
        cost = error * error
        cost = np.sum(cost) / num_samples
        return cost
    
    def gradient(self, x, y):
        z = self.forward(x)
        gradient_w = (z-y)*x
        gradient_w = np.mean(gradient_w, axis=0)
        gradient_w = gradient_w[:, np.newaxis]
        gradient_b = (z - y)
        gradient_b = np.mean(gradient_b)        
        return gradient_w, gradient_b
    
    def update(self, gradient_w, gradient_b, eta = 0.01):
        self.w = self.w - eta * gradient_w
        self.b = self.b - eta * gradient_b
        
    def train(self, x, y, iterations=100, eta=0.01):
        losses = []
        for i in range(iterations):
            z = self.forward(x)
            L = self.loss(z, y)
            gradient_w, gradient_b = self.gradient(x, y)
            self.update(gradient_w, gradient_b, eta)
            losses.append(L)
            if (i+1) % 10 == 0:
                print('iter {}, loss {}'.format(i, L))
        return losses

# 获取数据
train_data, test_data = load_data()
x = train_data[:, :-1]
y = train_data[:, -1:]
# 创建网络
net = Network(13)
num_iterations=1000
# 启动训练
losses = net.train(x,y, iterations=num_iterations, eta=0.01)

# 画出损失函数的变化趋势
plot_x = np.arange(num_iterations)
plot_y = np.array(losses)
plt.plot(plot_x, plot_y)
plt.show()

5.7 随机梯度下降法（ Stochastic Gradient Descent）

在上述程序中，每次损失函数和梯度计算都是基于数据集中的全量数据。对于波士顿房价预测任务数据集而言，样本数比较少，只有404个。但在实际问题中，数据集往往非常大，如果每次都使用全量数据进行计算，效率非常低，通俗地说就是“杀鸡焉用牛刀”。由于参数每次只沿着梯度反方向更新一点点，因此方向并不需要那么精确。一个合理的解决方案是每次从总的数据集中随机抽取出小部分数据来代表整体，基于这部分数据计算梯度和损失来更新参数，这种方法被称作随机梯度下降法（Stochastic Gradient Descent，SGD），核心概念如下：

• mini-batch：每次迭代时抽取出来的一批数据被称为一个mini-batch。
• batch_size：一个mini-batch所包含的样本数目称为batch_size。
• epoch：当程序迭代的时候，按mini-batch逐渐抽取出样本，当把整个数据集都遍历到了的时候，则完成了一轮训练，也叫一个epoch。启动训练时，可以将训练的轮数num_epochs和batch_size作为参数传入。

下面结合程序介绍具体的实现过程，涉及到数据处理和训练过程两部分代码的修改。

• 数据处理代码修改

数据处理需要实现拆分数据批次和样本乱序（为了实现随机抽样的效果）两个功能。

# 获取数据
train_data, test_data = load_data()

# 打乱样本顺序
np.random.shuffle(train_data)

# 将train_data分成多个mini_batch
batch_size = 10
n = len(train_data)
mini_batches = [train_data[k:k+batch_size] for k in range(0, n, batch_size)]

# 创建网络
net = Network(13)

# 依次使用每个mini_batch的数据
for mini_batch in mini_batches:
    x = mini_batch[:, :-1]
    y = mini_batch[:, -1:]
    loss = net.train(x, y, iterations=1)

• 训练过程代码修改

将每个随机抽取的mini-batch数据输入到模型中用于参数训练。训练过程的核心是两层循环：

1. 第一层循环，代表样本集合要被训练遍历几次，称为“epoch”，代码如下：

for epoch_id in range(num_epochs):

2. 第二层循环，代表每次遍历时，样本集合被拆分成的多个批次，需要全部执行训练，称为“iter (iteration)”，代码如下：

for iter_id,mini_batch in emumerate(mini_batches):

在两层循环的内部是经典的四步训练流程：前向计算->计算损失->计算梯度->更新参数，这与大家之前所学是一致的，代码如下：

            x = mini_batch[:, :-1]
            y = mini_batch[:, -1:]
            a = self.forward(x)  #前向计算
            loss = self.loss(a, y)  #计算损失
            gradient_w, gradient_b = self.gradient(x, y)  #计算梯度
            self.update(gradient_w, gradient_b, eta)  #更新参数

将两部分改写的代码集成到Network类中的train函数中，最终的实现如下。

import numpy as np

class Network(object):
    def __init__(self, num_of_weights):
        # 随机产生w的初始值
        # 为了保持程序每次运行结果的一致性，此处设置固定的随机数种子
        #np.random.seed(0)
        self.w = np.random.randn(num_of_weights, 1)
        self.b = 0.
        
    def forward(self, x):
        z = np.dot(x, self.w) + self.b
        return z
    
    def loss(self, z, y):
        error = z - y
        num_samples = error.shape[0]
        cost = error * error
        cost = np.sum(cost) / num_samples
        return cost
    
    def gradient(self, x, y):
        z = self.forward(x)
        N = x.shape[0]
        gradient_w = 1. / N * np.sum((z-y) * x, axis=0)
        gradient_w = gradient_w[:, np.newaxis]
        gradient_b = 1. / N * np.sum(z-y)
        return gradient_w, gradient_b
    
    def update(self, gradient_w, gradient_b, eta = 0.01):
        self.w = self.w - eta * gradient_w
        self.b = self.b - eta * gradient_b
            
                
    def train(self, training_data, num_epochs, batch_size=10, eta=0.01):
        n = len(training_data)
        losses = []
        for epoch_id in range(num_epochs):
            # 在每轮迭代开始之前，将训练数据的顺序随机打乱
            # 然后再按每次取batch_size条数据的方式取出
            np.random.shuffle(training_data)
            # 将训练数据进行拆分，每个mini_batch包含batch_size条的数据
            mini_batches = [training_data[k:k+batch_size] for k in range(0, n, batch_size)]
            for iter_id, mini_batch in enumerate(mini_batches):
                x = mini_batch[:, :-1]
                y = mini_batch[:, -1:]
                a = self.forward(x)
                loss = self.loss(a, y)
                gradient_w, gradient_b = self.gradient(x, y)
                self.update(gradient_w, gradient_b, eta)
                losses.append(loss)
                print('Epoch {:3d} / iter {:3d}, loss = {:.4f}'.
                                 format(epoch_id, iter_id, loss))
        
        return losses

# 获取数据
train_data, test_data = load_data()

# 创建网络
net = Network(13)
# 启动训练
losses = net.train(train_data, num_epochs=50, batch_size=100, eta=0.1)

# 画出损失函数的变化趋势
plot_x = np.arange(len(losses))
plot_y = np.array(losses)
plt.plot(plot_x, plot_y)
plt.show()

观察上述Loss的变化，随机梯度下降加快了训练过程，但由于每次仅基于少量样本更新参数和计算损失，所以损失下降曲线会出现震荡。由于房价预测的数据量过少，所以难以感受到随机梯度下降带来的性能提升。

六、模型保存

Numpy提供了save接口，可直接将模型权重数组保存为.npy格式的文件。

np.save('w.npy', net.w)
np.save('b.npy', net.b)

小结

本节我们详细介绍了如何使用NumPy实现梯度下降算法，构建并训练了一个简单的线性模型实现波士顿房价预测，可以总结出，使用神经网络建模房价预测有三个要点：

• 构建网络，初始化参数$w$和$b$，定义预测和损失函数的计算方法。
• 随机选择初始点，建立梯度的计算方法和参数更新方式。
• 从总的数据集中抽取部分数据作为一个mini_batch，计算梯度并更新参数，不断迭代直到损失函数几乎不再下降。