Hierarchical Clustering for Euclidean Data
Moses Charikar Vaggos Chatziafratis Rad Niazadeh Grigory Yaroslavtsev

























































 







all
arrows
chains,arrows




rnred     
vcblue    
msgreen   

plain


theoremTheorem[section]
lemma[theorem]Lemma
proposition[theorem]Proposition
corollary[theorem]Corollary
claim[theorem]Claim
fact[theorem]Fact
conj[theorem]Conjecture
definitionDefinition[section]

remark[theorem]Remark

proofof[1]









example[1][Example]
		
  #1






















𝔼
poly
Lap





































⌈⌉
⌊⌋





























arabic










Hierarchical Clustering for Euclidean Data
    Moses CharikarStanford University.
moses@cs.stanford.eduVaggos Chatziafratis[1]
vaggos@cs.stanford.edu Rad Niazadeh[1]
rad@cs.stanford.edu Grigory YaroslavtsevIndiana University, Bloomington. Research supported by NSF award CCF-1657477 and Facebook Faculty Research Award.
grigory@grigory.us
    December 30, 2023
=======================================================================================================================================================================================================================================================================================================





Recent works on Hierarchical Clustering (HC), a well-studied problem in exploratory data analysis, have focused on optimizing various objective functions for this problem under arbitrary similarity measures. In this paper we take the first step and give novel scalable algorithms for this problem tailored to Euclidean data in ℝ^d and under vector-based similarity measures, a prevalent model in several typical machine learning applications. We focus primarily on the popular Gaussian kernel and other related measures, presenting our results through the lens of the objective introduced recently by <cit.>. We show that the approximation factor in <cit.> can be improved for Euclidean data. We further demonstrate both theoretically and experimentally that our algorithms scale to very high dimension d, while outperforming average-linkage and showing competitive results against other less scalable approaches.




§ INTRODUCTION












Hierarchical clustering is a popular data analysis method, with various applications in data mining <cit.>, phylogeny <cit.>, and even finance <cit.>. In practice, simple agglomerative procedures like average-linkage, single-linkage and complete-linkage are often used for this task
(See the book by <cit.> for a comprehensive discussion). While the study of hierarchical clustering has focused on algorithms, there was a lack of objective functions to measure the quality of the output and compare the performance of algorithms.
To remedy this, <cit.> recently introduced and studied an interesting objective function for hierarchical clustering for similarity measures.



The input to the above problem is a set of points V with a similarity measure w_ij between i and j. 
Given a hierarchical clustering represented by a tree 𝒯 whose leaves correspond to points in V, Dasgupta's objective is defined as 
(𝒯) = ∑_(i,j) ∈ E w_ij |{x ∈ T(i,j)}|,
where T(i,j) is the subtree of 𝒯 rooted in the least common ancestor of i and j in 𝒯.   showed that solutions obtained from minimizing this objective has several desirable properties, which prompted a line of work on objective driven algorithms for hierarchical clustering, resulting in new algorithms as well as shedding light on the performance of classical methods <cit.>.





Recently, this 
viewpoint has been applied to understand the performance of average linkage agglomerative clustering, one of the most popular algorithms used in practice. 
<cit.> introduced a new objective, in some sense dual to the objective introduced by Dasgupta:
  
    (𝒯) = ∑_(i,j) ∈ E w_ij (|V| - |{x ∈ T(i,j)}|),

They showed that Average-Linkage obtains a 1/3-approximation for maximizing this objective function.





It turns out that a random solution 
also achieves a 1/3-approximation ratio for this problem.



Recently <cit.> showed that in the worst case, the approximation ratio achieved by Average-Linkage is no better than 1/3 for maximizing (𝒯).


They also gave an SDP based algorithm that achieves a (1/3+ϵ)-approximation for the problem, for small constant ϵ.






One drawback of the prior work on these hierarchical clustering objectives is the fact that they all consider arbitrary similarity scores (specified as an n × n matrix);
however, there is much more structure to such similarity scores in practice.
In this paper, we initiate the study of the commonly encountered case of Euclidean data, where the similarity score w_ij is computed by applying a monotone decreasing function to the Euclidean distance between i and j. Roughly speaking, we show how to exploit this structure to design improved approximation/faster algorithms for hierarchical clustering, and how to re-analyze algorithms commonly used in practice.
In this paper, we focus on the problem of maximizing (𝒯), introduced by <cit.>. 







Arguably the most common distance-based similarity measure used for Euclidean data is the Gaussian kernel. Here we use the spherical version w_ij∼exp(-v_i - v_j/2σ^2).  The parameter σ^2, referred to as the bandwidth, plays an important role for applications and a large body of literature exists on selection of this parameter <cit.>.





As discussed above, devising simple practical algorithms that improve on the 1/3-approximation for general similarity measures appears to be a major challenge as observed in <cit.>. 
In this paper, we show that this approximation can be improved for Euclidean data, through fast algorithms that can be scaled to very large datasets.

While it might seem that the case of Euclidean data with the Gaussian kernel is a very restricted class of inputs to the HC problem, we show that for suitably high dimensions and suitable small choice of bandwidth σ^2 it can simulate arbitrary similarity scores (scaled appropriately).
Thus any improvements to approximation guarantees for the Euclidean case (that do not apply to general similar scores) must necessarily involve assumptions about the dimension of the data (not very high) or on σ^2 (not  pathologically small).
Such assumptions on σ^2  are consistent with common methods for computing the bandwidth parameter based on data (e.g., <cit.>). 




  
Contributions.



