Adaptive Geometric Multiscale Approximations for
Intrinsically Low-dimensional Data
Wenjing Liao1

wliao@math.jhu.edu

Mauro Maggioni1,2,3

mauro.maggioni@jhu.edu
Departments of Mathematics, 2 Applied Mathematics and Statistics,
3
The Institute for Data Intensive Engineering and Science
Johns Hopkins University, 3400 N. Charles Street, Baltimore, MD 21218, USA

arXiv:1611.01179v2 [stat.ML] 18 Jul 2017

1

Editor:

Abstract
We consider the problem of efficiently approximating and encoding high-dimensional data
sampled from a probability distribution ρ in RD , that is nearly supported on a d-dimensional
set M - for example supported on a d-dimensional Riemannian manifold. Geometric MultiResolution Analysis (GMRA) provides a robust and computationally efficient procedure to
construct low-dimensional geometric approximations of M at varying resolutions. We introduce a thresholding algorithm on the geometric wavelet coefficients, leading to what
we call adaptive GMRA approximations. We show that these data-driven, empirical approximations perform well, when the threshold is chosen as a suitable universal function
of the number of samples n, on a wide variety of measures ρ, that are allowed to exhibit
different regularity at different scales and locations, thereby efficiently encoding data from
more complex measures than those supported on manifolds. These approximations yield
a data-driven dictionary, together with a fast transform mapping data to coefficients, and
an inverse of such a map. The algorithms for both the dictionary construction and the
transforms have complexity Cn log n with the constant linear in D and exponential in d.
Our work therefore establishes adaptive GMRA as a fast dictionary learning algorithm
with approximation guarantees. We include several numerical experiments on both synthetic and real data, confirming our theoretical results and demonstrating the effectiveness
of adaptive GMRA.

Keywords: Dictionary Learning, Multi-Resolution Analysis, Adaptive Approximation,
Manifold Learning, Compression

1. Introduction
We model a data set as n i.i.d. samples Xn := {xi }ni=1 from a probability measure ρ in
RD . We make the assumption that ρ is supported on or near a set M of dimension d  D,
and consider the problem, given Xn , of learning a data-dependent dictionary that efficiently
encodes data sampled from ρ.
In order to circumvent the curse of dimensionality, a popular model for data is sparsity: we say that the data is k-sparse on a suitable dictionary (i.e. a collection of vectors)
c Adaptive Geometric Multiscale Approximations.

Liao and Maggioni

D
d
Φ = {ϕi }m
i=1 ⊂ R if each data point x ∈ R may be expressed as a linear combination of
at most k elements of Φ. Clearly the case of interest is k  D. These sparse representations have been used in a variety of statistical signal processing tasks, compressed sensing,
learning (see e.g. Protter and Elad, 2007; Peyré, 2009; Lewicki et al., 1998; Kreutz-Delgado
et al., 2003; Maurer and Pontil, 2010; Chen et al., 1998; Donoho, 2006; Aharon et al., 2005;
Candes and Tao, 2007, among many others), and spurred much research about how to learn
data-adaptive dictionaries (see Gribonval et al., 2013; Vainsencher et al., 2011; Maurer and
Pontil, 2010, and references therein). The algorithms used in dictionary learning are often computationally demanding, being based on high-dimensional non-convex optimization
(Mairal et al., 2010). These approaches have the strength of being very general, with minimal assumptions made on geometry of the dictionary or on the distribution from which
the samples are generated. This “worst-case” approach incurs bounds dependent upon the
ambient dimension D in general (even in the standard case of data lying on one hyperplane).
In Maggioni et al. (2016) we proposed to attack the dictionary learning problem under
geometric assumptions on the data, namely that the data lies close to a low-dimensional set
M. There are of course various possible geometric assumptions, the simplest one being that
M is a single d-dimensional subspace. For this model Principal Component Analysis (PCA)
(see Pearson, 1901; Hotelling, 1933, 1936) is an effective tool to estimate the underlying
plane. More generally, one may assume that data lie on a union of several low-dimensional
planes instead of a single one. The problem of estimating multiple planes, called subspace
clustering, is more challenging (see Fischler and Bolles, 1981; Ho et al., 2003; Vidal et al.,
2005; Yan and Pollefeys, 2006; Ma et al., 2007, 2008; Chen and Lerman, 2009; Elhamifar
and Vidal, 2009; Zhang et al., 2010; Liu et al., 2010; Chen and Maggioni, 2011). This model
was shown effective in various applications, including image processing (Fischler and Bolles,
1981), computer vision (Ho et al., 2003) and motion segmentation (Yan and Pollefeys, 2006).
A different type of geometric model gives rise to manifold learning, where M is assumed
to be a d-dimensional manifold isometrically embedded in RD , see (Tenenbaum et al., 2000;
Roweis and Saul, 2000; Belkin and Niyogi, 2003; Donoho and Grimes, 2003; Coifman et al.,
2005a,b; Zhang and Zha, 2004) and many others. It is of interest to move beyond this model
to even more general geometric models, for example where the regularity of the manifold is
reduced, and data is not forced to lie exactly on a manifold, but only close to it.
Geometric Multi-Resolution Analysis (GMRA) was proposed in Allard et al. (2012), and
its finite-sample performance was analyzed in Maggioni et al. (2016). In GMRA, geometric
approximations of M are constructed with multiscale techniques that have their roots in
geometric measure theory, harmonic analysis and approximation theory. GMRA performs
a multiscale tree decomposition of data and build multiscale low-dimensional geometric
approximations to M. Given data, we run the cover tree algorithm (Beygelzimer et al.,
2006) to obtain a multiscale tree in which every node is a subset of M, called a dyadic cell,
and all dyadic cells at a fixed scale form a partition of M. After the tree is constructed,
we perform PCA on the data in each cell to locally approximate M by the d-dimensional
principal subspace so that every point in that cell is only encoded by the d coefficients in
principal directions. At a fixed scale M is approximated by a piecewise linear set. In Allard
et al. (2012) the performance of GMRA for volume measures on a C s , s ∈ (1, 2] Riemannian
manifold was analyzed in the continuous case (no sampling), and the effectiveness of GMRA
was demonstrated empirically on simulated and real-world data. In Maggioni et al. (2016),

2

Adaptive Geometric Multiscale Approximations for Intrinsically Low-dimensional Data

