topic Re: Getting product of matrix in SAS Programming

Getting product of matrix

rahul88888 — Wed, 06 Apr 2016 06:35:06 GMT

I am getting a wrong result for a bigger matrix however for a small matrix like this i am getting results correctly. Can anybody help me what i may be doing wrong

data test1;

input pol $ prod1-prod9;

datalines;

1abc 1 1 0 1 0 0 1 1 0

2abc 0 0 1 1 0 0 1 1 1

3abc 1 1 0 1 0 0 0 0 0

4abc 1 0 1 1 0 0 1 1 0

5abc 0 1 0 1 1 0 0 1 0

6abc 1 1 0 1 0 0 0 0 0

7abc 1 0 1 1 0 0 1 1 0

8abc 0 1 0 1 1 0 0 1 0

9abc 1 1 0 1 0 0 0 0 0

10abc 1 0 1 1 0 0 1 1 0

11abc 0 1 0 1 1 0 0 1 0

12abc 1 1 0 1 0 0 0 0 0

13abc 1 0 1 1 0 0 1 1 0

14abc 0 1 0 1 1 0 0 1 0

15abc 1 1 0 1 0 0 0 0 0

16abc 1 0 1 1 0 0 1 1 0

17abc 0 1 0 1 1 0 0 1 0

;

run;

options mprint mlogic merror symbolgen;

%macro gen_prod_matrix(ilib,ids,olib,oPds) ;

proc sql;

select num_numeric into :prod_size from dictionary.tables where upcase(memname)=upcase(&ids) and upcase(libname)=upcase(&ilib);

SELECT distinct name into :varb_list separated by ' ' FROM dictionary.columns where upcase(memname)=upcase(&ids) AND upcase(libname)=upcase(&ilib) and TYPE='num';

SELECT name into :char_var separated by ' ' FROM dictionary.columns where upcase(memname)=upcase(&ids) AND upcase(libname)=upcase(&ilib) and TYPE='char';

SELECT memname into :ds FROM dictionary.columns where upcase(memname)=upcase(&ids) AND upcase(libname)=upcase(&ilib) and TYPE='num';

SELECT libname into :lib FROM dictionary.columns where upcase(memname)=upcase(&ids) AND upcase(libname)=upcase(&ilib) and TYPE='num';

quit;

%let prod_squ=%eval(&prod_size*&prod_size);

%let diag_siz=%eval(&prod_size+1);

/************************************************/

data step1;

set &lib..&ds;

array prod [&prod_size] &varb_list;

array nmres[&prod_size,&prod_size];

do i=1 to &prod_size;

do j=1 to &prod_size;

val=prod[i]*prod[j];

nmres[i,j]+val;

end;

run;

/*************************************************/

data res1(keep=nmres: key);

set step1 end=last;

key=1;

if last;

run;

/**************************************************/

%macro diag;

%do i=1 %to &prod_squ %by &diag_siz;

nmres&i

%end;

%mend;

%let diag_var=%diag ;

/****************************************************/

data diagonal;

set res1;

keep &diag_var key;

run;

/************************************************/

data step2(keep=dmres: key);

set diagonal;

array diag[&prod_size] nmres:;

array dmres[&prod_size,&prod_size];

do i=1 to &prod_size;

do j=1 to &prod_size;

dmres[i,j]=sqrt(diag[i])*sqrt(diag[j]);

end;

key=1;

run;

/**********************************************************************/

proc sql;

create table dvd_ds as

select * from res1 t1 inner join step2 t2 on t1.key=t2.key;

quit;

/**********************************************************************/

data res3(keep=res:);

set dvd_ds;

array dm[&prod_squ] dmres:;

array nm[&prod_squ] nmres:;

array res[&prod_squ];

do i=1 to &prod_squ;

res[i]=nm[i]/dm[i];

end;

run;

/*******************************************************/

%macro prod;

%do i=1 %to &prod_size;

