ScaledMatrix
classThe ScaledMatrix provides yet another method of running
scale() on a matrix. In other words, these three operations
are equivalent:
## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] -1.86019233 2.0338753 -0.9085844 -0.5950883 1.21961525 0.5408903
## [2,] -0.77605190 0.8578318 -0.5294698 -1.2389670 -1.81567152 0.5036694
## [3,] 0.52924527 -0.2131484 1.8131660 -1.6654210 0.03559133 -0.5819311
## [4,] 0.01002606 0.4353355 1.3634738 -1.1645106 -0.83616436 2.0397293
## [5,] 0.83574863 1.4992263 0.1203298 -1.0979191 -0.01018285 1.1501193
## [6,] -0.84927624 1.6425179 1.4792750 -0.6131164 -0.35325179 0.2296403
## [,7] [,8] [,9] [,10]
## [1,] -1.0854147 -0.10832150 -0.6835120 0.7610132
## [2,] -0.4226883 1.28483414 -1.5062218 0.1392321
## [3,] -1.2613290 -0.23140165 -0.7502531 1.9715841
## [4,] -0.8171231 -0.04242887 1.9692675 -0.5650254
## [5,] -0.3702430 -2.53909477 0.9189684 -0.6471497
## [6,] -0.7013654 0.21727726 0.1022758 -0.2782957## <6 x 10> DelayedMatrix object of type "double":
## [,1] [,2] [,3] ... [,9] [,10]
## [1,] -1.86019233 2.03387528 -0.90858440 . -0.6835120 0.7610132
## [2,] -0.77605190 0.85783184 -0.52946980 . -1.5062218 0.1392321
## [3,] 0.52924527 -0.21314837 1.81316602 . -0.7502531 1.9715841
## [4,] 0.01002606 0.43533555 1.36347378 . 1.9692675 -0.5650254
## [5,] 0.83574863 1.49922631 0.12032976 . 0.9189684 -0.6471497
## [6,] -0.84927624 1.64251793 1.47927500 . 0.1022758 -0.2782957## <6 x 10> ScaledMatrix object of type "double":
## [,1] [,2] [,3] ... [,9] [,10]
## [1,] -1.86019233 2.03387528 -0.90858440 . -0.6835120 0.7610132
## [2,] -0.77605190 0.85783184 -0.52946980 . -1.5062218 0.1392321
## [3,] 0.52924527 -0.21314837 1.81316602 . -0.7502531 1.9715841
## [4,] 0.01002606 0.43533555 1.36347378 . 1.9692675 -0.5650254
## [5,] 0.83574863 1.49922631 0.12032976 . 0.9189684 -0.6471497
## [6,] -0.84927624 1.64251793 1.47927500 . 0.1022758 -0.2782957The 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 X with n columns, a vector of column centers c and a vector of column scaling values s, our scaled matrix can be written as:
Y = (X − c ⋅ 1nT)S
where S = diag(s1−1, ..., sn−1). If we wanted to right-multiply it with another matrix A, we would have:
YA = XSA − c ⋅ 1nTSA
The right-most expression is simply the outer product of c with the column sums of SA. More important is the fact that we can use the matrix multiplication operator for X with 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.017 0.000 0.015# The slower way with block processing.
da <- scale(DelayedArray(mat))
system.time(out2 <- da %*% blob)## user system elapsed
## 12.877 5.723 12.688The 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.
## user system elapsed
## 7.542 12.273 5.035Row and column sums are special cases of matrix multiplication and can be computed quickly:
## user system elapsed
## 0.011 0.002 0.007## user system elapsed
## 10.496 0.528 11.032Subsetting, transposition and renaming of the dimensions are all
supported without loss of the ScaledMatrix
representation:
## <20000 x 5> ScaledMatrix object of type "double":
## [,1] [,2] [,3] [,4] [,5]
## [1,] -0.001433462 -0.002912301 0.009737708 0.001718964 -0.014437939
## [2,] -10.044475920 -0.002912301 0.009737708 0.001718964 -0.014437939
## [3,] -0.001433462 -0.002912301 0.009737708 0.001718964 -0.014437939
## [4,] -0.001433462 -0.002912301 0.009737708 0.001718964 -0.014437939
## [5,] -0.001433462 -0.002912301 0.009737708 0.001718964 -0.014437939
## ... . . . . .
## [19996,] -0.001433462 -0.002912301 0.009737708 0.001718964 -0.014437939
## [19997,] -0.001433462 -0.002912301 0.009737708 0.001718964 -0.014437939
## [19998,] -0.001433462 -0.002912301 0.009737708 0.001718964 -0.014437939
## [19999,] -0.001433462 -0.002912301 0.009737708 0.001718964 -0.014437939
## [20000,] -0.001433462 -0.002912301 0.009737708 0.001718964 -0.014437939## <10000 x 20000> ScaledMatrix object of type "double":
## [,1] [,2] [,3] ... [,19999]
## [1,] -0.001433462 -10.044475920 -0.001433462 . -0.001433462
## [2,] -0.002912301 -0.002912301 -0.002912301 . -0.002912301
## [3,] 0.009737708 0.009737708 0.009737708 . 0.009737708
## [4,] 0.001718964 0.001718964 0.001718964 . 0.001718964
## [5,] -0.014437939 -0.014437939 -0.014437939 . -0.014437939
## ... . . . . .
## [9996,] 0.001377804 0.001377804 0.001377804 . 0.001377804
## [9997,] 0.006063338 0.006063338 0.006063338 . 0.006063338
## [9998,] -0.007030875 -0.007030875 -0.007030875 . -0.007030875
## [9999,] 0.010502583 0.010502583 0.010502583 . 0.010502583
## [10000,] 0.008426274 0.008426274 0.008426274 . 0.008426274
## [,20000]
## [1,] -0.001433462
## [2,] -0.002912301
## [3,] 0.009737708
## [4,] 0.001718964
## [5,] -0.014437939
## ... .
## [9996,] 0.001377804
## [9997,] 0.006063338
## [9998,] -0.007030875
## [9999,] 0.010502583
## [10000,] 0.008426274## <20000 x 10000> ScaledMatrix object of type "double":
## [,1] [,2] [,3] ... [,9999]
## GENE_1 -0.001433462 -0.002912301 0.009737708 . 0.010502583
## GENE_2 -10.044475920 -0.002912301 0.009737708 . 0.010502583
## GENE_3 -0.001433462 -0.002912301 0.009737708 . 0.010502583
## GENE_4 -0.001433462 -0.002912301 0.009737708 . 0.010502583
## GENE_5 -0.001433462 -0.002912301 0.009737708 . 0.010502583
## ... . . . . .
## GENE_19996 -0.001433462 -0.002912301 0.009737708 . 0.010502583
## GENE_19997 -0.001433462 -0.002912301 0.009737708 . 0.010502583
## GENE_19998 -0.001433462 -0.002912301 0.009737708 . 0.010502583
## GENE_19999 -0.001433462 -0.002912301 0.009737708 . 0.010502583
## GENE_20000 -0.001433462 -0.002912301 0.009737708 . 0.010502583
## [,10000]
## GENE_1 0.008426274
## GENE_2 0.008426274
## GENE_3 0.008426274
## GENE_4 0.008426274
## GENE_5 0.008426274
## ... .
## GENE_19996 0.008426274
## GENE_19997 0.008426274
## GENE_19998 0.008426274
## GENE_19999 0.008426274
## GENE_20000 0.008426274Other operations will cause the ScaledMatrix to collapse
to the general DelayedMatrix representation, after which
point block processing will be used.
## <20000 x 10000> DelayedMatrix object of type "double":
## [,1] [,2] [,3] ... [,9999] [,10000]
## GENE_1 0.9985665 0.9970877 1.0097377 . 1.010503 1.008426
## GENE_2 -9.0444759 0.9970877 1.0097377 . 1.010503 1.008426
## GENE_3 0.9985665 0.9970877 1.0097377 . 1.010503 1.008426
## GENE_4 0.9985665 0.9970877 1.0097377 . 1.010503 1.008426
## GENE_5 0.9985665 0.9970877 1.0097377 . 1.010503 1.008426
## ... . . . . . .
## GENE_19996 0.9985665 0.9970877 1.0097377 . 1.010503 1.008426
## GENE_19997 0.9985665 0.9970877 1.0097377 . 1.010503 1.008426
## GENE_19998 0.9985665 0.9970877 1.0097377 . 1.010503 1.008426
## GENE_19999 0.9985665 0.9970877 1.0097377 . 1.010503 1.008426
## GENE_20000 0.9985665 0.9970877 1.0097377 . 1.010503 1.008426For most part, the implementation of the multiplication assumes that
the A matrix and the
matrix product are small compared to X. It is also possible to
multiply two ScaledMatrixes 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.00256 -0.00160 -0.00096 -0.00304 0.00064 0.00352In 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.
## R version 4.4.2 (2024-10-31)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 24.04.1 LTS
##
## Matrix products: default
## BLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
## LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.26.so; LAPACK version 3.12.0
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_US.UTF-8 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: Etc/UTC
## tzcode source: system (glibc)
##
## attached base packages:
## [1] stats4 stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] BiocSingular_1.23.0 ScaledMatrix_1.15.0 DelayedArray_0.33.3
## [4] SparseArray_1.7.2 S4Arrays_1.7.1 abind_1.4-8
## [7] IRanges_2.41.2 S4Vectors_0.45.2 MatrixGenerics_1.19.0
## [10] matrixStats_1.4.1 BiocGenerics_0.53.3 generics_0.1.3
## [13] Matrix_1.7-1 BiocStyle_2.35.0
##
## loaded via a namespace (and not attached):
## [1] jsonlite_1.8.9 compiler_4.4.2
## [3] BiocManager_1.30.25 crayon_1.5.3
## [5] rsvd_1.0.5 Rcpp_1.0.13-1
## [7] DelayedMatrixStats_1.29.0 parallel_4.4.2
## [9] jquerylib_0.1.4 BiocParallel_1.41.0
## [11] yaml_2.3.10 fastmap_1.2.0
## [13] lattice_0.22-6 R6_2.5.1
## [15] XVector_0.47.0 knitr_1.49
## [17] maketools_1.3.1 bslib_0.8.0
## [19] rlang_1.1.4 cachem_1.1.0
## [21] xfun_0.49 sass_0.4.9
## [23] sys_3.4.3 cli_3.6.3
## [25] zlibbioc_1.52.0 digest_0.6.37
## [27] grid_4.4.2 irlba_2.3.5.1
## [29] sparseMatrixStats_1.19.0 lifecycle_1.0.4
## [31] evaluate_1.0.1 codetools_0.2-20
## [33] buildtools_1.0.0 beachmat_2.23.4
## [35] rmarkdown_2.29 tools_4.4.2
## [37] htmltools_0.5.8.1