Tuesday, December 16, 2014

The Mandelbrot Set in R

Introduction

I was having a conversation with someone I know about weather forecasts the other day and it went something like this:

"Yes, their forecasts are really not very good. They really need to work on improving them."
"Well, yes, I think they're okay, considering what they're trying to do is impossible."

The thought that a relatively simple set of equations can give rise to infinitely complex behavior is something which has intrigue and appeal beyond those mathematically minded and academic. Fractals became a part of pop culture, and the idea of chaos is not unknown to those outside mathematical research, as it was popularized by the term "The Butterfly Effect" (though unfortunately I can't say the movie starring Ashton Kutcher was anything to write home about).

When I was finishing high school and starting my undergraduate degree I got very interested in fractals and chaos. What really spurred this interest was James Gleick's Chaos, which outlined the history of the "mathematics of chaos" in a highly readable, narrative way. So much so that I later went on to take a pure mathematics course in chaos during my degree, and wrote a program in MATLAB for exploring the Mandelbrot set.

So I thought I'd give it a shot in R.

Background

Unquestionably the most well known fractal of all is the Mandelbrot fractal. Without going into laborious mathematical detail, the Mandelbrot set stems from using complex numbers, which are numbers of the form z = x + yi. The number i, the "imaginary unit", is the special quantity such that i2 = -1. As such complex numbers are unique in that multiplying two of them together can result in much different behavior than numbers on the real number line: the magnitude of their product is not guaranteed to be greater than the terms multiplied, even for two quantities with real and imaginary parts with magnitude greater than one.

The Mandelbrot set the set of numbers produced by the following:
1. Initialize a complex set of numbers in the complex plane that are all z = 0.
2. Iterate the formula

such where c is a set of complex numbers filling the complex plane.
3. The Mandelbrot set is the set where z remains bounded for all n.

And mathematically it can be shown that if under iteration z is greater than 2 than c is not in the set.

In practice one keeps track of the number of iterations it takes z to diverge for a given c, and then colours the points accordingly to their "speed of divergence".

Execution

In R this is straightforward code. Pick some parameters about the space (range and resolution for x and y, number of iterations, etc.) and then iterate. Naively, the code could be executed as below:

However, if you know anything about writing good code in R, you know that for loops are bad, and it's better to rely upon the highly optimized underpinnings of the language by using array-based functions like which and lapply.  So a more efficient version of the same code is below:

We can then plunk these into functions where the different criteria for the rendered fractal are parameters:


Let's compare the runtime between the two shall we? For the base settings, how does the naive version compare to the vectorized one? 

> compare_runtimes()
   user  system elapsed 
 37.713   0.064  37.773 
   user  system elapsed 
  0.482   0.094   0.576 

The results speak for themselves: a ~65x decrease in runtime by using R indexing functions instead of for loops! This definitely speaks to the importance of writing optimized code taking advantage of R's vectorized functions.

You can tweak the different parameters for the function to return different parts of the complex space. Below are some pretty example plots of what the script can do with different parameters:

Conclusion

What does all this have to do with my conversation about the weather? Why, everything, of course! It's where chaos theory came from. Enjoy the pretty pictures. It was fun getting R to churn out some nice fractals and good to take a trip down memory lane.

References and Resources

The Mandelbot Set: 

code on github: