# Aligning the Unseen in Attributed Graphs: Interplay between Graph Geometry and Node Attributes Manifold

Aldric Labarthe, Roland Bouffanais, Julien Randon-Furling

► To cite this version:

Aldric Labarthe, Roland Bouffanais, Julien Randon-Furling. Aligning the Unseen in Attributed Graphs: Interplay between Graph Geometry and Node Attributes Manifold. 2026. hal-05482618

HAL Id: hal-05482618

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

Preprint submitted on 29 Jan 2026

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 CC BY-SA 4.0 - Attribution - ShareAlike - International
License
# Aligning the Unseen in Attributed Graphs:
## Interplay between Graph Geometry and Node Attributes Manifold

Aldric Labarthe¹² Roland Bouffanais² Julien Randon-Furling¹

## Abstract

The standard approach to representation learning on attributed graphs—i.e., simultaneously reconstructing node attributes and graph structure—is geometrically flawed, as it merges two potentially incompatible metric spaces. This forces a destructive alignment that erodes information about the graph’s underlying generative process. To recover this lost signal, we introduce a custom variational autoencoder that separates manifold learning from structural alignment. By quantifying the metric distortion needed to map the attribute manifold onto the graph’s Heat Kernel, we transform geometric conflict into an interpretable structural descriptor. Experiments show our method uncovers connectivity patterns and anomalies undetectable by conventional approaches, proving both their theoretical inadequacy and practical limitations.

## 1. Introduction

Graphs and metric spaces share many fundamentally interwined properties. During the past two decades, numerous works have been conducted to construct a good embedding for networks (Perozzi et al., 2014; Grover & Leskovec, 2016; Hamilton et al., 2017; Liu et al., 2018; Tang et al., 2015). While early research focused on simple, unattributed graphs, modern network datasets carry significant auxiliary information in the form of node and edge labels (Yang et al., 2020; 2022). Representation learning algorithms now face a fundamentally different task than originally intended: they must not only capture network connectivity but also align it with the geometry of node attributes.

This challenge, often treated as a mere technical hurdle actually points to a critical theoretical oversight. Algorithms seeking such representations effectively assume that the graph and the attributes impose no conflicting geometric constraints. This view is disjoint from reality: random geometric graph theory (Duchemin & De Castro, 2023) demonstrates that graph connectivity implies geometry (especially for social networks), just as point clouds form intrinsic manifolds. By ignoring this, models default to an assumption of *geometric alignment*, presupposing that the graph is generated by a homophilic process strictly compatible with the attribute manifold.

We argue that this assumption is a potentially critical oversight. In many real-world networks exhibiting heterophily, the graph connectivity deviates significantly from the geometry of the node attributes. Current algorithms treat this deviation as reconstruction error or noise to be smoothed over. This central assumption has never been verified, and the interaction between the two geometries remains unexplored.

In this work, we reverse this perspective. Instead of forcing alignment, we formalize and measure the *geometric misalignment* between the attribute manifold and the graph structure. As a diagnostic tool, we design a custom framework based on a Variational Autoencoder (VAE) and a two-phase training approach: first, it learns the intrinsic manifold of the node attributes; second, it deforms this manifold to accommodate the graph’s geometry. Our framework should be understood as a probe: it intentionally restricts model flexibility to reveal geometric contradictions that standard architectures are instead designed to hide. The degree of deformation required to reconcile the attribute manifold with the graph geometry thus provides a powerful and interpretable signal.

We demonstrate that this misalignment is not random noise but a structural descriptor that highlights anomalies and non-homophilic connectivity.

**Our contributions are as follows:**

* **Formalization of Geometric Misalignment:** We introduce a theoretical framework that characterizes the

**Contributions:** AL conceived, formalized, conducted the experiments and wrote the manuscript. RB and JRF reviewed and supervised. ¹Université Paris Saclay, Université Paris Cité, ENS Paris Saclay, CNRS, SSA, INSERM, Centre Borelli, F-91190, Gif-sur-Yvette, France ²Department of Computer Science, University of Geneva, Route de Drize 7, CH-1227 Carouge, Switzerland. Correspondence to: Aldric Labarthe <aldric.labarthe@ens-paris-saclay.fr>.

Preprint. January 28, 2026.

Code and experiments available on a Github repository.
conflict between the Riemannian metric induced by node attributes and the diffusion metric induced by the graph kernel.

* **Diagnostic Probe:** We propose a VAE-based diagnostic protocol based, in which manifold learning is decoupled from structural alignment, thus allowing us to isolate the specific distortions required to match the graph topology. We implement it via novel approximations to Riemannian geodesics that allow gradients to flow through the metric tensor, enabling an end-to-end optimization of geometric misalignment constraints.

* **Misalignment as a Signal:** We demonstrate that the local curvature deformation required to align these geometries is a highly informative graph feature, improving community detection and identifying structural anomalies in both synthetic and real-world datasets.

# 2. Problem Statement

## 2.1. The geometric misalignment problem

Let $G = (V, E)$ be a weighted attributed graph with $N = |V|$ nodes each associated with an $\mathbb{R}^D$ attribute vector. The graph is classically described by two matrices:

* **Adjacency Matrix:** $A \in \mathbb{R}_+^{N \times N}$, where $A_{ij} > 0$ iff $(v_i, v_j) \in E$. We assume the graph is undirected and has no self-loops.

* **Attribute Matrix:** $X \in \mathcal{X}^N \subseteq \mathbb{R}^{N \times D}$, where the row $x_i$ represents the feature vector of node $v_i$.

Standard approaches for the attributed graph embedding problem seek to learn a low-dimensional latent representation $z_i \in \mathcal{Z} \subseteq \mathbb{R}^d$ of node $v_i$. The objective is to preserve information from both $X$ and $A$. This is typically formalized by assuming the existence of two decoding functions:

* An attribute decoder $\mathfrak{M}: \mathcal{Z} \to \mathbb{R}^D$ such that $\mathfrak{M}(z_i) = x_i$.

* A structural kernel $K: \mathcal{Z} \times \mathcal{Z} \to \mathbb{R}$ such that $K(z_i, z_j) = A_{ij}$.

**Assumption 1 (Manifold hypothesis).** Node attributes $x_i$ live on a lower dimensional manifold, i.e., $\mathcal{X}$ is a $d$-submanifold immersed in $\mathbb{R}^D$ (with $d \ll D$ and $d \ll N$).

Provided that $\mathfrak{M}$ is $C^1$ and $d\mathfrak{M}$ is injective, the attribute decoder $\mathfrak{M}$ induces a Riemannian metric $g^{\mathfrak{M}}$ on $\mathcal{Z}$ (pullback metric).

By the Moore-Aronszajn theorem, the structural kernel $K$ implies the existence of a unique reproducing kernel Hilbert space (RKHS) $\mathcal{H}_K$ and a feature map $\Phi_K: \mathcal{Z} \to \mathcal{H}_K$ from

which we can pull back the standard inner product of $\mathcal{H}_K$ onto the latent space $\mathcal{Z}$. This induces a Riemannian metric tensor $g^K$ on $\mathcal{Z}$:

$$g_{ij}^{K}(z) = \frac{\partial^2 K(u, v)}{\partial u_i \partial v_j} \Big|_{u=v=z}.$$

Consequently, the latent space $\mathcal{Z}$ is subject to two competing geometric and potentially conflicting structures.

**Problem 1.** We seek a single representation $\mathcal{Z}$ that faithfully embeds the data; however, the graph structure imposes a Riemannian metric $g^K$, while the node attributes impose a metric $g^{\mathfrak{M}}$ on the same manifold.

In this paper, we characterize this metric misalignment and investigate its utility as a structural descriptor of the graph.

## 2.2. Kernel influence

The choice of the structural kernel $K$ dictates the geometric nature of this misalignment and thus requires further attention.

Initial approaches (Kipf & Welling, 2016; Salha et al., 2019), building on a foundational intuition from the social network analysis (Hoff et al., 2002), used the canonical scalar product. This choice imposes a flat Euclidean geometry on the latent space ($g^K = \text{Id}$). However, this rigidity often conflicts with the manifold hypothesis: if the intrinsic geometry of the attributes requires a curved space to be represented faithfully in $d$ dimensions, forcing the latent space to be Euclidean induces significant distortion. Consequently, the model cannot reconcile the flat geometry imposed by the kernel with the intrinsic curvature of the data manifold, resulting in unavoidable reconstruction loss.

To address this persistent reconstruction gap, some authors have adopted bilinear kernels (Trouillon et al., 2016), such as $K(z_i, z_j) = z_i^T W z_j$. Consequently, $g^K = W$ implies that this kernel will induce a trade-off for a highly heterogeneously curved attribute manifold between properly reconstructing the attribute manifold or converging to a mean curvature manifold.

Another possible choice is the so-called Heat Kernel, which is the solution to the heat equation, often denoted $H(t, x, y)$, representing the temperature at point $x$ and time $t$ given an initial unit of heat at $y$ at time $t = 0$:

