Data segmentation based on the local intrinsic dimension

Michele Allegra; Elena Facco; Francesco Denti; Alessandro Laio,; Antonietta Mira

arXiv:1902.10459·stat.ML·July 14, 2020

Data segmentation based on the local intrinsic dimension

Michele Allegra, Elena Facco, Francesco Denti, Alessandro Laio,, Antonietta Mira

PDF

1 Repo

TL;DR

This paper introduces a robust, efficient method to segment high-dimensional data based on local intrinsic dimension variations, revealing heterogeneous regions with distinct structural properties across various real-world datasets.

Contribution

The authors develop a novel approach to identify and segment regions with different local intrinsic dimensions in large high-dimensional datasets.

Findings

01

Many real-world datasets contain regions with heterogeneous intrinsic dimensions.

02

Local ID segmentation can distinguish different structural states in data.

03

The method is computationally efficient and scalable to large datasets.

Abstract

One of the founding paradigms of machine learning is that a small number of variables is often sufficient to describe high-dimensional data. The minimum number of variables required is called the intrinsic dimension (ID) of the data. Contrary to common intuition, there are cases where the ID varies within the same data set. This fact has been highlighted in technical discussions, but seldom exploited to analyze large data sets and obtain insight into their structure. Here we develop a robust approach to discriminate regions with different local IDs and segment the points accordingly. Our approach is computationally efficient and can be proficiently used even on large data sets. We find that many real-world data sets contain regions with widely heterogeneous dimensions. These regions host points differing in core properties: folded vs unfolded configurations in a protein molecular…

Tables3

Table 1. Table 1: Performances of Hidalgo, restricting analysis to points with non-overlapping first and second neighbor For each data set we show the number of points in the original data set ( N 0 subscript 𝑁 0 N_{0} ), the number of points upon restriction N 𝑁 N , the estimated value of the d k subscript 𝑑 𝑘 d_{k} , and NMI betzeen the assignment and the ground truth. Qll results were obtained with ξ = 0.8 , q = 3 formulae-sequence 𝜉 0.8 𝑞 3 \xi=0.8,q=3 , repeating the sampling M = 10 𝑀 10 M=10 times and considering the sampling with the highest average log-posterior.

Dataset		Hidalgo d	NMI	Ref. [19] d	NMI	Ref. [10] d	NMI
2 Gaussians	$d_{1} = 5$ $d_{2} = 4$	$d_{1} = 4.7$ $d_{2} = 4.2$	0.91	$d_{1} = 5$ $d_{2} = 4$	0.96	$⟨ d ⟩ = 4.0$ $⟨ d ⟩ = 4.7$	0.44
2 Gaussians	$d_{1} = 6$ $d_{2} = 4$	$d_{1} = 5.8$ $d_{2} = 4.2$	0.93	$d_{1} = 6$ $d_{2} = 4$	0.94	$⟨ d ⟩ = 4.0$ $⟨ d ⟩ = 5.1$	0.95
2 Gaussians	$d_{1} = 7$ $d_{2} = 4$	$d_{1} = 6.9$ $d_{2} = 4.1$	0.94	$d_{1} = 6$ $d_{2} = 4$	0.97	$⟨ d ⟩ = 4.0$ $⟨ d ⟩ = 6.0$	0.96
2 Gaussians	$d_{1} = 8$ $d_{2} = 4$	$d_{1} = 7.6$ $d_{2} = 4.2$	0.92	$d_{1} = 7$ $d_{2} = 4$	0.94	$⟨ d ⟩ = 4.0$ $⟨ d ⟩ = 6.7$	0.93
2 Gaussians	$d_{1} = 9$ $d_{2} = 4$	$d_{1} = 9.0$ $d_{2} = 4.1$	0.92	$d_{1} = 7$ $d_{2} = 4$	0.95	$⟨ d ⟩ = 4.0$ $⟨ d ⟩ = 7.2$	0.92
5 Gaussians (linear)	$d_{1} = 1$ $d_{2} = 2$ $d_{3} = 4$ $d_{4} = 5$ $d_{5} = 9$	$d_{1} = 0.9$ $d_{2} = 1.9$ $d_{3} = 4.2$ $d_{4} = 4.6$ $d_{5} = 8.9$	0.89	$d_{1} = 7, d_{4} = 9, d_{5} = 10$ $d_{2} = 10$ $d_{3} = 9$ $d_{4} = 9$ $d_{5} = 10$	0.73	$⟨ d ⟩ = 1.0$ $⟨ d ⟩ = 2.0$ $⟨ d ⟩ = 3.5$ $⟨ d ⟩ = 4.3$ $⟨ d ⟩ = 7.3$	0.81
5 Gaussians (curved)	$d_{1} = 1$ $d_{2} = 2$ $d_{3} = 4$ $d_{4} = 5$ $d_{5} = 9$	$d_{1} = 0.9$ $d_{2} = 2.2$ $d_{2} = 2.2, d_{3} = 3.8$ $d_{4} = 5.2$ $d_{5} = 8.5$	0.84	$d_{1} = 6, d_{3} = 7$ $d_{2} = 10$ $d_{3} = 7, d_{4} = 8$ $d_{5} = 10$ $d_{5} = 10$	0.61	$⟨ d ⟩ = 1.0$ $⟨ d ⟩ = 2.0$ $⟨ d ⟩ = 2.0$ $⟨ d ⟩ = 4.6$ $⟨ d ⟩ = 7.2$	0.77
Protein Folding	$Q \leq 0.8$ $Q > 0.8$	$d_{1} = 11.8$ $d_{2} = 12.9$ $d_{3} = 13.2$ $d_{4} = 22.3$	0.60	$d_{1} = 11$	0	$⟨ d ⟩ = 11.4$ $⟨ d ⟩ = 13.7$	0.22

Table 2. Table S1: Independent points. We repeated the analysius of artifical data sets, using only indpedendent points with non-overlapping first and second neighbor For each scenario we show the number of points in the original data set ( N 0 subscript 𝑁 0 N_{0} ), the number of points upon restriction N 𝑁 N , the estimated value of the d k subscript 𝑑 𝑘 d_{k} , and NMI between the assignment and the ground truth. All results were obtained with ξ = 0.8 , q = 3 formulae-sequence 𝜉 0.8 𝑞 3 \xi=0.8,q=3 , repeating the sampling M = 10 𝑀 10 M=10 times, and considering the sampling with the highest average log-posterior.

Dataset	$𝐍_{𝟎}$	$𝐍$	$𝐝_{𝐢}$	NMI
2 Gaussians, $d_{1} = 5, d_{2} = 4$	2000	357	$d_{1} = 5.5, d_{2} = 3.9$	0.95
2 Gaussians, $d_{1} = 6, d_{2} = 4$	2000	360	$d_{1} = 6.2, d_{2} = 4.1$	0.93
2 Gaussians, $d_{1} = 7, d_{2} = 4$	2000	355	$d_{1} = 6.5, d_{2} = 4.3$	0.93
2 Gaussians, $d_{1} = 8, d_{2} = 4$	2000	365	$d_{1} = 9.1, d_{2} = 3.9$	0.96
5 Gaussians (linear), $d_{1} = 1, d_{2} = 2, d_{3} = 4, d_{4} = 5, d_{5} = 9$	5000	895	$d_{1} = 0.9, d_{2} = 1.7, d_{3} = 3.2, d_{4} = 4.4, d_{5} = 9.4$	0.88
5 Gaussians (curved), $d_{1} = 1, d_{2} = 2, d_{3} = 4, d_{4} = 5, d_{5} = 9$	5000	894	$d_{1} = 0.9, d_{2} = 1.7, d_{3} = 3.2, d_{4} = 4.4, d_{5} = 9.4$	0.61

Table 3. Table S2: Compustat variables used

#	Variable	#	Variable	#	Variable
1	Acquisitions	12	Liabilities - Total	23	Interest and Related Expense - Total
2	Assets - Total	13	Net Income (Loss)	24	Goodwill
3	Capital Expenditures	14	Operating Income Before Depreciation	25	Intangible Assets - Total
4	Cash	15	Property Plant and Equipment - Total (Net)	26	Pretax Income
5	Common Shares Outstanding	16	Purchase of Common and Preferred Stock	27	Pretax Income - Foreign
6	Common/Ordinary Shareholders	17	Sales/Turnover (Net)	28	Investment and Advances - Equity
7	Debt in Current Liabilities - Total	18	Stockholders Equity - Parent	29	Investment and Advances - Other
8	Long-Term Debt - Total	19	Income Taxes Paid	30	Increase in Investments
9	Cash Dividends (Cash Flow)	20	Research and Development Expense	31	Sale of Investments
10	Earnings Before Interest and Taxes	21	Price Close - Annual - Fiscal
11	Employees	22	Preferred/Preference Stock (Capital) - Total

