Fitting Distances by Tree Metrics Minimizing the Total Error within a Constant Factor

Typically our measured distances do not have an exact fit with any tree metric. We then have the following generic optimization problem for any \(L_p\)-norm:

Problem: \(L_p\)-fitting tree (ultra) metrics.

Input: A set S with a distance function \(\mathcal {D}:{S\choose 2}\rightarrow \mathbb {R}_{\gt 0}\).⁶

Desired Output: A tree metric (or ultrametric) T that spans S and fits \(\mathcal {D}\) in the sense of minimizing the \(L_p\)-normCavalli-Sforza and Edwards [11] introduced this numerical taxonomy problem for both tree and ultrametrics in the \(L_2\)-norm in 1967. Farris suggested using \(L_1\)-norm in 1972 [35, p. 662].

In this article, we focus on the \(L_1\)-norm, that is, the total sum of errors. The problem is APX-hard for both tree metrics and ultrametrics (see Section 9 and Reference [3]), so a constant approximation factor is the best we can hope for in polynomial time. The best previous approximation factor for both tree metrics and ultrametrics was \(O((\log n)(\log \log n))\) by Ailon and Charikar [3].

In this article, we present a deterministic polynomial time constant factor approximation both for tree metrics and for ultrametrics, that is, in both cases, we can find a tree T minimizing the \(L_1\)-norm within a constant factor of the best possible.

Thus, we will prove the following theorem.

Since Cavalli-Sforza and Edwards introduced the tree fitting problem, the problem has collected an extensive literature. In 1977 [54], it was shown that if there is a tree metric coinciding exactly with \(\mathcal {D}\), then it is unique and it can be found in time linear in the input size, i.e., \(O(|S|^2)\) time. The same then also holds trivially for ultrametrics. Unfortunately there is typically no tree metric coinciding exactly with \(\mathcal {D}\), and in 1987 [28] it was shown that for \(L_1\) and \(L_2\) the numerical taxonomy problem is NP-hard in the tree metric and the ultrametric cases. The problems are in fact APX-hard (see Section 9), which rules out the possibility of a polynomial-time approximation scheme. Thus, a constant factor, like ours for \(L_1\), is the best one can hope for from a complexity perspective for these problems.

For the \(L_\infty\) numerical taxonomy problem, there was much more progress. In 1993 [34], it was shown that for the ultrametric case an optimal solution can be found in \(\Theta (|S|^2)\) time. More recently, it was shown that when the points are embedded into \(\mathbb {R}^d\) and the distances are given by the pairwise Euclidean distances, the problem can be approximated in subquadratic time [22, 25]. For the general trees case (still in the \(L_\infty\)-norm objective), Reference [2] gave an \(O(|S|^2)\) algorithm that produces a constant factor approximation and proved that the problem is APX-hard (unlike the ultrametric case).

The technical result from Reference [2] was a general reduction from general tree metrics to ultrametrics. It modifies the input distance matrix and asks for fitting this new input with an ultrametric that can later be converted to a tree metric for the original distance matrix. The result states that for any p, if we can minimize the restricted ultrametric \(L_p\) error within a factor \(\alpha\) in polynomial-time, then there is a polynomial-time algorithm that minimizes the tree metric \(L_p\) error within a factor \(3\alpha\). The reduction from Reference [2] imposes a certain restriction on the ultrametric, but the restriction is not problematic, and in Section 8, we will even argue that the restriction can be completely eliminated with a slightly modified reduction. With n species, the reduction from tree metrics to ultrametrics can be performed in time \(O(n^2)\). Applying this to the optimal ultrametric algorithm from Reference [34] for the \(L_\infty\)-norm objective yielded a factor 3 for general metrics for the \(L_\infty\)-norm objective. The generic impact is that for any \(L_p\), later algorithms only had to focus on the ultrametric case to immediately get results for the often more important tree metrics case, up to losing a factor 3 in the approximation guarantee. Indeed, the technical result of this article is a constant factor approximation for ultrametric. Thus, letting TreeMetric and UltraMetric respectively denote the approximation factors of the tree metric and ultrametric problems, we have

For \(L_p\)-norms with constant p, the developments have been much slower. Ma et al. [44] considered the problem of finding the best \(L_p\) fit by an ultrametric where distances in the ultrametric are no smaller than the input distances. For this problem, they obtained an \(O(n^{1/p})\) approximation.

Later, Dhamdhere [29] considered the problem of finding a line metric to minimize additive distortion from the given data (measured by the \(L_1\)-norm) and obtained an \(O(\log n)\) approximation. In fact, his motivation for considering this problem was to develop techniques that might be useful for finding the closest tree metric with distance measured by the \(L_1\)-norm. Harb, Kannan, and McGregor [40] developed a factor \(O(\min \lbrace n,k\log n\rbrace ^{1/p})\) approximation for the closest ultrametric under the \(L_p\)-norm where k is the number of distinct distances in the input.

The best bounds known for the ultrametric variant of the problem are due to Ailon and Charikar [3]. They first focus on ultrametrics in \(L_1\) and show that if the distance matrix has only k distinct distances, then it is possible to approximate the \(L_1\) error within a factor \(k+2\). Next they obtain an LP-based \(O((\log n)(\log \log n))\) approximation for arbitrary distances matrices. Finally, they sketch how it can be generalized to an \(O(((\log n)(\log \log n))^{1/p})\) approximation of the \(L_p\) error for any p. Using the reduction from Reference [2], they also get an \(O(((\log n)(\log \log n))^{1/p})\) approximation for tree metrics under the \(L_p\)-norm objective. The \(O(((\log n)(\log \log n))^{1/p})\) approximation comes from an \(O((\log n)(\log \log n))\) approximation of the pth moment of the following function:Technically, Ailon and Charikar [3] present a simple LP relaxation for \(L_1\) ultrametric fitting—an LP that will also be used in our article. They get their \(O((\log n)(\log \log n))\) approximation using an LP rounding akin to the classic \(O(\log n)\) rounding of Leighton and Rao for multicut [43]. The challenge is to generalize the approach to deal with the hierarchical issues associated with ultrametric and show that this can be done paying only an extra factor \(O(\log \log n)\) in the approximation factor. Next they show that their LP formulation and rounding is general enough to handle different variants, including other \(L_p\)-norms as mentioned above, but also they can handle the weighted case, where for each pair of species \(i,j\), the error contribution to the overall error is multiplied by a value \(w_{ij}\). However, this weighted problem captures the multicut problem (and the weighted minimization correlation clustering problem) [3]. Since the multicut cannot be approximated within a constant factor assuming the unique games conjecture [19] and the best-known approximation bound remains \(O(\log n)\), it is beyond reach of current techniques to do much better in these more general settings.

