动手实现DNN,BP算法

：手动实现DNN网络

---

问题描述

本文以机器学习领域一个经典的数据集——Iris为例，手动实现一个4层的神经网络，以根据输入特征预测植物的种类。具体包括：数据的打乱和加载、初始化模型参数、前行计算、损失函数、后向计算、参数更新。在开始训练前以pytorch框架搭建相同的网络以验证手动过程是否正确。

---

以下是正文内容

一、搭建网络，实现各函数

1.导入python库，读入数据

需要的库如下：

from sklearn.datasets import load_iris
import random
import numpy as np
import torch
from torch import nn
from torch.nn import functional as F
from matplotlib import pyplot as plt

读入iris数据，分别得到特征和标签：

#下载数据集
iris=load_iris()

#获取feature数据
features=torch.tensor(iris.data,dtype=torch.float32) 

#获取label数据，注意通过sklearn下载的数据已经经过了预处理
#label里0、1、2分别代表三种类别
labels=torch.tensor(iris.target,dtype=torch.long)   

我们来看一下features和labels中的数据：

print('前三个样本特征:\n'+str(features[0:3]))
print('前三个样本标签:\n'+str(labels[0:3]))

output:
前三个样本特征:
tensor([[5.1000, 3.5000, 1.4000, 0.2000],
        [4.9000, 3.0000, 1.4000, 0.2000],
        [4.7000, 3.2000, 1.3000, 0.2000]])
前三个样本标签:
tensor([0, 0, 0])

---

2.定义数据加载函数

def data_iter(batch_size,which_iter='train_iter'):
    '''
    load_data函数第一个参数是每次用于计算的样本批量数目
    第二个参数指明用于训练集、验证集还是测试集
    '''
    
    all_examples=len(features)          #all_examples是所有样本数目，本例中为150
    iter_examples=0                     #iter_examples是每个训练集拥有的样本数
    
    indices=list(range(all_examples))   #生成一个列表[0:150]，起到下标的作用
    
    r=random.random                     #这里将上述下标列表打散
    random.seed(2)                      #指定随机种子，以复现试验结果
    random.shuffle(indices,r)
    
    
    #接下来划分数据集
    #我们规定训练集、验证集、训练集的样本分别取80%，10%，10%
    
    if which_iter=='train_iter':
        for i in range(0,120,batch_size):
            j=torch.tensor(indices[i:min(i+batch_size,120)])
            yield(torch.index_select(features,0,j),torch.index_select(labels,0,j))
        
    elif which_iter=='validation_iter':
        for i in range(120,135,batch_size):
            j=torch.tensor(indices[i:min(i+batch_size,135)])
            yield(torch.index_select(features,0,j),torch.index_select(labels,0,j))
        
    elif which_iter=='test_iter':
        for i in range(135,150,batch_size):
            j=torch.tensor(indices[i:min(i+batch_size,150)])
            yield(torch.index_select(features,0,j),torch.index_select(labels,0,j))

来测试一下data_iter

train_iter1=data_iter(3,'train_iter')
for x,y in train_iter1:
    print(x)
    print(y)
    break
    
output:
tensor([[6.7000, 2.5000, 5.8000, 1.8000],
        [7.0000, 3.2000, 4.7000, 1.4000],
        [5.0000, 3.4000, 1.5000, 0.2000]])
tensor([2, 1, 0])

嗯，没有问题

---

3.定义参数初始化函数

def initialize_parameters(layers_dims,seed=3):
    '''
    本函数用于初始化一个深层网络的参数
    layers_dims是一个列表，包含每层的节点数目
    在本实验中，layers_dims是一个四维向量，且layers_dims[0]=4,layers_dims[3]=3
    
    本函数的输出是一个网络的pearmeters，为字典类型
    字典的key分别是w1,b1,w2,b2...
    '''
    
    assert(len(layers_dims)==4)      #本次实验规定了层数为4层
    assert(layers_dims[0]==4)        #输入向量4维
    assert(layers_dims[3]==3)        #输出向量3维
    

    torch.manual_seed(seed)  #指定随机种子以复现结果
    
    parameters={}
    
    #本人习惯y=xw+b的形式这决定了权重参数的形状
    parameters['w1']=torch.randn(layers_dims[0],layers_dims[1],dtype=torch.float32)
    parameters['b1']=torch.zeros(1,layers_dims[1],dtype=torch.float32)
    
    parameters['w2']=torch.randn(layers_dims[1],layers_dims[2],dtype=torch.float32)
    parameters['b2']=torch.zeros(1,layers_dims[2],dtype=torch.float32)
    
    parameters['w3']=torch.randn(layers_dims[2],layers_dims[3],dtype=torch.float32)
    parameters['b3']=torch.zeros(1,layers_dims[3],dtype=torch.float32)
    
    return parameters        

