1 Introduction↩︎

Graph prediction is an important topic in statistics, machine learning, and signal processing, and is motivated by various applications, e.g., protein-protein interaction networks [1], wind power prediction [2], and user behavior modeling in social networks [3]. There can be various versions of graph prediction problems. We are interested in the scenario where one observes nodal (multi-dimensional) features \(X\), nodal responses \(Y\), and graph topology information. The prediction algorithm typically will then leverage graph topology information to facilitate the learning of predicting \(Y\) given \(X\). In practice, Graph Neural Network (GNN) models [4] have become a popular choice due to their expressiveness power and scalable computation. In the so-called inverse of a graph prediction problem, one would like to infer the input graph nodal features \(X\) given an outcome response \(Y\). Such a problem is of interest in various real-world applications, for instance, in molecular design [5], scientists want to infer features of the molecule that lead to certain outcomes; in power outage analysis [6], [7], we are interested in finding features (e.g., weather or past power output by sensors) that leads to outages for future prevention. The inverse graph prediction problem is the focus of the current work, for which we develop an invertible deep model that can be efficiently applied to graph data.

The prediction and inverse prediction problems are illustrated in Fig. 1. The left plot shows the case of general non-graph data, and the right plot shows the case with graph data \(X\) and graph label \(Y\). As an illustrative example, suppose the graph is a power grid, where each node \(v \in \{1,\ldots,5\}\) denotes a power generator, \(X_v \in \mathbb{R}^{d'}\) denotes historical power outputs from generator \(v\), and \(Y_v\) denotes the status of the generator \(v\) (e.g., \(Y_v=0\) for functioning normally, and 1 for anomaly). For system monitoring, it is useful to predict the status of generators given historical observations (the forward prediction, indicated by solid black arrows), as well as to generate unobserved possible circumstances given generator status (the inverse, indicated by dashed black arrows) for cause analyses.

To formalize the notion of inverse prediction, we adopt some probability terminologies. The forward prediction problem is to predict \(Y\) from input \(X\), which can be formulated as learning the conditional probability of \(p(Y|X)\). The inverse prediction problem is to learn the conditional probability of \(p(X|Y)\) and to generate samples \(X\) from it. Note that a discriminative task seeks to estimate the posterior \(p(X|Y)\) at a given (test) point \(X\), which can be done by a conventional classification model. The generative task is different in that one asks to sample \(X\) according to \(p(X|Y)\), that is, given a label \(Y\), to produce samples \(X\) which do not exist in the training data nor any provided test set. The problem is challenging when data \(X\) is in high dimensional space, where a grid of \(X\) can not be efficiently constructed. In particular, when \(X\) is graph nodal feature data, the dimensionality of \(X\) scales linearly with the graph size \(N\) (the number of nodes in the graph). In the case of categorical response, the graph label \(Y\) assigns one of the \(K\)-class labels to each node, which makes the total possible outcomes \(K^N\) many. Thus the inverse prediction problem on graph data poses both modeling and computational challenges when scalability to large graphs is needed.

In this work, we develop a deep generative model for the conditional generation task of the inverse prediction problem on graphs. While deep generative models like generative adversarial networks (GAN) [8], [9], variational auto-encoders (VAE) [10], [11] and normalizing flow networks [12] have been intensively developed in recent years, the conditional generation task given categorial prediction labels is less studied. We further review the literature and comment on related works in Section 2.1. Unlike previous conditional generative models [13]–[17], which typically concatenate the prediction label \(Y\) with the random code \(Z\) as input or rely on curated forms of neural network (NN) layers to ensure invertibility, we propose to encode the input data \(X\) one-to-one by an invertible network to an intermediate feature \(H\), from which the label \(Y\) can be predicted using a linear classifier, and in the other direction \(H\) can also be generated from \(Y\) by a parametric mixture model. The framework of our model is shown in Fig. 2. Because the general (non-graph) data case can be seen as a special case of graph data where the graph only has one node, we call the proposed model invertible Graph Neural Network (iGNN), as a unified name.

To efficiently handle the up to exponentially many outcomes of graph nodal labels, we introduce a factorized mixture distribution \(H | Y\) over graphs. We theoretically show that the capacity to learn the conditional distribution \(X|Y\) over a graph can be achieved when the per-node Gaussian mixture model \(H_v | Y_v\) has well-separated components, and the need separation increases with graph size \(N\) only mildly (an additive \({\log N}\) factor), cf. Proposition 1. Under this formulation, the \(Y\)-to-\(H\) sub-network is light in model and computation. The major part of the model capacity of the iGNN model is in the \(H\)-to-\(X\) sub-network, which is an invertible flow model based on Residual Network (ResNet) framework [18]. The computational scalability to large graphs is tackled by adopting GNN layers, which are not compatible with normalizing flow models having constrained formats in architecture [19], [20], and would also complicate the spectral normalization technique in [18] to ensure network invertibility.

To overcome these issues, we propose Wasserstein-2 regularization of the ResNet, which in the limit of a large number of residual blocks recovers the transport cost in the dynamic formulation of optimal transport (OT) [21]–[23]. Based on known results of OT theory, the Wasserstein-2 regularization leads to smooth trajectories of the transported densities. The invertibility of the induced flow map can then be guaranteed with a sufficient number of residual blocks. We empirically verify the invertibility of the flow model in experiments. This allows our model to use free-form neural network layers inside the ResNet blocks, including any GNN layer types. Theoretically, we study the existence and invertibility of flow mapping based on theories of OT and diffusion process. We also prove the expressiveness of GNN layers to approximate the velocity field in the theoretical flows when data is a Gaussian field on the graph.

We examine the performance of the proposed iGNN model on simulated graph data as well as data from real-world graph prediction applications. In summary, the contributions of the work are

\(\bullet\) We propose a two-step procedure, \(Y\)-to-\(H\) and \(H\)-to-\(X\), to tackle the generative task of inverse prediction problem viewed as a conditional generation problem and develop an invertible flow model consisting of two subnetworks accordingly. The model is made scalable to graph data by a factorized formulation of the parametric mixture model \(H|Y\) and the GNN layers in the invertible flow network between \(H\) and \(X\).

\(\bullet\) We introduce Wasserstein-2 regularization of the invertible flow network, which is computationally efficient and compatible with general free-form layer types, particularly the GNN layers in the flow model. The effect on preserving invertibility is backed by OT theory and verified in practice.

\(\bullet\) The existence of invertible flow maps is analyzed theoretically. Based on the theoretical flows, we analyze the expressiveness of spectral and spatial graph convolution layers applied to graph data.

\(\bullet\) The proposed iGNN model is applied to both simulated and real-data examples, showing improved generative performance over alternative conditional generation models.

Notation. We denote by \([n] = \{1,\cdots, n \}\) the integer set. \(\mathbb{E}_{x\sim p}\) denotes the expectation over \(x \sim p\), that is, \(\mathbb{E}_{x\sim p} f(x) = \int f(x) p(x) dx\). When a probability distribution has density, we use the same notation, e.g., \(p\), to denote the distribution and the density. For real-diagonalizable matrix \(A\), and function \(f: \mathbb{R}\to \mathbb{R}\), we denote by \(f(A)\) the function of the matrix \(A\). For \(T: \mathbb{R}^d \to \mathbb{R}^d\) and \(p\) a distribution on \(\mathbb{R}^d\), we denote by \(T_\# p\) the push-forward of distribution \(p\) by \(T\), i.e., \((T_\# p)(\cdot) = p( T^{-1} (\cdot))\).

2 Background↩︎

2.1 Related works↩︎

Generative deep models and normalizing flow. At present, generative adversarial networks (GAN) [8], [9] and variational auto-encoders (VAE) [10], [11] are two of the most popular frameworks that have achieved various successes [24]–[26]. However, they also suffer from clear limitations such as notable difficulties in training, such as mode collapse [27] and posterior collapse [28]. On the other hand, normalizing flows (see [12] for a comprehensive review) estimate arbitrarily complex densities via the maximum likelihood estimation (MLE), and they transport original random features \(X\) into distribution that are easier to sample from (e.g., standard multivariate Gaussian) through invertible neural networks. Flow-based models can be classified into two broad classes: the discrete-time models (some of which include coupling layers [19], autoregressive layers [20] and residual networks [18], [29]), and the continuous-time models as exemplified by neural ODE [30], [31]. Most normalizing flow methods focus on unconditional generation with little development in a conditional generation. In addition, to achieve numerically reliable training, regularization of the density transport trajectories in flow networks are necessary but remain a challenge. In this work, we propose Wasserstein-2 regularization motivated by the transport cost in the dynamic formula of optimal transport (OT) theory. In Remark 2, we compare with the original spectral normalization in iResNet [18] in more detail.

Conditional generation networks. The conditional generation versions of GAN (cGAN) have been studied in several places [13], [14], where the prediction outcome \(Y\) (one-hot encoded as \(e_Y\)) are typically concatenated with random noise \(Z\) and taken as input to the generator network, as illustrated in the left of Fig. 2. The one-hot encoding of categorical \(Y\) concatenated with Gaussian \(Z\) poses challenges in training cGAN models, in addition to known issues of their unconditional counterparts, such as mode collapse, posterior collapse, and failures to provide exact data likelihood. When the prediction label \(Y\) is lying on a graph having \(N\) nodes, the one-hot coding will increase up to \(O(N)\) more coordinates to the input \((Z,e_Y)\), which significantly increases the model complexity and computational load with large graphs. Conditional invertible neural network (cINN) model was developed in [15] for analyzing inverse problems. The model inherits the approach of cGAN models to concatenate one-hot encoded prediction label \(e_Y\) with normal code \(Z\) while using Real-NVP layers [19] to ensure neural network invertibility. In terms of training objective, [15] proposed to use maximum mean discrepancy (MMD) losses to encourage both the matching of the input data distribution and the independence between label \(Y\) and normal code \(Z\). We call the method in [15] cINN-MMD. Replacing the MMD losses with a flow-based objective, [16], [17] extended the invertible network approach in [15] and applied to image generation problems. In the models in [16] and [17], which we call cINN-Flowand cINN-Flow+respectively, the inputs to the Real-NVP layers contain encoded information of the prediction label \(Y\) so as to learn label-conditioned generation, see more in Fig. 10. In experiments, we compare with cGAN and cINN models on both non-graph and graph data.

2.2 Transport cost regularization↩︎

We review needed background on the transport cost and the dynamic formula of Wasserstein-2 optimal transport (OT). The transport cost has been proposed to regularize the normalizing flow neural network models, e.g., in [31]. Consider the trajectory \(x(t)\) satisfying \(\dot{x}(t) = v( x(t), t)\), that is, \(v(x,t)\) is the velocity field, where \(x(0) \sim p\) and \(t \in [0,1]\). The transport cost is defined as \[\label{eq:def-transport-cost} \mathcal{T} := \int_0^1 \mathbb{E}_{x \sim \rho(\cdot,t) } \| {v}(x,t) \|^2 dt,\tag{1}\] where \(\rho(x,t)\) is the marginal distribution of \(x(t)\) after pushing forward the initial distribution \(\rho(\cdot, 0) = p\) by \(v(x,t)\) for time \(t\). The cost 1 can be viewed as taking the kinetic energy in computing the Lagrangian action along the trajectory. A more general form of cost under the framework of Mean-field Games has been proposed in the current and independent work [32] to regularize trajectories of normalizing flows. It is also known that minimizing the cost \(\mathcal{T}\) under the constraint of transporting from fixed \(p = \rho(\cdot, 0)\) to \(q = \rho(\cdot, 1)\) leads to the dynamic formulation of OT, i.e., the Benamou-Brenier formula [21], [22]. Specifically, the solution of the minimization \[\label{eq:Benamou-Brenier} \begin{align} &\inf_{\rho, v} \int_0^1 \mathbb{E}_{x \sim \rho( \cdot, t)} \| {v}(x,t) \|^2 dt \\ & s.t. \quad \partial_t \rho + \nabla \cdot (\rho v) = 0, \quad \rho( \cdot, 0) = p, \quad \rho( \cdot, 1) = q, \end{align}\tag{2}\] recovers the Wasserstein-2 OT from \(p\) to \(q\), that is, at the minimizer \(v(x,t)\), \(\rho(\cdot, 1) = q\) and \(\mathcal{T} = W_2^2(p,q)\). The minimizer \(v(x,t)\) of 2 can be interpreted as the optimal control of the transport problem from \(p\) to \(q\).

In terms of neural network implementation, OT-Flow [31] used the potential model based on OT theory and parametrized the potential \(\Phi( x,t)\) by a neural network, where \(\nabla \Phi( x,t)\) gives the velocity \(v(x,t)\). Our method adopts a ResNet base model, where the invertibility is fulfilled by a Wasserstein-2 regularization, which can be viewed as a finite-step discrete-time counterpart of the transport cost \(\mathcal{T}\). The proposed Wasserstein-2 regularization recovers \(\mathcal{T}\) with a large number of steps, yet does not require integrating the continuous-time flow with accuracy for all time, cf. Remark 3.

3 Methods↩︎

Given data-label pairs \(\{X_i, Y_i\}\), we first describe our approach when \(X_i\) is a sample in \(\mathbb{R}^d\), and \(Y_i\) is the categorical label in \(K\) classes, i.e., \(Y_i \in [K]\). The case where both \(X_i\) and \(Y_i\) lie on a graph is addressed in Section 3.3 where we make the approach scalable to large graphs. In Section 3.4, we introduce Wasserstein-2 regularization to preserve the invertibility of the trained flow network with free-form residual block layer types. All proofs are in Appendix 8.

3.1 Inverse of prediction as conditional generation↩︎

The overall framework is to (end-to-end) train a network consisting of two sub-networks: the first sub-network maps invertibly from \(X\) to an intermediate representation \(H\), and the second sub-network maps from \(H\) to label \(Y\), which is a classifier and loses information. Specifically,

\(\bullet\) \(H\)-\(Y\) classification sub-network. We model \(H|Y\) by a Gaussian mixture model with “well-separated” means (detailed in Section 3.2). The parametric form of \(H|Y\) contains trainable parameter \(\theta_c^{g}\), and the generation of \(H|Y\) is by sampling the corresponding mixture component accordingly. The prediction of label \(Y\) from \(H\) can be conducted by a linear classifier parametrized by \(\theta_c^{c}\). The trainable parameters in the classification sub-network are denoted as \(\theta_c=(\theta_c^{g},\theta_c^{c})\).

\(\bullet\) \(X\)-\(H\) invertible sub-network. The invertible mapping from \(X\) to \(H\) is by a flow ResNet in \(\mathbb{R}^d\) (detailed in Section 3.4). The sub-network parameters are denoted as \(\theta = \{\theta_l, l=1,\cdots, L\}\), where \(L\) is the number of residual blocks, and the network mapping is denoted as \(F_\theta\). The prediction of label \(Y\) from input \(X\) is by first computing \(H = F_\theta(X)\) and then applying the \(H\)-to-\(Y\) classifier sub-network parametrized by \(\theta_c^c\). The generation of \(X|Y\) is by inversely mapping \(X = F_\theta^{-1}(H)\) once \(H\) is sampled according to \(p(H|Y)\) parametrized by \(\theta_c^g\).

The end-to-end training objective of the proposed network can be written as \[\label{end-to-end-obj} \underset{\{\theta, \theta_c\}}{\min} \;\mathcal{L}_g + \mu \mathcal{L}_c + \gamma \mathcal{W} ,\tag{3}\] where \(\mathcal{L}_g\), \(\mathcal{L}_c\) and \(\mathcal{W}\) are the generative loss, the classification loss and the Wasserstein-2 regularization, respectively. The scalars \(\mu, \gamma \ge 0\) are penalty factors. We will explain the choice of \(\mu\) later in this subsection, and the choice of \(\gamma\) is explained after the derivation in Section 3.4.

Given \(n_{\rm tr}\) many data-label pairs \(\{X_i, Y_i\}\), the generative loss is defined as \(\mathcal{L}_g = \frac{1}{n_{\rm tr}} \sum_{i=1}^{n_{\rm tr}} \ell_g ( X_i, Y_i)\) where \[\begin{align}- \ell_g ( X_i, Y_i) & =\log p_{H|Y_i}(F_{\theta}(X_i)) +\log |\det J_{F_{\theta}}(X_i)|.\label{Nflow95adaptive} \end{align}\tag{4}\] Because \(p_{H|Y}\) is a mixture model parametrized by \(\theta_c^g\), the per-sample loss \(\ell_g ( X_i, Y_i)\) is determined by both \(\theta\) and \(\theta_c^g\). We specify the construction and training of \(p_{H|Y}\) in Section 3.2. The mapping \(F_\theta\) is by an invertible ResNet. The construction of \(F_\theta\), including the computation of \(\log |\det J_{F_{\theta}}|\) and the regularization loss \(\mathcal{W}\), will be explained in Section 3.4.

Finally, the classification loss is defined as \(\mathcal{L}_c = \frac{1}{n_{\rm tr}} \sum_{i=1}^{n_{\rm tr}} \ell_c ( F_\theta(X_i), Y_i)\), where \(\ell_c( H_i, Y_i)\) is the per-sample \(K\)-class classification cross-entropy loss computed by softmax. We use a linear classifier to predict \(Y\) given \(H\), parametrized by \(\theta_c^c\). When \(F_\theta(X) | Y\) is close in distribution to the Gaussian mixture model \(H|Y\) specified by \(\theta_c^g\), one may also use \(\theta_c^g\) to construct the linear classifier from \(H\) to \(Y\), that is, to tie the parameters \(\theta_c^c\) and \(\theta_c^g\). Here we separately parametrize the forward and inverse prediction parameters in \(H\)-\(Y\) sub-network so as to facilitate optimization since, in both directions, the model is light. For the penalty factor \(\mu\), we find the experimental results insensitive to the choice. We use \(\mu = 1\) in all experiments.

3.2 Mixture model of \(H|Y\) and the shared flow↩︎

We parametrize \(p(H|Y)\) by a Gaussian mixture model as \[\label{eq:H124Y-one-node} H | Y=k \sim \mathcal{N}( \mu_k, \sigma^2 I_d), \quad k=1,\cdots, K,\tag{5}\] where \(\sigma > 0\) is a prefixed parameter. The \(K\) mean vectors \(\mu_k\) will be trainable, and we use isotropic Gaussian with the same covariance matrix for simplicity, which may be generalized. Note that we only parametrize the locations of means of the Gaussian components, and there is no need to specify the weights (fraction of samples in each component) because the fractions will be determined by the data. Fig. 3 provides an example where the fractions in the three classes differ.

We will initialize and train the mean vectors \(\mu_k\) to be sufficiently separated to ensure the Gaussian components have almost non-overlapping supports. Formally, we say that a collection of objects \(S\) in \(\mathbb{R}^d\) are \(\delta\)-separated if for any two distinct elements \(s, s' \in S\), \(d(s, s') \ge \delta\). When \(S\) is a collection of points in \(\mathbb{R}^d\), \(x, x' \in S\), then \(d(x,x') = \| x-x'\|\). When \(S\) is a collection of sets in \(\mathbb{R}^d\), and \(A, A' \in S\), we define \(d(A,A') = \inf_{x \in A, \, x' \in A'} \| x-x'\|\). A set \(S\subset \mathbb{R}^d\) is called an \(\epsilon\)-support of a probability distribution \(p\) if \(p(S^c) \le \epsilon\). When the means of Gaussian components are sufficiently separated, one would expect that there are \(\epsilon\)-supports of the \(K\) components that are mutually separated. Technically, we have that \(3\rho\)-separation of the means guarantees \(\rho\)-mutual separation of \(K\) \(\epsilon\)-supports, where \(\rho \sim \sigma \sqrt{\log (1/\epsilon)}\). The formal statement is proved in the following lemma.

The separation condition of the mixture model, especially the high dimensional counterpart in Proposition 1, may potentially be connected to the coding theory of Gaussian channel [33]. We further comment on this in the discussion section. The purpose of having separated \(\epsilon\)-supports of the \(K\) components in \(H|Y\) is for the construction of an invertible flow map from \(H\) to \(X\) that transports each class of samples correspondingly. Intuitively, the separated \(\epsilon\)-supports \(\Omega_k^{(H)}\) of the \(K\) components of \(p(H|Y)\) allow to partition \(\mathbb{R}^d\) into \(K\) domains with smooth boundary where each domain contains \(\Omega_k^{(H)}\). If the \(K\)-class distributions \(p(X|Y)\) also have separated \(\epsilon\)-supports \(\Omega_k^{(X)}\), then one can try to construct an invertible flow mapping in \(\mathbb{R}^d\) which transports \(p_{X|Y=k}\) to \(p_{H|Y=k}\) by transporting from \(\Omega_k^{(X)}\) to \(\Omega_k^{(H)}\) respectively. While the transport from each \(p_{X|Y=k}\) to \(p_{H|Y=k}\) can be constructed individually, the simultaneous transports of the \(K\) components by a shared flow would need the transported components to stay separated from each other along the flow, and this may not always be possible, e.g., due to topological constraints. For the shared invertible flow to exist, we introduce the following assumption:

While theoretical justification of Assumption 1 goes beyond the scope of the paper and is postponed here, we demonstrate the validity of the assumption empirically. An example of the shared invertible flow on a 3-class data in \(\mathbb{R}^2\) is illustrated in Fig. 3. In all our experiments, we find that the trained invertible ResNet successfully transports from \(p(X|Y)\) to \(p(H|Y)\) for all \(K\) classes. In the special case of having \(K=1\) class, the conditional generation problem reduces to an unconditional one and the shared flow is no longer an issue. In this case, one can set the distribution of \(H\) to be a standard normal, and the desired transport map is a normalizing flow. In Section 4 we further study the existence and theoretical properties of the flow mapping when \(K=1\), for general data and graph data.

Assumption 1 guarantees the existence of a shared flow mapping which transports the \(K\)-class conditional densities up to \(O(\epsilon)\) error as long as the components of the Gaussian mixture distribution \(p(H|Y)\) have \(\delta\)-separated \(\epsilon\)-supports. Meanwhile, by Lemma 1, \(p(H|Y)\) will have \(\rho\)-separated \(\epsilon\)-supports if the means \(\mu_k\) are separated as therein. Since \(\rho \sim \sigma \sqrt{\log (1/\epsilon)}\), it will exceed the \(O(1)\) constant \(\delta\) required by Assumption 1 when \(\epsilon\) is small. This suggests keeping the mean vectors \(\mu_k\) to be separated at the order of \(\sigma\sqrt{\log (1/\epsilon)}\). In practice, we initialize \(\mu_k\) to be separated at such a scale and preserve the separation during training via a barrier penalty on the distances \(\|\mu_k - \mu_{k'}\|\).

Note that Assumption 1 ensures the existence of a continuous-time flow, and we will use a ResNet model to represent the flow by finite-step composed transports which are regularized by summed discrete-time transport cost, to be detailed in Section 3.4.

3.3 Scalable conditional generation on graph data↩︎

The mixture modeling of \(H|Y\) described in the previous section is scalable to large graphs, as we explain in this subsection. Consider the graph data where both the data vector \(X\) and the label \(Y\) are defined on each node in a graph. In this subsection, the subscript \(_v\) indicates graph node \(v\). Suppose the graph \((V,E)\) has \(N\) nodes in \(V\), and the edge set \(E\) is specified by an adjacency matrix \(A\). We denote a graph data sample \(X\) and a graph label \(Y\) as \(X = [ X_1, \cdots, X_N]^T \in \mathbb{R}^{N\times d'}\), \(X_v \in \mathbb{R}^{d'}\), \(Y = [ Y_1, \cdots, Y_N]^T \in \mathbb{R}^N\), \(Y_v \in [K]\), \(v \in [N]\), where \(d'\) is the dimension of node feature. One may view \(X\) as a vector in \(\mathbb{R}^{d}\), \(d= d'N\), and label vector \(Y\) taking \(K^N\) many possibilities, and then apply the approach for non-graph data to the \(K^N\)-class prediction problem in \(\mathbb{R}^d\). However, this may lead to exponentially many classes when the graph size \(N\) is large, leading to difficulty in modeling \(H|Y\) by a mixture distribution as well as large model complexity.

To make our approach scalable to large graphs, we adopt two techniques: (i) we introduce a factorized form of \(H | Y\) which is independent and homogeneous over graph; (ii) we propose to use GNN layers in the invertible flow ResNet, which makes the neural network computation scalable. The scalability of our iGNN approach is demonstrated on a larger graph with \(N=500\) nodes, cf. Fig. 5.

(i) Factorized \(H|Y\) over graph. Recall that \(H\) is of same dimension as \(X\) and also defined on \(V\), \[H = [ H_1, \cdots, H_N]^T \in \mathbb{R}^{N\times d'}, \;H_v \in \mathbb{R}^{d'}.\] Suppose we have a \(K\)-component Gaussian mixture distribution of means \(\mu_k\) in \(\mathbb{R}^{d'}\), we specify the graph \(H|Y\) as \[\label{eq:H124Y-factored} p ( H | Y) = \prod_{v=1}^N p(H_v | Y_v), \quad H_v | Y_v \sim \mathcal{N}( \mu_{Y_v}, \sigma^2 I_{d'}),\tag{6}\] that is, the joint distribution of \(H|Y\) consists of independent and identical \(K\)-component Gaussian mixture distribution of \(H_v|Y_v\) in \(\mathbb{R}^{d'}\) across the graph. As a result, on the graph sample-label pair \(\{ X,Y \}\), the \(\log p(H|Y)\) term in the generative loss 4 can be computed as \[\label{Nflow95graph} \log p_{H|Y}( F_\theta(X)) = \sum_{v=1}^N \log p_{H_1| Y_1} ( (F_\theta(X))_v),\tag{7}\] where \(p_{H_1| Y_1}\) is specified by a Gaussian mixture model in \(\mathbb{R}^{d'}\).

The factorized form of \(H|Y\) reduces the complexity of modeling \(H|Y\) in \(\mathbb{R}^{d'N}\) to that of modeling a \(K\)-class mixture model in \(\mathbb{R}^{d'}\). At the same time, one would prefer that such a reduction does not prevent the mixture distribution \(H|Y\) from presenting any sufficiently separated conditional distribution \(X|Y\) via an invertible mapping in \(\mathbb{R}^{d'N}\). Under Assumption 1, the desired shared flow in \(\mathbb{R}^{d'N}\) exists as long as the \(O(1)\)-separation between some \(\epsilon\)-supports of the \(K^N\) components of \(p(H|Y)\) can be achieved. Due to that \(p(H|Y)\) is \(\mathcal{N}( \mu_Y, \sigma^2 I_{d' N})\) (\(\mu_Y\) is concatenated from \(\mu_{Y_v}\), \(v\in [N]\)), and the length of \(H-\mu_Y\) is of order \(\sigma \sqrt{d'N}\), this raises the question of how the separation of the Gaussian means of \(H_1|Y_1\) in \(\mathbb{R}^{d'}\) needs to scale with large \(N\). The following proposition shows that, theoretically, compared to Lemma 1 only an additional \(O(\sigma \sqrt{\log N })\) separation is needed.

The proposition suggests that the needed separation of Gaussian mixture means in \(\mathbb{R}^{d'}\) increases mildly when the graph size \(N\) increases. In experiments, we find that using the same algorithmic setting to preserve mean vectors’ separation in \(\mathbb{R}^{d'}\) suffices for graph data experiments where \(N\) is up to a few hundreds.

(ii) Scalable computation on the graph. As will be introduced in Section 3.4, the residual block in the proposed iGNN model can take a general form (the invertibility will be enforced by a Wasserstein transport regularization). This allows using any GNN layer type in the residual block, which can be made computationally scalable to large graphs [4].

Consider input graph data \(X\) in \(\mathbb{R}^{ N \times C}\), where \(N\) is the number of nodes, and \(C\) is the dimension of node features. We denote by \(Y\) the output graph data in \(\mathbb{R}^{ N \times C'}\) after the graph convolution, and we adopt spectral and spatial graph convolution layers in this work. Because \(X\) and \(Y\) can be hidden-layer features inside the residual block, the dimension \(C\) and \(C'\) may be larger than the data node feature dimension \(d'\). Let \(\tilde{L} \in \mathbb{R}^{N \times N}\) be the graph Laplacian matrix, possibly normalized so that the spectrum lines on a unit interval. Omitting the bias term and the non-linear node-wise activation function, the GNN layers used in this work are

\(\bullet\) ChebNet [34], where \(Y = \sum_{k=1}^K T_k(\tilde{L}) X \Theta_k\), parametrized by \(\Theta_k \in \mathbb{R}^{C \times C'}\), and \(T_k\) is the Chebyshev polynomial of degree \(k\).

\(\bullet\) L3Net [35], where \(Y = \sum_{k=1}^K B_kXA_k\), parametrized by \(B_k\) which are local filters on graph and \(A_k \in \mathbb{R}^{C \times C'}\).

The benefit of using GNN layer lies in the reduced model and computational complexities: the number of trainable parameters in a ChebNet layer is \(O(KCC')\), and that in an L3Net layer is \(O(K(CC'+Nv))\) where \(v\) denotes average local patch size [35]. In contrast, a fully-connected layer mapping from \(X\) to \(Y\) would have \(O(N^2CC')\) many parameters, which is not scalable to large graphs. We include the fully-connected layer baseline in experiments when the graph size is small and compare performance with the GNN layers.

3.4 Invertible flow network with Wasserstein-2 regularization↩︎

We use an invertible flow network to construct the encoding which maps one-to-one between \(X\) and \(H\). The base model follows the framework of [18]. The ResNet has \(L\) residual blocks, and the \(l\)-th block residual mapping \(f(x, \theta_l)\) is parametrized by \(\theta_l\). The overall ResNet mapping from \(X\) to \(H\) can be expressed as \(H=F_\theta(X) = x_L\) where \[\begin{align} \label{eq:resnet} x_{l} & = x_{l-1} + f(x_{l-1}, \theta_l), \quad l=1, \cdots, L, \quad x_0 = X. \end{align}\tag{8}\] In each of the \(L\) residual blocks, we use a shallow network with one or two hidden layers for \(f(x,\theta_l)\). The parameter \(\theta_l\) consists of the weights and bias vectors of the shallow network. The ResNet architecture is illustrated on the right of Fig. 2. The architecture of the residual block is free-form: one can use any layer type inside the residual blocks and even different layer types in different blocks. In this work, we use a fully-connected layer for Euclidean data and GNN layers for graph data. To compute the log determinant in 4 , we express the quantity as a power series and adopt the technique in [29] to obtain an unbiased estimator. Further details can be found in Appendix 9.1.

The invertibility of the ResNet is to be fulfilled by using sufficiently large \(L\) together with the Wasserstein-2 regularization in 3 , which we introduce here. The regularization \(\mathcal{W}\) takes the form as \[\mathcal{W} = \frac{1}{n_{\rm tr}} \sum_{i=1}^{n_{\rm tr}} \ell_w (X_i), \quad \ell_w (X_i)= \sum_{l=1}^L \| x_{l} - x_{l-1} \|_2^2,\] where \(x_0 = X_i\) and \(x_l\) is defined as in 8 . Replacing the empirical measure of \(n_{\rm tr}\) training samples by the data density \(p\) gives the population counterpart of \(\mathcal{W}\) as \(W := \mathbb{E}_{x_0 \sim p} \sum_{l=1}^L \| x_{l} - x_{l-1} \|_2^2\). We denote by \(T_l\) the transport map from \(x_{l-1}\) to \(x_l\), that is, \[\label{eq:def-Tl-resnet} T_{l}(x ) := x + f(x, \theta_l), \quad T_l: \mathbb{R}^d \to \mathbb{R}^d,\tag{9}\] and then \(x_l = T_l (x_{l-1})\). Define the composite of \(T_l\)’s as \(F_l\), which is the transport mapping from \(x_0\) to \(x_l\), namely, \(F_l := T_l\circ \cdots \circ T_1\), \(x_l = F_l( x_0)\), for \(l=1,\cdots, L\), and \(F_0 = \rm{Id}\). We also define \(\rho_l = (F_l)_{\#} p\) which is the marginal distribution of \(x_l\), and \(\rho_0 = p\). Then \(W\) can be equivalently written as \[\begin{align} \label{eq:W-reg-2} W = \sum_{l=1}^L \mathbb{E}_{x \sim \rho_{l-1}} \| T_{l}( x ) - x \|_2^2. \end{align}\tag{10}\] We will show in Section 4.1 that regularizing by \(W\) is equivalent to adding \[W':= \sum_{l=1}^L W_2( \rho_{l-1}, \rho_l )^2\] to the minimizing objective. We thus call the proposed regularization the “Wasserstein-2 regularization,” because \(W'\) sums over the step-wise (squared) Wasserstein-2 distance over the \(L\) steps.

In the limit of large \(L\), the proposed regularization by \(W\) serves to penalize the transport cost. When \(K=1\), the optimal flow induced by the minimizer gives the Wasserstein geodesic from \(p\) to the normal density, cf. Section 4.2. In our problem of conditional generation, when there are \(K>1\) classes, the limiting continuous time flow differs from the Wasserstein geodesic from \(\rho(\cdot, 0) = p(X|Y)\) to \(\rho(\cdot, 1) = p(H|Y)\) for fixed \(Y\), but is expected to give a shared flow in \(\mathbb{R}^d\) from \(X\) to \(H\), cf. Assumption 1 and Fig. 3. Under regularity conditions, one would expect the transport-cost regularized continuous-time flow to also be regular (the \(K=1\) case is proved in Proposition 3). As a result, the velocity field \(v(x,t)\) would have a finite \(x\)-Lipschitz constant \(B\) on any bounded domain. Then, the transport from \(x(t_{l-1})\) to \(x(t_l)\) induced by \(v(x,t)\) on the interval \([t_{l-1}, t_l]\) is invertible when \(B \Delta t < 1\), which holds when \(L > B\). In this case, also assuming that the solved discrete-time transport map \(T_l\) is close to that induced by \(v(x,t)\), one can expect the invertibility of \(T_l\) to hold for each \(l\), and then the composed transport \(F_L = F_\theta\) is also invertible.

We numerically verify the invertibility of the trained \(T_l\) in ResNet in experiments, cf. Table 1 and Appendix 9. Our analysis in Section 4.1 also suggests that the regularization factor \(\gamma\) should scale with \(L = 1/\Delta t\), the number of residual blocks. In experiments, we find that the invertibility of the ResNet can be guaranteed over a range of choices of \(\gamma\), and the quality of generation is insensitive to the choice, see more in Appendix 9.3. The number of blocks \(L\) to use depends on the complication of the data distribution, and in all our experiments, we find a few tens to be enough, including the large graph experiment.

4 Theory: Invertible flow↩︎

In this section, we first interpret the proposed Wasserstein-2 regularization in view of transport cost in Section 4.1. In Section 4.2, utilizing theories of diffusion process and the dynamic formulation of OT, we theoretically study the existence of invertible flow in the simplified case where \(K=1\). All proofs are in Appendix 8.

Throughout the section, we consider flow mapping induced by a velocity field \(v(x,t)\), that is, the continuous-time flow is represented by an initial value problem (IVP) of ODE \[\label{eq:ode-Xt} \dot{x} (t) = v( x(t), t), \quad x(0) \sim p.\tag{11}\] The transport in \(\mathbb{R}^d\) is the solution mapping from \(x(0)=x\) to \(x(t)\) at some time \(t > 0\).

4.1 Interpretation of Wasserstein-2 regularization↩︎

The proposed regularization by \(W\) in Section 3.4 can be interpreted using the dynamic formula of the Wasserstein-2 transport, where we recall the notations in Section 2.2. Suppose the time interval \([0,1]\) is divided into \(L\) time steps, \(\Delta t = t_{l+1} - t_l = {1}/{L}\). With the optimal velocity field \(v(x,t)\) and the corresponding \(\rho(x,t)\) that minimize the transport cost in 2 , on every time subinterval \([t_{l-1}, t_l]\), \(\mathbb{E}_{x(t_{l-1})\sim \rho(\cdot, t_{l-1})} \int_{t_{l-1}}^{t_l} \| v( x(t), t) \|^2 dt = { W_2( \rho( \cdot, t_{l-1} ), \rho( \cdot, t_{l} ) )^2 }/{\Delta t}\). By the definition of the transport cost 1 , we have \[\begin{align} \mathcal{T} = & \sum_{l=1}^L { W_2( \rho( \cdot, t_{l-1} ), \rho( \cdot, t_{l} ) )^2 }/{\Delta t}. \label{eq:T-by-rhol} \end{align}\tag{12}\] Note that the right hand side (r.h.s.) only depends on \(\rho(\cdot, t)\) at the time stamps \(t_l\). We define \(\rho_l: = \rho(\cdot, t_l)\), and denote by \(T_l\) the transport map from \(x(t_{l-1})\) to \(x(t_l)\) induced by the ODE of \(x(t)\) on the time interval \([t_{l-1}, t_l]\), i.e., \(T_l(x) = x + \int_{t_{l-1}}^{t_l} v(x(s),s) ds\) starting from \(x(t_{l-1}) = x\). About the notation, \(\rho_l\) and \(T_l\) here are determined by the continuous-time flow, and the notations coincide with those in 9 and 10 which are determined by the finite-step ResNet. We use the same notations here since \(T_l\) will be the variables to minimize the training objective, and thus ready to be parametrized by a residual block \(f(\cdot, \theta_l)\).

Specifically, when the other part in the loss, denoted as \(\mathcal{L}\), only depends on the terminal time density \(\rho(\cdot, 1)=\rho(\cdot, t_L)\), one can use the r.h.s. of 12 as the regularization term added to \(\mathcal{L}\), making the overall objective as \[\label{eq:L43T} \mathcal{L} + \mathcal{T} = \mathcal{L} [ \rho_L ] + \frac{1}{\Delta t} \sum_{l=1}^L W_2( \rho_{l-1}, \rho_l )^2,\tag{13}\] where \(\rho_l = (T_l)_{\#} \rho_{l-1}\), \(\rho_0 = p\). Note that for a given \(p\), the densities \(\rho_1, \cdots, \rho_L\) are determined by the \(L\) transport maps \(T_l\)’s, this means that \(T_1, \cdots, T_L\) can be used as the variables to minimize 13 .

The following proposition shows that this minimization is equivalent to using the proposed Wasserstein-2 regularization in Section 3.4 (specifically, using \(1/\Delta t\) times the r.h.s. of 10 ) in solving for \(T_l\)’s, and there is no need to solve for the OT distance \(W_2( \rho_{l-1}, \rho_l )^2\) with additional computation.

4.2 Existence and invertibility of flow↩︎

The transport mapping by the IVP 11 is invertible as long as 11 is well-posed. For the \(K\)-class case, Assumption 1 assumes the existence of the smooth velocity field \(v(x,t)\) as well as the shared flow mapping. Here we consider the special case of \(K=1\) (the unconditional generation problem), that is, the source density \(p\) is the data distribution in \(\mathbb{R}^d\) and the target density \(q\) is standard normal \(\mathcal{N}(0, I_d)\). Though the \(K\)-class conditional generation problem is of the main interest of the paper and has called for additional requirement of the share flow, our analysis of the \(K=1\) case provides theoretical insights, especially for the expressiveness of GNN layers in the flow model for graph data in Section 5.

The normalizing flow induced by \(v(x,t)\) from \(p\) to normal \(q\) is typically not unique. We introduce two constructions here that are related to the proposed Wasserstein-2 regularization of the flow model.

(i) Flow by Benamou-Brenier formula. Consider 11 on \(t\in [0,1]\), the flow \(F\) maps from \(x=x(0)\) to \(x(1)\). Let \(\rho(\cdot, t)\) be the density of \(x(t)\), \(\rho(\cdot, 0)= p\) and \(\rho(\cdot, 1) = F_\# p\). We consider the flow induced by the optimal \(v(x,t)\) in the Benamou-Brenier formula 2 , where the regularity of \(v\) and \(\rho\) follows from classical OT theory:

We now show that our optimization objective 3 is equivalent to the action minimization in the Benamou-Brenier formula when \(L\) is large so that the discrete-time transports approximate the continuous-time limit. When \(K=1\), the generalization loss (in population form) in 3 reduces to \[\mathcal{L}_g = - \mathbb{E}_{ x \sim p} \log ( (F^{-1})_\# q ( x ) ).\] By the relation [38] \[\label{eq:reverse-KL} {\rm KL}( F_\# p || q) = {\rm KL}( p || (F^{-1})_\# q),\tag{14}\] and that \({\rm KL}( p || (F^{-1})_\# q) = \mathcal{L}_g + \int p \log p\), we have that \(\mathcal{L}_g + \int p \log p = {\rm KL}( F_\# p || q)\), where \(\int p \log p\) is a constant independent of the model. The analysis in Section 3.4 gives that \(\mathcal{W}\) (in population form) is equivalent to \(\frac{1}{L} \mathcal{T}\) with large \(L\). Thus our minimizing objective \(\mathcal{L}_g + \gamma \mathcal{W}\) is equivalent to \[\label{eq:BB-loss} \mathcal{T} + \tilde{\gamma} {\rm KL}( \rho( \cdot, 1) || q )\tag{15}\] with some positive scalar \(\tilde{\gamma}\). Compared to the Benamou-Brenier formula 2 , the objective 15 relaxes the terminal condition that \(\rho(\cdot, 1) = q\) to be the KL divergence, which does not change the solution of the optimal \(v\).

(ii) Flow by Fokker-Planck equation. Because the transport cost \(\mathcal{T}\) equals squared Wasserstein-2 distance between \(\rho(\cdot,0)=p\) and \(\rho(\cdot,1)\) at optimal \(v\), the objective 15 is closely related to the problem of \(\min_{\tilde{q}} {\rm KL}( \tilde{q} || q ) + \frac{1}{\tilde{\gamma}} W_2( p, \tilde{q})^2\). When \(\tilde{\gamma}\) is small, the problem is the first step of the Jordan-Kinderleherer-Otto (JKO) scheme [39] to solve the Fokker-Planck equation of a stochastic diffusion process toward the equilibrium \(q\). Because \(q\) is standard normal, the stochastic process is an Ornstein-Uhlenbeck (OU) process. Without going to further connections to the JKO scheme, here we use the Fokker-Planck equation to provide another theoretical construction of the invertible normalizing flow. For the OU process in \(\mathbb{R}^d\), the Fokker-Planck equation can be written as \[\label{eq:fokker-planck} \partial_{t}\rho=\nabla\cdot(\rho\nabla V + \nabla \rho), \, V(x)={|x|^{2}}/{2}, \, \rho(x,0) = p(x),\tag{16}\] where \(\rho(x,t)\) represents the probability density of the OU process at time \(t\). The Liouville equation of 11 is \(\partial_{t}\rho=-\nabla\cdot(\rho v)\), where \(\rho(x,t)\) represents the density of \(x(t)\). Comparing to 16 , we see that the density evolution can be made the same if the velocity field \(v(x,t)\) is set to satisfy \[\label{eq:f-jko-flow} - v (x,t) = \nabla V(x) + \nabla \log \rho( x,t) = x + \nabla \log \rho( x,t).\tag{17}\] The smoothness of \(v\) follows from that of \(\rho(x,t)\), which has an explicit expression as the solution of 16 .

While theoretically the density \(\rho(x,t)\) in 16 converges to the normal equilibrium \(q\) in infinite time, the convergence is exponentially fast [40]. Thus the transported density \(F_\# p = \rho(\cdot, T)\) for a finite \(T \sim \log (1/\epsilon)\) can be \(\epsilon\)-close to \(q\), which then guarantees the closeness of \((F^{-1})_\# q\) to \(p\) by 14 .

5 Theory: Expressing flow of graph data↩︎

We analyze the expressiveness of graph convolution layers to approximate the theoretical flow maps identified in Section 4. All proofs are in Appendix 8.

5.1 Theoretical flows of graph Gaussian field data↩︎

Taking the continuous-time formulation, the goal is to represent the velocity field \(v(x,t)\) by a GNN layer where \(x\) is graph data. Assuming that the discrete-time flow can approximate the continuous-time counterpart when the number of residual blocks is large, the successful approximation of \(v (x,t)\) by GNN layers indicates that the flow can be constructed by the invertible GNN flow network.

For simplicity, we consider \(x(0) \sim \mathcal{N}( 0, \Sigma)\), which is a Gaussian field on the graph having \(N\) nodes (The node feature dimension \(d'=1\), and data dimension \(d=N\)). The \(N\)-by-\(N\) covariance matrix \(\Sigma\) is positive semi-definite, and we further assume it is invertible. Suppose PSD matrix \(A = VDV^T\) is the eigen-decomposition, where \(V\) is an orthogonal matrix, define \(A^{1/2} = V D^{1/2} V^T\). In this special case of Gaussian data, we have explicit expressions of the velocity field of the flow induced by Benamou-Brenier formula and by Fokker-Planck equation, as characterized by the following two lemmas.

In particular, at \(t=0\), \(T_0^{\rm FK} + I = \Sigma^{-1}\). As \(t \to \infty\), \(T_t^{\rm FK} \to 0\) exponentially fast.

In particular, \(T_0^{\rm BB} = \Sigma^{-1/2} - I\), and \(T_1^{\rm BB} = I - \Sigma^{1/2}\).

The velocity field \(v(x,t)\) in both cases is a linear mapping of \(x\) depending on \(t\), and thus it can always be exactly expressed by a fully-connected layer. For GNN layers, the problem is to represent the matrix \(T_t\) in both cases by a graph convolution. Because the desired operator itself is linear, we only consider the graph convolution in space (there is no channel mixing parameter because the node feature dimension \(d'=1\) here), omitting the bias vector and the non-linear activation.

5.2 Spectrum of \(\Sigma\) and an approximation lemma↩︎

The approximation error will depend on the condition number of \(\Sigma\), which reveals the fundamental difficulty of approximation when \(\Sigma\) is near singular (the operators \(T_t^{\rm FK}\) and \(T_t^{\rm BB}\) involves \(\Sigma^{-1}\) at \(\Sigma^{-1/2}\) at \(t=0\) respectively). Define \(S:= \Sigma^{-1}\) as the precision matrix of the distribution of \(x(0)\). Denote the condition number of \(\Sigma\) as \(\kappa := {\rm Cond}(\Sigma) < \infty\). \(\kappa = d_{\rm max}/d_{\rm min} \ge 1\), where \(d_{\rm max}\) and \(d_{\rm min}\) are the largest and smallest eigenvalues of \(\Sigma\). Suppose the data is properly normalized and without loss of generality, we divide \(\Sigma\) by the constant \(\sqrt{ d_{\rm max} d_{\rm min} }\), which makes the smallest and largest eigenvalue equal \(1/\sqrt{\kappa}\) and \(\sqrt{\kappa}\) respectively. Denote the spectrum of a real symmetric matrix \(A\) (i.e., the set of \(n\) eigenvalues) as \({\rm spec}(A)\). Because the eigenvalues of \(S\) are the reciprocal of those of \(\Sigma\), we have made \[\label{eq:spec-Sig-S} {\rm spec}( \Sigma ), {\rm spec}( S ) \subset [ 1/\sqrt{\kappa}, \sqrt{\kappa} ].\tag{18}\] Meanwhile, we also have \[\label{eq:spec-Sig-S-1472} {\rm spec}( \Sigma^{1/2} ), {\rm spec}( S^{1/2} ) \subset [ 1/\sqrt{\kappa'}, \sqrt{\kappa'} ], \quad \kappa' := \sqrt{\kappa}.\tag{19}\]

We introduce a function approximation lemma which will be used in the analysis of both \(T_t^{\rm FK}\) and \(T_t^{\rm BB}\). Define \[\label{eq:def-fs-gs} \begin{align} f_s(x) &: = \left( s + (1-s) x \right)^{-1},\\ g_s(x) &: = \left( s + (1-s) x^{-1} \right)^{-1}, \end{align} \quad 0 \le s \le 1.\tag{20}\]

In below, we first introduce the approximation of \(T_t^{\rm FK}\) by spectral and spatial graph convolutions. The approximation of \(T_t^{\rm BB}\) uses similar techniques and will be included afterward.

5.3 Approximating \(T_t^{\rm FK}\) by spectral graph convolution↩︎

By Lemma 2, with \(s =1 - e^{-2t} \in [0, 1)\), we have \(T_t^{\rm FK} + I = ( s I+ (1-s) \Sigma )^{-1} =: { T}_t'\), and it is equivalent to approximate \({ T}_t'\) by a graph convolution. By definition, \[\label{eq:Tt-by-fs-gs} { T}_t' = f_s(\Sigma) = g_s (S).\tag{21}\] Denote by \(\tilde{L}\) the (possibly normalized) graph Laplacian matrix, which is a real symmetric PSD matrix.

The assumption means that both \(\Sigma\) and \(S=\Sigma^{-1}\) can be approximated by some low degree polynomial of \(\tilde{L}\) up to relative error \(\delta\). This can be the case, for example, when \(\Sigma\) shares the eigenvectors with \(\tilde{L}\), and the polynomial \(p\) and \(q\) can be constructed to fit the eigenvalues. The condition (C1) may be relaxed by enlarging the interval \([ 1/\sqrt{\kappa}, \sqrt{\kappa} ]\) to be \([ 1/a, a ]\) where \(a\) is a multiple of \(\sqrt{\kappa}\), and then construct the polynomial approximation of \(g_s\) and \(f_s\) on \([ 1/a, a ]\) as in Lemma 4, which results in an \(O(1)\) constant factor in the bound. In this case, (C1) may also be induced by (C2) with \(O(1)\) constant \(\delta\). We keep the current form of (C1)(C2) for simplicity.

The following theorem constructs the throughout-time approximation of \({ T}_t^{\rm FK}\) by using the composed polynomials \(Q_s^{(n+1)} \circ q\) and \(P_s^{(n)}\circ p\) on short and long times respectively.

The theorem suggests that when the covariance matrix \(\Sigma\) is is well-conditioned, then small approximation error can be achieved using a low-degree \(n\).

5.4 Approximating \(T_t^{\rm FK}\) by spatial graph convolution↩︎

We consider approximating operator \(T_t'\) by spatially local graph filters, such as L3Net layer [35] as introduced in Section 3.3. We start by introducing the notion of locality on graph.

The neighborhood \(\mathcal{N}^{(v)}_i\) denotes all nodes accessible from node \(i\) within \(v\) steps along adjacent nodes. By definition, if a matrix \(B\) is \(v\)-local, then \(B^k\) is (\(kv\))-local.

The spatial graph convolution operator can be written as \(\sum_{r=1}^R c_r B_r\), where \(B_r\) are local filters on the graph, and \(c_r\) are coefficients. The following theorem proves the approximation of \(T_t'\), achieving the same bound as in Theorem 5. The construction is by using powers of the local filters \(\tilde{\Sigma}\) and \(\tilde{S}\) as \(B_r\).

Note that the construction has time-independent \(B_r\), which means that in the L3Net GNN layers, we can potentially share the basis filter \(B_r\) across residual blocks throughout the \(L\) blocks, which will further reduce model complexity.

5.5 Approximating \(T_t^{\rm BB}\) by graph convolutions↩︎

We take \(s = t\) the time on \([0,1]\). By Lemma 3, \(T_t^{\rm BB}\) has the following short and long time representation, \[\begin{align} T_t^{\rm BB} + \frac{1}{1-t} I &= \frac{1}{1-t} g_t( S^{1/2}), \quad t \in [0,1/2], \\ T_t^{\rm BB} - \frac{1}{t} I &= - \frac{1}{t} \bar{f}_t (\Sigma^{1/2}), \quad t \in [1/2,1], \end{align}\] where \(\bar{f}_s (x) : = x f_s(x)\), and \(f_s\), \(g_s\) are defined same as before 20 . Let \(\tilde{\Sigma}^{1/2}\) and \(\tilde{S}^{1/2}\) are the approximators by some low-degree polynomials of \(\tilde{L}\) in spectral convolution, and by \(v\)-local graph filters in spatial convolution, respectively. Following 19 , we replace \(\kappa\) to be \(\kappa'\) in the analysis, and the conditions (C1)(C2) in Assumptions 2-3 become

(C1’) \({\rm spec}( \tilde{\Sigma}^{1/2} ), {\rm spec}( \tilde{S}^{1/2} ) \subset [ 1/\sqrt{\kappa'}, \sqrt{\kappa'} ]\).

(C2’) \(\| \tilde{\Sigma}^{1/2} - {\Sigma}^{1/2} \|_{op}, \| \tilde{S}^{1/2} - S^{1/2} \|_{op} \le \delta \sqrt{\kappa'}\).

Noting that in the short and long time expression of \(T_t^{\rm BB}\), the constant factor \(1/(1-t)\) and \(1/t\) are bounded by 2 respectively, thus we aim to approximate \(g_t( S^{1/2})\) and \(\bar{f}_t (\Sigma^{1/2})\). In applying Lemma 4(ii), the bound becomes \(| \bar{f}_s(x) - x P_s^{(n)}(x) | \le a b (b/(b+1))^n\), which raises the polynomial degree to be \((n+1)\), and has another factor of \(a = \sqrt{\kappa'}\) in the bound. The analysis proceeds using the same technique, and the details are omitted. The final approximation bounds in both the short-time (\(t \in [0,1/2]\)) and long-time (\(t \in [1/2,1]\)) cases are \[4 \delta\sqrt{\kappa'} + 2 \kappa' \exp\{ - \frac{n}{\sqrt{\kappa'}+1} \} \sim \delta \kappa^{1/4} + \kappa^{1/2} \exp\{ - \frac{n}{ \kappa^{1/4}+1} \},\] which improves the dependence on condition number from \(\kappa^{1/2}\) to \(\kappa^{1/4}\).

5.6 Comparison of spectral and spatial graph convolutions↩︎

The analysis above shows that, for Gaussian field data \(X\) on graph, the ability of a GNN flow network to learn the normalizing flow is determined by the ability of the graph convolution filters to approximate the covariance matrix and precision matrix (or their square roots) of \(X\). In terms of the expressiveness of spectral and spatial filters, first note that when the spectral graph convolution filters are local (e.g., ChebNet filters of low polynomial degree), they become a special case of the spatial filters. In this case, spatial filters are always more expressive. Meanwhile, there can be local spatial filters that cannot be expressed by spectral filters [35]. Here we provide an example of covariance matrix \(\Sigma\) on a three-node graph that cannot be represented by spectral graph filter due to constructed symmetry.

This example will be empirically examined in Section 6.4 where we train iGNN models on simulated graph data, see Fig. 14. We find that for this example, using ChebNet layer fails to learn the generation flow, and switching to the L3Net layer resolves the issue. In other situations where spectral graph filters have sufficient expressiveness for graph data, the spectral GNN has the advantage of a lighter model size.

6 Experiment↩︎

We first examine the proposed iGNN model on simulated data, including large-graph data, in Section 6.2. We then apply the iGNN model to real graph data (solar ramping event data and traffic flow anomaly detection) in Section 6.3. In Section 6.4, we study the effect of using different GNN layers in the iGNN model. The experimental setup is introduced in Section 6.1, and further details and additional results are provided in Appendix 9.2. The code is available at https://github.com/hamrel-cxu/Invertible-Graph-Neural-Network-iGNN.

6.1 Experiment setup↩︎

Baselines and evaluation metrics. We consider three competing conditional generative models: (a) Conditional generative adversarial network (cGAN) [14]. (b) Conditional invertible neural network with maximum mean discrepancy (cINN-MMD) [15]. (c) Conditional invertible neural network using normalizing flow (cINN-Flow) [16]. To quantify performance, we measure the difference between two distributions (\(X|Y\) versus \(\hat{X}|Y\) at different \(Y\)) by kernel maximum mean discrepancy (MMD) [41] and energy statistics [42]. Details of the MMD and energy statistics metrics are contained in Appendix 9.1. We also provide qualitative comparison by visualizing the distribution of generated data.

Data and ResNet architecture. In the examples of graph data, the number of graph nodes ranges from 3 to 500. All graphs are undirected and unweighted, with inserted self-loops. Regarding the ResNet block layer type: for non-graph data, we use 2 fully-connected hidden layers of 64 neurons in all ResNet blocks. For graph data, in each residual block, we replace the first hidden layer to be a GNN layer (either ChebNet or L3Net layer); the second layer is a shared fully-connected layer that applies channel mixing across all the graph nodes (which can be viewed as a GNN layer with identity spatial convolution). The activation function is chosen as ELU [43] or LipSwish [29], which have continuous derivatives. The network is trained end-to-end with the Adam optimizer [44]. Hyperparameter selection and additional results with different choices are described in Appendix 9.3.

6.2 Simulated examples↩︎

We consider three simulated examples in this section. The first considers non-graph Euclidean data, and the second and the third consider graph data on small and large graphs. Detailed data-generating procedures can be found in Appendix 9.1.

1. Non-graph data in \(\mathbb{R}^2\). The dataset contains a Gaussian mixture of eight components with four classes, where each \(X|Y\) is further divided into two separated Gaussian distributions in \(\mathbb{R}^2\). Fig. 4 compares the generative results of the proposed iGNN model with cINN-MMD, where both methods can generate data that are reasonably close to \(X|Y\) at each \(Y\).

2. Data on a small graph. We consider the three-node graph introduced in Example 1. At each node \(v\), \(Y_v \in \{0,1\}\) and features \(X_v\in \mathbb{R}^2\), thus the graph node label vector \(Y\in \{0,1\}^3\). In this example, iGNN yields comparable generative performance with cINN-MMD, where both models use L3net GNN layers. See Fig. 13 and more results in appendix 9.2.

3. Data on a large graph. To demonstrate the scalability of our approach, we consider a 503-node chordal cycle graph [45], which is an expander graph. We design binary node labels and let node features \(X|Y \sim N(\mu_Y,\Sigma_Y)\), where \(X_v \in \mathbb{R}^2\) and the mean \(\mu_Y\) and covariance matrix \(\Sigma_Y\) contain graph information. Because enumerating all values of \(Y\) is infeasible, we randomly choose 50 values of outcome \(Y\), each of which has 50% randomly selected node labels to be 1. To visualize the generative performance of iGNN, we compare the covariance of true and generated data restricted to subgraphs. Specifically, we plot the covariance matrix of model-generated data \(\hat{X}|Y\) and true data \(X|Y\) on sub-graphs produced by 1 or 2-hop neighborhoods of a graph node. Fig. 5 shows the resemblance between learned and true covariance matrices on different neighborhoods on the graph.

6.3 Real-data examples↩︎

We apply the iGNN model to two graph prediction data in real applications, the solar ramping event data, and the traffic flow anomaly detection data. The inverse prediction problem is formulated as a conditional generation task.

1. Solar ramping events data. Consider the anomaly detection task on California solar data in 2017 and 2018, which were collected in ten downtown locations representing network nodes. Each node records non-negative bi-hourly radiation recordings measured in Global Horizontal Irradiance (GHI). After pre-processing, graph nodal features \(X_t \in \mathbb{R}^{10\times 2}\) denote the average of raw radiation recordings every 12 hours in the past 24 hours, and response vectors \(Y_t \in \{0,1\}^{10}\) contain the anomaly status of each city. Fig. 6 shows that the learned conditional distribution \(\hat{X}|Y\) by iGNN model closely resembles that of the true data \(X|Y\), and outperforms the generation of cINN-MMD. The quantitative evaluation is given in Table 2, which shows that iGNN has comparable or better performance than the alternative approaches (smaller test statistics indicate better generation). The table also shows that cINN-Flowperforms significantly worse than cINN-MMDand iGNN on this example, which is consistent with the visual comparison of \(\hat{X}|Y\) (not shown). Lastly, Fig. 7 shows the predictive capability of iGNN: given a test node feature matrix \(X\), we can compute \(\mathbb{P}(Y_i=1|X)\) for node \(i\) using the trained linear classifier on \(F_{\Theta}(X)\). The predicted probabilities learned by the model are consistent with the true nodal labels, and provide more information than the binary prediction output.

2. Traffic flow anomalies. We study the anomaly detection task on Los Angeles traffic flow data from April to September 2019. The whole network has 15 sensors with hourly recordings. Graph nodal features \(X_t \in \mathbb{R}^{15 \times 2}\) denote the raw hourly recording in the past two hours, and response vectors \(Y_t \in \{0,1\}^{15}\) contain the anomaly status of each traffic sensor. The graph topology is shown in Fig. [traffic95network], along with the raw input features in \(\mathbb{R}^2\) (over all graph nodes). The generated data distribution by iGNN resembles the ground truth, as shown in Fig. 8 (b), and the performance is better than that of cINN-MMDin (c). The quantitative evaluation metrics also reveal the better performance of iGNN over alternative baselines, cf. Table 2.

6.4 Comparison of GNN layers in iGNN models↩︎

We examine the empirical performance of different GNN layers in learning the normalizing flow of graph data, for which the theoretical analyses appeared in Section 5.

1. Simulated data on three-node graph. To validate our theory, we study two simulated datasets to show the possible insufficiency of spectral GNN layers, as has been explained in Section 5.6. The graph data are on a three-node graph: the first one is as in Example 1, where nodal feature dimension \(d'=1\), and \(K=1\); the second one has nodal feature dimension \(d'=2\), and \(K=2\), see more details in Appendix 9.1. We compare the generative performance of iGNN by using ChebNet and L3Net layers for both examples. The result on the \(d'=2\) example is shown in Fig. 9, where iGNN with ChebNet layers fails to learn the conditional distribution \(X|Y\), cf. plot (b). Meanwhile, iGNN with L3Net layers yields satisfactory performance by having sufficient model expressiveness for this example. The second experiment studies Example 1 considered in Section 5.6. The results are shown in Fig. 14, where, similarly, iGNN with ChebNet layers fails to generate the graph data \(X\) and switching to L3Net layers resolves the insufficiency of expressiveness.

2. Simulated large-graph data. We also compare the GNN layers on a larger graph using the same 503-node chordal cycle graph as in Fig. 5, and simulate \(X\sim \mathcal{N}(0,\Sigma)\), where \(\Sigma^{-1}=\sum_{k=0}^2 a_k T_k(\widetilde{L})\), \(T_k\) is the \(k\)-th Chebyshev polynomial, \(\widetilde{L} = I - D_A^{-1/2}AD_A^{-1/2}\), \(A\) being the adjacency matrix and \(D_A\) the degree matrix, and \(a_0={\color{black}0.6},a_1=0,a_2=0.5\). Fig. [Cheb95generation] (ChebNet) and [Local95generation] (L3Net) in appendix 9.2 show that the iGNN model learns to generate samples close in distribution to the ground truth, as reflected in the resemblance of the correlation matrices. For this dataset, both spectral and spatial GNN layers have sufficient expressiveness to learn the data distribution.

7 Discussion↩︎

In this work, we developed the iGNN model, a conditional generative deep model based on invertible normalizing flow ResNet for inverse graph prediction problems. The model encodes the conditional distribution \(X|Y\) into a parametric distribution \(H|Y\) by a one-to-one mapping from \(X\) to \(H\), which allows forward prediction and generation of input data \(X\) given an outcome \(Y\) (including the construction of uncertainty sets for \(X|Y\)). The scalability of the iGNN model to large graphs is achieved by taking a factorized formula of \(H|Y\) with a provable component separation guarantee, and the computational scalability is achieved via adopting GNN layers. In addition, the invertibility of the flow network is ensured by a Wasserstein-2 regularization, which can be computed efficiently in the forward pass of the flow network and is compatible with free-form neural network layer types in the residual blocks. Theoretically, we analyzed the existence of invertible flows and examined the expressiveness of GNN layers on graph data to express these normalizing flows. In experiments, we showed the improved performance of iGNN over alternatives on real data and the scalability of the iGNN model on large graphs.

There are several future directions to extend the work. On the theoretical side, first, under the framework of our problem, more analysis of the \(K >1\) case will be useful to go beyond the current Assumption 1. Second, it would be interesting to develop a full approximation result on how the theoretical continuous-time flow identified in Section 4 can be approximated by a deep residual network. Specifically, one can construct a discrete-time flow on \([0, T]\) with \(L\) time stamps \(t_l\) to approximate the continuous-time one. \(T=1\) for the Benamou-Brenier flow, and \(T \sim \log ({1}/{\epsilon})\) for the Fokker-Planck flow (so as to achieve \(\epsilon\)-closeness to the normal density at time \(T\)). The smoothness of the velocity field in both cases, cf. Propositions 3 and 4, implies that at each time stamp \(t_l\), \(v(x, t_l)\) may be approximated by a shallow residual block with finite many trainable parameters based on universal approximation theory. Compositing the \(L\) time steps constructs a neural network approximation of the flow with a provable approximate generation of the data density \(p\). The precise analysis is left to future work. At last, one may be curious about the connection between the encoding scheme of iGNN model and the Gaussian channel. If we view \(\mu_k\) as the input of the channel and the output as given by \(\mu_k + Z\), for noise \(Z\sim \mathcal{N}(0, \sigma^2 I_d)\), then the separation of \(K\) components can be viewed as a requirement to ensure the error probability of the linear decoder to be sufficiently small. A natural question arises in how many components can be encoded in this way and subsequently decoded (classified) successfully, for instance, by considering a compact encoding domain. To further develop the methodology, one interesting question is to consider the regression problem, where the outcome \(Y\) takes continuous (rather than categorical) values. Another related question is the inverse of the graph classification problem, that is, given an output label of a whole graph, generate the input graph data including changing graph topology and node/edge features. The current work tackles the node classification setting, and some of our techniques may extend to the graph classification setting.

Acknowledgement↩︎

The work is supported by NSF DMS-2134037. C.X. and Y.X. are supported by an NSF CAREER Award CCF-1650913, NSF DMS-2134037, CMMI-2015787, DMS-1938106, and DMS-1830210. X.C. is partially supported by NSF, NIH and the Alfred P. Sloan Foundation.

8 Proofs↩︎

8.1 Proofs in Section 3↩︎

8.2 Proofs in Section 4↩︎

8.3 Proofs in Section 5↩︎

9 Additional experimental details↩︎

9.1 Experimental set-up↩︎

9.1.1 Computation of \(\log\det\)↩︎

To compute the log determinant in 4 , we adopt the following unbiased log determinant approximation technique as proposed in [29]. Let \(F_{\theta}(x)=x+f_{\theta}(x)\) denote the output from a generic ResNet block with parameter \(\theta\). First, observe that for any input \(x\), \(\log|\det J_{F_{\theta}}(x)|=\text{tr}(\log J_{F_{\theta}}(x))\) because the matrix \(J_{F_{\theta}}(x)\) is non-singular. We thus have \(\text{tr}(\log J_{F_{\theta}}(x))=\text{tr}(\log (I+J_{f_{\theta}}(x)))\). As a result, the trace of the matrix logarithm can be expressed as \[\label{matrix95log} \text{tr}(\log (I+J_{f_{\theta}}(x))) = \text{tr}\left(\sum_{k=1}^{\infty} \frac{(-1)^{k+1}}{k}[J_{f_{\theta}}(x)]^k\right).\tag{41}\] Based on 41 , which takes infinite time to compute, we can obtain an unbiased estimator in finite time based on the “Russian roulette” estimator approach [29], [48]. For a ResNet as a concatenation of \(L\) ResNet blocks, the approximation is applied to each block and summed over all blocks. Lastly, to speed up gradient computation of the approximation, we further adopt memory efficient backpropagation through the early computation of gradients [29].

9.1.2 Model evaluation metrics↩︎

MMD statistics metric. Given two sets of samples \(\boldsymbol{X}=\{x_1,\ldots,x_n\},\boldsymbol{X'}=\{x'_1,\ldots,x'_n\}\) of same sample size \(n\), the MMD two-sample statistic between \(\boldsymbol{X}\) and \(\boldsymbol{X'}\) is defined as \[\label{MMD} \text{MMD}(\boldsymbol{X},\boldsymbol{X'}):=\frac{1}{n^2}\sum_{i=1}^n \sum_{j=1}^n k(x_i,x_j) + \frac{1}{n^2}\sum_{i=1}^n \sum_{j=1}^n k(x'_i,x'_j) - \frac{2}{n^2}\sum_{i=1}^n \sum_{j=1}^n k(x_i,x'_j),\tag{42}\] where we use the radial basis kernel \(k(x,x')=\exp(-\alpha\|x-x'\|^2)\) with \(\alpha=0.1\).

For the \(K\)-class conditional distribution \(\boldsymbol{X}|Y\), where we denote by \(\{ \boldsymbol{X}|Y=k \}\) the set of samples \(\{x_i:y_i=k\}_{i=1}^n\), the overall MMD statistic is defined using 42 as \[\label{weighted95MMD} {\rm MMD} = \sum_{k=1}^K w_k \text{MMD}( \{ \boldsymbol{X}|Y=k \} , \{ \boldsymbol{X'}|Y=k \}), \quad w_k = \frac{\sum_{i=1}^n \boldsymbol{1}(y_i=k)}{n}.\tag{43}\] Note that on graph data where \(Y\) concatenates all nodal labels, the summation is over all types of \(Y\) (up to \(K^N\) many).

Energy statistic metric. Given two sets of samples \(\boldsymbol{X}=\{x_1,\ldots,x_n\},\boldsymbol{X'}=\{x'_1,\ldots,x'_n\}\), The energy statistic under \(\ell_2\) norm is defined as \[\label{energy} \text{Energy}(\boldsymbol{X},\boldsymbol{X'}):=\frac{2}{n^2}\sum_{i=1}^n \sum_{j=1}^n \|x_i-x'_j\|_2 - \frac{1}{n^2}\sum_{i=1}^n \sum_{j=1}^n \|x_i-x_j\|_2- \frac{1}{n^2}\sum_{i=1}^n \sum_{j=1}^n \|x'_i-x'_j\|_2,\tag{44}\] For \(K\)-class conditional distribution \(\boldsymbol{X}|Y\), the weighted energy statistics is defined using 44 as \[\label{weighted95energy} {\rm Energy} = \sum_{k=1}^K w_k \text{Energy}( \{\boldsymbol{X}|Y=k\} , \{ \boldsymbol{X'}|Y=k\} ), \quad w_k = \frac{\sum_{i=1}^n \boldsymbol{1}(y_i=k)}{n}.\tag{45}\] Computation of the weighted statistics on graph data is identical to that of the weighted MMD statistics on graph.

Model invertibility error. We see from Fig. 11 that iGNN under the Wasserstein-2 regularization ensures model invertibility up to very high accuracy.

9.1.3 Construction of simulated graph data↩︎

We first describe the construction of simulated graph data on small graphs (corresponding to Figures 13 and 9). Each node has a binary label so that \(Y\in \{0,1\}^3\) is a binary vector. Conditioning on a specific binary vector \(Y\) out of the eight choices, the distribution of \(X|Y\) is defined as \[X|Y:=P_A\left(Z|Y+\begin{bmatrix} -4 & 0 & 4\\ 0 & 0 & 0 \end{bmatrix}\right),\] where \(P_A:=D_A^{-1}A\) is the graph averaging matrix.

The distribution of \(Z|Y\) is specified as:

• In Figure 9, we let \[\label{mixture95small} Z_v|Y_v \sim \begin{cases} \mathcal{N}((0,1.5)^T, 0.1I_2) & \text{if }Y_v=0 \\ \mathcal{N}((0,-1.5)^T, 0.1I_2) & \text{if } Y_v=1 \end{cases}.\tag{46}\]

• In Figure 13, \(Z_v|Y_v\) are non-identically distributed over nodes \(v=0,1,2\). When \(v=1\), \(Z_v|Y_v\) follow the same distribution as 46 . When \(v=0 \text{ or } 2\), \(Z_v|Y_v\) are rotated and shifted noisy two moons, whose positions depend on the binary value of \(Y_v\).

On large graphs (corresponding to Fig. 5), each node has a binary label so that \(Y\in \{0,1\}^{503}\) is a binary vector. Conditioning on a specific binary vector \(Y\) out of the eight choices, the distribution of \(X|Y\) is defined as \[X|Y:=R((1-\delta)I+\delta P_A)Z|Y \text{ s.t. } Z_v|Y_v \sim \begin{cases} \mathcal{N}((0,12)^T, I_2) & \text{if }Y_v=0 \\ \mathcal{N}((0,0)^T, I_2) & \text{if } Y_v=1 \end{cases},\] where \(R\) denotes a counter-clockwise rotation matrix for 90 degrees (applied node-wise to each two-dimensional nodal feature) and \(P_A:=D_A^{-1}A\) is the graph averaging matrix, where \(D_A\) is the degree matrix. We choose \(\delta=0.2\) so that a soft graph averaging is applied to the hidden variables \(Z|Y\).

9.2 Additional experimental results↩︎

This subsection contains experimental results to augment those in the main text. In particular,

• We show the confidence region on the three-moon non-graph conditional generation data in Fig. 12. The generative quality is presented in Fig. 3 in the main text.

• We compare the generative quality of iGNN and cINN-MMDon simulated three-node graph conditional generation data in Fig. 13, to augment the results in Section 6.2.

• We numerically verify Example 1 based on simulated three-node graph unconditional generation data. The generative results are presented in Fig. 14 for iGNN using Chebnet and L3net.

• We demonstrate the unconditional generation results of iGNN using GNN layers in Fig. [Cheb95generation] (spectral GNN) and Fig. [Local95generation] (spatial GNN). Empirical results are presented in Section 6.4 and theoretical analyses are presented in Section 5.

9.3 Hyperparameter selection↩︎

We first summarize the hyperparameters used in our experiments and then verify that iGNN is insensitive to alternative hyperparameter choices.

• Mixture distribution in \(\mathbb{R}^2\) with 8 components (cf. Fig. 4): we let iGNN contain 40 ResNet blocks, each of which is built with fully-connected layers with two hidden layers. We fix the learning rate at 5e-4 and train with a batch size of 1000. We fix the \(\gamma\) regularization factor of the \(W_2\) loss as 1.

• Conditional generation on small graph (cf. Fig. 9 and Fig. 13): we let iGNN contain 40 ResNet blocks, where the first hidden layer of each block is an L3Net layer. We fix the learning rate at 5e-4 and train with a batch size of 100. We fix the \(\gamma\) regularization factor of the \(W_2\) loss as 1.

• Conditional generation on large graph (cf. Fig. 5): we let iGNN contain 5 ResNet blocks, where the first hidden layer of each block is an L3Net layer. We fix the learning rate at 1e-3 and train with a batch size of 100. We fix the \(\gamma\) regularization factor of the \(W_2\) loss as 0.375/503.

• Solar ramping event (cf. Fig. 6): we let iGNN contain 40 ResNet blocks, where the first hidden layer of each block is a Chebnet layer. We fix the learning rate at 1e-4 and train with a batch size of 150. We fix the \(\gamma\) regularization factor of the \(W_2\) loss as 0.375/503.

• Traffic anomaly detection (cf. Fig. 8): we let iGNN contain 40 ResNet blocks, where the first hidden layer of each block is an L3Net layer. We fix the learning rate at 1e-4 and train with a batch size of 200. We fix the \(\gamma\) regularization factor of the \(W_2\) loss as 1.

• Unconditional generation on large graph (cf. Fig. [Cheb95generation] and Fig. [Local95generation]): we let iGNN contain 40 ResNet blocks, each of which contains a ChebNet or an L3Net layer. We fix the learning rate at 5e-4 and train with a batch size of 400. We fix the \(\gamma\) regularization factor of the \(W_2\) loss as 1.

Insensitivity of iGNN to \(\gamma\) for \(W_2\) regularization. We verify that as long as the ResNet is invertible under the choice of \(\gamma\), the value of \(\gamma\) only affects the training efficiency but not the final generative quality. This is because larger \(\gamma\) values restrict more the amount of movement by the ResNet, which thus potentially takes longer training epochs before transporting the distribution \(X|Y\) to \(H|Y\) at each \(Y\). Fig. 17 shows iGNN performance on conditional generation of data on large graph, where we increase the \(\gamma\) factor by 503 times. The performance is nearly identical to that in Fig. 5. On the other hand, Fig. 18 visualizes the generative loss of iGNN on unconditional generation of data on large graph, where we increase the \(\gamma\) factor by 5 times. Comparing to the generative loss in Fig. [Cheb95loss], the loss takes more epochs to converge, but eventual generative performance are nearly identical.

Insensitivity of iGNN to other hyperparameters. We verify that iGNN is insensitive to other hyperparameter choices such as learning rate, batch size, and the number of residual blocks. Figures 18 shows the loss trajectories of iGNN under alternative choices of these hyperparameters, where the setup is identical to that in Fig. [Cheb95generation]. We see that training losses under these choices all converge reasonably fast. We omit comparing the covariance matrices because the generative quality barely differs from those in Fig. [Cheb95generation].

References↩︎