Ускорьте расчет итеративного цикла с помощью R

Question

Я должен ускорить мой сценарий. У меня есть несколько циклов, как:

DT <- data.frame(Index=1:20, A=c(10:29))

cost1 <- 3
cost2 <- 0.05
cost3 <- 50

DT$S[1] <- cost1
for (j in 2:(20)) {
  DT$S[j] <- DT$S[j-1]-cost3+DT$S[j-1]*cost2/12
}

Где стоимость1 и стоимость2 являются константами. Можно ли избежать написания цикла?

Answer 1

Сумма может быть увеличена до

## S1 = c1
## S2 = S1 * c2 - c3 = c1 * c2 - c3
## S3 = S2 * c2 - c3 = c1 *c2^2 - c3 *c2 - c3
## S4 = S3 * c2 - c3 = c1 *c2^3 - c3 *c2^2 - c3 * c2 - c3

и реализовано как

f5 <- function(n) {
    c1 <- 3
    c2 <- 1 + 0.05 / 12
    c3 <- 50

    p <- cumprod(c(1, rep(c2, n - 1)))
    c1 * p - c3 * cumsum(c(0, p[-length(p)]))
}

По сравнению с my() реализовано как

my <- function(n) {
    cost1 <- 3
    cost2 <- 0.05
    cost3 <- 50

    cc <- (1 + cost2/12)      
    r <- vector('numeric', length = n)
    r[1] <- cost1
    for (j in 2:n)
        r[j] <-  r[j - 1] * cc - cost3
    r
}

у нас есть числовая эквивалентность и улучшенная производительность

> n <- 1e4
> all.equal(my(n), f5(n))
[1] TRUE
> microbenchmark(my(n), f5(n), times=5)
Unit: microseconds
  expr      min       lq      mean   median       uq      max neval
 my(n) 2495.459 2504.392 2516.5754 2505.541 2529.837 2547.648     5
 f5(n)  559.813  561.670  569.0204  563.739  565.325  594.555     5

но числовые проблемы (также для всех других реализаций) в целом n

> x = f5(1e6)
> x[which.min(x) + (-3):3]
[1] -1.778181e+308 -1.785590e+308 -1.793030e+308           -Inf           -Inf
[6]           -Inf           -Inf
> which.min(x)
[1] 168445

Answer 2

Ваш столбец S определяется линейным повторением первого порядка. I-й член может быть выражен в функции i, см., Например, эти слайды .

> DT <- data.frame(Index=1:20)
> cost1 <- 3; cost2 <- 0.05; cost3 <- 50
> DT$S[1] <- cost1
> for (j in 2:(20)) {
+   DT$S[j] <- DT$S[j-1]-cost3+DT$S[j-1]*cost2/12
+ }
> DT$S
 [1]    3.00000  -46.98750  -97.18328 -147.58821 -198.20316 -249.02901 -300.06663 -351.31691
 [9] -402.78073 -454.45898 -506.35256 -558.46236 -610.78929 -663.33424 -716.09814 -769.08188
[17] -822.28639 -875.71258 -929.36138 -983.23372
> s <- 1+cost2/12
> s_powers <- s^(1:(N-1))
> cost1*s_powers - cost3*(1-s_powers)/(1-s)
 [1]  -46.98750  -97.18328 -147.58821 -198.20316 -249.02901 -300.06663 -351.31691 -402.78073
 [9] -454.45898 -506.35256 -558.46236 -610.78929 -663.33424 -716.09814 -769.08188 -822.28639
[17] -875.71258 -929.36138 -983.23372

Давайте сравним четыре способа.

f1 <- function(){ # your way
  DT$S[1] <- cost1
  for (j in 2:N) {
    DT$S[j] <- DT$S[j-1]-cost3+DT$S[j-1]*cost2/12
  }
}
f2 <- function(){ # group the two DT$S[j-1] (cause DT$S[j-1] is slow)
  DT$S[1] <- cost1
  for (j in 2:N) {
    DT$S[j] <- (1+cost2/12)*DT$S[j-1]-cost3
  }
}
f3 <- function(){ # avoid DT$S[j-1] (@minem's answer)
  u <- numeric(N)
  u[1] <- cost1
  for (j in 2:N) {
    u[j] <- (1+cost2/12)*u[j-1]-cost3
  }
  DT$S <- u
}
f4 <- function(){ # express DT$S[j] in function of j
  s <- 1+cost2/12
  s_powers <- s^(1:(N-1))
  u2N <- cost1*s_powers - cost3*(1-s_powers)/(1-s)
  DT$S <- c(cost1, u2N)
}

Давайте сравним:

> library(microbenchmark)
> N <- 2000
> DT <- data.frame(Index=1:N)
> microbenchmark(
+   f1 = f1(),
+   f2 = f2(),
+   f3 = f3(),
+   f4 = f4()
+ )
Unit: microseconds
 expr       min        lq       mean    median         uq        max neval cld
   f1 65802.386 67920.918 73168.4472 69025.145 70347.8050 180938.153   100   c
   f2 52641.373 54790.698 58553.8418 55916.565 57021.0145 163660.112   100  b 
   f3   375.736   396.932   458.5317   418.798   459.6295    974.593   100 a  
   f4   220.890   235.170   266.3977   240.971   259.9360   1318.199   100 a  

Победитель - f4, тот, который не использует повторение.

Answer 3

Основная проблема вашего подхода заключается в том, что вы неоднократно вызываете элементы data.frame (DT$S), но в этих вычислениях это не требуется. Если мы заменим это вектором и добавим результаты в data.frame в конце, это будет намного быстрее. Также мы можем упростить формулу.

n <- 1e4
DT <- data.frame(Index = 1:n, A = seq(10, by = 1, length.out = n))

cost1 <- 3
cost2 <- 0.05
cost3 <- 50

your <- function() {
  DT$S[1] <- cost1
  for (j in 2:(n)) {
    DT$S[j] <- DT$S[j - 1] - cost3 + DT$S[j - 1]*cost2/12
  }
}
your()

Моя функция:

my <- function() {    
  cc <- (1 + cost2/12)      
  r <- vector('numeric', length = n)
  r[1] <- cost1
  for (j in 2:(n)) {
    # r[j] <- r[j - 1] - cost3 + r[j - 1] * cost2/12
    r[j] <-  r[j - 1] * cc - cost3
  }
  r
}

DT$S2 <- my()
all.equal(DT$S, DT$S2)
# [1] TRUE

microbenchmark::microbenchmark(your(), my(), times = 2)
# Unit: milliseconds
#   expr        min         lq      mean    median         uq        max neval cld
# your() 487.229621 487.229621 490.86917 490.86917 494.508715 494.508715     2   b
#   my()   1.515178   1.515178   1.59408   1.59408   1.672982   1.672982     2  a