Limitations of estimating local dimension and extremal index using exceedances in dynamical systems

In recent years, there has been a growing interest in extreme value theory based methods to obtain estimates of dynamical quantities. These methods stem from developments in the analysis of extreme events and their relationship with dynamical systems theory. Same as the sum of random variables has an asymptotic distributional limit described by the central limit theorem, less known theorems describe the distributional limits of extreme events. The asymptotic distributions of extremes allow the quantification of probabilities that relate the magnitude of an event with its frequency. These results have been extended to cover extremes that arise in observables computed along trajectories of chaotic dynamical systems. By defining an observable that depends on the distance between a reference point and the trajectory of the system in state space, one can obtain a representation of the extreme events as close recurrences to the reference point. These recurrences relate the dynamical and ergodic properties of the system to the statistics of the extremes. These methods have been extensively used in atmospheric data,^1–14 climate models,^15,16 and fractal dimension estimation.¹⁷ 

We focus our attention on two significant extreme value laws. The asymptotic limit of the maximum value in an observed data series partitioned in blocks follows a Generalized Extreme Value (GEV) distribution. Meanwhile, the asymptotic limit of the exceedances over a high threshold follows a Generalized Pareto Distribution (GPD). These laws have also been shown to apply to extremes of observables along orbits of dynamical systems under sufficiently fast decay of correlations.^18,19 The theory for estimating the local dimension through extreme value theory has been developed theoretically mostly for the GEV case,^20–23 with only a couple of papers dealing with the GPD estimator.^24,25 However, as typically the GPD approach requires less data, it is the one that has been mostly used in applications. However, this approach requires several abstract mathematical properties to hold, which are different from the ones required by the GEV approach. On one hand, there are some statistical properties of the dynamical system such as ergodicity and stationarity. On the other hand, there are some geometrical conditions on the invariant measure of the system, such as the existence of the local dimension and regular variation.

The dynamical quantities that are used to obtain information about the system are the local dimension and the extremal index. The setting required is a dynamical system whose evolution is contained in a bounded set of phase space, which possesses an invariant measure that describes the asymptotic frequency of visitation of sets. The local dimension is a geometrical quantity that describes the scaling of the invariant measure of a ball against its radius. This quantity is related naturally to the parameters of extreme value distributions. The extremal index is a statistical parameter that quantifies short-range correlations and appears when the extreme events cluster.

In this paper, we focus on numerically investigating the regular variation property, and the index-based estimation of the extremal index. Regular variation is a necessary property for the convergence of the exceedances to a GPD and the estimation of the local dimension. However, it is rarely mentioned and, to our knowledge, it has never been checked for any system or dataset in the literature dealing with the GPD method. Even some texts that explicitly intend to illustrate the use of this method do so in systems, which are not regularly varying. Index-based estimation is the most used method to estimate the extremal index from data. It is based on the indexed position of the extremes in the data series. It is as such a method that relies on the data series being of discrete time. However, it is routinely applied to continuous time processes discretized by sampling, a procedure that can be done in different ways and may impact the estimation.

In applications,^1–14 the output of the EBD algorithm is usually paired with an estimation of the extremal index $θ$⁠. The extremal index arises in the generalized extreme value distribution in the presence of short-range dependencies and can be seen as a number that describes the inverse of the mean size of clusters of extremes.²⁶ Note that the derivation of the exponential law for the exceedances above does not need a particular dependency structure and does not feature the extremal index in the limiting distribution. However, in most papers,^1–6,8–14 discussing applications it is included in the distribution function without a mathematical justification. The PoT approach erases the dependencies by considering the points over the threshold without regard for its position in the process. Shuffling the time series randomly will produce the exact same extremes and distribution. That is why, the extremal index should not appear in the distribution function. In contrast, the extremal index appears naturally in the block maxima approach since the clusters typically fall within a block, and thus only one extreme is chosen from each cluster, biasing the distribution.

While the extremal index can be seen as statistical quantity that is intrinsically linked to the process $U i$ and can be estimated from data, it does not appear in the asymptotic distribution of extremes through the PoT method.

Typically, the set of reference points in a chaotic invariant set that give raise to an extremal index different than 1 is the set of stable and periodic points (see, Ref. 27, Chap. 4), and thus a set of $μ$ measure zero, unlikely to be observed in most applications.

The index-based estimator for the extremal index is the most used in applications and infers the size of the clusters of exceedances by the index length between them. The following algorithm describes in detail how the computation works.

