# Effective estimation of the dimension of a manifold from random samples

Lucien Grillet  UNIV RENNES, IRMAR - UMR 6625, F-35000 RENNES, FRANCE
<EMAIL_ADDRESS>and  Juan Souto  UNIV RENNES, CNRS, IRMAR - UMR 6625,
F-35000 RENNES, FRANCE<EMAIL_ADDRESS>
###### Abstract.

We give explicit theoretical and heuristical bounds for how big does a data
set sampled from a reach-$1$ submanifold $M$ of euclidian space need to be, to
be able to estimate the dimension of $M$ with $90\%$ confidence.

##  1\. Introduction

The manifold hypothesis asserts that naturally ocurring data sets
$X\subset\mathbb{R}^{s}$ behave as if they had been sampled from a low-
dimensional submanifold $M$. There is then a large literature on manifold
learning, that is on understanding how properties of the underlying manifold
$M$ can be ”learned” from the data set $X$—see for example [3, 15, 22, 21, 5,
20, 7]—, but it is maybe fair to say that estimating the dimension $\dim(M)$
of the manifold, that is figuring out the intrinsic dimension of the data set,
seems to be of particular interest. This is maybe so because the dimension is
one of the most basic quantities associated to a manifold, but maybe also
because of the importance of understanding the applicability to the data set
$X$ of different dimension reduction schemes. In any case, there are very
numerous estimators for the intrinsic dimension of a data set. We refer to the
surveys [9, 10, 7] for a brief discussion of many of those estimators and to
the references therein for details—see also Chapter 3 in the monograph [20].
Anyways, all these dimension estimators are based on the idea that manifolds
look locally like their tangent spaces, that is like euclidean space. For
example, non-linear PCA aims at finding the dimension of tangent spaces. ANOVA
aims at computing the average angle between vectors in a tangent space,
exploiting the fact that for euclidean space itself one can read the dimension
from the expected value. The estimator (1.1) below, as well as closely related
algorithms due to Takens [25, 26], Theiler [27] or Grassberger-Procaccia [18],
exploits the fact that the correlation dimension of the manifold is just its
dimension. There are other such estimators where the correlation dimension is
replaced by the Box counting dimension or the Kolmogorov capacity. This list
of estimators is far from being all inclusive.

All those estimators work very well—in the sense that they yield with high
probability the desired result—as long as we apply them with common sense and
under suitable conditions, that is if

  * 1)

we work at a scale at which the manifold $M$ resembles euclidean space, and

  * 2)

the data set $X$ is large enough and the points therein have been obtained by
uniform sampling.

However, amazingle little is known when it comes to quantifying 1) or 2).
Indeed, there are classical estimates due to Grassberger [17], Procaccia [23],
and Eckmann-Ruelle [13] giving absolute lower bounds, in terms of the
dimension, for the number of measurements needed for the results to be
reliable: Grassberger argues that there are absolute lower bounds at all,
Procaccia argues that if the intrindic dimension is $\dim$ then one might need
a data set of at least $10^{\dim}$ measurements, while Eckmann-Ruelle suggest
a lower bound of the form $C\cdot 10^{\frac{1}{2}\dim}$ for some large but
undetermined $C$. One can summarize these results as asserting that, when the
dimension grows, the minimal cardinality of a data set allowing to compute the
dimension grows exponentially—this is indeed already the case (see [28]) when
one just wants to distinguish between spheres of consecutive dimensions.

The absolute lower bounds for the data size that we just mention do not help
however with the following less philosophical question: I suspect, or
hypothesize, that my data set is sampled out of a manifold with this or that
properties. To what extent can I trust the intrinsic dimension estimation
given by this or that estimator? The only result we know along those lines is
due to Niyogi-Smale-Weinberger [22], at least as long as one wants to allow
for variable curvature manifolds. In [22] the authors give namely an algorithm
to compute the homology of a closed submanifold $M\subset\mathbb{R}^{s}$ out
of a set $X$ sampled from $M$, and they estimated how large does the data size
have to be so that their algorithm has success in at least $90\%$ of the
cases—observe that knowing the homology we also know the dimension of the
manifold. A problem is that the needed data size is astronomical. For example,
the estimate in [22] for the number of points needed to be sampled to compute
with $90\%$ probability of success the homology of the 4-dimensional Clifford
torus
$\mathbb{T}^{4}=\mathbb{S}^{1}\times\mathbb{S}^{1}\times\mathbb{S}^{1}\times\mathbb{S}^{1}$
is 24.967.788 points, that is about 25 million points.

One might contend that the Niyosi-Smale-Weinberger algorithm computes
something much more sophisticated than the dimension, that they are basically
learning the whole manifold, and that being able to do that is an overkill if
what one wants to do is to estimate its dimension. We agree with that point of
view. We will consider the estimator

(1.1)
$\dim_{\operatorname{Corr}(\varepsilon_{1},\varepsilon_{2})}(X)=\operatorname{Round}\left(\frac{\log|PX(\varepsilon_{1})|-\log|PX(\varepsilon_{2})|}{\log(\varepsilon_{1})-\log(\varepsilon_{2})}\right)$

where

(1.2) $PX(\varepsilon)=\big{\\{}\\{x,y\\}\subset X\text{ with
}0<|x-y|\leqslant\varepsilon\big{\\}},$

and our goal will be to estimate how large does $X\subset M$ need to be for
the estimator (1.1) to have a $90\%$ success rate. It is definitively much
smaller. For example, applying (1.1) to randomly sampled data sets in
$\mathbb{T}^{4}$ we have a $90\%$ rate of sucess, as long as we sample
$18.262$ points, that is more than 1300 times less than than before.

Evidently, the value of
$\dim_{\operatorname{Corr}(\varepsilon_{1},\varepsilon_{2})}(X)$ does not only
depend on $X$ but also on the chosen scales. In some sense the goal of this
paper is to decide how to choose $\varepsilon_{1}$ and $\varepsilon_{2}$ in
such a way that one does not need to sample too many points to be $90\%$ sure
that (1.1) returns the correct value. This is the kind of results that we will
prove:

######  Theorem 1.1.

For $d=1,\cdots,10$ let $\varepsilon_{1}$ and $\varepsilon_{2}$ be scales as
in the table below. Also, given a closed $d$-dimensional manifold
$M\subset\mathbb{R}^{s}$ with reach $\tau(M)\geqslant 1$ let $n$ be also as in
the following table:

$d$ | $\varepsilon_{1}$ | $\varepsilon_{2}$ | n
---|---|---|---
$1$ | $1.5$ | $0.19$ | $9+21\cdot\operatorname{vol}(M)^{\frac{1}{2}}$
$2$ | $0.78$ | $0.2$ | $94+58\cdot\operatorname{vol}(M)^{\frac{1}{2}}$
$3$ | $0.63$ | $0.23$ | $635+146\cdot\operatorname{vol}(M)^{\frac{1}{2}}$
$4$ | $0.54$ | $0.23$ | $2786+392\cdot\operatorname{vol}(M)^{\frac{1}{2}}$
$5$ | $0.46$ | $0.22$ | $7013+1119\cdot\operatorname{vol}(M)^{\frac{1}{2}}$
$6$ | $0.4$ | $0.21$ | $13221+3366\cdot\operatorname{vol}(M)^{\frac{1}{2}}$
$7$ | $0.36$ | $0.21$ | $25138+10644\cdot\operatorname{vol}(M)^{\frac{1}{2}}$
$8$ | $0.33$ | $0.2$ | $50033+34890\cdot\operatorname{vol}(M)^{\frac{1}{2}}$
$9$ | $0.31$ | $0.19$ | $63876+119533\cdot\operatorname{vol}(M)^{\frac{1}{2}}$
$10$ | $0.29$ | $0.18$ | $139412+425554\cdot\operatorname{vol}(M)^{\frac{1}{2}}$

Then, if we sample independently and according to the riemannian volume form a
subset $X\subset M$ consisting of at least $n$ points, then we have

$\dim_{\operatorname{Corr}(\varepsilon_{1},\varepsilon_{2})}(X)=d$

with at least $90\%$ probability.

Here the reach $\tau(M)$ (see Definition 1), or in the terminology of [22] the
condition number, is taken as a measure for the local regularity of the
submanifold $M$. It is evident that in order to have specific bounds for the
needed data size, we do need to have some a priori control on the local
geometry: otherwise we could have a 1-dimesional submanifold so interwoven
that it looks as a $d$-dimensional manifold for some $d\geqslant 2$. Taking
the reach as a measure to quantitify to which extent does a submanifold
$M\subset\mathbb{R}^{s}$ resemble euclidean space seems to be actually pretty
common in the field of manifold learning [1, 22, 15, 6].

After discussing briefly the correlation dimension and the estimator (1.1) in
Section 2 we discuss some aspects of the geometry of the reach in Section 3.
We will mostly care about the volume of the thick diagonal

$DM(\varepsilon)=\\{(x,y)\in M\times M\text{ with
}|x-y|\leqslant\varepsilon\\}$

For example, in Theorem 3.1 we use the Bishop-Gromov theorem and the
properties of CAT(1)-spaces to give upper and lower bounds for the ratio
$\frac{\operatorname{vol}(DM(\varepsilon_{1}))}{\operatorname{vol}(DM(\varepsilon_{2}))}$
for sufficiently small $\varepsilon_{1}>\varepsilon_{2}$ positive:

######  Theorem 1.2.

Suppose that $M\subset\mathbb{R}^{s}$ is a $d$-dimensional ($d\geqslant 1$)
closed submanifold with reach $\tau(M)\geqslant 1$. Then we have

$\frac{\frac{\varepsilon_{1}}{2}}{\arcsin(\frac{\varepsilon_{2}}{2})}\\!\left(\\!\frac{\sin(\varepsilon_{1})}{\sin(2\arcsin\frac{\varepsilon_{2}}{2})}\\!\right)^{d-1}\\!\\!\\!\\!\leqslant\frac{\operatorname{vol}(DM(\varepsilon_{1}))}{\operatorname{vol}(DM(\varepsilon_{2}))}\\!\leqslant\frac{\int_{0}^{\sqrt{2}\cdot
2\cdot\arcsin(\frac{\varepsilon_{1}}{2})}\\!\sinh^{d-1}(t)dt}{\int_{0}^{\sqrt{2}\cdot\varepsilon_{2}}\sinh^{d-1}(t)dt}$

for any two $1>\varepsilon_{1}>\varepsilon_{2}>0$.

In Section 4 we come then to core of the present pamphlet. The reason why we
care about the volume of the thick diagonal is that for $X\subset M$ we have

$\frac{\log|PX(\varepsilon_{1})|-\log|PX(\varepsilon_{2})|}{\log(\varepsilon_{1})-\log(\varepsilon_{2})}\sim\frac{\log(\operatorname{vol}(DM(\varepsilon_{1})))-\log(\operatorname{vol}(DM(\varepsilon_{2})))}{\log(\varepsilon_{1})-\log(\varepsilon_{2})}$

with large probability, at least if $X$ has been obtained by independently
sampling a large number of points according to the riemannian measure on $M$.
Basically the goal of Section 4, or maybe even the goal of this paper, is to
find $\varepsilon_{1},\varepsilon_{2}>0$ so that a relatively small set $X$ is
such that with high probability the left side lies in the interval
$(d-\frac{1}{2},d+\frac{1}{2})$. We prove Theorem 1.1 above in Section 4. The
scales in the table in the theorem are obtained numerically—the process is
also described in Section 4, and a computer implementation can be found in
[19].

### Heuristical bounds

It is more or less evident that the bounds given by Theorem 1.1 are far from
being sharp. In Section 5 we add a heuristical discussion of what could be, in
practice, more realistic bounds. The starting point is that the number of
sampled points should not be what determines how reliable is the obtained
result, but rather the number $|PX(\varepsilon_{1})|$ of pairs at our larger
scale $\varepsilon_{1}>\varepsilon_{2}>0$. Arguing as if

  1. (1)

at our scales all balls in $M$ were euclidean, and

  2. (2)

the distances between the two points in pairs as in (1.2) were independent,

we get that, with $N$ as in Table 1, it would suffice to have a data set
$X\subset M$ with at least $N$ pairs as in (1.2) for (1.1) to give the correct
answer in about $90\%$ (resp. $70\%$) of the cases.

$d$ | $\varepsilon_{1}$ | $\varepsilon_{2}$ |  $N$ for 90% |  $N$ for 70%
---|---|---|---|---
1 | $1.5$ | $0.19$ | $30$ | $10$
2 | $0.78$ | $0.2$ | $122$ | $40$
3 | $0.63$ | $0.23$ | $249$ | $111$
4 | $0.54$ | $0.23$ | $516$ | $238$
5 | $0.46$ | $0.22$ | $878$ | $360$
6 | $0.4$ | $0.21$ | $1329$ | $554$
7 | $0.36$ | $0.21$ | $1719$ | $698$
8 | $0.33$ | $0.2$ | $2481$ | $1070$
9 | $0.31$ | $0.19$ | $3900$ | $1604$
10 | $0.29$ | $0.18$ | $5849$ | $2414$
Table 1. Heuristic bound $N$ for how large should the cardinality of
$PX(\varepsilon_{1})$ at least be to have $90\%$ (resp. $70\%$) rate of
success when using (1.1).