Equations52

P (μ ∣ d) = d^{N} i = 1 \prod N μ_{i}^{- (d + 1)} = d^{N} e^{- (d + 1) V},

P (μ ∣ d) = d^{N} i = 1 \prod N μ_{i}^{- (d + 1)} = d^{N} e^{- (d + 1) V},

P (μ_{i} ∣ d, p) ≐ k = 1 \sum K p_{k} d_{k} μ_{i}^{- (d_{k} + 1)} .

P (μ_{i} ∣ d, p) ≐ k = 1 \sum K p_{k} d_{k} μ_{i}^{- (d_{k} + 1)} .

P_{p os t} (z, d, p ∣ μ) \propto P (μ ∣ z, d) P_{p r} (z ∣ p) P_{p r} (d) P_{p r} (p) .

P_{p os t} (z, d, p ∣ μ) \propto P (μ ∣ z, d) P_{p r} (z ∣ p) P_{p r} (d) P_{p r} (p) .

P (N^{(q)} ∣ z) = i = 1 \prod N \frac{ξ ^{n_{i}^{in}} ( 1 - ξ ) ^{q - n_{i}^{in}}}{Z ( ξ , N _{z_{i}} )},

P (N^{(q)} ∣ z) = i = 1 \prod N \frac{ξ ^{n_{i}^{in}} ( 1 - ξ ) ^{q - n_{i}^{in}}}{Z ( ξ , N _{z_{i}} )},

P_{p os t} (z, d, p ∣ μ, N^{(q)}) \propto P (μ ∣ z, d) P (N^{(q)} ∣ z) P_{p r} (z ∣ p) P_{p r} (d) P_{p r} (p) .

P_{p os t} (z, d, p ∣ μ, N^{(q)}) \propto P (μ ∣ z, d) P (N^{(q)} ∣ z) P_{p r} (z ∣ p) P_{p r} (d) P_{p r} (p) .

\mathcal{N}_{ij}^{(q)}=\begin{cases}1&\mbox{if }j\neq i\ \mbox{\mbox{is among the first $q$ neighbors of }}i\\ 0&\mbox{otherwise, including }i=j\end{cases}

\mathcal{N}_{ij}^{(q)}=\begin{cases}1&\mbox{if }j\neq i\ \mbox{\mbox{is among the first $q$ neighbors of }}i\\ 0&\mbox{otherwise, including }i=j\end{cases}

P_{p os t} (z_{i} = z_{j} ∣ N_{ij}^{(q)} = 1, μ, p) > P_{p os t} (z_{i} = z_{j} ∣ N_{ij}^{(q)} = 0, μ, p) .

P_{p os t} (z_{i} = z_{j} ∣ N_{ij}^{(q)} = 1, μ, p) > P_{p os t} (z_{i} = z_{j} ∣ N_{ij}^{(q)} = 0, μ, p) .

L (N_{i}^{(q)} ∣ z) = \frac{ξ ^{n_{i}^{in} (z)} ( 1 - ξ ) ^{q - n_{i}^{in} (z)}}{Z ( ξ , N _{z_{i}} )},

L (N_{i}^{(q)} ∣ z) = \frac{ξ ^{n_{i}^{in} (z)} ( 1 - ξ ) ^{q - n_{i}^{in} (z)}}{Z ( ξ , N _{z_{i}} )},

n_{i}^{in} (z) = j \sum N_{ij}^{(q)} I_{z_{j} = z_{i}}

n_{i}^{in} (z) = j \sum N_{ij}^{(q)} I_{z_{j} = z_{i}}

q - n_{i}^{in} (z) = j \sum N_{ij}^{(q)} I_{z_{j} \neq = z_{i}}

q - n_{i}^{in} (z) = j \sum N_{ij}^{(q)} I_{z_{j} \neq = z_{i}}

Z (ξ, N_{z_{i}}) = {N_{i}^{(q)}} \sum ξ^{n_{i}^{in} (z)} (1 - ξ)^{q - n_{i}^{in} (z)} .

Z (ξ, N_{z_{i}}) = {N_{i}^{(q)}} \sum ξ^{n_{i}^{in} (z)} (1 - ξ)^{q - n_{i}^{in} (z)} .

Z (ξ, N_{z_{i}}) = (1 - ξ)^{q} (q N - N _{z_{i}})_{2} F_{1} (- q, 1 - N_{z_{i}}, N - N_{z_{i}} - q, \frac{ξ}{1 - ξ}),

Z (ξ, N_{z_{i}}) = (1 - ξ)^{q} (q N - N _{z_{i}})_{2} F_{1} (- q, 1 - N_{z_{i}}, N - N_{z_{i}} - q, \frac{ξ}{1 - ξ}),

L (N^{(q)} ∣ z, ξ) = i \prod L (N_{i}^{(q)} ∣ z, ξ) = k \prod \frac{ξ ^{n_{k}^{in}} ( 1 - ξ ) ^{q N_{k} - n_{k}^{in}}}{Z ( ξ , N _{k} ) ^{N_{k}}}

L (N^{(q)} ∣ z, ξ) = i \prod L (N_{i}^{(q)} ∣ z, ξ) = k \prod \frac{ξ ^{n_{k}^{in}} ( 1 - ξ ) ^{q N_{k} - n_{k}^{in}}}{Z ( ξ , N _{k} ) ^{N_{k}}}

n_{k}^{in} = ij \sum N_{ij}^{(q)} I_{z_{i} = k} I_{z_{j} = k}

n_{k}^{in} = ij \sum N_{ij}^{(q)} I_{z_{i} = k} I_{z_{j} = k}

n_{k}^{o u t} = ij \sum N_{ij}^{(q)} I_{z_{i} = k} (1 - I_{z_{j} = k}) = q N_{k} - n_{k}^{in}

n_{k}^{o u t} = ij \sum N_{ij}^{(q)} I_{z_{i} = k} (1 - I_{z_{j} = k}) = q N_{k} - n_{k}^{in}

\frac{P _{p os t} ( z _{i} = z _{j} ∣ N _{ij}^{(q)} = 1 , μ , p )}{P _{p os t} ( z _{i} = z _{j} ∣ μ , p )} = \frac{ξ}{1 - ξ} > 1/2.

\frac{P _{p os t} ( z _{i} = z _{j} ∣ N _{ij}^{(q)} = 1 , μ , p )}{P _{p os t} ( z _{i} = z _{j} ∣ μ , p )} = \frac{ξ}{1 - ξ} > 1/2.

L (N_{i}^{(q)} ∣ z) = (q N - 1)^{- 1}, \forall i .

L (N_{i}^{(q)} ∣ z) = (q N - 1)^{- 1}, \forall i .

(n _{b} N _{b}) (q - n _{b} N _{w})

(n _{b} N _{b}) (q - n _{b} N _{w})

\frac{ξ ^{n_{b}} ( 1 - ξ ) ^{q - n_{b}}}{Z}

\frac{ξ ^{n_{b}} ( 1 - ξ ) ^{q - n_{b}}}{Z}

_{2} F_{1} (- m, b, c, z) = n = 0 \sum m (-)^{n} (n m) \frac{( b ) _{n}}{( c ) _{n}} z^{n}

_{2} F_{1} (- m, b, c, z) = n = 0 \sum m (-)^{n} (n m) \frac{( b ) _{n}}{( c ) _{n}} z^{n}

Z = (1 - ξ)^{q} (q N _{w})_{2} F_{1} (- q, - N_{b}, N_{w} - q, \frac{ξ}{1 - ξ}) .

Z = (1 - ξ)^{q} (q N _{w})_{2} F_{1} (- q, - N_{b}, N_{w} - q, \frac{ξ}{1 - ξ}) .

L (N_{i}^{(q)} ∣ z, ξ) = \frac{ξ ^{n_{i}^{in} (z)} ( 1 - ξ ) ^{q - n_{i}^{in} (z)}}{( 1 - ξ ) ^{q} ( q N - N _{z_{i}} ) _{2} F _{1} ( - q , 1 - N _{z_{i}} , N - N _{z_{i}} - q , \frac{ξ}{1 - ξ} )} .

L (N_{i}^{(q)} ∣ z, ξ) = \frac{ξ ^{n_{i}^{in} (z)} ( 1 - ξ ) ^{q - n_{i}^{in} (z)}}{( 1 - ξ ) ^{q} ( q N - N _{z_{i}} ) _{2} F _{1} ( - q , 1 - N _{z_{i}} , N - N _{z_{i}} - q , \frac{ξ}{1 - ξ} )} .

S_{t} = τ = t_{0} + 1 \sum t (f_{τ} - μ)

S_{t} = τ = t_{0} + 1 \sum t (f_{τ} - μ)

