1 Introduction

The DelayedRandomArray package implements DelayedArray subclasses containing dynamically sampled random values. Specifically, the actual values are never fully held in memory but are generated when the relevant part of the array is accessed. This allows users to create very large arrays of random values that would not otherwise be possible by filling an ordinary matrix.

To install the package, follow the instructions on DelayedRandomArray landing page. Using the package is then as simple as:

library(DelayedRandomArray)
X <- RandomUnifArray(c(1e6, 1e6))
X
## <1000000 x 1000000> RandomUnifMatrix object of type "double":
##                  [,1]       [,2]       [,3] ...  [,999999] [,1000000]
##       [1,] 0.85328791 0.37369530 0.71554186   . 0.65866861 0.80340566
##       [2,] 0.05579900 0.43667441 0.25903451   . 0.36999319 0.99328220
##       [3,] 0.87722304 0.06033563 0.35807049   . 0.36755489 0.55061337
##       [4,] 0.04485846 0.82747437 0.19475783   . 0.80988399 0.96720696
##       [5,] 0.89834769 0.21547044 0.90036982   . 0.09178073 0.78401539
##        ...          .          .          .   .          .          .
##  [999996,] 0.82835042 0.95633107 0.97512839   .  0.6684754  0.2780033
##  [999997,] 0.21742800 0.86578729 0.07288855   .  0.2304744  0.7841290
##  [999998,] 0.45408800 0.09514295 0.54828125   .  0.5607946  0.6946765
##  [999999,] 0.50690367 0.37166847 0.55463962   .  0.1711145  0.4708021
## [1000000,] 0.36965499 0.37787916 0.26686422   .  0.3917955  0.6792303

The resulting array can be used in any pipeline that is compatible with DelayedArray objects. This object occupies only 64 MB in memory, whereas an ordinary matrix would require 8 PB instead.

2 Available distributions

Almost every distribution in stats is available here. To list a few:

RandomNormArray(c(100, 50))
## <100 x 50> RandomNormMatrix object of type "double":
##               [,1]        [,2]        [,3] ...       [,49]       [,50]
##   [1,]  -0.3959262  -0.1432747  -0.7105505   .  0.07387893  0.88515906
##   [2,]   1.2111220  -0.8693838  -1.5750164   .  2.33468300 -0.34168233
##   [3,]  -0.5087220   1.2074165  -0.5695149   . -0.10222970 -0.45476024
##   [4,]   1.6889490  -0.1994910   1.0103382   .  1.19312030 -0.57133766
##   [5,]  -0.8469510   0.2928810  -0.4797142   .  0.06591360  1.25419105
##    ...           .           .           .   .           .           .
##  [96,]  0.17208531  0.43518079  1.53910228   .  0.95817655  0.23899347
##  [97,]  1.06895435 -1.69531339 -1.91121637   .  1.85478098  0.89531145
##  [98,] -1.42676708  0.64426785  2.09879213   .  0.66999100  0.19304073
##  [99,] -0.05924496 -1.20812631 -0.67535752   . -0.84490851 -1.01294536
## [100,] -1.18939838 -0.64264774 -0.81627320   .  0.09319804 -1.98657689
RandomPoisArray(c(100, 50), lambda=5)
## <100 x 50> RandomPoisMatrix object of type "double":
##         [,1]  [,2]  [,3] ... [,49] [,50]
##   [1,]     5     6     5   .     5     2
##   [2,]     3     6     4   .     5     6
##   [3,]     6     4     2   .     6     3
##   [4,]     2     4     1   .     3     5
##   [5,]     2     5     3   .     5     2
##    ...     .     .     .   .     .     .
##  [96,]     4     4     3   .     6     3
##  [97,]     6     3     6   .     5     4
##  [98,]     5     6     5   .     8     5
##  [99,]     4     4     2   .     3     7
## [100,]     5     6     6   .     7     4
RandomGammaArray(c(100, 50), shape=2, rate=5)
## <100 x 50> RandomGammaMatrix object of type "double":
##              [,1]       [,2]       [,3] ...      [,49]      [,50]
##   [1,] 0.56041620 0.08550097 0.08339802   . 0.46876813 0.81976809
##   [2,] 0.20551171 0.30809741 0.11173671   . 0.07413771 0.62884006
##   [3,] 0.25105851 0.37852599 0.50584298   . 0.38250348 0.07156422
##   [4,] 0.52133533 0.38424874 0.18501654   . 0.12979928 0.59682565
##   [5,] 1.06267859 0.72513974 0.13791043   . 0.16378151 0.21418909
##    ...          .          .          .   .          .          .
##  [96,] 0.64384671 0.10249159 0.45598635   .  0.3291628  0.8177910
##  [97,] 0.17001010 0.33447570 0.39793379   .  0.5897948  1.0953518
##  [98,] 0.46702110 0.09881757 0.72320300   .  0.2250111  0.5275157
##  [99,] 0.47462797 0.33678091 0.45847482   .  0.6354723  0.3029273
## [100,] 0.25406533 0.01565724 0.22233539   .  0.5395694  0.5838045
RandomWeibullArray(c(100, 50), shape=5)
## <100 x 50> RandomWeibullMatrix object of type "double":
##             [,1]      [,2]      [,3] ...     [,49]     [,50]
##   [1,] 0.8986092 0.9897565 1.1105820   . 0.8292627 0.9547774
##   [2,] 0.4693174 1.0880235 0.9516790   . 1.2141667 0.7882257
##   [3,] 1.0392131 0.8846160 1.2421776   . 1.1774108 0.8239251
##   [4,] 1.0455797 0.9074672 1.1652238   . 1.5029199 0.7141953
##   [5,] 1.2097067 0.9138644 0.7434033   . 0.7779865 1.1213414
##    ...         .         .         .   .         .         .
##  [96,] 0.7639854 0.9900883 1.1913165   . 0.6529377 1.0908485
##  [97,] 0.7775855 1.0804989 0.8737506   . 0.7761219 0.8811226
##  [98,] 0.6099293 0.6950335 0.7615443   . 0.6160614 0.7050318
##  [99,] 0.9359207 0.6397807 0.9518008   . 0.5616216 0.9704564
## [100,] 1.0275823 1.5377616 1.1121941   . 0.5361636 0.8609717

Distributional parameters can either be scalars:

RandomNormArray(c(100, 50), mean=1)
## <100 x 50> RandomNormMatrix object of type "double":
##               [,1]        [,2]        [,3] ...      [,49]      [,50]
##   [1,]   2.9477003   1.6411427   0.3990734   .  1.3630894 -0.8013460
##   [2,]   0.3085601   1.7936422   0.9344531   .  1.5930408  0.3358776
##   [3,]  -0.1666180   1.0902513   3.8219420   . -1.2248101  1.9613019
##   [4,]  -0.1003908   0.5963145   1.5332615   .  0.5446926  2.0018224
##   [5,]   1.7973047   3.8472497   0.1980103   . -0.7056252  0.3092759
##    ...           .           .           .   .          .          .
##  [96,]  0.90172474  1.03487077  0.07031717   .  1.2388785 -0.3413436
##  [97,] -0.13545642  0.56302412  0.31622758   .  2.4236288  1.9835031
##  [98,]  3.19088233  1.24741552  1.96787901   .  0.9176761  0.1318759
##  [99,]  1.85738495 -0.16398195  2.16548827   .  1.7381401 -0.1457188
## [100,]  0.23595405  1.05399675  0.02967143   .  0.8805542  1.0494469

Or vectors, which are recycled along the length of the array:

RandomNormArray(c(100, 50), mean=1:100)
## <100 x 50> RandomNormMatrix object of type "double":
##             [,1]      [,2]      [,3] ...       [,49]       [,50]
##   [1,] 1.5494869 0.2016057 2.2277247   .  1.45340944 -0.06032189
##   [2,] 1.9768185 0.5691545 1.7998788   .  3.34157587  3.68873793
##   [3,] 3.1395949 1.9790573 1.9044215   .  2.12389053  4.13198945
##   [4,] 4.5274036 2.9318438 4.7491890   .  4.01933615  4.38185793
##   [5,] 4.8232035 4.5310373 6.9160595   .  3.96451882  4.75193388
##    ...         .         .         .   .           .           .
##  [96,]  96.76274  95.89193  95.57699   .    95.42683    96.30757
##  [97,]  95.98062  97.27188  97.18437   .    96.84343    97.26575
##  [98,]  97.24461  98.23586  96.52861   .    97.03885    97.48717
##  [99,]  99.74484 101.36906  97.29713   .    99.05148   100.37998
## [100,]  97.83604  99.28337  99.86997   .   102.91002    99.92274

Or other arrays of the same dimensions, which are used to sample the corresponding parts of the random array:

means <- RandomNormArray(c(100, 50))
RandomPoisArray(c(100, 50), lambda=2^means)
## <100 x 50> RandomPoisMatrix object of type "double":
##         [,1]  [,2]  [,3] ... [,49] [,50]
##   [1,]     1     2     1   .     0     0
##   [2,]     6     4     0   .     0     0
##   [3,]     2     0     1   .     0     1
##   [4,]     1     1     0   .     0     1
##   [5,]     0     1     3   .     2     1
##    ...     .     .     .   .     .     .
##  [96,]     6     0     0   .     1     2
##  [97,]     1     3     5   .     1     3
##  [98,]     0     1     2   .     0     0
##  [99,]     0     0     0   .     0     1
## [100,]     0     1     0   .     4     3

For example, a hypothetical simulation of a million-cell single-cell RNA-seq dataset might look like this:

ngenes <- 20000
log.abundances <- runif(ngenes, -2, 5)

nclusters <- 20 # define 20 clusters and their population means.
cluster.means <- matrix(2^rnorm(ngenes*nclusters, log.abundances, sd=2), ncol=nclusters)

ncells <- 1e6
clusters <- sample(nclusters, ncells, replace=TRUE) # randomly allocate cells
cell.means <- DelayedArray(cluster.means)[,clusters]

dispersions <- 0.05 + 10/cell.means # typical mean variance trend.

y <- RandomNbinomArray(c(ngenes, ncells), mu=cell.means, size=1/dispersions)
y
## <20000 x 1000000> RandomNbinomMatrix object of type "double":
##                [,1]       [,2]       [,3] ...  [,999999] [,1000000]
##     [1,]          2          0          0   .          0         41
##     [2,]          7          0          0   .          0         32
##     [3,]          1          0          0   .          0          0
##     [4,]         36          7          0   .          0         21
##     [5,]         24          6         46   .         35          5
##      ...          .          .          .   .          .          .
## [19996,]          0          1          0   .          0          5
## [19997,]          0          0          0   .          0          0
## [19998,]          0          0          0   .          0          0
## [19999,]          0          0          0   .          0          2
## [20000,]          4          4          0   .         25         29

3 Chunking

Each random DelayedArrays is broken into contiguous rectangular chunks of identical size and shape. Each chunk is assigned a seed at construction time that is used to initialize a random number stream (using the PCG32 generator from the dqrng package). When the user accesses any part of the array, we generate the random numbers in the overlapping chunks and return the desired values. This provides efficient random access to any subarray without the need to use any jump-ahead functionality.

The chunking scheme determines the efficiency of accessing our random DelayedArrays. Chunks that are too large require unnecessary number generation when a subarray is requested, while chunks that are too small would increase memory usage and book-keeping overhead. The “best” choice also depends on the downstream access pattern, if such information is known. For example, in a matrix where each column is a chunk, retrieval of a column would be very efficient while retrieval of a single row would be very slow. The default chunk dimensions are set to the square root of the array dimensions (or 100, whichever is larger), providing a reasonable compromise between all of these considerations. This can also be manually specified with the chunkdim= argument.

