# Re: [Rd] matrix bands

From: Martin Maechler <maechler_at_stat.math.ethz.ch>
Date: Fri, 26 Aug 2011 14:08:19 +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.

Yes, but did you not look at help(band) or not look carefully enough ?

--->

```       ‘bandSparse’ for the _construction_ of a banded sparse matrix
directly from its non-zero diagonals.

```

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

well, see above and let us know if you see anything lacking in bandSparse().
Till now we haven't got much feedback about it, and there may well be room for improvement.

Martin Maechler, ETH Zurich and co- maintainer("Matrix")

> 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 <= ij[,2] &
> ij[,2] <= dx,,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 <= ij[,2] &
> ij[,2] <= dx,,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 }

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

>  0.6452762 1.6994858 -0.9265627
```    >> band(A,-2)
```

>  -0.7762252

>> band(A,2)
>  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
```    >> 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
```    >> band(A,1)
```

> Error in band(A, 1) : cannot access this matrix band
```    >> band(A,-1)
```

>  2

>> band(A,-5)
>  6
```    >> for(i in -2:3){print(band(A,i));print(band2(A,i))}
```

>  -0.7762252  -0.7762252  0.7940685 -0.4824173 
> 0.7940685 -0.4824173  -1.556020 1.244182 -0.981055 
> -1.556020 1.244182 -0.981055  0.6452762 1.6994858
> -0.9265627  0.6452762 1.6994858 -0.9265627  2 1 
> 2 1  0.1923451  0.1923451
```    >> 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@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-devel

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

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

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 - 14:20: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.