e_{t} = ⎩ ⎨ ⎧ 1 \mbox i f S_{t} > S_{t - 1} \mbox an d S_{t} < S_{t - 1} 1 \mbox i f S_{t} < S_{t - 1} \mbox an d S_{t} > S_{t - 1} 0 \mbox o t h er w i se

e_{t} = ⎩ ⎨ ⎧ 1 \mbox i f S_{t} > S_{t - 1} \mbox an d S_{t} < S_{t - 1} 1 \mbox i f S_{t} < S_{t - 1} \mbox an d S_{t} > S_{t - 1} 0 \mbox o t h er w i se

E_{t} = τ = t_{0} + 1 \sum t - 1 e_{τ}

E_{t} = τ = t_{0} + 1 \sum t - 1 e_{τ}

1/2 - Z_{α /2} / t - t_{0} - 1 \leq E_{t} \leq 1/2 + Z_{α /2} / t - t_{0} - 1

1/2 - Z_{α /2} / t - t_{0} - 1 \leq E_{t} \leq 1/2 + Z_{α /2} / t - t_{0} - 1

Peer Reviews

No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.

Code & Models

Repositories

micheleallegra/Hidalgo
noneOfficial

Videos

No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.

Full text

Data segmentation based on the local intrinsic dimension

Michele Allegra

Institut de Neurosciences de la Timone UMR 7289, Aix Marseille Université, CNRS, 13385 Marseille, France

Scuola Internazionale Superiore di Studi Avanzati, Trieste, Italy

Elena Facco

Scuola Internazionale Superiore di Studi Avanzati, Trieste, Italy

Francesco Denti

University of California, Irvine, California, U.S.A.

Alessandro Laio

Scuola Internazionale Superiore di Studi Avanzati, Trieste, Italy

International Centre for Theoretical Physics, Trieste, Italy

[email protected]

Antonietta Mira

Università della Svizzera italiana, Lugano, Switzwerland

Università dell’Insubria, Como, Italy

[email protected]

Abstract

One of the founding paradigms of machine learning is that a small number of variables is often sufficient to describe high-dimensional data. The minimum number of variables required is called the intrinsic dimension (ID) of the data. Contrary to common intuition, there are cases where the ID varies within the same data set. This fact has been highlighted in technical discussions, but seldom exploited to analyze large data sets and obtain insight into their structure. Here we develop a robust approach to discriminate regions with different local IDs and segment the points accordingly. Our approach is computationally efficient and can be proficiently used even on large data sets. We find that many real-world data sets contain regions with widely heterogeneous dimensions. These regions host points differing in core properties: folded vs unfolded configurations in a protein molecular dynamics trajectory, active vs non-active regions in brain imaging data, and firms with different financial risk in company balance sheets. A simple topological feature, the local ID, is thus sufficient to achieve an unsupervised segmentation of high-dimensional data, complementary to the one given by clustering algorithms.

Introduction

From string theory to science fiction, the idea that we might be “glued” onto a low-dimensional surface embedded in a space of large dimensionality has tickled the speculations of scientists and writers alike. However, when it comes to multidimensional data such a situation is quite common rather than wild speculation: data often concentrate on hypersurfaces of low intrinsic dimension (ID). Estimating the ID of a dataset is a routine task in machine learning: it yields important information on the global structure of a dataset, and is a necessary preliminary step in several analysis pipelines.

In most approaches for dimensionality reduction and manifold learning, the ID is assumed to be constant in the dataset. This assumption is implicit in projection-based estimators, such as Principal Component Analysis (PCA) and its variants [1], Locally Linear Embedding [2], and Isomap [3]; and it also underlies geometric ID estimators [4, 5, 6], which infer the ID from the distribution of between-point distances. The hypothesis of a constant ID complies with simplicity and intuition but is not necessarily valid. In fact, many authors have considered the possibility of ID variations within a dataset [7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22], often proposing to classify the data according to this feature. However, the dominant opinion in the community is still that a variable ID is a peculiarity, or a technical detail, rather than a common feature to take into account before performing a data analysis. This perception is in part due to the lack of sufficiently general methods to track ID variations. Many of the methods developed in this field make assumptions that limit their applicability to specific classes of data. Refs. [7, 11, 12, 13, 14] use local ID estimators that implicitly or explicitly assume a uniform density. Refs. [8, 9, 10] jointly estimate the density and the ID from the scaling of neighbor distances, by approaches which work well only if the density varies slowly and is approximately constant in a large neighborhood of each point. Ref. [15] requires a priori knowledge on the number of manifolds and their IDs. Refs. [16, 17, 18, 20] all require that the manifolds on which the data lay are hyperplanes, or topologically isomorphic to hyperplanes. These assumptions (locally constant density, and linearity in a suitable set of coordinates) are often violated in the case of real-world data. Moreover, many of the above approaches [7, 11, 12, 13, 15, 16, 17, 18, 19, 20] work explicitly with the coordinates of the data, while in many applications only distances between pairs of data points are available. To our knowledge, only refs. [21, 22] do not make any assumption about the density, as they derive a parametric form of the distance distribution using extreme-value theory, which in principle is valid independently of the form of the underlying density. However, they assume that the low tail of the distance distribution is well approximated by its asymptotic form, an equally non-trivial assumption.

In this work, we propose a method to segment data based on the local ID which aims at overcoming the aforementioned limitations. Building on TWO-NN [23], a recently proposed ID estimator which is insensitive to density variations and uses only the distances between points, we develop a Bayesian framework that allows identifying, by Gibbs sampling, the regions in the data landscape in which the ID can be considered constant. Our approach works even if the data are embedded on highly curved and twisted manifolds, topologically complex and not isomorphic to hyperplanes, and if they are harvested from a non-uniform probability density. Moreover, it is specifically designed to use only the distance between data points, and not their coordinates. These features, as we will show, make our approach computationally efficient and more robust than other methods. Applying our approach to data of various origins, we show that ID variations between different regions are common. These variations often reveal fundamental properties of the data: for example, unfolded states in a molecular dynamics trajectory of a protein lie on a manifold of a lower dimension than the one hosting the folded states. Identifying regions of different dimensionality in a dataset can thus be a way to perform an unsupervised segmentation of the data. At odds with common approaches, we do not group the data according to their density but perform segmentation based on a geometric property defined on a local scale: the number of independent directions along which neighboring data are spread.

Methods

A Bayesian approach for discriminating manifolds with different IDs

The starting point of our approach is the recently proposed TWO-NN estimator [23], which infers the IDs from the statistics of the distances of the first two neighbors of each point. Let the data $(x_{1},x_{2},...,x_{N})$ , with $N$ indicating the number of points, be i) independent and identically distributed samples from a probability density $\rho(x)$ with support on a manifold with unknown dimension $d$ , such that ii) for each point $x_{i}$ , $\rho(x_{i})$ is approximately constant in the spherical neighborhood with center $i$ and radius given by the distance between $i$ and its second neighbor. Assumption i) is shared all methods that infer the ID from the distribution of distances (e.g., Refs. [5, 6, 9, 10, 11, 12, 14]. Assumption ii) is the key assumption of TWO-NN: while of course it may not be perfectly satisfied, other methods used to infer the ID from distances require uniformity of the density on a larger scale including the first $k$ neighbors, usually with $k\gg 2$ . Let $r_{i1}$ and $r_{i2}$ be the distances of the first and second neighbor of $x_{i}$ , then $\mu_{i}\doteq\frac{r_{i2}}{r_{i1}}$ follows the Pareto distribution $f(\mu_{i}|d)=d\mu_{i}^{-(d+1)}$ . This readily allows the estimation of $d$ from ${\mu_{i}},i=1,2,...,N$ . We can write the global likelihood of $\boldsymbol{\mu}\equiv\{\mu_{i}\}$ as

[TABLE]

where $V\doteq\sum_{i=1}^{N}\log(\mu_{i})$ . From Eq. (1), and upon specifying a suitable prior on $d$ , a Bayesian estimate of $d$ is immediately obtained. This estimate is particularly robust to density variations: in [23] it was shown that restricting the uniformity requirement only to the region containing the first two neighbors allows properly estimating the ID also in cases where $\rho$ is strongly varying.

TWO-NN can be extended to yield a heterogeneous-dimensional model with an arbitrarily high number of components. Let $\boldsymbol{x}=\{x_{i}\}$ be i.i.d samples from a density $\rho(x)$ with support on the union of $K$ manifolds with varying dimensions. This multi-manifold framework is common with many previous works investigating heterogeneous dimension in a dataset [8, 9, 15, 10, 14, 11, 24, 25, 16, 19, 20, 18]. Formally, let $\rho(x)=\sum_{k=1}^{K}p_{k}\rho_{k}(x)$ where each $\rho_{k}(x)$ has support on a manifold of dimension $d_{k}$ and $\boldsymbol{p}\doteq(p_{1},p_{2},...,p_{K})$ are the a priori probabilities that a point belongs to the manifolds $1,\dots,K$ . We shall first assume that $K$ is known, and later show how it can be estimated from the data. The distribution of the $\mu_{i}$ is simply a mixture of Pareto distributions:

[TABLE]

Here, like in ref. [23], we assume that the variables $\mu_{i}$ are independent. This condition could be rigorously imposed by restricting the product in eq. 1 to data points with no common first and second neighbors. In fact, one does not find significant differences between the empirical distribution of $\boldsymbol{\mu}$ obtained with all points and the one obtained by restricting the sample to the independent points (i.e., the empirical distribution of $\mu$ , obtained considering entire sample, is well fitted by the theoretical Pareto law). However, restricting the estimation to independent samples generally worsens the estimates, as one is left with a small fraction (typically, $\approx 15\%$ ) of the original points. Therefore, for the purpose of parametric inference, we decide to include all points in the estimation. In SI (Table S1) we show that restricting the estimation to independent points yields results that are qualitatively in good agreement with those obtained using all points, but with generally worse dimension estimates. Following the customary approach [26] we introduce latent variables $\boldsymbol{z}\doteq(z_{1},z_{2},...,z_{K})$ where $z_{i}=k$ indicates that point $i$ belongs to manifold $k$ . We have $P(\mu_{i}|\boldsymbol{d},\boldsymbol{p},\boldsymbol{z})=P(\mu_{i}|z_{i},\boldsymbol{d})P_{pr}(z_{i}|\boldsymbol{p})$ with with $P(\mu_{i}|z_{i},\boldsymbol{d})=d_{z_{i}}\mu_{i}^{-(d_{z_{i}}+1)}$ , $P_{pr}(z_{i}|\boldsymbol{p})=p_{z_{i}}$ . This yields the posterior

[TABLE]

We use independent Gamma priors on $\boldsymbol{d}$ , $d_{k}\sim Gamma(a_{k},b_{k})$ and a joint Dirichlet prior on $\boldsymbol{p}\sim Dir(c_{1},\dots,c_{K})$ . We fix $a_{k},b_{k},c_{k}=1\ \forall k$ , corresponding to a maximally non-informative prior on the $\boldsymbol{p}$ and an expectation of $\mathbb{E}\left[d_{k}\right]=1$ for all $k$ . If one has a different prior expectation on $\boldsymbol{d}$ and $\boldsymbol{p}$ , other choices of prior may be more convenient.

The posterior (3) does not have an analytically simple form, but it can be sampled by standard Gibbs sampling [27], allowing for the joint estimation of $\boldsymbol{d},\boldsymbol{p},\boldsymbol{z}$ . However, model (3) has a serious limitation: Pareto distributions with (even largely) different values of $d$ overlap to a great extent. Therefore, the method can not be expected to correctly estimate the $z_{i}$ : a given value $\mu_{i}$ may be compatible with several manifold memberships. This issue can be addressed by correcting an unrealistic feature of model (3), namely, the independence of the $z_{i}$ . We assume that the neighborhood of a point is more likely to contain points from the same manifold than from different manifolds. This requirement can be enforced with an additional term in the likelihood that penalizes local inhomogeneity (see also Ref. [14]). Consider the $q$ -neighbor matrix $\boldsymbol{\mathcal{N}}^{(q)}$ with nonzero entries $\mathcal{N}_{ij}^{(q)}$ only if $j\neq i$ is among the first $q$ neighbors of $i$ . Let $\xi$ be the probability to sample the neighbor of a point from the same manifold, and $1-\xi$ the probability to sample it from a different manifold, with $\xi>0.5.$ Define $n_{i}^{in}$ as the number of neighbors of $i$ with the same manifold membership ( $n_{i}^{in}\equiv\sum_{j}\mathcal{N}_{ij}^{(q)}\delta_{z_{i}z_{j}}$ ). Then we introduce a probability distribution for $\boldsymbol{\mathcal{N}}^{(q)}$ as:

[TABLE]

where $\mathcal{Z}_{q}$ is a normalization factor that depends also on the sizes of the manifolds (see SI for its explicit expression). This term favors homogeneity within a $q$ -neighborhood of each point. With this addition, the model now reads:

[TABLE]

The posterior (5) is sampled with Gibbs sampling starting from a random configuration of the parameters. We repeat the sampling $M$ times starting from different random configurations of the parameters, keeping the chain with the highest maximum log-posterior value. The parameters $\boldsymbol{d}$ and $\boldsymbol{p}$ can be estimated by their posterior averages. For all estimations, we only include the last $10\%$ of the points. This ensures that the initial burn-in period is excluded [28]. Moreover, visual inspection of the MCMC output is performed to ensure that we do not incur in any label switching issues [29]. Indeed, we observed that once the chains converge to one mode of the posterior distribution they do not transition to different modes. As for the $\boldsymbol{z}$ , we estimate the value of $\pi_{ik}\equiv P_{post}(z_{i}=k)$ . Point $i$ can be safely assigned to manifold $k$ if $\pi_{ik}>0.8$ , otherwise we will consider its assignment to be uncertain. Note that our definition of “manifold” requires compactness: two spatially separated regions with approximately the same dimension are recognized as distinct manifolds by our approach.

We name our method Hidalgo (Heterogeneous Intrinsic Dimension Algorithm). Hidalgo has three free parameters: the number of manifolds, $K$ ; the local homogeneity range, denoted by $q$ ; the local homogeneity level, denoted by $\xi$ .We fix $q$ and $\xi$ based on preliminary tests conducted on several artificial data sets, which show that the optimal “working point” of the method is at $q=3,\xi=0.8$ (SI, fig. S1). As is common for mixture models, the value of $K$ can be estimated by model selection. In particular, we can compare the models with $K$ and $K+1$ for increasing $K$ , starting from $K=1$ and stopping when there is no longer significant improvement, as measured by the average log-posterior.

The computational complexity of the algorithm depends on the total number of points $N$ , the number of manifolds $K$ , and the parameter $q$ . As input to the Gibbs sampling, the method requires the neighbor matrix and the $\boldsymbol{\mu}$ vector. To compute the latter, one must identify the first $q$ neighbors of each point, which, starting from the distance matrix, requires $O(qN^{2})$ steps (possibly reducible to $O(qN\log N)$ with efficient nearest-neighbor-search algorithms [30]). The computation of a distance matrix from a coordinate matrix takes as $O(DN^{2})$ steps, where $D$ is the total number of coordinates. In all practical examples with $N\lesssim 10^{6}$ points these preprocessing steps require much less computational time than the Gibbs sampling. The complexity of the Gibbs sampler is dominated by the sampling of the $\boldsymbol{Z}$ , which scales linearly with $N$ and $K$ , and quadratically with $q$ (in total, $KNq^{2}$ ). The number of iterations required for convergence cannot be simply estimated. In practice, we found this time to be roughly linear in $N$ (SI, Fig S2). Thus, in practice the total complexity of Hidalgo is roughly quadratic in $N$ . The code implementing Hidalgo is freely available at https://github.com/micheleallegra/Hidalgo.

Results

Validation of the method on artificial data

