Эффективный двойной цикл с максимальной операцией

Question

Какова эффективная реализация следующего двойного цикла for в R?

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]))
  }
}

Специально для очень длинных векторов оценка занимает довольно много времени, но мне интересно, есть ли способ векторизовать этодвойной цикл?Большое спасибо.

Answer 1

Похоже на тот, что задан @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

Benchmark

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 2

Мой быстрее.Он использует outer вместо циклов, вот для чего он предназначен.

Сначала функции, которые не нуждаются во внешних пакетах, OP, один в комментарии user20650 и мой.

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ответ и в ответ IceCreamToucan .

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 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