The minimal spanning tree as an estimator for generalized dimensions

PhysicsLettersA 169 (1992) 145—150 North-Holland

PHYSICS LETTERS A

The minimal spanning tree as an estimator for generalized dimensions Rien van de Weygaert a, Bernard J.T. Jones a b

a,b

and Vincent J. MartInez

Sterrewacht Leiden, Postbus 9513, 2300 RA Leiden, The Netherlands Niels Bohr Institutet, Blegdanisvej 17, DK-2100 Copenhagen, Denmark Departament deMatemàtica Aplicada i Astronomia, Universitat de València, Burjassot, 46100 Valencia, Spain

Received 29 October 1990; revised manuscript received 28 January 1992; accepted for publication 23 July 1992 Communicated by A.P. Fordy

Discrete random point sets can be characterized by their dimensionality. We show that the minimal spanning tree of the set provides a good estimator for the dimensionality. The method is shown to work by testing it on the Hénon attractor; advantages and disadvantages to other, more conventional, methods are discussed.

1. Introduction Describing point sets is a long standing problem. The aspect of point clustering we wish to discuss is that of the generalized dimension of the set. This is a simple extension of the concept of Hausdorff dimension introduced at the beginning of the century by Hausdorif [1,2]. A major problem in determining the dimensionality of point sets, especially in the case of experimental or observational data, is the fact that we only have a limited number of points in our set, so that it becomes difficult to get reliable estimates of generalized dimensions of these point sets from conventional methods. In this paper we wish to introduce the minimal spanning tree (MST) as a new, natural, tool for estimating a range of dimensions from a set consisting of a finite number of points,

2. The minimal spanning tree The MST is a graph-theoretical construct that was introduced by Kruskal [3] and Prim [4]. The mmimal spanning tree of a set of N points is the unique network of N— 1 edges (each linking two points) providing a route between any pair of points while minimizing the sum of the lengths of the edges. It

defines a set of disks whose diameters are the edge lengths, and which when taken together provide the minimal covering in the sense of the partition function definition (eqs. (1), (2)), so that it comes as close as possible to realizing the definition of the Hausdorif dimension for a discrete point set. Examples of two-dimensional minimal spanning trees are shown in fig. 1. Formally stated in graph-theoretical terms the minimal spanning tree can be formulated as follows: The data set is a graph G, consisting of a vertex set V (the points) and edge set E (an edge is a straight line connecting two points, E is a subset of V XV), each having a “length” or “weight” (in our case the Euclidian distance between the two vertices). A sequence of edges joining vertices is a path, a closed path is called a circuit. If there is a path between any pair of vertices the graph is called connected. A connected graph containing no circuits is called a tree. Ifthe tree of a connected graph contains all the vertices ofthe vertex set then it is called a spanning tree. The length of a tree is defined to be the sum of the weights of the component edges. The minimal spanning tree (MST) is the spanning tree of minimal length. The algorithm [4,5] we used for calculating the MST of a given point set starts by choosing an arbitrary point of the set and finding its nearest neigh-

0375-960 1/92/$ 05.00 © 1992 Elsevier Science Publishers B.V. All rights reserved.

145

Volume 169, number 3

PHYSICS LETTERS A

~ / ~ /

in the subtree the updating requires computation proportional to N—M, for Mgoing from ito N— 1, so that the total computation time is proportional to IN(N— 1). Moreover, the computation time is just of the point set. More efficient, dNlog N, algorithms linearly proportional to d, the embedding dimension are in principle possible, but not yet available. When looking at the behaviour of the minimal spanning tree in a Poissonian point set (fig. 1) in two dimensions we noticed its close link to the dimensionality ofthe point set, filling the part ofspace occupied by the points uniformly at every level of

V -

/ 7~ j-

__________

2

~ ‘\~(i~

4ç2. 5

~

~

~

~ _\~cK.~-

1

~

(~~<~S r/. ~<.,L’1

~

j~

rJ ~

~

J

~

~

~ ______ ______

—

________________ ~ ~‘

~‘y~’

I

~

5 j’~

I

~ ~

-

~

‘~--~

\J?-~ ~t

Fig. 1. Three Poisson point sets on the left, with increasing number density from top to bottom, and theircorresponding minimal spanning trees on the right.

bour. These two points and the corresponding edge form the subtree T1. For each isolated point (point not yet in the subtree) the identity and distance to its nearest neighbour within the subtree is stored; by definition this distance is called the distance of the isolated point to the subtree. These potential MST edges are called links. The Mth subtree, TM, is formed by adding to TM_I the nearest neighbour of the subtree, being the isolated point whose distance to the subtree is minimal, together with the corresponding link. After adding this new point to the subtree the new links are updated, i.e. the distance from each isolated point to the new subtree vertex is calculated to see whether it is smaller than the previous distance of the point to the subtree. The resulting TN_i is the minimal spanning tree. If there are M points 146

21 September 1992

resolution. In fact, this close resemblance had been noticed earlier, independently, by Steele [6] who suggested “that the asymptotic behaviour of the minimal spanning tree might be a useful adjunct to the concept of dimension in the modeling applications and analysis of fractals”. The MST construct has been used with some success in astronomy as a tool to identify filamentary structures in the distnbution of galaxies [7] In two earlier papers [8] we applied minimal spanning trees to the estimation of dimensions in the galaxy distri bution In particular, they gave a good estimation of the Hausdorif dimension of the galaxy distribution from a sample of a few hundred galaxies within an error of 0.1. 3. Generalized dimensions To understand the usefulness ofthe minimal spanning tree in a fractal context we first have to emphasize both the “infimum” in the partition sum, as well as the difference between the Hausdorif generalized dimensions and the Rényi dimensions. The generalizations of the Hausdorif dimension [1,21 of a set d embedded in Fl!’ can be defined by considering, for all r> 0, the family of coverings of d, T formed by sets B with diameters r, < r. Assigning a measure jL, to each subset and taking the moments of this measure we can define the partition function, ~,

q

F( q, r) = lim inf ~ r-.O

)~

if r~0 and q

.

F(q,r)=limsup~— r-.O

T~,.

~

1

(1)

,

,

r~

if

~0andq~1

.

(2)

Volume 169, number 3

PHYSICS LETTERS A

By demanding that f’(q, t) =const we get the function r( q). It is obvious that for q= 0, this demand and eq. (1) lead to the definition of the Hausdorff dimension, D~= x( 0). In a likewise manner we can define the Hausdorif generalized dimensions [9] from the function r( q), —

(3)

D(q)=(q—l)~’r(q).

A closely related, but somewhat different, quantity as the Hausdorif dimension is the capacity [10] of a set d c Fl”. If N( r) is the minimum number of ncubes of side r that we need to cover d, then the capacity of d is D~(~h1)= ~

logN(r) log(1 /r)~

(4)

Notice that the difference with the Hausdorif dimension is that now all the domains of the coverings are sets with the same diameter, The natural generalization of the concept of capacity are the R~nyidimensions [11] or generalized dimensions [121. Consider the measure p and a probability measure p ( r) corresponding to the coyering of d formed by disjoint cells of side r, {C1(r)}~c~, p(C~(r))=p1(r). Then for all qeR the Rényi dimensions are defined as the generalizations of the capacity (eq. (4)), ~ r-.O

log r

In the same way that D( q) for q= 0 is the Hausdorff dimension, Dq for q= 0 is the capacity D~,and it is possible to show [9,131 that D(q) ~Dq. It is important to appreciate the differences between the measures D4 and D(q) because the mmimal spanning tree formalism, as explained below, estimates the Hausdorifgeneralized dimensions while dimension estimators like box counting, correlation method and k nearest neighbours try to fit Dq, i.e. the Rényi generalized dimensions.

4. The MST estimator technique As is seen in eqs. (1) and (2), the partition sum leading to the generalized dimensions contains an “infimum” term (q~ 1), or a “supremum” term in

21 September 1992

the case q>~1. This means that one should try to make the partition sum as small as possible. In most methods developed to extract the dimensions of a point set this “infimum” (or “supremum”) term is not taken into account. However, when using the edge lengths of a minimal spanning tree, as the diameters of disks covering the set, one is quite naturally led to an extremum of the sum amongst all possible coyerings of a finite point set. The minimal spanning tree is the unique minimal covering of a finite point set. This means it is the natural generalization ofthe concept of an infimum covering of a set to a finite point set. Analyzing the scaling of the distribution of edge lengths of minimal spanning trees of finite point realizations of a process leads one directly to an estimate of the generalized dimensions. The minimal spanning tree with m edges produces a set of edge lengths {l~ } with m = NR 1; NR is the number of points in a subsample formed by random selection from the whole sample. In this case we define a partition function by

r=

S(r, m)=

—

m

~,

~ 11(m)’.

ji

—

(6)

If this sum scales as a power law of the number of edges m in the subsample we have S(t,m)=m~_i~const, (7) with the constant depending only on r. By fitting this scaling relation for different values of m (m+ 1 points selected randomly from the whole sample) we obtain q as a function of r. Notice that m has to be large enough for the branches of the corresponding MST to be small enough to be within the scaling region. This scaling relation is a generalization of the one analytically derived by Steele [6] for Poisson point processes in Euclidian spaces. The results can be somewhat unstable due to finite sampling effects caused by the random selection of m + 1 points. In order to surmount this problem we perform, for a given number of points m +1, the random selection of m + 1 points M times, and calculate for each of these M realizations the MST. From the set of these realizations we calculate q as a function of r by fitting the relationship

(q— 1) log m.

(8) 147

Volume 169, number 3

PHYSICS LETTERS A

The brackets denote the average taken over the M realizations. In practice we took M= 10. An additional advantage of this procedure is that it also reduces the intrinsic oscillations common to different scaling methods [14]. At both small and large edge lengths, where the l//~and 1, distributions have long and low tails one can expect “Poisson sampling noise” [151. This should make one especially careful in the case ofthe estimation of D(q)’s with q.<<0, because very long edges might distort the partition sum [16]. Careful study of the behaviour ofthe tails of the distribution function of the MST branch lengths as a function of m showed that for small values of m (m ~ 1000) the length of the largest branch of the tree is usually less than 10 times the mean length of all the branches, T( m), while the shortest branch is longer than 0.01 times T(m). However, the tails of the distribution get longer as m increases. The most reliable way of dealing with this problem is to take out these rarified tails for large m. For example, in our test (see below), we do not consider branches smaller than 0.OlT(m) or large than 10T(m). By applying these cut-offs nearly all the branches forming the MSTs are considered for m ~ 1000. For larger values of m only very few of the smallest and largest edges are removed. On average less than 0.09% of all the edges in the trees are taken out. By means ofthis procedure we have entered the region of the real distribution function.

5. The MST in practice: the Hénon attractor

21 September 1992 —_______________________

-

-. .-~.

_ . - -.

-. -

-----—-----~---—-

.

.

‘_~...

.

- -

_______________________

-SS~_

•-~~

,___.______,)

-. _____________

_________________

Fig. 2. A sample of 1000 points from the Hénon attractor (top) and the corresponding minimal spanning tree (bottom).

__________________________________ 8

\~

‘

I

I

‘• • 0 ~

.

a 0.5

1

1.5

~ 0

..~

+

.

* 4

~

—5

0

5

I

10

—T(q)

The use of the minimal spanning tree was extensively tested on the well-known Hénon attractor (a = 1.4, b= 0.3) [17]. The first 1000 points of the Hénon attractor are shown in fig. 2 (top). The bottorn plot offig. 2 depicts the resulting minimal spanning tree of the 1000 point Hénon attractor. Clearly the minimal spanning tree follows the attractor extremely well, In our test we used 20000 consecutive points of the Hénon attractor and took random samples from 100 to 7500 points. Minimal spanning trees were then calculated at nine different logarithmically equidistant values of m+ 1 points. Note that the number of data points required by the MST method presented 148

S

Fig. 3. A plot of q against — r(q) determined with the minimal S spanmng tree method for the Hénon attractor. The solid line shows the behaviour in the case ofa simple fractal. The inset shows the resultingf(a) spectrum.

here is considerably less than the number of points required by other, conventional, methods. For example, for the same attractor N= 4 x I 0~points were used with the box-counting method in ref. [18], while N= 108 points were used with the correlation method in ref. [19]. Then, using eqs. (6) and (8), the relation between q and t was determined. Figure 3 shows the plot of q against ~(q) for the Hénon attractor as determined from this minimal spanning —

Volume 169, number 3

PHYSICS LETTERS A

tree edge length scaling behaviour. Errors have been estimated by varying the scaling range where the fits of eq. (8) are evaluated [20]. Near q= 0 one sees that a straight line, with slope approximately 1 / D (0), fits the data points rather well, illustrating the fact that in the case of the Hénon attractor the generalized dimensions for I q I <1 are nearly equal to the Hausdorff dimension. The deviation from a straight line for both more negative and positive q’s is an indicator for multifractal behaviour [211. When applying a Legendre transform to the r( q) curve in fig. 3 we get an f( a) spectrum which is in good accordance with previously published results [18,19,22,23]. For the purpose of comparison we mention here that we found D(0)=l.268±0.021, and D( 1) = 1.244 ±0.014. We have analyzed some small regions of the attractor, in particular a region around the second turnback [23]. The MST algorithm still converges accurately towards thef( a) curve. However, the result is not equivalent to the f( a) spectrum of the whole attractor due to the nonuniformity of the fractal measure on the attractor. This fact implies that the local density scales differently for different regions [19]. —

6. Summary We have introduced a new method for estimating the generalized dimensions of a point set consisting of a limited number of points, based on the use of the minimal spanning tree, a uniquely defined graph within the point set. Its use has several advantages. Firstly, it comes as close as possible for a discrete point set to the definition ofgeneralized dimensions based on partition sums containing an “infimum” term. Secondly, it is an estimator of the Hausdorff generalized dimensions instead of the Rényi generalized dimensions as most estimators. Thirdly, it is a very good estimator in the case one has only a limited number of points as is mostly the case in sciences like e.g. experimental physics. Finally, there are some analytical indications found by Steele [6] that there is indeed a deep link between minimal spanning trees and fractal dimensions. The application to the well-known Hénon attractor shows that the dimensions and f( a) spectrum determined with

21 September 1992

the minimal spanning tree method compare quite well with those determined using more conventional estimators. Applications to continuous trajectories are also possible. For example, the multifractal character of the Lorenz attractor has been studied by means ofthe MST method in ref. [24]. All in all we may conclude from this that the minimal spanning tree is a valuable addition to the apparatus for determining the dimensionality of a point set.

Acknowledgement We are grateful to P. Grassberger, J. Barrow, M. Jensen, D. Umberger and R. DomInguez-Tenreiro for valuable comments and advice. We acknowledge the IFIC for providing us with computing facilities. This work was partially supported by the Dirección General de Investigacion Cientifica y Técnica (project number PB86-0292-C04-03-04), Spain.

References [l]F.Hausdorff,Math.Ann. 79(1919) 157. [2] K.J. Falconer, The geometry of fractal sets (Cambridge Univ. Press, Cambridge, 1985). [3] J.B. Kruskal, Proc. Am. Math. Soc. 7 (1956) 48. [4] R.C. Prim, Bell Syst. Tech. J. 36 (1957) 1389. [5] V.K.M. Whitney, Commun. ACM 15 (1972) 273. [6] J.M. Steele, Ann. Prob. 16 (1988) 1767. [7] J.D. Barrow, 5.P. Bhavsar and D.H. Sonoda, Mon. Not. R. Astron. Soc. 216 (1985) 17; S.P. Bhavsar and N. Ling, Astrophys. J. 331 (1988) L63; N. Ling, New statistical approaches to galaxy clustering, Ph.D. thesis, University of Sussex (1987). [8] V.J. Martinez and B.J.T. Jones, Mon. Not. R. Astron. Soc. 242 517; V.J. (1990) Martinez, B.J.T. Jones, R. Domfnguez-Tenreiro and R. van de Weygaert, Astrophys. J. 357 (1990) 50. [9) T.C. Halsey, M.H. Jensen, L.P. Kadanoff, I. Procaccia and B.!. Shraiman, Phys. Rev. A. 33 (1986) 1141. [10) A.N. Kolmogorov and V.M. Tihomirov, Usp. Mat. Nauk 14(1959)3. [11] A. Rényi, Probability theory (North-Holland, Amsterdam, 1970). [12] H.G.E. Hentschel and I. Procaccia, Physica D 8 (1983) 435. [13] V.J. Martinez, La textura de l’univers, Ph.D. thesis, University of Valencia (1989). [14] R. Badii and A. Politi, Phys. Lett. A 104 (1984) 303; J. Theiler, Phys. Rev. A 34 (1986) 2427. [15] 0. Broggi, J. Opt. Soc. Am. B 5 (1988)1020. [16) P. Grassberger, private communication.

149

Volume 169, number 3

PHYSICS LETTERS A

[17] M. Hénon, Commun. Math. Phys. 50 (1976) 69. [18] P. Grassberger, R. Badii and A. Politi, J. Stat. Phys. 51 (1988) 179. [19] A. Arneodo, 0. Grasseau and E.J. Kostelich, Phys. Lett. A 124 (1987) 426. [20] M.H. Jensen, L.P. Kadanoff, A. Libchaber, I. Procaccia and J. Stavans, Phys. Rev. Lett. 55 (1985) 2798. [21] R. Benzi, 0. Paladin, 0. Parisi and A. Vulpiani, J. Phys. A 17(1984)3521.

150

21 September 1992

[22) G.H. Gunaratne andI. Procaccia, Phys. Rev. Lett. 59 (1987) 1377; D.A. Rusell, J.D. Hansen and E. Ott, Phys. Rev. Lett. 45 (1980) 1175; P. Grassberger and I. Procaccia, Phys. Rev. Lett. 50 (1986) 346. [23] M.H. Jensen, Phys. Rev. Lett. 60 (1988) 1680. [24] R. Dominguez-Tenreiro, L.J. Roy and V.5. Martinez, Prog. Theor.Phys. (1992), in press.