Date: Wed, 12 Jan 2011 17:13:16 -0800
Reply-To: Dale McLerran <stringplayer_2@YAHOO.COM>
Sender: "SAS(r) Discussion" <SAS-L@LISTSERV.UGA.EDU>
From: Dale McLerran <stringplayer_2@YAHOO.COM>
Subject: Data step code for matrix inversion
Content-Type: text/plain; charset=us-ascii
On Monday, 1/10/11, I posted code for obtaining the inverse of
a square matrix. (See
http://listserv.uga.edu/cgi-bin/wa?A2=ind1101b&L=sas-l&F=&S=&P=15948)
A feature that falls out of the method of computing the
matrix inverse (by diagonalizing the input matrix without any
scaling of the input matrix until the very last step of the
algorithm) is that the product of the diagonal elements of the
diagonalized matrix is the determinant of the input matrix. The
cost of producing the determinant is, then, virtually nil once
the matrix has been diagonalized. Thus, I am posting a revised
version of the macro which adds the determinant as output
generated by the macro.
I would also note that in the previous version, I modified the
input matrix to fill in 0 values when the input matrix had missing
values. There are three problems with that code, all fixed here.
First, it should be standard programming practice not to modify
the input data unless one is expecting modified input data. The
code has been modified so that the input matrix returns unchanged.
Second, the algorithm that was employed did not always work the
way that one would expect. If one had on input
| 3 . 1 |
| . . . |
| 1 . 3 |
the value of the input matrix which was returned by the macro
was
| 0 0 0 |
| 0 0 0 |
| 0 0 3 |
As indicated previously, the current code does not even modify
the input matrix. Rather, it modifies a copy of the input
matrix and returns for the copy of the input matrix the values
| 3 0 1 |
| 0 0 0 |
| 1 0 3 |
The third modification to the macro also has to do with how the
scratch matrix is initialized if there are missing values.
Consider the situation where the value in row 1, column 2 is
missing and no other values in the same row/column are missing.
What is an appropriate value for that cell if we were to fill
in some value? It really is impossible to determine. So the
revised macro issues an error message identifying the problem
and then issues an ABORT statement to cease execution of the
data step or procedure.
The code below is not the complete code for computing the inverse
of a matrix. The MAT_GINV macro shown below calls another macro
which diagonalizes the input matrix (MAT_BIDIAG). The MAT_BIDIAG
macro is found in the previous post (see above). The previous
post also specifies properties of the inverse matrix which
is produced (when the input matrix is not full rank, it is a
generalized inverse with very specific properties) and also
had code to demonstrate how to use the macro.
Dale
---------------------------------------
Dale McLerran
Fred Hutchinson Cancer Research Center
mailto: dmclerra@nFoHsCpRaCm.org
Ph: (206) 667-2926
Fax: (206) 667-5977
---------------------------------------
%macro mat_ginv(input=, inverse=, dim=, det=det);
/**************************************************/
/* This macro computes the inverse of any general */
/* square matrix by solving a set of linear */
/* equations (like the SOLVE() function in IML.) */
/* */
/* Arguments: */
/* INPUT - name of a square data step array */
/* which contains the matrix that */
/* needs to be inverted */
/* INVERSE - name of a square data step array */
/* which will hold the inverse of */
/* the input matrix */
/* DIM - the dimension of the input matrix */
/* DET - the name of a scalar variable */
/* which will hold the determinant */
/* of the input matrix. If the */
/* input matrix is of rank p<&DIM, */
/* the value of DET is the */
/* determinant of the rank p */
/* submatrix which results from */
/* censoring rows and columns of the */
/* input matrix which are linearly */
/* dependent on previous rows and */
/* columns where we work from row 1, */
/* column 1 to row &DIM, column &DIM */
/**************************************************/
%global ginv_invoke;
/* Macro variable ginv_invoke is a counter for the */
/* number of times that this macro has been invoked. */
/* The purpose of this macro variable is to allow */
/* declaration of a scratch array by the MAT_GINV */
/* macro. By naming the array based on the number */
/* of times the MAT_GINV macro has been invoked */
/* within a single SAS session, it is possible to */
/* declare multiple scratch arrays in a single data */
/* step and for those scratch arrays to each have */
/* dimension necessary for the inversion problem. */
%if &ginv_invoke=%str() %then
%let ginv_invoke=1;
%else
%let ginv_invoke=%eval(&ginv_invoke+1);
%let invoke=&ginv_invoke;
array scratch_&invoke [&dim,&dim] _temporary_;
/* Operate on a copy of the original input matrix */
/* Also, initialize the inverse matrix as I(&DIM) */
do i=1 to &dim;
do j=1 to &dim;
scratch_&invoke[i,j] = &input[i,j];
&inverse[i,j] = (i=j);
end;
end;
/* Convert an entire row/column to 0 */
/* if the diagonal entry is missing. */
do i=1 to &dim;
if scratch_&invoke[i,i]=. then
do j=1 to &dim;
scratch_&invoke[i,j]=0;
scratch_&invoke[j,i]=0;
end;
end;
/* If there are still missing values in the input */
/* matrix, then the matrix is ill defined. Give */
/* an error message and generate the inverse as */
/* missing when ill defined. Otherwise, continue */
/* with matrix inversion. */
_continue=1;
do i=1 to &dim;
do j=1 to &dim;
if scratch_&invoke[i,j]=. then do;
_continue=0;
miss_row=i;
miss_col=j;
end;
end;
end;
if _continue=0 then do;
do i=1 to &dim;
do j=1 to &dim;
&inverse[i,j]=.;
end;
end;
put "ERROR: Matrix &input has a missing value in row " miss_row " and column " miss_col ". Missing values"/
" in the input matrix do not allow completion of the matrix inverse unless" /
" the diagonal entry at &input[" miss_row "," miss_row "] or &input[" miss_col "," miss_col "] is also missing" /
" (in which case the entire row and column with missing value" /
" at the diagonal entry is converted to all zero values.";
abort nolist;
end;
else do;
/* Make the input matrix lower diagonal. Operations performed */
/* on the input matrix to make it lower diagonal are also */
/* performed on the matrix which becomes the inverse matrix. */
%mat_bidiag(input=scratch_&invoke, supp=&inverse,
lower_upper=lower, dim=&dim);
/* Make the lower diagonal matrix also an upper diagonal */
/* matrix. Again, perform the same operations on the */
/* matrix which becomes the inverse matrix. */
%mat_bidiag(input=scratch_&invoke, supp=&inverse,
lower_upper=upper, dim=&dim);
/* Scale so all non-zero elements of diagonal matrix are 1. */
/* Perform the same scaling on all rows of the inverse matrix. */
/* Actually, we don't care about actually creating the identity */
/* matrix, so operations are only performed on the inverse mat. */
/* Note that the product of the entries of the diagonal matrix */
/* returned from the bidiag macro is the determinant of the */
/* rank p<=&DIM portion of the input matrix. */
&det=0;
do i=1 to &dim;
if scratch_&invoke[i,i]^=0 then do;
if &det=0 then &det=scratch_&invoke[i,i];
else &det = &det*scratch_&invoke[i,i];
do j=1 to &dim;
&inverse[i,j] = &inverse[i,j] / scratch_&invoke[i,i];
end;
end;
end;
end;
%mend;
