```Date: Wed, 12 Jan 2011 17:13:16 -0800 Reply-To: Dale McLerran Sender: "SAS(r) Discussion" From: Dale McLerran 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; ```

Back to: Top of message | Previous page | Main SAS-L page