The Shape of Data and Probability Measures
Diego H. Dı́az Martı́neza , Facundo Mémolib , Washington Mioa

arXiv:1509.04632v2 [stat.ML] 28 Feb 2017

a

Department of Mathematics, Florida State University, Tallahassee, FL 32306-4510 USA
b
Department of Mathematics, Ohio State University, Columbus, OH 43210-1174 USA

Abstract
We introduce the notion of multiscale covariance tensor fields (CTF) associated with Euclidean random variables as a gateway to the shape of their
distributions. Multiscale CTFs quantify variation of the data about every
point in the data landscape at all spatial scales, unlike the usual covariance tensor that only quantifies global variation about the mean. Empirical
forms of localized covariance previously have been used in data analysis and
visualization, for example, in local principal component analysis, but we develop a framework for the systematic treatment of theoretical questions and
mathematical analysis of computational models. We prove strong stability
theorems with respect to the Wasserstein distance between probability measures, obtain consistency results for estimators, as well as bounds on the
rate of convergence of empirical CTFs. These results show that CTFs are
robust to sampling, noise and outliers. We provide numerous illustrations of
how CTFs let us extract shape from data and also apply CTFs to manifold
clustering, the problem of categorizing data points according to their noisy
membership in a collection of possibly intersecting smooth submanifolds of
Euclidean space. We prove that the proposed manifold clustering method is
stable and carry out several experiments to illustrate the method.
Keywords: shape of data, multiscale data analysis, covariance fields,
Fréchet functions, manifold clustering
1. Introduction
Probing, analyzing and visualizing the shape of complex data are challenges that are magnified by the intricate dependence of their structural
properties, as basic as dimensionality, on location and scale (cf. [1]). As such,
resolving and integrating the geometry and topology of data across scales are
1

problems of foremost importance. In this paper, we develop the notion of
multiscale covariance tensor fields (CTF) associated with Euclidean random
variables and show that many properties of the shape of their distributions
become accessible through CTFs, which provide stable representations that
can be estimated reliably from data.
For a random vector y ∈ Rd , scale dependence is controlled by a kernel
function K(x, y, σ) > 0, where x, y ∈ Rd and σ > 0 is the scale parameter.
The idea is that from the standpoint of x, at scale σ > 0, the kernel masks
the distribution by attributing weight K(x, y, σ) to data located at y, creating a windowing effect. More simply put, K(x, y, σ) quantifies how well
an observer at x sees data at y at scale σ. Covariation of the weighted data
is measured relative to every point x ∈ Rd , not just about the mean as is
common practice, thus giving rise to a multiscale covariance field. Special
cases of these covariance fields were introduced in [2], targeting applications
to such problems as detection of local scales and feature rich points in shapes.
Here we present a more systematic treatment that includes a broader formulation of multiscale CTFs, stability theorems that ensure that properties of
probability measures derived from multiscale CTFs are robust, as well as
consistency results and convergence rates for empirical CTFs. We prove stability of CTFs with respect to the Wasserstein distance between probability
measures, a metric that is finding uses in an ever expanding landscape of
problems and whose origins are in optimal transport theory [3, 4]. Since
Wasserstein distance metrizes weak convergence of probability measures, we
obtain a strong stability result that ensures that if two probability distributions are similar in a weak sense, then their multiscale CTFs are uniformly
close over the entire domain. Convergence rates are derived from the stability theorems and results by Fournier and Guillin [5] and Garcı́a-Trillos
and Slepčev [6] on convergence of empirical measures. The standard covariance tensor of a random vector y ∈ Rd quantifies covariation of y about the
mean, but may be extended to a full covariance field by considering covariation about arbitrary points. Nonetheless, this field provides no information
about the organization of the data other than that already contained in the
covariance about the mean. Thus, a localized formulation is essential for
gaining additional insight into the shape of data.
The trace of a multiscale CTF is a scalar field that gives a multiscale
analogue of the classical Fréchet function V (x) = E [ky − xk2 ] of a random
variable y with finite second moment. The Fréchet function provides a more
geometric interpretation of the mean as the unique minimizer of V ; that is,
2

the point µ ∈ Rd with respect to which the spread of y is minimal. Similarly,
the local extrema and other properties of the multiscale Fréchet function
provide a wealth of information about the distribution of y. In fact, we show
that the distribution of any random vector may be fully recovered from the
multiscale Fréchet function associated with the Gaussian kernel.
Several variants of empirical localized or weighted covariance previously
have been used in data analysis, but we develop a framework for the formulation and systematic treatment of such problems. Allard et al. have
developed a computational model termed geometric multi-resolution analysis for multiscale data analysis based on covariance localized to hierarchies
of dyadic cubes [7]. In computer graphics, local principal component analysis (PCA) is commonly used in the estimation of normals to surfaces from
point-cloud data [8] for surface reconstruction; see also [9] and references
therein. In computer vision, tensor voting by Medioni et al. [10] has been
applied to multiple image analysis and processing tasks. Brox et al. have
used empirical covariance weighted by the isotropic Gaussian kernel in nonparametric density estimation targeting applications in motion tracking [11].
In the literature dealing with clustering, especially clustering of multiple possibly intersecting manifolds, local PCA ideas have been used in the works of
Kushnir et al. [12], Goldberg et al. [13], Gong et al. [14], Wang et al. [15],
and in a series of papers by Arias-Castro and collaborators [16, 17, 18].

(a)

(b)

(c)

Figure 1: Examples of data clustered along intersecting manifolds.

The special case of affine linear subspaces, known as subspace clustering,
has been addressed in the machine learning and computer vision literature by
many authors using a variety of techniques (cf. [19, 20, 21, 22, 23, 24, 25, 26]).
More general manifold clustering has been considered in [27, 28, 29, 12, 30,
13, 14, 15, 18, 31]. In our approach, we exploit the fact that localized covariance tensors encode rich information about the tangential structure of the
submanifolds that underlie the data. Combined with information about the
3

(relative) positions of the data points, they yield an effective data representation for manifold clustering. Although several different clustering techniques
could be applied to the “tensorized” data, we use the single linkage hierarchical method because it produces provenly stable dendrograms. In conjunction
with the stability and consistency results for covariance fields, this ensures
that the manifold clustering method is stable at all steps. Dendrogram stability is analyzed in the framework of [32].
Contributions and Organization of the paper. The paper includes several illustrations and applications of CTFs to data analysis. For example, to illustrate how geometric information can be extracted from CTFs, we show
that the curvature of plane curves and the principal curvatures of surfaces
in R3 can be calculated from the spectrum of multiscale CTFs. Thus, multiscale covariance tensors give a way of extending these infinitesimal measures
of geometric complexity to all scales and general probability distributions,
not just those supported on smooth submanifolds. We also apply multiscale
CTFs to manifold clustering, the problem of clustering Euclidean data that
are organized along a finite union of possibly intersecting smooth submanifolds. Fig. 1 shows three such examples.
The main goals of the paper are: (i) to establish the foundations for analysis, visualization and management of data with methods based on multiscale
covariance tensor fields, and (ii) to describe applications that characterize
the usefulness of CTFs in data analysis. In Section 2, we formulate the notion of multiscale CTFs for a broad class of kernels and give examples that
illustrate how CTFs reveal the geometry of data. In Section 3, we show that
the curvature of a plane curve and the principal curvatures of a surface in
R3 can be recovered from small-scale covariance. Section 4 is devoted to the
main theoretical developments. We prove stability and consistency theorems
for multiscale covariance tensor fields under mild regularity assumptions on
the kernel, and also analyze rates of convergence that are important for applications in data analysis. Since some discontinuous kernels are of practical
interest, we also investigate convergence results for such kernels, including a
pointwise central limit theorem. Multiscale Fréchet functions are discussed
in Section 5 and manifold clustering in Section 6. We close with a summary
and some discussion.

4

2. Covariance Tensor Fields
2.1. Preliminaries
To define covariance tensor fields, we introduce some notation. Elements
of the tensor product Rd ⊗ Rd may be identified with bilinear forms B : Rd ×
Rd → R through the Euclidean inner product. More precisely, a pure tensor
x ⊗ y corresponds to the bilinear form
x ⊗ y (u, v) = hx, ui · hy, vi ,

(1)

∀u, v ∈ Rd , where h , i denotes Euclidean inner product. Bilinear forms associated with more general elements of Rd ⊗ Rd can be described by linear
extension. In Euclidean coordinates, we abuse notation and also write the
coordinate vectors of x, y ∈ Rd as x and y. With this convention, letting A
be the d × d matrix A = xy T , we have
x ⊗ y (u, v) = hu, Avi ,

(2)

where the superscript T denotes transposition. In this manner, using Euclidean coordinates, an element of Rd ⊗ Rd also can be identified with a d × d
matrix by linear extension of the correspondence x ⊗ y ↔ A. Through these
identifications, we refer to an element Σ ∈ Rd ⊗ Rd interchangeably as a
tensor, a bilinear form or a matrix. We equip Rd ⊗ Rd with the inner product
defined on pure tensors by
hx1 ⊗ y1 , x2 ⊗ y2 i = hx1 , x2 i hy1 , y2 i

(3)

and extended linearly to Rd ⊗ Rd . Thus, the corresponding norm satisfies
kx ⊗ yk = kxkkyk ,

(4)

for any x, y ∈ Rd . In matrix representation, this is the Frobenius norm.
Throughout the paper, we view Rd as a measurable space equipped with
the Borel σ-algebra for the Euclidean metric. Let y be an Rd -valued random
variable distributed according to the probability measure α. Suppose that y
has expected value E [y] = µ ∈ Rd and finite second moment. As a motivation
for the definition of multiscale CTFs, recall that the covariance tensor of y
is defined as
Z
Σα (µ) = E [(y − µ) ⊗ (y − µ)] =
(y − µ) ⊗ (y − µ) α(dy) ∈ Rd ⊗ Rd . (5)
Rd

5

In matrix notation,
Z
Σα (µ) =

(y − µ)(y − µ)T α(dy) .

(6)

Rd

The bilinear form associated with Σα (µ) clearly is symmetric and positive
semi-definite.
Covariation of y may be measured with respect to any x ∈ Rd , not just µ.
Thus, Σα (µ) may be extended to a global covariance tensor field Σα : Rd →
Rd ⊗ Rd given by
Z
(y − x) ⊗ (y − x) α(dy) .
(7)
Σα (x) =
Rd

Note, however, that
Σα (x) = Σα (µ) + (µ − x) ⊗ (µ − x) ,

(8)

for any x ∈ Rd . Thus, for x 6= µ, Σα (x) does not reveal any information
about the distribution of y other than that already contained in Σα (µ). In
contrast, as we shall see below, multiscale analogues are rich in information
about the shape of α.
2.2. Multiscale Covariance Tensor Fields
We adopt the notation νd for the volume of the unit ball in Rd and ωd−1
for the “surface area” of the unit sphere Sd−1 ⊂ Rd , d ≥ 1. Recall that
ωd−1 = 2π d/2 /Γ(d/2), where Γ(·) is the Gamma function, and ωd−1 = d νd .
We make the convention that ν0 = 1.
Let y be an Rd -valued random variable with distribution α and let K be
a multiscale kernel; that is, a measurable function K : Rd × Rd × (0, ∞) → R
such that K(x, y, σ) > 0, for any x, y ∈ Rd and σ > 0.
Definition 1. The multiscale covariance tensor field (CTF) of y associated
with the kernel K is the one-parameter family of tensor fields, indexed by
σ ∈ (0, ∞), given by
Z
Σα (x, σ) :=
(y − x) ⊗ (y − x)K(x, y, σ) α(dy) ,
(9)
Rd