We start with the simplest case of 1-dimensional Euclidean data. 
Even in this seemingly simple setting, obtaining an efficient algorithm that produces an exactly optimum solution seems non-trivial, motivating the study of approximation algorithms.
In Section <ref> we prove that two algorithms – Random Cut (RC) and Average-Linkage (AL) – obtain 1/2-approximation of the optimal solution, i.e., tree maximizing (𝒯) (AL achieves this deterministically, while RC only in expectation). This beats the best known approximation for the general case, which is 0.336 <cit.>. Here RC is substantially faster than AL, with a running time of  O(n log n) vs. O(n^2).


We next consider the high-dimensional case with the Gaussian Kernel and show that Average-Linkage cannot beat the factor 1/3 even in poly-logarithmic dimensions. 
We propose the Projected Random Cut (PRC) algorithm that gets a constant improvement over 1/3, irrespective of the dimension d (the improvement is a function of the ratio of the diameter to σ^2, and drops as this ratio gets large). 


Furthermore, PRC runs in O(n(d+log n)) time while Average-Linkage runs in O(dn^2) time. Even single-linkage runs in almost-linear time only for constant d and has exponential dependence on the dimension (see e.g., <cit.>) and it is open whether it can be scaled to large d when n is large.







  
Experiments. Many existing algorithms have time efficiency shortcomings, and none of them can be used for really large datasets. On the contrary, our Projected Random Cut (see <Ref>) is a fast HC algorithm that scales to the largest ML datasets. The running time scales almost linearly and the algorithm can be implemented in one pass without needing to store similarities, so the memory is O(n). We also evaluate its quality on a small dataset (Zoo).










  
Further related work.
HC has been extensively studied across different domains (we refer the reader to <cit.> for a survey). One killer application in biology is in phylogenetics <cit.> and actually popular algorithms, like Average-Linkage are originating in this field. HC is also widely used to perform community detection in social networks <cit.>, is used in bioinformatics <cit.> and even in finance <cit.>. Other important applications include image and text classification <cit.>.

After a formal optimization HC objective was introduced, many works studied HC from an approximation algorithms perspective. The original top-down Sparsest Cut algorithm by Dasgupta, was later shown that it achieves an  O(log n) approximation by <cit.>; this was finally improved to O(√(log n)) by <cit.>. Subsequent work on beyond-worst-case analysis <cit.> proposes a hierarchical stochastic block model under which an older spectral algorithm of <cit.> combined with linkage steps gets a O(1)-approximation. Finally, alongside the objective of scaling HC to large graphs, the work of <cit.> presents two MapReduce implementations with a focus on minimum-spanning-tree based clusterings that scales to graphs with trillions of edges.




  
Open problems.


Here is a list of open questions that we remain unanswered in this paper:
(1) Can one get an improvement over 1/3 for the problem of maximizing (𝒯) as a function of d for small d?

(2) Can projection on low-dimensional subspaces be used to improve the approximation ratio for the high-dimensional case further?

(3) Does Average-Linkage achieve a 3/4-approximation for 1-dimensional data?





§ PRELIMINARIES



  
Euclidean data. We consider data sets represented as sets of d-dimensional feature vectors. Suppose these vectors are v_1,…,v_n∈ℝ^d. We focus on similarity measures between pairs of data points, denoted by [w_ij]_i,j∈[n], where the similarities only depend on the underlying vectors, w_ij=(v_i,v_j) for some function :ℝ^d×ℝ^d→ [0,1] and furthermore are fully determined by monotone functions of distances between them.


	A similarity measure w_ij = (v_i, v_j) is “distance-based” if (v_i, v_j) = (v_i - v_j_2) for some function ℝ→ [0,1], and is “monotone distance-based” if furthermore ℝ→ [0,1] is a monotone non-increasing function.


As an example of the monotone similarity measure it is natural to consider the Gaussian kernel similarity, 
    w_ij = (√(2 π)σ)^-n e^- v_i - v_j_2^2/2σ^2,Gaussian Kernel

where σ is a normalization factor determining the bandwidth of the Gaussian kernel <cit.>. For simplicity, we ignore the multiplicative factor (√(2 π)σ)^-n (unless noted otherwise), as our focus is on multiplicative approximations and scaling has no effects. 



  
Linkage-based hierarchical clustering. Among various algorithms popular in practice for HC, we focus on two very popular ones: single-linkage and average-linkage, which are two simple agglomerative clustering algorithms that recursively merge pairs of nodes (or super nodes) to create the final binary tree. Single-linkage picks the pair of super nodes with the maximum similarity weight between their data points, merges A and B maximizing max_i∈ A,j∈ Bw_ij.  On the contrary, average-linkage picks the pair of super-nodes with maximum average similarity weight at each step, merges  A and B maximizing w_AB/| A|·| B |, where w_AB:=∑_i∈ A, j∈ Bw_ij. Average-linkage has an approximation ratio of 1/3 for maximizing the  objective function <cit.>, and this factor is tight <cit.>.



  
General upper bounds. In order to analyze the linkage-based clustering algorithms and our proposed algorithms, we propose a natural upper bound on the value of the  objective function. The idea is to decompose the objective function  into contributions of triple of vertices: 

    (𝒯)=∑_i<j<k(w_ij^k_ij+w_jk^i_jk+w_ik^j_ik)