We first test Hidalgo on artificial data for which the true manifold partition of the data is known. We start from the simple case of two manifolds with different IDs, $d_{1}$ and $d_{2}$ . We consider several examples, varying the higher dimension $d_{1}$ from $5$ to $9$ while fixing the lower dimension $d_{2}$ to $4$ . On both manifolds $N=1000$ points are sampled from a multivariate Gaussian with a variance matrix given by $1/d_{i}$ , $i=1,2$ multiplied by the identity matrix of proper dimension. The two manifolds are embedded in a space with dimension corresponding to the highest dimension $d_{2}$ , with their centers at a distance of $0.5$ standard deviations, so they are partly overlapping. In Fig. 1a-b we illustrate the results obtained in the case of fixed $\xi=0.5$ , equivalent to the absence of any statistical constraint on neighborhood uniformity (note that for $\xi=0.5$ the parameter $q$ is irrelevant). The estimates of the two dimensions are shown together with the normalized mutual information (NMI) between the estimated assignment of points and the true one. The latter is defined as the MI of the assignment and the ground truth labeling divided by the entropy of the ground truth labeling. As expected, without a constraint on the assignment of neighbors, the method is not able to correctly separate the points and thus to estimate the dimensions of the two manifolds, even in the case of quite different IDs (NMI $<0.001$ ). As soon as we consider $\xi>0.5$ , results improve. A detailed analysis of the influence of the hyperparameters $q$ and $\xi$ is reported in SI. Based on such analysis, we identify the optimal parameter choice as $q=3,\xi=0.8$ . In Fig. 1c-d we repeat the same tests as in 1a-b but with $q=3$ and $\xi=0.8$ . Now the NMI between the estimated and ground truth assignment is almost $1$ in all cases and, correspondingly, the estimation of $d_{1}$ and $d_{2}$ is accurate. To verify whether our approach can discriminate between more than two manifolds ( $K>2$ ), we consider a more challenging scenario consisting of five Gaussians with unitary variance in dimensions $1,2,4,5,9$ respectively. Some of the Gaussians have similar IDs, as in the case of dimensions $1$ and $2$ , or $4$ and $5$ : Moreover, they can be very close to each other: the centers of those in dimensions $4$ and $5$ are only half a variance far from each other, and they are crossed by the Gaussian in dimension $1$ . To analyze such dataset we again choose the hyperparameters $q=3$ and $\xi=0.8$ . We do not fix the number of manifolds $K$ to its ground truth value $K=5$ , but we try to let the method estimate $K$ without relying on a priori information. We perform the analysis with different values of $K=1,\dots,6$ and compute an estimate of the average of the logarithm of the posterior value $\mathcal{L}$ for each $K$ . Results are shown in Fig. 1e. We see that $\mathcal{L}$ increases up to $K=5$ , and then decreases, from which we infer that the optimal number of manifolds is $K=5$ . In Fig. 1f we illustrate the final assignment of points to the respective manifolds together with the estimated dimensions, upon setting the number of manifolds to $K=5$ . The separation of the manifolds is very good. Only a few points of the manifold with dimension $1$ are incorrectly assigned to the one with dimension $2$ and vice versa. The values of normalized mutual information between the ground truth and our classification is 0.89.

Finally, we show that our method can handle data embedded in curved manifolds, or with complex topologies. Such data should not pose additional challenges to the method: Hidalgo is sensitive only to the local structure of data. In terms of the distribution of the first and second nearest-neighbor distances, a 2-d torus looks exactly like a 2-d plane, and thus we should be able to correctly classify these objects, that are so topologically different, as 2-d objects regardless of their shape. We generated a dataset with 5 Gaussians in dimension 1, 2, 4, 5, 9 (Fig. 2a). Each Gaussian is then embedded on a curved nonlinear manifold: a 1-dimensional circle, a 2-dimensional torus, a 4-dimensional Swiss roll, a 5-dimensional sphere, and a 9-dimensional sphere. None of these manifolds is topologically isomorphic to a hyperplane, except for the Swiss roll. Moreover, the manifolds are intersecting. As shown in Fig. 2b, our approach can distinguish the five manifolds. The dimensions are correctly retrieved (Fig. 2b), points are assigned to the right manifold with an accuracy corresponding to a value of the normalized mutual information (NMI) $=$ 0.84. Only some points of the 3d-torus are misassigned to the 4d-swiss roll. As in the previous example, $\mathcal{L}$ increases up to $K=5$ , and then decreases, from which we infer that the optimal number of manifolds is $K=5$ ( Fig. 2c).

ID variability in a protein folding trajectory

As a first real application of Hidalgo, we address ID estimation for a dynamical system. Asymptotically, dynamical systems are usually restrained to a low-dimensional manifold in phase space, called an attractor. Much effort has been devoted to characterizing the ID of such attractor [4]. However, in the presence of multiple metastable states an appropriate description of the visited phase space may require the use of multiple IDs. Here, we consider the dynamics of the villin headpiece (PDB entry: 2F4K). Due to its small size and fast folding kinetics, this small protein is a prototypical system for molecular dynamics simulations. Our analysis is based on the longest available simulated trajectory of the system from Ref. [33]. During the simulated $125$ $\mu s$ , the protein performs approximately 10 transitions between the folded and the unfolded state. We expect to find different dimensions in the folded and unfolded state, since these two states are metastable, and they would be considered as different attractors in the language of dynamical systems. Moreover, they are characterized by different chemical and physical features: the folded state is compact and dry in its core, while the unfolded state is swollen, with most of the residues interacting with a large number of water molecules. We extract the value of the 32 key observables (the backbone dihedral angles) for all the $N=31,000$ states in the trajectory and apply Hidalgo to this data of extrinsic dimension $D=32$ . We obtain a vector of estimated intrinsic dimensions $\boldsymbol{d}$ and an assignment of each point $i$ to one of the $K$ manifolds. We find four manifolds, three low-dimensional ones ( $d_{1}=11.8$ , $d_{2}=12.9$ , $d_{3}=13.2$ ) and a high-dimensional one ( $d_{4}=22.7$ ). We recall that two spatially separated regions with approximately the same dimension (in this case, $d_{2}$ and $d_{3}$ ) are recognized as distinct manifolds by our approach. To test whether this partition into manifolds is related to the separation between the folded and the unfolded state we relate the partition to the fraction of native contacts $Q$ , which can be straightforwardly estimated on each configuration of the system. $Q$ is close to one only if the configuration is folded, while it approaches zero when the protein is unfolded. In Fig. 3 we plot the probability distribution of $Q$ restricted to the four manifolds. We find that the vast majority of the folded configurations ( $Q>0.8$ ) are assigned to the high-dimensional manifold. Conversely, the unfolded configurations ( $Q\leq 0.8$ ) are most of the times assigned to one of the low-dimensional manifolds. This implies that a configuration belonging to the low dimensional manifolds is almost surely unfolded. The normalized mutual information between the manifold assignment and the folded/unfolded classification is 0.60. Thus, we can essentially identify the folded state using the intrinsic dimension, a purely topological observable unaware of any chemical detail.

ID variability in time-series from brain imaging

In the next example, we analyze a set of time-series from functional resonance imaging (fMRI) of the brain, representing the BOLD (blood oxygen-level dependent) signal of each voxel, which captures the activity a small part of the brain [34]. fMRI time-series are often projected on a lower dimension through linear projection techniques like PCA [35], a step that assumes a uniform ID. However, the gross features of the signal (e.g., power spectrum and entropy) are often highly variable in different parts of the brain, and also non-uniformities in the ID may well be present. Here, we consider a single-subject fMRI recording containing $D=202$ images collected while a subject was performing a visuomotor task [36, 31]. From the images we extract the $N=29851$ time series corresponding to the BOLD signals of each voxel. Applying Hidalgo, we find two manifolds with very different dimensions $d_{1}=16.2$ , $d_{2}=31.9$ . Again, we relate the identified manifolds to a completely independent quantity, the clustering frequency $\Phi$ introduced in [32, 31], which measures the temporal coherence of the signal of a voxel with the signals of other voxels in the brain. Voxels with non-negligible clustering frequency ( $\Phi>0.2$ ) are likely to belong to brain areas involved in the cognitive task at hand. In Fig. 4 we show the probability distribution of $\Phi$ restricted to the two manifolds. We find that the “task-related” voxels ( $\Phi>0.2$ ) almost invariably belong to the manifold with high dimensionality. These voxels appear concentrated in the occipital, parietal, and temporal cortex (Fig. 4 bottom), and belong to a task-relevant network of coherent activity [31]. The normalized mutual information between the manifold assignment and the task-relevant/task-irrelevant classification is 0.24; note that this value is lower than for the protein folding case, because the the low-dimensional manifold is “contaminated” by many points with low $\Phi$ . This result finds an interesting interpretation: the subset of “relevant” voxels gives rise to patterns that are not only coherent but also characterized by a larger ID than the remainder of the voxels. On the contrary, the incoherent voxels exhibit a lower ID, hence a reduced variability, which is consistent with the fact that the corresponding time series are dominated by low-dimensional noise. Again, this feature emerges from the global topology of the data, revealed by our ID analysis, without exploiting any knowledge of the task that the subject is performing.

ID variability in financial data

