R 提高功能的性能

Question

这个问题与我的previous one 有关。这是一个小样本数据。我同时使用了data.table 和data.frame 来找到更快的解决方案。

test.dt <- data.table(strt=c(1,1,2,3,5,2), end=c(2,1,5,5,5,4), a1.2=c(1,2,3,4,5,6), 
                   a2.3=c(2,4,6,8,10,12), a3.4=c(3,1,2,4,5,1), a4.5=c(5,1,15,10,12,10),
                   a5.6=c(4,8,2,1,3,9))

test.dt[,rown:=as.numeric(row.names(test.dt))]

test.df <- data.frame(strt=c(1,1,2,3,5,2), end=c(2,1,5,5,5,4), a1.2=c(1,2,3,4,5,6), 
                   a2.3=c(2,4,6,8,10,12), a3.4=c(3,1,2,4,5,1), a4.5=c(5,1,15,10,12,10),
                   a5.6=c(4,8,2,1,3,9))

test.df$rown <- as.numeric(row.names(test.df))

    > test.df
  strt end a1.2 a2.3 a3.4 a4.5 a5.6 rown
1    1   2    1    2    3    5    4    1
2    1   1    2    4    1    1    8    2
3    2   5    3    6    2   15    2    3
4    3   5    4    8    4   10    1    4
5    5   5    5   10    5   12    3    5
6    2   4    6   12    1   10    9    6

我想使用开始和结束列值来确定子集的列范围（从 a1.2 到 a5.6 的列）并获得平均值.比如第一行，由于strt=1和end=2，我需要得到a1.2和a2.3的均值；在第三行，我需要得到 a2.3、a3.4、a4.5 和 a5 的平均值.6

输出应该是这样的向量

> k
       1        2        3        4        5        6 
1.500000 2.000000 6.250000 5.000000 3.000000 7.666667 

这是我尝试过的：

解决方案 1：这使用 data.table 并在其上应用函数。

func.dt <- function(rown, x, y) {
      tmp  <- paste0("a", x, "." , x+1)
      tmp1 <- paste0("a", y, "." , y+1)
      rowMeans(test.dt[rown,get(tmp):get(tmp1), with=FALSE])
      }
    k <- test.dt[, func.dt(rown, strt, end), by=.(rown)]

解决方案 2：这使用 data.frame 并在其上应用一个函数。

func.df <- function(rown, x, y) {
  rowMeans(test.df[rown,(x+2):(y+2), drop=FALSE])
  }
k1 <- mapply(func.df, test.df$rown, test.df$strt, test.df$end)

解决方案 3：这使用 data.frame 并循环遍历它。

    test.ave <- rep(NA, length(test1$strt))
for (i in 1 : length(test.df$strt)) {
    test.ave[i] <- rowMeans(test.df[i, as.numeric(test.df[i,1]+2):as.numeric(test.df[i,2]+2), drop=FALSE])
    }

基准测试表明解决方案 2 是最快的。

test replications elapsed relative user.self sys.self user.child sys.child
1 sol1          100    0.67    4.786      0.67        0         NA        NA
2 sol2          100    0.14    1.000      0.14        0         NA        NA
3 sol3          100    0.15    1.071      0.16        0         NA        NA

但是，这对我来说还不够好。考虑到我的数据大小，这些函数需要运行几天才能得到输出。我确信我没有充分利用data.table 的力量，而且我也知道我的函数很糟糕（它们引用全局环境中的数据集而不传递它）。不幸的是，我的深度不够，不知道如何解决这些问题并使我的功能快速。我将非常感谢任何有助于改进我的功能或指出替代解决方案的建议。

Answer 1

我的解决方案是基准测试中的第一个

library(microbenchmark)
microbenchmark(
  lapply(
    apply(test.df,1, function(x){
        x[(x[1]+2):(x[2]+2)]}),
    mean),
  test.dt[, func.dt(rown, strt, end), by=.(rown)]
)

      min        lq      mean   median       uq      max neval
  138.654  175.7355  254.6245  201.074  244.810 3702.443   100
 4243.641 4747.5195 5576.3399 5252.567 6247.201 8520.286   100

