Manifold-adaptive dimension estimation revisited

The FSA and mFSA algorithm

There is a dataset with a sample size n, and sample points x_i ∈ ℝ^m. Then,

1. Compute distances: Calculate the distance of the kth and 2kth nearest neighbors (R_k, R_2k) for each data point (x_i). Here the neighborhood size is some positive integer k ∈ ℤ⁺.

2. Compute local estimates: Get local estimates δ_k(x_i) from the distances for each data point according to Eq. (5).

3. Calculate global estimate: Aggregate the local estimates into one global value. This last step is the only difference between the FSA and the mFSA method:

   1. FSA estimator: (7)${\displaystyle d_{FSA}^{\left(k\right)}=\frac{\sum {\delta }_k\left(x_i\right)}{n}}$

   2. mFSA estimator: (8)${\displaystyle d_{mFSA}^{\left(k\right)}=\text{M}\left[\left\{{\delta }_k\left(x_1\right),{\delta }_k\left(x_2\right),\dots ,{\delta }_k\left(x_n\right)\right\}\right]}$ where M stands for the sample median.

The cmFSA algorithm

There is a dataset with a sample size n, and sample points x_i ∈ ℝ^m. Then,

1. Compute mFSA estimate Apply the mFSA algorithm to get biased global estimate $d_{mFSA}^{\left(k\right)}$.

2. Model Calibration Fit a correction-model with the the given sample size n on uniform random hypercube calibration datasets consisting of various intrinsic dimension values, many instances each (at least N = 15 realizations). We used the following model: (9)${\displaystyle D\approx dexp\left(\sum _{l=1}^L{\alpha }_l{d}^l\right)}$ where D is the true dimension of the underlying manifold, α_l-s are sample size and k dependent coefficients, L is the order of the polynomial and $d=d_{mFSA}^{\left(k\right)}$ is a shorthand for the biased local estimate. This model is derived from heuristic reasoning, and simplifies to a linear model in the parameters, if the logarithm of the two sides is taken. First we calculate biased estimates on each test data. Second, we carry out the model fit by linear regression on the log–log values with the ordinary least squares method or with orthogonal distance regression.

3. Calculate cmFSA estimate Plug in the biased estimate into fitted the correction model to compute $d_{cmFSA}^{\left(k\right)}$.

A python implementation of the algorithms can be found at https://github.com/phrenico/cmfsapy along with the supporting codes for this article.

Simulations on D-hypercube datasets

The simulations were implemented in python3 (Van Rossum & Drake, 2009) using the numpy (Oliphant, 2006), scipy (Virtanen et al., 2020) and matplotlib (Hunter, 2007) packages, unless otherwise stated.

We generated test-datasets by uniform random sampling from the unit D-cube to demonstrate, that theoretical derivations fit to data. We measured distances with a circular boundary condition to avoid edge effects, hence the data is as close to the theoretical assumptions as possible.

To illustrate the probability density function of the FSA estimator, we computed the local FSA intrinsic dimension values (Fig. 2). We generated d-hypercubes (n = 10, 000, one realization) with dimensions of 2, 3, 5, 8, 10 and 12, then computed histograms of local FSA estimates for three neighborhood sizes: k = 1, 11, 50 respectively (Figs. 2A, 2F). We selected these specific neighborhoods because of didactic purposes: the k = 1 neighborhood is the smallest one, the k = 50 is a bigger neighborhood, which is still much smaller than the sample size, so the estimates are not affected by the finite sample effect. The k = 11 neighborhood represents a transition between the two “extremes”, the specific value is an arbitrary choice giving pleasing visuals suggesting the gradual change in the shape of the curve as a function of the k parameter. We drew the theoretically computed probability density function (pdf) to illustrate the fit.

To show that the theoretically computed sampling distribution of the mFSA fits to the hypercube datasets, we varied the sample size (n = 11, 101, 1001) with N = 5, 000 realizations from each. We computed the global mFSA for each realization and plotted the results for d = 2 (Fig. 3A) and d = 5 (Fig. 3B).

