Fractal Dimension Estimation with Persistent Homology: A Comparative Study

TL;DR

This paper introduces a new method using persistent homology to estimate fractal dimensions from point samples, demonstrating comparable or superior performance to classical techniques across various datasets.

Contribution

It presents an algorithm for persistent homology dimension estimation and compares its effectiveness with traditional methods on different fractal and empirical datasets.

Findings

Persistent homology dimension estimation performs comparably to correlation dimension.

It outperforms box-counting in fractal dimension estimation.

The method is practical for analyzing self-similar fractals and chaotic attractors.

Abstract

We propose that the recently defined persistent homology dimensions are a practical tool for fractal dimension estimation of point samples. We implement an algorithm to estimate the persistent homology dimension, and compare its performance to classical methods to compute the correlation and box-counting dimensions in examples of self-similar fractals, chaotic attractors, and an empirical dataset. The performance of the $0$-dimensional persistent homology dimension is comparable to that of the correlation dimension, and better than box-counting.

Figures

Tables

Table 1. Table 1. Data for the four self-similar fractals studied here. As discussed in Section  3.2 , lower values of α 𝛼 \alpha appear to yield better dimension estimates when comp 𝑃𝐻 i ​ ( X ) ≠ 0 . subscript comp subscript 𝑃𝐻 𝑖 𝑋 0 \text{comp}_{\mathit{PH}_{i}}\left(X\right)\neq 0.

Example	True Dim.	${\text{comp}}_{{\mathrm{𝑃𝐻}}_0}\left(X\right)$	${\text{comp}}_{{\mathrm{𝑃𝐻}}_1}\left(X\right)$	${\text{comp}}_{{\mathrm{𝑃𝐻}}_2}\left(X\right)$
$S$	$\frac{\log \left(3\right)}{\log \left(2\right)}$	$0$	$\frac{\log \left(3\right)}{\log \left(2\right)}$	$-$
$C\times I$	$1+\frac{\log \left(2\right)}{\log \left(3\right)}$	$\frac{\log \left(2\right)}{\log \left(3\right)}$	$0$	$-$
$C\times C$	$\frac{2\log \left(2\right)}{\log \left(3\right)}$	$\frac{2\log \left(2\right)}{\log \left(3\right)}$	$\frac{2\log \left(2\right)}{\log \left(3\right)}$	$-$
$M$	$\frac{\log \left(20\right)}{\log \left(3\right)}$	$0$	$\frac{\log \left(20\right)}{\log \left(3\right)}$	$\frac{\log \left(20\right)}{\log \left(3\right)}$

Table 2. Table 2. Dimension estimates for self-similar fractals, averaged over 10 10 10 trials of 10 6 superscript 10 6 10^{6} samples.

	True	Correl.	Box	${\mathrm{𝑃𝐻}}_0^{.5}$	${\mathrm{𝑃𝐻}}_0^1$	${\mathrm{𝑃𝐻}}_1^{.5}$	${\mathrm{𝑃𝐻}}_1^1$	${\mathrm{𝑃𝐻}}_2^1$	${\mathrm{𝑃𝐻}}_2^2$
$S$	$\approx 1.585$	$1.585$	$1.586$	$1.585$	$1.585$	$1.587$	$1.620$	$-$	$-$
$C\times I$	$\approx 1.631$	$1.633$	$1.618$	$1.629$	$1.623$	$1.634$	$1.642$	$-$	$-$
$C\times C$	$\approx 1.262$	$1.263$	$1.267$	$1.263$	$1.289$	$1.268$	$1.303$	$-$	$-$
$M$	$\approx 2.727$	$2.716$	$2.703$	$2.705$	$2.706$	2.878	$2.773$	$2.945$	$2.881$

Table 3. Table 3. Dimension estimates for chaotic attractors, averaged over 10 10 10 trials of 10 6 superscript 10 6 10^{6} samples.

	Correlation	Box-counting	${\mathrm{𝑃𝐻}}_0^{.5}$	${\mathrm{𝑃𝐻}}_0^1$	${\mathrm{𝑃𝐻}}_1^{.5}$	${\mathrm{𝑃𝐻}}_1^1$
Hénon	$1.21$	$1.22$	$1.28$	$1.28$	$1.52$	$1.36$
Ikeda	$1.68$	$1.72$	$1.71$	$1.72$	$1.71$	$1.72$
Rulkov	$1.01$	$1.52$	$1.62$	$1.87$	$2.02$
Lorenz	$2.04$	$>1.90$	$2.06$	$2.05$
MG	$3.04$	$>2.45$	$3.59$	$3.70$	$-$	$-$

Equations

$$E_{\alpha }^i\left(x_1,\dots ,x_n\right)=\sum _{I\in {PH}_i\left(x_1,\dots ,x_n\right)}{\left|I\right|}^{\alpha }.$$

$$X_{ϵ}=\left\{y\in {\mathbb{R}}^m:d(y,X)<ϵ\right\}.$$

$$E_{\alpha }^i\left(X\right)=\sum _{I\in {PH}_i\left(X\right)}{\left|I\right|}^{\alpha },$$

$${\text{dim}}_{{PH}_i^{\alpha }}\left(\mu \right)=\frac{\alpha }{1-\beta },$$

$$\beta =\underset{n\to \infty }{\mathrm{limsup}}\frac{\log \left(\mathbb{E}\left(E_{\alpha }^i\left(x_1,\dots ,x_n\right)\right)\right)}{\log \left(n\right)}.$$

$$F_i\left(X,ϵ\right)=\#\left\{I\in {PH}_i\left(X\right):\left|I\right|>ϵ\right\}$$

$${\text{comp}}_{{PH}_i}\left(X\right)=\underset{ϵ\to 0}{\lim }\frac{-\log \left(F_i\left(X,ϵ\right)\right)}{\log \left(ϵ\right)}.$$

$$C\left(ϵ\right)=\mathbb{P}\left(d\left(x,y\right)<ϵ\right).$$

$$\underset{ϵ\to 0}{\lim }\frac{\log \left(C\left(ϵ\right)\right)}{\log \left(ϵ\right)}$$

$$\underset{ϵ\to 0}{\lim }\underset{n\to \infty }{\lim }\frac{\log C\left(n,ϵ\right)}{\log ϵ}$$

$$C\left(n,ϵ\right)=\frac{\#\{\left(x_i,x_j\right):d\left(x_i,x_j\right)<ϵ,i<j\}}{n\left(n-1\right)},$$

$$C(N,{ϵ}_1)$$

$${\text{dim}}_{\overline{\text{box}}}\left(X\right)=\underset{\delta \to 0}{\mathrm{limsup}}-\frac{N_{\delta }\left(X\right)}{\delta }\text{and}{\text{dim}}_{\underline{\text{box}}}\left(X\right)=\underset{\delta \to 0}{\mathrm{liminf}}-\frac{N_{\delta }\left(X\right)}{\delta },$$

$$N_{j/10000}\left({\text{x}}_m\right)={\text{min}}_{1\le k\le 4}\left(N_{j/10000}\left({\text{x}}_{m-k}\right)\right).$$

$$\left(x,y\right)=\left(\sum _{i:a_i=1}{2}^{-i}+\sum _{i:a_i=2}{2}^{-i-1},\sum _{i:a_i=2}\sqrt{3}{2}^{-i-1}\right).$$

$$x=2\sum _{i=1}^{\infty }{a}_i{3}^{-i}.$$

$$\left(x,y,z\right)=\sum _{i=0}^{\infty }{3}^{-i}\left(x_i,y_i,z_i\right).$$

$$x_{n+1}$$

$$y_{n+1}$$

$$z_{n+1}=a+R\text{mbox}exp\left[i\left(\varphi -\frac{p}{1+|z_n{|}^2}\right)\right]z_n.$$

$$a$$

$$x_{n+1}$$

$$y_{n+1}$$

$$\mu$$

$$\dot{x}$$

$$\dot{y}$$

$$\dot{z}$$

$$\rho$$

$$y^{\prime }(t)$$

$$a$$

Full text

Research of the first author was supported in part by NIH T32 NS007292 and NSF DMS-1440140 while the author was in residence at the Mathematical Sciences Research Institute in Berkeley, California, during the Fall 2018 semester.

Research of the second author was supported in part by a NSF Mathematical Sciences Postdoctoral Research Fellowship under award number DMS-1606259.

1. Introduction

Loosely speaking, fractal dimension measures how local properties of a set depend on the scale at which they are measured. The Hausdorff dimension was perhaps the first precisely defined notion of fractal dimension [41, 23]. It is difficult to estimate in practice, but several other more computationally practicable definitions have been proposed, including the box-counting [10] and correlation [34] dimensions. These notions are in-equivalent in general.

Popularized by Mandelbrot [57, 58], fractal dimension has applications in a variety of fields including materials science [17, 92, 44], biology [6, 52, 66], soil morphology [68], and the analysis of large data sets [7, 86]. It is also important in pure mathematics and mathematical physics, in disciplines ranging from dynamics [81] to probability [8]. In some applications, it is necessary to estimate the dimension of a set from a point sample. These include earthquake hypocenters and epicenters [47, 39], rain droplets [54, 28], galaxy locations [37], and chaotic attractors [81, 34].

We propose that the recently defined persistent homology dimensions [1, 77] (Definition 1 below) are a practical tool for dimension estimation of point samples. Persistent homology [20] quantifies the shape of a geometric object in terms of how its topology changes as it is thickened; roughly speaking, $i$-dimensional persistent homology ($\mathit{PH}_{i}$) tracks $i$-dimensional holes that form and disappear in this process. Recently, it has found many applications in fields ranging from materials science [43, 73] to biology [31, 89]. There are further applications of persistent homology in machine learning, for which different methods of vectorizing persistent homology have been defined [2, 11]. In most applications, larger geometric features represented by persistent homology are of greatest interest. However, in the current context it is the smaller features — the “noise” — from which the dimension may be recovered.

For a finite subset of a metric space $\left\{x_{1},\ldots,x_{n}\right\}$ and a positive real number $\alpha$ define the $\alpha$ weight $E_{\alpha}^{i}\left(x_{1},\ldots,x_{n}\right)$ to be the sum of the lengths of the (finite) $\mathit{PH}_{i}$ intervals to the $\alpha$ power:

[TABLE]

If $\mu$ is a measure on bounded metric space, $\left\{x_{j}\right\}_{j\in\mathbb{N}}$ are i.i.d. samples from $\mu,$ and $i$ is a natural number the $\mathit{PH}_{i}$-dimension $\text{dim}_{\mathit{PH}_{i}^{\alpha}}\left(\mu\right)$ measures how $E_{\alpha}^{i}\left(x_{1},\ldots,x_{n}\right)$ scales as $n\rightarrow\infty.$ In particular, if $E_{\alpha}^{i}\left(x_{1},\ldots,x_{n}\right)\approx n^{\frac{d-\alpha}{d}},$ then $\text{dim}_{\mathit{PH}_{i}^{\alpha}}\left(\mu\right)=d.$ See below for a precise definition.

We estimate the persistent homology dimension of several examples, and compare its performance to classical methods to estimate the correlation and box-counting dimensions. We study the convergence of estimates as the sample size increases, and the variability of the estimate between different samples. The examples we consider are from three broad classes: self-similar fractals, chaotic attractors, and empirical data. The sets in the former class (shown in Figure 1– 1) have known dimensions, and are regular in a sense that implies that the various notions of fractal dimension agree for them. That is, they have a single well-defined “fractal dimension.” This makes it easy to compare the performance of the different dimension estimation techniques. In those cases, our previous theoretical results [77] imply that at least the zero-dimensional version of persistent homology dimension will converge to the true dimension.

The second class of fractals we study are strange attractors arising from chaotic dynamical systems, see e.g. Figure 1-1. Generically, these sets are not known to be regular, and the various definitions of fractal dimension may disagree for them. Finally, we apply the dimension estimation techniques to the Hauksson–Shearer Southern California earthquake catalog [40].

We also propose that a notion of persistent homology complexity due to MacPherson and Schweinhart [56] may be a good indicator of how difficult it is to estimate the correlation dimension or persistent homology dimension of a shape.

In summary, we conclude the following.

• •

Effectiveness. Based on our experiments, the $\mathit{PH}_{0}$ and correlation dimensions perform comparably well. In cases where the true dimension is known, they approach it at about the same rate. In most cases, the box-counting, $\mathit{PH}_{1},$ and $\mathit{PH}_{2}$ dimensions perform worse.

• •

Efficiency. Computation of the $\mathit{PH}_{0}$ dimension is fast and comparable with the correlation and box-counting dimensions. The $\mathit{PH}_{1}$ dimension is reasonably fast subsets of $\mathbb{R}^{2},$ but the $\mathit{PH}_{1}$ and $\mathit{PH}_{2}$ dimensions are quite slow for point clouds in $\mathbb{R}^{3},$ and computations for sets with higher ambient dimension are impractical.

• •

Equivalence. For a large class of regular fractals the $\mathit{PH}$ dimension coincides with various classical definitions of fractal dimension. However, there are fractal sets for which these definitions do not agree, with sometimes surprisingly large differences (e.g. the Rulkov and Mackey-Glass attractors). Within the class of $\mathit{PH}_{i}$ dimensions, there is variation between dimension estimates coming from different choices of the homological dimension or scaling weight. All this said, it can be difficult to determine whether dimensions truly disagree, or if the convergence is very slow.

• •

Error Error estimates, whether they are the “statistical test error” of the correlation dimension or the empirical standard deviation of estimates between trials, do not contain meaningful information about the difference between the dimension estimate and the true dimension. In general, it is difficult to tell whether a dimension estimate has approached its limiting value.

• •

Ease-of-use We found one simple rule for fitting a power law to estimate the $\mathit{PH}_{0}$ which worked well for all examples, in contrast to the correlation dimension and (especially) the box-counting dimension.

In the following, we briefly survey previous work comparing different methods for the estimation of fractal dimension (Section 1.1), outline the different methods considered here (Section 2), and discuss the results for each example (Sections 3, 4, and 5).

1.1. Background and Previous Work

Many previous studies have compared different methods for fractal dimension estimation. These include surveys focusing on applications to chaotic attractors [83], medical image analysis [52], networks [70], and time series and spatial data [32]. Fractal dimension estimation techniques have also been applied for intrinsic dimension analysis in applications where an integer-valued dimension is assumed [13, 63].

Several previous studies have observed a relationship between fractal dimension and persistent homology. Estimators based on [math]-dimensional persistent homology (minimum spanning trees) were proposed by Weygart et al. [87] and Martinez et al. [61]; we explain their relationship to the current work after Definition 1 below. The PhD thesis of Robins [69], arguably the first publication in the field of topological data analysis, studied persistent Betti numbers of fractals and proved results for the [math]-dimensional homology of totally disconnected sets. In 2012, MacPherson and Schweinhart introduced an alternate definition of $\mathit{PH}$ dimension, which we refer to here as “$\mathit{PH}_{i}$ complexity” to avoid confusion (Section 2.1.3). The two previous notions measure the complexity of a shape rather than a classical dimension; they are trivial for $\mathbb{R}^{m},$ for example. A 2019 paper of Adams et al. [1], proposed a $\mathit{PH}$ dimension that is a special case of the one we consider here, and performed computational experiments on self-similar fractals. However, they do not compare the estimates with those of other fractal dimension estimation techniques. In 2018, Schweinhart [78] proved a relationship between the upper box-counting dimension and the extremal properties of the persistent homology of a metric space, which was a stepping stone to the paper where the current definition was introduced [77]. The latter paper did not undertake a computational study of the persistent homology dimension; we mention relevant theoretical results in Section 2.1.1.

The second class of fractals we study arise in dynamical systems. A dynamical system describes the motion of trajectories $x(t)$ within a phase space, with a rule for describing how $x(t)$ changes as time $t$ increases, and are referred to as either discrete or continuous depending on whether $t\in\mathbb{Z}_{+}$ or $t\in\mathbb{R}_{+}$. An attractor $A$ is an invariant set (that is, if $x(t)\in A,$ then $x(s)\in A$ for all $s>t$) in the phase space such that if $x(t)$ is close to $A$, then its distance to $A$ will decrease to zero asymptotically as $t\to\infty$. Some attractors are simple, consisting of a collection of points or smooth periodic orbits, but some are complicated, such as those shown in Figure 1 (d-f). In dissipative systems strange attractors — attractors with a fractional dimension — can occur.

For self-similar fractals, dimension is easy to compute exactly and approximate computationally. The situation for strange attractors is quite different; in general, the fractal dimension of an attractor cannot be computed analytically. The measures on these sets may not be regular, and the different notions of dimension appear to disagree in many cases. Note that when estimates for different notions of fractal dimension disagree, it is difficult to tell whether this is because (1) one method converges much faster than the other (2) there is a significant systematic error in the estimates, or (3) the two definitions are genuinely different.

To offer a further cautionary tale, in [46] the authors study a family of quasi-periodically forced $1D$ maps with an attracting invariant curve. As a parameter changes, the curve becomes less stable and it becomes increasingly wrinkled; for any resolution, parameters exist for which the curve will appear to be a strange non-chaotic attractor [35]. However, until the invariant curve completely loses stability, the attractor remains a smooth 1-dimensional curve. Hence, no matter how fine a resolution we use to compute our dimension estimate, we can never be sure that our estimate is close to its convergent value.

The Lyapunov exponents of a dynamical system describe its exponential rates of expansion and contraction, and can be used to define the Lyapunov dimension of the attractor [91]. The Kaplan-Yorke conjecture claims that the Lyapunov dimension of a generic dynamical system will equal the information dimension of its attractor [48, 27]. For a class of strange attractors arising from 2D maps, Young showed this to be the case, and moreover that the Lyapunov, Hausdorff, box-counting, and Rényi dimensions all coincide [90]. Other attractors appear to exhibit multifractal properties in the sense that different notions of fractal dimension disagree. The multifractal properties of such an attractor can be studied with a one-parameter family of dimensions such as the generalized Hausdorff [33] or Rényi dimensions[5, 67]. We can also examine this with the $\mathit{PH}_{i}$ dimensions by varying the weight parameter $\alpha$ — see Figure 12.

The Takens embedding theorem also illustrates the importance of fractal dimension to the study of attractors. In an experimental setting, it is often impossible to record all of the relevant dynamic variables. However, one can reconstruct the entire attractor from the time series of a single observed quantity using a time delay embedding [81]. Namely, for a choice of time delay $\tau$ and embedding dimension $m$, a 1-dimensional time series $\{x_{1},\dots x_{N}\}\subseteq\mathbb{R}$ may be used to construct a $m$-dimensional time series consisting of points $\{x_{i},x_{i-\tau},\dots,x_{i-(m-1)\tau}\}\in\mathbb{R}^{m}$. If $m$ is at least twice the box-counting dimension of the attractor, then generically this reconstruction will be diffeomorphic to the original attractor [75]. These techniques have been widely applied and we briefly mention some references which use this in the context of topological data analysis [29, 62, 64].

2. Dimension Estimation Methods

2.1. Persistent Homology Dimension

We give a brief, informal introduction to persistent homology for the special case of a subset of Euclidean space. For a more in depth survey, see [30, 21, 22, 14]. Note that the definition of the $\mathit{PH}$ dimension also makes sense for subsets of an arbitrary metric space.

Persistent homology [20] quantifies the shape of a geometric object in terms of how the topology of the set changes as it is thickened. To be precise, if $X$ is a subset of $\mathbb{R}^{m},$ define the family of $\epsilon$-neighborhoods $\left\{X_{\epsilon}\right\}_{\epsilon\in\mathbb{R}^{+}}$ by

[TABLE]

Figure 2 shows $\epsilon$-neighborhoods of the Sierpinski triangle $S$ and a sample of 100 points from that set. The Sierpinski triangle $S$ contains infinitely many holes which disappear as we thicken it. The first homology group of $S_{\epsilon}$ is an algebraic invariance which essentially counts the number of holes in $S_{\epsilon}.$ The last hole to disappear is in the center of the triangle; it vanishes when we have thickened the triangle by $\epsilon=\frac{1}{4\sqrt{3}}.$ Persistent homology represents this largest hole by the single interval $\left(0,\frac{1}{4\sqrt{3}}\right).$ The next three largest holes disappear at $\epsilon=\frac{1}{8\sqrt{3}}$ and correspond to three intervals $\left(0,\frac{1}{8\sqrt{3}}\right).$ The following nine holes are represented by nine intervals $\left(0,\frac{1}{16\sqrt{3}}\right),$ and so on. These are the first dimensional persistent homology intervals of the Sierpinski triangle, which we denote by $\mathit{PH}_{1}\left(S\right).$

The zero-dimensional persistent homology of a metric space tracks the connected components that merge together as the geometric object is thickened. The Sierpinski triangle is already connected at $\epsilon=0,$ so $\mathit{PH}_{0}\left(S\right)$ is trivial. However, the finite point sample x shown in Figure 2 has 100 components. As we thicken x by an amount $\epsilon$ the first component disappears when $\epsilon$ equals $\delta/2,$ where $\delta$ is the smallest pairwise distance between the points. This corresponds to the interval $\left(0,\delta/2\right)$ in $\mathit{PH}_{0}\left(\textbf{x}\right).$ One can find all $\mathit{PH}_{0}$ intervals of x by increasing $\epsilon$ and forming an interval whenever there are two points $x_{1},x_{2}\in\textbf{x}$ so that $d\left(x_{1},x_{2}\right)<\epsilon$ and $x_{1}$ and $x_{2}$ are in different components of $\textbf{x}_{\epsilon_{0}}$ for all $\epsilon_{0}<\epsilon.$ This is essentially the same as Kruskal’s algorithm for computing the minimum spanning tree on $\textbf{x},$ which leads to a proof that there is a bijection between the edges of that tree and the intervals of $\mathit{PH}_{0}\left(\textbf{x}\right)$ where an interval corresponds to an edge of twice its length.

Roughly speaking, the higher dimensional homology groups $H_{i}(X_{\epsilon})$ count the number of higher dimensional “holes” in $X_{\epsilon}.$ The higher dimensional persistent homology is defined in terms of how these groups $H_{i}\left(X_{\epsilon}\right)$ change as $\epsilon$ increases. The structure of the persistent homology is captured by a unique set of intervals that track the birth and death of homology generators as $\epsilon$ changes [12, 14]. We denote this set of intervals by $\mathit{PH}_{i}\left(X\right).$

Traditionally, the information contained in persistent homology is plotted in a persistence diagram showing the scatter plot $\left(\text{birth},\text{death}\right)$ for each interval. The persistence diagrams of the Sierpinski triangle, and point samples from it with $100$ and $1000$ points are shown in Figure 3. Note that the $1$-dimensional persistent homology of the point sample with $1000$ roughly approximates that of the Sierpinski triangle, with $1$ interval that dies around $\epsilon=\frac{1}{4\sqrt{3}},$ three that die around $\epsilon=\frac{1}{8\sqrt{3}},$ and $9$ that die around $\epsilon=\frac{1}{16\sqrt{3}}.$ This is a consequence of the “bottleneck stability” of persistent homology [15].

It is a remarkable fact that if x is a finite subset of Euclidean space, then the intervals $\mathit{PH}_{i}\left(\textbf{x}\right)$ can be computed exactly and efficiently. This is done by replacing the infinite family of $\epsilon$-neighborhoods $\left\{\textbf{x}_{\epsilon}\right\}_{\epsilon\in\mathbb{R}^{+}}$ with a finite sequence of finite simplicial complexes, called the Alpha complex of x; these complexes are subcomplexes of the Delaunay triangulation on x [20].

Note that there is more than one way to define the persistent homology of a metric space. Here, we use the persistent homology of the Čech complex persistent homology of $X.$ If $X$ is a subset of Euclidean space, this is equivalent to the persistent homology of the $\epsilon$-neighborhood filtration of $X$ described above (as well as the persistent homology of the Alpha complex, if $X$ is finite). Another common (in-equivalent) notion is the persistent homology of the Vietoris–Rips Complex on $X$ [88, 18]. We use the Čech complex because there are efficient algorithms to compute its persistent homology for finite subsets of $\mathbb{R}^{2}$ and $\mathbb{R}^{3},$ as discussed below. A persistent homology dimension defined in terms of the Vietoris—Rips complex also has nice properties [77], and may work better in examples without an embedding into a small-dimensional Euclidean space.

2.1.1. Definition of the $\mathit{PH}$ dimension

Consider a sample $\left\{x_{1},\ldots,x_{n}\right\}$ of independent points from the natural measure on the Sierpinski triangle. As $n$ increases, the persistence diagram of the point sample will converge (in the bottleneck distance) to the persistence diagram of the Sierpinski triangle itself. As such, it might seem strange that we can recover the dimension of that set from the [math]-dimensional persistent homology of the samples ($\mathit{PH}_{0}(S)$ is trivial). However, one may notice in Figure 3 that there is a large cluster of points, both [math]- and $1$-dimensional, along the plotted diagonal line. These points are generally considered to be “noise” and do not significantly contribute to the larger features of interest in other applications of persistent homology. However, the rate at which this “noise” decays is linked to the dimension of the underlying object.

To track the growth of this “noise”, we define the power-weighted sum for $\alpha>0$ by:

[TABLE]

where the sum is taken over all finite intervals and $\left|I\right|$ denotes the length of an interval. The scaling properties of random variables of the form $E^{0}_{\alpha}\left(x_{1},\ldots,x_{n}\right)$ as $n\rightarrow\infty$ has been studied extensively in probabilistic combinatorics [3, 49, 80, 93], and the case $i>0$ has recently been of interest [19, 77]. Motivated by a theorem of Steele [80] and the computational work of Adams et al. [1], Schweinhart [77] introduced the following definition of the persistent homology dimension.

Definition 1**.**

Let $X$ be a bounded subset of a metric space and $\mu$ a measure defined on $X$. For each $i\in\mathbb{N}$ and $\alpha>0$ we define the persistent homology dimension:

[TABLE]

where

[TABLE]

We write this as the $\mathit{PH}_{i}^{\alpha}$ dimension, and sometimes omit the $i$ or $\alpha$ when making general statements.

That is, $\text{dim}_{\mathit{PH}_{i}^{\alpha}}\left(\mu\right)=d$ if $E_{\alpha}^{i}\left(x_{1},\ldots,x_{n}\right)$ scales as $n^{\frac{d-\alpha}{d}}.$ Larger values of $\alpha$ give relatively more weight to large intervals than to small ones. The case $\alpha=1$ is closely related to the dimension studied by Adams et al. [1], and agrees with it if defined. Weygaert et al. [87] defined a family of minimum spanning tree dimensions that are equivalent to the $\mathit{PH}_{0}$ dimensions, and used heuristic arguments to claim that they coincide with the generalized Hausdorff dimensions for chaotic attractors. Martinez et al. [61] asserted that the $\alpha\rightarrow 0$ of the $\mathit{PH}_{0}$ dimension gives the Hausdorff dimension for point samples from chaotic attractors.

It is a corollary of Steele [80] that if $\mu$ is a non-singular measure on $\mathbb{R}^{m},$ and $0<\alpha<m$ then $\text{dim}_{\mathit{PH}_{0}^{\alpha}}\left(\mu\right)=m.$ Schweinhart [77] proved that if $\mu$ satisfies a fractal regularity hypothesis called Ahlfors regularity, then $\text{dim}_{\mathit{PH}_{0}^{\alpha}}\left(\mu\right)$ equals the Hausdorff dimension of the support of $\mu$ (which coincides with the box-counting dimension under the regularity hypothesis). He also proved that if $d$ equals the upper box-counting dimension of the support of $\mu$ and $\alpha<d$ then $\text{dim}_{\mathit{PH}_{0}^{\alpha}}\left(\mu\right)\leq d,$ as well as weaker results about the cases where $i>0.$

Note that if $\mu$ is supported on a $k$-dimensional subspace of $\mathbb{R}^{m},$ then $H_{i}\left(X\right)$ is trivial for $i\geq k.$ It follows that $\text{dim}_{\mathit{PH}_{i}^{\alpha}}\left(\mu\right)=0.$ As such, even if $\mu$ is regular, its $\mathit{PH}_{i}$ dimension may not equal its Hausdorff dimension unless the Hausdorff dimension is sufficiently large. In particular, if $\mu$ is a $d$-Ahlfors regular measure supported on a $2$-dimensional subspace of $\mathbb{R}^{m},$ $d>1.5,$ and $\alpha<d$ then $\mathit{PH}_{1}^{\alpha}$–dimension of $\mu$ equals $d.$ [77]

Schweinhart’s results show that the $\mathit{PH}_{0}^{\alpha}$ and $\mathit{PH}_{1}^{\alpha}$ dimensions of the natural measures on the Sierpinski triangle, Cantor set cross an interval, and Cantor dust equal the Hausdorff dimensions of those sets when $\alpha$ is less than the true dimension.

2.1.2. Computation of the $\mathit{PH}$ dimension

We use different methods to compute persistent homology for the cases $i=0$ and $i>0.$ In the [math]-dimensional case, we use minimum spanning tree-based algorithms that work much more quickly, and are fast even for point samples with high ambient dimension. We use the implementation of the dual–tree Boruvka Euclidean minimum spanning tree algorithm [59] included in the mlpack library [16] to compute the edges of the minimum spanning tree on a point sample. This algorithm was sufficiently fast for point clouds with $10^{6}$ points in $\mathbb{R}^{8}.$ See [59] for a brief survey of other algorithms to compute minimum spanning trees for point sets in Euclidean space and abstract metric spaces.

For $i>0,$ we use GUDHI to compute the persistent homology of the Alpha complex of the point sample. [71, 60] Both this computation and the previous one could be optimized by re-writing the data structures to support insertion of vertices.

Given a sample of $n$ points $x_{1},\ldots,x_{n},$ we compute the $\alpha$-weighted sum $E_{\alpha}^{i}\left(x_{1},\ldots,x_{c_{i}}\right)$ for $100$ logarithmically spaced values of $c_{i}$ between $1,000$ and $n.$ Then, we use linear regression to fit a power law to the data $\left(c_{i},E_{\alpha}^{0}\left(x_{1},\ldots,x_{c_{i}}\right)\right).$ After some trial-and-error, we found that fitting the power law between $x=c_{n}$ and $x=c_{\lfloor n\rfloor/2}$ provided a reasonable estimate for all examples tested. Using a smaller range sometimes produced better convergence when $i>0,$ but also introduced more oscillations in the estimate. Alternate non-linear regression methods did not seem to result in better performance.

Figure 4 shows the power law fits this method produces for a sample from the Lorenz attractor.

2.1.3. $\mathit{PH}$ complexity

In some instances, the dimension of a metric space can be computed in terms of the persistent homology of the metric space itself. For example, as shown in Figure 3, the $1$-dimensional persistent homology of the Sierpinski triangle contains intervals that scale as its dimension. This is captured by an alternate notion of $\mathit{PH}$ dimension defined by MacPherson and Schweinhart [56], which measures the complexity of the connectivity of the shape rather than the classical dimension. Here, we refer to it as the “$\mathit{PH}$ complexity” of a shape to differentiate it from notions of dimension. As we will see below, this quantity may be an indicator of when the dimension is hard to estimate using any of the methods presented here.

If $X$ is a subset of a metric space, define the cumulative $\mathit{PH}_{i}$ curve $F_{i}$ by

[TABLE]

Then the $\mathit{PH}_{i}$ complexity of $X$ is

[TABLE]

Note that $\text{comp}_{\mathit{PH}_{i}}\left(\mathbb{R}^{n}\right)=0$ for any $i.$ Also, if $S$ is the Sierpinski triangle in Figure 2, our computation of the persistence diagram of $S$ shows that $\text{comp}_{\mathit{PH}_{0}}\left(S\right)=0$ and $\text{comp}_{\mathit{PH}_{1}}\left(S\right)=\frac{\log\left(3\right)}{\log\left(2\right)}.$

We can estimate $F_{i}\left(X,\epsilon\right)$ from samples; as the Hausdorff distance between $\left\{x_{1},\ldots,x_{n}\right\}$ and $X$ converges to zero, bottleneck stability implies that $F_{i}\left(X,\epsilon\right)$ will converge for values of $\epsilon$ large relative to the Hausdorff distance. In Figure 10 below we compute $\text{comp}_{\mathit{PH}_{1}}$ of the Ikeda attractor.

2.2. The Correlation Dimension

The correlation dimension [34] is commonly used in applications because it is easy to implement and provides reasonable answers even for relatively small sample sizes. A probability measure $\mu$ on a metric space $X$ induces a probability measure $\nu$ on the distance set of $X.$ Define the correlation integral of $X$ as the cumulative density function of $\nu$:

[TABLE]

The correlation dimension equals the limit

[TABLE]

if it exists. There is an extensive literature on the estimation and properties of the correlation dimension; see for example [65, 9, 24, 84, 79].

2.3. Computation of the Correlation Dimension

For a finite point sample $x_{1},\dots,x_{n}$ one can estimate the correlation dimension as

[TABLE]

where

[TABLE]

if $n$ is taken to $\infty$ appropriately as $\epsilon\rightarrow 0.$ That is, $C\left(n,\epsilon\right)$ measures the number of distances less than $\epsilon$ in proportion to the number of all inter-point distances.

As with other dimension estimates, this limiting expression converges logarithmically slowly. To accelerate the convergence, we can essentially apply l’Hôpital’s rule to find the limit by computing a slope. To do so, we fix a collection of logarithmically spaced values $\epsilon_{1}<\dots<\epsilon_{m}$, and compute a linear regression through the data $\left(\log\epsilon_{i},\log C(n,\epsilon_{i})\right)$. Our estimate for the correlation dimension is then given by the slope of the line of best fit. For a fixed point sample, the values of $\epsilon_{1}$ and $\epsilon_{m}$ are often chosen by hand to avoid outliers and edge effects.

While computing the $\mathcal{O}(n^{2})$ inter-point distances is prohibative for very large $n$, such a calculation is not needed, as many of these distances are large and do not factor into the dimension calculation [83]. Using a $kd$-tree, one can quickly compute the $\mathcal{O}(n)$ shortest distances with $\mathcal{O}(n\log n)$ effort. After some trial-and-error, we settled upon the heuristic for choosing $\epsilon_{1}$ and $\epsilon_{m}$ as below:

[TABLE]

See Figure 4. This heuristic provides a way to choose values of $\epsilon_{1}$ and $\epsilon_{m}$ consistently amongst different fractals and different sample sizes. It appears to give near-optimal convergence rates to the true dimension for self-similar fractals, and estimates that agree with the values in the literature for the Hénon, Ikeda, and Lorenz attractors [34, 79, 83].

As has been previously reported [83], the “statistical test error” from the linear regression calculation does not have much predictive value about the limiting dimension — see Figure 5. Also, it is much smaller than the empirical standard deviation of the dimension estimate between trials in some cases.

2.4. The Box-counting Dimension

The box-counting dimension [10] of a bounded subset of $X$ of $\mathbb{R}^{m}$ is defined in terms of the number of cubes of width $\delta$ needed to cover $X.$ Let $\left\{C_{i}^{\delta}\right\}_{i\in\mathbb{N}}$ be the cubes in the standard tiling of $\mathbb{R}^{m}$ by cubes of width $\delta,$ and let $N_{\delta}\left(X\right)$ be the number of cubes in $\left\{C_{i}^{\delta}\right\}_{i\in\mathbb{N}}$ that intersect $X.$ Define the upper and lower box-counting dimensions by

[TABLE]

respectively. If the upper and lower box-counting dimensions coincide, the shared value is called the box-counting dimension of $X$ and is denoted $\text{dim}_{\text{box}}\left(X\right).$ There are several equivalent definitions; see Falconer [25] for details. Many studies have investigated the properties and estimation of the box-counting dimension, including [74, 50, 82].

2.4.1. Computation of the Box-counting Dimension

We found it difficult to find a general method to compute the box-counting dimension that worked well for different examples and different numbers of samples. For example, computing box-counts of the form $N_{2^{-i}}\left(\textbf{x}\right)$ and fitting a power law between $2^{-i}$ and $2^{-j}$ resulted in estimates that were sensitive to $i,$ $j,$ and the number of samples. In some cases, it was easy to cherry-pick a specific choice just because it seemed to have the best convergence to the true dimension, though other choices resulted in power law fits that looked just as good. See Figure 6 for an example, which plots the error in the dimension estimate for many possible power law fits. Only very specific choices have errors as small as estimates of the correlation and $\mathit{PH}_{0}$ dimensions for the same sample.

We settled on the following method, which produced reasonably good results for planar examples. It is based on the observation that if $\left\{x_{j}\right\}_{j\in\mathbb{N}}$ are samples from $X\subset\mathbb{R}^{m}$ and $\delta>0$ then the box count $N_{\delta}\left(x_{1},\ldots,x_{m}\right)$ should converge to $N_{\delta}\left(X\right)$ as $n\rightarrow\infty.$ For a sample $\left\{x_{1},\ldots,x_{m}\right\}$, we estimate the box-counting dimension from the smallest boxes for which $N_{\delta}\left(x_{1},\ldots,x_{m}\right)$ appears to have stabilized.

For a a family of point samples $\textbf{x}_{m}$ in $\mathbb{R}^{n}$ (where the sizes of $\textbf{x}_{m}$ are logarithmically spaced as before) we rescale and translate the point samples so they are contained in a unit cube. Then, we compute box-counts of the form $N_{i/10000}\left(\textbf{x}_{m}\right)$ for $1\leq i\leq 1000.$ We fit a power law to the data $\left(i/10000,N_{i/10000}\left(\textbf{x}_{m}\right)\right)$ in the range $(\left\lceil{j/2}\right\rceil,j),$ where $j$ is the smallest index so that

[TABLE]

We used linear regression to fit the power law; non-linear regression did not produce substantially different results. See Figure 4.

We tried several variations. For example fitting the power law in the range $(\left\lceil{.9j}\right\rceil,j)$ produced estimates that converged faster with $n$ for some examples but exhibited large oscillations in others.

3. Results for Self-similar Fractals

We compare the performance of the fractal dimension estimation procedures for four different self-similar fractals: the Sierpinski triangle ($S$), the Cantor dust ($C\times C$), the Cantor set cross an interval ($C\times I$), and the Menger sponge ($M$). The first three are subsets of $\mathbb{R}^{2}$, and $M$ is contained in $\mathbb{R}^{3}.$ We chose these sets to illustrate the observed relationship between $\text{comp}_{\mathit{PH}_{0}}\left(X\right)$ and the performance of dimension estimation techniques; see Table 1. Definitions of the sets and sampling methods are contained in Appendix A.

For each set, we sample $50$ trials of $10^{6}$ points from the corresponding natural measure. We compute 100 dimension estimates for each trial, at logarithmically spaced numbers of points between $10^{3}$ and $10^{6}.$ To show the convergence of the dimension estimate to the true value, we plot number of points against the dimension estimate averaged across the trials, with thinner dotted lines one standard deviation above and below the estimate.

Dimension estimates for the four examples are plotted in Figure 7. In all examples, the $\mathit{PH}_{0}$ dimension and the correlation dimension perform better than the $\mathit{PH}_{1}$ or box-counting dimensions. The relatively poor performance of the $\mathit{PH}_{1}$ dimension can be ascribed to the fact that the weighted sums $E_{\alpha}^{1}\left(x_{1},\ldots,x_{n}\right)$ are smaller and noisier than the corresponding sums $E_{\alpha}^{0}\left(x_{1},\ldots,x_{n}\right)$ (see Figure 4). As mentioned previously, it was difficult to find an effective general method to produce box-counting estimates. This is illustrated for the Sierpinski triangle in Figure 6. While the plot of box width vs. box count looks linear for a large range on the log–log plot, only very specific choices of bounds will produce dimension estimates with errors on the same order as the $\mathit{PH}_{0}$ and correlation dimensions. These bounds vary unpredictably with the sample size and example.

The convergence of the $\mathit{PH}_{0}$ and correlation dimensions to the true dimension is very similar in all four examples. The correlation dimension estimates exhibit oscillations as the sample size increases, a phenomenon often ascribed to the fractal’s lacunarity [83]. The $\mathit{PH}_{0}$ dimension estimate appears less sensitive to oscillations, but this is likely due to how the power law fit is performed in the computations (as described in Section 2.1.2); choosing a narrower range over which to fit a power law results in a more oscillatory estimate. That said, for sizable samples, it is impractical to compute the correlation integral over a large range. The $\mathit{PH}_{0}$ dimension provides a computationally practicable way to use information from multiple lengthscales. It also has the advantage of having the parameter $\alpha$ which can be tuned to give a better estimate — see Section 3.2 for a discussion on the choice of $\alpha.$ The variance of the correlation dimension estimate is slightly lower, but this doesn’t mean much when oscillations are present.

3.1. Percentage Error

An interesting pattern emerges when we compare the percentage error of the dimension estimates across the three planar examples, as in Figure 8. Both $\mathit{PH}_{0}$ and correlation dimensions perform best for the Sierpinski triangle ($\text{comp}_{\mathit{PH}_{0}}(S)=0$), second best for the Cantor set cross an interval ($\text{comp}_{\mathit{PH}_{0}}(C\times I)=\text{dim}\left(C\times I\right)-1=\frac{\log\left(2\right)}{\log\left(3\right)}),$ and worst for the Cantor dust ($\text{comp}_{\mathit{PH}_{0}}(C\times C)=\text{dim}\left(C\times C\right)=\frac{2\log\left(2\right)}{\log\left(3\right)}$). That is, dimension estimation appears to be more difficult when the connectivity of the underlying set is more complex. Note that $\text{comp}_{\mathit{PH}_{0}}(X)$ can be computed (see [56]), and could be used as an indicator of whether additional caution is warranted when discussing dimension estimation results.

For the $\mathit{PH}_{1}$ dimension, the situation is different and the difficulty of dimension estimation appears to (unsurprisingly) be related to $\text{comp}_{\mathit{PH}_{1}}$ rather than $\text{comp}_{\mathit{PH}_{0}}.$ The rate of convergence was fastest for the Cantor set cross an interval ($\text{comp}_{\mathit{PH}_{1}}(C\times I)=0),$ but slower for the Cantor dust ($\text{comp}_{\mathit{PH}_{1}}(C\times C)=\text{dim}\left(C\times C\right)=\frac{2\log\left(2\right)}{\log\left(3\right)}$) and Sierpinski triangle ($\text{comp}_{\mathit{PH}_{1}}(S)=\text{dim}\left(S\right)=\frac{\log\left(3\right)}{\log\left(2\right)}$).

We exclude the Menger sponge from these figures, as the ambient dimension likely influences the difficulty of dimension estimation.

3.2. Dependence on $\alpha$

Estimates of the $\mathit{PH}_{0}^{\alpha}$ and $\mathit{PH}_{1}^{\alpha}$ dimensions for various choices of $\alpha$ are shown in Figures 9. In the cases where $\text{comp}_{\mathit{PH}_{i}}\left(X\right)=0,$ there are only small differences between dimension estimates for different choices of $\alpha,$ with perhaps a slight advantage for higher values of $\alpha$ (Figures 9 and 9). However, when $\text{comp}_{\mathit{PH}_{i}}\left(X\right)>0,$ dimension estimates for different values of $\alpha$ are substantially different. Lower values of $\alpha$ yield better estimates when $\text{comp}_{\mathit{PH}_{i}}(X)>0$ for planar examples, and middle values of $\alpha$ (i.e. equal to about half the true dimension) seem to provide the best convergence for $i=1,2$ (but convergence is slow, especially for $i=2.$) Lower values of $\alpha$ appear to give dimension estimates that have a higher variance between samples — see Figures 9 and 9 for the Sierpinski triangle, which is representative.

4. Results for Attractors

For each chaotic attractor we sample $50$ randomly selected initial conditions, and generate a time series of $10^{6}$ points after discarding an initial transient trajectory (except the Mackey-Glass attractor, where we sample $25$ initial conditions). For each time series, we compute 100 dimension estimates for each trial, at logarithmically spaced numbers of points between $10^{3}$ and $10^{6}.$ The various dynamical systems we studied are described in Appendix B and the dimension estimates we obtained for $N=10^{6}$ are summarized in Table 3 and Figure 11.

In all cases, we are skeptical that the box-counting dimension has approached a limiting value. The upper box dimension is known to be an upper bound for the $\mathit{PH}_{0}^{\alpha}$ [77] and correlation dimensions, but most of our estimates of these dimensions tend to be higher than our estimates of the box-counting dimension. This comports with previous observations that box-counting estimates do not seem to converge well for strange attractors [36].

Of the two-dimensional maps we studied, the Ikeda map appears to be most well-behaved. The box-counting, $\mathit{PH}_{0}$, and even the $\mathit{PH}_{1}$ dimension all seem to converge to a common value of $1.71$, with a correlation dimension of $1.68$ not that far off. The Ikeda attractor $I$ was the only chaotic attractor we studied with non-trivial persistent homology at multiple lengthscales (rather than the persistent homology “noise” of point samples). Here, we compute its $\mathit{PH}_{1}$ complexity, $\text{comp}_{\mathit{PH}_{1}}\left(I\right).$ Figure 10 shows the cumulative interval count function defined in 2.1.3, for a sample of $10^{6}$ points from the Ikeda attractor. The $\mathit{PH}_{0}$ curve only shows very small intervals, which come from the noise of point samples. This is as expected, because the Ikeda attractor is connected. However, the $\mathit{PH}_{1}$ curve shows two different regimes, one that corresponds to noise and one that appears to be picking up on the elongated holes visible in Figure 1. The cumulative length plot appears to follow a power law in the range $\epsilon=.0005$ to $\epsilon=.005.$ Fitting a power law to this range gives an estimate of $\text{comp}_{\mathit{PH}_{1}}\left(I\right)=.95.$ This is lower than the dimension estimates for the Ikeda attractor, indicating that the elongated holes in that set scale at a different rate than its other properties.

For the Hénon attractor, both $\mathit{PH}_{0}$ and $\mathit{PH}_{1}$ dimension estimates had large oscillations which make it hard to say whether they will converge to values different from either the box-counting or correlation dimensions. The correlation dimension exhibited smaller oscillations, which comports with previous observations by Theiler [85], who also warns of the danger of estimating the dimension of attractors with long-period oscillations. Previous studies [4] have claimed that the Hénon attractor is multifractal for the parameters chosen here, but we are not confident enough in our dimension estimates to make such an assertion.

The Rulkov attractor is non-homogeneous, as easily evidenced by Figure 1, and not surprisingly our fractal dimension estimates differed from each other, ranging from $1.00$ (correlation dimension) to $2.13$ ($\mathit{PH}_{1}^{1}$ dimension). The $\mathit{PH}_{1}^{\alpha}$ estimates began above $2$ for small samples sizes and decreased toward $2$ as the number of samples increased. For various values of $\alpha,$ the $\mathit{PH}_{0}^{\alpha}$ dimension estimates ranged from $1.38$ for $\alpha=1$ to $1.88$ for $\alpha=1.25$ (see Figure 12). As one may see in Figure 11, the variance of the correlation dimension estimate increases with sample size, suggesting the heuristic we chose for fitting the correlation integral requires further fine-tuning for this example. The Rulkov map is noninvertible, so it does not belong to the class of 2D maps considered in [90] where many of the fractal dimensions are known to coincide. While it is difficult to determine if numerical dimension estimates have converged, our results certainly suggest that the various fractal dimension definitions may not agree for the Rulkov attractor.

For the Lorenz attractor, the $\mathit{PH}_{0}$ and the correlation dimensions performed well and both appear to converge toward $2.05$. However, neither the $\mathit{PH}_{1}$ nor the $\mathit{PH}_{2}$ dimensions perform well. In fact, nearly $10^{5}$ are needed before the $\mathit{PH}_{2}$ dimension estimate less than 3, the ambient dimension in which the points reside! As in other cases, this is likely due to the fact that higher dimensional persistent homology contains less information — see Figure 4. The sum of the lengths of the $\mathit{PH}_{1}$ and $\mathit{PH}_{2}$ intervals of a point sample of $10^{6}$ points are $10$ and $400$ times smaller than the sum of the lengths of the $\mathit{PH}_{0}$ intervals, respectively. Also, the plot of $n$ versus $E_{1}^{i}\left(x_{1},\ldots,x_{n}\right)$ is much noisier for $\mathit{PH}_{2}$ than $\mathit{PH}_{0}.$

Lastly we consider the Mackey-Glass equation, a delay differential equation for which it is still an open conjecture whether there exist parameters for which the system exhibits mathematically provable chaos. The phase space for this system is infinite dimensional, and we use a projection into $\mathbb{R}^{8}$ for our fractal dimension calculations. In [26] Farmer reports the attractor’s Lyapunov dimension to be $3.58$, and we obtain a $\mathit{PH}_{0}^{.5}$ dimension estimate close to this value. However, the spectrum of $\mathit{PH}_{0}^{\alpha}$ dimension for various $\alpha$ had the largest spread of any of the fractals we studied, except the Rulkov attractor, ranging between $3.59$ for $\alpha=.5$ and $4.06$ for $\alpha=3.5$ (see Figure 12). While none of our dimension estimates appear to oscillate, all appear to converge slowly. The box-counting dimension estimation method described in Section 2.4.1 performed particularly poorly for this example, perhaps due to the high dimension and co-dimension. The estimate in Table 3 was computed by fitting a power law by hand. We did not attempt to compute either the $\mathit{PH}_{i}$ dimensions for $1\leq i\leq 7,$ because the algorithm we use is impractical if the ambient dimension is greater than $3.$

5. Earthquake Data

We estimate the dimension of the Hauksson–Shearer Waveform Relocated Southern California earthquake catalog [40, 51]. The hypocenter of an earthquake is the location beneath the earth’s surface where the earthquake originates. One can use the dimension estimate to study the geometry of earthquakes, i.e. by comparison to dimension estimates for fractures in rock surfaces or dislocations in crystals (see Section 6 of [47]). Previously, Harte [39] estimated the correlation dimension of earthquake hypocenters in New Zealand and Japan, and Kagan [47] studied the correlation dimension of earthquakes in southern California and developed extensive heuristics to correct for errors in that computation.

We downloaded the coordinates of the hypocenters of the earthquakes in the catalogue from the Southern California Earthquake Data Center Website [76], and selected the waveform-relocated earthquakes of magnitude greater than 2.0. This gave a data set of 74,929 earthquakes. The deepest earthquake in the data set is $33$ kilometers below sea level and the depth distribution has a single peak, so we did not split the data set into multiple samples as Harte did. The data set had small relocation errors, with the error in distance between nearby hypocenters estimated to be less than $.1$ kilometer for $90\%$ of such distances [40].

The presence of location errors means that a different methodology is required to estimate the correlation dimension than from the previously considered examples, where estimates at shorter lengthscales produce better results. To estimate the correlation dimension we plotted the correlation integral (Figure 13) and determined that the best range to fit a power law was between $2$ and $8.5$ kilometers. This yields a correlation estimate of about $1.66.$ This can be compared to Kagan’s estimate of 1.5 for earthquakes from an earlier version of the same earthquake catalog that was processed with different methodology (Figure 8 of [47]), and to Harte’s estimates of 1.8 and 1.5 for shallow earthquakes in Kanto, Japan and Wellington, New Zealand, respectively.

Figure 13 shows the plot of $n$ versus $E_{0}^{1}\left(x_{1},\ldots,x_{n}\right)$ used to estimate the $\mathit{PH}_{0}$ dimension. We estimated that the $\mathit{PH}_{0}^{1}$ dimension is approximately $1.76$ when $\alpha=1.$ Similar computations for $\alpha=.5$ and $\alpha=1.5$ yielded dimension estimates of $1.75$ and $1.83,$ respectively. From this, and the comparison with the correlation dimension, we have evidence that the earthquake probability distribution is not regular. Also, note that the $\mathit{PH}_{0}^{1}$ dimension plot (Figure 13) appears to follow a power law in a large range, and the dimension can be estimated without fiddling with example-specific parameters. This precludes difficulties that arose in Harte’s analysis where the correlation dimension estimate was very sensitive to the scale at which it was measured(Figure 20 in [39]). Of course, the correlation dimension has an advantage of interpretability at different lengthscales — the $x$-axis in Figure 13 is kilometers, while it is number of samples in Figure 13.

The box-counting, $\mathit{PH}_{1},$ and $\mathit{PH}_{2}$ dimensions did not produce good results. In the former case, it is unclear whether the plot of box width versus box count (Figure 13) follows a power law at any range, but certainly no range shorter than $5$ kilometers. Fitting between $5$ and $10$ kilometers yielded an estimate of $.91$, which is very different than the estimates computed with other methods. Fitting at larger length scales, between $20$ and $50$ kilometers produced an estimate of $1.22.$ It is not surprising that the $\mathit{PH}_{1}$ and $\mathit{PH}_{2}$ dimensions produced poor results, given the small sample size.

We also estimated that $\text{comp}_{\mathit{PH}_{i}}(X)$ likely equals zero for the earthquake data, for $i=0,1,2.$

6. Conclusion

Overall the performance of our $\mathit{PH}_{0}$ dimension calculations is comparable to the correlation dimension, and gives dimension estimates that are often greater for non-regular sets. Both the $\mathit{PH}_{0}$ and correlation dimensions were more reliable than either the box-counting, $\mathit{PH}_{1},$ or $\mathit{PH}_{2}$ dimensions. The correlation dimension and the $\mathit{PH}_{0}$ dimension are easy to implement and quick to compute. The $\mathit{PH}_{0}$ dimension may be more user-friendly in the sense that a single choice of power law range worked well for all examples. Of these two, we do not view one dimension estimation technique as “better” than the other. Rather, they are complementary methods that — if the underlying set is not regular — can provide complementary information.

In many applications of persistent homology, small intervals are often ignored in favor of larger features and discarded as “noise”. However, as we have demonstrated here, there is sometimes signal in that noise. Adding the $\mathit{PH}$ dimension to the current suite of $\mathit{PH}$-based techniques used in applications such as machine learning will offer an orthogonal data descriptor.

Appendix A Self-Similar Fractals

A.1. The Sierpinski Triangle

The Sierpinski triangle (Figure 1) is defined by iteratively removing equilateral triangles from a larger equilateral triangle. Begin by sub-dividing the equilateral triangle formed by $(0,0),(1,0),$ and $(\sqrt{3}/2,1/2)$ into four congruent triangles, and removing the center triangle. Repeat this process on the remaining three triangles, and continue ad infinitum. The Sierpinski triangle equals three copies of itself rescaled by the factor $1/2$ so the dimension of the resulting set is $\frac{\log\left(3\right)}{\log\left(2\right)}$ (for a self-similar set equal to $m$ copies of itself rescaled by a factor $r,$ the self-similarity dimension equals $\frac{\log\left(m\right)}{\log\left(1/r\right)}.$

We sampled points from the natural measure on the Sierpinski Triangle by sampling random integers $\left\{a_{1},a_{2},\ldots\right\}\in\left\{0,1,2\right\}$ and computing

[TABLE]

In practice, we end this procedure at $i=64.$

A.2. The Cantor Dust

The standard middle-thirds Cantor set $C$ is the set formed by removing the interval $\left(1/3,2/3\right)$ from the closed interval $\left[0,1\right],$ and iteratively removing the middle third of the remaining intervals. The dimension of $C$ is $\frac{\log\left(2\right)}{\log\left(3\right)}.$ The Cantor dust (Figure 1) is the product $C\times C;$ its dimension is $\frac{2\log\left(2\right)}{\log\left(3\right)}.$

We sampled points from the natural measure on $C$ by sampling random integers $\left\{a_{1},a_{2},\ldots\right\}\in\left\{0,1\right\}$ and computing

[TABLE]

In practice, we truncated the summation at $i=64.$ A point $(x_{1},x_{2})$ from the natural measure on the Cantor dust can then be sampled by independently sampling $x_{1}$ and $x_{2}$ as above.

A.3. The Cantor Set Cross an Interval

Consider the set $C\times\left[0,1\right],$ where $C$ is the Cantor set defined in the previous section. The dimension of this set is $1+\frac{\log\left(2\right)}{\log\left(3\right)},$ one greater than the dimension of the Cantor set.

We sampled points $(x,y)$ from the natural measure on the Cantor set cross an interval by sampling a random point $x$ on the Cantor set by the procedure described in the previous section, and a random real number $y$ from the uniform distribution on $\left[0,1\right].$

A.4. The Menger Sponge

The Menger sponge is defined by iteratively removing cubes from a larger cube. Begin by sub-dividing the unit cube in $\mathbb{R}^{3}$ into nine smaller sub-cubes, and remove nine of them: one from the center of the original cube one from the center of each of the eight faces. Repeat this process on each of the $20$ remaining sub-cubes and continue ad infinitum. The resulting structure is equal to $20$ copies of itself rescaled by a factor $1/3$ so the dimension of the resulting set is $\frac{\log\left(3\right)}{\log\left(20\right)}.$

For $j\in\mathbb{n}$ we sample a uniform three-tuple of integers $\left(x_{i},y_{i},z_{i}\right)\in\left\{0,1,2\right\}$ throwing out and re-sampling any three-tuple for which two or more of the coordinates equals one. To sample a point from the natural measure on the Menger sponge, we form the sum

[TABLE]

In practice, we end this procedure at $i=64.$

Appendix B Chaotic Attractors

B.1. Hénon Map

The Hénon map [42] is given by $(x_{n},y_{n})\mapsto(x_{n+1},y_{n+1})$ where:

[TABLE]

We used the parameters of $a=1.4$ and $b=0.3$. We generated $50$ time series using randomly chosen initial conditions, computing a trajectory of length $1.1\cdot 10^{7}$ and discarding the initial $10^{6}$ points in the series.

B.2. Ikeda Map

The complex Ikeda map [38] is given by $z_{n}\mapsto z_{n+1}$ where:

[TABLE]

We used the parameters:

[TABLE]

We generated 50 time series using randomly chosen initial conditions, computing a trajectory of length $1.1\cdot 10^{6}$ and discarding the initial $10^{5}$ points in the series.

B.3. Rulkov Map

The chaotic Rulkov map [72, 45] is given by $(x_{n},y_{n})\mapsto(x_{n+1},y_{n+1})$ where:

[TABLE]

We use the parameters

[TABLE]

We generated 50 time series using randomly chosen initial conditions, computing a trajectory of length $1.1\cdot 10^{6}$ and discarding the initial $10^{5}$ points in the series.

B.4. Lorenz System

The Lorenz system [53] given by the system of ordinary differential equations below:

[TABLE]

We use the parameters:

[TABLE]

We generated 50 time series using randomly chosen initial conditions, integrating forward for $1.01\times 10^{5}$ units of time, and discarding the initial $10^{3}$ units of time. We used MATLAB’s ode45 integrator with a relative error tolerance of $10^{-9}$ and an absolute error tolerance of $10^{-9}$. We sampled our trajectories at a rate of $10$Hz, producing time series with $10^{6}$ points each.

B.5. Mackey-Glass

The Mackey-Glass equation [55] is given by the scalar delay differential equation below:

[TABLE]

We take parameters:

[TABLE]

A natural phase space for this dynamical system is the Banach space of continuous functions $C=C^{0}([-\tau,0],\mathbb{R})$. For our dimension calculations we discretize this space using a projection map $\pi_{m}:C\to\mathbb{R}^{m}$ (for $m=8$) which evaluates a function $y\in C$ at $m$ maximally spaced points on the interval $[-\tau,0]$.

We generated 25 time series using randomly chosen initial conditions and integrating forward for $\tfrac{\tau}{m-1}\cdot 10^{6}+10^{4}$ units of time. We used MATLAB’s dde23 integrator with the default error tolerances for the first $5,000$ units of time, and the remainders of each trajectory were computed using a relative error tolerance of $10^{-5}$ and an absolute error tolerance of $10^{-9}$. The initial transient periods of $10^{4}$ units of time were discarded, and the remainders were sampled at a rate of $\tfrac{\tau}{m-1}$Hz, producing time series with $10^{6}$ points each.