Hierarchical clustering with dot products recovers hidden tree structure
Authors: Annie Gray
Alexander Modell
Patrick Rubin-Delanchy
Nick Whiteley
Venue: NeurIPS 2023 spotlight (2023)
Keywords: agglomerative clustering
generative model
graphical model
hierarchical clustering
high-dimensional data
TL;DR: A new perspective on recovery of tree structure from data.

Abstract: In this paper we offer a new perspective on the well established agglomerative clustering algorithm, focusing on recovery of hierarchical structure. We recommend a simple variant of the standard algorithm, in which clusters are merged by maximum average dot product and not, for example, by minimum distance or within-cluster variance. We demonstrate that the tree output by this algorithm provides a bona fide estimate of generative hierarchical structure in data, under a generic probabilistic graphical model. The key technical innovations are to understand how hierarchical information in this model translates into tree geometry which can be recovered from data, and to characterise the benefits of simultaneously growing sample size and data dimension. We demonstrate superior tree recovery performance with real data over existing approaches such as UPGMA, Ward's method, and HDBSCAN.

C, in which clusters of data points are merged Hierarchical clustering with dot products recovers hidden tree structure
Annie Gray
Alexander Modell Department of Mathematics Imperial College London, UK a. modeeollll fodeel@@ing 1elmmmmmiiaal.... ac uk
School of Mathematics
University Bristol, UK
annie. gray@bristol.ae uk
Patrick Rubin-Delanchy
Nick Whiteley School of Mathematies University of Bristol, UK
School Mathematics University of Bristol, UK
patrict.rutn..elnncyy0rrr...... uk
nick. whiteeeyybbrittl ac. uk
Abstract
In this paper offer perspective the well established agglomerative clus- tering algorithm, focusing on recovery of hierarchical structure. We recommend a simple variant of the standard algorithm, in which clusters merged by maximum average dot product and not, for example, by minimum distance or within-cluster variance. We demonstrate that the output by this algorithm provides bona fide estimate of generative hierarchical structure in data, under a generic probabilistic graphical model. The key technical innovations are to understand how hierarchical information in this model translates into tree geometry which can be recovered from data, and to characterise the benefits of simultancously growing sample size and data dimension. We demonstrate superior tree recovery performance with real data over existing approaches such as UPGMA, Ward's method, and HDBSCAN,
Introduction
Hierarchical structure is known occur in many natural and man-made systems 32 and the problem considered in this paper how to recover this structure from data. Hierarchical clustering algorithms, are very popular techniques which organise data into nested clusters, used routinely by data scientists and machine learning researchers, and are easily accessible through open source software packages such scikit-learn [38). We focus on perhaps the most popular family of such techniques: agglomerative clustering recursively.
[25
36 37 41
Agglomerative clustering methods are not model-based procedures, but rather simple algorithms. Nevertheless, in this work we uncover new perspective on agglomerative clustering by introducing a general form generative statistical model for the data, but without assuming specific parametric families of distributions (e.g., Gaussian). In our model (section 2.1). hierarchy takes the form of a tree defining the conditional independence structure of latent variables, using elementary concepts from probabilistic graphical modelling [30 28) In key innovation,
then augment this conditional independence tree to form what we call a dendrogram, whose geometry is related to population statistics of the data. The new insight which enables tree recovery in setting (made precise and explained section 3 that dot products between data vectors reveal heights of most recent common ancestors in the dendrogram.
We suggest an agglomerative algorithm which merges clusters according to highest sample average dot product (section 2.2p. This is in contrast to many existing approaches which quantify dissimilarity
37th Conference on Neural Information Processing Systems (NcurIPS 2023).

---

between data vectors using Euclidean distance. We also consider the case where data are preprocessed by reducing dimension using PCA. We mathematically analyse the performance of our dot product clustering algorithm and establish that under our model, with sample size and data dimension growing simultaneously at appropriate rates, the merge distortion 20 between the algorithm output and underlying tree vanishes. In numerical examples with real data (section + where we have access to labels providing a notion of true hierarchy, we compare performance of our algorithm against existing methods in terms of a Kendall correlation performance measure, which quantifies association between ground-truth and estimated tree structure. We examine statistical performance with and without dimension reduction by PCA, and illustrate how dot products versus Euclidean distances relate to semantic structure in ground-trum hierarchy.
Related work Agglomerative clustering methods combine some dissimilarity measure with 'linkage' function, determining which clusters are combined. Popular special cases include UPGMA (4] and Ward's method [48), against which we make numerical comparisons (section 4. Owing to the simple observation that, in general, finding the largest dot product between data vectors is not equivalent finding the smallest Euclidean distance, not correctly recover tree structure under our model (appendix E). Popular density- based clustering methods methods include CURE [22], OPTICS (7], BIRCH [52] and HDBSCAN [9]. In section 4 we discuss the extent to which these methods can and cannot be compared against ours in terms of tree recovery performance.
explain why these existing methods may
Existing theoretical treatments of hierarchical clustering involve different mathematical problem formulations and assumptions to ours. One common setup is assume underlying ultrametric space whose geometry specifies the unknown tree, and/or to study tree recovery as n + oo with An alternative problem formulation addresses recovery of the cluster tree of the probability density from which it is assumed data sampled 42 20 27]. The unknown tree in our problem formulation specifies conditional independence structure, and so it has a very different interpretation to the trees in all these cited works. Moreover, our data are vectors in I.," and p - oo is crucial aspect of our convergence arguments, but in the above cited works plays no role or is fixed. Our definition of dendrogram is different to that in e.g. [10] we do not require all leaf vertices to be equidistant from the root condition which arises as the "fixed molecular clock" hypothesis [8 34 in phylogenetics. We also allow data be sampled from non-leaf vertices. There is an enormous body of work on tree reconstruction methods in phylogenetics, e.g. listed at [4), but these are mostly not general-purposs solutions to the problem of inferring hierarchy. Associated theoretical convergence results are limited in scope, e.g, the famous work [88] is limited to a fixed molecular clock, five taxa and does not allow observation error.
respect to cost function, e.g, [TO 16
43
15 31 12]
2 Model and algorithm
2.1 Statistical model, tree and dendrogram
Where possible, we use conventional terminology from the field of probabilistic graphical models, e.g.. 28 to define our objects and concepts. Our model is built around an unobserved tree T - (V,E), that is directed acyelic graph with vertex and edge
V and &, with two properties: T is connected (ignoring directions of edges), and each vertex has at most one parent, where we say u is parent of there an edge from u to v. We observe data vectors Y, RP. 1,..., .n.,..,,, n, which we model as:
Y, = x + S(Z,),
(I)
comprising three independent sources of randomness:
Z...... Zn are i.i.d., discrete random variables, with distribution supported on subset of vertices Z C V, |Z| < oc;
X(v) = (X, ... (v)) is an IVPuud random vector for each vertex v E V; Ei...., E, are i.idd R valued random vectors. The elements of E, are i.i.d., zero mean and unit variance. For each E Z. S(z) RPx,, is a deterministic matrix.
For each u E Z, one can think of the vector X(w) as the random centre of "cluster", with correlation structure of the cluster determined by the matrix S()). The latent variable Z, indicates which cluster

---

the ith data vector Y; is associated with. The vectors X(v), v E V\ Z, correspond to unobserved vertices the underlying tree. As an example, Z could be the set of leaf vertices of T, but neither our methods nor theory require this to be the case. Throughout this paper, we assume that the tree TT determines two distributional properties of X. Firstly, we assume T is the conditional independence graph of the collection of random variables X, == (x) E V} for each j. that is, A1. for all j = .......P P, the marginal probability density or mass function of X, factorises as:
p(zz)=
II p
PEV
where Pa., denotes the parent of vertex
However we do not necessarily require that XI...., Xp are independent identically distributed. Secondly, we assume
.. A2. for each j = ......p. the following martingale-like property holds:
(PX, = (Pa.);
for all vertices u V except the root.
The conditions A and Anjuuee an additive hierarchical structure in X. For any distinct vertices u, E V and w an ancestor of both, A and Amimly that given (y(w), the increments X,()) (w) and (v) - (w) are conditionally independent and both conditionally mean zero.
To explain our algorithm we need to introduce the definition of a dendrogram, D = (T, (T. where h: V- 1 is function which assigns height each vertex of T, such that h(v) h(Pa,) for any vertex v E V other than the root. The term "dendrogram" is derived from the ancient Greek for "tree" and "drawing", and indeed the numerical values h(v), a drawing of T where height is measured with respect to some arbitrary baseline on the page, an example shown figure aa). With (:--) denoting the usual dot product between vectors, the function
V, can be used to construct
au, v) 1 (x) X())),, U,
(2)
will act as measure of affinity underlying our algorithm. The specific height function we consider is hvv)= == a(v, v)). The martingale property ensures that this height function satisfies h(v) 2 (Paa.) as required, see lemma 2 in appendix C In appendix EE] we also consider cosine similarity as an affinity measure, and analyse its performance under a multiplicative noise model, cf. the additive model
2.2 Algorithm
= = ,, Z2. for i [n]. Combining the model 1 with 2 we have a Zy)
[n] = {I, (I....,n], n... The imput to algorithm [ is an estimate .(.,)) of all the pairwise affinities aZZ, Zz), #jj E [n]. We consider two approaches to estimating 0r 74, Zj). with and without dimension reduction by uncentered PCA. For some r < min(p, let V E RRxxr denote the matrix whose columns are orthonormal eigenvectors of Σi= Y associated with its r largest eigenvalues. Then & = VT,, is r-dimensional vector of principal component scores for the ith data vector. The two possibilities for a(.,)) -) we consider are:
dddnaa (i. = .Y Y,,,
(Go, C3)-
(3)
In the case of dpc the dimension must be chosen. Our theory in section 3.2 assumes that v is chosen as the rank of the matrix with entries a (u, v), u, 0 E Z, which is at most [Z]. In practice r usually must be chosen based on the data, we discuss this in appendix[A]
Algorithm I returns& a dendrogram D (7, h),, comprising tree + (V,E) and height function h. Each vertex in V is a subset of [n). thus indexing subset of the data vectors Y YI...., Y,. The leaf vertices are the singleton sets (i), i E [n]). corresponding to the data vectors themselves. As algorithm I proceeds, vertices are appended to 12, edges are appended to É, and the domain of the function a(.,)) is extended as affinities between elements of 1/ are computed. Throughout the paper we simplify notation by writing j(i.j) as shorthand for a((i), (i)) for i, E [m). noting that each argument of a(,,)) -) is in fact a subset of [n).
3

---

(b)
(4)
X3
0
Y6
(5)
(6)
d
h(e) mla.b)
21 2z 2z b, Ze =
1 2(112,3,44) E.-12 a(3,4)
3 df
height of Y. = ||x1111 21, ......... 2ol
h((6)) a(6,6) 311661 -
Figure 1: (a) An example of the dendrogram D with V d, f.9), see lemma for interpretation of the merge height mn(--) and distance d(:,-). Horizontal distances in this diagram chosen arbitrarily. (b) DD augmented according to = {a, b, d} and the realization: 21, Zz. Z3 a. Z4 = b, 75 Ze = d, see section 3.3 for discussion. (c) The dendrogram D output from algorithm
in the case â Ódata D can be seen approximate the dendrogram in (b) and hence D.
Algorithm Dot product hierarchical clustering Input: pairwise affinities -(.,)) between n data points
1: Initialise partition Po {{1)...., {n}). vertex set V = Po and edge set to the empty set.
2: for m E f1..... n 1} do
Pm-1 with largest affinity 4: Update Pm-1 to Pm. by merging u, v to form w == U
3: Find distinct pair
5: Append vertex w. to 1 and directed edges w - u and w to
6: Define affinity between w and other members of Pm as
|u] W
|0
a(w,,)
a(u, :)
a(u,)) ))
and a(w, w) = a(u, v).
7: end for
Output: Dendrogram D = ( h)), comprising tree + (V,E) and heights h(v) = á(u, v) for E V Po. and h(v) max(h((a.)) o(v for v E Po-
Implementation using scikit-learn The presentation of algorithm 1 has been chosen to simplify its theoretical analysis, but alternative formulations the same method may much more computa- tionally efficient in practice. In appendix B we outline how algorithm [ can easily be implemented using the AggloomrrativeCluseriii class in scikit-learn (38]
Performance Analysis
3.1 Merge distortion upper bounded by affinity estimation error
In order to explain the performance of algorithm T we introduce the merge height functions: m(u, v) == fomost recent common ancestor of u and v).
u, v E V,
m (,, v) (most recent common ancestor of and v).
E 1.
To simplify notation we write ni(i..) .) as shorthand for mi(i), (i}). for i, E [n]]. The discrepancy between any two dendrograms whose vertices are in correspondence can quantified by merge distortion 20 the maximum absolute difference in merge height across all corresponding pairs of

---

