# A Phylogenetic Model of the Evolution of Discrete Matrices for the Joint
Inference of Lexical and Phonological Language Histories

Grégoire Clarté, Robin J. Ryder

###### Abstract

We propose a model of the evolution of a matrix along a phylogenetic tree, in
which transformations affect either entire rows or columns of the matrix. This
represents the change of both lexical and phonological aspects of linguistic
data, by allowing for new words to appear and for systematic phonological
changes to affect the entire vocabulary. We implement a Sequential Monte Carlo
method to sample from the posterior distribution, and infer jointly the
phylogeny, model parameters, and latent variables representing cognate births
and phonological transformations. We successfully apply this method to
synthetic and real data of moderate size.

Keywords: Bayesian inference; Historical linguistics; Sequential Monte Carlo;
Phylogenetics

###  1 Introduction

Languages change and diverge through time in a way that has often been
compared to the evolution of biological species. The history of related
languages, at a depth of a few millennia, is often represented by a dated
phylogenetic tree. These trees have traditionally been reconstructed by hand
by expert linguists, using the comparative method Campbell (2013) to
disentangle lexical, phonological and morphosyntactic evidence and deduce
which languages are close cousins. In their pioneering work, Gray and Atkinson
(2003) introduced quantitative methods to infer the phylogeny of the Indo-
European languages, using a model built from genetics. The data they used are
a binarized form of lexical data, recording presence or absence of _cognate
classes_ (defined below). This work has since led to the creation of many
models specifically tailored to cognate data (Nicholls and Gray, 2008; Chang
et al., 2015; Bouckaert et al., 2012; Rama, 2016), which have been applied to
Dravidian (Kolipakam et al., 2018), Sino-Tibetan (Sagart et al., 2019), Pama-
Nyungan (Bowern and Atkinson, 2012), or Austronesian languages (Gray et al.,
2009), amongst others. We refer the reader to Hoffmann et al. (2021) for a
review.

Most phylogenetic methods used in computational linguistics thus rely on
lexical datasets. The raw data are word forms attested in each language for a
number of pre-specified meanings. These raw data are then transformed into a
binary form, called cognacy classes: a cognacy class corresponds to a set of
words in different languages with provably shared history. For example, the
English _father_ , German _Vater_ and Latin _pater_ are cognates. In binary
format, the column corresponding to this cognacy class would bear a $1$ for
these three languages. The similarity between these, and other, words, is
evidence of shared history between these languages.

The transformation of raw word forms into cognacy classes is a time-consuming
process for expert linguists of the language family under study. It typically
transforms the raw data into a binary matrix recording the presence or absence
of each cognacy class in each language; the binary matrix is the only aspect
of the data used in the statistical inference. This allows for easier
stochastic modelling, but at the same time it loses an important part of the
information contained in the raw data, since the phonological aspects of the
data are lost.

Models of the evolution of other aspects of languages, including phonology,
remain scarce, even though this aspect of the data is also highly informative
of language history. List et al. (2017) describe a method to infer cognates
automatically, but this is done independently from the inference of the
phylogeny; Bouchard-Côté et al. (2013) reconstruct phonological changes, but
their method requires a previously reconstructed phylogeny as input. In a
similar vein, Machine Learning algorithms have been applied to the task of
cognate prediction, which proposes a plausible form for a word in a language
based on related languages; see Fourrier and Sagot (2022) and the recent
survey by Dekker and Zuidema (2020). These works do not attempt to jointly
infer the phylogeny of the languages.

A key aspect of cognate-only data is that they are intrinsically one-
dimensional: the observed data at a single language consist of a vector of
presence/absence of each cognate. (The whole data set with several languages
can then be seen as a matrix, with one language per row.) Most models then
assume that the vector elements evolve independently, allowing for likelihood
evaluations through a pruning algorithm Felsenstein (1981).

This paper describes the first attempt at modelling jointly the changes in
cognates and phonology, thus describing the evolution of the whole word form
and not just the cognacy classes. In our setting however, the data at a single
language are two-dimensional: they consist of a matrix, with each row
corresponding to a meaning, and each column corresponding to a phoneme. The
set of phonemes is finite and known; in this paper, we assume that the data
are aligned, in a way discussed later on. The whole data set with several
languages is thus a 3-dimensional array, and each array entry is a phoneme.
Indeed, the lexicon of a language is modified from its parent language via two
major processes: (1) a new word may arise, corresponding to the birth of a new
cognacy class; (2) existing words may be modified by a phonological change,
such as $/p/\to/f/$. Process (1) is well modelled by existing works; for
matrix data, it corresponds to the replacement of all values in a row of the
matrix. The difficulty with process (2) is that it typically affects
simultaneously many meanings in the language (Campbell, 2013): several values
in a column of the matrix will be modified.

By modelling both the arisal of new cognates and the phonological changes, we
are able to infer jointly the language phylogeny, cognate births and deaths,
and phonological changes. We present a stochastic model of the evolution of a
discrete matrix on a tree in Section 3, and describe in Section 4 numerical
methods to sample from the associated posterior. The application motivating
this work is the reconstruction of the history of Polynesian languages, which
we describe and analyse in Section 5, along with examples on simulated
datasets.

###  2 Background

,

meaning | eye | seven | hand | flower | to breathe | fruit
---|---|---|---|---|---|---
Samoan | mata | fitu | lima | fuŋa | mānava | fua
Tahitian | mata | hitu | rima | pua | huti | hua
Hawaiian | maka | hiku | rima | pua | hanu | hua
Table 1: Toy lexical dataset $D$ for three Austronesian languages: Samoan, Tahitian and Hawaiian. Cognates are displayed in the same colour. Data adapted from Greenhill et al. (2008). eye1 |  $\\{$mata, mata, maka$\\}$
---|---
seven1 |  $\\{$fitu, hitu, hiku$\\}$
hand1 |  $\\{$lima, rima, rima$\\}$
flower1 |  $\\{$fuŋa$\\}$
flower2 |  $\\{$pua, pua$\\}$
breathe1 |  $\\{$mānava$\\}$
breathe2 |  $\\{$huti$\\}$
breathe3 |  $\\{$hanu$\\}$
fruit1 |  $\\{$fua, hua, hua$\\}$
Table 2: Cognate sets $D/\hcrossing$ for the data of Table 1. |  eye1 |  seven1 |  hand1 |  flower1 |  flower2 |  breathe1 |  breathe2 |  breathe3 |  fruit1
---|---|---|---|---|---|---|---|---|---
Samoan | 1 | 1 | 1 | 1 | 0 | 1 | 0 | 0 | 1
Tahitian | 1 | 1 | 1 | 0 | 1 | 0 | 1 | 0 | 1
Hawaiian | 1 | 1 | 1 | 0 | 1 | 0 | 0 | 1 | 1
Table 3: Transformed data $\tilde{D}$, based on the data of Table 1 listing
presence/absence of each cognacy class.

The vast majority of quantitative studies for Historical Linguistics rely on
cognate data, which are a transformation of lexical data. A set of $L$
languages and of $M$ meanings are chosen in advance. Linguists collect the
word form(s) for each language $\ell$ and meaning $m$. This gives rise to raw
data $D$, in the form of a $L\times M$ matrix of words, exemplified in Table
1. We say two words $w=D_{\ell,m}$ and $w^{\prime}=D_{\ell^{\prime},m}$ in two
different languages are cognates, and write $w\hcrossing w^{\prime}$, if they
have the same meaning and have evolved from a common ancestral form. For
example, the Samoan, Tahitian and Hawaiian words for _eye_ are cognates:
"mata"$\hcrossing$"mata"$\hcrossing$"maka". Similarly, the three words for
_seven_ are cognates: "fitu"$\hcrossing$"hitu"$\hcrossing$"hiku". This is
shown by noticing not only the similarities between the words, but also that
differences occur in a systematic fashion: an $/f/$ in Samoan often
corresponds to an $/h/$ in Tahitian and Hawaiian (see _fruit_), and a $/t/$ in
Samoan and Tahitian often corresponds to a $/k/$ in Hawaiian (see _eye_). One
possible explanation would be that in the past, there was an $/h/\to/f/$
phonological transformation on the branch leading to Samoan, after it diverged
from Tahitian and Hawaiian, and a $/t/\to/k/$ transformation on the branch
leading to Hawaiian.

Since $\hcrossing$ is an equivalence relation, we can define equivalence
classes for $\hcrossing$, called cognacy classes, with $c_{m}$ classes for
meaning $m$; this corresponds to the quotient space $D/\hcrossing$ and an
example is given in Table 2. Column $m$ of $D$ is then transformed into a
$L\times c_{m}$ binary matrix, where a 1 in column $i$ for language $\ell$
means that cognacy class $i$ is represented in language $\ell$, and a 0 means
it is not represented in language $\ell$. We concatenate these binary matrices
to obtain the transformed data $\tilde{D}$, a $L\times\sum_{m}c_{m}$ binary
matrix, as shown in Table 3. It is $\tilde{D}$ which is then used in the
statistical analysis. Note that in $\tilde{D}$, we typically do not use the
original meanings, and in particular we ignore the fact that certain columns
of $\tilde{D}$ correspond to the same meaning; we then assume that the columns
of $\tilde{D}$ are iid.

The transformation $D\mapsto\tilde{D}$ is a time-consuming process performed
by hand by linguists expert of the language family under study. This
binarization also discards much potentially useful information. In particular,
the phonological transformations (which contain signal on the phylogeny) could
be recovered from $D$, but not from $\tilde{D}$.

In the format $\tilde{D}$, the data observed at a leaf language are assumed to
be a binary vector of cognate presence/absence, with vector elements assumed
to be independent. The models can thus be represented by a generative process
with state space $\\{0,1\\}$, of which we observe multiple independent
realizations. In our work however, the data at a single language are two-
dimensional: they consist of a matrix, in which each row corresponds to a
meaning, and each column corresponds to a phoneme. The set of phonemes is
finite and known. The whole data set with several languages is thus a
3-dimensional array, and each array entry is a phoneme.

Formalizing our notation, let $D_{\ell}$ be a matrix representing the data
observed at leaf language $\ell$. Table 4 shows how some of the raw data of
Table 1 are represented. The value of $D_{\ell}(m,k)$ corresponds to the $k$th
phoneme of the word for meaning $m$. It belongs to a known discrete set
identified with the integers $\\{1\ldots n_{k}\\}$. The number of possible
elements $n_{k}$ may vary between columns; this is the case in our application
of Section 5.

