Assign the diagonal elements of a matrix

3

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.

Share

About Author

Rick Wicklin

Distinguished Researcher in Computational Statistics

Rick Wicklin, PhD, is a distinguished researcher in computational statistics at SAS and is a principal developer of PROC IML and SAS/IML Studio. His areas of expertise include computational statistics, simulation, statistical graphics, and modern methods in statistical data analysis. Rick is author of the books Statistical Programming with SAS/IML Software and Simulating Data with SAS.

Back to Top