【吴恩达机器学习】练习4-神经网络

fishni

人工智能

发布于：2020年7月1日

对于这个练习，我们将再次处理手写数字数据集，这次使用反向传播的前馈神经网络。我们将通过反向传播算法实现神经网络代价函数和梯度计算的非正则化和正则化版本。我们还将实现随机权重初始化和使用网络进行预测的方法。

数据集加载

由于我们在练习3中使用的数据集是相同的，所以我们将重新使用代码来加载数据。

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.io import loadmat

1 2	data = loadmat('ex4data1.mat') data

{'__header__': b'MATLAB 5.0 MAT-file, Platform: GLNXA64, Created on: Sun Oct 16 13:09:09 2011',
 '__version__': '1.0',
 '__globals__': [],
 'X': array([[0., 0., 0., ..., 0., 0., 0.],
        [0., 0., 0., ..., 0., 0., 0.],
        [0., 0., 0., ..., 0., 0., 0.],
        ...,
        [0., 0., 0., ..., 0., 0., 0.],
        [0., 0., 0., ..., 0., 0., 0.],
        [0., 0., 0., ..., 0., 0., 0.]]),
 'y': array([[10],
        [10],
        [10],
        ...,
        [ 9],
        [ 9],
        [ 9]], dtype=uint8)}

由于我们以后需要这些（并将经常使用它们），我们先来创建一些有用的变量。

X = data['X']
y = data['y']

X.shape, y.shape#看下维度

((5000, 400), (5000, 1))

one-hot编码

one hot编

是将类别变量转换为机器学习算法易于利用的一种形式的过程。

One-Hot编码是分类变量作为二进制向量的表示。这首先要求将分类值映射到整数值。然后，每个整数值被表示为二进制向量，除了整数的索引之外，它都是零值，它被标记为1。
示例：

我们也需要对我们的y标签进行一次one-hot 编码。 one-hot 编码将类标签n（k类）转换为长度为k的向量，其中索引n为“hot”（1），而其余为0。 Scikitlearn有一个内置的实用程序，我们可以使用这个。

from sklearn.preprocessing import OneHotEncoder
# 初始化OneHotEncoder实例时，默认sparse参数为True，编码后返回的是一个稀疏矩阵的对象，
#如果要使用一般要调用toarray()方法转化成array对象。
# 若将sparse参数设置为False，则直接生成array对象，可直接使用。
encoder = OneHotEncoder(sparse=False,categories='auto')
y_onehot = encoder.fit_transform(y)
y_onehot.shape

(5000, 10)

1	y[0], y_onehot[0,:]

(array([10], dtype=uint8), array([0., 0., 0., 0., 0., 0., 0., 0., 0., 1.]))

我们要为此练习构建的神经网络具有与我们的实例数据（400 +偏置单元）大小匹配的输入层，25个单位的隐藏层（带有偏置单元的26个），以及一个输出层， 10个单位对应我们的一个one-hot编码类标签。有关网络架构的更多详细信息和图像，请参阅“练习”文件夹中的PDF。

我们需要实现的第一件是评估一组给定的网络参数的损失的代价函数。源函数在练习文本中（看起来很吓人）。以下是代价函数的代码。

sigmoid 函数

g 代表一个常用的逻辑函数（logistic function）为S形函数（Sigmoid function），公式为： $g\left( z \right)=\frac{1}{1+{{e}^{-z}}}$
合起来，我们得到逻辑回归模型的假设函数：
${{h}_{\theta }}\left( x \right)=\frac{1}{1+{{e}^{-{{\theta }^{T}}X}}}$

1 2	def sigmoid(z): return 1 / (1 + np.exp(-z))

前向传播函数

(400 + 1) -> (25 + 1) -> (10)

1 2	m = X.shape[0] np.insert(X,0,values=np.ones(m),axis=1)

array([[1., 0., 0., ..., 0., 0., 0.],
       [1., 0., 0., ..., 0., 0., 0.],
       [1., 0., 0., ..., 0., 0., 0.],
       ...,
       [1., 0., 0., ..., 0., 0., 0.],
       [1., 0., 0., ..., 0., 0., 0.],
       [1., 0., 0., ..., 0., 0., 0.]])

