# Invertible Neural Networks for Graph Prediction

Chen Xu111Emails<EMAIL_ADDRESS>H. Milton Stewart School of Industrial and
Systems Engineering, Georgia Institute of Technology. Xiuyuan
<EMAIL_ADDRESS>Department of Mathematics, Duke University Yao
<EMAIL_ADDRESS>H. Milton Stewart School of Industrial and
Systems Engineering, Georgia Institute of Technology.

###### Abstract

In this work, we address conditional generation using deep invertible neural
networks. This is a type of problem where one aims to infer the most probable
inputs $X$ given outcomes $Y$. We call our method invertible graph neural
network (iGNN) due to the primary focus on generating node features on graph
data. A notable feature of our proposed methods is that during network
training, we revise the typically-used loss objective in normalizing flow and
consider Wasserstein-2 regularization to facilitate the training process.
Algorithmic-wise, we adopt an end-to-end training approach since our objective
is to address prediction and generation in the forward and backward processes
at once through a single model. Theoretically, we characterize the conditions
for identifiability of a true mapping, the existence and invertibility of the
mapping, and the expressiveness of iGNN in learning the mapping.
Experimentally, we verify the performance of iGNN on both simulated and real-
data datasets. We demonstrate through extensive numerical experiments that
iGNN shows clear improvement over competing conditional generation benchmarks
on high-dimensional and/or non-convex data.

##  1 Introduction

We are interested in the problem of finding an invertible neural network $g$,
such that given a $Y$ subject to $Y=g(X)$ for an input $X$, we can find the
corresponding $X=g^{-1}(Y)$. In particular, this can be viewed as a type of
inverse problem. Inverse problem can be largely categorized as either one-to-
one (i.e., $Y=g(X)$ with an one-to-one mapping $g$) or can allow for one-to-
many mapping (i.e., $Y$ is categorical and $X\in\mathbb{R}^{d}$, where $d$
denotes the number of features). This work focuses on the latter (i.e., we do
not necessarily need one-to-one invertible mapping), which naturally arises in
many applications. In molecular design (Sánchez-Lengeling & Aspuru-Guzik,
2018), scientists want to infer features of the molecule that lead to certain
outcomes (e.g., type of disease). In anomaly detection (Zong et al., 2018)
such as power outage detection (Al-Shaalan, 2019), studying distributions of
normal vs. anomalous features (e.g., weathers or past power output by sensors)
is important for future prevention. In market analysis (Lu et al., 2015),
estimating patterns of customers based on their behaviors is crucial for the
accurate recommendation. Figure 1 illustrates how our iGNN applies to generic
Euclidean and non-Euclidean (e.g., graph) data, which subsume these
application examples above.

Figure 1: How iGNN applies to non-graph (left) and graph (right) observations.
The dashed lines indicate one-to-many mapping from $Y$ to either $X|Y$ or
$H|Y$. In our approach, $H|Y$ is first generated before mapping through an
invertible NN to $X|Y$. As a practical example, suppose the graph on the right
subfigure is a power grid, where each node $i\in\\{A,\ldots,E\\}$ denotes a
power generator, each $X_{i}\in\mathbb{R}^{C}$ denotes historical power
outputs from generator $i$, and each $Y_{i}\in\\{0,1\\}$ denotes whether
generator $i$ is functioning normally (e.g., $Y_{i}=0$) or abnormally (e.g.,
$Y_{i}=1$). For system monitoring, it is important to both predict (indicated
by solid black arrows) the status of generators given historical observations
and generate (indicated by dashed black arrows) possible circumstances (e.g.,
values on all generators) based on generator status for cause analyses.

This one-to-many task is synonymous with conditional generation, where the
goal is to generate $X|Y$ based on specific outcomes $Y$. Before introducing
common approaches in conditional generation, we acknowledge several
(unconditional) generation approaches, which address the simpler question of
generating samples as close to features $X$ as possible. At present,
generative adversarial networks (GAN) (Goodfellow et al., 2014; Gulrajani et
al., 2017) and variational auto-encoders (VAE) (Kingma & Welling, 2014, 2019)
are two of the most widely studied and popular frameworks. Both approaches
have created remarkable successes in many fields (Makhzani et al., 2015; Zhu
et al., 2017; Ledig et al., 2017). However, they also have clear limitations,
including the inability to provide an exact evaluation of data likelihood for
new samples and difficulties in training, such as mode collapse (Salimans et
al., 2016) and posterior collapse (Lucas et al., 2019).

For our purpose of conditional generation, there have also been several
streams of works, with an abstract description presented on the left of Figure
2. In the following, we first summarize the limitations of these methods and
the open questions in Section 1.1. Then, we introduce the general structure of
the invertible neural network used in our paper. Section 1.2 summarizes our
contributions and presents the roadmap for this work.

$Y$$H|Y$$X|Y$$W_{c}H+b_{c}$$W_{g}e_{Y}+b_{g}$$F_{\Theta}(X)$$F_{\Theta}^{-1}(H)$Generic
NNLinear generator &classifierOur
approach$X|Y$$H_{\Theta}(X|Y)$$G_{\Theta^{\prime}}(Y,Z)$Invertible
NN$\begin{bmatrix}Y\\\ Z\end{bmatrix}$Current Approaches

Figure 2: Abstract comparison of (several) current conditional generative
models (left) and ours (right). In the forward mapping through $H_{\Theta}$,
these methods either concatenates $X,Y$ (Isola et al., 2017; Ardizzone et al.,
2019b) or zero-pads $X$ (Ardizzone et al., 2019a). In the generative procedure
through $G_{\Theta^{\prime}}$, both the outcome variable $Y$ and random noise
$Z$ are taken as network inputs. If the networks are designed to be invertible
(Ardizzone et al., 2019b, a), only a single network is needed whereby
$G_{\Theta^{\prime}}=H^{-1}_{\Theta}$. These approaches inevitably increase
the input and output dimensions of the network, and because $Y$ and $Z$ are
inherently independent, practical training difficulties often arise. In
comparison, our approach maintains the dimensionality of the feature space
when training an invertible NN. During generation, we map one-hot-encoded $Y$
through a simple linear generator onto random vectors $H|Y$ (whose
distributions are disjoint over $Y$) before mapping through the inverse
mapping $F_{\Theta}^{-1}$. By disjointness, $H|Y$ allows a simple prediction
of $Y$.

###  1.1 Background

There have been many studies on conditional versions of GAN (cGAN) (Mirza &
Osindero, 2014; Isola et al., 2017), where the outcome $Y$ and random noise
$Z$ are both taken as inputs to the generator. However, these methods suffer
from similar issues as their unconditional counterparts; we would compare iGNN
with cGAN in our experiments. More recently, the work (Ardizzone et al.,
2019a) directly builds conditional invertible neural networks (CINN) for
analyzing inverse problems. A single invertible NN is trained by minimizing
maximum mean discrepancy (MMD) losses in the latent space $Z$ and the original
data space $X$; we call their method CINN-MMD for comparison. However, this
method also takes in the outcome $Y$ as an additional input to its invertible
NN, thus increasing the dimensionality of the original problem. In addition,
the scalability of MMD-based methods in high-dimensional spaces and choices of
kernel bandwidths can also limit its practical performance, as will be shown
in the experiments.

In this work, we tackle the question of conditional generation for inverse
problems by extending techniques in normalizing flow (Grathwohl et al., 2019;
Behrmann et al., 2019; Ziegler & Rush, 2019; Kobyzev et al., 2020; Ardizzone
et al., 2019b). Flow-based methods estimate arbitrarily complex densities via
the direct maximum likelihood estimation (MLE) approach–they transport
original random features $X$ into more tractable ones (e.g., standard
multivariate Gaussian) through invertible neural networks. However, most
methods only focus on unconditional generation. More recently, (Ardizzone et
al., 2019b) tackles the conditional generative task by combining conditional
invertible neural networks (Ardizzone et al., 2019a) with the normalizing
flow, a method which we call CINN-Nflow for comparison. Nevertheless, CINN-
Nflow suffers from a similar issue as CINN-MMD–by treating $Y$ as an input,
the data dimension is increased so that training is more difficult.

In addition, there have been many architecture designs of the invertible
neural network $F_{\Theta}$ (see Figure 2, right). Some are based on coupling
layers (Dinh et al., 2017) and autoregressive flows (Wehenkel & Louppe, 2019)
while others build on residual networks (Behrmann et al., 2019) and ordinary
differential equations (Grathwohl et al., 2019). At present, we choose
$F_{\Theta}$ to be an invertible residual network (iResnet), due to its
expressiveness and efficiently-generated inverses as proposed in (Behrmann et
al., 2019). Denote $F_{b,\Theta_{b}}:\mathbb{R}^{C}\rightarrow\mathbb{R}^{C}$
as the $b$-th residual block with its own trainable parameters $\Theta_{b}$,
which has the form

$\displaystyle F_{b,\Theta_{b}}(x)$ $\displaystyle:=x+g_{b,\Theta_{b}}(x),$
(1) $\displaystyle g_{b,\Theta_{b}}$
$\displaystyle:=W_{b_{l}}\circ\phi_{b_{l}}\circ\ldots\phi_{b_{1}}\circ
W_{b_{1}},$ (2)

where each pair of $\\{W_{b_{i}},\phi_{b_{i}}\\}$ in (2) denotes the weight
matrix and the activation function at the $i$-th layer of the $b$-th residual
block. Then, $F_{\Theta}:=F_{B,\Theta_{B}}\circ\ldots\circ F_{1,\Theta_{1}}$
is the concatenation of $B$ residual blocks. To ensure invertibility of each
block, (Behrmann et al., 2019) suggested using spectral normalizing so that
the Lipschitz constant of each $g_{b,\Theta_{b}}$ is less than one. However,
spectral normalization is computationally expensive and may not be sufficient
beyond fully-connected layers, so that we propose an alternative
regularization technique in Section 2.3.

Therefore, although conditional generation has been a long-standing and
actively researched area, there are still several open questions.

  1. 1.

How to design (graph) generative models that do not increase difficulties in
training? Treating $Y$ as additional inputs increase the data dimension. For
non-invertible networks such as cGAN and its variants, doing so requires even
more careful hyperparameter tuning and tends to cause many practical issues,
as described earlier. For invertible networks such as CINN-MMD or CINN-Nflow,
doing so causes more scalability issues in higher dimensions due to their
training objectives. In addition, conditional node feature generation on graph
without these training difficulties has not been well-addressed.

  2. 2.

How to reformulate the training objective in a natural way that also allows
for prediction? Most conditional generative models consider the same loss as
their unconditional counterparts and mainly differ by taking in $Y$ as network
inputs. In addition, for prediction, these models either require building
separate supervised networks or require $Y$ as an additional output by the
generative model (e.g., CINN-MMD). Both of these approaches increase
complexities in training.

###  1.2 Contribution

Compared to current approaches (see Figure 2, left), our approach (see Figure
2, right) circumvents the inevitable dimensionality issue: we first generate
hidden variables $H|Y$ that have disjoint supports over $H$ and then trains a
single invertible network that maintains the original data dimension. Due to
disjointness, the prediction of $Y$ can also be performed by simple linear
classifiers. The method is called invertible graph neural network (iGNN) due
to our primary focus on graph data. In summary, our technique contributions in
response to the open questions above include

  1. 1.

For the training objective, we formulate a conditional generative loss based
on the change-of-variable formula (see Eq. (3)), which naturally extends to
graph observations (see Eq. (7)). Predicting the outcome $Y$ is simple based
on our problem design and is built as a part of our end-to-end training. In
particular, the invertible neural network does not treat the conditional
variable $Y$ as an additional network input to maintain the original data
dimension.

  2. 2.

In network design, we impose the Wasserstein-2 regularization on each residual
block of the invertible residual network in iGNN. In contrast to the much more
expensive spectral normalization (Behrmann et al., 2019), this regularization
is easy-to-implement and facilitates smoother density transport.

  3. 3.

