9

我正在尝试采用两个“矩阵”类（稀疏矩阵）矩阵的元素最大值。我已经尝试过这个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
4

3 回答 3

8

正如您发现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.
于 2012-08-06T16:46:54.793 回答
7

尝试这个。它连接矩阵输出，然后在成对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
于 2012-08-06T17:13:21.593 回答
1

修改 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
于 2015-02-24T18:48:33.613 回答