\section{Introduction}

High-throughput sequencing is making large collections of sequences available to researchers at a low cost.
These genomic data represent a broad spectrum of life and motivates studies of the problem of reconstructing large phylogenetic trees using statistical methods.
Those data, however, also come from different sources, cover different genomic regions on which the evolutionary processes happen very differently, and may not be collected on the same set of species. 
Thus, combining trees on different, overlapping sets of species into a ``supertree" has become a popular approach for reconstructing large species trees.
Over the years, several methods of supertree reconstruction have been developed \citep{cotton2007majority}, and analyses of supertrees continue to play more and more important roles in the search for answers of many fundamental evolutionary questions.

While supertree methods have been of great interest in phylogenetics, little is known about their theoretical properties, especially in non-asymptotic settings. 
In cases when the individual trees are gene trees and the set of taxa are the same across the input trees, several statistical consistent methods have been derived \citep{ mossel2008incomplete, heled2009bayesian, kubatko2009stem, larget2010bucky, liu2010maximum, liu2011estimating, bryant2012inferring, mirarab2014astral, chifman2014quartet, vachaspati2015astrid}. 
However, most proofs of statistical consistency have been analyzed under the condition that each of the gene trees can be estimated accurately, and do not provide a guarantee of robustness to gene tree estimation errors  \citep{roch2015robustness}.
This is a critical concern because when the species tree is constructed using sequences taken from large genomic regions, those regions have a high chance of involving some recombination, which violates one of the main assumptions of the multi-species coalescent model.
On the other hand, limiting analyses to short regions increase gene tree estimation error (GTEE) and various summary methods had impaired accuracy when the error was high \citep{gatesy2014phylogenetic, mirarab2014statistical}.
It is later proved in \cite{roch2019long} that when the sequence length of each locus is bounded and the gene tree cannot be estimated reliably, most summary methods that estimate the species tree by combining gene trees are not statistically consistent.

Because some coalescent-based summary methods sometimes produce less accurate estimates than concatenation \citep{bayzid2013naive, patel2013error, mirarab2014statistical}, seemingly as a result of GTEE, the question of whether provable guarantees can be established in the presence of GTEE (for both species trees and supertree) naturally arises. 
In the context of species tree estimation, \cite{roch2015robustness} establish statistical consistency of the Rooted Triplet Consensus method and the Maximum Pseudo-likelihood for Estimating Species Trees method \citep{liu2010maximum} under GTEE and provide bounds on the sampling complexity for these methods to construct the correct species tree with high probability. 

The maximum likelihood (ML) supertree method is proposed by \cite{steel2008maximum} based on a probability model that permits ``errors" in gene tree topologies and allows the species tree to be estimated even if there is topological conflict amongst gene trees.
\cite{steel2008maximum} shows that ML estimate of the species tree is topologically consistent under fairly general conditions and also shows that the method of Matrix Representation with Parsimony \citep{baum1992combining} may be inconsistent under these same conditions.
However, while the ML estimate is topologically consistent, no results on the convergence rate of the estimator are obtained. 
In this paper, we propose a new analytic approach to study the convergence of the ML supertree method. 
Based on embedding trees into a metric space, we establish the conditions for which the convergence rate of the ML supertree can be obtained. 
We verify these conditions for the popular exponential error model of gene trees, thus obtain a polynomial convergence rate for the ML supertree to the true species tree under this model.

\section{Mathematical framework}

