如何按组加速子集
How to speed up subset by groups
我曾经使用 dplyr 来实现我的数据整理，但有些计算是 "slow"。特别是按组划分的子集，我读到当有很多组时 dplyr 很慢并且基于 this benchmark data.table 可能会更快所以我开始学习 data.table.
这里是如何使用 250k 行和大约 230k 组重现接近我的真实数据的东西。我想按 id1、id2 分组，并为每个组使用 max(datetime) 对行进行子集化。
数据
# random datetime generation function by Dirk Eddelbuettel
#
rand.datetime <- function(N, st = "2012/01/01", et = "2015/08/05") {
st <- as.POSIXct(as.Date(st))
et <- as.POSIXct(as.Date(et))
dt <- as.numeric(difftime(et,st,unit="sec"))
ev <- sort(runif(N, 0, dt))
rt <- st + ev
}
set.seed(42)
# Creating 230000 ids couples
ids <- data.frame(id1 = stringi::stri_rand_strings(23e4, 9, pattern = "[0-9]"),
id2 = stringi::stri_rand_strings(23e4, 9, pattern = "[0-9]"))
# Repeating randomly the ids[1:2000, ] to create groups
ids <- rbind(ids, ids[sample(1:2000, 20000, replace = TRUE), ])
# Adding random datetime variable and dummy variables to reproduce real datas
datas <- transform(ids,
datetime = rand.datetime(25e4),
var1 = sample(LETTERS[1:6], 25e4, rep = TRUE),
var2 = sample(c(1:10, NA), 25e4, rep = TRUE),
var3 = sample(c(1:10, NA), 25e4, rep = TRUE),
var4 = rand.datetime(25e4),
var5 = rand.datetime(25e4))
datas.tbl <- tbl_df(datas)
datas.dt <- data.table(datas, key = c("id1", "id2"))
我找不到使用 data.table 按组划分子集的直接方法，所以我问了这个问题：Filter rows by groups with data.table
我们建议我使用 .SD :
datas.dt[, .SD[datetime == max(datetime)], by = c("id1", "id2")]
但我有两个问题，它适用于日期但不适用于 POSIXct ("Error in UseMethod("as.data.table")：
没有适用于 'as.data.table' 的方法应用于 class "c('POSIXct', 'POSIXt')"" 的对象），这非常慢。例如日期：
> system.time({
+ datas.dt[, .SD[as.Date(datetime) == max(as.Date(datetime))], by = c("id1", "id2")]
+ })
utilisateur système écoulé
207.03 0.00 207.48
所以我找到了其他更快的方法来实现这个（并保持日期时间） data.table :
函数
f.dplyr <- function(x) x %>% group_by(id1, id2) %>% filter(datetime == max(datetime))
f.dt.i <- function(x) x[x[, .I[datetime == max(datetime)], by = c("id1", "id2")]$V1]
f.dt <- function(x) x[x[, datetime == max(datetime), by = c("id1", "id2")]$V1]
但后来我觉得data.table会快很多，和dplyr的时间差没有意义。
微基准测试
mbm <- microbenchmark(
dplyr = res1 <- f.dplyr(datas.tbl),
data.table.I = res2 <- f.dt.i(datas.dt),
data.table = res3 <- f.dt(datas.dt),
times = 50L)
Unit: seconds
expr min lq mean median uq max neval
dplyr 31.84249 32.24055 32.59046 32.61311 32.88703 33.54226 50
data.table.I 30.02831 30.94621 31.19660 31.17820 31.42888 32.16521 50
data.table 30.28923 30.84212 31.09749 31.04851 31.40432 31.96351 50
我 missing/misusing 和 data.table 有什么关系吗？您有加快计算速度的想法吗？
非常感谢任何帮助！谢谢
编辑：关于用于微基准测试的系统和包版本的一些精确度。 （电脑不是war机器，12Go i5）
系统
sessionInfo()
R version 3.1.3 (2015-03-09)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 7 x64 (build 7601) Service Pack 1
locale:
[1] LC_COLLATE=French_France.1252 LC_CTYPE=French_France.1252
[3] LC_MONETARY=French_France.1252 LC_NUMERIC=C
[5] LC_TIME=French_France.1252
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] readr_0.1.0 ggplot2_1.0.1 microbenchmark_1.4-2
[4] data.table_1.9.4 dplyr_0.4.1 plyr_1.8.2
loaded via a namespace (and not attached):
[1] assertthat_0.1 chron_2.3-45 colorspace_1.2-6 DBI_0.3.1
[5] digest_0.6.8 grid_3.1.3 gtable_0.1.2 lazyeval_0.1.10
[9] magrittr_1.5 MASS_7.3-39 munsell_0.4.2 parallel_3.1.3
[13] proto_0.3-10 Rcpp_0.11.5 reshape2_1.4.1 scales_0.2.4
[17] stringi_0.4-1 stringr_0.6.2 tools_3.1.3
> packageVersion("data.table")
[1] ‘1.9.4’
> packageVersion("dplyr")
[1] ‘0.4.1’
如何总结data.table和join原始数据
system.time({
datas1 <- datas.dt[, list(datetime=max(datetime)), by = c("id1", "id2")] #summarize the data
setkey(datas1, id1, id2, datetime)
setkey(datas.dt, id1, id2, datetime)
datas2 <- datas.dt[datas1]
})
# user system elapsed
# 0.083 0.000 0.084
正确过滤数据
system.time(dat1 <- datas.dt[datas.dt[, .I[datetime == max(datetime)], by = c("id1", "id2")]$V1])
# user system elapsed
# 23.226 0.000 23.256
all.equal(dat1, datas2)
# [1] TRUE
附录
如果您使用 data.table 的 devel version，setkey 参数是多余的（感谢@akrun 的指点）
system.time({
datas1 <- datas.dt[, list(datetime=max(datetime)), by = c("id1", "id2")] #summarize the data
datas2 <- datas.dt[datas1, on=c('id1', 'id2', 'datetime')]
})
好问题！
我假设 df 和 dt 是 easy/quick 键入的对象名称。
df = datas.tbl
dt = datas.dt
-O3级优化比较：
首先，这是我的系统在 dplyr 的当前 CRAN 版本和 data.table 的开发版本上的时间。 dplyr 的开发版本似乎受到性能回归的影响（并且正在被 Romain 修复）。
system.time(df %>% group_by(id1, id2) %>% filter(datetime == max(datetime)))
# 25.291 0.128 25.610
system.time(dt[dt[, .I[datetime == max(datetime)], by = c("id1", "id2")]$V1])
# 17.191 0.075 17.349
我运行这好几次了，力好像变了。但是，我使用 -O3 优化标志编译所有包（通过适当设置 ~/.R/Makevars）。而且我观察到 data.table 性能比我在 -O3.
比较过的其他软件包要好得多
分组速度比较
其次，了解这种缓慢的原因很重要。首先让我们将时间与组.
进行比较
system.time(group_by(df, id1, id2))
# 0.303 0.007 0.311
system.time(data.table:::forderv(dt, by = c("id1", "id2"), retGrp = TRUE))
# 0.002 0.000 0.002
即使总共有 250,000 行，您的数据大小也约为 38MB。在此大小下，分组速度不太可能出现明显差异。
data.table的分组在这里>100x更快，这显然不是这么慢的原因...
为什么这么慢？
那么是什么原因呢？让我们打开 datatable.verbose 选项并再次检查：
options(datatable.verbose = TRUE)
dt[dt[, .I[datetime == max(datetime)], by = c("id1", "id2")]$V1]
# Detected that j uses these columns: datetime
# Finding groups (bysameorder=TRUE) ... done in 0.002secs. bysameorder=TRUE and o__ is length 0
# lapply optimization is on, j unchanged as '.I[datetime == max(datetime)]'
# GForce is on, left j unchanged
# Old mean optimization is on, left j unchanged.
# Starting dogroups ...
# memcpy contiguous groups took 0.097s for 230000 groups
# eval(j) took 17.129s for 230000 calls
# done dogroups in 17.597 secs
所以仅 eval(j) 就占用了 ~97% 的时间！我们在 j 中提供的表达式针对 每个组 进行评估。由于您有 230,000 个组，并且 eval() 调用会受到惩罚，所以加起来。
避免eval()惩罚
自从我们意识到这种惩罚后，我们就开始实施一些常用函数的内部版本：sum、mean、min、max。这个 will/should 可以扩展到尽可能多的其他功能（当我们找到时间时）。
那么，我们先来计算一下刚刚获得max(datetime)的时间：
dt.agg = dt[, .(datetime = max(datetime)), by = .(id1, id2)]
# Detected that j uses these columns: datetime
# Finding groups (bysameorder=TRUE) ... done in 0.002secs. bysameorder=TRUE and o__ is length 0
# lapply optimization is on, j unchanged as 'list(max(datetime))'
# GForce optimized j to 'list(gmax(datetime))'
而且是即时的。为什么？因为 max() 在内部优化为 gmax() 并且没有 eval() 调用每个 230K 组。
那么为什么 datetime == max(datetime) 不是即时的？因为解析这样的表达式并在内部进行优化比较复杂，我们还没有做到。
解决方法
现在我们知道了这个问题，并且找到了解决它的方法，让我们使用它。
dt.agg = dt[, .(datetime = max(datetime)), by = .(id1, id2)]
dt[dt.agg, on = c("id1", "id2", "datetime")] # v1.9.5+
我的 Mac.
这需要大约 0.14 秒
请注意，这只是速度很快 ，因为 表达式已优化为 gmax()。比较：
dt[, .(datetime = base::max(datetime)), by = .(id1, id2)]
我同意优化更复杂的表达式以避免 eval() 惩罚是理想的解决方案，但我们还没有做到这一点。
我曾经使用 dplyr 来实现我的数据整理，但有些计算是 "slow"。特别是按组划分的子集，我读到当有很多组时 dplyr 很慢并且基于 this benchmark data.table 可能会更快所以我开始学习 data.table.
这里是如何使用 250k 行和大约 230k 组重现接近我的真实数据的东西。我想按 id1、id2 分组，并为每个组使用 max(datetime) 对行进行子集化。
数据
# random datetime generation function by Dirk Eddelbuettel
#
rand.datetime <- function(N, st = "2012/01/01", et = "2015/08/05") {
st <- as.POSIXct(as.Date(st))
et <- as.POSIXct(as.Date(et))
dt <- as.numeric(difftime(et,st,unit="sec"))
ev <- sort(runif(N, 0, dt))
rt <- st + ev
}
set.seed(42)
# Creating 230000 ids couples
ids <- data.frame(id1 = stringi::stri_rand_strings(23e4, 9, pattern = "[0-9]"),
id2 = stringi::stri_rand_strings(23e4, 9, pattern = "[0-9]"))
# Repeating randomly the ids[1:2000, ] to create groups
ids <- rbind(ids, ids[sample(1:2000, 20000, replace = TRUE), ])
# Adding random datetime variable and dummy variables to reproduce real datas
datas <- transform(ids,
datetime = rand.datetime(25e4),
var1 = sample(LETTERS[1:6], 25e4, rep = TRUE),
var2 = sample(c(1:10, NA), 25e4, rep = TRUE),
var3 = sample(c(1:10, NA), 25e4, rep = TRUE),
var4 = rand.datetime(25e4),
var5 = rand.datetime(25e4))
datas.tbl <- tbl_df(datas)
datas.dt <- data.table(datas, key = c("id1", "id2"))
我找不到使用 data.table 按组划分子集的直接方法，所以我问了这个问题：Filter rows by groups with data.table
我们建议我使用 .SD :
datas.dt[, .SD[datetime == max(datetime)], by = c("id1", "id2")]
但我有两个问题，它适用于日期但不适用于 POSIXct ("Error in UseMethod("as.data.table")： 没有适用于 'as.data.table' 的方法应用于 class "c('POSIXct', 'POSIXt')"" 的对象），这非常慢。例如日期：
> system.time({
+ datas.dt[, .SD[as.Date(datetime) == max(as.Date(datetime))], by = c("id1", "id2")]
+ })
utilisateur système écoulé
207.03 0.00 207.48
所以我找到了其他更快的方法来实现这个（并保持日期时间） data.table :
函数
f.dplyr <- function(x) x %>% group_by(id1, id2) %>% filter(datetime == max(datetime))
f.dt.i <- function(x) x[x[, .I[datetime == max(datetime)], by = c("id1", "id2")]$V1]
f.dt <- function(x) x[x[, datetime == max(datetime), by = c("id1", "id2")]$V1]
但后来我觉得data.table会快很多，和dplyr的时间差没有意义。
微基准测试
mbm <- microbenchmark(
dplyr = res1 <- f.dplyr(datas.tbl),
data.table.I = res2 <- f.dt.i(datas.dt),
data.table = res3 <- f.dt(datas.dt),
times = 50L)
Unit: seconds
expr min lq mean median uq max neval
dplyr 31.84249 32.24055 32.59046 32.61311 32.88703 33.54226 50
data.table.I 30.02831 30.94621 31.19660 31.17820 31.42888 32.16521 50
data.table 30.28923 30.84212 31.09749 31.04851 31.40432 31.96351 50
我 missing/misusing 和 data.table 有什么关系吗？您有加快计算速度的想法吗？
非常感谢任何帮助！谢谢
编辑：关于用于微基准测试的系统和包版本的一些精确度。 （电脑不是war机器，12Go i5）
系统
sessionInfo()
R version 3.1.3 (2015-03-09)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 7 x64 (build 7601) Service Pack 1
locale:
[1] LC_COLLATE=French_France.1252 LC_CTYPE=French_France.1252
[3] LC_MONETARY=French_France.1252 LC_NUMERIC=C
[5] LC_TIME=French_France.1252
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] readr_0.1.0 ggplot2_1.0.1 microbenchmark_1.4-2
[4] data.table_1.9.4 dplyr_0.4.1 plyr_1.8.2
loaded via a namespace (and not attached):
[1] assertthat_0.1 chron_2.3-45 colorspace_1.2-6 DBI_0.3.1
[5] digest_0.6.8 grid_3.1.3 gtable_0.1.2 lazyeval_0.1.10
[9] magrittr_1.5 MASS_7.3-39 munsell_0.4.2 parallel_3.1.3
[13] proto_0.3-10 Rcpp_0.11.5 reshape2_1.4.1 scales_0.2.4
[17] stringi_0.4-1 stringr_0.6.2 tools_3.1.3
> packageVersion("data.table")
[1] ‘1.9.4’
> packageVersion("dplyr")
[1] ‘0.4.1’
如何总结data.table和join原始数据
system.time({
datas1 <- datas.dt[, list(datetime=max(datetime)), by = c("id1", "id2")] #summarize the data
setkey(datas1, id1, id2, datetime)
setkey(datas.dt, id1, id2, datetime)
datas2 <- datas.dt[datas1]
})
# user system elapsed
# 0.083 0.000 0.084
正确过滤数据
system.time(dat1 <- datas.dt[datas.dt[, .I[datetime == max(datetime)], by = c("id1", "id2")]$V1])
# user system elapsed
# 23.226 0.000 23.256
all.equal(dat1, datas2)
# [1] TRUE
附录
如果您使用data.table 的 devel version，setkey 参数是多余的（感谢@akrun 的指点）
system.time({
datas1 <- datas.dt[, list(datetime=max(datetime)), by = c("id1", "id2")] #summarize the data
datas2 <- datas.dt[datas1, on=c('id1', 'id2', 'datetime')]
})
好问题！
我假设 df 和 dt 是 easy/quick 键入的对象名称。
df = datas.tbl
dt = datas.dt
-O3级优化比较：
首先，这是我的系统在 dplyr 的当前 CRAN 版本和 data.table 的开发版本上的时间。 dplyr 的开发版本似乎受到性能回归的影响（并且正在被 Romain 修复）。
system.time(df %>% group_by(id1, id2) %>% filter(datetime == max(datetime)))
# 25.291 0.128 25.610
system.time(dt[dt[, .I[datetime == max(datetime)], by = c("id1", "id2")]$V1])
# 17.191 0.075 17.349
我运行这好几次了，力好像变了。但是，我使用 -O3 优化标志编译所有包（通过适当设置 ~/.R/Makevars）。而且我观察到 data.table 性能比我在 -O3.
分组速度比较
其次，了解这种缓慢的原因很重要。首先让我们将时间与组.
进行比较system.time(group_by(df, id1, id2))
# 0.303 0.007 0.311
system.time(data.table:::forderv(dt, by = c("id1", "id2"), retGrp = TRUE))
# 0.002 0.000 0.002
即使总共有 250,000 行，您的数据大小也约为 38MB。在此大小下，分组速度不太可能出现明显差异。
data.table的分组在这里>100x更快，这显然不是这么慢的原因...
为什么这么慢？
那么是什么原因呢？让我们打开 datatable.verbose 选项并再次检查：
options(datatable.verbose = TRUE)
dt[dt[, .I[datetime == max(datetime)], by = c("id1", "id2")]$V1]
# Detected that j uses these columns: datetime
# Finding groups (bysameorder=TRUE) ... done in 0.002secs. bysameorder=TRUE and o__ is length 0
# lapply optimization is on, j unchanged as '.I[datetime == max(datetime)]'
# GForce is on, left j unchanged
# Old mean optimization is on, left j unchanged.
# Starting dogroups ...
# memcpy contiguous groups took 0.097s for 230000 groups
# eval(j) took 17.129s for 230000 calls
# done dogroups in 17.597 secs
所以仅 eval(j) 就占用了 ~97% 的时间！我们在 j 中提供的表达式针对 每个组 进行评估。由于您有 230,000 个组，并且 eval() 调用会受到惩罚，所以加起来。
避免eval()惩罚
自从我们意识到这种惩罚后，我们就开始实施一些常用函数的内部版本：sum、mean、min、max。这个 will/should 可以扩展到尽可能多的其他功能（当我们找到时间时）。
那么，我们先来计算一下刚刚获得max(datetime)的时间：
dt.agg = dt[, .(datetime = max(datetime)), by = .(id1, id2)]
# Detected that j uses these columns: datetime
# Finding groups (bysameorder=TRUE) ... done in 0.002secs. bysameorder=TRUE and o__ is length 0
# lapply optimization is on, j unchanged as 'list(max(datetime))'
# GForce optimized j to 'list(gmax(datetime))'
而且是即时的。为什么？因为 max() 在内部优化为 gmax() 并且没有 eval() 调用每个 230K 组。
那么为什么 datetime == max(datetime) 不是即时的？因为解析这样的表达式并在内部进行优化比较复杂，我们还没有做到。
解决方法
现在我们知道了这个问题，并且找到了解决它的方法，让我们使用它。
dt.agg = dt[, .(datetime = max(datetime)), by = .(id1, id2)]
dt[dt.agg, on = c("id1", "id2", "datetime")] # v1.9.5+
我的 Mac.
这需要大约 0.14 秒请注意，这只是速度很快 ，因为 表达式已优化为 gmax()。比较：
dt[, .(datetime = base::max(datetime)), by = .(id1, id2)]
我同意优化更复杂的表达式以避免 eval() 惩罚是理想的解决方案，但我们还没有做到这一点。