Ailon and Charikar [3] conclude that “determining whether an \(O(1)\) approximation can be obtained is a fascinating question. The LP formulation used in our [their] work could eventually lead to such a result.” For their main LP formulation for the (unweighted) \(L_1\) ultrametric fitting, the integrality gap was only known to be somewhere between 2 and \(O((\log n) (\log \log n))\). To break the \(\log n\)-barrier, we must come up with a radically different way of rounding this LP and free ourselves from the multicut-inspired approach.

For \(L_1\) ultrametric fitting, we give the first constant factor approximation, and we show this can be obtained by rounding the LP proposed by Ailon and Charikar, thus demonstrating a constant integrality gap for their LP. Our solution breaks the \(\log n\) barrier using the special combinatorial structure of the \(L_1\) problem.

Stepping a bit back, having different solutions for different norms should not come as a surprise. As an analogue, take the generic problem of placing k facilities in such a way that each of n cities is close to the nearest facility. Minimizing the vector of distances in the \(L_1\)-norm is called the k-median problem. In the \(L_2\)-norm, it is called the k-means problem and in the \(L_\infty\)-norm the k-center problem. Indeed, while the complexity of the k-center problem has been settled in the mid-1980s thanks to Gonzalez’s algorithm [38], it remained a major open problem for the next 15 years to obtain constant factor approximation for the k-median and the k-means problems. Similarly, our understanding of the k-means problem (\(L_2\)-objective) remains poorer than our understanding of the k-median problem, and the problem is in fact provably harder (no better than \(1+8/e\)-approximation algorithm [39], while k-median can be approximated within a factor 2.675 [9]).

For our tree fitting problem, the \(L_\infty\)-norm has been understood since the 1990s, and our result shows that the \(L_1\)-norm admits a constant factor approximation algorithm. The current status of affairs for tree and ultrametrics is summarized in Table 1. The status for \(L_p\) tree fitting is that we have a good constant factor approximation if we want to minimize the total error \(L_1\) or the maximal error \(L_\infty\). For all other \(L_p\)-norms, we only have the much weaker but general \(O(((\log n)(\log \log n))^{1/p})\) approximation from Reference [3]. In particular, we do not know if anything better is possible with \(L_2\). The difference is so big that even if we are in a situation where we would normally prefer an \(L_2\) approximation, our much better approximation guarantee with \(L_1\) might be preferable.

Computational Biology. Researchers have also studied reconstruction of phylogenies under stochastic models of evolution (see Farach and Kannan [33] or Mossel et al. [46] and the references therein; see also Henzinger et al. [41]).

Finally, related to the hierarchical correlation clustering problem that we introduce in this article is the hierarchical clustering problem introduced by Dasgupta [27] where the goal is, given a similarity matrix, to build a hierarchical clustering tree where the more similar two points are, the lower in the tree they are separated (formally, a pair \((u,v)\) induces a cost of similarity(\(u,v\)) times the size of the minimal subtree containing both u and v, the goal is to minimize the sum of the costs of the pairs). This has received a lot of attention in the past few years (References [5, 14, 15, 16, 18, 23, 24, 45, 47], see also References [1, 6, 13, 21]) but differs from our settings, since the resulting tree may not induce a metric space.

Metric Embeddings. There is a large body of work of metric embedding problems. For example, the metric violation distance problem asks to embed an arbitrary distance matrix into a metric space while minimizing the \(L_0\)-objective (i.e., minimizing the number of distances that are not preserved in the metric space). The problem is considerably harder and is only known to admit an \(O(\text{OPT}^{1/3})\)-approximation algorithm [31, 32, 36], while no better than a 2 hardness of approximation is known. More practical results on this line of work includes References [51] and [37]. Sidiropoulos et al. [48] also considered the problem of embedding into metric, ultrametric, and so on, while minimizing the total number of outlier points.

There is also a rich literature on metric embedding problems where the measure of interest is the multiplicative distortion, and the goal of the problem is to approximate the absolute distortion of the metric space (as opposed to approximating the optimal embedding of the metric space). Several such problems have been studied in the context of approximating metric spaces via tree metrics (e.g., References [8, 30]). The objective of these works is very different, since they are focused on the absolute expected multiplicative distortion over all input metrics while we aim at approximating the optimal expected additive distortion for each individual input metric.

While the embedding techniques developed in References [8, 30] and related works have been very successful for designing approximation algorithms for various problems in a variety of contexts, they are not aimed at numerical taxonomy. Their goal is to do something for general metrics. However, for our tree-fitting problem, the idea is that the ground truth is a tree, e.g., a phylogeny, and that the distances measured, despite noise and imperfection of the model, are close to the metric of the true tree. To recover an approximation to the true tree, we therefore seek a tree that compares well against the best possible fit of a tree metric.

We will now discuss the main idea of our algorithm. Our solution will move through several combinatorial problems that code different aspects of the \(L_1\)-fitting of ultrametrics but that do not generalize nicely to other norms.

Our result follows from a sequence of constant-factor-approximation reductions between problems. To achieve our final goal, we introduce several new problems that have a key role in the sequence of reductions. Some of the reductions and approximation bounds have already been extensively studied (e.g., Correlation Clustering). A roadmap of this sequence of results is given in Figure 1. While some the reductions and techniques we use for our new constant factor approximation were known, we do not make use of anything published since the previous \(O((\log n)(\log \log n))\) factor approximation of Ailon and Charikar [4].

