神经网络实现异或

什么是异或

---

在数字逻辑中，异或是对两个运算元的一种逻辑分析类型，符号为XOR或EOR或⊕。与一般的或（OR）不同，当两两数值相同时为否，而数值不同时为真。异或的真值表如下：

XOR truth table

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，基本上都是可以得到预期的结果的。其实在迭代次数少的情况下，也有可能得到预期的结果，这应该主要取决于初始的参数。