# Computing algebraic degrees
of phylogenetic varieties

Luis David Garcia Puente<EMAIL_ADDRESS>Department of
Mathematics and Computer Science, Colorado College, Colorado Springs, CO, USA
Marina Garrote-López<EMAIL_ADDRESS>Department of Mathematics,
University of British Columbia, Vancouver, Canada Elima Shehu elima.shehu@uni-
osnabrueck.de Osnabrück University, Osnabrück, Germany

(October 05, 2022)

###### Abstract

A phylogenetic variety is an algebraic variety parameterized by a statistical
model of the evolution of biological sequences along a tree. Understanding
this variety is an important problem in the area of algebraic statistics with
applications in phylogeny reconstruction. In the broader area of algebra
statistics, there have been important theoretical advances in computing
certain invariants associated to algebraic varieties arising in applications.
Beyond the dimension and degree of a variety, one is interested in computing
other algebraic degrees, such as the maximum likelihood degree and the
Euclidean distance degree. Despite these efforts, the current literature lacks
explicit computations of these invariants for the particular case of
phylogenetic varieties. In our work, we fill this gap by computing these
invariants for phylogenetic varieties arising from the simplest group-based
models of nucleotide substitution Cavender-Farris-Neyman model, Jukes-Cantor
model, Kimura 2-parameter model, and the Kimura 3-parameter model on small
phylogenetic trees with at most 5 leaves.

##  1 Introduction

One of the main goals of phylogenetics, is to reconstruct the evolutionary
history of a group of different biological entities, such as genes, bacterial
strains or species. In phylogenetics, we model these evolutionary
relationships with a phylogenetic tree, where the leaves represent the extant
entities of interest and the internal nodes represent their extinct common
ancestors. After some mild hypothesis, to model the substitution of
nucleotides in an evolutionary process we adopt a parametric statistical
model, that is, we model the substitution of nucleotides along a phylogenetic
tree as a Markov process on the tree. This allows to define a parameterization
map that describe the joint distributions at the leaves of the tree in terms
of the parameters of the model.

The resulting polynomial map
$\phi^{\mathcal{M}}_{T}:\mathbb{R}^{d}\rightarrow\mathbb{R}^{m}$ associated to
a tree $T$ and a model $\mathcal{M}$ maps any $d$-tuple of parameters to a
distribution vector of the $m$ possible observations at the leaves. The
Zariski closure of the image of this map, denoted
$\overline{Im(\phi^{\mathcal{M}}_{T}})$, is the _phylogenetic variety_
associated to $T$ and $\mathcal{M}$. For certain evolutionary models this
variety is a projective toric variety.

In applications, one is interested in finding the phylogenetic tree and the
parameters that best explains some given evolutionary data. Geometrically,
this can be interpreted as, given data point, computing the closest point on
the varieties associated to the different possible trees. It is known that,
for complex varieties, the number of complex non-singular critical points of
the squared Euclidean distance to a generic point outside the variety is
constant in a dense open subset, [DHO+16]. This number, known as the
_Euclidean distance degree_ (ED degree) of the variety, measures the algebraic
complexity of the variety providing an upper bound for the number of real
solutions. Another approach to recover the tree and parameters that best
explains some data is by maximizing the likelihood of the observed
distribution. That is, the maximum likelihood estimation is a similar
optimization problem where one tries to find the model parameters that best
explain the data. The _maximum likelihood degree_ (ML degree) is the number of
complex solutions to the critical equations of this optimization problem,
which is also shown to be constant [HKS05, CHKS06].

Computing the the ML degree of phylogenetic varieties is a classical problem,
and some degrees are already known for certain trees and models, see Section
3. On the other hand, the study ED degrees of these varieties is an emerging
and understudied problem. To the best of our knowledge, only the ED degree of
one tree of $4$ leaves with a specific evolutionary model has been studied so
far (see [CFSGL21]).

In this work we aim to extend the computations of ED and ML degrees for some
group-based models on small phylogenetic trees with at most 5 leaves. For very
small trees and the simplest models these degrees can be obtained using
symbolic computations. However, when we increase the number of leaves of the
trees or the complexity of the model, these tools are not sufficient. In this
paper, we use tools from numerical algebraic geometry (using the package
HomotopyContinuation.jl [BT18]) to compute lower bounds for the ED and ML
degrees. Our main contribution in this paper is Theorem 4.1, where we find
lower bound of these algebraic degrees for trees with up to 5 leaves and group
based models. In order to obtain these numbers, Lemma 2.1 is also necessary,
in which we present the cardinalities of the generic fibers of the
parametrization map $\phi^{\mathcal{M}}_{T}$ for the group based models.

The paper is organized as follows, in Section 2 we introduce the basic notions
on phylogenetic trees, phylogenetic varieties and group based models. We also
present the Fourier transform which allows to parametrize phylogenetic
varieties with a monomial map. In Section 3 we describe in more detail the ML
degree and the ED degree for phylogentic varieties. Section 4 is dedicated to
present the results we have obtained using numerical methods, and these
results are finally discussed in Section 5.

##  2 Phylogenetic models

In this section, we introduce phylogenetic models and the algebraic varieties
associated to them. We focus on group-based phylogenetic models and introduce
a change of coordinates called the discrete Fourier transformation, in which
the parameterization of the phylogenetic varieties becomes monomial. The
reader is referred to [Fel03, SS03, PS05, SS05] for a more detailed
introduction to the topic.

Evolutionary relationships between a set of biological entities, or taxa, are
commonly depicted in the form of phylogenetic trees. A phylogenetic tree $T$
is a connected acyclic graph whose leaves are in one-to-one correspondence
with a set of labels. That is, the leaves $L(T)$ of a phylogenetic tree $T$
represent a set of current taxa, the inner nodes $\text{In}(T)$ correspond to
common ancestors and the edges $E(T)$ represent the evolutionary relationships
among them. For simplicity, we assume leaves are in correspondence with the
set $\\{1,2,\ldots,n\\}$ where $n$ is the number of leaves. We choose an
internal vertex of $T$ as the root $r$ of $T$, which induces an orientation on
the set of edges, from the root to the leaves.

The substitution of nucleotides in an evolutionary process is modeled by
adopting a parametric statistical model as follows: associate a random
variable $X_{v}$ taking values on a set of states
$\Sigma:=\\{s_{1},\ldots,s_{k}\\}$111$\Sigma$ is usually assumed to be the set
of nucleotides $\\{\tt A,\tt C,\tt G,\tt T\\}$ or the binary set $\\{0,1\\}$
that can be identified with the two groups of nucleotide bases purines and
pyrimidines. at each node $v$ of the tree. Also, consider as parameters a
distribution $\pi=(\pi_{1},\ldots,\pi_{k})$ at the root $r$ of the tree, where
$p_{i}=Prob(X_{r}=s_{i})$ and such that $\sum_{i}\pi_{i}=1$; and a $k\times k$
transition matrix $M_{e}$ at each (oriented) edge $e:u\to v$ of $T$ given by
$(M_{e})_{ij}=Prob(X_{v}=s_{j}|X_{u}=s_{i})$. The transition matrices are
stochastic (or Markov) matrices, that is, all its entries are non-negative and
its rows sum up to $1$.

Assuming a hidden Markov process on a phylogenetic tree, the probability of
observing a particular set of nucleotides at the leaves of the tree can be
computed as a polynomial in terms of the model parameters as follows,

$p_{x_{1},\ldots,x_{n}}=\sum_{\begin{subarray}{c}x_{r},x_{v}\in\Sigma\\\
{\scriptsize v\in\text{In}(T)}\end{subarray}}\pi_{x_{r}}\prod_{e\in
E(T)}M_{e}\left(x_{u},x_{v}\right),$ (1)

where $x_{r}$ is the state of the root, $x_{u}$ is the state of the parent
node $u$ of $e$ , and $x_{v}$ is the state of the child node $v$ of $e$. If
$e$ is a terminal edge ending at leaf $i$ then $x_{v}=x_{i}$.

We consider the polynomial map
$\phi_{T}^{\mathcal{M}}:\mathbb{R}^{d}\to\mathbb{R}^{m}$ associated to a tree
$T$ and model $\mathcal{M}$ that maps any $d$-tuple to a distribution vector
of the $m:=k^{n}$ possible observations at the leaves,

$\begin{array}[]{rccc}\phi_{T}^{\mathcal{M}}:&\Theta\subset\mathbb{R}^{d}&\to&\bigtriangleup^{m-1}\subset\mathbb{R}^{m}\\\
&\theta=\left(\theta_{1},\ldots,\theta_{d}\right)&\mapsto&p=(p_{\tt
s_{1}s_{1}\dots s_{1}},p_{\tt s_{1}s_{1}\dots s_{2}},p_{\tt s_{1}s_{1}\dots
s_{3}},\dots,p_{\tt s_{k}s_{k}\dots s_{k}}),\end{array}$ (2)

