An easy way to generate a vector of letters

A little-known but useful feature of SAS/IML 12.3 (which was released with SAS 9.4) is the ability to generate a vector of lowercase or uppercase letters by using the colon operator (:).

Many SAS/IML programmers use the colon operator to generate a vector of sequential integers:

proc iml;
x = 1:5;     /* increasing sequence {1 2 3 4 5} */
y = 9:5;     /* decreasing sequence {9 8 7 6 5} */

The colon operator is also useful for generating variable names that start with a common prefix and are enumerated consecutively:

varNames = "var1":"var5";  /* {"var1" "var2" "var3" "var4" "var5"} */

The new feature of the colon operator is the ability to generate a sequence of consecutive English letters:

lower = "a":"z";    /* 26 lowercase letters */
upper = "A":"Z";    /* 26 uppercase letters */
s = "p":"k";        /* {"p" "o" "n" "m" "l" "k"} */
u = "U":"Z";        /* {U V W X Y Z} */

This feature is useful when you want to quickly name the levels of a categorical variable. For example, the following SAS/IML statements sample with replacement from four categories named A, B, C, and D:

call randseed(1234);
x = sample("A":"D", 10); /* 10 draw with replacement from {A B C D} */
print x;

Can you think of a way to use this feature in your work? Leave a comment.

Post a Comment

Using associativity can lead to big performance improvements in matrix multiplication

In a previous post, I stated that you should avoid matrix multiplication that involves a huge diagonal matrix because that operation can be carried out more efficiently. Here's another tip that sometimes improves the efficiency of matrix multiplication: use parentheses to prevent the creation of large matrices.

Matrix multiplication is associative, which means that a matrix expression like Z = A*B*C can be computed as Z = (A*B)*C or as Z = A*(B*C). The answer is the same, but sometimes the dimensions of the matrices are such that one way is more efficient than the other. Because matrix languages such as the SAS/IML language evaluate matrix expressions from left to right, the expression Z = (A*B)*C is equivalent to not using parentheses at all.

When do parentheses help? Parentheses help when the product A*B is large, but B*C is small. The following SAS/IML statements compute the same product, but the computation that uses parentheses is much faster because it avoids forming a huge intermediate result:

/* Compute Z = A*B*C where
   A is N x p matrix
   B is p x N matrix
   C is N x k matrix, and N is much large than p and k */
proc iml;
N=10000;  p = 500;  k = 300;
A = j(N,p,1); B = j(p,N,1); C = j(N,k,1);
t0 = time();
Z1 = A*B*C;
T1 = time()-t0;
t0 = time();
Z2 = A*(B*C);
T2 = time()-t0;
print T1 T2;

This computation was done in SAS 9.4, which includes support for multithreaded matrix multiplication in SAS/IML. If you are running SAS 9.3, use a smaller value for N.

You can see that the first computation required about 30 seconds, whereas the second required less than a second. Both expressions compute the same N x k matrix. Why is the second computation so much faster?

  • The first computation, must compute the temporary matrix A*B, which is an N x N matrix. After that huge matrix is formed, it is used to multiply C.
  • The second computation computes the temporary matrix B*C, which is only a small p x k matrix. That matrix is premultiplied by the medium-sized matrix A to form the final result.

So in the second computation there is no "inflation" of the matrix sizes. Each intermediate matrix is the same size or smaller than the original matrices. In contrast, the first computation forms a huge matrix that is much larger than any of the original matrices.

So next time you need to compute a product of several matrices, think about the dimensions of those matrices and choose an order for the multiplication that will keep intermediate matrices small. The effect can be dramatic!

Post a Comment

Never multiply with a large diagonal matrix

I love working with SAS Technical Support because I get to see real problems that SAS customers face as they use SAS/IML software. The other day I advised a customer how to improve the efficiency of a computation that involved multiplying large matrices. In this article I describe an important efficiency tip: Never multiply with a diagonal matrix.

