Convex clustering via `1 fusion penalization

arXiv:1412.0753v5 [stat.ME] 28 Dec 2016

P ETER R ADCHENKO AND G OURAB M UKHERJEE ∗

Abstract
We study the large sample behavior of a convex clustering framework, which minimizes the sample within cluster sum of squares under an `1 fusion constraint on the
cluster centroids. This recently proposed approach has been gaining in popularity,
however, its asymptotic properties have remained mostly unknown. Our analysis is
based on a novel representation of the sample clustering procedure as a sequence of
cluster splits determined by a sequence of maximization problems. We use this representation to provide a simple and intuitive formulation for the population clustering
procedure. We then demonstrate that the sample procedure consistently estimates its
population analog, and derive the corresponding rates of convergence. The proof conducts a careful simultaneous analysis of a collection of M-estimation problems, whose
cardinality grows together with the sample size. Based on the new perspectives gained
from the asymptotic investigation, we propose a key post-processing modification of
the original clustering framework. We show, both theoretically and empirically, that
the resulting approach can be successfully used to estimate the number of clusters in
the population. Using simulated data, we compare the proposed method with existing
number of clusters and modality assessment approaches, and obtain encouraging results. We also demonstrate the applicability of our clustering method for the detection
of cellular subpopulations in a single-cell virology study.
Some key words: Convex Clustering; Fusion Penalties; Number of Clusters; Rates of Convergence

1

Introduction

Clustering is one of the most popular statistical techniques for unsupervised classification
and taxonomy detection (Hartigan, 1975; Kaufman and Rousseeuw, 2009). One serious
∗

University of Southern California

1

limitation of the traditional methods, such as k-means, is the non-convexity of the corresponding optimization problems. Recently, several convex clustering algorithms have
been proposed (Xu et al., 2004; Bach and Harchaoui, 2008; Chi and Lange, 2013). Speed
and scalability of these algorithms make them increasingly popular for cluster analysis of
massive modern datasets. These approaches use convex relaxations of the traditional nonconvex clustering criteria, however, they do not naturally inherit the statistical properties
associated with the original methods. Here we study the large sample behavior of a popular
convex clustering framework that is based on an `1 fusion penalty (Hocking et al., 2011).
Consider the problem of clustering n observations, x1 , . . . , xn , which are sampled from
a Euclidean space, Rd . The well-studied k-means approach (MacQueen et al., 1967; Hartigan, 1978; Pollard, 1981, 1982; Jain, 2010) is based on minimizing the within cluster sum
P
of squares, ni=1 kxi − αi k22 , with respect to the cluster centroids, α1 , . . . , αn , under the
restriction that the number of distinct cluster centroids is at most k. This restriction can be
viewed as an `0 constraint on the centroids. Motivated by the Lasso and its variants (Tibshirani, 1996; Tibshirani et al., 2005), which successfully use the `1 constraint as a surrogate
for the NP-hard `0 constraint, Hocking et al. (2011) consider the following modification of
the k-means clustering criterion:
min

α1 ,...,αn

n
X

kxi − αi k22 subject to

i=1

X

αi − αj

1

≤ t.

(1)

1≤i<j≤n

When t = 0, the `1 penalty fuses all the cluster centroids together. Thus, all the observaP
tions are placed in the same cluster. When t ≥ i<j kxi − xj k1 , we have αi = xi for
all i, and, thus, each observation forms its own cluster. Varying t between the two extremes
creates a path of solutions to the regularized clustering problem. Note that the Lagrangian
P
P
form of the above criterion, minαi ni=1 kxi − αi k22 + λ i<j αi − αj 1 , is separable
across dimensions. Consequently, the corresponding optimization problem reduces to independently minimizing d univariate convex clustering criteria.
Thus, to understand the large sample behaviour of the multivariate solution, it is sufficient to focus on the analysis of the univariate clustering criterion,
n
X
X
min
(xi − αi )2 + λ

α1 ,...,αn

i=1

1≤i<j≤n

2

αi − αj .

(2)

As the penalty parameter λ varies from 0 to ∞, each corresponding solution determines
a cluster partition. We are interested in the asymptotics of the entire collection of such
partitions, which we view as the outcome of the sample clustering procedure.
Summary of the Main Contributions. We analyze the large sample behavior of the
sample clustering procedure determined by the solution path for criterion (2). We develop
a simple and intuitive formulation for the population clustering procedure, show that under some very mild regularity conditions the sample procedure consistently estimates its
population analog, and derive the corresponding rates of convergence.
More specifically, we first demonstrate that the path of solutions to (2) determines a
clustering tree, which can be formed by either successive merges of clusters, in a bottom
up fashion, or successive splits, in a top down approach. We then study the asymptotic behavior of the full clustering tree by representing each split as a solution to a maximization
problem. We define the corresponding population clustering procedure in a similar fashion, but replace sample averages with the corresponding expected values. The asymptotic
analysis is significantly complicated by the fact that, unlike in the standard M-estimation
setup (e.g. van der Vaart and Wellner 1996; van der Vaart 1998), the number of maximization problems at the sample level tends to infinity together with n, and the number of
the corresponding population problems is infinite. We establish consistency and the rates of
convergence of the sample clustering procedure through a careful analysis of the population
procedure and the corresponding empirical process.
Motivated by the results of our large sample investigation, we introduce a key postprocessing modification to the sample clustering procedure. We show, both theoretically and
empirically, that the resulting approach can be successfully used to estimate the number
of clusters in the population. We also compare the new methodology with a wide variety
of existing modality assessment and number of clusters approaches. Our results provide
strong support for the use of fusion penalization in clustering.
Connections to Related Work. Hocking et al. (2011), Chi and Lange (2013) and Tan
and Witten (2015) have studied modifications of optimization problem (1). These include
using `2 or `∞ regularization, as well as incorporating weights (Pelckmans et al., 2005;
Lindsten et al., 2011; Zhu et al., 2014). The large sample analysis in the papers listed
above focusses on showing that if the distance between clusters grows at a sufficiently fast
rate, then the corresponding method can separate the groups perfectly. Here we consider

3

a completely different perspective and investigate the asymptotics of a clustering approach
in the classical sense of Pollard (1981). We study a clustering procedure that is applied
to a random sample, and analyze its convergence to the outcome of the corresponding
population procedure, which is based on the underlying probability distribution. As we
point out in Section 5, the general framework of our theoretical analysis has the potential
to handle the aforementioned modifications of the optimization problem.
The criterion in (1) can be viewed (Hocking et al., 2011) as a convex relaxation of the
hierarchical clustering criterion (Hartigan, 1975). However, as clustering is a very mature
subject, approaches built on several other philosophies are also widely used in practice. A
detailed review of clustering methods can be found in Kaufman and Rousseeuw (2009).
One of the most popular methods is the k-means algorithm (MacQueen et al., 1967), which
follows a partitioning approach for making clusters. Other popular partitioning methods,
such as PAM (Kaufman and Rousseeuw, 1990) and CLARA (Kaufman and Rousseeuw,
1986), are based on the k-medoids algorithm. Density driven approaches, which include
mixture model based methods, such as Fraley and Raftery (2002) and Li (2005), as well as
non-parametric methods (see Li et al. 2007 and the references therein), provide a flexible
clustering framework, while spectral clustering methods, such as (Belkin and Niyogi, 2001;
Rohe et al., 2011; Shi et al., 2009) perform efficient dimension reduction before segmenting
the data. In our empirical analysis we compare the performance of the proposed approach
with most of the aforementioned clustering methods.
The `1 penalty, which is extensively used for variable selection (Tibshirani, 2011), also
finds its use in trend filtering (Tibshirani, 2013) and high-dimensional clustering problems
(Soltanolkotabi and Candés, 2012; Witten and Tibshirani, 2010). Another related approach,
the fused Lasso (Rinaldo et al., 2009; Tibshirani and Walther, 2005; Hoefling, 2010), deals
with applications having ordered features and checks for local constancy of their associated
coefficients. This approach penalizes the successive differences of the coefficients. Shen
and Huang (2010); Shen et al. (2012); Ke et al. (2013); Bondell and Reich (2008) have
proposed methods based on fusion penalties, which apply to all the pairwise differences
of coefficients. These approaches can successfully recover the grouping structure of predictors in a high-dimensional regression setup. However, the theory developed for these
methods focusses on the homogeneity of regression coefficients and cannot be applied in
the unsupervised clustering setup considered in this paper.

4

Organization of the Paper. In Section 2 we derive two equivalent algorithmic representations of the sample clustering procedure, which we use to formulate the corresponding
population procedure. Section 3 contains our main results, in which we establish consistency and the rates of convergence. Our asymptotic analysis reveals that an overwhelming
majority of the sample clusters are in some sense negligible. Motivated by this observation, we introduce a key post-processing modification to the clustering procedure. In
Section 4 we conduct a detailed empirical analysis of our approach. More specifically, we
use simulated data to show its strong performance relative to popular existing approaches
for assessing modality and estimating the number of clusters. We also illustrate the use
of our method in analysis of single-cell virology datasets. All the proofs, together with
additional technical details, are relegated to the Supplementary Material.

2

Sample and Population Clustering Procedures

In this section we derive two equivalent representations of the sample clustering procedure. First, we develop a computationally efficient merging algorithm for producing a path
of solutions to the clustering criterion (2). Then, in order to understand the large sample
behavior of the solution path, we introduce an equivalent splitting procedure, which can
recover all the corresponding cluster splits by solving a sequence of maximization problems. We use the splitting representation to define the population clustering procedure, and
describe its basic properties.

2.1

Equivalent Representations for the Sample Solution Path

Note that a solution path for problem (2) could be produced using the highly general fused
lasso algorithm in Hoefling (2010). Instead, we obtain a very simple and computationally efficient fitting procedure by analyzing our clustering criterion, (2), directly. The path
algorithm we describe here is a bottom up procedure, which starts at λ = 0, with each
observation forming its own cluster, and then gradually merges suitable clusters as λ increases. Fix λ, and suppose that C is one of the clusters identified by the solution to the
optimization problem (2). Write αC for the centroid of cluster C, and denote the corresponding cluster average by X C . As pointed out in Hocking et al. (2011), the first order

5

conditions for criterion (2) imply
αC = X C + λ

X

sign(αj − αC ).

(3)

j,αj 6=αC

Until the cluster partition or the ordering or the centroids are modified, parameter λ is the
only component on the right-hand side of the equation that can change. Thus, equation (3)
provides a simple way of tracking the piecewise linear paths of the centroids αi . Another
consequence of the first order conditions is that as λ increases, the only way the clusters get
modified is some of them get merged together (Hocking et al., 2011). Hence, we can store
the full cluster partition path by keeping track of the merges and the corresponding values
of the tuning parameter λ. Algorithm 1 makes this idea precise, and Theorem 1 provides a
rigorous justification. Here we use | · | to denote the cardinality of a set.
INITIALIZE:
Sort data in ascending order and store them as xn = {x1 , . . . , xn }.
Set K, the number of clusters, equal to n. For each i in 1, ..., n, set Ci = {xi }.
REPEAT:
Find the adjacent centroid distances standardized by cluster sizes:
d(j, j + 1) ← (X Cj+1 − X Cj ) /(|Cj | + |Cj+1 |).
Find the clusters that minimize this distance: j ∗ ← arg minj d(j, j + 1).
Merge the clusters that were found: Cj ∗ ← Cj ∗ ∪ Cj ∗ +1 .
Store the above merge and the corresponding λ value: λ = d(j ∗ , j ∗ + 1).
Relabel the remaining clusters: for j > j ∗ set Cj ← Cj+1 .
Reduce the total number of clusters: K ← K − 1.
UNTIL K = 1.
OUTPUT: Sequence of cluster merges and corresponding λ values.
Algorithm 1: Merging Algorithm
The following result shows that Algorithm 1 reproduces the sequence of cluster partitions and the corresponding λ values from the optimization problem (2). In the proof,
which is provided in the Supplementary Material, we also verify that the sequence of λ
values, corresponding to successive merges in Algorithm 1, is increasing.

6

Proposition 1. Suppose that the observations are generated from a continuous distribution.
Then, with probability one, the sequence of merges and λ values produced by the merging
algorithm is the same as the sequence corresponding to the optimization criterion (2).
For the asymptotic analysis, it is helpful to recover the sequence of cluster partitions
in a top down approach: we start with everything in one cluster and then split the clusters
iteratively. We call a representation of the cluster C as C = C1 ∪ C2 a split if max C1 <
min C2 . The full collection of splits corresponding to the optimization problem (2) is given
by the splitting procedure, described in Algorithm 2 below. Proposition 2, proved in the
Supplementary Material, provides theoretical justification. In particular, it shows that each
of the cluster splits is chosen to maximize the distance between the two sub-cluster means.
INITIALIZE:
Sort data in ascending order and store them as xn = {x1 , . . . , xn }.
Set the current partition of xn to xn .
REPEAT:
Select one cluster, C, with |C| > 1, from a current cluster partition of xn .
Find a split partition C = C1 ∪ C2 , that maximizes the distance X C2 − X C1 .
Store the split C = C1 ∪ C2 and the corresponding value λ = (X C2 − X C1 ) /|C|.
Replace C with C1 ∪ C2 in the current partition of xn .
UNTIL: All the clusters in the current partition of xn are of size one.
OUTPUT: Collection of cluster splits and corresponding λ values.
Algorithm 2: Splitting Procedure

Proposition 2. Suppose that the observations are generated from a continuous distribution.
Then, with probability one, the collection of splits and corresponding λ values produced by
the splitting procedure in Algorithm 2 exactly matches the sequence of merges and corresponding λ values produced by the merging algorithm .
Note that, unlike the merging algorithm, the splitting procedure does not provide a
computationally efficient way for producing the clustering tree. Instead, we use the splitting
procedure to understand the large sample behavior of the sample clustering procedure. It
is reasonable to expect that, as n tends to infinity, the collection of splits in the sample
procedure should resemble the collection of splits in an analogous procedure defined on
7

the population. The population procedure can be defined by replacing the averages with
the corresponding conditional means. The formal definition is given in the next section.

2.2

Population Clustering Procedure

For the remainder of the paper we assume that the underlying distribution has a finite first
moment and a real valued density, f . For concreteness, we focus on the case where the
support of the distribution is of the form (L0 , R0 ), where −∞ ≤ L0 < R0 ≤ ∞. Thus,
every open interval in (L0 , R0 ) contains positive probability. Given an interval (l, r) ⊆
(L0 , R0 ), we write µl,r for the population conditional mean on (l, r),


−1

Zr

µl,r = 

Zr

f (x)dx
l

xf (x)dx.

(4)

l

We set µr,r = r, by continuity. Given an interval (L, R) ⊆ (L0 , R0 ), we define
GL,R (a) = µa,R − µL,a ,

(5)

for a ∈ [L, R]. Note that GL,R (L) = µL,R − L and GL,R (R) = R − µL,R .
According to the results in Section 2.1, the sample clustering procedure determines the
split partition of a cluster by maximizing the distance between the empirical sub-cluster
means. We define the population clustering procedure by analogy. Given a cluster (L, R),
the population procedure chooses the split that maximizes the distance between the population sub-cluster means. In other words, it finds a point s that maximizes GL,R , then
partitions (L, R) into subintervals (L, s) and (s, R), on which the procedure is repeated.
If s is an interior point of (L, R), we call it a split point, and we call the corresponding
partition a split. Otherwise, the population procedure essentially wants to split off an endpoint, which forces the cluster to be truncated rather than split. More formally, given a
cluster (L, R), we distinguish three types of truncation, as specified below.
Definition 1.

(i) if arg max Gl,R = {l} for all l ∈ [L, L∗ ), and arg max GL∗ ,R 6= {L∗ },

then the interval (L, R) is truncated from the left to (L∗ , R∗ ), where R∗ = R;
(ii) if arg max GL,r = {r} for all r ∈ (R∗ , R], and arg max GL,R∗ 6= {R∗ }, then (L, R)
is truncated from the right to (L∗ , R∗ ), where L∗ = L;
8

(iii) if there exists a continuous decreasing function l 7→ Rl , satisfying RL = R, for which
arg max Gl,Rl = {l, Rl } for all l ∈ [L, L∗ ), and arg max GL∗ ,RL∗ 6= {L∗ , RL∗ }, then
(L, R) is truncated, in a two-sided fashion, to (L∗ , R∗ ), where R∗ = RL∗ .
Note that we incorporated a continuity requirement into the definition of a two-sided
truncation. In the next subsection we give regularity conditions under which this requirement is satisfied. We are now ready to formulate the full population clustering procedure.
INITIALIZE:
Set the current cluster collection, Σ, equal to {(L0 , R0 )}.
REPEAT:
Select one non-empty cluster, (L, R), from the current cluster collection, Σ.
− If the maximum value of GL,R is achieved at a point s in (L, R), then store s
as a split point and replace (L, R) in Σ with (L, s) and (s, R).
− Otherwise, replace (L, R) in Σ with the interval (L∗ , R∗ ) from Definition 1.
UNTIL: The current cluster collection, Σ, consists only of empty clusters.
OUTPUT: Set of split points.
Algorithm 3: Population Clustering Procedure
The collection of population split points determines the corresponding clusters. For
example, consider the symmetric mixture of two Gaussian distributions examined in Figure 1. The population procedure identifies one split point, located at zero. This specifies
the population cluster partition: (−∞, 0) ∪ (0, ∞).
Given an underlying distribution, Algorithm 3 determines the exact behaviour of the
population clustering procedure. In Section 4 of the Supplementary Material we document
the performance of the population procedure for a variety of Gaussian mixtures. In Section 2.3 we establish an important fact that the population procedure produces no splits for
unimodal distributions. We also provide conditions under which the population procedure
is well defined, by which we mean that it implements finitely many uniquely defined steps.

2.3

Properties of the Population Procedure

The proofs of the results established in this section are provided in the Supplementary
Material. We first consider an important special case, where the underlying distribution
9

5.0
4.5
4.0
3.5

GL, R(a)

−5

−3

−1

1

3

5

−5

−3

A

−1

1
B

3

5

−5

−3

−1

1

3

5

C

Figure 1: Population criterion function, GL,R , corresponding to the Gaussian mixture
0.5 N (−2, 1) + 0.5 N (2, 1), is plotted for three choices of (L, R), such that R = −L and
GL,R (L) = GL,R (R) = max GL,R . The population clustering procedure continuously truncates the support of the distribution, symmetrically in a two-sided fashion, until max GL,R
can be achieved at an interior point (plot C). Then, the procedure places a split point at
zero. Note that the resulting sub-clusters are then truncated down to empty sets according
to Proposition 3 in Section 2.3.
is unimodal. We suppose that the density f is either strictly monotone on its support,
(L0 , R0 ), or there exists a point c for which f is strictly increasing on (L0 , c) and strictly
decreasing on (c, R0 ). The following result shows that in this setting, under just a continuity
assumption on f , the population procedure is unique and does not reveal any clusters.
Proposition 3. If f is continuous and unimodal, then the population clustering procedure
is uniquely defined and produces no splits.
We now move to the general setting. The following simple regularity condition ensures
existence of a population clustering procedure with finitely many steps, as we demonstrate
in the proof of Proposition 4 below.
C1. Density f is nonzero and differentiable on (L0 , R0 ). It has finitely many modes and,
at each of its interior modes, admits a non-constant Taylor approximation.
Remark.

The last requirement means the following: for each interior mode c, there

exists a positive integer k, such that f is k times differentiable at c with f (k) (c) 6= 0.

Note that the differentiability assumption can be slightly relaxed: for example, the results that follow hold for continuous piece-wise linear densities with no constant segments.
However, we prefer to keep this assumption, as it simplifies the presentation of the results.
10

To address the question of uniqueness, consider the following counterexample. Suppose the underlying distribution is uniform on (L, R). Then, the criterion function GL,R is
constant on its domain, and the population procedure applied to cluster (L, R) may place
a split point anywhere in its interior. Thus, there are infinitely many versions of the population clustering procedure. The following regularity condition explicitly rules out such
settings, by requiring that the interior of (L, R) contains at most one maximizer of GL,R .
C2. When the population procedure performs a split, the location of the split point is
uniquely determined.
Note that condition C2 holds for each distribution with a continuous unimodal density, as
a direct consequence of Proposition 3. In the proof of Proposition 3 we also show that C2
holds for all bimodal densities, provided the smoothness condition, C1, is satisfied. The following result establishes existence and uniqueness of the population clustering procedure
in the general setting.
Proposition 4. If regularity conditions C1 and C2 are satisfied, then the population clustering procedure is uniquely defined and implements finitely many steps.
In Section 3 we show that under the same regularity conditions, C1 and C2, the sample
procedure consistently estimates its population counterpart.

3

Main Results

In this section we show that the clustering tree produced by criterion (2) consistently estimates the clustering tree produced by the population procedure defined in Section 2.2. We
also derive the corresponding rates of convergence and propose a novel post-processing
modification of the sample clustering procedure. Recall that we assume a finite first moment for the underlying distribution.

3.1

Consistency

We start with some useful notation. In both the sample and the population, each split is
characterized by a triple (L, s, R), where interval (L, R) is the cluster being split, and s
is the split point, located inside (L, R). We write PL,R for the probability assigned to the
11

