如何计算大型矩阵的行方差？答案

【问题标题】：How can I calculate the row variance for large matrices?如何计算大型矩阵的行方差？
【发布时间】：2021-11-28 14:37:16
【问题描述】：

我在 R 中有一个相当大的稀疏矩阵：

> dim(matrix)
[1] 60675 36807

现在，我想使用以下方法计算此矩阵的行方差：

apply(matrix,1,var)
Error in asMethod(object) : 
    Cholmod error 'problem too large' at file ../Core/cholmod_dense.c, line 102

我想我的矩阵太大了，无法正常工作。我的解决方案是使用多个内核来解决这个问题，如下所示：

mclapply(Matrix::t(matrix), var, mc.cores=16)
Error in asMethod(object) : 
    Cholmod error 'problem too large' at file ../Core/cholmod_dense.c, line 102

但正如您所见，我再次遇到同样的错误。你对如何处理这个问题有什么建议吗？也许对矩阵进行子集化然后计算方差？

【问题讨论】：

你能看看this article吗？我认为这可能会解决您的问题...
你能告诉我们class(matrix) 吗？
如果（1）您需要它比我当前的解决方案快得多并且（2）您给我们更多的上下文/minimal reproducible example（度稀疏/稀疏模式，速度要求...）
我有 RNAseq 表达矩阵数据，我想过滤掉低方差基因（行）。
我不能真正告诉你稀疏度（由于大小问题再次无法计算它），但有大约 42k 行的行和为 0。我还考虑过先删除它们，这会减小这个矩阵的大小。

标签： r matrix

【解决方案1】：

例子：

library(Matrix)
d <- Matrix(0, nrow=60675, ncol=36087)
d[5,5] <- 1

不幸的是，这不起作用：

library(matrixStats)
v <- rowVars(d)

rowVars(d) 中的错误：参数“x”必须是矩阵或向量。

但我们可以通过蛮力做到这一点：

v <- numeric(nrow(d))
for (i in seq(nrow(d))) {if (i %% 250 == 0) cat("."); v[i] <- var(d[i,]) }

这很慢（在老式 MacOS 桌面上大约需要 60 秒）但可以。并行化帮助：

library(parallel)
f <- function(i) var(d[i,])
v <- unlist(mclapply(seq(nrow(d)), f, mc.cores = 4))

（经过 22 秒）。

【讨论】：

没有考虑到这一点，但仅计算具有非零值的行是否可行，沿着 idx <- unique(d@i) 的谎言；然后遍历idx并计算var(d@x[d@i == i])（需要注意i索引）
当然，这样的东西应该可以。是否值得取决于您的稀疏模式：有多少行是空的？
可能还取决于存储模式——如果稀疏矩阵是面向行而不是面向列的，这会更容易一些，但我认为大多数Matrix 类都是面向列的。 ..
TBH 我不知道d[i,] 是否转换为非稀疏矩阵，或者是否可以/应该定义var 以对稀疏结构进行操作——您当然可以通过以下方式做得更好更加努力/多思考。我只是想提出一些可行的方法。
是的 ...（虽然我上面会忽略 var calc 中的零）

【解决方案2】：

Bioconductor 上的sparseMatrixStats 包为Matrix 对象实现了matrixStats API；

# BiocManager::install("sparseMatrixStats")
library(sparseMatrixStats)

mat <- matrix(0, nrow=10, ncol=6)
mat[sample(seq_len(60), 4)] <- 1:4
sparse_mat <- as(mat, "dgCMatrix")
class(sparse_mat)
#> [1] "dgCMatrix"
#> attr(,"package")
#> [1] "Matrix"

rowVars(sparse_mat)
#> [1] 0.0000000 0.1666667 1.5000000 0.0000000 0.0000000 0.0000000 0.0000000
#> [8] 0.0000000 2.8000000 0.0000000

【讨论】：

这在我的机器上需要 0.006 秒（与我的解决方案并行化的 22 秒相比！）