# [Rd] matrix bands

Jeremy David Silver jds at dmu.dk
Fri Aug 26 13:23:43 CEST 2011

```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

```