# Intrinsic Dimension Estimation: Relevant Techniques and a Benchmark Framework

Authors: P. Campadelli; E. Casiraghi; C. Ceruti; A. Rozza
Journal: Mathematical Problems in Engineering
 (2015)

---

## Abstract

When dealing with datasets comprising high-dimensional points, it is usually advantageous to discover some data structure. A fundamental information needed to this aim is the minimum number of parameters required to describe the data while minimizing the information loss. This number, usually called intrinsic dimension, can be interpreted as the dimension of the manifold from which the input data are supposed to be drawn. Due to its usefulness in many theoretical and practical problems, in the last decades the concept of intrinsic dimension has gained considerable attention in the scientific community, motivating the large number of intrinsic dimensionality estimators proposed in the literature. However, the problem is still open since most techniques cannot efficiently deal with datasets drawn from manifolds of high intrinsic dimension and nonlinearly embedded in higher dimensional spaces. This paper surveys some of the most interesting, widespread used, and advanced state-of-the-art methodologies. Unfortunately, since no benchmark database exists in this research field, an objective comparison among different techniques is not possible. Consequently, we suggest a benchmark framework and apply it to comparatively evaluate relevant state-of-the-art estimators.


---

## Body

## 1. Introduction

Since the 1950s, the rapid pace of technological advances allows us to measure and record increasing amounts of data, motivating the urgent need to develop dimensionality reduction systems to be applied on real datasets comprising high-dimensional points.To this aim, a fundamental information is provided by theintrinsic dimension (id) defined by Bennett [1] as the minimum number of parameters needed to generate a data description by maintaining the “intrinsic” structure characterizing the dataset, so that the information loss is minimized.More recently, a quite intuitive definition employed by several authors in the past has been reported by Bishop in [2], p. 314, where the author writes that “a set in D dimensions is said to have an id equal to d if the data lies entirely within a d-dimensional subspace of RD.”Though more specific and differentid definitions have been proposed in different research fields [3–5], throughout the pattern recognition literature the presently prevailing id definition views a point set as a sample set uniformly drawn from an unknown smooth (or locally smooth) manifold structure, eventually embedded in a higher dimensional space through a nonlinear smooth mapping; in this case, the id to be estimated is the manifold’stopological dimension.Due to the importance ofid in several theoretical and practical application fields, in the last two decades a great deal of research effort has been devoted to the development of effective id estimators. Though several techniques have been proposed in the literature, the problem is still open for the following main reasons.At first, it must be highlighted that though Lebesgue’s definition of topological dimension [6] (see Section 3.2) is quite clear, in practice its estimation is difficult if only a finite set of points is available. Therefore, id estimation techniques proposed in the literature are either founded on different notions of dimension (e.g., fractal dimensions, Section 3.2.1) approximating the topological one or on various techniques aimed at preserving the characteristics of data-neighborhood distributions, which reflect the topology of the underlying manifold. Besides, the estimated id value markedly changes as the scale used to analyze the input dataset changes [7] (see an example in Figure 1), and with the number of available points being practically limited, several methods underestimate id when its value is sufficiently high (namely, id⩾10). Other serious problems arise when the dataset is embedded in higher dimensional spaces through a nonlinear map. Finally, the too high computational complexity of most estimators makes them unpractical when the need is to process datasets comprising huge amounts of high-dimensional data.Figure 1
At very small scales (a) the dataset seems zero-dimensional; in this example, when the resolution is increased until including all the dataset (b) theid looks larger and seems to equal the embedding space dimension; the same effect happens when it is estimated at noise level (d); the correct id estimate is obtained at an intermediate resolution.In this work, after recalling the application domains of interest, we survey some of the most interesting, widespread used, and advancedid estimators. Unfortunately, since each method has been evaluated on different datasets, it is difficult to compare them by solely analyzing the results reported by the authors. This highlights the need of a benchmark framework, such as the one proposed in this work, to objectively assess and compare different techniques in terms of robustness with respect to parameter settings, high-dimensional datasets, datasets characterized by a high id, and noisy datasets.The paper is organized as follows: in Section2 the usefulness of the id knowledge is motivated and interesting id application domains profitably exploiting it are recalled; in Section 3 we survey notable state-of-the-art id estimators, grouping them according to the employed methods; in Section 4 we summarize mostly used experimental settings, we propose a benchmark framework, and we employ it to objectively assess and compare relevant id estimators; in Section 5 conclusions and open research problems are reported.


## 2. Application Domains

In this section we motivate the increasing research interest aimed at the development of automaticid estimators, and we recall different application contexts where the knowledge of the id of the available input datasets is a profitable information.In the field of pattern recognition, theid is one of the first and fundamental pieces of information required by several dimensionality reduction techniques [8–12], which try to represent the data in a more compact, but still informative, way to reduce the “curse of dimensionality” effects [13]. Furthermore, when using an autoassociative neural network to perform a nonlinear feature extraction, the id value d can suggest a reasonable value for the number of hidden neurons [14]. Indeed, a network with a single hidden layer of neurons with linear activation functions has an error function with a unique global minimum and, at this minimum, the network performs a projection on the subspace spanned by the first d principal components [15] estimated on the dataset (see 8.6.2 of [2]), with d being the number of hidden neurons. Besides, according to statistical learning theory [16], the capacity and generalization capability of a given classifier may depend on the id. More specifically, in the particular case of linear classifiers where the data are drawn from a manifold embedded through an identical map, the Vapnik-Chervonenkis (VC) dimension of the separation hyperplane is d+1 (see [16], pp. 156–158). Since the generalization error depends on the VC dimension, it follows that the generalization capability may depend on the id value d. Moreover, in [17] the authors mark that, in order to balance a classifier generalization ability and its empirical error, the complexity of the classification model should also be related to the id of the available dataset. Furthermore, since complex objects can be considered as structures composed by multiple manifolds that must be clustered to be processed separately, the knowledge of the local ids characterizing the considered object is fundamental to obtain a proper clustering [18].These observations motivate applications employing global or localid estimates to discover some structure within the data. In the following we summarize or simply recall some interesting examples [19–25].In [19] the authors introduce a fractal dimension estimator, called correlation dimension (CD) estimator (see Section 3.2.1), and show that the id estimate it computes is a reliable approximation of the strange attractor dimension in chaotic dynamical systems.In the field of gene expression analysis, the work proposed in [20] shows that the id estimate computed by the nearest neighbor estimator (described in [26] and Section 3.2.2) is a lower bound for the number of genes to be used in supervised and unsupervised class separation of cancer and other diseases. This information is important since generally used datasets contain large number of genes and the classification results strongly depend on the number of genes employed to learn the separation criteria.In [21], the authors show that id estimation methods being derived from the basis theory of fractal dimensions ([7, 19, 27], see Section 3.2.1) can be successfully used to evaluate the model order in signals and time series, which is the number of past samples required to model the time series adequately and is crucial to make reliable predictions. This comparative work employs fractal dimension estimators, since the domain of attraction of nonlinear dynamic systems has a very complex geometric structure, which could be captured by closely related studies on fractal geometry and fractal dimensions.A noteworthy research work in the field of crystallography [22] employs the fractal CD estimator [19] followed by a correction method [28] that, according to the authors, “is needed because the CD estimator, to give correct estimations of the id, requires an unrealistically large number of points.” Anyway, the experimental results show that id is a useful information to be exploited when analyzing crystal structures. This study not only proves that id estimates are especially useful when dealing with practical tasks concerning real data, but also underlines the need to compute reliable estimates on datasets drawn from manifolds characterized by high id and embedded in spaces of much greater dimensionality.The work of Carter et al. [23] is very interesting and notable because it is one of the first considering that the input data might be drawn from a multimanifold structure, where each submanifold has a (possibly) different id. To separate the manifolds, the authors compute local id estimates, by applying both a fractal dimension estimator (namely, MLE [27]; see Section 3.2.1) and a nearest neighbor-based estimator (described in [29, 30]; see Section 3.2.2) on properly defined data neighborhoods. The authors then show that the computed local ids might be helpful for the following interesting applications: (1) “Debiasing global id estimates”: the negative bias caused both by the limited number of available sample points and by thecurse of dimensionality is reduced by computing global id estimates through a weighted average of the local ones, which assign greater importance to the points away from the boundaries. However the authors themselves note that this method is only applicable for data with a relatively low id, since in high dimensions the points lie nearby the boundaries [31]. (2) “Statistical manifold learning”: the local id estimates are used to reduce the dimension of statistical manifolds [32], that is, manifolds whose points represent a pdf. When this step is applied as the first step of document classification applications, and analysis of patients’ samples acquired in the field of flow cytometry, it allows us to obtain lower dimensional points showing a good class separation. (3) “Network anomaly detection”: considering that the overall complexity of a router network is decreased when few sources account for a disproportionate amount of traffic, a decrease in the id of the entire network is searched for. (4) “Clustering”: problems of data clustering and image segmentation are dealt with by assuming that different clusters and image patches belong to manifold structures characterized by different complexity (and ids).In [24], to the aim of analyzing gene expression time series, the authors compute id estimates by comparing the fractal CD estimator and the nearest neighbor (NN) estimator [26]. The results on both simulated and real data show that NN seems to be more robust than CD with respect to nonlinear embedding and the underlying time-series model.In the field of geophysical signal processing, hyperspectral images, whose pixels represent spectra generated by the combination of an unknown set of independent contributions, called endmembers, often require estimating the number of endmembers. To this aim, the proposal in [25] is to substitute state-of-the-art algorithms specifically designed to solve this task with id estimators. After motivating the idea by describing the relation between the id of a dataset and the number of endmembers, the authors choose to experiment two fractal id estimators [7, 19] and a nearest neighbor-based one [33]. They obtain the most reliable results with the latest one after opportunely tuning the number of nearest neighbors to be considered.Finally, other noteworthy examples of research works that profitably exploitid and estimate it by usually applying fractal dimension estimators concern financial time series prediction [34], biomedical signal analysis [35–37], analysis of ecological time series [38], radar clutter identification [39], speech analysis [40], data mining and low dimensional representation of (biomedical) time series [41], and plant traits representation [42].


## 3. Intrinsic Dimension Estimators

In this section we survey some of the most notable, recent, and effective state-of-the-artid estimators, grouping them according to the main ideas they are based on.Specifically, in Section3.1 we describeprojectiveidestimators, which basically process a dataset PN≡{pi}i=1N⊆RD to identify a somehow appealing lower dimensional subspace where to project the data; the space dimension of the identified subspace is viewed as the id estimate.More recent projectiveid estimators exploit the assumption of datasets PN≡{pi}i=1N⊆RD being uniformly drawn from a smooth (or locally smooth) manifold M⊆Rd, embedded into a higher D-dimensional space through a nonlinear map; this is also the basic assumption of all the other groups of methods that will be referred to astopologicalidestimators (see Section 3.2) andgraph-basedidestimators (see Section 3.3).We note that the taxonomy we are using to group the reviewed methods is different from the one, commonly used by several authors in the past (as an example, see [43]), that viewed methods asglobal, when id estimation is performed by considering a dataset as a whole, orlocal, when all the data neighborhoods are analyzed separately and an estimate is computed by combining all the local results. All the recent methods have abandoned the global approach since it is now clear that analyzing a dataset at its biggest scale cannot produce reliable results. They thus estimate the global id by somehow combining local ids. This way of proceeding comes from the assumption that the underlying manifold is locally smooth.


### 3.1. Projectiveid Estimators

The first projectiveid estimators introduced in the literature explicitly compute the mapping that projects input points PN∈RD to the subspace M⊆Rd minimizing the information loss [27, 43] and therefore view the id as the minimal number of vectors linearly spanning the subspace M. It must be noted that, since these methods were originally designed for exploratory data analysis and dimensionality reduction, they often require the dimensionality of M (the id to be estimated) to be provided as input parameter. However, their extensions and variants include methodologies to automatically estimate id.Most of the projectiveid estimators can be grouped into two main categories: projection techniques based on Multidimensional Scaling (MDS) [44, 45] or its variants, which tend to preserve as much as possible pairwise distances among the data, and projection techniques based on Principal Component Analysis (PCA) [15, 46] and its variants that search for the best projection subspace M minimizing the projection error.Some of the best known examples ofMDS algorithms are MDSCAL [47–52], Bennett’s algorithm [1, 53], Sammon’s mapping [54], Curvilinear Component Analysis (CCA) [55], ISOMAP [56], and Local Linear Embedding (LLE) [57]. As shown by experiments reported in [56, 57] ISOMAP and variants of LLE compute the most reliable id estimates. We believe that their better performance is due to the fact that both ISOMAP and LLE have been the first projective methods based on the assumption that the input points are drawn from an underlying manifold, whose curvature might affect the precision of data neighborhoods computed by employing the Euclidean distance. However, as noted in [58, 59], these algorithms have shown that they suffer of all the major drawbacks affecting MDS-based algorithms, which are too much tied by the preservation of the pairwise distance values. Besides, as highlighted in [30], ISOMAP as an id estimator, as well as other spectral based methods like PCA, relies on a specific estimated eigenstructure that may not exist in real data. Regarding LLE, it either requires the id value to be known in advance or may automatically estimate it by analyzing the eigenvalues of the data neighborhoods [60]; however, as outlined in [7, 27], id estimates computed by means of eigenvalue analysis are as unreliable as those computed by most PCA-based approaches. Moreover, in [46] it is noted that methods such as LLE are based on the solution of a sparse eigenvalue problem under the unit covariance constraint; however, due to this imposed constraint, the global shape of the embedded data cannot reflect the underlying manifold.P
C
A [15, 46] is one of the most cited and well-known projective id estimators, often used as the first step of several pattern recognition problems, to compute low dimensional representations of the available datasets. When PCA is used for id estimation, the estimate is the number of “most relevant” eigenvectors of the sample covariance matrix, also called principal components (PCs). Due to the promising dimensionality reduction results, several PCA-based approaches, both deterministic and probabilistic, have been published. Among deterministic approaches, we recall the Kernel PCA (KPCA) [61] and the local PCA (LPCA) [62] and its extensions to automatically select the number of PCs [63, 64]. We observe that the work presented in [64] is one of the first works that estimates id by considering an underlying topological structure and therefore applies LPCA on data neighborhoods represented by an Optimally Topology Preserving Map (OTPM) built on clustered data (given an input dataset PN, its OTPM is usually computed through Topology Representing Networks (TRNs); these are unsupervised neural networks [65] developed to map PN to a set of neurons whose learnt connections define proximities in PN. These proximities correspond to the optimal topology preserving Voronoi tessellation and the corresponding Delaunay triangulation. In other words, TRNs compute connectivity structures that define and perfectly preserve the topology originally present in the data, forming a discrete path-preserving representation of the inner (topological) structure of the topological manifold underlying the dataset PN). The authors of this method state that their approach is more efficient and less sensitive to noise with respect to the PCA-based approaches. However, they do not show any experimental comparison and, besides, their algorithm employs critical thresholds and a data clustering technique whose result heavily influences the precision of the computed estimate [27].The usage of a probabilistic approach has been firstly introduced by Tipping and Bishop in [66]. Considering that deterministic methodologies lack an associated probabilistic model for the observed data, their Probabilistic PCA (PPCA) reformulates PCA as the maximum likelihood solution of a specific latent variable model. PPCA and its extensions to both mixture and hierarchical mixture models have been successfully applied to several real problems, but they still provide an id-estimation mechanism depending on critical thresholds. This motivates its subsequent variants [67] and developments, whose examples are Bayesian PCA (BPCA) [68] and two Bayesian model order selection methods introduced in [69, 70]. In [71] the asymptotic consistency of id estimation by (constrained) isotropic version of PPCA is shown with numerical experiments on simulated and real datasets.While the aforementioned methods have been simply recalled since theirid estimation results have shown to be unreliable [7, 27], in the following recent and promising proposals are described with more details.The Simple Exponential FamilyPCA (SePCA) [72] has been developed to overcome the assumption of Gaussian-distributed data that makes it difficult to handle all types of practical observations, for example, integers and binary values. SePCA achieves promising results by using exponential family distributions; however, it is highly influenced by critical parameter settings and it is successful only if the data distribution is known, which is often not the case, specially when highly nonlinear manifold structures must be treated.In [73] the authors propose the Sparse Probability PCA (SPPCA) as a probabilistic version of the Sparse PCA (SPCA) [74]. Precisely, SPCA selects id by forcing the sparsity of the projection matrix that is the matrix containing the PCs. However, based on the consideration that the level of sparsity is not automatically determined by SPCA, SPPCA employs a Bayesian formulation of SPCA, achieving sparsity by employing a different prior and automatically learning the hyperparameter related to the constraint weight through Evidence Approximation ([75] Section 3.5). The authors’ results and also the results of the comparative evaluation proposed in [76] show that this method seems to be less affected by the problems of the aforementioned projective schemes.An alternative method (MLSVD) [77] applies Singular Value Decomposition (SVD), basically a variant of PCA, locally and in a multiscale fashion to estimate the id characterizing D-dimensional datasets drawn from nonlinearly embedded d-dimensional manifolds M corrupted by Gaussian noise. Precisely, exploiting the same ideas of the theoretical PCA-based id estimator presented in [78], the authors note that the best way to avoid the effects of the curvature (induced by the nonlinearity of the embedding) is to apply SVD locally, that is, in hyperspheres B(p,r) centered on the data points p and having radius r. However, the choice of  r is constrained by the following considerations: (1) r must be big enough to have at least k≥d neighbors, (2) r must be small enough to ensure that M∩B is linear (or at least smooth), and (3) r must be big enough to ensure that the effects of noise are negligible. When these three constraints are met, the tangent space TMd(p,r), computed by applying SVD on the k neighbors, is a good approximation of the tangent space of M∩B and the number of its relevant eigenvalues corresponds to the (local) id of M. To find a proper value for r, the authors propose a multiscale approach that applies SVD on neighborhoods B(p,rs) whose radius varies in a range rs∈{rL,…,rH}. This allows us to compute D scale-dependent, local singular values λ1p,rs≥⋯≥λD(p,rs); using a least squares fitting procedure the SVs can be expressed as functions of r whose analysis allows us to identify the range of scales [rmin,…,rmax] not influenced by either noise or curvature. Finally, in the range rs=[rmin,…,rmax] the squared SVs are analyzed to get the id estimate d^ that maximizes the gap Δ(j)=λj(p,rs)-λj+1(p,rs) for the largest range of rs. The proposed algorithm has been evaluated on unit d-dimensional hyperspheres and cubes embedded in R100 and affected by Gaussian noise. The reported results are very good, while other ten well-known methods [19, 23, 27, 30, 79–81] show that the ids estimated on the same datasets are unreliable also in the absence of noise.


### 3.2. Topological Approaches

