loops - RcppArmadillo:for 循环中的负索引

Question

我是 Rcpp 的新手，正在尝试基于 for() 中的负索引执行计算- 使用 RcppArmadillo 循环。
我已经发现 RcppArmadillo 中的负索引不是那么简单，但是可以通过应该保留的元素向量来完成 (as I found here) .当要删除的元素是循环索引时，这对我来说似乎有点困难。我试图在 this answer 中实现最后一种方法，但没有成功。有没有一种简单的方法来指定不包括具有循环索引的元素的元素向量？
所以，我试图在 RcppArmadillo 中为 y[-i] 找到一个等价物。在以下 MWE 中:

# 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;
}

任何帮助是极大的赞赏!

Answer 1

对于这样的任务，最好跳过有问题的索引。标准中 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;
}

在上面，我们正在远离糖语法。 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++ 中的内部循环中有一个简单的条件语句不应该害怕(就像 base 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