1 Introduction
Animal classification systems are based on evolutionary relationships between different organisms, known as phylogeny. This approach allows us to organize species in a way that reflects our understanding of how they evolved from common ancestors. Phylogenetic trees are diagrams that graphically represent these evolutionary relationships between organisms. In these representations, the species of interest are placed at the tips of branches that emerge from a point representing a common ancestor. The natural mathematical object associated with this situation is a tree (a graph without cycles) with a root (the common ancestor). It is important to note that the hypotheses regarding how different species may have descended from a common ancestor are typically based on physical traits (which are therefore interpretable) or directly on DNA sequences.
In this context, one can attempt to infer phylogenetic distances based on phenotypic traits. For instance, in the case of birds, one might consider the possibility of inferring these distances based on their songs [Reference Beecher and Brenowitz2, Reference Cate8, Reference Chen and Wiens10, Reference Martens27, Reference Mason, Burns, Tobias, Claramunt, Seddon and Derryberry28, Reference Price and Lanyon37, Reference Tachibana, Oosugi and Okanoya42, Reference Wimberger and de Queiroz45]. If that is the case, what set of acoustic properties would allow us to predict phylogenetic distances? In this study, we propose a method based on dimensionality reduction (embeddings) and tree reconstruction to process sonograms (that is, graphical representations of bird songs) to test the hypothesis that it is possible to predict phylogenetic distances from bird songs without predefining the relevant properties.
On one hand, we address the general problem of reconstructing latent trees from images, and then we apply this method to a recently discussed example in biology, where pre-defined acoustic properties were used in conjunction with AI techniques to infer phylogenetic distances among a set of songbird species [Reference Moreira, Marques, Costeira and Hauptmann32, Reference Rivera, Edwards, Hauber and Woolley39].
Our approach consists of two main components:
-
(1) finding an embedding of the sonograms in a space where we can compute meaningful distances;
-
(2) reconstructing the tree based on these distances.
Each of the components is based on well-known methods: Poincaré embeddings [Reference Nickel and Kiela34] to obtain the desired distances and the neighbour-joining algorithm [Reference Compeau and Pevzner13, Reference Saitou and Nei40] to reconstruct the tree. Poincaré embeddings, which map data into hyperbolic space, have shown promise in representing hierarchical structures and obtaining simultaneously distances in the embedded spaces that reflect similarities. Hyperbolic geometry allows for the efficient representation of tree-like structures (see [Reference Atigh, Keller-Ressel and Mettes1, Reference Boguna, Krioukov and Klaffy4, Reference Chami, Ying, Ré and Leskovec9, Reference Ge, Mishra, Kornblith, Li and Jacobs15, Reference GhadimiAtigh, Schoep, Acar, van Noord and Mettes16, Reference Guo, Wang, Chen and Yu18, Reference Khrulkov, Mirvakhabova, Ustinova, Oseledets and Lempitsky20, Reference Krioukov, Papadopoulos, Kitsak, Vahdat and Boguna22–Reference Liu, Chen, Pan, Ngo, Chua and Jiang25, Reference Mathieu, Lan, Maddison, Tomioka and Teh29, Reference Mikolov, Chen, Corrado and Dean30, Reference Moreira, Marques, Costeira and Hauptmann32, Reference Nagano, Yamaguchi, Fujita and Koyama33, Reference Sala, De Sa, Gu and Ré41, Reference Yue, Lin, Yamada and Zhang47] and, in particular, the introduction in [Reference Moreira, Marques, Costeira and Hauptmann32] for an account of the state of the art), making it particularly suitable for our task. The neighbour-joining algorithm, a well-established method in phylogenetics, is then applied to the embedded data to reconstruct the tree.
Unlike previous work in hierarchical reconstruction from embeddings [Reference Nickel and Kiela34], our method deals with the challenge of having only samples of leaves from the latent tree. To address this, we employ Siamese networks, which are particularly well suited for learning similarity metrics between pairs of samples. This approach differs from the loss function considered in [Reference Nickel and Kiela34], which assumes access to nodes at all the levels of the latent hierarchy during training. Siamese networks have been previously used to embed sonograms in Euclidean space [Reference Bistel, Martinez and Mindlin3].
We demonstrate the effectiveness of our approach on both synthetic data with known hierarchical structure and real-world bird song spectrograms. For the latter, we focus on six species of estrildid finches: Stizoptera bichenovii, Poephila acuticauda, Lagonosticta senegala, Lonchura striata, Uraeginthus bengalus and Amandava subflava. These species were chosen based on a recent study comparing their phylogenetic organization with one emerging from pre-selected acoustic properties of their songs [Reference Rivera, Edwards, Hauber and Woolley39]. Our method offers several advantages over traditional approaches.
-
(1) It does not require pre-specification of acoustic or other types of properties, instead learning relevant features directly from the spectrograms.
-
(2) It provides a general framework for reconstructing trees from high-dimensional representations of their leaves, applicable beyond bird song analysis.
-
(3) The use of hyperbolic geometry allows for more efficient representation of hierarchical structures compared with Euclidean embeddings.
The rest of this paper is organized as follows. Section 2 provides the background on hyperbolic geometry and Poincaré embeddings. Section 3 describes our Siamese network architecture for learning embeddings. Section 4 gives details of the neighbour-joining algorithm and our tree reconstruction process. Section 5 presents our results on both synthetic and bird song data, comparing our reconstructed trees to known phylogenetic relationships. Finally, Section 6 discusses the implications of our findings and potential directions for future research.
This work contributes to the growing body of research at the intersection of artificial intelligence, machine learning and evolutionary biology, offering a novel method for inferring evolutionary relationships from complex, high-dimensional data such as spectrograms. By bridging the gap between behavioural traits and phylogenetic reconstruction, our approach has the potential to provide new insights into the evolution of animal communication and behaviour.
2 Poincaré embeddings
2.1 Hyperbolic geometry and the Poincaré ball model
Hyperbolic geometry is a non-Euclidean geometry in which the parallel postulate is replaced with the following one.
For any given line R and point P not on R, in the plane containing both line R and point P, there are at least two distinct lines through P that do not intersect R.
Hyperbolic spaces have begun to receive attention in the machine learning community, as they are suitable for modelling data with hierarchical latent structure. An introduction to the topic can be found in [Reference Cannon, Floyd, Kenyon and Parry7, Reference Greenberg17] with applications discussed in works such as [Reference Nickel and Kiela34, Reference Yue, Lin, Yamada and Zhang47]. Heuristically, they can be thought of as a continuous version of trees.
The Poincaré disk (ball) model is a hyperbolic space that uses the interior of a Euclidean unit disk (ball) as the representation space,
$\mathcal D = \{(x,y) \mid x^2+y^2 < 1\}$
. The boundary
$\partial \mathcal D = \{(x,y)\mid x^2 + y^2 = 1\}$
of
$\mathcal D$
does not belong to the disk. The distance between points
$u, v \in \mathcal D$
is defined by
$$ \begin{align} d (u, v) = \textrm{arcosh} \bigg( 1 + \frac{2 \|u - v\|^2}{(1 - \|u\|^2)(1 - \|v\|^2)}\bigg). \end{align} $$
Here,
$\|\cdot \|$
specifically refers to the Euclidean norm. This formula is valid for points
$u, v$
within the domain
$\mathcal {D}$
, where
$\|u\| < 1$
and
$\|v\| < 1$
. The d-dimensional Poincaré ball model is obtained similarly by replacing the two-dimensional disk with a d-dimensional ball. The Poincaré ball is a metric space and, moreover, it can be equipped with a Riemannian metric tensor to obtain a Riemannian manifold, but we will not require this.
When u is the centre of the circle and v is at Euclidean distance r from the origin, we get
$$ \begin{align*} d(0,v) = \textrm{arcosh} \bigg({\frac {1+r^2}{1-r^2}}\bigg) = \ln \bigg({\frac {(1+r)^2}{1-r^2}}\bigg)=2 \textrm{arctanh}( r). \end{align*} $$
Note that the distance grows unbounded as v approaches to the boundary and, consequently, the hyperbolic distance from a point in the interior of
$\partial \mathcal D$
to the boundary of the disk is infinite.
2.2 Poincaré embeddings
Adequate embeddings (representations) of symbolic data in different types of spaces and by different means have become a key point for many learning tasks. The most prominent examples are word embedding and their use to train Large Language Models [Reference Bojanowski, Grave, Joulin and Mikolov5, Reference Kenton Jacob Devlin and Toutanova19, Reference Liu, Lin, Shi and Zhao26, Reference Mikolov, Chen, Corrado and Dean31, Reference Pennington, Socher and Manning35, Reference Peters, Neumann, Iyyer, Gardner, Clark, Lee and Zettlemoyer36, Reference Radford, Narasimhan, Salimans and Sutskever38, Reference Vaswani, Shazeer, Parmar, Uszkoreit, Jones, Gomez, Kaiser and Polosukhin44].
In [Reference Nickel and Kiela34], the authors suggested to use the Poincaré ball as a suitable space to look for a good embedding for data with a latent hierarchical structure.
If we consider data with a hierarchy like a tree, then the Euclidean space is not ideal. For example, a tree with a branching factor b has
$b^{\ell }$
nodes at level
$\ell $
and
$\{(b^{\ell +1} -1)/({b-1})\}$
nodes at levels less than or equal to
$\ell $
. Thus, as the levels of the tree grow, the number of nodes grows exponentially. If we want to embed such a tree in Euclidean space in such a way that the distances between nodes is adequately represented by Euclidean distances, we need the dimension of the space to be as large as the number of levels in the tree minus one.
Remarkably, we can represent our tree within two-dimensional hyperbolic space by means of the following algorithm: we place the root of the tree at the centre and their children equidistantly. Then, we repeat this procedure recursively in each branch. In Figure 1, we can visualize this embedding. In this figure, all points are equidistant from their parents. This type of construction is possible because in hyperbolic space, we can accommodate as many equidistant points as desired by settling them close enough to the boundary.
Embedding of a tree with branching factor 2 into the Poincaré disk. Image is taken from [Reference Nickel and Kiela34].

