我的目标是使用 (并保留) 行和列名称, 对两个不兼容的矩阵(具有不同维度的矩阵)进行求和。
我已经想到了这种方法: 将矩阵转换为 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
)。
%ms%
是一个假设的运算符,它实现了所描述的行为。 - leodido