Topological approaches forid estimation consider a manifold M⊆Rd embedded in a higher dimensional space RD through a proper (locally) smooth map ϕ:M→RD and assume that the given dataset is PN={pi}i=1N={ϕ(xi)}i=1N⊂RD, where xi are independent identically distributed (i.i.d.) points drawn from M through a smooth probability density function (pdf) f:M→R+.Under this assumption theid to be estimated is the manifold’s topological dimension, defined either through the firstly proposed Brouwer Large Inductive Dimension [82] or the equivalent Lebesgue Covering Dimension [83]. Since Brouwer’s definition has been soon neglected by mathematicians for its difficult proof [83], the commonly adopted topological dimension definition is Lebesgue’s Covering Dimension, reported in the following.Definition 1 (cover).
Given a topological spaceX, a cover of a set Y⊆X is a countable collection C={Ci} of open sets such that each Ci⊂X and ⋃iCi⊇Y.Definition 2 (refinement of a cover).
A refinement of a coverC of a set Y is another cover C′ such that each set in C′ is contained in some sets of C.Definition 3 (topological dimension (Lebesgue Covering Dimension)).
Given the aforementioned definitions, the topological dimension of the topological spaceX, also called Lebesgue Covering Dimension, is d if every finite cover of X admits a refinement C′ such that no subset of X has more than d+1 intersecting open sets in C′. If no such minimal integer value exists, X is said to be of infinite topological dimension.To our knowledge, at the state of the art only two estimators have been explicitly designed to estimate the topological dimension.One of them, the Tensor Voting Framework (TVF) [84] and its variants [85], relies on the usage of an iterative information diffusion process based on Gestalt principles of perceptual organization [86]. TVF iteratively diffuses local information describing, for each pi∈PN, the tangent space approximating the underlying neighborhood of M. To this aim, the information diffused at each iteration is second-order symmetric positive definite tensors whose eigenvectors span the local tangent space. Practically, during the initialization step a ball tensor Ti0, which is an identity matrix representing the absence of orientation, is used to initialize a token pi for each point pi as {pi=(pi,Ti0)}i=1N. During iteration t each token pi “generates” the set of tensors Ti,jtj≠i that enact as votes cast to neighboring tokens; precisely, Ti,jt is sent to the jth neighbor, and it encodes information related to the local tangent space estimate in pi at time t and decays as the curvature and the distance from the jth neighbor increase. On the other side, at iteration t each token pj receives votes that are summed to update the pj’s tensor as Tjt+1=∑i≠j‍Ti,jt; this essentially refines the estimate of the local tangent space in pj. Based on the definition of topological dimension provided by Brouwer [82], in [87] it is noted that TVF can be employed to estimate the local ids by identifying the number of most relevant eigenvalues of the computed second-order tensors. Although interesting, this method has a too high computational cost, which makes it unfeasible for spaces of dimension D≥4.From the definition of Lebesgue Covering Dimension it can be derived [88] that the topological dimension of any M⊆Rd coincides with the affine dimension d of a finite simplicial complex (a simplicial complex in Rd has affine dimension d if it is a collection of affine simplexes in Rd, having at most dimension d or having at most d+1 vertices) covering M. This essentially means that a d-dimensional manifold could be approximated by a collection of d-dimensional simplexes (each having at most d+1 vertices); therefore, the topological dimension of M could be practically estimated by analyzing the number of vertices of the collection of simplexes estimated on PN. To this aim, in [89] a method is proposed to find the number of relevant positive coefficients that are needed to reconstruct each pi∈PN from a linear combination of its k neighbors, where k is a parameter to be manually set in the range d<k≤D+1. This algorithm is based on the fact that neighbors with positive reconstruction coefficients are the vertices of a simplex with dimension equal to the topological dimension of M. Practically, to ensure that k>d, its value is set to D, the reconstruction coefficients are calculated by means of an optimization problem constrained to be nonnegative, and the coefficients bigger than a user-defined threshold are considered as the relevant ones. The id estimate is then computed by employing two alternative approaches: the first one simply computes the mode of the number of relevant coefficients for each neighborhood and the second one sorts in descending order the coefficients computed for each neighborhood, computes the mean c¯ of the sorted coefficients, and estimates id as the number of relevant values in c¯. Note that, since k>d, this method is strongly affected by the curvature of the manifold when the id is big enough. Indeed, the results of the reported experimental evaluation make the authors assert that the method works well only on noisy-free data of low id (id≤6), under the assumption that the sampling process is uniform and the data points are sufficient.Though interesting, both approaches have shown to be effective only for manifolds of low curvature as well as lowid values.In the following we survey otherid estimators employing two different estimation approaches, which allow us to categorize them. More precisely, in Section 3.2.1fractalidestimators are described, which estimate different fractal dimensions since they are good approximations of the topological one; in Section 3.2.2nearest neighbors-based (NN) idestimators are recalled, which are often based on the statistical analysis of the distribution of points within small neighborhoods.


#### 3.2.1. Fractalid Estimators

Since topological dimension cannot be practically estimated, several authors implicitly assume thatM has a somehow fractal structure (see [90] for an exhaustive description of fractal sets) and estimate id by employing fractal dimension estimators, the most relevant of which are surveyed in this section.Very roughly, since the basis concept of all fractal dimensions is that the volume of ad-dimensional ball of radius r scales with its size s as rd [90, 91], all fractal dimension estimators are based on the idea of counting the number of observations in a neighborhood of radius r to (somehow) estimate the rate of growth of this number. If the estimated growth is rd, then the estimated fractal dimension of the data is considered to be equal to d.Note that all the derived estimators have the fundamental limitation that, in order to get an accurate estimation, the sizeN of the dataset with id  d has to satisfy the inequality proved by Eckmann and Ruelle in [92] for the correlation dimension estimator (CD [19], see below):(1)d<2log⁡1/ρ∗log⁡N,beingρ=rD≪1,12N2rDd≫1.This will lead to a large value of N, even for a dataset with lower id.Among fractal dimension estimators, one of the most cited algorithms is presented in [19] and will be referred to as CD in the following. It is an estimator of the correlation dimension (dimCorr), whose formal definition is as follows.Definition 4 (correlation dimension).
Given a finite sample setPN, let(2)CNr=2NN-1∑i=1,i<jNIr-pi-pj,where · is the Euclidean norm and I(·) is the step function used to simulate a closed ball of radius r centered on each pi (I(y)=0 if y<0, and I(y)=1 otherwise). Then, for a countable set, the correlation dimension dimCorr is defined as(3)dimCorr=limr→0limN→∞log⁡CNrlog⁡r.In practice CD computes an estimate, d^, of dimCorr by computing CN(r) for different ri and applying least squares to fit a line through the points log⁡ri;log⁡CN(ri). It has to be noted that, to produce correct id estimates, this estimator needs a very large number of data points [22], which is never available for practical applications; however, the computed unreliable estimations can be corrected by the correction method proposed in [28].The relevance of theCD estimator is shown by its several variants and extensions. An example is the work proposed in [91], where the authors propose a normalized CD estimator for binary data and achieve estimates approximating those computed by CD.SinceCD is heavily influenced by the setting of the scale parameters, in [93] the authors estimate the id by computing the expectation value of dimCorr through maximum likelihood estimate of the distribution of distances among points. The estimated d^ is computed as (4)d^=-1Q∑k=1Qrk-1,where Q is the set Q=rk∣rk<r and rk is the Euclidean distance between two generic data points and r is a real value, called cut-off radius.To develop an estimator more efficient thanCD, in [94] the authors choose a different notion of Fd, namely, the Information Dimension dimI:(5)dimI=-limδ→0∑i=1Nδprilog⁡prilog⁡δ,where N(δ) is the minimum number of δ-sized hypercubes covering a topological space and pri is the probability of finding a point in the ith hypercube. Noting that when the scale δ in (5) is big enough the different coverings used to estimate dimI could produce different values for N(δ), the authors look for the covering composed by the minimum number Nmin(δ) of nonempty sets. Similar to the CD algorithm, the id is the average slope of the curve obtained by fitting the points with coordinates log⁡δ;∑i=1Nmin(δ)prilog⁡pri.This algorithm is compared with theCD estimator, and the experimental tests show that both methods compute the same estimates. However, the achieved computation time is much lower than that of CD.Considering thatCD can severely underestimate the topological dimension if the data distribution on the manifold is nearly nonuniform, in [7] the author proposes the Packing Number (PN), a fractal dimension estimator that approximates the Capacity Dimension (dimCap). This choice is motivated by the fact that dimCap does not depend on the data distribution on the manifold and, if both dimCap and the topological dimension exist (which is certainly the case in machine learning applications), the two dimensions agree. To formally define dimCap, the ϵ-covering number N(ϵ) of a set S⊂X must be defined; N(ϵ) is the minimum number of open balls B(x0,ϵ)=x∈X:x0-x<ϵ whose union is a covering of S, where · is a distance metric. The definition of dimCap of S⊂X is based on the observation that the covering number N(ϵ) of a d-dimensional set is proportional to ϵ-d:(6)dimCap=-limϵ→0log⁡Nϵlog⁡ϵ.Since the estimation ofN(ϵ) is computationally expensive, based on the relation N(ϵ)≤NPack(ϵ)≤N(ϵ/2), the authors employ the ϵ-Packing Number NPack(ϵ), defined in [95] as the maximum cardinality of an ϵ-separated set. Employing a greedy algorithm to compute NPack(ϵ), the estimate, d^, of dimCap is computed as(7)d^ϵ1,ϵ2=-log⁡NPackϵ1-log⁡NPackϵ2log⁡ϵ1-log⁡ϵ2.To estimate d^ a greedy algorithm is used; however, as noted by the author, the dependency of d^ with respect to the order in which the points are visited by the greedy algorithm introduces a high variance. To avoid this problem, the algorithm iterates the procedure M times on random permutations of the data and considers the average as the final id estimate. The comparative evaluation with the CD estimator makes the authors assert that PN “seems more reliable if data contains noise or the distribution on the manifold is not uniform.” Unfortunately, also this method is scale-dependent.To avoid any scale-dependency in [79] the authors propose an estimator (Hein) based on the asymptotes of a smoothed version of (2), obtained by replacing the step function I(·) with a suitable kernel function. Precisely, they define(8)UN,h,d=2NN-1∑1≤i<j≤NN1hdKhpi-pjh2,where Kh is a kernel function with bandwidth h and d is the assumed dimensionality of the manifold from which the points are sampled. Note that, to guarantee the converge of (8), the bandwidth h has to fulfill the constraint limN→∞(Nhd)=∞. For this reason the authors formalize h as a function of N and, to achieve scale-independency, propose a method that estimates the id by analyzing the convergence of U(N,h,d) when varying the parameters N and d. Precisely, the dataset is subsampled to create sets of different cardinalities nsub∈Nsub={N,N/2,N/3,N/4,N/5} and the D curves whose points have coordinates (U(nsub,h(nsub),d),nsub) are considered. Employing this information the following id estimator is proposed:(9)Sloped=maxnsub∈Nsub⁡∂Unsub,hnsub,d∂n,d^=argmind∈1,…,D⁡Sloped.This work is notable since the empirical tests are performed on synthetic datasets specifically designed to study the influence of high curvature as well as noise on the proposed estimator. The usefulness of these datasets is confirmed by the fact that they have been also employed to assess several subsequent methods [76, 96].In [97] the authors present a fractal dimension estimator derived by the analysis of a vector quantizer applied to datasets PN⊆RD. Considering the codebook Y={y1,…,yk}⊂RD containing k code-vectors yi, ak-point quantizer is defined by a measurable function Qk:RD→Y, which brings each data point to one of the code-vectors in Y. This partitions the dataset into k so-calledquantizer cellsSi={pi∈PN:Qk(pi)=yi}, where log2(k) is called therate of the quantizer. Being X a random vector distributed according to a probability distribution ν, thequantization error is er(Qk∣ν)=(Eν[X-Qk(X)r])1/r, where r∈[1,∞) and · is the Euclidean norm in RD. Given the set Qk of all D-dimensional k-point quantizers, the performance achieved by an optimal k-point quantizer on X is er∗(Qk∣ν)=infQk∈Qk(er(Qk∣ν)). When the quantizer rate is high, the quantizer cells can be well approximated by D-dimensional hyperspheres with radius equal to ϵ and centered on each code-vector yi∈Y. In this case, the regularity of ν ensures that the probability of such balls is proportional to ϵ1/d, and it can be shown [98] that er∗(Qk∣ν)≈k-1/d. This is referred to as thehigh-rate approximation and motivates the definition of quantization dimension of order r:(10)drν=-limk→∞log⁡klog⁡er∗k∣ν.The theory of high-rate quantization [98] confirms that, for a regular ν supported on the manifold M, dr(ν) exists for each 1≤r≤∞ and equals the intrinsic dimension of M. Furthermore, the limit k→∞ allows us to motivate the relation between the quantization dimension and the Capacity Dimension. Indeed, according to the theory of high-rate quantization [98, 99], there exists a decreasing sequence ϵk, such that for sufficiently large values of k (i.e., in the high-rate regime that is when k→∞) the ratio -(log⁡k)/(log⁡er∗(k∣ν)) can be approximated increasingly finely, both from below and from above, by quantities converging to the common value dimCap. To practically compute an estimate of the quantization dimension, having fixed the value of r, the authors select a range k1≤k≤k2 of codebook sizes and design a set of quantizers {Qk}k=k1k2 giving good approximations e^r(k∣ν) of er∗(k∣p) over the chosen range of k. An id estimate is obtained by fitting the points with coordinates (log⁡(k);-log⁡e^r(k∣ν)) and measuring the average slope over the chosen range k. Though the authors mention that their algorithm is less affected by underestimation biases than neighborhood-based methods (see Section 3.2.2), in [23] this statement is confuted with theoretical arguments.


#### 3.2.2. Nearest Neighbors-Basedid Estimators

In this section we consider estimators, referred to asNN estimators in the following, that describe data-neighborhoods’ distributions as functions of d. They usually assume that close points are uniformly drawn from small d-dimensional balls (hyperspheres) Bd(x,r) having radius (a small radius r→0∈R+ guarantees that samples included into Bd(x,r), being less influenced by the curvature induced by the map ϕ, are approximating well enough the intrinsic structure of the underlying portion of M) r→0∈R+ and being centered on x∈M.Practically, given an input datasetPN, the value of functions f(d) is computed by approximating the sampling process related to Bd through the k-nearest neighbor algorithm (kNN).AmongNNid estimators, Trunk’s method [100] is often cited as one of the first methods. It formulates the distribution function, f(d), with an ad hoc statistic based on geometric considerations concerning angles; in practice, having fixed a threshold γ and a starting value for the parameter k, it applies kNN to find the neighbors of each pi∈PN and calculates the angle νi between the (k+1)th-nearest neighbor and the subspace spanned by the k-nearest neighbors. Considering a threshold parameter γ, if (1/N)∑i=1N‍νi≤γ, then k is considered as the id estimate; otherwise, k is incremented by 1 and the process is repeated. The main limitation of this method is the difficult choice of a proper value for γ.The work presented by Pettis et al. [26] is notable since it is one of the first works providing a mathematical motivation for the use of nearest-neighbor distances.Indeed, for an i.i.d. samplePN⊆RD drawn from a density distribution p(x) in Rd, the following approximation holds:(11)kN≃pxVdrd,where k is the number of nearest neighbors to x within the hypersphere Bd(x,r) of radius r and centered on x and V(d) is the volume of the (unit d-dimensional) ball in Rd.This means that the proportion of sample points falling inBd(x,r) is roughly approximated by p(x) times the volume of Bd(x,r). Since this volume grows as rd, assuming the density p(x) to be a constant, it follows that the number of samples in Bd(x,r) is proportional to rd. From the relationship in (11), and assuming that the samples are locally uniformly distributed, the authors derive an id estimator for d:(12)d^=r¯kkr¯k+1-r¯k,where r¯k is the average of the distances from each sample point to its kth nearest neighbors; defining ri(k) as the distance between xi and its kth-nearest neighbor, r¯k is expressed as r¯k=(1/N)∑i=1N‍ri(k).Since this algorithm is limited by the choice of a suitable value for parameterk, in [63] the authors propose a variant which considers a range of neighborhood sizes [kmin,kmax]. However, in the same work the authors themselves show that this technique generally yields an underestimate of the id when its value is high.Taking into account relation (11), in [101] the number NBd of data points in Bd(x,r) is described by a polynomial f(r)=∑s=0d‍βsrs of degree d. In practice, considering pi,pk∈PN, calling rik=pi-pk the interpoint distances, and being r=mini,k=1Nrik, R a parameter adaptively estimated (to estimate R by means of PN, the radius value corresponding to the first significant peak of the histogram of the rijs is found), a set of n radius values r=rj=r+j(R-r)/nj=1n is selected and used to calculate n pairs rj,f^(rj)j=1n, where f^rj=#rik<rji,k=1N is the number of interpoint distances strictly lower than rj. To estimate the coefficients βjj=1D, the computed pairs are fit by a least squares fitting procedure that estimates exactly D+1 coefficients. Since by hypothesis the degree of f is d, the significance test described in [17] is used to estimate the degree d^ of f^, which is considered as the id estimate. The comparative evaluation of this algorithm with the well-known Maximum Likelihood Estimator (MLE) [27] and its improved version [102], both described below, has shown that it is more robust than them when dealing with high-dimensional datasets.M
L
E [27], one of the most cited estimators, treats the neighbors of each point pi∈PN as events in a Poisson process and the distance r(j)(pi) between the query point pi and its jth nearest neighbor as the event’s arrival time. Since this process depends on d, MLE estimates id by maximizing the log-likelihood of the observed process. In practice a local id estimate is computed as(13)d^pi,k=1k∑j=1klog⁡rk+1pirjpi-1.Averaging thed^(pi,k)s, the global id estimate is d^(k)=(1/N)∑i=1N‍d^(pi,k).The theoretical stability of the proposedid estimator for data living in C1 submanifold of RD, d≤D, and for data in an affine subspace of RD has been proved, respectively, in [103, 104]. Though the authors’ comparative evaluation shows the superior performance of the proposed estimator with respect to the CD estimator [19] (see Section 3.2.1) and the NN estimator [26], they further improve it by removing its dependency from the parameter k; to this end, different values for k are adopted and the computed results are averaged to obtain the final id estimate: d^=(1/t)∑k∈{k1,…,kt}‍d^(k).Considering that, in practice,MLE is highly biased for both large and small values of k, a variant of MLE is proposed in [102], where the arithmetic mean is substituted with the harmonic average, leading to the following estimator: d^(k)=(1/N)∑i=1N(1/d^pi,k)-1.Though the proposal in [102] seems to achieve more accurate results, it is based on the assumption that neighbors surrounding each pi are independent, which is clearly incorrect. To cope with this problem, in [105] an interesting regularized version of MLE applies a regularized maximum likelihood technique to distances between neighbors. The comparative evaluation with the aforementioned MLE methods [27, 102] makes the authors state that, though the method might be the first to converge to the actual estimate given enough data points, its estimation accuracy is comparable to that achieved by the competing schemes.In [59, 106] a further improvement of MLE is presented; it achieves a better performance by substituting euclidean distances with geodesic ones.Despite the good results achieved byMLE-based approaches, these techniques have shown to be affected by the curvature induced by ϕ on the manifold neighborhoods approximated by kNN. To reduce this effect, various id estimators have been proposed in [76, 107]; here, we review those achieving the most promising experimental results.In [107] the authors firstly propose a family of id estimators (MiNDML∗), which exploit the pdfg(r;k,d) describing the distance r(1)(x) between the center x of Bd(x,r), x∈M, r→0+ and its nearest neighbor. Briefly, formulating g(r;k,d) as a function of the id value d (g(r;k,d)=kdrd-1(1-rd)k-1), the id estimator is computed by a maximum likelihood approach.After noting that this algorithm is still affected by a bias causing underestimations when the dataset dimensionality becomes sufficiently high (i.e.,id≥10), the authors present theoretical considerations which relate the bias to the fact that id estimators based on nearest-neighbor distances are often founded on statistics derived under the assumption that the amount of available data is unlimited, which is never the case in practical applications. Based on these considerations, two different estimators, named MiNDKL and IDEA, are presented.M
i
N
D
K
L compares the empirical pdf of the neighborhood distances computed on the dataset (gData) with the distribution of the neighborhood distances computed from points uniformly drawn from hyperspheres of known increasing dimensionality (gSphered). The id estimate is the dimensionality that minimizes their Kullback-Leibler divergence KL(gData,gSphered), which is evaluated by means of the data-driven technique proposed in [108].I
D
E
A relies on the authors’ observation that the quantities 1-r(j)(pi)/r(k+1)(pi) are distributed according to the beta distribution β1,d with parameters 1 and d, respectively. Therefore, since Eβ1,d=m=1/1+d, a consistent id estimator d^≃d equals(14)d^=m^1-m^≃d=m1-m,wherem^=1Nk∑i=1N∑j=1krjpirk+1pi≃m.To reduce the effect of the aforementioned bias, IDEA finally applies an asymptotic correction step that, inspired by the correction method presented in [28], models the underestimation error by considering both the base algorithm and the given dataset.Motivated by the promising results achieved byMiNDKL, in [76] the authors propose its extension, namely, DANCo; it further reduces the underestimation effect by combining an estimator employing normalized nearest-neighbor distances with one employing mutual angles. More precisely, DANCo compares the statistics estimated on PN with those estimated on (uniformly drawn) synthetic datasets of known id. The comparisons are performed by two Kullback-Leibler divergences applied to the distribution of normalized nearest-neighbor distances g(r;k,d), having g(r;k,d)=kdrd-1(1-rd)k-1, and the distribution of pairwise angles q(x;ν,τ), q(x;ν,τ) being the von Mises-Fisher distribution (VMF) [109] with parameters ν and τ.Theid estimate d^ is the one minimizing the sum of the two divergences:(15)d^=argmind∈1,…,D⁡KLgData,gSphered+KLqData,qSphered.A fast implementation of DANCo (Fast-DANCo) is also developed. Comparative evaluations show that this algorithm achieves promising results (as shown in [76] and Section 4).Another work, which is notable because the authors not only prove the consistency in probability of the presented estimators but also derive upper bounds (see (19) below) on the probability of the estimation-error for finite, and large enough, values of N, is proposed in [33]. More precisely, the authors introduce two estimators by firstly defining a function η:RD×R→R+ slowly varying near the origin (see [33] for a detailed description and motivation of this assumption). The function η is then used to express the logarithm of the probability of a point p of being in the hypersphere BD(pi,r): logPp∈BD(pi,r)=logη(p,r)+dlog(r), having Pp∈BD(pi,r)=η(p,r)rd.Considering thatP(p∈B(pi,rk(pi)))≈k/n for N big enough, the authors derive the following system of equations:(16)log⁡kn≈log⁡ηpi,r+d^pilog⁡rkpi,log⁡k2n≈log⁡ηpi,r+d^pilog⁡rk/2pi,and solve it for d^(pi) to obtain a local id estimate:(17)d^pi=log⁡2log⁡rkpi/rk/2pi.The two proposed estimators are then computed either by averaging (d^avg) or by voting (d^vote):(18)d^avg=1N∑i=1Nd^pi,d^vote=argmaxd′∈N+⁡#d^pi=d′i=1N,where #condi=1N denotes the number of points pi for which cond is verified.Under differentiability assumptions on the functionη and regularity assumptions on M the authors prove the consistency in probability of their estimators and provide upper bounds (see (19)) on the probability of the estimation-error for finite, and large enough, values of N. However, the derived bounds depend on unknown universal constants c,c′,c′′>0:(19)Pd^avg≠d≤exp⁡-c′NDcdk2,Pd^vote≠d≤exp⁡-c′′Ncdk2.


