rowsum 的 Rcpp 等效项 [关闭]

Posted 2023-02-17

【中文标题】rowsum 的 Rcpp 等效项 [关闭]

【英文标题】：Rcpp equivalent for rowsum [closed]

【发布时间】：2013-06-07 01:49:53

【问题描述】：

我正在为 C++ / Rcpp / Eigen 或 Armadillo 中的 R 函数 rowsum 寻找一个快速的替代方案。

目的是根据分组向量b得到向量a中元素的总和。例如：

> a
 [1] 2 2 2 2 2 2 2 2 2 2    
> b
 [1] 1 1 1 1 1 2 2 2 2 2
> rowsum(a,b)
  [,1]
1   10
2   10

在Rcpp 中编写一个简单的 for 循环非常慢，但也许我的代码效率低下。

我也尝试在Rcpp 中调用函数rowsum，但是，rowsum 不是很快。

【问题讨论】：

代码不使用提供的数据。当向量设计为与矩阵一起使用时，您在向量上使用rowsum。您没有提供任何 Cpp 代码。

rowsum 在上述情况下调度 rowsum.default 并且已经调用了 C 代码，因此它应该已经相当快了。您可以通过直接调用rowsum.default 或.Internal(rowsum_matrix(...)) 来获得一点点速度提升，尽管不鼓励使用后者并且不允许在 CRAN 上使用。

你查过这里的犰狳手册吗：arma.sourceforge.net/docs.html#sum至少有一些求和函数，符合你的目的吗？

听起来像是 data.table 擅长...

【参考方案1】：

不是答案，但可能有助于解决问题。似乎最坏情况的性能是对许多短组求和，这似乎与向量的大小成线性关系

> n = 100000; x = runif(n); f = sample(n/2, n, TRUE)
> system.time(rowsum(x, f))
   user  system elapsed 
  0.228   0.000   0.229 
> n = 1000000; x = runif(n); f = sample(n/2, n, TRUE)
> system.time(rowsum(x, f)) 
   user  system elapsed 
  1.468   0.040   1.514 
> n = 10000000; x = runif(n); f = sample(n/2, n, TRUE)
> system.time(rowsum(x, f))
   user  system elapsed 
 17.369   0.748  18.166 

似乎有两个捷径可用，避免重新排序

> n = 10000000; x = runif(n); f = sample(n/2, n, TRUE)
> system.time(rowsum(x, f, reorder=FALSE))
   user  system elapsed 
 16.501   0.476  17.025 

并避免对角色的内部强制

> n = 10000000; x = runif(n); f = as.character(sample(n/2, n, TRUE)); 
> system.time(rowsum(x, f, reorder=FALSE))
   user  system elapsed 
  8.652   0.268   8.949 

然后是似乎涉及的基本操作——找出分组因子的唯一值（预先分配结果向量）并求和

> n = 10000000; x = runif(n); f = sample(n/2, n, TRUE)
> system.time( t = tabulate(f); sum(x) )
   user  system elapsed 
  0.640   0.000   0.643 

所以是的，似乎有相当大的空间可以实现更快的单一用途。这对于data.table 来说似乎很自然，并且在 C 中实现起来并不难。这是一个混合解决方案，使用 R 进行制表并使用“经典”C 接口进行求和

library(inline)

rowsum1.1 <- function(x, f) 
    t <- tabulate(f)
    crowsum1(x, f, t)


