Estimating the Reach of a Manifold

Various problems in manifold estimation make use of a quantity called the reach, denoted by $\tau\_M$, which is a measure of the regularity of the manifold. This paper is the first investigation into the problem of how to estimate the reach. First, we study the geometry of the reach through an approximation perspective. We derive new geometric results on the reach for submanifolds without boundary. An estimator $\hat{\tau}$ of $\tau\_{M}$ is proposed in a framework where tangent spaces are known, and bounds assessing its efficiency are derived. In the case of i.i.d. random point cloud $\mathbb{X}\_{n}$, $\hat{\tau}(\mathbb{X}\_{n})$ is showed to achieve uniform expected loss bounds over a $\mathcal{C}^3$-like model. Finally, we obtain upper and lower bounds on the minimax rate for estimating the reach.


Background and Related Work
Manifold estimation has become an increasingly important problem in statistics and machine learning. There is now a large literature on methods and theory for estimating manifolds. See, for example, [29,24,23,10,31,8,25].
Estimating a manifold, or functionals of a manifold, requires regularity conditions. In nonparametric function estimation, regularity conditions often take the form of smoothness constraints. In manifold estimation problems, a common assumption is that the reach τ M of the manifold M is non-zero.
First introduced by Federer [21], the reach τ M of a set M ⊂ R D is the largest number such that any point at distance less than τ M from M has a unique nearest point on M . If a set has its reach greater than τ min > 0, then one can roll freely a ball of radius τ min around it [14]. The reach is affected by two factors: the curvature of the manifold and the width of the narrowest bottleneck-like structure of M , which quantifies how close M is from being self-intersecting.
Positive reach is the minimal regularity assumption on sets in geometric measure theory and integral geometry [22,33]. Sets with positive reach exhibit a structure that is close to being differential -the so-called tangent and normal cones. The value of the reach itself quantifies the degree of regularity of a set, with larger values associated to more regular sets. The positive reach assumption is routinely imposed in the statistical analysis of geometric structures in order to ensure good statistical properties [14] and to derive theoretical guarantees. For example, in manifold reconstruction, the reach helps formalize minimax rates [24,29]. The optimal manifold estimators of [1] implicitly use reach as a scale parameter in their construction. In homology inference [31,7], the reach drives the minimal sample size required to consistently estimate topological invariants. It is used in [15] as a regularity parameter in the estimation of the Minkowski boundary lengths and surface areas. The reach has also been explicitly used as a regularity parameter in geometric inference, such as in volume estimation [4] and manifold clustering [5]. Finally, the reach often plays the role of a scale parameter in dimension reduction techniques such as vector diffusions maps [32]. Problems in computational geometry such as manifold reconstruction also rely on assumptions on the reach [10].
In this paper we study the problem of estimating reach. To do so, we first provide new geometric results on the reach. We also give the first bounds on the minimax rate for estimating reach.
There are very few papers on this problem. When the embedding dimension is 3, the estimation of the local feature size (a localized version of the reach) was tackled in a deterministic way in [18]. To some extent, the estimation of the medial axis (the set of points that have strictly more than one nearest point on M ) and its generalizations [16,6] can be viewed as an indirect way to estimate the reach. A test procedure designed to validate whether data actually comes from a smooth manifold satisfying a condition on the reach was developed in [23]. The authors derived a consistent test procedure, but the results do not permit any inference bound on the reach.

Outline
In Section 2 we provide some differential geometric background and define the statistical problem at hand. New geometric properties of the reach are derived in Section 3, and their consequences for its inference follow in Section 4 in a setting where tangent spaces are known. We study minimax rates in Section 5. An extension to a model where tangent spaces are unknown is discussed in Section 6, and we conclude with some open questions in Section 7. For sake of readability, the proofs are given in the Appendix.

