Example 2014.4: Hilbert Matrix

April 14, 2014

[This article was first published on SAS and R, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

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

To leave a comment for the author, please follow the link and comment on their blog: SAS and R.

R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

If you got this far, why not subscribe for updates from the site? Choose your flavor: e-mail, twitter, RSS, or facebook...

Comments are closed.

Search R-bloggers


Never miss an update!
Subscribe to R-bloggers to receive
e-mails with the latest R posts.
(You will not see this message again.)

Click here to close (This popup will not appear again)