神经网络入门之矢量化

矢量化

---

这部分教程将介绍三部分：

• 矢量的反向传播
• 梯度检查
• 动量

在先前的教程中，我们已经使用学习了一个非常简单的神经网络：一个输入数据，一个隐藏神经元和一个输出结果。在这篇教程中，我们将描述一个稍微复杂一点的神经网络：包括一个二维的输入数据，三维的隐藏神经元和二维的输出结果，并且利用softmax函数来做最后的分类。在之前的网络中，我们都没有添加偏差项，但在这个网络模型中我们加入偏差项。网络模型如下：

神经网络模型

我们先导入教程需要使用的软件包。

import numpy as np  import sklearn.datasets import matplotlib.pyplot as plt  from matplotlib.colors import colorConverter, ListedColormap  from mpl_toolkits.mplot3d import Axes3D  from matplotlib import cm

定义数据集

在这个教程中，我们的分类目标还是二分类：红色类别（t=1）和蓝色类别（t=0）。红色样本是一个环形分布，位于环形中间的是蓝色样本。我们利用scikit-learn的make_circles的方法得到数据集。

这是一个二维的数据集，而且不是线性分割。如果我们利用教程2中的方法，那么我们不能正确将它们分类，因为教程2中的方法只能学习出线性分割的样本数据。但我们增加一个隐藏层，就能学习这个非线性分割了。

我们有N个输入数据，每个数据有两个特征，那么我们可以得到矩阵X如下：

输入数据

其中，xij表示第i个样本的第j个特征值。

经过softmax函数之后，该模型输出的最终结果T为：

输出结果

其中，当且仅当第i个样本属于类别j时，tij=1。因此，我们定义蓝色样本的标记是T = [0 1]，红色样本的标记是T = [1 0]。

# Generate the datasetX, t = sklearn.datasets.make_circles(n_samples=100, shuffle=False, factor=0.3, noise=0.1)T = np.zeros((100,2)) # Define target matrixT[t==1,1] = 1T[t==0,0] = 1# Separate the red and blue points for plottingx_red = X[t==0]x_blue = X[t==1]print('shape of X: {}'.format(X.shape))print('shape of T: {}'.format(T.shape))

shape of X: (100, 2)
shape of T: (100, 2)

# Plot both classes on the x1, x2 planeplt.plot(x_red[:,0], x_red[:,1], 'ro', label='class red')plt.plot(x_blue[:,0], x_blue[:,1], 'bo', label='class blue')plt.grid()plt.legend(loc=1)plt.xlabel('$x_1$', fontsize=15)plt.ylabel('$x_2$', fontsize=15)plt.axis([-1.5, 1.5, -1.5, 1.5])plt.title('red vs blue classes in the input space')plt.show()

样本分布

矢量的反向传播

1. 矢量的正向传播

计算隐藏层的激活函数

将二维的输入样本X转换成三维的输入层H，我们需要使用链接矩阵Wh（Whij表示第i个输入特征和第j个隐藏层神经元相连的权重）和偏差向量bh：

参数

之后计算结果如下：

计算结果

其中，σ是一个Logistic函数，H是一个N*3的输出矩阵。

整个计算流过程如下图所示。每个输入特征xij和权重参数whj1、whj2、whj3相乘。最后相乘的每行k（xij*whjk）结果进行累加得到zik，之后经过Logistic函数σ进行非线性激活得到hik。注意，偏差项Bh可以被整合在链接矩阵Wh中，只需要通过在输入参数xi中增加一个+1项。

在代码中，Bh和Wh使用bh和Wh表示。hidden_activations(X, Wh, bh)函数实现了隐藏层的激活输出。

计算输出层的激活结果

在计算输出层的激活结果之前，我们需要先定义3*2的链接矩阵Wo（Woij表示隐藏层第i个神经元和输出层第j个输出单元的链接权重）和2*1的偏差项矩阵bo：

参数

之后计算结果如下：

其中，ς表示softmax函数。Y表示最后输出的n*2的矩阵，Zod表示矩阵Zo的第d列，Wod表示矩阵Wo的第d列，bod表示向量bo的第d个元素。

在代码中，bo和Wo用变量bo和Wo来表示。output_activations(H, Wo, bo)函数实现了输出层的激活结果。

# Define the logistic functiondef logistic(z):     return 1 / (1 + np.exp(-z))# Define the softmax functiondef softmax(z):     return np.exp(z) / np.sum(np.exp(z), axis=1, keepdims=True)# Function to compute the hidden activationsdef hidden_activations(X, Wh, bh):    return logistic(X.dot(Wh) + bh)# Define output layer feedforwarddef output_activations(H, Wo, bo):    return softmax(H.dot(Wo) + bo)# Define the neural network functiondef nn(X, Wh, bh, Wo, bo):     return output_activations(hidden_activations(X, Wh, bh), Wo, bo)# Define the neural network prediction function that only returns#  1 or 0 depending on the predicted classdef nn_predict(X, Wh, bh, Wo, bo):     return np.around(nn(X, Wh, bh, Wo, bo))

