arXiv:1108.5781v1 [math.PR] 29 Aug 2011

Phase Transition in Distance-Based
Phylogeny Reconstruction∗
Sébastien Roch†
January 18, 2018

Abstract
We introduce a new distance-based phylogeny reconstruction technique
which provably achieves, at sufficiently short branch lengths, a logarithmic
sequence-length requirement—improving significantly over previous polynomial bounds for distance-based methods and matching existing results for
general methods. The technique is based on an averaging procedure that
implicitly reconstructs ancestral sequences.
In the same token, we extend previous results on phase transitions in phylogeny reconstruction to general time-reversible models. More precisely, we
show that in the so-called Kesten-Stigum zone (roughly, a region of the parameter space where ancestral sequences are well approximated by “linear
combinations” of the observed sequences) sequences of length O(log n) suffice for reconstruction when branch lengths are discretized. Here n is the
number of extant species.
Our results challenge, to some extent, the conventional wisdom that estimates of evolutionary distances alone carry significantly less information
about phylogenies than full sequence datasets.

Keywords: Phylogenetics, distance-based methods, phase transitions, reconstruction problem.

∗

The results detailed here were announced without proof in [Roc08, Roc10].
Department of Mathematics and Bioinformatics Program, University of California-Los Angeles,
Los Angeles, California 90095, USA. Work supported by NSF grant DMS-1007144.
†

1

Introduction

The evolutionary history of a group of organisms is generally represented by a
phylogenetic tree or phylogeny [Fel04, SS03]. The leaves of the tree represent the
current species. Each branching indicates a speciation event. Many of the most
popular techniques for reconstructing phylogenies from molecular data, e.g. UPGMA, Neighbor-Joining, and BIO-NJ [SS63, SN87, Gas97], are examples of what
are known as distance-matrix methods. The main advantage of these methods is
their speed, which stems from a straightforward approach: 1) the estimation of a
distance matrix from observed molecular sequences; and 2) the repeated agglomeration of the closest clusters of species. Each entry of the distance matrix is an
estimate of the evolutionary distance between the corresponding pair of species,
that is, roughly the time elapsed since their most recent common ancestor. This
estimate is typically obtained by comparing aligned homologous DNA sequences
extracted from the extant species—the basic insight being, the closer the species,
the more similar their sequences. Most distance methods run in time polynomial
in n, the number of leaves, and in k, the sequence length. This performance compares very favorably to that of the other two main classes of reconstruction methods, likelihood and parsimony methods, which are known to be computationally
intractable [GF82, DS86, Day87, MV05, CT06, Roc06].
The question we address in this paper is the following: Is there a price to
pay for this speed and simplicity? There are strong combinatorial [SHP88] and
statistical [Fel04] reasons to believe that distance methods are not as accurate as
more elaborate reconstruction techniques, notably maximum likelihood estimation
(MLE). Indeed, in a typical instance of the phylogenetic reconstruction problem,
we are given aligned DNA sequences {(ξli )ki=1 }l∈L , one sequence for each leaf
l ∈ L, from which we seek to infer the phylogeny on L. Generally, all sites
(ξli )l∈L , for i = 1, . . . , k, are assumed to be independent and identically distributed
according to a Markov model on a tree (see Section 1.1). For a subset W ⊆ L,
we denote by µW the distribution of (ξli )l∈W under this model. Through their
use of the distance matrix, distance methods reduce the data to pairwise sequence
correlations, that is, they only use estimates of µ2 = {µW : W ⊆ L, |W | = 2}.
In doing so, they seemingly fail to take into account more subtle patterns in the data
involving three or more species at a time. In contrast, MLE for example outputs
a model that maximizes the joint probability of all observed sequences. We call
methods that explicitly use the full dataset, such as MLE, holistic methods.
It is important to note that the issue is not one of consistency: when the sequence length tends to infinity, the estimate provided by distance methods—just
like MLE—typically converges to the correct phylogeny. In particular, under mild
assumptions, it suffices to know the pairwise site distributions µ2 to recover the
2

topology of the phylogeny [CH91, Cha96]. Rather the question is: how fast is
this convergence? Or more precisely, how should k scale as a function of n to
guarantee a correct reconstruction with high probability? And are distance methods significantly slower to converge than holistic methods? Although we do not
give a complete answer to these questions of practical interest here, we do provide
strong evidence that some of the suspicions against distance methods are based on
a simplistic view of the distance matrix. In particular, we open up the surprising
possibility that distance methods actually exhibit optimal convergence rates.
Context. It is well-known that some of the most popular distance-matrix methods actually suffer from a prohibitive sequence-length requirement [Att99, LC06].
Nevertheless, over the past decade, much progress has been made in the design
of fast-converging distance-matrix techniques, starting with the seminal work of
Erdös et al. [ESSW99a]. The key insight behind the algorithm in [ESSW99a], often dubbed the Short Quartet Method (SQM), is that it discards long evolutionary
distances, which are known to be statistically unreliable. The algorithm works by
first building subtrees of small diameter and, in a second stage, putting the pieces
back together. The SQM algorithm runs in polynomial time and guarantees the
correct reconstruction with high probability of any phylogeny (modulo reasonable
assumptions) when k = poly(n). This is currently the best known convergence rate
for distance methods. (See also [DMR06, DHJ+ 06, Mos07, GMS08, DMR09] for
faster-converging algorithms involving partial reconstruction of the phylogeny.)
Although little is known about the sequence-length requirement of MLE [SS99,
SS02], recent results of Mossel [Mos04], Daskalakis et al. [DMR06, DMR11], and
Mihaescu et al. [MHR09] on a conjecture of Steel [Ste01] indicate that convergence rates as low as k = O(log n) can be achieved when the branch lengths are
sufficiently short, using insights from statistical physics. We briefly describe these
results.
As mentioned above, the classical model of DNA sequence evolution is a
Markov model on a tree that is closely related to stochastic models used to study
particle systems [Lig85, Geo88]. This type of model undergoes a phase transition that has been extensively studied in probability theory and statistical physics:
at short branch lengths (in the binary symmetric case, up to 15% divergence per
edge), in what is called the reconstruction phase, good estimates of the ancestral sequences can be obtained from the observed sequences; on the other hand,
outside the reconstruction phase, very little information about ancestral states diffuses to the leaves. See e.g. [EKPS00] and references therein. The new algorithms
in [Mos04, DMR06, DMR11, MHR09] exploit this phenomenon by alternately 1)
reconstructing a few levels of the tree using distance-matrix techniques and 2) es-

3

timating distances between internal nodes by reconstructing ancestral sequences
at the newly uncovered nodes. The overall algorithm is not distance-based, however, as the ancestral sequence reconstruction is performed using a complex function of the observed sequences named recursive majority. The rate k = O(log n)
achieved by these algorithms is known to be necessary in general. Moreover, the
slower rate k = poly(n) is in fact necessary for all methods—distance-based or
holistic—outside the reconstruction phase [Mos03]. In particular, note that distance methods are in some sense “optimal” outside the reconstruction phase by the
results of [ESSW99a].
Beyond the oracle view of the distance matrix. It is an outstanding open problem to determine whether distance methods can achieve k = O(log n) in the reconstruction phase1 . From previous work on fast-converging distance methods, it
is tempting to conjecture that k = poly(n) is the best one can hope for. Indeed,
all previous algorithms use the following “oracle view” of the distance matrix, as
formalized by King et al. [KZZ03] and Mossel [Mos07]. As mentioned above, the
reliability of distance estimates depends on the true evolutionary distances. From
standard concentration inequalities, it follows that if leaves a and b are at distance
τ (a, b), then the usual distance estimate τ̂ (a, b) (see Section 1.1) satisfies:
if τ (a, b) < D + ε or τ̂ (a, b) < D + ε then |τ (a, b) − τ̂ (a, b)| < ε,

(1)

for ε, D such that k ∝ (1 − e−ε )−2 e2D . Fix ε > 0 small and k  poly(n). Let
T be a complete binary tree with log2 n levels. Imagine that the distance matrix
is given by the following oracle: on input a pair of leaves (a, b) the oracle returns
an estimate τ̂ (a, b) which satisfies (1). Now, notice that for any tree T 0 which is
identical to T on the first log2 n/2 levels above the leaves, the oracle is allowed to
return the same distance estimate as for T . That is, we cannot distinguish T and T 0
in this model unless k = poly(n). (This argument can be made more formal along
the lines of [KZZ03].)
What the oracle model ignores is that, under the assumption that the sequences
are generated by a Markov model of evolution, the distance estimates
(τ̂ (a, b))a,b∈[n]
are in fact correlated random variables. More concretely, for leaves a, b, c, d, note
that the joint distribution of (τ̂ (a, b), τ̂ (c, d)) depends in a nontrivial way on the
joint site distribution µW at W = {a, b, c, d}. In other words, even though the
distance matrix is—seemingly—only an estimate of the pairwise correlations µ2 ,
1

Mike Steel offers a 100$ reward for the solution of this problem.

4

it actually contains some information about all joint distributions. Note however
that it is not immediately clear how to exploit this extra information or even how
useful it could be.
As it turns out, the correlation structure of the distance matrix is in fact very
informative at short branch lengths. More precisely, we introduce in this paper a
new distance-based method with a convergence rate of k = O(log n) in the reconstruction phase (to be more accurate, in the so-called Kesten-Stigum phase; see
below)—improving significantly over previous poly(n) results. Note that the oracle model allows only the reconstruction of a o(1) fraction of the levels in that
case. Our new algorithm involves a distance averaging procedure that implicitly
reconstructs ancestral sequences, thereby taking advantage of the phase transition
discussed above. We also obtain the first results on Steel’s conjecture beyond the
simple symmetric models studied by Daskalakis et al. [DMR06, DMR11, MHR09]
(the so-called CFN and Jukes-Cantor models). In the next subsections, we introduce general definitions and state our results more formally. We also give an
overview of the proof.
Further related work. For further related work on efficient phylogenetic tree
reconstruction, see [ESSW99b, HNW99, CK01, Csu02].

1.1

Definitions