Here is the scenario. The customer needed to compute the matrix Z, which is the symmetric matrix product
     Z = W1/2 B R B′ W1/2

  • W = diag(d) is an N x N diagonal matrix
  • B is an N x p matrix
  • B′ is the transpose of B
  • R is a p x N symmetric matrix. (The symmetry of R isn't exploited in this article.)
In the customer's scenario, N and p were large. Roughly, N = 10,000 and p = 600. The matrix computation was taking a long time, and because the computation was inside a simulation loop, the entire program required many hours to run.

The brute force approach

The customer implemented the formula in the natural way. Let's time how long the straightforward computation takes. (I am using SAS 9.4, which uses multithreaded matrix multiplication. If you are using SAS 9.3, you might want to use N = 5000.) Since the contents of the matrices don't matter, I'll create random elements. As I've shown in a previous blog post, you can use the SQRVECH function to create the symmetric matrix R:

proc iml;
N = 10000; p = 600;
/* define matrices */
d = j(N,1,1);               /* d is N x 1 */
B = j(N,p,1);               /* B is N x p */
v = j(p*(p+1)/2, 1, 1);     /* allocate vector */
/* fill with random uniform numbers in (0,1) */
call randgen(d,"Uniform"); call randgen(B,"Uniform"); call randgen(v,"Uniform");
R = sqrvech(v);             /* create symmetric p x p matrix */
/* straightforward (but slow) computation */
t0 = time();
W = diag(d);                /* N x N diagonal matrix */
Z1 = sqrt(W) * B * R * B` * sqrt(W);
T1 = time() - t0;

On my computer, the naive computation with N = 10000 takes about 24 seconds. I think we can do better with a few small modifications.

Never multiply with a diagonal matrix

The time required to compute this matrix expression can be dramatically shortened by implementing the following improvements:

  • W is a diagonal matrix. Therefore computation sqrt(W) * B multiplies the ith row of B by the ith element of the diagonal of W1/2. You can compute this expression more efficiently by using elementwise multiplication (#) operator, as I showed in an article about converting a correlation matrix into a covariance matrix. The simpler expression is sqrt(d) # B, which also avoids forming the huge N x N diagonal matrix, W, and avoids taking the square-root of N2 elements, most of which are 0.
  • The expression sqrt(W) * B appears twice. The expression appears at the beginning of the formula, and the transpose of the expression appears at the end of the formula. Whenever you see a computation repeated twice, you should consider creating a matrix to hold the intermediate result, such as C = sqrt(d) # B.

If you implement these two improvements, the computation executes much quicker. On my computer it now takes less than a second:

free W;                               /* release the W memory */
/* avoid forming diag(d) and store temporary result */
t0 = time();
C = sqrt(d) # B;
Z2 = C * R * C`;
T2 = time() - t0;
print T1 T2;

When you use a huge N x N diagonal matrix to multiply B, most of the time is spent multiplying the off-diagonal elements, which are zero. The naive approach multiplies (and adds) about 100 million zeros! The elementwise multiplication does not multiply any zeros. Getting rid of the diagonal matrix makes a major difference in the speed of the computation, and leads to the following efficiency tip:
Tip: Never, ever, multiply with a large diagonal matrix! Instead, use elementwise multiplication of rows and columns.

Specifically, if d is a column vector:

  • Instead of diag(d) * A, use d # A to multiply the ith row of A by the ith element of d.
  • Instead of A * diag(d), use A # d` to multiply the jth column of A by the jth element of d.
Do you have any useful tips for computing with large matrices? Leave a comment.

Post a Comment

Tips for concatenating strings in SAS/IML

Last week, as part of an article on how spammers generate comments for blogs, I showed how to generate random messages by using the CATX function in the DATA step. In that example, the strings were scalar quantities, but you can also concatenate vectors of strings in the SAS/IML language. However, there are some problems that need to be handled when concatenating vectors. In this post I describe the following:
  • how to get rid of spaces when you concatenate strings
  • how to insert spaces (or other delimiters) between strings

A canonical application of concatenation is combining the first and last names of individuals to form the full name. For example, the following SAS/IML statements define vectors that contains the first and last names of three famous mathematicians. You can use the SAS/IML CONCAT function or the string concatenation operator (+) to concatenate the names. This example explicitly concatenates a space between the names:

proc iml;
first = {"C."    "Isaac"  "Leonhard"};
last  = {"Gauss" "Newton" "Euler"};
name = first + " " + last;   /* concatenate with space between first and last */
print name;

As I mentioned in my article on how to understand SAS/IML character vectors, there are actually several blank characters between the first and last names of Gauss and Newton. If you use the PrintChars module from my previous post, you see the following:

/* define or load the PrintChars module here... */
run PrintChars(name);

These blank characters appear because the array of first names (first) is a character array in which all elements have the same length. In this case, all elements have length 8, which is the number of characters in the longest name, "Leonhard." Shorter strings are padded with trailing blanks.

Vectorized approach to trim blanks

Although the SAS/IML language and the SAS DATA step language are similar, string concatenation in SAS/IML has some complexities that are not present in the DATA step. In the DATA step, you can use the TRIM function (or the STRIP function) to get rid of blanks. Unfortunately, when you apply these functions to a matrix, they don't solve the problem. The matrix that is returned by trim(first) is exactly the same as first because after TRIM strips off the trailing blanks, the trimmed strings are assembled into a matrix of length 8, which re-adds the trailing blanks!

So how can you get rid of the trailing blanks? One vectorized approach is to use the RIGHT function to right-align the first names prior to concatenating them with the last names. Of course, that will result in strings with leading spaces, so you then need to use the LEFT function to get rid of leading blank. This approach gets messy when you are concatenating many vectors.

I scratched my head over this problem for a long time. For a while I even abandoned trying to use a vectorized approach; I just iterated over the elements of the vectors and concatenated scalar quantities. (Mea culpa!) Then one day I remembered that you can call Base SAS functions from SAS/IML and pass in matrices as arguments. Can the CATX function, which solves the problem for scalar quantities, also solve the concatenation problem for vector quantities? Let's see:

name = catx(" ", first, last);   /* insert space delimiter between names */
run PrintChars(name);

Success! The concatenated strings have trailing blanks, and in every case the first name is separated from the last name by exactly one space. Once again, Base SAS functions help me to solve a problem that involves vectors! From now on I will use the CATX function to concatentate vectors of strings when I want to insert a space between strings.

The first argument to the CATX function specifies the delimiter that is inserted between strings. You can use that argument to insert commas, slashes, or any other delimiters between strings. You can even specify the "null string" (two quotation marks with no space betwen them) to concatenate strings when you do not want to insert a delimiter.

Post a Comment

How to create a string of a specified length in SAS/IML

In my recent post on how to understand character vectors in SAS/IML, I left out an important topic: How can you allocate a character vector of a specified length? In this article, "length" means the maximum number of characters in an element, not the number of elements in a vector.

In the SAS DATA step, you can use the LENGTH statement to specify the maximum number of characters that can be stored in a variable. Thus the length of a variable can be greater than the length of any string that it contains. However, the SAS/IML language does not support the LENGTH statement. This leads to an interesting problem: How can you create a vector in SAS/IML such that each element can contain up to k characters?

This problem comes up in statistical programming because sometimes you need to enforce a limit on the length of a character vector. For example, if you are creating strings that will be used to name SAS variables, the strings must not exceed 32 characters.

There are a couple of possible solutions. In SAS/IML 12.3, which was released with SAS 9.4, you can use the BLANKSTR function to create a blank string with a given number of characters. The following statements create a string that contains 10 space characters. You can then use the J function or the REPEAT function to create a character vector.

proc iml;
/* generate blank strings of an arbitrary size */
maxChars = 10;                /* maximum number of characters to store */
Blanks = BlankStr(maxChars);  /* string of length 10 (SAS/IML 12.3)    */
c = j(1, 5, Blanks);          /* allocate 1x5 vector of blank strings  */

If you do not have SAS/IML 12.3, use the SUBPAD function in Base SAS. The following statement duplicates the functionality of the BlankStr function.

Blanks = subpad(" ", 1, maxChars); /* return string of length 10 */

I'll leave it as an exercise to figure out why this call to the SUBPAD function creates a string with 10 blank characters!

Have you ever faced the problem of generating a string with a specified number of characters? What was your solution?

Post a Comment

This article is actually fastidious: How spammers generate random comments for blogs

Last week Chris Hemedinger posted an article about spam that is sent to SAS blogs and discussed how anti-spam software helps to block spam. No algorithm can be 100% accurate at distinguishing spam from valid comments because of the inherent trade-off between specificity and sensitivity in any statistical test. Therefore, some spam comments slip through the anti-spam filter, and I get the pleasure of reading the comments and deciding whether to allow them to appear on my blog, or whether to manually mark them as spam.

Why spammers submit comments to blogs

When I first started getting spam comments, I wondered what the spammers were trying to achieve. A typical spam comment seems fairly innocuous. Here are some actual spam comments that I have received:

  • Thanks for a marvelous posting! I certainly enjoyed reading it, you can be a great author.
  • Thank you for the good writeup. It in fact was a amusement account it.
  • If all bloggers made good content as you did, the internet will be much more useful than ever before.
  • Wow, this article is actually fastidious.

Yes, some of the grammar and word choices are strange, but these comments are not much different from some legitimate comments that I have received from legitimate readers whose native language is not English. What makes me sure that these are spam comments?

Along with each of these comments, the commenter included a URL link to some web site. The URL in a typical spam comment links to a web site that advertises cheap "name brand" merchandise, Russian brides, or get-rich-quick schemes. The spammers get paid for each link that they can successfully embed somewhere on the web, such as on my blog. As you might know, internet search engines use the number of "incoming links" as a measure of how important a web site is, and therefore how high it should appear in the search results. The goal of the blog spammer is to embed many links in many blog articles so that internet search engines rank their sponsoring web site highly when someone searches for something like "cheap viagra."

The link is not always embedded in the comment itself. When you comment on a blog, you have the option to include your name and to link to your personal web site. In a legitimate comment, the links points to the commenter's blog or business; spammers link to their sponsoring URL.

How spammers create random comments

As you can see from the sample comments, spammers try to construct complimentary but fairly generic message that they can submit regardless of an article's content.

Suppose that a spammer decides to construct the following generic message: Your blog is truly wonderful. He could write a program that submits this comment to a million blogs. However, he would not be very successful because anti-spam software can block this simple attack by applying simple logic: IF the comment is 'Your blog is truly wonderful' AND the URL field is filled in, THEN classify the comment as spam.

To attempt to defeat anti-spam software, spammers randomly generate comments by using synonyms for the nouns, verbs, and adjectives that appear in the comment. For example, synonyms for "blog" include "post" and "article." A synonym for "wonderful" is "marvelous," and so forth. Thus the spammer could modify his spam program to generate random messages according to the following grammatical template: Your NOUN is ADJECTIVE SUBJECT_COMPLEMENT. The result is like those Mad Libs® stories you and your sibling used to create on long car rides. You can create a huge number of possible sentences, but most of them sound silly.

It is easy to write a SAS DATA step program that generates comments that fit the grammatical template. The following program uses the RAND function to randomly generate elements of a character array and uses the CATX function to concatenate the elements into a sentence:

data SpamComments(keep=msg);
array noun{4} $12 _temporary_                /* nouns */
array adj{7}  $12 _temporary_                /* adjectives */
   (" ","simply","actually","truly","sincerely","honestly","really");
array sc{6}   $12 _temporary_                /* subject complements */
sp = " ";                                    /* delimiter between words */
call streaminit(12345);
length msg $ 120;
do i = 1 to 20;                              /* generate 20 messages   */
   noun_i = ceil(dim(noun)*rand("uniform")); /* random index into noun array */
   adj_i  = ceil(dim(adj)*rand("uniform"));
   sc_i   = ceil(dim(sc)*rand("uniform"));
   msg = catx(sp, "Your", noun[noun_i], "is", adj[adj_i], sc[sc_i]);
proc print; run;

As shown by the output on the left, these randomly generated comments can be amusing. It is often obvious when spammers use a thesaurus to automatically generate dozens of synonyms for each word in their message template. Words have subtle connotations; they cannot be mixed and matched like a verbal closet of Garanimals®. I call comments like this "grammaticals."

Do you have any experience dealing with spammers? Share your experience by leaving a comment. I'm sure it will be fastidious. If all readers make astonishing comments as you do, the web will be much more useful than ever before.

Post a Comment

Blanks and lengths: Understanding SAS/IML character vectors

SAS programmers are probably familiar with how SAS stores a character variable in a data set, but how is a character vector stored in the SAS/IML language?

Recall that a character variable is stored by using a fixed-width storage structure. In the SAS DATA step, the maximum number of characters that can be stored in a variable is determined when the variable is initialized, or you can use the LENGTH statement to specify the maximum number of characters. For example, the following statement specifies that the NAME variable can store up to 10 characters:

data A;
length name $ 10;   /* declare that a variable stores 10 characters */

The values in a character variable are left aligned. That is, values that have fewer than 10 characters are padded on the right with blanks (space characters).

SAS/IML character vectors

The same rules apply to character vectors in the SAS/IML language. A vector has a "length" that determines the maximum number of characters that can be stored in any element. (In this article, "length" means the maximum number of characters, not the number of elements in a vector.) Elements with fewer characters are blank-padded on the right. Consequently, the following two character vectors are equivalent. :

proc iml;
c  = {"A",      "B   C",  "  XZ",   "LMNOPQ"}; /* length set at initialization */
c2 = {"A     ", "B   C ", "  XZ  ", "LMNOPQ"}; /* all strings have length 6 */
if c=c2 then print "Character vectors are equal";
else print "Character vectors are not equal";

You can determine the maximum number of characters that can be stored in each element by using the NLENG function in SAS/IML. You can also discover the number of characters in each element of a vector (omitting any blank padding) by using the LENGTH function, as follows:

N = nleng(c);
trimLen = length(c);
print N trimLen c;

In this example, each element of the vector c can hold up to six characters. If you write the c variable to a SAS data set, the corresponding variable will have length 6. However, if you trim off the blanks at the end of the strings, most elements have fewer than six characters. Notice that the LENGTH function counts blanks at the beginning and the middle of a string but not at the end, so that the string "  XZ" counts as four characters.

Where are the blanks?

Notice that the ODS HTML destination is not ideal for visualizing blanks in strings. In HTML, multiple blank characters are compressed into a single blank when the string is rendered, so only one space appears on the displayed output. If you need to view the spaces in the strings, use the ODS LISTING destination, which uses a fixed-width font that preserves spaces. Alternatively, the following SAS/IML function prints each character (not including trailing blanks):

/* convert a string to a row vector of single characters (uses SAS/IML 12.1) */
start Str2Vec(s);
   return (substr(s, 1:length(s), 1));  /* row vector of characters */
/* print characters of all strings in a vector */
start PrintChars(v);
   L = length(v);     /* characters per name, not counting trailing blanks */
   do i = 1 to ncol(L)*nrow(L);
     c = char(1:L[i], 2);
     print (Str2Vec(v[i]))[colname=c];  /* print individual letters */
run PrintChars(c);

I think the Str2Vec function is very cool. It uses a feature of the SUBSTR function in SAS/IML 12.1 to convert a string into a vector of characters. The PrintChars function simply calls the Str2Vec function for each element of a character matrix and prints the characters with a column header. This makes it easy to see each character's position in a string.

This article provides a short overview of how strings are stored inside SAS/IML character vectors. For more details about SAS/IML character vectors and how you can manipulate strings, see Chapter 2 of Statistical Programming with SAS/IML Software.

Post a Comment

Simulating from the inverse gamma distribution in SAS

While at a conference recently, I was asked whether it was possible to use SAS to simulate data from an inverse gamma distribution. The SAS customer had looked at the documentation for the RAND function and did not see "inverse gamma" listed among the possible choices.

The answer is "yes." The distributions that are listed for the RAND function (and the RANDGEN subroutine in SAS/IML) are building blocks that enable you to build other less common distributions. I discuss this in Chapter 7 of Simulating Data with SAS where I show how to simulate data from a dozen distributions that are not directly supported by the RAND function.

When I want to simulate data from a distribution that is not directly supported by the RAND function, I first look at the documentation for the MCMC procedure, which lists the definitions for more than 30 distributions. If I can't find the distribution there, I walk to the library to browse the ultimate reference guide: the two volumes of Continuous Univariate Distributions by Johnson, Kotz, and Balakrishnan (2nd edition, 1994). The Wikipedia article on the inverse gamma distribution uses the same parameterization as these sources, although that is not always the case. (In fact, the Wikipedia article includes a link to an alternative parameterization of the inverse gamma distribution as a scaled inverse chi-squared distribution.)

The reference sources indicate that it is trivial to generate data from the inverse gamma distribution. You first draw x from the gamma distribution with shape parameter a and scale 1/b. Then the reciprocal 1/x is a draw from the inverse gamma distribution with shape parameter a and scale b. You can do this in the DATA step by using the RAND function, or in the SAS/IML language by using the following program:

proc iml;
call randseed(1);
a= 3;           /* shape parameter */
b= 0.5;         /* scale for inverse gamma */
x = j(1000,1);  /* 1000 draws */
call randgen(x, "Gamma", a, 1/b); /* X ~ gamma(a,1/b) */
x = 1/x;        /* reciprocal is inverse gamma */

The histogram to the left shows the distribution of 1000 draws from the inverse gamma distribution with parameters a=3 and b=0.5. The PDF of the inverse gamma distribution is overlaid on the histogram. (For details of this technique, see the article "How to overlay a custom density on a histogram in SAS.") The simulated data extend out past x=5, but the histogram and PDF have been truncated at x=2 to better show the density near the mode b/(1+a) = 1/8.

The lesson to learn is that it is often straightforward to use the elementary distributions in SAS to construct more complicated distributions. There are infinitely many distributions, so no software documentation can ever explicitly list them all! However, you can often use elementary arithmetic operations to express a complicated distribution in terms of simpler distributions.

Post a Comment

How much RAM do I need to store that matrix?

Dear Rick,
I am trying to create a numerical matrix with 100,000 rows and columns in PROC IML. I get the following error:
(execution) Unable to allocate sufficient memory.
Can IML allocate a matrix of this size? What is wrong?

Several times a month I see a variation of this question. The numbers vary. Sometimes it is 250,000 or one million rows. Sometimes the matrix is not square, but it is always big.

The IML procedure holds all matrices in RAM, so whenever I see this question I compute how much RAM is required for the specified matrix. Each element in a double-precision numerical matrix requires eight bytes. Therefore, if you know the size of a matrix, you can write a simple formula that computes the gigabytes (GB) of RAM required to hold the matrix in memory.

There are two definitions of a gigabyte. Disk drives and storage devices use 1 GB to mean 109 bytes. However, the historical definition in computer science is 10243 = 230 bytes, which is 1.07 x 109. These competing definitions are sufficiently close that you don't usually have to worry about the difference. For back-of-the envelope computations I use the simpler formula: an r x c matrix double-precision matrix requires r*c*8/109 gigabytes of RAM.

Because RAM is traditionally measured by using the second definition, the following SAS/IML module accurately calculates the number of gigabytes of RAM (using the computer science definition) required to hold a matrix of a given dimension:

proc iml;
/* Compute gigabytes (GB) of RAM required for matrix with r rows and c columns */
start HowManyGigaBytes(Rows, Cols);
   GB = 8 # Rows#Cols / 2##30;  /* 1024##3 bytes in a gigabyte */
   Fit2GB = choose(GB <= 2, "Yes", " No");
   print Rows[F=COMMA9.] Cols[F=COMMA9.] 
         GB[F=6.2] Fit2GB[L="Fits into 2GB"];
/* test:   rows   cols */
sizes = {250000   1000,
          10000  10000,
          16000  16000,
          40000  40000,
         100000 100000 };
run HowManyGigaBytes(sizes[,1], sizes[,2]);

You can see from the table that a square matrix with dimension 100,000 requires 74.5 GB of RAM when stored as a dense matrix. A typical "off-the-shelf" laptop or desktop computer might have 4 or 8 GB of RAM, but of course some of that is used by the operating system.

Usually you need to hold several matrices in RAM in order to add or multiply them together. The last column in the table indicates whether the specified matrix dimensions can fit into 2 GB of RAM, which is a convenient proxy for the practical question, "how big can my matrices be if I want to operate with several at once." Of course, the definition of "several" will depend on whether your computer has 8 GB, 16 GB, or more RAM. For square matrices, a matrix of size 16,000 will fit into 2 GB of RAM.

Because 2 GB is a useful benchmark and because not all matrices are square, the following statements show the dimensions of some non-square matrices that also fit into 2 GB:

rows = 1000 * ({5 10 16} || do(20, 100, 20));
colsPerGB = 2##30 / 8 / rows;
cols = floor( 2 * colsPerGB ); /* columns for 2 GB */
print (rows // cols)[F=comma10. R={"Rows" "Cols"} 
                     L="Matrices that Fit into 2 GB RAM"];

You can compute more values in this fashion and use PROC SGPLOT to create a graph that shows the dimensions of matrices that fit into 2 GB of RAM. The vertical axis of the following image has been truncated at 25,000 columns. (Click to enlarge.) The full graph is symmetric about the identity line because a matrix and its transpose contain the same number of columns.


In conclusion, it is useful to be able to calculate how many gigabytes of storage a matrix requires. In SAS/IML, if you intend to compute with several matrices, they must collectively fit into RAM. I have provided some simple computations that enable you to quickly check whether your matrices fit into 2 GB of RAM. You can, of course, modify the program to find matrices that fit into 1 GB, 0.5 GB, or any other size.

Post a Comment

The inverse of the Hilbert matrix

Just one last short article about properties of the Hilbert matrix. I've already blogged about how to construct a Hilbert matrix in the SAS/IML language and how to compute a formula for the determinant. One reason that the Hilbert matrix is a famous (some would say infamous!) example in numerical linear algebra is that the inverse matrix is known explicitly and is a matrix of integers. If H is an n x n Hilbert matrix, then the value of the (i,j)th element of the matrix is given by the following formula:


How could you compute the Hilbert inverse in a matrix language such as SAS/IML? The obvious way is to write nested DO loops and call the COMB function in Base SAS to compute the binomial coefficients that appear in the formula, as follows:

proc iml;
n = 5;
/* inverse of Hilbert matrix is a symmetric matrix of large integers */
invH = j(n,n);
do i = 1 to n;
   do j = i to n;
      b1 = comb(n+i-1, n-j);
      b2 = comb(n+j-1, n-i);
      b3 = comb(i+j-2, i-1);
      invH[i,j] = (-1)##(i+j) # (i+j-1) # b1 # b2 # b3##2;
      invH[j,i] = invH[i,j];
print invH;

For small matrices, nested loops are fine. However, you should assume that someday your program will be used to compute a large matrix. Therefore, you should vectorize the construction so that it avoids using double loops.

The COMB function is a scalar function that computes the number of combinations of r elements taken s at a time. In this definition, r and s are scalar values. But remember that from within the SAS/IML language you can pass matrix arguments to Base SAS functions. If you pass in matrix values for r and s, the COMB function will process the matrices elementwise and return a matrix of binomial coefficients.

Consequently, to vectorize the construction of the Hilbert matrix inverse you need to define i and j so that they are matrices instead of scalars. You can use the ROW function in SAS/IML 12.3 software to construct the i matrix and the COL function to construct the j matrix. (If you do not have SAS/IML 12.3 or later, the functions are defined in a previous blog.)

Therefore you can define the inverse of the Hilbert matrix as follows:

/* compute inverse of nxn Hilbert matrix */
invH = j(n,n);
i = row(invH); j = col(invH);    /* matrices */
b1 = comb(n+i-1, n-j);           /* matrix of binomial coefficients */
b2 = comb(n+j-1, n-i);
b3 = comb(i+j-2, i-1);
invH = (-1)##(i+j) # (i+j-1) # b1 # b2 # b3##2;

This examples is applicable to other matrices for which the (i,j)th element is defined by using a formula. The ROW and COL functions enable you to implement the formula and define the matrix without writing any loops over the rows and columns.

Post a Comment