下面测试一下initialize_parameters

layers_dims=[4,4,4,3]
parameters=initialize_parameters(layers_dims)

for i in range(1,len(layers_dims)):
    print('w'+str(i)+'.shape='+str(parameters['w'+str(i)].size()))
print('\nw1=\n'+str(parameters['w1']))
print('\nb1=\n'+str(parameters['b1']))

output:
w1.shape=torch.Size([4, 4])
w2.shape=torch.Size([4, 4])
w3.shape=torch.Size([4, 3])

w1=
tensor([[-0.8173, -0.5556, -0.8267, -1.2970],
        [-0.1974, -0.9643, -0.5133,  2.6278],
        [-0.7465,  1.0051, -0.2568,  0.4765],
        [-0.6652, -0.3627, -1.4504, -0.2496]])

b1=
tensor([[0., 0., 0., 0.]])

嗯，也不错！

---

4.定义损失函数

本实验由于是多分类问题，故输出单元采用SoftMax输出单元； y_hat是模型输出的概率质量函数，是由softmax输出单元得到的； y_hat=exp(z)/(sum(exp(z)))； y是样本的标签，取值为0或1或2；

def cost_function(y_hat,y):
    '''
    最大似然估计往往是大多数深度神经网络的最佳损失函数，负对数似然公式如下
    '''
    cost=-torch.log(torch.gather(y_hat,1,y.unsqueeze(1))) 
    #交叉墒损失等价于负对数似然，即样本的标签值作为模型概率质量函数的角标
    cost=torch.sum(cost)  
    
    return cost

---

5.定义sigmoid激活函数及其导数

def sigmoid(z):
    return 1/(1+torch.exp(-z))
def diff_sigmoid(z):
    return (1-sigmoid(z))*sigmoid(z)

---

6.定义前行计算

定义前向计算，输入一个四维向量x，输出一个字典cache，包含各个前向计算的中间值，以用于后向计算 我所采用的形式为y=xw+b 不同于pytorch的线性层采用的y=xw.T+b的形式

def forward_propagation(x,parameters,activation='sigmoid'):
    '''
    parameters是一个字典，包含用于前向计算的所有权重和偏差参数  
    '''
    
    cache={}
    cache['A0']=x
    
    cache['z1']=torch.matmul(x,parameters['w1'])+parameters['b1']
    cache['A1']=sigmoid(cache['z1'])
    cache['z2']=torch.matmul(cache['A1'],parameters['w2'])+parameters['b2']
    cache['A2']=sigmoid(cache['z2'])
    cache['z3']=torch.matmul(cache['A2'],parameters['w3'])+parameters['b3']
    
    #softmax输出单元通过对z指数化和归一化，得到有效的概率质量函数y_hat
    y_hat=torch.exp(cache['z3'])/(torch.sum(torch.exp(cache['z3']),1,keepdim=True))
    
    cache['y_hat']=y_hat
    
    return cache

---

下面测试一下forward_propagation

train_iter=data_iter(batch_size=3)
layers_dims=[4,4,4,3]
parameters=initialize_parameters(layers_dims)
for x,y in train_iter:
    break

cache=forward_propagation(x,parameters)
cost=cost_function(cache['y_hat'],y)
print('概率质量函数:\n'+str(cache['y_hat']))
print('损失:\n'+str(cost))

output:
概率质量函数:
tensor([[0.0597, 0.4179, 0.5224],
        [0.0649, 0.4120, 0.5231],
        [0.0698, 0.4206, 0.5096]])