Samoan | Tahitian | Hawaiian
---|---|---
| eye | m | a | t | a
---|---|---|---|---
seven | f | i | t | u
hand | l | i | m | a
| eye | m | a | t | a
---|---|---|---|---
seven | h | i | t | u
hand | r | i | m | a
| eye | m | a | k | a
---|---|---|---|---
seven | h | i | k | u
hand | r | i | m | a
Table 4: Representation of part of the data of Table 1 as a collection of
aligned matrices.

###  3 Model

We propose a generative model of the evolution of a discrete matrix on a
rooted and dated binary tree.

####  3.1 Data format

We consider data collected at $L$ leaves $\ell=1\ldots L$. At leaf language
$\ell$, the data are a $M\times K$ matrix, so that $D_{\ell}(m,k)$ is the
value for character $k$ of the meaning $m$ in language $\ell$. The set of
possible values for $D_{\ell}(m,k)$ is assumed known and finite; we identify
it with $1,\ldots,n_{k}$. Different characters may have different sets of
possible values.

####  3.2 Evolution model

We represent the vocabulary at point $x$ on the tree by a matrix $U_{x}$. At
the root, the word forms are drawn randomly from the uniform distribution
$\pi_{0}$ over all forms. At a branching event on the tree, the two child
branches inherit the vocabulary of the parent branch, and then evolve
independently. We model word change by two independent Poisson processes on a
phylogenetic tree; these processes modify $U$. These processes are:

  * •

a creation process: at a creation event, a new word is born. A meaning is
drawn at random and the entire corresponding row is replaced, with new values
drawn from the uniform distribution $\pi_{0}$. Creation events arise from a
Poisson Process with rate $\mu$, independently for each meaning. This process
is similar to the $M_{k}$ model (Pagel, 1994; Lewis, 2001).

  * •

a transformation process: at a transformation event, a phonological change $r$
takes place, of the form $r:/a/\to/a^{\prime}/$, affecting phoneme position
$k$. In column $k$, across all meanings, each occurrence of $/a/$ is replaced
by $/a^{\prime}/$ with probability $q_{k}$ independently. For each character
$k$, transformation events arise from a Poisson process with rate
$\lambda_{k}$. The change $r$ is then drawn from a catogerical distribution
over a predetermined set of possible transformations $\mathfrak{R}_{k}$. This
is a new adjunction to the language phylogeny literature, as we are aiming at
simulating the natural evolution of phonemes through time, as described by
phonological rules (Hock, 1991). This model, to our knowledge is not similar
to another one existing in the literature, because of the regularity of the
transformations.

We will identify a transition $r$ with its transition matrix, defined for
$a\rightarrow b$ as $r_{i,j}=\left\\{\begin{matrix}1&\text{if }i=j,i\neq a\\\
\beta_{k}&\text{if }i=a,j=b\\\ 1-\beta_{k}&\text{if }i=a,j=a\\\
0&\text{otherwise}\end{matrix}\right.$.

Figure 1 represents these processes and their effect on a small lexicon.

Let $x$ be the position of a creation event. We denote by $U_{x^{-}}$ and
$U_{x^{+}}$ the value of $U$ just before $x$ and just after $x$, respectively.
The creation event changes $U$ as follows:

$m\sim\mathcal{U}(1,\ldots,M)\qquad U_{x^{+}}(m,\cdot)\sim\pi\qquad\forall
m^{\prime}\neq m,U_{x^{+}}(m^{\prime},\cdot)=U_{x^{-}}(m^{\prime},\cdot).$ (1)

Now let $t$ be the position of a transformation event. We draw

$k\sim\mathcal{U}(1,\ldots,K)\qquad(a\rightarrow
a^{\prime})=r\sim\mathcal{U}(\mathfrak{R}_{k})$

so that the transformation is $/a/\to/a^{\prime}/$ in column $k$. The
transformation event changes $U$ as follows:

$\displaystyle\forall k^{\prime}\neq k,\forall m,$ $\displaystyle
U_{x^{+}}(m,k^{\prime})=U_{x^{-}}(m,k^{\prime})$ (2) $\displaystyle\forall
m,\text{ if }U_{x^{-}}(m,k)\neq a,$ $\displaystyle
U_{x^{+}}(m,k)=U_{x^{-}}(m,k)$ $\displaystyle\forall m,\text{ if
}U_{x^{-}}(m,k)=a,$ $\displaystyle
U_{x^{+}}(m,k)=\begin{cases}a^{\prime}&\text{ w.p. }\beta_{k}\\\ a&\text{ w.p.
}1-\beta_{k}.\end{cases}$

To account for unmodelled events, we allow for noise: at the end of an edge of
length $h$, we let

$\forall m,k,U_{x^{+}}(m,k)\begin{cases}\sim\mathcal{U}(1,\ldots,n_{k})&\text{
w.p. }1-e^{-h\nu}\\\ =U_{x^{-}}(m,k)&\text{ otherwise.}\end{cases}$ (3)

The noise level $\nu$ is useful to measure model fit. We also rely upon it in
the numerical implementation (Section 4) to define a tempering strategy.

This fully specifies the stochastic process $U$. The data are then observed at
the leaves: at leaf $\ell$, $U_{\ell}=D_{\ell}$.

Figure 1: Stochastic evolution processes. The first event corresponds to a
phonological transformation $/p/\to/f/$: in column 1 of $U_{x}$, all
occurrences of $/p$ may transform into $/f/$. The second event corresponds to
the creation of a new cognate: the second row of of $U_{x}$ is entirely
renewed.$\left[\begin{tabular}[]{ccc}p&a&n\\\
p&i&k\end{tabular}\right]$$\times$$\bullet$$\mathfrak{R}_{1}\ni
r:p\overset{\beta_{1}}{\rightarrow}f$$\left[\begin{tabular}[]{ccc}{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}f}&a&n\\\
{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}f}&i&k\end{tabular}\right]$$\times$$\bullet$$B^{1}_{e}$$\left[\begin{tabular}[]{ccc}f&a&n\\\
f&i&k\end{tabular}\right]$$\times$$\left[\begin{tabular}[]{ccc}{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}l}&{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}i}&{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}r}\\\
f&i&k\end{tabular}\right]$$\times$

To summarise, the whole process is described by:

  * •

$g$ the topology of the phylogenetic tree, with edges $E$.

  * •

the associated set of edges length $(h_{e})_{e\in E}$

  * •

$\mu$ the renewal rate of the meanings.

  * •

$\lambda_{k}$ the rate of the appearance of the transformations for each
traits; we write $\bm{\lambda}=\left(\lambda_{k}\right)_{k}$

  * •

$p_{k}$ a vector of probabilities for each transformation of each character;
we write $\bm{p}=\left(p_{k}\right)_{k}$

  * •

$\beta_{k}$ the probability that a transformation is applied; we write
$\bm{\beta}=\left(\beta_{k}\right)_{k}$

  * •

$\nu$ the intensity of the noise.

####  3.3 Questions of interest

As often in phylogenetic models, the dated tree is the main parameter of
interest. However several other question arise once the evolutive model is
set. The evolution rates $\lambda$ and $\rho$ in particular, as they inform on
evolution speed, can be compared between several subgroups of taxa. The
parameter $\beta$ is also interesting as it corresponds to different models
depending on its value. For $\beta$ close to $1$, the model is close to the
classical phonological model that assumes that phonological transformations
are _regular_ ; by opposition for $\beta$ small, transformations are not
regular.

####  3.4 Likelihood computation

$x$$y$$z$$Q(x,\cdot,\cdot)$$Q(y,\cdot,\cdot)=(S_{\langle
x,y\rangle}Q(y,\cdot,\cdot))\odot(S_{(x,z)}Q(z,\cdot,\cdot))$$Q(z,\cdot,\cdot)$$r_{5}$$r_{7}$$\bullet$$\bullet$$r_{2}$$r_{6}$$r_{3}$$\bullet$$\bullet$$\bullet$$r_{7}r_{5}S_{noise}(\ell_{\langle
x,y\rangle})=S_{\langle
x,y\rangle}$$S_{(x,z)}=r_{2}r_{6}r_{3}S_{noise}(\ell_{(x,z)})$ Figure 2:
Computation of the likelihood, without cognacy apparition, where $\odot$
denotes elementwise matrix multiplication.

The likelihood is not available in closed form because of the intricacy of the
two evolution processes. We thus add two sets of latent variables in order to
make the computations tractable:

  * •

$R^{k}_{e}$ the list of the dated transformations for character $k$ on branch
$e$, with associated position on the branch; we write
$\bm{R}^{k}=\left(R^{k}_{e}\right)_{e}$

  * •

$B_{e}^{m}$ the position of each apparitions for each meaning $m$ on branch
$e$; we write $\bm{B}=\left(B_{e}^{m}\right)_{e,m}$

Given these auxiliary variables, the likelihood can be factorized: denoting
$\theta=(g,h,\mu,\bm{\lambda},\bm{p},\bm{\beta},\nu,\bm{q})$ the set of all
parameters and auxiliary variables, and $D^{k}=D(\cdot,\cdot,k)$, we have

$\displaystyle\mathcal{L}(D\mid\theta,\bm{R},\bm{B})=\prod_{k}\mathcal{L}(D^{k}\mid\theta,\bm{R}^{k},\bm{B}).$
(4)

Then each term $\mathcal{L}(D^{k}\mid\theta,\bm{R}^{k},\bm{B})$ can be
computed separately by pruning (Felsenstein, 1981), as follows. In the
remainder of this subsection, we fix a character $k$. The quantity on which we
will prune is $\mathbb{P}(D_{desc(x)}(m,k)\mid U_{x}(m,k)=u,\theta)$ for a
node $x$ on the tree, where $desc(x)$ are the descendant leaves of $x$ and
$U_{x}(m,k)$ is the value for character $k$ and meaning $m$.For the sake of
simplicity we drop the dependency in $k$ and $\theta$, and define
$Q(x,m,u)=\mathbb{P}(D_{desc(x)}(m,k)\mid U_{x}(m,k)=u,\theta)$. The
likelihood for character $k$ is then:

$\mathcal{L}(D^{k}\mid\theta,\bm{R}^{k},\bm{B})=\prod_{i=1}^{M}\sum_{u}\pi_{0}(u)Q(root,m,u).$

The quantity $Q$ can be computed recursively through branching events of the
tree: let $x$ be an internal node with children nodes $y$ and $z$; we have

$\displaystyle Q(x,m,u)$ $\displaystyle=$
$\displaystyle\sum_{v}Q(y,m,v)P(U_{y}(m)=v\mid U_{x}(m)=u)\times$
$\displaystyle\qquad\sum_{w}Q(z,m,w)P(U_{z}(m)=w\mid U_{x}(m)=u).$

