Here, I demonstrate how to inverse a matrix using SAS/STAT.
There are several ways, invovling procedures such as PROC PRINCOMP, PROC REG, PROC ORTHOREG, each has their advantages and disadvantages.
1. PROC PRINCOMP conduct matrix inversion by relating SVD and Pseduo-Inverse. This is detailed @ here
2. We can also use PROC REG to Inverse a non-singular square matrix in SAS, in the following way:
%let dim=500;
data A;
array _X{&dim} X1-X&dim;
do i=1 to &dim;
do j=1 to &dim;
_X[j]=ranuni(0);
end;
output;
end;
run;
/* patch identity matrix */
data A_expanded;
array _I{&dim} I1-I&dim;
set A;
do i=1 to &dim;
if i=_n_ then _I[i]=1;
else _I[i]=0;
end;
run;
/* matrix inversion using formula: A%*% A^{-1}=I */
proc reg data=A_expanded outest=beta(keep=X1-X&dim) noprint;
model I1-I&dim=X1-X&dim /noint;
run;quit;
PROC REG output
PROC ORTHOREG output
3. More accurate approach requires the use of PROC ORTHREG, a SAS procedure designed for ill-conditioned matrices. Consider the following cases.
data A;
input X1-x5;
datalines;
0 0 1 0 1
1 0 0 1 0
0 1 1 0 1
1 0 0 0 1
0 1 0 1 0
;
run;
/* inv(A)=
1 1 1 0 1
1 0 1 0 0
0 1 1 1 1
1 0 1 0 1
1 1 1 1 1
*/
data A_i;
set A;
array _I{*} I1-I5;
do j=1 to 5;
if j=_n_ then _I[j]=1; else _I[j]=0;
end;
drop j;
run;
proc reg data=A_i outest=Ainv(keep=X1-X5) noprint;
model I1-I5=X1-X5/noint;
run;quit;
proc transpose data=Ainv out=Ainv_t; run;
%macro wrap;
%let n=5;
%do j=1 %to &n;
proc orthoreg data=A_i outest=A_inv_&j.(keep=X1-X&n)
noprint singular=1E-16;
model I&j=X1-X&n /noint;
run;quit;
%end;
data A_inverse;
set %do j=1 %to &n;
A_inv_&j
%end;;
run;
%mend;
%wrap;
proc transpose data=A_inverse out=A_inverse_t;
run;
At last but not least, I would like to emphasis that when the matrix is very ill-conditioned, PROC ORTHOREG will also fail to produce correct answer as well. Consider the Hilbert matrix of dimesion 8-by-8.
%let n=8;
data hilbert;
array _H{*} H1-H&n;
do i=0 to %eval(&n-1);
do k=1 to &n;
_H[k]=1/(i+k);
end;
output;
drop k i;
end;
run;
data hilbert_I;
set Hilbert;
array _I{*} I1-I&n;
do j=1 to dim(_I);
if j=_n_ then _I[j]=1; else _I[j]=0;
end;
drop j;
run;
%macro wrap;
%do j=1 %to &n;
proc orthoreg data=hilbert_i outest=Hilbert_inv_&j.(keep=H1-H&n)
noprint singular=1E-16;
model I&j=H1-H&n /noint;
run;quit;
%end;
data Hilbert_inverse;
set %do j=1 %to &n;
Hilbert_inv_&j
%end;;
run;
%mend;
%wrap;
proc reg data=hilbert_i outest=hilbert_inverse(keep=H1-H&n) noprint;
model I1-I&n=H1-H&n/noint;
run;quit;