the approximation error of M was estimated in the non-asymptotic regime with n i.i.d.
samples from a measure ρ, satisfying certain technical assumptions, supported on a tube
of a C 2 manifold of dimension d isometrically embedded in RD . The probability bounds in
Maggioni et al. (2016) depend on n and d, but not on D, successfully avoiding the curse of
dimensionality caused by the ambient dimension. The assumption that ρ is supported in a
tube around a manifold can account for noise and does not force the data to lie exactly on
a smooth low-dimensional manifold.
In Allard et al. (2012) and Maggioni et al. (2016), GMRA approximations are constructed on uniform partitions in which all the cells have similar diameters. However, when
the regularity, such as smoothness or curvature, weighted by the ρ measure, of M varies
at different scales and locations, uniform partitions are not optimal. Inspired by the adaptive methods in classical multi-resolution analysis (see Donoho and Johnstone, 1994, 1995;
Cohen et al., 2002; Binev et al., 2005, 2007, among many others), we propose an adaptive version of GMRA to construct low-dimensional geometric approximations of M on an
adaptive partition and provide finite sample performance guarantee for a much larger class
of geometric structures M in comparison with Maggioni et al. (2016).
Our main result (Theorem 8) in this paper may be summarized as follows: Let ρ be a
probability measure supported on or near a compact d-dimensional Riemannian manifold
M ,→ RD , with d ≥ 3. Let ρ admit a multiscale decomposition satisfying the technical
assumptions A1-A5 in section 2.1 below. Given n i.i.d. samples are taken from ρ, the
intrinsic dimension d, and a parameter κ large enough, adaptive GMRA outputs a dictionary
b n = {φbi }i∈Jn , an encoding operator D
bn : RD → RJn and a decoding operator D
b −1 : RJn →
Φ
n
D
D
b
R . With high probability, for every x ∈ R , kDn xk0 ≤ d + 1 (i.e. only d + 1 entries are
non-zero), and the Mean Squared Error (MSE) satisfies
b −1 D
bn xk2 ] .
MSE := Ex∼ρ [kx − D
n



log n
n



2s
2s+d−2

.

Here s is a regularity parameter of ρ as in definition 5, which allows us to consider M’s
and ρ’s with nonuniform regularity, varying at different locations and scales. Note that the
algorithm does not need to know s, but it automatically adapts to obtain a rate that depends
on s. We believe, but do not prove, that this rate is indeed optimal. As for computational
b n takes O((C d + d2 )Dn log n) and computing D
bn x only takes
complexity, constructing Φ
2
O(d(D + d ) log n), which means we have a fast transform mapping data to their sparse
encoding on the dictionary.
In adaptive GMRA, the dictionary is composed of the low-dimensional planes on adaptive partitions and the encoding operator transforms a point to the local d principal coefficients of the data in a piece of the partition. We state this results in terms of encoding
and decoding to stress that learning the geometry in fact yields efficient representations of
data, which may be used for performing signal processing tasks in a domain where the data
admit a sparse representation, e.g. in compressive sensing or estimation problems (see Iwen
and Maggioni, 2013; Chen et al., 2012; Eftekhari and Wakin, 2015). Adaptive GMRA is
designed towards robustness, both in the sense of tolerance to noise and to model error (i.e.
data not lying on a manifold). We assume d is given throughout this paper. If not, we refer
to Little et al. (2012, 2009a,b) for the estimation of intrinsic dimensionality.
3

Liao and Maggioni

The paper is organized as follows. Our main results, including the construction of
GMRA, adaptive GMRA and their finite sample analysis, are presented in Section 2. We
show numerical experiments in Section 3. The detailed analysis of GMRA and adaptive
GMRA is presented in Section 4. In Section 5, we discuss the computational complexity of
adaptive GMRA and extend our work to adaptive orthogonal GMRA. Proofs are postponed
till the appendix.
Notation. We will introduce some basic notation here. f . g means that there exists a
constant C independent on any variable upon which f and g depend, such that f ≤ Cg.
Similarly for &. f  g means that both f . g and f & g hold. The cardinality of a set
A is denoted by #A. For x ∈ RD , kxk denotes the Euclidean norm and Br (x) denotes the
Euclidean ball of radius r centered at x. Given a subspace V ∈ RD , we denote its dimension
by dim(V ) and the orthogonal projection onto V by ProjV . If A is a linear operator on RD ,
||A|| is its operator norm. The identity operator is denoted by I.

2. Main results
GMRA was proposed in Allard et al. (2012) to efficiently represent points on or near a
low-dimensional manifold in high dimensions. We refer the reader to that paper for details
of the construction, and we summarize here the main ideas in order to keep the presentation
self-contained. The construction of GMRA involves the following steps:
(i) construct a multiscale tree T and the associated decomposition of M into nested cells
{Cj,k }k∈Kj ,j∈Z where j represents scale and k location;
(ii) perform local PCA on each Cj,k : let the mean (“center”) be cj,k and the d-dim principal
subspace Vj,k . Define P j,k (x) := cj,k + ProjV j,k (x − cj,k ).
0

(iii) construct a “difference” subspace Wj+1,k0 capturing P j,k (Cj,k ) − P j+1,k (Cj+1,k0 ), for
each Cj+1,k0 ⊆ Cj,k (these quantities are associated with the refinement criterion in
adaptive GMRA).
M may be approximated, at each scale j, by its projection PΛj onto the family of linear
sets Λj := {Pj,k (Cj,k )}k∈Kj . For example, linear approximations of the S manifold at scale
6 and 10 are displayed in Figure 1. In a variety of distances, PΛj (M) → M. In practice
M is unknown, and the construction above is carried over on training data, and its result
is random with the training samples. Naturally we are interested in the performance of the
construction on new samples. This is analyzed in a setting of “smooth manifold+noise”
in Maggioni et al. (2016). When the regularity (such as smoothness or curvature) of M
varies at different locations and scales, linear approximations on fixed uniform partitions
are not optimal. Inspired by adaptive methods in classical multi-resolution analysis (see
Cohen et al., 2002; Binev et al., 2005, 2007), we propose an adaptive version of GMRA
which learns adaptive and near-optimal approximations.
We will start with the multiscale tree decomposition in Section 2.1 and present GMRA
and adaptive GMRA in Section 2.3 and 2.4 respectively.
4

Adaptive Geometric Multiscale Approximations for Intrinsically Low-dimensional Data

j=6

(a) S manifold

(b) scale 6

j = 10

(c) scale 10

Figure 1: (a) S manifold; (b,c) Linear approximations at scale 6, 10.
2.1 Multiscale partitions and trees
A multiscale set of partitions of M with respect to probability measure ρ is a family of sets
{Cj,k }k∈Kj ,j∈Z , called dyadic cells, satisfying Assumptions (A1-A5) below for all integers
j ≥ jmin :
(A1) for any k ∈ Kj and k 0 ∈ Kj+1 , either Cj+1,k0 ⊆ Cj,k or ρ(Cj+1,k0 ∩ Cj,k ) = 0. We
denote the children of Cj,k by C (Cj,k ) = {Cj+1,k0 : Cj+1,k0 ⊆ Cj,k }. We assume that
amin ≤ #C (Cj,k ) ≤ amax . Also for every Cj,k , there exists a unique k 0 ∈ Kj−1 such
that Cj,k ⊆ Cj−1,k0 . We call Cj−1,k0 the parent of Cj,k .
(A2) ρ(M \ ∪k∈Kj Cj,k ) = 0, i.e. Λj := {Cj,k }k∈Kj is a cover for M.
(A3) ∃θ1 > 0 : #Λj ≤ 2jd /θ1 .
(A4) ∃θ2 > 0 such that, if x is drawn from ρ|Cj,k , then a.s. kx − cj,k k ≤ θ2 2−j .
j,k
j,k
(A5) Let λj,k
1 ≥ λ2 ≥ . . . ≥ λD be the eigenvalues of the covariance matrix Σj,k of ρ|Cj,k ,
defined in Table 1. Then:
−2j /d,
(i) ∃θ3 > 0 such that ∀j ≥ jmin and k ∈ Kj , λj,k
d ≥ θ3 2

j,k
(ii) ∃θ4 ∈ (0, 1) such that λj,k
d+1 ≤ θ4 λd .

(A1) implies that the {Cj,k }k∈Kj ,j≥jmin are associated with a tree structure, and with some
abuse of notation we call the above tree decompositions. (A1)-(A5) are natural assumptions,
easily satisfied by natural multiscale decompositions when M is a d-dimensional manifold
isometrically embedded in RD : see the work (Maggioni et al., 2016) for a detailed discussion.
(A2) guarantees that the cells at scale j form a partition of M; (A3) says that there are
at most 2jd /θ1 dyadic cells at scale j. (A4) ensures diam(Cj,k ) . 2−j . When M is a ddimensional manifold, (A5)(i) is the condition that the best rank d approximation to Σj,k
is close to the covariance matrix of a d-dimensional Euclidean ball, while (A5)(ii) imposes
that the (d+1)-th eigenvalue is smaller that the d-th eigenvalue, i.e. the set has significantly
larger variances in d directions than in all the remaining ones.
5

Liao and Maggioni

Linear projection
on Cj,k
Linear projection
at scale j

GMRA

Empirical GMRA

Pj,k (x) := cj,k + ProjVj,k (x − cj,k )

bj,k (x) := b
P
cj,k + ProjVbj,k (x − b
cj,k )

Pj :=

Measure

P

k∈Kj

bj := P
b
P
k∈Kj Pj,k 1j,k

Pj,k 1j,k

ρb(Cj,k ) = n
bj,k /n

ρ(Cj,k )

Center

cj,k := Ej,k x

Principal
subspaces

Vj,k minimizes
Ej,k kx − cj,k − ProjV (x − cj,k )k2
among d-dim subspaces

Covariance
matrix

Σj,k := Ej,k (x − cj,k )(x − cj,k )T

hPX, QXi
kPXk

´

M hPx, Qxidρ

´

M kPxk



2 dρ

P

xi ∈Cj,k

xi

Vbj,k minimizes
1 P
cj,k − ProjV (x − b
cj,k )k2
xi ∈Cj,k kx − b
n
bj,k
among d-dim subspaces
b j,k := 1 P
cj,k )(xi − b
cj,k )T
Σ
xi ∈Cj,k (xi − b
n
bj,k
1/n

1
2

1
n
bj,k

b
cj,k :=

P

1/n

xi ∈Xn hPxi , Qxi i

P

xi ∈Xn

kPxi k2

1
2

Table 1: 1j,k is the indicator function on Cj,k (i.e.,1j,k (x) = 1 if x ∈ Cj,k and 0 otherwise).
Here Ej,k stands for expectation with respect to the conditional distribution dρ|Cj,k .
The measure of Cj,k is ρ(Cj,k ) and the empirical measure is ρb(Cj,k ) = n
bj,k /n where
n
bj,k is the number of points in Cj,k . Vj,k and Vbj,k are the eigen-spaces associated
b j,k respectively. Here P, Q: M → RD
with the largest d eigenvalues of Σj,k and Σ
are two operators.
We will construct such {Cj,k }k∈Kj ,j≥jmin in Section 2.6. In our construction (A1-A4) is
satisfied when ρ a doubling probability measure1 (see Christ, 1990; Deng and Han, 2008).
If we further assume that M is a d-dimensional C s , s ∈ (1, 2] closed Riemannian manifold
isometrically embedded in RD , then (A5) is satisfied as well (See Proposition 13).
Some notation: a master tree T is associated with {Cj,k }k∈Kj ,j≥jmin (using property
(A1)), constructed on M; since Cj,k ’s at scale j have similar diameters, Λj := {Cj,k }k∈Kj is
called a uniform partition at scale j. It may happen that at the coarsest scales conditions
(A3)-(A5) are satisfied but with very poor constants θ: it will be clear that in all that
follows we may discard a few coarse scales, and only work at scales that are fine enough
and for which (A3)-(A5) truly capture the local geometry of M.
A proper subtree T̃ of T is a collection of nodes of T with the properties: (i) the root
node is in T̃ , (ii) if Cj,k is in T̃ then the parent of Cj,k is also in T̃ . Any finite proper
subtree T̃ is associated with a unique partition Λ = Λ(T̃ ) which consists of its outer leaves,
by which we mean those Cj,k ∈ T such that Cj,k ∈
/ T̃ but its parent is in T̃ .
2.2 Empirical GMRA
In practice the master tree T is not given, nor can be constructed since M is not known:
we will construct one on samples by running a variation of the cover tree algorithm (see
Beygelzimer et al., 2006). From now on we denote the training data by X2n . We randomly
1. ρ is doubling if there exists C1 > 0 such that C1−1 rd ≤ ρ(M ∩ Br (x)) ≤ C1 rd for any x ∈ M and r > 0.
C1 is called the doubling constant of ρ.

6

Adaptive Geometric Multiscale Approximations for Intrinsically Low-dimensional Data

split the data into two disjoint groups such that X2n = Xn0 ∪ Xn where Xn0 = {x01 , . . . , x0n }
and Xn = {x1 , . . . , xn }, apply the cover tree algorithm (see Beygelzimer et al., 2006) on Xn0
to construct a tree satisfying (A1-A5) (see section 2.6). After the tree is constructed, we
assign points in the second half of data Xn , to the appropriate cells. In this way we obtain a
family of multiscale partitions for the points in Xn , which we truncate to the largest subtree
whose leaves contain at least d points in Xn . This subtree is called the data master tree,
denoted by T n . We then use Xn to perform local PCA to obtain the empirical mean b
cj,k
b
and the empirical d-dimensional principal subspace Vj,k on each Cj,k . Define the empirical
bj,k (x) := b
projection P
cj,k + ProjVbj,k (x − b
cj,k ) for x ∈ Cj,k . Table 1 summarizes the GMRA
objects and their empirical counterparts.
2.3 Geometric Multi-Resolution Analysis: uniform partitions
GMRA with respect to the distribution ρ associated with the multiscale tree T consists a
collection of piecewise affine projectors {Pj : RD → RD }j≥jmin on the multiscale partitions
{Λj := {Cj,k }k∈Kj }j≥jmin . At scale j, M is approximated by the piecewise linear sets
{Pj,k (Cj,k )}k∈Kj . In order to understand the empirical approximation error of M by the
bj,k (Cj,k )}k∈K at scale j, we split the error into a squared bias term
piecewise linear sets {P
j
and a variance term:
bj Xk2 = kX − Pj Xk2 + EkPj X − P
bj Xk2 .
EkX − P
|
{z
} |
{z
}

(1)

variance

bias2

bj Xk2 is also called the Mean Square Error (MSE) of GMRA. To bound the bias
EkX − P
term, we need regularity assumptions on the ρ, and for the variance term we prove concentration bounds of the relevant quantities around their expected values.
For a fixed distribution ρ, the approximation error of M at scale j, measured by kX −
Pj Xk := k(I − Pj )Xk, decays at a rate dependent on the regularity of M in the ρ-measure
(see Allard et al., 2012). We quantify the regularity of ρ as follows:
Definition 1 (Model class As ) A probability measure ρ supported on M is in As if
|ρ|As = sup inf{A0 : kX − Pj Xk ≤ A0 2−js , ∀ j ≥ jmin } < ∞ ,

(2)

T

where T varies over the set, assumed non-empty, of multiscale tree decompositions satisfying
Assumption (A1-A5).
We capture the case where the L2 approximation error is roughly the same on every cell
with the following definition:
∞
Definition 2 (Model class A∞
s ) A probability measure ρ supported on M is in As if
q
−js
|ρ|A∞
=
sup
inf{A
:
k(X
−
P
X)1
k
≤
A
2
ρ(Cj,k ), ∀ k ∈ Kj , j ≥ jmin } < ∞ (3)
0
0
j,k
j,k
s
T

where T varies over the set, assumed non-empty, of multiscale tree decompositions satisfying
Assumption (A1-A5).
7

Liao and Maggioni

−j
Clearly
A∞
s ⊂ As . Also, since diam(Cj,k ) ≤ 2θ2 2 , necessarily k(I − Pj,k )1j,k Xk ≤
p
θ2
ρ(Cj,k ), ∀ k ∈ Kj , j ≥ jmin , and therefore ρ ∈ A∞
1 in any case. Finally, these classes
contain suitable measures supported on manifolds:

2−j

Proposition 3 Let M be a closed manifold of class C s , s ∈ (1, 2] isometrically embedded
in RD , and ρ be a doubing probability measure on M with the doubling constant C1 . Then
our construction of {Cj,k }k∈Kj ,j≥jmin in Section 2.6 satisfies (A1-A5), and ρ ∈ A∞
s .
Proposition 3 is proved in Appendix A.2.
Example 1 We consider the d-dim S manifold whose x1 and x2 coordinates are on the
S curve and xi ∈ [0, 1], i = 3, . . . , d + 1. As stated above, the volume measure on the S
manifold is in A∞
2 . Numerically one can identify s from data sampled from ρ ∈ As as the
slope of the line approximating log10 kX − Pj Xk as a function of log10 rj where rj is the
average diameter of Cj,k ’s at scale j. Our numerical experiments in Figure 5 (b) give rise
to s ≈ 2.0239, 2.1372, 2.173 when d = 3, 4, 5 respectively.
Example 2 For comparison we consider the d-dimensional Z manifold whose x1 and x2
coordinates are on the Z curve and xi ∈ [0, 1], i = 3, . . . , d + 1. Volume measure on the Z
manifold is in A1.5 (see appendix B.2). Our numerical experiments in Figure 5 (c) give rise
to s ≈ 1.5367, 1.6595, 1.6204 when d = 3, 4, 5 respectively.
The squared bias in (1) satisfies kX − Pj Xk2 ≤ |ρ|2As 2−2js whenever ρ ∈ As . In Proposition 15 we will show that the variance term is estimated in terms of the sample size n and
the scale j:
!
2 #Λ log[αd#Λ ]
j(d−2)
d
j2
j
j
2
bj Xk ≤
EkPj X − P
=O
,
β22j n
n

where α, β are constants depending on θ2 , θ3 . In the case d = 1 both the squared bias and
the variance decrease as j increases, so choosing the finest scale of the data tree T n yields
the best rate of convergence. When d ≥ 2, the squared bias decreases but the variance
increases as j gets large as shown Figure 2, as a manifestation of the bias-variance tradeoff
in classical statistical and machine learning, except that it arises here in a geometric setting.
By choosing a proper scale j to balance these two terms, we obtain the following rate of
convergence for empirical GMRA:
Theorem 4 Suppose ρ ∈ As for s ≥ 1. Let ν > 0 be arbitrary and µ > 0. Let j ∗ be chosen
such that

µ logn n
for d = 1
∗
−j

 1
2
=
(4)
µ log n 2s+d−2 , for d ≥ 2
n
then there exists C1 := C1 (θ1 , θ2 , θ3 , θ4 , d, ν, µ) and C2 := C2 (θ1 , θ2 , θ3 , θ4 , d, µ) such that:


log n
s
b
P kX − Pj ∗ Xk ≥ (|ρ|As µ + C1 )
≤ C2 n−ν , for d = 1,
(5)
n
(

 s )
log n 2s+d−2
s
b
P kX − Pj ∗ Xk ≥ (|ρ|As µ + C1 )
≤ C2 n−ν , for d ≥ 2 .
(6)
n
8

Adaptive Geometric Multiscale Approximations for Intrinsically Low-dimensional Data

Total error

log 10(error) versus log10(partition size)

-0.4

Uniform
Adaptive

-0.6

-0.8

Optimal scale j *

log 10(L2 error)

-1

Variance

-1.2

-1.4

-1.6

-1.8

-2

Squared bias
-2.2
0

scale j

0.5

1

1.5

2

2.5

3

3.5

4

log 10(partition size)

(a) Bias and variance tradeoff

(b) Error versus the partition size

Figure 2: (a) Plot of the bias and variance estimates in Eq. (1), with s = 2, d = 5, n = 100.
(b) shows the approximation error on test data versus the partition size in GMRA
and adaptive GMRA for the 3-dim S manifold. When the partition size is between
1 and 102.8 , the bias dominates the error so the error decreases; after that, the
variance dominates the error, which becomes increasing.

Theorem 4 is proved in Section 4.2. In the perspective of dictionary learning, GMRA prod
vides a dictionary Φj ∗ of cardinality  dn/ log n for d = 1 and of cardinality  d(n/ log n) 2s+d−2
for d ≥ 2, so that every x sampled from ρ (and not just samples in the training data) may be
encoded with a vector with d + 1 nonzero entries: one entry encodes the location k of x on
the tree, e.g. (j ∗ , x) = (j ∗ , k) such that x ∈ Cj ∗ ,k , and the other d entries are VbjT∗ ,x (x−b
cj ∗ ,x ).
We also remind the reader that GMRA automatically constructs a fast transform mapping
points x to the vector representing Φj ∗ (See Allard et al. (2012); Maggioni et al. (2016) for
a discussion). Note that by choosing ν large enough,

 2s
2s+d−2
log
n
2
bj ∗ Xk .
(6) =⇒ MSE = EkX − P
,
n
and (5) implies MSE . ( logn n )2 for d = 1. Clearly, one could fix a desired MSE of size ε2 ,
and obtain from the above a dictionary of size dependent only of ε and independent of n,
for n sufficiently large, thereby obtaining a way of compressing data (see Maggioni et al.
(2016) for further discussion on this point). A special case of Theorem 4 with s = 2 was
proved in Maggioni et al. (2016).
2.4 Geometric Multi-Resolution Analysis: Adaptive Partitions
The performance guarantee in Theorem 4 is not fully satisfactory for two reasons: (i)
the regularity parameter s is required to be known to choose the optimal scale j ∗ , and
this parameter is typically unknown in any practical setting, and (ii) none of the uniform
partitions {Cj,k }k∈Kj will be optimal if the regularity of ρ (and/or M) varies at different
locations and scales. This lack of uniformity in regularity appears in a wide variety of
data sets: when clusters exist that have cores denser than the remaining regions of space,
9

Liao and Maggioni

Figure 3: Left: a master tree in which red nodes satisfy ∆j,k ≥ 2−j τn but blue nodes do not.
Right: the subtree of the red nodes is the smallest proper subtree that contains
all the nodes satisfying ∆j,k ≥ 2−j τn , i.e. were red in the figure on the left. Green
nodes form the adaptive partition.

when trajectories of a dynamical system are sampled that linger in certain regions of space
for much longer time intervals than others (e.g. metastable states in molecular dynamics
(Rohrdanz et al., 2011; Zheng et al., 2011)), in data sets of images where details exist at
different level of resolutions, affecting regularity at different scales in the ambient space, and
so on. To fix the ideas we consider again one simplest manifestations of this phenomenon
in the examples considered above: uniform partitions work well for the volume measure on
the S manifold but are not optimal for the volume measure on the Z manifold, for which the
ideal partition is coarse on flat regions but finer at and near the corners (see Figure 4). In
applications, for example to mesh approximation, it is often the case that the point clouds to
be approximated are not uniformly smooth and include different levels of details at different
locations and scales (see Figure 8). Thus we propose an adaptive version of GMRA that
will automatically adapts to the regularity of the data and choose a near-optimal adaptive
partition.
Adaptive partitions may be effectively selected with a criterion that determines whether
or not a cell should participate to the adaptive partition. The quantities involved in the
selection and their empirical version are summarized in Table 2.
We expect ∆j,k to be small on approximately flat regions, and large ∆j,k at many scales
b j,k to have the same behavior, at least when ∆
b j,k is
at irregular locations. We also expect ∆
with high confidence close to ∆j,k . We see this phenomenon represented in Figure 4 (a,b):
bj+1 xi − P
bj xi k decays uniformly at all points, while for
as j increases, for the S manifold kP
the Z manifold, the same quantity decays rapidly on flat regions but remains large at fine

Difference operator
Norm of difference

Definition (infinite sample)

Empirical version

Qj,k := (Pj − Pj+1 )1j,k
´
∆2j,k := Cj,k kQj,k xk2 dρ

bj,k := (P
bj − P
bj+1 )1j,k
Q
2
b 2 := 1 P
b
∆
xi ∈Cj,k kQj,k xi k
j,k
n

Table 2: Refinement criterion and the empirical version
10

Adaptive Geometric Multiscale Approximations for Intrinsically Low-dimensional Data

Algorithm 1 Empirical adaptive GMRA
Input: data X2n = Xn0 ∪ Xn , intrinsic dimension d, threshold κ
bb : multiscale tree, corresponding cells and adaptive piecewise
Output: T n , {Cj,k }, P
Λτn
linear projectors on adaptive partition.
1: Construct T n and {Cj,k } from Xn0
bj,k and ∆
b j,k on every node Cj,k ∈ T n .
2: Now use Xn . Compute P
b j,k ≥ 2−j τn where
smallest proper subtree of T n containing all Cj,k ∈ T n : ∆
3: Tbτn ← p
τn = κ (log n)/n.
b τn ← partition associated with outer leaves of Tbτn
4: Λ
b
bb ← P
5: P
b τ Pj,k 1j,k .
Λτ
Cj,k ∈Λ
n

n

scales around the corners. We wish to include in our approximation the nodes where this
quantity is large, since we may expect a large improvement in approximation by including
such nodes. However if too few samples exist in a node, then this quantity is not to be
trusted, because its variance is large. It turns out that it is enough to consider the following
criterion: let Tbτn be the smallestpproper subtree of T n that contains all C j,k ∈ T n for
b j,k ≥ 2−j τn where τn = κ (log n)/n. Crucially, κ may be chosen independently of
which ∆
the regularity index (see Theorem 8). Empirical adaptive GMRA returns piecewise affine
b τn , the partition associated with the outer leaves of Tbτn . Our algorithm is
projectors on Λ
summarized in Algorithm 1.
We will provide a finite sample performance guarantee of the empirical adaptive GMRA
for a model class that is more general than A∞
s . Given any fixed threshold η > 0, we let
T(ρ,η) be the smallest proper tree of T that contains all Cj,k ∈ T for which ∆j,k ≥ 2−j η.
The corresponding adaptive partition Λ(ρ,η) consists of the outer leaves of T(ρ,η) . We let
#j T(ρ,η) be the number of cells in T(ρ,η) at scale j.
Definition 5 (Model class Bs ) In the case d ≥ 3, given s > 0, a probability measure ρ
supported on M is in Bs if ρ satisfies the following regularity condition:
|ρ|pBs := sup sup η p
T

η>0

X

j≥jmin

2−2j #j T(ρ,η) < ∞, with p =

2(d − 2)
2s + d − 2

(7)

where T varies over the set, assumed nonempty, of multiscale tree decompositions satisfying
Assumption (A1-A5).
For elements in the model class Bs we have control on the growth rate of the truncated
tree T(ρ,η) as η decreases, namely it is O(η −p ). The key estimate on variance and sample
complexity in Lemma 14 indicates thatPthe natural measure of the complexity of T(ρ,η) is
the weighted tree complexity measure j≥jmin 2−2j #j T(ρ,η) in the definition above. First
of all, the class Bs is indeed larger than A∞
s (see appendix A.4 for a proof):
∞
Lemma 6 Bs is a more general model class than A∞
s . If ρ ∈ As , then ρ ∈ Bs and
|ρ|Bs . |ρ|A∞
.
s

11

Liao and Maggioni

log 10||Pj+1 x-P j x||2

log 10||Pj+1 x-P j x||2

2

0

5

2

d=3 Uniform:
d=3 Adaptive:
d=4 Uniform:
d=4 Adaptive:
d=5 Uniform:
d=5 Adaptive:

-0.5

0

5

log 10error versus log10(partition size)

0

slope= -0.78825 theory= -0.66667
slope= -0.76836 theory= -0.66667
slope= -0.59324 theory= -0.5
slope= -0.61705 theory= -0.5
slope= -0.48404 theory= -0.4
slope= -0.49943 theory= -0.4

-1

-2

-2

Scale

Scale

-4

-4

15

15
-6

log 10(L2 error)

10

10

-1.5

-2

-6

-2.5

20

20

-8

-8

25
-10

1

2

3

4

5

6

7

Point index

8

9

10

(a) S manifold

-3.5

-10

1

2

3

4

5

6

7

8

9

0

10

0.5

1

1.5

2

(b) Z manifold

2.5

3

3.5

4

4.5

log 10(partition size)

#10 4

Point index

#10 4

-3

(c) Error vs partition size, S manifold
log 10error versus log10(partition size)

0

d=3 Uniform:
d=3 Adaptive:
d=4 Uniform:
d=4 Adaptive:
d=5 Uniform:
d=5 Adaptive:

-0.5

slope= -0.58149 theory= -0.5
slope= -0.77297 theory= -0.75
slope= -0.4914 theory= -0.375
slope= -0.64198 theory= -0.5
slope= -0.40358 theory= -0.3
slope= -0.48847 theory= -0.375

log 10(L2 error)

-1

-1.5

-2

-2.5

-3

-3.5
0

0.5

1

1.5

2

2.5

3

3.5

4

log 10(partition size)

(d) Adaptive approx. of S

(e) Adaptive approx. of Z

(f) Error vs partition size, Z manifold

bj (xi ) − P
bj+1 (xi )|| from the coarsest scale (top) to the finest scale
Figure 4: (a,b): log10 ||P
(bottom), with columns indexed by points, sorted from left to right on the manifold. (d,e): adaptive approximations. (c,f): log-log plot of the approximation
error versus the partition size in GMRA and adaptive GMRA respectively. Theoretically, the slope is −2/d in both GMRA and adaptive GMRA for the S manifold. For the Z manifold, the slope is −1.5/d in GMRA and −1.5/(d − 1) in
adaptive GMRA (see appendix B).

Example 3 The volume measures on the d-dim (d ≥ 3) S manifold and Z manifold are in
B2 and B1.5(d−2)/(d−3) respectively (see appendix B). Numerically s is approximated by the
bb Xkd−2 versus the weighted complexity
negative of the slope in the log-log plot of kX − P
Λη
of the truncated tree according to Eq. (9). See numerical examples in Figure 5.

We also need a quasi-orthogonality condition which says that the operators {Qj,k }k∈Kj ,j≥jmin
applied on M are mostly orthogonal across scales and/or kQj,k Xk quickly decays.
Definition 7 (Quasi-orthogonality) There exists a constant B0 > 0 such that for any
proper subtree T̃ of any mater tree T satisfying Assumption (A1-A5),
X
X
k
Qj,k Xk2 ≤ B0
kQj,k Xk2 .
(8)
Cj,k ∈
/ T̃

Cj,k ∈
/ T̃

We postpone further discussion of this condition to Section 5.2. One can show (see appendix
D.1) that in the case d ≥ 3, ρ ∈ Bs along with quasi-orthogonality implies a certain rate of
12

Adaptive Geometric Multiscale Approximations for Intrinsically Low-dimensional Data

approximation of X by PΛ(ρ,η) X, as η → 0+ :


kX − PΛ(ρ,η) Xk2 ≤ Bs,d |ρ|pBs η 2−p ≤ Bs,d |ρ|2Bs 

X

j≥jmin

−

2−2j #j T(ρ,η) 

2s
d−2

,

(9)

where s = (d−2)(2−p)
and Bs,d := B0 2p /(1 − 2p−2 ).
2p
The main result of this paper is the following performance analysis of empirical adaptive
GMRA (see the proof in Section 4.3).
Theorem 8 Suppose ρ satisfies quasi-orthogonality and M isp
bounded: M ⊂ BM (0). Let
ν > 0. There exists κ0 (θ2 , θ3 , θ4 , amax , d, ν) such that if τn = κ (log n)/n with κ ≥ κ0 , the
following holds:
(i) if d ≥ 3 and ρ ∈ Bs for some s > 0, there are c1 and c2 such that
(
 s )

2s+d−2
log
n
bb Xk ≥ c1
≤ c2 n−ν .
P kX − P
Λτn
n

(ii) if d = 1, there exist c1 and c2 such that
(

1 )
2
log
n
bb Xk ≥ c1
P kX − P
≤ c2 n−ν .
Λτn
n

(iii) if d = 2 and

|ρ| := sup sup[− log η]−1
T

η>0

X

j≥jmin

(10)

(11)

2−2j #j T(ρ,η) < +∞ ,

then there exist c1 and c2 such that
(
 2  21 )
bb Xk ≥ c1 log n
≤ c2 n−ν .
P kX − P
Λτn
n

(12)

The dependencies of the constants in Theorem 8 on the geometric constants are as follows:
d≥3:
d=1:
d=2:

c1 = c1 (θ2 , θ3 , θ4 , amax , d, s, κ, |ρ|Bs , B0 , ν),
c1 = c1 (θ1 , θ2 , θ3 , θ4 , amax , d, κ, B0 , ν),
c1 = c1 (θ2 , θ3 , θ4 , amax , d, κ, |ρ|, B0 , ν),

c2 = c2 (θ2 , θ3 , θ4 , amin , amax , d, s, κ, |ρ|Bs , B0 )
c2 = c2 (θ1 , θ2 , θ3 , θ4 , amin , amax , d, κ, B0 )
c2 = c2 (θ2 , θ3 , θ4 , amin , amax , d, κ, |ρ|, B0 )

Notice that by choosing ν large enough, we have
(
 α β )
 α 2β
log n
log n
−ν
b
P kX − PΛb τ Xk ≥ c1
≤ c2 n ⇒ MSE ≤ c1
,
n
n
n
2s

so we also have the MSE bound: MSE . (log n/n) 2s+d−2 for d ≥ 3 and MSE . logd n/n for
d = 1, 2.
13

Liao and Maggioni

In comparison with Theorem 4, Theorem 8 is more satisfactory for two reasons: (i) when
2s
d ≥ 3, the same rate (log n/n) 2s+d−2 is proved for the model class Bs which is larger than
A∞
s . (ii) our algorithm is universal: it does not require a priori knowledge of the regularity
s, since the choice of κ is independent of s, yet it achieves the rate as if it knew the optimal
regularity parameter s.
In the perspective of dictionary learning, when d ≥ 3, adaptive GMRA provides a
d−2
dictionary ΦΛb τ associated with a tree of weighted complexity (n/ log n) 2s+d−2 , so that
n
every x sampled from ρ may be encoded by a vector with d + 1 nonzero entries, among
which one encodes the location of x in the adaptive partition and the other d entries are
the local principal coefficients of x.
For a given accuracy ε, in order to achieve MSE . ε2 , the number of samples we need
2s+d−2
is nε & (1/ε) s log(1/ε). When s is unknown, we can determine s as follows: we fix
a small n0 and run adaptive GMRA with n0 , 2n0 , 4n0 , . . . , 10n0 samples. For each sample
size, we evenly split data into the training set to construct adaptive GMRA and the test set
2s
to evaluate the MSE. According to Theorem 8, the MSE scales like [(log n)/n] 2s+d−2 where
n is the sample size. Therefore, the slope in the log-log plot of the MSE versus n gives an
approximation of −2s/(2s + d − 2).
The threshold τn in our adaptive algorithm is independent of s as κ0 does not depend on
s, which means our adaptive algorithm does not require s as a priori information but rather
will learn it from
also be natural to consider another stopping criterion:
´ data. It would
1
2 :=
2 dρ ≤ η 2 which suggests stopping refinement to finer scales
Ej,k
kP
x
−
xk
j
ρ(Cj,k ) Cj,k
if the approximation error is below certain threshold. The reason why we do not adopt
this stopping criterion is that in this case the threshold η would have to depend on s in
2s
order to guarantee the (adaptive) rate MSE . (log n/n) 2s+d−2 for d ≥ 3. More precisely,
E
for any threshold η > 0, let T(ρ,η)
be the smallest proper subtree of T whose leaves satisfy
2
E .
2
Ej,k ≤ η . The corresponding adaptive partition ΛE(ρ,η) consists of the leaves of T(ρ,η)
This stopping criterion guarantees kX − PΛE Xk ≤ η. It is natural to define the model
(ρ,η)
class Fs in the case d ≥ 3 to be the set of probability measures ρ supported on M such
d−2 P
−2j # ΛE
that supT supη>0 η s
j (ρ,η) < ∞ where T varies over the set of multiscale
j≥jmin 2
tree decompositions satisfying (A1-A5). One can show that A∞
s ( Fs . As an analogue of
Theorem 8, we can prove that, there exists κ0 > 0 such that if our adaptive algorithm adopts
s
the stopping criterion Ebj,k ≤ τnE where the threshold is chosen as τnE = κ (log n/n) 2s+d−2
b τ E satisfies
with κ ≥ κ0 , then the empirical approximation on the adaptive partition Λ
n
2s
bb Xk2 . (log n/n) 2s+d−2 . With this stopping criterion, the threshold τnE
MSE = kX − P
Λτ E
n

would depend on s, forcing us to know s as a priori information, unlike in Theorem 8.
Theorem 8 is stated when M is bounded. The assumption of the boundedness of M is
largely irrelevant, and may be replaced by a weaker assumption on the decay of ρ.
Theorem 9 Let d ≥ 3, s, δ, λ, µ > 0. Assume that there exists C1 such that
ˆ

BR (0)c

||x||2 dρ ≤ C1 R−δ , ∀R ≥ R0 .
14

Adaptive Geometric Multiscale Approximations for Intrinsically Low-dimensional Data

Suppose ρ satisfies quasi-orthogonality. If ρ restricted on BR (0), denoted by ρ|BR (0) , is in
2(d−2)
Bs for every R ≥ R0 and |ρ|BR (0) |pBs ≤ C2 Rλ for some C2 > 0, where p = 2s+d−2
. Then
p
there exists κ0 (θ2 , θ3 , θ4 , amax , d, ν) such that if τn = κ log n/n with κ ≥ κ0 , we have
)
(
δ

 s
2s+d−2 δ+max(λ,2)
log
n
bb Xk ≥ c1
≤ c2 n−ν
(13)
P kX − P
Λτn
n

bb X is obtained by adaptive GMRA
for some c1 , c2 independent of n, where the estimator P
Λτ
2s

n

within BRn (0) where Rn = max(R0 , µ(n/ log n) (2s+d−2)(δ+max(λ,2)) ), and is equal to 0 for the
points outside BRn (0).
2s

·

δ

Theorem 9 is proved at the end of Section 4.3. It implies MSE . (log n/n) 2s+d−2 δ+max(λ,2)
2s
As δ increases, i.e., δ → +∞, the MSE approaches (log n/n) 2s+d−2 , which is consistent with
Theorem 8 for bounded M. Similar results, with similar proofs, would hold under different
assumptions on the decay of ρ; for example for ρ decaying exponentially, or faster, only
log n terms in the rate would be lost compared to the rate in Theorem (8).
Remark 10 We claim that λ is not large in simple cases. For example, if ρ ∈ A∞
s and ρ
−jd
−(d+1+δ)
decays in the radial direction in such a way that ρ(Cj,k ) ≤ C2 kcj,k k
, it is easy
to show that ρ|BR (0) ∈ Bs for all R > 0 and |ρ|BR (0) |pBs ≤ Rλ with λ = d − (d+1+δ)(d−2)
(see
2s+d−2
the end of Section 4.3).
Remark 11 Suppose that ρ was supported in a tube of radius σ around a d-dimensional
manifold M, a model that can account both for (bounded) noise and situations where data
is not exactly on a manifold, but close to it, as in Maggioni et al. (2016). Then Theorem 8
and Theorem 9 apply in this case, provided one stops at scale j such that 2−j & σ.
Remark 12 In these Theorems we are assuming that d is given because it can be estimated
using existing techniques, see for example Little et al. (2012) and references therein.
2.5 Connection to previous works
The works by Allard et al. (2012) and Maggioni et al. (2016) are natural predecessors to
this work. In Allard et al. (2012), GMRA and orthogonal GMRA were proposed as datadriven dictionary learning tools to analyze intrinsically low-dimensional point clouds in a
high dimensions. The bias kX − Pj Xk were estimated for volume measures on C s , s ∈ (1, 2]
manifolds . The performance of GMRA, including sparsity guarantees and computational
costs, were systematically studied and tested on both simulated and real data. In Maggioni
et al. (2016) the finite sample behavior of empirical GMRA was studied. A non-asymptotic
bj Xk for the model class A2 (a special
probabilistic bound on the approximation error kX − P
case of Theorem 4 with s = 2) was established. It was further proved that if the measure
ρ is absolutely continuous with respect to the volume measure on a tube of a bounded C 2
manifold with a finite reach, then ρ is in A2 . Running the cover tree algorithm on data
gives rise to a family of multiscale partitions satisfying Assumption (A3-A5). The analysis
in Maggioni et al. (2016) robustly accounts for noise and modeling errors as the probability
15

Liao and Maggioni

measure is concentrated “near” a manifold. This work extends GMRA by introducing
adaptive GMRA, where low-dimensional linear approximations of M are built on adaptive
partitions at different scales. The finite sample performance of adaptive GMRA is proved
for a large model class. Adaptive GMRA takes full advantage of the multiscale structure of
GMRA in order to model data sets of varying complexity across locations and scales. We
also generalize the finite sample analysis of empirical GMRA from A2 to As , and analyze
the finite sample behavior of orthogonal GMRA and adaptive orthogonal GMRA.
In a different direction, a popular learning algorithm for fitting low-dimensional planes
to data is k-flats: let Fk be the collections of k flats (affine spaces) of dimension d. Given
data Xn = {x1 , . . . , xn }, k-flats solves the optimization problem
n

min

S∈Fk

1X
dist2 (xi , S)
n

(14)

i=1

where dist(x, S) = inf y∈S kx−yk. Even though a global minimizer of (14) exists, it is hard to
attain due to the non-convexity of the model class Fk , and practitioners are aware that many
local minima that are significantly worse than the global minimum exist. While often k is
considered given, it may be in fact chosen from the data: for example Theorem 4 in Canas
et al. (2012) implies that, given n samples from a probability measure that is absolutely
continuous with respect to the volume measure on a smooth d-dimensional manifold M, the
d
expected (out-of-sample) L2 approximation error of M by kn = C1 (M, ρ)n 2(d+4) planes is of
2
order O(n− d+4 ). This result is comparable with our Theorem 4 in the case s = 2 which says
1
that the L2 error by empirical GMRA at the scale j such that 2j  (n/ log n) d+2 achieves
2
a faster rate O(n− d+2 ). So we not only achieve a better rate, but we do so with provable
and fast algorithms, that are nonlinear but do not require non-convex optimization.
Multiscale adaptive estimation has been an intensive research area for decades. In the
pioneering works by Donoho and Johnstone (see Donoho and Johnstone, 1994, 1995), soft
thresholding of wavelet coefficients was proposed as a spatially adaptive method to denoise a
function. In machine learning, Binev et al. addressed the regression problem with piecewise
constant approximations (see Binev et al., 2005) and piecewise polynomial approximations
(see Binev et al., 2007) supported on an adaptive subpartition chosen as the union of dataindependent cells (e.g. dyadic cubes or recursively split samples). While the works above are
in the context of function approximation/learning/denoising, a whole branch of geometric
measure theory (following the seminal work by Jones (1990); David and Semmes (1993))
quantifies via multiscale least squares fits the rectifiability of sets and their approximability
by multiple images of bi-Lipschitz maps of, say, a d-dimensional square. We can the view
the current work as extending those ideas to the setting where data is random, possibly
noisy, and guarantees on error on future data become one of the fundamental questions.
Theorem 8 can be viewed as a geometric counterpart of the adaptive function approximation in Binev et al. (2005, 2007). Our results are a “geometric counterpart” of sorts.
We would like to point out two main differences between Theorem 8 and Theorem 3 in
Binev et al. (2005): (i) In Binev et al. (2005, Theorem 3), there is an extra assumption that
the function is in Aγ with γ arbitrarily small. This assumption takes care of the error at
the nodes in T \ T n where the thresholding criteria would succeed: these nodes should be
added to the adaptive partition but have not been explored by our data. This assumption
16

Adaptive Geometric Multiscale Approximations for Intrinsically Low-dimensional Data

is removed in our Theorem 8 by observing that the nodes below the data master tree have
small measure so their refinement criterion is smaller than 2−j τn with high probability.
b j,k ≥ 2−j τn unlike the criterion
(ii) we consider scale-dependent thresholding criterion ∆
in Binev et al. (2005, 2007) that is scale-independent. This difference arises because at
scale j our linear approximation is built on data within a ball of radius . 2−j and so the
variance of PCA on a fixed cell at scale j is proportional to 2−2j . For the same reason,
we measure the complexity of T(ρ,η) in terms of the weighted tree complexity instead of the
cardinality since the former one gives an upper bound of the variance in piecewise linear
approximation on partition via PCA (see Lemma 14). Using scale-dependent threshold and
measuring tree complexity in this way give rise to the best rate of convergence. In contrast, if we use scale-independent threshold and define a model class Γs for whose elements
2d
#T(ρ,η) = O(η − 2s+d ) (analogous to the function class in Binev et al. (2005, 2007)), we can
2s

2s+d . However
still show that A∞
s ⊂ Γs , but the estimator only achieves MSE . ((log n)/n)
0
2(d−2)
2d
2
∞
s
many elements of Γs not in As are in B with 2s0 +d−2 = 2s+d , and in Theorem 8 the
estimator based on scaled thresholding achieves a better rate, which we believe is optimal.
We refer the reader to Maggioni et al. (2016) for a thorough discussion of further related
work related to manifold and dictionary learning.

2.6 Construction of a multiscale tree decomposition
Our multiscale tree decomposition is constructed from a variation of the cover tree algorithm
(see Beygelzimer et al., 2006) applied on half of the data denoted by Xn0 . In brief the cover
tree T (Xn0 ) on Xn0 is a leveled tree where each level is a “cover” for the level beneath it. Each
level is indexed j and each node in T (Xn0 ) is associated with a point in Xn0 . A point can be
associated with multiple nodes in the tree but it can appear at most once at every level.
Let Tj (Xn0 ) ⊂ Xn0 be the set of nodes of T at level j. The cover tree obeys the following
properties for all j ∈ [jmin , jmax ]:
1. Nesting: Tj (Xn0 ) ⊂ Tj+1 (Xn0 );
2. Separation: for all distinct p, q ∈ Tj (Xn0 ), kp − qk > 2−j ;
3. Covering: for all q ∈ Tj+1 (Xn0 ), there is p ∈ Tj (Xn0 ) such that kp − qk < 2−j . The node
at level j associated with p is a parent of the node at level j + 1 associated with q.
In the third property, q is called a child of p. Each node can have multiple parents but
is only assigned to one of them in the tree. The properties above imply that for any q ∈ Xn0 ,
there exists p ∈ Tj such that kp − qk < 2−j+1 . The authors in Beygelzimer et al. (2006)
showed that cover tree always exists and the construction takes time O(C d Dn log n) .
We know show that from a set of nets {Tj (Xn0 )}j∈[jmin ,jmax ] as above we can construct a
set of Cj,k with desired properties. (see Appendix A for the construction of Cj,k ’s and the
f defined in (31) is equal to the union of the C 0 s constructed
proof of Proposition 13). M
j,k
above up to a set whose empirical measure is 0.
2. For these elements, the average cell-wise refinement is monotone such that: for every Cj,k and Cj+1,k0 ⊂
∆
0
∆
Cj,k , √ j+1,k
≤ √ j,k .
ρ(Cj+1,k0 )

ρ(Cj,k )

17

Liao and Maggioni

Proposition 13 Assume ρ is a doubling probability measure on M with the doubling constant C1 . Then {Cj,k }k∈Kj ,jmin ≤j≤jmax constructed above satisfies the Assumptions
1. (A1) with amax = C12 (24)d and amin = 1.
2. For any ν > 0,


f > 28ν log n
P ρ(M \ M)
3n



≤ 2n−ν ;

(15)

3. (A3) with θ1 = C1−1 4−d ;
4. (A4) with θ2 = 3.
5. If additionally
5a. if ρ satisfies the conditions in (A5) with Br (z), z ∈ M, replacing Cj,k with
constants θ̃3 , θ̃4 such that λd (Cov(ρ|Br (z) )) ≥ θ̃3 r2 /d and λd+1 (Cov(ρ|Br (z) )) ≤
θ̃4 λd (Cov(ρ|Br (z) )), then the conditions in (A5) are satisfied by the Cj,k ’s we
construct with θ3 := θ̃3 (4C1 )−2 12−d and θ4 := θ̃4 /θ̃3 122d+2 C14 .
5b. if ρ is the volume measure on a closed C s Riemannian manifold isometrically
embedded in RD , then the conditions in (A5) are satisfied by the Cj,k ’s when j
is sufficiently large.
Even though the {Cj,k } does not exactly satisfy Assumption (A2), we claim that (15) is
sufficient for our performance guarantees in the case that M is bounded by M and d ≥ 3,
f by 0 gives the error:
since simply approximating points on M \ M
P

(ˆ

28M 2 log n
kxk dρ ≥
3n
f
M\M
2

)

≤ 2n−ν .

(16)

The constants in Proposition 13 are extremely pessimistic, due to the generality of
the assumptions on the space M. Indeed when M is a nice manifold as in case (5b),
the statement in the Proposition says that the constants for the Cj,k ’s we construct are
similar to those of the ideal Cj,k ’s. In practice we use a much simpler and more efficient
tree construction method and we experimentally obtain the properties above with amax =
C12 4d and amin = 1, at least for the vast majority of the points, and θ{3,4} u θ̃{3,4} . We
describe this simpler construction for the multiscale partitions in Appendix A.3, together
with experiments suggesting that θ{3,4} u θ̃{3,4} .
Besides cover tree, there are other methods that can be used in practice for the multiscale
partition, such as METIS by Karypis and Kumar (1999) that is used in Allard et al. (2012),
iterated PCA (see some analysis in Szlam (2009)) or iterated k-means. These can be
computationally more efficient than cover trees, with the downside being that they may
lead to partitions not satisfying our usual assumptions.
18

Adaptive Geometric Multiscale Approximations for Intrinsically Low-dimensional Data

3. Numerical experiments
We conduct numerical experiments on both synthetic and real data to demonstrate the
performance of our algorithms. Given {xi }ni=1 , we split them to training data for the
constructions of empirical GMRA and adaptive GMRA and test data for the evaluation of
the approximation errors:

Absolute error
Relative error




1

1
ntest

ntest

L2 error
1
X
2
b i k2
kxi − Px

xi ∈test set

X

xi ∈test set

b i k2 /kxi k2
kxi − Px

L∞ error
max

1

xi ∈test set

2

max

xi ∈test set

b ik
kxi − Px

b i k/kxi k
kxi − Px

b is the piecewise linear projection given by
where ntest is the cardinality of the test set and P
empirical GMRA or adaptive GMRA. In our experiments we use absolute error for synthetic
data, 3D shape and relative error for the MNIST digit data, natural image patches.
3.1 Synthetic data
We take samples {xi }ni=1 on the d-dim S and Z manifold whose x1 , x2 coordinates are
on the S and Z curve and xi ∈ [0, 1], i = 3, 4, 5 and evenly split them to the training
set and the test set. In the noisy case, training data are corrupted by Gaussian noise:
x̃train
= xtrain
+ √σD ξi , i = 1, . . . , n2 where ξi ∼ N (0, ID×D ), but test data are noise-free.
i
i
Test data error below the noise level imply that we are denoising the data.
3.1.1 Regularity parameter s in the As and Bs model

We sample 105 training points on the d-dim S or Z manifold (d = 3, 4, 5). The measure
on the S manifold is not exactly the volume measure but is comparable with the volume
measure.
The log-log plot of the approximation error versus scale in Figure 5 (b) shows that
volume measures on the d-dim S manifold are in As with s ≈ 2.0239, 2.1372, 2.173 when
d = 3, 4, 5, consistent with our theory which gives s = 2. Figure 5 (c) shows that volume
measures on the d-dim Z manifold are in As with s ≈ 1.5367, 1.6595, 1.6204 when d = 3, 4, 5,
consistent with our theory which gives s = 1.5.
The log-log plot of the approximation error versus the weighted complexity of the adaptive partition in Figure 5 (d) and (e) gives rises to an approximation of the regularity
parameter s in the Bs model in the table.
3.1.2 Error versus sample size n
We take n samples on the 4-dim S and Z manifold. In Figure 6, we set the noise level σ = 0
(a) and σ = 0.05 (b), display the log-log plot of the average approximation error over 5
trails with respect to the sample size n for empirical GMRA at scale j ∗ which is chosen
1
∗
as per Theorem 4: 2−j = [(log n)/n] 2s+d−2 with d = 4 and s = 2 for the S manifold and
s = 1.5 for the Z manifold. For adaptive GMRA, the ideal κ increases as σ increases. We
19

Liao and Maggioni

log 10(L2 error)

-1.5

-2

-2.5

log1/δ (absolute error)
c

0

log1/δ (absolute error)

(a) Projection of S
and Z manifolds
c

0

-1

-0.5

-2

-1.5

-3

-2.5

-4
-5

-3.5

c

log
1/ δ
c

15
20

the S manifold
theoretical s numerical s
2
2.1 ± 0.2
2
2.1 ± 0.2
2
2.3 ± 0.1
the Z manifold
theoretical s numerical s
+∞
3.4 ± 0.3
3
3.1 ± 0.2
2.25
2.5 ± 0.2
d=3: num.= -2.0239 theor.= -2
d=4: num.= -2.1372 theor.= -2
d=5: num.= -2.173 theor.= -2

10

d=3
d=4
d=5

(error) versus scale

5

scale

(b) Error versus scale

d=3: num.= -1.5367 theor.= -1.5
d=4: num.= -1.6595 theor.= -1.5
d=5: num.= -1.6204 theor.= -1.5

0

log 1/ δ (error) versus scale

0

d=3
d=4
d=5

0

-1

-1.5

-2

-2.5

-3

0.2

0.4

0.6

0.8

1

1.2

1.4

-3.5
-0.4

-0.2

-1

-2

-3

-4

-5

(c) Z curve

-6
0.5

1

0

0.2

numeric scale (-log 10radii)

0.4

0.6

0.8

1

1.2

1.4

numeric scale (-log 10radii)

(b) S: error vs. scale

log 10(errord-2 )

-1

-0.5

-2

-1.5

-3

-2.5

-3.5

-3

d=3: slope= -1.5367 theory= -1.5
d=4: slope= -1.6595 theory= -1.5
d=5: slope= -1.6204 theory= -1.5

-0.5

(c) Z: error vs. scale

log 10(error d-2 ) versus log10(complexity)

log 10(error d-2 ) versus log10(complexity)

0

d=3 slope= -3.6798 theory= -Inf
d=4 slope= -3.2833 theory= -3
d=5 slope= -2.6695 theory= -2.25

d=3 slope= -2.4648 theory= -2
d=4 slope= -2.2121 theory= -2
d=5 slope= -2.3197 theory= -2

-1

-2

log 10(errord-2 )

(a) S curve

-1

log 10(error) versus numeric scale (-log10radii)

0

d=3: slope= -2.0239 theory= -2
d=4: slope= -2.1372 theory= -2
d=5: slope= -2.173 theory= -2

log 10(L2 error)

-0.5

-3.5
-0.2

Robust, Adaptive Geometric Multi-Resolution Analysis for data in high-dimensions

log 10(error) versus numeric scale (-log10radii)

0

-3

-4

-5

1.5

2

2.5

3

-6
0.6

0.8

1

1.2

1.4

1.6

1.8

2

2.2

2.4

2.6

log10(complexity)

log10(complexity)

(d) S: error vs. tree complexity

(e) Z: error vs. tree complexity

Figure 5: 105 training points are sampled on the d-dimensional S or Z manifold (d = 3, 4, 5).
bj Xkn , versus scale j. The negative of the
In (b) and (c), we display log10 kX − P
slope on the solid portion of the line approximates the regularity parameter s in
versus
the As model. In (d) and (e), we display the log-log plot of kX −PΛb η Xkd−2
n
the weighted complexity of the adaptive partition for the d-dimensional S and Z
manifold. The negative of the slope on the solid portion of the line approximates
the regularity parameter s in the Bs model. Our five experiments give the s in the
table. For the 3-dim Z manifold, while s = +∞ in the case of infinite samples,
we do obtain a large s with 105 samples.

let κ ∈ {0.05, 0.1} when σ = 0 and κ ∈ {1, 2} when σ = 0.05. We also test the Nearest
Neighbor (NN) approximation. The negative of the slope, determined by least squared fit,
gives rise to the rate of convergence: L2 error ∼ (1/n)−slope . When σ = 0, the convergence
rate for the nearest neighbor approximation should be 1/d = 0.25. GMRA gives rise to a
smaller error and a faster rate of convergence than the nearest neighbor approximation. For
the Z manifold, Adaptive GMRA yields a faster rate of convergence than GMRA. When
σ = 0.05, adaptive GMRA with κ = 0.5 and 1 gives rise to the fastest rate of convergence.
Adaptive GMRA with κ = 0.05 has similar rate of convergence as the nearest neighbor
approximation since the tree is almost truncated at the finest scales. We note a de-noising
20

Adaptive Geometric Multiscale Approximations for Intrinsically Low-dimensional Data

log 10(Error) versus log10(Sample size)

-0.5

log 10(Error) versus log10(Sample size)

-0.5

Nearest neighbor: slope= -0.25743 theory= -0.25
jstar: slope= -0.38219 theory= -0.33333
Adaptive κ=0.05: slope= -0.42188 theory= -0.33333
Adaptive κ=0.1: slope= -0.31197 theory= -0.33333

Nearest neighbor: slope= -0.17803
jstar: slope= -0.18289
Adaptive κ=1: slope= -0.25691
Adaptive κ=2: slope= -0.23912
log 10(σ)

-1

log 10(Absolute error)

log 10(Absolute error)

-1

-1.5

-2

-1.5

-2
-2.5

-3

-2.5
4

4.2

4.4

4.6

4.8

5

5.2

5.4

5.6

5.8

6

4

4.2

4.4

4.6

log 10(Sample size)

5

5.2

5.4

5.6

5.8

6

log 10(Sample size)

(a) the 4-dim S manifold, σ = 0

(b) the 4-dim S manifold, σ = 0.05

log 10(Error) versus log10(Sample size)

-0.5

4.8

log 10(Error) versus log10(Sample size)

-0.5

Nearest neighbor: slope= -0.25631 theory= -0.25
jstar: slope= -0.35185 theory= -0.3
Adaptive κ=0.05: slope= -0.39452 theory= -0.375
Adaptive κ=0.1: slope= -0.39127 theory= -0.375

Nearest neighbor: slope= -0.1858
jstar: slope= -0.1362
Adaptive κ=1: slope= -0.37071
Adaptive κ=2: slope= -0.33902
log 10(σ)

-1

log 10(Absolute error)

log 10(Absolute error)

-1

-1.5

-2

-1.5

-2
-2.5

-3

-2.5
4

4.2

4.4

4.6

4.8

5

5.2

5.4

5.6

5.8

6

4

4.2

4.4

log 10(Sample size)

4.6

4.8

5

5.2

5.4

5.6

5.8

6

log 10(Sample size)

(c) the 4-dim Z manifold, σ = 0

(d) the 4-dim Z manifold, σ = 0.05

Figure 6: L2 error versus the sample size n, for the 4-dim S and Z manifolds (d = 4), of
GMRA at the scale j ∗ chosen as per Theorem 4 and adaptive GMRA with varied
κ. We let κ ∈ {0.05, 0.1} when σ = 0 and κ ∈ {1, 2} when σ = 0.05.

effect when the approximation error falls below σ as n increases. In adaptive GMRA, when
κ is sufficiently large, i.e., κ = 0.5, 1 in this example, different values of κ do yield different
errors up to a constant, but the rate of convergence is independent of κ, as predicted by
Theorem 8.
3.1.3 Robustness of GMRA and adaptive GMRA
The robustness of the empirical GMRA and adaptive GMRA is tested on the 4-dim S and
Z manifold while σ varies but n is fixed to be 105 . Figure 3.1.3 shows that the average L2
approximation error in 5 trails increases linearly with respect to σ for both uniform and
adaptive GMRA with κ ∈ {0.05, 0.5, 1}.
21

Liao and Maggioni

0.4
S: uniform at j ⋆
S: adaptive κ =0.05
S: adaptive κ =0.5
S: adaptive κ =1
Z: uniform at j ⋆
Z: adaptive κ = 0.05
Z: adaptive κ = 0.5
Z: adaptive κ = 1
σ √
σ/ D

0.35

0.3

0.25

L2 error

Figure 7: The average L2 approximation error in 5
trails versus σ for GMRA and adaptive
GMRA with κ ∈ {0.05, 0.5, 1} on data sampled on the 4-dim S and Z manifolds. This
shows the error of approximation grows linearly with the noise size, suggesting robustness in the construction.

0.2

0.15

0.1

0.05

0
0

0.05

0.1

0.15

0.2

0.25

σ

3.2 3D shapes

(a) 41,472 points

(b) 165,954 points

(c) 437,645 points

(d) κ ≈ 0.18, partition size = 338 (e) κ ≈ 0.41, partition size = 749 (f) κ ≈ 0.63, partition size = 1141

Figure 8: Top line: 3D shapes;p
bottom line: adaptive partitions selected with refinement
b j,k ≥ 2−j κ (log n)/n. Every cell is colored by scale. In the adaptive
criterion ∆
partition, at irregular locations cells are selected at finer scales than at “flat”
locations.

We run GMRA and adaptive GMRA on 3D points clouds on the teapot, armadillo and
dragon in Figure 8. The teapot data are from the matlab toolbox and others are from the
Stanford 3D Scanning Repository http://graphics.stanford.edu/data/3Dscanrep/.
22

Adaptive Geometric Multiscale Approximations for Intrinsically Low-dimensional Data

Figure 8 shows that the adaptive partitions chosen by adaptive GMRA matches our
expectation that, at irregular locations, cells are selected at finer scales than at “flat”
locations.
In Figure 9, we display the absolute L2 /L∞ approximation error on test data versus
scale and partition size. The left column shows the L2 approximation error versus scale
for GMRA and the center approximation. While the GMRA approximation is piecewise
linear, the center approximation is piecewise constant. Both approximation errors decay
from coarse to fine scales, but GMRA yields a smaller error than the approximation by local
centers. In the middle column, we run GMRA and adaptive GMRA with the L2 refinement
criterion defined in Table 2 with scale-dependent (∆j,k ≥ 2−j τn ) and scale-independent
(∆j,k ≥ τn ) threshold respectively, and display the log-log plot of the L2 approximation error
versus the partition size. Overall adaptive GMRA yields the same L2 approximation error
as GMRA with a smaller partition size, but the difference is insignificant in the armadillo
and dragon, as these 3D shapes are complicated and the L2 error simply averages the error
at all locations. Then we implement adaptive GMRA with the L∞ refinement criterion:
bj+1 xi − P
bj xi k and display the log-log plot of the L∞ approximation
b ∞ = maxx ∈C kP
∆
i
j,k
j,k
error versus the partition size in the right column. In the L∞ error, adaptive GMRA saves
a considerable number (about half) of cells in order to achieve the same approximation
error as GMRA. In this experiment, scale-independent threshold is slightly better than
scale-dependent threshold in terms of saving the partition size.
3.3 MNIST digit data
We consider the MNIST data set from http://yann.lecun.com/exdb/mnist/, which contains images of 60, 000 handwritten digits, each of size 28 × 28, grayscale. The intrinsic
dimension of this data set varies for different digits and across scales, as it was observed in
Little et al. (2012). We run GMRA by setting the diameter of cells at scale j to be O(0.9j )
in order to slowly zoom into the data at multiple scales.
We evenly split the digits to the training set and the test set. As the intrinsic dimension
is not well-defined, we set GMRA to pick the dimension of Vbj,k adaptively, as the smallest
dimension needed to capture 50% of the energy of the data in Cj,k . As an example, we
display the GMRA approximations of the digit 0, 1, 2 from coarse scales to fine scales in
Figure 10. The histogram of the dimensions of the subspaces Vbj,k is displayed in (a). (b)
bj+1 xi − P
bj xi k from the coarsest scale (top) to the finest scale (bottom),
represents log10 kP
with columns indexed by the digits, sorted from 0 to 9. We observe that 1 has more fine
scale information than the other digits. In (c), we display the log-log plot of the relative
L2 error versus scale in GMRA and the center approximation. The improvement of GMRA
over center approximation is noticeable. Then we compute the relative L2 error for GMRA
and adaptive GMRA when the partition size varies. Figure 10 (d) shows that adaptive
GMRA achieves the same accuracy as GMRA with fewer cells in the partition. Errors
increase when the partition size exceeds 103 due to a large variance at fine scales. In
this experiment, scale-dependent threshold and scale-independent threshold yield similar
performances.
23

Liao and Maggioni

Teapot
log 10(error) versus scale (-log10radii)

0.5

log 10(error) versus log10(partition size)

0
Center: slope= -0.97936
GMRA: slope= -1.204

log 10(error) versus log10(partition size)

0.4

Uniform:
slope= -0.63399
Adaptive scaled: slope= -0.71067
Adaptive unscaled: slope= -0.67156

-0.2

Uniform:
slope= -0.49699
Adaptive scaled: slope= -0.51588
Adaptive unscaled: slope= -0.51359

0.2
0

-0.4

-0.6

0

log 10(error)

log 10(error)

log 10(error)

-0.5
-0.8

-1

-0.2

-1
-1.2

-0.4

-1.4

-1.5

-0.6
-1.6

-2
-0.6

-1.8
-0.4

-0.2

0

0.2

0.4

0.6

0.8

-0.8
0

0.5

1

1.5

-log 10(radii)

2

2.5

3

0

0.5

1

1.5

log 10(partition size)

2

2.5

3

log 10(partition size)

Armadillo
log 10(error) versus scale (-log10radii)

2

log 10(error) versus log10(partition size)

1.5
Center: slope= -0.94768
GMRA: slope= -1.2524

log 10(error) versus log10(partition size)

1.8

Uniform:
slope= -0.63738
Adaptive scaled: slope= -0.64301
Adaptive unscaled: slope= -0.63819

Uniform:
slope= -0.46372
Adaptive scaled: slope= -0.46232
Adaptive unscaled: slope= -0.45022

1.6

1.5
1
1.4
1
1.2

log 10(error)

log 10(error)

log 10(error)

0.5
0.5

0

1

0.8

0
0.6
-0.5
-0.5
0.4

-1
-2.2

-1
-2

-1.8

-1.6

-1.4

-1.2

-1

-0.8

-0.6

-0.4

-0.2

0.2
0

0.5

1

1.5

-log 10(radii)

2

2.5

3

3.5

4

0

0.5

1

1.5

log 10(partition size)

2

2.5

3

3.5

4

log 10(partition size)

Dragon
log 10(error) versus scale (-log10radii)

-1

log 10(error) versus log10(partition size)

-1.5
Center: slope= -0.9771
GMRA: slope= -1.1887

log 10(error) versus log10(partition size)

-1.2

Uniform:
slope= -0.54577
Adaptive scaled: slope= -0.6199
Adaptive unscaled: slope= -0.59983

Uniform:
slope= -0.34604
Adaptive scaled: slope= -0.40454
Adaptive unscaled: slope= -0.41434

-1.4

-1.5
-2

-1.6
-2

-2.5

log 10(error)

log 10(error)

log 10(error)

-2.5
-1.8

-2

-3
-3
-2.2
-3.5
-3.5

-4
0.8

-2.4

-4
1

1.2

1.4

1.6

1.8

2

2.2

2.4

2.6

2.8

-log 10(radii)

-2.6
0

0.5

1

1.5

2

log 10(partition size)

2.5

3

3.5

4

0

0.5

1

1.5

2

2.5

3

3.5

4

log 10(partition size)

Figure 9: Left column: log10 (L2 error) versus scalefor GMRA and center approximation;
Middle column: log-log plot of the L2 error versus partition size for GMRA and
adaptive GMRA with scale-dependent and scale-independent threshold under the
L2 refinement defined in Table 2; Right column: log-log plot of L∞ error versus
partition size for GMRA and adaptive GMRA with scale-dependent and scaleindependent threshold under the L∞ refinement.

3.4 Natural image patches
It was argued in Peyré (2009) that many sets of patches extracted from natural images
can be modeled a low-dimensional manifold. We use the Caltech 101 dataset from https:
//www.vision.caltech.edu/Image_Datasets/Caltech101/ (see F. Li and Perona, 2006),
take 40 images from four categories: accordion, airplanes, hedgehog and scissors and extract
multiscale patches of size 8 × 8 from these images. Specifically, if the image is of size m × m,
for ` = 1, . . . , log2 (m/8), we collect patches of size 2` 8, low-pass filter them and downsample
them to become patches of size 8 × 8 (see Gerber and Maggioni (2013) for a discussion
24

Adaptive Geometric Multiscale Approximations for Intrinsically Low-dimensional Data

1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

19

20

21

1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

19

20

21

1

2

3

4

5

6

7

8

9

10

11

12

(a) Scaling subspace

14

15

16

17

18

21

Uniform
Adaptive scaled
Adaptive unscaled

-0.1

-0.26

-0.15

-0.28

log 10(relative error)

log 10(error)

20

log 10(error) versus log10(partition size)

-0.24
Center
GMRA

-0.2

-0.3

-0.25

-0.32

-0.3

-0.34

-0.35
-3.5

19

bj+1 xi − P
bj xi k
(b) log10 kP

log 10(error) versus scale (-log10radii)

-0.05

13

-0.36
-3.45

-3.4

-3.35

-3.3

-3.25

-3.2

-3.15

-3.1

-3.05

0

-log 10(radii)

0.5

1

1.5

2

2.5

3

3.5

4

log 10(partition size)

(c) Error versus scale

(d) Error versus partition size

Figure 10: The top three rows: multiscale approximations of the digit 0, 1, 2 in the MNIST
data set, from the coarsest scale (left) to the finest scale (right). (a) the hisbj+1 xi − P
bj xi k from the
togram of dimensions of the subspaces Vbj,k ; (b) log10 kP
coarsest scale (top) to the finest scale (bottom), with columns indexed by the
digits, sorted from 0 to 9; (c) log-log plot of the relative L2 error versus scale
in GMRA and the center approximation; (d) log-log plot of the relative L2 error versus partition size for GMRA, adaptive GMRA with scale-dependent and
scale-independent threshold.

25

Liao and Maggioni

Learning image patches
log 10(error) versus scale (-log10radii)

0

log 10(error) versus log10(partition size)

-0.15
Center
GMRA

Uniform
Adaptive scaled
Adaptive unscaled

-0.16

-0.05

-0.17

log 10(relative error)

-0.18

log 10(error)

-0.1

-0.15

-0.19
-0.2
-0.21
-0.22

-0.2

-0.23
-0.24

-0.25
-0.1

-0.25
-0.05

0

0.05

0.1

0.15

0.2

0

0.5

1

1.5

-log 10(radii)

(a) Scaling subspace

2

2.5

3

3.5

4

4.5

log 10(partition size)

(b) Error versus scale

(c) Error versus partition size

Learning the Fourier magnitudes of image patches
log 10(error) versus scale (-log10radii)

-0.2

log 10(error) versus log10(partition size)

-0.4
Center
GMRA

-0.25

Uniform
Adaptive scaled
Adaptive unscaled

-0.45

-0.3

log 10(relative error)

log 10(error)

-0.35

-0.4

-0.45

-0.5

-0.55

-0.5

-0.55

-0.6

-0.6

-0.65

-0.65
-1

-0.9

-0.8

-0.7

-0.6

-0.5

-0.4

-0.3

0

-log 10(radii)

(d) Scaling subspace

(e) Error versus scale

0.5

1

1.5

2

2.5

3

3.5

4

4.5

log 10(partition size)

(f) Error versus partition size

Figure 11: Caltech 101 image patches
Figure 12: Top line: learning on 200, 000 image patches; bottom line: results of learning
the Fourier magnitudes of the same image patches. (a,d) histograms of the
dimensions of the subspaces Vbj,k ; (b,e) relative L2 error versus scale for GMRA
and the center approximation; (c,f) relative L2 error versus the partition size for
GMRA, adaptive GMRA with scale-dependent and scale-independent threshold.

about dictionary learning on patches of multiple sizes using multiscale ideas). Then we
randomly pick 200, 000 patches, evenly split them to the training set and the test set. In
the construction of GMRA, we set the diameter of cells at scale j to be O(0.9j ) and the
dimension of Vbj,k to be the smallest dimension needed to capture 50% of the energy of the
data in Cj,k . We also run GMRA and adaptive GMRA on the Fourier magnitudes of these
image patches to take advantage of translation-invariance of the Fourier magnitudes. The
results are shown in Figure 12. The histograms of the dimensions of the subspaces Vbj,k
are displayed in (a,e). Figure 12 (c) and (g) show the relative L2 error versus scale for
GMRA and the center approximation. We then compute the relative L2 error for GMRA
and adaptive GMRA when the partition size varies and display the log-log plot in (d) and
(h). It is noticeable that adaptive GMRA achieves the same accuracy as GMRA with a
smaller partition size. We conducted similar experiments on 200, 000 multiscale patches
from CIFAR 10 from https://www.cs.toronto.edu/~kriz/cifar.html (see Krizhevsky
and Hinton, 2009) with extremely similar results (not shown).
26

Adaptive Geometric Multiscale Approximations for Intrinsically Low-dimensional Data

4. Performance analysis of GMRA and adaptive GMRA
This section is devoted to the performance analysis of empirical GMRA and adaptive
GMRA. We will start with the following stochastic error estimate on any partition.
4.1 Stochastic error on a fixed partition
Suppose T̃ is a finite proper subtree of the data master tree T n . Let Λ be the partition
consisting the outer leaves of T̃ . The piecewise affine projector on Λ and its empirical
version are
X
X
bj,k 1j,k .
bΛ =
P
PΛ =
Pj,k 1j,k and P
Cj,k ∈Λ

Cj,k ∈Λ

bΛ Xk is given by:
A non-asymptotic probability bound on the stochastic error kPΛ X − P

Lemma 14 Let Λ be the partition associated a finite proper subtree T̃ of the data master
tree T n . Suppose Λ contains #j Λ cells at scale j. Then for any η > 0,
−

bΛ Xk ≥ η} ≤ αd · #Λ · e
P{kPΛ X − P
bΛ Xk2 ≤
EkPΛ X − P

d2 log(αd#Λ)

d2

βnη 2
P −2j
#j Λ
j2

P

j

(17)

2−2j #j Λ

βn

where α = α(θ2 , θ3 ) and β = β(θ2 , θ3 , θ4 ).

Lemma 14 and Proposition 15 below are proved in appendix C .
4.2 Performance analysis of empirical GMRA
According to Eq. (1), the approximation error of empirical GMRA is split into the squared
bias and the variance. A corollary of Lemma 14 with Λ = Λj results in an estimate of the
variance term.
Proposition 15 For any η ≥ 0,
− β2

2j nη 2

bj Xk ≥ η} ≤ αd#Λj e d2 #Λj
P{kPj X − P
2
bj Xk2 ≤ d #Λj log[αd#Λj ] .
EkPj X − P
β22j n

(18)
(19)

In Eq. (1), the squared bias decays like O(2−2js ) whenever ρ ∈ As and the variance
scales like O(j2j(d−2) /n). A proper choice of the scale j gives rise to Theorem 4 whose proof
is given below.
Proof of Theorem 4
Proof [Proof of Theorem 4]
bj Xk2
bj Xk2 ≤ kX − Pj Xk2 + EkPj X − P
EkX − P
≤ |ρ|2As 2−2sj +

d2 #Λj log[αd#Λj ]
d2 2j(d−2)
αd2jd
2
−2sj
≤
|ρ|
2
+
log
A
s
β22j n
θ1 βn
θ1
27

Liao and Maggioni

as #Λj ≤ 2jd /θ1 due to Assumption (A3).
Intrinsic dimension d = 1: In this case, both the squared bias and the variance decrease
as j increases, so we should choose the scale j ∗ as large as possible as long as most cells at
∗
scale j ∗ have d points. We will choose j ∗ such that 2−j = µ logn n for some µ > 0. After
log n
grouping Λj ∗ into light and heavy cells whose measure is below or above 28(ν+1)
, we can
3n
log n 2
show that the error on light cells is upper bounded by C( n ) and all heavy cells have at
least d points with high probability.
∗

Lemma 16 Suppose j ∗ is chosen such that 2−j = µ logn n with some µ > 0. Then


28(ν + 1)θ22 µ log n 2
,
k(X − P X)1{C ∗ :ρ(C ∗ )≤ 28(ν+1) log n } k ≤
j ,k
j ,k
3θ1
n
3n
n
o
log n
P each Cj ∗ ,k satisfying ρ(Cj ∗ ,k ) > 28(ν+1)
has at least d points ≥ 1 − n−ν .
3n
2

j∗

Lemma 16 is proved in appendix D. If j ∗ is chosen as above, The probability estimate in
(5) follows from


log n
log n s
∗
−sj
s
≤ |ρ|As µs
,
kX − Pj ∗ Xk ≤ |ρ|As 2
≤ |ρ|As µ
n
n


bj ∗ Xk ≥ C1 log n
P kPj ∗ X − P
n



≤

α
θ1 µ

2 log n
µβθ1 C1
−
d2

((log n)/n)−1 e

2
µβθ1 C1

n n− d2
≤ θα1 µ
log n

≤ C2 n−ν

provided that µβθ1 C12 /d2 − 1 > ν.
Intrinsic dimension d ≥ 2: When d ≥ 2, the squared bias decreases but the variance
1
∗
increases as j gets large. We choose j ∗ such that 2−j = µ ((log n)/n) 2s+d−2 to balance
these two terms. We use the same technique as d = 1 to group Λj ∗ into light and heavy
cells whose measure is below or above 28/3 · (ν + 1)(log n)/n, we can show that the error on
2s
light cells is upper bounded by C((log n)/n) 2s+d−2 and all heavy cells have at least d points
with high probability.
1

∗

Lemma 17 Let j ∗ be chosen such that 2−j = µ ((log n)/n) 2s+d−2 with some µ > 0. Then

 2s
28(ν + 1)θ22 µ2−d log n 2s+d−2
k(X − P X)1{C ∗ :ρ(C ∗ )≤ 28(ν+1) log n } k ≤
,
j ,k
j ,k
3θ1
n
3n


28(ν + 1) log n
P ∀ Cj ∗ ,k : ρ(Cj ∗ ,k ) >
, Cj ∗ ,k has at least d points ≥ 1 − n−ν .
3n
2

j∗

Proof of Lemma 17 is omitted since it is the same as the proof of Lemma 16. The probability
estimate in (6) follows from
kX − Pj ∗ Xk ≤ |ρ|As 2

−sj ∗

s

≤ |ρ|As µ

28



log n
n



s
2s+d−2

,

Adaptive Geometric Multiscale Approximations for Intrinsically Low-dimensional Data

(

bj ∗ Xk ≥ C1
P kPj ∗ X − P



log n
n



s
2s+d−2

)

αdµ−d
≤
θ1



log n
n

−

d
2s+d−2

e−

2 µd−2 log n
βθ1 C1
d2

≤ C2 n−ν

provided that βθ1 C12 µd−2 /d2 − 1 > ν.

4.3 Performance analysis of Adaptive GMRA
Proof [Proof of Theorem 8] In the case that M is bounded by M , the minimum scale
θ2
jmin = log2 M
. We first consider the case d ≥ 3. In our proof C stands for constants that
may vary at different locations, but it is independent of n and D. We will begin by defining
several objects of interest:
• T n : the data master tree whose leaf contains at least d points in Xn . It can be viewed
as the part of a multiscale tree that our data have explored.
• T : a complete multiscale tree containing T n . T can be viewed as the union T n and
some empty cells, mostly at fine scales with high probability, that our data have not
explored.
• T(ρ,η) : the smallest subtree of T which contains {Cj,k ∈ T : ∆j,k ≥ 2−j η}.
• Tη = T(ρ,η) ∩ T n .
b j,k ≥ 2−j η}.
• Tbη : the smallest subtree of T n which contains {Cj,k ∈ T n : ∆
• Λ(ρ,η) : the partition associated with T(ρ,η) .
• Λη : the partition associated with Tη .
b η : the partition associated with Tbη .
• Λ

• Suppose T 0 and T 1 are two subtrees of T . If Λ0 and Λ1 are two adaptive partitions
associated with T 0 and T 1 respectively, we denote by Λ0 ∨Λ1 and Λ0 ∧Λ1 the partitions
associated to the trees T 0 ∪ T 1 and T 0 ∩ T 1 respectively.
We also let b = 2amax + 5 where amax is the maximal number of children that a node has
in T ; κ0 = max(κ1 , κ2 ) where b2 κ21 /(21θ22 ) = ν + 1 and α2 κ22 /b2 = ν + 1 with α2 defined in
Lemma 19. In order the obtain the MSE bound, one can simply set
ν = 1.
bb = P
b
The empirical adaptive GMRA projection is given by P
b τ Pj,k 1j,k . Using
Λτn
Cj,k ∈Λ
n
the triangle inequality, we split the error as follows:
bb Xk ≤ e1 + e2 + e3 + e4
kX − P
Λτ
n

where
e1 := kX − PΛb τ
e3 := kPΛb τ

n ∨Λbτn

n ∧Λτn /b

e2 := kPΛb τ

Xk ,

bb
X −P
Λτ

n ∨Λbτn

n ∧Λτn /b

bb
e4 := kP
Λτ

Xk ,

X − PΛb τ

n ∧Λτn /b

29

n ∧Λτn /b

bb Xk.
X −P
Λτ
n

Xk

Liao and Maggioni

A similar split appears in the works of Binev et al. (2005, 2007). The partition built from
b j,k ≥ 2−j τn does not exactly coincide with the partition chosen
those Cj,k ’s satisfying ∆
based on those Cj,k satisfying ∆j,k ≥ 2−j τn . This is accounted by e2 and e4 , corresponding
b j,k is significantly larger or smaller than ∆j,k , which we will prove
to those Cj,k ’s whose ∆
to be small with high probability. The remaining terms e1 and e3 correspond to the bias
and variance of the approximations on the partition obtained by thresholding ∆j,k .
b τn ∨ Λbτ ⊇ Λbτ ,
Term e1 : The first term e1 is essentially the bias term. Since Λ
n
n
e21 = kX − PΛb τ

n ∨Λbτn

Xk2 ≤ kX − PΛbτn Xk2 ≤ kX − PΛ(ρ,bτn ) Xk2 + kPΛ(ρ,bτn ) X − PΛbτn Xk2 .
|
{z
} |
{z
}
e211

e212

e211 may be upper bounded deterministically from Eq. (9):
e211

≤

Bs,d |ρ|pBs (bτn )2−p

2(d−2)
2s+d−2

≤ Bs,d |ρ|Bs

(bκ)

4s
2s+d−2



log n
n



2s
2s+d−2

.

(20)

e12 encodes the difference between thresholding T and T n , but it is 0 with high probability:
Lemma 18 For any ν > 0, κ such that κ > κ1 , where b2 κ21 /(21θ22 ) = ν + 1,
P{e12 > 0} ≤ C(θ2 , amax , amin , κ)n−ν

(21)

The proof is postponed, together with those of the Lemmata that follow, to appendix D).
If M is bounded by M , then e212 ≤ 4M 2 and

 2s
log n 2s+d−2
2
2
2
−ν
2
Ee12 ≤ 4M P{e12 > 0} ≤ 4M Cn ≤ 4M C
(22)
n
if ν > 2s/(2s + d − 2), for example ν = 1.
b τn ∧ Λτ /b . For any η > 0,
Term e3 : e3 corresponds to the variance on the partition Λ
n
b τn ∧ Λτ /b )e
P{e3 > η} ≤ αd#(Λ
n

