具有最大操作的高效双 for 循环

Question

R 中以下双 for 循环的有效实现是什么？

set.seed(1)
u <- rnorm(100, 1)
v <- rnorm(100, 2)
x <- rnorm(100, 3)
y <- rnorm(100, 4)
sum = 0
for (i in 1:100){
  for (j in 1:100) {
    sum = sum + (1 - max(u[i], v[j])) * (1 - max(x[i], y[j]))
  }
}

特别是对于非常长的向量，评估需要相当长的时间，但我想知道是否有办法对这个双 for 循环进行向量化？非常感谢。

Answer 1

我的更快。它使用outer 而不是循环，这就是它的意义所在。

首先是不需要外部包的功能，OP的，user20650's comment和我的。

original <- function(u, v, x, y){
  sum1 = 0
  for (i in seq_along(u)){
    for (j in seq_along(v)) {
      sum1 = sum1 + (1 - max(u[i], v[j])) * (1 - max(x[i], y[j]))
    }
  }
  sum1
}

comment <- function(u, v, x, y){
  sum1 = 0
  for (i in seq_along(u)){
    sum1 = sum1 + (1 - pmax(u[i], v)) * (1 - pmax(x[i], y))
  }
  sum(sum1)
}

rui <- function(u, v, x, y){
  tmp1 <- outer(u, v, pmax)
  tmp2 <- outer(x, y, pmax)
  sum((1 - tmp1) * (1 - tmp2))
}

现在是www's answer 和IceCreamToucan's answer 中的函数。

library(tidyverse)

www <- function(u, v, x, y){
  dat <- data_frame(u, v, x, y)
  dat2 <- dat %>% complete(nesting(u, x), nesting(v, y))
  SUM2 <- sum(with(dat2, (1 - pmax(u, v)) * (1 - pmax(x, y))))
  SUM2
}

IceCream <- function(u, v, x, y){
  uv <- expand.grid(u, v)
  xy <- expand.grid(x, y)
  sum((1 - do.call(pmax, uv))*(1 - do.call(pmax, xy)))
}

对它们都进行测试，看看结果是否相同。请注意，存在浮点问题。

set.seed(1234)

u <- rnorm(1e2, 1)
v <- rnorm(1e2, 2)
x <- rnorm(1e2, 3)
y <- rnorm(1e2, 4)

o <- original(u, v, x, y)
c <- comment(u, v, x, y)
w <- www(u, v, x, y)
i <- IceCream(u, v, x, y)
r <- rui(u, v, x, y)

all.equal(o, c)
all.equal(o, w)
all.equal(o, i)
all.equal(o, r)

o - c
o - w
o - r
w - r
i - r
c - r

现在是速度测试。

library(microbenchmark)
library(ggplot2)

mb <- microbenchmark(
  loop = original(u, v, x, y),
  pmax = comment(u, v, x, y),
  tidy = www(u, v, x, y),
  ice = IceCream(u, v, x, y),
  outer = rui(u, v, x, y)
)

autoplot(mb)

Answer 2

与@www 给出的类似（但在基数 R 中）

uv <- expand.grid(u, v)
xy <- expand.grid(x, y)

sum((1 - do.call(pmax, uv))*(1 - do.call(pmax, xy)))

# [1] 37270.31

基准测试

library(microbenchmark)

microbenchmark(
  original = {
    SUM <- 0
    for (i in 1:100){
      for (j in 1:100) {
        SUM <- SUM + (1 - max(u[i], v[j])) * (1 - max(x[i], y[i]))
      }
    }
  }
  , tidyverse = {
      dat <- data_frame(u, v, x, y)
      dat2 <- dat %>% complete(nesting(u, x), nesting(v, y))

      sum(with(dat2, (1 - pmax(u, v)) * (1 - pmax(x, y))))
    }
  , expand = {
      uv <- expand.grid(u, v)
      xy <- expand.grid(x, y)

      sum((1 - do.call(pmax, uv))*(1 - do.call(pmax, xy)))
    }
  , outer = sum((1 - outer(u, v, pmax))*(1 - outer(x, y, pmax)))
)

# Unit: microseconds
#       expr       min         lq       mean     median        uq        max neval
#   original 12512.838 14315.3480 18210.6801 15189.9525 17504.480 217572.149   100
#  tidyverse  4373.285  4924.0305  5812.2483  5603.1585  6044.828  14461.375   100
#     expand   843.972   961.2120  1163.5428  1061.9080  1219.674   2865.911   100
#      outer   228.823   252.7905   301.5965   285.5315   322.832    686.055   100

Answer 3

这是您的代码的输出。

set.seed(1)

u <- rnorm(100, 1)
v <- rnorm(100, 2)
x <- rnorm(100, 3)
y <- rnorm(100, 4)
SUM <- 0
for (i in 1:100){
  for (j in 1:100) {
    SUM <- SUM + (1 - max(u[i], v[j])) * (1 - max(x[i], y[j]))
  }
}
SUM
# [1] 37270.31

使用tidyverse 和pmap 可以生成相同的输出。我们首先需要为每个向量创建正确的组合。然后我们可以使用pmap 来计算结果。

library(tidyverse)

dat <- data_frame(u, v, x, y)
dat2 <- dat %>% complete(nesting(u, x), nesting(v, y))

SUM2 <- sum(with(dat2, (1 - pmax(u, v)) * (1 - pmax(x, y))))
SUM2
# [1] 37270.31

tidyversse 和 pmap 方法比 for-loop 快。

library(microbenchmark)

microbenchmark(
  m1 = {SUM <- 0
for (i in 1:100){
  for (j in 1:100) {
    SUM <- SUM + (1 - max(u[i], v[j])) * (1 - max(x[i], y[i]))
  }
}},
  m2 = {
    dat <- data_frame(u, v, x, y)
    dat2 <- dat %>% complete(nesting(u, x), nesting(v, y))

    SUM2 <- sum(with(dat2, (1 - pmax(u, v)) * (1 - pmax(x, y))))
    SUM2
  })
# Unit: milliseconds
#  expr       min        lq      mean    median        uq      max neval cld
#    m1 13.983890 15.045932 17.579693 16.554175 18.267269 39.15417   100   b
#    m2  5.716827  6.226258  7.029025  6.735946  7.186002 14.09338   100  a