# Barycentric subspace analysis of network-valued data

Elodie Maignant, Xavier Pennec, Alain Trouvé, Anna Calissano

► To cite this version:

Elodie Maignant, Xavier Pennec, Alain Trouvé, Anna Calissano. Barycentric subspace analysis of network-valued data. 2025. hal-05191923

HAL Id: hal-05191923

https://hal.science/hal-05191923v1

Preprint submitted on 30 Jul 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# BARYCENTRIC SUBSPACE ANALYSIS OF NETWORK-VALUED DATA

BY ELODIE MAIGNANT<sup>2,3,1,a</sup> ID, XAVIER PENNEC<sup>2,b</sup> ID, ALAIN TROUVÉ<sup>3,c</sup>,
ANNA CALISSANO<sup>2,4,d</sup> ID

<sup>1</sup>Zuse Institute Berlin, Germany, <sup>a</sup>maignant@zib.de

<sup>2</sup>Université Côte d'Azur and Inria, France, <sup>b</sup>xavier.pennec@inria.fr

<sup>3</sup>Université Paris Saclay and ENS Paris-Saclay, France, <sup>c</sup>alain.trouve@ens-paris-saclay.fr

<sup>4</sup>University College London, UK, <sup>d</sup>a.calissano@ucl.ac.uk

Certain data are naturally modeled by networks or weighted graphs, be
they arterial networks or mobility networks. When there is no canonical la-
beling of the nodes across the dataset, we talk about unlabeled networks. In
this paper, we focus on the question of dimensionality reduction for this type
of data. More specifically, we address the issue of interpreting the feature sub-
space constructed by dimensionality reduction methods. Most existing meth-
ods for network-valued data are derived from principal component analysis
(PCA) and therefore rely on subspaces generated by a set of vectors, which
we identify as a major limitation in terms of interpretability. Instead, we pro-
pose to implement the method called barycentric subspace analysis (BSA),
which relies on subspaces generated by a set of points. In order to provide a
computationally feasible framework for BSA, we introduce a novel embed-
ding for unlabeled networks where we replace their usual representation by
equivalence classes of isomorphic networks with that by equivalence classes
of cospectral networks. We then illustrate BSA on simulated and real-world
datasets, and compare it to tangent PCA.

**1. Introduction.** In recent years, there has been a growing interest in the analysis of network-valued data, namely graphs with scalar weights on the edges, due to the various applications in which they arise. Examples of such data are brain connectivity networks [30, 9, 5], brain arterial networks [13], anatomical trees [37, 10], and mobility networks [36]. Different tasks have been addressed for network-valued data such as classification [33, 20], prediction [2, 27], hypothesis testing [30], data generation [29, 35], and dimensionality reduction [28, 13, 3].

1.1. *Principal component analysis versus barycentric subspace analysis.* Among the different possible tasks, dimensionality reduction techniques are crucial due to the high complexity of network-valued data. A significant proportion of the works addressing dimensionality reduction revolve around the extension of principal component analysis (PCA) from tangent approaches [28, 13, 3]. For interpreting the feature subspace, these methods rely on a set of one-dimensional principal components. It might be however misleading as these components consist in extrinsic Euclidean spans such that they blur the original and naturally discrete nature of data. To address this issue, we propose barycentric subspace analysis (BSA) as an alternative method [25, 26]. The general idea of BSA is the same as PCA: computing a subspace which minimizes the projection error of the dataset. The main difference is that such the feature subspace is formulated as a barycentric subspace, meaning it is generated by a set of points (called reference points) rather than a set of vectors.

MSC2020 subject classifications: Primary 62R30, 62H25, 05C50.

Keywords and phrases: Network-valued data, Dimensionality reduction, Riemannian manifolds.
<figure><img src="image_2.png" /><figcaption>FIG 1. A schematic visualization of barycentric subspace analysis in the Euclidean case.</figcaption></figure>

If the reference points are constrained to be selected among the data points (sample-limited BSA), the interpretation is straightforward. Motivated by these properties, the aim of the paper is to extend BSA to the case of a of set of network-valued data.

1.2. *Unlabeled networks.* When working with a set of network-valued data, there might be some uncertainty around the node labels and the correspondence between nodes across networks. In such cases, the networks are called unlabeled [18]. The ambiguity about the correspondence between nodes can vary from a totally unlabeled setting – two social networks connecting two different groups of individuals – to a partially labeled setting – two molecules involving different types of atoms. In this paper, we focus on defining BSA for unlabeled networks. The case of labeled networks can be addressed directly in the Euclidean setting adapting [25].

To extend BSA to the case of unlabeled networks, we need a suitable geometric framework. In the literature, the uncertainty around the node labels or node correspondence is often modeled via the action of a permutation group on the nodes of the networks [15, 18, 3]. For example, if weighted undirected networks are represented by their symmetric adjacency matrices, their unlabeled counterpart can be represented in the quotient space

$$ (1.1) \qquad \mathcal{G}_n = \mathrm{Sym}(n)/\mathfrak{S}_n. $$

However, performing BSA in this discrete setting is challenging as detailed in [4]. In [28], the authors propose another embedding that is inspired by shape analysis [8] and relies on a continuous group acting on networks. The main drawback of this approach is that the analysis of a network viewed as a shape does not always respect the symmetries that occur naturally in the network.

1.3. *Our proposal.* To overcome both the computational limitations of discrete quotient spaces and the difficulties in interpreting the last model, we propose a relaxed embedding based on the action of the orthogonal group by conjugation onto the set of weighted adjacency matrices

$$ (1.2) \qquad \Gamma_n = \mathrm{Sym}(n)/\mathrm{O}(n). $$

We show that this space is isometric to the space of sorted eigenvalues, automatically providing a full and computationally efficient geometric framework. In particular, we show that the corresponding quotient distance reduces to the well known spectral distance, which is naturally invariant under node relabeling [16, 7]. We call $\Gamma_n$ the spectral graph space of size $n$. Along with the distance function, we equip spectral graph spaces with the tools required to perform BSA, namely an exponential map and its inverse.

1.4. *Structure of the paper.* The paper is organized as follows. In Section 2, we introduce spectral graph spaces and we detail their geometric properties. In Section 3, we formulate barycentric subspace analysis and its sample-limited variant in the context of spectral graph spaces. We run a comparison with principal component analysis in Section 4. To further showcase the applicability of BSA, we study a simulated dataset in Section 5 and a real-world dataset in Section 6.
**2. Spectral graph spaces.** We consider undirected weighted graphs with a finite set of vertices (nodes). Such a graph or network will be given as a triplet $G = (V, E, w)$ consisting of the set $V = \{1, \dots, n\}$ of nodes, a set $E$ of edges i.e. pairs of nodes, and a scalar function $w$ defined on the edges. Treating non-existing edges as unweighted edges, we can encode the network $G$ by its weighted adjacency matrix $X = (x_{ij}) \in \mathbb{R}^{n \times n}$ defined by

$$ (2.1) \qquad x_{ij} = x_{ji} = \begin{cases} w\{i, j\} & \text{if } \{i, j\} \in E \\ 0 & \text{otherwise} \end{cases} $$

such that networks identify with the set $\text{Sym}(n)$ of symmetric matrices of size $n$. Note that we allow for networks to have self-loops, that is for a node to be connected to itself or equivalently for adjacency matrices to have non-zero diagonal entries.

