为了实践一下CNN运作的内部原理,加深算法的理解,从头开始写了遍CNN。因为刚开始学MATLAB还不太会用,所以拿C写的。写的有点难受,毕竟卷积等运算都得自己手写,所以只是实现了一个结构非常简单的CNN。等有时间了肯定好好学MATLAB~~~
网络框架:
为了方便起见只设置了一层卷积层。原始数据集的图片大小为28*28*1,一轮卷积操作与5个滤波器卷积,步长设置为1,得到24*24*5的卷积层。通过maxpooling后得到大小为12*12*5的池化层。展开后连接一个全连接层做分类器。因为该任务为多分类任务,输出层采用softmax函数,采用交叉熵作为目标函数。
下面罗列一下代码中用到的计算公式:
前向计算:
卷积层输出的feature map为:。
池化层采用maxpooling,在每个2*2的样本中取其最大值。
softmax层的输出为:。
反向传播:
反向传播的原理:反向传播基于链式求导法则实现,我们定义损失函数为,我们的目的是最小化损失函数。
采用梯度下降的法则,则我们要求出损失函数对每个权值的导数,更新的法则为:。
则我们的目标是求得,我们定义某一层的输出为
,敏感项为
,根据链式求导法则,
是
的函数,
是
的函数,所以有
,
。所以从后往前传递敏感项,再根据敏感项求出权值的梯度再更新即可。
我们采取交叉熵作为损失函数为:。
1.敏感项传播:
softmax层的敏感项。
为输出层的输出,
为标准输出。
当前层为全连接层的敏感项传递:
,其中
为隐藏层的输出。
当前层为池化层的敏感项传递:
对于max pooling,下一层的敏感项的值会原封不动的传递到上一层对应区块中的最大值所对应的神经元,而其他神经元的敏感项的值都是0。
当前层为卷积层的敏感项传递:
。因为采用Relu函数作为激活函数,所以输出值不为0的
,输出值为0的
。
2.权重的更新:
对全连接层权重的更新:,其中
为节点
对节点
的输出值。
对全连接层偏置项的更新:。
对filter权重的更新:,
。
对filter偏置项的更新:,
。
因为卷积结构设置不合理(只有一层卷积层)导致模型在测试集上的精确度较低,达到85%左右。
代码实现:
#include <bits/stdc++.h>
using namespace std ;
vector<double>labels;
vector<vector<double> >images;//训练集
vector<double>labels1;
vector<vector<double> >images1;//测试集
const int train_number=20000;//训练样本数
const int test_number=5000;//测试样本数
const int out=10;
const double alpha=0.01;//学习率
int step;
struct layer//定义卷积网络中的层
{
int L,W,H;
double m[30][30][5];
double b[30][30][5];
double delta[30][30][5];
void init()
{
for(int i=0;i<30;i++)
for(int j=0;j<30;j++)
for(int k=0;k<5;k++)
m[i][j][k]=0.01*(rand()%100);
}
};
struct fcnn_layer//定义全连接网络中的层
{
int length;
double m[1000];
double b[1000];
double delta[1000];
double w[20][1000];
void init()
{
for(int i=0;i<20;i++)
for(int j=0;j<1000;j++)
w[i][j]=0.01*(rand()%100);
}
};
struct Network//定义CNN
{
layer Input_layer;
layer conv_layer1;
layer pool_layer1;
layer filter1[5];
fcnn_layer fcnn_input;
fcnn_layer fcnn_w;
fcnn_layer fcnn_outpot;
}CNN;
void init()//权重初始化
{
CNN.Input_layer.init();
CNN.conv_layer1.init();
CNN.pool_layer1.init();
for(int i=0;i<5;i++) CNN.filter1[i].init();
CNN.fcnn_input.init();
CNN.fcnn_w.init();
CNN.fcnn_outpot.init();
}
layer conv(layer A,layer B[],int number,layer C,int flag);//卷积函数,表示卷积层A与number个filterB相卷积
layer CNN_Input(int num,layer A);//将图像输入卷积层
fcnn_layer Classify_input(layer A,fcnn_layer B);//将卷积提取特征输入到全连接神经网络
layer pool_input(layer A,fcnn_layer B);//全连接层的误差项传递到CNN中
layer pool_delta(layer A,layer B);//当前层为池化层的敏感项传递
fcnn_layer softmax(fcnn_layer A);//softmax函数
double Relu(double x);//Relu函数
layer Update(layer A,layer B,layer C);//filter更新
layer maxpooling(layer conv_layer,layer A);//池化前向输出
fcnn_layer fcnn_Mul(fcnn_layer A,fcnn_layer B,fcnn_layer C);//全连接层前向输出
double sum(layer A);//矩阵求和,此处用于敏感项求和
double sum1(layer A,layer B,int x,int y);
double sigmod(double x);
void test();
int test_out(int t);
/**************************此段为读取MNIST数据集模块**************/
int ReverseInt(int i)
{
unsigned char ch1, ch2, ch3, ch4;
ch1 = i & 255;
ch2 = (i >> 8) & 255;
ch3 = (i >> 16) & 255;
ch4 = (i >> 24) & 255;
return((int)ch1 << 24) + ((int)ch2 << 16) + ((int)ch3 << 8) + ch4;
}
void read_Mnist_Label(string filename, vector<double>&labels)
{
ifstream file;
file.open("train-labels.idx1-ubyte", ios::binary);
if (file.is_open())
{
int magic_number = 0;
int number_of_images = 0;
file.read((char*)&magic_number, sizeof(magic_number));
file.read((char*)&number_of_images, sizeof(number_of_images));
magic_number = ReverseInt(magic_number);
number_of_images = ReverseInt(number_of_images);
for (int i = 0; i < number_of_images; i++)
{
unsigned char label = 0;
file.read((char*)&label, sizeof(label));
labels.push_back((double)label);
}
}
}
void read_Mnist_Images(string filename, vector<vector<double> >&images)
{
ifstream file("train-images.idx3-ubyte", ios::binary);
if (file.is_open())
{
int magic_number = 0;
int number_of_images = 0;
int n_rows = 0;
int n_cols = 0;
unsigned char label;
file.read((char*)&magic_number, sizeof(magic_number));
file.read((char*)&number_of_images, sizeof(number_of_images));
file.read((char*)&n_rows, sizeof(n_rows));
file.read((char*)&n_cols, sizeof(n_cols));
magic_number = ReverseInt(magic_number);
number_of_images = ReverseInt(number_of_images);
n_rows = ReverseInt(n_rows);
n_cols = ReverseInt(n_cols);
for (int i = 0; i < number_of_images; i++)
{
vector<double>tp;
for (int r = 0; r < n_rows; r++)
{
for (int c = 0; c < n_cols; c++)
{
unsigned char image = 0;
file.read((char*)&image, sizeof(image));
tp.push_back(image);
}
}
images.push_back(tp);
}
}
}
void read_Mnist_Label1(string filename, vector<double>&labels)
{
ifstream file;
file.open("t10k-labels.idx1-ubyte", ios::binary);
if (file.is_open())
{
int magic_number = 0;
int number_of_images = 0;
file.read((char*)&magic_number, sizeof(magic_number));
file.read((char*)&number_of_images, sizeof(number_of_images));
magic_number = ReverseInt(magic_number);
number_of_images = ReverseInt(number_of_images);
for (int i = 0; i < number_of_images; i++)
{
unsigned char label = 0;
file.read((char*)&label, sizeof(label));
labels.push_back((double)label);
}
}
}
void read_Mnist_Images1(string filename, vector<vector<double> >&images)
{
ifstream file("t10k-images.idx3-ubyte", ios::binary);
if (file.is_open())
{
int magic_number = 0;
int number_of_images = 0;
int n_rows = 0;
int n_cols = 0;
unsigned char label;
file.read((char*)&magic_number, sizeof(magic_number));
file.read((char*)&number_of_images, sizeof(number_of_images));
file.read((char*)&n_rows, sizeof(n_rows));
file.read((char*)&n_cols, sizeof(n_cols));
magic_number = ReverseInt(magic_number);
number_of_images = ReverseInt(number_of_images);
n_rows = ReverseInt(n_rows);
n_cols = ReverseInt(n_cols);
for (int i = 0; i < number_of_images; i++)
{
vector<double>tp;
for (int r = 0; r < n_rows; r++)
{
for (int c = 0; c < n_cols; c++)
{
unsigned char image = 0;
file.read((char*)&image, sizeof(image));
tp.push_back(image);
}
}
images.push_back(tp);
}
}
}
/**************************************************************/
layer conv(layer A,layer B[],int number,layer C)
{
memset(C.m,0,sizeof(C.m));
for(int i=0;i<number;i++)
{
B[i].L=B[i].W=5;
B[i].H=1;
}
C.L=A.L-B[0].L+1;
C.W=A.W-B[0].W+1;
C.H=number;
for(int num=0;num<number;num++)
for(int i=0;i<C.L;i++)
for(int j=0;j<C.W;j++)
{
for(int a=0;a<B[0].L;a++)
for(int b=0;b<B[0].W;b++)
for(int k=0;k<A.H;k++)
C.m[i][j][num]+=A.m[i+a][j+b][k]*B[num].m[a][b][k];
C.m[i][j][num]=Relu(C.m[i][j][num]+C.b[i][j][num]);
}
return C;
}
layer CNN_Input(int num,layer A,int flag)
{
A.L=A.W=28;
A.H=1;
int x=0;
for(int i=0;i<28;i++)
{
for(int j=0;j<28;j++)
{
for(int k=0;k<1;k++)
{
if(flag==0) A.m[i][j][k]=images[num][x];
else if(flag==1) A.m[i][j][k]=images1[num][x];
x++;
}
}
}
return A;
}
fcnn_layer Classify_input(layer A,fcnn_layer B)
{
int x=0;
B.m[x++]=1.0;
for(int i=0;i<A.L;i++)
for(int j=0;j<A.W;j++)
for(int k=0;k<A.H;k++)
B.m[x++]=sigmod(A.m[i][j][k]);
B.length=x;
return B;
}
layer pool_input(layer A,fcnn_layer B)
{
int x=1;
for(int i=0;i<A.L;i++)
for(int j=0;j<A.W;j++)
for(int k=0;k<A.H;k++)
A.delta[i][j][k]=B.delta[x++];
return A;
}
layer pool_delta(layer A,layer B)
{
for(int k=0;k<A.H;k++)
for(int i=0;i<A.L;i+=2)
for(int j=0;j<A.W;j+=2)
{
if(fabs(A.m[i][j][k]-B.m[i/2][j/2][k])<0.01) A.delta[i][j][k]=B.delta[i/2][j/2][k];
else A.delta[i][j][k]=0;
}
return A;
}
layer maxpooling(layer conv_layer,layer A)
{
A.L=conv_layer.L/2;
A.W=conv_layer.W/2;
A.H=conv_layer.H;
for(int k=0;k<conv_layer.H;k++)
for(int i=0;i<conv_layer.L;i+=2)
for(int j=0;j<conv_layer.W;j+=2)
A.m[i/2][j/2][k]=max(max(conv_layer.m[i][j][k],conv_layer.m[i+1][j][k]),
max(conv_layer.m[i][j+1][k],conv_layer.m[i+1][j+1][k]));
return A;
}
fcnn_layer fcnn_Mul(fcnn_layer A,fcnn_layer B,fcnn_layer C)
{
memset(C.m,0,sizeof(C.m));
C.length=out;
for(int i=0;i<C.length;i++)
{
for(int j=0;j<A.length;j++)
{
C.m[i]+=B.w[i][j]*A.m[j];
}
C.m[i]+=B.b[i];
}
return C;
}
double sigmod(double x)
{
return 1.0/(1.0+exp(-x));
}
double sum(layer A)
{
double a=0;
for(int i=0;i<A.L;i++)
for(int j=0;j<A.W;j++)
for(int k=0;k<A.H;k++)
a+=A.delta[i][j][k];
return a;
}
double sum1(layer A,layer B,int x,int y)
{
double a=0;
for(int i=0;i<A.L;i++)
for(int j=0;j<A.W;j++)
for(int k=0;k<A.H;k++)
a+=A.delta[i][j][k]*B.m[i+x][j+y][k];
return a;
}
layer Update(layer A,layer B,layer C)
{
for(int i=0;i<A.L;i++)
for(int j=0;j<A.W;j++)
for(int k=0;k<A.H;k++)
{
A.m[i][j][k]-=alpha*sum1(A,B,i,j);
C.b[i][j][k]-=alpha*sum(A);
}
return A;
}
double Relu(double x)
{
return max(0.0,x);
}
fcnn_layer softmax(fcnn_layer A)
{
double sum=0.0;double maxx=-100000000;
for(int i=0;i<out;i++) maxx=max(maxx,A.m[i]);
for(int i=0;i<out;i++) sum+=exp(A.m[i]-maxx);
for(int i=0;i<out;i++) A.m[i]=exp(A.m[i]-maxx)/sum;
return A;
}
void forward_propagation(int num,int flag)//做一次前向输出
{
CNN.Input_layer=CNN_Input(num,CNN.Input_layer,flag);
CNN.conv_layer1=conv(CNN.Input_layer,CNN.filter1,5,CNN.conv_layer1);
CNN.pool_layer1=maxpooling(CNN.conv_layer1,CNN.pool_layer1);
//CNN.conv_layer2=conv(CNN.pool_layer1,CNN.filter2,3,CNN.conv_layer2,0);
//CNN.pool_layer2=maxpooling(CNN.conv_layer2,CNN.pool_layer2);
CNN.fcnn_input=Classify_input(CNN.pool_layer1,CNN.fcnn_input);
//for(int i=0;i<CNN.fcnn_input.length;i++) printf("%.5f ",CNN.fcnn_input.m[i]);
CNN.fcnn_outpot=fcnn_Mul(CNN.fcnn_input,CNN.fcnn_w,CNN.fcnn_outpot);
CNN.fcnn_outpot=softmax(CNN.fcnn_outpot);
for(int i=0;i<out;i++)
{
if(i==(int)labels[num]) CNN.fcnn_outpot.delta[i]=CNN.fcnn_outpot.m[i]-1.0;
else CNN.fcnn_outpot.delta[i]=CNN.fcnn_outpot.m[i];
//printf("%.5f ",CNN.fcnn_outpot.m[i]);
}
//printf(" %.0f\n",labels[num]);
}
void back_propagation()//反向传播算法
{
memset(CNN.fcnn_input.delta,0,sizeof(CNN.fcnn_input.delta));
for(int i=0;i<CNN.fcnn_input.length;i++)
{
for(int j=0;j<out;j++)
{
CNN.fcnn_input.delta[i]+=CNN.fcnn_input.m[i]*(1.0-CNN.fcnn_input.m[i])*CNN.fcnn_w.w[j][i]*CNN.fcnn_outpot.delta[j];
}
}
for(int i=0;i<CNN.fcnn_input.length;i++)
{
for(int j=0;j<out;j++)
{
CNN.fcnn_w.w[j][i]-=alpha*CNN.fcnn_outpot.delta[j]*CNN.fcnn_input.m[i];
CNN.fcnn_w.b[j]-=alpha*CNN.fcnn_outpot.delta[j];
//printf("%.5f ",CNN.fcnn_w.w[j][i]);
}
}
CNN.pool_layer1=pool_input(CNN.pool_layer1,CNN.fcnn_input);
CNN.conv_layer1=pool_delta(CNN.conv_layer1,CNN.pool_layer1);//pooling误差传递
//CNN.pool_layer1=conv(CNN.conv_layer1,CNN.filter2,3,CNN.pool_layer1,1);//conv误差传递
for(int i=0;i<5;i++) CNN.filter1[i]=Update(CNN.filter1[i],CNN.Input_layer,CNN.conv_layer1);
}
void train()
{
step=0;
for(int time=0;time<100;time++)
{
double err=0;
for(int i=0;i<train_number;i++)
{
forward_propagation(i,0);
err-=log(CNN.fcnn_outpot.m[(int)labels[i]]);
back_propagation();
}
printf("step: %d loss:%.5f\n",time,1.0*err/train_number);//每次记录一遍数据集的平均误差
test();
}
}
void test()
{
int sum=0;
for(int i=0;i<test_number;i++)
{
int ans=test_out(i);
int label=int(labels1[i]);
if(ans==label) sum++;
//printf("%d %d\n",ans,label);
}
//printf("\n");
//for(int i=0;i<out;i++) printf("%.5f ",data_out[i]);
//printf("\n");
printf("step:%d precision: %.5f\n",++step,1.0*sum/test_number);
}
int test_out(int t)
{
forward_propagation(t,1);
int ans=-1;
double sign=-1;
for(int i=0;i<out;i++)
{
if(CNN.fcnn_outpot.m[i]>sign)
{
sign=CNN.fcnn_outpot.m[i];
ans=i;
}
}
return ans;
}
int main()
{
init();
read_Mnist_Label("t10k-labels.idx1-ubyte", labels);
read_Mnist_Images("t10k-images.idx3-ubyte", images);
read_Mnist_Label1("t10k-labels.idx1-ubyte", labels1);
read_Mnist_Images1("t10k-images.idx3-ubyte", images1);//读取mnist数据集
for (int i = 0; i < images1.size(); i++)
{
for (int j = 0; j < images1[0].size(); j++)
{
images1[i][j]/=255.0;
}
}
for (int i = 0; i < images.size(); i++)
{
for (int j = 0; j < images[0].size(); j++)
{
images[i][j]/=255.0;
}
}//图像归一化
train();
return 0;
}
手造的轮子跟框架肯定没得比~~~,所以为了测试更多的性能,采用tensorflow框架也编写了一个CNN。
模型采用两个卷积层和全连接层构成,第一卷积层为32个大小为5*5*1的filter,然后通过2*2的maxpooling下采样。第二个卷积层为64个5*5*32的filter,2*2maxpooling下采样。最后连接到全连接层,通过softmax输出分类。
为了避免模型过拟合,加入dropout。dropout给隐藏层的神经元加上概率为keep_prob的失活率,从而同时训练出指数规模个共享权值的子网络用于分类,增加模型的鲁棒性。
采用随机批梯度下降法,取每100个训练数据的均值来迭代更新一次权值,有利于模型更好的收敛到最小值。
# -*- coding: utf-8 -*-
"""
Created on Wed Mar 6 15:04:23 2019
@author: yexiaohan
"""
import tensorflow as tf
from tensorflow.examples.tutorials.mnist import input_data
#权值初始化
def weight_value(shape):
init=tf.truncated_normal(shape,mean=0.0,stddev=0.1)
return tf.Variable(init)
#偏置初始化
def biase_value(shape):
init=tf.constant(0.1,shape=shape)
return tf.Variable(init)
#卷积
def conv2d(Inputs,Weight):
return tf.nn.conv2d(Inputs,Weight,strides=[1,1,1,1],padding='SAME')
#池化
def maxpooling(Inputs):
return tf.nn.max_pool(Inputs,ksize=[1,2,2,1],strides=[1,2,2,1],padding='SAME')
mnist=input_data.read_data_sets("MNIST_data",one_hot=True)
xs=tf.placeholder(tf.float32,[None,784])#读入图片
ys=tf.placeholder(tf.float32,[None,10])#读入标签
keep_prob=tf.placeholder(tf.float32)#dropout的失活率
Image=tf.reshape(xs,[-1,28,28,1])
#第一个卷积层的定义
conv1_W=weight_value([5,5,1,32])
conv1_b=biase_value([32])
conv1_h=tf.nn.relu(conv2d(Image,conv1_W)+conv1_b)
conv1_out=maxpooling(conv1_h)
#第二个卷积层的定义
conv2_W=weight_value([5,5,32,64])
conv2_b=biase_value([64])
conv2_h=tf.nn.relu(conv2d(conv1_out,conv2_W)+conv2_b)
conv2_out=maxpooling(conv2_h)
#全连接层的定义
fcnn_in=tf.reshape(conv2_out,[-1,49*64])
fcnn1_W=weight_value([49*64,49*64])
fcnn1_b=biase_value([49*64])
fcnn1_out=tf.nn.relu(tf.matmul(fcnn_in,fcnn1_W)+fcnn1_b)
fcnn1_dropout=tf.nn.dropout(fcnn1_out,keep_prob)
fcnn2_W=weight_value([49*64,10])
fcnn2_b=biase_value([10])
prediction=tf.nn.softmax(tf.matmul(fcnn1_dropout,fcnn2_W)+fcnn2_b)
#采用交叉熵做目标函数
cross_entropy=-tf.reduce_sum(ys*tf.log(prediction))
#精确度计算
num=tf.equal(tf.argmax(prediction,1),tf.argmax(ys,1))
accurate=tf.reduce_mean(tf.cast(num,tf.float32))
#训练模型
train_step=tf.train.AdamOptimizer(1e-4).minimize(cross_entropy)
sess=tf.Session()
sess.run(tf.initialize_all_variables())
for step in range(2000):
batch_x,batch_y=mnist.train.next_batch(100)
sess.run(train_step,feed_dict={xs:batch_x,ys:batch_y,keep_prob:0.5})
if step%50==0:
#没迭代50步输出在测试机上的精确度
print(sess.run(accurate,feed_dict={xs:mnist.test.images,ys:mnist.test.labels,keep_prob:1.0}))
用tensorflow搭建的CNN最终在测试集上达到98%以上的准确率。