Efficient Weingarten Map and Curvature Estimation on Manifolds

TL;DR

This paper introduces an efficient method for estimating the Weingarten map from point cloud data on manifolds, providing theoretical analysis and practical applications in curvature estimation and data simplification.

Contribution

It presents a novel estimator for the Weingarten map with proven convergence properties and demonstrates its effectiveness on real and simulated data.

Findings

Convergence rate of the estimator established theoretically.

Successful application to curvature estimation.

Effective point cloud simplification achieved.

Abstract

In this paper, we propose an efficient method to estimate the Weingarten map for point cloud data sampled from manifold embedded in Euclidean space. A statistical model is established to analyze the asymptotic property of the estimator. In particular, we show the convergence rate as the sample size tends to infinity. We verify the convergence rate through simulated data and apply the estimated Weingarten map to curvature estimation and point cloud simplification to multiple real data sets.

Equations

$$T_p\mathcal{M}\oplus T_p^{\perp }\mathcal{M}={\mathbb{E}}^d$$

$$A_{\xi }:T_p\mathcal{M}\to T_p\mathcal{M}$$

$$A_{\xi }(X)=-({\overline{\nabla }}_X\tilde{\xi }{)}^{\top }$$

$$A_{{\tilde{\xi }}_p}(X)=-\mathrm{d}g(X)$$

$$(\tilde{\xi}_{q}-\tilde{\xi}_{p})^{\top}=-A_{\tilde{\xi}_{p}}\big{(}(q-p)^{\top}\big{)}+O(\|q-p\|^{2})$$

$$\tilde{\xi }(\mathbf{u})=\tilde{\xi }(0)+\sum _{i=1}^m{u}^i\frac{\partial \tilde{\xi }}{\partial u^i}(0)+O(\Vert \mathbf{u}{\Vert }^2)$$

$$\mathbf{r}(\mathbf{u})=\mathbf{r}(0)+\sum _{i=1}^m{u}^i\frac{\partial \mathbf{r}}{\partial u^i}(0)+O(\Vert \mathbf{u}{\Vert }^2)$$

$$A_{\tilde{\xi}_{p}}\big{(}\frac{\partial\mathbf{r}}{\partial u^{i}}(0)\big{)}=-\big{(}\overline{\nabla}_{\frac{\partial\mathbf{r}}{\partial u^{i}}(0)}\tilde{\xi}\big{)}^{\top}=-\big{(}\frac{\partial\tilde{\xi}}{\partial u^{i}}(0)\big{)}^{\top}$$

$$\displaystyle\big{(}\tilde{\xi}(\mathbf{u})-\tilde{\xi}(0)\big{)}^{\top}=$$

$$\text{Cov}=\sum _{j\in I_i}(x_j-{\bar{x}}_i)(x_j-{\bar{x}}_i{)}^{\mathrm{t}}$$

$${\tilde{\xi }}_{x_j}^{\alpha }=\sum _{\beta =1}^{d-m}⟨{\xi }_i^{\alpha },{\xi }_j^{\beta }⟩{\xi }_j^{\beta }$$

$$K_h(y)=\frac{1}{h^m}K(\frac{\Vert y\Vert }{h})$$

$$({\tilde{\xi }}_{x_j}^{\alpha }-{\xi }_i^{\alpha }{)}^{\mathrm{t}}{E}_i=-(x_j-x_i{)}^{\mathrm{t}}{E}_i{A}_{{\xi }_i^{\alpha }}+O(\Vert x_j-x_i{\Vert }^2)$$

$$\sum _{j=1}^n\Vert ({\tilde{\xi }}_{x_j}^{\alpha }+{\xi }_i^{\alpha }{)}^{\mathrm{t}}{E}_i-(x_j-x_i{)}^{\mathrm{t}}{E}_i{\tilde{A}}_{{\xi }_i^{\alpha }}{\Vert }^2{K}_h(x_j-x_i)$$

$${\tilde{\Delta }}_i=E_i^{\mathrm{t}}\left[x_1-x_i,\cdots ,x_n-x_i\right]$$

$${\tilde{\Xi }}_i^{\alpha }=E_i^{\mathrm{t}}\left[{\tilde{\xi }}_{x_1}^{\alpha }-{\xi }_i^{\alpha },\cdots ,{\tilde{\xi }}_{x_n}^{\alpha }-{\xi }_i^{\alpha }\right]$$

$$W_i=\text{diag}\{K_h(x_1-x_i),\cdots ,K_h(x_n-x_i)\}$$

$${\tilde{A}}_{{\xi }_i^{\alpha }}=-{\tilde{\Xi }}_i^{\alpha }{W}_i{\tilde{\Delta }}_i^{\mathrm{t}}({\tilde{\Delta }}_i{W}_i{\tilde{\Delta }}_i^{\mathrm{t}}{)}^{-1}$$

$$(\xi_{P}-\xi_{p})^{\top}=-A\big{(}(P-p)^{\top}\big{)}+\eta(P)\|P-p\|^{2}$$

$${\mathrm{argmin}}_{\mathbf{A}\in {\mathbb{R}}^{m\times m}}{\mathbb{E}}_f\left[\Vert \mathbf{\Xi }+\mathbf{A\Delta }{\Vert }^2{\mathbf{K}}_P\right]$$

$$F(\mathbf{A})$$

$$=\text{Tr}({\mathbb{E}}_f[\mathbf{\Xi }{\mathbf{\Xi }}^{\mathrm{t}}{\mathbf{K}}_P])+2\text{Tr}({\mathbb{E}}_f[\mathbf{\Xi }{\mathbf{\Delta }}^{\mathrm{t}}{\mathbf{K}}_P]\mathbf{A})+\text{Tr}(\mathbf{A}{\mathbb{E}}_f[\mathbf{\Delta }{\mathbf{\Delta }}^{\mathrm{t}}{\mathbf{K}}_P]{\mathbf{A}}^{\mathrm{t}})$$

$$\mathbf{A}=-{\mathbb{E}}_f[\mathbf{\Xi }{\mathbf{\Delta }}^{\text{t}}{\mathbf{K}}_P]({\mathbb{E}}_f[\mathbf{\Delta }{\mathbf{\Delta }}^{\text{t}}{\mathbf{K}}_P]{)}^{-1}=-\mathbf{L}{\mathbf{D}}^{-1}$$

$$\mathbf{L}=\left[\begin{array}{ccc}{\mathbb{E}}_f[(\mathbf{\Xi }{)}_1(\mathbf{\Delta }{)}_1{\mathbf{K}}_P]& \cdots & {\mathbb{E}}_f[(\mathbf{\Xi }{)}_1(\mathbf{\Delta }{)}_m{\mathbf{K}}_P]\\ ⋮& \ddots & ⋮\\ {\mathbb{E}}_f[(\mathbf{\Xi }{)}_m(\mathbf{\Delta }{)}_1{\mathbf{K}}_P]& \cdots & {\mathbb{E}}_f[(\mathbf{\Xi }{)}_m(\mathbf{\Delta }{)}_m{\mathbf{K}}_P]\end{array}\right]$$

