I’m a fan of package {bigmemory}. It’s by far the most convenient solution I found for analyzing large genomic data in R on my computer. I’ve started using it in early 2016 and have also contributed some features.

At first, package {bigstatsr} was using the big.matrix objects of package {bigmemory}. Yet, at some point, I felt the need to become independent of package {bigmemory}. As package {bigstatsr} was a central tool of all my thesis work, I needed to add whatever feature I wanted whenever I wanted to. Thus, I reimplemented an object very similar to the filebacked big.matrix object, called “FBM” (Filebacked Big Matrix, very original) in package {bigstatsr}. These two formats are so similar that you can easily convert (without copying the data) between the two objects.

In this vignette, I explain the main differences between my package {bigstatsr} and the packages of the bigmemory family.

Formats and types

Format

Package {bigmemory} provides 3 types of big.matrix objects:

  • a “RAM” big.matrix, which is not shared between processes and use directly random access memory,
  • a shared big.matrix, which uses some shared memory (still a mystery for me),
  • a filebacked big.matrix (so, shared between processes), which stores the data on disk and access it via memory-mapping.

I placed a lot of interest for shared matrices (filebacked or not). Yet, I encountered memory limitations with the shared big.matrix (non-filebacked). So, at some point, I was using only filebacked big.matrix objects. So, in {bigstatsr}, you will found only the FBM format, which is very similar to the filebacked big.matrix format. To prove it, let us convert from one to the other (without copying the backingfile).

library(bigmemory)
library(bigstatsr)
## [1] "C:\\Users\\au639593\\AppData\\Local\\Temp\\Rtmpw3hTre\\file30ac2ebd53df.bk"
##       [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
##  [1,]    0    0    0    0    0    0    0    0    0     0
##  [2,]    0    0    0    0    0    0    0    0    0     0
##  [3,]    0    0    0    0    0    0    0    0    0     0
##  [4,]    0    0    0    0    0    0    0    0    0     0
##  [5,]    0    0    0    0    0    0    0    0    0     0
##  [6,]    0    0    0    0    0    0    0    0    0     0
##  [7,]    0    0    0    0    0    0    0    0    0     0
##  [8,]    0    0    0    0    0    0    0    0    0     0
##  [9,]    0    0    0    0    0    0    0    0    0     0
## [10,]    0    0    0    0    0    0    0    0    0     0
## [1] "C:/Users/au639593/AppData/Local/Temp/Rtmpw3hTre/file30ac2ebd53df.bk"
##       [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
##  [1,]    0    0    0    0    0    0    0    0    0     0
##  [2,]    0    0    0    0    0    0    0    0    0     0
##  [3,]    0    0    0    0    0    0    0    0    0     0
##  [4,]    0    0    0    0    0    0    0    0    0     0
##  [5,]    0    0    0    0    0    0    0    0    0     0
##  [6,]    0    0    0    0    0    0    0    0    0     0
##  [7,]    0    0    0    0    0    0    0    0    0     0
##  [8,]    0    0    0    0    0    0    0    0    0     0
##  [9,]    0    0    0    0    0    0    0    0    0     0
## [10,]    0    0    0    0    0    0    0    0    0     0
## [1] 2
##       [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
##  [1,]    2    0    1    0    0    0    0    0    0     0
##  [2,]    0    0    1    0    0    0    0    0    0     0
##  [3,]    0    0    1    0    0    0    0    0    0     0
##  [4,]    0    0    1    0    0    0    0    0    0     0
##  [5,]    0    0    1    0    0    0    0    0    0     0
##  [6,]    0    0    1    0    0    0    0    0    0     0
##  [7,]    0    0    1    0    0    0    0    0    0     0
##  [8,]    0    0    1    0    0    0    0    0    0     0
##  [9,]    0    0    1    0    0    0    0    0    0     0
## [10,]    0    0    1    0    0    0    0    0    0     0

Types

Package {bigmemory} handles many types:

  • unsigned char (1-byte unsigned integer)
  • char (1-byte signed integer)
  • short (2-byte signed integer)
  • integer (4-byte signed integer)
  • float (single precision floating-point numbers)
  • double (double precision floating-point numbers)
  • complex

For now, package {bigstatsr} handles the following types:

  • unsigned char
  • unsigned short
  • integer
  • float
  • double

