The inverse of the Hilbert matrix

0

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:

hilbertmatrixinveqn

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;
t_hilbertinv

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.

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.

Leave A Reply

Back to Top