1

I have an R function to calculate the logistic map, as below, using a for loop. But is there a way to change it (e.g. vectorise it) so that it doesn't use a loop?

logistic_map <- function(x,   # starting condition
                         r,   # rate parameter
                         N) { # number of iterations
    results <- numeric(length = N + 1)
    results[1] <- x
    for (i in seq_len(N)) {
        results[i + 1] <- r * results[i] * (1 - results[i])
    }

    data.frame(i = c(0, seq_len(N)), 
               x = results)
}

I've looked at the apply() family of functions and those in purrr, but I'm struggling to determine whether this is even possible. I'm tempted to conclude that it's not possible, just because each step relies entirely on the previous one, but it's entirely possible there's an elegant solution to this that I haven't been able to find.

Can I do this without a for loop?

Hamed
  • 206
  • 2
  • 8

2 Answers2

5
library(Rcpp)

cppFunction('NumericVector cpp_loop (NumericVector x, double r) {
  int n = x.size(), i = 0; n--;
  for (; i < n; i++) x[i + 1] = r * x[i] * (1 - x[i]);
  return x;
  }')

logistic_map <- function(x,   # starting condition
                         r,   # rate parameter
                         N) { # number of iterations
    results <- numeric(length = N + 1)
    results[1] <- x
    cpp_loop(results, r)

    data.frame(i = c(0, seq_len(N)), 
               x = results)
}

logistic_map(2, 0.2, 100)
Zheyuan Li
  • 71,365
  • 17
  • 180
  • 248
  • I can see that this would be faster as it compiles into C++, but this doesn't really vectorise the function in the sense I meant. I'm going to edit the question to make it a bit more explicit. – Hamed Oct 03 '18 at 07:31
  • 1
    @Hamed Read my answer https://stackoverflow.com/q/52597296/4891738 closely, understanding what "read-after-write" data dependency is and why SIMD-sense "vectorization" is impossible. However, implementing R-level loop as C/C++/FORTRAN level loop is R-sense "vectorization". And the latter is what is discussed in R for most of the time. – Zheyuan Li Oct 03 '18 at 07:37
  • Yes I've read it, thank you. So the fact that this has a read-after-write dependency means that it can't be vectorised in the sense I'm thinking of, but what you gave will achieve something similar: is that right? – Hamed Oct 03 '18 at 07:47
  • Ran `microbenchmark()` on both versions of the code and the C++ version is much faster, mean time of 2.5 milliseconds vs. 15.5 for the original version. – Hamed Oct 03 '18 at 08:23
2

Sure you can, here's a for-loop-free Reduce approach using only base R:

> v = {r=2.8;Reduce(function(a,b){xn=a[length(a)];b=r*xn*(1-xn);c(a,b)},rep(0,100),init=0.5)}
> v
  [1] 0.5000000 0.7000000 0.5880000 0.6783168 0.6109687 0.6655206 0.6232882
  [8] 0.6574401 0.6305953 0.6522456 0.6350996 0.6488947 0.6379250 0.6467347
 [15] 0.6397130 0.6453448 0.6408497 0.6444518 0.6415743 0.6438788 0.6420369

Whether you should is a different question. If you are doing this to try and get some speed then you should first learn to benchmark and then see if this is any faster. Using for loops is one of the things they tell you is a bad thing, but don't listen to them - sometimes for loops are faster than any of their packages.

To get more fundamental, one of the properties of recurrence relations as fractals like this is that they tend to not have closed form solutions. A closed form solution would let you compute x[i] without computing x[i-1] first, and hence be trivially vectorizable. For the logistic map, wikipedia tells us: https://en.wikipedia.org/wiki/Logistic_map that only for certain values of r do closed form solutions exist. Outside those values, you have to iteratively compute.

Spacedman
  • 92,590
  • 12
  • 140
  • 224
  • Yes I omitted the question of whether this is a good idea, but something like this gets closer to what I had in mind. Could I use `Reduce()` inside my function, so that you just pass the `x`, `r`, and `N` parameters and it would do its thing inside? – Hamed Oct 03 '18 at 07:50
  • Yep, I just hacked that on the command line. Replace the right numbers with your variables and it should all work. I expect it to be slower than your for loop though because it grows a vector. – Spacedman Oct 03 '18 at 07:57
  • Thank you. It seems like the answer to my original question is “Yeah, but probably don’t bother and just use the loop”! – Hamed Oct 03 '18 at 08:51