We are left to compute the quantity $P(U_{y}(m)=w\mid U_{x}(m)=u)$ going down
a branch $\langle x,y\rangle$. There are two cases to consider for this
computation.

###### Case 1: no cognate apparition

We first consider the case where there is no cognate apparition on branch
$\langle x,y\rangle$. In this case, we have to account for both the
transformations and the noise. We define $S$, the transition matrix of the
value for a given character, by:

$S_{\langle x,y\rangle}=\left(\prod_{r\in R_{\langle
x,y\rangle}}r\right)\times S_{noise}(h_{\langle x,y\rangle}),$ (6)

where $S_{noise}(h_{\langle x,y\rangle})$ is the matrix associated with the
noise at the end of the edge $\langle x,y\rangle$ from $x$ to $y$, as in eq.
3, and where the terms of the product are considered in the order they appear
along the branch. In this case, $P(U_{y}(m)=w\mid U_{x}(m)=u)$ is equal to the
entry at position $(u,w)$ of the matrix $S_{\langle x,y\rangle}$. Figure 2
summarizes the relationship between these objects.

$x$$y$$z$$Q(\cdot,y,\cdot)$$Q(x,m^{,}u)=\sum_{v}\sum_{w}\left[\pi_{\langle
x,y\rangle,\tau}(v)Q(m,y,v)\right]Q(m,z,w)S_{(x,z)}(u,w)$$Q(\cdot,z,\cdot)$$r_{5}$$r_{7}$$\bullet$$\bullet$$\pi_{\langle
x,y\rangle,\tau}=\pi_{0}r_{5}S_{noise}(\ell_{\langle
x,y\rangle})$$\bullet$$S_{(x,z)}$$\tau$ Figure 3: Computation of the
likelihood for meaning $m$, with an apparition on the edge $\langle
x,y\rangle$.

###### Case 2: at least one cognate apparition

The second case to consider is if there is apparition of a cognate on the edge
$\langle x,y\rangle$ for meaning $m$. In this case, the value at the child
node $U_{y}(m)$ is independent of the value at the parent node $U_{x}(m)$. Let
$\tau$ be the last apparition time of a cognate for meaning $m$ on branch
${\langle x,y\rangle}$. The distribution for character $k$ of meaning $m$ at
the end of the branch $y$ only depends on the transformations $\\{r\in
R_{{\langle x,y\rangle}}^{k}:t(r)<\tau\\}$, and is given by:

$\pi_{{\langle x,y\rangle},\tau}=\pi_{0}\times\left(\prod_{r\in R_{{\langle
x,y\rangle}}^{k},\ t(r)<\tau}r\right)\times S_{noise}(h_{{\langle
x,y\rangle}})$ (7)

and $P(U_{y}(m)=v\mid U_{x}(m)=u)$ is equal to the $v$th entry of vector
$\pi_{{\langle x,y\rangle},\tau}$. We represent these computations in Figure
3.

We initialize $Q$ at the leaves: $Q(x,m,u)=\mathbf{1}_{D_{x}(m)=u}$ if $x$ is
a leaf. Taken together, Equations 6 and 7 allow us to compute the quantity $Q$
on which we prune in equation 3.4. We are thus able to calculate the
likelihood.

####  3.5 Missing data

The computations above assume that all data have been observed. Our model
offers the advantage of handling missing directly, without any more cost.The
matrix $Q$ at a leaf $\ell$ is initialized with the observations; it
represents the probability, knowing the state at the leaf to observe the
observations at this leaf. In the case of missing data, since the pattern of
missingness is observed, it is straightforward to modify the leaf
initializations as

$\displaystyle Q(\ell,m,k)$ $\displaystyle=P(D_{\ell}(,m,k)\mid
U_{\ell}(m,k),\text{missingness})$
$\displaystyle=\left\\{\begin{matrix}1&\text{ if missing}\\\ 0&\text{ if not
missing and }D_{\ell}(m,k)\neq U_{\ell}(m,k)\\\ 1&\text{ if not missing and
}D_{\ell}(m,k)=U_{\ell}(m,k)\end{matrix}\right.$

####  3.6 Prior distributions

#####  3.6.1 On the evolution parameters

We take conjugate prior distributions on the evolution parameters. For
$\lambda$, $\mu$ and $\nu$ we chose gamma priors fitted so that on each branch
there are, in expectation, less than ten transformations or apparition, in
order to keep the model in reasonable.

We use conjugate priors for $\lambda$, $\mu$, $\nu$ and $p$ as these are the
parameters associated with Poisson and multinomial distributions. In our
applications below, we take priors $\beta_{k}\sim Beta(1,10)$ iid and
$\pi(\nu)\propto 1/\nu\mathbf{1}_{[a,b]}(\nu)$.

#####  3.6.2 On the tree

Following Nicholls and Gray (2008), we use a prior on the trees of the form:

$\pi(g)\propto\mathbf{1}_{t(root)\in[a,b]}t(root)^{L},$

where $L$ is the number of leaves and $a$ and $b$ are bounds on the root age.
We further constraints this prior by adding constraints on clades and their
ages. Notice that learning the root age is only possible if enough internal
constraints are provided, knowing this, we also propose to fix the root age to
an arbitrary value, so that we can focus on learning topologies and relative
ages — which does not depend on the root age. As appealing as learning the
root age is, we advise against putting too much faith in the posterior ; in
general this age is little informed, resulting only in the propagation of a
few constraints on the leaves. This was, anyway not a primary concern of our
study.

We also tried to add Gamma prior on the total length of the tree, that is

$\pi(g)\propto\Gamma_{c,d}\left(length(g)\right)\mathbf{1}_{t(root)\in[a,b]}t(root)^{L},$

where $\Gamma_{c,d}$ is the density of a Gamma distribution with parameters
$c$ and $d$. The results where similar compared to the uniform prior, which
indicates a low effect of the prior.

####  3.7 Inference of ancestral states

As stated in Section 3.4 the matrix $Q$ contains the probability of the
observations given the state of each character for each meaning. Given the
latent variables $R$ and $L$, it is possible to compute at any node $x$ of the
tree $\pi_{x}^{k}(m)$ the distribution of the values for meaning $m$, and
character $k$ at $x$, by descending the tree from the root, at which
$\pi_{r}^{k}(m)=\pi_{0}$, to $x$. We can thus compute exactly, with $n_{k}$
the number of possible values for character $k$:

$\pi(U(x,m,k)\mid
R,L,D)=\frac{Q(x,m,k)\odot\pi_{x}^{k}(m)}{\sum_{i=1}^{n_{k}}Q(x,m,k)(i)\pi_{x}^{k}(m)(i)},$

the renormalization is computed easily as a sum over a dozen elements. Notice
that in this case, when conditioning on $R$ and $L$, the value for each
meaning and each character is independent. If we are only interested in the
distribution of the ancestral value for a given meaning and character, this is
of little interest, but if we are interested in the _joint_ distribution of
several characters and meanings, care must be taken when renormalizing the
conditional posteriors.

The posterior distribution is then:

$\int_{R,L}\pi(U_{x}(m,k)\mid R,L,D)\pi(R,L\mid D)\mathrm{d}(R,L),$

which can be approximated by Monte Carlo, as we already have a sample from
$\pi(R,L\mid D)$.

###  4 Numerical methods

All the code is available, along with the datasets used, at
https://github.com/GClarte/PhylogenyFromMatrices.

We use a Sequential Monte Carlo (SMC, Del Moral et al. (2006)) algorithm to
sample from the posterior. SMC relies on importance sampling principles. This
method requires us to define a sequence of distributions $\pi_{t}$ that
creates a smooth transition between an initial known distribution and the
posterior. The usual method is to use a tempering of the form
$\pi_{t}(x)\propto\pi(x\mid D)^{1/t}$ or $\pi_{t}(x)\propto\pi(x)L(D\mid
x)^{t}$. This strategy is impractical in our case, as our mutation kernel
makes use of a Gibbs sampler, disallowing such tempered intermediate
distributions; none of our attempts in this direction were fruitful. Instead,
we use the noise parameter as a tempering parameter. For a high value of
$\nu$, on each branch the noise matrix $S_{noise}$ is close to the matrix
$(\pi_{0},\dots,\pi_{0})$, that is as if the model was a forest of trees
constituted of only one leaf. Conversely, for $\nu$ small we retrieve the
model we are looking for.

As the noise is a parameter that still needs to be learned, we propose the
following distributions. For a sequence $\nu_{1}>\dots>\nu_{T}$ we define :

$\pi_{t}(g,\mu,\lambda,p,\beta,q,\nu\mid
D)=\pi(g,\mu,\lambda,p,\beta,q\mid\nu=\nu_{t},D),$

with a final target

$\pi_{T+1}(g,\mu,\lambda,p,\beta,q,\nu\mid
D)=\pi(g,\mu,\lambda,p,\beta,q,\nu\mid D).$

That is we fix the noise parameter and we reduce it until it reaches a value
that will lie in a part of the space with non negligible posterior
probability.

This does not ensure a proper mixing of the SMC algorithm, as the number of
tempering steps needed can be quite high — and a high number of steps would
require a higher number of particles, which becomes quite heavy in terms of
memory costs — leading to a genealogy degeneracy. This is a severe limitation
to the use of our current version of the method, and this prevents us from
estimating the marginal likelihood of the model. Depsite the particle
degeneracy, we notice that the last particles do mix adequately. Our chosen
solution is thus to start with an SMC run, and use the final step to
initialize a Hastings-Metropolis-Rosenbluth-Teller algorithm (Metropolis et
al., 1953; Hastings, 1970). We describe thoroughly in the supplementary the
proposal kernel we used, the main difficulty lying in the handling of the
latent variables.

This method allows for parallelization, but it is at the cost of memory
efficiency, because of the size of the latent variables: we have to store the
value of $Q$ at each node (that is, a matrix of size $n_{k}\times M$ with
$n_{k}$ number of possible values for character $k$ and $M$ number of
cognates), for each character. In particular the number of particles is
problematic. On simulated data, 2000 particles are enough to provide stable
results. On real dataset, with larger inventories for each character the
number of particles cannot be higher than 1500, and the results are sometimes
unstable. We are aware of this issue, and have for now no other solution than
increasing the available memory.

####  4.1 Rooting the tree

This SMC scheme is able to adequately explore the posterior; with one
exception. Analysis of the output of the initial implementation of this
algorithm showed that the split closest to the root (_i.e._ the choice of
outgroup) mixed poorly in several examples. To circumvent this issue, we added
another taxon to the dataset with only missing data, constraining this taxon
to be an outgroup. The results following this method have been satisfying,
although we note that this a consequence of this choice is a small change in
the prior on the MRCA associated with the initial languages, that is the root
of the studied taxa.

###  5 Applications

To demonstrate the applicability of our method, we first apply it to synthetic
datasets in Section 5.1, and show that we are able to correctly reconstruct
parameters of interest. We then apply it to a subset of Polynesian languages
in 5.2.

We point out that our model is designed for particular datasets which present
the desired alignment properties. For many language families, this alignment
is not verified, and so our model would only be applicable with difficulty.
Nonetheless, there are numerous cases that are of interest. For example, many
Tai languages, which present monosyllabic structure, correspond to our
framework. Below, we apply our method to certain meanings from a set of well-
chosen Polysenian languages, which all fit a CVCV pattern.

Noticeably, these examples corresponds to closely related languages, that
often have only a few cognacy classes to differ, that is a very faint cognacy
signal. In such cases, using the phonological signal is essential.
Additionally, these data allow for an additional check: it is possible to
check the cognacy reconstruction by comparing with the linguists’
reconstruction.

####  5.1 On synthetic data

First, we can run numerical methods on synthetic datasets. Some experiments
with data from misspecified models are available in supplementary material.

For this experiment, we chose 4 characters with the following specifications :

Number of value | Possible transitions | $\lambda$ | $\pi_{0}$ | $p$ | $\beta$
---|---|---|---|---|---
5 | $1\leftrightarrow 2\leftrightarrow 3\leftrightarrow 4\leftrightarrow 5\leftrightarrow 1\leftrightarrow 4$ | $10^{-2}$ | $\mathcal{U}\\{1,\dots,5\\}$ | uniform | 1
5 | $1\leftrightarrow 2\leftrightarrow 3\leftrightarrow 4\leftrightarrow 5\leftrightarrow 1\leftrightarrow 4$ | $10^{-2}$ | $\mathcal{U}\\{1,\dots,5\\}$ | uniform | 1
2 | $1\leftrightarrow 2$ | $10^{-3}$ | $\mathcal{U}\\{1,2\\}$ | uniform | 1
2 | $1\leftrightarrow 2$ | $10^{-3}$ | $\mathcal{U}\\{1,2\\}$ | uniform | 1

The true tree is represented in top of Figure 4. We constrain leaves
$\\{8,12,13,15\\}$ and $\\{9,11,10\\}$ to form clades and we added the
following constraints on some internal node ages: the nearest common ancestor
of $10$ and $9$ in $[10,200]$, of $1$ and $2$ in $[10,200]$, of $8$ and $15$
in $[20,30]$. Different constraints on the internal nodes led to very similar
results.

We generated $100$ data points with these characters on a $15$-leaved tree
presented on the left in Figure 4. We ran the previously described SMC sampler
with a twelve step tempering on $\nu$ : $\\{1\cdot 10^{-1},\ 5\cdot 10^{-2},\
3\cdot 10^{-2},\ 2\cdot 10^{-2},\ 1\cdot 10^{-2},\ 9\cdot 10^{-3},\ 8\cdot
10^{-3},\ 7\cdot 10^{-3},\ 6\cdot 10^{-3},\ 5\cdot 10^{-3},\ 1\cdot
10^{-3}\\}$. We used $2\cdot 10^{3}$ particles, for the mutation step we ran
$10^{3}$ iterations of our Metropolis within Gibbs algorithm. After the last
tempering, we ran a last mutation step of $10^{4}$ iterations for each
particles. In total, on $30$ CPU Intel Xeon CPU E5-2630 v4 2.20GHz, it took
roughly 5 hours.

Figure 4: True tree and consensus tree of the final sample, simulated dataset.
Internal nodes of the true tree are tagged with their posterior probability.

The consensus tree is represented on the right in Figure 4. The overall
topology is well reconstructed, with uncertainty on the smaller branches,
which is expected as the small branches rarely host even a single
transformation. Noticeably each group of leaves supported by the consensus
tree corresponds to a true clade. Node ages do not seem accurately
reconstructed, in accordance with the root being to high.

The parameters are also well reconstructed; details are given in Supplementary
material. The existence of apparition events on the path between two leaves is
also quite well reconstructed, which has an interesting interpretation.

####  5.2 Polynesian languages

We now present a preliminary application of this method to data for 78
meanings from 14 Polynesian languages. The raw data form a subset of the
online Austronesian Basic Vocabulary Database (ABVD) (Greenhill et al., 2008).
We selected a subset of languages and meanings for which many word forms were
4 characters long, of the form Consonant-Vowel-Consonant-Vowel, to allow for
perfect alignment. This preprocessing step led us to exclude from the data all
words which are not of length 4 and replace them by missing data. This step
biases the data, and we therefore do not claim the inferred topology is
reliable. The cognacy judgments will also be biased, but we might hope
nonetheless that the inferred cognacies are somewhat trustworthy.

Despite this caveat, the results are interesting. In particular, the evolution
parameters in Figure 5 show a similarity between the evolution rate of the two
vowels, and between the two consonants. The noise parameter is relatively low,
indicating that the model is indeed detecting a signal from the dataset. The
topology is summarized by a majority-rule consensus tree in Figure 6). The
position of Hawaiian as an outlier in 54% of the posterior is unexpected, and
might be due to a bias in the way we selected the data. The other Eastern
Polynesian languages (Marquesan, Mangareva, Tuamotu, Rarotongan, Tahitian and
Rapanui) form a clade with posterior probability 99%, separating from the
other languages (which correspond to the Samoic languages). The Futunic
languages (Rennellese, Tikopia and Futuna Aniwa) also form a strongly
supported clade.

