arXiv:1804.02811v1 [math.ST] 9 Apr 2018

CONNECTING DOTS – FROM LOCAL COVARIANCE TO
EMPIRICAL INTRINSIC GEOMETRY AND LOCALLY LINEAR
EMBEDDING
JOHN MALIK, CHAO SHEN, HAU-TIENG WU, AND NAN WU

Abstract. Local covariance structure under the manifold setup has been
widely applied in the machine learning society. Based on the established theoretical results, we provide an extensive study of two relevant manifold learning
algorithms, empirical intrinsic geometry (EIG) and the locally linear embedding (LLE) under the manifold setup. Particularly, we show that without an
accurate dimension estimation, the geodesic distance estimation by EIG might
be corrupted. Furthermore, we show that by taking the local covariance matrix into account, we can more accurately estimate the local geodesic distance.
When understanding LLE based on the local covariance structure, its intimate relationship with the curvature suggests a variation of LLE depending
on the “truncation scheme”. We provide a theoretical analysis of the variation.
Keywords: local covariance matrix, empirical intrinsic geometry, locally linear embedding, geodesic distance

1. Introduction
Covariance is arguably one of the most important quantities in data analysis.
It has been widely studied in the past century and is still an active research topic
nowadays. In this paper, we focus on the local covariance structure under the
manifold setup, which has been widely applied, explicitly or implicitly, to various
applications in different fields; see, for example, a far-from-complete list [11, 19, 7, 4,
32, 13, 8, 20, 10, 23, 17, 14, 1]. In the past few years, its mathematical and statistical
properties have been well established [23, 5, 3, 12, 28, 30] for different purposes. In
this paper, based on the established theoretical foundation, we extensively discuss
two topics in the manifold learning society – empirical intrinsic geometry (EIG)
and locally linear embedding (LLE).
EIG [25, 26], or originally called non-linear independent component analysis [22],
is a technique aiming to deal with the distortion underlying the collected dataset
that is caused by the observation process. In many applications, the manifold
structure we have interest cannot be directly accessed, but only via an observation.
However, the observation process might nonlinearly deform the manifold we have
interest. As a result, the information inferred from the observed data point cloud
might not faithfully reflect the intrinsic properties. The goal of EIG is correcting
this deformation by taking the local covariance matrix into account. From the
statistical viewpoint, it is a generalization of the Mahalanobis distance. While it
has been successfully applied to different problems [29, 16, 31, 21, 15], to the best
of our knowledge, a systematic evaluation of how the algorithm works under the
manifold setup, and its sensitivity to the parameter choice, is missing. Due to its
1

2

JOHN MALIK, CHAO SHEN, HAU-TIENG WU, AND NAN WU

importance, the first contribution of this paper is providing a quantification of EIG
under the manifold setup, and discuss how the chosen parameter influences the
final result. In the special case that there is no deformation (that is, we can access
the manifold directly), we show a more accurate geodesic distance estimate by
correcting the Euclidean distance when the manifold is embedded in the Euclidean
space.
LLE [19] is a widely applied nonlinear dimension reduction technique in the manifold learning society. Despite its wide application, its theoretical properties were
studied only until recently. See [30] as an example. Based on the analysis, several peculiar behaviors of LLE have been better understood. While LLE depends
on the barycentric coordinate to determine the affinity between pairs of points, it
has a natural relationship with the local covariance matrix; the kernel associated
with LLE is not symmetric, which is different from the kernel commonly used in
graph Laplacian-based algorithms like Laplacian Eigenmaps [2] or Diffusion Maps
[6]. The regularization plays an essential role in the algorithm. Different regularizations lead to different embedding results. Based on the intimate relationship
between the curvature and regularization, the second contribution of this paper is
studying a variation of LLE by directly truncating the local covariance matrix.
In Section 3, we provide a thereotical argument of EIG under the manifold
setup; when the observation process is trivial, we provide a more accurate geodesic
estimator, called covariance-corrected geodesic distance estimator. In Section 4, we
discuss the relationship of EIG and LLE, and provide a variation of LLE. Numerical
results are shown in each section to support the theoretical findings. In Section 6,
discussion and conclusion are provided.
1.1. Notations and Mathematical Setup. Let X be a p-dimensional random
vector with the range supported on a d-dimensional, compact, smooth Riemannian
manifold (M, g) isometrically embedded in Rp via ι : M ,→ Rp . We assume that M
is boundary-free in this work. Let expx : Tx M → M be the exponential map at x.
Unless otherwise specified, we will carry out calculations using normal coordinates.
Let Tx M denote the tangent space at x ∈ M , and let ι∗ Tx M denote the embedded
tangent space in Rp . Write the normal space at y = ι(x) as (ι∗ Tx M )⊥ . Let IIx be the
second fundamental form of ι at x. Let P be the probability density function (p.d.f.)
associated with the random vector X [5, Section 4]. We assume that P ∈ C 5 (ι(M ))
and that there exists 0 < Pm ≤ PM such that Pm ≤ P (y) ≤ PM < ∞ for all
y ∈ ι(M ). Let ei ∈ Rp be the unit p-dimensional unitary vector with an 1 in the
p
ith entry. Let BhR (z) denote the p-dimensional Euclidean ball of radius h > 0
with centre z ∈ Rp , and let χA : Rp → {0, 1} denote the indicator function of the
set A ⊂ Rp . We define two main quantities, the truncated inverse and regularized
inverse, of a symmetric matrix that are related to EIG and LLE respectively.
Definition 1. Let A ∈ Rp×p be a real symmetric matrix. r = rank(A). Let
λ1 ≥ . . . ≥ λp be the eigenvalues of A, and let u1 , . . . , up be the corresponding
normalized eigenvectors. For 0 < α ≤ r, the α-truncated inverse of A is defined as
 −1
  >
u1
λ1




. 
.
..
(1)
Tα [A] = u1 · · · uα 
  ..  .
λ−1
α

u>
α

COVARIANCE MATRIX, EIG AND LLE

3

Choose a regularization constant c > 0. The c-regularized inverse of A is defined as

  >
(λ1 + c)−1
u1




. 
.
..
(2)
Ic [A] = u1 · · · ur 
  ..  ,
(λr + c)−1
u>
r
Note that if α = rank(A), then Tα [A] is the Penrose-Moore pseudo-inverse.
When c → 0, Ic [A] becomes Tr [A].
2. Local Covariance Matrix and Some Facts
For x ∈ M and a measurable set O ⊂ ι(M ), define
(3)

Cx,O := E[(X − ι(x))(X − ι(x))> χO (X)] ∈ Rp×p
p

the local covariance matrix at ι(x) ∈ ι(M ) associated with O. When O is BhR (ι(x))∩
ι(M ), where h > 0, we denote Ch (x) := Cx,B Rp (ι(x))∩ι(M ) and simply call Ch (x)
h
the local covariance matrix at ι(x). The local covariance matrix and its relationship with the embedded tangent space of the manifold have been widely studied
recently, including (but not exclusively) [23, 5, 28, 3, 12, 14] and the idea has been
widely applied. Recently, in order to systematically study the LLE algorithm, the
higher order structure of the local covariance matrix was explored in [30]. We now
summarize the result for our purpose. Since the local covariance matrix is invariant
up to translation and rotation and the analysis is local at one point, to simplify
the discussion, from now on, when we analyze the local covariance matrix at x, we
always assume that the manifold is translated and rotated in Rp so that ι∗ Tx M is
spanned by e1 , . . . , ed ∈ Rp .
Lemma 1. ([30, Proposition 3.2]) Fix x ∈ M and take h > 0. Write the eigendecomposition of Ch (x) = Uh (x)Λh (x)Uh (x)> . Then when h is sufficiently small,
the diagonal matrix Λh (x) ∈ Rp×p satisfies


|S d−1 |P (x)hd+2 Id×d + O(h2 )
0
Λh (x) =
(4)
,
0
O(h2 )
d(d + 2)
and the orthogonal matrix Uh (x) ∈ O(d) satisfies


U
0
(5)
Uh (x) = 1
(Ip×p + h2 S) + O(h4 ),
0 U2
where U1 ∈ O(d), U2 ∈ O(p − d), and S is an anti-symmetric matrix.
First, note that the local covariance matrix depends on the p.d.f.. Particularly,
when the sampling is non-uniform, the eigenvalues are deviated by the p.d.f.. Equation (4) in this lemma says that the first d eigenvalues of Ch (x) are of order hd+2 ,
while the rest of the eigenvalues are two orders higher; that is, of order O(hd+4 ).
Moreover, note that Ip×p + h2 S approximates a rotation, so equation (5) says that
the first d normalized eigenvectors of Ch (x) deviate from an orthonormal basis of
ι∗ Tx M by a rotational error of order O(h2 ), and the rest of the orthonormal eigenvectors of Ch (x) deviate from an orthonormal basis of (ι∗ Tx M )⊥ within a rotational
error of order O(h2 ).
In practice, we are given a finite sampling of points from the embedded manifold
and need to approximate the local covariance matrix. Since the finite convergence
argument of the sample local covariance matrix to the local covariance matrix is
standard (see, for example, [30, Propositions 3.1 and 3.2 and Lemma E.4]), to

4

JOHN MALIK, CHAO SHEN, HAU-TIENG WU, AND NAN WU

focus on the main idea and simplify the discussion, below we work directly on the
continuous setup; that is, we consider only the asymptotical case when n → ∞.
It has been well known that when a manifold M is isometrically embedded into
another manifold M 0 , then for two close points x, y ∈ M , the geodesic distance
between x, y in M could be well approximate by the geodesic distance between x, y
in M 0 , with the error depending on the second fundamental form of the embedding
(see, for example [24, Proposition 6]). We have the following lemma when M 0 is
Euclidean space.
Lemma 2. ([30, Lemma B.2]) Suppose that M is isometrically embedded in Rp
through ι. Fix x ∈ M and use polar coordinates (t, θ) ∈ [0, ∞)×S d−1 to parametrize
Tx M . For y = expx (θt) for sufficiently small t, we have:
IIx (θ, θ) 2 ∇θ IIx (θ, θ) 3
t +
t + O(t4 ).
2
6
Moreover, when h := kι(y) − ι(x)kRp is sufficiently small, we have
(6)