−

βnη 2
P
b τ ∧Λ
d2 j≥j
2−2j #j (Λ
n
τn /b )
min

b τn ∧ Λτ /b ⊂ Tτ /b , for any j ≥ 0, regardless of Λ
b τn , we have
according to Lemma 14. Since Λ
n
n
b τn ∧ Λτ /b ) ≤ #j Tτ /b ≤ #Tτ /b . Therefore
#j (Λ
n
n
n
P{e3 > η} ≤ αd#Tτn /b e

which implies
ˆ
ˆ +∞
Ee23 =
ηP {e3 > η} dη =
0

0

−

βnη 2
P
d2 j≥j
2−2j #j Tτ /b
n
min

+∞

≤ αd#Tτn /b e
−

η min 1, αd#Tτn /b e

−

βnη 2
p
d2 |ρ|
(τ /b)−p
Bs n

βnη 2
P
d2 j≥j
2−2j #j Tτ /b
n
min

!

,

(23)

dη



d2 log αd#Tτn /b X −2j
log n  τn −p
log n 2s+d−2
≤ C(θ2 , θ3 , d, κ, s, |ρ|Bs )
≤
2 #j Tτn /b ≤ C
.
βn
n
b
n
2s

j≥jmin

Term e2 and e4 : These terms account for the difference of truncating the master tree
b j,k ’s. We prove that ∆
b j,k ’s concentrate
based on ∆j,k ’s and its empirical counterparts ∆
near ∆j,k ’s with high probability if there are sufficient samples.
30