Our main technical contribution in this article is solving Hierarchical Correlation Clustering. Reductions from Ultrametric to Hierarchical Correlation Clustering, and from general Tree Metric to Ultrametric, are already known from References [3, 40] and [2], respectively. We discuss both reductions in Sections 7 and 8. This includes removing some restrictions in the reduction from Tree Metrics to Ultrametrics.

Focusing on Hierarchical Correlation Clustering, our input is the \(\ell\) weights \(\delta ^{(1)},\ldots ,\delta ^{(\ell)}\in \mathbb {R}_{\gt 0}\) and \(\ell\) edge sets \(E^{(1)},\ldots ,E^{(\ell)}~\subseteq ~{S\choose 2}\).

Step 1: Solve correlation clustering independently for each level. The first step in our solution is to solve the correlation clustering problem defined by \(E^{(t)}\) for each level \(t=1,\ldots ,\ell\) independently, thus obtaining an intermediate partitioning \(Q_t\). As we mentioned in Section 1.5.1, this can be done so that \(Q_t\) minimizes \(|E^{(t)}\triangle E(Q_t)|\) within a constant factor.

Step 2: Solve hierarchical cluster agreement. We now use the \(\ell\) weights \(\delta ^{(1)},\ldots ,\delta ^{(\ell)}\in \mathbb {R}_{\gt 0}\) and \(\ell\) partitions \(Q^{(1)},\ldots ,Q^{(\ell)}\) of S as input to the hierarchical cluster agreement problem, which we solve using an LP very similar to the one Ailon and Charikar [3] used to solve general hierarchical correlation clustering within a \(O((\log n)(\log \log n))\) factor. However, when applied to the special case of hierarchical cluster agreement, it allows an unusual kind of constant-factor LP rounding. We use the LP variables (which only indicate how much pairs of species should be together) to decide which sets from the input partitions are important to the hierarchy and which sets can be ignored. The important sets may not yet be hierarchical (or laminar), but we show that some small modifications suffice. This is done by a relatively simple combinatorial algorithm that modify the important sets bottom-up to generate the hierarchical output partitions \(P^{(1)},\ldots , P^{(\ell)}\). The result is a poly-time constant factor approximation for hierarchical cluster agreement, that is,The output partitions \(P^{(1)},\ldots , P^{(\ell)}\) are also returned as output to the original hierarchical correlation clustering problem.

We now provide a high level overview and the intuition behind the hierarchical cluster agreement algorithm. The algorithm can be broadly divided into two parts.

LP cleaning. We start by optimally solving the LP based on the weights \(\delta ^{(1)},\dots ,\delta ^{(\ell)}\) and partitions \(Q^{(1)},\dots ,Q^{(\ell)}\). For each level \(t\in \lbrace 1,\ldots , \ell \rbrace\), we naturally think of the relevant LP variables as distances and call them LP distances. That is because a small value means that the LP wants the corresponding species to be in the same part of the output partition at level t and vice versa, while the LP constraints also enforce the triangle inequality. The weights \(\delta ^{(1)},\dots ,\delta ^{(\ell)}\) impact the optimal LP variables but will otherwise not be used in the rest of the algorithm.

Using the LP distances, we clean each set in every partition \(Q^{(t)}\) independently. The objective of this step (LP-Cleaning - Algorithm 2) is to keep only the sets where most species are very close to each other and far away from the species not in the set. All other sets are disregarded. We process the sets of each level independently. The algorithm may only slightly modify sets that are not disregarded when processing their level. The property that we can clean each set independently to decide whether it is important or not, without looking at any other sets makes this part of our algorithm quite simple.

Omitting exact thresholds for simplicity, the algorithm works as follows. We process each set \(C_I\in Q^{(t)}\) by keeping only those species that are at very small LP distance from at least half of the other species in \(C_I\) and at large LP distance to almost all the species outside \(C_I\). Let us note that by triangle inequality and the pigeonhole principle, all species left in a set are at relatively small distance from each other. After this cleaning process, we only keep a set if at least \(90\%\) of its species are still intact, and we completely disregard it otherwise. The LP cleaning algorithm outputs the sequence \(L^{(*)}=(L^{(1)},\dots ,L^{(\ell)})\), where \(L^{(t)}\) is the family of surviving cleaned sets from \(Q^{(t)}\).

Derive hierarchy. Taking \(L^{(*)}\) as input, in the next step the algorithm Derive-Hierarchy (Algorithm 3) computes a hierarchical partition \(P^{(*)}=(P^{(1)},\dots ,P^{(\ell)})\) of S. This algorithm works bottom-up while initializing an auxiliary bottom most level of the hierarchy with \(|S|\) sets where each set is a singleton and corresponds to a species of S. Then the algorithm performs \(\ell\) iterations where at the tth iteration it processes all the disjoint sets in \(L^{(t)}\) and computes partition \(P^{(t)}\) while ensuring that at the end of the iteration \(P^{(1)},\dots ,P^{(t)}\) are hierarchical. An interesting feature of our algorithm is that while creating \(P^{(t)}\) we do not make any changes to the lower partitions \(P^{(1)},\dots ,P^{(t-1)}\). Next, we discuss how to compute \(P^{(t)}\) given \(L^{(t)}\) and the lower partitions \(P^{(1)},\dots ,P^{(t-1)}\).

Consider a set \(C_{LP}\in L^{(t)}\). Now if for each lower level set \(C^{\prime }\) either \(C^{\prime }\cap C_{LP}=\emptyset\) or \(C^{\prime }\subseteq C_{LP}\), then introducing \(C_{LP}\) at level t does not violate the hierarchy property. Otherwise, let \(C^{\prime }\) be a lower level set such that \(C^{\prime }\cap C_{LP}\ne \emptyset\) and \(C^{\prime }\not\subseteq C_{LP}\). Note that we already mentioned, once created, \(C^{\prime }\) is never modified while processing upper level sets. Thus, to ensure the hierarchy condition, the algorithm can either extend \(C_{LP}\) so that it completely covers \(C^{\prime }\) or can discard the common part from \(C_{LP}\).

