从感知机到神经网络：Python实现与测试

本文为吴恩达神经网络课程打的学习笔记，包含了本人自己的完整推导过程，并且为便于理解，将代码进行了部分重构。

简介

神经网络算法可以算是一种仿生学，其基本数据处理单元——感知机——模仿的是生物神经系统内的神经元，它能够接受来自多个源的信号输入，然后将信号转化为便于传播的信号在进行输出(在生物体内表现为电信号)。

感知机

感知机结构：



{z=wT∗x+bhw(x)=step(z)

其中z为神经元将接收到的信号进行整合，step()函数将整合后的信号进行转化并输出。

早期常用到的step函数有：

heaviside(x)={01x<0x≥0

sgn(x)=⎧⎩⎨−101x<0x=0x>0

多层感知机(Multi-Layer Perceptron)

多层感知机由输入层、输出层和隐含层组成：



实际上在机器学习发展的早期，感知机模型只能用来解决一些鸡毛蒜皮的问题，主要原因就是受step()函数的限制，因为早期的step()函数只能产生离散的输出，并且掩盖了原信号(即无法从感知机的输出来推断原信号)。下图为感知机模型解决异或问题的方案：

反向传播算法的出现

在很长一段时间里人们都没有找到训练MLP的方法，直到1986年D. E. Rumelhart等人发明了反向传播算法，MLP才具有学习能力。

反向传播算法具体过程是这样的：

对于一个训练用例，网络将其作为输出并计算输出，即预测值

测量网络的输出误差

计算在隐含层的最后一层中哪些单元对误差产生了贡献

再从那些产生误差的单元往回追溯，找出隐含层中所有导致误差的单元，直到回溯到输入层为止

调节那些单元的权重分布

为了使这个算法正常工作，作者还对MLP架构做出了一个至关重要的改进：将step函数替换为逻辑函数σ(x)=11+e(−x)，这个函数不仅是全局可导的，而且具有连续输出。

目前常用的激活函数有三种：

⎧⎩⎨⎪⎪σ(x)=11+e(−x)tanh(x)=2σ(2x)−1ReLU(x)=max(0,x)

下面是它们的图象与导数：

单感知机模型的实现

外部库

import numpy as np
import matplotlib.pyplot as plt
import h5py
import scipy
from PIL import Image
from scipy import ndimage
%matplotlib inline

np.random.seed(1)

引入数据

所用数据为h5格式的关于猫的图片数据，每张图片像素为64*64，RGB通道。其中猫的图片标记为y=1，不是猫的图片标记为y=0。

def load_dataset():
train_dataset = h5py.File('datasets/train_catvnoncat.h5', "r")
train_set_x_orig = np.array(train_dataset["train_set_x"][:]) # your train set features
train_set_y_orig = np.array(train_dataset["train_set_y"][:]) # your train set labels

test_dataset = h5py.File('datasets/test_catvnoncat.h5', "r")
test_set_x_orig = np.array(test_dataset["test_set_x"][:]) # your test set features
test_set_y_orig = np.array(test_dataset["test_set_y"][:]) # your test set labels

classes = np.array(test_dataset["list_classes"][:]) # the list of classes

train_set_y_orig = train_set_y_orig.reshape((1, train_set_y_orig.shape[0]))
test_set_y_orig = test_set_y_orig.reshape((1, test_set_y_orig.shape[0]))

return train_set_x_orig, train_set_y_orig, test_set_x_orig, test_set_y_orig, classes

train_set_x_orig, train_set_y, test_set_x_orig, test_set_y, classes = load_dataset()
print("X_train shape:\t{}\nY_train shape:\t{}\nX_test shape:\t{}\nY_test shape:\t{}\n".format
(train_set_x_orig.shape,train_set_y.shape,test_set_x_orig.shape,test_set_y.shape))

输出为：

X_train shape: (209, 64, 64, 3)

Y_train shape: (1, 209)

X_test shape: (50, 64, 64, 3)

Y_test shape: (1, 50)

index = 30        #更改index值以查看不同样本
plt.imshow(train_set_x_orig[index])
print ("y = " + str(train_set_y[0,index]))

输出为：

m_train=train_set_x_orig.shape[0]
m_test=test_set_x_orig.shape[0]
num_px=train_set_x_orig.shape[1]
print("m_train:{}\tm_test:{}\tnum_px:{}".format(m_train,m_test,num_px))

输出为：m_train:209 m_test:50 num_px:64

数据的矩阵化

在多样本输入条件下，对于每一个样本

xi=⎡⎣⎢⎢⎢feature1feature2...featuren⎤⎦⎥⎥⎥，都需进行计算wT∗x+b，最简单的方法就是使用循环让模型针对每一个样本进行训练，这样做的时间复杂度为O(样本数∗每个样本的训练次数)。

由于每个样本的计算之间没有依赖关系，可以将数据矩阵化实现并行计算：

X=[x1x2…xm]

，其中每个样本xi均为列向量。

感知机的最终输出

A=[σ(z1)σ(z2)…σ(zm)]，其中zi表示样本转换后的数字信号。

train_set_x_flatten = train_set_x_orig.reshape(m_train,num_px*num_px*3).T
test_set_x_flatten = test_set_x_orig.reshape(m_test,num_px*num_px*3).T

train_set_x = train_set_x_flatten/255
test_set_x = test_set_x_flatten/255

print("X_train_flatten shape:\t{}\nX_test_flatten shape:\t{}\n".format(train_set_x.shape,test_set_x.shape))

输出为：

X_train_flatten shape: (12288, 209)

X_test_flatten shape: (12288, 50)

实现激活函数

σ(x)=11+e(−x)

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

参数矩阵初始化

随即初始化权重矩阵w，而偏置常量b设为1。

#参数：输入层的节点数dim,同样本特征数
#返回：随机初始化的参数字典parameters
def initialize_parameters(dim):
parameters={
"w":np.random.randn(dim,1)*0.01,
"b":1
}

return parameters

前向传播

前向传播就是针对输入给出相应的输出并计算误差的过程，此处使用交叉熵L(a,y)=−ylna−(1−y)ln(1−a)来评估误差：

⎧⎩⎨⎪⎪⎪⎪⎪⎪Z=WT∗X+bA=σ(Z)Loss=1m∑i=1mL(ai,y)

#参数：包含W,b的参数字典parameters，样本数据X，样本标签Y
#返回：损失值cost，包含A的缓存cache
def forward_propagation(parameters, X, Y):
w=parameters['w']
b=parameters['b']

m = X.shape[1]

Z=np.dot(w.T,X)+b
A = sigmoid(Z)
cache=A
cost = ((-1/m)*(Y*np.log(A)+(1-Y)*np.log(1-A))).sum(axis=1)

cost = np.squeeze(cost)
assert(cost.shape == ())

return cost,cache

反向传播

反向传播算法实际上是一个优化问题，它需要找到一组能最小化Loss的解(W,b)。最简单的方法就是求导，直接计算满足一阶导为零且二阶导大于或小于零的参数即可，实际上在线性模型中就有这么一种方法称为正规方程。

另外一种逼近解的方法叫梯度下降法，对于某参数θ的计算过程如下：

θnext_step=θcur_step−α∂Loss∂θ

其中α称为学习率。

具体对于参数w跟b的计算为：

∂L∂A=1m∗(−YA+1−Y1−A)

∂L∂Z===∂L∂A∗∂A∂Z1m∗(−yA+1−Y1−A)∗A∗(1−A)1m(A−Y)

∂L∂W==∂L∂Z∗∂Z∂W1m(A−Y)∗XT

∂L∂b==∂L∂Z∗∂Z∂b1m(A−Y)

我们只需要∂L∂W=1mX∗(A−Y)T和∂L∂b=1m(A−Y)即可。为了减少计算量，可将前向传播时计算过的A缓存起来。

#参数：样本数据X，样本标签Y，包含A的缓存cache
#返回：包含dw,db的梯度字典grads
def backward_propagation(X,Y,cache):
m=X.shape[1]
A=cache

grads = {"dw": np.dot(X,(A-Y).T)/m,
"db": ((A-Y)/m).sum(axis=1)}

return grads

参数的迭代优化

对于给定的学习率与迭代次数，对参数W，b进行优化。optimize函数返回优化好的参数与记录下的损失值cost以供可视化。

θnext_step=θcur_step−α∂Loss∂θ

def optimize(parameters, X, Y, num_iterations, learning_rate):
costs = []

w=parameters['w']
b=parameters['b']

for i in range(num_iterations):
cost,cache=forward_propagation(parameters,X,Y)
grads=backward_propagation(X,Y,cache)

dw = grads["dw"]
db = grads["db"]

w = w-learning_rate*dw
b = b-learning_rate*db

parameters={
'w':w,
'b':b
}

#每一百次迭代记录一次cost
if i % 100 == 0:
costs.append(cost)

params = {"w": w,
"b": b}

return params, costs

预测函数

预测函数以给定的参数(W,b)与样本X为参数，会返回给定参数下的模型对样本X的预测值。

def predict(parameters, X):
w = parameters["w"]
b = parameters["b"]

m = X.shape[1]
Y_prediction = np.zeros((1,m))

Z=np.dot(w.T,X)+b
A = sigmoid(Z)

for i in range(A.shape[1]):
Y_prediction[0, i] = 1 if A[0,i]>0.5 else 0    #注意取矩阵元素与取数组元素的区别

return Y_prediction

评估模型

可视化模型在训练过程中的学习曲线。

def show_info(parameters,X_train,Y_train,X_test,Y_test,costs,learning_rate):
Y_prediction_test = predict(parameters, X_test)
Y_prediction_train = predict(parameters, X_train)

print("train accuracy: "  + str(np.sum((Y_prediction_train == Y_train)/X_train.shape[1])))
print("test accuracy: "  + str(np.sum((Y_prediction_test == Y_test)/X_test.shape[1])))

plt.plot(np.squeeze(costs))
plt.ylabel('cost')
plt.xlabel('iterations (per tens)')
plt.title("Learning rate =" + str(learning_rate))
plt.show()

整合模型

def model(X_train, Y_train, X_test, Y_test, num_iterations = 2000, learning_rate = 0.005):
#1.参数初始化
parameters = initialize_parameters(X_train.shape[0])

#2.迭代优化参数
parameters, costs = optimize(parameters, X_train, Y_train, num_iterations, learning_rate)

show_info(parameters,X_train,Y_train,X_test,Y_test,costs,learning_rate)
return parameters

训练模型

parameters = model(train_set_x, train_set_y, test_set_x, test_set_y, num_iterations = 2000, learning_rate = 0.005)

输出为：

多层网络实现

网络模型

此处实现一个L层的神经网络，最后一层的激活函数使用σ(x)，其余层均使用ReLU(x)。整个模型框图如下图所示：

参数矩阵维度分析

对于L层的网络模型，易得每一层(第i层)都有这么几个参数(i>0)：

神经元数，n[i]

权重参数，W[i]

偏置参数，b

线性输出，Z[i]=W[i]∗A[i−1]+b（注意，此处不再使用WT∗X）

激活函数，g[i](x)={σ(x)ReLU(x)i=Li<L

输出，A[i]=g(Z[i])

下面简单推导一下每一层的参数矩阵的维度。

令每个样本均为列向量，则易得

dim(A[i])=(n[i],m)

其中m为样本数量。

然后从输入层开始，输入层的矩阵维度为dim(A[0])=(n[0],m)，而

dim(W[1]∗A[0])=dim(A[1])=(n[1],m)

所以dim(W[1])=(n[1],n[0])；继续沿层数推下去会发现一个规律，得出

dim(W[i])=(n[i],n[i−1])

参数矩阵初始化

#三层模型
layers_dims = [12288, 7,1]         #输入层，隐含层，输出层

def initialize_parameters(layer_dims):
parameters = {}
L = len(layer_dims)            # number of layers in the network

for i in range(1, L):
parameters['W' + str(i)] = np.random.randn(layer_dims[i],layer_dims[i-1])*0.01
parameters['b' + str(i)] = np.ones((layer_dims[i],1))

assert(parameters['W' + str(i)].shape == (layer_dims[i], layer_dims[i-1]))
assert(parameters['b' + str(i)].shape == (layer_dims[i], 1))

return parameters

激活函数

def sigmoid(Z):
A = 1/(1+np.exp(-Z))
return A

def relu(Z):
A = np.maximum(0,Z)
assert(A.shape == Z.shape)
return A

每一层神经元的运算

完成每一层所有神经元的运算，支持两种激活函数。

#参数：上一层的输出A_prev，当前层的参数W,b，激活方式activation
#返回：当前层的输出A，缓存A_prev,W,Z
def perceptron_forward(A_prev, W, b, activation):
Z=np.dot(W,A_prev)

if activation == "sigmoid":
A = sigmoid(Z)

elif activation == "relu":
A = relu(Z)

assert (A.shape == (W.shape[0], A_prev.shape[1]))
cache=(A_prev,W,Z)

return A, cache

前向传播(L层的神经元运算)

def forward_propagation(X,Y, parameters):
caches = []
A = X        #第零层的输出就为X
L = len(parameters)//2                  # number of layers in the neural network

for l in range(1, L):
A_prev = A

A, cache = perceptron_forward(A_prev,parameters["W"+str(l)],parameters["b"+str(l)],"relu")

caches.append(cache)

AL, cache = perceptron_forward(A,parameters["W"+str(L)],parameters["b"+str(L)],"sigmoid")

caches.append(cache)

assert(AL.shape == (1,X.shape[1]))

m = Y.shape[1]
cost = -np.sum(np.multiply(Y,np.log(AL))+np.multiply(1-Y,np.log(1-AL)))/m

return AL, caches,cost

反向传播推导

首先从最后一层开始，

dA[L]==∂L∂A[L]1m∗(−YA[L]+1−Y1−A[L])

dZ[L]===∂L∂Z[L]dA[L]∗∂A[L]∂Z[L]dA[L]∗A[L]∗(1−A[L])

dW[L]===∂L∂W[L]dZ[L]∗∂Z[L]∂W[L]dZ[L]∗A[L−1].T

db[L]=1mdZ[L]

于是对于使用\sigma(x)的最后一层，得到了dW[L]=dZ[L]∗A[L−1].T和db[L]=1mdZ[L]，继续往前推导。

dA[L−1]===∂L∂A[L−1]dZ[L]∗∂Z[L]∂A[L−1]W[L].T∗dZ[L]

dZ[L−1]===∂L∂Z[L−1]dA[L−1]∗∂A[L−1]∂Z[L−1](0,dA[L−1])

dW[L−1]===∂L∂W[L−1]dZ[L−1]∗∂Z[L−1]∂W[L−1]dZ[L−1]∗A[L−2].T

db[L−1]=1mdZ[L−1]

可得对于使用ReLU(x)的L-1层，每层的梯度为dW[i]=dZ[i]∗A[i−1].T，db[i]=1mdZ[i]。

以上结论可变换为：

⎧⎩⎨⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪dA[L]dZ[L]dW[i]db[i]dA[i−1]dZ[i−1]=1m∗(−YA[L]+1−Y1−A[L])=dA[L]∗A[L]∗(1−A[L])=dZ[i]∗A[i−1].T=1mdZ[i]=W[i].T∗dZ[i]=(0,dA[i−1])(1)(2)(3)(4)(5)(6)

在多层神经网络的反向传播过程中，为快速求出第i层的dW[i]与db[i]，需要缓存前向传播过程中的A[i−1]、W[i]与Z[i]，还有反向传播过程中的dA[i]。

每一层神经元的运算

#参数：当前层的dA，包含A_prev,W,Z的cache，激活方式activation
#返回：上一层的梯度dA_prev，当前层的梯度dW,db
def perceptron_backward(dA, cache, activation):
A_prev,W,Z=cache

if activation == "relu":
dZ = np.array(dA, copy=True)
dZ[Z <= 0] = 0

elif activation == "sigmoid":
s = 1/(1+np.exp(-Z))
dZ = dA * s * (1-s)

assert (dZ.shape == Z.shape)

m = A_prev.shape[1]

dW = np.dot(dZ,A_prev.T)
db = np.sum(dZ,axis=1,keepdims=True)/m
dA_prev = np.dot(W.T,dZ)

return dA_prev, dW, db

反向传播

def backward_propagation(AL, Y, caches):
grads = {}
L = len(caches) # the number of layers
m = AL.shape[1]

dAL = - (np.divide(Y, AL) - np.divide(1 - Y, 1 - AL))/m
current_cache = caches[L-1]        #取L层的缓存数据
grads["dA" + str(L-1)], grads["dW" + str(L)], grads["db" + str(L)] = perceptron_backward(dAL,current_cache,"sigmoid")

for l in reversed(range(L - 1)):
current_cache = caches[l]
dA_prev, dW, db = perceptron_backward(grads["dA"+str(l+1)],current_cache,"relu")

grads["dA" + str(l)] = dA_prev
grads["dW" + str(l + 1)] = dW
grads["db" + str(l + 1)] = db

return grads

迭代优化

def update_parameters(parameters, grads, learning_rate):
L = len(parameters) // 2 # number of layers in the neural network

for l in range(L):
parameters["W" + str(l+1)] = parameters["W"+str(l+1)]-learning_rate*grads["dW"+str(l+1)]
parameters["b" + str(l+1)] = parameters["b"+str(l+1)]-learning_rate*grads["db"+str(l+1)]

return parameters

预测函数

若\sigma(x)的输出>0.5则判为正例，否则判为反例。

def predict(X, y, parameters):
m = X.shape[1]
n = len(parameters) // 2 # number of layers in the neural network
p = np.zeros((1,m))

probas, _ , _ = forward_propagation(X,y, parameters)

for i in range(0, probas.shape[1]):
if probas[0,i] > 0.5:
p[0,i]=1 if probas[0,i]>0.5 else 0

print("Accuracy: "  + str(np.sum((p == y)/m)))

return p

整合模型

def L_layer_model(X, Y, layers_dims, learning_rate = 0.005, num_iterations = 2000, print_cost=False):
costs = []

parameters = initialize_parameters(layers_dims)

for i in range(0, num_iterations):
AL, caches,cost = forward_propagation(X,Y, parameters)

grads = backward_propagation(AL, Y, caches)

parameters = update_parameters(parameters, grads, learning_rate)

if i % 100 == 0:
costs.append(cost)

# plot the cost
plt.plot(np.squeeze(costs))
plt.ylabel('cost')
plt.xlabel('iterations (per tens)')
plt.title("Learning rate =" + str(learning_rate))
plt.show()

return parameters

训练模型

parameters = L_layer_model(train_set_x, train_set_y, layers_dims, learning_rate=0.005,num_iterations = 2000, print_cost = True)

输出为：

评估模型

predict(train_set_x,train_set_y,parameters)
predict(test_set_x,test_set_y,parameters)

输出为：

Accuracy: 0.99043062201

Accuracy: 0.76

具有隐含层的神经网络比单感知机要提高了六个百分点的准确率，不过两个模型都存在明显的过拟合现象。