Normalization
In this section we run simple_normalize() on the three different
representations (TENxMatrix, HDF5Matrix, and “HDF5Matrix as sparse”)
of the smaller test dataset only (27,998 x 12,500), and we report the
time of each run.
See the Comprehensive timings obtained on various machines section
below in this document for simple_normalize() and simple_pca() timings
obtained on various machines on all our test datasets and using four different
block sizes: 16 Mb, 40 Mb, 100 Mb, and 250 Mb.
TENxMatrix (sparse)
dim(brain1_s)
## [1] 27998 12500
DelayedArray::setAutoBlockSize(250e6) # set block size to 250 Mb
## automatic block size set to 2.5e+08 bytes (was 1e+08)
system.time(norm_brain1_s <- simple_normalize(brain1_s))
## user system elapsed
## 64.627 5.202 69.778
dim(norm_brain1_s)
## [1] 1000 12500
HDF5Matrix (dense)
dim(brain1_D)
## [1] 27998 12500
DelayedArray::setAutoBlockSize(16e6) # set block size to 16 Mb
## automatic block size set to 1.6e+07 bytes (was 2.5e+08)
system.time(norm_brain1_D <- simple_normalize(brain1_D))
## user system elapsed
## 108.562 2.889 111.457
dim(norm_brain1_D)
## [1] 1000 12500
HDF5Matrix as sparse
dim(brain1_Ds)
## [1] 27998 12500
DelayedArray::setAutoBlockSize(250e6) # set block size to 250 Mb
## automatic block size set to 2.5e+08 bytes (was 1.6e+07)
system.time(norm_brain1_Ds <- simple_normalize(brain1_Ds))
## user system elapsed
## 75.464 4.375 79.779
dim(norm_brain1_Ds)
## [1] 1000 12500
On-disk realization of the normalized datasets
Note that the normalized datasets obtained in the previous section
are DelayedMatrix objects that carry delayed operations. These operations
can be displayed with showtree() e.g. for norm_brain1_s:
class(norm_brain1_s)
## [1] "DelayedMatrix"
## attr(,"package")
## [1] "DelayedArray"
showtree(norm_brain1_s)
## 1000x12500 double, sparse: DelayedMatrix object
## └─ 1000x12500 double, sparse: Unary iso op with args
## └─ 1000x12500 double, sparse: Stack of 1 unary iso op(s)
## └─ 1000x12500 double, sparse: Aperm (perm=c(2,1))
## └─ 12500x1000 double, sparse: Unary iso op with args
## └─ 12500x1000 integer, sparse: Aperm (perm=c(2,1))
## └─ 1000x12500 integer, sparse: Subset
## └─ 27998x1306127 integer, sparse: [seed] TENxMatrixSeed object
The other norm_brain1_* objects carry similar operations.
Before we proceed with PCA, we’re going to write the normalized datasets to
new HDF5 files. This introduces an additional cost, but, on the other hand,
it has the benefit of triggering on-disk realization of the object. This
means that all the delayed operations carried by the object will get realized
on-the-fly before the matrix data actually lands on the disk, making the new
object (TENxMatrix or HDF5Matrix) more efficient for PCA or whatever
block-processed operations will come next.
We will use blocks of 100 Mb for all the writing operations.
DelayedArray::setAutoBlockSize(1e8)
## automatic block size set to 1e+08 bytes (was 2.5e+08)
TENxMatrix (sparse)
dim(norm_brain1_s)
## [1] 1000 12500
system.time(norm_brain1_s <- writeTENxMatrix(norm_brain1_s, level=0))
## user system elapsed
## 5.238 0.458 5.697
The new norm_brain1_s object is a pristine TENxMatrix object:
class(norm_brain1_s)
## [1] "TENxMatrix"
## attr(,"package")
## [1] "HDF5Array"
showtree(norm_brain1_s) # "pristine" object (i.e. no more delayed operations)
## 1000x12500 double, sparse: TENxMatrix object
## └─ 1000x12500 double, sparse: [seed] TENxMatrixSeed object
HDF5Matrix (dense)
dim(norm_brain1_D)
## [1] 1000 12500
system.time(norm_brain1_D <- writeHDF5Array(norm_brain1_D, level=0))
## user system elapsed
## 5.182 0.587 5.771
The new norm_brain1_D object is a pristine HDF5Matrix object:
class(norm_brain1_D)
## [1] "HDF5Matrix"
## attr(,"package")
## [1] "HDF5Array"
showtree(norm_brain1_D) # "pristine" object (i.e. no more delayed operations)
## 1000x12500 double: HDF5Matrix object
## └─ 1000x12500 double: [seed] HDF5ArraySeed object
HDF5Matrix as sparse
dim(norm_brain1_Ds)
## [1] 1000 12500
system.time(norm_brain1_Ds <- writeHDF5Array(norm_brain1_Ds, level=0))
## user system elapsed
## 9.224 0.531 9.755
The new norm_brain1_Ds object is a pristine sparse HDF5Matrix object:
class(norm_brain1_Ds)
## [1] "HDF5Matrix"
## attr(,"package")
## [1] "HDF5Array"
showtree(norm_brain1_Ds) # "pristine" object (i.e. no more delayed operations)
## 1000x12500 double, sparse: HDF5Matrix object
## └─ 1000x12500 double, sparse: [seed] HDF5ArraySeed object
PCA
In this section we run simple_pca() on the normalized datasets obtained
in the previous section and report the time of each run.
See the Comprehensive timings obtained on various machines section
below in this document for simple_normalize() and simple_pca() timings
obtained on various machines on all our test datasets and using four different
block sizes: 16 Mb, 40 Mb, 100 Mb, and 250 Mb.
TENxMatrix (sparse)
DelayedArray::setAutoBlockSize(40e6) # set block size to 40 Mb
## automatic block size set to 4e+07 bytes (was 1e+08)
dim(norm_brain1_s)
## [1] 1000 12500
system.time(pca1s <- simple_PCA(norm_brain1_s))
## user system elapsed
## 83.267 4.645 81.493
HDF5Matrix (dense)
DelayedArray::setAutoBlockSize(1e8) # set block size to 100 Mb
## automatic block size set to 1e+08 bytes (was 4e+07)
dim(norm_brain1_D)
## [1] 1000 12500
system.time(pca1D <- simple_PCA(norm_brain1_D))
## user system elapsed
## 99.816 24.733 126.878
Sanity check:
stopifnot(all.equal(pca1D, pca1s))
HDF5Matrix as sparse
DelayedArray::setAutoBlockSize(40e6) # set block size to 40 Mb
## automatic block size set to 4e+07 bytes (was 1e+08)
dim(norm_brain1_Ds)
## [1] 1000 12500
system.time(pca1Ds <- simple_PCA(norm_brain1_Ds))
## user system elapsed
## 235.947 9.283 236.079
Sanity check:
stopifnot(all.equal(pca1Ds, pca1s))