向量化函数/提高 data.table 中的计算速度

Question

目前我有以下 data.table ：

item      city    dummyvar
A        Austin       1
A        Austin       1
A        Austin      100
B        Austin       2 
B        Austin       2
B        Austin      200
A          NY         1
A          NY         1
A          NY        100
B          NY         2 
B          NY         2
B          NY        200

我有一个名为ImbalancePoints 的用户定义函数，它应用于dummyvar，它返回检测到dummyvar 突然变化的行。我这样做的方式如下：

my.data.table[,
 .(item, city , imb.points = list(unique(try(ImbalancePoints(dummyvar), silent = T))) ),
 by = .(city, item)
]

对于NY 的情况，假设我得到了一个data.table 对象，如下所示：

 item   city   imb.points
   A     NY     3,449

其中imb.points 列是具有嵌套列表作为其元素的列，对于此示例，数字 3 和 449 表示对于 city = NY 和 item = A 的情况发生突然变化的行。然而，我面临的问题是我有大约。 12个不同城市的3000个不同的项目，计算这个需要很长时间。我想知道您是否可以告诉我如何矢量化/加速这个计算，因为我上次尝试这个花了将近 2 个小时但没有完成。

我不知道它是否有任何帮助，但我也附上了ImbalancePoints 函数：

library(pracma)

ImbalancePr <- function(eval.column) {
  n <- length(eval.column)
  imbalance <- rep(0, n)
  b_t = rep(0,n)
  elem_diff <- diff(eval.column)
  for(i in 2:n)
  {
    imbalance[i] <- sign(elem_diff[i-1]) * (elem_diff[i-1] != 0)
    + imbalance[i-1]*(elem_diff[i-1] == 0)
  }
  return(imbalance)
}

ImbalancePoints <- function(eval.column, w0 = 100, bkw_T = 10, bkw_b = 10){
  bv_t <- ImbalancePr(eval.column)
  w0 <- min(min(which(cumsum(bv_t) != 0)), w0)
  Tstar <- w0
  E0t <- Tstar
  repeat{
    Tlast <- sum(Tstar)
    nbt <- min(bkw_b, Tlast-1)
    P <- pracma::movavg(bv_t[1:Tlast], n = nbt, type = "e")
    P <- tail(P,1)
    bv_t_expected <- E0t * abs(P)
    bv_t_cumsum <- abs(cumsum(bv_t[-(1:Tlast)]))
    if(max(bv_t_cumsum) < bv_t_expected){break}else{
      Tnew <- min(which(bv_t_cumsum >= bv_t_expected))
    }
    Tlast <- Tlast + Tnew
    if(Tlast > length(eval.column)[1]){break}else{
      Tstar <- c(Tstar,Tnew)
      if(length(Tstar) <= 2){
        E0t <- mean(Tstar)
      }else{
        nt <- min(bkw_T,length(Tstar)-1)
        E0t <- pracma::movavg(Tstar[1:length(Tstar)], n = nt, type = "e")
        E0t <- tail(E0t,1)
      }
    }
  }
  return(sort(unique(Tstar)))
}

编辑：感谢 Paul 的洞察力，我的问题只是矢量化 ImbalancePoints 函数内的重复循环。但是，我不是很精通编码，也没有看到直接的解决方案。如果有人可以给我一个建议，或者如果您知道辅助功能/库，我将不胜感激。

Answer 1

这篇文章由几个部分组成，解决了不同的问题：

• 矢量化ImbalancePr()
• 分析ImbalancePoints()
• 将 movavg() 与 Rcpp 加速 4 倍

矢量化ImbalancePr()

我相信ImbalancePr()可以替换成

fImbalancePr <- function(x) c(0, sign(diff(x)))

至少，它返回相同的结果，经过基准测试（检查结果）：

library(bench)
library(ggplot2)
bm <- press(
  n = c(10, 100, 1000, 10000),
  {
    x <- rep(0, n)
    set.seed(123)
    x[sample(n, n/5)] <- 100
    print(table(x))
    mark(
      ImbalancePr(x),
      fImbalancePr(x)
    )
  }
  
)

Running with:
      n
1    10
x
  0 100 
  8   2 
2   100
x
  0 100 
 80  20 
3  1000
x
  0 100 
800 200 
4 10000
x
   0  100 
8000 2000

autoplot(bm)

fImbalancePr() 总是比 OP 的原始版本快。速度优势随着向量长度的增加而增加。

分析ImbalancePoints()

不过，这种改进对ImbalancePoints()的整体性能影响不大：

library(bench)
library(ggplot2)
bm <- press(
  n = c(10L, 100L, 1000L),
  {
    x <- replace(rep(0, n), n, 100)
    y <- c(rep(2, n), rep(-3, n), rep(5, n))
    mark(
      original = {
        list(
          ImbalancePoints(x),
          ImbalancePoints(y)
        )
      },
      modified = {
        list(
          fImbalancePoints(x),
          fImbalancePoints(y)
        )
      }
    )
  }
  
)

fImbalancePoint() 是ImbalancePoint() 的变体，其中ImbalancePr() 已替换为fImbalancePr()。

autoplot(bm)

有一个小的改进，但这无助于显着减少报告的 2 小时执行时间。

我们可以使用profvis来识别ImbalancePoints()内的时间花在哪里：

library(profvis)
x <- c(rep(0, 480L), rep(c(0:9, 9:0), 2L), rep(0, 480L))
profvis({
  for (i in 1:100) ffImbalancePoints(x)
})