In this paper, the term \emph{phylogenetic tree} refers to a tree $\mathbb{T}$ with leaves labeled by a set of species.
Each branch of $\mathbb{T}$ is associated with a non-negative branch length.
A tree is said to be resolved if it is bifurcating and all branch lengths are positive.
Given an unrooted phylogenetic tree $\mathbb{T}$ on a finite set $\mathscr{L}$ of species, any subset $\mathscr{L}'$ of $\mathscr{L}$ induces a phylogenetic tree on $\mathscr{L}'$, denoted $\mathbb{T} |_{\mathscr{L}'}$, which is the subtree of $\mathbb{T}$ that connects the species in $\mathscr{L}'$ only.

We regard the tree space as a metric space $(\mathcal{T}, d)$ where $\mathcal{T}$ is the set of all phylogenetic trees with branch lengths bounded from above by a positive constant $g_0$ and $d$ is a continuous metric. 
For simplicity of presentation, we will assume that $d$ is the BHV distance on the set of trees with $n$ species \citep{billera2001geometry}, but note that the analysis of the paper can be extended to any continuous and locally-Euclidean distance, including the branch-score distance \citep{kuhner1994simulation} and the AGPS distance \citep{amenta2007approximating}.

To describe trees that are ``near" to each other, the BHV distance use the class of nearest neighbor interchange (NNI) moves \citep{robinson1971comparison}.
An NNI move is defined as a transformation that collapses an interior branch to zero and then expands the resulting degree 4 vertex into a branch in a different way.
The BHV space models the set of trees $\mathcal{T}$ on $n$ species as a cubical complex consisting of a collection of orthants, each isomorphic to $\mathbb{R}_{\ge 0}^{2n-3}$. 
Each orthant of $\mathcal{T}$ corresponds uniquely to a tree topology, and the coordinates in each orthant parameterize the branch lengths for the corresponding tree.
The adjacent orthants of the complex with the same dimension correspond to NNI-adjacent trees.

The BHV space is equipped with a natural metric distance: the shortest path lying in the BHV space between the points. If two points lie in the same orthant, this distance is the usual Euclidean distance.
If two points are in different orthants, they can be joined by a sequence of straight segments, with each segment lying in a single orthant.
We can then measure the length of the path by adding up the lengths of the segments.
The distance between the two trees $\mathbb{T}_1$ and $\mathbb{T}_2$ on the BHV space is defined as the minimum of the lengths of such segmented paths joining the two points.

Throughout the paper, we assume that the evolution of the species has not involved reticulate processes, and there exists an underlying ``true species tree" in $\mathcal{T}$, denoted by $\mathbb{T}^*$.
Furthermore, $\mathbb{T}^*$ is a resolved tree with leaf set $\mathscr{L}^*$.
In a supertree reconstruction problem, we observe a sequence of gene trees $(\mathbb{T}_1, \mathbb{T}_2, \ldots, \mathbb{T}_k)$, where $\mathbb{T}_i$ has leaf set $\mathscr{L}_i$, and wish to combine these trees into a phylogenetic tree $\hat \mathbb{T}_k$ on the union of the leaf sets 
\[
\hat \mathscr{L}_k = \cup_{i=1}^k \mathscr{L}_i \subset \mathscr{L}^*.
\] 

In this paper, we consider a tree-generating probability model $\left( \mathcal{P}_{\mathbb{T}} \right )_{\mathbb{T} \in \mathcal{T}}$ such that for each set of species $\mathscr{L} \subset \mathscr{L}^*$, $\mathcal{P}^{\mathscr{L}}_{\mathbb{T}}(\cdot)$ is a distribution of trees forming by species in $\mathscr{L}$.
We assume that the observed gene trees $\{\mathbb{T}_i\}_{i=1}^k$ are independently distributed according to $\left \{ \mathcal{P}^{\mathscr{L}_i}_{\mathbb{T}^*}(\cdot) \right \}_{i=1}^k$ respectively.
In other words, the joint probability of the observed gene trees is
\[
\mathcal{P}_{\mathbb{T}^*}(\mathbb{T}_1, \mathbb{T}_2, \ldots, \mathbb{T}_k) = \prod_{i=1}^k{ \mathcal{P}^{\mathscr{L}_i}_{\mathbb{T}^*}(\mathbb{T}_i) }.
\]
\begin{Remark}
We note that the tree-generating probability model $\mathcal{P}_{\mathbb{T}^*}$ may not necessarily be nested. 
That is, given two nested sets of leaves $\mathscr{L}_1 \subset \mathscr{L}_2 \subset \mathscr{L}^*$, although the probability $\mathcal{P}^{\mathscr{L}_2}_{\mathbb{T}^*}$ (which is used to generate trees with leaf set $\mathscr{L}_2$) also induces a natural distribution on the set of trees with leaf set $\mathscr{L}_1$, this probability $\mathcal{P}^{\mathscr{L}_2}_{\mathbb{T}^*} |_{\mathscr{L}_1}$ may not be the same as $\mathcal{P}^{\mathscr{L}_1}_{\mathbb{T}^*}$.
\end{Remark}

Given the observed gene trees $(\mathbb{T}_1, \mathbb{T}_2, \ldots, \mathbb{T}_k)$, the ML supertree is defined as
\[
\hat \mathbb{T}_k = \argmax_{\mathbb{T} \in \mathcal{T}}  \ell_k(\mathbb{T}) 
\]
where $\ell_k(\mathbb{T}) $ denotes the log-likelihood function
\[
\ell_k(\mathbb{T}) = \sum_{i=1}^k \log \mathcal{P}^{\mathscr{L}_i}_{\mathbb{T}}(\mathbb{T}_i).
\]

To enable theoretical analyses of the ML supertree method, we make the following assumptions.

\begin{Assumption}[Weak covering property]
 The sequence of subsets of $\mathscr{L}^*$ satisfies the \emph{weak covering property}: there exists $c_1>0, \gamma>1/2$, and $K > 0$ such that for each subset $\mathscr{L}$ of taxa from $\mathscr{L}^*$ of size $4$,
\[
\frac{1}{k^\gamma}|\{i \le k: \mathscr{L} \subset \mathscr{L}_i\}| \ge c_1 \hspace*{.24in} \forall k \geq K.
\]
Here, $|\mathcal{A}|$ denotes the number of elements in the set $\mathcal{A}$.
\label{as:cover}
\end{Assumption}

We note that the weak covering property ensures that as the number of gene trees increases, all quartets (subtrees with $4$ leaves) of $\mathbb{T}^*$ are visited with enough frequency to enable a reliable estimate for each of the quartets. 
Assumption $\ref{as:cover}$ is a direct generalization of the \emph{covering property} introduced in \cite{steel2008maximum}, which can be obtained from Assumption $\ref{as:cover}$ by setting $\gamma=1$. 

\begin{Assumption}[Model identifiability]
For all $\mathscr{L} \subset \mathscr{L}^*$, the distribution $\mathcal{P}^{\mathscr{L}}_{\mathbb{T}^*}(\cdot)$ is identifiable.
That is, if $d(\mathbb{T}^* |_{\mathscr{L}}, \mathbb{T} |_{\mathscr{L}}) >  0$, then $\text{KL}(\mathcal{P}^{\mathscr{L}}_{\mathbb{T}^*}, \mathcal{P}^{\mathscr{L}}_{\mathbb{T}}) > 0$. 
\label{as:iden}
\end{Assumption}

Assumption $\ref{as:iden}$ guarantees that it is at least possible to reconstruct the restriction $\mathbb{T}^*|_{\mathscr{L}}$ of $\mathbb{T}^*$ to the subset $\mathscr{L}$ with a complete knowledge of $\mathcal{P}^{\mathscr{L}}_{\mathbb{T}^*}$. 
This assumption is similar to, but distinct from, the condition of \emph{basic centrality} introduced by \cite{steel2008maximum}, which requires that for all subsets of  $\mathscr{L} \subset \mathscr{L}^*$  of size 4, 
\[
\mathcal{P}_{\mathbb{T}^*}^{\mathscr{L}}[\mathbb{T}^* |_\mathscr{L}] \ge (1+\eta) \mathcal{P}_{\mathbb{T}^*}^{\mathscr{L}}[\mathbb{T}'] 
\]
for all trees $\mathbb{T}'$ on leaf set $\mathscr{L}$ that are different from $\mathbb{T}^* |_\mathscr{L}$, and where $\eta>0$:
\begin{itemize}
\item On one hand, the basic centrality condition implies that if $d(\mathbb{T}^* |_{\mathscr{L}}, \mathbb{T} |_{\mathscr{L}}) >  0$, the modes of the distributions $\mathcal{P}^{\mathscr{L}}_{\mathbb{T}^*}$ and $\mathcal{P}^{\mathscr{L}}_{\mathbb{T}}$ are different. 
From this, we can deduce that $\text{KL}(\mathcal{P}^{\mathscr{L}}_{\mathbb{T}^*}, \mathcal{P}^{\mathscr{L}}_{\mathbb{T}}) > 0$. 
Thus, our identifiability assumption is somewhat weaker than this condition.
Assumption $\ref{as:iden}$ also does not impose any assumption on the family of distribution themselves and thus can be applied to a wider class of probabilistic models. 

\item On the other hand, the basic centrality condition only concerns leaf sets of size 4, while Assumption $\ref{as:iden}$ impose restrictions on leaf sets of all sizes. 
However, we note that conditions on subset of leaves (such as the basic centrality condition) only work under the implicit assumption that there are some connection between $\mathcal{P}^{\mathscr{L}_2}_{\mathbb{T}}|_{\mathscr{L}_1}$ and  $\mathcal{P}^{\mathscr{L}_1}_{\mathbb{T}}$ for $\mathscr{L}_1 \subset \mathscr{L}_2$.
Since our framework does not assume any nested structure in the probability model, Assumption $\ref{as:iden}$ is more appropriate. 
\end{itemize}

Finally, we impose the following regularity conditions on the tree-generating probability model:

\begin{Assumption}
(Regularity)

\begin{itemize}

\item[(a)] For all $\mathbb{T} \in \mathcal{T}$, $\mathscr{L} \subset \mathscr{L}^*$, and any tree $\mathbb{T}'$ with leaf set $\mathscr{L}$, 
\begin{itemize}
\item[i. ] $\mathcal{P}^{\mathscr{L}}_{\mathbb{T}}(\mathbb{T}') > 0$,
\item[ii. ] $\log \mathcal{P}^{\mathscr{L}}_{\mathbb{T}}(\mathbb{T}')$ is a locally-Lipchitz function with respect to $\mathbb{T}$, and the Lipchitz constant does not depend on $\mathbb{T}'$.
\end{itemize}


\item[(b)] There exist $c_2, c_3 > 0$, $m \ge 2$ such that for any leaf set $\mathscr{L}$ and tree $\mathbb{T}$, if $d(\mathbb{T}^*, \mathbb{T}) \le c_2$, then 
\[
\text{KL}(\mathcal{P}^{\mathscr{L}}_{\mathbb{T}^*}, \mathcal{P}^{\mathscr{L}}_{\mathbb{T}}) \ge c_3 d(\mathbb{T}^* |_{\mathscr{L}}, \mathbb{T} |_{\mathscr{L}})^m.
\] 
\end{itemize}
\label{as:regular}
\end{Assumption}

\begin{Remark}
If for each tree $\mathbb{T}'$ with leaf set $\mathscr{L}$, the probability density function $\mathcal{P}^{\mathscr{L}}_{\mathbb{T}}(\mathbb{T}')$ is an analytic function with respect to $\mathbb{T}$ in a neighborhood of $\mathbb{T}^*$, then Assumption $\ref{as:regular}$(b) holds.
\label{rem:smooth}
\end{Remark}

\begin{proof}
Note that the function $\mathbb{E}_{\mathcal{P}_{\mathbb{T}^*}}\left [\log \mathcal{P}^{\mathscr{L}}_{\mathbb{T}}(\mathbb{T}')\right]$ (with respect to $\mathbb{T}$) is analytic in a neighborhood $\mathcal{U}$ of $\mathbb{T}^*$. 
Define 
\[
\mathcal{A} = \{ \mathbb{T} \in \mathcal{U}:  \mathbb{E}_{\mathbb{T}' \sim \mathcal{P}_{\mathbb{T}^*}}\left [\log \mathcal{P}^{\mathscr{L}}_{\mathbb{T}^*}(\mathbb{T}')\right] = \mathbb{E}_{\mathbb{T}' \sim \mathcal{P}_{\mathbb{T}^*}}\left [\log \mathcal{P}^{\mathscr{L}}_{\mathbb{T}}(\mathbb{T}')\right] \}.
\]
For all $\mathbb{T} \in \mathcal{A}$, we have 
\[
\text{KL}(\mathcal{P}^{\mathscr{L}}_{\mathbb{T}^*}, \mathcal{P}^{\mathscr{L}}_{\mathbb{T}}) = \mathbb{E}_{\mathbb{T}' \sim \mathcal{P}_{\mathbb{T}^*}}\left [\log \mathcal{P}^{\mathscr{L}}_{\mathbb{T}^*}(\mathbb{T}')\right] - \mathbb{E}_{\mathbb{T}' \sim \mathcal{P}_{\mathbb{T}^*}}\left [\log \mathcal{P}^{\mathscr{L}}_{\mathbb{T}}(\mathbb{T}')\right] = 0.
\]
By Assumption \ref{as:iden}, we conclude that $\mathbb{T} |_\mathscr{L} = \mathbb{T}^* |_\mathscr{L}$ for all $\mathbb{T} \in \mathcal{A}$.
Applying {\L}ojasiewicz inequality \citep[][Theorem 1]{ji1992global}, we deduce that there exists $c_3 >0$ and $m \ge 2$ such that for any $\mathbb{T}$ in the neighborhood,
\[
\mathbb{E}_{\mathcal{P}_{\mathbb{T}^*}}\left [\log \mathcal{P}^{\mathscr{L}}_{\mathbb{T}^*}(\mathbb{T}')\right] - \mathbb{E}_{\mathcal{P}_{\mathbb{T}^*}}\left [\log \mathcal{P}^{\mathscr{L}}_{\mathbb{T}}(\mathbb{T}')\right] \ge  c_3 d(\mathbb{T}, \mathcal{A})^m \geq c_3 d(\mathbb{T}^* |_{\mathscr{L}}, \mathbb{T} |_{\mathscr{L}})^m.
\]
This implies the result.
\end{proof}


Throughout this paper, we will assume that Assumptions $\ref{as:cover}$, $\ref{as:iden}$, and $\ref{as:regular}$ hold.
In the next section, we will establish the following convergence rate of the ML supertree.

\begin{thm}
Under Assumptions $\ref{as:cover}$, $\ref{as:iden}$, and $\ref{as:regular}$, for any $\delta>0$, there exist constants $m>0$, $ C_{\delta} > 0$ and $K_\delta>0$ such that for all $k \ge K_\delta$, we have
\[
d(\hat \mathbb{T}_k , \mathbb{T}^*)\leq C_{\delta} \left( \frac{\log k}{k^{(\gamma-1/2)}}\right)^{1/m}
\]
with probability at least $1-\delta$.
\label{thm:consistency}
\end{thm}

Theorem \ref{thm:consistency} establishes that under fairly mild regularity conditions, the ML supertree is consistent and has a polynomial convergence rate. 
Notably, the result holds for all identifiable family of locally-analytic distributions (Remark \ref{rem:smooth}). 
The degree of this (polynomial) convergence rate depends on the sampling scheme of the leaves (characterized by the covering coefficient $\gamma$) and on the geometry of the probabilistic model (characterized by the constant $m$ in Assumption $\ref{as:regular}$(b)). 
We further note that Assumption $\ref{as:regular}$(b)) is only required for establishing the convergence rate of the ML estimator and the absence of this condition does not affect the proof of consistency.

