R：矢量化循环以创建成对矩阵

Question

我想加快创建成对矩阵的函数，该矩阵描述了在一组位置中，在所有其他对象之前和之后选择对象的次数。

这是一个例子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] 被选中的次数。如果同时采摘成对水果（例如在ShopE），则不会记录。

查看预期输出：

>loop.function(df)
       apple orange pear
apple      0      0    2
orange     2      0    1
pear       1      2    0

您可以在这里看到pear 在apple 之前被选择了两次（在Shop C 和D 中），apple 在pear 之前被选择了一次（在Shop @987654346 @)。

我正在努力提高我对矢量化的了解，尤其是在循环的地方，所以我想知道如何对这个循环进行矢量化。

（我感觉可能有使用outer()的解决方案，但我对矢量化函数的了解仍然非常有限。）

更新

查看times = 10000loop.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 解决方案 - 我无法让它们处理真实数据。

感谢大家对这篇文章的兴趣和回答 - 我的目的是学习和提高性能，从给出的所有 cmets 和解决方案中肯定有很多东西可以学习。谢谢！

Answer 1

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 索引对于此示例不是必需的，但可能会提高更大数据集的性能。

由于这个问题引起了许多 cmets 的性能问题，我决定检查 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 可能值得付出努力。

Answer 2

这是一种通过简单修改使其速度提高 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个概念：

1. 原始的df[df$Fruits == fruits[m], ] 行效率低下，因为您将进行相同的比较length(Fruits)^2 次。相反，我们可以使用split()，这意味着我们只扫描水果一次。
2. df$var 被大量使用，它将在每个循环期间提取向量。在这里，我们将xm 的赋值放在内部循环之外，并尝试最小化我们需要子集/提取的内容。
3. 我将其更改为更接近 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)
``

Answer 3

似乎无法对原始数据框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))

Answer 4

好的，这里有一个解决方案：

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 到任何管道之前的代码以查看中间结果。