Our final example is in the realm of economics. We consider firms in the well-known Compustat database ( $N=8309$ ). For each firm, we consider $D=31$ balance sheet variables from the fiscal year 2016 (for details, see Table S2 in the SI). Applying Hidalgo we find four manifolds of dimensions $d_{1}=5.4$ , $d_{2}=6.3$ , $d_{3}=7.0$ and $d_{4}=9.1$ . To understand this result, we try to relate our classification with common indexes showing the type and financial stability of a firm. We start by relating our classification to the Fama-French classification [37], which assigns each firm to one of twelve categories depending on the firm’s trade. In Fig. 5(top) we separately consider firms belonging to the different Fama-French classes, and compute the fraction of firms assigned to each the four manifolds identified by Hidalgo. The two classifications are not independent, since the fractions for different Fama-French classes are highly non-uniform. More precisely, the normalized mutual information between the two classifications is $0.19$ , rejecting hypothesis of statistical independence ( $p$ -value $<10^{-5}$ ). In particular, firms in the utilities and energy sector show a preference for low dimensions ( $d_{1}$ and $d_{2}$ ), while firms purchasing products (nondurables, durables, manufacturing chemicals, equipment, wholesale) are concentrated in the manifold with the highest dimension $d_{4}$ . The manifold with intrinsic dimension $d_{3}$ mostly includes firms in the financial and health care sectors. Different dimensions are not only related to the classification of the firm, but also their financial robustness. We consider the S&P quality ratings for the firms assigned to each manifold (also from Compustat; ratings are available for $2894$ firms). In Fig. 5(bottom) we show the distribution of ratings for the different manifolds. These distributions appear to be different. We cannot predict the rating based on the manifold assignment (the normalized mutual information between the manifold assignment and the rating is 0.04), but companies with worse ratings show a preference to belong to low-dimensional manifolds: by converting ratings into numerical values scores from $1$ to $8$ , we found a Pearson correlation of $0.22$ between the local ID and the rating ( $p<10^{-10}$ ) We suggest a possible interpretation for this phenomenon: a low ID may imply a more rigid balance structure, which may entail a higher sensitivity to market shocks which, in turn, may trigger domino effects in contagion processes. This result shows that information about the S&P rating can be derived using only the topological properties of the data landscape, without any in-depth financial analysis. For example, no information on the commercial relationship between the firms or the nature of their business is used.

Comparison with other methods