\begin{cor}
Under Assumptions $\ref{as:cover}$, $\ref{as:iden}$, and $\ref{as:regular}(a)$, the ML supertree is consistent.
\end{cor}


\begin{Remark}
Theorem $\ref{thm:consistency}$ is still valid if the universal Assumption $\ref{as:regular}$(b) is replaced by the following condition:

\emph{
There exist $c_2, c_3 > 0$, $m \ge 2$ such that if $d(\mathbb{T}^*, \mathbb{T}) \le c_2$, then 
\[
\text{KL}(\mathcal{P}^{\mathscr{L}_i}_{\mathbb{T}^*}, \mathcal{P}^{\mathscr{L}_i}_{\mathbb{T}}) \ge c_3 d(\mathbb{T}^* |_{\mathscr{L}_i}, \mathbb{T} |_{\mathscr{L}_i})^m \hspace*{.24in} \forall i \in \mathbb{N}.
\] 
}
\label{rem:reg}
\end{Remark}

\section{Convergence of maximum likelihood supertree}

To enable the analysis of convergence, we define
\[
R_{k}(\mathbb{T}) = \frac{1}{k}\ell_k(\mathbb{T}^*) - \frac{1}{k}\ell_k(\mathbb{T}) 
\]
and 
\[
\mathcal{R}_k(\mathbb{T}) := E_{\mathcal{P}_{\mathbb{T}^*}}[R_{k}(\mathbb{T})] =  \frac{1}{k} \left(E_{\mathcal{P}_{\mathbb{T}^*}}[\ell_k(\mathbb{T}^*)]  - E_{\mathcal{P}_{\mathbb{T}^*}}[\ell_k(\mathbb{T})]  \right) = \frac{1}{k} \text{KL}(\mathcal{P}_{\mathbb{T}^*}, \mathcal{P}_{\mathbb{T}}).
\]
We refer to $R_{k}$ and $E_{\mathcal{P}_{\mathbb{T}^*}}[R_{k}(\mathbb{T})]$ as the \emph{empirical risk function} and the \emph{expected risk function}, respectively. 