Phylogenies. We define phylogenies and evolutionary distances more formally.
Definition 1 (Phylogeny) A phylogeny is a rooted, edge-weighted, leaf-labeled
tree T = (V, E, [n], ρ; τ ) where: V is the set of vertices; E is the set of edges;
L = [n] = {0, . . . , n − 1} is the set of leaves; ρ is the root; τ : E → (0, +∞) is a
positive edge weight function. We further assume that all internal nodes in T have
degree 3 except for the root ρ which has degree 2. We let Yn be the set of all such
phylogenies on n leaves and we denote Y = {Yn }n≥1 .
Definition 2 (Tree Metric) For two leaves a, b ∈ [n], we denote by Path(a, b) the
set of edges on the unique path between a and b. A tree metric on a set [n] is a
positive function d : [n] × [n] → (0, +∞) such that there exists a tree T = (V, E)
with leaf set [n] and an edge weight function w : E → (0, +∞) satisfying the
following: for all leaves a, b ∈ [n]
X
d(a, b) =
we .
e∈Path(a,b)

5

For convenience, we denote by (τ (a, b))a,b∈[n] the tree metric corresponding to the
phylogeny T = (V, E, [n], ρ; τ ). We extend τ (u, v) to all vertices u, v ∈ V in the
obvious way.
Example 1 (Homogeneous Tree) For an integer h ≥ 0, we denote by T (h) =
(V (h) , E (h) , L(h) , ρ(h) ; τ ) a rooted phylogeny where T (h) is the h-level complete
binary tree with arbitrary edge weight function τ and L(h) = [2h ]. For 0 ≤ h0 ≤ h,
(h)
(h)
we let Lh0 be the vertices on level h−h0 (from the root). In particular, L0 = L(h)
(h)
and Lh = {ρ(h) }. We let HY = {HYn }n≥1 be the set of all phylogenies with
homogeneous underlying trees.
Model of molecular sequence evolution. Phylogenies are reconstructed from
molecular sequences extracted from the observed species. The standard model of
evolution for such sequences is a Markov model on a tree (MMT).
Definition 3 (Markov Model on a Tree) Let Φ be a finite set of character states
with ϕ = |Φ|. Typically Φ = {+1, −1} or Φ = {A, G, C, T}. Let n ≥ 1 and let
T = (V, E, [n], ρ) be a rooted tree with leaves labeled in [n]. For each edge e ∈ E,
we are given a ϕ × ϕ stochastic matrix M e = (Mije )i,j∈Φ , with fixed stationary
distribution π = (πi )i∈Φ . An MMT ({M e }e∈E , T ) associates a state σv in Φ to
each vertex v in V as follows: pick a state for the root ρ according to π; moving
away from the root, choose a state for each vertex v independently according to the
distribution (Mσeu ,j )j∈Φ , with e = (u, v) where u is the parent of v.
The most common MMT used in phylogenetics is the so-called general timereversible (GTR) model.
Definition 4 (GTR Model) Let Φ be a set of character states with ϕ = |Φ| and
π be a distribution on Φ satisfying πi > 0 for all i ∈ Φ. For n ≥ 1, let T =
(V, E, [n], ρ; τ ) be a phylogeny. Let Q be a ϕ × ϕ rate matrix, that is, Qij > 0 for
all i 6= j and
X
Qij = 0,
j∈Φ

for all i ∈ Φ. Assume Q is reversible with respect to π, that is,
πi Qij = πj Qji ,
for all i, j ∈ Φ. The GTR model on T with rate matrix Q is an MMT on T =
(V, E, [n], ρ) with transition matrices M e = eτe Q , for all e ∈ E. By the reversibility assumption, Q has ϕ real eigenvalues
0 = Λ1 > Λ2 ≥ · · · ≥ Λϕ .
6

We normalize Q by fixing Λ2 = −1. We denote by Qϕ the set of all such rate
matrices. We let Gn,ϕ = Yn ⊗ Qϕ be the set of all ϕ-state GTR models on n
leaves. We denote Gϕ = {Gn,ϕ }n≥1 . We denote by ξW the vector of states on the
vertices W ⊆ V . In particular, ξ[n] are the states at the leaves. We denote by LT ,Q
the distribution of ξ[n] .
GTR models include as special cases many popular models such as the CFN model.
Example 2 (CFN Model) The CFN model is the GTR model with ϕ = 2, π =
(1/2, 1/2), and


−1/2 1/2
CFN
Q=Q
≡
.
1/2 −1/2
Example 3 (Binary Asymmetric Channel) More generally, letting Φ = {+, −}
and π = (π+ , π− ), with π+ , π− > 0, we can take