def forward_propagate(X, theta1, theta2):
    m = X.shape[0] # 行数
    
    a1 = np.insert(X, 0, values=np.ones(m), axis=1) # 0插入位置，np.ones(m) m维度，axis=1列的方向插入
    z2 = a1 * theta1.T
    a2 = np.insert(sigmoid(z2), 0, values=np.ones(m), axis=1)
    z3 = a2 * theta2.T
    h = sigmoid(z3)
    
    return a1, z2, a2, z3, h

代价函数

def cost(params, input_size, hidden_size, num_labels, X, y, learning_rate):
    # 传入：神经网络参数，输入层维度，隐含层维度，标签数，训练数据，学习率
    m = X.shape[0] # 样本个数
    X = np.matrix(X) # 转换为numpy型矩阵
    y = np.matrix(y)
    
    # reshape the parameter array into parameter matrices for each layer
    # 为每一层重塑参数数组为参数矩阵
    #从params中获取神经网络参数，并按照输入层维度和隐藏层维度重新定义参数的维度
    theta1 = np.matrix(np.reshape(params[:hidden_size * (input_size + 1)], (hidden_size, (input_size + 1))))
    theta2 = np.matrix(np.reshape(params[hidden_size * (input_size + 1):], (num_labels, (hidden_size + 1))))
    
    # run the feed-forward pass运行前向传播函数
    # 调用前面写好的前项传播函数
    a1, z2, a2, z3, h = forward_propagate(X, theta1, theta2)
    
    # compute the cost计算代价
    J = 0
    for i in range(m): #遍历每个样本
        first_term = np.multiply(-y[i,:], np.log(h[i,:]))
        second_term = np.multiply((1 - y[i,:]), np.log(1 - h[i,:]))
        J += np.sum(first_term - second_term)
    
    J = J / m
    
    return J

这个Sigmoid函数我们以前使用过。前向传播函数计算给定当前参数的每个训练实例的假设。它的输出形状应该与y的一个one-hot编码相同。

# 初始化设置
input_size = 400
hidden_size = 25
num_labels = 10
learning_rate = 1

# 随机初始化完整网络参数大小的参数数组
params = (np.random.random(size=hidden_size * (input_size + 1) + num_labels * (hidden_size + 1)) - 0.5) * 0.25

m = X.shape[0]
X = np.matrix(X)
y = np.matrix(y)

# 将参数数组解开为每个层的参数矩阵（在不更改数组数据的情况下为数组提供新形状）
theta1 = np.matrix(np.reshape(params[:hidden_size * (input_size + 1)], (hidden_size, (input_size + 1))))
theta2 = np.matrix(np.reshape(params[hidden_size * (input_size + 1):], (num_labels, (hidden_size + 1))))

theta1.shape, theta2.shape

((25, 401), (10, 26))

1 2	a1, z2, a2, z3, h = forward_propagate(X, theta1, theta2) a1.shape, z2.shape, a2.shape, z3.shape, h.shape

((5000, 401), (5000, 25), (5000, 26), (5000, 10), (5000, 10))

代价函数在计算假设矩阵h之后，应用代价函数来计算y和h之间的总误差。

1	cost(params, input_size, hidden_size, num_labels, X, y_onehot, learning_rate)

7.182002885834476

整合上述代码

def sigmoid(z):
    return 1/(1+np.exp(-z))
def forward_propagate(X,theta1,theta2):

    # a1 5000*401  theta1.T(25*401).T  = z1 5000*25 
    # a2 5000*26  theta2.T(10*26).T    =5000*10 
    m = X.shape[0]
    a1 = np.insert(X, 0, values=np.ones(m), axis=1)
    z2 = a1* theta1.T
    a2 = np.insert(sigmoid(z2),0,values=np.ones(m),axis=1)
    z3 = a2 * theta2.T
    h = sigmoid(z3)
    return a1,z2,a2,z3,h
