Filling an Upper Triangular Matrix from a Vector

8

The R You Ready blog posed an interesting problem. Essentially, you have a vector that contains n(n+1)/2 elements, and you want to pack those elements into the upper left triangular portion of a matrix. For example, if your data are

proc iml;
/** vector v is given: ncol(v) = n(n+1)/2 for some n **/
v = {6 6 6 6 6 6 5 5 5 5 5 4 4 4 4 3 3 3 2 2 1};  
n = (-1 + sqrt(1 + 8*ncol(v)))/2; /** find n using the quadratic formula **/

then you want to create the matrix M as follows:

        6         6         6         6         6         6
        5         5         5         5         5         0
        4         4         4         4         0         0
        3         3         3         0         0         0
        2         2         0         0         0         0
        1         0         0         0         0         0

In the SAS/IML language (which is similar to the MATLAB software mentioned in the blog), I would write the following statements:

M = j(n,n,0);       /** allocate the matrix and fill with zeros **/
vptr = 1;           /** this variable keeps track of our location in v **/
do i = 1 to n;
   j = n-i+1;       /** length of data for the i_th row **/
   M[i, 1:j] = v[ ,vptr:vptr+j-1]; /** copy relevant data **/
   vptr = vptr + j; /** increment the location variable **/
end;
print M;
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.

8 Comments

  1. Hi,

    I would like to take the opportunity to ask (rather off-topic) if anything can be said about the problems one might encounter when converting a Matlab programme tot SAS/IML.

    We are developing an application (using SAS9.2), and it involves, among other things, the incorporation of an existing algorithm written in Matlab. We haven't yet chosen between:
    A- exporting data to Matlab, run the Matlab code and import the result
    B- recoding the Matlab code in SAS/IML.

    In the long run we prefer B.
    But althoug we have extensive SAS expertise, our IML experience is very limited.
    I do know enough on matrix algebra (I think), and I think I roughly understand what the Matlab code is doing (there is also documentation, and the person that developed the algorithm is available).

    It seems that most Matlab code lines can be translated more or less 1-to-1 to IML code. Or might it be more complex?
    I have asked and googled around on some information but didn't really find anything.

    I wonder if you can say something on this?

    Thanks in advance!

  2. I have a few comments:
    1) Post this question on the SAS/IML Discussion Forum in order to solict advice from other SAS customers: http://support.sas.com/forums/index.jspa

    2) I have done (B) many times. It is my method of choice, but then I know both MATLAB and IML. I have also supervised a graduate student in a large-scale project of converting her advisor's MATLAB code to IML. She converted about 6,800 lines of code over the course of the semester. Yes, it is often 1-to-1 conversion, but it depend on what MATLAB functions are being called.

    3) The choice between (A) and (B) depends on the time you have, the expertise you have, and the complexity of the algorithm you are trying to convert.

    4) I once had bookmarked a "cheat sheet/reference sheet" that had a comparison of MATLAB and IML syntax, but I can no longer locate it!

  3. In Matlab you don't have to write loop, you could write 2 lines

    [i,j]=find(triu(ones(n)))
    M=full(sparse(i,j,v,n,n))

    and you are done. Way more compact

  4. Thanks for the feedback. I completely agree that when a language has a built-in function that can help you avoid a loop, then using it is often more efficient than writing an explicit loop. For this example, the loop simply keeps track of the indices and copies over some data, so it is very fast. As to compactness, if I plan to use this functionality often, I can wrap it up into a module for ease of use.

  5. When I type in Atul's MATLAB code I get

    M =

    6 6 6 5 5 3
    0 6 6 5 4 3
    0 0 6 5 4 3
    0 0 0 5 4 2
    0 0 0 0 4 2
    0 0 0 0 0 1

    which is not the correct answer.

  6. Here is the correct version of Atul's MATLAB code:

    [i,j]=find(fliplr(triu(ones(n))));
    M=full(sparse(i,j,v,n,n))';

    or more intuitively:

    M=zeros(n);
    b=fliplr(triu(~M));
    M(b)=v;
    M=M';

    Avoiding loops when doing operations over matrix elements tends to make things much more intuitive in MATLAB.

  7. In Matlab!
    For example, your vector for some n:
    n = 5;
    v = 1:n*(n+1)/2;
    An your matrix
    M = fliplr(triu(ones(n)));
    M(logical(M)) = v;
    Ciao

Leave A Reply

Back to Top