# Example 2014.4: Hilbert Matrix

April 14, 2014
By

(This article was first published on SAS and R, and kindly contributed to R-bloggers)

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/51/2    1/3    1/4    1/5    1/61/3    1/4    1/5    1/6    1/71/4    1/5    1/6    1/7    1/81/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 = 51/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=5h1 = 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