provided that the integral converges for each x ∈ Rd and σ > 0.
6

Remark 1. Note that Σα depends only on the probability measure α, not
on y. For this reason, we refer to Σα interchangeably as the multiscale CTF
of the random variable y or the probability measure α.
Σα (x, σ) measures the covariation of y about x with probability mass
at y weighted by K(x, y, σ). It is simple to verify that the bilinear form
Σα (x, σ) is symmetric and positive semi-definite. Note that if K is bounded
for each σ > 0, that is, ∃Mσ > 0 such that K(x, y, σ) ≤ Mσ , ∀x, y ∈ Rd ,
then Σα (x, σ) is well defined for any random variable y with finite second
moment. In particular if K ≡ 1, Σα (x, σ) = Σα (x), ∀x ∈ Rd . However, as
our primary goal is to study the organization of data and random variables
at scales ranging from local to global, we consider kernels in Rd that satisfy
additional decay conditions as they produce a windowing effect. The kernels
are constructed as follows.
Definition 2. Let d be a positive integer and f : [0, ∞) → R a bounded and
measurable function satisfying:
(a) f (r) > 0, ∀r ∈ [0, ∞);
R∞ d
(b) Md = 0 r 2 −1 f (r) dr < ∞;
(c) There is C > 0 such that rf (r) ≤ C, ∀r ∈ [0, ∞).
The multiscale kernel K : Rd × Rd × (0, ∞) → R associated with f is defined
as


ky − xk2
1
f
,
(10)
K(x, y, σ) :=
Cd (σ)
σ2
where Cd (σ) = 21 σ d Md ωd−1 .
Condition (b) in the definition implies that the normalizing
constant
R
Cd (σ) is well defined. The normalization is adopted so that K(x, y, σ) dy =
1, ∀x ∈ Rd and ∀σ > 0. Condition (c) guarantees that the integral in (9) is
convergent for any probability measure α. Henceforth, for convenience, we
assume that sup f = 1. This is not restrictive since scaling f does not change
the kernel K because of the normalization.
Whereas we investigate properties of multiscale CTFs in a more general
setting, our examples and experiments focus on two special kernels:

7

(i) The isotropic Gaussian kernel


1
ky − xk2
,
G(x, y, σ) =
exp −
(2πσ 2 )d/2
2σ 2

(11)

which is associated with the function f (x) = e−x/2 ;
(ii) The truncation kernel
1
T (x, y, σ) = d χ
σ νd



ky − xk2
σ2


(12)

associated with the characteristic function χ : [0, ∞) → R of the unit
interval [0, 1]. In measuring covariation of random variables about x,
the kernel T attributes a uniform weight to mass at points within the
closed ball of radius σ centered at x and weight zero to mass elsewhere.
Remark 2. The kernel K defined in (10) is homogeneous and isotropic; that
is, for any isometry ϕ : Rd → Rd , K(ϕ(x), ϕ(y), σ) = K(x, y, σ), ∀x, y ∈ Rd
and σ > 0. Moreover, if we write ϕ(x) = U x + b, with U ∈ O(d) and b ∈ Rd ,
then
U Σα (x, σ) U T = Σϕ∗ (α) (ϕ(x), σ),
(13)
for any (x, σ) ∈ Rd × (0, ∞). Here O(d) is the group of d × d orthogonal
matrices and ϕ∗ (α) is the pushforward of α under ϕ.
Remark 3. Multiscale covariance tensor fields can be defined for any positive
Borel measure α that satisfies
Z
kzk2 f (kzk2 ) α(dz) < ∞ ,
(14)
Rd

not just for probability measures. In particular, if f has compact support,
covariance fields are defined for any locally finite Borel measure α; that is,
measures for which every point p ∈ Rd has an open neighborhood Up such
that α(Up ) < ∞.
We conclude this section with examples that support our contention that
multiscale covariance tensor fields are rich in information about the shape of
data.

8

Example 1. This example shows that the spectrum of multiscale covariance
tensors allow us to estimate the dimensionality of data in a scale dependent
manner. We consider the data points y1 , . . . , yn in R2 , shown in Figure
2, and calculate Σαn centered at one of the data points for the Gaussian
kernel
Pnat scales σ = 0.1 and σ = 2. Here αn denotes the empirical measure
−1
n
i=1 δyi . The covariance tensors are depicted as ellipses whose principal
axes are in the direction of the eigenvectors
√ of the covariance matrix and
√
principal radii are proportional to λ1 and λ2 , where 0 ≤ λ1 ≤ λ2 are the
eigenvalues of the covariance. At scale σ = 0.1, λ1 /λ2 = 0.908, showing that
the covariance tensor is nearly isotropic, indicating that the “dimension” of
the data is 2. At σ = 2, the ratio of the eigenvalues is 0.025, giving a highly
anisotropic covariance tensor, from which we infer that the dimension is 1.

σ = 0.1

σ=2

Figure 2: Estimating data dimensionality at different scales through multiscale covariance.

Example 2 (A linear subspace of Rd ). Let v1 , . . . , vr ∈ Rd , 1 ≤ r ≤ d,
be orthonormal vectors and consider the subspace H =< v1 , . . . , vr > that
they span. Let α denote the singular measure supported on H induced by
the volume form on H. The measure α clearly is locally finite. We calculate
multiscale covariance fields at points x ∈ H to show that H may be recovered
from Σα (x, σ). By Remark 2, we may assume that x = 0. A calculation shows
that for the Gaussian kernel,
r

Σα (0, σ) = √
( 2π)d−r σ d−r−2

r
X

vi ⊗ vi .

(15)

i=1

For the truncation kernel,
Σα (0, σ) = λr

r
X

vi ⊗ vi ,

(16)

i=1

where
1

νr−1
λr = d−r−2
σ
νd

Z

π/2

−π/2

9

sin2 θ cosr θ dθ .

(17)

For r = 1, this expression simplifies to λ1 = 2/(3σ d−3 νd ). Thus, for both
kernels, the orthogonal complement of H is the null space of Σα (0, σ) and H
is the eigenspace associated with the positive eigenvalue λr .
Example 3 (Wedge of n segments). Consider the wedge (one-point union) W
of n segments L1 , . . . , Ln in Rd attached at the origin, as depicted in Fig. 3.
Each segment Li is determined by its length `i > 0 and a unit direction vector
vi . We assume that vi 6= vj , for any 1 ≤ i < j ≤ n. Let α be the singular
measure on Rd that is supported on W and agrees with the measure induced
by arc length on each segment Li . We consider the multiscale covariance
field of α associated with the truncation kernel. For x ∈ Li , x 6= 0, as in
the case r = 1 in Example 2, we have that Σα (x, σ) = (2/3σ d−3 νd ) vi ⊗ vi at
small enough scales. Thus, Σα (x, σ) has rank one. However, at the origin,
n
1 X
(min(σ, `i ))3 vi ⊗ vi .
Σα (0, σ) = d
3σ νd i=1

(18)