where $d$ is the number of free substitution parameters (that is, free entries
of the distribution at the root $\pi$ and the transition matrices $M_{e}$)
according to the model assumed and $\Theta$ denotes the space of stochastic
parameters. The coordinates $p_{x_{1},\ldots,x_{n}}$ correspond to the joint
distribution of characters at the leaves of $T$ given by a hidden Markov
process on $T$ governed by these parameters. Our assumption that the entries
of $\pi$ and rows of $M_{e}$ sum to $1$ implies $\sum
p_{x_{1},\ldots,x_{n}}=1$. Ignoring the stochastic restrictions on the
parameter space, we may regard
$\phi_{T}^{\mathcal{M}}:\mathbb{C}^{d}\to\mathbb{C}^{m}$ as a complex
polynomial map. The Zariski closure of the image of such complex map, denoted
$\mathcal{V}_{T}^{\mathcal{M}}=\overline{\phi_{T}^{\mathcal{M}}(\mathbb{C}^{d})}\in\mathbb{C}^{m}$,
is an algebraic variety known as the phylogenetic variety associated to $T$
and $\mathcal{M}$.

### Group-based models

Nucleotide substitution models $\mathcal{M}$ are specified by certain
restrictions, which may arise from biochemical properties, on the distribution
at the root and the transition matrices. The group-based model are certain
substitution models in which the set $\Sigma$ is given the additional
structure of an abelian group. That is, a group-based model relative to an
abelian group $G$ is a substitution model whose matrices can be defined in
terms of a function $f^{e}:G\to\mathbb{R}$ such that $M^{e}(g,h)=f^{e}(g-h)$,
for any $g,h\in G$ after the identification of $\Sigma$ with $G$.

In our work, we consider the Cavender-Farris-Neyman (CFN) model where
$\Sigma=\\{0,1\\}$ is identified with the group $\mathbb{Z}/2\mathbb{Z}$ and
the $3$-parameter Kimura (K81) model [Kim81] that identifies $\Sigma=\\{\tt
A,\tt C,\tt G,\tt T\\}$ with the group
$\mathbb{Z}/2\mathbb{Z}\times\mathbb{Z}/2\mathbb{Z}$. Both models assume a
uniform distribution at the root of the tree and respective transition
matrices have the form

$M_{e}=\left(\begin{array}[]{cc}a&b\\\ b&a\\\ \end{array}\right),\qquad
M_{e}=\left(\begin{array}[]{cccc}a&b&c&d\\\ b&a&d&c\\\ c&d&a&b\\\ d&c&b&a\\\
\end{array}\right).$

These two models are called general group-based models associated to a group
$G$. We also consider two submodels of the K81 model, the Jukes-Cantor
($JC69$) model [JC69] and the $2$-parameter Kimura (K80) model [Kim80].
Transition matrices for these models can be obtained from K81 matrices by
identifying parameters. The K80 model and the JC69 model assume $c=d$ and
$b=c=d$, respectively.

The parametrization map $\phi_{T}^{\mathcal{M}}$ in general is not one-to-one.
Then, understanding its generic fibers is important from both biological and
computational points of view. In what follows we present the cardinality of
such fibers, this corresponds to the number of parameters of the model that
are mapped by $\phi_{T}^{\mathcal{M}}$ to the same probability distribution at
the leaves of the tree.

######  Lemma 2.1.

The cardinalities $\delta_{T}^{\mathcal{M}}$ of the generic fibers of
$\phi_{T}^{\mathcal{M}}$ for the previously described group based models are

$\delta_{T}^{CFN}=2^{N},\quad\delta_{T}^{JC}=1,\quad\delta_{T}^{K80}=2^{N}\quad\mbox{and}\quad\delta_{T}^{K81}=4^{N}$
(3)

where $N$ is the number of interior nodes of $T$.

###### Proof.

The result is already known for any general group-based models, see
Proposition 2.1 in [CFSM15] where the authors prove that the cardinality of
the generic fiber of $\phi_{T}^{\mathcal{M}}$ is $\lvert G\rvert^{N},$ where
$G$ is the group corresponding to the model $\mathcal{M}$.

For the other models, the result follows from the the fact that certain
permutations of the rows and columns of the Markov matrices
$\left\\{M_{e}\right\\}_{e\in E(T)}$ can be applied without affecting the join
distribution $p$ (see [AR03] or [Cha96]). Moreover, according to [AR04], given
a set of parameters $\left\\{\pi,\\{M_{e}\\}_{e\in E(T)}\right\\}$ and a set
of permutations matrices $P_{v_{1}},\ldots,P_{v_{N}}$ for each interior node
$v_{i}\in Int(T)$, the distribution
$p=\phi_{T}\left(\pi,\left\\{M_{e}\right\\}_{e\in E(T)}\right)$ equals the
distribution $\phi_{T}\left(\pi
P_{r}^{T},\left\\{P_{u}M_{e}P_{v}^{T}\right\\}_{e:u\to v}\right)$, with
$P_{v}=Id$ if $v$ is a leaf of the tree. This implies that there are at most
$(|\Sigma|!)^{N}$ set of parameters that map to the same distribution.

However, if we fix a particular model, the number of permutations that one can
consider such that the structure of the transition matrices is preserved is
much smaller than $|\Sigma|!$. We start considering the K81 model and a tree
$T$. Let $e:u\to v$ be a terminal edge of $T$, that is, $v$ is a leaf, and
$M_{e}$ the transition matrix associated to $e$. Consider a permutation matrix
$P_{u}$ associated to the vertex $u$, then we can replace $M_{e}$ by the
matrix $P_{u}M_{e}$. It can be easily checked that, if $M_{e}$ is a K81
matrix, then $P_{u}M_{e}Id$ is a K81 matrix if and only if the permutations
$P_{u}$ is one of the following matrices:

$\Gamma_{1}=Id,\ \Gamma_{2}=\left(\begin{array}[]{cccc}0&1&0&0\\\ 1&0&0&0\\\
0&0&0&1\\\ 0&0&1&0\\\ \end{array}\right),\
\Gamma_{3}=\left(\begin{array}[]{cccc}0&0&1&0\\\ 0&0&0&1\\\ 1&0&0&0\\\
0&1&0&0\\\ \end{array}\right)\
\Gamma_{4}=\left(\begin{array}[]{cccc}0&0&0&1\\\ 0&0&1&0\\\ 0&1&0&0\\\
1&0&0&0\\\ \end{array}\right).$ (4)

Consider now and edge $e:u\to v$ such that $v$ is the parent node of a leaf of
$T$. Let $P_{u}$ and $P_{v}$ be two permutation matrices associated to $u$ and
$v$. As shown above, to keep the structure of the transition matrices, $P_{v}$
is one of the matrices $\\{\Gamma_{1},\Gamma_{2},\Gamma_{3},\Gamma_{4}\\}$.
Then, it is straightforward to see that $P_{u}M_{e}P_{v}^{T}$ is a K81 matrix
if and only if $P_{u}$ is also one of the permutation matrices in (4). This
same argument can be applied recursively to all edges of the tree, showing
that the only permutation matrices $P_{u}$ and $P_{v}$ that can be considered
for $PuM_{e}P_{v}^{T}$ to preserve the K81 structure are the matrices in
$\\{\Gamma_{1},\Gamma_{2},\Gamma_{3},\Gamma_{4}\\}$. This implies that the
number of parameters mapping to the same distribution by
$\phi_{T}^{\mathcal{M}}$ is $4^{N}$. Recursively along all edges of $T$ one
can see that all chosen $P_{u}$ for any inner node $u$ should be in
$\\{\Gamma_{1},\Gamma_{2},\Gamma_{3},\Gamma_{4}\\}$ so that $P_{u}M_{e}Pv^{T}$
are K81 matrices.

A similar argument can be applied to the K80 model, but in this case there are
only two matrices, $\Gamma_{1}$ and $\Gamma_{2}$ that preserve the structure
of the transition matrices, which implies that the cardinality of the fibers
of $\phi_{T}^{K80}$ is $2^{N}$. Finally, $\delta_{T}^{JC69}=1$ follows from
the fact that the only permutation $P_{u}$ that preserves the structure of the
JC69 matrices is the identity matrix. ∎

For more insight on group based models and their corresponding phylogenetic
varieties, we refer to [SS05], [Mic11] and [CFSM15].

### Fourier transform

The discrete Fourier transform is a powerful tool to analyse properties of
group-based models and the pattern distributions arising from them. This idea
was introduced by [Hen89] and [HP89] and was further explored by [SSE93]. The
discrete Fourier transform provides a linear change of coordinates that
reveals the irreducible variety of distributions of a group-based model
$\mathcal{V}_{T}^{\mathcal{M}}$ to be a toric variety. The following results
provide an explicit description of the parameterization of these varieties
after this linear change of coordinates.