ι(y) − ι(x) = (ι∗ θ)t +

(7)

h=t−

kIIx (θ, θ)k2 3 ∇θ IIx (θ, θ) · IIx (θ, θ) 4
t −
t + O(t5 ) .
24
24

t=h+

kIIx (θ, θ)k2 3 ∇θ IIx (θ, θ) · IIx (θ, θ) 4
h +
h + O(h5 ).
24
24

and hence
(8)

The proof of the Lemma could be found in, for example, in [24, Proposition 6]
when M 0 is a generic manifold or [30, Lemma B.2] when M 0 is Euclidean space.
3. Empirical Intrinsic Geometry and Local Covariance Matrix
EIG [25, 22] is a technique aiming to deal with the underlying distortion caused
by the observation process. The basic idea of the technique is that the local covariance matrix captures the distortion, under suitable assumptions, and we can
correct the distortion by manipulating the local covariance matrix and recover the
local geodesic distance. From the statistical viewpoint, it is a generalization of the
Mahalanobis distance. We now examine the intimate relationship between EIG and
the covariance-corrected geodesic estimator.
Suppose that M is a d dimensional closed Riemannian manifold that hosts the
information we have interest but we cannot directly access. We model the observation process as Φ : M → N , where Φ is a diffeomorphism and N is isometrically
embedded in Rq via ι. Under this setup, for the point x ∈ M that we cannot
access, there is a corresponding point ι(y) ∈ ι(N ) ⊂ Rq that we can access and
collect as the dataset, where y = Φ(x). The mission is estimating the geodesic
distance between two close points x ∈ M and w ∈ M through accessible data
points ι(y) = ι(Φ(x)) and ι(z) = ι(Φ(w)). This situation is commonly encountered
in data analysis; for example, the brain activity information contained in recorded
electroencephalogram signals might be deformed due to the recording process and
the electrophysiological properties. Clearly, due to the diffeomorphism, the pairwise
distance is no longer correct, and unsupervised learning algorithms, like diffusion
map, might lead to biased results. The interest of EIG is correcting this bias.
The main idea of EIG can be summarized in the following way. For each hidden
point x ∈ M , take the geodesic ball B (x) ⊂ M with the radius  and centered at
x. As is discussed in [22], if we can determine E(y) := Φ(B (x)) ⊂ N around y

COVARIANCE MATRIX, EIG AND LLE

5

E(y)
z" •+
y+

(Nd,+gN)

+

q+