$$\mathbf{D}=\left[\begin{array}{ccc}{\mathbb{E}}_f[(\mathbf{\Delta }{)}_1^2{\mathbf{K}}_P]& \cdots & {\mathbb{E}}_f[(\mathbf{\Delta }{)}_1(\mathbf{\Delta }{)}_m{\mathbf{K}}_P]\\ ⋮& \ddots & ⋮\\ {\mathbb{E}}_f[(\mathbf{\Delta }{)}_m(\mathbf{\Delta }{)}_1{\mathbf{K}}_P]& \cdots & {\mathbb{E}}_f[(\mathbf{\Delta }{)}_m^2{\mathbf{K}}_P]\end{array}\right]$$

$$\tilde{A}=-\tilde{L}{\tilde{D}}^{-1}$$

$$\tilde{L}=\left[\begin{array}{ccc}\frac{1}{n}{\sum }_{i=1}^n({\Xi }_i{)}_1({\Delta }_i{)}_1{\mathbf{K}}_{x_i}& \cdots & \frac{1}{n}{\sum }_{i=1}^n({\Xi }_i{)}_1({\Delta }_i{)}_m{\mathbf{K}}_{x_i}\\ ⋮& \ddots & ⋮\\ \frac{1}{n}{\sum }_{i=1}^n({\Xi }_i{)}_m({\Delta }_i{)}_1{\mathbf{K}}_{x_i}& \cdots & \frac{1}{n}{\sum }_{i=1}^n({\Xi }_i{)}_m({\Delta }_i{)}_m{\mathbf{K}}_{x_i}\end{array}\right]$$

$$\tilde{D}=\left[\begin{array}{ccc}\frac{1}{n}{\sum }_{i=1}^n({\Delta }_i{)}_1^2{\mathbf{K}}_{x_i}& \cdots & \frac{1}{n}{\sum }_{i=1}^n({\Delta }_i{)}_1({\Delta }_i{)}_m{\mathbf{K}}_{x_i}\\ ⋮& \ddots & ⋮\\ \frac{1}{n}{\sum }_{i=1}^n({\Delta }_i{)}_m({\Delta }_i{)}_1{\mathbf{K}}_{x_i}& \cdots & \frac{1}{n}{\sum }_{i=1}^n({\Delta }_i{)}_m^2{\mathbf{K}}_{x_i}\end{array}\right]$$

$${\mathrm{argmin}}_{\tilde{A}\in {\mathbb{R}}^{m\times m}}\sum _{i=1}^n\Vert {\Xi }_i-\tilde{A}{\Delta }_i{\Vert }^2{\mathbf{K}}_{x_i}$$

$$\text{MSE}=\mathbb{E}_{f}\left[\|\widetilde{A}-A\|_{F}^{2}\right]\leq 2\bigg{(}\underbrace{\mathbb{E}_{f}\left[\|\widetilde{A}-\mathbf{A}\|_{F}^{2}\right]}_{\text{Variance}}+\underbrace{\|A-\mathbf{A}\|_{F}^{2}}_{\text{Bias}^{2}}\bigg{)}$$

$$h\le R_p(h)\le h+C_{R_p}{h}^3$$

$$\mathbf{D}=h^2\left(f(p){\int }_{\Vert \mathbf{z}\Vert \le 1}(z^1{)}^{2}K(\Vert \mathbf{z}\Vert )\mathrm{d}\mathbf{z}\right)(I+o(1))$$

$${\mathbf{D}}^{-1}=\frac{1}{h^2\left(f(p){\int }_{\Vert \mathbf{z}\Vert \le 1}(z^1{)}^{2}K(\Vert \mathbf{z}\Vert )\mathrm{d}\mathbf{z}\right)}(I+o(1))$$

$${\mathbb{E}}_f[(\mathbf{\Delta }{)}_k(\mathbf{\Delta }{)}_l{\mathbf{K}}_P]={\int }_{\mathcal{M}}((P-p)\cdot e_k)((P-p)\cdot e_l)K_h(P-p)f(P)dv$$

Full text

1 Introduction

In the literature of manifold learning, data in high dimensional Euclidean space are always assumed to lie on a low dimensional manifold ${\mathcal{M}}$. The inference of geometry and topology of ${\mathcal{M}}$ becomes critical in the understanding and application of the data such as dimension reduction, clustering and visualization. Topological quantities often involve persistent homology [14, 45], homotopy groups [21] and fundamental groups [4], while geometric quantities include intrinsic dimension [24], tangent space [36], geodesics [23], Laplacian operators [6, 7] and curvature [1, 11]. There are also hot topics on the estimation of volume and support. This paper focuses on the estimation of the Weingarten map, or the second fundamental form, for point clouds sampled from submanifolds embedded in Euclidean spaces. The second fundamental form is a useful tool to study Riemannian manifold extrinsically. On the one hand, since there is no prior information about the manifold, it prohibits the intrinsic way to study the manifold from sampling data. On the other hand, the second fundamental form is closely related to other geometric quantities of the manifold. For example, in [30] the authors showed that the operator norm controls the reach of the manifold. In [2] they proposed a simple way to compute the Riemannian hessian based on the Weingarten map. Furthermore, it is well known that the second fundamental form measures how manifolds curve in the ambient space. Once the second fundamental form is known, we can compute all kinds of curvature from it.

Our motivation comes from the need to estimate the curvature for unstructured point cloud data. Efficient estimation of (Gaussian/mean/principal) curvature for point cloud data is an important but difficult problem. In geometry, curvature contains much information of the underlying space of an unordered point set. Therefore it provides prior information in many applications such as surface segmentation [43, 33], surface reconstruction [42, 8, 39], shape inference [38, 12], point cloud simplification [19, 32] and feature extraction [17]. However, methods used to estimate curvature are restricted. Direct computation on point clouds often requires local parametrization. One needs to fit a local parametric surface first. Then the curvature is obtained by substituting the coefficients into an analytical formula. Another way is to estimate curvature after surface reconstruction which turns point clouds into triangular meshes or level sets of some distance function [35, 16, 22, 27]. However, there is little theory about estimation error analysis like convergence rate. This is not surprising since these existing algorithms are not aiming at minimizing the estimation error of the curvature, but instead to minimize the error of the surface approximation. In addition, even when surfaces are close to each other under the Euclidean distance, their curvature might not be close to each other. For example, we can perturb a straight line a little bit so that the curvature is far away from zero. As a result, a more direct and efficient approach with theoretical guarantee on the estimation error is needed.