2. 矢量的反向传播

输出层的矢量反向传播

计算输出层的误差

最后一层的输出层使用的是softmax函数，该函数的交叉熵损失函数在这篇教程中已经有了很详细的描述。如果需要对N个样本进行C个分类，那么它的损失函数ξ是：

交叉熵损失函数

损失函数的误差梯度δo可以非常方便得到：

误差梯度

其中，Zo（Zo=H⋅Wo+bo）是一个n*2的矩阵，T是一个n*2的目标矩阵，Y是一个经过模型得到的n*2的输出矩阵。因此，δo也是一个n*2的矩阵。

在代码中，δo用Eo表示，error_output(Y, T)函数实现了该方法。

更新输出层的权重

对于N个样本，对输出层的梯度δwoj是通过∂ξ/∂woj计算的，具体计算如下：

其中，woj表示Wo的第j行，即是一个1*2的向量。因此，我们可以将上式改写成一个矩阵操作，即：

最后梯度的结果是一个3*2的Jacobian矩阵，如下：

在代码中，JWo表示上面的Jwo。gradient_weight_out(H, Eo)函数实现了上面的操作。

更新输出层的偏差项

对于偏差项bo可以采用相同的方式进行更新。对于批处理的N个样本，对输出层的梯度∂ξ/∂bo的计算如下：

最后梯度的结果是一个2*1的Jacobian矩阵，如下：

在代码中，Jbo表示上面的Jbo。gradient_bias_out(Eo)函数实现了上面的操作。

# Define the cost functiondef cost(Y, T):    return - np.multiply(T, np.log(Y)).sum()# Define the error function at the outputdef error_output(Y, T):    return Y - T# Define the gradient function for the weight parameters at the output layerdef gradient_weight_out(H, Eo):     return  H.T.dot(Eo)# Define the gradient function for the bias parameters at the output layerdef gradient_bias_out(Eo):     return  np.sum(Eo, axis=0, keepdims=True)

隐藏层的矢量化反馈

计算隐藏层的误差

在隐藏层上面的误差函数的梯度δh可以被定义如下：

其中，Zh是一个n*3的输入到Logistic函数中的输入矩阵，即：

其中，δh也是一个N*3的矩阵。

接下来，对于每个样本i和每个神经元j，我们计算误差梯度δhij的导数。通过前面一层的误差反馈，通过BP算法我们可以计算如下：

其中，woj表示Wo的第j行，它是一个1*2的向量。δoi也是一个1*2的向量。因此，维度是N*3的误差矩阵δh可以被如下计算：

其中，∘表示逐点乘积

在代码中，Eh表示上面的δh。error_hidden(H, Wo, Eo)函数实现了上面的函数。

更新隐藏层的权重

在N个样本中，隐藏层的梯度∂ξ/∂whj可以被如下计算：

其中，whj表示Wh的第j行，它是一个1*3的向量。我们可以将上面的式子写成矩阵相乘的形式如下：

梯度最后的结果是一个2*3的Jacobian矩阵，如下：

在代码中，JWh表示Jwh。gradient_weight_hidden(X, Eh)函数实现了上面的操作。

更新隐藏层的偏差项

偏差项bh可以按照相同的方式进行更新操作。在N个样本上面，梯度∂ξ/∂bh可以如下计算：

最后的梯度结果是一个1*3的Jacobian矩阵，如下：

在代码中，Jbh表示Jbh。gradient_bias_hidden(Eh)函数实现了上面的方法。

# Define the error function at the hidden layerdef error_hidden(H, Wo, Eo):    # H * (1-H) * (E . Wo^T)    return np.multiply(np.multiply(H,(1 - H)), Eo.dot(Wo.T))# Define the gradient function for the weight parameters at the hidden layerdef gradient_weight_hidden(X, Eh):    return X.T.dot(Eh)# Define the gradient function for the bias parameters at the output layerdef gradient_bias_hidden(Eh):     return  np.sum(Eh, axis=0, keepdims=True)

梯度检查

在编程计算反向传播梯度时，很容易产生错误。这就是为什么一直推荐在你的模型中一定要进行梯度检查。梯度检查是通过对于每一个参数进行梯度数值计算进行的，即检查这个数值与通过反向传播的梯度进行比较计算。

假设对于参数θi，我们计算它的数值梯度∂ξ/∂θi如下：

其中，f是一个神经网络的方程，X是输入数据，θ表示所有参数的集合。ϵ是一个很小的值，用来对参数θi进行评估。

