R：向量化循环以创建成对矩阵
R: Vectorize loop to create pairwise matrix
我想加速创建一个成对矩阵的函数，该矩阵描述了在一组位置中在所有其他对象之前和之后选择一个对象的次数。
这是一个例子df:
df <- data.frame(Shop = c("A","A","A","B","B","C","C","D","D","D","E","E","E"),
Fruit = c("apple", "orange", "pear",
"orange", "pear",
"pear", "apple",
"pear", "apple", "orange",
"pear", "apple", "orange"),
Order = c(1, 2, 3,
1, 2,
1, 2,
1, 2, 3,
1, 1, 1))
在每个 Shop 中，Fruit 是由客户在给定的 Order 中挑选的。
以下函数创建一个 m x n 成对矩阵：
loop.function <- function(df){
fruits <- unique(df$Fruit)
nt <- length(fruits)
mat <- array(dim=c(nt,nt))
for(m in 1:nt){
for(n in 1:nt){
## filter df for each pair of fruit
xm <- df[df$Fruit == fruits[m],]
xn <- df[df$Fruit == fruits[n],]
## index instances when a pair of fruit are picked in same shop
mm <- match(xm$Shop, xn$Shop)
## filter xm and xn based on mm
xm <- xm[! is.na(mm),]
xn <- xn[mm[! is.na(mm)],]
## assign number of times fruit[m] is picked after fruit[n] to mat[m,n]
mat[m,n] <- sum(xn$Order < xm$Order)
}
}
row.names(mat) <- fruits
colnames(mat) <- fruits
return(mat)
}
其中 mat[m,n] 是 fruits[m] 被选中的次数 在 fruits[n] 之后。而 mat[n,m] 是 fruits[m] 在 fruits[n] 之前 被选中的次数。如果同时采摘成对的水果则不记录（例如 Shop E）。
查看预期输出：
>loop.function(df)
apple orange pear
apple 0 0 2
orange 2 0 1
pear 1 2 0
这里可以看到pear在apple之前选择了两次（在ShopC和D），apple在pear之前选择了一次（在ShopA）。
我正在努力提高我对矢量化的了解，尤其是在代替循环方面，所以我想知道如何对这个循环进行矢量化。
（我感觉可能有使用outer()的解决方案，但我对向量化函数的了解仍然非常有限。）
更新
查看 times = 10000 的 loop.function()、tidyverse.function()、loop.function2()、datatable.function() 和 loop.function.TMS() 的基准测试：
Unit: milliseconds
expr min lq mean median uq max neval cld
loop.function(dat) 186.588600 202.78350 225.724249 215.56575 234.035750 999.8234 10000 e
tidyverse.function(dat) 21.523400 22.93695 26.795815 23.67290 26.862700 295.7456 10000 c
loop.function2(dat) 119.695400 126.48825 142.568758 135.23555 148.876100 929.0066 10000 d
datatable.function(dat) 8.517600 9.28085 10.644163 9.97835 10.766749 215.3245 10000 b
loop.function.TMS(dat) 4.482001 5.08030 5.916408 5.38215 5.833699 77.1935 10000 a
对我来说最有趣的结果可能是 tidyverse.function() 在真实数据上的表现。我将不得不在以后尝试添加 Rccp 解决方案 - 我无法让它们处理真实数据。
我感谢所有对此 post 的关注和回答 - 我的目的是学习和提高性能，从所有给出的评论和解决方案中肯定可以学到很多东西。谢谢！
好的，这是一个解决方案：
library(tidyverse)
# a dataframe with all fruit combinations
df_compare <- expand.grid(row_fruit = unique(df$Fruit)
, column_fruit = unique(df$Fruit)
, stringsAsFactors = FALSE)
df_compare %>%
left_join(df, by = c("row_fruit" = "Fruit")) %>%
left_join(df, by = c("column_fruit" = "Fruit")) %>%
filter(Shop.x == Shop.y &
Order.x < Order.y) %>%
group_by(row_fruit, column_fruit) %>%
summarise(obs = n()) %>%
pivot_wider(names_from = row_fruit, values_from = obs) %>%
arrange(column_fruit) %>%
mutate_if(is.numeric, function(x) replace_na(x, 0)) %>%
column_to_rownames("column_fruit") %>%
as.matrix()
apple orange pear
apple 0 0 2
orange 2 0 1
pear 1 2 0
如果您不知道第二个代码部分 (df_compare %>% ...) 中发生了什么，请将“管道”(%>%) 阅读为 'then'。 运行 从 df_compare 到任何管道之前的代码，以查看中间结果。
一个data.table解决方案：
library(data.table)
setDT(df)
setkey(df,Shop)
dcast(df[df,on=.(Shop=Shop),allow.cartesian=T][
,.(cnt=sum(i.Order<Order&i.Fruit!=Fruit)),by=.(Fruit,i.Fruit)]
,Fruit~i.Fruit,value.var='cnt')
Fruit apple orange pear
1: apple 0 0 2
2: orange 2 0 1
3: pear 1 2 0
此示例不需要 Shop 索引，但可能会提高更大数据集的性能。
由于这个问题引起了很多关于性能的评论，我决定检查一下 Rcpp 能带来什么：
library(Rcpp)
cppFunction('NumericMatrix rcppPair(DataFrame df) {
std::vector<std::string> Shop = Rcpp::as<std::vector<std::string> >(df["Shop"]);
Rcpp::NumericVector Order = df["Order"];
Rcpp::StringVector Fruit = df["Fruit"];
StringVector FruitLevels = sort_unique(Fruit);
IntegerVector FruitInt = match(Fruit, FruitLevels);
int n = FruitLevels.length();
std::string currentShop = "";
int order, fruit, i, f;
NumericMatrix result(n,n);
NumericVector fruitOrder(n);
for (i=0;i<Fruit.length();i++){
if (currentShop != Shop[i]) {
//Init counter for each shop
currentShop = Shop[i];
std::fill(fruitOrder.begin(), fruitOrder.end(), 0);
}
order = Order[i];
fruit = FruitInt[i];
fruitOrder[fruit-1] = order;
for (f=0;f<n;f++) {
if (order > fruitOrder[f] & fruitOrder[f]>0 ) {
result(fruit-1,f) = result(fruit-1,f)+1;
}
}
}
rownames(result) = FruitLevels;
colnames(result) = FruitLevels;
return(result);
}
')
rcppPair(df)
apple orange pear
apple 0 0 2
orange 2 0 1
pear 1 2 0
在示例数据集上，这比 data.table 解决方案运行 >500 倍，可能是因为它没有笛卡尔积问题。这在错误输入时不应该是健壮的，并且期望商店/订单按升序排列。
考虑到 data.table 解决方案的 3 行代码花费了几分钟时间，与更长的 Rcpp 解决方案/调试过程相比，我不建议选择 Rcpp 除非存在真正的性能瓶颈。
但有趣的是要记住，如果性能是必须的，Rcpp 可能值得付出努力。
这里有一种方法可以进行简单的修改，使其速度提高 5 倍。
loop.function2 <- function(df){
spl_df = split(df[, c(1L, 3L)], df[[2L]])
mat <- array(0L,
dim=c(length(spl_df), length(spl_df)),
dimnames = list(names(spl_df), names(spl_df)))
for (m in 1:(length(spl_df) - 1L)) {
xm = spl_df[[m]]
mShop = xm$Shop
for (n in ((1+m):length(spl_df))) {
xn = spl_df[[n]]
mm = match(mShop, xn$Shop)
inds = which(!is.na(mm))
mOrder = xm[inds, "Order"]
nOrder = xn[mm[inds], "Order"]
mat[m, n] <- sum(nOrder < mOrder)
mat[n, m] <- sum(mOrder < nOrder)
}
}
mat
}
有 3 个主要概念：
- 原来的
df[df$Fruits == fruits[m], ] 行效率低下，因为您要进行 length(Fruits)^2 次相同的比较。相反，我们可以使用 split()，这意味着我们只扫描水果一次。
df$var 有很多用途，它会在每个循环中提取向量。在这里，我们将 xm 的赋值放在内部循环之外，并尽量减少我们需要子集/提取的内容。- 我将其更改为更接近
combn，因为我们可以通过同时执行 sum(xmOrder > xnOrder) 然后将其切换为 sum(xmOrder < xnOrder) 来重新使用 match() 条件。
性能：
bench::mark(loop.function(df), loop.function2(df))
# A tibble: 2 x 13
## expression min median
## <bch:expr> <bch:tm> <bch:>
##1 loop.function(df) 3.57ms 4.34ms
##2 loop.function2(df) 677.2us 858.6us
我的直觉是，对于更大的数据集，@Waldi 的 data.table 解决方案会更快。但是对于较小的数据集，这应该是相当高效的。
最后，这是另一种 rcpp 方法，它似乎比@Waldi 慢：
#include <Rcpp.h>
using namespace Rcpp;
// [[Rcpp::export]]
IntegerMatrix loop_function_cpp(List x) {
int x_size = x.size();
IntegerMatrix ans(x_size, x_size);
for (int m = 0; m < x_size - 1; m++) {
DataFrame xm = x[m];
CharacterVector mShop = xm[0];
IntegerVector mOrder = xm[1];
int nrows = mShop.size();
for (int n = m + 1; n < x_size; n++) {
DataFrame xn = x[n];
CharacterVector nShop = xn[0];
IntegerVector nOrder = xn[1];
for (int i = 0; i < nrows; i++) {
for (int j = 0; j < nrows; j++) {
if (mShop[i] == nShop[j]) {
if (mOrder[i] > nOrder[j])
ans(m, n)++;
else
ans(n, m)++;
break;
}
}
}
}
}
return(ans);
}
loop_wrapper = function(df) {
loop_function_cpp(split(df[, c(1L, 3L)], df[[2L]]))
}
loop_wrapper(df)
``
似乎无法对原始数据帧进行矢量化 df。但是，如果您使用 reshape2::dcast() 对其进行转换，则每个商店一行：
require(reshape2)
df$Fruit <- as.character(df$Fruit)
by_shop <- dcast(df, Shop ~ Fruit, value.var = "Order")
# Shop apple orange pear
# 1 A 1 2 3
# 2 B NA 1 2
# 3 C 2 NA 1
# 4 D 2 3 1
# 5 E 1 1 1
...，那么您至少可以轻松地对 [m, n] 的每个组合进行矢量化：
fruits <- unique(df$Fruit)
outer(fruits, fruits,
Vectorize(
function (m, n, by_shop) sum(by_shop[,m] > by_shop[,n], na.rm = TRUE),
c("m", "n")
),
by_shop)
# [,1] [,2] [,3]
# [1,] 0 0 2
# [2,] 2 0 1
# [3,] 1 2 0
这可能是您希望对 outer 执行的解决方案。更快的解决方案是对水果 [m，n] 的所有组合进行真正的矢量化，但我一直在考虑它，但我没有看到任何方法可以做到这一点。所以我不得不使用 Vectorize 函数，这当然比真正的矢量化慢得多。
与您的原始函数的基准比较：
Unit: milliseconds
expr min lq mean median uq max neval
loop.function(df) 3.788794 3.926851 4.157606 4.002502 4.090898 9.529923 100
loop.function.TMS(df) 1.582858 1.625566 1.804140 1.670095 1.756671 8.569813 100
函数和基准代码（还添加了 dimnames 的保存）：
require(reshape2)
loop.function.TMS <- function(df) {
df$Fruit <- as.character(df$Fruit)
by_shop <- dcast(df, Shop ~ Fruit, value.var = "Order")
fruits <- unique(df$Fruit)
o <- outer(fruits, fruits, Vectorize(function (m, n, by_shop) sum(by_shop[,m] > by_shop[,n], na.rm = TRUE), c("m", "n")), by_shop)
colnames(o) <- rownames(o) <- fruits
o
}
require(microbenchmark)
microbenchmark(loop.function(df), loop.function.TMS(df))
我想加速创建一个成对矩阵的函数，该矩阵描述了在一组位置中在所有其他对象之前和之后选择一个对象的次数。
这是一个例子df:
df <- data.frame(Shop = c("A","A","A","B","B","C","C","D","D","D","E","E","E"),
Fruit = c("apple", "orange", "pear",
"orange", "pear",
"pear", "apple",
"pear", "apple", "orange",
"pear", "apple", "orange"),
Order = c(1, 2, 3,
1, 2,
1, 2,
1, 2, 3,
1, 1, 1))
在每个 Shop 中，Fruit 是由客户在给定的 Order 中挑选的。
以下函数创建一个 m x n 成对矩阵：
loop.function <- function(df){
fruits <- unique(df$Fruit)
nt <- length(fruits)
mat <- array(dim=c(nt,nt))
for(m in 1:nt){
for(n in 1:nt){
## filter df for each pair of fruit
xm <- df[df$Fruit == fruits[m],]
xn <- df[df$Fruit == fruits[n],]
## index instances when a pair of fruit are picked in same shop
mm <- match(xm$Shop, xn$Shop)
## filter xm and xn based on mm
xm <- xm[! is.na(mm),]
xn <- xn[mm[! is.na(mm)],]
## assign number of times fruit[m] is picked after fruit[n] to mat[m,n]
mat[m,n] <- sum(xn$Order < xm$Order)
}
}
row.names(mat) <- fruits
colnames(mat) <- fruits
return(mat)
}
其中 mat[m,n] 是 fruits[m] 被选中的次数 在 fruits[n] 之后。而 mat[n,m] 是 fruits[m] 在 fruits[n] 之前 被选中的次数。如果同时采摘成对的水果则不记录（例如 Shop E）。
查看预期输出：
>loop.function(df)
apple orange pear
apple 0 0 2
orange 2 0 1
pear 1 2 0
这里可以看到pear在apple之前选择了两次（在ShopC和D），apple在pear之前选择了一次（在ShopA）。
我正在努力提高我对矢量化的了解，尤其是在代替循环方面，所以我想知道如何对这个循环进行矢量化。
（我感觉可能有使用outer()的解决方案，但我对向量化函数的了解仍然非常有限。）
更新
查看 times = 10000 的 loop.function()、tidyverse.function()、loop.function2()、datatable.function() 和 loop.function.TMS() 的基准测试：
Unit: milliseconds
expr min lq mean median uq max neval cld
loop.function(dat) 186.588600 202.78350 225.724249 215.56575 234.035750 999.8234 10000 e
tidyverse.function(dat) 21.523400 22.93695 26.795815 23.67290 26.862700 295.7456 10000 c
loop.function2(dat) 119.695400 126.48825 142.568758 135.23555 148.876100 929.0066 10000 d
datatable.function(dat) 8.517600 9.28085 10.644163 9.97835 10.766749 215.3245 10000 b
loop.function.TMS(dat) 4.482001 5.08030 5.916408 5.38215 5.833699 77.1935 10000 a
对我来说最有趣的结果可能是 tidyverse.function() 在真实数据上的表现。我将不得不在以后尝试添加 Rccp 解决方案 - 我无法让它们处理真实数据。
我感谢所有对此 post 的关注和回答 - 我的目的是学习和提高性能，从所有给出的评论和解决方案中肯定可以学到很多东西。谢谢！
好的，这是一个解决方案：
library(tidyverse)
# a dataframe with all fruit combinations
df_compare <- expand.grid(row_fruit = unique(df$Fruit)
, column_fruit = unique(df$Fruit)
, stringsAsFactors = FALSE)
df_compare %>%
left_join(df, by = c("row_fruit" = "Fruit")) %>%
left_join(df, by = c("column_fruit" = "Fruit")) %>%
filter(Shop.x == Shop.y &
Order.x < Order.y) %>%
group_by(row_fruit, column_fruit) %>%
summarise(obs = n()) %>%
pivot_wider(names_from = row_fruit, values_from = obs) %>%
arrange(column_fruit) %>%
mutate_if(is.numeric, function(x) replace_na(x, 0)) %>%
column_to_rownames("column_fruit") %>%
as.matrix()
apple orange pear
apple 0 0 2
orange 2 0 1
pear 1 2 0
如果您不知道第二个代码部分 (df_compare %>% ...) 中发生了什么，请将“管道”(%>%) 阅读为 'then'。 运行 从 df_compare 到任何管道之前的代码，以查看中间结果。
一个data.table解决方案：
library(data.table)
setDT(df)
setkey(df,Shop)
dcast(df[df,on=.(Shop=Shop),allow.cartesian=T][
,.(cnt=sum(i.Order<Order&i.Fruit!=Fruit)),by=.(Fruit,i.Fruit)]
,Fruit~i.Fruit,value.var='cnt')
Fruit apple orange pear
1: apple 0 0 2
2: orange 2 0 1
3: pear 1 2 0
此示例不需要 Shop 索引，但可能会提高更大数据集的性能。
由于这个问题引起了很多关于性能的评论，我决定检查一下 Rcpp 能带来什么：
library(Rcpp)
cppFunction('NumericMatrix rcppPair(DataFrame df) {
std::vector<std::string> Shop = Rcpp::as<std::vector<std::string> >(df["Shop"]);
Rcpp::NumericVector Order = df["Order"];
Rcpp::StringVector Fruit = df["Fruit"];
StringVector FruitLevels = sort_unique(Fruit);
IntegerVector FruitInt = match(Fruit, FruitLevels);
int n = FruitLevels.length();
std::string currentShop = "";
int order, fruit, i, f;
NumericMatrix result(n,n);
NumericVector fruitOrder(n);
for (i=0;i<Fruit.length();i++){
if (currentShop != Shop[i]) {
//Init counter for each shop
currentShop = Shop[i];
std::fill(fruitOrder.begin(), fruitOrder.end(), 0);
}
order = Order[i];
fruit = FruitInt[i];
fruitOrder[fruit-1] = order;
for (f=0;f<n;f++) {
if (order > fruitOrder[f] & fruitOrder[f]>0 ) {
result(fruit-1,f) = result(fruit-1,f)+1;
}
}
}
rownames(result) = FruitLevels;
colnames(result) = FruitLevels;
return(result);
}
')
rcppPair(df)
apple orange pear
apple 0 0 2
orange 2 0 1
pear 1 2 0
在示例数据集上，这比 data.table 解决方案运行 >500 倍，可能是因为它没有笛卡尔积问题。这在错误输入时不应该是健壮的，并且期望商店/订单按升序排列。
考虑到 data.table 解决方案的 3 行代码花费了几分钟时间，与更长的 Rcpp 解决方案/调试过程相比，我不建议选择 Rcpp 除非存在真正的性能瓶颈。
但有趣的是要记住，如果性能是必须的，Rcpp 可能值得付出努力。
这里有一种方法可以进行简单的修改，使其速度提高 5 倍。
loop.function2 <- function(df){
spl_df = split(df[, c(1L, 3L)], df[[2L]])
mat <- array(0L,
dim=c(length(spl_df), length(spl_df)),
dimnames = list(names(spl_df), names(spl_df)))
for (m in 1:(length(spl_df) - 1L)) {
xm = spl_df[[m]]
mShop = xm$Shop
for (n in ((1+m):length(spl_df))) {
xn = spl_df[[n]]
mm = match(mShop, xn$Shop)
inds = which(!is.na(mm))
mOrder = xm[inds, "Order"]
nOrder = xn[mm[inds], "Order"]
mat[m, n] <- sum(nOrder < mOrder)
mat[n, m] <- sum(mOrder < nOrder)
}
}
mat
}
有 3 个主要概念：
- 原来的
df[df$Fruits == fruits[m], ]行效率低下，因为您要进行length(Fruits)^2次相同的比较。相反，我们可以使用split()，这意味着我们只扫描水果一次。 df$var有很多用途，它会在每个循环中提取向量。在这里，我们将xm的赋值放在内部循环之外，并尽量减少我们需要子集/提取的内容。- 我将其更改为更接近
combn，因为我们可以通过同时执行sum(xmOrder > xnOrder)然后将其切换为sum(xmOrder < xnOrder)来重新使用match()条件。
性能：
bench::mark(loop.function(df), loop.function2(df))
# A tibble: 2 x 13
## expression min median
## <bch:expr> <bch:tm> <bch:>
##1 loop.function(df) 3.57ms 4.34ms
##2 loop.function2(df) 677.2us 858.6us
我的直觉是，对于更大的数据集，@Waldi 的 data.table 解决方案会更快。但是对于较小的数据集，这应该是相当高效的。
最后，这是另一种 rcpp 方法，它似乎比@Waldi 慢：
#include <Rcpp.h>
using namespace Rcpp;
// [[Rcpp::export]]
IntegerMatrix loop_function_cpp(List x) {
int x_size = x.size();
IntegerMatrix ans(x_size, x_size);
for (int m = 0; m < x_size - 1; m++) {
DataFrame xm = x[m];
CharacterVector mShop = xm[0];
IntegerVector mOrder = xm[1];
int nrows = mShop.size();
for (int n = m + 1; n < x_size; n++) {
DataFrame xn = x[n];
CharacterVector nShop = xn[0];
IntegerVector nOrder = xn[1];
for (int i = 0; i < nrows; i++) {
for (int j = 0; j < nrows; j++) {
if (mShop[i] == nShop[j]) {
if (mOrder[i] > nOrder[j])
ans(m, n)++;
else
ans(n, m)++;
break;
}
}
}
}
}
return(ans);
}
loop_wrapper = function(df) {
loop_function_cpp(split(df[, c(1L, 3L)], df[[2L]]))
}
loop_wrapper(df)
``
似乎无法对原始数据帧进行矢量化 df。但是，如果您使用 reshape2::dcast() 对其进行转换，则每个商店一行：
require(reshape2)
df$Fruit <- as.character(df$Fruit)
by_shop <- dcast(df, Shop ~ Fruit, value.var = "Order")
# Shop apple orange pear
# 1 A 1 2 3
# 2 B NA 1 2
# 3 C 2 NA 1
# 4 D 2 3 1
# 5 E 1 1 1
...，那么您至少可以轻松地对 [m, n] 的每个组合进行矢量化：
fruits <- unique(df$Fruit)
outer(fruits, fruits,
Vectorize(
function (m, n, by_shop) sum(by_shop[,m] > by_shop[,n], na.rm = TRUE),
c("m", "n")
),
by_shop)
# [,1] [,2] [,3]
# [1,] 0 0 2
# [2,] 2 0 1
# [3,] 1 2 0
这可能是您希望对 outer 执行的解决方案。更快的解决方案是对水果 [m，n] 的所有组合进行真正的矢量化，但我一直在考虑它，但我没有看到任何方法可以做到这一点。所以我不得不使用 Vectorize 函数，这当然比真正的矢量化慢得多。
与您的原始函数的基准比较：
Unit: milliseconds
expr min lq mean median uq max neval
loop.function(df) 3.788794 3.926851 4.157606 4.002502 4.090898 9.529923 100
loop.function.TMS(df) 1.582858 1.625566 1.804140 1.670095 1.756671 8.569813 100
函数和基准代码（还添加了 dimnames 的保存）：
require(reshape2)
loop.function.TMS <- function(df) {
df$Fruit <- as.character(df$Fruit)
by_shop <- dcast(df, Shop ~ Fruit, value.var = "Order")
fruits <- unique(df$Fruit)
o <- outer(fruits, fruits, Vectorize(function (m, n, by_shop) sum(by_shop[,m] > by_shop[,n], na.rm = TRUE), c("m", "n")), by_shop)
colnames(o) <- rownames(o) <- fruits
o
}
require(microbenchmark)
microbenchmark(loop.function(df), loop.function.TMS(df))