Theoretically, we establish the identifiability of a true invertible mapping
and verify its existence and invertibility based on the well-posedness of the
initial-value problem for ordinary differential equations. We also show that
graph neural networks with enough (resp. insufficient) expressiveness can
(resp. cannot) theoretically learn the true invertible mapping.

  4. 4.

In experiments, we compare the performance of our model with competing models
on simulated and real data, with clear improvement over competing models on
high-dimensional and/or non-convex graph data. We also verify the importance
of model expressiveness through specific simulation designs.

The rest of the paper is organized as follows. Section 2 contains our method.
We describe how the training objective is formulated to include prediction and
conditional generation, how the framework is extended to graph data, and how
each residual block is regularized through the Wasserstein-2 penalty. Section
3 contains theoretical analyses with three parts. We first establish the
identifiability of the true invertible mapping, assuming its existence. Then,
we analyze when such invertible mappings exist upon connecting with ODE.
Lastly, we derive a closed-form expression for the invertible mapping in a
special case and consider the expressiveness of our iGNN in two examples.
Section 4 contains simulated and real-data experiments, which show the
improved performance of our method over other conditional generative models,
especially when data are high-dimensional and/or non-convex. In particular, we
also verify the model expressiveness statements in Section 3 with numerical
results. Section 5 concludes our work and discusses future directions. Section
6 contains complete derivation and proofs. Appendix A contains details
regarding experiment setups and additional numerical results.

##  2 Method

We first describe how the training objective is formulated in Section 2.1 for
generic Euclidean data. We then discuss the extension to non-Euclidean data,
such as graph observations in Section 2.2. Lastly, we propose the
regularization technique of residual blocks described in (2) in Section 2.3.

###  2.1 End-to-end training for conditional generation and prediction

We now present the concrete mathematical formulation of our training objective
based on the abstract representation on the right of Figure 2. Compared to
other conditional generative models, our iGNN performs a simple but effective
modification of the unconditional generative loss based on the change-of-
variable formula. As we will see, doing so has several benefits: (1) maintain
the original data dimension of $X$ instead of increasing it to include $Y$,
(2) facilitate model training as the network only needs to transport each
$X|Y$ to its corresponding $H|Y$ instead of all $X|Y$ to a standard Gaussian
variable $Z$, and (3) enable easy classification, as each $H|Y$ is designed to
be disjoint.

We first consider the case where the random feature $X\in\mathbb{R}^{d}$.
Given an invertible network $F_{\Theta}$ and additional categorical outcome
variables $Y\in[K]$ for $K$ classes, we consider the following conditional
generative loss:

$\displaystyle\log p_{X|Y}(X)$ $\displaystyle=\log
p_{H|Y}(F_{\Theta}(X))+\log|\det J_{F_{\Theta}}(X)|,$ (3) $\displaystyle
L_{g}$ $\displaystyle:=-\log p_{X|Y}(X),$

where $L_{g}$ denotes the generative loss to be minimized and each $H|Y$ is
disjoint over the value of $Y$. Note that conditional generation based on (3)
is simple—given different outcomes $Y$, we can sample from different
distributions $H|Y$ and then perform the inverse mapping $F_{\Theta}^{-1}$. In
implementation, we let distributions $H|Y$ be isotropic multivariate Gaussian
distributions, with mean parameters $\mu(H|Y):=W_{g}e_{Y}+b_{g}$ using
trainable parameters $\\{W_{g},b_{g}\\}$ and variances to be small enough for
disjointness. The term $e_{Y}$ denotes the $Y$-th standard basis vector in
$\mathbb{R}^{K}$, and the subscript $g$ denotes generation.

While the loss in (3) may seem simple, we highlight its benefit over the
current CINN-Nflow approach (Ardizzone et al., 2019b), whose loss is also
based on the change-of-variable formula. CINN-Nflow treats $Y$ as the network
input and writes

$\displaystyle\log p_{X|Y}(X)$ $\displaystyle:=\log
p_{Z}(G_{\Theta^{\prime}}(X,Y)|_{Z})+\log|\det J_{G_{\Theta^{\prime}}}(X,Y)|,$
(4) $\displaystyle G_{\Theta^{\prime}}(X,Y)$
$\displaystyle:=[G_{\Theta^{\prime}}(X,Y)|_{Y},G_{\Theta^{\prime}}(X,Y)|_{Z}].\quad\text{(Dimension
increase)}$

In (4), $Z$ denotes a standard multivariate Gaussian random vector and
$G_{\Theta}$ is the conditional invertible neural network. In contrast to our
objective (3), the output $G_{\Theta^{\prime}}(X,Y)$ in (4) contains
predictions for both $Y$ and $Z$, so that the dimension of the original
feature space is increased. In addition, because $Z$ and $Y$ are inherently
independent, the conditional generative model $G_{\Theta^{\prime}}^{-1}$ must
transport the uni-modal random vector $Z$ into multiple $X|Y$, each of which
may be multi-modal. Both procedures increase the difficulty in training.

In addition, our conditional formulation in (3) also allows one to perform
classification (e.g., estimate $Y|X$) using simple classifiers. For instance,
assuming $Y$ is binary but the original decision boundary of $Y|X$ is highly
nonlinear. Instead of relying on complex classifiers (e.g., another supervised
neural network), one may build a simple one (e.g., linear classifier) as long
as (1) distributions $H|Y$ do not overlap (e.g., disjoint multivariate
Gaussian distributions) (2) the invertible mapping can accurately transport
the distribution of $X|Y$ to match that of $H|Y$. Mathematically, we thus
minimize the classification loss using a simple linear classifier on
$F_{\Theta}(X)$:

$\mathcal{L}_{c}:=\mathcal{L}(Y,W_{c}F_{\Theta}(X)+b_{c})=-e_{Y}^{T}\text{softmax}(W_{c}F_{\Theta}(X)+b_{c}).$
(5)

Therefore, our end-to-end training objective can be written as

$\underset{\\{\Theta,\\{W_{c},b_{c}\\},\\{W_{g},b_{g}\\}\\}}{\min}\mathbb{E}_{X,Y}[\mathcal{L}_{g}+W_{2}(F_{\Theta})+\mu\mathcal{L}_{c}],$
(6)

which is estimated using training data. The term $W_{2}(F_{\Theta})$ will be
introduced in Section 2.3, Eq. (8) to ensure model invertibility and
smoothness in density transport. The term $\mu$ is a parameter that controls
the relative rate at which the generative model trains with respect to the
linear classifier. In practice, $\mu$ can be chosen small because of problem
formulation—classification is easy as long as estimated distributions using
$F_{\Theta}(X)|Y$ accurately match those of $H|Y$, which are disjoint. Figure
3 illustrates the result on a simple simulated two-moon data. For prediction,
both parts of the two-moon data are mapped separately through the invertible
map to disjoint parts of the Gaussian mixture, which can be easily classified
via the linear classifier $x=y$. For conditional generation, one simply
samples from the Gaussian mixture, and the invertible map $F_{\Theta}^{-1}$
generates the desired $X|Y$.

Figure 3: Illustration of iGNN on a toy example. We compare $X|Y$ vs.
$\hat{X|Y}:=F^{-1}_{\Theta}(H|Y)$ and $H|Y$ vs. $\hat{H}|Y:=F_{\Theta}(X|Y)$,
where the estimates are very close to the true conditional distributions.

###  2.2 Generative loss on graph data

A key emphasis of our paper is on non-Euclidean observation, such as graph
data, where we currently focus on node classification and node feature
generation. Notation-wise, given a graph $G=(V,E)$ with $|V|=V$ (with a slight
abuse of notation), let $A$ be its adjacency matrix, $D_{A}$ be the degree
matrix, $X^{(v)}\in\mathbb{R}^{C},v=1,\ldots,V$ be the node feature, and
$Y^{(v)}\in[K]$ be the node label. Denote $X\in\mathbb{R}^{V\times C}$ as the
matrix of node features and $Y\in\\{1,\ldots,K\\}^{V}$ as node labels.

For training, we still minimize the end-to-end objective in (6). However,
there is a key difference regarding the evaluation of the log-likelihood of
input node features $X$:

$\displaystyle\log p_{X|Y}(X)$ $\displaystyle=$ $\displaystyle\log
p_{H|Y}(F_{\Theta}(X))+\log\det J_{F_{\Theta}}(X)$ (Change of variable)
$\displaystyle\overset{(i)}{=}$ $\displaystyle\sum_{v=1}^{V}\left[\log
p_{H^{(v)}|Y^{(v)}}(F_{\Theta}(X)^{(v)})\right]+\log\det J_{F_{\Theta}}(X)$
(Independence) $\displaystyle\overset{(ii)}{=}$
$\displaystyle\sum_{v=1}^{V}\left[\log
p_{H^{(1)}|Y^{(1)}}(F_{\Theta}(X)^{(v)})\right]+\log\det J_{F_{\Theta}}(X),$
(Homogeneity) (7)

where in (i), we assume independence of node features upon transporting the
raw feature, and in (ii), we further assume that the conditional distribution
$H^{(v)}|Y^{(v)}$ at each node is homogenous, so the distribution only depends
on the label of $Y^{(v)}$ rather than the node $v$. In practice, requiring
(ii) is important as the number of nodes is often significantly larger than
the number of node classes. In that case, it is nearly infeasible to pre-
specify node-dependent disjoint conditional distributions $H^{(v)}|Y^{(v)}$.
Computationally, we would perform vectorize-then-unvectorize operations of
inputs $X$ and outputs $F_{\Theta}(X)$ to make sure
$J_{F_{\Theta}}(X)\in\mathbb{R}^{V\cdot C\times V\cdot C}$.

###  2.3 Invertible flow Resnet with Wasserstein-2 regularization

Recall that we choose $F_{\Theta}$ in our iGNN as an invertible residual
network (see right of Figure 2), which is a concatenation of residual blocks
defined in Eq. (2). To ensure network invertibility, instead of performing
spectral normalization of the weight matrices in each residual block, we
consider the following regularization on the movement of each residual block:

$W_{2}(F_{\Theta}):=\frac{1}{2}\sum_{b=1}^{B}\|g_{b,\Theta_{b}}(X)\|_{2}^{2}.$
(8)

Note that (8) is the Wasserstein-2 regularization (Bolley et al., 2012).
Figure 12 in the appendix shows that our model maintains invertibility on both
in-distribution and out-of-distribution data using the Wasserstein-2
regularization in (8).

There are several benefits of the Wasserstein-2 regularization compared to the
spectral normalization technique (Behrmann et al., 2019). Computationally, the
Wasserstein-2 penalty can be easily implemented as the network performs
forward passes on data. In contrast, spectral normalization relies on the
power iteration at each training epoch to normalize the spectral norm of
weight matrices in all residual blocks. Doing so can be very computationally
expensive when both the total number of residual blocks and the number of
hidden neurons are large, whereby the procedure must spectrally normalize many
high-dimensional weight matrices. In addition, the Wasserstein-2 penalty is
agnostic of the residual block design—each block can consist of only fully-
connected (FC) layers or generic graph filters (Defferrard et al., 2016; Cheng
et al., 2021). In contrast, spectral normalization is applicable when each
residual block concatenates fully-connected layers (FC) with contractive
nonlinearities (e.g., ReLU, ELU, Tanh). When each residual block contains more
than FC layers (e.g., general graph filters), it can be insufficient to
require the spectral norm of all weight matrix parameters to be less than one.

##  3 Theory

We provide several theoretical results in this section. First, we show in
Section 3.1 that under idealized assumptions, the optimal solution of the
conditional generative loss (3) identifies the true invertible mapping that
generates $X|Y$ for each $H|Y$. Then, we show in Section 3.2 that the true
invertible mapping exists under the well-posedness of an ordinary differential
equation, which can be verified under certain Lipschitz conditions. Lastly, we
discuss in Section 3.3 the expressiveness of our iGNN for Gaussian processes
in two examples that depend on graph topology. Throughout this section, we
assume that $X_{t}\in\mathbb{R}^{d}$, where graph features
$X_{t}\in\mathbb{R}^{V\times C}$ can be equivalently interpreted as random
vectors in $\mathbb{R}^{V\cdot C}$, with correlation among features according
to the graph structure. All derivatives and proofs appear in Appendix 6.

