一：神經網絡實現識別手寫數字
使用神經網絡再次實現處理手寫數字數據集。通過反向傳播算法實現神經網絡成本函數和梯度計算得非正則化和正則化版本。還將實現隨機權重初始化和使用網絡進行預測得方法。
（一）導入庫，並且讀取數據集
因為我們的數據集類型是.mat文件（是在matlab的本機格式），所以在使用python加載時，我們需要使用一個SciPy工具。
import numpy as np
import pandas as pd import matplotlib as plt from scipy.io import loadmat #用於載入matlab文件 data = loadmat('ex3data1.mat') #加載數據 print(data)
print(data['X'].shape)
print(data['y'].shape)
圖像在martix X中表示為400維向量（其中有5000個數據）。400維“特征”是原始圖像（20X20)圖像中每個像素的灰度強度。
類標簽在向量y中作為圖像中數字的數字類。
第一個任務是將我們的邏輯函數實現修改為完全向量化（即沒有for循環）。這是因為向量化代碼除了簡潔外，還能利用線性代數進行優化。
並且通常比迭代代碼快得多。
（二）部分數據可視化（5000中，隨機選取100個顯示）
def displayData(X,example_width=None): #顯示20x20像素數據在一個格子里面
if example_width is None: example_width = math.floor(math.sqrt(X.shape[1])) #X.shape[1]=400 print(X.shape) #計算行列 m,n = X.shape #返回的是我們隨機選取的100行樣本---（100，400） #獲取圖像高度 example_height = math.ceil(n/example_width) #400/20 = 20 # 計算需要展示的圖片數量 display_rows = math.floor(math.sqrt(m)) #10 display_cols = math.ceil(m / display_rows) #100/10=10 fig,ax = plt.subplots(nrows=display_rows,ncols=display_cols,sharey=True,sharex=True,figsize=(12,12)) #https://blog.csdn.net/qq_39622065/article/details/82909421 #拷貝圖像到子區域進行展示 for row in range(display_rows): for col in range(display_cols): ax[row,col].matshow(X[display_rows*row+col].reshape(example_width,example_height),cmap='gray') #希望將矩陣以圖形方式顯示https://www.cnblogs.com/shanger/articles/11872437.html plt.xticks([]) #用於設置每個小區域坐標軸刻度，設置為空，則不顯示刻度 plt.yticks([]) plt.show()
plt.rcParams['figure.dpi'] = 200
data = loadmat('ex3data1.mat') #加載數據 # print(data['X'].shape) # print(data['y'].shape) X = data.get('X') y = data.get('y') sample_idx = np.random.choice(np.arange(X.shape[0]),100) #隨機選取100個數（100個樣本的序號） # print(sample_idx) samples = X[sample_idx,:] #獲取上面選取的索引對應的數據集 smaples_y = y[sample_idx,:] #可以知道我們選取的數字標簽值 # print(smaples_y) displayData(samples) #進行展示
（三）對y標簽值進行one-hot編碼
舉例（1）：one-hot編碼
from sklearn.preprocessing import OneHotEncoder import numpy as np
a = np.array([[1,2,3,4,5，5，3，0,6,7,8,9,0]]) b = a.T
c = encoder.fit_transform(b)
可以看出下列對應關系：
舉例（2）：one-hot編碼
總之：one-hot編碼，會根據我們標簽值中的分類數（實驗中是10個），生成相應的維度。每個標簽值對應唯一一個行向量（該向量中只有一個1，其他全為0）
（詳細案例可以見--五：利用神經網絡解決多類別分類問題）
將輸出得多分類標簽值y（這里共10類），進行編碼，都轉換為一個10維列向量。
from sklearn.preprocessing import OneHotEncoder #預處理操作 encoder = OneHotEncoder(sparse=False) #sparse=True 表示編碼的格式，默認為 True，即為稀疏的格式，指定 False 則就不用 toarray() 了 y_onehot = encoder.fit_transform(y)
print(y.shape)
print(y_onehot.shape)
print(y[0]) #簡單查看編碼內容 print(y_onehot[0,:])
（四）神經網絡結構圖
我們要為實驗構建的神經網絡：
輸入層：包括400個特征（20X20），加一個偏置單元。 隱藏層：25個單位的隱藏層（詳細見下面），加一個偏置單元 輸出層：根據我們的分類，設置10個輸出單元
這里，我們已經提供了一系列的神經網絡參數：θ(1)和θ(2)---即全部的權重參數（已經提前為我們訓練過了，存儲在ex4weights.mat）。這些參數是取決於我們的第二層隱藏層單元數（由於ex4weights.mat文件提供了全部參數，所以同時也要配置隱藏層單元數---25）---這里是從文件中獲取θ的初始值，然后進行迭代下降（我們會在后面提到隨機初始化方法，用來放置參數對稱問題）
theta_data = loadmat('ex4weights.mat') #加載數據 theta_1 = theta_data.get('Theta1') theta_2 = theta_data.get('Theta2') print(theta_1.shape) print(theta_2.shape)
（五）實現sigmoid函數
def sigmoid(z): #實現sigmoid函數 return 1/(1+np.exp(-z))
（六）實現前向傳播函數---計算各個層（各個節點）的激活值
各層單元數變化如下：
神經網絡結構：
def forward_propagate(X, theta1, theta2): #X是我們輸入層訓練集，theta1是我們第一層的參數，theta2是我們的第二層參數 #首先，我們要為輸入數據集X，添加偏置單元 a1 = np.insert(X,0,1,axis=1) #插入一列全為1的列向量到X中 #計算第二層激活值 z2 = a1@theta1.T #theta1是一個(25，401)矩陣，a1是（5000，401）---z2是（5000，25） a2 = np.insert(sigmoid(z2),0,1,axis=1) #.之后添加偏置單元。變為（5000，26） #計算第三層激活值 z3 = a2@theta2.T #theta2是一個（10，26）矩陣，a2是（5000，26），所以z3是一個（5000，10）矩陣 h = sigmoid(z3) #獲取假設函數值 return a1,z2,a2,z3,h
theta_data = loadmat('ex4weights.mat') #加載數據 theta_1 = theta_data.get('Theta1') theta_2 = theta_data.get('Theta2') a1,z2,a2,z3,h = forward_propagate(X,theta_1,theta_2)
（七）代價函數---最好保持同下面梯度函數參數一致
def cost(theta_param,input_size,hidden_size,num_labels,X,y,lamda=1): #lamda是正則化參數 #獲取樣本個數 m = X.shape[0] theta_1 = theta_param[:(input_size + 1) * hidden_size].reshape(hidden_size, (input_size + 1)) theta_2 = theta_param[(input_size + 1) * hidden_size:].reshape(num_labels, (hidden_size + 1)) #將獲得神經元單元激活值等信息 a1,z2,a2,z3,h = forward_propagate(X,theta_1,theta_2) #初始化代價函數 J = 0 #根據公式計算代價函數 for i in range(m): #遍歷每一個樣本---m first_term = np.multiply(-y[i,:],np.log(h[i,:])) #標簽向量，乘以假設函數矩陣---K second_term = np.multiply((1-y[i,:]),np.log(1-h[i,:])) J += np.sum(first_term-second_term) J = J / m #加上正則化項 J += (lamda/(2*m))*(np.power(theta_1[:,1:],2).sum()+np.power(theta_2[:,1:],2).sum()) #不含θ_0項 return J
data = loadmat('ex4data1.mat') #加載數據 theta_data = loadmat('ex4weights.mat') # 加載數據 # 初始化操作 input_size = 400 #輸入層 hidden_size = 25 #隱藏層 num_labels = 10 #輸出層 theta_1 = theta_data.get('Theta1') theta_2 = theta_data.get('Theta2') theta_param = np.concatenate([np.ravel(theta_1),np.ravel(theta_2)])
1.當我們設置lamda為0，不含正則化時：
print(cost(theta_param,input_size,hidden_size,num_labels,X,y_onehot,0))
2.當設置lamda為1時：
print(cost(theta_param,input_size,hidden_size,num_labels,X,y_onehot,1))
（八）反向傳播 了解
1.輸入層不需要計算誤差
2.輸出層不需要計算sigmoid函數偏導：
3.其它層，我們需要以下計算：
（九）sigmoid函數求導
def sigmoid_gradient(z): return sigmoid(z)*(1-sigmoid(z))
（十）反向傳播算法---計算誤差、偏導
def backporp(theta_param,input_size,hidden_size,num_labels,X,y,lamda=1): # 獲取樣本個數 m = X.shape[0] theta_1 = theta_param[:(input_size + 1) * hidden_size].reshape(hidden_size, (input_size + 1)) theta_2 = theta_param[(input_size + 1) * hidden_size:].reshape(num_labels, (hidden_size + 1)) delta1 = np.zeros(theta_1.shape) #由代價函數公式可以知道delta和theta是同型矩陣 delta2 = np.zeros(theta_2.shape) # 先前向傳播，將獲得神經元單元激活值等信息 a1, z2, a2, z3, h = forward_propagate(X, theta_1, theta_2) # 由於theta_2是（10，26）矩陣，所以d3i*theta_2是（1，26）矩陣，為了可以相乘，則我們要將z2加上一列 z2 = np.insert(z2, 0, 1, axis=1) #求代價函數 J = cost(theta_param,input_size,hidden_size,num_labels,X,y,lamda) #實現反向傳播 for i in range(m): a1i = a1[i,:] #獲取a1中得第i個 z2i = z2[i,:] #是（1，25）矩陣 a2i = a2[i,:] hi = h[i,:] yi = y[i,:] #進行誤差計算 #先獲取輸出層3，的誤差δ d3i = hi - yi #是（1，10）矩陣 d3i = np.array([d3i]) #再獲取隱藏層誤差 d2i = np.multiply((d3i@theta_2),np.array([sigmoid_gradient(z2i)])) #(1,26) #獲取各層誤差 delta1 = delta1 + (d2i[:,1:]).T*a1i #獲取第一層誤差 (25,401) delta2 = delta2 + (d3i[:,:]).T*a2i #獲取第二層誤差 (10,26) delta1 = delta1 / m delta2 = delta2 / m #加入正則化操作 注意：我們這里不包含δ_0 delta1[:,1:] = delta1[:,1:] + (theta_1[:,1:]*lamda) / m delta2[:,1:] = delta2[:,1:] + (theta_2[:,1:]*lamda) / m #將梯度矩陣轉化為單個數組 grad = np.concatenate([np.ravel(delta1),np.ravel(delta2)]) return J,grad
J,grad = backporp(theta_param,X,y_onehot,1)
（十一） 預測函數
def predict_new(theta1,theta2,X): X = np.insert(X,0,1,axis=1) #插入一列全為1的列向量到X中 h_1 = sigmoid(X@theta1.T) h_1 = np.insert(h_1,0,1,axis=1) h_2 = sigmoid(h_1@theta2.T) p = np.argmax(h_2,axis=1)+1 return p
#保證輸入theta_param時一個向量，而不是矩陣。並且返回的梯度向量同theta_param同緯度---重點 fmin = minimize(fun=backporp,x0=theta_param,args=(input_size,hidden_size,num_labels,X,y_onehot,1),method='TNC',jac=True,options={'maxiter':250}) theta_param_new = fmin.x theta_1_new = theta_param_new[:(input_size + 1) * hidden_size].reshape(hidden_size, (input_size + 1)) theta_2_new = theta_param_new[(input_size + 1) * hidden_size:].reshape(num_labels, (hidden_size + 1)) y_pred = predict_new(theta_1_new,theta_2_new,X) 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 ))
（十二）全部代碼
import math import numpy as np import pandas as pd import matplotlib.pyplot as plt from scipy.io import loadmat #用於載入matlab文件 import scipy.optimize as opt from scipy.optimize import minimize from sklearn.preprocessing import OneHotEncoder #預處理操作 def displayData(X,example_width=None): #顯示20x20像素數據在一個格子里面 if example_width is None: example_width = math.floor(math.sqrt(X.shape[1])) #X.shape[1]=400 print(X.shape) #計算行列 m,n = X.shape #返回的是我們隨機選取的100行樣本---（100，400） #獲取圖像高度 example_height = math.ceil(n/example_width) #400/20 = 20 # 計算需要展示的圖片數量 display_rows = math.floor(math.sqrt(m)) #10 display_cols = math.ceil(m / display_rows) #100/10=10 fig,ax = plt.subplots(nrows=display_rows,ncols=display_cols,sharey=True,sharex=True,figsize=(12,12)) #https://blog.csdn.net/qq_39622065/article/details/82909421 #拷貝圖像到子區域進行展示 for row in range(display_rows): for col in range(display_cols): ax[row,col].matshow(X[display_rows*row+col].reshape(example_width,example_height),cmap='gray') #希望將矩陣以圖形方式顯示https://www.cnblogs.com/shanger/articles/11872437.html plt.xticks([]) #用於設置每個小區域坐標軸刻度，設置為空，則不顯示刻度 plt.yticks([]) plt.show() def sigmoid(z): # 實現sigmoid函數 return 1 / (1 + np.exp(-z)) def one_to_all(X,y,num_labels,learning_rate): #分類器構建 X訓練集 y標簽值 num_labels是構建的分類器數量（0-9標簽值，則10個），learning_rate是入值 rows = X.shape[0] cols = X.shape[1] #構建多維度θ值 all_theta = np.zeros((num_labels,cols+1)) #這里我們需要添加θ_0,相應的我們要在X訓練集中加入一列全為1的列向量 #為X訓練集添加一個全為1的列向量 X = np.insert(X,0,values=np.ones(rows),axis=1) #第二個參數是插入索引位置，第三個是插入數據，axis=1是列插入https://www.cnblogs.com/hellcat/p/7422005.html#_label2_9 #開始進行參數擬合 for i in range(1,num_labels+1): #從1-10代替0-9，因為數據集的標簽值是從1-10，print(np.unique(data['y']))可以知道 theta = np.zeros(cols+1) #初始化一個新的臨時的θ值 y_i = np.array([1 if label == i else 0 for label in y]) #進行二分類，對於對應的數設置為1，其他為0 y_i = np.reshape(y_i,(rows,1)) #進行轉置為列向量 #使用高級優化算法進行參數迭代 fmin = opt.minimize(fun=regularized_cost,x0=theta,args=(X, y_i), method='TNC',jac=regularized_gradient) all_theta[i-1,:]=fmin.x #將1-10變為0-9 return all_theta #返回多維度θ值 def forward_propagate(X, theta1, theta2): #X是我們輸入層訓練集，theta1是我們第一層的參數，theta2是我們的第二層參數 #首先，我們要為輸入數據集X，添加偏置單元 a1 = np.insert(X,0,1,axis=1) #插入一列全為1的列向量到X中 #計算第二層激活值 z2 = a1@theta1.T #theta1是一個(25，401)矩陣，a1是（5000，401）---z2是（5000，25） a2 = np.insert(sigmoid(z2),0,1,axis=1) #.之后添加偏置單元。變為（5000，26） #計算第三層激活值 z3 = a2@theta2.T #theta2是一個（10，26）矩陣，a2是（5000，26），所以z3是一個（5000，10）矩陣 h = sigmoid(z3) #獲取假設函數值 return a1,z2,a2,z3,h def cost(theta_param,X,y,lamda=1): #lamda是正則化參數 #獲取樣本個數 m = X.shape[0] theta_1 = theta_param[:(input_size + 1) * hidden_size].reshape(hidden_size, (input_size + 1)) theta_2 = theta_param[(input_size + 1) * hidden_size:].reshape(num_labels, (hidden_size + 1)) #將獲得神經元單元激活值等信息 a1,z2,a2,z3,h = forward_propagate(X,theta_1,theta_2) #初始化代價函數 J = 0 #根據公式計算代價函數 for i in range(m): #遍歷每一個樣本---m first_term = np.multiply(-y[i,:],np.log(h[i,:])) #標簽向量，乘以假設函數矩陣---K second_term = np.multiply((1-y[i,:]),np.log(1-h[i,:])) J += np.sum(first_term-second_term) J = J / m #加上正則化項 J += (lamda/(2*m))*(np.power(theta_1[:,1:],2).sum()+np.power(theta_2[:,1:],2).sum()) #不含θ_0項 return J def sigmoid_gradient(z): return sigmoid(z)*(1-sigmoid(z)) def backporp(theta_param,X,y,lamda=1): # 獲取樣本個數 m = X.shape[0] theta_1 = theta_param[:(input_size + 1) * hidden_size].reshape(hidden_size, (input_size + 1)) theta_2 = theta_param[(input_size + 1) * hidden_size:].reshape(num_labels, (hidden_size + 1)) delta1 = np.zeros(theta_1.shape) #由代價函數公式可以知道delta和theta是同型矩陣 delta2 = np.zeros(theta_2.shape) # 先前向傳播，將獲得神經元單元激活值等信息 a1, z2, a2, z3, h = forward_propagate(X, theta_1, theta_2) # 由於theta_2是（10，26）矩陣，所以d3i*theta_2是（1，26）矩陣，為了可以相乘，則我們要將z2加上一列 z2 = np.insert(z2, 0, 1, axis=1) #求代價函數 J = cost(theta_param,X,y,lamda) #實現反向傳播 for i in range(m): a1i = a1[i,:] #獲取a1中得第i個 z2i = z2[i,:] #是（1，25）矩陣 a2i = a2[i,:] hi = h[i,:] yi = y[i,:] #進行誤差計算 #先獲取輸出層3，的誤差δ d3i = hi - yi #是（1，10）矩陣 d3i = np.array([d3i]) #再獲取隱藏層誤差 d2i = np.multiply((d3i@theta_2),np.array([sigmoid_gradient(z2i)])) #(1,26) #獲取各層誤差 delta1 = delta1 + (d2i[:,1:]).T*a1i #獲取第一層誤差 (25,401) delta2 = delta2 + (d3i[:,:]).T*a2i #獲取第二層誤差 (10,26) delta1 = delta1 / m delta2 = delta2 / m #加入正則化操作 注意：我們這里不包含δ_0 delta1[:,1:] = delta1[:,1:] + (theta_1[:,1:]*lamda) / m delta2[:,1:] = delta2[:,1:] + (theta_2[:,1:]*lamda) / m #將梯度矩陣轉化為單個數組 grad = np.concatenate([np.ravel(delta1),np.ravel(delta2)]) return J,grad def predict_new(theta1,theta2,X): X = np.insert(X,0,1,axis=1) #插入一列全為1的列向量到X中 h_1 = sigmoid(X@theta1.T) h_1 = np.insert(h_1,0,1,axis=1) h_2 = sigmoid(h_1@theta2.T) p = np.argmax(h_2,axis=1)+1 return p data = loadmat('ex4data1.mat') #加載數據 theta_data = loadmat('ex4weights.mat') # 加載數據 # 初始化操作 input_size = 400 #輸入層 hidden_size = 25 #隱藏層 num_labels = 10 #輸出層 theta_1 = theta_data.get('Theta1') theta_2 = theta_data.get('Theta2') theta_param = np.concatenate([np.ravel(theta_1),np.ravel(theta_2)]) X = data.get('X') y = data.get('y') encoder = OneHotEncoder(sparse=False) #sparse=True 表示編碼的格式，默認為 True，即為稀疏的格式，指定 False 則就不用 toarray() 了 y_onehot = encoder.fit_transform(y) # J,grad = backporp(theta_param,X,y_onehot,1) #保證輸入theta_param時一個向量，而不是矩陣。並且返回的梯度向量同theta_param同緯度 fmin = minimize(fun=backporp,x0=theta_param,args=(X,y_onehot,1),method='TNC',jac=True,options={'maxiter':250}) theta_param_new = fmin.x theta_1_new = theta_param_new[:(input_size + 1) * hidden_size].reshape(hidden_size, (input_size + 1)) theta_2_new = theta_param_new[(input_size + 1) * hidden_size:].reshape(num_labels, (hidden_size + 1)) y_pred = predict_new(theta_1_new,theta_2_new,X) 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 ))