$$H(u, v) = \sum_{t \in \mathcal{T}} (4\pi t)^{-\frac{d}{2}} \exp \left( -\frac{d^2 g^{(\mathfrak{M})}(u, v)}{4t} \right). \quad (1)$$

The choice of time scales $\mathcal{T}$ explicitly governs the geometry induced by the kernel $H$. In the short-time limit ($t \to 0$), the Varadhan asymptotic formula (Varadhan, 1967) guarantees that the induced metric $g^K$ converges to $g^{\mathfrak{M}}$, effectively
unifying the two latent geometries. This construction enables us to address the misalignment problem by enforcing the structural kernel as a direct function of the attribute geometry.

When larger time scales are employed, $g^K$ transitions from the Riemannian metric to a diffusion metric. Crucially, on compact manifolds, the mapping from the original manifold to the diffusion geometry is a diffeomorphism. Since this transformation preserves topology and smoothness, the attribute decoder $\mathfrak{M}$ can effectively learn to invert this metric warping to reconstruct $x_i$, provided $t$ is not too large to cause spectral truncation (i.e., oversmoothing).

However, this geometric unification imposes a strong inductive bias, specifically equating structural connectivity with attribute proximity. By forcing the metric tensors to align ($g^K \approx g^{\mathfrak{M}}$), the model inherently assumes the graph is generated by a homophilic process on the attribute manifold. Consequently, a structural misalignment may still persist if the observed adjacency matrix $A$ exhibits heterophily or significant noise, as the model will treat edges that violate the attribute geometry as anomalies.

Therefore, this choice of kernel fundamentally reframes the misalignment from an optimization challenge to an informative structural descriptor. By fixing the structural kernel to be a function of the attribute geometry, the reconstruction error of the adjacency matrix no longer solely represents model incapacity, but rather quantifies the deviation of the observed graph from the homophilic ideal. A low reconstruction error implies that the graph connectivity is intrinsic to the attribute manifold, while a high residual signals that the graph structure is driven by external factors or heterophilic interactions not captured by the attribute diffusion. Thus, the magnitude of this residual misalignment serves as a direct, quantifiable proxy for the degree of attribute homophily in the graph. The remainder of this paper introduces an algorithm to evaluate this misalignment gap.

Note that in general settings, this geometric misalignment has been intuitively circumvented by introducing a learnable transformation $f: \mathcal{Z} \to \mathcal{Z}'$ prior to the inner product, such that the structural kernel becomes $K(z_i, z_j) = \langle f(z_i), f(z_j) \rangle$. This approach, ubiquitous in models employing MLP-based decoders, effectively absorbs the metric discrepancy by learning a diffeomorphism that warps the attribute-aligned latent space into a geometry compatible with the graph topology. However, this flexibility creates a fundamental identification challenge: the intrinsic geometry of the data becomes indistinguishable from the arbitrary distortions introduced by the decoder $f$. Furthermore, this resolution is often incomplete; since standard inner products impose a flat Euclidean signature, no continuous function $f$ can fully resolve the misalignment if the graph exhibits non-Euclidean topology (e.g., hyperbolicity) or if the attribute

manifold requires topological tearing to match the graph structure.

# 3. Methods

## 3.1. VAE for manifold learning

Under the manifold hypothesis, we seek to recover an underlying data manifold $\mathcal{X}$ using a variational autoencoder design (Kingma & Welling, 2013; Dilokthanakul et al., 2017). We provide a refresher on VAE and manifold learning through VAE in Appendix A. The encoder gets as inputs the nodes attributes vectors $x_i$ and outputs, for each individual vector a latent code $z_i$. The variational decoder is at the node level as its task is only to reconstruct each node's attribute vectors from its latent code. This approach allows us to perform our manifold learning objective.

We assume that the $x_i$ vectors are noisy samples drawn from the manifold $\mathcal{X}$, i.e., there exists $z_i \in \mathcal{Z}$ such that $x_i = \mathfrak{M}(z_i)$.

Through the VAE, we estimate a probabilistic model made of a generative process $p(x_i|z_i)$ and a posterior distribution $q(z_i|x_i)$. $p(x_i|z_i)$ is estimated using a variational decoder made of a multi-layer perceptron with two heads $\psi_{\theta}^{(\mu)}: \hat{\mathcal{Z}} \to \mathcal{X}$ and $\psi_{\theta}^{(\sigma)}: \hat{\mathcal{Z}} \to \mathbb{R}$ such that:

$$ \hat{p}_{\theta} = \mathcal{N} \left( \psi_{\theta}^{(\mu)}(z_i), \left( \psi_{\theta}^{(\sigma)}(z_i) \right)^2 \cdot \mathrm{Id}_D \right). $$

Provided that $\psi_{\theta}^{(\mu)}$ and $\psi_{\theta}^{(\sigma)}$ are differentiable on $\mathcal{Z}$ and their differential is of full rank, we have:

**Definition 1 (Estimated Attribute Geometry).** Let $\psi_{\theta}^{(\mu)}: \hat{\mathcal{Z}} \to \mathcal{X}$ and $\psi_{\theta}^{(\sigma)}: \hat{\mathcal{Z}} \to \mathbb{R}$ be the two heads of a smooth decoder parameterized by $\theta$. The latent space $\hat{\mathcal{Z}}$ is equipped with a Riemannian metric $\hat{g}_{\mathfrak{M}}$, defined as the pullback of the Euclidean metric of $\mathbb{R}^D$ via $\psi_{\theta}$. For any point $z \in \hat{\mathcal{Z}}$, the metric tensor is given by:

$$ \hat{g}_{\mathfrak{M}} = (J_{\theta}^{(\psi^{\mu})})^{\top} (J_{\theta}^{(\psi^{\mu})}) + (J_{\theta}^{(\psi^{\sigma})})^{\top} (J_{\theta}^{(\psi^{\sigma})}). $$

Our goal now is to *inject* the kernel geometry into this attribute geometry. Using the estimated attribute geometry through $\hat{g}_{\mathfrak{M}}$, we derive for any pair $(z_i, z_j) \in \hat{\mathcal{Z}}^2$ an approximate geodesic distance function $\hat{d}_{\hat{g}_{\mathfrak{M}}}(z_i, z_j|\theta)$, and estimate the heat kernel on the manifold $\hat{H}_{\mathcal{M}}(t|\theta)$ using the Varadhan asymptotic (Eq. (1)). The final heat kernel $\hat{H}_{\mathcal{M}}(\theta)$ is the sum of the heat kernels over the set of time scales $\mathcal{T}$, i.e., $\sum_{t \in \mathcal{T}} \hat{H}_{\mathcal{M}}(t|\theta)$.

**Definition 2 (Kernel Alignment Objective).** The Kernel Alignment Objective is defined as the minimization of the distance between the adjacency matrix $A$ of the graph and
<figure><img src="image_2.png" /><figcaption>Figure 1. Geometric Alignment VAE. (Black) Phase 1 learns the intrinsic attribute manifold $\mathcal{M}$ via standard reconstruction. (Blue) Phase 2 freezes the encoder and deforms the manifold geometry to align the geometry of $\hat{\mathcal{M}}$ with the kernel geometry.</figcaption></figure>

the kernel computed on the latent space $\hat{\mathcal{Z}}$:

$$ \mathcal{L}_2 = \| A - \hat{K}_{\hat{\mathcal{M}}}(\theta) \|_2^2. \qquad (2) $$

Minimizing this loss acts as a geometric regularization on the decoder's Jacobian, forcing $\hat{g}_{\mathcal{M}}$ to align with the diffusion geometry induced by the graph ($g_K$).

Since the Heat Kernel is a direct function of the attribute metric, the magnitude of the local distortion required to align $\hat{g}_{\mathcal{M}}$ with the kernel geometry $g^K$ serves as a proxy for attribute homophily in each neighborhood. Effectively, by enforcing this alignment, we recover the geometry that would define the manifold if the graph were generated by a purely homophilic process. Analyzing the resulting metric deformation thus reveals exactly where this assumption fails. This fundamentally reframes geometric misalignment: it is no longer a conceptual barrier, but an interpretable signal of homophily.

## 3.2. Training Algorithm

We propose a two-phase training procedure (Fig. 1) designed to decouple manifold learning from geometric alignment.

**Phase 1: Manifold Learning.** First, we learn a parametric approximation of the attribute manifold by training a standard VAE, ignoring graph topology. The encoder parameterization $\phi$ and decoder parameterization $\theta$ are optimized solely to reconstruct attributes $x_i$, minimizing the ELBO (see Appendix B.1). This yields a latent space $\mathcal{Z}$ with a Riemannian metric $g_{\mathcal{M}}$ induced exclusively by the attribute reconstruction task (and the variational prior).

**Phase 2: Geometric Alignment.** In the second phase, we freeze the encoder parameters $\phi$, fixing the latent coordinates $z_i$ and their variational distributions. We then replace

