Computing generalized cophenetic distances under all Lp norms: A near-linear time algorithmic framework

TL;DR

This paper introduces a fast algorithm to compute generalized cophenetic distances across all Lp norms, enabling efficient and versatile large-scale tree comparisons in biology.

Contribution

A near-linear time algorithmic framework for computing generalized cophenetic distances under all Lp norms is presented.

Findings

The new algorithm enables large-scale studies by overcoming the limitations of previous methods restricted to L1 and L2 norms.

The framework's scalability is demonstrated through practical performance evaluations.

The distributions of key components of generalized cophenetic distances are analyzed under various vector norms.

Abstract

The cophenetic distance is a well-established metric in biology used to compare pairs of trees represented in a vector format. This distance was introduced by Cardona and his co-authors, building on the foundational work of Sokal and Rohlf, which dates back over 60 years. It is widely recognized for its versatility since it can analyze trees with edge weights using various vector norms. However, when comparing large-scale trees, the quadratic runtime of the current best-known (i.e., naïve) algorithm for computing the cophenetic distance can become prohibitive. Recently, a new algorithmic framework with near-linear time complexity has been developed to calculate the distances of a generalized class of cophenetic distances, which are derived from the work of Sokal and Rohlf. This improvement not only allows the cophenetic distance to be utilized in large-scale studies but also enhances the versatility of these studies by incorporating generalized variants of the cophenetic distance. However, the framework is limited to applying only the L1 and L2 vector norms, which significantly restricts the versatility of generalized cophenetic distances in large-scale applications. To address this limitation, we present a near-linear time algorithmic framework for computing the generalized cophenetic distances across all Lp vector norms. In our scalability study, we showcase the practical performance of our unrestricted algorithmic framework. Furthermore, we investigate the applicability of the generalized cophenetic distances by analyzing the distributions of key components of these distances under various vector norms.

Biological research often relies on large-scale comparisons of evolutionary trees to extract valuable insights across different subfields of biology. To effectively compare trees on a large scale, sensitive metrics are needed to assess differences in topology and branch lengths alongside efficient computational algorithms. In this study, we focus on the classic cophenetic distance, a metric that compares pairs of trees represented in vector format. The cophenetic distance can analyze trees with edge weights using various vector norms, making it highly versatile. Recently, a fast algorithm was developed to calculate generalized cophenetic distances, including the cophenetic distance itself, which has enabled large-scale studies and expanded the types of cophenetic distances applicable for tree pair comparisons. However, this algorithm is limited to L1 and L2 norms, which greatly reduces the versatility of the cophenetic distance. To overcome this limitation, we introduce a fast algorithm that computes the generalized cophenetic distance for all vector norms, making it suitable for large-scale studies. We also conduct a scalability study to demonstrate the effectiveness of our algorithm in practice, and we analyze the distributions of key representatives of cophenetic distances across various vector norms.

Figures

Cophenetic distances under various Lp norms and contribution functions.The vectors have values computed for the least common ancestor (lca, see the bottom row) of the corresponding taxon pairs. For example, lcaT(a,b) is v, and the second position in cophenetic vectors for T reflects the contribution of v, that is, the depth of v is 1.

Scalability analysis comparing the near-linear algorithm and the naïve quadratic-time algorithm.The diagrams show the average runtime with standard deviation bands over 1000 runs for each tree size (n=25,50,…,600) on pairs of random trees, evaluated for Lp-norms (p=1,2,5,10,20,50,100). Crossing points indicate the tree sizes where the near-linear algorithm outperforms the naïve one.

Scalability analysis comparing the near-linear algorithm and the naïve quadratic-time algorithm for small norms.The diagrams show the average runtime with standard deviation bands over 1,000 runs for each tree size (n=200,400,…,2600) on pairs of random trees, evaluated for Lp-norms (p=1,2,7,8).

Scalability analysis comparing the near-linear algorithm (outlined bars) and the naïve quadratic-time algorithm on large trees.The diagrams show the average runtime with standard deviation error over 100 runs for each tree size (n=1000,2000,…,6000) on pairs of random trees, evaluated for Lp-norms (p=1,2,5,10,20,50,100).

Sampled distribution of three cophenetic distances under the Lp norms for p = 1 to 10 (rows).Pairs of trees of size 100 were sampled according to the uniform model (left columns) and the Yule model (right columns). The frequencies are grouped into 200 bins. Note that the y-axis scale and range are the same across all diagrams. To enhance visibility, low frequencies (e.g., at low distances close to 0) are omitted and the width of the diagrams is appropriately adjusted.

(Cont. from Fig 5). Sampled distribution of three cophenetic distances is shown for Lp norms with p = 20, 50, 100, and ∞ (rows). Similar to Fig 5, the frequencies for finite p are grouped into 200 bins. Since all three contribution functions return integers, the L∞ norms are also integers. As a result, the frequency values for L∞ are presented without binning.

Mean and standard deviation from the sampled distributions for the uniform and Yule models for Lp norms with p>1.

Equations

Full text

Introduction

The cophenetic distance, introduced by Cardona et al. [1], originates from the pioneering work of Sokal and Rohlf over 60 years ago [2] and has gained substantial recognition in biology due to its reputation for reliability in analysis. This distance is a vector-based metric that is more versatile than many other commonly used tree metrics [3–5], as it can be applied to tree pairs with edge weights and analyzed using various vector norms. This capability allows for a more in-depth analysis, particularly when comparing similar tree topologies. Consequently, the cophenetic distance has broad applicability across various fields, such as phylogenetics [6–8], genomics [9], ecology [10–12], epidemiology [13,14], and conservation biology [15].