Adaptive Geometric Multiscale Approximations for Intrinsically Low-dimensional Data

Lemma 19 For any η > 0 and any Cj,k ∈ T
n n
o n
oo
2j
2
b j,k ≤ η and ∆j,k ≥ bη , P ∆j,k ≤ η and ∆
b j,k ≥ bη
max P ∆
≤ α1 e−α2 2 nη (24)
for some constants α1 := α1 (θ2 , θ3 , amax , d) and α2 := α2 (θ2 , θ3 , θ4 , amax , d).

This Lemma enables one to show that e2 = 0 and e4 = 0 with high probability:
Lemma 20 Let α1 and α2 be the constants in Lemma 19. For any fixed ν > 0,
P{e2 > 0} + P{e4 > 0} ≤ α1 amin n−ν

(25)

when κ is chosen such that κ > κ2 , with α2 κ22 /b2 = ν + 1.
Since M is bounded by M , we have e22 ≤ 4M 2 so
Ee22

2

2

≤ 4M P{e2 > 0} ≤ 4M α1 amin n

−ν

2

≤ 4M α1 amin



log n
n



2s
2s+d−2

if ν > 2s/(2s + d − 2), for example ν = 1. The same bound holds for e4 .

2(d−2)

4s

(bκ) 2s+d−2 such
Finally, we complete the probability estimate (10): let c20 = Bs,d |ρ|B2s+d−2
s
s

that e11 ≤ c0 ((log n)/n) 2s+d−2 . We have
o
n
s
bb Xk ≥ c1 ((log n)/n) 2s+d−2
P kX − P
Λτn
n
o
s
≤ P e3 > (c1 − c0 ) ((log n)/n) 2s+d−2 + P{e12 > 0} + P{e2 > 0} + P{e4 > 0}
o
n
s
≤ P e3 > (c1 − c0 ) ((log n)/n) 2s+d−2 + Cn−ν ,

