我可以使用 Rcpp 加速我的 R 代码吗？

Question

我定义了一个包含矩阵、向量和参数的 R 函数 a。我需要针对 a 的不同值计算函数的结果。这在 R 中很容易编码，但是当矩阵为 "big" 并且参数值的数量很大时非常慢。

我可以在 R 中定义函数并在 Rcpp 中执行 for 循环吗？

它能加快计算速度吗？

R 中 foo 函数的一个最小示例是

f <- function(X,y,a){
  p = ncol(X)
  res = (crossprod(X) + a*diag(1,p))%*%crossprod(X,y)
  }

set.seed(0)
X <- matrix(rnorm(50*5),50,5)
y <- rnorm(50)
a <- seq(0,1,0.1)

result <- matrix(NA,ncol(X),length(a))

for(i in 1:length(a)){                  # Can I do this part in Rcpp?
  result[,i] <- f(X,y,a[i])
  }

result

Answer 1

我只是建议避免在循环X'X和X'y中重新计算，因为它们是循环不变的。

f <- function (XtX, Xty, a) (XtX + diag(a, ncol(XtX))) %*% Xty

set.seed(0)
X <- matrix(rnorm(50*5),50,5)
y <- rnorm(50)
a <- seq(0,1,0.1)

result1 <- matrix(NA, ncol(X), length(a))

XtX <- crossprod(X)
Xty <- crossprod(X, y)

for(i in 1:length(a)) {
  result1[,i] <- f(XtX, Xty, a[i])
  }

## compare with your `result`
all.equal(result, result1)
#[1] TRUE

小时后...

当我 return 我看到了更多需要简化的东西。

(XtX + diag(a, ncol(XtX))) %*% Xty = XtX %*% Xty + diag(a, ncol(XtX)) %*% Xty
                                   = XtX %*% Xty + a * Xty

所以实际上，XtX %*% Xty 也是循环不变的。

f <- function (XtX.Xty, Xty, a) XtX.Xty + a * Xty

set.seed(0)
X <- matrix(rnorm(50*5),50,5)
y <- rnorm(50)
a <- seq(0,1,0.1)

result2 <- matrix(NA, ncol(X), length(a))

XtX <- crossprod(X)
Xty <- c(crossprod(X, y))  ## one-column matrix to vector
XtX.Xty <- c(XtX %*% Xty)  ## one-column matrix to vector

for(i in 1:length(a)) {
  result2[,i] <- f(XtX.Xty, Xty, a[i])
  }

## compare with your `result`
all.equal(result, result2)
#[1] TRUE

事实证明我们可以摆脱循环：

## "inline" function `f`
for(i in 1:length(a)) {
  result2[,i] <- XtX.Xty + a[i] * Xty
  }

## compare with your `result`
all.equal(result, result2)
#[1] TRUE

## do it with recycling rule
for(i in 1:length(a)) {
  result2[,i] <- a[i] * Xty
  }
result2 <- XtX.Xty + result2

## compare with your `result`
all.equal(result, result2)
#[1] TRUE

## now use `tcrossprod`
result2 <- XtX.Xty + tcrossprod(Xty, a)

## compare with your `result`
all.equal(result, result2)
#[1] TRUE

我完全同意你的观点，你在问题中的示例代码只是一个 "foo"。而你发帖的时候可能没有仔细考虑过。然而，这足以表明在编写循环时，无论是 R 循环还是 C / C++ / FORTRAN 循环，我们都应该始终寻求将那些循环不变性从循环中拉出来以降低计算复杂度。

您关心的是使用 Rcpp（或任何编译语言）获得加速。您想要向量化一段不容易向量化的 R 代码。但是，"%*%"、crossprod 和 tcrossprod 映射到 BLAS（FORTRAN 代码），不是 R 级计算。您可以轻松地将所有内容矢量化。

不要总是将性能不佳归咎于 R 循环的解释开销（因为 R 是一种解释型语言）。如果每次迭代都进行一些 "heavy" 计算，例如大矩阵计算（使用 BLAS）或拟合统计模型（如 lm），那么这种开销是微不足道的。事实上，如果您确实想用编译语言编写这样的循环，请使用 lapply 函数。该函数在C层实现循环，每次迭代调用R函数。或者，是一个 Rcpp 等价物。在我看来，用R代码写的循环嵌套更有可能是低效的。

Answer 2

李哲源的回答正确地指出了在您的情况下应该做什么。至于你原来的问题，答案有两个：是的，你可以使用 Rcpp 将循环移动到 C++。不，您不会获得性能：

#include <Rcpp.h>

// [[Rcpp::export]]
Rcpp::NumericMatrix fillMatrix(Rcpp::NumericMatrix X,
                   Rcpp::NumericVector y,
                   Rcpp::NumericVector a,
                   Rcpp::Function f) {
  Rcpp::NumericMatrix result = Rcpp::no_init(X.cols(), a.length());
  for (int i = 0; i < a.length(); ++i) {
    result(Rcpp::_, i) = Rcpp::as<Rcpp::NumericVector>(f(X, y, a[i]));
  }
  return result;
}

/*** R
f <- function(X,y,a){
  p = ncol(X)
  res = (crossprod(X) + a*diag(1,p))%*%crossprod(X,y)
  }

X <- matrix(rnorm(500*50),500,50)
y <- rnorm(500)
a <- seq(0,1,0.01)

system.time(fillMatrix(X, y, a, f))
#   user  system elapsed 
#  0.052   0.077   0.075 
system.time({result <- matrix(NA,ncol(X),length(a))

for(i in 1:length(a)){
  result[,i] <- f(X,y,a[i])
  }
})
#   user  system elapsed 
#  0.060   0.037   0.049 
*/

所以在这种情况下，Rcpp 解决方案实际上比 R 解决方案慢。为什么？因为真正的工作是在函数 f 内完成的。这对于两种解决方案都是相同的，但 Rcpp 解决方案具有从 C++ 回调 R 的额外开销。注意for loops in R are not necessarily slow。顺便说一句，这里有一些基准数据：

          expr      min       lq     mean   median       uq      max neval
 fillMatrixR() 41.22305 41.86880 46.16806 45.20537 49.11250 65.03886   100
 fillMatrixC() 44.57131 44.90617 49.76092 50.99102 52.89444 66.82156   100

我可以使用 Rcpp 加速我的 R 代码吗？

Can I speedup my R code with Rcpp?

performance

for-loop

r

matrix

rcpp