where ^x_yz denotes the event that x is separated first among the vertices of triple {x,y,z} in tree 𝒯. Note that the final tree scores only one of the similarity weights between the triple {i,j,k}. Given this observation, we define the following benchmark:

    ≜∑_i<j<kmax(w_ij,w_jk,w_ik),

Clearly, for all trees 𝒯,  (𝒯)≤.



  
One-dimensional benchmarks. Consider the special case of 1D data points (with any monotone distance-based similarity measure as in <Ref>), and suppose v_1≤ v_2≤…≤ v_n∈ℝ. Now, for any triple i<j<k,  as a simple observation we have w_ik≤min(w_ij,w_jk). Hence we can modify the above benchmark to obtain two refined new benchmarks:

    ≜∑_i<j<kmax(w_ij,w_jk),
       ≜∑_i<j<k(w_ij+w_jk).

Again, clearly for all trees 𝒯 we have: 

    (𝒯)≤≤







§ HIERARCHICAL CLUSTERING IN 1D




In this section we look at the extreme case where the feature vectors have d=1, and we try to analyze the popular/natural algorithms existing in this domain by evaluating how well they approximate the objective function . We focus on average-linkage and Random Cut (will be formally defined later). In particular, Random Cut is a building block of our algorithm for high-dimensional data given in Section <ref>.

We show (1) average-linkage gives a 1/2-approximation, and can obtain no better than 3/4 fraction of the optimal objective value in the worst-case. We then show (2) Random Cut also is a 1/2-approximation (in expectation) and this factor is tight.  
In the Appendix <ref>, we discuss other simple algorithms: single-linkage and greedy cutting (will be defined formally later). We start by the observation that greedy cutting and single-linkage output the same tree (and so are equivalent). We finish by showing (3) there is an instance where single-linkage attains only 1/2 of the optimum objective value. We further show on this instance both average-linkage and random cutting are almost optimal.

For the rest of this section, suppose we have 1D points x_1≤…,≤ x_n∈ℝ, where w_ij=(‖ x_i-x_j ‖) for some monotone non-increasing function :ℝ→ [0,1].



  
Average-linkage. 
In order to simplify the notation, we define w_AB to be ∑_i∈ A,j∈ Bw_ij for any two sets A and B. We first prove a simple structural property of average-linkage in 1D.


For d=1, under any monotone distance-based similarity measure w_ij = (‖ x_i-x_j ‖),  average-linkage can always merge neighbouring clusters. 


We do the proof by induction. This property holds at the first step of average-linkage (base of induction). Now suppose up to a particular step of average-linkage, the algorithm has always merged neighbours. At this step, we have an ordered collection of super nodes (C_1,…,C_m), where for every i<j and for all (x,y)∈ C_i× C_j, x≤ y. If at this step average-linkage doesn't merge two neighbouring super nodes, then there exists i<j<k, where 

    w_C_iC_j/| C_i|| C_j| <w_C_iC_k/| C_i|| C_k|⇒∑_x∈ C_iw_xC_j/| C_j|<∑_x∈ C_iw_xC_k/| C_k|

But note that because of the monotone non-increasing distance-based similarity weights, for any triple (x,y,z)∈ C_i× C_j× C_k we have w_xy≥ w_xz, This is a contradiction to the above inequality, which finishes the inductive step.



	For d = 1, under any monotone distance-based similarity measure w_ij = (‖ x_i-x_j ‖), average-linkage obtains at least 12 of the , and hence is a 1/2-approximation for the objective .



 The proof uses a potential function argument. Given a partitioning of points x_1≤ x_2≤…≤ x_n into sets S_1,…,S_m, a triple of points i<j<k is called separated if no pair of these three points belong to the same set. 
Now, the potential function  gets {S_1,…,S_m} as input, and maps it to a summation over all separated triples by {S_i}_i∈[m] as below:



    (S_1,…,S_m)=∑_i<j<k: (i,j,k) is separated(w_ij+w_jk)

Note that ({x_1},…,{x_n})=, and ({x_1,…,x_n})=0. 

We now run average-linkage. Based on the definition of , every time that average-linkage merges two super nodes A and B it scores w_AB· (n-| A| -| B|), and sum over of all these per-step scores is equal to its final objective value. Let  denotes the variable that stores of the score of average-linkage over time. At every step we keep track of (1) the change in the potential function, denoted by Δ and (2) how much progress average-linkage had towards the final objective value, denoted by Δ. In order to prove 1/2-approximation, it suffices to show that at every step of average-linkage we have:

    Δ+1/2Δ≥ 0. ⋆

To see this note that average-linkage starts with all points separated, and ends with one cluster/super node with all the points. Therefore, by summing <ref> over all merging steps of average-linkage and canceling the terms in the telescopic sum, we have: 

    ((𝒯_AL)-0)+
       1/2(({x_1,…,x_n})-({x_1},…,{x_n}))≥ 0.

Plugging values of  at the start and the end, we get:

    (𝒯_AL)≥1/2(),

which implies the 1/2-approximation factor. 
	
	To prove <ref>, we focus on a single step of average-linkage where by <Ref> some two neighboring clusters denoted as A and B get merged. Let C denote the nodes on the left of A and let D denote the nodes on the right of B (<Ref>).[Note that C,D may be empty sets.] By merging two clusters A,B, the change in the score of average-linkage is:

    Δ = w_AB(| C|+| D|)
	
