加入和求和不兼容的矩阵答案

【问题标题】：Join and sum not compatible matrices加入和求和不兼容的矩阵
【发布时间】：2012-11-14 07:32:19
【问题描述】：

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

我想出了这种方法：将矩阵转换为 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.table 转换回 matrix 并修复行名和列名。

问题是：

如何概括这个过程？

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

总之，我需要一个 merge 过程，其行为与描述的一样。
还有其他实现相同目标的策略/实现同时更快、更通用吗？（希望有data.table怪物帮帮我）
这个过程可以同化什么样的加入（内部、外部等）？

提前致谢。

p.s.：我使用的是 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（使用 microbenchmark 包运行 100 次）

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）。

【问题讨论】：

%ms% 来自哪里？
%ms% 是一个 ipotetic 运算符，实现了所描述的行为
您的矩阵有多大，执行时间的差异很重要？他们总是有很多零吗？如果是这样，可能会有更快的替代方法使用与@MadScone 的解决方案在本质上更相似的稀疏矩阵。
是的，它们总是有很多零。对于那个简单的基准，我使用了上面发布的矩阵。稀疏矩阵是另一个很好的技巧 ;) ..
大小和稀疏性都会对首选解决方案产生巨大影响。请参阅下面的编辑。

标签： r join matrix merge data.table

【解决方案1】：

这是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

【讨论】：

【解决方案2】：

我只是把名字排好，然后带着基地 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 的解决方案，因为它易于编码且易于理解。当速度变慢时，我会像第一次尝试一样编写一个函数。当速度变慢时，我会像第二次尝试那样编写一个函数。

【讨论】：

出色的解决方案。谢谢！有人建议使用data.table 对象实现相同的行为？

【解决方案3】：

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

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()。

【讨论】：

哦，太好了。我认为您可以让 cast 函数进行聚合，因此不需要 aggregate 函数。使用reshape2 包，以便我们可以选择返回矩阵而不是数据框，它将是acast(rbind(melt(M1), melt(M2)), Var1~Var2, fun.aggregate=sum)。
这是一个很好的解决方案：谢谢！但是，我做了一些基准测试并在问题中报告了结果。
好东西@Aaron。由于某种原因，我仍然坚持使用reshape。是时候升级我的想法了。