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.482 6.485 70.923
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
## 95.425 3.701 99.135
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
## 70.411 5.101 75.437
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
## 6.346 0.452 6.800
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
## 4.508 0.650 5.160
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
## 7.832 0.655 8.488
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
## 85.151 4.493 80.283
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
## 74.567 20.184 94.858
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
## 265.007 12.893 271.821
Sanity check:
stopifnot(all.equal(pca1Ds, pca1s))