Figure 5: Parameter estimations from the ABVD dataset Figure 6: Consensus tree
from the ABVD dataset

The ABVD includes cognate sets by expert linguists for most pairs of words, so
we compare the inferred probability of cognacy with the true cognacy status.
There is overall agreement between the inferred and true cognacy status of
pairs of words, as shown in Figure 7. If we use the decision rule that two
words are inferred to be cognates iff the probability of cognacy is $>0.5$, we
correctly classify the cognacy status of 88% of word pairs.

In ongoing work not presented in this preprint, we apply this method to other
language families which correspond better to our alignment requirement.

Figure 7: Estimated probability of cognacy for pairs of words in the ABVD
dataset. Left: pairs of words which are not cognates according to expert
linguists; right: pairs of words which are true cognates.

###  6 Discussion

This works presents a new step towards an automatic version of the
quantitative method, which is arguably one of the ultimate goals of
computational phylolinguistics: learning cognacy classes along with the
phylogeny and evolution parameters.

This model is of course an imperfect representation of the true evolution of
lexicon and phonology. The first glaring simplification is that we do not
allow for the phonological transformations to depend on context. (Instead,
each transformation applies with probability $\beta_{k}\in]0,1[$, which allows
for transformations to not apply to all words.)

On the numerical side, the memory cost of our method — that have been
increased to decimate the time cost — remains the main limitation, computing
the marginal likelihood of our model is also out of reach for now.

An attainable extension of our model would be to allow for some or all cognacy
judgments to be known, constraining the latent variables $\bm{B}$, but still
inferring the phonological transformation to better inform the phylogeny.
Further down the road, future work will need to address modelling context-
dependent phonological changes, as well as insertions and deletions of
phonemes, and then inferring under such a model in an efficient manner.
Nonetheless, we believe this work to be an important step towards this
ultimate goal, and to already allow inference of language histories for
lesser-studied families for which no cognate data is available.

### References

  * Bouchard-Côté et al. [2013] Alexandre Bouchard-Côté, David Hall, Thomas L. Griffiths, and Dan Klein.  Automated reconstruction of ancient languages using probabilistic models of sound change.  _Proceedings of the National Academy of Sciences_ , 110(11):4224–4229, 2013.  ISSN 0027-8424.  doi: 10.1073/pnas.1204678110.
  * Bouckaert et al. [2012] Remco Bouckaert, Philippe Lemey, Michael Dunn, Simon J. Greenhill, Alexander V. Alekseyenko, Alexei J. Drummond, Russell D. Gray, Marc A. Suchard, and Quentin D. Atkinson.  Mapping the origins and expansion of the indo-european language family.  _Science_ , 337(6097):957–960, 2012.  ISSN 0036-8075.  doi: 10.1126/science.1219669.
  * Bowern and Atkinson [2012] Claire Bowern and Quentin Atkinson.  Computational phylogenetics and the internal structure of pama-nyungan.  _Language_ , pages 817–845, 2012.
  * Campbell [2013] Lyle Campbell.  _Historical linguistics_.  Edinburgh University Press, 2013.
  * Chang et al. [2015] Will Chang, David Hall, Chundra Cathcart, and Andrew Garrett.  Ancestry-constrained phylogenetic analysis supports the indo-european steppe hypothesis.  _Language_ , pages 194–244, 2015.
  * De Saussure [1989] Ferdinand De Saussure.  _Cours de linguistique générale_ , volume 1.  Otto Harrassowitz Verlag, 1989.
  * Dekker and Zuidema [2020] Peter Dekker and Willem Zuidema.  Word prediction in computational historical linguistics.  _Journal of Language Modelling_ , 8(2):295–336, 2020.
  * Del Moral et al. [2006] Pierre Del Moral, Arnaud Doucet, and Ajay Jasra.  Sequential monte carlo samplers.  _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ , 68(3):411–436, 2006.
  * Felsenstein [1981] Joseph Felsenstein.  _Inferring phylogenies_.  Sinauer associates Sunderland, MA, 1981.
  * Fourrier and Sagot [2022] Clémentine Fourrier and Benoît Sagot.  Probing multilingual cognate prediction models.  In _Findings of the Association for Computational Linguistics: ACL 2022_ , 2022.
  * Gray et al. [2009] R. D. Gray, A. J. Drummond, and S. J. Greenhill.  Language phylogenies reveal expansion pulses and pauses in pacific settlement.  _Science_ , 323(5913):479–483, 2009.  ISSN 0036-8075.  doi: 10.1126/science.1166858.
  * Gray and Atkinson [2003] Russell Gray and Quentin Atkinson.  Language-Tree Divergence Times Support the Anatolian Theory of Indo-European Origin.  _Journal of the Royal Statistical Society, series B_ , 426:435–9, 12 2003.
  * Greenhill et al. [2008] Simon J Greenhill, Robert Blust, and Russell D Gray.  The austronesian basic vocabulary database: from bioinformatics to lexomics.  _Evolutionary Bioinformatics_ , 4:EBO–S893, 2008.
  * Hastings [1970] W. Keith Hastings.  Monte Carlo sampling methods using Markov chains and their applications.  _Biometrika_ , 57(1):97–109, 1970.
  * Hock [1991] Hans Henrich Hock.  _Principles of Historical Linguistics_.  De Gruyter, 2 edition, 1991.  ISBN 0899258514.  doi: 10.1515/9783110219135.
  * Hoffmann et al. [2021] Konstantin Hoffmann, Remco Bouckaert, Simon J Greenhill, and Denise Kühnert.  Bayesian phylogenetic analysis of linguistic data using BEAST.  _Journal of Language Evolution_ , 6(2):119–135, 09 2021.  ISSN 2058-458X.  doi: 10.1093/jole/lzab005.  URL https://doi.org/10.1093/jole/lzab005.
  * Kolipakam et al. [2018] Vishnupriya Kolipakam, Fiona M. Jordan, Michael Dunn, Simon J. Greenhill, Remco Bouckaert, Russell D. Gray, and Annemarie Verkerk.  A bayesian phylogenetic study of the dravidian language family.  _Royal Society Open Science_ , 5(3):171504, 2018\.  doi: 10.1098/rsos.171504.
  * Lewis [2001] Paul Lewis.  A likelihood approach to estimating phylogeny from discrete morphological character data.  _Systematic biology_ , 50:913–25, 11 2001.  doi: 10.1080/106351501753462876.
  * List et al. [2017] Johann-Mattis List, Simon J. Greenhill, and Russell D. Gray.  The potential of automatic word comparison for historical linguistics.  _PLOS ONE_ , 12(1):1–18, 01 2017.  doi: 10.1371/journal.pone.0170046.
  * Metropolis et al. [1953] Nicholas Metropolis, Arianna W. Rosenbluth, Marshall N. Rosenbluth, Augusta H. Teller, and Edward Teller.  Equation of state calculations by fast computing machines.  _The Journal of Chemical Physics_ , 21(6):1087–1092, 1953.  doi: 10.1063/1.1699114.
  * Nicholls and Gray [2008] Geoff K. Nicholls and Russell D. Gray.  Dated ancestral trees from binary trait data and their application to the diversification of languages.  _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ , 70(3):545–566, 2008.
  * Pagel [1994] Mark Pagel.  Detecting correlated evolution on phylogenies: a general method for the comparative analysis of discrete characters.  _Proceedings of the Royal Society of London. Series B: Biological Sciences_ , 255(1342):37–45, 1994.
  * Rama [2016] Taraka Rama.  Ancestry sampling for indo-european phylogeny and dates.  2016\.
  * Ryder [2010] Robin J. Ryder.  _Phylogenetic Models of Language Diversification_.  DPhil thesis, University of Oxford, 2010.
  * Sagart et al. [2019] Laurent Sagart, Guillaume Jacques, Yunfan Lai, Robin J. Ryder, Valentin Thouzeau, Simon J. Greenhill, and Johann-Mattis List.  Dated language phylogenies shed light on the ancestry of sino-tibetan.  _Proceedings of the National Academy of Sciences_ , 116(21):10317–10322, 2019.  ISSN 0027-8424.  doi: 10.1073/pnas.1817972116.

