单纯形算法.txt

求多维函数极值的一种算法

求多维函数极值的一种算法，由Nelder和Mead提出，又叫单纯形算法，但和线性规划中的单纯形算法是不同的，由于未利用任何求导运算，算法比较简单，但收敛速度较慢，适合变元数不是很多的方程求极值，算法的基本思想如下： 
给定n个特征，可以构造一个具有n+1个顶点的单纯形，初始化时需(n+1)*n维矩阵(看成是有n+1个顶点的单纯形) ，矩阵的每一行为n元向量，x0为第一行，xi=x0+r*ei，r为对问题的特征长度大小的估计值，ei为单位向量，x0可初始化为全为1的向量，即认为每个特征权重是相同的，然后选取其余的，在选取过程中，r可以取相同的值也可以取不同的值（r可以看作是对第i个特征权重的调整） 。 
算法运行过程（以机器翻译中的rerank为例）： 
假定BLEU=f(特征的和)，对n+1个顶点（n维向量）分别计算BLEU值（取相反数），然后从中选出BLEU(相反数)最大，次大和最小的三个点，算法每次都是把其中的最大点对应的各权重进行调整，使其变小向最小点靠拢，调整完毕后，计算其对应的BLEU，再从这些BLEU中选出BLEU(相反数)最大，次大和最小的三个点，一直迭代下去，直到最高点到最低点的比率范围合适或达到最大迭代次数为止。 
源码： 
double famoeb(double x[],vector<double> feat) 
{//计算所有特征*权重的和 
double y=0.0; 
for(int i=0;i<FeatNum;i++) 
{ 
y+=x[i+1]*feat[i]; 
} 
return y; 
} 
//单纯形算法 
void amoeba(double p[],double y[],int mp,int np,int ndim,double ftol,int& iter) 
{ 
int i,j,ihi,inhi,mpts,nmax=20; 
double ypr,yprr,rtol,alpha=1.0; 
double beta=0.5; 
double gamma=2.0; 
int itmax=500; 
double pr[21],prr[21],pbar[21]; 
mpts=ndim+1; 
iter=0; 
do 
{ 
int ilo=1; 
if(y[1]>y[2]) 
{ 
ihi=1; 
inhi=2; 
} 
else 
{ 
ihi=2; 
inhi=1; 
} 
for(i=1;i<=mpts;i++) 
{//寻找函数值中的最大，最小和次大值 
if(y[i]<y[ilo]) 
{ 
ilo=i; 
} 
if(y[i]>y[ihi]) 
{ 
inhi=ihi; 
ihi=i; 
} 
else 
{ 
if(y[i]>y[inhi]) 
{ 
if(i!=ihi) 
{ 
inhi=i; 
} 
} 
} 
}//结束寻找各种函数极值 
rtol=2.0*fabs(y[ihi]-y[ilo])/(fabs(y[ihi])+fabs(y[ilo]));//计算从最高点到最低点的比率范围，如合适则返回 
if(rtol<ftol) 
{ 
erase(pbar,prr,pr); 
return; 
} 
if(iter==itmax)//如到了最大的迭代次数，则返回 
{ 
cout<<"amoeba exceeding maximum iterations."<<endl; 
return; 
} 
iter=iter+1;//进行下一次迭代 
for(j=1;j<=ndim;j++) 
{ 
pbar[j]=0.0; 
} 
for(i=1;i<=mpts;i++) 
{ 
if(i!=ihi) 
{ 
for(j=1;j<=ndim;j++) 
{ 
pbar[j]=pbar[j]+p[(i-1)*np+j]; 
} 
} 
} 
for(j=1;j<=ndim;j++) 
{ 
pbar[j]=pbar[j]/ndim; 
pr[j]=(1.0+alpha)*pbar[j]-alpha*p[(ihi-1)*np+j];//求反射点 
} 
vector<int> BestNo; 
ChooseOneBest(pr,numSentences,alldata,StartEndIndices,BestNo); 
//开始计算BLEU值 
vector<pairnum> initialScore(N_gram); 
double referenceLength=0.0;//参考翻译总长度 
for(int k=0;k<numSentences;k++) 
{ 
int sent=BestNo[k];//当前句子的最好候选翻译的序号 
for(int l=0;l<N_gram;l++) 
{ 
initialScore[l].left+=alldata[sent].ngram_data[l].left; 
initialScore[l].right+=alldata[sent].ngram_data[l].right; 
} 
referenceLength+=alldata[sent].closest_length; 
} 
ypr=-BLEU(initialScore,referenceLength);//计算本轮lamda所对应的bleu 
if(ypr<=y[ilo]) 
{//得到一个比最佳点稍好的结果，用gamma做一次外推 
for(j=1;j<=ndim;j++) 
{ 
prr[j]=gamma*pr[j]+(1.0-gamma)*pbar[j]; 
} 
vector<int> BestNo1; 
ChooseOneBest(prr,numSentences,alldata,StartEndIndices,BestNo1); 
//开始计算BLEU值 
vector<pairnum> initialScore1(N_gram); 
double referenceLength1=0.0;//参考翻译总长度 
for(int m=0;m<numSentences;m++) 
{ 
int sent=BestNo1[m];//当前句子的最好候选翻译的序号 
for(int n=0;n<N_gram;n++) 
{ 
initialScore1[n].left+=alldata[sent].ngram_data[n].left; 
initialScore1[n].right+=alldata[sent].ngram_data[n].right; 
} 
referenceLength1+=alldata[sent].closest_length; 
} 
yprr=-BLEU(initialScore1,referenceLength1);//计算本轮lamda所对应的bleu 
if(yprr<y[ilo]) 
{//以扩张点prr作为新的单纯形中的点 
for(j=1;j<=ndim;j++) 
{ 
p[(ihi-1)*np+j]=prr[j]; 
} 
y[ihi]=yprr; 
} 
else 
{//以反射点pr作为单纯形中得点 
for(j=1;j<=ndim;j++) 
{ 
p[(ihi-1)*np+j]=pr[j]; 
} 
y[ihi]=ypr; 
} 
} 
else 
{//反射点不如最佳点，同次高点比较 
if(ypr>=y[inhi]) 
{//反射点不如次高点，取一个中等程度低的点作一次一维收缩 
if(ypr<y[ihi]) 
{ 
for(j=1;j<=ndim;j++) 
{ 
p[(ihi-1)*np+j]=pr[j]; 
} 
} 
y[ihi]=ypr; 
for(j=1;j<=ndim;j++) 
{ 
prr[j]=beta*p[(ihi-1)*np+j]+(1.0-beta)*pbar[j]; 
} 
vector<int> BestNo2; 
ChooseOneBest(prr,numSentences,alldata,StartEndIndices,BestNo2); 
//开始计算BLEU值 
vector<pairnum> initialScore2(N_gram); 
double referenceLength2=0.0;//参考翻译总长度 
for(int s=0;s<numSentences;s++) 
{ 
int sent=BestNo2[s];//当前句子的最好候选翻译的序号 
for(int t=0;t<N_gram;t++) 
{ 
initialScore2[t].left+=alldata[sent].ngram_data[t].left; 
initialScore2[t].right+=alldata[sent].ngram_data[t].right; 
} 
referenceLength2+=alldata[sent].closest_length; 
} 
yprr=-BLEU(initialScore2,referenceLength2);//计算本轮lamda所对应的bleu 
if(yprr<y[ihi]) 
{//以prr作为新单纯形中的点 
for(j=1;j<=ndim;j++) 
{ 
p[(ihi-1)*np+j]=prr[j]; 
} 
y[ihi]=yprr;//更新当前最高点出的函数值 
} 
else 
{//单纯性太大，缩小原来的单纯形 
for(i=1;i<=mpts;i++) 
{ 
if(i!=ilo) 
{ 
for(j=1;j<=ndim;j++) 
{ 
pr[j]=0.5*(p[(i-1)*np+j]+p[(ilo-1)*np+j]); 
p[(i-1)*np+j]=pr[j]; 
} 
vector<int> BestNo3; 
ChooseOneBest(pr,numSentences,alldata,StartEndIndices,BestNo3); 
//开始计算BLEU值 
vector<pairnum> initialScore3(N_gram); 
double referenceLength3=0.0;//参考翻译总长度 
for(int u=0;u<numSentences;u++) 
{ 
int sent=BestNo3[u];//当前句子的最好候选翻译的序号 
for(int v=0;v<N_gram;v++) 
{ 
initialScore3[v].left+=alldata[sent].ngram_data[v].left; 
initialScore3[v].right+=alldata[sent].ngram_data[v].right; 
} 
referenceLength3+=alldata[sent].closest_length; 
} 
y[i]=-BLEU(initialScore3,referenceLength3);//计算本轮lamda所对应的bleu 
} 
} 
} 
} 
else 
{//反射点好于次高点，以反射点pr作为单纯形中得点 
for(j=1;j<=ndim;j++) 
{ 
p[(ihi-1)*np+j]=pr[j]; 
} 
y[ihi]=ypr; 
} 
} 
}while(1); 
}