prod&i

%end;

%mend;

%let prod_var=%prod ;

/****************************************************/

data matrix_op(keep=prod:);

set res3;

array results[&prod_squ] res:;

array products[&prod_size] &prod_var;

arrange=%eval(&prod_size-1);

do i=1 to &prod_squ by &prod_size;

do j=i to i+arrange;

x+1;

products[x]=results[j];

end;

output;

x=0;

end;

run;

/*************************************************************/

data &OLIB..&OPDS;

set matrix_op;

run;

%mend;

/*******************************************************************/

%macro colab_filter(in_A =,in_B=,ou_AB=);

/*determine number of rows and columns in the 2nd matrix*/

%let B_id=%sysfunc(open(&in_B));

%let B_rows=%sysfunc(attrn(&B_id,nobs));

%let B_cols=%sysfunc(attrn(&B_id,nvars));

%let rc=%sysfunc(close(&B_id));

/*transpose the 2nd matrix*/

proc transpose data=&in_B out=t&in_B(drop=_:);run;

/*making Cartesian product of the 1st and transposed 2nd matrices*/

data &ou_AB;

do until(eofA);

set &in_A end=eofA;

do i=1 to n;

set t&in_B nobs=n point=i;

output;

end;

run;

/*multiplication*/

data &ou_AB;

/*new columns for products, equal to number of columns in the 2nd matrix*/

array v[&B_cols];

do j=1 to &B_cols;

v[j]=0;

set &ou_AB;

array col _ALL_;

/*multiply corresponding pairs of columns*/

do i=&B_cols+2 to &B_cols+1+&B_rows;

v[j]+col[i]*col[i+&B_rows];

end;

output;

keep v:;

run;

%mend;

/**********************************/

%gen_prod_matrix('work','test1',work,PROD_MATRIX_test1);

%colab_filter(in_a=test1(DROP=f1:),in_b=PROD_MATRIX_test1,ou_AB=mat_prod_test1);

/**********************************/

Re: Getting product of matrix

Ksharp — Wed, 06 Apr 2016 09:40:21 GMT

If you are talking about Matrix Operator, That is a better question to IML. 
Post it at IML Forum, and Post some data and the output you want, and the logic behind that operator.

Re: Getting product of matrix

rahul88888 — Wed, 06 Apr 2016 10:24:36 GMT

Sorry I dont have proc IML i am using SAS 9.2 and it says

ERROR: Bad product ID for procedure IML.

Thats why rather than posting under IML i posted under base SAS and i wanted it to get done without IML. I am trying proc FCMP to create an user difined function hope this works out

Re: Getting product of matrix

FreelanceReinh — Wed, 06 Apr 2016 10:44:19 GMT

Have you tried the matrix CALL routines available through PROC FCMP?

Otherwise, I think you should clarify what you are trying to achieve and what exactly goes wrong.

You are talking about "product of matrix." So, you want to implement matrix multiplication in Base SAS?
If this is the case, how does input dataset TEST1 relate to this task? This dataset could be regarded as a 17x9 matrix with an additional character column (POL) of little use, but normally a product involves at least two factors.
You state that "for a small matrix like this i am getting results correctly." However, when I run your code, I get notes "Division by zero detected" and "Mathematical operations could not be performed," which are frowned upon, and a strange matrix with missing values in one row and one column (from the first macro call, %gen_prod_matrix). The second macro call (%colab_filter) fails with an obvious error message about the non-existent variable "f1:" -- I guess that variable POL should be used in the DROP= option.
The code of macro COLAB_FILTER (up to a renamed array) can be found in a 2014 Stackoverflow thread and is supposed to perform matrix multiplication. So, what's the purpose of macro GEN_PROD_MATRIX?
You wrote that you were "getting a wrong result for a bigger matrix." Why don't you share more details about that "wrong result" (and the expected correct result) with us?

Re: Getting product of matrix

