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];
   end;
end;
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.