The cophenetic distance between a pair of rooted binary trees is defined based on their representation as cophenetic vectors. A cophenetic vector for a rooted tree assigns a value to each pair of taxa within the tree, representing the depth of their least common ancestor in that particular tree. Cophenetic vectors encode their corresponding trees equivalently [1], allowing the distance between a pair of trees to be measured as a Lp norm of the vector difference of the corresponding trees. Vector norms are beneficial not only for biological data analysis [16] but also applicable in various fields, such as statistical analyses [17] and machine learning [18]. An example of various vector norms is depicted in Fig 1.

Cophenetic distances under various Lp norms and contribution functions.The vectors have values computed for the least common ancestor (lca, see the bottom row) of the corresponding taxon pairs. For example, lcaT(a,b) is v, and the second position in cophenetic vectors for T reflects the contribution of v, that is, the depth of v is 1.

With the advent of genome-scale data, the established reputation of the cophenetic distance and its analysis using various norms has led to increased interest among researchers in utilizing this metric for large-scale biological analyses. Such analyses include phylogenetics [19], median tree estimation [20,21], comparison of evolutionary process models [22], balancing of trees and networks [6,7], studies on theoretical properties of cophenetic $[eqn]$ norm [23], and representative sampling [10]. Calculating cophenetic distances using the naïve algorithm, which has a quadratic time complexity, is slow and impractical for large datasets. Fortunately, an efficient algorithmic framework now enables near-linear time computations for a generalized class of cophenetic distances for the L1 and L2 norms [24]. This advancement makes it feasible to conduct large-scale analyses. However, despite these advancements, algorithms for calculating cophenetic distances in sub-quadratic time for other Lp norms are still unknown.

In this paper, we present a near-linear time algorithmic framework for calculating the Lp norms for the generalized class of cophenetic distances. We demonstrate the performance of our algorithmic framework through a scalability study. Furthermore, we analyze the distributions of key representatives from the generalized class of cophenetic distances across various Lp norms. An implementation of the algorithmic framework can be found on GitHub [25].

Related work. The cophenetic distance is frequently used in phylogenetics to measure the similarity between two trees; both trees are encoded as vectors and are then compared in the corresponding vector space using different norms, such as the L1, L2, or $[eqn]$ norm. Another popular metric that utilizes vector encodings is the path-difference distance, and both metrics can be naïvely computed in O(pn^2^) time for a pair of trees with n taxa under the Lp norm. This quadratic runtime is highly restrictive in the realm of large-scale phylogenetic analysis. [26,27] addressed the quadratic barrier for the path-difference metric under the Lp norm by proposing near-linear time algorithms. For the cophenetic distances, [24] proposed a novel algorithmic framework to efficiently compute the distance between a pair of trees under the L1 and L2 norms in $[eqn]$ time and $[eqn]$ time, respectively [24]. Furthermore, [24] demonstrated that the framework can be applied to a broad class of cophenetic metrics constructed using path-monotonic mappings, which are monotonic on the paths of the input trees. For instance, the original cophenetic distance is derived from the depth of a node, which is a path-monotonic mapping. This condition enables the development of hybrid approaches, resulting in new cophenetic costs rather than metrics, by allowing distinct path-monotonic mappings for each input tree. Such an approach may be particularly suitable when the input trees originate from different sources.

Contribution. We introduce a near-linear time algorithmic framework for computing the Lp norms of a broad class of cophenetic distances. Our framework achieves the following time complexities: $[eqn]$ for even values of p, $[eqn]$ for odd values of p, and $[eqn]$ under the $[eqn]$ norm.

Our approach begins by partitioning the input trees into four approximately equal-sized subtrees. Based on the placement of taxon pairs within these subtrees, we identify several categories of taxon pair locations. For each category, we develop a specific method to compute their contributions to the cophenetic distance. The first category requires four recursive calls on controlled (smaller) subtrees extracted from the input trees. Other categories can be computed in linear time. For the final category, the computation depends on the parity of p, resulting in time complexities of $[eqn]$ and O(pn) for odd and even p, respectively. By integrating these methods, we develop a unified divide-and-conquer algorithm for computing cophenetic distances under Lp norms with finite p. Furthermore, we extend this approach to efficiently compute the $[eqn]$ norm of the cophenetic distance.

We demonstrate the efficiency of our algorithmic frameworks through experimental evaluations. We analyze the runtime of the proposed algorithm for tree pairs with varying numbers of taxa in a scalability study under multiple Lp norms. Our findings indicate that our implementation of the divide-and-conquer strategy significantly outperforms the runtime of the best-known quadratic algorithm for pairs of trees with more than 400 taxa across all tested norms $[eqn]$ .

Lastly, we investigated various types of cophenetic distances under different Lp norms, treating them as distributions of pairwise distances. The trees for each distribution were randomly generated using the uniform model and the Yule model [28]. Our analyses reveal that the class of cophenetic distances provides a high degree of diversity under larger Lp norms, which can benefit numerous applications in comparative phylogenetics.

Methods

This section presents algorithms for computing the Lp norms of cophenetic distances. We begin with the necessary definitions and outline the method for identifying a median vertex in a rooted tree. Next, we describe how input trees are partitioned into four approximately equal subtrees using median vertices. Based on the positions of taxon pairs within these subtrees, we classify them into several categories. For each category, we propose distinct algorithms to calculate their contributions to the cophenetic distance. Finally, we introduce a unified algorithm for computing cophenetic distances for all finite norms and the $[eqn]$ norm. We also analyze and present the time and space complexities of the proposed algorithms.