Notions of Differential Geometry
In what follows, D ≥ 2 and R D is endowed with the Euclidean inner product ·, · and the associated norm · . The associated closed ball of radius r and center x is denoted by B(x, r). We will consider compact connected submanifolds M of R D of fixed dimension 1 ≤ d < D and without boundary [19]. For every point p in M , the tangent space of M at p is denoted by T p M : it is the d-dimensional linear subspace of R D composed of the directions that M spans in the neighborhood of p. Besides the Euclidean structure given by R D ⊃ M , a submanifold is endowed with an intrinsic distance induced by the ambient Euclidean one, and called the geodesic distance. Given a smooth path c : [a, b] → M , the length of c is defined as Length(c) = b a c (t) dt. One can show [19] that there exists a path γ p→q of minimal length joining p and q. Such an arc is called geodesic, and the geodesic distance between p and q is given by d M (p, q) = Length(γ p→q ). We let B M (p, s) denote the closed geodesic ball of center p ∈ M and of radius s. A geodesic γ such that γ (t) = 1 for all t is called arc-length parametrized. Unless stated otherwise, we always assume that geodesics are parametrized by arc-length. For all p ∈ M and all unit vectors v ∈ T p M , we denote by γ p,v the unique arclength parametrized geodesic of M such that γ p,v (0) = p and γ p,v (0) = v. The exponential map is defined as exp p (vt) = γ p,v (t). Note that from the compactness of M , exp p : be the angle between u and v.

Reach
First introduced by Federer [21], the reach regularity parameter is defined as follows. Given a closed subset A ⊂ R D , the medial axis M ed(A) of A is the subset of R D consisting of the points that have at least two nearest neighbors on A. Namely, denoting by d(z, A) = inf p∈A p − z the distance function to A, The reach of A is then defined as the minimal distance from A to M ed(A).
Definition 2.1. The reach of a closed subset A ⊂ R D is defined as Some authors refer to τ −1 A as the condition number [31,32]. From the definition of the medial axis in (2.1), the projection π A (x) = arg min p∈A p − x onto A is well defined outside M ed(A). The reach is the largest distance ρ ≥ 0 such that π A is well defined on the ρ-offset x ∈ R D |d(x, A) < ρ . Hence, the reach condition can be seen as a generalization of convexity, since a set A ⊂ R D is convex if and only if τ A = ∞.
In the case of submanifolds, one can reformulate the definition of the reach in the following manner.
Theorem 2.2 (Theorem 4.18 in [21]). For all submanifolds M ⊂ R D , This formulation has the advantage of involving only points on M and its tangent spaces, while (2.2) uses the distance to the medial axis M ed(M ), which is a global quantity. The formula (2.3) will be the starting point of the estimator proposed in this paper (see Section 4). Figure 1: Geometric interpretation of quantities involved in (2.3).
The ratio appearing in 2.3 can be interpreted geometrically, as suggested in Figure 1. This ratio is the radius of an ambient ball, tangent to M at p and passing through q. Hence, at a differential level, the reach gives a lower bound on the radii of curvature of M . Equivalently, τ −1 M is an upper bound on the curvature of M . Proposition 2.3 (Proposition 6.1 in [31]). Let M ⊂ R D be a submanifold, and γ p,v an arc-length parametrized geodesic of M . Then for all t, In analogy with function spaces, the class M ⊂ R D |τ M ≥ τ min > 0 can be interpreted as the H older space C 2 (1/τ min ). In addition, as illustrated in Figure  2, the condition τ M ≥ τ min > 0 also prevents bottleneck structures where M is nearly self-intersecting. This idea will be made rigorous in Section 3.

