1 Overview

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.52573893  0.7225493  1.5262395  0.6702810  0.1757057  0.26399616
## [2,]  0.54502724 -0.2415321  0.6675121  0.5048663 -0.5948648  0.11466458
## [3,]  0.64147056 -1.2232704  2.2003185 -0.3917432  0.6180413  0.76202523
## [4,] -0.89598043  1.4593538 -0.3366019  1.0973632 -1.2881428 -0.13144343
## [5,] -0.06476385  0.2279120 -1.6414587 -0.6777206 -0.2154831 -1.30864938
## [6,]  0.06286008 -1.0902987 -0.8705417 -1.7642926 -0.5471422  0.03480639
##            [,7]        [,8]       [,9]       [,10]
## [1,] -0.4334852 -0.60154755 -0.4152499 -0.66802063
## [2,]  1.4673043  0.87789148 -0.6836570 -0.92698958
## [3,]  0.3733234  0.01965861  0.6769838 -0.26099460
## [4,]  0.7072324  2.59643984  0.7700574  0.32538791
## [5,]  0.8984243  0.36033247 -0.2546239  0.08718273
## [6,]  0.2741784 -1.12148144 -1.6819432 -0.51376923
library(DelayedArray)
smat2 <- scale(DelayedArray(mat))
head(smat2)
## <6 x 10> matrix of class DelayedMatrix and type "double":
##             [,1]        [,2]        [,3] ...        [,9]       [,10]
## [1,]  1.52573893  0.72254933  1.52623946   . -0.41524988 -0.66802063
## [2,]  0.54502724 -0.24153215  0.66751210   . -0.68365701 -0.92698958
## [3,]  0.64147056 -1.22327044  2.20031847   .  0.67698379 -0.26099460
## [4,] -0.89598043  1.45935379 -0.33660192   .  0.77005737  0.32538791
## [5,] -0.06476385  0.22791195 -1.64145870   . -0.25462386  0.08718273
## [6,]  0.06286008 -1.09029869 -0.87054172   . -1.68194325 -0.51376923
library(ScaledMatrix)
smat3 <- ScaledMatrix(mat, center=TRUE, scale=TRUE)
head(smat3)
## <6 x 10> matrix of class ScaledMatrix and type "double":
##             [,1]        [,2]        [,3] ...        [,9]       [,10]
## [1,]  1.52573893  0.72254933  1.52623946   . -0.41524988 -0.66802063
## [2,]  0.54502724 -0.24153215  0.66751210   . -0.68365701 -0.92698958
## [3,]  0.64147056 -1.22327044  2.20031847   .  0.67698379 -0.26099460
## [4,] -0.89598043  1.45935379 -0.33660192   .  0.77005737  0.32538791
## [5,] -0.06476385  0.22791195 -1.64145870   . -0.25462386  0.08718273
## [6,]  0.06286008 -1.09029869 -0.87054172   . -1.68194325 -0.51376923

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.

2 Matrix multiplication

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.022   0.001   0.021
# The slower way with block processing.
da <- scale(DelayedArray(mat))
system.time(out2 <- da %*% blob)
##    user  system elapsed 
##  32.407  10.462  46.573

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 
##  14.791   0.309  15.100

3 Other utilities

Row and column sums are special cases of matrix multiplication and can be computed quickly:

system.time(rowSums(smat))
##    user  system elapsed 
##   0.012   0.000   0.012
system.time(rowSums(da))
##    user  system elapsed 
##  20.335   7.944  28.280

Subsetting, transposition and renaming of the dimensions are all supported without loss of the ScaledMatrix representation:

