Monday, April 14, 2014

Example 2014.4: Hilbert Matrix

Rick Wicklin showed how to make a Hilbert matrix in SAS/IML. Rick has a nice discussion of these matrices and why they might be interesting; the value of H_{r,c} is 1/(r+c-1). We show how to make this matrix in the data step and in R. We also show that Rick's method for displaying fractions in SAS/IML works in PROC PRINT, and how they can be displayed in R.

In the SAS data step, we'll use an array to make the columns of the matrix and do loops to put values into each cell in a row and output the row into the data set before incrementing the row value. This is a little awkward, but does at least preserve the simple 1/(r+c-1) function of the cell values. We arrange the approach using a global macro to be more general.
%let n = 5;
data h;
array val [&n] v1 - v&n;
  do r = 1 to &n;
    do c = 1 to &n;
      val[c] = 1/(r+c-1);
To print the resulting matrix, we use the fract format, as Rick demonstrated. Pretty nice! The noobs option in the proc print statement suppresses the row number that would otherwise be shown.
proc print data = h noobs; 
var v1 - v5;
format _all_ fract.; 
 v1    v2     v3     v4     v5

  1    1/2    1/3    1/4    1/5
1/2    1/3    1/4    1/5    1/6
1/3    1/4    1/5    1/6    1/7
1/4    1/5    1/6    1/7    1/8
1/5    1/6    1/7    1/8    1/9

As is so often the case in R, a solution can be generated in one line using nesting. Also as usual, though, its a bit unnatural, and we don't deconstruct it here.
n = 5
1/sapply(1:n,function(x) x:(x+n-1))
A more straightforward approach follows. It's certainly less efficient, though efficiency would seem a non-issue in any application of this code. It's also the kind of code that R aficionados would scoff at. It does have the attractiveness of perfect clarity, however. We begin by defining an empty matrix, then simply loop through the cells of the matrix, assigning values one by one.
h1 = matrix(nrow=n,ncol=n)
for (r in 1:n) {
  for (c in 1:n)
    h1[r,c] = 1/(r+c-1)
To display the fractions, we use the fractions() function in MASS package that's distributed with R.
> library(MASS)
> fractions(h1)
     [,1] [,2] [,3] [,4] [,5]
[1,]   1  1/2  1/3  1/4  1/5 
[2,] 1/2  1/3  1/4  1/5  1/6 
[3,] 1/3  1/4  1/5  1/6  1/7 
[4,] 1/4  1/5  1/6  1/7  1/8 
[5,] 1/5  1/6  1/7  1/8  1/9 

An unrelated note about aggregators: We love aggregators! Aggregators collect blogs that have similar coverage for the convenience of readers, and for blog authors they offer a way to reach new audiences. SAS and R is aggregated by R-bloggers, PROC-X, and statsblogs with our permission, and by at least 2 other aggregating services which have never contacted us. If you read this on an aggregator that does not credit the blogs it incorporates, please come visit us at SAS and R. We answer comments there and offer direct subscriptions if you like our content. In addition, no one is allowed to profit by this work under our license; if you see advertisements on this page, the aggregator is violating the terms by which we publish our work.


Rick Wicklin said...

I think you overlooked the point of my last example. Both SAS and R support the row() and col() functions, so you can construct the Hilbert matrix using nearly identical statements. The SAS/IML statements are:

H = j(5,5, .); /* 5x5 matrix of missing */
H = 1 / (row(H) + col(H) - 1);

The R statements are similar:

H = matrix(NA, 5,5) # 5x5 matrix of missing
H = 1 / (row(H) + col(H) - 1)

Ken Kleinman said...

Yes, I did miss your point, Rick. I don't use IML, so I don't try to understand the IML parts of your posts. Apologies-- I know your focus is IML, but I find great value in the additional content.

The row() and col() functions in R return the row number or column number of each item in a matrix. Thus for the matrix
> matrix(1:4, nrow=2)
[,1] [,2]
[1,] 1 3
[2,] 2 4

we find
> row(matrix(1:4, nrow=2))
[,1] [,2]
[1,] 1 1
[2,] 2 2

This does suggest more elegant and readable one-line solutions than the sapply() version I supplied above.

anspiess said...

Nice post! Also a classical example for which to apply 'outer':

outer(1:5, 1:5, function(x, y) 1/(x + y - 1))


Ken Kleinman said...

Nice, Andrej, thanks.

For those not familiar, in R outer() generates the outer product of two vectors, by default. This is equivalent to the %o% operator:
> 1:2 %o% 1:2
[,1] [,2]
[1,] 1 2
[2,] 2 4
> outer(1:2,1:2)
[,1] [,2]
[1,] 1 2
[2,] 2 4
As Andrej demonstrates, you can supply a function on the vector inputs as an optional argument to outer().

(To get the inner product of two vectors, use the %*% operator:
> 1:2 %*% 1:2
[1,] 5