DNN网络的学习和封装实现（任意层全连接网络）

 这是之前自己实现的一个任意层全连接网络的模块，主要是各种正向、反向函数，主要依赖numpy

 基本功能已验证是没问题的，由于后面转而使用开源的框架，缺乏进一步验证，所以有些内容仅供参考

import numpy as np
import matplotlib.pyplot as plt
import sys
import copy

#@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
# activation functions and their derivative
#@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
def sigmoid(Z):
	A = 1 / (1 + np.exp(-Z))
	return A

def sigmoid_derivative(Z, A=None):
	if A is None:
		A = 1 / (1 + np.exp(-Z))
	dA_dZ = A * (1 - A)
	return dA_dZ

def tanh(Z):
	# 这个推荐直接使用np.tanh()
	# tanh = sinh/cosh, 下面的式子由欧拉公式推出
	ep = np.exp(Z)
	en = np.exp(-Z)
	A = (ep - en) / (ep + en)
	return A

def tanh_derivative(Z, A=None):
	if A is None:
		ep = np.exp(Z)
		en = np.exp(-Z)
		A = (ep - en) / (ep + en)

	dA_dZ = 1 - A**2
	return dA_dZ

def relu(Z):
	A = np.maximum(0, Z) 	# 0会被广播成相同的维度，逐一比较，取大值
							# 区别：max(a, axis=?)这个函数是取a的最大值，可以选择维度
	#a = np.where(z < 0, 0, z) # 这个也是可行的，不过不够效率（0广播，x<0为True取矩阵1值，否则取矩阵2值）
	return A

def relu_derivative(Z, A=None):
	dA_dZ = (Z > 0) * 1.0
	return dA_dZ

def leaky_relu(Z, leak=0.01):
	'''
	# 网上的一个实现
	f1 = 0.5 * (1 + leak)
	f2 = 0.5 * (1 - leak)
	return f1 * x + f2 * abs(x)
	'''
	'''
	# 采用下面更简洁的写法，显然网上的实现更有效率
	p = x > 0
	n = np.logical_not(p)
	s = x * p * 1.0 + x * n * leak
	'''
	A = Z * (Z > 0) * 1.0 + Z * (Z <= 0) * leak
	return A

def leaky_relu_derivative(Z, leak=0.2, A=None):
	dA_dZ = (Z > 0) * 1.0 + (Z <= 0) * leak
	return dA_dZ

def softmax(Z):
	''' x: (N, M) '''
	z_exp = np.exp(Z)
	ez_sum = np.sum(z_exp, axis=0, keepdims=True)
	A = z_exp / ez_sum
	return A