对于每个参数的数值梯度应该接近于反向传播梯度的参数。

# Initialize weights and biasesinit_var = 1# Initialize hidden layer parametersbh = np.random.randn(1, 3) * init_varWh = np.random.randn(2, 3) * init_var# Initialize output layer parametersbo = np.random.randn(1, 2) * init_varWo = np.random.randn(3, 2) * init_var# Compute the gradients by backpropagation# Compute the activations of the layersH = hidden_activations(X, Wh, bh)Y = output_activations(H, Wo, bo)# Compute the gradients of the output layerEo = error_output(Y, T)JWo = gradient_weight_out(H, Eo)Jbo = gradient_bias_out(Eo)# Compute the gradients of the hidden layerEh = error_hidden(H, Wo, Eo)JWh = gradient_weight_hidden(X, Eh)Jbh = gradient_bias_hidden(Eh)

# Combine all parameter matrices in a listparams = [Wh, bh, Wo, bo]# Combine all parameter gradients in a listgrad_params = [JWh, Jbh, JWo, Jbo]# Set the small change to compute the numerical gradienteps = 0.0001# Check each parameter matrixfor p_idx in range(len(params)):    # Check each parameter in each parameter matrix    for row in range(params[p_idx].shape[0]):        for col in range(params[p_idx].shape[1]):            # Copy the parameter matrix and change the current parameter slightly            p_matrix_min = params[p_idx].copy()            p_matrix_min[row,col] -= eps            p_matrix_plus = params[p_idx].copy()            p_matrix_plus[row,col] += eps            # Copy the parameter list, and change the updated parameter matrix            params_min = params[:]            params_min[p_idx] = p_matrix_min            params_plus = params[:]            params_plus[p_idx] =  p_matrix_plus            # Compute the numerical gradient            grad_num = (cost(nn(X, *params_plus), T)-cost(nn(X, *params_min), T))/(2*eps)            # Raise error if the numerical grade is not close to the backprop gradient            if not np.isclose(grad_num, grad_params[p_idx][row,col]):                raise ValueError('Numerical gradient of {:.6f} is not close to the backpropagation gradient of {:.6f}!'.format(float(grad_num), float(grad_params[p_idx][row,col])))print('No gradient errors found')

No gradient errors found

动量反向传播的更新

在前面几个例子中，我们使用最简单的梯度下降算法，根据损失函数来优化参数，因为这些损失函数都是凸函数。但是在多层神经网络，参数量非常巨大并且激活函数是非线性函数时，我们的损失函数极不可能是一个凸函数。那么此时，简单的梯度下降算法不是最好的方法去找到一个全局的最小值，因为这个方法是一个局部优化的方法，最后将收敛在一个局部最小值。

为了解决这个例子中的问题，我们使用一种梯度下降算法的改进版，叫做动量方法。这种动量方法，你可以想象成一只球从损失函数的表面从高处落下。在下降的过程中，这个球的速度会增加，但是当它上坡时，它的速度会下降，这一个过程可以用下面的数学公式进行描述：

其中，V(i)表示参数第i次迭代时的速度。θ(i)表示参数在第i次迭代时的值。∂ξ/∂θ(i)表示参数在第i次迭代时的梯度。λ表示速度根据阻力减小的值，μ表示学习率。这个数学描述公式可以被可视化为如下图：

速度VWh，Vbh，VWo和Vbo对应于参数Wh，bh，Wo和bo，在代码中，这些参数用VWh，Vbh，VWo和Vbo表示，并且被存储在一个列表Vs中。update_velocity(X, T, ls_of_params, Vs, momentum_term, learning_rate)函数实现了上面的方法。update_params(ls_of_params, Vs)函数实现了速度的更新方法。

# Define the update function to update the network parameters over 1 iterationdef backprop_gradients(X, T, Wh, bh, Wo, bo):    # Compute the output of the network    # Compute the activations of the layers    H = hidden_activations(X, Wh, bh)    Y = output_activations(H, Wo, bo)    # Compute the gradients of the output layer    Eo = error_output(Y, T)    JWo = gradient_weight_out(H, Eo)    Jbo = gradient_bias_out(Eo)    # Compute the gradients of the hidden layer    Eh = error_hidden(H, Wo, Eo)    JWh = gradient_weight_hidden(X, Eh)    Jbh = gradient_bias_hidden(Eh)    return [JWh, Jbh, JWo, Jbo]def update_velocity(X, T, ls_of_params, Vs, momentum_term, learning_rate):    # ls_of_params = [Wh, bh, Wo, bo]    # Js = [JWh, Jbh, JWo, Jbo]    Js = backprop_gradients(X, T, *ls_of_params)    return [momentum_term * V - learning_rate * J for V,J in zip(Vs, Js)]def update_params(ls_of_params, Vs):    # ls_of_params = [Wh, bh, Wo, bo]    # Vs = [VWh, Vbh, VWo, Vbo]    return [P + V for P,V in zip(ls_of_params, Vs)]