######  Theorem 2.2 ([ES93]).

Let $p$ be the joint distribution of a group based model for a phylogenetic
tree $T$ obtained as in (1). Then the Fourier transform of $p$ has the form

$q_{y_{1}\ldots y_{n}}=\begin{cases}\prod_{e\in E(T)}x_{k}^{e}&\text{if
$y_{n}=y_{1}+\cdots+y_{n-1}$,}\\\ 0&\text{otherwise}\end{cases}$ (5)

where $k=\sum_{i\in de(e)}y_{i}$, and $de(e)$ is the set of all leaves $l$
descendant of $e$222A leaf $l$ is a descendant of an edge $e$ if there is a
directed path from $e$ to $l$. and $x^{e}$ are the Fourier parameters
associated to edge $e$.

Fourier parameters are the Fourier transforms of the transition functions
$f^{e}:G\to\mathbb{R}$, which can be shown to be equivalent to the eigenvalues
of transition matrices $M_{e}$. For simplicity, we assume
$x_{0}^{e},\ldots,x_{l}^{e}$ are the $l+1$ different eigenvalues of $M_{e}$,
where $l$ coincides with the number of free parameters in $M_{e}$. As a
consequence of the previous theorem, we can consider a map
$\varphi_{T}^{\mathcal{M}}:\mathbb{C}^{d}\to\mathbb{C}^{m}$ parameterizing
phylogenetic varieties from Fourier parameters to Fourier coordinates. Since
the map $\varphi_{T}^{\mathcal{M}}$ is homogeneous we can consider the Zariski
closure of the image of such map as a projective variety in
$\mathbb{P}^{m-1}:$

$\begin{array}[]{rccc}\varphi_{T}^{\mathcal{M}}:&\mathbb{P}^{l}\times\cdots\times\mathbb{P}^{l}&\to&\mathbb{P}^{4^{n}-1}\\\
&\left([x_{0}^{e_{1}}:\dotsm:x_{l}^{e_{1}}],\ldots,[x_{0}^{e_{E}}:\dotsm:x_{l}^{e_{E}}]\right)&\mapsto&q=(q_{\tt
A\dots A},\dots,q_{\tt T\dots T}),\end{array}$ (6)

where $E$ denotes the number of edges $E(T)$ of $T$.

######  Remark 2.3.

The variety $\overline{\mathcal{V}_{T}^{\mathcal{M}}}=\overline{\text{Im
}\varphi_{T}^{\mathcal{M}}}$ is the projectivization of the affine
phylogenetic variety $\mathcal{V}_{T}^{\mathcal{M}}$. That is, the affine
variety $\mathcal{V}_{T}^{\mathcal{M}}$ corresponds to the affine chart given
by $x_{0}^{e_{1}}=\dots=x_{0}^{e_{E}}=1$.

Since the discrete Fourier transform is a linear change of coordinates, it
follows from Lemma 2.1 that the cardinallity of the fibers of the
parametrization map on Fourier parameters $\varphi_{T}^{\mathcal{M}}$ is
preserved.

######  Corollary 2.4.

The cardinalities of the generic fibers of the map $\varphi_{T}^{\mathcal{M}}$
for the models CFN, JC69, K80, K81 models are $2^{N}$, $1$, $2^{N}$ and
$4^{N}$ respectively.

##  3 Algebraic degrees for phylogenetic models

In this section, we discuss an important problem in mathematical
phylogenetics: given some data $u$ find a distribution $p$ in
$\mathcal{V}_{T}$ that makes $u$ most likely. Such distribution $p$ may depend
on the distance measure chosen. In this work we focus on two classical
optimization problems in algebraic statistics, the maximum likelihood estimate
and the closest point on a variety in the Euclidean distance. First, we
introduce the Maximum likelihood estimator problem and then the problem of
minimizing the Euclidean distance over an algebraic variety.

### Maximum likelihood estimator

Maximum likelihood estimation is a standard approach to parameter estimation,
and a fundamental computational task in statistics. In the case of
phylogenetics, given observed data $u$, a phylogenetic tree $T$, and a
nucleotide substitution model $\mathcal{M}$, the maximum likelihood estimate
is the set of parameters $\theta=(\theta_{1},\ldots,\theta_{d})$ that is most
likely to have produced the data $u$. We assume data is given by a non-
negative integer vector $u$ in $\mathbb{N}^{n}$, where $u_{i}$ is the number
of times the $i$-th event is observed. Then, the problem of maximum likelihood
estimation aims to find the parameters $\theta$ that maximizes the log-
likelihood function

$\ell_{u}(\theta)=\sum_{i=1}^{m}u_{i}\ log(p_{i}),$ (7)

where $p=(p_{1},\ldots,p_{m})=\phi_{T}^{\mathcal{M}}(\theta).$ Local and
global optimizers of the log-likelihood function are solutions of the critical
equations of $\ell_{u}(\theta)$, which are given by

$\frac{\partial\ell_{u}}{\partial\theta_{1}}=\dotsm=\frac{\partial\ell_{u}}{\partial\theta_{d}}=0$
(8)

For generic data $u$, the number of complex solutions
$\theta\in\mathbb{C}^{d}$ of the critical equations
$\frac{\partial\ell_{u}}{\partial\theta_{i}}=0$ such that
$\phi_{T}^{\mathcal{M}}(\theta)$ is a smooth point in $\mathcal{V}_{T}$ is
finite and constant. (see [HKS05, CHKS06]):

######  Definition 3.1.

The maximum likelihood degree (ML degree) of a tree $T$ and model
$\mathcal{M}$ is the number of solutions of the system of critical equations
of $\ell_{u}(\theta)$, for generic data $u$, divided by
$\delta_{T}^{\mathcal{M}}$, the cardinality of the fibers of
$\phi_{T}^{\mathcal{M}}$(see Remark 2.1).

ML degrees of phylogenetic varieties have been studied extensively, and they
are known for some specific trees and models. For instance, ML degree for the
$3$-leaf star tree evolving under Jukes Cantor model is $23$ and can be found
in [HKS05]. Known ML degrees of phylogenetic varieties are resumed in [CGS05]
and the companion website

https://www.coloradocollege.edu/aapps/ldg/small-trees/.

In section 4, we extend these computations, and present new ML degrees that
have been obtained using numerical methods.

### Minimizing the Euclidean distance

A different approach to find parameters that better explain a data point $u$
is to minimize the Euclidean distance function from $u$ to a given algebraic
variety $\mathcal{V}_{T}^{\mathcal{M}}$. In this case, we assume that the data
$u\in\mathbb{R}^{m}$ is a distribution. The aim is to minimize the squared
Euclidean distance function from $u$ to the variety
$\mathcal{V}_{T}^{\mathcal{M}}$. In order to accomplish this, we consider the
Fourier parameterization of the variety introduced in (6). Since the Fourier
transformation is an orthogonal linear change of coordinates, the Euclidean
distance between distributions can be computed using these new coordinates
(see Remark 2.5 in [CFSGL21] for more details). Then we consider the following
optimization problem

$\min_{x}d_{u}(x):=\min_{q\in\overline{\mathcal{V}_{T}^{\mathcal{M}}}}\sum_{i=1}^{m}(u_{i}-q_{i})^{2},$
(9)

where
$x=\left([x_{0}^{e_{1}}:\dotsm:x_{l}^{e_{1}}],\ldots,[x_{0}^{e_{E}}:\dotsm:x_{l}^{e_{E}}]\right)$
are the Fourier parameters and
$q=(q_{1},\ldots,q_{m})=\varphi_{T}^{\mathcal{M}}(x)$ the Fourier coordinates.
The critical points of this problem are the solutions of the system

$\frac{\partial d_{u}}{\partial x_{1}}=\dots=\frac{\partial d_{u}}{\partial
x_{d}}=0.$ (10)

For generic data $u$, the number of complex solutions $x\in\mathbb{C}^{d}$ of
the critical equations $\frac{\partial d_{u}}{\partial x_{i}}=0$ such that
$\phi_{T}^{\mathcal{M}}(x)$ is a smooth point in $\mathcal{V}_{T}$ is
constant, see [DHO+16].

######  Definition 3.2.

the Euclidean distance (ED) degree of $\mathcal{V}_{T}^{\mathcal{M}}$ is the
number of complex solutions of the critical equations $\frac{\partial
d_{u}}{\partial x_{i}}=0$, for generic data $u$, divided by
$\delta_{T}^{\mathcal{M}}$.

In this work we also consider a more general problem, that is, given a vector
$\lambda=(\lambda_{1},\ldots,\lambda_{m})$, consider the minimization problem

$\min_{x}d_{u,\lambda}(x)=\min_{q\in\overline{\mathcal{V}_{T}^{\mathcal{M}}}}\sum_{i=1}^{m}\lambda_{i}(u_{i}-q_{i})^{2}.$
(11)