def cost(params, input_size,hidden_size,num_lables,X,y,learning_rate):
    m = X.shape[0]
    X = np.matrix(X)
    y = np.matrix(y)
    
    theta1 = np.matrix(np.reshape(params[:hidden_size*(input_size+1)],(hidden_size,input_size+1)))
    theta2 = np.matrix(np.reshape(params[hidden_size*(input_size+1):],(num_lables,hidden_size+1)))

    a1,z2,a2,z3,h = forward_propagate(X,theta1,theta2)
    
    J=0
    for i in range(m):

        first_term = np.multiply(-y[i,:], np.log(h[i,:]))
        second_term = np.multiply((1 - y[i,:]), np.log(1 - h[i,:]))
        J += np.sum(first_term - second_term)

    J = J / m
    return J
input_size= 400
hidden_size= 25
num_labels= 10
learning_rate =1

params = (np.random.random(size=hidden_size * (input_size + 1) + num_labels * (hidden_size + 1)) - 0.5) * 0.25

cost(params,input_size,hidden_size,num_labels,X,y_onehot,learning_rate)

7.132323728323127

正则化代价函数

我们的下一步是增加代价函数的正则化。它实际上并不像看起来那么复杂 - 事实上，正则化术语只是我们已经计算出的代价的一个补充。下面是修改后的代价函数。

def cost(params, input_size, hidden_size, num_labels, X, y, learning_rate):
    m = X.shape[0]
    X = np.matrix(X)
    y = np.matrix(y)
    
    # reshape the parameter array into parameter matrices for each layer
    theta1 = np.matrix(np.reshape(params[:hidden_size * (input_size + 1)], (hidden_size, (input_size + 1))))
    theta2 = np.matrix(np.reshape(params[hidden_size * (input_size + 1):], (num_labels, (hidden_size + 1))))
    
    # run the feed-forward pass
    a1, z2, a2, z3, h = forward_propagate(X, theta1, theta2)
    
    # compute the cost
    J = 0
    for i in range(m):
        first_term = np.multiply(-y[i,:], np.log(h[i,:]))
        second_term = np.multiply((1 - y[i,:]), np.log(1 - h[i,:]))
        J += np.sum(first_term - second_term)
    
    J = J/m
    
    # add the cost regularization term
    J += (float(learning_rate) / (2 * m)) * (np.sum(np.power(theta1[:,1:], 2)) + np.sum(np.power(theta2[:,1:], 2)))
    
    return J

1	cost(params, input_size, hidden_size, num_labels, X, y_onehot, learning_rate)

7.137673923569304

反向传播

接下来是反向传播算法。反向传播参数更新计算将减少训练数据上的神经网络误差。我们需要的第一件事是计算我们之前创建的Sigmoid函数的梯度的函数。

1 2	def sigmoid_gradient(z): return np.multiply(sigmoid(z), (1 - sigmoid(z)))

现在我们准备好实施反向传播来计算梯度。由于反向传播所需的计算是代价函数中所需的计算过程，我们实际上将扩展代价函数以执行反向传播并返回代价和梯度。