###  3.1 Identifiability

Upon minimizing $L_{g}$ in (3), we first have the following proposition
regarding the identifiability of a true invertible mapping $F^{*}$.

######  Proposition 3.1 (Identifiability of the true mapping).

Let $Y\in\\{1,\ldots,K\\}$ be discrete labels. Assume that

  * •

For any $i,j\in[K]$ with $i\neq j$, the random vectors $H|Y=i$ and $H|Y=j$
follow disjoint tractable distributions (e.g., multivariate Gaussian).

  * •

There exists a bijective mapping
$F^{*}:\mathbb{R}^{d}\rightarrow\mathbb{R}^{d}$ such that for every $i\in[K]$,
$\hat{H}_{i}:=F^{*}(X|Y=i)$ follows the same distribution as $H|Y=i$.

  * •

There exist parameters $\Theta^{*}$ such that $F_{\Theta^{*}}=F^{*}$, where
$F_{\Theta}:\mathbb{R}^{d}\rightarrow\mathbb{R}^{d}$ is the pre-speficied
class of invertible models that depends on parameters $\Theta$.

Based on (3), denote

$\displaystyle\hat{\Theta}_{N}$
$\displaystyle:=\arg\max_{\Theta}\sum_{n=1}^{N}\log
p_{H|Y}(F_{\Theta}(X_{n}))+\log|\det J_{F_{\Theta}}(X_{n})|,$ (9)
$\displaystyle N_{K}$
$\displaystyle:=\min_{i=1,\ldots,K}\sum_{n=1}^{N}\textbf{1}(Y_{n}=i),$

where (9) is the maximum-likelihood estimator based on $N$ independent and
identically distributed (upon conditioning on values of $Y$) samples
$\\{(X_{n},Y_{n})\\}_{n=1}^{N}$. Then, under regularity conditions for maximum
likelihood estimation, we have that as $N_{K}\rightarrow\infty$,
$\hat{\Theta}_{N}\rightarrow\Theta^{*}$ with probability 1.

###  3.2 Existence and invertibility

Note that Proposition 3.1 relies on the idealized assumption that there exists
a true invertible mapping $F^{*}$ that transports $X|Y$ to $H|Y$ as desired.
In this subsection, we utilize the Fokker-Planck equation of an Ornstein-
Uhlenbeck (OU) process to show that such an invertible mapping exists and can
be approximated by a residual network with finite model complexity.

We first reduce the problem to the construction of invertible flow for the
unconditional case. The conditional case can be viewed as having a joint flow
which coincides with each individual flow upon restricting to the support of
$p_{X|Y}$. More precisely, this joint flow maps each $p_{X|Y}$ to its
corresponding $p_{H|Y}$, all of which are disjoint over $Y$; disjointness can
be achieved by re-scaling and shifting of a standard Gaussian (see explanation
earlier below Eq. (3)). As a result, once the existence of invertible flow
from some $p_{X}$ to a normal random vector can be constructed, then composed
with the shifting and re-scaling operations which are invertible transports,
one can construct an invertible joint flow for the conditional generating of
$X|Y$.

Below, we consider the unconditional case where the desired flow transports
from data distribution $p_{X}$ in $\mathbb{R}^{d}$ to a normal distribution
$p_{H}$ which is $\mathcal{N}(0,I_{d})$. Recall that the continuous-time flow
is represented by an initial value problem (IVP) of ODE

$\frac{dx_{t}}{dt}={\bf f}(x_{t},t),$ (10)

with initial value $x_{0}\sim p_{X}$. The solution mapping
$x_{t}:=P_{t}(x_{0})$ for $t>0$ is invertible as long as (10) is well-posed.
There may be more than one ${\bf f}(x,t)$ which induces such a flow towards
$p_{H}$, and to show existence, we construct ${\bf f}$ to correspond to the
following Fokker-Planck equation of an OU process ($\nabla$ denotes
$\nabla_{x}$)

$\partial_{t}\rho=\nabla\cdot(\rho\nabla V+\nabla\rho),\quad
V(x)=\frac{|x|^{2}}{2},\quad\rho(x,0)=p_{X}(x),$ (11)

where $\rho(x,t)$ represents the probability density of the OU process at time
$t$ starting from $\rho(x,0)$. Comparing (11) to the Liouville equation of
(10) which is $\partial_{t}\rho=-\nabla\cdot(\rho{\bf f})$, where $\rho(x,t)$
represents the distribution of $x_{t}$, we see that the density transportation
under the flow can be made the same as in (11) if we set the force ${\bf f}$
as

$-{\bf f}(x,t)=\nabla V(x)+\nabla\log\rho(x,t)=x+\nabla\log\rho(x,t).$ (12)

We first show the differentiability of ${\bf f}(x,t)$ in the following lemma.

######  Lemma 3.2.

Suppose $\rho(x,t)$ is the solution to (11) from $\rho(x,0)=p_{X}$, then ${\bf
f}(x,t)$ as in (12) is smooth on $\mathbb{R}^{d}\times(0,\infty)$.

The following proposition directly gives the well-posedness of the IVP, and
subsequently the invertibility of the continuous-time flow, assuming an
$x$-Lipschitz constant of ${\bf f}(x,t)$. The $x$-derivative of a vector field
${\bf v}(x)$ is denoted as $\nabla\otimes{\bf v}(x)$ which is a $d$-by-$d$
matrix and we use $\|\cdot\|$ to stand for matrix operator norm.

######  Proposition 3.3.

Suppose ${\bf f}(x,t)$ as in (12) with $\rho(x,t)$ being the solution to (11)
satisfies that

$\|\nabla\otimes{\bf f}(x,t)\|\leq L,\quad\forall t>0,\quad
x\in\mathbb{R}^{d},$

then the IVP (10) is well-posed.

The boundedness of the $x$-Lipschitz constant of ${\bf f}(x,t)$ can be shown
when the initial distribution $\rho(x,0)$ is sufficiently regular and decay
fast on $\mathbb{R}^{d}$, e.g., smooth and compactly supported. The precise
analysis is deferred to future work. Under generic conditions, the solution
$\rho(x,t)$ of Fokker-Planck equation converges to the equilibrium
distribution $\rho_{\infty}\propto e^{-V}$ exponentially fast under e.g.
$L^{2}$ or Wasserstein-2 distance (Bolley et al., 2012). This means that at
time $T\sim\log\frac{1}{\epsilon}$, the transported distribution of $x_{T}$ is
$\epsilon$-close to the normal density. One can construct a discrete-time flow
on $[0,T]$ to approximate the continuous-time one, and when the force ${\bf
f}(x,t)$ is sufficiently regular, at each $t$ it can be approximated by a
residual block with finite many trainable parameters. Furthermore, since the
$x$-Lipschitz constant of ${\bf f}(x,t)$ is bounded by $L$, invertibility of
the ResNet can be fulfilled as long as the discrete time-step $\Delta t<1/L$,
and the number of residual blocks equals $T/\Delta t$ which can be made $\sim
TL$.

###  3.3 Model expressiveness

Above in Section 3.2 that we verified the existence of an invertible mapping
and argued for approximation of the force ${\bf f}(x,t)$ in (12) using a
residual network with finite model complexity. In this subsection, we further
establish the exact form of ${\bf f}(x,t)$ for Gaussian processes (GP). In
particular, we discuss two graph cases where the covariance matrix $\Sigma$ of
the initial distribution of the GP directly relates to properties of the
graph, with an interesting implication regarding model expressivess of GNN.

We first provide the following analytic form of ${\bf f}(x,t)$ in (12), which
applies to general vector data in $\mathbb{R}^{d}$:

######  Lemma 3.4.

If $X_{0}\sim N(0,\Sigma)$ and $X_{\infty}\sim N(0,I_{d})$, then

${\bf f}(x,t)=-(I_{d}-\Sigma_{t}^{-1})x,$ (13)

where $\Sigma_{t}:=(1-\exp(-2t))I_{d}+\exp(-2t)\Sigma$ is the covariance
matrix that converges to $I_{d}$ exponentially fast.

Because the force ${\bf f}(x,t)$ in (13) is linear, it can always be
approximated by residual blocks of fully-connected layers, which are universal
approximators (Haykin, 1998; Lu et al., 2017).

On the other hand, one can interpret the random vector $X_{0}$ as the
collection of one-dimensional node features on a graph with $d$ nodes, with
$\Sigma$ encoding feature dependency. In such cases, it is not always true
that any graph filters can approximate $\Sigma$ and as a consequence, the
covariance $\Sigma_{t}^{-1}$ in (13) that depends on $\Sigma$; one example
that spectral-based GNN fails to do so will be given in Example 1. In
particular, we have the following propositions on the properties of
$\Sigma_{t}$ under different assumptions of $\Sigma$, which help us identify
networks that are able or unable to approximate the force ${\bf f}(x,t)$. More
precisely, we relate $\Sigma$ to the spectral-based GNN in one case and
spatial-based GNN in another.

For notation consistency with previous discussion on graph data, we change the
dimensionality $d$ to $V$ in statements on graph data from now on. For
simplicity, all analyses onward implicitly assume node features have dimension
$C=1$.

Case I. Spectral-based presentation of $\Sigma$: We first show that when
$\Sigma$ can be written as a polynomial expression of the graph Laplacian $L$,
then $\Sigma_{t}$ has a closed-form in terms of the polynomial evaluated on
eigenvalues of the graph Laplacian, whereby $\Sigma_{t}^{-1}$ can be
approximated up to arbitrary precision by another polynomial. Hence, spectral-
based GNN can in theory correctly approximate the true force in (13).

######  Proposition 3.5 (Spectral-based $\Sigma$).

Denote $L\in\mathbb{R}^{V\times V}$ as the graph Laplacian. Suppose there
exists a polynomial $P$ such that $\Sigma=P(L)$. Then

$\Sigma_{t}^{-1}=\sum_{i=1}^{V}((1-\exp(-2t))+\exp(-2t)P(\lambda_{i}(L)))^{-1}U_{i}U_{i}^{T},$

where $\lambda_{i}(L)$ denotes the $i$-th largest eigenvalue of $L$ and
$U=[U_{1},\ldots,U_{V}]$ with $L=U\Lambda U^{T}$ under eigen-decomposition.

The implication of Proposition 3.5 is that under certain regularity condition
on $P$, there exists another polynomial $P_{t}$ that can approximate
$\Sigma_{t}^{-1}$ up to arbitrary precision. For instance, when $P$ denotes
the Chebyshev polynomial of order $k$, then this approximation can be achieved
by (possibly) raising the degree of polynomial above $k$.

We also consider a corollary below where the covariance matrix $\Sigma$ is not
accurately represented by $P(L)$, so that estimation errors exist.

######  Corollary 3.6.

Under the same notations in Proposition 3.5, suppose

(1) $\Sigma$ and $L$ share the same eigen-bases.

(2) For some polynomial $\tilde{P},\min_{i}\tilde{P}(\lambda_{i}(L))>0$, where
$\lambda_{i}(A)$ denotes the $i$-th largest eigenvalue of a matrix.

(3) There exists $\delta>0$ such that
$\max_{i}|\lambda_{i}(\Sigma)-\tilde{P}(\lambda_{i}(L))|\leq\delta$.

Define $\tilde{\Sigma}_{t}:=(1-\exp(-2t))I_{V}+\exp(-2t)\tilde{P}(L)$. Then

$\rho(\Sigma_{t}^{-1}-\tilde{\Sigma}_{t}^{-1})=O(e^{-2t}\delta),$