In the process of modifying \(C_{LP}\) (where we can add or discard some species from it), at any point we define the core of \(C_{LP}\) to be the part that comes from the set created initially. Now to resolve the conflict between \(C_{LP}\) and \(C^{\prime }\), we work as follows. If the core of \(C_{LP}\) intersects the core of \(C^{\prime }\), then we extend \(C_{LP}\) so that \(C^{\prime }\) becomes a subset of it. Omitting technical details, there are two main ideas here: First, we ensure that the number of species in \(C^{\prime }\) (respectively, \(C_{LP})\) that are not part of its core is negligible with respect to the size of \(C^{\prime }\) (respectively, \(C_{LP})\). Furthermore, since the cores of \(C_{LP}, C^{\prime }\) have at least one common species, using triangle inequality we can claim that any pair of species from the cores of \(C^{\prime }, C_{LP}\) also have small LP distance; therefore, nearly all pairs of species in \(C_{LP}, C^{\prime }\) have small LP distance, meaning that the extension of \(C_{LP}\) is desirable (i.e., its cost is within a constant factor from the LP cost while it ensures the hierarchy).

Here we want to emphasize the point that because of the LP-cleaning, we can ensure that for any lower level set \(C^{\prime }\) at level t there exists at most one set whose core has an intersection with the core of \(C^{\prime }\). We call this the hierarchy-friendly property of the LP cleaned sets. This property is crucial for consistency, as it ensures that at level t there cannot exist more than one sets that are allowed to contain \(C^{\prime }\) as a subset.

In the other case, where the cores of \(C_{LP}\) and \(C^{\prime }\) do not intersect, the algorithm removes \(C_{LP}\cap C^{\prime }\) from \(C_{LP}\). The analysis of this part is more technical but follows the same arguments, namely using the aforementioned properties of LP-cleaning along with triangle inequality.

After processing all the sets in \(L^{(t)}\), the algorithm naturally combines these processed sets with \(P^{(t-1)}\) to generate \(P^{(t)}\), thus ensuring that \(P^{(1)}, \ldots , P^{(t)}\) are hierarchical.

High-level analysis. We will prove that the partitions \(P^{(1)},\ldots , P^{(\ell)}\) solves the original hierarchical clustering problem within a constant factor.

Using triangle inequality, we are going to show that the switch in Step 1, from the input edge sets \(E^{(1)},\ldots ,E^{(\ell)}\) to the partitions \(Q^{(1)},\ldots ,Q^{(\ell)}\) costs us no more than the approximation factor of correlation clustering used to generate each partition. This then becomes a multiplicative factor on our approximation factor for hierarchical cluster agreement, more specifically,We will even show that we can work with the LP from Reference [3] for the original hierarchical correlation clustering problem, and get a solution demonstrating a constant factor integrality gap.

In Section 2, we present the LP formulation and related definitions for Hierarchical Correlation Clustering. In Section 3, we show how to reduce Hierarchical Correlation Clustering to Hierarchical Cluster Agreement. In Section 4, we present the algorithm for Hierarchical Cluster Agreement, and in Section 5 we analyze it. In Section 6, we show that the LP formulation for Hierarchical Correlation Clustering has constant integrality gap. In Section 7, we show how \(L_p\)-fitting ultrametrics reduces to Hierarchical Correlation Clustering. In Section 8, we discuss the reduction from \(L_p\)-fitting tree metrics to \(L_p\)-fitting ultrametrics. In Section 9, we prove APX-Hardness of \(L_1\)-fitting ultrametrics and \(L_1\)-fitting tree metrics. We conclude in Section 10.

The IP-formulation of Hierarchical Cluster Agreement is akin to the IP-formulation of Hierarchical Correlation Clustering. Namely, the constraints are exactly the same for both problems. The only difference is in the objective function where we replace the general edge-sets \(E^{(1)}, \ldots , E^{(\ell)}\) with the disjoint clique edges from \(\mathcal {E}(Q^{(1)}), \ldots , \mathcal {E}(Q^{(\ell)})\) and similarly for the LP-relaxation of Hierarchical Cluster Agreement. Here each component in \(Q^{(t)}\) is called a level-t input-cluster.

To simplify our discussion, we use \(x^{(*)}\) to denote a fractional solution to the LP-relaxation of Hierarchical Cluster Agreement, that is, a vector containing all \(x_{i,j}^{(t)}, \lbrace i,j\rbrace \in {S\choose 2}, t\in [\ell ]\). One can think of \(x^{(*)}\) as the optimal fractional solution, but in principle it can be any solution.

We use \(x^{(t)}\), for some particular \(t\in [\ell ]\), to denote the vector containing all \(x_{i,j}^{(t)}, \lbrace i,j\rbrace \in {S\choose 2}\).

As previously, we use the term LP distances to refer to the entries of \(x^{(*)}\) and notice that for any particular \(t\in [l]\) the LP distances even satisfy the triangle inequality by the LP constraints.

Given \(x^{(*)}\), we define \(B_{\lt r}^{(t)}(i)\) to be the ball of species with LP-distance less than r from i at level t. More formally, \(B_{\lt r}^{(t)}(i) = \lbrace j\in S \mid x_{i,j}^{(t)}\lt r\rbrace\). Similarly, for a subset \(S^{\prime }\) of S we define the ball \(B_{\lt r}^{(t)}(S^{\prime }) = \lbrace j\in S \mid \exists i\in S^{\prime } \text{ s.t. }x_{i,j}^{(t)}\lt r\rbrace\).

We also define the LP cost of species \(i,j\) at level t as

as well as the LP cost of species in a set \(S^{\prime }\subseteq S\) at level t as

and in case \(S^{\prime }\) only contains a single species i, we write \(cost_i^{(t)}\) instead of \(cost_{\lbrace i\rbrace }^{(t)}\).

Then the LP cost at level t is denoted as \(cost^{(t)} = cost_{S}^{(t)}\).

Finally, the LP cost is simply \(cost^{(*)} = \sum _{t=1}^{\ell } cost^{(t)}\).

The pseudocode for our main algorithm for Hierarchical Cluster Agreement is given in Algorithm 1.