Now in an effort to relax the modeling of unlabeled networks by the action of node permutations, we consider the conjugation action of the orthogonal group $O(n)$ on the space $\text{Sym}(n)$ of symmetric matrices

$$ (2.2) \qquad (R, X) \mapsto RXR^T $$

which transforms one network into another cospectral network, that is a network whose adjacency matrix shares the same spectrum. The orbit of a network consists then exactly of all the other cospectral networks. As the spectrum of a network controls some of its connectivity properties, for example the degree of the nodes [31], the action of $O(n)$ tends to identify networks with a similar structure. Finally, because the permutation group is a subgroup of the orthogonal group, the equivalence class of a given network contains in particular all the other networks obtained by permuting its nodes.

Under such a new group action, unlabeled networks are then given as points of the quotient space $\Gamma_n$ defined by

$$ (2.3) \qquad \Gamma_n = \text{Sym}(n)/O(n). $$

It is a stratified differentiable manifold with singularities at the points where the action of the orthogonal group is not free (see [21] for an introduction to stratified spaces and [19] for another example of such differentiable manifold). We can show that these points correspond exactly to the networks with at least two equal eigenvalues. Let then $\text{Sym}(n)^*$ be the subset of symmetric matrices whose eigenvalues are pairwise distinct. On this subset, the action of the orthogonal group is proper and free such that the principal stratum $\Gamma_n^* = \text{Sym}(n)^*/O(n)$ is a connected differentiable manifold.

**DEFINITION 2.1.** The quotient space $\Gamma_n$ is called the spectral graph space of size $n$.

Let us now investigate the metric structure of the space. Firstly, we study the stratification as a metric space, and secondly, we focus on the Riemannian structure of the principal stratum. Mainly, we show that spectral graph spaces are isometric to convex polyhedral cones. From this result, we derive in particular that the Riemannian logarithm, central to the definition of the notion of barycenter, consists simply in a translation, and moreover that it extends to the whole stratification.

**2.1. Spectral graph spaces as metric spaces.** The graph space $\Gamma_n$ is naturally equipped with the quotient distance induced by the $O(n)$-invariant Frobenius distance on symmetric matrices. More precisely, if $\pi$ denotes the canonical projection, then the function

$$ (2.4) \qquad d(\pi(X), \pi(Y)) = \inf_{R \in O(n)} \|RYR^T - X\|_F $$
is well defined and satisfies all the axioms of a distance function. Such a construction is rather standard in the context of shape analysis [38] and makes $(\Gamma_n, d)$ a complete metric space. Since the permutation group is a subgroup of the orthogonal group as we already highlighted, the distance between two networks thus defined is always smaller than their Frobenius distance up to some permutation of their nodes. In other words, two networks with a similar structure have similar equivalence classes. Now we have the following result.

**PROPOSITION 2.2.** Let $X, Y \in \text{Sym}(n)$. Then the distance between the two corresponding networks is simply the $\ell^2$ spectral distance

$$ (2.5) \qquad d(\pi(X), \pi(Y)) = \left( \sum_{i=1}^{p} |\lambda_i(X) - \lambda_i(Y)|^2 \right)^{1/2} $$

where $\lambda_1(X) \le \dots \le \lambda_n(X)$ and $\lambda_1(Y) \le \dots \le \lambda_n(Y)$ denote the eigenvalues (with their multiplicity) of $X$ and $Y$ respectively.

Let us comment on this result. First, it is interesting to note that $\ell^p$ spectral distances are well known distances between networks [16, 7]. Then, what the proposition above actually tells us is that the map $\lambda = (\lambda_1, \dots, \lambda_n)$ defined accordingly, and constant over each equivalence class, descends to an isometry between the graph space $\Gamma_n$ and the convex cone of sorted vectors of $\mathbb{R}^n$ according to the commutative diagram below.

