require(knitr) opts_chunk$set(error=FALSE, message=FALSE, warning=FALSE)
The singular value decomposition (SVD) is an important procedure in general data analysis, primarily due to its use in the principal components analysis (PCA).
r Biocpkg("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
r Biocpkg("BiocParallel") framework.
To perform any SVD, we call the
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)
... but we can change the algorithm by specifying a different
BiocSingularParam class to the
For example, we could use an approximate approach with
r CRANpkg("irlba") instead:
set.seed(1000) i.out <- runSVD(dummy, k=10, BSPARAM=IrlbaParam())
... or we could use a randomized SVD via
set.seed(1000) r.out <- runSVD(dummy, k=10, BSPARAM=RandomParam())
By default, the exact algorithm is used when
The other algorithms are approximate but are much faster than the exact SVD for low
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.
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 cr.out <- runSVD(dummy, k=10, BSPARAM=epam)
fold=1 will compute the cross-product for all matrices, even for square matrices where the gains are negligible.
fold=Inf will never compute the cross-product for any matrix.
The default is to use
Centering and scaling of the input matrix can be performed with the
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
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))
Each individual algorithm has its own parameters that may need tuning.
These can be set in the constructors for the relevant
ipam <- IrlbaParam(tol=1e-8, extra.work=10) ipam
... which will be used in the
irlba() function when
runSVD() is called:
set.seed(3000) i.out <- runSVD(dummy, k=10, BSPARAM=ipam)
Internal matrix multiplications can be parallelized by setting the
This can speed up calculations when dealing with very large input matrices, e.g., from the
r Biocpkg("HDF5Array") package.
set.seed(4000) library(BiocParallel) i.out <- runSVD(dummy, k=10, BSPARAM=ipam, BPPARAM=bpparam())
The number of left or right singular vectors to return can be directly controlled with the
This can avoid wasting time in computing these vectors when they are not required.
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)
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.
Column centering is performed by default with
Standardization of each column can be enabled with
Alternatively, numeric vectors can be passed to either argument for subtraction from and scaling of each column.
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.