Definitions

We introduce needed notation and terminology partially following [24]. Let $[eqn]$ be a rooted binary tree, and let v and w be vertices in T. The root of T is denoted by $[eqn]$ . The least common ancestor of v and w in T is denoted by $[eqn]$ . We use $[eqn]$ to denote that w lies on the path between v and the root of T. Note that $[eqn]$ means $[eqn]$ and $[eqn]$ . A vertex v is called strictly internal in T, if v is neither a leaf nor the root. For any non-root vertex v, $[eqn]$ and $[eqn]$ denote the parent and sibling of v, respectively. The set of leaves in T is denoted by LT, and the number of leaves by $[eqn]$ . Similarly, $[eqn]$ denotes the set of all leaves reachable from v, and $[eqn]$ represents the size of $[eqn]$ . A weighted tree T is a rooted binary tree with an edge weight function $[eqn]$ .

We write that a function $[eqn]$ is path-monotonic if, for every v and w such that $[eqn]$ , either $[eqn]$ (descending) or for every v and w such that $[eqn]$ , $[eqn]$ (ascending). In this article, we distinguish the following three contribution functions, defined for unweighted trees, where for a vertex v from T its contribution is:

the depth of a v, i.e., the number of edges on the path from v to the root of T [1],the height of the subtree rooted at v, defined as the maximum number of edges in any path from v to a leaf in the subtree,the number of leaves in the subtree of T rooted at v.

Note that the first contribution function is descending, whereas the other two are ascending.

A similar definition applies to weighted depth and height in weighted trees, but instead of counting the number of edges, we sum the weights of edges. See examples in Fig 1. Note that our definition allows for the contribution of the root. Therefore, a tree may also include an additional rooting edge (with an associated weight), whose bottom vertex is the root, if necessary.

Let T and $[eqn]$ be trees with the same set of leaves $[eqn]$ , where the ordering of leaves is fixed. Let $[eqn]$ and $[eqn]$ (both either ascending or descending) be two path-monotonic contribution functions.

The cophenetic vector of T is defined as $[eqn]$ , and similarly for $[eqn]$ . The Lp-cophenetic distance with respect to the contribution functions $[eqn]$ and $[eqn]$ is defined for $[eqn]$ as:

$$[eqn]$$

where $[eqn]$ is the Lp norm. That is, dp is the Lp norm of the difference between the two cophenetic vectors induced by $[eqn]$ and $[eqn]$ .

Formally, in our case, if p is finite,

$$[eqn]$$

otherwise

$$[eqn]$$

It should be clear that dp is a metric [1] as long as both $[eqn]$ and $[eqn]$ are based on the same type of contribution functions (e.g., depth).

Median vertex in a rooted tree

A vertex t of a rooted tree T divides the tree into two parts: the subtree of T rooted at t, denoted by Tt and referred to as the lower tree with respect to t, and the tree T^t^, called the upper tree, which is obtained by replacing Tt with a leaf. The vertex t, as introduced in the next lemma, is called a median vertex and can be computed in O(n) time.

Lemma 1 (The existence of a median vertex; [24]). For every rooted binary tree T of size $[eqn]$ there is a vertex t such that $[eqn]$ and $[eqn]$ .

We will present the proof for Lemma 1, as it was not included in [24].

Proof: Let t be a vertex computed as follows: initialize t as the root of T, and repeatedly update t to its child with the largest subtree size until $[eqn]$ . We show that the resulting vertex t satisfies the condition.

Let $[eqn]$ be the sibling of t and p be the parent of t. We have |Tt| = |t|, $[eqn]$ , and $[eqn]$ . Hence, $[eqn]$ . If $[eqn]$ then $[eqn]$ . Thus, $[eqn]$ , which is a contradiction. Thus, $[eqn]$ . Finally, $[eqn]$ . $[eqn]$

A median vertex is usually non-unique, e.g., a rooted tree $[eqn]$ has two median vertices b and c. Note that the concept of a median vertex is similar to the centroid of a tree, which partitions an unrooted tree into subtrees, each with a size of at most $[eqn]$ . However, in our case, the median vertex divides a rooted binary tree into two subtrees, with the size of one subtree being between $[eqn]$ and $[eqn]$ , as shown in Lemma 1.

Classification of taxon pairs

In the remaining part of the article, we assume that T and $[eqn]$ are two trees that share the same fixed set of leaves. Let t and $[eqn]$ represent fixed median vertices of T and $[eqn]$ , respectively. Additionally, let $[eqn]$ and $[eqn]$ denote the contribution functions of T and $[eqn]$ , respectively. Without loss of generality, we further assume that all contribution functions $[eqn]$ and $[eqn]$ are descending. That is, for any pair of vertices v and w such that $[eqn]$ , it holds that $[eqn]$ .

Given two trees T and $[eqn]$ , we fix arbitrarily one median vertex in each tree. Then, the path connecting the median vertex with the root will be called a median path. We denote by A and B the sets of leaves excluding the median vertex in the upper and lower trees of T, respectively. Similarly, we denote $[eqn]$ and $[eqn]$ for $[eqn]$ . Note that $[eqn]$ .

We now have four possible classes for a pair of leaves $[eqn]$ , where x can be equal to y, from T, depending on their location:

AA — if both leaves are located in the upper tree (i.e., $[eqn]$ ),BB — if $[eqn]$ ,AB — if $[eqn]$ and $[eqn]$ ,BA — if $[eqn]$ and $[eqn]$ .