We investigated the dimensionality and sample-size effects on mFSA estimates (Figs. 4–5). We simulated the hypercube data in the 2–30 dimension-range, and applied various sample sizes: n = 10, 100, 1, 000, 2, 500, 10, 000; one hundred realizations each (N = 100). We computed the mFSA values with minimal neighborhood size (k = 1), and observed finite-sample-effects, and asymptotic convergence. We repeated the analysis with hard boundary conditions.

We fitted a correction formula on the logarithm of dimension values and estimates with the least squares method (Eq. 10), using all 100 realizations for each sample sizes separately (Fig. 6). (10)${\displaystyle \alpha =\frac{\sum \left(ln{E}_i\right)d^{\left(i\right)}}{\sum {\left(d^{\left(i\right)}\right)}^2}}$ where E_i = D_i/d⁽ⁱ⁾ is the relative error, D_i is the intrinsic dimension of the data, and d⁽ⁱ⁾ are the corresponding mFSA estimates. We carried out the model fit on the 2–30 intrinsic dimension range.

We also calibrated the cmFSA algorithm in a wider range of intrinsic dimension values (2–80) and applied more coefficients in the polynomial fit procedure (Fig. S1A). Also, we used orthogonal distance regression to fit the mean over realizations of lnE_i with the same D_i value (Fig. S1B). We utilized the mean and standard deviation of the regression error to compute the ideal error rate of cmFSA estimator, if the error-distributions are normal (Figs. S1C–S1F).

Simulations on customly sampled manifolds

We carried out simulations on datasets sampled from manifolds according to uniform, multivariate Gaussian, Cauchy distribution and on uniformly sampled D-spheres in the function of sample size as in Facco et al. (2017), Fig. 2.

The uniform sampling was carried out on D-hypercube data with periodic boundary conditions. The Gaussian datasets were sampled from a zero mean and unit variance and no covariance multivariate normal distribution. The Cauchy datasets were generated so that the probability density of the norms were a Cauchy distribution. We achieved this by the following procedure:

1. Generate n points according to D dimensional Gaussian distribution (ζ_i) and normalize the euclidean distance of the points from the origin. ${\displaystyle z_i=\frac{{\zeta }_i}{\left|{\zeta }_i\right|}\text{where}{\zeta }_i\sim \mathcal{N}\left(0,I\right)}$ and I is the D-dimensional identity matrix. Thus, the points z_i are uniformly distributed on the hyper-surface of a D − 1 dimensional hyper-sphere of unit radius.

2. Generate n positive real numbers u_i from a Cauchy distribution $f\left(u\right)=\frac{1}{\pi \left(1+u^2\right)}$ and multiply z_i by this to get a dataset: ${\displaystyle x_i=u_i\times z_i}$ Thus the norms of the resulting points are distributed according to a Cauchy distribution.

Finally, we produced the D-sphere data with the first step of the previous procedure.

We generated N = 200 instances of each dataset with the intrinsic dimension values D = 2, 5, 10, we estimated the global mFSA and cmFSA dimensions and plotted the mean and standard deviation in the function of sample size (Fig. 7).

Comparison on synthetic benchmark datasets

We simulated N = 100 instances of 15 manifolds (Table 1, M_i, n = 2, 500) with various intrinsic dimensions. We generated the datasets according to the first 15 manifolds proposed by Campadelli et al. (2015). More specifically, Table 1 contains the description manifold types, the first 15 manifolds of Table 2 are used in this work as synthetic benchmark, and the Table 4 shows the benchmark results in Campadelli et al. (2015), http://www.mL.uni-saarland.de/code/IntDim/IntDim.htm.

We applied the wide (D =2–80) calibration procedure (l₁ =  − 1, l₂ = 1, l₃ = 2, l₄ = 3) as in the previous subsection (n = 2, 500, k = 5) to compute cmFSA for the datasets. We used cmFSA in two modes, in integer and in fractal mode. In the former the global estimates are rounded to the nearest integer value, while in the latter case the estimates can take on real values.