损失:
tensor(4.1980)

不错嗷

7.定义后向计算

def backward_propagation(cache,parameters,y):
    '''
    定义后向计算函数，输入模型参数，样本标签，和缓存量
    输出一个字典，包含梯度
    '''
    
    y_hat=cache['y_hat']
    batch_size=cache['y_hat'].size()[0]
    
    A0,A1,A2=cache['A0'],cache['A1'],cache['A2']  #A0即是模型的输入，即x
    z1,z2=cache['z1'],cache['z2']
    w2,w3=parameters['w2'],parameters['w3']
    
    dz3=torch.zeros(y_hat.size())
    dz3.copy_(y_hat)
    for i in range(batch_size):
        dz3[i][y[i]]-=1
    
    
    dw3=torch.matmul(A2.transpose(0,1),dz3)/batch_size   #w3的导数
    
    dz2=torch.matmul(dz3,w3.transpose(0,1))*diff_sigmoid(z2)
    
    dw2=torch.matmul(A1.transpose(0,1),dz2)/batch_size   #w2的导数
    
    dz1=torch.matmul(dz2,w2.transpose(0,1))*diff_sigmoid(z1)
    
    dw1=torch.matmul(A0.transpose(0,1),dz1)/batch_size   #w1的导数
    
    db3=torch.sum(dz3,0,keepdim=True)/batch_size                    #b3的导数
    
    db2=torch.sum(dz2,0,keepdim=True)/batch_size                    #b2的导数
    
    db1=torch.sum(dz1,0,keepdim=True)/batch_size                    #b1的导数
    
    grad={'dw1':dw1,'dw2':dw2,'dw3':dw3,'db1':db1,'db2':db2,'db3':db3}
    return grad

---

8.用pytorch搭建神经网络

为了验证手动实现的梯度运算是否正确，这里用框架搭建同样参数的神经网络，然后自动求导以验证

class MyModule(nn.Module):
    def __init__(self,n_hidden1,n_hidden2,parameters):
        '''
        由于输入层维度为3，输出层维度为4，是定值
        故这里n_hidden1和n_hidden_2指定的是中间两层的维度
        
        parameters提供与手动BP算法同样的参数
        '''
        super(MyModule,self).__init__()
        
        self.hidden1=nn.Linear(4,n_hidden1)
        self.hidden2=nn.Linear(n_hidden1,n_hidden2)
        self.hidden3=nn.Linear(n_hidden2,3)  
        
		#为了得到相同的梯度，两个模型应该共用一批参数
        #前文已经提到，我采用的 y=xw+b 的形式
        #而pytorch线性层采用  y=xw.T+b 的形式
        #故对权重有一个转置操作
        self.hidden1.weight.data=parameters['w1'].transpose(0,1)
        self.hidden1.bias.data=parameters['b1']
        self.hidden2.weight.data=parameters['w2'].transpose(0,1)
        self.hidden2.bias.data=parameters['b2']
        self.hidden3.weight.data=parameters['w3'].transpose(0,1)
        self.hidden3.bias.data=parameters['b3']
        
    def forward(self,x):
        
        assert(x.dtype==torch.float32)
        
        x=torch.sigmoid(self.hidden1(x))
        x=torch.sigmoid(self.hidden2(x))
        x=self.hidden3(x)
        
        return x

---

二、验证梯度正确性

网络搭建完毕，开始验证手动BP算法的损失函数和梯度是否正确

train_iter=data_iter(4)            #加载数据集
layers_dims=[4,5,4,3]              #指定各层的神经元数目
parameters=initialize_parameters(layers_dims)     #初始化模型参数
for x,y in train_iter:          #取4个样本做验证即可
    break

用pytorch计算损失和梯度:

net=MyModule(5,4,parameters)
z3=net(x)
net_cost=F.cross_entropy(z3,y)
net_cost.backward()

用自定义模型计算损失和梯度

cache=forward_propagation(x,parameters)
cost=cost_function(cache['y_hat'],y)
grad=backward_propagation(cache,parameters,y)