vertices. (00 advocated merge distortion as a performance measure for cluster-tree recovery, which different to our model-based formulation, but merge distortion turns out be useful and tractable performance measure in our setting too. As preface to our main theoretical results, lemma lemma[]explainn exxlains how the geometry of the dendrogram D, in terms of merge heights, is related to population statistics of our model. Defining d(u, v) == hu)) h(w) + h(v) h(w) for * E V, where w is the most and v, we see du, v) is the vertical distance on the dendrogram from down then back up as illustrated in figure Ta).
recent common ancestor
Lemma 1. For any two vertices u, V V.
m(u, v) = 1 [X U), X (())) = a(u, d(u,v) E | ((
(4)
The proof is in appendix
C Considering the first two equalities in 44) it is natural to ask the estimated affinities a(--)) being close to the affinities a(:;)) implies small merge distortion
between ri (:))) and m(..). This the subject of our first main result, theorem below. In appendix E use the third equality 44) explain why popular agglomerative techniques such UPGMA 44 and Ward's method [8 which merge clusters based on proximity in Fuclidean distance may
enjoy limited success under our model, but in general do not correctly recover tree structure.
Let b denote the minimum branch length D, that is, = min(h(v) - h(Pa.)), where the minimum taken over all vertices in V except the root.
Theorem 1. Let the function á(:;") :) given as input to algorithm | real-valued and symmetric but otherwise arbitrary. For any , E 2, if
max [or(zs; mj) a (i. < b/2, k.E[m)),
then the dendrogram returned by algorithm satisfies
max
mm22, 24) mii,))) < max |a(2;; (((=,,)))a((1)))) A.j|m|lis)
(5)
(.jE[m],17]
The proof is in appendix
3.2
A Afiinity estimation error vanishes with increasing dimension and sample size
Our second main result, theorem 2 below, concerns the accuracy of estimating the affinitics a(::-)
using ådata or ácaaa as defined in B). We shall consider the following technical assumptions.
A3 (Mixing across dimensions). For mixing coefficients yp satisfying Σx²1 (13())) < oo and all
u. E Z, the sequence (u), (5))); 2 1} is 4p-mixing.
A4 (Bounded moments). For some q 2 supj21 maX [1]x)( Do and LE| E11 22] co. where En1 is the first element of the vector E..
A5 (Disturbance control). max,ez [[S(c)]lop O(1) p - where || llop is the spectral norm.
A6 (PCA rank). The dimension r chosen in definition of åpens see 3 equal to the rank of the matrix with entries a(u, v), u, Z.
The concept of çp-mixing is a classical weak deak-dependence condition, e.g. 19 39 implies that for each j 2 1, (X, (())) and (X)+o(u) (++5())). are asymptotically independent as + co. However, it is important to note that Gkaata and áipca in (3), and hence the operation of algorithm T with these inputs, are invariant to permutation of the data dimensions j = 1...., p. Thus our analysis under ABJoll requires there is some permutation of dimensions under which m-mmxxig holds. A a fairly mild integrability condition. Allallows control of magnitudes of the "disturbance" vectors Y, (2+)= 5(Z.),, Further background and discussion of assumptions is given in appendix
Theorem 2. Assume that A A5 hold and let be in A4 Then
2/4
max _
Jar Jar (i, 2z) Eaatta (i, E Op /p
(6)
If additionally A3 is strengthened from ip-mixing independence, a 2 and all v Z (in which case A5 holds), and Al6 holds, then
S(v) oI, for some constant
(z, Z2) 6pea (i, O (VTL ni {.E(m)),
max
+
5

---

The proof of theorem 2 is in appendix To prove we use recent cenfform-concentratinn result for principal component scores from [9]. Overall, theorem 2 says that affinity estimation error vanishes if the dimension p grows faster enough relative to n (and in the case of Öpea. noting that under A6 r I [2], so sensible to assume is much smaller than and p). The argument of Oz(-) is the convergence rate; if (Xp.)) is some collection of random variables indexed by p and n, Xp,n E (z/nn/1//1 means that for any e 0, there exists and M such that / 2 M implies P(I.Xp,.l > 6)
We give an original and self-contained proof of (6).
We note the convergence rate in (6) is decreasing in where as rate in 7 is not. It is an open mathematical question whether 7 can be improved in this regard; sharpening the results of Whiteley et al. [49 used
in the proof of 7 seems very challenging. However, when q = 2, ( gives Op (n/ vp) compared
to Og (V vr/n) in i.e. improvement from to /ar in the first term. We explore
empirical performance of ddata versus ppea in section 4.2
3.3 Interpretation
By combining theorems and 2 find that when b 0 is constant and & is either ddata dipca the merge distortion
max Im Z, 2z) fi((i (
(8)
converges to zero at rates given by the r.h.s of (6) and To gain intuition into what (8) tells
us about the resemblance between D and D, it is useful to consider intermediate dendrogram illustrated in figure [[)) which conveys the realized values of Zi....,Zy, Zn This dendrogram is constructed from 2> by adding leaf vertex corresponding to each observation Y, with parent Z; and height p -' Y | Zn] = h(Z,) + P 't (r(SS(2) (2,) TS(Z.)] and deleting any o E Z such that Z F for all E [ni] (e.g.. vertex in figure (b). The resulting merge height between the vertices corresponding to Y, and Y, is m Zz, Zz)- & being small implies this must be close to j) in D as in figure (c). Moreover, in the case á Odtam the height h(()) of leaf vertex {i} in D is upper bounded by 1 which under our statistical assumptions is concentrated around p- D-'Ei] Y
1.1. i... the heights of the leaves in figure (c) approximate those of the corresponding leaves in figure [b).
Overall we see that D in figure
(c) approximates the dendrogram in figure (b. and in turn D. However even if m Zz, Zn) ri for all i,j, the tree output from algorithm T, may not
be isomorphic (i.e., equivalent up to relabelling of vertices) to the tree in figure b): 7 is always binary and has 2n- vertices, whereas the tree in figure bb) may not be binary, depending on the underlying T and the realization of Zi..... Zn This reflects the fact that merge distortion, in general, is pseudometric of true dendrograms D, for instance binary trees with non-zero branch lengths, then asymptotically
dendrograms. However, if one restricts attention to specific classes
algorithm I can recover them exactly. We explain this point further in appendix
4 Numerical experiments
We explore the numerical performance of algorithm I in the setting of five data sets summarised below. The real datasets used are open source, and full details of data preparation and sources are given in appendix D
Simulated data. A simple tree structure with vertices V = 2,,,,3,455, 7,8). edge set E {6-1, 6-2,66666 6-3, 7-4,7 7-5, 8-6, 8-7) and Z = 3, (the Icaf vertices). Zi...., Zn are drawn from the uniform distribution on Z. The Xy(v) are Gaussian random variables, independent across j. Full details of how these variables are sampled are in appendix D The elements of E, are standard Gaussian, and S(r) olp with o =
20 Newsgroups. We used random subsample of 5000 documents from the well-known 20 Newsgroups data set (29]. Each data vector corresponds to one document, capturing its 12818 Term Frequency Inverse Document Frequency features. The value of n was chosen to put us in the regime p 2 to which our theory is relevant see section 3.2 Some ground-truth labelling of documents is known: each document is associated with of 20 newsgroup topics, organized two hierarchical levels.

---

Zebrafish gene counts. These data comprise gene counts in zebrafish embryo cells taken from their first day of development [47]. As embryos develop, cells differentiate into various types with specialised, distinct functions, so the data are expected to exhibit tree-like structure mapping these changes. We used a subsample such that n = 5079 and p = 5498 to put us in the p 2 n regime. Each cell has two labels: the tissue that the cell is from and a subcategory of this.
Amazon reviews. This dataset contains customer reviews on Amazon products M. A random sample of n = 5000 is taken and each data vector corresponds to one review with = 5594 Term Frequency Inverse Document Frequency features. Each product reviewed has labels which make up a three-level hierarchy of product types.
S&P 500 stock returns. The data are p 1259 daily returns between for n 368 stocks which were constituents of the S&P 500 market index (2) between 2013 to 2018. The two-level hierarchy of
stock sectors by industries and sub-industries follows the Global Industry Classification Standard (3.
4.1 Comparing algorithm to existing methods
We numerically compare algorithm I Jggainst three very popular variants of agglomerative clustering: UPGMA with Euclidean distance, Ward's method, and UPGMA with cosine distance. These are natural comparators because they work by iteratively merging clusters in manner similar to algorithm T. but using different criteria for choosing which clusters to merge. In appendix [+] we complement our numerical results with mathematical insights into how these methods perform under our modelling assumptions. Numerical results for other linkage functions and distances are given in appendix D Several popular density-based clustering methods use some hierarchical structure, such as CURE (22), OPTICS 77 and BIRCH (52) but these have limitations which prevent direct comparisons: they aren't equipped with a way to simplify the structure into a tree, which it is our aim to recover, and only suggest extracting a flat partition based on a density threshold. HDBSCAN [9] is density-based method that doesn't 't have these limitations, and report numerical comparisons against it.
Kendall Tr ranking correlation. For real data some ground-truth hierarchical labelling may be available but ground-truth merge heights usually are not. We need a performance measure to quantitatively compare methods operating on such data. Commonly used clustering performance measures such as the Rand index [40) and others [23 21] allow pairwise comparisons between partitions, but do not capture information about hierarchical structure. The cophenetic correlation coefficient [45] is commonly used to compare dendrograms, but relies on an assumption that points close in Euclidean distance should be considered similar which is incompatible with our notion of dot product affinity. To overcome these obstacles we formulate performance measure as follows. For each of n data points, we rank the other n data points according to the order in which they merge with it in the ground-truth hierarchy. We then compare these ground truth rankings to those obtained from a given hierarchical clustering algorithm using the Kendall 7. correlation coefficient (26). This outputs value in the interval [-1, 1]. with 11, 1 and 0 corresponding to negative, positive and lack of association between the ground-trutt -ruth and algorithm-derive association value across all n data points as the overall performance measure. Table I shows results with raw data vectors 1mm or PC scores Çim taken as input to the various algorithms. For all the data sets except S&P 500, algorithm 1 iis found to recover hierarchy more accurately than other methods. We include the results for the S&P 500 data to give a balanced scientific view, and in appendix [E] we discuss why our modelling assumptions may not be appropriate for these data, thus explaining the limitations of algorithm I
rankings. We report the mean
4.2 Simulation study of dot product estimation with and without PCA dimension reduction
For high-dimensional data, reducing dimension with PCA prior to clustering may reduce overall computational cost. Assuming Ç1.. are obtained from, e.g., a partial SVD, in time O(ppr), the time complexity of evaluating åepc is Onpr + n²r), versus O(npp) for ådma- although this ignores the cost of choosing r. In table [[] we see for algorithm [ the results for input Y are very similar to those for Çi.-- To examine this more closely and connect our findings to theorem [ we now compare ddata and ápen as estimates of a through simulation. The model is as described at the start of section 4 In figure (b))-bb, we see that when p is growing with n, and when p is constant, the åpea error is very slightly smaller than the d error. By contrast, in figure 2 c), when 10 is
7

---

Table 1: Kendall Tr, ranking performance measure. For the dot product method, i.e. algorithm I 1 as input corresponds using Odma. and (1:m corresponds ápea. The mean Kendall T. correlation coefficient reported alongside the standard error (numerical value shown is the standard errorx 10%).
Data
Input
Dot product
UPGMA w/ cos. dist.
HDBSCAN
UPGMA w/ Eucl. dist.
Ward
Y 1mm
0.26 (2.9)
0.26 (2.9)
-0.010 (0.65)
0.23 (2.7)
0.18 (2.5)
Newsgroups
Cim
0.24 (2.6)
0.18 (1.9)
-0.016 (1.9)
0.038 (1.5)
0.19 (2.7)
Y isn
0.34 (3.4)
0.25 (3.1)
0.023 (2.9)
0.27 (3.2)
0.30 (3.8)
Zebrafish
Sim
0.34 (3.4)
0.27 (3.2)
0.11 (2.8)
0.16 (2.5)
0.29 (3.8)
Y 1:n
0.15 (2.5)
0.12 (1.9)
0.014 (1.1)
0.070 (1.5)
0.10 (1.8)
Reviews
0.14 (2.4)
0.14 (2.4)
-0.0085 (0.78)
0.14 (2.6)
0.12 (2.4)
Y 1:n
0.34 (10)
0.34 (10)
0.14 (9.3)
0.34 (1)
0.35 (10)
S&P 500
0.36 (9.4)
0.42 (11)
0.33 (13)
(11)
0.39 (11)
Y Ion
0.86 (1)
0.81 (2)
0,52 (8)
0.52 (8)
0.52 (8)
Simulated
C1
0.86 (1)
0.81 (2)
0.52 (8)
0.52 (8)
0.52 (8)
fixed, we see that the Opea error is larger than that for ddata This inferior performance of Öpea for very small and fixed n is explained by n appearing in the denominator of the second term in the rate Or(ynri hur/p+ / /r/n) for Öpen in theorem 2 versus appearing only in the numerator of Op(n3/9 /VD) for Ordata- Since it is Gaussian, this simulation model has finite exponentiallof-qqarraii moments, which is a much stronger condition than A4 we conjecture the convergence rate in this Gaussian case is Op( Vlog n/p) for d which would be consistent with figure 2 a). These numerical results seem to suggest the rate for dpea is similar, thus the second result of theorem may not be sharp.
P growing with
(b)
P fixed
(c)
fixed
102
logn
101
10
Adata
maximum error
p
ápca
mooo
10
10
10a
n 10
P 10a
109
109
200
400
600
800
1000
200
400
600
800
1000
200
400
600
800
1000
p
Figure 2: Simulation study ōdaat and åpeaa as estimators of a. All three subplots display the maximum error jej],e) |a Zi, 2,) á(i, for â & (black in all subplots (a)-(c) and = åpa (red). Error bars showing the standard deviation from 100 simulations are present for all data points, but in some cases small they are barely visible.
4.3 Comparing dot product affinities and Euclidean distances for the 20 Newsgroups data
In this section we expand on the results in table I for the 20 Newsgroups data, by exploring how inter-topic and intra-topic dot product affinities and Euclidean distances relate to ground-truth labels. Most existing agglomerative clustering techniques quantify dissimilarity using Euclidean distance. To compare dot products and Euclidean distances, figures ala)-bb show, for each topic, the top five topics with the largest average dot product and smallest average Euclidean distance respectively. We see that clustering of semantically similar topic classes is apparent when using dot products but not when using Euclidean distance. Note that the average dot product affinity between comp. windows.x and itself is not shown in figure (b)) but is shown in (a). by the highest dark blue square in the comp. windows.x column. In appendix D we provide additional numerical results illustrating that with fixed, performance in terms of T6 correlation coefficient increases with p.

---

comp
comp.sys
misc.forsaie
rec
rec.sport
religion
sci
talk.aaltics
(a)
Average dot products
(b)
Average Euclidean distances
be netween documents of each topic
between documents of each topic
8
8-000
0000
199
comp gra.gicspccc.
comp ou wwnddwssmii
compo syst syt.bbm.. pi hardwair
campa ma..aaa ww..ma......... w.c.harddarn
micc-forsalo
rec auttos
roc m.c.molorrycin
recr Lbaeebal
rec cee-pport hockey
S0 crypt
sci ccleleecrrmii sclime
selispace
tal poliiic.guus guns
talk talk,po poliiccc mideast
tolk politics alt athelsm soc. relioiiion christian wwl.gioo.gniiiccccccc...........
comp.graphin graphicss
comp me rindowws miso
comp sys iemm pc. hardwa
comp.wind cerip/iiidd windows
misc miiccfoorss forsalo
rec autox
roc.spart baseetol
rec sport wwt.t.toooocccooooe hox Rey
scl criccrpp
ellecetooonccs
melime
scl. xcl.spacc
tal pow.iiicssssss
mideast
talk. alt atthelsm
religiio www.chrtttinn
talk politiiiic poliiic
Average dot product affinity between
(c)
windows.x documents and other topics
(d)
windows.x dacuments and other tapics
Average Euclidean distance between
Average Euclidean distences
Average oot product
products
and other sublect arge
digranes
winduws other tech.
suijecss are informaire
so crycpp
comp sys.bbm hardward
comp sys mar www.maa......a...mmmmmmm hardwaard
electrenicc
scl-space
micc micc..orsal wec.forsalcc
talk, alt atthelsmm
gane
sel celimee talk politic FDC autoo
roc motoroyctiici les
soc religion christian,
roc. spprpportttoo bascbslll
palik.olltic..mm meliccmmiiie deast talkk.p leice.
mec.metrryccles
nic. autos
sci elleelecronicc
sclapacc
rec.ppor ppprrtbbscclooo bascbbll
comp amppgrraphiss scumas
misc foorsatte
polit mise
talk rel plan miso
pec.spoottt hackeey
soc roligion christion.
comp, sys, mac mac.maaaaaaa.. hardda worthward
talk,po midddossmiis mise
guns
att.athiemm comp sys.ssssttte ibmm pr hardwars sci crypt
talk www.iiic.mmmdd...... midrast
comp comp.gre
comp osims wii--widowww
talk rollg
(e)
ff)
Dendrogram of tapics using dat products
Dendrogram of topics using UPGMA
Yeligion misc chrabann mideast
metorcycles
autos
electronics
machhardwaar
guns
forsale
medd
graahics
motoorcccls
space
windovss.mi
aults
pollo Cs. misss
bassebill
med
acchardware
christian
electronics
hockey
winddowss.iisc
mideast
graphics
crypt
Figure 3: Analysis the 20 Newsgroups data. Marker shapes correspond to newsgroup classes and marker colours correspond to topics within classes. The first/second columns show results for dot products/Euclidean distances respectively. First row: for each topic (æ-axis), the affinity/distance (y-axis) to the top five best-matching topies, calculated using average linkage of PC scores be- tween documents within topics. Second row: average affinity/distance between documents labelled 'comp.wine windows.xx and all other topies. Third row: dendrograms output from algorithm I and UPGMA applied to cluster topics.

---

For one topic ('comp. windows.x') the results are expanded in figures (d) to show the average dot products and average Euclidean distances to all other topics. Four out of the five topics with the largest dot product affinity belong to the same 'comp' topic class and other one is semantically similar 'sei.eryp'' topic. Whereas, the other topies in the same 'comp' class are considered dissimilar in terms of Euclidean distance.
(
In order display visually compact estimated dendrograms, we applied algorithm and UPGMA in a semi-supervised setting where each topic is assigned its own PC score, taken to be the average of the PC scores of the documents in that topic, and then the algorithms are applied cluster the topics. The results are shown in figures (fe)-(f) (for ease of presentation, leaf vertex 'heights' are fixed to be equal).
Limitations and opportunities
Our algorithm is motivated by modelling assumptions. these assumptions are not appropriate for the data at hand, then the algorithm cannot be expected to perform well. notable limitation our model is that (u, v) for all u, V (see lemma 3 in appendix C). This is inappropriate assumption when there are strong negative cross-correlations between some pairs of data vectors, and may explain why algorithm has inferior performance the S&P 500 data in table [ Further discussion given in appendix that their computational cost scales faster than On²²). Approximations to standard agglomerative
criticism agglomerative clustering algorithms in their basic form is
methods which improve computational scalability have been proposed [33 35). Future research could investigate analogous approximations and speed-up of our method. Fairness in hierarchical
clustering has been recently studied in cost function-based settings by [6] and in greedy algorithm
settings by [13. Future work could investigate versions of our algorithm which incorporate fairness measures.
10

---

References
[1] Amazon product reviews dataset. https///www..gggee..cooo//tttttt
Accesed:: srrrarciccal-eeet-claaassatt | 2023-05-16..
[2] S&P500 stock data. http:://www..eeeerrrrr/co/aaaatt cessed: 2023-05-16.
Ac-
[3] S&P500 hierarchy. companies Accessed: 2023-05-16.
[4] Wikipedia list phylogenetics software. https: http:://enn...iikie wwt.:/wwwii.eda......aaaaaaaaaaa byylloeeneticc_sfftaaaa Accessed: 2023-05-16.
[5] Amir Abboud, Vincent Cohen-Addad, and Hussein Houdrougé. Subquadratic high-dimensional hierarchical clustering. Advances Neural Information Processing Systems, 32, 2019.
[6] Sara Ahmadian, Alessandro Epasto, Marina Knittel, Ravi Kumar, Mohammad Mahdian, Ben- jamin Moseley, Philip Pham, Sergei Vassilvitskii, and Yuyan Wang. Fair hierarchical clustering. Advances Neural Information Processing Systems, 33:21050-21060, 2020.
[7] Mihael Ankerst, Markus M Breunig, Hans-Peter Kriegel, and Jörg Sander. Optics: Ordering points to identify the clustering structure. ACM Sigmod record, 28(2):49-60. 1999.
[8] Matteo Barigozzi, Haeran Cho, and Dom Owens. Fnets: Factor-adjusted network estimation and forecasting for high-dimensional time series. arXiv preprint arXiv:2201.061 2022.
[9] Ricardo JGB Campello, Davoud Moulavi, and Jörg Sander. Density-based clustering based on hierarchical density estimates. In Advances in Knowledge Discovery and Data Mining: 17th Pacific-Asia Conference, PAKDD 2013, Gold Coast, Australia, April 14-17, 2013, Proceedings, Part II 17, pages 160-172. Springer, 2013.
[10] Gunnar E Carlsson, Facundo Mémoli, al. Characterization, stability and convergence of hicrarchical clustering methods. J. Mach. Learn. Res., 11(Apr):125511470, 2010.
[11] Moses Charikar and Vaggos Chatziafratis. Approximate hierarchical clustering via sparsest cut and spreading metrics. Proceedings of the Twenty-Eighth Annual ACM-SIAM Symposium on Discrete Algorithms,, pages 8-1-854. SIAM, 2017.
[12] Vaggos Chatziafratis, Mohammad Mahdian, and Sara Ahmadian. Maximizing agreements for ranking, clustering and hierarchical clustering via max-cut. International Conference on Artificial Intelligence and Slutistics, pages 1657-1665. PMLR, 2021.
[13] Anshuman Chhabra and Prasant Mohapatra. Fair algorithms for hierarchical agglomerative clustering. In 2022 2s IEEE International Conference on Machine Learning and Applications (ICMLA), pages 206-211. IEEE, 2022.
[14] Vincent Cohen-Addad, Varun Kanade, and Frederik Mallmann-Trent Hicrarchical clustering beyond the worst-case. Advances Neural Information Processing Systems, 30, 2017.
[15] Vincent Cohen-Addad, Varun Kanade, Frederik Mallmann- Trenn, and Claire Mathieu, Hierar- chical clustering: Objective functions and algorithms. Journal of the ACM (JACM), 66(4):1-42, 2019.
[16] Sanjoy Dasgupta. A cost function for similarity-based hierarchical clustering. In Proceedings of the forty-eighth annual ACM symposium on Theory of Computing, pages 118-127, 2016.
[17] William HE Day and Herbert Edelsbrunner. Efficient algorithms for agglomerative hierarchical clustering methods. Journal of classification, 1(1):7-24, 1984.
[18] Ronald W DeBry. The consistency of several phylogeny-inference methods under varying evolutionary rates. Molecular Biology and Evolution, 93)):337-5511 1992.
[19] Roland I. Dobrushin. Central limit theorem for nonstationary markov chains. i. Theory of Probability & Its Applications, 1(1):65-88 80, 1956.

---

[20] Justin Eldridge, Mikhail Belkin, and Yusu Wang. Beyond hartigan consistency: Merge distortion metric for hierarchical clustering. In Conference on Learning Theory, pages 588-608 PMLR, 2015.
[21] Edward B Fowlkes and Colin L Mallows. A method for comparing two hierarchical clusterings. Journal of the American statistical association, 78(383):553-669 1983.
[22] Sudipto Guha, Rajeev Rastogi, and Kyuseok Shim. Cure: An efficient clustering algorithm for large databases. ACM Sigmod record, 27(2) 1998.
[23] Lawrence Hubert and Phipps Arabie. Comparing partitions. Journal of classification, 2:193-218, 1985.
[24] Anil K Jain and Richard C Dubes. Algorithms for clustering data. Prentice-Hal Hall, Inc., 1988.
[25] Stephen C Johnson. Hierarchical clustering schemes. Psychometrika, 32(3):241 1967.
[26] Maurice George Kendall. Rank correlation methods. Griffin, 1948.
[27] Jisu Kim, Yen-Chi Chen, Sivaraman Balakrishnan, Alessandro Rinaldo, and Larry Wasserman. Statistical inference for cluster trees. Advances in Neural Information Processing Systems, 29, 2016.
[28] Daphne Koller and Nir Friedman. Probabilistic graphical models: principles and techniques. MIT press, 2009.
[29] Ken Lang. Newsweeder: Learning to filter netnews. In Machine learning proceedings 1995, pages 331-339, Elsevier, 1995.
[30] Steffen Lauritzen. Graphical models, volume 17. Clarendon Press, 1996.
[31] Bogdan-Adrian Manghiuc and He Sun. Hierarchical clustering: appr-oooooinnaooonn for well-
clustered graphs. Advances in Neural Information Processing Systems, 14:9278-9289 2021. 2021.
[32] Henok Mengistu, Joost Huizinga, Jean-Baptiste Mouret, and Jeff Clune. The evolutionary origins of hierarchy. PLoS computational biology, 12(6):e 1004829, 2016.
[33] Nicholas Monath, Ari Kobren, Akshay Krishnamurthy, Michael R Glass, and Andrew McCallum. Scalable hierarchical clustering with tree grafting. In Proceedings of the 25th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining, pages 1438-1448, 2019.
[34] Priya Moorjani, Carlos Eduardo G Amorim, Peter F Amdt, and Molly Przeworski. Variation in the molecular clock of primates. Proceedings of the National Academy of Sciences, 113(38): 10607-10612, 2016.
[35] Benjamin Moseley, Sergei Vassilvtiskii, and Yuyan Wang. Hierarchical clustering in general metric spaces using approximate nearest neighbors. In International Conference on Artificial Intelligence and Statistics, pages 2440-2448. PMLR, 2021.
[36] Fionn Murtagh. A survey of recent advances in hierarchical clustering algorithms. The computer journal, 26(4):354-33 1983.
[37] Fionn Murtagh and Pedro Contreras. Algorithms for hierarchical clustering: an overview. Wiley Interdisciplinary Reviews: Data Mining and Knowledge Discovery, 2(1):86 97, 2012.
[38] Fabian Pedregosa, Gaël Varoquaux, Alexandre Gramfort, Vincent Michel, Bertrand Thirion, Olivier Grisel, Mathieu Blondel, Peter Prettenhofer, Ron Weiss, Vincent Dubourg, et al. Scikit- learn: Machine learning python. the Journal of machine Learning research, 12:2825-2830, 2011.
[39] Magda Peligrad. An invariance principle for ---ixing sequences. The Annals of Probability, pages 1304 1313, 1985.
[40] William M Rand. Objective criteria for the evaluation of clustering methods. Journal of the American Statistical association, 666336):846-85, 1971.
12