Statistical Model and Loss
Let us now describe the regularity assumptions we will use throughout. To avoid arbitrarily irregular shapes, we consider submanifolds M with their reach lower bounded by τ min > 0. Since the parameter of interest τ M is a C 2 -like quantity, it is natural -and actually necessary, as we shall see in Proposition 2.9 -to require an extra degree of smoothness. For example, by imposing an upper bound on the third order derivatives of geodesics.
Definition 2.4. We let M d,D τ min ,L denote the set of compact connected d-dimensional submanifolds M ⊂ R D without boundary such that τ M ≥ τ min , and for which every arc-length parametrized geodesic γ p,v is C 3 and satisfies It is important to note that any compact d-dimensional C 3 -submanifold M ⊂ R D belongs to such a class M d,D τ min ,L , provided that τ min ≤ τ M and that L is large enough. Note also that since the third order condition γ p,v (0) ≤ L needs to hold for all (p, v), we have in particular that γ p,v (t) ≤ L for all t ∈ R. To our knowledge, such a quantitative C 3 assumption on the geodesic trajectories has not been considered in the computational geometry literature.
Any submanifold M ⊂ R D of dimension d inherits a natural measure vol M from the d-dimensional Hausdorff measure H d on R D [22, p. 171]. We will consider distributions Q that have densities with respect to vol M that are bounded away from zero.
In order to focus on the geometric aspects of the reach, we will first consider the case where tangent spaces are observed at all the sample points. We let G d,D denote the Grassmanian of dimension d of R D , that is the set of all d-dimensional linear subspaces of R D .
Definition 2.6. For any distribution Q ∈ Q d,D τ min ,L,f min with support M we associate the distribution P of the random variable (X, T X M ) on R D × G d,D , where X has distribution Q. We let P d,D τ min ,L,f min denote the set of all such distributions. Formally, one can write P (dx dT ) = δ TxM (dT )Q(dx), where δ · denotes the Dirac measure. An i.i.d. n-sample of P is of the form (X 1 , T 1 ), . . . , (X n , T n ) ∈ R D × G d,D , where X 1 , . . . , X n is an i.i.d. n-sample of Q and T i = T X i M with M = supp(Q). For a distribution Q with support M and associated distribution P on R D × G d,D , we will write τ P = τ Q = τ M , with a slight abuse of notation. Note that the model does not explicitly impose an upper bound on τ M . Such an upper bound would be redundant, since the lower bound on f min does impose such an upper bound, as we now state in the following result. The proof relies on a volume argument (Lemma A.2 in the Appendix, leading to a bound on the diameter of M , and on a topological argument (Lemma A.5 in the Appendix to link the reach and the diameter. Proposition 2.7. Let M ⊂ R D be a connected closed d-dimensional manifold, and let Q be a probability distribution with support M . Assume that Q has a density f with respect to the Hausdorff measure on M such that inf x∈M f (x) ≥ f min > 0. Then, for some constant C d > 0 depending only on d.
To simplify the statements and the proofs, we focus on a loss involving the condition number. Namely, we measure the error with the loss In other words, we will consider the estimation of the condition number τ −1 M instead of the reach τ M .
Remark 2.8. For a distribution P ∈ P d,D τ min ,L,f min , Proposition 2.7 asserts that τ min ≤ τ P ≤ τ max := (C d /f min ) 1/d . Therefore, in an inference set-up, we can always restrict to estimatorsτ within the bounds τ min ≤τ ≤ τ max . Consequently, so that the estimation of the reach τ P is equivalent to the estimation of the condition number τ −1 P , up to constants. With the statistical framework developed above, we can now see explicitly why the third order condition γ ≤ L < ∞ is necessary. Indeed, the following Proposition 2.9 demonstrates that relaxing this constraint -i.e. setting L = ∞ -renders the problem of reach estimation intractable. Its proof is to be found in Section D.3 of the Appendix. Below, σ d stands for the volume of the ddimensional unit sphere S d .
Thus, one cannot expect to derive consistent uniform approximation bounds for the reach solely under the condition τ M ≥ τ min . This result is natural, since the problem at stake is to estimate a differential quantity of order two. Therefore, some notion of uniform C 3 regularity is needed.

Geometry of the Reach
In this section, we give a precise geometric description of how the reach arises. In particular, below we will show that the reach is determined either by a bottleneck structure or an area of high curvature (Theorem 3.4). These two cases are referred to as global reach and local reach, respectively. All the proofs for this section are to be found in Section B of the Appendix.
Consider the formulation (2.2) of the reach as the infimum of the distance between M and its medial axis M ed(M ). By definition of the medial axis (2.1), if the infimum is attained it corresponds to a point z 0 in M ed(M ) at distance τ M from M , which we call an axis point. Since z 0 belongs to the medial axis of M , it has at least two nearest neighbors q 1 , q 2 on M , which we call a reach attaining pair (see Figure 3b). By definition, q 1 and q 2 belong to B(z 0 , τ M ) and cannot be farther than 2τ M from each other. We say that (q 1 , q 2 ) is a bottleneck of M in the extremal case q 2 − q 1 = 2τ M of antipodal points of B(z 0 , τ M ) (see Figure  3a). Note that the ball B(z 0 , τ M ) meets M only on its boundary ∂B(z 0 , τ M ). • A pair of points (q 1 , q 2 ) in M is called reach attaining if there exists z 0 ∈ M ed(M ) such that q 1 , q 2 ∈ B(z 0 , τ M ). We call z 0 the axis point of (q 1 , q 2 ), and q 1 − q 2 ∈ (0, 2τ M ] its size. • A reach attaining pair (q 1 , q 2 ) ∈ M 2 is said to be a bottleneck of M if its size is 2τ M , that is q 1 − q 2 = 2τ M .
As stated in the following Lemma 3.2, if a reach attaining pair is not a bottleneck -that is q 1 − q 2 < 2τ M , as in Figure 3b -, then M contains an arc of a circle of radius τ M . In this sense, this "semi-local" case -when q 1 − q 2 can be arbitrarily small -is not generic. Though, we do not exclude this case in the analysis.
be their associated axis point, and write c z 0 (q 1 , q 2 ) for the arc of the circle with center z 0 and endpoints as q 1 and q 2 .
Then c z 0 (q 1 , q 2 ) ⊂ M , and this arc (which has constant curvature 1/τ M ) is the geodesic joining q 1 and q 2 .
In particular, in this "semi-local" situation, since τ −1 M is the norm of the second derivative of a geodesic of M (the exhibited arc of the circle of radius τ M ), the reach can be viewed as arising from directional curvature. Now consider the case where the infimum (2.2) is not attained. In this case, the following Lemma 3.3 asserts that τ M is created by curvature.
To summarize, there are three distinct geometric instances in which the reach may be realized: From now on, we will treat the first case separately from the other two. We are now in a position to state the main result of this section. It is a straightforward consequence of Lemma 3.2 and Lemma 3.3. • (Global Case) M has a bottleneck (q 1 , q 2 ) ∈ M 2 , that is, there exists z 0 ∈ M ed(M ) such that q 1 , q 2 ∈ ∂B(z 0 , τ M ) and q 1 − q 2 = 2τ M .
• (Local Case) There exists q 0 ∈ M and an arc-length parametrized geodesic γ 0 such that γ 0 (0) = q 0 and γ 0 (0) = 1 τ M . Let us emphasize the fact that the global case and the local case of Theorem 3.4 are not mutually exclusive. Theorem 3.4 provides a description of the reach as arising from global and local geometric structures that, to the best of our knowledge, is new. Such a distinction is especially important in our problem. Indeed, the global and local cases may yield different approximation properties and require different statistical analyses. However, since one does not know a priori whether the reach arises from a global or a local structure, an estimator of τ M should be able to handle both cases simultaneously.

Reach Estimator and its Analysis
In this section, we propose an estimatorτ (·) for the reach and demonstrate its properties and rate of consistency under the loss (2.4). For the sake of clarity in the analysis, we assume the tangent spaces to be known at every sample point. This assumption will be relaxed in Section 6.
We rely on the formulation of the reach given in (2.3) (see also Figure 1), and defineτ as a plugin estimator as follows: given a point cloud X ⊂ M , In particular, we haveτ (M ) = τ M . Since the infimum (4.1) is taken over a set X smaller than M ,τ (X) always overestimates τ M . In fact,τ (X) is decreasing in the number of distinct points in X, a useful property that we formalize in the following result, whose proof is immediate.
We now derive the rate of convergence ofτ . We analyze the global case (Section 4.1) and the local case (Section 4.2) separately. In both cases, we first determine the performance of the estimator in a deterministic framework, and then derive an expected loss bounds whenτ is applied to a random sample.
Respectively, the proofs for Section 4.1 and Section 4.2 are to be found in Section C.1 and Section C.2 of the Appendix.

Global Case
Consider the global case, that is, M has a bottleneck structure (Theorem 3.4). Then the infimum (2.3) is achieved at a bottleneck pair (q 1 , q 2 ) ∈ M 2 . When X contains points that are close to q 1 and q 2 , one may expect that the infimum over the sample points should also be close to (2.3): that is, thatτ (X) should be close to τ M . Proposition 4.2. Let M ⊂ R D be a submanifold with reach τ M > 0 that has a bottleneck (q 1 , q 2 ) ∈ M 2 (see Definition 3.1), and X ⊂ M . If there exist x, y ∈ X with q 1 − x < τ M and q 2 − y < τ M , then The error made byτ (X) decreases linearly in the maximum of the distances to the critical points q 1 and q 2 . In other words, the radius of the tangent sphere in Figure 1 grows at most linearly in t when we perturb by t < τ M its basis point p = q 1 and the point q = q 2 it passes through.
Based on the deterministic bound in Proposition 4.2, we can now give an upper bound on the expected loss under the model P d,D τ min ,L,f min . We recall that, throughout the paper, X n = {X 1 , . . . , X n } is an i.i.d. sample with common distribution Q associated to P (see Definition 2.6).

Local Case
Consider now the local case, that is, there exists q 0 ∈ M and v 0 ∈ T q 0 M such that the geodesic γ 0 = γ q 0 ,v 0 has second derivative γ 0 (0) = 1/τ M (Theorem 3.4). Estimating τ M boils down to estimating the curvature of M at q 0 in the direction v 0 .
We first relate directional curvature to the increment y−x 2 2d(y−x,TxM ) involved in the estimatorτ (4.1). Indeed, since the latter quantity is the radius of a sphere tangent at x and passing through y (Figure 1), it approximates the radius of curvature in the direction y − x when x and y are close. For x, y ∈ M , we let γ x→y denote the arc-length parametrized geodesic joining x and y, with the convention γ x→y (0) = x.
Let us now state how directional curvatures are stable with respect to perturbations of the base point and the direction. We let κ p denote the maximal directional curvature of M at p ∈ M , that is, . Let γ 0 be a geodesic such that γ 0 (0) = q 0 and γ 0 (0) = κ q 0 . Write θ x := ∠(γ 0 (0), γ q 0 →x (0)), θ y := ∠(γ 0 (0), γ q 0 →y (0)), and suppose that |θ x − θ y | ≥ π 2 . Then, γ x→y (0) In particular, geodesics in a neighborhood of q 0 with directions close to v 0 have curvature close to 1 τ M . A point cloud X sampled densely enough in M would contain points in this neighborhood. Hence combining Lemma 4.4 and Lemma 4.5 yields the following deterministic bound in the local case.
In other words, since the reach boils down to directional curvature in the local case,τ performs well if it is given as input a pair of points x, y which are close to the point q 0 realizing the reach, and almost aligned with the direction of interest v 0 . Note that the error bound in the local case (Proposition 4.6) is very similar to that of the global case (Proposition 4.2) with an extra alignment term sin 2 (|θ x − θ y |) . This alignment term appears since, in the local case, the reach arises from directional curvature τ M = γ q 0 ,v 0 (0) (Theorem 3.4). Hence, it is natural that the accuracy ofτ (X) depends on how precisely X samples the neighborhood of q 0 in the particular direction v 0 .
Similarly to the analysis of the global case, the deterministic bound in Proposition 4.6 yields a bound on the risk ofτ (X n ) when X n = {X 1 , . . . , X n } is random.
where C τ min ,d,L,f min ,p depends only on τ min , d, L, f min and p.
This statement follows from Proposition 4.6 together with the estimate of the probability of two points being drawn in a neighborhood of q 0 and subject to an alignment constraint. Proposition 4.3 and 4.7 yield a convergence rate ofτ (X n ) which is slower in the local case than in the global case. Recall that from Theorem 3.4, the reach pertains to the size of a bottleneck structure in the global case, and to maximum directional curvature in the local case. To estimate the size of a bottleneck, observing two points close to each point in the bottleneck gives a good approximation. However, for approximating maximal directional curvature, observing two points close to the curvature attaining point is not enough, but they should also be aligned with the highly curved direction. Hence, estimating the reach may be more difficult in the local case, and the difference in the convergence rates of Proposition 4.3 and 4.7 accords with this intuition.
Finally, let us point out that in both cases, neither the convergence rates nor the constants depend on the ambient dimension D.

Minimax Estimates
In this section we derive bounds on the minimax risk R n of the estimation of the reach over the class P d,D τ min ,L,f min , that is where the infimum ranges over all estimatorsτ n (X 1 , T X 1 ), . . . , (X n , T Xn ) based on an i.i.d. sample of size n with the knowledge of the tangent spaces at sample points.
The rate of convergence of the plugin estimatorτ (X n ) studied in the previous section leads to an upper bound on R n , which we state here for completeness.
Theorem 5.1. For all n ≥ 1, for some constant C τ min ,d,L,f min ,p depending only on τ min , d, L, f min and p.
We now focus on deriving a lower bound on the minimax risk R n . The method relies on an application of Le Cam's Lemma [34]. In what follows, let T V P, P = 1 2 |dP − dP | denote the total variation distance between P and P , where dP, dP denote the respective densities of P, P with respect to any dominating measure. Since |x − z| p + |z − y| p ≥ 2 1−p |x − y| p , the following version of Le Cam's lemma results from Lemma 1 in [34] and (1 − T V (P n , P n )) ≥ (1 − T V (P, P )) n .
Lemma 5.2 (Le Cam's Lemma). Let P, P ∈ P d,D τ min ,L,f min with respective supports M and M . Then for all n ≥ 1, Lemma 5.2 implies that in order to derive a lower bound on R n one needs to consider distributions (hypotheses) in the model that are stochastically close to each other -i.e. with small total variation distance -but for which the associated reaches are as different as possible. A lower bound on the minimax risk over P d,D τ min ,L,f min requires the hypotheses to belong to the class. Luckily, in our problem it will be enough to construct hypotheses from the simpler class Q d,D τ min ,L,f min . Indeed, we have the following isometry result between Q d,D τ min ,L,f min and P d,D τ min ,L,f min for the total variation distance, as proved in Section D.2 in the Appendix. In order to construct hypotheses in Q d,D τ min ,L,f min we take advantage of the fact that the class M d,D τ min ,L has good stability properties, which we now describe. Here, since submanifolds do not have natural parametrizations, the notion of perturbation can be well formalized using diffeomorphisms of the ambient space R D ⊃ M . Given a smooth map Φ : R D → R D , we denote by d i x Φ its differential of order i at x. Given a tensor field A between Euclidean spaces, let A op = sup x A x op , where A x op is the operator norm induced by the Euclidean norm. The next result states, informally, that the reach and geodesics third derivatives of a submanifold that is perturbed by a diffeomorphism that is C 3 -close to the identity map do not change much. The proof of Proposition 5.4 can be found in Section D.3 of the Appendix.
Hence, applying Lemma 5.2 with the hypotheses P, P associated to Q, Q of Proposition 5.5, and taking 12 ( /2τ min ) d = 1/n, together with Lemma 5.3, yields the following lower bound. Here, the assumptions on the parameters L and f min are necessary for the model to be rich enough. Roughly speaking, they ensure at least that a sphere of radius 2τ min belongs to the model. From Proposition 5.6, the plugin estimationτ (X n ) provably achieves the optimal rate in the global case (Theorem 4.3) up to numerical constants. In the local case (Theorem 4.7) the rate obtained presents a gap, yielding a gap in the overall rate. As explained above (Section 4.2), the slower rate in the local case is a consequence of the alignment required in order to estimate directional curvature. Though, let us note that in the one-dimensional case d = 1, the rate of Proposition 5.6 matches the convergence rate ofτ (X n ) (Theorem 5.1). Indeed, for curves, the alignment requirement is always fulfilled. Hence, the rate is exactly n −p for d = 1, andτ (X n ) is minimax optimal.
Here, again, neither the convergence rate nor the constant depend on the ambient dimension D.

Towards Unknown Tangent Spaces
So far, in our analysis we have used the key assumption that both the point cloud and the tangent spaces were jointly observed. We now focus on the more realistic framework where only points are observed. We once again rely on the formulation of the reach given in Theorem 2.3 and consider a new plug-in estimator in which the true tangent spaces are replaced by estimated ones. Namely, given a point cloud X ⊂ R D and a family T = {T x } x∈X of linear subspaces of R D indexed by X, the estimator is defined aŝ . (6.1) In particular,τ (X) =τ (X, T X M ), where T X M = {T x M } x∈X . Adding uncertainty on tangent spaces in (6.1) does not change drastically the estimator as the formula is stable with respect to T . We state this result quantitatively in the following Proposition 6.1, the proof of which can be found in Section E of the Appendix. In what follows, the distance between two linear subspaces U, V ∈ G d,D is measured with their principal angle π U − π V op . Proposition 6.1. Let X ⊂ R D and T = {T x } x∈X ,T = {T x } x∈X be two families of linear subspaces of R D indexed by X. Assume X to be δ-sparse, T andT to be θ-close, in the sense that In other words, the map T →τ (X, T ) −1 is smooth, provided that the basis point cloud X contains no zone of accumulation at a too small scale δ > 0. As a consequence, under the assumptions of Proposition 6.1, the bounds on τ (X) −1 − τ M −1 of Proposition 4.2 and Proposition 4.6 still hold with an extra error term 2 sin θ/δ if we replaceτ (X) byτ (X, T ).
For an i.i.d. point cloud X n asymptotic rates of tangent space estimation derived in C 3 -like models can be found in [13,32,2], yielding bounds on sin θ. In that case, the typical scale of minimum interpoint distance is δ n −2/d , as stated in the asymptotic result Theorem 2.1 in [27] for the flat case of R d . However, the typical covering scale of M used in the global case (Theorem 4.3) is ε (1/n) 1/d . It appears that we can sparsify the point cloud X n -that is, removing accumulation points -while preserving the covering property at scale ε = 2δ (log n/n) 1/d . This can be performed using the farthest point sampling algorithm [1, Section 3.3]. Such a sparsification pre-processing allows to lessen the possible instability ofτ (X n , ·) −1 . Though, whether the alignment property used in the local case (Theorem 4.7) is preserved under sparsification remains to be investigated.

Conclusion and Open Questions
In the present work, we gave new insights on the geometry of the reach. Inference results were derived in both deterministic and random frameworks. For i.i.d. samples, non-asymptotic minimax upper and lower bounds were derived under assumptions on the third order derivative of geodesic trajectories. Let us conclude with some open questions.
• The minimax upper and lower bounds given in Theorem 5.1 and Theorem 5.6 do not match. They are yet to be sharpened.
• In practice, since large reach ensures regularity, one may be interested with having a lower bound on the reach τ M . Giving the limiting distribution of the statisticτ (X n ) would allow to derive asymptotic confidence intervals for τ M .
• Other regularity parameters such as local feature size [10] and λ-reach [12] could be relevant to estimate, as they are used as tuning parameters in computational geometry techniques. Geodesic joining x to y with γ x→y (0 Operator norm A Some Technical Results on the Model

A.1 Metric Properties
This section garners geometric lemmas on embedded manifolds in the Euclidean space that are related to the reach, and that will be used several times in the proofs. (iv) For all p ∈ M , the map exp p : (vi) Let γ be a geodesic at p ∈ M , and P t the parallel transport operator along γ. Then for all t < πτ M and v ∈ T p M , Proof of Proposition A.1. (i) is stated in Proposition 2.1 in [31], yielding (ii) from Corollary 1.4 in [3]. (iii) follows using (i) again and the Gauss equation [19, p. 130]. (iv) is derived from (iii) by a direct application of Lemma 8 in [20]. (v) follows from (iv) and Lemma 6 in [5]. All that remain to be showed is (vi). For this, assume without loss of generality that v = 1. Let g : [0, t] → S d−1 be defined by g(s) = P s (v). Let u ∈ R D be a unit vector and denoting by∇ the ambient derivative. We may write g (s), u = ∇ γ (s) P s (w), u = II(γ (s), P s (w)), u .
Since g is a curve on S d−1 , this implies

A.2 Comparing Reach and Diameter
Let us prove Proposition 2.7. For this aim, we first state the following analogous bound on the (Euclidean) diameter.
for some constant C d > 0 depending only on d.
Proposition A.3. If K ⊂ R D is not homotopy equivalent to a point,

Proof of Proposition A.3. Follows from a straightforward combination of Lemma A.4 and Lemma A.5.
We recall that for two compact subsets A, B ⊂ R D , the Hausdorff distance [11, p. 252] between them is defined by We denote by conv(·) the closed convex hull of a set.
Proof of Lemma A.4. It is a straightforward corollary of Jung's Theorem 2.10.41 in [22], which states that K is contained in a (unique) closed ball with (minimal) radius at most Proof of Lemma A.5. Let us prove the contrapositive. For this, assume that τ K > d H (K, conv(K)). Then, Therefore, the map π K : conv(K) → K is well defined and continuous, so that K is a retract of conv(K) (see Chapter 0 in [26]). Therefore, K is homotopy equivalent to a point, since the convex set conv(K) is.
We are now in position to prove Proposition 2.7.
Proof of Proposition 2.7. From Theorem 3.26 in [26], M has a non trivial homology group of dimension d over Z/2Z, so that it cannot be homotopy equivalent to a point. Therefore (ii) Let 0 < r < q < ∞ be fixed. Let x, y / ∈ M ed(M ) ∪ M be such that d(x, M ) ∨ d(y, M ) ≤ r and Proof of Lemma B.1. The proof of (i) follows that of Theorem 4.8 (7) in [21]. Let v = x−a x−a and S = {t|π M (a + tv) = a}. As x − a > 0 belongs to S, sup S > 0 and from [21,Theorem 4.8 (6)] we get sup S ≥ τ M (a, v). Moreover, for 0 < t ∈ S, Developing and rearranging the square of previous inequality yields On the other hand, the proof of (ii) follows that of Theorem 4.8 (8) in [21]. Writing a = π M (x) and b = π M (y), the previous point yields As a consequence, hence the result.
Then for any fixed n ∈ N and t 1 , Noticing furthermore that For any fixed k and 0 ≤ j ≤ k, set t k,j = 2j−k k t 0 . The inequality above yields, Moreover, the γ n 's are curves joining q 1 to q 2 with images γ n ([−t 0 , t 0 ]) ⊂ R D \ • B(z 0 , τ M ), so that their lengths are at most that of the arc of great circle c z 0 (q 1 , q 2 ), that is Length (γ n ) ≥ Length (c z 0 (q 1 , q 2 )) = 2t 0 .
Lemma B.3. Let M be a compact manifold, and q 1 , q 2 ∈ M with q 1 = q 2 . Let (γ n ) n∈N be a sequence of curves on M joining q 1 and q 2 such that sup n Length(γ n ) < ∞ Then there exists a curve γ on M joining q 1 and q 2 such that Proof of Lemma B.3. Without loss of generality, we take the γ n 's to be arc length parametrized. For all n ∈ N, we let g n : [0, 1] → M be the reparametrization g n (t) = γ n (Length(γ n )t) . Notice that for all t ∈ [0, 1], the set (g n (t)) n∈N is contained in the compact set M , so that it is bounded uniformly in t. Moreover, writing K = sup n Length(γ n ) < ∞, we have that for all t 1 , t 2 ∈ [0, 1], Hence, the sequence (g n ) n∈N is pointwise bounded and equicontinuous. From Arzelà-Ascoli theorem [30,Theorem 45.4], there exists a curve γ : [0, 1] → M and subsequence (g n i ) i∈N converging uniformly to γ. For any fixed k and 1 ≤ j ≤ k, set t k,j = j/k. The (pointwise) convergence of (g n i ) i towards γ ensures that Furthermore, from the uniform convergence of (g n i ) i towards γ on the compact set [0, 1], hence the result. Then for all r ≤ τ M /2, To prove Lemma B.4 we need the following straightforward result.
Lemma B.5. Let U be a linear space and u ∈ U , n ∈ U ⊥ . If v = u + n + e, then Proof of Lemma B.4. First note that for all unit vector v ∈ T p M , γ p,v (r) belongs to B(p, r) ∩ M and, whenever 0 < r ≤ τ M 2 , Proposition A.1 (ii) ensures that γ p,v (r) = p. Therefore, it suffices to show that for all q ∈ B(p, r) ∩ M , there exists a unit tangent vector v ∈ T p M such that Let q ∈ B(p, r) ∩ M be different from p. Denoting t = d M (p, q) > 0, we let γ = γ p→q be the arc-length parametrized geodesic of minimal length such that γ(0) = p and γ(t) = q. γ exists from Proposition A.1 (ii) since r ≤ τ M 2 < πτ M . We will show that v = γ (0) provides the desired bound.
First, a Taylor expansion at zero of γ yields, Since γ (0) ∈ T p M ⊥ , Lemma B.5 shows that Therefore, This yields, Moreover, from q − p ≤ d M (p, q) and Proposition 6.3 in [31], we derive where the last two inequalities follow from elementary real analysis arguments. Therefore, we get t ≤ 2 q − p and Finally, using (2.3) we derive, Finally, the unit tangent bundle Proposition C.1. Let M ⊂ R D be a submanifold, and 0 < λ ≤ τ M . Assume that M has a reach attaining pair (q 1 , q 2 ) ∈ M 2 (see Definition 3.1) with q 1 − q 2 ≥ 2λ. Let X ⊂ M . If there exists x, y ∈ X with q 1 − x < λ and q 2 − y < λ, then depends only on the parameters τ M , λ, and is a decreasing function of τ M and λ when the other parameter is fixed.
Proof of Proposition C.1. The two left hand inequalities are a direct consequence of Corollary 4.1, let us then focus on the third one.
Without loss of generality, assume that q 1 − q 2 = 2λ. We set t to be equal to max {d M (q 1 , x), d M (q 2 , y)}, and z 1 := x + (q 2 − q 1 ). We have z 1 − x = q 2 − q 1 = 2λ and y − q 2 , q 1 − x ≤ t. Therefore, from the definition ofτ in (4.1) and the fact that the distance function to a linear space is 1-Lipschitz, we get The integration of the above bound gives where C τ M ,λ,f min ,d,p depends only on τ M , λ, f min , d, p, and is a decreasing function of τ M and λ when other parameters are fixed.