If the weight vector $\lambda$ is generic, then the number of complex critical
points of (11), that is, the solutions of

$\frac{\partial d_{u,\lambda}}{\partial x_{1}}=\dots=\frac{\partial
d_{u,\lambda}}{\partial x_{d}}=0,$ (12)

is independent of $\lambda$ and is called the generic ED degree (gED degree).
Since $\overline{\mathcal{V}_{T}^{\mathcal{M}}}$ is a complex toric variety,
then the parameterization $\varphi_{T}^{\mathcal{M}}$ is a monomial map and
can be described by a $d\times m$ matrix $A$ whose columns are the exponent
vectors of the monomials parameterizing each Fourier coordinate. Thus, the
generic ED degree can be computed in terms of the Chern-Mather volumes of all
faces of the polytope $P=\text{conv}(A)$. For the explicit expression of the
gED degree for toric varieties, the reader is referred to [HS18, Theorem 1.1].

######  Remark 3.3.

According to Theorem 6.11 in [DHO+16], the gED degree of
$\overline{\mathcal{V}_{T}^{\mathcal{M}}}$ is an upper bound of both the ED
degree of $\overline{\mathcal{V}_{T}^{\mathcal{M}}}$ and the ED degree of
${\mathcal{V}_{T}^{\mathcal{M}}}$.

The study of the optimization problem (9), i.e. minimizing the Euclidean
distance from a data point to a phylogenetic variety, as well as the
calculation of ED degrees for these varieties is still recent. As far as we
are aware, only the ED degree for the 4-leaf binary tree evolving under the
JC69 model is previously known, see [CFSGL21, Lemma 5.1].

##  4 Computational results

In this section we present the computational results that we obtained for the
ML and ED degrees for phylogenetic trees up to $5$ leaves and the group based
models previously introduced. We first discuss the computational approaches
used, then we present the computed degrees and finally the execution times for
such computations.

### Computational approaches

In order to compute the described algebraic degrees one could use symbolic
computations following the algorithms described in [HKS05, DHO+16]. However
these approaches have some limitations and are, in general, not feasible for
models such as the Kimura models, for example. However, we use the Macaulay2
package toricInvariants (see [HS18]) to compute the gED degree of the
projective toric varieties parameterized by (5), that has given us the
expected results in most cases.

Our approach in this paper is to use numerical methods to find out how much we
can extend the computations of these invariants for phylogenetic varieties. We
use the Julia package HomotopyContinuation.jl [BT18] for computing both, the
ED degrees and gED degrees for the toric phylogenetic varieties. To compute
these degrees we turn the computation of critical points of the optimization
problems (9) and (11) into the problem of solving the systems of polynomial
equations (10) and (12). We solve these systems of critical equations using
the solve function in HomotopyContinuation.jl. In order to solve the system, a
starting system, with easily computed solutions, is first derived using a
polyhedral homotopy algorithm developed by Hubert and Sturmfels [HS95]. After
that, the original system is subjected to a parameter homotopy, which along
with the solutions of the starting system define the solution path that we can
track.

To compute the ML degree we solve the system of critical equations given in
(8). This system can be converted into a polynomial system by adding new
variables as suggested in [PS05, Section 3.3]. In this way one can also use
HomotopyContinuation.jl. However, for our specific models, solving the
corresponding systems became computationally infeasible as the number of
leaves and parameters in the model increased. This is mostly due to the fact
that the number of paths to track increases at a significantly higher rate. As
a result, most computations did not finish due to time and memory constrains.

As an alternative approach, we use _monodromy_ to compute all the isolated
solutions with the specified parameter values and initial solutions. We search
for a pair of data $u_{0}$ and parameters $\theta_{0}$ satisfying the system
(8). New solutions for $\theta$ are found by making loops in the parameter
space of the system. In order to use monodromy one must compute a starting
solution in terms of variables $\theta_{0}$ and parameters $u_{0}$. One
approach to do this is to substitute some random variables $\theta_{0}$ into
the system (8) and then find initial parameters $u_{0}$ by solving the
resulting linear system. However, in our systems the number of parameters is
larger than the number of variables and therefore, we get an underdetermined
linear system. To solve it we assign random numbers to entries of $u_{0}$
until we get a full rank square linear system which we then solve for the
remaining entries of $u_{0}$.

Finally, we use a method based on interval arithmetic to certify solutions
obtained via numerical methods to the polynomial systems above, see [BRT21].
The method certify in the package HomotopyContinuation.jl uses interval
arithmetic to prove the correctness of an isolated nonsingular solution to a
(square) system of polynomial equations. Given a numerical approximation of a
complex zero, this package computes intervals for the real and imaginary parts
of each variable that contain a true solution to the system and certifies
whether the proposed numerical solution lies in this interval. This process
may not be able to certify all true solutions. Thus the number of certified
numerical solutions to our system of equations is a lower bound for the actual
algebraic degrees of the corresponding varieties, see [BRT21, Section 1.1].
This allows us to turn our numerical computations into rigorous proof of our
results in the following subsection.

The implementation of all methods discussed in this contribution is found at

https://github.com/marinagarrote/Computing-algebraic-degrees-of-phylogenetic-
varieties.

### Algebraic degrees of phylogenetic varieties for small trees

We display the results of our computations of the algebraic degrees of the
corresponding phylogenetic varieties in the tables Table 1, Table 2, and Table
3, . We also present the degree and dimension of the projective variety
$\overline{\mathcal{V}_{T}^{\mathcal{M}}}$. The latter invariant corresponds
to the number of free parameters needed to parameterize the variety plus one.
Observe that, by intersecting the corresponding variety with the hyperplane
given by $\left\\{q_{\tt A\dots\tt A}=1\right\\}$, the dimension drops by one
giving the dimension of ${\mathcal{V}_{T}^{\mathcal{M}}}$ which is the minimum
number of parameters needed to determine the model, see [CGS05].

### 3-leaf trees:

We start by considering the binary $3$-leaf tree.

Tree |  Model |  dim |  deg |  ML deg |  ED deg of ${\mathcal{V}_{T}^{\mathcal{M}}}$ |  gED deg of ${\mathcal{V}_{T}^{\mathcal{M}}}$ |  ED deg of $\overline{\mathcal{V}_{T}^{\mathcal{M}}}$ |  gED deg of $\overline{\mathcal{V}_{T}^{\mathcal{M}}}$
---|---|---|---|---|---|---|---|---
|  CFN |  $4$ |  $2$ |  $1$ |  $1$ |  $1$ |  $1$ |  $1$
|  JC69 |  $4$ |  $3$ |  $23$ |  $13$ |  $13$ |  $18$ |  $18$
|  K80 |  $7$ |  $12$ |  1964 |  $184$ |  $184$ |  $188$ |  $236$
|  K81 |  $10$ |  $96$ |  |  $12673$ |  $15553$ |  $10360$ |  $21112$
Table 1:  Dimension, degree, ML degree of ${\mathcal{V}_{T}^{\mathcal{M}}}$,
ED and gED degrees of ${\mathcal{V}_{T}^{\mathcal{M}}}$ and
$\overline{\mathcal{V}_{T}^{\mathcal{M}}}$ for the binary $3$-leaf tree
illustrated in the first column and the four group-based models presented in
Section 2.

### 4-leaf trees:

In the case of phylogenetic trees with four leaves, we consider both the star
tree and the binary tree.

Tree |  Model |  dim |  deg |  ML deg |  ED deg of ${\mathcal{V}_{T}^{\mathcal{M}}}$ |  gED deg of ${\mathcal{V}_{T}^{\mathcal{M}}}$ |  ED deg of $\overline{\mathcal{V}_{T}^{\mathcal{M}}}$ |  gED deg of $\overline{\mathcal{V}_{T}^{\mathcal{M}}}$
---|---|---|---|---|---|---|---|---
|  CFN |  $5$ |  $8$ |  $92$ |  $92$ |  $92$ |  $72$ |  $120$
|  JC69 |  $5$ |  $20$ |  $4315$ |  $124$ |  $220$ |  $243$ |  $267$
|  K80 |  $9$ |  $832$ |  |  $50972$ |  $88076$ |  $33176$ |  $106744$
|  CFN |  $6$ |  $4$ |  $14$ |  $10$ |  $30$ |  $4$ |  $36$
|  JC69 |  $6$ |  $34$ |  $17332$ |  $290$ |  $630$ |  $309$ |  $809$
Table 2:  Dimension, degree, ML degree of ${\mathcal{V}_{T}^{\mathcal{M}}}$,
ED and gED degrees of ${\mathcal{V}_{T}^{\mathcal{M}}}$ and
$\overline{\mathcal{V}_{T}^{\mathcal{M}}}$ for the star $4$-leaf tree evolving
under the CFN, JC69 and K80 models and also for the binary $4$-leaf tree with
the CFN and JC69 models.

### 5-leaf trees:

Finally, we consider three different $5$-leaf trees, the star tree, an
unresolved tree and the binary tree.

Tree |  Model |  dim |  deg |  ML deg |  ED deg of ${\mathcal{V}_{T}^{\mathcal{M}}}$ |  gED deg of ${\mathcal{V}_{T}^{\mathcal{M}}}$ |  ED deg of $\overline{\mathcal{V}_{T}^{\mathcal{M}}}$ |  gED deg of $\overline{\mathcal{V}_{T}^{\mathcal{M}}}$
---|---|---|---|---|---|---|---|---
|  CFN |  $6$ |  $52$ |  $4698$ |  $858$ |  $1098$ |  $484$ |  $1364$
|  JC69 |  $6$ |  $115$ |  |  $713$ |  $2313$ |  $1792$ |  $2792$
|  CFN |  $7$ |  $44$ |  $2098$ |  $334$ |  $1026$ |  $176$ |  $1256$
|  JC69 |  $7$ |  $315$ |  |  $3819$ |  $12787$ |  $4569$ |  $15977$
|  CFN |  $8$ |  $34$ |  $668$ |  $144$ |  $1044$ |  $62$ |  $1294$
|  JC69 |  $8$ |  $813$ |  |  $10073$ |  $62057$ |  $9010$ |  $78583$
Table 3:  Dimension, degree, ML degree of ${\mathcal{V}_{T}^{\mathcal{M}}}$,
ED and gED degrees of ${\mathcal{V}_{T}^{\mathcal{M}}}$ and
$\overline{\mathcal{V}_{T}^{\mathcal{M}}}$ for the star $5$-leaf tree, an
unresolved $5$-leaf tree and the binary $5$-leaf tree evolving under the CFN
and the JC69 model.

We have certified the numerical solutions of the systems (8), (10) and (12)
which allows us to state the following result:

######  Theorem 4.1.

The results presented in Tables 1, 2, and 3 are lower bounds for the
corresponding algebraic degrees of the respective phylogenetic varieties.

### Execution time.

The results presented in this section were computed on an iMac with the Apple
M1 chip and 16 GB of memory running macOS Monterey 12.6, Julia version 1.8.1
and Macaulay2 version 1.20. We used versions 2.6.4 and 3.01 of the packages
HomotopyContinuation.jl and toricInvariants respectively. In the following
Table 4 we display the times (in seconds) needed to compute the ML and ED
degrees of ${\mathcal{V}_{T}^{\mathcal{M}}}$ for all small trees and models
considered above. We report timings using both homotopy continuation and
monodromy methods and discuss the results in the next section. The blanks in
the table correspond to cases that we have not been able to compute, see
Section 5, while the entries marked with * in the last column correspond to
computational times with an extremely high standard deviation (a few seconds
to several hours) or computations that were halted after several hours.

Tree |  Model |  ML deg (HC) |  ML deg (Mon) |  ED deg of ${\mathcal{V}_{T}^{\mathcal{M}}}$ (HC) |  ED deg of ${\mathcal{V}_{T}^{\mathcal{M}}}$ (Mon)
---|---|---|---|---|---
|  CFN |  $0.03$ |  $0.02$ |  $0.04$ |  $0.03$
|  JC69 |  $3.2$ |  $1.8$ |  $0.04$ |  $0.05$
|  K80 |  $437$ |  $15$ |  $1$ |  $0.7$
|  K81 |  |  |  $92$ |  $190$
|  CFN |  $12$ |  $2.7$ |  $0.07$ |  $0.2$
|  JC69 |  $120$ |  $20$ |  $0.08$ |  $0.12$
|  K80 |  |  |  $140$ |  $3002$
|  CFN |  $89$ |  $3.7$ |  $0.5$ |  7.5
|  JC69 |  |  $358$ |  $0.5$ |  *
|  CFN |  |  $37$ |  $0.5$ |  $5$
|  JC69 |  |  |  $0.8$ |  $2$
|  CFN |  |  $78$ |  $10$ |  *
|  JC69 |  |  |  $13$ |  *
|  CFN |  |  $284$ |  $57$ |  *
|  JC69 |  |  |  $136$ |  *
Table 4:  Computational times to obtain the ML and ED degrees of
${\mathcal{V}_{T}^{\mathcal{M}}}$ using both homotopy continuation and
monodromy.

##  5 Discussion

In this work, we have highlighted the significant advantage of numerical
methods versus symbolic approaches to solve systems of polynomial equations
related to computing algebraic degrees of phylogenetic varieties. Most degrees
computed in this project were beyond the reach of symbolic methods, that
halted after long computational times due to machine memory problems. We have
verified that symbolic computations can only be carried out for trees under
the CFN model, some under the JC69 model and for the generic ED degree of
projective varieties. With that said, the approach described in [HS18] allowed
us to compute all but one of the gED degrees of the corresponding projective
varieties (the star tree with 3 leaves under the K81 model was the lone
exception). These symbolic computations were performed using the M2 package
ToricInvariants.m2 [Hel]. It is also worth noting that we have not been able
to verify numerically the generic ED degree for the projective variety of the
binary tree on 5 leaves under the JC69 model. The closest we have gotten is
78,577 certified solutions.

Although numerical methods allow us to present several new bounds for the
Euclidean distance degrees, we have not yet been able to compute the ML
degrees for many of the phylogenetic varieties we considered. Moreover, there
is evidence that we have reached the boundary of problems that are tractable
using the existing numerical techniques. These missing degrees, along with the
degrees of larger trees and more general models, remain an open question and
our next step in our studies.

It is also important to note the differences between homotopy continuation and
monodromy. In general, for all missing degrees in our tables, homotopy
continuation methods run out of memory during the computation of a starting
system and monodromy methods have excessively long computation times.
Furthermore, in our calculations we observe that although for most cases
monodromy is slower, it manages to complete the computation for some cases
that homotopy continuation is unable to do, as can be seen in Table 4.

In the first column of Table 4, all the missing ML degree computations failed
due to memory overflow while attempting to compute a start system for homotopy
continuation. In the second column, the missing entries are not due to memory
issues but prohibitively long running times. For example, the computation for
the star tree with 3 leaves under the K81 model halted after running for 6
days. Using monodromy, the system had found 10,080,070 solutions and exited
half way through the certification process. Similarly, for the binary tree
with 4-leaves under the K80 model, after 3 days of running time, using
monodromy, the system had computed 2,466,126 solutions. At that moment, we
paused the process and certified 62,141 solutions. This also demonstrates how
fast these degrees grow when we add more free parameters in the models.

In our work, we also kept track of the number of real solutions to our
systems. According to our simulations, generally the systems studied do not
have real and stochastic solutions. Although in a few occasions and for some
random initial data we have been able to find one or two stochastic solutions.
This shows that, from a biological point of view, if one wants to find the
point of the variety with stochastic (i.e. biologically meaningful) parameters
closest to a given data point, one should take into account the boundary of
the variety region defined by distributions that are the image of stochastic
parameters. This, therefore, adds considerable complexity to the problem of
finding the tree and parameters that best explain a given data point.

Finally, it would be interesting to extend our computations to include small
phylogenetic networks. Networks, are more general models that allow to
considering more complex evolutionary processes such as hybridization and gene
flow, see [GL18] for a nice introduction to the topic. Moreover, we believe
that it is possible to design tailored homotopy continuation methods that may
allow the computation of significant lower bounds for some of the missing
algebraic degrees and we propose to explore this direction further.

## Acknowledgments

The authors would like to thank Bernd Sturmfels for introducing the three of
us and encouraging us to work on this project. They also wish to thank Paul
Breiding, Sonja Petrović and Kemal Rose for helpful discussions on this topic.
The research of Elima Shehu was funded by the Deutsche Forschungsgemeinschaft
(DFG, German Research Foundation), Projektnummer 445466444.

