Re: [Rd] matrix bands

From: Martin Maechler <maechler_at_stat.math.ethz.ch>
Date: Fri, 26 Aug 2011 14:08:19 +0200

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

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

--->

  See Also:

       ‘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] & 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@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 ]

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

list of date sections of archive