Our LP relaxation has size polynomial in \(S, \ell\), and the two subroutines also run in polynomial time, as we show later. Therefore, the whole algorithm runs in polynomial time.

In Algorithm 2, we provide the pseudocode of the LP Cleaning step of our algorithm.

Intuitively, the aim of this algorithm is to clean the input sets so that (ideally) all species remaining in a set have small LP distances to each other and large LP distances to species not in the set.

Formally, Algorithm 2 takes a sequence \(Q^{(*)}=(Q^{(1)},\dots ,Q^{(\ell)})\) of partitions of S and a fractional solution \(x^{(*)}\) containing LP distances. It outputs a sequence of families of disjoint subsets of S, \(L^{(*)}=(L^{(1)},\dots ,L^{(\ell)})\). Here each component of \(L^{(t)}\) is called a level-t LP-cluster.

In the algorithm, for each input partition \(Q^{(t)}\) we process every level-t input-cluster \(C_I\in Q^{(t)}\) separately. For this, we remove all the species in \(C_I\) that do not have very small LP distance to at least half the species in \(C_I\) or that have small LP distance to many species not in \(C_I\). More formally, we remove all the species in \(C_I\) with LP distance less than 0.1 to at most half the species in \(C_I\) or with LP distance less than 0.6 to more than \(0.05|C_I|\) species not in \(C_I\).

After the cleaning step, we discard \(C_I\) if smaller than a \(9/10\) fraction of the species survive. Otherwise, we create an LP-cluster \(C_{LP}\) containing the species in \(C_I\) that survive. Next, we add the level-t LP-cluster \(C_{LP}\) to \(L^{(t)}\).

Of several properties that we prove concerning the output of the LP-Cleaning, we briefly mention the following one: The output sequence \(L^{(*)}\) is hierarchy-friendly in the sense that no two LP-clusters at the same level t can be intersected by the same LP-cluster at level \(t^{\prime }\lt t\). We formally prove this in Lemma 5.4.

The LP-Cleaning subroutine trivially runs in time polynomial in \(S,\ell\).

In this section, we introduce Derive-Hierarchy (Algorithm 3). It takes as input a hierarchy-friendly sequence \(L^{(*)}=(L^{(1)}, \ldots , L^{(\ell)})\) of families of disjoint subsets of S and outputs a sequence \({P^{(*)}=(P^{(1)}, \ldots , P^{(\ell)})}\) of hierarchical partitions of S. The execution of the algorithm can be seen, via a graphical example, in Figure 2.

The algorithm works bottom-up while performing \(\ell\) iterations for \(t=1,\dots ,\ell\). In the process, it incrementally builds a forest \(\mathcal {F}\). Throughout the algorithm, each non leaf node u in \(\mathcal {F}\) can be identified by an LP-cluster in \(L^{(*)}\). Moreover, for each node u the algorithm maintains two sets \(C(u) \subseteq S\), called the core-cluster of u, and \(C^{+}(u) \subseteq S\), called the extended-cluster of u.

The algorithm starts by initializing \(\mathcal {F}\) with \(|S|\) trees where each tree contains a single node \(u_i\) identified by a species \(i\in S\). Also, it initializes both sets \(C(u_i),C^{+}(u_i)\) with \(\lbrace i\rbrace\). Next, in iteration t the algorithm processes the LP-clusters in \(L^{(t)}\) and at the end of the iteration, the \(C^{+}()\) sets associated with the root nodes in \(\mathcal {F}\) define the partition \(P^{(t)}\). Precisely here, the \(C^{+}()\) set of a root node contains all the species descending from the respective root and will thus be the sets in the final partitions.

In the tth iteration, for each cluster \(C_{LP}\in L^{(t)}\) the algorithm adds a root node u to forest \(\mathcal {F}\) while initializing the set \(C(u)\) with \(C_{LP}\). Next for the root node u the algorithm decides on its children by processing the pre-existing roots in the following way. For consistency, first it detects all the pre-existing root nodes v such that \(C(u)\) does not intersect \(C(v)\). Then it removes from \(C(u)\) all the species that are descending from v, i.e., sets \(C(u)\leftarrow C(u)\setminus C^{+}(v)\). Last, it sets u as a parent of all other pre-existing root nodes v such that \(C(u)\) intersects \(C(v)\). Also, accordingly, it modifies the set of leaf nodes of the subtree rooted at u by setting \(C^+(u)\leftarrow C^ +(u)\cup C^+(v)\). Notice here that some of the root-nodes may correspond to sets from levels lower than t, in case no parent was assigned to them.

At the end of iteration t, the algorithm completes processing all the LP-clusters in \(L^{(1)},\dots ,L^{(t)}\) and constructs partitions \(P^{(1)},\dots ,P^{(t)}\). At the end of the \(\ell\) iterations it outputs the \(\ell\) partitions \(P^{(*)}=(P^{(1)},\dots ,P^{(\ell)})\).

The Derive-Hierarchy subroutine trivially runs in time polynomial in \(S,\ell\).

Notice that throughout the execution of the algorithm, \(\mathcal {F}\) is an incrementally updated graph (that is, no deletions occur). In fact, it is always a forest, as we start with \(|S|\) isolated nodes and only introduce new nodes as parents of roots of some of the existing trees. Moreover, this process implies that the subtree rooted at any specific node is never modified.

From now on, we use \(\mathcal {F}\) to refer to the final instance of the incrementally updated forest. We use \(\mathcal {F}(u)\) to refer to the state of this incrementally updated forest after introducing u; therefore, \(\mathcal {F}(u)\setminus \lbrace u\rbrace\) denotes the state of the forest exactly before introducing node u. We naturally identify the leaves of \(\mathcal {F}\) with the species of S.

For any node u in the forest \(\mathcal {F}\), the Derive-Hierarchy algorithm defines \(C(u)\), which we call the core-cluster of u, and \(C^+(u)\), which we call the extended-cluster of u. Furthermore, notice that each core-cluster \(C(u)\) is a subset of some LP-cluster \(C_{LP}\) (Lines 7–9 of Algorithm 3); we call this the LP-cluster of u and denote it by \(L(u)\). Moreover, each LP-cluster \(L(u)\) is a subset of an input-cluster \(C_I\) (Line 3 of Algorithm 2); we call this the input-cluster of u and denote it by \(I(u)\). These concepts are well defined for any new node u and never change throughout the algorithm. We remind the reader that LP-Cleaning discards some of the input-clusters, in the sense that they have no corresponding LP-cluster, and therefore they do not match \(I(u)\), for any node u.

