Bures-Wasserstein Flow Matching for Graph Generation

## Abstract

Abstract Graph generation has emerged as a critical task in fields ranging from molecule design to drug discovery. Contemporary approaches, notably diffusion and flow-based models, have achieved solid graph generative performance through constructing a probability path that interpolates between a reference distribution and the data distribution. However, these methods typically model the evolution of individual nodes and edges independently and use linear interpolations to build the path assuming that the data lie in Euclidean space. We show that this is suboptimal given the intrinsic non-Euclidean structure and interconnected patterns of graphs, and it poses risks to the sampling convergence. To build a better probability path, we model the joint evolution of the nodes and edges by representing graphs as connected systems parameterized by Markov random fields (MRF). We then leverage the optimal transport displacement between MRF objects to design the probability path for graph generation. Based on this, we introduce BWFlow, a flow-matching framework for graph generation that respects the underlying geometry of graphs and provides smooth velocities in the probability path. The novel framework can be adapted to both continuous and discrete flow-matching algorithms. Experimental evaluations in plain graph generation and 2D/3D molecule generation validate the effectiveness of BWFlow in graph generation with competitive performance, stable training, and guaranteed sampling convergence.

## 1 Introduction

Thanks to the capability of graphs in representing complex relationships, graph generation (Zhu et al., 2022; Liu et al., 2023a) has become an essential task in various fields such as protein design (Ingraham et al., 2019), drug discovery (Bilodeau et al., 2022), and social network analysis (Li et al., 2023). Among contemporary generative models, diffusion and flow models have emerged as two compelling approaches for their ability to achieve state-of-the-art performance in graph generation (Niu et al., 2020; Vignac et al., 2023a; Eijkelboom et al., 2024; Qin et al., 2024; Hou et al., 2024). In particular, these generative models can be unified under the framework of stochastic interpolation (Albergo & Vanden-Eijnden, 2023), which consists of four procedures (Lipman et al., 2024): 1) Drawing samples from the reference (source) distribution $p_{0}(\cdot)$ and/or the data (target) distribution $p_{1}(\cdot)$ for training set assembly; 2) Constructing a time-continuous probability path $p_{t}(\cdot),0\leq t\leq 1$ interpolating between $p_{0}$ and $p_{1}$; 3) Training a model to reconstruct the probability path by either approximating the score function or velocity fields (ratio matrix in the discrete case); and 4) sampling from $p_{0}$ and transforming it through the learned probability path to get samples that approximately follow $p_{1}$.

A core challenge in this framework is constructing the probability path $p_{t}$. Existing text and image generative models, operating either in the continuous (Ho et al., 2020; Song et al., 2021; Lipman et al., 2023; Liu et al., 2023b) or discrete (Campbell et al., 2022; Sun et al., 2023; Campbell et al., 2024; Gat et al., 2024; Minello et al., 2025) space, typically rely on linear interpolation between source and target distributions to construct the path. Graph generation models, including diffusion (Niu et al., 2020; Vignac et al., 2023a; Haefeli et al., 2022; Xu et al., 2024; Siraudin et al., 2024) and flow-based models (Eijkelboom et al., 2024; Qin et al., 2024; Hou et al., 2024), inherit this design by modeling every single node and edge independently and linearly build paths in the disjoint space. However, this approach is inefficient because it neglects the strong interactions and relational structure inherent in graphs, i.e., the significance of a node heavily depends on the configuration of its neighbors. While empirical success have been achieved via fine-grained searching on the training and sampling design (Qin et al., 2024) such as target guidance and time distortion, we argue that there remains a fundamental issue of the linear probability path construction, and these strategies only mitigate the problem by manipulating the probability path.

Figure: (a) Training path comparison
Refer to caption: /html/2506.14020/assets/x1.png

Motivation examples. The blue line in [Figure˜1(a)](#S1.F1.sf1) illustrates the probability path evolution through linear interpolation in plain graph generation, where the probability path remains flat until $t\approx 0.8$ before sharply dropping. This pattern provides a poor velocity(^1^11Formally, the per-sample velocity is the gradient of the corresponding sample-specific path. Here, we average among samples and view the path tangent as an approximation for the marginal velocity. [Section 2](#S2) gives a more accurate definition.) estimation, as ideally the velocity field should be smooth and consistently pointing to the data distribution as the green line in [Figure˜1(a)](#S1.F1.sf1)  (Kapusniak et al., 2024). The resultant velocity will make the sampling difficult to converge to the data distribution, as shown in [Figure˜1(c)](#S1.F1.sf3), where the termination point of the blue curve remains high. Though not explicitly mentioned, the superior performance achieved by Qin et al. (2024) is partially attributed to an extensive design space search to manipulate the path for smoother velocity estimation. The techniques they used, including target guidance, time distortion, and stochasticity injection, are conceptually visualized in [Figure˜1(a)](#S1.F1.sf1) and [1(b)](#S1.F1.sf2) with discussions in [Section˜F.1](#A6.SS1). More motivating examples in [Section˜I.2](#A9.SS2).

The limitations. The above examples reveal fundamental issues of the probability path construction in graph generation, attributable to two primary reasons:
1) The assumption of independence between nodes/edges and a linear interpolation in the locally disjoint space fails to capture the global co-evolution of graph components and properties such as community structure or spectrum (Haasler & Frossard, 2024). This potentially causes a sudden transition from reference to data distribution, which yields non-smooth probability paths. 2) The linear interpolation is derived through the optimal transport (OT) displacement (Tong et al., 2024) between distributions residing in the Euclidean space (Lipman et al., 2024). However, linear interpolation would stray away from the true data manifold when the underlying space is non-Euclidean (Chen & Lipman, 2024; Kapusniak et al., 2024). Since graphs naturally inhabit non-Euclidean geometries, linearly interpolating the nodes/edges neither guarantees an OT displacement nor respects the underlying geometry, making the constructed probability path suboptimal or even deviate from the valid graph domain.

Proposed solution. To address these limitations, we draw on statistical relational learning and model graphs using Markov Random Fields (MRFs) (Taskar et al., 2007; Qu et al., 2019). MRFs organize the nodes/edges as an interconnected system and interpolating between two MRFs captures the joint evolution of the whole graph system. Extending (Haasler & Frossard, 2024), we derive a closed-form Wasserstein distance between graph distributions and leverage it to construct the Bures-Wasserstein (BW) interpolation of two graphs that ensures the OT displacement in graph generation compared to linear interpolation. We then integrate these insights into a flow-matching framework(^2^22For conciseness, we focus on flow models, but the idea is generalizable to diffusions as in [Section F.2](#A6.SS2).) called Bures–Wasserstein Flow (BWFlow). Specifically, by defining a probability path via BW interpolation, we obtain smooth, globally coherent velocity fields at intermediate steps (see [Figure˜1(c)](#S1.F1.sf3)) that respect the non-Euclidean, interconnected structure of graphs. Crucially, BWFlow admits simulation-free computation of densities and velocities along the entire path, which translates into efficient, stable training and sampling.

Contributions. First, we theoretically and empirically show that the linear interpolation used in existing graph generation models gives suboptimal probability path construction and velocity estimation. Second, through parameterizing graphs as MRFs, we introduce BWFlow, a flow-matching model for graph generation that constructs probability paths respecting the graph geometry and develops smooth velocities. Third, BWFlow was tested on plain graph and 2D/3D molecule generation, demonstrating competitive performance without an excessive search for path manipulation techniques. We further show that BW interpolation consistently outperforms other interpolation methods in building flow matching models, leading to more stable training and sampling convergence.

## 2 Preliminaries

### 2.1 Flow matching for graph generation

#### Flow matching (FM).

Generative modeling considers fitting a mapping from state space ${\mathcal{S}}\rightarrow{\mathcal{S}}$ that transforms the samples from source distribution, $X_{0}\sim p_{0}$, to samples from target data distribution, $X_{1}\sim p_{1}$(^3^33For clarity, we denote the calligraphic style ${\mathcal{X}}$ being the random variable, the plain $X$ the relevant realizations and the bold symbol ${\bm{X}}$ the latent variables (parameters) that controls the distributons, i.e. $X\sim p({\mathcal{X}};{\bm{X}})$.). Continuous normalizing flow (CNF) (Chen et al., 2018) parameterizes the transformation through a push-forward equation that interpolates between $p_{0}$ and $p_{1}$ and constructs a probability path $p_{t}({\mathcal{X}})=\left[\psi_{t}p_{0}\right]({\mathcal{X}})$ through a time-dependent function $\psi_{t}$ (a.k.a flow). A vector field $u_{t}$, defined as $\frac{d}{dt}\psi_{t}\left({\mathcal{X}}\right)=u_{t}\left(\psi_{t}\left({\mathcal{X}}\right)\right)$ with $\psi_{0}\left({\mathcal{X}}\right)={\mathcal{X}}$, is said to generate $p_{t}$ if $\psi_{t}$ satisfies $X_{t}:=\psi_{t}\left(X_{0}\right)\sim p_{t}\text{ for }X_{0}\sim p_{0}$. The FM (Lipman et al., 2023) is designed to match the real velocity field through the loss:

$$ ${\mathcal{L}}_{\mathrm{FM}}(\theta)=\mathbb{E}_{t,X_{t}\sim p_{t}(\cdot)}\left\|v_{\theta}(X_{t})-u_{t}(X_{t})\right\|^{2}.$ (1) $$

where $v_{\theta}(\cdot):{\mathcal{S}}\rightarrow{\mathcal{S}}$ is the parameterized velocity field and $t\sim\mathcal{U}[0,1]$.

#### Conditional flow matching (CFM).

Given that the actual velocity field and the path are not tractable (Tong et al., 2024), one can construct the per-sample conditional flow. We condition the probability paths on variable $Z\sim\pi(\cdot)$ (for instance, a pair of source and target points $Z=(X_{0},X_{1})$) and re-write $p_{t}({\mathcal{X}})={\mathbb{E}}_{\pi(\cdot)}p_{t}({\mathcal{X}}\mid Z)$ and $u_{t}({\mathcal{X}})={\mathbb{E}}_{\pi(\cdot)}u_{t}({\mathcal{X}}\mid Z)$ where the conditional path and the velocity field are tractable. The CFM aims at regressing a velocity $v_{\theta}(\cdot)$ to $u_{t}({\mathcal{X}}\mid Z)$ by the loss,

$$ $\mathcal{L}_{\mathrm{CFM}}(\theta):=\mathbb{E}_{t,Z\sim\pi(\cdot),p_{t}(\cdot\mid Z)}\left\|v_{\theta}(X_{t})-u_{t}(X_{t}\mid Z)\right\|^{2},$ (2) $$

where it is shown that the CFM optimization has the same optimum as the FM objective (Tong et al., 2024).

#### Graphs as statistical objects.

When considering graph generation with CFM, the very first step is to model graphs as statistical objects. For notation, we let ${\mathcal{G}}=\{{\mathcal{V}},{\mathcal{E}},{\mathcal{X}}\}$ denote an undirected graph random variable with edges ${\mathcal{E}}=\{e_{uv}\}$, nodes ${\mathcal{V}}=\{v\}$, and node features ${\mathcal{X}}=\{x_{v}\}$. A graph realization is denoted as $G=\{V,E,X\}~\sim p({\mathcal{G}})$.
We consider a group of latent variables that controls the graph distribution, specifically the node feature mean ${\bm{X}}=\left[{\bm{x}}_{1},{\bm{x}}_{2},\ldots,{\bm{x}}_{|{\mathcal{V}}|}\right]^{\top}\in\mathbb{R}^{{|{\mathcal{V}}|}\times K}$, the weighted adjacency matrix ${\bm{W}}\in\mathbb{R}^{|{\mathcal{V}}|\times|{\mathcal{V}}|}$, and the Laplacian matrix ${\bm{L}}={\bm{D}}-{\bm{W}}\in\mathbb{R}^{|{\mathcal{V}}|\times|{\mathcal{V}}|}$, with ${\bm{D}}=\operatorname{diag}({\bm{W}}\mathbf{1})$ being the degree matrix (and $\mathbf{1}$ the all-one vector). In a nutshell, graphs are sampled from $G\sim p({\mathcal{G}};{\bm{G}})=p({\mathcal{X}},{\mathcal{E}};{\bm{X}},{\bm{W}})$.

#### Graph generation with CFM.

The CFM samples new graphs through iteratively building $G_{t+dt}=G_{t}+v^{\theta}_{t}(G_{t})\cdot dt$ with initial $G_{0}\sim p_{0}({\mathcal{G}})$ and a trained velocity field $v_{t}^{\theta}(G_{t})$, so that the medium points follows $G_{t}\sim p_{t}({\mathcal{G}})$ and terminates at $p_{1}({\mathcal{G}})$. We can parameterize $v^{\theta}_{t}(G_{t})$ as in (Gat et al., 2024),

$$ $v^{\theta}_{t}(G_{t})=\mathbb{E}_{G_{0}\sim p_{0}({\mathcal{G}}),G_{1}\sim p^{\theta}_{1\mid t}\left(\cdot\mid G_{t}\right)}\left[v_{t}\left(G_{t}\mid G_{0},G_{1}\right)\right]$ (3) $$

As such, training the velocity fields is replaced by training a denoiser $p^{\theta}_{1\mid t}\left(\cdot\mid G_{t}\right)$ to predict the clean datapoint, which is equivalent to maximizing the log-likelihood (Qin et al., 2024; Campbell et al., 2024),

$$ $\mathcal{L}_{\text{CFM}}=\mathbb{E}_{G_{1}\sim p_{1},G_{0}\sim p_{0},t\sim\mathcal{U}_{[0,1]},G_{t}\sim p_{t\mid 0,1}}\left[\log p_{1\mid t}^{\theta}\left(G_{1}\mid G_{t}\right)\right]$ (4) $$

where $t$ is sampled from a uniform distribution $\mathcal{U}_{[0,1]}$ and $G_{t}\sim p_{t\mid 0,1}$ can be obtained in a simulation-free manner. This framework avoids the evaluation of the conditional vector field at training time, which both increases the model robustness and training efficiency.

To proceed, a closed form of $p_{t}(\cdot\mid G_{0},G_{1})$ is required to construct both the probability path and the velocity field $v_{t}\left(G_{t}\mid G_{0},G_{1}\right)$. A common selection to decompose the probability density assumes independency for each node and edge (Hou et al., 2024; Qin et al., 2024; Eijkelboom et al., 2024) giving $p({\mathcal{G}})=p({\mathcal{X}})p({\mathcal{E}})=\prod_{v\in{\mathcal{V}}}p(x_{v})\prod_{e_{uv}\in{\mathcal{E}}}p(e_{uv})$. Choosing $\pi(\cdot)=p_{0}\left({\mathcal{G}}\right)p_{1}\left({\mathcal{G}}\right)$, the boundary conditions follow $p_{i}({\mathcal{G}})=\delta({\mathcal{X}}_{i}={\bm{X}}_{i})\cdot\delta({\mathcal{E}}_{i}={\bm{W}}_{i}),\forall i=\{0,1\}$ with $\delta$ the dirac function. This decomposition is further combined with linear interpolation to build the path, as introduced in (Tong et al., 2024), where,

$$ $\displaystyle p_{t}(x_{v}\mid G_{0},G_{1})=\mathcal{N}\left(t[X_{1}]_{v}+(1-t)[X_{0}]_{v},\sigma_{t}^{2}\right)$ (5) $\displaystyle p_{t}(e_{uv}\mid G_{0},G_{1})=\mathcal{N}\left(t[E_{1}]_{uv}+(1-t)[E_{0}]_{uv},\sigma_{t}^{2}\right).$ $$

Similarly, discrete flow matching frameworks for graph generation (Qin et al., 2024; Siraudin et al., 2024; Xu et al., 2024) is also based on linear interpolation, where the interpolant is sampled from a categorical distribution whose probabilities are simply linear interpolation between the boundary conditions.

### 2.2 Optimal transport and flow matching

Why linear interpolation? Existing literature (Liu et al., 2023b; Albergo & Vanden-Eijnden, 2023) argues that the probability path $p_{t}({\mathcal{X}}\mid Z)$ should be chosen to recover the optimal transport (OT) displacement interpolant (McCann, 1997). The (Kantorovich) optimal transport problem is to find the transport plan between two probability measures, $\eta_{0}$ and $\eta_{1}$, with the smallest associated transportation cost defined as follows.

###### Definition 1 (Wasserstein Distance) .

When the data follow Euclidean geometry and both boundary distributions, $p_{0}$ and $p_{1}$ are described by the Gaussian family, the probability path shown in [Equation˜5](#S2.E5) with $\sigma_{t}\rightarrow 0$ becomes a solution to [Equation˜6](#S2.E6).

###### Definition 1 (Wasserstein Distance) .

## 3 Methodology

In this paper, we introduce Bures–Wasserstein Flow Matching (BWFlow), a novel graph generation framework that is built upon the OT displacement when modeling graphs with Markov Random Fields (MRFs). We begin by casting graphs in an MRF formulation in [Section˜3.1](#S3.SS1). We then derive the BWFlow framework in [Section˜3.2](#S3.SS2) by formulating and solving the OT displacement problem on the MRF, thereby yielding the fundamental components, interpolations and velocity fields, for FM-based graph generation. Finally, in [Section˜3.3](#S3.SS3), we extend BWFlow to discrete FM regimes, enabling its application across a broad spectrum of graph-generation tasks. A schematic overview of the entire BWFlow is illustrated in [Figure˜2](#S3.F2).

### 3.1 Graph Markov random fields

We borrow the idea from MRF as a remedy to modeling the complex system organized by graphs, which intrinsically captures the underlying mechanism that jointly generates the nodes and edges. Mathematically, we assume the joint probability density distribution (PDF) of node features and graph structure as $p({\mathcal{G}};{\bm{G}})=p({\mathcal{X}},{\mathcal{E}};{\bm{X}},{\bm{W}})=p({\mathcal{X}};{\bm{X}},{\bm{W}})p({\mathcal{E}};{\bm{W}})$ where the node features and graph structure are interconnected through latent variables ${\bm{X}}$ and ${\bm{W}}$. For node features ${\mathcal{X}}$, we follow the MRF assumption in Zhu et al. (2003) and decompose the density into the node-wise potential $\varphi_{1}(v),\forall v\in{\mathcal{V}}$ and pair-wise potential $\varphi_{2}(u,v)$, $\forall e_{uv}\in{\mathcal{E}}$:

$$ $\displaystyle p({\mathcal{X}};{\bm{X}},{\bm{W}})\propto$ $\displaystyle\prod_{v}\underbrace{\exp\left\{-(\nu+d_{v})\|{\bm{V}}x_{v}-\bm{\mu}_{v}\|^{2}\right\}}_{\varphi_{1}(v)}$ (7) $\displaystyle\prod_{u,v}\underbrace{\exp\left\{w_{uv}\left[({\bm{V}}x_{u}-\bm{\mu}_{u})^{\top}({\bm{V}}x_{v}-\bm{\mu}_{v})\right]\right\}}_{\varphi_{2}(u,v)},$ $$

with $\|\cdot\|$ the $L_{2}$ norm, $(\cdot)^{\dagger}$ the pseudo-inverse, ${\bm{V}}$ the transformation matrix modulating the graph feature emission, and $\bm{\mu}_{v}$ the node-specific latent variable mean.  [Equation˜7](#S3.E7) can be expressed as a colored Gaussian distribution in [Equation˜8](#S3.E8) given that ${\bm{V}}x_{v}\sim{\mathcal{N}}(\bm{\mu}_{v},(\nu{\bm{I}}+{\bm{L}})^{-1})$. We further assume that edges are emitted via a Dirac delta, ${\mathcal{E}}\sim\delta({\bm{W}})$, yielding our definition of Graph Markov Random Fields (GraphMRF). The derivation can be found in [Section˜A.2](#A1.SS2).

###### Definition 2 (Graph Markov Random Fields) .

###### Definition 2 (Graph Markov Random Fields) .

### 3.2 Bures-Wasserstein flow matching for graph generation

Figure: Figure 2: Schematic overview of BWFlow, which consists of: a) Sample the marginal graph condition $G_{0}$ and $G_{1}$; b) Convert graphs to MRFs; c) Interpolate to get intermediate points; d) Convert back to get $G_{t}$; e) Train velocity based on $G_{t}$; and f) Generate new points with the trained velocity.
Refer to caption: /html/2506.14020/assets/x3.png

#### The optimal transport displacement between graph distributions.

Given that the joint probability of graphs decomposed as $p({\mathcal{G}})=p({\mathcal{X}};{\bm{X}},{\bm{W}})p({\mathcal{E}};{\bm{W}})$ and the measure factorized to $\eta_{{\mathcal{G}}_{j}}=\eta_{{\mathcal{X}}_{j}}\cdot\eta_{{\mathcal{E}}_{j}}$ with $j\in\{0,1\}$, the graph Wasserstein distance between $\eta_{{\mathcal{G}}_{0}}$ and $\eta_{{\mathcal{G}}_{1}}$ is written as,

$$ $d_{\text{BW}}({\mathcal{G}}_{0},{\mathcal{G}}_{1}):={\mathcal{W}}_{c}(\eta_{{\mathcal{X}}_{0}},\eta_{{\mathcal{X}}_{1}})+{\mathcal{W}}_{c}(\eta_{{\mathcal{E}}_{0}},\eta_{{\mathcal{E}}_{1}}).$ $$