−π− π−
Q=
.
π+ −π+
Phylogenetic reconstruction. A standard assumption in molecular evolution is
that each site in a sequence (DNA, protein, etc.) evolves independently according to a Markov model on a tree, such as the GTR model above. Because of the
reversibility assumption, the root of the phylogeny cannot be identified and we
reconstruct phylogenies up to their root.
e = {Y
e n }n≥1 be a
Definition 5 (Phylogenetic Reconstruction Problem) Let Y
e ϕ be a subset of rate matrices on ϕ states. Let T =
subset of phylogenies and Q
e If T = (V, E, [n], ρ) is the rooted tree underlying T , we
(V, E, [n], ρ; τ ) ∈ Y.
denote by T− [T ] the tree T where the root is removed: that is, we replace the
two edges adjacent to the root by a single edge. We denote by Tn the set of all
leaf-labeled trees on n leaves with internal degrees 3 and we let T = {Tn }n≥1 .
A phylogenetic reconstruction algorithm is a collection of maps A = {An,k }n,k≥1
i )k
[n] k
from sequences (ξ[n]
i=1 ∈ (Φ ) to leaf-labeled trees T ∈ Tn . We only consider
algorithms A computable in time polynomial in n and k. Let k(n) be an increasing function of n. We say that A solves the phylogenetic reconstruction problem
e ⊗Q
e ϕ with sequence length k = k(n) if for all δ > 0, there is n0 ≥ 1 such
on Y
e n, Q ∈ Q
e ϕ,
that for all n ≥ n0 , T ∈ Y
h


i
i k(n)
P An,k(n) (ξ[n]
)i=1 = T− [T ] ≥ 1 − δ,
k(n)

i )
where (ξ[n]
i=1 are i.i.d. samples from LT ,Q .

7

An important result of this kind was given by Erdos et al. [ESSW99a].
Theorem 1 (Polynomial Reconstruction [ESSW99a]) Let 0 < f ≤ g < +∞
and denote by Yf,g the set of all phylogenies T = (V, E, [n], ρ; τ ) satisfying f ≤
τe ≤ g, ∀e ∈ E. Then, for all ϕ ≥ 2 and all 0 < f ≤ g < +∞, the phylogenetic
reconstruction problem on Yf,g ⊗ Qϕ can be solved with k = poly(n).
This result was recently improved by Daskalakis et al. [DMR06, DMR11] (see
also [MHR09])
in the so-called Kesten-Stigum reconstruction phase, that is, when
√
g < ln 2.
Definition 6 (∆-Branch Model) Let 0 < ∆ ≤ f ≤ g < +∞ and denote by Yf,g
∆
the set of all phylogenies T = (V, E, [n], ρ; τ ) satisfying f ≤ τe ≤ g where τe is an
integer multiple of ∆, for all e ∈ E. For ϕ ≥ 2 and Q ∈ Qϕ , we call Yf,g
∆ ⊗ {Q}
the ∆-Branch Model (∆-BM).
√
Let g ∗ = ln 2.
Theorem 2 (Logarithmic Reconstruction [DMR06, DMR11, MHR09]) For
CFN }
0 < ∆ ≤ f ≤ g < g ∗ , the phylogenetic reconstruction problem on Yf,g
∆ ⊗{Q
2
can be solved with k = O(log n) .
Distance methods. The proof of Theorem 1 uses distance methods, which we
now define formally.
Definition 7 (Correlation Matrix) Let Φ be a finite set with ϕ ≥ 2. Let
(ξai )ki=1 , (ξbi )ki=1 ∈ Φk
be the sequences at a, b ∈ [n]. For υ1 , υ2 ∈ Φ, we define the correlation matrix
between a and b by
k
1X
ab
b
Fυ1 υ2 =
1{ξai = υ1 , ξbi = υ2 },
k
i=1

and Fbab = (Fbυab1 υ2 )υ1 ,υ2 ∈Φ .
Definition 8 (Distance Method) A phylogenetic reconstruction algorithm A =
i )k
{An,k }n,k≥1 is said to be distance-based if A depends on the data (ξ[n]
i=1 ∈
[n]
k
ab
b
(Φ ) only through the correlation matrices {F }a,b∈[n] .
2

The correct statement of this result appears in [DMR11]. Because of different conventions, our
edge weights are scaled by a factor of 2 compared to those in [DMR11]. The dependence of k in ∆
is ∆−2 .

8

The previous definition takes a very general view of distance-based methods: any
method that uses only pairwise sequence comparisons. In practice, most distancebased approaches actually use a specific distance estimator, that is, a function of
Fbab that converges to τ (a, b) in probability as n → +∞. We give two classical
examples below.
Example 4 (CFN Metric) In the CFN case with state space Φ = {+, −}, a standard distance estimator (up to a constant) is


D(Fb) = − ln 1 − 2(Fb+− + Fb−+ ) .
Example 5 (Log-Det Distance [BH87, Lak94, LSHP94, Ste94]) More generally,
a common distance estimator (up to scaling) is the so-called log-det distance
D(Fb) = − ln | det Fb|.
Loosely speaking, the log-det distance can be thought as a generalization of the
CFN metric. We will use a different generalization of the CFN metric. See section 1.3.

1.2

Results

In our main result, we prove that phylogenies under GTR models of mutation can
be inferred using a distance-based method from k = O(log n) sequence length.
Theorem 3 (Main Result) For all ϕ ≥ 2, 0 < ∆ ≤ f ≤ g < g ∗ and Q ∈ Qϕ ,
there is a distance-based method solving the phylogenetic reconstruction problem
3
on Yf,g
∆ ⊗ {Q} with k = O(log n).
Note that this result is a substantial improvement over Theorem 1—at least, in a
certain range of parameters—and that it matches the bound obtained in Theorem 2.
The result is also novel in two ways over Theorem 2: only the distance matrix is
used; the result applies to a larger class of mutation matrices. A weaker version
of the result stated here was first reported without proof in [Roc08]. Note that
in [Roc08] the result was stated without the discretization assumption which is in
fact needed for the final step of the proof. This is further explained in Section 7.3
of [DMR11]. The new proofs presented here rely on recent joint work with Yuval
Peres [PR11] on exponential moment bounds for quantities such as σ̄a .
In an attempt to keep the paper as self-contained as possible we first give a
proof in the special case of homogeneous trees. This allows to keep the algorithmic details to a minimum. The proof appears in Section 3. We extend the result
to general trees in Section 4. The more general result relies on a combinatorial
algorithm of [DMR11].
3

As in Theorem 2, the dependence of k in ∆ is ∆−2 [Roc10].

9

1.3

Proof Overview

Distance averaging. The basic insight behind Steel’s conjecture is that the accurate reconstruction of ancestral sequences in the reconstruction phase can be harnessed to perform a better reconstruction of the phylogeny itself. For now, consider
the CFN model with character space {+1, −1} and assume that our phylogeny is
homogeneous with uniform branch lengths ω. Generate k i.i.d. samples (σVi )ki=1 .
Let a, b be two internal vertices on level h − h0 < h (from the root). Suppose we
seek to estimate the distance between a and b. This estimation cannot be performed
directly because the sequences at a and b are not known. However, we can try to
estimate these internal sequences. Denote by A, B the leaf set below a and b respectively. An estimate of the sequence at a is the (properly normalized) “site-wise
average” of the sequences at A
1 X σai 0
σ̄ai =
,
(2)
|A| 0 e−ωh0
a ∈A

for i = 1, . . . , k, and similarly for b. It is not immediately clear how such a
site-wise procedure involving simultaneously a large number of leaves can be performed using the more aggregated information in the correlation matrices {Fbuv }u,v∈[n] .
Nevertheless, note that the quantity we are ultimately interested in computing is the
following estimate of the CFN metric between a and b
!
k
1X i i
τ̄ (a, b) = − ln
σ̄a σ̄b .
k
i=1

Our results are based on the following observation:
!
k
1 X 1 X σai 0
τ̄ (a, b) = − ln
k
|A| 0 e−ωh0
i=1

= − ln

a ∈A

X X
1
0
|A||B|e−2ωh 0
0

a ∈A b ∈B

= − ln

1 X σbi 0
|B| 0 e−ωh0
b ∈B
!!
k
1X i i
σa0 σb0
k
i=1
!

X X
1
0 0
e−τ̂ (a ,b )
0
−2ωh
|A||B|e
0
0

!!

,

a ∈A b ∈B

where note that the last line depends only on distance estimates τ̂ (a0 , b0 ) between
leaves a0 , b0 in A, B. In other words, through this procedure, which we call exponential averaging, we perform an implicit ancestral sequence reconstruction using
only distance estimates. One can also think of this as a variance reduction technique. When the branch lengths are not uniform, one needs to use a weighted
version of (2). This requires the estimation of path lengths.
10

GTR models. In the case of GTR models, the standard log-det estimator does not
lend itself well to the exponential averaging procedure described above. Instead,
we use an estimator involving the right eigenvector ν corresponding to the second
eigenvalue Λ2 of Q. For a, b ∈ [n], we consider the estimator


τ̂ (a, b) = − ln ν > Fbab ν .
(3)
This choice is justified by a generalization of (2) introduced in [MP03]. Note that
ν may need to be estimated.
Concentration. There is a further complication in that to obtain results with
high probability, one needs to show that τ̄ (a, b) is highly concentrated. However,
one cannot directly apply standard concentration inequalities because σ̄a is not
bounded. Classical results on the reconstruction problem imply that the variance
of σ̄a is finite—which is not quite enough. To show concentration, we bound the
moment generating function of σ̄a .

1.4

Organization

In Section 2, we provide a detailed account of the connection between ancestral
sequence reconstruction and distance averaging. We then give a proof of our main
result in the case of homogeneous trees in Section 3. In Section 4, we conclude
with a sketch of the proof in the general case.
In the Appendix, we provide a few complimentary results. In Section A, we
show that the distance matrix is not in general a sufficient statistic. In Section B,
we analyze a standard algorithm, known as WPGMA, in the so-called molecular
clock case, that is, when the mutation rate is the same on all branches of the tree. In
particular, in the latter case we note that the discretized branch length assumption
is not needed.

2

Ancestral Reconstruction and Distance Averaging

√
Let ϕ ≥ 2, 0 < ∆ ≤ f ≤ g < g ∗ = ln 2, and Q ∈ Qϕ with corresponding stationary distribution π > 0. In this section we restrict ourselves to the homogeneous
case T = T (h) = (V, E, [n], ρ; τ ) where we take h = log2 n and f ≤ τe ≤ g and
τe is an integer multiple of ∆, ∀e ∈ E. (See Examples 1 and 2 and Theorem 2.)4
4

Note that, without loss of generality, we can consider performing ancestral state reconstruction
on a homogeneous tree as it is always possible to “complete” a general tree with zero-length edges.
We come back to this point in Section 4.

11

Throughout this section, we use a sequence length k > κ log(n) where κ is
a constant to be determined later. We generate k i.i.d. samples (ξVi )ki=1 from the
GTR model (T , Q) with state space Φ.

2.1

Distance Estimator

The standard log-det estimator does not lend itself well to the averaging procedure
discussed above. For reconstruction purposes, we instead use an estimator involving the right eigenvector ν corresponding to the second eigenvalue Λ2 of Q. For
a, b ∈ [n], consider the estimator


(4)
τ̂ (a, b) = − ln ν > Fbab ν ,
where the correlation matrix Fbab was introduced in Definition 7. We first give a
proof that this is indeed a legitimate distance estimator. For more on connections
between eigenvalues of the rate matrix and distance estimation, see e.g. [GL96,
GL98, GMY09].
Lemma 1 (Distance Estimator) Let τ̂ be as above. For all a, b ∈ [n], we have
E[e−τ̂ (a,b) ] = e−τ (a,b) .

Proof: Note that E[Fbijab ] = πi e−τ (a,b)Q ij . Then
h
i

X X 
E ν > Fbab ν =
νi
πi e−τ (a,b)Q νj
i∈Φ

=

X

ij

j∈Φ

νi (πi e−τ (a,b) νi )

i∈Φ

= e−τ (a,b)

X

πi νi2

i∈Φ

= e

−τ (a,b)

.


For a ∈ [n] and i = 1, . . . , k, let
σai = νξai .
Then (4) is equivalent to
τ̂ (a, b) = − ln

k
1X i i
σa σb
k

!
.

(5)

i=1

Note that in the CFN case, we have simply ν = (1, −1)> and hence (5) can be
interpreted as a generalization of the CFN metric.
12

2.2

Ancestral Sequence Reconstruction

Let e = (x, y) ∈ E and assume that x is closest to ρ (in topological distance). We
define Path(ρ, e) = Path(ρ, y), |e|ρ = |Path(v, e)|, and

Rρ (e) = 1 − θe2 Θ−2
ρ,y ,
where Θρ,y = e−τ (ρ,y) and θe = e−τ (e) .
Proposition 1 below is a variant of Lemma 5.3 in [MP03]. For completeness,
we give a proof.
Proposition 1 (Weighted Majority: GTR Version) Let ξ[n] be a sample from LT ,Q
(see Definition 4) with corresponding σ[n] . For a unit flow Ψ from ρ to [n], consider
the estimator
X Ψ(x)σx
S=
.
Θρ,x
x∈[n]

Then, we have
E[S] = 0,
E[S | ξρ ] = σρ ,
and
Var[S] = 1 + KΨ ,
where
KΨ =

X

Rρ (e)Ψ(e)2 .

e∈E

Proof: We follow the proofs of [EKPS00, MP03]. Let ēi be the unit vector in
direction i. Let x ∈ [n], then
> τ (ρ,x)Q
.
E[ē>
ξx | ξρ ] = ēξρ e

Therefore,
τ (ρ,x)Q
E[σx | ξρ ] = ē>
ν = σρ e−τ (ρ,x) ,
ξρ e

and
E[S | ξρ ] =

X Ψ(x)σρ e−τ (ρ,x)
X
= σρ
Ψ(x) = σρ .
Θρ,x

x∈[n]

x∈[n]

In particular,
E[S] =

X
ι∈Φ

13

πi νi = 0.

For x, y ∈ [n], let x ∧ y be the meeting point of the paths between ρ, x, y. We
have
X
E[σx σy ] =
P[ξx∧y = ι]E[σx σy | ξx∧y = ι]
ι∈Φ

=

X

πι E[σx | ξx∧y = ι]E[σy | ξx∧y = ι]

ι∈Φ

=

X

πι e−τ (x∧y,x) νι e−τ (x∧y,y) νι

ι∈Φ

= e−τ (x,y)

X

πι νι2

ι∈Φ
−τ (x,y)

= e

.

Then
Var[S] = E[S 2 ]
X Ψ(x)Ψ(y)
=
E[σx σy ]
Θρ,x Θρ,y
x,y∈[n]
X
=
Ψ(x)Ψ(y)e2τ (ρ,x∧y) .
x,y∈[n]

For e ∈ E, let e = (e↑ , e↓ ) where e↑ is the vertex closest to ρ. Then, by a telescoping sum, for u ∈ V
X
X
X
Rρ (e) =
e2τ (ρ,e↓ ) −
e2τ (ρ,e↑ )
e∈Path(ρ,u)

e∈Path(ρ,u)

e∈Path(ρ,u)

= e2τ (ρ,u) − 1,
and therefore
E[S 2 ] =

X

Ψ(x)Ψ(y)e2τ (v,x∧y)

x,y∈[n]


=

X



Ψ(x)Ψ(y) 1 +

x,y∈[n]

= 1+
= 1+

X

X

Rρ (e)

e∈Path(ρ,x∧y)

Rρ (e)

1{e ∈ Path(ρ, x ∧ y)}Ψ(x)Ψ(y)

X

e∈E

x,y∈[n]

X

Rρ (e)Ψ(e)2 .

e∈E

14


Let Ψ be a unit flow from ρ to [n]. We will use the following multiplicative
decomposition of Ψ: If Ψ(x) > 0, we let
ψ(e) =

Ψ(y)
,
Ψ(x)

and, if instead Ψ(x) = 0, we let ψ(y) = 0. Denoting x↑ the immediate ancestor of
−τ
x ∈ V and letting θx = e (x↑ ,x) , it will be useful to re-write
KΨ =

h−1 X
X

Y

(1 − θx2 )

h0 =0 x∈L(h)

e∈Path(ρ,x)

h0

ψ(e)2
,
θe2

(6)

and to define the following recursion from the leaves. For x ∈ [n],
Kx,Ψ = 0.
Then, let u ∈ V − [n] with children v1 , v2 with corresponding edges e1 , e2 and
define


X
ψ(eα )2
2
Ku,Ψ =
((1 − θvα ) + Kvα ,Ψ )
.
θe2α
α=1,2

Note that, from (6), we have Kρ,Ψ = KΨ .
Because of our use of short sequences, bounds on the variance are not enough
for our purposes: We need exponential concentration on our distance estimates. To
obtain such concentration, we give bounds on the exponential moment of S. Our
proof generalizes a recent argument of Peres and Roch [PR11].
Proposition 2 (Weighted Majority: Exponential Bound) For ζ ∈ R, let
Γi (ζ) = ln E[exp(ζS) | ξρ = νi ].
Then, there exists c > 0 depending only on Q and f such that for all ζ ∈ R, we
have
1
Γi (ζ) ≤ νi ζ + cζ 2 KΨ .
2
Proof: We prove the claim by induction, moving away from the leaves. We begin
with an analytical lemma inspired by the proof of [PR11].

15

Lemma 2 (Recursion Step) Let M = eτ Q with second right eigenvector ν and
corresponding eigenvalue λ = e−τ satisfying τ ≥ f . Then there is c > 0 depending on Q and f such that for all i ∈ Φ
F (x) ≡

X
j∈Φ

1
Mij exp(νj x) ≤ exp(λνi x + c(1 − λ2 )x2 ) ≡ G(x),
2

(7)

for all x ∈ R.
Proof: Let c0 = c(1 − λ2 ). Note that
X
F 0 (x) =
Mij νj exp(νj x),
j∈Φ

F 00 (x) =

X

Mij νj2 exp(νj x),

j∈Φ

1
G0 (x) = (λνi + c0 x) exp(λνi x + c0 x2 ),
2
and

1
G00 (x) = ((λνi + c0 x)2 + c0 ) exp(λνi x + c0 x2 ).
2

Hence,
F (0) = G(0) = 1,
F 0 (0) = G0 (0) = λνi .
Let
π̄ = min πι ,
ι

and

1
ν̄ ≡ max |νi | ≤ √ .
i
π̄

Note that
F 00 (x) ≤ ν̄ 2 exp(ν̄|x|) ≡ F (x),
and

1
G00 (x) ≥ c0 exp(−ν̄|x| + c0 x2 ) ≡ G(x).
2

Choose c0 = c∗ > 0 such that F (x) < G(x) for all x ∈ R. Note in particular that
taking

c∗ > max 4ν̄, ν̄ 2 exp(2ν̄) ,

16

is enough. Indeed, for |x| > 1 we have c∗ > ν̄ 2 and exp(−ν̄|x| + 12 c∗ x2 ) >
exp(ν̄|x|) so that F (x) < G(x). For |x| ≤ 1, we have
G(x) > c∗ exp(−ν̄) > ν̄ 2 exp(ν̄) ≥ F (x).
Now choose c = c∗ (1 − e−2f )−1 in (7) (which implies c0 ≥ c∗ by τ ≥ f ).
Then,
G00 (x) ≥ G(x) > F (x) ≥ F 00 (x),
and therefore
G(x) ≥ F (x),
for all x ∈ R. 
Going back to the proof of Proposition 2, let Sx = σx for all x ∈ [n] and
Su =

X

Svα

α=1,2

ψ(eα )
,
θeα

where u ∈ V − [n] with children v1 , v2 with corresponding edges e1 , e2 . Note that
Sρ = S. Let
Γiu (ζ) = ln E[exp(ζSu ) | ξu = i].
Take c > 0 as in Lemma 2. The main claim is clearly true at the leaves, that is,
for all x ∈ [n]
Γix (ζ) = ln E[exp(ζSx ) | ξx = i]
= ln E[exp(ζσx ) | ξx = i]
= νi ζ
1
≤ νi ζ + cζ 2 Kx,Ψ .
2

17

For u ∈ V −[n] as above, we have by the Markov property, induction, and Lemma 2




X
ψ(e
)
α 
Svα
Γiu (ζ) = ln E exp ζ
| ξu = i
θe α
α=1,2




X
ψ(eα )
ln E exp ζSvα
| ξu = i
=
θeα
α=1,2






X
X
ψ(eα )
ln 
Mijeα E exp ζSvα
| ξvα = j 
=
θe α
α=1,2
j∈Φ





X
X
ψ(eα )

=
ln 
Mijeα exp Γjvα ζ
θe α
α=1,2
j∈Φ





2 !
X
X
ψ(eα )
1
ψ(eα )

+ cKvα ,Ψ ζ
≤
ln 
Mijeα exp νj ζ
θeα
2
θeα
α=1,2
j∈Φ


ψ(eα ) 2
1 2 X
=
cζ
Kvα ,Ψ
2
θe α
α=1,2


 

X
X
ψ(eα )

+
ln 
Mijeα exp νj ζ
θeα
α=1,2
j∈Φ


ψ(eα ) 2
1 2 X
cζ
Kvα ,Ψ
≤
2
θe α
α=1,2




X
1
ψ(eα ) 2
ψ(eα )
2
+
θeα νi ζ
+ c(1 − θvα ) ζ
θeα
2
θeα
α=1,2


X
1
ψ(eα ) 2
= νi ζ + cζ 2
((1 − θv2α ) + Kvα ,Ψ )
2
θeα
α=1,2

1
= νi ζ + cζ 2 Ku,Ψ .
2


2.3

Distance Averaging

The input to our tree reconstruction algorithm is the matrix of all estimated distances between pairs of leaves {τ̂ (a, b)}a,b,∈[n] . For short sequences, these estimated distances are known to be accurate for leaves that are close enough. We now
18

show how to compute distances between internal nodes in a way that involves only
{τ̂ (a, b)}a,b,∈[n] (and previously computed internal weights) using Proposition 2.
(h)

Let 0 ≤ h0 < h. For v ∈ Lh0 , let Tv = (Vv , Ev ) be the subtree of T = T (h)
(h)
rooted at v with leaf set denoted Lv . Let a, b ∈ Lh0 . For x ∈ {a, b}, denote by X
the leaves of T = T (h) below x. Assume that we are given θe , for all e below a, b.
We estimate τ (a, b) as follows
!
X X
1
−1
−1 −τ̂ (a0 ,b0 )
τ̄ (a, b) ≡ − ln
Θa,a0 Θb,b0 e
.
|A||B| 0
0
a ∈A b ∈B

This choice of estimator is suggested by the following observation
X X
0
−1 −τ̂ (a0 ,b0 )
e−τ̄ (a,b) ≡
2−2h Θ−1
a,a0 Θb,b0 e
a0 ∈A b0 ∈B

=

k
1X
k
i=1

X 2−h0 σ i 0
a
0
Θ
a,a
0

!

a ∈A

X 2−h0 σ i 0
b
0
Θ
b,b
0

!
.

b ∈B

Note that the first line depends only on estimates (τ̂ (u, v))u,v∈[n] and {Θv,· }v∈Va ∪Vb .
The last line is the empirical distance between the reconstructed states at a and b
when the flow is chosen to be homogeneous in Proposition 1.
(h)

Lemma 3 (Large Deviations) Let 0 ≤ h0 < h and let a, b ∈ Lh0 . For x = a, b,
let
X 2−h0 σx0
Sx =
.
Θx,x0
0
x ∈X

It holds that
E[e−τ̄ (a,b) ] = e−τ (a,b) ,
and there exists ζ ∗ > 0 small enough such that
E[exp(ζSa Sb )] < +∞,
for all |ζ| < |ζ ∗ |. In particular, for all ε > 0 there exists 0 < χ < 1 such that
h
i
P e−τ̄ (a,b) − E[e−τ̄ (a,b) ] > ε ≤ χk .
Moreover, χ is a constant independent of h0 .

19

Proof: We first prove the expectation formula. Note that
" k
!
!#
X X 2−h0 σ i 0
X 2−h0 σ i 0
1
b
a
E[e−τ̄ (a,b) ] = E
0
0
k
Θ
Θ
a,a
b,b
i=1 a0 ∈A
b0 ∈B
!
!#
"
0
X 2−h0 σb0
X 2−h σa0
= E
Θa,a0
Θb,b0
b0 ∈B
a0 ∈A
" "
!
!
##
X 2−h0 σa0
X 2−h0 σb0
= E E
| ξa , ξb
Θa,a0
Θb,b0
a0 ∈A
b0 ∈B
" "
# "
##
X 2−h0 σa0
X 2−h0 σb0
= E E
| ξa E
| ξb
Θa,a0
Θb,b0
0
0
a ∈A

b ∈B

= E [σa σb ]
= e−τ (a,b) ,
0

where we used that |A| = |B| = 2h .
To prove the large deviation result, it suffices by standard arguments [Dur96]
to bound the exponential moment of
!
!
X 2−h0 σ i 0
X 2−h0 σ i 0
a
b
Sa Sb =
.
Θa,a0
Θb,b0
0
0
a ∈A

b ∈B

2

Let N be Normal(0, 1) and recall that E[eζN ] = eζ /2 . By applying Proposition 2
twice and using Fubini’s Theorem for positive random variables (see also [PR11]),
we get (letting Ψ be the homogeneous flow on T )
1
E[exp(ζSa Sb ) | ξa , ξb ] ≤ E[exp(σa ζSb + cζ 2 Sb2 Ka,Ψ ) | ξa , ξb ]
2
p
= E[exp(σa ζSb + cKa,Ψ ζSb N ) | ξa , ξb ]
p
= E[exp(Sb (σa ζ + cKa,Ψ ζN )) | ξa , ξb ]
p
≤ E[exp(σb (σa ζ + cKa,Ψ ζN )
p
1
+ c(σa ζ + cKa,Ψ ζN )2 Kb,Ψ ) | ξa , ξb ]
2
< +∞,
uniformly in σa , σb for |ζ| > 0 small enough, where we used |σa |, |σb | ≤ ν̄ < +∞,
Cauchy-Schwarz, and
1/2

1
2 2
2
E[ec ζ Ka,Ψ Kb,Ψ N ] =
< +∞,
1 − 2(c2 ζ 2 Ka,Ψ Kb,Ψ )
20

for small enough ζ. Above we used the moment-generating function of the chisquare distribution.5
To prove that the large deviation result is independent of the level h0 , we show
that Ka,Ψ is uniformly bounded in h0 . From (6), we have
Ka,Ψ ≤

0 −1
hX

(1 − e

−2g

h0 −i e

)2

i=0
0

≤

h
X

e2jg e−(2 ln

√

2(h0 −i)g

22(h0 −i)

2)j

j=1
0

=

h
X

e2j(g−g

∗)

j=1

≤

+∞
X

(e−2(g

∗ −g)

)j

j=0

=

1
1−

e−2(g∗ −g)

< +∞,

(8)

√
where recall that g ∗ = ln 2 and g < g ∗ . 
In the next section, we use the previous lemma in two situations: 1) to estimate
the distance between two close vertices on the same level; 2) to detect that two
vertices on the same level are “far apart.” These specializations of Lemma 3 are
stated below. We only sketch the proofs, which are straightforward.
Proposition 3 (Deep Distance Computation: Small Diameter) Let D > 0, γ >
(h)
0, and ε > 0. Let a, b ∈ Lh0 as above. There exist κ > 0 such that if the following
conditions hold:
• [Small Diameter] τ (a, b) < D,
• [Sequence Length] k > κ log(n),
then
|τ̄ (a, b) − τ (a, b)| < ε,
with probability at least 1 − O(n−γ ).
5