the reconstruction loss with the Kernel Alignment Objective (Eq. (2)).

By freezing the encoder, we establish a fixed coordinate system determined solely by the attribute manifold. This prevents the model from minimizing the loss trivially by rearranging latent codes to cluster connected nodes (as in standard graph embedding). Instead, the decoder is forced to minimize $\mathcal{L}_2$ by locally warping the metric tensor $g_{\mathcal{M}}$ to match the graph's Heat Kernel. This design avoids the latent space entanglement observed in simultaneous multi-decoder designs (Yang et al., 2020), where the model arbitrarily balances structural and attribute constraints. Here, the coordinate system remains stable, allowing the geometric deformation to be measured explicitly.

**Differentiable Geodesic Computation.** Minimizing Eq. (2) requires computing pairwise geodesic distances $\hat{d}_g(z_i, z_j)$ at every training step, a process typically intractable for auto-differentiation due to the computational cost of solving the boundary value problem ($O(d^3)$ for tensor inversion, $O(d^2D)$ for Jacobian computation). Unlike post-hoc analysis methods (Arvanitidis et al., 2018; Kalatzis et al.), our loss requires gradients to flow through these distances to update the decoder.

We resolve this by introducing a differentiable discrete approximation. For low-dimensional spaces ($d \le 3$), we discretize the latent domain into a grid, caching the metric tensor for each cell, and compute shortest paths via Dijkstra's algorithm. Crucially, we decouple the path selection from the weight evaluation: we fix the discrete path found by Dijkstra's algorithm and backpropagate solely through the sum of edge weights along that path. To our knowledge, this is the first approach to enable end-to-end geometric latent space regularization during VAE training.

However, as grid complexity scales exponentially with di-
mension ($O(P^d)$), we employ a linear interpolation estimator for high-dimensional settings. This method approximates the geodesic by integrating the metric along the linear segment connecting $z_i$ and $z_j$. While theoretically an upper bound, we find this estimator to be highly correlated with true geodesics in practice, likely because VAE latent manifolds typically exhibit simple curvature obstacles that do not require complex path deviations. Full details are provided in **Appendix B.2.** In **Section 4.1**, we present results obtained with synthetic data using both estimators, thus highlighting their practical equivalence despite the theoretical superior fidelity of the grid-based approach. For instance, we substantiated this claim on the synthetic manifold, where the linear-interpolation estimator achieved near-perfect alignment with an exact numerical solver ($r = 0.999, \rho = 0.999$ and a percentage error of only 1.0%, with a bias of +0.015).

# 4. Results

Standard graph representation learning literature typically aims to maximize predictive performance on downstream tasks (e.g., node classification), effectively treating the divergence between graph topology and node attributes as noise to be smoothed over. Our work reverses this paradigm: we posit that this divergence is not noise, but a meaningful structural signal—the geometric misalignment. Consequently, the following experiments are designed to demonstrate that this misalignment is not only a theoretical issue but can also lead to valuable information being missed. Our architecture is not designed to minimize a classification error rate but to recover and interpret this geometric tension.

## 4.1. Synthetic Dataset

### 4.1.1. EXPERIMENTAL SETUP

We generate 500 vectors $x_i$ of dimension $D = 20$ that live on a Swiss-roll manifold of intrinsic dimension $d = 2$. From these vectors, we generate a preferential attachment-based social network using an existing generator (Labarthe & Kerzreho, 2026).

The generator creates preferential attachment via a compatibility matrix. We construct it using a Gaussian similarity kernel on the true $z_i$ (the inner-coordinates of the $x_i$). This incorporates the manifold's geometry into the network's generating process.

We introduce a perturbation in our similarity matrix. Our objective is to reproduce settings in which a small subset of vertices is corrupted (e.g., not homogeneous with other vertices in their underlying homophily). To do so, we select 70 nodes whose attribute vectors are close on the attribute manifold, and we shuffle their intra-group pairwise scores. This way, we preserve the distribution of scores, but break the link between the scores and the manifold. Following our

previous assumption, if we were not to perturb the similarity matrix, the information in the network topology would be redundant to some extent with that in the node attributes manifold. Indeed, similarity scores are fully determined by manifold structure and topology. Perturbing the matrix would induce neighborhoods where the network connectivity is not related to the attribute manifold. If our assumption is correct, we should observe a discrepancy in the effect of phase 2 between neighborhoods with fully homophilic connectivity and perturbed neighborhoods. This perturbation strategy creates two communities with similar connectivity structures that are nearly indistinguishable to standard algorithms.

### 4.1.2. MANIFOLD RECOVERY (PHASE 1)

First, we verify the quality of the manifold learning task in Phase 1. Our estimated geodesic distances show a strong correspondence with the theoretical geodesic distances ($R^2 = 0.894$ with Dijkstra estimation, $R^2 = 0.889$ with linear interpolation estimation), and Spearman rank correlation shows that they are highly correlated ($\rho = 0.85$ with Dijkstra estimation, $\rho = 0.84$ with linear interpolation estimation). The correlation between angles is also significant ($\rho = 0.832$).

This indicates that the manifold learned in Phase 1 is highly consistent with the theoretical manifold under investigation. Panels A and B of Fig. 2 show a consistency between the theoretical curvature and the estimated one, despite the manifold being estimated up to a diffeomorphism.

### 4.1.3. DETECTING MISALIGNMENT (PHASE 2)

In Phase 2, the graph constraint forces a deformation of the latent geometry. We quantify this deformation by computing the link-level distortion $Z_{ij}$, defined as the modified Z-score of the pairwise distance variation between Phase 1 and Phase 2:

$$Z_{ij} = \frac{0.6745 (\delta_{ij} - \tilde{\delta})}{\text{median}_{kl} (|\delta_{kl} - \tilde{\delta}|)}, \quad (3)$$

where $\delta_{ij} = \log |\hat{d}_{\text{Phase1}}(z_i, z_j) - \hat{d}_{\text{Phase2}}(z_i, z_j)|$ and $\tilde{\delta} = \text{median}_{kl} (\delta_{kl})$. We also define the *Node Geometric Distortion Score* $S_i = \sum_{j \neq i} Z_{ij}$.

The effects of Phase 2 correction on the neighborhoods of the manifold are heterogeneous. Panel D of Fig. 2 illustrates the relative variation of the Gaussian curvature at each point, while Panel C shows the latent space curvature after introducing kernel geometry. We observe an overall smoothing of the manifold, with the removal of the characteristic curvature drop line that was also a feature of the theoretical manifold. In contrast, two regions have been consistently altered: the hollow portion of the C-shaped point cloud and
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Synthetic dataset latent space. Panel A: theoretical latent space $\mathcal{Z}$; Panel B: estimated latent space after Phase 1; Panel C: estimated latent space after Phase 2; Panel D: estimated curvature changes between Phases 1 and 2. Nodes are plotted in dark purple for normal nodes, and in yellow for perturbed ones.</figcaption></figure>

the neighborhood of the two classes' separation. The first change carries no useful information: curvature was due to missing data points, leading to a high curvature region due to the variance in the head of the decoder.

<figure><img src="image_4.png" /><figcaption>**Figure 3.** Node Geometric Distortion Scores $S_i$ (from eq. 3). Non-perturbed nodes are plotted in blue, perturbed nodes in red.</figcaption></figure>

However, the second effect suggests that our hypothesis is validated: an alteration in the underlying process that generates the network does induce a geometrical alteration of the alignment between graph and attribute geometry. Using $S_i$ as an anomaly score (Fig. 3), we find that this simple metric distinguishes two separate classes with an F1-score of 82.50%. This significantly outperforms naïve baselines such as spectral clustering (F1 36.65%) and Louvain clustering (F1 44.87%). Even more sophisticated benchmarks, such as a standard VAE and VGAE (Kipf & Welling, 2016), struggle to match this performance because they cannot extract a signal as clear as ours, reaching a peak F1 of 45.71% using an Isolation Forest and VAE. Finally, using the VGAE's aggregated reconstruction error (similar to the state-of-the-art anomaly detection approach DOMINANT (Ding et al., 2019)) proved effective, though it did not reach our performance, yielding an F1 of 74.29%. Full performance metrics are summarized in Table 1. This achievement is consistent with our experimental setting: while reconstruction error flags random noise, our geometric misalignment score specifically isolates structural contradictions—edges that

exist in the graph but violate the manifold logic.