def backprop(params, input_size, hidden_size, num_labels, X, y, learning_rate):
    # 步骤一：获取样本个数
    m = X.shape[0]
    
    # 步骤二：将矩阵X,y转换为numpy型矩阵
    X = np.matrix(X)
    y = np.matrix(y)
    
    # 步骤三：从params中获取神经网络参数，并按照输入层维度和隐含层维度重新定义参数的维度
    # reshape the parameter array into parameter matrices for each layer
    theta1 = np.matrix(np.reshape(params[:hidden_size * (input_size + 1)], (hidden_size, (input_size + 1))))
    theta2 = np.matrix(np.reshape(params[hidden_size * (input_size + 1):], (num_labels, (hidden_size + 1))))
    
    # 步骤四：调用前面写好的前向传播函数
    # run the feed-forward pass
    a1, z2, a2, z3, h = forward_propagate(X, theta1, theta2)
    
    # 步骤五：初始化
    # initializations
    
    J = 0
    delta1 = np.zeros(theta1.shape)  # (25, 401)
    delta2 = np.zeros(theta2.shape)  # (10, 26)
    
    # 步骤六：计算代价函数
    # compute the cost
    for i in range(m):
        first_term = np.multiply(-y[i,:], np.log(h[i,:]))
        second_term = np.multiply((1 - y[i,:]), np.log(1 - h[i,:]))
        J += np.sum(first_term - second_term)
    
    J = J / m
    
    # 步骤七：正则化
    # add the cost regularization term
    J += (float(learning_rate) / (2 * m)) * (np.sum(np.power(theta1[:,1:], 2)) + np.sum(np.power(theta2[:,1:], 2)))
    
    # 步骤八：实现反向传播
    # perform backpropagation执行反向传播
    
    # a1 （5000,401），theta1 （25,401）,theta2 (10,26)
    # z2= （5000,401）*（25，401）.T=(5000,25) , a2 = (5000,26)
    # z3 =(5000,26)*(10,26).T= (5000,10)   a3 = h = (5000,10)
    for t in range(m):
        a1t = a1[t,:]  # (1, 401)
        z2t = z2[t,:]  # (1, 25)
        a2t = a2[t,:]  # (1, 26)
        ht = h[t,:]  # (1, 10)
        yt = y[t,:]  # (1, 10)
        
        
        d3t = ht - yt  # (1, 10)
        # theta2 （10,26）
        z2t = np.insert(z2t, 0, values=np.ones(1))  # (1, 26)
        # sigmoid_gradient(z2t) （1,26）
        # (theta2.T (26,10) d3t（1,10）.T ).T= (26,1).T = (1,26)
        d2t = np.multiply((theta2.T * d3t.T).T, sigmoid_gradient(z2t))  # (1, 26)
        
        delta1 = delta1 + (d2t[:,1:]).T * a1t  # 输入层到隐含层的误差
        delta2 = delta2 + d3t.T * a2t     # 隐含层到输出层的误差
        
    delta1 = delta1 / m
    delta2 = delta2 / m
    
    # 步骤九：将梯度矩阵转换为单个数组
    # unravel the gradient matrices into a single array
    grad = np.concatenate((np.ravel(delta1), np.ravel(delta2)))
    
    return J, grad

梯度矩阵转换为单个数组

a = np.array([[1, 2], [3, 4]])
b = np.array([[5, 6]])
print("111:",np.ravel(a))
print("222:",np.concatenate((a,b))) # 默认axis=0
print("ddd:",np.concatenate((a, b.T), axis=1))
print("--------")
np.concatenate((np.ravel(a),np.ravel(b)))

111: [1 2 3 4]
222: [[1 2]
 [3 4]
 [5 6]]
ddd: [[1 2 5]
 [3 4 6]]
--------





array([1, 2, 3, 4, 5, 6])

反向传播计算的最难的部分（除了理解为什么我们正在做所有这些计算）是获得正确矩阵维度。顺便说一下，你容易混淆了A * B与np.multiply（A，B）使用。基本上前者是矩阵乘法，后者是元素乘法（除非A或B是标量值，在这种情况下没关系）。
无论如何，让我们测试一下，以确保函数返回我们期望的。

1 2	J, grad = backprop(params, input_size, hidden_size, num_labels, X, y_onehot, learning_rate) J, grad.shape

(7.137673923569304, (10285,))

矩阵乘法和元素乘法

矩阵乘法

a = np.matrix(np.array([[1,2,3,3]])) # numpy型矩阵
b = np.matrix(np.array([[2],[3],[4],[5]]))

print(a.shape,b.shape)
print(a*b)
print(np.dot(a,b))

(1, 4) (4, 1)
[[35]]
[[35]]

元素乘法示例

a = np.matrix(np.array([[1,2,3,3]]))
b = np.matrix(np.array([[2,4,6,8]]))

np.multiply(b,a)

matrix([[ 2,  8, 18, 24]])

反向传播正则化

我们还需要对反向传播函数进行一个修改，即将梯度计算加正则化。最后的正式版本如下。