### 3.3. Graph-Basedid Estimators

As noted in [96], the work of [110] has cleared up the fact that theories underlying graphs can be applied to solve a variety of statistical problems; indeed, also in the field of id estimation various types of graph structures have been proposed [29, 96, 111, 112] and used for id estimation. Examples are the kNN graph (kNNG) [30], the Minimum Spanning Tree (MST) [113] and its variation, the geodesic MST (GMST) [29]; and the sphere of influence graph (SIG) [114] and its generalization, the k-sphere of influence graph (kSIG) [96].Given a sample setPN={pi}i=1N a graph built on PN, usually denoted by G(PN)=({pi}i=1N,{ei,j}i,j∈{1,…,N}), employs the sample points pi as nodes (vertices) of the graph and connects them with weighted arcs (edges) {ei,j}i,j∈{1,…,N}.AkNNGN(PN) is built by employing a distance function, which commonly is the Euclidean one, to weight the arcs connecting each pi to its kNNs.AMST(PN) is the spanning tree minimizing the sum of the edge weights. When the weights approximate Geodesic distances [56], a GMSTN(PN) is obtained.ASIGN(PN) is defined by connecting nodes pi and pj iff pi-pj≤ρ(i)+ρ(j), where ρ(i) is the distance between pi and its nearest neighbor in PN. Essentially, two vertices pi and pj are connected if the corresponding NN hyperspheres intersect. A generalization of SIGN(PN) is kSIG(PN), where nodes pi and pj are connected iff pi-pj≤ρk(i)+ρk(j), ρk(i) being the distance between pi and its kNN in PN. This means that the kNN hyperspheres centered on pi and pj intersect.In the following we recall interestingid estimators based on GMST(PN), kNNG(PN), and kSIG(PN).In [29, 30], after defining the length functional L(GN(PN))=∑‍ei,jγ, γ∈(0,d), to build either the GMST(PN) or the MST(PN) of kNNG(PN), graph theories are exploited to estimate both the id of the underlying manifold structure M and its intrinsic Rènyi α-entropy HM. To this aim, the authors derive the linear model: logLMSTPN=alog⁡d+b, a=(d-γ)/d, b=logc+HM, c being an unknown constant, and exploit it to define an estimator of both d and HM. Briefly, a set of cardinalities {nk}k=1K is chosen and, for each nk, the MST(Pnk) is constructed on the set Pnk, which contains nk points randomly sampled from PN, to obtain a set of K pairs logLMST(Pnk),nk. Fitting them with a least squares procedure the estimates a^≃a and b^≃b are computed. Recalling that a=(d-γ)/d, the id is calculated as d^=round{γ/(1-a^)}≃d. This process is iterated to produce the final estimate as the average of the obtained results.The aforementionedkNNG based algorithm [29, 30] is exploited in [112], where the authors consider datasets sampled from a union of disjoint manifolds with possibly different ids. To estimate the local ids, the authors propose a heuristic, which is not described here, to automatically determine the local neighborhoods with similar geometric structures without any prior knowledge on the number of manifolds, their ids, and their sampling distributions.In [96] the authors present three id estimation approaches, defined as “graph theoretic methods” since the statistics they compute are functions only of graph properties (such as vertex degrees and vertex eccentricities) and do not directly depend on the interpoint distances.The first statistic, denoted bySN1(PN)=r¯j(kNNG(PN)) in the following, is based on the reach (the reach rj,i(pi,G), in j steps of a node pi∈G, is the total number of vertices which are connected to pi by a path composed of j arcs or less in G) of vertices in the kNNG(PN). Considering that the reach of each vertex pi∈kNNG(PN) grows as the id increases, in [115] the average reach r¯j(kNNG) in j steps of vertices in kNNG(PN) is employed: SN1(PN)=r¯j(kNNG(PN))=(1/N)∑i=1Nrj,i(pi,kNNG(PN)).The second statistic, denoted bySN2(PN)=MN(MST(PN)), is computed by considering the degree of vertices in the MST(PN). Recalling that, for datasets PN obtained from a continuous distribution on Rd, the ratio of nodes with a given degree j in MSTN(PN) converges a.s. to a limit depending only on j and d [116] and that the average degree in a tree is a constant depending only on the number of vertices, the authors empirically observe a dependency between the average degree and the id. This leads to the definition of an id estimator employing the statistic SN2=MN(MST(PN))=(1/N)∑i=1N(degMST(PN)(pi))2.The third statistic, denoted bySN3(PN)=UNk(kSIG(PN)), is motivated by studies in the literature [117] showing that the expected number of neighbors shared by a given pair of points depends on the id of the underlying manifold. Accordingly, calling Ni,j the number of samples in the intersection of the two kNN hyperspheres centered on pi and pj, intuitions similar to those considered for r¯j(kNNG) lead to defining of SN3(PN)=UNk(kSIG(PN))=(1/n)∑i≤j‍Ni,j.Based on their theoretical results and empirical tests on synthetically generated datasets characterized byid values dj in a finite range F⊆N+ (where F={dj}dj=212 in the reported experiments), the authors propose an approximate Bayesian estimator that could indistinctly employ each of the three statistics SN1, SN2, and SN3, denoted by SN∗ in the following. To this aim, they assume that each statistic can be approximated by a Gaussian density fdj(·)=N(μ(dj),σ2(dj)); to estimate μ(dj) and σ2(dj), for each dj∈F, L datasets of large size are synthetically generated by random sampling from the uniform distribution on the unit dj-cube. These datasets are then used to estimate the parameters μ~(dj)≃μ(dj) and σ~2(dj)≃σ2(dj) that define the approximation f~dj(·), computed on a generic sample set with size N and id=dj, of the Gaussian density fdj(·) of SN∗.At this stage, given a new input datasetPN having unknown id, the statistic SN∗(PN)=sP is computed and used to calculate the approximated value f~dj(sP)=N(μ~2(dj), σ~2(dj)/N)≃fdj(sP). Assuming equal a priori probability for all the dj∈F, the posterior probability P[dj∣SN∗] is given by(20)PdjSN∗=f~djsP∑dj∈Ff~djsP,dj∈F,and employed to compute an “a posteriori expected value” of the id: (21)d^=round∑dj∈FdjPdjSN∗.The authors evaluate the performance of their methods on synthetic datasets, some of which have been used by similar studies in the literature [79], while the others (challenging ones) are proposed by the authors to have manifolds with nonconstant curvature. The comparison of the achieved results with those obtained by the estimators proposed in [27, 33, 112, 118] has led to the conclusion that none of the methods has a good performance on all the tested datasets. However, graph theoretic approaches would appear to behave better when manifolds of nonconstant curvature are processed.This interesting comparison strengthens the need of defining a benchmark framework to allow an objective and reproducible comparative evaluation ofid estimators. For this reason, in Section 4 we describe our proposal in this direction.


## 3.1. Projectiveid Estimators

The first projectiveid estimators introduced in the literature explicitly compute the mapping that projects input points PN∈RD to the subspace M⊆Rd minimizing the information loss [27, 43] and therefore view the id as the minimal number of vectors linearly spanning the subspace M. It must be noted that, since these methods were originally designed for exploratory data analysis and dimensionality reduction, they often require the dimensionality of M (the id to be estimated) to be provided as input parameter. However, their extensions and variants include methodologies to automatically estimate id.Most of the projectiveid estimators can be grouped into two main categories: projection techniques based on Multidimensional Scaling (MDS) [44, 45] or its variants, which tend to preserve as much as possible pairwise distances among the data, and projection techniques based on Principal Component Analysis (PCA) [15, 46] and its variants that search for the best projection subspace M minimizing the projection error.Some of the best known examples ofMDS algorithms are MDSCAL [47–52], Bennett’s algorithm [1, 53], Sammon’s mapping [54], Curvilinear Component Analysis (CCA) [55], ISOMAP [56], and Local Linear Embedding (LLE) [57]. As shown by experiments reported in [56, 57] ISOMAP and variants of LLE compute the most reliable id estimates. We believe that their better performance is due to the fact that both ISOMAP and LLE have been the first projective methods based on the assumption that the input points are drawn from an underlying manifold, whose curvature might affect the precision of data neighborhoods computed by employing the Euclidean distance. However, as noted in [58, 59], these algorithms have shown that they suffer of all the major drawbacks affecting MDS-based algorithms, which are too much tied by the preservation of the pairwise distance values. Besides, as highlighted in [30], ISOMAP as an id estimator, as well as other spectral based methods like PCA, relies on a specific estimated eigenstructure that may not exist in real data. Regarding LLE, it either requires the id value to be known in advance or may automatically estimate it by analyzing the eigenvalues of the data neighborhoods [60]; however, as outlined in [7, 27], id estimates computed by means of eigenvalue analysis are as unreliable as those computed by most PCA-based approaches. Moreover, in [46] it is noted that methods such as LLE are based on the solution of a sparse eigenvalue problem under the unit covariance constraint; however, due to this imposed constraint, the global shape of the embedded data cannot reflect the underlying manifold.P
C
A [15, 46] is one of the most cited and well-known projective id estimators, often used as the first step of several pattern recognition problems, to compute low dimensional representations of the available datasets. When PCA is used for id estimation, the estimate is the number of “most relevant” eigenvectors of the sample covariance matrix, also called principal components (PCs). Due to the promising dimensionality reduction results, several PCA-based approaches, both deterministic and probabilistic, have been published. Among deterministic approaches, we recall the Kernel PCA (KPCA) [61] and the local PCA (LPCA) [62] and its extensions to automatically select the number of PCs [63, 64]. We observe that the work presented in [64] is one of the first works that estimates id by considering an underlying topological structure and therefore applies LPCA on data neighborhoods represented by an Optimally Topology Preserving Map (OTPM) built on clustered data (given an input dataset PN, its OTPM is usually computed through Topology Representing Networks (TRNs); these are unsupervised neural networks [65] developed to map PN to a set of neurons whose learnt connections define proximities in PN. These proximities correspond to the optimal topology preserving Voronoi tessellation and the corresponding Delaunay triangulation. In other words, TRNs compute connectivity structures that define and perfectly preserve the topology originally present in the data, forming a discrete path-preserving representation of the inner (topological) structure of the topological manifold underlying the dataset PN). The authors of this method state that their approach is more efficient and less sensitive to noise with respect to the PCA-based approaches. However, they do not show any experimental comparison and, besides, their algorithm employs critical thresholds and a data clustering technique whose result heavily influences the precision of the computed estimate [27].The usage of a probabilistic approach has been firstly introduced by Tipping and Bishop in [66]. Considering that deterministic methodologies lack an associated probabilistic model for the observed data, their Probabilistic PCA (PPCA) reformulates PCA as the maximum likelihood solution of a specific latent variable model. PPCA and its extensions to both mixture and hierarchical mixture models have been successfully applied to several real problems, but they still provide an id-estimation mechanism depending on critical thresholds. This motivates its subsequent variants [67] and developments, whose examples are Bayesian PCA (BPCA) [68] and two Bayesian model order selection methods introduced in [69, 70]. In [71] the asymptotic consistency of id estimation by (constrained) isotropic version of PPCA is shown with numerical experiments on simulated and real datasets.While the aforementioned methods have been simply recalled since theirid estimation results have shown to be unreliable [7, 27], in the following recent and promising proposals are described with more details.The Simple Exponential FamilyPCA (SePCA) [72] has been developed to overcome the assumption of Gaussian-distributed data that makes it difficult to handle all types of practical observations, for example, integers and binary values. SePCA achieves promising results by using exponential family distributions; however, it is highly influenced by critical parameter settings and it is successful only if the data distribution is known, which is often not the case, specially when highly nonlinear manifold structures must be treated.In [73] the authors propose the Sparse Probability PCA (SPPCA) as a probabilistic version of the Sparse PCA (SPCA) [74]. Precisely, SPCA selects id by forcing the sparsity of the projection matrix that is the matrix containing the PCs. However, based on the consideration that the level of sparsity is not automatically determined by SPCA, SPPCA employs a Bayesian formulation of SPCA, achieving sparsity by employing a different prior and automatically learning the hyperparameter related to the constraint weight through Evidence Approximation ([75] Section 3.5). The authors’ results and also the results of the comparative evaluation proposed in [76] show that this method seems to be less affected by the problems of the aforementioned projective schemes.An alternative method (MLSVD) [77] applies Singular Value Decomposition (SVD), basically a variant of PCA, locally and in a multiscale fashion to estimate the id characterizing D-dimensional datasets drawn from nonlinearly embedded d-dimensional manifolds M corrupted by Gaussian noise. Precisely, exploiting the same ideas of the theoretical PCA-based id estimator presented in [78], the authors note that the best way to avoid the effects of the curvature (induced by the nonlinearity of the embedding) is to apply SVD locally, that is, in hyperspheres B(p,r) centered on the data points p and having radius r. However, the choice of  r is constrained by the following considerations: (1) r must be big enough to have at least k≥d neighbors, (2) r must be small enough to ensure that M∩B is linear (or at least smooth), and (3) r must be big enough to ensure that the effects of noise are negligible. When these three constraints are met, the tangent space TMd(p,r), computed by applying SVD on the k neighbors, is a good approximation of the tangent space of M∩B and the number of its relevant eigenvalues corresponds to the (local) id of M. To find a proper value for r, the authors propose a multiscale approach that applies SVD on neighborhoods B(p,rs) whose radius varies in a range rs∈{rL,…,rH}. This allows us to compute D scale-dependent, local singular values λ1p,rs≥⋯≥λD(p,rs); using a least squares fitting procedure the SVs can be expressed as functions of r whose analysis allows us to identify the range of scales [rmin,…,rmax] not influenced by either noise or curvature. Finally, in the range rs=[rmin,…,rmax] the squared SVs are analyzed to get the id estimate d^ that maximizes the gap Δ(j)=λj(p,rs)-λj+1(p,rs) for the largest range of rs. The proposed algorithm has been evaluated on unit d-dimensional hyperspheres and cubes embedded in R100 and affected by Gaussian noise. The reported results are very good, while other ten well-known methods [19, 23, 27, 30, 79–81] show that the ids estimated on the same datasets are unreliable also in the absence of noise.


## 3.2. Topological Approaches

Topological approaches forid estimation consider a manifold M⊆Rd embedded in a higher dimensional space RD through a proper (locally) smooth map ϕ:M→RD and assume that the given dataset is PN={pi}i=1N={ϕ(xi)}i=1N⊂RD, where xi are independent identically distributed (i.i.d.) points drawn from M through a smooth probability density function (pdf) f:M→R+.Under this assumption theid to be estimated is the manifold’s topological dimension, defined either through the firstly proposed Brouwer Large Inductive Dimension [82] or the equivalent Lebesgue Covering Dimension [83]. Since Brouwer’s definition has been soon neglected by mathematicians for its difficult proof [83], the commonly adopted topological dimension definition is Lebesgue’s Covering Dimension, reported in the following.Definition 1 (cover).
Given a topological spaceX, a cover of a set Y⊆X is a countable collection C={Ci} of open sets such that each Ci⊂X and ⋃iCi⊇Y.Definition 2 (refinement of a cover).
A refinement of a coverC of a set Y is another cover C′ such that each set in C′ is contained in some sets of C.Definition 3 (topological dimension (Lebesgue Covering Dimension)).
Given the aforementioned definitions, the topological dimension of the topological spaceX, also called Lebesgue Covering Dimension, is d if every finite cover of X admits a refinement C′ such that no subset of X has more than d+1 intersecting open sets in C′. If no such minimal integer value exists, X is said to be of infinite topological dimension.To our knowledge, at the state of the art only two estimators have been explicitly designed to estimate the topological dimension.One of them, the Tensor Voting Framework (TVF) [84] and its variants [85], relies on the usage of an iterative information diffusion process based on Gestalt principles of perceptual organization [86]. TVF iteratively diffuses local information describing, for each pi∈PN, the tangent space approximating the underlying neighborhood of M. To this aim, the information diffused at each iteration is second-order symmetric positive definite tensors whose eigenvectors span the local tangent space. Practically, during the initialization step a ball tensor Ti0, which is an identity matrix representing the absence of orientation, is used to initialize a token pi for each point pi as {pi=(pi,Ti0)}i=1N. During iteration t each token pi “generates” the set of tensors Ti,jtj≠i that enact as votes cast to neighboring tokens; precisely, Ti,jt is sent to the jth neighbor, and it encodes information related to the local tangent space estimate in pi at time t and decays as the curvature and the distance from the jth neighbor increase. On the other side, at iteration t each token pj receives votes that are summed to update the pj’s tensor as Tjt+1=∑i≠j‍Ti,jt; this essentially refines the estimate of the local tangent space in pj. Based on the definition of topological dimension provided by Brouwer [82], in [87] it is noted that TVF can be employed to estimate the local ids by identifying the number of most relevant eigenvalues of the computed second-order tensors. Although interesting, this method has a too high computational cost, which makes it unfeasible for spaces of dimension D≥4.From the definition of Lebesgue Covering Dimension it can be derived [88] that the topological dimension of any M⊆Rd coincides with the affine dimension d of a finite simplicial complex (a simplicial complex in Rd has affine dimension d if it is a collection of affine simplexes in Rd, having at most dimension d or having at most d+1 vertices) covering M. This essentially means that a d-dimensional manifold could be approximated by a collection of d-dimensional simplexes (each having at most d+1 vertices); therefore, the topological dimension of M could be practically estimated by analyzing the number of vertices of the collection of simplexes estimated on PN. To this aim, in [89] a method is proposed to find the number of relevant positive coefficients that are needed to reconstruct each pi∈PN from a linear combination of its k neighbors, where k is a parameter to be manually set in the range d<k≤D+1. This algorithm is based on the fact that neighbors with positive reconstruction coefficients are the vertices of a simplex with dimension equal to the topological dimension of M. Practically, to ensure that k>d, its value is set to D, the reconstruction coefficients are calculated by means of an optimization problem constrained to be nonnegative, and the coefficients bigger than a user-defined threshold are considered as the relevant ones. The id estimate is then computed by employing two alternative approaches: the first one simply computes the mode of the number of relevant coefficients for each neighborhood and the second one sorts in descending order the coefficients computed for each neighborhood, computes the mean c¯ of the sorted coefficients, and estimates id as the number of relevant values in c¯. Note that, since k>d, this method is strongly affected by the curvature of the manifold when the id is big enough. Indeed, the results of the reported experimental evaluation make the authors assert that the method works well only on noisy-free data of low id (id≤6), under the assumption that the sampling process is uniform and the data points are sufficient.Though interesting, both approaches have shown to be effective only for manifolds of low curvature as well as lowid values.In the following we survey otherid estimators employing two different estimation approaches, which allow us to categorize them. More precisely, in Section 3.2.1fractalidestimators are described, which estimate different fractal dimensions since they are good approximations of the topological one; in Section 3.2.2nearest neighbors-based (NN) idestimators are recalled, which are often based on the statistical analysis of the distribution of points within small neighborhoods.