When considering the pair $[eqn]$ in both trees, there are 16 possible types of locations, denoted in the form $[eqn]$ , where $[eqn]$ and $[eqn]$ . We say that the pair $[eqn]$ has type $[eqn]$ if $[eqn]$ and $[eqn]$ .

Since some of the types due to symmetry represent the same sets of taxon pairs, we introduce categories for the joint representation of types as indicated in Table 1. There are 4 non-mixed categories $[eqn]$ that correspond uniquely to the types $[eqn]$ , $[eqn]$ , $[eqn]$ and $[eqn]$ . Taking symmetry into account, there are 4 single-mixed categories $[eqn]$ , e.g., in S1 the type $[eqn]$ is equivalent to the $[eqn]$ type, since they involve the same taxon pairs, and 2 double-mixed categories: D1 (which includes $[eqn]$ and $[eqn]$ ) and D2 (which includes $[eqn]$ and $[eqn]$ ).

Table 1: Categories of the types of taxon pair sets. The representative type of each category is marked by a star.

To compute $[eqn]$ , it is sufficient to demonstrate how to compute it for one representative type $[eqn]$ from each category, where $[eqn]$ and $[eqn]$ . Specifically, we calculate the partial distances as:

$$[eqn]$$

Then, the cophenetic distance is the $[eqn]$ -th power of the sum of the partial distances computed for each representative of the category. In the next sections, we show algorithms to compute partial distances for each representative type.

The partial distance of non-mixed types

According to Table 1, there are four types of non-mixed taxon pairs: $[eqn]$ , $[eqn]$ , $[eqn]$ , and $[eqn]$ . To calculate the partial distance for each non-mixed type, we begin by contracting the trees T and $[eqn]$ to the set of leaves defined by the corresponding pair. For instance, for the type $[eqn]$ , the trees T and $[eqn]$ are contracted to the set $[eqn]$ .

We then compute the partial distance dp for these contracted trees recursively. See also lines 14-15 in Algorithm 5.

If f(n) represents the complexity of the algorithm for computing the distance between trees of size n, the total time to calculate these four partial distances is $[eqn]$ , where $[eqn]$ , $[eqn]$ , $[eqn]$ , and $[eqn]$ (note that $[eqn]$ .), plus the linear cost of performing the contractions.

The partial distance of double-mixed types

Here, we present algorithms for computing partial distances for the categories D1 and D2, represented by the two double-mixed types $[eqn]$ and $[eqn]$ , respectively.

Double-mixed types AB|A′B′.

The type $[eqn]$ denotes pairs $[eqn]$ , where x belongs to the upper trees and y to the lower trees of both T and $[eqn]$ . In this case, the $[eqn]$ of $[eqn]$ is positioned along the median path in both trees. Furthermore, $[eqn]$ is the same as $[eqn]$ . Similarly, we have $[eqn]$ .

Having this, the partial distance is

$$[eqn]$$

which can be computed in O(pn) time.

Double-mixed types AB|B′A′.

For $[eqn]$ the naïve approach requires $[eqn]$ steps. Below, we show an O(pn) time solution. We begin with the following problem.

Problem 1. Given two sequences of numbers: $[eqn]$ and $[eqn]$ . Compute: $[eqn]$ .

Algorithm 1 Function $[eqn]$ (partially adapted from [24]).

1: Function $[eqn]$ , where g is a median vertex of G

2:     Set $[eqn]$ for every vertex v on the median path of G

Init counters

3:     For every leaf l in X: $[eqn]$

4:     Initialize an empty sequence $[eqn]$ (a list)

5:     For every vertex v on the median path: append $[eqn]$ to

$[eqn]$ , repeated $[eqn]$ times

6:     Return $[eqn]$

Lemma 2. $[eqn]$ from Algorithm 2 computes $[eqn]$ in $[eqn]$ time and space.

Proof: Correctness: Note that the swap in the 5-th line ensures that $[eqn]$ . Let $[eqn]$ be the number of elements from $[eqn]$ that are smaller or equal to $[eqn]$ . In particular we have $[eqn]$ . The algorithm of $[eqn]$ in the main loop computes $[eqn]$ as $[eqn]$ , i.e., it is the sum of the l-th powers of all elements from $[eqn]$ that are smaller or equal to $[eqn]$ . In particular, $[eqn]$ . Then, for a fixed j,

Algorithm 2 Partial distances: type $[eqn]$ .

1: Input: T and $[eqn]$ with median vertices t and $[eqn]$ , resp.

2: Output: $[eqn]$

3: Function $[eqn]$ :

4:     $[eqn]$

5:     If $[eqn]$ Then swap $[eqn]$ and $[eqn]$

6:     Let $[eqn]$ for all $[eqn]$ and $[eqn]$

7:     While $[eqn]$ :

8:        If $[eqn]$ and $[eqn]$

9:        Then i = i + 1; For $[eqn]$

10:        Else j = j + 1; For $[eqn]$

11:     Return $[eqn]$

12: Let $[eqn]$ and $[eqn]$

13: Return $[eqn]$

$$[eqn]$$

Now, the formula in the 11-th line is derived by applying above the following identities $[eqn]$ and $[eqn]$ .

Time and space complexity: The main loop requires $[eqn]$ steps, while computing the formula in the 11-th line requires O(pm) steps plus precomputing all values of binomial coeficients $[eqn]$ for $[eqn]$ , which can be done once in O(p) time. In total, the overall time complexity is $[eqn]$ , the same applies to space complexity. $[eqn]$

