Bivariate Linear Approximation

Robert P. Munafo, 2023 Jul 3.

"Bivariate linear approximation" (BLA), also called "bilinear approximation", is a method to extrapolate "deltas" from two sets of reference point calculations performed using the derivative extrapolation method.

It is effectively a "series" approximation of degree 1 (linear) in two variables, using coefficients that can be iterated based on the iterated "deltas" of each of the two input iterations.

The model is "bivariate" because two deltas are used: one for the parameter c which remains constant through all iterations of a pixel (the ∂ in the cubic deltas method), and a delta representing the difference added to the reference (high-precision) z value when "rebasing".

A "rebasing" test (actually two tests, one for a c delta and one for a z delta) was introduced by "Zhuoran" on Fractal Forums in 2021 Dec.

BLA is most useful for deep zooming tasks that require a very high dwell limit, typically for dwell limits in the millions or higher. It requires an existing derivative extrapolation implementation, but need not be implemented if high dwell limits are not considered important. It is used in the more recent versions of Kalles Fraktaler and contemporaneous deep zooming programs.

The Model

Given existing iterates Z_n of a reference point (preperiod plus limit cycle) that is known to be periodic (a nucleus or Misiurewicz point) and relatively close to the rest of the pixels we want to draw, the z and c values for the points to be iterated via approximation are called z_n=Z_n+Δz_n and c_n=c+Δc respectively.

A bivariate linear model (a polynomial of degree 1 in two variables) is used to approximate each Δz_n in terms of an earlier Δz_m (with m=0 initially and increased when needed to refresh Δz_n. This "rebasing" repairs the errorneous low bits in Δz_n that have occurred because of accumulating roundoff error, by recalculating it from the known Z_n that was derived from high-precision iteration; this is necessary only when Z_n+Δz_n has gotten close to the origin.

Because of this use of two subscripts m and n we define a third amount l which is their difference:

n = m + l

Then the bivariate linear model is:

Δz_n = A_l Δz_m + B_l Δc (1)

The iteration process for iterating (A_l, B_l) is modeled by expanding the initial conditions "z₀=0" with (1) giving A₀=1, B₀=0; and expanding the recurrence "z_n+1 = z_n²+c_n" by substituting in z_n=Z_n+Δz_n and c_n=C+Δc, and then substituting (1) making sure to notice that the subscripts change from n to m (Δz_n changes to Δz_m).

In that process, terms get generated for Δz_n², Δz_nΔc, and Δc²; but there are no coefficients for those in the bilinear model so we ignore them.

The result is:

A₀ = 1
B₀ = 0
A_l+1 = 2 Z_n A_l (3)
B_l+1 = 2 Z_n B_l + 1 (4)

Converting the Modeled Iteration to Code

If we didn't have the bilinear model and just wanted to perform iterations on (Z+Δz), it would just be:

f(Z+Δz) = (Z+Δz)²+C+Δc = Z² + 2ZΔz + Δz² + C + Δc

The "Z²" and "C" have been precomputed at high precision for the reference point, so we just need to compute the new residual, which is:

f(Z+Δz)-f(Z) = 2ZΔz + Δz² + Δc (5)

In the bilinear model above there is no coefficient "C_l" for an Δz² term, so we didn't try to derive the expression to calculate for this "C_l" and therefore no way to translate it into code.

But it turns out that the coefficients A_l and B_l are not needed either, as we can show this with some algebraic manipulation.

Start with (1), but substitute in n+1 for n on the left, and l+1 for l on the right:

Δz_n+1 = A_l+1 Δz_m + B_l+1 Δ_c

Use (3) and (4) to replace A_l+1 and B_l+1:

Δz_n+1 = 2 Z_n A_l Δz_m + (2 Z_n B_l + 1) Δ_c

Distribute the Δ_c:

Δz_n+1 = 2 Z_n A_l Δz_m + 2 Z_n B_l Δ_c + Δ_c

Un-distribute the 2 Z_n:

Δz_n+1 = 2 Z_n (A_l Δz_m + B_l Δ_c) + Δ_c

Use (1) to replace the A and B part with Δz_n:

Δz_n+1 = 2 Z_n Δz_n + Δ_c

Here we have agreement with (5), except for the Δz² part which wasn't being modeled. So we didn't need coefficients A_l and B_l and so there's no problem simply using (5) for the iteration of Δz_n.

Therefore the code for the main iteration step is simply:

/* Calculate new dz */ dz := 2 * Z[n] * dz + dz * dz + c;

This matches the two published implementations (computer code). Zhuoran's code has:

dz = 2 * dz * Reference[RefIteration] + dz * dz + dc;

Claude Heiland-Allen's code is refactored slightly to save one multiplication (here with my comments added):

Z = Z_ptr[m]; // Z_ptr is the precomputed array of reference iterates z = (2 * Z + z) * z + c; // "z" is their name for Δz

BLA is discussed in this thread on FractalForums. Read through the whole thread for links to earlier discussions that are relevant to BLA.

revisions: 20230620 first version; 20230622 add model equations and derivations of iteration 20230703 rewrite and move differential estimation content to its own article

Mu-ency main page — index — recent changes — DEMZ

This page was written in the "embarrassingly readable" markup language RHTF, and was last updated on 2023 Jul 03.

s.27