RcppArmadillo：for循环中的负索引

Posted 2023-02-16

【中文标题】RcppArmadillo：for循环中的负索引

【英文标题】：RcppArmadillo: Negative indexing in for-loop

【发布时间】：2020-09-02 09:14:06

【问题描述】：

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

所以，我正在尝试在 RcppArmadillo 中为以下 MWE 中的 y[-i] 找到等效项：

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

非常感谢任何帮助！

【参考方案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++ 的内部循环中有一个简单的条件语句（就像基础 @987654331 中的情况一样@)。据我所知，在 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

【讨论】：

好的，谢谢！事实上，我需要对一个小矩阵进行子集化，所以这样一个带有 if 语句的 for 循环来填充矩阵是有效的。然而，这感觉有点像一种变通方法，对于大型矩阵来说有点耗时。对于更大的问题，您还会建议此解决方案吗？

@Suzanne，请参阅更新的解决方案。一般来说，C++ 中的 if 语句不如R 中的类似语句那么密集。 C++ 解决方案比我提供的优化 R 解决方案慢得多的原因是它们实现了完全不同的算法。希望这会有所帮助！