GPU Enables Obsession with Fractals

I am a fractals addict. And I now have a GPU that is enabling my addiction. GPUs were originally intended to speed up graphics, but MATLAB^® uses them to speed up computation. Fractals¹ are graphics that require extensive computation. Perfect candidates for a GPU.

The GPU has renewed my interest in the Mandelbrot set itself and facilitated work with three fractal variants: the “Burning Ship,” the “tower of powers,” and the global convergence behavior of Newton’s method.

My GPU is an NVIDIA^® Titan V, housed in a separate peripheral enclosure that is bigger than the laptop and that provides separate power and cooling. It is roughly 300 times faster than the CPU for computing these fractals, and completely changes how I interact with my mandelbrot program. I introduced gpuArray into the program and can now use a grid resolution comparable to my screen resolution and iterate until fine details in the image become visible.

Inside the Mandelbrot Set

What happens in Mandelbrot stays in Mandelbrot. The Mandelbrot involves a simple iteration with complex numbers, starting at an initial point \(z_0\). The Mandelbrot set is the region in the complex plane consisting of the values \(z_0\) for which the trajectories defined by

\[z_{k+1} = z_{k}^2 + z_{0}, k=1, 2, ...\]

remain bounded. That's it. That's the entire definition. It's amazing that such a simple definition can produce such fascinating complexity.

The color image shown above is a bird’s eye view of the Mandelbrot set. It does not have the resolution to show the richly detailed structure of the fringe just outside the boundary of the set. In fact, the Mandelbrot set has tiny filaments reaching into the fringe region, even though the fringe is hard to see in this view.

Figure 1 is a sketch of the geometry of the Mandelbrot set. The largest component is a heart-shaped cardioid, bounded by a curve with the parametric equation

\[z = e^{it}/2 - e^{2it}/4, -\pi \le t \le \pi \]

To the left of the cardioid is a circular disc of radius ¼. The cardioid and the disc touch at the point marked by the red dot. More about that red dot later.

Surrounding these two components are infinitely many nearly circular discs, and discs upon discs, with ever-decreasing diameters. The exact locations and shapes of the smallest discs can only be determined computationally.

It has recently been proved that the Mandelbrot set is mathematically connected, but the connected region is sometimes so thin that we cannot resolve it on a graphics screen or even compute it in a reasonable amount of time.

Computing Mandelbrot

The fascinating patterns that we associate with Mandelbrot come from the fringe region just outside the set. 

Here is the function that is the core of the computation. The input is a complex scalar starting point z0 and a real scalar parameter that I call depth. The output is a count kz. If the iteration is terminated because z is outside the disc of radius of two, then z was destined to become infinite and the count kz can be used as an index into a colormap. On the other hand, if z survives depth iterations, then it is declared to be within the Mandelbrot set.

function kz = mandelbrot(z0,depth)
    z = z0;
    kz = 0;
    while (z*conj(z) <= 4) & (kz <= depth)
        kz = kz + 1;
        z = z*z + z0;
    end
end

Let’s generate a grid of points in the complex plane covering a square of width w, centered at the complex point zc.

grid = 512;
s = w*(-1:1/grid:1);
[u,v] = meshgrid(s+real(zc),s+imag(zc));

Now comes the only difference between a program that runs on the CPU and one that runs on the GPU. To generate the grid on a CPU,

z0 = u + i*v;

For the GPU,

z0 = gpuArray(u + i*v);

Then, for either processor, the statement

kz = arrayfun(@mandelbrot,z0);

applies the scalar function mandelbrot to all the elements of the array z0. On the CPU, this is a vectorized double for loop over the grid. On the GPU, the statement is broken down into hundreds of individual tasks that run in parallel.

Now let’s look at two extraordinary images derived from the Mandelbrot fringe.

On the Fringe

The region known as the Valley of the Seahorses lies between the central cardioid and the disc to its left. There are actually two valleys, one above and one below the real axis. They meet at the red dot that we saw in Figure 1. With a little imagination, you can picture small marine creatures living in Figure 2. As we shall see later, the Valley also contains buried π.

Because the Mandelbrot set is self-similar, it contains an infinite number of miniature Mandelbrots, each with the same shape as the big one. The one shown in Figure 3 has a magnification factor of 10¹⁰.

Buried π

A remarkable result was discovered in 1991 by Dave Boll, then a graduate student at Colorado State University. Boll was investigating the behavior of the Mandelbrot iteration in the Valley of the Seahorses. The valleys become narrower as they approach the axis, which they meet at (-3/4,0), the red dot in Figure 1.

