Fractal basin boundary dimension estimation: final state uncertainty

[Datseris]

This paper: https://www.sciencedirect.com/science/article/pii/0375960183909453?via%3Dihub describes a method to estimate fractal dimensions of fractal basins boundary of attraction, which is useful because it is very hard to get the precise form of the boundary and then use stuff like genentropy (because you typically have all states that converge to an attractor, not just their boundary).

Implementation sketch:

The way I imagine it is with either by initializing several integrators in parallel (what EnsembleProblem does) or by doing reinit! for every new initial condition (what I typically do in DynamicalSystems.jl, because you typically parallelize at a higher level, e.g. for different parameters).

The user can provide a terminate function that checks the state of the integrator and terminates integration on request (e.g. terminate(integ) = norm(integ.u[3:4]) < 1e-3 a typical condition that terminates when velocity becomes too small). Once terminated, you check the final state of the integrator with another user provided function that returns true if the integrator has reached the fixed point you care about, or false for something else (and thus the point is "uncertain").

Regarding optimization, the function can take a series of sizes ε as input, because for ε2 > ε1 you don’t have to recompute what initial conditions starting in ball ε1 do; you have already done that in the ε1 step.

---

Want to back this issue? Post a bounty on it! We accept bounties via Bountysource.

[LnBn]

Another optimization may be how I envisioned my example in https://discourse.julialang.org/t/ensemble-simulations-with-paired-trajectories/45194. I.e, that instead of integrating N trajectories for each value of ε and taking the ratio of "uncertain" points to N, you instead integrate until you find a certain number of uncertain points (a colleague suggests ~10) and take 10/N* as the estimate for f, where N* is the number of trajectories you needed.

In problems where the basin geometry is very complicated (the magnetic pendulum may be a good example), this would ensure the integration doesn't spend too long on "large" values of ε where the perturbed point is more likely to be in a different basin. It would also, I think, have the bonus of better controlling the error in the estimator.

[Datseris]

The magnetic pendulum is implemented https://github.com/JuliaDynamics/DynamicalSystemsBase.jl/blob/master/src/continuous_famous_systems.jl#L736 , here is a figure I made with it:

you instead integrate until you find a certain number of uncertain points (a colleague suggests ~10) and take 10/N* as the estimate for f, where N* is the number of trajectories you needed.

I do not understand this. Can you please rephrase it or explain it in more detail? Feel free to refer to the above figure. Specifically, what qualifies as an "uncertain point"?

[LnBn]

I do not understand this. Can you please rephrase it or explain it in more detail? Feel free to refer to the above figure.

The idea of the algorithm in the paper is that you take a point x0 (which belongs to, say, the red basin), and another point on a circle of size ε centered at x0. If the new point belongs to a different basin then the original point has "final state uncertainty" (is an "uncertain point") at the resolution ε. (In the example in the paper they take two points, a "left" and "right" perturbation in the x coordinate.) The claim is that the volume f of uncertain points in the phase space should scale as ε^α, and in the paper they estimate this volume by integrating N initial conditions and their perturbations and counting the number k of those initial conditions that have final state uncertainty (for that ε). The estimate of f is then k/N. Then the value of ε is decreased and the integration runs again.

I argue that performing this many computations is unnecessary, and that you could instead set some target number T of uncertain points to find and observe the growth in the number N* of initial conditions needed to find T uncertain points as ε varies.

[Datseris]

So in your case you pick random initial conditions, and you integrate until you find T uncertain initial conditions. Then you claim that f is instead T/N with N the total you integrated?

[LnBn]

That's right, the estimate of f would then be T/N*. I believe this approach has been used in publications before. I can look for an example.

[Datseris]

Okay. I don't find any fundamental difference between the two algorithms, in the end of the day they do the same thing with more or less accuracy. Your version, even if faster, is less accurate, I believe. A uniform partition of a sphere of size ε will give a more accurate representation of the fraction of uncertain numbers. Correct?

Another point is how do you decide on the value of T? It doesn't appear obvious to me.

[LnBn]

A uniform partition of a sphere of size ε will give a more accurate representation of the fraction of uncertain numbers. Correct?

I'm not sure I understand this. How exactly will you choose the points?

Your version, even if faster, is less accurate, I believe.

Of course, that's the trade-off. Larger T would be more accurate but take longer.

Another point is how do you decide on the value of T? It doesn't appear obvious to me.

This is going to be problem specific but I expect one could get away with choosing something at least an order of magnitude less than N, perhaps more. I've been playing around with the algorithm as is from the paper so I'd be keen to compare my current results if I can find a nice way to implement my proposal here.

[Datseris]

Uniformly inside the ε-sphere.

You first choose a subset of the state space, say S, and a point of the state space, say x_0, near the boundary of attraction of S. Then you pick an ε and for a sphere of size ε you uniformly generate initial conditions in the sphere of radius ε centered at x_0.

[LnBn]

Ok I think I understand. So you're comparing N perturbations (of size ≤ ε) of a single initial condition x0, and seeing what fraction of those are in a different basin to x0?

[Datseris]

Yeap.

[awage]

Hi!

I think I may have a solution to this one. The challenging part is to obtain the basin. Once there is an estimation on a finite grid it is easy to extract the boundary. The last step consists in computing its box counting dimension with generalized_dim. It takes only a few seconds to compute. The two measures are related: https://en.wikipedia.org/wiki/Uncertainty_exponent

I have a unsupervised algorithm to compute the basins written in Julia. It works fine but still requires some fine tuning.

[Datseris]

Awesome, looking forwards to hearing more about this!