[Rd] matrix bands

[Jeremy David Silver]

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