We can repeat Boll’s computation on a small grid centered just off the axis at the point -3/4 + yi for a tiny value of y and imaginary unit i.

y = 1.0e-07
zc = -3/4 + y*i

Make the grid just large enough to touch the axis.

width = 2*y
grid = 4

Choose depth to be inversely proportional to y, which makes it huge.

depth = 4/y

With these parameters, run the code above. It produces

kz =

40000000    40000000    40000000    40000000    40000000
40000000    40000000    40000000    40000000    40000000
40000000    40000000    31415926    40000000    40000000
40000000    40000000    20943951    40000000    40000000
40000000    40000000    15707963    40000000    40000000

Look at the iteration count in the center of the grid. See a familiar value?

This isn’t a fluke. A 2001 paper by Aaron Klebanoff, “π in the Mandelbrot Set,” analyzes a similar computation in the cusp at the front of the cardioid.

Next, a curious variant of the Mandelbrot iteration.

The Burning Ship

The Burning Ship comes from a strange iteration:

\[z_{k+1} = F(z_k) + z_{0}, k=1, 2, ...\]

where

\[F(z) = (|Re(z)| + i |Im(z)|)^2\]

I say this is strange because the function \(F(z)\) is not analytic. I am interested in this iteration because of the uncanny similarities in the following pictures.

The initial domain, shown in Figure 4, is a square of width 3.5 centered at -0.5-0.5i. I’ve inserted an arrow pointing to the region of interest in the wake of the ship.

Zoom in on the ship’s wake by a factor of 500 to the point -1.861-.002i. Apply the bone colormap to make it appear cold instead of burning. Figure 5 shows the resulting fractal next to a 1915 photograph of Antarctic explorer Ernest Shackleton’s ship Endurance frozen in the ice in the Weddell Sea.

Tower of Powers

Start with any complex number, \(z\), and repeatedly exponentiate it.

\[z, z^z, z^{z^z}, . . .\]

We can express this as an iteration. Start with \(y_0 = 1\) and let

\[y_{k+1} = z^{y_k}\]

If \(z\) is too big, this iteration will blow up to infinity. For some \(z\), it will converge to a finite limit. For example, if \(z = \sqrt 2\), the \(y_k\) will converge to 2. The most interesting case is when \(y_k\) approaches a cycle. For example, if \(z\) is near 2.5+4\(i\), the cycle has length 7.

 2.4684 + 4.0754i
-0.6216 + 0.3634i
 0.2603 - 0.0184i
 1.4868 + 0.3613i
-3.4877 + 6.1054i
 7.7632e-06 - 2.6617e-06i
 1.0000 + 0.0000i
 2.4684 + 4.0755i

This cycle length is the basis for the “tower of powers” fractal. Figure 6 shows the overall fractal.

Zooming in by a factor of 10⁵ and changing the color map produces the image shown in Figure 7.

Newton’s Method

When the starting point of Newton's method is not close to a zero of the function, the global behavior can appear to be unpredictable. Contour plots of iteration counts to convergence from a region of starting points in the complex plane generate thought-provoking fractal images.

The iteration is familiar. Pick a function \(f(z)\) with derivative \(f '(z)\). Start at \(z_0\) and let

\[z_{k+1} = z_{k} – f(z_{k})/f '(z_{k})\]

This will eventually converge to a zero of \(f\). Count the number of iterations it takes to get close.

There are many functions (and colormaps) to choose from. My favorite cubic polynomial is

\[f(z) = z^3 – 2z – 5\]

Figure 8 shows the complex plane divided into three regions where the iteration converges to one of the three zeroes of the cubic. Between these regions are areas of intense fractal action, shown in black in the figure.

Figure 9 shows the global behavior of Newton’s method seeking the zeroes of

\[f(z) = tan(sin(z)) - sin(tan(z))\]

The function has infinitely many zeroes and an unbounded first derivative.

The function for Figure 10 is

\[f(z) = z\;sin(1/z) \]

The most prominent blue regions surround the zeroes at ±1/π.

Many of the images described here were new to me—I’d never seen them before. Interactive experiments with the GPU made them possible. 

¹ The term fractal was coined by Benoit Mandelbrot (1924–2010), a Polish/French/American mathematician who spent most of his career at the IBM Watson Research Center in Yorktown Heights, N.Y.