topic Vectorizing in SAS/IML Software and Matrix Computations

Vectorizing

CerditoSalvaje — Sat, 29 Jan 2022 18:17:30 GMT

Hello to everyone,

I have asked so many questions about this subject, I am trying to code an algorithm for the Particle Swarm Optimization. I have asked in this forum some days ago and Rick_SAS told me that it would be interesting to vectorize my code. So I have read some articles about the subject and I have taken a sheet of paper and a pen, for trying to achieve the result. I think that I don't have the logic behind and it's a little bit frustrating. So I don't know if someone can show me a little example with this code. In order to try to see the logic with my own code. Here you have a little fragment of my code, I have tried so many things and I failed. Do not hesitate to ask me if you need some additional explanation. Thank you all.

xmin=-100;
xmax=100;
/***parameters and variables***/
n=20;						/*number of particules*/
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);

do run=1 to maxrun; /*this line is not important in this example*/
	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;
end;

Re: Vectorizing

sbxkoenk — Sat, 29 Jan 2022 19:35:38 GMT

Hello,

As this is SAS/IML code, I have moved your post to the

SAS/IML Software and Matrix Computations

board.

No absolute need to do it in SAS/IML by the way.
See here :
Particle Swarm Optimization in SAS
Gregory Watson, UCLA Department of Biostatistics, Los Angeles, CA
https://www.lexjansen.com/wuss/2014/63_Final_Paper_PDF.pdf

Cheers,

Koen

Re: Vectorizing

Rick_SAS — Sat, 29 Jan 2022 21:12:57 GMT

I applaud you for taking on this challenge. It might be frustrating in the short term, but it will pay off in the long term. See these tips: https://blogs.sas.com/content/iml/2013/11/18/vectorizing-the-construction-of-a-structured-matrix.html

1. Assuming that the objective function involves fractional powers (as in your last post), I have adjusted xmin=0.

2. I have also vectorized the objective function.

3. If you are not familiar with the '#' multiplication operator, the operation v#Y is shorthand notation for multiplying the i_th column of the matrix Y by i_th element of the row vector, v. See "Shorthand notation for row and column operations."

proc iml;
xmin=0;     /* min has to be >= 0 or else x^(1/4) and x^(3/4) are undefined */
xmax=100;
/***parameters and variables***/
n=20;			   			/*number of particles*/
p=2;
maxt=1000;					/*maximum iteration*/
maxrun=1;					/*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);

/* vectorization of f: R^2 --> R */
start f(x);
   con=100*x[,1]+200*x[,2]-1500;
   fc = 3*(x[,1]##0.25) # (x[,2]##0.75);
   pen=10**9;
   idx = loc(con > 0);
   if ncol(type)>0 then 
      fc[idx] = pen*con[idx];
   return(fc);
finish f;

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

print X[c={'x1' 'x2' 'F(x1,x2)'}];

I understand that vectorized operations can be hard to learn, so write back if you need anything clarified. When you vectorize, the resulting code is very compact. It is like reading poetry: there is a lot of information in each line.

Re: Vectorizing

CerditoSalvaje — Sat, 29 Jan 2022 22:55:48 GMT

Thank you so much. I will start right now. If I get it for the whole code I'll tell you. Thanks.