Series Approximation and the Mandelbrot set

Series approximation has been implemented in Very Plotter's Mandelbrot set plot! After some initial calculations, Very Plotter is now able to skip some amount (often a large percentage) of the calculations needed to render a Mandelbrot set image. This vastly speeds up rendering at deep locations.

 

Here's a performance comparison with the "Flake" location discovered by Dinkydau:

No series approx: 293s With series approx: 129s

Both images were rendered on my M1 Mac mini with 7 worker threads:

How was it that much faster on the same hardware? Using 16 series approximation terms (explanation below), Very Plotter skipped 24,961 iterations for each pixel in the image. With the average pixel having 30,691 total iterations, that means we skipped 81% of the overall number of iterations in the image!

Introduction

This series approximation method (which we'll call SA) has been widely known to the Mandelbrot set community for about 9 years now, ever since "mrflay" (K. I. Martin) announced his application of perturbation theory and series approximation to the Mandelbrot set. (See my previous post on perturbation theory).

As K. I. Martin's original PDF shows, by using relatively straightforward math, we see that repeated invocations of perturbation theory calculations can be approximated with just one equation:

ฮ” n = A n ๐›ฟ + B n ๐›ฟ 2 + C n ๐›ฟ 3 + ...

(where lower-case delta, ๐›ฟ, is the initial distance between our point and the reference orbit)

In other words, the offset or delta (upper-case delta ฮ”) for any point (distance from reference orbit) at iteration n can be computed with the above equation. This is huge. Instead of having to calculate every single interation for a pixel, like we do for the basic escape time algorithm and for perturbation theory, we can just calculate where the point is after n iterations! Next, we'll see how to actually go about computing this.

The Coefficients

In the equation above, those coefficients An, Bn, ... are calculated with:

X n = reference orbit location at  n th iteration A 0 = 1 , i B 0 = 0 , i C 0 = 0 A n+1 = 2 X n A n + 1 B n+1 = 2 X n B n + A n 2 C n+1 = 2 X n C n + 2 A n B n

Following the methodology of "Botond Kosa" on fractalforums.com, where they describe the operation of the now-defunct Mandel Machine software, we can create an arbitrarily-long series approximation equation.

Again from sft_maths.pdf (original K. I. Martin perturbation theory and series approximation document):

A n+1 = 2 X n A n + 1 B n+1 = 2 X n B n + A n 2 C n+1 = 2 X n C n + 2 A n B n

Then "Botond Kosa" shows the 4th and 5th (D and E) terms in this fractalforums.com post:

D n+1 = 2 X n D n + 2 A n C n + B n 2 E n+1 = 2 X n E n + 2 A n D n + 2 B n C n

If we continue this pattern, we can calculate as many terms as we want. For example, here are the next two terms that I came up with:

F n+1 = 2 X n F n + 2 A n E n + 2 B n D n + C n 2 G n+1 = 2 X n G n + 2 A n F n + 2 B n E n + 2 C n D n

We can generalize this for any number of terms. With more terms, we can generally skip more iterations, but at the cost of many more computations required to calculate the coefficients.

Stopping Criteria

To calculate these coefficients, we need to step through the reference orbit, iteration by iteration, and calculate each term for each iteration. But how do we know when to stop and how many iterations to skip? Again, "Botond Kosa", the developer of Mandel Machine, made a telling post on fractalforums.com.

Let's look at our SA formula, where we draw a dividing line | between two terms:

ฮ” n = A n ๐›ฟ + B n ๐›ฟ 2 | C n ๐›ฟ 3 + ...

The "stopping criteria" is to stop iterating the coefficients once we can no longer draw a line between any two terms where all the terms (An๐›ฟ,Bn๐›ฟ2) on the left side of the line are not several orders of magnitude (hypotenuse length of the position of the complex number on the complex plane back to the (0+0i) origin) larger than all the terms (Cn๐›ฟ3,Dn๐›ฟ4,...) on the right side.

I believe we can make that "several orders of magnitude" number 1,000 or 10,000 but the larger it is, the sooner we'll stop. And depending on the location, we may need to stop earlier to produce an accurate image.

These stopping criteria can only be tested with some value for ๐›ฟ. So we'll need to pick an actual location in the image somewhere, calculate the initial distance from the reference orbit for that position (๐›ฟ) and use that to test the stopping criteria. Very Plotter uses a configurable number of these testing points, where at least 4 (the four corners of the image) are used (the default is 16 points, evenly spaced over the image). Since the SA is applied to all points in the image (all points in the image skip the same number of iterations), once any of the test points meets the above stopping criteria, we stop iterating our coefficients and use that nth iteration as our number of iterations to skip.

Skipping Iterations

To apply the SA, before starting perturbation iterations, we first run our SA equation with our computed coefficients and our initial ๐›ฟ value for the point. Then we simply increment our iteration number and reference orbit iteration number. That's it. We've skipped all those iterations! From there, we continue with the "Zhuoran" perturbation algorithm (see the previous perturbation theory post for details on that).