interval (L, R) by the underlying distribution. For a split sp = (L, s, R) we define its
size as size(sp) = min {PL,s , Ps,R }. When the probabilities in the above definition are
c
replaced with the corresponding sample frequencies, we write size(sp)
for the resulting
quantity, and refer to it as the empirical size.
The set of all the population splits, denoted by S, defines the population clustering tree.
b defines the sample tree. The cardinality of Sb
Similarly, the set of all sample splits, S,
tends to infinity as the sample size grows. Alternatively, according to Proposition 4, under
mild regularity conditions the population procedure produces finitely many splits, together
with some truncations. To establish consistency, we divide the sample splits into “big”
and “small”, based on their empirical size, then show that the first group converges to the
population splits, while the second is asymptotically negligible. The formal definition is
given below. We write dH for the Hausdorff distance between subsets of a Euclidean space.
c
Definition 2. Write Sbα for the set {sp ∈ Sb : size(sp)
> α} and let α∗ be the smallest
split size in the population procedure. We call the sample clustering procedure strongly
consistent if, for each α in (0, α∗ ), the following statements hold almost surely,
|Sbα | → |S|

(6)

dH (Sbα , S) → 0

and

max{size(sp) : sp ∈ Sb \ Sbα } → 0.

(7)
(8)

If we replace almost sure convergence with convergence in probability, we have a weaker
notion of consistency, which holds automatically when the sample procedure is strongly
consistent. In particular, displays (6) and (7) imply that, except on a set of probability tending to zero, there is a one to one correspondence between Sbα and the set of all population
splits, such that each split in Sbα converges to its population counterpart with respect to the
usual Euclidean distance. The next result, which is proved in the Supplementary Material,
establishes consistency of the sample procedure.
Theorem 1. Suppose that regularity conditions C1 and C2, given in Section 2.3, are satisfied. Then, the sample clustering procedure is strongly consistent.
Remark.

It follows from the proof that the result continues to hold if we replace
c with size in the definition of Sbα and/or replace size with size
c in display (8).
size

12

If condition C2 is violated, and the locations of the population splits are not uniquely determined, a modification of Theorem 1 continues to hold. More specifically, suppose that
the number of versions of the population procedure is finite. Then, the sample procedure
converges to the set of population versions, rather than a specific one. In other words, the
number and the locations of the big sample splits approach the corresponding quantities for
an appropriately chosen population version, where the choice depends on the sample.
Consider the important case of a unimodal underlying distribution. Proposition 3 in
Section 2.3 and the proof of Theorem 1 imply that in this case condition C1 is not required
for consistency. Note also that the population procedure produces no splits, by Proposition 3. It follows that all of the sample splits are uniformly asymptotically negligible.
Corollary 1. If f is continuous and unimodal, then the maximum size of all the splits in
the sample clustering procedure goes to zero almost surely.
In the next section we extend the results in Theorem 1 by establishing the rates of
convergence for the sample clustering procedure.

3.2

Rates of Convergence