An intuitive argument for the consistency of the ML estimator can be described as follows. 
As the number of gene trees increases, we have
\[
\left |R_{k}(\mathbb{T})  - \mathcal{R}_k(\mathbb{T}) \right|  \to 0 
\]
with sufficiently high probability. 
Therefore, the ML estimator will converge to the optimal value of the risk function, which is attained at the true species tree $\mathbb{T}^*$.
This simple argument is formalized by a lower bound of the expected risk function (Section 3.1) and a uniform concentration bound on the deviation of the empirical risk function and its expectation (Section 3.2). 

\subsection{Lower bound of the expected risk}

\begin{Lemma}

There exist a neighborhood $\mathcal{U}$ of $\mathbb{T}^*$ and $c_4>0$ such that 
\[
\mathbb{E}_{\mathcal{P}_{\mathbb{T}^*}}\left [\ell_{k}(\mathbb{T}^*) - \ell_{k}(\mathbb{T})  \right] \ge c_4 k^{\gamma} d(\mathbb{T}^*, \mathbb{T})^m \hspace*{.24in} \forall \mathbb{T} \in \mathcal{U}.
\]
\label{info}
\end{Lemma}

\begin{proof}

We note that
\[
\mathbb{E}_{\mathcal{P}_{\mathbb{T}^*}}\left [\ell_{k}(\mathbb{T}^*) - \ell_{k}(\mathbb{T})  \right]  = \text{KL}(\mathcal{P}_{\mathbb{T}^*}, \mathcal{P}_{\mathbb{T}}).
\]
Consider an arbitrary leaf set $\mathscr{L} \subset \mathscr{L}^*$ of size 4 and define $\mathcal{I}= \{i \le k: \mathscr{L} \subset \mathscr{L}_i\}$.
By the weak covering property, we have $|\mathcal{I}| \ge c_1 k^{\gamma}$.
Thus,
\[
 \text{KL}(\mathcal{P}_{\mathbb{T}^*}, \mathcal{P}_{\mathbb{T}})  \ge \sum_{i \in \mathcal{I}}{\text{KL}(\mathcal{P}^{\mathscr{L}_i}_{\mathbb{T}^*}, \mathcal{P}^{\mathscr{L}_i}_{\mathbb{T}})} \geq \sum_{i \in \mathcal{I}}{c_3 d(\mathbb{T}^*|_{\mathscr{L}_i}, \mathbb{T}|_{\mathscr{L}_i})^m} \ge c_1 c_3 k^{\gamma} (2n-3)^{m/2} d(\mathbb{T}^*|_\mathscr{L}, \mathbb{T}|_\mathscr{L})^m
\]
for all $\mathbb{T}$ in some neighborhood $U$ around $\mathbb{T}^*$. 
Here, the second inequality comes from Assumption $\ref{as:regular}$(b).

