加法和求和不兼容的矩阵

Question

加法和求和不兼容的矩阵

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

然后我可以将此数据表转换回一个矩阵，并修复行和列名称。

问题是：

- 如何推广此过程？我需要一种自动创建列表的方法，例如list(X1 = X1 + X1.1, X3 = X3 + X3.1, X4 = X4 + X4.1, X5 = X5 + X5.1, X7, X8 = X8 + X8.1)，因为我想将此函数应用于其行列名称和维度未知的矩阵。总之，我需要一个行为如上述描述的合并过程。

- 是否有其他实现相同目标且更快且更通用的策略/实现？（希望一些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)
}

基准测试（使用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）。

- leodido

“%ms%” 是从哪里来的？ - GSee

%ms% 是一个假设的运算符，它实现了所描述的行为。 - leodido

你的矩阵有多大，才会影响执行时间差异？它们总是有很多零吗？如果是这样，可能会有更快的替代方法，使用稀疏矩阵，这些方法更类似于@MadScone的解决方案。 - Aaron left Stack Overflow

是的，它们总是有很多个零。对于这个简单的基准测试，我使用了上面发布的矩阵。稀疏矩阵是另一个好的提示 ;) .. - leodido

1

无论大小还是稀疏度都会对选择哪种解决方案产生巨大影响。请参见下面的编辑。 - Aaron left Stack Overflow

3个回答

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

- mnel

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 没有第七列/行。

最后，cast 转换矩阵，使其以 aggregate 的第一列作为行，第二列作为列的结果编写。有效地撤销了之前的融合。[,-1] 取消了 cast 中多余的一列（我认为可能有更好的方法来做到这一点，但我不知道如何）。

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

- Ciarán Tobin

1

哦，很好。我认为你可以让cast函数执行聚合操作，这样aggregate函数就不再需要了。使用reshape2包，我们可以选择获取矩阵而不是数据框，代码如下：acast(rbind(melt(M1), melt(M2)), Var1~Var2, fun.aggregate=sum)。 - Aaron left Stack Overflow

这是一个不错的解决方案：谢谢！但是，我进行了一些基准测试并在问题中报告了结果。 - leodido

很棒的东西@Aaron。由于某些原因，我仍然在使用“reshape”。该升级了。 - Ciarán Tobin

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- Aaron left Stack Overflow · Accepted Answer

我只需使用基本的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的解决方案，因为它易于编码和理解。当这变慢时，我会编写类似于我的第一次尝试的函数。当那变慢时，我会编写类似于我的第二次尝试的函数。