A tensorial approach to the inversion of group-based phylogenetic models

Abstract

Background

Hadamard conjugation is part of the standard mathematical armoury in the analysis of molecular phylogenetic methods. For group-based models, the approach provides a one-to-one correspondence between the so-called “edge length” and “sequence” spectrum on a phylogenetic tree. The Hadamard conjugation has been used in diverse phylogenetic applications not only for inference but also as an important conceptual tool for thinking about molecular data leading to generalizations beyond strictly tree-like evolutionary modelling.

Results

For general group-based models of phylogenetic branching processes, we reformulate the problem of constructing a one-one correspondence between pattern probabilities and edge parameters. This takes a classic result previously shown through use of Fourier analysis and presents it in the language of tensors and group representation theory. This derivation makes it clear why the inversion is possible, because, under their usual definition, group-based models are defined for abelian groups only.

Conclusion

We provide an inversion of group-based phylogenetic models that can implemented using matrix multiplication between rectangular matrices indexed by ordered-partitions of varying sizes. Our approach provides additional context for the construction of phylogenetic probability distributions on network structures, and highlights the potential limitations of restricting to group-based models in this setting.

Background

Fundamental to evolutionary biology is the development and implementation of molecular phylogenetic methods [1]. These methods provide the means to reconstruct the past evolutionary history of biological entities given present-day molecular data, such as DNA. Considering Kimura’s neutral theory of molecular evolution, it is logical to apply a stochastic model at the level of DNA substitutions to construct probabilistic description of what molecular alignments are expected to be observed, given a proposed evolutionary history (tree topology and edge lengths). is commonly implemented assuming an IID (across sites in the alignment) and Markov process for DNA substitution, leading to a model that has a continuous-time Markov chain at its core (see Semple and Steel [2] for an introduction to the mathematics underlying modern phylogenetic methodology).

In a series of papers, Hendy and colleagues introduced the Hadamard conjugation as a novel tool for phylogenetic analyses [3]–[5]. They found an invertible relationship between a phylogenetic tree, as characterized by its edge length spectrum, and the probability distribution of site patterns (referred to as the sequence spectrum). Originally introduced only for the 2-state symmetric model, the Hadamard conjugation was later extended to the K3ST model [6]–[8] and further to any of the so-called “group-based” models [9]. Hadamard conjugation has been used as both a tool for simulation [10] and to look at statistical properties of methods, exploring the inconsistency of parsimony under a molecular clock [5],[11]. For these sorts of applications, following the notation in Felsenstein [1], we can use the Hadamard transform H to start with an edge length spectrum γ and calculate the sequence spectrum s=H ⁻¹ log(H γ). The beauty of Hadamard conjugations is that one can also begin with an observed sequence spectrum $ŝ$ and perform the inverse of the conjugation to empirically obtain an edge length spectrum $\widehat{\gamma }=H^{-1}log\left(\mathrm{Hŝ}\right)$. Although it is not expected that the $\widehat{\gamma }$ spectrum will precisely match a tree, Hendy [12] proposed using an optimisation criterion to map from $\widehat{\gamma }$ to the “closest tree”.

Several authors have commented that it is potentially a useful feature of Hadamard conjugation that data isn’t forced onto a fixed tree. The conflicting information can be retained and interpreted in the form of a “lentoplot” [13] or a splits-graph [14], with both of these methods implemented in Spectronet [15]. Schliep [16] gives some more statistical justification for such an approach by making a link to modern statistical techniques such as the Lasso and Ridge regression.

von Haeseler and Churchill [17] seems to be the first paper that explicitly suggests using Hadamard conjugation to provide a likelihood framework for networks. The principle idea in this work was to start with an edge length spectrum that encodes a set of incompatible splits, use the Hadamard transformation to get site probabilities and use these to determine a likelihood. This idea was further explored by Bryant [18], and Bryant [19] followed this through defining the “n-taxon process” for group-based models. It should be noted that likelihoods calculated via Hadamard are not equivalent to likelihoods calculated by taking a mixture of trees. Indeed, Matsen and Steel [20], Matsen et al. [21] used Hadamard methods in combination with phylogenetic invariants to show that mixtures of trees with the same topology can exactly mimic another tree under the 2-state model. Considering biological applications, thinking in terms of mixtures of trees or partitions where the data can be thought of as arising on a set of trees [22]–[24] seems more reasonable than the Hadamard conjugation. Strimmer and Moulton [25] suggested using split networks as a spring board to likelihood-based analyses on DAGs, but later identified several problems with the approach [26]; most notably, in split-networks internal nodes do not have a biological interpretation as an ancestor.

In Sumner et al. [27], we gave some additional insight into the interpretation of applying the Hadamard conjugation in a network setting. We showed that permutation group structure inherent to the Hadamard transformation – as for any group-based model – restricts the resulting process from being capable of reproducing truly convergent processes. This is a serious limitation, as one of the biological motivations for explicit network models is the ability to model convergent processes. We also presented an alternative algebraic formalism for the general Markov model, analogous to the n-taxon process, but capable of reproducing convergent processes.

From the point of view of group representation theory, the inversion of group-based models relies on the fact that the irreducible representations of an abelian group are one-dimensional, and the model structure essentially reduces to analysing group characters – hence the standard presentation of a Fourier inversion. In this article, we make this connection concrete. For the general Markov model, it is then immediately apparent that an analogous inversion is not possible because the algebraic structure underlying the model is not abelian and hence the irreducible representations are not one-dimensional. In fact, to obtain one-dimensional representations for the general Markov model, it is necessary to apply higher-degree polynomial maps (beyond the degree 1, linear case), and define “Markov invariants” [28]. These invariants present one-dimensional representations but at the cost of the higher degree – degree 5 in the case of the general Markov model with four states on quartet trees [29],[30]. This connection between Hadamard transformation and Markov invariants is an interesting one, but we do not discuss it further here.

In this paper we approach the inversion of group-based phylogenetic models by taking a representation-theoretic perspective and working explicitly with tensor indices. Our approach rests heavily on the formalism of “phylogenetic tensors”, as presented in Bashford et al. [31], for the binary-symmetric and K3ST model, and Sumner et al. [27],[28], for the general Markov model.

Although the main inversion results presented here are not more general than those in in Székely et al. [7], we think it is important to reformulate them using the language of tensors and representation theory. This viewpoint has already led to new approaches for modeling convergent evolution [27] and for studying non-group-based models [28]. However, in none of our previous work was the link to Hadamard conjugation explicitly discussed. By presenting an old technique (Hadamard conjugation) in a new light we hope to introduce other researchers to the viewpoint of tensor analysis and representation theory.

Methods

Group-based models

We consider the continuous-time formulation of Markov processes, and show how to implement the inversion of a group-based phylogenetic model based on any abelian group G. We note that such an inversion requires a map from tensor product space (where elements are indexed by ordered-n-partitions) to phylogenetic splits (where elements are indexed by bipartitions). We achieve this by finding canonical maps from bipartitions to ordered-n-partitions.

For a group G (not necessarily abelian) with order |G|=d, we write G={σ ₁,σ ₂,…,σ _d }, and, when necessary, write ε∈G to specify the identity element of G. We will discuss the “regular representation” of G shortly, but skipping ahead we find that any rate matrix Q occurring in a group-based Markov model can be written in the form

$$\begin{array}{l}\begin{array}{c}Q=-\lambda \mathbf{1}+\sum _{\epsilon \ne \sigma \in G}{\alpha }^{\sigma }{K}_{\sigma },\end{array}\end{array}$$

(1)

where each $0\le {\alpha }^{\sigma }\in ℝ$, $\lambda =\sum _{\epsilon \ne \sigma \in G}{\alpha }^{\sigma }$ and the K _σ are the permutation matrices corresponding to the (non-identity) group elements σ∈G.

For the reader interested in deriving this result, consider the d-dimensional vector space ${〈G〉}_{ℂ}\equiv {〈{\sigma }_1,{\sigma }_2,\dots ,{\sigma }_d〉}_{ℂ}=\{v=v_1{\sigma }_1+v_2{\sigma }_2+\dots +v_d{\sigma }_d:v_i\in ℂ\}$, with scalar multiplication and vector addition defined via

$$\begin{array}{l}\begin{array}{cc}v+\lambda v^{\prime }& =(v_1{\sigma }_1+v_2{\sigma }_2+\dots +v_d{\sigma }_d)\\ +\lambda (v_1^{\prime }{\sigma }_1+v_2^{\prime }{\sigma }_2+\dots +v_d^{\prime }{\sigma }_d)\\ =(v_1+\lambda v_1^{\prime }){\sigma }_1+(v_2+\lambda v_2^{\prime }){\sigma }_2+\dots \\ +(v_d+\lambda v_d^{\prime }){\sigma }_d,\end{array}\end{array}$$

for all $v,v^{\prime }\in {〈G〉}_{ℂ}$ and $\lambda \in ℂ$. The regular representation, ${\rho }_{\text{reg}}:G\to \mathit{\text{GL}}(d,ℂ)$, is then defined by setting the group action

$$\sigma :v\mapsto \mathrm{\sigma v}=v_1\left(\sigma {\sigma }_1\right)+v_2\left(\sigma {\sigma }_2\right)+\dots +v_d\left(\sigma {\sigma }_d\right),$$

for all $v\in {〈G〉}_{ℂ}$ and σ∈G. If we fix {σ ₁,σ ₂,…,σ _d } as an ordered basis for ${〈G〉}_{ℂ}$, it is then clear – via Cayley’s theorem – that each group element σ gets mapped to a permutation matrix K _σ :=ρ _reg(σ), with $K_{\sigma }{\sigma }_i=\sum _j{\left[K_{\sigma }\right]}_i^j{\sigma }_j:=\sigma {\sigma }_i$. Thus K _σ has matrix elements