We measured the performance of the mFSA and corrected-mFSA estimators on the benchmark datasets, and compared them with the performance of ML (Levina & Bickel, 2004) DANCo (Ceruti et al., 2014) and the 2NN (Facco et al., 2017) (Table 2) estimators. We used the Matlab (MATLAB, 2020; Lombardi, 2020) (see on github) and an R package (Johnsson, Soneson & Fontes, 2015) implementation of DANCo. In the case of DANCo, we also investigated the results for integer and for fractal mode just as for the cmFSA algorithm.

To quantify the performance we adopted the Mean Percentage Error (MPE, Eq. 11) metric (Campadelli et al., 2015): (11)${\displaystyle \mathrm{MPE}=\frac{100}{MN}\sum _{j=1}^M\sum _{i=1}^N\frac{|D_j-d_{ij}|}{D_j}}$ where there is N realizations of M types of manifolds, D_j are the true dimension values, d_ij are the dimension estimates.

Also, we used the error rate—the fraction of cases, when the estimator did not find (missed) the true dimensionality—as an alternative metric (Fig. 8). We used this metric to compare the performace of DANCo and cmFSA in integer mode, we simply counted the cases, when the estimator missed the true dimension value: (12)${\displaystyle H_j=\frac{1}{N}\sum _{i=1}^{N}I\left(D_j\ne d_{ij}\right)}$ where H_j is the error rate for a manifold computed from N realizations and I = 1 if D_j ≠ d_ij is the indicator function for the error. We computed the mean error rate H by averaging the manifold specific values.

Dimension estimation of interictal and epileptic dynamics

We used data of intracranial field potentials from two subdural grids positioned –parietofrontally (6*8 channels, Gr A-F and 1–8) and frontobasally (2*8 channels, Fb A-B and 1–8) –on the brain surface and from three strips located on the right temporal cortex (8 channels, JT 1–8), close to the hippocampal formation and two interhemispheric strips, located within the fissura longitudinalis, close to the left and right gyrus cinguli (8 channels BIH 1–8 and 8 channels JIH 1–8) as part of presurgical protocol for a subject with drug resistant epilepsy (Fig. 9A). The participants signed a written consent form and the study was approved by the relevant institutional ethical committee (Medical Research Council, Scientific and Research-Ethics Committee TUKEB, Ref number: 20680-4/2012/EKU (368/PI/2012)). This equipment recorded extracellular field potentials at 88 neural channels at a sampling rate of 2048 Hz. Moreover, we read in—using the neo package (Garcia et al., 2014)—selected 10 second long chunks of the recordings from interictal periods (N = 16) and seizures (N = 18) to further analysis.

We standardised the data series and computed the Current Source Density (CSD) as the second spatial derivative of the recorded potential. We rescaled the 10 second long signal chunks by subtracting the mean and dividing by the standard deviation. Then, we computed the CSD of the signals by applying the graph Laplacian operator on the time-series. The Laplacian contains information about the topology of the electrode grids, to encode this topology, we used von Neumann neighborhood in the adjacency matrix. After CSD computation, we bandpass-filtered the CSD signals (Gramfort et al., 2013) (1–30 Hz, fourth order Butterworth filter) to improve signal to noise ratio.

We embedded CSD signals and subsampled the embedded time series. We used an iterative manual procedure to optimize embedding parameters (Fig. S2). Since the fastest oscillation is (30 Hz) in the signals, a fixed value with one fourth period (2048/120 ≈ 17 samples) were used as embedding delay. We inspected the average space–time separation plots of CSD signals to determine a proper subsampling, with the embedding dimension of D = 2 (Fig. S2A). We found, that the first local maximum of the space–time separation was at around 5 ms: 9–10, 10–11, 11–12 samples for the 1%, 25%, 50% percentile contour-curves respectively. Therefore, we divided the embedded time series into 10 subsets to ensure the required subsampling. Then, we embedded the CSD signal up to D = 12 and measured the intrinsic dimensionality for each embeddings (Figs. S2B and S2C). We found that intrinsic dimension estimates started to show saturation at D ≥ 3, therefore we chose D = 7 as a sufficiently high embedding dimension (averaged over k =10–20 neighborhood sizes).

