Perturbation theory

  • 211 Replies
  • 8031 Views

0 Members and 1 Guest are viewing this topic.

Offline hapf

  • *
  • Fractal Friend
  • **
  • Posts: 17
« Reply #30 on: November 08, 2017, 12:38:24 AM »
Here's a clearer counterexample to that proposition...
As I said, it depends where you go for the highest period. Have a look at the dominant periods in this region. I think you can figure out where the highest period should be searched for after some testing with different locations.

Offline gerrit

  • *
  • 3f
  • ******
  • Posts: 1666
« Reply #31 on: November 13, 2017, 07:32:13 AM »
I am confused about the formula's for perturbation theory for the distance estimator using the derivative \( dz = \frac{\partial z}{\partial c} \). Can somebody check my math?

Our iteration is \( z_{k+1} = z_k^2+c \) and \( dz_{k+1} = 2z_k dz_k+1 \). We write \( c=c_0+\delta \) and the reference orbit \( X_k \) and \( dX_k \) is like z and dz but with \( c_0 \).

We then write \( z_k = w_k + X_k \) and \( dz_k = dw_k + dX_k \) which gives for the perturbation w
\( w_{k+1} = w_k^2+2X_kw_k + \delta \) and \( dw_{k+1} =2(X_k+w_k)dw_k+2w_kdX_k \).

See for example https://mathr.co.uk/mandelbrot/perturbation.pdf.

Initial conditions are \( z_1=c_0+\delta \) and \( dz_1=1 \), same for X but with \( c_0 \) and for w \( w_1=\delta \) and \( dw_1 = 0 \).

The last formula confuses me, as \( w_1=\delta = c-c_0 \), should we not have \( dw_1 = \frac{\partial w_1}{\partial c}=1 \)?? In fact my perturbation code works only if I init dw as  \( dw_1 = 1 \) which however does not seem correct according to the algebra.

Any help would be appreciated.


Offline claude

  • *
  • 3e
  • *****
  • Posts: 1082
    • mathr.co.uk
