Ultrametric Fitting by Gradient Descent

TL;DR

This paper introduces a novel gradient descent-based optimization framework for fitting ultrametrics to data, enabling flexible hierarchical clustering with improved performance over traditional methods.

Contribution

It presents a new unconstrained optimization approach for ultrametric fitting, allowing the use of various cost functions and improving upon standard hierarchical clustering techniques.

Findings

Effective ultrametric fitting on synthetic datasets

Competitive results against state-of-the-art algorithms

Flexible framework for hierarchical clustering

Abstract

We study the problem of fitting an ultrametric distance to a dissimilarity graph in the context of hierarchical cluster analysis. Standard hierarchical clustering methods are specified procedurally, rather than in terms of the cost function to be optimized. We aim to overcome this limitation by presenting a general optimization framework for ultrametric fitting. Our approach consists of modeling the latter as a constrained optimization problem over the continuous space of ultrametrics. So doing, we can leverage the simple, yet effective, idea of replacing the ultrametric constraint with a min-max operation injected directly into the cost function. The proposed reformulation leads to an unconstrained optimization problem that can be efficiently solved by gradient descent methods. The flexibility of our framework allows us to investigate several cost functions, following the classic paradigm of combining a data fidelity term with a regularization. While we provide no theoretical guarantee to find the global optimum, the numerical results obtained over a number of synthetic and real datasets demonstrate the good performance of our approach with respect to state-of-the-art agglomerative algorithms. This makes us believe that the proposed framework sheds new light on the way to design a new generation of hierarchical clustering methods. Our code is made publicly available at \url{https://github.com/PerretB/ultrametric-fitting}.

Equations

$$(\forall (x,y,z)\in V^3)d(x,y)\le \max \left\{d(x,z),d(z,y)\right\}.$$

$$\text{u∈W is an ultrametric on G}\iff (\forall C\in \mathcal{C},\forall e\in C)u(e)\le \underset{e^{\prime }\in C\backslash \left\{e\right\}}{\max }u(e^{\prime }).$$

$$(\forall (x,y)\in V^2)d_u(x,y)=\underset{P\in {\mathcal{P}}_{xy}}{\min }\underset{e^{\prime }\in P}{\max }u(e^{\prime }),$$

$$\bordermatrix{&x_{1}&x_{2}&x_{3}&x_{4}\cr x_{1}&0&r_{1}&r_{3}&r_{3}\cr x_{2}&r_{1}&0&r_{3}&r_{3}\cr x_{3}&r_{3}&r_{3}&0&r_{2}\cr x_{4}&r_{3}&r_{3}&r_{2}&0\cr}$$

$$\underset{u\in \mathcal{W}}{\mathrm{minimize}}J(u;w)\text{s.t.}\text{u is an ultrametric on G}.$$

$$(\forall \tilde{w}\in \mathcal{W},\forall e_{xy}\in E){\Phi }_{\mathcal{G}}(\tilde{w})(e_{xy})=\underset{P\in {\mathcal{P}}_{xy}}{\min }\underset{e^{\prime }\in P}{\max }\tilde{w}(e^{\prime }),$$

$$\operatorname*{minimize}_{\tilde{w}\in\mathcal{W}}\;J\big{(}\Phi_{\mathcal{G}}(\tilde{w});\,w\big{)}.$$

$$J_{\rm closest}(u;w)=\sum_{e\in E}\Big{(}u(e)-w(e)\Big{)}^{2}.$$

$$J_{\mathrm{size}}(u)=\sum _{e_{xy}\in E}\frac{u(e_{xy})}{{\gamma }_u({\text{lca}}_u(x,y))}.$$

$$(\forall n\in T_u){\gamma }_u(n)=\min \left\{|c|,c\in \mathrm{Childre}{\mathrm{n}}_{\mathrm{u}}(\mathrm{n})\right\},$$

$$\mathcal{T}=\big{\{}({\rm ref},{\rm pos},{\rm neg})\in V^{3}\;|\;\mathcal{L}_{\rm ref}=\mathcal{L}_{\rm pos}\quad\textrm{and}\quad\mathcal{L}_{\rm ref}\neq\mathcal{L}_{\rm neg}\big{\}}.$$

$$J_{\mathrm{triplet}}(u)=\sum _{(\mathrm{ref},\mathrm{pos},\mathrm{neg})\in \mathcal{T}}\max \{0,\alpha +d_u(\mathrm{ref},\mathrm{pos})-d_u(\mathrm{ref},\mathrm{neg})\}.$$

$$J_{\mathrm{Dasgupta}}(u;w)=\sum _{e_{xy}\in E}\frac{|{\text{lca}}_u(x,y)|}{w(e_{x,y})}.$$

$$|n|=\sum _{y\in V}H({\text{alt}}_u(n)-d_u(x,y)),$$

$$\bordermatrix{&e_{13}&e_{12}&e_{23}&e_{35}&e_{25}&e_{24}&e_{45}\cr e_{13}&0&0&0&0&0&0&0\cr e_{12}&1&1&0&0&0&0&0\cr e_{23}&0&0&1&0&0&0&0\cr e_{35}&0&0&0&0&0&0&0\cr e_{25}&0&0&0&0&0&0&0\cr e_{24}&0&0&0&1&1&1&0\cr e_{45}&0&0&0&0&0&0&1}$$

$$(\forall \tilde{w}\in \mathcal{W},\forall e_{xy}\in E){\Phi }_{\mathcal{G}}(\tilde{w})(e_{xy})=\tilde{w}(e_{xy}^{\mathrm{mst}})\mathrm{with}{e}_{xy}^{\mathrm{mst}}=\sigma ({\text{lca}}_{\text{SL}(\tilde{w})}(x,y)).$$

$$\frac{\partial \Phi (\tilde{w})}{\partial \tilde{w}}=\left[\text{mathbbm}{1}_{{\Phi }_{\mathcal{G}}^{\ast }({\tilde{w}}_1)},\dots ,\text{mathbbm}{1}_{{\Phi }_{\mathcal{G}}^{\ast }({\tilde{w}}_M)}\right],$$

$$J_{\rm triplet}\big{(}{u}\big{)}=\sum_{({\rm ref},{\rm pos},{\rm neg})\in\mathcal{T}}\max\{0,\alpha+{u}\big{(}\sigma(\textrm{lca}_{u}({\rm ref},{\rm pos}))\big{)}-{u}\big{(}\sigma(\textrm{lca}_{u}({\rm ref},{\rm neg}))\big{)}\}.$$

$$\text{card}(n)$$

$$\displaystyle=\frac{1}{2}\sum_{x\in\sigma(n)}\Big{(}\ell(\textrm{alt}_{u}(n)-d_{u}(x,x))+\sum_{y\in V\backslash{\{x\}}}\ell(\textrm{alt}_{u}(n)-d_{u}(x,y))\Big{)}$$

$$\displaystyle=\frac{1}{2}\sum_{x\in\sigma(n)}\Big{(}\ell(\textrm{alt}_{u}(n))+\sum_{y\in\mathcal{A}(x)}\sum_{z\in c_{\hat{x}}(y)}\ell(\textrm{alt}_{u}(n)-d_{u}(x,z))\Big{)}$$

$$\displaystyle=\frac{1}{2}\sum_{x\in\sigma(n)}\Big{(}\ell(\textrm{alt}_{u}(n))+\sum_{y\in\mathcal{A}(x)}\sum_{z\in c_{\hat{x}}(y)}\ell(\textrm{alt}_{u}(n)-\textrm{alt}_{u}(y))\Big{)}$$

$$\displaystyle=\frac{1}{2}\sum_{x\in\sigma(n)}\Big{(}\ell(\textrm{alt}_{u}(n))+\sum_{y\in\mathcal{A}(x)}|c_{\hat{x}}(y)|\ell(\textrm{alt}_{u}(n)-\textrm{alt}_{u}(y))\Big{)}$$

$$\displaystyle=\frac{1}{2}\sum_{x\in\sigma(n)}\Big{(}\ell\big{(}{u}(\sigma(n))\big{)}+\sum_{y\in\mathcal{A}(x)}|c_{\hat{x}}(y)|\ell\big{(}{u}(\sigma(n))-{u}(\sigma(y))\big{)}\Big{)},$$

$$(\forall e\in E)u(e)=\sum _{k=1}^K{r}_k{\mathcal{L}}_k(e),$$

$$\widetilde{J}_{\rm closest}(\mathcal{L},r;w)=\sum_{e\in E}\sum_{k=1}^{K}\mathcal{L}_{k}(e)\big{(}r_{k}-w(e)\big{)}^{2},$$

$${\bar{r}}_k=\frac{{\sum }_{e\in E}{\bar{\mathcal{L}}}_k(e)w(e)}{{\sum }_{e\in E}{\bar{\mathcal{L}}}_k(e)}.$$

Full text

1 Introduction

Ultrametrics provide a natural way to describe a recursive partitioning of data into increasingly finer clusters, also known as hierarchical clustering [1]. Ultrametrics are intuitively represented by dendrograms, i.e., rooted trees whose leaves correspond to data points, and whose internal nodes represent the clusters of its descendant leaves. In topology, this corresponds to a metric space in which the usual triangle inequality is strengthened by the ultrametric inequality, so that every triple of points forms an isosceles triangle, with the two equal sides at least as long as the third side. The main question investigated in this article is: << How well can we construct an ultrametric to fit the given dissimilarity data? >> This is what we refer to as ultrametric fitting.

Ultrametric fitting can be traced back to the early work on numerical taxonomy [2] in the context of phylogenetics [3]. Several well-known algorithms originated in this field, such as single linkage [4], average linkage [5], and Ward method [6]. Nowadays, there exists a large literature on ultrametric fitting, which can be roughly divided in four categories: agglomerative and divisive greedy heuristics [7, 8, 9, 10, 11, 12, 13], integer linear programming [14, 15, 16], continuous relaxations [17, 18, 19, 20], and probabilistic formulations [21, 22, 23]. Our work belongs to the family of continuous relaxations.

The most popular methods for ultrametric fitting probably belong to the family of agglomerative heuristics. They follow a bottom-up approach, in which the given dissimilarity data are sequentially merged through some specific strategy. But since the latter is specified procedurally, it is usually hard to understand the objective function being optimized. In this regard, several recent works [9, 16, 19, 11, 24] underlined the importance to cast ultrametric fitting as an optimization problem with a well-defined cost function, so as to better understand how the ultrametric is built.

Recently, Dasgupta [9] introduced a cost function for evaluating an ultrametric, and proposed an heuristic to approximate its optimal solution. The factor of this approximation was later improved by several works, based on a linear programming relaxation [16], a semidefinite programming relaxation [19], or a recursive $\phi$-sparsest cut algorithm [24]. Along similar lines, it was shown that average linkage provides a good approximation of the optimal solution to Dasgupta cost function [25, 26]. Closer to our approach, a differentiable relaxation inspired by Dasgupta cost function was also proposed [20]. Moreover, a regularization for Dasgupta cost function was formulated in the context of semi-supervised clustering [23, 27], based on triplet constraints provided by the user.

More generally, the problem of finding the closest ultrametric to dissimilarity data was extensively studied through linear programming relaxations [18] and integer linear programming [15]. A special case of interest arises when the dissimilarities are specified by a planar graph, which is a natural occurrence in image segmentation [28, 29, 30, 31]. By exploiting the planarity of the input graph, a tight linear programming relaxation can be derived from the minimum-weight multi-cut problem [14]. There exist many other continuous relaxations of discrete problems in the specific context of image segmentation [32, 33, 34, 35, 36, 37], but they typically aim at a flat representation of data, rather than hierarchical.

Contribution. We propose a general optimization framework for ultrametric fitting based on gradient descent. Our approach consists of optimizing a cost function over the continuous space of ultrametrics, where the ultrametricity constraint is implicitly enforced by a min-max operation. We demonstrate the versatility of our approach by investigating several cost functions:

the closest-ultrametric fidelity term, which expresses that the fitted ultrametric should be close to the given dissimilarity graph; 2. 2.

the cluster-size regularization, which penalizes the presence of small clusters in the upper levels of the associated hierarchical clustering; 3. 3.

the triplet regularization for semi-supervised learning, which aims to minimize the intra-class distance and maximize the inter-class distance; 4. 4.

the Dasgupta fidelity term, which is a continuous relaxation of Dasgupta cost function expressing that the fitted ultrametric should associate large dissimilarities to large clusters.

We devise efficient algorithms with automatic differentiation in mind, and we show that they scale up to millions of vertices on sparse graphs. Finally, we evaluate the proposed cost functions on synthetic and real datasets, and we show that they perform as good as Ward method and semi-supervised SVM.

2 Ultrametric fitting

Central to this work is the notion of ultrametric, a special kind of metric that is equivalent to hierarchical clustering [38]. Formally, an ultrametric $d\colon V\times V\to\mathbb{R}_{+}$ is a metric on a space $V$ in which the triangle inequality is strengthen by the ultrametric inequality, defined as

[TABLE]

The notion of ultrametric can also be defined on a connected (non-complete) graph $\mathcal{G}=(V,E)$ with non-negative edge weights $w\in\mathcal{W}$, where $\mathcal{W}$ denotes the space of functions from $E$ to $\mathbb{R}_{+}$. In this case, the distance is only available between the pairs of vertices in $E$, and the ultrametric constraint must be defined over the set of cycles $\mathcal{C}$ of $\mathcal{G}$ as follows:

[TABLE]

Note that an ultrametric $u$ on $\mathcal{G}$ can be extended to all the pairs of vertices in $V$ through the min-max distance on $u$, which is defined as

[TABLE]

where $\mathcal{P}_{xy}$ denotes the set of all paths between the vertices $x$ and $y$ of $\mathcal{G}$. This observation allows us to compactly represent ultrametrics as weight functions $u\in\mathcal{W}$ on sparse graphs $\mathcal{G}$, instead of more costly pairwise distances. Figure 1 shows an example of ultrametric and its possible representations.

Notation.

The dendrogram associated to an ultrametric $u$ on $\mathcal{G}$ is denoted by $T_{u}$ [38]. It is a rooted tree whose leaves are the elements of $V$. Each tree node $n\in T_{u}$ is the set composed by all the leaves of the sub-tree rooted in $n$. The altitude of a node $n$, denoted by $\textrm{alt}_{u}(n)$, is the maximal distance between any two elements of $n$: i.e., $\textrm{alt}_{u}(n)=\max\left\{u(e_{xy})\mid x,y\in n\textrm{ and }e_{xy}\in E\right\}$. The size of a node $n$, denoted by $|n|$, is the number of leaves contained in $n$. For any two leaves $x$ and $y$, the lowest common ancestor of $x$ and $y$, denoted $\textrm{lca}_{u}(x,y)$, is the smallest node of $T_{u}$ containing both $x$ and $y$.

2.1 Optimization framework

Our goal is to find the ultrametric that "best" represents the given edge-weighted graph. We propose to formulate this task as a constrained optimization problem involving an appropriate cost function $J\colon\mathcal{W}\to\mathbb{R}$ defined on the (continuous) space of distances $\mathcal{W}$, leading to

[TABLE]

The ultrametricity constraint is highly nonconvex and cannot be efficiently tackled with standard optimization algorithms. We circumvent this issue by replacing the constraint with an operation injected directly into the cost function. The idea is that the ultrametricity constraint can be enforced implicitly through the operation that computes the subdominant ultrametric, defined as the largest ultrametric below the given dissimilarity function. One way to compute the subdominant ultrametric is through the min-max operator $\Phi_{\mathcal{G}}\colon\mathcal{W}\to\mathcal{W}$ defined by

[TABLE]

where $\mathcal{P}_{xy}$ is defined as in (3). Then, Problem (4) can be rewritten as

[TABLE]

Since the min-max operator is sub-differentiable (see (15) in Section 3), the above problem can be optimized by gradient descent, as long as $J$ is sub-differentiable. This allows us to devise Algorithm 1.

Note that the mix-max operator already proved useful in image segmentation to define a structured loss function for end-to-end supervised learning [28, 31]. The goal was however to find a flat segmentation rather than a hierarchical one. To the best of our knowledge, we are the first to use the mix-max operator within an optimization framework for ultrametric fitting.

2.2 Closest ultrametric

A natural goal for ultrametric fitting is to find the closest ultrametric to the given dissimilarity graph. This task fits nicely into Problem (4) by setting the cost function to the sum of squared errors between the sought ultrametric and the edge weights of the given graph, namely

[TABLE]

Although the exact minimization of this cost function is a NP-hard problem [40], the proposed optimization framework allows us to compute an approximate solution. Figure 2 shows the ultrametric computed by Algorithm 1 with $J_{\rm closest}$ for an illustrative example of hierarchical clustering.

A common issue with the closest ultrametric is that small clusters might branch very high in the dendrogram. This is also true for average linkage and other agglomerative methods. Such kind of ultrametrics are undesirable, because they lead to partitions containing very small clusters at large scales, as clearly shown in Figures 2(b)-2(c). We now present two approaches to tackle this issue.

2.3 Cluster-size regularization

To fight against the presence of small clusters at large scales, we need to introduce a mechanism that pushes down the altitude of nodes where such incorrect merging occurs. This can be easily translated in our framework, as the altitude of a node corresponds to the ultrametric distance between its children. Specifically, we penalize the ultrametric distance proportionally to some non-negative coefficients that depend on the corresponding nodes in the dendrogram, yielding

[TABLE]

Here above, the $\gamma$ coefficients play an essential role: they must be small for the nodes that need to be pushed down, and large otherwise. We thus rank the nodes by the size of their smallest child, that is

[TABLE]

where $\rm Children_{u}(n)$ denotes the children of a node $n$ in the dendrogram $T_{u}$ associated to the ultrametric $u$ on $\mathcal{G}$. Figure 2(d) shows the ultrametric computed by Algorithm 1 with $J_{\rm closest}+J_{\rm size}$. The positive effect of this regularization can be appreciated by observing that small clusters are no longer branched very high in the dendrogram.

2.4 Triplet regularization

Triplet constraints [23, 27] provide an alternative way to penalize small clusters at large scales. Like in semi-supervised classification, we assume that the labels $\mathcal{L}_{v}$ of some data points $v\in V$ are known, and we build a set of triplets according to the classes they belong to:

[TABLE]

These triplets provide valuable information on how to build the ultrametric. Intuitively, we need a mechanism that reduces the ultrametric distance within the classes, while increasing the ultrametric distance between different classes. This can be readily expressed in our framework with a regularization acting on the altitude of nodes containing the triplets, leading to

[TABLE]

Here above, the constant $\alpha>0$ represents the minimum prescribed distance between different classes. Figure 2(e) shows the ultrametric computed by Algorithm 1 with $J_{\rm closest}+J_{\rm triplet}$.

2.5 Dasgupta cost function

Dasgupta cost function [9] has recently gained traction in the seek of a theoretically grounded framework for hierarchical clustering [16, 11, 41, 24, 27]. However its minimization is known to be a NP-hard problem [9]. The intuition behind this function is that large clusters should be associated to large dissimilarities. The idea it then to minimize, for each edge $e\in E$, the size of the dendrogram node associated to $e$ divided by the weight of $e$, yielding

[TABLE]

However, we cannot directly use (12) in our optimization framework, as the derivative of $|\textrm{lca}_{u}(x,y)|$ with respect to the underlying ultrametric $u$ is equal to 0 almost everywhere. To solve this issue, we propose a soft-cardinal measure of a node that is differentiable w.r.t. the associated ultrametric $u$. Let $n$ be a node of the dendrogram $T_{u}$, and let $\{x\}\subseteq n$ be a leaf of the sub-tree rooted in $n$. We observe that the cardinal of $n$ is equal to the number of vertices $y\in V$ such that the ultrametric distance $d_{u}(x,y)$ between $x$ and $y$ is strictly lower than the altitude of $n$, namely

[TABLE]

where $H$ is the Heaviside function. By replacing $H$ with a continuous approximation, such as a sigmoid function, we provide a soft-cardinal measure of a node of $T_{u}$ that is differentiable with respect to the ultrametric $u$. Figure 2(f) shows the ultrametric computed by Algorithm 1 with $J_{\rm Dasgupta}$.

Note that a differentiable cost function inspired by Dasgupta cost function was proposed in [20]. This function replaces the node size by a parametric probability measure which is optimized over a fixed tree. This is fundamentally different from our approach, where the proposed measure is a continuous relaxation of the node size, and it is directly optimized over the ultrametric distance.

3 Algorithms

In this section, we present a general approach to efficiently compute the various terms appearing in the cost functions introduced earlier. All the proposed algorithms rely on some properties of the single linkage (agglomerative) clustering, which is a dual representation of the subdominant ultrametric. We perform a detailed analysis of the algorithm used to compute the subdominant ultrametric. The other algorithms can be found in the supplemental material.

Single-linkage clustering can be computed similarly to a minimum spanning tree (m.s.t.) with Kruskal’s algorithm, by sequentially processing the edges of the graph in non decreasing order, and merging the clusters located at the extremities of the edge, when a m.s.t. edge is found. One consequence of this approach is that each node $n$ of the dendrogram representing the single linkage clustering is canonically associated to an edge of the m.s.t. (see Figure 3(a)), which is denoted by $\sigma(n)$.

In the following, we assume that we are working with a sparse graph $\mathcal{G}=(V,E)$, where $\mathcal{O}(|E|)=|V|$. The number of vertices (resp. edges) of $\mathcal{G}$ is denoted by $N$ (resp. $M$). For the ease of writing, we denote the edge weights as vectors of $\mathbb{R}^{M}$. The dendrogram corresponding to the single-linkage clustering of the graph $\mathcal{G}$ with edge weights $\tilde{w}\in\mathcal{W}$ is denoted by $\textrm{SL}(\tilde{w})$.

3.1 Subdominant ultrametric

To obtain an efficient and automatically differentiable algorithm for computing the subdominant ultrametric, we observe that the min-max distance between any two vertices $x,y\in V$ is given by the weight of the pass edge between $x$ and $y$. This is the edge holding the maximal value of the min-max path from $x$ to $y$, and an arbitrary choice is made if several pass edges exist. Moreover, the pass edge between $x$ and $y$ corresponds to the l.c.a. of $x$ and $y$ in the single linkage clustering of $(\mathcal{G},\tilde{w})$ (see Figure 3(a)). Equation (5) can be thus rewritten as

[TABLE]

The single-linkage clustering can be computed in time $\mathcal{O}(N\log N)$ with a variant of Kruskal’s minimum spanning tree algorithm [4, 42]. Then, a fast algorithm allows us to compute the l.c.a. of two nodes in constant time $\mathcal{O}(1)$, thanks to a linear time $\mathcal{O}(N)$ preprocessing of the tree [43]. The subdominant ultrametric can thus be computed in time $\mathcal{O}(N\log N)$ with Algorithm 2. Moreover, the dendrogram associated to the ultrametric returned by Algorithm 2 is the tree computed on line 2.

Note that Algorithm 2 can be interpreted as a special max pooling applied to the input tensor $w$, and can be thus automatically differentiated. Indeed, a sub-gradient of the min-max operator $\Phi$ at a given edge $e_{xy}$ is equal to 1 on the pass edge between $x$ and $y$ and 0 elsewhere. Then, the Jacobian of the min-max operator $\Phi$ can be written as the matrix composed of the indicator column vectors giving the position of the pass edge associated to the extremities of each edge in $E$:

[TABLE]

where $\Phi^{*}_{\mathcal{G}}(\tilde{w}_{i})$ denotes the index of the pass edge between the two extremities of the $i$-th edge, and $\mathbbm{1}_{j}$ is the column vector of $\mathbb{R}^{M}$ equals to 1 in position $j$, and 0 elsewhere (see Figure 3(b)).

3.2 Regularization terms

The cluster-size regularization defined in (8) can be implemented through the same strategy used in Algorithm 2, based on the single-linkage clustering and the fast l.c.a. algorithm, leading to a time complexity in $\mathcal{O}(N\log N)$. See supplemental material.

Furthermore, thanks to Equation (14), the triplet regularization defined in (11) can be written as

[TABLE]

This can be implemented with a time complexity in $\mathcal{O}(|\mathcal{T}|+N\log N)$. See supplemental material.

3.3 Dasgupta cost function

The soft cardinal of a node $n$ of the tree $T_{u}$ as defined in (13) raises two issues: the arbitrary choice of the reference vertex $x$, and the quadratic time complexity $\Theta(N^{2})$ of a naive implementation. One way to get rid of the arbitrary choice of $x$ is to use the two extremities of the edge $\sigma(n)$ canonically associated to $n$. To efficiently compute (13), we can notice that, if $c_{1}$ and $c_{2}$ are the children of $n$, then the pass edge between any element $x$ of $c_{1}$ and any element $y$ of $c_{2}$ is equal to edge $\sigma(n)$ associated to $n$. This allows us to define $\operatorname{card}(n)$ as the relaxation of $|n|$ in (13), by replacing the Heaviside function $H$ with the sigmoid function $\ell$, leading to

[TABLE]

where $\mathcal{A}(x)$ is the set of ancestors of $x$, and $c_{\hat{x}}(y)$ is the child of $y$ that does not contain $x$. The time complexity to evaluate (3.3) is dominated by the sum over the ancestors of $n$ which, in the worst case, is in the order of $\mathcal{O}(N)$, leading also to worst case time complexity of $\mathcal{O}(N^{2})$. In practice, dendrograms are usually well balanced, and thus the number of ancestors of a node is in the order of $\mathcal{O}(\log N)$, yielding an empirical complexity in $\mathcal{O}(N\log N)$.

4 Experiments

In the following, we evaluate the proposed framework on two different setups. The first one aims to assess our continuous relaxation of the closest ultrametric problem with respect to the (almost) exact solution provided by linear programming on planar graphs [14]. The second one aims to compare the performance of our cost functions to the classical hierarchical and semi-supervised clustering methods. The implementation of our algorithms, based on Higra [44] and PyTorch [45] libraries, is available at https://github.com/PerretB/ultrametric-fitting.

Framework validation

As Problem (6) is non-convex, there is no guarantee that the gradient descent method will find the global optimum. To assess the performance of the proposed framework, we use the algorithm proposed in [14], denoted by CUCP (Closest Ultrametric via Cutting Plane), as a baseline for the closest ultrametric problem defined in (7). Indeed, CUCP can provides an (almost) exact solution to the closest ultrametric problem for planar graphs based on a reformulation as a set of correlation clustering/multi-cuts problems with additional hierarchical constraints. However, CUCP requires to define a priori the set of levels which will compose the final hierarchy.

We generated a set of superpixels adjacency graphs of increasing scale from a high-resolution image (see Figure 4). The weight of the edge linking two superpixels is defined as the mean gradient value, obtained with [46], along the frontier between the two superpixels. The results presented in Figure 5 shows that the proposed approach is able to provide solutions close to the optimal ones (Figure 5(a)) using only a fraction of the time needed by the combinatorial algorithm (Figure 5(b)), and without any assumption on the input graph. The complete experimental setup is in the supplemental material.

The computation time of some combinations of cost terms are presented in Figure 5(c). Note that, Algorithm 1 usually achieves convergence in about one hundred iterations (see supplemental material). Closest and Closest+Size can handle graphs with millions of edges. Dasgupta relaxation is computationally more demanding, which decreases the limit to a few hundred thousands of edges.

Hierarchical clustering

We evaluate the proposed optimization framework on five datasets downloaded from the LIBSVM webpage,222https://www.csie.ntu.edu.tw/$\sim$cjlin/libsvmtools/datasets/ whose size ranges from $270$ to $1500$ samples. For each dataset, we build a 5-nearest-neighbor graph, to which we add the edges of a minimum spanning tree to ensure the connectivity. Then, we perform hierarchical clustering on this graph, and we threshold the resulting ultrametric at the prescribed number of clusters. We divide our analysis in two sets of comparisons: hierarchical clustering (unsupervised), and semi-supervised clustering. To be consistent among the two types of comparisons, we use the classification accuracy as a measure of performance.

Figure 6(a) compares the performance of three hierarchical clustering methods. The baseline is "Ward" agglomerative method, applied to the pairwise distance matrix of each dataset. Average linkage and closest ultrametric are not reported, as their performance is consistently worst. The "Dasgupta" method refers to Algorithm 1 with $J_{\rm Dasgupta}+\lambda J_{\rm size}$ and $\lambda=1$. The "Closest+Size" method refers to Algorithm 1 with the cost function $J_{\rm closest}+\lambda J_{\rm size}$ and $\lambda=10$. In both cases, the regularization is only applied to the top-10 dendogram nodes (see supplemental material). The results show that the proposed approach is competitive with Ward method (one of the best agglomerative heuristics). On the datasets Digit1 and Heart, "Dasgupta" performs slightly worse than "Closest+Size": this is partly due to the fact that our relaxation of the Dasgupta cost function is sensible to data scaling.

Figure 6(b) compares the performance of two semi-supervised clustering methods, and an additional unsupervised method. The first baseline is "Spectral" clustering applied to the Gaussian kernel matrix of each dataset. The second baseline is "SVM" classifier trained on the fraction of labeled samples, and tested on the remaining unlabeled samples. Between $10\%$ and $40\%$ of training samples are drawn from each dataset using a 10-fold scheme, and the cross-validated performance is reported in terms of mean and standard deviation. The "Closest+Triplet" method refers to Algorithm 1 with $J_{\rm closest}+\lambda J_{\rm triplet}$, $\lambda=1$ and $\alpha=10$. The results show that the triplet regularization performs comparably to semi-supervised SVM, which in turn performs better than spectral clustering.

5 Conclusion

We have presented a general optimization framework for fitting ultrametrics to sparse edge-weighted graphs in the context of hierarchical clustering. We have demonstrated that our framework can accommodate various cost functions, thanks to efficient algorithms that we have carefully designed with automatic differentiation in mind. Experiments carried on simulated and real data allowed us to show that the proposed approach provides good approximate solutions to well-studied problems.

The theoretical analysis of our optimization framework is beyond the scope of this paper. Nonetheless, we believe that statistical physics modelling [47] may be a promising direction for future work, based on the observation that ultrametricity is a physical property of spin-glasses [48, 49, 50, 51, 52, 53]. Other possible extensions include the end-to-end learning of neural networks for hierarchical clustering, possibly in the context of image segmentation.

6 Acknowledgment

This work was partly supported by the INS2I JCJC project under grant 2019OSCI. We are deeply grateful to Julian Yarkony and Charless Fowlkes for sharing their code, and to Fred Hamprecht for many insightful discussions.

Appendix A Supplemental material

A.1 Algorithms

This section describes in detail the algorithms proposed to compute the cost terms associated to cluster-size regularization, triplet regularization, and Dapgusta cost function.

Cluster-size regularization

This regularization penalizes small clusters at large scales (see Figure 7(a)), and the associated cost is computed by Algorithm 3. It proceeds by first computing the size of the smallest child of each node of the tree. The size of each node can be trivially computed recursively from the leaves to the root in linear time. Then, it identifies the pass node associated to each edge thanks to he fast l.c.a. algorithm, and it deduces the individual cost for this edge. Note that we assume the weights $\gamma_{u}$ are constants, even though they depend on the variable $u$ being optimized. This allows us to simplify the gradient evaluation. Moreover, the algorithm presents an additional hyper-parameter $k$ for applying the regularization only to the top-$k$ dendrogram nodes. In the algorithm we denote by $\textrm{rank}(n)$, the rank of a node $n$ according to the ordering given by their altitudes (from highest to lowest values). The root of the tree is thus ranked 1, the second highest node is ranked 2 and so on. Note that in practice, the single linkage algorithm (as any agglomerative clustering method) naturally orders nodes according to this ordering and no extra computation is required. We set $k=10$ in all our numerical experiments.

Triplet regularization

This regularization for semi-supervised learning enforces triplet constraints (see Figure 7(b)), and the associated cost is computed by Algorithm 4, which is very similar to the subdominant ultrametric algorithm. For every triplet $({\rm ref},{\rm pos},{\rm neg})$, it searches for $n_{1}$ and $n_{2}$, the smallest clusters containing the pairs $({\rm ref},{\rm pos})$ and $({\rm ref},{\rm neg})$, with the fast l.c.a. algorithm. It then introduces a penalization if the altitude of $n_{1}$ (i.e., the distance between ${\rm ref}$ and ${\rm pos}$) is not small enough compared to the altitude of $n_{2}$ (i.e., the distance between ${\rm ref}$ and ${\rm neg}$).

Dasgupta cost function

The difficulty in implementing the proposed relaxation of Dasgupta cost term lies in the soft-cardinal function defined in (3.3). The main function described in Algorithm 5 is similar to previously presented algorithm. The soft-cardinal function is computed by algorithm 6. As with Algorithm 3, the size of the nodes of the tree can be computed recursively from leaves to root in linear time. Note that, on line 9, the child of $y$ that does not contain $x$ can easily be determined by remembering the previous node of the "for each" loop: it is the sibling of the latter.

A.2 Framework validation

All the tests in the comparison with the CUCP algorithm (Section 4) were conducted on a computer with an Intel I7 4 cores processor and 16 GB of memory. For the optimization, we use the AMSGrad variation [54] of the ADAM method with step-size $0.01$. Our implementation of the proposed algorithms are all single threaded.

For each of the 52 test instances, the values $J_{\rm closest}$ obtained at each iteration of Algorithm 1 were normalized between 0 (lowest achieved cost for this instance) and 1 (highest cost reached for this instance). Figure 8 shows the mean-normalized convergence curve (with its standard deviation). We can see, that the convergence rate appears to be very good and smooth in practice. The convergence is usually reached within a bit more than a hundred iterations.

A.3 Illustrative examples

Figure 9 shows more illustrative examples of hierarchical clustering. For each dataset, we build a 5-nearest-neighbor graph, to which we add the edges of a minimum spanning tree to ensure the connectivity. Then, we perform hierarchical clustering on this graph, and we threshold the resulting ultrametric at the prescribed number of clusters. The column "Closest" is the solution to Algorithm 1 with the cost function $J_{\rm closest}$. The column "Closest+Size" is the solution to Algorithm 1 with the cost function $J_{\rm closest}+\lambda J_{\rm size}$ and $\lambda=0.1$, where the regularization is only applied to the top-10 dendogram nodes. The column "Closest+Triplet" is the solution to Algorithm 1 with the cost function $J_{\rm closest}+\lambda J_{\rm triplet}$, $\lambda=1$, and $\alpha=3$. The column "Dasgupta" is the solution to Algorithm 1 with the cost function $J_{\rm Dasgupta}$. For the optimization, we use the AMSGrad variation [54] of the ADAM method with step-size $0.1$.

A.4 Average linkage approximates the closest ultrametric problem

To help understand how the closest ultrametric problem relates to average linkage, note that an ultrametric takes a finite set of non-negative values $\{r_{1},\dots,r_{K}\}$ with $K<|V|$. Hence, it can be represented as

[TABLE]

where $\mathcal{L}=(\mathcal{L}_{1},\dots,\mathcal{L}_{K})$ are functions from $E$ to $\{0,1\}$ defining a hierarchical partition. In this setting, the cost function $J_{\rm closest}$ boils down to

[TABLE]

and, for a fixed hierarchical clustering $\bar{\mathcal{L}}$, the optimal altitudes are given by

[TABLE]

This is exactly the criterion used by average linkage to build a hierarchical clustering. We can thus argue that the latter provides an approximate solution to the closest ultrametric problem. As a matter of fact, average linkage and Algorithm 1 with $J_{\rm closest}$ produce structurally similar ultrametrics, as shown in Figures 2 and 9 for illustrative examples of hierarchical clustering.