where $\rho$ denotes the spectral norm of the matrix difference, and the
constant in big-O depends on the minimum eigenvalues of $\Sigma$ and
$\tilde{P}(L)$. In particular, note that $\tilde{\Sigma}_{t}$ can be
approximated up to arbitrary precision based on the implication of Proposition
3.5.

Case II. Spatial-based presentation of $\Sigma$: We next show that under
certain spatial properties of $\Sigma$, $\Sigma_{t}^{-1}$ satisfies similar
spatial properties and can thus be approximated by spatial-based GNN filters
such as the low-rank locally learnable network (L3net) (Cheng et al., 2021).
However, it may not be approximated by spectral-based GNN such as Chebnet
(Defferrard et al., 2016) as shown in Example 1. We first define the desired
spatial property in terms of locality.

######  Definition 3.7 ($v$-locality).

Denote $A\in\mathbb{R}^{V\times V}$ as the graph adjacency matrix for $d$
nodes. We define that the covariance matrix $\Sigma\in\mathbb{R}^{V\times V}$,
which encodes correlation among node features, is $v$-local for some integer
$v\geq 1$ if

$\Sigma_{ij}=0\text{ if }j\notin\mathcal{N}^{(v)}_{i},$

where $\mathcal{N}^{(v)}_{i}$ denotes the set of $v$-th neighbors of node $i$
(e.g., all nodes accessible from node $i$ within $v$ steps along adjacent
nodes).

######  Proposition 3.8 (Local $\Sigma$ and $\Sigma^{-1}$).

Suppose $\Sigma$ and $\Sigma^{-1}$ are $v_{1}$-local and $v_{2}$-local based
on Definition 3.7. Then depending on values of $t$, $\Sigma_{t}^{-1}$ in (13)
can be expressed by power series of $\Sigma$ or $\Sigma^{-1}$. If the power
series are truncated at order $k$, $\Sigma_{t}^{-1}$ can be approximated by a
$\max(kv_{1},kv_{2})$-local covariance matrix.

On the other hand, the following simple example show that spectral-based GNN
can never learn $\Sigma$ and therefore, the covariance $\Sigma_{t}^{-1}$
required by the force ${\bf f}(x,t)$ in (13). We will revisit this example in
Section 4.3.

######  Example 1 ($\Sigma_{t}^{-1}$ cannot be learned by spectral-based graph
filters).

Consider a simple graph with three nodes $\\{1,2,3\\}$ and two edges
$\\{(1,2),(2,3)\\}$ between nodes. Self-loops at each node are also inserted.
Note that the adjacency matrix $A$ satisfies the property $\pi A\pi^{T}=A$,
where $\pi\in S_{3}$ is any permutation matrix over graph nodes. Thus, any
graph filter $f[A]$ based on $A$ satisfies $\pi f(A)\pi^{T}=f(A)$. However, if
the covariance matrix $\Sigma$ and permutation matrix $\pi$ take the form

$\Sigma:=\begin{bmatrix}1&\rho&0\\\ \rho&0&\rho_{1}\\\
0&\rho_{1}&1\end{bmatrix},\pi=\begin{bmatrix}0&0&1\\\ 0&1&0\\\
1&0&0\end{bmatrix},$

then it is easy to see that $\pi\Sigma\pi^{T}\neq\Sigma$ if
$\rho\neq\rho_{1}$. Therefore, any $f(A)$ cannot approximate $\Sigma$ up to
arbitrary accuracy.

In summary, Propositions 3.5 and 3.8 state the existence of approximations for
$\Sigma_{t}$ under correct model specification. In experiments, we will
explore the importance of architecture design—whether $\Sigma_{t}^{-1}$
can/cannot indeed be approximated by our iGNN under different designs of
residual blocks. Results are shown in Section 4.1 Figures 10 (for Proposition
3.5) and 11 (for Proposition 3.8).

##  4 Numerical Examples

The experiments are organized as follows. In Section 4.1, we first compare our
iGNN with competing methods on simulated examples, where one of the
competitors performs similarly as iGNN. In Section 4.2, we then consider more
challenging simulated data and real-data examples. In particular, the
competitors all perform worse than iGNN in terms of conditionally generating
the high-dimensional and/or non-convex node features. In Section 4.3, we
lastly verify model expressiveness statements in Section 3.3. Code to produce
our results can be found at https://github.com/hamrel-cxu/Invertible-Graph-
Neural-Network-iGNN.

We consider three competing conditional generative models and several two-
sample testing metrics for quantitative assessment. Regarding competing
methods, we compare iGNN with the conditional generative adversarial network
(cGAN) (Isola et al., 2017), the conditional invertible neural network trained
with maximum mean discrepancy (CINN-MMD) (Ardizzone et al., 2019a), and the
same CINN trained with the normalizing flow (CINN-Nflow) (Ardizzone et al.,
2019b). Details and hyperparameters of the competitors in each setting are
described in Appendix A.1. Regarding quantitative assessments, two-sample
testing statistics are used to compare the similarity between two
distributions (e.g., $X|Y$ vs. $\hat{X}|Y$ at different $Y$). In particular,
we use kernel MMD under the Gaussian kernel with various bandwidth choices
(Gretton et al., 2012) and the energy-based test (Székely & Rizzo, 2013),
which utilizes the energy distance or E-statistics for comparisons of
distributions. We remark that since CINN-MMD is based on the MMD test
(although the kernel functions are different), the kernel MMD statistics can
be biased towards this method.

The following settings are fixed in all experiments; other hyperparameters
will be explained separately in each experiment. All graphs are assumed
undirected and unweighted, with inserted self-loops. The activation function
is chosen as ELU (Clevert et al., 2015) due to its continuous derivatives and
the hidden dimension is fixed at 64. During training, we use the Adam
optimizer (Kingma & Ba, 2015) to minimize the end-to-end objective (6). We
terminate the training of all methods once the percentage decrease of training
loss is smaller than 0.01%; for compactness and clarity in visualization, we
show the training and test losses for competing methods in each experiment in
Appendix A.2, Figure 14.

###  4.1 Simple examples

We consider two simulated examples in this section. The first considers non-
graph Euclidean vector data, and the second deals with graph node features. In
particular, we find that CINN-MMD yields competitive performances as iGNN, but
the other two other competing methods perform much worse.

1\. Simulation: 8 Gaussians. We consider one non-graph toy example that is
widely used by conditional generative methods for assessment. Each conditional
distribution $X|Y$ comprises two disjoint Gaussian distributions for four
possible choices of $Y$. The task is to estimate $X|Y$ as closely as possible
for each $Y$. We build the iGNN to have 40 residual blocks, where each block
has four FC layers with 64 hidden nodes. During training, we generate 2500
i.i.d. training and test samples for each of the four labels of $Y$, fix the
learning rate at 5e-4, and train with a batch size of 1000. Figure 4 compares
our generative results with those of CINN-MMD, where both methods can
reasonably and correctly estimate each $X|Y$. Table 1 provides the two-sample
testing statistics, where the values by CINN-Nflow and cGAN are much larger
than those by ours and by CIIN: MMD. These larger values indicate a greater
discrepancy between the true and estimated conditional distributions.

a iGNN loss

b iGNN $\hat{X}|Y$

c CINN loss

d CINN $\hat{X}|Y$

Figure 4: 8 Gaussian visual comparison, iGNN vs. CINN-MMD. Both methods can
estimate the true conditional densities $X|Y$ reasonably well. In addition, we
note that while CINN losses increase after approximately 120 epochs, we did
not pause the training because generative performances before the loss
increase are visually poorer than those after the entire 200 epochs.

Table 1: 8 Gaussian two-sample testing statistics. We considered both the MMD
statistics under the Gaussian kernel with various bandwidths (Gretton et al.,
2012) and the Energy-based statistics (Székely & Rizzo, 2013). Statistics are
weighted by the number of observations at each potential value of $Y$. Both
CINN-Nflow and cGAN yield much larger test statistics than iGNN and CINN-MMD,
indicating a worse estimation of the true conditional distribution.

| MMD: $\alpha=0.1$ | MMD: $\alpha=1.0$ | MMD: $\alpha=5.0$ | MMD: $\alpha=10.0$ | Energy
---|---|---|---|---|---
iGNN | 0.001 | 0.002 | 0.004 | 0.003 | 0.004
CINN-MMD | 0.000 | 0.002 | 0.004 | 0.004 | 0.004
CINN-Nflow | 0.367 | 0.947 | 0.620 | 0.460 | 1.582
cGAN | 0.148 | 0.299 | 0.273 | 0.266 | 0.601

2\. Simulation: 2D convex node features. We now consider the graph introduced
in Example 1. Each node has a binary label so that $Y\in\\{0,1\\}^{3}$ is a
three-dimensional binary vector. Conditioning on a specific binary vector $Y$
out of the eight choices, the random matrix $X|Y$ is defined as

$X|Y:=P_{A}\left(H|Y+\begin{bmatrix}-4&0&4\\\ 0&0&0\end{bmatrix}\right)\text{
s.t. }H^{(v)}|Y^{(v)}\sim\begin{cases}N((1.5,0)^{T},0.1I_{3})&\text{if
}Y^{(v)}=0\\\ N((-1.5,0)^{T},0.1I_{3})&\text{if }Y^{(v)}=1\end{cases},$

where $P_{A}:=D_{A}^{-1}A$ is the graph average matrix. In particular, we
purposefully design the example so that the marginal distribution of each node
feature does not yield enough information for classification (i.e. $Y^{(v)}$
can be either 1 or zero for the same $X^{(v)}$, when features
$X^{(v^{\prime})},v^{\prime}\neq v$ remain unobserved). We build the iGNN to
have 40 residual blocks, where each block has 3 layers in the form of
$\text{FC}\circ\text{FC}\circ\text{L3net}$ (same neighborhood orders as
section 4.1.2). During training, we generate 4000 i.i.d. training and test
samples for each of the eight binary vectors $Y$, fix the learning rate at
5$\times 10^{-4}$, and train with a batch size of 100. Figure 5 shows that the
generate and classification performances by our method are both satisfactory,
and CINN-MMD can yield similar performances. Meanwhile, Table 2 shows the
similar issue that neither CINN-Nflow nor cGAN can yield satisfactory
performances.

######  Remark 4.1 (Using FC layers on graph data).

Note that we can always build each residual block using only FC layers and
ignore the relationship among graph nodes. The main difference from using a
graph filter is that one has to first vectorize the node feature matrix
$X\in\mathbb{R}^{V\times C}$ as
$\widetilde{X}:=\text{vec}(X)\in\mathbb{R}^{V\cdot C}$ and then transport the
high-dimensional density. Doing so can significantly increase the number of
parameters in the model if $V\cdot C$ is large. Generative performance may
also deteriorate in some instances since graph information is ignored.

a iGNN losses

b iGNN $\hat{X}|Y$

c CINN $\hat{X}|Y$

Figure 5: 2D convex node feature visual comparison, iGNN vs. CINN-MMD. We
connect in (b) and (c) two-dimensional node features belonging to the same
3-by-2 feature matrix $X$ by blue lines to reveal the high-dimensional
distribution of $X$. The conditional generative performances are satisfactory
at different values of the random vector $Y\in\\{0,1\\}^{3}$. In addition,
generative or classification losses and classification errors converge to
desired values in a few epochs.

Table 2: 2D convex node feature two-sample testing statistics. Statistics are
weighted by the number of observations at each potential value of $Y$. Similar
to before, neither CINN-Nflow nor cGAN can yield satisfactory performances.
Although iGNN yields larger test statistics than that by CINN-MMD, the
generative performance in Figure 5 shows that both methods yield similar
generative performances.

3 node, convex | MMD: $\alpha=0.1$ | MMD: $\alpha=1.0$ | MMD: $\alpha=5.0$ | MMD: $\alpha=10.0$ | Energy
---|---|---|---|---|---
iGNN | 0.007 | 0.053 | 0.049 | 0.027 | 0.062
CINN-MMD | 0.003 | 0.021 | 0.032 | 0.024 | 0.028
CINN-Nflow | 0.749 | 0.628 | 0.205 | 0.087 | 5.570
cGAN | 0.499 | 1.325 | 1.181 | 1.084 | 2.813