Recently several methods have been proposed to estimate the second fundamental form. In [11], the authors proposed the notions of generalized second fundamental form and mean curvature based on the theory of varifolds. The generalized mean curvature descends to the classical mean curvature if the varifold is the standard measure on a submanifold. Under specific conditions, they proved that the optimal convergence rate is $O(n^{-1/2})$ where $n$ is the sample size, regardless of the dimension of the varifold. In [1], the authors used polynomial fitting to estimate the tangent space, the second fundamental form and the support of manifold simultaneously. Under some assumptions for the regularity of the manifold and sampling, they proved the convergence rate for the second fundamental form is $O(n^{-(k-2)/m})$ where $k$ is regularity parameter and $m$ is the dimension of manifold. In computer vision, several approaches are based on triangular meshes rather than point clouds. For example, S. Rusinkiewicz approximates the second fundamental form per triangle for meshes [34]. J. Berkmann and T. Caelli proposed two covariance matrices to approximate the second fundamental form [9]. For curvature estimation, there are methods proposed without the need of estimation of second fundamental form. In [28], the authors introduced the definition of Voronoi covariance measure (VCM) to estimate the curvature for noisy point clouds. Though a stability theorem is proved, there is no assertion on convergence rate. In the classic paper [40], G. Taubin defined a matrix by an integration formula. As illustrated in [20], this is nothing but the Weingarten map. This formula is adopted by [20] to estimate the principal curvature on point sets. The authors then proposed a method for anisotropic fairing of a point sampled surface using mean curvature flow.

1.1 Our Contribution

In this paper, we proposed a two-step procedure to estimate the Weingarten map for point clouds sampled from submanifolds embedded in Euclidean spaces. Firstly, we use local PCA to estimate the tangent spaces and normal spaces. Secondly, we apply a least-square method to compute the matrix representing the Weingarten map under the estimated tangent basis. The algorithm is general for point clouds in any dimension, and is efficient to implement due to low complexity.

A statistical model is set up to analyze the convergence rate for the Weingarten map estimator. Under the assumption of exact tangent spaces and given normal vector field, we proved that if the bandwidth is chosen to be $O(n^{-1/(m+4)})$, then the optimal convergence rate will be $O(n^{-4/(m+4)})$. Other than kernel method, we also discussed the $k$-nearest-neighbor method, which is more convenient to use in practice. Compared with the method proposed in [11], our method converges faster in low dimension. In comparison with the estimator proposed in [1], our method gives a closed form and easier to compute.

The convergence rate is verified by numerical experiments on two synthetic data sets. We also compare WME with the traditional quadratic fitting method, the state-of-art algorithm in this literature. Our method yields better results than the quadratic fitting in both MSE and robustness. As an application, we propose a curvature-based clustering method in point cloud simplification. Furthermore, we reconstruct surfaces based on the simplified point clouds to give a visible comparison. Three real data sets are tested to show the gain of our WME algorithm.

Outlines

This paper is organized as follows. We introduce our WME method in section 2 followed by a statistical model to analyze the convergence rate in section 3. In section 4, we verify the convergence rate and compare our method with quadratic surface fitting method using synthetic data. Applications to brain cortical data and experiments on point cloud simplification are given in section 5. In section 6, we discuss the possible applications of WME algorithm in future works.

2 Algorithm Description

Let ${\mathcal{M}}\subseteq\mathbb{E}^{d}$ be an $m$-dimensional submanifold in a $d$-dimensional Euclidean space with induced Riemannian metric. At each point $p$ we have the following decomposition

[TABLE]

Let $\overline{\nabla}$ be the standard connection in $\mathbb{E}^{d}$. For any normal vector $\xi\in T_{p}^{\perp}{\mathcal{M}}$, extend $\xi$ to be a normal vector field $\tilde{\xi}$ on ${\mathcal{M}}$. The Weingarten map or the shape operator at $p$ with respect to $\xi$ is defined as

[TABLE]

where ${}^{\top}:E^{d}\to T_{p}{\mathcal{M}}$ is the orthogonal projection to the tangent space. It can be verified that the definition of $A_{\xi}$ is independent of the extension of $\xi$ (see appendix A).

Example 1**.**

Let ${\mathcal{H}}^{d-1}\subseteq\mathbb{E}^{d}$ be a hypersurface, and $\tilde{\xi}$ be a unit normal vector field on ${\mathcal{H}}^{d-1}$. The Gauss map $g:{\mathcal{H}}^{d-1}\to\mathbb{S}^{d-1}$ sending any point on the hypersurface to a point on the unit sphere is defined by $g(p)=\tilde{\xi}_{p}$. For any $X\in T_{p}{\mathcal{H}}^{d-1}$, we have

[TABLE]

that is, $-A_{\tilde{\xi}_{p}}$ is the tangent map of the Gauss map.

The Weingarten map measures the variation of the normal vector field. In fact, from the following ‘Taylor expansion’ of a normal vector field, we can see that the Weingarten map plays the role similar to the derivative of a function.

Proposition 1**.**

Let $\tilde{\xi}$ be a normal vector field on the submanifold ${\mathcal{M}}\subseteq\mathbb{R}^{d}$. Suppose that $p$ is a point on ${\mathcal{M}}$ and $q$ is any point within the geodesic neighborhood of $p$. Denote ⊤ to be the orthogonal projection to the tangent space at $p$. We have

[TABLE]

Proof.

Let $\mathbf{r}:U^{m}\subseteq\mathbb{R}^{m}\to{\mathcal{M}}$ be the exponential map such that $\mathbf{r}(0)=p$. Denote $\mathbf{u}=(u^{1},\cdots,u^{m})\in U^{m}$. Then the vector fields $\{\frac{\partial\mathbf{r}}{\partial u^{i}}\}_{i=1}^{m}$ form a local tangent frame. Denote $\tilde{\xi}(\mathbf{r}(\mathbf{u}))$ by $\tilde{\xi}(\mathbf{u})$. Consider the following expansions

[TABLE]

By definition of $A_{\tilde{\xi}_{p}}$, we have

[TABLE]

Substituting to (2.5), we obtain that

[TABLE]

Note that $\|\mathbf{u}\|$ represents the geodesic distance on the manifold by the property of the exponential map. According to the proposition 2 in appendix B, the geodesic distance can be approximated by the Euclidean distance in the same order. Therefore, (2.4) follows. ∎

Assume that $n$ points $x_{1},\cdots,x_{n}$ viewed as points in $\mathbb{R}^{d}$ are independently sampled from some distribution on ${\mathcal{M}}$. The object is to estimate the Weingarten maps at each $x_{i}$ ($i=1,2,\cdots,n$). Thus, firstly, we need to estimate the tangent and normal spaces. Otherwise, the estimation of Weingarten maps does not make sense without specifying the normal directions. Local PCA is an extensively used method in manifold learning to estimate tangent and normal spaces, whose effectiveness and consistency is well understood (see [36]). Secondly, after tangent space estimation, proposition 1 indicates a simple linear model to estimate the Weingarten maps from the data points, which can be resolved by a least-square method. Therefore, our WME algorithm is a two-step procedure where each step can be proceeded efficiently.

Local PCA: For every data point $x_{i}$ ($i=1,2,\cdots,n$) we try to estimate a basis $e_{i}^{1},\cdots,e_{i}^{m}$ to the tangent space and a basis $\xi_{i}^{1},\cdots,\xi_{i}^{d-m}$ to the normal space. Fix a parameter $h_{\text{PCA}}>0$ and define the following set$I_{i}=\{j\in\mathbb{N}|\|x_{j}-x_{i}\|\leq h_{\text{PCA}}\}$. The local covariance matrix is defined as

[TABLE]

