Intrinsic dimension estimation for discrete metrics

Data produced by experiments and observations are very often high-dimensional, with each data-point being defined by a sizeable number of features. To the please of modelers, real-world datasets seldom occupy this high-dimensional space uniformly, as strong regularities and constraints emerge. Such a property is what allows for low-dimensional descriptions of these high-dimensional data, ultimately making science possible.

In particular, data-points are often effectively contained in a manifold which can be described by a relatively small amount of coordinates. The number of such coordinates is called Intrinsic Dimension (ID). More formally, the ID is defined as the minimum number of variables needed to describe the data without significant information loss. Its knowledge is of paramount importance in unsupervised learning[1, 2, 3] and has found applications across disciplines. In solid state physics and statistical physics, the ID can be used as a proxy of an order parameter describing phase transitions[4, 5]; in molecular dynamics it can be used to quantify the complexity of a trajectory[6]; in deep learning theory the ID indicates how information is compressed throughout the various layers of a network [7, 8, 9]. During the last three decades much progress has been made in the development of sophisticated tools to estimate the ID[10, 11] and most estimators have been formulated (and are supposed to work) in spaces where distances can vary continuously. However, many datasets are characterised by discrete features and, consequently, discrete distances. For instance, categorical datasets like satisfaction questionnaires, clinical trials, unweighted networks, spin systems, protein, and DNA sequences fall into this category.

Two main methods are usually employed in these cases. The Box Counting (BC) estimator [12, 13, 14] — which is defined by measuring the scaling between the number of boxes needed to cover a dataset and the boxes size — provides good results for 2-3-dimensional datasets but is computationally demanding for higher-dimensional datasets. The second popular method is the Fractal Dimension (FD) estimator[12, 15, 16] and it is based on the assumption of a power law relationship $N\sim r^{d}$ for the number $N$ of neighbors within a sphere of radius $r$ from a given point, where $d$ is the fractal dimension of the data.This estimator has been successfully applied, on discrete datasets, to model the phenomena of dielectric breakdown [17] and Anderson Localization [18]. For non-fractal objects, both methods are reliable only in the limit of small boxes and small radii, since the manifold containing the data can be curved and the data points can be distributed non-uniformly[19]. However, in discrete spaces such a limit is not well defined due to the minimum distance induced by any discrete lattice, and this can lead to systematic errors [20, 21].

In this letter, we introduce an ID estimator explicitly formulated for spaces with discrete features. In discrete spaces, the ID can be thought of as the dimension of a (hyper)cubic lattice where the original data-points can be (locally) projected without a significant information loss. The key challenge in dealing with the discrete nature of the data lies in the proper definition of volumes on lattices. To this end, we introduce a novel method that makes use of Ehrhart’s theory of polytopes [22], which allows to enumerate the lattice points of a given region. By measuring a suitable statistics, depending on the number of data-points observed within a given (discrete) distance, one can infer the value of the dimension of the region, which we interpret as the ID of the dataset. The statistics we use is defined in such a way that density of points is required to be constant only locally and not in the whole dataset. Importantly, our estimator allows to explicitly select the scale at which the ID is computed.

Methods - We assume data points to be uniformly distributed on a generic domain, and that their density is $\rho$. In such domain, we consider a region $A$ with volume $V\!(A)$. Since we are assuming points to be independently generated, the probability of observing $n$ points in $A$ is given by the Poisson distribution [23]

so that ${\langle n\rangle}=\rho\,V\!(A)$. Consider now a data-point $i$ and two regions $A$ and $B$, one containing the other, and both containing the data-point: $i\in A\subset B$. Then the number of points $n$ and $k-n$ falling, respectively, in $A$ and $B\setminus A$ are Poisson distributed with rates $\lambda_{1}=\rho\,V\!(A)$ and $\lambda_{2}=\rho\,V\!(B\setminus A)$. The conditional probability of having $n$ points in $A$ given that there are $k$ points in $B$ is

with