### 3.2.1. Fractalid Estimators

Since topological dimension cannot be practically estimated, several authors implicitly assume thatM has a somehow fractal structure (see [90] for an exhaustive description of fractal sets) and estimate id by employing fractal dimension estimators, the most relevant of which are surveyed in this section.Very roughly, since the basis concept of all fractal dimensions is that the volume of ad-dimensional ball of radius r scales with its size s as rd [90, 91], all fractal dimension estimators are based on the idea of counting the number of observations in a neighborhood of radius r to (somehow) estimate the rate of growth of this number. If the estimated growth is rd, then the estimated fractal dimension of the data is considered to be equal to d.Note that all the derived estimators have the fundamental limitation that, in order to get an accurate estimation, the sizeN of the dataset with id  d has to satisfy the inequality proved by Eckmann and Ruelle in [92] for the correlation dimension estimator (CD [19], see below):(1)d<2log⁡1/ρ∗log⁡N,beingρ=rD≪1,12N2rDd≫1.This will lead to a large value of N, even for a dataset with lower id.Among fractal dimension estimators, one of the most cited algorithms is presented in [19] and will be referred to as CD in the following. It is an estimator of the correlation dimension (dimCorr), whose formal definition is as follows.Definition 4 (correlation dimension).
Given a finite sample setPN, let(2)CNr=2NN-1∑i=1,i<jNIr-pi-pj,where · is the Euclidean norm and I(·) is the step function used to simulate a closed ball of radius r centered on each pi (I(y)=0 if y<0, and I(y)=1 otherwise). Then, for a countable set, the correlation dimension dimCorr is defined as(3)dimCorr=limr→0limN→∞log⁡CNrlog⁡r.In practice CD computes an estimate, d^, of dimCorr by computing CN(r) for different ri and applying least squares to fit a line through the points log⁡ri;log⁡CN(ri). It has to be noted that, to produce correct id estimates, this estimator needs a very large number of data points [22], which is never available for practical applications; however, the computed unreliable estimations can be corrected by the correction method proposed in [28].The relevance of theCD estimator is shown by its several variants and extensions. An example is the work proposed in [91], where the authors propose a normalized CD estimator for binary data and achieve estimates approximating those computed by CD.SinceCD is heavily influenced by the setting of the scale parameters, in [93] the authors estimate the id by computing the expectation value of dimCorr through maximum likelihood estimate of the distribution of distances among points. The estimated d^ is computed as (4)d^=-1Q∑k=1Qrk-1,where Q is the set Q=rk∣rk<r and rk is the Euclidean distance between two generic data points and r is a real value, called cut-off radius.To develop an estimator more efficient thanCD, in [94] the authors choose a different notion of Fd, namely, the Information Dimension dimI:(5)dimI=-limδ→0∑i=1Nδprilog⁡prilog⁡δ,where N(δ) is the minimum number of δ-sized hypercubes covering a topological space and pri is the probability of finding a point in the ith hypercube. Noting that when the scale δ in (5) is big enough the different coverings used to estimate dimI could produce different values for N(δ), the authors look for the covering composed by the minimum number Nmin(δ) of nonempty sets. Similar to the CD algorithm, the id is the average slope of the curve obtained by fitting the points with coordinates log⁡δ;∑i=1Nmin(δ)prilog⁡pri.This algorithm is compared with theCD estimator, and the experimental tests show that both methods compute the same estimates. However, the achieved computation time is much lower than that of CD.Considering thatCD can severely underestimate the topological dimension if the data distribution on the manifold is nearly nonuniform, in [7] the author proposes the Packing Number (PN), a fractal dimension estimator that approximates the Capacity Dimension (dimCap). This choice is motivated by the fact that dimCap does not depend on the data distribution on the manifold and, if both dimCap and the topological dimension exist (which is certainly the case in machine learning applications), the two dimensions agree. To formally define dimCap, the ϵ-covering number N(ϵ) of a set S⊂X must be defined; N(ϵ) is the minimum number of open balls B(x0,ϵ)=x∈X:x0-x<ϵ whose union is a covering of S, where · is a distance metric. The definition of dimCap of S⊂X is based on the observation that the covering number N(ϵ) of a d-dimensional set is proportional to ϵ-d:(6)dimCap=-limϵ→0log⁡Nϵlog⁡ϵ.Since the estimation ofN(ϵ) is computationally expensive, based on the relation N(ϵ)≤NPack(ϵ)≤N(ϵ/2), the authors employ the ϵ-Packing Number NPack(ϵ), defined in [95] as the maximum cardinality of an ϵ-separated set. Employing a greedy algorithm to compute NPack(ϵ), the estimate, d^, of dimCap is computed as(7)d^ϵ1,ϵ2=-log⁡NPackϵ1-log⁡NPackϵ2log⁡ϵ1-log⁡ϵ2.To estimate d^ a greedy algorithm is used; however, as noted by the author, the dependency of d^ with respect to the order in which the points are visited by the greedy algorithm introduces a high variance. To avoid this problem, the algorithm iterates the procedure M times on random permutations of the data and considers the average as the final id estimate. The comparative evaluation with the CD estimator makes the authors assert that PN “seems more reliable if data contains noise or the distribution on the manifold is not uniform.” Unfortunately, also this method is scale-dependent.To avoid any scale-dependency in [79] the authors propose an estimator (Hein) based on the asymptotes of a smoothed version of (2), obtained by replacing the step function I(·) with a suitable kernel function. Precisely, they define(8)UN,h,d=2NN-1∑1≤i<j≤NN1hdKhpi-pjh2,where Kh is a kernel function with bandwidth h and d is the assumed dimensionality of the manifold from which the points are sampled. Note that, to guarantee the converge of (8), the bandwidth h has to fulfill the constraint limN→∞(Nhd)=∞. For this reason the authors formalize h as a function of N and, to achieve scale-independency, propose a method that estimates the id by analyzing the convergence of U(N,h,d) when varying the parameters N and d. Precisely, the dataset is subsampled to create sets of different cardinalities nsub∈Nsub={N,N/2,N/3,N/4,N/5} and the D curves whose points have coordinates (U(nsub,h(nsub),d),nsub) are considered. Employing this information the following id estimator is proposed:(9)Sloped=maxnsub∈Nsub⁡∂Unsub,hnsub,d∂n,d^=argmind∈1,…,D⁡Sloped.This work is notable since the empirical tests are performed on synthetic datasets specifically designed to study the influence of high curvature as well as noise on the proposed estimator. The usefulness of these datasets is confirmed by the fact that they have been also employed to assess several subsequent methods [76, 96].In [97] the authors present a fractal dimension estimator derived by the analysis of a vector quantizer applied to datasets PN⊆RD. Considering the codebook Y={y1,…,yk}⊂RD containing k code-vectors yi, ak-point quantizer is defined by a measurable function Qk:RD→Y, which brings each data point to one of the code-vectors in Y. This partitions the dataset into k so-calledquantizer cellsSi={pi∈PN:Qk(pi)=yi}, where log2(k) is called therate of the quantizer. Being X a random vector distributed according to a probability distribution ν, thequantization error is er(Qk∣ν)=(Eν[X-Qk(X)r])1/r, where r∈[1,∞) and · is the Euclidean norm in RD. Given the set Qk of all D-dimensional k-point quantizers, the performance achieved by an optimal k-point quantizer on X is er∗(Qk∣ν)=infQk∈Qk(er(Qk∣ν)). When the quantizer rate is high, the quantizer cells can be well approximated by D-dimensional hyperspheres with radius equal to ϵ and centered on each code-vector yi∈Y. In this case, the regularity of ν ensures that the probability of such balls is proportional to ϵ1/d, and it can be shown [98] that er∗(Qk∣ν)≈k-1/d. This is referred to as thehigh-rate approximation and motivates the definition of quantization dimension of order r:(10)drν=-limk→∞log⁡klog⁡er∗k∣ν.The theory of high-rate quantization [98] confirms that, for a regular ν supported on the manifold M, dr(ν) exists for each 1≤r≤∞ and equals the intrinsic dimension of M. Furthermore, the limit k→∞ allows us to motivate the relation between the quantization dimension and the Capacity Dimension. Indeed, according to the theory of high-rate quantization [98, 99], there exists a decreasing sequence ϵk, such that for sufficiently large values of k (i.e., in the high-rate regime that is when k→∞) the ratio -(log⁡k)/(log⁡er∗(k∣ν)) can be approximated increasingly finely, both from below and from above, by quantities converging to the common value dimCap. To practically compute an estimate of the quantization dimension, having fixed the value of r, the authors select a range k1≤k≤k2 of codebook sizes and design a set of quantizers {Qk}k=k1k2 giving good approximations e^r(k∣ν) of er∗(k∣p) over the chosen range of k. An id estimate is obtained by fitting the points with coordinates (log⁡(k);-log⁡e^r(k∣ν)) and measuring the average slope over the chosen range k. Though the authors mention that their algorithm is less affected by underestimation biases than neighborhood-based methods (see Section 3.2.2), in [23] this statement is confuted with theoretical arguments.


### 3.2.2. Nearest Neighbors-Basedid Estimators

In this section we consider estimators, referred to asNN estimators in the following, that describe data-neighborhoods’ distributions as functions of d. They usually assume that close points are uniformly drawn from small d-dimensional balls (hyperspheres) Bd(x,r) having radius (a small radius r→0∈R+ guarantees that samples included into Bd(x,r), being less influenced by the curvature induced by the map ϕ, are approximating well enough the intrinsic structure of the underlying portion of M) r→0∈R+ and being centered on x∈M.Practically, given an input datasetPN, the value of functions f(d) is computed by approximating the sampling process related to Bd through the k-nearest neighbor algorithm (kNN).AmongNNid estimators, Trunk’s method [100] is often cited as one of the first methods. It formulates the distribution function, f(d), with an ad hoc statistic based on geometric considerations concerning angles; in practice, having fixed a threshold γ and a starting value for the parameter k, it applies kNN to find the neighbors of each pi∈PN and calculates the angle νi between the (k+1)th-nearest neighbor and the subspace spanned by the k-nearest neighbors. Considering a threshold parameter γ, if (1/N)∑i=1N‍νi≤γ, then k is considered as the id estimate; otherwise, k is incremented by 1 and the process is repeated. The main limitation of this method is the difficult choice of a proper value for γ.The work presented by Pettis et al. [26] is notable since it is one of the first works providing a mathematical motivation for the use of nearest-neighbor distances.Indeed, for an i.i.d. samplePN⊆RD drawn from a density distribution p(x) in Rd, the following approximation holds:(11)kN≃pxVdrd,where k is the number of nearest neighbors to x within the hypersphere Bd(x,r) of radius r and centered on x and V(d) is the volume of the (unit d-dimensional) ball in Rd.This means that the proportion of sample points falling inBd(x,r) is roughly approximated by p(x) times the volume of Bd(x,r). Since this volume grows as rd, assuming the density p(x) to be a constant, it follows that the number of samples in Bd(x,r) is proportional to rd. From the relationship in (11), and assuming that the samples are locally uniformly distributed, the authors derive an id estimator for d:(12)d^=r¯kkr¯k+1-r¯k,where r¯k is the average of the distances from each sample point to its kth nearest neighbors; defining ri(k) as the distance between xi and its kth-nearest neighbor, r¯k is expressed as r¯k=(1/N)∑i=1N‍ri(k).Since this algorithm is limited by the choice of a suitable value for parameterk, in [63] the authors propose a variant which considers a range of neighborhood sizes [kmin,kmax]. However, in the same work the authors themselves show that this technique generally yields an underestimate of the id when its value is high.Taking into account relation (11), in [101] the number NBd of data points in Bd(x,r) is described by a polynomial f(r)=∑s=0d‍βsrs of degree d. In practice, considering pi,pk∈PN, calling rik=pi-pk the interpoint distances, and being r=mini,k=1Nrik, R a parameter adaptively estimated (to estimate R by means of PN, the radius value corresponding to the first significant peak of the histogram of the rijs is found), a set of n radius values r=rj=r+j(R-r)/nj=1n is selected and used to calculate n pairs rj,f^(rj)j=1n, where f^rj=#rik<rji,k=1N is the number of interpoint distances strictly lower than rj. To estimate the coefficients βjj=1D, the computed pairs are fit by a least squares fitting procedure that estimates exactly D+1 coefficients. Since by hypothesis the degree of f is d, the significance test described in [17] is used to estimate the degree d^ of f^, which is considered as the id estimate. The comparative evaluation of this algorithm with the well-known Maximum Likelihood Estimator (MLE) [27] and its improved version [102], both described below, has shown that it is more robust than them when dealing with high-dimensional datasets.M
L
E [27], one of the most cited estimators, treats the neighbors of each point pi∈PN as events in a Poisson process and the distance r(j)(pi) between the query point pi and its jth nearest neighbor as the event’s arrival time. Since this process depends on d, MLE estimates id by maximizing the log-likelihood of the observed process. In practice a local id estimate is computed as(13)d^pi,k=1k∑j=1klog⁡rk+1pirjpi-1.Averaging thed^(pi,k)s, the global id estimate is d^(k)=(1/N)∑i=1N‍d^(pi,k).The theoretical stability of the proposedid estimator for data living in C1 submanifold of RD, d≤D, and for data in an affine subspace of RD has been proved, respectively, in [103, 104]. Though the authors’ comparative evaluation shows the superior performance of the proposed estimator with respect to the CD estimator [19] (see Section 3.2.1) and the NN estimator [26], they further improve it by removing its dependency from the parameter k; to this end, different values for k are adopted and the computed results are averaged to obtain the final id estimate: d^=(1/t)∑k∈{k1,…,kt}‍d^(k).Considering that, in practice,MLE is highly biased for both large and small values of k, a variant of MLE is proposed in [102], where the arithmetic mean is substituted with the harmonic average, leading to the following estimator: d^(k)=(1/N)∑i=1N(1/d^pi,k)-1.Though the proposal in [102] seems to achieve more accurate results, it is based on the assumption that neighbors surrounding each pi are independent, which is clearly incorrect. To cope with this problem, in [105] an interesting regularized version of MLE applies a regularized maximum likelihood technique to distances between neighbors. The comparative evaluation with the aforementioned MLE methods [27, 102] makes the authors state that, though the method might be the first to converge to the actual estimate given enough data points, its estimation accuracy is comparable to that achieved by the competing schemes.In [59, 106] a further improvement of MLE is presented; it achieves a better performance by substituting euclidean distances with geodesic ones.Despite the good results achieved byMLE-based approaches, these techniques have shown to be affected by the curvature induced by ϕ on the manifold neighborhoods approximated by kNN. To reduce this effect, various id estimators have been proposed in [76, 107]; here, we review those achieving the most promising experimental results.In [107] the authors firstly propose a family of id estimators (MiNDML∗), which exploit the pdfg(r;k,d) describing the distance r(1)(x) between the center x of Bd(x,r), x∈M, r→0+ and its nearest neighbor. Briefly, formulating g(r;k,d) as a function of the id value d (g(r;k,d)=kdrd-1(1-rd)k-1), the id estimator is computed by a maximum likelihood approach.After noting that this algorithm is still affected by a bias causing underestimations when the dataset dimensionality becomes sufficiently high (i.e.,id≥10), the authors present theoretical considerations which relate the bias to the fact that id estimators based on nearest-neighbor distances are often founded on statistics derived under the assumption that the amount of available data is unlimited, which is never the case in practical applications. Based on these considerations, two different estimators, named MiNDKL and IDEA, are presented.M
i
N
D
K
L compares the empirical pdf of the neighborhood distances computed on the dataset (gData) with the distribution of the neighborhood distances computed from points uniformly drawn from hyperspheres of known increasing dimensionality (gSphered). The id estimate is the dimensionality that minimizes their Kullback-Leibler divergence KL(gData,gSphered), which is evaluated by means of the data-driven technique proposed in [108].I
D
E
A relies on the authors’ observation that the quantities 1-r(j)(pi)/r(k+1)(pi) are distributed according to the beta distribution β1,d with parameters 1 and d, respectively. Therefore, since Eβ1,d=m=1/1+d, a consistent id estimator d^≃d equals(14)d^=m^1-m^≃d=m1-m,wherem^=1Nk∑i=1N∑j=1krjpirk+1pi≃m.To reduce the effect of the aforementioned bias, IDEA finally applies an asymptotic correction step that, inspired by the correction method presented in [28], models the underestimation error by considering both the base algorithm and the given dataset.Motivated by the promising results achieved byMiNDKL, in [76] the authors propose its extension, namely, DANCo; it further reduces the underestimation effect by combining an estimator employing normalized nearest-neighbor distances with one employing mutual angles. More precisely, DANCo compares the statistics estimated on PN with those estimated on (uniformly drawn) synthetic datasets of known id. The comparisons are performed by two Kullback-Leibler divergences applied to the distribution of normalized nearest-neighbor distances g(r;k,d), having g(r;k,d)=kdrd-1(1-rd)k-1, and the distribution of pairwise angles q(x;ν,τ), q(x;ν,τ) being the von Mises-Fisher distribution (VMF) [109] with parameters ν and τ.Theid estimate d^ is the one minimizing the sum of the two divergences:(15)d^=argmind∈1,…,D⁡KLgData,gSphered+KLqData,qSphered.A fast implementation of DANCo (Fast-DANCo) is also developed. Comparative evaluations show that this algorithm achieves promising results (as shown in [76] and Section 4).Another work, which is notable because the authors not only prove the consistency in probability of the presented estimators but also derive upper bounds (see (19) below) on the probability of the estimation-error for finite, and large enough, values of N, is proposed in [33]. More precisely, the authors introduce two estimators by firstly defining a function η:RD×R→R+ slowly varying near the origin (see [33] for a detailed description and motivation of this assumption). The function η is then used to express the logarithm of the probability of a point p of being in the hypersphere BD(pi,r): logPp∈BD(pi,r)=logη(p,r)+dlog(r), having Pp∈BD(pi,r)=η(p,r)rd.Considering thatP(p∈B(pi,rk(pi)))≈k/n for N big enough, the authors derive the following system of equations:(16)log⁡kn≈log⁡ηpi,r+d^pilog⁡rkpi,log⁡k2n≈log⁡ηpi,r+d^pilog⁡rk/2pi,and solve it for d^(pi) to obtain a local id estimate:(17)d^pi=log⁡2log⁡rkpi/rk/2pi.The two proposed estimators are then computed either by averaging (d^avg) or by voting (d^vote):(18)d^avg=1N∑i=1Nd^pi,d^vote=argmaxd′∈N+⁡#d^pi=d′i=1N,where #condi=1N denotes the number of points pi for which cond is verified.Under differentiability assumptions on the functionη and regularity assumptions on M the authors prove the consistency in probability of their estimators and provide upper bounds (see (19)) on the probability of the estimation-error for finite, and large enough, values of N. However, the derived bounds depend on unknown universal constants c,c′,c′′>0:(19)Pd^avg≠d≤exp⁡-c′NDcdk2,Pd^vote≠d≤exp⁡-c′′Ncdk2.