$$\begin{array}{l}\begin{array}{c}{\left[K_{\sigma }\right]}_i^j=\left\{\begin{array}{c}1,\text{if}{\sigma }_j=\sigma {\sigma }_i,\\ 0,\text{otherwise.}\end{array}\right.\end{array}\end{array}$$

(2)

Consider the unit column vectors

$$\begin{array}{cc}{\xi }_1& ={(1,0,0,\dots ,0)}^T,{\xi }_2={(0,1,0,0,\dots ,0)}^T,\dots \\ {\xi }_d& ={(0,0,\dots ,0,1)}^T;\end{array}$$

and identify each ${\sigma }_i\in {〈G〉}_{ℂ}$ with ${\xi }_i\in {ℂ}^d$, so that the group action becomes σ:ξ _i ↦K _σ ξ _i =ξ _j where σ _j =σ σ _i . Thus the matrix elements ${\left[K_{\sigma }\right]}_i^j$ have i as the column label and j as the row label.

A group-based Markov model is then obtained by taking a continuous-time Markov chain with state space G={σ ₁,σ ₂,…,σ _d } and using the group multiplication in G to assign a rate α _σ to all substitutions σ ₁↦σ ₂ where σ σ ₁=σ ₂. Following this through (as is done in detail in [32]) we are led to the formula (1) for rate matrices in any group-based model.

The regular representation is one example of the general concept of a representation of G on a vector space V, defined as a homomorphism ρ:G→G L(V) satisfying ρ(g ₁g ₂)=ρ(g ₁)ρ(g ₂) for all g ₁,g ₂∈G. A representation is said to be reducible if there exists a proper subspace U⊂V satisfying ρ(g)U⊂U, i.e. the set of matrices ρ(G) send vectors in U back to U. In this case, U is called an invariant subspace. The representation ρ is then called irreducible if V does not contain any invariant subspaces.

The reader should note that the usual construction of a “group-based” model [2] stipulates that G be abelian. Although the construction just given using the regular representation allows for non-abelian G, we will nonetheless only consider the abelian case in this paper, because, as discussed in the introduction, it is only in the abelian case that a (linear) inversion of phylogenetic models is possible. In this case the irreducible representations of G are all one-dimensional [33], and hence the analysis reduces to computations with group characters, as is exploited in the previous approaches using Fourier analysis [9],[34].

Phylogenetic tensors

We denote [d]:={1,2,…,d} as the state space for a continuous-time Markov chain. Consider an n-taxa phylogenetic tree and a d-state phylogenetic pattern distribution ${\left\{p_{i_1{i}_2\dots i_n}\right\}}_{i_1,i_2,\dots ,i_n\in \left[d\right]}$ with the interpretation that $p_{i_1{i}_2\dots i_n}$ is the probability that the observed state at the k ^th leaf on the tree is i _k . As is shown in Sumner and Jarvis [35] and in more detail in Sumner et al. [27], such phylogenetic pattern distributions can be represented abstractly as tensors in the n-fold tensor product space ${\otimes }^n{ℂ}^d:={ℂ}^d\otimes {ℂ}^d\otimes \dots \otimes {ℂ}^d$, as follows. If we choose {ξ ₁,ξ ₂,…,ξ _d } as an ordered basis for ${ℂ}^d$, and ordered basis ${\{{\xi }_{i_1}\otimes {\xi }_{i_2}\otimes \dots \otimes {\xi }_{i_d}\}}_{i_1,i_2,\dots ,i_n\in \left[d\right]}$ for the tensor product space, a “phylogenetic tensor” $P\in {\otimes }^n{ℂ}^d$ is then defined as

$$\begin{array}{l}\begin{array}{c}P=\sum _{i_1,i_2,\dots ,i_n\in \left[d\right]}{p}_{i_1{i}_2\dots i_n}{\xi }_{i_1}\otimes {\xi }_{i_2}\otimes \dots \otimes {\xi }_{i_n}.\end{array}\end{array}$$

For readers who are unfamiliar with tensor products, it is possible to understand the general concept via the definition of the “Kronecker” product of a n×m matrix A and a n ^′×m ^′ matrix B as the n n ^′×m m ^′ matrix given by

We can index the matrix A⊗B with row indicies i ₁j ₁=11,12,…,n n ^′ and column indices j ₁j ₂=11,12,…,m m ^′, i.e. generically ${(A\otimes B)}_{i_1{j}_1,i_2{j}_2}=A_{i_1{i}_2}{B}_{j_1{j}_2}$ and specifically (A⊗B)_12,32=A ₁₃B ₂₂. This point of view is useful if one wants to write out specific matrix representations of tensors, however, in the development that follows will focus heavily on the indexing of tensor components in the various cases discussed.

Suppose $\pi =\sum _{i\in \left[d\right]}{\pi }_i{\xi }_i\in {ℂ}^d$ represents the state distribution of a single taxa, i.e. π _i is the probability that a randomly chosen site in the sequence will be in state i. Now suppose a phylogenetic branching event occurs and the sequence is copied. The corresponding phylogenetic tensor $P=\sum _{i_1,i_2\in \left[d\right]}{p}_{i_1{i}_2}{\xi }_{i_1}\otimes {\xi }_{i_2}$ representing the joint distribution of the two-taxa just after the branching event then has the property that $p_{i_1{i}_2}={\pi }_{i_1}$ if i ₂=i ₁ and is zero otherwise. Thinking in terms of tensor operations, we find that phylogenetic branching events can be generated by a linear operator $\delta :{ℂ}^d\to {ℂ}^d\otimes {ℂ}^d$ determined by δ(π)=P and defined in general using our chosen basis as

$$\begin{array}{l}\begin{array}{cc}\delta \left({\xi }_i\right):={\xi }_i\otimes {\xi }_i,\delta \left(\pi \right)& =\delta \left(\sum _i{\pi }_i{\xi }_i\right)=\sum _i{\pi }_i\delta \left({\xi }_i\right)\\ =\sum _i{\pi }_i{\xi }_i\otimes {\xi }_i.\end{array}\end{array}$$

The remarkable fact for group-based models, central to the present article, is that the permutation matrices “intertwine” particularly simply with the branching operator:

$$\begin{array}{l}\begin{array}{c}\delta \left(K_{\sigma }{\xi }_i\right)=\delta \left({\xi }_{\sigma \left(i\right)}\right)={\xi }_{\sigma \left(i\right)}\otimes {\xi }_{\sigma \left(i\right)}=K_{\sigma }\otimes K_{\sigma }·\delta \left({\xi }_i\right).\end{array}\end{array}$$

Thus, for any rate matrix Q arising from a group-based model, we have (via the linearity of δ):

$$\begin{array}{l}\begin{array}{c}\delta ·Q=\left(-\lambda \mathbf{1}\otimes \mathbf{1}+\sum _{\epsilon \ne \sigma \in G}{\alpha }^{\sigma }{K}_{\sigma }\otimes K_{\sigma }\right)·\mathrm{\delta .}\end{array}\end{array}$$

(3)

We also note that, since Q can be expressed a linear combination of permutation matrices representing elements in a group G, the matrix powers Q ²,Q ³,Q ⁴… will also be expressible as linear combinations of the same permutation matrices (although precise expressions for the relevant coefficients may or may not be easily computable). Together with (3), this implies that, for any substitution matrix e ^Qt arising from matrix exponentiation,

$$\begin{array}{l}\begin{array}{c}\delta ·e^{\mathit{\text{Qt}}}=e^{-\lambda }exp\left(\sum _{\epsilon \ne \sigma \in G}{\alpha }^{\sigma }{K}_{\sigma }\otimes K_{\sigma }\right)·\mathrm{\delta .}\end{array}\end{array}$$

(4)

This relation shows that mathematically, and hence conceptually, “Markov evolution on a single followed by a branching event” can be replaced with “Branching event on a single taxon followed by (correlated) Markov evolution of two taxa.” This equivalence is illustrated in Figure 1, and should be compared to the equivalent discussion of the “n-taxa process” given in [18] and [19].

Figure 1

Markov evolution on a single followed by a branching event (illustrated on the left), is equivalent to a branching event on a single taxon followed by correlated Markov evolution of two taxa (illustrated on the right). Mathematically, this equivalence can be implemented by exploiting the equality given in (4).

In Sumner et al. [27] we showed how to generalise this intertwining action to the case of the general Markov model. Interestingly, for the general Markov model the appropriate intertwining has quite a different structure from what occurs in group-based models, and hence the simplicity of (4) is somewhat misleading in general. We refer the reader to Sumner et al. [27] for more discussion on this point.

Returning to the case of group-based models, for each subset A⊆[n], we define a linear map on ${\otimes }^n{ℂ}^d$ as the tensor product $K_{\sigma }^{\left(A\right)}:=K_{\sigma }^{a_1}\otimes K_{\sigma }^{a_2}\otimes \dots \otimes K_{\sigma }^{a_n}$ where a _i =1 if i∈A and 0 otherwise. For example, if n=5, we have

$$\begin{array}{l}\begin{array}{c}{K}_{\sigma }^{\left(\right\{1,2,4\left\}\right)}=K_{\sigma }\otimes K_{\sigma }\otimes \mathbf{1}\otimes K_{\sigma }\otimes \mathbf{1}.\end{array}\end{array}$$

To develop a phylogenetic tensor on a tree, we root the phylogenetic tree at taxon n, and label edges by subsets ∅≠e⊆[n−1], where i∈e if the path from taxon n to taxon i crosses the edge labelled by e. A five taxon tree with this labelling, is presented in Figure 2. To each edge labelled by ∅≠e⊆[n−1], we assign the rate matrix

$$\begin{array}{l}\begin{array}{c}{Q}_e:=-{\lambda }_e\mathbf{1}+\sum _{\epsilon \ne \sigma \in G}{\alpha }_e^{\sigma }{K}_{\sigma },\end{array}\end{array}$$

Figure 2

A six taxa tree rooted at taxon 6 with edges labelled by subsets of {1,2,3,4,5}.

where each ${\alpha }_e^{\sigma }\ge 0$ is the rate of substitution for all states σ ₁ to σ ₂ satisfying $\sigma ={\sigma }_2{\sigma }_1^{-1}$, and ${\lambda }_e=\sum _{\sigma \in G}{\alpha }_e^{\sigma }$. Each edge is then assigned substitution matrix $M_e=e^{Q_e}$, so that the time parameter for each edge is absorbed into the definition of Q _e .

Now iterating (4) multiple times, Bashford et al. [27],[31] show that any phylogenetic tensor can be written as

$$\begin{array}{l}\begin{array}{c}P=e^{-\lambda }exp\left(\sum _{\varnothing \ne e\subseteq \left[n-1\right],\sigma \in G}\underset{e}{\overset{\sigma }{\alpha }}\underset{\sigma }{\overset{\left(e\right)}{K}}\right)·{\delta }^{n-1}\mathrm{\pi .}\end{array}\end{array}$$

(5)

where $\lambda =\sum _{\varnothing \ne e\subseteq \left[n-1\right]}{\lambda }_e=\sum _{\varnothing \ne e\subseteq \left[n-1\right],\epsilon \ne \sigma \in G}{\alpha }_e^{\sigma }$, and δ ⁿ⁻¹π is the d×d×…×d tensor that represents the “zero edge-length star tree” distribution on n taxa. It is this form of phylogenetic tensors that will do a lot of the heavy lifting in the discussion that follows. The reader should note that under this representation, there is no need for the edge parameters $\left\{{\alpha }_e^{\sigma }:\varnothing \ne e\subseteq \left[n-1\right],\sigma \in G\right\}$ to be chosen to be compatible with a particular tree, hence the possibilities for generalising to non-tree-like or network models, as discussed in the introduction.

The stationary distribution for group-based models is uniform (because the rate matrices are doubly stochastic). In this paper we always assume a stationary distribution, so that:

$$\begin{array}{l}\begin{array}{c}\pi =\frac{1}{d}{(1,1,\dots ,1)}^T,\end{array}\end{array}$$

and δ ⁿ⁻¹π has tensor components

$$\begin{array}{l}\begin{array}{c}{\left[{\delta }^{n-1}\pi \right]}_{i_1{i}_2\dots i_n}=\left\{\begin{array}{cc}\frac{1}{d},& \text{if}{i}_1=i_2=\dots =i_n,\\ 0,& \text{otherwise.}\end{array}\right.\end{array}\end{array}$$

This concludes our discussion of the tensor presentation of phylogenetic probability distributions under group-based models. It is important to note that everything discussed so far works for any group-based model, with no requirement that the underlying group G be abelian.

In what follows, we discuss the inversion of abelian group-based models. We present the simplest case with $G={\mathbb{Z}}_2$; the $G={\mathbb{Z}}_3$ case; the $G={\mathbb{Z}}_2\times {\mathbb{Z}}_2$ case; the general $G={\mathbb{Z}}_r$ case; and finally we discuss the case of any abelian group.

Results

The binary-symmetric case

We begin with the inversion of the so-called “binary-symmetric” model. Consider ${ℂ}^2$ with standard basis

$$\begin{array}{l}\begin{array}{c}\left\{{\xi }_0=\left(\begin{array}{c}1\\ 0\end{array}\right),{\xi }_1=\left(\begin{array}{c}0\\ 1\end{array}\right)\right\}.\end{array}\end{array}$$

As a group-based model, the binary-symmetric model arises by taking the group

$$\begin{array}{l}\begin{array}{c}G:={\mathbb{Z}}_2={\{0,1\}}_{+\text{(mod 2)}}\cong 〈\sigma |{\sigma }^2=\epsilon 〉,\end{array}\end{array}$$

with a generic rate matrix given by

$$\begin{array}{l}\begin{array}{c}Q=\left(\begin{array}{cc}-1& 1\\ 1& -1\end{array}\right)=-\mathbf{1}+K,\end{array}\end{array}$$

where $K=\left(\begin{array}{cc}0& 1\\ 1& 0\end{array}\right)$ is the permutation matrix representing σ in the standard basis.

Now ${\rho }_{\text{reg}}:{\mathbb{Z}}_2\to {\mathbf{M}}_2\left(ℂ\right)$, with σ↦K, is the regular representation of ${\mathbb{Z}}_2$, and the character table of ${\mathbb{Z}}_2$ given in Table 1 is easily recognised to be the Hadamard matrix

$$\begin{array}{l}\begin{array}{c}h=\left(\begin{array}{cc}1& 1\\ 1& -1\end{array}\right).\end{array}\end{array}$$

Table 1 The character table of ${\mathbb{Z}}_2$

As ${\mathbb{Z}}_2$ is an abelian group, the irreducible representations are one-dimensional.

The corresponding projection operators can be read off from the columns of the character table. That is, the operators

$$\begin{array}{l}\begin{array}{c}{\Theta }_{\mathit{\text{id}}}:=\frac{1}{2}\left(\epsilon +\sigma \right),{\Theta }_{\mathit{\text{sgn}}}:=\frac{1}{2}\left(\epsilon -\sigma \right);\end{array}\end{array}$$

project ρ _reg=i d⊕s g n onto the id and sgn representations of ${\mathbb{Z}}_2$, respectively.

This observation prompts us to work in the alternative basis:

$$\begin{array}{l}\begin{array}{cc}{f}_0:& ={\Theta }_{\mathit{\text{id}}}·{\xi }_0={\Theta }_{\mathit{\text{id}}}·{\xi }_1=h{\xi }_0={\xi }_0+{\xi }_1,\\ f_1:& ={\Theta }_{\mathit{\text{sgn}}}·{\xi }_0=-{\Theta }_{\mathit{\text{sgn}}}·{\xi }_1=h{\xi }_1={\xi }_0-{\xi }_1.\end{array}\end{array}$$

In this basis the permutation matrix is diagonal:

$$\begin{array}{l}\begin{array}{c}\hat{K}:=\mathit{\text{hK}}{h}^{-1}=\left(\begin{array}{ll}1& 0\\ 0& -1\end{array}\right),\hat{Q}:=-\mathbf{1}+\hat{K}=\left(\begin{array}{ll}0& 0\\ 0& -2\end{array}\right).\end{array}\end{array}$$

The representation-theoretic perspective on $\hat{K}$ is to observe that i d(σ)=1 and s g n(σ)=−1.

Referring to (5), we know that we can write a generic phylogenetic tensor as

$$\begin{array}{l}\begin{array}{c}P=e^{-\lambda }exp\left(\sum _{\varnothing \ne e\subseteq \left[n-1\right]}{\alpha }_e{K}^{\left(e\right)}\right)·{\delta }^{n-1}\pi ,\end{array}\end{array}$$

where $\lambda =\sum _{\varnothing \ne e\subseteq \left[n-1\right]}{\alpha }_e$.

We index matrix and tensor indices by using $i,j,k=0,1\in {\mathbb{Z}}_2$ and allow multiplication × in the ring of integers . The Hadamard matrix then has matrix elements ${\left[h\right]}_i^j={(-1)}^{i\times j}$ where j is the row index and i is the column index. Observe that in the diagonal basis, the permutation matrix has elements

$$\begin{array}{l}\begin{array}{c}{\left[\hat{K}\right]}_i^j={\delta }_{\mathit{\text{ij}}}{(-1)}^i.\end{array}\end{array}$$

Thus we have expressions such as

$$\begin{array}{l}\begin{array}{c}{\left[{\hat{K}}^{\left(\right\{2,3\left\}\right)}\right]}_{i_1{i}_2{i}_3}^{j_1{j}_2{j}_3}={\delta }_{i_1{j}_1}{\delta }_{i_2{j}_2}{\delta }_{i_3{j}_3}{(-1)}^{i_2+i_3},\end{array}\end{array}$$

where ${\hat{K}}^{\left(\right\{2,3\left\}\right)}=\mathbf{1}\otimes \hat{K}\otimes \hat{K}$.

As we are dealing with tensors of arbitrary size, it is convenient to represent a string such as i ₁i ₂…i _n as an ordered-bipartition μ=μ ₀: μ ₁ of the set [n], where μ ₀,μ ₁⊆[n] with j∈μ _k if and only if i _j =k. For example we have the following equivalences:

$$\begin{array}{l}\begin{array}{cc}00110& \equiv \{1,2,5\}:\{3,4\},01111\equiv \left\{1\right\}:\{2,3,4,5\},\\ 10001& \equiv \{2,3,4\}:\{1,5\}\end{array}\end{array}$$

and inequivalence:

$$\begin{array}{l}\begin{array}{c}01010\equiv \{1,3,5\}:\{2,4\}\ne \{2,4\}:\{1,3,5\}\equiv 10101.\end{array}\end{array}$$

We then have

$$\begin{array}{l}\begin{array}{cc}{\left[{\hat{K}}^{\left(e\right)}\right]}_{i_1{i}_2\dots i_n}^{j_1{j}_2\dots j_n}& ={\left[{\hat{K}}^{\left(e\right)}\right]}_{\mu }^{\nu }={\left[{\hat{K}}^{\left(e\right)}\right]}_{{\mu }_0:{\mu }_1}^{{\nu }_0:{\nu }_1}\\ ={\delta }_{{\mu }_0{\nu }_0}{\delta }_{{\mu }_1{\nu }_1}{(-1)}^{|e\cap {\mu }_1|}.\end{array}\end{array}$$

Defining h ⁽ⁿ⁾:=h ⁽ⁿ⁻¹⁾⊗h where h ⁽¹⁾:=h, in the diagonal basis $\hat{P}:=h^{\left(n\right)}·P$ and using our notation h ⁽ⁿ⁾ has tensor components

$$\begin{array}{l}\begin{array}{cc}{\left[h^{\left(n\right)}\right]}_{\mu }^{\nu }={\left[h^{\left(n\right)}\right]}_{{\mu }_0:{\mu }_1}^{{\nu }_0:{\nu }_1}& ={\left[h^{\left(n\right)}\right]}_{i_1{i}_2\dots i_n}^{j_1{j}_2\dots j_n}\\ ={(-1)}^{i_1\times j_1+i_2\times j_2+\dots +i_n\times j_n}\\ ={(-1)}^{|{\mu }_1\cap {\nu }_1|}.\end{array}\end{array}$$

The zero edge-length star-tree initial distribution has tensor components

$$\begin{array}{l}\begin{array}{c}{\left[{\delta }^{n-1}\pi \right]}_{i_1{i}_2\dots i_n}=\frac{1}{2}{\delta }_{i_1{i}_2}{\delta }_{i_1{i}_3}\dots {\delta }_{i_1{i}_n},\end{array}\end{array}$$

(where, although it seems we have given preference to taxon 1 in this expression, there are many ways that this distribution can be expressed using the δ _ij ). In the diagonal basis with $\widehat{{\delta }^{n-1}\pi }:=h^{\left(n\right)}·{\delta }^{n-1}\pi$, we have components

$$\begin{array}{l}\begin{array}{c}{\left[\widehat{{\delta }^{n-1}\pi }\right]}_{i_1{i}_2\dots i_n}\\ =\frac{1}{2}\sum _{j_1,j_2,\dots ,j_n}{(-1)}^{i_1\times j_1+i_2\times j_2+\dots +i_n\times j_n}{\delta }_{j_1{j}_2}{\delta }_{j_1{j}_3}\dots {\delta }_{j_1{j}_n}\\ =\frac{1}{2}\sum _{j_1}{(-1)}^{\left(i_1+i_2+\dots +i_n\right)\times j_1}=\frac{1}{2}\left(1+{(-1)}^{i_1+i_2+\dots +i_n}\right),\end{array}\end{array}$$

which is exactly the statement

$$\begin{array}{l}\begin{array}{c}{\left[\widehat{{\delta }^{n-1}\pi }\right]}_{\mu }={\left[\widehat{{\delta }^{n-1}\pi }\right]}_{{\mu }_0:{\mu }_1}=\frac{1}{2}\left(1+{(-1)}^{\left|{\mu }_1\right|}\right).\end{array}\end{array}$$

Since $\hat{K}$ is diagonal in the transformed basis, we can conclude that

$$\begin{array}{l}\begin{array}{cc}{\left[\hat{P}\right]}_{\mu }& ={\left[\hat{P}\right]}_{{\mu }_0:{\mu }_1}\\ =e^{-\lambda }exp\left(\sum _{\varnothing \ne e\subseteq \left[2,n\right]}{\alpha }_e\underset{{\mu }_0:{\mu }_1}{\overset{{\mu }_0:{\mu }_1}{\left[{\hat{K}}^{\left(e\right)}\right]}}\right)\frac{1}{2}\left(1+{(-1)}^{\left|{\mu }_1\right|}\right).\end{array}\end{array}$$

Of course many of these tensor components will be zero and we would like to ignore these.

Take u=u ₀: u ₁ as an ordered bipartition of the reduced set [n−1], so that u≡i ₁i ₂…i _n−1 where j∈u _k if and only if i _j =k, and define

$$\begin{array}{l}\begin{array}{cc}\gamma \left(u\right)& =\left\{\begin{array}{l}0,\text{if}\left|u_1\right|\text{is even,}\\ 1,\text{if}\left|u_1\right|\text{is odd;}\end{array}\right.\\ =2-\left(0\left|u_0\right|+1\left|u_1\right|\right)\text{(mod 2)},\end{array}\end{array}$$

and interpret u·γ(u) as a string: u·γ(u) = i ₁i ₂…i _n−1γ(u).

If we make the definitions

$$\begin{array}{l}\begin{array}{c}{\mathcal{P}}_u:={\left[\hat{P}\right]}_{u·\gamma \left(u\right)},{\eta }_u:=\frac{1}{2}\sum _{\varnothing \ne e\subseteq \left[n-1\right]}{\alpha }_e{\left[{\hat{K}}^{\left(e\right)}\right]}_{u·\gamma \left(u\right)}^{u·\gamma \left(u\right)},\end{array}\end{array}$$

then we can write the non-zero components as

$$\begin{array}{l}\begin{array}{c}{\mathcal{P}}_u=e^{-\lambda }exp\left({\eta }_u\right),\end{array}\end{array}$$

with inverses

$$\begin{array}{l}\begin{array}{c}{\eta }_u=ln\left({\mathcal{P}}_u\right)+\mathrm{\lambda .}\end{array}\end{array}$$

(6)

This is the first part of the inversion.

We would like to go further and actually recover the individual edge weights α _e . To do this we define the (square) 2ⁿ⁻¹×2ⁿ⁻¹ matrix F with components

$$\begin{array}{l}\begin{array}{c}{\left[F\right]}_u^e:={\left[{\hat{K}}^{\left(e\right)}\right]}_{u·\gamma \left(u\right)}^{u·\gamma \left(u\right)}={(-1)}^{|e\cap u|}={\left[h^{(n-1)}\right]}_u^e,\end{array}\end{array}$$

with e a subset and u an ordered-bipartition of [n−1]. As ${\left(h^{(n-1)}\right)}^2=\frac{1}{2^{n-1}}\mathbf{1}$, we see that F provides its own inverse F ⁻¹ with components

$$\begin{array}{l}\begin{array}{c}{\left[F^{-1}\right]}_e^u:=\frac{1}{2^{n-1}}{\left[F\right]}_u^e.\end{array}\end{array}$$

Defining the column vectors $\overrightarrow{\alpha }=\left\{{\alpha }_e\right\}$ and $\overrightarrow{\eta }=\left\{{\eta }_u\right\}$, we can write the matrix equations

$$\begin{array}{l}\begin{array}{c}\overrightarrow{\eta }=F\overrightarrow{\alpha },\overrightarrow{\alpha }=F^{-1}\overrightarrow{\eta }.\end{array}\end{array}$$

Together with the first part of the inversion (6), these equations give a one-one map between pattern probabilities and edge weights for the binary-symmetric model.

Inversion of the ${\mathbb{Z}}_3$model

Taking confidence from the previous case we now discuss the inversion of the group-based phylogenetic model with $G={\mathbb{Z}}_3$. We take

$${\mathbb{Z}}_3={\{0,1,2\}}_{+\text{(mod 3)}}\cong 〈\sigma |{\sigma }^3=\epsilon 〉,$$

and, by analogy to the ${\mathbb{Z}}_2$ case, index tensors with indices i,j=0,1,2 and allow multiplication × by extending ${\mathbb{Z}}_3$ to the ring ${\mathbf{F}}_3={\{0,1,2\}}_{+,\times \text{(mod 3)}}$.

In this case a generic rate matrix is given by

$$\begin{array}{l}\begin{array}{cc}Q& =\left(\begin{array}{lll}-(\alpha +\beta )& \beta & \alpha \\ \alpha & -(\alpha +\beta )& \beta \\ \beta & \alpha & -(\alpha +\beta )\end{array}\right)\\ =-\left(\alpha +\beta \right)\mathbf{1}+\alpha K_1+\beta K_2,\end{array}\end{array}$$

where

$$\begin{array}{l}\begin{array}{c}{K}_1=\left(\begin{array}{lll}0& 0& 1\\ 1& 0& 0\\ 0& 1& 0\end{array}\right),K_2=\left(\begin{array}{lll}0& 1& 0\\ 0& 0& 1\\ 1& 0& 0\end{array}\right),\end{array}\end{array}$$

are the matrices representing the permutations σ≅(123) and σ ²≅(132) under the regular representation, respectively.

We define ω=e ^2πi/3, and present the character table of ${\mathbb{Z}}_3$ in Table 2. The decomposition of the regular representation is ρ _reg=i d⊕ω⊕ω ², and the columns of the character table give the projection operators onto the (one-dimensional) irreducible subspaces:

$$\begin{array}{l}\begin{array}{cc}{\Theta }_{\mathit{\text{id}}}:& =\frac{1}{3}\left(\epsilon +\sigma +{\sigma }^2\right),\\ {\Theta }_{\omega }:& =\frac{1}{3}\left(\epsilon +\mathrm{\omega \sigma }+{\omega }^2{\sigma }^2\right),\\ {\Theta }_{{\omega }^2}:& =\frac{1}{3}\left(\epsilon +{\omega }^2\sigma +\omega {\sigma }^2\right).\end{array}\end{array}$$

Table 2 The character table of ${\mathbb{Z}}_3$

Therefore, the matrix

$$\begin{array}{l}\begin{array}{c}f=\left(\begin{array}{lll}1& 1& 1\\ 1& \omega & {\omega }^2\\ 1& {\omega }^2& \omega \end{array}\right),\end{array}\end{array}$$

diagonalizes the generic rate matrix for this model:

$$\begin{array}{l}\begin{array}{c}\hat{Q}=\mathit{\text{fQ}}{f}^{-1}=\left(\begin{array}{lll}0& 0& 0\\ 0& \mathrm{\alpha \omega }+\beta {\omega }^2& 0\\ 0& 0& \alpha {\omega }^2+\mathrm{\beta \omega }\end{array}\right),\end{array}\end{array}$$

or, equivalently,

$$\begin{array}{l}\begin{array}{cc}{\hat{K}}_1& ={\mathit{\text{fK}}}_1{f}^{-1}=\left(\begin{array}{lll}1& 0& 0\\ 0& \omega & 0\\ 0& 0& {\omega }^2\end{array}\right),\\ {\hat{K}}_2& ={\mathit{\text{fK}}}_2{f}^{-1}=\left(\begin{array}{lll}1& 0& 0\\ 0& {\omega }^2& 0\\ 0& 0& \omega \end{array}\right).\end{array}\end{array}$$

We recall our basic result (5) that for group-based models, a generic phylogenetic tensor can be expressed as

$$\begin{array}{l}\begin{array}{c}P=e^{-\lambda }exp\left(\sum _{\varnothing \ne e\subseteq \left[n-1\right]}\left({\alpha }_e\underset{1}{\overset{\left(e\right)}{K}}+{\beta }_e\underset{2}{\overset{\left(e\right)}{K}}\right)\right)·{\delta }^{n-1}\pi ,\end{array}\end{array}$$

where $\lambda =\sum _{\varnothing \ne e\subseteq \left[n-1\right]}\left({\alpha }_e+{\beta }_e\right)$. We take the stationary distribution as initial distribution, so $\pi ={\left(\frac{1}{3},\frac{1}{3},\frac{1}{3}\right)}^T$.

The matrix elements of f can be expressed as ${\left[f\right]}_i^j={\omega }^{i\times j}$, where we extend $i,j\in {\mathbb{Z}}_3$ to include multiplication × from the ring of integers . Similarly,

$$\begin{array}{l}\begin{array}{c}{\left[{\hat{K}}_1\right]}_i^j={\delta }_{\mathit{\text{ij}}}{\omega }^i,{\left[{\hat{K}}_2\right]}_i^j={\delta }_{\mathit{\text{ij}}}{\left({\omega }^2\right)}^i.\end{array}\end{array}$$

More generally, tensorial components can be expressed as

$$\begin{array}{l}\begin{array}{c}{\left[\mathbf{1}\otimes \widehat{K_1}\otimes \widehat{K_1}\right]}_{i_1{i}_2{i}_3}^{j_1{j}_2{j}_3}={\delta }_{i_1{j}_1}{\delta }_{i_2{j}_2}{\delta }_{i_3{j}_3}{\omega }^{i_2+i_3}.\end{array}\end{array}$$

We represent a string i ₁i ₂…i _n as an ordered-tripartition, i ₁i ₂…i _n ≡μ=μ ₀: μ ₁: μ ₂, of the set [n], where j∈μ _k if and only if i _j =k. For example, if we take n=5, we have:

$$\begin{array}{l}\begin{array}{cc}00000& \equiv \{1,2,3,4,5\}\mathrm{:\varnothing :\varnothing },20120\equiv \{2,5\}:\left\{3\right\}:\{1,4\},\\ 01122& \equiv \left\{1\right\}:\{2,3\}:\{4,5\}.\end{array}\end{array}$$

Taking n = 3, we have

$$\begin{array}{l}\begin{array}{cc}{\left[{\hat{K}}_1^{\left(\right\{2,3\left\}\right)}\right]}_{\mu }^{\nu }={\left[\mathbf{1}\otimes {\hat{K}}_1\otimes {\hat{K}}_1\right]}_{\mu }^{\nu }& ={\left[\mathbf{1}\otimes {\hat{K}}_1\otimes {\hat{K}}_1\right]}_{{\mu }_0:{\mu }_1:{\mu }_2}^{{\nu }_0:{\nu }_1:{\nu }_2}\\ ={\delta }_{\mathrm{\mu \nu }}{\omega }^{|{\mu }_1\cap \{2,3\left\}\right|+2|{\mu }_2\cap \{2,3\left\}\right|},\end{array}\end{array}$$

and in general:

$$\begin{array}{l}\begin{array}{cc}{\left[{\hat{K}}_1^{\left(e\right)}\right]}_{\mu }^{\nu }& ={\delta }_{\mathrm{\mu \nu }}{\omega }^{|e\cap {\mu }_1|+2|e\cap {\mu }_2|},\\ {\left[{\hat{K}}_2^{\left(e\right)}\right]}_{\mu }^{\nu }& ={\delta }_{\mathrm{\mu \nu }}{\omega }^{|e\cap {\mu }_2|+2|e\cap {\mu }_1|}.\end{array}\end{array}$$

Taking the uniform distribution as initial distribution, the initial star-tree distribution can be written as

$$\begin{array}{l}\begin{array}{c}{\left[{\delta }^{n-1}\pi \right]}_{i_1{i}_2\dots i_n}=\frac{1}{3}{\delta }_{i_1{i}_2}{\delta }_{i_1{i}_3}\dots {\delta }_{i_1{i}_n}.\end{array}\end{array}$$

Defining f ⁽ⁿ⁾=f ⁽ⁿ⁻¹⁾⊗f where f ⁽¹⁾=f, we have

$$\begin{array}{l}\begin{array}{cc}{\left[f^{\left(n\right)}\right]}_{\mu }^{\nu }& ={\left[f^{\left(n\right)}\right]}_{i_1{i}_2\dots i_n}^{j_1{j}_2\dots j_n}={\left[f\right]}_{i_1}^{j_1}{\left[f\right]}_{i_2}^{j_2}\dots {\left[f\right]}_{i_n}^{j_n}\\ ={\omega }^{i_1\times j_1+i_2\times j_2+\dots +i_n\times j_n},\end{array}\end{array}$$

and in the transformed basis, where $\widehat{{\delta }^{n-1}\pi }:=f^{\left(n\right)}·{\delta }^{n-1}\pi$, we have

$$\begin{array}{l}\begin{array}{cc}{\left[\widehat{{\delta }^{n-1}\pi }\right]}_{i_1{i}_2\dots i_n}& =\frac{1}{3}\sum _{j_1,j_2,\dots ,j_n}{\omega }^{i_1\times j_1+i_2\times j_2+\dots +i_n\times j_n}\\ \times {\delta }_{j_1{j}_2}{\delta }_{j_1{j}_3}\dots {\delta }_{j_1{j}_n}\\ =\frac{1}{3}\sum _{j_1}{\omega }^{j_1\times (i_1+i_2+\dots +i_n)}\\ =\frac{1}{3}\left(1+{\omega }^{i_1+i_2+\dots +i_n}+{\left({\omega }^2\right)}^{i_1+i_2+\dots +i_n}\right).\end{array}\end{array}$$

Indexing by ordered-tripartitions, we conclude that

$$\begin{array}{l}\begin{array}{cc}{\left[\widehat{{\delta }^{n-1}\pi }\right]}_{\mu }& =\frac{1}{3}\left(1+{\omega }^{i_1+i_2+\dots +i_n}+{\left({\omega }^2\right)}^{i_1+i_2+\dots +i_n}\right)\\ =\frac{1}{3}\left(1+{\omega }^{\left|{\mu }_1\right|+2\left|{\mu }_2\right|}+{\left({\omega }^2\right)}^{\left|{\mu }_1\right|+2\left|{\mu }_2\right|}\right).\end{array}\end{array}$$

Now suppose |μ ₁|+2|μ ₂|=0 (mod 3), then

$$\begin{array}{l}\begin{array}{c}{\left[\widehat{{\delta }^{n-1}\pi }\right]}_{\mu }=\frac{1}{3}\left(1+1+1\right)=1.\end{array}\end{array}$$

If |μ ₁|+2|μ ₂|=1 (mod 3), then

$$\begin{array}{l}\begin{array}{c}{\left[\widehat{{\delta }^{n-1}\pi }\right]}_{\mu }=\frac{1}{3}\left(1+\omega +{\omega }^2\right)=0,\end{array}\end{array}$$

and if |μ ₁|+2|μ ₂|=2 (mod 3), then

$$\begin{array}{l}\begin{array}{c}{\left[\widehat{{\delta }^{n-1}\pi }\right]}_{\mu }=\frac{1}{3}\left(1+{\omega }^2+\omega \right)=0.\end{array}\end{array}$$

Thus we have found a basis where all the elements of the initial star-tree tensor are zero unless the tripartion μ satisfies |μ ₁|+2|μ ₂|=0 (mod 3). Crucially, this statement also holds for the phylogenetic tensor $\hat{P}$ because in this basis the rate matrices of this model are diagonal:

$$\begin{array}{l}\begin{array}{cc}{\left[\hat{P}\right]}_{\mu }& ={\left[\hat{P}\right]}_{{\mu }_0:{\mu }_1:{\mu }_2}\\ =e^{-\lambda }exp\left(\frac{1}{2}\sum _{\varnothing \ne e\subseteq \left[n-1\right]}\underset{{\mu }_0:{\mu }_1:{\mu }_2}{\overset{{\mu }_0:{\mu }_1:{\mu }_2}{\left[{\alpha }_e\underset{1}{\overset{\left(e\right)}{K}}+{\beta }_e\underset{2}{\overset{\left(e\right)}{K}}\right]}}\right)\\ \times \frac{1}{3}\left(1+{\omega }^{1\left|{\mu }_1\right|}+{\omega }^{2\left|{\mu }_2\right|}\right).\end{array}\end{array}$$

We deal with this condition on μ by taking u=u ₀: u ₁: u ₂ as an ordered-tripartion of the reduced set [n−1] and setting μ=u·γ(u) (considered as the concatenation of strings) where

$$\begin{array}{l}\begin{array}{cc}\gamma \left(u\right)& =\left\{\begin{array}{ll}0,& \text{if}\left|u_1\right|+2\left|u_2\right|=0\\ 1,& \text{if}\left|u_1\right|+2\left|u_2\right|=2\\ 2;& \text{if}\left|u_1\right|+2\left|u_2\right|=1\end{array}\right.\\ =3-\left(0\left|u_0\right|+1\left|u_1\right|+2\left|u_2\right|\right)\text{(mod 3)}.\end{array}\end{array}$$

If we make the definitions

$$\begin{array}{l}\begin{array}{cc}{\mathcal{P}}_u& :={\left[\hat{P}\right]}_{u·\gamma \left(u\right)},\\ {\eta }_u& :={\left[\sum _{\varnothing \ne e\subseteq \left[n-1\right]}{\alpha }_e{K}_1^{\left(e\right)}+{\beta }_e{K}_2^{\left(e\right)}\right]}_{u·\gamma \left(u\right)}^{u·\gamma \left(u\right)},\end{array}\end{array}$$

we then have the first part of the inversion

$$\begin{array}{l}\begin{array}{c}{\mathcal{P}}_u=e^{-\lambda }exp\left({\eta }_u\right),{\eta }_u=ln\left({\mathcal{P}}_u\right)+\mathrm{\lambda .}\end{array}\end{array}$$

(7)

As in the ${\mathbb{Z}}_2$ case, we would like to use η _u to recover the rate parameters α _e ,β _e for all ∅≠e⊆[n−1] and thus complete the full inversion for this model. Of course, it is little bit more difficult this time.

Recall that μ=μ ₀: μ ₁: μ ₂ with μ _i ⊆[n], whereas u=u ₀: u ₁: u ₂ with u _i ⊆[n−1], and ∅≠e⊆[n−1]. Considering

$$\begin{array}{l}\begin{array}{cc}{\left[K_1^{\left(e\right)}\right]}_{\mu }^{\mu }& ={\omega }^{|e\cap {\mu }_1|+2|e\cap {\mu }_2|},\end{array}\end{array}$$

it follows that

$$\begin{array}{l}\begin{array}{cc}{\left[K_1^{\left(e\right)}\right]}_{u·\gamma \left(u\right)}^{u·\gamma \left(u\right)}& ={\omega }^{|e\cap u_1|+2|e\cap u_2|},\end{array}\end{array}$$

and similarly

$$\begin{array}{l}\begin{array}{c}{\left[K_2^{\left(e\right)}\right]}_{u·\gamma \left(u\right)}^{u·\gamma \left(u\right)}={\omega }^{|e\cap u_2|+2|e\cap u_1|}.\end{array}\end{array}$$

We make the observation that

$$\begin{array}{l}\begin{array}{c}{\left[F_1\right]}_u^e:={\left[f^{(n-1)}\right]}_{u_0:u_1:u_2}^{e^c\mathrm{:e:\varnothing }}={\omega }^{|u_1\cap e|+2|u_2\cap e|}={\left[K_{\alpha }^{\left(e\right)}\right]}_{u·\gamma \left(u\right)}^{u·\gamma \left(u\right)},\end{array}\end{array}$$

and

$$\begin{array}{l}\begin{array}{c}{\left[F_2\right]}_u^e:={\left[f^{(n-1)}\right]}_{u_0:u_1:u_2}^{e^c\mathrm{:\varnothing :e}}={\omega }^{|u_2\cap e|+2|u_1\cap e|}={\left[K_{\beta }^{\left(e\right)}\right]}_{u·\gamma \left(u\right)}^{u·\gamma \left(u\right)},\end{array}\end{array}$$

where F ₁ and F ₂ are 2ⁿ⁻¹×3ⁿ⁻¹ matrices.

Thus we may write

$$\begin{array}{l}\begin{array}{c}{\eta }_u=\sum _{\varnothing \ne e\subseteq \left[n-1\right]}{\alpha }_e{\left[F_1\right]}_u^e+{\beta }_e{\left[F_2\right]}_u^e.\end{array}\end{array}$$

Defining the column vectors $\overrightarrow{\alpha }=\left\{{\alpha }_e\right\},\overrightarrow{\beta }=\left\{{\beta }_e\right\}$ and $\overrightarrow{\eta }=\left\{{\eta }_u\right\}$, we can write

$$\begin{array}{l}\begin{array}{c}\overrightarrow{\eta }=F_1\overrightarrow{\alpha }+F_2\overrightarrow{\beta },\end{array}\end{array}$$

and define two 3ⁿ⁻¹×2ⁿ⁻¹ matrices G ₁ and G ₂ as

$$\begin{array}{l}\begin{array}{c}{\left[G_1\right]}_e^u:={\left[{f^{-1}}^{(n-1)}\right]}_{e^c\mathrm{:e:\varnothing }}^u,{\left[G_2\right]}_e^u:={\left[{f^{-1}}^{(n-1)}\right]}_{e^c\mathrm{:\varnothing :e}}^u,\end{array}\end{array}$$

where

$$\begin{array}{l}\begin{array}{c}{f}^{-1}=\left(\begin{array}{lll}1& 1& 1\\ 1& \omega & {\omega }^2\\ 1& {\omega }^2& \omega \end{array}\right),\end{array}\end{array}$$

with f f ⁻¹=1.

Considering that

$$\begin{array}{l}\begin{array}{c}\sum _v{\left[{f^{-1}}^{(n-1)}\right]}_u^v{\left[f^{(n-1)}\right]}_v^w={\delta }_{\mathit{\text{uw}}},\end{array}\end{array}$$

for all ordered-triparitions u,w of [n−1], we have the matrix products

$$\begin{array}{l}\begin{array}{c}\begin{array}{ll}{G}_1{F}_1=\mathbf{1},& G_1{F}_2=0,\\ G_2{F}_2=\mathbf{1},& G_2{F}_1=0.\end{array}\end{array}\end{array}$$

Thus the second part of the inversion for this model is

$$\begin{array}{l}\begin{array}{c}\overrightarrow{\alpha }=G_1\overrightarrow{\eta },\overrightarrow{\beta }=G_2\overrightarrow{\eta }.\end{array}\end{array}$$

Together with (7), these equations give a one-one map between pattern probabilities and edge weights for the group-based model with $G={\mathbb{Z}}_3$.

Inversion of the K3ST model

We now consider the K3ST model [36] which occurs as the group-based model with

$$\begin{array}{cc}G& ={\mathbb{Z}}_2\times {\mathbb{Z}}_2={\left\{\right(0,0),(0,1),(1,0),(1,1\left)\right\}}_{+\text{(mod 2)}}\\ \cong 〈\left(12\right)\left(34\right),\left(13\right)\left(24\right)〉.\end{array}$$

In this model a generic rate matrix is given by

$$\begin{array}{l}\begin{array}{c}Q=-\left(\alpha +\beta +\gamma \right)\mathbf{1}+\alpha K_{01}+\beta K_{10}+\gamma K_{11},\end{array}\end{array}$$

where

$$\begin{array}{l}\begin{array}{cc}{K}_{01}& =\mathbf{1}\otimes K=\left(\begin{array}{llll}0& 1& 0& 0\\ 1& 0& 0& 0\\ 0& 0& 0& 1\\ 0& 0& 1& 0\end{array}\right),\\ K_{10}& =K\otimes \mathbf{1}=\left(\begin{array}{llll}0& 0& 1& 0\\ 0& 0& 0& 1\\ 1& 0& 0& 0\\ 0& 1& 0& 0\end{array}\right),\\ K_{11}& =K\otimes K=\left(\begin{array}{llll}0& 0& 0& 1\\ 0& 0& 1& 0\\ 0& 1& 0& 0\\ 1& 0& 0& 0\end{array}\right).\end{array}\end{array}$$

(8)

We already know that the 2×2 Hadamard matrix h diagonalizes K, so we see immediately that H=h⊗h diagonalizes this model:

$$\begin{array}{l}\begin{array}{cc}{\hat{K}}_{01}:& ={\mathit{\text{HK}}}_{01}{H}^{-1}=\mathbf{1}\otimes \mathit{\text{hK}}{h}^{-1}=\left(\begin{array}{llll}1& 0& 0& 0\\ 0& -1& 0& 0\\ 0& 0& 1& 0\\ 0& 0& 0& -1\end{array}\right),\\ {\hat{K}}_{10}:& ={\mathit{\text{HK}}}_{10}{H}^{-1}=\left(\begin{array}{llll}1& 0& 0& 0\\ 0& 1& 0& 0\\ 0& 0& -1& 0\\ 0& 0& 0& -1\end{array}\right),\\ {\hat{K}}_{11}:& ={\mathit{\text{HK}}}_{11}{H}^{-1}=\left(\begin{array}{llll}1& 0& 0& 0\\ 0& -1& 0& 0\\ 0& 0& -1& 0\\ 0& 0& 0& 1\end{array}\right).\end{array}\end{array}$$

Of course H is the character table of ${\mathbb{Z}}_2\times {\mathbb{Z}}_2$ and the permutation matrices (8), together with K ₀₀:=1, give the regular representation ρ _reg≅i d⊗i d⊕i d⊗s g n⊕s g n⊗i d⊕s g n⊗s g n, where we recall the basic result that the tensor product of two irreducible representations of a group G gives an irreducible representation of G×G.

Simplifying notation, for this model we index tensors with indices given as pairs: $i,j=00,01,10,11\in {\mathbb{Z}}_2\times {\mathbb{Z}}_2$; and we express the individual parts using lower case Roman characters. For example, we write i:=a b=01, with a=0 and b=1. This gives matrix elements:

$$\begin{array}{l}\begin{array}{cc}{\left[{\hat{K}}_{01}\right]}_{\mathit{\text{ab}}}^{\mathit{\text{cd}}}& ={\delta }_{\mathit{\text{ac}}}{\delta }_{\mathit{\text{bd}}}{(-1)}^b,{\left[{\hat{K}}_{10}\right]}_{\mathit{\text{ab}}}^{\mathit{\text{cd}}}={\delta }_{\mathit{\text{ac}}}{\delta }_{\mathit{\text{bd}}}{(-1)}^a,\\ {\left[{\hat{K}}_{11}\right]}_{\mathit{\text{ab}}}^{\mathit{\text{cd}}}& ={\delta }_{\mathit{\text{ac}}}{\delta }_{\mathit{\text{bd}}}{(-1)}^{a+b};\end{array}\end{array}$$

and more complicated tensor products such as

$$\begin{array}{l}\begin{array}{c}{\left[{\hat{K}}_{01}\otimes {\hat{K}}_{01}\otimes \mathbf{1}\right]}_{a_1{b}_1{a}_2{b}_2{a}_3{b}_3}^{c_1{d}_1{c}_2{d}_2{c}_3{d}_3}\\ ={\delta }_{a_1{c}_1}{\delta }_{b_1{d}_1}{\delta }_{a_2{c}_2}{\delta }_{b_2{d}_2}{\delta }_{a_3{c}_3}{\delta }_{b_3{d}_3}{(-1)}^{b_1+b_2}.\end{array}\end{array}$$

Again we interpret strings such as μ≡a ₁a ₂…a _n and ν≡b ₁b ₂…b _n as ordered-bipartitions μ=μ ₀: μ ₁ and ν=ν ₀: ν ₁ of the set [n]. We can then write matrix elements of tensor products as

$$\begin{array}{l}\begin{array}{cc}{\left[{\hat{K}}_{01}^{\left(e\right)}\right]}_{\mu ,\nu }^{{\mu }^{\prime },{\nu }^{\prime }}& ={\delta }_{\mu {\mu }^{\prime }}{\delta }_{\nu {\nu }^{\prime }}{(-1)}^{|e\cap {\nu }_1|},\\ {\left[{\hat{K}}_{10}^{\left(e\right)}\right]}_{\mu ,\nu }^{{\mu }^{\prime },{\nu }^{\prime }}& ={\delta }_{\mu {\mu }^{\prime }}{\delta }_{\nu {\nu }^{\prime }}{(-1)}^{|e\cap {\mu }_1|},\\ {\left[{\hat{K}}_{11}^{\left(e\right)}\right]}_{\mu ,\nu }^{{\mu }^{\prime },{\nu }^{\prime }}& ={\delta }_{\mu {\mu }^{\prime }}{\delta }_{\nu {\nu }^{\prime }}{(-1)}^{|e\cap {\mu }_1|+|e\cap {\nu }_1|}.\end{array}\end{array}$$

Taking the stationary distribution $\pi =\frac{1}{4}{(1,1,1,1)}^T$ as initial distribution, the zero edge-length star-tree distribution is given by

$$\begin{array}{l}\begin{array}{c}{\left[{\delta }^{n-1}\pi \right]}_{i_1{i}_2\dots i_n}=\frac{1}{4}{\delta }_{i_1{i}_2}{\delta }_{i_1{i}_3}\dots {\delta }_{i_1{i}_n},\end{array}\end{array}$$

which in the finer index representation is

$$\begin{array}{l}\begin{array}{c}{\left[{\delta }^{n-1}\pi \right]}_{a_1{b}_1{a}_2{b}_2\dots a_n{b}_n}\\ =\frac{1}{4}{\delta }_{a_1{a}_2}{\delta }_{a_1{a}_3}\dots {\delta }_{a_1{a}_n}{\delta }_{b_1{b}_2}{\delta }_{b_1{b}_3}\dots {\delta }_{b_1{b}_n}.\end{array}\end{array}$$

Recall that elements of the Hadamard matrix can be written as ${\left[h\right]}_b^a={(-1)}^{a\times b}$, where $a,b\in {\mathbb{Z}}_2$ and we allow multiplication × by extending to the ring of integers . In the transformed basis, we have

$$\begin{array}{l}\begin{array}{cc}{\left[\widehat{{\delta }^{n-1}\pi }\right]}_{a_1{b}_1{a}_2{b}_2\dots a_n{b}_n}& ={\left[\widehat{{\delta }^{n-1}\pi }\right]}_{\mu ,\nu }\\ =\frac{1}{4}\sum _{c_1,c_2,\dots ,c_n}^{d_1,d_2,\dots ,d_n}{\left[h\right]}_{c_1}^{a_1}{\left[h\right]}_{c_2}^{a_2}\dots {\left[h\right]}_{c_n}^{a_n}{\left[h\right]}_{d_1}^{b_1}{\left[h\right]}_{d_2}^{b_2}\dots \\ \times {\left[h\right]}_{d_n}^{b_n}{\delta }_{a_1{a}_2}{\delta }_{a_1{a}_3}\dots {\delta }_{a_1{a}_n}{\delta }_{b_1{b}_2}{\delta }_{b_1{b}_3}\dots {\delta }_{b_1{b}_n}\\ =\frac{1}{4}\sum _{c_1,d_1}{(-1)}^{(a_1+a_2+\dots a_n)\times c_1+(b_1+b_2+\dots +b_n)\times d_1}\\ =\frac{1}{4}\left(1+{(-1)}^{a_1+\dots +a_n}+{(-1)}^{b_1+\dots +b_n}\right.\\ \left.+{(-1)}^{a_1+\dots +a_n+b_1+\dots +b_n}\right)\\ =\left\{\begin{array}{ll}0,& \text{if either}\left|{\mu }_1\right|\text{or}\left|{\nu }_1\right|\text{is odd};\\ 1,& \text{if}\left|{\mu }_1\right|\text{and}\left|{\nu }_1\right|\text{are both even}.\end{array}\right.\end{array}\end{array}$$

We recall (5), so under this model we can express a generic phylogenetic tensor as

$$\begin{array}{l}\begin{array}{c}P=e^{-\lambda }exp\left(\sum _{\varnothing \ne e\subseteq \left[n-1\right]}{\alpha }_e\underset{01}{\overset{\left(e\right)}{K}}+{\beta }_e\underset{10}{\overset{\left(e\right)}{K}}+{\gamma }_e\underset{11}{\overset{\left(e\right)}{K}}\right)·{\delta }^{n-1}\mathrm{\pi .}\end{array}\end{array}$$

To exclude the vanishing components we define, for all ordered bipartitions u=u ₀: u ₁ of the reduced set [n−1],

$$\begin{array}{l}\begin{array}{cc}\gamma \left(u\right)& =\left\{\begin{array}{l}0,\text{if}\left|u_1\right|\text{is even},\\ 1,\text{if}\left|u_1\right|\text{is odd};\end{array}\right.\\ =2-\left(0\right|u_0|+1|u_1\left|\right)\text{(mod 2)},\end{array}\end{array}$$

and intepret u·γ(u) as the string u·γ(u)=a ₁a ₂…a _n−1γ(u). Then, for each pair u,v of ordered-bipartitions of [n−1], we define

$$\begin{array}{l}\begin{array}{c}{\eta }_{u,v}:={\left[\sum _{\varnothing \ne e\subseteq \left[n-1\right]}{\alpha }_e{K}_{01}^{\left(e\right)}+{\beta }_e{K}_{10}^{\left(e\right)}+{\gamma }_e{K}_{11}^{\left(e\right)}\right]}_{u·\gamma \left(u\right),v·\gamma \left(v\right)}^{u·\gamma \left(u\right),v·\gamma \left(v\right)},\end{array}\end{array}$$

and

$$\begin{array}{l}\begin{array}{c}{\mathcal{P}}_{u,v}:={\left[P\right]}_{u·\gamma \left(u\right),v·\gamma \left(v\right)},\end{array}\end{array}$$

This gives the inversion

$$\begin{array}{l}\begin{array}{c}{\mathcal{P}}_{u,v}=e^{-\lambda }exp\left({\eta }_{u,v}\right),{\eta }_{u,v}=\lambda +ln\left({\mathcal{P}}_{u,v}\right).\end{array}\end{array}$$

Consider the 2ⁿ×2ⁿ⁻¹ rectangular matrices F ₀₁, F ₁₀ and F ₁₁ with components

$$\begin{array}{l}\begin{array}{c}\begin{array}{l}{\left[F_{01}\right]}_{u,v}^e={\left[K_{01}^{\left(e\right)}\right]}_{u,v}^{u,v}={(-1)}^{|e\cap v_1|},\\ {\left[F_{10}\right]}_{u,v}^e={\left[K_{10}^{\left(e\right)}\right]}_{u,v}^{u,v}={(-1)}^{|e\cap u_1|},\\ {\left[F_{11}\right]}_{u,v}^e={\left[K_{11}^{\left(e\right)}\right]}_{u,v}^{u,v}={(-1)}^{|e\cap u_1|+|e\cap v_1|};\end{array}\end{array}\end{array}$$

where e⊆[n−1] and u=u ₀: u ₁ and v=v ₀: v ₁ are ordered-bipartitions of [n−1]. If we define the column vector $\overrightarrow{\eta }:=\left\{{\eta }_{u,v}\right\}$ indexed by pairs of ordered-bipartitions and the column vectors $\overrightarrow{\alpha }:=\left\{{\alpha }_e\right\}$, $\overrightarrow{\beta }:=\left\{{\alpha }_e\right\}$ and $\overrightarrow{\gamma }:=\left\{{\alpha }_e\right\}$ indexed by subsets of [n−1], we then have the matrix equation

$$\begin{array}{l}\begin{array}{c}\overrightarrow{\eta }=F_{01}\overrightarrow{\alpha }+F_{10}\overrightarrow{\beta }+F_{11}\overrightarrow{\gamma }.\end{array}\end{array}$$

Writing H ⁽ⁿ⁾=H ⁽ⁿ⁻¹⁾⊗H with H ⁽¹⁾=H, we note that

$$\begin{array}{l}\begin{array}{cc}{\left[F_{01}\right]}_{u,v}^e& ={\left[H^{(n-1)}\right]}_{u,v}^{\varnothing ,e},{\left[F_{10}\right]}_{u,v}^e={\left[H^{(n-1)}\right]}_{u,v}^{e,\varnothing },\\ {\left[F_{11}\right]}_{u,v}^e& ={\left[H^{(n-1)}\right]}_{u,v}^{e,e};\end{array}\end{array}$$

and define the 2ⁿ⁻¹×2ⁿ rectangular matrices G ₀₁,G ₁₀ and G ₁₁ as

$$\begin{array}{l}\begin{array}{cc}{\left[G_{01}\right]}_e^{u,v}& ={\left[{H^{-1}}^{(n-1)}\right]}_{\varnothing ,e}^{u,v},{\left[G_{10}\right]}_e^{u,v}={\left[{H^{-1}}^{(n-1)}\right]}_{e,\varnothing }^{u,v},\\ {\left[G_{11}\right]}_e^{u,v}& ={\left[{H^{-1}}^{(n-1)}\right]}_{e,e}^{u,v}.\end{array}\end{array}$$

Noting that

$$\begin{array}{l}\begin{array}{c}\sum _{w,x}{\left[{H^{-1}}^{(n-1)}\right]}_{u,v}^{w,x}{\left[H^{(n-1)}\right]}_{w,x}^{y,z}={\delta }_{u,y}{\delta }_{v,z},\end{array}\end{array}$$

for all u,v,y,z ordered-bipartitions of [n−1], we then have the matrix identities

$$\begin{array}{l}\begin{array}{c}{G}_{01}{F}_{01}=\mathbf{1},G_{10}{F}_{10}=\mathbf{1},G_{11}{F}_{11}=\mathbf{1},\end{array}\end{array}$$

and

$$\begin{array}{l}\begin{array}{cc}{G}_{01}{F}_{10}& =0=G_{01}{F}_{11}=G_{\beta }{F}_{11}=G_{10}{F}_{01}\\ =G_{11}{F}_{01}=G_{11}{F}_{10}.\end{array}\end{array}$$

Writing

$$\begin{array}{l}\begin{array}{c}\overrightarrow{\alpha }=G_{01}\overrightarrow{\eta },\overrightarrow{\beta }=G_{10}\overrightarrow{\eta },\overrightarrow{\gamma }=G_{11}\overrightarrow{\eta },\end{array}\end{array}$$

completes the inversion for the K3ST model.

Inversion of the ${\mathbb{Z}}_r$model

We now consider the group based model for ${\mathbb{Z}}_r={\left\{0,1,2,\dots (r-1)\right\}}_{+\left(\text{mod r}\right)}\cong 〈\sigma :{\sigma }^r=e〉$. For this model the generic rate matrix has the form

$$\begin{array}{l}\begin{array}{c}Q=-\lambda \mathbf{1}+\sum _{i=1}^r{\alpha }^i{K}_{{\sigma }^i},\end{array}\end{array}$$

where $\lambda =\sum _{i=1}^r{\alpha }^i$ and

$$\begin{array}{l}\begin{array}{c}{K}_{\sigma }=\left(\begin{array}{lllll}0& 0& \dots & 0& 1\\ 1& 0& 0& 0& 0\\ 0& 1& \dots & 0& 0\\ ⋮& ⋮& \ddots & ⋮& ⋮\\ 0& 0& \dots & 1& 0\end{array}\right),\end{array}\end{array}$$

so that $K_{{\sigma }^i}=K_{\sigma }^i$.

Defining ω=e ^2πi/r, we have ω ^r=1 and 1+ω+ω ²+…+ω ^r−1=0 and ${\left[f\right]}_i^j={\omega }^{\mathit{\text{ij}}}$ where i,j=0,1,2,…,r−1. Of course, f is the character table of ${\mathbb{Z}}_r$ and ${\left[f^{-1}\right]}_j^i=\frac{1}{r}{\omega }^{-\mathit{\text{ij}}}$.

Lemma 1.

$$\begin{array}{l}\begin{array}{c}\sum _{\nu }{\left[f\otimes f\otimes \dots \otimes f\right]}_{\mu }^{\nu }{\left[f^{-1}\otimes f^{-1}\otimes \dots \otimes f^{-1}\right]}_{\nu }^{{\mu }^{\prime }}={\delta }_{\mu {\mu }^{\prime }},\end{array}\end{array}$$

where μ,ν,μ ^′ are ordered-r-partitions of the set [n] defined by the strings i ₁i ₂…i _n , j ₁j ₂…j _n and k ₁k ₂…k _n , respectively.

Proof.

The result is obvious by the definition of tensor product. However, explicitly we have

$$\begin{array}{l}\begin{array}{cc}\sum _{\nu }{\left[f\otimes f\otimes \dots \otimes f\right]}_{\mu }^{\nu }& {\left[f^{-1}\otimes f^{-1}\otimes \dots \otimes f^{-1}\right]}_{\nu }^{{\mu }^{\prime }}\\ =\frac{1}{r^n}\sum _{0\le j_1,j_2,\dots ,j_{r-1}\le (r-1)}\\ \times {\omega }^{i_1{j}_1+i_2{j}_2+\dots i_{r-1}{j}_{r-1}}{\omega }^{-\left(j_1{k}_1+j_2{k}_2+\dots +j_n{k}_n\right)}\\ =\frac{1}{r^n}\sum _{0\le j_1,j_2,\dots ,j_{r-1}\le (r-1)}\\ \times {\omega }^{j_1(i_1-k_1)+j_2(i_2-k_2)+\dots +j_n(i_n-k_n)}\end{array}\end{array}$$

which clearly equals 1 if i _ℓ −k _ℓ =0 for all ℓ, and, by repeatedly applying 1+ω+ω ²+…+ω ^r−1=0, equals 0 otherwise.

The regular representation contains exactly one copy of every irreducible representation and the irreducible representations of ${\mathbb{Z}}_r$ are given by the powers of ω:

$${\rho }_i:\begin{array}{l}{\mathbb{Z}}_r\to ℂ\\ \sigma \mapsto {\omega }^i\end{array}.$$

Thus the change of basis $K_{{\sigma }^i}\mapsto {\hat{K}}_{{\sigma }^i}={\mathit{\text{fK}}}_{{\sigma }^i}{f}^{-1}$ will give diagonal matrices ${\hat{K}}_{{\sigma }^i}$. Additionally,

Lemma 2.

In the diagonal basis, the matrices ${\hat{K}}_{{\sigma }^i}:={\mathit{\text{fK}}}_{{\sigma }^i}{f}^{-1}$ have matrix elements given by ${\left[{\hat{K}}_{{\sigma }^s}\right]}_i^j={\omega }^{\mathit{\text{is}}}{\delta }_{\mathit{\text{ij}}}$.

Proof.

Consider the matrix elements ${\left[K_{{\sigma }^s}\right]}_i^j={\delta }_{i{\sigma }^s\left(j\right)}$. Thus

$$\begin{array}{l}\begin{array}{cc}{\left[{\mathit{\text{fK}}}_{{\sigma }^s}{f}^{-1}\right]}_j^i& =\sum _{k,l}{\omega }^{\mathit{\text{ik}}}{\delta }_{k{\sigma }^s\left(l\right)}{\omega }^{-\mathit{\text{lj}}}=\sum _l{\omega }^{i{\sigma }^s\left(l\right)-\mathit{\text{lj}}}\\ =\sum _l{\omega }^{i(l+s)-\mathit{\text{lj}}}\\ ={\omega }^{\mathit{\text{is}}}\sum _l{\omega }^{l(i-j)}={\omega }^{\mathit{\text{is}}}{\delta }_{\mathit{\text{ij}}},\end{array}\end{array}$$

where we have used ${\omega }^{{\sigma }^s\left(m\right)}={\omega }^{m+s}$.

Now

$$\begin{array}{l}\begin{array}{c}{\left[{\delta }^{n-1}\pi \right]}_{i_1{i}_2\dots i_n}=\frac{1}{r}{\delta }_{i_1{i}_2}{\delta }_{i_1{i}_3}\dots {\delta }_{i_1{i}_n},\end{array}\end{array}$$

and

$$\begin{array}{l}\begin{array}{cc}{\left[\widehat{{\delta }^{n-1}\pi }\right]}_{i_1{i}_2\dots i_n}& =\frac{1}{r}\sum _{j_1,j_2,\dots ,j_r}{\omega }^{i_1{j}_1+i_2{j}_2+\dots +i_n{j}_n}\\ \times {\delta }_{j_1{j}_2}{\delta }_{j_1{j}_3}\dots {\delta }_{j_1{j}_n}\\ =\frac{1}{r}\sum _{j_1}{\omega }^{j_1(i_1+i_2+\dots +i_n)}\\ =\left\{\begin{array}{ll}1& \text{if}{i}_1+i_2+\dots +i_n=0\text{(mod r)}\\ 0,& \text{otherwise.}\end{array}\right.\end{array}\end{array}$$

Translating this result using the ordered-r-partitions for indices, we have

Lemma 3.

In the diagonal basis, the uniform initial distribution on the star tree has components

$$\begin{array}{l}\begin{array}{c}{\left[\widehat{{\delta }^{n-1}\pi }\right]}_{\mu }\\ =\left\{\begin{array}{l}1\text{if}0\left|{\mu }_0\right|+1\left|{\mu }_1\right|+2\left|{\mu }_2\right|+\dots +(r-1)\left|{\mu }_{r-1}\right|=0\mathit{\text{(mod r)}}\\ 0,\text{otherwise.}\end{array}\right.,\end{array}\end{array}$$

where μ=μ ₀: μ ₁: μ ₂:…: μ _r−1 is an ordered-r-partition of the set [n].

Again recall that for this model a generic phylogenetic tensor can be written as

$$\begin{array}{l}\begin{array}{c}P=e^{-\lambda }exp\left(\sum _{\varnothing \ne e\subseteq \left[n-1\right],s\in \left[r-1\right]}\underset{e}{\overset{s}{\alpha }}\underset{{\sigma }^s}{\overset{\left(e\right)}{K}}\right){\delta }^{n-1}\pi ,\end{array}\end{array}$$

where $\pi =\frac{1}{r}{(1,1,\dots ,1)}^T$. In the diagonal basis $\hat{P}:=f^{\left(n\right)}·P$ and as a consequence of Lemma 3$\hat{P}$ will have many vanishing components. To avoid these we take u=u ₀: u ₁: u ₂:…: u _r−1 as an ordered-r-partition of [n−1] and set

$$\begin{array}{l}\begin{array}{cc}\gamma \left(u\right)& =r-\left(0\left|u_0\right|+1\left|u_1\right|+2\left|u_2\right|+\dots \right.\\ \left.+(r-1)\left|u_{r-1}\right|\right)\text{(mod r)}.\end{array}\end{array}$$

If we define ${\mathcal{P}}_u:={\left[\hat{P}\right]}_{u·\gamma \left(u\right)}$ and

$$\begin{array}{l}\begin{array}{c}{\eta }_u:={\left[\sum _{\varnothing \ne e\subseteq \left[n-1\right],s\in \left[r-1\right]}{\alpha }_e^s{\hat{K}}_{{\sigma }^s}^{\left(e\right)}\right]}_{u·\gamma \left(u\right)}^{u·\gamma \left(u\right)},\end{array}\end{array}$$

we then have the first part of the inversion for the ${\mathbb{Z}}_r$ model:

$$\begin{array}{l}\begin{array}{cc}{\mathcal{P}}_u& =e^{-\lambda }exp\left({\eta }_u\right),\\ {\eta }_u& =ln\left({\mathcal{P}}_u\right)+\mathrm{\lambda .}\end{array}\end{array}$$

For each i∈[r−1], we define the column vectors ${\overrightarrow{\alpha }}_i:={\left\{{\alpha }_e^i\right\}}_{\varnothing \ne e\subseteq \left[n-1\right]}$, and, for each ∅≠e⊆[n−1] and u an ordered- (r−1)-partition of [n−1], we define the rectangular r ⁿ⁻¹×2ⁿ⁻¹ matrices

$$\begin{array}{l}\begin{array}{c}\begin{array}{lll}{\left[F_1\right]}_u^e& :={\left[K_{\sigma }^{\left(e\right)}\right]}_{u·\gamma \left(u\right)}^{u·\gamma \left(u\right)},& {\left[F_2\right]}_u^e:={\left[K_{{\sigma }^2}^{\left(e\right)}\right]}_{u·\gamma \left(u\right)}^{u·\gamma \left(u\right)},\dots \\ {\left[F_{r-1}\right]}_u^e& :={\left[K_{{\sigma }^{r-1}}^{\left(e\right)}\right]}_{u·\gamma \left(u\right)}^{u·\gamma \left(u\right)},\end{array}\end{array}\end{array}$$

so we have the vector equation

$$\begin{array}{l}\begin{array}{c}\eta =F_1\overrightarrow{{\alpha }_1}+F_2\overrightarrow{{\alpha }_2}+\dots +F_{r-1}{\overrightarrow{\alpha }}_{r-1}.\end{array}\end{array}$$

We claim that

Lemma 4.

$$\begin{array}{l}\begin{array}{cc}{\left[F_1\right]}_u^e& ={\left[f^{(n-1)}\right]}_u^{e^c\mathrm{:e:\varnothing :\varnothing :}\dots \mathrm{:\varnothing }},\\ {\left[F_2\right]}_u^e& ={\left[f^{(n-1)}\right]}_u^{e^c\mathrm{:\varnothing :e:\varnothing :}\dots \mathrm{:\varnothing }},\cdots \\ \cdots {\left[F_{r-1}\right]}_u^e& ={\left[f^{(n-1)}\right]}_u^{e^c\mathrm{:\varnothing :\varnothing :\varnothing :}\dots \mathrm{:e}}.\end{array}\end{array}$$

Proof.

We recall that ${\left[{\hat{K}}_{{\sigma }^s}\right]}_i^j={\omega }^{\mathit{\text{is}}}{\delta }_{\mathit{\text{ij}}}$, so, for μ=μ ₀: μ ₁: μ ₂:…: μ _r−1 an ordered-r-parition of [n], and e a subset of [n−1] we have

$$\begin{array}{l}\begin{array}{c}{\left[{\hat{K}}_{{\sigma }^s}^{\left(e\right)}\right]}_{\mu }^{\mu }={\omega }^{s\left(0|{\mu }_0\cap e|+1|{\mu }_1\cap e|+\dots +(r-1)|{\mu }_{r-1}\cap e|\right)},\end{array}\end{array}$$

so

$$\begin{array}{l}\begin{array}{c}{\left[{\hat{K}}_{{\sigma }^s}^{\left(e\right)}\right]}_{u·\gamma \left(u\right)}^{u·\gamma \left(u\right)}={\omega }^{s\left(0|u_0\cap e|+1|u_1\cap e|+\dots +(r-1)|u_{r-1}\cap e|\right)},\end{array}\end{array}$$

because e⊆[n−1]. On the other hand ${\left[f\right]}_i^j={\omega }^{\mathit{\text{ij}}}$, so

$$\begin{array}{l}\begin{array}{c}{\left[f^{(n-1)}\right]}_u^{e^c\mathrm{:\varnothing :}\dots \mathrm{:\varnothing :e:\varnothing :}\dots \mathrm{:\varnothing }}={\omega }^{s\left(0|u_0\cap e|+1|u_1\cap e|+\dots +(r-1)|u_{r-1}\cap e|\right)},\end{array}\end{array}$$

where e appears in the s ^th position.

Define, for i∈[r−1], the rectangular 2ⁿ⁻¹×r ⁿ⁻¹ matrices

$$\begin{array}{l}\begin{array}{cc}{\left[G_1\right]}_e^u:& ={\left[{f^{-1}}^{(n-1)}\right]}_{u·\gamma \left(u\right)}^{e^c·\gamma \left(u\right)\mathrm{:e:\varnothing :\varnothing :}\dots \mathrm{:\varnothing }}\\ {\left[G_2\right]}_e^u:& ={\left[{f^{-1}}^{(n-1)}\right]}_{u·\gamma \left(u\right)}^{e^c·\gamma \left(u\right)\mathrm{:\varnothing :e:\varnothing :}\dots \mathrm{:\varnothing }}\\ ⋮\\ {\left[G_{r-1}\right]}_e^u& ={\left[{f^{-1}}^{(n-1)}\right]}_{u·\gamma \left(u\right)}^{e^c·\gamma \left(u\right)\mathrm{:\varnothing :\varnothing :\varnothing :}\dots \mathrm{:e}}.\end{array}\end{array}$$

Of course G _i F _j =δ _ij 1, so we now have the second part of the inversion:

$$\begin{array}{l}\begin{array}{c}\overrightarrow{{\alpha }_i}=G_i\mathrm{\eta .}\end{array}\end{array}$$

Inversion of any abelian group-based model

Lemma 5.

Any (finitely generated) abelian group G is isomorphic to a direct product of cyclic groups of prime-power order, ie. $G\cong {\mathbb{Z}}_{r_1}\times {\mathbb{Z}}_{r_2}\times \dots \times {\mathbb{Z}}_{r_q}$ where each $r_i=p_i^{n_i}$ where p _i is prime and n _i is a positive integer.

Lemma 6.

The group-based model arising from the G is defined only up to group isomorphisms of G.

Proof.

A generic rate matrix for the group-based model arsing from G is given by

$$\begin{array}{l}\begin{array}{c}Q=-\lambda \mathbf{1}+\sum _{e\ne \sigma \in G}{\alpha }^{\sigma }{K}_{\sigma }.\end{array}\end{array}$$

Under a group isomorphism ϕ:G→G ^′, we have ϕ(σ _i σ _j )=ϕ(σ _i )ϕ(σ _j ).

Recall (2), so that the matrix elements ${\left[K_{\sigma }\right]}_i^j$ is set via the action σ _i ↦σ σ _i =σ _j . If we consider the regular representation of G ^′ we then have ${\left[K_{\varphi (\sigma })\right]}_i^j$ defined by ϕ(σ _i )↦ϕ(σ)ϕ(σ _i ). Now ϕ(σ)ϕ(σ _i )=ϕ(σ σ _i )=ϕ(σ _j ) and, because ϕ is a group isomorphism, this occurs if and only if σ σ _i =σ _j . Thus ${\left[K_{\varphi (\sigma })\right]}_i^j={\left[K_{\sigma }\right]}_i^j$ for all i and j.

This means that we can restrict attention to a single representitive in the isomorphism class of G. Of course, for this purpose we choose the representative guaranteed by Lemma 5.

Thus, for any abelian group G, with generators σ ₁,σ ₂,…,σ _q the corresponding group-based model has rate generators given by

$$\begin{array}{l}\begin{array}{c}{L}_{\sigma }=-\mathbf{1}+K_{{\sigma }_1^{m_1}}\otimes K_{{\sigma }_2^{m_2}}\otimes \dots \otimes K_{{\sigma }_q^{m_q}},\end{array}\end{array}$$

for all $e\ne \sigma =\left({\sigma }_1^{m_1},{\sigma }_2^{m_2},\dots ,{\sigma }_q^{m_q}\right)\in G$, where $K_{{\sigma }_i}$ is the permutation matrix representing the generator ${\sigma }_i\in {\mathbb{Z}}_{r_i}$. The character table f of G is simply the tensor product of the individual character tables of the ${\mathbb{Z}}_{r_i}$:

$$\begin{array}{l}\begin{array}{c}f=f_1\otimes f_2\otimes \dots \otimes f_q.\end{array}\end{array}$$

In the diagonal basis we have matrix elements

$$\begin{array}{l}\begin{array}{c}{\left[f_k{K}_{{\sigma }_k^s}{f}_k^{-1}\right]}_i^j={\left[{\widehat{K}}_{{\sigma }_k^s}\right]}_i^j={\left({\omega }_k\right)}^{\mathit{\text{is}}}{\delta }_{\mathit{\text{ij}}},\end{array}\end{array}$$

where ω _k is a k ^th root of unity. Thus

$$\begin{array}{l}\begin{array}{c}{\left[{\widehat{K}}_{{\sigma }_1^{m_1}}\otimes {\widehat{K}}_{{\sigma }_2^{m_2}}\otimes \dots {\widehat{K}}_{{\sigma }_q^{m_q}}\right]}_{i_1{i}_2\dots i_q}^{j_1{j}_2\dots j_q}\\ ={\delta }_{i_1{j}_1}{\delta }_{i_2{j}_2}\dots {\delta }_{i_q{j}_q}{\left({\omega }_1\right)}^{i_1{m}_1}{\left({\omega }_2\right)}^{i_2{m}_2}\dots {\left({\omega }_q\right)}^{i_q{m}_q}.\end{array}\end{array}$$

We write phylogenetic tensors for this model in the form

$$P_{i_{11}{i}_{12}\dots i_{1n},i_{21}{i}_{22}\dots i_{2n}\dots \dots i_{q1}{i}_{q2}\dots i_{\mathit{\text{qn}}}},$$

where 0≤i _sj ≤r _s for all 0≤s≤q. We simplify notation by writing each group of indices as μ ^(s):=i _s1i _s2…i _sn where μ ^(s) is an ordered- r _s -partition of [n].

Lemma 7.

In the diagonal basis, the uniform initial distribution on the star tree has components

$$\begin{array}{l}\begin{array}{c}{\left[\widehat{{\delta }^{n-1}\pi }\right]}_{{\mu }^{\left(1\right)}{\mu }^{\left(2\right)}\dots {\mu }^{\left(q\right)}}\\ =\left\{\begin{array}{ll}1,& \text{if}0\left|{\mu }_0^{\left(i\right)}\right|+1\left|{\mu }_1^{\left(i\right)}\right|+\dots +(r_i-1)\left|{\mu }_{r-1}^{\left(i\right)}\right|=0\forall i;\\ 0,& \text{otherwise.}\end{array}\right.\end{array}\end{array}$$

A generic phylogenetic tensor for this model can be expressed as

$$\begin{array}{l}\begin{array}{c}P=e^{-\lambda }exp\left(\sum _{\varnothing \ne e\subseteq \left[n-1\right]}^{s_i\in \left[r_i-1\right]}\underset{e}{\overset{s_1{s}_2\dots s_q}{\alpha }}\underset{\underset{1}{\overset{s_1}{\sigma }}}{\overset{\left(e\right)}{K}}\otimes \underset{\underset{2}{\overset{s_2}{\sigma }}}{\overset{\left(e\right)}{K}}\otimes \dots \otimes \underset{\underset{q}{\overset{s_q}{\sigma }}}{\overset{\left(e\right)}{K}}\right)·{\delta }^{n-1}\pi ,\end{array}\end{array}$$

where π is the unifrom distribution on $\sum _{i=1}^q{r}_i$ states, i.e.

$$\pi ={\left(\sum _{i=1}^q{r}_i\right)}^{-1}{\left(1,1,\dots ,1\right)}^T.$$

In the diagonal basis $\hat{P}={(f_1\otimes f_2\otimes \dots \otimes f_q)}^{\left(n\right)}·P$, and, as a consequence of the previous lemma, P has many vanishing components. To avoid these, for each i∈[q] we take $u^{\left(i\right)}=u_0^{\left(i\right)}:u_1^{\left(i\right)}:u_2^{\left(i\right)}:\dots :u_{r_i-1}^{\left(i\right)}$ as an ordered- r _i -partition of [n−1] and set

$$\begin{array}{l}\begin{array}{cc}{\gamma }_i\left(u^{\left(i\right)}\right)& =r_i-\left(0\left|u_0^{\left(i\right)}\right|+1\left|u_1^{\left(i\right)}\right|+2\left|u_2^{\left(i\right)}\right|+\dots \right.\\ \left.+(r_i-1)\left|u_{r-1}^{\left(i\right)}\right|\right)\left(\text{mod}r\right).\end{array}\end{array}$$

We then define

$$\begin{array}{l}\begin{array}{c}{\mathcal{P}}_{u^{\left(1\right)}{u}^{\left(2\right)}\dots u^{\left(q\right)}}:={\left[\hat{P}\right]}_{u^{\left(1\right)}·{\gamma }_1\left(u^{\left(1\right)}\right)u^{\left(2\right)}·{\gamma }_2\left(u^{\left(2\right)}\right)\dots u^{\left(q\right)}·{\gamma }_1\left(u^{\left(q\right)}\right)},\end{array}\end{array}$$

and

$$\begin{array}{l}\begin{array}{c}{\eta }_{u^{\left(1\right)}\dots u^{\left(q\right)}}\\ :={\left[\sum _{\varnothing \ne e\subseteq \left[n-1\right]}^{s_i\in \left[r_i-1\right]}{\alpha }_e^{s_1\dots s_q}{\hat{K}}_{{\sigma }_1^{s_1}}^{\left(e\right)}\otimes \dots \otimes {\hat{K}}_{{\sigma }_q^{s_q}}^{\left(e\right)}\right]}_{u^{\left(1\right)}·{\gamma }_1\left(u^{\left(1\right)}\right)\dots u^{\left(q\right)}·{\gamma }_1\left(u^{\left(q\right)}\right)}^{u^{\left(1\right)}·{\gamma }_1\left(u^{\left(1\right)}\right)\dots u^{\left(q\right)}·{\gamma }_1\left(u^{\left(q\right)}\right)},\end{array}\end{array}$$

so that we have the first part of the inversion

$$\begin{array}{l}\begin{array}{cc}{\mathcal{P}}_{u^{\left(1\right)}{u}^{\left(2\right)}\dots u^{\left(q\right)}}& =e^{-\lambda }exp\left({\eta }_{u^{\left(1\right)}{u}^{\left(2\right)}\dots u^{\left(q\right)}}\right),\\ {\eta }_{u^{\left(1\right)}{u}^{\left(2\right)}\dots u^{\left(q\right)}}& =\lambda +ln\left({\mathcal{P}}_{u^{\left(1\right)}{u}^{\left(2\right)}\dots u^{\left(q\right)}}\right).\end{array}\end{array}$$

We define the column vectors ${\overrightarrow{\alpha }}^{s_1{s}_2\dots s_q}:={\left\{{\alpha }_e^{s_1{s}_2\dots s_q}\right\}}_{\varnothing \ne e\subseteq \left[n-1\right]}$ and $\overrightarrow{\eta }:=\left\{{\eta }_{u^{\left(1\right)}{u}^{\left(2\right)}\dots u^{\left(q\right)}}\right\}$ where u _i is an ordered- r _i -partition of [n−1], and we define the (r ₁r ₂…r _q )ⁿ⁻¹×2ⁿ⁻¹ matrices

$$\begin{array}{l}\begin{array}{c}{\left[F_{s_1{s}_2\dots s_q}\right]}_{u_1{u}_2\dots u_q}^e\\ :={\left[K_{{\sigma }_1^{s_1}}^{\left(e\right)}\right]}_{u_1·\gamma \left(u_1\right)}^{u_1·\gamma \left(u_1\right)}{\left[K_{{\sigma }_2^{s_2}}^{\left(e\right)}\right]}_{u_2·\gamma \left(u_2\right)}^{u_2·\gamma \left(u_2\right)}\dots {\left[K_{{\sigma }_q^{s_q}}^{\left(e\right)}\right]}_{u_q·\gamma \left(u_q\right)}^{u_q·\gamma \left(u_q\right)}\\ ={\left[f_1^{(n-1)}\right]}_{u_1}^{e^c\mathrm{:\varnothing :}\dots \mathrm{:\varnothing :e:\varnothing :}\dots \mathrm{:\varnothing }}{\left[f_2^{(n-1)}\right]}_{u_2}^{e^c\mathrm{:\varnothing :}\dots \mathrm{:\varnothing :e:\varnothing :}\dots \mathrm{:\varnothing }}\dots \\ \dots {\left[f_q^{(n-1)}\right]}_{u_q}^{e^c\mathrm{:\varnothing :}\dots \mathrm{:\varnothing :e:\varnothing :}\dots \mathrm{:\varnothing }},\end{array}\end{array}$$

where in each term e appears in the $s_i^{\mathit{\text{th}}}$ position and the equality follows from Lemma 4.

We can then write the vector equation

$$\begin{array}{l}\begin{array}{c}\overrightarrow{\eta }=\sum _{s_1{s}_2\dots s_q:1\le s_i\le r_i-1}{F}_{s_1{s}_2\dots s_q}{\overrightarrow{\alpha }}^{s_1{s}_2\dots s_q}.\end{array}\end{array}$$

If we define the 2ⁿ⁻¹×(r ₁r ₂…r _q )ⁿ⁻¹ matrices

$$\begin{array}{l}\begin{array}{c}{\left[G_{s_1{s}_2\dots s_q}\right]}_e^{u_1{u}_2\dots u_q}\\ ={\left[{f_1^{-1}}^{(n-1)}\right]}_{u_1}^{e^c\mathrm{:\varnothing :}\dots \mathrm{:\varnothing :e:\varnothing :}\dots \mathrm{:\varnothing }}& {\left[{f_2^{-1}}^{(n-1)}\right]}_{u_2}^{e^c\mathrm{:\varnothing :}\dots \mathrm{:\varnothing :e:\varnothing :}\dots \mathrm{:\varnothing }}\dots \\ \dots {\left[{f_q^{-1}}^{(n-1)}\right]}_{u_q}^{e^c\mathrm{:\varnothing :}\dots \mathrm{:\varnothing :e:\varnothing :}\dots \mathrm{:\varnothing }},\end{array}\end{array}$$

where in each term e appears in the $s_i^{\mathit{\text{th}}}$ position, we have the orthogonality relations

$$\begin{array}{l}\begin{array}{c}{G}_{s_1{s}_2\dots s_q}{F}_{s_1^{\prime }{s}_2^{\prime }\dots s_q^{\prime }}={\delta }_{s_1{s}_1^{\prime }}{\delta }_{s_2{s}_2^{\prime }}\dots {\delta }_{s_q{s}_q^{\prime }}\mathbf{1}.\end{array}\end{array}$$

This gives us the second part of the inversion of any group-based model:

$$\begin{array}{l}\begin{array}{c}{\overrightarrow{\alpha }}^{s_1{s}_2\dots s_q}=G_{s_1{s}_2\dots s_q}\overrightarrow{\eta }.\end{array}\end{array}$$

Conclusion

In this article we have given an alternative derivation of the inversion of group-based phylogenetic models. Primarily our method relies on the remarkable intertwining relation between branching events and Markov evolution (4), and the resulting simplified expression of phylogenetic tensors given in (5). From there we took a representation theoretic approach concentrating on the structure of tensor indices.

Authors’ information

Jeremy G. Sumner is a ARC Research Fellow, Peter D. Jarvis is a Alexander von Humboldt Fellow and Barbara R. Holland is a ARC Future Fellow.