where each data point is regarded as a column vector in $\mathbb{R}^{d}$ and $\bar{x}_{i}=\frac{1}{|I_{i}|}\sum_{j\in I_{i}}x_{j}$ is the mean of neighboring points. The first $m$ eigenvectors of Cov corresponding to the first $m$ largest eigenvalues constitute the basis to the tangent space at $x_{i}$, while the last $d-m$ eigenvectors corresponding to the last $d-m$ smallest eigenvalues form the basis to the normal space at $x_{i}$.

Normal Vector Extension: Suppose that we pick up the normal vector $\xi_{i}^{\alpha}$ ($i=1,2,\cdots,m$, $\alpha=1,2,\cdots,d-m$) and want to estimate the Weingarten map with respect to $\xi_{i}^{\alpha}$. Extend $\xi_{i}^{\alpha}$ to be a normal vector field $\tilde{\xi}^{\alpha}$ on ${\mathcal{M}}$ by setting

[TABLE]

for $j=1,2,\cdots,n$. That is, we project $\xi_{i}^{\alpha}$ to the normal space at $x_{j}$. Since projection varies smoothly with respect to the points on the manifold, it results in a smooth normal vector field.

Assume $K:\mathbb{R}\to\mathbb{R}$ is a twice differentiable function supported on $[0,1]$. For example, the truncated Gaussian kernel on $[0,1]$ or Epanechnikov kernel. Let $y\in\mathbb{R}^{d}$. Given $h>0$, define $K_{h}:\mathbb{R}^{d}\to\mathbb{R}$ to be

[TABLE]

Let $E_{i}=[e_{i}^{1},\cdots,e_{i}^{m}]$ be the matrix consisting of the basis to the tangent space at $x_{i}$. According to proposition 1, if $x_{j}$ is close to $x_{i}$,

[TABLE]

where $A_{\xi_{i}^{\alpha}}$ is understood as $m\times m$ matrix. Therefore, we want to find a matrix $\widetilde{A}_{\xi_{i}^{\alpha}}$ which minimizes the following residual

[TABLE]

Set

[TABLE]

Then the solution of (2.12) is given in the following closed form

[TABLE]

Remark 1**.**

The above method gives the estimation of the Weingarten maps at each point with respect to all normal basis. Since each Weingarten map is an $m\times m$ matrix and the normal basis consists of $d-m$ vectors, the whole procedure gives $m\times m\times(d-m)$ coefficients (in fact $\frac{m(m+1)}{2}\times(d-m)$ coefficients since the Weingarten map is symmetric) at each point. Using these coefficients we can give the estimation of second fundamental form and mean curvature and sectional curvature at each point. If we assume the underlying manifold is of low dimension, then all the coefficients are about the size $O(dn)$ where $d$ is the dimension of the ambient space and $n$ is the size of the point cloud.

Remark 2**.**

Since the Weingarten map is a self-adjoint operator on the tangent space, the matrix $\widetilde{A}_{\xi_{i}^{\alpha}}$ should be symmetric. It is natural to solve (2.12) on the space of symmetric matrices. However, as we will prove later, the solution of (2.12) converges to the true matrix. It is not necessary to solve a more complex optimization problem on the space of symmetric matrices. In cases where symmetry is important we can always use the symmetrization $\frac{1}{2}(\widetilde{A}_{\xi_{i}^{\alpha}}^{\rm t}+\widetilde{A}_{\xi_{i}^{\alpha}})$.

3 Convergence Rate

3.1 Statistical Modeling

Let $\xi$ be a normal vector field on ${\mathcal{M}}$, and $P$ be a random vector valued in ${\mathcal{M}}$ with smooth positive density function $f$. Fix a point $p\in{\mathcal{M}}$. To avoid complexity on notations, we drop the subscripts and $A:T_{p}\mathcal{M}\to T_{p}\mathcal{M}$ is always understood as the Weingarten map (or its matrix representation if a basis is specified) at $p$ with respect to $\xi_{p}$. We rewrite proposition 1 as follows: when $P$ is within the normal neighborhood of $p$, we have

[TABLE]

where $\eta:{\mathcal{M}}\to\mathbb{R}^{m}$ is assumed to be a bounded smooth function in the neighborhood of $p$. In addition, assume that a basis $e_{1},e_{2},\cdots,e_{m}$ for the tangent space $T_{p}\mathcal{M}$ is given. Set $\mathbf{\Xi}=(\xi_{P}-\xi_{p})^{\top}$, $\mathbf{\Delta}=(P-p)^{\top}$ to be the vectors representing the coordinates under the basis, and set $\mathbf{K}_{P}=K_{h}(P-p)$. Consider the following optimization problem

[TABLE]

That is, we want to find the minimizer of the function

[TABLE]

By setting $dF/d\mathbf{A}=0$, the population solution can be given in the following closed form

[TABLE]

where we have set $\mathbf{L}=\mathbb{E}_{f}[\mathbf{\Xi}\mathbf{\Delta}^{\textrm{t}}\mathbf{K}_{P}]$ and $\mathbf{D}=\mathbb{E}_{f}[\mathbf{\Delta}\mathbf{\Delta}^{\textrm{t}}\mathbf{K}_{P}]$. Denote the coordinate components $\mathbf{\Xi}\cdot e_{j}$ by $(\mathbf{\Xi})_{j}$ and $\mathbf{\Delta}\cdot e_{j}$ by $(\mathbf{\Delta})_{j}$ for $j=1,2,\cdots,m$. In matrix form we have

[TABLE]

Let $x_{1},\cdots,x_{n}$ be i.i.d. samples from $P$. Let $\Xi_{i},\Delta_{i}$ be the quantities obtained by replacing the random vector in $\mathbf{\Xi},\mathbf{\Delta}$ with samples for $i=1,2,\cdots,n$. Substituting the expectation by empirical mean, we obtain the empirical solution

[TABLE]

where

[TABLE]

It is easy to check that the empirical solution given here is the same as the one given in (2.14). That is, it is the solution of the following empirical optimization problem

[TABLE]

Finally, the mean square error (MSE) is defined as

[TABLE]

where $\|\cdot\|_{F}$ denotes the Frobenius norm.

3.2 Mean Square Error

In the estimation of either variance or bias, the main obstacle is how to control the norm of inverse matrices $\mathbf{D}^{-1}$ and $\widetilde{D}^{-1}$. The method is somewhat similar to that in kernel density estimation. When the bandwidth $h$ is small, the integration is taken over a normal neighborhood near $p$. Using the parametrization of the exponential map, the integral domain is the tangent space and we can use Taylor expansion to find the leading terms. Before that we notice the following fact concerning Euclidean distance and manifold distance, which can be found in [31].

Lemma 1**.**

Let $d_{p}$ and $\rho_{p}$ be the Riemannian distance and Euclidean distance to $p$, respectively. There exists a function $R_{p}(h)$ and positive constants $\delta_{R_{p}}$, $C_{R_{p}}$ such that when $h<\delta_{R_{p}}$, $d_{p}(y)\leq R_{p}(h)$ for all $y$ with $\rho_{p}(y)\leq h$, and furthermore,

[TABLE]

Lemma 1 indicates that for $h$ small enough if $\|P-p\|\leq h$ then there exists a function such that the geodesic distance between $P$ and $p$ is controlled by $R_{p}(h)$. First we give the estimation of $\mathbf{D}^{-1}$. Some properties of exponential map are presented in appendix B.

