BiocSingular 1.0.0
The singular value decomposition (SVD) is an important procedure in general data analysis, primarily due to its use in the principal components analysis (PCA). The BiocSingular package provides methods to perform both exact and approximate SVD with a single wrapper function. Our aim is to provide a standard framework for use in other Bioconductor packages, where users can easily switch between different SVD algorithms. We also support parallelization throughout various parts of the SVD calculations via the BiocParallel framework.
# SVD algorithm choices
To perform any SVD, we call the runSVD() function.
The code below performs an exact SVD to obtain the first k singular vectors and values (in this case, 10):
dummy <- matrix(rnorm(10000), ncol=25)
e.out <- runSVD(dummy, k=10, BSPARAM=ExactParam())
str(e.out)## List of 3
## $ d: num [1:10] 24.4 23.9 23.1 22.9 22.5 ...
## $ u: num [1:400, 1:10] -0.09372 -0.02707 0.00729 -0.00641 0.01066 ...
## $ v: num [1:25, 1:10] 0.065061 0.465998 -0.338473 -0.029838 -0.000806 ...… but we can change the algorithm by specifying a different BiocSingularParam class to the BSPARAM= argument.
For example, we could use an approximate approach with irlba instead:
set.seed(1000)
i.out <- runSVD(dummy, k=10, BSPARAM=IrlbaParam())… or we could use a randomized SVD via rsvd:
set.seed(1000)
r.out <- runSVD(dummy, k=10, BSPARAM=RandomParam())By default, the exact algorithm is used when BSPARAM=NULL.
The other algorithms are approximate but are much faster than the exact SVD for low k.
They also involve randomization, so note the use of a random seed.
For tall or fat matrices, it may be faster to perform the SVD on the cross-product rather than on the matrix itself. The cross-product is usually much smaller than the original matrix and can be decomposed rapidly, with time savings offsettting the cost of the cross-product in the first place. The results of the cross-product decomposition is used to obtain the decomposition results for the original matrix.
The definition of “tall” or “fat” is based on the fold difference in dimensions.
All BiocSingularParam objects have a fold slot that specifies the minimum fold difference between dimensions for a cross-product to be computed.
This provides an automated mechanism for deciding whether or not a cross-product should be used.
epam <- ExactParam(fold=10)
epam## class: ExactParam
## cross-product fold-threshold: 10.00
## deferred centering/scaling: offcr.out <- runSVD(dummy, k=10, BSPARAM=epam)Setting fold=1 will compute the cross-product for all matrices, even for square matrices where the gains are negligible.
Setting fold=Inf will never compute the cross-product for any matrix.
The default is to use fold=5.
Centering and scaling of the input matrix can be performed with the center= and scale= arguments, respectively.
These should be numeric vectors of length equal to the number of rows of the input matrix.
The SVD will then be effectively performed on t(t(dummy)-center)/scale).
cs.out <- runSVD(dummy, k=10, scale=runif(ncol(dummy)),
center=rnorm(ncol(dummy)))Many of the approximate algorithms (and computation of the cross-product) involve matrix multiplications. By default, any centering and scaling is applied before any multiplication. However, it is possible to “defer” the centering and scaling such that the multiplication can make use of features of the underlying matrix representation (e.g., sparsity). This enables faster calculations at the expense of numerical stability.
set.seed(2000)
def.out <- runSVD(dummy, k=10, scale=runif(ncol(dummy)),
center=rnorm(ncol(dummy)), BSPARAM=IrlbaParam(deferred=TRUE))The main practical purpose of the SVD is to perform PCAs.
This is achieved with the runPCA() wrapper function, which also accepts a BSPARAM= argument specifying the type of algorithm to use:
pcs.out <- runPCA(dummy, rank=10, BSPARAM=ExactParam())
str(pcs.out)## List of 3
## $ sdev : num [1:10] 1.22 1.2 1.15 1.15 1.12 ...
## $ rotation: num [1:25, 1:10] 0.065141 0.465864 -0.33853 -0.029862 -0.000804 ...
## ..- attr(*, "dimnames")=List of 2
## .. ..$ : NULL
## .. ..$ : chr [1:10] "PC1" "PC2" "PC3" "PC4" ...
## $ x : num [1:400, 1:10] -2.288 -0.662 0.176 -0.157 0.26 ...
## ..- attr(*, "dimnames")=List of 2
## .. ..$ : NULL
## .. ..$ : chr [1:10] "PC1" "PC2" "PC3" "PC4" ...The output is equivalent to that of prcomp() - a matrix of principal component scores, the rotation vectors, and the standard deviation explained by each component.
head(pcs.out$x)## PC1 PC2 PC3 PC4 PC5 PC6
## [1,] -2.28836543 -0.005088504 -0.1215437 1.14325491 2.4615798 -0.3795913
## [2,] -0.66231190 -0.593552869 0.8889813 -0.08355081 -1.6417589 0.5222694
## [3,] 0.17610681 -0.581214302 -1.3314832 -0.40977421 -1.4605159 -0.4954922
## [4,] -0.15681334 -0.783137887 0.6542274 0.11570220 0.8610293 -0.3942107
## [5,] 0.25961665 0.174974164 1.2359999 -0.74731536 -0.5603019 0.9505845
## [6,] 0.05301977 0.539412924 0.4282860 -0.79711774 -1.3153335 -0.1527365
## PC7 PC8 PC9 PC10
## [1,] -1.3398241 1.2252040 -0.43533736 -1.08645020
## [2,] -0.2877867 -0.3340612 0.09030705 0.26171719
## [3,] -1.4815112 -0.1693165 3.53207643 -0.52140348
## [4,] 1.3050660 0.2166619 -0.20253499 -0.70411648
## [5,] 0.4708046 0.1893536 -1.42183396 -0.00935545
## [6,] -2.6508875 -0.9834861 0.04690857 -0.29286364Column centering is performed by default with center=TRUE.
Standardization of each column can be enabled with scale=TRUE.
Alternatively, numeric vectors can be passed to either argument for subtraction from and scaling of each column.
sessionInfo()## R version 3.6.0 (2019-04-26)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows Server 2012 R2 x64 (build 9600)
##
## Matrix products: default
##
## locale:
## [1] LC_COLLATE=C
## [2] LC_CTYPE=English_United States.1252
## [3] LC_MONETARY=English_United States.1252
## [4] LC_NUMERIC=C
## [5] LC_TIME=English_United States.1252
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] BiocParallel_1.18.0 BiocSingular_1.0.0 knitr_1.22
## [4] BiocStyle_2.12.0
##
## loaded via a namespace (and not attached):
## [1] Rcpp_1.0.1 magrittr_1.5 BiocGenerics_0.30.0
## [4] IRanges_2.18.0 lattice_0.20-38 stringr_1.4.0
## [7] tools_3.6.0 parallel_3.6.0 grid_3.6.0
## [10] rsvd_1.0.0 xfun_0.6 snow_0.4-3
## [13] irlba_2.3.3 htmltools_0.3.6 matrixStats_0.54.0
## [16] yaml_2.2.0 digest_0.6.18 bookdown_0.9
## [19] Matrix_1.2-17 BiocManager_1.30.4 S4Vectors_0.22.0
## [22] evaluate_0.13 rmarkdown_1.12 DelayedArray_0.10.0
## [25] stringi_1.4.3 compiler_3.6.0 stats4_3.6.0