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:
## <1000000 x 1000000> RandomUnifMatrix object of type "double":
## [,1] [,2] [,3] ... [,999999] [,1000000]
## [1,] 0.5445787 0.5013758 0.9593754 . 0.55048973 0.47121212
## [2,] 0.1206542 0.3429983 0.8109840 . 0.22519468 0.14338738
## [3,] 0.3549180 0.4064891 0.5481541 . 0.19236118 0.11929859
## [4,] 0.2673351 0.1761746 0.5083625 . 0.37233067 0.69728049
## [5,] 0.8718528 0.8326842 0.3711016 . 0.08662414 0.89830577
## ... . . . . . .
## [999996,] 0.90603576 0.72201695 0.07483648 . 0.1342257 0.2856341
## [999997,] 0.81573934 0.91197400 0.67199090 . 0.9863887 0.3270062
## [999998,] 0.85010868 0.09198367 0.35933124 . 0.5665313 0.9851832
## [999999,] 0.80506488 0.41376440 0.70219421 . 0.6818469 0.3175612
## [1000000,] 0.97441563 0.95324582 0.75664554 . 0.5042735 0.4325741The 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.
Almost every distribution in stats is available here. To list a few:
## <100 x 50> RandomNormMatrix object of type "double":
## [,1] [,2] [,3] ... [,49] [,50]
## [1,] 0.57437591 1.67080608 -0.53787926 . 0.1646838 -0.9897245
## [2,] 0.10584850 -1.05566982 0.10263087 . 0.5451388 -0.3861510
## [3,] -0.91415884 -0.57786326 0.53600812 . -0.2588984 0.8100572
## [4,] 2.28042921 -0.60543510 0.06041909 . 0.3238388 -0.7579666
## [5,] -0.34474789 0.07346831 -0.34864905 . 0.8305572 -0.6314178
## ... . . . . . .
## [96,] -0.72680845 -0.04132511 1.52039419 . -0.1871196 -0.3499974
## [97,] -0.86614684 -0.56351662 -0.02220191 . 0.4015388 1.9945382
## [98,] 1.89481048 0.50924323 -0.86620603 . -1.6157598 0.5618521
## [99,] -1.56281632 -0.65826974 -1.00543651 . 0.1782101 2.0409314
## [100,] -0.77624862 -0.69562326 2.97296592 . -0.6861892 -1.3428163## <100 x 50> RandomPoisMatrix object of type "double":
## [,1] [,2] [,3] ... [,49] [,50]
## [1,] 10 5 4 . 5 5
## [2,] 4 6 4 . 6 4
## [3,] 4 5 5 . 8 6
## [4,] 3 10 5 . 2 3
## [5,] 6 4 6 . 4 6
## ... . . . . . .
## [96,] 4 6 5 . 8 7
## [97,] 8 9 7 . 1 7
## [98,] 5 5 3 . 5 6
## [99,] 5 3 7 . 6 7
## [100,] 5 3 2 . 7 2## <100 x 50> RandomGammaMatrix object of type "double":
## [,1] [,2] [,3] ... [,49] [,50]
## [1,] 0.1355294 0.2343008 0.1641899 . 0.63559270 0.42821422
## [2,] 0.4469419 0.6154198 0.5858003 . 0.67435912 0.06900228
## [3,] 0.1749913 0.2482699 0.1526097 . 0.69493046 0.54426192
## [4,] 1.4195444 0.3566785 0.4352261 . 0.13279141 0.41597769
## [5,] 0.3136003 0.1432691 0.4226424 . 0.66323925 0.25595516
## ... . . . . . .
## [96,] 0.12987178 0.19360116 0.33723164 . 0.9904717 0.3949928
## [97,] 0.14106275 0.43181119 0.28056473 . 0.2376704 0.2414953
## [98,] 0.77273936 0.13447534 0.22027570 . 0.1657394 0.9850449
## [99,] 0.35527022 0.62109959 0.44513435 . 0.1815611 0.3162159
## [100,] 0.08816819 0.37331893 0.22763814 . 0.2451349 1.3988476## <100 x 50> RandomWeibullMatrix object of type "double":
## [,1] [,2] [,3] ... [,49] [,50]
## [1,] 0.6471215 0.5961889 1.1183776 . 1.0872129 0.5931966
## [2,] 1.1619233 0.9815352 0.9575827 . 1.0390676 1.3080007
## [3,] 0.5961542 0.9142290 1.1167185 . 0.9082715 0.9623295
## [4,] 1.0196618 0.8637012 0.6962438 . 0.5459533 0.9240985
## [5,] 0.1813089 1.0313057 0.9515449 . 0.7621557 0.9751172
## ... . . . . . .
## [96,] 1.0242366 1.0849748 0.7081675 . 1.0595052 1.1615160
## [97,] 0.7477094 0.7676396 1.1236418 . 0.7726990 0.5697889
## [98,] 1.4910909 0.9744276 0.9659669 . 1.0247492 1.2238684
## [99,] 0.4744775 0.7998682 1.1220066 . 0.9514502 1.1244156
## [100,] 1.2502716 0.9682002 0.4750009 . 1.0116810 0.9886629Distributional parameters can either be scalars:
## <100 x 50> RandomNormMatrix object of type "double":
## [,1] [,2] [,3] ... [,49] [,50]
## [1,] 1.77633411 0.26882013 0.99587797 . 0.7136975 0.5967379
## [2,] 2.02807868 1.67746852 0.03493045 . 0.6866012 2.0694251
## [3,] 0.11996393 2.09587326 2.42702008 . 0.7173492 2.3213342
## [4,] -0.82854808 1.05556343 -0.39507086 . 1.7204409 1.8914305
## [5,] 0.27391543 -0.07783539 0.68493848 . 1.7894565 2.0576551
## ... . . . . . .
## [96,] 0.1111025 1.4696861 0.6048629 . 1.45296611 -0.16728257
## [97,] 1.7744261 1.6817285 -0.2823688 . -0.14609373 0.78914931
## [98,] 0.7847522 1.3235311 0.4760542 . 0.20328228 1.10552910
## [99,] -0.2383344 1.6071420 0.6655284 . 0.47511918 0.07580588
## [100,] 2.0011128 1.8275233 -0.3120024 . 0.13127603 1.72557726Or vectors, which are recycled along the length of the array:
## <100 x 50> RandomNormMatrix object of type "double":
## [,1] [,2] [,3] ... [,49] [,50]
## [1,] 0.67068825 0.05229021 1.20272691 . -1.1133348 -0.0962316
## [2,] -1.09449266 0.85571226 2.59110731 . 3.1688850 3.0241547
## [3,] 2.24150331 1.37904411 2.05603003 . 4.7254942 3.1903058
## [4,] 5.52574239 3.87825200 4.43360149 . 3.2002917 5.4460445
## [5,] 4.22573699 5.17437919 6.54013134 . 4.8659081 3.9445034
## ... . . . . . .
## [96,] 97.42638 98.06552 96.28356 . 94.60765 95.53379
## [97,] 97.26434 95.26832 96.93035 . 96.08934 94.96657
## [98,] 98.26621 98.86981 98.81143 . 98.69824 97.97102
## [99,] 98.19509 99.18958 98.68624 . 100.39596 99.38022
## [100,] 98.10809 101.12615 100.31765 . 98.44123 99.67917Or other arrays of the same dimensions, which are used to sample the corresponding parts of the random array:
## <100 x 50> RandomPoisMatrix object of type "double":
## [,1] [,2] [,3] ... [,49] [,50]
## [1,] 2 0 1 . 1 1
## [2,] 2 1 0 . 0 2
## [3,] 3 1 1 . 2 1
## [4,] 1 2 0 . 2 1
## [5,] 4 0 0 . 3 2
## ... . . . . . .
## [96,] 0 0 0 . 1 0
## [97,] 1 0 1 . 1 0
## [98,] 1 0 6 . 0 3
## [99,] 2 2 2 . 5 5
## [100,] 0 2 0 . 1 0For 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,] 5 5 14 . 30 66
## [2,] 42 1 107 . 4 0
## [3,] 59 0 7 . 0 0
## [4,] 0 12 0 . 16 0
## [5,] 29 1 0 . 0 15
## ... . . . . . .
## [19996,] 17 8 74 . 6 12
## [19997,] 159 12 16 . 18 64
## [19998,] 0 4 0 . 0 2
## [19999,] 66 12 0 . 9 12
## [20000,] 7 0 3 . 3 22Each 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.
## <1000 x 500> RandomUnifMatrix object of type "double":
## [,1] [,2] [,3] ... [,499] [,500]
## [1,] 0.680440166 0.883694901 0.556259171 . 0.6428653 0.5355112
## [2,] 0.254891545 0.042005806 0.636162331 . 0.1360661 0.1167859
## [3,] 0.306360161 0.062161709 0.005092454 . 0.4023786 0.7455437
## [4,] 0.557970699 0.350718441 0.938491546 . 0.6834968 0.7700486
## [5,] 0.973388026 0.994772896 0.265654398 . 0.8018585 0.7067080
## ... . . . . . .
## [996,] 0.16995222 0.91547886 0.98937035 . 0.9990970 0.9448759
## [997,] 0.90338664 0.20239602 0.14016361 . 0.9989606 0.8051184
## [998,] 0.78123588 0.23824850 0.99869329 . 0.6578197 0.9473315
## [999,] 0.28654407 0.72677683 0.02594512 . 0.4078759 0.3573185
## [1000,] 0.90073651 0.33013362 0.08393459 . 0.1116284 0.7170677## <1000 x 500> RandomUnifMatrix object of type "double":
## [,1] [,2] [,3] ... [,499] [,500]
## [1,] 0.78042143 0.24797643 0.94144263 . 0.923906075 0.645603517
## [2,] 0.18896119 0.35186949 0.42286269 . 0.725465333 0.906551875
## [3,] 0.98204860 0.78749705 0.11713305 . 0.739652304 0.004846214
## [4,] 0.50163820 0.14143287 0.14520502 . 0.487257476 0.786116496
## [5,] 0.88308795 0.26678325 0.01130176 . 0.835532834 0.724486915
## ... . . . . . .
## [996,] 0.7862784 0.6114389 0.4585025 . 0.41859961 0.06005039
## [997,] 0.3676248 0.4207759 0.4913526 . 0.74138445 0.03816334
## [998,] 0.2295269 0.8975206 0.4563791 . 0.48240336 0.04912318
## [999,] 0.7062255 0.4162533 0.6312912 . 0.66140499 0.34596018
## [1000,] 0.7040103 0.3162627 0.8549438 . 0.52937940 0.35269592Unlike 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:
## <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## <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.64178865Like 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.
## <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## <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.17849232For 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.
## <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## <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## R version 4.4.2 (2024-10-31)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 24.04.1 LTS
##
## Matrix products: default
## BLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
## LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.26.so; LAPACK version 3.12.0
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=C
## [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: Etc/UTC
## tzcode source: system (glibc)
##
## attached base packages:
## [1] stats4 stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] DelayedRandomArray_1.15.0 DelayedArray_0.33.2
## [3] SparseArray_1.7.2 S4Arrays_1.7.1
## [5] IRanges_2.41.1 abind_1.4-8
## [7] S4Vectors_0.45.2 MatrixGenerics_1.19.0
## [9] matrixStats_1.4.1 BiocGenerics_0.53.3
## [11] generics_0.1.3 Matrix_1.7-1
## [13] BiocStyle_2.35.0
##
## loaded via a namespace (and not attached):
## [1] jsonlite_1.8.9 compiler_4.4.2 BiocManager_1.30.25
## [4] crayon_1.5.3 Rcpp_1.0.13-1 jquerylib_0.1.4
## [7] yaml_2.3.10 fastmap_1.2.0 lattice_0.22-6
## [10] R6_2.5.1 XVector_0.47.0 knitr_1.49
## [13] maketools_1.3.1 bslib_0.8.0 rlang_1.1.4
## [16] cachem_1.1.0 xfun_0.49 sass_0.4.9
## [19] sys_3.4.3 cli_3.6.3 zlibbioc_1.52.0
## [22] digest_0.6.37 grid_4.4.2 dqrng_0.4.1
## [25] lifecycle_1.0.4 evaluate_1.0.1 buildtools_1.0.0
## [28] rmarkdown_2.29 tools_4.4.2 htmltools_0.5.8.1