We have compared our method with two state-of-the-art methods in the literature, considering the synthetic datasets as well as the protein folding dataset (for which the folded/unfolded classification yields an approximate “ground truth”). The first method we consider is the “sparse manifold clustering and embedding” (SMCE, Ref. [19]). The method creates a graph by connecting neighboring points, supposedly lying on the same manifold, and then uses spectral clustering on this graph to retrieve the manifolds. As an output, it returns a manifold assignment for each point, together with an ID estimate. We resort to the implementation provided by Ehsan Elhamifar (http://khoury.neu.edu/home/eelhami/codes.htm). For SMCE, we can compute the NMI between the manifold assignment and the ground-truth manifold label. Furthermore, we can report the estimated local ID for points with a given ground-truth assignment, identified with the ID of the manifold to which they are assigned. SMCE requires to specify the number of manifolds $K$ as input and depends on a free parameter $\lambda$ . To present a fair comparison, we fixed $K$ at its ground truth value and explored different values in the range $0.01\leq\lambda\leq 20$ , and we report the results corresponding to the highest NMI with the ground truth. Furthermore, we repeated the estimation $M=10$ times and kept the results with the highest NMI with the ground truth. The second method is the local ID (LID) estimation by Carter et al. [10], which combines ID estimation restricted to a (large) neighborhood and local smoothing to produce an integer ID estimate for each point. We used the implementation provided by Kerstin Johnsson [38] (https://github.com/kjohnsson/intrinsicDimension) using a neighborhood size $k=50$ which is a standard value. Note that the LID estimation is deterministic, so we do not need to repeat the estimation several times. For LID, we can compute the NMI between the (integer) label given by the local ID estimate and the ground truth; moreover, we can compute the average ID, $\langle d\rangle$ for the points with a given ground truth assignment.

The results of the comparison are presented in Table 1. For the examples with two manifolds, SMCE [19] is able to correctly retrieve the manifolds with high accuracy (NMI $\geq 0.94$ ). The estimate of the lower dimension $d_{2}=4$ is always correct, while the highest dimension is underestimated for $d_{1}\geq 7$ . Overall, SMCE performs very well (even better than Hidalgo) on the datasets with two manifolds, as the two manifolds are well separated, matching well the assumptions of the model. On the examples with five manifolds, the performance of SMCE is worse, as the method cannot correctly retrieve intersecting manifolds. In the example with five linear manifolds, SMCE cannot retrieve the manifold of dimension $d_{1}=1$ , which is intersecting the manifolds of dimension $d_{3}=4$ and $d_{4}=5$ . Moreover, the ID estimates are quite poor. In the example with five curved manifolds, SMCE merges the manifolds with $d_{4}=5$ and $d_{5}=9$ , and cannot retrieve the manifolds with $d_{1}=1$ and $d_{3}=4$ , which get split. Again, the dimension estimates are poor. For the protein folding data, where regions with different IDs are probably not well separated, SMCE finds a single manifold with $d=11$ , and thus it is not able to detect any ID variation.

For the examples with two manifolds, the LID method [10] is generally able to correctly retrieve regions with different IDs with high accuracy (NMI $\geq 0.92$ ), except for the case with $d_{1}=5,d_{2}=4$ where results are inaccurate. The estimate of the lower dimension $d_{2}=4$ is always precise, while the higher dimension is always underestimated. The degree of underestimation increases with increasing $d_{1}$ . This is expected, since the ID estimates are produced by assuming a uniform density in a large neighborhood of each point, an assumption that fails in these data, especially for points on the manifold with higher dimension. Overall, LID performs quite well (comparably with Hidalgo) in separating regions with different IDs in the datasets with two manifolds, but gives worse ID estimates than Hidalgo. In the examples with five manifolds, the performance of LID is slightly worse. In the example with five linear manifolds LID merges the manifolds of dimension $d_{3}=4$ and $d_{=}5$ , and misassigns some of their points to the manifold of dimension $1$ . This is because relatively large neighborhoods are used for the estimation of the local ID: this leads to a difficulty in discriminating close regions with similar ID ( $d_{3}=4$ and $d_{4}=5$ ) and to a “contamination” of results by points of the manifold with $d_{1}=1$ intersecting the two. In the example with five curved manifolds, LID can correctly identify the manifolds with $d_{2}=2$ , $d_{4}=5$ , and $d_{5}=9$ , even if it underestimates the ID in the latter two cases (high dimension, large density variations). Many points in the manifold of dimension $d_{3}=4$ are mis-assigned to the manifold of dimension $d_{1}=1$ , due to the fact that the two manifolds are highly intersecting. Finally, for the protein folding data, where the density is highly variable, LID shows a tendency for unfolded configurations to have a higher local ID than unfolded ones. However, it cannot discriminate between folded and unfolded configurations as Hidalgo. The ID of the folded configurations is highly underestimated compared to the Hidalgo estimate, probably because of density variations.

In Figure 6 we compare the results of the three methods (Hidalgo, SMCE, and LID) for the dataset with five curved manifolds.

Discussion

The increasing availability of a large amount of data has considerably expanded the opportunities and challenges for unsupervised data analysis. Often data come in the form of a completely uncharted “point cloud” for which no model is at hand. The primary goal of the analyst is to uncover some structure within the data. For this purpose, a typical approach is dimensionality reduction, whereby the data are simplified by projecting them onto a low-dimensional space.

The appropriate intrinsic dimension (ID) of the space onto which one should project the data is not constant everywhere. In this work, we developed an algorithm (Hidalgo) to find manifolds of different IDs in the data. Applying Hidalgo, we observed large variations of the ID in datasets of diverse origin (a molecular dynamics simulation, a set of time series from brain imaging, a dataset of firm balance sheets). This finding suggests that a highly non-uniform ID is not an oddity, but a rather common feature. Generally speaking, ID variations can be expected whenever the system under study can be in various “macrostates” characterized by a different number of degrees of freedom due, for instance, to differences in the constraints. As an example, the folded state the protein is able to locally explore the phase space in many independent directions, while in the unfolded state it only performs wider movements in fewer directions. In the case of companies, a financially stable company may have more degrees of freedom to adjust its balance sheet.

In the cases we have analyzed, regions characterized by different dimensions were found to host data points differing in important properties. Thus, variations of the ID within a dataset can be used to classify the data into different categories. It is remarkable that such classification is achieved by looking at a simple topological property, the local ID. Let us stress that ID-based segmentation should not be considered an alternative to clustering methods. In fact, in most cases (e.g., most classical benchmark sets for clustering) different clusters do not correspond to significantly different intrinsic dimensions - rather, they correspond to regions of high density of points within a single manifold of well-defined ID. In such cases, clusters cannot be identified as different manifolds by Hidalgo. Conversely, when the ID is variable, regions of different IDs do not necessarily correspond to clusters in the standard meaning of the word: they may or may not correspond to regions of high density of points. A typical example could be that of protein folding: while the folded configurations are quite similar to one another, and hence constitute a cluster in the traditional sense, the unfolded configurations may be very heterogeneous, hence quite far in data space and then not a “cluster” in the standard meaning of the word.

The idea that ID may vary in the same data is not new. In fact, many works have discussed the possibility of a variable ID and developed methods to estimate multiple IDs [7, 8, 9, 10, 11, 12, 13, 21, 22, 14, 15, 16, 17, 18, 20]. Our method builds on these previous contributions but is designed with the specific goal of overcoming technical limitations of other available approaches and make ID-based segmentation a general-purpose tool. Our scheme uses only the distances between the data points, and not their coordinates, which significantly enlarges its scope of applicability. Moreover, the scheme uses only the distances between a point and its $q$ nearest neighbors, with $q\leq 5$ . We thus circumvent the notoriously difficult problem of defining a globally meaningful metric [3], only needing a consistent metric on a small local scale. Finally, Hidalgo is not hindered by density variations or curvature. For these reasons, Hidalgo is competitive with other manifold learning and variable ID estimation methods and, in particular, can yield better ID estimates and manifold reconstruction. We have compared our method with two state-of-the-art methods in the literature, the local ID method (LID, [10]), and the sparse manifold clustering and embedding (SMCE, [19]) method. Both methods show issues in the cases of intersecting manifolds and variable densities and yield worse ID estimates than Hidalgo, especially for large IDs.

Hidalgo is computationally efficient and therefore suitable for the analysis of large data sets. For example, it takes $\approx$ 30 mins to perform the analysis of the neuroimaging data, which includes 30000 data points, on a standard computer using a single core. Implementing the algorithm with parallel Gibbs sampling [39] may considerably reduce computing time and then yield a fully scalable method. Obviously, Hidalgo has some limitations. Some are intrinsic to the way the data are modeled. Hidalgo is not suitable to cover cases in which the ID is a continuously varying parameter [21], or in which sparsity is so strong that points cannot be assumed to be sampled from a continuous distribution. In particular, Hidalgo is not suitable for discrete data for which the basic assumptions of the method are violated. For example, ties in the sample could lead to null distances between neighbors, jeopardizing the computation of $\boldsymbol{\mu}$ . Others are technical issues related with the estimation procedure, and, and least in principle, susceptible to improvement in refined versions of the algorithm: for instance, one may improve the estimation with suitable enhanced sampling techniques. Finally, let us point out some further implications of our work. Our findings suggest a caveat with respect to common practices of dimensionality reduction, which assume a uniform ID. In the case of significant variations, a global dimensionality reduction scheme may become inaccurate. In principle, the partition in manifolds obtained with Hidalgo may be the starting point for using standard dimensionality reduction schemes. For example, one can imagine to apply PCA[1] or Isomap[3], or sketchmap[40] separately to each manifold. However, we point out that a feasible scheme to achieve this goal does not come as an immediate byproduct of our method. Once a manifold with given ID is identified, it is highly nontrivial to provide a suitable parameterization thereof, especially because the manifolds may be highly nonlinear, and even topologically non-trivial. How to suitably integrate our approach with a dimensionality reduction scheme remains a topic for further research. Another implication is that a simple topological invariant, the ID, can be very a powerful tool for unsupervised data analysis, lending support to current efforts at characterizing topological properties of the data [41, 42].

Data availability

Figure 1, 2, 6 and 7 have raw simulated data. Figure 3,4,5 use real data. All data used are publicly available in the following github repository: https://github.com/micheleallegra/Hidalgo

Code availability

Code is publicly available in the following github repository: https://github.com/micheleallegra/Hidalgo

Acknowledgements

MA was supported by Horizon 2020 FLAG-ERA (European Union), grant ANR-17-HBPR-0001-02 during completion of this work. We thank Giovanni Barone Adesi and Julia Reynolds (Institute of Finance, USI, Lugano, Switzerland) for helping us with the Compustat dataset and offering precious suggestions for its analysis. We thank Giulia Sormani (SISSA, Trieste, Italy) for suggesting, and providing us with the protein dynamics data. We thank Alex Rodriguez (SISSA, Trieste, Italy) for helping with the analysis of protein data. M.A. thanks the SISSA community at large for the substantial moral and intellectual support received, which has been critical for completion of this work.

The authors declare that they have no competing financial interests.

Correspondence and requests for materials should be addressed to A.L. (email: [email protected]) and A.M.. (email: [email protected]).

Supplementary Information

Enforcing neighborhood uniformity. In our model, we wish to obtain well separated manifolds. We do not wish to impose this condition in the form of a rigid constraint, since in real cases regions with different IDs are not completely separated, but only as a “soft constraint”, privileging configurations of $\boldsymbol{z}$ such that the first neighbors of each point are preferentially assigned to the same manifold. In a Bayesian framework, this means that given that $j$ is among the first neighbors of $i$ , the probability that $z_{i}=z_{j}$ is increased. Consider the neighbor matrix $\mathcal{N}_{ij}^{(q)}$ defined as:

[TABLE]

Intuitively, we would like to impose

[TABLE]

However, Eq. S2 is a relation between posterior probabilities, hence it cannot be directly embedded in the likelihood. What we can specify in the likelihood is the probability of observing the data $\mathcal{N}_{ij}^{(q)}$ , given an assignment $\boldsymbol{z}$ of the points. The way to enforce Eq. S2 is assuming that the first neighbors of each point are preferentially points of the same manifold. Consider the $i$ -th row of the neighbor matrix, $\boldsymbol{\mathcal{N}}_{i}^{(q)}\equiv\{\mathcal{N}_{ij}^{(q)},j=1,\dots,N\}$ . $\boldsymbol{\mathcal{N}}_{i}^{(q)}$ is a vector containing $q$ ones and $N-q$ zeros. Without any assumption, all configurations of $q$ zeros and $N-q$ ones are equally likely. Instead, we assume that neighbors are preferentially points from the same manifold. Formally, we assume that neighbors are selected forn the points of the same manifold with probability $\xi$ and from a different manifold with probability $1-\xi$ , with $\xi>1/2$ . Correspondingly, we introduce a new term in the likelihood:

[TABLE]

where

[TABLE]

is the number of neighbors of $i$ sampled from the same manifold, and

[TABLE]

is the number of neighbors of $i$ sampled from a different manifold. Function $\mathcal{Z}$ is a normalization factor that depends on $\xi$ :

[TABLE]

and can be expressed in a compact way as

[TABLE]

where ${}_{2}F_{1}(a,b,c,x)$ is the hypergeometric function. The derivation of this expression and the details about the likelihood term in (S3) are presented below. By considering all points $i$ , we obtain the global likelihood

[TABLE]

where

[TABLE]

is the total number of “internal” neighbors of points from manifold $k$ and

[TABLE]

is the total number of “external” neighbors of points from $k$ . Note that since $\mathcal{Z}$ depends on $i$ only through the hidden variables $\boldsymbol{z}$ we are able to split the product into $K$ components.

With this additional term in the likelihood, we obtain

[TABLE]

Derivation of the neighborhood uniformity term. With reference to $\boldsymbol{\mathcal{N}}_{i}^{(q)}$ , without any assumption, all configurations containing $q$ zeros and $N-q$ ones are equally likely. It is easy to compute the number of such configurations. The problem is analogous to the problem of selecting $q$ balls from a box containing $N-1$ balls: we have to choose $q$ neighbors among $N-1$ points, point $i$ being excluded. The number of possible choices is $\binom{N-1}{q}$ . Hence, all configurations of $\boldsymbol{\mathcal{N}}_{i}^{(q)}$ being equally likely we would have

[TABLE]

Instead, we assume that the neighbors of a point are preferentially points from the same manifold. Formally, we assume that neighbors are selected with probability $\xi$ among the $N_{z_{i}}$ points assigned to the same manifold of $i$ , and with probability $1-\xi$ among the $N-N_{z_{i}}$ points assigned to a different manifold. Here $\xi>1/2$ , so that configurations with neighbors assigned to the same manifold are more likely. Now, the problem is analogous to the problem where we have to select $q$ balls from two boxes, a black box containing $N_{b}$ balls and a white one containing $N_{w}$ balls. Before selecting each ball, we choose the box, the black one with probability $\xi$ and the white one with probability $1-\xi$ . Clearly, the probability of a choice with $n_{b}$ black and $q-n_{b}$ white balls is then proportional to $\xi^{n_{b}}(1-\xi)^{q-n_{b}}$ . For a given $n_{b}$ , the number of possible choices of balls is

[TABLE]

One can easily verify that $\sum_{n_{b}=0}^{q}\binom{N_{b}}{n_{b}}\binom{N_{w}}{q-n_{b}}=\binom{N_{b}+N_{w}}{q}$ . The probability of a given choice is then

[TABLE]

where $\mathcal{Z}=\sum_{n_{b}=0}^{q}\binom{N_{b}}{n_{b}}\binom{N_{w}}{q-n_{b}}\xi^{n_{b}}(1-\xi)^{q-n_{b}}$ . By using the formula (Abramowitz and Stegun, 15.4.1)

[TABLE]

where $(a)_{n}=a(a+1)\dots(a+n-1)$ is the Pochammer symbol and doing some simple algebra, $\mathcal{Z}$ can be compactly expressed as

[TABLE]

Replacing $N_{b}$ with $N_{z_{i}}-1$ (the number of points assigned to the same manifold as $i$ , excluding $i$ ), $N_{w}$ with $N-N_{z_{i}}$ (the number of points assigned to a different manifold), and $n_{b}$ with $n_{i}^{in}$ , we obtain the likelihood of a given configuration of $\boldsymbol{\mathcal{N}}_{i}^{(q)}$ as

[TABLE]

Choice of the free parameters. In order to find a good configuration for the parameters $(q,\xi)$ , we perform tests with several values of $q\in\{2,3,\dots,10\}$ and $\xi\in[0.5,1)$ for all the artificial datasets used for validation: i) two manifolds of dimension $4$ and $5$ , with a Gaussian density of points; ii) five manifolds of dimension $1,2,4,5,9$ , with a Gaussian density of points; iii) five manifolds of dimension $1,2,4,5,9$ , with a Gaussian density of points, embedded on a 1-dimensional circle, a 2-dimensional torus, a 4-dimensional Swiss roll, a 5-dimensional sphere and a 9-dimensional sphere.