---

[41] Chandan K Reddy and Bhanukiran Vinzamuri. A survey of partitional and hierarchical clustering algorithms. Data clustering, pages 87-110. Chapman and Hall/CRC, 2018.
[42] Alessandro Rinaldo, Aarti Singh, Rebecca Nugent, and Larry Wasserman. Stability density- based clustering. Journal of Machine Learning Research, 13:905, 2012.
[43] Aurko Roy and Sebastian Pokutta. Hierarchical clustering via spreading metrics. Advances in Neural Information Processing Systems, 29, 2016.
[44] R. R. Sokal and C. D. Michener. A statistical method for evaluating systematic relationships. University of Kansas Science Bulletin, 38: 1409-1434, 1958.
[45] Robert Sokal and James Rohlf. The comparison of dendrograms by objective methods. Taxon, pages 33 33-40, 1962.
[46] Baris Sumengen, Anand Rajagopalan, Gui Citovsky, David Simcha, Olivier Bachem, Pradipta Mitra, Sam Blasiak, Mason Liang, and Sanjiv Kumar. Scaling hierarchical agglomerative clustering to billion-sized datasets. arXiv preprint arXiv:2005.11653 2021.
[47] Daniel E Wagner, Caleb Weinreb, Zach M Collins, James A Briggs, Sean G Megason, and Allon M Klein. Single-cell mapping of gene expression landscapes and lineage in the zebrafish embryo. Science, 36066392)9981-98877 2018.
[48] Joe Ward Jr. Hierarchical grouping to optimize an objective function, Journal of the American statistical association, 58(301):236-244 1963.
[49] Nick Whiteley, Annie Gray, and Patrick Rubin- Subin-Delannch. Discovering latent topology and geometry data: law of large dimension, arXiv preprint arXiv:2208. 1665, 2022.
[50] Alexander Wolf, Philipp Angerer, and Fabian Theis. Scanpy: large-scale single-cell gene
expression data analysis. Genome biology, 19:1-5, 2018.
[51] Wang Xuejun, Hu Shuhe, Yang Wenzhi, and Shen Yan. On complete convergence for weighted sums of-mixing random variables. Journal of Inequalities and Applications, 2010:1-13, 2010.
[52] Tian Zhang, Raghu Ramakrishnan, and Miron Livny. Birch: an efficient data clustering method for very large databases. ACM sigmod record, 103-)1003-114 1996.
13

---

Appendices
A Wasserstein dimension selection
In order compute the principal component scores Ç1, Si...., Sn- must choose the dimension Traditionally and somewhat heuristically, this is done by finding the "elbow" in the scree plot of eigenvalues. The bias-variance tradeoff associated with choosing was explored by [49], who suggested a data-splitting method of dimension selection for high-dimensional data using Wasserstein distances. Whilst this method may be costly than the traditional "elbow" approach, was empirically demonstrated in [9] to have superior performance.
Algorithm 2 Wasserstein PCA dimension selection 49 Input: data vectors I....., Y, I..
1: for {1,..., min,(,)) do
2: Let V E IR"x" denote the matrix whose columns are orthonormal eigenvectors associated [m/2]
with the r largest eigenvalues of | Y.Y.
3: Orthogonally project I....., [m/22] onto the column space of Y, VV
4: Compute Wasserstein distance d, between Y....., Y(n/2) and Y [n/2]+1,-.. Y, (as point
sets in R²)
5: end for
Output: selected dimension argmin {d, }.
We note that in practice, the eigenvectors appearing in this procedure could be computed sequentially as grows, and to limit computational cost one might consider r only up to some Tmax < min(n, p).
B Implementation using scikit-learn
Algorithm [ can be implemented using the Python module scikit-learn 38 via their class, using the standard 'average' linkage criterion and custom metric function to compute the desired affinities. However, AgglomerativeClue ustering merges clusters based minimum distance metric, whereas algorithm merges according to maximum dot product. Therefore, the custom metric function we used in our implementation calculates all the pairwise dot products and subtracts them from the maximum. This transformation needs to then be rectified if accessing the merge heights.
AgglomerativeCluster
All code released part of this paper is under the MIT License and can be found at https
//githut www.//ooomnnnnnnieeee wttp///ninan//aa.oo.rrroeea http://t..ooouccceeearraaaaiiaae hierarchical
C Proofs and supporting theoretical results
Lemma 2. The height function h(v) = a(v, v) satisfies h(b)) 2 h(Pav) for all v E V except the root.
Proof.
1 [b]l[((v) X(Pa.) ) X(Pa.)
h(w))
(() X(Pa + 2 X(Pa.). X(Pa.)) + h(Pa.) 2 h(Pap),
p
where A2 combined with the tower property and linearity of conditional expectation implies
E X (v) X(Pa,,, X(Pa)))
Lemma 3. For all u. V, a(u,e) > 0.
14

---

Proof. = = a(u, v) 2 0 holds immediately from the definition of in 2 For u * with most recent common ancestor w,
B (
X(M Σ ( ==1
20 0,
where the second equality uses A logether with standard conditional independence aguments, and the third equality uses A2
Proof of lemmal 7 Let be the most recent common ancestor of and For each 1, ..... the property together with standard conditional independence arguments imply that Xú and (v( conditionally independent given X, (w), and the property A2 implies that (u)x(3 (w)] [[) = (2()). Therefore, by the tower property conditional expectation,
E (x)x))() = (x))
E
(w)?)
Hence, using the definitions of the merge height m. the height h and the affinity a,
P
m(u, v) h((w) at w. w) (
Σ (x) = a(u =((
which proves the first equality in the statement. The second equality is the definition of a. For the third equality the statement, we have
d(u, v) = (u)) + h(v) - 2h(w)
axu, u) + a(v, 2a(u,)) v)
E (() X (())) + 1 [ X(w)) - -2- 2'E(((( (xx)
P
E [||(u)) X())]
p
where the first equality uses the definition of d, and the second equality uses the definition of h and h(w) = m(u, v) = a(u, v).
C.1 Proof of Theorem
The following lemma establishes an identity concerning the affinities computed in algorithm I which will be used in the proof of theorem
Lemma 4. Let P.. m 2 0, be the sequence of partitions of [n]] constructed in algorithm Then for
any m 2 0,
a(u,u) |ull
1
Σ
a(i,)), for all distinct pairs u, v E Pm-
SEajER
Proof. The proof is by induction on m. With m = 0, 9 holds immediately since Po [1). ..... {n]}. Now suppose holds at step m. Then for any distinct pair w, w/ Pm+11
15

---