as long as κ is chosen such that κ > max(κ1 , κ2 ) where b2 κ21 /(21θ22 ) = ν + 1 and α2 κ22 /b2 =
ν + 1 according to (21) and (25). Applying (23) gives rise to
2s
n
o
− p βn −p (c1 −c0 )2 ((log n)/n) 2s+d−2
s
|ρ|
(τn /b)
P e3 > (c1 − c0 ) ((log n)/n) 2s+d−2 ≤ αd#Tτn /b e Bs

≤ αd#Tτn /b n

−

β(c1 −c0 )2 κp
p
bp |ρ|
Bs

≤ αdamin n

if c1 is taken large enough such that



β(c1 −c0 )2 κp
−
−1
p
p
b |ρ|
Bs

β(c1 −c0 )2 κp
bp |ρ|pBs

≤ αdamin n−ν

≥ ν + 1.

We are left with the cases d = 1, 2. When d = 1, for any distribution ρ satisfying
quasi-orthogonality (8) and any η > 0, the tree complexity may be bounded as follows:
X
X
2−2j #j T(ρ,η) ≤
2−2j 2j /θ1 = 2/θ1 2−jmin = 2M/(θ1 θ2 ) ,
j≥jmin

j≥jmin

so kX − PΛ(ρ,η) Xk2 ≤ 8M B0 η 2 /(3θ1 θ2 ). Hence
e211 ≤

2
8M B0
3θ1 θ2 (bτn )

≤

8M B0 b2 κ2
3θ1 θ2 (log n)/n,

31

P{e3 > η} ≤ αd#Tτn /b e−

θ1 θ2 βnη 2
2M d2

,

Liao and Maggioni

which yield Ee23 ≤ 2M d2 log αd#Tτn /b /(θ1 θ2 βn) ≤ C(log n)/n and estimate (11).
WhenPd = 2, for any distribution satisfying quasi-orthogonality and given any η > 0,
we have j≥jmin 2−2j #j T(ρ,η) ≤ −|ρ|−1 log η, whence kX − PΛ(ρ,η) Xk2 ≤ − 43 B0 |ρ|η 2 log η.
Therefore
e211 ≤ − 43 B0 |ρ|(bτn )2 log(bτn ) ≤ C(log2 n)/n

,

P{e3 > η} ≤ αd#Tτn /b e

−

2βnη 2
d2 |ρ| log n

,

which yield Ee23 ≤ d2 |ρ| log αd#Tτn /b (log n)/(2βn) ≤ C(log2 n)/n and the probability estimate (12).
Proof [Proof of Theorem 9] Let R > 0. If we run adaptive GMRA on BR (0), and approximate points outside BR (0) by 0, the MSE of the adaptive GMRA in BR (0) is
2s

2s

bb )1{kxk≤R} Xk2 . (|ρ|B (R) |p + R2 ) ((log n)/n) 2s+d−2 . Rmax(λ,2) ((log n)/n) 2s+d−2 .
k(I − P
0
Λτ
n

The squared error outside BR (0) is

ˆ
2

k1{kxk≥R} Xk =

BR (0)c

||x||2 dρ ≤ CR−δ .

(26)

The total MSE is
2s

MSE . Rmax(λ,2) ((log n)/n) 2s+d−2 + R−δ .
2s
− (2s+d−2)(δ+max(2,λ))

Minimizing over R suggests taking R = Rn = max(R0 , µ(log n/n)
2s
δ
·
2s+d−2 δ+max(λ,2)

ing MSE . ((log n)/n)
and Eq. (10) in Theorem 8.

), yield-

. The probability estimate (13) follows from Eq. (26)

In Remark 10, we claim that λ is not large in simple cases. If ρ ∈ A∞
s and ρ decays such that ρ(Cj,k ) ≤ 2−jd kcj,k k−(d+1+δ) , we have ∆j,k ≤ 2−js 2−jd/2 kcj,k k−(d+1+δ)/2 .
Roughly speaking, for any η > 0, the cells of distance r to 0 satisfying ∆j,k ≥ 2−j η
d+1+δ

2

will satisfy 2−j ≥ (ηr 2 ) 2s+d−2 . In other words, the cells of distance r to 0 are trun2
d+1+δ
cated at scale jmax such that 2−jmax = (ηr 2 ) 2s+d−2 , which gives rise to complexity
2(d−2)

(d+1+δ)(d−2)

≤ 2−2jmax rd−1 2jmax d ≤ η − 2s+d−2 rd−1− 2s+d−2 . If we run adaptive GMRA with threshold η on BR (0), the weighted complexity of the truncated tree is upper bounded by
2(d−2)

η − 2s+d−2 rd−
λ=d−

(d+1+δ)(d−2)
2s+d−2

(d+1+δ)(d−2)
.
2s+d−2

. Therefore, ρ|BR (0) ∈ Bs for all R > 0 and |ρ|BR (0) |pBs ≤ Rλ with

5. Discussions and extensions
5.1 Computational complexity
The computational costs of GMRA and adaptive GMRA are summarized in Table 3.
32

Adaptive Geometric Multiscale Approximations for Intrinsically Low-dimensional Data

Operations

Computational cost

Multiscale tree construction

C d Dn log n
−jd
Ddn2
| {z }

Randomized PCA at scale j

·

2jd
|{z}

= Ddn

PCA cost at one node number of nodes

Randomized PCA at all nodes

Ddn
|{z}

·

1/d log n
| {z }

= Dn log n

cost at a fixed scale number of scales

Computing ∆j,k ’s

Ddn
|{z}

·

1/d log n
| {z }

= Dn log n

cost for single scale number of scales

Compute Pj (x) for a new sample x

D log n
| {z }

find Cj,k containing x

+

Dd
|{z}

= D(log n + d)

compute Pj,k (x)

Table 3: Computational cost

5.2 Quasi-orthogonality
A main difference between GMRA and orthonormal wavelet bases (see Daubechies, 1992;
Mallat, 1998) is that Vj,x * Vj+1,x where (j, x) = (j, k) such that x ∈ Cj,k . Therefore the
geometric wavelet subspace ProjVj,x ⊥ Vj+1,x which encodes the difference between Vj+1,x and
Vj,x is in general not orthogonal across scales.
Theorem 8 involves a quasi-orthogonality condition (8), which is satisfied if the operators
{Qj,k } applied on Mpare rapidly decreasing in norm or are orthogonal. When ρ ∈ A∞
1 such
−j
that kQj,k Xk ∼ 2
ρ(Cj,k ), quasi-orthogonality is guaranteed. In this case, for any node
p
p
0
Cj,k andPCj 0 ,k0 ⊂ Cj,k , we have kQj 0 ,k0 Xk/ ρ(Cj 0 ,k0 ) . 2−(j −j) kQj,k Xk/ ρ(Cj,k ), which
implies C 0 0 ⊂Cj,k hQj,k X, Qj 0 ,k0 Xi . 2kQj,k Xk2 . Therefore B0 . 2. Another setting is
j ,k
when Qj 0 ,k0 and Qj,k are orthogonal whenever Cj 0 ,k0 ⊂ Cj,k , as guaranteed in orthogonal
GMRA in Section 5.3, in which case exact orthogonality is automatically satisfied.
Quasi-orthogonality enters in the proof of Eq. (9). If quasi-orthogonality is violated, we
still have a convergence
result in Theorem 8 but the convergence rate will be worse: MSE
s
1
2s+d−2
. [(log n)/n]
when d ≥ 3 and MSE . [(logd n)/n] 2 when d = 1, 2.
5.3 Orthogonal GMRA and adaptive orthogonal GMRA
A different construction, called orthogonal geometric multi-resolution analysis in Section
5 of Allard et al. (2012), follows the classical wavelet theory by constructing a sequence
of increasing subspaces and then the corresponding wavelet subspaces exactly encode the
orthogonal complement across scales. Exact orthogonality is therefore satisfied.
5.3.1 Orthogonal GMRA
In the construction, we build the sequence of subspaces {Sbj,k }k∈Kj ,j≥jmin with a coarse-tofine algorithm in Table 4. For fixed x and j, (j, x) denotes (j, k) such that x ∈ Cj,k . In
orthogonal GMRA the sequence of subspaces Sj,x is increasing such that S0,x ⊂ S1,x ⊂
· · · Sj,x ⊂ Sj+1,x · · · and the subspace Uj+1,x exactly encodes the orthogonal complement of
33

Liao and Maggioni

Orthogonal GMRA

Empirical orthogonal GMRA

S0,x = V0,x

Sb0,x = Vb0,x

U1,x = ProjS ⊥ V1,x , S1,x = S0,x ⊕ U1,x
0,x

...

Subpaces

j,x

projectors

0,x

...
b
Uj+1,x = ProjSb⊥ Vbj+1,x

Uj+1,x = ProjS ⊥ Vj+1,x

j,x

Sj+1,x = Sj,x ⊕ Uj+1,x
P
Sj := k∈Kj Sj,k 1j,k

Affine

b1,x = Proj b⊥ Vb1,x , Sb1,x = Sb0,x ⊕ U
b1,x
U
S

Sj,k (x) := cj,k + ProjSj,k (x − cj,k )

bj+1,x
Sbj+1,x = Sbj,x ⊕ U
P
Sbj := k∈Kj Sbj,k 1j,k

Sbj,k (x) := b
cj,k + ProjSbj,k (x − b
cj,k )

Table 4: Orthogonal GMRA

Sj,x in Sj+1,x . Orthogonal GMRA with respect to the distribution ρ corresponds to affine
projectors onto the subspaces {Sj,k }k∈Kj ,j≥jmin .
For a fixed distribution ρ, the approximation error kX − Sj Xk decays as j increases.
We will consider the model class Aos where kX − Sj Xk decays like O(2−js ).
Definition 21 A probability measure ρ supported on M is in Aos if
|ρ|Aos = sup inf{Ao0 : kX − Sj Xk ≤ Ao0 2−js , ∀ j ≥ jmin } < ∞ ,

(27)

T

where T varies over the set, assumed non-empty, of multiscale tree decompositions satisfying
Assumption (A1-A5).
Notice that As ⊂ Aos . We split the MSE into the squared bias and the variance as: EkX−
Sbj Xk2 = kX − Sj Xk2 +EkSj X − Sbj Xk2 . The squared bias kX −Sj Xk2 ≤ |ρ|2Aos 2−2js when 5 j(d−2) 
d2 j 4 #Λj log[αj#Λj ]
j 2
ever ρ ∈ Ao . In Lemma 33 we show EkSj X − Sbj Xk2 ≤
=
O
2j
s

β2 n

n

where α and β are the constants in Lemma 14. A proper choice of the scale yields the
following result:

Theorem 22 Assume that ρ ∈ Aos , s ≥ 1. Let ν > 0 be arbitrary and µ > 0. If j ∗ is
properly chosen such that

µ logn n
for d = 1
∗
−j
 5  1
2
,
=
µ log n 2s+d−2
for
d
≥
2
n

then there exists a constant C1 (θ1 , θ2 , θ3 , θ4 , d, ν, µ, s) such that

5 
log
n
s
≤ C2 (θ1 , θ2 , θ3 , θ4 , d, µ)n−ν
P kX − Sbj ∗ Xk ≥ (|ρ|Aos µ + C1 )
n
(
)
s
 5  2s+d−2
log
n
P kX − Sbj ∗ Xk ≥ (|ρ|Aos µs + C1 )
≤ C2 (θ1 , θ2 , θ3 , θ4 , d, µ, s)n−ν
n
Theorem 22 is proved in appendix E.1.
34

for d = 1,
for d ≥ 2.
(28)

Adaptive Geometric Multiscale Approximations for Intrinsically Low-dimensional Data

Algorithm 2 Empirical Adaptive Orthogonal GMRA
Input: data X2n = Xn0 ∪ Xn , intrinsic dimension d, threshold κ
Output: SbΛb o : adaptive piecewise linear projectors
τn

1:

2:
3:

4:
5:

Construct T n and {Cj,k } from Xn0
b o on every node Cj,k ∈ T n .
Compute Sbj,k and ∆
j,k
b o ≥ 2−j τ o where
Tbτno ← smallest proper subtree of T n containing all Cj,k ∈ T n : ∆
n
j,k
q
5
o
τn = κ (log n)/n.
b τ o ← partition associated with the outer leaves of Tbτ o
Λ
n
n
P
b
Sbb ←
S
1
.
j,k j,k
b
Cj,k ∈Λτno

Λτno

5.3.2 Adaptive Orthogonal GMRA
Definition (infinite sample)
∆oj,k

Empirical version

:= k(Sj − Sj+1 )1j,k Xk

1
2
= k(I − Sj )1j,k Xk2 − k(I − Sj+1 )1j,k Xk2

b o := (Sbj − Sbj+1 )1j,k X
∆
j,k

2
=
(I − Sbj )1j,k X − (I − Sbj+1 )1j,k X

2

1
2

Table 5: Refinement criterion in adaptive orthogonal GMRA

Orthogonal GMRA can be constructed adaptively to the data with the refinement cri1
terion defined in Table 5. We let τno := κ(log5 n/n) 2 where κ is a constant, truncate the
data master tree T n to the smallest proper subtree that contains all C j,k ∈ T n satisfying
b o ≥ 2−j τno , denoted by Tbτ o . Empirical adaptive orthogonal GMRA returns piecewise
∆
j,k
n
b τ o consisting of the outer leaves of Tbτ o . Our
affine projectors on the adaptive partition Λ
n
n
algorithm is summarized in Algorithm 2.
If ρ is known, given any fixed threshold η > 0, we let T(ρ,η) be the smallest proper tree of
T that contains all Cj,k ∈ T for which ∆oj,k ≥ 2−j η. This gives rise to an adaptive partition
Λ(ρ,η) consisting the outer leaves of T(ρ,η) . We introduce a model class Bso for whose elements
we can control the growth rate of the truncated tree T(ρ,η) as η decreases.
Definition 23 In the case d ≥ 3, given s > 0, a probability measure ρ supported on M is
in Bso if the following quantity is finite
|ρ|pBso := sup sup η p
T

η>0

X

j≥jmin

2−2j #j T(ρ,η) with p =

2(d − 2)
2s + d − 2

(29)

where T varies over the set, assumed non-empty, of multiscale tree decompositions satisfying
Assumption (A1-A5).
35

Liao and Maggioni

Notice that exact orthogonality is satisfied for orthogonal GMRA. One can show that,
as long as ρ ∈ Bso ,


o
o
|ρ|2Bso 
kX − SΛ(ρ,η) Xk2 ≤ Bs,d
|ρ|pBso η 2−p ≤ Bs,d

X

j≥jmin

−

2s
d−2

2−2j #j T(ρ,η) 

,

o := 2p /(1 − 2p−2 ). We can prove the following performance guarantee of the
where Bs,d
empirical adaptive orthogonal GMRA (see Appendix E.2):

Theorem 24 Suppose M is bounded: M ⊂ BM (0) and the multiscale tree satisfies ρ(Cj,k ) ≤
θ0 2−jd for some θ0 > 0. Let d ≥ 3 and ν > 0. There exists κ0 (θ0 , θ2 , θ3 , θ4 , amax , d, ν) such

1
that if ρ ∈ Bso for some s > 0 and τno = κ (log5 n)/n 2 with κ ≥ κ0 , then there is a c1 and
c2 such that
(
)
s
 5  2s+d−2
log
n
P kX − SbΛb o Xk ≥ c1
≤ c2 n−ν .
(30)
τn
n
In Theorem 24, the constants are c1 := c1 (θ0 , θ2 , θ3 , θ4 , amax , d, s, κ, |ρ|Bso , ν) and c2 :=
2s

5

c2 (θ0 , θ2 , θ3 , θ4 , amin , amax , d, s, κ, |ρ|Bso ). Eq. (30) implies that MSE . ( logn n ) 2s+d−2 for
orthogonal adaptive GMRA when d ≥ 3. In the case of d = 1, 2, we can prove that
4+d
MSE . log n n .

Acknowledgments
This research was partially funded by ONR N00014-12-1-0601, NSF-DMS-ATD-1222567,
1708553, 1724979 and AFOSR FA9550-14-1-0033.

Appendix A. Tree construction, regularity of geometric spaces
A.1 Tree construction
We now show that from a set of nets {Tj (Xn0 )}j∈[jmin ,jmax ] from the cover tree algorithm we
N (j)

can construct a set of Cj,k with desired properties. Let {aj,k }k=1 be the set of points in
Tj (Xn0 ). Given a set of points {z1 , . . . , zm } ⊂ RD , the Voronoi cell of z` with respect to
{z1 , . . . , zm } is defined as
Voronoi(z` , {z1 , . . . , zm }) = {x ∈ RD : kx − z` k ≤ kx − zi k for all i 6= `}.
Let
f=
M

