Paper deep dive
Stability of Flow Models for Graph Signals
Martin Schmidt, Gonzalo Mateos
Intelligence
Status: succeeded | Model: Gemma-4-26B-A4B | Prompt: intel-v1 | Confidence: 93%
Last extracted: 7/18/2026, 8:38:06 AM
Summary
This paper analyzes the stability of continuous normalized flow (CNF) models parameterized by Graph Neural Networks (GNNs) for generating graph signals. It proves that permutation equivariance is preserved in the continuous-time dynamics and discrete approximations. The authors derive explicit stability bounds quantifying how structural perturbations in the graph topology affect the generated probability distributions. Based on these bounds, they propose a stability-promoting regularized flow matching strategy that penalizes the spatial Lipschitz constant of the vector field, demonstrating improved robustness to structural noise in synthetic and real-world fMRI experiments.
Entities (10)
Relation Signals (6)
Continuous Normalized Flow (CNF) → parameterizedby → Graph Neural Network (GNN)
confidence 98% · we analyze continuous normalized flow models parameterized by GNNs
Continuous Normalized Flow (CNF) → exhibitsproperty → Permutation Equivariance
confidence 97% · show that permutation equivariance is preserved for both the resulting continuous-time ordinary differential equations
Graph Neural Network (GNN) → providesinductivebias → Permutation Equivariance
confidence 96% · Permutation equivariance is a non-negotiable desideratum for graph-structured data... infuse inductive biases directly into the network architecture
Stability Bounds → quantifieseffectof → Structural Perturbations
confidence 95% · derive explicit stability bounds on the generated probability distributions, which quantify how relative graph perturbations affect the final sampled signals
Flow Matching → incorporatesregularization → Spatial Lipschitz Constant
confidence 94% · introduce a stability-promoting regularized flow matching strategy that actively penalizes the spatial Lipschitz constant
Wasserstein Distance → usedtomeasure → Stability
confidence 90% · quantified in terms of Wasserstein distance under perturbations of the graph topology
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:Generating signals on graphs requires permutation-equivariant models that exhibit stability with respect to relative structural perturbations. While favorable stability properties of Graph Neural Networks (GNNs) have been well documented, it is unclear how structural errors propagate through the dynamics of continuous generative flow models that are gaining traction for graph signal generation. In this paper, we analyze continuous normalized flow models parameterized by GNNs and show that permutation equivariance is preserved for both the resulting continuous-time ordinary differential equations and their discrete numerical approximations used as graph signal samplers. Our primary contribution is to derive explicit stability bounds on the generated probability distributions, which quantify how relative graph perturbations affect the final sampled signals. Motivated by these theoretical bounds, we introduce a stability-promoting regularized flow matching strategy that actively penalizes the spatial Lipschitz constant of the vector field during model training. Experiments using synthetic smooth signals on stochastic block model graphs and real-world fMRI signals on brain connectomes demonstrate that this bound-oriented approach yields generative models that are more robust to structural noise, without sacrificing output quality.
Tags
Links
- Source: https://arxiv.org/abs/2607.07510v1
- Canonical: https://arxiv.org/abs/2607.07510v1
PDF not stored locally. Use the link above to view on the source site.
Full Text
108,076 characters extracted from source content.
Expand or collapse full text
Stability of Flow Models for Graph Signals Martin Schmidt and Gonzalo Mateos This work was supported in part by NSF under Grant ECCS 2231036. (Corresponding author: Gonzalo Mateos.) Martin Schmidt and Gonzalo Mateos are with the Dept. of Electrical and Computer Engineering, University of Rochester, Rochester, NY 14627, USA (e-mails: mschmi21@ur.rochester.edu; gmateosb@ece.rochester.edu). Abstract Generating signals on graphs requires permutation-equivariant models that exhibit stability with respect to relative structural perturbations. While favorable stability properties of Graph Neural Networks (GNNs) have been well documented, it is unclear how structural errors propagate through the dynamics of continuous generative flow models that are gaining traction for graph signal generation. In this paper, we analyze continuous normalized flow models parameterized by GNNs and show that permutation equivariance is preserved for both the resulting continuous-time ordinary differential equations and their discrete numerical approximations used as graph signal samplers. Our primary contribution is to derive explicit stability bounds on the generated probability distributions, which quantify how relative graph perturbations affect the final sampled signals. Motivated by these theoretical bounds, we introduce a stability-promoting regularized flow matching strategy that actively penalizes the spatial Lipschitz constant of the vector field during model training. Experiments using synthetic smooth signals on stochastic block model graphs and real-world fMRI signals on brain connectomes demonstrate that this bound-oriented approach yields generative models that are more robust to structural noise, without sacrificing output quality. I Introduction Generative models are powerful tools for learning underlying data distributions from finite samples. Recent advances in continuous-time generative modeling, including flow matching [29], diffusion models [16], and neural ordinary differential equations (ODEs) [9], have achieved remarkable success generating realistic images, audio and video. However, their application to non-Euclidean domains is still developing. Specifically, while considerable effort has gone into generating graph topologies [31, 34, 23, 30, 47, 8], significantly less attention has been paid to generating signals supported on the nodes of a given graph structure [39, 46, 52, 49], despite their ubiquity in real-world applications [35]. Such signals arise in diverse contexts, including neural activity over structural brain connectomes [18], traffic flows on transportation networks [50], sensor measurements in distributed systems [22], and power injections in electrical grids [38]. Learning to faithfully generate such signals can then be useful for various reasons, e.g., to produce synthetic samples in low-data regimes, to serve as implicit priors in inverse or reconstruction problems, and to enable simulation of realistic system behavior. For instance, recent works have leveraged graph signal diffusion models for probabilistic forecasting of stock prices to capture uncertainties and tail events [46], modeling correlated user-item interactions for collaborative filtering in recommender systems [52], and optimizing wireless resource allocation [45]. While said application-driven impetus has fueled exciting architectural advances [39, 46, 52, 49], the present paper studies fundamental equivariance and stability properties of flow models for graph signal generation. In addition to addressing theoretical questions left unanswered by prior work, the upshot of our analysis has practical implications to model training. Permutation equivariance and stability. Because high-quality data can be both limited and expensive, neural networks must leverage structural priors to learn effectively. To this end, one can e.g., incorporate prior knowledge through data augmentation [1] or penalties in the loss function [37]. When possible, arguably the preferred approach is to infuse inductive biases directly into the network architecture [4]. Permutation equivariance is a non-negotiable desideratum for graph-structured data, which guarantees that the generative dynamics natively respect the graph’s structural symmetries. Beyond achieving good performance on training data, machine learning models are expected to generalize well to unseen samples. A key principle underlying such behavior is stability: models that respond smoothly to small perturbations in their input tend to generalize better and exhibit improved robustness and transferability [44, 3]; learning with graphs is no exception [53]. Graph-structured data are inherently irregular and often noisy, and in many applications the observed topology is only an imperfect reflection of the underlying system [27, Ch. 7]. For instance, graphs constructed from empirical measurements—such as correlation networks [51] derived from finite-length time series—can vary significantly with the sampling window, leading to noisy graph representations [32, 11]. Crucially, measuring a model’s robustness to such topological noise requires comparing predictive outputs across different graph structures. Because node indexing is often arbitrary, any formal measure of structural stability must be defined modulo permutation [13]. Consequently, permutation equivariance is not merely a useful inductive bias, but a strict mathematical prerequisite for analyzing the structural stability of these graph-parametrized learning models. Innovations in context. While permutation equivariance, structural stability, and transferability have been rigorously established for Graph Neural Networks (GNNs) [13, 40, 24, 28], these properties have primarily been studied in the context of static, discriminative tasks. In contrast, continuous-time generative models evolve a state trajectory that maps a simple prior distribution to the data distribution. Although related equivariance results exist for flow models [26], in our graph signal processing (GSP) setting equivariance must account for coupled permutations of both the graph signal and the underlying topology. Recent works in graph signal generation have primarily focused on model design rather than theoretical guarantees. In particular, [39] incorporates the graph structure into a diffusion-style forward process via a heat equation, while [46] proposes a U-Net–inspired GNN architecture. However, these approaches do not provide a formal analysis of how structural perturbations propagate through the induced dynamical system. In continuous-time graph signal flows, the evolution is initialized from a random prior and driven by a time-dependent vector field parameterized by a GNN, so that graph-dependent operations are repeatedly applied along the entire trajectory. As a result, even small errors in the underlying graph can accumulate over time and affect the generated samples. Consequently, extending the stability properties of GNNs to this generative setting requires characterizing error propagation through the flow’s dynamics. Beyond pointwise stability of the vector field, we are interested in stability of the induced distributions, quantified in terms of Wasserstein distance under perturbations of the graph topology. In this context, our main contribution is to derive a Wasserstein stability bound that quantifies how relative graph perturbations affect the learnt distribution over graph signals. In support of this main result, we also: (i) formally establish that a continuous normalized flow (CNF) parametrized by a GNN vector field inherits the permutation equivariance of the base architecture; and (i) propose a regularized flow matching strategy informed by our theoretical findings, which encourages robustness to graph perturbations. Crucially, our bounds demonstrate that the stability of the vector field is intrinsically tied to the graph spectrum, enabling a fundamentally graph-aware regularization strategy. Paper outline. The remainder of the paper is organized as follows. Section I reviews the necessary preliminaries. Section I formulates the GNN-parametrized CNF, formally proves their permutation equivariance, and derives the continuous-time stability bounds. Because these continuous models must be evaluated via numerical integration, Section IV extends our theoretical framework to discrete samplers, deriving explicit stability bounds for the Euler and Heun methods. Section V details the practical implications of these theoretical results, explaining how to leverage our bounds for judicious regularization of the training loss. Section VI numerically illustrates the stability of these models on a synthetic stochastic block model (SBM) graph with smooth signals and a test case of generating fMRI signals supported on a functional brain connectome. Finally, Section VII concludes the paper by summarizing our findings and discussing limitations of our analysis, with an outlook towards future work. All technical details and proofs are deferred to the appendices. I Preliminaries This section reviews the preliminary background required to formalize the proposed generative model. We outline the key GSP concepts and survey GNN stability results, followed by the fundamentals of CNFs and flow matching. Throughout the paper, ∥⋅∥\|·\| denotes the Euclidean norm for vectors, while for matrices, ∥⋅∥2\|·\|_2 and ∥⋅∥F\|·\|_F represent the spectral and Frobenius norms, respectively. I-A Graph Signal Processing Let =(,ℰ,)G=(V,E,W) be a known undirected graph with node set V of cardinality N, edge set ℰ⊆×E ×V, and edge weight map :ℰ→ℝW:E . Associated with G is a matrix representation ∈ℝN×NS ^N× N, called a graph shift operator (GSO), which respects the sparsity structure of the graph; specifically, Sij=0S_ij=0 if (i,j)∉ℰ(i,j) for i≠ji≠ j. Common choices for the GSO include the adjacency matrix A, the graph Laplacian L, and their normalized variants. Because G is undirected then ∈S is a symmetric matrix, where =∈ℝN×N:=⊤S=\S ^N× N:S=S \. Denote the GSO eigendecomposition by =⊤S=V , where ∈ℝN×NV ^N× N contains the orthonormal eigenvectors of S and is the diagonal matrix of corresponding eigenvalues λii=1N\ _i\_i=1^N. Beyond the graph structure itself, we also consider signals defined over the nodes in V. A graph signal is a function x:→ℝx:V , mapping each node to a real value. This function can be represented as a vector ∈ℝNx ^N, where xix_i denotes the signal value at node i. The matrix ⊤V defines the graph Fourier transform (GFT), so that ¯=⊤ x=V x is the GFT of graph signal x [35]. The eigenvectors V thus serve as the Fourier basis, while the eigenvalues in are the graph frequencies. Graph convolutional filters are defined as polynomial functions of S parameterized by coefficients θpp=0P−1\ _p\_p=0^P-1. We define a graph filter operator of order P as ()≔∑p=0P−1θppH(S) _p=0^P-1 _pS^p, which acts on an input graph signal x as =()z=H(S)x [21]. Leveraging the GFT definition and the diagonal frequency response of the filter ()=∑p=0P−1θppH( )= _p=0^P-1 _p ^p, the graph convolution can be written in the spectral domain as ¯=()¯ z=H( ) x. Equivalently, the i-th spectral component satisfies z¯i=h(λi)x¯i, z_i=h( _i)\, x_i, where h(λ):=∑p=0P−1θpλph(λ):= _p=0^P-1 _pλ^p is the scalar frequency response of the graph and ()=diag(h(λ1),…,h(λN))H( )= diag(h( _1),…,h( _N)). A GNN [14, 25, 10] can be constructed from these filters via a cascade of L layers, where each layer entails a graph convolution followed by a pointwise nonlinearity σ(⋅)σ(·) such as ReLU. At layer ℓ , let ℓ∈ℝN×FℓX_ ^N× F_ denote the matrix whose columns are the graph signal features, with 0X_0 and LX_L representing the input and output signals, respectively. For each filter tap p, let ℓp∈ℝFℓ−1×Fℓ _ p ^F_ -1× F_ denote the filter coefficient matrix at layer ℓ . The output of layer ℓ in a GNN is then recursively given by ℓ=σ(∑p=0P−1pℓ−1ℓp). _ =σ ( _p=0^P-1S^pX_ -1 _ p ). (1) Henceforth assume, for simplicity, that the input and output features are one-dimensional (F0=FL=1F_0=F_L=1), while the hidden layers have dimension F1=⋯=FL−1=F_1=…=F_L-1=F. For a given input vector ∈ℝNx ^N, we set 0=X_0=x and define the overall GNN mapping as u(;)≔Lu_ θ(x;S) _L, representing the outcome of applying (1) sequentially. The GNN is parameterized by a tensor of learnable filter coefficients :=ℓpℓ,p θ:=\ _ p\_ ,p and a fixed GSO S. The latter endows the GNN with useful inductive biases about the signal’s relational structure encoded in G. I-B Stability of Graph Neural Networks Let =∈0,1N×N:=,⊤=P=\P∈\0,1\^N× N:P1=1,P 1=1\ be the set of permutation matrices, where ∈ℝN1 ^N is the all-ones vector. A fundamental property of GNNs constructed via (1) is permutation equivariance. Namely, for any permutation matrix ∈P , if we define the permuted GSO ^=⊤ S=P SP and the permuted input signal ^=⊤ x=P x, the GNN output satisfies u(^;^)=⊤u(;). u_ θ( x; S)=P u_ θ(x;S). (2) Beyond exact permutations, it is also prudent to study how the GNN behaves under structural noise. To this end, we consider the relative perturbation model from [13], where a perturbed GSO ~ S is related to the nominal GSO S modulo permutation. Specifically, we evaluate the symmetric error matrix ∈E at the permutation 0∈P_0 that affords the smallest error norm, i.e., ⋆,0 \E ,P_0\ =argmin∈,∈‖2 = argmin_E ,P \|E\|_2 (3) s.t.⊤~ .t. SP =+12(+). =S+ 12(ES+SE). Unlike absolute perturbations where edge corruptions are uniform, this relative perturbation model ensures that regions of G with weaker connectivity experience proportionally smaller changes compared to others with stronger links; see also [41, Sec. IV-A] for a detailed discussion further justifying the practical relevance of this widely adopted model. We assume that the pointwise non-linearity σ(⋅)σ(·) is normalized Lipschitz, i.e., |σ(b)−σ(a)|≤|b−a||σ(b)-σ(a)|≤|b-a|. About the graph filter, we require that the scalar frequency response h(λ)h(λ) is integral Lipschitz with constant C (see [13, Def. 4]) and it is uniformly bounded, i.e., ‖()‖2≤B\|H( )\|_2≤ B. Then, GNN outputs are stable, meaning for ‖⋆‖2≤ε\|E \|_2≤ it holds (see [13, Thm. 4]) ‖0⊤u(;~)−u(0⊤;)‖≤(Γε+(ε2))‖, \|P_0 u_ θ(x; S)-u_ θ(P_0 x;S)\|≤( +O( ^2))\|x\|, (4) for all ∈ℝNx ^N, where Γ:=C(1+δN)L(BF)L−1 :=C(1+δ N)L(BF)^L-1, and δ≤8δ≤ 8 represents the eigenvector misalignment between ⋆E and S [13, Thm. 1]. The quantity Γ elucidates the effect of GNN depth L and width F on stability, as well as graph filter spectral response summaries C and B that also play a role. We leverage this foundational bound to derive the stability of the GNN-parametrized continuous-time flow we introduce in Section I for graph signal generation. However, we note that our analytical framework is general and naturally extends to any base GNN satisfying a stability condition like (4). I-C Continuous Normalizing Flows and Flow Matching Given a dataset of signals (1),…,(n)x^(1),…,x^(n) drawn from an unknown data distribution pdatap_ data over ℝNR^N, the goal of continuous-time generative modeling is to construct a dynamical system that maps a simple prior distribution pinitp_ init to pdatap_ data. A CNF is described by an ODE for times t∈[0,1]t∈[0,1], dtt=ut(t),0∼pinit, ddtx_t=u_t θ(x_t), _0 p_ init, where the time-dependent vector field ut:ℝN×[0,1]→ℝNu_t θ:R^N×[0,1] ^N is parameterized by a neural network with learnable weights θ. We let Φt(0) _t θ(x_0) denote the continuous flow map, defined as the solution to this initial value problem at time t evaluated from 0x_0. The objective is to train the vector field such that the distribution of the endpoint 1x_1 closely approximates pdatap_ data. Let ptp_t be a target probability density path satisfying p0=pinitp_0=p_ init and p1=pdatap_1=p_ data, and let uttarget:ℝN×[0,1]→ℝNu_t^target:R^N×[0,1] ^N be the underlying target vector field that generates this path. The flow matching (FM) loss is defined as ℒFM()=t∼(0,1)∼pt[‖ut()−uttarget()‖2], _ FM( θ)=E_ subarraylt (0,1)\\ x p_t subarray [\|u_t θ(x)-u_t^target(x)\|^2 ], where (0,1)U(0,1) stands for the uniform distribution in (0,1)(0,1). In practice, both the marginal probability path ptp_t and the target vector field uttargetu_t^target are intractable. To circumvent this, [29] introduces conditional flow matching (CFM). Instead of targeting the marginal vector field directly, we define a conditional probability path pt(⋅|)p_t(·|z) for each data point ∼pdataz p_ data. This path satisfies p0(⋅|)=pinitp_0(·|z)=p_ init, while the endpoint p1(⋅|)p_1(·|z) is a distribution concentrated around the data point z, such as the multivariate Gaussian (;,σ2)N(x;z,σ^2I) for a sufficiently small σ, where I denotes the N×N× N identity matrix. If uttarget(⋅|):ℝN×[0,1]→ℝNu_t^target(·|z):R^N×[0,1] ^N is the corresponding conditional vector field that generates pt(⋅|)p_t(·|z), the CFM loss is ℒCFM()=t∼(0,1)∼pdata∼pt(⋅|)[∥ut()−uttarget(|)∥2]. _ CFM( θ)=E_ subarraylt (0,1)\\ z p_ data\\ x p_t(·|z) subarray [\|u_t θ(x)-u_t^target(x|z)\|^2 ]. (5) Crucially, it follows that ℒCFM()L_ CFM( θ) and ℒFM()L_ FM( θ) are equal up to a constant independent of θ [29], allowing us to train the model purely on the conditional paths. For example, a standard choice is the optimal transport path pt(|)=(;t,(1−t)2)p_t(x|z)=N(x;tz,(1-t)^2I). I-D Objective and Scope Our objective is to learn the distribution pdatap_ data of signals supported on the nodes of a graph G, starting from a simple prior pinitp_ init. To achieve this, we design a CNF model where the time-dependent vector field utu_t θ is parameterized by a GNN. By embedding the graph topology into the generative dynamics, we prove that the resulting flow is permutation equivariant and derive explicit bounds on its stability under structural errrors adhering to the relative perturbation model in (3). I Graph Flow Models To construct a generative flow model for a given G with GSO S, we define a time-dependent vector field ut(⋅;):ℝN×[0,1]→ℝNu_t θ(·;S):R^N×[0,1] ^N. Since the data we aim to generate are signals supported on a graph, it is natural to parameterize this continuous vector field using a GNN. Doing so allows the flow dynamics to inherently leverage the structural inductive biases of the underlying graph topology. Because a GNN natively processes only graph signals, we introduce a continuous temporal conditioning function g:ℝN×[0,1]→ℝN×dg:R^N×[0,1] ^N× d defined by the mapping (,t)↦g(,t)(x,t) g(x,t) to incorporate the time variable t∈[0,1]t∈[0,1]. We require g to be permutation equivariant, meaning that for all permutation matrices ∈P , it holds that g(^,t)=⊤g(,t). g( x,t)=P g(x,t). (6) This condition is naturally satisfied by the standard time-embedding techniques commonly used in generative models [43, 33]. For instance, given a time embedding emb:ℝ→ℝd−1emb:R ^d-1, g can be implemented via concatenation as g(,t)=[,emb(t)⊤]∈ℝN×dg(x,t)=[x,1emb(t) ] ^N× d; or via addition as g(,t)=+emb(t)∈ℝNg(x,t)=x+1emb(t) ^N using a scalar emb:ℝ→ℝemb:R . Because an identical time embedding is applied across all nodes, the operation remains entirely independent of the node indexing and (6) holds. Without loss of generality, to simplify the notational burden in the subsequent theoretical analysis, we will henceforth assume d=1d=1, such that g(,t)∈ℝNg(x,t) ^N. If a multi-dimensional mapping is used in practice, the stability bounds derived in Sections I-B and IV-B naturally extend across the feature dimensions via the triangle inequality. The vector field is then defined by evaluating a GNN on the time-augmented graph signal g(,t)g(x,t), namely ut(;)≔u(g(,t);). u_t θ(x;S) u_ θ (g(x,t);S ). (7) Given an initial condition 0∼pinitx_0 p_ init, the continuous-time dynamics of our generative model are governed by the ODE dtt=ut(t;). ddtx_t=u_t θ(x_t;S). (8) We denote by Φt(0;) _t θ(x_0;S) the associated graph-aware flow map. To guarantee the existence and uniqueness of Φt _t θ via the Picard-Lindelöf theorem, we assume the vector field is Lipschitz continuous in the state variable. This is a standard property of neural networks, including GNNs (see Proposition 3). Specifically, we assume there exists a constant M>0M>0 such that for all graph signals ,∈ℝNx,y ^N and all t∈[0,1]t∈[0,1], ‖ut(;)−ut(;)‖≤M‖−‖. \|u_t θ(x;S)-u_t θ(y;S)\|≤ M\|x-y\|. (9) I-A Permutation Equivariance of the Continuous Flow Having defined our GNN-parametrized flow model, the foremost property to verify is permutation equivariance. We first establish the permutation equivariance of the vector field. Definition 1 (Permutation equivariant vector field). We say a vector field utu_t θ is permutation equivariant if it satisfies ut(^;^)=⊤ut(;), u_t θ( x; S)=P u_t θ(x;S), for all ∈P , ∈ℝNx ^N, and t∈[0,1]t∈[0,1]. By the equivariance of the temporal conditioning g in (6) and the base GNN operator u_ θ (2), we can explicitly verify this property holds for the vector field (7) as well. For any permutation matrix ∈P , evaluating the vector field on the permuted state ^=⊤ x=P x and graph ^=⊤ S=P SP yields ut(^;^) u_t θ( x; S) =u(g(^,t);^)=u(⊤g(,t);^) =u_ θ (g( x,t); S )=u_ θ (P g(x,t); S ) (10) =⊤u(g(,t);)=⊤ut(;). =P u_ θ (g(x,t);S )=P u_t θ(x;S). Thus, the vector field satisfies Definition 1. Ultimately, our goal is to ensure that the generated signals themselves exhibit this structural symmetry. The following proposition states that indeed the resulting flow map Φt(0;) _t θ(x_0;S) inherits the permutation equivariance of the vector field. Proposition 1 (Permutation equivariance of the flow). Assume the vector field utu_t θ is permutation equivariant and Lipschitz continuous. Then, the associated flow map satisfies Φt(^0;^)=⊤Φt(0;), _t θ( x_0; S)=P _t θ(x_0;S), (11) for all t∈[0,1]t∈[0,1] and 0∈ℝNx_0 ^N. Proof. See Appendix A. In particular, evaluating (11) at t=1t=1 yields the equivariance of the generated graph signal sample. If we consider the initial state as a random vector 0∼pinitx_0 p_ init, Proposition 1 implies that the random vectors Φt(^0;^) _t θ( x_0; S) and ⊤Φt(0;)P _t θ(x_0;S) are equal almost surely. Therefore, they are equal in distribution for any pinitp_ init, which we denote as Φt(^0;^)=⊤Φt(0;). _t θ( x_0; S) D=P _t θ(x_0;S). I-B Stability of the Continuous Flow With the permutation equivariance of the flow established in Proposition 1, we can now move on to assess its structural stability. We reiterate that because node indexing is often arbitrary, comparing outputs on a nominal GSO S and a perturbed one ~ S must be carried out modulo permutation. Recall that (4) asserts the stability of the base GNN, u_ θ, under relative perturbations. To determine whether this desired property translates to the flow, we first bound the perturbation error of the vector field ut(;)=u(g(,t);)u_t θ(x;S)=u_ θ(g(x,t);S) evaluated at an arbitrary state. To this end, note that (4) holds for all ∈ℝNx ^N. Thus, we can plug the signal g(,t)g(x,t) into it, so that the perturbation error of our parameterized vector field is immediately bounded by ‖0⊤ut(;~)−ut(0⊤;)‖≤(Γε+(ε2))‖g(,t)‖, -7.0pt\|P_0 u_t θ(x; S)-u_t θ(P_0 x;S)\|≤( +O( ^2))\|g(x,t)\|, (12) for all ∈ℝNx ^N and t∈[0,1]t∈[0,1]. Because the vector field is Lipschitz continuous [cf. (9)], the flow is continuous with respect to time and initial conditions. Let t≔Φt(⊤0;)y_t^P _t θ(P x_0;S) denote the trajectory starting from a permuted initial state. Since g is continuous and the set of permutation matrices P is finite, the signal g(t,t)g(y_t^P,t) is strictly bounded across all possible permuted trajectories. This guarantees the existence of a finite constant Cg(0)≔max∈supτ∈[0,1]‖g(τ,τ)‖<∞C_g(x_0) _P _τ∈[0,1]\|g(y_τ^P,τ)\|<∞. Since the vector field is Lipschitz continuous with constant M, it is also one-sided Lipschitz with a global constant m≤Mm≤ M. We define an integrable time-dependent one-sided Lipschitz coefficient mt∈L1([0,1])m_t∈ L^1([0,1]) such that ⟨ut(;)−ut(;),−⟩≤mt‖−‖2, u_t θ(x;S)-u_t θ(y;S),x-y ≤ m_t\|x-y\|^2, (13) for all ,∈ℝNx,y ^N. Equipped with this one-sided Lipschitz condition, we can formalize the stability of the flow model. Theorem 1 (Stability of the flow model). Assuming the vector field utu_t θ is one-sided Lipschitz with coefficient mtm_t, and the base GNN, u_ θ, satisfies the stability conditions yielding (4), then it holds that min∈∥⊤Φt(0;~)−Φt(⊤0;)∥≤ΩtCg(0)(Γε+(ε2)), split _P \|P _t θ(x_0; S)&- _t θ(P x_0;S)\|\\ &≤ _tC_g(x_0) ( +O( ^2) ), split (14) for all t∈[0,1]t∈[0,1] and 0∈ℝNx_0 ^N, where Ωt _t ≔∫0texp(∫τtmss)τ, _0^t ( _τ^tm_sds )dτ, (15) Cg(0) C_g(x_0) ≔max∈supτ∈[0,1]‖g(Φτ(⊤0;),τ)‖. _P _τ∈[0,1]\|g( _τ θ(P x_0;S),τ)\|. Proof. See Appendix B. To obtain more interpretable closed-form expressions, we can evaluate the integral Ωt _t using the global constants m and M introduced previously. Corollary 1 (Simplified bounds of the flow model). Under the assumptions of Theorem 1, substituting the global Lipschitz constants simplifies the integral Ωt _t in (15) as follows: (i) For a global one-sided Lipschitz constant m, then Ωt=emt−1m _t= e^mt-1m. (i) For the standard Lipschitz constant M, then Ωt=eMt−1M _t= e^Mt-1M. Furthermore, evaluating the bound (14) of Theorem 1 at time t=1t=1 yields the stability of the generated sample: min∈∥⊤Φ1(0;~)−Φ1(⊤0;)∥≤Ω1Cg(0)(Γε+(ε2)), split _P \|P _1 θ(x_0; S)&- _1 θ(P x_0;S)\|\\ &≤ _1C_g(x_0) ( +O( ^2) ), split for all 0∈ℝNx_0 ^N. So far, our results characterize the stability of individual trajectories starting from a fixed, deterministic initial state 0x_0. However, CNFs ultimately model probability distributions. By treating the initial condition as a random vector 0∼pinitx_0 p_ init, we can broaden these pointwise results to bound the 2-Wasserstein distance [48] between the generated distributions. Corollary 2 (Wasserstein stability of the flow model). Under the assumptions of Theorem 1, let the initial state be a random vector 0∼pinitx_0 p_ init. Let ⊤Φt(0;~)∼π~tP _t θ(x_0; S) π_t^P and Φt(⊤0;)∼πt _t θ(P x_0;S) _t^P. Then, it holds that min∈W2(π~t,πt)≤ΩtCg′(Γε+(ε2)), _P W_2( π_t^P, _t^P)≤ _tC_g ( +O( ^2) ), (16) for all t∈[0,1]t∈[0,1], where Cg′=0[Cg2(0)]C_g = E_x_0 [C_g^2(x_0) ]. Proof. See Appendix C. In summary, these stability results are governed by three terms: (i) Γ , which captures the inherent stability of the base GNN; (i) Ω1 _1, which is exponential in the Lipschitz constant and dictates how the continuous dynamics amplify errors; and (i) Cg(0)C_g(x_0), the trajectory supremum across all initial condition permutations. Notably, the expectation Cg′C_g becomes independent of P when pinit()=(;,)p_ init(x)=N(x;0,I), since the initial and permuted pushforward distributions are identical. IV Discrete Graph Signal Samplers The flow map Φt _t θ associated with the ODE in (8) defines the continuous-time dynamics of our generative model. Even though this provides the exact model for transporting the prior pinitp_ init to the desired data distribution pdatap_ data, the resulting trajectories cannot be computed analytically in practice. To actually obtain a generated sample, we must approximate the continuous flow using discrete samplers, such as the Euler and Heun methods [5]. We discretize the time interval [0,1][0,1] into K steps of size h=1/Kh=1/K, corresponding to times tk=kht_k=kh. To streamline notation, we will henceforth use the iteration index k interchangeably with the time tkt_k, writing the vector field evaluated at step k simply as uk(⋅;)≔utk(⋅;)u_k θ(·;S) u_t_k θ(·;S). Under this discretization, the final iterate Kx_K approximates the terminal sample of the continuous flow, Φ1(0;) _1 θ(x_0;S). Formally, given a GSO S, a deterministic K-step sampler is defined by a sequence of update maps Tk(⋅;)k=0K−1\T_k θ(·;S)\_k=0^K-1 that generates a trajectory kk=0K\x_k\_k=0^K according to k+1=Tk(k;)x_k+1=T_k θ(x_k;S), for k=0,…,K−1k=0,…,K-1. From this, we define the associated discrete flow map Ψk(0;):ℝN→ℝN _k θ(x_0;S):R^N ^N at step k∈0,…,Kk∈\0,…,K\ recursively as: Ψ0(0;) _0 θ(x_0;S) ≔0, _0, Ψk(0;) _k θ(x_0;S) ≔Tk−1(Ψk−1(0;);),for k≥1. T_k-1 θ ( _k-1(x_0;S);S ), k≥ 1. (17) Specific numerical methods provide distinct instantiations of the maps TkT_k θ in terms of the parameterized graph-aware vector field uku_k θ. For instance, the Euler sampler yields Tk(;)=+huk(;), T_k θ(x;S)=x+h\,u θ_k(x;S), (18) while the Heun sampler is defined by the update map Tk(;)=+h2(uk(;)+uk+1(+huk(;);)). T_k θ(x;S)=x+ h2 (u_k θ(x;S)+u_k+1 θ(x+hu_k θ(x;S);S) ). (19) IV-A Permutation Equivariance of the Discrete Flow As for the continuous-time model, the foremost property to verify for any viable numerical approximation is that it preserves the structural symmetries in graph G with GSO S. Definition 2 (Permutation equivariant sampler). We say a deterministic K-step sampler is permutation equivariant if its update maps satisfy Tk(^;^)=⊤Tk(;), T_k θ( x; S)=P T_k θ(x;S), (20) for all k∈0,…,K−1k∈\0,…,K-1\, ∈ℝNx ^N, and ∈P . Proposition 2 (Permutation equivariance of the discrete flow). Assume the underlying sampler Tk(⋅;)k=0K−1\T_k θ(·;S)\_k=0^K-1 is permutation equivariant. Then, the discrete flow map in (IV) satisfies Ψk(^0;^)=⊤Ψk(0;), _k θ( x_0; S)=P _k θ(x_0;S), for all k∈0,…,Kk∈\0,…,K\, 0∈ℝNx_0 ^N and ∈P . Proof. See Appendix D. In particular, Proposition 2 guarantees the equivariance of the generated terminal sample, ΨK(^0;^)=⊤ΨK(0;) _K θ( x_0; S)=P _K θ(x_0;S). This condition is naturally satisfied by many widely used discrete samplers. Notably, if the underlying vector field uku_k θ is permutation equivariant (cf. Definition 1), the Euler and Heun samplers defined in (18) and (19) are permutation equivariant as well. And by Proposition 2, the discrete flow maps induced by both methods are permutation equivariant. The full details to establish the equivariance of these samplers are provided in Appendices E and F. IV-B Stability of the Discrete Flow Building on the continuous case in Section I-B, our goal is to bound the deviation between a trajectory generated on the perturbed GSO ~ S, where ~ S is defined as in (3), and the appropriately permuted trajectory on the nominal GSO S. Arguably, this is the most important stability consideration when it comes to generating graph signals in practice. By expanding the recursive definition of the discrete flow, the error at any step k+1k+1 can be decomposed using the triangle inequality into two distinct sources of error. This motivates two key characterizations of the numerical method: (i) how robust the sampler’s update map is to graph perturbations (graph stability); and (i) how sensitive it is to input variations (state stability). Definition 3 (Graph stable sampler). We say a permutation equivariant sampler Tk(⋅;)k=0K−1\T_k θ(·;S)\_k=0^K-1 is graph stable if there exists a permutation invariant bounding function β:ℝN×[0,1]→ℝ≥0β:R^N×[0,1] _≥ 0 (i.e., β(⊤,t)=β(,t)β(P x,t)=β(x,t) for all ∈P ) and a constant Γ>0 >0, such that ‖0⊤Tk(;~)−Tk(0⊤;)‖≤β(,tk)(Γε+(ε2)), \|P_0 T_k θ(x; S)-T_k θ(P_0 x;S)\|≤β(x,t_k) ( +O( ^2) ), for all ∈ℝNx ^N and k∈0,…,K−1k∈\0,…,K-1\. Definition 4 (State stable sampler). We say a sampler Tk(⋅;)k=0K−1\T_k θ(·;S)\_k=0^K-1 is state stable if there exists a constant α≥0α≥ 0, such that ‖Tk(;)−Tk(;)‖≤α‖−‖, \|T_k θ(x;S)-T_k θ(y;S)\|≤α\|x-y\|, for all ,∈ℝNx,y ^N and k∈0,…,K−1k∈\0,…,K-1\. Equipped with these two characterizations, we can establish a general stability bound for the resulting discrete flow map by tracking the accumulated error across the K integration steps. Theorem 2 (Stability of the discrete flow). Assuming the underlying permutation equivariant sampler Tk(⋅;)k=0K−1\T_k θ(·;S)\_k=0^K-1 is graph stable with bounding function β, and state stable with constant α. Then, it holds that min∈∥⊤Ψk(0;~) _P \|P _k θ(x_0; S) −Ψk(⊤0;)∥ - _k θ(P x_0;S)\| ≤Υk(α)βmax(Γε+(ε2)), ≤ _k(α)\, _ ( +O( ^2) ), for all k∈1,…,Kk∈\1,…,K\ and 0∈ℝNx_0 ^N, where Υk(α) _k(α) ≔1−αk1−αif α≠1,kif α=1, cases 1-α^k1-α&if α≠ 1,\\ k&if α=1, cases βmax _ ≔max∈maxj∈0,…,K−1β(j,tj), _P _j∈\0,…,K-1\β(y_j^P,t_j), j _j^P ≔Ψj(⊤0;). _j θ(P x_0;S). Proof. See Appendix G. As established previously, both the Euler and Heun methods are permutation equivariant, allowing us to directly apply this general stability framework. By evaluating these specific numerical methods under the standard assumption of Lipschitz continuity, we obtain the following explicit stability bounds. Corollary 3 (Stability of Euler and Heun samplers). Under the assumptions of Theorem 1, if the vector field utu_t θ is Lipschitz continuous with constant M, then the discrete flow maps induced by the respective samplers satisfy the following bounds for all k∈1,…,Kk∈\1,…,K\ and 0∈ℝNx_0 ^N: (i) Euler sampler: min∈‖⊤Ψk(0;~)−Ψk(⊤0;)‖≤(1+hM)k−1MCg(0)(Γε+(ε2)), split _P &\|P _k θ(x_0; S)- _k θ(P x_0;S)\|\\ &≤ (1+hM)^k-1MC_g(x_0) ( +O( ^2) ), split where Cg(0) C_g(x_0) ≔max∈maxj<k‖g(j,tj)‖. _P _j<k\|g(y_j^P,t_j)\|. (i) Heun sampler: min∈‖⊤Ψk(0;~)−Ψk(⊤0;)‖≤(1+hM+h2M22)k−1MCg(0)(Γε+(ε2)), split _P &\|P _k θ(x_0; S)- _k θ(P x_0;S)\|\\ &≤ (1+hM+ h^2M^22 )^k-1MC_g(x_0) ( +O( ^2) ), split where Cg(0)≔max∈maxj<k C_g(x_0) _P _j<k max(∥g(j,tj)∥, (\|g(y_j^P,t_j)\|, ∥g(j+huj(j;),tj+1)∥). \|g(y_j^P+hu_j θ(y_j^P;S),t_j+1)\| ). Both bounds are evaluated along the unperturbed trajectories defined by j≔Ψj(⊤0;)y_j^P _j θ(P x_0;S). Proof. See Appendix H. A few noteworthy observations can be made regarding these discrete stability bounds. First, their structure perfectly mirrors the bounds established in Theorem 1, consisting of a trajectory dependent supremum scaled by a time dependent growth factor. Notably, in the limit as the step size h→0h→ 0 (and the number of steps k→∞k→∞ such that kh=tkh=t), the growth factors for both Euler and Heun samplers converge to their continuous-time counterpart [cf. (i) in Corollary 1]. Second, just as we did for the continuous flow, we can evaluate these discrete bounds at the final integration step k=Kk=K (which corresponds to time t=1t=1). This yields the explicit stability bound for the generated sample. Again, our ultimate robustness objective is to measure the stability of the generated distributions. By treating the initial condition as a random vector 0∼pinitx_0 p_ init, we can lift the pointwise guarantees of our numerical samplers to a Wasserstein stability bound for the discrete flow. Corollary 4 (Wasserstein stability of the discrete flow). Under the assumptions of Theorem 2, let the initial state be a random vector 0∼pinitx_0 p_ init. Let the resulting distributions at step k be denoted by ⊤Ψk(0;~)∼π~kP _k θ(x_0; S) π_k^P and Ψk(⊤0;)∼πk _k θ(P x_0;S) _k^P. Then, it holds that min∈W2(π~k,πk)≤Υk(α)βmax′(Γε+(ε2)), _P W_2( π_k^P, _k^P)≤ _k(α)\,β _ ( +O( ^2) ), (21) for all k∈1,…,Kk∈\1,…,K\, where βmax′≔0[βmax2(0)]β _ E_x_0 [ _ ^2(x_0) ]. Proof. See Appendix I. V Implications for Training Regularization Beyond providing theoretical guarantees, the stability bounds derived in the previous sections offer actionable information on how to regularize our generative model to promote robustness. By inspecting these expressions, we can isolate the factors inherent to the data from those governed by our design choices. The graph size N is fixed, and the eigenvector misalignment δ depends entirely on the perturbed GSO ~ S. On the other hand, the GNN parameters L, F, B, C, and M are determined by our architecture and learned filter weights. The number of layers L (GNN depth) and hidden features F (GNN width) dictate the architectural capacity of the vector field. As evidenced by the factor (BF)L−1(BF)^L-1, increasing depth and width exponentially amplifies instability [13]. As expected, this implies the architecture should be kept as simple as possible. While acknowledging this structural trade-off is standard practice, the true potential for principled regularization lies in actively controlling the filter-dependent constants B, C, and M during training. Recall that the filter’s frequency response is bounded such that ‖()‖2≤B\|H( )\|_2≤ B. Expressed in terms of the filter coefficients ℓp _ p, the infimum B that satisfies this condition across all layers ℓ=1,…,L =1,…,L is given by B=maxℓmaxi∈1,…,n‖∑p=0P−1ℓpλip‖2. B= _ _i∈\1,…,n\ \| _p=0^P-1 _ p _i^p \|_2. (22) Similarly, the integral Lipschitz condition implies that |λdh(λ)dλ|≤C |λ dh(λ)dλ |≤ C. Expressed in terms of the filter coefficients, the infimum C satisfying this inequality is bounded by C≤maxℓmaxi∈1,…,n‖∑p=0P−1pℓpλip‖2. C≤ _ _i∈\1,…,n\ \| _p=0^P-1p _ p _i^p \|_2. (23) Beyond bounding the filter conditions of the base GNN, next we establish an explicit bound for the spatial Lipschitz constant M of the vector field introduced in (9). Proposition 3 (Lipschitz constant of the vector field). The infimum spatial Lipschitz constant M of the vector field (7) is bounded by M≤∏ℓ=1Lmaxi∈1,…,n‖∑p=0P−1ℓpλip‖2.M≤ _ =1^L _i∈\1,…,n\ \| _p=0^P-1 _ p _i^p \|_2. (24) Proof. See Appendix J. We note that if the temporal conditioning function g decouples the state tx_t and time t, for example, via concatenation g(,t)=[,emb(t)⊤]∈ℝN×dg(x,t)=[x,1emb(t) ] ^N× d, the time component does not affect the spatial Lipschitz constant. This follows as the norm of the difference perfectly isolates the state ‖g(,t)−g(,t)‖F=‖−‖.\|g(x,t)-g(y,t)\|_F=\|x-y\|. Consequently, the weights of the first layer operating exclusively on the time embedding do not contribute to the bound in Proposition 3. All in all, the bounds in (22), (23) and (24), provide a principled basis to regularize B, C, and M, respectively. A key observation is that these bounds depend directly on the spectrum λii=1N\ _i\_i=1^N of the given GSO S, making them fundamentally graph-aware. Furthermore, since the underlying graph topology is fixed, this spectrum only needs to be computed once prior to training. Related works have explored regularization strategies based on similar quantities for static GNNs [2, 6]. However, highlighting the specific theoretical results of our work, the derived bounds suggest that the spatial Lipschitz constant M acts as the dominating term since it drives the exponential growth factor eMM e^MM (see Corollary 1). Consequently, for the scope of this paper, our proposed regularization strategy solely penalizes M to ensure stable generative dynamics across CNF training paradigms, though alternative approaches could be chosen to jointly regularize both the spatial and integral Lipschitz constants. To demonstrate these ideas, we consider a regularized training objective of the form ℒ()=ℒCFM()+μℛ(), ( θ)=L_ CFM( θ)+ ( θ), (25) where the ℒCFML_CFM is defined as in (5) and the regularization term ℛR is defined as the bound in Proposition 3, while μ>0μ>0 is a tuning hyperparameter. VI Numerical Experiments To illustrate the derived stability bounds in a practical setting, we consider two test cases: (i) a synthetic setting based on a SBM graph with smooth signals; and (i) a real-world example using fMRI data [20], where brain connectivity defines the underlying graph structure. In both cases, given the respective graphs the goal is to generate signals that follow the data distribution. For the SBM case, we introduce a controlled synthetic perturbation governed by a parameter ε . For the fMRI case, perturbations arise naturally from errors in estimating the graph structure, specifically from the empirical correlation matrix used to construct the functional connectivity graph. The main objective of this numerical evaluation is to illustrate how loose the derived stability bounds can be in practice and to show that, despite this looseness, the proposed bound-informed regularization leads to more stable models in the presence of graph perturbations. Figure 1: Performance and stability on a SBM graph with smooth signals under a controlled synthetic perturbation of magnitude ε . The aMMD (left) shows that Regularized Flow Matching (RFM) consistently outperforms standard Flow Matching (FM) across all ε values. The empirical stability (center), measured by the change in the generated output signals due to ε , demonstrates that RFM is more stable than FM. The theoretical stability bounds (right) correctly bound the empirical errors but are not tight. VI-A Experimental Setup Datasets. In the synthetic setting, we use an SBM graph [17] with two communities of 10 nodes each (N=20N=20). Signals in each community are drawn from a Gaussian distribution with standard deviation 1. One community has mean 1, while the other has mean −1-1. We generate 500 training, 500 validation, and 500 test graph signals. To induce smoothness, the signals are processed using a low-pass graph filter. For the fMRI data, we consider a single subject from the HCP dataset, with data in ℝ360×1190R^360× 1190, where 360 represents brain regions (nodes) and 1190 is the number of time points. We treat each time point as an individual graph signal, i.e., a vector (i)∈ℝ360x^(i) ^360. Thus, the dataset consists of (1),…,(1190)\x^(1),…,x^(1190)\. This provides a well-motivated practical scenario, as fMRI signals are known to inherently exhibit graph signal smoothness over functional brain connectomes [19, 36]. We first split the dataset into 80% training/validation data and 20% test data. Subsequently, 20% of the training/validation set is held out for validation, while the remaining 80% is used for training. Architecture. We parameterize the continuous vector field ut(t;)u_t θ(x_t;S) defined in (7) using a two-layer GNN (L=2L=2). The temporal conditioning function is defined as g(,t)=[,emb(t)⊤]g(x,t)=[x,1emb(t) ], where emb(t)∈ℝ64emb(t) ^64 is a standard 64-dimensional sinusoidal time embedding. Because the input consists of the 1-dimensional node feature and the 64-dimensional time embedding, the input feature dimension is F0=65F_0=65. The GNN hidden layer has a width of F1=4F_1=4 features, and the output layer F2=1F_2=1. Both graph convolutional layers utilize P=4P=4 filter taps and use the normalized graph Laplacian as the GSO. We employ a normalized SiLU [12] as the pointwise nonlinearity σ(⋅)σ(·). While we opt for this relatively simple architecture to balance expressiveness and stability, we have empirically observed consistent results across a range of alternative architectural configurations. We generate samples using an Euler sampler with K=100K=100 steps. Loss functions. We consider two training alternatives. For the first case, standard Flow Matching (FM), we directly minimize the conditional flow matching loss in (5) using the optimal transport path [29], defined as pt(|)=(;t,(1−t)2)p_t(x|z)=N(x;tz,(1-t)^2I). For the second case, Regularized Flow Matching (RFM), we use the new regularized loss proposed in (25), which combines (5) with the regularization term introduced in Proposition 3, with μ=0.01μ=0.01 selected via a grid search. Training and evaluation. We use an ADAM optimizer with learning rate 0.0010.001. We train for 500 epochs with early stopping of 50 epochs and a batch size of 64. In all experiments, the model is trained using the unperturbed graph, i.e., ε=0 =0, and for the fMRI setting the graph is constructed using the full training and validation data. At test time, we evaluate robustness by introducing graph perturbations. In both cases, we run 10 random datasets; for the SBM this corresponds to drawing new samples, while for the fMRI case we generate new train/validation/test partitions. For each dataset draw, we sample 10 perturbed GSOs ~ S. We report two evaluation metrics: (i) a performance-oriented metric based on the average Maximum Mean Discrepancy (aMMD) [15]; and (i) the variation in generated outputs induced by perturbations of the graph support, quantified by ‖Ψk(0;~)−Ψk(0;)‖\| _k θ(x_0; S)- _k θ(x_0;S)\|. We follow the evaluation protocol of prior work for aMMD [39]. Specifically, we compute MMD between generated and test distributions over three signal descriptors—quadratic variation, spectral centroid, and degree correlation—capturing smoothness, spectral content, and graph alignment, respectively, and average the resulting MMDs. In all resulting plots, the curves represent the median across these independent runs, and the shaded regions indicate the 25th and 75th percentiles. Synthetic perturbation. At test time, we evaluate robustness under controlled perturbations of the graph structure in the SBM setting. Specifically, we generate a random diagonal matrix E such that ‖2≤ε\|E\|_2≤ . Following the framework of [13], the entries of E are drawn uniformly from the interval [(1−ε)ε,ε][(1- ) , ]. The perturbed GSO is then defined as ~=++ S=S+ES+SE. We simulate logarithmically-spaced ε values between 10−310^-3 and 10−0.510^-0.5. Functional connectivity perturbation. At test time, we further evaluate robustness under data-driven perturbations of the graph in the fMRI setting. We construct the GSO ~ S from an empirical Pearson correlation matrix estimated on a subset of the available data. Specifically, we vary the fraction of samples used to estimate the correlation network from the training and validation sets, ranging from 1%1\% to 90%90\% in logarithmic scale. VI-B Results and Discussion Synthetic perturbation on SBM graphs. Fig. 1 summarizes the results for the synthetic SBM experiment. In terms of generative quality (Fig. 1, left), the RFM model achieves a lower aMMD than standard FM across all noise levels ε . The empirical stability evaluation (Fig. 1, center) demonstrates that RFM is consistently more stable than FM. By actively penalizing the spatial Lipschitz constant during training, RFM effectively limits the accumulation of structural errors along the trajectory. When comparing the empirical output variations against the theoretical bounds (Fig. 1, right), we observe that while the bounds correctly limit the maximum error, they are not tight. This looseness is expected and stems from two primary sources. First, computing the exact Lipschitz constant of a neural network is NP-hard [42], so we rely on the upper bound established in Proposition 3. This approximation is inherently loose because it is a conservative upper bound, and because it bounds the global Lipschitz constant M instead of the tighter one-sided Lipschitz constant m. Substituting this worst-case approximation into Theorem 1 yields a pessimistic estimate for the exponential growth factor Ω1 _1. Second, the stability bound of the base GNN in (4), as established in [13], is already a conservative upper bound. This native lack of tightness naturally propagates through our flow bounds. Realistic perturbation on fMRI data. Fig. 2 depicts the fMRI experiments, reporting how these models behave under real-world, data-driven perturbations. As the ratio of training samples decreases, the empirical correlation matrix becomes a noisier estimate of the brain connectome used to train the model. For high sample ratios (Fig. 2, left), RFM and FM yield comparable aMMD, indicating that the regularization does not degrade performance when the graph structure is reliable. For low sample ratios, RFM achieves better aMMD. The empirical divergence (Fig. 2, center) shows that RFM is consistently more stable than FM when the graph estimate is degraded. The bounds for the fMRI test case (Fig. 2, right) exhibit the same conservative behavior discussed previously, but they are notably looser than in the SBM case. This increased looseness is driven by two factors: the brain connectome has a much higher number of nodes, and the complexity of the fMRI signals requires the learned models (both FM and RFM) to have a larger spatial Lipschitz constant M. Nevertheless, the same overall general conclusions hold. Crucially, despite the looseness of the theoretical bounds, explicitly controlling the architectural constants derived from them yields a model that is more robust to graph perturbations. Figure 2: Performance and stability on fMRI signals over a functional brain connectome under data-driven graph perturbations. The x-axis represents the ratio of training samples used to estimate the empirical correlation matrix. The aMMD (left) shows that Regularized Flow Matching (RFM) performs comparably to standard Flow Matching (FM) at high sample ratios, but achieves better performance at lower ratios. The empirical stability (center), measured by the change in the generated output signals due to graph estimation errors, demonstrates that RFM is more stable than FM. The theoretical stability bounds (right) correctly bound the empirical errors but are not tight. VII Conclusion In this paper, we established the equivariance and stability properties relevant to the generation of graph signals using GNN-parametrized CNFs. We proved that both the continuous-time ODEs and their discrete numerical approximations, specifically Euler and Heun samplers, preserve permutation equivariance inherited from the underlying GNN. By analyzing how structural errors propagate through the generative flow, we derived explicit Wasserstein stability bounds for the generated distributions under relative graph perturbations. Beyond theoretical guarantees, these bounds informed a practical training strategy. Having identified the spatial Lipschitz constant of the vector field as the main source of error accumulation, we introduced RFM to actively penalize this constant. Our experiments on synthetic SBMs and empirical fMRI data confirmed that explicitly controlling this constant yields generative models that are more robust to structural noise, without sacrificing the quality of the generated signals. While our practical implementation focused specifically on regularizing the spatial Lipschitz constant, exploring the theoretical and empirical compromises of jointly regularizing the spatial and frequency-response Lipschitz constants represents one direction to further optimize the trade-off between model expressiveness and structural robustness. The analytical framework presented here establishes a foundation for evaluating broader classes of continuous generative models. Because the continuous-time and discrete stability guarantees are grounded in the properties of the parameterized vector field, this analysis can be extended to stochastic differential equations, provided the injected noise satisfies permutation equivariance. Furthermore, our results provide the basis to establish transferability bounds, addressing the problem of training a generative model on smaller (sampled) graphs and transferring them to large-scale graphs [7]. Appendix A Permutation Equivariance of the Continuous Flow Proof of Proposition 1. Let t=Φt(0;)x_t= _t θ(x_0;S) and define t≔⊤ty_t x_t. Differentiating with respect to t and using the equivariance of the vector field utu_t θ in (10) yields dtt ddty_t =⊤dtt=⊤ut(t;) =P ddtx_t=P u_t θ(x_t;S) =ut(⊤t;⊤)=ut(t;^). =u_t θ(P x_t;P SP)=u_t θ(y_t; S). Since 0=⊤0=^0y_0=P x_0= x_0, the path ty_t satisfies the initial value problem defining Φt(^0;^) _t θ( x_0; S). By the uniqueness of solutions guaranteed by the Picard-Lindelöf theorem, it follows that t=Φt(^0;^)y_t= _t θ( x_0; S), which concludes the proof. ∎ Appendix B Stability of the Flow Model Proof of Theorem 1. Let 0∈P_0 be the permutation matrix from Section I-B that defines the relative perturbation in (3). Let ~t≔Φt(0;~) x_t _t θ(x_0; S) and t≔Φt(0⊤0;)x_t _t θ(P_0 x_0;S). First, we upper bound the minimum over all permutations by evaluating the error specifically at 0P_0, to obtain min∈ _P ‖⊤Φt(0;~)−Φt(⊤0;)‖ \|P _t θ(x_0; S)- _t θ(P x_0;S)\| ≤‖0⊤Φt(0;~)−Φt(0⊤0;)‖ ≤\|P_0 _t θ(x_0; S)- _t θ(P_0 x_0;S)\| =‖0⊤~t−t‖. =\|P_0 x_t-x_t\|. Let t≔0⊤~t−tz_t _0 x_t-x_t. Differentiating yields 12dt‖t‖2 12 ddt\|z_t\|^2 =12dt⟨t,t⟩=⟨t,dtt⟩ = 12 ddt _t,z_t = _t, ddtz_t =⟨0⊤~t−t,0⊤ut(~t;~)−ut(t;)⟩. = _0 x_t-x_t,P_0 u_t θ( x_t; S)-u_t θ(x_t;S) . By adding and subtracting the term 0⊤ut(0t;~)P_0 u_t θ(P_0x_t; S), we can split the inner product 12dt‖t‖2 12 ddt\|z_t\|^2 =⟨0⊤~t−t,0⊤(ut(~t;~)−ut(0t;~)) = _0 x_t-x_t,P_0 (u_t θ( x_t; S)-u_t θ(P_0x_t; S) ) +0⊤ut(0t;~)−ut(t;)⟩ +P_0 u_t θ(P_0x_t; S)-u_t θ(x_t;S) =⟨0⊤~t−t,0⊤(ut(~t;~)−ut(0t;~))⟩ = _0 x_t-x_t,P_0 (u_t θ( x_t; S)-u_t θ(P_0x_t; S) ) +⟨0⊤~t−t,0⊤ut(0t;~)−ut(t;)⟩. + _0 x_t-x_t,P_0 u_t θ(P_0x_t; S)-u_t θ(x_t;S) . To bound this expression, we handle each term separately. For the first term, we use the fact that the permutation matrix 0P_0 is orthogonal (i.e., ⟨0,0⟩=⟨,⟩ _0a,P_0b = ,b ) and apply the one-sided Lipschitz condition (13). For the second term, we apply the Cauchy-Schwarz inequality and we use the stability of the vector field by substituting the input signal with 0tP_0x_t in (12). This yields 12dt‖t‖2 12 ddt\|z_t\|^2 =⟨~t−0t,ut(~t;~)−ut(0t;~)⟩ = x_t-P_0x_t,u_t θ( x_t; S)-u_t θ(P_0x_t; S) +⟨t,0⊤ut(0t;~)−ut(t;)⟩ + _t,P_0 u_t θ(P_0x_t; S)-u_t θ(x_t;S) ≤mt‖~t−0t‖2 ≤ m_t\| x_t-P_0x_t\|^2 +‖t‖(Γε+(ε2))‖g(t,t)‖. +\|z_t\| ( +O( ^2) )\|g(x_t,t)\|. Since ‖~t−0t‖=‖0⊤~t−t‖=‖t‖\| x_t-P_0x_t\|=\|P_0 x_t-x_t\|=\|z_t\|, the inequality simplifies to 12dt‖t‖2≤mt‖t‖2+‖t‖(Γε+(ε2))‖g(t,t)‖. 12 ddt\|z_t\|^2≤ m_t\|z_t\|^2+\|z_t\| ( +O( ^2) )\|g(x_t,t)\|. By the definition of Cg(0)C_g(x_0), which takes the supremum over the time interval and the maximum over all permutation matrices in P, it directly follows that ‖g(t,t)‖≤Cg(0)\|g(x_t,t)\|≤ C_g(x_0), and hence 12dt‖t‖2≤mt‖t‖2+‖t‖(Γε+(ε2))Cg(0). 12 ddt\|z_t\|^2≤ m_t\|z_t\|^2+\|z_t\| ( +O( ^2) )C_g(x_0). Let rt=‖t‖r_t=\|z_t\|. From the chain rule 12dtrt2=rtdtrt 12 ddtr_t^2=r_t ddtr_t, and we thus obtain rtdtrt≤mtrt2+rt(Γε+(ε2))Cg(0). r_t ddtr_t≤ m_tr_t^2+r_t ( +O( ^2) )C_g(x_0). For t such that rt=0r_t=0, the result is trivially true. For rt>0r_t>0, we can divide by rtr_t to obtain dtrt≤mtrt+(Γε+(ε2))Cg(0). ddtr_t≤ m_tr_t+ ( +O( ^2) )C_g(x_0). Using Grönwall’s inequality and the fact that the initial error r0=‖0⊤0−0⊤0‖=0r_0=\|P_0 x_0-P_0 x_0\|=0, we have rt≤Cg(0)(Γε+(ε2))(∫0texp(∫τtmss)τ), r_t≤ C_g(x_0) ( +O( ^2) ) ( _0^texp ( _τ^tm_sds )dτ ), which concludes the proof. ∎ Appendix C Wasserstein Stability of the Flow Model Proof of Corollary 2. Let 0∈P_0 be the specific permutation matrix from Section I-B that minimizes the relative perturbation error in (3). Evaluating the flow on the exact same initial random variable 0x_0 under both the nominal graph S and the perturbed graph ~ S forms a valid coupling. Thus, we can upper bound the Wasserstein-2 distance using the bound derived in Theorem 1 min∈W2(π~t,πt)2 _P W_2( π_t^P, _t^P)^2 ≤0[‖0⊤Φt(0;~)−Φt(0⊤0;)‖2] _x_0 [\|P _0 _t θ(x_0; S)- _t θ(P _0x_0;S)\|^2 ] ≤0[(ΩtCg(0)(Γε+(ε2)))2]. _x_0 [ ( _tC_g(x_0) ( +O( ^2) ) )^2 ]. Because Ωt _t and (Γε+(ε2)) ( +O( ^2) ) are deterministic and independent of 0x_0, we can factor them out of the expectation 0 _x_0 [(ΩtCg(0)(Γε+(ε2)))2] [ ( _tC_g(x_0) ( +O( ^2) ) )^2 ] =Ωt2(Γε+(ε2))20[Cg(0)2] = _t^2 ( +O( ^2) )^2E_x_0 [C_g(x_0)^2 ] =Ωt2(Γε+(ε2))2(Cg′)2. = _t^2 ( +O( ^2) )^2(C_g )^2. Taking the square root of both sides directly yields (16). ∎ Appendix D Permutation Equivariance of the Discrete Flow Proof of Proposition 2. We proceed by induction on k. For the base case k=0k=0, the definition of the discrete flow in (IV) and our initial condition yields Ψ0(^0;^)=^0=⊤0=⊤Ψ0(0;). _0 θ( x_0; S)= x_0=P x_0=P _0 θ(x_0;S). Assume the claim holds for some k∈0,…,K−1k∈\0,…,K-1\, meaning Ψk(^0;^)=⊤Ψk(0;) _k θ( x_0; S)=P _k θ(x_0;S). For step k+1k+1, we apply the recursive definition of the discrete flow and the permutation equivariance of the update map TkT_k in (20) to obtain Ψk+1(^0;^) _k+1 θ( x_0; S) =Tkθ(Ψk(^0;^);^) =T_k^θ ( _k θ( x_0; S); S ) =Tk(⊤Ψk(0;);^) =T_k θ (P _k θ(x_0;S); S ) =⊤Tk(Ψk(0;);) =P T_k θ ( _k θ(x_0;S);S ) =⊤Ψk+1(0;). =P _k+1 θ(x_0;S). By induction, the property holds for all k∈0,…,Kk∈\0,…,K\. ∎ Appendix E Permutation Equivariance of the Euler Sampler Proof. Assuming the vector field uku_k θ is permutation equivariant, we evaluate the Euler update map (18) on the permuted state ^=⊤ x=P x and permuted graph ^=⊤ S=P SP. This yields Tk(^;^) T_k θ( x; S) =⊤+huk(⊤;⊤) =P x+h\,u θ_k(P x;P SP) =⊤+h⊤uk(;) =P x+h\,P u θ_k(x;S) =⊤(+huk(;)) =P (x+h\,u θ_k(x;S) ) =⊤Tk(;). =P T_k θ(x;S). So, the Euler sampler is permutation equivariant. ∎ Appendix F Permutation Equivariance of the Heun Sampler Proof. Assuming the vector field uku_k θ is permutation equivariant, we evaluate the Heun update map (19) on the permuted inputs ^=⊤ x=P x and ^=⊤ S=P SP. Using the equivariance property on both evaluations of the vector field yields Tk(^;^)=^+h2(uk(^;^)+uk+1(^+huk(^;^);^)) T_k θ( x; S)= x+ h2 (u_k θ( x; S)+u_k+1 θ ( x+hu_k θ( x; S); S ) ) =⊤+h2(⊤uk(;)+uk+1(⊤(+huk(;));^)) =P x+ h2 (P u_k θ(x;S)+u_k+1 θ (P (x+hu_k θ(x;S)); S ) ) =⊤+h2(⊤uk(;)+⊤uk+1(+huk(;);)) =P x+ h2 (P u_k θ(x;S)+P u_k+1 θ (x+hu_k θ(x;S);S ) ) =⊤(+h2(uk(;)+uk+1(+huk(;);))) =P (x+ h2 (u_k θ(x;S)+u_k+1 θ (x+hu_k θ(x;S);S ) ) ) =⊤Tk(;). =P T_k θ(x;S). Thus, the Heun sampler is permutation equivariant as well. ∎ Appendix G Stability of the Discrete Flow Proof of Theorem 2. Let 0∈P_0 be the specific permutation matrix from Section I-B that minimizes the relative perturbation error in (3). Define ~k≔Ψk(0;~) x_k _k θ(x_0; S) and k≔Ψk(0⊤0;)x_k _k θ(P_0 x_0;S). By definition of the discrete flow min∈ _P ‖⊤Ψk+1(0;~)−Ψk+1(⊤0;)‖ \|P _k+1 θ(x_0; S)- _k+1 θ(P x_0;S)\| ≤‖0⊤Ψk+1(0;~)−Ψk+1(0⊤0;)‖ ≤\|P_0 _k+1 θ(x_0; S)- _k+1 θ(P_0 x_0;S)\| =‖0⊤~k+1−k+1‖ =\|P_0 x_k+1-x_k+1\| =‖0⊤Tk(~k;~)−Tk(k;)‖. =\|P_0 T_k θ( x_k; S)-T_k θ(x_k;S)\|. By adding and subtracting the term 0⊤Tk(0k;~)P_0 T_k θ(P_0x_k; S), applying the triangle inequality, and using the orthogonality of the permutation matrix 0⊤P_0 , we can bound the error as ∥0⊤Tk(~k;~) \|P_0 T_k θ( x_k; S) −Tk(k;)∥ -T_k θ(x_k;S)\| ≤‖0⊤Tk(~k;~)−0⊤Tk(0k;~)‖ ≤\|P_0 T_k θ( x_k; S)-P_0 T_k θ(P_0x_k; S)\| +‖0⊤Tk(0k;~)−Tk(k;)‖ +\|P_0 T_k θ(P_0x_k; S)-T_k θ(x_k;S)\| =‖Tk(~k;~)−Tk(0k;~)‖ =\|T_k θ( x_k; S)-T_k θ(P_0x_k; S)\| +‖0⊤Tk(0k;~)−Tk(k;)‖. +\|P_0 T_k θ(P_0x_k; S)-T_k θ(x_k;S)\|. Because the sampler is state stable (Definition 4), we have ‖Tk(~k;~)−Tk(0k;~)‖ \|T_k θ( x_k; S)-T_k θ(P_0x_k; S)\| ≤α‖~k−0k‖ ≤α\| x_k-P_0x_k\| =α‖0⊤~k−k‖. =α\|P_0 x_k-x_k\|. By substituting 0kP_0x_k into the definition of a graph stable sampler (Definition 3) and using the permutation equivariance of β (i.e., β(0k,tk)=β(k,tk)β(P_0x_k,t_k)=β(x_k,t_k)), yields ‖0⊤Tk(0k;~)−Tk(k;)‖≤β(k,tk)(Γε+(ε2)). \|P_0 T_k θ(P_0x_k; S)-T_k θ(x_k;S)\|≤β(x_k,t_k) ( +O( ^2) ). By definition, k=Ψk(0⊤0;)x_k= _k θ(P_0 x_0;S) is the unperturbed trajectory starting from the permuted initial condition 0⊤0P_0 x_0. Because βmax _ is defined as the maximum over all ∈P , it directly follows that β(k,tk)≤βmaxβ(x_k,t_k)≤ _ . Thus, the error bound satisfies the recurrence relation ‖0⊤~k+1−k+1‖≤α‖0⊤~k−k‖+βmax(Γε+(ε2)). \|P_0 x_k+1-x_k+1\|≤α\|P_0 x_k-x_k\|+ _ ( +O( ^2) ). Since ‖0⊤~0−0‖=0\|P_0 x_0-x_0\|=0, iterating this recurrence relation yields a geometric series ‖0⊤~k−k‖ \|P_0 x_k-x_k\| ≤βmax(Γε+(ε2))∑i=0k−1αi ≤ _ ( +O( ^2) ) _i=0^k-1α^i =1−αk1−αβmax(Γε+(ε2))if α≠1,kβmax(Γε+(ε2))if α=1, = which concludes the proof. ∎ Appendix H Stability of Euler and Heun Samplers Proof of Corollary 3. As established in Appendices E and F, both the Euler and Heun samplers are permutation equivariant. We will proceed by showing that each sampler is state stable (Definition 4) with a specific constant α, and graph stable (Definition 3) with a bounding function β. Then, applying Theorem 2, we will arrive at the respective bounds. Euler sampler. By the triangle inequality and Lipschitz continuity of the vector field with constant M, we have ∥Tk(;) \|T_k θ(x;S) −Tk(;)∥ -T_k θ(y;S)\| =‖−+h(uk(;)−uk(;))‖ =\|x-y+h (u_k θ(x;S)-u_k θ(y;S) )\| ≤‖−‖+h‖uk(;)−uk(;)‖ ≤\|x-y\|+h\|u_k θ(x;S)-u_k θ(y;S)\| ≤‖−‖+hM‖−‖ ≤\|x-y\|+hM\|x-y\| =(1+hM)‖−‖. =(1+hM)\|x-y\|. Thus, the Euler sampler is state stable with α=1+hMα=1+hM. For graph stability, we evaluate the error between the perturbed and nominal update maps using the stability of the continuous vector field (12), which yields ∥0⊤Tk(;~) \|P_0 T_k θ(x; S) −Tk(0⊤;)∥ -T_k θ(P_0 x;S)\| =h‖0⊤uk(;~)−uk(0⊤;)‖ =h\|P_0 u_k θ(x; S)-u_k θ(P_0 x;S)\| ≤h‖g(,tk)‖(Γε+(ε2)). ≤ h\|g(x,t_k)\| ( +O( ^2) ). We identify and define the bounding function β(,tk)≔h‖g(,tk)‖β(x,t_k) h\|g(x,t_k)\|. To verify that β is permutation invariant, we use the permutation equivariance of F (6) and the orthogonality of ⊤P , namely β(⊤,tk) β(P x,t_k) =h‖g(⊤,tk)‖=h‖⊤g(,tk)‖ =h\|g(P x,t_k)\|=h\|P g(x,t_k)\| =h‖g(,tk)‖=β(,tk). =h\|g(x,t_k)\|=β(x,t_k). Thus, the Euler sampler is graph stable. Substituting α and βmax=hCg(0) _ =hC_g(x_0) into Theorem 2 concludes the proof. Heun sampler. Applying the exact same logic, we evaluate the state stability of the Heun update map. Using the triangle inequality and applying the Lipschitz constant M sequentially ‖Tk(;)−Tk(;)‖ \|T_k θ(x;S)-T_k θ(y;S)\| ≤‖−‖+h2‖uk(;)−uk(;)‖ ≤\|x-y\|+ h2\|u_k θ(x;S)-u_k θ(y;S)\| +h2‖uk+1(+huk(;);)−uk+1(+huk(;);)‖ + h2 \|u_k+1 θ (x+hu_k θ(x;S);S )-u_k+1 θ (y+hu_k θ(y;S);S ) \| ≤‖−‖+hM2‖−‖ ≤\|x-y\|+ hM2\|x-y\| +hM2(‖−‖+h‖uk(;)−uk(;)‖) + hM2 (\|x-y\|+h\|u_k θ(x;S)-u_k θ(y;S)\| ) ≤‖−‖+hM2‖−‖+hM2(‖−‖+hM‖−‖) ≤\|x-y\|+ hM2\|x-y\|+ hM2 (\|x-y\|+hM\|x-y\| ) =(1+hM+h2M22)‖−‖. = (1+hM+ h^2M^22 )\|x-y\|. Thus, the Heun sampler is state stable with α=1+hM+h2M22α=1+hM+ h^2M^22. For graph stability, let ~≔+huk(;~) z +hu_k θ(x; S) and ≔0⊤+huk(0⊤;)z _0 x+hu_k θ(P_0 x;S). The error in the update map is bounded by ∥0⊤Tk(;~) \|P_0 T_k θ(x; S) −Tk(0⊤;)∥ -T_k θ(P_0 x;S)\| ≤h2‖0⊤uk(;~)−uk(0⊤;)‖ ≤ h2\|P_0 u_k θ(x; S)-u_k θ(P_0 x;S)\| +h2‖0⊤uk+1(~;~)−uk+1(;)‖. + h2\|P_0 u_k+1 θ( z; S)-u_k+1 θ(z;S)\|. The first term is bounded directly by h2‖g(,tk)‖(Γε+(ε2)) h2\|g(x,t_k)\|( +O( ^2)). For the second term, we add and subtract 0⊤uk+1(0;~)P_0 u_k+1 θ(P_0z; S) and use the Lipschitz continuity of the vector field and its stability evaluated at z to obtain ∥0⊤ \|P_0 uk+1(~;~)−uk+1(;)∥ u_k+1 θ( z; S)-u_k+1 θ(z;S)\| ≤‖0⊤uk+1(~;~)−0⊤uk+1(0;~)‖ ≤\|P_0 u_k+1 θ( z; S)-P_0 u_k+1 θ(P_0z; S)\| +‖0⊤uk+1(0;~)−uk+1(;)‖ +\|P_0 u_k+1 θ(P_0z; S)-u_k+1 θ(z;S)\| ≤M‖~−0‖+‖g(,tk+1)‖(Γε+(ε2)). ≤ M\| z-P_0z\|+\|g(z,t_k+1)\| ( +O( ^2) ). Since ‖~−0‖=h‖0⊤uk(;~)−uk(0⊤;)‖≤h‖g(,tk)‖(Γε+(ε2))\| z-P_0z\|=h\|P_0 u_k θ(x; S)-u_k θ(P_0 x;S)\|≤ h\|g(x,t_k)\|( +O( ^2)), we can combine everything to bound the total graph error by h2(‖g(,tk)‖+hM‖g(,tk)‖+‖g(,tk+1)‖)(Γε+(ε2)). h2 (\|g(x,t_k)\|+hM\|g(x,t_k)\|+\|g(z,t_k+1)\| ) ( +O( ^2) ). Hence, the Heun sampler is graph stable with bounding function β(,tk)≔h2((1+hM)‖g(,tk)‖+‖g(,tk+1)‖).β(x,t_k) h2 ((1+hM)\|g(x,t_k)\|+\|g(z,t_k+1)\| ). (26) By the definition of Cg(0)C_g(x_0), evaluated along the unperturbed trajectories, we have ‖g(k,tk)‖≤Cg(0)\|g(x_k,t_k)\|≤ C_g(x_0) and ‖g(,tk+1)‖≤Cg(0)\|g(z,t_k+1)\|≤ C_g(x_0). Therefore, evaluating the bounding function along this trajectory yields βmax _ ≤h2((1+hM)Cg(0)+Cg(0)) ≤ h2 ((1+hM)C_g(x_0)+C_g(x_0) ) =h(1+hM2)Cg(0). =h (1+ hM2 )C_g(x_0). Making the required substitutions in the stability bound of Theorem 2, we find min∈‖⊤Ψk(0;~)−Ψk(⊤0;)‖ _P \|P _k θ(x_0; S)- _k θ(P x_0;S)\| ≤αk−1α−1βmax(Γε+(ε2)) ≤ α^k-1α-1 _ ( +O( ^2) ) ≤(1+hM+h2M22)k−1MCg(0)(Γε+(ε2)). ≤ (1+hM+ h^2M^22 )^k-1MC_g(x_0) ( +O( ^2) ). ∎ Appendix I Wasserstein Stability of the Discrete Flow Proof of Corollary 4. Let 0∈P_0 be the specific permutation matrix from Section I-B that minimizes the relative perturbation error. Evaluating the discrete flow on the exact same initial random vector 0x_0 under both the nominal graph S and the perturbed graph ~ S forms a valid coupling. Thus, we can upper bound the Wasserstein-2 distance using the bound derived in Theorem 2, to obtain min∈W2(π~k,πk)2 _P W_2( π_k^P, _k^P)^2 ≤0[‖0⊤Ψk(0;~)−Ψk(0⊤0;)‖2] _x_0 [\|P _0 _k θ(x_0; S)- _k θ(P _0x_0;S)\|^2 ] ≤0[(Υk(α)βmax(0)(Γε+(ε2)))2]. _x_0 [ ( _k(α) _ (x_0) ( +O( ^2) ) )^2 ]. Because Υk(α) _k(α) and (Γε+(ε2)) ( +O( ^2) ) are deterministic and independent of 0x_0, we can factor them out of the expectation 0 _x_0 [(Υk(α)βmax(0)(Γε+(ε2)))2] [ ( _k(α) _ (x_0) ( +O( ^2) ) )^2 ] =Υk(α)2(Γε+(ε2))20[βmax(0)2] = _k(α)^2 ( +O( ^2) )^2E_x_0 [ _ (x_0)^2 ] =Υk(α)2(Γε+(ε2))2(βmax′)2. = _k(α)^2 ( +O( ^2) )^2(β _ )^2. Taking the square root of both sides directly yields (21). ∎ Appendix J Lipschitz constant of the vector field Proof of Proposition 3. For a normalized Lipschitz pointwise nonlinearity σ(⋅)σ(·), i.e., |σ(b)−σ(a)|≤|b−a||σ(b)-σ(a)|≤|b-a|, we have ‖ℓ−ℓ‖F2≤‖∑p=0P−1p(ℓ−1−ℓ−1)ℓp‖F2.\|X_ -Y_ \|^2_F≤ \| _p=0^P-1S^p(X_ -1-Y_ -1) _ p \|_F^2. Let ℓ−1≔ℓ−1−ℓ−1Z_ -1 _ -1-Y_ -1 and define its GFT as ¯ℓ−1≔⊤ℓ−1 Z_ -1 Z_ -1. Substituting =⊤S=V yields ‖ℓ−ℓ‖F2≤‖∑p=0P−1p¯ℓ−1ℓp‖F2.\|X_ -Y_ \|^2_F≤ \| _p=0^P-1 ^p Z_ -1 _ p \|_F^2. Let ¯≔∑p=0P−1p¯ℓ−1ℓp Y _p=0^P-1 ^p Z_ -1 _ p and denote its i-th row as ¯i⊤ y_i . Because =diag(λ1,…,λN) = diag( _1,…, _N), the operations decouple across the rows in the spectral domain ¯i⊤=∑p=0P−1λip¯i⊤ℓp=¯i⊤ℓ(λi), y_i = _p=0^P-1 _i^p z_i _ p= z_i H_ ( _i), where ℓ(λi)≔∑p=0P−1λipℓpH_ ( _i) _p=0^P-1 _i^p _ p. Then ∥ℓ \|X_ −ℓ∥F2≤‖¯‖F2=∑i=1N‖¯i‖2=∑i=1N‖¯i⊤ℓ(λi)‖2 -Y_ \|^2_F≤ \| Y \|_F^2= _i=1^N\| y_i\|^2= _i=1^N\| z_i H_ ( _i)\|^2 ≤∑i=1N‖¯i‖2‖ℓ(λi)‖22≤Mℓ2∑i=1N‖¯i‖2 ≤ _i=1^N\| z_i\|^2\|H_ ( _i)\|^2_2≤ M_ ^2 _i=1^N\| z_i\|^2 =Mℓ2‖¯ℓ−1‖F2=Mℓ2‖ℓ−1−ℓ−1‖F2, =M_ ^2\| Z_ -1\|^2_F=M_ ^2\|X_ -1-Y_ -1\|_F^2, where MℓM_ is the layer-wise Lipschitz constant defined as Mℓ≔maxi∈1,…,N‖∑p=0P−1ℓpλip‖2. M_ _i∈\1,…,N\ \| _p=0^P-1 _ p _i^p \|_2. Applying this relation recursively across all L layers bounds the global spatial Lipschitz constant M≤∏ℓ=1LMℓM≤ _ =1^LM_ , which concludes the proof. ∎ References [1] K. Alomar, H. I. Aysel, and X. Cai (2023) Data augmentation in classification and segmentation: a survey and new strategies. J. Imaging 9 (2), p. 1–26. Cited by: §I. [2] R. Arghal, E. Lei, and S. S. Bidokhti (2022) Robust graph neural networks via probabilistic lipschitz constraints. Proc. Mach. Learn. Res. 168, p. 1–13. Cited by: §V. [3] O. Bousquet and A. Elisseeff (2002) Stability and generalization. J. Mach. Learn. Res. 2, p. 499–526. Cited by: §I. [4] M. M. Bronstein, J. Bruna, T. Cohen, and P. Veličković (2021) Geometric Deep Learning: Grids, Groups, Graphs, Geodesics, and Gauges. arXiv preprint arXiv:2104.13478. Cited by: §I. [5] J. C. Butcher (2016) Numerical methods for ordinary differential equations. John Wiley & Sons. Cited by: §IV. [6] J. Cerviño, L. Ruiz, and A. Ribeiro (2022) Training stable graph neural networks through constrained learning. In Proc. IEEE Intl. Conf. Acoustics, Speech and Signal Process. (ICASSP), p. 4223–4227. Cited by: §V. [7] J. Cerviño, L. Ruiz, and A. Ribeiro (2023) Learning by transference: training graph neural networks on growing graphs. IEEE Trans. Signal Process. 71 (), p. 233–247. Cited by: §VII. [8] H. Chen, N. Zilberstein, and S. Segarra (2026) Prior-informed flow matching for graph reconstruction. In ICML Wrksp. Struct. Prob. Inference and Gen. Modeling, p. 1–33. Cited by: §I. [9] R. T. Q. Chen, Y. Rubanova, J. Bettencourt, and D. Duvenaud (2018) Neural ordinary differential equations. In Proc. Adv. Neural. Inf. Process. Syst., p. 1–13. Cited by: §I. [10] M. Defferrard, X. Bresson, and P. Vandergheynst (2016) Convolutional neural networks on graphs with fast localized spectral filtering. In Proc. Adv. Neural. Inf. Process. Syst., p. 1–9. Cited by: §I-A. [11] X. Dong, D. Thanou, M. Rabbat, and P. Frossard (2019) Learning graphs from data: a signal representation perspective. IEEE Signal Process. Mag. 36 (3), p. 44–63. Cited by: §I. [12] S. Elfwing, E. Uchibe, and K. Doya (2018) Sigmoid-weighted linear units for neural network function approximation in reinforcement learning. Neural Networks 107, p. 3–11. Cited by: §VI-A. [13] F. Gama, J. Bruna, and A. Ribeiro (2020) Stability properties of graph neural networks. IEEE Trans. Signal Process. 68, p. 5680–5695. Cited by: §I, §I, §I-B, §I-B, §I-B, §V, §VI-A, §VI-B. [14] F. Gama, A. G. Marques, G. Leus, and A. Ribeiro (2019) Convolutional neural network architectures for signals supported on graphs. IEEE Trans. Signal Process. 67 (4), p. 1034–1049. Cited by: §I-A. [15] A. Gretton, K. M. Borgwardt, M. J. Rasch, B. Schölkopf, and A. Smola (2012) A kernel two-sample test. J. Mach. Learn. Res. 13 (25), p. 723–773. Cited by: §VI-A. [16] J. Ho, A. Jain, and P. Abbeel (2020) Denoising diffusion probabilistic models. In Proc. Adv. Neural. Inf. Process. Syst., p. 1–12. Cited by: §I. [17] P. W. Holland, K. B. Laskey, and S. Leinhardt (1983) Stochastic blockmodels: first steps. Social Networks 5 (2), p. 109–137. Cited by: §VI-A. [18] W. Huang, T. A. Bolton, J. D. Medaglia, D. S. Bassett, A. Ribeiro, and D. Van De Ville (2018) A graph signal processing perspective on functional brain imaging. Proc. IEEE 106 (5), p. 868–885. Cited by: §I. [19] W. Huang, L. Goldsberry, N. F. Wymbs, S. T. Grafton, D. S. Bassett, and A. Ribeiro (2016) Graph frequency analysis of brain signals. IEEE J. Sel. Topics Signal Process. 10 (7), p. 1189–1203. Cited by: §VI-A. [20] Human Connectome Project Note: https://w.humanconnectome.org/ Cited by: §VI. [21] E. Isufi, F. Gama, D. I. Shuman, and S. Segarra (2024) Graph filters for signal processing and machine learning on graphs. IEEE Trans. Signal Process. 72 (), p. 4745–4781. Cited by: §I-A. [22] I. Jabłoński (2017) Graph signal processing in applications to sensor networks, smart grids, and smart cities. IEEE Sensors J. 17 (23), p. 7659–7666. Cited by: §I. [23] J. Jo, S. Lee, and S. J. Hwang (2022) Score-based generative modeling of graphs via the system of stochastic differential equations. In Proc. Int. Conf. Mach. Learn., p. 10362–10383. Cited by: §I. [24] H. Kenlay, D. Thanou, and X. Dong (2021) Interpretable stability bounds for spectral graph filters. In Proc. Int. Conf. Mach. Learn., p. 5388–5397. Cited by: §I. [25] T. N. Kipf and M. Welling (2017) Semi-supervised classification with graph convolutional networks. In Proc. Int. Conf. Learn. Representations, p. 1–14. Cited by: §I-A. [26] J. Köhler, L. Klein, and F. Noé (2020) Equivariant flows: exact likelihood generative learning for symmetric densities. In Proc. Int. Conf. Mach. Learn., p. 5361–5370. Cited by: §I. [27] E. D. Kolaczyk (2009) Statistical analysis of network data: methods and models. Springer, New York, NY. Cited by: §I. [28] R. Levie, W. Huang, L. Bucci, M. Bronstein, and G. Kutyniok (2021) Transferability of spectral graph convolutional neural networks. J. Mach. Learn. Res. 22, p. 1–59. Cited by: §I. [29] Y. Lipman, R. T. Chen, H. Ben-Hamu, M. Nickel, and M. Le (2023) Flow matching for generative modeling. In Proc. Int. Conf. Learn. Representations, p. 1–28. Cited by: §I, §I-C, §I-C, §VI-A. [30] C. Liu, W. Fan, Y. Liu, J. Li, H. Li, H. Liu, J. Tang, and Q. Li (2023) Generative diffusion models on graphs: methods and applications. In Proc. Int. Joint Cont. Artificial Intell., p. 6702–6711. Cited by: §I. [31] J. Liu, A. Kumar, J. Ba, J. Kiros, and K. Swersky (2019) Graph normalizing flows. In Proc. Adv. Neural. Inf. Process. Syst., p. 1–11. Cited by: §I. [32] G. Mateos, S. Segarra, A. G. Marques, and A. Ribeiro (2019-05) Connecting the dots: identifying network structure via graph signal processing. IEEE Signal Process. Mag. 36 (3), p. 16–43. Cited by: §I. [33] A. Q. Nichol and P. Dhariwal (2021) Improved denoising diffusion probabilistic models. In Proc. Int. Conf. Mach. Learn., p. 8162–8171. Cited by: §I. [34] C. Niu, Y. Song, J. Song, S. Zhao, A. Grover, and S. Ermon (2020) Permutation invariant graph generation via score-based generative modeling. In Proc. Int. Conf. Artif. Intell. Statist., p. 1–10. Cited by: §I. [35] A. Ortega, P. Frossard, J. Kovačević, J. M. Moura, and P. Vandergheynst (2018) Graph signal processing: overview, challenges, and applications. Proc. IEEE 106 (5), p. 808–828. Cited by: §I, §I-A. [36] M. G. Preti and D. Van De Ville (2019) Decoupling of brain function from structure reveals regional behavioral specialization in humans. Nat. Commun. 10, p. 1–7. Cited by: §VI-A. [37] M. Raissi, P. Perdikaris, and G. E. Karniadakis (2019) Physics-informed neural networks: a deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations. J. Comput. Phys. 378, p. 686–707. Cited by: §I. [38] R. Ramakrishna and A. Scaglione (2021) Grid-graph signal processing (Grid-GSP): a graph signal processing framework for the power grid. IEEE Trans. Signal Process. 69, p. 2725–2739. Cited by: §I. [39] S. Rozada, Vimal K B, A. Cavallo, A. G. Marques, H. Jamali-Rad, and E. Isufi (2026) Graph-aware diffusion for signal generation. In Proc. IEEE Intl. Conf. Acoustics, Speech and Signal Process. (ICASSP), p. 461–465. Cited by: §I, §I, §VI-A. [40] L. Ruiz, L. F. O. Chamon, and A. Ribeiro (2023) Transferability properties of graph neural networks. IEEE Trans. Signal Process. 71, p. 3474 – 3489. Cited by: §I. [41] L. Ruiz, F. Gama, and A. Ribeiro (2021) Graph neural networks: architectures, stability and transferability. Proc. IEEE 109 (5), p. 660–682. Cited by: §I-B. [42] K. Scaman and A. Virmaux (2018) Lipschitz regularity of deep neural networks: analysis and efficient estimation. In Proc. Adv. Neural. Inf. Process. Syst., p. 1–10. Cited by: §VI-B. [43] Y. Song, J. Sohl-Dickstein, D. P. Kingma, A. Kumar, S. Ermon, and B. Poole (2021) Score-based generative modeling through stochastic differential equations. In Proc. Int. Conf. Learn. Representations, p. 1–36. Cited by: §I. [44] C. Szegedy, W. Zaremba, I. Sutskever, J. Bruna, D. Erhan, I. Goodfellow, and R. Fergus (2014) Intriguing properties of neural networks. In Proc. Int. Conf. Learn. Representations, p. 1–10. Cited by: §I. [45] Y. B. Uslu, S. Hadou, S. S. Bidokhti, and A. Ribeiro (2026; see also arXiv:2604.05175 [eess.SP]) Graph signal diffusion models for wireless resource allocation. In spawc, p. . Cited by: §I. [46] Y. B. Uslu, S. Hadou, S. Rozada, S. S. Bidokhti, and A. Ribeiro (2026) Graph signal generative diffusion models. In Proc. IEEE Intl. Conf. Acoustics, Speech and Signal Process. (ICASSP), p. 626–630. Cited by: §I, §I. [47] C. Vignac, I. Krawczuk, A. Siraudin, B. Wang, V. Cevher, and P. Frossard (2023) DiGress: discrete denoising diffusion for graph generation. In Proc. Int. Conf. Learn. Representations, p. 1–22. Cited by: §I. [48] C. Villani et al. (2009) Optimal transport: old and new. Vol. 338, Springer. Cited by: §I-B. [49] M. Yang (2025) Topological schrödinger bridge matching. In Proc. Int. Conf. Learn. Representations, p. 1–42. Cited by: §I. [50] B. Yu, H. Yin, and Z. Zhu (2018) Spatio-temporal graph convolutional networks: a deep learning framework for traffic forecasting. In Proc. Int. Joint Cont. Artificial Intell., p. 3634–3640. Cited by: §I. [51] A. Zalesky, A. Fornito, and E. Bullmore (2012) On the use of correlation as a measure of network connectivity. Neuroimage 60 (4), p. 2096–2106. Cited by: §I. [52] Y. Zhu, C. Wang, Q. Zhang, and H. Xiong (2024) Graph signal diffusion model for collaborative filtering. In Proc. ACM SIGIR Conf. Res. Develop. Inf. Retrieval, p. 1380–1390. Cited by: §I. [53] D. Zügner, A. Akbarnejad, and S. Günnemann (2018) Adversarial attacks on neural networks for graph data. Intl. Conf. Knowledge Discovery Data Mining, p. 2847–2856. Cited by: §I.