Now, $[eqn]$ and similarly, $[eqn]$ is computed by calling $[eqn]$ , where $[eqn]$ is the sequence of contributions of $[eqn]$ ’s for x in $[eqn]$ and $[eqn]$ is the sequence of contributions of $[eqn]$ ’s for y in $[eqn]$ . Such sequences are inferred in O(n) steps by the function $[eqn]$ in Algorithm 1.

Lemma 3. Algorithm 2 computes $[eqn]$ in O(pn) time.

Proof: The proof is similar to the proof of Lemma 3 from [24]. The difference is in the more general computation of the sum from Lemma 2. We leave out simple details. $[eqn]$

The partial distance of single-mixed types

There are four single-mixed categories represented by the types $[eqn]$ , $[eqn]$ , $[eqn]$ , and $[eqn]$ (see Table 1). Each of these variants can be solved similarly. Thus, we present only the algorithm for computing the partial distance for $[eqn]$ .

Assume that the pair of taxa $[eqn]$ has the type $[eqn]$ , meaning $[eqn]$ and $[eqn]$ . In this case, $[eqn]$ is a vertex from the upper part of tree T, while $[eqn]$ lies on the median path of $[eqn]$ . To apply our algorithm, we consider the tree to be ordered. For non-leaf vertices, the right child of a vertex v is denoted $[eqn]$ , and the left child is denoted $[eqn]$ .

Let $[eqn]$ . Our solution is divided depending on the parity of p. We start with the solution to odd norms.

Single-mixed types under odd Lp norms.

We start with the following definitions.

$[eqn]$ and $[eqn]$ .

In the next two lemmas, we prove several properties of vertex attributes $[eqn]$ from Algorithm 3.

Algorithm 3 Partial distances: type $[eqn]$ for odd p’s.

1: Input: T and $[eqn]$ with median verticest and $[eqn]$ , respectively;   p is odd

2: Output: $[eqn]$

3: For every $[eqn]$ in the upper tree of T:

4:     $[eqn]$ ; $[eqn]$ # the zero vector

5: For $[eqn]$ : # The preprocessing loop

6:     $[eqn]$ ; $[eqn]$

7:     If $[eqn]$

8:        Then $[eqn]$ ;

$[eqn]$

9:     If $[eqn]$ Then $[eqn]$ Else $[eqn]$

10: For every non-root $[eqn]$ in the upper tree of T in postfix

order: # The main loop

11:     $[eqn]$ ;

$[eqn]$

12:     $[eqn]$

13:     $[eqn]$

14: Return $[eqn]$

Lemma 4. If v is strictly internal, then after the preprocessing loop of Algorithm 3, we have

$$[eqn]$$

where $[eqn]$ .

Proof: In line 8 of Algorithm 2, $[eqn]$ is defined as the vertex such that $[eqn]$ , and this inequality is not satisfied by the parent of $[eqn]$ ; that is, $[eqn]$ . Combining these observations, we can conclude that $[eqn]$ consists of vertices $[eqn]$ for which $[eqn]$ as specified in line 8. Consequently, from the assignment in line 8, we have

$[eqn]$ . $[eqn]$

Lemma 5. If v is strictly internal, then after the main loop of Algorithm 3,

$[eqn]$ ,and $[eqn]$ .

Proof: We present an inductive proof for a fixed $[eqn]$ and for $[eqn]$ , which denotes the $[eqn]$ -th element of the vector $[eqn]$ . If v is a leaf, then the equality follows directly from line 9 of Algorithm 3 and the definition of $[eqn]$ . Now, let us consider the case where v is strictly internal. Note that $[eqn]$ .

$$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$

The proof for $[eqn]$ is similar. We omit details. $[eqn]$

Lemma 6. If v is strictly internal in T, then after the main loop of Algorithm 3 we have

$$[eqn]$$

Proof: We can prove (2) by using Lemmas 4 and 5. Let v be strictly internal in the upper tree of T. Then, for a leaf $[eqn]$ and a leaf y from $[eqn]$ (i.e., from $[eqn]$ , see line 12), we have $[eqn]$ and $[eqn]$ . Let R be the right side of (2). Then,

$$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$

The last equation is obtained by identities from Lemma 5 and binomial expansions. R equals the right side of the assignment from line 13. $[eqn]$

Lemma 7. Algorithm 3 computes $[eqn]$ in $[eqn]$ time.

Proof: We show the algorithm’s correctness followed by the stated time complexity.

Correctness: Let I be the set of non-root vertices from the upper tree of T. Then, every pair of leaves $[eqn]$ of type $[eqn]$ uniquely determines $[eqn]$ such that $[eqn]$ and $[eqn]$ . It also follows that $[eqn]$ . Let $[eqn]$ denote such vertex v. For a given v, $[eqn]$ is the set of all pairs $[eqn]$ such that $[eqn]$ . Hence,

$$[eqn]$$

By Lemma 6 the above sum equals the value returned in the last line of Algorithm 3.

Time complexity: The key aspect of the time complexity is found in line 8. We show that $[eqn]$ can be found by a binary search in $[eqn]$ time that seeks the value in an ordered array composed of vertices on the path connecting a given leaf x with the root of T. An infix traversal of T can construct such an array. Then, a vertex is inserted into the array when visited for the first time. When a vertex is visited for the last time, it is removed from the array. Due to the monotonic ordering of paths, the array is always sorted, and its size is limited by n. The time complexity of the remaining loops is O(pn), which gives $[eqn]$ time of Algorithm 3. $[eqn]$