Thus, if the root of a tree is placed at the origin, the relative distance to other nodes is small. However, the leaves can be located near the boundary to guarantee that the hyperbolic distance coincides with the graph distance in the tree.
In view of this, we are going to look for an embedding of sonograms from bird songs into the Poincaré disk (ball) with the hope that such embeddings can give information about a latent hierarchical structure determined by phylogenetics. Previously, to address the performance of the method, we will deal with synthetic data composed of images with a latent hierarchical structure.
The goal of computing embeddings for images coming from both synthetic data and spectrograms is to obtain a representation in a space in which distances between points are meaningful. They simplify the representation of images by reducing dimensionality drastically while preserving relationships between them.
The Poincaré ball is used as the embedding space for both our synthetic data and spectrograms data because they both have the desired hierarchical structure and are able to accommodate easily such hierarchy.
In the synthetic data, they represent the leaves of a tree constructed by branching areas of a square into zeros or ones. Spectrograms are the songs of different species that can be related through a phylogenetic tree. Remember that the hypothesis is that spectrograms are the leaves of the evolutionary tree.
We note that our inspiration for using Poincaré embeddings comes from [Reference Nickel and Kiela34], but both our context and, consequently, our algorithm differ. In [Reference Nickel and Kiela34], the algorithm relies solely on information about hyponymy relations between pairs of symbols. Additionally, it is implicitly assumed that the data conform to the Poincaré disk, with data present at all levels of the hierarchy. In our case, given a pair of sonograms, in addition to information on whether they come from the same species, we have a (high-dimensional and poor) notion of similarity between the sonograms themselves. Also, our data correspond only to leaves of the tree, so that we do not expect to fulfil the whole disk but to be close to the boundary.
For these reasons, this approach differs from that in [Reference Nickel and Kiela34]. In that case, embeddings were obtained by accessing nodes at all levels of the latent hierarchy. However, in our case, we only had access to the leaves. Consequently, we needed to use a tailored loss function during training to address this difference.
3 Siamese networks
A Siamese network is a type of neural network that consists of two or more identical sub-networks with shared weights and the same architecture. These sub-networks process different inputs independently, and their outputs are combined through a separate layer to produce a final result. Effectively, the network
$G_{w}$
is the same for all inputs, but it is evaluated multiple times with different input values. They are used to find the similarity between two comparable objects. Among other tasks, they have been used successfully for face recognition and signature verification [Reference Bromley, Guyon, LeCun, Säckinger and Shah6, Reference Chicco11, Reference Chopra, Hadsell and LeCun12, Reference Taigman, Yang, Ranzato and Wolf43].
The architecture of a Siamese network typically consists of a convolutional backbone followed by a dense layer. The input images are passed through the backbone and the output features are passed to a dense layer, which calculates a vector of dimension n. This vector is then passed to the output layer, which calculates the Euclidean distance between the feature vectors of the input images.
A typical loss function used for training a Siamese network is the contrastive loss function,
$$ \begin{align} L(w \mid x, y) &= \sum_{(i,j) \in \mathcal T}\frac{(1 - y_{ij}) }{2} d_{ij}^{2} + \frac{y_{ij}}{2} [ \max(0, m - d_{ij}) ]^{2} \nonumber\\ &\text{with} \quad d_{ij} = d(G_{w}(x_{i}), G_{w}(x_{j})), \end{align} $$
where
$y_{ij}$
is the label indicating whether the images
$x_i$
and
$x_j$
are similar (
$y_{ij} = 0$
) or dissimilar (
$y_{ij} = 1$
), and
$d_{ij}$
is the distance between the feature vectors of the input images, which is usually an
$L^p$
-norm in Euclidean space. The training set
$\mathcal T$
is a suitably chosen subset of the product space of all possible pairs. The variable w stands for the weights of the network.
The term
${(1 - y) } d^{2}/2$
penalizes similar images that have a large distance between them, and
${y}/{2} [ \max (0, m - d) ]^{2}$
penalizes dissimilar images whose distance is less than a margin m.
Siamese networks are trained with pairs of similar and dissimilar images that constitute the set
$\mathcal T$
. Their output can be used to measure the similarity between two images.
Since we are looking for an embedding in hyperbolic space, this work proposes replacing the standard choice of
$d_{ij}$
with the distance in the Poincaré ball. While such distances have been used in other contexts, to the best of our knowledge, they have not been applied to generating embeddings for bird songs.
4 Tree reconstruction
In this section, we review the neighbour-joining algorithm for tree reconstruction. The input of the algorithm is a set of (compatible) distances between leaves of the (latent) tree. The output is a tree in which the distance between leaves are those given by the input.
4.1 Phylogenetic tree
Phylogenetic trees allow visualizing the evolutionary development of different species, subspecies or populations called taxa, and the relationships among them. Each of the leaves of these trees represents different taxa and the interior vertices represent common ancestors from which they derive. If there is a common ancestor to all taxa, it is located at the root of the tree.
The most recent taxa are represented as leaves of the trees. Given a leaf j, the only node to which it is connected is called the parent of j. Also, the edge connecting a leaf to its parent is called a branch.
In the reconstruction of phylogenetic trees, the aim is typically to relate taxa that are currently present. Information about these is available, but not about those represented by internal nodes. The goal is to reconstruct evolutionary relationships from the distances between the terminal taxa. Distances information can be obtained in various ways. In the case of the bird species analysed in this work, distance information arises from differences and similarities between their songs and is measured through the embedding of the sonograms in the Poincaré ball.
4.2 Distance matrices
A matrix
$D \in \mathbb R^{n \times n}$
with nonnegative entries is used to represent the distances between different taxa. For such a matrix D, we will say that it is a distance matrix if it is symmetric and satisfies the triangular inequality,
$D_{i,k} \leq D_{i,j} + D_{j,k}$
for all
$1\le i,j,k \le n$
.
To represent the evolutionary distance between taxa, a length is assigned to each edge. The length of a path within the tree is defined as the sum of the distances of each edge. Thus, for a tree T, the evolutionary distance between two present taxa i and j, called
$d_{i,j}(T)$
, corresponds to the length of the shortest path from i to j.
4.3 Neighbour-joining algorithm
One of the most commonly used methods for reconstructing evolutionary trees is the neighbour-joining (NJ) algorithm proposed by Saitou and Masatoshi [Reference Saitou and Nei40]. We briefly review the algorithm following [Reference Compeau and Pevzner13]. Given a distance matrix D that represents distances between leaves of an unknown tree, the algorithm selects pairs of leaves that are inferred to be neighbours and substitutes them with a single leaf, thus reconstructing a tree by repeating the process recursively. A natural way to do this is to choose a pair of leaves with minimum distance and to assume that they have a common parent. However, finding a minimum in the distance matrix does not guarantee finding a pair of siblings in the tree. Hopefully, it is possible to transform this matrix into another one whose minimum element is indeed a pair of siblings. Proceeding recursively, we can reconstruct the whole tree.
The NJ algorithm produces an unrooted phylogenetic tree. However, knowing the root is of great importance for interpreting evolutionary relationships, as it represents the common ancestor to all taxa of the tree and provides the tree with an evolutionary path of the species.
There are several methods to place the root of a tree, the most used one being outgroup rooting. In this case, we need external (biological) knowledge to identify subgroups in the tree. Since we want to have an algorithm which is free of external knowledge, we opt to use midpoint rooting [Reference Farris14], which is only based on the distances calculated by the tree reconstruction algorithm.
5 Results
In this section, we show how to apply our method. First, we will work with synthetic data, generated with a controlled “evolutionary” algorithm. Then, we will apply our method to sonograms computed from the public data base Xeno-canto [46]. Our method is implemented with the following algorithm.

