Rcpp equivalent for rowsum
Rcpp equivalent for rowsum
我正在寻找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不是很快。
为了补充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
不是答案，但可能有助于构建问题。似乎最坏情况下的性能是对许多短组求和，这似乎与向量
的大小成线性关系。> 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相同的结果，需要将其塑造为具有适当dim名称的矩阵
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值)。是的，存在效率低下，例如，分配没有计数的空间级别(尽管这避免了对名称的字符向量进行昂贵的强制转换)。
这是我使用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的示例数据来说，这是相当快的，但只有当因子由非负整数组成时才会起作用，并且会按照因子向量中最大整数的顺序消耗内存(对上述方法的一个明显改进是从max中减去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
还需要注意的是，很多的减速(与tabulate相比)发生在R代码中，所以如果你把它移到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])
}
}
在@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 has
> 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
相关文章：
- "error: no matching function for call to"构造函数错误
- 表示"accepting anything for this template argument" C++概念的通配符
- 如何在C++中从两个单独的for循环中添加两个数组
- 在Linux for Windows上编译C++代码时出错
- 调用专用模板时出错"no matching function for call to [...]"
- 为什么我的for循环不能正确获取argv
- 为什么我不能在 FOR LOOP 中使用 i/10，C++？
- Arduino:for/while/if在void setup(）或void loop(）之前?——错误:之前需要不合格
- 在基于范围的for循环中使用结构化绑定声明
- 通过for循环使用用户输入填充列表
- 使用for循环检查数组中的重复项
- 在for循环中使用auto vs decltype(vec.size(））来处理字符串的向量
- 为什么 const std：:p air<K，V>& 在 std：：map 上基于范围的 for 循环不起作用？
- 正在使用for循环创建QScatterSerie
- Python中的for循环与C++有何不同
- std：：memory_order for std：：atomic：<T>：wait
- 在更改for循环的第三部分后，未使用for循环结果
- 在 for 循环中查找问题时遇到困难
- 嵌套for循环C++的问题(初学者）
- Rcpp equivalent for rowsum
