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

$\mathrm{BFGS:\ Initialize\ population\ } X$

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

$\mathrm{BFGS:\ for\ } i=1 \mathrm{\ to\ } |X|\ H_i=I$

H=zeros(2*n,dim,dim);for i=1:n,H(i,:,:)=eye(dim);end

$\mathrm{BFGS:\ }k=0$

k=0;

$\mathrm{BFGS:\ while\ } f_{best}\ne f_{opt}$

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;

$\mathrm{BFGS:\ for\ } i=1 \mathrm{\ to\ } |X|$

  for i=1:n

$\mathrm{BFGS:\ }d=-H_i\nabla f_i$

    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

$\mathrm{BFGS:\ }\alpha \mathrm{\ is\ computed\ from\ a\ line\ search\ procedure\ to\ satisfy\ the\ Wolfe\ conditions}$

      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

$\mathrm{BFGS:\ else,\ }X_i^+=X_i+\alpha d$

        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

$\mathrm{DE:\ for\ } i=1 \mathrm{\ to\ } |X|$

  for i=n+1:2*n

$\mathrm{DE:\ } O_i^+=\mathrm{DE\_candidate}(X,i)$

    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);

$\mathrm{DE:\ }X^+=\mathrm{best\ }|X|\mathrm{\ individuals\ from\ }\{X^+,O^+\}$

  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;

$\mathrm{BFGS:\ for\ } i=1 \mathrm{\ to\ } |X|$

  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

initialize a list of secant equations E to the previous parent found

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

$\mathrm{if\ }K(S_{1,2,\dots,i})>K_{max} | K(S_{1,i})>K_{max} | K(S_{2,i})>K_{max},\dots,| K(S_{i-1,i})>K_{max}$

$\mathrm{\ eliminate\ }S_i\mathrm{\ and\ }Y_i$

    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

$\mathrm{BFGS:\ }S_j=X_i^+-X_j\mathrm{\ and\ }Y_j=\nabla f_i^+-\nabla f_j$

      if mxcond<threshold,S=[S P(i,:)'-P0(j,:)'];Y=[Y G(i,:)'-G0(j,:)'];end
    end
    I=permute(H(i,:,:),[2 3 1]);

$\mathrm{Schnabel:\ }(S,Y) = \mathrm{Schnabel\_Algorithm\_3.31}(X,X_i^+)$

    [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

$\mathrm{BFGS:\ }B_i^+=B_i-B_i S(S^T B_i S)^{-1}S^T B_i+Y(Y^T S)^{-1}Y^T$

    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

$X_n^+=\mathrm{\ rotation\ in\ span}(X^+)^\perp\mathrm{\ of\ }X_1^+\mathrm{\ around\ }\Sigma_{i=1}^{|X|/2} w_i X_i^+$

  [~,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';

$\mathrm{BFGS:\ }k=k+1$

  k=k+1;
end
x=fg('xbest');