A more careful analysis gives the dependence of χ in ∆ as χ = 1 − O(∆2 ) [Roc10].

21

Proof: Let
ε0 = min{(eε − 1)e−D , (1 − e−ε )e−D },
and observe that
τ̄ (a, b) − τ (a, b) < −ε
=⇒ e−τ̄ (a,b) > e−τ (a,b)+ε
=⇒ e−τ̄ (a,b) − e−τ (a,b) > (eε − 1)e−D ≥ ε0 .
A similar implication holds in the other direction. The result now follows from
Lemma 3. 
Proposition 4 (Deep Distance Computation: Diameter Test) Let D > 0, W >
(h)
5, and γ > 0. Let a, b ∈ Lh0 as above. There exists κ > 0 such that if the
following conditions hold:
• [Large Diameter] τ (a, b) > D + ln W ,
• [Sequence Length] k > κ log(n),
then

W
,
2
with probability at least 1 − n−γ . On the other hand, if the first condition above is
replaced by
τ̄ (a, b) > D + ln

• [Small Diameter] τ (a, b) < D + ln W
5 ,
then
τ̄ (a, b) ≤ D + ln

W
,
4

with probability at least 1 − n−γ .
Proof: The proof is similar to the proof of Proposition 3. 

3

Reconstructing Homogeneous Trees