On the other hand, \cite{dinh2018consistency} (Lemma 6.2 (i)) proved that for some leaf set $\mathscr{L} \subset \mathscr{L}^*$ of size 4, we have
\[
d(\mathbb{T}^*|_\mathscr{L}, \mathbb{T}|_\mathscr{L}) \ge \frac{1}{(2n-3)} d(\mathbb{T}^*, \mathbb{T}).
\]
We deduce that
\[
\text{KL}(\mathcal{P}_{\mathbb{T}^*}, \mathcal{P}_{\mathbb{T}}) \ge \frac{1}{(2n-3)^{m/2}} c_1 c_3 k^{\gamma} d(\mathbb{T}^*, \mathbb{T})^m
\]
\end{proof}

\begin{Lemma}
For any leaf set $\mathscr{L}$ and any $x>0$, there exists $C_x>0$ such that
\[
 \text{KL}(\mathcal{P}^{\mathscr{L}}_{\mathbb{T}^*}, \mathcal{P}^{\mathscr{L}}_{\mathbb{T}}) \ge C_x 
\]
when $d(\mathbb{T}^*|_\mathscr{L}, \mathbb{T}|_\mathscr{L}) \ge x$.
\label{prelimc}
\end{Lemma}

\begin{proof}
Assume that there exists a sequence of tree $\{\mathbb{T}_i\}$ such that $d(\mathbb{T}^*|_\mathscr{L}, \mathbb{T}_i|_\mathscr{L}) \ge x$ and $ \text{KL}(\mathcal{P}^{\mathscr{L}}_{\mathbb{T}^*}, \mathcal{P}^{\mathscr{L}}_{\mathbb{T}_i}) \to 0$. 
Since the tree space is compact, we can extract a sub-sequence $\{\mathbb{T}_{i_j}\}$ of $\{\mathbb{T}_i\}$ that converges to some tree $\mathbb{T}_0$. 
We deduce that $ \text{KL}(\mathcal{P}^{\mathscr{L}}_{\mathbb{T}^*}, \mathcal{P}^{\mathscr{L}}_{\mathbb{T}_0}) =0$ and $\mathbb{T}^*|_\mathscr{L} \ne \mathbb{T}_0|_\mathscr{L}$. 
This contradicts Assumption $\ref{as:iden}$.
\end{proof}


\begin{Lemma}
Let $\mathcal{V}$ be a neighborhood of $\mathbb{T}^*$. 
There exist $c_{\mathcal{V}} > 0$ such that for any $\mathbb{T} \not \in \mathcal{V}$, we have 
\[
\text{KL}(\mathcal{P}_{\mathbb{T}^*}, \mathcal{P}_{\mathbb{T}}) \ge c_{\mathcal{V}} k^{\gamma}.
\] 
\label{prelimb}
\end{Lemma}

\begin{proof}
By \cite{dinh2018consistency} (Lemma 6.2 (i)), there exists $\mathscr{L} \subset \mathscr{L}^*$ of size 4 such that
\[
d(\mathbb{T}^*|_\mathscr{L}, \mathbb{T}|_\mathscr{L}) \ge \frac{1}{(2n-3)} d(\mathbb{T}^*, \mathbb{T}) \ge C_{\mathcal{V}} >0
\]
since $\mathbb{T} \not \in \mathcal{V}$. 
Define $\mathcal{I}= \{i \le k: \mathscr{L} \subset \mathscr{L}_i\}$, we note that for all $i \in \mathcal{I}$
\[
d(\mathbb{T}^*|_{\mathscr{L}_i}, \mathbb{T}|_{\mathscr{L}_i}) \ge (2n-3)^{1/2} d(\mathbb{T}^*|_\mathscr{L}, \mathbb{T}|_\mathscr{L}) \ge C'_{\mathcal{V}},
\]
which implies $ \text{KL}(\mathcal{P}^{\mathscr{L}_i}_{\mathbb{T}^*}, \mathcal{P}^{\mathscr{L}_i}_{\mathbb{T}}) \ge C''_{\mathcal{V}}$ for some $ C''_{\mathcal{V}} >0$ by Lemma $\ref{prelimc}$. 

By the weak covering property, we have $|\mathcal{I}| \ge c_1 k^{\gamma}$.
Thus,
\[
 \text{KL}(\mathcal{P}_{\mathbb{T}^*}, \mathcal{P}_{\mathbb{T}})  \ge \sum_{i \in \mathcal{I}}{\text{KL}(\mathcal{P}^{\mathscr{L}_i}_{\mathbb{T}^*}, \mathcal{P}^{\mathscr{L}_i}_{\mathbb{T}})} \ge c_1 k^{\gamma} C''_{\mathcal{V}}
\]
which completes the proof.
\end{proof}



\subsection{Uniform concentration bound}

\begin{Lemma}[Concentration bound]
For any $\delta>0$, $k \ge 3$, there exists $c_5(\delta)$ such that
\[
\left |R_{k}(\mathbb{T})  - \mathbb{E}_{\mathcal{P}_{\mathbb{T}^*}}\left [R_{k}(\mathbb{T}) \right] \right| \le \frac{c_5 \log k} {\sqrt{k}} \hspace*{.24in} \forall \mathbb{T}, 
\]
with probability at least $1-\delta$.
\label{lem:variance}
\end{Lemma}