5.1 Synthetic data
We generate synthetic images with a hierarchical tree structure. The algorithm to produce the images is as follows: we split the image in two equal parts. The left half of the image is coloured either black or white representing a first branch that splits into a zero or a one. Then, we split the right half of the image into two equal parts, the upper half is coloured either white or black. This defines a new branching of the tree into zero or one. One quarter of the image is still left undecided about its colour at this stage. In the next step, the left half of the lower right quarter takes either white or black colour. Repeating the procedure for five steps, with the sector of the image chosen in a spiral manner at each step, different classes representing the leaves of a tree are obtained. After five branchings, the remaining area of the image is coloured in black and white in equal proportions. The set of all images that can be produced with this algorithm can be thought of as leaves of a tree. In this way, images representing different classes are formed. To generate distinct elements for each class, noise is added to the images by altering the colour of 10
$\%$
of the pixels, changing them to various shades of grey.
Figure 2 shows six leaves chosen at random. For those leaves, we generate 80 images to train the model and 20 for validation. The tree structure representing them is shown in Figure 3. The tree in panel (b) is the desired result to be reconstructed.
Six images with implicit hierarchy.

Implicit hierarchical structure of the synthetic images.

The model, which aims to find embeddings to represent the data is a Siamese neural network. Each of the twin networks is itself a convolutional net. Using the notation
$N @ k \times k$
for a convolutional layer with N filters and a kernel of size
$k \times k$
, the structure of the convolutional network is as follows:
-
• Convolution + ReLU
$64 @ 10 \times 10$
; -
• Max Pool
$2 \times 2$
; -
• Convolution + ReLU
$128 @ 7\times 7$
; -
• Max Pool
$2\times 2$
; -
• Convolution + ReLU
$256 @ 4 \times 4$
; -
• Max Pool
$2 \times 2$
; -
• Convolution + ReLU
$512 @ 4 \times 4$
; -
• Max Pool
$2 \times 2$
; -
• Dense Layer 512 + ReLU;
-
• Dense Layer 2 + Linear;
-
•
$x^* \gets 0.99 \tanh (|x|) {x}/{|x|}$
.
Thus, the Siamese neural network is defined using two convolutional networks with ReLU activation function; each of which takes an image and produces a vector within the Poincaré disk in two dimensions. The network follows the typical structure of Siamese networks, after each convolutional layer, which doubles the number of filters of the previous one, a pooling layer is placed to reduce the dimension.
To compute the hyperbolic distance between these vectors using the Poincaré distance defined in (2.1), we need the output vectors to belong to the (open) ball of radius one; then we apply the following transformation to the vector x obtained as the two-dimensional result of the last dense layer:
$x^* = .99 \tanh (|x|) {x}/{|x|}$
.
Also, by using the
$\tanh $
function, points close to the centre remain close and those with large norm are closer to the boundary of the ball. The new point
$x^*$
maintains the direction of x. Once both vectors are obtained within the ball of radius one, it is possible to calculate the Poincaré distance between them, (see Figure 4).
Siamese network architecture. The distance between the outputs of each of the networks is computed in the Poincaré ball. The loss function imposes the embedding of different classes to be separated in the Poincaré ball while embeddings of elements of the same class are forced to be close in the Poincaré ball.