« Reply #32 on: November 13, 2017, 05:36:27 PM »
(Minor point, all this arbitrary naming of variables is really confusing to me, that's why I use angle-brackets in my PDF to distinguish deltas - maybe not so easy in the forum LaTeX, I think I use \langle and \rangle but that might need a LaTeX package not supported by MathJAX?)

You don't need to use perturbation for the derivative, because it isn't "small".  In mandelbrot-perturbator I use perturbation iterations for Z and regular iterations for dZ.

You can use series approximation to initialize the derivative, knighty told me (but I haven't implemented it yet) that the series for the derivative deltas is just the derivative of the series for the deltas, so you don't need to calculate two sets of coefficients: \[ \sum k A_k \delta^{k-1} \]

Offline gerrit

  • *
  • 3f
  • ******
  • Posts: 1666
« Reply #33 on: November 13, 2017, 06:22:31 PM »
Thanks Claude, I never thought of that. Problem fixed.

I agree with notation issue, but find \( \langle \langle z \rangle \rangle \) too cumbersome. Common notation for a small but finite disturbance in numerical literature is \( z +\Delta z \) or \( z +\delta z \).

Offline quaz0r

  • *
  • Fractal Friar
  • *
  • Posts: 106
« Reply #34 on: November 13, 2017, 08:43:58 PM »
Quote from: claude
You can use series approximation to initialize the derivative, knighty told me (but I haven't implemented it yet)

im confused.  there is more than one way to initialize the derivative is that what we're saying?  i implemented whatever you wrote about in your perturbation document back in the day and it works good.  is there a better way now?

edit:  ok looking at my code again i see how we are producing two sets of coefficients.  so in this other approach you would do nothing different/extra to initialize the derivative than to simply take the derivative of the series at the point of initialization?  it sounds too simple to be true..  >:D

edit 2:  instead of crowdfunding a professor, i think we just need to crowdfund knighty to dedicate more time to this stuff   8)
« Last Edit: November 13, 2017, 09:19:51 PM by quaz0r »

Offline 3DickUlus

  • *
  • 3e
  • *****
  • Posts: 1142
    • Digilantism
« Reply #35 on: November 14, 2017, 12:40:42 AM »
@quaz0r I agree with edit2  :thumbs:
« Last Edit: November 14, 2017, 01:25:03 AM by 3dickulus »
Resistance is fertile... you will be illuminated!

https://en.wikibooks.org/wiki/Fractals/fragmentarium

Offline gerrit

  • *
  • 3f
  • ******
  • Posts: 1666
« Reply #36 on: December 02, 2017, 03:02:29 AM »
Here's something I came up with to possibly give the Pauldelbrot glitch condition some theoretical justification. Not sure if it's sound, but here goes.

Notation: Let \( Z \leftarrow Z^2+C \) be the usual iteration (all iterations starting at 0 unless otherwise indicated). We split by writing \( C = c_0+c \) with \( c_0 \) a fixed reference point and c ranges in very small values over a rectangle. Next we write \( Z = X+z \) (X reference orbit, z perturbed orbit) which satisfy \( X \leftarrow X^2+c_0 \) and \[ z \leftarrow z^2+2Xz+c\ \ (1) \] z is in native precision, X is extended precision (but dowgraded when used in eq. (1)).

Eq.(1) seems to fail as soon as z becomes so big that the tiny c is lost in the addition, but it is not an option to then give up and continue in full precision as this happens rather soon. In practice we can just continue (1) and there is no problem unless there is one, a glitch. Pauldelbrots glitch trigger is \( |z+X|/|X|<tol \) with \( tol \approx 10^{-3} \) and works great.

So far nothing new, now for my idea. Let's say at iteration N z gets big so that the c term drops out of eq. (1) and let \( z^* = z_N, X^*=X_N, Z^* = Z_N \). and reset the iteration counter to zero (so we continue N+1 etc. but count 1,2,3..). Now look at two iterations, the correct one, and the one with the dropped c term. Let's call the correct one z and the other one v. We have
\[ z\leftarrow z^2+2Xz+c,  z_0 = z^*\ \ \ (2a) \]\[ v\leftarrow v^2+2Xv,  v_0 = z^*\ \ \ (2b) \]
Now change variables again by writing \( z = Z-X \) and \( v = V-X \) (V is the new variable) we get
\[ Z\leftarrow Z^2+c_0+c,  Z_0 = Z^*\ \ (3a) \]\[ V\leftarrow V^2+c_0,  V_0 = Z^*\ \ (3b) \]

These are just the equations for Julia sets \( J_1=J(c_0+c) \) and \( J_2=J(c_0) \), evaluated (iterated) at \( Z^* \).
Now if \( c_0 \) is on a nucleus (as it should be) J2 will be a blob centered on 0, and the orbit will remain bounded if \( Z^* \in J_2 \). Assume our image is mostly filaments, most likely \( c_0+c \) will not sit on a hyperbolic component, so the orbit will diverge in (3a) which is the behavior we want.
So it seems problems will arise if \( Z^* \) is in the Julia set of \( c_0 \). If \( R_0 \) is the size of the Julia set the error signal then becomes
\[ |w_n+X_n|<R_0 \]
which would be the Pauldelbrot measure if it was true that \( R_0 \approx tol|X_n| \).
Unfortunately I see no reason why this would be so, so I'm not sure if this makes sense or can lead to anything but I can make no further progress.

Edit: Looks like the size of the Julia set is about the size of the nucleus, so very tiny so as-is it will never work. Still, the glitches must be related to the difference in (2a) and (2b) (where the c has dropped out) and this difference can be described by these Julia sets but maybe not as simplistic as I first thought, or maybe this is a total dead-end...

Edit2: Another thing that occurred to me is that the Julia set of a periodic point is not just a connected center but has large filaments spreading out, which at finite iteration depth have an area. We know that \( X_n \in J_2 \) so if we think of \( |X_n| \) as an effective Julia set size overestimate, \( Z^* \in J_2 \) becomes \( |w_n+X_n| < tol |X_n| \) with tol = 1e-3 the "correct" overestimation factor. Maybe...
« Last Edit: December 03, 2017, 09:41:55 AM by gerrit »

Offline gerrit

  • *
  • 3f
  • ******
  • Posts: 1666
« Reply #37 on: December 18, 2017, 05:51:39 AM »
About Claude's post on the old forum: http://www.fractalforums.com/index.php?topic=25458.msg103793#msg103793

Quote
An example might help:  say you are using reference c0 and there is a glitch at iteration 12345, then the new reference c1 at the semantic center of the glitch (where z = 0, a minibrot nucleus) has period 12345 and has z = 0 at iteration 0, 12345, 24690, ...  As we are at iteration 12345, we can just adjust the delta c values by (c1 - c0), and leave the delta z values the same- they are near 0  (we found a glitch, remember) and surrounding the new reference c1 which has its z 12345 at 0.  By periodicity, we don't need to restart the iterations of the reference, its z is 0 after 12345 iterations anyway. So split the glitched pixels out into a separate batch with the new reference, the old reference is still used for the non-glitched pixels.  So the glitch correction happens at the same time as rendering, no separate passes afterwards.

Before asking questions, here's my notation: \( X_n(c_i) \) is the high precision (reference) orbit of point \( c_i \), \( z_n(c) \) is the actual computed native precision PT orbit of \( c \), \( w_n(\delta) \) is the perturbation in low precision s.t. \( z_n(c) = w_n(\delta) + X_n(c_i) \) for reference point \( c_i \), and \( c = c_i + \delta \).

We start with reference \( c_0 \) until we detect an error ("glitch") in \( w \) at a \( \delta_1 \) corresponding to \( c_1 = c_0 + \delta_1 \) at iteration N=12345. Using Pauldelbrot's criterium we have at the glitch \( |z_n(c1)/X_n(c0)|<tol\approx 10^{-3} \).

We then promote \( c_1 \) (actually a point close to it, call it \( c_1^* \)) to a reference point by searching around it for a mini with period N, and compute its orbit from N on in full precision, starting at 0. Question: do you always find this mini? As we have an "error" in \( w_n \) should we not distrust what it tells us (period is N)? Esp. since the condition is only triggered after \( w \) gets big, and the \( \delta \) term has dropped out a while ago already.

Then a new reference orbit \( X_n(c_1^*) \) computed from N on is used for "the glitched pixels".  Question: are those the pixels that glitched at the same iteration N that \( w_N(\delta_1) \) glitched?

If so, would it not be better to "hand over" more pixels that may glitch later on to \( c_1^* \), by dividing all the pixels according to whether the are closer to \( c_0 \) or \( c_1 \) using the distance metric \( d(c_1,c_2)=|z_n(c1)-z_n(c_2)| \)?

Finally, could you get away with saying at iteration N we discovered another period mini with period N, just use \( c_1 \) as new reference point, and set \( X_n(c1) = z_n(c_1) \) for \( n<N \), and \( X_N(c_1) = 0  \) (was small) and simply reuse this sequence without recomputing anything? I think not, but am not sure why.

Hope those are not too many question but I intend to implement this, and it's better to figure out things before coding, rather than after.

If you could point me to the right files on mathr.com where this is implemented that would help too. I searched but got lost...

Offline claude

  • *
  • 3e
  • *****
  • Posts: 1082
    • mathr.co.uk
« Reply #38 on: December 19, 2017, 12:30:50 AM »
Question: do you always find this mini?

I just do one iteration of Newton's method using the low precision derivative.  I don't check for convergence or do any more iterations.  Whether it improves or makes things worse sometimes I haven't checked.

Quote
As we have an "error" in \( w_n \) should we not distrust what it tells us (period is N)? Esp. since the condition is only triggered after \( w \) gets big, and the \( \delta \) term has dropped out a while ago already.

My understanding is that the irreparable damage is done on the next iteration, when the squaring pulls w to 0 while X jumps somewhere else.  In any case the shifting as I described works for me...

Quote
Then a new reference orbit \( X_n(c_1^*) \) computed from N on is used for "the glitched pixels".  Question: are those the pixels that glitched at the same iteration N that \( w_N(\delta_1) \) glitched?

Yes, "the glitched pixels" are 1. all at the same iteration of the glitch, 2. all with the same parent reference.  You end up with a tree of references and subreferences at glitches within glitches...

Quote
If so, would it not be better to "hand over" more pixels that may glitch later on to \( c_1^* \), by dividing all the pixels according to whether the are closer to \( c_0 \) or \( c_1 \) using the distance metric \( d(c_1,c_2)=|z_n(c1)-z_n(c_2)| \)?

Possibly, it's worth a try!

Quote
Finally, could you get away with saying at iteration N we discovered another period mini with period N, just use \( c_1 \) as new reference point, and set \( X_n(c1) = z_n(c_1) \) for \( n<N \), and \( X_N(c_1) = 0  \) (was small) and simply reuse this sequence without recomputing anything? I think not, but am not sure why.
Yes! The whole point is that you don't need to know X_n(c_1) for n < N, because X_N(c_1) = 0, so you can just carry on with the summed iterate as new delta.  Refining X_n(c_1) with Newton's method with the low precision derivative adds a few bits. which may be as good as "perfect" if you don't use super-high precision for the references.

Quote
If you could point me to the right files on mathr.com where this is implemented that would help too. I searched but got lost...

I don't know the .com, mine is .co.uk.  Probably these are the most relevant:
https://code.mathr.co.uk/fractal-bits/blob/HEAD:/mandelbrot-series-approximation/m.cpp#l1185
https://code.mathr.co.uk/mandelbrot-perturbator/blob/HEAD:/c/lib/perturbator.c#l559

Offline gerrit

  • *
  • 3f
  • ******
  • Posts: 1666
« Reply #39 on: December 19, 2017, 01:06:14 AM »
Quote
Finally, could you get away with saying at iteration N we discovered another period mini with period N, just use c1 as new reference point, and set Xn(c1)=zn(c1) for n<N, and XN(c1)=0 (was small) and simply reuse this sequence without recomputing anything? I think not, but am not sure why.
Yes! The whole point is that you don't need to know X_n(c_1) for n < N, because X_N(c_1) = 0, so you can just carry on with the summed iterate as new delta.  Refining X_n(c_1) with Newton's method with the low precision derivative adds a few bits. which may be as good as "perfect" if you don't use super-high precision for the references.
What I meant was instead of refining \( c_1 \), just keep it as is, set \( X_N(c1)=0 \) and use as reference orbit \( X_i(c_1) = z_{(i\ mod\ N)}(c_1) \) for \( i>N \), which we already calculated. Sounds like a free lunch.

Offline claude

  • *
  • 3e
  • *****
  • Posts: 1082
    • mathr.co.uk
« Reply #40 on: December 19, 2017, 01:35:15 AM »
Sounds like a free lunch.

Storing and retrieving the orbits (all of them, because we don't know yet which c1 is the glitch) may take longer than recalculating.  Using the single reference and low period delta as reference may be risky as it is at a glitch, which will break at the next iteration... but maybe it would be posssible to do the "glitch iteration" at high precision, then go back to native precision deltas for the later iterations?

Offline gerrit

  • *
  • 3f
  • ******
  • Posts: 1666
« Reply #41 on: December 19, 2017, 02:23:41 AM »
Storing and retrieving the orbits (all of them, because we don't know yet which c1 is the glitch) may take longer than recalculating.  Using the single reference and low period delta as reference may be risky as it is at a glitch, which will break at the next iteration... but maybe it would be posssible to do the "glitch iteration" at high precision, then go back to native precision deltas for the later iterations?
Duh, I didn't think of that, stupid.

But could you not go back to iterations 1 -- N-1 and repeat the PT calculation of \( z(c_1) = X(c_0)+w(\delta_1 \), set \( X_i(c_1) =z_i(c_1) \) i<N and \( X_N(c1)=0 \) and repeat the period, all in low precision?

Offline claude

  • *
  • 3e
  • *****
  • Posts: 1082
    • mathr.co.uk
« Reply #42 on: December 19, 2017, 02:55:59 AM »
I guess so.  Try it!  (mandelbrot-perturbator can't use that approach, as it computes in chunks, not storing the whole reference orbit at any one time)

Offline gerrit

  • *
  • 3f
  • ******
  • Posts: 1666
« Reply #43 on: December 19, 2017, 03:29:17 AM »
I guess so.  Try it!  (mandelbrot-perturbator can't use that approach, as it computes in chunks, not storing the whole reference orbit at any one time)
Easier said than done, esp. in MATLAB with no for loops, and no thread communication, but I'm working on it.

Easy to try for me is your suggestion to just switch to full precision (everywhere) for a few iterations after a glitch is detected which I did.
It looks a lot better than just ignoring the glitch, but not close to acceptable.

Thanks for the help.

Offline knighty

  • *
  • Fractal Feline
  • **
  • Posts: 185
« Reply #44 on: December 19, 2017, 10:06:29 AM »
Here's something I came up with to possibly give the Pauldelbrot glitch condition some theoretical justification. Not sure if it's sound, but here goes.

Notation: Let \( Z \leftarrow Z^2+C \) be the usual iteration (all iterations starting at 0 unless otherwise indicated). We split by writing \( C = c_0+c \) with \( c_0 \) a fixed reference point and c ranges in very small values over a rectangle. Next we write \( Z = X+z \) (X reference orbit, z perturbed orbit) which satisfy \( X \leftarrow X^2+c_0 \) and \[ z \leftarrow z^2+2Xz+c\ \ (1) \] z is in native precision, X is extended precision (but dowgraded when used in eq. (1)).

Eq.(1) seems to fail as soon as z becomes so big that the tiny c is lost in the addition, but it is not an option to then give up and continue in full precision as this happens rather soon. In practice we can just continue (1) and there is no problem unless there is one, a glitch. Pauldelbrots glitch trigger is \( |z+X|/|X|<tol \) with \( tol \approx 10^{-3} \) and works great.

What makes you think that? :)

In fact (IMHO) the glitches appear for two main complementary reasons (beside when the reference diverges prematurely):
- z being a function of c, when the magnitude of dz/dc is very small you get a kind of magnification effect that may become huge.
- cancellation:
The perturbation formula can be written like this,
\[ z \leftarrow z((z+X)+X)+c = z(Z+X)+c \]
If Z is "almost" opposite to X we get a cancellation, we loose most of the precision or all of it. When c is small (w.r.t. z and/or X) the next iterate will be almost 0. That's why The glitch detection works. That also means that one can also detect the glitch one iteration ahead.
 


xx
Resource for Learning Perturbation Theory

Started by Byte11 on Programming

0 Replies
264 Views
Last post September 10, 2018, 06:49:35 PM
by Byte11
clip
Actual Implementation of Perturbation Theory for the Mandelbrot Set

Started by LionHeart on Fractal Mathematics And New Theories

5 Replies
302 Views
Last post February 20, 2019, 07:52:36 AM
by LionHeart
xx
42 - an ironic theory of everything

Started by taurus on Fractal Mathematics And New Theories

0 Replies
119 Views
Last post March 05, 2019, 11:05:15 PM
by taurus
xx
String Theory

Started by FractalDave on Fractal Image Gallery

5 Replies
203 Views
Last post September 19, 2018, 05:39:39 PM
by gerrit
xx
algebraic number theory and the Mandelbrot set

Started by claude on Fractal Mathematics And New Theories

0 Replies
137 Views
Last post August 05, 2018, 06:32:41 PM
by claude