Moreover, any separated triple i<j<k such that either i∈ C, j∈ A,k∈ B or i∈ A,j∈ B,k∈ D will not be separated anymore after this merge. For each such triple, the potential function drops by w_ij+w_jk. Therefore:
	
    -Δ = (w_AB|C| + w_AC|B|)+ (w_AB|D|+ w_BD|A| )

		
To compare the two, we show that 	w_AB(|C|+|D|) ≥ w_AC|B| + w_BD|A|, and hence:

    Δ   = w_AB(|C|+|D|) ≥ w_AC|B| + w_BD|A|
       = -Δ-Δ ,

which implies <ref> as desired. To prove the last claim, note that average-linkage picks the pair (A,B) over both (A,C) and (B,D). Therefore, by definition of average-linkage: 
	
    w_AB|A||B|≥w_AC|A||C| w_AB|C| ≥ w_AC|B|

	
    w_AB|A||B|≥w_BD|B||D| w_AB|D| ≥ w_BD|A|

By summing  above inequalities, we get w_AB(|C|+|D|) ≥ w_AC|B| + w_BD|A|, which finishes the proof.


As a final note, in <Ref> we discuss hard instances for average-linkage under Gaussian kernels when d=1, where we essentially show the result of <Ref> is tight when comparing against , and there is no hope to get an approximation factor better than 3/4 for average-linkage in general for d=1. 



  
Random Cut. The following algorithm (termed as Random Cut) picks a uniformly random point r in the range [x_1, x_n] and divides the set of points into left and right using r as the splitter. The same process is applied recursively until the leaves are reached.





For d = 1 under any monotone distance-based similarity measure w_ij = g(x_i,x_j) the algorithm Random Cut obtains at least 1/2 fraction of the , and  hence gives a 1/2-approximation for the objective  in expectation.


For every triple i < j < k conditioned on partitioning the interval [i,k] for the first time the longer edge amongst (x_i, x_j) and (x_j, x_k) gets cut with probability p_1 ≥ 1/2 and the shorter with probability p_2 ≤ 1/2 so that p_1 + p_2 = 1.
W.l.o.g let's assume that (x_i, x_j) is the longer edge. 
Then the algorithm Random Cut scores w_jk p_1 + w_ij(1 - p_1) in expectation for the (i, j,k) triple. 
Note that: 

    w_jk p_1 + w_ij(1 - p_1) = (w_jk  - w_ij)(p_1 - 1/2) + 1/2(w_ij + w_jk)≥1/2(w_ij + w_jk)

By the linearity of expectation taking the sum over all triples i < j < k Random Cut scores at least 1/2 of  in expectation and hence gives 1/2-approximation for the objective .












§ HIERARCHICAL CLUSTERING IN HIGH DIMENSIONS

We now describe an algorithm Projected Random Cut which we use for high-dimensional data. This algorithm is given as Algorithm <ref>. It first projects on a random spherical Gaussian vector and then clusters the resulting projections using Random Cut.





For any input set of vectors v_1, …, v_n ∈ℝ^d the algorithm Projected Random Cut gives an α-approximation (in expectation) for the objective  under the Gaussian kernel similarity measure w_ij∼ e^-v_i - v_j_2^2/2σ^2 where α = (1 + δ)/3 for δ = min_i,jexp( - v_i - v_j_2^2/2σ^2).


Recall an upper bound on the optimum: 

    OPT ≤ = ∑_i < j < kmax(w_ij, w_ik, w_jk).

Fix any triple (i,j,k) where i < j < k. Note that the objective value achieved by the algorithm Projected Random Cut can also be expressed as ALG = ∑_i < j < k ALG_i,j,k where ALG_i,j,k is the contribution to the objective from the triple (i,j,k) defined as follows.
Consider the tree constructed by the algorithm. If v_k is the first vector in the triple (v_i, v_j, v_k) to be separated from the other two in the hierarchical partition (starting from the root) then A_i,j,k is defined to be w_ij (in the other two cases when i or j are separated first the definition is analogous).
Note that since Projected Random Cut is a randomized algorithm ALG_i,j,k is a random variable. By the linearity of expectation we have 𝔼[ALG] = ∑_i < j < k𝔼[ALG_i,j,k].
Thus in order to complete the proof it suffices to show that for every i < j < k it holds that:

    𝔼[ALG_i,j,k] ≥αmax(w_ij, w_ik, w_jk).



Fix any triple (v_i, v_j, v_k) which forms a triangle in ℝ^d. Conditioned on cutting this triangle for the first time let (p_ij,p_ik,p_jk) be the vector of probabilities corresponding to the events that the corresponding edge is not cut. I.e. this is the probability that we score the contribution of this edge in the objective. Note that p_ij + p_ik + p_jk = 1.

Consider any triangle whose vertices are v_i, v_j, v_k. To simplify presentation we set i = 1, j = 2, k = 3. We can assume that v_2 - v_1≥v_2 - v_3≥v_1 - v_3. Let θ_1 = v_1 - v_3 v_2 - v_3, θ_2 = v_2 - v_1v_3 - v_1 and θ_3 = v_1 - v_2v_3 - v_2 so that θ_1 ≥θ_2 ≥θ_3. See Figure <ref>. Note that the probability that the i-th longest side of the triangle has the longest projection is then θ_i/π. 



