SAS/IML programmers know that the VECDIAG matrix can be used to extract the diagonal elements of a matrix. For example, the following statements extract the diagonal of a 3 x 3 matrix:
proc iml; m = {1 2 3, 4 5 6, 7 8 9}; v = vecdiag(m); /* v = {1,5,9} */ |
But did you know that you can also assign the diagonal elements without using a loop? Because SAS/IML matrices are stored in row-major order, the elements on the diagonal of an n x p matrix have the indices 1, p+1, 2p+2, ... np. In other words, the following statements assign the diagonal elements of a matrix:
start SetDiag(A, v); diagIdx = do(1,nrow(A)*ncol(A), ncol(A)+1); A[diagIdx] = v; /* set diagonal elements */ finish; run SetDiag(m, 3:1); /* set diagonal elements to {3,2,1} */ run SetDiag(m, .); /* set diagonal elements to missing */ |
I've used this trick in several past blog posts (for example, in my post about how to compute the log-determinant of a matrix), but recently I saw some SAS/IML code that used loops to set the diagonal elements, as follows:
do i = 1 to ncol(m); m[i,i] = v[i]; /* set diagonal elements (not efficient) */ end; |
As I have said before, to increase efficiency, replace loops over matrix elements with an equivalent vector operation. Assigning the diagonal elements of a matrix by using this trick is another example of how to avoid looping in the SAS/IML language.
4 Comments
Pingback: An efficient way to increment a matrix diagonal - The DO Loop
Pingback: Compute nearest neighbors in SAS - The DO Loop
Pingback: Difference operators as matrices - The DO Loop
Pingback: Define or extract the diagonals of a matrix - The DO Loop