Learning is done with a series of mini-batches of 32 pairs of images. The choice of these pairs is not completely random: half of them are positive, meaning that they are chosen so that they belong to the same class, and the other half are negative, meaning that they belong to different classes. The model is trained for 150 epochs or iterations. The parameter m for the contrastive loss is set to
$m=25$
. Adam optimizer [Reference Kingma and Ba21] is used, along with a learning rate parameter of
$10^{-6}$
.
We compute the embedding ten times. For each of them, the initialization of the network weights is varied following a Gaussian distribution. To measure the performance of the neural network, 10 classes are randomly chosen, allowing repetition. For each class, one image is selected, called the support image. Then, a test image is chosen, and to compare it with the support images, the convolutional model is applied to each of them. If the Poincaré distance between the test image and the support image of the same class is the lowest, then we say that the task is successful. Repeating the process 100 times, we obtain a measure of the accuracy of the model as the average of the successful tasks.
The embeddings obtained for this data can be seen in Figure 5. The neural network successfully separates the different classes and achieved
$98.2\%$
accuracy in training and
$97\%$
in validation. The results for each run are shown in Table 1.
Accuracy of the model applied to synthetic data.

Embeddings in the Poincaré disk. Top: runs 1, 2 and 7. Bottom: runs 8, 9 and 10.