<table><caption>**Table 1.** Performance benchmark for community detection comparing our misalignment correction node-level effect ($S_i$) against clustering and VAE and VGAE based baselines. Values are in %, † is the VGAE variant with GCN decoder, and ‡ is our method computed with the linear interpolation geodesic distance estimator.</caption><thead><tr><th>Method</th><th>F1</th><th>ROC AUC</th><th>ARI</th></tr></thead><tbody><tr><td>Proposed Metric (S<sub>i</sub>)</td><td><strong>82.50</strong></td><td><strong>96.48</strong></td><td><strong>73.13</strong></td></tr><tr><td>Proposed Metric (S<sub>i</sub>, ‡)</td><td><strong>83.75</strong></td><td><strong>97.25</strong></td><td><strong>70.95</strong></td></tr><tr><td>Cosine Assortativity</td><td>34.90</td><td>64.29</td><td>4.54</td></tr><tr><td>Spectral Clustering</td><td>36.65</td><td>71.86</td><td>-3.28</td></tr><tr><td>Louvain Clustering</td><td>44.87</td><td>80.00</td><td>9.60</td></tr><tr><td>VAE K-Means</td><td>45.60</td><td>80.58</td><td>10.81</td></tr><tr><td>VAE Isolation Forest</td><td>45.71</td><td>68.44</td><td>29.45</td></tr><tr><td>VAE Recon. Error</td><td>28.16</td><td>54.87</td><td>-9.03</td></tr><tr><td>VGAE K-Means</td><td>46.51</td><td>81.28</td><td>12.33</td></tr><tr><td>VGAE Isolation Forest</td><td>28.57</td><td>58.47</td><td>12.42</td></tr><tr><td>VGAE Recon. Error</td><td>74.29</td><td>94.90</td><td>68.43</td></tr><tr><td>VGAE Recon. Error (†)</td><td>34.56</td><td>46.53</td><td>0.00</td></tr></tbody></table>

This novel and powerful approach demonstrates that geometric misalignment conceals a signal about the graph's underlying structure that can be leveraged to improve existing architectures for downstream tasks, such as anomaly detection, especially when the perturbation is very subtle. In particular, our $Z_{ij}$ coefficients shed a new light on homophily than assortativity does (see Table 1, despite being correlated as they both relate to homophily $r = 0.3193$, $\rho = 0.3288$, MI = 0.2275) by offering a true link-level measure. This also provides a practical application and a new motivation for analyzing graph-attribute manifold geometry alignment.
## 4.2. Empirical Dataset

### 4.2.1. DATASET AND HYPOTHESIS

Beyond the validation on synthetic data, we investigate the utility of our framework on a novel real-world graph constructed for this study. We analyze the transportation network of Île-de-France (Paris region, France), where nodes represent municipalities (N = 1,266) and edges represent public transit connectivity (see **Appendix C.2** for construction details). Node attributes consist of high-dimensional socio-demographic indicators.

This dataset provides an ideal testbed for examining geometric misalignment. First, the attributes satisfy the manifold hypothesis: socio-demographic indicators are highly correlated (e.g., housing density and population), implying a lower-dimensional intrinsic structure. Second, the graph topology (transit links) is partially, but not entirely, driven by these attributes. While transit often connects similar dense urban centers (homophily), it also creates “long-range” shortcuts between socio-economically distinct areas to facilitate commuting. We hypothesize that our alignment score can identify non-homophilic links, distinguishing purely geographic or demographic connections from those driven by external urban planning decisions.

### 4.2.2. LATENT STRUCTURE AND INTERPRETATION

<figure><img src="image_5.png" /><figcaption>**Figure 4.** Phase 1 Latent manifold clustering Points are colored according to the log. robust Z-score of the variation between Phase 1 and Phase 2 of the total pairwise distances with respect to other nodes.</figcaption></figure>

We train our VAE architecture following our two-phase methodology on this dataset. Figure 4 shows that at the end of Phase 1, the latent space of the VAE is structured around 5 regions, separated by curvature walls. To interpret

these latent structures, we trained a surrogate decision-tree classifier on the input features, which were pre-normalized with a standard normal distribution (Z-scores). This analysis reveals that the manifold encodes a distinct socioeconomic hierarchy: separating metropolitan centers (e.g., Paris, Cluster 3) and dense urban hubs (Cluster 2) from the residential periphery (Cluster 1) and specialized Productive Nodes (Cluster 4), while the majority of communes form a rural baseline (Cluster 0). Detailed decision rules, threshold quantiles, and representative medoids for each cluster are provided in **Appendix C.2.3**.

Following the methodology we developed for the synthetic dataset, we compute the $Z_{ij}$ scores (from eq. 3), and find that the five identified regions respond differently to our Phase 2 correction.

<figure><img src="image_6.png" /><figcaption>**Figure 5.** Cluster Interaction Network. Nodes represent the five socioeconomic functional regions identified in Phase 1 (0: Rural, 1: Residential, 2: Dense Urban, 3: Metro Center, 4: Productive). Edges are colored by the Average Manifold-Graph Distortion (modified Z-Score), measuring the tension between the attribute manifold and the transport topology. Red edges (positive Z) identify structural “shortcuts” where transit reduces the distances between points placed according to socio-demographic proximity.</figcaption></figure>

Figure 5 aggregates these pairwise variations into a cluster-level interaction map. The metropolitan center (Cluster 3) acts as the primary source of geometric tension, exhibiting extreme positive deviations both internally (avg. Z-score of 16.2) and in its connections to dense urban hubs (Cluster 2) and productive nodes (Cluster 4). This validates our hypothesis: the public transit system creates efficient “shortcuts” into the capital that bridge socioeconomically distinct areas, a structural feature the attribute manifold alone does not
capture. Conversely, the rural baseline (Cluster 0) and the residential periphery (Cluster 1) are dominated by negative Z-scores (e.g., -3.8 and -2.7 self-interaction), indicating that transit connectivity in these outer rings aligns closely with—or is even sparser than—the intrinsic attribute similarity and geographic proximity.

These geometric distortions provide a quantitative signature of the region's monocentric transport topology, corroborating well-established critiques in urban geography regarding the “radial” nature of the so-called Francilien network (Halbert, 2008; 2015). The extreme positive variations associated with the metropolitan center (Cluster 3)—particularly the thick red edges connecting it to Clusters 2 and 4—capture the efficiency of the hub-and-spoke system: the network effectively shrinks the distance between the capital and its productive satellites, rendering them “closer” in the transport graph than their socio-spatial attributes would suggest. In sharp contrast, the deep blue self-loops of the rural baseline (Cluster 0) and the residential Periphery (Cluster 1) (Z-score < 2.7) reveal the systemic failure of lateral connectivity. While these municipalities are often geographically adjacent and socio-demographically indistinguishable, the lack of direct banlieue-to-banlieue transit links forces disjointed, inefficient paths. Our metric thus successfully recovers the accessibility gap driven by the historic star-shaped topology of the Greater Paris RER/Metro system, effectively quantifying the phenomenon where the center exerts a powerful gravitational pull while peripheral towns remain structurally isolated from their neighbors (Wenglen-ski & Orfeuil, 2004).

Finally, we validate the utility of our framework by contrasting it with direct statistical baselines. We attempted to regress the graph transportation convenience directly on node attributes ($R^2$ = 0.003) and on pairwise similarities ($R^2$ = 0.0031). The lack of correlation underscores the core intuition of our methodology: some information lies in the heterogeneous misalignment between attribute geometry and network connectivity. Here, this enabled us to make explicit the structural inequities of the transport network that otherwise remain effectively unprovable from our data alone.

# 5. Discussion

**Towards Geometry-Aware Architectures.** Most attributed graph representation learning methods optimize for the alignment of graph geometry and node attributes, implicitly assuming that the corresponding metric structures are compatible. When this assumption fails, the optimization does not resolve the incompatibility but removes it, systematically discarding geometric information. In this sense, existing joint embedding architectures are mathematically mis-specified: they enforce metric equivalence rather than

model metric tension. This work exposes this modeling error and shows that the suppressed discrepancies correspond to structured, interpretable signals, motivating a shift toward architectures that explicitly reason about geometric incompatibility rather than smoothing it away under homophilic assumptions.

**Geometric misalignment as a signal.** In this work, we hypothesize that the discrepancy between node attributes and graph structure is a feature, not a bug. While standard representation learning algorithms typically smooth over this discrepancy, our framework explicitly measures the deformation required to map one modality onto the other. By quantifying the deformation needed to align graph diffusion distances with the attribute manifold, the resulting misalignment scores ($Z_{ij}$) identify structural shortcuts—connections that bridge distant regions of feature space. Beyond the transport network studied here, this perspective applies broadly: it distinguishes connectivity driven by intrinsic similarity from connectivity driven by functional or relational constraints, as in biological networks with non-chemical interactions or social graphs where ties cross demographic boundaries. Therefore, this paper argues that not only does a conceptual issue exist in existing architectures, but that this issue also carries useful information, making it relevant to practical applications.

**Scope and scalability.** Our VAE framework is deliberately not proposed as a scalable or predictive architecture. Instead, its primary role is diagnostic: to isolate geometric incompatibilities in a controlled setting and to quantify the information lost when standard models enforce joint alignment. While our geodesic computation step is computationally intensive, this cost is the necessary price of interpretability. Specifically, we trade the ability to process millions of nodes for the ability to quantify the precise structural signals encoded within. As such, the framework is best viewed as an analytical tool, applicable to moderate-sized graphs or subsampled regions of larger systems, in the same spirit as spectral or diffusion-based exploratory methods.