Note that the assumptions (1) and (2) are not that outlandish, at least if
$\varepsilon_{1}$ is small and if the number of pairs is small when compared
to the cardinality of $X$. In any case, for what it is worth, numerical
simulations (see Table 5) seem to support the values given in Table 1. In
particular, our numerical simulations also indicate that if
$|PX(\varepsilon_{1})|$ is less or equal than the value in the right column in
Table 1, then we should count with about a $30\%$ failure rate when we use
(1.1).

######  Remark.

As a bigger value of $|PX(\varepsilon_{1})|$ results in better performance, it
is a natural idea to try to maximize it. This can be done in two ways. First,
we can try to obtain more data to increase the number $n$, but as this is not
always possible in concrete situations, the second solution is the increase
the scale $\varepsilon_{1}$. The downside of this second method is that
considering bigger balls increases the effect of the curvature. One must
carefully balance those two effects.

Thinking on the reliability of (1.1) in terms of the number of pairs has the
huge advantage that one does not need to have any a priori knowledge of the
volume of the manifold. It is however not clear how different are the
heuristic bound from Table 1 and the formal bound given in Theorem 1.1. To be
able to compare both bounds we note that, always under the assumption that all
$\varepsilon_{1}$ balls in our $d$-manifold $M$ are euclidean, then if we
sample a set $X\subset M$ with $n$ points then we expect to have

$|PX(\varepsilon_{1})|=\frac{n(n-1)}{2}\cdot\frac{\operatorname{vol}(B^{\mathbb{R}^{d}}(\varepsilon_{1}))}{\operatorname{vol}(M)}$

pairs of points. Using this relation we get for example that, using the
heuristic bound, it would suffice to sample $1958$ points from the torus
$\mathbb{T}^{4}$ for (1.1) to be $90\%$ of the time correct. See Table 4 for
more on the comparisson between the heuristic bound and the bound in Theorem
1.1 for the number of points that suffice to have $90\%$ success rate when
using (1.1). This is, once again, supported by numerical experiments—see Table
6.

Note now that thinking of the applicability of (1.1) in terms of the
cardinality of $PX(\varepsilon_{1})$ leads to an implementation of (1.1) which
could be applied to data sets $X$ sampled from a manifold $M$ whose reach we
ignore. The basic idea is the following:

> If we want to check if $M$ has dimension $d$ then we take $\varepsilon_{1}$
> minimal so that $|PC(\varepsilon_{1})|$ is as large as Table 1 asks for in
> dimension $d$, then we choose $\varepsilon_{2}$ so that the ratio
> $\frac{\varepsilon_{1}}{\varepsilon_{2}}$ is as in Table 1), and then
> compute $\dim_{\operatorname{Corr}(\varepsilon_{1},\varepsilon_{2})}(X)$.

The preceeding discussion implies that for sufficiently rich synthetic data
sets this algorithm has a reliability of about 90%. For the sake of
completeness we decided to test it also on data sets each consisting in 200
grayscale pictures of the 3D-model Suzanne (see Figure 1) obtained by randomly
choosing 1, 2 and 3 Euler angles—the pictures are 64 by 64 pixels large and
can thus be represented as points in $\mathbb{R}^{4096}$. See Table 7 for the
obtained results.

Figure 1. 3 sample pictures from the data sets.

### Comparisson with other estimators

In Section 6 we compare briefly the performance of (1.1) with that of a few
other dimension estimators: reading the correlation dimension via a log-log
plot, ANOVA, and local PCA. The basic observation is that, at least if the
volume of $M$ greater or equal to that of the $d$-dimensional Clifford torus
and if we take sets of cardinality close to the heuristic bound proposed
either in Table 1 or in Table 4, then (1.1) seems to perform better than ANOVA
and local PCA, at least with the specific implementations we proposed. It is
however harder to compare any of these estimators with the incarnation of
(1.1) using a log-log plot: it is namely unclear how to meaningfully quantify
what does it mean for the non-constant slope of a curve to look constant. It
thus only seems to make sense to test the log-log technique numerically, but
it is also unclear how to formally do this: while one might well print
diagrams for 100 randomly chosen random subsets of a given manifold, it is not
clear how can one formally decide what any given diagram is suggesting to us.
Still, the limited experiments we conducted lead to the conclussion that the
log-log plot method is pretty reliable for relatively small data sets.

### Acknowledgements

The second author would like to thank Kaie Kubjas for getting him interested
in this topic.

##  2\. The estimator

In this section we recall what is the correlation dimension and how it can be
approximated to obtain dimension estimators, but first we introduce some
notation that we will use through out the paper. First, distances in eulidean
space will be denoted either by $d_{\mathbb{R}^{n}}(x,y)$ or by $|x-y|$. If
$x$ and $y$ are points in a Riemannian manifold $M$ then $d_{M}(x,y)$ is their
distance with respect to the Riemannian metric. The ball in $M$ of radius $R$
centered at some point $x\in M$ is denoted by $B^{M}(x,R)$, although
sometimes, when there is no risk of confussion, we drop the subscript writing
thus simply $B(x,R)$. Again, if there is no risk of confussion with the
dimension, balls in the ambient euclidean space $\mathbb{R}^{s}$ will be
denoted by blackboard bold, that is $\mathbb{B}(x,R)=B^{\mathbb{R}^{s}}(x,R)$.
Finally, the intersection of balls in the ambient space $\mathbb{R}^{s}$ with
the manifold $M$ will be denoted by

$\mathbb{B}^{M}(x,R)=M\cap\mathbb{B}(x,R)$

With this notation in place, we turn our attention to the correlation
dimension and the estimator (1.1).

As we mentioned earlier, there are plenty of dimension estimators (see [9, 10]
and [20, Chapter 3]). In the language of [10], (1.1) is a fractal-based
estimator. Fractal-based because what one is aiming at, is computing a
dimension which makes sense for fractal objects—in this case the correlation
dimension. Recall that the correlation integral at scale $\varepsilon$ of a
Borel measure $\mu$ on $\mathbb{R}^{s}$ is the integral

$C(\mu,\varepsilon)=\int\mu(\mathbb{B}(x,\varepsilon))\,d\mu(x)$

The upper correlation dimension and lower correlation dimension are then
defined as

$D^{+}(\mu)=\limsup_{\varepsilon\to
0}\frac{\log(C(\mu,\varepsilon))}{\log(\varepsilon)}\text{ and
}D^{-}(\mu)=\liminf_{\varepsilon\to
0}\frac{\log(C(\mu,\varepsilon))}{\log(\varepsilon)}.$

When both of them agree, then one refers to

$D(\mu)=D^{+}(\mu)=D^{-}(\mu)$

as the correlation dimension of $\mu$. There are numerous situations of
interest in dynamics [24] in which the correlation dimension of a measure
exists. It is also trivial that it exists if $\mu$ is a smooth measure whose
support is a submanifold of euclidean space. Since this is the case in which
we will find ourselves, we state this fact as a lemma:

######  Lemma 2.1.

If $\mu$ is a smooth finite measure of full support of a $d$-dimensional
submanifold $M\subset\mathbb{R}^{s}$, then the correlation dimension exists
and we have $D(\mu)=d$.∎

In [18] Grassberger and Procaccia note that if $\varepsilon>0$ is small enough
and if we have a set $X=(x_{1},\dots,x_{n})$ consisting of $n$ points ($n$
very large) sampled independently with respect to the measure $\mu$ then
$D(\mu)$ is approximated by

(2.3)
$\dim_{GP}(X,\varepsilon)=\frac{\log\left(\frac{|DX(\varepsilon)|}{n(n-1)}\right)}{\log(\varepsilon)}$

where

(2.4) $DX(\varepsilon)=\\{(x_{i},x_{j})\in X\times X\text{ with }i\neq j\text{
and }|x_{i}-x_{j}|\leqslant\varepsilon\\}$

is the set of (ordered) pairs of points in $X$ within $\varepsilon$ of each
other. Grassberger and Procaccia propopose (2.3) as an estimator for the
intrinsic dimension of the data set $X$—variants are discussed by Takens [25,
26] and Theiler [27]. We will however focus here on the version from [20], or
more precisely the quantity

(2.5)
$\dim_{\operatorname{Corr}(\varepsilon_{1},\varepsilon_{2})}(X)=\operatorname{Round}\left(\frac{\log|DX(\varepsilon_{1})|-\log|DX(\varepsilon_{2})|}{\log\varepsilon_{1}-\log\varepsilon_{2}}\right)$

Note that this expression is nothing other than (1.1) above.

A problem when using (2.5) as an estimator, or for that matter when we use
(2.3), is that we have to pick up the appropriate scales. In practial
implementations, this is often by-passed by taking many possible scales
$\varepsilon_{1}<\varepsilon_{2}<\dots<\varepsilon_{k}$, computing
$|DX(\varepsilon_{i})|$ for each one of those scales, plotting the result in a
log-log-plot and choosing what looks to us as the slope at some region where
the slope looks constant. In practice, at least when tested on synthetic data,
the log-log plot method is surprisingly effective, but it is very unclear how
one can get out of that a formal statistical test, or how can one evaluate how
much confidence can one have on the obtained number. Theorem 1.1 from the
introduction gives such bounds, at least if the manifold we are working with
is sufficiently regular, in the sense of having reach $\tau(M)\geqslant 1$. We
discuss a few aspects of the geometry of the reach in the next section.

##  3\. Some Geometry

In this section we recall a few facts about the geometry of reach-1
submanifolds $M$ of euclidean space. Combining these facts with standard
arguments from Riemannian geometry we give upper and lower bounds for the
volume of the $\varepsilon$ fat diagonal in $M\times M$.

###  3.1. Reach-1 manifolds

We start recalling the definition of the reach of a closed subset of euclidean
space:

######  Definition 1.

The reach $\tau(S)$ of a closed subset $S\subset\mathbb{R}^{s}$ is the
supremum of those $T\geqslant 0$ with the property that for every
$x\in\mathbb{R}^{s}$ with $d_{\mathbb{R}^{s}}(x,S)\leqslant T$ there is a
unique point in $S$ closest to $x$.

Since the reach was introduced by Federer in [14], it has proved to be a
useful notion. Indeed, it follows from the very definition that positive reach
sets, that is sets $S$ with $\tau(S)>0$, have some neighborhood
$\mathcal{N}(S)$ on which the closest point projection $\pi:\mathcal{N}(S)\to
S$ is well-defined—the existence of such a projection is enough to show that
sets of positive reach share many regularity properties with convex sets.

Here we will be working from the very beginning with very regular objects,
namely closed smooth submanifolds $M\subset\mathbb{R}^{s}$ of euclidean space.
In this setting the reach, or rather a lower bound for the reach, helps to
quantify how distorted is the inner geometry of $M$ with respect to that of
the ambient euclidean space. We summarize what we will need in the following
proposition:

######  Proposition 3.1.

Let $M\subset\mathbb{R}^{s}$ be a closed submanifold with reach
$\tau(M)\geqslant 1$. Then we have:

  1. (1)

$d_{M}(x,y)\leqslant 2\arcsin\left(\frac{|x-y|}{2}\right)$ for any two $x,y\in
M$ with $|x-y|<2$.

  2. (2)

The set $\mathbb{B}^{M}(x,r)=M\cap\mathbb{B}(x,r)$ is geodesically convex for
any $x\in\mathbb{R}^{s}$ and any $r<1$.

  3. (3)

We have $\angle\big{(}\gamma^{\prime}(0),\gamma^{\prime}(t)\big{)}\leqslant
d_{M}(\gamma(0),\gamma(\ell))$ for every geodesic $\gamma:[0,\ell]\to M$
parametrized by arc length.

  4. (4)

The manifold $M$ has sectional curvature pinched by
$-2\leqslant\kappa_{M}\leqslant 1$.

See, in that order Lemma 3, Corollary 1 and Lemma 5 in [6] for the first three
claims of Proposition 3.1. See then Proposition A.1 (iii) in [1] for the final
claim. In any case, we refer to [1] and [6], and to the references therein for
general facts about submanifolds $M\subset\mathbb{R}^{s}$ of positive reach.

###  3.2. Injectivity radius

Armed with Proposition 3.1 we can now derive a lower bound for the injectivity
radius of those submanifolds $M\subset\mathbb{R}^{s}$ with reach
$\tau(M)\geqslant 1$. Recall that the injectivity radius of a geodesically
complete Riemannian manifold $M$ at a point $x\in M$ is defined as

$\operatorname{inj}(M,x)=\sup\\{t>0\ |\ \exp_{x}:T_{x}M\to M\text{ is
injective on }B^{T_{x}M}(0,t)\\}$

where $B^{T_{x}M}(0,t)=\\{v\in T_{x}M\text{ with }|v|\leqslant r\\}$ and where
$\exp_{x}$ is the Riemannian exponential map. The injectivity radius of $M$
itself is then defined to be

$\operatorname{inj}(M)=\inf_{x\in M}\operatorname{inj}(M,x).$

Recall also that the systole $\operatorname{syst}(M)$ of a closed manifold $M$
is the length of the shortest non-trivial closed geodesic. The importance of
the systole now is that, together with an upper bound
$\kappa_{M}\leqslant\kappa$ for the sectional curvature, the systole yields
very a usable bound [4, Thm. 89] for the injectivity radius:

$\operatorname{inj}(M)\geqslant\min\left\\{\frac{\pi}{\sqrt{\kappa}},\frac{1}{2}\operatorname{syst}(M)\right\\}.$

Altogether, see [4] for basic facts and definitions from Riemannian geometry.