From the representation of the data in this space, the Chebyshev centre of each class is computed using the Poincaré distances. Thus, each class has a representative in the embedding space. Using these representatives, the Poincaré distance matrix between the classes is computed. Then, we apply the NJ algorithm to reconstruct the trees. The result is that in
$90\%$
of the runs, the tree in Figure 6(a) is obtained; therefore, each pair of images is correctly identified. The algorithm is robust for this dataset. After applying midpoint rooting to choose the root, in
$70\%$
of the runs, the tree in panel (b) is obtained. Therefore,
$70\%$
of the times, the entire hierarchical structure of Figure 3 is recovered.
Reconstructed trees with NJ. (a) Tree without root obtained in
$90\%$
of the runs; (b) rooted tree (midpoint-rooting) obtained in
$70\%$
of the runs.

5.2 Sonograms
We work with sonograms of syllables from the songs of each of the following species: Stizoptera bichenovii (ST), Poephila acuticauda (PO), Lagonosticta senegala (LA), Lonchura striata (LO), Amandava subflava (AM), Uraeginthus bengalus (UR). The songs from which the representative syllables are extracted have been taken from the Xeno-canto database [46]. We selected these species because a recent study compared their phylogenetic organization with that emerging from comparing the songs by pre-selecting certain acoustic properties [Reference Rivera, Edwards, Hauber and Woolley39]. The method we present here does not require a pre-specification of properties to carry out the reconstruction.
The complete set of audio files is first converted to mono and down-sampled to 22 050 Hz. A 12th-order Butterworth FIR-type band-pass filter is applied with cut-off frequencies between 1 and 8 kHz. Noise reduction is then performed using spectral gating, a technique that applies a noise gate by setting a threshold for each frequency band, suppressing noise components below that threshold. To implement this, we estimate stationary noise thresholds using a signal and noise threshold set at 1.5 standard deviations above the mean, followed by noise suppression at a reduction ratio of 95%.
Next, the audio files are segmented into 1-second chunks with a 0.5-second overlap, padding with zeros when the chunk is shorter than 1 second, but discarding chunks that required more than 0.25 seconds of padding. Spectrograms for each chunk are generated using the short-time Fourier transform (STFT) with a Hamming window. The window size and hop-size are adapted to obtain spectrograms with a shape of (128, 256). A clipping operation is applied to limit the values to a range between 1.5 times the median (values below this are set to zero) and the maximum. The logarithm of the values is then computed for each element in the spectrogram matrix. Finally, the spectrograms are rescaled to values between 0 and 1. The resulting spectrograms are exported as images at 50 dpi resolution for use as inputs to the Siamese networks.
Although the files from Xeno-canto are labelled with the presence of a particular species, the audio recordings contain a variety of noise and vocalizations, which may belong to the labelled species or to others present in the recordings. Consequently, the data are considered weakly labelled. The challenge is to extract the specific vocalizations corresponding to the labelled species. To address this, we manually select audio segments that contained vocalizations of the species of interest. We cross-referenced these selections with audiovisual material from YouTube to increase confidence in the species-specific vocalizations and identified recurring patterns in the Xeno-canto files. We then select audio segments corresponding to one unique syllable type for each species, both from the Xeno-canto database and the YouTube recordings, yielding a total of 568 segments (poephila + acuticauda, 51; stizoptera + bichenovii, 133; lonchura + striata, 159; lagonosticta + senegala, 69; uraeginthus + bengalus, 69; amandava + subflava, 87). These selected segments are used for training and evaluating the network, chosen in advance to avoid data leakage.
Sonograms of the syllables are obtained for each of these species. An example from each class is shown in Figure 7. The dataset consists of 568 images, 80
$\%$
of them are used to train the network and the remaining 20
$\%$
for validation. Additionally, the data are augmented by applying transformations: rotations, translations and zoom.
Six sonograms of representatives of different bird species.