rahul88888 — Wed, 06 Apr 2016 11:25:39 GMT

For this matrix the results is correctly displayed across all the outputs however if we increase the number of rows to lets say 50000 or so the results are not correct.

I also checked the matrix code and then found out that the matrix multiplication is not getting dynamic

Re: Getting product of matrix

rahul88888 — Wed, 06 Apr 2016 17:02:44 GMT

I am attaching the excel file and the logic that i am trying to achieve

Thanks

Re: Getting product of matrix

Ksharp — Thu, 07 Apr 2016 01:30:35 GMT

You can run IML code under University Edition. It is totally free.
Can you elaborate How you get Matrix numerator and  denominator ?

Re: Getting product of matrix

rahul88888 — Thu, 07 Apr 2016 01:46:36 GMT

Numerotor is nothing just the sum of the columns which can be obtained also by proc corr /Cov and creating a daataset for the same

while the dinomintor too is just square root of the combinations of the first matrix

Like Square root for user 1 - it is SQR rooot P1*P1, P1*p2 and so on

Re: Getting product of matrix

Ksharp — Thu, 07 Apr 2016 02:47:14 GMT


IML is the best tool for matrix operation.
I really suggest you to use IML code , which is a lot simpler and easier:



proc iml;
x={
1 0 1 1,
0 0 1 1,
0 1 0 1,
1 0 1 0,
1 1 0 0
};