Suppose now that $M\subset\mathbb{R}^{s}$ is a closed submanifold of euclidean
space with $\tau(M)\geqslant 1$, and note the that assumption that $M$ is
closed implies that it is metrically complete and thus geodesically complete
by the Hopf-Rinow theorem. On the other hand if $\gamma:\mathbb{S}^{1}\to
M\subset\mathbb{R}^{s}$ is any smooth curve and if $t$ is such that
$|\gamma(t)-\gamma(0)|$ is maximal, then
$\langle\gamma^{\prime}(0),\gamma^{\prime}(t)\rangle=0$. We thus get from (3)
in Proposition 3.1 that every closed geodesic reaches at least distance
$\frac{\pi}{2}$, and hence that

$\operatorname{syst}(M)\geqslant\pi.$

Now, this fact together with the upper bound for the sectional curvature from
(4) in Proposition 3.1 and with the lower bound for the injectivity radius
yields the following:

######  Corollary 3.2.

If $M\subset\mathbb{R}^{s}$ is a closed submanifold of euclidean space with
reach $\tau(M)\geqslant 1$, then we have
$\operatorname{inj}(M)\geqslant\frac{\pi}{2}$.∎

It should be noted that if $\dim(M)=1$ and $\tau(M)\geqslant 1$ then
$\operatorname{syst}(M)\geqslant 2\pi$ and hence
$\operatorname{inj}(M)\geqslant\pi$. The standard circle shows that this is
optimal. We do not know by how much can one improve the bound given in
Corollary 3.2 in other dimensions.

###  3.3. Volumes of balls

Although we will be mostly interested in the volume of sets
$\mathbb{B}^{M}(x,r)=M\cap\mathbb{B}(x,r)$ we start by considering the volume
of actual metric balls in the manifold. First note that having bounds on the
curvature, we also get bounds on the volumes of balls, at least as long as the
radius remains below the injectivity radius. More precisely, if $M$ is a
Riemannian manifold of dimension $d\geqslant 2$ and with curvature pinched in
$[-2,1]$ then we have

(3.6)
$\operatorname{vol}(B^{\mathbb{S}^{d}}(r))\leqslant\operatorname{vol}(B^{M}(x,r))\leqslant\operatorname{vol}(B^{\frac{1}{\sqrt{2}}\mathbb{H}^{d}}(r))$

for all $r<\operatorname{inj}_{M}(x)$ (see [4, Thm. 103] and [4, Thm. 107]).
Here, the sphere $\mathbb{S}^{d}$ and the scaled hyperbolic space
$\frac{1}{\sqrt{2}}\mathbb{H}^{d}$ are respectively the simply connected
complete $d$-dimensional manifolds of constant curvature $1$ and $-2$. Having
bounds for the volumes of balls we also have bounds for the ratios between
volumes of balls of different radius. We can get however somewhat better
results:

######  Proposition 3.3.

Suppose that $M$ is a $d$-dimensional ($d\geqslant 2$) Riemannian manifold
with sectional curvature pinched in $[-2,1]$. We then have

$\frac{R}{r}\cdot\left(\frac{\sin(R)}{\sin(r)}\right)^{d-1}\leqslant\frac{\operatorname{vol}(B^{M}(R,x))}{\operatorname{vol}(B^{M}(r,x))}\leqslant\frac{\int_{0}^{\sqrt{2}\cdot
R}\sinh^{d-1}(t)\ dt}{\int_{0}^{\sqrt{2}\cdot r}\sinh^{d-1}(t)\ dt}$

for any two $0<r<R<\min\\{\pi,\operatorname{inj}(M)\\}$ and any $x\in M$.

###### Proof.

We begin with the upper bound. We get from the curvature bound
$\kappa\geqslant-2$ and the Bishop-Gromov comparison theorem [4, p.310] that

$\frac{\operatorname{vol}(B^{M}(x,R))}{\operatorname{vol}(B^{M}(x,r))}\leqslant\frac{\operatorname{vol}(B^{\frac{1}{\sqrt{2}}\mathbb{H}^{d}}(R))}{\operatorname{vol}(B^{\frac{1}{\sqrt{2}}\mathbb{H}^{d}}(r))}=\frac{\int_{0}^{\sqrt{2}\cdot
R}\sinh^{d-1}(t)\ dt}{\int_{0}^{\sqrt{2}\cdot r}\sinh^{d-1}(t)\ dt}$

for all $x\in M$, as we wanted.

Let us now deal with the lower bound. Well, the fact that $M$ has curvature
pinched from above by $1$ implies that $M$ is is locally a CAT($1$)-space—see
[8, 11] for facts about CAT($\kappa$)-spaces and comparisson geometry. Recall
now that we are working at a scale smaller than $\pi$ and the injectivity
radius. In particular, the CAT($1$) property implies that geodesic triangles
we encounter are thinner in our manifold than in $\mathbb{S}^{d}$. This
implies in particular that, for $0<r<R<\min\\{\pi,\operatorname{inj}(M)\\}$,
the radial projection

$\operatorname{proj}:S^{M}(x,R)\to S^{M}(x,r),\ \
\operatorname{proj}(y)=\exp_{x}(r\cdot R^{-1}\cdot\exp_{x}^{-1}(y))$

contracts distances more (expands distances less) than the corresponding map
in the sphere, meaning that $\operatorname{proj}$ is
$\frac{\sin(r)}{\sin(R)}$-Lipschitz. We deduce thus that

$\mathrm{Area}(S^{M}(x,r))\leqslant\left(\frac{\sin(r)}{\sin(R)}\right)^{d-1}\cdot\mathrm{Area}(S^{M}(x,R))$

where $S^{M}(x,r)=\partial B^{M}(x,r)$ is the distance-$r$-sphere in $M$ and
$\mathrm{Area}(\cdot)$ stands for the $(d-1)$-dimensional volume. Anyways, if
we set $T=\frac{r}{R}$ then we get from the co-area formula that

$\displaystyle\operatorname{vol}(B^{M}(x,r))$
$\displaystyle=\int_{0}^{r}\mathrm{Area}(S^{M}(x,t))\ dt$
$\displaystyle\stackrel{{\scriptstyle
s=\frac{1}{T}t}}{{=}}\int_{0}^{R}\mathrm{Area}(S^{M}(x,T\cdot s))\cdot T\ ds$
$\displaystyle\leqslant\frac{r}{R}\cdot\int_{0}^{R}\left(\frac{\sin(Ts)}{\sin(s)}\right)^{d-1}\mathrm{Area}(S^{M}(x,s))\
ds$

The function

$(0,\pi)\to\mathbb{R},\ s\to\frac{\sin(Ts)}{\sin(s)}$

is monotonically increasing (because $T=\frac{r}{R}\in(0,1)$). This means that

$\displaystyle\operatorname{vol}(B^{M}(x,r))$
$\displaystyle\leqslant\frac{r}{R}\cdot\left(\frac{\sin(r)}{\sin(R)}\right)^{d-1}\cdot\int_{0}^{R}\mathrm{Area}(S^{M}(x,s))\
ds$
$\displaystyle=\frac{r}{R}\cdot\left(\frac{\sin(r)}{\sin(R)}\right)^{d-1}\cdot\operatorname{vol}(B^{M}(x,R))$

And we are done. ∎

We come now to the result we really care about:

######  Corollary 3.4.

Suppose that $M\subset\mathbb{R}^{s}$ is a closed $d$-dimensional ($d\geqslant
1$) submanifold with reach $\tau(M)\geqslant 1$. We then have

$\operatorname{vol}(B^{\mathbb{S}^{d}}(r))\leqslant\operatorname{vol}(\mathbb{B}^{M}(x,r))\leqslant\operatorname{vol}\left(B^{\frac{1}{\sqrt{2}}\mathbb{H}^{d}}\left(2\arcsin\left(\frac{r}{2}\right)\right)\right)$

and

$\frac{\frac{R}{2}}{\arcsin(\frac{r}{2})}\\!\left(\\!\frac{\sin(R)}{\sin(2\arcsin\frac{r}{2})}\\!\right)^{d-1}\\!\\!\\!\\!\leqslant\frac{\operatorname{vol}(\mathbb{B}^{M}(x,R))}{\operatorname{vol}(\mathbb{B}^{M}(x,r))}\\!\leqslant\frac{\int_{0}^{\sqrt{2}\cdot
2\cdot\arcsin(\frac{R}{2})}\\!\sinh^{d-1}(t)dt}{\int_{0}^{\sqrt{2}\cdot
r}\sinh^{d-1}(t)dt}$

for any two $0<r<R<1$ and any $x\in M$.

###### Proof.

Let’s assume for the time being that $d\geqslant 2$. From part (1) in
Proposition 3.1 we get for any $t<2$ that

(3.7) $B^{M}(x,t)\subset\mathbb{B}^{M}(x,t)\subset
B^{M}\left(x,2\arcsin\left(\frac{t}{2}\right)\right)$

The first claim then follows directly from (3.6).

On the other hand, if we combine (3.7) with the upper bound in Proposition 3.3
we get that

$\displaystyle\frac{\operatorname{vol}(\mathbb{B}^{M}(x,R))}{\operatorname{vol}(\mathbb{B}^{M}(x,r))}$
$\displaystyle\leqslant\frac{\operatorname{vol}\left(B^{M}\left(x,2\arcsin\left(\frac{R}{2}\right)\right)\right)}{\operatorname{vol}(B^{M}(r,x))}$
$\displaystyle\leqslant\frac{\int_{0}^{\sqrt{2}\cdot
2\arcsin\left(\frac{R}{2}\right)}\sinh^{d-1}(t)\ dt}{\int_{0}^{\sqrt{2}\cdot
r}\sinh^{d-1}(t)\ dt},$

and we are done with the upper bound of the second claim. The lower bound is
obtained analogously and we leave the details to the reader.

So far we have been focusing on the case of dimension $d\geqslant 2$. In
dimension $d=1$ we actually get from (3.7) that

$R-r\leqslant\operatorname{vol}(\mathbb{B}^{M}(x,R))-\operatorname{vol}(\mathbb{B}^{M}(x,R))\leqslant
2\arcsin\left(\frac{R-r}{2}\right)$

This implies directly that

$\frac{R}{r}\leqslant\frac{\operatorname{vol}(\mathbb{B}^{M}(x,R))}{\operatorname{vol}(\mathbb{B}^{M}(x,r))}\leqslant
1+\frac{2\arcsin\left(\frac{R-r}{2}\right)}{r}$

for all $0<r<R<1$. The so-obtained bound for $d=1$ is slightly better than the
one we had claimed. We also note that this bound still works for $0<r<R<2$. ∎

###  3.4. Volume of thick diagonal

Our goal here is to estimate how the volume
$\operatorname{vol}(DM(\varepsilon))$ of the $\varepsilon$-thick diagonal

$DM(\varepsilon)=\\{(x,y)\in M\times M\text{ with
}|x-y|\leqslant\varepsilon\\}$

of a reach-1 submanifold $M\subset\mathbb{R}^{s}$ varies when we replace
$\varepsilon$ by something else. Here the volume is computed as a subset of
the Riemannian manifold $M\times M$, but can be expressed as an integral over
$M$ as follows:

(3.8)
$\operatorname{vol}(DM(\varepsilon))=\int_{M}\operatorname{vol}(\mathbb{B}^{M}(x,\varepsilon))\
dx$

In other words, $\operatorname{vol}(DM(\varepsilon))$ is nothing other than
the correlation integral at scale $\varepsilon$ of the riemannian measure of
$M$ when considered as a measure on the ambient euclidean space
$\mathbb{R}^{s}$. We stress that the thick diagonal is defined in terms of the
ambient distance in euclidean space, not in terms of the intrinsic distance of
$M$. Anyways, now we prove the following:

######  Theorem 3.1.

Suppose that $M\subset\mathbb{R}^{s}$ is a $d$-dimensional ($d\geqslant 1$)
closed submanifold with reach $\tau(M)\geqslant 1$. Then we have

$\frac{\frac{\varepsilon_{1}}{2}}{\arcsin(\frac{\varepsilon_{2}}{2})}\\!\left(\\!\frac{\sin(\varepsilon_{1})}{\sin(2\arcsin\frac{\varepsilon_{2}}{2})}\\!\right)^{d-1}\\!\\!\\!\\!\leqslant\frac{\operatorname{vol}(DM(\varepsilon_{1}))}{\operatorname{vol}(DM(\varepsilon_{2}))}\\!\leqslant\frac{\int_{0}^{\sqrt{2}\cdot
2\cdot\arcsin(\frac{\varepsilon_{1}}{2})}\\!\sinh^{d-1}(t)dt}{\int_{0}^{\sqrt{2}\cdot\varepsilon_{2}}\sinh^{d-1}(t)dt}$

for any two $0<r<R<1$.

###### Proof.

From the expression (3.8) we get that

$\min_{x\in
M}\frac{\operatorname{vol}(\mathbb{B}^{M}(R,x))}{\operatorname{vol}(\mathbb{B}^{M}(r,x))}\leqslant\frac{\operatorname{vol}(DM(R))}{\operatorname{vol}(DM(r))}\leqslant\max_{x\in
M}\frac{\operatorname{vol}(\mathbb{B}^{M}(R,x))}{\operatorname{vol}(\mathbb{B}^{M}(r,x))}$

Now the claim follows from Corollary 3.4. ∎

###  3.5. The gap