Directly from the algorithm we get that

To help with our discussion, we also define the following variables related to the Derive-Hierarchy algorithm (Algorithm 3):

For a node \(u\in \mathcal {F}\), we define its level \(t(u)\) to be the value of iteration t in Algorithm 3 when internal node u was introduced, and 0 when u is a leaf node.

We start with some observations that are heavily used in proving structural results regarding the core and the extended-clusters. These are in turn used for lower-bounding the LP cost.

The most important reason we are using the LP-Cleaning subroutine is so that any two species belonging in the same LP-Cluster at level t have small LP-distance.

For the analysis, it is convenient that our relations involve the LP-clusters instead of the input-clusters. Therefore, we rephrase Line 4 of Algorithm 2 in terms of LP-clusters, effectively proving that few species outside of an LP-cluster \(L(u)\) have small LP-distances to \(L(u)\).

The following lemma is just a convenient application of the triangle inequality of our LP, that is heavily used in subsequent proofs. Informally, it states that, under certain mild conditions, the LP-distance is small not only if \(i,j\) belong in the same LP-cluster (or core-cluster), but even if they happen to be in different clusters that are both intersected by the same third cluster.

We are now ready to prove the hierarchy-friendly property of the output of LP-Cleaning, as we informally claimed when introducing the algorithm. We claim that two LP-clusters of the same level cannot be intersected by the same lower level LP-cluster.

We finally present a simple lower bound on the LP cost. Recall thatand

In this section, we present several structural results related to our algorithm.

We start with pointing out that our algorithm ends up with the same output, no matter the order in which we process LP-clusters of the same level. This is due to the input sequence \(L^{(*)}\) being hierarchy-friendly.

Next, we prove two claims that we have already mentioned informally while introducing Algorithm 3 (Derive-Hierarchy). First we claim that the incrementally built graph is always a forest and also for any node u in F, its extended-cluster contains exactly the species that are identified with leaves descending from u in \(\mathcal {F}\). We call these species the descending species of u. We note that the results of Section 5.2 do not require these properties.

This simple lemma alone is enough to prove the following corollaries:

We also need the following result.

Using the developed toolkit of structural results, we are ready to show that for any node \(u\in \mathcal {F}\), all three of the LP-cluster \(L(u)\), the core-cluster \(C(u)\), and the extended-cluster \(C^+(u)\) are similar; more than that, we show lower bounds of the LP cost related to \(\Delta ^-(u) = L(u)\setminus C(u)\) and \(\Delta ^+(u) = C^+(u) \setminus C(u)\).

In particular, we claim that the following inequality holds for every \(u\in \mathcal {F}\),

We prove this claim inductively, based on the order in which nodes are added in \(\mathcal {F}\). As a base case, we initially create a node \(u_i\) for each species \(i\in S\) with \(C(u_i)=C^+(u_i)=\lbrace i\rbrace\), meaning that \(\Delta ^+(u_i) = \emptyset\).

For any other node u, we argue about the size of its extended-cluster \(C^+(u)\) in relation with the core-clusters of its descendants, as suggested by Lemma 5.12. We now partition the descendants of u in three parts and argue about each one of them.

For a node u, we define its top-non-intersecting descendants J as the set of highest level descendant nodes in \(\mathcal {F}\) whose core-clusters do not intersect \(C(u)\) (the reader is encouraged to consult Figure 3 before proceeding). More formally, using \(v \prec _{\mathcal {F}} u\) to denote that v is a descendant of u in forest \(\mathcal {F}\), we have

Notice that, by definition, if two nodes \(v,w\) belong in u’s top-non-intersecting descendants, then none is an ancestor of the other. Therefore u’s top-non-intersecting descendants J naturally partitions the proper descendants of u in three parts: J itself, the set \(J^{+}\) of proper descendants of nodes in J, and R containing the rest of the proper descendants of u (i.e., the proper descendants of u that are not descendants of any node in J). We also define sets of species related to these sets,

The apparent asymmetry of not excluding \(C(u)\) from \(S_J\) follows from the definition of u’s top-non-intersecting descendants J; the core-clusters of nodes in J are disjoint from \(C(u)\), meaning that \(S_J\) would be the same even if we excluded species in \(C(u)\). Note that this is not the case for the core-clusters in \(J^{+}\), as proper descendants of nodes in J might still intersect \(C(u)\), as in Figure 3.

Notice that by Lemma 5.12 we have \(C^+(u) = C(u) \cup (S_J\cup S_{J^+} \cup S_R)\), and thus

If \(v\in J \cup J^+ \cup R\), and its core-cluster does not intersect the core-cluster of u, then by definition of J we have that v is in descendants (not necessarily proper) of J. Therefore, \(v\in J\cup J^+\), meaning that nodes in R have core-clusters that intersect \(C(u)\). Furthermore, by definition, each node v in u’s top-non-intersecting descendants has a parent whose core-cluster intersects \(C(u)\). Therefore, for any species \(i\in C(u)\) and any species \(j\in S_J \cup S_R\), by Lemma 5.3 we have that their LP-distance is small, that is,

Species in \(S_J\cup S_R\) are not in \(C(u)\), and so by Corollary 5.10 they are not in \(L(u)\), as they belong in \(\Delta ^+(u)\). Then Lemma 5.2 gives

We are left to argue about species in \(J^+\), that is, in core-clusters of the descendants of u’s top-non-intersecting descendants. By Lemma 5.12, these species all belong in the extended-clusters of u’s top-non-intersecting descendants, \(\bigcup _{v\in J} C^{+}(v) \supseteq S_J \bigcup S_{J^+}\). By Corollary 5.11, these extended-clusters are disjoint, and thus \(S_{J^+} \subseteq \bigcup _{v\in J} \Delta ^+(v)\). By the inductive hypothesis (9), we get

