RcppArmadillo：for循环中的负索引 更新

如何解决RcppArmadillo：for循环中的负索引 更新

我是Rcpp的新手，正在尝试使用RcppArmadillo在for()循环中基于负索引进行计算。 我已经发现RcppArmadillo中的负索引并不是那么简单，但是可以通过应保留(as I found here)的元素向量来完成。当要删除的元素是循环索引时，这对我来说似乎有点困难。我尝试在this answer中实现最后一种方法，但是没有成功。是否有一种简单的方法来指定除具有循环索引的元素以外的元素向量？

因此，我试图在以下MWE中为y[-i]在RcppArmadillo中找到等效项：

# In R:
# Input 
n <- 10
y <- 1:20

# Computation
x <- rep(NA,n)
for(i in 1:n){
x[i] <- sum(y[-i])
}

到目前为止，我的代码是Rcpp代码：

// [[Rcpp::depends(RcppArmadillo)]]
#include <RcppArmadillo.h>   

// [[Rcpp::export]]
arma::vec rcpp_sum (arma::vec y,int n){

arma::vec x(n);
arma::vec ind(n);

for (i=0; i<n; i++){
ind[i] = /*no idea...*/
x[i] = sum(y[ind]);
}

return x;
}

非常感谢您的帮助！

解决方法

对于这样的任务，最好只是跳过有问题的索引。在标准C++中，我们将检查索引并跳过它。像这样：

// [[Rcpp::export]]
arma::vec rcpp_sum (arma::vec y,int n){
    
    arma::vec x(n);
    
    for (int i = 0; i < n; i++) {
        x[i] = 0; // Initialize value
        
        for (int j = 0; j < y.size(); ++j) {
            if (i != j) {
                x[i] += y[j];
            }
        }
    }
    
    return x;
}

在以上内容中，我们不再使用Sugar语法。在这种情况下，IMO没关系，因为替代方案并不过于复杂。在简化过程中，不必依赖RcppArmadillo，因为我们可以只使用纯Rcpp

#include <Rcpp.h>
using namespace Rcpp;

// [[Rcpp::export]]
NumericVector pure_rcpp_sum (NumericVector y,int n){
    
    NumericVector x(n);
    
    for (int i = 0; i < n; i++) {
        for (int j = 0; j < y.size(); ++j) {
            if (i != j) {
                x[i] += y[j];
            }
        }
    }
    
    return x;
}

验证输出：

all.equal(as.vector(rcpp_sum(y,n)),x)
[1] TRUE

all.equal(pure_rcpp_sum(y,n),x)
[1] TRUE

更新

根据OP的要求，我们在基础R中针对此特定目的采用了一种优化方法。上面的示例演示了如何解决非常具体的问题，即仅对向量中的值求和，而在C++中忽略一个值。这本质上是教学法，不一定是完成此特定任务的最佳方法（如下所示）。

在显示简单的R代码之前，我想指出，应该避免OP关于在C++的内部循环中具有简单条件语句的担忧。在基础R中）。据我所知，OP链接中所演示的子集为O（n），并具有额外逻辑向量的额外开销。我们上面提出的解决方案应该更有效，因为它在没有额外对象的情况下基本上可以完成相同的事情。

现在，获取更新的代码：

baseR <- function(y,n) {
    mySum <- sum(y)
    vapply(1:n,function(x) mySum - y[x],FUN.VALUE = 1)    
}

## Here is the OP code for reference
OP <- function(y,n) {
    x <- rep(NA,n)
    for(i in 1:n) {x[i] <- sum(y[-i])}
    x
}

就是这样。闪电般的快：

huge_y <- rnorm(1e6)
huge_n <- 1e3

system.time(t1 <- baseR(huge_y,huge_n))
 user  system elapsed 
 0.003   0.000   0.003

system.time(t2 <- pure_rcpp_sum(huge_y,huge_n))
 user  system elapsed 
2.776   0.003   2.779 

system.time(t3 <- OP(huge_y,huge_n))
 user  system elapsed 
9.555   1.248  10.805

all.equal(t1,t2)
[1] TRUE

all.equal(t1,t3)
[1] TRUE