Additionally, the unsigned char type is used in the FBM.code256 format, which instead of accessing integer values ranging from 0 to 255, it uses some code to access 256 arbitrary different values. I use this format in my other R package {bigsnpr} to store genotype dosages.

Class

A big.matrix is basically an S4 class object that stores a pointer to a C++ object (an external pointer). When you restart your R session, this pointer becomes Nil and it may make your R session crash. You need to use a different object, a big.matrix.descriptor (using describe()) which stores enough information to make it possible to create this external pointer again (using attach.big.matrix()). Therefore, one has to often switch between descriptors and big.matrix objects.

For FBMs, I use the nice idea of package bigmemoryExtras. Basically, I use a Reference Class (RC) object with active binding. In this object, I store the external pointer and the information needed to create the pointer to the C++ object. The active binding makes this automatic so that the user never need to use attach.big.matrix() or describe() anymore (and no more session crash!).

What this also means is that you can now serialize a FBM (for example, saving it in an rds file with saveRDS() or using it in a parallel algorithm). For instance, with a standard big.matrix object, you’ll need to pass the descriptor object in parallel algorithms:

X <- FBM(10, 10); X[] <- rnorm(length(X))
bm <- X$bm()
library(foreach)
## Warning: package 'foreach' was built under R version 3.6.2
## Error in unserialize(socklist[[n]]): error reading from connection
## Error in serialize(data, node$con): error writing to connection
##  [1] -4.5987802  5.5223872 -2.2085774  1.0856417  0.6811090  4.4952763
##  [7]  2.9756624  3.8597190  0.2475205  4.7507939

In contrast, you can directly use FBMs in parallel algorithms, as the address will be reattached automatically thanks to the active bindings.

foreach(j = 1:10, .combine = 'c') %dopar% {
  sum(X[, j])
}
##  [1] -4.5987802  5.5223872 -2.2085774  1.0856417  0.6811090  4.4952763
##  [7]  2.9756624  3.8597190  0.2475205  4.7507939

C++ accessors

Let us compute the column sums of a big.matrix object in Rcpp.

##  [1] -4.5987802  5.5223872 -2.2085774  1.0856417  0.6811090  4.4952763
##  [7]  2.9756624  3.8597190  0.2475205  4.7507939

Now, let us do it for an FBM.

##  [1] -4.5987802  5.5223872 -2.2085774  1.0856417  0.6811090  4.4952763
##  [7]  2.9756624  3.8597190  0.2475205  4.7507939

So, the main difference is that {bigmemory} uses macc[j][i] whereas FBM objects use the same accessor in C++ as standard Rcpp matrices, macc(i, j). So, it is easier to adapt existing Rcpp algorithms to be used for FBM objects, e.g. using templates. Note that there is also a sub-FBM accessor, so that you can also use the same algorithms on a subset of the FBM object. For example:

## [1] -4.598780  5.522387 -2.208577  1.085642  0.681109  4.495276

Apply an R function

##    user  system elapsed 
##    0.02    0.01    0.03
##    user  system elapsed 
##    1.73    0.03    1.77
all.equal(test1, true)
## [1] TRUE

The {biganalytics} strategy is to make a loop, which is slow because there are a lot of elements to loop through.

Package {bigstatsr} uses a trade-off between accessing all the matrix at once and accessing only one column/row at each iteration. You can access blocks of the big matrix and apply efficient vectorized R functions to each block, and then combine the results.

##    user  system elapsed 
##    0.03    0.05    0.08
all.equal(test2, true)
## [1] TRUE
##    user  system elapsed 
##    0.05    0.00    0.05
all.equal(sqrt(test3), true)
## [1] TRUE

Matrix operations

##    user  system elapsed 
##    0.21    0.01    0.23
##    user  system elapsed 
##    0.19    0.00    0.18
identical(test[], true)
## [1] TRUE
##    user  system elapsed 
##    0.27    0.06    0.33
identical(test2, true)
## [1] TRUE
##    user  system elapsed 
##    0.21    0.00    0.21
identical(test3[], true)
## [1] TRUE
##    user  system elapsed 
##    0.20    0.02    0.22
identical(test4[], true)
## [1] TRUE

Making functions (not operators) makes it possible to use subsetting and scaling.