自适应遗传算法

上回文说到基于误差梯度下降的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;
}

  

posted @ 2011-12-02 11:35  张朝阳  阅读(8721)  评论(1编辑  收藏  举报