**Limitations and Future Work.** While our grid-based and linear interpolation estimators suffice for the experiments presented, they highlight a critical gap in the geometric deep learning toolkit: the need for scalable, differentiable geodesic estimators. Unlike recent developments in post-hoc analysis, which rely on non-differentiable solvers (e.g., discrete graph searches or ODE solvers), training geometric regularizers end-to-end requires gradients to propagate through the distance computation itself.

# Impact Statement

This paper presents work whose goal is to advance the field of Machine Learning. There are many potential societal
consequences of our work, none of which we feel must be
specifically highlighted here.

References

Arvanitidis, G., Hansen, L. K., and Hauberg, S. Latent space oddity: On the curvature of deep generative models. In *International Conference on Learning Representations (ICLR)*, 2018. URL https://openreview.net/forum?id=SJzRZ-WCZ.

Davidson, T. R., Falorsi, L., Cao, N. D., Kipf, T., and Tomczak, J. M. Hyperspherical variational auto-encoders. *arXiv preprint arXiv:1804.00891*, (arXiv:1804.00891), 2018. doi: 10.48550/arXiv.1804.00891. URL http://arxiv.org/abs/1804.00891.

Dilokthanakul, N., Mediano, P. A. M., Garnelo, M., Lee, M. C. H., Salimbeni, H., Arulkumaran, K., and Shanahan, M. Deep unsupervised clustering with gaussian mixture variational autoencoders. *arXiv preprint*, (arXiv:1611.02648), 2017. doi: 10.48550/arXiv.1611.02648. URL http://arxiv.org/abs/1611.02648.

Ding, K., Li, J., Bhanushali, R., and Liu, H. Deep anomaly detection on attributed networks. In *Proceedings of the 2019 SIAM International Conference on Data Mining (SDM)*, pp. 594–602. SIAM, 2019. doi: 10.1137/1.9781611975673.67. URL https://doi.org/10.1137/1.9781611975673.67.

Duchemin, Q. and De Castro, Y. Random geometric graph: Some recent developments and perspectives. *High Dimensional Probability IX: The Ethereal Volume*, pp. 347–392, 2023.

Grover, A. and Leskovec, J. node2vec: Scalable feature learning for networks. In *Proceedings of the 22nd ACM SIGKDD international conference on Knowledge discovery and data mining*, pp. 855–864, 2016.

Halbert, L. Examining the mega-city-region hypothesis: Evidence from the paris city-region/bassin parisien. *Régional studies*, 42(8):1147–1160, 2008.

Halbert, L. *L'avantage métropolitain*. PUF, 2015.

Hamilton, W., Ying, Z., and Leskovec, J. Inductive representation learning on large graphs. *Advances in neural information processing systems*, 30, 2017.

Hoff, P. D., Raftery, A. E., and Handcock, M. S. Latent space approaches to social network analysis. *Journal of the american Statistical association*, 97(460):1090–1098, 2002.

Kalatzis, D., Eklund, D., Arvanitidis, G., and Hauberg, S. Variational autoencoders with riemannian brownian motion priors. *arXiv preprint*, (arXiv:2002.05227). doi: 10.48550/arXiv.2002.05227. URL http://arxiv.org/abs/2002.05227.

Kingma, D. P. and Welling, M. Auto-encoding variational bayes. *arXiv preprint*, (arXiv:1312.6114), 2013. doi: 10.48550/arXiv.1312.6114. URL http://arxiv.org/abs/1312.6114.

Kipf, T. N. and Welling, M. Variational graph auto-encoders, 2016. URL http://arxiv.org/abs/1611.07308.

Labarthe, A. and Kerzreho, Y. Generating social networks with static and dynamic utility-maximization approaches. *Social Network Analysis and Mining*, 2026. URL https://arxiv.org/abs/2411.16464. To appear.

Liu, J., He, Z., Wei, L., and Huang, Y. Content to node: Self-translation network embedding. In *Proceedings of the 24th ACM SIGKDD international conference on knowledge discovery & data mining*, pp. 1794–1802, 2018.

Miolane, N. and Holmes, S. Learning weighted submanifolds with variational autoencoders and riemannian variational autoencoders. In *Proceedings of the IEEE/CVF conference on computer vision and pattern recognition*, pp. 14503–14511, 2020.

Perozzi, B., Al-Rfou, R., and Skiena, S. DeepWalk: Online learning of social representations. In *Proceedings of the 20th ACM SIGKDD international conference on Knowledge discovery and data mining*, pp. 701–710, 2014. doi: 10.1145/2623330.2623732. URL http://arxiv.org/abs/1403.6652.

Que, Q. and Belkin, M. Back to the future: Radial basis function networks revisited. In *Artificial intelligence and statistics*, pp. 1375–1383. PMLR, 2016.

Salha, G., Hennequin, R., and Vazirgiannis, M. Keep it simple: Graph autoencoders without graph convolutional networks. *arXiv preprint arXiv:1910.00942*, 2019.

Shao, H., Kumar, A., and Thomas Fletcher, P. The riemannian geometry of deep generative models. In *Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition Workshops*, pp. 315–323, 2018.

Syrota, S., Zainchkovskyy, Y., Xi, J., Bloem-Reddy, B., and Hauberg, S. Identifying metric structures of deep latent variable models. *arXiv preprint*, (arXiv:2502.13757), 2025.

Tang, J., Qu, M., and Mei, Q. Pte: Predictive text embedding through large-scale heterogeneous text networks.
In *Proceedings of the 21st ACM SIGKDD international conference on knowledge discovery and data mining*, pp. 1165–1174, 2015.

Tosi, A., Hauberg, S., Vellido, A., and Lawrence, N. D. *Metrics for probabilistic geometries*. *arXiv preprint*, (arXiv:1411.7432), 2014.

Trouillon, T., Welbl, J., Riedel, S., Gaussier, É., and Bouchard, G. *Complex embeddings for simple link prediction*. In *International conference on machine learning*, pp. 2071–2080. PMLR, 2016.

Varadhan, S. R. S. On the behavior of the fundamental solution of the heat equation with variable coefficients. *Communications on Pure and Applied Mathematics*, 20 (2):431–455, 1967.

Wenglenski, S. and Orfeuil, J.-P. *Differences in accessibility to the job market according to social status and place of residence in the paris area*. *Built Environment*, 30(2): 116–126, 2004.

Yang, C., Zhang, J., Wang, H., Li, S., Kim, M., Walker, M., Xiao, Y., and Han, J. *Relation learning on social networks with multi-modal graph edge variational autoencoders*. In *Proceedings of the 13th International Conference on Web Search and Data Mining*, pp. 699–707, 2020.

Yang, C., Xiao, Y., Zhang, Y., Sun, Y., and Han, J. *Heterogeneous network representation learning: A unified framework with survey and benchmark volume 34*, pp. 4854–4873, 2022. doi: 10.1109/TKDE. 2020.3045924. URL https://ieeexplore.ieee. org/document/9300240. Conference Name: IEEE Transactions on Knowledge and Data Engineering.

Zheng, Y., He, T., Qiu, Y., and Wipf, D. P. *Learning manifold dimensions with conditional variational autoencoders*. *Advances in Neural Information Processing Systems*, 35:34709–34721, 2022.
## A. Prior on Variational Autoencoders (VAE)

A VAE is a generative model that approximates the data-generating distribution $p(x)$ by introducing a latent variable $z \in \mathcal{Z}$, and learning a probabilistic mapping between observations $x \in \mathcal{X}$ and the latent space. Formally, it consists of:

1. Prior distribution over latent variables $p(z)$

2. Generative model (decoder): A likelihood model $p_{\theta}(x | z)$, parameterized by $\theta$, which defines the conditional distribution of data given the latent variables.

3. Variational posterior (encoder): An approximate posterior $q_{\phi}(z | x)$, parameterized by $\phi$, typically modeled as a Gaussian with diagonal covariance:

$$q_{\phi}(z | x) = \mathcal{N}(z | \mu_{\phi}(x), \text{diag}(\sigma_{\phi}^{2}(x)))$$

where $\mu_{\phi}(x)$ and $\sigma_{\phi}(x)$ are outputs of the encoder network.

The objective of a VAE it to train an encoder network to output the mean and variance of a distribution on the latent space, and a decoder network to reconstruct the input by outputting a mean and variance for a distribution in the data space. Therefore, we try to maximize $p_{\theta}(x) = \int_{\mathcal{Z}} p_{\theta}(x | z)p(z)dz$, but it is intractable due to $p_{\theta}(x | z)$. For this reason, we will try to maximize $p_{\theta}(x) = \int_{\mathcal{Z}} q_{\psi}(x | z)p(z)dz$ and to make our approximation $q_{\psi}(x | z)$ as close as possible of $p_{\theta}(x | z)$. This will be accomplished by maximizing the evidence lower bound (ELBO).

