Exactly what I wanted. Thanks much.

      -thomas


On Sat, 21 Feb 2009, Martin Maechler wrote:

"KK" == Ken Knoblauch <ken.knobla...@inserm.fr>
    on Fri, 20 Feb 2009 14:11:55 +0000 (UTC) writes:

   KK> Thomas Lumley <tlumley <at> u.washington.edu> writes:
   >> I want to construct a symmetric band matrix in the Matrix
   >> package from a matrix where the first column contains
   >> data for the main diagonal, the second column has data
   >> for the first subdiagonal/superdiagonal and so on.
   >>
   >> Since the Matrix will be 10^5 x 10^5 or so, with perhaps
   >> 10-20 non-zero elements above the diagonal per row, I
   >> can't do it by constructing a full matrix and then using
   >> the band() function to subset it to a band matrix.
   >>
   >> Any suggestions?
   >>
   >> -thomas
   >>
   >> Thomas Lumley Assoc. Professor, Biostatistics tlumley
   >> <at> u.washington.edu University of Washington, Seattle
   >>
   KK> I'm not expert but just looking at the package, perhaps
   KK> the spMatrix function (constructor) would do it.


   KK> spMatrix(5, 5, c(2:5, 1:5, 1:4), c(1:4, 1:5, 2:5), rnorm(13))

   KK> 5 x 5 sparse Matrix of class "dgTMatrix"

   KK> [1,] -1.7109379 -0.5576860  .            .          .
   KK> [2,] -0.3219899 -0.9294558  0.001323253  .          .
   KK> [3,]  .          0.4099070  0.646068453  0.5388224  .
   KK> [4,]  .          .         -1.007848357 -0.1117796 -0.5555834
   KK> [5,]  .          .          .           -0.6816979 -0.6052127

Yes, very good, Ken!

My solution uses the successor of spMatrix(), called
sparseMatrix(), unless in the symmetric case, where it class
new("dsTMatrix", ...) directly:

bandSparse <- function(n, m = n, k, listDiag, symmetric = FALSE)
{
   ## Purpose: Compute a band-matrix by speciyfying its (sub-)diagonal(s)
   ## ----------------------------------------------------------------------
   ## Arguments: (n,m) : Matrix dimension
   ##                 k : integer vector of "diagonal numbers",  with identical
   ##                    meaning as in  band(*, k)
   ##          listDiag: list of (sub)diagonals
   ##         symmetric: if TRUE, specify only upper or lower triangle;
   ## ----------------------------------------------------------------------
   ## Author: Martin Maechler, Date: 20 Feb 2009, 22:42
   use.x <- !missing(listDiag)
   len.k <- length(k)
   stopifnot(is.list(listDiag),
              k == as.integer(k), n == as.integer(n), m == as.integer(m))
   if(use.x && length(listDiag) != len.k)
        stop(sprintf("'listDiag' must have the same length (%d) as 'k'", len.k))
   k <- as.integer(k)
   n <- as.integer(n)
   m <- as.integer(m)
   stopifnot(n >= 0, m >= 0, -n+1 <= k, k <= m - 1)
   if(symmetric && any(k < 0) && any(k > 0))
        stop("for symmetric band matrix, only specify upper or lower triangle",
             "\n hence, all k must have the same sign")
   dims <- c(n,m)

   k.lengths <- ## This is a bit "ugly"; I got the cases "by inspection"
        if(n >= m) {
            ifelse(k >= m-n,  m - pmax(0,k), n+k)
        } else { ## n < m
            ifelse(k >= -n+1, n + pmin(0,k), m-k)
        }
   i <- j <- integer(sum(k.lengths))
   if(use.x)
        x <- if(len.k > 0) # carefully getting correct type/mode
            rep.int(listDiag[[1]][1], length(i))
   off.i <- 0L
   for(s in seq_len(len.k)) {
        kk <- k[s] ## *is* integer
        l.kk <- k.lengths[s] ## == length of (sub-)diagonal kk
        ii1 <- seq_len(l.kk)
        ind <- ii1 + off.i
        if(kk >= 0) {
            i[ind] <- ii1
            j[ind] <- ii1 + kk
        } else { ## k < 0
            i[ind] <- ii1 - kk
            j[ind] <- ii1
        }
        if(use.x) {
            xx <- listDiag[[s]]
            if(length(xx) < l.kk)
                warning(sprintf("the %d-th (sub)-diagonal (k = %d) is %s",
                                s, kk, "too short; filling with NA's"))
            x[ind] <- xx[ii1]
        }
        off.i <- off.i + l.kk
   }
   if(symmetric) { ## we should have smarter sparseMatrix()
       UpLo <- if(min(k) >= 0) "U" else "L"
        as(if(use.x) {
            if(is.integer(x)) x <- as.double(x)
            new("dsTMatrix", i= i-1L, j = j-1L, x = x, Dim= dims, uplo=UpLo) } 
else
           new("nsTMatrix", i= i-1L, j = j-1L, Dim= dims, uplo=UpLo),
           "CsparseMatrix")
   }
   else { ## general, not symmetric
        if(use.x)
            sparseMatrix(i=i, j=j, x=x, dims=dims)
        else
            sparseMatrix(i=i, j=j, dims=dims)
   }
}

if(FALSE) { ## Examples

   dList <- list(1:30, 10*(1:20), 100*(1:20))

   s1 <- bandSparse(13, k = -c(0:2, 6), lis = c(dList, dList[2]), symm=TRUE)
   s2 <- bandSparse(13, k =  c(0:2, 6), lis = c(dList, dList[2]), symm=TRUE)
   stopifnot(identical(s1, t(s2)), is(s1,"dsCMatrix"))

   n <- 1e6 ## is already biggish; MM sees top memory up to 1.1 GB;
            ## and it (the bandSparse() call) takes 1 or 2 minutes (slow NB)
   bk <- c(0:5, 7,11)
   bLis <- as.data.frame(matrix(1:8, 1e6, 8, byrow=TRUE))
   B <- bandSparse(1e6, k = bk, list = bLis)
   object.size(B) ## 100 mio Bytes
   Bs <- bandSparse(1e6, k = bk, list = bLis, symmetric=TRUE)
   object.size(Bs)## 100 mio Bytes
   gc()
   show(Bs[1:30, 1:15])
}





Thomas Lumley                   Assoc. Professor, Biostatistics
tlum...@u.washington.edu        University of Washington, Seattle

______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Reply via email to