For measures defined on invariant sets of dynamical systems, it is complicated to obtain an analytical result regarding its regular variation properties. Sometimes, analytic calculations can be done for measures supported on self-similar attractors, but this applies only on a case-by-case basis. Here, we show through a series of selected examples how the regular variation of the measure relates to the geometry of the attractor. We put special emphasis on those which have been previously used in the literature related to the EBD algorithm.

A summary of the dynamical systems used, their dynamical rules, and some of their properties can be found in Table I.

We start our exposition with the Cantor shift map for two reasons. One is that, it is such a simple system that it can be shown that its invariant measure is not regularly varying. The second one is that it has been used in the literature (see, Ref. 27, Chap. 9) to illustrate the EBD method, recognizing difficulties to apply it but without identifying it as a system with a non-regularly varying invariant measure. The dynamics in the Cantor set are given by a left shift in a symbolic coding of the Cantor set. See  Appendix B, for more details and a proof of the statement that the invariant measure on the Cantor set is not regularly varying.

Figure 1 shows the result of taking the exceedances from the Cantor shift map and fitting a GPD distribution to them according to the EBD algorithm. The holes in the measure, in turn, produce holes in the distribution of exceedances as certain values are not possible.

The top panel of Fig. 2(a) shows that the quantity described in (9) does not converge to $2 − log ⁡ 2 / log ⁡ 3≈0.646$⁠. Instead, it oscillates strongly for any of the orders of magnitude reached. Since the measure is not regularly varying, we see a non-vanishing factor that prevents the convergence of the asymptotic distribution to an exponential. The middle panel shows that the estimate of the local dimension is highly correlated with the concentration of the measure and oscillates as well. This implies that the estimates of the local dimension are resolution dependent due to the lack of regular variation. The correlation estimate for the dimension shows reasonable values, and the final value is relatively close to the actual analytical solution $log⁡2/log⁡3≈0.631$⁠.

An interesting edge case is given by the so-called “fat fractals.” These are sets of positive Lebesgue measure and, thus, of integer dimension but that might have a very intricate fine structure similar to a fractal set. The example that we show here is a Cantor set, a nowhere dense set of isolated points and thus topologically equivalent to the middle third Cantor set, but it has positive Lebesgue measure. The set is produced via a dynamic construction, using a non-autonomous modification of the tent map, for more details see  Appendix C.

Despite the holes in the set and the topological equivalence to the Cantor set, the $(1/2,1/2)$ Bernoulli measure here is regularly varying. This can be seen by the convergence of the quantity $R(r)$ [Eq. (9)] in the top graph of the bottom panel of Fig. 2(c) to the expected value of $2 − 1$⁠. The estimates of the local dimension also seem to converge to the right answer, which is 1. Note that the similarities both in terms of geometry and topology between the two Cantor set examples leave as main difference the Lebesgue size in which the $(1/2,1/2)$ Bernoulli measure is concentrated.

The Cantor shift map serves as an example of the consequences of the lack of regular variation. Cantor sets are also models for the sections of attractors such as the Hénon map. In Ref. 30, it is stated that the EBD algorithm is adequate for non-Axiom A chaotic dynamical systems with multifractal properties. The Hénon map is a counterexample of this statement.

Figure 3 shows the same computation as in the previous figure but performed on a randomly selected point in the Hénon attractor. In the bottom panels, it also shows the exceedances and their geometrical disposition on four different stages of the zooming in process around $ζ$⁠. As the threshold becomes larger, the algorithm zooms more into the local structure of the attractor. A branch that looks like a line from afar, doubles, and becomes separate lines and, eventually, the one not containing $ζ$ disappears and the one containing it starts dividing into more. This process causes the quantity $R(r)$ and the estimation of the local dimension to oscillate strongly. Due to the structure of the attractor, we expect this phenomenon to arise at all sufficiently small scales, for all points.

Figure 2(b) shows the result of zooming and averaging over 1000 randomly selected reference points. The quantity $R(r)$ oscillates for any given point, at different values of the radius. There is a clear cluster around the value $2 − 1$ and corresponding $Δ ζ l=1$⁠. This means that at any small scale, the exceedances around most reference points look like a line at that resolution.