crowsum1 = cfunction(c(x_in="numeric", f_in="integer", t_in = "integer"), "
    SEXP res_out;
    double *x = REAL(x_in), *res;
    int len = Rf_length(x_in), *f = INTEGER(f_in);

    res_out = PROTECT(Rf_allocVector(REALSXP, Rf_length(t_in)));
    res = REAL(res_out);
    memset(res, 0, Rf_length(t_in) * sizeof(double));
    for (int i = 0; i < len; ++i)
        res[f[i] - 1] += x[i];
    UNPROTECT(1);
    return res_out;
")

与

> system.time(r1.1 <- rowsum1.1(x, f))
   user  system elapsed 
  1.276   0.092   1.373 

要实际返回与rowsum 相同的结果，需要将其成形为具有适当暗名称的矩阵

rowsum1 <- function(x, f) 
    t <- tabulate(f)
    r <- crowsum1(x, f, t)
    keep <- which(t != 0)
    matrix(r[keep], ncol=1, dimnames=list(keep, NULL))


> system.time(r1 <- rowsum1(x, f))
   user  system elapsed 
  9.312   0.300   9.641

所以对于所有这些工作，我们的速度只有 2 倍（而且更不通用——x 必须是数字，f 必须是整数；没有 NA 值）。是的，存在效率低下的问题，例如，分配没有计数的空间级别（尽管这避免了对名称的字符向量进行昂贵的强制转换）。

【讨论】：

哇，我印象深刻！这几乎正​​是我想要的！我的框架需要rowsum1的规格，x是数字，f是整数。无论如何，你的函数比我电脑上的 rowsum 快 4 倍。 :)

马丁。 IIRC，Rf_allocVector 没有将初始值设置为 0.0。我认为res[i] 的初始值是垃圾

@RomainFrancois 哎呀，已修复，谢谢。

另外，再考虑一下。 rowsum1 中的 which 测试可能是错误的（我猜我的答案也是同样的问题），例如x = c(1,-1,1,-1) 和 f = c(1L,1L,2L,2L)。

which(t != 0) 更好，再次感谢 Romain。更令人沮丧的是——大量的工作、不完整的实现、新的错误，而且速度仍然没有那么快。

【参考方案2】：

为了补充 Martin 的代码，这里有一些基于 Rcpp 的版本。

int increment_maybe(int value, double vec_i)
    return vec_i == 0 ? value : ( value +1 ) ;  


// [[Rcpp::export]]
NumericVector cpprowsum2(NumericVector x, IntegerVector f)
    std::vector<double> vec(10) ;
    vec.reserve(1000); 
    int n=x.size(); 
    for( int i=0; i<n; i++)
        int index=f[i]; 
        while( index >= vec.size() )
            vec.resize( vec.size() * 2 ) ;    
        
        vec[ index ] += x[i] ;
    
    // count the number of non zeros
    int s = std::accumulate( vec.begin(), vec.end(), 0, increment_maybe) ; 
    NumericVector result(s) ;
    CharacterVector names(s) ;

    std::vector<double>::iterator it = vec.begin() ;
    for( int i=0, j=0 ; j<s; j++ ,++it, ++i )
        // move until the next non zero value
        while( ! *it ) i++ ; ++it ;
        result[j] = *it ;
        names[j]  = i ;
    
    result.attr( "dim" ) = IntegerVector::create(s, 1) ;
    result.attr( "dimnames" ) = List::create(names, R_NilValue) ; 
    return result ;

C++ 代码处理所有事情，包括格式化为rowsum 给出的矩阵格式，并显示（稍微）更好的性能（至少在示例中）。

# from Martin's answer
> system.time(r1 <- rowsum1(x, f))
   user  system elapsed
  0.014   0.001   0.015

> system.time(r3 <- cpprowsum2(x, f))
   user  system elapsed
  0.011   0.001   0.013

> identical(r1, r3)
[1] TRUE

【讨论】：

干得好。想为 Rcpp Gallery 编写它吗？ ;-)

也许吧。无论如何，我必须熟悉整个过程。

没什么——不要让 git 的东西迷惑你。查看任何一篇文章，查看“来源”链接并查看其来源。带有 cmets 的 .cpp 或 .Rmd。如果您愿意，我们可以通过 Google 聊天/环聊继续。

@Dirk，你应该这样做！ :)

【参考方案3】：

这是我尝试使用Rcpp 执行此操作（第一次使用该软件包，所以请指出我的低效率）：

library(inline)
library(Rcpp)

