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

Question

我定义了一个 R 函数，它包含一个矩阵、一个向量和一个参数 a。我需要针对 a 的不同值计算函数的结果。这在 R 中编码很简单，但是当矩阵“大”并且参数值的数量很大时速度很慢。

我可以在 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

---

几小时后...

当我回来时，我看到了更多需要简化的东西。

(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 是一种解释型语言)。如果每次迭代都在进行一些“繁重”的计算，例如大矩阵计算(使用 BLAS)或拟合统计模型(如 lm)，那么这种开销是微不足道的。事实上，如果您确实想用编译语言编写这样的循环，请使用 lapply 函数。该函数在C层实现循环，每次迭代调用R函数。或者，Ralf's answer是一个 Rcpp 等价物。在我看来，用 R 代码编写的循环嵌套更有可能是低效的。