##  Part I Supplementary Materials

####  6.1 Metropolis-Hastings steps

For the mutation steps, we run many iterations of a simple Metropolis-Hastings
within Gibbs sampler. That is we update one by one each coordinate of the
parameter and latent variables conditionally on the value of the data and the
other coordinates. More precisely we adopt the following update schemes, the
probability to choose each of these update is chosen by the user so that the
most important and most difficult updates are the most frequent. More
precisely, we favour the update of the latent variables $R$ and $B$, and the
topology $g$.

We write $L$ the full likelihood computed as
$\sum_{m}\pi_{0}U_{root}(m,\cdot)$, where $U_{root}$ is computed through
pruning as described in the main text.

#####  6.1.1 Update of $R$

We sample the value of the transformations from the prior. We also tried to
add or remove one by one the transformations but that lead to poor acceptance
rate, this is due to the high correlation of the transformations, especially
if all the transformations between two values do not exist. The prior on the
transformations is induced by $k$, $p$ and $\lambda$.

The conditional writes

$\pi(R^{k}\mid B,\lambda,\mu,p,\beta,\nu,D)\propto
L(D\mid\theta)\frac{(\lambda h_{e})^{|R|}}{|R|!}p^{R},$

where $p^{R}=\prod_{r\in R}p[r]$ and $|R|=\sum_{e}|R_{e}|$.

The density of the kernel cancels out the term $\frac{(\lambda
h_{e})^{|R|}}{|R|!}p^{R}$ in the Metropolis acceptance ratio.

#####  6.1.2 Update of $B$

We choose an edge and we change the whole vector of apparition for each
cognate, we propose a new value of $B$ from the prior and accept it
independently. This may seem too simple, but it provides better results than
the addition and removal of single transformation in $B$ — as the
transformations are highly correlated.

As previously, the proposal density cancels out a part of the acceptance
ratio, only leaving the likelihood ratio.

Noticeably for this latent variable, we did not store the position of each
apparition but only the total number and the position of the last, which makes
the update arguably more difficult in the SPR moves.

#####  6.1.3 Update of the parameters $\lambda$, $\mu$, $p$, $\beta$ and $\nu$

We chose, when possible, conjugated priors — Poisson and multinomial
distributions. This makes the update of the parameters quite simple. The other
parameters are updated according to a Metropolis-Hastings step.

For $\lambda_{k}$, the update is exact, from the conjugated posterior
$\Gamma(\alpha+n,\beta+l)$, where $n$ is the number of transformations on the
whole tree for character $k$, and $l$ the length of the tree.

For $p$ the update is also exact, as the conditional is merely
$\mathcal{D}ir(\alpha+\hat{R})$, where $\hat{R}$ is the vector of counts of
the nature of the transformations.

For the other parameters, we rely on a Metropolis-Hastings step, with a
proposal from the prior for $\beta$, and Gaussian for $\nu$ and $\mu$. An
exact sampling could be use for $\mu$ but we tried non conjugate priors as
well.

#####  6.1.4 Update of $t$

Concerning the update of the length of the edges, we use a simple MH step,
where we propose to move a node in the possible range of time constrained by
the parent and children nodes. We propose also to rescale the tree, or part of
the tree. If the root age is fixed, some of those moves are impossible.

The conditional writes:

$\pi(t\mid D,B,R,\lambda,\rho,\nu)\propto L(D\mid
D,B,R,t,\nu)\prod_{k,e}\frac{(\lambda_{k}h_{e})^{|R^{k}|}}{|R^{k}|!}\frac{(\mu
h_{e})^{|B|}}{|B|!},$

the proposal is symmetric which cancels the jacobian.

#####  6.1.5 Update of the topology

This last move is the most important and consequently the most difficult to
choose. In the previous works [Gray and Atkinson, 2003, Ryder, 2010] a new
tree is built by an SPR (sample-prune-regraft) move, in practice a subtree is
always regraft at the nearest possible position — its cousin, meaning that the
(sub)-tree $((a,b),c)$ becomes $(a,(b,c))$.

The question that remains is how to transport the latent variables of the
changed parts of the tree. In order to preserve the intuition behind the
latent variables, it seems reasonable to either keep all the ages of the nodes
constant, or to change the ages while only moving to the nearest possible
position — it is not possible to keep both as it would be incompatible most of
the time. We describe in algorithm 1 the main method we used, a graphical
representation can be found in Figure 8, the case of change with the nearest
neighbour corresponds to $i=h$.

Let $g$ be a tree, $X$ the edges dependent latent variables.

Choose $f\in g$ a non root node, let $e$ be its parent, $h$ its grandparent
and $f^{\prime}$ its brother.

Choose $j$ a node, whose parent is denoted $i$, such that $t(e)\in[t(j),t(i)]$
where $t$ is the height of node.

Set $j$ as the new brother of $f$, by splitting the edge $\langle i,j\rangle$
into $\langle i,e\rangle$ and $\langle e,f\rangle$, such that $t(e)$ remains
constant. Merge the edges $\langle h,e\rangle$ and $\langle
e,f^{\prime}\rangle$.

Merge the latent variables on $\langle h,e\rangle$ and $\langle
e,f^{\prime}\rangle$ into latent variables on $\langle h,f^{\prime}\rangle$.

Split the latent variables on $\langle i,j\rangle$ into latent variables on
$\langle i,e\rangle$ and $\langle e,j\rangle$ proportionally to the length of
the edges.

Sample from the prior latent variables on $\langle e,f\rangle$.

Algorithm 1 Sample Prune and Regraft algorithm with latent variables. See
Figure 8.
$h$$e$$f^{\prime}$$f$$i$$j$$k$$h$$e$$f^{\prime}$$f$$i$$j$$k$$\dots$$\longrightarrow$$\dots$
Figure 8: Move of the subtree starting at $f$ so that $f$ has $j$ as brother.
The colors represent the edge-attached latent variables.

####  6.2 Further results on the simulated dataset

On the side of the parameters, we provide in Figure 9 the posterior
distributions for the main evolution parameters. They seem quite efficiently
reconstructed — particularly the transformation rate $\beta$ which is close to
$1$, the true value — with an underestimation overall the edge length
parameters, which is linked to the overestimation in edge length — the number
of events per branch is well reconstructed in the lower parts of the tree —.
This induces a poor performance in the reconstruction of cognacy when studying
distantly related leaves, as the number of apparitions on the long edges
starting from the root is underestimated. Maybe the number of character should
be increased so that the creation of new words becomes cheaper compared to the
evolution. In Figure 10 we represent the posterior cognacy probability. The
number of false positive is quite high. However, in general, meanings with low
cognacy probability do corresponds to words that are non cognates in truth.
The root age is overestimated in accordance to what has been said before.

Figure 9: Posterior estimations of $\lambda$, $\rho$, $\nu$, $\beta$ ant the
root age. The true values are indicated by the vertical lines. Figure 10:
Posterior cognacy probability for each meaning between $t4$ and $t6$. The
color indicates the true cognacy status.

Concerning the reconstructed ancestral state, the result highly depends on the
age of the common ancestor. Inferring the state of the common ancestor to $9$
and $10$ is fairly simple. More interestingly, we present in Figure 11 the
posterior probability for the values of characters $1$ and $4$ at the common
ancestor between $2$ and $6$, this node is correctly reconstructed only with
$87\%$. If the majority of the values are correctly reconstructed, some are
completely off. Noticeably, the first character, that evolves faster, is more
uncertain than the fourth. We advise to use this reconstruction method with
caution, as the efficiency of the method can greatly vary between meanings and
characters.

Figure 11: Posterior distribution of internal ancestral values for the first
and fourth character at the nearest common ancestor of $6$ and $2$. The color
represents the posterior probability for each value, the true value
corresponds to the black dot.