## References

  * [AR03] Elizabeth S. Allman and John A. Rhodes.  Phylogenetic invariants for the general Markov model of sequence mutation.  Mathematical Biosciences, 186(2):113–144, dec 2003.
  * [AR04] Elizabeth S. Allman and John A. Rhodes.  Quartets and Parameter Recovery for the General Markov Model of Sequence Mutation.  Applied Mathematics Research eXpress, 2004(4):107–131, 2004.
  * [BRT21] Paul Breiding, Kemal Rose, and Sascha Timme.  Certifying zeros of polynomial systems using interval arithmetic, 2021\.
  * [BT18] Paul Breiding and Sascha Timme.  HomotopyContinuation.jl: A Package for Homotopy Continuation in Julia.  In Mathematical Software – ICMS 2018, pages 458–465, Cham, 2018\. Springer International Publishing.
  * [CFSGL21] Marta Casanellas, Jesús Fernández-Sánchez, and Marina Garrote-López.  Distance to the stochastic part of phylogenetic varieties.  Journal of Symbolic Computation, 104:653–682, may 2021.
  * [CFSM15] Marta Casanellas, Jesús Fernández-Sánchez, and Mateusz Michałek.  Low degree equations for phylogenetic group-based models.  Collectanea Mathematica, 66(2):203–225, may 2015.
  * [CGS05] Marta Casanellas, Luis David Garcia, and Seth Sullivant.  Catalog of Small Trees.  In L Pachter and B Sturmfels, editors, Algebraic Statistics for Computational Biology, chapter 15, pages 305–321. Cambridge University Press, aug 2005.
  * [Cha96] Joseph T. Chang.  Full reconstruction of Markov models on evolutionary trees: Identifiability and consistency.  Mathematical Biosciences, 137(1):51–73, oct 1996.
  * [CHKS06] Fabrizio Catanese, Serkan Hoşten, Amit Khetan, and Bernd Sturmfels.  The maximum likelihood degree.  American Journal of Mathematics, 128(3):671–697, 2006.
  * [DHO+16] Jan Draisma, Emil Horobeţ, Giorgio Ottaviani, Bernd Sturmfels, and Rekha R. Thomas.  The Euclidean Distance Degree of an Algebraic Variety.  Foundations of Computational Mathematics, 16(1):99–149, feb 2016\.
  * [ES93] Steven N. Evans and T. P. Speed.  Invariants of Some Probability Models Used in Phylogenetic Inference.  The Annals of Statistics, 21(1):355–377, mar 1993.
  * [Fel03] Joseph Felsenstein.  Inferring phylogenies.  Sinauer Associates, 2003.
  * [GL18] Elizabeth Gross and Colby Long.  Distinguishing phylogenetic networks.  SIAM Journal on Applied Algebra and Geometry, 2(1):72–93, 2018\.
  * [Hel] Martin Helmer.  Toricinvariants: A _Macaulay2_ package. version 3.01.
  * [Hen89] Michael D. Hendy.  The Relationship Between Simple Evolutionary Tree Models and Observable Sequence Data.  Systematic Zoology, 38(4):310, dec 1989.
  * [HKS05] Serkan Hosten, Amit Khetan, and Bernd Sturmfels.  Solving the likelihood equations.  Foundations of Computational Mathematics, 5(4):389–407, 2005.
  * [HP89] Michael D. Hendy and David Penny.  A Framework for the Quantitative Study of Evolutionary Trees.  Systematic Zoology, 38(4):297, dec 1989.
  * [HS95] Birkett Huber and Bernd Sturmfels.  A polyhedral method for solving sparse polynomial systems.  Math. Comp., 64(212):1541–1555, 1995.
  * [HS18] Martin Helmer and Bernd Sturmfels.  Nearest points on toric varieties.  Mathematica Scandinavica, 122(2):213–238, Apr. 2018.
  * [JC69] Thomas H. Jukes and Charles R. Cantor.  Evolution of protein molecules.  Mammalian protein metabolism, 3:21–132, 1969.
  * [Kim80] Motoo Kimura.  A simple method for estimating evolutionary rates of base substitutions through comparative studies of nucleotide sequences.  Journal of Molecular Evolution, 16(2):111–120, jun 1980.
  * [Kim81] Motoo Kimura.  Estimation of evolutionary distances between homologous nucleotide sequences.  Proceedings of the National Academy of Sciences, 78(1):454–458, jan 1981.
  * [Mic11] Mateusz Michałek.  Geometry of phylogenetic group-based models.  Journal of Algebra, 339(1):339–356, aug 2011.
  * [PS05] Lior Pachter and Bernd Sturmfels, editors.  Algebraic statistics for computational biology.  Cambridge University Press, New York, 2005.
  * [SS03] Charles Semple and Mike Steel.  Phylogenetics.  Oxford University Press, 2003.
  * [SS05] Bernd Sturmfels and Seth Sullivant.  Toric Ideals of Phylogenetic Invariants.  Journal of Computational Biology, 12(2):204–228, mar 2005.
  * [SSE93] László A.A. Székely, Mike A. Steel, and Péter L. Erdös.  Fourier Calculus on Evolutionary Trees.  Advances in Applied Mathematics, 14(2):200–216, jun 1993.


# Attention Based Fully Convolutional Network for Speech Emotion Recognition

Yuanyuan Zhang, Jun Du, Zirui Wang, Jianshu Zhang, Yanhui Tu National
Engineering Laboratory for Speech and Language Information Processing,
University of Science and Technology of China, Hefei, China
{zyuan, cs211, xysszjs<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

Speech emotion recognition is a challenging task for three main reasons: 1)
human emotion is abstract, which means it is hard to distinguish; 2) in
general, human emotion can only be detected in some specific moments during a
long utterance; 3) speech data with emotional labeling is usually limited. In
this paper, we present a novel attention based fully convolutional network for
speech emotion recognition. We employ fully convolutional network as it is
able to handle variable-length speech, free of the demand of segmentation to
keep critical information not lost. The proposed attention mechanism can make
our model be aware of which time-frequency region of speech spectrogram is
more emotion-relevant. Considering limited data, the transfer learning is also
adapted to improve the accuracy. Especially, it’s interesting to observe
obvious improvement obtained with natural scene image based pre-trained model.
Validated on the publicly available IEMOCAP corpus, the proposed model
outperformed the state-of-the-art methods with a weighted accuracy of 70.4%
and an unweighted accuracy of 63.9% respectively.

##  1 Introduction

Emotions play an important role in human communications [1] and successfully
detecting the emotion states is helpful to improve the efficiency of human-
computer interaction. For instance, in call centers, tracking customers’
emotion states can be useful for quality measurement [2] and the calls from
angry customers can therefore be assigned to experienced agents. Speech is one
of the communication channels that emotions could have serious influence on.
Technically, emotions affect both the voice characteristics and linguistic
content. In this study, we focus on the change of voice characteristics to
recognize the underlying emotions in speech.

Figure 1: The traditional speech emotion recognition system.

Speech emotion recognition (SER) has been an active research field for decades
[3, 4, 5, 6]. We demonstrate the architecture of traditional approaches for
SER in Figure 1. First, acoustic features which are believed to incorporate
the information of human emotions are extracted from raw speech waveform frame
by frame. The features include pitch, voicing probability, energy, etc. Then
various statistical functions (e.g. mean, max, linear regression coefficients,
etc.) are applied to the frame-level features. And the outputs are
concatenated as a feature vector to represent the whole utterance. Finally,
the utterance feature vector is fed to the classifier. There are many
classification models that have been used [3, 4, 5, 6], with support vector
machine (SVM) being one of the most popular choices.

Recently, deep learning methods have been introduced to this field. In [7],
deep neural network (DNN) was used on the top of traditional utterance-level
features and achieved a significant improvement on the accuracy compared with
conventional classifiers. [8] used DNN to learn the short-term acoustic
features, followed by traditional statistical functions to construct
utterance-level features and the extreme learning machine (ELM) was used as
the classifier. In [9], the state-of-the-art result was reported by using both
convolutional and recurrent layers to directly learn the mapping from speech
spectrogram to the corresponding emotion state. In [9], the speech spectrogram
must be segmented into pieces or zero-padded to a fixed size to satisfy the
requirement of convolutional neural network (CNN). Each sub-utterance was
assigned the emotion label of the corresponding whole utterance. And during
the testing procedure, the prediction of the whole utterance was evaluated by
averaging the posterior probabilities of all sub-utterances. However, it is
not quite reasonable to assume that each sub-sentence within a whole sentence
represents the overall emotion. In addition, the speech continuity could be
destroyed by segmentation which made the system more difficult to catch the
whole process of emotion changing from rise to fall.

To solve this problem, in this study, the fully convolutional neural network
(FCN) is adopted to handle variable-length speech, free of the demand of
segmentation to keep critical information not lost. In addition, attention
mechanism has shown its efficiency especially in encoder-decoder models [10,
11, 12], which is employed to guide the decoder to know which parts of the
outputs of the encoder are more important. Specific to classification models,
a self-attention mechanism has been proposed and it is designed to tell the
classifier which parts of the input are more relevant to the output classes.
In [13, 14], the self-attention was used to extract sentence embedding for
semantic analysis. In [15], the authors used the self-attention mechanism on
SER, enabling the network to focus on emotional salient of an utterance. The
encoder they adopted is long short-term memory (LSTM). Considering many
irrelevant signals are mixed with speech signals, we adopt attention mechanism
with FCN to achieve 2D attention visualization on top of spectrograms rather
than 1D attention visualization only on the time axis in [15].

Another problem is the that speech data with emotional labeling is usually
hard to collect. Transfer learning is a useful method to solve the current
task with the help of the knowledge obtained from related problems [10], which
can be operated as finetuning the parameters of network from a pre-trained
model, It has been widely used when the training data is insufficient [16, 17,
18, 19], especially when the model is based on CNN. In this paper, we present
an interesting observation, i.e., an obvious improvement on SER can be
obtained with natural scene image based pre-trained model. It worth noting
that speech signal is very different from image.