Single-mixed types under even Lp norms.

Algorithm 4 Partial distances: type $[eqn]$ for even p.

1: Input: T and $[eqn]$ with median vertices t and $[eqn]$ , resp.

2: Output: $[eqn]$

3: For every v in the upper tree of T: $[eqn]$ . $[eqn]$

(the zero vector)

4: For $[eqn]$ : # Preprocessing

5:     $[eqn]$ $[eqn]$

6: For every non-root v in the upper tree of T in postfix

order:

7:     $[eqn]$

8:     $[eqn]$

9: Return $[eqn]$

The advantage of even Lp norms is given by the relation $[eqn]$ for p even. This fact allows us to circumvent the extra complexity of Algorithm 3. Algorithm 4 outlines the computation of $[eqn]$ when p is even.

Lemma 8. For each non-root vertex v in the upper tree of T after Algorithm 4 we have

$[eqn]$ , $[eqn]$ .

Proof: The relation for $[eqn]$ follows directly from the bottom-up nature of the algorithm and the relation for $[eqn]$ can be observed using the following equation:

$$[eqn]$$

$[eqn]$

Lemma 9. Algorithm 4 computes $[eqn]$ in O(pn) time for even p.

Proof: The correctness of the algorithm follows from Lemmas 8 and 7. Finally, the time complexity of Algorithm 4 is O(pn) since each line within the loops requires O(p) time and is executed at most 2n times. $[eqn]$

Partial distances of the remaining single-mixed types.

As previously mentioned, Algorithm 3 and Algorithm 4 for computing the partial distances of single mixed type $[eqn]$ can be adapted to solve the other types as follows. For the type $[eqn]$ , replace the term “upper” with “lower” and A with B in both algorithms. For the type $[eqn]$ , swap the input trees and execute both algorithms. For the type $[eqn]$ , swap the input trees and run the algorithm designed for the type $[eqn]$ . These modifications are detailed in lines 14 through 16 of Algorithm 5.

Algorithm to compute the cophenetic distance under the Lp norm

The pseudo-code in Algorithm 5 summarizes the complete procedure for computing the cophenetic distance for finite p. The correctness of the algorithm follows from the results presented in the previous section. Below, we analyze the time complexity in two scenarios: when p is constant, and when p is considered as a parameter in the asymptotic analysis.

Algorithm 5 Computing cophenetic distance.

1: Input: T and $[eqn]$ with the same set of leaves of size n, two

ascending contribution functions $[eqn]$ and $[eqn]$ ,

and an integer $[eqn]$

2: Output: The Lp-cophenetic distance between T and $[eqn]$ with

respect to $[eqn]$ and $[eqn]$

3: Compute all binomial coeficients $[eqn]$ , for all i, by

$[eqn]$ and a0 = 1

4: Function $[eqn]$

5:     Compute t and $[eqn]$ the median vertices of T and $[eqn]$ ,

resp.

6:     Let A and $[eqn]$ be the set of all leaves in the upper

tree of T and $[eqn]$ , resp.

7:     $[eqn]$ ; $[eqn]$ ; $[eqn]$

8:     For X, in $[eqn]$ : For $[eqn]$ in $[eqn]$

9:        Let $[eqn]$ and $[eqn]$

10:        $[eqn]$ # Compute partial distances

for type $[eqn]$

11:     $[eqn]$ see Eq. (1) # *Double-mixed $[eqn]$ *

12:     $[eqn]$ by Alg. 2 # *Double-mixed $[eqn]$ *

13:     # All single-mixed types

14:     $[eqn]$ using Alg. 3 if p is odd, or Alg. 4

otherwise

15:     Similarly to line 14 , $[eqn]$ , but replace

“upper” with “lower” and A

with B in Alg. 3/4

16:     For $[eqn]$ and $[eqn]$ swap the input trees and

repeat lines 14 and 15

17:     Return s # Return the sum of partial distances

18: Return $[eqn]$ # Return the cophenetic distance

First, consider the case where p is constant. Since $[eqn]$ , it follows from Lemmas 3 and 7 that computing the partial distances of mixed types requires $[eqn]$ time when p is odd, and O(n) time when p is even. Consequently, the overall time complexity of Algorithm 5, as derived in [24], is $[eqn]$ for odd p, and $[eqn]$ for even p.

If p is a parameter, $[eqn]$ is the worst-case time complexity of the complete algorithm. Then, computation of partial distances of non-mixed types (see lines 8-10 in Algoritm 5) requires $[eqn]$ time where $[eqn]$ and $[eqn]$ for each i (by Lemma 1), while the computation of mixed types requires $[eqn]$ time if p is odd or O(pn) time if p is even, (Lemmas 3 and 7). Therefore, for some $[eqn]$ we can write that, $[eqn]$ if $[eqn]$ , $[eqn]$ $[eqn]$ , if n>5 and p is odd, and $[eqn]$ $[eqn]$ , otherwise.

Theorem 10. The time complexity of the algorithm is $[eqn]$ , when p is odd and $[eqn]$ when p is even.

Proof: The proof for odd norms. Assume that p is odd. We show that there are constants $[eqn]$ and d>0 such that for every n>0 and p>0, $[eqn]$ . The proof is by induction on n. For $[eqn]$ we have $[eqn]$ and the inequality is satisfied. For n>5, $[eqn]$ .

Let $[eqn]$ . Then, for n>5, $[eqn]$ . Also, $[eqn]$ and $[eqn]$ $[eqn]$ . Finally, for n>5, $[eqn]$ .