either w or w' is the result of merging two elements of Pm. or to and w' are both elements of Pm. In the latter case the induction hypothesis immediately implies:
d(w, w')
Σ
a(i,j).
AEW.JEW
In the case that or w' is the result of merge, suppose w.Loog. that u U for some u, E Pm and E Pm. Then by definition in algorithm
|u|
|e|
a(w, w/)
á
w")
+
(
w')
|u|
|0|
á(i,j)
Σ
a(i,j)
|u||2o|
[wl
iEN.JEW'
iEV.jEm
1
Σ
a(i,j);
iew.jEw
where the final equality uses w U... v. The induction hypothesis thus holds at step m 1.
The following proposition establishes the validity of the height function constructed in algorithm Some of the arguments used this proof are qualitatively similar those used to study reducible linkage functions by, e.g., Sumengen et al. [46), see also historical references therein.
Proposition 1. With v the vertex set and h the height function constructed in algorithm / with any symmetric, real-valued input a((;-). holds that h((v) h(Paa)) for all vertices E V except the root.
Proof. The required inequality h(v) 2 h(Pa.) holds immediately for all the leaf vertices v Po {[I)...., {n}} } by the definition of h in algorithm All the remaining vertices the output tree, i.e., those in V Po. are formed by merges over the course the algorithm. For m 2 0 let Wm = um Umn denote the vertex formed by merging some Um, Um E P Then Wm Pa and wm. Pa.m. Each um is either member of Po or equal to Wm' for some m' < m. The same is true of each Um. It therefore suffices to show that h(ww)) 2 h(wm+11 for m 2 O, where by definition in the algorithm, ((wm)) â(Um: Um). Also by definition in the algorithm, h(w,, h(wm+1) is the largest pairwise affinity between elements of Pm+1. Our objective therefore to upper-bound this largest affinity and compare it to M(wm) = á( fmm )m)
The affinity between Wm Um U Um and any other element w' of Pm+1 (which must also be an element of Pm) is, by definition the algorithm,
Um
[uml]
@(wm, w')
Ôl (uum w')
Ô(Un, w')
Umm + Fmml
[mmm + Pm
mmxxfc ( Um w') alVms w')3
6(Um: Um))
where the last inequality holds because Um, Um. by definition, have the largest affinity amongst all elements of Pm. For the same reason, the affimity between any two distinct elements of Pm+1 neither of which is wmn. (and therefore both of which are elements of P..) is upper-bounded by á(U V..). We have therefore established h(w,)) á( (um: vm) 2 h(wm+1) as required, and this completes the proof.
Proof of theorem] 7 Let us introduce some definitions used throughout the proof.
M max ((ii.) :) ar(2,, 2,)l- [
(10)
We arbitrarily chose and then fix j E [n] with i * j. and define
H
A = mi.jj).
(1I)
16

---

Let denote the most recent common ancestor of the leaf vertices { and {5} in D and let m 2 1 denote the step of the algorithm at which u is created by a merge, that m min 2 1: e
Pmn]) J. We note that by construction, u is equal to the union of all leaf vertices with ancestor and by definition of h in algorithm
(u)) = H.
Let denote the most recent common ancestor of 2; and 2j in D, which has height h(v) = H.
There is no partition of u into two non-empty sets A, B S [n]] such that á(k; 1) < A for all &s E A and I e B. We prove this by contradiction. Suppose that such partition exists. There must be a step m' S m at which some A' C A is merged some B' - B. The vertex formed by this merge would have height
Lower bound on m(2, -
h(w) - á(A',, B')
Σ
a(k, < 11 = hhw),
B'
REA'
where the first equality is the definition of h(w) in the algorithm and the second equality holds by lemma 4 However, in this construction u. is an ancestor of and h(w) < h(()) therefore contradicts the result of proposition !
As a device to be used in the next step of the proof, consider an undirected graph with vertex set u, in which there is an edge between two vertices h and I if and only if b(k,1) > H. Then, because there no partition as established above, this graph must be connected. Now consider a second undirected graph, also with vertex set u, in which there is any edge between two vertices is and l if and only if a12R 2) 2 Ā M. Due to the definition of M in (10). any edge in the first graph is an edge in the second, so the second graph is connected too. Let k, L, and É be any distinct members of u. Using the fact established in lemma I that azz,, 21) and a< ze) are respectively the merge heights in D between 2y and 2,, and 2 and žg, can be seen that if there are edges between A: and and between I and f in the second graph, there must also be an edge in that graph between & and f. Combined with the connectedness, this implies that the second graph is complete, so that a 2)) 2 H M for all distinct k, 2 E u. In particular a(z,, 23) 2 A M, and since m(z,, 22) = a(z,, 2y), we find
m( 2z) (i, 2 -M.
(12)
Upper bound on m(2,, 23) - m(.. Let S. = fi E [m] = e or has ancestor v in D). For k, E Sy, lemma I tells us (= 21) is the merge height between 2 and æ1, so a( #)) 2 H. Using 10). we therefore have
a(k, 2 H M, Vk, E S..
(13)
It follows from the definition of h in the algorithm that if S. [n],, the heights of all vertices in D greater than or equal to H M. This implies i 2 H - MM In summary, we have shown that when S, = [n]].
m( 22) m(.j)) S M.
(14)
remains to consider the case S. * [n]. The proof of the same upper bound (14 in this case is more involved. In summary, we need to establish that the most recent common ancestor of {i) and {f) in
D has height least .M. The main idea of the proof is to consider the latest step of the algorithm at which a vertex with height at least H M is formed by merge, and show the partition formed by this merge contains the most recent common ancestor of (i) and (j), or an ancestor thereof.
To this end let m² denote the latest step in algorithm | at which the vertex formed, w", has height greater than or equal to H M. To see that m² must exist, notice
max @(k, I) 2 alz, 2z) M, ky]E[[]]
(15)
by definition of M in (10). Combined with the definition of h in algorithm 1 the vertex formed by the merge at step 1 of the algorithm therefore has height greater than or equal to H M. Therefore m is indeed well-defined.
Our next objective to show that the partition Pm- formed at step m² contains an element which itself contains both and j. We proceed by establishing some facts about S, and Pm..
17

---

Let S, = [n] S.. For & E Se, E Se, U cannot be an ancestor of 21, by lemma (22,, 21) is the merge height of 2i and 24, and b the minimum branch length in D, we have at (Zk, 21) - b. From 10 we then find
a(k, S b E
(16)
We claim that no element of Pm. can contain both an element of S, and an element of S,,. We prove this claim by contradiction. If such element of P P.-- did exist, there would be step m' S m² at which some A' S Sy is merged with some B' C S,. But the vertex w' formed by this merge would be assigned height h(w')= á(A', B") < H b + M < H M, where the first inequality uses lemma and (T6). and the second inequality uses the assumption of the theorem that M b/2. Recalling the definition wn we have h(w*) 2 - M. We therefore see that w" is an ancestor of w' with h(w") h(w"), contradicting the result of proposition] !
Consider the elements of Pm. denoted A and B, which contain and j respectively. We claim that A = B. We prove this claim by contradiction. Suppose A o B. As established in the previous paragraph, neither A nor B can contain element of Sp. Therefore, using lemma and [3
a(A, B)
1
Σ KEAJEB
H
||||B]
Again using the established fact that no element Pm can contain both an element of S, and an element of S,,, m² cannot be the final step of the algorithm, since that would require Pm- {[]]] Therefore á( A, B) is one of the affinities which algorithm | would maximise over at step m + 1, so the height of the vertex formed by merge step m would be greater than equal H M, which contradicts the definition m Thus we have proved there exists an element of Pm- which contains both i and j. This element must be the most recent common ancestor of {i) and [i}. or an ancestor thereof. Also, this element must have been formed by merge at step less than or equal to m² and must have height greater than equal to H M. Invoking proposition I we have thus
established H H M. In summary, in the case S, + [n]] we have shown
m(2,,2 24) mi,j) < M.
(17)
Combining the lower bound (12) with the upper bounds 14 17 and the fact that i. chosen arbitrarily, completes the proof.
C.2 Supporting material and proof for Theorem 2
Definitions and interpretation for assumptions A3 and A5
We recall the definition of ç-mmixin from, e.g., define:
39 For sequence of random variables {Eji 2 1}.
(k)) = sup 121
sup
IPBBAA) F(B)).
where Fo the o-algebra generated by Ei,..., Then {5ir 2 1} is said to be yp-mixing (kk)) Y 0 as
To interpret assumption A5 notice
Z.1 < [S(Z)) E.1
max [S()) P E 1), #E2
where the first inequality uses the independence of E; and Z; and the second inequality uses the
fact that the elements of the vectors E, are i.i.d. Since Y. - x(Z,) = S(Z,) (2,,),, E,, A5 thus implies E X [ O(p) as p oo, which can be viewed as natural growth rate since is the dimension of the disturbance vector Y, ((Z)). the proof of proposition 2 below, AB is used in
similar manner to control dot products of the form (Y: X, X,) and (Y: X, Y X₃).
18

---

Proof of Theorem 2 For the first claim of the theorem, proposition combined with the tower property of conditional expectation imply that for any 8 >
P
max
P YY Y)) al2,, 2();, 26)
Isisisn
1 1 n 1)
1 89
C(g, E,,,), (18)
po/2
2
from which (6) follows.
The second claim of the theorem in essence corollary to [49] 1J. little work is needed to map the setting of the present work on to the setting of Whiteley et al. [Thm 1j. To see the connection, we endow the finite set Z in the present work with the discrete metric: dz (u, v) = 0 for v, and dz u) 0. Then (Z, dz) is compact metric space, and in the setting specified in the statement of theorem 2 where A3 strengthened to independence, 8 = p and S(v) = I for all Z, the variables Y
Y,. 7 (()), 2); E E, exactly follow the Latent Metric
Model of Whiteley et al. [49]
Zn are i.idd. according to probability distribution supported on 2. As [49], by Mercer's theorem there exists feature map Φ 2 -> IR'" associated with this probability distribution, such that (Φ(u), x(v)) a(u, v), for
Moreover, according to the description in section 2.1 the variables ,
u, 2. Here r, as A6 is the rank of the matrix with elements au,, v), which most 2.
Theorem of 49 in this context implies there exists a random orthogonal matrix Q E R such that
max *E
1/2Q o(zi) E Op
(VT
nr +
(19)
Consider the bound:
1 -
pápea (i, (x(2x,
(cr, S) - ((2)))
D
p 1/-2 o(2.), (0 1/20))
(0(Zz), x- 1/2 --/00-10000000 Q2Q - -d(Z,))
1 1/2 QG: ((z) (11-1/2Q6 - - +|1((2)11)
[() |00 1/2 1/2QC: - -ZZ)
where orthogonality of has been used, and the final inequality uses Cauchy-Schwarz and the triangle inequality for the || norm. Combining the above estimate with (19). the bound:
max [l((2)) < max [1)()]] = max auv *E(n]
DE2
vE2
S sup max [(IX,(())3]
< sup max j21 vEZ
(1)x (20)
j21
and completes the proof of the second claim of the theorem.
Proposition 2. Assume the model in section 2.1 satisfies assumptions A3 and let ip and g be as in
A3 and A4 Then there exists a constant C( sp) depending only on q and sp such that for any >
P
max Isi<JSn
1 (Y, Y,) aZ7, 23)l 2 8 7 71....; 7m
1 1 n(n -
((, (21)
89 pn/2
19

---

where
Mgg. X, E) = sup max E [
121 MEZ
E 1111 9]]
sup max |S(())lbpp p21 nEZ
mup
sup max (1))1 121 NEZ
E
E1
2q
sup max [S(e)lopp p21 REZ
Proof. Fix any i.j such that <i<js n. Consider the decomposition:
p Y, - Zz)
A
where
A = p-' -X1)) - aa Ej)
A₂ = p- (X(2),
A3 = p'' (x(2)),
A1 = p-' E, S(2,)
The proof proceeds by bounding A (1 7 Zn] for ls 1, 4. Writing A as
A
11
ALA Xxx=Xx(2)Xxx3))) (X
Xx
Zn]
we see that A₁ a sum p random variables each of which conditionally mean zero given Zz,..., Zn- Noting that the two collections of random variables (Zi.... Zn} and (X(v): E V) are independent (as per the description of the model in section 2.1p. under assumption AB we may apply a moment inequality for se-mixing random variables [Sm]ma 1.7] to show that there exists constant Ci (q, s) depending only on ip such that
E 11 E1
Zn]
0/2
1
P
S <i((p)
E
IZ...............
19 Zi....
Zn] +
| 1,4
|| Zi.... Zn|
7
n12
P
p
<i((4))))
A 1,k E
21..... Za] +
E [| 1.& 1 21 ....... Znl
109
pa/2
k=1
200i(9.44)
sup E
A 1
91
Z1
Zul
pA/2
tz1
C1 19)
sup max E k21 vEZ
[x)
(22)
2ala
where second inequality holds by two applications of Jensen's inequality and 2 2, and the final inequality uses the fact that for a, b 2 0. (a + b)d < 29-1(a" + 69), the Cauchy- -Sehhy--- inequality, and the independence of [Z,.... Zn} and {X(
For A2. we have
A = P k=1
A = X(Za))
where H) denotes the kth element of vector. Since the three collections of random variables, {Zi,...; Znh. (X(v); E 2) and (EI,..., E.,) are mutually independent, and the elements of each vector E, RRP are mean zero and independent, we see that given (Z,..., Zn) and (X(e); e VJ. Ag is simple average of conditionally independent and conditionally mean-zero
20

---

random variables. Applying the Marcinkiewicz-Zynnne Zygmunn inequality we find there exists a constant C2(4) depending only on such that
E 21,,,1922,.... Und X 2]
2,,..,,
9/2
S C2(g)
2
2
7 Zn X(u)0
(23)
Noting that and applying Minkowski's inequality to the r...s. of (2 then using the indepen- dence of [Zi...., Zn (X(e);u E 2} and the vector Ej.
and the i.dd. nature of the elements of
19/2
2/4
E
A 2, 1
Zi...., Zno X(u);UE
2, 91
X,, 278/9
E [[S(2,) (x( X(2)){a |E Zı, 21,...,Zu, X(v)uue 2/9
Σ
p
1 pa
E 2// n
E
[S(2,) T x(2)))1
E
2// ||(())) X(Z) 1 [
pa
max 2/2
2EZ
Substituting into
23 and using the tower property of conditional expectation we obtain:
[
19
Z
Zp]
1
E
E11
max DEZ
|S())
X(Z))l
ZI...,,Z! Zn]
pa/2
pa/2
9/2)
P (x,(2))
E11
9
max DE
|S( E
Z Zi...., On
pa/2
k=1
1
P
E1
max 2
|S
E
X 21
ZI.... Un
p//22
1
E.
E1
|S()) sup max TE (+) 121 vez
(24)
pA/2
vEZ
where the second inequality holds by Jensen's inequality (recall 2). Since the r.h.s. of (24 does not depend on i or j, the same bound holds with Ag on the l.h.s. replaced by Ag.
Turning to A., we have
1 A p
((Z))),
Σ ISRRPSP
= S(Z)) S(()))kee
P
Noting that * and that the elements of E, and E, are independent, identically distributed, and mean zero, we that A is sum of p² random variables which are all conditionally mean zero and conditionally independent given ......,y The Marcinkiewiez-Zyemmnn inequality gives:
9/2
E E| ZI,
Z. S E2(g)
1
Σ
Zi.... Zn
(25)
21

---