(" X) #+=+Q(x)dVN(x)+
":+diﬀeomorphsim+
(Md,+gM)
ε+

X+
($,+%,+#)+

•+
x+

X #+=+P(x)dVM(x)+

w"

Bε(x)

Figure 1. The illustration of the theoretical framework of EIG.
that is associated with B (x), then up to a constant, the geodesic distance between
two close points w and x can be well approximated by ι(y), ι(z), and the local
covariance matrix Cy,ι(E(y)) associated with ι(E(y)) as defined in equation (3). To
simplify the notation, denote C (y) := Cy,ι(E(y)) . See Figure 1 for an illustration
of the setup. Note that when M is Euclidean and Φ is linear, E(y) is an ellipsoid.
In the following, although we consider the manifold model that is in general not
Euclidean, we abuse the terminology and call E(y) an ellipsoid. Numerically, C (y)
is estimated by taking points in Y ∩E(y) into account. Again, since the convergence
proof is standard, to simplify the discussion, we skip the finite convergence step and
focus on the continuous setup.
This main idea is carried out by noting that the local covariance matrix at
y associated with C (y) captures the Jacobian of the diffeomorphism associated
with the observation, that is, C (y) ≈ ∇Φ|x Φ|>
x [22]. Under the assumption that
the ellipsoid E(y) is known, C (y) ≈ ∇Φ|x Φ|>
x is true and the dimension of the
manifold d is known, authors in [22, 25] consider the following quantity, called the
EIG distance, to estimate the geodesic distance between x and w:


Tα [C (y)] + Tα [C (z)]
>
(9)
(ι(z) − ι(y)) ,
(ι(z) − ι(y))
2
where α is chosen to be the dimension of d.
3.1. Challenges of EIG and analysis. However, there are two challenges of
applying the EIG idea: how to determine the ellipsoid E(y) associated with B (x),
and how to estimate the dimension d of the intrinsic manifold. In [22] and most of its
citations, the state space model combined with the stochastic differential equation
(SDE) is considered, and the ellipsoid can be determined simply by taking the
temporal relationship into account under this model. We refer readers with interest
to [22, 25] for more details about this state-space model. If the state-space model

6

JOHN MALIK, CHAO SHEN, HAU-TIENG WU, AND NAN WU

and SDE cannot be directly applied, this task is most of time a big challenge and to
the best of our knowledge not too much is known. This challenge is however out of
the scope of this paper. On the other hand, most of time we do not have an access
to the dimension of d, and a dimension estimation is needed. Although theoretically
we can count on the spectral gap to estimate d, in practice it depends on the tryand-error process and little can be guaranteed. Particularly, when the manifold is
highly distorted by the observation, faithfully estimating the dimension is another
challenging task. To the best of our knowledge, although EIG has obtained several
successes in different applications [29, 16, 31, 21], a systematic exploration of the
approximation under the general manifold setup is lacking. It is also not clear what
may happen when the dimension is not estimated correctly.
Below, we assume that we have had the knowledge of the ellipsoid associated
with B (x) for all x, and hence we can evaluate the local covariance matrix at
y = Φ(x) associated with E(y). We show how C (y) ≈ ∇Φ|x Φ|>
x holds under the
manifold setup, and the influenced of an erroneously estimated dimension. We start
with the following definitions.
Definition 2. For y ∈ N , define the normalized local covariance matrix at y
associated with E(y) ⊂ N as
C¯ (y) :=

(10)

C (y)
2 E[χE(y) (Y )]

The normalization step in (10) is introduced to remove the impact of the nonuniform sampling. As we will see, without this normalization, when the sampling is
non-uniform, the p.d.f. will play a role in the final analysis.
Lemma 3. Let x ∈ M and B (x) ⊂ M be the geodesic ball around x. The local
covariance matrix at y = Φ(x) associated with the ellipsoid E(y) = Φ(B (x)) ⊂ N
satisfies
(11)

C (y) =

|S d−1 ||∇Φ(x)|P (x)d+2
d+4
[ι∗ |y ∇Φ(x)][∇Φ(x)> ι∗ |>
)
y ] + O(
d(d + 2)

and the normalized local covariance matrix at y associated with E(y) satisfies
(12)

C¯ (y) =

1
2
[ι∗ |y ∇Φ(x)][∇Φ(x)> ι∗ |>
y ] + O( ) .
d+2

If v1 , v2 ∈ (ι∗ Ty N )⊥ , then
v1> C¯ (y)v2 =

d2
d−1
4|S
|(d + 4)

Z
S d−1

v1> IIy (∇θ Φ(x), ∇θ Φ(x))(IIy (∇θ Φ(x), ∇θ Φ(x)))> v2 dθ + O(4 ) .

By this lemma, we see that the expansion of C (y) depends not only on the p.d.f.
but also the Jacobian of the deformation, and the normalization step cancels this
dependence. The proof is postponed to Appendix, page 17. Based on this lemma
and the considered setup in [22, 25], we consider the following quantity.
Definition 3. For y, z ∈ N , define the EIG distance of order α between y, z, where
1 ≤ α ≤ q, as


Tα [C¯ (y)] + Tα [C¯ (z)]
2
>
(13)
EIGα (y, z) = (ι(y) − ι(z))
(ι(y) − ι(z)) .
2

COVARIANCE MATRIX, EIG AND LLE

7

Note that this definition is slightly different from that considered in (9). Below
we show that the EIG distance between y = Φ(x) ∈ N and z = Φ(w) ∈ N is a good
estimator of the geodesic distance between x ∈ M and w ∈ M only when some
conditions are satisfied.
Theorem 1. Suppose two d-dimensional smooth closed Riemannian manifolds M
and N are diffeomorphic via Φ : M → N , and suppose N is isometrically embedded
in Rq via ι. Take x ∈ M and w ∈ B (x) ⊂ M , where  > 0, and denote t to
be the geodesic distance between x and w. Let y = Φ(x), z = Φ(w). Take α ≤
min(rank[C¯ (y)], rank[C¯ (z)]). When  is sufficiently small, we have the following:
(1) When α = d, we have
√
EIGd (y, z) = d + 2t + O(t2 ) + O(t3 ) ;
(2) When 1 ≤ α < d, we have
(14)

EIGα (y, z) =

√

d + 2 t + O(t) .

Let Vα (y) and Vα (z) be the subspaces of ι∗ Ty N and ι∗ Tz N generated by the
first α eigenvectors of C¯ (y) and C¯ (z) respectively. Suppose z = expy ϑ(y)
for ϑ(y) ∈ Ty N and y = expz ϑ(z) for ϑ(z) ∈ Tz N . If ι∗ |y ϑ(y) ∈ Vα (y) and
ι∗ |z ϑ(z) ∈ Vα (z), we have
√
(15)
EIGα (y, z) = d + 2 t + O(t2 ) + O(t3 ) .
(3) When α > d, assume that the smallest nonzero eigenvalues of C¯ (y) and
C¯ (z) are of order 4 . Let t = β for β > 1, when  is sufficiently small, we
have
√
(16)
EIGα (y, z) = d + 2 t + O(t(t/ + )2 + t(t2 /2 + 1)2 ) ,
and when t = , we have
(17)

EIGα (y, z) =

√

d + 2 t + O(t) .

The proof is postponed to Section A. We now have discussion of the theorem.
In general, rank[C¯ (y)] and rank[C¯ (z)] are different, thus we need the condition
α ≤ min(rank[C¯ (y)], rank[C¯ (z)]). First of all, if we know the dimension of the
manifold, then the EIG distance of order d between y = Φ(x) ∈ N and z = Φ(w) ∈
N is an accurate estimator
of the geodesic distance between x ∈ M and w ∈ M , up
√
to a global constant d + 2. This result coincides with the claim in [22, 25] when
the sampling on N is uniform. However, when the sampling on N is non-uniform,
the result is deviated by the p.d.f. if we replace the normalized local covariance
matrix C¯ (y) in (13) by the local covariance matrix C (y). This deviation can be
seen from comparing (11) and (12) in Lemma 3. We mention that in practice,
if the purpose of applying EIG is finding the k nearest neighbors of x ∈ M but
not determining the precise geodesic distance, then this normalization step is not
critical.
Second, if the dimension is unknown and wrongly estimated, the result depends
on the situation. When the dimension is under-estimated, in general the estimator
is wrong. In a non-generic situation that the geodesic direction from y to z and
that from z to y are both located on the first α eigenvectors of the associated local
covariant matrices, we may still obtain an accurate estimator. Note that due to the
curvature, there is no guarantee that the first α eigenvectors of C¯ (y) will contain
ϑ(y); even if it does, there is no guarantee that the first α eigenvectors of C¯ (z)

8

JOHN MALIK, CHAO SHEN, HAU-TIENG WU, AND NAN WU

will contain ϑ(z). When the dimension is over-estimated, then we still obtain a
reasonably well estimate of the intrinsic geodesic distance when t/ is sufficiently
small, but with a slow convergence rate.
The assumption “the smallest nonzero eigenvalues of C¯ (y) and C¯ (z) are of order
4
 ” needs some discussion. The rule of thumb is that eigenvalues of C¯ (y) associated
with the tangent directions are of order 1, but the eigenvalues associated with the
normal directions are of order 2 or higher. The assumption that the smallest
eigenvalues are of order 4 means that when the principle curvature is zero, the
higher order curvature is not too small. As is shown in the proof of this theorem,
under this assumption, we can control the error induced by the interaction between
ι(y)−ι(z) and Tα [C¯ (y)]. However, when there is an eigenvalue of even higher order,
the interaction between ι(y) − ι(z) and Tα [C¯ (y)] becomes more complicated, and a
more delicate analysis involving a systematic recursive formula of Taylor expansion
of ι(y) − ι(z) and Tα [C¯ (y)] is needed. We will explore this issue in our future work.
To sum up, in practice if we are not confident about the estimated dimension, we
may want to choose a larger α if the dimension estimation is not accurate. Although
we do not discuss it in this paper, we mention that this conservative way might not
be that preferred if noise exists, since the noise will contaminate the eigenvector
associated with the small eigenvalues. A more statistical approach to solve these
issues will be studied in our future work.

3.2. A special EIG setup – covariance-corrected geodesic estimator. When
Φ is identity, EIG is reduced to the ordinary manifold learning setup; that is, we
have samples from the manifold and we want to estimate its geometric structure.
By Lemma 2, when two points on the smooth manifold are close enough, we can
accurately estimate their geodesic distance t by the ambient Euclidean distance
h = kι(y) − ι(x)kRp up to the third order O(t3 ), and the third order term is essentially the second fundamental form. We call h the Euclidean-distance-based geodesic
distance estimator. In general, it is not an easy task to directly estimate the second
fundamental form from the point cloud if M 0 is not Euclidean space. However,
when M 0 is Euclidean, like the setup in Lemma 2, the second fundamental form
information could be well approximated by the local covariance matrix. We define
the following projection operator associated with the local covariance matrix.
Definition 4. Suppose M is a closed, smooth d dimensional Riemannian manifold
isometrically embedded in Rp through ι. Fix x ∈ M . For h̄ > 0, let Ch̄ (x) :=
Cx,B Rp (ι(x))∩ι(M ) as defined in equation (3). Let λ1 ≥ . . . ≥ λp be the eigenvalues
h̄
of Ch̄ (x), and let u1 , . . . , up be the corresponding normalized eigenvectors. Then,
for any z ∈ Rp , we define the orthogonal projection associated with Ch̄ (x), Ph̄⊥ (z),
to be the projection of z onto the subspace spanned by ud+1 , . . . , up .
The name orthogonal projection follows from Lemma 1, where we prove that
ud+1 , . . . , up deviate from an orthonormal basis of (ι∗ Tx M )⊥ within a rotational
error of order O(h̄2 ). Inspired by Lemmas 1 and 2, we have the following theorem.
The proof is postponed to Section B.
p

Theorem 2. Fix x ∈ M . Suppose that y ∈ M and ι(y) ∈ Bh̄R (ι(x)). We use polar
coordinates (t, θ) ∈ [0, ∞) × S d−1 to parametrize Tx M and we have y = expx (θt)

COVARIANCE MATRIX, EIG AND LLE

9

for t > 0. Denote h := kι(y) − ι(x)kRp ≤ h̄. When h̄ is small enough,
!
kPh̄⊥ (ι(y) − ι(x))k2Rp
t− h+
= O(h2 h̄2 ) .
(18)
6h
Based on this Theorem, we have the following definition.
Definition 5. Following the notation in Theorem 2, the covariance-corrected geodesic distance estimator between points x, y ∈ M is defined as
(19)

h+

kPh̄⊥ (ι(y) − ι(x))k2Rp
.
6h

Clearly, the covariance-corrected geodesic distance estimator accurately estimates the geodesic distance t up to the fourth order. It is worth mentioning that
a better estimation for the geodesic distance cannot be directly achieved only with
the local covariance matrix. Indeed, to estimate d(x, y) within an error of O(h5 ),
by Lemma 2, we need to estimate ∇θ IIx (θ, θ)⊥ within an error of order O(h). Thus,
more information is needed if a more accurate geodesic distance estimate is needed.
4. Locally Linear Embedding and a Variation by Truncation
LLE [19] is a widely applied dimension reduction technique in the manifold
learning society. Some of its theoretical properties have been explored in [30]. For
example, under the manifold setup, the barycentric coordinate is intimately related
to the local covariance structure of the manifold. By reformulate the barycentric coordinate in the local covariance structure, the natural kernel associated with
LLE is discovered, and it is very different from the kernel in graph Laplacian-based
algorithms like Laplacian Eigenmaps [2], Diffusion Maps [6] or commute time embedding [18]. The regularization step in LLE is crucial in the whole algorithm,
and it might provide different embeddings if the regularization is not chosen correctly. In this section, we connect LLE and EIG, and take the geometric structure
to explore a variation of LLE via truncation.
4.1. A summary of LLE algorithm. We start with recalling the LLE algorithm in the finite sampling setup. Let X = {xi }ni=1 ⊂ Rp denote a set of point

p
cloud. Fix one point xk ∈ X . For h > 0, assume Nxk := X ∩ BhR (xk ) \ {xk } =
{xk,1 , . . . , xk,Nk }, where h > 0. Here we take the h-radius ball to determine neighbors to be consistent with the following discussion. In practice we can use the
k-nearest neighbors, where k ∈ N is determined by the user. The relationship
between the h-radius ball and k-nearest neighbors is discussed in [30, Section 5].
Define the local data matrix associated with xk as


|
|
Gn,h (xk ) := xk,1 − xk · · · xk,Nk − xk  ∈ Rp×Nk .
|
|
With the local data matrix, the LLE algorithm is composed of three steps. First,
for each xk ∈ X , find its barycentric coordinate associated with Nxk by solving
(20)

wxk =

argmin
w∈RNk ,w> 1Nk =1

kxk −

Nk
X
j=1

w(j)xk,j k2 ∈ RNk .

10

JOHN MALIK, CHAO SHEN, HAU-TIENG WU, AND NAN WU

Solving (20) is equivalent to minimizing w> Gn,h (xk )> Gn,h (xk )w over w ∈ RNk
under the constraint w> 1Nk = 1. Since in general Gn,h (xk )> Gn,h (xk ) is singular,
it is recommended in [19] to solve (Gn,h (xk )> Gn,h (xk ) + cINk ×Nk )y = 1Nk , where
c > 0 is the regularizer chosen by the user, and hence obtain
(21)

w xk =

(Gn,h (xk )> Gn,h (xk ) + cINk ×Nk )−1 1Nk
.
((Gn,h (xk )> Gn,h (xk ) + cINk ×Nk )−1 1Nk )> 1N

By viewing wxk as the “affinity” of xk and its neighbors, the second step is defining
a LLE matrix W ∈ Rn×n by

wxk (j) if xl = xk,j ∈ Nxk ;
(22)
Wk,l =
0
otherwise.
Finally, to reduce the dimension of X or to visualize X , it is suggested in [19] to
embed X into a low dimension Euclidean space
xk 7→ [v1 (k), · · · , v` (k)]> ∈ R` ,

(23)

for each xk ∈ X , where ` is the dimension of the embedded points chosen by the
user, and v1 , · · · , v` ∈ Rn are eigenvectors of (I − W )> (I − W ) corresponding to
the ` smallest eigenvalues. Note that in general W is not symmetric, and this is
why the spectrum of (I − W )> (I − W ) is suggested but not that of I − W .
It is shown in [30, Section 2] that by taking the relationship between the sample covariance matrix Cn,h (xk ) := Gn,h (xk )Gn,h (xk )> and the Gramian matrix
Gn,h (xk )> Gn,h (xk ) into account, we can rewrite (21) as
(24)

wx>k =

>
>
1>
Nk − 1Nk Gn,h (xk ) Ic (Cn,h (xk ))Gn,h (xk )
>
Nk − 1>
Nk Gn,h (xk ) Ic (Cn,h (xk ))Gn,h (xk )1Nk

,

>
>
and view 1>
Nk −1Nk Gn,h (xk ) Ic (Cn,h (xk ))Gn,h (xk ) as the “kernel” associated with
the LLE algorithm. It is clear that in general this kernel has negative values, and
the LLE matrix is not a transition matrix. By defining

Tn,xk := Ic (Cn,h (xk ))Gn,h (xk )1Nk

(25)

we obtain the result claimed in [30, Proposition 2.1].
4.2. LLE and Mahalanobis distance in the continuous setup. The above is
for the general dataset. When the dataset X = {ι(xi )}ni=1 ⊂ ι(M ) ⊂ Rp is sampled
from a manifold, with (24) and (25), the asymptotical behavior of LLE under the
manifold setup is explored in [30]. As is shown in [30, (3.12)], in the continuous
setup, the unnormalized kernel associated with LLE at x ∈ M is
(26)

p
KhLLE (x, y) = [1 − T>
ι(x) (ι(y) − ι(x))]χB R (ι(x))∩ι(M ) (ι(y)),
h

where y ∈ M and
(27)



Tι(x) := Ic (Ch (x)) E(X − ι(x))χB Rp (x) (X) ∈ Rp ,
h

and hence the normalized kernel is defined by
(28)

PhLLE (x, y) := R

KhLLE (x, y)
.
KhLLE (x, y)P (y)dV (y)
M

Clearly, (24) and (25) are discretization of PhLLE (x, y) and Tι(x) . The convergence
behavior could be found in [30, (3.7) and Theorem 3.1].

COVARIANCE MATRIX, EIG AND LLE

11

There are two geometric facts we should mention. First, rewriting (21) as (24)
is equivalent to representing (21) in the coordinate of the frame bundle, and hence
an explicit form of the kernel underlying LLE. Second, the unnormalized kernel
at x involves projecting ι(y) − ι(x) onto the vector space spanned by u1 , . . . , ur ,
and scaling the resulting vector component-wise by (λ1 + c)−1/2 , . . . , (λr + c)−1/2 .
When r > d, the involvement of the regularizer c and the eigenvalues λd+1 , . . . , λr
are associated with the normal bundle. An important intuition is that when c →
0, if we put aside the expectation in T>
ι(x) (ι(y) − ι(x)) and replace E(X − ι(x))
by (ι(y) − ι(x)), we obtain the term (ι(y) − ι(x))> Tr (Ch (x))(ι(y) − ι(x)), which
can be understood as a variation of the Mahalanobis distance. Therefore, we can
interpret LLE as taking the neighbor information (the 0-1 kernel in (26), that is,
1 × χB Rp (ι(x))∩ι(M ) (ι(y))), and the Mahalanobis distance of neighboring points into
h
account to design the kernel.
Lemma 1 guarantees that the first d eigenvalues of Ch (x) are sufficiently large.
However, the remaining (r − d) non-trivial eigenvalues are small. Moreover, the
remaining (r − d) non-trivial eigenvalues are related to the curvature of the manifold. The importance of choosing the regularization constant c has been extensively
discussed in [30, Theorem 3.2], and it is known that c is to adjust the curvature
information in these scaling factors (λd+1 + c)−1/2 , . . . , (λr + c)−1/2 for different
purposes. To be more precise, the regularization c plays a role of “radio tuner”.
The larger the c is, the more enhanced the p.d.f. information will be; the smaller the
c is, the more dominated the curvature information will be [30, Theorem 3.2]. It is
shown that if we want to obtain the Laplace-Beltrami operator of the Riemannian
manifold, c should be chosen properly so as to suppress the curvature information
contained in the (d + 1)-th to the r-th eigenvalues.
4.3. A variation of LLE and EIG. Under the manifold setup, the above discussion suggests that the curvature information associated with the d + 1, . . . , r eigenvalues and eigenvectors is not needed if we want to obtain the Laplace-Beltrami
operator. Therefore, an alternative to choosing a regularization parameter is a direct truncation by taking only the largest d eigenvalues into account; that is, replace
Ic [Ch (x)] in (26) by Td [Ch (x)], and define the normalized truncated LLE kernel at
x ∈ M by
(29)

PhtLLE (x, y) := R

KhtLLE (x, y)
tLLE
Kh (x, y)P (y)dV
M

(y)

,

where
(30)

p
KhtLLE (x, y) := [1 − T̃>
ι(x) (ι(y) − ι(x))]χB R (ι(x))∩ι(M ) (ι(y)),
h

y ∈ M and
(31)



T̃ι(x) := Td (Ch (x)) E(X − ι(x))χB Rp (x) (X) ∈ Rp .
h

Numerically when we have only finite data points, we run the same LLE algorithm,
but replace wx>k and Tn,xk by the following terms:
(32)

w̃x>k :=

>
1>
Nk − T̃n,xk Gn,h (xk )

Nk − T̃>
n,xk Gn,h (xk )1Nk

,

12

JOHN MALIK, CHAO SHEN, HAU-TIENG WU, AND NAN WU

where
T̃n,xk := Td (Gn,h (xk )Gn,h (xk )> )Gn,h (xk )1Nk ;

(33)

that is, instead of running a regularized pseudo-inversion, we run a truncated inversion. LLE with low-dimensional neighborhood representation (LDR-LLE) proposed
in [9] is an algorithm related to this truncated idea. While the geometric structure
of local covariance matrix is not specifically discussed in LDR-LLE, it can be systematically studied in our framework. To simplify the nomination, we also call the
above truncation scheme LDR-LLE.
Geometrically, (ι(y) − ι(x))> Td [Ch (x)](ι(y) − ι(x)) evaluates geodesic distances
between points in the h-neighbourhood of x using a method directly related to
both EIG in (9) when Φ is identity. With the normalized truncated LLE kernel,
we can proceed with the standard LLE dimension reduction step. Therefore, the
LDR-LLE as taking the neighbor information, 1 × χB Rp (ι(x))∩ι(M ) (ι(y))), and the
h
EIG distance of neighboring points into account to design the kernel.
Under the same condition specified in [30, Theorem 3.2], LDR-LLE has the
same asymptotical behavior as that in [30, Theorem 3.3] with the properly chosen
regularization; that is,
4
Theorem
√ 3. Assume the setup in Section 1.1. Take f ∈ C (M ). If h = h(n) so
log(n)
that n1/2 hd/2+1 → 0 and h → 0 when n → ∞, with probability greater than 1 − n−2 ,
for all xk ∈ X we have
!
p
Z
Nk
X
log(n)
w̃n,xk (j)f (xk,j ) =
PhtLLE (xk , y)f (y)P (y)dV (y) + O
n1/2 d/2−1
M
j=1

and
Z
M

PhtLLE (x, y)f (y)P (y)dV (y) = f (y) +

h2
∆f (x) + O(h3 ) .
2(d + 2)

As a result, with probability greater than 1 − n−2 , for all xk ∈ X we have
!
p
Nk
i
log(n)
1 hX
1
∆f (x) + O(h) + O
.
w̃n,xk (i)f (xk,j ) − f (xk ) =
h2 j=1
2(d + 2)
n1/2 d/2+1
In other words, by the LDR-LLE, we recover the Laplace-Beltrami operator of the
manifold even if the second fundamental form is non-trivial and the sampling is
non-uniform. While the proof is a direct modification of the proof of [30, Theorem
3.2], we provide the proof in Appendix C for the sake of self-containedness. We
mention that the finite sample convergence has the same statement and follows
exactly the same line as that of [30, Theorem 3.1], so we skip it.
Based on the theoretical development, we identify two benefits of LDR-LLE.
First, the normal bundle information is automatically removed by the truncation,
and we obtain the intrinsic geometric quantity ∆. Second, the non-uniform sampling effect is automatically removed, and no kernel density estimation is needed.
This comes from the fact that the EIG part of the kernel T̃ι(x) (or the Mahalanobis
distance part of the kernel Tι(x) ) automatically provides the density information
(see Lemma A.7). Compared with DM, since we do not need normalization as is
proposed in [6] to eliminate the non-uniform sampling effect, the convergence rate
is faster, as is stated in [30, Theorem 3.3].

COVARIANCE MATRIX, EIG AND LLE
#10!5

#10!3

Euclidean Distance
Covariance Correction

Euclidean Distance
Covariance Correction

1

Absolute Error

Absolute Error

1.5

13

1

0.5

0.8
0.6
0.4
0.2

0

0

0.05

0.1

0.15

0

0

5

Geodesic Distance to xk

10

15

20

25

Geodesic Distance to xk

Figure 2. Comparison of the Euclidean-distance-based geodesic
distance estimator (black) and the covariance-corrected geodesic
distance estimator (red) in the logarithmic spiral L. Left: local
distance for  = 0.2. Right: Global geodesic distance estimation by
applying Dijkstra’s algorithm with different local geodesic distance
estimators.
However, we emphasize that although this approach leads to aesthetic asymptotical properties under the manifold assumption, it may not work well when we
do not have a manifold structure. For a general graph or point cloud, we found it
more reliable to apply the original LLE, and a model or analysis explaining why
is needed. Last but not least, note that this action requires that we know the
dimension of the manifold, so even if we know the point cloud is sampled from a
manifold, we need to estimate its dimension before applying the LDR-LLE. This
requirement might render LDR-LLE less applicable in non-manifold data.
5. Numerical Results
We demonstrate some numerical results associated with the studied algorithms
in this paper. Uniformly sample n = 8000 points from the logarithmic spiral
L ⊂ R2 , a 1-dim manifold parametrized as x(s) = ( √s2 + 1) cos log( √s2 + 1) and
y(s) = ( √s2 + 1) sin log( √s2 + 1), and choose  = 0.2. We fix xk ∈ L and plot the
absolute error between the true geodesic distance and the Euclidean and covariancecorrected distances. The result is shown in Figure 2. Note that the two “branches”
correspond to differences in curvature to either side of xk . In some manifold learning
algorithms, knowing the local geodesic distance is not enough, and we need to
know the global geodesic distance between any two points. ISOMAP [27] is a
typical example. Usually, we estimate the global geodesic distance by applying
Dijkstra’s algorithm to find the shortest-path distance between all pairs, where
the geodesic distance between neighboring points is estimated by the Euclidean
distance. By combining the covariance-corrected geodesic distance estimator and
Dijkstra’s algorithm, we achieve a more accurate global geodesic distance, as is
shown in Figure 2.
We show in Figure 3 that on S 1 ⊂ R2 with a non-uniform sampling, the eigenvalues of the Laplace-Beltrami operator can be accurately recovered using the LDRLLE. When compared with the eigenvalues obtained from the 1-normalized DM,

14

JOHN MALIK, CHAO SHEN, HAU-TIENG WU, AND NAN WU

Truncated LLE
Di,usion Map

Absolute Error

150

100

50

0

0

500

1000

1500

Eigenvalue of "

Figure 3. The eigenvalues of the Laplace-Beltrami operator on
S 1 are well-recovered by the LDR-LLE in the presence of a nonuniform sampling, compared with the α-normalized DM.
we see that the LDR-LLE performs better, which comes from the faster convergence rate of LDR-LLE. We use n = 8000 points and we choose h = 0.03. The
non-uniform sampling is obtained as follows. Sample n points {θi }ni=1 uniformly
from [0, 1] and set xi = (cos(2π(θi + 0.3 sin(θi ))), sin(2π(θi + 0.3 sin(θi )))) ∈ S 1 ,
where i = 1, . . . , n.
6. Conclusion
In this paper, we extend the knowledge of local covariance matrix, particularly
the higher order expansion of the local covariance matrix, to study two commonly
applied manifold learning algorithms, EIG and LLE. We provide a theoretical analysis of EIG under different situations under the manifold setup, and emphasize
the importance of correctly estimating the dimension and its sensitivity to chosen
parameters. The curvature effect is specifically carefully discussed. Under the trivial EIG setup when there is no deformation, that is, Φ is identity, we show that
the local covariance matrix structure allows us to obtain a more accurate geodesic
distance for neighboring points. The geometric relationship between local covariance matrix and LLE leads to a natural generalization of LLE by taking EIG into
account. We provide a theoretical justification of LDR-LLE and compare its pros
and cons with the original LLE.
7. Acknowledgement
Hau-Tieng Wu thanks Professor Ronald Coifman, Professor Ronen Talmon, and
Professor Stefan Steinerberger for various discussion of the topic.
References
[1] E. Arias-Castro, G. Lerman, and T. Zhang. Spectral Clustering Based on Local PCA. JMLR,
18:1–57, 2017.
[2] M. Belkin and P. Niyogi. Laplacian eigenmaps for dimensionality reduction and data representation. Neural computation, 15(6):1373–1396, 2003.
[3] A. Bernstein and A. Kuleshov. Data-based Manifold Reconstruction via Tangent Bundle
Manifold Learning. ICML-2014, Topological Methods for Machine Learning Workshop, pages
1–6, 2014.

COVARIANCE MATRIX, EIG AND LLE

15

[4] M. Brand. Charting a manifold. In NIPS, volume 15, pages 961–968, 2003.
[5] M.-Y. Cheng and H.-T. Wu. Local linear regression on manifolds and its geometric interpretation. J. Am. Stat. Assoc., 108:1421–1434, 2013.
[6] R. R. Coifman and S. Lafon. Diffusion maps. Appl. Comput. Harmon. Anal., 21(1):5–30,
2006.
[7] D. L. Donoho and C. Grimes. Hessian eigenmaps: Locally linear embedding techniques for
high-dimensional data. P. Natl. Acad. Sci. USA, 100(10):5591–5596, 2003.
[8] A. B. Goldberg, X. Zhu, A. Singh, Z. Xu, and R. Nowak. Multi-Manifold Semi-Supervised
Learning. In NIPS, pages 169–176, 2009.
[9] Y. Goldberg and Y. Ritov. LDR-LLE: LLE with low-dimensional neighborhood representation. In ISVC, pages 43–54, 2008.
[10] D. Gong, X. Zhao, and G. Medioni. Robust Multiple Manifolds Structure Learning. In ICML,
number 1, 2012.
[11] N. Kambhatla and T. K. Leen. Dimension reduction by local principal component analysis.
Neural Computation, 9(7):1493–1516, 1997.
[12] D. N. Kaslovsky and F. G. Meyer. Non-asymptotic analysis of tangent space perturbation.
Information and Inference, 3(2):134–187, 2014.
[13] D. Kushnir, M. Galun, and A. Brandt. Fast multiscale clustering and manifold identification.
Pattern Recognition, 39(10):1876–1891, 2006.
[14] A. V. Little, M. Maggioni, and L. Rosasco. Multiscale geometric methods for data sets I:
Multiscale SVD, noise and curvature. Appl. Comput. Harmon. Anal., 43(3):504–567, 2017.
[15] G.-R. Liu, Y.-L. Lo, Y.-C. Sheu, and H.-T. Wu. Diffuse to fuse EEG spectra – intrinsic
geometry of sleep dynamics for classification. ArXiv 1803.01710, 2018.
[16] G. Mishne, R. Talmon, and I. Cohen. Graph-based supervised automatic target detection.
IEEE Transactions on Geoscience and Remote Sensing, 53(5):2738–2754, 2015.
[17] S. Pedagadi, J. Orwell, S. Velastin, and B. Boghossian. Local fisher discriminant analysis
for pedestrian re-identification. In Proceedings of the IEEE Computer Society Conference on
Computer Vision and Pattern Recognition, pages 3318–3325, 2013.
[18] H. J. Qiu and E. R. Hancock. Clustering and embedding using commute times. IEEE Trans.
Patt. Anal. Mach. Intell., 29(11):1873–90, 2007.
[19] S. T. Roweis and L. K. Saul. Nonlinear dimensionality reduction by locally linear embedding.
Science, 290(5500):2323–2326, 2000.
[20] M. Salhov, G. Wolf, and A. Averbuch. Patch-to-tensor embedding. Appl. Comput. Harmon.
Anal., 33(2):182–203, 2012.
[21] A. Shemesh, R. Talmon, O. Karp, I. Amir, M. Bar, and Y. J. Grobman. Affective response
to architecture investigating human reaction to spaces with different geometry. Architectural
Science Review, 60(2):116–125, 2017.
[22] A. Singer and R. Coifman. Non-linear independent component analysis with diffusion maps.
Appl. Comput. Harmon. Anal., 25(2):226–239, 2008.
[23] A. Singer and H.-T. Wu. Vector diffusion maps and the connection Laplacian. Comm. Pure
Appl. Math., 65(8):1067–1144, 2012.
[24] O. Smolyanov, H.v. Weizsacker, and O. Wittich. Chernoff’s theorem and discrete time approximations of brownian motion on manifolds. Potential Anal., 26(1):1–29, 2007.
[25] R. Talmon and R. R. Coifman. Empirical intrinsic modeling of signals and information geometry. Yale-CS TR, 1467, 2012.
[26] R. Talmon and R. R. Coifman. Empirical intrinsic geometry for nonlinear modeling and time
series filtering. P. Natl. Acad. Sci. USA, 110(31):12535–12540, 2013.
[27] J. B. Tenenbaum, V. De Silva, and J. C. Langford. A global geometric framework for nonlinear
dimensionality reduction. Science, 290(5500):2319–2323, 2000.
[28] H. Tyagi, E. Vural, and P. Frossard. Tangent space estimation for smooth embeddings of
Riemannian manifolds. Information and Inference, 2(1):69–114, 2013.
[29] H.-T. Wu, R. Talmon, and Y.-L. Lo. Assess sleep stage by modern signal processing techniques. IEEE Trans. Biomed. Eng., 62(4):1159–1168, 2015.
[30] H.-T. Wu and N. Wu. Think globally, fit locally under the manifold setup: Asymptotic
analysis of locally linear embedding. Annals of Statistics, In press, 2018.
[31] O. Yair and R. Talmon. Local Canonical Correlation Analysis for Nonlinear Common Variables Discovery. IEEE Transactions on Signal Processing, 65(5):1101–1115, 2017.

16

JOHN MALIK, CHAO SHEN, HAU-TIENG WU, AND NAN WU

[32] Z. Zhang and H. Zha. Principal manifolds and nonlinear dimensionality reduction via tangent
space alignment. SIAM J. Sci. Comput., 26:313 – 338, 2004.

Appendix A. Proof of Theorem 1
Suppose that M and N are d dimensional closed Riemannian manifolds. Φ :
M → N is a diffeomorphsim. Moreover, suppose that N is isometrically embedded
in Rq via ι. In this section, all the derivatives of Φ are calculated in normal coordinates. For any y ∈ Rq , we identify Ty Rq with Rq . The proof of the theorem is
composed of three lemmas, including Lemma 3.
Lemma A.4. Fix x ∈ M and use polar coordinates (t, θ) ∈ [0, ∞) × S d−1 to
parametrize Tx M . For w = expx (θt) for sufficiently small t, let y = Φ(x) and
z = Φ(w), then we have:
(34)

ι(z)−ι(y) = ι∗ |y ϑt +
+


1
IIy (ϑ, ϑ) + ι∗ |y ∇2θθ Φ(x) t2
2


1
∇ϑ IIy (ϑ, ϑ) + 3IIy (∇2θθ Φ(x), ϑ) + ι∗ |y ∇3θθθ Φ(x) t3 + O(t4 ),
6

where ϑ := ∇θ Φ(x) ∈ Ty N , and
(35)

2t2 = (ι(z) − ι(y))> ι∗ |y [∇Φ(x)∇Φ(x)> ]−1 ι∗ |>
y (ι(z) − ι(y))
>

4
+ ι(y) − ι(z) ι∗ |z [∇Φ(w)∇Φ(w)> ]−1 ι∗ |>
z ι(y) − ι(z) + O(t ).

Above lemma can be regarded as a generalization of Lemma 2. Indeed, when
N = M and Φ is the identity, ∇Φ is identity and we recover Lemma 2.
Proof. Equation (34) follows from a direct Taylor expansion. Since t is small
enough, let γ : [0, t] 7→ M to be the unique unit speed geodesic from x to w;
that is, γ(0) = x and γ(t) = w. By equation (34), we have
ι(z) − ι(y) =ι∗ |y ϑt +


1
IIy (ϑ, ϑ) + ι∗ |y ∇2γ 0 (0)γ 0 (0) Φ(x) t2 + O(t3 ) ,
2

where ϑ := ∇θ Φ(x) = ∇Φ|x (θ). Since the calculation is made in normal coordinates, we have ι∗ |y = Oy Jp,d , where Oy ∈ O(q) and Jq,d ∈ Rq×d with 1 on
diagonal entries and 0 on others. Therefore, we have ι∗ |>
y ι∗ |y = Id×d . Moreover,
d−1
ι∗ | >
II
(ν,
ν)
=
0
for
any
ν
∈
S
.
Hence,
if
we
multiply
[∇Φ(x)]−1 ι∗ |>
y y
y to both
sides of above equation, we have
1
0
−1 2
[∇Φ(x)]−1 ι∗ |>
∇γ 0 (0)γ 0 (0) Φ(x)t2 + O(t3 ).
y (ι(z) − ι(y)) = γ (0)t + [∇Φ(x)]
2
Consequently,
(36)

(ι(z) − ι(y))> ι∗ |y [∇Φ(x)∇Φ(x)> ]−1 ι∗ |>
y (ι(z) − ι(y))
= (ι(z) − ι(y))> ι∗ |y [∇Φ(x)−1 ]> ∇Φ(x)−1 ι∗ |>
y (ι(z) − ι(y))
= γ 0 (0)> γ 0 (0)t2 + γ 0 (0)> [∇Φ(x)]−1 ∇2γ 0 (0)γ 0 (0) Φ(x)t3 + O(t4 )
= t2 + γ 0 (0)> [∇Φ(x)]−1 ∇2γ 0 (0)γ 0 (0) Φ(x)t3 + O(t4 ).

Note that the last step follows from γ 0 (0)> γ 0 (0) = 1.

COVARIANCE MATRIX, EIG AND LLE

17

Let c : [0, t] 7→ M be the unique geodesic from w to x so that c(0) = w and
c(t) = x. Similarly, we obtain
(ι(y) − ι(z))> ι∗ |z [∇Φ(w)∇Φ(w)> ]−1 ι∗ |>
z (ι(y) − ι(z))

(37)

= t2 + c0 (0)> [∇Φ(w)]−1 ∇2c0 (0)c0 (0) Φ(w)t3 + O(t4 )
= t2 − γ 0 (t)> [∇Φ(γ(t))]−1 ∇2γ 0 (t)γ 0 (t) Φ(γ(t))t3 + O(t4 )
= t2 − γ 0 (0)> [∇Φ(x)]−1 ∇2γ 0 (0)γ 0 (0) Φ(x)t3 + O(t4 ),
where the second last equality comes from w = γ(t), c0 (0) = −γ 0 (t) and the last
equality comes from the Taylor expansion of γ 0 (t)> [∇Φ(γ(t))]−1 ∇2γ 0 (t)γ 0 (t) Φ(γ(t))
at t = 0. The conclusion follows from adding equations (36) and (37) together. 

Proof of Lemma 3. Fix x ∈ M and use polar coordinates (t, θ) ∈ [0, ∞) × S d−1 to
parametrize Tx M . Set y 0 = Φ(x0 ) and x0 = expx (θt). By definition we have
Z
C (y) =

(38)

(ι(y 0 ) − ι(y))(ι(y 0 ) − ι(y))> Q(y 0 )dVN (y 0 ) ∈ Rq×q ,

E(y)

where Q is the probability density function of the random vector Y := Φ ◦ X
associated with the Riemannian measure defined on N . Let P (x0 ) = Q ◦ Φ(x0 ),
then by a direct Taylor expansion we have
Z
C (y) =

(ι ◦ Φ(x0 ) − ι ◦ Φ(x))(ι ◦ Φ(x0 ) − ι ◦ Φ(x))> P (x0 )|∇Φ(x0 )|dVM (x0 )

B (x)

Z

Z

=
S d−1



(ι∗ |y ∇Φ(x)θt + O(t2 ))(ι∗ |y ∇Φ(x)θt + O(t2 ))> (P (x) + O(t))

0

× |∇Φ(x)|(1 + O(t))(td−1 + O(td+1 ))dtdθ
Z
Z 
= |∇Φ(x)|P (x)
(ι∗ |y ∇Φ(x)θ)(ι∗ |y ∇Φ(x)θ)> td+1 + O(td+2 )dtdθ
S d−1 0
nZ
o
|∇Φ(x)|P (x)d+2
=
[ι∗ |y ∇Φ(x)]
θθ> dθ [ι∗ |y ∇Φ(x)]> + O(d+4 ).
d+2
S d−1
Note that all terms
td+2 above contain a factor θθ> θ. Due to the symmetry
R of order
>
of the sphere, S d−1 θθ θdθ = 0. Hence, the error in the last step above is of
R
d−1
order O(d+4 ) rather than O(d+3 ). Moreover, we have S d−1 θθ> dθ = |S d | Id×d ,
therefore

(39)

C (y) =

|S d−1 ||∇Φ(x)|P (x)d+2
>
>
[ι∗ |y ∇Φ(x)][∇Φ(x) ι∗ |y ] + O(d+4 ).
d(d + 2)

18

JOHN MALIK, CHAO SHEN, HAU-TIENG WU, AND NAN WU

By a similar calculation, we have
Z
E[χE (y) (Y )] =
(40)
Q(y 0 )dVN (y 0 )
E(y)
Z
0
=
P (x )|∇Φ(x0 )|dVM (x0 )
B (x)

Z

Z



(P (x) + O(t))|∇Φ(x)|(1 + O(t))(td−1 + O(td+1 ))dtdθ

=
S d−1

Z

0

Z

=
S d−1



[P (x)|∇Φ(x)|td−1 + O(td )]dtdθ

0

|S d−1 ||∇Φ(x)|P (x)d
+ O(d+2 )
=
d
d
Note that all the
R order t terms above contain a factor θ. Due to the symmetry
of the sphere, S d−1 θdθ = 0. Hence, the error in the last step above is of order
O(d+2 ) rather than O(d+1 ). The first result follows from taking the quotient of
the above two equations.
By Lemma A.4, we have

ι ◦ Φ(x0 ) − ι ◦ Φ(x) = ι∗ |y ∇θ Φ(x)t +


1
IIy (∇θ Φ(x), ∇θ Φ(x)) + ι∗ |y ∇2θθ Φ(x) t2 + O(t3 ).
2

If v1 , v2 ∈ (ι∗ Ty N )⊥ , since ι∗ |y ∇θ Φ(x) and ι∗ |y ∇2θθ Φ(x) are in ι∗ Ty N , then
1 >
v IIy (∇θ Φ(x), ∇θ Φ(x))t2 + O(t3 ),
2 1
1
v2> (ι ◦ Φ(x0 ) − ι ◦ Φ(x)) = v2> IIy (∇θ Φ(x), ∇θ Φ(x))t2 + O(t3 ),
2

v1> (ι ◦ Φ(x0 ) − ι ◦ Φ(x)) =

(41)
(42)

and we have
v1> C (y)v2
Z
=
v1> (ι ◦ Φ(x0 ) − ι ◦ Φ(x))(ι ◦ Φ(x0 ) − ι ◦ Φ(x))> v2 P (x0 )|∇Φ(x0 )|dVM (x0 )
B (x)

Z

Z

=
S d−1

0



1 >
(v IIy (∇θ Φ(x), ∇θ Φ(x))t2 + O(t3 ))(v2> IIy (∇θ Φ(x), ∇θ Φ(x))t2 + O(t3 ))>
4 1

× (P (x) + O(t))|∇Φ(x)|(1 + O(t))(td−1 + O(td+1 ))dtdθ
Z
Z 
|∇Φ(x)|P (x)
=
v1> IIy (∇θ Φ(x), ∇θ Φ(x))(IIy (∇θ Φ(x), ∇θ Φ(x)))> v2 td+3 + O(td+4 )dtdθ
4
d−1
S
0
Z
|∇Φ(x)|P (x)d+4
=
v1> IIy (∇θ Φ(x), ∇θ Φ(x))(IIy (∇θ Φ(x), ∇θ Φ(x)))> v2 dθ + O(d+6 ).
4(d + 4)
d−1
S
Due to the symmetry of the sphere, the error in the last step above is of order
O(d+6 ) rather than O(d+5 ). The result follows from taking the quotient.

To further control the influence of the truncated pseudo-inverse, we need the
following lemma to further explore the higher order eigenvalue structure of the
local covariance matrix.

COVARIANCE MATRIX, EIG AND LLE

Lemma A.5. Up to a rotation of Rq , we have
"
>
1
∇Φ(x)∇Φ(x) + O(2 )
(43)
C¯ (y) = d+2
O(2 )

19

#
O(2 )
,
d
2
4
4(d+4) M  + O( )

where M is a (q − d) × (q − d) diagonal matrix satisfying
(44)
Mi,i =

1
|S d−1 |

Z
S d−1

>
e>
q−d+i IIy (∇θ Φ(x), ∇θ Φ(x))(IIy (∇θ Φ(x), ∇θ Φ(x))) eq−d+i dθ.

Moreover, if Mi,i = 0, e>
q−d+i IIy (∇θ Φ(x), ∇θ Φ(x)) = 0.
Proof. First, we can rotate
 ι(N ) so that ι∗ Ty N is generated by {e1 , e2 , · · · , ed }. In
O1
other words, ι∗ |y =
∈ Rq×d , where O1 ∈ Od . Then, by the previous lemma
0
we have

 1
>
O ∇Φ(x)∇Φ(x) O1> + O(2 )
O(2 )
(45)
,
C¯ (y) = d+2 1
O(2 )
M̄ 2 + O(4 )
where
M̄ is a q − d by q − d symmetric matrix. If O2 diagonalizes M̄ , let O =

O1 0
, and we have
0 O2
#
"
>
1
∇Φ(x)∇Φ(x) + O(2 )
O(2 )
>¯
d+2
,
O C (y)O =
d
2
4
O(2 )
4(d+4) M  + O( )
where M is a (q − d) × (q − d) diagonal matrix. And by the previous lemma
Z
1
Mi,i = d−1
e>
IIy (∇θ Φ(x), ∇θ Φ(x))(IIy (∇θ Φ(x), ∇θ Φ(x)))> eq−d+i dθ.
|S
| S d−1 q−d+i
Pd
d−1
Note that if θ =
⊂ Ty N , then e>
∇θ Φ(x)) =
q−d+i IIy (∇θ Φ(x),
j=1 θj ∂j ∈ S
R
pi (θ1 , θ2 , · · · , θd ) where pi is a quadratic form. If Mi,i = 0, then S d−1 p2i dθ = 0.
Hence. pi = 0. The conclusion follows.

We are now ready to finish the proof of Theorem 1.
Proof of Theorem 1. Note that (ι(z)−ι(y))> Tα [C¯ (y)](ι(z)−ι(y)) is invariant under
any rotation of Rq . Therefore, by Lemma A.5 we can assume that
#
"
>
1
∇Φ(x)∇Φ(x) + O(2 )
O(2 )
d+2
¯
(46)
,
C (y) =
d
2
4
O(2 )
4(d+4) M  + O( )
where M is a (q − d) × (q − d) diagonal matrix shown in (44). Suppose that M has
l ≥ 0 zero diagonal entries. By assumption, we have p − d − l eigenvalues of order
2 that are related to principle curvatures, and l eigenvalues of order 4 . Write the
eigendecomposition of C¯ (y) as
(47)
C¯ (y) = U (x)Λ (x)U (x)> ,
where U (x) ∈ O(p) and Λ (x) is a p × p diagonal matrix. By the perturbation
argument (see Appendix A of [30] for details), Λ (x) and U (x) satisfy
 1

2
0
0
d+2 Λ1 (x) + O( )
d
2
4
0
0 ,
(48)
Λ (x) = 
4(d+4) Λ2 (x) + O( )
0
0
O(4 )

20

JOHN MALIK, CHAO SHEN, HAU-TIENG WU, AND NAN WU
>

where Λ1 (x) is the eigenvalue matrix of ∇Φ(x)∇Φ(x) , which is of order 1, and
Λ2 (x) is a (q − d − l) × (q − d − l) diagonal matrix which consists of nonzero diagonal
entries of M , which is also of order 1, and


U1 (x)
0
0
U2 (x)
0  (Ip×p + 2 S(x)) + O(4 ),
(49)
U (x) =  0
0
0
U3 (x)
>

where U1 (x) ∈ O(d) is the orthornormal eigenvector matrix of Φ(x)∇Φ(x) , U2 (x) ∈
O(q − d − l), U3 (x) ∈ O(l) and S(x) is anti-symmetric.1 We now finish the proof.
Case 1: α = d. In this case,




U (x) + O(2 )
Td [C̄ (y)] = (d + 2) 1
[Λ−1
(x) + O(2 )] U1 (x)> + O(2 ) O(2 )
2
1,
O( )


>
[∇Φ(x)∇Φ(x) ]−1 0
= (d + 2)
+ O(2 )
0
0
2
= (d + 2)ι∗ |y [∇Φ(x)∇Φ(x)> ]−1 ι∗ |>
y + O( ).

Therefore,
(50)
(ι(z) − ι(y))> Td [C¯ (y)](ι(z) − ι(y))
2
2
= (d + 2)(ι(z) − ι(y))> ι∗ |y [∇Φ(x)∇Φ(x)> ]−1 ι∗ |>
y (ι(z) − ι(y)) + O( kι(y) − ι(z)kRq ),
2 2
= (d + 2)(ι(z) − ι(y))> ι∗ |y [∇Φ(x)∇Φ(x)> ]−1 ι∗ |>
y (ι(z) − ι(y)) + O( t ) ,

where the last step follows from Lemma A.4 that O(kι(y) − ι(z)k2Rq ) = O(t2 ).
Similarly,
(51)
(ι(y) − ι(z))> Td [C¯ (z)](ι(y) − ι(z))
2 2
= (d + 2)(ι(y) − ι(z))> ι∗ |z [∇Φ(w)∇Φ(w)> ]−1 ι∗ |>
z (ι(y) − ι(z)) + O( t ).

By Lemma A.4, we have
h T [C¯ (y)] + T [C¯ (z)] i
d 
d 
(ι(z) − ι(y))
2
= (d + 2)t2 + O(t2 2 ) + O(t4 ) .

EIG2d (y, z) = (ι(z) − ι(y))>

Preparation for α 6= d. To show the case when α 6= d, we need some preparations. Let λ1 ≥ λ2 ≥ · · · λq denote the eigenvalues of C¯ (y), and let {ui } be
the corresponding
normalized eigenvectors. Note that ui is the i-th column of


U1 (x)
 0  + O(2 ) when 1 ≤ i ≤ d in (49). Based on (49), for 1 ≤ i ≤ d, the first d
0
entries of ui are of order 1 and the other entries are of order 2 . Note that the first
d entries of ui are associated with the tangent space and the others are associated
with the normal space. For d + 1 ≤ i ≤ q − l, the tangent components of ui (the
first d entries) are of order 2 and the normal components of ui (the remaining
1We mention that if the eigenvalues of Φ(x)∇Φ(x)> repeat, U (x) may not be uniquely deter1
mined, or if the nonzero diagonal terms of Λ2 (x) repeat, U2 (x) may not be identity matrix. In this
case, if one wants to fully determine U1 (x), U2 (x) and U3 (x), he/she needs to further explore the
higher order expansion of C¯ (y). Since it is irrelevant to this proof, we refer readers with interest
to Appendix A of [30].

COVARIANCE MATRIX, EIG AND LLE

21

p − d entries) are of order 1. On the other hand, based on Lemma A.4, the tangent
component of ι(z) − ι(y) is of order t and the normal component of ι(z) − ι(y) is
of order t2 . Following the notation in Lemma A.4 and putting the above together,
for 1 ≤ i ≤ p − l, we have
(ι(z) − ι(y))> ui = O(t)

(52)

>

2

for 1 ≤ i ≤ d,

2

(ι(z) − ι(y)) ui = O(t + t )

(53)

for d + 1 ≤ i ≤ q − l.
2

For q −l + 1 ≤ i ≤ q, we have ui = ūi + O( ), where ūi is (i − q − l)’s column of

0
 0  shown in (49). By Lemma A.4, ū>
i (ι(z) − ι(y)) becomes
U3 (x)
2
1 >
2
3
ūi IIy (∇θ Φ(x), ∇θ Φ(x)) + ū>
ū>
i ι∗ |y ∇θθ Φ(x) t + O(t ).
i ι∗ |y ∇θ Φ(x)t +
2
>
2
Note that ū>
i ι∗ |y ∇θ Φ(x) = 0 and ūi ι∗ |y ∇θθ Φ(x) = 0, since ūi is in the nor>
mal direction. By Lemma A.5, we have ūi IIy (∇θ Φ(x), ∇θ Φ(x)) = 0. Therefore
3
2
ū>
i (ι(z) − ι(y)) = O(t ). Since ui − ūi is of order O( ), ι(z) − ι(y) is of order O(t),
>
2
we have (ui − ūi ) (ι(z) − ι(y)) = O(t ). Putting the above together gives
(ι(z) − ι(y))> ui = O(t3 + t2 )

(54)

for q − l + 1 ≤ i ≤ q.

With (52), (53) and (54), we can finish the proof for α 6= d.
Case 2: 1 ≤ α < d. In this case, we have

(55)
(ι(z) − ι(y))> Tα [C¯ (y)] − Td [C¯ (y)] (ι(z) − ι(y))
= (ι(z) − ι(y))>

d
i
h X
uj u>
j
(ι(z) − ι(y)).
λj
j=α+1

Recall that λj is of order 1 for 1 ≤ j ≤ d and also (52). As a result,
(ι(z) − ι(y))> Tα [C¯ (y)](ι(z) − ι(y)) = (ι(z) − ι(y))> Td [C¯ (y)](ι(z) − ι(y)) + O(t2 ).
Similarly,
(ι(z) − ι(y))> Tα [C¯ (z)](ι(z) − ι(y)) = (ι(z) − ι(y))> Td [C¯ (z)](ι(z) − ι(y)) + O(t2 ).
We then achieve the claim
EIG2α (y, z) = (d + 2)t2 + O(t2 ) .

(56)

We now show a special case when 1 ≤ α < d. Let the geodesic distance between y
and z be d(y, z), then by Lemma A.4, we have d(y, z) = kι(y) − ι(z)kRq + O(kι(y) −
ι(z)k3Rq ) = O(t). In a special case when ι∗ |y ϑ(y) ∈ Vα (y), where Vα (y) is the
subspace of ι∗ Ty N generated by the first α eigenvectors of C¯ (y), we have

(57)
(ι(z) − ι(y))> Tα [C¯ (y)] − Td [C¯ (y)] (ι(z) − ι(y))
d
h X
i
uj u>
j
= (ι(z) − ι(y))
(ι(z) − ι(y)) = O(d(y, z)4 ) = O(t4 ).
λ
j
j=α+1
>

If furthermore ι∗ |z ϑ(z) ∈ Vα (z), where Vα (z) is the subspace of ι∗ Tz N generated
by the first α eigenvectors of C¯ (z), we have

(58)
(ι(z) − ι(y))> Tα [C¯ (z)] − Td [C¯ (z)] (ι(z) − ι(y)) = O(t4 )

22

JOHN MALIK, CHAO SHEN, HAU-TIENG WU, AND NAN WU

and hence
EIG2α (y, z) = (d + 2)t2 + O(t2 2 ) + O(t4 ) .

(59)

Case 3: α > d. In this case, without lose of generality, we assume 1 < l ≤ q − d
and α ≥ q − l. Similarly, we have

(ι(z) − ι(y))> Tα [C¯ (y)] − Td [C¯ (y)] (ι(z) − ι(y))
= (ι(z) − ι(y))>

q−l
h X
uj u>
j
+
λj
j=d+1

α
X
j=q−l+1

i
uj u>
j
(ι(z) − ι(y))
λj

Since λj is of order 2 for d + 1 ≤ j ≤ q − l, by (53)
(ι(z) − ι(y))

>

q−l
X
uj u>
j
(ι(z) − ι(y)) = O((t2 / + t)2 ).
λj

j=d+1

When q − l + 1 ≤ j ≤ α, λj is of order 4 by assumption. By (54)
(ι(z) − ι(y))>

α
X
j=q−l+1

uj u>
j
(ι(z) − ι(y)) = O((t3 /2 + t)2 ).
λj

Hence,

(ι(z) − ι(y))> Tα [C¯ (y)] − Td [C¯ (y)] (ι(z) − ι(y)) = O((t2 / + t)2 + (t3 /2 + t)2 ).
By the same argument, we have

(ι(z) − ι(y))> Tα [C¯ (z)] − Td [C¯ (z)] (ι(z) − ι(y)) = O((t2 / + t)2 + (t3 /2 + t)2 ).
Summing above two equations together leads to
(60)

EIG2α (y, z) = (d + 2)t2 + O((t2 / + t)2 + (t3 /2 + t)2 ) .

Since t ≤  by assumption, we set t = β , where β ≥ 1. By a straightforward
calculation, t2 dominates O((t2 / + t)2 + (t3 /2 + t)2 ) asymptotically when β > 1;
otherwise t2 and O((t2 / + t)2 + (t3 /2 + t)2 ) are asymptotically of the same order.
We thus conclude the claim.

Appendix B. Proof of Theorem 2
By Lemma 1 we have
Uh̄ (x) =


U1
0


0
+ O(h̄2 ),
U2

where U1 ∈ O(d) and U2 ∈ O(p − d). By Equation (6),
ι(y) − ι(x) = (ι∗ θ)t +

IIx (θ, θ) 2
t + O(t3 ).
2

By Equation (7),
t=h+

kIIx (θ, θ)k2 3
h + O(h4 ).
24

Hence
ι(y) − ι(x) = (ι∗ θ)h +

IIx (θ, θ) 2
h + O(h3 ).
2

COVARIANCE MATRIX, EIG AND LLE

23

Note that (ι∗ θ)h is in ι∗ Tx M . Therefore
Ph⊥ (ι(y) − ι(x)) =

IIx (θ, θ) 2
IIx (θ, θ) 2
h + O(h3 + hh̄2 ) =
h + O(hh̄2 ),
2
2

where we use the fact that h < h̄ in the last step. Hence,
kPh⊥ (ι(y) − ι(x))k2Rp
kIIx (θ, θ)k2 3
=
h + O(h2 h̄2 ).
6h
24
The conclusion follows.



Appendix C. Proof of Theorem 3
The proof of Theorem 3 consists of two steps, the variance analysis and bias
analysis. The variance analysis of the LDR-LLE is similar to Case 0 in [30, Theorem
3.1], so we only provide the result and the proof is omitted.
√
log(n)
Proposition 1. Fix f ∈ C(ι(M )). Suppose h = h(n) so that n1/2 hd/2+1 → 0 and
h → 0 as n → ∞. With probability greater than 1 − n−2 , for all xk ∈ X ,
Z
N
 plog(n) 
X
tLLE
w̃n,xk (i)f (xk,j ) =
Ph (xk , y)f (y)P (y)dV (y) + O 1/2 d/2−1 .
n h
M
j=1
The bias analysis part of Theorem 3 depends on the following two technical
lemmas. We use the following notations to simplify the proof. For p, d ∈ N so that
d ≤ p, denote Jp,d ∈ Rp×d so that the (i, i) entry is 1 for i = 1, . . . , d, and zeros
elsewhere. For v ∈ Rp ,
v = [[v1 , v2 ]] ∈ Rp ,

(61)

where v1 ∈ Rd forms the first d coordinates of v and v2 ∈ Rp−d forms the last
p − d coordinates of v. Thus, by a proper translation and rotation of ι(M ) in Rp so
that ι(x) = 0 and ι∗ Tx occupies the first d axes of Rp , for v = [[v1 , v2 ]] ∈ Tι(x) Rp ,
>
v is tangential to ι∗ Tx M and [[0, v2 ]] is normal to ι∗ Tx M .
v1 = Jp,d
Lemma A.6. ([30, Lemma B.5]) Fix x ∈ M and f ∈ C 3 (M ). When h > 0 is
sufficiently small, the following expansions hold.
(1) The scalar E[χB Rp (ι(x)) (X)] satisfies
h

|S d−1 | d
h + O(hd+2 ).
h
d
(2) The scalar E[(f (X) − f (x))χB Rp (ι(x)) (X)] satisfies
E[χB Rp (ι(x)) (X)] = P (x)
h

E[(f (X) − f (x))χ
(X)]
i
d−1 h
|S
| 1
=
P (x)∆f (x) + ∇f (x) · ∇P (x) hd+2 + O(hd+3 ).
d(d + 2) 2
Rp (ι(x))
Bh

(3) The vector E[(X − ι(x))χB Rp (ι(x)) (X)] satisfies
h

E[(X − ι(x))χ
(X)]
hh |S d−1 |
ii
>
=
Jp,d
ι∗ ∇P (x)hd+2 + O(hd+3 ), O(hd+2 ) .
d(d + 2)
Rp (ι(x))
Bh

24

JOHN MALIK, CHAO SHEN, HAU-TIENG WU, AND NAN WU

(4) The vector E[(X − ι(x))(f (X) − f (x))χB Rp (ι(x)) (X)] satisfies
h

(X)]
E[(X − ι(x))(f (X) − f (x))χ
hh |S d−1 |
ii
>
=
Jp,d
P (x)ι∗ ∇f (x)hd+2 + O(hd+3 ), O(hd+4 ) .
d(d + 2)
Rp (ι(x))
Bh

The proof of the above Lemma can be found in [30, Lemma B.5]. The next
lemma describes the vector T when h is sufficiently small. The main significance is
that the vector T̃ι(x) almost recovers ∇ log P (x) in the tangent direction.
Lemma A.7. Fix x ∈ M . When h is sufficiently small, we have
ii
hh J > ι∗ ∇P (x)
p,d
+ O(h2 ), O(h2 ) .
T̃ι(x) =
P (x)
Proof. Note that


d
ui u>
Rp (ι(x))

 X
i E(X − ι(x))χBh
T̃ι(x) := Td [Ch (x)] E(X − ι(x))χB Rp (ι(x)) =
,
h
λi
i=1


>
ei + O(h2 )
U1 Jp,d
,
where ui and λi form the i-th eigenpair of Ch (x). By Lemma 1, ui =
O(h2 )
i = 1, . . . , d and U1 ∈ O(d), and λi =
we have for 1 ≤ i ≤ d

|S d−1 |P (x) d+2
d(d+2) h

+ O(hd+4 ). By Lemma A.6 ,

u>
Rp (x ) (X)]
i E[(X − ι(x))χBh
k
ii hh
ii
d−1 hh J > ι ∇P (x)
|S
| p,d ∗
>
hd+2 + O(hd+4 ), O(hd+2 ) · U1 Jp,d
ei + O(h2 ), O(h2 )
=
d+2
d
d−1
|S
|
>
=
(ι∗ ∇P (x))> Jp,d U1 Jp,d
ei hd+2 + O(hd+4 ) .
d(d + 2)
Thus, for 1 ≤ i ≤ d,
p
u>
i E[(X − ι(x))χB R (xk ) (X)]
h

λi

=

|S d−1 |
>
>
d+2
+
d(d+2) (ι∗ ∇P (x)) Jp,d U1 Jp,d ei h
|S d−1 |P (x) d+2
+ O(hd+4 )
d(d+2) h

=

>
(ι∗ ∇P (x))> Jp,d U1 Jp,d
ei
+ O(h2 ) .
P (x)

O(hd+4 )

Hence, we have
T̃ι(x) =

d
p
X
u>
i E[(X − ι(x))χB R (xk ) (X)]
h

i=1

=

λi

d 
>
X
(ι∗ ∇P (x))> Jp,d U1 Jp,d
ei
i=1

P (x)

ui

hh
ii
>
+ O(h2 ) U1 Jp,d
ei + O(h2 ), O(h2 )

hh J > ι∗ ∇P (x)
ii
p,d
=
+ O(h2 ), O(h2 ) ,
P (x)
where the last step follows from the fact that U1 ∈ O(d).
We are ready to prove Theorem 3.



COVARIANCE MATRIX, EIG AND LLE

25

Proof of Theorem 3. Note that
Z
PtLLE (x, y)f (y)P (y)dV (y) − f (x)
M

=

p
E[(f (X) − f (x))χB Rp (ι(x)) (X)] − T̃>
ι(x) E[(X − ι(x))(f (X) − f (x))χB R (ι(x)) (X)]
h

h

p
E[χB Rp (ι(x)) (X)] − T̃>
ι(x) E[(X − ι(x))χB R (ι(x)) (X)]

.

h

h

By Lemma A.6 and Lemma A.7, we have
d+2
p
),
T̃>
ι(x) E[(X − ι(x))χB R (ι(x)) (X)] = O(h

(62)

h

and
p
T̃>
ι(x) E[(X − ι(x))(f (X) − f (x))χB R (ι(x)) (X)] =
h

|S d−1 |
∇P (x) · ∇f (x)hd+2 + O(hd+3 ).
d(d + 2)

Hence,
p
E[χB Rp (ι(x)) (X)] − T>
ι(x) E[(X − ι(x))χB R (ι(x)) (X)] = P (x)
h

h

|S d−1 | d
h + O(hd+2 ),
d

and
E[(f (X) − f (x))χB Rp (ι(x)) (X)] − T>
Rp (ι(x)) (X)]
ι(x) E[(X − ι(x))(f (X) − f (x))χBh
h
i
|S d−1 |
|S d−1 | h 1
P (x)∆f (x) + ∇f (x) · ∇P (x) hd+2 −
∇P (x) · ∇f (x)hd+2 + O(hd+3 )
=
d(d + 2) 2
d(d + 2)
|S d−1 |
=
P (x)∆f (x)hd+2 + O(hd+3 ).
2d(d + 2)
If we combine the above two equations, we conclude the claim that
Z
1
∆f (x)h2 + O(h3 ).
PtLLE (x, y)f (y)P (y)dV (y) − f (x) =
2(d + 2)
M

Department of Mathematics, Duke University, Durham, NC, USA
E-mail address: john.malik@duke.edu
Department of Mathematics, Duke University, Durham, NC, USA
E-mail address: chao.shen@duke.edu
Department of Mathematics and Department of Statistical Science, Duke University, Durham, NC, USA; Mathematics Division, National Center for Theoretical Sciences, Taipei, Taiwan
E-mail address: hauwu@math.duke.edu
Department of Mathematics, University of Toronto, Toronto, ON, Canada
E-mail address: n.wu@mail.utoronto.ca