## 3.2.1. Fractalid Estimators

Since topological dimension cannot be practically estimated, several authors implicitly assume thatM has a somehow fractal structure (see [90] for an exhaustive description of fractal sets) and estimate id by employing fractal dimension estimators, the most relevant of which are surveyed in this section.Very roughly, since the basis concept of all fractal dimensions is that the volume of ad-dimensional ball of radius r scales with its size s as rd [90, 91], all fractal dimension estimators are based on the idea of counting the number of observations in a neighborhood of radius r to (somehow) estimate the rate of growth of this number. If the estimated growth is rd, then the estimated fractal dimension of the data is considered to be equal to d.Note that all the derived estimators have the fundamental limitation that, in order to get an accurate estimation, the sizeN of the dataset with id  d has to satisfy the inequality proved by Eckmann and Ruelle in [92] for the correlation dimension estimator (CD [19], see below):(1)d<2log⁡1/ρ∗log⁡N,beingρ=rD≪1,12N2rDd≫1.This will lead to a large value of N, even for a dataset with lower id.Among fractal dimension estimators, one of the most cited algorithms is presented in [19] and will be referred to as CD in the following. It is an estimator of the correlation dimension (dimCorr), whose formal definition is as follows.Definition 4 (correlation dimension).
Given a finite sample setPN, let(2)CNr=2NN-1∑i=1,i<jNIr-pi-pj,where · is the Euclidean norm and I(·) is the step function used to simulate a closed ball of radius r centered on each pi (I(y)=0 if y<0, and I(y)=1 otherwise). Then, for a countable set, the correlation dimension dimCorr is defined as(3)dimCorr=limr→0limN→∞log⁡CNrlog⁡r.In practice CD computes an estimate, d^, of dimCorr by computing CN(r) for different ri and applying least squares to fit a line through the points log⁡ri;log⁡CN(ri). It has to be noted that, to produce correct id estimates, this estimator needs a very large number of data points [22], which is never available for practical applications; however, the computed unreliable estimations can be corrected by the correction method proposed in [28].The relevance of theCD estimator is shown by its several variants and extensions. An example is the work proposed in [91], where the authors propose a normalized CD estimator for binary data and achieve estimates approximating those computed by CD.SinceCD is heavily influenced by the setting of the scale parameters, in [93] the authors estimate the id by computing the expectation value of dimCorr through maximum likelihood estimate of the distribution of distances among points. The estimated d^ is computed as (4)d^=-1Q∑k=1Qrk-1,where Q is the set Q=rk∣rk<r and rk is the Euclidean distance between two generic data points and r is a real value, called cut-off radius.To develop an estimator more efficient thanCD, in [94] the authors choose a different notion of Fd, namely, the Information Dimension dimI:(5)dimI=-limδ→0∑i=1Nδprilog⁡prilog⁡δ,where N(δ) is the minimum number of δ-sized hypercubes covering a topological space and pri is the probability of finding a point in the ith hypercube. Noting that when the scale δ in (5) is big enough the different coverings used to estimate dimI could produce different values for N(δ), the authors look for the covering composed by the minimum number Nmin(δ) of nonempty sets. Similar to the CD algorithm, the id is the average slope of the curve obtained by fitting the points with coordinates log⁡δ;∑i=1Nmin(δ)prilog⁡pri.This algorithm is compared with theCD estimator, and the experimental tests show that both methods compute the same estimates. However, the achieved computation time is much lower than that of CD.Considering thatCD can severely underestimate the topological dimension if the data distribution on the manifold is nearly nonuniform, in [7] the author proposes the Packing Number (PN), a fractal dimension estimator that approximates the Capacity Dimension (dimCap). This choice is motivated by the fact that dimCap does not depend on the data distribution on the manifold and, if both dimCap and the topological dimension exist (which is certainly the case in machine learning applications), the two dimensions agree. To formally define dimCap, the ϵ-covering number N(ϵ) of a set S⊂X must be defined; N(ϵ) is the minimum number of open balls B(x0,ϵ)=x∈X:x0-x<ϵ whose union is a covering of S, where · is a distance metric. The definition of dimCap of S⊂X is based on the observation that the covering number N(ϵ) of a d-dimensional set is proportional to ϵ-d:(6)dimCap=-limϵ→0log⁡Nϵlog⁡ϵ.Since the estimation ofN(ϵ) is computationally expensive, based on the relation N(ϵ)≤NPack(ϵ)≤N(ϵ/2), the authors employ the ϵ-Packing Number NPack(ϵ), defined in [95] as the maximum cardinality of an ϵ-separated set. Employing a greedy algorithm to compute NPack(ϵ), the estimate, d^, of dimCap is computed as(7)d^ϵ1,ϵ2=-log⁡NPackϵ1-log⁡NPackϵ2log⁡ϵ1-log⁡ϵ2.To estimate d^ a greedy algorithm is used; however, as noted by the author, the dependency of d^ with respect to the order in which the points are visited by the greedy algorithm introduces a high variance. To avoid this problem, the algorithm iterates the procedure M times on random permutations of the data and considers the average as the final id estimate. The comparative evaluation with the CD estimator makes the authors assert that PN “seems more reliable if data contains noise or the distribution on the manifold is not uniform.” Unfortunately, also this method is scale-dependent.To avoid any scale-dependency in [79] the authors propose an estimator (Hein) based on the asymptotes of a smoothed version of (2), obtained by replacing the step function I(·) with a suitable kernel function. Precisely, they define(8)UN,h,d=2NN-1∑1≤i<j≤NN1hdKhpi-pjh2,where Kh is a kernel function with bandwidth h and d is the assumed dimensionality of the manifold from which the points are sampled. Note that, to guarantee the converge of (8), the bandwidth h has to fulfill the constraint limN→∞(Nhd)=∞. For this reason the authors formalize h as a function of N and, to achieve scale-independency, propose a method that estimates the id by analyzing the convergence of U(N,h,d) when varying the parameters N and d. Precisely, the dataset is subsampled to create sets of different cardinalities nsub∈Nsub={N,N/2,N/3,N/4,N/5} and the D curves whose points have coordinates (U(nsub,h(nsub),d),nsub) are considered. Employing this information the following id estimator is proposed:(9)Sloped=maxnsub∈Nsub⁡∂Unsub,hnsub,d∂n,d^=argmind∈1,…,D⁡Sloped.This work is notable since the empirical tests are performed on synthetic datasets specifically designed to study the influence of high curvature as well as noise on the proposed estimator. The usefulness of these datasets is confirmed by the fact that they have been also employed to assess several subsequent methods [76, 96].In [97] the authors present a fractal dimension estimator derived by the analysis of a vector quantizer applied to datasets PN⊆RD. Considering the codebook Y={y1,…,yk}⊂RD containing k code-vectors yi, ak-point quantizer is defined by a measurable function Qk:RD→Y, which brings each data point to one of the code-vectors in Y. This partitions the dataset into k so-calledquantizer cellsSi={pi∈PN:Qk(pi)=yi}, where log2(k) is called therate of the quantizer. Being X a random vector distributed according to a probability distribution ν, thequantization error is er(Qk∣ν)=(Eν[X-Qk(X)r])1/r, where r∈[1,∞) and · is the Euclidean norm in RD. Given the set Qk of all D-dimensional k-point quantizers, the performance achieved by an optimal k-point quantizer on X is er∗(Qk∣ν)=infQk∈Qk(er(Qk∣ν)). When the quantizer rate is high, the quantizer cells can be well approximated by D-dimensional hyperspheres with radius equal to ϵ and centered on each code-vector yi∈Y. In this case, the regularity of ν ensures that the probability of such balls is proportional to ϵ1/d, and it can be shown [98] that er∗(Qk∣ν)≈k-1/d. This is referred to as thehigh-rate approximation and motivates the definition of quantization dimension of order r:(10)drν=-limk→∞log⁡klog⁡er∗k∣ν.The theory of high-rate quantization [98] confirms that, for a regular ν supported on the manifold M, dr(ν) exists for each 1≤r≤∞ and equals the intrinsic dimension of M. Furthermore, the limit k→∞ allows us to motivate the relation between the quantization dimension and the Capacity Dimension. Indeed, according to the theory of high-rate quantization [98, 99], there exists a decreasing sequence ϵk, such that for sufficiently large values of k (i.e., in the high-rate regime that is when k→∞) the ratio -(log⁡k)/(log⁡er∗(k∣ν)) can be approximated increasingly finely, both from below and from above, by quantities converging to the common value dimCap. To practically compute an estimate of the quantization dimension, having fixed the value of r, the authors select a range k1≤k≤k2 of codebook sizes and design a set of quantizers {Qk}k=k1k2 giving good approximations e^r(k∣ν) of er∗(k∣p) over the chosen range of k. An id estimate is obtained by fitting the points with coordinates (log⁡(k);-log⁡e^r(k∣ν)) and measuring the average slope over the chosen range k. Though the authors mention that their algorithm is less affected by underestimation biases than neighborhood-based methods (see Section 3.2.2), in [23] this statement is confuted with theoretical arguments.


## 3.2.2. Nearest Neighbors-Basedid Estimators

In this section we consider estimators, referred to asNN estimators in the following, that describe data-neighborhoods’ distributions as functions of d. They usually assume that close points are uniformly drawn from small d-dimensional balls (hyperspheres) Bd(x,r) having radius (a small radius r→0∈R+ guarantees that samples included into Bd(x,r), being less influenced by the curvature induced by the map ϕ, are approximating well enough the intrinsic structure of the underlying portion of M) r→0∈R+ and being centered on x∈M.Practically, given an input datasetPN, the value of functions f(d) is computed by approximating the sampling process related to Bd through the k-nearest neighbor algorithm (kNN).AmongNNid estimators, Trunk’s method [100] is often cited as one of the first methods. It formulates the distribution function, f(d), with an ad hoc statistic based on geometric considerations concerning angles; in practice, having fixed a threshold γ and a starting value for the parameter k, it applies kNN to find the neighbors of each pi∈PN and calculates the angle νi between the (k+1)th-nearest neighbor and the subspace spanned by the k-nearest neighbors. Considering a threshold parameter γ, if (1/N)∑i=1N‍νi≤γ, then k is considered as the id estimate; otherwise, k is incremented by 1 and the process is repeated. The main limitation of this method is the difficult choice of a proper value for γ.The work presented by Pettis et al. [26] is notable since it is one of the first works providing a mathematical motivation for the use of nearest-neighbor distances.Indeed, for an i.i.d. samplePN⊆RD drawn from a density distribution p(x) in Rd, the following approximation holds:(11)kN≃pxVdrd,where k is the number of nearest neighbors to x within the hypersphere Bd(x,r) of radius r and centered on x and V(d) is the volume of the (unit d-dimensional) ball in Rd.This means that the proportion of sample points falling inBd(x,r) is roughly approximated by p(x) times the volume of Bd(x,r). Since this volume grows as rd, assuming the density p(x) to be a constant, it follows that the number of samples in Bd(x,r) is proportional to rd. From the relationship in (11), and assuming that the samples are locally uniformly distributed, the authors derive an id estimator for d:(12)d^=r¯kkr¯k+1-r¯k,where r¯k is the average of the distances from each sample point to its kth nearest neighbors; defining ri(k) as the distance between xi and its kth-nearest neighbor, r¯k is expressed as r¯k=(1/N)∑i=1N‍ri(k).Since this algorithm is limited by the choice of a suitable value for parameterk, in [63] the authors propose a variant which considers a range of neighborhood sizes [kmin,kmax]. However, in the same work the authors themselves show that this technique generally yields an underestimate of the id when its value is high.Taking into account relation (11), in [101] the number NBd of data points in Bd(x,r) is described by a polynomial f(r)=∑s=0d‍βsrs of degree d. In practice, considering pi,pk∈PN, calling rik=pi-pk the interpoint distances, and being r=mini,k=1Nrik, R a parameter adaptively estimated (to estimate R by means of PN, the radius value corresponding to the first significant peak of the histogram of the rijs is found), a set of n radius values r=rj=r+j(R-r)/nj=1n is selected and used to calculate n pairs rj,f^(rj)j=1n, where f^rj=#rik<rji,k=1N is the number of interpoint distances strictly lower than rj. To estimate the coefficients βjj=1D, the computed pairs are fit by a least squares fitting procedure that estimates exactly D+1 coefficients. Since by hypothesis the degree of f is d, the significance test described in [17] is used to estimate the degree d^ of f^, which is considered as the id estimate. The comparative evaluation of this algorithm with the well-known Maximum Likelihood Estimator (MLE) [27] and its improved version [102], both described below, has shown that it is more robust than them when dealing with high-dimensional datasets.M
L
E [27], one of the most cited estimators, treats the neighbors of each point pi∈PN as events in a Poisson process and the distance r(j)(pi) between the query point pi and its jth nearest neighbor as the event’s arrival time. Since this process depends on d, MLE estimates id by maximizing the log-likelihood of the observed process. In practice a local id estimate is computed as(13)d^pi,k=1k∑j=1klog⁡rk+1pirjpi-1.Averaging thed^(pi,k)s, the global id estimate is d^(k)=(1/N)∑i=1N‍d^(pi,k).The theoretical stability of the proposedid estimator for data living in C1 submanifold of RD, d≤D, and for data in an affine subspace of RD has been proved, respectively, in [103, 104]. Though the authors’ comparative evaluation shows the superior performance of the proposed estimator with respect to the CD estimator [19] (see Section 3.2.1) and the NN estimator [26], they further improve it by removing its dependency from the parameter k; to this end, different values for k are adopted and the computed results are averaged to obtain the final id estimate: d^=(1/t)∑k∈{k1,…,kt}‍d^(k).Considering that, in practice,MLE is highly biased for both large and small values of k, a variant of MLE is proposed in [102], where the arithmetic mean is substituted with the harmonic average, leading to the following estimator: d^(k)=(1/N)∑i=1N(1/d^pi,k)-1.Though the proposal in [102] seems to achieve more accurate results, it is based on the assumption that neighbors surrounding each pi are independent, which is clearly incorrect. To cope with this problem, in [105] an interesting regularized version of MLE applies a regularized maximum likelihood technique to distances between neighbors. The comparative evaluation with the aforementioned MLE methods [27, 102] makes the authors state that, though the method might be the first to converge to the actual estimate given enough data points, its estimation accuracy is comparable to that achieved by the competing schemes.In [59, 106] a further improvement of MLE is presented; it achieves a better performance by substituting euclidean distances with geodesic ones.Despite the good results achieved byMLE-based approaches, these techniques have shown to be affected by the curvature induced by ϕ on the manifold neighborhoods approximated by kNN. To reduce this effect, various id estimators have been proposed in [76, 107]; here, we review those achieving the most promising experimental results.In [107] the authors firstly propose a family of id estimators (MiNDML∗), which exploit the pdfg(r;k,d) describing the distance r(1)(x) between the center x of Bd(x,r), x∈M, r→0+ and its nearest neighbor. Briefly, formulating g(r;k,d) as a function of the id value d (g(r;k,d)=kdrd-1(1-rd)k-1), the id estimator is computed by a maximum likelihood approach.After noting that this algorithm is still affected by a bias causing underestimations when the dataset dimensionality becomes sufficiently high (i.e.,id≥10), the authors present theoretical considerations which relate the bias to the fact that id estimators based on nearest-neighbor distances are often founded on statistics derived under the assumption that the amount of available data is unlimited, which is never the case in practical applications. Based on these considerations, two different estimators, named MiNDKL and IDEA, are presented.M
i
N
D
K
L compares the empirical pdf of the neighborhood distances computed on the dataset (gData) with the distribution of the neighborhood distances computed from points uniformly drawn from hyperspheres of known increasing dimensionality (gSphered). The id estimate is the dimensionality that minimizes their Kullback-Leibler divergence KL(gData,gSphered), which is evaluated by means of the data-driven technique proposed in [108].I
D
E
A relies on the authors’ observation that the quantities 1-r(j)(pi)/r(k+1)(pi) are distributed according to the beta distribution β1,d with parameters 1 and d, respectively. Therefore, since Eβ1,d=m=1/1+d, a consistent id estimator d^≃d equals(14)d^=m^1-m^≃d=m1-m,wherem^=1Nk∑i=1N∑j=1krjpirk+1pi≃m.To reduce the effect of the aforementioned bias, IDEA finally applies an asymptotic correction step that, inspired by the correction method presented in [28], models the underestimation error by considering both the base algorithm and the given dataset.Motivated by the promising results achieved byMiNDKL, in [76] the authors propose its extension, namely, DANCo; it further reduces the underestimation effect by combining an estimator employing normalized nearest-neighbor distances with one employing mutual angles. More precisely, DANCo compares the statistics estimated on PN with those estimated on (uniformly drawn) synthetic datasets of known id. The comparisons are performed by two Kullback-Leibler divergences applied to the distribution of normalized nearest-neighbor distances g(r;k,d), having g(r;k,d)=kdrd-1(1-rd)k-1, and the distribution of pairwise angles q(x;ν,τ), q(x;ν,τ) being the von Mises-Fisher distribution (VMF) [109] with parameters ν and τ.Theid estimate d^ is the one minimizing the sum of the two divergences:(15)d^=argmind∈1,…,D⁡KLgData,gSphered+KLqData,qSphered.A fast implementation of DANCo (Fast-DANCo) is also developed. Comparative evaluations show that this algorithm achieves promising results (as shown in [76] and Section 4).Another work, which is notable because the authors not only prove the consistency in probability of the presented estimators but also derive upper bounds (see (19) below) on the probability of the estimation-error for finite, and large enough, values of N, is proposed in [33]. More precisely, the authors introduce two estimators by firstly defining a function η:RD×R→R+ slowly varying near the origin (see [33] for a detailed description and motivation of this assumption). The function η is then used to express the logarithm of the probability of a point p of being in the hypersphere BD(pi,r): logPp∈BD(pi,r)=logη(p,r)+dlog(r), having Pp∈BD(pi,r)=η(p,r)rd.Considering thatP(p∈B(pi,rk(pi)))≈k/n for N big enough, the authors derive the following system of equations:(16)log⁡kn≈log⁡ηpi,r+d^pilog⁡rkpi,log⁡k2n≈log⁡ηpi,r+d^pilog⁡rk/2pi,and solve it for d^(pi) to obtain a local id estimate:(17)d^pi=log⁡2log⁡rkpi/rk/2pi.The two proposed estimators are then computed either by averaging (d^avg) or by voting (d^vote):(18)d^avg=1N∑i=1Nd^pi,d^vote=argmaxd′∈N+⁡#d^pi=d′i=1N,where #condi=1N denotes the number of points pi for which cond is verified.Under differentiability assumptions on the functionη and regularity assumptions on M the authors prove the consistency in probability of their estimators and provide upper bounds (see (19)) on the probability of the estimation-error for finite, and large enough, values of N. However, the derived bounds depend on unknown universal constants c,c′,c′′>0:(19)Pd^avg≠d≤exp⁡-c′NDcdk2,Pd^vote≠d≤exp⁡-c′′Ncdk2.