\begin{proof}
Since the functions $\log \mathcal{P}^{\mathscr{L}_i}_{\mathbb{T}}(\mathbb{T}_i)$ is locally Lipschitz with respect to $\mathbb{T}$ (Assumption $\ref{as:regular}$(a)) and number of species is finite, there exists $C_1 > 0$ such that
\[
|\log \mathcal{P}^{\mathscr{L}_i}_{\mathbb{T}}(\mathbb{T}_i) -  \log \mathcal{P}^{\mathscr{L}_i}_{\mathbb{T}'}(\mathbb{T}_i)| \le C_1 d(\mathbb{T}, \mathbb{T}') \hspace*{.24in} \forall i, \mathbb{T}_i, \mathbb{T}, \mathbb{T}'.
\]
Since the tree space is compact, $C_1 d(\mathbb{T}, \mathbb{T}')$ is bounded by a constant $C_2$.
Using Hoeffding's inequality \citep{hoeffding1963probability}, we obtain
\[
\mathbb{P}\left[ \left|R_{k}(\mathbb{T})  - \mathbb{E}_{\mathcal{P}_{\mathbb{T}^*}}\left [R_{k}(\mathbb{T}) \right]   \right|  \ge \frac{x}{\sqrt{k}} \right] \le 2\exp\left(-\frac{2x^2}{kC_2^2}\right).
\]
On the other hand, we have $|R_{k}(\mathbb{T}) - R_k(\mathbb{T}')| \le C_1 d(\mathbb{T}, \mathbb{T}')$ and $|\mathbb{E}_{\mathcal{P}_{\mathbb{T}^*}}\left [R_{k}(\mathbb{T}) \right] - \mathbb{E}_{\mathcal{P}_{\mathbb{T}^*}}\left [R_{k}(\mathbb{T}') \right]| \leq C_1 d(\mathbb{T}, \mathbb{T}')$.
Thus, if we define the events
\[
\mathcal{A}(x, k, \mathbb{T})=\left\{ \left|R_{k}(\mathbb{T})  - \mathbb{E}_{\mathcal{P}_{\mathbb{T}^*}}\left [R_{k}(\mathbb{T}) \right]   \right|  \ge \frac{x}{2\sqrt{k}} \right\}
\]
and
\[
\mathcal{B}(x, k, \mathbb{T})=\left\{ \exists \mathbb{T}': d(\mathbb{T}, \mathbb{T}') \le \frac{x}{4C_1\sqrt{k}} ~ \text{and} ~ \left|R_{k}(\mathbb{T}')  - \mathbb{E}_{\mathcal{P}_{\mathbb{T}^*}}\left [R_{k}(\mathbb{T}') \right]   \right|  \ge \frac{x}{\sqrt{k}} \right\}
\]
then $\mathcal{B}(x, k, \mathbb{T}) \subset \mathcal{A}(x, k, \mathbb{T})$.
Note that the total number of balls of radius $x/(4C_1\sqrt{k})$ required to cover the tree space is bounded above by
\[
C_n \left(\frac{4g_0C_1 \sqrt{k}}{x}\right)^{2n-3} (2n-3)!!
\]
where $g_0$ is the upper bound of branch lengths and $C_n$ is a constant depending on the number of species.
To obtain the desired inequality, we will chose $x$ such that
\[
C_n \left(\frac{4g_0C_1 \sqrt{k}}{x}\right)^{2n-3} (2n-3)!! \times 2\exp\left(-\frac{2x^2}{kC_2^2}\right) \le \delta 
\] 
which can be done with $x = C(\delta, n, C_2, g_0) \log k$.
\end{proof}

\subsection{Proof of Theorem \ref{thm:consistency}} 



First, we establish that the ML estimator is consistent. 

By Lemma $\ref{lem:variance}$, we have
\[
\mathcal{R}(\hat \mathbb{T}_k) \le  R_k(\hat \mathbb{T}_k) + c_5 \frac{ \log(k)}{\sqrt{k}} \le  c_5 \frac{ \log(k)}{\sqrt{k}} 
\]
with probability at least $1-\delta$, since $\hat \mathbb{T}_k$ is the maximizer of the empirical risk function. 

On the other hand, let $\mathcal{V}$ be a neighborhood of $\mathbb{T}^*$, by Lemma \ref{prelimb},  we have
\[
\mathbb{E}_{\mathcal{P}_{\mathbb{T}^*}} [R_k(\mathbb{T}')] \ge c_{\mathcal{V}} \frac{1}{k^{1-\gamma}}\hspace*{.24in} \forall \mathbb{T}' \not \in \mathcal{V}.
\]
Since $\gamma>1/2$, there exists $K_{\delta, \mathcal{V}}$ such that for $k \ge K_{\delta, \mathcal{V}}$, we have 
\[
c_5 \frac{ \log(k)}{\sqrt{k}}   \le c_{\mathcal{V}} \frac{1}{k^{1-\gamma}}
\]
We deduce that for $k \ge K_{\delta, \mathcal{V}}$, $\hat \mathbb{T}_k \in \mathcal{V}$ with probability at least $1 -\delta$. 
This proves that the ML estimator is consistent. 

Next, we will derive the convergence rate of the ML supertree. 
From Lemma \ref{info}, there exist a neighborhood $\mathcal{U}$ of $\mathbb{T}^*$ and $c_4, m>0$ such that 
\[
\mathbb{E}_{\mathcal{P}_{\mathbb{T}^*}}[R_k(\mathbb{T})] = \frac{1}{k}\mathbb{E}_{\mathcal{P}_{\mathbb{T}^*}}\left [\ell_{k}(\mathbb{T}^*) - \ell_{k}(\mathbb{T})  \right] \ge c_4 \frac{1}{k^{1-\gamma}}\ d(\mathbb{T}^*, \mathbb{T})^m \hspace*{.24in} \forall \mathbb{T} \in \mathcal{U}.
\]
For $k \ge K_{\delta, \mathcal{U}}$, we have $\hat \mathbb{T}_k \in \mathcal{U}$ with probability at least $1 -\delta$.
By Lemma \ref{lem:variance}, with probability $1 - 2\delta$, we obtain
\begin{align*}
c_4 \frac{1}{k^{1-\gamma}}\ d(\mathbb{T}^*, \hat \mathbb{T}_k)^m & \le \mathcal{R}_k(\hat \mathbb{T}_k) -  R_k(\hat \mathbb{T}_k)   \\
& \le c_5 \frac{ \log(k)}{\sqrt{k}}.
\end{align*}
Here, 
\[
R_k(\hat \mathbb{T}_k) = \frac{1}{k}\ell_k(\mathbb{T}^*) - \frac{1}{k}\ell(\hat \mathbb{T} _k) \le 0
\]
because $\hat \mathbb{T}_k$ is the maximizer of $R_k$. 
This completes the proof.

\section{Applications: convergence under the exponential model}

The \emph{exponential model} is a simple model of gene tree estimation errors in which the probability of observing a given tree decreases exponentially with its distance from the species tree \citep{steel2008maximum}. 
Suppose $d$ is some metric on the set of trees, in the exponential model, the probability of reconstructing any tree $\mathbb{T}'$ with a leaf set $\mathscr{L}$, when $\mathbb{T}^*$ is the generating tree is proportional to an exponentially decaying function of the distance from $\mathbb{T}'$ to $\mathbb{T}^*|_{\mathscr{L}}$:
\[
\mathcal{P}^{\mathscr{L}}_{\mathbb{T}^*}(\mathbb{T}') = \alpha_{\mathbb{T}^*, \mathscr{L}} \exp (-\beta_{\mathscr{L}} d(\mathbb{T}', \mathbb{T}^* |_{\mathscr{L}})).
\] 
Here $\beta_{\mathscr{L}}$ is a constant that depends only on the set of leaves $\mathscr{L}$, while $\alpha_{\mathbb{T}^*, \mathscr{L}}$ is the normalizing constant to ensure that $\mathcal{P}^{\mathscr{L}}_{\mathbb{T}^*}$ is a density function.

In this section, we verify the identifiability and regularity conditions for the exponential model when $d$ is any continuous tree distance such that if $\mathbb{T}$ and $\mathbb{T}'$ have the same topology, then $d(\mathbb{T}, \mathbb{T}')$ is the Euclidean distance.