Thus $n|k\sim\mathrm{Binomial}(n;k,p)$. As far as the density $\rho$ is constant within $A$ and $B$, $p$ is simply equal to the ratio of the volumes of the considered regions and, remarkably, density independent. This is a key property which, as we will show, allows using the estimator even when the density is approximately constant only locally, and varies, even substantially, across larger distance scales. One can then write a conditional probability of the observations $n_{i}$ (one for each data point), given the parameters $k_{i}$ and $p_{i}$, which can possibly be point-dependent:

Such formulation assumes all the observations to be statistically independent. Strictly speaking this is typically not true, since the regions $A$ and $B$ of different points can be overlapping. We will address this issue in Supplementary Information (SI), demonstrating that neglecting correlations does not induce significant errors.

The next step consists in defining the volumes in Eq. (3) according to the nature of the embedding manifold. We now assume our space to be a lattice where the $L^{1}$ metric is a natural choice. In this space the volume $V(A)$ is the number of lattice points contained in $A$. According to Ehrhart theory of polytopes [24], the number of lattice points within distance $t$ in dimension $d$ from a given point amounts to [25]

where ${}_{2}F_{1}(a,b,c,z)$ is the ordinary hypergeometric function. At a given $t$, the above expression is a polynomial in $d$ of order $t$. As a consequence, the ratio of volumes defining the value of $p$ in Eq. (3) becomes a ratio of two polynomials in $d$. Given a dataset, the choice of $t_{1}$ and $t_{2}$ fixes the values of $n_{i}$ and $k_{i}$ in the likelihood (4). Its maximization with respect to $d$ allows to infer the data-manifold’s ID, which is simply given the root of equation (see SI for more details on the derivation)

where the mean value over $n$ and $k$ is intended over all the points of the dataset. The root can be easily found with standard optimization libraries. This procedure defines an ID estimator that, for brevity, we will call I3D (Intrinsic Dimension estimator for Discrete Datasets).

Very importantly, the ID estimate is density independent as such factor cancels out (see Eq. (3)). The error on the estimator has a theoretical lower bound, given by the Cramer-Rao inequality, which has an explicit analytic expression. As an alternative, the ID can be estimated by a Bayesian approach as the mean value of its posterior distribution, and the error estimated via the posterior variance (details in SI).

The estimation of the ID depends on the choice of the volumes of the smaller and larger regions, which are parametrised by the “radii” $t_{1}$ and $t_{2}$. By varying $t_{2}$, the radius of the largest probe region, one can explore the behaviour of the ID at different scales. The proper range of $t_{2}$ is dataset dependent and should be chosen by plotting the value of the ID as a function of it, as we will illustrate in the following. If the dataset has a well-defined ID, one will observe an (approximate) plateau in this plot. This leaves the procedure with one free parameter: the ratio $r=t_{1}/t_{2}$ and its choice influences the statistical error. In continuous space the ratio between volumes in Eq. (3) is simply $p=r^{d}$ and the Cramer-Rao variance has a simple dependence on the parameter $r$. By minimising it with respect to $r$, one obtains that the optimal value for the ratio is $r_{opt}\sim 0.2032^{\frac{1}{d}}$ (see SI).

In order to check the goodness of the estimator, we test whether the number of points $n$ contained within the internal shells are actually distributed as a mixture of binomials, as our model assumes:

where $P(k)$ is the the empirical probability distribution of the $k$ found by fixing $t_{2}$. In the following we will compare the empirical cumulative distribution of $n$ to the cumulative distribution of $P(n)$.

Results: Uniform distribution - We tested the I3D estimator on artificial datasets, and compared it against the two aforementioned methods: the Box Counting (BC) and the Fractal Dimension (FD). The BC estimate of the ID is obtained by a linear fit between the logarithm of the number of occupied covering boxes and the logarithm of the boxes side. Seemingly, for the FD, the linear fit is computed among the logarithm of the average number of neighbours within a given radius and the logarithm of the radius. In both cases, the scale reported in the figures is given by the largest box or radius included in the fit. We started by analysing uniformly-distributed points in 2d and 6d square lattices. We adopted periodic boundary conditions in order to reduce boundary effects as much as possible. For the I3D estimator, in this and all following cases, we set $t_{1}/t_{2}=r=0.5$. Results are shown in Fig. 1.