#####  6.2.1 Other prior

We tried a more classical prior on the tree, with the same specifications as
previously. Leading to the results presented in Figure 12

Figure 12: Consensus tree and root age posterior distribution for the Gamma
prior.

####  6.3 Misspecifications of the model

Our model is clearly an oversimplification of the reality. We propose here to
test the resilience of the model to some of the most common

#####  6.3.1 Correlated traits

We ran the experiment on the same tree, with the same characters and values —
although with a newly generated dataset, with the particularity that the
values for the second character do not come from the above described dynamic.
At the end of each edge, the value for the second character is the same as the
first one with some noise.

The resulting consensus trees are presented in Figure 13. Overall the
inference does not seem very much affected, the tree posterior is strikingly
concentrated around the consensus topology. The parameters associated with the
second character seem a little bit off, which may be strange as _both_
characters are correlated.

Figure 13: Consensus tree for the correlated dataset test.

Figure 14: Posterior estimations of $\lambda$, $\rho$, $\nu$, $\beta$ ant the
root age. The true values are indicated by the vertical lines.

#####  6.3.2 Selection of transformations

If the set of possible transformations in the implementation is larger than
the true one, we can assume that the inference will lose in efficiency. In the
previous model, it is clear that the likelihood of the model is drastically
changed if the possible transformations are changed. The identifiability of
the model is even unclear, as there are several parameters $p$ that can lead
to a same likelihood.

We ran our numerical methods on such a model, the parameter of the experiment
are the same as the synthetic exact experiment. The only change is that we did
not reduce the number of transitions possible for the two first characters (we
have 20 possible transitions instead of the 12 truly present).

Figure 15: True and consensus tree for the unknown transformations sets

However, this is clearly not efficient if the size of $\mathfrak{R}_{k}$ is
too high. For example, if we allow all the pairwise transformations in a
character with $30$ value we get more than a thousand pairs ; the convergence
of the numerical methods are then vastly slower and we cannot reach a
satisfying equilibrium, and the data for this character can only be explained
with noise or apparitions, which is not optimal. We recommend to reduce to
less than 30 the size of $\mathfrak{R}_{k}$.

In these experiments the inference of the other parameter is of lesser
quality, with parameter overall underestimated, see Figure 16. The inference
of the probability of the transformations is not quite satisfying, as we can
see in Figure 17, there is no significant difference in the reconstituted
probability, this comes certainly from the small number of transformations in
total on the tree — around a dozen.

Figure 16: Posterior estimations of $\lambda$, $\rho$, $\nu$, $\beta$ ant the
root age for the experiments with unknown $\mathfrak{R}_{k}$. The true values
are indicated by the vertical lines. Figure 17: Posterior of the
transformation repartition $p$ for the first characters the real
transformations are indicated by dots. The color represents the posterior
intensity.

#####  6.3.3 Data simulated from a forest

In this last example, maybe the most interesting for our application, we
concatenate the dataset coming from a $11$-leaved tree and a $4$-leaved tree
represented in Figure 18. That is we try to infer a single tree where there
should be two. This is a quite important case as we are never sure that some
language are indeed related.

The results are presented in Figure 19. Interestingly enough, the two group of
language are reconstructed as subtree with rather long edges. Clearly, if the
evolution parameter were higher this would not be the case, especially for a
small number of characters and a small number of possible value. The evolution
parameter however are quite off, as they must be quite important to explain
this unique tree. The root is not so old, which is an effect of the prior.

Figure 18: Two trees used for the forest example

Figure 19: Results of the inference for the topology and the apparition
parameter $\rho$ on the forest dataset.


11institutetext: NepAL Applied Mathematics and Informatics Institute for
Research (NAAMII), Kathmandu, Nepal 22institutetext: Department of Radiology,
Feinberg School of Medicine, Northwestern University, Chicago, USA
33institutetext: Institute of Biomedical Engineering, Department of
Engineering Science, University of Oxford, OX3 7DQ, Oxford, UK
44institutetext: NIHR Oxford Biomedical Research Centre, Oxford, UK

# TGANet: Text-guided attention for improved polyp segmentation

Nikhil Kumar Tomar 11 Debesh Jha 22 Ulas Bagci 22 Sharib Ali 3344

###### Abstract

Colonoscopy is a gold standard procedure but is highly operator-dependent.
Automated polyp segmentation, a precancerous precursor, can minimize missed
rates and timely treatment of colon cancer at an early stage. Even though
there are deep learning methods developed for this task, variability in polyp
size can impact model training, thereby limiting it to the size attribute of
the majority of samples in the training dataset that may provide sub-optimal
results to differently sized polyps. In this work, we exploit size-related and
polyp number-related features in the form of text attention during training.
We introduce an auxiliary classification task to weight the text-based
embedding that allows network to learn additional feature representations that
can distinctly adapt to differently sized polyps and can adapt to cases with
multiple polyps. Our experimental results demonstrate that these added text
embeddings improve the overall performance of the model compared to state-of-
the-art segmentation methods. We explore four different datasets and provide
insights for size-specific improvements. Our proposed text-guided attention
network (TGANet) can generalize well to variable-sized polyps in different
datasets. Codes are available at https://github.com/nikhilroxtomar/TGANet.

###### Keywords:

Label embedding polyp multi-scale features attention

##  1 Introduction

Colorectal cancer (CRC) is one of the leading causes of cancer-related deaths
[16] worldwide. However, high operator-dependence and subjectivity during gold
standard colonoscopic procedures remain high. This is also due to the complex
topology of organ, severe artefacts, constant deformation of organ, debris and
stool etc. Even though cleansing of the bowel is done to improve detection
rates of cancer and cancer precursor lesions such as polyps, the missed rate
is still high that accounts for 26.8% for polyps located on the right colon
and 21.4% to polyps on the left colon [10, 12]. In addition, the missed rate
for flat or sessile polyps (diminutive polyps) is grim (nearly 32.7%). An
automated system is thus clearly needed to minimize the operator subjectivity
and missed rate. Semantic segmentation can classify each pixel into a class
category, allowing the opportunity to learn meaningful semantic
representations of polyps and their complex surroundings. Several methods do
exist in the literature [3, 7, 14] but most them focus on exploiting only
localized spatial context. However, the nature and occurrence of polyps in
colonic mucosa can be confused with colonic folds. Exploiting associated
attributes such as size and occurrence (one or a few) could be used to infer
and improve segmentation for hard samples.

Encoder-decoder networks has been widely used for polyp segmentation using
various modifications to boost network performance [3, 7, 14]. PraNet [3]
applied area and boundary cues in reverse attention to focus on the polyp
boundary regions. The high-level feature aggregation and boundary attention
blocks in the network help to calibrate some of the misaligned predictions and
improve the segmentation accuracy. Similarly, HRENet [14] designed an
informative context enhancement (ICE) technique and adaptive feature
aggregation (AFA) module and trained the model on their edge and structure
consistency aware loss (ESCLoss), and obtained superior performance. Other
works such as PolypSeg [18] and MSRFNet [15] uses modules incorporating
multiple-scale information. An adaptive scale context module (ASCM) and
semantic global context module (SGCM) was used in PolypSeg [18]. The ASCM
tackles the size variations among the polyp and improves the better feature
representation capability, while SGCM enhances the feature fusion between the
high-level and low-level features and remove noise in the low-level features
to improve the segmentation accuracy. Similarly, MSRFNet [15] integrated
cross-scale fusion modules to propagate both high resolution and low-
resolution features and an added shape stream network to prune polyp
boundaries.

Most of these works in the literature [3, 14, 15] focuses on size variation,
boundary curves, background regions, dense skip connections and dense residual
scale fusions that can boost performance. However, these adjustments are made
using additional layers and explicit extensions of networks and their
connections. This adds to the complexity of the model that can adversely
affect the generalization of test samples coming from a similar distribution
and require a large dataset. In addition, it can also affect images with
underrepresented polyp sizes. In this work, we propose incorporating a text
guided attention mechanism using a simple byte-pair encoding for the
attributes comprising polyp number and its size. In addition, we use the same
encoder layer of the network to provide weights for each of these attributes.

The main contributions of the presented work include - 1) text guided
attention to learn different features in the context of the number of polyps
presence (one or many) and size (small, medium and large), 2) feature
enhancement module to strengthen the features of the encoder and pass them to
the decoder, and 3) multi-scale feature aggregation to capture features
learned by different decoder blocks. We have evaluated our TGANet on four
publicly available polyp datasets and compared it with five SOTA medical image
segmentation methods.

Figure 1: Block diagram of the proposed TGANet

##  2 Method

The proposed TGANet is a polyp segmentation architecture with text guided
attention that enables to enhance feature representations such that the polyps
present in images are segmented optimally independent of their size
variability and occurrence. Our TGANet architecture consists of various
components that are shown in Figure 1 and elaborated below.

###  2.1 Encoder module

TGANet is built upon a pre-trained ResNet50 [4] as backbone encoder network
for which we use its four different encoding blocks, $e_{i},i\in{1,2,3,4}$.
These blocks are consecutively used for our auxiliary attribute classification
task and for main polyp segmentation task. For the text-attribute
classification, we use the output from the fourth encoder block as two
classification task modules separately, i.e., number of polyps (one or many)
and their size (small, medium and large). Here, softmax probabilities
$\mathrm{\sigma}_{polyp}^{no}(.)$ and $\mathrm{\sigma}_{polyp}^{sz}(.)$ are
predicted. For the main segmentation task, we take the output from each
ResNet50 block and passes them through the feature enhancement module (FEM,
$f_{i},i\in{1,2,3,4}$) that is responsible for strengthening the features by
applying multiple dilated convolutions and an attention mechanism.

###  2.2 Feature enhancement module

Feature enhancement module (FEM) (see Figure 1 (b)) begins with four parallel
dilated convolutions $\mathrm{Conv}$ with a dilation rate $r=\\{1,6,12,18\\}$.
Each dilation is followed by a batch normalization $\mathrm{BN}$ and a
rectified linear unit $\mathrm{ReLU}$ which we refer as $\mathrm{CBR}$. The
output features are passed through a channel-attention module $\mathrm{CAM}$
[17] to capture the explicit relationship between the feature channels. The
highlighted features from these four dilated convolutions are then
concatenated and passed through a $\mathrm{Conv}_{3\times 3}$ followed by
$\mathrm{BN}$ layer and added with the original input features through a
$\mathrm{Conv}_{1\times 1}$. The resulting features are then followed by a
$\mathrm{ReLU}$ activation function, and a spatial attention mechanism
$\mathrm{SAM}$ [17] is applied to suppress the irrelevant regions.