Lemma 2**.**

There exists a positive constant $h_{0}$ such that when $h<h_{0}$,

[TABLE]

where $I$ is the identity matrix. Thus, the inverse is given by

[TABLE]

Proof.

Pick $h_{0}$ such that for $h<h_{0}$ the exponential map is well defined on the geodesic ball of radius $R_{p}(h)$. For an arbitrary element $\mathbf{D}_{kl}$, let $\mathbf{r}$ be the exponential map. We have

[TABLE]

where $g(\mathbf{u})$ denotes the Riemannian metric matrix. In normal coordinates, $g_{ij}(\mathbf{u})=\delta_{ij}+o(\|\mathbf{u}\|)$. Changing the variable of integration to $\mathbf{u}=h\mathbf{z}$, we obtain that

[TABLE]

The integral domain can be divided into two parts: the first is the unit ball $Q_{1}=\{\mathbf{z}|\|\mathbf{z}\|\leq 1\}$ and the second is $Q_{2}=\{\mathbf{z}|1\leq\|\mathbf{z}\|\leq R_{x}(h)/h\}$. On $Q_{1}$, the integration is

[TABLE]

where we have used the symmetry of integration and $\delta_{kl}$ denotes the Kronecker delta. On $Q_{2}$, note that the integrand is a quantity of $O(h^{2})$ whereas the area of integral domain is $O(h^{2})$ by lemma 1. Overall, we have shown that

[TABLE]

which proves equation (3.11). The inverse matrix is given by the identity $(I-D)^{-1}=\sum_{k=0}^{\infty}D^{k}$ for small $D$. ∎

Using the presented results, we are able to give the convergence order of bias.

Lemma 3**.**

As $h\to 0$, $\|\mathbf{A}-A\|_{F}^{2}=O(h^{4})$.

Proof.

According to the model assumption (3.1), if we set $\mathbf{\Theta}=\mathbf{\Xi}-A\mathbf{\Delta}=\eta(P)\|P-p\|^{2}$, then we have

[TABLE]

It suffices to estimate $\mathbb{E}_{f}[\mathbf{\Theta}\mathbf{\Delta}^{\textrm{t}}\mathbf{K}_{P}]$. As in the proof of lemma 2, for $h<h_{0}$, the integration is

[TABLE]

After the changing of variable $\mathbf{u}=h\mathbf{z}$, note that by symmetry we have

[TABLE]

Thus the coefficient of $h^{3}$ vanishes. Set

[TABLE]

The leading term is

[TABLE]

Utilizing the estimation for $\mathbf{D}$, we obtain that

[TABLE]

∎

For the matrix $\widetilde{D}$ we have error coming from random sampling. Therefore, the estimation only holds with high probability (w.h.p.). That is, for any $\epsilon>0$ the estimation holds with probability greater than $1-\epsilon$. For the simplicity of statements, we omit the specific computation of quantities involving $\epsilon$.

Lemma 4**.**

Suppose that when $h\to 0$ and $nh^{m}\to\infty$ as $n\to\infty$, w.h.p. the following equality holds

[TABLE]

Thus the inverse is given by

[TABLE]

Proof.

For an arbitrary element $\widetilde{D}_{kl}$, note that

[TABLE]

The variance is

[TABLE]

Using the same method as in the proof in lemma 2, we can find the leading term is

[TABLE]

By Chebyshev’s inequality, for any $\epsilon>0$,

[TABLE]

as $h\to 0$ and $nh^{m}\to\infty$. Set the right side to be $O(1)$, we can show that w.h.p. $|\widetilde{D}_{kl}-\mathbf{D}_{kl}|\leq h^{2}/\sqrt{nh^{m}}$, which proves the equation (3.23). The inverse is given by the following identity

[TABLE]

∎

Using the presented results, now we can estimate the variance.

Lemma 5**.**

Suppose that when $h\to 0$ and $nh^{m}\to\infty$ as $n\to\infty$, the variance is $O(\frac{1}{nh^{m}})$ w.h.p.

Proof.

First by interpolating a mixing term we have

[TABLE]

For the first term we have

[TABLE]

For an arbitrary term in $\mathbf{L}-\widetilde{L}$ we have

[TABLE]

When $h<h_{0}$, by the model assumption (3.1), the above quantity (3.32) is bounded by

[TABLE]

where the leading term is

[TABLE]

Together with the estimation for $\mathbf{D}^{-1}$, we see that $(*)=O(\frac{1}{nh^{m}})$.

For the second term we have

[TABLE]

Now the order of each quantity in (3.35) is

• (1)

Similar to the proof of lemma 4, note that $\mathbb{E}_{f}[\widetilde{L}_{kl}]=\mathbf{L}_{kl}$ and the variance is controlled by a quantity of order $O(h^{4}/nh^{m})$ w.h.p. Thus,

[TABLE]

• (2)

Using the same method as in the estimation for $\|\mathbf{L}-\widetilde{L}\|_{F}^{2}$, we see that it is controlled by $O(1/nh^{m-4})$ w.h.p.

• (3)

By lemma 4, this term is $O(1/h^{4})$ w.h.p.

• (4)

By lemma 2, this term is $O(1/h^{4})$.

Hence, the order of $(**)$ is $O(\frac{1}{nh^{m}})$. Overall, the rate for variance is proved. ∎

Altogether we can give the convergence rate of MSE.

Theorem 1**.**

Let $\xi$ be a normal vector field on ${\mathcal{M}}$, and $P$ be a random vector valued in ${\mathcal{M}}$ with smooth positive density function $f$. Assume that $K:\mathbb{R}\to\mathbb{R}$ is a twice differentiable function supported on $[0,1]$. In addition, assume that a basis $e_{1},e_{2},\cdots,e_{m}$ for the tangent space $T_{p}\mathcal{M}$ is given. Let $x_{1},\cdots,x_{n}$ be i.i.d. samples from $P$. When $h\to 0$ and $nh^{m}\to\infty$ as $n\to\infty$, the mean square error defined as in (3.9) is

[TABLE]

If $h$ is chosen to be proportional to $n^{-1/(m+4)}$, the optimal convergence rate is given by $O(n^{-4/(m+4)})$.

3.3 $k$-Nearest-Neighbor Method

$k$-nearest-neighbor method is widely used in many settings since it is intuitively simple, easy to implement and computationally efficient. For simplicity, we focus on disk neighbors with fixed radius $h$ in this article. However, the disk neighbor can be replaced by $k$-nearest-neighbor without any other changes. According to Theorem 1 and the well known relation $\frac{k}{n}\sim h^{m}$ (citation), we suggest to set $k=O(n^{4/(m+4)})$ in practice to reach the optimal convergence rate $O(n^{-4/(m+4)})$.

3.4 Comparison with Other Estimators

In [11], the authors proposed generalized second fundamental form and mean curvature for varifolds based on geometric measure theory. Therefore, even if the underlying space has singularities (for example, self-intersecting points), they are able to compute the generalized mean curvature. Basically, a $m$-varifold is a Radon measure on the space $\mathbb{R}^{d}\times\mathbb{G}_{m,d}$ where $\mathbb{G}_{m,d}$ is the Grassmannian manifold of $m$-planes in $\mathbb{R}^{d}$. For an $m$-varifold $V$, they define the generalized mean curvature field using the first variation

