Supervised Learning of Labeled Pointcloud Differences via Cover-Tree Entropy Reduction
Abraham Smith Paul Bendich John Harer Alex Pieloch Jay Hineman





CDER]Supervised Learning of Labeled Pointcloud Differences via Cover-Tree Entropy Reduction




A. Smith

Mathematics, Statistics, and Computer Science

University of Wisconsin-Stout

Menomonie, WI, 54751, USA

Geometric Data Analytics, Inc. 

Durham, NC, 27707, USA
smithabr@uwstout.edu
P. Bendich 

Mathematics 

Duke University 

Durham, NC, 27708, USA 

Geometric Data Analytics, Inc. 

Durham, NC, 27707, USA 

United States
bendich@math.duke.edu
J. Harer 

Mathematics, Electrical and Computer Engineering

Duke University

Durham, NC, 27707, USA 

Geometric Data Analytics, Inc. 

Durham, NC, 27707, USA
harer@math.duke.edu
A. Pieloch

Mathematics

Columbia University

New York, NY, 10027, USA
pieloch@math.columbia.edu
J. Hineman 

Geometric Data Analytics, Inc.

Durham, North Carolina, 27707, USA
jay.hineman@geomdata.com





We introduce a new algorithm, called CDER, for supervised machine learning
that merges the multi-scale geometric properties of Cover Trees with the
information-theoretic properties of entropy.  CDER applies to a training set of 
labeled pointclouds embedded in a common Euclidean space.  If typical pointclouds corresponding to distinct labels tend to differ
at any scale in any sub-region, CDER can identify these differences in linear time, creating a set of distributional
coordinates which act as a feature extraction mechanism for supervised learning. We describe theoretical properties and implementation details
of CDER, and illustrate its benefits on several synthetic examples.

MSC 62H30 MSC 60G55