## 3.3. Graph-Basedid Estimators

As noted in [96], the work of [110] has cleared up the fact that theories underlying graphs can be applied to solve a variety of statistical problems; indeed, also in the field of id estimation various types of graph structures have been proposed [29, 96, 111, 112] and used for id estimation. Examples are the kNN graph (kNNG) [30], the Minimum Spanning Tree (MST) [113] and its variation, the geodesic MST (GMST) [29]; and the sphere of influence graph (SIG) [114] and its generalization, the k-sphere of influence graph (kSIG) [96].Given a sample setPN={pi}i=1N a graph built on PN, usually denoted by G(PN)=({pi}i=1N,{ei,j}i,j∈{1,…,N}), employs the sample points pi as nodes (vertices) of the graph and connects them with weighted arcs (edges) {ei,j}i,j∈{1,…,N}.AkNNGN(PN) is built by employing a distance function, which commonly is the Euclidean one, to weight the arcs connecting each pi to its kNNs.AMST(PN) is the spanning tree minimizing the sum of the edge weights. When the weights approximate Geodesic distances [56], a GMSTN(PN) is obtained.ASIGN(PN) is defined by connecting nodes pi and pj iff pi-pj≤ρ(i)+ρ(j), where ρ(i) is the distance between pi and its nearest neighbor in PN. Essentially, two vertices pi and pj are connected if the corresponding NN hyperspheres intersect. A generalization of SIGN(PN) is kSIG(PN), where nodes pi and pj are connected iff pi-pj≤ρk(i)+ρk(j), ρk(i) being the distance between pi and its kNN in PN. This means that the kNN hyperspheres centered on pi and pj intersect.In the following we recall interestingid estimators based on GMST(PN), kNNG(PN), and kSIG(PN).In [29, 30], after defining the length functional L(GN(PN))=∑‍ei,jγ, γ∈(0,d), to build either the GMST(PN) or the MST(PN) of kNNG(PN), graph theories are exploited to estimate both the id of the underlying manifold structure M and its intrinsic Rènyi α-entropy HM. To this aim, the authors derive the linear model: logLMSTPN=alog⁡d+b, a=(d-γ)/d, b=logc+HM, c being an unknown constant, and exploit it to define an estimator of both d and HM. Briefly, a set of cardinalities {nk}k=1K is chosen and, for each nk, the MST(Pnk) is constructed on the set Pnk, which contains nk points randomly sampled from PN, to obtain a set of K pairs logLMST(Pnk),nk. Fitting them with a least squares procedure the estimates a^≃a and b^≃b are computed. Recalling that a=(d-γ)/d, the id is calculated as d^=round{γ/(1-a^)}≃d. This process is iterated to produce the final estimate as the average of the obtained results.The aforementionedkNNG based algorithm [29, 30] is exploited in [112], where the authors consider datasets sampled from a union of disjoint manifolds with possibly different ids. To estimate the local ids, the authors propose a heuristic, which is not described here, to automatically determine the local neighborhoods with similar geometric structures without any prior knowledge on the number of manifolds, their ids, and their sampling distributions.In [96] the authors present three id estimation approaches, defined as “graph theoretic methods” since the statistics they compute are functions only of graph properties (such as vertex degrees and vertex eccentricities) and do not directly depend on the interpoint distances.The first statistic, denoted bySN1(PN)=r¯j(kNNG(PN)) in the following, is based on the reach (the reach rj,i(pi,G), in j steps of a node pi∈G, is the total number of vertices which are connected to pi by a path composed of j arcs or less in G) of vertices in the kNNG(PN). Considering that the reach of each vertex pi∈kNNG(PN) grows as the id increases, in [115] the average reach r¯j(kNNG) in j steps of vertices in kNNG(PN) is employed: SN1(PN)=r¯j(kNNG(PN))=(1/N)∑i=1Nrj,i(pi,kNNG(PN)).The second statistic, denoted bySN2(PN)=MN(MST(PN)), is computed by considering the degree of vertices in the MST(PN). Recalling that, for datasets PN obtained from a continuous distribution on Rd, the ratio of nodes with a given degree j in MSTN(PN) converges a.s. to a limit depending only on j and d [116] and that the average degree in a tree is a constant depending only on the number of vertices, the authors empirically observe a dependency between the average degree and the id. This leads to the definition of an id estimator employing the statistic SN2=MN(MST(PN))=(1/N)∑i=1N(degMST(PN)(pi))2.The third statistic, denoted bySN3(PN)=UNk(kSIG(PN)), is motivated by studies in the literature [117] showing that the expected number of neighbors shared by a given pair of points depends on the id of the underlying manifold. Accordingly, calling Ni,j the number of samples in the intersection of the two kNN hyperspheres centered on pi and pj, intuitions similar to those considered for r¯j(kNNG) lead to defining of SN3(PN)=UNk(kSIG(PN))=(1/n)∑i≤j‍Ni,j.Based on their theoretical results and empirical tests on synthetically generated datasets characterized byid values dj in a finite range F⊆N+ (where F={dj}dj=212 in the reported experiments), the authors propose an approximate Bayesian estimator that could indistinctly employ each of the three statistics SN1, SN2, and SN3, denoted by SN∗ in the following. To this aim, they assume that each statistic can be approximated by a Gaussian density fdj(·)=N(μ(dj),σ2(dj)); to estimate μ(dj) and σ2(dj), for each dj∈F, L datasets of large size are synthetically generated by random sampling from the uniform distribution on the unit dj-cube. These datasets are then used to estimate the parameters μ~(dj)≃μ(dj) and σ~2(dj)≃σ2(dj) that define the approximation f~dj(·), computed on a generic sample set with size N and id=dj, of the Gaussian density fdj(·) of SN∗.At this stage, given a new input datasetPN having unknown id, the statistic SN∗(PN)=sP is computed and used to calculate the approximated value f~dj(sP)=N(μ~2(dj), σ~2(dj)/N)≃fdj(sP). Assuming equal a priori probability for all the dj∈F, the posterior probability P[dj∣SN∗] is given by(20)PdjSN∗=f~djsP∑dj∈Ff~djsP,dj∈F,and employed to compute an “a posteriori expected value” of the id: (21)d^=round∑dj∈FdjPdjSN∗.The authors evaluate the performance of their methods on synthetic datasets, some of which have been used by similar studies in the literature [79], while the others (challenging ones) are proposed by the authors to have manifolds with nonconstant curvature. The comparison of the achieved results with those obtained by the estimators proposed in [27, 33, 112, 118] has led to the conclusion that none of the methods has a good performance on all the tested datasets. However, graph theoretic approaches would appear to behave better when manifolds of nonconstant curvature are processed.This interesting comparison strengthens the need of defining a benchmark framework to allow an objective and reproducible comparative evaluation ofid estimators. For this reason, in Section 4 we describe our proposal in this direction.


## 4. A Benchmark Proposal

At the present, an objective comparison of differentid estimators is not possible due to the lack of a standardized benchmark framework; therefore, in this section, after recalling experimental datasets and evaluation procedures introduced in the literature (see Sections 4.1 and 4.2), we choose some of them to propose a benchmark framework (see Section 4.3) that allows for reproducible and comparable experimental setups. The usefulness of the proposed benchmark is then shown by employing it to compare relevant state-of-the-art id estimators whose code is publicly available (see Section 4.4).


### 4.1. Datasets

The datasets employed in the literature are both synthetically generated datasets and real ones. In the following sections we describe those we choose to use in our benchmark study.


#### 4.1.1. Synthetic Datasets