###  4.2 Complex examples

We now consider more challenging node feature generation examples beyond ones
in Section 4.1. We first build upon the three-node graph by considering non-
convex feature distributions $X|Y$. We then consider two high-dimensional
real-data graph examples in anomaly detection, where it is essential to
generate graph features jointly based on the label at each node (recall Figure
1). In all these cases, we will see that CINN-MMD performs worse than iGNN
when node features are non-convex or high-dimensional and that neither CINN-
Nflow nor cGAN yields satisfactory performances as iGNN.

1\. Simulation: 2D non-convex node features. We still use the graph introduced
in Example 1 under a more complex node feature setting. At each node $v$,
features $X^{(v)}$ can be either convex (e.g., linear transform of Gaussian
vectors) or non-convex (e.g., a part of the two-moon dataset as in Figure 3).
The response vector $Y$ is still binary. We ignore details regarding the
generation of $X|Y$, which is very similar to before, except for the way
$H^{(v)}|Y^{(v)}$ is chosen. We train the same iGNN under the same
hyperparameter settings. Figure 6 shows that our method tends to yield better
performances when the true conditional distribution is non-convex, which is
also verified in Table 3 in terms of the smaller two-sample testing
statistics.

a iGNN losses

b iGNN $\hat{X}|Y$

c CINN $\hat{X}|Y$

Figure 6: 2D non-convex node feature visual comparison, iGNN vs. CINN-MMD. The
setup is nearly identical to Figure 5. The conditional generative performance
by iGNN tends to outperform that by CINN-MMD, which fails to capture the non-
convexity of feature distributions in some instances.

Table 3: 2D non-convex node feature two-sample testing statistics. Statistics
are weighted by number of observations at each potential value of $Y$. The
smaller test statistics by iGNN further verify the visual comparisons in
Figure 6.

| MMD: $\alpha=0.1$ | MMD: $\alpha=1.0$ | MMD: $\alpha=5.0$ | MMD: $\alpha=10.0$ | Energy
---|---|---|---|---|---
iGNN | 0.001 | 0.006 | 0.005 | 0.003 | 0.008
CINN-MMD | 0.003 | 0.012 | 0.011 | 0.007 | 0.020

a Raw data

b iGNN loss

c iGNN $\hat{X}|Y$

d CINN $\hat{X}|Y$

Figure 7: Solar ramping event, visual comparison. We show the network and raw
features in (a), loss and error metrics in (b), and generative performances of
the 20-dimensional node feature matrices (10 nodes $\times$ 2D features) in
(c) or (d), where rows are ordered by values of $Y\in\\{0,1\\}^{10}$ with top
three occurrences. We see that the generative performances by ours are
satisfactory in general, but CINN-MMD can fail to generate the non-convex
distribution $X|Y$ in the second row of (d).

2\. Real-data: solar ramping events. We then consider the anomaly detection
task on California solar data in 2017 and 2018. The whole network has ten
sensors with bi-hourly recordings. 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. The first 75% (resp. res 25%) data denote the training
(resp. test) data, which have 1094 (resp. 365) graph observations. Figure 7
visualizes the graph and the raw features colored by anomaly status (i.e.,
black means normal). During training, we fix the learning rate at 1e-4 and
train with a batch size of 150. We build iGNN to have 40 residual blocks,
where we considered in each block three different architectures: (a)
$\text{FC}\circ\text{FC}\circ\text{Chebnet}$ of order 2, (b)
$\text{FC}\circ\text{FC}\circ\text{L3net}$ of neighborhood order (0, 1, 2),
and (c) $\text{FC}\circ\text{FC}\circ\text{FC}$. In terms of generative
performance, architecture (a) performs the best, whose results are shown in
Figure 7 upon comparing with CINN-MMD; we show results for architecture (b)
and (c) in Appendix A.2, Figure 15 and 16. In particular, Figure 7 shows that
when iGNN captures the entire distributional pattern for different $X|Y$,
CINN-MMD fails to do so when the underlying distribution is non-convex (see
the second row in (c) or (d)). Meanwhile, Table 4 shows that both CINN-Nflow
and cGAN yield much larger test statistics than iGNN and CINN-MMD, indicating
much poorer generative performance. Although CINN-MMD may yield smaller
training statistics than iGNN, the reason is likely attributed to how we
weight the statistics for different $X|Y$ and high variance under insufficient
samples. Details can be found in Table 6, Appendix A.2 where we show the test
statistics at each $X|Y$ without weighting.

Table 4: Solar ramping event two-sample testing statistics. Statistics are
weighted by the number of observations at each potential value of $Y$. Similar
to simulation results, CINN-Nflow and cGAN yield much larger test statistics
than iGNN and CINN-MMD, indicating much poorer generative performance. The
quantitative results by ours and CINN-MMD can be similar, but Figure 7 shows
that CINN-MMD cannot accurately predict the distribution $X|Y$ well when it is
non-convex.

Training |  MMD: $\alpha=0.1$ |  MMD: $\alpha=1.0$ |  MMD: $\alpha=5.0$ |  MMD: $\alpha=10$ |  Energy
---|---|---|---|---|---
iGNN |  0.053 |  0.037 |  0.010 |  0.005 |  0.257
CINN-MMD |  0.037 |  0.039 |  0.011 |  0.005 |  0.197
CINN-Nflow |  0.348 |  0.067 |  0.012 |  0.005 |  2.907
cGAN |  0.453 |  0.874 |  0.988 |  0.997 |  2.683

Test |  MMD: $\alpha=0.1$ |  MMD: $\alpha=1.0$ |  MMD: $\alpha=5.0$ |  MMD: $\alpha=10$ |  Energy
---|---|---|---|---|---
Ours |  0.062 |  0.063 |  0.014 |  0.006 |  0.341
CINN-MMD |  0.061 |  0.056 |  0.014 |  0.006 |  0.344
CINN-Nflow |  0.402 |  0.091 |  0.015 |  0.006 |  3.488
cGAN |  0.572 |  0.938 |  0.997 |  1.000 |  3.422

3\. Real-data: traffic flow anomalies. We lastly consider the anomaly
detection task on Los Angeles traffic flow data from April to September 2019.
The whole network has 15 sensors with hourly recordings. 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 first 75% (resp. res 25%) data denote the
training (resp. test) data, which have 3070 (resp. 1320) graph observations.
Figure 8 visualizes the graph and the raw training features, which have
different anomaly patterns from the solar ramping events in Figure 7. During
training, we fix the learning rate at 1e-4 and train with a batch size of 200.
We build iGNN to have 40 residual blocks, where each block contains four FC
layers; empirical evidence suggests that combining FC layers with graph
filters yields poorer generative performance, which we will examine further in
the future. Figure 8 shows that our generative performance is satisfactory in
general, whereas CINN-MMD tends to generate distributions much more dispersed
than the actual ones. Table 5 further verifies this through the much larger
two-sample energy-based test statistics by CINN-MMD.

######  Remark 4.2 (Further implications).

For many real-world graph datasets, test data may contain node label vectors
that lie outside the set of node label vectors in training data, but we can
still estimate the distribution of $X|Y$ upon sampling from $H|Y$ and applying
$F_{\Theta}^{-1}$. In words, this tells us “what circumstances (i.e., $X$)
most likely result in (previously unseen) outcomes (i.e., $Y$)”. In the
context of anomaly detection on graph, specific nodes may only contain normal
observations in training data, but for prognosis, we can examine what features
would lead to anomaly outcomes at these nodes.

Table 5: Traffic anomaly detection two-sample testing statistics. Statistics
are weighted by the number of observations at each potential value of $Y$.
Similar to before, both CINN-Nflow and cGAN yield much larger test statistics
than iGNN and CINN-MMD, indicating much poorer generative performance. In
terms of the energy-based metric, CINN-MMD is also significantly larger than
iGNN, coinciding with our observation in Figure 8. Even under the MMD metric
on which the objective of CINN-MMD is based, the iGNN can still sometimes
reach smaller values.

Training |  MMD: $\alpha=0.1$ |  MMD: $\alpha=1.0$ |  MMD: $\alpha=5.0$ |  MMD: $\alpha=10$ |  Energy
---|---|---|---|---|---
iGNN |  0.058 |  0.037 |  0.010 |  0.005 |  0.214
CINN-MMD |  0.140 |  0.053 |  0.011 |  0.005 |  1.327
CINN-Nflow |  0.316 |  0.055 |  0.011 |  0.005 |  6.487
cGAN |  0.823 |  1.052 |  1.011 |  1.005 |  3.856

Test |  MMD: $\alpha=0.1$ |  MMD: $\alpha=1.0$ |  MMD: $\alpha=5.0$ |  MMD: $\alpha=10$ |  Energy
---|---|---|---|---|---
iGNN |  0.128 |  0.059 |  0.008 |  0.003 |  0.537
CINN-MMD |  0.152 |  0.056 |  0.008 |  0.003 |  1.484
CINN-Nflow |  0.281 |  0.058 |  0.008 |  0.003 |  6.183
cGAN |  0.916 |  1.056 |  1.008 |  1.003 |  4.132

###  4.3 iGNN expressiveness

We focus on verifying model expressiveness analyses in Section 3.3 with
empirical evidence. We first revisit the three-node example in Figure 5, where
we show that models that lack expressiveness can never correctly generate the
features. Then, we design two examples according to Propositions 3.5 and 3.8.
In the first case, the covariance matrix $\Sigma$ is similar to the output of
a spectral-based graph filter. In the second case, $\Sigma$ satisfies certain
spatial locality properties. We show graph filters that have enough
expressiveness can accurately estimate the covariance matrix $\Sigma$ of $X$
and $I_{d}$ of $H$, hence transporting the densities correctly.

a Raw data

b iGNN loss

c iGNN $\hat{X}|Y$

d CINN $\hat{X}|Y$

Figure 8: Traffic anomaly detection, visual comparison. We show the network
and raw features in (a), loss and error metrics in (b), and generative
performances of the 30-dimensional node feature matrices (15 nodes $\times$ 2D
features) in (c) and (d), where rows are ordered by values of
$Y\in\\{0,1\\}^{15}$ with top three occurrences. Our generative performances
are generally satisfactory, but CINN-MMD tends to yield much more dispersed
distributions even when $X|Y$ are convex, possibly due to the scalability
issue in high-dimensional spaces.

1\. GNN lacks expressiveness. Recall Figure 5 on the three-node graph in
Example 1, where we showed the satisfactory performance of our iGNN. In
particular, each residual block takes the form
$\text{FC}\circ\text{FC}\circ\text{L3net}$. Figure 9 shows the results where
we replace the first L3net in each residual block by a Chebnet of order 2 or
by a FC layer. In particular, iGNN under the Chebnet variant cannot learn the
distribution $X|Y$, due to the symmetry issue inherent in spectral-based graph
filters mentioned in Example 1.

a Cheb

b Cheb generation

c FC

d FC generation

Figure 9: Results for 2D convex node features. We replace the L3net layer with
the Chebnet layer under an identical setup as Figure 5. The conditional
generative performances are satisfactory by using FC layers (which have enough
model expressiveness in theory) but not by using a Chebnet layer (spectral-
based with symmetric graph Laplacian).

a Network

b Loss $L_{g}$ in (7)

c True vs. estimated correlation matrices.

Figure 10: Results for the spectral-based $\Sigma$ as in Proposition 3.5.
Losses converge in a few epochs, and correlation matrices are estimated with
high accuracy, revealing that iGNN with enough expressiveness can be trained
to estimate the true mapping.

a Network

b Loss $L_{g}$ in (7)

c True vs. estimated correlation matrices.

Figure 11: Results for the local $\Sigma$ as in Proposition 3.8. Losses
converge reasonably fast, and correlation matrices are estimated with high
accuracy, revealing that iGNN with enough expressiveness can be trained to
estimate the true mapping.

