R - 三元函数的向量化实现答案

【问题标题】：R - Vectorized implementation of ternary functionR - 三元函数的向量化实现
【发布时间】：2015-06-13 20:28:13
【问题描述】：

我有三个向量 X、Y 和 Z 等长 n。我需要创建一个函数f(X[i],Y[j],Z[k]) 的n x n x n 数组。执行此操作的直接方法是依次循环遍历 3 个向量中每个向量的每个元素。但是，计算数组所需的时间随着n 呈指数增长。有没有办法使用矢量化操作来实现这一点？

编辑：正如 cmets 中所述，我添加了一个简单示例来说明所需内容。

set.seed(1)
X = rnorm(10)
Y = seq(11,20)
Z = seq(21,30)

F = array(0, dim=c( length(X),length(Y),length(Z) ) )
for (i in 1:length(X))
  for (j in 1:length(Y))
    for (k in 1:length(Z))
      F[i,j,k] = X[i] * (Y[j] + Z[k])

谢谢。

【问题讨论】：

一个可重现的例子可能会有所帮助。

标签： r vectorization

【解决方案1】：

你可以使用嵌套的outer：

set.seed(1)
X = rnorm(10)
Y = seq(11,20)
Z = seq(21,30)

F = array(0, dim = c( length(X),length(Y),length(Z) ) )
for (i in 1:length(X))
  for (j in 1:length(Y))
    for (k in 1:length(Z))
      F[i,j,k] = X[i] * (Y[j] + Z[k])

F2 <- outer(X, outer(Y, Z, "+"), "*")

> identical(F, F2)
[1] TRUE

包含 Nick K 提出的expand.grid 解决方案的微基准测试：

X = rnorm(100)
Y = seq(1:100)
Z = seq(101:200)

forLoop <- function(X, Y, Z) {
  F = array(0, dim = c( length(X),length(Y),length(Z) ) )
  for (i in 1:length(X))
    for (j in 1:length(Y))
      for (k in 1:length(Z))
        F[i,j,k] = X[i] * (Y[j] + Z[k])
  return(F)
}

nestedOuter <- function(X, Y, Z) {
  outer(X, outer(Y, Z, "+"), "*")
}

expandGrid <- function(X, Y, Z) {
  df <- expand.grid(X = X, Y = Y, Z = Z)
  G <- df$X * (df$Y + df$Z)
  dim(G) <- c(length(X), length(Y), length(Z))
  return(G)
}

library(microbenchmark)
mbm <- microbenchmark(
  forLoop = F1 <- forLoop(X, Y, Z), 
  nestedOuter = F2 <- nestedOuter(X, Y, Z), 
  expandGrid = F3 <- expandGrid(X, Y, Z), 
  times = 50L)

> mbm
Unit: milliseconds
expr         min         lq        mean      median          uq        max neval
forLoop 3261.872552 3339.37383 3458.812265 3388.721159 3524.651971 4074.40422    50
nestedOuter    3.293461    3.36810    9.874336    3.541637    5.126789   54.24087    50
expandGrid   53.907789   57.15647   85.612048   88.286431  103.516819  235.45443    50

【讨论】：

很好的答案，虽然它不能推广到任意三元函数 f(X, Y, Z)
但正如您刚刚指出的那样，它比使用 expand.grid 快得多！
谢谢。这要快得多，但是可以按照尼克 K 的上述评论进行概括吗？我的代码中的函数比给出的示例更复杂。具体来说，它是F[i,j,k] = X[i] + c1*X[i]*c2 + X[i]*sqrt(V[j]*c2)*Z[k]，其中c1 和c2 是任意常量。这可以使用 Nick K 的expand.grid 方法轻松实现。
outer(X, outer(sqrt(Y * c2), Z, "*"), function(a, b) a + a * c1 * c2 + a * b) 虽然 expand.grid 方法更清晰
您可以在outer 中使用匿名函数，但是是的，也许在扩展网格上使用简单的矢量化是最好的解决方案

【解决方案2】：

这是一个额外的选项，一个可能的 Rcpp 实现（如果你喜欢你的循环）。虽然我无法超越@Juliens 解决方案（也许有人可以），但它们或多或少有相同的时间

library(Rcpp)
cppFunction('NumericVector RCPP(NumericVector X,  NumericVector Y, NumericVector Z){

             int nrow = X.size(), ncol = 3, indx = 0;
             double temp(1) ;
             NumericVector out(pow(nrow, ncol)) ;
             IntegerVector dim(ncol) ;

             for (int l = 0; l < ncol; l++){
               dim[l] = nrow;
             }             

            for (int j = 0; j < nrow; j++) {
               for (int k = 0; k < nrow; k++) {
                     temp = Y[j] + Z[k] ;
                   for (int i = 0; i < nrow; i++) {
                         out[indx] = X[i] * temp ;
                         indx += 1 ;
                   }
               }
            }

            out.attr("dim") = dim;
            return out;
}')

验证

identical(RCPP(X, Y, Z), F)
## [1] TRUE

快速基准测试

set.seed(123)
X = rnorm(100)
Y = 1:100
Z = 101:200

nestedOuter <- function(X, Y, Z) outer(X, outer(Y, Z, "+"), "*")

library(microbenchmark)
microbenchmark( 
  nestedOuter = nestedOuter(X, Y, Z),  
  RCPP = RCPP(X, Y, Z),
  unit = "relative",
  times = 1e4)

# Unit: relative
#        expr      min       lq     mean   median       uq       max neval
# nestedOuter 1.000000 1.000000 1.000000 1.000000 1.000000 1.0000000 10000
#        RCPP 1.164254 1.141713 1.081235 1.100596 1.080133 0.7092394 10000

【讨论】：

【解决方案3】：

您可以按如下方式使用 expand.grid：

df <- expand.grid(X = X, Y = Y, Z = Z)
G <- df$X * (df$Y + df$Z)
dim(G) <- c(length(X), length(Y), length(Z))
all.equal(F, G)

如果你有一个矢量化的函数，这也可以。如果没有，您可以使用 plyr::daply。

【讨论】：