Matrix multiplication in SAS (Without PROC IML)

[Jhoony]

Hello,

I am trying to perform a simple matrix multiplication in SAS without PROC IML.

So, the matrix multiplication is in the form of: 

A * transpose(A)

I have found a following code using a simple input-data: 

* Input-data;
data matrix_A;                   
   input a b c;
   datalines;        
4 5 6
7 8 9
10 11 12
;          


* Matrix multiplication;
proc transpose data=matrix_A out=matrix_A_t(drop=_name_);
run;

%macro prod_mat(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;
    end;
run;

/* multiplication*/
data &ou_AB;
    /* new columns for products, equal to number of columns in the 2nd matrix*/
    array p[&B_cols];
    do j=1 to &B_cols;
        p[j]=0;
        set &ou_AB;
        array col _ALL_;
        /* multiply corresponding pairs of columns*/
        do i=&B_cols+2 to &B_cols+1+&B_rows;
            p[j]+col[i]*col[i+&B_rows];
        end;
    end;
    output;
    keep p:;
run;
%mend prod_mat;
%prod_mat(in_A =matrix_A,in_B=matrix_A_t,ou_AB=work.matrix_AA_t)

However, this code is quite slow on larger datasets. 

Is there a way to optimize this code (without using PROC IML)? 

Or is there another way of doing A*transpose(A) matrix multiplication (without PROC IML)? 

[FreelanceReinh]

Hello @Jhoony,

I know that there are routines and methods for matrix operations in

1. PROC FCMP: see CALL TRANSPOSE Routine and CALL MULT Routine for your example
2. PROC DS2: see TRANS Method and MULT Method for your example

but I have to admit that I have little to no experience with them (still on my to-do list ...).

Here's how your example can be computed with PROC FCMP:

proc fcmp;
array A[3,3] (4 5 6 7 8 9 10 11 12);
array A_t[3,3];
array AA_t[3,3];
call transpose(A, A_t);
call mult(A, A_t, AA_t);
put AA_t=;
run;

For the DS2 methods see Example 2: Multiply Two Matrices That Are Read from External Data in the MULT method documentation.

[Ksharp]

data matrix_A;                   
   input a b c;
   datalines;        
4 5 6
7 8 9
10 11 12
;     


proc transpose data=matrix_A out=A(drop=_name_) ;
run;
proc transpose data=A(obs=0) out=vnames ;
run;
proc sql noprint;
create table cartesian as
select a._name_ as a,b._name_ as b  from vnames as a,vnames as b;

select cat(cats('sum(',a,'*',b),') as ',cats(a,'_',b)) into : list separated by ','
 from cartesian;

create table temp as
select &list from A;
quit;
proc transpose data=temp out=temp2;
run;
data temp3;
 set temp2;
 id1=scan(_name_,1,'_');
 id2=scan(_name_,2,'_');
run;
proc transpose data=temp3 out=want(drop=_name_);
by id1;
id id2;
var col1;
run;

[mkeintz]

Once you've done a proc transpose on B, then you can use PROC SCORE to do the multiplication.  PROC SCORE (see https://support.sas.com/documentation/onlinedoc/stat/131/score.pdf ) is primarily meant to generate factor scores, principal component scores, or predicted y values from a set of regression coefficients and a set of regressor variables.

But it can also be used for ordinary matrix multiplication.

data a;
  input idLTR :$1. var1 var2 var3;
datalines;
A 11 12 13
B 21 22 23
C 31 32 33
D 41 42 43
E 51 52 53
run;
data b;
  input var1 var2 var3;
datalines;
1 2 3
10 20 30
100 200 300
run;

data bview / view=bview;
  set b;
  _TYPE_='SCORE';
run;
PROC TRANSPOSE data=bview out=b_t;
  copy _type_;
run;
proc score data=A (rename=(var1=col1 var2=col2 var3=col3))
           score=B_T
           out=WANT; 
   var col1 col2 col3;
   id IDLTR;
run; 

Note you have to rename the VAR names from dataset A, but that's just to match the output of the proc transpose (which might be tweaked to avoid this rename requirement, but I didn't take time to confirm).  The WANT dataset will still get the original dataset A variable names, (due to the _NAME_ variable also produced by proc transpose.

[Tom]

You don't need the RENAME= with such simple sequential variable names.  Just tell PROC TRANSPOSE to use the same prefix when generating variable names so they match those used in dataset A.

proc tranpose data=bview out=b_t prefix=var;
  copy _type_;
run;
proc score data=A score=B_T out=WANT ; 
   var var1-var3;
   id IDLTR;
run; 

[FreelanceReinh · Accepted solution]

---

@Jhoony wrote:

A * transpose(A)

---

This particular matrix product can also be obtained with PROC CORR (after PROC TRANSPOSE):

proc transpose data=matrix_A out=A_t;
run;

ods output sscp=AA_t;
proc corr data=A_t sscp;
var col:;
run;

Edit: The remaining output from PROC CORR can be suppressed with ods select none;.