What is Mandelbrot Set?

Mandelbrot set is the set of complex numbers \(c\) for which the sequence \(\{0, f_{c}(0), f_{c}(f_{c}(0)), \ldots, f^{(n)}_{c}(0), \ldots\}\) remains bounded in absolute value where function \(f_{c}(z)=z^{2}+c\). For more detailed introduction, you could go to Wikipedia page for Mandelbrot set.

The calculation for Mandelbrot set involves heavy computation, and here we will use JuliaCall to speed up the computation by rewriting R function in Julia, and we will also see some easy and useful tips in writing higher performance Julia code. In this example, we will see that JuliaCall actually brings more than 100 times speedup of the calculation.

Pure R Implementation

The R implementation for the calculation of Mandelbrot set is quite straightforward. Note that it could be proved that if absolute value of one item in the sequence is greater than 2, the sequence \(\{0, f_{c}(0), f_{c}(f_{c}(0)), \ldots, f^{(n)}_{c}(0), \ldots\}\) will be divergent.

In the mandelbrot function, we keep the iteration number that the absolute value of the item is greater than 2 until the maximum iteration times.

And in the mandelbrotImage function, for two given sequences of \(x\) and \(y\), we calculate the value of mandelbrot function divided by the maximum iteration times on every grid point.

Performance Comparison Between R and Julia Implementation

This is the setting we use to compare the performance between R and Julia implementations.

Time for Pure R Implementation

system.time(zR <- mandelbrotImage(xs, ys, iterate_max))
#>    user  system elapsed 
#>  20.312   0.160  20.762

Time for Julia Implementation using JuliaCall

# A little warm up for Julia implementation
invisible(julia_call("mandelbrotImage", xs, ys, 2L))
system.time(zJL <- julia_call("mandelbrotImage", xs, ys, iterate_max))
#>    user  system elapsed 
#>   0.223   0.001   0.225

We could see that JuliaCall brings a lot of times speedup of the calculation, actually, we could see more speedup with larger problem scale, like 100 times speedup or even more. I won’t show the result here because I don’t want to wait minutes for this RMarkdown document to be knitted.

Tips for Julia Performance

Write Small Functions

A general advice in writing Julia (as well as R) is that you should write small functions which target at doing one thing. For example, it is possible to write mandelbrot and mandelbrotImage function together, but it is not a good practice. And the function call is also very cheap in Julia.

Type Stability

If you want to write high performance Julia code, you should write type stability functions. It means the variable in the functions should be of only one type.

For example, if you change the first line of mandelbrot functions like this:

And we do the timing again:

# A little warm up for Julia implementation
invisible(julia_call("mandelbrotImage1", xs, ys, 2L))
system.time(zJL <- julia_call("mandelbrotImage1", xs, ys, iterate_max))
#>    user  system elapsed 
#>   0.244   0.001   0.247

We could see the function becomes much slower, because in the mandelbrot1 function, z is an integer at the beginning, but becomes a complex number in the iteration. We could use @code_warntype or code_warntype and other tools provided by Julia to check about this problem, see https://docs.julialang.org/en/stable/manual/performance-tips/ for more information.


For the end of this article, let us have a look at our Mandel set!

image(xs, ys, zR, col = topo.colors(12))

image(xs, ys, zJL, col = topo.colors(12))