Synthetic datasets are generated by drawing samples from manifolds (M) linearly or nonlinearly embedded in higher dimensional spaces.The publicly available tool (http://www.mL.uni-saarland.de/code/IntDim/IntDim.htm) proposed by Hein and Audibert in [79] allows us to generate 13 kinds of synthetic datasets by uniformly drawing samples from 13 manifolds of known id; they are schematically described in Table 1, where they are referred to as M∗H. These manifolds are embedded in higher dimensional spaces through both linear and nonlinear maps and are characterized by different curvatures. We note that manifold M8H is particularly challenging for its high curvature; indeed, when it is used for testing, most relevant id estimators compute pronounced id overestimates (see also the results reported in [107]).Table 1
The13 types of synthetic datasets generated with the tool proposed in [79].
Dataset
Underlying   manifold name
Description
d
D
Synthetic
M
1
H
d-dimensional sphere  linearly embedded
D
-
1
User-defined
M
2
H
Affine space
3
5
M
3
H
Concentrated figure, mistakable with a3-dimensional one
4
6
M
4
H
Nonlinear manifold
4
8
M
5
H
2-dimensional helix
2
3
M
6
H
Nonlinear manifold
6
36
M
7
H
Swiss-Roll
2
3
M
8
H
Nonlinear  (highly curved) manifold
12
72
M
9
H
Affine space
D
User-defined
M
10
H
d-dimensional hypercube
D
-
1
User-defined
M
11
H
Möebius band10-times twisted
2
3
M
12
H
Isotropic multivariate Gaussian
D
User-defined
M
13
H
1-dimensional helix curve
1
User-definedAnother interesting dataset [96] is generated by sampling a d-dimensional paraboloid, MPd, nonlinearly embedded in a higher (3(d+1)) dimensional space, according to a multivariate Burr distribution with parameter α=1. Tests on this dataset are particularly challenging since the underlying manifold is characterized by a nonconstant curvature.To perform tests on datasets generated by employing a smooth nonuniformpdf, we propose the dataset Mbeta, obtained as follows: we sample N points in [0,1)10, according to a beta distribution β0.5,10 with parameters 0.5 and 10, respectively (high skewness), and store them in a matrix XN∈RN×10; multiply each point of XN (XN(i,j)) by sin⁡(cos⁡(2πXN(i,j))), thus obtaining a matrix D1∈RN×10; multiply each point of XN by cos⁡(sin⁡(2πXN(i,j))), thus obtaining another matrix D2∈RN×10; append D1 and D2 to generate a matrix D3∈R2500×20; append D3 to its duplicate to finally generate a test dataset containing N points in R40.To further test estimators’ performance on nonlinearly embedded manifolds of highid, we propose to generate two datasets, referred to as MN1 and MN2 in the following (a tool to generate the datasets sampled from d-dimensional paraboloids, the Mbeta dataset, the MN1 dataset, and the MN2 dataset, is available at http://security.di.unimi.it/~fox721/dataset_generator.m). Precisely, to generate MN1 we uniformly draw N points in [0,1]18, we transform each point by means of tan(xicos(x18-i+1)) where i=1,…,18, we obtain points in R36 by appending each transformed x to arctan(x18-i+1sin(xi)), and we duplicate the coordinates of each point to finally generate points in R72. The id of MN1 is 18, and its points are drawn from a manifold nonlinearly embedded in R72. To generate MN2 containing N points in R96, we applied the same procedure on vectors sampled in [0,1]24.


#### 4.1.2. Real Datasets

Real datasets employed in the literature generally concern problems in the fields of image analysis, signal processing, time series prediction, and biochemistry. Among them, the most known and used datasets areISOMAP face database [56], MNIST database [119], Isolet dataset [120], D2  SantaFe [121] dataset, and DSVC1 time series [21]. Recently, the Crystal Fingerprint space for the chemical compound silicon dioxide dataset has also been proposed [22].I
S
O
M
A
P face database consists in 698 gray-level images of size 64×64 depicting the face of a sculpture. This dataset has three degrees of freedom: two for the pose and one for the lighting direction (see Figure 2(a)).Figure 2
(a) Samples fromISOMAP face database. (b) Samples from digit “0” to digit “9” in MNIST database.
(a)
(b)M
N
I
S
T database consists in 70000 gray-level images of size 28×28 of hand-written digits (see Figure 2(b)). The real id of this database is not actually known, but some works [79, 122] propose similar estimates for the different digits; as an example, the proposed id values for the digit “1” are in the range {8,…,11}.I
s
o
l
e
t dataset has been generated as follows: 150 subjects spoke the name of each letter of the alphabet twice, thus producing about 52 training examples from each speaker, for a total of 7797 samples. The speakers are grouped into 5 sets of 30 speakers each and are referred to as isolet1, isolet2, isolet3, isolet4, and isolet5. The real id value characterizing this dataset is not actually known, but a study reported in [123] shows that the correct estimate could be in the range {16,…,22}.The versionD2 of SantaFe dataset is a time series of 50000 one-dimensional points having nine degrees of freedom (id=9) and being generated by a simulation of particle motion. In order to estimate the attractor dimension of this time series, it is possible to employ the method of delays described in [124], which generates D-dimensional vectors by partitioning the original dataset in blocks containing D consecutive values; as an example, by choosing D=50 a dataset containing 1000 points in R50 is obtained.D
S
V
C
1 is a time series composed by 5000 samples measured from a hardware realization of Chua’s circuit [125]. Employing the method of delays with D=20, a dataset containing 250 points in R20 is obtained. The id characterizing this dataset is ~2.26 [21].Crystal Fingerprint spaces, or Crystal Finger spaces, have been recently proposed in crystallography [22] with the aim of representing crystalline structures; these spaces are built starting from the measured distances between atoms in the crystalline structure. The theoretical id of one Crystal Finger space consists in 3Na+3 crystal degrees of freedom, where Na is the number of atoms in the crystalline unitary cell.


### 4.2. Experimental Procedures and Evaluation Measures

At the state of the art, two approaches have been mainly used to assessid estimators on datasets of known id.The first one subsamples the test dataset to obtainT subsets of fixed cardinality and computes the percentage of correct estimations. To analyze estimators’ behavior with respect to the cardinality of input datasets, this procedure may be repeated by using different cardinality values [29, 30, 79, 122], thus obtaining a distinct performance evaluation measure for each cardinality.The second approach estimates theid on T permutations of the same dataset and averages the T  id estimates to obtain the final one [27, 76, 107, 126]. This value is then compared with the real one to assess the id estimator.To also test the estimator’s robustness with respect to its parameter settings, in [27, 107, 126] the authors apply a further test, originally proposed by Levina and Bickel in [27]. Precisely, sample sets with different cardinalities are drawn from the standard Gaussian pdf in R5 and, for each set, the estimator is applied varying the values of its parameters in fixed ranges; this allows us to analyze the behavior of the id estimate as a function of both the dataset’s cardinality and the parameter settings.Note that, sinceid estimators are usually tested on different datasets to evaluate their reliability when confronted by different dataset structures and configurations, in [126] an overall evaluation measure is proposed. This indicator, called Mean Percentage Error (MPE), summarizes all the obtained results in a unique value computed as MPE=100/#M∑M(d^M-dM/dM), where #M is the number of tested datasets, d^M is the id estimated on the dataset M, and dM is the real id of M. To apply this technique to real datasets whose id belongs to the range {dmin,…,dmax}, the same authors propose to calculate the associated MPE’s term as mind∈{dmin,…,dmax}d^M-d/dM, where dM is the mean of the range.


### 4.3. Benchmark

In this section we propose an evaluation approach which can be used as a standard framework to assess estimators performance, comparing it to relevantid estimators whose code is publicly available. In this benchmark, we suggest to use the following estimators (see Section 3): Hein, MLE, kNNG, MLSVD, BPCA, CD, MiNDKL, and DANCo (the source code of the mentioned methods is available at Hein: http://www.mL.uni-saarland.de/code.shtml, MLE: http://dept.stat.lsa.umich.edu/~elevina/mledim.m, kNNG: http://web.eecs.umich.edu/~hero/IntrinsicDim/, MLSVD: http://www.math.duke.edu/~mauro/code.html#MSVD, BPCA: http://research.microsoft.com/en-us/um/cambridge/projects/infernet/blogs/bayesianpca.aspx, CD: http://cseweb.ucsd.edu/lvdmaaten/dr/download.php, MiNDKL, and DANCo: http://www.mathworks.it/matlabcentral/fileexchange/40112-intrinsic-dimensionality-estimation-techniques). Note that these estimators cover all the groups described in Section 3, that is,Projective,Fractal,Nearest-Neighbors-based, andGraph-based estimators.The benchmark is composed by following steps:(1)
Test all the considered estimators on both the synthetic and real datasets described below. We highlight that the synthetic datasets whoseid is a user-defined parameter should be created with sufficiently high id values (id≥10).(2)
Comparative evaluation steps are as follows:
(a)
compute theMPE indicator both for synthetic and real datasets,(b)
compute a ranking test with control methods; to this aim, we suggest the Friedman test with Bonferroni-Dunn post hoc analyses [127],(c)
perform the tests proposed in [27] to evaluate the robustness, with respect to different cardinalities and parameter settings.The21 synthetic datasets used in the benchmark, referred to as M∗ in the following, are listed in Table 2 with their relevant characteristics (N, d, and D). The first 15 datasets are generated with the tool proposed in [79]; they include 4 instances, M10∗, of dataset M10, which are drawn from M10H after its embedding in RD by setting D={11,18,25,71}. Note that we did not include the dataset sampled from M8H (see Table 1) since relevant and recent id estimators have similarly produced highly overestimated results when tested on it [107]. Indeed, dealing with highly curved manifolds is still a quite challenging problem in the field.Table 2
Synthetic datasets and real datasets suggested by the benchmark;N is the dataset cardinality, d is the id, and D is the embedding space dimension.
Dataset
Dataset name
N
d
D
Synthetic
M
1
2500
10
11
M
2
2500
3
5
M
3
2500
4
6
M
4
2500
4
8
M
5
2500
2
3
M
6
2500
6
36
M
7
2500
2
3
M
9
2500
20
20
M
10
a
2500
10
11
M
10
b
2500
17
18
M
10
c
2500
24
25
M
10
d
2500
70
71
M
11
2500
2
3
M
12
2500
20
20
M
13
2500
1
13
M
N
1
2500
18
72
M
N
2
2500
24
96
M
beta
2500
10
40
M
P
3
2500
3
12
M
P
6
2500
6
21
M
P
9
2500
9
30
Real
M
D
S
C
V
1
250
2.26
20
M
I
S
O
M
A
P
698
3.00
4096
M
S
a
n
t
a
F
e
1000
9.00
50
M
M
N
I
S
T
1
70000
8.00–11.00
784
M
S
i
O
2
4738
12.00
1800
M
I
s
o
l
e
t
7797
16.00–22.00
617The last six synthetic datasets areMN1, MN2, Mbeta, and 3 instances of dataset MP∗, which are sampled from paraboloids MPd whose id is, respectively, d={3,6,9}.To perform multiple tests,20 instances of each dataset have been generated, and the achieved results have been averaged.Regarding the real datasets we used theDSVC1 time series [21] (MDSVC1, id~2.26), the ISOMAP face database [56] (MISOMAP, id=3), the SantaFe dataset [121] (MSantaFe, id=9), the MNIST database [119] (MMNIST1, id∈{8,…,13}), the Isolet dataset [120] (MIsolet, id∈{16,…,22}), and the Crystal Fingerprint space for the chemical compound silicon dioxide SiO2 structure with 3 atoms (this allows us to obtain the MSiO2 dataset containing 4738 points embedded in R1800 and being characterized by an id equal to 12).To run multiple tests also onMMNIST1, MSiO2, and MIsolet, for each of them we generated 5 instances by extracting random subsets containing 2500 points each and we averaged the achieved results.Table3 summarizes the parameter values we employed for different estimators. Note that, to relax the dependency of the kNNG algorithm from the setting of its parameter k, we performed multiple runs with k1≤k≤k2 and we averaged the achieved results. Furthermore, we tested two versions of the algorithm (referred to as kNNG1 and kNNG2) obtained by varying the parameters M and N.Table 3
Parameter settings for the different estimators:k represents the number of neighbors, γ represents the edge weighting factor for kNN, M represents the number of Least Square (LS) runs, N represents the number of resampling trials per LS iteration, α and π represent the parameters (shape and rate) of the Gamma prior distributions, which describe the hyperparameters and the observation noise model of BPCA, and μ contains the mean and the precision of the Gaussian prior distribution describing the bias inserted in the inference of BPCA.
Dataset
Method
Parameters
Synthetic
MLE
k
1
=
6
,
k
2
=
20
DANCo
k
=
10
k
N
N
G
1
k
1
=
6
,
k
2
=
20
,
γ
=
1
,
M
=
1
,
N
=
10
k
N
N
G
2
k
1
=
6
,
k
2
=
20
,
γ
=
1
,
M
=
10
,
N
=
1
BPCA
iters
=
2000, α=(2.0,2.0),π=(2.0,2.0),μ=(0.0,0.01)
Hein
N
o
n
e
CD
N
o
n
e
MLSVD
N
o
n
e
M
i
N
D
K
L
k
=
10
Real
MLE
k
1
=
3
,
k
2
=
8
DANCo
k
=
5
k
N
N
G
1
k
1
=
3
,
k
2
=
8
,
γ
=
1
,
M
=
1
,
N
=
10
k
N
N
G
2
k
1
=
3
,
k
2
=
8
,
γ
=
1
,
M
=
10
,
N
=
1
BPCA
iters
=
2000, α=(2.0,2.0),π=(2.0,2.0),μ=(0.0,0.01)
Hein
N
o
n
e
CD
N
o
n
e
MLSVD
N
o
n
e
M
i
N
D
K
L
k
=
5


### 4.4. Experimental Results

Table4 summarizes the results obtained by the compared estimators on the synthetic datasets, while in Table 5 the results obtained on the real datasets are reported.Table 4
Results achieved on the synthetic datasets. The bottom row reports theMPE achieved by each algorithm; anyhow, for each test dataset the best approximation results are highlighted in boldface.
Dataset
d
M
LE
k
N
N
G
1
k
N
N
G
2
BPCA
Hein
CD
M
i
N
D
K
L
DANCo
MLSVD
M
1
10.00
9.10
9.16
9.89
5.45
9.45
9.12
10.30
10.09
10.00
M
2
3.00
2.88
2.95
3.03
3.00
3.00
2.88
3.00
3.00
3.00
M
3
4.00
3.83
3.75
3.82
4.00
4.00
3.23
4.00
4.00
2.08
M
4
4.00
3.95
4.05
4.76
4.25
4.00
3.88
4.15
4.00
8.00
M
5
2.00
1.97
1.96
2.06
2.00
2.00
1.98
2.00
2.00
2.00
M
6
6.00
6.39
6.46
11.24
12.00
5.95
5.91
6.50
7.00
12.00
M
7
2.00
1.96
1.97
2.09
2.00
2.00
1.93
2.07
2.00
2.35
M
9
20.00
14.64
15.25
10.59
13.55
15.50
13.75
19.15
19.71
20.00
M
10
a
10.00
8.26
8.62
10.21
5.20
8.90
8.09
9.85
9.86
10.00
M
10
b
17.00
12.87
13.69
15.38
9.46
13.85
12.30
16.25
16.62
17.00
M
10
c
24.00
16.96
17.67
21.42
13.3
17.95
15.58
22.55
24.28
24.00
M
10
d
70.00
36.49
39.67
40.31
71.00
38.69
31.4
64.38
70.52
70.00
M
11
2.00
2.21
1.95
2.03
1.55
2.00
2.19
2.00
2.00
1.00
M
12
20.00
15.82
16.40
24.89
13.7
15.00
11.26
19.35
19.90
20.00
M
13
1.00
1.00
0.97
1.07
5.70
1.00
1.14
1.00
1.00
1.00
M
N
1
18.00
12.25
14.26
19.8
36.00
14.10
10.40
17.76
18.76
18.00
M
N
2
24.00
14.72
17.62
26.87
48.00
17.76
12.43
23.76
25.76
24.00
M
beta
10.00
6.36
6.45
14.77
19.7
4.00
3.05
7.00
7.00
10.00
M
P
3
3.00
2.89
2.93
3.12
7.00
2.00
2.43
3.00
3.00
1.00
M
P
6
6.00
4.96
4.98
5.82
7.00
2.66
3.58
5.04
6.00
1.00
M
P
9
9.00
6.35
6.89
8.04
10.95
2.85
4.55
7.00
8.00
1.00
MPE
17.29
14.50
16.79
62.62
19.92
25.96
5.55
3.70
26.34Table 5
Results achieved on the real datasets by the compared approaches. The bottom row reports theMPE achieved by each algorithm; anyhow, for each test dataset the best approximation results are highlighted in boldface (when the real id takes values in a range, we highlighted the results that best approximate the mean value of the range).
Dataset
id
MLE
k
N
N
G
1
k
N
N
G
2
BPCA
Hein
CD
M
i
N
D
K
L
DANCo
MLSVD
M
D
S
C
V
1
2.26
2.03
1.77
1.86
6.00
3.00
1.92
2.50
2.26
1.75
M
I
S
O
M
A
P
3.00
4.05
3.60
4.32
4.00
3.00
3.37
3.9
4.00
1.00
M
S
a
n
t
a
F
e
9.00
7.16
7.28
7.43
18.00
6.00
4.39
7.60
8.19
1.00
M
M
N
I
S
T
1
8.00–11.00
10.29
10.37
9.58
11.00
8.00
6.96
11.00
9.98
1.00
M
S
i
O
2
12.00
39.28
10.24
10.36
3.00
4.80
1.05
17.20
12.60
1.00
M
I
s
o
l
e
t
16.00–22.00
15.78
6.50
8.32
19.00
3.00
3.65
20.00
19.00
1.00
MPE
53.83
27.41
26.76
71.68
34.50
43.34
27.00
15.14
75.17Looking at the number of correct estimations computed by each algorithm (highlighted in boldface), we have the following ranking:MLSVD is correct on 13 out of 21 synthetic datasets, DANCo (correct on 10 out of 21 datasets), Hein (correct on 6 out of 21), MiNDKL (6 out of 21), BPCA (4 out of 21), and MLE (1 out of 21). It can be further noted that kNNG∗, CD, MLE, and Hein obtain good estimates only for low id manifolds, while they produce underestimated values when processing datasets of high id.By observing theMPE indicator, which accounts for the precision of the achieved estimates, we obtain a different ranking: DANCo, MiNDKL, and kNNG1 and kNNG2, MLE, Hein, CD, and MLSVD. This difference is due to the fact that algorithms, such as kNNG1 and kNNG2, MLE, and Hein, most of the times produce results close to the correct value.Regarding the real datasets, all the algorithms achieve a much worseMPE indicator, and again DANCo is best performing method.Furthermore, we compute the Friedman ranking test with the Bonferroni-Dunn post hoc analysis as proposed in Section4.3 to state the quality of the achieved results on both the synthetic and real datasets. Tables 6 and 7 summarize the ranking results.Table 6
Friedman ranking results achieved on all the datasets. The null hypothesis that the algorithms perform comparably is rejected withp value <0.00001.
Method
Ranking
D
ANCo
2.40
M
i
N
D
K
L
3.46
Hein
4.67
k
N
N
G
2
5.11
MLSVD
5.17
k
N
N
G
1
5.17
MLE
5.70
CD
6.63
BPCA
6.68Table 7
Hypothesis testing of significance between techniques. Bonferroni-Dunn’s procedure rejects those hypotheses that have ap value ≤0.0125.
M
i
N
D
K
L
Hein
k
N
N
G
1
k
N
N
G
2
MLE
CD
MLSVD
BPCA
D
ANCo
0.1567
0.0024
0.0003
0.0002
0.0002
0.0000
0.0000
0.0000
M
i
N
D
K
L
∗
∗
∗
0.0801
0.0303
0.0244
0.0055
0.0020
0.0000
0.0000
Hein
∗
∗
∗
∗
∗
∗
0.7528
0.6366
0.1564
0.1474
0.0034
0.0018
k
N
N
G
1
∗
∗
∗
∗
∗
∗
∗
∗
∗
0.8557
0.3443
0.2301
0.0164
0.0071
k
N
N
G
2
∗
∗
∗
∗
∗
∗
∗
∗
∗
∗
∗
∗
0.9314
0.3894
0.1113
0.0282
MLE
∗
∗
∗
∗
∗
∗
∗
∗
∗
∗
∗
∗
∗
∗
∗
0.3428
0.1876
0.0307
CD
∗
∗
∗
∗
∗
∗
∗
∗
∗
∗
∗
∗
∗
∗
∗
∗
∗
∗
0.7337
0.1961Finally, we performed the tests proposed in [27] to evaluate the robustness of MiNDKL, MLE, DANCo, and kNNG∗ with respect to the settings of their k parameter. Precisely, these tests employ synthetic datasets subsampled from the standard Gaussian pdf in R5 (id=5). As proposed in Section 4.2, we repeated the tests for datasets with cardinalities N∈{200,500,1000,2000} varying the parameter k in the range {5,…,100}.As shown in Figure3 many of the tested techniques are strongly influenced by the parameter settings; therefore, studying the variability of the algorithms’ behavior when changing their parameter settings is of utmost importance.Figure 3
Behavior of (a)MLE, (b) DANCo, (c) kNNG1, (d) kNNG2, and (e) MiNDKL applied to points drawn from a 5-dimensional standard normal distribution; in this test N∈{200,500,1000,2000} and k∈{5,…,100}.
(a)
(b)
(c)
(d)
(e)


## 4.1. Datasets

The datasets employed in the literature are both synthetically generated datasets and real ones. In the following sections we describe those we choose to use in our benchmark study.


### 4.1.1. Synthetic Datasets

Synthetic datasets are generated by drawing samples from manifolds (M) linearly or nonlinearly embedded in higher dimensional spaces.The publicly available tool (http://www.mL.uni-saarland.de/code/IntDim/IntDim.htm) proposed by Hein and Audibert in [79] allows us to generate 13 kinds of synthetic datasets by uniformly drawing samples from 13 manifolds of known id; they are schematically described in Table 1, where they are referred to as M∗H. These manifolds are embedded in higher dimensional spaces through both linear and nonlinear maps and are characterized by different curvatures. We note that manifold M8H is particularly challenging for its high curvature; indeed, when it is used for testing, most relevant id estimators compute pronounced id overestimates (see also the results reported in [107]).Table 1
The13 types of synthetic datasets generated with the tool proposed in [79].
Dataset
Underlying   manifold name
Description
d
D
Synthetic
M
1
H
d-dimensional sphere  linearly embedded
D
-
1
User-defined
M
2
H
Affine space
3
5
M
3
H
Concentrated figure, mistakable with a3-dimensional one
4
6
M
4
H
Nonlinear manifold
4
8
M
5
H
2-dimensional helix
2
3
M
6
H
Nonlinear manifold
6
36
M
7
H
Swiss-Roll
2
3
M
8
H
Nonlinear  (highly curved) manifold
12
72
M
9
H
Affine space
D
User-defined
M
10
H
d-dimensional hypercube
D
-
1
User-defined
M
11
H
Möebius band10-times twisted
2
3
M
12
H
Isotropic multivariate Gaussian
D
User-defined
M
13
H
1-dimensional helix curve
1
User-definedAnother interesting dataset [96] is generated by sampling a d-dimensional paraboloid, MPd, nonlinearly embedded in a higher (3(d+1)) dimensional space, according to a multivariate Burr distribution with parameter α=1. Tests on this dataset are particularly challenging since the underlying manifold is characterized by a nonconstant curvature.To perform tests on datasets generated by employing a smooth nonuniformpdf, we propose the dataset Mbeta, obtained as follows: we sample N points in [0,1)10, according to a beta distribution β0.5,10 with parameters 0.5 and 10, respectively (high skewness), and store them in a matrix XN∈RN×10; multiply each point of XN (XN(i,j)) by sin⁡(cos⁡(2πXN(i,j))), thus obtaining a matrix D1∈RN×10; multiply each point of XN by cos⁡(sin⁡(2πXN(i,j))), thus obtaining another matrix D2∈RN×10; append D1 and D2 to generate a matrix D3∈R2500×20; append D3 to its duplicate to finally generate a test dataset containing N points in R40.To further test estimators’ performance on nonlinearly embedded manifolds of highid, we propose to generate two datasets, referred to as MN1 and MN2 in the following (a tool to generate the datasets sampled from d-dimensional paraboloids, the Mbeta dataset, the MN1 dataset, and the MN2 dataset, is available at http://security.di.unimi.it/~fox721/dataset_generator.m). Precisely, to generate MN1 we uniformly draw N points in [0,1]18, we transform each point by means of tan(xicos(x18-i+1)) where i=1,…,18, we obtain points in R36 by appending each transformed x to arctan(x18-i+1sin(xi)), and we duplicate the coordinates of each point to finally generate points in R72. The id of MN1 is 18, and its points are drawn from a manifold nonlinearly embedded in R72. To generate MN2 containing N points in R96, we applied the same procedure on vectors sampled in [0,1]24.


### 4.1.2. Real Datasets

Real datasets employed in the literature generally concern problems in the fields of image analysis, signal processing, time series prediction, and biochemistry. Among them, the most known and used datasets areISOMAP face database [56], MNIST database [119], Isolet dataset [120], D2  SantaFe [121] dataset, and DSVC1 time series [21]. Recently, the Crystal Fingerprint space for the chemical compound silicon dioxide dataset has also been proposed [22].I
S
O
M
A
P face database consists in 698 gray-level images of size 64×64 depicting the face of a sculpture. This dataset has three degrees of freedom: two for the pose and one for the lighting direction (see Figure 2(a)).Figure 2
(a) Samples fromISOMAP face database. (b) Samples from digit “0” to digit “9” in MNIST database.
(a)
(b)M
N
I
S
T database consists in 70000 gray-level images of size 28×28 of hand-written digits (see Figure 2(b)). The real id of this database is not actually known, but some works [79, 122] propose similar estimates for the different digits; as an example, the proposed id values for the digit “1” are in the range {8,…,11}.I
s
o
l
e
t dataset has been generated as follows: 150 subjects spoke the name of each letter of the alphabet twice, thus producing about 52 training examples from each speaker, for a total of 7797 samples. The speakers are grouped into 5 sets of 30 speakers each and are referred to as isolet1, isolet2, isolet3, isolet4, and isolet5. The real id value characterizing this dataset is not actually known, but a study reported in [123] shows that the correct estimate could be in the range {16,…,22}.The versionD2 of SantaFe dataset is a time series of 50000 one-dimensional points having nine degrees of freedom (id=9) and being generated by a simulation of particle motion. In order to estimate the attractor dimension of this time series, it is possible to employ the method of delays described in [124], which generates D-dimensional vectors by partitioning the original dataset in blocks containing D consecutive values; as an example, by choosing D=50 a dataset containing 1000 points in R50 is obtained.D
S
V
C
1 is a time series composed by 5000 samples measured from a hardware realization of Chua’s circuit [125]. Employing the method of delays with D=20, a dataset containing 250 points in R20 is obtained. The id characterizing this dataset is ~2.26 [21].Crystal Fingerprint spaces, or Crystal Finger spaces, have been recently proposed in crystallography [22] with the aim of representing crystalline structures; these spaces are built starting from the measured distances between atoms in the crystalline structure. The theoretical id of one Crystal Finger space consists in 3Na+3 crystal degrees of freedom, where Na is the number of atoms in the crystalline unitary cell.


## 4.1.1. Synthetic Datasets

Synthetic datasets are generated by drawing samples from manifolds (M) linearly or nonlinearly embedded in higher dimensional spaces.The publicly available tool (http://www.mL.uni-saarland.de/code/IntDim/IntDim.htm) proposed by Hein and Audibert in [79] allows us to generate 13 kinds of synthetic datasets by uniformly drawing samples from 13 manifolds of known id; they are schematically described in Table 1, where they are referred to as M∗H. These manifolds are embedded in higher dimensional spaces through both linear and nonlinear maps and are characterized by different curvatures. We note that manifold M8H is particularly challenging for its high curvature; indeed, when it is used for testing, most relevant id estimators compute pronounced id overestimates (see also the results reported in [107]).Table 1
The13 types of synthetic datasets generated with the tool proposed in [79].
Dataset
Underlying   manifold name
Description
d
D
Synthetic
M
1
H
d-dimensional sphere  linearly embedded
D
-
1
User-defined
M
2
H
Affine space
3
5
M
3
H
Concentrated figure, mistakable with a3-dimensional one
4
6
M
4
H
Nonlinear manifold
4
8
M
5
H
2-dimensional helix
2
3
M
6
H
Nonlinear manifold
6
36
M
7
H
Swiss-Roll
2
3
M
8
H
Nonlinear  (highly curved) manifold
12
72
M
9
H
Affine space
D
User-defined
M
10
H
d-dimensional hypercube
D
-
1
User-defined
M
11
H
Möebius band10-times twisted
2
3
M
12
H
Isotropic multivariate Gaussian
D
User-defined
M
13
H
1-dimensional helix curve
1
User-definedAnother interesting dataset [96] is generated by sampling a d-dimensional paraboloid, MPd, nonlinearly embedded in a higher (3(d+1)) dimensional space, according to a multivariate Burr distribution with parameter α=1. Tests on this dataset are particularly challenging since the underlying manifold is characterized by a nonconstant curvature.To perform tests on datasets generated by employing a smooth nonuniformpdf, we propose the dataset Mbeta, obtained as follows: we sample N points in [0,1)10, according to a beta distribution β0.5,10 with parameters 0.5 and 10, respectively (high skewness), and store them in a matrix XN∈RN×10; multiply each point of XN (XN(i,j)) by sin⁡(cos⁡(2πXN(i,j))), thus obtaining a matrix D1∈RN×10; multiply each point of XN by cos⁡(sin⁡(2πXN(i,j))), thus obtaining another matrix D2∈RN×10; append D1 and D2 to generate a matrix D3∈R2500×20; append D3 to its duplicate to finally generate a test dataset containing N points in R40.To further test estimators’ performance on nonlinearly embedded manifolds of highid, we propose to generate two datasets, referred to as MN1 and MN2 in the following (a tool to generate the datasets sampled from d-dimensional paraboloids, the Mbeta dataset, the MN1 dataset, and the MN2 dataset, is available at http://security.di.unimi.it/~fox721/dataset_generator.m). Precisely, to generate MN1 we uniformly draw N points in [0,1]18, we transform each point by means of tan(xicos(x18-i+1)) where i=1,…,18, we obtain points in R36 by appending each transformed x to arctan(x18-i+1sin(xi)), and we duplicate the coordinates of each point to finally generate points in R72. The id of MN1 is 18, and its points are drawn from a manifold nonlinearly embedded in R72. To generate MN2 containing N points in R96, we applied the same procedure on vectors sampled in [0,1]24.


## 4.1.2. Real Datasets

Real datasets employed in the literature generally concern problems in the fields of image analysis, signal processing, time series prediction, and biochemistry. Among them, the most known and used datasets areISOMAP face database [56], MNIST database [119], Isolet dataset [120], D2  SantaFe [121] dataset, and DSVC1 time series [21]. Recently, the Crystal Fingerprint space for the chemical compound silicon dioxide dataset has also been proposed [22].I
S
O
M
A
P face database consists in 698 gray-level images of size 64×64 depicting the face of a sculpture. This dataset has three degrees of freedom: two for the pose and one for the lighting direction (see Figure 2(a)).Figure 2
(a) Samples fromISOMAP face database. (b) Samples from digit “0” to digit “9” in MNIST database.
(a)
(b)M
N
I
S
T database consists in 70000 gray-level images of size 28×28 of hand-written digits (see Figure 2(b)). The real id of this database is not actually known, but some works [79, 122] propose similar estimates for the different digits; as an example, the proposed id values for the digit “1” are in the range {8,…,11}.I
s
o
l
e
t dataset has been generated as follows: 150 subjects spoke the name of each letter of the alphabet twice, thus producing about 52 training examples from each speaker, for a total of 7797 samples. The speakers are grouped into 5 sets of 30 speakers each and are referred to as isolet1, isolet2, isolet3, isolet4, and isolet5. The real id value characterizing this dataset is not actually known, but a study reported in [123] shows that the correct estimate could be in the range {16,…,22}.The versionD2 of SantaFe dataset is a time series of 50000 one-dimensional points having nine degrees of freedom (id=9) and being generated by a simulation of particle motion. In order to estimate the attractor dimension of this time series, it is possible to employ the method of delays described in [124], which generates D-dimensional vectors by partitioning the original dataset in blocks containing D consecutive values; as an example, by choosing D=50 a dataset containing 1000 points in R50 is obtained.D
S
V
C
1 is a time series composed by 5000 samples measured from a hardware realization of Chua’s circuit [125]. Employing the method of delays with D=20, a dataset containing 250 points in R20 is obtained. The id characterizing this dataset is ~2.26 [21].Crystal Fingerprint spaces, or Crystal Finger spaces, have been recently proposed in crystallography [22] with the aim of representing crystalline structures; these spaces are built starting from the measured distances between atoms in the crystalline structure. The theoretical id of one Crystal Finger space consists in 3Na+3 crystal degrees of freedom, where Na is the number of atoms in the crystalline unitary cell.


## 4.2. Experimental Procedures and Evaluation Measures

At the state of the art, two approaches have been mainly used to assessid estimators on datasets of known id.The first one subsamples the test dataset to obtainT subsets of fixed cardinality and computes the percentage of correct estimations. To analyze estimators’ behavior with respect to the cardinality of input datasets, this procedure may be repeated by using different cardinality values [29, 30, 79, 122], thus obtaining a distinct performance evaluation measure for each cardinality.The second approach estimates theid on T permutations of the same dataset and averages the T  id estimates to obtain the final one [27, 76, 107, 126]. This value is then compared with the real one to assess the id estimator.To also test the estimator’s robustness with respect to its parameter settings, in [27, 107, 126] the authors apply a further test, originally proposed by Levina and Bickel in [27]. Precisely, sample sets with different cardinalities are drawn from the standard Gaussian pdf in R5 and, for each set, the estimator is applied varying the values of its parameters in fixed ranges; this allows us to analyze the behavior of the id estimate as a function of both the dataset’s cardinality and the parameter settings.Note that, sinceid estimators are usually tested on different datasets to evaluate their reliability when confronted by different dataset structures and configurations, in [126] an overall evaluation measure is proposed. This indicator, called Mean Percentage Error (MPE), summarizes all the obtained results in a unique value computed as MPE=100/#M∑M(d^M-dM/dM), where #M is the number of tested datasets, d^M is the id estimated on the dataset M, and dM is the real id of M. To apply this technique to real datasets whose id belongs to the range {dmin,…,dmax}, the same authors propose to calculate the associated MPE’s term as mind∈{dmin,…,dmax}d^M-d/dM, where dM is the mean of the range.


## 4.3. Benchmark

In this section we propose an evaluation approach which can be used as a standard framework to assess estimators performance, comparing it to relevantid estimators whose code is publicly available. In this benchmark, we suggest to use the following estimators (see Section 3): Hein, MLE, kNNG, MLSVD, BPCA, CD, MiNDKL, and DANCo (the source code of the mentioned methods is available at Hein: http://www.mL.uni-saarland.de/code.shtml, MLE: http://dept.stat.lsa.umich.edu/~elevina/mledim.m, kNNG: http://web.eecs.umich.edu/~hero/IntrinsicDim/, MLSVD: http://www.math.duke.edu/~mauro/code.html#MSVD, BPCA: http://research.microsoft.com/en-us/um/cambridge/projects/infernet/blogs/bayesianpca.aspx, CD: http://cseweb.ucsd.edu/lvdmaaten/dr/download.php, MiNDKL, and DANCo: http://www.mathworks.it/matlabcentral/fileexchange/40112-intrinsic-dimensionality-estimation-techniques). Note that these estimators cover all the groups described in Section 3, that is,Projective,Fractal,Nearest-Neighbors-based, andGraph-based estimators.The benchmark is composed by following steps:(1)
Test all the considered estimators on both the synthetic and real datasets described below. We highlight that the synthetic datasets whoseid is a user-defined parameter should be created with sufficiently high id values (id≥10).(2)
Comparative evaluation steps are as follows:
(a)
compute theMPE indicator both for synthetic and real datasets,(b)
compute a ranking test with control methods; to this aim, we suggest the Friedman test with Bonferroni-Dunn post hoc analyses [127],(c)
perform the tests proposed in [27] to evaluate the robustness, with respect to different cardinalities and parameter settings.The21 synthetic datasets used in the benchmark, referred to as M∗ in the following, are listed in Table 2 with their relevant characteristics (N, d, and D). The first 15 datasets are generated with the tool proposed in [79]; they include 4 instances, M10∗, of dataset M10, which are drawn from M10H after its embedding in RD by setting D={11,18,25,71}. Note that we did not include the dataset sampled from M8H (see Table 1) since relevant and recent id estimators have similarly produced highly overestimated results when tested on it [107]. Indeed, dealing with highly curved manifolds is still a quite challenging problem in the field.Table 2
Synthetic datasets and real datasets suggested by the benchmark;N is the dataset cardinality, d is the id, and D is the embedding space dimension.
Dataset
Dataset name
N
d
D
Synthetic
M
1
2500
10
11
M
2
2500
3
5
M
3
2500
4
6
M
4
2500
4
8
M
5
2500
2
3
M
6
2500
6
36
M
7
2500
2
3
M
9
2500
20
20
M
10
a
2500
10
11
M
10
b
2500
17
18
M
10
c
2500
24
25
M
10
d
2500
70
71
M
11
2500
2
3
M
12
2500
20
20
M
13
2500
1
13
M
N
1
2500
18
72
M
N
2
2500
24
96
M
beta
2500
10
40
M
P
3
2500
3
12
M
P
6
2500
6
21
M
P
9
2500
9
30
Real
M
D
S
C
V
1
250
2.26
20
M
I
S
O
M
A
P
698
3.00
4096
M
S
a
n
t
a
F
e
1000
9.00
50
M
M
N
I
S
T
1
70000
8.00–11.00
784
M
S
i
O
2
4738
12.00
1800
M
I
s
o
l
e
t
7797
16.00–22.00
617The last six synthetic datasets areMN1, MN2, Mbeta, and 3 instances of dataset MP∗, which are sampled from paraboloids MPd whose id is, respectively, d={3,6,9}.To perform multiple tests,20 instances of each dataset have been generated, and the achieved results have been averaged.Regarding the real datasets we used theDSVC1 time series [21] (MDSVC1, id~2.26), the ISOMAP face database [56] (MISOMAP, id=3), the SantaFe dataset [121] (MSantaFe, id=9), the MNIST database [119] (MMNIST1, id∈{8,…,13}), the Isolet dataset [120] (MIsolet, id∈{16,…,22}), and the Crystal Fingerprint space for the chemical compound silicon dioxide SiO2 structure with 3 atoms (this allows us to obtain the MSiO2 dataset containing 4738 points embedded in R1800 and being characterized by an id equal to 12).To run multiple tests also onMMNIST1, MSiO2, and MIsolet, for each of them we generated 5 instances by extracting random subsets containing 2500 points each and we averaged the achieved results.Table3 summarizes the parameter values we employed for different estimators. Note that, to relax the dependency of the kNNG algorithm from the setting of its parameter k, we performed multiple runs with k1≤k≤k2 and we averaged the achieved results. Furthermore, we tested two versions of the algorithm (referred to as kNNG1 and kNNG2) obtained by varying the parameters M and N.Table 3
Parameter settings for the different estimators:k represents the number of neighbors, γ represents the edge weighting factor for kNN, M represents the number of Least Square (LS) runs, N represents the number of resampling trials per LS iteration, α and π represent the parameters (shape and rate) of the Gamma prior distributions, which describe the hyperparameters and the observation noise model of BPCA, and μ contains the mean and the precision of the Gaussian prior distribution describing the bias inserted in the inference of BPCA.
Dataset
Method
Parameters
Synthetic
MLE
k
1
=
6
,
k
2
=
20
DANCo
k
=
10
k
N
N
G
1
k
1
=
6
,
k
2
=
20
,
γ
=
1
,
M
=
1
,
N
=
10
k
N
N
G
2
k
1
=
6
,
k
2
=
20
,
γ
=
1
,
M
=
10
,
N
=
1
BPCA
iters
=
2000, α=(2.0,2.0),π=(2.0,2.0),μ=(0.0,0.01)
Hein
N
o
n
e
CD
N
o
n
e
MLSVD
N
o
n
e
M
i
N
D
K
L
k
=
10
Real
MLE
k
1
=
3
,
k
2
=
8
DANCo
k
=
5
k
N
N
G
1
k
1
=
3
,
k
2
=
8
,
γ
=
1
,
M
=
1
,
N
=
10
k
N
N
G
2
k
1
=
3
,
k
2
=
8
,
γ
=
1
,
M
=
10
,
N
=
1
BPCA
iters
=
2000, α=(2.0,2.0),π=(2.0,2.0),μ=(0.0,0.01)
Hein
N
o
n
e
CD
N
o
n
e
MLSVD
N
o
n
e
M
i
N
D
K
L
k
=
5


## 4.4. Experimental Results

Table4 summarizes the results obtained by the compared estimators on the synthetic datasets, while in Table 5 the results obtained on the real datasets are reported.Table 4
Results achieved on the synthetic datasets. The bottom row reports theMPE achieved by each algorithm; anyhow, for each test dataset the best approximation results are highlighted in boldface.
Dataset
d
M
LE
k
N
N
G
1
k
N
N
G
2
BPCA
Hein
CD
M
i
N
D
K
L
DANCo
MLSVD
M
1
10.00
9.10
9.16
9.89
5.45
9.45
9.12
10.30
10.09
10.00
M
2
3.00
2.88
2.95
3.03
3.00
3.00
2.88
3.00
3.00
3.00
M
3
4.00
3.83
3.75
3.82
4.00
4.00
3.23
4.00
4.00
2.08
M
4
4.00
3.95
4.05
4.76
4.25
4.00
3.88
4.15
4.00
8.00
M
5
2.00
1.97
1.96
2.06
2.00
2.00
1.98
2.00
2.00
2.00
M
6
6.00
6.39
6.46
11.24
12.00
5.95
5.91
6.50
7.00
12.00
M
7
2.00
1.96
1.97
2.09
2.00
2.00
1.93
2.07
2.00
2.35
M
9
20.00
14.64
15.25
10.59
13.55
15.50
13.75
19.15
19.71
20.00
M
10
a
10.00
8.26
8.62
10.21
5.20
8.90
8.09
9.85
9.86
10.00
M
10
b
17.00
12.87
13.69
15.38
9.46
13.85
12.30
16.25
16.62
17.00
M
10
c
24.00
16.96
17.67
21.42
13.3
17.95
15.58
22.55
24.28
24.00
M
10
d
70.00
36.49
39.67
40.31
71.00
38.69
31.4
64.38
70.52
70.00
M
11
2.00
2.21
1.95
2.03
1.55
2.00
2.19
2.00
2.00
1.00
M
12
20.00
15.82
16.40
24.89
13.7
15.00
11.26
19.35
19.90
20.00
M
13
1.00
1.00
0.97
1.07
5.70
1.00
1.14
1.00
1.00
1.00
M
N
1
18.00
12.25
14.26
19.8
36.00
14.10
10.40
17.76
18.76
18.00
M
N
2
24.00
14.72
17.62
26.87
48.00
17.76
12.43
23.76
25.76
24.00
M
beta
10.00
6.36
6.45
14.77
19.7
4.00
3.05
7.00
7.00
10.00
M
P
3
3.00
2.89
2.93
3.12
7.00
2.00
2.43
3.00
3.00
1.00
M
P
6
6.00
4.96
4.98
5.82
7.00
2.66
3.58
5.04
6.00
1.00
M
P
9
9.00
6.35
6.89
8.04
10.95
2.85
4.55
7.00
8.00
1.00
MPE
17.29
14.50
16.79
62.62
19.92
25.96
5.55
3.70
26.34Table 5
Results achieved on the real datasets by the compared approaches. The bottom row reports theMPE achieved by each algorithm; anyhow, for each test dataset the best approximation results are highlighted in boldface (when the real id takes values in a range, we highlighted the results that best approximate the mean value of the range).
Dataset
id
MLE
k
N
N
G
1
k
N
N
G
2
BPCA
Hein
CD
M
i
N
D
K
L
DANCo
MLSVD
M
D
S
C
V
1
2.26
2.03
1.77
1.86
6.00
3.00
1.92
2.50
2.26
1.75
M
I
S
O
M
A
P
3.00
4.05
3.60
4.32
4.00
3.00
3.37
3.9
4.00
1.00
M
S
a
n
t
a
F
e
9.00
7.16
7.28
7.43
18.00
6.00
4.39
7.60
8.19
1.00
M
M
N
I
S
T
1
8.00–11.00
10.29
10.37
9.58
11.00
8.00
6.96
11.00
9.98
1.00
M
S
i
O
2
12.00
39.28
10.24
10.36
3.00
4.80
1.05
17.20
12.60
1.00
M
I
s
o
l
e
t
16.00–22.00
15.78
6.50
8.32
19.00
3.00
3.65
20.00
19.00
1.00
MPE
53.83
27.41
26.76
71.68
34.50
43.34
27.00
15.14
75.17Looking at the number of correct estimations computed by each algorithm (highlighted in boldface), we have the following ranking:MLSVD is correct on 13 out of 21 synthetic datasets, DANCo (correct on 10 out of 21 datasets), Hein (correct on 6 out of 21), MiNDKL (6 out of 21), BPCA (4 out of 21), and MLE (1 out of 21). It can be further noted that kNNG∗, CD, MLE, and Hein obtain good estimates only for low id manifolds, while they produce underestimated values when processing datasets of high id.By observing theMPE indicator, which accounts for the precision of the achieved estimates, we obtain a different ranking: DANCo, MiNDKL, and kNNG1 and kNNG2, MLE, Hein, CD, and MLSVD. This difference is due to the fact that algorithms, such as kNNG1 and kNNG2, MLE, and Hein, most of the times produce results close to the correct value.Regarding the real datasets, all the algorithms achieve a much worseMPE indicator, and again DANCo is best performing method.Furthermore, we compute the Friedman ranking test with the Bonferroni-Dunn post hoc analysis as proposed in Section4.3 to state the quality of the achieved results on both the synthetic and real datasets. Tables 6 and 7 summarize the ranking results.Table 6
Friedman ranking results achieved on all the datasets. The null hypothesis that the algorithms perform comparably is rejected withp value <0.00001.
Method
Ranking
D
ANCo
2.40
M
i
N
D
K
L
3.46
Hein
4.67
k
N
N
G
2
5.11
MLSVD
5.17
k
N
N
G
1
5.17
MLE
5.70
CD
6.63
BPCA
6.68Table 7
Hypothesis testing of significance between techniques. Bonferroni-Dunn’s procedure rejects those hypotheses that have ap value ≤0.0125.
M
i
N
D
K
L
Hein
k
N
N
G
1
k
N
N
G
2
MLE
CD
MLSVD
BPCA
D
ANCo
0.1567
0.0024
0.0003
0.0002
0.0002
0.0000
0.0000
0.0000
M
i
N
D
K
L
∗
∗
∗
0.0801
0.0303
0.0244
0.0055
0.0020
0.0000
0.0000
Hein
∗
∗
∗
∗
∗
∗
0.7528
0.6366
0.1564
0.1474
0.0034
0.0018
k
N
N
G
1
∗
∗
∗
∗
∗
∗
∗
∗
∗
0.8557
0.3443
0.2301
0.0164
0.0071
k
N
N
G
2
∗
∗
∗
∗
∗
∗
∗
∗
∗
∗
∗
∗
0.9314
0.3894
0.1113
0.0282
MLE
∗
∗
∗
∗
∗
∗
∗
∗
∗
∗
∗
∗
∗
∗
∗
0.3428
0.1876
0.0307
CD
∗
∗
∗
∗
∗
∗
∗
∗
∗
∗
∗
∗
∗
∗
∗
∗
∗
∗
0.7337
0.1961Finally, we performed the tests proposed in [27] to evaluate the robustness of MiNDKL, MLE, DANCo, and kNNG∗ with respect to the settings of their k parameter. Precisely, these tests employ synthetic datasets subsampled from the standard Gaussian pdf in R5 (id=5). As proposed in Section 4.2, we repeated the tests for datasets with cardinalities N∈{200,500,1000,2000} varying the parameter k in the range {5,…,100}.As shown in Figure3 many of the tested techniques are strongly influenced by the parameter settings; therefore, studying the variability of the algorithms’ behavior when changing their parameter settings is of utmost importance.Figure 3
Behavior of (a)MLE, (b) DANCo, (c) kNNG1, (d) kNNG2, and (e) MiNDKL applied to points drawn from a 5-dimensional standard normal distribution; in this test N∈{200,500,1000,2000} and k∈{5,…,100}.
(a)
(b)
(c)
(d)
(e)


## 5. Conclusions and Open Problems

This work presents the base theories of state-of-the-artid estimators and surveys the most relevant and recent among them, highlighting their strengths and their drawbacks.Unfortunately, performing an objective comparative evaluation among the surveyed methods is difficult because, to our knowledge, no benchmark framework exists in this research field; therefore, in Section4 we propose an evaluation approach that employs both real and synthetic datasets and suggests experiments to evaluate the estimators’ robustness with respect to their parameter settings. Note that, the benchmark is designed to evaluate the performance achieved by id estimators when both low and high id data must be processed; this consideration is due to the fact that, to our knowledge, only few methods [28, 76, 107, 126] have empirically investigated the problem of datasets drawn from manifolds nonlinearly embedded in higher dimensional spaces and characterized by a sufficiently high id (i.e., id⩾10). However, due to the continuous technological advances, high id datasets are becoming more and more common, and the construction of a theoretically well-formed and robust id estimator able to deal with high id data and limited amount of points remains one of the open research challenges in machine learning. Besides, id estimators should be developed by also considering datasets drawn through nonuniform smooth pdfs from manifolds M characterized by a nonconstant curvature; indeed, most of the algorithms are tested by only employing data drawn by means of uniform pdf.We further note that, though the aforementioned problems still need further investigations, most researches in this field are presently focusing on tasks that require to estimate theid as the first step. Examples are “multimanifold learning,” whose aim is to process datasets drawn from multiple manifolds, each characterized by different id, to identify the underlying structures (see [128] for an example); “nonlinear dimensionality reduction”; or “manifold reconstruction,” whose aim is to find the mapping that projects the data (embedded in a higher D-dimensional space) on the lower d^-dimensional subspace, d^ being the id estimated on the input dataset (as examples, see [9, 11, 129]).



---

*Source: 759567-2015-10-19.xml*