MATLAB Examples

## Contents

function x=BFGSE(FUN,dim,ftarget,maxfunevals,range) 
n=dim;%population size m=n-1;%extra secant equations Cr=0.1;%crossover probability rangemn=range(1)+fg('os');rangemx=range(2)+fg('os');%domain 
P=zeros(2*n,dim);P(1:n,:)=(rangemx-rangemn)*rand(n,dim)+rangemn; f=zeros(2*n,1);G=zeros(n,dim);for i=1:n,[g,f(i)]=fg_AD_BBOB(FUN,P(i,:)',[],P,f);G(i,:)=g';end%evaluation 
H=zeros(2*n,dim,dim);for i=1:n,H(i,:,:)=eye(dim);end 
k=0; 
while fg('fbest')>=ftarget & fg('evaluations')<maxfunevals 
 Gradbeta=zeros(2*n,dim);Q=[1:n 1:n];%line search return values P0=P; G0=G;%old gradient P00=P;G00=G;H00=H;%all previous generation information lsgrad=zeros(2*n,1); P(n+1:2*n,:)=P(1:n,:);G0(n+1:2*n,:)=G0(1:n,:);H(n+1:2*n,:,:)=H(1:n,:,:);%inherited information P0(n+1:2*n,:)=P0(1:n,:);%old X f0=f; 
 for i=1:n 
 p=-permute(H(i,:,:),[2 3 1])*G(i,:)'; if norm(p)<1e-7 P(i,:)=(rangemx-rangemn)*rand(1,dim)+rangemn; lsgrad(i)=1; else 
 fg('toogle_queue') [alpha,~,Falpha,gradbeta,Q(i),fail]=linesch_armijo_find_parent(P(i,:)',... f0(i),G(i,:)*p,p,[],0,1,0,P0(1:n,:),G0(1:n,:),f0,i,permute(H(i,:,:),[2 3 1])); fg('toogle_queue') if fail==0 & Q(i)==-1,Q(i)=i;end%if no parent is adecuate, set the i-th individual as parent Gradbeta(i,:)=gradbeta; 

## if line search fails or d is small, resample

 if fail | alpha==-Inf | (norm(p)<1e-7 & alpha<1 & alpha>0)% | norm(p*alpha)<1e-5%| norm(p)<1e-5% & flag>10^5) P(i,:)=(rangemx-rangemn)*rand(1,dim)+rangemn; lsgrad(i)=1; else 
 P(i,:)=P(i,:)+alpha*p'; 
 end 
 end if lsgrad(i),f(i)=feval(FUN,P(i,:)');else f(i)=Falpha;end 

## if f overflows, resample

 if isnaninf(f(i)) P(i,:)=(rangemx-rangemn)*rand(1,dim)+rangemn; H(i,:,:)=eye(dim);f(i)=feval(FUN,P(i,:)'); lsgrad(i)=1;continue end 
 end 
 for i=n+1:2*n 
 p=DE_candidate_component(P0(1:n,:),i-n); jrand=randi(dim); for q=1:dim,if rand>=Cr & q~=jrand,p(q)=P(i,q);end,end p=p-P(i,:)'; 

## if mutant vector is small or an ascent direction, discard DE individual

 if norm(p)<1e-5 | G(i-n,:)*p>=0,f(i)=Inf;continue,end P(i,:)=P(i,:)+p';[f(i),P(i,:)]=feval(FUN,P(i,:)'); 

## if f overflows, resample

 if isnaninf(f(i)) P(i,:)=(rangemx-rangemn)*rand(1,dim)+rangemn; H(i,:,:)=eye(dim);f(i)=feval(FUN,P(i,:)');continue end 

## if the sufficient descent condition from the target vector is not met, discard DE individual

 if f(i)>f0(i-n)+.0001*G(i-n,:)*p,f(i)=Inf;continue,end 
 end queue=fg('queue'); if ~isempty(queue) [~,idx]=max(f(n+1:2*n)); if queue.f<f(n+idx) f(n+idx)=queue.f; P(n+idx,:)=queue.x; end end next=n+1;%next points to an individual to replace inadequate individuals [~,arindex]=sort(f); 
 P=P(arindex,:);P0=P0(arindex,:);G0=G0(arindex,:);f=f(arindex);H=H(arindex,:,:); lsgrad=lsgrad(arindex);Gradbeta=Gradbeta(arindex,:);Q=Q(arindex); rep=zeros(n,1);%marks DE individuals whose target vector is in the population for i=1:n,if arindex(i)>n && any(arindex(1:n)==arindex(i)-n),rep(i)=1;end,end i=1; 
 while i<=n 
 if arindex(i)<=n & ~lsgrad(i),G(i,:)=Gradbeta(i,:);else,G(i,:)=fg_AD_BBOB(FUN,P(i,:)',f(i),P,f);end 

## if f or the gradient overflows, resample and skip the hessian update

 if isnaninf(f(i)) | isnaninf(G(i,:)) P(i,:)=(rangemx-rangemn)*rand(1,dim)+rangemn;H(i,:,:)=eye(dim); [g,f(i)]=fg_AD_BBOB(FUN,P(i,:)',[],P,f);G(i,:)=g'; i=i+1;continue end if k>0 & lsgrad(i),H(i,:,:)=eye(dim);i=i+1;continue,end 

## if a BFGS individual has no line search parent, skip the hessian update

 if Q(i)==-1,i=i+1;continue,end 

## choose previouos parent found if possible or the nearest new parent found

 q=find_nearest_parents(P00(1:n,:),G00(1:n,:),G(i,:)',P(i,:)',i); if any(q==Q(i)),q=Q(i);else q=q(1);end if q==-1 

## if no parent is adequate for a DE individual, resample and skip the hessian update

 if arindex(i)>n if next<=2*n && ~isnaninf(f(arindex(next))) P(i,:)=P(arindex(next),:);f(i)=f(arindex(next));next=next+1;continue else P(i,:)=(rangemx-rangemn)*rand(1,dim)+rangemn;H(i,:,:)=eye(dim); [f(i),P(i,:)]=feval(FUN,P(i,:)');i=i+1;continue end else 

## if no parent is adequate for a BFGS individual, set parent as previous parent found

 q=Q(i); 
 end 
 end P0(i,:)=P00(q,:);G0(i,:)=G00(q,:);H(i,:,:)=H00(q,:,:);%information inherited by DE offspring 

## make a list L of individuals sorted by distance to the i-th individual

 dist=zeros(n,1);for q=1:n,dist(q)=norm(P(i,:)-P0(q,:));end,[~,idx]=sort(dist); 

## make sqrt(n) clusters

 warning('off','stats:kmeans:EmptyCluster'),warning('off','stats:kmeans:FailedToConverge') kmeansk=round(sqrt(n/2));IDX=kmeans(P(1:n,:)/max(max(abs(P(1:n,:)))),kmeansk,'emptyaction','drop'); 

## send individuals in the cluster of the i-th individual to the end of L

 for q=2:n qq=0; while IDX(idx(q))==IDX(idx(1)) qq=qq+1;if qq==n-q+2,break,end,idx(q:n)=circshift(idx(q:n),-1); end end 

## for DE individuals E is the target vector equation

 if arindex(i)<=n,ind=Q(i);else ind=i;end 

## append to E the secant equations in the order of L which satisfy the curvature condition

 if arindex(i)<=n for q=1:n if any(ind==idx(q))==0 & (G(i,:)-G0(idx(q),:))*(P(i,:)'-P0(idx(q),:)')>0%2nd Wolfe if ~rep(idx(q)),ind=[ind idx(q)];end%avoid YT nan because of repeated P0 end end end 

## insert the elite secant equation in E after the original secant equation

 if arindex(i)<=n [~,idx]=min(f(1:n));if ~any(ind==idx),ind=[ind idx];end for q=length(ind)-1:-1:2,if ind(q+1)==idx,ind(q+1)=ind(q);ind(q)=idx;end,end end 
 S=zeros(dim,0);Y=zeros(dim,0); threshold=10^4; for q=1:length(ind) 
 j=ind(q); if isnaninf(P(i,:)'-P0(j,:)'),continue,end mxcond=cond([S P(i,:)'-P0(j,:)']); for qq=1:size(S,2) if mxcond>threshold,break,end mxcond=max(mxcond,cond([S(:,qq) P(i,:)'-P0(j,:)'])); end 
 if mxcond<threshold,S=[S P(i,:)'-P0(j,:)'];Y=[Y G(i,:)'-G0(j,:)'];end 
 end I=permute(H(i,:,:),[2 3 1]); 
 [YT,ST]=perturbation(Y,S,m,'3.31',I); 

## if Y overflows or is empty, set (S,Y) to original secant pair

 if isnaninf(YT) | isempty(YT) S=P(i,:)'-P0(i,:)';YT=G(i,:)'-G0(i,:)'; else 
 S=ST; 

## if S'Y is nonpositive, resample and skip the hessian update

 if S'*YT<=0 P(i,:)=(rangemx-rangemn)*rand(1,dim)+rangemn;H(i,:,:)=eye(dim); [g,f(i)]=fg_AD_BBOB(FUN,P(i,:)',[],P,f);G(i,:)=g';i=i+1;continue end 
 end 
 warning('off','MATLAB:singularMatrix') H(i,:,:)=I+(S-I*YT)*inv(S'*YT)*S'+S*inv(S'*YT)*(S-I*YT)'-... S*inv(S'*YT)*(S-I*YT)'*YT*inv(S'*YT)*S'; 

## if hessian overflows, resample

 if isnaninf(H(i,:,:)) P(i,:)=(rangemx-rangemn)*rand(1,dim)+rangemn;H(i,:,:)=eye(dim); [g,f(i)]=fg_AD_BBOB(FUN,P(i,:)',[],P,f);G(i,:)=g';i=i+1;continue end i=i+1; 
 end 
 [~,arindex]=sort(f(1:n));mu=floor(n/2);weights=log(mu+1/2)-log(1:mu)';xmean=P(arindex(1:mu),:)'*weights; weights=[weights;zeros(n-mu,1)];Z=Null(P(1:n,:));z=size(Z,2); if z==0 sigma=sqrt(var(P(arindex,:),weights));V=normrnd(0,sigma); else inverse=P(arindex(1),:)';inverse=fxISO(xmean,P(arindex(1),:)'); V=Z*(2*rand(z,1)-1);V=V'/norm(V)*norm(inverse-xmean); end if randi(2)==1 V=zeros(1,dim); end P(arindex(n),:)=P(arindex(1),:)+V;H(arindex(n),:,:)=eye(dim); [g,f(arindex(n))]=fg_AD_BBOB(FUN,P(arindex(n),:)',[],P,f);G(arindex(n),:)=g'; 
 k=k+1; 
end x=fg('xbest');