If (v_1, v_3) is the shortest edge in the triangle (v_1, v_2, v_3) then it holds that p_13≥1/3


Suppose that 𝐠 forms an angle π/2 - θ with (v_3 - v_1), i.e. the vector 𝐠^⊥ orthogonal to 𝐠 forms angle θ with v_3 - v_1. 
We define three auxiliary points x, y, z as follows (see also Figure <ref>.). Let ℓ_1 be a line parallel to 𝐠^⊥ going through v_3, let ℓ_2 to be a line parallel to 𝐠 going through v_2 and let ℓ_3 be the line parallel to ℓ_1 going through v_1.
We then let x be the intersection of ℓ_1 and (v_1, v_2), y be the intersection of ℓ_1 and ℓ_2 and z be the intersection of ℓ_2 and ℓ_3 (see Fig <ref>).

	
Thus the projections of v_3 - v_1, v_2 - v_1 and v_2 - v_3 on 𝐠 are y - z, v_2 - z and v_2 - y.
Hence conditioned on (v_1, v_2) having the longest projection the probability of scoring the contribution of (v_2, v_3) is given as p_23^12(θ) = y - z/y - v_2 since we are applying the Random Cut algorithm after the projection.
Note that by Thales's theorem we have p_23^12(θ) = y - z/y - v_2 = x - v_1/v_2 - v_1.
Applying the law of sines to the triangles (v_1, v_2, v_3) and (x, v_1, v_3) we have:

    sinθ_1/v_1 - v_2 = sin (θ_1 + θ_2)/v_1 - v_3,   sinθ/v_1 - v_3 = sin (θ + θ_2)/v_1 - x,

where we used the fact that sinθ_3 = sin (π - θ_1 - θ_2) = sin (θ_1 + θ_2). Similarly, we use the fact that sin (v_3 - x v_1 - x) = sin (θ + θ_2).
Using the above we can express p_23^12(θ) as:

    p_23^12(θ) = sinθ/sin (θ + θ_2)sin (θ_1 + θ_2)/sinθ_1


Thus the overall probability of scoring the contribution of the edge (v_1, v_3) conditioned on (v_1, v_2) having the longest projection which we denote as p_13^12 is given as:

    p_23^12   = 1/θ_1∫_0^θ_1 p_23^12 dθ


Similarly, consider the probability p_13^12 of scoring the contribution of (v_1,v_3) conditioned on (v_1,v_2) having the longest projection.
We can express it as:

    p_13^12   =  1/θ_1∫_0^θ_1 p_13^12(θ) dθ,

where 
    p_13^12(θ) = sinθ/sin (θ + θ_3)sin (θ_1 + θ_3)/sinθ_1.
 

Below we will show that p_13^12≥ p_23^12.
In fact, we will show that for any fixed θ∈ [0, θ_1] it holds that p_13^12(θ) ≥ p_23^12(θ).
Comparing the expressions for both it suffices to show that sin(θ_1 + θ_3)/sin (θ + θ_3)≥sin (θ_1 + θ_2)/sin (θ + θ_2) for all θ∈ [0, θ_1].
Since θ_1 = π - θ_2 - θ_3 this is equivalent to:

    sinθ_3/sin (θ + θ_2)≤sinθ_2/sin(θ + θ_3)

It suffices to show that:

    sinθ_3 sin(θ + θ_3) ≤sinθ_2 sin (θ + θ_2)

Using the formula sinαsinβ = 1/2(cos(α - β) - cos(α + β)) it suffices to show that:

    cos(θ + 2 θ_3) ≥cos(θ + 2 θ_2).

The above inequality follows for all θ∈ [0, π - θ_2 - θ_3] since θ_3 ≤θ_2 .
This shows that p_13^12≥ p_23^12. Since the probability that (v_1, v_2) has the longest projection is θ_1 / π we have that the probability of scoring (v_1, v_3) and (v_1, v_2) having the longest projection is at least θ_1/2π.
An analogous argument shows that the probability of scoring (v_1, v_3) and (v_2, v_3) having the longest projection is at least θ_2/2π.


Putting things together: 

    p_13≥1/2θ_1 + θ_2/π = 1/2π - θ_3/π≥1/3,

where we used that θ_3 ≤π/3 since θ_3 ≤θ_2 ≤θ_1.


We are now ready to complete the proof of Theorem <ref>.
Let γ = 2/3δ/(1 + δ) and note that γ≥δ / 3 since δ≤ 1.
If p_13≥ 1/3 + γ then the desired guarantee follows immediately.
Otherwise, if p_13≤ 1/3 + γ then we have:

    𝔼[ALG_1,2,3]/OPT_1,2,3   = p_13 w_13 + p_12 w_12 + p_23w_23/w_13
       ≥1/3 + (2/3 - γ)w_12/w_13≥1/3 + (2/3 - γ) δ
        = 1/3 + 2/3δ/1 + δ≥1 + δ/3,