Another use of this algorithm is to compute the information dimension of the attractor (Refs. 17 and 27, Chap. 9) by averaging the local dimension estimates along a trajectory. This typically gives values around $1.35$ for the Hénon attractor, which is only one decimal point away for the value $1.26±0.02$ estimated by other authors.³¹ This is possibly a reason why the this system has been used to illustrate the EBD method in some texts (Refs. 24 and 27, Chap. 9), despite the lack of regular variation.

The Solenoid map is an example of a uniformly hyperbolic system, also called Axiom A. This class of systems have been presented in the literature as a valid setting to apply this algorithm (Refs. 24 and 27, Chap. 9). In Refs. 32 and 33, it is stated that the exceedances of the process (1) around a point in the attractor of an Axiom A system converge to an exponential distribution with parameters equal to the inverse of the local dimension around that point. Here, we produce a counterexample.

The Solenoid map is a discrete time dynamical system on a torus, which stretches the torus to be twice as long and $0<a≤1/4$ times thinner and wraps it twice inside the original torus. A detailed description of the map and a estimation procedure for its invariant measure can be found in  Appendix D.

Figure 2(d) shows the quantity $R(r)$ and the estimated local dimension for 1000 points in the Solenoid attractor, with parameter $a=0.076$⁠. The quotient of the measure of the balls becomes independent of the chosen reference point due to the uniformity of the attracting set. The result is an oscillation which synchronizes for any reference point, and the same for the estimated local dimension. As a result, there is no convergence for the algorithm even on average.

Using the approximation procedure described in  Appendix D, one can plot the measure of a ball of radius $r$ centered around a point in the attractor. This behaves approximately as a power law $r 1 − log ⁡ 2 / log ⁡ a$ but with a series of kinks spaced by powers of the parameter $a$⁠. These kinks induce the logarithmic oscillations of period $a$ and are consistent with the geometrical phenomenon observed in Fig. 3.

To compute the true recurrences in continuous time systems, we have implemented an interpolation scheme. This makes the dimension estimates 1 unit lower than in reality, see  Appendix E for an explanation of the interpolation procedure.

The Lorenz 63 system has also been used to illustrate the EBD method.^2,7,30 Figure 4(a) shows the analysis as performed with the discrete time systems but with the interpolation scheme. As the radius decreases the estimates of the measure ratio cluster around 0.5 with some downward excursions. Similarly, the estimates of the dimension cluster around 1 but have oscillations that push them upward for some range of scales. Similar to the Hénon map, the attractor is locally like a surface around most points at any given sufficiently small scale but not for any point through all the scales. The average of the local dimensions stabilizes around the value 2.054, very close to the value $2.06±0.01$ reported in Ref. 29.

An example of a less flat attractor is given by the 4D Lorenz 96 model. Figure 4(c) shows the results for this system which results in chaotic time evolution. The quantities plotted do not seem so converge or cluster in any way. A sufficiently small scale seems to be reached around $r= 10 0.5$⁠, after that the distribution of estimates appears to be stable. The range of the estimations remains very broad and with a big dispersion. Our estimation of the information dimension of the attractor over a trajectory with $10 7$ points is $2.834$⁠, which differs slightly from the average answer of $3.010$ of the run.

Figure 4(b) shows a similar plot with 25 points randomly chosen from the attractor of the Lorenz96 model in 8D. The lines corresponding to the lowest local dimension (dotted) and highest (dashed) at the end of the estimation have been highlighted. In the bottom panel of Fig. 4(c), the difference between the dotted and the dashed curves is plotted. These lines show that the oscillations do not preserve relative order or magnitude and thus inferring any property based on the value of the local dimension is an artifact of the data length. This is directly relevant for studies that associate events to particular atmospheric configurations depending on if its local dimension is high or low, such as Refs. 2,3,5–8,15. We show clearly that which event would get which classification depends on the data length if the underlying dynamics that generated the data are not regularly varying. This is because, the same initial conditions can have larger or smaller local dimension based on the data length.

The Hénon–Heiles system is a Hamiltonian system that was introduced to model the orbit of stars around galactic centers. Its orbits display a fat fractal behavior.³⁴ The initial condition provided to this system determines the total energy available to the orbits. For our simulations we have used far apart initial conditions belonging to the trajectory with initial condition (0.0, $−$0.25, 0.42, 0.0).

The results for this system are not conclusive. A problem that arises is that the numerical integration produces round up errors that do not conserve the energy of the system on the long run. This means in our case that the reference point is in a manifold from which the dynamics drift away. Based in the rest of examples, we would expect to see convergence in this case. The graphs of Fig. 4(d) could be converging around $r≈ 10 − 0.6$⁠, but they start diverging significantly for smaller radii.

