7

我的目标是使用（并保留）行名和列名来“求和”两个不兼容的矩阵（具有不同维度的矩阵）。

我想出了这种方法：将矩阵转换为data.table对象，将它们连接起来，然后对列向量求和。

一个例子：

> M1
  1 3 4 5 7 8
1 0 0 1 0 0 0
3 0 0 0 0 0 0
4 1 0 0 0 0 0
5 0 0 0 0 0 0
7 0 0 0 0 1 0
8 0 0 0 0 0 0
> M2
  1 3 4 5 8
1 0 0 1 0 0
3 0 0 0 0 0
4 1 0 0 0 0
5 0 0 0 0 0
8 0 0 0 0 0
> M1 %ms% M2
  1 3 4 5 7 8
1 0 0 2 0 0 0
3 0 0 0 0 0 0
4 2 0 0 0 0 0
5 0 0 0 0 0 0
7 0 0 0 0 1 0
8 0 0 0 0 0 0

这是我的代码：

M1 <- matrix(c(0,0,1,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0), byrow = TRUE, ncol = 6)
colnames(M1) <- c(1,3,4,5,7,8)
M2 <- matrix(c(0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0), byrow = TRUE, ncol = 5)
colnames(M2) <- c(1,3,4,5,8)
# to data.table objects
DT1 <- data.table(M1, keep.rownames = TRUE, key = "rn")
DT2 <- data.table(M2, keep.rownames = TRUE, key = "rn")
# join and sum of common columns
if (nrow(DT1) > nrow(DT2)) {
    A <- DT2[DT1, roll = TRUE]
    A[, list(X1 = X1 + X1.1, X3 = X3 + X3.1, X4 = X4 + X4.1, X5 = X5 + X5.1, X7, X8 = X8 + X8.1), by = rn]
}

输出：

   rn X1 X3 X4 X5 X7 X8
1:  1  0  0  2  0  0  0
2:  3  0  0  0  0  0  0
3:  4  2  0  0  0  0  0
4:  5  0  0  0  0  0  0
5:  7  0  0  0  0  1  0
6:  8  0  0  0  0  0  0

然后我可以将其转换回data.tableamatrix并修复行名和列名。

问题是：

  • 如何概括这个过程？

    我需要一种自动创建的方法，list(X1 = X1 + X1.1, X3 = X3 + X3.1, X4 = X4 + X4.1, X5 = X5 + X5.1, X7, X8 = X8 + X8.1)因为我想将此函数应用于事先不知道维度（和行/列名称）的矩阵。

    总之，我需要一个行为如所述的合并过程。

  • 是否有其他实现相同目标的策略/实施同时更快、更通用？（希望一些data.table怪物帮助我）

  • 这个过程可以吸收什么样的连接（内部、外部等）？

提前致谢。

ps：我使用的是data.table 1.8.2版

编辑 - 解决方案

@Aaron 解决方案。没有外部库，只有基础 R。它也适用于矩阵列表。

add_matrices_1 <- function(...) {
  a <- list(...)
  cols <- sort(unique(unlist(lapply(a, colnames))))
  rows <- sort(unique(unlist(lapply(a, rownames))))
  out <- array(0, dim = c(length(rows), length(cols)), dimnames = list(rows,cols))
  for (m in a) out[rownames(m), colnames(m)] <- out[rownames(m), colnames(m)] + m
  out
}

@MadScone 解决方案。使用reshape2包。它每次调用仅适用于两个矩阵。

add_matrices_2 <- function(m1, m2) {
  m <- acast(rbind(melt(M1), melt(M2)), Var1~Var2, fun.aggregate = sum)
  mn <- unique(colnames(m1), colnames(m2))
  rownames(m) <- mn
  colnames(m) <- mn
  m
}

@Aaron 解决方案。使用Matrix包。它仅适用于稀疏矩阵，也适用于它们的列表。

add_matrices_3 <- function(...) {
  a <- list(...)
  cols <- sort(unique(unlist(lapply(a, colnames))))
  rows <- sort(unique(unlist(lapply(a, rownames))))
  nrows <- length(rows)
  ncols <- length(cols)
  newms <- lapply(a, function(m) {
    s <- summary(m)
    i <- match(rownames(m), rows)[s$i]
    j <- match(colnames(m), cols)[s$j]
    ilj <- i < j
    sparseMatrix(
      i         = ifelse(ilj, i, j),
      j         = ifelse(ilj, j, i),
      x         = s$x,
      dims      = c(nrows, ncols),
      dimnames  = list(rows, cols),
      symmetric = TRUE
    )
  })
  Reduce(`+`, newms)
}