[TABLE]

The generalized mean curvature descends to the classical mean curvature in the sense that $\delta V=-H\mathcal{H}^{m}_{|{\mathcal{M}}}$ if $V$ is the standard measure on ${\mathcal{M}}$. However, in general, the first variation cannot be given in closed form. For a point cloud varifold $V=\sum_{j}=1^{n}m_{j}\delta_{x_{j}}\otimes\delta_{P_{j}}$, the first variation $\delta V$ is not a measure. Therefore, they put forward the following quantity to approximate the mean curvature for point clouds

[TABLE]

where $\alpha,\beta$ are functions supported on $[-1,1]$ and such that $\int_{\mathbb{R}^{d}}\alpha(x){\rm d}x=\int_{\mathbb{R}^{d}}\beta(x){\rm d}x=1$, and $C_{\alpha},C_{\beta}$ are constants related to $\alpha,\beta$. They proved that $H^{V}_{\alpha,\beta,\epsilon}$ converges to $H$ under certain conditions (Theorem 3.6 in [11]). Specifically, if the points are uniformly sampled and the tangents are exact, the convergence rate is $\frac{1}{n\epsilon}+\epsilon$ given that $n\epsilon\to\infty$. Therefore, the optimal bandwidth is $\epsilon=n^{-1/2}$ and the optimal convergence rate is also $n^{-1/2}$. Compared with WME algorithm, this method yields faster convergence when the underlying manifold is of relatively high dimension. However, for curves and surfaces, our algorithm provides faster convergence.

In [1], the authors proposed estimators for the tangent space, the second fundamental form and the support of manifold based on local polynomial fitting. Finite sample rates are derived. Let ${\mathcal{M}}$ be a submanifold subject to some regularity conditions, and in addition assume that the data are sampled uniformly on the manifold. Then the convergence rate is $O(n^{-(k-2)/m})$ where $k$ is regularity parameter and $m$ is the dimension of manifold. However, the estimator cannot be given in a closed form. Nonconvex optimization techniques should be involved in order to approximate the estimators.

4 Numerical Experiments

We apply WME algorithm on various synthetic data sets and verify the optimal convergence rate proved in theorem 1. Both kernel method and $k$-Nearest-Neighbor method are used to show the consistency of WME algorithm. Then the algorithm is applied to curvature estimation. A comparison with classical local quadratic fitting method is carried out to demonstrate the efficiency and robustness.

4.1 Kernel Method

We verify the convergence rates on three synthetic data sets:

**Conical Spiral: ** A conical spiral is a space curve given by

[TABLE]

where $r,a>0$ are parameters. Here we set $r=1$ and $a=1.5$. It is a 1-dimensional manifold embedded in 3-Euclidean space. If we choose $h$ proportional to $n^{-1/5}$, then the optimal convergence rate will be $O(n^{-4/5})$.

**Torus: ** A 2-dimensional torus is given by

[TABLE]

where $R>2r>0$ are parameters. Here we set $R=4$ and $r=0.5$. Thus if we choose $h$ proportional to $n^{-1/6}$, the optimal rate will be $O(n^{-2/3})$.

**Ellipsoid: ** A 3-dimensional ellipsoid is given by

[TABLE]

where $a,b,c,d>0$ are parameters. Here we set $a=1$, $b=1.15$, $c=1.31$ and $d=1.7$. If the bandwidth is chosen proportional to $n^{-1/7}$, the convergence rate will be $O(n^{-4/7})$.

We uniformly sample $1,000\sim 20,000$ points using the parametrizations given as above. The bandwidth is chosen so that the convergence rate will be optimal. We perform leave-one-out estimates on the data points. The $\log(\text{MSE})$ is drawn with respect to $\log(\text{number of points})$. Therefore the slope corresponds to the order of convergence rate. The results are shown in figure 1.

4.2 $k$-Nearest-Neighbor Method

We test the optimal convergence rate on a 2-dimensional torus and a 2-dimensional ellipsoid. The sample size for the torus with major radius $5$ and minor radius $2$ ranges from $1000$ to $20000$. The sample size for the ellipsoid with length of principal axes $6,6,8$ ranges from $7000$ to $30000$. $k$ is chosen to be $n^{\frac{2}{3}}$. Points are uniformly sampled. Figure 2 is the $\log(n)$-$\log(\text{MSE})$ plot. The slopes match the optimal convergence rate.

4.3 Curvature Estimation

We compare our method with a local quadratic surface fitting method. This is chosen for two reasons. On one hand, quadratic fitting is a commonly used method. Other complicated fitting algorithms involve extra scaling parameters which are difficult to tune in practice. On the other hand, quadratic fitting is studied by many scientists. In [26], the authors compared five methods in computing Gaussian and mean curvature for triangular meshes of 2 dimensional surfaces. The result turns out that quadratic fitting exceeds other methods in computing mean curvature.

The method of quadratic surface fitting is illustrated as follows. First translate and rotate the $k$-nearest neighbors of a point so that its normal vector coincides with $z$-axis. Then fit the paraboloid $z=ax^{2}+bxy+cy^{2}$ by least square. The Gaussian curvature and mean curvature at origin $P$ are given by

[TABLE]

The MSE of Gaussian curvature and mean curvature are compared as follows. We sample $1000\sim 20000$ points on a torus with major radius 5 and minor radius 2. The number of $k$-nearest neighbors is set to be 100 for each iteration. The result in Figure 2 shows that our method excels the quadratic fitting method without introducing any computational complexity.

The robustness is compared as follows. Again we sample $1000,2000,\cdots,20000$ points on the same torus with multivariate Gaussian noise with zero mean and covariance $\sigma^{2}I_{3}$ where $\sigma^{2}\in\{0.01,0.05,0.001\}$. The MSE of Gaussian curvature and mean curvature for different noises in Figure 4 show that our method is more robust.

5 Applications

We apply our WME algorithm to real data sets. The first application is curvature estimation for brain cortical surfaces. We demonstrate the robustness of WME method on a real cortical surface data set. The second application is point cloud simplification which is a hot topic in computer vision. We propose a new method based on curvature to simplify large point cloud data sets. The results of following experiments show that our algorithm is also practical for real point clouds.

5.1 Brain Cortical Surface Data

