什么是異或
在數字邏輯中，異或是對兩個運算元的一種邏輯分析類型，符號為XOR或EOR或⊕。與一般的或（OR）不同，當兩兩數值相同時為否，而數值不同時為真。異或的真值表如下：
| Input | Output | |
|---|---|---|
| A | B | |
| 0 | 0 | 0 |
| 0 | 1 | 1 |
| 1 | 0 | 1 |
| 1 | 1 | 0 |
- 0, false
- 1, true
據說在人工神經網絡（artificial neural network, ANN）發展初期，由於無法實現對多層神經網絡（包括異或邏輯）的訓練而造成了一場ANN危機，到最后BP算法的出現，才讓訓練帶有隱藏層的多層神經網絡成為可能。因此異或的實現在ANN的發展史是也是具有里程碑意義的。異或之所以重要，是因為它相對於其他邏輯關系，例如與（AND）, 或（OR）等，異或是線性不可分的。如下圖：
在實際應用中，異或門（Exclusive-OR gate, XOR gate）是數字邏輯中實現邏輯異或的邏輯門，這一函數能實現模為2的加法。因此，異或門可以實現計算機中的二進制加法。
異或的神經網絡結構
在【機器學習】課程中，使用了AND（與），NOR（或非）和OR（或）的組合實現了XNOR（同或），與我們要實現的異或（XOR）正好相反。因此還是可以采用課程中的神經網絡結構，如下圖：
如果算上輸入層我們的網絡共有三層，如下圖所示，其中第1層和第2層中的1分別是這兩層的偏置單元。連線上是連接前后層的參數。
- 輸入：我們一共有四個訓練樣本，每個樣本有兩個特征，分別是(0, 0), (1, 0), (0, 1), (1, 1);
- 理想輸出：參考上面的真值表，樣本中兩個特征相同時為0，相異為1
- 參數：隨機初始化，范圍為(-1, 1)
- 關於神經網絡的基礎知識以及前向傳播、反向傳播的實現請參考下面兩篇文章，寫的非常精彩：
機器學習公開課筆記(4)：神經網絡(Neural Network)——表示
機器學習公開課筆記(5)：神經網絡(Neural Network)——學習
代碼
- 原生態的代碼：
下面的實現是完全根據自己的理解和對【機器學習】課程中作業題的模仿而寫成的，雖然代碼質量不是非常高，但是算法的所有細節都展示出來了。
在66, 69, 70行的注釋是我之前沒有得到正確結果的三個原因，其中epsilon確定的是隨機初始化參數的范圍，例如epsilon=1，參數范圍就是(-1, 1)
1 # -*- coding: utf-8 -*-
2 """
3 Created on Tue Apr 4 10:47:51 2017
4
5 @author: xin
6 """
7 # Neural Network for XOR
8 import numpy as np
9 import matplotlib.pyplot as plt
10
11 HIDDEN_LAYER_SIZE = 2
12 INPUT_LAYER = 2 # input feature
13 NUM_LABELS = 1 # output class number
14 X = np.array([[0, 0], [0, 1], [1, 0], [1, 1]])
15 y = np.array([[0], [1], [1], [0]])
16
17
18 def rand_initialize_weights(L_in, L_out, epsilon):
19 """
20 Randomly initialize the weights of a layer with L_in
21 incoming connections and L_out outgoing connections;
22
23 Note that W should be set to a matrix of size(L_out, 1 + L_in) as
24 the first column of W handles the "bias" terms
25 """
26 epsilon_init = epsilon
27 W = np.random.rand(L_out, 1 + L_in) * 2 * epsilon_init - epsilon_init
28 return W
29
30
31 def sigmoid(x):
32 return 1.0 / (1.0 + np.exp(-x))
33
34
35 def sigmoid_gradient(z):
36 g = np.multiply(sigmoid(z), (1 - sigmoid(z)))
37 return g
38
39
40 def nn_cost_function(theta1, theta2, X, y):
41 m = X.shape[0] # m=4
42 # 計算所有參數的偏導數（梯度）
43 D_1 = np.zeros(theta1.shape) # Δ_1
44 D_2 = np.zeros(theta2.shape) # Δ_2
45 h_total = np.zeros((m, 1)) # 所有樣本的預測值, m*1, probability
46 for t in range(m):
47 a_1 = np.vstack((np.array([[1]]), X[t:t + 1, :].T)) # 列向量, 3*1
48 z_2 = np.dot(theta1, a_1) # 2*1
49 a_2 = np.vstack((np.array([[1]]), sigmoid(z_2))) # 3*1
50 z_3 = np.dot(theta2, a_2) # 1*1
51 a_3 = sigmoid(z_3)
52 h = a_3 # 預測值h就等於a_3, 1*1
53 h_total[t,0] = h
54 delta_3 = h - y[t:t + 1, :].T # 最后一層每一個單元的誤差, δ_3, 1*1
55 delta_2 = np.multiply(np.dot(theta2[:, 1:].T, delta_3), sigmoid_gradient(z_2)) # 第二層每一個單元的誤差（不包括偏置單元）, δ_2, 2*1
56 D_2 = D_2 + np.dot(delta_3, a_2.T) # 第二層所有參數的誤差, 1*3
57 D_1 = D_1 + np.dot(delta_2, a_1.T) # 第一層所有參數的誤差, 2*3
58 theta1_grad = (1.0 / m) * D_1 # 第一層參數的偏導數，取所有樣本中參數的均值，沒有加正則項
59 theta2_grad = (1.0 / m) * D_2
60 J = (1.0 / m) * np.sum(-y * np.log(h_total) - (np.array([[1]]) - y) * np.log(1 - h_total))
61 return {'theta1_grad': theta1_grad,
62 'theta2_grad': theta2_grad,
63 'J': J, 'h': h_total}
64
65
66 theta1 = rand_initialize_weights(INPUT_LAYER, HIDDEN_LAYER_SIZE, epsilon=1) # 之前的問題之一，epsilon的值設置的太小
67 theta2 = rand_initialize_weights(HIDDEN_LAYER_SIZE, NUM_LABELS, epsilon=1)
68
69 iter_times = 10000 # 之前的問題之二，迭代次數太少
70 alpha = 0.5 # 之前的問題之三，學習率太小
71 result = {'J': [], 'h': []}
72 theta_s = {}
73 for i in range(iter_times):
74 cost_fun_result = nn_cost_function(theta1=theta1, theta2=theta2, X=X, y=y)
75 theta1_g = cost_fun_result.get('theta1_grad')
76 theta2_g = cost_fun_result.get('theta2_grad')
77 J = cost_fun_result.get('J')
78 h_current = cost_fun_result.get('h')
79 theta1 -= alpha * theta1_g
80 theta2 -= alpha * theta2_g
81 result['J'].append(J)
82 result['h'].append(h_current)
83 # print(i, J, h_current)
84 if i==0 or i==(iter_times-1):
85 print('theta1', theta1)
86 print('theta2', theta2)
87 theta_s['theta1_'+str(i)] = theta1.copy()
88 theta_s['theta2_'+str(i)] = theta2.copy()
89
90 plt.plot(result.get('J'))
91 plt.show()
92 print(theta_s)
93 print(result.get('h')[0], result.get('h')[-1])
下面是輸出結果：
# 隨機初始化得到的參數
('theta1', array([[ 0.18589823, -0.77059558, 0.62571502], [-0.79844165, 0.56069914, 0.21090703]])) ('theta2', array([[ 0.1327994 , 0.59513332, 0.34334931]]))
# 訓練后得到的參數 ('theta1', array([[-3.90903729, -7.44497437, 7.20130773], [-3.76429211, 6.93482723, -7.21857912]])) ('theta2', array([[ -6.5739346 , 13.33011993, 13.3891608 ]]))
# 同上，第一次迭代和最后一次迭代得到的參數 {'theta1_0': array([[ 0.18589823, -0.77059558, 0.62571502], [-0.79844165, 0.56069914, 0.21090703]]), 'theta2_9999': array([[ -6.5739346 , 13.33011993, 13.3891608 ]]), 'theta1_9999': array([[-3.90903729, -7.44497437, 7.20130773], [-3.76429211, 6.93482723, -7.21857912]]), 'theta2_0': array([[ 0.1327994 , 0.59513332, 0.34334931]])}
# 預測值h: 第1個array里是初始參數預測出來的值，第2個array中是最后一次得到的參數預測出來的值
(array([[ 0.66576877], [ 0.69036552], [ 0.64994307], [ 0.67666546]]), array([[ 0.00245224], [ 0.99812746], [ 0.99812229], [ 0.00215507]]))
下面是隨着迭代次數的增加，代價函數值J(θ)的變化情況：
- 更加精煉的代碼
下面這段代碼是我在排除之前自己的代碼中的問題時，在Stack Overflow上發現的，發帖的人也碰到了同樣的問題，但原因不一樣。他的代碼里有一點小問題，已經修正。這段代碼，相對於我自己的原生態代碼，有了非常大的改進，沒有限定層數和每層的單元數，代碼本身也比較簡潔。
說明：由於第44行，傳的參數是該層的a值，而不是z值，所以第11行需要做出一點修改，其實直接傳遞a值是一種更方便的做法。
1 # -*- coding: utf-8 -*-
2
3 import numpy as np
4 import matplotlib.pyplot as plt
5
6
7 def sigmoid(x):
8 return 1/(1+np.exp(-x))
9
10 def s_prime(z):
11 return np.multiply(z, 1.0-z) # 修改的地方
12
13 def init_weights(layers, epsilon):
14 weights = []
15 for i in range(len(layers)-1):
16 w = np.random.rand(layers[i+1], layers[i]+1)
17 w = w * 2*epsilon - epsilon
18 weights.append(np.mat(w))
19 return weights
20
21 def fit(X, Y, w):
22 # now each para has a grad equals to 0
23 w_grad = ([np.mat(np.zeros(np.shape(w[i])))
24 for i in range(len(w))]) # len(w) equals the layer number
25 m, n = X.shape
26 h_total = np.zeros((m, 1)) # 所有樣本的預測值, m*1, probability
27 for i in range(m):
28 x = X[i]
29 y = Y[0,i]
30 # forward propagate
31 a = x
32 a_s = []
33 for j in range(len(w)):
34 a = np.mat(np.append(1, a)).T
35 a_s.append(a) # 這里保存了前L-1層的a值
36 z = w[j] * a
37 a = sigmoid(z)
38 h_total[i, 0] = a
39 # back propagate
40 delta = a - y.T
41 w_grad[-1] += delta * a_s[-1].T # L-1層的梯度
42 # 倒過來，從倒數第二層開始到第二層結束，不包括第一層和最后一層
43 for j in reversed(range(1, len(w))):
44 delta = np.multiply(w[j].T*delta, s_prime(a_s[j])) # 這里傳遞的參數是a，而不是z
45 w_grad[j-1] += (delta[1:] * a_s[j-1].T)
46 w_grad = [w_grad[i]/m for i in range(len(w))]
47 J = (1.0 / m) * np.sum(-Y * np.log(h_total) - (np.array([[1]]) - Y) * np.log(1 - h_total))
48 return {'w_grad': w_grad, 'J': J, 'h': h_total}
49
50
51 X = np.mat([[0,0],
52 [0,1],
53 [1,0],
54 [1,1]])
55 Y = np.mat([0,1,1,0])
56 layers = [2,2,1]
57 epochs = 5000
58 alpha = 0.5
59 w = init_weights(layers, 1)
60 result = {'J': [], 'h': []}
61 w_s = {}
62 for i in range(epochs):
63 fit_result = fit(X, Y, w)
64 w_grad = fit_result.get('w_grad')
65 J = fit_result.get('J')
66 h_current = fit_result.get('h')
67 result['J'].append(J)
68 result['h'].append(h_current)
69 for j in range(len(w)):
70 w[j] -= alpha * w_grad[j]
71 if i == 0 or i == (epochs - 1):
72 # print('w_grad', w_grad)
73 w_s['w_' + str(i)] = w_grad[:]
74
75
76 plt.plot(result.get('J'))
77 plt.show()
78 print(w_s)
79 print(result.get('h')[0], result.get('h')[-1])
下面是輸出的結果：
# 第一次迭代和最后一次迭代得到的參數
{'w_4999': [matrix([[ 1.51654104e-04, -2.30291680e-04, 6.20083292e-04],
[ 9.15463982e-05, -1.51402782e-04, -6.12464354e-04]]), matrix([[ 0.0004279 , -0.00051928, -0.00042735]])],
'w_0': [matrix([[ 0.00172196, 0.0010952 , 0.00132499],
[-0.00489422, -0.00489643, -0.00571827]]), matrix([[-0.02787502, -0.01265985, -0.02327431]])]}
# 預測值h: 第1個array里是初始參數預測出來的值，第2個array中是最后一次得到的參數預測出來的值
(array([[ 0.45311095],
[ 0.45519066],
[ 0.4921871 ],
[ 0.48801121]]),
array([[ 0.00447994],
[ 0.49899856],
[ 0.99677373],
[ 0.50145936]]))
觀察上面的結果，最后一次迭代得到的結果並不是我們期待的結果，也就是第1、4個值接近於0, 第2、3個值接近於1。下面是代價函數值J(θ)隨着迭代次數增加的變化情況：
從上圖可以看到，J(θ)的值從2000以后就一直停留在0.35左右，因此整個網絡有可能收斂到了一個局部最優解，也有可能是迭代次數不夠導致的。
將迭代次數改成10000后， 即epochs = 10000，基本上都是可以得到預期的結果的。其實在迭代次數少的情況下，也有可能得到預期的結果，這應該主要取決於初始的參數。