A continuous time notion of the extremal index exists.³⁵ However, previous works^1–10,14 use the discrete time version of the extremal index under the assumption that it is adequate when the data are a discrete time sampling of a continuous time process. The papers use the index-based estimator given by Ref. 28, which infers the value of $θ$ through the indexes of the extremes as in a discrete time process $X 1, X 2,…$⁠, see Algorithm 2. Clustered exceedances correspond to a lower value of $θ$⁠, meanwhile longer gaps between them may increase it. Thus, the method used to discretize the continuous time process affects the estimation. An unambiguous method of discretizing is, for example, the return map to a Poincaré surface or a suspension flow. However, if one uses a numerical scheme as temporal discretization, the extremal index estimator will depend on the sampling time, the local speed of the flow, the length of the data series and the threshold choice. Denser sampling will make clusters bigger but also make the gaps between clusters bigger in the index sense. Longer trajectories may decrease the sizes of the clusters, but they will introduce more clusters and more gaps of unknown length. Increasing the threshold will make the clusters smaller and the gaps bigger, therefore increasing the value of $θ$⁠. If an adaptive time step is used, additional parameters may affect the estimation.

To study the effects of varying these parameters we have made a simple experiment with the Lorenz 63 model, the result of which can be seen in Figure 5. On the top-left panel the extremal index is computed and averaged over the process (8) using points on a trajectory while varying the time step $Δt$⁠. The temporal length is fixed at $t l=1000$ time units. Two curves have been computed, one with a fixed quantile $q=0.99$ and one with a varying quantile $q=1−1/ t l / Δ t$⁠. This quantile corresponds at taking as extremes approximately the $N$ biggest values, where $N≈ t l/Δt$ is the total amount of data in the process. For both curves the extremal index decreases strongly as the time step decreases, given that $0≤θ≤1$⁠. The bottom left panel shows a similar experiment but with a fixed time step $Δt≈0.0198$ and varying the temporal length of the trajectory. The quantiles are the same as in the top-left panel. Now, there is a difference in the behavior of the two curves. While the extremal index estimation does not change with length when the quantile stays fixed, the extremal index grows with length when the quantile also increases.

This quantity, averaged over the trajectory, is shown in the right panels of Fig. 5. In the top right panel, we see that the time duration of the clusters decreases with denser sampling, although this time the behavior is more affected by the quantile choice. When keeping the quantile fixed the t0ime spent in the cluster decreases but seems to stabilize as $Δt→0$⁠. When the quantile changes with the amount of data, the average time in clusters decreases almost linearly. This shows that in practice we cannot compare results from index-based estimators over time series sampled with different frequencies. The bottom right panel shows how the average time in the cluster changes with the length of the trajectory. If the quantile is fixed, the average time spent in clusters remains approximately constant, but if the quantile increases the average time spent in clusters decreases with the length of the trajectory.

The choice of quantile has a significant impact on the results. Generally, there is a trade-off between the quantile and the number of data points. While higher quantiles deliver a better approximation to the asymptotic limit, they also result in smaller sample sizes from which to perform statistical inference. The changing quantile chosen leaves $N$ points as extremes. The reason for this is that increasing $N$ simultaneously increases the quantile and the size of the sample, so it should give a better approximation in both fronts. This shows how sensitive the parameters are to this estimation procedure. The index-based estimator still quantifies a sort of persistence in phase space, but it is related to the spurious correlations induced by the speed of the flow around a point. This is not what the extremal index is meant to quantify. Instead is meant to quantify correlations between distinct recurrences, such as those experimented by a discrete time system when the reference point is a periodic point.

This raises the question of whether the values of extremal indices calculated using the index-based approach say something meaningful about the underlying continuous time dynamical systems or are simply numerical artifacts. Indeed, an unambiguous discretization method through a Poincaré return map assigns an extremal index of $θ=1$ to almost all points for the Lorenz 63 system attractor,³⁶ which is what should be expected given the theory of the extremal index (Ref. 27, Chap. 4). This also becomes obvious from the top-left panel of Fig. 5 as $θ→1$⁠, irrespectively of reference point, in the limit $Δt→T$⁠, where $T$ is a timescale comparable to the typical variability timescales of the dynamical system.