where we used the fact that 

    w_12/w_13 = e^v_1 - v_3_2^2 - v_1 - v_2_2^2/2 σ^2≥ e^- v_1 - v_2_2^2/2 σ^2≥δ.





 §.§ Gaussian Kernels with small  

  <Ref> only provides an improved approximation guarantee for Projected Random Cut compared to the factor 1/3 (the tight approximation guarantees of average-linkage in high dimensions; see <Ref>) if δ is not too small, where δ =  min_i,jexp( - v_i - v_j_2^2/2σ^2). In particular, we get constant improvement if δ =Ω(1). Is this a reasonable assumption? Interestingly, we answer this question in the affirmative by showing that if we have δ = exp(- Ω̃(n)), then the Gaussian kernel can encode arbitrary similarity weights (up to scaling, which has no effects on multiplicative approximations). For simplicity, we only prove this result for {ϵ,1} weights here, while it can be generalized to arbitrary weights.

Given any undirected graph G=(V,E) on n nodes and ϵ>0, there exist unit vectors {k_v}_v∈ V in ℝ^d and bandwidth parameter σ∈ℝ_+, such that  d=O(n^2), 1/σ^2=Ω(nlog(1/)), and for some α>0 we have:

    ∀ (u,v)∈ E: e^-‖ k_u-k_v‖_2^2/2σ^2=α,
       ∀ (u,v)∉ E: e^-‖ k_u-k_v‖_2^2/2σ^2=αϵ.


 Our proof is constructive. Let d=n2. Pick orthogonal vectors {x_e}_e∈ E in ℝ^d such that ‖ x_e‖_2 =1/d_u d_v, where d_u is the degree of node u in graph G. For each v∈ V, define y_v∈ℝ^d as follows:

    y_v≜∑_(u,v) ∈ E(d_ud_v)x_uv.

Note that ‖ y_v ‖_2^2=∑_(u, v) ∈ Ed_u^2d_v^21/d_u^2 d_v^2=d_v
As the next step, pick a set of n orthonormal vectors {z_v} in the null space of {y_v}_v∈ V. Finally, for each v∈ V, define the final vector k_v∈ℝ^d as follows:

    k_v≜√(1-d_v/n)z_v+√(1/n)y_v

First, note that these vectors have unit length:

    ‖ k_v ‖_2^2=1-d_v/n+‖ y_v‖_2^2/n=1-d_v/n+d_v/n=1

Now, pick any two vertices u and v. If (u,v)∉ E, then:

    ⟨ k_u,  k_v ⟩   = ⟨√(1-d_v/n)z_u+√(1/n)y_u, √(1-d_v/n)z_v+√(1/n)y_v⟩
       =1/n⟨ y_u, y_v ⟩=0 ,

where we used the fact that ⟨ z_u, z_v⟩ = 0 , ⟨ z_u, y_v⟩ = ⟨ z_v, y_u⟩ =0 and the fact that ⟨ x_e,  x_e'⟩ =0 for every edge e incident to u and every edge e' incident to v, as e≠ e' when (u,v) ∉ E. Similarly, when (u,v)∈ E:

    ⟨ k_u, k_v ⟩ =1/n⟨ y_u, y_v⟩=1/n(d_u^2d_v^2‖ x_uv‖_2^2)=1/n

Now, consider a Gaussian kernel with bandwidth σ=(nlog(1/ϵ))^-1/2 and vectors {k_v}_v∈ V. Since k_u - k_v_2^2 = 2(1 - ⟨ k_u, k_v ⟩) from the above calculations of the inner products it follows that:

    ∀ (u,v)∈ E: e^-‖ k_u-k_v‖_2^2/2σ^2=e^1-1/n/σ^2,
       ∀ (u,v)∉ E:e^-‖ k_u-k_v‖_2^2/2σ^2=e^1/σ^2.

Thus the ratio is e^-1/nσ^2=ϵ, as desired.




§ HARD INSTANCES FOR AVERAGE-LINKAGE UNDER GAUSSIAN KERNEL SIMILARITY






  
High-dimensional case.

We embed the construction of <cit.> shown in Figure <ref> into vectors with similarities given by the Gaussian kernel.



There exists a set of vectors v_1, …, v_n ∈ℝ^d for d = poly(log n) for which the average-linkage clustering algorithm achieves an approximation at most 1/3 + o(1) for  under the Gaussian kernel similarity measure.


We start by this theorem:


For any constant ϵ∈ (0,1) the instance ℐ average-linkage clustering achieves the value of  at most 1/6 n^8/3 + O(n^7/3) while the optimum is at least 1/2 n^8/3 - O(n^7/3).




 Let Δ, τ > 0 be real-valued parameters to be chosen later.
We use indices i and j to index our set of vectors.
For i ∈{1, 2, …, n^1/3}, j ∈{1, 2, …, n^2/3} let

    v_i,j = Δ (e_i + (1 + τ)e_k + j),
 
where k = n^1/3 and e_i is the t-th standard unit vector e_t = (0 …, 1, …, 0) with the 1 in the t-th entry.
Then it is easy to see that for any fixed i ∈ [n^1/3] and j_1 ≠ j_2 ∈ [n^2/3] it holds that:

    v_i,j_1 - v_i, j_2_2^2 = 2 (1 + τ)^2 Δ^2

Similarly, for any fixed j ∈ [n^2/3] and i_1 ≠ i_2 ∈ [n^1/3] it holds that:

    v_i_1, j - v_i_2,j_2^2 = 2 Δ^2.