它似乎快了 25 倍，但这是一个小数据集。我确信有比我所做的更好的方法来做到这一点。

Answer 2

除非你能想出一个聪明的子集方法来做到这一点，否则我认为你已经达到了 R 的速度障碍。对于这个问题，您需要使用像 C++ 这样的低级语言。幸运的是，Rcpp 包使得与R 中的C++ 的接口变得简单。免责声明：我一生中从未编写过一行 C++ 代码。这段代码可能效率很低。

library(Rcpp)

cppFunction('NumericVector MYrcpp(NumericMatrix x) {
  int nrow = x.nrow(), ncol = x.ncol();
  NumericVector out(nrow);

  for (int i = 0; i < nrow; i++) {
    double avg = 0;

    int start = x(i,0);
    int end = x(i,1);

    int N = end - start + 1;

    while(start<=end){

      avg += x(i, start + 1); 
    
      start = start + 1;

    }

    out[i] = avg/N;

  }
  return out;
}')

对于这段代码，我将把data.frame 作为matrix 传递（即testM <- as.matrix(test.df)）

让我们看看它是否有效......

 MYrcpp(testM)
[1] 1.500000 2.000000 6.250000 5.000000 3.000000 7.666667

有多快？

Unit: microseconds
          expr      min        lq      mean   median       uq       max neval
          f2() 1543.099 1632.3025 2039.7350 1843.458 2246.951  4735.851   100
          f3() 1859.832 1993.0265 2642.8874 2168.012 2493.788 19619.882   100
          f4()  281.541  315.2680  364.2197  345.328  375.877  1089.994   100
 MYrcpp(testM)    3.422   10.0205   16.7708   19.552   21.507    56.700   100

其中f2()、f3() 和f4() 定义为

f2 <- function(){
  func.df <- function(rown, x, y) {
    rowMeans(test.df[rown,(x+2):(y+2), drop=FALSE])
  }
  k1 <- mapply(func.df, test.df$rown, test.df$strt, test.df$end)
}

f3 <- function(){
  test.ave <- rep(NA, length(test.df$strt))
  for (i in 1 : length(test.df$strt)) {
    test.ave[i] <- rowMeans(test.df[i,as.numeric(test.df[i,1]+2):as.numeric(test.df[i,2]+2), drop=FALSE])
  }
}

f4 <- function(){
  lapply(
    apply(test.df,1, function(x){
      x[(x[1]+2):(x[2]+2)]}),
    mean)
}

这大约是最快的 20 倍。

注意，要实现上述代码，您需要R 可以访问的C 编译器。对于 Windows，请查看 Rtools。有关Rcpp 的更多信息，请阅读this

现在让我们看看它是如何扩展的。

N = 5e3
test.df <- data.frame(strt = 1, 
                 end = sample(5, N, replace = TRUE), 
                 a1.2 = sample(3, N, replace = TRUE), 
                 a2.3 = sample(7, N, replace = TRUE), 
                 a3.4 = sample(14, N, replace = TRUE),
                 a4.5 = sample(8, N, replace = TRUE),
                 a5.6 = sample(30, N, replace = TRUE))
test.df$rown <- as.numeric(row.names(test.df))


test.dt <- as.data.table(test.df)

microbenchmark(f4(), MYrcpp(testM))
Unit: microseconds
          expr       min         lq        mean     median          uq       max neval
          f4() 88647.256 108314.549 125451.4045 120736.073 133487.5295 259502.49   100
 MYrcpp(testM)   196.003    216.533    242.6732    235.107    261.0125    499.54   100

使用 5e3 行 MYrcpp 现在快 550 倍。这部分是因为f4() 不会像 Richard 在评论中讨论的那样很好地扩展。 f4() 本质上是通过在 lapply 中调用 apply 来调用嵌套的 for 循环。有趣的是，C++ 代码还通过在 for 循环中利用 while 循环来调用嵌套循环。速度差异在很大程度上是由于C++ 代码已经被编译，不需要被打断成机器在运行时可以理解的东西。