BENCHMARK（100 次带microbenchmark包运行）

Unit: microseconds
   expr                min         lq    median         uq       max
1 add_matrices_1   196.009   257.5865   282.027   291.2735   549.397
2 add_matrices_2 13737.851 14697.9790 14864.778 16285.7650 25567.448

无需评论基准：@Aaron 解决方案获胜。

细节

有关性能的见解（取决于矩阵的大小和稀疏性），请参阅@Aaron 的编辑（以及稀疏矩阵的解决方案：）add_matrices_3。

4

3 回答 3

6

我只是把名字排好，然后带着基地 R 去镇上。

这是一个简单的函数，它采用未指定数量的矩阵并将它们按行/列名称相加。

add_matrices_1 <- function(...) {
  a <- list(...)
  cols <- sort(unique(unlist(lapply(a, colnames))))
  rows <- sort(unique(unlist(lapply(a, rownames))))
  out <- array(0, dim=c(length(rows), length(cols)), dimnames=list(rows,cols))
  for(M in a) { out[rownames(M), colnames(M)] <- out[rownames(M), colnames(M)] + M }
  out
}

然后它像这样工作：

# giving them rownames and colnames
colnames(M1) <- rownames(M1) <- c(1,3,4,5,7,8)
colnames(M2) <- rownames(M2) <- c(1,3,4,5,8)

add_matrices_1(M1, M2)
#   1 3 4 5 7 8
# 1 0 0 2 0 0 0
# 3 0 0 0 0 0 0
# 4 2 0 0 0 0 0
# 5 0 0 0 0 0 0
# 7 0 0 0 0 1 0
# 8 0 0 0 0 0 0

然而，对于更大的矩阵，它的效果并不好。这是一个制作矩阵的函数，从可能性中选择n列，并用非零值N填充点。k（这假设是对称矩阵。）

makeM <- function(N, n, k) {
  s1 <- sample(N, n)
  M1 <- array(0, dim=c(n,n), dimnames=list(s1, s1))
  r1 <- sample(n,k, replace=TRUE)
  c1 <- sample(n,k, replace=TRUE)
  M1[cbind(c(r1,c1), c(c1,r1))] <- sample(N,k)
  M1
}

然后这是另一个使用稀疏矩阵的版本。

add_matrices_3 <- function(...) {
  a <- list(...)
  cols <- sort(unique(unlist(lapply(a, colnames))))
  rows <- sort(unique(unlist(lapply(a, rownames))))
  nrows <- length(rows)
  ncols <- length(cols)
  newms <- lapply(a, function(m) {
    s <- summary(m)
    i <- match(rownames(m), rows)[s$i]
    j <- match(colnames(m), cols)[s$j]
    ilj <- i<j
    sparseMatrix(i=ifelse(ilj, i, j),
                 j=ifelse(ilj, j, i),
                 x=s$x,
                 dims=c(nrows, ncols),
                 dimnames=list(rows, cols), symmetric=TRUE)
  })
  Reduce(`+`, newms)
}

当矩阵大而稀疏时，这个版本肯定更快。（请注意，我没有安排转换为稀疏对称矩阵的时间，希望如果这是一种合适的格式，您将在整个代码中使用该格式。）

set.seed(50)
M1 <- makeM(10000, 5000, 50)
M2 <- makeM(10000, 5000, 50)
mm2 <- Matrix(M2)
mm1 <- Matrix(M1)
system.time(add_matrices_1(M1, M2))
#   user  system elapsed 
#  2.987   0.841   4.133 
system.time(add_matrices_3(mm1, mm2))
#   user  system elapsed 
#  0.042   0.012   0.504

但是当矩阵很小时，我的第一个解决方案仍然更快。

set.seed(50)
M1 <- makeM(100, 50, 20)
M2 <- makeM(100, 50, 20)
mm2 <- Matrix(M2)
mm1 <- Matrix(M1)
microbenchmark(add_matrices_1(M1, M2), add_matrices_3(mm1, mm2))
# Unit: microseconds
#                       expr      min       lq   median        uq       max
# 1   add_matrices_1(M1, M2)  398.495  406.543  423.825  544.0905  43077.27
# 2 add_matrices_3(mm1, mm2) 5734.623 5937.473 6044.007 6286.6675 509584.24