In this section, we prove our main result in the case of homogeneous trees. More
precisely, we prove the following.

22

Theorem 4 (Main Result: Homogeneous Case) Let 0 < ∆ ≤ f ≤ g < +∞
and denote by HYf,g
∆ the set of all homogeneous phylogenies T = (V, E, [n], ρ;
√τ )
satisfying f ≤ τe ≤ g and τe is an integer multiple of ∆, ∀e ∈ E. Let g ∗ = ln 2.
Then, for all ϕ ≥ 2, 0 < ∆ ≤ f ≤ g < g ∗ and Q ∈ Qϕ , there is a distancebased method solving the phylogenetic reconstruction problem on HYf,g
∆ ⊗ {Q}
with k = O(log n).
In the homogeneous case, we can build the tree level by level using simple
“four-point” techniques [Bun71]. See e.g. [SS03, Fel04] for background and details. See also Section 3.2 below. The underlying combinatorial algorithm we use
here is essentially identical to the one used by Mossel in [Mos04]. From Propositions 3 and 4, we get that the “local metric” on each level is accurate as long as we
compute adequate weights. We summarize this fact in the next proposition. For
∆ > 0 and z ∈ R+ , we let [z]∆ be the closest multiple of ∆ to z (breaking ties
arbitrarily). For D > 0, W > 5, we define


W
SD(a, b) = 1 [τ̄ (a, b)]∆ ≤ D + ln
,
3
and we let


d(a, b) =

[τ̄ (a, b)]∆ , if SD(a, b) = 1,
+∞,
o.w.

Proposition 5 (Deep Distorted Metric) Let D > 0, W > 5, and γ > 0. Let
(h)
0
∗
T = (V, E, [n], ρ; τ ) ∈ HYf,g
∆ with g < g . Let a, b ∈ Lh0 for 0 ≤ h < h.
Assume we are given, for x = a, b, θe for all e ∈ Vx . There exists κ > 0, such that
if the following condition holds:
• [Sequence Length] The sequence length is k > κ log(n),
then we have, with probability at least 1 − O(n−γ ),
d(a, b) = τ (a, b)
under either of the following two conditions:
1. [Small Diameter] τ (a, b) < D, or
2. [Finite Estimate] d(a, b) < +∞.
Proof: We let ε < ∆/2. The first part of the proposition follows immediately from
Proposition 3 and the second part of Proposition 4. For the second part, choose κ
so as to satisfy the conditions of Proposition 3 with diameter D + ln W and apply
the first part of Proposition 4. 
It remains to show how to compute the weights, which is the purpose of the next
section.
23

3.1

Estimating Averaging Weights

Proposition 5 relies on the prior computation of the weights θe for all e ∈ Vx , for
x = a, b. In this section, we show how this estimation is performed.
(h)
Let a, b, c ∈ Lh0 . Denote by z the meeting point of the paths joining a, b, c.
We define the “three-point” estimate


1
θ̂z,a = O(a; b, c) ≡ exp − [d(a, b) + d(a, c) − d(b, c)] .
2
Note that the expression in parenthesis is an estimate of the distance between a and
z.
(h)

Proposition 6 (Averaging Weight Estimation) Let a, b, c ∈ Lh0 as above. Assume that the assumptions of Propositions 3, 4, 5 hold. Assume further that the
following condition hold:
• [Small Diameter] τ (a, b), τ (a, c), τ (b, c) < D + ln W ,
then
θ̂z,a = θz,a ,
with probability at least 1 − O(n−γ ) where θ̂z,a = O(a; b, c).
Proof: The proof follows immediately from Proposition 5 and the remark above
the statement of Proposition 6. 

3.2

Putting it All Together
(h)

Let 0 ≤ h0 < h and Q = {a, b, c, d} ⊆ Lh0 . The topology of T (h) restricted
to Q is completely characterized by a bi-partition or quartet split q of the form:
ab|cd, ac|bd or ad|bc. The most basic operation in quartet-based reconstruction
algorithms is the inference of such quartet splits. In distance-based methods in
particular, this is usually done by performing the so-called four-point test: letting
1
F(ab|cd) = [τ (a, c) + τ (b, d) − τ (a, b) − τ (c, d)],
2
we have


 ab|cd if F(a, b|c, d) > 0
ac|bd if F(a, b|c, d) < 0
q=

ad|bc o.w.

Of course, we cannot compute F(a, b|c, d) directly unless h0 = 0. Instead we use
Proposition 5.
24

Deep Four-Point Test. Assume we have previously computed weights θe for all
e ∈ Vx , for x = a, b, c, d. We let
1
F(ab|cd) = [d(a, c) + d(b, d) − d(a, b) − d(c, d)],
2

(9)

and we define the deep four-point test
FP(a, b|c, d) = 1{F(ab|cd) > f /2},
with FP(a, b|c, d) = 0 if any of the distances in (9) is infinite. Also, we extend
the diameter test SD to arbitrary subsets by letting SD(S) = 1 if and only if
SD(x, y) = 1 for all pairs x, y ∈ S.
Algorithm. Fix D > 4g, W > 5, γ > 3. Choose κ so as to satisfy Propositions 5
and 6. Let Z0 be the set of leaves. The algorithm—a standard cherry picking
algorithm—is detailed in Figure 1.
Proof of Theorem 4: The proof of Theorem 4 follows from Propositions 5 and 6.
Indeed, at each level h0 , we are guaranteed by the above to compute a distorted
metric with a radius large enough to detect all cherries on the next level using
four-point tests. The proof follows by induction. 

4

Extension to General Trees

It is possible to generalize the previous arguments to general trees, using a combinatorial algorithm of [DMR11], thereby giving a proof of Theorem 3. To apply
the algorithm of [DMR11] we need to obtain a generalization of Proposition 5 for
disjoint subtrees in “general position.” This is somewhat straightforward and we
give a quick sketch in this section.

4.1

Basic Definitions

The algorithm in [DMR11] is called Blindfolded Cherry Picking. We refer the
reader to [DMR11] for a full description of the algorithm, which is somewhat involved. It is very similar in spirit to the algorithm introduced in Section 3.2, except
for complications due to the non-homogeneity of the tree. The proof in [DMR11]
is modular and relies on two main components: a distance-based combinatorial argument which remains unchanged in our setting; and a statistical argument which
we now adapt. The key to the latter is [DMR11, Proposition 4]. Note that [DMR11,
Proposition 4] is not distance-based as it relies on a complex ancestral reconstruction function—recursive majority. Our main contribution in this section is to show
25

Algorithm
Input: Distance estimates {τ̂ (a, b)}a,b∈[n] ;
Output: Tree;
• For h0 = 1, . . . , h − 1,
1. Four-Point Test. Let
Rh0 = {q = ab|cd : ∀a, b, c, d ∈ Zh0 distinct such that FP(q) = 1}.
2. Cherry Picking. Identify the cherries in Rh0 , that is, those pairs of vertices that only appear on the same side of the quartet splits in Rh0 . Let
(h0 +1)

Zh0 +1 = {a1

(h0 +1)

, . . . , a2h−(h0 +1) },

be the parents of the cherries in Zh0
3. Weight Estimation. For all z ∈ Zh0 +1 ,
(a) Let x, y be the children of z. Choose w to be any other vertex in Zh0
with SD({x, y, w}) = 1.
(b) Compute
θ̂z,x = O(x; y, w).
(c) Repeat the previous step interchanging the role of x and y.

Figure 1: Algorithm.
how this result can be obtained using the techniques of the previous sections—
leading to a fully distance-based reconstruction algorithm.
In order to explain the complications due to the non-homogeneity of the tree
and state our main result, we first need to borrow a few definitions from [DMR11].
Basic Definitions. Fix 0 < ∆ ≤ f ≤ g < g ∗ as in Theorem 3. Let T =
(V, E, [n], ρ; τ ) ∈ Yf,g
∆ be a phylogeny with underlying tree T = (V, E). In this
section, we sometimes refer to the edge set, vertex set and leaf set of a tree T 0 as
E(T 0 ), V(T 0 ), and L(T 0 ) respectively.
Definition 9 (Restricted Subtree) Let V 0 ⊆ V be a subset of the vertices of T .
The subtree of T restricted to V 0 is the tree T 0 obtained by 1) keeping only nodes
and edges on paths between vertices in V 0 and 2) by then contracting all paths
composed of vertices of degree 2, except the nodes in V 0 . We sometimes use the
notation T 0 = T |V 0 . See Figure 2 for an example.

26

Figure 2: Restricting the top tree to its white nodes.

27

u3
u2
u4

u1

u5
u7

u8

u6