我不确定您的数据集有多大，但是当我在 data.frame 和 1e7 行上运行 MYrcpp 时，这是我可以在我破旧的笔记本电脑上分配的最大 data.frame，它运行在 500 毫秒内。

更新：C++ 代码的 R 等效项

MYr <- function(x){
  nrow <- nrow(x)
  ncol <- ncol(x)
  out <- matrix(NA, nrow = 1, ncol = nrow)

  for(i in 1:nrow){
    avg <- 0

    start <- x[i,1]
    end <- x[i,2]

    N <- end - start + 1

    while(start<=end){
  
      avg <- avg + x[i, start + 2]
  
      start = start + 1
  
    }

    out[i] <- avg/N

  }

  out

}

MYrcpp 和 MYr 在很多方面都相似。让我讨论几个不同之处

1. MYrcpp 的第一行与MYr 不同。换句话说，MYrcpp 的第一行 NumericVector MYrcpp(NumericMatrix x) 意味着我们正在定义一个函数，其名称为 MYrcpp，它返回类 NumericVector 的输出并接受类 NumericMatrix 的输入 x。
2. 在C++中，引入变量时必须定义其类，即int nrow = x.row()是一个变量，其名称为nrow，其类为int（即整数）并赋值为x.nrow()即x的行数。 （如果您不知所措，请忽略，nrow() 是类 `NumericVector 实例的方法。就像在 Python 中一样，您通过将方法附加到实例来调用方法。R 等效于 S3 和 S4 方法）
3. 当您在 C++ 中设置子集时，您使用 () 而不是 R 中的 []。此外，索引从零开始（如在 Python 中）。例如，C++ 中的x(0,1) 等价于R 中的x[1,2]
4. ++ 是一个运算符，表示加 1，即 j++ 与 j + 1 相同。 += 是一个运算符，表示加在一起并赋值，即a += b 与a = a + b 相同

Answer 3

我很好奇在不编写自定义 C 或 C++ 代码的情况下我能以多快的速度完成此操作。我能想到的最好的如下。请注意，使用mean.default 将提供更高的精度，因为它会第二次传递数据以进行纠错。

f_jmu <- compiler::cmpfun({function(m) {
  # remove start/end columns from 'm' matrix
  ma <- m[,-(1:2)]
  # column index for each row in 'ma' matrix
  cm <- col(ma)
  # logical index of whether we need the column for each row
  i <- cm >= m[,1L] & cm <= m[,2L]
  # multiply the input matrix by the index matrix and sum it
  # divide by the sum of the index matrix to get the mean
  rowSums(i*ma) / rowSums(i)
}})

Rcpp 函数仍然更快（不足为奇），但上面的函数相当接近。这是我的笔记本电脑上 5000 万次观察的示例，笔记本电脑配备 i7-4600U 和 12GB 内存。

set.seed(21)
N <- 5e7
test.df <- data.frame(strt = 1L, 
                 end = sample(5, N, replace = TRUE), 
                 a1.2 = sample(3, N, replace = TRUE), 
                 a2.3 = sample(7, N, replace = TRUE), 
                 a3.4 = sample(14, N, replace = TRUE),
                 a4.5 = sample(8, N, replace = TRUE),
                 a5.6 = sample(30, N, replace = TRUE))
test.df$strt <- pmax(1L, test.df$end - sample(3, N, replace = TRUE) + 1L)
test.m <- as.matrix(test.df)

还要注意，我注意确保test.m 是一个整数 矩阵。这有助于减少内存占用，从而加快速度。

R> system.time(st1 <- MYrcpp(test.m))
   user  system elapsed 
  0.900   0.216   1.112 
R> system.time(st2 <- f_jmu(test.m))
   user  system elapsed 
  6.804   0.756   7.560 
R> identical(st1, st2)
[1] TRUE