Real data example

For this post, I will use a big.matrix which represents genotypes for 15,283 individuals, corresponding to the number of mutations (0, 1 or 2) at 287,155 different loci. Here, I will use only the first 10,000 loci (columns).

What you need to know about the big.matrix format:

• you can easily and quickly access matrice-like objects stored on disk,
• you can use different types of storage (I use type char to store each element on only 1 byte),
• it is column-major ordered as standard R matrices,
• you can access elements of a big.matrix using X[i, j] in R,
• you can access elements of a big.matrix using X[j][i] in Rcpp,
• you can get a RcppEigen or RcppArmadillo view of a big.matrix (see Appendix).
• for more details, go to the GitHub repo.

Peek at the data:

print(dim(X))

## [1] 15283 10000

print(X[1:10, 1:12])

##       [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12]
##  [1,]    2    0    2    0    2    2    2    1    2     2     2     2
##  [2,]    2    0    1    2    1    1    1    1    2     1     2     2
##  [3,]    2    0    2    2    2    2    1    1    2     1     2     2
##  [4,]    2    2    0    2    0    0    0    2    2     2     0     2
##  [5,]    2    1    2    2    2    2    2    1    2     2     2     2
##  [6,]    2    1    2    1    2    2    1    1    2     2     2     2
##  [7,]    2    0    2    0    2    2    2    0    2     1     2     2
##  [8,]    2    1    1    2    1    1    1    1    2     1     2     2
##  [9,]    2    1    2    2    2    2    2    2    2     2     2     2
## [10,]    2    0    2    1    2    2    2    0    2     1     2     1

What I needed

I needed a fast matrix-vector multiplication between a big.matrix and a vector. Moreover, I could not use any RcppEigen or RcppArmadillo multiplication because I needed some options of efficiently subsetting columns or rows in my matrix (see Appendix).

Writing this multiplication in Rcpp is no more than two loops:

// [[Rcpp::depends(RcppEigen, bigmemory, BH)]]
#include <RcppEigen.h>
#include <bigmemory/MatrixAccessor.hpp>

using namespace Rcpp;

// [[Rcpp::export]]
NumericVector prod1(XPtr<BigMatrix> bMPtr, const NumericVector& x) {
  
  MatrixAccessor<char> macc(*bMPtr);

  int n = bMPtr->nrow();
  int m = bMPtr->ncol();

  NumericVector res(n);
  int i, j;
  
  for (j = 0; j < m; j++) {
    for (i = 0; i < n; i++) {
      res[i] += macc[j][i] * x[j];
    }
  }

  return res;
}

One test:

y <- rnorm(ncol(X))

print(system.time(
  test <- prod1(X@address, y)
))

##    user  system elapsed 
##   0.664   0.004   0.668

What comes next should be transposable to other applications and other types of data.