ScaledMatrix
classScaledMatrix 1.14.0
The ScaledMatrix
provides yet another method of running scale()
on a matrix.
In other words, these three operations are equivalent:
mat <- matrix(rnorm(10000), ncol=10)
smat1 <- scale(mat)
head(smat1)
## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] -1.0189944 -0.04021421 0.2410524 0.658419549 2.13301741 0.2325958
## [2,] -0.5656357 -1.21521339 -0.6694777 0.008762209 0.23088765 1.7735139
## [3,] 0.2048440 -0.82088235 2.2136665 -2.408495993 -0.07910219 0.1703811
## [4,] -0.1005533 -0.66432255 0.5229881 1.452529571 -0.01343078 0.9846661
## [5,] -0.4576884 0.86536764 0.1756861 -1.685649304 -0.56254312 -0.1955790
## [6,] -0.6138258 -0.96674005 -1.0665189 0.742141596 -0.69392880 0.1944565
## [,7] [,8] [,9] [,10]
## [1,] 0.02809433 1.09892225 -0.22192860 1.5728416
## [2,] -0.75715462 1.57109129 0.07064569 2.2285042
## [3,] -2.63107861 -0.04690459 -0.43939053 0.4494341
## [4,] 0.68847868 0.83362869 0.10802644 0.3029970
## [5,] -0.25035579 -1.15887145 1.00366338 0.4223062
## [6,] 1.53104668 1.88558681 1.38476936 0.4658005
library(DelayedArray)
smat2 <- scale(DelayedArray(mat))
head(smat2)
## <6 x 10> DelayedMatrix object of type "double":
## [,1] [,2] [,3] ... [,9] [,10]
## [1,] -1.01899436 -0.04021421 0.24105241 . -0.22192860 1.57284155
## [2,] -0.56563566 -1.21521339 -0.66947773 . 0.07064569 2.22850416
## [3,] 0.20484395 -0.82088235 2.21366645 . -0.43939053 0.44943408
## [4,] -0.10055334 -0.66432255 0.52298814 . 0.10802644 0.30299704
## [5,] -0.45768842 0.86536764 0.17568613 . 1.00366338 0.42230622
## [6,] -0.61382576 -0.96674005 -1.06651892 . 1.38476936 0.46580047
library(ScaledMatrix)
smat3 <- ScaledMatrix(mat, center=TRUE, scale=TRUE)
head(smat3)
## <6 x 10> ScaledMatrix object of type "double":
## [,1] [,2] [,3] ... [,9] [,10]
## [1,] -1.01899436 -0.04021421 0.24105241 . -0.22192860 1.57284155
## [2,] -0.56563566 -1.21521339 -0.66947773 . 0.07064569 2.22850416
## [3,] 0.20484395 -0.82088235 2.21366645 . -0.43939053 0.44943408
## [4,] -0.10055334 -0.66432255 0.52298814 . 0.10802644 0.30299704
## [5,] -0.45768842 0.86536764 0.17568613 . 1.00366338 0.42230622
## [6,] -0.61382576 -0.96674005 -1.06651892 . 1.38476936 0.46580047
The biggest difference lies in how they behave in downstream matrix operations.
smat1
is an ordinary matrix, with the scaled and centered values fully realized in memory.
Nothing too unusual here.smat2
is a DelayedMatrix
and undergoes block processing whereby chunks are realized and operated on, one at a time.
This sacrifices speed for greater memory efficiency by avoiding a copy of the entire matrix.
In particular, it preserves the structure of the original mat
, e.g., from a sparse or file-backed representation.smat3
is a ScaledMatrix
that refactors certain operations so that they can be applied to the original mat
without any scaling or centering.
This takes advantage of the original data structure to speed up matrix multiplication and row/column sums,
albeit at the cost of numerical precision.Given an original matrix \(\mathbf{X}\) with \(n\) columns, a vector of column centers \(\mathbf{c}\) and a vector of column scaling values \(\mathbf{s}\), our scaled matrix can be written as:
\[ \mathbf{Y} = (\mathbf{X} - \mathbf{c} \cdot \mathbf{1}_n^T) \mathbf{S} \]
where \(\mathbf{S} = \text{diag}(s_1^{-1}, ..., s_n^{-1})\). If we wanted to right-multiply it with another matrix \(\mathbf{A}\), we would have:
\[ \mathbf{YA} = \mathbf{X}\mathbf{S}\mathbf{A} - \mathbf{c} \cdot \mathbf{1}_n^T \mathbf{S}\mathbf{A} \]
The right-most expression is simply the outer product of \(\mathbf{c}\) with the column sums of \(\mathbf{SA}\). More important is the fact that we can use the matrix multiplication operator for \(\mathbf{X}\) with \(\mathbf{SA}\), as this allows us to use highly efficient algorithms for certain data representations, e.g., sparse matrices.
library(Matrix)
mat <- rsparsematrix(20000, 10000, density=0.01)
smat <- ScaledMatrix(mat, center=TRUE, scale=TRUE)
blob <- matrix(runif(ncol(mat) * 5), ncol=5)
system.time(out <- smat %*% blob)
## user system elapsed
## 0.023 0.002 0.024
# The slower way with block processing.
da <- scale(DelayedArray(mat))
system.time(out2 <- da %*% blob)
## user system elapsed
## 28.837 7.718 40.748
The same logic applies for left-multiplication and cross-products.
This allows us to easily speed up high-level operations involving matrix multiplication by just switching to a ScaledMatrix
,
e.g., in approximate PCA algorithms from the BiocSingular package.
library(BiocSingular)
set.seed(1000)
system.time(pcs <- runSVD(smat, k=10, BSPARAM=IrlbaParam()))
## user system elapsed
## 10.307 0.166 10.474
Row and column sums are special cases of matrix multiplication and can be computed quickly:
system.time(rowSums(smat))
## user system elapsed
## 0.008 0.000 0.008
system.time(rowSums(da))
## user system elapsed
## 17.400 4.093 21.496
Subsetting, transposition and renaming of the dimensions are all supported without loss of the ScaledMatrix
representation:
smat[,1:5]
## <20000 x 5> ScaledMatrix object of type "double":
## [,1] [,2] [,3] [,4] [,5]
## [1,] 0.0011816179 0.0055379141 0.0027777196 -0.0026284908 0.0008401287
## [2,] 0.0011816179 0.0055379141 0.0027777196 -0.0026284908 0.0008401287
## [3,] 0.0011816179 0.0055379141 0.0027777196 -0.0026284908 0.0008401287
## [4,] 0.0011816179 0.0055379141 0.0027777196 -0.0026284908 0.0008401287
## [5,] 0.0011816179 0.0055379141 0.0027777196 -0.0026284908 0.0008401287
## ... . . . . .
## [19996,] 0.0011816179 0.0055379141 0.0027777196 -0.0026284908 0.0008401287
## [19997,] 0.0011816179 0.0055379141 0.0027777196 -0.0026284908 0.0008401287
## [19998,] 0.0011816179 0.0055379141 0.0027777196 -0.0026284908 0.0008401287
## [19999,] 0.0011816179 0.0055379141 0.0027777196 -0.0026284908 0.0008401287
## [20000,] 0.0011816179 0.0055379141 0.0027777196 -0.0026284908 0.0008401287
t(smat)
## <10000 x 20000> ScaledMatrix object of type "double":
## [,1] [,2] [,3] ... [,19999]
## [1,] 0.0011816179 0.0011816179 0.0011816179 . 0.0011816179
## [2,] 0.0055379141 0.0055379141 0.0055379141 . 0.0055379141
## [3,] 0.0027777196 0.0027777196 0.0027777196 . 0.0027777196
## [4,] -0.0026284908 -0.0026284908 -0.0026284908 . -0.0026284908
## [5,] 0.0008401287 0.0008401287 0.0008401287 . 0.0008401287
## ... . . . . .
## [9996,] -0.005793043 -0.005793043 -0.005793043 . -0.005793043
## [9997,] 0.008322339 0.008322339 0.008322339 . 0.008322339
## [9998,] -0.019809196 -0.019809196 -0.019809196 . -0.019809196
## [9999,] 0.006689966 0.006689966 0.006689966 . 0.006689966
## [10000,] -0.005910080 -0.005910080 -0.005910080 . -0.005910080
## [,20000]
## [1,] 0.0011816179
## [2,] 0.0055379141
## [3,] 0.0027777196
## [4,] -0.0026284908
## [5,] 0.0008401287
## ... .
## [9996,] -0.005793043
## [9997,] 0.008322339
## [9998,] -0.019809196
## [9999,] 0.006689966
## [10000,] -0.005910080
rownames(smat) <- paste0("GENE_", 1:20000)
smat
## <20000 x 10000> ScaledMatrix object of type "double":
## [,1] [,2] [,3] ... [,9999] [,10000]
## GENE_1 0.001181618 0.005537914 0.002777720 . 0.006689966 -0.005910080
## GENE_2 0.001181618 0.005537914 0.002777720 . 0.006689966 -0.005910080
## GENE_3 0.001181618 0.005537914 0.002777720 . 0.006689966 -0.005910080
## GENE_4 0.001181618 0.005537914 0.002777720 . 0.006689966 -0.005910080
## GENE_5 0.001181618 0.005537914 0.002777720 . 0.006689966 -0.005910080
## ... . . . . . .
## GENE_19996 0.001181618 0.005537914 0.002777720 . 0.006689966 -0.005910080
## GENE_19997 0.001181618 0.005537914 0.002777720 . 0.006689966 -0.005910080
## GENE_19998 0.001181618 0.005537914 0.002777720 . 0.006689966 -0.005910080
## GENE_19999 0.001181618 0.005537914 0.002777720 . 0.006689966 -0.005910080
## GENE_20000 0.001181618 0.005537914 0.002777720 . 0.006689966 -0.005910080
Other operations will cause the ScaledMatrix
to collapse to the general DelayedMatrix
representation, after which point block processing will be used.
smat + 1
## <20000 x 10000> DelayedMatrix object of type "double":
## [,1] [,2] [,3] ... [,9999] [,10000]
## GENE_1 1.001182 1.005538 1.002778 . 1.0066900 0.9940899
## GENE_2 1.001182 1.005538 1.002778 . 1.0066900 0.9940899
## GENE_3 1.001182 1.005538 1.002778 . 1.0066900 0.9940899
## GENE_4 1.001182 1.005538 1.002778 . 1.0066900 0.9940899
## GENE_5 1.001182 1.005538 1.002778 . 1.0066900 0.9940899
## ... . . . . . .
## GENE_19996 1.001182 1.005538 1.002778 . 1.0066900 0.9940899
## GENE_19997 1.001182 1.005538 1.002778 . 1.0066900 0.9940899
## GENE_19998 1.001182 1.005538 1.002778 . 1.0066900 0.9940899
## GENE_19999 1.001182 1.005538 1.002778 . 1.0066900 0.9940899
## GENE_20000 1.001182 1.005538 1.002778 . 1.0066900 0.9940899
For most part, the implementation of the multiplication assumes that the \(\mathbf{A}\) matrix and the matrix product are small compared to \(\mathbf{X}\).
It is also possible to multiply two ScaledMatrix
es together if the underlying matrices have efficient operators for their product.
However, if this is not the case, the ScaledMatrix
offers little benefit for increased overhead.
It is also worth noting that this speed-up is not entirely free.
The expression above involves subtracting two matrix with potentially large values, which runs the risk of catastrophic cancellation.
The example below demonstrates how ScaledMatrix
is more susceptible to loss of precision than a normal DelayedArray
:
set.seed(1000)
mat <- matrix(rnorm(1000000), ncol=100000)
big.mat <- mat + 1e12
# The 'correct' value, unaffected by numerical precision.
ref <- rowMeans(scale(mat))
head(ref)
## [1] -0.0025584703 -0.0008570664 -0.0019225335 -0.0001039903 0.0024761772
## [6] 0.0032943203
# The value from scale'ing a DelayedArray.
library(DelayedArray)
smat2 <- scale(DelayedArray(big.mat))
head(rowMeans(smat2))
## [1] -0.0025583534 -0.0008571123 -0.0019226040 -0.0001039539 0.0024761618
## [6] 0.0032943783
# The value from a ScaledMatrix.
library(ScaledMatrix)
smat3 <- ScaledMatrix(big.mat, center=TRUE, scale=TRUE)
head(rowMeans(smat3))
## [1] -0.00480 0.00848 0.00544 -0.00976 -0.01056 0.01520
In most practical applications, though, this does not seem to be a major concern, especially as most values (e.g., log-normalized expression matrices) lie close to zero anyway.
sessionInfo()
## R version 4.4.1 (2024-06-14)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 24.04.1 LTS
##
## Matrix products: default
## BLAS: /home/biocbuild/bbs-3.20-bioc/R/lib/libRblas.so
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.12.0
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_GB 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: America/New_York
## tzcode source: system (glibc)
##
## attached base packages:
## [1] stats4 stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] BiocSingular_1.22.0 ScaledMatrix_1.14.0 DelayedArray_0.32.0
## [4] SparseArray_1.6.0 S4Arrays_1.6.0 abind_1.4-8
## [7] IRanges_2.40.0 S4Vectors_0.44.0 MatrixGenerics_1.18.0
## [10] matrixStats_1.4.1 BiocGenerics_0.52.0 Matrix_1.7-1
## [13] BiocStyle_2.34.0
##
## loaded via a namespace (and not attached):
## [1] jsonlite_1.8.9 compiler_4.4.1
## [3] BiocManager_1.30.25 crayon_1.5.3
## [5] rsvd_1.0.5 Rcpp_1.0.13
## [7] DelayedMatrixStats_1.28.0 parallel_4.4.1
## [9] jquerylib_0.1.4 BiocParallel_1.40.0
## [11] yaml_2.3.10 fastmap_1.2.0
## [13] lattice_0.22-6 R6_2.5.1
## [15] XVector_0.46.0 knitr_1.48
## [17] bookdown_0.41 bslib_0.8.0
## [19] rlang_1.1.4 cachem_1.1.0
## [21] xfun_0.48 sass_0.4.9
## [23] cli_3.6.3 zlibbioc_1.52.0
## [25] digest_0.6.37 grid_4.4.1
## [27] irlba_2.3.5.1 sparseMatrixStats_1.18.0
## [29] lifecycle_1.0.4 evaluate_1.0.1
## [31] codetools_0.2-20 beachmat_2.22.0
## [33] rmarkdown_2.28 tools_4.4.1
## [35] htmltools_0.5.8.1