[Rd] matrix bands

From: Jeremy David Silver <jds_at_dmu.dk>
Date: Fri, 26 Aug 2011 13:23:43 +0200


Dear R developers,

I was looking for a function analogous to base::diag() for getting and setting bands of a matrix. The closest I could find was Matrix::band(), but this was not exactly what I wanted for two reasons. Firstly, Matrix::band() returns a matrix rather than just the specified band. Secondly, Matrix::band() cannot be used for setting the values for a matrix band.

Setting or getting a matrix band could be of interest for sufficiently many users that you might consider including it in base R. Alternatively, something like this could be included in the Matrix package.

I have included two versions of these functions, a simple and naive version, and a more efficient version. The band index, n, is positive for bands above the diagonal and negative for bands below the diagonal.

Regards,
Jeremy Silver

###############################################

## less clear formulation - more efficient band <- function(x,n = 0){

   dx <- dim(x)
   if(length(dx) != 2L)
     stop("only matrix bands can be accessed")

   max.dx <- max(dx)
   n <- as.integer(n)

   ij <- cbind(i = seq(1,max.dx) - n,

               j = seq(1,max.dx))
   ij <- ij[1 <= ij[,1] & ij[,1] <= dx[1] & 1 <= ij[,2] & ij[,2] <= dx[2],,drop=FALSE]

   if(nrow(ij) == 0)
     stop('cannot access this matrix band')

   x[ij]
}

'band<-' <- function(x,n = 0, value){

   dx <- dim(x)
   if(length(dx) != 2L)
     stop("only matrix bands can be replaced")

   max.dx <- max(dx)
   n <- as.integer(n)

   ij <- cbind(i = seq(1,max.dx) - n,

               j = seq(1,max.dx))
   ij <- ij[1 <= ij[,1] & ij[,1] <= dx[1] & 1 <= ij[,2] & ij[,2] <= dx[2],,drop=FALSE]

   if(nrow(ij) == 0)
     stop('cannot replace this matrix band')

   x[ij] <- value

   x
}

## simple, clear formulation - not very efficient band2 <- function(x, n = 0) {

   x[col(x) - row(x) == as.integer(n)]
}

'band2<-' <- function(x, n = 0, value) {

   x[which(col(x) - row(x) == as.integer(n))] <- value    x
}

## here are some examples to show that it works

## define a test matrix
 > A <- matrix(rnorm(12),3,4)
 > A

            [,1] [,2] [,3] [,4]
[1,] -1.5560200 0.6452762 1.072565 0.1923451
[2,] 0.7940685 1.2441817 1.699486 -0.2998814

[3,] -0.7762252 -0.4824173 -0.981055 -0.9265627

## access some of the bands

 > band(A,1)
[1] 0.6452762 1.6994858 -0.9265627

 > band(A,-2)
[1] -0.7762252

 > band(A,2)
[1] 1.0725649 -0.2998814

## set one of the bands

 > band(A,2) <- 2:1
 > A

            [,1] [,2] [,3] [,4]
[1,] -1.5560200 0.6452762 2.000000 0.1923451
[2,] 0.7940685 1.2441817 1.699486 1.0000000

[3,] -0.7762252 -0.4824173 -0.981055 -0.9265627

## another example - a single column

 > A <- matrix(1:10)
 > A

       [,1]

  [1,]    1
  [2,]    2
  [3,]    3
  [4,]    4
  [5,]    5
  [6,]    6
  [7,]    7
  [8,]    8
  [9,]    9

[10,] 10

 > band(A,0)
[1] 1

 > band(A,1)
Error in band(A, 1) : cannot access this matrix band  > band(A,-1)
[1] 2

 > band(A,-5)
[1] 6

## compare the results from the two versions of the function

 > for(i in -2:3){print(band(A,i));print(band2(A,i))}
[1] -0.7762252
[1] -0.7762252
[1] 0.7940685 -0.4824173
[1] 0.7940685 -0.4824173
[1] -1.556020 1.244182 -0.981055
[1] -1.556020 1.244182 -0.981055

[1] 0.6452762 1.6994858 -0.9265627
[1] 0.6452762 1.6994858 -0.9265627
[1] 2 1
[1] 2 1
[1] 0.1923451
[1] 0.1923451

## show that the naive version is very slow for large matrices

 > N <- 1e4
 > M <- matrix(0,N,N)
 > system.time(band(M,2))

    user system elapsed
   0.005 0.003 0.007
 > system.time(band2(M,2))

    user system elapsed
  18.509 2.121 20.754



R-devel_at_r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel Received on Fri 26 Aug 2011 - 11:41:24 GMT

This quarter's messages: by month, or sorted: [ by date ] [ by thread ] [ by subject ] [ by author ]

All messages

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.2.0, at Fri 26 Aug 2011 - 13:10:23 GMT.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-devel. Please read the posting guide before posting to the list.

list of date sections of archive