The crucial figure of merit to assess the performance of the method is the normalized mutual information (NMI) between the estimated and the true assignment of $\boldsymbol{z}$ , which measures the quality of the assignment of points to manifolds. Indeed, once the manifolds are correctly separated, the problem is essentially reduced to a dimension estimation within the single manifolds (which is successfully solved by TWO-NN). In Fig. S1 we show the NMI as a function of $(q,\xi)$ . In all cases, we performed $10^{5}$ iterations of the Gibbs sampling and repeated the sampling $M=10$ times starting from different random configurations of the parameters. We kept the sampling with the highest maximum average of the log-posterior value.

For all values of $q$ , the MI first increases and then decreases with $\xi$ . This can be expected based on the following considerations. When $\xi$ is close to $0.5$ , as we discussed above, the method cannot discriminate between different manifolds. When $\xi$ is increased, the posterior distribution starts to prefer configurations that approximately satisfy the neighborhood homogeneity constraint. For sufficiently high $\xi$ , the posterior distribution is sharply peaked at the configuration that optimally satisfies this constraint; correspondingly, if the Gibbs sampler can explore the parameter space exhaustively, it will eventually find this peaked region and remain trapped there. Hence, the NMI achieves average values close to 1. However, for $\xi$ close to $1$ , the posterior distribution is very likely to also have pronounced local maxima. Therefore, depending on the initial configuration, the sampler may remain trapped in one of them. Hence, one can observe a drop in the NMI. In general, these sampling issues can be worsened when $q$ is increased since the local maxima become more and more pronounced. In principle, this problem may be dealt with by resorting to well established enhanced sampling techniques. For simplicity, in the present work, we prefer to verify that there is a region of the parameter space where the results appear optimal and restrict to these regions for subsequent analyses.

In general, $q=2$ yields poor results ( $NMI\lesssim 0.4$ ) in all cases. This means that is general $q=2$ is too low to effectively enforce the uniformity constraint. Analogously, $q=10$ too high yields poor results ( $NMI\lesssim 0.5$ ), except in the simple case of two manifolds. There are two independent reasons for this behavior. First, $q=10$ can lead to sharp peaks in the posterior, and hence to sampling issues; second, $q=10$ gives issues in the case of intersecting manifolds, as it enforces uniformity on too large a scale. Good results ( $NMI\gtrsim 0.75$ ) are found for $q=2,0.75\leq\xi\leq,0.80,0.85$ and for $q=5,0.65\leq\xi\leq 0.70$ . Based on these findings, we identify the optimal “working point" of the method at $q=3,\xi=0.8$ , which yields $NMI>0.80$ ) for all data sets considered.

Convergence time as a function of sample size.

In order to test the speed of convergence of the Gibbs sampling as a function of the number of points $N$ , we have generated random data sampled from two manifolds of dimension $d_{1}=9$ and $d_{2}=4$ . As in Results, $N/2$ points are sampled from a multivariate Gaussian with a variance matrix given by $1/d_{i}$ , $i=1,2$ , times the identity matrix of proper dimension. The two manifolds are embedded in a space with dimension corresponding to the higher dimension $d_{2}$ , with their centers at a distance of $0.5$ standard deviations. We consider different sample sizes $N=2^{i},\{i=6,15\}$ . In order to test convergence, we use the criterion in [43]. Given independent samples $Y_{t}$ (with $t=1,\dots,T$ and $T$ the total number of samples), and given a function $f_{t}=f(Y_{t})$ , convergence can be verified from the behavior of the cumulative sums

[TABLE]

where $\mu=\frac{1}{T-t_{0}}\sum_{\tau=t_{0}+1}^{T}f_{\tau}$ , and $t_{0}\geq 1$ is any time, usually chosen as higher to the putative burn-in time. At convergence, the behavior of $F_{t}$ is highly non-smooth; In particular, consider the indicator variable

[TABLE]

which indicates whether $S_{t}$ has a local maximum/minimum at $t$ . At convergence, approximately half of the points should be local minima or maxima. Formally, $e_{t}$ behaves as a Bernoulli variable $\mathcal{B}(1/2)$ , and consequently

[TABLE]

behaves as a binomial variable $B(t-t_{0}-1,1/2)$ . Therefore, one should have

[TABLE]

for at least $(1-\alpha)\%$ of the time (for large $t$ , the binomial will be approximately normal). In other words, when convergence is reached, $E_{t}$ should fluctuate around $1/2$ with a variance compatible with that of a binomial variable.

To adopt this test to the data, we studied the point-wise dimension estimates $d_{1t}$ and $d_{2t}$ . We subsampled the $d_{it}$ at $t=1,\theta,2\theta,\dots$ with $\theta=10$ to ensure independence of the samples. Let $D_{1t}$ and $D_{2t}$ be the resulting samples. We used $f_{1\tau}=d_{1t}$ and $f_{2\tau}=d_{2t}$ and first verified whether the corresponding $E_{1t}$ , $E_{2t}$ (Eq. S13) converged to $1/2$ . While an excessive variance around $1/2$ is a sign of residual non-stationarity, failure of $E_{t}$ to approach the value $1/2$ for increasing $t$ is a likely sign of the non-independence of the samples used [43]. Therefore, if $E_{it}$ did not approach to $1/2$ , we increased the subsampling at $\theta=20,\theta=50,\theta=100,\theta=200,\theta=500$ until $\langle E_{it}\rangle\simeq 1/2$ for large $t$ . We divided the $D_{1t}$ and $D_{2t}$ in non-overlapping windows of $20$ successive points and checked whether $E_{it}$ was within the bounds (S14) with $\alpha=0.05$ for 95% of the times (at 19 times out of 20) within each window. We then defined as convergence time $T_{conv}$ the time corresponding to the first window where both $E_{1t}$ , $E_{2t}$ were satisfying the criterion.

For each dataset size $N$ , we repeated the assessment of $T_{conv}$ for $M=20$ independent Gibbs samplings. The results are shown in Fig. S2: one can see that the convergence time increases roughly linearly with data size $N$ .

Restricting the analysis to points with no common first and second neighbor.

Compustat variables used.

Bibliography43

The reference list from the paper itself. Each links out to its DOI / PubMed record.

1[1] Jolliffe, I. T. Principal component analysis and factor analysis. In Principal component analysis , 115–128 (Springer, 1986).
2[2] Roweis, S. T. & Saul, L. K. Nonlinear dimensionality reduction by locally linear embedding. \Journal Title science 290 , 2323–2326 (2000).
3[3] Tenenbaum, J. B., De Silva, V. & Langford, J. C. A global geometric framework for nonlinear dimensionality reduction. \Journal Title science 290 , 2319–2323 (2000).
4[4] Grassberger, P. & Procaccia, I. Measuring the strangeness of strange attractors. In The Theory of Chaotic Attractors , 170–189 (Springer, 2004).
5[5] Levina, E. & Bickel, P. J. Maximum likelihood estimation of intrinsic dimension. In Advances in neural information processing systems , 777–784 (2005).
6[6] Rozza, A., Lombardi, G., Ceruti, C., Casiraghi, E. & Campadelli, P. Novel high intrinsic dimensionality estimators. \Journal Title Machine learning 89 , 37–65 (2012).
7[7] Barbará, D. & Chen, P. Using the fractal dimension to cluster datasets. In Proceedings of the sixth ACM SIGKDD international conference on Knowledge discovery and data mining , 260–264 (ACM, 2000).
8[8] Gionis, A., Hinneburg, A., Papadimitriou, S. & Tsaparas, P. Dimension induced clustering. In Proceedings of the eleventh ACM SIGKDD international conference on Knowledge discovery in data mining , 51–60 (ACM, 2005).