Fibonacci Sequence in R with Memoization

This is in response to Andrew Z’s post on R-Bloggers Friday about using recursion to calculate numbers in the Fibonacci sequence. http://heuristicandrew.blogspot.com/2014/12/fibonacci-sequence-in-r-and-sas.html I’ve re-written the author’s Fibonacci function here. The only really change is that this one is extended to negative integers.

fib <- function(n) {

    # Handle "vectors" by element
    if (length(n) > 1) {
        return(sapply(n, fib))
    }

    # Base cases
    if (n == 0) 
        return(0)
    if (n == 1) 
        return(1)

    # Check to see if n is an integer Do not use is.integer as that is very
    # strict
    if (round(n, 0) != n) 
        return(NA)

    # Negative numbers
    if (n < 0) 
        return(fib(-1 * n) * ((-1)^((n + 1)%%2)))

    # Everything else
    return(fib(n - 1) + fib(n - 2))

}

This is a great example of recursion, but not very useful in practice.  It is pretty slow because it has to go to the base case several times for every top-level call.

system.time(fib(20))

##    user  system elapsed 
##    0.19    0.00    0.19

system.time(fib(25))

##    user  system elapsed 
##    1.54    0.00    1.54

system.time(fib(30))

##    user  system elapsed 
##   15.98    0.00   15.97

Using fib to calculate Fibonacci numbers over 10 is really just too slow. This function is a perfect candidate for memoization. The memoise package is great for this, but gets tricky when the function being memoized is recursive. This example uses the principles of the memoise package and even steals a little code.

fibM <- (function() {

    # The code here related to the cache *mostly* comes from the memoise
    # package's object new_cache.

    cache <- NULL

    cache_reset <- function() {
        cache <<- new.env(TRUE, emptyenv())
        cache_set('0', 0)
        cache_set('1', 1)
    }

    cache_set <- function(key, value) {
        assign(key, value, envir = cache)
    }

    cache_get <- function(key) {
        get(key, envir = cache, inherits = FALSE)
    }

    cache_has_key <- function(key) {
        exists(key, envir = cache, inherits = FALSE)
    }

    # Initialize the cache
    cache_reset()


    # This is the function that gets returned by the anonymous function and
    # becomes fibM.
    function(n) {

        nc <- as.character(n)

        # Handle "vectors" by element
        if (length(n) > 1) {
            return(sapply(n, fibM))
        }

        # Check to see if n is an integer Do not use is.integer as that is very
        # strict
        if (round(n, 0) != n) 
            return(NA)

        # Cached cases
        if (cache_has_key(nc)) 
            return(cache_get(nc))

        # Negative numbers
        if (n < 0) 
            return(fibM(-1 * n) * ((-1)^((n + 1)%%2)))

        # Everything else
        out <- fibM(n - 1) + fibM(n - 2)
        cache_set(nc, out)
        return(out)

    }

})()

This function demonstrates the use of memoization and using an anonymous function to create a closure. I wrote it this way so that the cache appears to be embedded in the function instead of leaving it in the Global Environment. This keeps it “safe” and promotes tidiness. Every time you call fibM, whether at the top level or through recursion, it checks its cache to see if it has already solved this problem. If it has it just returns the result. If not, it solves it and stores the answer in the cache for later use. Before it has been called the first time the cache is empty except for the “seed” values.

ls(environment(fibM)$cache)

## [1] '0' '1'

Call it once and see how fast it fills up.

fibM(30)

## [1] 832040

ls(environment(fibM)$cache)

##  [1] '0'  '1'  '10' '11' '12' '13' '14' '15' '16' '17' '18' '19' '2'  '20'
## [15] '21' '22' '23' '24' '25' '26' '27' '28' '29' '3'  '30' '4'  '5'  '6' 
## [29] '7'  '8'  '9'

Even on the first run this is faster than fib but that speed improvement intensifies once the cache starts to fill up. To demonstrate, I’ll clear the cache first.

environment(fibM)$cache_reset()
ls(environment(fibM)$cache)

## [1] '0' '1'

Now for some speed tests.

nums <- -25:25

system.time(res.fib <- fib(nums))

##    user  system elapsed 
##    7.11    0.00    7.11

system.time(res.fibM.1 <- fibM(nums))

##    user  system elapsed 
##       0       0       0

system.time(res.fibM.2 <- fibM(nums))

##    user  system elapsed 
##       0       0       0

Let me know if you have any questions or comments on the techniques used here.