$$\begin{aligned} \log p_{\theta}(x) &\ge \log p_{\theta}(x) - D_{\text{KL}}(q_{\phi}(z|x) \| p_{\theta}(z|x)) \\ &= \mathbb{E}_{q_{\phi}(z|x)}[\log p_{\theta}(x | z)] \\ &\quad - D_{\text{KL}}(q_{\phi}(z|x) \| p(z)) \end{aligned}$$

The decoder can be seen as a smooth immersion recovering the Riemannian submanifold upon which live the data (Shao et al., 2018; Arvanitidis et al., 2018; Syrota et al., 2025). Formally, if $\mathcal{M}$ is a $d$-dimensional manifold immersed in $\mathbb{R}^D$ defined with a smooth $(\mathcal{C}^{\infty})$ immersion $\mathfrak{M}: \mathcal{Z} \subseteq \mathbb{R}^d \to \mathcal{X} \subseteq \mathbb{R}^D$ such that $\mathcal{M} = \mathfrak{M}(\mathcal{M})$, the objective of the VAE is to estimate $\mathfrak{M}$ through its decoder. We denote the latent space of the VAE as $\mathcal{Z}$, which is usually called the parametric space of the parametrization $\mathfrak{M}$ of the manifold $\mathcal{M}$.

The smooth immersion $\mathfrak{M}$ induces a Riemannian metric on $\mathcal{M}$, with $g_{\mathcal{M}} = (D\mathfrak{M})^{\top}(D\mathfrak{M})$ being the metric tensor. The Riemannian distance, or inner distance $d_{\mathcal{M}}$, between two points on the manifold $\mathcal{M}$ is the length of the shortest geodesic between the two:

$$d_{\mathcal{M}}(X, Y) = \begin{cases} \inf_{\gamma} \int_{0}^{1} \sqrt{g_{\mathcal{M}}(\dot{\gamma}(t), \dot{\gamma}(t))} dt & \text{s.t. } \gamma \in \mathcal{C}^{1}, \gamma(0) = X, \gamma(1) = Y \\ \inf_{\gamma} \int_{0}^{1} \sqrt{g_{\mathcal{M}}(\dot{\gamma}(t), \dot{\gamma}(t))} dt & \text{s.t. } \gamma \in \mathcal{C}^{\infty}, \gamma(0) = X, \gamma(1) = Y \end{cases}$$

with $\dot{\gamma}(t)$ the tangent vector of the curve $\gamma$ at point $t$.

A MLP decoder in a VAE design defines a smooth immersion if (i) its activation functions are $\mathcal{C}^{\infty}$; and (ii) each weight matrix has maximal rank. In this case, the metric tensor can be derived from the decoder Jacobian (Shao et al., 2018). When the decoder is variational, the expected Riemannian metric of the manifold the stochastic generating process induce is:

$$\mathbb{E}(g_{\mathcal{M}}) = (J^{(\mu)})^{\top} (J^{(\mu)}) + (J^{(\sigma)})^{\top} (J^{(\sigma)})$$

where $J^{(\mu)}$ and $J^{(\sigma)}$ are the Jacobian matrices of the $\mu$ and $\sigma$ heads of the decoder (Arvanitidis et al., 2018). The variance of the metric under the L2 measure vanishes when the data dimension increases (Tosi et al., 2014).

While standard VAEs are known to introduce regularization biases that may distort the learned isometry (Miolane & Holmes, 2020), they nonetheless recover the topological dimension (Zheng et al., 2022) and provide a differentiable parametric approximation of the manifold geometry. We acknowledge that the choice of prior can affect the shape of the Riemannian latent space (Davidson et al., 2018; Kalatzis et al.). However, in our preliminary experiments, we implemented the same algorithm with more expressive priors (notably learnable ones based on Riemannian Brownian motions (Kalatzis et al.)), and observed that performance gains were negligible despite a dramatic increase in computational cost. This suggests that for
the specific task of detecting geometric misalignment, the structural tension is robustly captured even by standard Gaussian priors.

One main issue with this framework comes from the behavior of neural networks: as both of our mean and variance functions are estimated with neural networks, they will behave unexpectedly outside their training region. Following (Arvanitidis et al., 2018), we estimate the inverse variance estimate called precision $\beta_{\psi}(z) = \frac{1}{\sigma_{\psi}^2(z)}$ which is specially implemented with a neural network that extrapolates towards 0 (which gives a guarantee for the variance to be large outside of the training regions). This extrapolation towards 0 is produced by the use of a radial basis function (RBF) neural network (Que & Belkin, 2016) trained afterwards on the outputs of the variational decoder. As this did not bring any gains nor changes to our results, we decided to ignore this step in the results we presented here for the sake of clarity and conciseness.

# B. Implementation details

## B.1. Losses

In phase 1 of the training of our VAE, we use a three-part reconstruction loss that combines a standard reconstruction likelihood with two variance-matching penalties to stabilize training and avoid any posterior collapse. The first term is the Gaussian negative log-likelihood, ensuring that the decoder's mean $\mu$ and variance $\sigma^2$ explain the observed data $x_i$. The second term compares the variance of the decoder's means $\mu$ across dimensions with the empirical variance of the target $x_i$, penalizing deviations so that the deterministic component of the decoder captures the true signal variability rather than collapsing. The third term enforces consistency for the stochastic part of the decoder: the average predicted variance $\mathbb{E}(\sigma^2)$ must match the residual variance in the data not explained by $\mu$. Together, these terms encourage the model to (i) reconstruct the data faithfully, (ii) attribute variance to either structured signal ($\mu$) or noise ($\sigma^2$) in proportion to what is present in the data, and (iii) avoid degenerate solutions where either the mean or variance dominates improperly. As we will use the manifold learned at this phase for the second phase, we need to guarantee good training performance for the algorithm at this stage, as every error at this stage would snowball on the second phase.

$$ \mathcal{L}_{\text{phase 1}} = \sum_i \left[ \underbrace{-\log p_{\theta}(x_i | \mu_i, \sigma_i^2) + \lambda_1 \| \text{Var}[\mu_i] - \text{Var}[x_i] \|^2 + \lambda_2 \| \mathbb{E}[\sigma_i^2] - (\text{Var}[x_i] - \text{Var}[\mu_i]) \|^2}_{\text{Reconstruction loss}} \right] + \beta D_{\text{KL}}(q_{\phi}(z | x_i) \| p(z)). $$

Phase 2 loss is trivially implemented from 2. We computed the heat times using the eigenvalues of the Graph Laplacian by defining a temporal range based on the non-trivial spectral limits. Specifically, we set $t_{\min} = 1/\lambda_{\max}$ and $t_{\max} = 4/\lambda_2$, where $\lambda_2$ and $\lambda_{\max}$ denote the smallest and largest non-zero eigenvalues, respectively. We then extract $k = 15$ log-spaced diffusion scales $T = \{t_i\}_{i=1}^k$ within this interval. This multiscale sampling ensures the geometric alignment captures structural information across the full spectral decomposition, from local neighborhoods to global graph topology.

## B.2. Online geodesic computation during Phase 2 training

Keeping our $\mathcal{M}$ bounded smooth manifold equipped with a Riemannian metric tensor $g$, our goal is to evaluate $g$ efficiently at arbitrary points while maintaining compatibility with automatic differentiation frameworks such as PyTorch.

We define a uniform Cartesian grid over the bounded domain of $\mathcal{Z}$ (the parametric space of $\mathcal{M}$). Let $[a_i, b_i]$ denote the bounds along the $i$-th coordinate axis. The grid points along dimension $i$ are

$$ x_i^k = a_i + k \Delta_i, \quad k = 0, 1, \dots, P, \quad \Delta_i = \frac{b_i - a_i}{P}. $$

The full grid $\mathcal{G}$ is

$$ \mathcal{G} = \{(x_1^{k_1}, \dots, x_d^{k_d}) \mid k_i \in \{0, \dots, P\}\}. $$

Let $\mathcal{G} \in \mathbb{R}^{|\mathcal{G}| \times d \times d}$ store the precomputed metric tensors at all grid points. To approximate geodesic distances on $\mathcal{M}$, we transform the discretized grid $\mathcal{G}$ into a weighted graph. Each grid point $x \in \mathcal{G}$ is represented as a node in the graph, and edges are introduced between neighboring grid points according to the Cartesian connectivity (e.g. 2d-connectivity for axis-aligned neighbors, optionally augmented with diagonal neighbors for improved accuracy).
The weight assigned to an edge $(i, j)$ corresponds to the Riemannian length of the infinitesimal displacement vector $v = x_j - x_i$ under the local metric tensor. Since the exact length is defined by the line integral

$$ \ell_{ij} = \int_{0}^{1} \sqrt{v^{\top} g(x_i + tv) v} dt, $$