Therefore, by Inequality (13) we get that

By Equations (13) and (15), we bound the size of \(\Delta ^+(u)\),

We are only left with bounding \(|L(u)|\). For this, we prove that

which, combined with Equation \((16)\), proves our initial claim, as we have

Before proving Equation \((17),\) we make some definitions (see Figure 4). Roughly speaking, we want to identify an appropriate set K of nodes such that the union of their extended-clusters both contains \(\Delta ^-(u)\) and its cardinality is reasonably boundable. In fact, the nodes in K are descendants of roots of \(\mathcal {F}(u)\setminus \lbrace u\rbrace\) that satisfy the condition in Line 8 of Algorithm 3 (i.e., nodes v such that \(C(v) \cap C(u) = \emptyset\), and \(C^{+}(v) \cap L(u) \not= \emptyset\)).

We now give a formal constructive definition of the set K. Let M be the set containing all the non-descendants of u at level at most \(t(u)\) whose core-clusters intersect \(L(u)\). We define \(K^{\prime }\) to be the set of parents of the nodes in M. Finally, K is obtained from \(K^{\prime }\) by removing the nodes who have a proper ancestor in \(K^{\prime }\). Notice by Corollary 5.11 that their extended-clusters are disjoint. We also define sets of species associated with K as follows:

Note \(\Delta ^-(u)\subseteq S_{K^+}\). Next, we claim for each node \(v\in K\), \(C(v)\cap L(u)=\emptyset\). Now if we can prove this claim, then it implies \(S_K\cap L(u)=\emptyset\), and thus we can write \(\Delta ^-(u)\subseteq S_{K^+}\setminus S_K\). Next, we prove the claim. Notice that for any node \(v \in M\), v is not a descendant of u but \(C(v)\cap L(u)\ne \emptyset\); thus there always exists a node \(w\in K\) such that w is an ancestor of v and \(C(w)\cap L(u)=\emptyset\).

Now, for the sake of contradiction, assume there exists a node \(w\in K\) such that \(C(w)\cap L(u)\ne \emptyset\). But then \(w\in M\), and following the previous argument there exists a node \(w^{\prime }\in K\) such that \(w^{\prime }\) is an ancestor of w and \(C(w^{\prime })\cap L(u)=\emptyset\). This is a contradiction, as by construction both w and \(w^{\prime }\) cannot be present in K.

Furthermore, notice that no node \(w\in K\) is at level \(t(w)= t(u)\), as that would imply a child \(w^{\prime }\in M\) of w, but \(C(w^{\prime })\) intersects \(C(w)\) (and therefore \(L(w)\)) as \(w^{\prime }\) is a child of w, and \(C(w^{\prime })\) intersects \(L(u)\), since \(w\in M\). By Lemma 5.4, this is a contradiction.

We conclude that K contains nodes at level at most \(t(u)-1\), which allows us to apply the inductive hypothesis \(|\Delta ^+(w)| \le 0.3|C(w)|\) for nodes \(w\in K\). Thus, from \(\Delta ^-(u) \subseteq S_{K^+} \setminus S_{K}\), we get

Furthermore, all nodes in K have a child whose core-cluster intersects \(C(u)\), and so by Lemma 5.3 the LP-distance between a species \(i\in L(u)\) and a species \(j\in S_K\) is small, \(x_{i,j}^{t(u)} \lt 0.6\). By Lemma 5.2, we get \(|S_K| \lt \frac{1}{6}|L(u)|\), which gives us \(|\Delta ^-(u)| \le \frac{1}{6}\cdot 0.3|L(u)|\).

By the definition of \(\Delta ^-(u)=L(u)\setminus C(u)\), we get \(|C(u)| \ge (1-\frac{1}{6}\cdot 0.3)|L(u)|\), by whichwhich concludes the proof of claim (17) and, as previously argued, the proof of claim (9).

As a by-product of this analysis, we can also give some lower bounds on the LP cost.

In this section, we prove that Algorithm 1 is an \(O(1)\) approximation of the LP cost.

We first make some definitions. Let \(t\in [\ell ]\). An input-cluster \(C_I\in Q^{(t)}\) is strong if there exists a level-t node \(u\in \mathcal {F}\) such that \(I(u)=C_I\). Similarly, a part P of the output partition \(P^{(t)}\) is strong if there exists a level-t node \(u\in \mathcal {F}\) such that \(C^+(u)=P\). In both cases, we say that u is the corresponding node. We characterize an input-cluster as weak if it is not strong, and similarly a part of the output partition \(P^{(t)}\) as weak if it is not strong.

We start with upper bounding the cost of Algorithm 1. The upper bound is related to the input-clusters (distinguishing between strong and weak) and the parts of the output partitions (again distinguishing between strong and weak). Informally, for weak input-clusters and weak parts, the cost of our algorithm is proportional to the sum of squares of their size. For a strong input-cluster with corresponding node u, the cost of our algorithm is proportional to its size times the number of species of the input-cluster that did not end up in u’s core-cluster. For a strong part of the output partitions, the cost of our algorithm is proportional to its size times the number of its species that did not end up in u’s core-cluster.

For each term of Lemma 5.15, we show a matching lower bound for the LP cost. First, we give a lower bound of the LP cost related to the weak input-clusters.

Next, we give a lower bound of the LP cost related to the strong input-clusters.

The following lemma lower bounds the LP cost in relation to the strong parts of the output partition \(P^{(t)}\).

The following lemma lower bounds the LP cost in relation to the weak parts of the output partition \(P^{(t)}\).

We are now ready to prove the following lemma:

Lemma 5.20 directly provesas Algorithm 1 simply picks \(x^{(*)}\) to be an optimal fractional solution to the LP relaxation.

Combining this with Inequality C concludes the proof of Theorem 3.1.

The claim from Reference [2] is that approximating a certain restricted ultrametric within a factor \(\alpha\) can be used to approximate tree pseudometric within a factor \(3\alpha\). Here we completely lift these restrictions for \(L_1\) and show that they can be lifted for all \(L_p\) with \(p\in \lbrace 2,3,\ldots \rbrace\) with an extra factor of at most 2.

