从 Rcpp 中的列表中提取元素似乎有点慢

It seems it is a bit slow to extract element from a List in Rcpp

我刚刚用三个相同大小的输入向量写了一个 Rcpp 函数,x(数字)y(数字)和 category(字符)。然后我要 return 一个列表,列表的大小等于唯一类别值的长度。此列表中的每个元素都是基于 xy 具有相应类别的相同大小矩阵(相等的行和列)。

但是,当 n 的大小很大时,我发现我的代码不够快。我认为原因是我需要从列表中提取一些东西,进行一些计算并每次都将其插入。有没有人对如何加快这个过程有什么建议。

Rcpp代码

#include <Rcpp.h>
using namespace Rcpp;

//[[Rcpp::export]]
List myList(NumericVector x, NumericVector y, CharacterVector category) {

  int n = x.size();
  CharacterVector levels = unique(category);
  int levels_size = levels.size();
  List L(levels_size);

  int plot_width = 600;
  int plot_height = 600;

  // Each element in the list L has the same size Matrix
  for(int j = 0; j < levels_size; j++) {
    NumericMatrix R(plot_height, plot_width);
    L[j] = R;
  }
  int id = 0;

  double xmax = max(x);
  double ymax = max(y);
  double xmin = min(x);
  double ymin = min(y);

  for(int i=0; i < n; i++) {

    for(int j = 0; j < levels_size; j++) {
      if(category[i] == levels[j]) {
        id = j;
        break;
      }
    }

    int id_x = floor((x[i] - xmin)/(xmax - xmin) * (plot_width - 1));
    int id_y = floor((y[i] - ymin)/(ymax - ymin) * (plot_height - 1));

    NumericMatrix M = L[id];
    // some computation in M
    M(id_y, id_x) += 1;
    L[id] = M;
  }
  return(L);
}

R码

n <- 1e8
class <- 20

x <- rnorm(n)
y <- rnorm(n)
category <- sample(as.factor(1:class), size = n, replace = TRUE)

start_time <- Sys.time()
L <- myList(x = x, y = y, category = category)
end_time <- Sys.time()
end_time - start_time
# Time difference of 35.3367 secs

我怀疑有关性能的两个主要问题:

  • 很多字符串比较(1e9 的顺序)
  • 矩阵有很多缓存未命中,因为通常两个连续的 xy 对不会来自同一类别,因此需要不同的矩阵

两者都指向同一个方向:不要尝试实现自己的 GROUP BY 操作。像 data.table 这样的数据库引擎和包更清楚如何做到这一点。例如,当使用 data.table 时,我们需要一个更简单的函数,它期望 x 和 y 用于一个类别 并输出单个矩阵:

#include <Rcpp.h>
using namespace Rcpp;

//[[Rcpp::export]]
NumericMatrix getMat(NumericVector x, NumericVector y,
                     double xmin, double xmax, double ymin, double ymax,
                     int plot_width = 600, int plot_height = 600) {
    int n = x.size();
    NumericMatrix M(plot_height, plot_width);

    for(int i=0; i < n; i++) {
        int id_x = floor((x[i] - xmin)/(xmax - xmin) * (plot_width - 1));
        int id_y = floor((y[i] - ymin)/(ymax - ymin) * (plot_height - 1));
        M(id_y, id_x) += 1;
    }
    return M;
}

/***R
n <- 1e8
class <- 20

library("data.table")
foo <- data.table(x = rnorm(n),
                  y = rnorm(n),
                  category = sample(as.factor(1:class), size = n, replace = TRUE))

xmin <- min(foo$x)
xmax <- max(foo$x)
ymin <- min(foo$y)
ymax <- max(foo$y)

system.time(bar <- foo[,
                       list(baz = list(getMat(x, y, xmin, xmax, ymin, ymax))),
                       by = category])
*/

备注:

  • 在我的系统上,聚合时间不到 6 秒。
  • 如果在聚合之前做一个 setkey(foo, category) 会更快。不过,这实际上改变了行的顺序。小心使用!
  • data.table 语法有点简洁,但是习惯了...
  • 输出的结构不同,但如果需要可以转换。