def softmax_derivative(Z, dA_input, A):
	'''
	返回dJ/dx的导数，目前测试和dx = A-Y效果相同
	'''
	'''
	# 其它激活函数都是单输出，这个是多输出，有些差异
	# 较复杂，一般和损失函数一起提前计算出来
	#https://blog.csdn.net/sunlanchang/article/details/88822573
	#https://www.zhihu.com/question/360140818/answer/1006428058

	目前已知：dJ/dZ = A - y（sigmoid+交叉熵，softmax+交叉熵求导化归都是这个通式），推导过程如下：
	# 推导（单个样本，三个类别）:
	z1 --exp--> ez1 	-->		a1 = ez1 / (ez1 + ez2 + ez3)
	z2 --exp--> ez2 	-->		a2 = ez2 / (ez1 + ez2 + ez3)
	z3 --exp--> ez3 	-->		a3 = ez2 / (ez1 + ez2 + ez3)
	可见：z1会影响a1, a2, a3，同样z2, z3也影响所有输出a
	因此，根据链式法则法则：dz1 = da1/dz1 + da2/dz1 + da3/dz1，dz2, dz3类似
	具体求法：
		令(ez1 + ez2 + ez3) = sum，要求da1/dz1 = (ez1 / sum)'
		根据导数除法公式：(u(x)/v(x))' = (u'v - v'u) / (v**2)
		这里u(x) = exp(z1), v(x) = ez1 + ez2 + ez3
		（注意这里exp(z1)并不是复合函数，因此结果后面不需要再乘exp(z1)
		复合函数指的是log(exp(x))这种）
		容易求得：
		da1/dz1 = (ez1 * sum - ez1**2) / (sum**2) = a1 - a1**2
		da2/dz1 = -(ez2 * ez1) / (sum**2) = -a2 * a1
		da3/dz1 = -(ez3 * ez1) / (sum**2) = -a3 * a1

		如果a1,a2,a3是F(z)最终输出，则不能直接求导（多输出，不可导），即以下式子没有意义
		dz1 = a1 - a1**2 - a2 * a1 - a3 * a1
			= a1 * (1 - a1 - a2 - a3) = 0
		因为这是一个“三个输入产生三个输出”的问题，将三个导数直接相加没有直接意义
		实际上，softmax后面还要接损失函数，最终构成单输出的函数，则可以运用链式法则求dJ/d?的偏导数。
	此时（单个样本）:
		loss = J = -y1 * log(a1) - y2 * log(a2) - y3 * log(a3)， y1,y2,y3中只有一个1，a1,a2,a3是三个类别的概率
		dJ/da1 = -y1/a1
		dJ/da2 = -y2/a2
		dJ/da3 = -y3/a3
		偏导数：
		dJ/dz1 = dJ/da1 * da1/dz1 + dJ/da2 * da2/dz1 + dJ/da3 * da3/dz1，链式法则
			   = a1 - y1
		同样，dJ/dz2 = a2 - y2, dJ/dz3 = a3 - y3
		即：dJ/dZ = A - Y
		由于softmax各个样本间是独立的，因此多样本也满足dJ/dZ = A - Y
	'''

	# 以下是softmax求导过程的实现，效率较低，但是正确
	N = Z.shape[0]
	M = Z.shape[1]
	K = np.zeros((N, N, M)) # 导数矩阵，第一个N对应N个x，第二个N表示各个局部的分量
	dZ = np.zeros((N, M))
	for i in range(N):		# 对每一个z(i)求导
		for j in range(N):	# N个方向上的导数分量
			if j == i:
				K[i][j] = A[i] - A[i]**2
			else:
				K[i][j] = -A[j] * A[i]
		dZ[i] = np.sum(dA_input * K[i], axis=0) # 列向求和，链式法则
		# 直接返回K矩阵的话，可其它激活函数更加一致，但是会导致函数外面多一个for循环
		# 所以此函数直接要求输出导数dJ/dA，返回dJ/dZ

	return dZ

#@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
# Loss functions
# 由于我们一般手动计算出Loss和final activation的最简形式，所以这里不考虑其导数
#@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
def loss_l1(A, Y):
	'''Y = (No, M)'''
	loss = np.sum(np.abs(A - Y)) / Y.shape[1]
	return loss

def loss_l2(A, Y):
	loss = np.sum((Y - A) ** 2) / Y.shape[1]
	return loss

def loss_cross_entropy(A, Y, activation):
	if activation == 'sigmoid':
		# 对二分类(2, M)，它是：
		loss = np.sum(-Y * np.log(A) - (1-Y) * np.log(1-A)) / Y.shape[1]
	else:
		# 对多分类(No, M)
		loss = np.sum(-Y * np.log(A)) / Y.shape[1]
	return loss

def loss_cross_entropy_derivative(A, Y, activation):
	if activation == 'sigmoid':
		# 对sigmoid分类(2, M)，它是：
		dA = -np.divide(Y, A) + np.divide(1-Y, 1-A)
		# dA = -np.divide(Y, A) - np.divide(1-Y, 1-A) 如果写成这种错误形式，会导致梯度爆炸
		# 为什么第二项的前面是加呢？因为-(1-Y)log(1-A)对A求导，log里面是-A+1，属于复合函数kx+b，所以需要乘以-1
	else:
		# 对softmax分类(No, M)
		dA = -np.divide(Y, A)
	return dA

#@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
# initial functions and others utils
#@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
def normalize_universal(X, normalize_paras=None):
	'''数据集归一化'''
	if normalize_paras is None:
		'''通用的初始化：均值归0，方差归1'''
		'''在调整好维度的训练集x = (Ni, Mepoch)上进行这个调用，对每个特征进行归一化'''
		u = np.sum(X, axis=1, keepdims=True) / X.shape[1]
		X = X - u

		# 本来方差是x**2 - u**2，但是此时u已经归0了
		d2 = np.sum(X**2, axis=1, keepdims=True) / X.shape[1]
		X = X / d2
	else:
		u = normalize_paras['u']
		d2 = normalize_paras['d2']
		X = X - u
		X = X / d2

	# 验证集、测试集也要用同样的参数（训练集得到的u,d2）进行归一化
	normalize_paras = {'u': u, 'd2': d2}
	return X, normalize_paras