The structure of the convolutional network used to obtain the embeddings is similar to that used with the synthetic data. The only difference is that the final output layer has three dimensions instead of two. The following layers outline the architecture of the network, showing the sequence of operations applied to the input data.
-
• Convolution + ReLU
$64 @ 10 \times 10$
; -
• Max Pool
$2 \times 2$
; -
• Convolution + ReLU
$128 @ 7\times 7$
; -
• Max Pool
$2\times 2$
; -
• Convolution + ReLU
$256 @ 4 \times 4$
; -
• Max Pool
$2 \times 2$
; -
• Convolution + ReLU
$512 @ 4 \times 4$
; -
• Max Pool
$2 \times 2$
; -
• Dense Layer 512 + ReLU;
-
• Dense Layer 3 + Linear;
-
•
$x^* \gets 0.99 \tanh (|x|) {x}/{|x|}$
.
The Siamese neural network uses two convolutional networks, each of which takes a sonogram and produces a point
$x^*$
within the Poincaré ball in dimension three.
Similar to the synthetic data, learning is performed with a series of mini-batches of 32 pairs of images where half of them correspond to sonograms of the same class and the other half to sonograms of different classes. The model required 2500 epochs to converge.
The average result obtained is
$88.3\%$
accuracy for the training set and
$86.3\%$
for the validation set. In Table 2, the accuracy values for each of the ten runs are shown.
Accuracy for the model applied to sonogram data.