While the BC and FD proved to be reliable in finding the fractal dimension of repeating, self-similar lattices [12, 17], they do not manage to assess the proper dimension of randomly distributed points, especially at small scales. The I3D estimator, instead, returns accurate values for the ID at all scales and, importantly, provides the correct estimate also on self-similar lattices (see SI). Remarkably the I3D estimator allows to select the scale explicitly by varying the radius $t_{2}$. In Fig. 1, lower panels, we also report a first example of model validation for I3D. The two cdfs (empirical and theoretical one, according to Eq. (7)) perfectly match, meaning that the ID estimation is reliable.

Gaussian distribution - Secondly, we tested the estimators on Gaussian distributed points in 5 dimensions, analysing a case in which the data are uncorrelated and a case in which a correlation is induced by a non-diagonal covariance matrix. In both cases, we set diagonal elements of the covariance matrix to $\sigma=5$ (implying an effective standard deviation of the distribution of $\sigma_{eff}=\sqrt{d}\sigma$), while off-diagonal terms –for correlated data– were uniformly extracted in the interval (0,2). The values were chosen in order to keep the dimension of the dataset under control, as correlations of the same order of the diagonal would reduce the dimensionality of the dataset. The points were projected on a lattice by taking the nearest integer in each coordinate.

As one can observe in panel A of Fig.2, I3D is accurate as far as it explores a neighborhood where the density does not vary too much (namely, as far as $t_{2}/\sigma_{eff}\lesssim 1$). Correspondingly, empirical and model cdfs in panel B are superimposed. Beyond such distance, neighborhoods are characterised by non-constant density; consequently, estimates gets less precise and, accordingly, the two cdfs show inconsistencies (panel C: $t_{2}/\sigma_{eff}\sim 1.5$). On the other hand, the BC and FD estimations are far from desired values at any scales, for both correlated and uncorrelated cases.

Spin dataset - As a third test, we created synthetic Ising-like spin systems with a tunable ID, which is given by the number of independent parameters used to generate the dataset. The 1d ensemble is obtained by generating a set of points belonging to a line embedded in $\mathbb{R}^{D}$ with the process $\bm{\varphi}_{i}=\bm{\varphi}_{0}+\bm{\alpha}\epsilon(i)$. Here, $\bm{\alpha}$ is a fixed random vector of unitary norm with uniformly distributed components and $\bm{\varphi}_{0}=-0.5$ is the y-intercept that, for simplicity, is equal for all the components; $\epsilon_{i}$ are gaussian-distributed: $\epsilon\sim\mathcal{N}(0,10)$ and independently drawn for each sample $i$. We then proceed to the discretization by extracting the $\bm{z}_{i}=\text{sign}(\bm{\varphi}_{i})$, an ensemble of $N$ states of $D$ discrete spins. The pipeline is summarized in Fig.3. The role of $\bm{\varphi}_{0}$ is to introduce an offset in order to enhance the number of the reachable discrete states. In fact, for $\bm{\varphi}_{0}=0$, we would obtain only two different states, given by $\bm{z}=\text{sign}(\bm{\alpha}\epsilon)=\pm\text{sign}(\bm{\alpha})$, since the spins would change sign synchronously. An offset $\neq 0$ allows the angles $\bm{\varphi}_{i}$ and the spins $\bm{z}_{i}$ to shift sign in an asynchronous way. The extension to higher dimensions is straightforward and consists in generating the initial points as $\bm{\varphi}_{i}=\bm{\varphi}_{0}+\sum_{j=1}^{id}\bm{\alpha}_{j}\epsilon_{j}(i)$, with $\bm{\alpha}_{j}\cdot\bm{\alpha}_{k}\sim\delta_{jk}$. Due to the nature of data domain (a $D$-dimensional hypercube with side 1), the BC cannot be applied, as boxes with side larger than 1 would include the whole data set. FD and I3D estimates for the 1d system are very close. This is not surprising as both continuous and discrete volumes (and, consequently, the neighbors) scale linearly with the radius. In the 2d case, I3D clearly outperforms the other methods, although even the best estimate remains slightly lower than the true value. This effect, due to non-uniform density, is relatively small and indeed the empirical and theoretical cdfs are rather consistent (panels B and C). Such an effect becomes more important as the dimension rises (see SI for examples in $d=3$ and $d=4$).