To establish the rates of convergence for the sample splitting procedure, we need an additional regularity condition. We use the term population cluster to refer to all intervals that
appear along the path of the population procedure.
C3. For each population cluster (L, R) and each t ∈ arg max[L,R] GL,R , if t ∈ (L, R),
then G00 (t) 6= 0, otherwise G0L,R (t) 6= 0.
The requirement on GL,R (t)00 , imposed for each population split (L, t, R), is the standard
M-estimation assumption that requires the second derivative of the population criterion
function to be nonsingular at the population maximum (e.g. van der Vaart and Wellner
1996; van der Vaart 1998). The requirement on G0L,R is the analog of the aforementioned
M-estimation assumption in the case where the population criterion function, GL,R , is maximized at an endpoint of [L, R], rather than in the interior. In this case, the behaviour of
GL,R near its maximum is characterized by the first derivative, rather than by the second.
b ) contain all the sample splits sp = (L, b
Let S(τ
t, R), for which the sample frequencies
of (L, R), (−∞, b
t) and (b
t, ∞) are greater than or equal to τ . In Theorem 2 we restrict our
13

b ), for arbitrarily small but positive τ . Without this reattention to the sample splits in S(τ
striction, the rate of convergence in (10) would change. In particular, split sizes larger than
Op (log n/n) are produced when the sample procedure is applied to intervals whose widths
tend to zero. Also, larger split sizes may appear near the boundary of the support of the distribution. The general approach used in the proof of Theorem 2 would also establish these
slower rates of convergence. However, the exact form of the new rates depends on the behavior of the density near the boundary of the support and on the aforementioned intervals
of negligible width. Instead, we present a clean result, with just one rate of convergence
for all the small sample splits, while only imposing some simple regularity conditions.
Theorem 2. Suppose that regularity conditions C1-C3 are satisfied. Let α∗ be the smallest
split size in the population procedure. Then, for each α in (0, α∗ ),
dH (Sbα , S) = Op n−1/3



and

b )} = Op (log n/n) .
max{size(sp) : sp ∈ (Sb \ Sbα ) ∩ S(τ
Remark.

(9)
(10)

As we point out in the proof, if the domain, (L0 , R0 ), of the underlying

distribution is bounded, then we can remove the lower bounds on the sample frequenb ) by assuming, instead, that
cies of (−∞, b
t) and (b
t, ∞) from the definition of S(τ
f (L0 ) and f (R0 ) are nonzero. It also follows from the proof that the result continues
c in display (10).
to hold if we replace size with size

The proof of Theorem 2 is provided in the Supplementary Material. The intuition for the
presented rates of convergence can be described, informally, as follows. We bound the distance between the sample split point, b
t, and its population counterpart, t, by characterizing
bL,R (a), is
the behaviour of the sample criterion function near t. The sample criterion, G
the empirical analog of the population criterion, GL,R (a), and is defined as the difference
between the averages of the observations in [a, R] and [L, a], respectively. We examine the
decrease in GL,R (t) that occurs when t is perturbed by a small amount, δ. Then, we contrast
this decrease with the stochastic term given by the difference between the corresponding
p
bL,R and GL,R . The order of this term is roughly δ/n. When t is a populadeviations in G
tion split point and, thus, lies in the interior of (L, R), the corresponding decrease in GL,R
is quadratic in δ. Balancing out the two terms yields the cube root asymptotic behaviour
(cf. Kim and Pollard 1990) for the sample split point. When function GL,R is maximized
14

at an endpoint of the interval [L, R], as in the case of truncation, the decrease in GL,R is
p
linear in δ. Balancing this decrease with the stochastic term of order δ/n suggests that
the sample split point is a Op (n−1 ) amount away from the boundary of (L, R). Uniformity
of the rate over all the small sample splits requires an additional log n factor.
In the next section we take advantage of our asymptotic results to propose a key modification to the sample clustering procedure.

3.3

Big Merge Tracker: Post-processing the Sample Procedure

Theorem 1 demonstrates that the sample analog of the truncation operation is peeling a
large number of tiny clusters off the ends of a large cluster. It follows that when recording
sample splits we should distinguish between those that correspond to splits in the population procedure and those that correspond to truncations. Based on this observation, we
propose to post-process the sample clustering procedure by only keeping the splits with
significant empirical sizes. More specifically, given a threshold α, if the cardinality of one

0.0

0.2

0.4

0.6

0.8

1.0

of the sub-clusters is below αn, the corresponding split is removed from the final output.

0.0

0.2

0.4

0.6

0.8

1.0

−4

−2

0

2

4

6

−4

−2

0

2

4

6

Figure 2: The plots illustrate the path of Algorithm 1 on a sample of 1000 observations
from a symmetric Gaussian mixture, 0.5 N (−2, 1)+0.5 N (2, 1). The scatter plot on the left
displays the sample frequencies for each pair of clusters merged along the path. The next
two plots show the cluster memberships before and after the big merge. For the leftmost
plot, the x and y axes denote the proportions of the points in the two merging clusters. For
the other two plots, only the x axis, which marks the locations of the points, is informative.
Figure 2 illustrates the path of Algorithm 1 on a sample of 1000 observations, generated
independently from the symmetric Gaussian mixture distribution used in Figure 1. The
scatter plot on the left displays the sample frequencies for each pair of clusters merged
along the path. We found only one merge in which both clusters pass the α = 0.1 threshold.
The big merge occurs at a point where the current number of clusters is 32. The two
15

rightmost plots display cluster memberships before and after the merge. The non-shaded
points belong to clusters with non-appreciable size.
The equivalence between the splitting procedure and the merging algorithm implies that
in the post-processing step of Algorithm 1 we only keep the merges with the cardinality of
each of the merging clusters above αn. For any such merge, we place the split point midway
between the two closest representatives of the two clusters being merged. We also replace
the stored merges with the corresponding split points. The resulting sequence of split points
can then be reinterpreted as a sequence of splits, or a sequence of merges, using the full
sample. For example, if the final output contains no split points, then all of the observations
in the sample are placed in the same cluster. We call this modified approach the Big Merge
Tracker (BMT) with threshold α. As a direct consequence of Theorems 1 and 2, under the
respective regularity conditions, the BMT consistently estimates the number of population
clusters, and its split points converge to their population counterparts at the Op (n−1/3 ) rate.
In the next section we analyze the empirical performance of the BMT approach.

4

Simulation Study and Real Data Analysis

In this section we show strong performance of the proposed BMT approach relative to
popular existing methods for assessing modality and estimating the number of clusters.
We also illustrate the use of our methodology in analysis of single-cell virology datasets.
In addition, in Section 5.3 of the Supplementary Material we apply BMT on very large
simulated data sets and demonstrate its superior scalability properties.
When the separation between two clusters is very small, the population splitting procedure can still be successful at finding a split point by massively truncating the support.
This zooming-in effect may result in larger sizes of the small sample splits, as we discussed,
from a theoretical standpoint, in the paragraph above the statement of Theorem 2. To counteract this phenomenon, we propose an adjustment to the Big Merge Tracker. If the sum of
the sample frequencies for the two merging clusters in the last big merge is less than 50%,
we do not report any merges. Preventing the corresponding splitting procedure from truncating more than 50% of the data, while searching for the first split, slightly reduces its
efficiency, but makes it more robust to sampling fluctuations. Throughout this section we
use the adjustment described above and set the BMT threshold, α, equal to 0.1 (Algorithm 1

16

in Section 5 of the Supplementary Material contains the full pseudocode). Note that a large
number of additional simulation results, corresponding to a wider range of sample sizes,
together with an analysis regarding the choice of α, are provided in Section 5.4 of the
Supplementary Material.

4.1

Modality Assessment

Testing for homogeneity of a population is an important statistical problem (Aitkin and
Rubin, 1985; Müller and Sawitzki, 1991; Roeder, 1994). Here, we use the BMT to detect
the presence of two or more dominant modes in the density. In Table 1 we compare our approach with two popular modality assessment procedures: (i) kernel density estimate based
test of Silverman (1981) (ii) histogram based Diptest proposed by Hartigan and Hartigan
(1985). P-values of the Silverman test are calculated using the R-package referenced in
Vollmer et al. (2013). R-packge of Maechler (2013) is used for implementing the Dip test.
Detailed descriptions of these procedures are given in Section 5.1 of the Supplementary
Material.
We consider 5 different simulation scenarios, in which 100 independent samples of
size 10000 were generated and subjected to modality analysis. Table 1 reports the percentage of cases in which multi-modality was detected. P-values for the Dip and Silverman
tests were computed based on 1000 MCMC simulations, and decision on the null hypothesis of unimodality was made at the 5% level of significance. The mean and the standard
deviation of the p-values from these tests are also reported. In the two unimodal scenarios
the BMT is on par with the Silverman and the Dip tests in confirming unimodality of the
population distribution with high certainty. In the three non-unimodal cases, which include
normal and beta mixtures, the BMT shows better performance in detecting multi-modality.

4.2

Estimating the Number of Clusters

We study the potency of the BMT in detecting the true number of clusters. We compare
its performance with the following number of clusters estimation methods: (i) the CH index of Caliński and Harabasz (1974) (ii) the KL index of Krzanowski and Lai (1988) (iii)
the H measure of Hartigan (1975) (iv) the Silhouette statistic based KR index of Kaufman and Rousseeuw (2009), (v) the Gap statistic of Tibshirani et al. (2001) (vi) the Jump

17

Population Density

Silverman Test P-value (S)

BMT

Std(D)

% multi-mode

Mean (S)

Std(S)

% multi-mode

% multi-mode

N(0,1)

0.99

0.04

0.00

0.48

0.25

0.00

0.00

Beta(2,4)

0.98

0.04

0.00

0.54

0.28

2.00

0.00

N (−1.1, 1) + N (1.1, 1) /2

0.81

0.22

0.00

0.22

0.21

29.00

69.00

Beta(4, 6) + Beta(7, 3) /2

0.84

0.22

0.00

0.31

0.25

21.00

49.00

N (−2.5, 1) + N (0, 1) + N (2.5, 1) /3

0.10

0.14

52.00

0.03

0.03

79.00

96.00





Dip Test P-value (D)
Mean (D)

Table 1: Simulation study to compare multi-modality detection methods
statistic of Sugar and James (2003) (vii) the clustering prediction strength criterion of Tibshirani and Walther (2005), and (viii) the bootstrap based cluster instability minimizing
criterion of Fang and Wang (2012), which is inspired by the stable clusters selection approach of Wang (2010a). Detailed descriptions of these procedures are provided in Section 5.2 of the Supplementary Material. We consider one multivariate and five univariate
regimes. 100 independent replications with the sample size of 5000 are used in each simulation setting, and the distribution of the number of clusters detected by each method
is reported. The eight competitor methods are implemented using a number of different clustering approaches via the NbClust R-package of Charrad et al. (2014) and the
fpc package of Hennig (2014). More specifically, we use k-means clustering with the Euclidean distance metric (the corresponding results are reported in Table 2), as well as Ward’s
method (Ward Jr, 1963), Centroid-based clustering (Kaufman and Rousseeuw, 2009), PAM
(Kaufman and Rousseeuw, 1990), CLARA (Kaufman and Rousseeuw, 1986), clustering by
merging Gaussian mixture components (Hennig, 2010) and hierarchical clustering initialized Gaussian mixture model based clustering method of Fraley and Raftery (2002). All
the results for the non-k-means clustering approaches are reported in the tables provided in
Section 5.2 of the Supplementary Material.
In our first univariate example, we consider a non-symmetric mixture of two normal
densities. Each of them has unit variance and their means are fairly well-separated. We
observe that the CH, KL and Hartigan methods struggle to recover the bimodal structure,
while the others successfully detect the two clusters (we note that in this setting, the CH
index performs better when it uses the centroid based clustering algorithm). The next
three simulation scenarios correspond to non-symmetric tri-modal population densities that
18

True Population Density

Methods

0.3 N (−4, 1) + 0.7 N (4, 1)

0.3 N (−3, 1) + 0.35 N (0, 1) + 0.35 N (3, 1)

0.3 t1 (−3) + 0.35 t1 (0) + 0.35 t1 (3)

0.3 dexp(−3) + 0.35 dexp(0) + 0.35 dexp(3)



Beta(8, 2) + Beta(5, 5) + Beta(2, 8)

{0.5 N (−2, 1) + 0.5 N (2, 1)}
⊗ N (0, 1)
⊗ N (0, 1)
⊗ χ21
⊗ χ21


3

CH
KL
Hartigan
Silhouette
Gap
Jump
Pred Str.
Stability
BMT
CH
KL
Hartigan
Silhouette
Gap
Jump
Pred Str.
Stability
BMT
CH
KL
Hartigan
Silhouette
Gap
Jump
Pred Str.
Stability
BMT
CH
KL
Hartigan
Silhouette
Gap
Jump
Pred Str.
Stability
BMT
CH
KL
Hartigan
Silhouette
Gap
Jump
Pred Str.
Stability
BMT
CH
KL
Hartigan
Silhouette
Gap
Jump
Pred Str.
Stability
BMT

1
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
17
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
99
0
0
0
0
0
0
0
0
0
0
0
0

2
29
39
0
100
100
99
100
100
93
0
11
0
8
100
0
0
11
5
7
18
0
71
36
11
24
69
1
0
13
0
39
100
0
0
0
0
0
13
0
0
100
0
78
0
0
0
18
0
0
100
14
100
85
96

3
12
10
32
0
0
0
0
0
7
0
8
69
92
0
100
100
89
95
9
23
46
26
64
9
53
30
99
0
8
55
61
0
100
100
100
100
2
3
57
100
0
1
17
25
100
1
11
52
87
0
0
0
0
4

Number of Clusters
4
5
6
7
8
4
2
2
4 11
5
5
8
8
8
16 12 10 11
5
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
1
0
6 13
11
9
6 14 12
13
7
4
2
4
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
7 11
7
6 10
15 14
5
4
9
24 11
5
4
5
3
0
0
0
0
0
0
0
0
0
13 12
5 11
9
19
4
0
0
0
0
1
0
0
0
0
0
0
0
0
0
1
0
6 14
18
9 11
7 20
19
9 10
3
1
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
3
3
6 16
11
9 10 15 16
14
9
6
5
5
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
4
1
0
0
0
60 15
0
0
0
0
0
0
0
0
0
4
5 11 10
14
7 13
9 11
16 13
4
7
4
12
0
1
0
0
0
0
0
0
0
62
0
0
0
0
0
0
0
0
0
15
0
0
0
0
0
0
0
0
0

9
14
9
7
0
0
0
0
0
0
28
19
1
0
0
0
0
0
0
14
6
1
0
0
7
0
0
0
23
10
1
0
0
0
0
0
0
22
9
0
0
0
0
0
0
0
23
11
4
0
0
8
0
0
0

Table 2: Number of clusters detected in 100 trials for six simulation scenarios
19

10+
22
8
7
0
0
1
0
0
0
52
10
0
0
0
0
0
0
0
29
6
4
0
0
6
0
0
0
56
4
2
0
0
0
0
0
0
48
14
4
0
0
0
0
0
0
46
6
0
0
0
16
0
0
0

are mixtures of standard normals, non-central t-densities with one degree of freedom and
double exponential densities with the unit rate parameter, respectively. The medians of the
mixing densities and the mixture weights match across the three settings, and the separation
between the adjacent medians is not large. In these three simulation settings, the Gap
statistic approach, as well as the CH, KL and Hartigan measures, has difficulty detecting
the true number of clusters. The Silhouette method performs well for Gaussian mixtures
but has difficulties in the other two cases. Jump statistic, prediction strength and bootstrap
stability approaches do well in the normal and the double exponential cases, however, they
do not show good performance in the considerably thicker-tailed case of the mixture of
t-densities. For our fifth simulation setting we consider a bounded population density that
is a mixture of three Beta densities. Here, only the BMT and the Silhouette do well in
recovering the true number of clusters. In our last example we consider a 5 dimensional
data set, which is generated from a product density. The first dimension is generated from
a symmetric mixture of two Gaussians; the next two dimensions contain white noise, while
the forth and fifth dimensions are generated from a central chi-square distribution with one
degree of freedom. We observe that, together with the CH, KL and Hartigan measures, the
Silhouette and the Jump approaches do not perform well in detecting the two clusters in
this data.
BMT does consistently well across each of the six simulation scenarios, outperforming
all the other approaches overall. The prediction strength and the cluster stability methods,
which are modern state of the art approaches, deliver the best results among the competitors. However, these two methods have significant trouble in the cases of the beta and
the non-central t mixtures. When implemented with various non-k-means clustering approaches (see Section 5.2 of the Supplementary Material), neither of the competitors considerably improves the performance reported in Table 2. Figure 3 in Section 5.2 of the
Supplementary Material provides plots of the densities used in the numerical experiments.

4.3

Sub-population Analysis in Single Cell Virology

We demonstrate an application of our clustering approach in a single-cell Mass Cytometry
(Bendall et al., 2011) based virology study. We analyze the data reported in Sen et al.
(2014), where the effect of Varicella Zoster Virus (VZV) on human tonsil T cell is studied. VZV is a human herpesvirus and causes varicella and zoster (Zerboni et al., 2014).
20

We study protein expressions from five independent experiments, each containing an Uninfected (UN) and a Bystander (BY) populations. Bystanders are cells in the VZV infected
population, which are not directly infected by the virus, but are influenced by neighboring virus infected cells. Protein expression values are studied on the arcsinh scale. Nonexpressed values are uniformly distributed between [−1, 1]. Cellular sub-populations are
detected by clustering the populations based on the expressions of “core-proteins”, which
are associated with T cell activation (Newell et al., 2012). Most of the samples have large
sizes, usually on the order of ∼ 105 . Traditional clustering techniques fail to accommodate
such large sample sizes and resort to sub-sampling based approaches (Qiu et al., 2011;
Linderman et al., 2012). The BMT, on the other hand, has the advantage of being scalable
enough to conduct clustering analysis on the entire sample.
SUB-POPULATIONS

Experiment I
UN
BY

Experiment II
UN
BY

Experiment III
UN
BY

Experiment IV
UN
BY

Experiment V
UN
BY

DUAL POSITIVE

8411
(8.8%)

6596
(7.3%)

5253
(5.8%)

4169
5.7%)

4971
(6.0%)

2703
(5.4%)

3795
(5.0%)

1510
(4.6%)

8047
(8.5%)

5225
(8.0%)

DUAL NEGATIVE

2723
(2.8%)

2973
(3.3%)

3537
(3.9%)

2631
(3.6%)

4433
(5.3%)

2935
(5.9%)

4354
(5.8%)

2196
(6.7%)

5012
(5.3%)

2881
(4.4%)

CD4 NON-NAIVE

7993
(8.4%)

10636
(11.8%)

15144
(16.7%)

11556
(15.9%)

21444
(25.9%)

12429
(25.0%)

22149
(29.7%)

8508
(26.1%)

30034
(31.9%)

20469
(31.3%)

CD4 NAIVE

69977
(73.7%)

64119
(71.1%)

57744
(63.7%)

47374
(65.1%)

45764
(55.3%)

27987
(56.3%)

35458
(47.5%)

16390
(50.4%)

40398
(43.0%)

28524
(43.7%)

CD8 NAIVE

5654
(6.0%)

5671
(6.3%)

8599
(9.5%)

6571
(9.0%)

5798
(7.0%)

3490
(7.0%)

8271
(11.1%)

3774
(11.6%)

9869
(10.5%)

7829
(12.0%)

POPULATION SIZE

94837

90157

90641

72699

82637

49672

74540

32497

93878

65244

Table 3: Sizes and Proportions of dominant clusters detected by BMT across 5 independent
Virology experiments
We treat three proteins, CD4, CD8 and CD45RA (naive), as core-proteins, as they are
typically used to classify T cells. For each of the 10 samples (UN and BY from experiments I-V), based on the expressions of the above three proteins, we performed automated
clustering by using BMT in the three dimensional space. Figure 4 and 5 in Section 5.5
of the Supplementary Material show that in all the cases, the BMT detects unimodality
for CD4 and CD45RA and bimodality for CD8 expression values. Using the bi-modality
of CD8 and the BMT detected splits, we classify cells as CD8-high and CD8-low. Also,
considering the expression and non-expressions of the other two markers we simultane21

ously classify cells into the following clusters, or sub-populations: (i) Dual positive: CD4
expressed and CD8 High (ii) Dual-negative: CD4 non-expressed and CD8 low (iii) CD4
Non-Naive: CD4 expressed and CD45RA non-expressed (iv) CD4 Naive: CD4 expressed
and CD45RA non-expressed (v) CD8 Naive: CD8 high and CD45RA expressed.
Table 3 reports the sizes and proportional representations of these sub-populations,
across the five experiments. BMT sub-populations resemble the T-cell biology based phenotypic classification in Sen et al. (2014). They also revalidate that the sub-population
distribution in the Bystander cells is not much different from that of the uninfected, though
the UN sub-population distribution varies across experiments. Using this BMT based categorization of the T cells, sub-population level cell-signaling patterns can be subsequently
studied. Figure 6 in Section 5.5 of the Supplementary Material shows the heatmap of the
protein expressions (core + signaling proteins) of the sub-populations from Experiment I.

5

Discussion

In this paper we focus on the analysis of a popular convex clustering approach that is based
on the `1 fusion regularization. However, we believe that our general theoretical framework can be extended to handle other types of fusion penalties. In particular, based on a
preliminary analysis of the corresponding `2 approach, we conjecture that the asymptotic
results in this paper also hold in the `2 case, under appropriate regularity conditions. The
corresponding population procedure can be defined by analogy, as a collection of cluster
splits and truncations, where each operation seeks to maximize the Euclidian distance between the corresponding sub-cluster means. The proofs require a rigorous formulation and
a thorough analysis of a multivariate analog of the truncation operation.
High computational efficiency of the proposed BMT approach allows it to be applied to
massive high-dimensional data sets. In particular, BMT can be used for high-dimensional
feature screening, to rule out predictors that do not reveal any clusters in the data. Moreover,
we can take further advantage of BMT’s computational efficiency and apply the screening
procedure to several linear combinations for each pair of variables. This way we can move
beyond marginal screening, similarly to how regression models with interaction terms move
beyond the simple additive structure.

22

6

Supplementary Material

6.1

Proof of Theorem 1

Preliminaries
Here we provide three lemmas that make important contributions to the proof. The first
allows us to focus on a compact support, the second established uniform convergence of
the sample criterion functions on compact sets, and the third derives useful continuity properties of the population solution. The proofs of the lemmas are provided below, after the
main argument for Theorem 1. All of the results take advantage of regularity conditions
C1 and C2, which are given in Section 2.3 of the main paper. In the proof of the lemmas
we verify that these conditions can be removed for densities f that are unimodal and continuous. First, we handle the case where the support of density is unbounded. By analogy
with the term population cluster (page 13 of the main paper), we use sample cluster to refer
to intervals that appear along the path of the sample clustering procedure. We follow the
traditional approach of using ω to denote a general element of the sample space.
Lemma 1. For every positive , δ and δ 0 , and almost all ω, there exist a sample cluster
b R),
b a population cluster (L, R), and a bounded interval B , which depends only on ,
(L,
such that, for all sufficiently large n,
≤
(a) P−∞,Lb ≤  and PR,∞
b
b − R| ≤ δ and |L
b − L| ≤ δ 0
(b) |R
b R)
b ⊆ B .
(c) (L,
Remark.

b and R
b depend on ω, while B does not.
Note that L, R, L

To draw connections between the sample and the population clustering procedures,
bL,R , as the empirical analog of the population
we define the sample criterion function, G
criterion, GL,R . More formally, we write µ
bl,r for the the sample average on [l, r] and let
bL,R (a) = µ
G
ba,R − µ
bL,a ,

(11)

for a ∈ [L, R], with the following important caveat. If either of the intervals [L, a] and
bL,R (a) to be undefined. We now summarize
[a, R] contains no observations, we consider G
23

the sample clustering procedure using the newly defined criterion function. Given a cluster
bL,R , then splits (L, R) into
(L, R), the sample procedure finds a point s that maximizes G
sub-clusters (L, s) and (s, R). Note that s is guaranteed to be located strictly between two
of the observations.
To understand the behavior of the sample clustering procedure on a bounded support,
we establish uniform convergence of the of the sample criterion functions to their population counterparts.
Lemma 2. Define AM = {(L, R), −M ≤ L < R ≤ M }. Then, as n → ∞,
sup
(L,R)∈AM

bL,R (a) − GL,R (a)| → 0,
max |G

a∈[L,R]

(12)

almost surely, for each positive M .
Remark.

b L,R (a) the maximum is taken over the
In the expression maxa∈[L,R] G

closed interval [ min{xn ∩ [L, R]} , max{xn ∩ [L, R]} ]. By the classical GlivenkoCantelli theorem, the endpoints of this interval converge to L and R almost surely, and
the convergence is uniform over (L, R) in AM , for each fixed M .

The next result derives some important properties of the population procedure, which
we then relate to the sample procedure with the help of Lemma 2. To state the result,
we need to modify RL , which is a function of L that appears in the definition of the two
sided truncation. Suppose that the population procedure truncates interval (Ls , Rs ) to interval (L∗ , R∗ ) in a two-sided fashion. We will extend function RL to be defined in open
neighborhoods of Ls and L∗ , using the equation GL,RL (L) = GL,RL (RL ). The proof of
Proposition 4 shows that, provided the neighborhoods are small enough, the extension is
uniquely determined if we require that function RL remains continuous and decreasing. To
simplify the notation we write max GL,R for maxa∈[L,R] GL,R (a).
Lemma 3. Let (Ls , Rs ) be a bounded cluster, produced by the population procedure, such
that arg max GLs ,Rs = {Ls , Rs }. Suppose that the first split of the population procedure
applied to (Ls , Rs ) is given by (L∗ , s∗ , R∗ ). For every sufficiently small positive  there
exist positive 0 , δ and δ 0 , such that 0 ≤ , δ 0 ≤ δ ≤ /2 and |RL − RL+δ0 | ≤ δ for all L

24

in [Ls − δ 0 , L∗ + 0 ], and the following sets of inequalities hold as long as |R − RL | ≤ 2δ:
max GL,R >

max

[L+δ 0 ,R−δ]

max GL,R >

max

[L+δ 0 ,R−δ]\(s∗ −,s∗ +)

max GL,R >

max GL,R

>

[R−δ,R]

max

[L,R]\(s∗ −,s∗ +)

max GL,R >
Remark.

for L ∈ [Ls − δ 0 , L∗ − 0 ],

GL,R

for L ∈ [L∗ − 0 , L∗ + 0 /2],

GL,R

for L ∈ [L∗ + 0 /2, L∗ + 0 ].

max0 GL,R

for L ∈ [Ls − δ 0 , L∗ + 0 ], R ∈ [RL + δ, RL + 2δ],

max GL,R

for L ∈ [Ls − δ 0 , L∗ + 0 ], R ∈ [RL − 2δ, RL − δ].

[L,L+δ ]

[L,L+δ 0 ]

GL,R

[R−δ,R]

A minor modification is required when Ls = L0 or Rs = R0 , where

(L0 , R0 ) is the support of f : each interval appearing in the statement needs to be
replaced by its intersection with [L0 , R0 ].

b in the five sets of inWe will now use Lemma 2 to show that if we replace G with G
equalities above, then for each ω in a set of probability one the resulting sample inequalities simultaneously hold for all sufficiently large n. For example, the first set becomes
bL,R > max G
bL,R . We also allow the cluster (Ls , Rs ) to depend on ω, unmax G
[L+δ,R−δ]

der an additional assumption that there exists a finite deterministic M , for which bound
max{|Ls − δ 0 |, |Rs + δ|} ≤ M holds almost surely. To establish the sample inequalities,
bL,R − GL,R | and
we define ∆n = supA max |G
M

τ =

min

L∈[Ls −δ 0 ,L∗ −0 ], |R−RL |≤2δ


max GL,R −


max

[L+δ 0 ,R−δ]

GL,R .

Continuity of GL,R (a) with respect to a, L an R, together with the first set of inequalities
in Lemma 3, implies τ > 0. Taking into account Lemma 2, as well as the continuity
of GL,R (a), we derive the following inequalities, which hold for all ω in a set of probability
one, and all sufficiently large n:
bL,R − max G
bL,R > max GL,R − max GL,R − 2∆n − o(1) ≥ τ − o(1) > 0.
max G
[L+δ,R−δ]

[L+δ,R−δ]

The argument for the rest of the inequalities is analogous.
25

Main Body of the Proof
We restrict our attention to the set of probability one that remains after we cast out the negligible sets during the application of the almost sure results from Lemma 1 and Lemma 2.
We focus on a fixed ω, but suppress the dependence on ω to simplify the notation.
Consider an arbitrary positive . We will investigate the behavior of the sample probs , R
bs ), that is close to
cedure, starting with an appropriately bounded sample cluster, (L
bs , R
bs ) ⊆ B , where B
some population cluster, (Ls , Rs ). More specifically, we assume (L
bs − Rs | ≤ δ and
is a bounded interval that does not depend on ω. We also assume |R
bs − Ls | ≤ δ 0 , for all sufficiently large n, where δ and δ 0 are to be specified later. We will
|L
track the sample and population procedures up to the first population split and show, under
only the above assumptions, that the two procedures are close in an appropriate sense. We
need to consider the following options for the population procedure applied to (Ls , Rs ).
First: two-sided truncation, then a split; second: one-sided truncation, then a two sided
truncation, followed by a split; third: one-sided truncation, followed by a split. There are
three more possibilities, where in each of the initial options the procedure does not produce
a split. Finally, the cluster (Ls , Rs ) may be split right away.
We first focus on the important case of a two-sided truncation, followed by a split.
The rest of the cases can be handled with only minor modifications to the argument. Let
sp = (L∗ , s∗ , R∗ ) be the first split of the population procedure applied to (Ls , Rs ). The
population procedure performs a two-sided truncation along (L, RL ), where RL is a continuous decreasing function of L, such that RLs = Rs and RL∗ = R∗ . Given our values
of , Ls , Rs and sp, we let 0 , δ and δ 0 be the quantities from Lemma 3. Recall that 0 ≤ 
and δ 0 ≤ δ ≤ /2. As we demonstrate in the paragraph below the statement of Lemma 3,
the five sets of population inequalities given in the lemma have natural sample counterparts.
The first three sets of inequalities hold for all R with |R − RL | ≤ 2δ and all sufficiently

26

large n:
bL,R >
max G

bL,R >
max G

bL,R >
max G

max

[L+δ 0 ,R−δ]

max

[L+δ 0 ,R−δ]\(s∗ −,s∗ +)

max

[L,R]\(s∗ −,s∗ +)

bL,R
G

for L ∈ [Ls − δ 0 , L∗ − 0 ]

bL,R
G

for L ∈ [L∗ − 0 , L∗ + 0 /2]

bL,R
G

for L ∈ [L∗ + 0 /2, L∗ + 0 ].

(13)

b R)
b for the current cluster in the sample procedure.
To simplify the notation, we write (L,
b L
b + δ 0 ] or [R
b − δ, R]
b and big if
We call a sample split small if the split point is placed in [L,
the split point is placed in [s∗ − , s∗ + ]. Inequalities (13) directly imply that the following
b − R b | ≤ 2δ is satisfied:
properties hold for all sufficiently large n, provided |R
L

b ≤ L∗ − 0 , then the next split of the sample procedure is small;
(i) if L
b ∈ [L∗ − 0 , L∗ + 0 /2], then the next sample split is either small or big;
(ii) if L
b ∈ [L∗ + 0 /2, L∗ + 0 ], then the next sample split is big.
(iii) if L
b
We refer to |R−R
b | ≤ 2δ as the requirement. Note that it is satisfied for the starting cluster,
L
bs , R
bs ), by assumption, and because |RL − RL+δ0 | ≤ δ for all L in [Ls − δ 0 , L∗ + 0 ], by
(L
Lemma 3. We will verify that the requirement remains valid until the sample procedure
b − R b | ≤ δ. By properties (i)-(ii), together with the
makes a big split. First, suppose |R
L

b b
bound |RLb − RL+δ
b 0 | ≤ δ, either the next sample split is big or it reduces (L, R) to a cluster
b − R b | ∈ [δ, 2δ]. Consider the
that still satisfies the requirement. Now consider the case |R
L

last two inequalities in Lemma 3, and, again, derive the sample sample counterparts,
bL,R
max G

>

[R−δ,R]

bL,R >
max G

[L,L+δ 0 ]

bL,R
max G

for L ∈ [Ls − δ 0 , L∗ + 0 ], R ∈ [RL + δ 0 , RL + 2δ]

bL,R
max G

for L ∈ [Ls − δ 0 , L∗ + 0 ], R ∈ [RL − 2δ 0 , RL − δ].

[L,L+δ 0 ]

[R−δ,R]

Taking inequality |RLb − RL+δ
b 0 | ≤ δ into account, we deduce that the next sample split,
b − R b |. Thus, when the sample procedure is applied
if small, reduces the distance |R
L
b
b
to (Ls , Rs ), the requirement remains satisfied until the big split, and properties (i)-(iii)
remain valid. Consequently, we can use these properties repeatedly to establish that, for
27

bs , R
bs ) makes a number of small
all sufficiently large n, the sample procedure starting at (L
b∗ , sb∗ , R
b∗ ), such that ||sp
splits, followed by a big split, sp
b = (L
b − sp||∞ < .
Now consider the case where the population procedure truncates, in a two-sided fashion, the cluster (Ls , Rs ) all the way to an empty set. To establish that the sample procedure,
bs , R
bs ), makes only small splits, we can apply the same reasoning as in the preapplied to (L
vious case, but using only the first of the three inequality sets in display (13) and only the
first of the properties (i)-(iii). The rest of the cases can be handled analogously, with only
minor modifications to the original argument. Thus, for all sufficiently large n, the sample
bs , R
bs ) either makes small splits until the end or makes small splits
procedure starting at (L
until a split sp,
b such that ||sp
b − sp||∞ < . Note that the width of the smaller sub-cluster
produced by a small split is less than . Because  can be chosen arbitrarily small, we conclude that sp
b converges to sp, while the maximum size of the splits prior to sp
b goes to zero.
Note that size(sp)
b converges to size(sp) = α∗ , by continuity of the function size. Then,
c sp)
by the classical Glivenko-Cantelli theorem, size(
b converges to α∗ as well.
Because the number of population splits is finite, we can repeat all of the preceding
arguments sequentially. For example, to analyze the clustering procedures between the
first and, potentially, second population split, we first consider the left sub-cluster and set
bs and R
bs to L∗ , s∗ , L
b∗ and sb∗ , respectively. Note that (Ls , Rs ) is appropriately
Ls , Rs , L
bs − Rs | ≤ δ and |L
b s − Ls | ≤ δ 0
bounded, and, for every positive δ and δ 0 , inequalities |R
hold for all sufficiently large n. This verifies the assumptions stated in the first paragraph of
the proof, which are the only ones needed for the preceding argument. The same is true for
bs and R
bs to s∗ , R∗ , sb∗ and R
b∗ , respectively.
the right sub-cluster, on which we set Ls , Rs , L
To complete the proof of the theorem, it is only left to show that the maximum size of all
bs , R
bs ), converges to zero. Recall
the sample splits leading to the very first starting cluster, (L
bs , R
bs ), with the required properties, is guaranteed by Lemma 1. The
that the existence of (L
same result also tells us that P−∞,Lbs and PRbs ,∞ are bounded above by . As  can be chosen
arbitrarily small, we have established the desired convergence.
Proof of Lemma 1
The proof takes advantage of several applications of Lemma 2 and the strong law of large
numbers. We restrict our attention to the set of probability one that remains after casting out
the negligible sets associated with the aforementioned convergence results. We conduct a

28

pointwise argument, at a fixed ω, but suppress the dependence on ω to simplify the notation.
When (L0 , R0 ), the support of the underlying distribution, is bounded, the conclusion
of the lemma holds for L = L0 and R = R0 . Consider an unbounded case of R0 = ∞
and |L0 | < ∞. Note that the population procedure starts with a right-sided truncation.
Let T ∗ denote the right endpoint of the cluster at which the population procedure either
transitions to a two-sided truncation or makes a split. Fix a K1 , such that PK1 ,∞ ≤ ,
with a further requirement that K1 > T ∗ , if T ∗ exists. Take K2 large enough to satisfy
K2 > µL0 ,∞ + µK1 ,∞ − L0 . By the law of large numbers, K2 − µ
bL0 ,∞ > µ
bK1 ,∞ − L0 ,
for all sufficiently large n. Suppose that n is large enough for the above lower bound to be
satisfied. Then, the following inequalities hold for all R > K2 and r ∈ [L0 , K1 ]:
bL0 ,R (R) = R − µ
bL0 ,R (r).
G
bL0 ,R > K2 − µ
bL0 ,∞ > µ
bK1 ,∞ − L0 > µ
br,R − µ
bL0 ,r = G
b with R
b > K2 , the next sample split
This guarantees that for every sample cluster (L0 , R),
point is placed to the right of K1 , which implies that the sample procedure will eventually
b with R
b ∈ [K1 , K2 ]. Because the population procedure reduces
produce a cluster (L0 , R)
(L0 , ∞) to (L0 , K1 ) using the right-sided truncation, interval (L0 , R) is a bounded population cluster for each R in [K1 , K2 ]. This completes the proof for the case R0 = ∞,
|L0 | < ∞. The case L0 = −∞, R < ∞ can be handled analogously.
Now consider the case (L0 , R0 ) = (−∞, ∞). Let K1− and K1+ satisfy P−∞,K1− ≤ 
and PK1+ ,∞ ≤ . If the population procedure makes a split, and the first one is applied
to the cluster (L∗ , R∗ ), we further require (K1− , K1+ ) ⊃ [L∗ , R∗ ]. Additional conditions
are placed on K1− and K1+ below. Take K2+ sufficiently large to ensure K2+ > µK1− ,∞ +
µK1+ ,∞ − µ−∞,K1− . By the law of large numbers, K2 − µ
bK1− ,∞ > µ
bK1+ ,∞ − µ
b−∞,K1− for all
sufficiently large n. Consequently, for all r ∈ [K1− , K1+ ], R > K2+ , L ≤ K1− we have:
bL,R (R) = R − µ
bL,R (r).
G
bL,R > K2+ − µ
bK1− ,∞ > µ
bK1+ ,∞ − µ
b−∞,K1− ≥ µ
br,R − µ
bL,r = G
b R),
b with L
b ≤ K1− and R
b > K2+ , the next
This guarantees that for every sample cluster (L,
sample split point is placed outside of [K1− , K1+ ].
bL,R (L) >
If K2− lies below µ−∞,K1+ + µ−∞,K1− − µK1+ ,∞ , then, by similar arguments, G
bL,R (r) for all r ∈ [K1− , K1+ ], L < K2− and R ≥ K1+ . This guarantees that for every
G
b R),
b with L
b < K2− and R
b ≥ K1+ , the next sample split point is placed
sample cluster (L,
29

outside of [K1− , K1+ ].
b R)
b with L
b ∈ [K2− , K1− ] and
Consequently, the sample procedure produces a cluster (L,
b ∈ [K1+ , K2+ ]. Because K1− < L∗ , interval (L,
b R b ) is a population cluster. It is only left
R
L
b
to show that we can find an R that is located within δ of R b . We will use the following
L

result, which is proved further below.
Lemma 4. Let (L, RL ) be a population cluster, achieved via a two-sided truncation of the
real line. If K1− is sufficiently small, then for each positive δ we have
for all L ∈ [K2− , K1− ] and R ∈ [RL + δ, K2+ ].

max GL,R > max GL,R

[RL +δ,R]

[L,RL ]

Applying the argument in the paragraph below the statement of Lemma 3, and replacb b for all R ∈ [R b +δ, K2+ ],
b we deduce that max[R +δ,R] G
b b > max b
G
ing L with L,
b
L

L,R

[L,RL
b]

L

L,R

b > R b + δ, the sample procedure sequentially
as long as n is sufficiently large. Hence, if R
L
b
b
moves the right endpoint of cluster (L, R) further left, until it falls in [R b , R b + δ]. The
L

L

b < R b − δ can be handled using analogous arguments.
case R
L
Proof of Lemma 2
Let Pn denote the empirical measure associated with observations x1 , ..., xn , and let P
be the corresponding population distribution. To simplify the presentation, we replace
R
R
expressions h(x)dPn (x) and h(x)dP (x) by Pn h and P h, respectively.
The classical Glivenko-Cantelli theorem gives the following uniform convergence,
sup |Pn (l, r) − P (l, r)| → 0,

(14)

l<r

almost surely, as n goes to infinity. Note that the collection of functions hl,r (x) = x1{l<x<r} ,
defined for all real l and r, forms a VC class with an integrable envelope, H(x) = |x|. Consequently, by a functional generalization of the Glivenko-Cantelli theorem (e.g. Theorem
19.13 in van der Vaart 1998), we have
sup |Pn hl,r − P hl,r | → 0,

(15)

l,r

almost surely, as n goes to infinity. For the rest of the proof we cast out the negligible sets
on which (14) and (15) break down, and conduct a pointwise argument on the remaining
30

set of probability one. We suppress the dependence on ω to simplify the notation.
Fix an arbitrarily small positive . Define AM () = AM ∩ {(l, r) : r − l ≥ }. Let
c = inf AM () P (l, r), and note that c is positive. Taking advantage of the convergence
in (14) and (15), we can bound supAM () |b
µl,r − µl,r | above by
sup
AM ()

Pn hl,r − P hl,r
1
1
−
+ sup |Pn hl,r |
P (l, r)
Pn (l, r) P (l, r)
AM ()

−1
≤ c−1
 sup |Pn hl,r − P hl,r | + c
l,r

−1
= c−1
 o(1) + c

supl,r Pn |hl,r |
sup |Pn (l, r) − P (l, r)|
inf AM () Pn (l, r) l<r

P H + o(1)
o(1) = o(1),
c + o(1)

as n tends to infinity. Consequently,
sup
(L,R)∈AM

bL,R (a) − GL,R (a)| ≤ 2 + 2 sup |b
max |G
µl,r − µl,r | = 2 + o(1).

a∈[L,R]

AM ()

This completes the proof, because  can be chosen arbitrarily small.
Proof of Lemma 3
For concreteness we assume Ls > L0 throughout the proof. The case Ls = L0 follows
with only minor modifications. Suppose that positive 0 and δ 0 are small enough for RL
to be defined on [Ls − δ 0 , L∗ + 0 ]. Then, maxL∈[Ls −δ0 ,L∗ +0 ] |RL − RL+δ0 | ≤ δ for every
positive δ, as long as δ 0 is sufficiently small, by uniform continuity of RL on a compact set.
Differentiating GL,R (s) with respect to s, we derive the following formulas:
G0L,R (L) =

f (L)(µL,R − L) 1
−
PL,R
2

and G0L,R =

1 f (R)(R − µL,R )
−
.
2
PL,R

(16)

Note that for L ∈ [Ls , L∗ ], function GL,R is maximized at the endpoints. Consequently,
G0L,RL (L) ≤ 0 and G0L,RL (RL ) ≥ 0, when L ∈ [Ls , L∗ ]. For the clarity of the exposition
we first focus on the regular case, where both inequalities are strict for L ∈ [Ls , L∗ ]. By
continuity of G0L,R , there exist positive δ1 and 1 , such that G0L,R (L) < 0 and G0L,R (R) > 0
for all L ∈ [Ls − δ1 , L∗ + 1 ] and |R − RL | ≤ δ1 . Using uniform continuity of G0L,R on

31

compact sets, we conclude that there exists a positive δ3 , such that
GL,R (L) ∨ GL,R (R) > GL,R (s)

for all s ∈ (L, L + δ3 ] ∪ [R − δ3 , R),

(17)

when L ∈ [Ls − δ1 , L∗ + 1 ] and |R − RL | ≤ δ1 . Because arg max GL,RL = {L, RL } for
L ∈ [Ls , L∗ ), continuity of GL,R implies that, given an arbitrarily small positive 0 , we can
find a positive δ2 , such that
GL,R (L) ∨ GL,R (R) >

max
[L+δ3 ,R−δ3 ]

GL,R ,

(18)

when L ∈ [Ls − δ2 , L∗ − 0 ] and |R − RL | ≤ δ2 . Inequalities (17) and (18) yield
GL,R (L) ∨ GL,R (R) > GL,R (s)

for all

s ∈ (L, R),

when L ∈ [Ls − δ1 ∧ δ2 , L∗ − 0 ] and |R − RL | ≤ δ1 ∧ δ2 . This implies the first set of
inequalities in Lemma 3.
By a similar argument, for every positive  there exist positive 2 and δ4 , such that
GL,R (L) ∨ GL,R (R) >

max

[L+δ3 ,R−δ3 ]\(s∗ −,s∗ +)

GL,R ,

(19)

when L ∈ [Ls −δ4 , L∗ +2 ] and |R−RL | ≤ δ4 . Combining this bound with inequality (17),
we derive the second set of inequalities in Lemma 3.
In the proof of Proposition 3 we establish inequalities µL∗ ,s∗ ≤ (L∗ + s∗ )/2 and
µs∗ ,R∗ ≥ (s∗ + R∗ )/2. In our case they are equalities, because GL∗ ,R∗ (s∗ ) = GL∗ ,R∗ (L∗ ) =
GL∗ ,R∗ (R∗ ). Taking this fact into account when solving G0L∗ ,R∗ (s∗ ) = 0 we derive a useful
equality, (R∗ − s∗ )/Ps∗ ,R∗ = (s∗ − L∗ )/PL∗ ,s∗ . It follows that
f (L∗ )(R∗ − s∗ )/Ps∗ ,R∗ < 1

and

f (L∗ )(s∗ − L∗ )/PL∗ ,s∗ < 1,

(20)

because otherwise f (L∗ )(R∗ − L∗ )/PL∗ ,R∗ ≥ 1, contradicting the fact that G0L,RL (L) < 0
for L ∈ [Ls , L∗ ], which holds in the regular setting that we now focus on.
Consider the function h(L) = GL,RL (s∗ ) − GL,RL (L), and note that h(L∗ ) = 0. Inequalities (20) guarantee that h0 (L∗ ) > 0. Hence, there exists a positive 3 , such that for

32

every 0 ≤ 3 ,
GL,RL (s∗ ) > GL,RL (L) = GL,RL (RL ),
when L ∈ [L∗ + 0 /2, L∗ + 0 ]. Continuity of GL,R (s), together with compactness of the
intervals involved, implies that for every 0 ≤ 3 there exists a positive δ5 , such that
GL,R (s∗ ) > GL,R (L) ∨ GL,R (R),

(21)

when L ∈ [L∗ +0 /2, L∗ +0 ] and |R−RL | ≤ δ5 . Because arg max GL∗ ,RL∗ = {L∗ , RL∗ , s∗ },
we can also find positive 4 and δ6 , for which
GL,R (s∗ ) >

max

[L+δ3 ,R−δ3 ]\(s∗ −,s∗ +)

GL,R ,

(22)

when L ∈ [L∗ , L∗ + 4 ] and |R − RL | ≤ δ6 . Combining bounds (17), (21) and (22), we
establish the third set of inequalities in Lemma 3.
Consider function hL (R) = GL,R (R) − GL,R (L). Differentiation gives h0L (RL ) =
1 − 2f (RL )(RL − µL,RL )/PL,RL , which is strictly positive for L ∈ [Ls , L∗ ] in the regular
case that we now focus on. Taking advantage of equality hL (RL ) = 0, continuity of h0L and
compactness of the intervals involved, we deduce that, provided 0 , δ and δ 0 are sufficiently
small, inequalities GL,R (R) > max[L,L+δ0 ] GL,R are satisfied for all L ∈ [Ls − δ 0 , L∗ + 0 ]
and R ∈ [RL + δ, RL + 2δ]. This establishes the fourth set of inequalities in Lemma 3. The
fifth set can be derived using analogous arguments.
We now examine the irregular setting, where G0L,RL (L) = 0 or G0L,RL (RL ) = 0 for
some L in [Ls , L∗ ]. For concreteness, we examine the case G0Li ,Ri (Ri ) = 0 6= G0Li ,Ri (Li ),
where Li ∈ [Ls , L∗ ] and Ri = RLi . Suppose that Ri < Rs ; the case Ri = R0 can be
handled with minor modifications. Differentiation gives G00Li ,Ri (Ri ) = −f 0 (Ri )/[6f (Ri )],
which, together with arg max GLi ,Ri ⊇ {Li , Ri }, implies f 0 (Ri ) ≥ 0. However, a strict
inequality, f 0 (Ri ) > 0 implies that 2f (RL )(RL − µL,RL ) > PL,RL for all L that are sufficiently close to Li and satisfy L < Li . Because this contradicts inequality G0L,RL (RL ) ≥ 0,
we conclude that f 0 (Ri ) = 0, and thus, G0Li ,Ri (Ri ) = G00Li ,Ri (Ri ) = 0. Further differeni
00
i
i
00
i
tiation yields G000
Li ,Ri (R ) = −f (R )/[2f (R )], which implies f (R ) ≤ 0. This stepwise

argument can be continued if f 00 (Ri ) = 0, however we will focus on the case f 00 (Ri ) < 0,
for concreteness. Note that Ri is an interior mode of f .
Let l = L − Li , t = s − Ri and r = R − Ri . After deriving the third order Taylor series
33

expansion for the function (L, s, R) 7→ GL,R (s) at (Li , Ri , Ri ), we establish the following
approximation, which holds for (l, t, r) near zero:
GLi +l,Ri +r (Ri + t) − GLi +l,Ri +r (Ri ) =

−f 00 (Ri ) 3
(3t − t2 r + 3tr2 ) + o(|l|3 + |t|3 + |r|3 ).
6f (Ri )

Analysis of the above expression reveals that, because f 00 (Ri ) < 0, the approximating
cubic function is increasing in t for every fixed r and l.
We now revisit the first set of inequalities in Lemma 3; the rest of the sets can be established using similar arguments. For the corresponding proof given in the regular setting to
still go through, it is sufficient to establish that for every small enough positive δ1 ,
GL,R (L) ∨ GL,R (R) > GL,R (s)

for all

s ∈ (L, R − δ1 /4],

when |L − Li | ≤ δ1 and |R − Ri | ≤ δ1 .
By the established monotonicity of the cubic approximation, we can find a positive δ3 ,
such that for each positive δ < 4δ3 ,
GL,R (R) > GL,R (s)

for all s ∈ [R − δ3 , R − δ1 /4],

(23)

when |L − Li | ≤ δ1 and |R − Ri | ≤ δ1 . Because arg max GLi ,Ri = {Li , Ri }, we can also
find a positive δ2 , such that
GL,R (R) ∨ GL,R (L) >

max
[L+δ3 ,R−δ3 ]

GL,R ,

(24)

when |L − Li | ≤ δ2 and |R − Ri | ≤ δ2 . Combining inequalities (23) and (24) establishes
GL,R (R) ∨ GL,R (L) >

max
[L+δ3 ,R−δ1 /4]

GL,R ,

when |L − Li | ≤ δ1 , |R − Ri | ≤ δ1 , and δ1 is sufficiently small. Using the regular case
argument in the beginning of the proof we also establish the left side bound,
GL,R (L) > GL,R (s)

for all s ∈ (L, L + δ3 ],

when |L − Li | ≤ δ2 , |R − Ri | ≤ δ2 , and δ2 is sufficiently small, which completes the proof.

34

Proof of Lemma 4
Let Rm be the maximum of R∗ and the right-most mode of f . Let Lm be the minimum of
L∗ and the left-most mode of f . As we point out in the proof of Proposition 3, any interior
maximizer of GL,R must lie in (Lm , Rm ). Thus, it is sufficient to verify
(i) GL,R (R) > max[Lm ,Rm ] GL,R

and

(ii) GL,R (R) > GL,R (L),
for all L ∈ [K2− , K1− ] and R ∈ [RL + δ, K2+ ]. If we choose K1− sufficiently small, so that
K1− < Lm , RK1− > Rm and RK1− > µ−∞,Rm + µRm ,∞ − µ−∞,Lm , then (i) follows from
GL,R (R) > RK1− − µ−∞,Rm > µRm ,∞ − µ−∞,Lm > max
GL,R .
m
m
[L ,R ]

The existence argument in the proof of Proposition 4 implies that that the derivative
of the function hL (R) = GL,R (R) − GL,R (L) is strictly positive on [RL , ∞), for every
L < Lm . Statement (ii) then follows from hL (RL ) = 0.

6.2

Proof of Theorem 2

Note that for the purposes of this proof we can replace “sample frequency” with “underb ). By the classical Glivenko-Cantelli
lying probability” in the definition of the set S(τ
theorem, the corresponding modified statement of the theorem implies the original one.
Write J (τ ) for the collection of all intervals (L, R), such that such that P (L, R) ≥ τ ,
P (−∞, L) ≥ τ /2 and P (R, ∞) ≥ τ /2. By Theorem 1, all intervals (L, R) corresponding
b ) lie in the set J (τ ), with probability tending to
to the splits in the redefined collection S(τ
one. We will only use the last two inequalities in the definition of J (τ ) to ensure that, with
probability tending to one, all (L, R) are contained in a fixed bounded interval, on which f
is bounded away from zero.
Main Body of the Proof
The following two lemmas give us appropriate control over the difference between small
bL,R and GL,R .
perturbations in G

35

Lemma 5. For every positive  and τ , there exists a Op (1) random sequence Mn , such that
bL,R (R) − G
bL,R (a)) − (GL,R (R) − GL,R (a))
(G

≤ (R − a) + (log n/n)Mn

bL,R (L) − G
bL,R (a)) − (GL,R (L) − GL,R (a))
(G

≤ (a − L) + (log n/n)Mn ,

bL,R (a) is well-defined.
for all (L, R) ∈ J (τ ) and a ∈ (L, R), such that G
The proof of Lemma 5 is given further below. The next result can be established using a
similar type of argument, however, it is a direct corollary of Lemma 4.1 in Kim and Pollard
(1990). Note that the log n factor, which appears in Lemma 5 but not in Lemma 6, is due
bL,R is approximated in Lemma 5, are
to the fact that points L and R, near which function G
allowed to vary, while s, the corresponding point in Lemma 6 is fixed.
Lemma 6. Given positive τ and , and a point s ∈ (L0 , R0 ), there exists a Op (1) sequence
Mn , such that
bL,R (s) − G
bL,R (a)) − (GL,R (s) − GL,R (a)) ≤ |s − a|2 + n−2/3 Mn ,
(G
bL,R (a) is
for all L,a and R, such that (L, s) ∈ J (τ ), (s, R) ∈ J (τ ), a ∈ (L, R) and G
well-defined.
We start by deriving the uniform rate of convergence for the small sample splits. As
in the proof of consistency, we focus on the case of the population procedure performing
a two-sided truncation, followed by a split. The rest of the cases can be handled using
analogous arguments. Let (L∗ , s∗ , R∗ ) be the first population split. The population procedure truncates (L0 , R0 ), along (L, RL ), down to (L∗ , R∗ ). As we showed in the consistency proof, the sample procedure reduces the cluster consisting of all the observations
down to, approximately, (L∗ , R∗ ), by repeatedly splitting off small clusters near the boundary. The split point for the big split is placed near s∗ . Given a positive δ, we define
A(δ) = {(L, R) : L ≤ L∗ + δ, |R − RL | ≤ δ}. In the consistency proof we showed that
the following property holds with probability tending to one for the starting cluster, (L, R),
bL,R ⊂ [L, L + δ] ∪ [R − δ, R].
of each small sample split: (L, R) ∈ A(δ) and arg max G
bL,R . Then, to establish bound (10) in the statement of
Let b
aL,R = arg max[L,L+δ]∪[R−δ,R] G
Theorem 2 it is sufficient to prove that for every positive τ there exists a positive δ and a
36

Op (log n/n) sequence Bn , such that inequality (R − b
aL,R ) ∧ (b
aL,R − L) ≤ Bn holds for all
(L, R) ∈ A(δ) ∩ J (τ ).
To characterize small perturbations in the population criterion function, GL,R , we will
use the derivations in the proof of Lemma 3 (we focus only on the regular part of the
proof, because of the new regularity condition, C3). Those derivations imply that there
exist positive c0 and δ, such that, for all (L, R) ∈ A(δ) ∩ J (τ ),
GL,R (R) − GL,R (a) ≥ c0 (R − a),

for a ∈ [R − δ, R], and

(25)

GL,R (L) − GL,R (a) ≥ c0 (a − L),

for a ∈ [L, L + δ].

(26)

Define
∆n (L, a, R) =

bL,R (R) − G
bL,R (a)) − (GL,R (R) − GL,R (a)) +
(G

bL,R (L) − G
bL,R (a)) − (GL,R (L) − GL,R (a)) .
(G
Combining Lemma 5, applied for  = c0 /2, with inequalities (25) and (26), we deduce
bL,R (R) − G
bL,R (b
0 ≥ G
aL,R )

≥ c0 (R − b
aL,R ) ∧ (b
aL,R − L) − ∆n (L, b
aL,R , R)
≥ (c0 /2)(R − b
aL,R ) ∧ (b
aL,R − L) − (log n/n)Mn .
Hence, (R − b
aL,R ) ∧ (b
aL,R − L) ≤ (c0 /2)−1 (log n/n)Mn for all (L, R) ∈ A(δ) ∩ J (τ ),
which is what we needed to prove. Thus, up to the first big split, the sizes of all the small
sample splits are uniformly Op (log n/n). In the consistency theorem we showed that, with
probability tending to one, the number of the big sample splits equals the (finite) number of
the population splits. Consequently, the behavior of the sample procedure after the first big
split can be handled by repeating the argument given above. This establishes bound (10) in
the statement of Theorem 2.

37

We now focus on deriving the rate of convergence in (9). Recall that (L∗ , s∗ , R∗ ) and
b∗ , sb∗ , R
b∗ ) denote the first population split and the first big sample split, respectively.
(L
Given a positive δ, define B(δ) = {(L, s, R) : |L − L∗ | ≤ δ, |s − s∗ | ≤ δ, |R − R∗ | ≤ δ}.
b∗ , sb∗ , R
b∗ ) ∈ B(δ) with probability tending to one. Derivations in the proof
Note that (L
of Lemma 3, together with the assumption G00L∗ ,R∗ (s∗ ) 6= 0, imply that there exist positive
constants δ, c1 and c2 , such that
GL,R (s∗ ) − GL,R (a) ≥ c1 |a − s∗ |2 − c2 |a − s∗ | (|L − L∗ | + |R − R∗ |) ,
b b∗ b∗ and let
for all (L, s, R) ∈ B(δ). Define sb = arg max[s∗ −δ,s∗ +δ] G
L ,R
∗
b b∗ b∗ (s∗ ) − G
b b∗ b∗ (b
s) .
Dn = G
b ∗ ,R
b∗ (s ) − GL
b ∗ ,R
b∗ (b
L ,R
L ,R s) − GL

We can handle Dn by applying Lemma 6, with sufficiently small  and τ . It follows that,
with probability tending to one,
b b∗ b∗ (b
b b∗ b∗ (s∗ ) − G
0 ≥ G
L , R s)
L ,R
≥ GLb∗ ,Rb∗ (s∗ ) − GLb∗ ,Rb∗ (b
s) − Dn


b∗ − L∗ | + |R
b∗ − R∗ | − n−2/3 Mn .
≥ c1 |b
s − s∗ |2 − c2 |b
s − s∗ | |L
Consequently,


b∗ − L∗ | + |R
b ∗ − R∗ | .
|b
s − s∗ | = Op n−1/3 + |L
b∗ , sb∗ , R
b∗ ) it is only left to show |L
b∗ −L∗ | = Op (n−1/3 )
Thus, to establish the rate in (9) for (L
b∗ − R∗ | = Op (n−1/3 ). In the remainder of the proof we establish the last two
and |R
b∗ , sb∗ , R
b∗ ) to the
stochastic bounds, and then extend the rate of convergence derived for (L
subsequent big sample splits.
We again use the derivations in the proof of Lemma 3, from which it follows that, given
an arbitrarily small positive τ , there exist positive constants δ, c1 and c2 , such that
GL,R (R) − GL,R (a) ≥ c1 (R − RL ) − c2 (a − L),

38

for all (L, R) ∈ A(δ) ∩ J (τ ) and a ∈ [L, L + δ]. Taking advantage of a maximal inequality for the empirical process indexed by a VC class of functions with a square integrable
envelope (e.g. Lemma 19.38 in van der Vaart 1998), Lemma 5, and the already established
rate of convergence for the small sample splits, we deduce that there exists a Op (n−1/2 )
sequence Cn , such that
bL,R (R) − max G
bL,R (a) ≥ c1 (R − RL ) − Cn ,
G
a∈[L,L+δ]

for all (L, R) ∈ A(δ) ∩ J (τ ). Consequently, when R − RL > c−1
1 Cn , the maximizer of
bL,R is at least a δ away from L. Because the split points of the small sample splits are
G
uniformly within Op (log n/n) of the cluster boundary, it follows that R − RL is bounded
above by a 2c−1
1 Cn for all the sample clusters, (L, R), produced up to the first big sample
split. A similar argument that focuses on L instead of R derives a lower bound for R − RL .
b∗ − L∗ | = Op (n−1/3 ).
Thus, |R − RL | = Op (n−1/2 ), and it is only left to establish |L
It follows from the derivations in the proof of Lemma 3 that there exist positive constants δ, c1 , c2 and c3 , such that
GL,R (s∗ ) − GL,R (a) ≥ c1 (L − L∗ ) − c2 (R − RL ) − c3 [(a − L) ∧ (R − a)],

(27)

provided L ∈ [L∗ , L∗ + δ], |R − R∗ | ≤ δ and |a − L| ∧ |a − R| ≤ δ. Handling the
stochastic term the same way we did in the paragraph above, we conclude that there exists
a Op (n−1/2 ) sequence Tn , such that inequality
bL,R (s∗ ) ≤
G

max
a∈[L,L+δ]∪[R−δ,R]

bL,R (a)
G

implies L − L∗ ≤ Tn . Note that the same Tn is chosen for all intervals (L, R) satisfying
L ∈ [L∗ , L∗ + δ] and |R − R∗ | ≤ δ. The companion lower bound on L − L∗ follows
from an analogous argument, which replaces the lower bound in display (27) with an upper
b∗ − L∗ | = Op (n−1/2 ) then follows directly after taking into
bound. The relationship |L
consideration the rate of convergence already established for the small sample splits.
This completes the derivation of the rate of convergence for the first big sample split.
For the subsequent big splits the same argument can be repeated with the following negligible modification. Due to the fact that endpoints of the clusters resulting from the first

39

big sample split are within Op (n−1/3 ) of their population counterparts, all the Op (n−1/2 )
expressions should be replaced with Op (n−1/3 ).
Proof of Lemma 5
We will only establish the first inequality. The second can be proved by using the same
argument with the appropriate adjustment of the notation. For points a that are bounded
away from R, the first inequality holds by the functional generalization of the GlivenkoCantelli theorem (e.g. Theorem 19.13 in van der Vaart 1998). Hence, from here on we only
focus on a ∈ [(L + R)/2, R]. We will use “.” to mean that inequality “≤” holds when
the right hand side is multiplied by a positive constant, which is chosen independently
from the involved parameters, such as n, L and R. Note that the statement of the lemma
restricts L, a and R to a bounded interval, on which the infumum of f is positive, while the
supremum is finite. In particular, we have R − a . P (a, R) and P (a, R) . R − a. Define
dl,r (x) = (r − x)(r − l)−1 1{l<x<r} and note that
|b
µa,R − µa,R | =

Pn da,R (R − a) P da,R (R − a)
−
Pn (a, R)
P (a, R)

≤ |Pn da,R − P da,R |

(R − a)
(R − µ
ba,R )
+ |Pn (a, R) − P (a, R)|
P (a, R)
P (a, R)

. |Pn da,R − P da,R | + |Pn (a, R) − P (a, R)| =: E1 + E2 .
As before, we let hl,r (x) = x1l<x<r . Observe that |Pn ha,R − P ha,R | . E2 , and stochastic
bound supl<r (|Pn hl,r − P hl,r |) + (|Pn (l, r) − P (l, r)|) = op (1) holds by the functional
generalization Glivenko-Cantelli theorem. It follows that
|(µL,R − µL,a ) − (b
µL,R − µ
bL,a )|

=

P hL,a P (a, R) + P ha,R P (L, a) Pn hL,a Pn (a, R) + Pn ha,R Pn (L, a)
−
P (L, R)P (L, a)
Pn (L, R)Pn (L, a)

. E2 + [(R − a) + E2 ] op (1),

(28)

40

where the op (1) term comes from the functional generalization of the Glivenko-Cantelli
theorem and is chosen uniformly for all L, R and a in consideration.
It is only left to show that the bound in the statement of the lemma holds for E1 and E2 .
We start with E2 . We need to show that there exists a Op (1) sequence of random variables
Mn , such that
|Pn (l, r) − P (l, r)| ≤ (r − l) + (log n/n)Mn ,

(29)

for all (l, r) contained within a given bounded set. Define Mn as the infimum of those
values for which the above inequality holds. Write Aj,n for the set of intervals (l, r) that
satisfy [2j−1 − 1] log n/n < P (l, r) ≤ 2j log n/n. Recall that we have restricted our
attention to a bounded interval, on which kf k∞ is finite. In the argument that follows the
bounded interval is not explicitly present, however we do use the fact that c1 := kf k−1
∞ is
positive. Observe that
P (Mn > t) ≤

∞
X

P ({∃(l, r) ∈ Aj,n s.t. |Pn (l, r) − P (l, r)| > (r − l) + t log n/n}

j=1

≤

∞
X

P ({sup |Pn (l, r) − P (l, r)| > (log n/n)[c1 (2j−1 − 1) + t]}
Aj,n

j=1

≤

∞
X
E supAj,n |Pn (l, r) − P (l, r)|
j=1

(log n/n)[c1 (2j−1 − 1) + t]

:=

∞
X

Sj .

j=1

We will bound the summands in the above expression using a maximal inequality from
the empirical process theory. Note that the class of indicator functions for all intervals
has polynomial bracketing numbers with respect to L2 (P ), and its envelope function is
identically equal to one. Consequently, an application of Lemma 19.36 in van der Vaart
(1998) gives us the following bound:
E

sup |Pn (l, r) − P (l, r)| . n−1/2 θ[1 +

p

− log(θ ∧ 1)] − log(θ ∧ 1)/n.

(30)

P (l,r)<θ2

We can bound each Sj by applying inequality (30) with θ2 = 2j log n/n:
Sj .

2j/2 log n/n
2j/2
=
.
(log n/n)[(2j−1 − 1) + t]
(2j−1 − 1) + t

It follows that the sum bounding P (Mn > t) can be made arbitrarily small, uniformly over
41

all n, by increasing t. Consequently, Mn = Op (1).
The proof of the bound for E1 is essentially identical to the one above. The class
of functions dl,r , for l and r in a fixed bounded interval, also has polynomial bracketing
numbers with respect to L2 (P ), and its envelope function is also identically equal to one.
Hence, the argument above works for E1 after each (l, r) is replaced with dl,r .

6.3

Proof of Propositions 1-4

Proposition 1
Let λ0 = 0. Suppose that the first merge happens at λ = λ1 , the second at λ = λ2 , and so
on. We will first show that, with probability one, values λk form an increasing sequence.
Consider two merges, C3 = C1 ∪ C2 and C = C3 ∪ C4 . For concreteness, we will focus on
the case where cluster C4 exists at the time of the first merge, and establish
X C2 − X C1
X C4 − X C3
≤
.
|C3 |
|C|

(31)

This will complete the proof because of the continuity of the underlying distribution. The
complementary case, where C4 is formed after the first merge can be analyzed analogously.
Suppose that the above inequality does not hold. Then, taking into account representation
|C3 |X C3 = |C1 |X C1 + |C2 |X C2 , we can derive
X C4 − X C2 = X C4 − X C3 − |C1 | · |C2 |−1 · (X C3 − X C1 )
= X C4 − X C3 − |C1 | · |C3 |−1 · (X C2 − X C1 )


|C|
|C1 |
|C2 | + |C4 |
−
=
(X C2 − X C1 ).
< (X C2 − X C1 )
|C3 | |C3 |
|C3 |
The resulting inequality contradicts the merge C3 = C1 ∪ C2 .
We will now verify that the KKT conditions for optimization problem (2) hold for the
solutions produced by Algorithm 1. The KKT conditions are satisfied if there exist βij with
|βij | ≤ 1 and βij = −βji , such that for every i:
α i − xi + λ

X

sign(αi − αj ) + λ

j6=i,αj 6=αi

X

βij = 0.

(32)

j6=i,αj =αi

Write C(i) for the current cluster containing xi . Taking into account equations (3), the
42

KKT conditions can be rewritten as follows,
X

X C(i) − xi + λ

βij = 0.

(33)

j6=i,j∈C(i)

We will argue by induction over the number of merges and establish the KKT conditions,
(33), for all i and for all λ = λk . If (33) holds for a particular tuning parameter value λk ,
it also holds for λ > λk , provided the cluster C(i) is not modified, if we shrink all the
corresponding βij by a factor of λk /λ.
For λ = λ0 , conditions (33) hold trivially, with each observation forming its own cluster. Suppose we are able to verify the KKT conditions up to the merge k − 1. Suppose that
the k-th merge, at λ = λk , is C = C1 ∪ C2 . By the discussion above, conditions (33) hold
at λ = λk for all i ∈
/ C, and there exist βij with |βij | ≤ 1 and βij = −βji , such that
X

X C1 − xi + λk

βij = 0

for all i ∈ C1

(34)

βij = 0

for all i ∈ C2 .

(35)

j6=i,j∈C1

X

X C2 − xi + λk

j6=i,j∈C2

We will set βij = −1 for each i ∈ C1 and j ∈ C2 , and we will keep the remaining βij
values intact. We need to show that
X C − xi + λ k

X

βij = 0 for all i ∈ C.

(36)

j6=i,j∈C

Consider an i ∈ C1 . Equations (34) imply λk

P

j6=i,j∈C1

βij = xi − X C1 . Recall that

λk = (X C2 − X C1 )/|C|. It follows that
X C − xi + λk

X

βij = X C − X C1 + (X C1 − X C2 )|C2 |/|C| = 0,

(37)

j6=i,j∈C

as required. The argument for i ∈ C2 is analogous, but uses equations (35) instead of (34).
Proposition 2
We will prove the result by induction over the number of merges. For each merge we will
establish the following claim: the splitting procedure applied to the last formed cluster

43

matches the merging procedure that formed that cluster. It follows that the sequence of
clusters formed by the splitting procedure matches the sequence of clusters formed by the
merging algorithm.
The claim for merge number one is trivial: the only possible split exactly matches the
first merge. Suppose that the claim has been established for the first k merges. Let C
be the cluster formed by the merge k + 1, which combines clusters C1 and C2 , with C1
being the left one. It is only left to show that the first splitting procedure applied to C
produces clusters C1 and C2 . Consider a possible alternative split: C = C3 ∪ C4 , with C3
being the left cluster that is different from C1 . To verify the claim, we need to establish
X C2 − X C1 > X C4 − X C3 . For concreteness, we will focus on the case C3 ⊂ C1 . The case
C1 ⊂ C3 can be handled analogously, taking advantage of the inclusion C4 ⊂ C2 .
Define C5 = C1 \ C3 . Using representations X C1 = X C3 |C3 |/|C1 | + X C5 |C5 |/|C1 | and
X C4 = X C5 |C5 |/|C4 | + X C2 |C2 |/|C4 |, we can rewrite the desired inequality as
X C2 − X C5
X C5 − X C3
>
.
|C4 |
|C1 |

(38)

Suppose that the cluster C1 was formed by the merge C11 ∪ C12 . The induction claim for
merges one through k implies that X C12 − X C11 maximizes the corresponding difference of
the averages over all partitions of C1 . By the monotonicity of the λ values in the merging
algorithm, we have (X C2 − X C1 )/|C| > (X C12 − X C11 )/|C1 |, which yields
X C2 − X C1
X C5 − X C3
>
.
|C|
|C1 |

(39)

Consequently, if we can establish that
X C2 − X C5
X C2 − X C1
>
,
|C4 |
|C|

(40)

then the required inequality (38) is satisfied. Using representation X C1 = X C3 |C3 |/|C1 | +
X C5 |C5 |/|C1 | we can rewrite (40) as (X C2 − X C1 )/|C| > (X C5 − X C3 )/|C1 |. The last
inequality is true by (39), which completes the proof.

44

Proposition 3
We will first establish that the population procedure produces no splits. More specifically,
we will show that for every s that is an extremum of GL,R (·) on the interval (L, R), we
have GL,R (s) < GL,R (L) ∨ GL,R (R). First, consider the case where s is less than or equal
to the mode of the density, f , restricted to (L, R). Differentiating function GL,R gives
G0L,R (s) =

f (s)PL,R
[µL,s + µs,R − µL,R − s].
PL,s Ps,R

Thus, if s is an extremum, then µs,R − µL,R = s − µL,s (see Figure 3).

Figure 3: Illustration of the proof of Proposition 3 on a truncated normal density.

Consequently,
GL,R (s) = µs,R − µL,s = 2(s − µL,s ) + µL,R − s
< s − L + µL,R − s = GL,R (L).
The last inequality follows from the fact that f is strictly increasing on (L, s), which implies
s − µL,s < (s − L)/2. Hence, GL,R (s) < GL,R (L). In the case where s is greater than the
mode of f on (L, R), an analogous argument establishes GL,R (s) < GL,R (R).
We will now prove that the population clustering procedure is uniquely defined. As we
have established in the first part of the proof, the population procedure has the following
options: (i) reduce the support to an empty set using a one-sided truncation; (ii) do the
same using two-sided truncation; and (iii) start with a one-sided truncation, followed by
the two-sided truncation, which reduces the support to an empty set. Because one-sided
truncation is always well defined, we can focus on the two-sided truncation exclusively.
45

Consider an interval (L1 , R1 ), such that arg max GL1 ,R1 = {L1 , R1 }. We need to show
that there exists, one and only one, continuous and decreasing function l 7→ Rl , defined
on [L1 , L2 ], for some L2 > L1 , such that R1 = RL1 and L2 = RL2 . Recalling the more
general existence and uniqueness argument in the proof of Proposition 4, we note that it is
enough to verify that
[f (L) ∨ f (R)](R − L) < PL,R ,

(41)

for all L < R with L ≥ L1 , R ≤ R1 and arg max GL,R = {L, R}. Suppose that inequality (41) is violated. Then, unimodality of f implies f (L) 6= f (R). Another consequence of
the violation is that µL,R < (L + R)/2 if f (L) > f (R) and µL,R > (L + R)/2, otherwise.
This contradicts the equality µL,R = (L + R)/2, implied by GL,R (L) = GL,R (R). We
conclude that [f (L) ∨ f (R)](R − L) < PL,R . This completes the proof.
Note that in the first part of the proof we establish the following general property,
which holds for any continuous f : if the population procedure splits (L∗ , R∗ ) at point s∗ ,
then µL∗ ,s∗ ≤ (L∗ + s∗ )/2 and µs∗ ,R∗ ≥ (s∗ + R∗ )/2. One consequence of this is that
a population split point must lie to the right of left-most mode of f and to the left of the
right-most mode.
This completes the proof of Proposition 3. We will now show that C2 holds for all
bimodal densities, provided the smoothness condition, C1, is satisfied.
Let m denote the interior local minimum of density f . Consider an interval (L, R) that
is a subset of the support of f . Suppose there exist points s1 < s2 , such that max GL,R =
GL,R (s1 ) = GL,R (s2 ). We will establish that for every s in (s1 , s2 ) equality G0L,R (s) = 0
implies G00L,R (s) < 0. Consequently, on the interval (s1 , s2 ) function GL,R is not flat and
has no local minima, which results in a contradiction.
Define gR (s) = f (s)(µs,R − s)/Ps,R . Note that
G00L,R (s)


= f (s)[2gR (s) − 1]

1
1
+
Ps,R PL,s



when G0L,R (s) = 0.

(42)

Because G0 (s1 ) = 0 and G00L,R (s1 ) ≤ 0, we have 2gR (s1 ) − 1 ≤ 0.
In the proof of Proposition 3 we showed that each interior maximizer of GL,R must lie in
between the two modes of f . Suppose that s1 is located to the left of m. It follows that f is
decreasing on (s1 , m). Define κ(s) = Ps,R (2gR (s) − 1) and recall that κ(s1 ) ≤ 0. Because
κ0 (s) = 2f 0 (s)(µs,R − s) + f (s)(2gR (s) − 1), we have κ0 (s1 ) < 0, and, consequently,

46

κ0 (s) < 0 and κ(s) < 0 for all s in (s1 , m]. Thus, it follows from formula (42) that, for s in
(s1 , m], equality G0L,R (s) = 0 implies G00L,R (s) < 0.
Define gL (s) = f (s)(s − µL,s )/PL,s , and note that gR (s) = gL (s) when G0L,R (s) = 0.
An argument analogous to the one in the previous paragraph, but with s2 and gL used
instead of s1 and gR , establishes that, for s in (m, s2 ), equality G0L,R (s) = 0 implies
G00L,R (s) < 0. This completes the proof.
Proposition 4
We will first focus on the case where the two-sided truncation starts at some bounded
interval (L1 , R1 ). Thus, RL1 = R1 . We will show that there exists a point L∗ > L1
and a small  > 0, such that the function RL , with the required properties, is uniquely
defined either for L ∈ (L1 − , L∗ ], where L∗ = RL∗ , or for L ∈ (L1 − , L∗ + ), where
the population procedure splits cluster (L∗ , RL∗ ). Recall that the required properties of
RL are that it is continuous and decreasing, arg max GL,RL = {L, RL } for L ∈ [L1 , L∗ ),
arg max GL∗ ,RL∗ ⊇ {L∗ , RL∗ }, and GL,RL (L) = GL,RL (RL ) for L ∈ (L1 − , L∗ + ).
Write F (L, R) for GL,R (L)/2−GL,R (R)/2, and note that F (L, R) = µL,R −(R+L)/2.
Observe that arg max GL,R = {L, R} implies F (L, R) = 0. To simplify the notation we
write F1 (L, R) for the derivative ∂F (l, r)/∂l, evaluated at (L, R), and write F2 (L, R) for
the corresponding value of ∂F (l, r)/∂r. Differentiation of F gives us
F1 (L, R) =

f (L)(µL,R − L) 1
−
PL,R
2

and F2 (L, R) =

f (R)(R − µL,R ) 1
− .
PL,R
2

(43)

According to display (16) in the proof of Lemma 3, we have F1 (L, R) = G0L,R (L) and
F2 (L, R) = −G0L,R (R). Because arg max GL1 ,R1 = {L1 , R1 }, we then also have inequalities F1 (L1 , R1 ) ≤ 0 and F2 (L1 , R1 ) ≤ 0. Suppose, first, that both of the above
inequalities are strict. Then, both F1 (L1 , R1 ) and F2 (L1 , R1 ) are negative. By the classical implicit function theorem, function L 7→ RL , for which F (L, RL ) = 0, and hence
GL,RL (L) = GL,RL (RL ), is uniquely defined on (L1 − , L1 + ), for some small positive
. By the same theorem RL is also continuous and decreasing.
As L increases from the value L1 , we define RL as the largest R below R1 , for which
F (L, R) = 0. Note that the smallest RL can be is L. We continuously decrease L until
F1 (L, RL ) = 0, or F2 (L, RL ) = 0, or L = RL , and let L2 be the value of L corresponding

47

to the stopping point. Note that, by the implicit function theorem, function RL is continuous
and decreasing on (L1 − , L2 ), as well as uniquely defined. Because RL is continuous,
decreasing and bounded below by L2 , we can define it at L2 by continuity, also preserving
the monotonicity and the lower bound. We write R2 for RL2 .
If L2 = R2 or L∗ < L2 , then the proof of existence and uniqueness of the two-sided
truncation is complete. Otherwise, we need to handle the case where F1 (L2 , R2 ) = 0 or
F1 (L2 , R2 ) = 0. Note that GL2 ,R2 (L2 ) = GL2 ,R2 (R2 ) = max GL2 ,R2 . We will now focus
on the case F1 (L2 , R2 ) = F2 (L2 , R2 ) = 0. The remaining two cases follow from analogous
arguments.
The second order Taylor series approximation to F (L, R) at (L2 , R2 ) is
F (L, R)−F (L2 , R2 ) =



(R2 − L2 )  0
f (L2 )(∆L)2 + f 0 (R2 )(∆R)2 +o |∆L|2 + |∆R|2 ,
2PL2 ,R2

where ∆R = R −R2 and ∆L = L−L2 . As we point out in the proof of Lemma 3, equality
G0L2 ,R2 (R2 ) = F2 (L2 , R2 ) = 0 implies f 0 (R2 ) ≥ 0. Similarly, we have f 0 (L2 ) ≤ 0.
Suppose, first, that both inequalities are strict; the remaining case can be handled using
analogous arguments that use higher order Taylor expansions. Then, given an arbitrarily
small neighborhood N of (L2 , R2 ) we can find (L, R) ∈ N , such that L > L2 and R > R2
and F (L, R) = 0. Also note that F1 (L, R) and F2 (L, R) are both negative. Consequently,
by the implicit function theorem, the map l 7→ Rl is uniquely defined in a neighborhood
of L, such that R = RL , and the function Rl is continuous and decreasing. Thus, the map
is uniquely defined on (L2 , L2 + ), for some positive , it is continuous and decreasing,
and Rl ↑ R2 as l ↓ L2 . Thus, the continuity and monotonicity of the map is preserved
when we glue together the functions defined on (L1 , L2 ] and (L2 , L2 + ). The above
Taylor approximation implies that when  is sufficiently small, derivatives F1 (L, R) and
F2 (L, R) are both negative for L ∈ (L2 , L2 + ). Thus, we are back in the regular setting
of the previous paragraph. Applying the above arguments sequentially establishes that the
function RL is uniquely defined, in a continuous and decreasing fashion, either for L ∈
(L1 − , L∗ ], where L∗ = RL∗ , or for L ∈ (L1 − , L∗ + ), where the population procedure
splits cluster (L∗ , RL∗ ). The remaining case (L1 , R1 ) = (−∞, ∞) can be handled using
an additional, but similar, argument, where L is decreased, rather than increased from the
value L1 .
Thus, we have shown that the two-sided truncation of the population procedure is
48

uniquely defined. Note that the one sided truncation exists and is unique directly by definition, while the uniqueness of the split is ensured by the regularity condition C2. It is only
left to show that the population procedure will implement finitely many steps, for which
we just need to check that the procedure will make finitely many splits.
In the proof of Proposition 3 we establish the following general property: a population
split point must lie to the right of left-most mode of the truncated version of f , and to the
left of the corresponding right-most mode, where the truncation is done to the cluster that
is being split. We will assume that the procedure makes infinitely many splits and conduct
an argument to reach a contradiction. The assumption implies that there exists an interval
(Li , Ri ), such that for every positive  there exists an infinitely large nested sequence of
distinct population clusters, of the form (L, R), sandwiched between (Li − , Ri + ) and
[Li , Ri ], such that either the left or right endpoint of each cluster is a population split point
for the previous cluster in the sequence. Without loss of generality, we focus on the case
where each right endpoint of the considered cluster sequence is a split point for the previous
cluster. By continuity of GL,R (s) we must then have GLi ,Ri (Ri ) ≥ GLi ,Ri (Li ). Below we
examine the case where the last inequality is strict; the case of equality can be handled
using an analogous argument. It follows that when  is sufficiently small, the population
procedure cannot place a split point in (L, Li ] or do a left sided truncation of the cluster
(L, R). Consequently, we have L = Li for each cluster in the considered sequence.
Recall that G0L,R (R) = 1/2−f (R)(R−µL,R )/PL,R . If G0 (Li , Ri )(Ri ) is negative, then,
for a sufficiently small , G0L,R (R) is negative and bounded away from zero. Consequently,
when R gets sufficiently close to Ri , the corresponding split point must lie below Ri , which
is a contradiction. Similarly, we reach a contradiction when G0 (Li , Ri )(Ri ) is positive. For
the rest of the proof we focus on the remaining case G0 (Li , Ri )(Ri ) = 0.
The general property given in the beginning of the paragraph before the previous one
implies that f 0 (R) > 0 for R ∈ (Ri , Ri + ) and f 0 (Ri ) ≥ 0. We first consider the case
where the last inequality is strict. Let t = s − Ri and r = R − Ri . After deriving the second
order Taylor series expansion for the function (s, R) 7→ GLi ,R (s) at (Ri , Ri ), we establish
approximation
GLi ,Ri +r (Ri + t) − GLi ,Ri +r (Ri ) =

−f 0 (Ri ) 2
(t + tr) + o(|t|2 + |r|2 ),
6f (Ri )

which holds for (t, r) near zero. Because we now focus on the case f 0 (Ri ) > 0, we must
49

have G(Li , Ri +r)(Ri −r/2) > maxs∈[Ri ,Ri +r] G(Li , Ri +r)(s), for every sufficiently small
positive r. Consequently, when  is sufficiently small, the population procedure applied to
an interval (Li , R) satisfying R ∈ (Ri , Ri + ) must place the split point below Ri . This is
again a contradiction of our assumption.
Now consider the remaining case f 0 (Ri ) = 0. Let t = s − Ri and r = R − Ri .
After deriving the third order Taylor series expansion for the function (s, R) 7→ GLi ,R (s)
at (Ri , Ri ), we establish approximation
GLi ,Ri +r (Ri + t) − GLi ,Ri +r (Ri ) =

−f 00 (Ri ) 3
(3t − t2 r + 3tr2 ) + o(|t|3 + |r|3 ),
6f (Ri )

for (t, r) near zero. We will focus on the case f 00 (Ri ) 6= 0, as the case f 00 (Ri ) = 0 can
be handled using analogous arguments that use higher order Taylor approximation. Because Ri is a local minimum of f , we must have f 00 (Ri ) > 0. Analysis of the above Taylor
expansion reveals that G(Li , Ri + r)(Ri − r) > maxs∈[Ri ,Ri +r] G(Li , Ri + r)(s), for every sufficiently small positive r. Consequently, when  is sufficiently small, the population
procedure applied to an interval (Li , R) satisfying R ∈ (Ri , Ri + ) must place the split
point below Ri . This contradiction completes the proof.

6.4

Analysis of the Population Procedure for Gaussian Mixtures

Consider the case where f is a mixture of two Gaussian densities on the real line. Note that
regularity condition C1 is satisfied. Proposition 3 imply that condition C2 is satisfied as
well. Thus, the population procedure is well defined, and, according to the results in Section 3, the sample clustering procedure is consistent. In Table 4 we document the behaviour
of the population clustering procedure for various mixtures of two Gaussian distributions
on the real line. For 7 different levels of separation between the two normal means we
consider 9 different mixing proportions, from the symmetric case of 50 : 50 mixing to
the highly non-symmetric 10 : 90 mixing. The behavior of the population splitting procedure in other cases can be interpolated from the table using continuity arguments. We
present, where applicable, the location of the split point, s∗ (“NO” denotes that no splits
was detected), as well as the endpoints, L∗ and R∗ , of the corresponding cluster. The local
minimum of the density, m, and the split point minimizing the expected misclassification
error, sMC , are also provided. We also report, under Excess MCE, how much the mis50

classification error of the population clustering procedure exceeds that of the sMC based
oracle rule. Even though the population procedure focusses on cluster separation, rather
than classification, its Excess MCE is well controlled when the number of splits is detected
correctly. In all cases where s∗ is not reported, we checked that the two-sided truncation
takes the left endpoint of the cluster all the way to m, and thus reduces the density to a unimodal one. Proposition 3 then implies that the population procedure does not produce any
splits. Similarly, in the cases where s∗ is provided, we verified that the bimodal sub-cluster
is truncated down to a unimodal one, which means that the procedure does not produce a
second split.
6.4.1

Additional Technical Details

Note that the population procedure for a distribution supported by the real line starts with
a two-sided truncation of the support. The truncation proceeds along the interval (L, RL ),
as L is increased, and stops when either the support is reduced all the way to an empty
set or the maximum of GL,RL is achieved at an interior point. In the latter case, define L∗
as the smallest L for which there exists an s in (L, RL ), such that µL,s = (L + s)/2 and
µs,RL = (s + RL )/2. Write s∗ for the corresponding point s, and let R∗ stand for RL∗ . The
population procedure truncates the support to the interval (L∗ , R∗ ), which is then split at s∗ .
Note that L∗ must be smaller than m, the local minimum of the density, by Proposition 3.
Because analytical solutions are not available, we find L∗ , R∗ and s∗ numerically. For
each L on a dense grid we locate RL using equation µL,RL = (L + RL )/2. The symmetry
of Gaussian distributions implies that in the search for L∗ we only need to consider the
values of L that satisfy inequality L ≥ 2(µ1 − µ2 ) + RL . If the above inequality holds,
then we focus on s, the local maximum of GL,RL , and compute δ1 = µL,s − (L + s)/2
and δ2 = µs,RL − (s + RL )/2. The proof of Proposition 4 shows that at L = L∗ we have
δ1 = δ2 = 0, and, as L crosses L∗ , inequalities δ1 < 0 and δ2 > 0 are satisfied for the first
time. Thus, as L is increased, L∗ can be taken as the first point where the above inequalities
hold. To check for the second split we apply the same approach on the bimodal sub-cluster.
In all the reported cases, we verify that the two-sided truncation takes an endpoint of the
sub-cluster all the way to m, and thus reduces the density to a unimodal one. Proposition 3
then implies that the population procedure does not produce a second split.
We now provide closed form expressions for µL,R , GL,R and G0L,R . Figure 4 displays
51

CASE

|µ2 − µ1 | = 9

|µ2 − µ1 | = 8

|µ2 − µ1 | = 7

|µ2 − µ1 | = 6

|µ2 − µ1 | = 5

|µ2 − µ1 | = 4

|µ2 − µ1 | = 3

p1
0.50
0.45
0.40
0.35
0.30
0.25
0.20
0.15
0.10
0.50
0.45
0.40
0.35
0.30
0.25
0.20
0.15
0.10
0.50
0.45
0.40
0.35
0.30
0.25
0.20
0.15
0.10
0.50
0.45
0.40
0.35
0.30
0.25
0.20
0.15
0.10
0.50
0.45
0.40
0.35
0.30
0.25
0.20
0.15
0.10
0.50
0.45
0.40
0.35
0.30
0.25
0.20
0.15
0.10
0.50
0.45
0.40
0.35
0.30
0.25
0.20
0.15
0.10

p2
0.50
0.55
0.60
0.65
0.70
0.75
0.80
0.85
0.90
0.50
0.55
0.60
0.65
0.70
0.75
0.80
0.85
0.90
0.50
0.55
0.60
0.65
0.70
0.75
0.80
0.85
0.90
0.50
0.55
0.60
0.65
0.70
0.75
0.80
0.85
0.90
0.50
0.55
0.60
0.65
0.70
0.75
0.80
0.85
0.90
0.50
0.55
0.60
0.65
0.70
0.75
0.80
0.85
0.90
0.50
0.55
0.60
0.65
0.70
0.75
0.80
0.85
0.90

µ1
-4.50
-4.50
-4.50
-4.50
-4.50
-4.50
-4.50
-4.50
-4.50
-4.00
-4.00
-4.00
-4.00
-4.00
-4.00
-4.00
-4.00
-4.00
-3.50
-3.50
-3.50
-3.50
-3.50
-3.50
-3.50
-3.50
-3.50
-3.00
-3.00
-3.00
-3.00
-3.00
-3.00
-3.00
-3.00
-3.00
-2.50
-2.50
-2.50
-2.50
-2.50
-2.50
-2.50
-2.50
-2.50
-2.00
-2.00
-2.00
-2.00
-2.00
-2.00
-2.00
-2.00
-2.00
-1.50
-1.50
-1.50
-1.50
-1.50
-1.50
-1.50
-1.50
-1.50

µ2
4.50
4.50
4.50
4.50
4.50
4.50
4.50
4.50
4.50
4.00
4.00
4.00
4.00
4.00
4.00
4.00
4.00
4.00
3.50
3.50
3.50
3.50
3.50
3.50
3.50
3.50
3.50
3.00
3.00
3.00
3.00
3.00
3.00
3.00
3.00
3.00
2.50
2.50
2.50
2.50
2.50
2.50
2.50
2.50
2.50
2.00
2.00
2.00
2.00
2.00
2.00
2.00
2.00
2.00
1.50
1.50
1.50
1.50
1.50
1.50
1.50
1.50
1.50

m
0.00
-0.02
-0.05
-0.07
-0.10
-0.13
-0.16
-0.20
-0.26
0.00
-0.03
-0.05
-0.08
-0.11
-0.15
-0.18
-0.23
-0.29
0.00
-0.03
-0.06
-0.10
-0.13
-0.17
-0.22
-0.27
-0.34
0.00
-0.04
-0.08
-0.12
-0.16
-0.21
-0.26
-0.33
-0.41
0.00
-0.05
-0.10
-0.15
-0.20
-0.26
-0.33
-0.41
-0.53
0.00
-0.07
-0.14
-0.21
-0.28
-0.37
-0.47
-0.58
-0.74
0.00
-0.12
-0.24
-0.38
-0.53
-0.71
-1.50
-1.50
-1.50

s?
0.00
-0.45
-0.90
-1.36
-1.82
-2.31
-2.90
-3.82
NO
0.00
-0.40
-0.80
-1.22
-1.64
-2.12
-2.72
NO
NO
0.00
-0.35
-0.71
-1.09
-1.49
-1.97
-2.66
NO
NO
0.00
-0.32
-0.64
-0.99
-1.39
-1.91
NO
NO
NO
0.00
-0.30
-0.61
-0.96
-1.41
NO
NO
NO
NO
0.00
-0.32
-0.67
-1.12
NO
NO
NO
NO
NO
0.00
-0.50
NO
NO
NO
NO
NO
NO
NO

L?
-8.99
-8.53
-8.08
-7.62
-7.17
-6.67
-6.09
-5.09

R?
8.99
9.43
9.88
10.33
10.79
11.24
11.70
12.16

2nd split
NO
NO
NO
NO
NO
NO
NO
NO

-7.99 7.99
-7.58 8.38
-7.17 8.78
-6.77 9.19
-6.34 9.59
-5.86 9.99
-5.25 10.40

NO
NO
NO
NO
NO
NO
NO

-6.98
-6.63
-6.27
-5.90
-5.49
-5.01
-4.32

6.98
7.34
7.69
8.04
8.39
8.75
9.12

NO
NO
NO
NO
NO
NO
NO

-5.99
-5.66
-5.34
-4.99
-4.59
-4.07

5.99
6.28
6.59
6.89
7.20
7.52

NO
NO
NO
NO
NO
NO

-4.97
-4.68
-4.37
-4.01
-3.56

4.97
5.23
5.49
5.75
6.02

NO
NO
NO
NO
NO

-3.89
-3.62
-3.28
-2.85

3.89
4.15
4.38
4.62

NO
NO
NO
NO

-2.68
-2.29

2.68
2.97

NO
NO

sMC
0.00
-0.02
-0.04
-0.07
-0.09
-0.12
-0.15
-0.19
-0.24
0.00
-0.03
-0.05
-0.08
-0.11
-0.14
-0.17
-0.22
-0.28
0.00
-0.03
-0.06
-0.09
-0.12
-0.16
-0.20
-0.25
-0.31
0.00
-0.03
-0.07
-0.10
-0.14
-0.18
-0.23
-0.29
-0.37
0.00
-0.04
-0.08
-0.12
-0.17
-0.22
-0.28
-0.35
-0.44
0.00
-0.05
-0.10
-0.15
-0.21
-0.28
-0.35
-0.43
-0.55
0.00
-0.07
-0.14
-0.21
-0.28
-0.37
-0.46
-0.58
-0.73

Excess MCE
0.000
0.000
0.000
0.000
0.001
0.004
0.011
0.037
0.100
0.000
0.000
0.000
0.001
0.003
0.008
0.020
0.150
0.100
0.000
0.000
0.001
0.003
0.007
0.016
0.040
0.150
0.100
0.000
0.001
0.004
0.008
0.016
0.034
0.200
0.150
0.100
0.000
0.005
0.011
0.021
0.041
0.250
0.200
0.149
0.100
0.000
0.015
0.034
0.065
0.298
0.248
0.198
0.148
0.097
0.000
0.057
0.396
0.344
0.293
0.241
0.190
0.139
0.090

Table 4: Finding the population splits for 2-normal mixtures: p1 N (µ1 , 1) + p2 N (µ2 , 1).
52

GL,R (a), f and G0L,R (normalized), corresponding to the 0.35 N (−4, 1) + 0.65 N (4, 1)
distribution, for three values of L: (a) L < L∗ , (b) L = L∗ , (c) L > L∗ . Locations of
µL,s , (L + s)/2, µs,RL and (s + RL )/2 are also provided.
If the population density f is a mixture of k normals: f (x) :=

Pk

i=1

wi φ(x − µi ) then

using the following inequality,
Z
x φ(x − µ) dx = µ Φ(x − µ) − φ(x − µ),
we get closed from expressions of GL,R which can be subsequently optimized:
GL,R (a) =

−1
Pa,R

Z

R

xf (x) dx −

−1
PL,a

Z

a

Z

R

xf (x) dx =
a

Z

a

xf (x) dx =
L

k
X
i=1
k
X

a

xf (x), dx

where

L

wi µi {Φ(R − µi ) − Φ(a − µi )} − f (R) + f (a),
wi µi {Φ(a − µi ) − Φ(L − µi )} − f (a) + f (L), and

i=1

Pa,R =

k
X



wi Φ(R − µi ) − Φ(a − µi ) ,

PL,a =

i=1

k
X


wi Φ(a − µi ) − Φ(L − µi ) .

i=1

The following alternative expression for GL,R (a) in terms of the conditional means and
−1
−1
density fL,R (x) = PL,R
f (x) I{L, R}, where P̃a1 ,a2 = PL,R
Pa1 ,a2 , is also useful:

−1
P̃L,a

Z

R
−1
P̃a,R

Z

a

x fL,R (x) dx −
x fL,R (x) dx
a
L


Z a
−1
= {P̃L,a P̃a,R }
P̃L,a µL,R −
x fL,R (x) dx
∞

Z
−2
2
= fL,R (a){P̃L,a P̃a,R } µL,R P̃a,L + (1 − 2 P̃a,L )

GL,R (a) =

a


x fL,R (x) dx − P̃a,L P̃a,R .

−∞

Differentiating the above with respect to a, we arrive at: G0L,R (a) = κL,R (a) × HL,R (a),

53

9
7
5

GL, R(a)

3

0

5

10

−10

−5

0

5

10

−10

−5

0

5

10

−10

−5

0

5

10

−10

−5

0

5

10

−10

−5

0

5

10

−10

−5

0

5

10

−10

−5

0

5

10

−10

−5

0

5

10

0.6

−5

−0.2
−0.4
−0.6

dGL, R da

0.0

0.2

0.0

density

−10

Figure 4: Across rows we have the plots of GL,R (a), the truncated density and the normalized G0L,R (a), as L varies across columns over the 3 cases: (i) L < L? (ii) L = L? (iii)
L > L? . The population density used here is 0.35 N (−4, 1) + 0.65N (4, 1). The dotted
blue line and the stars denote the position of the zero of G0 . In the last row, the conditional means, µLs , µLR and µsR , are denoted by circles and the corresponding mid-points,
(L + s)/2, (L + RL )/2 and (s + RL )/2, by squares. The signs differences δ1 and δ2 , defined
in Section 3, vary as follows (i) (δ1 > 0, δ2 < 0) (ii) (δ1 = 0, δ2 = 0) (iii) (δ1 < 0, δ2 > 0).

54

where
−2 −2
−1
κL,R (a) = f (a) PL,a
Pa,R PL,R
, and

HL,R (a) =

2
µL,R PL,a

Z

a

+ {Pa,R − PL,a }

xf (x) dx − a Pa,L Pa,R .
L

Note that κL,R (a) > 0 for all a ∈ (L, R). Hence, to track the the extremas of GL,R , it is
enough to search for the zeros of HL,R (a). We call HL,R the normalized G0L,R . Figure 4
shows the plots of GL,R (a), the truncated density and HL,R (a), when f = 0.35 N (−4, 1) +
0.65N (4, 1). The plot of GL,R appears flat in the neighborhood of the split. But, the plot of
HL,R clearly shows only one zero-crossing and demonstrates uniqueness of the maximum
of GL,R in the case of interest.

6.5

Further Details on the Simulation Study & Real Data Analysis

For the numerical experiments in Section 4 of the main paper, we implemented the BMT
with threshold α uniformly set at 10%, and with the adjustment of at most 50% truncation
for the first split. The details are provided in Algorithm 4.
6.5.1

Modality Assessment

We compare the performance of the BMT with the following two popular modality assessment procedures:
Silverman Test is based on a kernel density estimate. It uses the idea that if the population density is non-unimodal, a large value of the bandwidth will be required to
smooth the data to a unimodal density estimate (Silverman, 1981) . The test uses the
minimum bandwidth that produces a unimodal kernel estimator. Large values of the
minimum bandwidth based test-statistic provide evidence to support the alternative
hypothesis of multi-modality. To conduct the Silverman test, we use the R-package
referenced in Vollmer et al. (2013). It is based on Gaussian kernels and incorporates
Hall and York (2001) adjustment for calculating the p-value.
The Dip Test proposed by Hartigan and Hartigan (1985) is a histogram based method,
which does not require estimating the density. The Dip-statistic is the minimum
Kolmogorov-Smirnov distance between the empirical distribution and the class of
55

INITIALIZE:
K = number of clusters = n.
Sort data in ascending order and store them as: x = {x1 , . . . , xn }.
Assign cluster mean {a1 , a2 , . . . , an } to them: ai = xi for i = 1, . . . , n.
Cluster size: si = 1, i = 1, . . . , n.
Cluster Membership Indices of x: I(x) = {1, . . . , n}.
WHILE K > 1:
Find the consecutive adjacent centroid distance standardized by cluster sizes:
d(j, j + 1) ← (aj+1 − aj )/(sj + sj+1 )
Find the clusters with minimum merging distance:
j ? ← arg min1≤j≤K−1 d(j, j + 1)
Check if it is a Big Merge: min{sj ? , sj ? +1 } > dnαe
IF Big Merge: Find and Store Mass after merge = (sj ? + sj ? +1 )/n and
New Split = max{x[I(x) is j ? } sj ? + min{x[I(x) is (j ? + 1)]} sj ? +1




(sj ? + sj ? +1 );

Merge the j ? and (j ? + 1) clusters and update the centroid and size of the new cluster:
aj ? ← (sj ? aj ? + sj ? +1 aj ? +1 )/(sj ? + sj ? +1 )
sj ? ← (sj ? + sj ? +1 )
Reduce the number of clusters in path K ← (K − 1)
Change cluster indices & cluster member indices of data according to the above reduction:
FOR k in (j ? + 1) : K, sk ← sk+1 ; ak ← ak+1
FOR ALL I(x) > j ? : reduce index by 1 ,i.e., I(x) = I(x) − 1
ADJUSTMENT: IF Mass after merge in the TOP SPLIT < 50%, Stored.Splits = NULL;
OUTPUT Stored Splits.
Algorithm 4: α-thresholded BMT algorithm with truncation adjustment

56

unimodal distributions. Larger values of the Dip-statistic signify departure from
the null hypothesis of unimodality. P-values are calculated using the R-package of
Maechler (2013). The p-value of this test is quite conservative.
6.5.2

Estimating the number of clusters

We compare the performance of the BMT with eight statistical methods that are popularly
used for estimating the number of clusters in a dataset. A comparison study of 30 different
approaches in Milligan and Cooper (1985) reports the approach in Caliński and Harabasz
(1974) as being one of the best performing global methods. It prescribes maximizing the
following index over k:
CH(k) =

B(k)/(k − 1)
,
W (k)/(n − k)

where B(k) and W (k) are respectively the between and the within clusters sum of squares
for k clusters. Another popular approach, due to Krzanowski and Lai (1988), is based on
the changes in the within clusters sum of squares as new clusters are formed, and seeks to
maximize the following ratio over k:
KL(k) =

DIFF(k)
where DIFF(k) = (k − 1)2/p Wk−1 − k 2/p Wk .
DIFF(k + 1)

Both these approaches are not defined for k = 1 and can not be used for testing population
unimodality. Hartigan (1975) proposed using the smallest k for which the following ratio
of the within cluster sum of squares is greater than 10:

H(k) =


W (k)
−1
(n − k − 1).
W (k + 1)

It can be used for testing presence of only one cluster. Theoretical thresholds based on the
F distribution can also be used. Gordon (1996) further sub-divides these approaches into
local and global methods. Local methods consider individual pairs of clusters and check
whether they should be merged. On the other hand, global methods incorporate the entire
data in evaluating measures that are subsequently optimized as a function of the number of
clusters. Note that the BMT is a local method. In addition to the above three methods, we
compare BMT with the following five methods.

57

Scenarios 2,3 & 4

Scenario 5

−6

−4

−2

0

2

4

6

16, 4)
0.0

0.00

0.00

0.05

0.05

0.5

0.10

0.10

0.15

1.0

0.20

0.15

0.25

1.5

0.20

Scenario 1

−10

−5

0

5

10

0.0

Scenario 6 : dimensions 2 & 3

0.2

0.4

0.6

0.8

1.0

Scenario 6 : dimensions 4 & 5

−6

−4

−2

0

2

4

6

0.5
0.0

0.0

0.00

0.1

0.05

0.2

0.10

1.0

0.3

0.15

1.5

0.4

0.20

Scenario 6 : dimension 1

−6

−4

−2

0

2

4

6

−6

−4

−2

0

2

4

6

Figure 5: Plot of the univariate densities used in the different numerical experiments of
Table 2 of the main paper.
Given k clusters, for each data-point xi the silhouette statistic of Kaufman and Rousseeuw
(2009) uses: a(i) - the average distance of xi to other points in its cluster, and b(i) the average distance of xi from points in its nearest neighboring cluster. It is given by
shk (i) = (b(i) − a(i))/max{a(i), b(i)}. Large values of shk (i) signify good clustering.
A popular estimate of the optimal number of clusters is based on maximizing the average
silhouette statistic,
−1

KR(k) = n

n
X

shk (i)

over k ≥ 2.

i=1

The Gap statistic of Tibshirani et al. (2001) uses the ‘elbow phenomenon’ (Thorndike,
1953) by estimating the number of clusters at the transition point, where the decline in the
within cluster dispersion first slackens. The goodness of clustering for k clusters is defined

58

as:
Gapn (k) = E?n {log(W̃ (k))} − log(W̃ (k)),
where W̃ (k) is the size-normalized intra-cluster sums of squares. The expectation is over
reference datasets and can be estimated by the mean of log W̃ ∗ (k) over B i.i.d. datasets
that are generated by sampling uniformly from the original dataset’s range. The standard
deviation, std(k) of log W̃ ∗ (k), is also recorded, and an estimate of the optimal number of
clusters in the datasets is given by the smallest k for which the following holds:
Gap(k) ≥ Gap(k + 1) − (1 + B −1 )−1/2 std(k + 1).
Another method that uses the ‘elbow phenomenon’ is the Jump statistic of Sugar and
James (2003). It again is a global method and involves computation of the mean square
error (mse) associated with k clusters for different choices of k. The ‘jumps’ in the p/2
moments of the mean square errors are subsequently calculated: Jk = {mse(k)}−p/2 −
{mse(k − 1)}−p/2 for k = 1, . . . , n. The estimate of the number of clusters is the value of
k that maximizes the jump Jk .
The prediction strength criterion of Tibshirani and Walther (2005) is computed by repeatedly dividing the data-set into two halves: train and test. For each value of k ≥ 1,
we cluster the test and training data into k groups and measure how well the training set
cluster centroids predict co-memberships in the test set. For each of the k test clusters, the
proportion of observation pairs that are also assigned to the same cluster by the training
set centroids is computed. The prediction strength is the minimum of this value over the k
test clusters. We have used the R-package of Hennig (2014) for computing the prediction
strength criterion. 50 divisions of the data were used and the estimate of the number of
clusters is the maximum the value of k which has prediction strength above 0.8.
The number of clusters estimate by Fang and Wang (2012) involves drawing 2 bootstrap samples from the data several times. For each of these bootstrap sample pairs and
for each value of k, a global clustering instability measure (Wang, 2010b), which reflects
the distance between clusterings in paired samples, is computed. The number of clusters
estimate minimizes the clustering instability aggregated over the bootstrap samples.
In Tables 5, 6 and 7 we report the performance of the aforementioned seven number of

59

clusters estimation approaches (all except for the BMT and Jump), across the six simulation
experiments, for the following clustering algorithms:
• Ward’s algorithm (Ward Jr, 1963), which was implemented through the NbCLust
package.
• Centroid-based clustering algorithm (Kaufman and Rousseeuw, 2009).
• PAM (Partitioning Around Medoids,) which is the most common realization of kmedoid clustering algorithm (Kaufman and Rousseeuw, 1990). It uses a greedy
search, which may not find the optimum solution but is faster than exhaustive search.
• CLARA (Clustering LARge Applications) algorithm of Kaufman and Rousseeuw
(1986), which extends the k-medoids approach for large sample sizes. It first conducts a down-sampling and then clusters the down-sampled observations. Thereafter,
it assigns all objects in the dataset to the clusters that were formed.
• GMM-merge, which conducts clustering by merging Gaussian mixture components
from an initial mclust clustering (Hennig, 2010). The mixture components are merged
in a hierarchical fashion. The merging criterion is computed for all pairs of current
clusters and the two clusters with the extremal criterion value are merged. Then,
criterion values are recomputed for the merged cluster and the merging process is
continued up to the length dictated by number of clusters selection methods.
• GMM-BIC, which is based on Fraley and Raftery (2002). It fits a parameterized
Gaussian mixture models by EM algorithm which is initialized by model-based hierarchical clustering.
Not all the above clustering algorithms were compatible with all of the number of clusters methods that we consider. Some of the combinations returned frequent computational
errors due to non-convergence and related issues. We restricted ourselves to those pairs
which were implementable in all our simulation set-ups. Also, as most of these clustering
algorithms were very computationally expensive, we could not conduct 100 repetitions of
the simulation scenarios as we had done for Table 2 in the main paper with the much faster
k-means algorithm. We limited ourselves to 20 independent trails here, and most of the
performance patterns were clearly visible by then. Figure 5, provided in this subsection,
60

True Population Density

0.3 N (−4, 1) + 0.7 N (4, 1)

0.3 N (−3, 1) + 0.35 N (0, 1) + 0.35 N (3, 1)

0.3 t1 (−3) + 0.35 t1 (0) + 0.35 t1 (3)

0.3 dexp(−3) + 0.35 dexp(0) + 0.35 dexp(3)



Beta(8, 2) + Beta(5, 5) + Beta(2, 8) 3

{0.5 N (−2, 1) + 0.5 N (2, 1)}
⊗ N (0, 1)
⊗ N (0, 1)
⊗ χ21
⊗ χ21

Methods
CH
KL
Hartigan
Silhouette
Gap
CH
KL
Hartigan
Silhouette
Gap
CH
KL
Hartigan
Silhouette
Gap
CH
KL
Hartigan
Silhouette
Gap
CH
KL
Hartigan
Silhouette
Gap
CH
KL
Hartigan
Silhouette
Gap

1
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0

2
3
11
0
20
20
1
2
0
9
20
0
2
0
16
10
0
2
0
7
19
0
4
0
11
17
19
5
0
9
20

3
1
1
7
0
0
0
4
11
11
0
0
5
7
4
10
3
3
14
13
1
0
1
11
8
3
0
3
5
5
0

Number of Clusters
4 5 6 7 8
1 0 2 1 0
0 1 1 2 2
4 6 3 0 0
0 0 0 0 0
0 0 0 0 0
0 0 2 0 1
1 1 2 1 6
1 3 1 2 1
0 0 0 0 0
0 0 0 0 0
0 1 0 1 1
2 0 0 2 5
4 1 0 3 2
0 0 0 0 0
0 0 0 0 0
0 1 0 3 1
3 4 1 6 0
3 0 0 1 1
0 0 0 0 0
0 0 0 0 0
0 0 0 2 2
5 4 3 0 1
5 0 1 0 0
1 0 0 0 0
0 0 0 0 0
1 0 0 0 0
3 1 2 2 0
13 0 2 0 0
6 0 0 0 0
0 0 0 0 0

9
5
1
0
0
0
5
1
0
0
0
3
2
1
0
0
6
0
0
0
0
6
0
2
0
0
0
3
0
0
0

10+
7
1
0
0
0
11
2
1
0
0
14
2
2
0
0
6
1
1
0
0
10
2
1
0
0
0
1
0
0
0

Table 5: Number of clusters detected in 20 trials for the Ward’s clustering algorithm

61

True Population Density

Methods

0.3 N (−4, 1) + 0.7 N (4, 1)

0.3 N (−3, 1) + 0.35 N (0, 1) + 0.35 N (3, 1)

0.3 t1 (−3) + 0.35 t1 (0) + 0.35 t1 (3)

0.3 dexp(−3) + 0.35 dexp(0) + 0.35 dexp(3)



Beta(8, 2) + Beta(5, 5) + Beta(2, 8)

{0.5 N (−2, 1) + 0.5 N (2, 1)}
⊗ N (0, 1)
⊗ N (0, 1)
⊗ χ21
⊗ χ21


3

CH
KL
Hartigan
Silhouette
Gap
CH
KL
Hartigan
Silhouette
Gap
CH
KL
Hartigan
Silhouette
Gap
CH
KL
Hartigan
Silhouette
Gap
CH
KL
Hartigan
Silhouette
Gap
CH
KL
Hartigan
Silhouette
Gap

1
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0

2
17
1
0
20
20
1
2
0
5
17
1
3
0
20
6
0
1
0
11
19
0
2
0
5
19
7
5
0
20
20

Number of Clusters
3 4 5 6 7 8
0 0 0 0 0 1
0 1 2 1 1 1
0 1 2 4 0 4
0 0 0 0 0 0
0 0 0 0 0 0
1 0 0 0 1 3
1 4 2 2 0 4
12 2 2 1 2 1
3 8 3 1 0 0
3 0 0 0 0 0
1 0 0 5 3 1
4 3 5 2 1 2
12 4 2 0 0 0
0 0 0 0 0 0
14 0 0 0 0 0
0 0 1 2 1 2
4 2 2 3 1 4
8 2 5 4 0 0
0 0 3 3 3 0
1 0 0 0 0 0
1 0 0 0 2 1
1 5 1 1 3 0
11 6 2 1 0 0
12 3 0 0 0 0
1 0 0 0 0 0
4 3 3 1 0 1
4 2 1 0 3 2
4 5 2 2 3 0
0 0 0 0 0 0
0 0 0 0 0 0

9
1
9
7
0
0
6
3
0
0
0
2
0
0
0
0
5
2
1
0
0
5
5
0
0
0
1
1
2
0
0

10+
1
4
2
0
0
8
2
0
0
0
7
0
2
0
0
9
1
0
0
0
11
2
0
0
0
0
2
2
0
0

Table 6: Number of clusters detected in 20 trials for the centroid clustering algorithm

62

True Population Density

Methods

0.3 N (−4, 1) + 0.7 N (4, 1)

0.3 N (−3, 1) + 0.35 N (0, 1) + 0.35 N (3, 1)

0.3 t1 (−3) + 0.35 t1 (0) + 0.35 t1 (3)

0.3 dexp(−3) + 0.35 dexp(0) + 0.35 dexp(3)



Beta(8, 2) + Beta(5, 5) + Beta(2, 8)

{0.5 N (−2, 1) + 0.5 N (2, 1)}
⊗ N (0, 1) ⊗ N (0, 1)
⊗ χ21 ⊗ χ21


3

BS and GMM-BIC
BS and GMM-merge
BS and CLARA
PS and CLARA
BS and PAM
PS and PAM
BS and GMM-BIC
BS and GMM-merge
BS and CLARA
PS and CLARA
BS and PAM
PS and PAM
BS and GMM-BIC
BS and GMM-merge
BS and CLARA
PS and CLARA
BS and PAM
PS and PAM
BS and GMM-BIC
BS and GMM-merge
BS and CLARA
PS and CLARA
BS and PAM
PS and PAM
BS and GMM-BIC
BS and GMM-merge
BS and CLARA
PS and CLARA
BS and PAM
PS and PAM
BS and GMM-BIC
BS and GMM-merge
BS and CLARA
PS and CLARA
BS and PAM
PS and PAM

1
0
0
0
0
0
0
0
0
0
0
0
5
0
0
0
6
0
20
0
0
0
0
0
11
0
0
0
1
0
20
0
0
0
0
0
15

2
20
20
20
20
20
6
0
0
1
0
2
15
17
20
1
7
3
0
1
0
0
0
0
9
0
20
12
9
7
0
0
20
20
20
5
5

Number of Clusters
3 4 5 6 7 8
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
14 0 0 0 0 0
19 0 0 0 0 0
20 0 0 0 0 0
19 0 0 0 0 0
20 0 0 0 0 0
18 0 0 0 0 0
0 0 0 0 0 0
3 0 0 0 0 0
0 0 0 0 0 0
0 3 8 3 4 1
2 4 1 0 0 0
0 0 0 2 3 6
0 0 0 0 0 0
5 0 4 5 2 0
20 0 0 0 0 0
20 0 0 0 0 0
20 0 0 0 0 0
20 0 0 0 0 0
0 0 0 0 0 0
1 2 1 5 2 3
0 0 0 0 0 0
5 0 1 0 0 0
9 1 0 0 0 0
12 0 1 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 1
0 0 0 0 0 0

9
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
6
0
3
0
0
0
0
0
1
0
2
0
0
0
0
0
0
0
0
0

10+
0
0
0
0
0
0
1
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
5
0
0
0
0
0
20
0
0
0
14
0

Table 7: Number of clusters detected in 20 trials by the prediction strength (PS) and bootstrap stability (BS) methods for selecting the number of clusters. Four clustering algorithms, CLARA, PAM, GMM-BIC and GMM-merge, were used in combination with each
of the two aforementioned number of clusters methods whenever the two parts were compatible in the R-package fpc.

63

contains the plots of the densities used in our numerical experiments. In the second sub-plot
of the top row we have the true population densities for the second, third and fourth simulation set-ups. In red, black and blue we have mixture of normals, t and double-exponential
densities, respectively. Bootstrap stability with CLARA or with the GMM-merge algorithm
and Prediction Strength with CLARA produced performances that are similar to those seen
with the k-means clustering algorithm. The other combinations were not as effective.
6.5.3

Performance on Large Data Sets

In Table 8 we report the performance of the BMT across 4 different simulation examples
involving large samples generated from Gaussian mixtures. For each example, 100 independent data sets are generated from the population density and the distribution of the
number of clusters detected by the BMT is reported (the frequencies are in parenthesis).
We also report the average and the standard deviation of the Mean Square Error (MSE)
over the cases where the the true number of clusters is correctly detected. The oracle MSE
is calculated based on the partition that uses the minima of the true population density.
Population Density

Sample
Size

Time in Sec
per Replicate

Number of Clusters


N (−2.5, 1) + N (0, 1) + N (2.5, 1) /3

104

0.70

0.5 N (−5, 1) + 0.25 N (0, 1) + 0.25 N (5, 1)

5 × 104

N (−1.1, 1) + N (1.1, 1) /2


N (0, 1) ± N (4, 1) ± N (8, 1) /5



MSE
Mean
SD

Oracle
MSE

2 (3), 3 (96), 4 (1)

0.6817

0.0405

0.6564

11.58

3 (100)

1.1152

0.0125

0.9769

105

44.07

1 (25), 2 (70), 3 (5)

0.6905

0.0258

0.6789

105

38.98

5 (100)

0.8909

0.0036

0.8909

Table 8: Performance of the BMT on simulated datasets of large sample sizes
Overall, the BMT correctly detects the true number of clusters with high certainly.
Also, the average MSE is observed to be very close to the Oracle one. To demonstrate
the scalability of the proposed method, we report the average elapsed time (in seconds)
per replication. The numerical experiments were performed at the Center for High Performance Computing (http://hpcc.usc.edu) of the University of Southern California.
The computations were done in R version 3.1.1 on Dual Quadcore Intel Xeon 2.33 GHz,
16GB Memory nodes. We used the Snowfall package of Knaus (2013) to distribute
computations over 100 CPUs.
64

Most of the competitor methods discussed above fail to accommodate the large sample
sizes of the data sets in Table 8. To get a better understanding of the relative scalibility of
the BMT, we implemented the HMAC algorithm of Li et al. 2007, which can handle large
sample sizes. HMAC is a non-parametric mode identification based clustering approach,
which produces clusters in hierarchical levels. Compared to HMAC, the BMT was found
to be substantially faster in producing the entire hierarchical path. On an iMac desktop with
2.9 GHz Intel Core i5 processor and 8 GB memory, for grids of 20 hierarchical levels, the
R-code associated with the HMAC paper required approximate run times of 3 minutes, 70
minutes and 4 hours for the sample sizes of 10000, 50000 and 100000, respectively. On the
other hand, the corresponding average run times for the BMT were 10 seconds, 5 minutes
and 20 minutes, respectively.
6.5.4

Performance for smaller sample sizes and the choice of α

For the numerical experiments in Section 4 of the main paper we used relatively large sample sizes and implemented the BMT with the threshold α uniformly set at 10%. The choice
of the threshold α is, however, important in smaller sample sizes. Our theoretical results
suggest that α should, in general, be increased when the sample size decreases. For this
purpose, we have conducted an extensive simulation study, where, for a variety of sample
sizes (100, 500, 1000, 2000, 5000 and 10000), and for a number of different threshold
values (1%, 2%, 5%, 7.5%, 10%, 15%, 20% and 25%), we record the percentage of cases
where the BMT approach produces splits, despite the fact that the true population density
is unimodal. We would like the threshold α to be sufficiently high, so that the instances of
splits in unimodal cases (which represent false discoveries) are well controlled. We tried 4
different population densities: (a) standard normal, (b) t with 1 df, (c) exponential with
rate parameter 1, and (d) cauchy. The results are reported in Table 9. We found that the
Gaussian case is the most difficult, requiring the highest threshold. Based on the simulation
results reported in Table 9, we suggest threshold sizes of 20% and 15% for sample sizes
100 and 500, respectively, and 10% or, possibly, lower for sample sizes 1000 and larger.
We also report the results for modality detection and number of cluster estimation for
smaller sample sizes of 500 and 1000. The comparisons with modality detection methods
are reported in tables 10 and 11, and the results on estimating the number of clusters by
BMT are presented in tables 12, 13, 14 and 15. The underlying distributions in these

65

tables exactly match those reported in the simulation study in the main paper. We report
the performance of our method for a number of different threshold sizes. In the modality
detection experiments, our proposed BMT approach performed better than the Dip and
Silverman tests in detecting modes, though it made a few more false discoveries. We
found that even in low sample sizes, BMT generally outperformed other number of clusters
approaches for finding the number of Beta and t mixture components (cases III and V).
Table 9: The table reports the percentage of cases where splits were detected by the BMT
algorithm, as the sample size, n, and the threshold size, α, are varied. We consider 4
unimodal densities: standard normal, t with 1 degrees of freedom, exponential with a unit
rate parameter and Cauchy.
Density

NORMAL

t with 1 df

sample size
2

5

7.5

10

15

20

25

100

100

99

97

81

42

16

5

5

500

100

75

44

28

14

6

1

1

1000

99

59

36

28

16

6

2

2

2000

80

20

9

5

2

1

1

1

5000

28

2

1

1

1

1

0

0

10000

3

0

0

0

0

0

0

0

100

100

64

39

23

14

5

2

2

500

54

5

3

0

0

0

0

0

1000

19

1

0

0

0

0

0

0

2000

0

0

0

0

0

0

0

0

5000

0

0

0

0

0

0

0

0

10000

EXP

CAUCHY

Threshold size (in %)
1

0

0

0

0

0

0

0

0

100

100

97

82

63

38

23

10

10

500

99

42

25

9

3

0

0

0

1000

78

17

3

0

0

0

0

0

2000

30

2

0

0

0

0

0

0

5000

6

0

0

0

0

0

0

0

10000

0

0

0

0

0

0

0

0

100

100

64

37

17

7

4

1

1

500

59

7

3

0

0

0

0

0

1000

15

0

0

0

0

0

0

0

2000

1

0

0

0

0

0

0

0

5000

0

0

0

0

0

0

0

0

10000

0

0

0

0

0

0

0

0

66

Table 10: Results comparing various multi-modality detection methods, across the five
different simulation scenarios considered in Table 1, are reported for the sample size of
500. The results for the BMT algorithm with threshold sizes 5%, 10%, 15% and 20% are
reported.
Dip Test P-value (D)

Silverman Test P-value (S)

BMT with threshold size
% multi-mode

Mean (D)

Std (D)

% multi-mode

Mean (S)

Std (S)

% multi-mode

5%

10%

15%

20%

Case I

0.9

0.16

0

0.51

0.25

0

71

35

12

7

Case II

0.86

0.18

0

0.59

0.26

1

85

46

27

15

Case III

0.8

0.24

0

0.39

0.27

7

97

74

55

40

Case IV

0.72

0.26

1

0.32

0.24

19

100

89

65

50

Case V

0.6

0.28

0

0.23

0.19

14

100

96

90

68

Table 11: Results comparing various multi-modality detection methods, across the five different simulation scenarios considered in Table 1, are reported for the sample size of 1000.
The results for the BMT algorithm with threshold sizes 5%, 10%, 15% and 20% are reported.
Dip Test P-value (D)

Silverman Test P-value (S)

BMT with threshold size
% multi-mode

Mean (D)

Std (D)

% multi-mode

Mean (S)

Std (S)

% multi-mode

5%

10%

15%

20%

Case I

0.94

0.11

0

0.48

0.22

0

45

12

9

4

Case II

0.92

0.12

0

0.61

0.25

0

65

26

14

4

Case III

0.8

0.22

0

0.33

0.28

20

97

70

57

47

Case IV

0.76

0.25

1

0.35

0.26

16

99

83

57

41

Case V

0.5

0.31

2

0.2

0.19

19

100

98

92

74

67

Table 12: The table reports the number of clusters detected by the BMT algorithm, with 3 different
prefixed threshold choices, in 100 trials for six simulation scenarios of Table 2, with sample sizes
being kept fixed at 500 for all the experiments.

Threshold

Simulation

size

Scenario

1

2

3

4

5

6

7

8

9

10+

I

0

53

41

6

0

0

0

0

0

0

II

0

9

65

26

0

0

0

0

0

0

III

0

12

84

4

0

0

0

0

0

0

IV

0

0

100

0

0

0

0

0

0

0

V

0

3

77

20

0

0

0

0

0

0

VI

0

27

15

20

0

19

0

6

6

7

I

0

87

13

0

0

0

0

0

0

0

II

1

35

64

0

0

0

0

0

0

0

III

1

36

63

0

0

0

0

0

0

0

IV

0

8

92

0

0

0

0

0

0

0

V

0

17

83

0

0

0

0

0

0

0

VI

0

67

1

26

0

1

0

5

0

0

I

0

100

0

0

0

0

0

0

0

0

II

7

61

32

0

0

0

0

0

0

0

III

18

63

19

0

0

0

0

0

0

0

IV

2

41

57

0

0

0

0

0

0

0

V

0

34

66

0

0

0

0

0

0

0

VI

0

83

16

0

0

0

0

1

0

0

10%

15%

20%

Number of clusters

68

Table 13: The table reports the number of clusters detected by the BMT algorithm, with 3 different
prefixed threshold choices, in 100 trials for six simulation scenarios of Table 2, with sample sizes
being kept fixed at 1000 for all the experiments.

Threshold

Simulation

Size

Scenario

1

2

3

4

5

6

7

8

9

10+

I

0

68

32

0

0

0

0

0

0

0

II

0

9

83

8

0

0

0

0

0

0

III

0

8

88

4

0

0

0

0

0

0

IV

0

1

99

0

0

0

0

0

0

0

V

0

6

78

16

0

0

0

0

0

0

VI

0

51

14

15

0

14

0

1

1

4

I

0

94

6

0

0

0

0

0

0

0

II

0

34

65

1

0

0

0

0

0

0

III

1

36

63

0

0

0

0

0

0

0

IV

0

5

95

0

0

0

0

0

0

0

V

0

11

89

0

0

0

0

0

0

0

VI

0

82

0

14

0

2

0

2

0

0

I

0

100

0

0

0

0

0

0

0

0

II

0

62

38

0

0

0

0

0

0

0

III

15

67

18

0

0

0

0

0

0

0

IV

0

36

64

0

0

0

0

0

0

0

V

0

27

73

0

0

0

0

0

0

0

VI

0

94

6

0

0

0

0

0

0

0

10%

15%

20%

Number of clusters

69

Table 14: Number of clusters detected in 100 trials for simulation scenarios of Table 2,
with the sample size of 500.
Simulation

Methods

Scenario

I

II

III

IV

V

VI

Number of clusters
1

2

3

4

5

6

7

8

9

10+

CH

0

33

9

3

6

4

4

8

9

24

KL

0

58

4

8

3

3

4

11

4

5

Hartigan

0

0

38

10

12

9

8

8

5

10

Silhouette

0

100

0

0

0

0

0

0

0

0

Gap

0

100

0

0

0

0

0

0

0

0

Jump

0

100

0

0

0

0

0

0

0

0

Pred Str.

0

100

0

0

0

0

0

0

0

0

Stability

0

100

0

0

0

0

0

0

0

0

CH

0

0

4

1

1

4

6

14

28

42

KL

0

10

18

13

17

10

9

8

8

7

Hartigan

0

0

58

17

9

5

2

4

2

3

Silhouette

0

37

63

0

0

0

0

0

0

0

Gap

0

100

0

0

0

0

0

0

0

0

Jump

6

0

76

0

0

6

3

7

1

1

Pred Str.

0

2

98

0

0

0

0

0

0

0

Stability

0

19

81

0

0

0

0

0

0

0

CH

0

8

7

1

3

5

6

14

20

36

KL

0

26

13

12

11

6

8

10

9

5

Hartigan

0

0

39

12

11

13

8

6

7

4

Silhouette

0

72

26

0

0

0

1

0

1

0

Gap

0

60

40

0

0

0

0

0

0

0

Jump

18

12

6

5

4

5

9

13

17

11

Pred Str.

10

53

36

1

0

0

0

0

0

0

Stability

0

47

9

0

0

0

6

5

21

12

CH

0

2

4

5

5

5

7

8

20

44

KL

0

17

6

9

8

11

13

11

13

12

Hartigan

0

0

54

18

12

6

4

4

1

1

Silhouette

0

8

92

0

0

0

0

0

0

0

Gap

0

100

0

0

0

0

0

0

0

0

Jump

67

0

27

0

0

0

1

0

2

3

Pred Str.

0

0

100

0

0

0

0

0

0

0

Stability

0

4

96

0

0

0

0

0

0

0

CH

0

0

2

3

2

4

5

15

28

41

KL

0

13

18

10

14

13

9

6

8

9

Hartigan

0

0

57

17

10

7

2

3

2

2

Silhouette

0

54

46

0

0

0

0

0

0

0

Gap

0

100

0

0

0

0

0

0

0

0

Jump

0

0

62

2

0

6

8

8

4

10

Pred Str.

13

77

10

0

0

0

0

0

0

0

Stability

0

43

57

0

0

0

0

0

0

0

CH

0

1

4

1

3

9

16

24

42

KL

0

13

14

14

14

10

8

7

14

6

Hartigan

0

0

41

24

13

9

4

1

7

1

Silhouette

0

5

53

19

0

9

4

5

3

2

Gap

0

100

0

0

0

0

0

0

0

0

Jump

0

2

0

13

0

0

2

9

24

50

Pred Str.

0

100

0

0

0

0

0

0

0

0

Stability

0

68

27

5

0

0

0

0

0

0

70

Table 15: Number of clusters detected in 100 trials for simulation scenarios of Table 2,
with the sample size of 1000.
Simulation

Methods

Scenario

I

II

III

IV

V

VI

Number of clusters
1

2

3

4

5

6

7

8

9

10+

CH

0

31

4

3

4

4

6

5

15

28

KL

0

46

11

6

4

5

3

4

12

9

Hartigan

0

0

25

15

14

10

2

15

7

12

Silhouette

0

100

0

0

0

0

0

0

0

0

Gap

0

100

0

0

0

0

0

0

0

0

Jump

0

100

0

0

0

0

0

0

0

0

Pred Str.

0

100

0

0

0

0

0

0

0

0

Stability

0

100

0

0

0

0

0

0

0

0

CH

0

0

2

0

3

5

4

14

25

47

KL

0

11

3

12

10

11

19

14

13

7

Hartigan

0

0

55

15

15

7

1

4

3

0

Silhouette

0

27

73

0

0

0

0

0

0

0

Gap

0

100

0

0

0

0

0

0

0

0

Jump

4

0

96

0

0

0

0

0

0

0

Pred Str.

0

0

100

0

0

0

0

0

0

0

Stability

0

14

86

0

0

0

0

0

0

0

CH

0

10

10

6

0

5

12

11

12

34

KL

0

19

14

11

11

9

12

8

6

10

Hartigan

0

0

39

20

9

10

8

5

4

5

Silhouette

0

72

28

0

0

0

0

0

0

0

Gap

0

60

40

0

0

0

0

0

0

0

Jump

22

12

11

4

6

2

19

4

10

10

Pred Str.

3

40

51

6

0

0

0

0

0

0

Stability

0

64

14

0

0

0

0

0

9

13

CH

0

2

4

0

1

3

5

16

20

49

KL

0

6

10

6

13

13

11

10

14

17

Hartigan

0

0

53

9

15

6

8

3

3

3

Silhouette

0

0

100

0

0

0

0

0

0

0

Gap

0

100

0

0

0

0

0

0

0

0

Jump

80

0

20

0

0

0

0

0

0

0

Pred Str.

0

0

100

0

0

0

0

0

0

0

Stability

0

2

98

0

0

0

0

0

0

0

CH

0

0

1

1

1

5

12

22

16

42

KL

0

12

18

14

12

12

14

3

8

7

Hartigan

0

0

70

10

4

5

5

3

2

1

Silhouette

0

53

47

0

0

0

0

0

0

0

Gap

0

100

0

0

0

0

0

0

0

0

Jump

0

0

90

0

0

0

2

1

4

3

Pred Str.

0

67

32

1

0

0

0

0

0

0

Stability

0

38

62

0

0

0

0

0

0

0

CH

0

0

3

2

0

7

8

15

21

44

KL

0

13

15

7

11

13

13

9

11

8

Hartigan

0

0

42

26

8

9

9

2

1

3

Silhouette

0

1

69

13

1

13

2

0

1

0

Gap

0

100

0

0

0

0

0

0

0

0

Jump

0

8

0

14

0

0

0

9

21

48

Pred Str.

0

100

0

0

0

0

0

0

0

0

Stability

0

79

10

11

0

0

0

0

0

0

71

6.5.5

BMT & Sub-population Analysis in Single Cell Virology

We demonstrate an application of our clustering method in an immunology study conducted
at single cell level. Emerging technologies (Wang and Bodivitz, 2011) have recently enabled us to collect proteomic data sets at single cell resolution. These data sets reflect the
variations of protein expressions across cells and need clustering techniques for detection
of cellular sub-populations. Typically, sub-populations are detected by core-protein expressions based cluster analysis of the samples, and the signaling expressions of the resultant
sub-populations are subsequently studied. In Figures 6 and 7 we display the results of the
BMT induced clustering on the virology datasets of Sen et al. (2014). Figure 8 shows the
post-clustering, sub-population level signaling expressions.

Acknowledgement
Radchenko’s research was partially supported by NSF Grant DMS-1209057. Computing
resources were provided by the USC Center for High Performance Computing. The R code
and the data sets used for producing the numerical results of Section 4 can be downloaded
from http://www-bcf.usc.edu/˜gourab/code-bmt.

References
Aitkin, M. and Rubin, D. B. (1985). Estimation and hypothesis testing in finite mixture
models. Journal of the Royal Statistical Society. Series B (Methodological) 67–75.
Bach, F. R. and Harchaoui, Z. (2008). Diffrac: a discriminative and flexible framework for
clustering. In Advances in Neural Information Processing Systems, 49–56.
Belkin, M. and Niyogi, P. (2001). Laplacian eigenmaps and spectral techniques for embedding and clustering. In NIPS, vol. 14, 585–591.
Bendall, S. C., Simonds, E. F., Qiu, P., Amir, E., Krutzik, P. O., Finck, R., Bruggner, R. V.,
Melamed, R., Trejo, A., Ornatsky, O. I., Balderas, R. S., Plevritis, S. K., Sachs, K., Pe’er,
D., Tanner, S. D., and Nolan, G. P. (2011). Single-cell mass cytometry of differential
immune and drug responses across a human hematopoietic continuum. Science (New
York, N.Y.) 332, 6030, 687–696.
Bondell, H. D. and Reich, B. J. (2008). Simultaneous regression shrinkage, variable selection, and supervised clustering of predictors with oscar. Biometrics 64, 1, 115–123.

72

Figure 6: Across columns we have histograms of the expression values of the proteins CD4,
CD8 and CD45RA, respectively. Along rows, from top to bottom, we have the histograms
of the Uninfected and Bystander population, respectively, for the independent experiments
I-III. The shaded gray region denotes unexpressed values. Splits in the expression values
(if any) detected by BMT are shown by vertical red lines.

73

Figure 7: Across columns we have histograms of the expression values of CD4, CD8 and
CD45RA. Across rows are the histograms of UN and BY populations for Experiments
IV-V. Splits detected by BMT (if any) are shown by red lines.

Figure 8: The above plot shows the heatmaps of the protein expression values (in order of decreasing intensity: Red, Yellow, Green and Black) of the Uninfected and Bystander populations in Experiment I. The horizontal white lines demarcate the five major
sub-population detected by BMT algorithm, based on the expression of the three surface
markers on the left of the vertical white line. The proteins on the right of the vertical line
are associated with cell-signaling. The heatmaps are standardized separately for the two
populations.

74

Caliński, T. and Harabasz, J. (1974). A dendrite method for cluster analysis. Communications in Statistics-theory and Methods 3, 1, 1–27.
Charrad, M., Ghazzali, N., Boiteau, V., and Niknafs, A. (2014). NbClust: An R package for
determining the relevant number of clusters in a data set. Journal of Statistical Software
61, 6, 1–36.
Chi, E. C. and Lange, K. (2013). Splitting methods for convex clustering. arXiv preprint
arXiv:1304.0499 .
Fang, Y. and Wang, J. (2012). Selection of the number of clusters via the bootstrap method.
Computational Statistics & Data Analysis 56, 3, 468–477.
Fraley, C. and Raftery, A. E. (2002). Model-based clustering, discriminant analysis, and
density estimation. Journal of the American statistical Association 97, 458, 611–631.
Gordon, A. D. (1996). Null models in cluster validation. In From data to knowledge,
32–44. Springer.
Hall, P. and York, M. (2001). On the calibration of silverman’s test for multimodality.
Statistica Sinica 11, 2, 515–536.
Hartigan, J. (1978). Asymptotic distributions for clustering criteria. The Annals of Statistics
117–131.
Hartigan, J. A. (1975). Clustering algorithms. Wiley.
Hartigan, J. A. and Hartigan, P. (1985). The dip test of unimodality. The Annals of Statistics
70–84.
Hennig, C. (2010). Methods for merging gaussian mixture components. Advances in data
analysis and classification 4, 1, 3–34.
Hennig, C. (2014). fpc: Flexible procedures for clustering. R package version 2.1-9.
Hocking, T., Vert, J.-P., Bach, F., and Joulin, A. (2011). Clusterpath: an algorithm for
clustering using convex fusion penalties. In L. Getoor and T. Scheffer, eds., ICML, 745–
752. Omnipress.
Hoefling, H. (2010). A path algorithm for the fused lasso signal approximator. Journal of
Computational and Graphical Statistics 19, 4, 984–1006.
Jain, A. K. (2010). Data clustering: 50 years beyond k-means. Pattern Recognition Letters
31, 8, 651–666.
Kaufman, L. and Rousseeuw, P. (1986). Clustering large data sets. Elsevier.
Kaufman, L. and Rousseeuw, P. J. (1990). Partitioning around medoids (program pam).
Finding groups in data: an introduction to cluster analysis 68–125.

75

Kaufman, L. and Rousseeuw, P. J. (2009). Finding groups in data: an introduction to
cluster analysis, vol. 344. John Wiley & Sons.
Ke, T., Fan, J., and Wu, Y. (2013).
arXiv:1303.7409 .

Homogeneity in regression.

arXiv preprint

Kim, J. and Pollard, D. (1990). Cube root asymptotics. The Annals of Statistics 191–219.
Knaus, J. (2013). snowfall: Easier cluster computing (based on snow). R package version
1.84-6.
Krzanowski, W. J. and Lai, Y. (1988). A criterion for determining the number of groups in
a data set using sum-of-squares clustering. Biometrics 23–34.
Li, J. (2005). Clustering based on a multilayer mixture model. Journal of Computational
and Graphical Statistics 14, 3, 547–568.
Li, J., Ray, S., and Lindsay, B. G. (2007). A nonparametric statistical approach to clustering
via mode identification. Journal of Machine Learning Research 8, 8, 1687–1723.
Linderman, M. D., Bjornson, Z., Simonds, E. F., Qiu, P., Bruggner, R. V., Sheode, K.,
Meng, T. H., Plevritis, S. K., and Nolan, G. P. (2012). Cytospade: high-performance
analysis and visualization of high-dimensional cytometry data. Bioinformatics 28, 18,
2400–2401.
Lindsten, F., Ohlsson, H., and Ljung, L. (2011). Clustering using sum-of-norms regularization: With application to particle filter output computation. In Statistical Signal
Processing Workshop (SSP), 2011 IEEE, 201–204. IEEE.
MacQueen, J. et al. (1967). Some methods for classification and analysis of multivariate
observations. In Proceedings of the fifth Berkeley symposium on mathematical statistics
and probability, vol. 1, 281–297. California, USA.
Maechler, M. (2013). diptest: Hartigan’s dip test statistic for unimodality - corrected code.
R package version 0.75-5.
Milligan, G. W. and Cooper, M. C. (1985). An examination of procedures for determining
the number of clusters in a data set. Psychometrika 50, 2, 159–179.
Müller, D. W. and Sawitzki, G. (1991). Excess mass estimates and tests for multimodality.
Journal of the American Statistical Association 86, 415, 738–746.
Newell, E., Sigal, N., Bendall, S., Nolan, G., and Davis, M. (2012). Cytometry by timeof-flight shows combinatorial cytokine expression and virus-specific cell niches within a
continuum of cd8+ t cell phenotypes. Immunity 36, 1, 142 – 152.
Pelckmans, K., De Brabanter, J., Suykens, J., and De Moor, B. (2005). Convex clustering
shrinkage. In PASCAL Workshop on Statistics and Optimization of Clustering Workshop.

76

Pollard, D. (1981). Strong consistency of k-means clustering. The Annals of Statistics 9,
1, 135–140.
Pollard, D. (1982). A central limit theorem for k-means clustering. The Annals of Probability 919–926.
Qiu, P., Simonds, E. F., Bendall, S. C., Jr., K. D. G., Bruggner, R. V., Linderman, M. D.,
Sachs, K., Nolan, G. P., and Plevritis, S. K. (2011). Extracting a cellular hierarchy from
high-dimensional cytometry data with spade. Nature Biotechnology .
Rinaldo, A. et al. (2009). Properties and refinements of the fused lasso. The Annals of
Statistics 37, 5B, 2922–2952.
Roeder, K. (1994). A graphical technique for determining the number of components in a
mixture of normals. Journal of the American Statistical Association 89, 426, 487–495.
Rohe, K., Chatterjee, S., and Yu, B. (2011). Spectral clustering and the high-dimensional
stochastic blockmodel. The Annals of Statistics 1878–1915.
Sen, N., Mukherjee, G., Sen, A., Bendall, S., Sung, P., Nolan, G., and Arvin, A. (2014).
Single-cell mass cytometry analysis of human tonsil t cell remodeling by varicella zoster
virus. Cell Reports 8, 2, 633 – 645.
Shen, X. and Huang, H.-C. (2010). Grouping pursuit through a regularization solution
surface. Journal of the American Statistical Association 105, 490.
Shen, X., Huang, H.-C., and Pan, W. (2012). Simultaneous supervised clustering and
feature selection over a graph. Biometrika 99, 4, 899–914.
Shi, T., Belkin, M., and Yu, B. (2009). Data spectroscopy: Eigenspaces of convolution
operators and clustering. The Annals of Statistics 3960–3984.
Silverman, B. W. (1981). Using kernel density estimates to investigate multimodality.
Journal of the Royal Statistical Society. Series B (Methodological) 97–99.
Soltanolkotabi, M. and Candés, E. J. (2012). A geometric analysis of subspace clustering
with outliers. Ann. Statist. 40, 4, 2195–2238.
Sugar, C. A. and James, G. M. (2003). Finding the number of clusters in a dataset. Journal
of the American Statistical Association 98, 463.
Tan, K. M. and Witten, D. (2015). Statistical properties of convex clustering. arXiv preprint
arXiv:1503.08340 .
Thorndike, R. L. (1953). Who belongs in the family? Psychometrika 18, 4, 267–276.
Tibshirani, R. (1996). Regression shrinkage and selection via the lasso. Journal of the
Royal Statistical Society. Series B (Methodological) 267–288.

77

Tibshirani, R. (2011). Regression shrinkage and selection via the lasso: a retrospective.
Journal of the Royal Statistical Society: Series B (Statistical Methodology) 73, 3, 273–
282.
Tibshirani, R., Saunders, M., Rosset, S., Zhu, J., and Knight, K. (2005). Sparsity and
smoothness via the fused lasso. Journal of the Royal Statistical Society: Series B (Statistical Methodology) 67, 1, 91–108.
Tibshirani, R. and Walther, G. (2005). Cluster validation by prediction strength. Journal of
Computational and Graphical Statistics 14, 3, 511–528.
Tibshirani, R., Walther, G., and Hastie, T. (2001). Estimating the number of clusters in a
data set via the gap statistic. Journal of the Royal Statistical Society: Series B (Statistical
Methodology) 63, 2, 411–423.
Tibshirani, R. J. (2013). Adaptive piecewise polynomial estimation via trend filtering.
arXiv preprint arXiv:1304.2986 .
van der Vaart, A. W. (1998). Asymptotic Statistics. Cambridge University Press.
van der Vaart, A. W. and Wellner, J. A. (1996). Weak Convergence and Empirical Processes: With Applications to Statistics. Springer-Verlag.
Vollmer, S., Holzmann, H., and Schwaiger, F. (2013). Peaks vs components. Review of
Development Economics 17, 2, 352–364.
Wang, D. J. and Bodivitz, S. (2011). Single cell analysis: the new frontier in omics. Nature
Methods .
Wang, J. (2010a). Consistent selection of the number of clusters via crossvalidation.
Biometrika 97, 4, 893–904.
Wang, J. (2010b). Consistent selection of the number of clusters via crossvalidation.
Biometrika 97, 4, 893–904.
Ward Jr, J. H. (1963). Hierarchical grouping to optimize an objective function. Journal of
the American statistical association 58, 301, 236–244.
Witten, D. M. and Tibshirani, R. (2010). A framework for feature selection in clustering.
Journal of the American Statistical Association 105, 490.
Xu, L., Neufeld, J., Larson, B., and Schuurmans, D. (2004). Maximum margin clustering.
In Advances in neural information processing systems, 1537–1544.
Zerboni, L., Sen, N., Oliver, S. L., and Arvin, A. M. (2014). Molecular mechanisms of
varicella zoster virus pathogenesis. Nature Reviews Microbiology .
Zhu, C., Xu, H., Leng, C., and Yan, S. (2014). Convex optimization procedure for clustering: Theoretical revisit. In Advances in Neural Information Processing Systems, 1619–
1627.
78

