Creating a tridiagonal matrix

3

I was recently asked how to create a tridiagonal matrix in SAS/IML software. For example, how can you easily specify the following symmetric tridiagonal matrix without typing all of the zeros?

proc iml;
m = {1 6 0 0 0,
     6 2 7 0 0,
     0 7 3 8 0,
     0 0 8 4 9,
     0 0 0 9 5 };

This matrix can be created in two steps. The first step is to use the DIAG function to create a matrix that contains specific values on the diagonal and zeros elsewhere:

/** create diagonal matrix **/
d  = 1:5;
m = diag(d);

The second step relies on the fact that SAS/IML matrices are stored in row-major order. Therefore, the indices of the upper diagonal of m are 2, 8, 14, and 20. In general, for an n x n matrix, the following statements compute the indices of the upper diagonal and assign values to those locations:

n = ncol(m);                /** size of matrix **/
d1 = 6:9;                   /** values for upper diag **/
uppIdx = do(2, n*n, n+1);   /** indices of upper diag **/
m[uppIdx]= d1;              /** set upper diag **/

Similarly, the following statements assign the lower diagonal:

lowIdx = do(n+1, n*n, n+1); /** indices of lower diag **/
m[lowIdx] = d1;             /** set lower diag **/
print m;

Of course, you can also use the indices of the upper or lower diagonal to extract the values, rather than to assign the values.

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.

3 Comments

  1. Hi Rick,

    Need to find out how to calculate the tri-diagonal matrix H for a sample logit[y|X] = theta0 +(h1(X1))`*( theta1) +(h2(X2))`*( theta2)+........ +(hn(Xn))`*( thetan) for natural cubic splines with knots at quantiles?* I'm struggling to get such plus I'm new to iml programmin...

    Kind regards,
    Mxo

Leave A Reply

Back to Top