j[
max

[

j=jmin aj,k ∈Tj (Xn0 )

B 1 2−j (aj,k ) .

(31)

4

Our Cj,k ’s are constructed in Algorithm 3. These Cj,k ’s form a multiscale tree decomposition
f We will prove that M\M
f has a negligible measure and {Cj,k }k∈K ,j∈[j ,j ] satisfies
of M.
j

36

min

max

Adaptive Geometric Multiscale Approximations for Intrinsically Low-dimensional Data

Algorithm 3 Construction of a multiscale tree decomposition {Cj,k }
Input: data Xn0
Output: A multiscale tree decomposition {Cj,k }
1: Run cover tree on Xn0 to obtain a set of nets {Tj (Xn0 )}j∈[jmin ,jmax ]
f defined in (31)
2: j = jmin : Cjmin ,0 = M
3: for j = jmin + 1, . . . , jmax : For every Cj−1,k0 at scale j − 1, Cj−1,k0 has
# (Tj (Xn0 ) ∩ Cj−1,k0 ) children indexed by aj,k ∈ Tj (Xn0 ) ∩ Cj−1,k0 with corresponding
Cj,k ’s constructed as follows:
(j)
f
Cj,k = M

and for i = j + 1, . . . , jmax

(i)

Cj,k =

(j



\

Voronoi(aj,k , Tj (Xn0 ) ∩ Cj−1,k0 )
[

B 1 2−i (ai,k0 )
4

(i−1)

ai,k0 ∈Cj,k

[

(i−1)

Cj,k

)

Finally, let Cj,k = Cj,kmax .
Assumptions (A1-A5). The key is that every Cj,k is contained in a ball of radius 3 · 2−j and
also contains a ball of radius 2−j /4.
Lemma 25 Every Cj,k constructed in Algorithm 3 satisfies B 2−j (aj,k ) ⊆ Cj,k ⊆ B3·2−j (aj,k )
4

Proof For any x ∈ RD and any set C ∈ RD , the diameter of C with respect to x is defined
as diam(C, x) := supz∈C kz −xk. First, we prove that, for every j, Cj,k1 ∩Cj,k2 = ∅ whenever
k1 6= k2 . Take any aj+1,k10 ∈ Cj,k1 and aj+1,k20 ∈ Cj,k2 . Our construction guarantees that
diam(Cj+1,k10 , aj+1,k10 ) ≤ 41 2−(j+1) + 41 2−(j+2) + . . . < 12 2−(j+1)
and similarly for diam(Cj+1,k20 , aj+1,k20 ). Since kaj+1,k10 − aj+1,k20 k ≥ 2−(j+1) , this implies
that Cj+1,k10 ∩ Cj+1,k20 = ∅. In our construction,
Cj,k1 =



[

aj+1,k0 ∈Cj,k1
1

Cj+1,k10

[

B 2−j (aj,k1 ), Cj,k2 =
4



[

aj+1,k0 ∈Cj,k2
2

Cj+1,k20

[

B 2−j (aj,k2 ).
4

Since kaj,k1 − aj,k2 k ≥ 2−j , we observe that B 1 2−j (aj,k1 ) ∩ B 1 2−j (aj,k2 ) = ∅, Cj+1,k10 ∩
4
4
B 1 2−j (aj,k2 ) = ∅ for every aj+1,k10 ∈ Cj,k1 , and Cj+1,k20 ∩ B 1 2−j (aj,k1 ) = ∅ for every aj+1,k20 ∈
4
4
Cj,k2 . Therefore Cj,k1 ∩ Cj,k2 = ∅.
Our construction of Cj,k ’s guarantees that every Cj,k contains a ball of radius 41 · 2−j .
Next we prove that every Cj,k is contained in a ball of radius 3 · 2−j . The cover tree structure guarantees that Xn0 ⊂ ∪aj,k ∈Tj (Xn0 ) B2·2−j (aj,k ) for every j. Hence, for every aj,k and
every aj+1,k0 ∈ Cj,k , we obtain kaj+1,k0 − aj,k k ≤ 2 · 2−j and the computation above yields
diam(Cj+1,k0 , aj+1,k0 ) ≤ 2−j /4, and therefore diam(Cj,k , aj,k ) ≤ 2 · 2−j + 2−j /4 ≤ 3 · 2−j . In
37

Liao and Maggioni

summary Cj,k is contained in the ball of radius 3 · 2−j centered at aj,k .
The following Lemma will be useful when comparing comparing covariances of sets:
Lemma 26 If B ⊆ A, then we have λd (cov(ρ|A )) ≥

ρ(B)
ρ(A) λd (cov(ρ|B )).

Proof Without loss of generality, we assume both A and B are centered at x0 . Let V be
the eigenspace associated with the largest d eigenvalues of cov(ρ|B ). Then
uT cov(ρ|A )u
v T cov(ρ|A )v
λd (cov(ρ|A )) = max min
≥
min
v∈V
dimU =d u∈U
uT u
vT v

´
T
T
v
A (x − x0 )(x − x0 ) dρ v
≥ min
v∈V
ρ(A)v T v
´
 


´
T
T dρ v
T
T
v
(x
−
x
)(x
−
x
)
0
0
v
A\B
B (x − x0 )(x − x0 ) dρ v

= min 
+
T
T
v∈V
ρ(A)v v
ρ(A)v v

´
v T B (x − x0 )(x − x0 )T dρ v
ρ(B)
≥ min
=
λd (cov(ρ|B )).
v∈V
ρ(A)v T v
ρ(A)

A.2 Regularity of geometric spaces
To fix the ideas, consider the case where M is a manifold of class C s , s ∈ R+ \ Z, i.e.
around every point x0 there is a neighborhood Ux0 that is parametrized by a function
f : V → Ux0 , where V is an open connected set of Rd , and f ∈ C s , i.e. f is bsc times
continuously differentiable and the bsc-th derivative f bsc is Hölder continuous of order s −
bsc, i.e. ||f bsc (x) − f bsc (y)|| ≤ ||f bsc ||C s−bsc ||x − y||s−bsc . In particular, for s ∈ (0, 1), f is
simply a Hölder function of order s. For simplicity we assume x = f (xd ) where xd ∈ V .
If M is a manifold of class C s , s ∈ (0, 1), a constant approximation of f on a set I by
the value x0 := f (xd0 ) on such set yields
ˆ
ˆ
1
1
d
d 2
|f (x ) − f (x0 )| dρ(x) ≤
||xd − xd0 ||2s ||f ||2C s dρ(x) ≤ ||f ||2C s diam(I)2s
ρ(I) I
ρ(I) I
where we used continuity of f . If I was a ball, we would obtain a bound which would
be better by a multiplicative´ constant no larger than 1/d. Moreover, the left hand side is
1
minimized by the mean ρ(I)
I f (y)dρ(y) of f on I, and so the bound on the right hand side
holds a fortiori by replacing f (xd0 ) by the mean.
Next we consider the linear approximation of M on I ⊂ M. Suppose there exits θ0 , θ2
such that I is contained in a ball of radius θ2 r and contains a ball of radius θ0 r. Let x0 ∈ I be
the closest point on I to the mean. Then I is the graph of a C s function f : PTx0 (I) → PTx⊥ (I)
0

where Tx0 (I) is the plane tangent to I at x0 and Tx⊥0 (I) is the orthogonal complement of
Tx0 (I). Since all the quantities involved are invariant under rotations and translations, up
to a change of coordinates, we may assume xd = (x1 , . . . , xd ) and f = (f1 , . . . , fD−d ) where
38

Adaptive Geometric Multiscale Approximations for Intrinsically Low-dimensional Data

fi := fi (xd ), i = d+1, . . . , D. A linear approximation of f = (fd+1 , . . . , fD ) based on Taylor
expansion and an application of the mean value theorem yields the error estimates.
• Case 1: s ∈ (1, 2)
ˆ
2
1
f (xd ) − f (xd0 ) − ∇f (xd0 ) · (xd − xd0 ) dρ
ρ(I) I
ˆ
D
X
2
1
=
sup
∇fi (ξi )(xd − xd0 ) − ∇fi (xd0 ) · (xd − xd0 ) dρ
ρ(I) ξi ∈domain(fi ) I
i=d+1

≤

D
X

i=d+1

≤D

1
sup
ρ(I) ξi ∈domain(fi )

max

i=1,...,D−d

ˆ

Cj,k

||xd − xd0 ||2 kξi − xd0 k2(s−bsc) ||∇fi ||2C s−bsc dρ

||∇fi ||2C s−bsc diam(I)2s .

• Case 2: s = 2
ˆ
2
1
f (xd ) − f (xd0 ) − ∇f (xd0 ) · (xd − xd0 ) dρ
ρ(I) I
ˆ
D
X
1
=
fi (xd ) − fi (xd0 ) − ∇fi (xd0 ) · (xd − xd0 )
ρ(I) I
i=d+1

≤
≤

D
X

i=d+1

D
2

1
sup
ρ(I) ξi ∈domain(fi )

max

i=1,...,D−d

ˆ

I

2

dρ

1
(ξi − xd0 )T D2 fi |xd (ξi − xd0 ) + o(kξi − xd0 k2 )
0
2

2

dρ

kD2 fi kdiam(I)4 + o(2−4j ).

M does not have boundaries, so the Taylor expansion in the computations above can
be performed on the convex hull of PTx0 (I) , whose diameter is no larger than diam(I). Note
that this bound then holds for other linear approximations which are at least as good, in
L2 (ρ|I ), as Taylor expansion. One such approximation is, by definition, the linear least
square fit of f in L2 (ρ|I ). Let LI be the least square fit to the function x 7→ f (x). Then
D
X

i=d+1

≤



1
λi (cov(ρ|I )) =
ρ(I)

ˆ

2

I

||f (x) − LI (x)||2 dρ(x)

D maxi=1,...,D−d ||∇fi ||2C s−bsc diam(I)2s , s ∈ (1, 2)
.
4
D
2
s=2
2 maxi=1,...,D−d kD fi kdiam(I) ,

(32)

Proof [Proof of Proposition 13] Claim (A1) follows by a simple volume argument: Cj,k is
contained in a ball of radius 3 · 2−j , and therefore has volume at most C1 (3 · 2−j )d , and each
child contains a ball of radius 2−(j+1) /4, and therefore volume at least C1−1 (2−(j+1) /4)d . It
follows that amax ≤ C12 (3 · 2−j /2−(j+1) · 4)d . Clearly amin ≥ 1 since every aj,k belongs to
both Tj (Xn0 ) and Tj 0 (Xn0 ) with j 0 ≥ j. (A1),(A3), (A4) are straightforward consequences of
39

Liao and Maggioni

the doubling assumption and Lemma 25. As for (A2), for any ν > 0, we have




f = 0 and ρ(M \ M)
f > 28ν log n
f > 28ν log n = P ρb(M \ M)
P ρ(M \ M)
3n
3n


1
28ν log n
f
f
f
f
≤ P |b
ρ(M \ M) − ρ(M \ M)| > ρ(M \ M) and ρ(M \ M) >
2
3n
3

f

≤ 2e− 28 nρ(M\M) ≤ 2n−ν .

In order to prove the last statement about property (A5) in the case of 5a, observe that
B2−j /4 (aj,k ) ⊆ Cj,k ⊆ B3·2−j (aj,k ). By Lemma 26 we have
C1−1 (2−j /4)d
C1 (3 · 2−j )d
λd (cov(ρ|B2−j /4 (aj,k ) ) ≤ λd (cov(ρ|Cj,k ) ≤
λd (cov(ρ|B3·2−j (aj,k ) )
ρ(Cj,k )
ρ(Cj,k )
and therefore λd (cov(ρ|Cj,k ) ≥ C1−2 (1/12)d λd (cov(ρ|B2−j /4 (aj,k ) ) ≥ C1−2 (1/12)d θ̃3 (2−j /4)2 /d,

so that (A5)-(i) holds with θ3 = θ̃3 (4C1 )−2 (1/12)d . Proceeding similarly for λd+1 , we obtain
from the upper bound above that
λd+1 (cov(ρ|Cj,k ) ≤

C1 (3 · 2−j )d
λd+1 (cov(ρ|B3·2−j (aj,k ) ) ≤ (12d )2 · 144C14 θ̃4 /θ̃3 λd (cov(ρ|Cj,k )
C1−1 (2−j /4)d

so that (A5)-(ii) holds with θ4 = (12d )2 · 144C14 θ̃4 /θ̃3 .
In order to prove (A5) in the case of 5b, we use calculations as in Little et al. (2012);
Maggioni et al. (2016) where one obtains that the first d eigenvalues of the covariance matrix
of ρ|Br (z) with z ∈ M, is lower bounded by θ̃3 r2 /d for some θ̃3 > 0. Then (A5)-(i) holds for
Cj,k with θ3 = θ̃3 (4C1 )−2 (1/12)d . The estimate of λd+1 (cov(ρ|Cj,k )) follows from (32) such
that

D
X
D maxi=1,...,D−d ||∇fi ||2C s−bsc (6 · 2−j )2s , s ∈ (1, 2)
2
λi (cov(ρ|Cj,k )) ≤
.
D
2
−j 4
s=2
2 maxi=1,...,D−d kD fi k(6 · 2 ) ,
i=d+1

Therefore, there exists j0 such that λd+1 (cov(ρ|Cj,k )) < θ4 λd (cov(ρ|Cj,k )) when j ≥ j0 . The
2
calculation above also implies that ρ ∈ A∞
s if maxi=1,...,D−d ||∇fi ||C s−bsc for s ∈ (1, 2) or
maxi=1,...,D−d kD2 fi k for s = 2 is uniformly upper bounded.

A.3 An alternative tree construction method
The {Cj,k } constructed by Algorithm 3 is proved to satisfy Assumption (A1-A5). In numerical experiments, we use a much simpler algorithm to construct {Cj,k } as follows:
Cjmax ,k = Voronoi(ajmax ,k , Tjmax (Xn0 )) ∩ B2−jmax (ajmax ,k ),
S
and for any j < jmax , we define Cj,k =
Cj−1,k0 .
aj−1,k0 child of aj,k

40

Adaptive Geometric Multiscale Approximations for Intrinsically Low-dimensional Data

We observe that the vast majority of Cj,k ’s constructed above satisfy Assumption (A1A5) in our numerical experiments. In Fig. 13, we will show that (A5) is satisfied when
we experiment on volume measures on the 3-dim S and Z manifold. Here we sample 105
training data, perform multiscale tree decomposition as stated above, and compute θ3j,k , θ4j,k
at every Cj,k . In Fig. 13, we display the mean of {θ3j,k }k∈Kj or {θ4j,k }k∈Kj versus scale j,
with a vertical error bar representing the standard deviation of {θ3j,k }k∈Kj or {θ4j,k }k∈Kj at
each scale. We observe that θ3 = minj,k θ3j,k ≥ 0.05 at all scales and θ4 = maxj,k θ4j,k ≤ 1/2
except at very coarse scales, which demonstrates Assumption (A5) is satisfied here. Indeed
θ4 is not only bounded, but also decreases from coarse scales to fine scales.
θ3

0.3

θ4

1

θ3

0.25

θ4

1.2

1
0.8
0.25
0.2
0.8
0.6
0.2
0.6
0.4

0.15
0.4

0.15
0.2
0.2
0.1
0.1
0
0

0.05
-2

-1

0

1

2

3

-0.2
-2

4

scale

-1

0

1

2

3

0.05
-2

4

scale

-1

0

1

scale

2

3

4

-0.2
-2

-1

0

1

2

3

4

scale

(a) {θ3j,k }k∈Kj versus j of (b) {θ4j,k }k∈Kj versus j of (c) {θ3j,k }k∈Kj versus j of (d) {θ4j,k }k∈Kj versus j of
the 3-dim S manifold
the 3-dim S manifold
the 3-dim Z manifold
the 3-dim Z manifold

Figure 13: The mean of {θ3j,k }k∈Kj or {θ4j,k }k∈Kj versus scale j, with a vertical error bar
representing the standard deviation of {θ3j,k }k∈Kj or {θ4j,k }k∈Kj at each scale.
We observe that, every Cj,k constructed above is contained in a ball of radius θ2 2−j and
contains a ball of radius θ0 2−j , with θ2 /θ0 ∈ [1, 2] for the majority of Cj,k ’s. In Fig. 14, we
take the volume measures on the 3-dim S and Z manifold, and plot log2 of the outer-radius
and the statistics a lower bound for the in-radius3 versus the scale of cover tree. Notice
that the in-radius is a fraction of the outer-radius at all scales, and log2 θ2 − log2 θ0 ≤ 1 for
the majority of cells.
A.4 A∞
s ⊂ Bs

Proof [proof of Lemma 6] Assumepρ(Cj,k )  2−jd . According to Definition 2, ρ ∈ A∞
s
−js ρ(C ), ∀k ∈ K , j ≥ j
2
if k(X − Pj,k X)1j,k k ≤ |ρ|A∞
,
which
implies
∆
≤
j
min
j,k
j,k
s
p
d
−j(s+
)
−js
2
2|ρ|A∞
2
ρ(Cj,k ) . |ρ|A∞
2
.
s
s
Let η > 0 and T(ρ,η) be the smallest proper subtree of T that contains all Cj,k for which
∆j,k

d

2−j(s+ 2 ) & 2−j η
≥ 2−j η. All the nodes satisfying ∆j,k ≥ 2−j η will satisfy |ρ|A∞
s
2

which implies 2−j & (η/|ρ|A∞
) 2s+d−2 . Therefore, the truncated tree T(ρ,η) is contained in
s
∗

2

) 2s+d−2 , so the entropy of T(ρ,η) is upper bounded by
Tj ∗ = ∪j≤j ∗ Λj with 2−j  (η/|ρ|A∞
s
3. The in-radius of Cj,k is approximately computed as follows: we randomly pick a center, and evaluate the
largest radius with which the ball contains at least 95% points from Cj,k . This procedure is repeated for
two centers, and then we pick the maximal radius as an approximation of the in-radius.

41

Liao and Maggioni

log2(radii) versus scale j

1

log2(radii) versus scale j

1

Ratio between outer-radii and in-radii
4.5

outer-radii:
slope = -0.88764
in-radii 0.95 quantile: slope = -1.0374
in-radii mean:
slope = -1.0001
in-radii median:
slope = -1.0001
in-radii 0.05 quantile: slope = -0.93474

0

outer-radii:
slope = -0.86723
in-radii 0.95 quantile: slope = -1.1825
in-radii mean:
slope = -1.0797
in-radii median:
slope = -1.0797
in-radii 0.05 quantile: slope = -0.90587

0

4

3.5

3

2.5

2

-1

-1

1.5

log2(radius)

log2(radius)

1

-2

0.5
-2

-1

0

1

2

3

4

5

4

5

scale

-2

Ratio between outer-radii and in-radii
4

3.5

-3

-3

3

2.5

2

-4

-4
1.5

1

-5
0

0.5

1

1.5

2

2.5

-5
-0.5

3

0.5
-2

0

scale

0.5

1

1.5

2

2.5

-1

0

1

2

3

scale

scale

Figure 14: From left to right: the in-radius and outer-radius of a pentagon; log2 of the
outer-radius and the statistics of the in-radius versus the scale of cover tree for
the 3-dim S manifold, and then the same plot for the 3-dim Z manifold; ratio
between outer-radii and in-radii, for the 3-dim S manifold top) and the 3-dim Z
manifold (bottom).

2(d−2)
P
∗
the entropy of Tj ∗ , which is j≤j ∗ 2−2j #Λj  2j (d−2)  (η/|ρ|A∞
)− 2s+d−2 . Then ρ ∈ Bs
s
and |ρ|Bs . |ρ|A∞
according to Definition 5.
s

Appendix B. S manifold and Z manifold
We consider volume measures on the d dimensional S manifold and Z manifold whose x1
and x2 coordinates are on the S curve and Z curve in Figure 5 (a) and xi , i = 3, . . . , d + 1
are uniformly distributed in [0, 1].
B.1 S manifold
Since S manifold is smooth and has a bounded curvature, the volume measure on the S
manifold is in A∞
2 . Therefore, the volume measure on the S manifold is in A2 and B2 when
d ≥ 3.
B.2 Z manifold
B.2.1 The volume on the Z manifold is in A1.5
The uniform distribution on the d dimensional Z manifold is in A1 at two corners and
satisfies k(X − Pj,k X)1j,k k = 0 when p
Cj,k is away from the corners. There exists A0 > 0
such that k(X − Pj,k X)1j,k k ≤ A0 2−j ρ(Cj,k ) when Cj,k intersects with the corners. At
scale j, there are about 2jd cells away from the corners and there are about 2j(d−1) cells
which intersect with the corners. As a result,
kX − Pj Xk ≤ O

p

2jd · 0 · 2−jd + 2j(d−1) · 2−2j · 2−jd = O(2−1.5j ),

so the volume measure on Z manifold is in A1.5 .
42

Adaptive Geometric Multiscale Approximations for Intrinsically Low-dimensional Data

B.2.2 Model class Bs

Assume ρ(Cj,k )  2−jd . We compute the regularity parameter s in the Bs model class
when d ≥ 3. It is easy to see that ∆j,k = 0 when Cj,k is away from the corners and
p
d
∆j,k ≤ 2A0 2−j ρ(Cj,k ) . 2−j( 2 +1) when Cj,k intersects with the corners. Given any fixed
threshold η > 0, in the truncated tree T(ρ,η) , the parent of the leaves intersecting with the
d

corners satisfy 2−j( 2 +1) & 2−j η. In other words, at the corners the tree is truncated at a
2
∗
scale coarser than j ∗ such that 2−j = O(η d ). Since ∆j,k = 0 when Cj,k is away from the
corners, the entropy of T(ρ,η) is dominated by the nodes intersecting with the corners whose
cardinality is 2j(d−1) at scale j. Therefore
Entropy of T(ρ,η) .

X

j≤j ∗

 2(d−3) 
2−2j 2j(d−1) = O η − d
,

and s ≥ 3(d−2)
which implies that p ≤ 2(d−3)
d
2(d−3) > 1.5.
Then we study the relation between the error kX − PΛ(ρ,η) Xk and the partition size
#Λ(ρ,η) , which is numerically verified in Figure 4. Since all the nodes in T(ρ,η) that intersect
with corners are at a scale coarser than j ∗ , #Λ(ρ,η) ≈ 2j
η . [#Λ(ρ,η) ]

d
− 2(d−1)

∗ (d−1)

 η−

and

kX − PΛ(ρ,η) Xk . η

2−p
2

2sd
− 2(d−1)(2s+d−2)

2s

= η 2s+d−2 . [#Λ(ρ,η) ]

2(d−1)
d

. Therefore,

3
− 2(d−1)

= [#Λ(ρ,η) ]

.

Appendix C. Proofs of Lemma 14 and Proposition 15
C.1 Concentration inequalities
We first recall a Bernstein inequality from Tropp (2014) which is an exponential inequality
to estimate the spectral norm of a sum independent random Hermitian matrices of size
D × D. It features the dependence on an intrinsic dimension parameter which is usually
much smaller than the ambient dimension D. For a positive-semidefinite matrix A, the
intrinsic dimension is the quantity
intdim(A) =

trace(A)
.
kAk

Proposition 27 (Theorem 7.3.1 in Tropp (2014)) Let ξ1 , . . . , ξn be D × D independent random Hermitian matrices that satisfy
Eξi = 0 and kξi k ≤ R, i = 1, . . . , n.

P
Form the mean ξ = n1 ni=1 ξi . Suppose E(ξ 2 )  Φ. Introduce the intrinsic dimension
parameter din = intdim(Φ). Then, for nt ≥ nkΦk1/2 + R/3,
P{kξk ≥ t} ≤ 4din e
43

nt2 /2

− nkΦk+Rt/3

.

Liao and Maggioni

We use the above inequalities to estimate the deviation of the empirical mean from the
mean and the deviation of the empirical covariance matrix from the covariance matrix when
the data Xj,k = {x1 , . . . , xn } (with a slight abuse of notations) are i.i.d. samples from the
distribution ρ|Cj,k .
Lemma 28 Suppose x1 , . . . , xn are i.i.d. samples from ρ|Cj,k . Let
ˆ

n

cj,k =
Cj,k

ˆ
Σj,k =
Cj,k

, b
cj,k

xdρ|Cj,k

1X
:=
xi
n

b j,k := 1
, Σ
n

(x − cj,k )(x − cj,k )T dρ|Cj,k

Then
−

3nt2
2 −2j

i=1
n
X
i=1

(xi − b
cj,k )(xi − b
cj,k )T

−j

P{kb
cj,k − cj,k k ≥ t} ≤ 8e 6θ2 2 +2θ2 2 t ,
 2

3nt2
−
4θ2
4 2−4j +8θ 2 2−2j t
24θ2
b
2
P{kΣj,k − Σj,k k ≥ t} ≤
d+8 e
.
θ3

(33)
(34)

Proof We start by proving (33). We will apply Bernstein inequality with ξi = xi − cj,k ∈
RD . Clearly
Eξi = 0, and kξi k ≤ θ2 2−j due to Assumption (A4). We form the mean
P
n
cj,k − cj,k and compute the variance
ξ = n1 i=1 ξi = b
2

2

T

σ = n kEξ ξk = E

n
X
i=1

xi − cj,k

!T

n
X
i=1

xi − cj,k

!

n
X

=

i=1

E(xi − cj,k )T (xi − cj,k ) ≤ nθ22 2−2j .

Then for nt ≥ σ + θ2 2−j /3,
P{kb
cj,k − cj,k k ≥ t} ≤ 8e

−

n2 t2 /3
σ 2 +θ2 2−j nt/3

−

≤ 8e

3nt2
2 2−2j +2θ 2−j t
6θ2
2

.

Pn

We now prove (34). Define the intermediate matrix Σ̄j,k = n1 i=1 (xi − cj,k )(xi − cj,k )T .
b j,k − Σj,k = Σ̄j,k − Σj,k − (b
Since Σ
cj,k − cj,k )(b
cj,k − cj,k )T , we have
b j,k − Σj,k k ≤ kΣ̄j,k − Σj,k k + kb
kΣ
cj,k − cj,k k2 ≤ kΣ̄j,k − Σj,k k + θ2 2−j kb
cj,k − cj,k k.

b j,k −Σj,k k < t is kΣ̄j,k −Σj,k k < t/2 and kb
A sufficient condition for kΣ
cj,k −cj,k k < 2j t/(2θ2 ).
We apply Proposition 27 to estimate P{kΣ̄j,k − Σj,k k ≥ t/2}: let ξi = (xi − cj,k )(xi −
cj,k )T −
Σj,k ∈ RD×D . One can verify that Eξi = 0 and kξi k ≤ 2θ22 2−2j . We form the mean
P
ξ = n1 ni=1 ξi = Σ̄j,k − Σj,k , and then
!
n
n
n
n
X
X
1
1 X 2
1 X 2 −2j
θ22 2−2j
2
Eξ = E
ξ
ξ
=
Eξ

θ
2
Σ

Σj,k ,
i
i
j,k
i
2
n2
n2
n2
n
i=1

which satisfies

θ22 2−2j
n Σj,k

i=1

i=1

i=1

≤ θ24 2−4j /n. Meanwhile

din = intdim(Σj,k ) =

trace(Σj,k )
θ2 2−2j
θ2
≤ 2 −2j = 2 d.
kΣj,k k
θ3 2 /d
θ3
44

Adaptive Geometric Multiscale Approximations for Intrinsically Low-dimensional Data

Then, Proposition 27 implies
−nt2 /8

2

−3nt
2 −2j t
4θ2
4θ2
4 −4j θ2 2
4 −4j
2 −2j
P{kΣ̄j,k − Σj,k k ≥ t/2} ≤ 2 de θ2 2 + 3
= 2 de 24θ2 2 +8θ2 2 t .
θ3
θ3

Combining with (33), we obtain


2j t
b
P{kΣj,k − Σj,k k ≥ t} ≤ P{kΣ̄j,k − Σj,k k ≥ t/2} + P kb
cj,k − cj,k k ≥
2θ2
 2

3nt2
−
4θ2
4 −4j
2 −2j
≤
d + 8 e 24θ2 2 +8θ2 2 t .
θ3

In Lemma 28 data are assumed to be i.i.d. samples from the conditional distribution
ρ|Cj,k . Given Xn = {x1 , . . . , xn } which contains i.i.d. samples from ρ, we will show that the
empirical measure ρb(Cj,k ) = n
bj,k /n is close to ρ(Cj,k ) with high probability.

Lemma 29 Suppose x1 , . . . , xn are i.i.d. samples from ρ. Let ρ(Cj,k ) =
ρb(Cj,k ) = n
bj,k /n where n
bj,k is the number of points in Cj,k . Then

´

Cj,k

1dρ and

2

P{|b
ρ(Cj,k ) − ρ(Cj,k )| ≥ t} ≤ 2e

for all t ≥ 0. Setting t = 12 ρ(Cj,k ) gives rise to


1
P |b
ρ(Cj,k ) − ρ(Cj,k )| ≥ ρ(Cj,k )
2



− 6ρ(C3nt )+2t
j,k

3

≤ 2e− 28 nρ(Cj,k ) .

(35)

(36)

Combining Lemma 28 and Lemma 29 gives rise to probability bounds on kb
cj,k −cj,k k and
b j,k − Σj,k k where cj,k , b
b j,k are the conditional mean, empirical conditional
kΣ
cj,k , Σj,k and Σ
mean, conditional covariance matrix and empirical conditional covariance matrix on Cj,k ,
respectively.
b j,k as
Lemma 30 Suppose x1 , . . . , xn are i.i.d. samples from ρ. Define cj,k , Σj,k and b
cj,k , Σ
Table 1. Then given any t > 0,
3
− 28
nρ(Cj,k )

−

3nρ(Cj,k )t2
2 2−2j +4θ 2−j t
12θ2
2

P {kb
cj,k − cj,k k ≥ t} ≤ 2e
+ 8e
,


3nρ(Cj,k )t2
n
o
−
3
4θ22
4 2−4j +16θ 2 2−2j t
− 28
nρ(Cj,k )
96θ2
b
2
P kΣj,k − Σj,k k ≥ t
≤ 2e
+
d+8 e
.
θ3

(37)
(38)

P
Proof The number of samples on Cj,k is n
bj,k = ni=1 1j,k (xi ). Clearly E[b
nj,k ] = nρ(Cj,k ).
Let I ⊂ {1, . . . , n} and |I| = s. Conditionally on the event AI := {xi ∈ Cj,k for i ∈
45

Liao and Maggioni

I and xi ∈
/ Cj,k for i ∈
/ I}, the random variables {xi , i ∈ I} are i.i.d. samples from ρ|Cj,k .
According to Lemma 29,
P{kb
cj,k − cj,k k ≥ t | n
bj,k = s} =

X

I⊂{1,...,n}
|I|=s

= P{kb
cj,k − cj,k k ≥ t | A{1,...,s} } ≤ 8e
and
b j,k − Σj,k k ≥ t | n
P{kΣ
bj,k = s} ≤



−

P{kb
cj,k − cj,k k ≥ t | AI }

3st2
2 2−2j +2θ 2−j t
6θ2
2

1

n
s

,


3st2
−
4θ22
4 2−4j +8θ 2 2−2j t
24θ2
2
d+8 e
.
θ3

bj,k ≥ 21 nρ(Cj,k ) and then
Furthermore |b
ρ(Cj,k ) − ρ(Cj,k )| ≤ 12 ρ(Cj,k ) yields n


P kb
cj,k − cj,k k ≥ t



b j,k − Σj,k k ≥ t
P kΣ


3nρ(Cj,k )t2
−
1
2 2−2j +4θ 2−j t
12θ2
2
,
|b
ρ(Cj,k ) − ρ(Cj,k )| ≤ ρ(Cj,k ) ≤ 8e
2

(39)

  2

3nρ(Cj,k )t2
−
1
4θ2
4 2−4j +16θ 2 2−2j t
48θ2
2
.(40)
|b
ρ(Cj,k ) − ρ(Cj,k )| ≤ ρ(Cj,k ) ≤
d+8 e
2
θ3

Eq. (39) (40) along with Lemma 29 gives rise to
P {kb
cj,k − cj,k k ≥ t} ≤ 2e

3
− 28
nρ(Cj,k )

−

+ 8e

3nρ(Cj,k )t2
2 2−2j +4θ 2−j t
12θ2
2

,

 2

3nρ(Cj,k )t2
n
o
−
3
4θ
4
2 −2j
−4j
2
b j,k − Σj,k k ≥ t ≤ 2e− 28 nρ(Cj,k ) +
P kΣ
d + 8 e 48θ2 2 +16θ2 2 t .
θ3
b j,k − Σj,k k, we can estimate the angle between the eigenspace of Σ
b j,k and Σj,k
Given kΣ
with the following proposition.

Proposition 31 (Davis and Kahan (1970) or Theorem 3 in Zwald and Blanchard (2006))
j,k
1
b
Let δd (Σj,k ) = 12 (λj,k
d − λd+1 ). If kΣj,k − Σj,k k ≤ 2 δd (Σj,k ), then
ProjVj,k − ProjVbj,k ≤

b j,k − Σj,k k
kΣ
.
δd (Σj,k )

According to Assumption (A4) and (A5), δd (Σj,k ) ≥ θ3 2−2j /(4d). An application of Proposition 31 yields


n
θ3 (1 − θ4 )t
b
P kProjVj,k − ProjVbj,k ≥ t} ≤ P kΣj,k − Σj,k k ≥
2d22j
2
2
2
 2

3θ3 (1−θ4 ) nρ(Cj,k )t
−
3
4θ2
4 d2 +32θ 2 θ (1−θ )dt
− 28
nρ(Cj,k )
384θ2
4
2 3
≤ 2e
+
d+8 e
.
(41)
θ3
46

Adaptive Geometric Multiscale Approximations for Intrinsically Low-dimensional Data

Proof [Proof of Lemma 14] Since
X ˆ

bΛ Xk2 =
kPΛ X − P

Cj,k ∈Λ Cj,k

bj,k xk2 dρ =
kPj,k x − P

X
j

X

ˆ

bj,k xk2 dρ,
kPj,k x − P

k:Cj,k ∈Λ Cj,k

we obtain the estimate

n
o X  X
bΛ Xk ≥ η ≤
P kPΛ X − P
P

j

ˆ

k:Cj,k ∈Λ Cj,k

bj,k xk2 dρ ≥ P
kPj,k x − P

P

ˆ



k:Cj,k ∈Λ Cj,k

j

Λη 2




−2j # Λ 
j
j≥jmin 2

Next we prove that, for any fixed scale j,

 X

2−2j #

bj,k xk2 dρ ≥ t2
kPj,k x − P





.(42)

2j nt2

≤ α#j Λe

− β2#

jΛ

.

(43)

P
Then Lemma 14 is proved by setting t2 = 2−2j #j Λη 2 /( j≥0 2−2j #j Λ).
The proof of (43) starts with the following calculation:
X

ˆ

bj,k xk2 dρ
kPj,k x − P

k:Cj,k ∈Λ Cj,k

X

=

ˆ

k:Cj,k ∈Λ Cj,k

≤

X

ˆ

k:Cj,k ∈Λ Cj,k

≤2
≤2

ˆ

X

kcj,k + ProjVj,k (x − cj,k ) − b
cj,k − ProjVbj,k (x − b
cj,k )k2 dρ

k(I − ProjVbj,k )(cj,k − b
cj,k ) + (ProjVj,k − ProjVbj,k )(x − cj,k )k2 dρ

k:Cj,k ∈Λ Cj,k

ˆ

X

k:Cj,k ∈Λ Cj,k

h

h

i
kcj,k − b
cj,k k2 + k(ProjVj,k − ProjVbj,k )(x − cj,k )k2 dρ
i
kcj,k − b
cj,k k2 + θ22 2−2j kProjVj,k − ProjVbj,k k2 dρ

+
For any fixed j and given t > 0, we divide Λ into light cells Λ−
j,t and heavy cells Λj,t , where

Λ−
j,t

:=



Cj,k ∈ Λ : ρ(Cj,k ) ≤

t2
20θ22 2−2j #j Λ



−
and Λ+
j,t := Λ \ Λj,t .

h
i
´
Since Cj,k kcj,k − b
cj,k k2 + θ22 2−2j kProjVj,k − ProjVbj,k k2 dρ ≤ 5θ22 2−2j ρ(Cj,k ), for light sets
we have
2

X

k:Cj,k ∈Λ−
j,t

ˆ

h
i
t2
kcj,k − b
cj,k k2 + θ22 2−2j kProjVj,k − ProjVbj,k k2 dρ ≤ .
2
Cj,k
47

(44)

Liao and Maggioni

Next we consider Cj,k ∈ Λ+
j,t . We have
(

t
P kb
cj,k − cj,k k ≥ p
8#j Λρ(Cj,k )


3
≤ 2 exp − nρ(Cj,k ) + 8e
28

)

t2
3nρ(Cj,k )
8#j Λρ(Cj,k )
−
t
2 2−2j +4θ 2−j
12θ2
2
8#j Λρ(Cj,k )

√

2j

≤ C1 e

2

−C2 2 # nt
Λ
j

,

(45)

and
(

P kProjVj,k

2j t
− ProjVbj,k k ≥ p
θ2 8#j Λρ(Cj,k )

3
− 28
nρ(Cj,k )

≤ 2e

+




4θ22
d+8 e
θ3

)

2 (1−θ )2 nρ(C
22j t2
3θ3
4
j,k ) 8θ 2 # Λρ(C
j
j,k )
2
−
2j t
4 d2 +32θ 2 θ (1−θ )d
384θ2
3
4
2
θ2
8#j Λρ(Cj,k )

√

2j nt2
d #j Λ

−C4 22

≤ C3 de

(46)

where positive constants C1 , C2 , C3 , C4 depend on θ2 and θ3 . Combining (44), (45) and (46)
gives rise to (43) with α = max(C1 , C3 ) and β = min(C2 , C4 ).

Proof [Proof of Proposition 15] The bound (18) follows directly from Lemma 14 applied
to Λ = Λj ; (19) follows from (18) by integrating the probability over η:

bj Xk2 =
EkPj X − P
(
ˆ
+∞

≤

ˆ
0

+∞

n
o
bj Xk ≥ η dη
ηP kPj X − P
)
ˆ
ˆ
2j
2

η min 1, αd#Λj e

0

where αd#Λj e

−

2
β22j nη0
d2 #Λj

− β22

nη
d #Λj

η0

dη =

+∞

ηdη +

0

2j nη 2
d #Λj

− β22

αdη#Λj e

η0

= 1. Then
2j

2

#Λ2j −β 2 #Λnη0
d2 #Λj log[αd#Λj ]
bj Xk2 = 1 η02 + α ·
j
EkPj X − P
e
≤
.
2
2β 22j n
β22j n

48

dη

Adaptive Geometric Multiscale Approximations for Intrinsically Low-dimensional Data

Appendix D. Proof of Eq. (9), Lemma 16, 18, 19, 20
D.1 Proof of Eq. (9)
+(n−1)

+n
Let Λ+0
(ρ,η) = Λ(ρ,η) and Λ(ρ,η) be the partition consisting of the children of Λ(ρ,η)
n = 1, 2, . . .. Then

kX − PΛ(ρ,η) Xk = k
=k
≤k

n−1
X
`=0
∞
X
`=0

for

(PΛ+` X − PΛ+(`+1) X) + PΛ+n (ρ,η) X − Xk
(ρ,η)

(ρ,η)

(PΛ+` X − PΛ+(`+1) X) + lim PΛ+n X − Xk
(ρ,η)

X

Cj,k ∈T
/ (ρ,η)

n→∞

(ρ,η)

(ρ,η)

Qj,k Xk + k lim PΛ+n X − Xk.
n→∞

(ρ,η)

We have k limn→∞ PΛ+n X − Xk = 0 due to Assumption (A4). Therefore,
(ρ,η)

kX − PΛ(ρ,η) Xk2 ≤ k
≤ B0
≤ B0

X

X

Cj,k ∈T
/ (ρ,η)

X

Qj,k Xk2 ≤

`≥0 Cj,k ∈T
\T(ρ,2−` η)
(ρ,2−(`+1) η)

X

2−2` η 2

X

j≥jmin

`≥0

X

Cj,k ∈T
/ (ρ,η)

∆2j,k ≤ B0

X X

B0 kQj,k Xk2 = B0

`≥0 j≥jmin

2−2j #j T(ρ,2−(`+1) η) ≤ B0

X
`≥0

X

∆2j,k

Cj,k ∈T
/ (ρ,η)

(2−j 2−` η)2 #j T(ρ,2−(`+1) η)

2−2` η 2 |ρ|pBs [2−(`+1) η]−p



X
≤ B0 2p 
2−`(2−p)  |ρ|pBs η 2−p ≤ Bs,d |ρ|pBs η 2−p .
`≥0

D.2 Proof of Lemma 16

2

(X − Pj ∗ X)1nC


j ∗ ,k :ρ(Cj ∗ ,k )≤

≤ # Cj ∗ ,k : ρ(Cj ∗ ,k ) ≤
≤

28(ν+1) log n
3n

o



≤

n

X

Cj ∗ ,k :ρ(Cj ∗ ,k )≤

ˆ
o Cj ∗ ,k
28(ν+1) log n

kx − Pj ∗ ,k k2 dρ

3n

28(ν + 1) log n 2 −2j ∗ 28(ν + 1) log n
θ2 2
3n
3n

28(ν+1)θ22 j ∗ (d−2)
2
(log n)/n
3θ1

≤

28(ν+1)θ22 µ
3θ1

((log n)/n)2 .

For every Cj ∗ ,k , we have

P ρ(Cj ∗ ,k ) > 28
b(Cj ∗ ,k ) < d/n
3 (ν + 1)(log n)/n and ρ

≤ P |b
ρ(Cj ∗ ,k ) − ρ(Cj ∗ ,k )| > ρ(Cj ∗ ,k )/2 and ρ(Cj ∗ ,k ) >
for n so large that 14(ν + 1) log n > 3d
3

≤ 2e− 28 nρ(Cj ∗ ,k ) ≤ 2n−ν−1 .
49

28
3 (ν

+ 1)(log n)/n

Liao and Maggioni

Then

P each Cj ∗ ,k satisfying ρ(Cj ∗ ,k ) > 28
3 (ν + 1)(log n)/n has at most d points

28
≤ # Cj ∗ ,k : ρ(Cj ∗ ,k ) < 3 (ν + 1)(log n)/n 2n−ν−1 ≤ #Λj ∗ 2n−ν−1 ≤ 2n−ν /(θ1 µ log n) < n−ν ,
when n is so large that θ1 µ log n > 2.
D.3 Proof of Lemma 18
Since Tbτn ⊂ T(ρ,bτn ) , P{e12 > 0} if and only if there exists Cj,k ∈ T(ρ,bτn ) \ Tbτn . In other
words, P{e12 > 0} if and only if there exists Cj,k ∈ T(ρ,bτn ) such that ρb(Cj,k ) < d/n and
∆j,k > 2−j bτn . Therefore,
X
P{e12 > 0} ≤
P{b
ρ(Cj,k ) < d/n and ∆j,k > 2−j bτn }
≤
≤

X

Cj,k ∈T(ρ,bτn )

X

Cj,k ∈T(ρ,bτn )

Cj,k ∈T(ρ,bτn )

n
P ρb(Cj,k ) < d/n and ρ(Cj,k ) >

4b2 τn2
9θ22



o

since ∆j,k ≤

n
P |b
ρ(Cj,k ) − ρ(Cj,k )| > ρ(Cj,k )/2 and ρ(Cj,k ) >

4b2 τn2
9θ22

o

−j
3
2 θ2 2


q
ρ(Cj,k )

(for n large enough so that 2b2 κ2 log n > 9θ22 d)
2 2
2 2
3
X
− 28
n· 4b κ 2log n
− b κ2
9θ2 n
21θ2
#T(ρ,bτn ) .
≤
2e
≤ 2n
Cj,k ∈T(ρ,bτn )

The leaves of T(ρ,bτn ) satisfy ρ(Cj,k ) > 4b2 τn2 /(9θ22 ). Since ρ(M) = 1, there are at most
9θ22 /(4b2 τn2 ) leaves in T(ρ,bτn ) . Meanwhile, since every node in T has at least amin children,
#T(ρ,bτn ) ≤ 9θ22 amin /(4b2 τn2 ). Then for a fixed but arbitrary ν > 0,
b2 κ2

P{e12 > 0} ≤

18θ22 amin − 21θ2
2
n
4b2 τn2

b2 κ2

≤

18θ22 amin 1− 21θ2
2
n
4b2 κ2

≤ C(θ2 , amax , amin , κ)n−ν ,

if κ is chosen such that κ > κ1 where b2 κ21 /(21θ22 ) = ν + 1.
D.4 Proof of Lemma 19
We first prove (24). Introduce the intermediate variable
¯ j,k := kQj,k kn = k(Pj − Pj+1 )1j,k Xk
∆
n
and then observe that
n
o
n
o
b j,k ≤ η and ∆j,k ≥ bη ≤ P ∆
b j,k ≤ η and ∆
¯ j,k ≥ (amax + 2)η
P ∆

¯ j,k ≤ (amax + 2)η and ∆j,k ≥ (2amax + 5)η . (47)
+P ∆

The bound in Eq. (24) is proved in the following three nsteps. In Step One, we show that o
b j,k ≤ η and ∆
¯ j,k ≥ (amax + 2)η
∆j,k ≥ bη implies ρ(Cj,k ) ≥ O(22j η 2 ). Then we estimate P ∆

¯ j,k ≤ (amax + 2)η and ∆
¯ j,k ≥ (2amax + 5)η in Step Three.
in Step Two and P ∆
50

Adaptive Geometric Multiscale Approximations for Intrinsically Low-dimensional Data

Step One: Notice that ∆j,k ≤ 32 θ2 2−j

p
ρ(Cj,k ). As a result, ∆j,k ≥ bη implies

ρ(Cj,k ) ≥

4b2 22j η 2
.
9θ22

(48)

Step Two:
n
o
n
o
b j,k ≤ η and ∆
¯ j,k ≥ (amax + 2)η ≤ P |∆
b j,k − ∆
¯ j,k | ≥ (amax + 1)η .
P ∆

(49)

We can write

b j,k − ∆
¯ j,k | ≤ (Pj,k − P
bj,k )1j,k X
|∆

n

X

+

Cj+1,k0 ∈C (Cj,k )

bj+1,k0 − Pj+1,k0 )1j+1,k0 X
(P


q
≤ kcj,k − b
cj,k k + θ2 2−j kProjVj,k − ProjVbj,k k
ρb(Cj,k )
|
{z
}
X

+

Cj+1,k0 ∈C (Cj,k )

|



e1

kcj+1,k0 − b
cj+1,k0 k + θ2 2−(j+1) kProjVj+1,k0 − ProjVb

j+1,k

{z
e2

k
0

n

q
ρb(Cj+1,k0 ) .
}

(50)

Term e1 : We will estimateP{e1 > η}. Conditional on the event that{|b
ρ(Cj,k ) − ρ(Cj,k )| ≤
p
1
3
−j
cj,k k + θ2 2 kProjVj,k − ProjVbj,k k
ρ(Cj,k ). A similar
2 ρ(Cj,k )}, we have e1 ≤ 2 kcj,k − b
argument to the proof of Lemma 14 along with (48) give rise to
 

q
3
2j
2
−j
P
ρ(Cj,k ) > η ≤ γ̃1 e−γ̃2 2 nη
kcj,k − b
cj,k k + θ2 2 kProjVj,k − ProjVbj,k k
2


where γ̃1 := γ̃1 (θ2 , θ3 , d) and γ̃2 := γ̃2 (θ2 , θ3 , θ4 , d); otherwise P |b
ρ(Cj,k ) − ρ(Cj,k )| > 12 ρ(Cj,k ) ≤
3
− 28
nρ(Cj,k )

2e

−b

≤ 2e

2 22j nη 2
2
21θ2

. Therefore

P{e1 > η} ≤ max(γ̃1 , 2)e

− min(γ̃2 ,

b2
2j
2
2 )2 nη
21θ2

(51)

n
Term e2 : We will estimate P{e2 > amax η}. Let Λ− = Cj+1,k0 ∈ C (Cj,k ) : ρ(Cj+1,k0 ) ≤

22j η 2
8θ22

+
−
−
and
n Λ = C (C2jj,k2) \ Λ . For every C2jj+1,k
o0 ∈ Λ , when we condition on the event that
2
ρ(Cj+1,k0 ) ≤ 28θη2 and ρb(Cj+1,k0 ) ≤ 24θη2 , we obtain
2

X

Cj+1,k0 ∈Λ−

≤

X



2

kc

Cj+1,k0 ∈Λ−

j+1,k0

θ2 2

−j

−b
c

j+1,k0

k + θ2 2

−(j+1)

q
ρb(Cj,k ) ≤ amax η/2;

kProjVj+1,k0 − ProjVb

51

j+1,k

q
k
ρb(Cj+1,k0 )
0

(52)

o

Liao and Maggioni

otherwise,
n
o
2j 2
2j 2
P ρ(Cj+1,k0 ) ≤ 28θη2 and ρb(Cj+1,k0 ) > 24θη2
2
2
n
22j η 2
ρ(Cj+1,k0 ) − ρ(Cj+1,k0 )| ≥
≤ P ρ(Cj+1,k0 ) ≤ 8θ2 and |b
2

2 !

22j η 2
22j η 2
2j
2
≤ 2e

− 3n

2
8θ2

6ρ(Cj+1,k0 )+2

2
8θ2

nη
2
64θ2

− 3·2

≤ 2e

22j η 2
8θ22

o

.

(53)

For Cj+1,k0 ∈ Λ+ , a similar argument to e1 gives rise to





q
−(j+1)
0
0
0
kcj+1,k − b
cj+1,k k + θ2 2
kProjVj+1,k0 − ProjVb
k
P
ρb(Cj+1,k ) > amax η/2
j+1,k0


Cj+1,k0 ∈Λ+


q
X
−(j+1)
≤
P
kcj+1,k0 − b
cj+1,k0 k + θ2 2
kProjVj+1,k0 − ProjVb
k
ρb(Cj+1,k0 ) ≥ η/2
0


X

j+1,k

Cj+1,k0 ∈Λ+

≤ γ̃3 e−γ̃4 2

2j nη 2

(54)

where γ̃3 := γ̃3 (θ2 , θ3 , amax , d) and γ̃4 := γ̃4 (θ2 , θ3 , θ4 , amax , d).
Finally combining (49), (50), (51), (52), (53) and (54) yields
n
o
n
o
b j,k ≤ η and ∆
¯ j,k ≥ (amax + 2)η ≤ P |∆
b j,k − ∆
¯ j,k | ≥ (amax + 1)η
P ∆
≤ P{e1 > η} + P{e2 > amax η} ≤ γ̃5 e−γ̃6 2

2j nη 2

for some constants γ̃5 := γ̃5 (θ2 , θ3 , a
max , d) and γ̃6 := γ̃6 (θ2 , θ3 , θ4 , amax , d).
¯ j,k ≤ (amax + 2)η and ∆j,k ≥ (2amax + 5)η
Step Three: The probability P ∆
mated as follows. For a fixed Cj,k , we define the function
f (x) = k(Pj − Pj+1 ) 1j,k xk , x ∈ M.
Observe
that |f (x)| ≤ 23 θ2 2−j for any x ∈ M. We define kf k2 =
P
n
1
2
i=1 f (xi ). Then
n

¯ j,k ≤ (amax + 2)η and ∆j,k ≥ (2amax + 5)η
P ∆

´

Mf

2 (x)dρ

2j

(55)

is esti-

and kf k2n =

2


− 2 nη
¯ j,k ≥ η = P {kf k − 2kf kn ≥ η} ≤ 3e 648θ22
≤ P ∆j,k − 2∆

(56)

where the last inequality follows from Györfi et al. (2002, Theorem 11.2). Combining (47),
(55) and (56) yields (24).
Next we turn to the bound in Eq. (24), which corresponds to the case that ∆j,k ≤ η
p
b j,k ≥ bη. In this case we have ∆
b j,k ≤ 3 θ2 2−j ρb(Cj,k ) which implies
and ∆
2
ρb(Cj,k ) ≥

4b2 22j η 2
,
9θ22
52

(57)

Adaptive Geometric Multiscale Approximations for Intrinsically Low-dimensional Data

instead of (48). We shall use the fact that ρ(Cj,k ) ≥ (2b2 22j η 2 )/(9θ22 ) given (57) with high
probability, by writing
n
o
n
o
b j,k ≥ bη ≤P ∆j,k ≤ η and ∆
b j,k ≥ bη ρ(Cj,k ) ≥ 2b2 22j2 η2
P ∆j,k ≤ η and ∆
9θ2
n
o
2b2 22j η 2
4b2 22j η 2
+ P ρ(Cj,k ) ≤ 9θ2
and ρb(Cj,k ) ≥ 9θ2
(58)
2

2

where the first term is estimated as above and the second one is estimated through Eq. (35)
in Lemma 29:
n
o
2b2 22j η 2
4b2 22j η 2
P ρ(Cj,k ) ≤ 9θ2
and ρb(Cj,k ) ≥ 9θ2
2
2


2
2j
2
2b2 22j η 2
2b 2 η
and
|b
ρ
(C
)
−
ρ(C
)|
≥
≤P ρ(Cj,k ) ≤
j,k
j,k
9θ22
9θ22



2b2 22j η 2
2b2 22j η 2
3b2 22j nη 2
≤2e

− 3n(

)2

2
9θ2

6ρ(Cj,k )+2

2
9θ2

−

≤ 2e

2
36θ2

.

Using the estimate in (58), we obtain the bound (24) which concludes the proof.
D.5 Proof of Lemma 20

b τn ∨ Λbτ = Λ
b τn ∧ Λτ /b , or
We will show how Lemma 19 implies Eq. (25). Clearly e2 = 0 if Λ
n
n
b
b
b
b
equivalently Tτn ∪Tbτn = Tτn ∩Tτn /b . In the case e2 > 0, the inclusion Tτn ∩Tτn /b ⊂ Tτn ∪Tbτn is
strict, i.e. there exists Cj,k ∈ T n such that either Cj,k ∈ Tbτn and Cj,k ∈
/ Tτn /b , or Cj,k ∈ Tbτn
and Cj,k ∈
/ Tbτn . In other words, there exists Cj,k ∈ T n such that either ∆j,k < 2−j τn /b and
b j,k ≥ 2−j τn , or ∆j,k ≥ b2−j τn and ∆
b j,k < 2−j τn . As a result,
∆
n
o
X
b j,k < 2−j τn and ∆j,k ≥ b2−j τn
P{e2 > 0} ≤
P ∆
(59)
Cj,k ∈T n

+

X

Cj,k ∈T n

P{e4 > 0} ≤

X

Cj,k ∈T n

n
o
b j,k ≥ 2−j τn .
P ∆j,k < 2−j τn /b and ∆

n
o
b j,k ≥ 2−j τn .
P ∆j,k < 2−j τn /b and ∆

(60)

We apply (24) in Lemma 19 to estimate the first term in (59):
n
o
X
X
2j −2j 2 log n
b j,k < 2−j τn and ∆j,k ≥ b2−j τn ≤
P ∆
α1 e−α2 n2 ·2 κ n
Cj,k ∈T n

= α1 #T n n−α2

κ2

≤ α1 amin nn−α2

κ2

Cj,k ∈T n
2

2 −1)

≤ α1 amin n1−α2 κ = α1 amin n−(α2 κ

.

Using (24), we estimate the second term in (59) and (60) as follows
n
o
X
X
2j 2−2j 2 log n
2 2
b j,k ≥ 2−j τn ≤
P ∆j,k < 2−j τn /b and ∆
α1 e−α2 n2 · b2 κ n ≤ α1 amin n−(α2 κ /b −1) .
Cj,k ∈T n

Cj,k ∈T n

We therefore obtain (25) by choosing κ > κ2 with α2 κ22 /b2 = ν + 1.
53

Liao and Maggioni

Appendix E. Proofs in orthogonal GMRA
E.1 Performance analysis of orthogonal GMRA
The proofs of Theorem 22 and Theorem 24 are resemblant to the proofs of Theorem 4 and
Theorem 8. The main difference lies in the variance term, which results in the extra log
factors in the convergence rate of orthogonal GMRA. Let Λ be the partition associated with
a finite proper subtree T̃ of the data master tree T n , and let
X
X
Sbj,k 1j,k .
SΛ =
Sj,k 1j,k and SbΛ =
Cj,k ∈Λ

Cj,k ∈Λ

Lemma 32 Let Λ be the partition associated with a finite proper subtree T̃ of the data
master tree T n . Suppose Λ contains #j Λ cells at scale j. Then for any η > 0,


βnη 2
X
− 2P
d
j 4 2−2j #j Λ
j≥j
min
P{kSΛ X − SbΛ Xk ≥ η} ≤ αd 
j#j Λ e
(61)
j≥jmin

where α and β are the constants in Lemma 14.

Proof [Proof of Lemma 32] The increasing subspaces {Sj,x } in the construction of orthogonal GMRA may be written as
S0,x = V0,x
⊥
S1,x = V0,x ⊕ V0,x
V1,x

⊥
⊥ ⊥
S2,x = V0,x ⊕ V0,x
V1,x ⊕ V1,x
V0,x V2,x

···

⊥
⊥
⊥ ⊥
Sj,x = V0,x ⊕ V0,x
V1,x ⊕ . . . ⊕ Vj−1,x
· · · V1,x
V0,x Vj,x .

Therefore kProjSj,x − ProjSbj,x k ≤

Pj

`=0 (j

+ 1 − `)kProjV`,x − ProjVb`,x k, and then

j
o
n
o X
n
P kProjSj,x − ProjSbj,x k ≥ t ≤
P kProjV`,x − ProjVb`,x k ≥ t/j 2 .

(62)

`=0

The rest of the proof is almost the same as the proof of Lemma 14 in appendix C with a
slight modification of (41) substituted by (62).
The corollary of Lemma 32 with Λ = Λj results in the following estimate of the variance
in empirical orthogonal GMRA.
Lemma 33 For any η ≥ 0,
2j

P{kSj X − Sbj Xk ≥ η} ≤ αdj#Λj e

EkSj X − Sbj Xk2 ≤

d2 j 4 #Λ

nη
− β2
2 4

2

d j #Λj

log[αdj#Λj ]
.
β22j n

54

j

,

(63)
(64)

Adaptive Geometric Multiscale Approximations for Intrinsically Low-dimensional Data

Proof [Proof of Theorem 22]
EkX − Sbj Xk2 ≤ kX − Sj Xk2 + EkSj X − Sbj Xk2
≤ |ρ|2Aos 2−2sj +

d2 j 4 #Λj log[αdj#Λj ]
d2 j 4 2j(d−2)
αdj2jd
−2sj
2
2
+
.
≤
|ρ|
log
o
As
β22j n
θ1 βn
θ1

 1
∗
When d ≥ 2, We choose j ∗ such that 2−j = µ (log5 n)/n 2s+d−2 . By grouping Λj ∗ into
5
light and heavy cells whose measure is below or above 28
3 (ν + 1) log n/n, we can show that
2s

the error on light cells is upper bounded by C((log5 n)/n) 2s+d−2 and all heavy cells have at
least d points with high probability.
∗

Lemma 34 Suppose j ∗ is chosen such that 2−j = µ
k(X − Pj ∗ X)1

5

n
P ∀ Cj ∗ ,k : ρ(Cj ∗ ,k ) >

k2 ≤

28(ν+1) log n
}
3n
28(ν+1) log5 n
, Cj ∗ ,k
3n

{Cj ∗ ,k :ρ(Cj ∗ ,k )≤



log5 n
n



1
2s+d−2

28(ν+1)θ22 µ2−d
3θ1

with some µ > 0. Then


log5 n
n



2s
2s+d−2

,

o
has at least d points ≥ 1 − n−ν .

Proof of Lemma 34 is omitted since it is the same as the proof of Lemma 16. Lemma 34
guarantees that a sufficient amount of cells at scale j ∗ has at least d points. The probability
estimate in (28) follows from

 5  s 
 5 − d
2s+d−2
2s+d−2 −βθ1 µd−2 C 2 (2s+d−2)4 /d2 log n
b
1
P kSj ∗ X − Sj ∗ Xk ≥ C1 logn n
≤ C2 log n logn n
e
d

d−2 C 2 (2s+d−2)4 /d2
1

≤ C2 (log n) n 2s+d−2 n−βθ1 µ

≤ C2 n1−βθ1 µ

d−2 C 2 (2s+d−2)4 /d2
1

≤ C2 n−ν

provided C1 is chosen such that βθ1 µd−2 C12 (2s + d − 2)4 /d2 − 1 > ν. The proof when d = 1
is completely analogous to that of Theorem 4.

E.2 Performance analysis of adaptive orthogonal GMRA
Proof [Proof of Theorem 24] Empirical adaptive orthogonal GMRA is given by SbΛb o =
τn
P
b
b o Sj,k 1j,k . Using triangle inequality, we have
Cj,k ∈Λ
τn

kX − SbΛb o Xk ≤ e1 + e2 + e3 + e4
τn

with each term given by

e3 := kSΛb

e1 := kX − SΛb

o ∧Λτ o /b
τn
n

X − SbΛb

o ∨Λbτ o
τn
n

o ∧Λτ o /b
τn
n

e2 := kSΛb

Xk

o ∨Λbτ o
τn
n

e4 := kSbΛb

Xk

X − SΛb

o ∧Λτ o /b
τn
n

o ∧Λτ o /b
τn
n

X − SbΛb o Xk

Xk

τn

where b = 2amax + 5. We will prove the case d ≥ 3. Here one proceeds in the same way as
in the proof of Theorem 8. A slight difference lies in the estimates of e3 , e2 and e4 .
55

Liao and Maggioni

Term e3 : Ee23 is the variance. One can verify that T(ρ,τno /b) ⊂ Tj0 := ∪j≤j0 Λj where j0 is the
p
largest integer satisfying 2j0 d ≤ 9b2 θ0 θ22 /(4τno 2 ). The reason is that ∆oj,k ≤ 23 θ2 2−j θ0 2−jd
so ∆oj,k ≥ 2−j τno /b implies 2j0 d ≤ 9b2 θ0 θ22 /(4τno 2 ). For any η > 0,
−

P{e3 > η} ≤ αdj0 #Tτno /b e

4
j0

P

βnη 2
−2j # T o
j τn /b
j≥jmin 2

−

≤ αdj0 #Tτno /b e

βnη 2
4 |ρ|p (τ o /b)−p
j0
o n
Bs

The estimate of Ee23 follows from
ˆ
Ee23 =
≤

0

ˆ

+∞

ηP {e3 > η} dη =

j04 log αj0 #Tτno /b
βn

X

j≥jmin

2

−2j

0

+∞

η min 1, αdj0 #Tτno /b e

#j Tτno /b ≤

5
C logn n (τno /b)−p

−

βnη 2
4P
−2j # T o
j0
j τn /b
j≥jmin 2

!

≤ C(θ0 , θ2 , θ3 , amax , κ, d, s)

dη


log5 n
n



2s
2s+d−2

Term e2 and e4 : These two terms are analyzed with Lemma 35 stated below such that
for any fixed but arbitrary ν > 0,
P{e2 > 0} + P{e4 > 0} ≤ β1 amin /dn−ν
if κ is chosen such that κ > κ2 with d4 β2 κ22 /b2 = ν + 1.
Lemma 35 b = 2amax + 5. For any η > 0 and any Cj,k ∈ T
o
o n
 n
2j 2 4
b o ≥ bη
b o ≤ η and ∆o ≥ bη , P ∆o ≤ η and ∆
≤ β1 je−β2 n2 η /j ,
max P ∆
j,k
j,k
j,k
j,k
with positive constants β1 := β1 (θ2 , θ3 , θ4 , amax , d) and β2 := β2 (θ2 , θ3 , θ4 , amax , d).

References
M. Aharon, M. Elad, and A. Bruckstein. K-SVD: Design of dictionaries for sparse representation. In PROCEEDINGS OF SPARS 05’, pages 9–12, 2005.
W. K. Allard, G. Chen, and M. Maggioni. Multi-scale geometric methods for data sets II:
Geometric multi-resolution analysis. Applied and Computational Harmonic Analysis, 32
(3):435–462, 2012.
M. Belkin and P. Niyogi. Laplacian eigenmaps for dimensionality reduction and data representation. Neural Computation, 15(6):1373–1396, 2003.
A. Beygelzimer, S. Kakade, and J. Langford. Cover trees for nearest neighbor. In International Conference on Machine Learning, pages 97–104, 2006.
P. Binev, A. Cohen, W. Dahmen, R.A. DeVore, and V. Temlyakov. Universal algorithms for
learning theory part i: piecewise constant functions. J. Mach. Learn. Res., 6:1297–1321,
2005.
56

.

Adaptive Geometric Multiscale Approximations for Intrinsically Low-dimensional Data

P. Binev, A. Cohen, W. Dahmen, R.A. DeVore, and V. Temlyakov. Universal algorithms for
learning theory part ii: piecewise polynomial functions. Constr. Approx., 26(2):127–152,
2007.
G. Canas, T. Poggio, and L. Rosasco. Learning manifolds with k-means and k-flats. In
Advances in Neural Information Processing Systems, 2012.
E. Candes and T. Tao. The Dantzig selector: statistical estimation when p is much larger
than n. Annals of Statistics, pages 2313–2351, 2007. /0506081.
G. Chen and G. Lerman. Foundations of a multi-way spectral clustering framework for
hybrid linear modeling. Foundations of Computational Mathematics, 9:517–558, 2009.
DOI 10.1007/s10208-009-9043-7.
G. Chen and M. Maggioni. Multiscale geometric and spectral analysis of plane arrangements. In Conference on Computer Vision and Pattern Recognition, 2011.
Guangliang Chen, M. Iwen, Sang Chin, and M. Maggioni. A fast multiscale framework
for data in high-dimensions: Measure estimation, anomaly detection, and compressive
measurements. In Visual Communications and Image Processing (VCIP), 2012 IEEE,
pages 1–6, 2012.
Scott Shaobing Chen, David L. Donoho, and Michael A. Saunders. Atomic decomposition
by basis pursuit. SIAM Journal on Scientific Computing, 20(1):33–61, 1998. doi: 10.
1137/S1064827596304010. URL http://link.aip.org/link/?SCE/20/33/1.
M. Christ. A T (b) theorem with remarks on analytic capacity and the Cauchy integral.
Colloquium Mathematicae, 60/61(2):601–628, 1990.
A. Cohen, I. Daubechies, O. G. Guleryuz, and M. T. Orchard. On the importance of combining wavelet-based nonlinear approximation with coding strategies. IEEE Transactions
on Information Theory, 48(7):1895–1921, 2002.
R. R. Coifman, S. Lafon, A. B. Lee, M. Maggioni, B. Nadler, F. Warner, and S. W. Zucker.
Geometric diffusions as a tool for harmonic analysis and structure definition of data:
Diffusion maps. Proceedings of the National Academy of Sciences of the United States of
America, 102(21):7426–7431, 2005a.
R. R. Coifman, S. Lafon, A. B. Lee, M. Maggioni, B. Nadler, F. Warner, and S. W. Zucker.
Geometric diffusions as a tool for harmonic analysis and structure definition of data:
Multiscale methods. Proceedings of the National Academy of Sciences of the United
States of America, 102(21):7432–7438, 2005b.
I. Daubechies. Ten lectures on wavelets. Society for Industrial and Applied Mathematics,
1992. ISBN 0-89871-274-2.
G. David and S. Semmes. Analysis of and on uniformly rectifiable sets, volume 38 of
Mathematical Surveys and Monographs. American Mathematical Society, Providence,
RI, 1993. ISBN 0-8218-1537-7.
57

Liao and Maggioni

C. Davis and W. M. Kahan. The rotation of eigenvectors by a perturbation. iii. SIAM
Journal on Numerical Analysis, 7(1):1–46, 1970.
D. Deng and Y. Han. Harmonic analysis on spaces of homogeneous type, volume 19. Springer
Science and Business Media, 2008.
D. Donoho. Compressed sensing. IEEE Tran. on Information Theory, 52(4):1289–1306,
April 2006.
D. L Donoho and C. Grimes. Hessian eigenmaps: new locally linear embedding techniques
for high-dimensional data. Proceedings of the National Academy of Sciences of the United
States of America, pages 5591–5596, March 2003.
D. L. Donoho and J. M. Johnstone.
Biometrika, 81(3):425–455, 1994.

Ideal spatial adaptation by wavelet shrinkage.

D. L. Donoho and J. M. Johnstone. Adapting to unknown smoothness via wavelet shrinkage.
Journal of the American Statistical Association, 90(432):1200–1224, 1995.
A. Eftekhari and M. B. Wakin. Sparse subspace clustering. Applied and Computational
Harmonic Analysis, 39(1):67–109, 2015.
E. Elhamifar and R. Vidal. Sparse subspace clustering. In Conference on Computer Vision
and Pattern Recognition, pages 2790–2797, June 2009.
R. Fergus F. Li and P. Perona. One-shot learning of object categories. IEEE transactions
on pattern analysis and machine intelligence, 28(4):594–611, 2006.
M. Fischler and R. Bolles. Random sample consensus: A paradigm for model fitting with
applications to image analysis and automated cartography. Comm. of the ACM, 24(6):
381–395, 1981.
S. Gerber and M. Maggioni. Multiscale dictionaries, transforms, and learning in highdimensions. In Proc. SPIE conference Optics and Photonics, 2013.
Rémi Gribonval, Rodolphe Jenatton, Francis Bach, Martin Kleinsteuber, and Matthias
Seibert. Sample complexity of dictionary learning and other matrix factorizations. 2013.
L. Györfi, M. Kohler, A. Krzyżak, and H. Walk. A distribution-free theory of nonparametric
regression. New York: Springer, 2002.
J. Ho, M. Yang, J. Lim, K. Lee, and D. Kriegman. Clustering appearances of objects
under varying illumination conditions. In Conference on Computer Vision and Pattern
Recognition, volume 1, pages 11–18, 2003.
H. Hotelling. Analysis of a complex of statistical variables into principal components.
Journal of Educational Psychology, 24(4):17–441,498–520, 1933.
H. Hotelling. Relations between two sets of variates. Biometrika, 28(3/4):321–377, 1936.
58

Adaptive Geometric Multiscale Approximations for Intrinsically Low-dimensional Data

M. A. Iwen and M. Maggioni. Approximation of points on low-dimensional manifolds via random linear projections. Inference and Information, 2(1):1–31, 2013.
arXiv:1204.3337v1, 2012.
P. W. Jones. Rectifiable sets and the traveling salesman problem. Inventiones Mathematicae,
102(1):1–15, 1990.
G. Karypis and V. Kumar. A fast and high quality multilevel scheme for partitioning
irregular graphs. SIAM Journal on Scientific Computing, 20(1):359–392, 1999.
K. Kreutz-Delgado, J. F. Murray, B. D. Rao, Kjersti Engan, T.-W. Lee, and T. J. Sejnowski.
Dictionary learning algorithms for sparse representation. Neural Comput., 15(2):349–396,
February 2003.
A. Krizhevsky and G. Hinton. Learning multiple layers of features from tiny images. 2009.
M.S. Lewicki, T.J. Sejnowski, and H. Hughes. Learning overcomplete representations. Neural Computation, 12:337–365, 1998.
A. V. Little, M. Maggioni, and L. Rosasco. Multiscale geometric methods for data
sets I: Multiscale SVD, noise and curvature. Technical report, MIT-CSAIL-TR-2012029/CBCL-310, MIT, Cambridge, MA, September 2012. URL http://dspace.mit.
edu/handle/1721.1/72597.
A.V. Little, Y.-M. Jung, and M. Maggioni. Multiscale estimation of intrinsic dimensionality
of data sets. In Proc. A.A.A.I., 2009a.
A.V. Little, J. Lee, Y.-M. Jung, and M. Maggioni. Estimation of intrinsic dimensionality of
samples from noisy low-dimensional manifolds in high dimensions with multiscale SV D.
In Proc. S.S.P., 2009b.
G. Liu, Z. Lin, and Y. Yu. Robust subspace segmentation by low-rank representation.
Proceedings of the 27th International Conference on Machine Learning, 2010.
Y. Ma, H. Derksen, W. Hong, and J. Wright. Segmentation of multivariate mixed data via
lossy data coding and compression. IEEE Transactions on Pattern Analysis and Machine
Intelligence, 29(9):1–17, 2007.
Y. Ma, A. Y. Yang, H. Derksen, and R. Fossum. Estimation of subspace arrangements
with applications in modeling and segmenting mixed data. SIAM Review, 50(3):413–458,
2008.
Mauro Maggioni, Stanislav Minsker, and Nate Strawn. Multiscale dictionary learning: Nonasymptotic bounds and robustness. J. Mach. Learn. Res., 17(1):43–93, January 2016.
ISSN 1532-4435. URL http://dl.acm.org/citation.cfm?id=2946645.2946647.
J. Mairal, F. Bach, J. Ponce, and G. Sapiro. Online learning for matrix factorization and
sparse coding. Journ. Mach. Learn. Res., 11:19–60, 2010.
S. G. Mallat. A wavelet tour in signal processing. Academic Press, 1998.
59

Liao and Maggioni

Andreas Maurer and Massimiliano Pontil. K-dimensional coding schemes in hilbert spaces.
IEEE Transactions on Information Theory, 56(11):5839–5846, 2010.
K. Pearson. On lines and planes of closest fit to systems of points in space. Philosophical
Magazine, 2(11):559–572, 1901.
G. Peyré. Manifold models for signals and images. Computer Vision and Image Understanding, 113(2):249–260, 2009.
G. Peyré. Sparse modeling of textures. Journal of Mathematical Imaging and Vision, 34
(1):17–31, 2009.
M. Protter and M. Elad. Sparse and redundant representations and motion-estimation-free
algorithm for video denoising, 2007.
M. A. Rohrdanz, W. Zheng, M. Maggioni, and C. Clementi. Determination of reaction
coordinates via locally scaled diffusion map. J. Chem. Phys., (134):124116, 2011.
S. T Roweis and LK Saul. Nonlinear dimensionality reduction by locally linear embedding.
Science, 290:2323–2326, 2000.
A. Szlam. Asymptotic regularity of subdivisions of euclidean domains by iterated PCA and
iterated 2-means. Applied and Computational Harmonic Analysis, 27(3):342–350, 2009.
J. B. Tenenbaum, V. De Silva, and J. C. Langford. A global geometric framework for
nonlinear dimensionality reduction. Science, 290(5500):2319–2323, 2000.
J. A. Tropp. User-friendly tools for random matrices: An introduction. NIPS version, 2014.
URL http://users.cms.caltech.edu/~jtropp/.
Daniel Vainsencher, Shie Mannor, and Alfred M Bruckstein. The sample complexity of
dictionary learning. The Journal of Machine Learning Research, 12:3259–3281, 2011.
R. Vidal, Y. Ma, and S. Sastry. Generalized principal component analysis (GPCA). IEEE
Transactions on Pattern Analysis and Machine Intelligence, 27(12), 2005.
J. Yan and M. Pollefeys. A general framework for motion segmentation: Independent,
articulated, rigid, non-rigid, degenerate and nondegenerate. In European conference on
computer vision, volume 4, pages 94–106, 2006.
T. Zhang, A. Szlam, Y. Wang, and G. Lerman. Randomized hybrid linear modeling by
local best-fit flats. In CVPR, San Francisco, CA, June 2010.
Z. Zhang and H. Zha. Principal manifolds and nonlinear dimensionality reduction via
tangent space alignment. Journal of Shanghai University, 8(4):406–424, 2004.
W. Zheng, M. A. Rohrdanz, M. Maggioni, and C. Clementi. Polymer reversal rate calculated
via locally scaled diffusion map. J. Chem. Phys., (134):144108, 2011.
L. Zwald and G. Blanchard. On the convergence of eigenspaces in kernel principal component analysis. In Advances in Neural Information Processing Systems, pages 1649–1656,
2006.
60

