Obsidian | Level 7

Vectorizing the Particle Swarm Optimization

Hi all. As you can see in my profile, I am writing code about Particle Swarm Optimization in SAS. A few days ago, Rick recommended that I vectorize the algorithm to make it more optimal and perform better, as well as being much more aesthetically pleasing and easier to understand.

Here is the original code 100% functional and without vectorizing:

``````/*creating a template so as to display our results*/
proc template;
/*the main structure of this template comes from the SAS documentation (https://documentation.sas.com/doc/en/pgmsascdc/v_021/imlug/imlug_tables_sect015.htm)*/
define table PSO_Results;
column Bestfun bx1 bx2 average gradx1 gradx2;     /* names and order of columns */
text "Best results of the simulation";
end;
text "Best variables";     /* title of spanning header */
start = x1;                  /* span starts at second column */
end   = x2;                  /* span ends at third column */
end;
end;
run;

/****************************************BASIC PSO****************************************/
proc iml;
/*****initialization*****/
/*function we want to minimize*/
start f(x);
y= x[1]**2 + x[2]**2 -100 ;
return(y);
finish f;

xmin=-100;
xmax=100;

/***parameters and variables***/
n=20;						/*number of particules*/
p=2;						/*number of dimensions*/
c1=2.05;						/*acceleration factor*/
c2=2.05;						/*acceleration factor*/
maxt=1000;					/*maximum iteration*/
maxrun=100;					/*trials of the algorithm*/
X=j(n,p+1,0);				/*initialization of the particle's matrix*/
lb=j(1,p,xmin);
ub=j(1,p,xmax);

/*other useful matrixes*/
ffmin=j(maxt,maxrun,0);
fft=j(1,maxrun,0);
fff=j(1, maxrun,0);
rgbest=j(maxrun,p,0);
nb_t=j(2, maxrun,0);

do run=1 to maxrun;
nb_t[1,run]=run;
/***initialization of the pso***/
do i=1 to n;
do j=1 to p;
/*sample in a uniform distribution to initialize the particles position*/
call randgen(u,'uniform');
X[i,j]=lb[j]+u*(ub[j]-lb[j]);
end;
/*to fill the third column of X (which corresponds to the output of the function), we evaluate the fitting for each particule*/
X[i,p+1]=f(X[i,]);
end;

pbest=X;			/*initialization of pbest so we can iterate then*/
V=0.1*X;			/*initialization of the velocity (10%) */

/*we find the current minimum of f, given the current values of X (if we want to minimize the function)*/
titi=min(pbest[,p+1]);

/*then we iterate to assign the coordinates of this minimum to gbest*/
do i=1 to n;
if pbest[i,p+1]=titi then gbest=pbest[i,];
else e=0;
end;

*	print X, pbest, gbest;
/**********************************/
t=1;							/*resetting the iteration counter*/
tolerance=1;					/*resetting the tolerance(stopping criteria)*/
do until(tolerance < (10**-(12))|t=maxt);

/*for all the particles in all dimensions : */
do i=1 to n;
do j=1 to p;

call randgen(u1,'uniform');
call randgen(u2,'uniform');
V[i,j]=V[i,j]+c1*u1*(pbest[i,j]-X[i,j])+c2*u2*(gbest[1,j]-X[i,j]);

X[i,j]= X[i,j]+V[i,j];

/*handling boundary violations*/
if X[i,j]<lb[j] then X[i,j]=lb[j];
else e=0;
if X[i,j]>ub[j] then X[i,j]=ub[j];
else e=0;

end;
end;

/*computing fitness for all the particules + updating pbest*/
do i=1 to n;

X[i,p+1]=f(X[i,]);

if X[i,p+1]<pbest[i,p+1] then pbest[i,]=X[i,];
else e=0;

end;

ffmin[t,run]=gbest[,p+1];
fft[run]=t;

/* updating gbest*/
titi=min(pbest[,p+1]);
do i=1 to n;
if pbest[i,p+1]=titi then gbest=pbest[i,];
else e=0;
end;

/*calculating tolerance*/
if t>100 then tolerance=abs(ffmin[(t-100),run]-gbest[,p+1]);
else e=0;

t=t+1;

end;
*print gbest;
nb_t[2,run]=t;					/*storing the number of iterations it took the algorithm to converge*/
fff[run]=gbest[,p+1];
rgbest[run,]=gbest[,1:p];
end;

