上回文说到基于误差梯度下降的BP网络算法容易陷入局部极小,通常的改进方法先使用遗传算法生成比较好的权重值,再交给神经网络训练。
遗传算法随着进化的进行,其选择率、交叉算子、变异率应该是动态改变的。
编码方式
在使用BP网络进行文本分类时,大都是采用实数编码,把权值设为[0,1]上的实数,这是因为要使用权值调整公式要求权值是实数。但是在使用遗传算法优化这些权值时,完全可以把它们编码为整数。比如设为[1,64]上的整数,一个权值只有64种选择,而[0,1]上的实数有无穷多个,这样既可以缩小搜寻的范围,同时也加大了搜寻的步长。毕竟BP网络中很多个极小点,使用遗传的目的只是在全局找个一个比较优的解,进一步的精确寻优交给BP神经网络来做。
选择算子
在进化初期我们应该使用较小的选择压力,以鼓励种群向着多样化发展;在进化后期个体差异不大,适应度都很高,这时应增大选择压力以刺激进化速度。可以使用模拟退火(SA,simulated annealing)来决定选择率,即我们以一定的概率来接收不好的个体:
这是模拟退火的原始表达式,意思是说在金属退火的过程中,其能量在降低(
<0),我们以
的概率接收本次变化,显然当温度T越低时,接收概率越大,温度T越高时,接收概率越小,k是常数。对应到遗传算法,就是当种群平均适应值越低时,接收劣等个体的概率越高,当种群平均适应值越高时,接收劣等个体的概率越小。
另外M.Srinivas提出当群体适应度比较集中时,使交叉概率增大;当群体适应度比较分散时,使交叉概率减小。种群适应度分散与否通过最大、最小和平均适应度来衡量。
选择算子是保证遗传算法能找到近优解的唯一手段,当染色体唯度很高时,遗传算法很难找到较好的解。这是因为最开始生成的初始种群适应度都极其的低,个体之间(适应度)差异不大,如果使用锦标赛选择法则跟随机选择无异,即使使用赌轮法选择到最优个体的概率会大增加,但是最优个体也不比最劣个体好到哪儿去,最优个体也不含有优良的基因片段。所以对于高维数据,在进化初期主要靠交叉进行全局搜索来搜寻较优的个体。
交叉算子
交叉实际上就是在进行全局搜索,所以遗传算法不过是穷举算法的一个变种。在进化初期,种群多样性高,采用单点交叉就可以获得较广的搜索空间;在进化后期,个体差异不大,需要采用多点交叉,或者有人采用变异交叉点的方法。
当发现种群的适应度操持不变时,可能已经进入了局部最优,应该变异交叉点,大步跨出当前的小山峰。
由于要保留精英个体,所以交叉要以一定的概率进行。随着进化的进行,交叉率应逐渐降低趋于某个值,以避免算法不收敛。
变异算子
直观上对好的个体应施以较小的变异率,对劣等个体应施以较大的变异率。
当发现种群的适应度在降低时应增大变异率。
另外M.Srinivas提出当群体适应度比较集中时,使变异概率增大;当群体适应度比较分散时,使变异概率减小。种群适应度分散与否通过最大、最小和平均适应度来衡量。
下面的代码是用遗传算法来为BP网络寻找比较好的初始解。但是遗传算法根本就没有起到作用,因为我的神经网络输入层1000个节点,隐藏层20个节点,这就2万个权值了,也就是说染色体的长度在2万维以上,用遗传算法根本就找不到较优解,它始终是在随机地遍历,一点儿没有想“进化”的意思。
#include<iostream>
#include<cmath>
#include<cstdlib>
#include<ctime>
#include<cassert>
#include<vector>
#include<fstream>
#include<sstream>
#include<limits>
#include<algorithm>
using namespace std;
inline double sigmoid(double activation,double response);
void initIO();
void initWeight();
void printWeight();
void getOutput(vector<double> &input,vector<double> &Y,vector<double> &output);
void adjustWeight(vector<vector<double> > &input,vector<vector<double> > &Y,vector<vector<double> > &output,vector<vector<double> > &D);
double getMSE(vector<vector<double> > &output,vector<vector<double> > &D);
void initPopulation();
void select();
template <typename Comparable>void InsertSort(vector<Comparable> &vec);
template<typename T>void multipointcross(vector<T> &vec1,vector<T> &vec2,int p);
void cross();
void mutation();
void runGA();
/*************************************/
/** 语料库信息 **/
/*************************************/
const int dim=1000; //样本向量的维度
const int catnum=9; //9大分类
const string classes[]={"C3","C7","C11","C19","C31","C32","C34","C38","C39"}; //类别名称(训练集和测试集都按照这个顺序)
const int test_total=7193; //测试样本总数
const int testnum[catnum]={508,447,463,1066,741,822,1404,758,1066}; //测试样本每个类别的文本数
/*************************************/
/** BP网络参数 **/
/*************************************/
const int train_per_cat=40; //每个类别训练文本的个数
const int P=catnum*train_per_cat; //总的训练样本数
const int in_count=dim+1; //输入层节点数
const int hidden_count=20; //隐藏层节点数
const int out_count=catnum; //输出层节点数等于类别数
const int iter_lim=200; //最大迭代次数
const double Epsilon=10; //允许误差
const double flat=1; //判断误差曲面是否进入平坦区域的依据
double Eta=0.2; //学习率
const double beta=0.8; //调整无效时减小学习率
const double theta=1.2; //调整有效时增大学习率
double lambda=1; //陡度因子
double W[hidden_count][out_count]={0}; //从隐藏层到输出层的权值
double V[in_count][hidden_count-1]={0}; //从输入层到隐藏层的权值
double alpha=0.7; //动量系数
double pre_delte_W[hidden_count][out_count]={0}; //上一次的权值调整量
double pre_delte_V[in_count][hidden_count]={0}; //数组必须显式地赋0,否则它的初始值是一个随机的数
double pre_E=RAND_MAX; //上一次的系统误差
vector<vector<double> > X; //神经网络的输入
vector<vector<double> > D; //网络的期望输出
/*************************************/
/** 遗传算法参数 **/
/*************************************/
const int genomelen=in_count*(hidden_count-1)+hidden_count*out_count; //基因长度
//基因结构体
struct genome{
vector<double> gene;
double fitness;
genome(){
gene.assign(genomelen,0);
}
bool operator < (const genome & rhs) const{
return fitness<rhs.fitness;
}
//计算基因适应度
double calFitness(){
//先创建两个变量Y和O,后面要用
vector<vector<double> > Y;
vector<vector<double> > O;
Y.reserve(P);
O.reserve(P);
for(int i=0;i<P;++i){
vector<double> y(hidden_count);
y[0]=-1;
Y.push_back(y);
vector<double> o(out_count);
O.push_back(o);
}
//把基因赋给神经网络的V和W
int index=0;
for(int i=0;i<in_count;++i){
for(int j=0;j<hidden_count-1;++j){
V[i][j]=gene[index++];
}
}
for(int i=0;i<hidden_count;++i){
for(int j=0;j<out_count;++j){
W[i][j]=gene[index++];
}
}
//计算每个样本经过网络后的输出
for(int p=0;p<P;++p)
getOutput(X.at(p),Y.at(p),O.at(p));
//计算一批样本的均方误差
double err=getMSE(O,D);
//将误差转化为适应度
fitness=pow(M_E,(0-err)/1000);
return fitness;
}
};
const int population_size=20; //种群大小
vector<genome> population(population_size); //当代种群
genome best_individual; //保存全局最优个体
double pre_maxfit=0; //上一代的最高适应度
double pre_avgfit=0; //上一代平均适应值
double maxfit; //上一代的平均适应度
double avgfit; //当代的平均适应度
int ga_iteration_limit=100; //GA最大迭代次数
double pc=0.8; //交叉率
double pm=0.05; //变异率
/**
* 初始化种群
*/
void initPopulation(){
srand(time(0));
double avg=0.0;
for(int i=0;i<population_size;++i){
genome g;
for(int j=0;j<genomelen;++j){
g.gene[j]=rand()/(double)RAND_MAX;
}
double ff=g.calFitness();
if(ff>maxfit)
maxfit=ff;
avg+=ff;
population[i]=g;
}
avg/=population_size;
avgfit=avg;
}
/**
* 选择算子
*/
void select(){
//对适应度按照模拟退火进行缩放
vector<pair<genome,double> > vec;
//cout<<"选择前的适应度:";
for(int i=0;i<population_size;i++){
double f=population[i].fitness;
//double sf=pow(M_E,(f-maxfit)/(0.5*avgfit));
double sf=f;
//cout<<sf<<"\t";
pair<genome,double> p(population[i],sf);
vec.push_back(p);
}
//cout<<endl;
//对sf进行累加
for(int i=1;i<population_size;i++){
vec[i].second+=vec[i-1].second;
}
//生成随机数进行选择
vector<genome> newpopulation;
//cout<<"选择后的适应度:";
for(int i=0;i<population_size;i++){
double r=rand()/(double)RAND_MAX*vec[population_size-1].second;
//逆着找到第一个比r小
int j=population_size-1;
for(;j>=0;){
if(vec[j].second<r)
break;
--j;
}
//cout<<vec[j+1].first.fitness<<"\t";
newpopulation.push_back(vec[j+1].first);
}
//cout<<endl;
//用新种群替换旧种群
population=newpopulation;
}
/*插入排序*/
template <typename Comparable>
void InsertSort(vector<Comparable> &vec){
for(int i=1;i<vec.size();i++){
if(vec[i]<vec[i-1]){
int num=vec[i];
int k=i-1;
for(;k>=0;k--){
if(num<vec[k]){
vec[k+1]=vec[k];
}
else
break;
}
vec[k+1]=num;
}
}
}
/**
* 对两个基因进行p点交叉
*/
template<typename T>
void multipointcross(vector<T> &vec1,vector<T> &vec2,int p){
assert(vec1.size()==vec2.size());
//随机生成p个断开点
vector<int> points;
while(points.size()<p){
int point=rand()%vec1.size();
if(point==0) continue;
vector<int>::iterator itr=find(points.begin(),points.end(),point);
if(itr==points.end())
points.push_back(point);
}
InsertSort(points); //对断开的位置进行排序
int i=0;
for(;i<p-1;i+=2){
vector<T> tmp;
tmp.assign(vec1.begin()+points[i],vec1.begin()+points[i+1]);
int index=0;
for(int j=points[i];j<points[i+1];++j){
vec1[j]=vec2[j];
vec2[j]=tmp[index++];
}
}
vector<T> tmp;
tmp.assign(vec1.begin()+points[i],vec1.end());
int index=0;
for(int j=points[i];j<vec1.size();++j){
vec1[j]=vec2[j];
vec2[j]=tmp[index++];
}
}
/**
* 交叉算子
*/
void cross(){
int slicep=1; //一般采用单点交叉
//如果平均适应度没有达到0.8,前后两代平均适应度变化不大,且本代个体的适应度都比较集中,则采用3点交叉,跳出局部最优
if(fabs(avgfit-pre_avgfit)<0.05 && fabs(maxfit-avgfit)<0.1 && avgfit<0.8)
slicep=3;
for(int i=0;i<population_size-2;i+=2){
if(rand()/(double)RAND_MAX<pc){
multipointcross(population[i].gene,population[i+1].gene,slicep);
}
}
}
/**
* 变异算子
*/
void mutation(){
pm/=(maxfit-avgfit); //当种群适应度比较集中时,选择较大的变异率
for(int i=0;i<population_size;++i){
int count=population_size/100; //在基因长度的百分之一个位置点上进行变异
while(count-->0){
if(rand()/(double)RAND_MAX<pm){
int position=rand()%population[i].gene.size(); //随机选择一个变异位置
population[i].gene[position]=rand()/(double)RAND_MAX; //变异就是重新赋它一个随机数
}
}
}
}
/**
* 运行遗传算法
*/
void runGA(){
int itr=ga_iteration_limit;
while(itr-->0){
cout<<"遗传第"<<ga_iteration_limit-itr<<"代,最高适应度:";
initPopulation();
select();
cross();
mutation();
pre_avgfit=avgfit;
pre_maxfit=maxfit;
for(int i=0;i<population_size;i++)
population[i].calFitness();
sort(population.begin(),population.end()); //先对当代种群按适应度从小到大排序
maxfit=population[population_size-1].fitness; //计算当代的最大适应度
cout<<maxfit;
if(maxfit>best_individual.fitness) //检查是否是全局最优个体
best_individual=population[population_size-1];
if(best_individual.fitness>0.8){
cout<<"遗传算法得到的最优个体适应度已超过0.8,遗传退出。"<<endl;
break;
}
//计算当代的平均适应度
for(int i=0;i<population_size;i++)
avgfit+=population[i].fitness;
avgfit/=population_size;
cout<<",平均适应度:"<<avgfit<<endl;
}
cout<<"遗传算法达到最大的迭代次数,得到的最优个体适应度为"<<best_individual.fitness<<",遗传退出。"<<endl;
//对全局最优个体进行解码
int index=0;
for(int i=0;i<in_count;++i){
for(int j=0;j<hidden_count-1;++j){
V[i][j]=best_individual.gene[index++];
}
}
for(int i=0;i<hidden_count;++i){
for(int j=0;j<out_count;++j){
W[i][j]=best_individual.gene[index++];
}
}
}
/**
* 单极性Sigmoid函数
*/
inline double sigmoid(double activation,double response){
double ex=-activation/response;
return 1.0/(pow(M_E,ex)+1);
}
/**
* 初始化网络的输入和期望输出
*/
void initIO(){
X.reserve(P);
D.reserve(P);
ifstream infile("/home/orisun/master/fudan_corpus/tc_ann.txt");
for(int i=0;i<catnum;++i){
for(int j=0;j<train_per_cat;++j){
vector<double> in;
in.reserve(in_count);
in.push_back(-1);
string fn;
getline(infile,fn);
ifstream input(fn.c_str());
if(!input) //打开文件失败
break;
string s;
while(input>>s){
in.push_back(atof(s.c_str()));
}
//cout<<in.size()<<endl;
assert(in.size()==in_count);
input.close();
X.push_back(in);
vector<double> d(out_count);
d[i]=1;
D.push_back(d);
}
}
infile.close();
}
/**
* 初始网络权值W和V,赋予[0,1]上的随机数
*/
void initWeight(){
srand(time(0));
for(int i=0;i<hidden_count;++i){
for(int j=0;j<out_count;++j)
W[i][j]=rand()/(double)RAND_MAX;
}
for(int i=0;i<in_count;++i){
for(int j=0;j<hidden_count-1;++j)
V[i][j]=rand()/(double)RAND_MAX;
}
}
/**
* 打印权重
*/
void printWeight(){
cout<<"W="<<endl;
for(int i=0;i<hidden_count;++i){
for(int j=0;j<out_count;++j)
cout<<W[i][j]<<"\t";
cout<<endl;
}
cout<<"V="<<endl;
for(int i=0;i<in_count;++i){
for(int j=0;j<hidden_count-1;++j)
cout<<V[i][j]<<"\t";
cout<<endl;
}
}
/**
* 给定输入,求网络的输出
*/
void getOutput(vector<double> &input,vector<double> &Y,vector<double> &output){
//cout<<input.size()<<endl;
//cout<<Y.size()<<endl;
//cout<<output.size()<<endl;
assert(input.size()==in_count && Y.size()==hidden_count && output.size()==out_count);
assert(input[0]==-1);
assert(Y[0]==-1);
for(int j=1;j<hidden_count;++j){
double net=0.0; //隐藏层的净输入
for(int i=0;i<in_count;++i)
net+=input[i]*V[i][j];
Y[j]=sigmoid(net,lambda); //把净输入抛给S形函数,得到隐藏层的输出
}
for(int k=0;k<out_count;++k){
double net=0.0; //输出层的净输入
for(int j=0;j<hidden_count;++j)
net+=Y[j]*W[j][k];
output[k]=sigmoid(net,lambda); //把净输入抛给S形函数,得到输出层的输出
//cout<<output[k]<<"\t";
}
//cout<<endl;
}
/**
* 批训练法根据样本总体误差调整权重W和V
*/
void adjustWeight(vector<vector<double> > &input,vector<vector<double> > &Y,
vector<vector<double> > &output,vector<vector<double> > &D){
double delte_W[hidden_count][out_count]={0}; //数组必须显式地赋0,否则它的初始值是一个随机的数
double delte_V[in_count][hidden_count]={0};
for(int j=0;j<hidden_count;++j){
for(int k=0;k<out_count;++k){
for(int p=0;p<P;++p){
delte_W[j][k]+=(D[p][k]-output[p][k])*output[p][k]*(1-output[p][k])*Y[p][j];
}
delte_W[j][k]*=Eta;
W[j][k]+=delte_W[j][k]+alpha*pre_delte_W[j][k]; //加入动量项
pre_delte_W[j][k]=delte_W[j][k];
}
}
for(int i=0;i<in_count;++i){
for(int j=0;j<hidden_count;++j){
for(int p=0;p<P;++p){
double tmp=0.0;
for(int k=0;k<out_count;++k){
tmp+=(D[p][k]-output[p][k])*output[p][k]*(1-output[p][k])*W[j][k];
}
delte_V[i][j]+=tmp*Y[p][j]*(1-Y[p][j])*input[p][i];
}
delte_V[i][j]*=Eta;
V[i][j]+=delte_V[i][j]+alpha*pre_delte_V[i][j]; //加入动量项
pre_delte_V[i][j]=delte_V[i][j];
}
}
}
/**
* 计算所有样本的均方误差和
*/
double getMSE(vector<vector<double> > &output,vector<vector<double> > &D){
double error=0.0;
for(int p=0;p<P;++p){
for(int k=0;k<out_count;++k){
error+=pow((D[p][k]-output[p][k]),2);
}
}
error/=2;
return error;
}
int main(){
initIO();
initWeight();
runGA();
vector<vector<double> > Y;
vector<vector<double> > O;
Y.reserve(P);
O.reserve(P);
for(int i=0;i<P;++i){
vector<double> y(hidden_count);
y[0]=-1;
Y.push_back(y);
vector<double> o(out_count);
O.push_back(o);
}
int iteration=iter_lim;
//printWeight();
while(iteration-->0){
for(int p=0;p<P;++p)
getOutput(X.at(p),Y.at(p),O.at(p));
double err=getMSE(O,D);
cout<<"第"<<iter_lim-1-iteration<<"次迭代误差:"<<err<<endl;
//printWeight();
if(err<Epsilon){ //如果误差小于允许的误差,则退出迭代
cout<<"误差小于允许误差,迭代退出。"<<endl;
break;
}
//动态调整学习率
if(err>pre_E) //误差比上次大,减小学习率
Eta*=beta;
else if(err<pre_E) //误差比上次小,增大学习率
Eta*=theta;
pre_E=err;
//动态调整陡度因子
if(err-pre_E<flat && pre_E-err<flat && err>5*Epsilon){ //误差变化量接近于0(进入平坦区),而误差仍很大
lambda*=theta;
}
else if(err-pre_E>flat || pre_E-err>flat){ //退出平坦区后
lambda*=beta;
}
adjustWeight(X,Y,O,D);
}
//使用原样本进行测试
vector<double> Out(out_count);
int errnum=0;
for(int p=0;p<P;++p){
getOutput(X[p],Y[p],Out);
int max_index=0;
double max=numeric_limits<double>::min();
for(int k=0;k<out_count;k++){
if(Out.at(k)>max){
max=Out.at(k);
max_index=k;
}
}
if(max_index!=p/train_per_cat){
errnum++;
cout<<p/train_per_cat<<"--->"<<max_index<<endl;
}
}
cout<<"错误总数:"<<errnum<<endl;
return 0;
}