[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
More information about the R-devel
mailing list