\begin{thm}
Under the exponential model with Assumptions $\ref{as:cover}$, for any $\delta>0$, $ C_{\delta} > 0$ and $K_\delta > 0$ such that for all $k \ge K_\delta$,
\[
d(\hat \mathbb{T}_k , \mathbb{T}^*)\leq C_{\delta} \left( \frac{\log k}{k^{(\gamma-1/2)}}\right)^{1/m}
\]
with probability at least $1-\delta$, where 
\[
m = 4\sup_{i}{|\mathscr{L}_i|} -4. 
\]
\label{thm:exp}
\end{thm}

\begin{proof}

First, we note that for all $\mathbb{T}'$ with leaf set $\mathscr{L}$, $\mathcal{P}^{\mathscr{L}}_{\mathbb{T}}(\mathbb{T}')$ is a positive and locally Lipschitz function. 
To obtain an upper bound on the convergence rate for the ML supertree, we need to verify Assumption $\ref{as:iden}$ and Assumption $\ref{as:regular}(b)$/Remark $\ref{rem:reg}$.

\smallskip

\textbf{Identifiability}. Since it is sufficient to prove either $\text{KL}(\mathcal{P}^{\mathscr{L}}_{\mathbb{T}^*}, \mathcal{P}^{\mathscr{L}}_{\mathbb{T}}) > 0$ or $\text{KL}(\mathcal{P}^{\mathscr{L}}_{\mathbb{T}}, \mathcal{P}^{\mathscr{L}}_{\mathbb{T}^*}) > 0$, we can assume that $\alpha_{\mathbb{T}^*, \mathscr{L}} \geq \alpha_{\mathbb{T}, \mathscr{L}}$. 