The remainder of the paper is organized as follows. In Section 2, we first
introduce the proposed architecture. In Section 3, we report and analyze
experiment results. Finally we summarize our work and present conclusions in
Section 4.

##  2 The Proposed Architecture

Figure 2: The overall architecture of an attention based fully convolutional
neural network.

In this paper, we propose a novel attention based fully convolutional neural
network. The input of the model is also the spectrogram. But inherently unlike
[9], we do not need to segment spectrograms into pieces or pad them to a fixed
shape. The FCN is able to handle spectrogram with variable sizes. The overall
architecture is shown in Figure 2. The FCN encodes the spectrogram into a
high-level representation while the attention mechanism impels the remaining
sub-layers of the model to focus on specific time-frequency regions of the
input spectrogram. All components of the system can be optimized jointly.

Figure 3: The AlexNet based FCN configurations. The convolutional layer
parameters are denoted as “Conv(kernel size)-[stride size]-[number of
channels]”. The maxpooling layer parameters are denoted as “Maxpool-[kernel
size]-[stride size]”. For brevity, the local response normalization layer and
ReLU activation function is not shown.

###  2.1 Fully convolutional network

CNN has been widely used for deep learning, which does not require traditional
handcrafted feature extraction and it has been proved that CNN based system
can obtain a comparable or even better accuracy compared with the traditional
systems on the SER task [15, 20, 21]. The basic components of CNN are
convolution, pooling and activation layers. The convolutional layer is
determined by the number of input channels, the number of output feature maps,
the kernel size and stride. Each kernel can be considered as a filter whose
size is usually much smaller than the input. Hence, a kernel operates on a
local region of input rather than the whole feature map. The locations that
connect to higher layers are called receptive fields. On a given feature map,
the kernel weights are shared to detect certain feature in different locations
and to reduce the complexity of network. The pooling layers usually conduct an
average or max pooling operation to remove noise and extract robust features.
The activation layers are actually element-wise nonlinear functions [12].

The typical CNNs, including AlexNet [22], Oxford VGGNet [23], and ResNet [24]
take fixed-size input. Inspired by [18], we turn the AlexNet into a fully
convolutional network by simply removing its fully connected layers. And then
it is used as our encoder, which is shown in Figure 3. All the convolution
layers are followed by a ReLU activation function, and the first two
convolution layers are equipped with a local response normalization layer. We
also directly adapt the VGGnet to classify emotion states, but it yields a
lower accuracy than AlexNet due to the limited training data.

Assuming that the output of FCN encoder is a 3-dimensional array of size
$F\times T\times C$,where the F and T correspond to the frequency and time
domains of spectrogram and C is channel size. We can consider the output as a
variable-length grid of $L$ elements, $L=F\times T$. Each of the elements is a
C-dimensional vector corresponding to a region of speech spectrogram,
represented as $\boldsymbol{a}_{i}$.

$\boldsymbol{A}=\\{\boldsymbol{a}_{1},\cdots,\boldsymbol{a}_{L}\\},\boldsymbol{a}_{i}\in\mathbb{R}^{C}$
(1)

###  2.2 Attention layer

Intuitively, not all time-frequency units contribute equally to the emotion
state of the whole utterance, i.e., not all the element vectors of set
$\boldsymbol{A}$ contribute equally to the emotion state. Hence, we introduce
attention mechanism to extract the elements that are important to the emotion
of the utterance and aggregate those element arrays to form an utterance
emotion vector. We use the following formulas to realize this idea:

$e_{i}=\boldsymbol{u}^{T}\tanh(\boldsymbol{W}\boldsymbol{a}_{i}+\boldsymbol{b})$
(2) $\alpha_{i}=\frac{\exp(\lambda e_{i})}{\sum_{k=1}^{L}\exp(\lambda e_{k})}$
(3) $\boldsymbol{c}=\sum_{i=1}^{L}\alpha_{i}\boldsymbol{a}_{i}$ (4)

That is, first we feed the annotation $\boldsymbol{a}_{i}$ through a multi-
layer-perceptron (MLP) layer with the $\tanh$ as the non-linear activation
function to obtain a new representation of $\boldsymbol{a}_{i}$. Then we
measure the importance weight, $e_{i}$, of the $\boldsymbol{a}_{i}$ by the
inner product between this new vector and the learnable vector
$\boldsymbol{u}$. After that, the normalized importance weight $\alpha_{i}$ is
calculated through the softmax function. Finally, the utterance emotion vector
$\boldsymbol{c}$ is computed as the weighted sum of set $\boldsymbol{A}$ with
importance weights. $\lambda$ is a scale factor which controls the uniformity
of the importance weights of the annotation vectors. $\lambda$ ranges between
0 and 1. If $\lambda=1$, the scaled-softmax becomes the commonly used softmax
function. If $\lambda=0$, the importance weights will be a uniform
distribution on the set $\boldsymbol{A}$, which means all the time-frequency
units have the same importance weights for the final utterance emotion vector.
In this study, $\lambda=0.3$ is used for the balance.

##  3 Experiments

###  3.1 Database and feature extraction

We validate our systems on the IEMOCAP database [25], one of the widely used
databases on speech emotion recognition. The IEMOCAP corpus comprises five
sessions, each of which includes labeled emotional speech utterances from
recordings of dialogs between two actors. There is no actor overlapping
between these sessions. To be comparable with [9], we utilize the database in
the same way:

  * •

The IEMOCAP database contains scripted and improvised dialogs. We only use
improvised data.

  * •

We use the speech utterances from four emotion categories, i.e., happy, sad,
angry and neutral.

  * •

We implement a five-fold cross validation. In each fold, the data from four
sessions is used for model training, and the data from the remaining session
is splited: one actor for validation and the other one as the testing set.

The experiments only apply the raw spectrogram as the input, the spectrogram
extraction process is consistent with [9]: First, a sequence of overlapping
Hamming windows are applied to the speech waveform, with window shift set to
10 msec, and window size set to 40 msec. Then, for each frame we calculate a
discrete Fourier transform (DFT) of length 800. Finally the 200-dimensional
low-frequency part of the spectrogram is used as the input. Please note that
in [9], the Hamming window size of 20 msec is used, and the authors concluded
that the size of 20 msec is better. In our study, we set the window size to 40
msec and achieve a higher accuracy.

###  3.2 Evaluation metric

The IEMOCAP database is imbalanced with respect to the emotional classes. So
we adopt both the weighted accuracy (WA) and the unweighted accuracy (UA) as
the metric:

  * •

Weighted accuracy - the overall accuracy across all utterances of the testing
set.

  * •

Unweighted accuracy - the average of accuracies across all the classes.

###  3.3 Experiment results and analysis

Table 1: The accuracy comparison of AlexNet and VGGNet-16 with random initialization or finetuning. System |  | Weighted
---
Accuracy
| Unweighted
---
Accuracy
AlexNet Random-init | $66.5\%$ | $54.8\%$
AlexNet Finetuning | $\boldsymbol{67.9\%}$ | $\boldsymbol{57.3\%}$
VGGNet-16 Random-init | $65.3\%$ | $54.8\%$
VGGNet-16 Finetuning | $66.8\%$ | $56.7\%$

First, we directly adapt AlexNet and VGGNet-16 to classify. The only
difference is that there are 4 nodes in softmax layer. The utterances are
split or padded to fixed-length sub-utterances by using the same method in
[9]. During the testing procedure, the posterior probabilities are the average
of the all sub-utterances respectively. Considering the limited data, we
compare the networks with random initialization and the pre-trained networks
based on ImageNet dataset [26].

Table 1 summarizes the results of AlexNet and VGGNet-16 with different
initializations. It’s interesting to observe the pre-trained neural networks
(NNs) always outperform the NNs with random initialization. It’s worth noting
that the speech signal is very different from image. The only explanation is
the pre-trained NNs have been empowered to detect some certain structures so
that they can be more easily trained. By comparing the first row and third row
(or comparing the second row and fourth row), we demonstrate that the AlexNet
outperforms the VGGNet-16 on this task. We think the lack of sufficient
training speech data is one main reason. Based on these results, our FCN model
directly uses AlexNet, excluding its full connected layers. And we initialize
FCN by using the pre-trained parameters.

Table 2: The accuracy comparison between FCN based attention model and the other systems. System |  | Weighted
---
Accuracy
| Unweighted
---
Accuracy
Our FCN+Attention | $\boldsymbol{70.4\%}$ | $\boldsymbol{63.9\%}$
CNN+LSTM Model1 in [9] | $68.8\%$ | $59.4\%$
CNN+LSTM Model2 in [9] | $67.3\%$ | $62.0\%$
Figure 4: The 2D-attention weights of FCN model for 4 examples in different
emotion categories. Top: The spectrogram. Bottom: The 2D-attention weights
figure of spectrogram. Each point of the figure corresponds to the point of
spectrogram in the same location and the brighter color represents the larger
weight.

