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

Question

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

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

它可以加快计算速度吗？

一个foo函数的最小例子R是

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-loop 的解释开销(因为R是解释语言)归咎于性能不佳.如果每次迭代都进行一些"重"计算,如大矩阵计算(使用BLAS)或拟合统计模型(如lm),则这种开销是微不足道的.实际上,如果您确实想用编译语言编写这样的循环,请使用lapply函数.此函数在C级实现循环,并为每次迭代调用R函数.或者,拉尔夫的回答是Rcpp等价物.在我看来,用R代码编写的循环嵌套更有可能效率低下.