2\. Sufficient GNN expressiveness for spectral-based $\Sigma$, see Proposition
3.5. We first consider the case in which
$\Sigma=\sum_{k=0}^{2}a_{k}T_{k}(\widetilde{L})$, where $T_{k}$ denotes the
$k$-th degree Chebyshev polynomial and $\widetilde{L}$ is the scaled and
normalized graph Laplacian. We let $a_{0}=0.5,a_{1}=0.1,a_{2}=0.5$ so that
$\Sigma$ is positive definite with a significant number of non-zero and large
off-diagonal entries. We consider a seven-node chordal cycle graph (Lubotzky,
1994), which is one example of expander graphs. Figure 10 visualizes the raw
graph. To verify the existence statement, we design iGNN to have 40 residual
blocks, where each residual block is a single Chebnet filter of order two
(Defferrard et al., 2016) mapping from $\mathbb{R}$ to $\mathbb{R}$. During
training, we generate 4000 i.i.d. training and test samples following the GP,
fix the learning rate at 5e-4, and train with a batch size of 400. Figure 10
shows that both the training and test losses converge in a small number of
training epochs. Figure 10 shows that the correlation matrices $\Sigma$ and
$I_{V}$ are both estimated with high accuracy.

3\. Sufficient GNN expressiveness for spatial-based $\Sigma$ and
$\Sigma^{-1}$, see Proposition 3.8. We then consider the three-node graph in
Example 1, where node 1 is connected to both node 0 and 2, neither of which
connects each other. In particular, we specify the correlation between node 0
and 1 (resp. 2 and 1) as 0.6 (resp. $-0.4$), so that $\Sigma$ (resp.
$\Sigma^{-1}$) is one-neighbor local (resp. two-neighbor local). Figure 11
visualizes the raw graph. We also use a 40-block iGNN for estimation, where
the only difference is that an L3net layer replaces each Chebnet layer with
neighborhood orders $(d_{1},d_{2},d_{3})=(0,1,2)$. All other hyperparameters
are the same as those in section 4.1.1 above. Figure 11 shows that the losses
converge reasonably fast. Figure 11 shows that the correlation matrices are
estimated with high accuracy. Lastly, note that we intentionally design the
correlation matrix $\Sigma$ to be asymmetric as in Example 1, whereas the
graph Laplacian is symmetric due to graph design. As a result, any spectral-
based graph filter is still unable to express the matrix $\Sigma$ and,
therefore, perform the generative task. Results are shown in Appendix A.2,
Figures 13 and 13, which show that even if losses converge, neither
correlation matrices is estimated correctly.

##  5 Conclusion

We proposed a general framework based on normalizing flow for conditional
generation in this work. Although our primary focus is on graph data, the
framework is general for any Euclidean data, which may be viewed as
observations on graphs with a single node, so iGNN naturally applies to these
problems. This framework benefits one-to-many inverse problems that aim to
infer $X|Y$. In particular, our techniques address both prediction and
conditional generation through a single invertible residual neural work at
once, thus allowing for implementation simplicity and clarity in training.
Theoretically, we analyzed the identifiability of a true invertible mapping
assuming its existence, verified its existence and invertibility utilizing the
Fokker-Planck equation of an Ornstein-Uhlenbeck process, and examined the
expressiveness of iGNN in learning the mapping. Experimentally, we compared
iGNN with competing conditional generative models on simulated and real-data
examples, verifying the improved performance of iGNN, especially in high-
dimensional and non-convex graph cases.

There are many open avenues for research. Method-wise, we want to consider
continuous outcomes $Y$ and other generative tasks such as edge feature
generation and graph generation. Regarding network training, we want to
consider the continuous-time variant of residual networks (Grathwohl et al.,
2019) and utilize techniques by monotone variational inequalities in each
residual block for performance guarantees (Xu et al., 2022). Theory-wise, we
will further the analyses regarding the invertibility of the true force ${\bf
f}(x,t)$ in (10) based on the Jordan-Kinderlehrer-Otto (JKO) proximal scheme
regularization. We also want to study the expressiveness of iGNN in estimating
the true force more deeply.

## 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.

##  6 Proofs

###### Proof of Proposition 3.1.

Recall that $F^{*}$ is bijective in $\mathbb{R}^{d}$ and maps $X|Y$ to $H|Y$
for each potential value of $Y$. Hence, for each value $i\in[K]$, the random
vector $\hat{X}_{i}:=F^{*-1}(H|Y=i)$ follows the same distribution as $X|Y=i$.
Therefore, the random variables generated by the inverse of $F^{*}$ follow the
same conditional distributions $X|Y$. In addition, since we assume that there
exists a $\Theta^{*}$ at which the model $F_{\Theta^{*}}$ is identical to
$F^{*}$, $F_{\Theta^{*}}^{-1}(H|Y)$ also follows the true conditional
distribution.

By the assumption that $N_{K}$ approaches infinity as $N\rightarrow\infty$, we
know that infinitely many samples from each $X|Y=i$, $i\in[K]$ are available.
Under additional regularity conditions, we know that the MLE estimate
$\hat{\Theta}_{N}$ is also consistent (Newey & McFadden, 1986)—it converges to
the true parameter $\Theta^{*}$ with probability one. Therefore, the model
$F_{\hat{\Theta}_{N}}^{-1}\rightarrow F_{\Theta^{*}}^{-1}$ with probability
one. ∎

###### Proof of Lemma 3.2.

The solution $\rho(x,t)$ of (11) has the explicit expression as

$\rho(x,t)=\int K_{t}(x,y)\rho_{0}(y)dy,\quad
K_{t}(x,y):=\frac{1}{(2\pi\sigma_{t}^{2})^{d/2}}e^{-\frac{|x-e^{-t}y|^{2}}{2\sigma_{t}^{2}}},\quad\sigma_{t}^{2}:=1-e^{-2t}.$
(14)

By the smoothness of $K_{t}(x,y)$, $\rho(x,t)$ is smooth over
$\mathbb{R}^{d}\times(0,\infty)$. By (12), for any $t>0$, ${\bf
f}(x,t)=-x-\frac{\nabla\rho(x,t)}{\rho(x,t)}$, and thus ${\bf f}(x,t)$ is also
smooth over $\mathbb{R}^{d}\times(0,\infty)$. ∎

###### Proof of Proposition 3.3.

The well-posedness of the IVP can be guaranteed as long as ${\bf f}(x,t)$ is
continuous and satisfies a Lipschitz condition w.r.t $x$ (Sideris, 2013). The
former is by Lemma 3.2, and the latter is the assumption of the proposition. ∎

###### Proof of Lemma 3.4.

We only need to derive the distribution of $X_{t}$ in order to find its
density function $\rho(x,t)$. Because the Ornstein–Uhlenbeck process is a
Gaussian process, it is evident that $X_{t}\sim N(0,\Sigma_{t})$, where
$\Sigma_{t}$ is the time-dependent covariance matrix.

By the transition probability, also known as the Green’s function

$K(X,Y;t):=(2\pi\sigma_{t}^{2})^{-d/2}\exp\left(-(2\sigma_{t}^{2})^{-1|}\|X-\exp(-t)Y\|^{2}_{2}\right),\sigma_{t}^{2}:=1-\exp(-2t),$

we know that $X_{t}|X_{0}=Y\overset{d}{=}N(0,\sigma_{t}^{2}I_{d})$. When
$X_{0}\sim N(0,\Sigma)$, we thus have

$X_{t}=\exp(-t)X_{0}+\mathcal{E},\mathcal{E}\sim
N(0,\sigma_{t}^{2}I_{d}),X_{0}\perp\mathcal{E}.$

Hence, $\Sigma_{t}:=Var(X_{t})=\exp(-2t)\Sigma+(1-\exp(-2t))I_{d}$.

As a result, $\rho(x,t)\propto\exp(-1/2X_{t}^{T}\Sigma_{t}^{-1}X_{t})$, with
$\nabla_{X}\ln\rho(x,t)=-\Sigma_{t}^{-1}X_{t}$. ∎

###### Proof of Proposition 3.5.

Note that by the spectral decomposition $L=U\Lambda U^{T}$ and the assumption
$\Sigma=P(L)$, we have

$\displaystyle\Sigma$
$\displaystyle=\sum_{i=1}^{V}P(\lambda_{i}(L))U_{i}U_{i}^{T},$
$\displaystyle\Sigma^{-1}$
$\displaystyle=\sum_{i=1}^{V}P(\lambda_{i}(L))^{-1}U_{i}U_{i}^{T}.$

Thus,

$\displaystyle\Sigma_{t}$ $\displaystyle=(1-\exp(-2t))I_{V}+\exp(-2t)\Sigma$
$\displaystyle=U[(1-\exp(-2t))I_{V}+\exp(-2t)P(\Lambda)]U^{T}.$

As a result,

$\Sigma_{t}^{-1}=\sum_{i=1}^{V}\left((1-\exp(-2t))+\exp(-2t)P(\lambda_{i}(L))\right)^{-1}U_{i}U_{i}^{T}.$

∎

###### Proof of Corollary 3.6.

By the assumption that $\Sigma$ and $L$ share the same eigen-bases, we have

$\displaystyle\Sigma_{t}$
$\displaystyle=\sum_{i=1}^{V}\left((1-\exp(-2t))+\exp(-2t)\lambda_{i}(\Sigma)\right)U_{i}U_{i}^{T},$
$\displaystyle\tilde{\Sigma}_{t}$
$\displaystyle=\sum_{i=1}^{V}\left((1-\exp(-2t))+\exp(-2t)\tilde{P}(\lambda_{i}(L))\right)U_{i}U_{i}^{T}.$

Under diagonalization, the inverse of $\Sigma_{t}$ and $\tilde{\Sigma}_{t}$
can be expressed as

$\displaystyle\Sigma_{t}^{-1}$
$\displaystyle=\sum_{i=1}^{V}\left((1-\exp(-2t))+\exp(-2t)\lambda_{i}(\Sigma)\right)^{-1}U_{i}U_{i}^{T},$
$\displaystyle\tilde{\Sigma}_{t}^{-1}$
$\displaystyle=\sum_{i=1}^{V}\left((1-\exp(-2t))+\exp(-2t)\tilde{P}(\lambda_{i}(L))\right)^{-1}U_{i}U_{i}^{T}.$

Thus, we have

$\displaystyle\rho(\Sigma_{t}^{-1}-\tilde{\Sigma}_{t}^{-1})$
$\displaystyle\leq$
$\displaystyle\max_{i}\mid\left((1-\exp(-2t))+\exp(-2t)\lambda_{i}(\Sigma)\right)^{-1}-\left((1-\exp(-2t))+\exp(-2t)\tilde{P}(\lambda_{i}(L))\right)^{-1}\mid$
$\displaystyle\leq$
$\displaystyle\left|\frac{(1-c_{t})\delta}{(c_{t}+(1-c_{t})\lambda_{d}(\Sigma))(c_{t}+(1-c_{t})\min_{i}\tilde{P}(\lambda_{i}(L)))}\right|,\quad
c_{t}:=1-\exp(-2t).$

Note that the denominator
$(c_{t}+(1-c_{t})\lambda_{d}(\Sigma))(c_{t}+(1-c_{t})\min_{i}\tilde{P}(\lambda_{i}(L)))$
is the product of two convex combinations of 1 and either
$\lambda_{V}(\Sigma)$ or $\min_{i}\tilde{P}(\lambda_{i}(L))$. Because $\Sigma$
is positive definite, its minimum eigenvalue is bounded away from zero.
Therefore, as long as $\min_{i}\tilde{P}(\lambda_{i}(L))>0$, both
$c_{t}+(1-c_{t})\lambda_{V}(\Sigma)$ and
$c_{t}+(1-c_{t})\min_{i}\tilde{P}(\lambda_{i}(L))$ are bounded above and away
from zero for all $t\geq 0$. Hence, the denominator is bounded above, leading
to