这里grad是由手动模型得到的梯度字典，将grad[‘w1’]与net.hidden1的权重做差，若每个元素都小于1e(-7)，说明二者相等。其他梯度亦然。

print('手动BP网络的损失:'+str(cost))
print('由框架实现的损失:'+str(net_cost))

assert((torch.sum(grad['dw1']-net.hidden1.weight.grad.transpose(0,1)))<10**(-7))
assert((torch.sum(grad['dw2']-net.hidden2.weight.grad.transpose(0,1)))<10**(-7))
assert((torch.sum(grad['dw3']-net.hidden3.weight.grad.transpose(0,1)))<10**(-7))

assert((torch.sum(grad['db1']-net.hidden1.bias.grad))<10**(-7))
assert((torch.sum(grad['db2']-net.hidden2.bias.grad))<10**(-7))
assert((torch.sum(grad['db3']-net.hidden3.bias.grad))<10**(-7))

output:
手动BP网络的损失:tensor(1.2027)
由框架实现的损失:tensor(1.2027, grad_fn=<NllLossBackward>)

可见两种方式实现的网络的损失函数和梯度是相等的，这说明上述手动实现的BP算法是正确的

---

三、开始训练

在开始训练之前，还要定义参数更新函数和评价函数

def update_parameters(parameters,grad,lr):
    for i in range(1,4):
        parameters['w'+str(i)]=parameters['w'+str(i)]-lr*grad['dw'+str(i)]
        parameters['b'+str(i)]=parameters['b'+str(i)]-lr*grad['db'+str(i)]
    return parameters

def evaluate_accuracy(test_set,parameters):
    
    num_correct,len_set=0,0   #num_correct是正确分类的个数
                              #len_set是测试集的长度
    
    for x,y in test_set:
        cache=forward_propagation(x,parameters)
        predict=torch.argmax(cache['y_hat'])
        if y==predict:
            num_correct+=1
        len_set+=1
        
    return num_correct/len_set*100

接下来可以正式训练了

def train_model(batch_size,parameters,num_epochs,lr=0.005):
    cost_on_train=[]
    cost_on_validation=[]
    for n in range(num_epochs):
        train_iter=data_iter(batch_size,'train_iter',seed=7)
        validation_iter=data_iter(15,'validation_iter')
        
        train_cost_cache=[]
        
        for x,y in train_iter:
            cache=forward_propagation(x,parameters)
            train_cost_cache.append(cost_function(cache['y_hat'],y))
            grad=backward_propagation(cache,parameters,y)
            parameters=update_parameters(parameters,grad,lr)
        cost_on_train.append(sum(train_cost_cache)/len(train_cost_cache))
        
        for x,y in validation_iter:
            cache=forward_propagation(x,parameters)
            cost_on_validation.append(cost_function(cache['y_hat'],y))
            
    fig,ax = plt.subplots()
    
    ax.plot(range(num_epochs),cost_on_train,label='average cost on train set')
    ax.plot(range(num_epochs),cost_on_validation,label='average cost on validation set')
    ax.set_xlabel('num_epoch')
    ax.set_ylabel('average cost')
    ax.legend()
    plt.show()
    return parameters,cost_on_train

layers_dims=[4,6,6,3]

parameters=initialize_parameters(layers_dims,seed=3)

parameters,cost_on_train=train_model(1,parameters,num_epochs=300,lr=0.02)

来看看在训练集、验证集和测试集上的准确率：

训练集准确率

train_iter=data_iter(1,'train_iter')
train_accuracy=evaluate_accuracy(train_iter,parameters)
print('训练集上的准确率:'+str(train_accuracy)+'%')

output:
训练集上的准确率:98.33333333333333%

验证集准确率

validation_iter=data_iter(1,'validation_iter')
validation_accuracy=evaluate_accuracy(validation_iter,parameters)
print('验证集上的准确率:'+str(validation_accuracy)+'%')

output:
验证集上的准确率:93.33333333333333%

测试集准确率

test_iter=data_iter(1,'test_iter')
test_accuracy=evaluate_accuracy(test_iter,parameters)
print('测试集上的准确率:'+str(test_accuracy)+'%')

output:
测试集上的准确率:100.0%

---