Neural networks
Visualizing the data
在這一部分,首先需要加載數據並隨機輸出幾個圖像。
加載的數據有 \(15000\) 個訓練樣本(training examples),每一個訓練樣本是一個 \(64\times 64\) 像素的灰度圖。每一個像素代表了一個 \(8\) 位的無符號整數,代表了每個位置的灰度強度。這 \(64\times 64\) 的像素網格將被展開成為一個 \(4096\) 維的向量。這些訓練樣本將會成為矩陣 \(X\) 的每一行。所以維度為 \(15000\times 4096\) 的矩陣。
第二部分,是一個 \(15000\) 維的向量 \(y\) 其包含了所有訓練集的標簽。其值是從 \(1\) 到 \(15\) 分別表示漢字:零、一、二、三、四、五、六、七、八、九、十、百、千、萬、億。
並且將加載的數據集按照 \(7:3\) 的比例分成兩個集合:訓練集(training set)、測試集(test set)。
from matplotlib import pyplot as plt
from matplotlib import image as img
import numpy as np
from os import listdir
from scipy import optimize
import pandas as pd
import openpyxl
import time
from scipy import io
import re
n = 64 * 64 # 特征數量
m = 15000 # 數據樣本數量
input_layer_size = n # 64 * 64 的灰度圖像
hidden_layer_size = 256 # 256 個隱藏單元
num_labels = 15 # 15 個分類
Lambda = 1 # 正則化系數
def read_images(dir_str, n):
files_list = listdir(dir_str)
files_str = str(files_list)
# 樣本數量
m = len(files_list)
data = np.zeros((num_labels, m // num_labels, n))
for label in range(1, num_labels + 1):
pattern = 'input_\d*_\d*_' + str(label) + '\.jpg'
label_list = re.findall(pattern, files_str)
count = 0
for file in label_list:
data[label - 1, count, :] = np.ravel(img.imread(dir_str + file), order='F')
count += 1
return data
dir_str = './data/'
data = read_images(dir_str, n)
%matplotlib inline
fig, axes = plt.subplots(nrows=3, ncols=5) # 子圖為 3 行,5 列
i = 0
for axe in axes:
for ax in axe:
ax.imshow(data[i, 2, :].reshape(int(np.sqrt(n)), int(np.sqrt(n)), order='F'), cmap='gray')
i += 1
Xy = np.empty((int(m * 0.7), n + 1))
Xytest = np.empty((int(m * 0.3), n + 1))
for i in range(num_labels):
np.random.shuffle(data[i])
Xy[i * int(m // num_labels * 0.7) : i * int(m // num_labels * 0.7) + int(m // num_labels * 0.7), 0:-1] = data[i, 0:int(m // num_labels * 0.7), :]
Xy[i * int(m // num_labels * 0.7) : i * int(m // num_labels * 0.7) + int(m // num_labels * 0.7), -1] = np.ones(int(m // num_labels * 0.7)) * (i + 1)
Xytest[i * int(m // num_labels * 0.3) : i * int(m // num_labels * 0.3) + int(m // num_labels * 0.3), 0:-1] = data[i, int(m // num_labels * 0.7):, :]
Xytest[i * int(m // num_labels * 0.3) : i * int(m // num_labels * 0.3) + int(m // num_labels * 0.3), -1] = np.ones(int(m // num_labels * 0.3)) * (i + 1)
np.random.shuffle(Xy)
np.random.shuffle(Xytest)
X = Xy[:, 0:n]
y = Xy[:, -1]
Xtest = Xytest[:, 0:n]
ytest = Xytest[:, -1]
mat = io.loadmat('ex4data2.mat')
X = mat.get('X')
y = mat.get('y')
y = np.ravel(y)
Xtest = mat.get('Xtest')
ytest = mat.get('ytest')
ytest = np.ravel(ytest)
Model representation
神經網絡如下圖所示:
共有 \(3\) 層, \(1\) 個輸入層, \(1\) 個隱層, \(1\) 個輸出層。
每張圖片共有 \(64\times 64\) 個像素,所以輸入層有 \(4096\) 個輸入單元,每個隱層有 \(256\) 個單元,輸出層為 \(15\) 個輸出單元,表示 \(15\) 個類別的個數。
Feedforward and cost function
神經網絡的帶有正則化項的損失函數(cost function with regularization)表示為:
其中,\(h_\theta(x^{(i)})\) 表示向神經網絡輸入計算后的輸出, \(K=15\) 表示所有可能的標簽的個數, \(h_\theta(x^{(i)})_k=a_k^{(3)}\) 表示第 \(k\) 個輸出單元的激勵值(activation),並且不對偏置單元(bias unit)進行正則化。
除此之外,模型的原始標簽值是 \(1,2,\dots,15\) 的整數,為了訓練神經網絡,需要將這些十進制的整數值重新編碼為只有 \(0\) 和 \(1\) 的標簽向量:
假設 \(x^{(i)}\) 表示漢字“五”的圖片,那么對應的 \(y^{(i)}\) 應該是一個 \(15\) 維的向量,並且 \(y_5=1\) 其它元素都等於 \(0\) 。
def nnCostFun(nn_params, X, y, input_layer_size=input_layer_size, hidden_layer_size=hidden_layer_size, num_labels=num_labels, Lambda=Lambda):
Theta1 = nn_params[0:hidden_layer_size * (input_layer_size + 1)].reshape(hidden_layer_size, input_layer_size + 1, order='F')
Theta2 = nn_params[hidden_layer_size * (input_layer_size + 1) : ].reshape(num_labels, hidden_layer_size + 1, order='F')
m = X.shape[0] # 樣本的個數
# 將樣本標簽轉換為 15 維的 0-1 向量
yy = np.zeros((m, num_labels))
for i in range(m):
yy[i, int(y[i] - 1)] = 1
# 執行前向傳播
# 給 X 增加 1 列,全為 1
a1 = np.ones((m, 1), np.uint)
a1 = np.hstack((a1, X))
# 前向傳播
z2 = np.matmul(a1, Theta1.T)
a2 = sigmoid(z2)
# 增加第二層(隱藏層)的偏置單元
ones = np.ones((m, 1), dtype=np.int)
a2 = np.hstack((ones, a2))
z3 = np.matmul(a2, Theta2.T)
a3 = sigmoid(z3) # m * num_labels
# 計算 J
s = -yy * np.log(a3)
s -= (1 - yy) * np.log(1 - a3)
s = np.sum(s) # 計算每個輸出單元的累加
J = (1 / m) * s;
# 正則化 J
t1 = np.power(Theta1, 2)
t2 = np.power(Theta2, 2)
s = np.sum(t1[:, 1:]) + np.sum(t2[:, 1:])
J += (Lambda / (2 * m)) * s
# 反向傳播計算梯度
delta3 = a3 - yy
delta2 = (np.matmul(delta3, Theta2)[:, 1:] * sigmoid_gradient(z2))
Theta1_grad = (1 / m) * np.matmul(delta2.T, a1)
Theta2_grad = (1 / m) * np.matmul(delta3.T, a2)
# 正則化 D
Theta1_grad[:, 1:] += (Lambda / m) * Theta1[:, 1:]
Theta2_grad[:, 1:] += (Lambda / m) * Theta2[:, 1:]
gradient = [np.ravel(Theta1_grad, order='F').tolist(), np.ravel(Theta2_grad, order='F').tolist()]
gradient = gradient[0] + gradient[1]
return J, np.array(gradient)
def costFun(params):
return nnCostFun(params, X, y, input_layer_size, hidden_layer_size, num_labels, Lambda)
Backpropagation
在這一部分,我們實現后 BP 算法計算神經網絡的損失函數的梯度(gradient)。
Sigmoid gradient
首先,我們實現 Sigmoid 函數的梯度:
對於預測的越准確的激勵值, \(z\) 的值也就越大,梯度就越趨向於 \(0\) 。當 \(z=0\) 時,梯度的值應該正好是 \(0.25\) 。
def sigmoid(z):
return 1.0 / (1 + np.exp(-z))
def sigmoid_gradient(z):
return sigmoid(z) * (1 - sigmoid(z))
Random initialization
當開始訓練網絡之前,隨機初始化每一個參數用來對稱破缺(symmetry breaking)是非常重要的。一個隨機初始化非常有效的策略是,隨機地為 \(\Theta^{(l)}\) 選擇一個在統一 \(\left[-\epsilon_{init},\epsilon_{init}\right]\) 范圍內的值。我們可以使用 \(\epsilon_{init}=0.12\) 。這個范圍保證了參數保持在一個非常小的值,並使得學習更加有效。
def randInitialize(L_in, L_out):
epsilon = 0.12
W = np.random.rand(L_out, 1 + L_in) * 2 * epsilon - epsilon
return W
Backpropagation
給定一個訓練樣本 \((x^{(t)}, y^{(t)})\) ,首先輸入神經網絡並前向傳播計算神經網絡所有的激勵值,包括輸出值 \(h_\Theta(x)\) 。然后,對於第 \(l\) 層的結點 \(j\) ,將通過測量結點對於任何誤差的“負責度”(responsible),計算誤差項(error item) \(\delta_j^{(l)}\) 。
對於一個輸出結點,可以直接通過測量神經網絡的激勵值與真實的目標值的差值,並使用其定義 \(\delta_j^{(3)}\) 。對於隱藏單元,將基於第 \((l+1)\) 層的結點的誤差項的加權平均值去計算 \(\delta_j^{(l)}\) 。
- 對於第 \(t\) 個訓練樣本 \(x^{(t)}\) ,設置輸入層的值 \(a^{(1)}\) 。執行前向傳播,計算每層的激勵值 \((z^{(2)},a^{(2)},z^{(3)},a^{(3)})\) 。每一層都需要增加一個 \(+1\) 項的偏置單元。
- 對於第 \(3\) 層的輸出單元 \(k\) ,設
- 對於隱藏層 \(l=2\) ,設
- 使用如下公式積累神經網絡的梯度值:
注意,此步驟不對偏置單元 \(\delta_0^{(l)}\) 進行計算。
5. 將通過使用累積的梯度除以 \(m\) 得到神經網絡的損失函數的梯度:
Gradient checking
可以將參數矩陣 \(\Delta\) 展開成為一個非常長的向量 \(\theta\) 。通過這樣做,可以將損失函數看作 \(J(\theta)\) 並使用下面的步驟進行梯度檢查。
假設有一個函數 \(f_i(\theta)\) 去計算 \(\frac{\partial}{\partial\theta_i}J(\theta)\) 檢查 \(f_i\) 是否輸出了正確的導數值。設
因此, \(\theta^{(i+)}\) 除了第 \(i\) 個元素的值被增加了 \(\epsilon\) 其它值都與 \(\theta\) 的值相同,與之類似地, \(\theta^{(i-)}\) 也是除了第 \(i\) 個元素被減少了 \(\epsilon\) 其余的值都與 \(\theta\) 相同。可以使用數值驗證對於每一個 \(i\) 的 \(f_i(\theta)\) 的正確性:
這兩個值彼此的接近程度將取決於 \(J\) 的細節。假設 \(\epsilon=1e^{-4}\) ,通常將會發現上面左手和右手邊的式子至少接受 \(4\) 位有效數字。
def numericalGradient(theta, X, y, in_lay, hidden_lay, num, Lambda):
numgrad = np.zeros(theta.shape)
perturb = np.zeros(theta.shape)
e = 1e-4
for p in range(np.size(theta)):
# 設置擾動向量
perturb[p] = e
loss1 = nnCostFun(theta - perturb, X, y, in_lay, hidden_lay, num, Lambda)
loss1 = loss1[0]
loss2 = nnCostFun(theta + perturb, X, y, in_lay, hidden_lay, num, Lambda)
loss2 = loss2[0]
# 計算數值解梯度
numgrad[p] = (loss2 - loss1) / (2 * e)
perturb[p] = 0
return numgrad
def checkGradient(Lambda=0):
# 測試數據
in_lay = 3
hidden_lay = 5
num = 3
mm = 5
# 生成參數和樣本
Theta1 = np.arange(hidden_lay * (in_lay + 1)) + 1
Theta2 = np.arange(num * (hidden_lay + 1)) + 1
Theta1 = np.sin(Theta1).reshape(hidden_lay, in_lay + 1)
Theta2 = np.sin(Theta2).reshape(num, hidden_lay + 1)
X = np.arange(mm * in_lay) + 1
X = np.sin(X).reshape(mm, in_lay)
y = 1 + np.mod(np.arange(mm) + 1, num)
# 展開參數
pars = [np.ravel(Theta1).tolist(), np.ravel(Theta2).tolist()]
pars = pars[0] + pars[1]
pars = np.array(pars)
# 分別計算梯度
grad = nnCostFun(pars, X, y, in_lay, hidden_lay, num, Lambda)
grad = grad[1]
numgrad = numericalGradient(pars, X, y, in_lay, hidden_lay, num, Lambda)
# 計算差的范數
diff = np.linalg.norm(numgrad - grad) / np.linalg.norm(numgrad + grad)
return diff
checkGradient()
1.6291376831737205e-10
Regularized Neural Networks
增加一個正則化項到梯度上,在使用反向傳播計算完 \(\Delta_{ij}^{(l)}\) 之后,使用如下公式對梯度正則化:
注意不會對 \(\Theta^{(l)}\) 的第一列的偏置項進行正則化。
Learning parameters
在成功地實現了神經網絡的損失函數和梯度的計算之后,下一步走將會使用優化函數學習一個良好的參數集。
Minibatch Stochasitc Gradient Descent
按照數據生成分布抽取 \(m\) 個小批量(獨立同分布的)樣本,通過計算它們的梯度均值,可以得到梯度的無偏估計。
\(SGD\) 及相關的小批量亦或更廣義的基於梯度優化的在線學習算法,一個重要的性質是每一步更新的計算時間不依賴訓練樣本數目的多寡。即使訓練樣本非常大時,它們也能收斂。對於足夠大的數據集, \(SGD\) 可能會在處理整個訓練集之前就收斂到最終測試誤差的某個固定容差范圍內。
在求數值解的優化算法中,先選取一組模型參數的初始值,如隨機選取;接下來對參數進行多次迭代,使每次迭代都可能降低損失函數的值。在每次迭代中,先隨機均勻采樣一個由固定數目訓練數據樣本所組成的小批量(mini-batch)\(\mathcal{B}\),然后求小批量中數據樣本的平均損失有關模型參數的導數(梯度),最后用此結果與預先設定的一個正數的乘積作為模型參數在本次迭代的減小量。
模型的每個參數將作如下迭代:
在上式中,\(|\mathcal{B}|\)代表每個小批量中的樣本個數(批量大小,batch size),\(\alpha\)稱作學習率(learning rate)並取正數。
Vectorized
在模型訓練或預測時,常常會同時處理多個數據樣本並用到矢量計算。
廣義上講,當數據樣本數為\(m\),特征數為\(n\)時,矢量計算表達式為
其中模型輸出 \(\hat{y}\in\mathbb{R}^{m\times1}\) , 批量數據樣本特征 \(X\in\mathbb{R}^{m\times n}\) ,權重 \(w\in\mathbb{R}^{n\times 1}\) , 偏差 \(b \in \mathbb{R}\) 。相應地,批量數據樣本標簽 \(\boldsymbol{y} \in \mathbb{R}^{\mathcal{B} \times 1}\) 。設模型參數 \(\boldsymbol{\theta}\) ,我們可以重寫損失函數為
小批量隨機梯度下降的迭代步驟將相應地改寫為
def unroll(params, input_layer_size, hidden_layer_size, num_labels):
Theta1 = params[0:hidden_layer_size * (input_layer_size + 1)].reshape(hidden_layer_size, input_layer_size + 1, order='F')
Theta2 = params[hidden_layer_size * (input_layer_size + 1) : ].reshape(num_labels, hidden_layer_size + 1, order='F')
return Theta1, Theta2
def roll(Theta1, Theta2):
l = [np.ravel(Theta1, order='F').tolist(), np.ravel(Theta2, order='F').tolist()]
l = l[0] + l[1]
return np.array(l)
def predict(params, X):
Theta1 = params[0:hidden_layer_size * (input_layer_size + 1)].reshape(hidden_layer_size, input_layer_size + 1, order='F')
Theta2 = params[hidden_layer_size * (input_layer_size + 1) : ].reshape(num_labels, hidden_layer_size + 1, order='F')
m = X.shape[0] # 樣本的個數
# 執行前向傳播
# 給 X 增加 1 列,全為 1
a1 = np.ones((m, 1))
a1 = np.hstack((a1, X))
# 前向傳播
z2 = np.matmul(a1, Theta1.T)
a2 = sigmoid(z2)
# 增加第二層(隱藏層)的偏置單元
ones = np.ones((m, 1), dtype=np.int)
a2 = np.hstack((ones, a2))
z3 = np.matmul(a2, Theta2.T)
a3 = sigmoid(z3) # m * num_labels
return a3
def predict_accuracy(params, X, y):
a3 = predict(params, X)
# 預測、計算准確率
pre = np.argmax(a3, axis=1) + 1
return np.mean(pre == y) * 100
def predict_graph(params, X, y):
a3 = predict(params, X)
pre = np.argmax(a3, axis=1) + 1
labels = list(range(1, 16))
accuracies = []
for label in labels:
accuracy = np.sum((pre == label) & (y == label)) / np.sum(y == label)
accuracies.append(accuracy)
plt.bar(labels, accuracies)
def mgd(params, learning_rate=0.1, batch_size=100, maxepoch=10):
epoch_list = []
cost_list = []
test_accuracy_list = []
train_accuracy_list = []
m = X.shape[0]
# cost
fig1, axes1 = plt.subplots()
# accuracy
fig2, axes2 = plt.subplots()
for epoch in range(maxepoch):
# 小批量進行更新全部訓練集一次
for i in range(0, m, batch_size):
XX = X[i:i + batch_size, :]
yy = y[i:i + batch_size]
J, grad = nnCostFun(params, XX, yy)
# 更新參數
params -= learning_rate * grad
# 記錄
epoch_list.append(epoch)
cost_list.append(J)
test_accuracy = predict_accuracy(params, Xtest, ytest)
test_accuracy_list.append(test_accuracy)
train_accuracy = predict_accuracy(params, X, y)
train_accuracy_list.append(train_accuracy)
print(epoch, J, train_accuracy, test_accuracy)
axes1.plot(epoch_list, cost_list)
axes1.set_xlabel('epoch')
axes1.set_ylabel('cost')
axes1.text(0, J, 'learning rate: %f' %learning_rate)
axes2.plot(epoch_list, test_accuracy_list)
axes2.plot(epoch_list, train_accuracy_list)
axes2.set_xlabel('epoch')
axes2.set_ylabel('accuracy')
axes2.text(0, test_accuracy, 'learning rate: %f' %learning_rate)
axes2.legend(['test accuracy', 'train accuracy'], loc=0) # 圖例
fig1.savefig('cost_' + str(int(time.time())))
fig2.savefig('accuracy_' + str(int(time.time())))
return cost_list, test_accuracy_list, train_accuracy_list
不同的學習速率、正則化項對學習曲線(learning curve)的表現形式都有所不同
下面對三個不同的學習速率進行簡單的測試,可以看到不同的學習速率對梯度下降的性能是有着非常大的影響的;一個經常的方法是在迭代中逐漸修改學習速率,可以較好避免噪聲對迭代的影響。
mat = io.loadmat('ex4weights2.mat')
Theta1 = mat.get('Theta1')
Theta2 = mat.get('Theta2')
params = roll(Theta1, Theta2)
在執行完小批量的梯度下降后,可以看到訓練集的准確率為 \(96\%\) , 每個數字的准確率如條形圖所示
predict_graph(params, X, y)
predict_accuracy(params, X, y)
96.86666666666667
測試集的准確率為 \(70\%\) ,對於非二分類的問題,這里是十五個類別,\(70\%\) 的准確率不算太低
predict_graph(params, Xtest, ytest)
predict_accuracy(params, Xtest, ytest)
70.44444444444444
下面從測試集中隨機抽取 \(10\) 張圖片進行識別,查看效果
label_map = {1:'零', 2:'一', 3:'二', 4:'三', 5:'四', 6:'五', 7:'六', 8:'七', 9:'八', 10:'九', 11:'十', 12:'百', 13:'千', 14:'萬', 15:'億'}
fig, axes = plt.subplots(nrows=2, ncols=5)
i = 0
for axe in axes:
for ax in axe:
i = np.random.randint(0, 4500)
ax.imshow(Xtest[i, :].reshape((64, 64), order='F'), cmap='gray')
print(label_map.get((predict(params, Xtest[i, :].reshape((1, 4096), order='F')).argmax() + 1)), end=' ')
九 零 千 六 七 千 十 四 三 七
可以看到仍然有些誤差,這是非常正常的,因為此模型選用的 \(BP\) 神經網絡,對於這種圖像識別的,使用卷積神經網絡才是非常有力的模型。
下面可以使用 Photoshop 軟件創建一個 \(64\times 64\) 的灰度圖,並在圖里手寫一個漢字數字,然后放到神經網絡里進行識別
ttt = img.imread('tt.jpg')
plt.imshow(ttt, cmap='gray')
pp = predict(params, ttt.reshape((1, 4096), order='F'))
idx = np.argmax(pp, axis=1) + 1
print(label_map[idx[0]])
八