We measured the intrinsic dimensionality of the embedded CSD signals using the mFSA method during interictal and epileptic episodes (Fig. 9). We selected the neighborhood size between k = 10 and k = 20 and averaged the resulting estimates over the neighborhoods and subsampling realizations. We investigated the dimension values (Figs. 9C and 9D) and differences (Fig. 9E) between interictal and epileptic periods.

We also compared the mFSA estimates with the original –mean based –FSA estimates in the function of neighborhood size on a recording in the k =1–12 neighborhood range and plotted the estimates against each other to visualize differences (Fig. 9B).

Manifold adaptive dimension estimator revisited

Results on randomly sampled hypercube datasets

Theoretical probability density function of the local FSA estimator fits to empirical observations (Eq. (15), Fig. 2). We simulated hypercube datasets with fixed sample size (n = 10, 000) and of various intrinsic dimensions (D = 2, 3, 5, 8, 10, 12). We measured the local FSA estimator at each sample point with three different k parameter values (k = 1, 11, 50). We visually confirmed that the theoretical pdf fits perfectly to the empirical histograms.

The empirical sampling distribution of mFSA fits to the theoretical curves for small intrinsic dimension values (Fig. 3). To demonstrate the fit, we drew the density of mFSA on two hypercube datasets D = 2 and D = 5 with the smallest possible neighborhood (k = 1), for different sample sizes (n = 11, 101, 1,001). At big sample sizes the pdf is approximately a Gaussian (Laplace, 1986), but for small samples the pdf is non-Gaussian and skewed.

The mFSA estimator underestimates intrinsic dimensionality in high dimensions. This phenomena is partially a finite sample effect (Fig. 4), but edge effects make this underestimation even more severe. This phenomenon was pronounced at low sample sizes and high dimensions, but we experienced convergence to the real dimension value as we increased sample size.

We graphically showed that mFSA estimator asymptotically converged to the real dimension values for hypercube-datasets, when we applied periodic boundary conditions (Fig. 5). We found, that the convergence is much slower for hard boundary conditions, where edge effects make systematic estimation errors higher.

From the shape of the curves in Fig. 4, we heuristically derived a correction formula for finite sample size and edge effects (Eq. 9). The heuristics is as follows. We tried to find a formula, which intuitively describes the true intrinsic dimension in the function (C) of the estimated values. One can see on Fig. 4, that at small values the error converges to zero and also the curve lies approximately on the diagonal, so it’s derivative goes to one. (30)${\displaystyle \underset{d\to 0}{lim}C\left(d\right)=D\underset{d\to 0}{lim}{C}^{\prime }\left(d\right)=1}$ where C is the correction function and d is the biased estimate. Equation (9) satisfies these conditions and gives good fit to empirical data (Fig. 6).

From an other point of view, Eq. (9) means that one could estimate the logarithm of relative error with an L-order polynomial: (31)${\displaystyle log\left(E_{rel}\right)=log\left(\frac{D}{d}\right)=\sum _{l=1}^L{\alpha }_l{d}^l}$

The order of the polynomial was different for the two types of boundary conditions. When we applied hard boundary, the order was L = 1, however in the periodic case higher order polynomials fit the data. Thus, in the case of hard-boundary, we could make the empirical correction formula: (32)${\displaystyle D\approx C\left(d\right)=d{e}^{{\alpha }_{n}d}}$ where α_n is a sample size dependent coefficient that we could fit with the least squares method. This simple model described well the data in the 2–30 intrinsic dimension range (Figs. 6A–6F).

Results on customly sampled manifolds

We investigated the case when the assumption of uniform sampling or flatness is violated through gaussian, Cauchy and hypersphere datasets (Fig. 7) with various intrinsic dimensions and sample sizes. We added hypercube datasets with periodic boundary conditions as a control with the same parameter setting respectively (k = 5).