# Row-wise chunks:
RandomUnifArray(c(1000, 500), chunkdim=c(1, 500))
## <1000 x 500> RandomUnifMatrix object of type "double":
##                [,1]        [,2]        [,3] ...     [,499]     [,500]
##    [1,]  0.37226634  0.39919838  0.28537464   . 0.52218746 0.78883913
##    [2,]  0.47511242  0.21105794  0.06129817   . 0.59366773 0.27559308
##    [3,]  0.22525454  0.23468316  0.11345375   . 0.36871533 0.71581683
##    [4,]  0.47630730  0.59396691  0.65288421   . 0.03670637 0.68712554
##    [5,]  0.02879088  0.54073843  0.20558288   . 0.25747474 0.67049646
##     ...           .           .           .   .          .          .
##  [996,] 0.932079839 0.768391092 0.710301048   . 0.59157251 0.09392468
##  [997,] 0.806321804 0.461381833 0.423722171   . 0.28614326 0.17036559
##  [998,] 0.106257328 0.117255454 0.119953389   . 0.47170863 0.17724409
##  [999,] 0.633041224 0.008374148 0.584569189   . 0.46489931 0.51009300
## [1000,] 0.648421318 0.060406642 0.057043200   . 0.22036006 0.63447768
# Column-wise chunks:
RandomUnifArray(c(1000, 500), chunkdim=c(1000, 1))
## <1000 x 500> RandomUnifMatrix object of type "double":
##               [,1]       [,2]       [,3] ...     [,499]     [,500]
##    [1,]  0.3357943  0.7578555  0.4854952   .  0.5670060  0.7300359
##    [2,]  0.7523783  0.2882169  0.4567683   .  0.4489082  0.1831905
##    [3,]  0.4652075  0.9506465  0.9113736   .  0.8924244  0.7680325
##    [4,]  0.6515398  0.1281461  0.7269326   .  0.4632325  0.5288809
##    [5,]  0.6225496  0.8378159  0.7527625   .  0.2598345  0.4756020
##     ...          .          .          .   .          .          .
##  [996,] 0.03366255 0.37670137 0.14619793   . 0.89467090 0.96757530
##  [997,] 0.70469184 0.36156388 0.63886486   . 0.92991270 0.79307689
##  [998,] 0.49098190 0.77164663 0.60820953   . 0.60973434 0.84366857
##  [999,] 0.08713370 0.68073395 0.47815019   . 0.37399044 0.12821038
## [1000,] 0.64075908 0.01144930 0.30457419   . 0.69384228 0.07908277

Unlike other chunk-based DelayedArrays, the actual values of the random DelayedArray are dependent on the chunk parameters. This is because the sampling is done within each chunk and any alteration to the chunk shape or size will rearrange the stream of random numbers within the array. Thus, even when the seed is set, a different chunkdim will yield different results:

set.seed(199)
RandomUnifArray(c(10, 5), chunkdim=c(1, 5))
## <10 x 5> RandomUnifMatrix object of type "double":
##               [,1]         [,2]         [,3]         [,4]         [,5]
##  [1,] 0.4432813900 0.7119991907 0.7257159729 0.9557158216 0.2104000037
##  [2,] 0.3134128384 0.0670862596 0.3807734565 0.1526811279 0.1971345709
##  [3,] 0.5525409468 0.7807079460 0.9905007351 0.0898367104 0.8541555854
##  [4,] 0.8914577304 0.3764661429 0.6028462166 0.2576166289 0.0145520803
##  [5,] 0.5042253651 0.2618340398 0.2405915416 0.6800763716 0.5603335327
##  [6,] 0.3953960796 0.1314800435 0.9008538304 0.1165704445 0.9035755624
##  [7,] 0.6270407308 0.0001234491 0.3542078000 0.3546805161 0.9760325255
##  [8,] 0.2882098067 0.9133890264 0.8018615395 0.7615402588 0.6458653482
##  [9,] 0.2233627134 0.2957882064 0.7371763836 0.3001928469 0.8730950402
## [10,] 0.6481108814 0.5491673953 0.6821353873 0.4434014931 0.1461723484
set.seed(199)
RandomUnifArray(c(10, 5), chunkdim=c(10, 1))
## <10 x 5> RandomUnifMatrix object of type "double":
##             [,1]       [,2]       [,3]       [,4]       [,5]
##  [1,] 0.44328139 0.31341284 0.55254095 0.89145773 0.50422537
##  [2,] 0.71199919 0.06708626 0.78070795 0.37646614 0.26183404
##  [3,] 0.72571597 0.38077346 0.99050074 0.60284622 0.24059154
##  [4,] 0.95571582 0.15268113 0.08983671 0.25761663 0.68007637
##  [5,] 0.21040000 0.19713457 0.85415559 0.01455208 0.56033353
##  [6,] 0.63126383 0.10118984 0.76830283 0.17583353 0.62183470
##  [7,] 0.87478047 0.26076972 0.58018989 0.92146566 0.90680388
##  [8,] 0.55180537 0.20464839 0.38008429 0.56606812 0.75101891
##  [9,] 0.40323090 0.20083487 0.28559824 0.67769322 0.11943279
## [10,] 0.13404760 0.94959186 0.61432837 0.67675354 0.64178865

4 Further comments

Like any other random process, the seed should be set to achieve reproducible results. We stress that the R-level seed only needs to be set before construction of the random DelayedArray; it is not necessary to set the seed during its use. This is because the class itself will define further seeds (one per chunk) and store these in the object for use in per-chunk sampling.

set.seed(999)
RandomNormArray(c(10, 5))
## <10 x 5> RandomNormMatrix object of type "double":
##              [,1]        [,2]        [,3]        [,4]        [,5]
##  [1,]  0.07299897 -0.62967402  0.48686043 -1.00406271  0.09839211
##  [2,] -1.87549362  0.11653501  0.98496315 -0.85515988 -1.44932474
##  [3,]  0.62150210 -0.24704984 -0.69766033 -1.60089511  0.94752723
##  [4,] -1.68035676  0.20939894 -3.29663662  0.57270346  0.94413669
##  [5,]  0.22933664  1.12333362 -0.77725398  1.59886411  1.07042538
##  [6,]  0.43082191 -0.22740387  1.25381398  0.32842239  0.20448699
##  [7,] -0.71424097  0.27773134  1.53748664  0.01761569  1.15470498
##  [8,]  0.35908991  0.35618136 -0.49375892 -0.38652760  0.95591577
##  [9,]  0.06472627  0.79333115 -0.73803086 -0.56434493  0.22473943
## [10,] -2.41389122  0.76273286  0.47567035 -0.08524682 -1.17849232
set.seed(999)
RandomNormArray(c(10, 5))
## <10 x 5> RandomNormMatrix object of type "double":
##              [,1]        [,2]        [,3]        [,4]        [,5]
##  [1,]  0.07299897 -0.62967402  0.48686043 -1.00406271  0.09839211
##  [2,] -1.87549362  0.11653501  0.98496315 -0.85515988 -1.44932474
##  [3,]  0.62150210 -0.24704984 -0.69766033 -1.60089511  0.94752723
##  [4,] -1.68035676  0.20939894 -3.29663662  0.57270346  0.94413669
##  [5,]  0.22933664  1.12333362 -0.77725398  1.59886411  1.07042538
##  [6,]  0.43082191 -0.22740387  1.25381398  0.32842239  0.20448699
##  [7,] -0.71424097  0.27773134  1.53748664  0.01761569  1.15470498
##  [8,]  0.35908991  0.35618136 -0.49375892 -0.38652760  0.95591577
##  [9,]  0.06472627  0.79333115 -0.73803086 -0.56434493  0.22473943
## [10,] -2.41389122  0.76273286  0.47567035 -0.08524682 -1.17849232

