我的目标是“求和”两个 不兼容的矩阵 (具有不同维度的矩阵)使用(并保留)行和列名称。
我想出了这种方法:将矩阵转换为 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)
因为我想要 将此函数应用于事先未知维度(和行/列名称)的矩阵 .总之,我需要一个 合并 行为如所描述的过程。
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)
}
基准 (使用
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
)。
最佳答案
我只是把名字排好,然后带着基地 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 的解决方案,因为它易于编码且易于理解。当它变慢时,我会写一个像我第一次尝试一样的函数。当它变慢时,我会写一个函数,就像我第二次尝试一样。
关于r - 连接和求和不兼容的矩阵,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/13571359/