data.table 分组操作的性能
Performance of Grouped Operations with data.table
我将从描述我正在执行的任务开始。我必须重复计算分组总和，通常在 5 到 10 次之间。我正在执行分组总和的列中的值随每次迭代而变化，但我分组所依据的列却没有。下面是一个例子table，其中w、x、y一起构成了分组，z是要求和的列。
DT <- data.table(w = sample(1:10, size = 1000000, replace = TRUE),
x = sample(1:100, size = 1000000, replace = TRUE),
y = sample(1:1000, size = 1000000, replace = TRUE),
z = runif(n = 1000000, min = 0, max = 1))
setkey(DT, w, x, y)
我认为我最初的解决方案是最明显的：
DT[, Group_Total := sum(z), keyby = list(w, x, y)]
microbenchmark(DT[, Group_Total := sum(z), keyby = list(w, x, y)])
Unit: milliseconds
min lq mean median uq max neval
447.8674 467.3366 481.1989 479.3057 490.5499 691 100
我发现这个解决方案并没有我希望的那么快，特别是因为我不得不多次执行这个计算。但是，据我所知，data.table 中没有一种方法可以提供比这更好的性能。
最近，我开始尝试使用 Rccp 并编写了以下函数来执行分组求和：
cppFunction('NumericVector Group_Total(NumericVector Grouping, NumericVector x) {
int n = Grouping.size();
NumericVector out(n);
int curr_start_index = 0;
double grp_total = x[0];
for(int i = 1; i < (n-1); ++i) {
if(Grouping[i] == 1) {
for(int j = curr_start_index; j < i; ++j) {
out[j] = grp_total;
}
curr_start_index = i;
grp_total = x[i];
} else {
grp_total = grp_total + x[i];
}
}
if(Grouping[(n-1)] == 1) {
for(int j = curr_start_index; j < (n-1); ++j) {
out[j] = grp_total;
}
out[(n-1)] = x[(n-1)];
} else {
grp_total = grp_total + x[(n-1)];
for(int j = curr_start_index; j < n; ++j) {
out[j] = grp_total;
}
}
return out;
}')
为了使用这个函数，我在数据排序后计算一个组内索引列：
DT[, Within_Group_Index := 1][, Within_Group_Index := cumsum(Within_Group_Index), keyby = list(w, x, y)]
由于我只需要计算一次这一列，所以增加的时间并不重要。创建此列后，可以通过以下方式计算分组总和：
DT[, Group_Total := Group_Total(Within_Group_Index, z)]
microbenchmark(DT[, Group_Total := Group_Total(Within_Group_Index, z)])
Unit: milliseconds
min lq mean median uq max neval
5.3286 6.9879 9.078267 7.4814 8.0171 161.6406 100
这大大加快了速度，现在我们来回答我的问题：为什么？
我的理解是分组操作在data.table中非常高效。我可以做些什么来加速 data.table 分组操作，或者这是基本速度与灵活性权衡的示例，即 data.table 分组操作的灵活性是以成本为代价的考虑特定用例时的速度？
我认为你的怀疑是正确的。
至于为什么，我们可以做出有根据的猜测。
编辑：以下信息忽略了 data.table 对 GForce 优化的作用，
GForce 支持的函数可能会避免类似于 C++ 代码的数据副本，但是，
正如弗兰克所说，
:= 在 1.14.3 之前没有尝试检测 GForce-optimizable 表达式，
而下面的addrs肯定没有被GForce优化。
您编写的 C++ 代码不会修改其 输入 数据的任何内容，
它只需要分配 out。
正如你所说的那样，
data.table 旨在更加灵活，
并且必须支持用户提供的任何有效 R 代码。
这让我觉得，
对于分组操作，
它必须根据组索引为 z 的每个子集分配新的 R 向量，
多次有效地复制一些输入。
我想尝试验证我的假设，
所以我写了一些 C++ 代码来查看内存地址：
set.seed(31L)
DT <- data.table(w = sample(1:3, size = 20, replace = TRUE),
x = sample(1:3, size = 20, replace = TRUE),
y = sample(1:3, size = 20, replace = TRUE),
z = runif(n = 20, min = 0, max = 1))
setkey(DT, w, x, y)
cppFunction(plugins = "cpp11", includes = "#include <sstream>", '
StringVector addrs(NumericVector z, IntegerVector indices, IntegerVector group_ids) {
StringVector ans(indices.size());
for (auto i = 0; i < indices.size(); ++i) {
std::ostringstream oss;
oss << "Group" << group_ids[i] << " " << &z[indices[i]];
ans[i] = oss.str();
}
return ans;
}
')
idx <- DT[, .(indices = .I[1L] - 1L, group_ids = .GRP), keyby = list(w, x, y)]
DT[, list(addrs(z, idx$indices, idx$group_ids))]
# V1
# 1: Group1 0x55b7733361f8
# 2: Group2 0x55b773336200
# 3: Group3 0x55b773336210
# 4: Group4 0x55b773336220
# 5: Group5 0x55b773336228
# 6: Group6 0x55b773336230
# 7: Group7 0x55b773336238
# 8: Group8 0x55b773336248
# 9: Group9 0x55b773336250
# 10: Group10 0x55b773336258
# 11: Group11 0x55b773336260
# 12: Group12 0x55b773336270
# 13: Group13 0x55b773336280
# 14: Group14 0x55b773336288
# 15: Group15 0x55b773336290
正如预期的那样，
如果我们把 z 作为一个整体来看，
没有复制发生，不同元素的地址彼此接近，
例如 0x55b773336200 = 0x55b7733361f8 + 0x8。
您可以多次执行前面代码的最后一行，它总是会显示相同的地址。
我部分没想到的是：
DT[, list(addrs(z, 0L, .GRP)), keyby = list(w, x, y)]
# w x y V1
# 1: 1 1 2 Group1 0x55b771b03440
# 2: 1 1 3 Group2 0x55b771b03440
# 3: 1 2 1 Group3 0x55b771b03440
# 4: 1 2 2 Group4 0x55b771b03440
# 5: 1 3 2 Group5 0x55b771b03440
# 6: 1 3 3 Group6 0x55b771b03440
# 7: 2 1 1 Group7 0x55b771b03440
# 8: 2 2 1 Group8 0x55b771b03440
# 9: 2 2 2 Group9 0x55b771b03440
# 10: 2 2 3 Group10 0x55b771b03440
# 11: 2 3 1 Group11 0x55b771b03440
# 12: 3 2 1 Group12 0x55b771b03440
# 13: 3 2 2 Group13 0x55b771b03440
# 14: 3 2 3 Group14 0x55b771b03440
# 15: 3 3 3 Group15 0x55b771b03440
一方面，
内存中的地址确实改变了，
所以有些东西被复制了，
如果你 运行 多次，你每次都会看到不同的地址。
但是，似乎 data.table 以某种方式重用了具有相同地址的缓冲区，
也许分配一个长度最大的 group-subset 的数组并将不同的组值复制到它？
我想知道他们是如何管理的。
或者也许我的代码是错误的¯\_(ツ)_/¯
编辑：我在 addrs 的第一行添加了以下内容
（由于 R 解析的双重转义）
Rcout << "input length = " << z.size() << "\n";
如果我 运行 上面的最后一个 DT[...] 代码，
即使地址相同，它也会打印不同的长度。
我发现初始代码在开发分支版本 1.14.3 上更快。来自 the NEWS:
- := is now optimized by group
我相信这适用于 GForce 工作的任何地方。参见 ?GForce。 （感谢@Alexis 指出这一点）
我看到的时间和你的差不多，但是升级后发现原来的代码还是挺快的：
microbenchmark(
direct = DT[, Group_Total := sum(z), keyby = list(w, x, y)],
indirect = DT[, idx := rowid(w,x,y)][, gt := Group_Total(idx, z)]
)
Unit: milliseconds
expr min lq mean median uq max neval
direct 42.92247 44.61097 48.77803 45.48076 46.93751 72.29237 100
indirect 30.49422 31.95343 34.69138 33.18616 35.32934 58.37454 100
备注：
rowid 构建 within-group 索引。DT[, Group_Total := sum(z), keyby = list(w, x, y), verbose=TRUE] 将显示有关是否使用 GForce 优化的信息。
我将从描述我正在执行的任务开始。我必须重复计算分组总和，通常在 5 到 10 次之间。我正在执行分组总和的列中的值随每次迭代而变化，但我分组所依据的列却没有。下面是一个例子table，其中w、x、y一起构成了分组，z是要求和的列。
DT <- data.table(w = sample(1:10, size = 1000000, replace = TRUE),
x = sample(1:100, size = 1000000, replace = TRUE),
y = sample(1:1000, size = 1000000, replace = TRUE),
z = runif(n = 1000000, min = 0, max = 1))
setkey(DT, w, x, y)
我认为我最初的解决方案是最明显的：
DT[, Group_Total := sum(z), keyby = list(w, x, y)]
microbenchmark(DT[, Group_Total := sum(z), keyby = list(w, x, y)])
Unit: milliseconds
min lq mean median uq max neval
447.8674 467.3366 481.1989 479.3057 490.5499 691 100
我发现这个解决方案并没有我希望的那么快，特别是因为我不得不多次执行这个计算。但是，据我所知，data.table 中没有一种方法可以提供比这更好的性能。 最近，我开始尝试使用 Rccp 并编写了以下函数来执行分组求和：
cppFunction('NumericVector Group_Total(NumericVector Grouping, NumericVector x) {
int n = Grouping.size();
NumericVector out(n);
int curr_start_index = 0;
double grp_total = x[0];
for(int i = 1; i < (n-1); ++i) {
if(Grouping[i] == 1) {
for(int j = curr_start_index; j < i; ++j) {
out[j] = grp_total;
}
curr_start_index = i;
grp_total = x[i];
} else {
grp_total = grp_total + x[i];
}
}
if(Grouping[(n-1)] == 1) {
for(int j = curr_start_index; j < (n-1); ++j) {
out[j] = grp_total;
}
out[(n-1)] = x[(n-1)];
} else {
grp_total = grp_total + x[(n-1)];
for(int j = curr_start_index; j < n; ++j) {
out[j] = grp_total;
}
}
return out;
}')
为了使用这个函数，我在数据排序后计算一个组内索引列：
DT[, Within_Group_Index := 1][, Within_Group_Index := cumsum(Within_Group_Index), keyby = list(w, x, y)]
由于我只需要计算一次这一列，所以增加的时间并不重要。创建此列后，可以通过以下方式计算分组总和：
DT[, Group_Total := Group_Total(Within_Group_Index, z)]
microbenchmark(DT[, Group_Total := Group_Total(Within_Group_Index, z)])
Unit: milliseconds
min lq mean median uq max neval
5.3286 6.9879 9.078267 7.4814 8.0171 161.6406 100
这大大加快了速度，现在我们来回答我的问题：为什么？ 我的理解是分组操作在data.table中非常高效。我可以做些什么来加速 data.table 分组操作，或者这是基本速度与灵活性权衡的示例，即 data.table 分组操作的灵活性是以成本为代价的考虑特定用例时的速度？
我认为你的怀疑是正确的。 至于为什么，我们可以做出有根据的猜测。
编辑：以下信息忽略了 data.table 对 GForce 优化的作用，
GForce 支持的函数可能会避免类似于 C++ 代码的数据副本，但是，
正如弗兰克所说，
:= 在 1.14.3 之前没有尝试检测 GForce-optimizable 表达式，
而下面的addrs肯定没有被GForce优化。
您编写的 C++ 代码不会修改其 输入 数据的任何内容，
它只需要分配 out。
正如你所说的那样，
data.table 旨在更加灵活，
并且必须支持用户提供的任何有效 R 代码。
这让我觉得，
对于分组操作，
它必须根据组索引为 z 的每个子集分配新的 R 向量，
多次有效地复制一些输入。
我想尝试验证我的假设， 所以我写了一些 C++ 代码来查看内存地址：
set.seed(31L)
DT <- data.table(w = sample(1:3, size = 20, replace = TRUE),
x = sample(1:3, size = 20, replace = TRUE),
y = sample(1:3, size = 20, replace = TRUE),
z = runif(n = 20, min = 0, max = 1))
setkey(DT, w, x, y)
cppFunction(plugins = "cpp11", includes = "#include <sstream>", '
StringVector addrs(NumericVector z, IntegerVector indices, IntegerVector group_ids) {
StringVector ans(indices.size());
for (auto i = 0; i < indices.size(); ++i) {
std::ostringstream oss;
oss << "Group" << group_ids[i] << " " << &z[indices[i]];
ans[i] = oss.str();
}
return ans;
}
')
idx <- DT[, .(indices = .I[1L] - 1L, group_ids = .GRP), keyby = list(w, x, y)]
DT[, list(addrs(z, idx$indices, idx$group_ids))]
# V1
# 1: Group1 0x55b7733361f8
# 2: Group2 0x55b773336200
# 3: Group3 0x55b773336210
# 4: Group4 0x55b773336220
# 5: Group5 0x55b773336228
# 6: Group6 0x55b773336230
# 7: Group7 0x55b773336238
# 8: Group8 0x55b773336248
# 9: Group9 0x55b773336250
# 10: Group10 0x55b773336258
# 11: Group11 0x55b773336260
# 12: Group12 0x55b773336270
# 13: Group13 0x55b773336280
# 14: Group14 0x55b773336288
# 15: Group15 0x55b773336290
正如预期的那样，
如果我们把 z 作为一个整体来看，
没有复制发生，不同元素的地址彼此接近，
例如 0x55b773336200 = 0x55b7733361f8 + 0x8。
您可以多次执行前面代码的最后一行，它总是会显示相同的地址。
我部分没想到的是：
DT[, list(addrs(z, 0L, .GRP)), keyby = list(w, x, y)]
# w x y V1
# 1: 1 1 2 Group1 0x55b771b03440
# 2: 1 1 3 Group2 0x55b771b03440
# 3: 1 2 1 Group3 0x55b771b03440
# 4: 1 2 2 Group4 0x55b771b03440
# 5: 1 3 2 Group5 0x55b771b03440
# 6: 1 3 3 Group6 0x55b771b03440
# 7: 2 1 1 Group7 0x55b771b03440
# 8: 2 2 1 Group8 0x55b771b03440
# 9: 2 2 2 Group9 0x55b771b03440
# 10: 2 2 3 Group10 0x55b771b03440
# 11: 2 3 1 Group11 0x55b771b03440
# 12: 3 2 1 Group12 0x55b771b03440
# 13: 3 2 2 Group13 0x55b771b03440
# 14: 3 2 3 Group14 0x55b771b03440
# 15: 3 3 3 Group15 0x55b771b03440
一方面，
内存中的地址确实改变了，
所以有些东西被复制了，
如果你 运行 多次，你每次都会看到不同的地址。
但是，似乎 data.table 以某种方式重用了具有相同地址的缓冲区，
也许分配一个长度最大的 group-subset 的数组并将不同的组值复制到它？
我想知道他们是如何管理的。
或者也许我的代码是错误的¯\_(ツ)_/¯
编辑：我在 addrs 的第一行添加了以下内容
（由于 R 解析的双重转义）
Rcout << "input length = " << z.size() << "\n";
如果我 运行 上面的最后一个 DT[...] 代码，
即使地址相同，它也会打印不同的长度。
我发现初始代码在开发分支版本 1.14.3 上更快。来自 the NEWS:
- := is now optimized by group
我相信这适用于 GForce 工作的任何地方。参见 ?GForce。 （感谢@Alexis 指出这一点）
我看到的时间和你的差不多，但是升级后发现原来的代码还是挺快的：
microbenchmark(
direct = DT[, Group_Total := sum(z), keyby = list(w, x, y)],
indirect = DT[, idx := rowid(w,x,y)][, gt := Group_Total(idx, z)]
)
Unit: milliseconds
expr min lq mean median uq max neval
direct 42.92247 44.61097 48.77803 45.48076 46.93751 72.29237 100
indirect 30.49422 31.95343 34.69138 33.18616 35.32934 58.37454 100
备注：
rowid构建 within-group 索引。DT[, Group_Total := sum(z), keyby = list(w, x, y), verbose=TRUE]将显示有关是否使用 GForce 优化的信息。