def normalize_norm(X, normalize_paras=None):
	# 除范数形式初始化，目前不使用这种形式
	# linalg是liner + algebra，即线性代数, norm即范数
	# norm()的ord参数表示范数类型，默认为l2范数
	# axis=1表示水平方向（第二个维度，0是第一个维度即竖直方向）
	'''
	norm_row_x = np.linalg.norm(x, axis=1, keepdims=True)
	x_normalized = X / norm_row_x
	return x_normalized
	'''
	return X, normalize_paras

def normalize_rgb(X):
	'''rgb图像经常用到的归一化方式'''
	return X / 255

def init_weight(networks, Wl=None, Bl=None):
	if Wl is not None and Bl is not None:
		networks['weights'] = {'Wl': Wl, 'Bl': Bl}
		return networks

	'''
	Nl: 各层的神经元数，包括输入层（特征数）
	activations: 各层激活函数string, 第0层赋''即可
	'''
	Nl = networks['Nl']
	activations = networks['activations']

	layers = len(Nl)
	Wl = [0] * layers
	Bl = [0] * layers

	for l in range(1, layers, 1):
		k = 1.0
		if activations[l] in ['relu', 'leaky_relu']:
			k = 2.0
		Wl[l] = np.random.randn(Nl[l], Nl[l-1]) * np.sqrt(k / Nl[l-1]) # Xavier初始化
		Bl[l] = np.zeros((Nl[l], 1))

	networks['weights'] = {'Wl': Wl, 'Bl': Bl}
	return networks

def flatten(x, y):
	'''扁平化：经常用到的一个函数，定义成虚函数，具体问题具体处理'''
	'''
	M = x.shape[0] 				# x = (M, ?, ?, ?), y = (M, 1)
	X = x.reshape(M, -1).T 		# (N, M)
	Y = y.T 					# (1, M)
	return X, Y
	'''
	return x, y


def softmax_y2mat(Y, No):
	'''
	将实验测试结果（向量）转换成神经网络输出的矩阵形式
	Y：二维数组(1,M)，每个元素的取值范围[0,1,2,3,4..]，
	分别表示"优良中差"No=5个等级
	'''
	M = Y.shape[1]
	Yr = np.zeros((No, M))

	for i in range(No):
		Yr[i] = (i == Y[0])
	return Yr
	'''
	# 以下算法可以达到同样目的
	# np.eye(C)产生CxC的对角方阵，根据索引进行取出，构成一个新的矩阵
	def convert_to_one_hot(Y, C):
		Y = np.eye(C)[Y.reshape(-1)].T
		# 或者Y = np.eye(C)[:,Y.reshape(-1)]也是一样的
		return Y
	'''

def softmax_mat2y(A, No):
	'''
	将神经网络输出（矩阵）转换成实验结果的形式（向量）
	Ar: all interger
	'''
	Ar = np.argmax(A, axis=0).reshape(1, -1)
	return Ar

def learning_rate_decay_ip(alpha0, stage, decay_rate=1):
	'''inverse proportional function反比例函数衰减'''
	'''stage是一个阶段性参数，常用的有epoch'''
	alpha_decayed = alpha0 * 1./(1 + decay_rate * stage)
	return alpha_decayed

#@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
# 前向传播过程函数
#@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
def step_forward(A_input, W, B, activation):
	'''
	一个层网络的前向传播，输入需要的参数，输出前向计算参数cache（反向传播需要）、A_output
	'''
	Z = np.dot(W, A_input) + B 		# (Nl, Nl-1) x (Nl-1, M) + (Nl, 1) = (Nl, M)

	if activation == "sigmoid": A_output = sigmoid(Z)
	elif activation == "tanh": A_output = np.tanh(Z)
	elif activation == "relu": A_output = relu(Z)
	elif activation == "leaky_relu": A_output = leaky_relu(Z)
	elif activation == "softmax": A_output = softmax(Z)
	else:
		print("Error, unknown activation function: ", activation, " default: g(z) = z")
		A_output = Z

	forward_cache = (Z,)
	return A_output, forward_cache