The published state-of-the-art results using the IEMOCAP corpus are given in
[9]. We list their two best models, i.e., CNN+LSTM Model1 and CNN+LSTM Model2
in Table 2. Model1 is a CNN-LSTM model while Model2 is trained based on Model1
in order to improve the unweighted accuracy. The attention based FCN model is
trained just by one step. And compared with the best results in both Model1
and Model2, our attention based FCN model achieves $1.6\%$ and $1.9\%$
absolute improvements on WA and UA, respectively.

To explain why the improvement can be gained from our proposed approach, we
plot the 2D-attention weights of FCN model for 4 test examples of different
emotion categories in Figure 4. The Figure 4 illustrates that the attention
weights of the non-speech frames are quite small which indicates that the
voice activation detection is implicitly implemented and the information from
non-speech frames are ignored by the attention mechanism automatically.
Besides, the time-frequency units of spectrogram are assigned different
weights based on the degrees they are relevant to emotion states. That
explains why the attention weights are also small on parts of the voice
frames. And the attention weights always have small values in high frequency
areas, which is consistent with the common sense that the information of
speech is mainly contained in the low frequency area. Actually, the bright
area extends from low frequency to high frequency with a decreasing
brightness. This indicates our 2D-attention mechanism has detected the
emotional segment successfully and assigned decreasing weights from low to
high frequency bands. The 2D-attention mechanism is able to scan the
spectrogram not only in the time domain but also in the frequency domain.

##  4 Conclusions

We demonstrated that the CNN architectures designed for visual recognition can
be directly adapted for speech emotion recognition. Besides, it’s interesting
to see the transfer learning can build a solid bridge between natural image
and speech signal. Finally, we proposed an attention based FCN model. Our
model is able to handle utterances with variable lengths and the attention
mechanism empowers the network to focus on emotionally salient regions of
spectrogram. Our system achieves beyond the state-of-the-art accuracy on the
benchmark dataset IEMOCAP.

## References

  * [1] R. Cowie, E. Douglas-Cowie, N. Tsapatsoulis, G. Votsis, S. Kollias, W. Fellenz, and J. G. Taylor, “Emotion recognition in human-computer interaction,” _IEEE Signal processing magazine_ , vol. 18, no. 1, pp. 32–80, 2001.
  * [2] F. Burkhardt, J. Ajmera, R. Englert, J. Stegmann, and W. Burleson, “Detecting anger in automated voice portal dialogs,” in _Ninth International Conference on Spoken Language Processing_ , 2006.
  * [3] C. Vinola and K. Vimaladevi, “A survey on human emotion recognition approaches, databases and applications,” _ELCVIA Electronic Letters on Computer Vision and Image Analysis_ , vol. 14, no. 2, pp. 24–44, 2015.
  * [4] M. El Ayadi, M. S. Kamel, and F. Karray, “Survey on speech emotion recognition: Features, classification schemes, and databases,” _Pattern Recognition_ , vol. 44, no. 3, pp. 572–587, 2011.
  * [5] P. Chandrasekar, S. Chapaneri, and D. Jayaswal, “Automatic speech emotion recognition: A survey,” in _Circuits, Systems, Communication and Information Technology Applications (CSCITA), 2014 International Conference on_. IEEE, 2014, pp. 341–346.
  * [6] S. G. Koolagudi and K. S. Rao, “Emotion recognition from speech: a review,” _International journal of speech technology_ , vol. 15, no. 2, pp. 99–117, 2012.
  * [7] A. Stuhlsatz, C. Meyer, F. Eyben, T. Zielke, G. Meier, and B. Schuller, “Deep neural networks for acoustic emotion recognition: raising the benchmarks,” in _Acoustics, speech and signal processing (ICASSP), 2011 IEEE international conference on_. IEEE, 2011, pp. 5688–5691.
  * [8] K. Han, D. Yu, and I. Tashev, “Speech emotion recognition using deep neural network and extreme learning machine,” in _Fifteenth Annual Conference of the International Speech Communication Association_ , 2014.
  * [9] A. Satt, S. Rozenberg, and R. Hoory, “Efficient emotion recognition from speech using deep learning on spectrograms,” _Proc. Interspeech 2017_ , pp. 1089–1093, 2017.
  * [10] D. Bahdanau, K. Cho, and Y. Bengio, “Neural machine translation by jointly learning to align and translate,” _arXiv preprint arXiv:1409.0473_ , 2014\.
  * [11] M.-T. Luong, H. Pham, and C. D. Manning, “Effective approaches to attention-based neural machine translation,” _arXiv preprint arXiv:1508.04025_ , 2015.
  * [12] J. Zhang, J. Du, S. Zhang, D. Liu, Y. Hu, J. Hu, S. Wei, and L. Dai, “Watch, attend and parse: An end-to-end neural network based approach to handwritten mathematical expression recognition,” _Pattern Recognition_ , vol. 71, pp. 196–206, 2017.
  * [13] Z. Yang, D. Yang, C. Dyer, X. He, A. Smola, and E. Hovy, “Hierarchical attention networks for document classification,” in _Proceedings of the 2016 Conference of the North American Chapter of the Association for Computational Linguistics: Human Language Technologies_ , 2016, pp. 1480–1489.
  * [14] Z. Lin, M. Feng, C. N. d. Santos, M. Yu, B. Xiang, B. Zhou, and Y. Bengio, “A structured self-attentive sentence embedding,” _arXiv preprint arXiv:1703.03130_ , 2017.
  * [15] S. Mirsamadi, E. Barsoum, and C. Zhang, “Automatic speech emotion recognition using recurrent neural networks with local attention,” in _Acoustics, Speech and Signal Processing (ICASSP), 2017 IEEE International Conference on_. IEEE, 2017, pp. 2227–2231.
  * [16] R. Girshick, J. Donahue, T. Darrell, and J. Malik, “Rich feature hierarchies for accurate object detection and semantic segmentation,” in _Proceedings of the IEEE conference on computer vision and pattern recognition_ , 2014, pp. 580–587.
  * [17] J. Donahue, Y. Jia, O. Vinyals, J. Hoffman, N. Zhang, E. Tzeng, and T. Darrell, “Decaf: A deep convolutional activation feature for generic visual recognition,” in _International conference on machine learning_ , 2014, pp. 647–655.
  * [18] J. Long, E. Shelhamer, and T. Darrell, “Fully convolutional networks for semantic segmentation,” in _Proceedings of the IEEE conference on computer vision and pattern recognition_ , 2015, pp. 3431–3440.
  * [19] O. Koller, O. Zargaran, H. Ney, and R. Bowden, “Deep sign: hybrid cnn-hmm for continuous sign language recognition,” in _Proceedings of the British Machine Vision Conference 2016_ , 2016.
  * [20] G. Trigeorgis, F. Ringeval, R. Brueckner, E. Marchi, M. A. Nicolaou, B. Schuller, and S. Zafeiriou, “Adieu features? end-to-end speech emotion recognition using a deep convolutional recurrent network,” in _Acoustics, Speech and Signal Processing (ICASSP), 2016 IEEE International Conference on_. IEEE, 2016, pp. 5200–5204.
  * [21] Z. Aldeneh and E. M. Provost, “Using regional saliency for speech emotion recognition,” in _Acoustics, Speech and Signal Processing (ICASSP), 2017 IEEE International Conference on_. IEEE, 2017, pp. 2741–2745.
  * [22] A. Krizhevsky, I. Sutskever, and G. E. Hinton, “Imagenet classification with deep convolutional neural networks,” in _Advances in neural information processing systems_ , 2012, pp. 1097–1105.
  * [23] K. Simonyan and A. Zisserman, “Very deep convolutional networks for large-scale image recognition,” _arXiv preprint arXiv:1409.1556_ , 2014.
  * [24] K. He, X. Zhang, S. Ren, and J. Sun, “Deep residual learning for image recognition,” in _Proceedings of the IEEE conference on computer vision and pattern recognition_ , 2016, pp. 770–778.
  * [25] C. Busso, M. Bulut, C.-C. Lee, A. Kazemzadeh, E. Mower, S. Kim, J. N. Chang, S. Lee, and S. S. Narayanan, “Iemocap: Interactive emotional dyadic motion capture database,” _Language resources and evaluation_ , vol. 42, no. 4, p. 335, 2008.
  * [26] O. Russakovsky, J. Deng, H. Su, J. Krause, S. Satheesh, S. Ma, Z. Huang, A. Karpathy, A. Khosla, M. Bernstein _et al._ , “Imagenet large scale visual recognition challenge,” _International Journal of Computer Vision_ , vol. 115, no. 3, pp. 211–252, 2015.