$\rho(\Sigma_{t}^{-1}-\tilde{\Sigma}_{t}^{-1})=O(e^{-2t}\delta),$

which denotes the tail bound of the spectral norm of the matrix difference. ∎

###### Proof of Proposition 3.8.

Recall that $\Sigma_{t}:=(1-\exp(-2t))I_{V}+\exp(-2t)\Sigma$. We thus analyze
the locality of $\Sigma_{t}$ based on the magnitude of $t$.

When $t=0$, $\Sigma_{t}=\Sigma$ so that $\Sigma_{t}=\Sigma^{-1}$.

When $t>0$, we consider two cases:

(1) Long-time: Assume $\exp(-2t)<=1/2$ and denote
$c_{t}:=\exp(-2t)/(1-\exp(-2t))$. Then,

$\Sigma_{t}=(1-\exp(-2t)[I_{V}+c_{t}\Sigma],$

which yields

$\Sigma_{t}^{-1}=(1-\exp(-2t))[\sum_{p=0}^{\infty}(-1)^{p}(c_{t}\Sigma)^{p}]$
(15)

by the power series expansion $(I+A)^{-1}=\sum_{p=0}^{\infty}(-1)^{p}A^{p}$ if
the spectral norm $\rho(A)<1$. We may assume that $\rho(c_{t}\Sigma)$ is small
enough so that this hold, especially given that $c_{t}$ converges to zero
exponentially fast.

Therefore, if $\Sigma$ is local, we can truncate (15) up to the first $k$
summands. Meanwhile, note that when $\Sigma$ is $v_{1}$-local, $\Sigma^{k}$ is
$kv_{1}$-local, as each matrix power yields a non-zero entry between node $i$
and node $k$ if and only if there is at least one edge between $i$ and $j$ and
one edge between $j$ and $k$, where $j$ is another node.

(2) Short-time: Assume $\exp(-2t)>1/2$. The proof is nearly the same as in the
previous case, where we have

$\Sigma_{t}=\exp(-2t)\Sigma[I_{V}+c_{t}^{-1}\Sigma^{-1}],$

which yields

$\Sigma_{t}^{-1}=\exp(-2t)[\sum_{p=0}^{\infty}-c_{t}^{-p}\Sigma^{-(p+1)}].$

The remaining steps are identical to above ones, where we use the
$v_{2}$-locality of $\Sigma^{-1}$ to establish the $kv_{2}$-locality of
$\Sigma_{t}^{-1}$. ∎

## References

  * Al-Shaalan (2019) Al-Shaalan, A. M.  Reliability evaluation of power systems.  2019\.
  * Ardizzone et al. (2019a) Ardizzone, L., Kruse, J., Wirkert, S. J., Rahner, D., Pellegrini, E., Klessen, R. S., Maier-Hein, L., Rother, C., and Köthe, U.  Analyzing inverse problems with invertible neural networks.  _ArXiv_ , abs/1808.04730, 2019a.
  * Ardizzone et al. (2019b) Ardizzone, L., Lüth, C., Kruse, J., Rother, C., and Köthe, U.  Guided image generation with conditional invertible neural networks.  _arXiv preprint arXiv:1907.02392_ , 2019b.
  * Behrmann et al. (2019) Behrmann, J., Grathwohl, W., Chen, R. T., Duvenaud, D., and Jacobsen, J.-H.  Invertible residual networks.  In _International Conference on Machine Learning_ , pp. 573–582. PMLR, 2019.
  * Bolley et al. (2012) Bolley, F., Gentil, I., and Guillin, A.  Convergence to equilibrium in wasserstein distance for fokker–planck equations.  _Journal of Functional Analysis_ , 263(8):2430–2457, 2012.
  * Cheng et al. (2021) Cheng, X., Miao, Z., and Qiu, Q.  Graph convolution with low-rank learnable local filters.  In _International Conference on Learning Representations_ , 2021.
  * Clevert et al. (2015) Clevert, D.-A., Unterthiner, T., and Hochreiter, S.  Fast and accurate deep network learning by exponential linear units (elus).  _arXiv preprint arXiv:1511.07289_ , 2015.
  * Defferrard et al. (2016) Defferrard, M., Bresson, X., and Vandergheynst, P.  Convolutional neural networks on graphs with fast localized spectral filtering.  _Advances in neural information processing systems_ , 29, 2016.
  * Dinh et al. (2017) Dinh, L., Sohl-Dickstein, J., and Bengio, S.  Density estimation using real nvp.  _ArXiv_ , abs/1605.08803, 2017.
  * Goodfellow et al. (2014) Goodfellow, I. J., Pouget-Abadie, J., Mirza, M., Xu, B., Warde-Farley, D., Ozair, S., Courville, A. C., and Bengio, Y.  Generative adversarial nets.  In _NIPS_ , 2014.
  * Grathwohl et al. (2019) Grathwohl, W., Chen, R. T. Q., Bettencourt, J., Sutskever, I., and Duvenaud, D. K.  Ffjord: Free-form continuous dynamics for scalable reversible generative models.  _ArXiv_ , abs/1810.01367, 2019.
  * Gretton et al. (2012) Gretton, A., Borgwardt, K. M., Rasch, M. J., Schölkopf, B., and Smola, A.  A kernel two-sample test.  _J. Mach. Learn. Res._ , 13:723–773, 2012.
  * Gulrajani et al. (2017) Gulrajani, I., Ahmed, F., Arjovsky, M., Dumoulin, V., and Courville, A. C.  Improved training of wasserstein gans.  In _NIPS_ , 2017.
  * Haykin (1998) Haykin, S.  Neural networks: A comprehensive foundation.  1998\.
  * Ioffe & Szegedy (2015) Ioffe, S. and Szegedy, C.  Batch normalization: Accelerating deep network training by reducing internal covariate shift.  _ArXiv_ , abs/1502.03167, 2015.
  * Isola et al. (2017) Isola, P., Zhu, J.-Y., Zhou, T., and Efros, A. A.  Image-to-image translation with conditional adversarial networks.  _2017 IEEE Conference on Computer Vision and Pattern Recognition (CVPR)_ , pp. 5967–5976, 2017.
  * Kingma & Ba (2015) Kingma, D. P. and Ba, J.  Adam: A method for stochastic optimization.  _CoRR_ , abs/1412.6980, 2015.
  * Kingma & Welling (2014) Kingma, D. P. and Welling, M.  Auto-encoding variational bayes.  _CoRR_ , abs/1312.6114, 2014.
  * Kingma & Welling (2019) Kingma, D. P. and Welling, M.  An introduction to variational autoencoders.  _Foundations and Trends® in Machine Learning_ , 12(4):307–392, 2019.  doi: 10.1561/2200000056.
  * Kobyzev et al. (2020) Kobyzev, I., Prince, S. J., and Brubaker, M. A.  Normalizing flows: An introduction and review of current methods.  _IEEE transactions on pattern analysis and machine intelligence_ , 43(11):3964–3979, 2020.
  * Ledig et al. (2017) Ledig, C., Theis, L., Huszár, F., Caballero, J., Aitken, A. P., Tejani, A., Totz, J., Wang, Z., and Shi, W.  Photo-realistic single image super-resolution using a generative adversarial network.  _2017 IEEE Conference on Computer Vision and Pattern Recognition (CVPR)_ , pp. 105–114, 2017.
  * Lu et al. (2015) Lu, J., Wu, D., Mao, M., Wang, W., and Zhang, G.  Recommender system application developments: A survey.  _Decis. Support Syst._ , 74:12–32, 2015.
  * Lu et al. (2017) Lu, Z., Pu, H., Wang, F., Hu, Z., and Wang, L.  The expressive power of neural networks: A view from the width.  In _NIPS_ , 2017.
  * Lubotzky (1994) Lubotzky, A.  Discrete groups, expanding graphs and invariant measures.  In _Progress in mathematics_ , 1994.
  * Lucas et al. (2019) Lucas, J., Tucker, G., Grosse, R. B., and Norouzi, M.  Understanding posterior collapse in generative latent variable models.  In _DGS@ICLR_ , 2019.
  * Makhzani et al. (2015) Makhzani, A., Shlens, J., Jaitly, N., and Goodfellow, I. J.  Adversarial autoencoders.  _ArXiv_ , abs/1511.05644, 2015.
  * Mirza & Osindero (2014) Mirza, M. and Osindero, S.  Conditional generative adversarial nets.  _ArXiv_ , abs/1411.1784, 2014.
  * Newey & McFadden (1986) Newey, W. and McFadden, D.  Large sample estimation and hypothesis testing.  _Handbook of Econometrics_ , 4:2111–2245, 1986.
  * Salimans et al. (2016) Salimans, T., Goodfellow, I. J., Zaremba, W., Cheung, V., Radford, A., and Chen, X.  Improved techniques for training gans.  _ArXiv_ , abs/1606.03498, 2016.
  * Sánchez-Lengeling & Aspuru-Guzik (2018) Sánchez-Lengeling, B. and Aspuru-Guzik, A.  Inverse molecular design using machine learning: Generative models for matter engineering.  _Science_ , 361:360 – 365, 2018.
  * Sideris (2013) Sideris, T. C.  Ordinary differential equations and dynamical systems.  2013\.
  * Székely & Rizzo (2013) Székely, G. J. and Rizzo, M. L.  Energy statistics: A class of statistics based on distances.  _Journal of Statistical Planning and Inference_ , 143:1249–1272, 2013.
  * Wehenkel & Louppe (2019) Wehenkel, A. and Louppe, G.  Unconstrained monotonic neural networks.  _ArXiv_ , abs/1908.05164, 2019.
  * Xu et al. (2022) Xu, C., Cheng, X., and Xie, Y.  Training neural networks using monotone variational inequality.  _arXiv preprint arXiv:2202.08876_ , 2022.
  * Zhu et al. (2017) Zhu, J.-Y., Park, T., Isola, P., and Efros, A. A.  Unpaired image-to-image translation using cycle-consistent adversarial networks.  _2017 IEEE International Conference on Computer Vision (ICCV)_ , pp. 2242–2251, 2017.
  * Ziegler & Rush (2019) Ziegler, Z. M. and Rush, A. M.  Latent normalizing flows for discrete sequences.  In _ICML_ , 2019.
  * Zong et al. (2018) Zong, B., Song, Q., Min, M. R., Cheng, W., Lumezanu, C., ki Cho, D., and Chen, H.  Deep autoencoding gaussian mixture model for unsupervised anomaly detection.  In _ICLR_ , 2018.

##  Appendix A Additional experiments

###  A.1 Setup details

1\. Verify model invertibility We see from Figure 12 that iGNN under the
Wasserstein-2 regularization ensures model invertibility up to very high
accuracy.

a In-distribution sample $X$

b Out-of-distribution sample $X$

Figure 12: iGNN invertibility assessment. Given $X$ and a fixed iResnet
architecture, examine $F^{-1}_{\Theta}(F_{\Theta}(X))$ (forward then invert)
vs. $F_{\Theta}(F^{-1}_{\Theta}(X))$ (invert then forward) on in-distribution
data (two-moon) and out-of-distribution data (random $U[0,1)$ samples).

2\. Comparison methods We compare iGNN with three conditional generative
models, which are representative of different approaches to this task: the
first builds on Generative Adversarial Network (GAN) (Goodfellow et al.,
2014), the second uses invertible neural networks to analyze inverse problems
(Ardizzone et al., 2019a), and the last extends the idea of normalizing flow
(Ardizzone et al., 2019b). These works were originally developed for image
generation tasks, but we apply them here due to the generality of these
methods.

2.1. Conditional Generative Adversarial Network (cGAN) (Isola et al., 2017):
Denote $Y\in[K]$ as the one-hot encoded categorical response,
$X_{Y}\in\mathbb{R}^{n}$ as the corresponding true feature, and
$Z\in\mathbb{R}^{n}$ as the random noise. Given a generator
$G:(Z,Y)\rightarrow\mathbb{R}^{n}$ and a discriminator
$D:(X_{Y},Y)\rightarrow[0,1]$, the objective of the cGAN loss is

$\displaystyle\arg\min_{G}\max_{D}\mathcal{L}_{cGAN}(G,D)+\lambda\mathcal{L}_{2}(G)$
(16) $\displaystyle\mathcal{L}_{cGAN}(G,D):=\mathbb{E}_{X_{Y},Y}[\log
D(X_{Y},Y)]+\mathbb{E}_{Z,Y}[\log(1-D(G(Z,Y),Y))]$ (17)
$\displaystyle\mathcal{L}_{2}(G):=\mathbb{E}_{X_{Y},Z,Y}[\|X_{Y}-G(Z,Y)\|_{2}]$
(18)

We note that the discriminator $D$ aims to maximize (17) in order to
distinguish fake features $G(Z,Y)$ from the true ones $X_{Y}$, whereas the
generator $G$ tries to fool $D$ by minimizing (17). Therefore, this min-max
game leads to objective (16), which also includes an $\ell_{2}$ penalty to
force $G$ producing features as similar to $X_{Y}$ as possible.

During maximization, we adopt the recommendation in (Isola et al., 2017;
Goodfellow et al., 2014) by training $G$ to maximize $\mathbb{E}_{Z,Y}[\log
D(G(Z,Y),Y)]$. We also divide the objective by 2 when optimizing $D$, so as to
slow down the rate it learns comparing to $G$.

We adopt the following hyperparameters and model architectures upon visual
cross-validation (i.e., comparing the generative performances). Regarding
architecture, we aim to use networks as similar to our iGNN as possible for
fair comparisons. In particular, we concatenate blocks of layers together but
do not restrict the output of each block, as we need not worry about model
invertibility; a batch normalization layer (Ioffe & Szegedy, 2015) is used
after each block. Because cGAN is not flow-based, we keep the hidden dimension
of all except the last layer in the last block at 64, rather than forcing the
output of each block to lie within the original dimension (as when training
our iGNN). In practice, the generative task seems harder than the
discriminative task, so we always make $G$ twice deeper as $D$, with the
former having 40 blocks. Regarding hyperparameters, we fix $\lambda=5$ in (16)
and always use the Adam optimizer under weight decay equal to 0.5 (i.e.,
$l_{2}$ penalty on model parameters). The learning rates for network $G$ and
network $D$ are fixed as (a) 5e-4 and 1e-4 for generating 2D convex node
features and the solar anomaly detection (b) 1e-1 and 1e-3 for the traffic
anomaly detection. All other hyperparameters (e.g., structure of each block,
batch size, termination criteria, etc.) are the same as we did in our training
of (16).

2.2. Conditional Invertible Neural Networks with Maximum Mean Discrepancy
(CINN-MMD) (Ardizzone et al., 2019a): This method constructs conditional
coupling blocks as follows. Denote $Y$ as one-hot encoded categorical feature
and $X_{b}$ as the input to the $b-$th block $f_{b}$, which first splits
$X_{b}$ into two parts $[X_{b,1},X_{b,2}]$. Then, the output by $f_{b}$ is
$[X_{b+1,1},X_{b+1,2}]$, where

$\displaystyle X_{b+1,1}$
$\displaystyle=X_{b,1}\odot\exp(s_{b,1}(X_{b,2},Y))+t_{b,1}(X_{b,2},Y)),$ (19)
$\displaystyle X_{b+1,2}$
$\displaystyle=X_{b,2}\odot\exp(s_{b,2}(X_{b+1,1},Y))+t_{b,2}(X_{b+1,1},Y)),$