Otherwise if i_1 ≠ i_2 ∈ [n^1/3] and j_1 ≠ j_2  ∈ [n^2/3] then:

    v_i_1, j_1 - v_i_2, j_2_2^2 = 2 Δ^2 + 2 (1 + τ)^2 Δ^2 ≥ 4 Δ^2.


By setting Δ^2 = 2 σ^2 c log n for a sufficiently large constant c the contribution of pairs of vectors with i_1 ≠ i_2 and j_1 ≠ j_2 can be made negligible. Let 2 (1 + τ)^2 Δ^2 = α and 2 Δ^2 = β.
The rest of the pairs correspond to an hard instance ℐ for which average-linkage only achieves a 1/3 + o(1)-approximation compared to the optimum.
By setting τ = 1 / poly(log(n)) we have e^(β^2 - α^2)/2σ^2  = Ω(1) and hence by Theorem <ref> it follows that average-linkage clustering can't achieve better than 1/3 + o(1) approximation for this instance.

Finally, note that by applying the Johnson-Lindenstrauss transform we can reduce the dimension required for the above reduction to d = poly(log n). Indeed, projecting on a random subspace of dimension O(log n/ξ^2) would preserve ℓ_2^2-distances between all pairs of vectors up to a multiplicative factor of (1 ±ξ) setting ξ = τ/10 = 1/poly(log n) it follows that our hard instance can be embedded in dimension d = poly(log n).


































  
Low-dimensional case.

For d = 1 a hard instance for average-linkage clustering can also be constructed.


For points x_1, …, x_n ∈ℝ average-linkage clustering achieves approximation at most 3/4 for  under the Gaussian kernel similarity measure.


Consider the instance consisting of four equally spaced points on a line, i.e. 0, Δ, 2Δ, 3 Δ.
Then (after shifting the two middle points slightly) the average-linkage clustering algorithm might first connect the two middle points and then connect the two other points in arbitrary order. We denote the cost of this solution as AVG.
An alternative solution would be to create two groups (0, Δ) and (2Δ, 3Δ) and then merge them together.
We denote the cost of this solution as OPT.
By making Δ sufficiently large the contribution of pairs at distance more than Δ from each other can be ignored. We thus have  AVG ≤ (√(2 π)σ)^-n (3 e^-Δ^2/2σ^2 + o(e^-Δ^2/2σ^2)) and OPT ≥  (√(2 π)σ)^-n  4 e^- Δ^2/2σ^2, which gives the ratio of 3/4 for sufficiently large Δ.


In the four-point instance of the proof of <Ref>, the  evaluates to 

(√(2 π)σ)^-n 6 e^- Δ^2/2σ^2, and hence average-linkage cannot obtain more than 1/2 of .








§ EXPERIMENTAL RESULTS



In this section we demonstrate the quality of the solution returned by Projected Random Cut (PRC) on a small dataset, and highlight that running time only scales linearly on a large dataset.
PRC does not compute the similarity weights (it only takes one pass over the feature vectors with dimension-free memory requirements) and then sorts the projected points in time O(nlog n). Hence, it is fast even in use-cases with over millions of datapoints (in contrast to average-linkage, spectral clustering, or even single-linkage which all have superlinear running times in high dimensions).

We run PRC algorithm on two real datasets from the UCI ML repository <cit.>: (1) the Zoo dataset <cit.> (the small dataset) that contains 100 animals given as 16D feature vectors (this dataset comes from applications in biology), and (2) the SIFT10M dataset <cit.> (the large dataset) that contains around 10M datapoints, where each datapoint is a 128D Scale Invariant Feature Transform (SIFT) vector (this dataset comes from applications in computer vision).  For more details, refer to the Appendix <ref>.




  
Small dataset (Zoo). We compare PRC to (1) the recursive spectral clustering using the second eigenvector of the normalized Laplacian of the weight matrix (Spectral) <cit.>, (2) average-linkage (AL), and (3) , an upper-bound on the objective value of the optimum tree; see <Ref>. In contrast to PRC, these three benchmarks need to compute the weights and are slow on large datasets. 

 <Ref> summarizes the results of this experiment for various choices of the parameter σ (first column). The second, third and fourth columns report the objective values () of the trees returned by PRC, Spectral and AL respectively, and the fifth column gives an empirically observed approximation guarantee for PRC by comparing it against the upper bound on optimum . 
 We observe that PRC attains high approximation factors of ≈ [0.73,…,0.92] compared to  . Moreover, PRC obtains competitive objective values compared to Spectral and AL, although it runs faster (in almost linear-time). 
 



On the choice of bandwidth parameter σ, as a rule of thumb, we find an interval [α,β] such that (1) if σ>β, a considerable portion of pairs of datapoints with very different weights under the cosine similarities <cit.> have almost equal weights under the Gaussian kernel, and (2) if σ<α,  a considerable portion of pairs of datapoints with almost equal weights under the cosine similarities have very different weights under the Gaussian kernel. We end up with α=1.5,β=5. 
 




  
Large dataset (SIFT10M) The focus of this experiment is measuring the running time of PRC, and showing that it only scales linearly with the dataset size. Note that evaluating the performance of any other algorithm or upper bound (or even one pass over the similarity matrix) would be prohibitive.

We run PRC on truncated versions of SIFT10M of sizes 10K, 100K, 500K, 1M and 10M. σ is set to 450 [Note that the σ value does not affect the running time.]. We emphasize that our PRC algorithm runs extremely fast on a 2014 MacBook[System specs: 8 GB 1600 MHz DDR3 RAM, 2,5 GHz Intel Core i5 CPU.]. The running times are summarized in <Ref>. Observe that PRC scales almost linearly with the data and has almost the same running time as just a single pass over the datapoints.