故事的寓意：大小和稀疏很重要。

此外，正确处理比节省几微秒更重要。几乎总是最好使用简单的功能，除非遇到麻烦，否则不要担心速度。所以在小情况下，我更喜欢 MadScone 的解决方案，因为它易于编码且易于理解。当速度变慢时，我会像第一次尝试一样编写一个函数。当速度变慢时，我会编写一个类似于我第二次尝试的函数。

于 2012-11-26T21:20:19.453 回答
3

这是一个data.table解决方案。神奇的是添加.SD组件（两者都具有相同的名称），然后通过引用分配剩余的列。

# a function to quickly get the non key columns
nonkey <- function(DT){ setdiff(names(DT),key(DT))}
# the columns in DT1 only
notinR <- setdiff(nonkey(DT1), nonkey(DT2))

#calculate; .. means "up one level"
result <- DT2[DT1, .SD + .SD, roll= TRUE][,notinR := unclass(DT1[, ..notinR])]

# re set the column order to the original (DT1) order
setcolorder(result, names(DT1))

# voila!
result

   rn 1 3 4 5 7 8
1:  1 0 0 2 0 0 0
2:  3 0 0 0 0 0 0
3:  4 2 0 0 0 0 0
4:  5 0 0 0 0 0 0
5:  7 0 0 0 0 1 0
6:  8 0 0 0 0 0 0

我不相信这是一个特别稳定的解决方案，因为我不确定它不是侥幸获得答案，因为M1并且M2是彼此的子集

编辑，一个丑陋的方法使用eval

这变得更加困难，因为您有非语法名称（`1`等）

inBoth <- intersect(nonkey(DT1), nonKey(DT2))

 backquote <- function(x){paste0('`', x, '`')}
 bqBoth <- backquote(inBoth)

 charexp <- sprintf('list(%s)',paste(c(paste0( bqBoth,'=',  bqBoth, '+ i.',inBoth), backquote(notinR)), collapse = ','))

result2 <- DT2[DT1,eval(parse(text = charexp)), roll = TRUE]
 setcolorder(result2, names(DT1))

# voila!
result2


   rn 1 3 4 5 7 8
1:  1 0 0 2 0 0 0
2:  3 0 0 0 0 0 0
3:  4 2 0 0 0 0 0
4:  5 0 0 0 0 0 0
5:  7 0 0 0 0 1 0
6:  8 0 0 0 0 0 0
于 2012-11-27T04:40:14.290 回答
1

我想我设法用这条恶心的线做到了：

cast(aggregate(value ~ X1 + X2, rbind(melt(M1), melt(M2)), sum), X1 ~ X2)[,-1]

这利用了reshape包。作为数据帧返回，因此根据需要转换为矩阵。

如果您希望采用示例中建议的格式，请尝试以下操作：

"%ms%" <- function(m1, m2) {
  m <- as.matrix(cast(aggregate(value ~ X1 + X2, rbind(melt(m1), melt(m2)), sum), X1 ~ X2)[,-1])
  mn <- unique(colnames(m1), colnames(m2))
  rownames(m) <- mn
  colnames(m) <- mn
  return (m)
}

然后你可以这样做：

M1 %ms% M2

编辑：

解释

显然应该有一些解释抱歉。

melt(M1)

M1从它的原始形式转换成这样的格式（行、列、值）。例如

    1 3 4 5 7 8
  1 0 0 1 0 0 0
  3 0 0 0 0 0 0
  4 1 0 0 0 0 0
  5 0 0 0 0 0 0
  7 0 0 0 0 1 0
  8 0 0 0 0 0 0

转换为：

  X1 X2 value 
1  1  1     0
2  3  1     0
3  4  1     1

等。将两个矩阵中所有可能的（行、列、值）组合M1并M2列出到一个矩阵中。现在这个：

aggregate(value ~ X1 + X2, rbind(melt(M1), melt(M2)), sum)

对行和列相同的值求和。例如，它将在两个矩阵中求和 (1, 1)。和 (3, 1) 等。它不会做任何不存在M2的事情，例如没有第 7 列/行。

最后cast转换矩阵，使其aggregate第一列的结果为行，第二列的结果为列。有效地消除了之前的融化。[,-1]正在删除一个不必要的列（cast我认为可能有更好的方法，但我不知道如何）。

正如我所说，它是作为数据框返回的as.matrix()，因此如果您愿意，请在结果上使用。

于 2012-11-26T21:25:32.847 回答