16S Genomics strands - Lastly, we present the application of our methodology to a real-world dataset in the field of genomics. The dataset consists of DNA sequences of $\sim\!100\!-\!300$ nucleotides. We selected a dataset downloaded from the Qiita server (https://qiita.ucsd.edu/study/description/13596)[26]. In such study, they sequenced the v4 region of the 16S ribosomal RNA of the microbiome associated with sponges and algal water blooms. This small-subunit of rRNA genes is widely used to study the composition of microbial communities [27, 28, 29, 30]. Hamming distance and the binary mapping A:11, T:00, C:10 and G:01 were used to compute sequences’ distance. The canonical letter representation leads to almost identical results (see SI). To avoid dealing with isolated sequences, we kept only sequences having at least 10 neighbors within a distance of 10. Sequences come with their associated multiplicity, related to the number of times the same read has been found in the samples. We ignore such degeneracy and compute an ID which describes just the distribution of the points regardless of their abundance.

To begin with, we estimated the ID on a subset of sequences that are similar to each other. In order to find such sets, we perform a k-means clustering and calculate the ID separately for each of them. Panel A in Fig. 4 shows the ID at small to medium scale for one of such clusters. The empirical and reconstructed cdfs, performed at $t_{2}=20$ (see inset), are fairly compatible. Panel B shows the average and the standard deviation of the ID of all clusters (weighted according to the respective populations). One can appreciate that the ID is always between 1 and 3 in a wide range of distances, showing a plateau around 2 for $15<t_{2}<40$.

Such a low value for the ID is an interesting and unexpected feature, as it suggests that, despite the high-dimensionality of sequences’ space, evolution effectively operates in a low-dimensional space. Qualitatively, an ID $\sim$ 2 on a scale of $\sim$ 20 means that if one considers all the sequences differing by approximately 20 mutations from a given sequence, these mutations cannot be regarded as independent one from each other, but are correlated in such a way that approximately 18 degrees of freedom are effectively forbidden. The “direction” of these correlated mutations can be, at least approximately, measured by performing PCA in the space of sequences with the binary mapping. The first two dominant eigenvectors, shown in panel C, were estimated using all the sequences within a distance of 20 (top) and 30 (bottom) from the center of the cluster of Panel A. Remarkably, the eigenvectors do not change significantly on this distance range, indicating that, consistently with the low value of the ID, the data manifold on this scale can be approximately described by a two-dimensional plane. In order to provide an interpretation of the vectors defining this plane, we repeated this same analysis on the previously mentioned spin model. In this case, if the generative model is defined by two vectors $\bm{\alpha}_{1}$ and $\bm{\alpha}_{2}$, the first two dominant eigenvectors of a PCA performed on $\sim$1000 points are contained in the span of the two generating vectors, with a residual of 0.04 (see SI for details). The components of a vector $\bm{\alpha}$ can then be qualitatively interpreted as proportional to the mutation probabilities of the associated nucleotide for a collective mutation process. In the genomics dataset this reasoning can applied only locally: the direction of correlated mutation is significantly different in different clusters, indicating that the data manifold is highly curved.

Conclusions - We presented an ID estimator formulated to analyze discrete datasets. Our method relies on few mathematical hypotheses and is asymptotically correct if the density is constant within the probe radius $t_{2}$. In order to prove the estimator’s effectiveness, we tested the algorithm against three different artificial datasets and compared it to the well known Box Counting and Fractal Dimension estimators. While the last two performed poorly, the new one achieved good results in all cases, providing reliable ID estimations corroborated by the comparison of empirical and model cumulative distribution functions for one of the observables. We finally applied the estimator on a genomics dataset, finding an unexpectedly low ID which hints at strong constraints in the sequences’ space, and then exploited such information to give a qualitative interpretation of such ID. The newly developed method paves the way to push the investigation even further, towards the extension to discrete metrics of distance-based algorithms and routines that are, nowadays, consolidated in the continuum, such as density estimation methods or clustering algorithms.

The code implementing the algorithm is available in open source within the DADApy [31] software.

Intrinsic dimension estimation for discrete metrics
— Supplementary Information —