We extend Haasler & Frossard (2024) and analytically derive the graph Wasserstein distance using the OT formula between Gaussians (Dowson & Landau, 1982; Olkin & Pukelsheim, 1982; Takatsu, 2010) (see [Lemma˜2](#Thmlemma2) proved in [Section˜B.1](#A2.SS1)) as follows.

###### Proposition 1 (Bures-Wasserstein Distance) .

Based on the Bures-Wasserstein (BW) distance, we then derive the OT interpolant for two graphs, which is the solution of the displacement minimization problem described as,

$$ ${\mathcal{G}}_{t}=\operatorname*{arg\,min}_{\tilde{{\mathcal{G}}}}~(1-t)d_{\text{BW}}({\mathcal{G}}_{0},\tilde{{\mathcal{G}}})+td_{\text{BW}}(\tilde{{\mathcal{G}}},{\mathcal{G}}_{1}).$ (10) $$

###### Proposition 1 (Bures-Wasserstein Distance) .

#### The probability path.

The interpolation is obtained through solving [Equation˜10](#S3.E10) with the BW distance defined in [˜1](#Thmprop1), we prove the minimizer of the above problem has the form in [˜2](#Thmprop2). The proof can be found in [Section˜C.1](#A3.SS1).

###### Proposition 2 (Bures-Wasserstein interpolation) .

The interpolant provides a closed form for the induced probability path $p({\mathcal{G}}_{t}\mid G_{0},G_{1})$ and the velocity $v(G_{t}\mid G_{0},G_{1})$ that is easy to access without any simulation.

###### Proposition 2 (Bures-Wasserstein interpolation) .

#### The velocity.

We consider the reparameterization as in [Equation˜3](#S2.E3) and derive the conditional velocity $v_{t}\left(G_{t}\mid G_{1},G_{0}\right)$ as in [˜3](#Thmprop3).

###### Proposition 3 (Bures-Wasserstein velocity) .

With [˜2](#Thmprop2) and [˜3](#Thmprop3), we are now able to formally construct the algorithms for Bures-Wasserstein flow matching. Taking continuous flow matching as an example,  [Algorithms˜1](#alg1) and [2](#alg2) respectively introduce the training and sampling pipelines for our BWFlow.

###### Proposition 3 (Bures-Wasserstein velocity) .

### 3.3 Discrete Bures-Wasserstein flow matching for graph generation

Up to now we are working on the scenario when $p({\mathcal{X}}\mid{\bm{X}},{\bm{W}})$ is a Gaussian and $p({\mathcal{E}}\mid{\bm{W}})$ is a Dirac distribution. However, previous studies have observed a significant improvement of the discrete counterpart of the continuous graph generation models (Vignac et al., 2023a; Xu et al., 2024; Qin et al., 2024). To benefit our model from such a nature, we derive the discrete Bures-Wasserstein flow matching for graph generation.

#### The discrete probability path.

We design the probability path as discrete distributions,

$$ $\displaystyle p_{t}(x_{v}\mid G_{0},G_{1})=\operatorname{Categorical}([{\bm{X}}_{t}]_{v}),$ (13) $\displaystyle p_{t}(e_{uv}\mid G_{0},G_{1})=\operatorname{Bernoulli}([{\bm{W}}_{t}]_{uv})$ $\displaystyle\text{ s.t. }p_{0}({\mathcal{G}})=\delta(G_{0},\cdot),p_{1}({\mathcal{G}})=\delta(G_{1},\cdot)$ $$

where ${\bm{W}}_{t}={\bm{D}}_{t}-{\bm{L}}_{t}$ with ${\bm{X}}_{t}$ and ${\bm{L}}_{t}$ defined the same in [Equation˜11](#S3.E11). We consider the fact that the Dirac distribution is a special case when the Categorical/Bernoulli distribution has probability 1 or 0, so the boundary condition $p_{0}({\mathcal{G}})=\delta(G_{0},\cdot),p_{1}({\mathcal{G}})=\delta(G_{1},\cdot)$ holds. Even though we are not sampling from Gaussian distributions anymore, it is possible to approximate the Wasserstein distance between two multivariate discrete distributions with the Gaussian counterpart so the conclusions, such as optimal transport displacements, still hold. More discussions in [Section˜D.2](#A4.SS2).

#### The discrete velocity fields.

The path of node features ${\mathcal{X}}_{t}$ can be re-written as $p_{t}({\mathcal{X}})=(1-t)\delta(\cdot,{\bm{X}}_{0})+t\delta(\cdot,{\bm{X}}_{1})$ so the conditional velocity can be accessed through $v_{t}(X_{t}\mid G_{0},G_{1})=[\delta(\cdot,{\bm{X}}_{1})-\delta(\cdot,{\bm{X}}_{t})]/(1-t)$. However, the probability path of edges ${\mathcal{E}}_{t}$, shown in [Equations˜11](#S3.E11) and [13](#S3.E13), cannot be written as a mixture of two boundary conditions given the non-linear interpolation. To this end, we derive in [Section˜D.3](#A4.SS3) that the discrete velocity follows,

$$ $v_{t}(E_{t}\mid G_{1},G_{0})=\left(1-2E_{t}\right)\frac{\dot{{\bm{W}}_{t}}}{{\bm{W}}_{t}\circ\left(1-{\bm{W}}_{t}\right)},$ (14) $$

where ${\bm{W}}_{t}={\bm{D}}_{t}-{\bm{L}}_{t}$, $\dot{{\bm{W}}_{t}}=\text{diag}(\dot{{\bm{L}}_{t}})-\dot{{\bm{L}}_{t}}$ with ${\bm{L}}_{t}$, $\dot{{\bm{L}}_{t}}$ defined in [Equations˜11](#S3.E11) and LABEL:eq:bw_velocity respectively. With the interpolation and velocity defined, the discrete flow matching is built in [Algorithms˜3](#alg3) and [4](#alg4).

## 4 Experiments

We evaluate the BWFlow algorithms through both the plain graph generation and real-world molecule generation tasks. We first outline the experimental setup in [Section˜4.1](#S4.SS1), followed by a general comparison in [Section˜4.2](#S4.SS2). Next, we conduct behavior analysis on the impact of interpolation methods and the corresponding velocity construction on graph generation performance in [Section˜4.3](#S4.SS3), which demonstrates the effectiveness and benefit of flow along Bures-Wasserstein interpolation.

Figure: (a) The evolution of graph statistics ratio along the probability path.
Refer to caption: /html/2506.14020/assets/figs/png_output/flow_step_ratio.png

**Table 1: Plain Graph generation performance. We sampled 5 times (each run generates 40 graphs) to calculate the mean and standard deviation. We only keep the main diffusion/flow model for comparison, while other models are included in the full version at [Table˜8](#A9.T8).**
|  |  | Planar | Tree | SBM |  |  |  |
| --- | --- | --- | --- | --- | --- | --- | --- |
| Model | Class | V.U.N. $\uparrow$ | A.Ratio $\downarrow$ | V.U.N. $\uparrow$ | A.Ratio $\downarrow$ | V.U.N. $\uparrow$ | A.Ratio $\downarrow$ |
| Train set | — | 100 | 1.0 | 100 | 1.0 | 85.9 | 1.0 |
| DiGress (Vignac et al., 2023a) | Diffusion | 77.5 | 5.1 | 90.0 | 1.6 | 60.0 | 1.7 |
| HSpectre (Bergmeister et al., 2024) | Diffusion | 95.0 | 2.1 | 100.0 | 4.0 | 75.0 | 10.5 |
| GruM (Jo et al., 2024) | Diffusion | 90.0 | 1.8 | — | — | 85.0 | 1.1 |
| CatFlow (Eijkelboom et al., 2024) | Flow | 80.0 | — | — | — | 85.0 | — |
| DeFoG (Qin et al., 2024) | Flow | 99.5±1.0 | 1.6±0.4 | 96.5±2.6 | 1.6±0.4 | 90.0±5.1 | 4.9±1.3 |
| BWFlow | Flow | 97.5±2.5 | 1.3±0.4 | 75.5±2.4 | 1.4±0.3 | 90.5±4.0 | 3.8±0.9 |
| DeFoG (EMA) | Flow | 77.5±8.37 | 3.5±1.7 | 73.1±11.4 | 1.50±0.3 | 85.0±7.1 | 3.7±0.9 |
| BWFlow (EMA) | Flow | 84.8±6.44 | 2.4±0.9 | 68.1±12.2 | 1.24±0.2 | 83.5±6.0 | 2.4±0.6 |

2D molecule graph generation. The model performance is illustrated in [Table˜2](#S4.T2). In both datasets, BWFlow can achieve competitive results near the state-of-the-art (SOTA) flow matching models (Qin et al., 2024), and outperforms the diffusion models. Given that MOSES and GUACAMOL benchmarks are approaching saturation, the fact that BWFlow achieves results on par with the SOTA models serves as strong evidence of its effectiveness.

**Table 2: Large molecule generation results. [Table˜15](#A9.T15) gives further experiments with binary edge types.**
|  | Guacamol | MOSES |  |  |  |  |  |  |  |  |  |
| --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- |
| Model | Val. $\uparrow$ | V.U. $\uparrow$ | V.U.N.$\uparrow$ | FCD $\uparrow$ | Val. $\uparrow$ | Unique. $\uparrow$ | Novelty $\uparrow$ | Filters $\uparrow$ | FCD $\downarrow$ | SNN $\uparrow$ | Scaf $\uparrow$ |
| Training set | 100.0 | 100.0 | 0.0 | 92.8 | 100.0 | 100.0 | 0.0 | 100.0 | 0.01 | 0.64 | 99.1 |
| DiGress (Vignac et al., 2023a) | 85.2 | 85.2 | 85.1 | 68.0 | 85.7 | 100.0 | 95.0 | 97.1 | 1.19 | 0.52 | 14.8 |
| DisCo (Xu et al., 2024) | 86.6 | 86.6 | 86.5 | 59.7 | 88.3 | 100.0 | 97.7 | 95.6 | 1.44 | 0.50 | 15.1 |
| Cometh (Siraudin et al., 2024) | 98.9 | 98.9 | 97.6 | 72.7 | 90.5 | 99.9 | 92.6 | 99.1 | 1.27 | 0.54 | 16.0 |
| DeFoG (Qin et al., 2024) | 99.0 | 99.0 | 97.9 | 73.8 | 92.8 | 99.9 | 92.1 | 98.9 | 1.95 | 0.55 | 14.4 |
| BWFlow (Ours) | 98.8 | 98.9 | 97.4 | 69.2 | 92.0 | 100.0 | 94.5 | 98.4 | 1.32 | 0.56 | 15.3 |

### 4.1 Experiment settings

Dataset. For plain graph generation, we evaluate the quality of generated graphs on three benchmark datasets following previous works (Martinkus et al., 2022a; Vignac et al., 2023a; Bergmeister et al., 2024), specifically, planar graphs, tree graphs, and stochastic blocking models (SBM). Two datasets, namely MOSES (Polykovskiy et al., 2018) and GUACAMOL (Brown et al., 2019), are benchmarked to test the model performance on 2D molecule generation. For 3D molecule generation with coordinate data, we test the model on QM9 (Ramakrishnan et al., 2014) and GEOM-DRUGS (Axelrod & Gómez-Bombarelli, 2020).

Metrics. In plain graph generation, the evaluation metrics include the percentage of Valid, Unique, and Novel (V.U.N.) graphs, and the average maximum mean discrepancy ratio (A.Ratio) of graph statistics between the set of generated graphs and the test set are reported (details in [Section˜I.3](#A9.SS3)). For molecule generation, we test two scenarios with and without bond type information, where the latter validates the capacity of our methods in generating the graph structures. To this end, we develop a new relaxed metric to measure the stability and validity of atoms and molecules when bond types are not available. Specifically, the atom-wise stability is relaxed as (Stability of Atom $i$ is defined as $s_{i}$):

$$ $\displaystyle s_{i}=\mathbb{I}[\exists\{b_{ij}\}_{j\in{\mathcal{N}}_{i}}\in\prod_{j\in{\mathcal{N}}_{i}}B_{ij}:\sum_{j\in{\mathcal{N}}_{i}}b_{ij}=\text{EV}_{i}],$ $$

with the identity function $\mathbb{I}$. This means atom $i$ is “relaxed-stable” if there is at least one way to pick allowed bond types ($B_{ij}$) to its neighbors ${\mathcal{N}}_{i}$ so that their total exactly matches the expected valences $\text{EV}_{i}$. Such a relaxed stability of atoms (Atom.Stab.) inherently defines molecule stability (Mol.Stab.) and the validity of a molecule, which are the shared metrics for both 2D/3D molecule generation. In addition to these metrics, distribution metrics are also used for 2D molecules (FCD, Scaf, etc.), and 3D generations (charge distributions, atom total variation, angles, etc.). Details in [Section˜I.3](#A9.SS3).

Setup. To isolate the impact from model architecture, we follow Qin et al. (2024) to fix the backbone model as the same graph transformers. It is shown that sampling/training distortion and target guidance have a significant impact on the performance of graph generation tasks (Qin et al., 2024). In our experiment, the best model performance is obtained with these technologies, but in behavior analysis, we disabled time distortion and target guidance for a fair comparison. In molecular generation, two scenarios with and without bond types information are considered to better evaluate the ability of generating graph structures. More experimental details can be found in [Section˜I.1](#A9.SS1).

### 4.2 Main results for graph generation

Plain graph generation. In [Table˜1](#S4.T1), we report both V.U.N. and A.Ratio. As performance on these benchmarks continues to fluctuate significantly even after convergence and the results are near saturated, we present not only the best scores but also the exponentially moving averaged (EMA) results on last 5 checkpoints and decay 0.999. BWFlow outperforms most competitors on Planar and SBM graph generation. The lone exception is the tree graphs, where our model falls short. We attribute this gap to the fundamentally different generation process for tree graphs (which reside in hyperbolic space (Yang et al., 2022)) and thus violate our MRF assumptions.

#### 3D molecule generation.

[Table˜3](#S4.T3) gives the results on the 3D molecule generation task with explicit hydrogen, where we ignore the bond type but just view the adjacency matrix as a binary one for validating the power of generating graph structures. Interestingly, the empirical results show that even without edge type, the 3D graph generation model already can capture the molecule data distribution. And our BWFlow significantly outperforms the SOTA models, including MiDi (Vignac et al., 2023b) and FlowMol (Dunn & Koes, 2024). We believe a promising future direction is to incorporate the processing of multiple bond types into our framework, which would potentially raise the performance by a margin.

**Table 3: Quantitative experimental results on 3D Molecule Generation with explicit hydrogen.**
| Dataset | Interpolation | Metrics |  |  |  |  |  |  |  |  |
| --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- |
| $\mu$ | V.U.N(%) | Mol.Stab. | Atom.Stab. | Connected(%) | Charge($10^{-2}$) | Atom($10^{-2}$) | Angles(^∘) |  |  |  |
| QM9<br>(with h) | MiDi | 1.01 | 93.13 | 93.98 | 99.60 | 99.21 | 0.2 | 3.7 | 2.21 |  |
| FlowMol | 1.01 | 87.53 | 88.45 | 99.13 | 99.09 | 0.4 | 4.2 | 2.72 |  |  |
| BWFlow | 1.01 | 96.45 | 97.84 | 99.84 | 99.24 | 0.1 | 2.3 | 1.96 |  |  |
| GEOM<br>(with h) | Midi | 1.34 | 78.23 | 32.42 | 89.61 | 79.15 | 0.6 | 11.2 | 9.6 |  |
| FlowMol | 1.34 | 82.20 | 36.90 | 94.60 | 59.98 | 0.4 | 8.8 | 6.5 |  |  |
| BWFlow | 1.20 | 87.75 | 46.80 | 95.08 | 73.53 | 0.1 | 6.5 | 3.96 |  |  |

### 4.3 Behavior analysis

BWFlow provides smooth velocity in probability paths. To illustrate how BWFlow models the smooth evolution of graphs, we compute the A.Ratio on SBM datasets (the figures for the others are in [Figure˜6](#A9.F6)) between generated graph interpolants and test data for $t\in[0,1]$, as shown in [Figure˜3(a)](#S4.F3.sf1). In contrast to the linear (arithmetic) interpolation, BW interpolation initially exposes the model to more out-of-distribution samples with increased A.Ratio. After this early exploration, the A.Ratio monotonously converges, yielding a smooth interpolation between the reference graphs and the data points. This behavior enhances both the model robustness and velocity estimation, which helps in covering the convergence gap in the generation stage as in [Figure˜1(c)](#S1.F1.sf3). In comparison, harmonic and geometric interpolations step outside the valid graph domain, making the learning ill-posed.

The impact of interpolation metrics on the model performance. [Figure˜3(b)](#S4.F3.sf2) illustrated a bar plot that compares interpolation methods on the ability of generating valid plain graphs measured by V.U.N., which shows the superiority of BW interpolation in capturing graph distributions.  [Figure˜3(c)](#S4.F3.sf3) illustrated an example (in planar graph generation) of the convergence curve at the training stage (full results in [Table˜10](#A9.T10)), which suggests that BWFlow can bring a faster convergence speed compared to FM methods constructed with linear (arithmetic) interpolations.

## 5 Discussion and future work

In this paper, we introduce BWFlow, a flow matching model that captures the non-Euclidean and interconnected properties of graphs. While we show BWFlow exhibits outstanding performance in various graph generation tasks, it faces the following limitations that motivate solid future work.

Extension to multiple relation types. As our framework is built upon the interpolation parameterized by the Graph Laplacian, it is not easily generalizable to the graph generation with multiple edge types. We made preliminary attempts at the extension but a comprehensive design is still required.

Lower computational complexity. While constructing the probability path and the velocity, our BW interpolation suffers from an extra cost due to its request to compute the pseudo-inverse of the Laplacian. Compared to linear interpolation, this brings $O(N^{3})$ extra complexity theoretically, and empirically 2x training time and inference time, which is non-negligible in large graph generation. We aim to develop iterative optimization methods to make training faster in future work.

More universal interpolation that accommodates the geometry. In our experiments on the tree dataset, performance was unsatisfactory. We attribute this issue to the unique geometry of tree-structured graphs. A promising future work includes selecting the adaptive interpolation schemes that accommodate the intrinsic geometry of a graph.

## 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 which we feel must be specifically highlighted here.

## References

- Albergo & Vanden-Eijnden (2023)
Albergo, M. S. and Vanden-Eijnden, E.
Building normalizing flows with stochastic interpolants.
In *ICLR*. OpenReview.net, 2023.
- Axelrod & Gómez-Bombarelli (2020)
Axelrod, S. and Gómez-Bombarelli, R.
GEOM: energy-annotated molecular conformations for property prediction and molecular generation.
*CoRR*, abs/2006.05531, 2020.
- Bach et al. (2020)
Bach, E., Rogers, S., Williamson, J., and Rousu, J.
Probabilistic framework for integration of mass spectrum and retention time information in small molecule identification.
*Bioinformatics*, 37(12):1724–1731, 11 2020.
ISSN 1367-4803.
doi: 10.1093/bioinformatics/btaa998.
URL [https://doi.org/10.1093/bioinformatics/btaa998](https://doi.org/10.1093/bioinformatics/btaa998).
- Bergmeister et al. (2024)
Bergmeister, A., Martinkus, K., Perraudin, N., and Wattenhofer, R.
Efficient and scalable graph generation through iterative local expansion.
In *ICLR*. OpenReview.net, 2024.
- Bhatia et al. (2019)
Bhatia, R., Jain, T., and Lim, Y.
On the bures–wasserstein distance between positive definite matrices.
*Expositiones Mathematicae*, 37(2):165–191, 2019.
ISSN 0723-0869.
doi: https://doi.org/10.1016/j.exmath.2018.01.002.
URL [https://www.sciencedirect.com/science/article/pii/S0723086918300021](https://www.sciencedirect.com/science/article/pii/S0723086918300021).
- Bilodeau et al. (2022)
Bilodeau, C., Jin, W., Jaakkola, T., Barzilay, R., and Jensen, K. F.
Generative models for molecular discovery: Recent advances and challenges.
*Wiley Interdisciplinary Reviews: Computational Molecular Science*, 12(5):e1608, 2022.
- Brown et al. (2019)
Brown, N., Fiscato, M., Segler, M. H. S., and Vaucher, A. C.
Guacamol: Benchmarking models for de novo molecular design.
*J. Chem. Inf. Model.*, 59(3):1096–1108, 2019.
- Campbell et al. (2022)
Campbell, A., Benton, J., Bortoli, V. D., Rainforth, T., Deligiannidis, G., and Doucet, A.
A continuous time framework for discrete denoising models.
In *NeurIPS*, 2022.
- Campbell et al. (2024)
Campbell, A., Yim, J., Barzilay, R., Rainforth, T., and Jaakkola, T. S.
Generative flows on discrete state-spaces: Enabling multimodal flows with applications to protein co-design.
In *ICML*. OpenReview.net, 2024.
- Cao & Kipf (2018)
Cao, N. D. and Kipf, T.
Molgan: An implicit generative model for small molecular graphs.
*CoRR*, abs/1805.11973, 2018.
- Chen & Lipman (2024)
Chen, R. T. Q. and Lipman, Y.
Flow matching on general geometries.
In *ICLR*. OpenReview.net, 2024.
- Chen et al. (2018)
Chen, T. Q., Rubanova, Y., Bettencourt, J., and Duvenaud, D.
Neural ordinary differential equations.
In *NeurIPS*, pp.  6572–6583, 2018.
- Chen et al. (2023)
Chen, X., He, J., Han, X., and Liu, L.
Efficient and degree-guided graph generation via discrete diffusion modeling.
In *ICML*, volume 202 of *Proceedings of Machine Learning Research*, pp.  4585–4610. PMLR, 2023.
- Dai et al. (2020)
Dai, H., Nazi, A., Li, Y., Dai, B., and Schuurmans, D.
Scalable deep generative modeling for sparse graphs.
In *ICML*, volume 119 of *Proceedings of Machine Learning Research*, pp.  2302–2312. PMLR, 2020.
- Diamant et al. (2023)
Diamant, N. L., Tseng, A. M., Chuang, K. V., Biancalani, T., and Scalia, G.
Improving graph generation by restricting graph bandwidth.
In *ICML*, volume 202 of *Proceedings of Machine Learning Research*, pp.  7939–7959. PMLR, 2023.
- Dowson & Landau (1982)
Dowson, D. and Landau, B.
The fréchet distance between multivariate normal distributions.
*Journal of Multivariate Analysis*, 12(3):450–455, 1982.
ISSN 0047-259X.
doi: https://doi.org/10.1016/0047-259X(82)90077-X.
URL [https://www.sciencedirect.com/science/article/pii/0047259X8290077X](https://www.sciencedirect.com/science/article/pii/0047259X8290077X).
- Dunn & Koes (2024)
Dunn, I. and Koes, D. R.
Mixed continuous and categorical flow matching for 3d de novo molecule generation.
*ArXiv*, pp.  arXiv–2404, 2024.
- Eijkelboom et al. (2024)
Eijkelboom, F., Bartosh, G., Naesseth, C. A., Welling, M., and van de Meent, J.
Variational flow matching for graph generation.
In *NeurIPS*, 2024.
- Gat et al. (2024)
Gat, I., Remez, T., Shaul, N., Kreuk, F., Chen, R. T. Q., Synnaeve, G., Adi, Y., and Lipman, Y.
Discrete flow matching.
In *NeurIPS*, 2024.
- Goyal et al. (2020)
Goyal, N., Jain, H. V., and Ranu, S.
Graphgen: A scalable approach to domain-agnostic labeled graph generation.
In *WWW*, pp.  1253–1263. ACM / IW3C2, 2020.
- Grover & Leskovec (2016)
Grover, A. and Leskovec, J.
node2vec: Scalable feature learning for networks.
In *KDD*, pp.  855–864. ACM, 2016.
- Haasler & Frossard (2024)
Haasler, I. and Frossard, P.
Bures-wasserstein means of graphs.
In *AISTATS*, volume 238 of *Proceedings of Machine Learning Research*, pp.  1873–1881. PMLR, 2024.
- Haefeli et al. (2022)
Haefeli, K. K., Martinkus, K., Perraudin, N., and Wattenhofer, R.
Diffusion models for graphs benefit from discrete state spaces.
In *The First Learning on Graphs Conference*, 2022.
URL [https://openreview.net/forum?id=CtsKBwhTMKg](https://openreview.net/forum?id=CtsKBwhTMKg).
- Ho et al. (2020)
Ho, J., Jain, A., and Abbeel, P.
Denoising diffusion probabilistic models.
In *NeurIPS*, 2020.
- Hou et al. (2024)
Hou, X., Zhu, T., Ren, M., Bu, D., Gao, X., Zhang, C., and Sun, S.
Improving molecular graph generation with flow matching and optimal transport.
*CoRR*, abs/2411.05676, 2024.
- Ingraham et al. (2019)
Ingraham, J., Garg, V. K., Barzilay, R., and Jaakkola, T. S.
Generative models for graph-based protein design.
In *NeurIPS*, pp.  15794–15805, 2019.
- Jiang et al. (2025)
Jiang, K., Tang, B., Dong, X., and Toni, L.
Heterogeneous graph structure learning through the lens of data-generating processes.
In *The 28th International Conference on Artificial Intelligence and Statistics*, 2025.
URL [https://openreview.net/forum?id=JHK0QBKdYY](https://openreview.net/forum?id=JHK0QBKdYY).
- Jo et al. (2024)
Jo, J., Kim, D., and Hwang, S. J.
Graph generation with diffusion mixture.
In *ICML*. OpenReview.net, 2024.
- Kapusniak et al. (2024)
Kapusniak, K., Potaptchik, P., Reu, T., Zhang, L., Tong, A., Bronstein, M. M., Bose, A. J., and Giovanni, F. D.
Metric flow matching for smooth interpolations on the data manifold.
In *NeurIPS*, 2024.
- Kipf & Welling (2016)
Kipf, T. N. and Welling, M.
Variational graph auto-encoders.
*CoRR*, abs/1611.07308, 2016.
- Li et al. (2023)
Li, M., Kreacic, E., Potluru, V. K., and Li, P.
Graphmaker: Can diffusion models generate large attributed graphs?
*CoRR*, abs/2310.13833, 2023.
- Liao et al. (2019)
Liao, R., Li, Y., Song, Y., Wang, S., Hamilton, W. L., Duvenaud, D., Urtasun, R., and Zemel, R. S.
Efficient graph generation with graph recurrent attention networks.
In *NeurIPS*, pp.  4257–4267, 2019.
- Lipman et al. (2023)
Lipman, Y., Chen, R. T. Q., Ben-Hamu, H., Nickel, M., and Le, M.
Flow matching for generative modeling.
In *The Eleventh International Conference on Learning Representations*, 2023.
URL [https://openreview.net/forum?id=PqvMRDCJT9t](https://openreview.net/forum?id=PqvMRDCJT9t).
- Lipman et al. (2024)
Lipman, Y., Havasi, M., Holderrieth, P., Shaul, N., Le, M., Karrer, B., Chen, R. T. Q., Lopez-Paz, D., Ben-Hamu, H., and Gat, I.
Flow matching guide and code.
*CoRR*, abs/2412.06264, 2024.
- Liu et al. (2023a)
Liu, C., Fan, W., Liu, Y., Li, J., Li, H., Liu, H., Tang, J., and Li, Q.
Generative diffusion models on graphs: Methods and applications.
In *IJCAI*, pp.  6702–6711. ijcai.org, 2023a.
- Liu et al. (2023b)
Liu, X., Gong, C., and Liu, Q.
Flow straight and fast: Learning to generate and transfer data with rectified flow.
In *ICLR*. OpenReview.net, 2023b.
- Martinkus et al. (2022a)
Martinkus, K., Loukas, A., Perraudin, N., and Wattenhofer, R.
SPECTRE: spectral conditioning helps to overcome the expressivity limits of one-shot graph generators.
In *ICML*, volume 162 of *Proceedings of Machine Learning Research*, pp.  15159–15179. PMLR, 2022a.
- Martinkus et al. (2022b)
Martinkus, K., Loukas, A., Perraudin, N., and Wattenhofer, R.
SPECTRE: spectral conditioning helps to overcome the expressivity limits of one-shot graph generators.
In *ICML*, volume 162 of *Proceedings of Machine Learning Research*, pp.  15159–15179. PMLR, 2022b.
- McCann (1997)
McCann, R. J.
A convexity principle for interacting gases.
*Advances in Mathematics*, 128(1):153–179, 1997.
ISSN 0001-8708.
doi: https://doi.org/10.1006/aima.1997.1634.
URL [https://www.sciencedirect.com/science/article/pii/S0001870897916340](https://www.sciencedirect.com/science/article/pii/S0001870897916340).
- Minello et al. (2025)
Minello, G., Bicciato, A., Rossi, L., Torsello, A., and Cosmo, L.
Generating graphs via spectral diffusion.
In *The Thirteenth International Conference on Learning Representations*, 2025.
URL [https://openreview.net/forum?id=AAXBfJNHDt](https://openreview.net/forum?id=AAXBfJNHDt).
- Niu et al. (2020)
Niu, C., Song, Y., Song, J., Zhao, S., Grover, A., and Ermon, S.
Permutation invariant graph generation via score-based generative modeling.
In *AISTATS*, volume 108 of *Proceedings of Machine Learning Research*, pp.  4474–4484. PMLR, 2020.
- Olkin & Pukelsheim (1982)
Olkin, I. and Pukelsheim, F.
The distance between two random vectors with given dispersion matrices.
*Linear Algebra and its Applications*, 48:257–263, 1982.
ISSN 0024-3795.
doi: https://doi.org/10.1016/0024-3795(82)90112-4.
URL [https://www.sciencedirect.com/science/article/pii/0024379582901124](https://www.sciencedirect.com/science/article/pii/0024379582901124).
- Polykovskiy et al. (2018)
Polykovskiy, D., Zhebrak, A., Sánchez-Lengeling, B., Golovanov, S., Tatanov, O., Belyaev, S., Kurbanov, R., Artamonov, A., Aladinskiy, V., Veselov, M., Kadurin, A., Nikolenko, S. I., Aspuru-Guzik, A., and Zhavoronkov, A.
Molecular sets (MOSES): A benchmarking platform for molecular generation models.
*CoRR*, abs/1811.12823, 2018.
- Pooladian et al. (2023)
Pooladian, A., Ben-Hamu, H., Domingo-Enrich, C., Amos, B., Lipman, Y., and Chen, R. T. Q.
Multisample flow matching: Straightening flows with minibatch couplings.
In *ICML*, volume 202 of *Proceedings of Machine Learning Research*, pp.  28100–28127. PMLR, 2023.
- Qin et al. (2024)
Qin, Y., Madeira, M., Thanou, D., and Frossard, P.
Defog: Discrete flow matching for graph generation.
*CoRR*, abs/2410.04263, 2024.
- Qu et al. (2019)
Qu, M., Bengio, Y., and Tang, J.
GMNN: graph markov neural networks.
In *ICML*, volume 97 of *Proceedings of Machine Learning Research*, pp.  5241–5250. PMLR, 2019.
- Ramakrishnan et al. (2014)
Ramakrishnan, R., Dral, P. O., Rupp, M., and Von Lilienfeld, O. A.
Quantum chemistry structures and properties of 134 kilo molecules.
*Scientific data*, 1(1):1–7, 2014.
- Siraudin et al. (2024)
Siraudin, A., Malliaros, F. D., and Morris, C.
Cometh: A continuous-time discrete-state graph diffusion model, 2024.
URL [https://arxiv.org/abs/2406.06449](https://arxiv.org/abs/2406.06449).
- Song et al. (2021)
Song, Y., Sohl-Dickstein, J., Kingma, D. P., Kumar, A., Ermon, S., and Poole, B.
Score-based generative modeling through stochastic differential equations.
In *ICLR*. OpenReview.net, 2021.
- Stärk et al. (2024)
Stärk, H., Jing, B., Wang, C., Corso, G., Berger, B., Barzilay, R., and Jaakkola, T. S.
Dirichlet flow matching with applications to DNA sequence design.
In *ICML*. OpenReview.net, 2024.
- Sun et al. (2023)
Sun, H., Yu, L., Dai, B., Schuurmans, D., and Dai, H.
Score-based continuous-time discrete diffusion models.
In *ICLR*. OpenReview.net, 2023.
- Takatsu (2010)
Takatsu, A.
On wasserstein geometry of gaussian measures.
*Probabilistic approach to geometry*, 57:463–472, 2010.
- Taskar et al. (2007)
Taskar, B., Abbeel, P., Wong, M.-F., and Koller, D.
Relational markov networks.
*Introduction to statistical relational learning*, 175:200, 2007.
- Tong et al. (2024)
Tong, A., Fatras, K., Malkin, N., Huguet, G., Zhang, Y., Rector-Brooks, J., Wolf, G., and Bengio, Y.
Improving and generalizing flow-based generative models with minibatch optimal transport.
*Trans. Mach. Learn. Res.*, 2024, 2024.
- Vignac et al. (2023a)
Vignac, C., Krawczuk, I., Siraudin, A., Wang, B., Cevher, V., and Frossard, P.
Digress: Discrete denoising diffusion for graph generation.
In *ICLR*. OpenReview.net, 2023a.
- Vignac et al. (2023b)
Vignac, C., Osman, N., Toni, L., and Frossard, P.
Midi: Mixed graph and 3d denoising diffusion for molecule generation.
In *ECML/PKDD (2)*, volume 14170 of *Lecture Notes in Computer Science*, pp.  560–576. Springer, 2023b.
- Villani & Society (2003)
Villani, C. and Society, A. M.
*Topics in Optimal Transportation*.
Graduate studies in mathematics. American Mathematical Society, 2003.
ISBN 9781470418045.
URL [https://books.google.co.uk/books?id=MyPjjgEACAAJ](https://books.google.co.uk/books?id=MyPjjgEACAAJ).
- Wang et al. (2024)
Wang, F., Yang, L., Huang, Z., Wang, M., and Li, H.
Rectified diffusion: Straightness is not your need in rectified flow.
*CoRR*, abs/2410.07303, 2024.
- Weigt et al. (2009)
Weigt, M., White, R. A., Szurmant, H., Hoch, J. A., and Hwa, T.
Identification of direct residue contacts in protein–protein interaction by message passing.
*Proceedings of the National Academy of Sciences*, 106(1):67–72, 2009.
doi: 10.1073/pnas.0805923106.
URL [https://www.pnas.org/doi/abs/10.1073/pnas.0805923106](https://www.pnas.org/doi/abs/10.1073/pnas.0805923106).
- Xu et al. (2024)
Xu, Z., Qiu, R., Chen, Y., Chen, H., Fan, X., Pan, M., Zeng, Z., Das, M., and Tong, H.
Discrete-state continuous-time diffusion for graph generation.
*arXiv preprint arXiv:2405.11416*, 2024.
- Yang et al. (2022)
Yang, M., Zhou, M., Li, Z., Liu, J., Pan, L., Xiong, H., and King, I.
Hyperbolic graph neural networks: A review of methods and applications.
*CoRR*, abs/2202.13852, 2022.
- You et al. (2018)
You, J., Ying, R., Ren, X., Hamilton, W. L., and Leskovec, J.
Graphrnn: Generating realistic graphs with deep auto-regressive models.
In *ICML*, volume 80 of *Proceedings of Machine Learning Research*, pp.  5694–5703. PMLR, 2018.
- Yu & Zhan (2025)
Yu, M. and Zhan, K.
Bias mitigation in graph diffusion models.
In *The Thirteenth International Conference on Learning Representations*, 2025.
URL [https://openreview.net/forum?id=CSj72Rr2PB](https://openreview.net/forum?id=CSj72Rr2PB).
- Zhu et al. (2003)
Zhu, X., Lafferty, J., and Ghahramani, Z.
*Semi-supervised learning: From Gaussian fields to Gaussian processes*.
School of Computer Science, Carnegie Mellon University, 2003.
- Zhu et al. (2022)
Zhu, Y., Du, Y., Wang, Y., Xu, Y., Zhang, J., Liu, Q., and Wu, S.
A survey on deep graph generation: Methods and applications.
In *LoG*, volume 198 of *Proceedings of Machine Learning Research*, pp.  47. PMLR, 2022.

## Appendix A Graph Markov random fields: background and theory

### A.1 Background of Markov random fields

Markov random fields (MRFs) were originally developed to describe the dynamics of interconnected physical systems such as molecules and proteins (Weigt et al., 2009; Bach et al., 2020). For a given graph $G=\{V,E\}$, MRFs are energy-based models that describe the graph with the following probability density:

$$ $p(G)=\frac{1}{Z}\prod_{\xi\in\mathrm{cl}(G)}\varphi_{\xi}\left({\mathcal{V}}_{\xi}\right)=\frac{1}{Z}e^{-U(G)/kT},$ (15) $$

where the energy $U(G)$ is used to describe the whole connected system, $\operatorname{cl}(G)$ is the set of cliques of $G$, ${\mathcal{V}}_{\xi}$ is the subset of nodes related to $\xi$, and $Z$ is the partition function. In our paper, we follow Zhu et al. (2003) and parameterize the energy function with node-wise potential $\varphi_{1}(v),\forall v\in{\mathcal{V}}$ and edge-wise potential $\varphi_{2}(u,v)$, $\forall e_{uv}\in{\mathcal{E}}$:

$$ $U(G)=\sum_{v\in V}\underbrace{\log\varphi_{1}(v)}_{\text{Node-wise potential}}+\sum_{\{u,v\}\in E}\underbrace{\log\varphi_{2}(u,v)}_{\text{Edge-wise potential}}$ $$

As a concrete example, in the molecule system that consists of atoms and bonds, the node-wise potential can be the kinetic energy for an atom and the edge-wise potential is determined by the interatomic forces, such as the Electrostatic force. MRFs thus serve as a natural and elegant way to describe general graph systems.

The energy-based models have an intrinsic relationship with generative models. As an example, Song et al. (2021) derived the relationship between diffusion models and the Langevin dynamics, which is used to describe the evolution of an energy-based model. It is shown that the diffusion models are trying to approximate the score function $\nabla_{\mathcal{X}}\operatorname{log}p({\mathcal{X}})$. In the energy-based models, the score function is just the gradient of energy, $\nabla_{\mathcal{X}}\operatorname{log}p({\mathcal{X}})=-\nabla_{\mathcal{X}}U({\mathcal{X}})$, and the Langevin dynamics samples data points towards a lower energy profile state. Thus, interpolating between two graph distributions, respectively, reference and data distributions, can be viewed as transitioning between two Markov random fields with different energy landscapes. As illustrated in  [Figure˜4](#A1.F4), the random molecule graph sampled from the reference distributions corresponds to the high energy state and the data distribution the low energy state.

Figure: Figure 4: Molecule graphs correspond to MRFs with different energies.
Refer to caption: /html/2506.14020/assets/x5.png

The idea of our paper originated from the two facts: MRFs are energy-based model describing connected systems, and the energy-based models have their intrinsic relationship with the diffusion/flow models. Thus, if a model is required to describe the evolution of the whole graph system, we believe it is natural to consider constructing a probability path for two graph distributions with MRFs as the backbone.

### A.2 Derivation of Graph Markov random fields

We show the derivation of [Definition˜2](#Thmdefinition2), which is restated here:

###### Definition 3 (Graph Markov Random Fields) .

Derivation:

We start from

$$ $\displaystyle p({\mathcal{X}};{\bm{X}},{\bm{W}})\propto\prod_{v}\exp\left\{-(\nu+d_{v})\|{\bm{V}}x_{v}-\bm{\mu}_{v}\|^{2}\right\}\prod_{u,v}\exp\left\{w_{uv}\left[({\bm{V}}x_{u}-\bm{\mu}_{u})^{\top}({\bm{V}}x_{v}-\bm{\mu}_{v})\right]\right\}.$ (17) $$

We assume that the linear transformation matrix has dimension ${\bm{V}}\in\mathbb{R}^{K^{\prime}\times K}$ given that $x_{v}\in\mathbb{R}^{K}$ and define a transformed variable

$$ $h_{v}\equiv{\bm{V}}\,x_{v}-\bm{\mu}_{v}\in\mathbb{R}^{K^{\prime}},\text{ stacking as }{\mathcal{H}}\in\mathbb{R}^{|{\mathcal{V}}|\times K^{\prime}}.$ (18) $$

The probability becomes

$$ $P({\mathcal{H}};{\bm{X}},{\bm{W}})\propto\prod_{v}\exp\Bigl{\{}-(\nu+d_{v})\|h_{v}\|^{2}\Bigr{\}}\prod_{u,v}\exp\Bigl{\{}w_{uv}\,h_{u}^{\top}h_{v}\Bigr{\}}.$ (19) $$

Then, the terms inside the exponent in [Equation˜19](#A1.E19) become

$$ $\displaystyle-\sum_{v}(\nu+d_{v})\|h_{v}\|^{2}+\sum_{u,v}w_{uv}\,h_{u}^{\top}h_{v}$ $\displaystyle=-\sum_{v}(\nu+d_{v})\,h_{v}^{\top}h_{v}+\sum_{u,v}w_{uv}\,h_{u}^{\top}h_{v}$ $\displaystyle=-\sum_{u,v}h_{u}^{\top}\Bigl{[}(\nu+d_{u})\delta_{uv}-w_{uv}\Bigr{]}h_{v},$ $$

where the Kronecker delta $\delta_{uv}=1$ if $u=v$ and 0 else. We define a squared matrix $\Lambda^{\prime}$ to arrange the inner term, which can be written as,

$$ $\bm{\Lambda}^{\prime}=\nu{\bm{I}}+{\bm{L}}\quad\text{with}\quad\bm{\Lambda}^{\prime}_{uv}=(\nu+d_{u})\delta_{uv}-w_{uv}.$ (20) $$

${\bm{I}}$ is the identity matrix. Thus, the exponent in compact matrix form gives

$$ $-\frac{1}{2}\operatorname{Tr}({\mathcal{H}}^{\top}\bm{\Lambda}^{\prime}{\mathcal{H}}),\text{ where }{\mathcal{H}}=\begin{pmatrix}h_{1}\\ h_{2}\\ \vdots\\ h_{|{\mathcal{V}}|}\end{pmatrix}.$ (21) $$

It is possible to rearrange the exponent as

$$ $\operatorname{Tr}({\mathcal{H}}^{\top}\bm{\Lambda}^{\prime}{\mathcal{H}})=\operatorname{vec}({\mathcal{H}})^{\top}(\bm{\Lambda}\otimes{\bm{I}})\operatorname{vec}({\mathcal{H}}),$ (22) $$

where $\otimes$ denotes the Kronecker product. This is exactly in the form of a multivariate colored Gaussian. Thus, the joint distribution of $\operatorname{vec}({\mathcal{H}})$ (of dimension $|{\mathcal{V}}|K^{\prime}$) is given by

$$ $\operatorname{vec}({\mathcal{H}})\sim\mathcal{N}\Bigl{(}0,\;\,(\nu{\bm{I}}+{\bm{L}}))^{-1}\otimes{\bm{I}}_{K^{\prime}}\Bigr{)},$ (23) $$

Recall that $h_{v}={\bm{V}}\,x_{v}-\bm{\mu}_{v}$, we obtain

$$ $\operatorname{vec}({\mathcal{H}})=({\bm{I}}\otimes{\bm{V}})\,\operatorname{vec}({\mathcal{X}})-\operatorname{vec}(\bm{\mu}).$ (24) $$

Since the transformation is linear, the distribution over ${\mathcal{X}}$ remains Gaussian. By the properties of linear transformations of Gaussians, if

$$ $\operatorname{vec}({\mathcal{H}})\sim\mathcal{N}(\operatorname{vec}(\bm{\mu}),\,{\bm{\Sigma}}_{h}),\operatorname{vec}({\mathcal{X}})=({\bm{I}}\otimes{\bm{V}}^{\dagger})\,\operatorname{vec}({\mathcal{H}}),$ (25) $$

then

$$ $\operatorname{vec}({\mathcal{X}})\sim\mathcal{N}\Bigl{(}({\bm{I}}\otimes{\bm{V}}^{\dagger})\,\operatorname{vec}(\bm{\mu}),\;({\bm{I}}_{n}\otimes{\bm{V}}^{\dagger})\,{\bm{\Sigma}}_{\mathcal{H}}\,({\bm{I}}_{n}\otimes{\bm{V}}^{\dagger})^{\top}\Bigr{)}.$ (26) $$

Thus, using the mixed-product property of the Kronecker product,

$$ $({\bm{I}}\otimes{\bm{V}}^{\dagger})((\nu{\bm{I}}+{\bm{L}})^{-1}\otimes{\bm{I}})({\bm{I}}_{n}\otimes{\bm{V}}^{\dagger})^{\top}=({\bm{L}}+\nu{\bm{I}})^{-1}\otimes({\bm{V}}^{\dagger}{\bm{V}}^{\dagger\top})$ (27) $$

Finally, the joint distribution over ${\mathcal{X}}$ is

$$ $\displaystyle\operatorname{vec}({\mathcal{X}})\sim{\mathcal{N}}({\bm{X}},\Sigma),$ (28) with $\displaystyle{\bm{X}}=({\bm{I}}\otimes{\bm{V}}^{\dagger})\,\operatorname{vec}(\bm{\mu})=\operatorname{vec}({\bm{V}}^{\dagger}\bm{\mu})$ and $\displaystyle\Sigma=\,(\nu{\bm{I}}+{\bm{L}})^{-1}\otimes\left({\bm{V}}^{\dagger}{\bm{V}}^{\dagger\top}\right),$ $$

We use the following lemma:

###### Lemma 1 .

Given two invertible matrices ${\bm{A}}$ and ${\bm{B}}$, their Kronecker product satisfies $({\bm{A}}\otimes{\bm{B}})^{-1}={\bm{A}}^{-1}\otimes{\bm{B}}^{-1}$.

So that we get

$$ $\operatorname{vec}({\mathcal{X}})\sim\mathcal{N}\left({\bm{X}},\Lambda^{\dagger}\right),\text{ with }{\bm{X}}=\operatorname{vec}({\bm{V}}^{\dagger}\bm{\mu}),\Lambda=(\nu{\bm{I}}+{\bm{L}})\otimes\bm{V}^{\top}\bm{V}.$ (29) $$

which ends the derivation.

###### Definition 3 (Graph Markov Random Fields) .

###### Lemma 1 .

### A.3 The usage scope of graph Markov random fields

Given that our Graph Markov random fields (GraphMRF) have an explicit form to constrain the graph distribution, it inherits certain inductive biases and we have to properly understand the usage scenarios.

To understand the scenarios which GraphMRF could be used, we start by stating the plain form of MRF by melting the linear transformer matrix and the mean term, i.e., giving $h_{v}={\bm{V}}x_{v}-\bm{\mu}$ which gives the probability density as,

$$ $P({\mathcal{H}}\mid{\bm{L}})\propto\exp(-\operatorname{trace}({\mathcal{H}}^{\top}({\bm{L}}+\nu{\bm{I}}){\mathcal{H}}))=\exp(-{\bm{w}}_{uv}\|h_{u}-h_{v}\|^{2}-\nu\sum_{u}\|h_{u}\|^{2}_{\text{F}})$ (30) $$

#### Plain Graph Markov random fields.

To understand the scenarios which we can utilize MRF to model graphs, we first consider the simplest case when ${\bm{V}}$ is rectangular orthogonal (semi-orthogonal) matrix such that ${\bm{V}}^{\top}{\bm{V}}={\bm{I}}$ and the mean $\bm{\mu}=0$, the probability density becomes,

$$ $P({\bm{X}},{\bm{L}})\propto\exp(-{\bm{X}}^{\top}({\bm{L}}+\nu{\bm{I}}){\bm{X}})=\exp(-\sum_{\{u,v\}\in{\mathcal{E}}}{\bm{W}}_{uv}({\bm{x}}_{u}-{\bm{x}}_{v})^{2}-\nu\sum_{u}{\bm{x}}_{u}^{2})$ (31) $$

As $\nu\rightarrow 0$, the exponent term inside becomes

$$ ${\mathcal{S}}({\bm{X}},{\bm{L}})=-\sum_{\{u,v\}\in{\mathcal{E}}}{\bm{W}}_{uv}({\bm{x}}_{u}-{\bm{x}}_{v})^{2}=\operatorname{trace}({\bm{X}}^{\top}{\bm{L}}{\bm{X}}),$ (32) $$

where we name ${\mathcal{S}}({\bm{X}},{\bm{L}})$ as the smoothness of the graph features. The smoothness measures how similar the neighbors connected are. For instance, if there exists an edge between node $u$ and $v$ weighted as ${\bm{W}}_{uv}$, the likelihood will be higher if ${\bm{x}}_{u}$ and ${\bm{x}}_{v}$ be similar, so that $\|{\bm{x}}_{u}-{\bm{x}}_{v}\|^{2}$ are small. This suggests that the probability will be higher if the ${\mathcal{S}}({\bm{X}},{\bm{L}})$ is small.

This vanilla form captures the first type of graphs the GraphMRF can model - the homophily graphs, i.e., similar nodes (measured by the node attribute) may be more likely to attach to each other than dissimilar ones. This includes the social networks with friendship, traffic networks, etc.

#### Graph Markov random fields with embeddings.

Now we move one step further to consider the graphs with linear transformer matrix ${\bm{V}}$. Linear transformation provides a map from the feature space to the latent space, which can be considered as an embedding method to empower the models with better expressiveness. As a simple example, when the ${\bm{V}}$ provides a negative projection, the mapping can capture the heterophily relationships, which means the nodes connected are dissimilar.

Coincidently, this aligns well with the famous embedding method Node2Vec as in Grover & Leskovec (2016), where the edge weights are proportional to the negative distance, or the inner product of the embeddings. i.e.,

$$ ${\bm{W}}_{uv}\propto\exp(-\|{\bm{V}}{\bm{x}}_{u}-{\bm{V}}{\bm{x}}_{v}\|_{\text{F}}^{2})$ (33) $$

In (Jiang et al., 2025) it is derived that learning the parameters of MRFs is intrinsically equivalent to learning embeddings similar to Node2Vec. As such, the expressiveness of MRFs are as good as Node2Vec, which grants its usage to molecule graphs, protein interaction networks, social networks, and knowledge graphs. In our paper we make the assumption is that the linear mapping from ${\bm{X}}$ the observation is shared. This requirement translates to that the two graphs should have the same embedding space and feature space, which is practical if the reference distribution and data distributions share the same space.

#### Graphs without features.

We wish to emphasize that even though the GraphMRF is constructed under the assumption that graph features exist, it is capable of modeling the non-attributed graphs, such as planar and SBM graphs. To do so, we consider the optimization over the Rayleigh function: It is shown that, if $v_{1},\ldots,v_{k-1}$ are orthonormal eigenvectors for $\lambda_{1},\ldots,\lambda_{k-1}$, then the eigenvalues satisfy,

$$ $\lambda_{k}=\min_{\begin{subarray}{c}{\bm{x}}\neq 0\\ {\bm{x}}\perp{\bm{v}}_{1},\ldots,{\bm{v}}_{k-1}\end{subarray}}R({\bm{x}}),\text{ with }R({\bm{x}})=\frac{{\bm{x}}^{T}{\bm{L}}{\bm{x}}}{{\bm{x}}^{T}{\bm{x}}}$ (34) $$

In such a scenario, the graphs are no longer related to the actual node features, but instead, the eigenvectors ${\bm{v}}_{k}$ serve as an intrinsic graph feature. Interestingly, $R({\bm{x}})$ is the normalized form of the smoothness as in [Equation˜31](#A1.E31), which is exactly equivalent to the eigenvalues. This means that if the spectrum of the graph Laplacian focuses on the low-frequency components, the corresponding graphs would give a higher probability in terms of MRFs if we view the eigenvectors as the node features. Such a pattern could capture the graph distributions almost all the synthetic datasets, such as Planar, SBM, TLS, COMM20 datasets, satisfy such a property. There is no surprise that our models capture better global patterns in those datasets. But it is also worth pointing out that there exists exceptions, such as the tree graph, which does not have a clear clustering pattern, thus the spectrum does not follow our GraphMRF.

#### Future work.

A limitation of our method is that it cannot easily capture the generation of graphs with multiple relation types, which we name heterogeneous graphs. Even though we utilize an intuitive solution in the experiment to produce [Table˜2](#S4.T2): we first sample the pure graph structure without edge types to produce the graph backbone, and then sample the edge types via liner interpolated probability on top of the backbone. The solution provides preliminary results for the graph generation in multi-relational graphs, but still requires improvements. Fortunately, there exists a few ways to extend the GraphMRF to heterogeneous graphs (Jiang et al., 2025). An interesting future work can be generalizing our model to heterogeneous graphs by considering GraphMRF variants, such as the H2MN proposed in Jiang et al. (2025).

## Appendix B Proofs

### B.1 Wasserstein Distance between two colored gaussian distributions

We first prove the lemma that captures the Wasserstein distance between two colored Gaussians, which will be used in deriving our Bures-Wasserstein distances in graph generations.

###### Lemma 2 .

Proof.
We first state the following proposition.

###### Proposition 4 .

(Translation Invariance of the 2-Wasserstein Distance for Gaussian Measures) Consider two measures $\eta_{0}\sim\mathcal{N}\left(\bm{\mu}_{0},\Sigma_{0}\right)$ and $\eta_{1}\sim\mathcal{N}\left(\bm{\mu}_{1},\Sigma_{1}\right)$ and their centered measure as $\tilde{\eta}_{0}=\mathcal{N}\left(0,\Sigma_{0}\right)$ and $\tilde{\eta}_{1}=\mathcal{N}\left(0,\Sigma_{1}\right)$, the squared Wasserstein distance decomposes as

$$ ${\mathcal{W}}_{2}^{2}\left(\eta_{0},\eta_{1}\right)=\left\|\bm{\mu}_{0}-\bm{\mu}_{1}\right\|_{2}^{2}+{\mathcal{W}}_{2}^{2}\left(\tilde{\eta}_{0},\tilde{\eta}_{1}\right)$ $$

Proof:

Consider two random vectors ${\mathcal{X}},{\mathcal{Y}}$ distributed as $\eta_{0},\eta_{1}$,

$$ ${\mathcal{X}}=\bm{\mu}_{0}+\tilde{{\mathcal{X}}},{\mathcal{Y}}=\mu_{1}+\tilde{{\mathcal{Y}}},\text{ with }\tilde{{\mathcal{X}}}\sim\tilde{\eta}_{0},\tilde{{\mathcal{Y}}}\sim\tilde{\eta}_{1}.$ $$

For any coupling $({\mathcal{X}},{\mathcal{Y}})$, we consider the expected squared Euclidean distance,

$$ $\displaystyle\mathbb{E}_{{\mathcal{X}},{\mathcal{Y}}}\|{\mathcal{X}}-{\mathcal{Y}}\|^{2}$ $\displaystyle=\mathbb{E}_{{\mathcal{X}},{\mathcal{Y}}}\left\|\bm{mu}_{0}-\bm{\mu}_{1}+(\tilde{{\mathcal{X}}}-\tilde{{\mathcal{Y}}})\right\|^{2}.$ (35) $\displaystyle=\left\|\bm{\mu}_{0}-\bm{\mu}_{1}\right\|^{2}+2\left\langle\bm{\mu}_{0}-\mu_{1},\tilde{{\mathcal{X}}}-\tilde{{\mathcal{Y}}}\right\rangle+\mathbb{E}_{\tilde{{\mathcal{X}}},\tilde{{\mathcal{Y}}}}\|\tilde{{\mathcal{X}}}-\tilde{{\mathcal{Y}}}\|^{2}$ $$

Since $\tilde{{\mathcal{X}}}$ and $\tilde{{\mathcal{Y}}}$ both have zero mean, we have $\mathbb{E}[\tilde{{\mathcal{X}}}-\tilde{{\mathcal{Y}}}]=0$ so the cross-term vanishes. Thus,

$$ $\mathbb{E}\|{\mathcal{X}}-{\mathcal{Y}}\|^{2}=\left\|\bm{\mu}_{0}-\bm{\mu}_{1}\right\|^{2}+\mathbb{E}\|\tilde{{\mathcal{X}}}-\tilde{{\mathcal{Y}}}\|^{2}$ (36) $$

Take the definition of 2-Wasserstein distance, the infimum over all couplings directly yields

$$ $\displaystyle\left({\mathcal{W}}_{2}(\eta_{0},\eta_{1})\right)^{2}$ $\displaystyle=\inf_{\pi\in\Pi(\eta_{0},\eta_{1})}\int\|{\mathcal{X}}-{\mathcal{Y}}\|^{2}\,d\pi({\mathcal{X}},{\mathcal{Y}}).$ (37) $\displaystyle=\left\|\bm{\mu}_{0}-\bm{\mu}_{1}\right\|^{2}+{\mathcal{W}}_{2}^{2}\left(\tilde{\eta}_{0},\tilde{\eta}_{1}\right)$ $$

This completes the proof of Proposition [4](#Thmprop4).

Now we prove the flowing proposition, which will give us our lemma.

###### Proposition 5 .

Given two centered measures as $\tilde{\eta}_{0}=\mathcal{N}\left(0,\Sigma_{0}\right)$ and $\tilde{\eta}_{1}=\mathcal{N}\left(0,\Sigma_{1}\right)$

$$ ${\mathcal{W}}_{2}^{2}\left(\tilde{\eta}_{0},\tilde{\eta}_{1}\right)=\operatorname{Tr}\left(\Sigma_{0}+\Sigma_{1}-2\left(\Sigma_{1}^{1/2}\Sigma_{0}\Sigma_{1}^{1/2}\right)^{1/2}\right).$ (38) $$

proof. The coupling $\pi$ of $\tilde{\eta}_{0}$ and $\tilde{\eta}_{1}$ is a joint Gaussian measure with zero mean and covariance matrix

$$ $\Sigma_{c}=\begin{pmatrix}\Sigma_{0}&{\bm{C}}\\ {\bm{C}}^{T}&\Sigma_{1}\end{pmatrix}\succeq 0,$ (39) $$

where ${\bm{C}}$ is the cross-covariance and $\succeq$ means the matrix is positive semi-definitive (PSD). The expected squared distance between the two random vectors $({\mathcal{X}},{\mathcal{Y}})$ drawn from $\pi$ is then described as,

$$ $\displaystyle\mathbb{E}\|{\mathcal{X}}-{\mathcal{Y}}\|^{2}$ $\displaystyle=\operatorname{Tr}(\mathbb{E}[({\mathcal{X}}-{\mathcal{Y}})({\mathcal{X}}-{\mathcal{Y}})^{\top}])$ (40) $\displaystyle=\operatorname{Tr}(\Sigma_{0})+\operatorname{Tr}(\Sigma_{1})-2\,\operatorname{Tr}({\bm{C}}).$ $$

The definition of Wasserstein distance gives,

$$ ${\mathcal{W}}_{c}(\eta_{0},\eta_{1})=\inf_{\pi\in\Pi(\eta_{0},\eta_{1})}\mathbb{E}\|{\mathcal{X}}-{\mathcal{Y}}\|^{2}$ (41) $$

Thus, minimizing the Wasserstein distance is equivalent to maximizing $\operatorname{Tr}({\bm{C}})$ over all ${\bm{C}}$ subject to the joint covariance is positive semi-definite (PSD).
It turns out (see Dowson & Landau (1982); Olkin & Pukelsheim (1982); Takatsu (2010) ) that the condition in [Equation˜39](#A2.E39) is equivalent to,

$$ $\Sigma_{1}-{\bm{C}}^{\top}\Sigma_{0}^{-1}{\bm{C}}\succeq 0\leftrightarrow\Sigma_{0}^{-1/2}{\bm{C}}\Sigma_{1}^{-1/2}\text{ has operator norm }\leq 1$ (42) $$

So we denote $\,{\bm{K}}:=\Sigma_{0}^{-1/2}{\bm{C}}\Sigma_{1}^{-1/2}$ with $\|{\bm{K}}\|_{\text{op}}\leq 1$. Then

$$ $\operatorname{Tr}({\bm{C}})=\operatorname{Tr}\left(\Sigma_{0}^{1/2}{\bm{K}}\Sigma_{1}^{1/2}\right)=\operatorname{Tr}\left({\bm{K}}\Sigma_{1}^{1/2}\Sigma_{0}^{1/2}\right).$ $$

Using von Neumann trace inequality, its trace inner-product with ${\bm{K}}$ is maximized by choosing ${\bm{K}}={\bm{I}}$ on the support.

$$ $\max_{\|{\bm{K}}\|_{op}\leq 1}\operatorname{Tr}({\bm{K}}{\bm{A}})=\operatorname{Tr}\left({\bm{M}}^{1/2}\right),\quad{\bm{M}}=\sqrt{{\bm{A}}{\bm{A}}^{\top}}=\Sigma_{1}^{1/2}\Sigma_{0}\Sigma_{1}^{1/2}$ $$

Hence the optimal value of $\operatorname{Tr}({\bm{C}})$ is

$$ $\operatorname{Tr}({\bm{C}}^{*})=\operatorname{Tr}\left[\left(\Sigma_{1}^{1/2}\Sigma_{0}\,\Sigma_{1}^{1/2}\right)^{1/2}\right]$ $$

Substituting this optimal value into the expression of Wasserstein distance, we obtain

$$ $\displaystyle{\mathcal{W}}_{2}^{2}\left(\tilde{\eta}_{0},\tilde{\eta}_{1}\right)$ $\displaystyle=\operatorname{Tr}(\Sigma_{0})+\operatorname{Tr}(\Sigma_{1})-2\,\operatorname{Tr}\left[\left(\Sigma_{1}^{1/2}\Sigma_{0}\,\Sigma_{1}^{1/2}\right)^{1/2}\right].$ (43) $$

This completes the proof of proposition [5](#Thmprop5). Taking Proposition [4](#Thmprop4) and Proposition [5](#Thmprop5) together, we proved [Lemma˜2](#Thmlemma2).

###### Lemma 2 .

###### Proposition 4 .

###### Proposition 5 .

### B.2 Derivation of the graph Wasserstein distance under MRF

We then prove the Bures-Wasserstein distance for two graph distributions. We restate [˜1](#Thmprop1),

###### Proposition 6 (Bures-Wasserstein Distance) .

Specifically, [Definition˜2](#Thmdefinition2) uses graph Markov random fields to describe a graph as

$$ $\displaystyle p({\mathcal{G}};{\bm{G}})=p({\mathcal{X}},{\mathcal{E}};{\bm{X}},{\bm{W}})=p({\mathcal{X}};{\bm{X}},{\bm{W}})\cdot p({\mathcal{E}};{\bm{W}})\text{ where }{\mathcal{E}}\sim\delta({\bm{W}})\text{ and }$ (45) $\displaystyle\operatorname{vec}({\mathcal{X}})\sim\mathcal{N}\left({\bm{X}},\bm{\Lambda}^{\dagger}\right),\text{ with }{\bm{X}}=\operatorname{vec}({\bm{V}}^{\dagger}\bm{\mu}),\bm{\Lambda}=(\nu{\bm{I}}+{\bm{L}})\otimes\bm{V}^{\top}\bm{V}.$ $$

With the graph Wasserstein distance defined as,

$$ $\text{(Graph Wasserstein Distance)}\quad d_{\text{BW}}({\mathcal{G}}_{0},{\mathcal{G}}_{1}):={\mathcal{W}}_{c}\left(\eta_{{\mathcal{G}}_{0}},\eta_{{\mathcal{G}}_{1}}\right)={\mathcal{W}}_{c}(\eta_{{\mathcal{X}}_{0}},\eta_{{\mathcal{X}}_{1}})+{\mathcal{W}}_{c}(\eta_{{\mathcal{E}}_{0}},\eta_{{\mathcal{E}}_{1}}).$ $$

We first consider calculating ${\mathcal{W}}_{c}(\eta_{{\mathcal{X}}_{0}},\eta_{{\mathcal{X}}_{1}})$. Specifically, this is the distance between two colored Gaussian measures where

$$ $\begin{gathered}\eta_{i}\sim\mathcal{N}\Bigl{(}\bm{\mu}_{i}^{\prime},\Sigma_{i}\Bigr{)},\quad i=0,1,\[2.84526pt] \text{where }\bm{\mu}_{i}^{\prime}={\bm{V}}_{i}\otimes\bm{\mu}_{i}\text{ and }\Sigma_{i}^{-1}=\Lambda_{i}=(\nu{\bm{I}}+{\bm{L}}_{i})\otimes({\bm{V}}_{i}^{\top}{\bm{V}}_{i}).\end{gathered}$ (46) $$

where we first assume that these two Gaussians are emitted from different linear transformation matrices ${\bm{V}}_{0}$ and ${\bm{V}}_{1}$. This will bring us the most general and flexible form that could be universally applicable, and potentially can bring more insights to future work. Next, we will inject a few assumptions to arrive at a more practical form for building the flow matching models.

An important property of Kronecker product: Given two invertible matrices ${\bm{A}}$ and ${\bm{B}}$, their Kronecker product satisfies $({\bm{A}}\otimes{\bm{B}})^{-1}={\bm{A}}^{-1}\otimes{\bm{B}}^{-1}$.
Using such a property, in the limit as $\nu\to 0$, we have

$$ $\Lambda_{i}\to{\bm{L}}_{i}\otimes({\bm{V}}_{i}^{\top}{\bm{V}}_{i})\quad\Longrightarrow\quad\Sigma_{i}={\bm{L}}_{i}^{-1}\otimes({\bm{V}}_{i}^{\top}{\bm{V}}_{i})^{-1}.$ (47) $$

According to Propensity [2](#Thmlemma2), the squared 2-Wasserstein distance between two Gaussian measures is given by

$$ $\mathcal{W}_{2}^{2}(\eta_{0},\eta_{1})=\underbrace{\|\bm{\mu}_{0}^{\prime}-\bm{\mu}_{1}^{\prime}\|^{2}}_{\text{Mean term}}+\underbrace{\operatorname{Tr}\Bigl{(}\Sigma_{0}+\Sigma_{1}-2\Bigl{(}\Sigma_{0}^{1/2}\Sigma_{1}\,\Sigma_{0}^{1/2}\Bigr{)}^{1/2}\Bigr{)}}_{\text{Covariance Term}}.$ (48) $$

###### Proposition 6 (Bures-Wasserstein Distance) .

#### Mean Term.

Since $\bm{\mu}_{i}^{\prime}={\bm{V}}\otimes\bm{\mu}_{i}$,
the mean difference becomes

$$ $\|\bm{\mu}_{0}^{\prime}-\bm{\mu}_{1}^{\prime}\|^{2}=\|{\bm{V}}_{0}\bm{\mu}_{0}-{\bm{V}}_{1}\bm{\mu}_{1}\|_{F}^{2}=\|{\bm{X}}_{0}-{\bm{X}}_{1}\|^{2}_{\text{F}}$ (49) $$

#### Covariance term.

Using the property of the Kronecker product, the square root of [Equation˜47](#A2.E47) factors in as

$$ $\Sigma_{i}^{1/2}={\bm{L}}_{i}^{-1/2}\otimes({\bm{V}}_{i}^{\top}{\bm{V}}_{i})^{-1/2}.$ (50) $$

and

$$ $\displaystyle\Sigma_{0}^{1/2}\Sigma_{1}\,\Sigma_{0}^{1/2}$ $\displaystyle=\Bigl{(}{\bm{L}}_{0}^{-1/2}{\bm{L}}_{1}^{-1}{\bm{L}}_{0}^{-1/2}\Bigr{)}\otimes\Bigl{(}({\bm{V}}_{0}^{\top}{\bm{V}}_{0})^{-1/2}({\bm{V}}_{1}^{\top}{\bm{V}}_{1})^{-1}({\bm{V}}_{0}^{\top}{\bm{V}}_{0})^{-1/2}\Bigr{)}$ (51) $$

We first look into the term related to ${\bm{V}}_{0}$ and ${\bm{V}}_{1}$, which is,

$$ $\displaystyle\operatorname{Tr}\Bigl{(}({\bm{V}}_{0}^{\top}{\bm{V}}_{0})^{-1/2}({\bm{V}}_{1}^{\top}{\bm{V}}_{1})^{-1}({\bm{V}}_{0}^{\top}{\bm{V}}_{0})^{-1/2}\Bigr{)}$ $\displaystyle=\operatorname{Tr}\Bigl{(}({\bm{V}}_{1}^{\top}{\bm{V}}_{1})^{-1}({\bm{V}}_{0}^{\top}{\bm{V}}_{0})^{-1/2}({\bm{V}}_{0}^{\top}{\bm{V}}_{0})^{-1/2}\Bigr{)}$ (52) $\displaystyle=\operatorname{Tr}\Bigl{(}({\bm{V}}_{1}^{\top}{\bm{V}}_{1})^{-1}({\bm{V}}_{0}^{\top}{\bm{V}}_{0})^{-1}\Bigr{)}$ $$

As $\operatorname{Tr}({\bm{A}}+{\bm{B}})=\operatorname{Tr}({\bm{A}})+\operatorname{Tr}({\bm{B}})$ the covariance term becomes

$$ Covariance Term (53) $\displaystyle=\operatorname{Tr}\Bigl{(}\Sigma_{0}+\Sigma_{1}-2\Bigl{(}\Sigma_{0}^{1/2}\Sigma_{1}\,\Sigma_{0}^{1/2}\Bigr{)}^{1/2}\Bigr{)}$ $\displaystyle=\operatorname{Tr}(\Sigma_{0})+\operatorname{Tr}(\Sigma_{1})-2\operatorname{Tr}\left((\Sigma_{0}^{1/2}\Sigma_{1}\,\Sigma_{0}^{1/2})^{1/2}\right)$ $\displaystyle=\operatorname{Tr}\Bigl{(}{\bm{L}}_{0}^{-1}\otimes({\bm{V}}_{0}^{\top}{\bm{V}}_{0})^{-1}+{\bm{L}}_{1}^{-1}\otimes({\bm{V}}_{1}^{\top}{\bm{V}}_{1})^{-1}-2\,\Bigl{(}{\bm{L}}_{0}^{-1/2}{\bm{L}}_{1}^{-1}{\bm{L}}_{0}^{-1/2}\Bigr{)}^{1/2}\otimes({\bm{V}}_{1}^{\top}{\bm{V}}_{1})^{-1/2}({\bm{V}}_{0}^{\top}{\bm{V}}_{0})^{-1/2}\Bigr{)}$ $$

Given that $\operatorname{Tr}({\bm{A}}\otimes{\bm{B}})=\operatorname{Tr}({\bm{A}})\operatorname{Tr}({\bm{B}})$ and $\operatorname{Tr}({\bm{V}}^{\top}{\bm{V}})=\|{\bm{V}}\|_{\text{F}}^{2}$ for any real-valued matrix ${\bm{V}}$, we can further derive,

$$ Covariance Term $\displaystyle=\operatorname{Tr}[({\bm{V}}_{0}^{\top}{\bm{V}}_{0})^{-1}]\operatorname{Tr}({\bm{L}}_{0}^{\dagger})+\operatorname{Tr}[({\bm{V}}_{1}^{\top}{\bm{V}}_{1})^{-1}]\operatorname{Tr}({\bm{L}}_{1}^{\dagger})$ (54) $\displaystyle-2\operatorname{Tr}\Bigl{(}{\bm{L}}_{0}^{\dagger/2}{\bm{L}}_{1}^{\dagger}{\bm{L}}_{0}^{\dagger/2}\Bigr{)}^{1/2}\cdot\operatorname{Tr}[({\bm{V}}_{1}^{\top}{\bm{V}}_{1})^{-1/2}({\bm{V}}_{0}^{\top}{\bm{V}}_{0})^{-1/2}].$ $$

Unfortunately, to simplify this equation, we have to make the two gram matrix, $({\bm{V}}_{0}^{\top}{\bm{V}}_{0})^{-1}$ and $({\bm{V}}_{1}^{\top}{\bm{V}}_{1})^{-1}$ agree, i.e., $({\bm{V}}_{1}^{\top}{\bm{V}}_{1})^{-1}=({\bm{V}}_{0}^{\top}{\bm{V}}_{0})^{-1}$. This will be satisfied if and only if there exists an orthogonal matrix ${\bm{Q}}$ such that

$$ ${\bm{V}}_{1}^{\dagger}={\bm{V}}_{0}^{\dagger}{\bm{Q}}.$ $$

Thus, to further process, we simply consider the case when ${\bm{V}}_{1}$ and ${\bm{V}}_{0}$ are exactly the same, i.e., ${\bm{V}}_{1}={\bm{V}}_{0}={\bm{V}}$ (we have already discussed how realistic this assumption is in [Section˜A.3](#A1.SS3)). So that we work under the assumptions that $\|{\bm{V}}_{0}^{\dagger}\|_{F}^{2}=\|{\bm{V}}_{1}^{\dagger}\|_{F}^{2}=\beta$, which simplify the trace as

$$ Covariance Term $\displaystyle=\beta\cdot\operatorname{Tr}\Bigl{(}{\bm{L}}_{0}^{\dagger}+{\bm{L}}_{1}^{\dagger}-2\,\Bigl{(}{\bm{L}}_{0}^{\dagger/2}{\bm{L}}_{1}^{\dagger}{\bm{L}}_{0}^{\dagger/2}\Bigr{)}^{1/2}\Bigr{)}.$ (55) $$

Combining the mean term and the covariance term, we obtain the Wasserstein distance of ${\mathcal{W}}_{c}(\eta_{{\mathcal{X}}_{0}},\eta_{{\mathcal{X}}_{1}})$

For calculating ${\mathcal{W}}_{c}(\eta_{{\mathcal{E}}_{0}},\eta_{{\mathcal{E}}_{1}})$, we have the freedom to choose the cost function when obtaining the Wasserstein distance. Note that ${\bm{W}}$ serves as the prior for the Gaussian covariance matrix $\Sigma$, where the covariance has to be positive-semi definite. Thus, according to (Bhatia et al., 2019), a proper distance between two positive semi-definite matrices is measured by

$$ ${\mathcal{W}}(\eta_{{\mathcal{E}}_{0}},\eta_{{\mathcal{E}}_{1}})=\left\|\Sigma_{0}^{1/2}-\Sigma_{1}^{1/2}\right\|_{F}^{2}.$ (56) $$

Coincidently, this is another usage case when the Bures-Wasserstein metric is utilized. Putting everything together, the Wasserstein distance in the limit $\nu\to 0$ is

$$ $\displaystyle d_{\text{BW}}({\mathcal{G}}_{0},{\mathcal{G}}_{1})$ $\displaystyle=\|{\bm{V}}_{0}\bm{\mu}_{0}-{\bm{V}}_{1}\bm{\mu}_{1}\|_{F}^{2}+(\beta+1)\cdot\operatorname{Tr}\Bigl{(}{\bm{L}}_{0}^{\dagger}+{\bm{L}}_{1}^{\dagger}-2\,\Bigl{(}{\bm{L}}_{0}^{\dagger/2}{\bm{L}}_{1}^{\dagger}{\bm{L}}_{0}^{\dagger/2}\Bigr{)}^{1/2}\Bigr{)}.$ (57) $\displaystyle=\underbrace{\|{\bm{X}}_{0}-{\bm{X}}_{1}\|_{F}^{2}}_{d_{\bm{X}}({\bm{X}}_{0},{\bm{X}}_{1})}+(\beta+1)\underbrace{\cdot\operatorname{Tr}\Bigl{(}{\bm{L}}_{0}^{\dagger}+{\bm{L}}_{1}^{\dagger}-2\,\Bigl{(}{\bm{L}}_{0}^{\dagger/2}{\bm{L}}_{1}^{\dagger}{\bm{L}}_{0}^{\dagger/2}\Bigr{)}^{1/2}\Bigr{)}}_{{}_{d_{\bm{L}}({\bm{L}}_{0},{\bm{L}}_{1})}}.$ $$

This expression separates the contribution of the mean difference (transformed by ${\bm{V}}$) and the discrepancy between the covariance structures (encoded in ${\bm{L}}_{0}$ and ${\bm{L}}_{1}$). This could be further used to derive BW interpolation, which we will show in [Section˜C.1](#A3.SS1). In the main body, constant $\beta$ actually corresponds to $\beta+1$ here. This complete our derivation in [˜1](#Thmprop1).

## Appendix C Derivation of Bures-Wasserstein flow matching

In order to build the flow matching framework, we need to derive the optimal interpolation and the corresponding velocities for the probability path $p(G_{t}\mid G_{0},G_{1})$. This is achieved via the OT displacement between two graph distributions.

### C.1 The Bures-Wasserstein graph interpolation

We aim to recover the proposition stated as follows.

###### Proposition 7 (Bures-Wasserstein interpolation) .

The interpolation is an extension of the concept of mean, where in the optimal transport world, the Wasserstein barycenter (mean) of measures $\eta_{0},\ldots\eta_{m-1}$ under weights $\lambda_{0},\ldots\lambda_{m-1}$ can be derived over the following optimization problem:

$$ $\bar{\eta}=\underset{\eta}{\arg\min}\sum_{j=0}^{m-1}\lambda_{j}\left({\mathcal{W}}_{2}\left(\eta,\eta_{j}\right)\right)^{2}$ (59) $$

When $m=2$, based on the Bures-Wasserstein (BW) distance, we can define the OT displacement minimization problem on graphs described as,

$$ ${\mathcal{G}}_{t}=\operatorname*{arg\,min}_{\tilde{{\mathcal{G}}}}~(1-t)d_{\text{BW}}({\mathcal{G}}_{0},\tilde{{\mathcal{G}}})+td_{\text{BW}}(\tilde{{\mathcal{G}}},{\mathcal{G}}_{1}).$ (60) $$

where $d_{\text{BW}}({\mathcal{G}}_{0},{\mathcal{G}}_{1})$ is described in [˜1](#Thmprop1). The optimal graph interpolation is the solution to the problem.

In the setting of graph, this becomes a two-variable optimization problem, where

$$ ${\mathcal{X}}_{t},{\mathcal{E}}_{t}=\operatorname*{arg\,min}_{\tilde{{\mathcal{X}}},\tilde{{\mathcal{E}}}}~(1-t)d_{\text{BW}}({\mathcal{G}}_{0},\tilde{{\mathcal{G}}})+td_{\text{BW}}(\tilde{{\mathcal{G}}},{\mathcal{G}}_{1}).$ (61) $$

Fortunately, recall in [Equation˜57](#A2.E57) that our distance measurement $d_{\text{BW}}({\mathcal{G}}_{0},{\mathcal{G}}_{1})$ is decomposed into $d_{{\bm{X}}}({\bm{X}}_{0},{\bm{X}}_{1})$ and $d_{{\bm{L}}}({\bm{L}}_{0},{\bm{L}}_{1})$, then the optimization over node and edges are disentangleable into solving the two sub optimization problem,

$$ $\displaystyle\text{Sub-question 1: }\quad\bar{{\bm{X}}_{t}}=\underset{\tilde{{\bm{X}}}}{\arg\min}\quad(1-t)\|{\bm{X}}_{0}-\tilde{{\bm{X}}}\|_{F}^{2}+t\|\tilde{{\bm{X}}}-{\bm{X}}_{1}\|_{F}^{2}$ (62) $\displaystyle\text{Sub-question 2: }\quad\bar{{\bm{L}}_{t}}=\underset{\tilde{{\bm{L}}}}{\arg\min}\quad(1-t)d_{\bm{L}}({\bm{L}}_{0},\tilde{{\bm{L}}})+td_{\bm{L}}({\bm{L}}_{1},\tilde{{\bm{L}}})$ $$

This two problems are completely disentangled thus we can solve them separately.

###### Proposition 7 (Bures-Wasserstein interpolation) .

#### Sub-question 1

For the first problem, we simply set the derivate to 0 and get,

$$ $(1-t)(\tilde{{\bm{X}}}-{\bm{X}}_{0})+t(\tilde{{\bm{X}}}-{\bm{X}}_{1})=0\rightarrow{\bm{X}}_{t}=(1-t){\bm{X}}_{0}+t({\bm{X}}_{1})$ (63) $$

#### Subquestion 2

The second subproblem is equivalent in deriving the covariance of Bures-Wasserstein interpolation between two Gaussian measures, $\eta_{0}\sim\mathcal{N}\left(0,{\bm{L}}_{0}^{\dagger}\right)$ and $\eta_{1}\sim\mathcal{N}\left(0,{\bm{L}}_{1}^{\dagger}\right)$. This problem has been properly addressed in Haasler & Frossard (2024) and here we just verbose their results. For more details we refer the reader to  Haasler & Frossard (2024) for a further discussion.

The optimal transport geodesic between $\eta_{0}\sim\mathcal{N}\left(0,{\bm{L}}_{0}^{\dagger}\right)$ and $\eta_{1}\sim\mathcal{N}\left(0,{\bm{L}}_{1}^{\dagger}\right)$ is defined by $\eta_{t}=((1-t)I+t{\bm{T}})_{\#}\eta_{0}$, where the symbol “$\#$” denotes the push-forward of a measure by a mapping, ${\bm{T}}$ is a linear map that satisfies ${\bm{T}}{\bm{L}}_{0}^{\dagger}{\bm{T}}={\bm{L}}_{1}^{\dagger}$.

We define a new matrix ${\bm{M}}$ and do normalization, which leads to,

$$ ${\bm{T}}={\bm{L}}_{0}^{1/2}{\bm{M}}{\bm{L}}_{0}^{1/2}$ (64) $$

Plug in gives,

$$ $\displaystyle{\bm{T}}{\bm{L}}_{0}^{\dagger}{\bm{T}}^{\top}$ $\displaystyle={\bm{L}}_{0}^{1/2}{\bm{M}}{\bm{L}}_{0}^{1/2}{\bm{L}}_{0}^{\dagger}\left({\bm{L}}_{0}^{1/2}{\bm{M}}{\bm{L}}_{0}^{1/2}\right)^{\top}$ (65) $\displaystyle={\bm{L}}_{0}^{1/2}{\bm{M}}{\bm{M}}^{\top}{\bm{L}}_{0}^{1/2}.$ $$

So that we obtain

$$ ${\bm{L}}_{1}^{\dagger}={\bm{L}}_{0}^{1/2}{\bm{M}}{\bm{M}}^{\top}{\bm{L}}_{0}^{1/2}\rightarrow{\bm{M}}=({\bm{L}}_{0}^{\dagger/2}{\bm{L}}_{1}^{\dagger}{\bm{L}}_{0}^{\dagger/2})^{1/2}$ (66) $$

Replace ${\bm{T}}$ and we get,

$$ ${\bm{T}}={\bm{L}}_{0}^{1/2}\left({\bm{L}}_{0}^{\dagger/2}{\bm{L}}_{1}^{\dagger}{\bm{L}}_{0}^{\dagger/2}\right)^{1/2}{\bm{L}}_{0}^{1/2}$ (67) $$

Given that the geodesic $\eta_{t}=((1-t){\bm{I}}+t{\bm{T}})_{\#}\eta_{0}$ which also has a Gaussian form $\eta_{t}\sim\mathcal{N}\left(0,\Sigma_{t}\right)$, We can then write the covariance matrix and obtain

$$ $\displaystyle{\bm{L}}_{t}^{\dagger}$ $\displaystyle=\Sigma_{t}=((1-t)I+t{\bm{T}}){\bm{L}}_{0}^{\dagger}((1-t)I+t{\bm{T}})$ (68) $\displaystyle={\bm{L}}_{0}^{1/2}\left((1-t){\bm{L}}_{0}^{\dagger}+t\left({\bm{L}}_{0}^{\dagger/2}{\bm{L}}_{1}^{\dagger}{\bm{L}}_{0}^{\dagger/2}\right)^{1/2}\right){\bm{L}}_{0}^{1/2}{\bm{L}}_{0}^{\dagger}{\bm{L}}_{0}^{1/2}\left((1-t){\bm{L}}_{0}^{\dagger}+t\left({\bm{L}}_{0}^{\dagger/2}{\bm{L}}_{1}^{\dagger}{\bm{L}}_{0}^{\dagger/2}\right)^{1/2}\right){\bm{L}}_{0}^{1/2}$ $\displaystyle={\bm{L}}_{0}^{1/2}\left((1-t){\bm{L}}_{0}^{\dagger}+t\left({\bm{L}}_{0}^{\dagger/2}{\bm{L}}_{1}^{\dagger}{\bm{L}}_{0}^{\dagger/2}\right)^{1/2}\right)^{2}{\bm{L}}_{0}^{1/2}$ $$

Which ends the derivation.

### C.2 Deriving the velocity of BW interpolation

We first show the general form of the velocity term for the Gaussian and Dirac measures.

#### Gaussian Measure.

For a time-parametrized Gaussian density $p_{t}(x)=\mathcal{N}\left(x;\bm{\mu}_{t},\Sigma_{t}\right)$,
the velocity field $v_{t}(x)$ satisfies the continuity equation

$$ $\partial_{t}p_{t}+\nabla\cdot\left(p_{t}v_{t}\right)=0,$ $$

is an affine function of $x$. And the instantaneous velocity field follows,

$$ $v_{t}({\mathcal{X}})=\dot{\bm{\mu}}_{t}+\frac{1}{2}\dot{\Sigma}_{t}\Sigma_{t}^{-1}\left({\mathcal{X}}-\bm{\mu}_{t}\right).$ $$

#### Dirac Measure.

When the measure is a Dirac function,

$$ $p_{t}(x)=\delta\left(\cdot,\bm{\mu}_{t}\right).$ $$

We can just consider it as the limited case of the Gaussian measure, when $\Sigma_{t}\rightarrow 0$. So that the velocity at simply takes

$$ $v_{t}(x)=\dot{\bm{\mu}}_{t}.$ $$

We then move to prove the following propensity for Bures-wasserstein velocity.

###### Proposition 8 (Bures-Wasserstein velocity) .

Proof:

###### Proposition 8 (Bures-Wasserstein velocity) .

#### The graph structure velocity.

As we assume the edges, $E_{t}\sim\delta(\cdot,{\bm{W}}_{t}),$ following a dirac distribution, the velocity is defined as

$$ $v_{t}(E_{t})=\dot{{\bm{W}}}_{t}.$ $$

Given that, $\dot{{\bm{W}}}_{t}=\text{diag}(\dot{{\bm{L}}_{t}})-\dot{{\bm{L}}_{t}}$, we transit fo deriving the derivative of the Laplacian matrix, $\dot{{\bm{L}}_{t}}$.
Using the fact that,

$$ $\frac{d}{dt}\left({\bm{A}}^{-1}\right)=-{\bm{A}}^{-1}\frac{d{\bm{A}}}{dt}{\bm{A}}^{-1}$ $$

we obtain the derivate of Laplacian matrix,

$$ $\dot{{\bm{L}}_{t}}=\frac{d(\Sigma_{t}^{\dagger})}{dt}=\Sigma_{t}^{\dagger}\frac{d\Sigma_{t}}{dt}\Sigma_{t}^{\dagger}={\bm{L}}_{t}\frac{d\Sigma_{t}}{dt}{\bm{L}}_{t}$ (70) $$

According to [Equation˜68](#A3.E68) and  [Equation˜67](#A3.E67), the covariance matrix is defined through the interpolation,

$$ $\Sigma_{t}=((1-t){\bm{I}}+t{\bm{T}}){\bm{L}}_{0}^{\dagger}((1-t){\bm{I}}+t{\bm{T}}):={\bm{R}}_{t}{\bm{L}}_{0}^{\dagger}{\bm{R}}_{t}$ (71) $$

where ${\bm{R}}_{t}=(1-t){\bm{I}}+t{\bm{T}}$. Taking the derivative, we get,

$$ $\dot{\Sigma}_{t}=\frac{d}{dt}\left({\bm{R}}_{t}\Sigma_{0}{\bm{R}}_{t}\right)={\bm{R}}_{t}^{\prime}\Sigma_{0}{\bm{R}}_{t}+{\bm{R}}_{t}\Sigma_{0}{\bm{R}}_{t}^{\prime}=({\bm{T}}-{\bm{I}})\Sigma_{0}{\bm{R}}_{t}+{\bm{R}}_{t}\Sigma_{0}({\bm{T}}-{\bm{I}})$ (72) $$

Using the fact that $\Sigma_{0}{\bm{R}}_{t}={\bm{R}}_{t}\Sigma_{0}=\Sigma_{t}$,
we obtain the covariance gradient

$$ $\dot{\Sigma}_{t}=({\bm{T}}-{\bm{I}})\Sigma_{t}+\Sigma_{t}({\bm{T}}-{\bm{I}})$ (73) $$

So that,

$$ $\displaystyle-\dot{{\bm{L}}_{t}}$ $\displaystyle=\frac{d(\Sigma_{t}^{\dagger})}{dt}=\Sigma_{t}^{\dagger}\frac{d\Sigma_{t}}{dt}\Sigma_{t}^{\dagger}={\bm{L}}_{t}\frac{d\Sigma_{t}}{dt}{\bm{L}}_{t}$ (74) $\displaystyle={\bm{L}}_{t}(({\bm{T}}-{\bm{I}}){\bm{L}}^{\dagger}_{t}+{\bm{L}}^{\dagger}_{t}({\bm{T}}-{\bm{I}})){\bm{L}}_{t}$ $\displaystyle={\bm{L}}_{t}({\bm{T}}-{\bm{I}})+({\bm{T}}-{\bm{I}}){\bm{L}}_{t}$ $\displaystyle={\bm{L}}_{t}{\bm{T}}+{\bm{T}}{\bm{L}}_{t}-2{\bm{L}}_{t}$ $$

Thus, $\dot{{\bm{L}}_{t}}=2{\bm{L}}_{t}-{\bm{L}}_{t}{\bm{T}}-{\bm{T}}{\bm{L}}_{t}$.

Given that ${\bm{L}}_{t}={\bm{D}}_{t}-{\bm{W}}_{t}$ so that ${\bm{W}}_{t}=\text{diag}({\bm{L}}_{t})-{\bm{L}}_{t}$, taking the derivative gives $\dot{{\bm{W}}_{t}}=\text{diag}(\dot{{\bm{L}}_{t}})-\dot{{\bm{L}}_{t}}$. As we assume the edges, $E_{t}\sim\delta(\cdot,{\bm{W}}_{t})$, the derivate directly yields the velocity,

$$ $v_{t}(E_{t}\mid G_{0},G_{1})=\dot{{\bm{W}}_{t}}=\text{diag}(\dot{{\bm{L}}_{t}})-\dot{{\bm{L}}_{t}}.$ $$

#### The node feature velocity.

The instantaneous velocity field follows,

$$ $v_{t}({\mathcal{X}}\mid G_{0},G_{1})=\dot{\bm{\mu}}_{t}+\frac{1}{2}\dot{\Sigma}_{t}\Sigma_{t}^{-1}\left({\mathcal{X}}-\bm{\mu}_{t}\right).$ $$

The mean gradient interpolating $\eta_{0}$ and $\eta_{1}$ can be written as $\dot{\bm{\mu}}_{t}={\bm{X}}_{1}-{\bm{X}}_{0}$ and ${\bm{X}}_{t}=(1-t){\bm{X}}_{0}+t{\bm{X}}_{1}$. So that the velocity leads to,

$$ $v_{t}({\mathcal{X}}\mid G_{0},G_{1})={\bm{X}}_{1}-{\bm{X}}_{0}+\frac{1}{2}\dot{{\bm{L}}^{\dagger}_{t}}{\bm{L}}_{t}\left({\mathcal{X}}-{\bm{X}}_{t}\right).$ $$

However, in practice, we do not need such a complicated velocity term. We wish to avoid the estimation of complex gradient-inverse term so that we can escape from the complicated computation. Under the assumption that the amplitude of covariance is much smaller than the mean difference, we can omit the second term and just keep the mean difference. Hence the instantaneous velocity is simply described as

$$ $v_{t}(X_{t}\mid G_{0},G_{1})={\bm{X}}_{1}-{\bm{X}}_{0}=X_{1}-X_{0}=\frac{1}{1-t}({\bm{X}}_{1}-{\bm{X}}_{t})$ (75) $$

Figure: Algorithm 1 BWFlow Training

## Appendix D Discrete Bures-Wasserstein flow matching for graph generation

### D.1 Probability path construction for discrete Bures-Wasserstein flow matching

#### The discrete probability path.

We design the probability path as discrete distributions,

$$ $\displaystyle p_{t}(x_{v}\mid G_{0},G_{1})=\operatorname{Categorical}([{\bm{X}}_{t}]_{v}),\quad p_{t}(e_{uv}\mid G_{0},G_{1})=\operatorname{Bernoulli}([{\bm{W}}_{t}]_{uv})$ (76) $\displaystyle\text{ s.t. }p_{0}({\mathcal{G}})=\delta(G_{0},\cdot),p_{1}({\mathcal{G}})=\delta(G_{1},\cdot)$ $$

where ${\bm{W}}_{t}={\bm{D}}_{t}-{\bm{L}}_{t}$ with ${\bm{X}}_{t}$ and ${\bm{L}}_{t}$ defined the same in [Equation˜11](#S3.E11). We consider the fact that the Dirac distribution is a special case when the Categorical/Bernoulli distribution has probability 1 or 0, so the boundary condition $p_{0}({\mathcal{G}})=\delta(G_{0},\cdot),p_{1}({\mathcal{G}})=\delta(G_{1},\cdot)$ holds. As such, ${\bm{X}}_{t}=(1-t){\bm{X}}_{0}+t{\bm{X}}_{1}\in[0,1]^{|{\mathcal{V}}|\times K}$. Since the boundary condition for each entry, $[{\bm{X}}_{0}]_{v}$ and $[{\bm{X}}_{1}]_{v}$ are two one-hot embeddings, $[{\bm{X}}_{t}]_{v}=t[{\bm{X}}_{0}]_{v}+(1-t)[{\bm{X}}_{1}]_{v}$ would sum to one, which works as a valid probability vector. Thus, $\operatorname{Categorical}([{\bm{X}}_{t}]_{v})$ is a K-class categorical distribution.

For the edge distribution, we just consider $e_{uv}$ is conditionally independent of the other given $[{\bm{W}}_{t}]_{uv}$. One thing to emphasize is that, given the nature of Bures-Wasserstein interpolation, the yielded ${\bm{W}}_{t}$ is not always bounded by $[0,1]$ thus we have to hard-clip the boundary.

### D.2 Approximating Wasserstein distance in Bernoulli distributions

To make sure that the individual nodes are structured and developed jointly while doing flow matching, we assume that the $\operatorname{vec}({\mathcal{X}})$ still maintains a covariance matrix similar to  [Equation˜8](#S3.E8), which gives $\bm{\Lambda}=(\nu{\bm{I}}+{\bm{L}})\otimes\bm{V}^{\top}\bm{V}$ given that ${\mathcal{X}}$ is emitted from a latent variable ${\mathcal{H}}$ through an affine transformation and the latent variable has a covariance matrix $(\nu{\bm{I}}+{\bm{L}})^{-1}$. Different from the Gaussian case, the latent variable would still be a discrete distribution, so that the affine transformation carries the covariance matrix out.

Unfortunately, the Wasserstein distance between two discrete graph distributions that follow [Equation˜13](#S3.E13) does not have a closed-form solution given the complex interwined nature. However, it is possible to use the central limit theorem applied to ${\mathcal{X}}$ so that we can approximate the Wasserstein distance of two Bernoulli distributions with the Gaussian counterpart. This approximation works when we are in high-dimensional case (high dimension means $|{\mathcal{V}}|d$ is moderately large.), and the OT-distance between two such Bernoulli distributions is well-captured by the corresponding Gaussian formula, which we already introduced in [Equation˜57](#A2.E57).

With such nature, even though we are not sampling from Gaussian distributions anymore, it is possible to approximate the Wasserstein distance between two multivariate discrete distributions with the Gaussian counterpart, so the conclusions, such as optimal transport displacements, still hold. And we can similarly derive the Bures-Wasserstein velocity as in the next section.

### D.3 Velocity for discrete Bures-Wasserstein flow matching

#### Node Velocity.

For node-wise, the path of node features ${\mathcal{X}}_{t}$ can be re-written as $p_{t}({\mathcal{X}})=(1-t)\delta(\cdot,{\bm{X}}_{0})+t\delta(\cdot,{\bm{X}}_{1})$ so the conditional velocity can be accessed through $v_{t}(X_{t}\mid G_{0},G_{1})=[\delta(\cdot,{\bm{X}}_{1})-\delta(\cdot,{\bm{X}}_{t})]/(1-t)$ similar as the derivation in (Gat et al., 2024).

#### Edge Velocity.

For edge-wise, we look into each entry of the adjacency matrix ${\bm{W}}$, and consider a time-dependent Bernoulli distribution, the probability density function is:

$$ $p_{t}(e_{uv})=[{\bm{W}}_{t}]_{uv}^{e_{uv}}(1-[{\bm{W}}_{t}]_{uv})^{1-e_{uv}},\qquad e_{uv}\in\{0,1\}.$ (77) $$

To properly define a velocity $v(x,t)$, it should follow the continuity equation

$$ $\frac{\partial}{\partial t}p_{t}(e_{uv})+\nabla\cdot(p\,v)_{t}(e_{uv})\;=\;0.$ (78) $$

We use $x$ and $y$ to denote two states of $e_{uv}$ ($p(e_{uv}=x):=p(x),p(e_{uv}=y):=p(y)$), then the divergence term is

$$ $\nabla\cdot(p\,v)(e_{uv}=x)=\sum_{y\neq x}[p_{t}(y)\,v_{t}(y\to x)-p_{t}(x)\,v_{t}(x\to y)].$ (79) $$

As we are working on a Bernoulli distribution, then the forward equations become

$$ $\begin{cases}\partial_{t}p(0)=p(1)\,v_{t}(1\to 0)-p(0)\,v_{t}(0\to 1),\\ \partial_{t}p(1)=p(0)\,v_{t}(0\to 1)-p(1)\,v_{t}(1\to 0).\end{cases}$ (80) $$

Since $p_{t}(1)=[{\bm{W}}_{t}]_{uv},$, we have $\partial_{t}p(1)=[\dot{{\bm{W}}_{t}}]_{uv}$ and
$\partial_{t}p(0)=-[\dot{{\bm{W}}_{t}}]_{uv}$.
Hence

$$ $p(0)\,v_{t}(0\to 1)\;-\;p(1)\,v_{t}(1\to 0)=[\dot{{\bm{W}}_{t}}]_{uv}.$ $$

There are many solutions to the above equation. We chose a symmetric solution so that the transition of $e_{uv}\to 1-e_{uv}$ with

$$ $v_{t}(0\to 1)=\frac{[\dot{{\bm{W}}_{t}}]_{uv}}{1-[{\bm{W}}_{t}]_{uv}},\quad v_{t}(1\to 0)=-\frac{[\dot{{\bm{W}}_{t}}]_{uv}}{[{\bm{W}}_{t}]_{uv}}.$ $$

Finally for concise, we can write write it as a velocity field on states $e_{uv}\in\{0,1\}$, note $1-2e_{uv}$ is $+1$ at $e_{uv}=0$ and $-1$ at $e_{uv}=1$. Thus, we have

$$ $v(e_{uv},t)=(1-2e_{uv})\,\frac{[\dot{{\bm{W}}}_{t}]_{uv}}{[{\bm{W}}_{t}]_{uv}\,(1-{\bm{W}}_{t}]_{uv})},e_{uv}\in\{0,1\},$ $$

which in matrix form gives

$$ $v_{t}(E_{t}\mid G_{1},G_{0})=\left(1-2E_{t}\right)\frac{\dot{{\bm{W}}_{t}}}{{\bm{W}}_{t}\circ\left(1-{\bm{W}}_{t}\right)}.$ (81) $$

Combine the node velocity and the edge velocity, we can now introduce the Discrete Bures-Wasserstein Flow matching algorithm, with the training and inference part respectively introduced in [Algorithm˜3](#alg3) and [Algorithm˜4](#alg4).

Figure: Algorithm 3 Discrete BWFlow Training

## Appendix E Design space for Bures Wasserstein interpolation and velocity

In the introduction part, we have already compared different probability paths and how they are impacting the inference time sampling. While the Bures-Wasserstein flow path is shown to produce a better probability path for the model to learn, as we illustrated in [Figure˜1(a)](#S1.F1.sf1), we have to point out that linear interpolation and the corresponding probability path can still converge to the data distribution with sufficiently large flow steps. As if we conduct sampling with infinite flow steps during the later stage of flow, the samples are still able to arrive at the target distributions. A similar pattern exists in diffusion models when they are considered as a Monte-Carlo Markov Chain, and they need sufficiently large steps to converge. We emphasize that the convergence gap in [Figure˜1(c)](#S1.F1.sf3) would be slowly recovered as the number of flow steps increases.

Given that different sampling algorithms can all bring the samples to the data distributions under certain conditions, we wish to understand the huge design space of Bures Wasserstein interpolation. We list the advantages and disadvantages in different techniques and discuss further when each techniques should be used.

In general, we consider two important steps to construct the flow matching for graph generation, specifically, training and sampling. In training, the main challenge is to obtain a valid real velocity $u(G_{t})$ to be regressed to, so we listed a few strategies that can help us with that. In sampling, the challenge becomes how to reconstruct the probability path through the velocity estimated.

### E.1 The Training Design

In general, the learning objective in flow matching depends on regressing the velocity term. There are several way to obtain the velocity.

- 1.
Exact velocity estimation. Use [Equation˜3](#S2.E3) as the parameterization and learn $p_{\theta}({\mathcal{G}}_{1}\mid G_{t})$
- 2.
Numerical Approximation. In the implementation of (Stärk et al., 2024), the derivative is calculated through numerical approximation. To achieve better efficiency in calculating velocity, we simply consider a numerical estimation as in (Stärk et al., 2024), where the velocity term is obtained as, $\dot{{\bm{L}}_{t}}=({\bm{L}}_{t+\Delta t}-{\bm{L}}_{t})/\Delta t$. Regressing on the numerical difference can provide an estimation for the velocity.
- 3.
AutoDiff. In (Chen & Lipman, 2024), the derivative of the probability path is evaluated through Pytorch AutoDiff. However, in practice we find this method unstable.

We summarized the training stage model parameterization in [Table˜4](#A5.T4)

**Table 4: The model parameterization for flow matching in training stage**
|  | Continuous Flow Matching | Discrete Flow Matching |
| --- | --- | --- |
| $x$-prediction | $v^{\theta}_{t}\left(G_{t}\right)=\frac{1}{1-t}\left[\tilde{G}^{\theta}_{t}(G_{t})-G_{t}\right]$ | $v^{\theta}_{t}\left(G_{t}\right)=\frac{1}{1-t}\left[p_{1\mid t}^{\theta}\left({\mathcal{G}}_{1}\mid G_{t}\right)-\delta\left(\cdot,G_{t}\right)\right]$ |
| Numerical Approximation | $v^{\theta}_{t}(G_{t})\approx G_{t+dt}-G_{t}$ | $v^{\theta}_{t}(G_{t})\approx p(G_{t+dt}\mid G_{0},G_{1})-p(G_{t}\mid G_{0},G_{1})$ |
| AutoDiff | $v^{\theta}_{t}(G_{t})\approx\dot{G_{t}}$ | Discrete velocity introduced in [Equation˜14](#S3.E14) |

### E.2 The Sampling Design

As we described in the [Equation˜3](#S2.E3), in our training framework, we actually train a denoised $p_{\theta}({\mathcal{G}}_{1}\mid{\mathcal{G}}_{t})$. With such a parameterization and taking discrete flow matching as an example, the sampling can be done through one of the following design choices:

- 1.
Target Guided Velocity Sampling. The velocity is designed as,
$v_{\theta}(G_{t})=\frac{1}{1-t}(p_{\theta}({\mathcal{G}}_{1}\mid G_{t})-\delta(G_{t},\cdot)).$
This design directly moves the current point $G_{t}$ towards the direction pointing to the predicted $G_{1}$. The target-guided velocity is guaranteed to converge to the data distribution, but the interpolant might lie outside the valid graph domain.
- 2.
BW velocity sampling. We use [Equation˜14](#S3.E14) to directly estimate the velocity and flow the Bures-Wasserstein probability path to generate new data points. This path is smooth in the sense of graph domain. However, this path requires more computational cost.
- 3.
Probability Path Reconstruction. The third option is directly reconstructing the probability path, i.e., we first obtain an estimated point,
$\tilde{G_{1}}\sim p_{\theta}({\mathcal{G}}_{1}\mid G_{t})$
and then construct the data point at $t+dt$, which gives
$G_{t+dt}\sim p({\mathcal{G}}_{t}\mid\tilde{G_{1}},G_{0})$
through LABEL:eq:bw_velocity. This is the most computationally costly method, which is obtained through the diffusion models. But this method also provide accurate probability path reconstruction.

In [Section˜4](#S4), we show BW velocity follows a path that minimizes the Wasserstein distance thus provides better performance, but sampling following linear velocity also provides convergence with much lower computational cost. So it is a trade-off to be considered in the real-world application.

**Table 5: Reconstructing probability path choices in flow matching during inference**
|  | Continuous Flow Matching | Discrete Flow Matching |
| --- | --- | --- |
| BW Velocity | LABEL:eq:bw_velocity | [Equation˜14](#S3.E14) |
| Target-guided Velocity | $v_{\theta}(G_{t})=\frac{1}{1-t}(\tilde{G_{1}}-G_{t})$ | $v_{\theta}(G_{t})=\frac{1}{1-t}(p_{\theta}({\mathcal{G}}_{1}\mid G_{t})-\delta(G_{t},\cdot))$ |
| Path Reconstruction | $G_{t+dt}=\delta(G_{t+dt}\mid\tilde{G_{1}},G_{0})$ | $G_{t+dt}\sim p(G_{t+dt}\mid\tilde{G_{1}},G_{0})$ |

## Appendix F Discussion and Limitations

### F.1 The implicit manipulation of probability path

Figure: (a) Linear training path
Refer to caption: /html/2506.14020/assets/x6.png

Though not explicitly mentioned, Qin et al. (2024) makes huge efforts to manipulate the probability path for better velocity estimation by extensively searching the design space, and their finding aligns well with the statement that the velocity should be smooth and consistently directing to the data points: 1) Time distortion: (The oragne line in [Figure˜5(b)](#A6.F5.sf2)) the polynomial distortion of training and sampling focus on the later stage of the probability trajectory, providing better velocity estimation in this area. This uneven sampling strategy is equivalent to pushing the probability path left to make it smooth. 2) Target guidance: (The orange line in [Figure˜5(a)](#A6.F5.sf1)) the target guidance directly estimate the direction from a point along the path towards the termination graph, so that the manipulated probability could smoothly pointing to the data distribution. and 3) Stochasticity injection: (The green line in [Figure˜5(a)](#A6.F5.sf1)) Stochasticity explores the points aside from the path, which avoid the path to be stuck in the platform area.

### F.2 Potential extension to diffusion models

In order to extend the flow matching algorithms with diffusion models, one important thing is to convert the pair-conditioned probability path and velocity to single boundary conditions. For instance, the probability path in flow matching has the form $p(G_{t}\mid G_{1},G_{0})$ and the velocity follows $v(G_{t}\mid G_{1},G_{0})$. As suggested in (Siraudin et al., 2024; Campbell et al., 2022; Xu et al., 2024), the discrete graph diffusion models require a velocity (which is equivalent to a ratio matrix) to perturb the data distribution conditioned on the data points, which we denote as $v(G_{t}\mid G_{1})$. As long as the unilateral conditional velocity has a tractable form, one can first sample a $G_{1}$ and get $G_{t}$ through iteratively doing to:

$$ $G_{t-dt}=G_{t}-v(G_{t}\mid G_{1})dt$ $$

starting from $G_{1}$. So that one can easily construct the probability path $p(G_{t}\mid G_{1})$ to fit into the diffusion model framework. In practice, given that we know the explicit form of $v(G_{t}\mid G_{1},G_{t^{\prime}})$ (just replace $G_{0}$ in the expression), the unilateral conditional velocity can be obtained through taking the limitation,

$$ $v(G_{t}\mid G_{1})=v(G_{t}\mid G_{1},G_{t})=\lim_{t^{\prime}\to t}v(G_{t}\mid G_{1},G_{t^{\prime}}).$ $$

Both linear interpolation and our Bures-Wasserstein interpolation can achieve this easily. We just provide a discussion here and will leave this as future work as this paper does not focus on diffusion models but on flow matching models.

### F.3 Permutation invariance

The Bures-Wasserstein distance between two graph distributions is not permutation invariant, and the minimal value is obtained through the graph alignment. So ideally, to achieve optimal transport, graph alignment and mini-batch matching could provide a better probability path. However, permutation invariance is not always a desired property since we only want to find a path that better transforms from the reference distribution to the data distributions. As an illustration, the widely used linear interpolation to construct graph flow (Qin et al., 2024) does not guarantee permutation invariance as well. And it is proved that, if the measurement is based on Wasserstein distance between two Gaussian distributions.

$$ $\displaystyle d_{\text{BW}}(\eta_{{\mathcal{G}}_{0}},\eta_{{\mathcal{G}}_{1}})\leq d_{\text{Arithmetic}}(\eta_{{\mathcal{G}}_{0}},\eta_{{\mathcal{G}}_{1}})$ (82) with $\displaystyle d_{\text{BW}}(\eta_{{\mathcal{G}}_{0}},\eta_{{\mathcal{G}}_{1}})=\left\|{\bm{X}}_{0}-{\bm{X}}_{1}\right\|_{F}^{2}+\beta\operatorname{trace}\left({\bm{L}}^{\dagger}_{0}+{\bm{L}}^{\dagger}_{1}-2\left({\bm{L}}_{0}^{\dagger/2}({\bm{P}}^{\top}{\bm{L}}_{1}{\bm{P}})^{\dagger}{\bm{L}}_{0}^{\dagger/2}\right)^{1/2}\right),$ and $\displaystyle d_{\text{Arithmetic}}(\eta_{{\mathcal{G}}_{0}},\eta_{{\mathcal{G}}_{1}})=\left\|{\bm{X}}_{0}-{\bm{X}}_{1}\right\|_{F}^{2}+\beta\|{\bm{L}}_{0}-{\bm{P}}^{\top}{\bm{L}}_{1}{\bm{P}}\|^{2},\forall{\bm{P}}\in\text{ potential permutation set}$ $$

## Appendix G Related Works

### G.1 Diffusion and Flow Models

Among contemporary generative models, diffusion (Ho et al., 2020) and flow models (Lipman et al., 2023) have emerged as two compelling approaches for their superior performance in generating text and images. In particular, these generative models can be unified under the framework of stochastic interpolation (Albergo & Vanden-Eijnden, 2023), which consists of four procedures (Lipman et al., 2024) as we introduced in [Section˜1](#S1). These contemporary generative models rely on constructing a probability path between data points of an easy-to-sample reference distribution and of the data distribution, and training a machine learning model to simulate the process (Lipman et al., 2024). So that one can sample from the reference (a.k.a source) distribution and iteratively transform it to approximate data samples from the target distribution. Diffusion models construct the probability path with a unilateral path conditioned on the data distribution, where one start sampling from a data point $X_{1}$ and construct the path $p({\mathcal{X}}_{t}\mid X_{1})$. While flow models can condition on either two boundary conditions, $\{X_{1},X_{0}\}$ or just one-side boundary condition $X_{1}$.

Depending on the space that the algorithm operates on, both models can be categorized into continuous or discrete models. The continuous generative models assume the data distributions are themself lying in continuous space (such as Gaussian) and build models, with examples in diffusion (Ho et al., 2020; Song et al., 2021; Wang et al., 2024) and flow (Lipman et al., 2023; Liu et al., 2023b). The discrete generateive models assume the data follows a discrete distribution, for instance categorical or Bernoulli distributions. Examples include discrete diffusion (Campbell et al., 2022; Sun et al., 2023) and discrete flow models (Campbell et al., 2024; Gat et al., 2024; Minello et al., 2025).

Under the stochastic interpolation framework, the interpolation methods are commonly selected through optimal transport (OT) displacement interpolant (Liu et al., 2023b; Albergo & Vanden-Eijnden, 2023; McCann, 1997). Optimal transport is a classical topic in mathematics that was originally used in economics and operations research (Villani & Society, 2003), and has now become a popular tool in generative models. OT aims for finding the best transport plan between two probability measures with the smallest associated transportation cost. It has been shown that generative models can be combined with technologies such as iterative matching (Tong et al., 2024) and mini batching (Pooladian et al., 2023) to approximate the OT cost, and get a significant boost in their performance in generative modeling.

### G.2 Graph Generation Models

Thanks to the capability of graphs in representing complex relationships, graph generation (Zhu et al., 2022; Liu et al., 2023a) has become an essential task in various fields such as protein design (Ingraham et al., 2019), drug discovery (Bilodeau et al., 2022), and social network analysis (Li et al., 2023). The initial attempt at graph generation is formalized through autoregression. For instance, GraphRNN (You et al., 2018) organizes the node interactions into a series of connection events and conducts autoregressive prediction for generation. Later, one shot generation methods such as Variational Graph Auto-Encoder were proposed (Kipf & Welling, 2016; Cao & Kipf, 2018).

Among various generative models, diffusion models and flow-based models have emerged as two compelling approaches for their ability to achieve state-of-the-art performance in graph generation tasks (Niu et al., 2020; Vignac et al., 2023a; Eijkelboom et al., 2024; Qin et al., 2024; Hou et al., 2024).
In the early stage, continuous diffusion models were first extended to the task of graph generation (Niu et al., 2020), where they just view the adjacency matrix as a special signal living on the $\mathbb{R}^{|{\mathcal{V}}|\times|{\mathcal{V}}|}$ domain. However, these methods fail to capture the natural discreteness of graphs, and  Vignac et al. (2023a) first brings discrete diffusion into graph generation. After that, more work (Siraudin et al., 2024; Xu et al., 2024) starts to focus on designing better discrete diffusion models for graph generation.

On the other hands, with the development of flow matching techniques, a few works have been developed to utilize flow models for graph generation and they have achieved huge success.  Eijkelboom et al. (2024) utilizes variational flow matching to process categorical data and Qin et al. (2024) developed discrete flow matching for graph generation tasks.

In parallel, there are a number of work that have managed to respect the intrinsic nature of graphs, such as global patterns. For instance, Jo et al. (2024) brings a mixture of graph technique to enhance the performance by explicitly learning final graph structures; Yu & Zhan (2025) mitigates exposure bias and reverse-start bias in graph generation; Hou et al. (2024) improves graph geneartion through optimal transport flow matching techniques but they still assume the independence between nodes and edges and use hamming distance to measure the transport cost; and  Li et al. (2023) gives the large-scale attributed graph generation framework through batching edges.

However, there remain a core challenge: constructing the probability path $p_{t}$. Existing text and image generative models, operating either in the continuous (Ho et al., 2020; Song et al., 2021; Lipman et al., 2023; Liu et al., 2023b) or discrete (Campbell et al., 2022; Sun et al., 2023; Campbell et al., 2024; Gat et al., 2024; Minello et al., 2025) space, typically rely on linear interpolation between source and target distributions to construct the path. Graph generation models, including diffusion (Niu et al., 2020; Vignac et al., 2023a; Haefeli et al., 2022; Xu et al., 2024; Siraudin et al., 2024) and flow-based models (Eijkelboom et al., 2024; Qin et al., 2024; Hou et al., 2024), inherit this design by modeling every single node and edge independently and linearly build paths in the disjoint space. However, this approach is inefficient because it neglects the strong interactions and relational structure inherent in graphs, i.e., the significance of a node heavily depends on the configuration of its neighbors. While empirical success have been achieved via fine-grained searching on the training and sampling design (Qin et al., 2024) such as target guidance and time distortion, we argue that there remains a fundamental issue of the linear probability path construction, and these strategies only mitigate the problem by manipulating the probability path.

## Appendix H Comparison with other interpolation methods

In the experimental part, we compare our methods with arithmetic (linear) interpolation, geometric interpolation and harmonic interpolation. We state the equation of them respectively as follows.

We consider the boundary graph $G_{0}$ and $G_{1}$ with $\bm{X}_{0},\bm{X}_{1}\in\mathbb{R}^{|\mathcal{V}|\times d}$ and $\bm{W}_{0},\bm{W}_{1}\in\mathbb{R}^{|\mathcal{V}|\times|\mathcal{V}|}$. Let $t\in[0,1]$, we fixed the feature interpolation as,

$$ $\bm{X}_{t}=(1-t)\,\bm{X}_{0}\;+\;t\,\bm{X}_{1},$ $$

the graph structure interpolation can be expressed as,

Linear interpolation:

$$ ${\bm{W}}_{t}=(1-t)\,{\bm{W}}_{0}\;+\;t\,{\bm{W}}_{1}.$ $$

Geometric interpolation:

$$ ${\bm{W}}_{t}={\bm{W}}_{0}^{1/2}\,\bigl{(}{\bm{W}}_{0}^{-1/2}\,{\bm{W}}_{1}\,{\bm{W}}_{0}^{-1/2}\bigr{)}^{t}\,{\bm{W}}_{0}^{1/2},$ $$

Harmonic interpolation:

$$ ${\bm{W}}_{t}=\Bigl{(}(1-t)\,{\bm{W}}_{0}^{-1}\;+\;t\,{\bm{W}}_{1}^{-1}\Bigr{)}^{-1}.$ $$

Each interpolation methods actually handle each special manifold assumption, which should be designed under a comprehensive understanding of the task. In our experimental part, we conduct intensive analysis on the impact of interpolation methods to the graph generation quality.

## Appendix I Additional Experiment Results

### I.1 Experiment setups and computational cost

The training and sampling computation time are provided in [Table˜6](#A9.T6). The experiments were run on a single NVIDIA A100-SXM4-80GB GPU. The hyperparameter configuration in producing [Tables˜1](#S4.T1), [2](#S4.T2) and [3](#S4.T3) is reported in [Table˜7](#A9.T7).

**Table 6: Training and sampling time on each dataset. TG means using target-guided velocity; BW means using BW velocity.**
| Dataset | Min Nodes | Max Nodes | Training Time (h) | Graphs Sampled | Sampling Time (h) |
| --- | --- | --- | --- | --- | --- |
| Planar | 64 | 64 | 45 (1.55x) | 40 | 0.07(TG); 0.13 (BW) |
| Tree | 64 | 64 | 10 (1.25x) | 40 | 0.07(TG);0.14(BW) |
| SBM | 44 | 187 | 74 (0.98x) | 40 | 0.07(TG) 0.14(BW) |
| Moses | 8 | 27 | 35 (0.76x) | 25000 | 5(TG); 6(BW) |
| Guacamol | 2 | 88 | 251 (1.8x ) | 10000 | 7(TG); 21(BW) |
| QM9 | 3 | 29 | 15 | 25000 | 5(TG) 6(BW) |
| GEOM | / | 181 | 141 | 10000 | 7(TG) 14(BW) |

**Table 7: Best Configuration for Training and Sampling when producing [Tables˜1](#S4.T1), [2](#S4.T2) and [3](#S4.T3).**
|  | Train | Sampling |  |  |
| --- | --- | --- | --- | --- |
| Dataset | Initial Distribution | Train Distortion | Sample Distortion | Sampling steps |
| Planar | Marginal | Identity | Identity | 1000 |
| Tree | Marginal | Polydec | Polydec | 1000 |
| SBM | Absorbing | Identity | Identity | 1000 |
| MOSES | Marginal | Identity | Identity | 500 |
| GUACAMOL | Marginal | Identity | Identity | 500 |
| QM9 | Marginal | Identity | Identity | 500 |
| GEOM | Marginal | Identity | Identity | 500 |

### I.2 Additional results for the training paths

Figure: (a) Probability Path for Planar Graphs
Refer to caption: /html/2506.14020/assets/figs/png_output/prob_path_planar.png

[Figure˜6](#A9.F6) gives the training probability path construction for planar graphs and tree graphs. While planar graphs have a similar pattern as the SBM datasets as in [Figure˜3(a)](#S4.F3.sf1), the probability path constructed for tree graphs does not follow a similar pattern. We attribute this to the different geometry of tree graphs that reside in hyperbolic space (Yang et al., 2022).

### I.3 More experiments on plain graph generations

#### Additional results for sampling paths.

We then give the sampling path construction in  [Figure˜7](#A9.F7). To better illustrate the advantage of BWFlow, we fix the sampling steps to be as small as 50. It is clear that in planar and SBM dataset, the BW velocity can still provide a smooth probability and stable convergence towards the data distribution. While the linear velocity does not give a good probability path and fails to converge to the optimal value, especially when the sampling size is small.

The maximum mean discrepancy (MMD) of four graph statistics between the set of generated graphs
and the test set is measured, including degree (Deg.), clustering coefficient (Clus.),
count of orbits with 4 nodes (Orbit), the eigenvalues of
the graph Laplacian (Spec.), wavelet ratio (Wavelet.). To verify that the model learns to generate graphs with valid topology, we gives the percentage of valid,
unique, and novel (V.U.N.) graphs for where a valid graph satisfies the corresponding property of each dataset (Planar, Tree, SBM, etc.).

Figure: (a) Sampling path comparision on Planar dataset
Refer to caption: /html/2506.14020/assets/x7.png

Figure: (a) QM9 with explicit hydrogen
Refer to caption: /html/2506.14020/assets/figs/png_output/training_curve_qm9.png

**Table 8: Graph Generation Performance on Synthetic Graphs. Results are obtained through tuning the probability path manipulation techniques. The remaining values are obtained from Qin et al. (2024).**
|  | Planar Dataset |  |  |  |  |
| --- | --- | --- | --- | --- | --- |
| Model | Ratio $\downarrow$ | Valid $\uparrow$ | Unique  $\uparrow$ | Novel $\uparrow$ | V.U.N.  $\uparrow$ |
| Train set | 1.0 | 100 | 100 | 0.0 | 0.0 |
| GRAN (Liao et al., 2019) | 2.0 | 97.5 | 85.0 | 2.5 | 0.0 |
| SPECTRE (Martinkus et al., 2022b) | 3.0 | 25.0 | 100 | 100 | 25.0 |
| DiGress (Vignac et al., 2023a) | 5.1 | 77.5 | 100 | 100 | 77.5 |
| EDGE (Chen et al., 2023) | 431.4 | 0.0 | 100 | 100 | 0.0 |
| BwR (Diamant et al., 2023) | 251.9 | 0.0 | 100 | 100 | 0.0 |
| BiGG (Dai et al., 2020) | 16.0 | 62.5 | 85.0 | 42.5 | 5.0 |
| GraphGen (Goyal et al., 2020) | 210.3 | 7.5 | 100 | 100 | 7.5 |
| HSpectre (one-shot) (Bergmeister et al., 2024) | 1.7 | 67.5 | 100 | 100 | 67.5 |
| HSpectre (Bergmeister et al., 2024) | 2.1 | 95.0 | 100 | 100 | 95.0 |
| GruM (Jo et al., 2024) 1.8 | — | — | — | 90.0 |  |
| CatFlow (Eijkelboom et al., 2024) | — | — | — | — | 80.0 |
| DisCo (Xu et al., 2024) | — | 83.6 ±2.1 | 100.0 ±0.0 | 100.0 ±0.0 | 83.6 ±2.1 |
| Cometh - PC (Siraudin et al., 2024) | — | 99.5 ±0.9 | 100.0 ±0.0 | 100.0 ±0.0 | 99.5 ±0.9 |
| DeFoG | 1.6 ±0.4 | 99.5 ±1.0 | 100.0 ±0.0 | 100.0 ±0.0 | 99.5 ±1.0 |
| BWFlow | 1.3 ±0.4 | 97.5 ±2.5 | 100.0 ±0.0 | 100.0 ±0.0 | 97.5±2.5 |
|  | Tree Dataset |  |  |  |  |
| Train set | 1.0 | 100 | 100 | 0.0 | 0.0 |
| GRAN (Liao et al., 2019) | 607.0 | 0.0 | 100 | 100 | 0.0 |
| DiGress (Vignac et al., 2023a) | 1.6 | 90.0 | 100 | 100 | 90.0 |
| EDGE (Chen et al., 2023) | 850.7 | 0.0 | 7.5 | 100 | 0.0 |
| BwR (Diamant et al., 2023) | 11.4 | 0.0 | 100 | 100 | 0.0 |
| BiGG (Dai et al., 2020) | 5.2 | 100 | 87.5 | 50.0 | 75.0 |
| GraphGen (Goyal et al., 2020) | 33.2 | 95.0 | 100 | 100 | 95.0 |
| HSpectre (one-shot) (Bergmeister et al., 2024) | 2.1 | 82.5 | 100 | 100 | 82.5 |
| HSpectre (Bergmeister et al., 2024) | 4.0 | 100 | 100 | 100 | 100 |
| DeFoG | 1.6 ±0.4 | 96.5 ±2.6 | 100.0 ±0.0 | 100.0 ±0.0 | 96.5 ±2.6 |
| BWFlow | 1.4±0.3 | 75.5 ±2.4 | 100.0±0 | 100.0±0 | 75.5±2.4 |
|  | Stochastic Block Model ($n_{\text{max}}=187$, $n_{\text{avg}}=104$) |  |  |  |  |
| Model | Ratio $\downarrow$ | Valid $\uparrow$ | Unique  $\uparrow$ | Novel $\uparrow$ | V.U.N.  $\uparrow$ |
| Training set | 1.0 | 85.9 | 100 | 0.0 | 0.0 |
| GraphRNN (You et al., 2018) | 14.7 | 5.0 | 100 | 100 | 5.0 |
| GRAN (Liao et al., 2019) | 9.7 | 25.0 | 100 | 100 | 25.0 |
| SPECTRE (Martinkus et al., 2022b) | 2.2 | 52.5 | 100 | 100 | 52.5 |
| DiGress (Vignac et al., 2023a) | 1.7 | 60.0 | 100 | 100 | 60.0 |
| EDGE (Chen et al., 2023) | 51.4 | 0.0 | 100 | 100 | 0.0 |
| BwR (Diamant et al., 2023) | 38.6 | 7.5 | 100 | 100 | 7.5 |
| BiGG (Dai et al., 2020) | 11.9 | 10.0 | 100 | 100 | 10.0 |
| GraphGen (Goyal et al., 2020) | 48.8 | 5.0 | 100 | 100 | 5.0 |
| HSpectre (one-shot) (Bergmeister et al., 2024) | 10.5 | 75.0 | 100 | 100 | 75.0 |
| HSpectre (Bergmeister et al., 2024) | 10.2 | 45.0 | 100 | 100 | 45.0 |
| GruM (Jo et al., 2024) | 1.1 | — | — | — | 85.0 |
| CatFlow (Eijkelboom et al., 2024) | — | — | — | — | 85.0 |
| DisCo (Xu et al., 2024) | — | 66.2 ±1.4 | 100.0 ±0.0 | 100.0 ±0.0 | 66.2 ±1.4 |
| Cometh (Siraudin et al., 2024) | — | 75.0 ±3.7 | 100.0 ±0.0 | 100.0 ±0.0 | 75.0 ±3.7 |
| DeFoG | 4.9 ±1.3 | 90.0 ±5.1 | 100.0 ±0.0 | 100.0 ±0.0 | 90.0 ±5.1 |
| BWFlow | 3.8±0.9 | 90.5 ±4.0 | 100.0 ±0.0 | 100.0 ±0.0 | 90.5 ±4.0 |

#### Full results for plain graph generation.

[Table˜8](#A9.T8) gives the full results with other generative models aside from the diffusion and flow models. [Table˜9](#A9.T9) gives the results on smaller datasets, i.e., comm20

**Table 9: Quantitative experimental results on COMM20 (smaller dataset).**
| FM type | Interpolation | comm-20 |  |  |  |
| --- | --- | --- | --- | --- | --- |
| Deg. | Clus. | Orbit. |  |  |  |
| Discrete | Arithmetic | 0.071 | 0.115 | 0.037 |  |
| Harmonic | 0.011 | 0.036 | 0.027 |  |  |
| Geometric | 0.047 | 0.083 | 0.02 |  |  |
| BW | 0.009 | 0.013 | 0.017 |  |  |

### I.4 Full table for the synthetic graph generation

In [Table˜10](#A9.T10), we illustrate the numerical results for comparing the interpolation methods in plain graph generation without node features. It is clear that BWFlow outperforms other methods in planar and SBM graphs. But the performance was not good in tree graph generations.

**Table 10: Ablation study on interpolation methods when probability path manipulation techniques are all disabled. The clustering and orbit ratios in tree graphs are omitted, given that in the training set, the corresponding statistics are 0. The results go over Exponential Moving Average (decay 0.999) for the last 5 checkpoints. The table is produced with Marginal boundary distributions, without time distortion.**
| Dataset | FM type | Interpolation | V.U.N metrics | Spectral Metrics |  |  |  |  |  |  |  |
| --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- |
| Novelty | Uniqueness | Validity | Orbit | Spec | Clustering | Degree | Wavelet | Avg. Ratio |  |  |  |
| Planar | Discrete | DeFog | 100 | 100 | 78.25 | 8.98 | 1.45 | 2.09 | 2.65 | 2.38 | 3.51 |
| Arithmetic | 100 | 100 | 73.25 | 10.83 | 1.33 | 1.74 | 2.24 | 2.39 | 3.70 |  |  |
| harmonic | 100 | 100 | 0.00 | 4519.63 | 2.57 | 17.01 | 25.10 | 42.41 | 921.35 |  |  |
| Geometric | 100 | 100 | 23.25 | 655.66 | 1.61 | 10.17 | 13.25 | 6.68 | 137.47 |  |  |
|  |  | BW | 100 | 100 | 84.75 | 5.14 | 1.27 | 1.69 | 1.78 | 2.02 | 2.38 |
| Tree | Discrete | Defog | 100 | 100 | 61.53 | / | 1.17 | / | 1.27 | 1.51 | 1.32 |
| Arithmetic | 100 | 100 | 56.91 | / | 1.16 | / | 1.04 | 1.45 | 1.22 |  |  |
| harmonic | 100 | 100 | 0.53 | / | 2.32 | / | 1.93 | 3.31 | 2.52 |  |  |
| Geometric | 100 | 100 | 48.38 | / | 1.62 | / | 2.13 | 2.10 | 1.94 |  |  |
|  |  | BW | 100 | 100 | 51.45 | / | 1.58 | / | 2.56 | 2.13 | 2.09 |
| SBM | Discrete | Arithmetic | 100 | 100 | 44.63 | 2.57 | 1.40 | 1.46 | 15.55 | 7.88 | 5.77 |
| harmonic | 100 | 100 | 9.73 | 3.10 | 10.23 | 1.59 | 172.10 | 103.04 | 58.10 |  |  |
| Geometric | 100 | 100 | 0.0 | 3.11 | 4.45 | 1.80 | 150.41 | 60.60 | 54.65 |  |  |
| BW | 100 | 100 | 58.70 | 2.03 | 1.50 | 1.50 | 9.04 | 8.41 | 4.51 |  |  |

### I.5 3D Molecule Generation: QM9 without explicit Hydrogen

In [Table˜11](#A9.T11) we report the results of QM9 without explicit hydrogen. This task is relatively easy compared to the generation task with explicit hydrogen, and both Midi and our BWFlow have achieved near-saturated performance with validity near to 100%.

**Table 11: Quantitative experimental results on QM9 datasets without explicit hydrogen in 3D molecule generation.**
| Dataset | Interpolation | Metrics |  |  |  |  |  |
| --- | --- | --- | --- | --- | --- | --- | --- |
| $\mu$ | V.U.N(%) | Connected(%) | Charges($10^{-2}$) | Atom($10^{-2}$) | Angles(^∘) |  |  |
| QM9<br>(w/o h) | MiDi | 1.00 | 98.0 | 100.0 | 0.4 | 5.1 | 1.49 |
| Linear Flow | 1.60 | 79.33 | 52.3 | 0.7 | 14.0 | 8.77 |  |
| BWFlow | 1.02 | 99.8 | 100.0 | 0.4 | 4.8 | 1.53 |  |

### I.6 Convergence Analysis

[Figure˜8](#A9.F8) are the training convergence analysis on Planar and QM9 dataset, showing that BWFlow provides a fast convergence speed than others.

**Table 12: Graph generation performance on the synthetic datasets: Planar, Tree, and SBM. Given that the synthetic datasets are usually unstable in evaluation, we applied an exponential moving average to stabilize the results and sample 5 times (each run generates 40 graphs) to calculate the mean and standard deviation. The experiment settings are in [Table˜7](#A9.T7)**
|  |  | Planar | Tree | SBM |  |  |  |
| --- | --- | --- | --- | --- | --- | --- | --- |
| Model | Class | V.U.N. $\uparrow$ | A.Ratio $\downarrow$ | V.U.N. $\uparrow$ | A.Ratio $\downarrow$ | V.U.N. $\uparrow$ | A.Ratio $\downarrow$ |
| Train set | — | 100 | 1.0 | 100 | 1.0 | 85.9 | 1.0 |
| DiGress (Vignac et al., 2023a) | Diffusion | 77.5 | 5.1 | 90.0 | 1.6 | 60.0 | 1.7 |
| EDGE (Chen et al., 2023) | Diffusion | 0.0 | 431.4 | 0.0 | 850.7 | 0.0 | 51.4 |
| HSpectre (Bergmeister et al., 2024) | Diffusion | 95.0 | 2.1 | 100.0 | 4.0 | 75.0 | 10.5 |
| GruM (Jo et al., 2024) | Diffusion | 90.0 | 1.8 | — | — | 85.0 | 1.1 |
| CatFlow (Eijkelboom et al., 2024) | Flow | 80.0 | — | — | — | 85.0 | — |
| DisCo (Xu et al., 2024) | Diffusion | 83.6±2.1 | — | — | — | 66.2±1.4 | — |
| Cometh (Siraudin et al., 2024) | Diffusion | 99.5±0.9 | — | — | — | 75.0±3.7 | — |
| DeFoG (Qin et al., 2024) | Flow | 99.5±1.0 | 1.6±0.4 | 96.5±2.6 | 1.6±0.4 | 90.0±5.1 | 4.9±1.3 |
| BWFlow | Flow | 97.5±2.5 | 1.3±0.4 | 75.5±2.4 | 1.4±0.3 | 90.5±4.0 | 3.8±0.9 |
| DeFoG (EMA) | Flow | 77.5±8.37 | 3.5±1.7 | 73.1±11.4 | 1.50±0.3 | 85.0±7.1 | 3.7±0.9 |
| BWFlow (EMA) | Flow | 84.8±6.44 | 2.4±0.9 | 68.1±12.2 | 1.24±0.2 | 83.5±6.0 | 2.4±0.6 |

**Table 13: Comparison of interpolation methods on 3D Molecule Generation with explicit hydrogen in QM9 dataset.**
| Flow Type | Interpolation | Metrics |  |  |  |  |  |  |  |
| --- | --- | --- | --- | --- | --- | --- | --- | --- | --- |
| $\mu$ | V.U.N(%) | Mol Stable | Atom Stable | Connected(%) | Charge($10^{-2}$) | Atom($10^{-2}$) | Angles(^∘) |  |  |
| Discrete | MiDi | 1.01 | 93.13 | 93.98 | 99.60 | 99.21 | 0.2 | 3.7 | 2.21 |
| Arithmetic | 1.01 | 87.53 | 88.45 | 99.13 | 99.09 | 0.4 | 4.2 | 2.72 |  |
| harmonic | 1.01 | 94.91 | 94.54 | 99.65 | 99.03 | 0.6 | 6.4 | 2.21 |  |
| Geometric | 1.01 | 91.26 | 91.29 | 99.42 | 98.42 | 0.1 | 4.4 | 3.63 |  |
| BW | 1.01 | 96.45 | 97.84 | 99.84 | 99.24 | 0.1 | 2.3 | 1.96 |  |
| Continuous | Arith | 2.15 | 25.45 | 10.23 | 76.85 | 28.82 | 0.7 | 5.6 | 14.47 |
| harmonic | 1.01 | 11.38 | 11.64 | 73.48 | 99.65 | 1.2 | 17.2 | 15.04 |  |
| Geometric | 1.00 | 42.07 | 46.08 | 91.13 | 99.87 | 1.0 | 12.7 | 8.03 |  |
| BW | 1.02 | 62.02 | 61.76 | 95.99 | 97.72 | 0.6 | 8.7 | 7.80 |  |

**Table 14: Large molecule generation results. Only iterative denoising-based methods are reported.**
|  | Guacamol | MOSES |  |  |  |  |  |  |  |  |  |  |
| --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- |
| Model | Val. $\uparrow$ | V.U. $\uparrow$ | V.U.N.$\uparrow$ | KL div $\uparrow$ | FCD $\uparrow$ | Val. $\uparrow$ | Unique. $\uparrow$ | Novelty $\uparrow$ | Filters $\uparrow$ | FCD $\downarrow$ | SNN $\uparrow$ | Scaf $\uparrow$ |
| Training set | 100.0 | 100.0 | 0.0 | 99.9 | 92.8 | 100.0 | 100.0 | 0.0 | 100.0 | 0.01 | 0.64 | 99.1 |
| DiGress (Vignac et al., 2023a) | 85.2 | 85.2 | 85.1 | 92.9 | 68.0 | 85.7 | 100.0 | 95.0 | 97.1 | 1.19 | 0.52 | 14.8 |
| DisCo (Xu et al., 2024) | 86.6 | 86.6 | 86.5 | 92.6 | 59.7 | 88.3 | 100.0 | 97.7 | 95.6 | 1.44 | 0.50 | 15.1 |
| Cometh (Siraudin et al., 2024) | 98.9 | 98.9 | 97.6 | 96.7 | 72.7 | 90.5 | 99.9 | 92.6 | 99.1 | 1.27 | 0.54 | 16.0 |
| DeFoG (Qin et al., 2024) | 99.0 | 99.0 | 97.9 | 97.7 | 73.8 | 92.8 | 99.9 | 92.1 | 98.9 | 1.95 | 0.55 | 14.4 |
| BWFlow (Ours) | 98.8 | 98.9 | 97.4 | / | 69.2 | 92.0 | 100.0 | 94.5 | 98.4 | 1.32 | 0.56 | 15.3 |

**Table 15: Large molecule generation results. Only comparing the representative diffusion and flow models. B.E. is the scenario that only considers binary edge types. The results are almost saturated, thus not very informative.**
|  | Guacamol | MOSES |  |  |  |  |
| --- | --- | --- | --- | --- | --- | --- |
| Model | Val. $\uparrow$ | V.U. $\uparrow$ | V.U.N.$\uparrow$ | Val. $\uparrow$ | Unique. $\uparrow$ | Novelty $\uparrow$ |
| Digress (B.E.) | 96.0 | 98.9 | 97.4 | 96.1 | 100 | 100 |
| Defog (B.E.) | 98.4 | 98.4 | 97.9 | 99.3 | 100 | 100 |
| BWFlow (B.E.) | 98.0 | 98.0 | 97.7 | 99.6 | 100 | 100 |

## Peer Reviews

confidence: 2
contribution: 3
initial_rating: 8
presentation: 3
questions: Please check the above questions.
soundness: 3
strengths: S1. Good motivation. "MRFs organize the nodes/edges as an interconnected system and interpolating between two MRFs captures the joint evolution of the graph system" sounds reasonable to me.

S2. This paper's presentation is nice. Even without checking the detailed mathematical derivations, I can still follow the story in the main content.

S3. The experimental results seem impressive.
summary: This paper studies the graph generation problem, via flow (matching) model. The notable point is that this paper formulates the graph with Markov Random Fields (MRFs) so that the interpolation on GraphMRFs can include the interaction between nodes and edges, greatly different from the previous work, whose node and edge interpolations are independent (Eq. 5). Some followup results, e.g., what is the closed-form interpolation, and the velocity at the interpolation point, are presented.
weaknesses: Generally, i am satisfied about the quality of this paper. Here are some questions/suggestion which might be able to improve this paper.

Q1. From lines 66 to 82, it introduces the motivating example of this paper. It is a bit confusing why t=0.8 is such an important point. If this is purely empirical observation, or there are some mathematical reasons?

Q2. Can we understand the proposed method as a kind of latent diffusion/flow method, just like the stable diffusion did? I understand that the proposed method does not include an explicit encoder/decoder, and the I like the interpolations/velocity derived based on the colored
Gaussian distribution. Asking because the Figure 2 somewhat shows such a connection between flows in raw graph space and latent space.

Q3. About the OT paths. To be specific,

Q3.1. According to the Algorithm 1, the interpolation is on the OT path between a pair of G_0 and G_1, but not the in-batch OT coupling between a batch of {G_0} and {G_1}, is that right? But I think it should be easy to generalize to the mini-batch OT coupling setting as [1] did, is that right?

Q3.2 I can buy the story that "capture global co-evolution (Haasler & Frossard, 2024) of the graph components" for the interpolation is important, but whether the OT interpolation/path really matters? Asking because in many recent flow model studies, the straight line interpolation's best advantage is more on the efficiency (so fewer sampling steps), but not on the effectiveness.

[1] Improving and generalizing flow-based generative models with minibatch optimal transport

---

confidence: 4
contribution: 3
initial_rating: 6
presentation: 3
questions: Given the computational overhead, would a hybrid schedule be feasible—one that leverages linear interpolation for early stages and switches to the proposed BWFlow when most beneficial?


In Section A.3 there are cases where GraphMRF may be suboptimal for certain graph structures. Can structural information be injected to improve performance—for example, via higher-order graph features, subgraph/ motif features, or structural role embeddings?

line 479 "When scaled up to large
but sparse graphs, the complexity can be reduced to O(TN2)", what is T
soundness: 3
strengths: The paper is well-written, clear, and methodologically sound, presenting complex ideas with precision and rigor.

It establishes a theoretical framework for probability-path construction directly within Markov Random Field (MRF) space, enabling coherent joint node–edge evolution that respects graph-structured dependencies.

It introduces a novel integration of MRFs with optimal transport, using MRFs as the ambient representation space for probability paths in a way that naturally aligns with graph topology and conditional dependencies.
summary: The paper proposes BWFlow, a flow-matching framework for graph generation that constructs a smooth, joint probability path by modeling graphs as Markov random fields (MRFs) and leveraging optimal transport between MRFs. It argues that existing diffusion/flow methods use independent, linear node/edge interpolations that disrupt graph dependencies, leading to irregular paths and unstable training. By enforcing co-evolution of nodes and edges along the derived optimal path, BWFlow aims to improve training dynamics and sampling convergence. Experiments on plain graph and molecular generation report competitive or better performance, faster convergence, and more efficient sampling.
weaknesses: While acknowledged in the limitations, BWFlow incurs substantial computational overhead, and reducing this cost appears non-trivial.

The experimental evaluation omits several GNN-based baselines, including DisCo (Xu et al., 2024), limiting the completeness of the comparison.

---

confidence: 2
contribution: 3
initial_rating: 6
presentation: 2
questions: 1. In Fig. 1(c), when flow steps = 1.0, the Time Distortion method yields a slightly smaller value than the BW interpolation. Does this imply that Time Distortion is comparably efficient in this setting?

2. Could the authors elaborate on the modeling assumption $p(\mathcal{G};\mathbf{G})=p(\mathcal{X};\mathbf{X},\mathbf{W})\cdot p(\mathcal{E};\mathbf{W})$? Why are the node features and edge structures conditionally independent, except for the coupling through 𝑊? How critical is this assumption for the GraphMRF formulation?

3. Does the model scale to large graphs?
soundness: 3
strengths: 1. The paper presents a mathematically coherent framework connecting graph optimal transport and flow matching. The closed-form derivations of the BW interpolation and velocity fields are valuable theoretical contributions.

2. Using the BW metric on graph Laplacians provides a geometrically meaningful way to measure and interpolate between graph distributions, which addresses a long-standing issue of linear interpolation violating graph manifold constraints.

3. The experimental section is relatively comprehensive, covering various types of graph data, including plain graphs, molecular graphs, and 3D molecular structures.
summary: This paper proposes Bures–Wasserstein Flow Matching (BWFlow), a generative modeling framework for graphs that integrates flow matching with the Bures–Wasserstein (BW) geometry on a Graph Markov Random Field (GraphMRF) representation. The authors first formalize graphs as random fields on nodes and edges (GraphMRF), establishing a Gaussian joint distribution whose covariance is tied to the graph Laplacian. They then derive a closed-form expression of the Bures–Wasserstein distance between two graphs, together with an analytical interpolation path and conditional velocity field, which enable flow matching without simulation. The method is further extended to a discrete version (with categorical node features and Bernoulli edges), and training/inference procedures are given. Comprehensive experiments on graph generation benchmarks are reported.
weaknesses: 1. The methodological novelty is mainly compositional: the work combines existing Bures–Wasserstein distance (already applied to graph covariance learning) with the established flow-matching framework under a Graph MRF formulation, resulting in limited originality.

2. The empirical evaluation lacks consistency with standard benchmarks (e.g., DiGress, GruM, DeFoG) and omits common molecular-graph metrics such as Validity, Uniqueness, and Novelty. The results would be more convincing if experiments used the same datasets and evaluation metrics as prior baselines.

3. The paper does not isolate the effect of replacing Euclidean interpolation with the BW distance, nor does it provide a systematic comparison between the continuous and discrete variants. Without such analysis, it is unclear whether BW geometry truly contributes to the observed performance.

4. The notation system is heavy, and some equations are symbolically dense without intuitive explanation or visualization.

5. The paper does not mention a code release to ensure reproducibility.

---

confidence: 3
contribution: 3
initial_rating: 8
presentation: 3
questions: 1. I would like to clarify the experimental setting details for Figure 1a in the Introduction.
2. I would like the authors to address the question regarding the discussion of the results in Figure 3a, which is the second point under Weaknesses.
soundness: 3
strengths: 1. (Originality) In graph generation tasks using flow matching, incorporating graph structure is non-trivial. This paper solves this issue by modeling a graph using the Markov Random Field. To my knowledge, this is the first paper to propose such a solution.
2. (Quality) There is design flexibility in how to set the velocity field in flow matching. This paper provides a theoretical background by employing the Optimal Transport. Furthermore, while optimal transport computations are typically high, the proposed method overcomes this problem using the Markov Random Field. This makes the approach practically implementable.
3. (Clarity) The writing is clear. It explains the necessary background knowledge (Markov Random Fields, Optimal Transport) to understand the proposed method. The mathematical descriptions are also appropriate. Therefore, I had no significant difficulties in understanding the paper's main points.
4. (Significance) Numerical experiments demonstrate generation performance of the proposed method surpasses existing methods, which strengthens the significance of the work in the area of graph generation research.
summary: This paper proposes a novel graph generation method using flow matching. First, it is pointed out that methods generating graphs via linear interpolation fail to smoothly change graph statistics, leading to poor flow velocity estimation accuracy and convergence gaps. Building upon this, the proposed method models a graph using Markov Random Field, estimates a vector field derived by the optimal transport on this model using the flow matching. Numerical experiments applied the proposed method to three graph generation tasks, 2D and 3D molecular graph generation tasks, verifying its generation performance.
weaknesses: 1. The paper points out that the problem with the existing methods is that they ignore graph structure and generate nodes and edges independently. However, the experiments verify the limitations of linear interpolations. Therefore, it seems unclear whether the problem lies in ignoring graph structure or in using linear interpolations.
2. I have a question about the interpretation of the results in Figure 3a. Looking at Figure 1a of the introduction, the authors assume that ideal interpolation changes graph statistics linearly. However, the discussion on Figure 3a suggests that it is desirable that statistics first increase and then decay. It seemingly contradicts the claim in Figure 1a and it would be more natural to claim that linear interpolation is preferable.


Decision: Accept (Poster)