#@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
# 反向传播过程函数
#@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
def step_backward(dA_input, forward_cache, A, Apre, W, B, activation):
	'''
	# Apre参与计算dW，W参数预算dA_output
	# 输入需要的参数
	# 输出下一层的dA、本层计算的梯度
	'''
	Z = forward_cache[0]

	if activation == "sigmoid": dZ = dA_input * sigmoid_derivative(Z, A)
	elif activation == "tanh": dZ = dA_input * tanh_derivative(Z, A)
	elif activation == "relu": dZ = dA_input * relu_derivative(Z, A)
	elif activation == "leaky_relu": dZ = dA_input * leaky_relu_derivative(Z, A)
	elif activation == "softmax": dZ = softmax_derivative(Z, dA_input, A)
	elif activation == "dZ": dZ = dA_input # 表示外面已经提前算了dZ
	else:
		print("Error, unknown activation function: ", activation, " default: g(z) = z")
		dZ = dA_input

	M = dZ.shape[1]
	dW = np.dot(dZ, Apre.T) / M
	dB = np.sum(dZ, axis=1, keepdims=True) / M
	dApre_output = np.dot(W.T, dZ)

	gradient = (dW, dB)
	backward_cache = (dZ,)
	return dApre_output, gradient, backward_cache

#@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
# 训练过程函数
#@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
def propagation_forward(networks, X):
	weights = networks['weights']
	Wl = weights['Wl']
	Bl = weights['Bl']
	activations = networks['activations']
	layers = len(activations)

	forward_caches = [0] * layers
	Al = [0] * layers
	Al[0] = X

	for l in range(1, layers, 1):
		Al[l], forward_caches[l] = step_forward(Al[l-1], Wl[l], Bl[l], activations[l])

	return Al, forward_caches

def propagation_backward(networks, dA_input, forward_caches, Al):
	weights = networks['weights']
	Wl = weights['Wl']
	Bl = weights['Bl']
	activations = networks['activations']
	layers = len(activations)

	dAl = [0] * layers
	gradients = [0] * layers
	backward_caches = [0] * layers

	dAl[-1] = dA_input

	for l in range(layers-1, 0, -1):
		dAl[l-1], gradients[l], backward_caches[l] = step_backward(dAl[l], forward_caches[l],\
			Al[l], Al[l-1], Wl[l], Bl[l], activations[l])

	return dAl, gradients, backward_caches

def update_weight(weights, gradients, alpha):
	layers = len(gradients)
	Wl = weights['Wl']
	Bl = weights['Bl']
	for l in range(1, layers, 1):
		Wl[l] -= alpha * gradients[l][0]
		Bl[l] -= alpha * gradients[l][1]

#@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
# 训练、预测接口
#@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
def train(networks, X, Y, alpha, iterations):
	J = [0] * iterations
	activations = networks['activations']

	for i in range(iterations):
		Al, forward_caches = propagation_forward(networks, X)
		J[i] = loss_cross_entropy(Al[-1], Y, activations[-1])
		dA = loss_cross_entropy_derivative(Al[-1], Y, activations[-1])
		dAl, gradients, backward_caches = propagation_backward(networks, dA, forward_caches, Al)

		update_weight(networks['weights'], gradients, alpha)

	return J, dAl[0]

def predict(networks, X):
	Al, forward_caches = propagation_forward(networks, X)
	return Al[-1]

 我的一个使用demo，鸢尾花数据集分类

import dnn_manual as dnn
import numpy as np
import matplotlib.pyplot as plt
import sys
import copy
import dataset_utils as utils

Ni = 4 											# 输入层神经元数, 特征数
No = 3 											# 输出层神经元数
Nl = [Ni, Ni+1, Ni+1, No] 						# 各层含有的神经元个数
activations = ['', 'relu', 'relu', 'softmax']
alpha = 0.1
iterations = 10000
networks = {'Nl': Nl, 'activations': activations, 'weights': None}


'''初始化权重'''
dnn.init_weight(networks)

'''加载数据集'''
Xtrain, Ytrain, Xtest, Ytest = utils.load_iris(".\\iris\\iris_training.csv", ".\\iris\\iris_test.csv")

'''归一化输入'''
Xtrain, normalize_paras = dnn.normalize_universal(Xtrain)

'''训练'''
J, dX = dnn.train(networks, Xtrain, dnn.softmax_y2mat(Ytrain, No),\
	alpha=alpha, iterations=iterations)

'''预测训练集'''
Ytrain_predict = dnn.softmax_mat2y(dnn.predict(networks, Xtrain), No)
accuracy = np.sum(Ytrain == Ytrain_predict)
print('train set accuracy = ', accuracy)

'''预测测试集'''
# 用同样的参数归一化测试集
Xtest, _ = dnn.normalize_universal(Xtest, normalize_paras)
Ytest_predict = dnn.softmax_mat2y(dnn.predict(networks, Xtest), No)
accuracy = np.sum(Ytest == Ytest_predict)
print('test set accuracy = ', accuracy)

plt.plot(range(iterations), J)
plt.show()