Figure 3: The subtrees T |{u1 ,u2 ,u3 ,u8 } and T |{u4 ,u5 ,u6 ,u7 } are edge-disjoint. The
subtrees T |{u1 ,u5 ,u6 ,u8 } and T |{u2 ,u3 ,u4 ,u7 } are edge-sharing.

Definition 10 (Edge Disjointness) Denote by PathT (x, y) the path (sequence of
edges) connecting x to y in T . We say that two restricted subtrees T1 , T2 of T are
edge disjoint if
PathT (x1 , y1 ) ∩ PathT (x2 , y2 ) = ∅,
for all x1 , y1 ∈ L(T1 ) and x2 , y2 ∈ L(T2 ). We say that T1 , T2 are edge sharing if
they are not edge disjoint. See Figure 3 for an example.
Definition 11 (Legal Subforest) We say that a tree is a rooted full binary tree if
all its internal nodes have degree 3 except the root which has degree 2. A restricted
subtree T1 of T is a legal subtree of T if it is also a rooted full binary tree. We say
that a forest
F = {T1 , T2 , . . .},
is legal subforest of T if the Tι ’s are edge-disjoint legal subtrees of T . We denote
by ρ(F) the set of roots of F.
Definition 12 (Dangling Subtrees) We say that two edge-disjoint legal subtrees
T1 , T2 of T are dangling if there is a choice of root for T not in T1 or T2 that is
consistent with the rooting of both T1 and T2 . See Figure 4 below for an example
where two legal, edge-disjoint subtrees are not dangling.
Definition 13 (Basic Disjoint Setup (General)) Let T1 = Tx1 and T2 = Tx2 be
two restricted subtrees of T rooted at x1 and x2 respectively. Assume further that
T1 and T2 are edge-disjoint, but not necessarily dangling. Denote by yι , zι the
28

x2
z2

y2

x1
y1

z1
u2
w2
v2

u1
w1
v1

T2

T1

Figure 4: Basic Disjoint Setup (General). The rooted subtrees T1 , T2 are edgedisjoint but are not assumed to be dangling. The white nodes may not be in the
restricted subtrees T1 , T2 . The case w1 = x1 and/or w2 = x2 is possible. Note that
if we root the tree at any node along the dashed path, the subtrees rooted at y1 and
y2 are edge-disjoint and dangling (unlike T1 and T2 ).

children of xι in Tι , ι = 1, 2. Let wι be the node in T where the path between T1
and T2 meets Tι , ι = 1, 2. Note that wι may not be in Tι since Tι is restricted,
ι = 1, 2. If wι 6= xι , assume without loss of generality that wι is in the subtree
of T rooted at zι , ι = 1, 2. We call this configuration the Basic Disjoint Setup
(General). See Figure 4. Let τ (T1 , T2 ) be the length of the path between w1 and
w2 in the metric τ .

4.2

Deep Distorted Metric

Our reconstruction algorithm for homogeneous trees (see Section 3) builds the tree
level by level and only encounters situations where one has to compute the distance
between two dangling subtrees (that is, the path connecting the subtrees “goes
above them”). However, when reconstructing general trees by growing a subforest
from the leaves, more general situations such as the one depicted in Figure 4 cannot
be avoided and have to be dealt with carefully.
Hence, our goal in this subsection is to compute the distance between the internal nodes x1 and x2 in the Basic Disjoint Setup (General). We have already shown
29

how to perform this computation when T1 and T2 are dangling, as this case is handled easily by Proposition 5 (after a slight modification of the distance estimate; see
below). However, in the general case depicted in Figure 4, there is a complication.
When T1 and T2 are not dangling, the reconstructed sequences at x1 and x2 are not
conditionally independent. But it can be shown that for the algorithm Blindfolded
Cherry Picking to work properly, we need: 1) to compute the distance between
x1 and x2 correctly when the two subtrees are close and dangling; 2) detect when
the two subtrees are far apart (but an accurate distance estimate is not required in
that case). This turns out to be enough because the algorithm Blindfolded Cherry
Picking ensures roughly that close reconstructed subtrees are always dangling. We
refer the reader to [DMR11] for details.
The key point is the following: if one computes the distance between y1 and y2
rather than the distance between x1 and x2 , then the dangling assumption is satisfied (re-root the tree at any node along the path connecting w1 and w2 ). However,
when the algorithm has only reconstructed T1 and T2 , we cannot tell which pair
in {y1 , z1 } × {y2 , z2 } is the right one to use for the distance estimation. Instead,
we compute the distance for all pairs in {y1 , z1 } × {y2 , z2 } and the following then
holds: in the dangling case, all these distances will agree (after subtracting the
length of the edges between x1 , x2 and {y1 , z1 , y2 , z2 }); in the general case, at
least one is correct. This is the basic observation behind the routine D ISTORTED M ETRIC in Figure 5 and the proof of Proposition 7 below. We slightly modify the
definitions of Section 3.
Using the notation of Definition 13, fix (a, b) ∈ {y1 , z1 } × {y2 , z2 }. For x =
a, b, denote by X the leaves of Tx and let |`|x be the graph distance (that is, the
number of edges) between x and leaf ` ∈ X. Assume that we are given θe for all
e ∈ E(Ta ) ∪ E(Tb ). We estimate τ (a, b) as follows
!
X X
0
0
0
0
−1 −τ̂ (a ,b )
τ̄ (a, b) ≡ − ln
2−|a |a −|b |b Θ−1
.
a,a0 Θb,b0 e
a0 ∈A b0 ∈B

Note that, because the tree is binary, it holds that
X X
X
X
0
0
0
0
2−|a |a −|b |b =
2−|a |a
2−|b |b = 1,
a0 ∈A b0 ∈B

a0 ∈A

b0 ∈B

and we can think of the weights on A (similarly for B) as resulting from a homogeneous flow Ψa from a to A. Then, the bounds on the variance and the exponential
moment of
X
0
0
Sa ≡
2−|a |a Θ−1
a,a0 σa ,
a0 ∈A

30

in Propositions 1 and 2 still hold with
X
Ra (e)Ψ(e)2 .
Ka,Ψa =
e∈E(Ta )

Moreover Ka,Ψa is uniformly bounded following an argument identical to (8) in the
proof of Lemma 3. In particular, the same large deviations result hold for τ̄ (a, b).
For D > 0, W > 5, we define


W
SD(a, b) = 1 [τ̄ (a, b)]∆ ≤ D + ln
,
3
and we let


d(a, b) =

[τ̄ (a, b)]∆ , if SD(a, b) = 1,
+∞,
o.w.

Algorithm D ISTORTED M ETRIC
Input: Rooted forest F = {T1 , T2 } rooted at vertices x1 , x2 ; weights τe , for all e ∈
E(T1 ) ∪ E(T2 );
Output: Distance Υ;
• [Children] Let yι , zι be the children of xι in F for ι = 1, 2 (if xι is a leaf, set
zι = yι = xι );
• [Distance Computations] For all pairs (a, b) ∈ {y1 , z1 } × {y2 , z2 }, compute
D(a, b) := d(a, b) − τ (a, x1 ) − τ (b, x2 );
• [Multiple Test] If
(1)

(1)

(2)

(2)

max{ D(r1 , r2 ) − D(r1 , r2 ) :
(ι)

(ι)

(r1 , r2 ) ∈ {y1 , z1 } × {y2 , z2 }, ι = 1, 2} = 0,
return Υ := D(z1 , z2 ), otherwise return Υ := +∞ (return Υ := +∞ if any of
the distances above is +∞).

Figure 5: Routine D ISTORTED M ETRIC.
Proposition 7 (Accuracy of D ISTORTED M ETRIC) Let D > 0, W > 5, γ > 0
and g < g 0 < g ∗ . Consider the Basic Disjoint Setup (General) with F = {T1 , T2 }
and Q = {y1 , z1 , y2 , z2 }. Assume we are given θe for all e ∈ E(T1 ) ∪ E(T2 ). Let
Υ denote the output of D ISTORTED M ETRIC in Figure 5. There exists κ > 0, such
that if the following condition holds:
31

• [Edge Length] It holds that τ (e) ≤ g 0 , ∀e ∈ E(Tx ), x ∈ Q6 ;
• [Sequence Length] The sequence length is k > κ log(n),
then we have, with probability at least 1 − O(n−γ ),
Υ = τ (x1 , x2 )
under either of the following two conditions:
1. [Dangling Case] T1 and T2 are dangling and τ (T1 , T2 ) < D, or
2. [Finite Estimate] Υ < +∞.
Proof: The proof, which is a simple combination of the proof of Proposition 5 and
the remarks above the statement of Proposition 7, is left out. 
Full Algorithm. The rest of the Blindfolded Cherry Picking algorithm is unchanged except for an additional step to compute averaging weights as in the algorithm of Section 3. This concludes our sketch of the proof of Theorem 3.

6

For technical reasons explained in [DMR11], we allow edges slightly longer than the upper
bound g.

32

Acknowledgments
This work was triggered by a discussion with Elchanan Mossel about lower bounds
for distance methods, following a talk of Joseph Felsenstein. Elchanan pointed out
that the distance matrix has a potentially useful correlation structure. I am also
indebted to Yuval Peres [PR11].

References
[Att99]

K. Atteson. The performance of neighbor-joining methods of phylogenetic reconstruction. Algorithmica, 25(2-3):251–278, 1999.

[BH87]

Daniel Barry and J. A. Hartigan. Statistical analysis of hominoid
molecular evolution. Statist. Sci., 2(2):191–210, 1987. With comments by Stephen Portnoy and Joseph Felsenstein and a reply by the
authors.

[Bun71]

P. Buneman. The recovery of trees from measures of dissimilarity.
In Mathematics in the Archaelogical and Historical Sciences, pages
187–395. Edinburgh University Press, Edinburgh, 1971.

[CH91]

Joseph T. Chang and John A. Hartigan. Reconstruction of evolutionary trees from pairwise distributions on current species, 1991.

[Cha96]

Joseph T. Chang. Full reconstruction of Markov models on evolutionary trees: identifiability and consistency. Math. Biosci., 137(1):51–
73, 1996.

[CK01]

Miklós Csurös and Ming-Yang Kao. Provably fast and accurate recovery of evolutionary trees through harmonic greedy triplets. SIAM
Journal on Computing, 31(1):306–322, 2001.

[Csu02]

M. Csurös. Fast recovery of evolutionary trees with thousands of
nodes. J. Comput. Biol., 9(2):277–97, 2002.

[CT06]

Benny Chor and Tamir Tuller. Finding a maximum likelihood tree is
hard. J. ACM, 53(5):722–744, 2006.

[Day87]