We will need the well-known characterization of ultrametrics discussed in the Introduction that U is an ultrametric iff it is a metric and \(\forall \lbrace i,j,k\rbrace \in {S \choose 3}: U(i,j) \le \max \lbrace U(i,k),U(k,j)\rbrace\).

For simplicity, we prove the theorem only for \(p\lt \infty\), as for \(L_{\infty }\) the 3 approximation [2] cannot be improved by our theorem.

In this section, we prove that to find a good tree metric, it suffices to find a good tree pseudometric. This is a minor detail that we add for completeness. Informally, the construction simply replaces 0 distances with some parameter \(\epsilon\) and accordingly adapts the whole metric. By making the parameter \(\epsilon\) very small, the cost is not significantly changed.

Technically, our main lemma is the following.

Therefore, for any \(p\ge 1\), we can approximate \(L_p\)-fitting tree metrics by using an approximation to \(L_p\)-fitting tree pseudometrics. The error is at most \((1+\alpha)\) times the approximation factor of the tree pseudometric, as any tree metric is also a tree pseudometric.

Setting \(\alpha =\frac{1}{|S|}\) and using the result from Reference [2], we conclude that

This concludes the proof of Theorem 1.1.

As a final note, in the case of \(L_0\) (that is, we count the number of disagreements between \(\mathcal {D}\) and \(T^{\prime }\)) one cannot hope for a similar result. To see this, let S be a set of species, and let \(c_1,c_2\in S\) be two special species. The distance between any pair of species is 2, except if the pair contains either \(c_1\) or \(c_2\), in which case the distance is 1. The optimal tree pseudometric simply sets the distance between \(c_1\) and \(c_2\) to 0 and preserves everything else (1 disagreement).

Any tree metric requires at least \(|S|-3\) disagreements: We say that a non-special species is good if it has tree-distance 1 to both \(c_1\) and \(c_2\) and bad otherwise. Bad species have distance different than 1 to at least one special species, while good species have distance less than 2 with each other; the disagreements minimize at \(|S|-3\), when there is either one or two good species.

The correlation clustering problem has been shown to be APX-Hard in Reference [17]. As noted in References [3, 40] correlation clustering is the same as the \(L_1\)-fitting ultrametrics in case both the input and the output are only allowed to have distances in \(\lbrace 1,2\rbrace\). We refer to this problem as \(L_1\)-fitting \(\lbrace 1,2\rbrace\)-ultrametrics. Therefore, the \(L_1\)-fitting \(\lbrace 1,2\rbrace\)-ultrametrics is also APX-Hard.

For completeness, we sketch this relation here. Let \(E\subseteq {S\choose 2}\) be an instance of correlation clustering; then \(\mathcal {D}(i,j)\) is an instance to \(L_1\)-fitting \(\lbrace 1,2\rbrace\)-ultrametrics, where \(\mathcal {D}(i,j)=1\) if \(\lbrace i,j\rbrace \in E\), and \(\mathcal {D}(i,j)=2\) otherwise. Similarly, given \(\mathcal {D}\) we can obtain E by setting \(\lbrace i,j\rbrace \in S\) iff \(\mathcal {D}(i,j)=1\). Given any solution to correlation clustering (permutation P of S), we get a solution T to \(L_1\)-fitting \(\lbrace 1,2\rbrace\)-ultrametrics with \(T(i,j)=1\) if \(i,j\) are in the same part of P and \(T(i,j)=2\) otherwise. As T is an ultrametric, we are guaranteed that \(T(i,j)\le \max \lbrace T(i,k),T(j,k)\rbrace\); therefore, if \(T(i,k)=T(j,k)=1\), then \(T(i,j)=1\) as only distances in \(\lbrace 1,2\rbrace\) are allowed. Thus distance-1 is a transitive relation and P can be obtained by the equivalence classes of species with distance 1 in T. The observation from Reference [3] is that \(|E \triangle \mathcal {E}(P)| = \Vert T-\mathcal {D}\Vert _1\), which follows by trivial calculations.

The bird’s eye view of our approach for showing APX-Hardness of \(L_1\)-fitting ultrametrics is the following. For the sake of contradiction, we assume that \(L_1\)-fitting ultrametrics is not APX-Hard. We then show how to solve the \(L_1\)-fitting \(\lbrace 1,2\rbrace\)-ultrametrics problem in polynomial time within any constant factor greater than 1, contradicting the fact that it is APX-Hard. The main idea is that we first solve the general \(L_1\)-fitting ultrametrics problem. Then we apply a sequence of local transformations that converts the general ultrametric to an ultrametric with distances in \(\lbrace 1,2\rbrace\) without increasing the error. To achieve this, we first eliminate distances smaller than 1, then eliminate distances larger than 2, and then eliminate distances in \((1,2)\).

We first prove the following result concerning the local transformations. We remind the reader that an ultrametric T is defined as a metric with the property that for \(i,j,k \in S\) we have \(T(i,j) \le \max \lbrace T(i,k),T(j,k)\rbrace\).

We note that the ideas in Lemma 1 of Reference [40] can be directly used to prove Lemma 9.1. We only include Lemma 9.1, because we need these ideas applied to any ultrametric, not just to optimal ones, as Lemma 1 of Reference [40] does.

In this section, we show that \(L_1\)-fitting tree metrics is APX-Hard. Our reduction is based on the techniques used in Reference [28] to prove NP-Hardness of the same problem. The bird’s eye view of our approach is that we solve \(L_1\)-fitting ultrametrics by solving \(L_1\)-fitting tree metrics on a modified instance. In this instance, we introduce new species having small distance to each other and large distance to the original species. Through a sequence of local transformations, we show that we can modify the tree describing the obtained tree metric so as to consist of a star connecting the new species and an ultrametric tree connecting the original species (the center of the star and the root of the ultrametric tree are connected by a large edge). This ultrametric would refute APX-Hardness of \(L_1\)-fitting ultrametrics in case \(L_1\)-fitting tree metrics was not APX-Hard.