n=ncol(x);
numerator=j(n,n,.);
do i=1 to n;
 do j=1 to n;
  numerator[i,j]=sum(x[,i]#x[,j]);
 end;
end;

vecdiag=vecdiag(numerator);
denominator=j(n,n,.);
do i=1 to n;
 do j=1 to n;
  denominator[i,j]=sqrt(vecdiag[i]#vecdiag[j]);
 end;
end;

divided=numerator/denominator;

Multiplicated=xmult(x,divided);

print numerator ,denominator,divided,Multiplicated;
quit;

Re: Getting product of matrix

rahul88888 — Thu, 07 Apr 2016 02:57:39 GMT

I wish I had IML I got the error though the array of the first matrix is not becoming dynamic, for some reason will look into the same again. However when I try with fixed product and say let's say for 5 products and fix it Then it works perfectly

Re: Getting product of matrix

Ksharp — Thu, 07 Apr 2016 03:14:09 GMT

I don't know what you are talking about.
If you want read data from table as your sample data. You can do this.
But that lead to ERROR information due to P6 value are all zero , and make Matrix Product meaning nothing.



data test1;
input pol $ prod1-prod9;
datalines;
1abc 1 1 0 1 0 0 1 1 0
2abc 0 0 1 1 0 0 1 1 1
3abc 1 1 0 1 0 0 0 0 0
4abc 1 0 1 1 0 0 1 1 0
5abc 0 1 0 1 1 0 0 1 0
6abc 1 1 0 1 0 0 0 0 0
7abc 1 0 1 1 0 0 1 1 0
8abc 0 1 0 1 1 0 0 1 0
9abc 1 1 0 1 0 0 0 0 0
10abc 1 0 1 1 0 0 1 1 0
11abc 0 1 0 1 1 0 0 1 0
12abc 1 1 0 1 0 0 0 0 0
13abc 1 0 1 1 0 0 1 1 0
14abc 0 1 0 1 1 0 0 1 0
15abc 1 1 0 1 0 0 0 0 0
16abc 1 0 1 1 0 0 1 1 0
17abc 0 1 0 1 1 0 0 1 0
;
run;
proc iml;
use test1;
read all var _num_ into x;
close;

n=ncol(x);
numerator=j(n,n,.);
do i=1 to n;
 do j=1 to n;
  numerator[i,j]=sum(x[,i]#x[,j]);
 end;
end;

vecdiag=vecdiag(numerator);
denominator=j(n,n,.);
do i=1 to n;
 do j=1 to n;
  denominator[i,j]=sqrt(vecdiag[i]#vecdiag[j]);
 end;
end;

divided=numerator/denominator;

Multiplicated=xmult(x,divided);

print numerator ,denominator,divided,Multiplicated;
quit;




OUTPUT:
numerator
11	6	5	11	0	0	6	6	0
6	11	0	11	5	0	1	6	0
5	0	6	6	0	0	6	6	1
11	11	6	17	5	0	7	12	1
0	5	0	5	5	0	0	5	0
0	0	0	0	0	0	0	0	0
6	1	6	7	0	0	7	7	1
6	6	6	12	5	0	7	12	1
0	0	1	1	0	0	1	1	1
denominator
11	11	8.1240384	13.674794	7.4161985	0	8.7749644	11.489125	3.3166248
11	11	8.1240384	13.674794	7.4161985	0	8.7749644	11.489125	3.3166248
8.1240384	8.1240384	6	10.099505	5.4772256	0	6.4807407	8.4852814	2.4494897
13.674794	13.674794	10.099505	17	9.2195445	0	10.908712	14.282857	4.1231056
7.4161985	7.4161985	5.4772256	9.2195445	5	0	5.9160798	7.7459667	2.236068
0	0	0	0	0	0	0	0	0
8.7749644	8.7749644	6.4807407	10.908712	5.9160798	0	7	9.1651514	2.6457513
11.489125	11.489125	8.4852814	14.282857	7.7459667	0	9.1651514	12	3.4641016
3.3166248	3.3166248	2.4494897	4.1231056	2.236068	0	2.6457513	3.4641016	1
divided
1	0.5454545	0.6154575	0.8043997	0	.	0.6837635	0.522233	0
0.5454545	1	0	0.8043997	0.6741999	.	0.1139606	0.522233	0
0.6154575	0	1	0.5940885	0	.	0.9258201	0.7071068	0.4082483
0.8043997	0.8043997	0.5940885	1	0.5423261	.	0.6416889	0.8401681	0.2425356
0	0.6741999	0	0.5423261	1	.	0	0.6454972	0
.	.	.	.	.	.	.	.	.
0.6837635	0.1139606	0.9258201	0.6416889	0	.	1	0.7637626	0.3779645
0.522233	0.522233	0.7071068	0.8401681	0.6454972	.	0.7637626	1	0.2886751
0	0	0.4082483	0.2425356	0	.	0.3779645	0.2886751	1


You notice that P6 are all missing, thus lead to Matrix 'divided' has missing value and can't perform matrix product.

Re: Getting product of matrix

rahul88888 — Thu, 07 Apr 2016 03:20:00 GMT

I accept what you said however I don't have IML on the system then I am working on so will have to try out other ways possible .. I am looking into my code to find the flaws and actually I am trying to get the same answer as you have got . For this dataset my code will work without having to use IML don't know why its failing in a large dataset .. let me try

Re: Getting product of matrix

rahul88888 — Thu, 07 Apr 2016 03:20:36 GMT

Thanks very much for your patience to look through

Re: Getting product of matrix

Ksharp — Thu, 07 Apr 2016 03:54:10 GMT

I am NOT  going to do it by data step. That is a very nasty work.
Maybe that is loss of precision problem.
When you have some very very small value ,that would generate some error value.


Example :

a = 1e13;
b = 1e13;
c = 100 * a;
a = a+1;
x = c || a || b || c;
y = c || a || (-b) || (-c);
z = xmult(x,y`); / * correct answer * /
print z [format=16.0];
wrong = x * y`; / * loss of precision * /
print wrong [format=16.0];
Figure 24.433 Extended-Precision Multiplication
z
20000000000001
wrong
0

Re: Getting product of matrix

rahul88888 — Thu, 07 Apr 2016 06:14:06 GMT

Yes well this may be the case however I am trying in data step again and storing it as a function and then run the same to get the product of matrix.