We denote 
\[
\mathcal{A}_c = \{ \mathbb{T}' : \beta_{\mathscr{L}} d(\mathbb{T}', \mathbb{T}^* |_{\mathscr{L}}) \leq  c \}
\] 
and pick $c$ small enough such that $\beta_{\mathscr{L}} d(\mathbb{T}', \mathbb{T} |_{\mathscr{L}}) \ge 2c$.
We note that in $A_c$,
\[
 \exp \left(-\beta_{\mathscr{L}} d(\mathbb{T}', \mathbb{T}^* |_{\mathscr{L}}) \right)  \ge \exp \left(-\beta_{\mathscr{L}} d(\mathbb{T}', \mathbb{T} |_{\mathscr{L}}) \right).
\]
Since $\alpha_{\mathbb{T}^*, \mathscr{L}} \geq \alpha_{\mathbb{T}, \mathscr{L}}$, we have
\begin{align*}
 \alpha_{\mathbb{T}^*, \mathscr{L}} \exp \left(-\beta_{\mathscr{L}} d(\mathbb{T}', \mathbb{T}^* | \mathscr{L}) \right) &-  \alpha_{\mathbb{T}, \mathscr{L}} \exp \left(-\beta_{\mathscr{L}} d(\mathbb{T}', \mathbb{T} | \mathscr{L}) \right)\\
 &\ge  \alpha_{\mathbb{T}, \mathscr{L}} \left[ \exp \left(-\beta_{\mathscr{L}} d(\mathbb{T}', \mathbb{T}^* | \mathscr{L}) \right)  -\exp \left(-\beta_{\mathscr{L}} d(\mathbb{T}', \mathbb{T}^* | \mathscr{L}) \right)  \right].
\end{align*}
By Pinsker's inequality,
\begin{align*}
\text{KL}(\mathcal{P}^{\mathscr{L}}_{\mathbb{T}^*}, \mathcal{P}^{\mathscr{L}}_{\mathbb{T}}) &\ge 2d_{TV}(\mathcal{P}^{\mathscr{L}}_{\mathbb{T}^*}, \mathcal{P}^{\mathscr{L}}_{\mathbb{T}}) ^2 \\
& \geq 2  \alpha_{\mathbb{T}, \mathscr{L}}^2  \left [ \int_{\mathcal{A}_c}{ \exp \left(-\beta_{\mathscr{L}} d(\mathbb{T}', \mathbb{T}^* | \mathscr{L}) \right)- \exp \left(-\beta_{\mathscr{L}} d(\mathbb{T}', \mathbb{T} | \mathscr{L}) \right)}  \right ]^2 \\
& \geq  2 \alpha_{\mathbb{T}, \mathscr{L}}^2 \mu(\mathcal{A}_c)^2 (e^c - e^{2c})^2>0,
\end{align*}
which establishes the identifiability of the exponential model.
Here, $\mu$ is the Lebesgue measure.


\bigskip

\textbf{Regularity}. 
Consider a fixed tree $\mathbb{T} \in \mathcal{T}$ and leaf set $\mathscr{L}$, we first assume that $\alpha_{\mathbb{T}^*, \mathscr{L}} \geq \alpha_{\mathbb{T}, \mathscr{L}}$.
If we define
\[
\mathcal{A} = \{ \mathbb{T}' : \beta_{\mathscr{L}} d(\mathbb{T}', \mathbb{T}^* |_{\mathscr{L}}) \leq  \frac{1}{3}d(\mathbb{T} |_{\mathscr{L}}, \mathbb{T}^* |_{\mathscr{L}})) \},
\] 
then
\[
\beta_{\mathscr{L}} d(\mathbb{T}', \mathbb{T} |_{\mathscr{L}}) \ge  \frac{2}{3}d(\mathbb{T} |_\mathscr{L}, \mathbb{T}^* |_{\mathscr{L}}) ~~~ \forall \mathbb{T}' \in \mathcal{A}.
\] 
Using the same argument as above, we have
\begin{align*}
d_{TV}(\mathcal{P}^{\mathscr{L}}_{\mathbb{T}^*}, \mathcal{P}^{\mathscr{L}}_{\mathbb{T}}) & \geq  \alpha_{\mathbb{T}, \mathscr{L}_i} \mu(A) \left [e^{-\frac{1}{3} d(\mathbb{T} |_\mathscr{L}, \mathbb{T}^* |_{\mathscr{L}})} - e^{-\frac{2}{3}d(\mathbb{T} |_\mathscr{L}, \mathbb{T}^* |_{\mathscr{L}})} \right].
\end{align*}
It can be verified that there exists $C > 0$ such that for all $x \in [0, C]$, we have
\begin{equation}
e^{-x} - e^{-2x} \ge \frac{x}{2}.
\label{eq:ex}
\end{equation}
Moreover since $d$ is Euclidean inside each orthant, if $d(\mathbb{T}, \mathbb{T}^*)$ is small enough, we have
\begin{equation}
\mu(A) = C_{2|\mathscr{L}| - 3} \left(\frac{d(\mathbb{T} |_\mathscr{L}, \mathbb{T}^* |_{\mathscr{L}})}{ 3 \beta_{\mathscr{L}}}\right)^{2 |\mathscr{L}| -3}.
\label{eq:mu}
\end{equation}
Thus, let $\mathcal{W}$ be a neighborhood in the same topology of $\mathbb{T}^*$ such that
\[
\alpha_{\mathbb{T}, \mathscr{L}}  \ge \frac{1}{2} \alpha_{\mathbb{T}^*, \mathscr{L}} \hspace*{.24in} \forall \mathbb{T} \in \mathcal{W},
\]
and both Equations $\eqref{eq:ex}$ and $\eqref{eq:mu}$ are satisfied, we have
\[
\text{KL}(\mathcal{P}^{\mathscr{L}}_{\mathbb{T}^*}, \mathcal{P}^{\mathscr{L}}_{\mathbb{T}}) \ge C d(\mathbb{T}, \mathbb{T}^* |_{\mathscr{L}})^{4|\mathscr{L}| - 4} \hspace*{.24in}  \forall \mathbb{T} \in \mathcal{W},
\]
for some constant $C>0$ independent of $\mathbb{T}$. 

For the case when $\alpha_{\mathbb{T}^*, \mathscr{L}} < \alpha_{\mathbb{T}, \mathscr{L}}$ we define
\[
\mathcal{A}'= \{ \mathbb{T}' : \beta_{\mathscr{L}} d(\mathbb{T}', \mathbb{T}^* |_{\mathscr{L}}) \leq  \frac{1}{3}d(\mathbb{T}'|_{\mathscr{L}}, \mathbb{T}|_{\mathscr{L}})) \},
\] 
and the argument proceeds similarly. 
This validates Assumption $\ref{as:regular}(b)$ for exponential model with $m=4n-4$. 
However, if we use the regularity condition in Remark $\ref{rem:reg}$, we can obtain the result with
\[
m = 4\sup_{i}{|\mathscr{L}_i|} -4. 
\]
\end{proof}


\section{Discussion and Conclusion}

In this paper, we propose a novel analytic approach to analyze the convergence of the ML supertree method. 
Instead of focusing on reconstructing the correct discrete topology of the species tree as in previous studies \citep[e.g.][]{roch2015robustness,steel2008maximum}, we employ a continuous model of the tree space and analyze the ML estimator on this metric space, aiming at recovering both the topology and the branch lengths of the species tree. 
This framework enables us to use tools from statistical learning and information theory to establish the convergence rate of the ML estimator and at the same time, to weaken the conditions to obtain consistency and convergence of the estimator.
Our \emph{weak covering property} is an extension of the classical covering property \citep{steel2008maximum} and provides a considerable relaxation on the sampling schemes for supertree estimation.
Our identifiability condition is also more intuitive and generalizable than the well-known basic centrality condition and does not impose constraints on the shape of the probabilistic model of gene tree estimation errors. 
Our information-theoretical approach to analyze statistical estimator on tree spaces is of independent interest and can be extended to other problems in phylogenetics.

There are several avenues for future directions for this work. 
The first direction is extending our results to other practical models of phylogenetic errors, including the multiple-coalescent model (along with a detailed model of the effects of short sequence length on the accuracy in estimating the individual gene trees).
Second, while our result provides a polynomial bound on the convergence rate, the power of the convergence (characterized by the geometric constant $m$) is not sharp. 
A sharper bound of the convergence rate would be of great interest to the field (from both theoretical and applied perspective) and would require further understanding of the tree-generating probabilistic model. 

\section*{Acknowledgement}
LSTH was supported by startup funds from Dalhousie University, the Canada Research Chairs program, the
NSERC Discovery Grant RGPIN-2018-05447, and the NSERC Discovery Launch Supplement DGECR-2018-00181.
VD was supported by a startup fund from University of Delaware and National Science Foundation grant DMS-1951474.

\bibliographystyle{chicago}