smat[,1:5]
## <20000 x 5> matrix of class ScaledMatrix and type "double":
##                  [,1]         [,2]         [,3]         [,4]         [,5]
##     [1,] -0.002890835 -0.006727507 -0.004947143  0.006352725 -0.008927862
##     [2,] -0.002890835 -0.006727507 -0.004947143  0.006352725 -0.008927862
##     [3,] -0.002890835 -0.006727507 -0.004947143  0.006352725 -0.008927862
##     [4,] -0.002890835 -0.006727507 -0.004947143  0.006352725 -0.008927862
##     [5,] -0.002890835 -0.006727507 -0.004947143  0.006352725 -0.008927862
##      ...            .            .            .            .            .
## [19996,] -0.002890835 -0.006727507 -0.004947143  0.006352725 -0.008927862
## [19997,] -0.002890835 -0.006727507 -0.004947143  0.006352725 -0.008927862
## [19998,] -0.002890835 -0.006727507 -0.004947143  0.006352725 -0.008927862
## [19999,] -0.002890835 -0.006727507 -0.004947143  0.006352725 -0.008927862
## [20000,] -0.002890835 -0.006727507 -0.004947143  0.006352725 -0.008927862
t(smat)
## <10000 x 20000> matrix of class ScaledMatrix and type "double":
##                  [,1]         [,2]         [,3] ...     [,19999]     [,20000]
##     [1,] -0.002890835 -0.002890835 -0.002890835   . -0.002890835 -0.002890835
##     [2,] -0.006727507 -0.006727507 -0.006727507   . -0.006727507 -0.006727507
##     [3,] -0.004947143 -0.004947143 -0.004947143   . -0.004947143 -0.004947143
##     [4,]  0.006352725  0.006352725  0.006352725   .  0.006352725  0.006352725
##     [5,] -0.008927862 -0.008927862 -0.008927862   . -0.008927862 -0.008927862
##      ...            .            .            .   .            .            .
##  [9996,] -0.002191261 -0.002191261 -0.002191261   . -0.002191261 -0.002191261
##  [9997,] -0.008541684 -0.008541684 -0.008541684   . -0.008541684 -0.008541684
##  [9998,]  0.005359300  0.005359300  0.005359300   .  0.005359300  0.005359300
##  [9999,]  0.007243546  0.007243546  0.007243546   .  0.007243546  0.007243546
## [10000,]  0.007491605  0.007491605  0.007491605   .  0.007491605  0.007491605
rownames(smat) <- paste0("GENE_", 1:20000)
smat
## <20000 x 10000> matrix of class ScaledMatrix and type "double":
##                    [,1]         [,2]         [,3] ...     [,9999]    [,10000]
##     GENE_1 -0.002890835 -0.006727507 -0.004947143   . 0.007243546 0.007491605
##     GENE_2 -0.002890835 -0.006727507 -0.004947143   . 0.007243546 0.007491605
##     GENE_3 -0.002890835 -0.006727507 -0.004947143   . 0.007243546 0.007491605
##     GENE_4 -0.002890835 -0.006727507 -0.004947143   . 0.007243546 0.007491605
##     GENE_5 -0.002890835 -0.006727507 -0.004947143   . 0.007243546 0.007491605
##        ...            .            .            .   .           .           .
## GENE_19996 -0.002890835 -0.006727507 -0.004947143   . 0.007243546 0.007491605
## GENE_19997 -0.002890835 -0.006727507 -0.004947143   . 0.007243546 0.007491605
## GENE_19998 -0.002890835 -0.006727507 -0.004947143   . 0.007243546 0.007491605
## GENE_19999 -0.002890835 -0.006727507 -0.004947143   . 0.007243546 0.007491605
## GENE_20000 -0.002890835 -0.006727507 -0.004947143   . 0.007243546 0.007491605

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> matrix of class DelayedMatrix and type "double":
##                 [,1]      [,2]      [,3] ...  [,9999] [,10000]
##     GENE_1 0.9971092 0.9932725 0.9950529   . 1.007244 1.007492
##     GENE_2 0.9971092 0.9932725 0.9950529   . 1.007244 1.007492
##     GENE_3 0.9971092 0.9932725 0.9950529   . 1.007244 1.007492
##     GENE_4 0.9971092 0.9932725 0.9950529   . 1.007244 1.007492
##     GENE_5 0.9971092 0.9932725 0.9950529   . 1.007244 1.007492
##        ...         .         .         .   .        .        .
## GENE_19996 0.9971092 0.9932725 0.9950529   . 1.007244 1.007492
## GENE_19997 0.9971092 0.9932725 0.9950529   . 1.007244 1.007492
## GENE_19998 0.9971092 0.9932725 0.9950529   . 1.007244 1.007492
## GENE_19999 0.9971092 0.9932725 0.9950529   . 1.007244 1.007492
## GENE_20000 0.9971092 0.9932725 0.9950529   . 1.007244 1.007492

4 Caveats

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 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.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.

Session information

sessionInfo()
## R version 4.1.0 (2021-05-18)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 20.04.2 LTS
## 
## Matrix products: default
## BLAS:   /home/biocbuild/bbs-3.13-bioc/R/lib/libRblas.so
## LAPACK: /home/biocbuild/bbs-3.13-bioc/R/lib/libRlapack.so
## 
## 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       
## 
## attached base packages:
## [1] parallel  stats4    stats     graphics  grDevices utils     datasets 
## [8] methods   base     
## 
## other attached packages:
##  [1] BiocSingular_1.8.0   ScaledMatrix_1.0.0   DelayedArray_0.18.0 
##  [4] IRanges_2.26.0       S4Vectors_0.30.0     MatrixGenerics_1.4.0
##  [7] matrixStats_0.58.0   BiocGenerics_0.38.0  Matrix_1.3-3        
## [10] BiocStyle_2.20.0    
## 
## loaded via a namespace (and not attached):
##  [1] Rcpp_1.0.6                bslib_0.2.5.1            
##  [3] compiler_4.1.0            BiocManager_1.30.15      
##  [5] jquerylib_0.1.4           tools_4.1.0              
##  [7] DelayedMatrixStats_1.14.0 digest_0.6.27            
##  [9] jsonlite_1.7.2            evaluate_0.14            
## [11] lattice_0.20-44           rlang_0.4.11             
## [13] yaml_2.2.1                xfun_0.23                
## [15] stringr_1.4.0             knitr_1.33               
## [17] sass_0.4.0                grid_4.1.0               
## [19] R6_2.5.0                  BiocParallel_1.26.0      
## [21] rmarkdown_2.8             bookdown_0.22            
## [23] irlba_2.3.3               magrittr_2.0.1           
## [25] htmltools_0.5.1.1         sparseMatrixStats_1.4.0  
## [27] rsvd_1.0.5                beachmat_2.8.0           
## [29] stringi_1.6.2