# 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}=0" with (1) giving A_{0}=1,
B_{0}=0; and expanding the recurrence "z_{n+1} = z_{n}^{2}+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}^{2}, Δz_{n}Δc, and
Δc^{2}; but there are no coefficients for those in the bilinear model
so we ignore them.

The result is:

A_{0} = 1

B_{0} = 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)^{2}+C+Δc = Z^{2} + 2ZΔz + Δz^{2} + C + Δc

The "Z^{2}" 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^{2} + Δc (5)

In the bilinear model above there is no coefficient "C_{l}" for an
Δz^{2} 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^{2} 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 ΔzBLA 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

From the Mandelbrot Set Glossary and Encyclopedia, by Robert Munafo, (c) 1987-2024.

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