The embeddings obtained in the Poincaré space in the first run can be observed in Figure 8.
Sonograms embedding.

Based on the data representation in this space, we compute the Chebyshev centre of each class with respect to the Poincaré distance. Using this matrix of Poincaré distances between the centres, the NJ algorithm is applied.
The unrooted trees obtained from the ten runs can be seen in Figure 9. Subsequently, each of the trees is rooted using midpoint rooting. The results can be found in Figure 10. In
$60\%$
of the runs, the same tree (b) is obtained. Its structure presents similarities with the phylogenetic tree. First, Poephila acuticauda and Stizoptera bichenovii are neighbours, as in the phylogenetic tree. Amandava subflava and Lagonosticta senegala derive from the same common ancestor for both trees. The main difference between these two trees lies in the placement of Lonchura striata and Uraeginthus bengalus. The former shares a parent with Lagonosticta senegala in the song tree, while in the phylogenetic tree, the closest common ancestor is the parent of Poephila acuticauda and Stizoptera bichenovii. In contrast, for the song tree, Uraeginthus bengalus occupies the place that Lonchura striata has in the phylogenetic tree.
Trees generated by NJ: (a) unrooted tree generated from sonograms for 70% of the runs; (b) unrooted tree generated from sonograms for 20% of the runs; (c) unrooted tree generated from sonograms for 10% of the runs.

(a) Phylogenetic tree of bird species. (b)–(d) Trees with roots using midpoint rooting and NJ. The tree in panel (b) was obtained in 60% of the runs. The tree in panel (c) was obtained in 30% of the runs. The tree in panel (d) was obtained in 10% of the runs.