William H. E. Day. Computational complexity of inferring phylogenies from dissimilarity matrices. Bull. Math. Biol., 49(4):461–467,
1987.

33

[DHJ+ 06]

Constantinos Daskalakis, Cameron Hill, Alexander Jaffe, Radu Mihaescu, Elchanan Mossel, and Satish Rao. Maximal accurate forests
from distance matrices. In RECOMB, pages 281–295, 2006.

[DMR06]

Constantinos Daskalakis, Elchanan Mossel, and Sébastien Roch. Optimal phylogenetic reconstruction. In STOC’06: Proceedings of the
38th Annual ACM Symposium on Theory of Computing, pages 159–
168, New York, 2006. ACM.

[DMR09]

Constantinos Daskalakis, Elchanan Mossel, and Sébastien Roch.
Phylogenies without branch bounds: Contracting the short, pruning
the deep. In RECOMB, pages 451–465, 2009.

[DMR11]

Constantinos Daskalakis, Elchanan Mossel, and Sébastien Roch.
Evolutionary trees and the ising model on the bethe lattice: a proof of
steel’s conjecture. Probability Theory and Related Fields, 149:149–
189, 2011. 10.1007/s00440-009-0246-2.

[DS86]

William H. E. Day and David Sankoff. Computational complexity of
inferring phylogenies by compatibility. Syst. Zool., 35(2):224–229,
1986.

[Dur96]

Richard Durrett. Probability: theory and examples. Duxbury Press,
Belmont, CA, second edition, 1996.

[EKPS00]

W. S. Evans, C. Kenyon, Y. Peres, and L. J. Schulman. Broadcasting
on trees and the Ising model. Ann. Appl. Probab., 10(2):410–433,
2000.

[ESSW99a] P. L. Erdös, M. A. Steel, L. A. Székely, and T. A. Warnow. A few
logs suffice to build (almost) all trees (part 1). Random Struct. Algor.,
14(2):153–184, 1999.
[ESSW99b] P. L. Erdös, M. A. Steel, L. A. Székely, and T. A. Warnow. A few
logs suffice to build (almost) all trees (part 2). Theor. Comput. Sci.,
221:77–118, 1999.
[Fel04]

J. Felsenstein. Inferring Phylogenies. Sinauer, Sunderland, MA,
2004.

[Gas97]

O. Gascuel. BIO-NJ: An improved version of the NJ algorithm based
on a simple model of sequence data. Mol. Biol. Evol., 14(7):685–695,
1997.
34

[Geo88]

H. O. Georgii. Gibbs measures and phase transitions, volume 9 of
de Gruyter Studies in Mathematics. Walter de Gruyter & Co., Berlin,
1988.

[GF82]

R. L. Graham. and L. R. Foulds. Unlikelihood that minimal phylogenies for a realistic biological study can be constructed in reasonable
computational time. Math. Biosci., 60:133–142, 1982.

[GL96]

X Gu and W H Li. A general additive distance with time-reversibility
and rate variation among nucleotide sites. Proceedings of the
National Academy of Sciences of the United States of America,
93(10):4671–4676, 1996.

[GL98]

Xun Gu and Wen-Hsiung Li. Estimation of evolutionary distances
under stationary and nonstationary models of nucleotide substitution.
Proceedings of the National Academy of Sciences of the United States
of America, 95(11):5899–5905, 1998.

[GMS08]

Ilan Gronau, Shlomo Moran, and Sagi Snir. Fast and reliable reconstruction of phylogenetic trees with very short edges. In SODA
’08: Proceedings of the nineteenth annual ACM-SIAM symposium on
Discrete algorithms, pages 379–388, Philadelphia, PA, USA, 2008.
Society for Industrial and Applied Mathematics.

[GMY09]

Ilan Gronau, Shlomo Moran, and Irad Yavneh. Towards optimal distance functions for stochastic substitutions models. Preprint, 2009.

[HNW99]

D. H. Huson, S. H. Nettles, and T. J. Warnow. Disk-covering, a fastconverging method for phylogenetic tree reconstruction. J. Comput.
Biol., 6(3–4), 1999.

[KZZ03]

Valerie King, Li Zhang, and Yunhong Zhou. On the complexity of
distance-based evolutionary tree reconstruction. In SODA ’03: Proceedings of the fourteenth annual ACM-SIAM symposium on Discrete
algorithms, pages 444–453, Philadelphia, PA, USA, 2003. Society for
Industrial and Applied Mathematics.

[Lak94]

JA Lake. Reconstructing Evolutionary Trees from DNA and Protein Sequences: Paralinear Distances. Proceedings of the National
Academy of Sciences, 91(4):1455–1459, 1994.

[LC06]

Michelle R. Lacey and Joseph T. Chang. A signal-to-noise analysis
of phylogeny estimation by neighbor-joining: insufficiency of polynomial length sequences. Math. Biosci., 199(2):188–215, 2006.
35

[Lig85]

Thomas M. Liggett. Interacting particle systems, volume 276
of Grundlehren der Mathematischen Wissenschaften [Fundamental
Principles of Mathematical Sciences]. Springer-Verlag, New York,
1985.

[LSHP94]

PJ Lockhart, MA Steel, MD Hendy, and D Penny. Recovering Evolutionary Trees under a More Realistic Model of Sequence. Mol Biol
Evol, 11(4):605–612, 1994.

[MHR09]

R. Mihaescu, C. Hill, and S. Rao. Fast phylogeny reconstruction
through learning of ancestral sequences. Preprint, 2009.

[Mos03]

E. Mossel. On the impossibility of reconstructing ancestral data and
phylogenies. J. Comput. Biol., 10(5):669–678, 2003.

[Mos04]

E. Mossel. Phase transitions in phylogeny. Trans. Amer. Math. Soc.,
356(6):2379–2404, 2004.

[Mos07]

E. Mossel. Distorted metrics on trees and phylogenetic forests.
IEEE/ACM Trans. Comput. Bio. Bioinform., 4(1):108–116, 2007.

[MP03]

E. Mossel and Y. Peres. Information flow on trees. Ann. Appl.
Probab., 13(3):817–844, 2003.

[MV05]

Elchanan Mossel and Eric Vigoda. Phylogenetic MCMC Algorithms
Are Misleading on Mixtures of Trees. Science, 309(5744):2207–
2209, 2005.

[PR11]

Yuval Peres and Sébastien Roch. Reconstruction on trees: Exponential moment bounds for linear estimators. Electron. Comm. Probab.,
16:251–261 (electronic), 2011.

[Roc06]

Sébastien Roch. A short proof that phylogenetic tree reconstruction
by maximum likelihood is hard. IEEE/ACM Trans. Comput. Biology
Bioinform., 3(1):92–94, 2006.

[Roc08]

Sébastien Roch. Sequence-length requirement for distance-based
phylogeny reconstruction: Breaking the polynomial barrier. In
FOCS, pages 729–738, 2008.

[Roc10]

Sebastien Roch. Toward extracting all phylogenetic information from
matrices of evolutionary distances. Science, 327(5971):1376–1379,
2010.

36

[RS96]

A Rzhetsky and T Sitnikova. When is it safe to use an oversimplified
substitution model in tree- making? Mol Biol Evol, 13(9):1255–1265,
1996.

[SHP88]

M. A. Steel, M. D. Hendy, and D. Penny. Loss of information in
genetic distances. Nature, 336(6195):118, 1988.

[SN87]

N. Saitou and M. Nei. The neighbor-joining method: A new method
for reconstructing phylogenetic trees. Mol. Biol. Evol., 4(4):406–425,
1987.

[SS63]

R. Sokal and P. Sneath. Principles of Numerical Taxonomy. W. H.
Freeman and Co., San Francisco, Calif., 1963.

[SS99]

Michael A. Steel and László A. Székely. Inverting random functions.
Ann. Comb., 3(1):103–113, 1999. Combinatorics and biology (Los
Alamos, NM, 1998).

[SS02]

M. A. Steel and L. A. Székely. Inverting random functions. II. Explicit bounds for discrete maximum likelihood estimation, with applications. SIAM J. Discrete Math., 15(4):562–575 (electronic), 2002.

[SS03]

C. Semple and M. Steel. Phylogenetics, volume 22 of Mathematics
and its Applications series. Oxford University Press, 2003.

[Ste94]

M. Steel. Recovering a tree from the leaf colourations it generates
under a Markov model. Appl. Math. Lett., 7(2):19–23, 1994.

[Ste01]

M. Steel. My Favourite Conjecture. Preprint, 2001.

[Was04]

Larry Wasserman. All of statistics. Springer Texts in Statistics.
Springer-Verlag, New York, 2004. A concise course in statistical inference.

37

A

The Distance Matrix is Not Sufficient

A statistic (i.e., a function of the full data) is called sufficient if, conditioned on the
value of the statistic, the distribution of the full data does not depend on the parameters of the generating model. Roughly speaking, a sufficient statistic encapsulates
all the information about the data. See e.g. [Was04]. In this section, we show that
the pairwise correlation matrices do not constitute a sufficient statistic for the full
Markov model of evolution. Hence, there is in principle more information in the
full sequence dataset than there is in the matrix of evolutionary distances.
We give a simple example of non-sufficiency. Consider a four-leaf tree with
leaf set L = {a, b, c, d} and split ab|cd. Assume we use a CFN model with purines
denoted “0” and pyrimidines denoted “1” with equal mutation probabilities p. Consider the following correlation matrices
1
Fbυij1 υ2 = ,
4
for all i 6= j ∈ L and υ1 , υ2 ∈ {0, 1}.
correlation matrices are

0 0
 0 1
Data1 = 
 0 1
1 1
and



0
 1
Data2 = 
 1
0

0
1
0
1

Two different datasets consistent with these
0
1
0
0

0
0
1
0

1
1
1
0

1
0
0
0

1
0
1
1


1
1 
,
0 
1

0
0
1
1

0
0
0
0

1
0
0
1

1
0
1
0

1
1
0
0


1
1 
,
1 
1

where the columns are the sites and the rows are the leaves in the order a, b, c, d.
We compare the probability of observing the two datasets under two different
values of p: p = ε and = 1/2 − ε for ε > 0 small. In the first case, in a first
approximation it suffices to compute the parsimony scores and we have
Pε [Data1 ] =
and

 ε 8
2

+ O(ε9 ) =

ε8
+ O(ε9 ),
256

 2  
1
ε 2 2 4
ε10
Pε [Data2 ] =
ε
+ O(ε11 ) =
+ O(ε11 ).
2
2
16

38