$$ (2.6) \qquad \begin{tikzcd}[column sep=2.8em, row sep=2.8em] 
    \text{Sym}(n) \arrow[r, "\lambda"] & C_n = \{\mu \in \mathbb{R}^n \mid \mu_1 \le \dots \le \mu_n\} \\
    \pi \arrow[uur, "\cong"'] \arrow[dl, "\pi_n"] & \\
    \Gamma_n & 
    \end{tikzcd} $$

An important consequence of this statement is that the geodesics (distance-minimizing curves) of the graph space $\Gamma_n$ identify with that of the cone, such that there exists a unique minimizing geodesic segment between any two networks described by the straight line segment

$$ (2.7) \qquad \gamma(t) = (1-t)\lambda(X) + t\lambda(Y). $$

For further insight into these different results, we refer to Alekseevsky et al. [1]. In fact, the map $\lambda$ descends to a diffeomorphism between the principal stratum $\Gamma_n^*$, which we recall consists of the quotient by the orthogonal group of symmetric matrices with pairwise distinct eigenvalues, and the interior of the open cone $\text{int}(C_n)$. Let us now extend the diagram 2.6 by equipping the manifold $\Gamma_n^*$ with a Riemannian structure compatible with the existing metric structure.

**2.2. The Riemannian structure of spectral graph spaces.** One way to endow the quotient space $\Gamma_n^*$ with a Riemannian metric is to leverage the Riemannian submersion theorem as detailed in [34]. Namely, that there exists a unique Riemannian metric $g$ on $\Gamma_n^*$ that makes the surjective map $d\pi$ an isometry between the horizontal subbundle

$$ (2.8) \qquad H \text{Sym}(n)^* = (\ker d\pi)^{\perp} $$

and the tangent bundle $T\Gamma_n^*$. Moreover, the corresponding Riemannian distance coincides with the quotient distance as introduced in 2.4. Note that such a construction is also precisely the one developed in Kendall and Le [19] to equip shape spaces with a Riemannian metric. Let us then construct the metric $g$ explicitly. First, we prove the following.
**PROPOSITION 2.3.** Let $x \in \text{Sym}(n)^*$. Let us write $x = Q_X \text{diag}(\lambda_1(X), \dots, \lambda_n(X))Q_X^T$ for $Q_X \in \text{O}(n)$ unique. Then the horizontal subspace at $x$ is

$$ (2.9) \qquad H_X \text{Sym}(n)^* = \{Q_X \text{diag}(h_1, \dots, h_n) Q_X^T \mid (h_1, \dots, h_n) \in \mathbb{R}^n\}. $$

In particular, it means that the horizontal bundle identifies with the vector space $\mathbb{R}^n$.

Now let us consider two tangent vectors at $\pi(X)$ and let $U = Q_X \text{diag}(h_1, \dots, h_n)Q_X^T$ and $V = Q_X \text{diag}(k_1, \dots, k_n)Q_X^T$ be their horizontal lifts at $X$, that is their preimages in the horizontal subspace at $X$. Note that the diagonal entries of $u$ and $v$ are independent of the choice of the representative $X$. By construction, we have

$$ (2.10) \qquad g_{\pi(X)}(d_X \pi(U), d_X \pi(V)) = \langle U, V \rangle_F = \sum_{i=1}^{n} h_i k_i $$

which is non-other than the Euclidean metric of $\mathbb{R}^n$ such that $\Gamma^*$ is isometric to the open cone $\text{int}(C_n)$ of $\mathbb{R}^n$ as a Riemannian manifold, or in other words such that the map $\lambda$ descends to a Riemannian isometry. Another way to prove such a result is thanks to Lee Myers Steerod's theorem that states that an isometry between two Riemannian manifolds is always a Riemannian isometry.

**2.3. Geodesics and extension of the Riemannian logarithm.** Now in particular, this result tells us that on the principal stratum, the Riemannian geodesic joining two networks identify exactly with the geodesic segment described in 2.7 joining their respective spectra, and whose expression can also be derived from that of horizontal geodesics, which describe exactly the geodesics of the quotient space $\Gamma_n^*$ given that $\pi$ is a Riemannian submersion [23]. In this respect, the geodesics of $\Gamma_n$ as a metric space extend the Riemannian geodesics of the main stratum to all the strata. Accordingly, we extend the Riemannian logarithm to $\Gamma_n$ as the inverse of the corresponding unique geodesic segment

$$ (2.11) \qquad \log_{\pi(X)}(\pi(Y)) = \lambda(Y) - \lambda(X). $$

From the notion of Riemannian logarithm, one derives that of a barycenter, or weighted Fréchet mean. With such a tool now in our hands, we then move on to describing the barycentric geometry of the spectral graph spaces and setting up the barycentric subspace analysis.

**3. Barycentric subspace analysis on spectral graph spaces.** Essentially, barycentric subspace analysis (BSA) consists in the approximation of a collection of data points by a lower-dimensional barycentric subspace. In this section, we review the mathematical construction of BSA in the specific case of spectral graph spaces, starting with the definition of a barycentric subspace. We refer to the original work by Pennec [25] for a general formulation of the method on a Riemannian manifold.

**3.1. Barycentric subspaces of spectral graph spaces.** A barycentric subspace is defined for a finite family of points as the set of all their barycenters. Concretely, the barycentric subspace of the points $\pi(A_0), \dots, \pi(A_k) \in \Gamma_n$ is defined as

$$ (3.1) \qquad \text{BS}(\pi(A_0), \dots, \pi(A_k)) = \bigcup_{\substack{w_0, \dots, w_k \in \mathbb{R} \\ w_0 + \dots + w_k = 1}} \left\{ \pi(X) \in \Gamma_n \mid \sum_{i=0}^{k} w_i \log_{\pi(X)}(\pi(A_i)) = 0 \right\}. $$

The points $\pi(A_0), \dots, \pi(A_k)$ are then called the reference points of the subspace. From the definition of the logarithm in 2.11, we derive an explicit formulation of the barycentric subspaces of spectral graph spaces.
**THEOREM 3.1.** Let $A_0, \dots, A_k \in \text{Sym}(n)$. Then the barycenters of $\pi(A_0), \dots, \pi(A_k)$ identify exactly with barycenters of $\lambda(A_0), \dots, \lambda(A_k)$ and the corresponding barycentric subspace is isometric to the convex set

$$ (3.2) \quad \text{BS}(\pi(A_0), \dots, \pi(A_k)) \approx \left\{ \sum_{i=0}^{k} w_i \lambda(A_i) \;\middle|\; \sum_{i=0}^{k} w_i = 1 \right. \\ \qquad \left. \text{and } \sum_{i=0}^{k} w_i \lambda_r(A_i) \le \sum_{i=0}^{k} w_i \lambda_{r+1}(A_i) \text{ for } 1 \le r \le n-1 \right\}. $$

It is a convex polytope with at most $n-1$ facets, and it corresponds to the cross-section of the cone $C_n$ passing through the vectors $\lambda(A_0), \dots, \lambda(A_k)$. In particular, the barycentric subspace of $\pi(A_0), \dots, \pi(A_k)$ always contains the convex hull of $\lambda(A_0), \dots, \lambda(A_k)$.

Note that if the spectra $\lambda(A_0), \dots, \lambda(A_k)$ are affinely independent, then the barycentric subspace of $\pi(A_0), \dots, \pi(A_k)$ is a subspace of $\Gamma_n$ of dimension $k$. Now, one can also write the polytope $\text{BS}(\pi(A_0), \dots, \pi(A_k))$ as the intersection of the following $n-1$ half-spaces

$$ (3.3) \qquad \alpha_{11} \lambda_1(X) + \dots + \alpha_{1n} \lambda_n(X) \ge \beta_1 \\ \qquad \qquad \vdots \\ \qquad \alpha_{n-11} \lambda_1(X) + \dots + \alpha_{n-1n} \lambda_n(X) \ge \beta_{n-1} $$

where for $1 \le r \le n-1$, the pair $(\alpha_r, \beta_r)$ solves

$$ (3.4) \qquad \begin{bmatrix} \lambda_1(A_0) & \cdots & \lambda_n(A_0) & 1 \\ \vdots & \ddots & \vdots & \vdots \\ \lambda_1(A_k) & \cdots & \lambda_n(A_k) & 1 \end{bmatrix} \begin{bmatrix} \alpha_{r1} \\ \vdots \\ \alpha_{rn} \\ -\beta_r \end{bmatrix} = \begin{bmatrix} \lambda_{r+1}(A_0) - \lambda_r(A_0) \\ \vdots \\ \lambda_{r+1}(A_k) - \lambda_r(A_k) \end{bmatrix}. $$

The derivations leading to 3.3 and 3.4 as well the proof of the theorem are detailed in Appendix A. We illustrate these results for $k=2$ in Appendix B.

BSA then consists in finding the barycentric subspace that minimizes the distance to the dataset, in our case the spectral distance. Optimizing for such a subspace is a combination of the search for the reference points and the projection of the data onto the corresponding barycentric subspace. Since the latter is a projection onto a closed convex set, it has a unique minimum which can be achieved by usual gradient descent methods on $\mathbb{R}^n$. Without additional constraint, optimal reference points are not well defined, however, as the same $k$-dimensional barycentric subspace can be generated by any family of $k+1$ affinely independent points of the subspace. Here, we consider more specifically what is referred to as sampled-limited barycentric Subspace Analysis (sample-limited BSA) in the original paper by Pennec [25].

**3.2. Sample-limited barycentric subspace analysis.** In this BSA variant, the optimal reference points are selected from among the data points so they are directly interpretable as such, following the same approach as in $k$-medoids clustering [17]. In particular, rather than reference spectra, we obtain reference networks. This in turn increases the interpretability of the low-dimensional subspace, whether used for visualization purposes or as a support for lower-dimensional analyses of the dataset. Mathematically speaking, sample-limited BSA
consists of the following optimization problem

$$
\begin{equation} \tag{3.5}
\begin{array}{ll}
\underset{\substack{1 \le i_0 < \dots < i_k \le N \\ w_{ij} \in \mathbb{R}}} \text{minimise} & \displaystyle \sum_{i=1}^{N} \sum_{r=1}^{n} \left| \lambda_r(X_i) - \sum_{j=0}^{k} w_{ij} \lambda_r(X_{i_j}) \right|^2 \\[1.5ex]
\forall 1 \le i \le N, w_{i_0} + \dots + w_{i_k} = 1 & \text{subject to} \quad \sum_{j=0}^{k} w_{ij} \lambda_r(X_{i_j}) < \sum_{j=0}^{k} w_{ij} \lambda_{r+1}(X_{i_j}) \quad \forall i, r.
\end{array}
\end{equation}
$$

The search for reference points is now a well defined but combinatorial optimization problem. This poses a major computational challenge. One way of reducing the computational cost of the problem would be to restrict the search of the reference networks to pre-assigned clusters. This approach seems particularly relevant in the case of a mixed distribution as for example for the dataset studied in Section 5. Another way of doing so could be to enforce some minimal distance between any two reference networks, thereby selecting more contrasting structures. Additionally, note that the datasets we are interested in the context of geometric data analysis are in practice rather small, since larger ones are already well-covered by linear techniques.

Finally, to further increase the interpretability of the projection, we may enforce the ref-
erence networks to be extreme points of the dataset. More precisely, we restrict the projec-
tion to the convex hull of the reference networks and solve Problem 3.5 with the additional
constraint that the barycentric weights must be positive. We refer to this method as sample-
limited convex BSA. It is particularly suited for the identification of archetypal points [6].
Another aspect is that sample-limited convex BSA should be more robust to noise than stan-
dard sample-limited BSA as the volume of the convex hull of given reference networks is
significantly smaller than that of their barycentric subspace (see Appendix B).

3.3. *Spectral network reconstruction.* In terms of interpretability, the use of spectral graph spaces to implement BSA raises its own issues. Precisely, the projection of a network onto the optimal barycentric subspace is no longer a network but only the spectrum of one. Still, it might be interesting for some further analysis to recover such a projection as a network, for instance when combining dimension reduction with another method specifically designed for network-valued data. Given the projected spectrum of a data point, a canonical candidate is the network with this spectrum that is the closest to the data point (with respect to the Frobenius distance on the space of symmetric matrices). Based on a similar proof to that of Proposition 2.2, we can see that such an optimum corresponds exactly to the conjugation of the projected spectrum by the orthogonal matrix that achieves the eigendecomposition of the original data point

$$
(3.6) \quad \hat{x}_i = R_i \operatorname{diag} \left( \sum_{j=1}^{k} w_{ij} \lambda(X_{i_j}) \right) R_i^T \quad \text{where} \quad X_i = R_i \operatorname{diag}(\lambda(X_i)) R_i^T.
$$

The reconstruction error is then the same as the projection error in the spectral graph spaces. Note that the orthogonal matrix $R_i$ is not unique when $X_i$ has multiple equal eigenvalues. Whenever this occurs, we simply choose one of the possible eigendecompositions. This reconstruction algorithm can be seen as a low-rank approximation of $X_i$. An artifact introduced by the procedure is the occurrence of self-loops within the reconstructed networks, despite the data not having any in the first place (see Figure 2). This problem can be avoided by reconstructing explicitly a network with no diagonal entries instead, for example the closest such network to the original data point. This is possible because every zero trace matrix is
<figure><img src="image_3.png" /><figcaption>FIG 2. Spectral network reconstruction of the projection of a random network onto the barycentric subspace of a complete network and a star. The networks here are drawn in such a way that the length of the edge {$i, j$} satisfies $x_{ij} \propto \exp[-\ell\{i, j\}^2]$. To control the edge length, we rely on the Kamada-Kawai layout available in the Python library NetworkX [14]. We keep this visualization throughout the paper.</figcaption></figure>

orthogonally similar to at least one zero diagonal matrix [11]. However, we observe in the context of BSA that self-loops reconstructed by the first algorithm carry a marginal weight, so that the two algorithms output close solutions. It should be the case as long as the projection error is small enough.

**4. A comparison with tangent principal component analysis.** In this first experiment, we wish to underline the better interpretability of the low-dimensional subspace computed by sample-limited BSA over that computed by existing dimensionality reduction methods such as geodesic PCA [3] or tangent PCA [13]. We focus the comparison on tangent PCA for it is straightforward to implement. Through our experiment, we show in particular that when it comes to analyzing the structural variability of networks, reference points, that is networks, are more intuitive quantities than vector components.

To this end, we consider a two-dimensional set of networks featuring changes of topology, in the sense of the edge sets $E$. We generate the dataset as a two-parameter dataset $X(s_i, t_i)$ in the following fashion. The network $X(s, t)$ consists of $n = 6$ nodes, two of which are always connected by a central edge of weight $s$, while the connectivity of the others depends on the value of $t$. If $t = 0$, then $X(s, t)$ consists of only five edges and has a tree topology. Else, two edges of weight $|t|$ are added so that $X(s, t)$ has either an "eight" topology if $t < 0$ or a "hourglass" topology if $t > 0$ as illustrated in Figure 3. Overall, the dataset consists in $N = 16$ networks sampled from the uniform distribution over $[\frac{1}{2}, \frac{3}{2}] \times [-\frac{1}{2}, \frac{1}{2}]$. We apply sample-limited BSA and its convex variant for $k = 3$ on one hand, and tangent PCA with 2 components on the other hand. Let us recall that tangent PCA is defined for a dataset on a Riemannian manifold as the PCA of the corresponding Riemannian logarithms in the tangent space at the Fréchet mean [12]. In the paper by Guo et al. [13], the method is redefined on the graph space introduced in Section 1 – where the permutation matrices act instead of the orthogonal matrices. Although such a quotient space is not a manifold, one can still define a
<figure><img src="image_4.png" /><figcaption>FIG 3. A two-parameters dataset. The first parameter *s* can be read directly as the weight of the central edge. Now if the networks represent molecules, the parameter *t* can be interpreted as the angle between the central bond and the extreme ones. It generates three different topologies: an "eight" topology, a tree topology, and a "hourglass" topology (from left to right).</figcaption></figure>

transformation analog to the logarithms by optimally aligning the networks of the dataset to their mean and then subtracting such. Note that the alignments may be computed exactly in this experiment because the networks have a small enough numbers of nodes. The two first components of tangent PCA define then a 2-dimensional subspace. In order to interpret such a subspace, one relies in practice on the geodesics spanned by the components, in our case two straight lines of the space of symmetric matrices. As for sample-limited (convex) BSA, we base our interpretation on the three reference networks selected.

A first observation is that the projections output by the three methods provide planar repre-
sentations of the dataset which are comparable both visually and numerically. Precisely, the
squared projection error is 4, 7 × 10⁻² for tangent PCA while it is 3, 4 × 10⁻² for sample-
limited BSA and 5, 2 × 10⁻² for sample-limited convex BSA. Now, although the two first
components of tangent PCA cover most variability i.e. 88% of the dataset, they display in
Figure 4 rather misleading deformations of the mean network for they generate topologies
that are not observed in the original dataset. This can be explained partly by the mean it-
self having an unobserved topology, but most importantly because the topological changes
characterizing the dataset occur along non linear deformations. More generally, due to the
discrete nature of topology, structural variations within of a population of networks are dif-
ficult to reproduce from continuous objects such as components. Instead, the three reference
<figure><img src="image_5.png" /><figcaption>FIG 4. Tangent PCA of the two-parameters dataset. The straight lines spanned by the two first components (in red and blue) can be read as deformations of the Fréchet mean (framed). As expected, neither the mean network nor the deformed networks belong to the two-parameters family generating the dataset. Instead, the two deformations interpolate linearly the "eight" topology and the "hourglass" one, eventually resulting in negative edges. As a result, they provide a misleading interpretation of the feature subspace, and by extension of the dataset itself, if one would rely on the observation that the projection explains almost 90% of the data variance.</figcaption></figure>

networks picked by sample-limited BSA and its convex variant (see Figures 5 and 6) recover by definition existing topologies. Additionally, we observe that the convex constraint allows to retrieve reference networks with more different structures for they are more distant to each other, though without coping with a greater projection error.

**5. Reference points of a clustered dataset.** In this second experiment, we investigate the case of a clustered dataset. In particular, we are interested in how the reference networks selected by sample-limited BSA localize with respect to the clusters.

We consider a data set of size *N* = 15 generated from three networks with a clearly identified topology: the star, the 2-fully meshed star [22], and the complete network, all with same number of nodes *n* = 10. The generative model simply consists in adding Gaussian noise of standard deviation *σ* to the weight of the edge between any two different nodes, eventually adjusting the resulting weight to its absolute value so that it is positive. We fix *σ* = 5 × 10⁻² such that the clusters those formed are well separated in the spectral domain. Because the spectral distance is less than the Frobenius distance for any node permutation, the clusters
<figure><img src="image_6.png" /><figcaption>FIG 5. Sample-limited BSA of the two-parameters dataset. The barycentric coordinates thus assigned to the data points define their projection onto the plane only up to a choice of coordinates for the reference networks. We fix the position of those in such a way that their spectral distance to each other is preserved. The projected distribution is similar to that achieved with PCA, except that it appears more stretched. This is likely a consequence of the bias of PCA towards Gaussian distributions. The reference networks computed alongside the coordinates reflect the two prevailing topologies within the dataset, yet two of the three reference networks are fairly similar and therefore redundant in terms of interpretability. This issue is overcome by adding the convex constraint.</figcaption></figure>

are then also well separated in the adjacency domain as we observe in Figure 7. We perform sample-limited backward BSA, a reformulation of sample-limited BSA as the search for the nested subspaces $BS(X_{i_0}, \dots, X_{i_k}) \subset BS(X_{i_0}, \dots, X_{i_{k-1}}) \subset \dots \subset BS(X_{i_0})$ such that the mean squared projection error increases as less as possible each time one reference point is removed [25]. Note that the first subset is the optimal subspace found by sample-limited BSA. We set $k = N - 1$ such that the mean squared projection error is initially 0, and we analyze how the error increases as we remove reference points.

The main observation in Figure 7 is that the three last reference networks are precisely located in one cluster each, regardless of whether we apply the convex or the non-convex variant of sample-limited backward BSA. In fact, removing one of those three points results in a sudden increase of the mean squared projection error, suggesting that the two-dimensional barycentric subspace is the smallest suitable for projecting the dataset. Additionally, though sample-limited backward BSA should theoretically achieve a lower projection error than its convex variant, there is no striking discrepancy for this dataset. Still, we note that the reference networks selected by the convex variant have a more accentuated structure. While this can be attributed to them being constrained to be distant from each other as we already pointed out in the previous experiment, another possible explanation comes from the geometry of spectral graph spaces. Precisely, the three template networks are singular networks and lie on the boundary of the cone $C_{10}$. Consequently, they also lie outside of the open convex
<figure><img src="image_7.png" /><figcaption>FIG 6. Sample-limited convex BSA of the two-parameters dataset. Under the convex constraint, the three reference networks selected are extreme enough for us to witness three very different structure – and in fact, almost all three possible topologies. They therefore provide a strong intuition of the data structural variability. The other side of the coin is that the variability of the projected distribution is less than with standard sample-limited BSA. In particular, several points lie on the boundary of the triangular projection domain delimited by the three reference networks. An even more critical manifestation of this phenomenon is, for example, that the projection of the network X₀ is superimposed on that of the reference network X₁.</figcaption></figure>

hull of the dataset such that the data points with the closer structure should be among the extreme points of the dataset.

**6. A small case study on European airlines route maps.** As a concluding experiment, we run a sample-limited backward convex BSA of a real dataset derived from the OpenFlight database [24], a database collecting airline traffic data. The dataset consists of 12 networks representing the route maps of the main European airlines as follows: each node stands for a macro-region in and around Europe (Northern Europe, Eastern Europe, Southern Europe, Western Europe, North Africa, Western Asia) and the edge between two nodes is weighted by the proportion of routes operated between the two corresponding macro-regions out of all routes operated by the airline worldwide. Although the nodes are canonically labeled by the macro-regions, we choose to place ourselves in the unlabeled framework so as to compare the strategy of the airlines independently of their geographical establishment. For example, while Air France and Aeroflot base their traffic in two different regions, Western Europe on the one hand and Eastern Europe on the other, both airlines nonetheless share the same strategy, namely a centralized strategy (see Figure 9). Following a similar heuristic as in the previous experiment, we rely on the backward implementation of sample-limited convex BSA to select the dimension for projecting the dataset, and we interpret the lower-dimensional representation thus obtained in the light of the original data.
<figure><img src="image_8.png" /><figcaption>FIG 7. A data set with three clusters. The networks are clustered according to their underlying topology: star, 2-fully meshed star or complete network (of edge weight 1/2).</figcaption></figure>

The evolution of the mean squared projection error along the sequence of nested barycentric subspaces recovered suggests that the one-dimensional subspace offers the best trade-off between accuracy and dimension reduction for projecting the dataset, although the threshold is not as clear-cut as for a well-clustered dataset (see Figure 10). The two reference airlines that define the one-dimensional subspace have two opposite topologies: a well-connected one (three-meshed star) and a sparse one (star). Schematically, these two topologies exemplify the strategy of low-cost airlines and national airlines respectively. The projection of the dataset onto the one-dimensional barycentric subspace can then be interpreted as the level of centralization of the airlines.

**7. Conclusions and Further Developments.** Motivated by the need to define interpretable techniques for unlabeled network-valued data, we adapt to this particular case the method called barycentric subspace analysis (BSA) and a priori defined for manifold-valued data. Typically, unlabeled networks are studied as elements of a quotient space resulting from the action of node permutations. However, this kind of group action is not compatible with BSA in terms of computational cost. To overcome such an issue, we introduce spectral graph spaces, a novel embedding obtained by relaxing the action – by conjugation – of the permutation group onto the set of weighted adjacency matrices into that of the orthogonal group. While the former identifies isomorphic networks, the latter identifies cospectral networks.
<figure><img src="image_9.png" /><figcaption>FIG 8. Sample-limited backward BSA (left) and sample-limited backward convex BSA (right) of the clustered dataset. For both methods, we plot the evolution of the mean squared projection error with the dimension of the barycentric subspace (top). Note that the y-axis is cut below the value of projection error on the zero-dimensional subspace in order to be able to observe more clearly the first increments. The key information to take from the plots is that the projection error suddenly increases – in ratio – when jumping from dimension 2 to dimension 1. Coincidentally, The two-dimensional barycentric subspace is defined by one reference network per cluster (bottom). Coming back to the plots, it is interesting to check that in the standard variant, the error does not increase before the dimension of the barycentric subspace becomes lower than that of the data space, in this case 9. It is one less than the number of eigenvalues due to all the networks having no self-loops and therefore a trace zero, a property preserved by the action of the orthogonal group.</figcaption></figure>

Because spectral graph spaces are isometric to convex linear cones, namely that defined by sorted eigenvalues, they allow for efficient computations, solving the aforementioned issue. We implement BSA on such spectral graph spaces. BSA consists in computing the barycentric subspace minimizing the projection error of the dataset, where barycentric subspaces are

<figure><img src="image_10.png" /><figcaption>FIG 9. Route maps of 12 European airlines linking macro-regions in and around Europe. The topology of the networks reflects the strategy of the airlines and in particular the difference in this respect between low-cost airlines (three-meshed star) and most national airlines (star).</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>FIG 10. Sample-limited backward convex BSA of the European airline dataset. We plot the evolution of the mean squared projection error (MSE) with the dimension of the barycentric subspace (left). The most sudden increase of the projection error occurs here when jumping from dimension 1 to dimension 0, suggesting dimension 1 as the lowest suitable dimension for projection the dataset. The two corresponding reference networks (right) illustrate strategies that are opposite: easyJet covers the largest number of routes, while Swiss is centralized with a focus on one particular route. Then, the projection of the dataset onto the one-dimensional subspace (right) ranks the airlines from Easy Jet to Swiss, and can therefore be interpreted as a measure of how equitably distributed their traffic is. Towards the center of the projected distribution, we have on one hand AirFrance and Aeroflot, implementing a centralized but more equally weighted strategy, and on the other hand Ryanair and Alitalia, for they are not strictly centralized but cover fewer routes than easyJet. In fact, Alitalia's topology looks similar to both that of EasyJet and Swiss, namely the three-meshed star topology and the star one.</figcaption></figure>

subspaces generated by a set of reference points. Thanks to the encoding of the feature sub-
space by reference points, that is networks, its interpretation is straightforward. To illustrate
our point, we compare BSA with tangent PCA and apply it both to simulated data sets and to
a real-world data set consisting of flight connectivity networks. In summary, spectral graph
space offer an efficient embedding for analyzing unlabeled networks. They represent a solid
trade-off between computational efficiency and interpretability of the model. Furthermore,
BSA is also easier to interpret than component-based dimensionality reduction methods such
as PCA.

Both the geometric and the methodological aspects offer further possible developments. With regard to geometric aspects, since the spectrum of the Laplacian matrix of a network is more widely studied in the literature than that of its adjacency matrix, it would be interesting to investigate an extension of spectral graph spaces to the set of Laplacian matrices. Such an extension is however not straightforward as the space of Laplacian matrices is not closed under the action of the orthogonal group. Regarding the methodological aspects, we could further explore the interplay between the optimal reference networks defined by convex BSA and the singular points of spectral graph spaces, which often correspond to networks with a very clear structure. Indeed, on one hand, the convex constraint enforces reference networks to be extreme points of the sampled distribution. On the other hand, singular networks lie at the boundary of the spectral graph spaces. Finally, we believe our framework could provide interesting insights about real-world network-valued data such as structural brain connectiv- ity matrices and electroencephalography (EEG) data [5, 32, 30].

**Acknowledgment.** We would like to thank Peter Michor for the fruitful discussion on the geometric aspects discussed in this paper. The work was supported by the ERC grant \#786854 G-Statistics from the European Research Council, by the Chapman Fellowship scheme by the Department of Mathematics at Imperial College and by the DFG funded Cluster of Excellence EXC 2046 MATH+ (project ID: AA1-20).
# APPENDIX A: PROOFS

## PROOF OF PROPOSITION 2.2. Consider the functional

$$ (A.1) \qquad d_{X,Y} : R \mapsto \|RYR^T - X\|^2. $$

Its gradient at $R \in O(n)$ with respect to the Frobenius metric on $T_R O(n) = \text{Skew}(n) \cdot R$ is

$$ (A.2) \qquad \text{grad } d_{X,Y}(R) = 2RYR^T X - 2XRYR^T \in \text{Skew}(n). $$

Now the functional $d_{X,Y}$ is minimized only if the first order condition

$$ (A.3) \qquad \text{grad } d_{X,Y}(R) \in \text{Skew}(n)^\perp = \text{Sym}(n) $$

applies. This condition is satisfied if and only if $X$ and $RYR^T$ commute with each other. Let then $R \in O(n)$ such that this is the case. Then $X$ and $RYR^T$ can be diagonalized in a common orthonormal basis

$$ (A.4) \qquad X = Q_{X,Y} D_X Q_{X,Y} \quad \text{and} \quad RYR^T = Q_{X,Y} D_Y Q_{X,Y} $$

and we have

$$ (A.5) \qquad d_{X,Y}(R) = \|D_X - D_Y\|^2. $$

Up to a permutation $\sigma \in \mathfrak{S}(n)$, this amounts exactly to

$$ (A.6) \qquad d_{X,Y}(R) = \sum_{i} |\lambda_{i}(X) - \lambda_{\sigma(i)}(Y)|^2. $$

Finally, the rearrangement inequality states that

$$ (A.7) \qquad \sum_{i=1}^{n} \lambda_{i} \mu_{\sigma(i)} \leq \sum_{i} \lambda_{i} \mu_{i} $$

for any permutation $\sigma$. Therefore, the functional $d_{X,Y}$ is minimal only if $\sigma = id$ and its minimum is the quantity

$$ (A.8) \qquad d_{X,Y}(R) = \sum_{i} |\lambda_{i}(X) - \lambda_{i}(Y)|^2. $$

Moreover, the reader may want to check that such a choice of $\sigma$ corresponds exactly to taking $R = Q_X Q_Y^T$ where $Q_X$ and $Q_Y$ are such that $X = Q_X \text{diag}(\lambda_1(X), \dots, \lambda_n(X)) Q_X^T$ and $Y = Q_Y \text{diag}(\lambda_1(Y), \dots, \lambda_n(Y)) Q_Y^T$. Conversely, if $R = Q_X Q_Y^T$ then Equation A.8 is satisfied. $\square$

## PROOF OF PROPOSITION 2.3. Let $X \in \text{Sym}(n)^*$. The horizontal subspace at $X$ is defined by the following construction

$$ (A.9) \qquad H_X \text{Sym}(n)^* + (\text{ker } d_X \pi)^\perp = T_X \text{Sym}(n)^*. $$

The subspace $\text{ker } d_X \pi$ consists exactly of the tangent vectors at $X$ which are also tangent to the fiber $F$ of $\pi$ at $\pi(X)$. Consider then a curve $\gamma$ in the fiber $F$ such that $\gamma(0) = X$. It can be written as

$$ (A.10) \qquad \gamma(s) = R(s)XR(s)^T $$

where $R$ is a curve in $O(n)$ with $R(0) = \text{Id}$. Then let us compute the derivative of $\gamma$ at 0

$$ (A.11) \qquad \dot{\gamma}(0) = \dot{R}(0)XR(0)^T + R(0)X\dot{R}(0)^T. $$
Since $O(n)$ is a Lie group of Lie algebra $\text{Skew}(n)$, then there exists $A \in \text{Skew}(n)$ such that $\dot{R}(0) = A$ and we have

$$ (A.12) \qquad \dot{\gamma}(0) = AX + XA^T. $$

Now let $U \in \text{Hor}_X \text{Sym}(n)^*$. Then for all $A \in \text{Skew}(n)$ we have

$$ (A.13) \qquad \langle U, AX - XA \rangle = 0, $$

that is

$$ (A.14) \qquad \langle XU - UX, A \rangle = 0. $$

Since $XU - UX$ is skew-symmetric, this is equivalent to

$$ (A.15) \qquad XU - UX = 0. $$

Finally, if $Q_X \in O(n)$ is such that $x = Q_X \text{diag}(\lambda_1(X), \dots, \lambda_n(X))Q_X^T$, then $U$ commutes with $X$ if and only if $Q_X^T U Q_X$ is a diagonal matrix. $\square$

**PROOF OF THEOREM 3.1.** Let $X \in \text{Sym}(n)$. Then $\pi(X)$ belongs to the barycentric subspace of $\pi(A_0), \dots, (A_k)$ if and only if there exists $w_0, \dots, w_n \in \mathbb{R}$ summing to 1 such that

$$ (A.16) \qquad \sum_{i=0}^{k} w_i \log_{\pi(X)} (\pi(A_i)) = 0. $$

Following the definition of the logarithm, this condition translates directly into

$$ (A.17) \qquad \sum_{i=0}^{k} w_i (\lambda(A_i) - \lambda(X)) = 0. $$

Adding the constraint that $w_0, \dots, w_k$ have to sum to 1, we get

$$ (A.18) \qquad \lambda(X) = \sum_{i=0}^{k} w_i \lambda(A_i). $$

Now we recall that by definition of the map $\lambda$, we must have $\lambda(X) \in C_n$. Therefore, $w_0, \dots, w_k$ should satisfy

$$ (A.19) \qquad \forall 1 \le r < n, \quad \sum_{i=0}^{k} w_i \lambda_r(A_i) \le \sum_{i=0}^{k} w_i \lambda_{r+1}(A_i). $$

Let us now derive from Equations A.18 and A.19 a simple geometric constraint on $\lambda(X)$. Let us denote

$$ (A.20) \qquad \Lambda = \begin{bmatrix} \lambda_1(A_0) \cdots \lambda_1(A_k) \\ \vdots & \vdots \\ \lambda_n(A_0) \cdots \lambda_n(A_k) \\ 1 & \cdots & 1 \end{bmatrix} \quad \text{and} \quad \theta_r = \begin{bmatrix} \lambda_{r+1}(A_0) - \lambda_r(A_0) \\ \vdots \\ \lambda_{r+1}(A_k) - \lambda_r(A_k) \end{bmatrix} $$

On one hand, Equation A.19 may be rearranged in

$$ (A.21) \qquad \forall 1 \le q < p, \quad [w_0 \cdots w_k]^T \theta_r \ge 0. $$

On the other hand, Equation A.18 together with the constraint $\sum_{i=0}^{k} w_i = 1$ is expressed as

$$ (A.22) \qquad [\lambda_1(X) \cdots \lambda_n(X) 1] = \Lambda [w_0 \cdots w_k]. $$

Notice that $\theta_r \in \text{ker}(\Lambda)^\perp$. It follows that $\theta_r \in \text{im}(\Lambda^T)$. For all $1 \le r < n$, let then $\alpha_{r1}, \dots, \alpha_{rn}$ and $\beta_r$ solve

$$ (A.23) \qquad \Lambda^T [\alpha_{r1} \cdots \alpha_{rn} - \beta_r] = \theta_r. $$

Then Equation A.19 is equivalent to $\lambda(X)$ satisfying $\alpha_r^T \lambda(X) \ge \beta_r$. $\square$
APPENDIX B: ILLUSTRATING THEOREM 3.1 FOR $k = 2$

The barycentric subspace of three reference networks $\pi(A_0), \pi(A_1)$ and $\pi(A_2) \in \Gamma_n$ is a two-dimensional polytope, that is a polygon. It can be drawn isometrically in the plane in the following way. Let $a_0, a_1$ and $a_2 \in \mathbb{R}^2$ be such that $\|a_i - a_j\|_{\mathbb{R}^2} = \|\lambda(A_i) - \lambda(A_j)\|_{\mathbb{R}^n}$. Then for $1 \le r \le n-1$ the equation in the plane of the r-th facet of the polygon is

$$ (B.1) \qquad \alpha_{r1}x + \alpha_{r2}y \ge \beta_r $$

where the three coefficients $\alpha_{r1}, \alpha_{r2}, \beta_r$ solve

$$ (B.2) \qquad \begin{bmatrix} a_{01} & a_{02} & 1 \\ a_{11} & a_{12} & 1 \\ a_{21} & a_{22} & 1 \end{bmatrix} \begin{bmatrix} \alpha_{r1} \\ \alpha_{r2} \\ -\beta_r \end{bmatrix} = \begin{bmatrix} \lambda_{r+1}(A_0) - \lambda_r(A_0) \\ \lambda_{r+1}(A_1) - \lambda_r(A_1) \\ \lambda_{r+1}(A_2) - \lambda_r(A_2) \end{bmatrix}. $$

In Figure 11, we visualize such a polygon for different choices of $n$ and reference networks.

<figure><img src="image_12.png" /><figcaption>FIG 11. Barycentric subspace of 3 networks. (Top) The barycentric subspace of the 3 networks is isometric to a closed tetragon. It corresponds to the situation where the subspace has the maximum number of sides, that is one less than the number of nodes of the reference networks. (Center) The second example illustrates a situation where the barycentric subspace is isometric to a polygon that is not closed. (Bottom) In the last example, where we doubled the number of nodes, the intersection of the half-planes is not minimal and the corresponding polygon is an hexagon.</figcaption></figure>
REFERENCES

[1] ALEKSEEVSKY, D., KRIEGL, A., LOSIK, M. and MICHOR, P. W. (2003). The Riemannian geometry of orbit spaces – the metric, geodesics, and integrable systems. *Publicationes Mathematicae Debrecen* **62** 247–276. https://doi.org/10.5486/PMD.2003.2821

[2] CALISSANO, A., FERAGEN, A. and VANTINI, S. (2022). Graph-valued regression: Prediction of unlabelled networks in a non-Euclidean graph space. *Journal of Multivariate Analysis* **190** 104950. https://doi.org/10.1016/j.jmva.2022.104950

[3] CALISSANO, A., FERAGEN, A. and VANTINI, S. (2024). Populations of unlabelled networks: graph space geometry and generalized geodesic principal components. *Biometrika* **111** 147–170. https://doi.org/10.1093/biomet/asad024

[4] CALISSANO, A., MAIGNANT, E. and PENNEC, X. (2023). Towards Quotient Barycentric Subspaces. In *Geometric Science of Information* (F. NIELSEN and F. BARBARESCO, eds.) 366–374. Springer Nature Switzerland. https://doi.org/10.1007/978-3-031-38271-0_36

[5] CALISSANO, A., PAPADOPOULO, T., PENNEC, X. and DESLAURIERS-GAUTHIER, S. (2024). Graph alignment exploiting the spatial organization improves the similarity of brain networks. Human brain mapping. *Human brain mapping* **45** e26554. https://doi.org/10.1002/hbm.26554

[6] CUTLER, A. and BREIMAN, L. (1994). Archetypal Analysis. *Technometrics* **36** 338–347. https://doi.org/10.2307/1269949

[7] DONNAT, C. and HOLMES, S. (2018). Tracking network dynamics: A survey using graph distances. *The Annals of Applied Statistics* **12** 971–1012. https://doi.org/10.1214/18-AOAS1176

[8] DRYDEN, I. L. and MARDIA, K. V. (2016). *Statistical Shape Analysis, with Applications in R*. *Wiley Series in Probability and Statistics*. John Wiley & Sons, Ltd. https://doi.org/10.1002/9781119072492

[9] DURANTE, D., DUNSON, D. B. and VOGELSTEIN, J. T. (2017). Nonparametric Bayes Modeling of Populations of Networks. *Journal of the American Statistical Association* **112** 1516–1530. https://doi.org/10.1080/01621459.2016.1219260

[10] FERAGEN, A., OWEN, M., PETERSEN, J., WILLE, M. M. W., THOMSEN, L. H., DIRKSEN, A. and DE BRUIJNE, M. Tree-Space Statistics and Approximations for Large-Scale Analysis of Anatomical Trees. In *Information Processing in Medical Imaging* (J. C. GEE, S. JOSHI, K. M. POHL, W. M. WELLS and L. ZÖLLEI, eds.) 74–85. Springer Berlin, Heidelberg. https://doi.org/10.1007/978-3-642-38868-2_7

[11] FILLMORE, P. A. (1969). On Similarity and the Diagonal of a Matrix. *The American Mathematical Monthly* **76** 167–169. https://doi.org/10.1080/00029890.1969.12000162

[12] FLETCHER, P. T., LU, C., PIZER, S. M. and JOSHI, S. (2004). Principal geodesic analysis for the study of nonlinear statistics of shape. *IEEE Transactions on Medical Imaging* **23** 995-1005. https://doi.org/10.1109/TMI.2004.831793

[13] GUO, X., SRIVASTAVA, A. and SARKAR, S. (2021). A Quotient Space Formulation for Generative Statistical Analysis of Graphical Data. *Journal of Mathematical Imaging and Vision* **63** 735–752. https://doi.org/10.1007/s10851-021-01027-1

[14] HAGBERG, A. A., SCHULT, D. A. and SWART, P. J. (2008). Exploring Network Structure, Dynamics, and Function using NetworkX. In *Proceedings of the 7th Python in Science Conference* (G. VAROQUAUX, T. VAUGHT and J. MILLMAN, eds.) 11 - 15.

[15] JAIN, B. J. and OBERMAYER, K. Structure Spaces. *Journal of Machine Learning Research* **10** 2667–2714. https://doi.org/10.5555/1577069.1755876

[16] JURMAN, G., VISINTAINER, R. and FURLANELLO, C. (2011). An introduction to spectral distances in networks. In *Proceedings of the 2011 Conference on Neural Nets WIRN10: Proceedings of the 20th Italian Workshop on Neural Nets* 227–234. IOS Press. https://doi.org/10.5555/1940632.1940662

[17] KAUFMAN, L. and ROUSSEEUW, P. J. (1990). *Partitioning Around Medoids (Program PAM)* In *Finding Groups in Data* 2, 68-125. John Wiley and Sons, Ltd. https://doi.org/10.1002/9780470316801.ch2

[18] KOLACZYK, E. D., LIN, L., ROSENBERG, S., WALTERS, J. and XU, J. (2020). Averages of unlabeled networks: Geometric characterization and asymptotic behavior. *The Annals of Statistics* **48** 514–538. https://doi.org/10.1214/19-AOS1820

[19] LE, H. and KENDALL, D. G. (1993). The Riemannian Structure of Euclidean Shape Spaces: A Novel Environment for Statistics. *The Annals of Statistics* **21** 1225 – 1271. https://doi.org/10.1214/aos/1176349259

[20] LIU, H., CHEN, X., WASSERMAN, L. and LAFFERTY, J. (2010). Graph-Valued Regression. In *Advances in Neural Information Processing Systems* (J. LAFFERTY, C. WILLIAMS, J. SHAWE-TAYLOR, R. ZEMEL and A. CULOTTA, eds.) **23**. Curran Associates, Inc. Available at https://papers.nips.cc/paper_files/paper/2010/hash/821fa74b50ba3f7cba1e6c53e8fa6845-Abstract.html.
[21] MATHER, J. N. (1973). Stratifications and Mappings. In *Dynamical Systems* (M. M. Peixoto, ed.) 195-232. Academic Press. https://doi.org/10.1016/B978-0-12-550350-1.50021-7

[22] MIEGHEM, P. V. (2010). *Graph Spectra for Complex Networks*. Cambridge University Press. https://doi.org/10.1017/CBO9780511921681

[23] O'NEILL, B. (1966). The fundamental equations of a submersion. **13** 459–469. https://doi.org/10.1307/mmj/1028999604

[24] OpenFlights Database. Available at https://openflights.org/data.html. Accessed the 2023-05-16.

[25] PENNEC, X. (2018). Barycentric subspace analysis on manifolds. **46** 2711–2746. https://doi.org/10.1214/17-AOS1636

[26] ROHÉ, M.-M., SERMESANT, M. and PENNEC, X. (2018). Low-dimensional representation of cardiac motion using Barycentric Subspaces: A new group-wise paradigm for estimation, analysis, and reconstruction. *Medical Image Analysis* **45** 1–12. https://doi.org/10.1016/j.media.2017.12.008

[27] SEVERN, K. E., DRYDEN, I. L. and PRESTON, S. P. (2021). Non-parametric regression for networks. *Stat* **10** e373. https://doi.org/10.1002/sta4.373

[28] SEVERN, K. E., DRYDEN, I. L. and PRESTON, S. P. (2022). Manifold valued data analysis of samples of networks, with applications in corpus linguistics. *The Annals of Applied Statistics* **16** 368–390. https://doi.org/10.1214/21-AOAS1480

[29] SIMONOVSKY, M. and KOMODAKIS, N. (2018). GraphVAE: Towards Generation of Small Graphs Using Variational Autoencoders. In *Artificial Neural Networks and Machine Learning – ICANN 2018* (V. KÚRKOVÁ, Y. MANOLOPOULOS, B. HAMMER, L. ILIADIS and I. MAGLOGIANNIS, eds.) 412–422. pringer International Publishing. https://doi.org/10.1007/978-3-030-01418-6_41

[30] SIMPSON, S. L., LYDAY, R. G., HAYASAKA, S., MARSH, A. P. and LAURIENTI, P. J. A permutation testing framework to compare groups of brain networks. *Frontiers in Computational Neuroscience* **7**. https://doi.org/10.3389/fncom.2013.00171

[31] SPIELMAN, D. A. Spectral and Algebraic Graph Theory. Available at http://cs-www.cs.yale.edu/homes/spielman/sagt/.

[32] SPORNS, O. (2018). Graph theory methods: applications in brain networks. *Dialogues in Clinical Neuroscience* **20** 111–121. https://doi.org/10.31887/DCNS.2018.20.2/osporns

[33] TSUDA, K. and SAIGO, H. Graph Classification. In *Managing and Mining Graph Data* (C. C. Aggarwal and H. Wang, eds.) 337–363. Springer US. https://doi.org/10.1007/978-1-4419-6045-0_11

[34] TUMPACH, A. B. and PRESTON, S. C. (2023). Three Methods to Put a Riemannian Metric on Shape Space. In *Geometric Science of Information* (F. NIELSEN and F. BARBARESCO, eds.) 366–374. Springer Nature Switzerland. https://doi.org/10.1007/978-3-031-38271-0_1

[35] VIGNAC, C., KRAWCZUK, I., SIRAUDIN, A., WANG, B., CEVHER, V. and FROSSARD, P. (2022). Digress: Discrete denoising diffusion for graph generation. Available at arXivpreprintarXiv:2209.14734.

[36] VON FERBER, C., HOLOVATCH, T., HOLOVATCH, Y. and PALCHYKOV, V. (2009). Public transport networks: empirical analysis and modeling. **68** 261–275. https://doi.org/10.1140/epjb/e2009-00090-x

[37] WANG, H. and MARRON, J. S. (2007). Object oriented data analysis: Sets of trees. *The Annals of Statistics* **35** 1849–1873. https://doi.org/10.1214/009053607000000217

[38] YOUNES, L. (2010). *Shapes and Diffeomorphisms* **171**. Springer Berlin, Heidelberg. https://doi.org/10.1007/978-3-642-12055-8