Generating random sparse matrix

Question

My goal is to generate a large sparse matrix with majority (~99%) zeros and ones. Ideally, I would be working with 10,000 rows and 10,000,000 columns. Additionally, each column is generated as a sequence of Bernoulli samples with a column-specific probability. So far, I've implemented 3 ways to generate the data:

Function 1

Creating basic dense matrix of 0/1:

spMat_dense <- function(ncols,nrows,col_probs){
  matrix(rbinom(nrows*ncols,1,col_probs),
         ncol=ncols,byrow=T)
}

Function 2

Using Rcpp:

#include <RcppArmadillo.h>

// [[Rcpp::depends(RcppArmadillo)]]

using namespace std;
using namespace Rcpp;
using namespace arma;

// [[Rcpp::export]]
arma::sp_mat spMat_cpp(const int& ncols, const int& nrows, const NumericVector& col_probs){

  IntegerVector binom_draws = no_init(nrows);
  IntegerVector row_pos;
  IntegerVector col_pos;
  int nz_counter=0;

  //Generate (row,cell)-coordinates of non-zero values
  for(int j=0; j<ncols; ++j){
    binom_draws = rbinom(nrows,1,col_probs[j]);
    for(int i=0; i<nrows; ++i){
      if(binom_draws[i]==1){ 
        row_pos.push_back(i);
        col_pos.push_back(j);
        nz_counter += 1;
      }
    }
  }

  //Create a 2 x N matrix - indicates row/col positions for N non-zero entries
  arma::umat loc_mat(2,nz_counter);

  for(int i=0;i<nz_counter; ++i){
    loc_mat(0,i) = row_pos[i];
    loc_mat(1,i) = col_pos[i];
  }

  IntegerVector x_tmp = rep(1,nz_counter);
  arma::colvec x = Rcpp::as<arma::colvec>(x_tmp);

  //sparse matrix constructor
  arma::sp_mat out(loc_mat,x);
  return out;
}

Function 3

Using dgCMatrix construction in Matrix package:

spMat_dgC <- function(ncols,nrows,col_probs){
  #Credit to Andrew Guster (https://stackoverflow.com/a/56348978/4321711)
  require(Matrix)
  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_probs)            #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)
}

Benchmarking

ncols = 100000
nrows = 1000
col_probs = runif(ncols, 0.001, 0.002)

microbenchmark::microbenchmark(generate_SpMat1(ncols=ncols,nrows=nrows,col_probs=col_probs),
                               generate_SpMat2(ncols=ncols,nrows=nrows,col_probs = col_probs),
                               generate_spMat(ncols=ncols,nrows=nrows,col_probs=col_probs),
                               times=5L)

Unit: seconds
                                                          expr
      spMat_dense(ncols = ncols, nrows = nrows, col_probs = col_probs)
 spMat_cpp(ncols = ncols, nrows = nrows, col_probs = col_probs)
     spMat_dgC(ncols = ncols, nrows = nrows, col_probs = col_probs)
       min        lq      mean   median        uq       max neval
  6.527836  6.673515  7.260482  7.13241  7.813596  8.155053     5
 56.726238 57.038976 57.841693 57.24435 58.325564 59.873333     5
  6.541939  6.599228  6.938952  6.62452  7.402208  7.526867     5

Interestingly, my Rcpp code is not as optimal as I thought it would be. I'm not entirely sure why it's not as efficient as the basic, dense construction. The advantage however in the Rcpp and dgCMatrix construction is that they don't create a dense matrix first. The memory used is much less:

ncols = 100000
nrows = 1000
col_probs = runif(ncols, 0.001, 0.002)

mat1 <- spMat_dense(ncols=ncols,nrows=nrows,col_probs=col_probs)
mat2 <- spMat_cpp(ncols=ncols,nrows=nrows,col_probs = col_probs)
mat3 <- spMat_dgC(ncols=ncols,nrows=nrows,col_probs=col_probs)

object.size(mat1)
object.size(mat2)
object.size(mat3)

> object.size(mat1)
400000216 bytes
> object.size(mat2)
2199728 bytes
> object.size(mat3)
2205920 bytes

Question

What is it about my Rcpp code that makes it slower than the other two? Is it possible to optimize or is the well-written R code with dgCMatrix as good as it gets?

Answer 1

This seems faster in r:

minem2 <- function(ncols, nrows, col_probs){
  r <- lapply(1:ncols, function(x) {
    p <- col_probs[x]
    i <- sample.int(2L, size = nrows, replace = T, prob = c(1 - p, p))
    which(i == 2L)
  })
  rl <- lengths(r)
  nc <- rep(1:ncols, times = rl) # col indexes
  nr <- unlist(r) # row index
  ddims <- c(nrows, ncols)
  sparseMatrix(i = nr, j = nc, dims = ddims)
}

Andrew Guster commented on this approach (link)

Maybe using this logic Rcpp code could be written faster...

Generally, we do not need to generate all values, but just get the indexes where value is 1.