plainnat
















§ GREEDY CUTTING AND SINGLE-LINKAGE.

Consider a simple algorithm, denoted by Greedy Cut, that picks the interval  with maximum length among {[x_i,x_i+1]}_i=1:n-1 (lets say [x_m,x_m+1]), and repeats the same operation recursively on (x_1,…,x_m) and (x_m+1,…,x_n) until the leaves are reached. 

For d=1 and under any monotone distance-based similarity measure  Greedy Cut and single-linkage return the same HC tree. Moreover, the edges picked by Greedy Cut are exactly the same edges picked by single-linkage, picked in reverse order. 


It is known that single-linkage is essentially the Kruskal algorithm (and hence edges picked by single-linkage form a Maximum Spanning Tree (MST)). Clearly, for any monotone distance-based measure the line connecting x_1 to x_n is the unique MST (as any tree can be shortcut-ed with this line), and hence single-linkage picks intervals {[x_i,x_i+1]}_i=1:n-1 in increasing order of their lengths. At the same time, Greedy Cut picks also the same intervals, but in decreasing order of their length. Moreover, because single-linkage merges the edges picked by Greedy Cut in the reverse ordering (and it creates the HC tree from bottom to top), it returns the same HC tree as Greedy Cut.



 As a simple observation, Greedy Cut is equivalent to reverse-Kruskal in 1D; It starts from all edges,  goes over them in increasing order of weights (here, in decreasing order of lengths), and only keeps an edge when its removal makes the graph disconnected. We claim the first edge picked by reverse-kurskal corresponds to the interval [x_m,x_m+1] with the maximum length, and hence the equivalence between the two algorithms by induction. To prove the claim, the line between x_1 and x_n keeps the graph connected, so the first picked edge corresponds to an interval (x_i,x_i+1). Also, it should be the interval [x_m,x_m+1] with maximum length. Moreover, removing this edge makes the graph disconnected, because otherwise there is was cross edge between the left-side (x_1,...,x_m) and the right-side (x_m+1,…,x_n). The length of such an edge is more than the length of (x_m,x_m+1), so it should have been removed before, a contradiction.


We finish the section by demonstrating the lack of performance of single-linkage (and hence Greedy Cut) through an example, which justifies using both average-linkage and Random Cut as smooth versions of single-linkage for 1D data points.
 For d=1, single-linkage can obtain at most 12 of the optimum objective, and 12 of the objective values of average-linkage and Random Cut, under Gaussian kernels.


Consider an example with n equally spaced points on a line, where the distance between any two adjacent node is Δ. We now slightly move the points so that 
weight of (x_i,x_i+1) is (√(2 π)σ)^-n (1-(i-1)ϵ) e^- Δ^2/2σ^2 for i=1:n-1. Now, single-linkage peels off points in the order x_1,…,x_n. Roughly speaking, we let Δ/σ^2 to be large enough so we can ignore the similarity weights between any two non-adjacent points. Hence single-linkage gets an objective value of ≈(∑_i=1^n-1 (n-i))(√(2 π)σ)^-n e^- Δ^2/2σ^2, which evaluates to (n^2/2+o(n^2))(√(2 π)σ)^-n e^- Δ^2/2σ^2. Average-linkage returns the symmetric binary tree (and hence the objective value is 
    ((n-2)(n-1)-nlog(n))(√(2 π)σ)^-n e^- Δ^2/2σ^2=(n^2+o(n^2))(√(2 π)σ)^-n e^- Δ^2/2σ^2
. Random returns a random binary tree, with (roughly speaking) a similar objective value of 
    (n^2+o(n^2))(√(2 π)σ)^-n e^- Δ^2/2σ^2
 in expectation. The fact that optimum objective value is as large as these two quantities finishes the proof.






























































§ DEFERRED DISCUSSIONS IN 




  
Datasets.
We use real data in our experiments. To be able to measure the performance of produced Hierarchical Clusterings, we need to compute , and hence we use a small dataset. To be able to show linear scaling of running time even for tens of millions of datapoints, we will use large datasets of sizes 10K to 10M.


The small dataset: The Zoo dataset contains 100 animals given as 16-dimensional vectors, forming 7 different classes (mammals, amphibians, etc.). The features contain information about the animal's characteristics (if it has tail, hair, fins, etc.). Here, we want to showcase the quality of the solution produced by PRC; in the case of the small dataset, we can afford keeping track of two benchmarks: the performance of the widely-used Spectral algorithm and the [ is based on triples of datapoints, thus taking O(n^3) time.]. 

The large dataset: In the SIFT10M dataset, each data point is a SIFT feature, extracted from Caltech-256 <cit.> by the open source VLFeat library <cit.>.
Caltech-256 is used as a computer vision benchmark image data set, that features a total of 256 different classes with high intra-class variations for each category. Each datapoint is a 128-dimensional vector and similar to the Zoo dataset, we use this information as features to perform fast Hierarchical Clustering. Here, we want to test the scalability of our algorithm so we use successively larger datasets from the SIFT10M dataset of 10K, 100K, 500K, 1M and finally 10M datapoints by truncating the data file as necessary.