Even norms. If p is even, it suffices to show that $[eqn]$ , for some constants d>0 and $[eqn]$ . The proof is similar and simpler compared to the odd case and follows the same reasoning. We omit details. $[eqn]$

The last theorem shows that, despite the higher complexity of odd norms, the factor p only appears in the asymptotically minor term of $[eqn]$ . This suggests that the computational effort for cophenetic distances is more influenced by the size of the input trees than by the chosen norm level p. Additionally, this theoretical result highlights that the binary search algorithm - required for computing double-mixed types and representing the main distinction between the two algorithmic variants - introduces $[eqn]$ more steps in the worst-case than the algorithm for even norms. See the experimental section for a detailed discussion on scalability and a comparison of the algorithms for odd and even norms.

Cophenetic distances under L∞ norm.

Algorithm 5 can be easily adapted to solve the last remaining case of cophenetic distances. The algorithm for the $[eqn]$ norm is similar to the case p = 1. The only difference is taking the maximum value instead of adding partial distances as follows. First, fix $[eqn]$ . Then in line 7, replace $[eqn]$ with $[eqn]$ . In lines 10-16 instead of assignments s = s + r, where r is the right-hand expression, write $[eqn]$ . Now, the time complexity is $[eqn]$ , which follows from Theorem 10 for the odd case p = 1.

Results

In the following sections, we present the scalability study and distribution analysis results. The near-linear time algorithm, i.e., Algorithm 5, for computing the Lp norm cophenetic distance and the naïve (quadratic time) algorithm were implemented in Rust 1.79.0. The implementation of both algorithms, along with instructions on reproducing the simulation results below, is available on GitHub [25]. Note that the distances computed by the naïve and our algorithm are the same.

Scalability analysis

We investigate the scalability of our algorithm by comparing its runtime to that of the best-known quadratic-time algorithm for computing cophenetic distance.

Datasets. We generated three datasets consisting of pairs of random binary trees with n leaves using the classic Yule model, where for each n we generated q pairs. The first dataset was generated using $[eqn]$ and q = 1000, the second with $[eqn]$ and q = 1000, and the last one with $[eqn]$ and q = 100. Algorithms were executed on each pair of trees for various norms Lp with p values ranging from 1 to 100. The comparative runtime analysis for each dataset is illustrated in the diagrams of Figs 2, 3, and 4.

Scalability analysis comparing the near-linear algorithm and the naïve quadratic-time algorithm.The diagrams show the average runtime with standard deviation bands over 1000 runs for each tree size (n=25,50,…,600) on pairs of random trees, evaluated for Lp-norms (p=1,2,5,10,20,50,100). Crossing points indicate the tree sizes where the near-linear algorithm outperforms the naïve one.

Scalability analysis comparing the near-linear algorithm and the naïve quadratic-time algorithm for small norms.The diagrams show the average runtime with standard deviation bands over 1,000 runs for each tree size (n=200,400,…,2600) on pairs of random trees, evaluated for Lp-norms (p=1,2,7,8).

Scalability analysis comparing the near-linear algorithm (outlined bars) and the naïve quadratic-time algorithm on large trees.The diagrams show the average runtime with standard deviation error over 100 runs for each tree size (n=1000,2000,…,6000) on pairs of random trees, evaluated for Lp-norms (p=1,2,5,10,20,50,100).

Results. Our algorithm shows a significant improvement in efficiency compared to the quadratic solution for trees with more than 250 taxa and norms below 20, as illustrated in Fig 2. Furthermore, the crossing point at which our algorithm outperforms the quadratic solution increases gradually as p increases, reaching 366 at p = 100. This finding indicates the potential for developing a hybrid algorithm that utilizes the divide-and-conquer method for larger trees while reverting to the quadratic time algorithm for trees with fewer than 230–370 leaves during recursion.

Despite the differences in asymptotic worst-case time complexity between odd and even norms, these differences are not evident in the average runtime diagrams presented in Figs 3 and 4. Many of these differences can be attributed to the binary search step in Algorithm 3, as explained in Lemma 7, where the runtime is relatively straightforward to estimate. Notably, the computation of odd norms accounted for only 0.01% of the total average runtime in the binary search steps. In practical terms, this indicates that the time complexity of Algorithm 3 is closer to O(pn) rather than the more conservative estimate of $[eqn]$ provided in Lemma 7. Consequently, the overall runtime complexity for calculating the cophenetic distance can be approximated as $[eqn]$ for any given norm. Since the trees were generated randomly, we propose that $[eqn]$ is also a valid estimator of the average time complexity of our algorithm. However, as shown in Fig 4, there are no significant differences in runtime across various norms for a fixed tree size n. This indicates that the runtime is primarily affected by the overhead associated with maintaining supporting data structures, rather than by the computation of values, which involves loops over the range from 0 to p. We conjecture that for extremely large values of p, the computational cost of these loops will become increasingly noticeable in the runtime; however, such scenarios were not tested in this study.

Cophenetic distance distributions

We investigate the distributions of three key representatives from the generalized class of cophenetic distances. These representatives are (i) the original cophenetic distance metric defined by depths, (ii) the metric defined by the heights of the subtrees, and (iii) the metric defined by the number of taxa in the subtrees. For ease of reference, we refer to these metrics as the depth, subtree-height, and subtree-size cophenetic metrics, respectively.