we approximate this quantity using quadrature. In practice, we adopt the midpoint rule

$$ w_{ij} = \sqrt{v^{\top} \left( \frac{1}{2}(g_i + g_j) \right) v}, $$

which yields a symmetric weight $w_{ij} = w_{ji}$, ensuring that the resulting graph is undirected.

The outcome is a weighted undirected graph $G = (\mathfrak{G}, E, w)$, where $E$ denotes the set of neighbor connections and $w: E \to \mathbb{R}^+$ assigns edge lengths. Standard shortest-path algorithms (e.g. Dijkstra, A*) can then be applied to compute approximate geodesic distances between any two grid points. For off-grid query points, the metric tensor $g$ is recovered by interpolation from $G$, preserving compatibility with automatic differentiation frameworks such as PyTorch.

For $N$ query points, solving the geodesic equation requires $O(N^2)$ initial value problems to be solved when computing all-pairs distances, with each solve itself involving iterative numerical integration steps whose cost grows with the stiffness of the metric. This quickly becomes prohibitive.

In contrast, once the grid and graph are constructed, geodesic distances reduce to a graph shortest-path problem. For example, Dijkstra's algorithm computes single-source distances in $O(|E| \log |\mathfrak{G}|)$ time, and all-pairs distances can be obtained in $O(|\mathfrak{G}| |E| \log |\mathfrak{G}|)$, which is far more efficient in practice for large grids. Moreover, the graph structure is sparse and can be reused across multiple queries.

Since shortest-path algorithms such as Dijkstra's are inherently discrete and non-differentiable, directly backpropagating through the path computation would break the gradient flow. To address this, we decouple the discrete path selection from the continuous evaluation of metric-dependent edge weights. In practice, we first compute the shortest path without gradients, treating the sequence of edges as fixed with respect to automatic differentiation. Once the path is obtained, we recover the geodesic length by summing the precomputed differentiable edge weights along that path. In this way, the gradients still propagate correctly through the metric tensor evaluations and interpolations stored in $G$, while avoiding the non-differentiability of the path-finding procedure itself.

While the proposed grid-based graph approach offers significant speedups for geodesic computations on low-dimensional manifolds, its scalability is inherently constrained by the dimension $d$ of the parametric space $\mathcal{Z}$. The total number of grid points scales as $|\mathfrak{G}| = (P + 1)^d$, and the memory required to store the precomputed metric tensor $G$ grows as $O(d^2 P^d)$. Consequently, the graph construction and the shortest-path search–with a time complexity of approximately $O(dP^d \log P)$ for sparse connectivity–suffer from the *curse of dimensionality*. This exponential growth limits the applicability of the uniform grid method primarily to low-dimensional manifolds (e.g., $d = 2$ or $d = 3$), as the computational cost and memory footprint become prohibitive for higher-dimensional latent spaces.

To bypass this scalability bottleneck, we introduce an estimator of geodesic distances with a linear interpolation scheme between latent points $u, v \in \mathbb{R}^{d_z}$. We discretize the straight line in latent space,

$$ x_t = (1 - t)u + tv, \quad t \in [0, 1], $$

into $T$ points and evaluate the length of the curve with respect to the induced metric:

**Definition 3 (Linear interpolation distance approximation).** For two points $(u, v)$ on a Riemannian manifold $M$, we define our Linear interpolation distance approximation $(\hat{d}_{LI}(u, v))$ as:

$$ \hat{d}_{LI}(u, v) \triangleq \int_{0}^{1} \sqrt{\dot{\gamma}(t)^{\top} g(x_t) \dot{\gamma}(t)} dt, \quad \dot{\gamma}(t) = v - u. $$

In practice, the integral is approximated with trapezoidal integration over $T$ discretization steps. This reduces the cost to evaluating the metric tensor at $O(T)$ points along the segment, which is tractable and parallelizable. Our implementation computes batched distances in float precision, making it feasible to evaluate approximate manifold distances even for large graphs.
This estimator is designed as a fallback for high-dimensional settings (e.g., $d \gg 3$) where grid-based graph construction is infeasible due to the curse of dimensionality. Despite the theoretical sub-optimality of assuming a straight path in latent space, we empirically found $\hat{d}_{LI}$ to be extremely correlated with the grid-based geodesic estimator in the context of Variational Autoencoder (VAE) latent spaces. We hypothesize that this effectively captures the topology of many learned manifolds, in which the path typically needs to bypass only a single dominant curvature obstacle (e.g., the low-density region near the origin of the prior). Consequently, the straight-line latent interpolation provides a highly efficient yet representative proxy for the true geodesic distance in these regimes.

# C. Datasets

## C.1. Synthetic dataset

### C.1.1. DATASET CONSTRUCTION

**Manifold data** Our synthetic dataset is a collection of points that lie on or near a 2d manifold embedded within a higher-dimensional ambient space ($D = 20$). We sample a collection of $N = 500$ 2d-vectors using coordinates $(u, v) \in \mathbb{R}^2$ and map them through a smooth immersion. Specifically, we generate a soft swiss-roll-like surface without hard boundaries or singularities. The construction defines a smooth radial component $r(u) = r_0 + \text{softplus}(\text{spread} \cdot u)$ and an angular component $\theta(v) = \text{twist} \cdot v$, yielding the embedding

$$f(u, v) = (r(u) \cos(\theta(v)), r(u) \sin(\theta(v)), v).$$

Unlike the standard Swiss roll, this formulation is homeomorphic to $\mathbb{R}^2$ (since $u, v \in \mathbb{R}$) and has full rank everywhere, ensuring that the resulting manifold is smooth and without boundary effects.

To immerse the data in a higher-dimensional ambient space $\mathbb{R}^D$, we extend this base immersion with additional nonlinear features of $(u, v)$. For each extra dimension, we generate a random linear combination $au + bv$ (with $a, b$ drawn from a Gaussian distribution) and apply a smooth nonlinear transformation selected from a family of functions (e.g., sin, cos, tanh, $\exp(-0.1x^2)$, or $x/(1+x^2)$). This guarantees that the higher-dimensional embedding retains smoothness while introducing additional curvature and variability across dimensions. As a result, the final dataset forms a nontrivial 2d manifold in $\mathbb{R}^{20}$, suitable for evaluating generative models and representation learning methods, as this approach ensures that every dimension of $x$ is a continuous function of both $u$ and $v$.

<figure><img src="image_7.png" /><figcaption>Figure 6. Our synthetic data set generated using manifold generation.</figcaption></figure>

**Compatibility matrix and network generation** As we want to work on a graph with realistic properties, we need to find a way to generate from the attributes a network that depends on them, but not in a totally trivial way. We use an existing
generator (Labarthe & Kerzreho, 2026): following a matrix of pairwise compatibility scores, the algorithm returns a network that features similar properties as a social network (power law degree distribution for instance) using an agent-based model. This also fits our requirement of generating a graph that will rely on node attributes, but not too trivially, as the generating process is based on homophily, but the resulting graph is not fully homophilic.

To generate this matrix of pairwise compatibility scores in a way that incorporates the manifold structure, we rely on an exponential similarity matrix: for each pair of vectors $z_a, z_b \in [0, 2\pi]^2$, we take as similarity score $s_{a,b}$:

$$s_{a,b} = \exp \left( - \frac{2 \|z_1 - z_2\|^2}{\underset{i,j \le N}{\operatorname{med}} (\|z_i - z_j\|^2)} \right)$$

To increase sparsity, we remove all values of $s_{a,b}$ that are below 0.2. We illustrate our corruption process in Figure 7.

<figure><img src="image_8.png" /><figcaption>Figure 7. Our perturbation process: the original adjacency matrix, the matrix after placing the 70 selected nodes at the end, and the final poisoned matrix (left to right).</figcaption></figure>

### C.1.2. ARCHITECTURE

For this experiment, the encoder is implemented as a feedforward MLP with two hidden layers and an ELU activation function. No dropout regularization is applied to the decoder, as we overfit the manifold we are trying to learn (we use dropout = 0.2 on the encoder). The node attribute decoder is variational, with two hidden layers, also using ELU activations. Table 2 summarizes the layer-wise architecture.

<table><caption>Table 2. Layers configurations of the encoder and decoder.</caption><thead><tr><th>Module</th><th>Layer type</th><th>Dimensions</th></tr></thead><tbody><tr><td rowspan="4">Encoder</td><td>Input</td><td>d<sub>in</sub> = 20</td></tr><tr><td>Dense</td><td>d<sub>in</sub> → 16</td></tr><tr><td>Dense</td><td>16 → 16</td></tr><tr><td>Latent projection</td><td>16 → 2</td></tr><tr><td rowspan="4">Decoder</td><td>Input</td><td>z ∈ ℝ<sup>2</sup></td></tr><tr><td>Dense</td><td>2 → 16</td></tr><tr><td>Output (mean head)</td><td>16 → d<sub>in</sub></td></tr><tr><td>Output (variance head)</td><td>16 → d<sub>in</sub></td></tr></tbody></table>