###  2.3 Label attention

Label attention module is designed to provide learned text-based attention to
the output features of the decoder blocks in our TGANet. Here, we use three
label attention modules, $l_{i},i\in{1,2,3}$, as soft channel-wise attention
to the three decoder outputs that enables larger weights to the representative
features and suppress the redundant ones. The first label attention module
uses the output of the embedding fusion $\mathcal{E}(.)$ obtained by element-
wise dot product between the softmax probability concatenation
$\\{\mathrm{\sigma}_{one},\mathrm{\sigma}_{many},\mathrm{\sigma}_{small},\mathrm{\sigma}_{medium},\mathrm{\sigma}_{large}\\}$
with the encoded text embedding. Say,
$\mathcal{A}=\\{{\mathrm{one},\mathrm{many},\mathrm{small},\mathrm{medium},\mathrm{large}}\\}$
be the attributes that are encoded using byte-pair encoding [5] and denoted by
$\mathcal{A}_{encode}$ with $\\{a_{j}^{k}\\}$ as vector embedding for each
attribute $j$ of length $|k|$, then $\mathcal{E}(.)$ that is given by:

$\mathcal{E}=\sigma_{j}\odot{a_{j}^{k}},\quad\forall{k}.$ (1)

The output of the label attention module is referred to as label features
$l_{f}$ in this paper.

###  2.4 Decoder

The decoder in the proposed TGANet is comprised of three different decoder
blocks $d_{i},i\in{1,2,3}$, of which each takes the input features to upsample
it and pass it through some convolutional layers to produce the output. This
output is refined using the label attention module $l_{i}$ and passed to the
subsequent decoder blocks $d_{i}$ (see Figure 1 (c)). The first decoder block
takes the output of the fourth FEM $f_{4}$ to upsample it using bilinear
interpolation by a factor of two, and then it is concatenated with the output
features from the third FEM $f_{3}$. The resulting concatenated feature is
passed through a $\mathrm{Conv}_{1\times 1}$-$\mathrm{BN}$-$\mathrm{ReLU}$
referred as $\mathrm{CBR}$ followed by a sequence of three
$\mathrm{Conv}_{3\times 3}$-$\mathrm{BN}$, further accompanied by their
multiple residual connections and a $\mathrm{ReLU}$ activation function with
subsequent convolutional block attention module represented as $d_{i}^{cbam}$.
An element-wise multiplication is done to allow additional soft-attention from
the computed label features $l_{f}$ using a sigmoid function for each decoder
block output $d_{i}^{out},i\in{1,2,3}$ given by:

$d_{i}^{out}=d_{i}^{cbam}\odot\sigma([\mathrm{Conv}-\mathrm{ReLU}-\mathrm{Conv}]l_{f}),\quad\forall{i\in{1,2,3}}$
(2)

###  2.5 Multi-scale feature aggregation

Multi-scale feature aggregation (MSFA) module (see supplementary Figure 1) is
used to fuse multi-scale feature representations at various decoder outputs
$d_{i}^{out},i\in{1,2,3}$ that allows to capture learned features. We take the
first two features $\\{d_{1}^{out},d_{2}^{out}\\}$ and pass them through a
bilinear upsampling to ensure that all three features have the exact spatial
dimensions followed by linear $1\times 1$ convolution layers, $\mathrm{BN}$
and $\mathrm{ReLU}$ activation before concatenation. To boost the capture of
non-linear features we further apply a series of convolutional layers,
$\mathrm{BN}$ and $\mathrm{ReLU}$ together with multiple residual connections
for improved flow of information. The output is represented as $m_{f}$ which
is responsible for our predicted segmentation map $I^{pred}_{mask}$ given by:
$I^{pred}_{mask}=\sigma(\mathrm{Conv}_{1\times 1}(m_{f}))$.

Joint loss optimization: We jointly minimize loss for both the auxiliary
classification tasks (cross-entropy losses, $\mathrm{CE}_{loss1}$,
$\mathrm{CE}_{loss2}$) and the segmentation task (binary cross entropy,
$\mathrm{BCE}_{loss3}$ and dice loss, $\mathrm{DSC}_{loss4}$) with equal
weights.

##  3 Experiments and results

Table 1: Quantitative results on the experimented polyp datasets. Method | Backbone | mIoU | mDSC | Recall | Precision | F2
---|---|---|---|---|---|---
Dataset: Kvasir-SEG [9]
U-Net [13] | - | 0.7472 | 0.8264 | 0.8504 | 0.8703 | 0.8353
HarDNet-MSEG [6] | HardNet68 | 0.7459 | 0.8260 | 0.8485 | 0.8652 | 0.8358
ColonSegNet [7] | - | 0.6980 | 0.7920 | 0.8193 | 0.8432 | 0.7999
DeepLabV3+ [2] | ResNet50 | 0.8172 | 0.8837 | 0.9014 | 0.9028 | 0.8904
PraNet [3] | Res2Net | 0.8296 | 0.8942 | 0.9060 | 0.9126 | 0.8976
TGANet (Ours) | ResNet50 | 0.8330 | 0.8982 | 0.9132 | 0.9123 | 0.9029
Dataset: CVC-ClinicDB [1]
U-Net [13] | - | 0.8428 | 0.8978 | 0.9001 | 0.9209 | 0.8981
HarDNet-MSEG [6] | HardNet68 | 0.8388 | 0.8967 | 0.8929 | 0.9216 | 0.8938
ColonSegNet [7] | - | 0.8248 | 0.8862 | 0.8828 | 0.9017 | 0.8826
DeepLabV3+ [6] | ResNet50 | 0.8973 | 0.9391 | 0.9441 | 0.9442 | 0.9389
PraNet [3] | Res2Net | 0.8866 | 0.9318 | 0.9347 | 0.9479 | 0.9333
TGANet (Ours) | ResNet50 | 0.8990 | 0.9457 | 0.9437 | 0.9519 | 0.9439
Dataset: BKAI [11]
U-Net [13] | - | 0.7599 | 0.8286 | 0.8295 | 0.8999 | 0.8264
HarDNet-MSEG | HardNet68 | 0.6734 | 0.7627 | 0.7532 | 0.8344 | 0.7528
ColonSegNet [7] | - | 0.6881 | 0.7748 | 0.7852 | 0.8711 | 0.7746
DeepLabV3+ [2] | ResNet50 | 0.8314 | 0.8937 | 0.8870 | 0.9333 | 0.8882
PraNet [3] | Res2Net | 0.8264 | 0.8904 | 0.8901 | 0.9247 | 0.8885
TGANet (Ours) | ResNet50 | 0.8409 | 0.9023 | 0.9026 | 0.9208 | 0.9002
Dataset: Kvasir-Sessile [8]
U-Net [13] | - | 0.2472 | 0.3688 | 0.7237 | 0.3264 | 0.4635
HarDNet-MSEG | HardNet68 | 0.1565 | 0.2558 | 0.5403 | 0.2236 | 0.3298
ColonSegNet [7] | - | 0.2113 | 0.3278 | 0.5234 | 0.3336 | 0.3868
DeepLabV3+ [2] | ResNet50 | 0.5927 | 0.7078 | 0.7085 | 0.8225 | 0.7009
PraNet [3] | Res2Net | 0.6671 | 0.7736 | 0.8069 | 0.8244 | 0.7871
TGANet (Ours) | ResNet50 | 0.6910 | 0.7980 | 0.7925 | 0.8588 | 0.7879
Training dataset: Kvasir-SEG – Test dataset: CVC-ClinicDB
U-Net [13] | - | 0.5433 | 0.6336 | 0.6982 | 0.7891 | 0.6563
HarDNet-MSEG [6] | HardNet68 | 0.6058 | 0.6960 | 0.7173 | 0.8528 | 0.7010
ColonSegNet [7] | - | 0.5090 | 0.6126 | 0.6564 | 0.7521 | 0.6246
DeepLabV3+ [2] | ResNet50 | 0.7388 | 0.8142 | 0.8331 | 0.8735 | 0.8198
PraNet [3] | Res2Net | 0.7286 | 0.8046 | 0.8188 | 0.8968 | 0.8077
TGANet (Ours) | ResNet50 | 0.7444 | 0.8196 | 0.8290 | 0.8879 | 0.8207

###  3.1 Datasets

To evaluate the performance of our TGANet, we have used four publicly
available polyp segmentation benchmark datasets including Kvasir-SEG [9], CVC-
ClinicDB [1], BKAI[11], and Kvasir-Sessile [8] (details are presented in
supplementary Table 1). Relevant to our experiment, Kvasir-Sessile [8]
contains 196 small, diminutive, sessile and flat polyps that are less than 10
mm in size.

###  3.2 Implementation Details

All models are trained on NVIDIA GeForce RTX 3090 GPU, and images are resized
to $256\times 256$ pixels with 80:10:10 training, validation, and testing
splits except for Kvasir-SEG, where we adopted the official split of 880/120
for training and testing. Simple data augmentation strategy including random
rotation, vertical flipping, horizontal flipping, and coarse dropout are used.
All models are trained on a similar hyperparameters configuration with a
learning rate of $1e^{-4}$, batch size of 16, and optimized with Adam
optimizer. An early stopping mechanism and ReduceLROnPlateau is used to
prevent models from overfitting.

Standard medical image segmentation metrics such as mean intersection over
union (mIoU), mean Sørensen–dice coefficient (mDSC), recall, precision,
F2-score and frame per second (FPS) are used.

###  3.3 Results

Table 2: mDSC for different sizes and polyp counts on Kvasir-SEG

[9] Method small medium large one many DeepLabV3+ [2] 0.8776 0.9003 0.8633
0.8922 0.8289 PraNet [3] 0.8826 0.9079 0.8900 0.9071 0.8106 TGANet 0.8869
0.9203 0.8769 0.9075 0.8378

Figure 2: Qualitative results comparison on the Kvasir-SEG dataset.

We have compared our results with five SOTA methods that include UNet [13],
HarDNet-MSEG [6], ColonSegNet [7], DeepLabv3+ [2], and PraNet [3]. These
algorithms are widely used baselines in both polyp segmentation and general
medical image segmentation. The quantitative results are presented in Table 1.

#### Results on Kvasir-SEG:

Table 1 shows that TGANet outperforms all the SOTA methods with a mIoU of
0.8330 and mDSC of 0.8982. Our TGANet outperforms most competitive PraNet [3]
by 1.58% in mIoU and 1.45% mDSC.

#### Results on CVC-ClinicDB:

For CVC-ClinicDB dataset, TGANet outperforms all SOTA methods reporting the
highest mIoU and mDSC of 0.8990 and 0.9457, respectively. Our method
outperformed the most competitive DeepLabV3+ [2] with a mIoU of 0.17% and mDSC
of 0.66%.