To the best of our knowledge, the specific distributions for any selected metrics remain unpublished. Therefore, we present here the sampled distributions of cophenetic distances based on two standard models of phylogenetic tree sampling: the uniform model and the Yule model [28,29]. The sampled distributions for the classical depth-based cophenetic distance were previously discussed in [1] for both the L1 norm and the L2 norm under the uniform model.

Data. We generated a dataset containing 10^6^ pairs of trees, each with 100 leaves. Each tree was independently generated under the uniform model. Similarly, we generated another dataset based on the Yule model. It is important to note that the generated trees do not include edge lengths; therefore, both depths and subtree heights are measured in terms of the number of edges.

Results. The sampled distributions are illustrated in Figs 5 and 6, with the corresponding mean and standard deviation statistics depicted in Fig 7.

Sampled distribution of three cophenetic distances under the Lp norms for p = 1 to 10 (rows).Pairs of trees of size 100 were sampled according to the uniform model (left columns) and the Yule model (right columns). The frequencies are grouped into 200 bins. Note that the y-axis scale and range are the same across all diagrams. To enhance visibility, low frequencies (e.g., at low distances close to 0) are omitted and the width of the diagrams is appropriately adjusted.

(Cont. from Fig 5). Sampled distribution of three cophenetic distances is shown for Lp norms with p = 20, 50, 100, and ∞ (rows). Similar to Fig 5, the frequencies for finite p are grouped into 200 bins. Since all three contribution functions return integers, the L∞ norms are also integers. As a result, the frequency values for L∞ are presented without binning.

Mean and standard deviation from the sampled distributions for the uniform and Yule models for Lp norms with p>1.

The distribution of subtree-height and depth cophenetic metrics exhibit similar shapes under the uniform model for smaller values of p as illustrated in Fig 5. A similar trend is observed for the depth under the Yule model. Furthermore, these distributions are positively skewed, which is consistent with the findings for the L1 and L2 norms of the depth cophenetic metric reported in [1]. On the contrary, the subtree-height distributions exhibit positive skewness. It is also worth noting that our study involved a significantly larger number of pairs for sampling. As a result, the diagrams in Fig 5 appear smoother compared to the sampled distributions of depth metric from [1].

Since cophenetic vectors are finite, it is straightforward to prove that $[eqn]$ under any contribution function. Consequently, Lp cophenetic distributions converge to $[eqn]$ distributions as p tends to infinity. Furthermore, in our case, all three contribution functions return integer values, which implies that the $[eqn]$ distributions have an integer domain, as illustrated in Fig 6. This property is evident in diagrams, where the distributions become increasingly irregular as p grows and more and more similar to the corresponding discrete distributions under $[eqn]$ . This effect is particularly noticeable for the depth and subtree-height metrics. By contrast, the distributions for the subtree-size metric remain relatively smooth. This smoothness arises from the fact that, under the $[eqn]$ norm, almost all frequencies under uniform and Yule models are concentrated at 98 (see low standard deviation in Fig 7). This value is derived from the pair of leaf labels that form a cherry in one tree but are separated by the root in the other tree, yielding subtree sizes of 2 and 100, respectively, and a corresponding distance of 98.

The sampled distributions for the subtree-size cophenetic metric are notably distinct from the analysis mentioned above. The histograms for this metric are negatively-skewed, and the mean value under the Yule model is, in this case, larger than the mean value under the uniform model.

In comparing the depth and subtree-height cophenetic metrics, we find that the mean values under the Yule model are significantly lower than those under the uniform model. A similar bias towards the Yule model was previously noted in the sampled distributions for the path-difference distance [30,31].

Additionally, the mean and standard deviation of the depth metric under the L1 and L2 norms align with the exact values reported in [32].

Results and discussion

We introduced a novel algorithmic framework for computing the Lp norm cophenetic distance, achieving a time complexity of $[eqn]$ for odd values of p, $[eqn]$ for even values of p, and $[eqn]$ for $[eqn]$ norm. This represents a substantial advancement compared to the previously best-known naïve algorithm, which requires $[eqn]$ time.

Additionally, our scalability studies suggest that the estimated runtime of our algorithm approaches $[eqn]$ time under all Lp norms with finite p, contrasting with the larger upper asymptotic bound observed for odd values of p. These advancements greatly improve the usability of the cophenetic distance for large-scale phylogenetic studies and the median-tree inference of species trees from gene trees using this metric.

Distribution analyses of these three key representative metrics from the cophenetic class further enhance this work, offering practitioners valuable guidance in selecting appropriate metrics for their specific needs.

The framework demonstrates broad practical applicability by generalizing to all metrics that rely on the path-monotonic property, here referred to as the class of generalized cophenetic distances. This generalization can be achieved by either designing contribution functions that satisfy basic monotonicity properties or by forming linear, positively weighted combinations of existing contribution functions.

As a result, the class of generalized cophenetic distances includes the original cophenetic distance based on depth, as well as other metrics, particularly the subtree height and subtree size cophenetic distances. Additionally, this class incorporates the more recently proposed metric from [19], which combines both weighted and unweighted contribution functions. Consequently, this metric can be computed in near-linear time using our algorithm under any Lp norm.

Furthermore, the framework can be applied to mixed scenarios where different types of cophenetic vectors are used for the trees — for example, depth in one tree and height in another. Although these mixed scenarios may not fully satisfy metric properties, they can still be useful for comparing trees from different origins by allowing asymmetry. For instance, this approach may be more suitable when comparing a gene tree, which represents the evolutionary history of a gene, to a species tree that illustrates the evolutionary history of the species from which the genes were sampled.