As we mentioned in the introduction, if we sample more and more points from a
manifold $M$ and we apply the our algorithm, then what we are doing is
computing the quantity

$\frac{\log(\operatorname{vol}(DM(\varepsilon_{1})))-\log(\operatorname{vol}(DM(\varepsilon_{2})))}{\log(\varepsilon_{1})-\log(\varepsilon_{2})}.$

Now, armed with Theorem 3.1 we could analyse what happens when one of the
scales $\varepsilon_{1}$ and $\varepsilon_{2}$ tends to $0$, or when the gap
between them tends to $0$, or when the dimension grows. All of this would be
nice and well, but what we actually care about is to find scales that on the
one hand keep our estimator reliable while working with as few points as
possible. Implementing numerically a procedure described in Section 4.4 below,
we find convenient scales for the dimensions we are mostly interested in. Note
that the condition $R<1$ can be replaced by $R<2$ for dimension 1, as
explained in the proof of Corollary 3.4.

######  Corollary 3.5.

Let $M\subset\mathbb{R}^{s}$ be a submanifold of dimension $d=1,2,\dots,10$
and with reach $\tau(M)\geqslant 1$, and let $\varepsilon_{1},\varepsilon_{2}$
and $\operatorname{gap}_{d}$ be as in the table below. Then we have

$d-\frac{1}{2}+\operatorname{gap}_{d}\leqslant\frac{\log\left(\frac{\operatorname{vol}(\mathbb{B}^{M}(\varepsilon_{1},x))}{\operatorname{vol}(\mathbb{B}^{M}(\varepsilon_{2},x))}\right)}{\log\left(\frac{\varepsilon_{1}}{\varepsilon_{2}}\right)}\leqslant
d+\frac{1}{2}-\operatorname{gap}_{d}$

$d$ | $\varepsilon_{1}$ | $\varepsilon_{2}$ | $\operatorname{gap}_{d}$
---|---|---|---
$1$ | $1.5$ | $0.19$ | $0.463241$
$2$ | $0.78$ | $0.2$ | $0.387573$
$3$ | $0.63$ | $0.23$ | $0.307476$
$4$ | $0.54$ | $0.23$ | $0.249891$
$5$ | $0.46$ | $0.22$ | $0.223958$
$d$ | $\varepsilon_{1}$ | $\varepsilon_{2}$ | $\operatorname{gap}_{d}$
---|---|---|---
$6$ | $0.4$ | $0.21$ | $0.208521$
$7$ | $0.36$ | $0.21$ | $0.178814$
$8$ | $0.33$ | $0.2$ | $0.166892$
$9$ | $0.31$ | $0.19$ | $0.155560$
$10$ | $0.29$ | $0.18$ | $0.152528$

∎

##  4\. Sampling the thick diagonal

In this section we will be still assuming that $M\subset\mathbb{R}^{s}$ is a
closed $d$-dimensional submanifold ($d\geqslant 1$) with reach
$\tau(M)\geqslant 1$. Basically our goal is to bound the number of points that
we have to sample from $M$ to get a decent result when we use (1.1).

###  4.1. Some probability

Suppose that we have a symmetric, say bounded, function

$f:M\times M\to\mathbb{R}$

We are interested in the sequence of random variables

(4.9) $X_{n}^{f}:M^{\mathbb{N}}\to\mathbb{R},\ \
X_{n}(x_{1},x_{2},\dots)\mapsto\sum_{i,j\leqslant n,\ i\neq j}f(x_{j},x_{j})$

when $n$ tends to $\infty$. Here we have endowed $M$ with the probability
measure

$\operatorname{Prob}=\frac{1}{\operatorname{vol}(M)}\operatorname{vol}$

proportional to the Riemannian measure. Accordingly, $M^{k}$ and
$M^{\mathbb{N}}$ are all endowed with the corresponding product measure, again
denoted by $\operatorname{Prob}$.

Being the sum of random variables, the expectation and variance of $X_{n}^{f}$
are easy to get. Here they are:

$\displaystyle E(X_{n}^{f})$ $\displaystyle=n(n-1)\cdot E(f)$
$\displaystyle\operatorname{Var}(X_{n}^{f})$ $\displaystyle=2\cdot
n(n-1)\cdot\operatorname{Var}(f)+4\cdot n(n-1)(n-2)\cdot\operatorname{cov}(f)$

where $\operatorname{cov}(f)$ is the co-variance of
$(x_{1},\dots,x_{n})\mapsto f(x_{1},x_{2})$ and $(x_{1},\dots,x_{n})\mapsto
f(x_{1},x_{3})$, or in a formula

$\operatorname{cov}(f)=\int_{M\times M\times M}f(x,y)\cdot f(x,z)\
d\operatorname{Prob}(x,z,y)-E(f)^{2}.$

Besides the expectation and the variance, what we will need to estimate is the
quantity $\frac{\operatorname{Var}}{E^{2}}$ for the random variables
$X_{n}^{f}$. Well, this is what we get if we just use our expressions for the
expectation and the variance:

(4.10)
$\frac{\operatorname{Var}(X_{n}^{f})}{E(X_{n}^{f})^{2}}=\frac{2}{n(n-1)}\cdot\frac{\operatorname{Var}(f)}{E(f)^{2}}+\frac{4(n-2)}{n(n-1)}\cdot\frac{\operatorname{cov}(f)}{E(f)^{2}}$

The reason why we will care about this last quantity is the following surely
standard consequence of the Bienaymé-Chebyshev inequality:

######  Lemma 4.1.

For any integrable random variable $X$ in a probability space $X$ we have

$\operatorname{Prob}\left(\left|\log\left(\frac{X}{E(X)}\right)\right|>\delta\right)\leqslant\frac{1}{(1-e^{-\delta})^{2}}\cdot\frac{\operatorname{Var}(X)}{E(X)^{2}}$

###### Proof.

Well, let us compute

$\displaystyle\operatorname{Prob}\left(\left|\log\left(\frac{X}{E(X)}\right)\right|>\delta\right)$
$\displaystyle=\operatorname{Prob}\left(\frac{X}{E(X)}\notin[e^{-\delta},e^{\delta}]\right)$
$\displaystyle\leqslant\operatorname{Prob}\left(\left|\frac{X}{E(X)}-1\right|>1-e^{-\delta}\right)$
$\displaystyle\leqslant\operatorname{Prob}\left(\left|X-E(X)\right|>(1-e^{-\delta})\cdot
E(X)\right)$

Setting

$k=\frac{(1-e^{-\delta})\cdot E(X)}{\sqrt{\operatorname{Var}(X)}}$

in the standard Bienaymé-Chebysheff inequality

$\operatorname{Prob}\left(|X-E(X)|\geqslant
k\cdot\sqrt{\operatorname{Var}(X)}\right)\leqslant k^{-2}$

we get

$\operatorname{Prob}\left(\left|\log\left(\frac{X}{E(X)}\right)\right|>\delta\right)\leqslant\frac{\operatorname{Var}(X)}{(1-e^{-\delta})^{2}\cdot
E(X)^{2}}$

as we had claimed. ∎

###  4.2. The function we care about

We are going to be interested in all of this in the case that
$f=f_{\varepsilon}$ is the characteristic function of $DM(\varepsilon)$, that
is