[
    Abraham Smith Paul Bendich John Harer Alex Pieloch Jay Hineman
    October 24, 2016
==================================================================






§ ACKNOWLEDGMENTS


All authors partially supported by OSD contract number N00024-13-D6400 via
JHU-APL subcontract number 131753.  (Distribution Statement A - Approved for
public release; distribution is unlimited.) Bendich and Harer were also
partially supported by NSF award BIGDATA 1444791. 
Harer was also partially supported by the DARPA MoDyL program, number HR0011-16-2-0033.
Pieloch was supported by the National Science Foundation Graduate Student Fellowship Program through grant DGE 16-44869.
We thank Christopher J.
Tralie for many helpful insights about cover trees, and we thank David Porter
and Michael Williams at JHU-APL for many motivational discussions.

Our code is available at <https://github.com/geomdata/gda-public/> with
documentation at <https://geomdata.github.io/gda-public/>



§ OVERVIEW AND ASSUMPTIONS

We propose a new supervised machine-learning method for classification, where each object to be
classified is a weighted pointcloud, and the classification task is to learn which of a finite set of labels should be applied
to the pointcloud. The method is fast, data-driven,
multi-scale, and requires no tuning.  Additionally, its details are transparently
geometric; it does not suffer from the “black box” interpretation
difficulties that arise in many machine-learning methods.  We call it
Cover-tree Differencing via Entropy Reduction [CDER, pronounced “cedar”].

A high-level sketch of CDER is as follows. We construct a partial cover tree
<cit.> on the union of all labeled pointclouds that are given in some
training set, and we search the cover tree for convex regions that are likely
to be local minima of entropy <cit.>.  For each such region, we build
distributional coordinates from the dominant labels' data.  An ordered sequence
of these distributional coordinates allows us to determine the likelihood of
each label for an unlabeled test pointcloud.

Section <ref> explains the notion of distributional coordinates
defined on pointclouds as a front-end for supervised learning.
Section <ref> details the cover-tree algorithm and our
enhancements to it.  Section <ref> explains our simple approach to
entropy minimization in the labeled pointcloud context. Section <ref> gives a formal
description of CDER. Section <ref> illustrates CDER with several
synthetic examples.  The remainder of this section establishes context and
notation for this problem.


Fix a Euclidean[See the Discussion section for further thoughts on this requirement.] space ℝ^D.  A
pointcloud is a finite[In this document, |X| indicates the
cardinality of a finite set X, whereas v indicates a norm.] set X = {
x_1, x_2, …, x_|X|}⊂ℝ^D.  We allow weighted
pointclouds, where each x ∈ X has a positive weight w(x); the utility of
weights is explained in Section <ref>.  A cloud
collection is a finite set of pointclouds 
    𝒳 = { X_1, X_2, …,
    X_N},
 with N=|𝒳|.  Let 𝒳=
⋃_i=1^N  X_i ⊂ M denote the union of all pointclouds in a cloud
collection, which is itself a pointcloud with weights and labels of each point
x ∈𝒳 inherited from the respective X_i  x.

For our supervised learning task, the training data is a cloud collection
𝒳, where each pointcloud is labeled from a finite set of labels,
Λ.  For the sake of discussion and visualization, we usually interpret
Λ as a set of colors, but for mathematical purposes, Λ = { 1,
…, L} where L = |Λ|.  Let λ: 𝒳→Λ
denote the label function, and we also use λ to indicate a generic
label in Λ.

It is acceptable that the pointclouds X_i have unequal sizes across X_i ∈𝒳.  It is also acceptable that the labeled sub-collections
𝒳|_λ = { X_i  : λ(X_i) = λ} =
λ^-1(λ) have unequal sizes across labels λ∈Λ.
Pointwise weights can be assigned to compensate for these sizes, as in Section <ref>.
Our only structural hypothesis is that, for each labeled sub-collection, the
pointclouds are sampled from an underlying density function—or from several
density functions as chosen by a random process—on ℝ^D.  For each
label λ, let f_λ denote this density function.

We aim to efficiently find regions in ℝ^D where one label has greater density than the
other labels.  That is, we seek convex regions Ω⊂ℝ^D that are
“characteristic” of a particular label λ because f_λ is
unusually prominent there, while ignoring regions where the various density
functions are similar.  See Figure <ref>.



We analyze the relative density of labels over regions Ω using the
information-theoretic notion of entropy (see Section <ref>), locating subregions
of low entropy.
Because of the use of cover trees, our method emphasizes regions Ω that
are very large, very dense, or small-yet-remote (see
Section <ref>).   For each of these regions Ω, we
construct a function g that approximates f_λ near Ω; such a
g is an example of what we call a distributional coordinate (see
Section <ref>).



§ WEIGHTS

Typically, pointclouds are experimental samples that do not come with
pre-assigned weights, and we assign weights as follows: We have the
prior assumption that each label λ∈Λ is equally likely
among the test data.[See the discussion section for more general thoughts about this assumption.]
We give the training set total weight 1.  That is, each label
is alloted total weight 1/L, where L = |Λ|.  Moreover, we make the
assumption that each sample X_i is equally representative of its underlying
density function f_λ_i, regardless of |X_i|.  That is, each
pointcloud X_i is alloted total weight 1/L N_i, where N_i is the
cardinality of X|_λ_i, the number of training pointclouds
with label λ_i.  Finally, we assume all points x ∈ X_i have equal
weight, so w(x) = 1/L N_i |X_i|.  Let 𝒳
denote the labeled and weighted pointcloud given by the union of all X_i with
these labels and weights.

For example, suppose that the training set consists of X_1, X_2, X_3 where
|X_1| = 4, |X_2|=10, and |X_3|=3 with labels red, blue, and red,
respectively.  Suppose we want 𝒳 to have unit weight.
Then each x ∈ X_1 has weight w(x) = 1/2·2·4 =
1/16, each x ∈ X_2 has weight w(x)= 1/2·1·10 =
1/20, and each x ∈ X_3 has weight w(x) = 1/2·2·3=
1/12.

These weights are chosen to model relative likelihoods in the following way:
If each label λ has an underlying density function f_λ, and
each X_i with label λ is a sample from f_λ, then
∫_Ω f_λ should be proportional to
∑_X_i ∈𝒳|_λ∑_x ∈ X_i∩Ω w(x)
for any measurable set Ω.




§ DISTRIBUTIONAL COORDINATES


Given a cloud collection 𝒳 as above, we would like to map each
pointcloud X to a vector v(X) ∈ℝ^k for some fixed
dimension k. While this is an interesting topic in its own right, we are
mainly concerned with viewing this as a feature extraction method: that
is,
the vectors v_i = v(X_i) ∈ℝ^k and corresponding
labels λ_i = λ(X_i) will be used as input for any supervised
learning algorithm that requires collections of vectors in a common Euclidean
space, along with class labels.

One could imagine a supervised learning approach to labeled pointclouds that
does not perform feature extraction; for example, by first defining a
(non-Euclidean) metric between pointclouds, and then applying a
nearest-neighbor or clustering method.  However, most such metrics are unstable
to noise (such as Hausdorff distance) or are slower-than-linear to compute
(such as various Kullback-Leibler Divergence <cit.> or the
Wasserstein metric <cit.>).

This section frames the problem of feature extraction in the common language of
distributional coordinates. We start with a definition of this intuitive
concept, and then describe the advantages and disadvantages of some standard
examples. The approach we advocate comes in the next section.

Let g be an integrable function on Euclidean space.  For any
pointcloud X = {x_1, …, x_|X|}, with pointwise weights w_i = w(x_i),
define 
    ∫_X g = w_1 g(x_1) + … + w_|X|g(x_|X|).
 We
refer to g as a distributional coordinate.  Any ordered set (g_1,
…, g_k) of distributional coordinates transforms X into a vector
𝐯_X = (∫_X g_1, …, ∫_X g_k) ∈ℝ^k.
It is computationally important to keep the dimension k—or at least the intrinsic
dimension of {𝐯_X :  X ∈𝒳}—as small as possible,
as we hope to minimize the negative effect of the curse of dimensionality on the supervised learning
problem.

We use the cloud collection 𝒳 shown on the bottom of Figure <ref> as a
notional working example in the Euclidean plane. There are two class labels,
“magenta” and “green.” We imagine that each magenta pointcloud (a typical example on the top-left of the figure) contains a large
number of points sampled from a large central blob and a few points sampled
from much smaller blobs along the horizontal axis. Each green pointcloud (a typical example on the top-right of the figure) contains a
large number of points sampled from the same large central blob and a few
points sampled from smaller, different blobs along the horizontal axis.

In a typical supervised learning context, 𝒳 would be split (say 80:20) into
training and testing sets, and we would build a set of distributional
coordinates from the labeled training set which could then infer  the correct
color labels on the testing set.



 §.§ Binning

A simple approach would be to let each g_i be the indicator function on some
region R_i; that is, g_i takes the value 1 on some region R_i and 0
elsewhere. Then ∫_X g_i counts the total weight of the points X ∩ R_i. Assuming the regions R_1, … R_k cover the domain with overlap only on boundaries (for example, the R_i could be the squares shown in Figure <ref>)
this is what is commonly referred to as binning, and it has been used (for example, <cit.>) successfully for supervised-learning on
labeled pointclouds. On the other hand, binning is obviously unstable: consider a point x ∈ X that lies very close to the boundary between two adjacent bins.
It also can suffer from the curse of dimensionality: a larger number of bins than are shown in Figure <ref> would be needed to distinguish the class labels.




 §.§ Gaussian Mixtures

To address instability, one could define each g_i to be a Gaussian with mean
μ_i and covariance matrix Σ_i. The obvious questions are: how many
Gaussians (modes) should be used, what are their means, and what are their
covariances?  A data-driven approach should be used to fit a Gaussian Mixture Model
(<cit.>, VI.8) to 𝒳.

It is apparent to the human eye that nine Gaussians (1 large green and magenta, 4
small magenta, and 4 small green) model the example in Figure <ref>.  That
is, we could “by hand” define a good set of distributional coordinates to establish a feature vector
𝐯 = (g_1, …, g_9) : 𝒳→ℝ^9.  Moreover,
for the purposes of classification, the central “large” distributional
coordinate g_1 is useless.  And perhaps the two innermost “small”
Gaussians are less useful, as the remaining six distributional coordinates
are sufficient for accurate classification.  Our goal is to replicate this
observation: we construct Gaussian modes (g_1, …, g_k) in a
data-driven manner from 𝒳 that are appropriate for
use as labeled distributional coordinates for the label classification
problem.  We now proceed to the computational details.




§ COVER TREES

Cover trees were originally conceived <cit.> as a way to accelerate
the nearest-neighbors problem.  More recently, they have been used in
dictionary learning <cit.> and for speeding up the analysis of the
topology of pointclouds <cit.>, to pick just a few applications. This section gives basic
definitions for cover trees as well as a fast algorithm for their construction;
our account is far more detailed than usually appears in the literature and is
meant to accompany our publicly available code. Sections <ref> and <ref>
use cover trees to define a set of label-driven distributional coordinates on
a cloud collection with class labels.

For any finite pointcloud[For the moment we forget that
𝒳 is a union of pointclouds and simply treat it as a
single pointcloud to be studied.] 𝒳⊂ M, consider a filtration

    A_0 ⊂⋯⊂ A_ℓ⊂ A_ℓ+1⊂⋯⊂𝒳.

An element a ∈ A_ℓ is called
an adult[To reduce confusion (hopefully) by maintaining a
consistent metaphor, we diverge from
prior works' notation in two places.  Our adult is called a
center elsewhere.  Our guardian is usually called a
parent.  We avoid the word parent because there are many distinct tree
relations embedded in this algorithm.] at level ℓ.  The set A_ℓ∖ A_ℓ-1 is the
cohort at ℓ.

A cover tree builds a filtration by covering 𝒳 with balls of smaller and
smaller radius centered at the points in A_ℓ.
See Figure <ref>.



Specifically, a cover tree is a filtration of 𝒳 with the following
additional properties:


  * |A_0|=1.  Let A_0 = {a_0};

  * There is a radius r_0 such that 𝒳⊂ B(a_0, r_0);

  * There is a real number 0< θ < 1, called the shrinkage
ratio such that, for every ℓ,
𝒳⊂∪_a_i ∈ A_ℓ B(a_i, r_ℓ) where r_ℓ =
θ^ℓ r_0.

  * For each ℓ, if a_i, a_j ∈ A_ℓ, then d(a_i, a_j) > r_ℓ.

  * For each ℓ, each point x ∈𝒳 is assigned to a
guardian a_i ∈ A_ℓ such that x lies in the ball B(a_i, r_ℓ); we say x is
a child of a_i.  Each a_i ∈ A_ℓ is its own guardian and its own
child.

  * There is a tree structure on the (level, adult) pairs of the filtration
(ℓ, a_i), where there is a tree relation (ℓ, a_i) → (ℓ+1, a_k) if
a_k was a child of a_i at level ℓ.  We say a_k is a
successor of a_i, and a_i is a predecessor of a_k.
Note that (ℓ, a_i) → (ℓ+1, a_i) for all a_i ∈ A_ℓ.




 §.§ Vocabulary

Cover trees can be generated using a fast algorithm relying on the notion of
friends, first introduced in <cit.>.  The rough algorithm is
familiar to many practitioners, but is not widely published.  We enhance this
algorithm by giving each point a weight and a label.   We
describe our method here because the details of the construction are crucial to
the supervised-learning task of this article.

Extending the maturation/reproduction metaphor of adults,
children, and guardians above, a child x with guardian a_i at
level ℓ is called a teen if 1/2 r_ℓ < d(a_i, x), and it is
called a youngin if d(a_i, x) ≤1/2 r_ℓ.  See
Figure <ref>.
When changing from level ℓ to level ℓ+1, the radius of each ball shrinks to
r_ℓ+1 = θ r_ℓ.  Children farther than r_ℓ+1 from their guardians
become orphans.  We must decide whether these orphans should be
adopted by other adults at level ℓ+1, or if the orphans should be
emancipated as new adults at level ℓ+1.
That is, the newly emancipated adults at level ℓ+1 comprise the
cohort at level ℓ+1.




Our algorithm uses several notions of friends, set by the following bounds
at level ℓ:

    T_1(ℓ)    = (2 + θ)r_ℓ,
    
    T_2(ℓ)    = (2 + 2θ)r_ℓ,
    
    T_3(ℓ)    = 2( 1 + θ + θ^2 + ⋯ ) r_ℓ =
    2/1-θr_ℓ.

For each k=1,2,3, we say that adults a_i and a_j at level ℓ are
type-k friends if d (a_i, a_j) ≤ T_k(ℓ).

It is easy to verify the following recursion relations:

    T_1(ℓ) ≤ T_2(ℓ) ≤ T_3(ℓ),
    
    T_1(ℓ+1) ≤ T_1(ℓ),
    
    T_2(ℓ+1) ≤ T_2(ℓ),
    
    T_3(ℓ+1) ≤ T_3(ℓ).













Here are the essential facts about friends, which are immediate consequences of
the triangle inequality:



  T1 If x is a child of a_i at level ℓ that is orphaned at level
ℓ+1, then x can be adopted only by either the type-1 friends of a_i at
level ℓ or the newly-emancipated children at level ℓ+1 of those type-1
friends.  See Figure <ref>.


  T2 If we choose to re-assign children to the nearest adult, then teen x
with guardian a_k can be adopted by adult a_h at level ℓ+1, where the
level-ℓ predecessors of a_k and a_h are type-2 friends.
See Figure <ref>.


  T3 If a_k and a_h are type-3 friends at level ℓ+1, then their
predecessors a_i and a_j are type-3 friends at level ℓ.
See Figure <ref>.





Suppose a_k is an adult at level ℓ+1 with predecessor a_i.
If x has guardian a_k at level ℓ+1 and guardian of a_j at level
ℓ,
then a_j is a type-1 friend of a_i, and d(a_k, a_j) ≤ r_ℓ +
r_ℓ+1= 1/2 T_2(ℓ).  See Figure <ref>.


Based on this scenario, we say a_j ∈ A_ℓ is an elder of a_k
∈ A_ℓ+1 if d(a_k, a_j) ≤1/2T_2(ℓ).
Moreover, if a_i = a_k is its own predecessor, then a_j = a_i = a_k is its
only elder, so the notion of elder is nontrivial only if a_k belongs to
cohort ℓ.


For the purpose of computing elders, there are three particularly good choices
of θ which require no additional computation.
If θ = 1/2(√(5)-1), then 1/2
T_2(ℓ) = T_1(ℓ+1).
If θ = 1/2, then 1/2 T_2(ℓ) = T_2(ℓ+1).
If θ = √(2)-1, then 1/2 T_2(ℓ) = T_3(ℓ+1).
In these cases, the elders of a_k are the friends (of the respective type)
a_j that have cohort ℓ or less.  See Figure <ref>.



These facts allow the following algorithm.



 §.§ Algorithm

At level 0, the cover tree has a single adult, a_0.  All points are its
children.  Its only friend is itself.  Every x ∈𝒳 has some weight and some
label, which we assume as given input.  For each label λ, compute the total weight and weighted mean of all children of
a_0 with that label.

Suppose at level ℓ, the friends, children, per-label weights, and per-label means are
known for all adults a_i ∈ A_ℓ.  We compute these for level ℓ+1 as
follows:



  * Advance. Shrink the radius by a factor of θ to r_ℓ+1 =
r_ℓθ.   Every adult a_i ∈ A_ℓ remains an adult a_i ∈
A_ℓ+1.  Mark a_i as its own elder.


  * Orphan. For each a_i ∈ A_ℓ, find its orphaned children: those x
such that d(a_i, x) > r_ℓ+1.  Sort the orphans so that the first orphan
is the one nearest the per-label mean for the label of greatest weight, the
second orphan is the one nearest the per-label mean of the label of
second-greatest weight, and so on, via a two-dimensional argsort without
repetition.  The reason for this sorting will become apparent in
Section <ref>.


  * Adopt or emancipate. For each a_i ∈ A_ℓ, consider its orphans in
order.  For each x orphaned by a_i, check whether x can
be adopted by another adult using the type-1 friends of a_i:


  * If so, the adoptive guardian will be a adult a_h such that d(a_h, x)
≤ r_ℓ+1.  The predecessor a_j of a_h has d(a_j, a_i) ≤ T_1(l).


  * If not, then x is a new adult at level ℓ+1.  Its only child is
itself.  Its friends are initialized to be itself and its predecessor a_i (but
these will be updated in the Befriending stage).  Note the set A_ℓ+1
is ordered by cohort, and within a cohort is ordered by predecessor, and within a
predecessor block is ordered by proximity to the predecessor's per-label means.



  * Exchange Teens.[This step is optional for general cover
tree, but it is particularly useful for the computation of elders, which are
needed in Section <ref>.]
If x has guardian a_i, then check whether it is closer to some other
guardian using type-2 friends.  That is, find a_h = _A_ℓ+1 d(x, ·).
The predecessor a_j of a_h has d(a_j, a_i) ≤ T_2(ℓ).  Note that
a_h ≠ a_i implies d(a_h, x) > 1/2 r_ℓ+1, meaning that only teens (not
youngins) can be exchanged.


  * Befriend: Now, all (ℓ+1)-adults and their children are known.  We need to
find new friends using the recursion relation and the predecessors.

    
  * A pair of ℓ+1 adults, a_k and a_h, are type-3
friends if their predecessors a_i and a_j are type-3 friends at level ℓ
and d(a_k, a_h) ≤ T_3(ℓ+1).

    
  * A pair of ℓ+1 adults, a_k and a_h, are type-2 friends if they are
type-3 friends and d(a_k, a_h) ≤ T_2(ℓ+1).

    
  * A pair of ℓ+1 adults, a_k and a_h, are type-1 friends if they are
type-2 friends and d(a_k, a_h) ≤ T_1(ℓ+1).

    
  * For each newly emancipated adult a_k of cohort ℓ+1, find its
    elders by searching its friends for elements of cohort ℓ or less.


  * Weigh: For each adult a_k at level ℓ+1, compute the per-label
weights and per-label weighted means for its children.
Sort the labels by their weights from greatest to least.

At level ℓ+1, the friends, children, per-label weights, and per-label
means are all known.



 §.§ Algorithmic Efficiency

Suppose n = |𝒳|, and the cover-tree filtration stops at
level ℓ_max.  We are interested in the storage and computational
complexity of this algorithm in terms of n, for a fixed ambient
dimension D.

The desired output is the successor/predecessor tree and the children/guardian
relationships at each level.  Thus, storage of friend lists is the only issue
that might lead us to another algorithm.
The main question is: given an adult a_i at level ℓ, how many other
adults are within T_3(ℓ) = 2/(1-θ)r_ℓ of a_i?

Since no two adults at level ℓ are within r_ℓ of one another, this
is related to the idea of “doubling dimension.”   Define 𝔻(D, α)
to be the maximum number of points in ℝ^D that can be placed in the
ball of radius α, including the center, so that no two points have
distance ≤ 1.  This number grows with the ambient dimension D.  In
Euclidean space, it is independent of scale; that is, the same number works for
a ball of radius rα using points of distance at least r.  The number
can be bounded by the number of simplices with edge length at least 1 that
can be put into a ball of radius α.  Since these simplices have a
minimum volume √(D+1)/D!2^D/2, we get a bound like:

    𝔻≤(B_α)D! 2^D/2/√(D+1)∼D! (α)^D (2 π)^D/2/(D/2)!√(D+1)∼ O(D! 2^D α^D)
    .



Setting α = 2/1-θ, we get an upper bound on the expected
size of the friends list, which is a constant independent of the level ℓ.
Thus, the total storage size is O(n); in particular, it is at most 3
𝔻n, where the 3 accounts for the lists of type-1, type-2, and
type-3 friends.

Setting α = 1/θ, we get an upper bound on the expected
number of successors of an adult from level ℓ to ℓ+1, which is a
constant independent of the level ℓ.  The size of A_0 is |A_0|=1. The
maximum expected size of A_1 is 𝔻.  The maximum expected size of
A_2 is 𝔻^2.  Generally, the maximum expected size of A_ℓ is
𝔻^ℓ.  Moreover, the larger the size of each cohort A_ℓ­
A_ℓ-1, the shallower the cover tree.  That is, ℓ_max is
expected to be of the order log_𝔻(n).

For each adult a_i ∈ A_ℓ, we must:


  * Sum and sort the weights of its children, by label.  The expected number of
children is n/𝔻^D.  Sorting the labels is O(L log(L)) for L labels.


  * Compute the distance (an O(D) operation) to each of its roughly
θ^ℓ n children. These distances yield the orphans and teens.  The
expected number of orphans is (1-θ^D)n/𝔻^ℓ, and the expected
number of teens is ( 1 - (1/2)^D)n/𝔻^ℓ.


  * For each orphan, determine whether it should be
adopted or liberated.  Compute the distance (an O(D) operation) to the T_1
friends or successors
of the T_1 friends of a_i.  This set has expected size 𝔻^2,
which is constant in n.


  * For each teen, determine whether it should be
exchanged.  Compute the distance (an O(D) operation) to the T_2
friends or successors of the T_2 friends of a_i.  This set has expected size 𝔻^2,
which is constant in n.


  * For each successor—a list of size 𝔻—update its friends by
computing the distance (an O(D) operation) to the successors of the T_3
friends of a_i.  This is expected to be 𝔻^2 distance computations.


For each adult at level ℓ, this is roughly n/𝔻^ℓ
operations, albeit with an additional penalty from D.  Level ℓ has at most 𝔻^ℓ adults, in
which case the total number of computations is roughly

    ∑_ℓ = 0^ℓ_max𝔻^ℓn/𝔻^ℓ
    O(𝔻^2)
    ∼ O(n log_𝔻n)O(𝔻^2).


It is interesting that the rapid growth from the ambient dimension D is
counteracted by the logarithm term.   In our experience with applications, this
is effectively linear in n, as the sum terminates earlier than this estimate.
This is because the estimate of ℓ_max by log_𝔻(n) is
typically irrelevant; the maximum level ℓ_max of the cover tree
is actually controlled by the ratio of the maximum to the minimum distances in
the pointcloud.  In applications, these are typically within several orders of
magnitude, so with θ = 1/2, we expect ℓ_max∼ 14.
Moreover, in our supervised-learning application in Section <ref>, we
abort the cover tree much earlier, once the low-entropy regions have been
found.



§ ENTROPY REDUCTION

For each label λ, let 𝒳|_λ denote the
(weighted) points in 𝒳 that have label λ.  For
any convex, compact Ω, let Ω_λ = Ω∩𝒳|_λ.  The total weight in Ω of each label is
w_λ = ∑{ w(x)  :  x ∈Ω_λ}, which by
Section <ref> corresponds to the relative likelihood of a point with
label λ being sampled from Ω.
Let W = ∑_λ w_λ.

To quantify the relative prominence of the various labels, we use
the information-theoretic notion of entropy.
Entropy is characterized by several nice properties, detailed in <cit.>.
The formula is

    S(Ω)    =
    - ∑_λw_λ/Wlog_L(
    w_λ/W)
       =
    -1/W∑_λ w_λlog_L(
    w_λ) + log_L(W).


If all weights are roughly equal, then entropy is 1.  If one label has high
weight and the others have weight 0, then the entropy is 0.  Generally, entropy
decreases as a smaller subset of labels becomes more likely.  Thus, low entropy
is a good indicator that the region Ω is particularly prominent for a
small subset of labels.
The question remains: how can we locate low-entropy regions?

Our search for low-entropy regions is motivated by the following elementary
theorems; however, we do not insist that these strict regularity hypotheses
hold for any given application.


Fix a function f:ℝ^D →ℝ that is continuous on
a compact, convex set Ω.
Fix a radius r>0.  Let μ_r(x) denote the centroid of f on the compact,
convex set Ω∩ B_r(x) (the closed ball).
There exists x ∈Ω such that μ_r(x) = x.

Proof: Note that μ_r(x) varies continuously with x, and apply Brouwer's fixed-point
theorem (<cit.>, II.22).


Suppose that f:ℝ^D→ℝ is an analytic function with local
max at x.   For any ϵ >0, there exists a radius r such that the
centroid μ_r(x) of f on B_r(x) has d(μ_r(x) ,x ) < ϵ.

Proof: Expand f as a power series around x, and approximate f by
truncating at the second-order term.  The level sets are concentric ellipsoids about x, as given
by the positive-definite symmetric matrix f”.


Suppose f and f̃ are positive, continuous functions defined on a
compact, convex region Ω⊂ℝ^D.
Suppose that f achieves a local max f(x) = M on Ω and that
f̃(x) ≤ε on Ω, with ε < M.  Then,
there is a radius r such that the average value of f on Ω∩
B_s(x) is greater than the average value of f̃ on Ω∩
B_s(x) for all s ≤ r.

Proof:   The mean value theorem for integrals.

Putting these together, we try to find a region Ω and a ball of radius
r_ℓ in that region whose center is near the mean of a particular label.
We ask whether the entropy of that ball is non-decreasing with radius;  that
is, does the entropy become lower (or remain the same) when the ball shrinks?
If true, this is consistent with the hypothesis that a local maximum of a
particular label occurs near the center of that ball, while other labels remain
roughly constant.  If false, then we subdivide the region and search around the
boundary of the original ball for other smaller regions of low entropy.

Specifically, the compact, convex regions corresponds to the children of an
adult a_i ∈ A_ℓ in the cover tree, and the sub-balls correspond to the
children of successors of a_i at level ℓ+1 in the cover tree.   This is
the heart of the CDER algorithm in Section <ref>.




§ COVER-TREE DIFFERENCING VIA ENTROPY REDUCTION

The input data is a ratio θ and a training set of pointclouds
X_1, X_2, …, X_N with labels λ_1, λ_2, …, λ_N,
respectively.    Set weights as in Section <ref>.




 §.§ Region Selection

Here we describe the heart of the CDER algorithm, wherein we search the adults
of the cover tree as a data-driven method to locate regions on which we want to
build distributional coordinates.  The construction and weighting of those
distributional coordinates occurs in Section <ref>.

For each level of the cover tree, define a subset Â_ℓ⊂ A_ℓ
of adults that are still potential candidates for distributional coordinates.
Set Â_0 = A_0 = { a_0}.  Once Â_ℓ = {}, we
.

The cover tree of 𝒳 is constructed level-by-level as in
Section <ref>.  At each level ℓ of the cover tree (until ), we
perform the following computation:

For each adult a_i ∈ A_ℓ, let β_ℓ(i) denote the set of children
of a_i at level ℓ. Let α_ℓ(i) = β_ℓ+1(i) denote the
set of children of a_i at level ℓ+1.  Let γ_ℓ(i) denote the
union of the children (at level ℓ-1) of the elders of a_i.
(If a_i belongs to cohort ℓ-1 or less, then γ_ℓ(i) =
β_ℓ-1(i).)
Note that α_ℓ(i) ⊂β_ℓ(i) ⊂γ_ℓ(i); see
Figure <ref>.  Then, we decide what to do with a_i:





  * if Â_ℓ = ∅, then .  Nothing more is needed.

  * for each adult a_i ∈Â_ℓ:


  * if a_i has only one child (itself) at level ℓ+1, then .  This
    region is useless.

  * elif a_i has only one child (itself) at level ℓ, then .  This
    region is useless.

  * elif S(α) ≤ S(β) ≤ S(γ) < 1, then
     a distributional coordinate using the children of a_i of
    each of the dominant labels.[A label is dominant if its points represent more
    than 1/L of the total weight of the children.  For a two-label task,
    this is simply the majority label.]

  * elif S(γ) ≤ S(α) ≤ S(β) < 1, then  a_i to
    Â_ℓ+1.  The smaller ball B(a_i, r_ℓ+1) may contain a
    lower-entropy region.

  * elif S(β) ≤ S(α) ≤ S(γ) < 1, then  the successors of
    a_i—other than itself!—to Â_ℓ+1.  The annulus B(a_i,
    r_ℓ) ∖ B(a_i, r_ℓ+1) may contain a lower-entropy region.

  * elif S(β) ≤ S(γ) ≤ S(α) < 1, then  the successors of
    a_i—other than itself!—to Â_ℓ+1. The annulus B(a_i,
    r_ℓ) ∖ B(a_i, r_ℓ+1) may contain a lower-entropy region.

  * elif S(γ) ≤ S(α) ≤ S(β) < 1, then .  The region B(a_i,
    r_ℓ) is useless.

  * elif if S(γ) ≤ S(β) ≤ S(α) < 1, then . The region B(a_i,
    r_ℓ) is useless.

  * else (meaning max{ S(α), S(β), S(γ)} = 1)  the
    successors of a_i to Â_ℓ+1.  There is too little
    information to decide whether B(a_i, r_ℓ) contains a low-entropy
    region, so re-analyze everything at the next level.



Recall that our method in Section <ref> produces the list
A_ℓ which is ordered by cohort, and within a cohort is ordered by
predecessor, and within a predecessor block is ordered by proximity to the
predecessor's per-label means.  Therefore, if a_i < a_j in this ordering,
then a_i tends to be nearer the mean of the most prominent label of a larger
region than a_j.  Hence, the distributional coordinates are sorted by
“granularity,” from coarse to fine.

In some scenarios, detailed analysis of clusters is more important than
efficient classification.  Therefore, one can implement a “non-parsimonious”
version by replacing all  operations with .  For
machine learning applications, the “parsimonious” version is far faster, as
the cover tree can be stopped much earlier than the limit suggested in
Section <ref>.



 §.§ Build Distributional Coordinates

For each a_i ∈ A_ℓ selected to , we construct a
distributional coordinate g_i in the following way.

For each of the dominant labels λ among the children of a_i, use
PCA/SVD to build a Gaussian model g_i,λ of mass 1.   Let Δ S
denote the entropy difference on β_ℓ caused by erasing the
non-dominant labels.
We amplify (or attenuate) each
Gaussian model by the
coefficient

    m_i = w_λ  (1 - Δ S) (r_ℓ)^D.

The distributional coordinate is m_i g_i,λ(x).

One can experiment with different weighting systems; we have found these to be
generally effective, and we justify this choice for the following geometric
reasons:


  * The weight represents the relative likelihood that a point of the given
label was selected from the region on which this coordinate was detected.
Perhaps, by chance in our selection process, a single region was detected at a
later level in the form of several smaller regions.  The weight of the smaller
regions will sum to the weight of the larger region.


  * The entropy term is the “certainty” with which this region was
selected.  The entropy Δ S is the information that is lost if the
other labels are removed from the region, so 1-Δ S is the information
that remains. In other words, this term penalizes regions according to the
impurity of their labels.


  * The (r_ℓ)^D = (r_0 θ^ℓ)^D term accounts for the relative
remoteness of the sample region.  All else being equal, remote regions are
particularly distinctive. Because of the cover-tree construction, remoteness
can be measured by the size of the cover-tree ball.
For example,  suppose that two different
regions are detected by CDER, each with the same entropy, weight, and
configuration of points; the Gaussian coordinates are translates of
one-another.  Suppose that the first region is very distant from the rest of
the pointcloud, so it is detected at an early level ℓ.  Suppose that the second region
is surrounded closely by other data, so it is not detected by the cover tree
until a later level ℓ' > ℓ. The first region's relative volume is
θ^(ℓ-ℓ')D.





§ EXAMPLES

This section illustrates the promise of CDER on several synthetic examples
of varying levels of complexity that are designed to stress certain features of
the algorithm.



 §.§ Blobs

This example shows that CDER ignores high-density regions unless they
have low entropy.

Consider two labels[Labels are indicated by colors, and we have intentionally chosen colors that seem to show up non-catastrophically in black-and-white printouts. Nonetheless, we recommend printing this paper in color or at least reading this section on the computer screen!]: 0/magenta and 1/green.
The cloud collection has 25 pointclouds of each label.
Each magenta pointcloud consists of:


  * 100 points sampled from a standard normal
distribution

  * 2 points sampled from a normal distribution
with σ = 0.2 and mean (4,0)

  * 2 points sampled from a normal distribution
with σ = 0.2 and mean (5,0)

  * 2 points sampled from a normal distribution
with σ = 0.2 and mean (-3,0)

  * 2 points sampled from a normal distribution
with σ = 0.2 and mean (-6,0)

Each green pointcloud consists of


  * 100 points sampled from a standard normal
distribution

  * 2 points sampled from a normal distribution
with σ = 0.2 and mean (-4,0)

  * 2 points sampled from a normal distribution
with σ = 0.2 and mean (-5,0)

  * 2 points sampled from a normal distribution
with σ = 0.2 and mean (3,0)

  * 2 points sampled from a normal distribution
with σ = 0.2 and mean (6,0)

This produces the cloud collection in Figure <ref>.
The output of CDER appears in Figure <ref>


In a typical run, the algorithm terminates at cover-tree level 8, even though
the complete cover tree of 𝒳 ends at level 14.  A total
of 26 distributional coordinates are produced, and their weights vary by about
3 orders of magnitude.
Table <ref> shows how conservative the entropy-reduction process is at
selecting regions to generate distributional coordinates.



To judge accuracy on testing data in a supervised-learning scenario, we need a
method to label a test pointcloud using the distributional coordinates.  Each
test pointcloud X_test is mapped to a point in ℝ^26.
Many sophisticated methods are possible, but for simplicity, we simply ask: For a
given pointcloud X_test, which is bigger:  the Euclidean norm
of the magenta Gaussian coordinates evaluated on X_test,

    √(∑_g magenta(∫_X_test g)^2)

or the Euclidean norm of the green Gaussian coordinates evaluated on
X_test,

    √(∑_g green(∫_X_test g)^2)?


With this simple comparison, the algorithm achieves 100% accuracy in a 5-fold
cross-validation of this cloud collection with a 80/20 training/testing split.
More precisely, for each run of the cross-validation, we take 80 percent of the point clouds,
turn these into a cloud collection, build distributional coordinates using this cloud collection, and then
test the results on the remaining 20 percent of the point clouds. This entire procedure is then repeated five times.



Moreover, the relative masses of the 28 distributional coordinates vary over
four orders of magnitude, so for this sort of comparison, one could dispose of
many of them in the interest of speed while preserving accuracy.
Note that these distributions are contrived to have the same mass (0-moment)
and mean (1-moment) and variance.  Elementary statistical tests would not
distinguish them; 2-moments or skewness tests would be necessary.



 §.§ Blocks

This example shows that CDER is fairly robust against background noise that
prevents strict separation.  It also demonstrates that smoothness of the
underlying distributions is not necessary for good results.

Consider two labels: 0/magenta and 1/green.
The cloud collection consists of 100 magenta pointclouds and 100 green
pointclouds.
Each magenta pointcloud in generated by sampling 30 points uniformly in the
unit square, as well as 2 extra points sampled uniformly in a 0.1× 0.1 square at
(1/4, 1/4), as well as
2 extra points sampled uniformly in a 0.1× 0.1 square at
(1/2, 1/2).
Each green pointcloud in generated by sampling 30 points uniformly in the
unit square, as well as 2 extra points sampled uniformly in a 0.1× 0.1 square at
(3/4, 3/4), as well as
2 extra points sampled uniformly in a 0.1× 0.1 square at
(1/2, 1/2).
See Figures <ref> and <ref>.

Using the same simple comparison as in Section <ref>, the algorithm
achieves 88% accuracy despite the high background noise.



 §.§ Deep Field

In this example, we demonstrate that Gaussian mixtures are (nearly)
“fixed points” of CDER.  It also demonstrates that the algorithm can handle
unequal sample sizes via the weighting system in Section <ref>.




Consider two labels: 0/magenta and 1/green.
The underlying distribution is build from 50 Gaussian distributions.
For each of these, we choose a random label.  We choose a mean point, uniformly
on the 10 × 10 square.  We choose a pair of orthonormal covariances uniformly
from 0 to 0.5.  We choose a random rotation angle.  We choose an
amplification factor (that is, a sample size) uniformly between 50 and 5000.
For a particular random seed, we arrive at some combined density function f_0 and
f_1, as in Figure <ref>

Then, for each label, we produce between 20 and 40 pointclouds,
each with between 50 and 500 points.  For a particular random seed, we arrive
at the cloud collection in Figure <ref>.

The CDER algorithm builds a Gaussian mixture model focused on regions of low
entropy (Figure <ref>), so it should not be surprising that
it builds Gaussians near the original Gaussians, as in
Figure <ref>.






 §.§ Three Labels

For simplicity, the previous examples have involved only two labels.    In a
two-label system, a low-entropy region has exactly one dominant label.
However, the algorithm is sensible for any number of labels, and it may be that
the low-entropy regions are dominated by multiple labels.  Hence, an ensemble
of regions may be necessary to distinguish pointclouds.


Consider three labels:  0/magenta, 1/green, and 2/orange.
Let f_0 be a triple of standard normal distributions: one each at (0,0),
(4,0), and (-2, 2√(3)).
Let f_1 be a triple of standard normal distributions: one each at
(0,0),(-2,2√(3)), and (-2, -2√(3)).
(-2,2√(3)), and (-2, -2√(3)).
Let f_2 be a triple of standard normal distributions: one each at (0,0),
(-2,-2√(3)), and (4, 0).  See Figure <ref>.

The CDER algorithm detects these shared regions perfectly, achieving 100% in a
5-fold cross-validation of this cloud collection with a 80/20 training/testing
split, using the same simple comparison method as in Section <ref>.




§ STABILITY

Ideally, one would like to prove that CDER is stable or backward-stable,
in the sense of <cit.>.
To formulate stability or backward-stability, CDER must be expressed as a
numerical approximation of a formal function between normed vector spaces.

Let C denote the vector space of piecewise continuous
functions on
our compact domain.  Define a norm on C × C by
(f_1,f_2) = √(f_1_2^2 + f_2_2^2), where each term is the L_2 norm
on our (compact) domain.
Given two “input ”functions (f_1, f_2) ∈ V we consider them as
density functions for each of two labels.  Consider the
pointwise entropy,

    H(x) =
    -f_1(x)/f_1(x)+f_2(x)log_2(f_1(x)/f_1(x)+f_2(x))
    -f_2(x)/f_1(x)+f_2(x)log_2(f_2(x)/f_1(x)+f_2(x)).


Let T denote the formal function T: C
× C → C × C by
T(f_1, f_2) = (g_1, g_2)
where

    g_1(x)= 
    H(x)f_1(x),   if f_1(x) > f_2(x),
    
    0,   otherwise.

and

    g_2(x)= 
    H(x)f_2(x),   if f_2(x) > f_1(x),
    
    0,   otherwise.

A numerical algorithm T̃ for T can be obtained roughly as “sample
equally from f_1 and f_2, and apply CDER to generate weighted
distributional coordinates.  The weighted sum of the distributional coordinates
for label 1 is g_1, and the sum of the distributional coordinates for label 2
is g_2.”  A stability result such as this is a subject of future work. We
have not yet succeeded in proving either stability statement for CDER, but the
high cross-validation of the examples above is promising.




§ DISCUSSION
 

This paper introduced CDER, a data-driven, label-driven feature extraction method for collections of labeled pointclouds. CDER is fast, both in terms of theoretical complexity
and in initial tests on examples. It does not require the user to choose any tuning parameters, and the geometric meaning of its output features are transparently clear.
This section outlines some future directions and makes a few generalizing remarks.

The attentive reader of Section <ref> may notice that cover trees can be defined in an arbitrary metric space, and indeed this is how they
were originally defined <cit.>.
We use them to construct distributional coordinates, and our algorithm for doing so (Section <ref>) demands that we be able to
quickly compute means of sets of points. While the Fréchet Mean (<cit.>, IX) can be defined in any metric space, there are not always fast algorithms for its computation (it also need not exist in all cases).
And so while more generality could be achieved in the definition
of CDER, we only make our complexity statements for cloud collections in a common Euclidean space.

All examples in Section <ref> were artificial, and simply intended to help visualize the novel and fairly-technical CDER
algorithm and emphasize its key properties.
Future work will involve applications of CDER to real data from a variety of real fields, including large administrative datasets of interest to social scientists
and to vehicle tracking <cit.>.

We also hope that CDER will prove useful as a feature-extraction method in topological data analysis <cit.>,
since persistence diagrams can be thought of as point clouds in the plane. Future work will compare the performance of CDER against other such feature-extraction methods (for example, <cit.> and <cit.>).

Finally, we recall the weighting discussion in Section <ref>, where we used the simplifying assumptions that each color/label was equally likely and that each point
within a single point cloud should be given equal weight. We note that CDER can be easily adapted to accommodate other prior assumptions about relative likelihoods of labels
or even prior assumptions about outlier status of certain points in a cloud, say as part of a Bayesian learning process.
















plain