for any σ > 0. Thus, for σ ≤ min{`i , 1 ≤ i ≤ n},
Σα (0, σ) =

1

n
X

3σ d−3 νd

i=1

vi ⊗ vi .

(19)

Figure 3: Covariance at the one-point union of line segments.

3. Geometry of Curves and Surfaces
In this section, we show how multiscale CTFs associated with the truncation kernel extract precise local geometric information from plane curves
and surfaces in R3 .
10

3.1. Plane Curves
Example 4. We begin with the special case of a circle. Let CR ⊂ R2 be the
circle of radius R centered at the origin in R2 and α the singular measure
supported on CR induced by arc length. For any x ∈ R2 , we denote r = kxk.
If x is such that |r − R| > σ then Σα (x, σ) = 0. Assume that x ∈ R2 and
0 < σ < R are such that r ∈ [R − σ, R + σ]. In this case, in the coordinate
system given by the directions n = x/kxk and t = n⊥ , a calculation shows
that Σα (x, σ) is diagonal with entries


1 
2
2
2
+
R
(R
cos
φ
−
4r)
sin
φ
Rφ
R
+
2r
λn (x, σ) =
πσ 2
(20)
R3
λt (x, σ) =
(φ − sin φ cos φ) ,
πσ 2
 2 2 2
−σ
where φ = arccos R +r
. Thus, the normal and tangential vectors, n
2rR
and t, are eigenvectors with eigenvalues λn and λt , respectively. Fig. 4 shows
the eigenvalues as functions of r, 0.9 ≤ r ≤ 1.1, for σ = 0.1.

Figure 4: Tangential (blue) and normal (red) eigenvalues as a function of r, 0.9 ≤ r ≤ 1.1,
at σ = 0.1, of the multiscale CTF associated with the truncation kernel for the singular
measure induced by arc length, supported on the unit circle in R2 .

Now we consider a general smooth curve C ⊂ R2 , that is, a 1-dimensional,
smooth, properly embedded submanifold of R2 . Let α be the singular measure on R2 supported on C and induced by arc length. This measure is
locally finite because the embedding is proper. We calculate the small-scale
covariance at points on C for the truncation kernel and show that the curvature can be recovered from the eigenvalues of Σα . Let x ∈ C be fixed. The
arc-length parametrization of C near x may be written as
X(s) = s −

κ2 s3
κs2 κs s3
+ O(s4 ) and Y (s) =
+
+ O(s4 ) ,
6
2
6
11

(21)

where X(s) and Y (s) are coordinates along the tangent and normal to C at
x, respectively [33]. Here, the curvature κ and its derivative κs are evaluated
at x. A calculation yields:
Proposition 1. Let σ > 0 be small. If C is a smooth plane curve and x ∈ C,
then in the coordinates specified above we have

 2σ κ2 σ3
κs σ 3
− 20π + O(σ 4 )
+ O(σ 4 )
3π
15π
.
(22)
Σα (x, σ) =
κs σ 3
κ2 σ 3
+ O(σ 4 )
+ O(σ 4 )
15π
10π
Proposition 1 implies that, for σ > 0 small, the eigenvalues of Σα are
λ1 =

2σ κ2 σ 3
κ2 σ 3
−
+ O(σ 4 ) and λ2 =
+ O(σ 4 ) ,
3π
20π
10π

(23)

so that

2σ κ2 σ 3
+
+ O(σ 4 ) .
3π
20π
Thus, the curvature at x ∈ C may be recovered, up to a sign, as
√
1/2

20π
2σ
.
κ = ± lim 3/2 tr Σα (x, σ) −
σ→0 σ
3π
tr Σα (x, σ) =

(24)

(25)

3.2. Surfaces in R3
Example 5. Let SR be the sphere of radius R centered at the origin in R3 .
For x ∈ R3 , we let r = kxk. If x is such that |r − R| > σ, then Σα (x, σ) = 0.
Assume that x 6= (0, 0, 0) and σ > 0 are such that r ∈ [R − σ, R + σ]. In
the coordinate system given by the vector n = x/kxk, and any orthonormal
basis {t1 , t2 } of the orthogonal complement n⊥ , a direct calculation shows
that Σα (x, σ) is a 3 × 3 diagonal matrix with entries
 
φ
R4
4
(cos φ + 2)
λt1 (x, σ) = λt2 (x, σ) = 3 sin
σ
2

R2
(26)
λn (x, σ) = 3 (1 − cos φ) R2 + R cos φ(R cos φ + R − 3r)
2σ

R2
+ 3 (1 − cos φ) −3Rr + 3r2
2σ
 2 2 2
−σ
. In particular, this means that λn is the einwhere φ = arccos R +r
2Rr
genvalue corresponding to the eigenvector n along the normal direction to
12

Figure 5: Tangential (blue) and normal eigenvalues (red) of the multiscale CTF associated
with the truncation kernel as a function of r, 0.9 ≤ r ≤ 1.1, at σ = 0.1, for the singular
measure induced by surface area, supported on the unit sphere in R3 .

the sphere at x/kxk, and {t1 , t2 } span the eigenspace along the tangent directions with eigenvalue λt1 = λt2 . Fig. 5 shows a plot of the eigenvalues as
a function of r, 0.9 ≤ r ≤ 1, 1, for σ = 0.1.
Now we consider a general smooth compact surface S ⊂ R3 . Let α be the
singular measure on R3 supported on S and induced by the area measure on
S. We calculate the small-scale covariance at points on S for the truncation
kernel and show that the principal curvatures may indeed be recovered from
the spectrum of Σ. Given a non-umbilic point p ∈ S, one can choose a
Cartesian coordinate system centered at p so that the x-axis is along the
direction of maximal curvature at p, the y-axis is along the direction of
minimal curvature at p, and the z-axis is along the normal to S at p.
Proposition 2. Let σ > 0 be small, p ∈ S be non-umbilic, and α be the
surface area measure on S. In the coordinate system described above, the
covariance tensor for the truncation kernel is given by


At1 O(σ 4 ) O(σ 5 )
Σα (p, σ) = O(σ 4 ) At2 O(σ 5 ) ,
(27)
5
5
O(σ ) O(σ ) An
where

13

3σ
1
+
(−3κ21 − 6κ1 κ2 + κ22 )σ 3 + O(σ 4 ),
16 256
3σ
1
At2 =
+
(κ2 − 6κ1 κ2 − 3κ22 )σ 3 + O(σ 4 ),
16 256 1
3κ21 + 2κ1 κ2 + 3κ22 3
An =
σ + O(σ 4 ),
128

At1 =

(28)

and κ1 > κ2 are the principal curvatures of S at p.
It follows from this result that, for σ > 0 small,
1
3
σ + (κ1 − κ2 )2 σ 3 + O(σ 4 ) and
16
64
 π 2 11
det Σα (p, σ) = 3κ21 + 2κ1 κ2 + 3κ22
σ + O(σ 12 ) .
2048
tr Σα (p, σ) =

(29)

As a consequence, κ1 and κ2 can be recovered from the spectrum of Σα (p, σ)
as a function of σ. Indeed, from the small scale asymptotics of the trace
and determinant of Σα (p, σ), we can extract the values of (κ1 − κ2 )2 and
3κ21 + 2κ1 κ2 + 3κ22 from which we can determine the values of κ1 and κ2 .
Proof of Proposition 2. Using cylindrical coordinates in the chosen reference
system, we can parametrize the patch S ∩ B(p, σ) as (ρ cos φ, ρ sin φ, z(ρ, φ)),
2
for φ ∈ [0, 2π], ρ ∈ [0, ρσ (φ)], where ρσ (φ) = σ− 81 κ1 (cos φ)2 +κ2 (sin φ)2 σ 3 +

2
O(σ 4 ), and z(ρ, φ) = ρ2 κ1 (cos φ)2 + κ2 (sin φ)2 + O(σ 3 ). The area element
on the patch is given by


ρ3 2
2
2
2
5
(30)
dA = ρ + (κ1 (cos φ) + κ2 (sin φ) ) + O(ρ ) dρ dφ.
2
Now we have all the ingredients neededRR
to compute Σα (p, σ). For example,
to calculate the (1, 1)-entry, we express S∩B(0,σ) x2 dA as
2π

ρσ (φ)



ρ5 cos2 φ 2 2
2
2
6
ρ cos φ +
κ2 sin φ + κ1 cos φ + O(ρ ) dρ dφ ,
2
0
0
(31)
which after a simple but tedious calculation yields the desired result. The
computation of other entries of the matrix follows similar steps.
Z

Z



3

2

14

4. Stability and Consistency
For each p ∈ [1, ∞), let Pp (Rd ) denote the collection
Borel probaR of all
d
p
bility measures α on R whose pth moment Mp (α) = kzk α(dz) is finite.
1/p
We adopt the notation mp (α) = Mp (α). For p = ∞, we let P∞ (Rd ) be the
collection of all Borel probability measures on Rd with bounded support and
m∞ (α) = sup{kzk, z ∈ supp [α]}. By Jensen’s inequality, if 1 ≤ q ≤ p ≤ ∞,
then Pp (Rd ) ⊂ Pq (Rd ) and mq (α) ≤ mp (α), for any α ∈ Pp (Rd ).
Definition 3. For p ∈ [1, ∞] and λ > 0, we define Ppλ (Rd ) ⊂ Pp (Rd ) as the
subset of all α ∈ Pp (Rd ) such that α(A) ≤ λ L(A), for all measurable sets A,
where L stands for Lebesgue measure.
Example 6. If α ∈ Pp (Rd ) is absolutely continuous with respect to the
Lebesgue measure with density function f ∈ L∞ (Rd ) satisfying kf k∞ ≤ λ,
then α ∈ Ppλ (Rd ).
Let us recall the definition of the p-Wasserstein distance Wp (α, β) between
α, β ∈ Pp (Rd ). Let Γ(α, β) be the collection of all couplings of α and β; that
is, probability measures µ on Rd × Rd such that (π1 )∗ µ = α and (π2 )∗ µ = β,
where π1 , π2 : Rd × Rd → Rd denote projections onto the first and second
components, respectively.
Definition 4. For p ∈ [1, ∞), the p-Wasserstein distance between α, β ∈
Pp (Rd ) is given by
Z Z
1/p
p
Wp (α, β) := inf
kz1 − z2 k µ(dz1 × dz2 )
,
µ∈Γ(α,β)

and the ∞-Wasserstein distance between α, β ∈ P∞ (Rd ) by
W∞ (α, β) :=

inf
µ∈Γ(α,β)

sup {kz1 − z2 k, (z1 , z2 ) ∈ supp [µ]} .

Remark 4.
(i) For any α, β ∈ Pp (Rd ), p ∈ [1, ∞], there exists a coupling that realizes
the infimum in the definition of Wp (α, β) (cf. [34]).
(ii) It is a standard result that, for each p ∈ [1, ∞), Wp defines a metric
on Pp (Rd ) that is compatible with weak convergence of probability
measures [3].
(iii) If ϕ : Rd → Rd is an isometry, then Wp (α, β) = Wp (ϕ∗ (α), ϕ∗ (β)), for
any α, β ∈ Pp (Rd ).
15

4.1. Smooth Kernels
Theorem 1 (Stability for Smooth Kernels). Let f : [0, ∞) → R be as in
Definition 2 with multiscale kernel K. Suppose that f is differentiable and
there exists a constant A1 > 0 such that r3/2 |f 0 (r)| ≤ A1 , ∀r ≥ 0. Then,
there is a constant Af > 0, that depends only on f , such that
sup kΣα (x, σ) − Σβ (x, σ)k ≤
x∈Rd

σAf
W1 (α, β),
Cd (σ)

for any α, β ∈ P1 (Rd ) and any σ > 0. Here k · k is the norm associated with
the inner product defined in (3).
Theorem 1 shows that multiscale covariance fields yield a robust representation of probability measures that make their geometric properties more
readily accessible, as illustrated in our examples. In Section 5, we show that
not only is Σα (·, σ) stable, but all the information contained in the probability measure α is fully absorbed into the multiscale CTF associated with the
Gaussian kernel. In fact, α may be recovered from the multiscale scalar field
given by Vα (x, σ) = tr Σα (x, σ), x ∈ Rd and σ > 0.
The following lemma will be used in the proof of the stability theorem
for smooth kernels. To simplify notation, we define Kσ : Rd → R by


kzk2
1
f
(32)
Kσ (z) = K(z, 0, σ) =
Cd (σ)
σ2
and Qσ : Rd → Rd ⊗ Rd by
Qσ (z) = (z ⊗ z)Kσ (z) .

(33)

Lemma 1. Let f be as in Definition 2 and suppose that f is differentiable
and there is a constant A1 > 0 such that r3/2 |f 0 (r)| ≤ A1 , ∀r ≥ 0. Then,
there is a constant Af > 0, that depends only on f , such that
kQσ (z1 ) − Qσ (z2 )k ≤

Af σ
kz1 − z2 k,
Cd (σ)

for any z1 , z2 ∈ Rd and all σ > 0.
Proof. Let z(t) = tz1 + (1 − t)z2 , 0 ≤ t ≤ 1. Then,
Z 1
d
Qσ (z1 ) − Qσ (z2 ) =
Qσ (z(t)) dt .
0 dt
16

(34)

Since




d
(z1 − z2 ) ⊗ z
kzk2
z ⊗ (z1 − z2 )
kzk2
Qσ (z) =
f
+
f
dt
Cd (σ)
σ2
Cd (σ)
σ2


kzk2
2
+ (z ⊗ z) 2
f0
(z · (z1 − z2 )) ,
σ Cd (σ)
σ2

(35)

it follows that


d
2kzk
kzk2
Qσ (z) ≤
f
kz1 − z2 k
dt
Cd (σ)
σ2


2kzk3
kzk2
0
+ 2
f
kz1 − z2 k
σ Cd (σ)
σ2


kzk2
2σ kzk
=
f
kz1 − z2 k
Cd (σ) σ
σ2


2σ kzk3 0 kzk2
f
kz1 − z2 k .
+
Cd (σ) σ 3
σ2

(36)

Since f is smooth, condition
√ (c) of Definition 2 ensures that there is a constant A2 > 0 such that rf (r) ≤ A2 , ∀r ≥ 0. Moreover, by hypothesis
r3/2 |f 0 (r)| ≤ A1 . Thus, (36) implies that
σAf
d
Qσ (z) ≤
kz2 − z1 k ,
dt
Cd (σ)

(37)

where Af = 2(A1 + A2 ). The lemma follows from (34) and (37).
Proof of Theorem 1. Without loss of generality, we may assume that x = 0.
We express the covariance fields as
Z
Z
Σα (x, σ) = Qσ (z1 ) α(dz1 ) and Σβ (x, σ) = Qσ (z2 ) β(dz2 ).
(38)
Given η > 0 satisfying W1 (α, β) < η, let µ ∈ Γ(α, β) be a coupling such that
ZZ
kz1 − z2 kµ(dz1 × dz2 ) < η .
(39)
We may write
Z
Z
Σα (x, σ) = Qσ (z1 ) µ(dz1 × dz2 ) and Σβ (x, σ) = Qσ (z2 ) µ(dz1 × dz2 ) .
(40)
17

Therefore,
ZZ
kΣα (x, σ) − Σβ (x, σ)k ≤

kQσ (z1 ) − Qσ (z2 )k µ(dz1 × dz2 ) .

(41)

ZZ
σAf
kΣα (x, σ) − Σβ (x, σ)k ≤
kz1 − z2 k µ(dz1 × dz2 )
Cd (σ)
σAf
η.
≤
Cd (σ)

(42)

Lemma 1 and (41) imply that

Since (42) holds for any η > W1 (α, β), we can conclude that
kΣα (x, σ) − Σβ (x, σ)k ≤

σAf
W1 (α, β) ,
Cd (σ)

(43)

as claimed.
d
PnIn what follows, given random vectors yi ∈ R , i ∈ N, we let αn =
i=1 δyi /n.

Corollary 1 (Consistency for Smooth Kernels). Let K be a multiscale kernel
as in Theorem 1 and σ > 0. If α ∈ P1 (Rd ) and yi ∈ Rd , i ∈ N, are i.i.d.
random variables with distribution α, then
n↑∞

sup kΣαn (x, σ) − Σα (x, σ)k −−→ 0
x∈Rd

almost surely.
Proof. Theorem 1 implies that
sup kΣαn (x, σ) − Σα (x, σ)k ≤
x∈Rd

σAf
W1 (αn , α).
Cd (σ)

(44)

The conclusion follows from the fact that W1 metrizes weak convergence
of probability measures in P1 (Rd ) and Varadarajan’s Theorem [35] about
convergence of empirical measures on Polish spaces that ensures that αn
converges weakly to α almost surely.

18

Corollary 1 guarantees the asymptotic consistency of empirical CTFs.
However, in applications, it is important to have estimates of the rate of
convergence, which we derive from the stability theorem and a result of
Fournier and Guillin [5, Theorem 1].
Theorem 2 (Fournier and Guillin, [5]). Let α ∈ Ps (Rd ), where s > 1. If
y1 , . . . , yn are i.i.d. random variables with distribution α and p ∈ [1, s), then
there exists a constant b > 0 that depends only on p, s and d such that
E [Wp (α, αn )] ≤
 s−p
−1
−

if p > d/2 and s 6= 2p;
n s + n 2
s−p
1
p
−
−
b ms (α) · n s + n 2 log(1 + n) if p = d/2 and s 6= 2p;

p
 − s−p
n s + n− d
if p ∈ [1, d/2) and s 6= d/(d − p),
for any n ≥ 1.
Corollary 2. Let f be as in Theorem 1 and σ > 0. Suppose that α ∈ P3 (Rd )
and yi , i ∈ N, are i.i.d. random variables with distribution α. Then, there is
a constant b > 0, that depends only on d, such that


E sup kΣαn (x, σ) − Σα (x, σ)k ≤
x∈Rd
 2
1

n− 3 + n− 2
if d = 1;

σAf b
2
1
m3 (α) · n− 3 + n− 2 log(1 + n) if d = 2;

Cd (σ)
1
 − 23
n + n− d
if d ≥ 3.
Proof. Since P3 (Rd ) ⊆ P1 (Rd ), α is also an element of P1 (Rd ). The conclusion follows by invoking Theorem 1 and the result of Fournier and Guillin
with p = 1 and s = 3. The constant b depends only on d because we are
fixing p and s.
Theorem 1 ensures that multiscale covariance fields are stable. However,
the results do not apply to some discontinuous kernels of practical interest.
Nonetheless, we prove a stability theorem for the truncation kernel, as well
as pointwise convergence results for more general kernels.

19

4.2. The Truncation Kernel
We begin our discussion of covariance fields associated with the truncation
kernel with a stability theorem with respect to the ∞-Wasserstein metric. In
preparation for the proof of the theorem, we introduce some notation. For
0 ≤ a < b, let
Rd (a, b) = {y ∈ Rd : a < kyk ≤ b} .
(45)
and

Z
sd (a, b) =

kyk2 dy

(46)

Rd (a,b)

be its radial moment of inertia. We will use the fact that for any B ≥ b, the
inequality
ωd−1 B d+2
(47)
sd (a, b) ≤ (b − a)
d+2 B−a
holds. Indeed,


ωd−1 d+2
ωd−1 d+2
sd (a, b) =
b
− ad+2 =
a
(b/a)d+2 − 1
d+2
d+2
 
 d+1 !
b
b
ωd−1 d+2 (b − a)
a
1+
=
+ ··· +
d+2
a
a
a
 
 d+1 !
(48)
B
b − a ωd−1 d+2 B − a
B
a
1+
≤
+ ··· +
B−a d+2
a
a
a
≤ (b − a)

ωd−1 B d+2 − ad+2
ωd−1 B d+2
≤ (b − a)
.
(d + 2)
B−a
d+2 B−a

λ
(Rd ) be as in Definition 3.
Let P∞ (Rd ) and P∞

Theorem 3 (Stability for the Truncation Kernel). Let σ > 0 and λ > 0.
Suppose Ω ⊂ Rd is a compact set and let c > 0 satisfy diam(Ω) ≤ c. There is
λ
a constant A = A(σ, d, c) > 0 such that if α ∈ P∞
(Rd ) and β ∈ P∞ (Rd ) have
their supports contained in Ω, then the multiscale covariance tensor fields of
α and β associated with the truncation kernel T satisfy
sup kΣα (x, σ) − Σβ (x, σ)k ≤ λAW∞ (α, β) .
x∈Rd

Proof. Without loss of generality, we may assume that x ∈ Rd is the origin.
We abbreviate η = W∞ (α, β) and let µ ∈ Γ(β, α) be a coupling that realizes
η. Clearly, η ≤ diam(Ω) ≤ c.
20

Using the notation introduced in (45), let Aσ,η = Rd (σ, σ + η). We write
(σ + η)d
Σα (x, σ) − Σβ (x, σ) = Σα (x, σ) −
Σα (x, σ + η)
d
σ{z
|
}
T1

d

(σ + η)
Σα (x, σ + η) − Σβ (x, σ) .
+
d
| σ
{z
}

(49)

T2

λ
(Rd ), we bound the norm of T1 as follows:
Using the fact that α ∈ P∞
Z
Z
1
1
kT1 k = d
y ⊗ y dα(y) ≤ d
kyk2 dα(y)
σ νd Aσ,η
σ νd Aσ,η
Z
λ
λ
≤ d
kyk2 dy = d sd (σ, σ + η) .
σ νd Aσ,η
σ νd

(50)

Using (47) with a = σ, b = σ + η and B = σ + c, we have that
sd (σ, σ + η) ≤ η

ωd−1 (σ + c)d+2
d (σ + c)d+2
= η νd
.
d+2
c
d+2
c

(51)

Thus,
kT1 k ≤

λd (σ + c)d+2
W∞ (α, β) .
d+2
cσ d

(52)

Now we examine T2 . Let I : Rd → R the characteristic function of the closed
ball of radius 1 centered at the origin. Using the coupling µ, we write



ZZ
 y 
1
y1
2
T2 = d
(y1 ⊗ y1 ) I
− (y2 ⊗ y2 ) I
µ(dy1 × dy2 ) .
σ νd
σ
+
η
σ
d
d
R ×R
(53)
Note that the integrand in (53) vanishes on (Rd \ Bσ+η ) × (Rd \ Bσ ) and
the integral also vanishes over Rd \ Bσ+η × Bσ because this subdomain is
disjoint from supp [µ]. Combining these remarks with
y1 ⊗ y1 − y2 ⊗ y2 = (y1 − y2 ) ⊗ y1 + y2 ⊗ (y1 − y2 ) ,

21

(54)

we may rewrite (53) as
ZZ
1
T2 = d
[(y1 ⊗ y1 ) − (y2 ⊗ y2 )] µ(dy1 × dy2 )
σ νd
Bσ+η ×Bσ
ZZ
1
+ d
y1 ⊗ y1 µ(dy1 × dy2 )
σ νd
Bσ+η ×(Rd \Bσ )
ZZ
1
= d
(y1 − y2 ) ⊗ y1 µ(dy1 × dy2 )
σ νd
Bσ+η ×Bσ
ZZ
1
y2 ⊗ (y1 − y2 ) µ(dy1 × dy2 )
+ d
σ νd
Bσ+η ×Bσ
ZZ
1
+ d
y1 ⊗ y1 µ(dy1 × dy2 ) .
σ νd
Bσ+η ×Aσ,2η

(55)

For the last equality, we used again the fact that
(y1 , y2 ) ∈ supp [µ] =⇒ ky1 − y2 k ≤ η .

(56)

From (55) and (56), using the facts that ky1 k ≤ σ + η for y1 ∈ Bσ+η and
ky2 k ≤ σ for y2 ∈ Bσ , we may conclude that
ZZ
σ+η
ky1 − y2 k µ(dy1 × dy2 )
kT2 k ≤ d
σ νd
Bσ+η ×Bσ
ZZ
σ
ky1 − y2 k µ(dy1 × dy2 )
+ d
σ νd
Bσ+η ×Bσ
ZZ
1
+ d
ky1 k2 µ(dy1 × dy2 )
(57)
σ νd
Bσ+η ×Aσ,2η
Z
Z
σ+c
σ
≤ d η
α(dy1 ) + d η
α(dy1 )
σ νd
σ νd Bσ+η
Bσ+η
Z
1
+ d
ky1 k2 α(dy1 ) ,
σ νd A(σ−η)+ ,2η
λ
where (σ − η)+ = max{σ − η, 0}. Since α ∈ P∞
(Rd ) and η ≤ c,
Z
Z
α(dy1 ) ≤ λ
dy1 = λ(σ + η)d νd ≤ λ(σ + c)d νd
Bσ+η

Bσ+η

22

(58)

and
Z

Z

2

ky1 k α(dy1 ) ≤ λ
Aσ−η,2η

ky1 k2 dy1

A(σ−η)+ ,2η

= λsd ((σ − η)+ , σ + η)
ωd−1 (σ + c)d+2
d+2 η+c
d
(σ + c)d+2
νd
,
≤ 2λη
d+2
c

(59)

≤ 2λη

where we used (47) with a = (σ − η)+ , b = σ + η and B = σ + c. Combining
(57), (58) and (59), we obtain
kT2 k ≤

2σ + c
2λ d (σ + c)d+2
d
λ(σ
+
c)
η
+
η
σd
σd d + 2
c

(60)

From (49), (52) and (60), it follows that
d (σ + c)d+2
W∞ (α, β)
d+2
cσ d
(2σ + c)(σ + c)d
+λ
W∞ (α, β)
σd
2d (σ + c)d+2
+λ
W∞ (α, β)
(d + 2)
cσ d
= λA(σ, d, c)W∞ (α, β) .

kΣα (x, σ) − Σβ (x, σ)k ≤ λ

(61)

Since x ∈ Rd is arbitrary, the claim follows.
We now derive a consistency result and estimates for the rate of convergence of empirical approximations to multiscale covariance fields. The
following result is a W∞ -counterpart to the theorem by Fournier and Guillin
stated above.
Theorem 4 (Garcı́a-Trillos and Slepčev [6]). Let Ω ⊂ Rd be a bounded
connected open subset with Lipschitz boundary. Let α be a probability measure
on Ω with density fα : Ω → (0, ∞) such that there exists λ ≥ 1 with λ−1 ≤
fα (x) ≤ λ, for all x ∈ Ω, and let y, . . . , yn be i.i.d. random variables with
distribution α. Then, there exist constants c1 , C1 , C2 > 0, depending only on
Ω and λ, such that for all n ∈ N and p > 1,

√
P W∞ (α, αn ) ≤ (C1 + C2 p) rd (n) ≥ 1 − c1 n−p ,
23

where r2 (n) =

ln(n)3/4
n1/2

and rd (n) =

ln(n)1/d
,
n1/d

for d ≥ 3.

Corollary 3 (Consistency for the Truncation Kernel). Let α be a probabilty
measure on Rd with density fα and let Ωα be the interior of the support of
α. Assume that Ωα is bounded and connected with Lipschitz boundary ∂Ωα .
Furthermore, assume that there exists λ ≥ 1 such that λ−1 ≤ fα (z) ≤ λ, for
all z ∈ Ωα . If yi , i ∈ N, are i.i.d. random variables with distribution α, then,
for any p > 1, there are constants C = C(Ωα , λ, p) > 0 and c1 = c1 (Ωα , λ) >
0 such that


P sup kΣαn (x, σ) − Σα (x, σ)k ≤ C rd (n) ≥ 1 − c1 n−p .
x∈Rd

Here, αn =

Pn

i=1 δyi /n.

√
Proof. We use Theorem 4 and write C 0 = C1 + p C2 . Theorem 3 implies
that there is a constant C 00 = C 00 (Ωα , λ) > 0 such that
sup kΣαn (x, σ) − Σα (x, σ)k ≤ C 00 W∞ (α, αn ) .

(62)

x∈Rd

Thus,



P sup kΣαn (x, σ) − Σα (x, σ)k ≤ C C rd (n) ≥
x∈Rd

≥ P W∞ (α, αn ) ≤ C 0 rd (n) ≥ 1 − c1 n−p .
0

00

(63)

The claim follows by setting C = C 0 C 00 .
Corollary 4. Let σ > 0 and p > 1. Under the assumptions of Corollary 3,
for the truncation kernel, there exist N = N (σ, Ωα , λ) ∈ N and a constant
A = A(σ, Ωα , λ) > 0 such that


E sup kΣαn (x, σ) − Σα (x, σ)k ≤ Ard (n),
x∈Rd

for all n ≥ N.
R∞
Proof. We apply to W∞ (α, αn ) the identity E(Z) = 0 P(Z > t) dt that is
valid for any non-negative random variable Z with finite first moment. Since
W∞ (α, αn ) ≤ D = diam(Ωα ), we get
Z D
E [W∞ (α, αn )] =
P (W∞ (α, αn ) > t) dt .
(64)
0

24

Theorem 4 implies that
P (W∞ (α, αn ) > (C1 +
Let t0 = min{D, (C1 +
Z

√

√
p C2 ) rd (n)) ≤ n−p .

p C2 ) rd (n)}. From (64) and (65),

t0

Z

D

P (W∞ (α, αn ) > t) dt

P (W∞ (α, αn ) > t) dt +

E [W∞ (α, αn )] =

(65)

t0

0

≤ t0 + D P (W∞ (α, αn ) > (C1 +
√
≤ (C1 + p C2 ) rd (n) + Dn−p .

√

p C2 ) rd (n))

(66)

Fixing p, say p = 2, for n sufficiently large, the dominant term on this last
expression is the one involving rd (n). Thus, the claim follows from (66) and
Theorem 3 applied to α and β = αn .
Remark 5. We carry out an experiment to test the convergence rates obtained in Corollary 4. We consider the probability measure α supported
on the unit circle S1 ⊂ R2 induced by the normalized arc length element
(2π)−1 ds. In this case, for the truncation kernel, Σα was calculated explicitly
in Example 4. We consider sets of i.i.d. samples of size n, 10 ≤ n ≤ 106 .
For each n, thirty sets of samples are taken. For each such set, we compute
Σαn and estimate the “error” as max kΣαn (x, σ) − Σα (x, σ)k, for σ = 0.6,
where the maximum is taken over gridpoints on a 24 × 24 grid on the square
[−1.5, 1.5] × [−1.5, 1.5]. We let εn be the average error over all thirty sets of
samples. Figure 6 shows a plot (in blue) of εn in log-log scale. To compare
εn with the predicted rates, we use a least-squares fit, in log-log scale, of
3/4
, also shown in Figure 6 (in red). The disthe form ε = Cr2 (n) = C ln(n)
n1/2
crepancy between the predicted and observed rates suggests that Corollary
4 might not be optimal. A curve of the form ε = Cn−1/2 , shown in green,
produces a tighter fit to the data, suggesting that the optimal bound might
be O(n−1/2 ).
4.3. General Kernels
We conclude the discussion of convergence of empirical CTFs with a
pointwise central limit theorem (CLT) that holds for kernels in the full generality of Definition 2. One may think of it as a CLT for each entry of the
matrix Σα (x, σ). If e1 , . . . , ed is an orthonormal basis of Rd , the (i, j)-entry of
the covariance matrix in this coordinate system is given by Σα (x, σ)(ei , ej ),
25

Figure 6: Log-log plots of experimental error rates (in blue) for empirical covariance fields,
rates predicted by Corollary 4 (in red), and a least-squares fit of order n−1/2 (in green).

the bilinear form Σα (x, σ) evaluated at (ei , ej ). In matrix notation, this is
the same as hei , Σα (x, σ)ej i. More generally, for fixed u, v, x ∈ Rd and σ > 0,
we consider
Z
Σα (x, σ)(u, v) = (y − x) ⊗ (y − x)(u, v)K(x, y, σ) α(dy).
(67)
Consider the random variable
zuv (y) = (y − x) ⊗ (y − x) (u, v)K(x, y, σ),

(68)

where y has distribution α. Clearly,
E [zuv ] = Σα (x, σ)(u, v) .

(69)

Theorem 5 (Central Limit). If f is as in Definition 2, then zuv has finite
2
variance σuv
. Moreover, if zi , i ∈ N, are i.i.d. random variables with the
same distribution as zuv , then
!
n
√
1X
d
2
n
zi − Σα (x, σ)(u, v) →
− N (0, σuv
),
n i=1
2
as n → ∞, where convergence is in distribution and N (0, σuv
) is normally
2
distributed with mean zero and variance σuv .

26

Proof. We show that zuv has finite second moment. From (68) and (4),
Z
Z
2
zuv α(dy) 6 kukkvk ky − xk4 K 2 (x, y, σ) α(dy)


Z
ky − xk4 2 ky − xk2
σ 4 kukkvk
α(dy)
≤
f
(70)
Cd2 (σ)
σ4
σ2
σ 4 kukkvk 2
≤
C .
Cd2 (σ)
The last inequality follows from condition (c) in Definition 2 that ensures
that r2 f 2 (r) < C 2 , for any r > 0. The theorem now follows from a direct
application of the classical CLT.
Remark 6. Note that (70) implies that if kuk = kvk = 1, then
Z
C 2σ4
2
2
,
σuv = zuv
α(dy) − (E [zuv ])2 ≤ 2
Cd (σ)

(71)

giving a uniform bound on the variance of zuv over x ∈ Rd and u, v ∈ Sd−1 .
5. Multiscale Fréchet Functions
The mean of a random vector y ∈ Rd is a simple and yet oftentimes informative, “one-element” summary of the distribution of y. If y has finite second
moment and is distributed according to the probability measure α, then the
mean may be characterized more geometrically as the unique minimizer of
the Fréchet function
Z


2
Fα (x) = E ky − xk = ky − xk2 α(dy) ,
(72)
which measures the spread of y about x ∈ Rd . The mean, however, is not as
effective for complex distributions of practical interest such as multimodal
distributions or those supported in nonlinear subspaces. In this section,
we introduce a multiscale analogue of the Fréchet function that is rich in
information about the shape of the distribution of y. At each fixed scale,
the local minima of the function may be viewed as localized analogues of the
mean, as illustrated in examples below. However, instead of just focusing on
the local extrema, we take the view that it is more informative to investigate
the behavior of the full multiscale Fréchet function, as this lets us uncover
more information about the distribution of y.
27

Definition 5. Let f : [0, ∞) → R be as in Definition 2 with associated kernel
K : Rd ×Rd ×(0, ∞) → R. The multiscale Fréchet function Vα : Rd ×(0, ∞) →
R is defined as
Z
Vα (x, σ) := ky − xk2 K(x, y, σ) α(dy) .
Proposition 3. For each σ > 0, the multiscale Fréchet function satisfies
Vα (x, σ) = tr Σα (x, σ) .
Proof. Let {e1 , . . . , ed } ⊂ Rd be an orthonormal basis. Then,
ky − xk2 =

d
X

hy − x, ei i2 =

d
X

i=1

(y − x) ⊗ (y − x)(ei , ei ) .

(73)

i=1

Hence,
Vα (x, σ) =
=

=

d Z
X

(y − x) ⊗ (y − x)(ei , ei )K(x, y, σ) α(dy)

i=1
d Z
X
i=1
d
X


(y − x) ⊗ (y − x)K(x, y, σ) α(dy) (ei , ei )

(74)

Σα (x, σ)(ei , ei ) = tr Σα (x, σ) ,

i=1

as claimed.
Corollary 5 (Stability). Let f : [0, ∞) → R be as in Definition 2 with multiscale kernel K. Suppose that f is differentiable and there exists a constant
A > 0 such that r3/2 |f 0 (r)| ≤ A, ∀r ≥ 0. Then, there is a constant Af > 0,
that depends only on f , such that
sup |Vα (x, σ) − Vβ (x, σ)| ≤
x∈Rd

σdAf
W1 (α, β), ,
Cd (σ)

for any α, β ∈ P1 (Rd ) and any σ > 0.
Proof. The result follows from Proposition 3, Theorem 1 and the fact that
for any d × d matrix X, |tr X| ≤ d kXk, where kXk is the Frobenius norm
of X.
28

Similarly, Corollary 1 and Proposition 3 yield the following consistency
result for multiscale Fréchet functions.
Corollary 6 (Consistency). Suppose that α ∈ P1 (Rd ). Let yi ∈ Rd , i ∈ N,
be i.i.d. random variables with distribution α and K a multiscale kernel as
in Theorem 1. Then, for each fixed σ > 0,
n↑∞

sup |Vα (x, σ) − Vαn (x, σ)| −−→ 0
x∈Rd

almost surely.
The following result about convergence of multiscale Fréchet functions
are immediate consequences of Corollary 2 and Corollary 5.
Corollary 7. Let f be as in Theorem 1 and σ > 0. Suppose that α ∈ P3 (Rd )
and yi , i ∈ N, are i.i.d. random variables with distribution α. Then, there is
a constant β > 0, that depends only on d, such that


E sup |Vα (x, σ) − Vαn (x, σ)| ≤
x∈Rd
 2
1

n− 3 + n− 2
if d = 1;

σdAf β
2
1
m3 (α) · n− 3 + n− 2 log(1 + n) if d = 2;

Cd (σ)
1
 − 32
n + n− d
if d ≥ 3.
Remark 7. Analogous stability and consistency results for the truncation
kernel follow from Theorem 3, Corollary 3 and Corollary 4.
For more general kernels, the following pointwise central limit theorem
holds. For fixed x ∈ Rd and σ > 0, let
t(y) = ky − xk2 K(x, y, σ),

(75)

whose expected value is E [t] = Vα (x, σ). As in Theorem 5, the variance of t
is finite and denoted σt2 .
Theorem 6 (Central Limit). Let f be as in Definition 2. If ti ∈ R, i ∈ N,
are i.i.d. random variables with the same distribution as t, then
!
n
√
1X
d
n
ti − Vα (x, σ) →
− N (0, σt2 ) ,
n i=1
as n → ∞, where convergence is in distribution and N (0, σV2 ) is normally
distributed with mean zero and variance σt2 .
29

Multiscale Fréchet functions not only give stable representations of probability measures, but any probability measure α may be fully recovered from
its multiscale Fréchet function associated with the Gaussian kernel, as the
following result shows.
Proposition 4. Let σ > 0 be fixed. Any probability measure α is completely
determined by the Fréchet function Vα (·, σ) associated with the Gaussian kernel at scale σ.
Proof. Let hσ : Rd → R be given by


kxk2
kxk2
exp − 2 .
hσ (x) =
(2πσ 2 )d/2
2σ

(76)

Then, for the Gaussian kernel, we may express the multiscale Fréchet function
as the convolution Vα (x, σ) = (hσ ∗ α)(x). Under Fourier transform, for each
fixed σ > 0, we obtain
Vbα (ξ, σ) = b
hσ (ξ) φα (−2πξ) ,
where φα is the characteristic function of α defined as φα (ξ) =
Therefore,
φα (ξ) = Vbα (−ξ/2π, σ)/ b
hσ (−ξ/2π)

(77)
R

eihx,ξi α(dx).
(78)

provided that b
hσ (−ξ/2π) 6= 0. A calculation shows that

 2


σ kξk2
σ 2 kξk2
2
b
exp −
,
(79)
hσ (−ξ/2π) = σ d −
π
2π
√
which only vanishes at points ξ on the sphere of radius ρσ = πd/σ about
b
the origin.
√ Thus, (78) implies that we can recover φα (ξ) from Vα (·, σ), if
kξk 6= πd/σ. By continuity, we can recover φα (ξ), for any ξ. The claim
now follows from the fact that the characteristic function φα determines α
[35].
The following examples illustrate how information about the shape of
data can be extracted from multiscale Fréchet functions.
Example 7. We consider n = 400 data points distributed into two clusters
of 200 points, each sampled from a Gaussian of variance 0.36 centered at
different points. The data points are plotted in blue in Fig. 7(a), which
30

also shows the empirical Fréchet function Vn at scale σ = 3. The local
minima of Vn captures what is perceived as the “centers” of the two clusters
at that scale. However, more information about the data distribution can
be uncovered from Vn . For example, the local minima may be viewed as
attractors of the (negative) gradient field −∇Vn , indicated by the arrows in
the figure. The stable manifold of each attractor, which comprises points
that move toward the attractors under the associated flow may be viewed
as clusters inferred from the data at that scale. These clusters are delimited
by the repellers of the system, which correspond to the local maxima of Vn .
Fig. 7(b) shows how Vn varies across scales, highlighting the bifurcation of
the attractors (in red) and repellers (in green) as σ changes. In data analysis,

(a)

(b)

Figure 7: (a) Fréchet function for data on the line (highlighted in blue) computed with
the Gaussian kernel at scale σ = 3; (b) Fréchet function across scales.

such bifurcation diagrams may find several applications. For example, if the
data represent the distribution of some phenotypic trait for two species that
have evolved from a single group, the multiscale Fréchet function and the
associated bifurcation diagram let us create an evolutionary model for the
trait from the observed data.
Example 8. Here we consider the dataset in R2 shown in panel (a) of Fig. 8.
Panels (b)–(h) show the Fréchet function for the Gaussian kernel calculated
at increasing scales. The gradient field −∇Vn at scale σ = 2.25 is depicted
in panel (a) of Fig. 9 along with the two attractors p1 and p2 , and their
stable manifolds that were estimated numerically.
The stable manifolds
may be viewed as estimations at scale σ = 2.25 of clusters of the underlying
probability measure α from which the data was sampled. Panel (b) shows
the gradient field and the covariance tensors at the attractors depicted as
31

data

σ =1.0

σ = 1.5

σ = 2.0

σ = 3.0

σ = 4.5

σ = 5.0

σ = 6.0

σ = 7.0

σ = 9.0

Figure 8: Heat maps of the multiscale Fréchet function for 2D data at increasing scales
computed with the Gaussian kernel.

(a)

(b)

Figure 9: (a) 2D data, attractors and their stable manifolds at a fixed scale (σ = 2.25);
(b) gradient vector field and covariance tensors at the attractors.

ellipses with principal radii proportional to the square root of the eigenvalues
of the covariance matrix. This may be viewed as a localized analogue of
principal component analysis (PCA) that is able to uncover geometry that
is not detectable with standard PCA. Analysis of the spectra of Σn (pi , σ),
i = 1, 2, suggests that the data is organized around two one-dimensional
clusters, whereas standard PCA is not sensitive to the local dimensionality
because of the orientation of the clusters.
32

These examples are intended as proof-of-concept illustrations. Topological and other methods will be explored in forthcoming work for extraction
of structural information from Vn .
6. Hierarchical Manifold Clustering
Clustering is a central theme in pattern analysis with a rich history; cf.
[36]. One of the most studied forms of the problem is that of partitioning a
dataset into various subsets if there is some form of spatial separation of the
data into subgroups. Motivated by problems in such areas as computer vision
and video analysis, cf. [37], there has been a growing interest in clustering
data that are organized as a finite union of possibly intersecting subspaces
that have some special geometric structure [38, 16, 18, 26]. As illustrated in
Fig. 1, the data may consist of noisy samples from an arrangement of (affine)
linear subspaces of a Euclidean space such as a collection of lines in a plane,
or an arrangement of lines and planes in R3 . More generally, the clusters may
comprise a finite collection of possibly non-linear, smooth submanifolds of a
Euclidean space that intersect transversely. Here we propose an approach
to manifold clustering based on CTFs. The basic idea is to use covariance
fields to incorporate directional information at each data point. Formally,
this is achieved via a section of the tensor bundle Rd × (Rd ⊗ Rd ) over Rd , as
follows. Given a probability measure α and a multiscale kernel, let Σα (x, σ)
be the associated CTF. For each σ > 0, consider the section ια;σ : Rd →
Rd × (Rd ⊗ Rd ) given by x 7→ (x, Σα (x, σ)). On the total space of the tensor
bundle, define the metric
1/2
k(x, Σ) − (x0 , Σ0 )kγ = kΣ − Σ0 k2 + γ 2 kx − x0 k2
,
(80)
where x, x0 ∈ Rd , Σ, Σ0 ∈ Rd ⊗ Rd , and γ ≥ 0 is a parameter that balances
the contributions of the spatial and tensor components. Note that k · k0 only
defines a pseudo-metric since k · k0 disregards “horizontal” distances.
For any subset X ⊆ Rd , we denote by Xα;γ,σ the metric space (X, dα;γ,σ ),
where
dα;γ,σ (x, x0 ) = kια;σ (x) − ια;σ (x0 )kγ .
(81)
For a dataset A = {a1 , . . . , an } ⊂ Rd , the proposed clustering method is
based on the single-linkage method [39] applied to the finite
Pn metric space
−1
Aαn ;γ,σ associated with the empirical measure αn = n
i=1 δai . Equivalently, clustering is based on the n × n affinity matrix D whose (i, j)-entry
33

is
dij = dαn ;γ,σ (ai , aj ) .

(82)

Recall that single linkage on a finite metric space A = (A, dA ) starts from n
clusters, each a singleton {ai }, 1 6 i 6 n, sequentially merging the closest
clusters until all data points coalesce into a single cluster. Closeness of two
clusters, say A1 , A2 ⊂ A, is measured by the inter-cluster distance
dsl (A1 , A2 ) =

min

a∈A1 ,a0 ∈A2

dA (a, a0 ) ,

(83)

We choose single linkage because it yields stable dendrograms, as expounded
below, under assumptions on the probability measure from which the data
is sampled that are not very restrictive. Combined with our stability and
consistency results for covariance fields, this guarantees that the manifold
clustering method is stable at all stages.
6.1. Dendrogram Stability
We denote a metric space by X = (X, dX ). An ultrametric space is a pair
(X, uX ), where uX : X × X → R+ is a metric on X that satisfies the strong
triangle inequality
uX (x, x0 ) ≤ max {uX (x, x00 ), uX (x00 , x0 )} ,

(84)

for all x, x0 , x00 ∈ X. Any such function uX is called an ultrametric on X.
As proved in [32], dendrograms over a finite set X are in structurepreserving, bijective correspondence with ultrametrics on X. In this formulation, a hierarchical clustering method can be regarded as a map H : M → U
from finite metric spaces into finite ultrametric spaces. Henceforward, H will
denote the map given by single linkage hierarchical clustering. It is known
[32] that if X = (X, dX ) ∈ M, then H(X) = (X, uX ) is given by
uX (x, x0 ) =

min

x=x0 ,...,xr =x0

max dX (xi , xi+1 ) .
i

(85)

The minimum above is taken over all finite ordered sequences x0 , x1 , . . . , xr
of points in X such that x0 = x and xr = x0 . If x, x0 ∈ X, then uX (x, x0 )
may be interpreted as the dendrogram level at which the clusters containing
x and x0 first merge. This is known as the cophenetic distance between x
and x0 .

34

The main goal of this section is to formulate and prove stability of the
map M 3 X 7→ H(X) ∈ U. The question of stability of single linkage
clustering can be approached using ideas related to the Gromov-Hausdorff
distance [40], as follows. A correspondence R between two sets X and Y is
a subset of X × Y such that π1 (R) = X and π2 (R) = Y , where π1 and π2
denote projections onto the first and second factors. Given X and Y , we
denote by R(X, Y ) the set of all correspondences between X and Y .
Definition 6. Let X and Y be compact metric spaces.
(i) The distortion of a correspondence R between X and Y is defined by
dis (R; X, Y) :=

max

(x,y),(x0 ,y 0 )∈R

|dX (x, x) − dY (y, y 0 )| .

(ii) The Gromov-Hausdorff distance between X and Y is given by
dGH (X, Y) :=

1
inf dis (R; X, Y),
2 R

where the infimum is taken over all correspondences between X and Y.
The following stability result is a generalization of [32, Proposition 26].
Proposition 5. For any X, Y ∈ M and any correspondence R ∈ R(X, Y ),
dis(R; H(X), H(Y)) ≤ dis(R; X, Y) .
As a consequence, dGH (H(X), H(Y)) ≤ dGH (X, Y).
Remark 8. The claim of the proposition may be written, equivalently, as
follows. If uX and uY denote the ultrametrics produced by single linkage
hierarchical clustering on X and Y, then
|uX (x, x0 ) − uY (y, y 0 )| ≤

max
0 0

(x,y),(x ,y )∈R

|dX (x, x0 ) − dY (y, y 0 )| ,

(86)

for any correspondence R between X and Y and all (x, y), (x0 , y 0 ) ∈ R.
Proof of Proposition 5. We prove (86). Given a correspondence R ∈ R(X, Y )
and (x, y), (x0 , y 0 ) ∈ R, let x = x0 , x1 , . . . , xn = x0 in X be such that
maxi dX (xi , xi+1 ) = uX (x, x0 ). Let y0 = y, yn = y 0 and choose y1 , . . . , yn−1 ∈
35

Y such that (xi , yi ) ∈ R for all i = 1, . . . , n − 1. This is possible since any
correspondence R satisfies π1 (R) = X. Notice that
uY (y, y 0 ) ≤ max dY (yi , yi+1 )
i

≤ max (dX (xi , xi+1 ) + |dX (xi , xi+1 ) − dY (yi , yi+1 )|)
i

≤ max dX (xi , xi+1 ) +
i

0

= uX (x, x ) +

max
0 0

(x,y),(x ,y )∈R

max

(x,y),(x0 ,y 0 )∈R

|dX (x, x0 ) − dY (y, y 0 )|

(87)

|dX (x, x0 ) − dY (y, y 0 )|.

The claim follows since (87) also holds if we reverse the roles of X and Y .
Lemma 2. Let α, β ∈ P∞ (Rd ) and σ > 0. If a kernel satisfies the conditions
of Lemma 1, then
sup kΣα (a, σ) − Σβ (b, σ)k ≤
(a,b)∈Rµ

2Af σ
Cd (σ)

sup ka − bk ,
(a,b)∈Rµ

for any coupling µ ∈ Γ(α, β), where Rµ := supp [µ] and Af > 0 is as in
Lemma 1.
Proof. Set ζ = sup(a,b)∈supp [µ] ka − bk. Let µ ∈ Γ(α, β) and (y, y 0 ), (a, b) ∈
supp [µ]. In the notation of Lemma 1, setting z1 = y − a and z2 = y 0 − b, we
have
Af σ
2Af σ
kQσ (y − a) − Qσ (y 0 − b)k ≤
(ky − y 0 k + ka − bk) ≤
ζ , (88)
Cd (σ)
Cd (σ)
where in the last inequality we used ky − y 0 k ≤ ζ and ka − bk ≤ ζ. Since
ZZ
kΣα (a, σ) − Σβ (b, σ)k ≤
kQσ (y − a) − Qσ (y 0 − b)k µ(dy × dy 0 ) , (89)
the lemma follows.
Lemma 3 (Lemma 2.2 of [41]). Let α, β ∈ P∞ (Rd ). Then, for any coupling
µ ∈ Γ(α, β), Rµ = supp [µ] gives a correspondence between A = supp [α] and
B = supp [β].
Theorem 7. Let α, β ∈ P∞ (Rd ), A = supp [α], B = supp [β], σ > 0 and
γ ≥ 0. Then, for any kernel satisfying the conditions of Lemma 1,


2Af σ
dGH ((A, dα;σ,γ ), (B, dβ;σ,γ )) ≤
+ γ W∞ (α, β) ,
Cd (σ)
with Af > 0 as in Lemma 1.
36

Proof. Let µ ∈ Γ(α, β) be a coupling that realizes W∞ (α, β). By Lemma 3,
Rµ = supp [µ] is a correspondence between A and B. Thus,
 1
dGH Aα;σ,γ , Bβ;σ,γ ≤ dis(Rµ ; A, B)
2
1
=
|dα;σ,γ (a, a0 ) − dβ;σ,γ (b, b0 )|
sup
2 (a,b),(a0 ,b0 )∈Rµ

1
kΣα (a, σ) − Σβ (b, σ)k +
≤
sup
2 (a,b),(a0 ,b0 )∈Rµ
0

0

0

0

+ Σα (a , σ) − Σβ (b , σ) + γka − bk + γka − b k



(90)

≤ sup kΣα (a, σ) − Σβ (b, σ)k + γ sup ka − bk
(a,b)∈Rµ


≤

2Af σ
+γ
Cd (σ)

(a,b)∈Rµ


sup ka − bk,
(a,b)∈Rµ

where the last step follows from Lemma 2. The conclusion follows since
W∞ (α, β) = sup(a,b)∈Rµ ka − bk.
Combining Proposition 5 and Theorem 7, we obtain:
Corollary 8 (Stability of Hierarchical Manifold Clustering). Let α, β ∈
P∞ (Rd ) be probability measures with finite support, A = supp [α], B =
supp [β], σ > 0 and γ ≥ 0. Then,


2Af σ
+ γ W∞ (α, β).
dGH (H(Aα;γ,σ ), H(Bα;γ,σ )) ≤
Cd (σ)
6.2. Comments About Consistency of Hierarchical Manifold Clustering
A question that our paper leaves open is whether, under a reasonable
generative model for the sampling from a collection of intersecting manifolds one may be able to prove that the empirical dendrogram converges in
probability to a dendrogram that represents the spatial organization of the
underlying manifolds. In the simpler context of clustering n i.i.d. samples
{x1 , . . . , xn } from a compact metric space (X, dX ) endowed with a Borel
probability measure αX , it has been established in [32] that single linkage hierarchical clustering converges to a dendrogram whose hierarchical structure
depends on the support of αX in a precise way.

37

In the case of flat clustering, i.e. when the goal is to obtain a single partition of the dataset, consistency results for some multi-manifold clustering
methods based on local PCA are given in [16, 17, 18].
6.3. Examples and Applications
Let X = {x1 , . . . , xn } be a dataset in Rd . For γ, σ > 0, we apply the
single linkage method to the metric space Xαn ;γ,σ = (X, dαn ;γ,σ ), where αn is
the empirical measure associated with X and dαn ;γ,σ is the distance defined
in (81). The ultrametric associated with H(Xαn ;γ,σ ) is abbreviated uαn ;γ,σ .
In this setting, analyzing informative dendrogram cutoff levels often is
an important task, which can be approached in different ways, depending on
the nature of the problem. For example, a cutoff level h may be based on a
pre-assigned number of clusters, be learned from training data, or be more
exploratory. We give examples that illustrate all three viewpoints.
Example 9 (Lines and planes). In this experiment we consider the unlabeled
point cloud in Fig. 1(a) that represents an arrangement of two parallel planes
and two lines that intersect the planes transversely. Each plane contains 225
points on a uniform grid and each line contains 30 equally spaced points.
Cutting the dendrogram at four clusters, our method finds the four affine
linear subspaces accurately with the Gaussian kernel at σ = 0.6. In this
case, it is important to choose γ 6= 0 since the spatial component of (82) is
needed to discriminate the parallel planes. In Fig. 1(a), the points are colored
according to cluster membership. In this case, the covariance tensors at data
points on the planes that are away ’from the cluster intersections have two
dominating eigenvalues, whereas for points on the lines they have only one
such eigenvalue. Thus, an analysis of the spectrum of the covariance tensors
at data points let us infer the dimension of each cluster.
Example 10 (A Line Arrangement). In this example, the point-cloud data
represents three intersecting lines in R2 , as shown in Fig. 10(a). Each line
segment is sampled at 200 equally spaced points. Since the slopes of the
lines are different, we expect the covariance matrices to be able to cluster the
points without the aid of additional spatial information. Thus, we set γ = 0
in (82) and σ = 0.4. The number of clusters was set to six to test the ability
of the algorithm to detect not only the lines, but also the three intersections.
Fig. 10(b) shows the single-linkage dendrogram, highlighting each of the six
clusters. The data points are colored according to cluster membership. As
expected, well delineated clusters are detected away from the intersection
38

(a)

(b)

(c)

(d)

Figure 10: (a) an arrangement of three lines and (b) clustering dendrogram; (c) noisy lines
with outliers and (d) clustering dendrogram.

points because the covariance matrices are highly anisotropic with principal
axes that align well with the corresponding line segments. Although the
covariance matrices are not as anisotropic near the intersection points, there
are enough differences in their behavior near the three intersection loci for
the algorithm to be able to place them into different clusters.
The next two examples are of a more exploratory nature in that dendrogram cutoff was chosen through experimentation with the data.
Example 11 (Noisy Lines with Outliers). This is a noisy version of Example
10, as shown in Fig. 10(c). As before, each line is represented by 200 points,
but we have added Gaussian noise of width 0.015 to the data, as well as 180
outliers sampled from the uniform distribution on a rectangle containing the
lines. Because of the nature of the data, the number of clusters was set to
m = 80 so that the three main clusters did not get merged because of the
outliers. The figure also shows a line fitted to each of the three largest clusters
using principal component analysis. The method was able to sharply recover
the three lines, even in the presence of noise and outliers. The majority of
the 80 clusters are singletons of outliers and these are colored black in the
figure. We remark that the choice of σ = 0.51 is crucial when dealing with
data contaminated by noise. In this case, it was also important to set γ 6= 0
to better cope with noise.
Example 12 (Floor cracks). We apply the clustering method to segmentation
of two images of concrete floor cracks. Panel (a) of Fig. 11 shows the original images, whereas panel (b) shows binary images obtained from an edge
detection algorithm. We cluster the foreground pixels of the binary images.
39

As in Example 11, it is important to allow a fairly large number of clusters
so that the clusters that detect the main cracks do not get merged because
of the noisy pixels. Panels (c) and (d) show the outputs (not to scale) of the
clustering algorithm.

(a)

(b)

(c)

(d)

Figure 11: (a) original and (b) processed images of floor cracks; (c) and (d) show clustering
based on CTFs.

To further test the ability of the method to cluster intersecting manifolds,
we experimented with synthetic data comprising multiple arrangements such
as the intersecting lines in Fig. 1(b).
Example 13. We consider three syntehtic datasets of point clouds representing random arrangements of: (i) three line segments in R2 ; (ii) four curves in
R2 that are either line segments or arcs of parabolas; and (iii) three patches
of planes in R3 . Each of these datasets contains a total of 250 point clouds,
50 used for training the algorithm and 200 test samples. The points in each
point cloud are labeled to allow quantification of the accuracy of the output
of the algorithm. Fig. 12 shows a few samples from each of these datasets.
Parameter values that optimize classification performance are learned from
the training samples. Note, however, that even though the number k of
clusters is known, specifying a height h that yields precisely k clusters may
yield undesirable results. For example, important clusters representing different components of an arrangement of manifolds may get merged due to
the presence of outliers or the behavior near the intersections. Thus, it is
often preferable to choose a lower cutoff level before this phenomenon occurs
at the expense of getting a larger number ` of clusters. In such situations,
we select the largest k clusters and assign each point in the remaining (` − k)
clusters to the closest of the top k clusters. Experiments indicate that a good
40

Figure 12: Random arrangements of line segments (row 1), segments of lines and parabolas
(row 2), and plane patches (row 3).

baseline for the cutoff level h is the mean cophenetic distance, which for a
point cloud {x1 , . . . , xn } ⊂ Rd is given by
h0 =

X
2
uα ;γ,σ (xi , xj ) .
n(n − 1) i<j n

(91)

In the learning phase, we typically search for h in a neighborhood of h0 whose
width is determined by the variance of the distribution of the cophenetic
distances.
With the learned parameter values, the algorithm performs well in all
three cases. For each point cloud we count the number of misclassified points
and calculate the average error (AE) and the mean error (ME) rates over all
test samples obtaining:
(i) arrangements of lines: 9.59% (AE) and 4.17% (ME);
(ii) lines and parabolas: 9.93% (AE) and 3.38% (ME);
(iii) arrangements of planes: 7.00% (AE) and 2.42% (ME).
As expected, a closer inspection of the results reveals that most of the errors
occur at points near the intersections of the clusters, where the covariance
tensors are not as informative for clustering purposes.
41

7. Concluding Remarks
We introduced the notion of multiscale covariance tensor fields associated with Euclidean random variables and developed a framework for the
systematic study of the shape of data using localized covariance tensors. We
investigated foundational questions such as stability and consistency of multiscale CTFs, provided illustrations of how CTFs let us uncover geometry
underlying data, and applied the methods to manifold clustering. We also
introduced multiscale Fréchet functions, which are scalar fields derived from
CTFs that fully capture the distributions of random vectors. Multiscale
Fréchet functions are particularly well suited for extension of the methods
of this paper to non-Euclidean random variables, a problem that is receiving
ever increasing attention in data science. In this setting, the goal is to devise
methods that can cope with random variables taking values in spaces such as
Riemannian manifolds and more general metric spaces. Unless restrictive assumptions are imposed on the sample space and the distributions, CTFs may
be difficult to define in this nonlinear realm. In contrast, the Fréchet function
formulation can be easily extended to metric spaces supporting a diffusion
kernel [42]. In forthcoming work, we will investigate theoretical and computational aspects of such extensions, including the accessibility of information
residing in multiscale Fréchet functions, a problem that poses computational
challenges even in the case of high-dimensional Euclidean random variables.
In this paper, we only considered radial basis kernels; however, many
results extend easily to more general kernels. We emphasized the multiscale formulation largely because of the questions that motivated this work.
Nevertheless, the majority of the results apply to kernels that are not scale
dependent.
Covariance tensor fields also suggest ways of formalizing the notion of
shape of Euclidean data and probability measures. For example, for a distribution α with the property that the covariance tensor field Σα (·, σ) associated
with a smooth kernel (such as the Gaussian kernel) is non-singular for every x ∈ Rd , Σ−1
α (·, σ) defines a metric tensor with close ties to α. This poses
the problem of uncovering relationships between Riemannian metrics derived
from CTFs, such as Σ−1
α (·, σ), and the shape of α.
In a different direction, for a fixed point x ∈ Rd , an interesting problem
is that of capturing the values of σ for which Σα (x, σ) exhibits a “jump”
in behavior. This study, in the context of images, gives rise to notions of
local scales. Knowledge of local scales for each point x leads to criteria
42

for selecting important, salient points in the spirit of SIFT [43, 2]. The
concept of local scales arose first in the context of images [44] and was later
extended to probabilty distributions [45]. The notions of local scales in [44,
45] were isotropic. Thus, future developments related to characterizing shape
using CTFs are suggested by the possibility of constructing notions of local
scales on general shapes [45, 2] which —by exploiting the tensorial nature of
Σα (x, σ)— become sensitive to direction.
Data Accessibility
The synthetic data used in the manifold clustering experiments is available at https://bitbucket.org/diegodiaz-math/ctf-files/.
Acknowledgements
This research was supported in part by NSF grants IIS-1422400, DMS1418007 and DBI-1262351, and by the Erwin Schrödinger Institute in Vienna.
We thank Dejan Slepčev for a discussion about the results of [6].
References
[1] A. Little, Y.-M. Jung, M. Maggioni, Estimation of intrinsic dimensionality of samples from noisy low-dimensional manifolds in high dimensions
with multiscale SVD, in: IEEE/SP 15th Workshop on Statistical Signal
Processing, 2009, pp. 85–88.
[2] D. Dı́az Martı́nez, F. Mémoli, W. Mio, Multiscale covariance fields, local
scales, and shape transforms, in: F. Nielsen, F. Barbaresco (Eds.), Geometric Science of Information, Vol. 8085 of Lecture Notes in Computer
Science, Springer, 2013, pp. 794–801.
[3] C. Villani, Topics in Optimal Transportation, Vol. 58 of Graduate Studies in Mathematics, American Mathematical Society, Providence, RI,
2003.
[4] C. Villani, Optimal Transport, Old and New, Springer, 2009.
[5] N. Fournier, A. Guillin, On the rate of convergence in Wasserstein distance of the empirical measure, Probability Theory and Related Fields.
43

[6] N. Garcı́a-Trillos, D. Slepčev, On the rate of convergence of empirical
measures in ∞-transportation distance, Canadian Journal of Mathematics, doi:10.4153/CJM-2014-044-6.
[7] W. Allard, G. Chen, M. Maggioni, Multi-scale geometric methods for
data sets II: Geometric multi-resolution analysis, Appl Comput Harmonic Anal 32 (2012) 435–462.
[8] J. Berkmann, T. Caelli, Computation of surface geometry and segmentation using covariance techniques, IEEE Trans Pattern Anal Mach Intell
16 (1994) 1114–1116.
[9] R. Rusu, Semantic 3D object maps for everyday manipulation in human living environments, Ph.D. thesis, Technische Universität München
(2009).
[10] G. Medioni, M.-S. Lee, C.-H. Tang, A Computational Framework for
Segmentation and Grouping, Elsevier Science, 2000.
[11] T. Brox, B. Rosenhahn, D. Cremers, H.-P. Seidel, Nonparametric density estimation with adaptive, anisotropic kernels for human motion
tracking, in: Human Motion – Understanding, Modeling, Capture and
Animation, Vol. 4814 of Lecture Notes in Computer Science, Springer,
2007, pp. 152–165.
[12] D. Kushnir, M. Galun, A. Brandt, Fast multiscale clustering and manifold identification, Pattern Recognition 39 (10) (2006) 1876–1891.
[13] A. B. Goldberg, X. Zhu, A. Singh, Z. Xu, R. D. Nowak, Multi-manifold
semi-supervised learning., in: AISTATS, 2009, pp. 169–176.
[14] D. Gong, X. Zhao, G. Medioni, Robust multiple manifolds structure
learning, arXiv preprint arXiv:1206.4624.
[15] Y. Wang, Y. Jiang, Y. Wu, Z.-H. Zhou, Spectral clustering on multiple
manifolds, IEEE Transactions on Neural Networks 22 (7) (2011) 1149–
1161.
[16] E. Arias-Castro, Clustering based on pairwise distances when the data is
of mixed dimensions, IEEE Transactions on Information Theory 57 (3)
(2011) 1692–1706.
44

[17] E. Arias-Castro, G. Chen, G. Lerman, et al., Spectral clustering based
on local linear approximations, Electronic Journal of Statistics 5 (2011)
1537–1587.
[18] E. Arias-Castro, G. Lerman, T. Zhang, Spectral clustering based on
local pca, arXiv preprint arXiv:1301.2007.
[19] R. Vidal, Y. Ma, S. Sastry, Generalized principal component analysis
(GPCA), IEEE Trans Pattern Anal Mach Intell 27 (12) (2005) 1945–
1959.
[20] V. Govind, A tensor decomposition for geometric grouping and segmentation, in: IEEE Conference on Computer Vision and Pattern Recognition, Vol. 1, 2005, pp. 1150–1157.
[21] G. Chen, G. Lerman, Spectral curvature clustering (SCC), Int J Comput
Vis 81 (3) (2009) 317–330.
[22] G. Lerman, T. Zhang, Robust recovery of multiple subspaces by geometric `p minimization, Ann Stat 39 (5) (2011) 2686–2715.
[23] T. Zhang, A. Szlam, Y. Wang, G. Lerman, Hybrid linear modeling via
local best-fit flats, Int J Comput Vis 100 (3) (2012) 217–240.
[24] G. Liu, Z. Lin, S. Yan, J. Sun, Y. Ma, Robust recovery of subspace
structures by low-rank representation, IEEE Trans Pattern Anal Mach
Intell 35 (1) (2013) 171–184.
[25] E. Elhamifar, R. Vidal, Sparse subspace clustering: algorithm, theory,
and applications, IEEE Trans Pattern Anal Mach Intell 35 (11) (2013)
2765–2781.
[26] M. Soltanolkotabi, E. Elhamifar, E. J. Candes, et al., Robust subspace
clustering, The Annals of Statistics 42 (2) (2014) 669–699.
[27] M. Polito, P. Perona, Grouping and dimensionality reduction by locally
linear embedding, in: NIPS, 2001, pp. 1255–1262.
[28] A. Gionis, A. Hinneburg, S. Papadimitriou, P. Tsaparas, Dimension
induced clustering, in: Proceedings of the eleventh ACM SIGKDD international conference on Knowledge discovery in data mining, ACM,
2005, pp. 51–60.
45

[29] R. Souvenir, R. Pless, Manifold clustering, in: International Conference
on Computer Vision, 2005, pp. 648–653.
[30] G. Haro, G. Randall, G. Sapiro, Stratification learning: Detecting mixed
density and dimensionality in high dimensional point clouds, in: NIPS,
2006, pp. 553–560.
[31] X. Wang, K. Slavakis, G. Lerman, Riemannian multi-manifold modeling,
arXiv preprint arXiv:1410.0095.
[32] G. Carlsson, F. Mémoli, Characterization, stability and convergence of
hierarchical clustering methods, J Mach Learn Res 11 (2010) 1425–1470.
[33] J. Bruce, P. Giblin, Curves and singularities: a geometrical introduction
to singularity theory, Cambridge University Press, 1992.
[34] C. Givens, R. Shortt, A class of Wasserstein metrics for probability
distributions, Michigan Math J 31 (2) (1984) 231–240.
[35] R. M. Dudley, Real Analysis and Probability in Cambridge Studies in
Advanced Mathematics, Vol. 74, Cambridge University Press, Cambridge, England, 2002.
[36] A. Jain, R. Dubes, Algorithms for Clustering Data, Prentice Hall, Englewood Cliffs, NJ, 1988.
[37] R. Tron, R. Vidal, A benchmark for the comparison of 3-D motion
segmentation algorithms, in: IEEE Conference on Computer Vision and
Pattern Recognition, 2007, pp. 1–8.
[38] G. Chen, G. Lerman, Foundations of a multi-way spectral clustering
framework for hybrid linear modeling, Foundations of Computational
Mathematics 9 (5) (2009) 517–558.
[39] R. Sibson, SLINK: an optimally efficient algorithm for the single-link
cluster method, Comp J 16 (1) (1973) 30–34.
[40] D. Burago, Y. Burago, S. Ivanov, A Course in Metric Geometry, Graduate Studies in Mathematics, American Mathematical Society, Providence, RI, 2001.

46

[41] F. Mémoli, Gromov-Wasserstein distances and the metric approach to
object matching, Foundations of Computational Mathematics (2011) 1–
71.
[42] R. Coifman, S. Lafon, Diffusion maps, Appl Comput Harmonic Anal
21 (1) (2006) 5–30.
[43] D. G. Lowe, Distinctive image features from scale-invariant keypoints,
International Journal of Computer Vision 60 (2) (2004) 91–110.
[44] P. W. Jones, T. M. Le, Local scales and multiscale image decompositions, Applied and Computational Harmonic Analysis 26 (3) (2009)
371–394.
[45] T. Le, F. Mémoli, Local scales on curves and surfaces, Applied and
Computational Harmonic Analysis 33 (3) (2012) 401 – 437.

47

