[Rd] matrix bands

[Martin Maechler]

>>>>> 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 at r-project.org mailing list
    > https://stat.ethz.ch/mailman/listinfo/r-devel