# Run backpropagation# Initialize weights and biasesinit_var = 0.1# Initialize hidden layer parametersbh = np.random.randn(1, 3) * init_varWh = np.random.randn(2, 3) * init_var# Initialize output layer parametersbo = np.random.randn(1, 2) * init_varWo = np.random.randn(3, 2) * init_var# Parameters are already initilized randomly with the gradient checking# Set the learning ratelearning_rate = 0.02momentum_term = 0.9# define the velocities Vs = [VWh, Vbh, VWo, Vbo]Vs = [np.zeros_like(M) for M in [Wh, bh, Wo, bo]]# Start the gradient descent updates and plot the iterationsnb_of_iterations = 300  # number of gradient descent updateslr_update = learning_rate / nb_of_iterations # learning rate update rulels_costs = [cost(nn(X, Wh, bh, Wo, bo), T)]  # list of cost over the iterationsfor i in range(nb_of_iterations):    # Update the velocities and the parameters    Vs = update_velocity(X, T, [Wh, bh, Wo, bo], Vs, momentum_term, learning_rate)    Wh, bh, Wo, bo = update_params([Wh, bh, Wo, bo], Vs)    ls_costs.append(cost(nn(X, Wh, bh, Wo, bo), T))

# Plot the cost over the iterationsplt.plot(ls_costs, 'b-')plt.xlabel('iteration')plt.ylabel('$\\xi$', fontsize=15)plt.title('Decrease of cost over backprop iteration')plt.grid()plt.show()

Decrease of cost over backprop iteration

可视化训练分类结果

在下图中，我们利用输入数据X和目标结果T，利用动量方法对BP算法进行更新得到的分类边界。我们利用红色和蓝色去区分输入数据的分类域。从图中，我们可以看出，我们把所有的样本都正确分类了。

# Plot the resulting decision boundary# Generate a grid over the input space to plot the color of the#  classification at that grid pointnb_of_xs = 200xs1 = np.linspace(-2, 2, num=nb_of_xs)xs2 = np.linspace(-2, 2, num=nb_of_xs)xx, yy = np.meshgrid(xs1, xs2) # create the grid# Initialize and fill the classification planeclassification_plane = np.zeros((nb_of_xs, nb_of_xs))for i in range(nb_of_xs):    for j in range(nb_of_xs):        pred = nn_predict(np.asmatrix([xx[i,j], yy[i,j]]), Wh, bh, Wo, bo)        classification_plane[i,j] = pred[0,0]# Create a color map to show the classification colors of each grid pointcmap = ListedColormap([        colorConverter.to_rgba('b', alpha=0.30),        colorConverter.to_rgba('r', alpha=0.30)])# Plot the classification plane with decision boundary and input samplesplt.contourf(xx, yy, classification_plane, cmap=cmap)# Plot both classes on the x1, x2 planeplt.plot(x_red[:,0], x_red[:,1], 'ro', label='class red')plt.plot(x_blue[:,0], x_blue[:,1], 'bo', label='class blue')plt.grid()plt.legend(loc=1)plt.xlabel('$x_1$', fontsize=15)plt.ylabel('$x_2$', fontsize=15)plt.axis([-1.5, 1.5, -1.5, 1.5])plt.title('red vs blue classification boundary')plt.show()

red vs blue classification boundary

输入域的转换

存在两个理由去解释为什么这个神经网络可以将这个非线性的数据进行分类。第一，因为隐藏层中使用了非线性的Logistic函数来帮助数据分类。第二，隐藏层使用了三个维度，即我们将数据从低维映射到了高维数据。下面的图绘制除了在隐藏层中的三维数据分类。

# Plot the projection of the input onto the hidden layer# Define the projections of the blue and red classesH_blue = hidden_activations(x_blue, Wh, bh)H_red = hidden_activations(x_red, Wh, bh)# Plot the error surfacefig = plt.figure()ax = Axes3D(fig)ax.plot(np.ravel(H_blue[:,0]), np.ravel(H_blue[:,1]), np.ravel(H_blue[:,2]), 'bo')ax.plot(np.ravel(H_red[:,0]), np.ravel(H_red[:,1]), np.ravel(H_red[:,2]), 'ro')ax.set_xlabel('$h_1$', fontsize=15)ax.set_ylabel('$h_2$', fontsize=15)ax.set_zlabel('$h_3$', fontsize=15)ax.view_init(elev=10, azim=-40)plt.title('Projection of the input X onto the hidden layer H')plt.grid()plt.show()

Projection of the input X onto the hidden layer H