#### Results on BKAI:

Table 1 shows the comparison of the result on the BKAI dataset that show that
our proposed TGANet obtains mIoU of 0.8409 and mDSC of 0.9023 and outperforms
the best performing DeepLabV3+ [2] by 0.95% on mIoU and 0.86% on mDSC.

#### Results on Kvasir-Sessile:

Kvasir-Sessile dataset is clinically most relevant as it has flat and sessile
polyps. On this dataset, it can be observed (see Table 1) that our TGANet
surpasses all the other methods in all the evaluation metrics. It outperforms
the best performing PraNet [3] by a large margin of 2.39% on mIoU and 2.44% on
mDSC. Similarly, almost 10% improvement is observed compared to the DeepLabV3+
[2] in this case which is a significant improvement.

#### Results on cross dataset:

To explore the generalization capability of our proposed TGSNet, we train the
model on Kvasir-SEG and test it on the CVC-ClinicDB. The cross-dataset test
(Table 1) also suggested improvements compared to all SOTA methods and
obtained an increment of 0.56% on mIoU and 0.54% on mDSC compared to the SOTA
DeepLabv3+ [2].

#### Results on size and number-based sampled polyps:

To show the effectiveness of our proposed TGANet, we evaluated test samples of
Kvasir-SEG-based on the attributes used in training. It can be observed in
Table 2 that our model outperforms the best SOTA methods for almost all cases.
For the ‘small’, ‘medium’ and ‘many cases’, the improvement ranges from nearly
1-2%.

Our qualitative results (see Figure 2) demonstrate a clear improvement of our
text-based attention method for different sizes and number polyp samples. It
is evident that both PraNet [3] and DeepLabV3+ [2] failed to capture sample
with two polyps (4th row) and also provided over segmentation for the small
(1st row) and medium polyps (2nd row). Additionally, we have provided the
total number of parameters, flops and FPS in supplementary material Table 2.

###  3.4 Ablation study

Table 3: Ablation study of TGANet on Kvasir-SEG No | Method | mIoU | mDSC | Recall | Precision | F2
---|---|---|---|---|---|---
#1 | TGANet w/o label and classifier | 0.8104 | 0.8786 | 0.8987 | 0.8970 | 0.8850
#2 | TGANet w/o MSFA | 0.8151 | 0.8832 | 0.9061 | 0.8999 | 0.8907
#3 | TGANet w/o FEM | 0.8084 | 0.8766 | 0.8968 | 0.9010 | 0.8838
#4 |  TGANet w/o (label+classifier+ MSFA+FEM) | 0.8063 | 0.8747 | 0.8963 | 0.8971 | 0.8798
#5 | TGANet (Ours) | 0.8330 | 0.8982 | 0.9132 | 0.9123 | 0.9029

To validate the effectiveness and importance of the core components used in
the network, we compare TGANet with its five variants, which is presented in
Table 3. The results suggest that the introduction of the text guided
attention along with the label boosts the performance of the network. The
results show that TGANet improves the baseline without the label and
classifier (#1) by 2.26% on mIoU and 1.96% on mDSC.

##  4 Conclusion

We proposed a text-guided attention architecture (TGANet) to tackle polyps’
variable size and number for robust polyp segmentation. We have used multiple
feature enhancement modules connected with different encoder blocks to achieve
this. An auxiliary task is learned together with the main task to compliment
both the size-based and number-based feature representations and used as label
attentions in the decoder blocks. Additionally, the multi-scale fusion of the
features at the decoder enabled our network to deal with these attribute
changes. Our experimental results demonstrated the effectiveness of our TGANet
outperformed and provided higher segmentation performance on flat and sessile
polyps that are clinically important.

## References

  * [1] Bernal, J., Sánchez, F.J., Fernández-Esparrach, G., Gil, D., Rodríguez, C., Vilariño, F.: WM-DOVA maps for accurate polyp highlighting in colonoscopy: Validation vs. saliency maps from physicians. Computerized Medical Imaging and Graphics 43, 99–111 (2015)
  * [2] Chen, L.C., Zhu, Y., Papandreou, G., Schroff, F., Adam, H.: Encoder-decoder with atrous separable convolution for semantic image segmentation. In: Proceedings of the European conference on computer vision (ECCV). pp. 801–818 (2018)
  * [3] Fan, D.P., Ji, G.P., Zhou, T., Chen, G., Fu, H., Shen, J., Shao, L.: Pranet: Parallel reverse attention network for polyp segmentation. In: Proceedings of the International conference on medical image computing and computer-assisted intervention (MICCAI). pp. 263–273 (2020)
  * [4] He, K., Zhang, X., Ren, S., Sun, J.: Deep residual learning for image recognition. In: Proceedings of the IEEE conference on computer vision and pattern recognition (CVPR). pp. 770–778 (2016)
  * [5] Heinzerling, B., Strube, M.: BPEmb: Tokenization-free Pre-trained Subword Embeddings in 275 Languages. In: Proceedings of the International Conference on Language Resources and Evaluation (LREC 2018) (2018)
  * [6] Huang, C.H., Wu, H.Y., Lin, Y.L.: HarDNet-MSEG A Simple Encoder-Decoder Polyp Segmentation Neural Network that Achieves over 0.9 Mean Dice and 86 FPS. arXiv preprint arXiv:2101.07172 (2021)
  * [7] Jha, D., Ali, S., Tomar, N.K., Johansen, H.D., Johansen, D., Rittscher, J., Riegler, M.A., Halvorsen, P.: Real-time polyp detection, localization and segmentation in colonoscopy using deep learning. IEEE Access 9, 40496–40510 (2021)
  * [8] Jha, D., Smedsrud, P.H., Johansen, D., de Lange, T., Johansen, H.D., Halvorsen, P., Riegler, M.A.: A comprehensive study on Colorectal Polyp Segmentation With ResUNet++, Conditional Random Field and Test-Time Augmentation. IEEE journal of Biomedical and Health Informatics 25(6), 2029–2040 (2021)
  * [9] Jha, D., Smedsrud, P.H., Riegler, M.A., Halvorsen, P., Lange, T.d., Johansen, D., Johansen, H.D.: Kvasir-SEG: a segmented polyp dataset. In: Proceedings of the International Conference on Multimedia Modeling (MMM). pp. 451–462 (2020)
  * [10] Kim, N.H., Jung, Y.S., Jeong, W.S., Yang, H.J., Park, S.K., Choi, K., Park, D.I.: Miss rate of colorectal neoplastic polyps and risk factors for missed polyps in consecutive colonoscopies. Intestinal research 15(3), 411 (2017)
  * [11] Lan, P.N., An, N.S., Hang, D.V., Van Long, D., Trung, T.Q., Thuy, N.T., Sang, D.V.: NeoUNet: towards accurate colon polyp segmentation and neoplasm detection. arXiv preprint arXiv:2107.05023 (2021)
  * [12] Rex, D.K., Cutler, C.S., Lemmel, G.T., Rahmani, E.Y., Clark, D.W., Helper, D.J., Lehman, G.A., Mark, D.G.: Colonoscopic miss rates of adenomas determined by back-to-back colonoscopies. Gastroenterology 112(1), 24–28 (1997)
  * [13] Ronneberger, O., Fischer, P., Brox, T.: U-Net: Convolutional Networks for Biomedical Image Segmentation. In: Proceedings of the International Conference on Medical image computing and computer-assisted intervention (MICCAI). pp. 234–241 (2015)
  * [14] Shen, Y., Jia, X., Meng, M.Q.H.: Hrenet: A hard region enhancement network for polyp segmentation. In: Proceedings of the International Conference on Medical Image Computing and Computer-Assisted Intervention (MICCAI). pp. 559–568 (2021)
  * [15] Srivastava, A., Jha, D., Chanda, S., Pal, U., Johansen, H.D., Johansen, D., Riegler, M.A., Ali, S., Halvorsen, P.: MSRF-Net: A Multi-scale Residual Fusion Network for Biomedical Image Segmentation. IEEE Journal of biomedical imaging and health informatics (2021)
  * [16] Sung, H., Ferlay, J., Siegel, R.L., Laversanne, M., Soerjomataram, I., Jemal, A., Bray, F.: Global cancer statistics 2020: Globocan estimates of incidence and mortality worldwide for 36 cancers in 185 countries. CA: a cancer journal for clinicians 71(3), 209–249 (2021)
  * [17] Woo, S., Park, J., Lee, J.Y., Kweon, I.S.: Cbam: Convolutional block attention module. In: Proceedings of the European conference on computer vision (ECCV). pp. 3–19 (2018)
  * [18] Zhong, J., Wang, W., Wu, H., Wen, Z., Qin, J.: Polypseg: An efficient context-aware network for polyp segmentation from colonoscopy videos. In: Proceedings of the International Conference on Medical Image Computing and Computer-Assisted Intervention (MICCAI). pp. 285–294 (2020)

## Supplementary material

List of suppfigures 1 Block diagram of the proposed multi-scale fusion aggregation module that takes different scale features from three decoder blocks and concatenates them using upsampling, convolution, batch normalization and $\mathrm{ReLU}$ activations. Also, we apply residual-connections for different layers as shown. List of suppTables 1 Polyp datasets used in our experiments with number of images, input size and their availability. For some indicated dataset test data are not available (e.g. BKAI) for which we used training data and split it into 80:10:10 for train, validation and test. Also, we use small, flat and sessile polyps in Kvasir-SEG. Dataset |  Images | Input size | Availability
---|---|---|---
Kvasir-SEG | 1000 | Variable | Public
CVC-ClinicDB | 612 | 384 × 288 | Public
BKAI | 1000† | Variable | Public
Kvasir-Sessile⋄ | 196 | Variable | Public
† test data unavailable ⋄sessile polyps from Kvasir-SEG
List of suppTables 2 Algorithm complexity of methods used in the study. Number of model parameters, flops, image size and frame-per-second is provided for both state-of-the-art (SOTA) methods and proposed TGANet. Method | Parameters (Million) | Flops (GMac) | Image Size | FPS
---|---|---|---|---
U-Net | 31.04 | 54.75 | $256\times 256$ | 164.93
HarDNet-MSEG | 33.34 | 6.02 | $256\times 256$ | 42.95
ColonSegNet | 5.01 | 62.16 | $256\times 256$ | 127.95
DeepLabV3+ | 39.76 | 43.31 | $256\times 256$ | 99.79
PraNet | 32.55 | 6.93 | $256\times 256$ | 36.90
TGANet (Ours) | 19.84 | 41.88 | $256\times 256$ | 37.64