To further illustrate the robustness, we test our method on the real brain cortical surface data. A point cloud of left brain cortical surface is obtained from Human Connectome Projects (s3://hcp-openaccess/HCP_1200/100206/MNINonLinear/Native/), consisting of 166,737 position vectors. This data is noisy and there is no information about the true curvature of the surface so there is no ground truth and the error can’t be calculated. Instead we propose an indirect way to evaluate the performance. Firstly, we estimate the Gaussian and mean curvature for the point cloud based on the entire dataset. The results are regarded as the true curvature for the underlying cortical surface. Then, the data is divided into training and testing sets. We recalculate the Gaussian and mean curvature for training data. For each testing data, the curvature is inferred to be the mean of curvature for its $k$ nearest neighbors in training data. Finally, we compute the mean square error of the curvature for testing data. The same procedure is also carried out using quadratic surface fitting method. From Figure 5, the mean square error obtained from WME is monotonically decreasing as the number of testing data increases but the error from quadratic surface fitting method is fluctuating, which means that WME is more robust on this real and complicated dataset.

5.2 Point Cloud Simplification

Point clouds are often converted into a continuous surface representation such as polygonal meshes and splines. This process is called surface reconstruction [8]. The reconstruction algorithms require large amounts of memory and do not scale well with data size. Hence before further processed, the complexity of point cloud data should be reduced first. In [32], the authors proposed three types of simplification algorithms: clustering methods, iterative simplification and particle simulation. These methods are based on a quantity defined by the covariance of local data. As claimed by [32], this quantity reflects the curving of point cloud. However, the clear relation between this quantity and the curvatures needs to be further studied. Here we propose a curvature-adaptive clustering simplification algorithm and compare with uniform clustering simplification algorithm.

The uniform clustering method is described as follows. Starting from a random seed point, a cluster $C_{0}$ is built by successively adding the nearest neighbors. The process is terminated when the size of clusters reaches the previously set threshold. The next cluster $C_{1}$ is built in the same procedure with all points in $C_{0}$ excluded. Each cluster is represented by its mean as a representative. The simplified point cloud is given by the representatives.

Intuitively, to preserve the geometric details of point cloud, the points in highly curved region should be kept. Therefore, a seed point with larger curvature should generate smaller cluster. Let $\Omega$ represent any kind of (Gaussian, mean or principal) curvature. Suppose that $|\Omega|_{\max}$ is the largest absolute curvature of the entire surface. Starting from a random seed point $p$, with absolute curvature $|\Omega|_{p}$, a cluster $C_{p}$ is built by successively adding the nearest neighbors. The process is terminated when the size of cluster reaches

[TABLE]

where $0<c<1$ is the scaling constant and $T$ is the preset threshold. $\lceil\cdot\rceil$ denotes the ceiling function. The cluster and its curvature are represented by the mean of its points and mean of corresponding curvature. This yields a non-uniform clustering method.

The algorithms are applied to three scanned data sets: the Duke dragon, the Stanford bunny and the Armadillo. Here we adopt absolute mean curvature for curvature-based simplification. After simplification, we apply the Moving Least Square (MLS) method for surface reconstruction [8]. The visualized surfaces in Figure 6, 7 and 8 give a direct comparison of two algorithms. In each of these figures, the first subfigure represents the simplified point cloud using uniform clustering and the second subfigure is the reconstructed surface from the point cloud. Similarly, the third subfigure represents result using curvature-based clustering and the fourth subfigure shows the surface reconstructed from it. Results show that WME preserves more geometric information than the uniform method, especially for the region with larger curvature.

6 Discussions and Future Works

This paper introduced a new algorithm to estimate the Weingarten map for point clouds sampled from some distribution on a submanifold embedded in Euclidean space. A statistical model was also established to investigate the optimal convergence rate of the Weingarten map estimator. Numerical experiments were carried out to validate the consistency of the algorithm. Compared with other methods concerning the second fundamental form and curvature estimation, our method showed great robustness and efficiency. Applications to real data sets indicated that our method worked well in practice.

The Weingarten map is a vital tool which is indispensable in the study for submanifolds. In the literature of manifold learning, where data sets from high dimensional Euclidean space are always assumed to be sampled from some underlying low dimensional manifold, the Weingarten map helps researchers to investigate the underlying space from an extrinsic view. Furthermore, the Weingarten map is also in close relation with many other concepts in Riemannian geometry. Curvature is mentioned as one of them. Therefore, we expect many possible applications in the future. For example, the followings are of our concern.

Optimization on Manifolds. In [2] the authors presented a useful relationship between the Weingarten map and Riemannian hessian of functions, and the latter plays an important role in Riemannian optimization. Especially, Riemannian hessian is critical in all kinds of second-order optimization methods on Riemannian manifolds [10, 18]. We expect our method can be applied in Riemannian optimization problems such as low-rank matrix completion [41], high-dimensional tensor completion [37] and independent subspace analysis [29].

Laplacian-based Methods. Laplacian-based methods are extensively studied in manifold learning. Theories and applications about graph Laplacian are popular in the realm of computer vision [13], network analysis[15], and statistical learning [5, 7, 6]. It is known that the Weingarten map is related to the Beltrami-Laplacian operator on Riemannian manifolds. We expect to find other ways to construct the Laplacian estimator on point clouds.

Statistical Inference on Point Clouds. In the application for brain cortical data our method showed great efficiency in curvature estimation. In recent researches, many kinds of curvature appeared to be useful indexes in statistical inference and regression [25, 44]. We expect to find applications in fields such as biostatistics, medical imaging and neuroscience.

\appendixpage

Appendix A A The Second Fundamental Form, Weingarten Map and Curvature

Let ${\mathcal{M}}\subseteq\mathcal{N}$ be a submanifold with induced Riemannian metric. The connection in $\mathcal{N}$ is denoted by $\overline{\nabla}$ and the Riemannian connection on ${\mathcal{M}}$ is denoted by $\nabla$. Let $\mathfrak{X}({\mathcal{M}})$ be the set of vector fields on ${\mathcal{M}}$. For any vector fields $X,Y\in\mathfrak{X}({\mathcal{M}})$, we have the Gauss formula

[TABLE]

where

[TABLE]

The map defined by

[TABLE]

has the following properties:

[TABLE]

Therefore, $h$ is a second-order symmetric covariant tensor field on ${\mathcal{M}}$.

Definition 1**.**

The tensor field $h$ is called the second fundamental form of submanifold ${\mathcal{M}}$. Especially, at every point $p\in{\mathcal{M}}$, the second fundamental form defines a symmetric, bilinear map

[TABLE]

Let $\xi\in\Gamma(T^{\bot}{\mathcal{M}})$ be a normal vector field. We have the Weingarten formula

[TABLE]

where

[TABLE]

Theorem 2**.**

For every $\xi\in\Gamma(T^{\bot}{\mathcal{M}})$, the map defined by

[TABLE]

is a smooth $(1,1)$-form tensor field. At every point $p\in{\mathcal{M}}$, there is a self-adjoint transformation

[TABLE]

and satisfies the identity

[TABLE]

Therefore, at every point $p\in{\mathcal{M}}$, there is a bilinear map $A:T_{p}^{\bot}{\mathcal{M}}\times T_{p}{\mathcal{M}}\to T_{p}{\mathcal{M}}$ such that

[TABLE]

Definition 2**.**

At $p\in{\mathcal{M}},\xi\in T_{p}^{\bot}{\mathcal{M}}$, the linear map $A\xi:T_{p}{\mathcal{M}}\to T_{p}{\mathcal{M}}$ is called the shape operator or the Weingarten map at $p$ with respect to $\xi$.

Let $e_{1},\cdots,e_{m}$ be a basis of the tangent space and $\xi_{1},\cdots,\xi_{d-m}$ be a basis of the normal space. Assume

[TABLE]

Then

[TABLE]

Definition 3**.**

Let

[TABLE]

then $H$ is independent to the choice of orthonormal fields, and is called the mean curvature vector field.

According to the definition, we have

[TABLE]

where $H^{\xi}=\tilde{g}(H,\xi)$ is called the mean curvature along $\xi$. The value $\|H\|=(\sum_{\alpha}\|H^{\alpha}\|^{2})^{1/2}=\frac{1}{m}(\sum_{\alpha=1}^{d-m}\text{trace}(A_{\xi_{\alpha}})^{2})^{1/2}$ is called the mean curvature. If $M$ is a hypersurface, $\|H\|=\frac{1}{m}|\text{trace}(A_{\xi_{\alpha}})|$. For a surface in 3-dimensional Euclidean space, this coincides with the absolute mean curvature.

Let $R$ and $\overline{R}$ be the Riemann curvature tensor of ${\mathcal{M}}$ and $\mathcal{N}$ respectively. From the equations (A.1) and (A.6) we can derive the Gauss equation

[TABLE]

Especially the sectional curvature can be expressed as

[TABLE]

If the ambient space is Euclidean, then $\overline{R}$ vanishes identically. Let $e_{i},e_{j}$ span a two-plane $\pi_{ij}$, then the sectional curvature of $\pi_{ij}$ is

[TABLE]

From $A_{\xi_{\alpha}}$ we extract the $2\times 2$ submatrix $A_{\xi_{\alpha}}|_{\pi_{ij}}$ with $i$th and $j$th row and column. Then

[TABLE]

Example 2** (Hypersurfaces).**

If ${\mathcal{M}}$ is a hypersurface of $\mathbb{E}^{d}$ with a unit normal vector field $\xi$, the sectional curvature is $K(\pi_{ij})=\det(A_{\xi}|_{\pi_{ij}})$. If $e_{1},\cdots,e_{d-1}$ is a basis that diagonalizes the Weingarten map with eigenvalues $\lambda_{1},\cdots,\lambda_{d-1}$, then $K(\pi_{ij})=\lambda_{i}\lambda_{j}$. The eigenvalues are called principal curvature and eigenvectors are called principal directions. The determinant of $A_{\xi}$ is called the Gauss-Kronecker curvature.

Example 3** (Planar Curves).**

A planar curve is a 1-dimensional manifold embedded in the 2-plane. Let $\mathbf{t}$ be the tangent vector field and $\mathbf{n}$ be the normal vector field. We have

[TABLE]

where $\kappa$ is the curve curvature. Thus the estimation can be simply given by

[TABLE]

Example 4** (Space Curves).**

A space curve is a 1-dimensional manifold embedded in 3-space. Let $\mathbf{t}$ be the tangent vector field, $\mathbf{n}$ be the normal vector field, and $\mathbf{b}$ be the binormal vector field. We have the Frenet formula

[TABLE]

where $\kappa$ is curvature and $\tau$ is torsion. Let $\xi=\cos(\theta)\mathbf{n}+\sin(\theta)\mathbf{b}$ be a unit normal vector field. We have

[TABLE]

Similarly, let $\xi_{\perp}=-\sin(\theta)\mathbf{n}+\cos(\theta)\mathbf{b}$ be the unit normal vector perpendicular to $\xi$. Then

[TABLE]

By definition the mean vector field is $H=\cos(\theta)\kappa\xi-\sin(\theta)\kappa\xi_{\perp}=\kappa\mathbf{n}$, and the mean curvature is $\|H\|=\kappa$. As a byproduct we obtain the normal vector $\mathbf{n}=H/\|H\|$.

Appendix B B The Exponential Map

The second fundamental form is closely related to the injective radius $\iota({\mathcal{M}})$, as proved in [3]: If the operator norm $\|h\|\leq C$, then $\iota({\mathcal{M}})\geq\pi/C$. Thus the exponential map can be defined on the ball of radius $\pi/C$ in $T_{p}\mathcal{M}$ for all $p$.

The exponential map for submanifolds embedded in Euclidean spaces can be expressed in special forms. Fix $p\in\mathcal{M}$. Without loss of generality we transform $\mathcal{M}$ isometrically in $\mathbb{R}^{d}$ so that $p$ is the origin and the first $m$ standard basis $\{e_{1},\cdots,e_{m}\}$ spans the tangent space $T_{p}\mathcal{M}$. Consider the coordinate components of the exponential map

[TABLE]

We can express the derivatives of component functions using the geometric terms defined on $\mathcal{M}$. In fact, we have the following proposition.

Proposition 2**.**

Let $\mathbf{u}=(u^{1},\cdots,u^{m},0,\cdots,0)\in T_{p}{\mathcal{M}}\subseteq\mathbb{R}^{d}$. The exponential map for ${\mathcal{M}}$ can be expressed by

[TABLE]

Furthermore, the derivatives of component functions can be expressed by

[TABLE]

for $i,j,k,l=1,\cdots,m;\alpha=m+1,\cdots,n$.

Proof.

Firstly, since the tangent map of $\mathbf{r}$ at origin is the identity, we have

[TABLE]

We obtain that $\partial_{i}r^{j}(0)=\delta^{j}_{i}$ for $i,j=1,\cdots,m$ and $\partial_{i}r^{\alpha}(0)=0$ for $\alpha=m+1,\cdots,n$.

For the second derivatives, let $\mathbf{u}=\sum_{i=1}^{m}u^{i}e_{i}$ and $\gamma(s)=\mathbf{r}(s\mathbf{u})$. we have

[TABLE]

Thus, $\partial_{i}\partial_{j}r^{k}(0)=0$ for $i,j,k=1,\cdots,m$. Let $h_{\alpha}(\mathbf{u},\mathbf{v})=\langle h(\mathbf{u},\mathbf{v}),e_{\alpha}\rangle=\langle A_{e_{\alpha}}(\mathbf{u}),\mathbf{v}\rangle$. We see that \big{[}\partial_{i}\partial_{j}r^{\alpha}(0)=\langle A_{e_{\alpha}}e_{i},e_{j}\rangle\big{]}_{i,j} is the matrix representation for the Weingarten map $A_{e_{\alpha}}$ for $\alpha=m+1,\cdots,n$.

By (A.6), we have

[TABLE]

If we define $\nabla h(\mathbf{u},\mathbf{v},\mathbf{w})=\nabla_{\mathbf{w}}^{\perp}h(\mathbf{u},\mathbf{v})-h(\nabla_{\mathbf{w}}\mathbf{u},\mathbf{v})-h(\mathbf{u},\nabla_{\mathbf{w}}\mathbf{v})$, then we have

[TABLE]

Comparing the components on both sides we obtain that

[TABLE]

for $i,j,k,l=1,\cdots,m$ and $\alpha=m+1,\cdots,n$. We have deduced the Taylor expansion for the exponential map. ∎

Let $q=\mathbf{r}(s\mathbf{u})$ where $\|\mathbf{u}\|=1$. We can compare the Euclidean distance $\|q-p\|$ and the Riemannian distance $s=d_{\mathcal{M}}(p,q)$. In fact, from the Taylor expansion of the exponential map, we have

[TABLE]

Therefore, we have

[TABLE]