The regular variation property is necessary for the convergence of the exceedances to an exponential distribution. In the previous literature, this property is assumed to hold for the invariant measure of large classes of dynamical systems. Our study shows that the lack of regular variation makes the local dimension estimates resolution dependent. Thus, the estimation depends strongly on the number of data. The bias introduced by this effect can be small relative to the dimension of the system, like in the Solenoid map case, or large as in the 4D Lorenz 96 case. This bias affects differently each data point, or equivalently, each initial condition. In higher dimensional systems, another important source of bias is the one introduced by the concentration of the Euclidean norm,³³ but this affects all points equally.

There is an apparent relationship between the geometry of the invariant set and the regular variation property. It appears that when the invariant measure is singular with respect to Lebesgue measure, then it is not regularly varying. Thus, any measure whose supporting set has a dimension different from an integer would not be regularly varying. If true, this would imply that regularly varying systems might be the exception rather than the rule in high dimensional chaotic dynamics.

Index-based estimation of the extremal index is performed by constructing a method designed for discrete time systems. For continuous time systems sampled at a fixed time intervals, the index-based extremal index estimation depends heavily on parameters irrelevant to the dynamics of the system, such as the sampling frequency and the threshold (quantile). This shows that the estimator measures the speed of the flow around a point more than anything else. Thus, its deviation from the theoretically expected value of 1 for most systems is likely a numerical artifact rather than a true characterization of the dynamics. This problem is not fixed by renormalizing by the time step, as previous literature suggested.

These findings are directly relevant for the application of dynamical indicators (local dimensional and extremal index) in really complex systems, such as the output of climate models or observational data. In some other applications, analog states are used instead of recurrences in a phase space, but the theory of dynamical systems is applied to them nevertheless. It is reasonable to think that the same assumptions that apply to recurrences should apply to the analogs as well. The estimation procedures are hard to falsify.^17,30 However, they are routinely used to analyze data of enormous complexity where checking properties such that of regular variation is unfeasible, yet the findings of the literature are heavily reliant on it. We believe that additional efforts are needed to clarify the suitability of these methods for real (i.e., continuous time, high dimensional, and complex) data.

We want to thank the anonymous reviewers for helpful suggestions that greatly improved the quality of the paper. This research is supported by the European Union’s Horizon 2020 research and innovation program under the Marie Sklodowska-Curie Actions Grant No. 956170 CriticalEarth. For the purpose of open access, the author has applied a Creative Commons Attribution (CC BY) license to any Author Accepted Manuscript version arising from this submission.

The authors have no conflicts to disclose.

Ignacio del Amo: Formal analysis (equal); Investigation (equal); Methodology (equal); Software (equal); Writing – original draft (equal). George Datseris: Methodology (equal); Software (supporting); Supervision (equal); Writing – review & editing (equal). Mark Holland: Conceptualization (equal); Funding acquisition (lead); Supervision (equal); Writing – review & editing (equal).

The code and the data used to generate the plots in this paper are fully available in GitHub at https://github.com/PythagoreanCult/FractalDimensionsQ, Ref. 41.

The measure of ball a radius $r<1$ and centered in a point $ζ∈ C ∞$ can be computed in the following way (see Fig. 6). Note that for some $N∈ N$⁠, we have that $1/ 3 N + 1≤r<1/ 3 N$⁠. Regardless of the position of the center, we have that $1/ 2 N≥ μ b( B r(ζ))≥1/ 2 N + 1$ because the ball will always contain a segment of $C N + 1$⁠, with equality holding for all $x$ if $r$ is exactly $1/ 3 N + 1$⁠. When $r≤2/ 3 N + 1$⁠, the measure of the ball will depend on the position of the point, since there will be a subset of positive measure of $C ∞$ which is at $1/ 3 N + 1$ distance, but only in one of the sides, and there might be an intersection. Let $x e$ be the end point of the subset $C N + 1$ such that $ζ∈ C N + 1$ which is at distance $1/ 3 N + 1$ of the closest subset of $C N + 1$⁠. Then, the measure of the ball is $1/ 2 N + 1(1+C( 3 N + 1(r−1/ 3 N + 1− |ζ− x e |)))$⁠, where $C(x)$ is the Cantor ternary function (also called Devil’s staircase) whenever $r−1/ 3 N + 1− |ζ− x e |$ is positive, or 0 otherwise. When $r≥2/ 3 N + 1$⁠, the ball will always contain some positive measure subset of the nearest subset. The measure of the ball as before can be expressed again as $1/ 2 N + 1(1+C( 3 N + 1(r−1/ 3 N + 1− |ζ− x e |)))$ with $C(x)+1$ when the $x≥1$⁠.

