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.

SAS
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);
   end;
    output;
  end;
run;

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.; 
run;

 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

R
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.

n=5
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.