On the hypercube datasets with periodic boundary conditions the mFSA algorithm produced a massive underestimation of intrinsic dimension for low sample sizes for D = 10, but cmFSA corrects this bias caused by finite sample size (Fig. 7A). For the small-dimensional cases. When D = 2 and D = 5 both cmFSA and FSA estimated well the true intrinsic dimension values. On the gaussian datasets with non-periodic boundary conditions mFSA produced even more severe underestimation for D = 10 or D = 5, but cmFSA overestimated the intrinsic dimensions (Fig. 7B). On the heavy tailed Cauchy datasets mFSA showed a non-monotonic behaviour in the function of sample size: for fewer points it had low values with a maximum at mid sample sizes and exhibited slow decline convergence to true dimension value for big samples (Fig. 7C). This shape of the curve resulted in underestimation for small samples followed by an overestimation part depressing towards the true dimension values as N goes to infinity (D = 5, 10). For D = 2 the first underestimation section was missing. cmFSA produced severe overestimation for these Cauchy datasets. The hypersphere dataset is an example when the point density is approximately uniform, but the manifold is curved (Fig. 7D). On this datset mFSA produced underestimation for D = 5, 10 and good estimates for D = 2. cmFSA overestimated the dimension value.

Results on synthetic benchmarks

We tested the mFSA estimator and its corrected version on synthetic benchmark datasets (Hein & Audibert, 2005; Campadelli et al., 2015). We simulated N = 100 instances of 15 manifolds (Table 1, M_i, n = 2, 500) with various intrinsic dimensions.

We estimated the intrinsic dimensionality of each sample and computed the mean, the error rate and Mean Percentage Error (MPE) for the estimators. We compared the mFSA, cmFSA, the R and the Matlab implementation of DANCo, the Levina-Bickel and the 2NN estimator (Table 2). cmFSA and DANCo was evaluated in two modes, in a fractal-dimension mode and in an integer dimension mode.

The mFSA estimator underestimated intrinsic dimensionality, especially in the cases when the data had high dimensionality. The Levina-Bickel estimator overestimated low intrinsic dimensions and underestimated the high ones. The 2NN estimator produced underestimation on most test manifolds it reached the best average result on the M₆ and M₁₃ manifolds.

In contrast, the cmFSA estimator found the true intrinsic dimensionality of the datasets, it reached the best overall error rate (0.277) and 2nd best MPE (Fig. 8, Table 2). In some cases, it slightly over-estimated the dimension of test datasets. Interestingly, DANCo showed implementation-dependent performance, the Matlab algorithm showed the 2nd best error rate (0.323) and the best MPE value (Table 2). The R version overestimated the dimensionality of datasets in most cases.

Analysing epileptic seizures

To show how mFSA works on real-world noisy data, we applied it to human neural recordings of epileptic seizures.

We acquired field potential measurements from a patient with drug-resistant epilepsy by 2 electrode grids and 3 electrode strips. We analyzed the neural recordings during interictal periods and during epileptic activity to map possible seizure onset zones (see Methods).

We found several characteristic differences in the dimension patterns between normal and control conditions. In interictal periods (Fig. 9C), we found the lowest average dimension value at the FbB2 position on the fronto-basal grid. Also, we observed gradually increasing intrinsic dimensions on the cortical grid (Gr) between the F1 and D6 channels. In contrast, we observed the lowest dimension values at the right interhemispherial strip (JIH 1–2) and on the temporo-basal electrode strip (JT 3–5) close to the hippocampus, and the gradient on the cortical grid altered during seizures (Fig. 9D). Comparing the dimensions between seizure and control periods, the majority of the channels showed lower dimensions during seizures. This decrease was most pronounced close to the hippocampal region (strip JT) and the parietal region mapped by the main electrode grid (GrA-C). Curiously, the intrinsic dimensionality became higher at some frontal (GrE1-F2) and fronto-basal (FbA1-B3) recording sites during seizure (Figs. 9A and 9E).

Comparison of the original FSA and the mFSA dimension estimators on the seizure data series showed characteristic difference similar to the one observed in the simulated data: mFSA resulted in lower dimension estimates than FSA and the difference between the two methods decreases as the k neighbourhood increases (Fig. 9B, compare it with Figs. 1C and 1D).