(4.11) $f_{\varepsilon}(x,y)=\left\\{\begin{array}[]{cl}1&\text{ if
}|x-y|\leqslant\varepsilon\\\ 0&\text{ otherwise}\end{array}\right.$

This function satisfies that

$E(f_{\varepsilon})=\frac{\operatorname{vol}(DM(\varepsilon))}{\operatorname{vol}(M)^{2}}$

and hence we get that

$E(X_{n}^{f_{\varepsilon}})=n(n-1)\cdot\frac{\operatorname{vol}(DM(\varepsilon))}{\operatorname{vol}(M)^{2}}$

for all $n\geqslant 2$. Besides knowing the expectation, to be able to use
Lemma 4.1 when we need to know, or at least estimate, is the quantity
$\frac{\operatorname{Var}}{E^{2}}$. To apply (4.10) we need first to be able
to estimate the variance and covariance of $f_{\varepsilon}$. Well, since
$f_{\varepsilon}$ only takes the values $0$ and $1$, the variance is easily
calculated:

$\operatorname{Var}(f_{\varepsilon})=E(f_{\varepsilon})-E(f_{\varepsilon})^{2}=\frac{\operatorname{vol}(DM(\varepsilon))}{\operatorname{vol}(M)^{2}}-\left(\frac{\operatorname{vol}(DM(\varepsilon))}{\operatorname{vol}(M)^{2}}\right)^{2}$

When it comes to the covariance we have

$\displaystyle\operatorname{cov}(f_{\varepsilon})$
$\displaystyle=\int\left(\frac{\operatorname{vol}(\mathbb{B}^{M}(x,\varepsilon))}{\operatorname{vol}(M)}\right)^{2}d\operatorname{Prob}(x)-\frac{\operatorname{vol}(DM(\varepsilon))^{2}}{\operatorname{vol}(M)^{4}}$
$\displaystyle=\int\left(\frac{\operatorname{vol}(\mathbb{B}^{M}(x,\varepsilon))}{\operatorname{vol}(M)}\right)^{2}d\operatorname{Prob}(x)-\left(\int_{M}\frac{\operatorname{vol}(DM(\varepsilon))}{\operatorname{vol}(M)}d\operatorname{Prob}(x)\right)^{2}$
$\displaystyle=\operatorname{Var}\left(x\mapsto\frac{\operatorname{vol}(\mathbb{B}^{M}(x,\varepsilon))}{\operatorname{vol}(M)}\right)$

This means that when the volume of
$\mathbb{B}^{M}(x,\varepsilon)=M\cap\mathbb{B}(x,\varepsilon)$ is constant
then the covariance vanishes. This is for example the case for
$M=\mathbb{S}^{d}\subset\mathbb{R}^{d+1}$ or for the Clifford torus
$M=\mathbb{T}^{d}\subset\mathbb{R}^{2d}$. However, in general we do not get
anything better than the bound coming from Popoviciu’s inequality, that is

$\operatorname{cov}(f_{\varepsilon})\leqslant\frac{\left(V^{M}_{\max}(\varepsilon)-V^{M}_{\min}(\varepsilon)\right)^{2}}{4\cdot\operatorname{vol}(M)^{2}}$

where we have set

$V_{\max}^{M}(\varepsilon)=\max_{x\in
M}\operatorname{vol}(\mathbb{B}^{M}(x,\varepsilon))\text{ and
}V_{\min}^{M}(\varepsilon)=\min_{x\in
M}\operatorname{vol}(\mathbb{B}^{M}(x,\varepsilon))$

Now, using (4.10), the bound for $\operatorname{cov}(f_{\varepsilon})$, as
well as the bound
$\operatorname{vol}(DM(\varepsilon))\geqslant\operatorname{vol}(M)\cdot
V_{\min}(\varepsilon)$ we get that

$\frac{\operatorname{Var}(X_{n}^{f_{\varepsilon}})}{E(X_{n}^{f_{\varepsilon}})^{2}}\leqslant\frac{2}{(n-1)^{2}}\cdot\frac{\operatorname{vol}(M)}{V^{M}_{\min}(\varepsilon)}+\frac{1}{n-1}\cdot\left(\frac{V_{\max}^{M}(\varepsilon)}{V_{\min}^{M}(\varepsilon)}-1\right)^{2}$

To get bounds that only depend on the dimension and on $\varepsilon<2$ recall
that from Corollary 3.4 we get that

(4.12)
$\begin{split}V^{M}_{\min}(\varepsilon)&\geqslant\operatorname{vol}(B^{\mathbb{S}^{d}}(\varepsilon))\stackrel{{\scriptstyle\text{def}}}{{=}}\mathcal{V}(\varepsilon)\\\
\frac{V_{\max}^{M}(\varepsilon)}{V_{\min}^{M}(\varepsilon)}&\leqslant\frac{\operatorname{vol}\left(B^{\frac{1}{\sqrt{2}}\mathbb{H}^{d}}\left(2\arcsin\left(\frac{\varepsilon}{2}\right)\right)\right)}{\operatorname{vol}(B^{\mathbb{S}^{d}}(\varepsilon))}\stackrel{{\scriptstyle\text{def}}}{{=}}\mathcal{R}(\varepsilon)\end{split}$

Using these bounds we get

$\frac{\operatorname{Var}(X_{n}^{f_{\varepsilon}})}{E(X_{n}^{f_{\varepsilon}})^{2}}\leqslant\frac{2}{(n-1)^{2}}\cdot\frac{\operatorname{vol}(M)}{\mathcal{V}(\varepsilon)}+\frac{1}{n-1}\cdot\left(\mathcal{R}(\varepsilon)-1\right)^{2}$

We record what we have so far:

######  Lemma 4.2.

Let $M\subset\mathbb{R}^{s}$ be a closed submanifold with dimension
$\dim(M)=d$ and reach $\tau(M)\geqslant 1$, and for some $\varepsilon<2$ and
$n\in\mathbb{N}$ consider $f_{\varepsilon}$ and $X_{n}^{f_{\varepsilon}}$ as
in (4.11) and (4.9). Then we have

$\displaystyle E(X_{n}^{f_{\varepsilon}})=$ $\displaystyle
n(n-1)\cdot\frac{\operatorname{vol}(DM(\varepsilon))}{\operatorname{vol}(M)^{2}}$
$\displaystyle\frac{\operatorname{Var}(X_{n}^{f_{\varepsilon}})}{E(X_{n}^{f_{\varepsilon}})^{2}}\leqslant$
$\displaystyle\frac{2}{(n-1)^{2}}\cdot\frac{\operatorname{vol}(M)}{\mathcal{V}(\varepsilon)}+\frac{1}{n-1}\cdot\left(\mathcal{R}(\varepsilon)-1\right)^{2}$

where $\mathcal{V}(\varepsilon)$ and $\mathcal{R}(\varepsilon)$ are as in
(4.12). ∎

Before moving any further let us give explicit formulas for
$\mathcal{V}(\varepsilon)$ and $\mathcal{R}(\varepsilon)$:

(4.13)
$\begin{split}\mathcal{V}(\varepsilon)&=\operatorname{vol}(\mathbb{S}^{d-1})\cdot\int_{0}^{\varepsilon}\sin^{d-1}(t)dt\\\
\mathcal{R}(\varepsilon)&=\frac{2^{-\frac{d}{2}}\int_{0}^{2\sqrt{2}\arcsin\frac{\varepsilon}{2}}\sinh(t)^{d-1}dt}{\int_{0}^{\varepsilon}\sin^{d-1}(t)dt}\end{split}$

Note also that the two summands in the bound for
$\frac{\operatorname{Var}(X_{n}^{f_{\varepsilon}})}{E(X_{n}^{f_{\varepsilon}})^{2}}$
in Lemma 4.2 are rather different. Assume for example that $d$ is fixed. Then
the weight of the second factor decreases when $\varepsilon$ decreases. On the
other hand the value of the first one explodes. Recall also that the second
factor can be ignored if $\operatorname{cov}$ vanishes, that is if all balls
$\mathbb{B}^{M}(x,\varepsilon)$ have the same volume.

###  4.3. Some technical results

Recall that to estimate the dimension of $M$ via (1.1), or equivalently via
(2.5), what we do is to take, for two scales $\varepsilon_{1}>\varepsilon_{2}$
random values of $X^{\varepsilon_{1}}_{n}=X^{f_{\varepsilon_{1}}}_{n}$ and
$X^{\varepsilon_{2}}_{n}=X^{f_{\varepsilon_{2}}}_{n}$, compute

$\frac{\log\frac{X^{\varepsilon_{1}}_{n}}{n(n-1)}-\log\frac{X^{\varepsilon_{2}}_{n}}{n(n-1)}}{\log\varepsilon_{1}-\log\varepsilon_{2}}$

and hope the that obtained value has something to do with $\dim(M)$. Well,
what we get from Lemma 4.1 is an estimate of the probability that this value
is far from the expectation. Indeed, since we have

$\displaystyle\operatorname{Prob}\left(\left|\frac{\log
X^{\varepsilon_{1}}_{n}-\log
X^{\varepsilon_{2}}_{n}}{\log\varepsilon_{1}-\log\varepsilon_{2}}-\frac{\log\operatorname{vol}(DM(\varepsilon_{1}))-\log\operatorname{vol}(DM(\varepsilon_{2}))}{\log\varepsilon_{1}-\log\varepsilon_{2}}\right|>\rho\right)$
$\displaystyle\phantom{blablabla}=\operatorname{Prob}\left(\left|\frac{\log
X^{\varepsilon_{1}}_{n}-\log
X^{\varepsilon_{2}}_{n}}{\log\varepsilon_{1}-\log\varepsilon_{2}}-\frac{\log
E(X^{\varepsilon_{1}}_{n})-\log
E(X^{\varepsilon_{2}}_{n})}{\log\varepsilon_{1}-\log\varepsilon_{2}}\right|>\rho\right)$
$\displaystyle\phantom{blablabla}=\operatorname{Prob}\left(\left|\log\left(\frac{X^{\varepsilon_{1}}_{n}}{E(X^{\varepsilon_{1}}_{n})}\right)-\log\left(\frac{X^{\varepsilon_{2}}_{n}}{E(X^{\varepsilon_{2}}_{n})}\right)\right|>\log\left(\left(\frac{\varepsilon_{1}}{\varepsilon_{2}}\right)^{\rho}\right)\right)$
$\displaystyle\phantom{blablabla}\leqslant\sum_{i=1,2}\operatorname{Prob}\left(\left|\log\left(\frac{X^{\varepsilon_{i}}_{n}}{E(X^{\varepsilon_{i}}_{n})}\right)\right|>\log\left(\left(\frac{\varepsilon_{1}}{\varepsilon_{2}}\right)^{\frac{\rho}{2}}\right)\right)$
$\displaystyle\phantom{blablabla}\leqslant\frac{1}{\left(1-\left(\frac{\varepsilon_{2}}{\varepsilon_{1}}\right)^{\frac{\rho}{2}}\right)^{2}}\sum_{i=1,2}\frac{\operatorname{Var}(X_{n}^{\varepsilon_{i}})}{E(X_{n}^{\varepsilon_{i}})^{2}}$

Plugging in the statement of Lemma 4.2 we get:

######  Theorem 4.1.

Let $M\subset\mathbb{R}^{s}$ be a closed submanifold with reach
$\tau(M)\geqslant 1$, pick two scales $0<\varepsilon_{2}<\varepsilon_{1}$.
Also, for $n\geqslant 2$ set

$\rho=\sum_{i=1,2}\left(\frac{2}{(n-1)^{2}}\cdot\frac{\operatorname{vol}(M)}{\mathcal{V}(\varepsilon_{i})}+\frac{1}{n-1}\cdot\left(\mathcal{R}(\varepsilon_{i})-1\right)^{2}\right)$

Then we have

$\operatorname{Prob}\left(\left|\frac{\log X^{\varepsilon_{1}}_{n}-\log
X^{\varepsilon_{2}}_{n}}{\log\varepsilon_{1}-\log\varepsilon_{2}}-\frac{\log\left(\frac{\operatorname{vol}(DM(\varepsilon_{1}))}{\operatorname{vol}(DM(\varepsilon_{2}))}\right)}{\log\varepsilon_{1}-\log\varepsilon_{2}}\right|>\Delta\right)\leqslant\frac{\rho}{\left(1-\left(\frac{\varepsilon_{2}}{\varepsilon_{1}}\right)^{\frac{\Delta}{2}}\right)^{2}}$

for any $\delta>0$. ∎

Let us get an slightly more user friendly version:

######  Corollary 4.3.

With the same assumptions and notation as in Theorem 4.1 suppose that for some
positive $\alpha_{i}$’s with $\alpha_{1}+\alpha_{2}=1$ and for some $\rho>0$
we have

$n\geqslant
1+\frac{1}{\alpha_{i}\cdot\rho}\cdot\left(\mathcal{R}(\varepsilon_{i})-1\right)^{2}+\sqrt{\frac{2}{\alpha_{i}\cdot\rho}\cdot\frac{\operatorname{vol}(M)}{\mathcal{V}(\varepsilon_{i})}}$

for $i=1,2$. Then we also have

$P\left(\left|\frac{\log X^{\varepsilon_{1}}_{n}-\log
X^{\varepsilon_{2}}_{n}}{\log\varepsilon_{1}-\log\varepsilon_{2}}-\frac{\log\left(\frac{\operatorname{vol}(DM(\varepsilon_{1}))}{\operatorname{vol}(DM(\varepsilon_{2}))}\right)}{\log\varepsilon_{1}-\log\varepsilon_{2}}\right|>\Delta\right)\leqslant\rho\cdot\left(1-\left(\frac{\varepsilon_{2}}{\varepsilon_{1}}\right)^{\frac{\Delta}{2}}\right)^{-2}$

###### Proof.

In terms of Theorem 4.1 what we have to do is to guarantee for $i=1,2$ that

$\alpha_{i}\rho\geqslant\frac{2}{(n-1)^{2}}\cdot\frac{\operatorname{vol}(M)}{\mathcal{V}(\varepsilon_{i})}+\frac{1}{n-1}\cdot\left(\mathcal{R}(\varepsilon_{i})-1\right)^{2}$

It thus suffices to ensure that $n-1$ is larger than the solution $X$ of the
equation

$a\stackrel{{\scriptstyle
def}}{{=}}\alpha_{i}\rho=\frac{1}{X^{2}}\cdot\frac{2\cdot\operatorname{vol}(M)}{\mathcal{V}(\varepsilon_{i})}+\frac{1}{X}\cdot\left(\mathcal{R}(\varepsilon_{i})-1\right)^{2}\stackrel{{\scriptstyle
def}}{{=}}\frac{1}{X^{2}}c+\frac{1}{X}b$

This is now a quadratic equation with positive solution

$X=\frac{b}{a}+\sqrt{\frac{c}{a}}$

The claim follows. ∎

Again, if $M$ is such that all balls $\mathbb{B}^{M}(x,\varepsilon)$ have
constant volume, then one can replace the first displayed equation in
Corollary 4.3 by

$n\geqslant
1+\sqrt{\frac{2}{\alpha_{i}\cdot\rho}\cdot\frac{\operatorname{vol}(M)}{\mathcal{V}(\varepsilon_{i})}}.$

###  4.4. Searching decent scales

Given the dimension $d$ and the volume $\operatorname{Vol}(M)$, what are the
optimal scales to run (1.1) so that we have $\geqslant 90\%$ success
probability? Let’s see how we could find, if not the optimal scales, at least
decent ones. First, for $1>\varepsilon_{1}>\varepsilon_{2}>0$ consider the
quantity

$\Delta=\Delta_{\varepsilon_{1},\varepsilon_{2}}=\max\\{\Delta_{1},\Delta_{2}\\}$

where

$\displaystyle\Delta_{1}$
$\displaystyle\leqslant\frac{1}{2}-\frac{\log\left(\frac{\int_{0}^{\sqrt{2}\cdot
2\cdot\arcsin(\frac{\varepsilon_{1}}{2})}\sinh^{d-1}(t)\
dt}{\int_{0}^{\sqrt{2}\cdot{\varepsilon_{2}}}\sinh^{d-1}(t)\
dt}\right)}{\log\frac{\varepsilon_{1}}{\varepsilon_{2}}}+d$
$\displaystyle\Delta_{1}$
$\displaystyle\leqslant\frac{1}{2}+\frac{\log\left(\frac{\varepsilon_{1}}{2\cdot\arcsin(\frac{\varepsilon_{2}}{2})}\cdot\left(\frac{\sin(\varepsilon_{1})}{\sin(2\cdot\arcsin\frac{\varepsilon_{2}}{2})}\right)^{d-1}\right)}{\log\frac{\varepsilon_{1}}{\varepsilon_{2}}}-d$

From Lemma 4.2 and Theorem 3.1 we get for all $n$ that

$\left|\frac{\log\frac{E(X_{n}^{\varepsilon_{1}})}{E(X_{n}^{\varepsilon_{2}})}}{\log\varepsilon_{1}-\log\varepsilon_{2}}-d\right|=\left|\frac{\log\frac{\operatorname{vol}(DM(\varepsilon_{1}))}{\operatorname{vol}(DM(\varepsilon_{2}))}}{\log\varepsilon_{1}-\log\varepsilon_{2}}-d\right|\leqslant\frac{1}{2}-\Delta$

Note that Theorem 3.1 has the condition $R<1$, but for dimension 1, it be
replaced by the condition $R<2$, as explained in the proof of Corollary 3.4.

It follows that, as long as $\Delta>0$, if we take a very large number of
points $n$ then we get that it is very likely that the (1.1) returns the value
$d$. Now, how many points we do actually need if we want to guarantee a $90\%$
rate of success? Well, with notation as in Theorem 4.1 we start by setting

(4.14)
$\rho=\rho_{\varepsilon_{1},\varepsilon_{2}}=\frac{1}{10}\cdot\left(1-\left(\frac{\varepsilon_{2}}{\varepsilon_{1}}\right)^{\frac{\Delta}{2}}\right)^{2}$

Now, once we have $\rho$ we get from Corollary 4.3 that if we take
$\alpha\in(0,1)$, set $\alpha_{1}=\alpha$ and $\alpha_{2}=1-\alpha$, and if we
take at least

$\displaystyle
n(\varepsilon_{1},\varepsilon_{2},\alpha,\operatorname{vol}(M))$
$\displaystyle=\max_{i=1,2}\left(1+\frac{1}{\alpha_{i}\cdot\rho}\cdot\left(\mathcal{R}(\varepsilon_{i})-1\right)^{2}+\sqrt{\frac{2}{\alpha_{i}\cdot\rho}\cdot\frac{\operatorname{vol}(M)}{\mathcal{V}(\varepsilon_{i})}}\right)$
$\displaystyle\leqslant
1+\max_{i=1,2}\left(\frac{1}{\alpha_{i}\cdot\rho}\cdot\left(\mathcal{R}(\varepsilon_{i})-1\right)^{2}\right)+$
$\displaystyle\phantom{BLABLABLA}+\left(\max_{i=1,2}\sqrt{\frac{2}{\alpha_{i}\cdot\rho\cdot\mathcal{V}(\varepsilon_{i})^{2}}}\right)\cdot\operatorname{vol}(M)^{\frac{1}{2}}$

points, then

(4.15)
$\operatorname{Prob}\left(\left|\frac{\log\frac{X_{n}^{{\varepsilon_{1}}}}{X_{n}^{{\varepsilon_{2}}}}}{\log\varepsilon_{1}-\log\varepsilon_{2}}-d\right|<\frac{1}{2}d\right)>90\%.$

If we are interested in manifolds with $\operatorname{Vol}(M)\leqslant V$ then
(4.15) holds as long as we take at least

(4.16) $\min_{\tiny\begin{array}[]{l}1>\varepsilon_{1}>\varepsilon_{2}>0\\\
\text{with
}\Delta_{\varepsilon_{1},\varepsilon_{2}}>0\end{array}}\min_{\alpha\in(0,1)}n(\varepsilon_{1},\varepsilon_{2},\alpha)$

points and we use (1.1) with constants $1>\varepsilon_{1}>\varepsilon_{2}>0$.
Now, to find decent scales we can now minimize (4.16). A program which
numerically approximates that is available at [19].

In fact, running also the program in each $d=1,2,\dots,10$ for the volume of
the corresponding $d$-dimensional torus we get that
$(\varepsilon_{1},\varepsilon_{2},\alpha)$ as in Table 2 give smallish values
for (4.16). If we plug these constants in the formula for
$n(\varepsilon_{1},\varepsilon_{2},\alpha,\operatorname{vol}(M))$ that we gave
above we recover the statement of Theorem 1.1 stated in the introduction:

$d$ | $\varepsilon_{1}$ | $\varepsilon_{2}$ | $\alpha_{1}$ |  | $d$ | $\varepsilon_{1}$ | $\varepsilon_{2}$ | $\alpha_{1}$
---|---|---|---|---|---|---|---|---
1 | $1.5$ | $0.19$ | $0.15$ |  | 6 | $0.4$ | $0.21$ | $0.03$
2 | $0.78$ | $0.2$ | $0.11$ |  | 7 | $0.36$ | $0.21$ | $0.03$
3 | $0.63$ | $0.23$ | $0.09$ |  | 8 | $0.33$ | $0.2$ | $0.02$
4 | $0.54$ | $0.23$ | $0.06$ |  | 9 | $0.31$ | $0.19$ | $0.02$
5 | $0.46$ | $0.22$ | $0.04$ |  | 10 | $0.29$ | $0.18$ | $0.01$
Table 2. Decent scales for
$\operatorname{vol}=\operatorname{vol}(\mathbb{T}^{d})$ in dimension
$d=1,2,\dots,10$. It is evident that the values in Table 2 can change if
instead of using $\operatorname{vol}(\mathbb{T}^{d})$ as an input one chooses
any other value. However, for whatever it is worth, if instead one chooses
$10\cdot\operatorname{vol}(\mathbb{T}^{d})$ or even
$100\cdot\operatorname{vol}(\mathbb{T}^{d})$ then nor much changes: in small
dimensions (that is, up to dimension 3) the scales increase a bit, but for
dimensions at least $4$ nothing changes.

######  Theorem 1.1.

For $d=1,\cdots,10$ let $\varepsilon_{1}$ and $\varepsilon_{2}$ be scales as
in the table below. Also, given a closed $d$-dimensional manifold
$M\subset\mathbb{R}^{s}$ with reach $\tau(M)\geqslant 1$ let $n$ be also as in
the following table:

$d$ | $\varepsilon_{1}$ | $\varepsilon_{2}$ | n
---|---|---|---
$1$ | $1.5$ | $0.19$ | $9+21\cdot\operatorname{vol}(M)^{\frac{1}{2}}$
$2$ | $0.78$ | $0.2$ | $94+58\cdot\operatorname{vol}(M)^{\frac{1}{2}}$
$3$ | $0.63$ | $0.23$ | $635+146\cdot\operatorname{vol}(M)^{\frac{1}{2}}$
$4$ | $0.54$ | $0.23$ | $2786+392\cdot\operatorname{vol}(M)^{\frac{1}{2}}$
$5$ | $0.46$ | $0.22$ | $7013+1119\cdot\operatorname{vol}(M)^{\frac{1}{2}}$
$6$ | $0.4$ | $0.21$ | $13221+3366\cdot\operatorname{vol}(M)^{\frac{1}{2}}$
$7$ | $0.36$ | $0.21$ | $25138+10644\cdot\operatorname{vol}(M)^{\frac{1}{2}}$
$8$ | $0.33$ | $0.2$ | $50033+34890\cdot\operatorname{vol}(M)^{\frac{1}{2}}$
$d$ | $\varepsilon_{1}$ | $\varepsilon_{2}$ | n
---|---|---|---
$9$ | $0.31$ | $0.19$ | $63876+119533\cdot\operatorname{vol}(M)^{\frac{1}{2}}$
$10$ | $0.29$ | $0.18$ | $139412+425554\cdot\operatorname{vol}(M)^{\frac{1}{2}}$

Then, if we sample independently and according to the riemannian volume form a
subset $X\subset M$ consisting of at least $n$ points, then we have

$\dim_{\operatorname{Corr}(\varepsilon_{1},\varepsilon_{2})}(X)=d$

with at least $90\%$ probability.∎

In the next section we discuss some (much smaller) heuristic bounds, discuss
some numerical experiments, and compare with the performance of other
estimators.

##  5\. Heuristics

Theorem 4.1 gives us a bound for the number of points needed in a data set to
be able to get from (1.1) at least $90\%$ of the time its dimension. In
concrete examples, we expect that this confidence level can be achieved with
significantly less points. We will discuss this difference between theory and
practice, suggesting a simpler heuristic model supported by computational
examples.

### Heuristic bound

We begin by discussing a heuristic model representing an ideal situation
without curvature. More concretely we will be running (1.1) at some scales
$\varepsilon_{1}$ and $\varepsilon_{2}$ at which we can ignore curvature
effects. In effect, we will act as if all balls in $M$ of radius at most
$\varepsilon_{1}$ were euclidean and totally geodesic.

The first observation is that the statistic (1.1) is computed from information
extracted from the distances $|x-y|$ for those (unordered) pairs

$\\{x,y\\}\in PX(\varepsilon)=\big{\\{}\\{x,y\\}\subset X\text{ with }x\neq
y\text{ with }|x-y|\leqslant\varepsilon\big{\\}}$

rather than from the points themselves. This means that the performance of the
algorithm should depend on $|PX(\varepsilon_{1})|$ and on the dimension $d$,
instead of directly on the total number $n$ of points.

We think of the distance $|x-y|$ for $\\{x,y\\}\in PX(\varepsilon_{1})$ as a
random variable, and from now on, we will take the point of view that we have
$N=|PX(\varepsilon_{1})|$ random variables $(X_{i})_{1\leqslant i\leqslant N}$
given by taking the distance between pairs of points at distance at most
$\varepsilon_{1}$ and sampled uniformly on $M$. We then consider the variables
$Y_{i}$ equal to 1 if $X_{i}$ is smaller than $\varepsilon_{2}$ and 0
otherwise. With this notation in place, the estimator (1.1) becomes

$\dim_{\operatorname{Corr}(\varepsilon_{1},\varepsilon_{2})(X)}=\operatorname{Round}\left(\frac{\log(\frac{1}{N}\sum_{i=1}^{N}Y_{i})}{\log(\varepsilon_{2})-\log(\varepsilon_{1})}\right)$

Note that since we are assuming that all the balls are euclidean, the mean
value of $Y_{i}$ is then $E_{d}=(\varepsilon_{2}/\varepsilon_{1})^{d}$ and as
the variables $Y_{i}$ only take the values 0 and 1, their variance is
$\sigma_{d}^{2}=E_{d}-E_{d}^{2}$.

In reality, the variables $X_{i}$, and thus the variables $Y_{i}$ have no
reason to be independent. Still, most of them are when the volume is large
when compared to the size of the data set. So, from now on, we will put
ourselves in the ideal situation that the $N$ variables $Y_{i}$ are
independent. Independence implies that the distribution of the sample mean
$Z=\frac{1}{N}\sum_{i}Y_{i}$ is binomial of parameters $N$ and $E_{d}$, which
can be approximated by a normal distribution of mean value $E_{d}$ and of
variance $\frac{1}{N}\sigma_{d}^{2}$ using the central limit theorem. It is
then known that the probability of $Z$ to be in the interval
$[E_{d}-1.64\cdot\sigma_{d}/\sqrt{N},E_{d}+1.64\cdot\sigma_{d}/\sqrt{N}]$ is
about 90%. If we set
$\operatorname{gap}_{d}=\min\left(E_{d-0.5}-E_{d},E_{d}-E_{d+0.5}\right)$, we
want to find $N$ so that
$1.64\cdot\sigma_{d}/\sqrt{N}=\operatorname{gap}_{d}$. This number gives a
number of pairs sufficient to obtain the right dimension with a confidence of
90%.

Suppose for example that the manifold $M$ has dimension 4. For
$\varepsilon_{1}=0.54$ and $\varepsilon_{2}=0.23$, the scales coming from
Theorem 1.1, we can then compute the values of $\operatorname{gap}_{4}$ and of
$\sigma_{4}$ and deduce the required value of $N$.

$\displaystyle\operatorname{gap}_{4}$
$\displaystyle=\min(E_{3.5}-E_{4},E_{4}-E_{4.5})\simeq 0.01143$
$\displaystyle\sigma_{4}$ $\displaystyle=\sqrt{E_{4}-E_{4}^{2}}\simeq 0.1784$
$\displaystyle N$
$\displaystyle=\left(1.64\cdot\sigma_{4}/\operatorname{gap}_{4}\right)^{2}\simeq
655\text{ pairs}$

This is an approximation using the central limit theorem (that is, replacing
the binomial distribution by the normal distribution), but more precise
computation can be done working directly with the binomial distribution. Doing
this, we can take $N$ down to 516 (see [19]). This reasoning can be applied to
obtain the required number of pairs for each dimension—the results are
summerized in Table 3.

$d$ | $\varepsilon_{1}$ | $\varepsilon_{2}$ |  $N$ for 90% |  $N$ for 70%
---|---|---|---|---
1 | $1.5$ | $0.19$ | $30$ | $10$
2 | $0.78$ | $0.2$ | $122$ | $40$
3 | $0.63$ | $0.23$ | $249$ | $111$
4 | $0.54$ | $0.23$ | $516$ | $238$
5 | $0.46$ | $0.22$ | $878$ | $360$
6 | $0.4$ | $0.21$ | $1329$ | $554$
7 | $0.36$ | $0.21$ | $1719$ | $698$
8 | $0.33$ | $0.2$ | $2481$ | $1070$
9 | $0.31$ | $0.19$ | $3900$ | $1604$
10 | $0.29$ | $0.18$ | $5849$ | $2414$
Table 3. Heuristic bounds for the size of $PX(\varepsilon_{1})$ needed to have
$90\%$ and $70\%$ rate of success when applying (1.1) to data sets sampled
from a reach $1$ manifold.

We compare next these bounds with those in Theorem 1.1 and then discuss a few
numerical experiments.

### Comparisson between theoretical and heuristic bounds

A problem when comparing the heurestic bounds in Table 3 and those in Theorem
1.1 is that the former ones are given interms of the cardinality of
$PX(\varepsilon_{1})$ while the latter ones are given in terms of the
cardinality of $X$.

To connect these two quantities recall that we are acting as if all
$\varepsilon_{1}$-balls in $M$ were euclidean. Now, we get for example from
Lemma 4.2 that the number $|PX(\varepsilon_{1})|$ of unordered pairs of points
at distance at most $\varepsilon_{1}$ is approximately given by the formula

(5.17)
$|PX(\varepsilon_{1})|\simeq\dfrac{n(n-1)}{2}\dfrac{\operatorname{vol}(B^{\mathbb{R}^{d}}(\varepsilon_{1}))}{\operatorname{vol}(M)}.$

Acting as if (5.17) were to give a perfect relation between the number of
points and that of pairs, we get that to get the estimated $516$ pairs in the
case that $M=\mathbb{T}^{4}$ is the 4-dimensional torus we need $1958$ data
points. In comparison, Theorem 1.1 gives an upper bound of $18262$ points.

Arguing like this we can convert the heuristic bounds in Table 3 to bounds for
the needed cardinality of a data set in terms of the dimension and the volume
of the underlying manifold—see Table 4.

$d$ | heuristic $n$ |  $n$ from Theorem 1.1
---|---|---
1 | $5\cdot\operatorname{vol}(M)^{\frac{1}{2}}$ | $9+21\cdot\operatorname{vol}(M)^{\frac{1}{2}}$
2 | $12\cdot\operatorname{vol}(M)^{\frac{1}{2}}$ | $94+58\cdot\operatorname{vol}(M)^{\frac{1}{2}}$
3 | $22\cdot\operatorname{vol}(M)^{\frac{1}{2}}$ | $635+146\cdot\operatorname{vol}(M)^{\frac{1}{2}}$
4 | $50\cdot\operatorname{vol}(M)^{\frac{1}{2}}$ | $2786+392\cdot\operatorname{vol}(M)^{\frac{1}{2}}$
5 | $128\cdot\operatorname{vol}(M)^{\frac{1}{2}}$ | $7013+1119\cdot\operatorname{vol}(M)^{\frac{1}{2}}$
6 | $355\cdot\operatorname{vol}(M)^{\frac{1}{2}}$ | $13221+3366\cdot\operatorname{vol}(M)^{\frac{1}{2}}$
7 | $964\cdot\operatorname{vol}(M)^{\frac{1}{2}}$ | $25138+10644\cdot\operatorname{vol}(M)^{\frac{1}{2}}$
8 | $2949\cdot\operatorname{vol}(M)^{\frac{1}{2}}$ | $50033+34890\cdot\operatorname{vol}(M)^{\frac{1}{2}}$
9 | $9458\cdot\operatorname{vol}(M)^{\frac{1}{2}}$ | $63876+119533\cdot\operatorname{vol}(M)^{\frac{1}{2}}$
10 | $33021\cdot\operatorname{vol}(M)^{\frac{1}{2}}$ | $139412+425554\cdot\operatorname{vol}(M)^{\frac{1}{2}}$
Table 4. Comparisson between the heuristic bound and the bound in Theorem 1.1
for the number of points that suffice to have $90\%$ success rate when
applying (1.1) to data sets sampled from a reach $1$ manifold.

One should keep in mind that entries the middle column in Table 4 are only
meningful for $\operatorname{vol}(M)$ large. Still, there is a very clear
difference between both bounds, the heuristic bound and that from Theorem 1.1.
This difference is at least in part due to the curvature of the submanifold
$M$, but things are not helped by either all the nested inequalities leading
to Theorem 1.1 or the fact that the Bienaymé-Chebyshev inequality is not very
precise.

### Numerical evidence

Experimental examples seem to confirm the heuristical bounds presented on
Table 3. By numerically sampling points on different manifolds, we obtain
results close to what was expected. Here is the procedure we followed:

  1. (1)

Choose a manifold $M$ of known dimension $d$ and consider the scales and
number of pairs $N$ given by Table 3.

  2. (2)

Uniformly and independently sample points in $M$ until we obtain $N$ pairs at
distance $\varepsilon_{1}$. The sampling is done by repeating a program
specific to the desired manifold that samples a single point randomly and
uniformly on it.

  3. (3)

Estimate the dimension using estimator (1.1) with scales $\varepsilon_{1}$ and
$\varepsilon_{2}$.

  4. (4)

Repeat steps 2 and 3 one hundred times and count the number of success.

We ran this experiment a variety of manifolds of reach 1 and always obtained a
rate of success in a range of $\pm 6\%$ of the target rate. This difference
between the actual rate and the target would be totally normal even in an
ideal situation. Indeed, repeating an experiment 100 time with a probability
of success of 90% gives a result with a standard deviation of 3 successes.
Every experimental result then falls into the usual range of two times the
standard deviation.

The manifolds were chosen to observe different situations:

  * •

Worms: $0$-level set of a randomly produced function on $\mathbb{R}^{2}$. The
precise algorithm used for generating these manifolds and sampling from them
is available in [19], as well as the algorithms for the other manifolds.

  * •

Rotation torus: Rotate around the $z$-axis the circle in the $xy$-plane of
radius 1 and center $(2,0,0)$. This surface has reach $1$ and there is a mix
of positive and negative curvature.

  * •

Clifford torus: The product
$\mathbb{T}^{d}=\mathbb{S}^{1}\times\dots\times\mathbb{S}^{1}\subset\mathbb{R}^{2d}$
of $d$-circles of radius $1$. The Clifford torus is curved in Euclidean space
but is flat as a Riemannian manifold.

  * •

Flat torus: This is an ideal situation. We consider namely the abstract
manifold $\mathbb{R}^{d}/2\pi\cdot\mathbb{Z}^{d}$ with its inner distance—it
is not embedded in some larger Euclidean space.

  * •

Swiss roll: The Swiss roll is one of the standard objects on which manifold
learning argorithms seem to be tested, but it also adds a manifold with
boundary to our list. For the sampling, we used the function make_swiss_roll
from the library scikit-learn.

  * •

Schwarz $P$ surface: This is the triply periodic surface in $\mathbb{R}^{3}$
with equation $\cos(x)+\cos(y)+\cos(z)=0$—it approximates one of Schwartz’s
triply periodic minimal surfaces and thus shows features of negative
curvature, both intrinsic and extrinsic. To be able to deal with a finite
volume surface we consider it as a submanifold of the 3-dimension flat torus.

  * •

Spheres: This is the standard sphere
$\mathbb{S}^{d}=\\{x\in\mathbb{R}^{d+1}\text{ with }\|x\|=1\\}$. Spheres have
reach $1$ and are positively curved.

  * •

Gaussian distribution: The standard Gaussian distribution in $\mathbb{R}^{d}$.
The reason why we test this in particular is to include a non-uniform
distribution in our list.

The way we sample points depends on the concrete manifold under consideration,
but we stress that we are sampling each point independently. More precisely,
we are not aiming at getting point at some uniform distance of each other.

Figure 2. Examples of samplings from manifolds. 500 points on a ”Worms” manifold (Left) and 2 000 points on a piece of the Schwarz surface (Right). Manifold | $d$ | 90% target | 70% target
---|---|---|---
Worms | 1 | 88% | 66%
Rotation torus | 2 | 92% | 70%
Clifford torus | 2 | 89% | 69%
Flat torus | 2 | 88% | 66%
Swiss roll | 2 | 93% | 69%
Schwarz surface | 2 | 88% | 66 %
3-sphere | 3 | 92% | 76%
4-sphere | 4 | 89% | 75%
Product of two rotation tori | 4 | 92% | 70%
Clifford torus | 4 | 93% | 72%
Flat torus | 4 | 90% | 74%
Product of two Schwarz surfaces | 4 | 92% | 72%
Gaussian distribution in $\mathbb{R}^{4}$ | 4 | 90% | 76%
5-sphere | 5 | 93% | 74%
Table 5. Experimental rates of success for different manifolds by considering
the ammount of pairs of points suggested by Table 3.

These results tend to confirm the conclusions of the heuristic model. We can
test this model a bit more by testing the estimator with the number of points
given by Table 4. This is what we did in Table 6, and we still obtain results
close to what was expected.

Manifold | $d$ | rate of success
---|---|---
Clifford torus | 2 | 91%
3-sphere | 3 | 91%
Flat torus | 4 | 91%
Product of two tori | 4 | 94%
Table 6. Experimental rates of success for different manifolds by sampling the
number of points given by the ”heuristic” column in Table 4.

### Reach free estimator

A problem with all the results we have been discussing so far is that in
practice one has little clue what the reach of the underlying manifold could
be. However, as we already mentioned in the introduction, one can actually
derive from Table 3 an estimator which does not need any a priori bound on the
reach.

> Assumption: We have a data set $X\subset\mathbb{R}^{s}$ of which we think
> that it has been sampled from some mysterious submanifold
> $M\subset\mathbb{R}^{s}$. We trust however that our data set is good enough
> and we want to test if $M$ could plaussibly have dimension $d$.

> Test: For the chosen $d$, let $\varepsilon_{1}$, $\varepsilon_{2}$ and $N$
> be as in Table 3 (say from the 90% column). Now take $R>0$ to be minimal
> with $|PX(R)|\geqslant N$ and set
> $r=\frac{\varepsilon_{2}}{\varepsilon_{1}}R$. Now check if
> $\dim_{\operatorname{Corr}(R,r)}(X)=d$.

The heuristic discussion above, as well as the numerical experiments, suggest
that as long as our data set is rich enough so that $R\leqslant
0.54\cdot\tau(M)$ then we should get $\dim_{\operatorname{Corr}(R,r)}(X)=d$
with a 90% probability.

Wondering what would happen if we run this algorithm on more ”real” data sets,
we chose 3 different data sets of respective dimension 1, 2 and 3 (see [19]).
Each data sets consist in 200 grayscale pictures of the 3D-model Suzanne with
random rotations (see Figure 1). We respectively randomized 1, 2 and 3 Euler
angles to obtain the desired dimensions. The pictures are 64 by 64 pixels
large and can thus be represented as points in $\mathbb{R}^{4096}$. The
results are presented in Table 7. For each data set, we test dimension 1, 2, 3
and 4 and we show the estimated dimension (before rounding it to the closest
integer).

Note that when using the parameters of dimension 1 for estimating the
dimension of a higher dimensional set, we usually get no points at distance
$\varepsilon/2$, as $(1.5/0.19)^{2}$ is higher than $30$. When this happens,
the estimator (1.1) cannot be computed, we can only conclude that the
dimension is probably bigger than 1.

Hypothesis | Data set 1 | Data set 2 | Data set 3
---|---|---|---
dimension = 1 | 1.13 | $>1$ | $>1$
dimension = 2 | 1.59 | 2.02 | 3.54
dimension = 3 | 1.64 | 2.03 | 3.33
dimension = 4 | 1.45 | 2.08 | 3.40
Table 7. Testing the dimension of the 3 ”real” data sets. To test dimension=1
(resp. dimension=2, resp. dimension=3) we set $\varepsilon_{1}$ so that
$|PX(\varepsilon_{1})|=30$ (resp. 122, resp. 249). Cells with a pink
background represent the tests in which the estimated dimension is consistent
with the tested dimension. The cells with the value ”$>d$” mean that there
were no pairs at distance $\varepsilon_{2}$ for the corresponding number of
pairs at distance $\varepsilon_{1}$.

We can see in Table 7 that there were no Type I errors (meaning that we never
rejected the true hypothesis). On the other hand we found a Type II error,
when the data set of dimension 1 passed the test for dimension 2.

##  6\. Comparison with other estimators

In this paper we study the estimator (1.1), but a variety of other algorithms
exist. For example, instead of doing statistics using only the distances
between points, we could also use the angles between points, or other more
complex features. These different approaches can be compared using the
previous heuristic model. We refer to [9, 10] and specially to [20, Chapter 3]
for a review of different dimension estimators. We will compare (1.1) with the
estimators ANOVA, local PCA, and with the implementation of (1.1) where one
tries to avoid picking scales, reading instead the dimension from a log-log
chart.

### ANOVA

Diaz, Quiroz and Velasco propose in [12] a method based on the angles refered
as ANOVA in the literature. Their idea is to estimate the local dimension
around a point $x$ by considering its $k$ nearest neighbors and the ${k\choose
2}$ angles at $x$ formed by these points. From the variance of these
${k\choose 2}$ angles, we can identify the closest $\beta_{d}$ and deduce the
dimension $d$, where $\beta_{d}$ is defined as follows.

$\beta_{d}=\frac{1}{\operatorname{vol}(\mathbb{S}^{d-1})^{2}}\int_{\mathbb{S}^{d-1}\times\mathbb{S}^{d-1}}\left(\angle(\theta,\eta)-\dfrac{\pi}{2}\right)^{2}d\theta
d\eta$

The global dimension can then be recovered by taking the median, the mode or
the mean of the local dimensions.

To be able to compare ANOVA to the estimator (1.1) we will instead take the
variance of all angles. More precisely, consider every ordered triple of
points in which all three points are at distance at most $\varepsilon_{1}$
from each other. Each such triple $(x,y,z)$ determines an angle
$\angle(y-x,z-x)$. We compute the variance of the so-obtained angles we can
locate the closest $\beta_{d}$. We take that $d$ to be the ANOVA dimension of
our data set.

Putting ourselves again in the ideal situation that we are working in a scale
at which curvature can be ignored and assuming (and that is a lot of assuming
in this case) that the angles we find are independent of each other, we can
argue as earlier in the discussion of the heuristic bound and we get that, in
dimension 4, we would need at least $652$ angles to achieve a confidence level
of 90% with the ANOVA estimator.

Suppose now that we sample 1958 points from the 4-dimensional Clifford torus
$M=\mathbb{T}^{4}$. That number was chosen so that we expect to have $516$
pairs within $\varepsilon_{1}=0.54$ of each other, the heuristic bound for
(1.1) in dimension $4$. On the other hand we expect to have

$\text{number of unordered triples}={1958\choose
3}\cdot\left(\frac{\operatorname{vol}(B^{\mathbb{R}^{4}}(\varepsilon_{1}))}{\operatorname{vol}(M)}\right)^{2}\simeq
91.$

Each unordered pair gives 3 angles, meaning can expect to find about $273$
angles. This is much less than what we estimated that would needed.

And this phenomenon gets worse when the volume of the underlying manifold
grows: if the manifold has volume $64\pi^{4}$ and if we have a data set with
$516$ pairs of points within $\varepsilon_{1}=0.54$ then we expect to only
have $45$ triples, that is about $135$ angles. The reason for this is that the
number of points needed to have a given number of pairs grows with
$\operatorname{vol}(M)^{\frac{1}{2}}$ while it grows as
$\operatorname{vol}(M)^{\frac{2}{3}}$ when we fix the number of triples
instead. This means that for large volumes the algorithm will only become
worse as finding triples of points will become more and more difficult.

We numerically compared the rates of success of estimator (1.1) and of ANOVA
for the Clifford torus with the same scales. The results are presented in
Table 8. This experiment shows that, for this example and for these scales,
estimator (1.1) gives significantly better results than ANOVA. However, we
cannot conclude that estimator (1.1) has better performance in general. In
particular, for manifolds of smaller volume or for different choices of
scales, ANOVA could in principle give better results.

Manifold | $d$ | number of points | estimator (1.1) | ANOVA
---|---|---|---|---
Clifford torus | 2 | 76 | 93% | 65%
Clifford torus | 3 | 347 | 93% | 67%
Table 8. Comparison of the experimental rates of success between estimator
(1.1) and ANOVA on the Clifford torus.

######  Remark.

Recall that we have considered a variation of ANOVA—our conclusions should
also apply to the original algorithm as well.

### Local PCA

Principal Component Analysis aims to find the best linear space containing a
given data set. According to [7], PCA is the gold standard of dimension
estimation. It works as follows. To our given data set
$X=(x_{1},\dots,x_{n})\subset\mathbb{R}^{N}$ we associate first the mean

$\bar{x}=\frac{1}{n}\sum_{i=1}^{n}x_{n}$

and then the $n\times N$ matrix $A$ whose rows are the vectors
$u_{i}-\bar{u}$, and one computes then the singular values $s_{1}\geqslant
s_{2}\geqslant\dots\geqslant s_{\min\\{n,N\\}}$. If $X$ is contained in a
linear subspace of dimension $d$ then $s_{k}=0$ for all $k\geqslant d+1$.
Accordingly, one can declare that $X$ has PCA-dimension $k$ if the gap
$s_{k}-s_{k+1}$ is maximal. Another possibility would be to fix a threshold
$\varepsilon$ and declare the PCA-dimension of $X$ to be the largest $k$ with
$s_{k}\geqslant\varepsilon$. For example, if what one wants to do is to test
the hypothesis that $X$ is contained in a $d$-dimensional subspace then one
can check for example if $s_{d+1}$ is below some threshold
$\varepsilon<\frac{1}{d+2}$—the number $\frac{1}{d+2}$ arises because it is
the expected value for $s_{d}$ if $X$ is uniformly sampled out of the unit
ball in $\mathbb{R}^{d}$.

In any case, that was PCA. The idea of local PCA, or Nonlinear PCA, is to
apply PCA to certain subsets of the data set and then, for good measure,
average the so obtained numbers. As we see, one does not only need to agree on
what one calls PCA, but also on what subsets does one wants to subject to the
PCA treatement. For example, as in [7] one can cluster the data set using
single linkage clustering111In single linkage clustering the clusters are, for
some $\varepsilon$, the connected components of the graph with vertex set $X$
and where two vertices are joined by an edge if they are with $\varepsilon$ of
each other. and then apply PCA to each custer. This is, in our humble and
uneducated opinion, a very reasonable choice if, for data sets $X$ uniformly
sampled out of a submanifold $M\subset\mathbb{R}^{N}$, what one wants to do is
to find linear spaces (approximately) containing each connected component of
$M$. On the other hand, if what one wants to do is to recover the dimension of
$M$ then it seems reasonable to rather apply PCA to (some of) the sets
$X\cap\mathbb{B}(x,\varepsilon)$ for $x\in X$ and for some $\varepsilon$
chosen so that $M\cap\mathbb{B}(x,\varepsilon)$ is well-approximated by its
tangent space.

Note now that every set consisting of $d+1$ points is contained in a
$d$-dimensional affine subspace of $\mathbb{R}^{s}$. This means that if we
want to distinguish dimension $d$ from dimension $d+1$ we need at the very
least to find $d+2$ tuples of nearby points. Now, if we once again sample
$1958$ points out of the Clifford torus $M=\mathbb{T}^{4}$ so that we can
expect the magic number of $516$ paris of points at most at distance $0.54$
from each other, then we can expect to find

$\begin{array}[]{c}\text{number of unordered}\\\ \text{6-tuples at scale
}0.54\end{array}={1958\choose
6}\cdot\left(\frac{\operatorname{vol}(B^{\mathbb{R}^{4}}(0.54))}{\operatorname{vol}(M)}\right)^{5}\simeq
0.110373....$

In other words, the expectation is to not have any $6$-tuples, meaning that,
if we are working at the same scales as we were implementing (1.1), then we
have no chance of distinguishing our 4-manifold from a 5-manifold using local
PCA. Note that even if we work at a scale of $2$, scale at which we are acting
as if the round sphere were totally flat, then we can expect to only have 2
6-tuples

$\begin{array}[]{c}\text{number of unordered}\\\ \text{6-tuples at scale
}2\end{array}={1958\choose
6}\cdot\left(\frac{\operatorname{vol}(B^{\mathbb{R}^{4}}(2))}{\operatorname{vol}(M)}\right)^{5}\simeq
2.043944....$

As we see, even if we work at scales at which we are flat earthers, we do not
have by far enough sufficiently populated clusters for local PCA to be
meaningful.

### log-log plots

Another method to derive the estimator (1.1) when the reach is unknown, is to
take a lot of scales

$\varepsilon_{1}<\varepsilon_{2}<\dots<\varepsilon_{k},$

and to consider the graph of the piecewise linear function with corners

$\big{(}\log(\varepsilon_{i}),\log(|PX(\varepsilon_{i})|)\big{)}$

and try to read out some sort of meaningful slope of that function (see [20,
Chapter 3]). We implemented this procedure for 1 000 points uniformly sampled
from the product of two 2-dimensional rotation tori of reach 1. We obtained
the graph presented on Figure 3.

$\displaystyle{-2}$$\displaystyle{-1}$$\displaystyle{0}$$\displaystyle{1}$$\displaystyle{2}$$\displaystyle\log(\varepsilon_{i})$$\displaystyle{-5}$$\displaystyle{0}$$\displaystyle{5}$$\displaystyle{10}$$\displaystyle{15}$$\displaystyle{20}$$\displaystyle\log(|PX(\varepsilon_{i})|)$

Figure 3. Log-log plot for a set of 1000 points sampled out of the product of
two rotation tori. The graph also features three lines of slope 3, 4 and 5,
for comparison.

Three different parts can be observed on this graph:

  1. (1)

A flat part, when the scale is smaller than the minimal distance between the
points (represented with the value -1 on the y-axis).

  2. (2)

A mostly linear part whose slope should approximate the dimension of $M$.

  3. (3)

A flat part when the scale becomes greater than the diameter of $M$ (which
plateaus at $\log(1000\cdot 999/2)\simeq 13.12$).

It is notable that in this example the second part looks linear even after
twice the reach, even if the behavior of $\log(|PX(\varepsilon_{i})|)$ beyond
this scale is unpredictable. Around the middle of this graph, the slope of the
linear part seems to be close to 4, as predicted.

We also applied this method from the 3 ”real” datasets of pictures of Suzanne
with random rotations. The results are presented in Figure 4. For dimension 1
and 2, it worked surprisingly well and the log-log plots clearly show what is
the right dimension. The plot for the data set of dimension 3 is not as clear,
but the result seems consistent.

Let us conclude with two further comments on the log-log procedure:

(1) The log-log procedure seems to allow us to estimate the dimension of $M$
without any assumption on its reach. It is however difficult to define a
precise algorithm for it. Consequently, we cannot really compare its
performance with that of the other estimators.

(2) Additionally, this method gives us insight on the choice of scales that
could be relevant. Indeed, the scales determined in Section 4.4 have no reason
to be optimal for the heuristic model. And as this model does not take the
curvature nor the reach into account, it is impossible to use it to find
optimal scales. When considering Figure 3, it seems natural to choose two
scales that are far apart to measure the slope as precisely as possible, but
close enough so that they do no get too close to the extremities. Here, two
scales that would seem relevant could be $\exp(0.5)\simeq 1.6$ and
$\exp(-0.5)\simeq 0.6$. In any case, the choice of scales
$0.54\simeq\exp(-0.6)$ and $0.23\simeq\exp(-1.5)$ seem to be far from optimal.
Using the estimator (1.1) with scales 1.6 and 0.6 on computational examples,
we obtain a rate of success of 99% with only 1 000 points for the product of
two tori. For the scales 0.54 and 0.23, the heuristic model gives us a number
of points of 3916 for a rate of success of 90%. This rate falls to 41% with
only 1 000 points on computational examples. However, these scales have been
chosen after the study of the results to obtain the desired answer on a given
manifold. But this does not give us a method to choose optimal scales a
priori. For example, the scales 1.6 and 0.6 one give a rate of success of 45%
for 100 points on the 4-sphere, while the scale 1 and 0.6 give a rate of 92%.

$\displaystyle{-6}$$\displaystyle{-4}$$\displaystyle{-2}$$\displaystyle{0}$$\displaystyle{2}$$\displaystyle{4}$$\displaystyle\log(\varepsilon_{i})$$\displaystyle{-5}$$\displaystyle{0}$$\displaystyle{5}$$\displaystyle{10}$$\displaystyle{15}$$\displaystyle\log(|PX(\varepsilon_{i})|)$

$\displaystyle{-3}$$\displaystyle{-2}$$\displaystyle{-1}$$\displaystyle{0}$$\displaystyle{1}$$\displaystyle{2}$$\displaystyle{3}$$\displaystyle{4}$$\displaystyle\log(\varepsilon_{i})$$\displaystyle{-5}$$\displaystyle{0}$$\displaystyle{5}$$\displaystyle{10}$$\displaystyle{15}$$\displaystyle\log(|PX(\varepsilon_{i})|)$

$\displaystyle{-1}$$\displaystyle{0}$$\displaystyle{1}$$\displaystyle{2}$$\displaystyle{3}$$\displaystyle\log(\varepsilon_{i})$$\displaystyle{-5}$$\displaystyle{0}$$\displaystyle{5}$$\displaystyle{10}$$\displaystyle{15}$$\displaystyle\log(|PX(\varepsilon_{i})|)$

Figure 4. Log-log plots the 3 ”real” datasets of pictures of Suzanne with
random rotations of dimension 1, 2 and 3. The blue lines have slopes of $d-1$,
$d$ and $d+1$.

## References

  * [1] E. Aamari, J. Kim, F. Chazal, B. Michel, A. Rinaldo, and L. Wasserman, Estimating the reach of a manifold, Electron. J. Stat. 13 (2019).
  * [2] H. Adams, M. Aminian, E. Farnell, M. Kirby, J. Mirth, R. Neville, C. Peterson and C. Shonkwiler, A fractal dimension for measures via persistent homology, in Topological data analysis—the Abel Symposium 2018, Abel Symp., 15, Springer 2020.
  * [3] M. Belkin, P. Niyogi, and V. Sindhwani. Manifold regularization: a geometric framework for learning from labeled and unlabeled examples, J. Mach. Learn. Res., (2006).
  * [4] M. Berger, A panoramic view of Riemannian geometry, Springer-Verlag, Berlin, 2003.
  * [5] J.-D. Boissonnat and A. Ghosh, Manifold reconstruction using tangential Delaunay complexes, Discrete Comput. Geom., 51, (2014).
  * [6] J.-D. Boissonnat, A. Lieutier, and M. Wintraecken, The reach, metric distortion, geodesic convexity and the variation of tangent spaces, J. Appl. Comput. Topol. 3, (2019).
  * [7] P. Breiding, S. Kališnik, B. Sturmfels and M. Weinstein, Learning algebraic varieties from samples, Rev. Mat. Complut. 31, (2018).
  * [8] M. Bridson and A. Haefliger, Metric spaces of non-positive curvature, Grundlehren der mathematischen Wissenschaften 319. Springer-Verlag, 1999.
  * [9] F. Camastra, Data dimensionality estimation methods: a survey, Pattern Recognition, Volume 36, (2003).
  * [10] F. Camastra and A. Staiano, Intrinsic dimension estimation: advances and open problems, Information Sciences 328, (2016).
  * [11] J. Cheeger and D. Ebin, Comparison theorems in Riemannian geometry, AMS Chelsea Publishing, 2008.
  * [12] M. Díaz, A. Quiroz, and M. Velasco, Local angles and dimension estimation from data on manifolds, J. Multivariate Anal. 173, (2019).
  * [13] J.P. Eckmann and D. Ruelle, Fundamental liminations for estimating dimensions and lyapounov exponents in dynamical systems, Physica D-56, (1992).
  * [14] H. Federer, Curvature measures, Trans. Am. Math. Soc. 93, (1959).
  * [15] C. Fefferman, S. Mitter, and H. Narayanan, Testing the manifold hypothesis, J. Amer. Math. Soc. 29, (2016).
  * [16] K. Fukunaga, Intrinsic dimensionality extraction, in Handbook of Statistics, Volume 2, 1982.
  * [17] P. Grassberger, Do climatic attractors exist?, Nature 323, (1986).
  * [18] P. Grassberger and I. Procaccia, Measuring the strangeness of strange attractors, Physica D9, (1983).
  * [19] L. Grillet, https://github.com/lgrillet/dim-estimation
  * [20] J. Lee and M. Verleysen, Nonlinear dimensionality reduction, Information Science and Statistics, Springer 2010.
  * [21] H. Narayanan and S. Mitter, Sample complexity of testing the manifold hypothesis, in Advances in Neural Information Processing Systems 23, Neural Information Processing Systems, 2011.
  * [22] P. Niyogi, S. Smale and S. Weinberger, Finding the homology of submanifolds with high confidence from random samples, Discrete Comput. Geom. 39, (2008).
  * [23] I. Procaccia, Complex of just complicated?, in Nature 333, (1998).
  * [24] D. Simpelaere, Correlation Dimension, Journal of Statistical Physics volume 90, (1998).
  * [25] F. Takens, Invariants related to dimension and entropy, in Atas do 13 coloqkio brasileiro de matématica, Rio de Janeiro, 1983.
  * [26] F. Takens, On the numerical determination of the dimension of an attractor, in Dynamical systems and bifurcations (Groningen, 1984), Lecture Notes in Math., 1125, Springer, 1985.
  * [27] J. Theiler, Estimating fractal dimension, J. Opt. Soc. Am. A 7, (1990).
  * [28] S. Weinberger, The complexity of some topological interference problems, Found. Comput. Math. 14, (2014).