R 矩阵中哪些 rows/columns 与其他哪些重复？
Which rows/columns are duplicates of which others in R matrices?
我有一个矩阵，有很多行和列，性质
x <- matrix(c(1, 1, 3, 3, 55, 55, 1, 3, 3, 1,
1, 1, 3, 3, 55, 55, 1, 3, 9, 1), ncol = 2)
我的问题
在每组重复行中（即每组相同的行），我希望确定第一行索引并将其分配给该组中的所有事件。例如，在两列（第 1、2、7、10 行）中有几个重复行 1。在每一行中，我都想要 第一行索引 ，即 1.
x
# [,1] [,2]
# [1,] 1 1 # first row of 1-1. Assign its row index, 1, to all 1-1 rows
# [2,] 1 1
# [3,] 3 3 # first row of 3-3. Assign its row index, 3, to all 3-3 rows
# [4,] 3 3
# [5,] 55 55 # first row of 55-55. Assign its row index, 5, to all 55-55 rows
# [6,] 55 55
# [7,] 1 1
# [8,] 3 3
# [9,] 3 9 # first (and only) row of 3-9; row index 9
# [10,] 1 1
想要的结果：
1 1 3 3 5 5 1 3 9 1
我的尝试
我想出的最好办法是基于 duplicated 和 for 循环的复杂方法，既不高效也不优雅。我也知道可能 ；那些涉及将行连接成字符串的也非常耗费资源。
# Identify duplicates
duplicate <- duplicated(x, MARGIN = 1)
# Identify first occurrence of each duplicate
firstDup <- duplicated(x, MARGIN = 1, fromLast = TRUE) & !duplicate
indices <- which(firstDup)
# Initialize index for unique rows
index <- seq_len(dim(x)[1])
cf <- duplicate
for (i in indices) {
# Duplicates must occur after first occurrence
cf[seq_len(i)] <- FALSE
dups <- apply(x[cf, , drop = FALSE], 1L, identical, x[i, ])
index[which(cf)[dups]] <- i
}
index
是否有使用 base R 的优雅解决方案？
这个怎么样？
l <- asplit(x, 1L)
match(l, l)
[1] 1 1 3 3 5 5 1 3 9 1
这里，我们使用asplit获取x行的列表l，match获取每行第一次出现的索引.
如果矩阵很大，则以下解决方案可能就足够了：
l <- do.call(paste, data.frame(x))
match(l, l)
[1] 1 1 3 3 5 5 1 3 9 1
TL;DR
对于大小相同但形状不同的整数矩阵（5e+06×2、5e+05×20、5000×2000），包含从 1 到 10 的整数，测试最快的 base 答案是 grouping/match，@alexis_laz 在 中建议。最快的非 base 答案是 data.table::frank/match，尽管 grouping/match 在所有情况下都具有可比性，甚至优于 data.table 答案在 5000 x 2000 的情况下。
请注意，对于范围更大的双精度矩阵或整数矩阵，结果可能会有所不同，具体取决于 data.table 可用的线程数。 [TODO?]
背景
@MikaelJagan 的 asplit/match(<list>, <list>) 似乎是“使用 base R 的优雅解决方案”。但是，?match 警告：
Matching for lists is potentially very slow and best avoided except in simple cases.
鉴于 OP 具有“具有许多行和列的矩阵”，我们想将 asplit/match(<list>, <list>) 答案的性能与另一个 base 的答案进行比较答案：
- @Onyambu 的
paste/match(<chr>, <chr>) ;
- @ThomasIsCoding 的
interaction/match(<int>, <int>) ;
- @alexis_laz的
grouping/match(<int>, <int>).
我们将这些与一些非 base 答案作为基准（我们认识到 OP 仅要求 base）：
- @MikaelJagan 的
Rcpp ;
- @Henrik 的
data.table 回答：
- 将
which = TRUE 和 mult = "first" 传递给 [.data.table 的自连接；
- 两种基于行排名的方法，根据处理关系的方式而有所不同：
frank(ties.method = "average")/match(<dbl>, <dbl>),frank(ties.method = "dense")/match(<int>, <int>).
设置
library(microbenchmark)
library(data.table)
getDTthreads() # 4
f_asplit <- function(x) {
l <- asplit(x, 1L)
match(l, l) }
f_paste <- function(x) {
s <- do.call(paste, as.data.frame(x))
match(s, s) }
f_interaction <- function(x) {
z <- as.integer(interaction(as.data.frame(x)))
match(z, z) }
f_grouping <- function(x) {
g <- do.call(grouping, as.data.frame(x))
o <- order(g, method = "radix")
e <- attr(g, "ends")
z <- rep.int(seq_along(e), c(e[1L], e[-1L] - e[-length(e)]))[o]
match(z, z) }
f_join <- function(x) {
d <- as.data.table(x)
d[d, on = names(d), mult = "first", which = TRUE] }
f_frank_average <- function(x) {
d <- as.data.table(x)
r <- frank(d, ties.method = "average")
match(r, r) }
f_frank_dense <- function(x) {
d <- as.data.table(x)
r <- frank(d, ties.method = "dense")
match(r, r) }
Rcpp::sourceCpp('<copy source code from @MikaelJagan\'s answer here>')
基准测试
行多，列少
我们首先使用 5e+06×2 整数矩阵评估性能：
set.seed(1L)
x <- matrix(sample(10L, size = 1e+07L, replace = TRUE), ncol = 2L)
microbenchmark(
f_asplit(x),
f_paste(x),
f_interaction(x),
f_grouping(x),
f_join(x),
f_frank_average(x),
f_frank_dense(x),
f_rcpp(x),
times = 10L,
check = "identical",
setup = gc(FALSE)
)
Unit: milliseconds
expr min lq mean median uq max neval
f_asplit(x) 17369.93905 18861.91195 19070.21298 19013.0180 19207.29194 22420.71085 10
f_paste(x) 502.63884 507.35077 509.01823 509.2443 511.72301 515.10083 10
f_interaction(x) 234.19311 236.52494 241.80098 238.7392 242.32923 259.75644 10
f_grouping(x) 182.25226 182.89358 187.09642 184.6124 187.10444 208.15532 10
f_join(x) 119.43460 120.86829 123.16607 122.9332 125.07169 128.44722 10
f_frank_average(x) 104.40150 107.53607 111.00268 108.5597 116.80375 121.83675 10
f_frank_dense(x) 86.60926 88.29555 91.42976 90.4716 92.32413 99.30659 10
f_rcpp(x) 459.02304 464.79855 472.43669 468.2492 470.25508 523.06734 10
f_asplit 比 base 替代方案慢两个数量级。 f_grouping 是最快的 base 答案，但 f_frank_dense 快了大约 2 倍（总体上最快）。
更少的行，更多的列
以上结果并未推广到所有整数矩阵输入。例如，f_interaction 与 ncol(x) 的比例非常差：如果 x 的每一列都有 u 个唯一元素，则可能的交互数量为 u^ncol(x)。
出于这个原因，我们进行了第二次基准测试，这次考虑的是行数较少 (5e+05) 列数较多 (20) 的矩阵。
set.seed(1L)
x <- matrix(sample(10L, size = 1e+07L, replace = TRUE), ncol = 20L)
f_interaction 的初始测试导致内存分配错误，因此它被排除在基准测试之外。
system.time(f_interaction(x))
Error: cannot allocate vector of size 7.5 Gb
Timing stopped at: 173.2 6.05 200.4
microbenchmark(
f_asplit(x),
f_paste(x),
## f_interaction(x),
f_grouping(x),
f_join(x),
f_frank_average(x),
f_frank_dense(x),
f_rcpp(x),
times = 10L,
check = "identical",
setup = gc(FALSE)
)
Unit: milliseconds
expr min lq mean median uq max neval
f_asplit(x) 5416.08762 5681.23523 5731.89246 5732.31779 5905.44517 5913.77141 10
f_paste(x) 592.92990 604.15083 629.31101 623.78679 637.81814 724.83871 10
f_grouping(x) 63.89522 64.14134 65.42723 65.11530 66.00557 68.06045 10
f_join(x) 340.73722 342.18096 353.35774 352.08861 359.88480 382.13480 10
f_frank_average(x) 69.90496 70.81840 72.29819 72.04409 73.11977 77.44347 10
f_frank_dense(x) 52.58033 53.33760 54.42029 54.01672 55.63532 56.99664 10
f_rcpp(x) 184096.21999 184816.36584 185774.76817 186218.58335 186696.31674 186781.24972 10
f_grouping 仍然是最快的 base 答案。值得注意的是，它现在比 f_paste 快一个数量级，只比 f_frank_dense.
慢一点点
更少的行，更多的列
我们执行了最终基准测试，排除了上一轮中最慢的答案（f_asplit 和 f_rcpp），现在考虑一个 5000×2000 整数矩阵：
set.seed(1L)
x <- matrix(sample(10L, size = 1e+07L, replace = TRUE), ncol = 2000L)
microbenchmark(
## f_asplit(x),
f_paste(x),
## f_interaction(x),
f_grouping(x),
f_join(x),
f_frank_average(x),
f_frank_dense(x),
## f_rcpp(x),
times = 10L,
check = "identical",
setup = gc(FALSE)
)
Unit: milliseconds
expr min lq mean median uq max neval
f_paste(x) 1067.47994 1075.45148 1083.17391 1080.72997 1089.74027 1102.45249 10
f_grouping(x) 19.24007 19.50026 19.86404 19.79002 20.25302 20.60127 10
f_join(x) 616.66706 621.29854 630.61460 628.16315 636.39097 650.16180 10
f_frank_average(x) 59.82007 61.41706 62.68610 62.99318 64.56520 64.88463 10
f_frank_dense(x) 58.03648 60.59857 63.50526 61.99278 66.03694 71.30638 10
现在 f_grouping 总体上最快，比 f_frank_dense 快大约 3 倍。
如果您使用的是 base R
，我们可以使用 ave
> ave(1:nrow(x), x[, 1], x[, 2], FUN = function(v) v[1])
[1] 1 1 3 3 5 5 1 3 9 1
如果你有多个列，你可以试试
> z <- as.integer(interaction(as.data.frame(m2)))
> ave(seq_along(z), z, FUN = function(x) x[1])
[1] 1 1 3 3 5 5 1 3 9 1
或
> z <- as.integer(interaction(as.data.frame(x)))
> match(z, z)
[1] 1 1 3 3 5 5 1 3 9 1
基准测试
set.seed(1)
v <- sample(1:10, 1e7, replace = TRUE)
m2 <- matrix(v, ncol = 2)
microbenchmark(
ave1 = {
ave(1:nrow(m2), m2[, 1], m2[, 2], FUN = function(v) v[1])
},
ave2 = {
z <- as.integer(interaction(as.data.frame(m2)))
ave(seq_along(z), z, FUN = function(x) x[1])
},
match = {
z <- as.integer(interaction(as.data.frame(m2)))
match(z, z)
},
times = 10L
)
我们会看到
Unit: milliseconds
expr min lq mean median uq max neval
ave1 648.0755 655.9521 715.8848 701.1927 747.4759 885.9838 10
ave2 785.4868 883.2935 913.3867 899.1789 929.6571 1050.9020 10
match 417.1598 447.3718 507.0462 495.8791 551.9436 625.0841 10
不是基础 R，而是另一个比较点：这是一个天真的 Rcpp 实现。最坏的情况是O(n * m^2 / 2)，其中m = nrow(x)和n = ncol(x)，这实际上是相当糟糕的。在行上并行化循环可能会有很大帮助。
它的优点是不复制矩阵x，所以根本不需要太多内存。目前给出的答案调用asplit、as.data.frame或as.data.table，全部复制（分区）x.
Rcpp::sourceCpp(code = '
#include <Rcpp.h>
using namespace Rcpp;
// [[Rcpp::export]]
IntegerVector f_rcpp(IntegerMatrix x)
{
int m = x.nrow();
int n = x.ncol();
IntegerVector res(m);
if (n == 0) {
res.fill(1);
} else {
int i, ki, p, kp, j;
for (i = 0; i < m; ++i) {
for (p = 0; p < i; ++p) {
for (j = 0, ki = i, kp = p; j < n; ++j) {
if (x[kp] != x[ki]) {
break;
}
ki += m;
kp += m;
}
if (j == n) {
res[i] = p + 1;
break;
}
}
if (p == i) {
res[i] = i + 1;
}
}
}
return res;
}
')
f_rcpp(x)
[1] 1 1 3 3 5 5 1 3 9 1
我有一个矩阵，有很多行和列，性质
x <- matrix(c(1, 1, 3, 3, 55, 55, 1, 3, 3, 1,
1, 1, 3, 3, 55, 55, 1, 3, 9, 1), ncol = 2)
我的问题
在每组重复行中（即每组相同的行），我希望确定第一行索引并将其分配给该组中的所有事件。例如，在两列（第 1、2、7、10 行）中有几个重复行 1。在每一行中，我都想要 第一行索引 ，即 1.
x
# [,1] [,2]
# [1,] 1 1 # first row of 1-1. Assign its row index, 1, to all 1-1 rows
# [2,] 1 1
# [3,] 3 3 # first row of 3-3. Assign its row index, 3, to all 3-3 rows
# [4,] 3 3
# [5,] 55 55 # first row of 55-55. Assign its row index, 5, to all 55-55 rows
# [6,] 55 55
# [7,] 1 1
# [8,] 3 3
# [9,] 3 9 # first (and only) row of 3-9; row index 9
# [10,] 1 1
想要的结果：
1 1 3 3 5 5 1 3 9 1
我的尝试
我想出的最好办法是基于 duplicated 和 for 循环的复杂方法，既不高效也不优雅。我也知道可能
# Identify duplicates
duplicate <- duplicated(x, MARGIN = 1)
# Identify first occurrence of each duplicate
firstDup <- duplicated(x, MARGIN = 1, fromLast = TRUE) & !duplicate
indices <- which(firstDup)
# Initialize index for unique rows
index <- seq_len(dim(x)[1])
cf <- duplicate
for (i in indices) {
# Duplicates must occur after first occurrence
cf[seq_len(i)] <- FALSE
dups <- apply(x[cf, , drop = FALSE], 1L, identical, x[i, ])
index[which(cf)[dups]] <- i
}
index
是否有使用 base R 的优雅解决方案？
这个怎么样？
l <- asplit(x, 1L)
match(l, l)
[1] 1 1 3 3 5 5 1 3 9 1
这里，我们使用asplit获取x行的列表l，match获取每行第一次出现的索引.
如果矩阵很大，则以下解决方案可能就足够了：
l <- do.call(paste, data.frame(x))
match(l, l)
[1] 1 1 3 3 5 5 1 3 9 1
TL;DR
对于大小相同但形状不同的整数矩阵（5e+06×2、5e+05×20、5000×2000），包含从 1 到 10 的整数，测试最快的 base 答案是 grouping/match，@alexis_laz 在 base 答案是 data.table::frank/match，尽管 grouping/match 在所有情况下都具有可比性，甚至优于 data.table 答案在 5000 x 2000 的情况下。
请注意，对于范围更大的双精度矩阵或整数矩阵，结果可能会有所不同，具体取决于 data.table 可用的线程数。 [TODO?]
背景
@MikaelJagan 的 asplit/match(<list>, <list>) ?match 警告：
Matching for lists is potentially very slow and best avoided except in simple cases.
鉴于 OP 具有“具有许多行和列的矩阵”，我们想将 asplit/match(<list>, <list>) 答案的性能与另一个 base 的答案进行比较答案：
- @Onyambu 的
paste/match(<chr>, <chr>) - @ThomasIsCoding 的
interaction/match(<int>, <int>) - @alexis_laz的
grouping/match(<int>, <int>)
我们将这些与一些非 base 答案作为基准（我们认识到 OP 仅要求 base）：
- @MikaelJagan 的
Rcpp - @Henrik 的
data.table回答：- 将
which = TRUE和mult = "first"传递给[.data.table的自连接； - 两种基于行排名的方法，根据处理关系的方式而有所不同：
frank(ties.method = "average")/match(<dbl>, <dbl>),frank(ties.method = "dense")/match(<int>, <int>).
- 将
设置
library(microbenchmark)
library(data.table)
getDTthreads() # 4
f_asplit <- function(x) {
l <- asplit(x, 1L)
match(l, l) }
f_paste <- function(x) {
s <- do.call(paste, as.data.frame(x))
match(s, s) }
f_interaction <- function(x) {
z <- as.integer(interaction(as.data.frame(x)))
match(z, z) }
f_grouping <- function(x) {
g <- do.call(grouping, as.data.frame(x))
o <- order(g, method = "radix")
e <- attr(g, "ends")
z <- rep.int(seq_along(e), c(e[1L], e[-1L] - e[-length(e)]))[o]
match(z, z) }
f_join <- function(x) {
d <- as.data.table(x)
d[d, on = names(d), mult = "first", which = TRUE] }
f_frank_average <- function(x) {
d <- as.data.table(x)
r <- frank(d, ties.method = "average")
match(r, r) }
f_frank_dense <- function(x) {
d <- as.data.table(x)
r <- frank(d, ties.method = "dense")
match(r, r) }
Rcpp::sourceCpp('<copy source code from @MikaelJagan\'s answer here>')
基准测试
行多，列少
我们首先使用 5e+06×2 整数矩阵评估性能：
set.seed(1L)
x <- matrix(sample(10L, size = 1e+07L, replace = TRUE), ncol = 2L)
microbenchmark(
f_asplit(x),
f_paste(x),
f_interaction(x),
f_grouping(x),
f_join(x),
f_frank_average(x),
f_frank_dense(x),
f_rcpp(x),
times = 10L,
check = "identical",
setup = gc(FALSE)
)
Unit: milliseconds
expr min lq mean median uq max neval
f_asplit(x) 17369.93905 18861.91195 19070.21298 19013.0180 19207.29194 22420.71085 10
f_paste(x) 502.63884 507.35077 509.01823 509.2443 511.72301 515.10083 10
f_interaction(x) 234.19311 236.52494 241.80098 238.7392 242.32923 259.75644 10
f_grouping(x) 182.25226 182.89358 187.09642 184.6124 187.10444 208.15532 10
f_join(x) 119.43460 120.86829 123.16607 122.9332 125.07169 128.44722 10
f_frank_average(x) 104.40150 107.53607 111.00268 108.5597 116.80375 121.83675 10
f_frank_dense(x) 86.60926 88.29555 91.42976 90.4716 92.32413 99.30659 10
f_rcpp(x) 459.02304 464.79855 472.43669 468.2492 470.25508 523.06734 10
f_asplit 比 base 替代方案慢两个数量级。 f_grouping 是最快的 base 答案，但 f_frank_dense 快了大约 2 倍（总体上最快）。
更少的行，更多的列
以上结果并未推广到所有整数矩阵输入。例如，f_interaction 与 ncol(x) 的比例非常差：如果 x 的每一列都有 u 个唯一元素，则可能的交互数量为 u^ncol(x)。
出于这个原因，我们进行了第二次基准测试，这次考虑的是行数较少 (5e+05) 列数较多 (20) 的矩阵。
set.seed(1L)
x <- matrix(sample(10L, size = 1e+07L, replace = TRUE), ncol = 20L)
f_interaction 的初始测试导致内存分配错误，因此它被排除在基准测试之外。
system.time(f_interaction(x))
Error: cannot allocate vector of size 7.5 Gb
Timing stopped at: 173.2 6.05 200.4
microbenchmark(
f_asplit(x),
f_paste(x),
## f_interaction(x),
f_grouping(x),
f_join(x),
f_frank_average(x),
f_frank_dense(x),
f_rcpp(x),
times = 10L,
check = "identical",
setup = gc(FALSE)
)
Unit: milliseconds
expr min lq mean median uq max neval
f_asplit(x) 5416.08762 5681.23523 5731.89246 5732.31779 5905.44517 5913.77141 10
f_paste(x) 592.92990 604.15083 629.31101 623.78679 637.81814 724.83871 10
f_grouping(x) 63.89522 64.14134 65.42723 65.11530 66.00557 68.06045 10
f_join(x) 340.73722 342.18096 353.35774 352.08861 359.88480 382.13480 10
f_frank_average(x) 69.90496 70.81840 72.29819 72.04409 73.11977 77.44347 10
f_frank_dense(x) 52.58033 53.33760 54.42029 54.01672 55.63532 56.99664 10
f_rcpp(x) 184096.21999 184816.36584 185774.76817 186218.58335 186696.31674 186781.24972 10
f_grouping 仍然是最快的 base 答案。值得注意的是，它现在比 f_paste 快一个数量级，只比 f_frank_dense.
更少的行，更多的列
我们执行了最终基准测试，排除了上一轮中最慢的答案（f_asplit 和 f_rcpp），现在考虑一个 5000×2000 整数矩阵：
set.seed(1L)
x <- matrix(sample(10L, size = 1e+07L, replace = TRUE), ncol = 2000L)
microbenchmark(
## f_asplit(x),
f_paste(x),
## f_interaction(x),
f_grouping(x),
f_join(x),
f_frank_average(x),
f_frank_dense(x),
## f_rcpp(x),
times = 10L,
check = "identical",
setup = gc(FALSE)
)
Unit: milliseconds
expr min lq mean median uq max neval
f_paste(x) 1067.47994 1075.45148 1083.17391 1080.72997 1089.74027 1102.45249 10
f_grouping(x) 19.24007 19.50026 19.86404 19.79002 20.25302 20.60127 10
f_join(x) 616.66706 621.29854 630.61460 628.16315 636.39097 650.16180 10
f_frank_average(x) 59.82007 61.41706 62.68610 62.99318 64.56520 64.88463 10
f_frank_dense(x) 58.03648 60.59857 63.50526 61.99278 66.03694 71.30638 10
现在 f_grouping 总体上最快，比 f_frank_dense 快大约 3 倍。
如果您使用的是 base R
，我们可以使用ave
> ave(1:nrow(x), x[, 1], x[, 2], FUN = function(v) v[1])
[1] 1 1 3 3 5 5 1 3 9 1
如果你有多个列，你可以试试
> z <- as.integer(interaction(as.data.frame(m2)))
> ave(seq_along(z), z, FUN = function(x) x[1])
[1] 1 1 3 3 5 5 1 3 9 1
或
> z <- as.integer(interaction(as.data.frame(x)))
> match(z, z)
[1] 1 1 3 3 5 5 1 3 9 1
基准测试
set.seed(1)
v <- sample(1:10, 1e7, replace = TRUE)
m2 <- matrix(v, ncol = 2)
microbenchmark(
ave1 = {
ave(1:nrow(m2), m2[, 1], m2[, 2], FUN = function(v) v[1])
},
ave2 = {
z <- as.integer(interaction(as.data.frame(m2)))
ave(seq_along(z), z, FUN = function(x) x[1])
},
match = {
z <- as.integer(interaction(as.data.frame(m2)))
match(z, z)
},
times = 10L
)
我们会看到
Unit: milliseconds
expr min lq mean median uq max neval
ave1 648.0755 655.9521 715.8848 701.1927 747.4759 885.9838 10
ave2 785.4868 883.2935 913.3867 899.1789 929.6571 1050.9020 10
match 417.1598 447.3718 507.0462 495.8791 551.9436 625.0841 10
不是基础 R，而是另一个比较点：这是一个天真的 Rcpp 实现。最坏的情况是O(n * m^2 / 2)，其中m = nrow(x)和n = ncol(x)，这实际上是相当糟糕的。在行上并行化循环可能会有很大帮助。
它的优点是不复制矩阵x，所以根本不需要太多内存。目前给出的答案调用asplit、as.data.frame或as.data.table，全部复制（分区）x.
Rcpp::sourceCpp(code = '
#include <Rcpp.h>
using namespace Rcpp;
// [[Rcpp::export]]
IntegerVector f_rcpp(IntegerMatrix x)
{
int m = x.nrow();
int n = x.ncol();
IntegerVector res(m);
if (n == 0) {
res.fill(1);
} else {
int i, ki, p, kp, j;
for (i = 0; i < m; ++i) {
for (p = 0; p < i; ++p) {
for (j = 0, ki = i, kp = p; j < n; ++j) {
if (x[kp] != x[ki]) {
break;
}
ki += m;
kp += m;
}
if (j == n) {
res[i] = p + 1;
break;
}
}
if (p == i) {
res[i] = i + 1;
}
}
}
return res;
}
')
f_rcpp(x)
[1] 1 1 3 3 5 5 1 3 9 1