The $(1/2,1/2)$ Bernoulli measure is not regularly varying. Note that because of the gaps in the support of the measure, there are values of $b$ for which the measure does not change value at all. Let $b=1/2$ and $ζ− x e=1/ 3 N + 1$ meaning $ζ$ is the opposite end point. Note that the end points will always belong to $C ∞$ by construction. If $r=2/ 3 N + 1$⁠, then Eq. (5) can only be true if the index $γ=0$⁠. Values of $b$ and $r$ such that this holds can be found for any point $ζ∈ C ∞$⁠, so decreasing sequences of $r$ such that this holds can be found.

The measure we effectively use is again the $(1/2,1/2)$ Bernoulli measure. Starting with a uniform distribution and because the maps are symmetrical with respect to $x=1/2$⁠, each iteration removes the same amount of measure from each side. Thus, the measure of each $n$-cylinder set (corresponding to a $(i/ 2 n,(i+1)/ 2 n)$ interval) is $2 − n$⁠, same as for the middle third Cantor set described above.

The proof of formula (D3) in Ref. 37, Chap. 13, is made by finding a suitable covering by balls of an iterate $f S k( T S)$⁠. This idea can be modified to estimate the volume of a small ball intersecting the Solenoid attractor. The solenoid attractor embedded in $R 3$ looks like a line coiled infinitely many times within a torus, hence the name, and the intersection of the solenoid attractor with a small ball centered at a point in the attractor has the structure of a series of branches of variable length, slightly curved, one of which goes through the center. We make the following simplifications:

1. The measure is uniform along the length of the branches.

2. We assume that the branches within the ball are approximately straight.

3. All branches are of total length $2π$ when they go around the torus.

Equation (D5) gives an approximation of where the set $f S k( T S)$ intersects $P φ k$⁠, since $v 0∈V$⁠, we have $v k Γ 1 , k∈ f S k( T S)$⁠, and thus its distance to the nearest point of $A S∩ P φ$ is bounded by $2 a k≤2/ 4 k$⁠, corresponding to the maximal diameter of a connected component of $f S k( T S)∩ P φ k$⁠.

One thing that becomes clear in the derivation of the approximation of the measure is that the point chosen to do the approximation does not matter much. This is, in formula (D4), we see that the initial angle at which we take the Poincaré section is exponentially dampened, so the angles will be almost the same for any Poincaré section taken. The initial condition within the disk $v 0∈V$ is also exponentially dampened because $a<1/4$⁠, so again it does not contribute much for high $k$⁠. This is consistent with the uniform hyperbolicity of the system, which implies that the invariant set should have locally the same structure everywhere.

Using this approximation, the measure of a ball of radius $r$ centered around a point $φ,v$ can be computed. It behaves approximately as a power law $μ S( B r(φ,v))≈ r 1 − log ⁡ 2 / log ⁡ a$⁠, but with a series of kinks spaced according to the geometrical series $r n= a n$⁠, as Fig. 7 shows. This is consistent with the geometrical phenomenon described in Fig. 3 since clusters of points in the Poincaré section $A S∩ P φ$ are separated by powers of $a$⁠.

The approximation follows closely the numerical average of the quantity (9) computed over all the points and allows us to extend it with much lower computational effort. This approximation can yield good values for very small radii, and as a verification, estimating the local dimension as the slope of the plot $log⁡ μ S( B r(φ,v))$ vs $log⁡r$ over 16 orders of magnitude gives $1.2695…$⁠, which is very close to the exact value $1−log⁡2/log⁡0.076=1.2689…$⁠, which is an accuracy hard to reach when estimating dimensions of fractal sets.

When computing recurrences, numerical schemes introduce a source of error by producing points in the trajectory, which may be far from the reference point even if the continuous time trajectory goes very close to it. As the radius of the ball decreases, the numerical approximation may even miss to hit the ball and omit a true recurrence. These problems can be avoided by performing an interpolation and solving a minimization problem that produces only the closest recurrence of the trajectory within the ball. This means that for every trajectory that enters the ball, only one point is kept. Thus, the set of points that remain in the ball has a dimension which is the dimension of the intersection between the attractor and ball minus 1. Figure 8 shows a schema of an attractor intersecting a ball centered around a reference point and the points that remain inside it after the optimization procedure takes place.