def backprop(params, input_size, hidden_size, num_labels, X, y, learning_rate):
    m = X.shape[0]
    X = np.matrix(X)
    y = np.matrix(y)
    
    # reshape the parameter array into parameter matrices for each layer
    theta1 = np.matrix(np.reshape(params[:hidden_size * (input_size + 1)], (hidden_size, (input_size + 1))))
    theta2 = np.matrix(np.reshape(params[hidden_size * (input_size + 1):], (num_labels, (hidden_size + 1))))
    
    # run the feed-forward pass
    a1, z2, a2, z3, h = forward_propagate(X, theta1, theta2)
    
    # initializations
    J = 0
    delta1 = np.zeros(theta1.shape)  # (25, 401)
    delta2 = np.zeros(theta2.shape)  # (10, 26)
    
    # compute the cost
    for i in range(m):
        first_term = np.multiply(-y[i,:], np.log(h[i,:]))
        second_term = np.multiply((1 - y[i,:]), np.log(1 - h[i,:]))
        J += np.sum(first_term - second_term)
    
    J = J / m
    
    # add the cost regularization term
    J += (float(learning_rate) / (2 * m)) * (np.sum(np.power(theta1[:,1:], 2)) + np.sum(np.power(theta2[:,1:], 2)))
    
    # perform backpropagation
    for t in range(m):
        a1t = a1[t,:]  # (1, 401)
        z2t = z2[t,:]  # (1, 25)
        a2t = a2[t,:]  # (1, 26)
        ht = h[t,:]  # (1, 10)
        yt = y[t,:]  # (1, 10)
        
        d3t = ht - yt  # (1, 10)
        
        z2t = np.insert(z2t, 0, values=np.ones(1))  # (1, 26)
        d2t = np.multiply((theta2.T * d3t.T).T, sigmoid_gradient(z2t))  # (1, 26)
        
        delta1 = delta1 + (d2t[:,1:]).T * a1t
        delta2 = delta2 + d3t.T * a2t
        
    delta1 = delta1 / m
    delta2 = delta2 / m
    
    # add the gradient regularization term
    delta1[:,1:] = delta1[:,1:] + (theta1[:,1:] * learning_rate) / m
    delta2[:,1:] = delta2[:,1:] + (theta2[:,1:] * learning_rate) / m
    
    # unravel the gradient matrices into a single array
    grad = np.concatenate((np.ravel(delta1), np.ravel(delta2)))
    
    return J, grad

1 2	J, grad = backprop(params, input_size, hidden_size, num_labels, X, y_onehot, learning_rate) J, grad.shape

(7.137673923569304, (10285,))

我们终于准备好训练我们的网络，并用它进行预测。这与以往的具有多类逻辑回归的练习大致相似。

from scipy.optimize import minimize

# minimize the objective function
fmin = minimize(fun=backprop, x0=params, args=(input_size, hidden_size, num_labels, X, y_onehot, learning_rate), 
                method='TNC', jac=True, options={'maxiter': 250})
fmin

    fun: 0.36936848574990994
    jac: array([-1.94682911e-04,  2.79038111e-06, -8.05264111e-06, ...,
      -7.06122696e-04, -9.89911805e-04, -9.76751347e-04])
message: 'Max. number of function evaluations reached'
   nfev: 250
    nit: 16
 status: 3
success: False
      x: array([-0.13037322,  0.01395191, -0.04026321, ..., -1.50678966,
      -0.13882294, -0.24985413])

由于目标函数不太可能完全收敛，我们对迭代次数进行了限制。我们的总代价已经下降到0.5以下，这是算法正常工作的一个很好的指标。让我们使用它发现的参数，并通过网络转发，以获得一些预测。

让我们使用它找到的参数，并通过网络前向传播以获得预测。

X = np.matrix(X)
theta1 = np.matrix(np.reshape(fmin.x[:hidden_size * (input_size + 1)], (hidden_size, (input_size + 1))))
theta2 = np.matrix(np.reshape(fmin.x[hidden_size * (input_size + 1):], (num_labels, (hidden_size + 1))))

a1, z2, a2, z3, h = forward_propagate(X, theta1, theta2)
y_pred = np.array(np.argmax(h, axis=1) + 1)
y_pred

array([[10],
       [10],
       [10],
       ...,
       [ 9],
       [ 9],
       [ 9]], dtype=int64)

np.argmax测试

1
2
3

a=np.array([[10, 11, 12],
           [13, 14, 15]])
np.argmax(a, axis=1)

array([2, 2], dtype=int64)

最后，我们可以计算准确度，看看我们训练完毕的神经网络效果怎么样。

1
2
3

correct = [1 if a == b else 0 for (a, b) in zip(y_pred, y)]
accuracy = (sum(map(int, correct)) / float(len(correct)))
print ('accuracy = {0}%'.format(accuracy * 100))