时间是通过采样收集的，因此需要足够的重复次数才能获得良好的覆盖率。

RStudio 的屏幕截图显示了一次运行的结果：

• movavg() 消耗了 ImbalancePoints() 所花费时间的 25%。
• 根据分析，pracma::movavg() 中的双冒号运算符需要另外 20%。预先使用 library(pracma) 加载 pracma 包是否可以加快速度可能值得测试。
• 10% 用于调用tail()。 tail(x, 1) 可以替换为 x[length(x)]，速度快了一个数量级以上。

如果我们通过键入pracma::movavg（不带括号）来查看movavg() 的代码，我们会看到有一个无法向量化的迭代循环：

...
else if (type == "e") {
    a <- 2/(n + 1)
    y[1] <- x[1]
    for (k in 2:nx) y[k] <- a * x[k] + (1 - a) * y[k - 1]
}
...

此外，仅使用调用movavg() 创建的时间序列的最后一个值。因此，这里可能有两种性能改进选项：

• 选择其他加权均值函数，该函数仅使用有限窗口内的数据点。
• 使用 Rcpp 在 C++ 中重新实现 movavg()。

用Rcpp 加速movavg()

用Rcpp 函数替换对pracma::movavg() 的调用以及对tail() 的后续调用，我们可以将ImbalancePoints() 的整体速度提高到4 倍。

EMA_last_cpp(x, n) 替换 tail(pracma::movavg(x, n, type = "e"), 1)

library(Rcpp)
cppFunction("
double EMA_last_cpp(const NumericVector& x, const int n) {
  int nx = x.size(); 
  double a = 2.0 / (n + 1.0);
  double b = 1.0 - a;
  double y;

  y = x[0];
  for(int k = 1; k < nx; k++){
    y = a * x[k] + b * y;
  }
  
  return y;
}
")

现在，我们可以相应地修改ImbalancePoints()。另外，替换了ImbalancePr()的调用，另外两个地方修改了代码（见cmets）：

fImbalancePoints <-
  function(eval.column, 
           w0 = 100,
           bkw_T = 10,
           bkw_b = 10) {
    # bv_t <- ImbalancePr(eval.column)
    bv_t <- c(0, sign(diff(eval.column)))
    # w0 <- min(min(which(cumsum(bv_t) != 0)), w0)
    w0 <- min(which(bv_t != 0)[1L], w0) # pick first change point
    Tstar <- w0
    E0t <- Tstar
    repeat {
      Tlast <- sum(Tstar)
      # remove warning: 
      # In max(bv_t_cumsum) : no non-missing arguments to max; returning -Inf
      if (Tlast >= length(bv_t)) break
      nbt <- min(bkw_b, Tlast - 1)
      # P <- movavg(bv_t[1:Tlast], n = nbt, type = "e")
      # P <- tail(P, 1)
      P <- EMA_last_cpp(bv_t[1:Tlast], n = nbt)
      bv_t_expected <- E0t * abs(P)
      bv_t_cumsum <- abs(cumsum(bv_t[-(1:Tlast)]))
      if (max(bv_t_cumsum) < bv_t_expected) {
        break
      } else{
        Tnew <- min(which(bv_t_cumsum >= bv_t_expected))
      }
      Tlast <- Tlast + Tnew
      if (Tlast > length(eval.column)[1]) {
        break
      } else{
        Tstar <- c(Tstar, Tnew)
        if (length(Tstar) <= 2) {
          E0t <- mean(Tstar)
        } else{
          nt <- min(bkw_T, length(Tstar) - 1)
          # E0t <- movavg(Tstar[1:length(Tstar)], n = nt, type = "e")
          # E0t <- tail(E0t, 1)
          E0t <- EMA_last_cpp(Tstar[1:length(Tstar)], n = nt)
        }
      }
    }
    return(sort(unique(Tstar)))
  }

基准

library(bench)
library(ggplot2)
bm <- press(
  n = c(10L, 100L, 1000L),
  {
    x <- replace(rep(0, n), n, 100)
    y <- c(rep(2, n), rep(-3, n), rep(5, n))
    mark(
      original = {
        list(
          ImbalancePoints(x),
          ImbalancePoints(y)
        )
      },
      modified = {
        list(
          fImbalancePoints(x),
          fImbalancePoints(y)
        )
      },
      min_time = 1
    )
  }
)
bm

  expression     n     min median `itr/sec` mem_alloc `gc/sec` n_itr  n_gc total_time
  <bch:expr> <int> <bch:t> <bch:>     <dbl> <bch:byt>    <dbl> <int> <dbl>   <bch:tm>
1 original      10 315.1us  369us   2318.      2.66KB     4.16  2231     4   962.49ms
2 modified      10   120us  136us   6092.    195.11KB     5.31  5733     5   940.99ms
3 original     100 583.4us  671us   1283.     55.09KB     4.16  1234     4   961.78ms
4 modified     100 145.4us  167us   5146.     47.68KB     4.19  4916     4   955.29ms
5 original    1000 438.1ms  469ms      2.17  157.37MB     4.33     3     6      1.38s
6 modified    1000  97.1ms  103ms      9.53  152.09MB    17.1     10    18      1.05s

显示修改后的版本比原始版本快大约 3 到 5 倍。这可能有助于 OP 将他的生产数据集的计算时间从 2 个多小时减少一个重要因素。