/*displaying our best results by putting them in a template*/
fff2=fff`;
bestfun=mean(fff2);			/*taking the best run among all of the runs of the simulation*/
bx1=mean(rgbest[1]);
bx2=mean(rgbest[2]);

*print bestrun, bestfun, best_variables;

titi=nb_t`;
average=mean(titi[,2]);

best_variables= bx1||bx2;
opt={1 1};

toto=bestfun || bx1 || bx2 || average || gradx;			/*we concatenate the main results*/

call TablePrint(best_results) template="PSO_Results";

quit;``````

Here all my progress made but I have several problems that I mention in the code and that I would like to tell you about here :

``````/*creating a template so as to display our results*/
proc template;
/*the main structure of this template comes from the SAS documentation (https://documentation.sas.com/doc/en/pgmsascdc/v_021/imlug/imlug_tables_sect015.htm)*/
define table PSO_Results;
column Bestfun bx1 bx2 average gradx1 gradx2;     /* names and order of columns */
text "Best results of the simulation";
end;
text "Best variables";     /* title of spanning header */
start = x1;                  /* span starts at second column */
end   = x2;                  /* span ends at third column */
end;
end;
run;
proc iml;

/* vectorization of f: R^2 --> R */
start f(x);
ros= 100*(x[,2]-x[,1]##2)##2+(x[,1]-1)##2;
return(ros);
finish f;

xmin=-5;     /* min has to be >= 0 or else x^(1/4) and x^(3/4) are undefined */
xmax=10;
/***parameters and variables***/
n=20;			   			/*number of particles*/
p=2;
maxt=1000;					/*maximum iteration*/
maxrun=100;					/*trials of the algorithm*/
X=j(n,p+1,0);				/*initialization of the particle's matrix*/
lb=j(1,p,xmin);
ub=j(1,p,xmax);
c1=2.05;						/*acceleration factor*/
c2=2.05;						/*acceleration factor*/

phi=c1+c2;
chi= 2/abs(2-phi-sqrt(phi**2 -4*phi));

/*other useful matrixes*/
ffmin=j(maxt,maxrun,0);
fft=j(1,maxrun,0);
fff=j(1, maxrun,0);
rgbest=j(maxrun,p,0);
nb_t=j(2, maxrun,0);

call randseed(1234);
Y = j(n,P);/* allocate for random values */
do run=1 to maxrun; /*this line is not important in this example*/

/***initialization of the pso***/
call randgen(Y, 'uniform');

/*sample in a uniform distribution to initialize the particles position*/
X[,1:p] = lb + (ub-lb)#Y;
/*to fill the third column of X (which corresponds to the output of the function), we evaluate the fitting for each particule*/
X[,p+1]=f(X);

pbest=X;			/*initialization of pbest so we can iterate then*/
V=0.1*X[,1:p];

/*then we iterate to assign the coordinates of this minimum to gbest*/
titi=min(pbest[,p+1]);
minIdx = loc(pbest[,p+1]=titi);
if ncol(minIdx)>0 then gbest = pbest[minIdx, ];
else e=0;
t=1;							/*resetting the iteration counter*/
tolerance=1;					/*resetting the tolerance(stopping criteria)*/

do until(tolerance < (10**-(12))|t=maxt);
/*for all the particles in all dimensions : */
call randseed(1234);
u1 = j(n,p);
u2 = j(n,p);
call randgen(u1,'uniform');
call randgen(u2,'uniform');

V=chi*(V+c1*u1#(pbest[,1:p]-X[,1:p]) + c2*u2#(gbest[,1:p]-X[,1:p]));

X[,1:p] = X[,1:p] + V;

/*handling boundary violations*/
X[,1:p]=X[,1:p]<>lb;
X[,1:p]=X[,1:p]><ub;

/*computing fitness for all the particules + updating pbest*/

X[,p+1]=f(X[,1:p]);

inferior_Idx=loc(X[,p+1]<pbest[,p+1]);
if ncol(inferior_Idx)>0 then pbest[inferior_Idx,]=X[inferior_Idx,];
else e=0;

ffmin[t,run]=gbest[,p+1];
fft[run]=t;

/* updating gbest*/
titi=min(pbest[,p+1]);

minIdx_2 = loc(pbest[,p+1]=titi);
if ncol(minIdx_2)>0 then
gbest = pbest[minIdx, ];

/*calculating tolerance*/
if t>100 then tolerance=abs(ffmin[(t-100),run]-gbest[,p+1]);
else e=0;

t=t+1;

end;
*print gbest;
nb_t[2,run]=t;					/*storing the number of iterations it took the algorithm to converge*/
fff[run]=gbest[,p+1];
rgbest[run,]=gbest[,1:p];
end;
/*displaying our best results by putting them in a template*/
fff2=fff`;
bestfun=mean(fff2);			/*taking the best run among all of the runs of the simulation*/
bx1=mean(rgbest[,1]);
bx2=mean(rgbest[,2]);

titi=nb_t`;
print titi;
average=mean(titi[,2]);

best_variables= bx1||bx2;
opt={1 1};

toto=bestfun || bx1 || bx2 || average || gradx;			/*we concatenate the main results*/

call TablePrint(best_results) template="PSO_Results";

quit;``````

To begin with, the definition of bounds is complicated. In this case they work but only for 2 columns or p=2. I would like to generalize it but I don't know how.

Another problem lies in the fact of speed. Gbest is a 1x3 format vector. while X in this case is nx2. subtracting from each other is tricky. Actually what I want is for it to be subtracted line by line. I've thought about concatenating gbest n times so I can subtract it. But it seems to me that there has to be another way to do it.

At the moment everything is ready, although I still have many pending things. Like for example, how do you remove the do until loop in vector format?

Thank you so much everyone. I am learning more than ever with all your advice. And being a bit heavy.

All the best.

1 ACCEPTED SOLUTION

Accepted Solutions
SAS Super FREQ

Re: Vectorizing the Particle Swarm Optimization

Here's another hint. There is n error on the line

inferior_Idx=loc(X[,p+1]<pbest[,p+1]);
pbest=X[inferior_Idx,];

ERROR: (execution) Matrix has not been set to a value.

operation : [ at line 2026 column 14
operands : X, inferior_Idx,
X 20 rows 3 cols (numeric)

inferior_Idx 0 row 0 col (type ?, size 0)

The error is because no element satisfies the condition you are searching for. Anytime you use LOC to search for elements that satisfy a criterion, you need to ask yourself whether it is possible that no elements exist.  If so, the LOC function will return the empty matrix.  Thus, you shouldn't use the results form LOC without checking that some result was found. See the article, "Beware the naked LOC."

7 REPLIES 7
SAS Super FREQ

Re: Vectorizing the Particle Swarm Optimization

I don't have time to fix all your problems, but

1. For the lower bounds, study this example

``````
proc iml;
p = 2;
X = {1 2  3,
-1 3  5,
3 0  7,
2 4  9};
lb = {1 2};

T = (X[ ,1:p] <> lb);
print T;

X[ ,1:p] = (X[ ,1:p] <> lb);
print X;``````

2. You should not try to remove loops that are controlling the convergence of an iterative algorithm. Is that what your loop is doing?

Obsidian | Level 7

Re: Vectorizing the Particle Swarm Optimization

Yes, it is. Don't worry Rick. You have already helped me a lot. Thanks a lot. It's okay like that.
SAS Super FREQ

Re: Vectorizing the Particle Swarm Optimization

Here's another hint. There is n error on the line

inferior_Idx=loc(X[,p+1]<pbest[,p+1]);
pbest=X[inferior_Idx,];

ERROR: (execution) Matrix has not been set to a value.

operation : [ at line 2026 column 14
operands : X, inferior_Idx,
X 20 rows 3 cols (numeric)

inferior_Idx 0 row 0 col (type ?, size 0)

The error is because no element satisfies the condition you are searching for. Anytime you use LOC to search for elements that satisfy a criterion, you need to ask yourself whether it is possible that no elements exist.  If so, the LOC function will return the empty matrix.  Thus, you shouldn't use the results form LOC without checking that some result was found. See the article, "Beware the naked LOC."

Obsidian | Level 7

Re: Vectorizing the Particle Swarm Optimization

Yes, Rick. Thank you so much. Problem solved.
Obsidian | Level 7

Re: Vectorizing the Particle Swarm Optimization

I have put some if statement for solving the naked loc.
SAS Super FREQ

Re: Vectorizing the Particle Swarm Optimization

A general technique for debugging a program is to run the program, then look for the FIRST error in the log. Figure out how to fix that error. Rerun. Fix the new first error. Repeat until the program runs all the way through.

That will take care of the syntax and run time errors. It might still contain logical errors. You can read about the different kinds of programming errors.

Calcite | Level 5

Re: Vectorizing the Particle Swarm Optimization

it worked for me, thanks for the answers.   roblox tutorials mobdro

From The DO Loop