In [Reference Rivera, Edwards, Hauber and Woolley39], the authors construct dendrograms based on clustering with
$k=3$
using acoustic features, with the most effective dendrogram based on the primary principal component derived from frequency-related traits. This approach results in a dendrogram that is aligned with the phylogenetic tree. In contrast, although our method is not as successful at reconstructing the phylogenetic tree, it offers the advantage of not relying on pre-selected features or specifying parameters like the number of clusters.
It is noteworthy that the phylogenetic tree differs from the tree generated from the behaviour in the position of a species: Uraenginthus bengalus. This species has some characteristics that warrant a review of the data used. In particular, the female sings, although with a less complex song than the male. There is a possibility that our dataset included examples of female songs for this species, which would compromise the classification. It would be interesting to assess this hypothesis.
The remaining two song trees have structures with more differences compared with the phylogenetic one. In all runs, Poephila acuticauda and Stizoptera bichenovii are neighbours, and Lagonosticta senegala and Lonchura striata are also neighbours.
6 Conclusion
In recent years, we have witnessed the impressive growth of a new paradigm in science. One that, in some way, contrasts with the “interpretative” attitude towards phenomena that has been characteristic of natural sciences in the last centuries. This new approach, based on data, establishes complex correlations between data without aiming to elucidate mechanisms or employ easily interpretable variables. In the context of this work, we explore what kind of phylogenetic inferences can be drawn from data without the need for prior interpretative processing of the characteristics from which the task is performed. We discuss two examples of how to unveil phylogenia from data conveyed as images. One where the “evolutionary” process is simulated computationally, and another where the data come from the natural world (songs produced by individuals of closely related species). We take the second example, because neural networks were recently employed to make phylogenetic inferences from birdsong. Opposite to our case, in that case, the variables used are “interpretable”. In fact, they are intuitive acoustic features easily recognizable in the song spectrograms. This leads us to the second example.
The first example is synthetic and is used to test the method since in this case, the “phylogeny” is known in advance. We observe that the algorithm is capable to find the hierarchical structure in the data in addition to having only information about the leaves of the tree. That is, the algorithm reconstructs the sub-tree determined by the leaves for which we have examples in the dataset. The only input required by the algorithm is images that encode, in some way, relevant characteristics to determine the hierarchy. We remark that we do not need to know what those characteristics are. They just need to be present in the data. In view of this, we can expect the algorithm to be able to reconstruct, in real situations, phylogenetic trees departing from images that contain relevant phylogenetic characteristics. Moreover, this leads us to the natural problem of inferring if a set of characteristics contains phylogenetic information or not and even more, if there is an additional hierarchical structure besides the phylogeny or if it is the case that it is just that.
The second example is very subtle. Songs, in the species that are studied in this work, are learned. Therefore, some features are genetically conditioned, while others are learned from other individuals. Then, choosing the correct interpretable variables is challenging: we have to choose those that are more likely to be genetically programmed. A data-driven approach might identify the similarities without predetermination of the features.
In this work, we train a Siamese neural network with a metric which is designed to unveil hierarchies. The result is that, statistically, it is capable to infer the tree structure behind the species that is pretty similar to the phylogenetic one.
In any data-driven study, data preparation is a key aspect of the process and might induce specific biases. In our application, the selection of audio segments was performed manually, using a single syllable type as the representative song pattern for each species. Yet, the Siamese network demonstrated strong performance in species recognition, yielding nontrivial results. The generated trees exhibited consistent robustness across multiple runs, further underscoring the method’s effectiveness. Incorporating laboratory-controlled data or automated methods for labelling audio segments could expand the applicability of this approach.
Future work could investigate how the method performs on datasets with a large number of species or more complex acoustic data, as well as its scalability to other phenotypic traits or species. Notably, the computational cost associated with Poincaré embeddings and Siamese network training is a key consideration. Training Siamese networks involves generating and processing input pairs, which significantly increases the dataset size compared with standard supervised learning. This makes scalability to large datasets or extensive feature sets an important area for further exploration.
Nevertheless, the current results highlight both the potential and robustness of the method, as it was evaluated on data that were selected prior to model training, ensuring it remained agnostic and unaffected by data leakage.









