CODEWITHSUNDEEP
rperformancefor-loop
stefanodv
stefanodv

Reputation: 533

Speed up iterative loop calculation with R

I have to speed up my script. I have some cycles like:

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
}

Where cost1 and cost2 are constants. Is it possible to avoid writing a cycle?

Upvotes: 5

Views: 228

Answers (3)

Martin Morgan
Martin Morgan

Reputation: 46856

The sum can be expanded as

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

and implemented as

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

Compared to my() implemented as

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
}

we have numerical equivalence and improved performance

> 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

but numerical problems (also for all other implementations) at large 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

Upvotes: 2

St&#233;phane Laurent
St&#233;phane Laurent

Reputation: 84529

Your column S is defined by a first-order linear recurrence. The i-th term can be expressed in function of i, see e.g. these slides.

> 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

Let's compare four ways. 

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

Let's compare:

> 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

The winner is f4, the one which does not use recurrence.

Upvotes: 2

minem
minem

Reputation: 3650

The main problem with your approach is that you are repeatedly calling elements of data.frame (DT$S), but that is not needed in this calculations. If we replace that with vector and add the results to data.frame at the end, it is much faster. Also we can simplify the formula.

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:

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

Upvotes: 5

Related Questions