In particular, we get the ratio
Pε [Data2 | Fb]
Pε [Data2 ]
=
= ε2 + O(ε3 ).
b
P
[Data
]
ε
1
Pε [Data1 | F ]
On the other hand, if p = 1/2 − ε then the state distribution is almost uniform
and we get
P1/2−ε [Data2 | Fb]
P1/2−ε [Data2 ]
=
= 1 + O(ε).
P1/2−ε [Data1 ]
P1/2−ε [Data1 | Fb]
Since the ratios are different, we have shown that the distribution of the data conditioned on the correlation matrices depends on the parameters of the model. Therefore, the distance matrix is not a sufficient statistic.

B

Probabilistic Analysis of WPGMA

Let 0 < f < g < +∞ and denote by UYf,g the set of all phylogenies T =
(V, E, [n], ρ; τ ) ∈ Yf,g where we have further that τ is ultrametric, that is, for all
v ∈ V it holds that τ (v, x) = τ (v, y) ≡ τ (v), for all leaves x, y below v. This
is known as the molecular clock assumption, that is, the case where the mutation
rate is equal on all edges. In that case, there are particularly simple clustering
algorithms. We recall the WPGMA algorithm in Figure 6. In the molecular clock
case, it is enough to consider “uncorrected” distances [RS96]. Therefore, we run
WPGMA with the uncorrected distance estimates
1 − ω̂(a, b)
,
2

τ̂u (a, b) =
where

ω̂(a, b) = ν > Fbab ν,
for a, b ∈ [n]. We call a subset of leaves A a clade if it corresponds to all
leaf descendants of an internal node a∗ called the most recent common ancestor
(MRCA). For a clade A with MRCA a∗ and a leaf a ∈ A, we let |a|A = |a|a∗ and
ΘA = Θa∗ ,a . For disjoint clades A and B, we let
τ̂u (A, B) =

XX

2−|a|A 2−|b|B τ̂u (a, b) =

a∈A b∈B

1 − ω̂(A, B)
,
2

where
ω̂(A, B) =

XX

2−|a|A 2−|b|B ω̂(a, b).

a∈A b∈B

39

We define
ω(a, b) = e−τ (a,b) ,
and
1 − e−τ (a,b)
.
2
And similarly for ω(A, B) and τu (A, B).
Throughout this section, we use a sequence length k > κ log(n) where κ is a
constant to be determined later.
τu (a, b) =

Algorithm WPGMA
Input: Distance estimates {τ̂u (a, b)}a,b∈[n] ;
Output: Tree;
• Initialization. Let Z0 be the set of leaves as clusters, that is,
Z0 = {{l} : l ∈ [n]},
and for all a, b ∈ [n] let
τ̂u ({a}, {b}) = τ̂u (a, b).
• Main Loop. For i = 1, . . . , n − 1,
– Selection Step. Let
(A∗ , B ∗ ) ∈ arg min{τ̂u (A, B) : A, B ∈ Zi−1 distinct}.
Merge clusters A∗ , B ∗ to obtain Zi .
– Reduction Step. For all C ∈ Zi − {A∗ ∪ B ∗ }, compute
τ̂u (C, A∗ ∪ B ∗ ) =

1
[τ̂u (C, A∗ ) + τ̂u (C, B ∗ )].
2

(10)

• Output. Output tree implied by the successive clusterings Z0 , . . . , Zn−1 .

Figure 6: Algorithm WPGMA.
Theorem 5 (Analysis of WPGMA) For all 0 < f < g < g ∗ , WPGMA solves the
phylogenetic reconstruction problem on UYf,g ⊗ {Q} with k = O(log n).
Proof: Fix D > 3g + 2f , 2g + 2f < D < D, and
 2f

e − 1 eD−2g−2f − 1
0
ε < min 2f
,
.
e + 1 eD−2g−2f + 1
40

This choice ensures that
e2f

1 − ε0
> 1,
1 + ε0

and
eD−2g−2f

1 − ε0
> 1,
1 + ε0

which will be needed later. Let
ε = min{ε0 e−D , ε0 e−D },
and let χ be as in Lemma 3 for this choice of ε. Taking κ large enough, assume the
conclusion of Lemma 3 holds for all pairs of clades in the tree, an event we denote
by (?).
By definition, we have
τ̂u (A, B) ≤ τ̂u (A0 , B 0 )

ω̂(A, B) ≥ ω̂(A0 , B 0 ).

⇐⇒

For convenience, in the rest of the proof we work with ω̂ rather than τ̂u . If A, B
are disjoint clades with respective MRCA a∗ and b∗ satisfying τ (a∗ , b∗ ) < D, we
have
ω̂(A, B) < ω(A, B) + ΘA ΘB ε
∗ ,b∗ )

+ ε0 e−D )

−τ (a∗ ,b∗ )

+ ε0 e−τ (a

≤ ΘA ΘB (e−τ (a
< ΘA ΘB (e

∗ ,b∗ )

)

0

= ω(A, B)(1 + ε ),
and similarly
ω̂(A, B) > ω(A, B)(1 − ε0 ).
On the other hand, if τ (a∗ , b∗ ) > D, we have
ω̂(A, B) < ω(A, B) + ΘA ΘB ε
∗ ,b∗ )

≤ ΘA ΘB (e−τ (a

+ ε0 e−D )

< ΘA ΘB (e−D + ε0 e−D )
= ΘA ΘB e−D (1 + ε0 ).
By (?) these inequalities hold for all such pairs of clades.
Two clades A, B are sister clades if their MRCA is their immediate ancestor.
We use the following convention. Recall that the leaves are denoted {1, . . . , n}. We
let min A be the smallest label in A. When denoting a pair of sister clades (A, B),
41

we always assume min A < min B. There are n−1 pairs of sister clades. Order the
sister pairs by decreasing value of ω̂(A, B), breaking ties by lexicographic order
over (min A, min B):
(A1 , B1 ), . . . , (An−1 , Bn−1 ).
We assume that WPGMA uses the same tie-breaking rule. We let Ci = Ai ∪ Bi .
We prove the following basic claim. For all i = 1, . . . , n − 1, at Selection
Step i we choose (A∗ , B ∗ ) = (Ai , Bi ). The result then follows. We work by
induction. For i = 0, there is nothing to prove. Assume the claim holds up to some
1 ≤ i < n − 1. We make a series of observations:
1. All the current clusters in Zi−1 are clades. This follows from the induction hypothesis. By the induction hypothesis, we also get that the values
τ̂u (A, B) computed at the Reduction Steps indeed correspond to our original definition:
τ̂u (A, B) =

XX

2−|a|A 2−|b|B τ̂u (a, b) =

a∈A b∈B

1 − ω̂(A, B)
,
2

where
ω̂(A, B) =

XX

2−|a|A 2−|b|B ω̂(a, b).

a∈A b∈B

2. We show that for all C ∈ Zi−1 , we have
ω(Ai , Bi )e−2f < Θ2C ≤ ω(Ai , Bi )e2g+2f .
Let C ∈ Zi−1 such that C = A ∪ B for sister clades A, B. By (?), we have
that
Θ2C

= ω(A, B)
> ω̂(A, B)(1 + ε0 )−1
> ω̂(Ai , Bi )(1 + ε0 )−1
1 − ε0
> ω(Ai , Bi )
1 + ε0
> ω(Ai , Bi )e−2f .

Conversely, if a clade C = A ∪ B with sister clades A, B satisfies
Θ2C = ω(A, B) > ω(Ai , Bi )e2f ,

42

(11)

then
ω̂(A, B) > (1 − ε0 )ω(A, B)
> (1 − ε0 )ω(Ai , Bi )e2f
1 + ε0
> (1 − ε0 )ω(Ai , Bi )
1 − ε0
0
> (1 + ε )ω(Ai , Bi )
> ω̂(Ai , Bi ),

(12)

so that C must be included in a cluster of Zi by our induction hypothesis. In
particular, if two sister clades A, B are such that Θ2A , Θ2B > ω(Ai , Bi )e2g+2f
then (11) is satisfied, that is, ω(A, B) > ω(Ai , Bi )e2f . By (12), (A, B)
would have been selected in a previous iteration by induction. That implies,
for all C ∈ Zi−1 ,
Θ2C ≤ ω(Ai , Bi )e2g+2f .
3. We claim that Ai , Bi ∈ Zi−1 . Indeed, by the previous paragraph all clades
with Θ2 -value at least ω(Ai , Bi )e2f have been constructed in a previous
iteration. In particular, the clade Ai has been constructed in a previous step
as it satisfies
p
ΘAi e−f > ΘCi = ω(Ai , Bi ).
The same holds for Bi . Moreover, Ai and Bi being sister clades of each
other (and no other clades), they cannot have been selected inside another
pair by our induction hypothesis.
4. By construction, (Ai , Bi ) is chosen over all other sister clades present in
Zi−1 . So it remains to show that (Ai , Bi ) is selected over all other pairs.
Pairs of clades that are far enough will not be selected. That is, if A, B with
MRCA a∗ , b∗ is such that
τ (a∗ , b∗ ) ≥ D,
then
ω̂(A, B) < ΘA ΘB e−D (1 + ε0 )
< ω(Ai , Bi )e2g+2f e−D (1 + ε0 )
< ω̂(Ai , Bi )(1 − ε0 )−1 e2g+2f e−D (1 + ε0 )
< ω̂(Ai , Bi ),
by assumption on ε0 .
43

5. Finally, non-sister clades that are closer than D cannot be selected. Indeed,
assume by contradiction that (A∗ , B ∗ ) is such a pair. Since (A∗ , B ∗ ) are not
sister clades, at least one of them, say A∗ without loss of generality, has an
immediate ancestor u that is stricly lower than the MRCA of A∗ and B ∗ .
Take C ∗ to be any clade in Zi−1 below u that is different than A∗ . There
must be such a clade because otherwise A∗ would have been merged with its
sister already. The MRCA of A∗ and C ∗ is u. Moreover, we must have
Θ2A∗ > ω(Ai , Bi )e−2f ,
and
Θ2C ∗ ≤ ω(Ai , Bi )e2g+2f ,
so that
τ (a∗ , c∗ ) < 2g + g + 2f < 3g + 2f < D,
where a∗ and c∗ are the MRCA of A∗ and C ∗ respectively. Finally by (?)
ω̂(A∗ , C ∗ ) > ω(A∗ , C ∗ )(1 − ε0 )
> ω(A∗ , B ∗ )e2f (1 − ε0 )
> ω̂(A∗ , B ∗ )(1 + ε0 )−1 e2f (1 − ε0 )
> ω̂(A∗ , B ∗ ).
This is a contradiction.


44

