R中对稀疏矩阵进行逐元素最大值运算

Question

我正在尝试对两个类别为“矩阵”的矩阵（稀疏矩阵）进行逐元素最大值计算。我尝试使用pmax（...）函数，该函数似乎可以在两个“普通”矩阵上运行, 但是当我传入两个稀疏矩阵时，在R 2.15上会出现以下错误。

library(Matrix)
# Loading required package: lattice
v=Matrix(0,100,100); v[1,1]=1; 
x=v
pmax(v,x)
# Error in pmax(v, x) : (list) object cannot be coerced to type 'logical'
# In addition: Warning message:
# In any(nas) : coercing argument of type 'list' to logical

- user1449378

---

只是确认一下：执行 as.matrix(v) 不是一个选项，因为它是一个大型的稀疏矩阵？ - David Robinson

另外，我们谈论的是多大/多稀疏？ - David Robinson

翻译：太大了，无法转换为标准矩阵——它是100,000乘以100,000。 - user1449378

Answer 1

你发现 pmax 不支持稀疏矩阵。原因是因为 cbind 不支持稀疏矩阵。 Matrix 的作者编写了 cBind，它是 cbind 的等效替代品。如果你在 pmax 函数中更改一行代码，它便可以正确运行：

pmax.sparse=function (..., na.rm = FALSE) 
{
    elts <- list(...)
    if (length(elts) == 0L) 
        stop("no arguments")
    if (all(vapply(elts, function(x) is.atomic(x) && !is.object(x), 
        NA))) {
        mmm <- .Internal(pmax(na.rm, ...))
    }
    else {
        mmm <- elts[[1L]]
        attr(mmm, "dim") <- NULL
        has.na <- FALSE
        for (each in elts[-1L]) {
            attr(each, "dim") <- NULL
            l1 <- length(each)
            l2 <- length(mmm)
            if (l2 < l1) {
                if (l2 && l1%%l2) 
                  warning("an argument will be fractionally recycled")
                mmm <- rep(mmm, length.out = l1)
            }
            else if (l1 && l1 < l2) {
                if (l2%%l1) 
                  warning("an argument will be fractionally recycled")
                each <- rep(each, length.out = l2)
            }
            # nas <- cbind(is.na(mmm), is.na(each))
            nas <- cBind(is.na(mmm), is.na(each)) # Changed row.
            if (has.na || (has.na <- any(nas))) {
                mmm[nas[, 1L]] <- each[nas[, 1L]]
                each[nas[, 2L]] <- mmm[nas[, 2L]]
            }
            change <- mmm < each
            change <- change & !is.na(change)
            mmm[change] <- each[change]
            if (has.na && !na.rm) 
                mmm[nas[, 1L] | nas[, 2L]] <- NA
        }
    }
    mostattributes(mmm) <- attributes(elts[[1L]])
    mmm
}

pmax.sparse(x,v)
# Works fine.

- nograpes

---

当我在大型稀疏矩阵上测试您的方法时，我得到了“Cholmod错误'问题太大'”。即使是较小的问题，我也收到了“<sparse>[ <logic> ] : .M.sub.i.logical()可能效率低下”的警告，并且输出不是预期的“pmax”。有什么想法吗？ - flodel

这个答案实际上是不正确的。正如你所指出的，该函数无法工作。第一个问题是 ! 运算符将稀疏矩阵转换为密集矩阵，通过将所有空单元格转换为 FALSE。我解决了这个问题，但下一个问题是该包中没有实现高效的索引（通过另一个逻辑稀疏矩阵）。看起来这应该是一个容易解决的问题，虽然我还没有尝试过。 - nograpes

Answer 2

尝试这个。它将summary输出的矩阵连接起来，然后在按(i, j)对进行分组后取最大值。它还具有通用性，您可以执行任何类型的元素间操作，只需将max替换为所选函数（或编写一个接受FUN参数的通用函数）。

pmax.sparse <- function(..., na.rm = FALSE) {

   # check that all matrices have conforming sizes
   num.rows <- unique(sapply(list(...), nrow))
   num.cols <- unique(sapply(list(...), ncol))
   stopifnot(length(num.rows) == 1)
   stopifnot(length(num.cols) == 1)

   cat.summary <- do.call(rbind, lapply(list(...), summary))
   out.summary <- aggregate(x ~ i + j, data = cat.summary, max, na.rm)

   sparseMatrix(i = out.summary$i,
                j = out.summary$j,
                x = out.summary$x,
                dims = c(num.rows, num.cols))
}

如果你的矩阵太大，非常稀疏，导致这段代码运行速度过慢，我会考虑使用类似的方法来使用 data.table。以下是一个应用示例：

N <- 1000000
n <- 10000
M1 <- sparseMatrix(i = sample(N,n), j = sample(N,n), x = runif(n), dims = c(N,N))
M2 <- sparseMatrix(i = sample(N,n), j = sample(N,n), x = runif(n), dims = c(N,N))
M3 <- sparseMatrix(i = sample(N,n), j = sample(N,n), x = runif(n), dims = c(N,N))
system.time(p <- pmax.sparse(M1,M2,M3))
#   user  system elapsed 
#   2.58    0.06    2.65

另一个提出的解决方案出现了以下错误：

Error in .class1(object) : 
  Cholmod error 'problem too large' at file ../Core/cholmod_dense.c, line 106

- flodel

Answer 3

对flodel的回答进行修改（无法直接评论回答），以使用data.table包加速处理较大矩阵的计算。

使用原始的flodel版本运行：

> object.size(m1)
# 131053304 bytes
> dim(m1)
# [1] 8031286      39
> object.size(m2)
# 131053304 bytes
> dim(m2)
# [1] 8031286      39
> system.time(pmax.sparse(m1, m2))
# user  system elapsed 
# 326.253  21.805 347.969

修改cat.summary、out.summary和结果矩阵的计算方法为：

cat.summary <- rbindlist(lapply(list(...), summary)) # that's data.table
out.summary <- cat.summary[, list(x = max(x)), by = c("i", "j")]

sparseMatrix(i = out.summary[,i],
             j = out.summary[,j],
             x = out.summary[,x],
             dims = c(num.rows, num.cols))

运行修改后的版本：

> system.time(pmax.sparse(m1, m2))
# user  system elapsed 
# 21.546   0.049  21.589 

- elanmi