Generating random sparse matrix

Question

I am dealing with large matrices - on the order of 10^8 columns and 10^3-10^4 rows. Since these matrices are only ones and zeros (over 99% zeros), I think the sparse construction in the Matrix package is appropriate. However, I don't see a way to generate a random matrix like in the example below. Note that non-zero entries are defined by the column probabilities col_prob.

set.seed(1) #For reproducibility
ncols <- 20
nrows <- 10
col_prob <- runif(ncols,0.1,0.2)
rmat <- matrix(rbinom(nrows*ncols,1,col_prob),
       ncol=ncols,byrow=T)

Certainly I can convert rmat into a sparse matrix:

rmat_sparse <- Matrix(rmat, sparse=TRUE)

However, I would like to generate the sparse matrix in one step. I'm not sure that the function Matrix::rsparsematrix can accomplish this.

Answer 1

The following function will generate a sparse matrix of the type you are looking for, by manipulating the values of a blank dgCMatrix object. It basically creates the rbinom rows one at a time and populates the @i and @p values accordingly.

library(Matrix)    
randsparse <- function(nrows, ncols, col_prob) {
  mat <- Matrix(0, nrows, ncols, sparse = TRUE)  #blank matrix for template
  i <- vector(mode = "list", length = ncols)     #each element of i contains the '1' rows
  p <- rep(0, ncols)                             #p will be cumsum no of 1s by column
  for(r in 1:nrows){
    row <- rbinom(ncols, 1, col_prob)            #random row
    p <- p + row                                 #add to column identifier
    if(any(row == 1)){
      for (j in which(row == 1)){
        i[[j]] <- c(i[[j]], r-1)                 #append row identifier
      }
    }
  }
  p <- c(0, cumsum(p))                           #this is the format required
  i <- unlist(i)
  x <- rep(1, length(i))
  mat@i <- as.integer(i)
  mat@p <- as.integer(p)
  mat@x <- x
  return(mat)
}

set.seed(1)
randsparse(10, 20, runif(20, 0.1, 0.2))

10 x 20 sparse Matrix of class "dgCMatrix"

 [1,] 1 . . . . . . . 1 . . . . . 1 . . . . .
 [2,] . . . . . . . . . . . . . . . . . . . .
 [3,] 1 . . . . . . . . . . . . . . 1 1 . . 1
 [4,] . . . . . . . . . . . . . 1 . . . . . .
 [5,] . . . 1 . . . . 1 . 1 . . . . . . . . .
 [6,] 1 . . . . . . . . . . . . . 1 . . . 1 .
 [7,] . . . . . . . . . . . . . . . . . . . .
 [8,] . 1 . . 1 . . . . . . . 1 . . 1 . . . 1
 [9,] . . 1 . . . . . 1 . . . . 1 . . . 1 . .
[10,] . . . . . . . . . . 1 . . 1 . . . 1 1 .