Applying Minkowski's inequality to the r...s. of
25
9/2
2/9
E
Z1
Zn
p
I<k
E E 1 Zi.... :
2/1
p2
2/9
E
Eu
Ejel" [[((2)) S(2))loo))
Z1, 2.y.... Zn
ISk,ESp
1 EE
[[11
19/4//
E | k(2))
pd
PSP
1 E
Em
20
2/9
max S S(e)l 0
where the final inequality holds by Jensen's inequality. Substituting back into 22 and using
[S( S(w)]le $ sp' (S() S(u))lop S // (/]|[((
we obtain:
E A Z1, ,,,--.. ... Zml C2(2)
E
En
24
max WEZ
|S(())l
(26)
p9/2
Combining 22
24 and
26 using the fact that for b 2 O, (a + 6)S
S ( + b9), we find
that there exists a constant
Clq. (p) depending only on and ip such that
E
( Y) 2(),, Zz Zn
$ C(g. p)) pa/2 (q,X,,
where (q. X, E) is defined in the statement of the proposition and is finite by assumptions A4 and A5 By Markov's 's inequality, for any 2
P ((p Y,, Y) - at Z4, 2z) 2 8 Z1, Zn) S _ C(g.s) M(g.X,E,5 E,
(27)
50
Daf2
and the proof is completed by union bound:
P
max 1SicsSn
Y Y,) a(Z2, < Zi,... 2m
p
YY, Y,) al 2,) < Zi...., ... Zn
Izicjšn
U P 1 (Y,, - (21,, EI...., Un
1 P
ISieiSn
[2n
1)
2 1
Clg. 4)
X. E.S).
89
pA/2
C.3 Interpretation of merge heights and exact tree recovery
Here we expand on the discussion in section 3.3 and provide further interpretation of merge heights and algorithm[] [ In particular our aim is to clarify in what circumstances algorithm I will asymptoti- cally correctly recover underlying tree structure. For ease exposition throughout section we assume that 2 are the leaf vertices of T.
As a preliminary we note the following corollary theorem D assuming > 0. if one takes as input algorithm [ the true merge heights, i.c. (up bijective relabelling of leaf vertices)
22

---

á(:.-) = m(..)) a(+,), ), where = 12 then theorem implies that algorithm 1 outputs dendrogram D whose merge heights ( are equal to m(..) (up to bijective relabeling over vertices). This clarifies that with knowledge of m(::)) algorithm constructs dendrogram which has m(:)) )) as its merge heights.
We now ask for more: if once again -) is taken as input to algorithm/ | under what conditions is the output tree I equal to T (upto bijective relabelling of vertices)? We claim this holds when T is binary tree and that all its non-leaf nodes have different heights. We provide sketch proof of this claim, since complete proof involves many tedious and notationally cumbersome details.
To remove the need for repeated considerations of relabelling, suppose T (V,E) E is given, then w.l.o.g. relabel the leaf vertices of 7 as {1}.... 2 and relabel each non-leaf vertex to be the union of its children. Thus each vertex is some subset of [|]].
Now assume that T is a binary tree and that all its non-leaf nodes have different heights. Note that |V 2| - i..., there are | non-leaf vertices. The tree T is uniquely characterized by sequence of partitions Po...., Piz -1 where Po ((1).... {121), and for
....., 2
Pm is constructed from 1 1 by merging the two elements of Pm-1 whose most recent common ancestor is the mtth highest non-leaf vertex (which is uniquely defined since are assuming two non-leaf vertices have equal heights).
To see that in this situation algorithm with a(:,)) m(::)) and |2] as input, per- forms exact recovery of the tree, i..., T T, it suffices to notice that the sequence of partitions Po,... PIzI-1 constructed by algorithm uniquely characterizes T, and moreover (Po...., -1) =
P1))-)). The details of this last equality involve simple but tedious substitutions of m(.)) in place of a(:,.) algorithm so are omitted.
D Further details of numerical experiments and data preparation
All real datasets used are publicly available under the CCO: Public domain license. Further, all experiments were run locally on laptop with an integrated GPU (Intel UHD Graphics 620).
D.1 Simulated data
For each v E V, (,(v) are independent and identically distributed Gaussian random variables with:
X,(1) ~ ((X,()),)
Xy(2) (6),(
X,(3) ~ (6)2)
X,(4) (7).5)
X,(5) ~ ~N(X((7)
X,(6)~N(,) (8),
~ 18),
X,(8) 1(0,
for I.....p.
D.2 20 Newsgroups
The dataset originates from 29]. however, the version used is the one available in the Python package 'scikit-leam' 38 Each document is pre-processed in the following way: generic stopwords and e-mail addresses are removed, and words are lemmatised. The processed documents are then converted into a matrix of TF-IDF features. Labels be found the 20 Newsgroups website http:://qont
jasoaa/llelssgropsss..t but are mainly intuitive from the title of labels, with full stops separating levels of hierarchy. When using PCA dimension of = 34 was selected by the method described in appendix [
The following numerical results complement those in the main part of the paper.
23

---

Table 2: Kendall т. ranking performance measure, for Algorithm 1 and the 20 Newsgroups data set. The mean Kendall 7b correlation coefficient reported alongside the standard error (numerical value shown is the standard emrox 109). This numerical results are plotted in figure below.
Data
Input
500
p 1000
p = 5000
7500
p 12818
Y Iim
0.026 (0.55)
0.016 (1.0)
0.13 (2.2)
0.17 (2.5)
0.26 (2.9)
Newsgroups
Sis
0.017 (0.72)
0.047 (1.2)
0.12 (1.9)
0.15 (2.5)
0.24 (2.6)
Performance of Alg. 1 for the 20 Newsgroups data, with n raxed and varying
0.20
T. correlatio
0.15
0.10
0.05
2000
4000
6000
0000
10000
12000
Figure 4: Performance of Algorithm for the 20 Newsgroups data set as a function of number of TF-IDF features, p, with fixed. See table for numerical values and standard errors.
D.3 Zebrafish gene counts
These data were collected over a 24-hour period (timestamps available in the data), however, the temporal aspect of the data was ignored when selecting a sub-sample. To process the data, we
followed the steps in [47] which are the standard steps the popular SCANPY (50 Python package for analysing single-cell gene-expression data to process the data. This involves filtering out
genes that are present in less than cells or are highly variable, taking the logarithm, scaling and regressing out the effects of the total counts per cell. When using PCA a dimension of r 29 was selected by the method described in appendix A
D.4 Amazon reviews
The pre- pre-processing of this dataset is similar that the Newsgroup data except e-mail addresses are no longer removed. Some data points had a label of 'unknown' in the third level of the hierarchy, these were removed from the dataset. In addition, reviews that are two words or less are not included. When using PCA dimension of r = 22 was selected by the method described in appendix A
D.5 S&P 500 stock data
This dataset was found through the paper 8 with the authors code used to process the data, The labels follow the Global Industry Classification Standard and can be found here: [3]. The return on day of each stock is calculated by
op Pt
return;
P;
where p."h and pel is the respective opening and closing price on day i. When using PCA a dimension of 10 was used selected by the method described in appendix
24

---

D.6 Additional method comparison
Table T reports additional method comparison results, complementing those in table which con- cerned only the average linkage function (noting UPGMA is equivalent to using the average linkage function with Euclidean distances). In table we also compare Euclidean and cosine distances paired with complete and single linkage functions. To aid comparison, the first column (average linkage with the dot product) is the same in T In general, using complete or single linkage performs worse for both Euclidean and cosine distances. The only notable exception being slight improvement on the simulated dataset.
Table 3: Kendall тà ranking performance measure. For the dot product method, i... algorithm I Y as input corresponds using ddata and Çi:m corresponds ápea. The mean Kendall т. correlation coefficient is reported alongside the standard error (numerical value shown is the standard error x10%).
Average linkage Dot product
Complete linkage
Single linkage
Data
Input
Euclidean
Cosine
Fuclidean
Cosine
0.26 (2.9)
0.022 (0.87)
-0.010 (0.88)
-0.0025 (0.62)
-0.0225 (0.62)
Newsgroups
<1m
0.24 (2.6)
0.0041 (1.2)
0.036 (1.2)
-0,016 (2.0)
0.067 (1.5)
11
0.34 (3.4)
0.15 (2.2)
0.24 (3.2)
0.023 (3.0)
0.032 (2.9)
Zebrafish
Sim
0.34 (3.4)
0.17 (2.0)
0.30 (3.4)
0.12 (2.8)
0.15 (3.2)
Reviews
1im
0.15 (2.5)
0.019 (0.90)
0.023 (1.0)
0.0013 (0.81)
0.0013 (0.81)
0.14 (2.4)
0.058 (1.5)
0.063 (1.8)
0.015 (1.2)
0.038 (1.0)
Y 1:o
0.34 (10)
0.33 (10)
0.33 (10)
0.17 (10)
0.17 (10)
S&P 500
0.36 (9.4)
0.32 (10)
0.31 (10)
0.36 (13)
0.39 (12)
Y 1om
0.86 (1)
0.55 (8.7)
0.84 (2.0)
0.55 (8.7)
0.84 (2.0)
Simulated
mm
0.86 (1)
0.55 (8.7)
0.84 (2.0)
0.55 (8.7)
0.84 (2.0)
Table 4: Kendall ranking performance measure. The mean Kendall т correlation coefficient is reported alongside the standard error (numerical value shown is the standard error 109).
Data
Input
UPGMA with dot product *dissimilarity*
UPGMA with Manhattan distance
Y Im
-0.0053 (0.24)
-0,0099 (1.3)
Newsgroups
Smm
0.0029 (0.33)
0.052 (1.6)
0.0012 (3.13)
0.16 (2.4)
Zebrafish
S1:n
0.00046 (0.12)
0.050 (2.8)
Y 1:n.
-0.0005 (0.29)
0.0018 (0.44)
Reviews
Sim
-0.0015 (0.41)
0.061 (1.3)
Y 1m
0.0026 (7.7)
0.37 (9.4)
S&P 500
Sim
0.0028 (7.5)
0.39 (11)
1m
-0.0026 (1.6)
0.55 (8.7)
Simulated
Slon.
-0.0023 (1.8)
0.84 (2)
E
Understanding agglomerative clustering with Euclidean or cosine distances in our framework
E.1
Quantifying dissimilarity using Euclidean distance
The first step of many standard variants of aggiomerative clustering such UPGMA and Ward's method is to find and merge the pair of data vectors which are closest to each other in Euclidean distance. From the elementary identity:
YY: Y 2 YY 2 + |2 2(Y, Yj),
we see that, general, this not equivalent finding the pair with largest dot product, because of the presence of the terms Y I|2 and || Y For some further insight in to how this relates to
25

---

aa) m(a, d) < m(c,d) d) & d(a, d) > dec d)
(b) mía, d) < mmfed) d) & d(a, d) < d(c, d)
b d
g
Figure 5: Illustration of how maximising merge height m(::-) may or may not be equivalent to minimising distance d(::,)) -), depending the geometry of the dendrogram. (a) Equivalence holds (b) Equivalence does not hold.
modelling and theoretical analysis framework, it is revealing to the consider the idealised case of choosing to merge by maximising merge height m(...) versus minimising d(:::) (recall the identities for and d established in lemma D Figure shows two simple scenarios in which geometry of the dendrogram has an impact on whether or not maximising merge height m(::)) :) is equivalent to minimising -. From this example see that, in situations where some branch lengths are disproportionately large, minimising d(..-) will have different results to maximising m(...).
UPGMA and Ward's method differ in their linkage functions, and so differ in the clusters they create practice after their respective first algorithmic steps. UPGMA uses average linkage to combine Euclidean distances, and there does not seem to be mathematically simple connection between this and algorithm| D except to say that in general will return different results. Ward's method merges the pair of clusters which results in the minimum increase in within-cluster variance. When clusters contain equal numbers of samples, this increase is equal to the squared Fuclidean distance between the clusters' respective centroids.
E.2 Agglomerative clustering with cosine distance is equivalent to an instance of algorithm 1
The cosine 'distance' between Y, and Y, is:
(i,)) Y Y
In table [ we show results for standard agglomerative clustering using deas as a measure of dissimi- larity, combined with average linkage. At each iteration, this algorithm works by merging the pair of data-vectors/ccussers which are closest with respect to dess say and merged to form w, with dissimilarities between and the existing lata-vectors/clusters computed according to the average linkage function:
|0|| |w|
doos (w. -)
deas (u, .) doos (v. -).
(28)
This procedure can be seen to be equivalent to algorithm ! with input a(.,)) 1- doos (..)). Indeed maximising 1 - deos (., -) is clearly equivalent to minimizing deos(-,)) and with a(: = 1 doos (-:))
the affinity computation at line 6 of algorithm D ii::
[2]
ā (w, -)
ā(u, -) + [cc|
a(v, -)
[xc]
|0||
[1 doos (u, }} +
[] dece(u -)]
|w|
|u]]
|00|
[
|2| deos (v,
deos (2, -)
|w]]
|wo|
deos(u, -)
doos
deos (1, :)
|w]
w
where the r.h.s. of the final equality recognise
(28
26

---

E.3
Using cosine similarity an affinity measure removes multiplicative noise
Building from the algorithmic equivalence identified section E.2 we now address the theoretical
performance of agglomerative clustering with cosine distance in our modelling framework.
Intuitively, cosine distance is used in situations where the magnitudes of data vectors are thought not to convey useful information about dissimilarity. To formalise this idea, we consider a variation of statistical model from section 2. in which:
Σ the leaf vertices of T, | n,, and take these vertices to be labelled 2 {1, ........ ,.....) properties AT and A2 hold, and it is assumed that for all v E 2,
X is as in section 2.1
p X
the additive model
is replaced by multiplicative noise model:
Y, (i),
(29)
where % > 0 are all strictly positive random scalars, independent of other variables, but otherwise arbitrary.
The interpretation of this model is that the expected square magnitude of the data vector Y; is entirely determined by indeed we have
- >f ?
[Bll
Y Y
p
where h is as in section
2.1
We note that in this multiplicative model, the random vectors E, and 2.1 play role, and can view the random variables Zi,..., Zn as
matrices (()) from section
being replaced by constants Z, i, rather than being i.i.d.
Now define:
(Y Y,) =1
coo (i.
Theorem 3. Assume that the model specified in section sup max,ez (V))( < oc. Then
E.3 satisfies
3
and for some q
nz/9
max <jE(n]aa
|a(i, &ens (1,5)) OF
Vp
Proof. The main ideas of the proof are that the multiplicative factors Y. in the numerator dccs (i, cancel out with those in the denominator, and combined with the condition "='[[[[(())I1 1 we may then establish that a j) and (coe(i, are probabilistically close using similar arguments to those in the proof of proposition
Consider the decomposition:
p (X(())),
a(i,j) acos (i,)) aii.))
p 1/2 1-1/21 X
1 (X(1),
= aii,,)
(30)
(X(I)X()))
(1)lp 1120011111110000) X(2xX( P 1-11
(31)
Applying the Cauchy-Sehie Schwartz inequality. and adding and subtracting p 1/2 [x2 and in the final term of this decomposition leads to:
max kJE(n],
(a(i.)) &cau (i.5)
< mAX 1.)
((i.j) 1
((XI).X(M)
(32)
p
+ max sE/n]
1/2 [((i)] -
(33)
+ 1 max iE|n]
p 1/1/21[X(()11
-
-1l)
max je(mi
P
1/2 1112 | x()iI - -1|-
(34)
27

---

The proof proceeds by arguing that the term 32 is in Op (n2/0/ /VD) and the terms 33 and 34
in Oş (n'// /p), where we note that this asymptotic concerns the limit n / //D 0.
In order to analyse the term (2) 3. let denote the probability law of the additive model for ......; Y,,
in equation expectation. Then for any 6 0 and i, E [m]. ė 7 J.
T in section 2.1 in the case that (v) for all v E Z. Let denote the associated
(a(l.))
1 (X1), x()) > 8)
P
p
=P
- YY. ,) - a a2 2)) > 8 21 = 1a,....... Nn. -....,,n. =n
1
E 1121 1.... Zn =njpp (Y, Y,) - -a(Zi, 2z) > Z1 Zi.... 2n) -6,.....,Z)) 7, =)]
11
The conditional probability on the r.h.s the final equality can be upper bounded using the inequal- ity (27) in the proof of proposition 2 and combined with the same union bound argument used immediately after (27). this establishes (32 is in Op (n2/4/ vP) as required.
To show that (33 and are in Oj (T// it suffices to show that maxte[]] |--1/ 1
in Oe (n (n ///))) By re-arranging the equality: (la) ()a 1) |a|2
have
1/2 X (1)l - 1| 1 | =
(1
(35)
where 4,(1) = (x,(1) Thus under the model from section E.3
Σラ= (1) is an average
of mean-zero random variables, and by the same arguments as in the proof of proposition 2 under the mixing assumption A3 and and the assumption of the theorem that sup,21 max EE [ 29 oc, combined with a union bound, we have
max | [X()) - 1| Op(n /p)) #E[n]
Together with 35 this implies maxie[n] |p11/2| - 1 is in Oz (n// required, and that completes the proof.
E,4 Limitations of our modelling assumptions and failings of dot-product affinities
As noted in section 5 algorithm is motivated and theoretically justified by our modelling assump- tions, laid out in sections and E.3 If these assumptions are not well matched to data in practice, then algorithm D may not perform well.
The proof of lemma CommmaBshone shows that, as consequence of the conditional independence and martingale- like assumptions, Alland AQ au,,)) 2 0 for all u, By theorems 2 and 3] ådata åpea approximate a (at the vertices 2) under the additive model from section 2.1 and óeess approximates a under the multiplicative model from section F.3 Therefore if in practice (i, j). (pea(f, âcou(i, j)
are found to take non-negligible negative values for some pairs i, j, that is an indication that our modelling assumptions may not be appropriate.
Even if the values taken by (d (i, (i, j) or dooa (i, are all positive in practice, there could be other failings of our modelling assumptions. As an academic but revealing example, if instead of or 29 the data were to actually follow a combined additive and multiplicative noise model, i...
Y, X(,) + E,,
then in general neither dettaa ppea nor Ôcas would approximate p Co.
28


Reviews:
[{"review_id": "oc8yLXdn1s", "reviewer": "Reviewer_dt3f", "summary": "The paper discusses a new perspective on hierarchical clustering using the dot product as similarity measure instead of some distance. Under mild conditions on the probabilistic graphical model that generates the data, the proposed algorithm is shown to faithfully recover the underlying tree geometry. Surprisingly, the theoretical results show that the performance of the approach does not suffer, but instead benefits from high-dimensional data.", "questions": "- In line 102, it is not clear why $\\mathbb{E}[< Y_i, Y_j>|Z_1,\\dots,Z_n]$ is conditioned on all $Z_n$\n- What is the meaning of the operator $O_{\\mathbb{P}}$?\n- In Table 1, the results for simulated data do not depend on whether PCA was applied or not. Is this caused by the fact that the dimensions are independent?\n- Since the PCA version of the algorithm operates on fixed dimension $r$, it is not clear to my why in (9) an increase of $p$ should bring any benefit. Could you expand on that?", "limitations": "The authors clearly outline the limitations of their study, which is highly appreciated. A possible further limitation of the method is that, apparently, good performance guarantees can only be given for high-dimensional data, which is slightly counterintuitive. Societal implications are not to be expected.", "rating": 8, "confidence": 3, "soundness": 4, "presentation": 4, "contribution": 3, "strengths": "The paper is well written and all theoretical results and definitions are accompanied by intuitive explanations (e.g., lines 84-90, lines 159-166, Fig. 1, or Sec. 3.3). This makes the paper accessible to an audience not versed in graphical models, for example. The results seem quite promising; it is noteworthy that the authors also included \"negative\" results for the S&P500 dataset, and that they made an effort in investigating the limitations of the proposed method.", "weaknesses": "I see very little weaknesses in the paper (see below for a list of typos). There are only few points unclear (see Questions below). The only concern I have is that the results claim that the performance of the method improves if the dimensionality of the data increases. This is probably caused by the mixing condition, which ensures that more dimensions deliver more relevant information. This assumption is quite questionable in my opinion, as usually more dimensions do not reveal substantially more information about the underlying distribution. As an example, consider natural images: going from low resolutions to higher resolutions may at first improve clustering performance, but at some point increasing the resolution will bring only diminishing returns. The reason is that, as the resolution increases, also the correlation between dimensions increases, which is in contrast to A1. I would appreciate seeing results based on, e.g., lower numbers of TF-IDF features to understand how this affects clustering performance.\n\nAnother concern is that Fig. 3a and 3c seem to be inconsistent. In 3a, at the position of comp.windows.x, there are three dark blue squares, one light blue square, and (maybe?) a hidden pink x. In 3c, the top five topics contain two dark blue and two light blue squares. Please check if there was a mistake in generating the figure.\n\n## Minor:\n- line 56: \"interpretation to in all...\"\n- line 78: \"distributional properties OF $\\mathbf{X}$\"\n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": ">The only concern I have is that the results claim that the performance of the method improves if the dimensionality of the data increases. This is probably caused by the mixing condition ... I would appreciate seeing results based on, e.g., lower numbers of TF-IDF features to understand how this affects clustering performance.\n\nThanks for this insightful comment. As requested, in the .pdf attached to this response, table 1 and figure 1, we give additional numerical results showing the effect of the number of TF-IDF features, i.e. varying $p$. These results were obtained by randomly choosing features to exclude, to give the desired $p$. We find that, as $p$ grows performance improves, without evidence of \"diminshing returns\". We note a slight increase in the standard error with $p$, but the standard errors are roughly 100 times smaller than the values of the performance measure (note the factor of $10^3$), so we think this is not a serious concern. \n\nThe mixing condition is indeed a key assumption. However, it is not the only reason that performance improves with dimension. To see why, note that from eq. (1), \n\n$$\np^{-1}\\langle\\mathbf{Y}_i,\\mathbf{Y}_j\\rangle = p^{-1}\\langle\\mathbf{X}(Z_i),\\mathbf{X}(Z_j)\\rangle + p^{-1}\\langle\\mathbf{X}(Z_i),\\mathbf{S}(Z_j)\\mathbf{E}_j\\rangle + p^{-1}\\langle\\mathbf{X}(Z_j),\\mathbf{S}(Z_i)\\mathbf{E}_i\\rangle + p^{-1} \\langle \\mathbf{S}(Z_i)\\mathbf{E}_i,\\mathbf{S}(Z_j)\\mathbf{E}_j\\rangle.\n$$\n\nThe mixing assumption allows us to prove $p^{-1}\\langle\\mathbf{X}(Z_i),\\mathbf{X}(Z_j)\\rangle-\\alpha(Z_i,Z_j)\\to 0$ (convergence in probability) as $p\\to\\infty$, whilst the properties that the \"noise\" vectors $\\mathbf{E}_i$ are zero-mean, have independent entries, are independent of each other, $\\mathbf{X}$ and $Z_1,\\ldots,Z_n$, allow us to deduce that the second to fourth terms on the r.h.s. of the displayed equation above converge to zero as $p\\to\\infty$.  In this latter sense, high-dimensionality can help remove noise, irrespective of whether the mixing condition A1 holds.\n\nWe agree that for image data as you describe, where $p\\to\\infty$ corresponds to increasingly high resolution, the mixing assumption is not realistic. Future work would be needed to investigate alternative assumptions in complete detail, but for example, in such situations it could possibly be realistic that $p^{-1}\\langle\\mathbf{X}(Z_i),\\mathbf{X}(Z_j)\\rangle-\\alpha(Z_i,Z_j)$ is \"close\" to zero with some high probability for all $p$ large enough, without actually converging to zero. Another interesting possibility is that if one does _not_ assume the mixing condition A1, then one might consider $p^{-1}\\langle\\mathbf{X}(Z_i),\\mathbf{X}(Z_j)\\rangle$ as an affinity measure rather than $\\alpha(Z_i,Z_j)$. Of course these are just hypotheses, but we believe the work in the manuscript could be a first step towards understanding such scenarios.\n\nFor some types of data, such as time series on increasingly long time scales, or geospatial data collected over increasingly large areas, $p\\to\\infty$ does _not_ correspond higher resolution, and the mixing assumption is arguably realistic. For other data types, such as the document and biological data sets in our paper, it is perhaps less easy to firmly rule the mixing assumption in or out, but we believe our numerical results are encouraging.\n\n>Another concern is that Fig. 3a and 3c seem to be inconsistent. ... Please check if there was a mistake in generating the figure.\n\nThe explanation here is that in 3c we do not plot the average dot product affinity between comp.windows.x and itself. However, this is plotted in 3a; the highest dark blue square in the comp.windows.x column corresponds to dot-product with itself. We could update this, or add a comment to clarify.\n\n>**Minor:**...\n\nThanks for catching these typo's. \n\n>**Questions**\n> - In line 102, it is not clear why  $\\mathbb{E}[\\langle Y_i,Y_j\\rangle|Z_1,\\ldots,Z_n]$ is conditioned on all $Z_1,\\ldots,Z_n$\n\nUnder our model indeed $\\mathbb{E}[\\langle Y_i,Y_j\\rangle|Z_1,\\ldots,Z_n]=\\mathbb{E}[\\langle Y_i,Y_j\\rangle|Z_i,Z_j]$, we would happily amend line 102 to the latter.\n\n> - What is the meaning of the operator $O_\\mathbb{P}$?\n\nIntuitively, the argument of $O_{\\mathbb{P}}(\\cdot)$ is the convergence rate. We would happily add a precise defintion of this \"big Oh in probability\" notation to the manuscript: if $X_{p,n}$ is random variable indexed by $p$ and $n$, then e.g. $X_{p,n}\\in O_{\\mathbb{P}}\\left(n^{2/q} / \\sqrt{p}\\right)$ means that for any $\\epsilon>0$ there exists finite $\\delta$ and $M$ such that $n^{2/q} / \\sqrt{p}\\geq M$ implies:\n\n$$\n\\mathbb{P}(|X_{p,n}|>\\delta)<\\epsilon.\n$$\n\n> - In Table 1, the results for simulated data do not depend on whether PCA was applied or not. Is this caused by the fact that the dimensions are independent?\n\nFrom a rigorous mathematical point of view, we cannot answer this question definitively: it is possible that a conclusion similar to that in Theorem 2 for PCA holds with the mixing assumption A1 rather than independence assumption, but we haven't proved it. In numerical experiments, we have found that some dependence does not adversely effect PCA, but a detailed examination of this matter would necessitate a longer paper.  \n\n> - Since the PCA version of the algorithm operates on fixed dimension, it is not clear to my why in (9) an increase of $p$ should bring any benefit. ...\n\nThe proof of theorem 2 and eq. (9) involves analysing the top $r$ eigenvalues and eigenvectors of the $n\\times n$ matrix whose $(i,j)$-th entry is $p^{-1}\\langle\\mathbf{Y}_i,\\mathbf{Y}_j\\rangle$. As discussed above,  $p^{-1}\\langle\\mathbf{Y}_i,\\mathbf{Y}_j\\rangle - \\alpha(Z_i,Z_j) \\to 0$ when $p\\to\\infty$. So $p\\to\\infty$ helps us prove the above mentioned eigenvalues/vectors are close to those of the matrix with elements $\\alpha(Z_i,Z_j)$.  For this we need a stronger form of convergence than element-wise matrix convergence, this is where the additional assumption of independence is used."}, {"review_id": "CduF6QvQpA", "reviewer": "Reviewer_s6Xk", "summary": "The paper studies hierarchical agglomerative clustering under a specified generative process for the underlying data vectors. The main focus of the paper is similarity based clustering that deals with inner products between vectors rather than computing pairwise distances that has been studied before.\n\nThe goal of the paper is to provide recovery guarantees for the underlying tree structure. The paper provides such guarantees under their specified generative model and they bound the maximum merge distortion  by the affinity estimation error (see Theorem 1). Moreover, they provide tradeoffs for the estimation error based on the dimension of the data vectors and the sample size. Finally, the authors test their method in terms of runtime and quality against several other methods from the literature.\n\n", "questions": "From above: \nQ1: Do we learn something interesting by this dot-product analysis and how are the dimension/sample/estimation tradeoffs different that those other papers that did distance-based analyses?\nQ2: Perhaps easy, but giving more evidence to motivate your particular hierarchical tree model would be appreciated.\nQ3: Technical question: Perhaps to help illustrate some of your techniques/ideas, it would be helpful to instantiate your model to a small-depth hierarchy. Say there are only 2 or 3 levels like you had in some of your experimental datasets (S&P500, 20Newsgroups etc). Can your theorems be more interpretable in this scenario?", "limitations": "I like the paper and I find that the authors addressing the concerns mentioned above will significantly improve their paper.\nI don't find any major limitations, apart from better motivating their model, comparing it with the existing ones, and better illustrating their Theorems and what we learn from them (in comparison to distance-based methods).", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "+hierarchical tree recovery is not a well-understood problem so the overall research direction is interesting\n+dot products are used a lot in real-world applications so having analyses based on them is important\n+the flavor of the guarantees seem to be towards the right direction for what we would mean \"hierarchical tree recovery\"\n+the technical aspects of the proofs", "weaknesses": "-The generative model should be better explained and compared to other previously studied models for hierarchical clustering. In particular, it was not clear to me how to think of your model given that Hierarchical Stochastic Block Models and well-clusterable graphs are defined and seem more natural to me, see e.g. papers below:\n\nCohen-Addad, Kanade, Mallmann-Trenn, Mathieu: \"Hierarchical clustering: Objective functions and algorithms\" \nChatziafratis, Mahdian, Ahmadian: \"Maximizing Agreements for Ranking, Clustering and Hierarchical Clustering via MAX-CUT\"\nManghiuc, Sun: \"HierarchicalClustering: O(1)-Approximation for Well-Clustered Graphs\" \n\nSo overall, having a better exposition and examples for why this particular model you proposed is a natural one, would help the readers a lot.\n\n-The difference from other papers that focus on distance-based hierarchical clustering is important for the mathematics, but conceptually the algorithm is very related. So here the innovation in terms of conceptual contribution is slightly weakened. However, from technical viewpoint there are interesting ideas the paper introduces.\n\n-In the statement of the Theorems 1 and 2, I missed some interpretaion of the results and especially some comparison to related work on distance-based methods. Do we learn something interesting by this dot-product analysis? How are the dimension/sample/estimation tradeoffs different that those other papers?\n\n-The numerical experiments: I am not sure they convey the benefits of your algorithm in terms of estimation error and recovery. In Table 1, I would clarify what are the most significant advantages of your algorithm as currently it is hard to grasp.\n\n\n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "> The generative model should be better explained and compared to other previously studied models for hierarchical clustering. ... see e.g. papers below ... having a better exposition and examples for why this particular model you proposed is a natural one, would help the readers a lot.\n\nThanks for these references. In fact we already cite Cohen-Addad et al. \"Hierarchical clustering: Objective functions and algorithms\" in the last paragraph of section 1, it is item 14 in our bibliography. The paper of Manghiuc et al. falls within the cost-function based framework of Dasgupta (bibliograpy item 15) cited and discussed in section 1. There we note that this existing modelling approach involves assuming an underlying ultrametric space whose geometry specifies the unknown tree. The paper of  Chatziafratis et al. you point to is closely related to the paper of Charikar and Chatziafratis \"Approximate Hierarchical Clustering via Sparsest Cut and Spreading Metrics\" (bibliography item 11) also cited in section 1. There we highlight that high-dimensionality $p\\to\\infty$ is a key consideration in our analysis, whereas in existing works dimension is either fixed, or not considered at all. \n\nRegarding the matter of how to think of our model, our entire setup is really quite different and new compared to existing works, so direct comparisons may not be very meaningful. However, one key feature is that the \"ground truth\" tree structure underlying our model is the conditional independence graph (see eq. (2) in the manuscript), as opposed to e.g. the geoemtry of an ultrametric as mentioned above. This is really a new perspective, and we think it is important because of the fundamental role that conditional independence plays in hierarchical statistical modelling. \n\nRegarding exposition, we intended fig. 1 to help the reader understand how this tree structure relates to data, via the specific notion of dendrogram which we introduce in the paper. Regarding examples, one way we could add some content here is to expand upon the simulation example described in appendix D, which we use in our numerical results in section 4. \n\nWe would like to clarify a somewhat subtle point about our objectives: it is _not_ our main objective to specify a model which we necessarily believe is more natural than those in other works. Rather, our intention is to uncover the fact that a fairly standard algorithm, albeit using a dot-product affinity, allows the conditional independence tree underlying data to be recovered, under assumptions which we believe are quite general. \n\n> - The difference from other papers that focus on distance-based hierarchical clustering is important for the mathematics, but conceptually the algorithm is very related. So here the innovation in terms of conceptual contribution is slightly weakened.\n\nIndeed we are not claiming significant methodological innovation, but rather a new perspective on a well-known type of algorithm.\n\n>  - In the statement of the Theorems 1 and 2, I missed some interpretaion ... especially some comparison to related work on distance-based methods. Do we learn something interesting by this dot-product analysis? How are the dimension/sample/estimation tradeoffs different that those other papers?\n\nAs noted above and in section 1 of the manuscript, in existing works, dimension $p$ is either fixed or does not appear at all, and various forms of convergence are driven by increasing sample size $n\\to\\infty$, whereas in our setup $p\\to\\infty$ drives convergence - see eq. (8) and (9). To our knowledge, this is the first time it has been established that high-dimensionality might have a beneficial effect on performance. The convergence rate $n^{2/q}/\\sqrt{p}$ in theorem 2 is unusual. At first glance it may appear strange that sample size $n$ appears in the numerator, but this is a reflection of the very conservative $\\max$-form of error we are considering on the r.h.s. of (8). Moreover, we see from $n^{2/q}$ that the polynomial moment parameter $q$  from assumption A2 ameliorates the effect of $n$ being large. These are all new insights. \n\nRegarding comparison to distance-based methods, we note that supplementary material appendix E \"Understanding agglomerative clustering with Euclidean or cosine distances in our framework\" is dedicated to helping make such connections. In E.3 we establish a convergence result (theorem 3) similar to theorem 2 but concerning cosine distance, under a modified model in which noise is multiplicative rather than additive as in eq. (1).  \n\n> - The numerical experiments: I am not sure they convey the benefits of your algorithm... . In Table 1, I would clarify what are the most significant advantages of your algorithm as currently it is hard to grasp.\n\nThanks for this suggestion, the basic message here is that our algorithm does a better job of recovering ground-truth hierarchy, except for with the S/&P500 data - shortcomings of our modelling approach which may explain why are discussed in appendix E.4.\n\n>**Questions**\n> - Q1: Do we learn something interesting by this dot-product analysis and how are the dimension/sample/estimation tradeoffs different that those other papers that did distance-based analyses?\n\nSee above.\n\n> - Q2: Perhaps easy, but giving more evidence to motivate your particular hierarchical tree model would be appreciated.\n\nAlso see above.\n\n> - Q3: Technical question: Perhaps to help illustrate some of your techniques/ideas, it would be helpful to instantiate your model to a small-depth hierarchy. ... Can your theorems be more interpretable in this scenario?\n\nThanks for this thought-provoking question. Having considered this, we think that the number of levels as you suggest does enter very directly in to our theoretical framework, and so we can't see an obvious simplification here. This may seem suprising, but we believe it is related to the fact we quantfy performance with merge heights, pair-wise between data points.\n"}, {"review_id": "EUGG21nbtW", "reviewer": "Reviewer_FShC", "summary": "The authors discuss a phylogenetic reconstruction problem (I'm not 100% clear on the exact problem, though) and suggest to use the dot product as a measure of similarity (or affinity). In particular, they seem to use the UPGMA algorithm (Alg. 1) where similarity is defined by dot product (scaled by $1/p$).\n\nThe proposed method is shown to recover the underlying tree structure with guaranteed accuracy under reasonable theoretical assumptions. Empirical results indicate good performance in a number of cases where the ground-truth hierarchy is known.", "questions": "I am familiar with phylogenetics and probabilistic graphical models, but despite reasonable efforts, I fail to be able to comprehend the modeling setup. The authors should clearly explain what the model in Eq. (1) means. In particular, the meaning of the latent variables $Z_i$, and the two mappings, $X(Z_i)$ and $S(Z_i)$, that generate the observed data from the latent variables and the random noise, $E_i$, should be explained.\n\nAnother general thing I was left wondering is the logic behind the dot product as an affinity measure. It is common to use the cosine similarity $\\cos \\theta = \\frac{\\langle x, y\\rangle}{||x|| ||y||}$, where $\\theta$ is the angle between the vectors $x$ and $y$, and $\\langle x, y\\rangle$ is the dot product between them. In other words, the cosine similarity is the dot product scaled by the product of the Euclidean norms of the two vectors. Such scaling can be motivated by noting that without normalization, the dot product has the peculiar property that, for instance, $\\langle Cx, x\\rangle = C \\langle x, x\\rangle$, i.e., for $C>1$, the vectors $Cx$ and $x$ are more similar (in terms of the dot product) to each other than the vectors $x$ and $x$ (so the vector $x$ to itself). For normalized vectors, the dot product and the cosine similarity are identical. (Further, for normalized vectors, the Euclidean distance equals $\\sqrt{2 - 2{\\langle x, y\\rangle}}$ so it too is a monotonic function of the dot product.)  I would appreciate some comments on the motivation of using the dot product as opposed to commonly used metrics such as the Euclidean distance or the cosine similarity, and on when it is likely to be appropriate (and when not).\n\ndetailed comments:\n- I suppose some would prefer to use the $x \\cdot y$ notation for dot product since $\\langle x,y \\rangle$ often denotes the (more general) inner product\n- p. 2: \"distributional properties $\\mathbf{X}$\": should this be \"... of $\\mathbf{X}$\"?\n- p. 7: When you say you compare against the UPGMA method, you should say what distance metric you use in UPGMA. In fact, it would be interesting to see the results of UPGMA (and other methods) based on various distance metrics (dot product, cosine, Euclidean, ...). And indeed, why not include more phylogenetic methods such as Neighbor-Joining (NJ), etc?\n- appendix C, proof of Lemma 2: the second line of the displayed equation should have a factor 2 in front of the dot product term (as in $||x+y||^2 = ||x||^2 + 2 \\langle x,y \\rangle + ||y||^2$)", "limitations": "See above", "rating": 7, "confidence": 3, "soundness": 3, "presentation": 2, "contribution": 3, "strengths": "- the theoretical guarantees appear very strong", "weaknesses": "- presentation should be improved so that the problem under study is clearer\n- it seems that methodologically, there is no novelty beyond proposing to use the dot product as a similarity measure in the UPGMA method\n- limited set of benchmark methods in the empirical part", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": ">**Weaknesses**\n> - presentation should be improved so that the problem under study is clearer\n\nSee response to questions below.\n\n> - it seems that methodologically, there is no novelty beyond proposing to use the dot product as a similarity measure in the UPGMA method\n\nIndeed beyond the point of using the dot-product affinity measure, we are not claiming methodological novelty. Rather our contribution is to put forward an entirely new perspective on a very well known and very widely used type of algorithm. \n\n> - limited set of benchmark methods in the empirical part\n\nPlease see below for clarifications - especially regarding UPGMA and the additional numerical results in appendix D.\n\n>**Questions**\n> - I am familiar with phylogenetics and probabilistic graphical models, but despite reasonable efforts, I fail to be able to comprehend the modeling setup. The authors should clearly explain what the model in Eq. (1) means.\n\nWe would happily add some commentary. For example,  for each $v\\in\\mathcal{Z}$ one can think of the vector $\\mathbf{X}(v)\\in\\mathbb{R}^p$ as the random centre of a \"cluster\", with correlation structure of the cluster determined by the matrix $\\mathbf{S}(v)$. The latent variable $Z_i$ indicates which cluster the ith data vector $\\mathbf{Y}_i$ is associated with. The vectors $\\mathbf{X}(v)$, $v\\in\\mathcal{V}\\setminus\\mathcal{Z}$, correspond to unobserved vertices in the underlying tree.\n\n> -Another general thing I was left wondering is the logic behind the dot product as an affinity measure.  It is common to use the cosine similarity ... I would appreciate some comments on the motivation of using the dot product as opposed to commonly used metrics such as the Euclidean distance or the cosine similarity, and on when it is likely to be appropriate (and when not).\n\nWe wonder if the reviewer may have missed appendix E in the supplementary material, \"Understanding agglomerative clustering with Euclidean or cosine distances in our framework\". We refer to appendix E on lines 46, 139, 232, 256 and 299 of the main part of the manuscript.\n\nSection E.1 discusses use of Euclidean distance as a similarity meausre and section E.2 makes the connection between our algorithm and using cosine distance. Theorem 3 in section E.3 shows that using cosine distance as a similarity measure can work well under an alternative model in which errors are multiplicative, rather than additive as in eq. (1) in the main part of the manuscript, and the \"error-free\" data vectors $\\mathbf{X}(Z_i)$ all have the same expected square magnitude. Contrasting this with theorem 2 highlights that the dot-product affinity is a good choice to help remove additive noise, and can cope with data where expected square magnitudes vary across data points - see also figure 1 of the main manuscript for interpretation of dendrogram height.\n\nIn appendix E.4 we discuss the limitations of our modelling assumptions and failings of the dot-product affinities.\n\nIn connection with what the reviewer writes about $\\langle Cx,x\\rangle,$ we had thoughts along the same lines: we chose the word dot-product \"affinity\" rather than \"similarity\" purposefully;  indeed it is possible that $\\langle Cx,x\\rangle >\\langle x,x \\rangle$, whereas it would be confusing to suggest that $Cx$ and $x$ are more \"similar\" to each other than $x$ is to itself. We would happily add a note to explain this in the manuscript.\n\n> - I suppose some would prefer to use the $\\cdot$ notation for dot product since often denotes the (more general) inner product.\n\nWe see your point here and have experimented with other notation, but found the \"$\\cdot$\"\" notation to be less visually clear in many equations.\n\n> - p. 2: \"distributional properties \": should this be \"... of \"?\n\nThanks for catching this.\n\n> - p. 7: When you say you compare against the UPGMA method, you should say what distance metric you use in UPGMA. In fact, it would be interesting to see the results of UPGMA (and other methods) based on various distance metrics (dot product, cosine, Euclidean, ...). And indeed, why not include more phylogenetic methods such as Neighbor-Joining (NJ), etc?\n\nWe had used \"UPGMA\" to refer specifically to the combination of average linkage and Euclidean distance, in retrospect we realise this is confusing. In fact results for average linkage with cosine distance (i.e. UPGMA with cosine distance) are already given in table 1 in the main manuscript (there labelled \"Cosine distance\"). We could clarify this. \n\nFollowing the reviewers suggestion, in the .pdf attached to this response we give numerical results for UPGMA with Manhattan distance, and UPGMA with dot-product as a dissimilarity measure (hence \"opposite\" to our Algorithm 1). Performance of these algorithms is worse than that of Algorithm 1. \n\nWe wonder if the reviewer may not have seen the additional numerical results in table 2 in appendix D of the supplementary material, where we compare our algorithm against other combinations of cosine and Euclidean distances, with complete and single linkage functions. In total across the main manuscript and appendix D there are comparisons of our method against 8 others:\n\nmain manuscript, table 1:\n- \"Cosine distance\" (i.e. UPGMA with cosine distance)\n- HDBSCAN\n- \"UPGMA\" (i.e. UPGMA with Euclidean distance)\n- Ward's method\n\nappendix D, table 2:\n- Complete linkage with Euclidean distance\n- Complete linkage with Cosine distance\n- Single linkage with Euclidean distance\n- Single linkage with Cosine distance\n\nso including the new UPGMA results, the total number of algorithms we compare against is 10.\n\nThanks for the suggestion about neighbour-joining. Our theory and performance measure concern hierarchy recovery, and in the time available we haven't found an implementation of neighbour-joining which allows hierarchy (rather than a \"flat\" clustering) to be extracted. With more time we could investigate this further.\n\n> - appendix C, proof of Lemma 2: the second line...\nThanks for catching this."}, {"review_id": "YYaFnecscV", "reviewer": "Reviewer_wZ8j", "summary": "This paper presents new analysis and perspectives on one of the most widely used clustering algorithms, hierarchical agglomerative clustering (HAC). In particular, the authors consider the relationship between a particular generative process of data and a dot-product based linkage of HAC. Empirical and theoretical results are presented.", "questions": "* Q1. How could the presentation be modified to address W1?\n* Q2. (W2) How does reducibility related to the proofs of Lemma 4 / Prop 1?\n* Q3. (W3) did you consider using any other metrics for evaluation? If you were to add another metric to provide a new slice of information, which metrics/measurements would you want to add?\n* Q4. Please see W4\n\n", "limitations": "Yes.", "rating": 5, "confidence": 3, "soundness": 3, "presentation": 2, "contribution": 3, "strengths": "This paper presents an interesting perspective on a widely studied clustering algorithm. \n\nI don't think that the intention of the paper is to yield surprise that for certain kinds of data there might be one linkage function (dot product) that works better than others. Rather I think the intention of the paper is to draw connections between clustering models and HAC.\n\nIn particular, I see the strengths as: \n\n* Connections between the model for latent tree structures in Eq. 1 and Theorem 2 and dot-product based average linkage\n* Connections between models for trees with heights and HAC.\n* Empirical analysis of some of the theoretical scaling properties: Figure 2.\n", "weaknesses": "I think that the paper could be improved in the following ways: \n\n* W1. I think that more explicit treatment of Eq. (2) and (3) would improve the presentation; e.g., showing where/why these hold under the model in Eq. (1). \n* W2. My understanding is that Lemma 4 & Prop 1 looks quite similar to standard HAC proofs about reducibility? I think it would greatly improvement the treatment of the result to explain the differences / distinctions. Apologies if I have missed something.\n* W3. While the evaluation metric seems quite reasonable, might also be interesting to show results across other metrics such as dendrogram purity or similar metrics. I see this as a minor point though.\n* W4. In the advent of distributed (typically relatively low) dimensional representations from deep neural encoder models. I wonder if the authors could provide perspective on \"theorem 2 says that affinity estimation error vanishes if the dimension p grows faster enough relative to n\". E.g., should we think of this method as appropriate in such circumstances where p is likely in range of 128-1024 or so and n might be in the millions or billions? ", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "> **Weaknesses**\n>> - W1. I think that more explicit treatment of Eq. (2) and (3) would improve the presentation; e.g., showing where/why these hold under the model in Eq. (1).\n\n(2) and (3) are not a consequence of (1), but rather distributional assumptions we make about the ingredients of (1), we would happily clarify this in the manuscript. We could make (2) and (3) numbered assumptions, of the form (Ax), if that would be clearer.\n\n\n>> - W2. My understanding is that Lemma 4 & Prop 1 looks quite similar to standard HAC proofs about reducibility? I think it would greatly improvement the treatment of the result to explain the differences / distinctions. Apologies if I have missed something.\n\nThe arguments involved in the proofs of Lemma 4 & Prop 1 are indeed similar in flavour to those concerning reducible linkage functions in e.g.:\n\nSumengen, Baris, et al. \"Scaling hierarchical agglomerative clustering to billion-sized datasets.\" arXiv preprint arXiv:2105.11653 (2021).\n\nand historical references therein. However we are not aware of any results in the literature which give us precisely what we need for Lemma 4 and Prop 1, and indeed we derived them from scratch. We could of course clarify this in the manuscript with appropriate references.\n\n>> - W3. While the evaluation metric seems quite reasonable, might also be interesting to show results across other metrics such as dendrogram purity or similar metrics. I see this as a minor point though.\n\nThanks for this suggestion. Looking up dendrogram purity in:\n\nHeller, Katherine A., and Zoubin Ghahramani. \"Bayesian hierarchical clustering.\" Proceedings of the 22nd international conference on Machine learning. 2005.\n\nit appears to be applicable in cases where the ground-truth labels specify a partition of the data points, i.e. a \"flat\" clustering, but are not hierachical in nature. By contrast, the tau-B correlation measure in our manuscript is constructred to compare estimated with ground-truth hierarchical labelling, and indeed our theoretical results concern the quality of estimated hierarchy.  We agree it could be interesting to consider other metrics, but the connection of such metrics to our theory might be quite unclear if they don't directly quantify hierarchy recovery. This could be an avenue for future investigation.\n\n>> - W4. In the advent of distributed (typically relatively low) dimensional representations from deep neural encoder models. I wonder if the authors could provide perspective on \"theorem 2 says that affinity estimation error vanishes if the dimension p grows faster enough relative to n\". E.g., should we think of this method as appropriate in such circumstances where p is likely in range of 128-1024 or so and n might be in the millions or billions?\n\nThe polynomial moment parameter $q$, appearing and the convergence rate $n^{2/q} / p^{1/2}$ in theorem 2 is key here. Under assumption A2, $q$ quantifies how light/heavy the tails of the distributions of the data are -- higher $q$ corresponding to lighter tails. In the numerator of $n^{2/q} / p^{1/2}$ we see increasing $q$ acts to ameliorate the effect of increasing $n$. Futhermore, as mentioned on line 272 in section 4.2, in line with our empirical findings (fig. 2) we conjecture that when A2 is strengthened from polynomial to exponential-of-quadratic (i.e. sub-Gaussian) moments, the convergence rate improves from $n^{2/q} / p^{1/2}$ to $(\\log_e  n)^{1/2}  / p^{1/2}$. This means that, theoretically, convergence would occur as long as $p$ grows faster than $\\log_e n$, in which case $p$ in the range $128-1024$ you mention could be enough to counter-act $n$ in the range $10^6 - 10^9$.\n\n>**Questions:**\n>> - Q1. How could the presentation be modified to address W1?\n\nSee above.\n\n>> - Q2. (W2) How does reducibility related to the proofs of Lemma 4 / Prop 1?\n\nSee above.\n\n>> - Q3. (W3) did you consider using any other metrics for evaluation? If you were to add another metric to provide a new slice of information, which metrics/measurements would you want to add?\n\nFurther to our response on this topic above, to illustrate the appliability of our theory in greater detail, it would be desirable to have ground-truth merge heights available, so that estimated versus ground-truth merge heights could be used as a metric. Such ground-truth merge heights might potentially be available in some phylogenetic application domains, and in future we hope to engage with experts in such domains to investigate this further.  If in future we could extend out theoretical results from merge-height estimation to quantifying the quality of \"flat\" clusterings derived from those merge heights, then the dendrogram purity measure you suggest could be an interesting addition.\n\n>> - Q4. Please see W4\n\nSee above."}, {"review_id": "hJukoIV8hz", "reviewer": "Reviewer_Emfa", "summary": "For agglomerative clustering, this paper proposes a method where the affinity of clusters at each hierarchy is naturally visualized as their height. Section 2 describes its statistical model (a nonparametric model based on first order moments) and the clustering algorithm derived from it (Algorithm 1) is presented. Algorithm 1 is extremely simple (needless to say, this simplicity is a virtue of the proposed method), and is only a slight modification of many of the standard agglomerative clustering methods. However, despite its simplicity, the deep insights of the algorithm are explained in Section 3. Remarkably, very roughly speaking, the affinity of each hierarchical cluster is expressed as the height of its dendrogram, which, under certain assumptions, asymptotically approaches its true value. Finally, Section 4 describes the effectiveness of the proposed algorithm on real data, as well as its limitations, and provides code for follow-up studies. \n\nNote: I was not able to understand the details (techniques in the supplementary material) of some of the theoretical analysis in this paper within the initial review period. For Theorem 1, I was able to follow the outline of its key mathematical induction. For Theorem 2, I have not been able to grasp the details. However, the statements in the text are very reliable, and the authors have made an effort to add intuitive explanations to these theoretical results.", "questions": "I am very grateful to authors for sharing very insightful and interesting ideas. I have some questions because this paper contains simple methodology but very surprising (and perhaps intuitively non-trivial) results. These questions are the reason I am not giving this paper a higher score at this time for its overall rating.\nThese may simply be my lack of insight, but perhaps a brief mention in the text or in the supplementary material would make the advantages of the proposed method clearer to some readers.\nI would be very grateful if a response from the authors could clear up some of my questions.\n\n(1) Another aspect of the model assumptions (Section 2.1).\n\nIt is very interesting to note that the only model mentioned in this paper is the chain rule by Equation 3, and everything else can be expressed nonparametrically. However, does Equation 3 imply that \"all data are i.i.d. samples from a single cluster (a single nonparametric distribution that assumes only first order moments)\"? In other words, is there model misspecification occurring for data sets with multiple clusters (e.g., a 2-cluster Gaussian mixture model)? If we use the chaining rule in Equation 3 to marginalize the random variables at each node in turn from the terminal node to the root node of the tree structure, this may be a model in which all data are sampled i.i.d. from some nonparametric distribution with identical mean values (only the first order moments are specified). Is this intuition wrong? If this intuition is correct, then surely the \"PCA does not do anything wrong\" in the second half of Theorem 2 (Equation 9) makes a lot of sense. Conversely, if the proposed statistical model is such that it has multiple clusters (like a Gaussian mixture model, where each cluster tends to gather data at its center), then it is highly counterintuitive that \"still PCA does not do anything bad\". In summary, my question comes down to whether the proposed model assumes a one-cluster i.i.d. model.\n\n(2) The validity of greedily searching for the closest affinity pair in row 3 of Algorithm 1. \n\nFrom the viewpoint of an agglomerative algorithm, it seems quite reasonable to greedily look for the closest affinity pair in line 3 of the algorithm. On the other hand, from the perspective of a statistical model, a probabilistic selection (i.e., Markov chain selection), such as finding a single pair from a categorical distribution weighted by affinity proportions, also seems natural. Does the fact that this row 3 is a greedy selection play an important role in the theoretical analysis in section 3? The actual intent of this question can be attributed to something that is also closely related to question (1), as follows.\nFor example, if the true data distribution is like a Gaussian mixture model with multiple clusters, is it possible to recover the true data distribution (mixture model structure) by, for example, truncating the tree at some hierarchy (resolution) of the tree structure obtained with the proposed model? In other words, would there not be a misspecification problem as a statistical model? My rough intuition tells me that it is not easy to recover the structure of a statistical model in the case of a greedy merge.\nI also know that such a topic is a bit off the author's actual intent for this paper. However, the question arises because I am very curious as to why the statistical model behind the proposed method leads to very good properties (Theorem 1 and Theorem 2) in Algorithm 1.\n\n(3) (Very minor.) Reason why DAG is first introduced. \n\nIs there a reason why when the tree model is introduced in lines 68-69, it is broadly defined as directed acyclic graphs (DAGs) rather than restricted to binary trees?\nDAGs are often used as the structure of Bayesian networks. Are there any indications that the theoretical results of this study for binary trees may be extended to DAGs or Bayesian networks?\n\nFinally, once again, I appreciate your sharing of some very interesting ideas.", "limitations": "The paper carefully questions the weaknesses of the proposed method in Section 5. It is well explained what kind of data and in what cases the proposed method is unlikely to show its true value.\n\nOn another point, I am slightly concerned whether the statistical model in the proposed method implicitly assumes i.i.d. data with a single cluster, as I have asked the authors in my question. However, this concern may be addressed by a response from the authors.", "rating": 7, "confidence": 3, "soundness": 4, "presentation": 4, "contribution": 3, "strengths": "- This paper is an excellent combination of a simple methodology and deep insights behind it. Needless to say, methodological simplicity is a virtue of the proposed method at this time.\n- The phenomenon of affinities between clusters at each level of hierarchy appearing visualized as heights on the dendrogram will be the attractive results for a wide range of readers, both theoretically and in terms of application.\n- It is excellent that the theoretical results (Theorem 1 and Theorem 2) are carefully interpreted in Section 3.3 as detailed intuitions.\n- Authors also share the code (Jupyter Notebook) for re-experimentation with the community, which is a great contribution to subsequent researchers.\n- In Section 5, the paper explicitly mentions the slightly stronger assumption of the proposed model and clarifies for what cases (what data) it has a negative impact. This is an essential finding in the development of science and technology.", "weaknesses": "I have not been able to find any convincing weaknesses for this paper in the initial peer review. However, I have some concern about whether the statistical model assumed by the proposed method is causing model speciﬁcation for data with so-called multiple clusters. I have inquired about this concern with the author in the following question. If the concern is resolved, I have not found any major weakness in this paper.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "> **Weaknesses**\n>> I have not been able to find any convincing weaknesses for this paper in the initial peer review. However, I have some concern about whether the statistical model assumed by the proposed method is causing model speciﬁcation for data with so-called multiple clusters. I have inquired about this concern with the author in the following question. If the concern is resolved, I have not found any major weakness in this paper.\n\nPlease see below. \n\n> **Questions:**\n>> - Another aspect of the model assumptions (Section 2.1). ... It is very interesting to note that the only model mentioned in this paper is the chain rule by Equation 3, and everything else can be expressed nonparametrically. However, does Equation 3 imply that \"all data are i.i.d. samples from a single cluster (a single nonparametric distribution that assumes only first order moments)\"? In other words, is there model misspecification occurring for data sets with multiple clusters (e.g., a 2-cluster Gaussian mixture model)? \n\nThanks for this thought-provoking question. \n\nWe can cast a 2-cluster Gaussian mixture model as a special case of our model, but subject to the Gaussian cluster centers being random and themselves having the same mean. To set this up we need a tree with three vertices $\\mathcal{V}=\\\\{0,1,2\\\\}$ with $0$ being the root, edges $(0,1)$ and $(0,2)$, and $\\mathcal{Z}=\\\\{1,2\\\\}$. For simplicity, let $\\mathbf{X}(0)$ be a constant (that is, not random), and let $\\mathbf{X}(1)$ and $\\mathbf{X}(2)$ be iid random vectors, each with mean $\\mathbf{X}(0)$. If the vectors $\\mathbf{E}_i$, for $=1,\\ldots,n$ are Gaussian, then overall $Y_1,\\ldots,Y_n$ are draws from a 2-component, Gaussian mixture with random centres, as claimed.\n\nWe need $\\mathbf{X}(1)$ and $\\mathbf{X}(2)$ to have mean $\\mathbf{X}(0)$ in order for eq (2) to hold, which is in turn needed to establish the relationship between affinity $\\alpha$ and merge height $m$ in Lemma 1. This is really fundamental to our approach.\n\nSo, mathematically speaking, if data were generated from a model with deterministic cluster centers, then technically our model of the latent cluster centers would be misspecified. However, if we are only given data $Y_1,\\ldots,Y_n$ we believe it is not possible to identify whether cluster centers are random or not, i.e. this technical model misspecification for latent cluster centers does not imply misspecification of the induced model of observed data. Hence, overall we think there is no problem here.\n\n>> - (2) The validity of greedily searching for the closest affinity pair in row 3 of Algorithm 1. From the viewpoint of an agglomerative algorithm, it seems quite reasonable to greedily look for the closest affinity pair in line 3 of the algorithm. On the other hand, from the perspective of a statistical model, a probabilistic selection (i.e., Markov chain selection), such as finding a single pair from a categorical distribution weighted by affinity proportions, also seems natural. ... The actual intent of this question can be attributed to something that is also closely related to question (1), as follows. For example, if the true data distribution is like a Gaussian mixture model with multiple clusters, is it possible to recover the true data distribution (mixture model structure) by, for example, truncating the tree at some hierarchy (resolution) of the tree structure obtained with the proposed model? In other words, would there not be a misspecification problem as a statistical model? \n\nThe topic of recovering the data distribution is a fascinating one, but not one we have considered before. We note in passing, it is true that if one knows all the pairwise merge heights of a dendrogram (in out setup) then one can reconstruct the entire tree - this is discussed in appendix 3.C. What you suggest about an alternative to the greedy merge sounds somewhat like the integration performed in the EM algorithm for mixture models, or Bayesian approaches to hierarchical clustering, where the associations between data points and clusters are integrated out. We have not attempted to analyse such approaches in our framework, but this could be an interesting  and new avenue of future research -- perhaps some kind of consistency result for likelihood-based or Bayesian model fitting -- but this is well beyond the scope of the present paper.\n\n>> - (3) (Very minor.) Reason why DAG is first introduced. Is there a reason why when the tree model is introduced in lines 68-69, it is broadly defined as directed acyclic graphs (DAGs) rather than restricted to binary trees? DAGs are often used as the structure of Bayesian networks. Are there any indications that the theoretical results of this study for binary trees may be extended to DAGs or Bayesian networks?\n\nThe point of not restricting the binary trees is just to emphasise that our theory tells us the algorithm can produce useful output, without assuming such binary structure, even though the algorithm itself involves binary merges. We haven't yet considered more general Bayesian networks, but again this could be an interesting and new avenue for future research.\n\nThanks very much for these stimulating questions!\n"}]

Author rebuttal:
The authors are very grateful for the effort which the reviewers have put in to reading our manuscript and providing feedback.

We are pleased to see that, on the whole, the reviewers have recognised and engaged with the new perspective on hierarchical clustering and tree recovery which we report in this paper. Since this perspective is novel in various ways (e.g., in terms of our probabilistic model-based formulation involving conditional independence, emphasis on dot-product affinities, and theoretical investigation of high-dimensional behaviour) it necessitates the use of some concepts and notation which are non-standard in the literature on hierarchical clustering. The authors are very grateful for the various suggestions made by the reviewers which will help us improve the presentation in this regard. 

We would like to highlight appendices D and E in the supplementary material of the original submission, in case they were missed. Appendix D.6 reports additional comparative numerical results. Appendix E, entitled "Understanding agglomerative clustering with Euclidean or cosine distances in our framework", provides methodological and theoretical perspectives to help make connections between our framework and agglomerative clustering using more standard distance measures techniques.  

Attached to this rebuttal is a .pdf containing additional numerical results in response to the comments from reviewers dt3f and FShC.

Meta-review: This paper studies hierarchical clustering via the probabilistic graphical model. The authors propose a simple agglomerative hierarchical clustering algorithm, where clusters are merged by the maximum average dot product, and show its interesting theoretical properties based on the formulation. 

This paper is overall well-written, and all the reviewers are happy with the contribution of the paper and recommend acceptance. I also agree with it and like the simplicity of the algorithm with its interesting theoretical analysis.
In particular, this paper gives a new insight into hierarchical clustering using the probabilistic graphical model, which is inspiring and can trigger further development of theoretical studies of hierarchical clustering. 

Since the authors made a significant effort on rebuttals, which contain relevant information and successfully address concerns raised by the reviewers, I strongly recommend integrating such contents in the camera-ready version of the paper.

To summarize, I would appreciate the contribution of the paper, which is sound and innovative. Therefore, I recommend the spotlight presentation of the paper.
Decision: Accept (spotlight)