To encourage a gradual regularization of the latent distribution, we employ a sigmoid-based KL annealing schedule. The KL weight increases linearly from 0.0 to 2 over 1500 steps (our epochs number for the Phase 1 training). This prevents premature collapse of the latent space. The learning rate during Phase 1 is fixed at $5 \times 10^{-3}$, optimized with Adam. No learning rate decay or warmup is used, ensuring stable convergence throughout the long training horizon.
## C.2. Empirical dataset

### C.2.1. DATASET CONSTRUCTION

**Data** We integrate two complementary sources of information to construct our mobility-socioeconomic graph over the Île-de-France region in France. First, we use the *Comparateur de territoires* dataset published by the French National Institute of Statistics and Economic Studies (INSEE), which provides city-level demographic and socioeconomic indicators (population, housing, employment, enterprises, median income, etc.) together with geospatial geometries of commune (city) boundaries. Second, we exploit the *General Transit Feed Specification* (GTFS) dataset published by Île-de-France Mobilités (IDFM), which encodes the region's public transport offer in terms of stops, stop times, trips, and routes. The GTFS format allows us to reconstruct the operational connectivity of the public transit network (metro, RER, tramway, bus, etc.) with precise temporal information.

**Preprocessing** From the GTFS feed, we extract stops (with geographic coordinates), trips (vehicle journeys following a sequence of stops), and routes (the transit lines grouping trips of the same service). We ensure consistency by enforcing a one-to-one mapping between each trip and its parent route, and enrich every stop-to-stop movement with its associated vehicle category (route type). Each route type is mapped to a canonical transit mode (e.g., bus, metro, tramway, regional train), and we assign an approximate passenger capacity to capture heterogeneous vehicle sizes (e.g., 60 for buses, 1, 600 for regional trains). Temporal information is standardized by correcting for wraparound at midnight and applying a minimum threshold to prevent unrealistically short travel times, ensuring robust, interpretable inter-stop travel times.

**Stop-level network construction** For each trip, we build ordered stop pairs by matching consecutive stops along its sequence. Each directed stop-stop edge is assigned an effective contribution score defined as

$$ \text{contrib}(u, v) = \frac{\text{capacity(trip)}}{\tau + \text{transport_time}(u, v)} $$

where $\tau$ is a smoothing hyperparameter preventing division by zero while controlling the relative weight of capacity versus travel time. Edges are aggregated across trips to yield undirected stop-stop links, with attributes such as the number of trips, total capacity, and average travel time. The contribution score provides a monotonic measure of effective connectivity that favors high-capacity and short-travel-time links. This contribution design make the addition of bus lines or inefficient lines (in travel times) always an improvement for the connectivity score between two stops, even if it increases the average transport time. Indeed, we consider that a higher number of transportation possibilities increases the convenience of transportation.

**City-level graph aggregation** We then embed the stop network into the INSEE geospatial layer. Each GTFS stop is converted into a point geometry and spatially joined to the commune polygon in which it lies. This yields a stop-to-commune mapping. For each pair of communes, we sum contribution scores across all inter-commune stop pairs, thus producing weighted edges reflecting the potential mobility supply between municipalities. Additional attributes (trips, capacity, average travel time) are stored as auxiliary edge features.

**City-level Socioeconomic indicators** From the INSEE Comparateur des territoires, we retain a reduced set of 20 commune-level descriptors covering population, households, housing, income, employment, and enterprises. For demographics, we include population counts in 2020 (p20_pop) and 2014 (p14_pop), together with the commune's surface area (superf), household counts (p20_men), domiciled births (naisd22) and deaths (decesd22). Housing is represented by the number of dwellings in 2020 (p20_log), broken down into main residences (p20_rp), secondary or occasional residences (p20_rsecocc), and vacant dwellings (p20_logvac). Economic well-being is captured by the median standard of living in 2020 (med20). Employment indicators include the total number of jobs at the workplace (p20 Emplt), salaried jobs specifically (p20 Emplt sal), and their historical value in 2014 (p14 Emplt). Finally, productive structure is characterized by the number of establishments in 2021 (ettot21), further disaggregated into agriculture, forestry and fishing (etaz21), industry (etbe21), construction (etfz21), and commerce/transport/services (etgu21), along with the number of medium-to-large establishments employing ten or more workers (ettefp1021). Many of these descriptors are redundant-for instance, total population across years, or dwellings and their subcategories-so the feature space is highly correlated. **This redundancy justifies the manifold hypothesis,** namely that socioeconomic profiles of communes lie on a lower-dimensional smooth structure embedded in the high-dimensional indicator space.
**Node labelling** Commune-level socioeconomic attributes are preprocessed by quantile transformation (with Gaussian output distribution), thereby normalizing heterogeneous scales and mitigating the effect of heavy-tailed distributions (e.g., income, population, enterprise counts). Nodes inherit their INSEE socioeconomic descriptors, standardized as described above. To ensure robustness, communes with missing socioeconomic attributes are removed, as are isolated communes with no incident edges.

<figure><img src="image_9.png" /><figcaption>Figure 8. Final Transportation network extracted from GTFS data. The darker and thicker the link is, the easier it is to go from one city to another.</figcaption></figure>

### C.2.2. ARCHITECTURE

For this experiment, the encoder is implemented as a feedforward MLP with two hidden layers and ReLU activation functions. No dropout regularization is applied to the decoder as we overfit the manifold we try to learn (we use dropout = 0.2 on the encoder). The node attribute decoder is variational, with two hidden layers, also using ReLU activations. Table 3 summarizes the layer-wise architecture.

<table><caption>Table 3. Layers configurations of the encoder and decoder.</caption><thead><tr><th>Module</th><th>Layer type</th><th>Dimensions</th></tr></thead><tbody><tr><td rowspan="4">Encoder</td><td>Input</td><td>d<sub>in</sub> = 20</td></tr><tr><td>Dense</td><td>d<sub>in</sub> → 96</td></tr><tr><td>Dense</td><td>96 → 64</td></tr><tr><td>Latent projection</td><td>64 → 2</td></tr><tr><td rowspan="4">Decoder</td><td>Input</td><td>z ∈ ℝ<sup>2</sup></td></tr><tr><td>Dense</td><td>2 → 128</td></tr><tr><td>Output (mean head)</td><td>128 → d<sub>in</sub></td></tr><tr><td>Output (variance head)</td><td>128 → d<sub>in</sub></td></tr></tbody></table>

To encourage a gradual regularization of the latent distribution, we employ a linear KL annealing schedule. The KL weight increases linearly from 0.0 to 1 over 1200 steps (our epochs number for the phase 1 training). This prevents premature collapse of the latent space. The learning rate during Phase 1 is fixed at $5 \times 10^{-3}$, optimized with Adam. No learning rate decay or warmup is used, ensuring stable convergence throughout the long training horizon.

### C.2.3. CITY CLUSTERS
<table><caption>Table 4. Cluster definitions derived from surrogate decision tree rules after VAE Phase 1 training. Feature thresholds are reported as standard normal quantiles ($\sigma$). Representative Communes listed are the medoids (most central instances) of each cluster in the learned latent space.</caption><thead><tr><th>ID</th><th>Key Distinguishing Rules</th><th>Socioeconomic Profile</th><th>Representative Medoids</th></tr></thead><tbody><tr><td>3</td><td>Pop'14 &gt; +2.91σ</td><td>Metropolitan Centers<br>Extreme population density.</td><td>Paris 20<sup>e</sup>, Paris 18<sup>e</sup>, Paris 8<sup>e</sup></td></tr><tr><td>2</td><td>Housing'20 ∈ (0.82, 2.91]<sup>σ</sup><br>Constraint: Pop'14 ≤ 2.91</td><td>Dense Urban<br>High residential density without being a global outlier.</td><td>Le Plessis-Robinson, Houilles, Chatou, Brétigny-sur-Orge, Limeil-Brévannes</td></tr><tr><td>4</td><td>Large Est. &gt; -2.92σ<br>OR Industry &gt; -2.79σ</td><td>Productive Nodes<br>Defined by specific industrial presence or large establishments.</td><td>Monthyon, Thoiry, Bennecourt, Maincy, Janville-sur-Juine</td></tr><tr><td>1</td><td>Pop'20 &gt; 0.15σ<br>Constraint: Low Ind./Agri.</td><td>Residential Periphery<br>Above-average population but low economic activity.</td><td>Lisses, Montlhéry, Le Pecq, Fleury-Mérogis, Aubergenville</td></tr><tr><td>0</td><td>All other communes<br>(Fails all criteria above)</td><td>Rural & Low-Activity<br>The baseline majority; low density and establishment counts.</td><td>Everly, Marchémoret, Coulommes, Boutervilliers, Fouju</td></tr></tbody></table>