For certain distributions, it is possible to indicate that the array is sparse. This does not change the result or efficiency of the sampling process, but can still be useful as it allows downstream functions to use more efficient sparse algorithms. Of course, this is only relevant if the distributional parameters are such that sparsity is actually observed.

RandomPoisArray(c(1e6, 1e6), lambda=0.5) # dense by default
## <1000000 x 1000000> RandomPoisMatrix object of type "double":
##                  [,1]       [,2]       [,3] ...  [,999999] [,1000000]
##       [1,]          0          0          1   .          1          0
##       [2,]          0          1          1   .          0          1
##       [3,]          0          1          1   .          0          0
##       [4,]          0          1          0   .          1          0
##       [5,]          0          0          0   .          1          0
##        ...          .          .          .   .          .          .
##  [999996,]          1          0          0   .          1          1
##  [999997,]          0          1          0   .          1          1
##  [999998,]          2          0          1   .          0          0
##  [999999,]          1          0          1   .          0          1
## [1000000,]          0          0          0   .          1          0
RandomPoisArray(c(1e6, 1e6), lambda=0.5, sparse=TRUE) # treat as sparse
## <1000000 x 1000000> sparse RandomPoisMatrix object of type "double":
##                  [,1]       [,2]       [,3] ...  [,999999] [,1000000]
##       [1,]          0          1          1   .          1          0
##       [2,]          0          1          1   .          1          0
##       [3,]          0          1          1   .          0          0
##       [4,]          1          2          0   .          0          0
##       [5,]          0          1          0   .          0          2
##        ...          .          .          .   .          .          .
##  [999996,]          0          0          1   .          0          0
##  [999997,]          0          0          0   .          1          0
##  [999998,]          3          0          0   .          1          1
##  [999999,]          1          0          0   .          1          0
## [1000000,]          0          2          1   .          1          0

Session information

sessionInfo()
## R version 4.4.0 beta (2024-04-15 r86425)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 22.04.4 LTS
## 
## Matrix products: default
## BLAS:   /home/biocbuild/bbs-3.19-bioc/R/lib/libRblas.so 
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.10.0
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## time zone: America/New_York
## tzcode source: system (glibc)
## 
## attached base packages:
## [1] stats4    stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
##  [1] DelayedRandomArray_1.12.0 DelayedArray_0.30.0      
##  [3] SparseArray_1.4.0         S4Arrays_1.4.0           
##  [5] abind_1.4-5               IRanges_2.38.0           
##  [7] S4Vectors_0.42.0          MatrixGenerics_1.16.0    
##  [9] matrixStats_1.3.0         BiocGenerics_0.50.0      
## [11] Matrix_1.7-0              BiocStyle_2.32.0         
## 
## loaded via a namespace (and not attached):
##  [1] crayon_1.5.2        cli_3.6.2           knitr_1.46         
##  [4] rlang_1.1.3         xfun_0.43           jsonlite_1.8.8     
##  [7] dqrng_0.3.2         htmltools_0.5.8.1   sass_0.4.9         
## [10] rmarkdown_2.26      grid_4.4.0          evaluate_0.23      
## [13] jquerylib_0.1.4     fastmap_1.1.1       yaml_2.3.8         
## [16] lifecycle_1.0.4     bookdown_0.39       BiocManager_1.30.22
## [19] compiler_4.4.0      Rcpp_1.0.12         XVector_0.44.0     
## [22] lattice_0.22-6      digest_0.6.35       R6_2.5.1           
## [25] bslib_0.7.0         tools_4.4.0         zlibbioc_1.50.0    
## [28] cachem_1.0.8