rowsum_helper = cxxfunction(signature(x = "numeric", y = "integer"), '
  NumericVector var(x);
  IntegerVector factor(y);

  std::vector<double> sum(*std::max_element(factor.begin(), factor.end()) + 1,
                          std::numeric_limits<double>::quiet_NaN());
  for (int i = 0, size = var.size(); i < size; ++i) 
    if (sum[factor[i]] != sum[factor[i]]) sum[factor[i]] = var[i];
    else sum[factor[i]] += var[i];
  

  return NumericVector(sum.begin(), sum.end());
', plugin = "Rcpp")

rowsum_fast = function(x, y) 
  res = rowsum_helper(x, y)
  elements = which(!is.nan(res))
  list(elements - 1, res[elements])

Martin 的示例数据非常快，但仅当因子由非负整数组成时才有效，并且将消耗因子向量中最大整数的内存（对上述的一个明显改进是减去 min从最大到减少内存使用 - 这可以在 R 函数或 C++ 中完成）。

n = 1e7; x = runif(n); f = sample(n/2, n, T)

system.time(rowsum(x,f))
#    user  system elapsed 
#   14.241  0.170  14.412

system.time(tabulate(f); sum(x))
#    user  system elapsed 
#   0.216   0.027   0.252

system.time(rowsum_fast(x,f))
#    user  system elapsed 
#   0.313   0.045   0.358

另请注意，R 代码中出现了很多减速（与 tabulate 相比），因此如果您将其移至 C++，您应该会看到更多改进：

system.time(rowsum_helper(x,f))
#    user  system elapsed 
#   0.210   0.018   0.228

---

这是一个可以处理几乎所有y 的概括，但会慢一些（我实际上更喜欢在 Rcpp 中执行此操作，但不知道如何在那里处理任意 R 类型）：

rowsum_fast = function(x, y) 
  if (is.numeric(y)) 
    y.min = min(y)
    y = y - y.min
    res = rowsum_helper(x, y)
   else 
    y = as.factor(y)
    res = rowsum_helper(x, as.numeric(y))
  

  elements = which(!is.nan(res))

  if (is.factor(y)) 
    list(levels(y)[elements-1], res[elements])
   else 
    list(elements - 1 + y.min, res[elements])
  

【讨论】：

很棒的代码！非常感谢！

【参考方案4】：

在@Ben 已删除的评论和“答案”中，f 是有序且不断增加的。

n = 1e7; x = runif(n);
f <- cumsum(c(1L, sample(c(TRUE, FALSE), n - 1, TRUE)))

所以

rowsum3 <- function(x, f)

    y <- cumsum(x)
    end <- c(f[-length(f)] != f[-1], TRUE)
    diff(c(0, y[end]))

是一种常见的 R 解决方案（如果不太关心精度的话），并且

crowsum3 <- cfunction(c(x_in="numeric", f_in="integer"), "
    int j = 0, *f = INTEGER(f_in), len = Rf_length(f_in), 
        len_out = len == 0 ? 0 : f[len - 1];
    SEXP res = Rf_allocVector(REALSXP, len_out);
    double *x = REAL(x_in), *r = REAL(res);
    memset(r, 0, len_out * sizeof(double));
    for (int i = 0; i < len; ++i) 
        if (i != 0 && f[i] != f[i-1]) ++j;
        r[j] += x[i];
    
    return res;
")

可能是 C 解决方案。这些都有时间

> system.time(r3 <- rowsum3(x, f))
   user  system elapsed 
  1.116   0.120   1.238 
> system.time(c3 <- crowsum3(x, f))
   user  system elapsed 
  0.080   0.000   0.081 

R 实现中的精度损失很明显

> all.equal(r3, c3)
[1] TRUE
> identical(r3, c3)
[1] FALSE

rowsum_helper有

> system.time(r2 <- rowsum_helper(x, f))
   user  system elapsed 
  0.464   0.004   0.470 

但也假设基于 0 的索引所以

> head(rowsum_helper(x, f))
[1]       NaN 0.9166577 0.4380485 0.7777094 2.0866507 0.7300764
> head(crowsum3(x, f))
[1] 0.9166577 0.4380485 0.7777094 2.0866507 0.7300764 0.7195091

【讨论】：

这里也一样，注意res的非初始化值。

@RomainFrancois 没有，原版初始化r正确；不过crowsum3(numeric(), integer()) 没做好事！