where the subnetworks $\\{s_{b,1},s_{b,2},t_{b,1},t_{b,2}\\}$ can be arbitrary
neural networks. The full network $f:=f_{B}\circ\ldots\circ f_{1}$ thus
concatenates $B$ such blocks and takes in input $X_{0}=X$ and outputs $X_{B}$.

Given a CINN as the concatenation of $B$ blocks, the method optimizes the
following three losses where $Z$ follows a standard multivariate Gaussian:

$\displaystyle\mathcal{L}_{Y}$ $\displaystyle:=\mathbb{E}[(Y-f_{Y}(X))^{2},]$
$\displaystyle\mathcal{L}_{Z}$ $\displaystyle:=D(q(Y,Z),p(Y)p(Z)),$
$\displaystyle\mathcal{L}_{X}$ $\displaystyle:=D(q(X),p(X)).$

Note that (1) $\mathcal{L}_{y}$ measures the prediction error between true $Y$
and $f_{Y}(X)$, where $f(X)=[f_{Y}(X),f_{Z}(X)]$ denotes the full prediction.
(2) $\mathcal{L}_{Z}$ (resp. $\mathcal{L}_{X}$) measures the distance between
the joint distribution of $(Y,Z)$ based on prediction (resp. the distribution
based on generated samples) and the product of their true marginal
distributions (resp. true distribution of $X$). The notation $q$ denotes the
estimated distribution based on $f(X)$ and $p$ denotes the true distribution
for variables $Y$ and $Z$. The distance function $D$ is chosen as the MMD
distance under the Inverse Multiquadratic kernel. In addition to this setup,
(Ardizzone et al., 2019a) suggested several tricks during optimization and
model construction, including blocking the gradient of $\mathcal{L}_{Z}$ with
respect to $Y$, balancing the effect of three losses, and random shuffling of
hidden variables.

We adopt the following hyperparameters and model architectures upon visual
cross-validation (i.e., comparing the generative performances). Regarding
architecture, 40 blocks are used, and each subnetwork has three FC layers with
hidden dimensions fixed at 128; we did not use a graph filter because of the
splitting of features at each block. Splitted features are always assumed to
have the same dimension. Regarding hyperparameters, we fix $\lambda=0.5$ in
(16) and always use the Adam optimizer with learning rate at 5e-4. All other
hyperparameters (e.g., batch size, termination criteria, etc.) are the same as
we did in our training of (16).

2.3. Conditional Invertible Neural Networks with Normalizing Flow (CINN-Nflow)
(Ardizzone et al., 2019b): Given a CINN as the concatenation of $B$ blocks,
the method directly optimizes the unconditional generative loss in (4) using
$X_{B}$. It also adds an $\ell_{2}$ weight regularization for better practical
performance. The construction of the conditional coupling blocks is identical
to CINN-MMD, so that we omit the descriptions. The same hyperparameters as
used earlier are also adopted.

###  A.2 Further results

a Network

b Loss $L_{g}$ in (7)

c True vs. estimated correlation matrices.

Figure 13: IGNN results for the local $\Sigma$ as in Proposition 3.8. Each
residual block in iGNN is the Chebnet. We see that losses converge reasonably
fast, but correlation matrices are not estimated well, revealing the effect of
GNNs lacking expressiveness in generative tasks. This example contrasts with
iGNN under L3net in Figure 11, which shows satisfactory generative
performances when GNN has enough expressiveness.

a 3 node, convex

b 3 node, non-convex

c Solar

d Traffic

e 8 Gaussians

f 3 node, non-convex

g Solar

h Traffic

i 8 Gaussians

j 3 node, convex

k Solar

l Traffic

Figure 14: Losses by competing methods on simulated and real-data experiments.
We see that losses all converges but generative performance differs. First
row: losses for CINN-MMD. Second row: losses for CINN-Nflow. Last row: losses
for cGAN

a Raw data

b Metrics

c True vs. estimated $X$ and $H$ at the two most frequent choices of $Y$.

Figure 15: CA solar ramping event detection. We replace in iGNN the Chebnet
layer with an L3net layer and show the network and raw features in (a), loss
and error metrics in (b), and generative performances of the 20-dimensional
node feature matrices in (c). The generative performances are similar to the
ones in Figure 7.

a Raw data

b Metrics

c True vs. estimated $X$ and $H$ at the two most frequent choices of $Y$.

Figure 16: CA solar ramping event detection by using FC layers in iGNN only.
We show the network and raw features in (a), loss and error metrics in (b),
and generative performances of the 20-dimensional node feature matrices in
(c). The generative performances are poorer than the earlier ones in Figure 7
and 15.

Table 6: Solar ramping event, two-sample testing statistics at each $X|Y$
before weighting. The first column indicates the outcome of $Y$ on the graph
and the number of graph observations in the data. We see that our method has
very similar test statistics as CINN-MMD, but the statistics can have too much
variance when samples are too few, where it is essential to examine the
generative performances qualitatively as we did in Figure 7.

Training |  | iGNN | CINN-MMD | CINN-Nflow | cGAN
---|---|---|---|---|---
$Y=$(0,0,0,0,0,0,0,0,0,0), 679 obs | MMD: $\alpha$=0.1 | 0.041 | 0.021 | 0.304 | 0.267
| MMD: $\alpha$=1.0 | 0.081 | 0.086 | 0.582 | 1.297
| MMD: $\alpha$=5.0 | 0.231 | 0.210 | 0.197 | 0.668
| MMD: $\alpha$=10.0 | 0.019 | 0.017 | 0.040 | 0.800
| Energy | 0.123 | 0.146 | 0.201 | 1.201
$Y=$(1,1,1,1,1,1,1,1,1,1), 144 obs | MMD: $\alpha$=0.1 | 0.038 | 0.031 | 0.027 | 1.026
| MMD: $\alpha$=1.0 | 0.004 | 0.004 | 0.004 | 0.974
| MMD: $\alpha$=5.0 | 0.041 | 0.046 | 0.047 | 1.047
| MMD: $\alpha$=10.0 | 0.009 | 0.009 | 0.009 | 1.009
| Energy | 0.002 | 0.002 | 0.002 | 0.992
$Y=$(0,0,0,0,0,0,0,0,1,0), 22 obs | MMD: $\alpha$=0.1 | 0.020 | 0.021 | 0.021 | 1.021
| MMD: $\alpha$=1.0 | 0.004 | 0.004 | 0.004 | 1.004
| MMD: $\alpha$=5.0 | 0.189 | 0.075 | 2.418 | 1.895
| MMD: $\alpha$=10.0 | 0.387 | 0.632 | 5.221 | 6.288
| Energy | 1.485 | 1.106 | 2.850 | 3.425

Test |  | iGNN | CINN-MMD | CINN-Nflow | cGAN
---|---|---|---|---|---
$Y=($0,0,0,0,0,0,0,0,0,0) 222 obs | MMD: $\alpha$=0.1 | 0.029 | 0.014 | 0.286 | 0.334
| MMD: $\alpha$=1.0 | 0.150 | 0.172 | 0.774 | 1.380
| MMD: $\alpha$=5.0 | 0.281 | 0.446 | 0.798 | 1.084
| MMD: $\alpha$=10.0 | 0.013 | 0.010 | 0.032 | 0.831
| Energy | 0.207 | 0.199 | 0.274 | 1.273
$Y=($1,1,1,1,1,1,1,1,1,1) 59 obs | MMD: $\alpha$=0.1 | 0.311 | 0.239 | 0.331 | 1.327
| MMD: $\alpha$=1.0 | 0.002 | 0.002 | 0.003 | 0.980
| MMD: $\alpha$=5.0 | 0.048 | 0.050 | 0.051 | 1.051
| MMD: $\alpha$=10.0 | 0.066 | 0.063 | 0.066 | 1.066
| Energy | 0.001 | 0.001 | 0.001 | 0.993
$Y=($1,0,1,1,1,1,1,1,1,1) 10 obs | MMD: $\alpha$=0.1 | 0.021 | 0.021 | 0.021 | 1.021
| MMD: $\alpha$=1.0 | 0.025 | 0.025 | 0.025 | 1.025
| MMD: $\alpha$=5.0 | 0.200 | 0.145 | 2.646 | 2.422
| MMD: $\alpha$=10.0 | 0.677 | 0.799 | 6.374 | 6.877
| Energy | 1.483 | 2.072 | 5.146 | 5.230