# MULTIRESOLUTION EQUIVARIANT GRAPH VARIA## TIONAL AUTOENCODER

**Anonymous authors**
Paper under double-blind review

ABSTRACT

In this paper, we propose Multiresolution Equivariant Graph Variational Autoen_coders (MGVAE), the first hierarchical generative model to learn and generate_
graphs in a multiresolution and equivariant manner. At each resolution level, MGVAE employs higher order message passing to encode the graph while learning
to partition it into mutually exclusive clusters and coarsening into a lower resolution that eventually creates a hierarchy of latent distributions. MGVAE then
constructs a hierarchical generative model to variationally decode into a hierarchy
of coarsened graphs. Importantly, our proposed framework is end-to-end permutation equivariant with respect to node ordering. MGVAE achieves competitive results with several generative tasks including general graph generation, molecular
generation, unsupervised molecular representation learning to predict molecular
properties, link prediction on citation graphs, and graph-based image generation.

1 INTRODUCTION

Understanding graphs in a multiscale and multiresolution perspective is essential for capturing the
structure of molecules, social networks, or the World Wide Web. Graph neural networks (GNNs)
utilizing various ways of generalizing the concept of convolution to graphs (Scarselli et al., 2009)
(Niepert et al., 2016b) (Li et al., 2016) have been widely applied to many learning tasks, including
modeling physical systems (Battaglia et al., 2016), finding molecular representations to estimate
quantum chemical computation (Duvenaud et al., 2015) (Kearnes et al., 2016) (Gilmer et al., 2017b)
(Hy et al., 2018), and protein interface prediction (Fout et al., 2017). One of the most popular
types of GNNs is message passing neural nets (MPNNs) that are constructed based on the message
passing scheme in which each node propagates and aggregates information, encoded by vectorized
messages, to and from its local neighborhood. While this framework has been immensely successful
in many applications, it lacks the ability to capture the multiscale and multiresolution structures that
are present in complex graphs (Rustamov & Guibas, 2013) (Chen et al., 2014) (Cheng et al., 2015)
(Xu et al., 2019).

Ying et al. (2018) proposed a multiresolution graph neural network that employs a differential pooling operator to coarsen the graph. While this approach is effective in some settings, it is based on
_soft assigment matrices, which means that (a) the sparsity of the graph is quickly lost in higher layers_
(b) the algorithm isn’t able to learn an actual hard clustering of the vertices. The latter is important
in applications such as learning molecular graphs, where clusters should ideally be interpretable as
concrete subunits of the graphs, e.g., functional groups.

In contrast, in this paper we propose an arhictecture called Multiresolution Graph Network (MGN),
and its generative cousin, Multiresolution Graph Variational Autoencoder (MGVAE), which explicitly learn a multilevel hard clustering of the vertices, leading to a true multiresolution hierarchy. In
the decoding stage, to “uncoarsen” the graph, MGVAE needs to generate local adjacency matrices,
which is inherently a second order task with respect to the action of permutations on the vertices,
hence MGVAE needs to leverage the recently developed framework of higher order permutation
equivariant message passing networks (Hy et al., 2018; Maron et al., 2019b).

Learning to generate graphs with deep generative models is a difficult problem because graphs are
combinatorial objects that typically have high order correlations between their discrete substructures (subgraphs) (You et al., 2018a) (Li et al., 2018) (Liao et al., 2019) (Liu et al., 2019) (Dai
et al., 2020). Graph-based molecular generation (Gmez-Bombarelli et al., 2018) (Simonovsky &


-----

Komodakis, 2018) (Cao & Kipf, 2018) (Jin et al., 2018) (Thiede et al., 2020) involves further challenges, including correctly recognizing chemical substructures, and importantly, ensuring that the
generated molecular graphs are chemically valid. MGN allows us to extend the existing model of
variational autoencoders (VAEs) with a hierarchy of latent distributions that can stochastically generate a graph in multiple resolution levels. Our experiments show that having a flexible clustering
procedure from MGN enables MGVAE to detect, reconstruct and finally generate important graph
substructures, especially chemical functional groups.

2 RELATED WORK

There have been significant advances in understanding the invariance and equivariance properties of
neural networks in general (Cohen & Welling, 2016a) (Cohen & Welling, 2016b), of graph neural
networks (Maron et al., 2019b), of neural networks learning on sets (Zaheer et al., 2017) (Serviansky
et al., 2020) (Maron et al., 2020), along with their expressive power on graphs (Maron et al., 2019c)
(Maron et al., 2019a). Our work is in line with group equivariant networks operating on graphs and
sets. Multiscale, multilevel, multiresolution and coarse-grained techniques have been widely applied
to graphs and discrete domains such as diffusion wavelets (Coifman & Maggioni, 2006); spectral
wavelets on graphs (Hammond et al., 2011); finding graph wavelets based on partitioning/clustering
(Rustamov & Guibas, 2013); graph clustering and finding balanced cuts on large graphs (Dhillon
et al., 2005) (Dhillon et al., 2007) (Chiang et al., 2012) (Si et al., 2014); and link prediction on social
networks (Shin et al., 2012). Prior to our work, some authors such as (Zhou et al., 2019) proposed
a multiscale generative model on graphs using GAN (Goodfellow et al., 2014), but the hierarchical
structure was built by heuristics algorithm, not learnable and not flexible to new data that is also
an existing limitation of the field. In general, our work exploits the powerful group equivariant
networks to encode a graph and to learn to form balanced partitions via back-propagation in a datadriven manner without using any heuristics as in the existing works.

In the field of deep generative models, it is generally recognized that introducing a hierarchy of
latents and adding stochasticity among latents leads to more powerful models capable of learning
more complicated distributions (Blei et al., 2003) (Ranganath et al., 2016) (Ingraham & Marks,
2017) (Klushyn et al., 2019) (Wu et al., 2020) (Vahdat & Kautz, 2020). Our work combines the
hierarchical variational autoencoder with learning to construct the hierarchy that results into a generative model able to generate graphs at many resolution levels.

3 MULTIRESOLUTION GRAPH NETWORK

3.1 CONSTRUCTION

An undirected weighted graph _,_ _,_ _,_ _v,_ _e_ with node set and edge set is represented
by an adjacency matrix N[∣V∣×∣V∣], where _ij_ 0 implies an edge between node vi and vj
with weight _ij (e.g.,_ _ij_ G = (V0, 1 in the case of unweighted graph); while node features areE A F F ) V E
represented by a matrix A ∈v R[∣V∣×][d][v], and edge features are represented by a tensor A > _e_ R[∣V∣×∣V∣×][d][e].
The second-order tensor representation of edge features is necessary for our higher-order message A A ∈{ }
passing networks described in the next section. Indeed, F ∈ _v can be encoded in the diagonal of F_ ∈ _e._
**Definition 1. A K-cluster partition of graph** _is a partition of the set of nodes_ _into K mutually_
F F
_exclusive clusters_ 1,.., _K. Each cluster corresponds to an induced subgraph_ _k_ _k_ _._
G V

**Definition 2. A coarsening of** _is a graph_ [˜] of K nodes defined by a K-cluster partition in which

V V G = G[V ]

_node ˜vk of_ [˜] corresponds to the induced subgraph _k. The weighted adjacency matrix_ [˜] N[K][×][K]

_of_ [˜] is G G

G ˜ 21 _i[,v]j_ [∈V]k G _if k_ _k[′],_ A ∈

G _kk_

′ _vi_ _k,vj_ _k_ _ij,_ _if k_ _k[′],_

[∑][v] [A]′ _[ij][,]_ =

_where the diagonal of_ [˜] denotes the number of edges inside each cluster, while the off-diagonalA = { ∈V ∈V

∑ A ≠

_denotes the number of edges between two clusters._
A

Fig. 3.1 shows an example of Defs. 1 and 2: a 3-cluster partition of the Aspirin C9H8O4 molecular
graph and its coarsening graph. Def. 3 defines the multiresolution of graph in a bottom-up manner
G


-----

2 3

1 3

1 1

V V

2
V

6


Figure 1: Aspirin C9H8O4, its 3-cluster partition and the corresponding coarsen graph

in which the bottom level is the highest resolution (e.g., itself) while the top level is the lowest
resolution (e.g., is coarsened into a single node).
G

**Definition 3. An L-level coarsening of a graph** _is a series of L graphs_ _,..,_ _in which_

G

_1._ _is_ _itself._ G G[(][1][)] G[(][L][)]

_2. For G[(][L] 1[)]_ Gℓ _L_ 1, _is a coarsening graph of_ _as defined in Def. 2. The number_
_of nodes in_ _is equal to the number of clusters in_ _._
≤ ≤ − G[(][ℓ][)] G[(][ℓ][+][1][)]

_3. The top level coarsening G[(][ℓ][)]_ _is a graph consisting of a single node, and the corresponding G[(][ℓ][+][1][)]_
_adjacency matrix_ N[1][×][1].
G[(][1][)]

**Definition 4. An L-level Multiresolution Graph Network (MGN) of a graph** _consists of L_ 1 tu
A[(][1][)] ∈

_ples of five network components_ **_c[(][ℓ][)],_** **_e[(]local[ℓ][)]_** _[,]_ **_[d][(]local[ℓ][)]_** _[,]_ **_[d][(]global[ℓ][)]_** _[,]_ **_[p][(][ℓ][)][)}]ℓ[L]_** 2[. The][ ℓ][-th tuple encodes][ G][(][ℓ][)]
G −

_and transforms it into a lower resolution graph_ _in the higher level. Each of these network=_

{(

_components has a separate set of learnable parameters_ **_θ1[(][ℓ][)][,]_** **_[θ]2[(][ℓ][)][,]_** **_[θ]3[(][ℓ][)][,]_** **_[θ]4[(][ℓ][)][,]_** **_[θ]5[(][ℓ][)]_**
_ity, we collectively denote the learnable parameters as G[(][ℓ][−][1][)]_ **_θ, and drop the superscript. The network_**
_components are defined as follows:_ ( [)][. For simplic-]

_1. Clustering_ _procedure_ **_c_** ; **_θ_** _,_ _which_ _partitions_ _graph_ _into_ _K_ _clusters_
1 _[,...,]_ [V]K[(][ℓ][)][. Each cluster is an induced subgraph][ G]k[(][ℓ][)] _of_ _with adjacency matrix_
_k_ _[.]_ (G[(][ℓ][)] ) G[(][ℓ][)]
V [(][ℓ][)] G[(][ℓ][)]

_2. Local encoderA[(][ℓ][)]_ **_elocal_** _k_ [;] **_[θ][)][, which is a permutation equivariant (see Defs. 5, 6) graph]_**
_neural network that takes as input the subgraph_ _k_ _[, and outputs a set of node latents]_
_k_ _represented as a matrix of size(G[(][ℓ][)]_ _k_
G[(][ℓ][)]

_3. Local decoderZ_ [(][ℓ][)] **_dlocal_** _k_ [;] **_[θ][)][, which is a permutation equivariant neural network that tries] ∣V_** [(][ℓ][)][∣×][ d][z][.]
_to reconstruct the subgraph adjacency matrix_ _k_ _for each cluster from the local encoder’s_
_output latents._ (Z [(][ℓ][)]
A[(][ℓ][)]

_4. (Optional) Global decoder dglobal_ ; **_θ_** _, which is a permutation equivariant neural net-_
_work that reconstructs the full adjacency matrix_ _from all the node latents of K clusters_
_k_ _k_ _represented as a matrix of size(Z_ [(][ℓ][)] ) _dz._
A[(][ℓ][)]

_5. Pooling networkZ_ [(][ℓ][)] = ⊕ Z [(][ℓ][)] **_p_** _k_ [;] **_[θ][)][, which is a permutation invariant (see Defs. 5, 6) neural net-] ∣V_** [(][ℓ][)]∣×
_work that takes the set of node latents_ _k_ _and outputs a single cluster latent ˜zk[(][ℓ][)]_ _dz._
_The coarsening graph(Z_ [(][ℓ][)] _has adjacency matrix_ _built as in Def. 2, and the cor-_
_responding node features_ _k ˜ Zzk[(][ℓ][(][)][ℓ][)]represented as a matrix of size K_ _dz._ ∈
G[(][ℓ][−][1][)] A[(][ℓ][−][1][)]

Algorithmically, MGN works in a bottom-up manner as a tree-like hierarchy starting from the high- Z [(][ℓ][−][1][)] = ⊕ ×
est resolution graph, going to the lowest resolution (see Fig. 3.1). Iteratively, at ℓ-th level,
MGN partitions the current graph into K clusters by running the clustering procedure c[(][ℓ][)]. Then,
the local encoder e[(]local G[ℓ][)] [(][L][and local decoder][)] **_[ d]global[(][ℓ][)]_** [operate on each of the] G[(][1][)] _[ K][ subgraphs separately, and]_
can be executed in parallel. This encoder/decoder pair is responsible for capturing the local structures. Finally, the pooling network p[(][ℓ][)] shrinks each cluster into a single node of the next level.


-----

0 1 0 0 0 1
1 0 1 0 0 0

2 ⎛⎜00 10 01 10 01 00⎞⎟ 21 16 01
A[(][3][)] = ⎜⎜⎜⎜⎜01 00 00 10 01 10⎟⎟⎟⎟⎟ G[(][1][)] A[(][2][)] = ⎛⎜0 1 3⎞⎟

0⎝ 1 1 ⎠ ⎝ ⎠ 0 1 0 0

1 0 1 1

A1[(][3][)] = ⎛⎜11 00 00⎞⎟ V1[(][2][)] V2[(][2][)] V3[(][2][)] A3[(][3][)] = ⎛⎜⎜00 11 00 00
⎝ ⎠ ⎝

_v0_ _v1_ _v2_ _v3_ _v4_ _v5_ _v6_ _v7_ _v8_ _v9_ _v10_ _v11_ _v12_


Figure 2: Hierarchy of 3-level Multiresolution Graph Network on Aspirin molecular graph

Optionally, the global decoder d[(]global[ℓ][)] [makes sure that the whole set of node latents][ Z] [(][ℓ][)][ is able to]
capture the inter-connection between clusters.

In terms of time and space complexity, MGN is more efficient than existing methods in the field.
The cost of global decoding the highest resolution graph is proportional to . For example, while
the encoder can exploit the sparsity of the graph and has complexity, a simple dot-product
decoder dglobal sigmoid has both time and space complexity of ∣V∣[2] which is infeasible for large graphs. In contrast, the cost of running K local dot-product decoders is O(∣E∣) _K_,
which is approximately(Z) = _K times more efficient.(ZZ_ _[T]_ ) O(∣V∣[2])
O(∣V∣[2]/ )

3.2 HIGHER ORDER MESSAGE PASSING

In this paper we consider permutation symmetry, i.e., symmetry to the action of the symmetric
group, Sn. An element σ Sn is a permutation of order n, or a bijective map from 1,...,n to
1,...,n . The action of Sn on an adjacency matrix R[n][×][n] and on a latent matrix R[n][×][d][z]
are ∈ { }
{ }σ _i1,i2_ _σ−1_ _i1_ _,σ−1_ _i2_ _,_ _σ_ A ∈i,j _σ−1_ _i_ _,j,_ _σ_ Sn Z ∈.
Here, the adjacency matrix ( is a second order tensor with a single feature channel, while the latent) ( ) ( )
matrix is a first order tensor with[ ⋅A] = A _dz feature channels. In general, the action of[_ ⋅Z] = Z Sn on a ∈ _k-th order_
tensor R[n][k][×][d] (the last index denotes the feature channels) is defined similarly as: A
Z
X ∈ _σ_ _i1,..,ik,j_ _σ_ 1 _i1_ _,..,σ_ 1 _ik_ _,j,_ _σ_ Sn.

− ( ) − ( )

Network components of MGN (as defined in Sec. 3.1) at each resolution level must be either[ ⋅X] = X ∈ _equiv-_
_ariant, or invariant with respect to the permutation action on the node order of_ . Formally, we
define these properties in Def. 5.
G[(][ℓ][)]
**Definition** **5.** _An_ Sn-equivariant _(or_ _permutation_ _equivariant)_ _function_ _is_ _a_ _function_
_f_ R[n][k][×][d] R[n][k][′] [×][d][′] _that satisfies f_ _σ_ _σ_ _f_ _for all σ_ Sn and R[n][k][×][d]. Similarly,
_we say that f is Sn-invariant (or permutation invariant) if and only if f_ _σ_ _f_ _._
∶ → ( ⋅X) = ⋅ (X) ∈ X ∈

**Definition 6. An Sn-equivariant network is a function f** R[n][k][×][d] R[n][k][′]([×][d][′] ⋅X) =defined as a composition(X)
_of Sn-equivariant linear functions f1,..,fT and Sn-equivariant nonlinearities γ1,..,γT :_
∶ →

_f_ _γT_ _fT_ _.._ _γ1_ _f1._

_On the another hand, an Sn-invariant network is a function =_ ○ ○ ○ ○ f R[n][k][×][d] R defined as a composition
_of an Sn-equivariant network f_ _and an Sn-invariant function on top of it, e.g., f_ _f_ _f_ _._
∶ →

[′] [′′] [′]

In order to build higher order equivariant networks, we revisit some basic tensor operations: the = ○
tensor product A _B and tensor contraction A_ _x1,..,xp (details and definitions are Sec. A). It can be_
shown that these tensor operations respect permutation equivariance (Kondor et al., 2018). Based
↓
on these tensor contractions and Def. 5, we can construct the second-order ⊗ Sn-equivariant networks
as in Def. 6 (see Example 1): f _γ_ _T_ _.._ _γ_ 1. The second-order networks are particularly essential for us to extend the original variational autoencoder (VAE) (Kingma & Welling,
= ○M ○ ○ ○M


-----

2014) model that approximates the posterior distribution by an isotropic Gaussian distribution with
a diagonal covariance matrix and uses a fixed prior distribution 0, 1 . In constrast, we generalize
by modeling the posterior by _µ,_ Σ in which Σ is a full covariance matrix, and we learn an adaptive parameterized prior _µ,ˆ_ Σ[ˆ] instead of a fixed one. Only the second-order encoders can output N( )
a permutation equivariant full covariance matrix, while lower-order networks such as MPNNs are N( )
unable to. See Sec. 4.2, B and C for details. N( )
**Example 1. The second order message passing has the message** 0 R[∣V∣×∣V∣×(][d][v][+][d][e][)] _initialized_
_by promoting the node features_ _v to a second order tensor (e.g., we treat node features as self-loop_
_edge features), and concatenating with the edge features_ _e. Iteratively, H_ ∈
F

_t_ _γ_ _t_ _,_ _t_ _t_ F _t_ 1 _i,j_ _,_
_i,j_
− ↓

_where_ _t_ 1 results in a fourth order tensor whileH = (M ) M = W [⊕(A ⊗Hi,j contracts it down to a second order) ]
_tensor along the i-th and j-th dimensions,_ _denotes concatenation along the feature channels, and_
−
_t denotes a multilayer perceptron on the feature channels. We remark that the popular MPNNs A ⊗H_ ↓
_(Gilmer et al., 2017b) is a lower-order one and a special case in which ⊕_ _t_ _t_ 1 _t_ 1
Wwhere _ii_ _j_ _ij is the diagonal matrix of node degrees. The message_ _T of the last iteration is_
− −
_still second order, so we contract it down to the first order latent_ _i_ MT _i. = D[−][1]AH_ W
D = ∑ A H
↓

3.3 LEARNING TO CLUSTER Z = ⊕ H

**Definition 7. A clustering of n objects into k clusters is a mapping π** 1,..,n 1,..,k _in which_
_π_ _i_ _j if the i-th object is assigned to the j-th cluster. The inverse mapping π[−][1]_ _j_ _i_ 1,n
_π_ _i_ _j_ _gives the set of all objects assigned to the j-th cluster. The clustering is represented by an ∶{_ } →{ }
_assignment matrix(_ ) = Π 0, 1 _such that Πi,π_ _i_ 1. ( ) = { ∈[ ] ∶
( ) = }
**Definition 8. The action of Sn on a clustering π(** _of)_ _n objects into k clusters and its corresponding_
_assignment matrix Π are ∈{_ }[n][×][k] =

_σ_ _π_ _i_ _π_ _σ[−][1]_ _i_ _,_ _σ_ Π _i,j_ Πσ−1 _i_ _,j,_ _σ_ Sn.

**Definition 9. Let** _be a neural network that takes as input a graph(_ ) _of n nodes, and outputs a_
_clustering π of k[ clusters. ⋅_ ]( ) = _is said to be equivariant if and only if(_ ( )) [ ⋅ ] = _σ_ ∈ _σ_ _for all_
_σ_ Sn. N G
N N( ⋅G) = ⋅N(G)

From Def. 9, intuitively the assignement matrix ∈ Π still represents the same clustering if we permute
its rows. The learnable clustering procedure c ; **_θ_** is built as follows:

1. A graph neural network parameterized by θ encodes graph into a first order tensor of K

(G[(][ℓ][)] )

feature channels ˜p[(][ℓ][)] R[∣V] [(][ℓ][)][∣×][K].

2. The clustering assignment is determined by a row-wise maximum pooling operation: G[(][ℓ][)]

∈

_π[(][ℓ][)]_ _i_ arg max _p˜[(]i,k[ℓ][)]_ (1)
_k_ 1,K

that is an equivariant clustering in the sense of Def. 9.( ) = ∈[ ]

A composition of an equivariant function (e.g., graph net) and an equivariant function (e.g., maximum pooling given in Eq. 1) is still an equivariant function with respect to the node permutation.
Thus, the learnable clustering procedure c ; **_θ_** is permutation equivariant.

In practice, in order to make the clustering procedure differentiable for backpropagation, we replace

(G[(][ℓ][)] )

the maximum pooling in Eq. 1 by sampling from a categorical distribution. Let π[(][ℓ][)] _i_ be a categorical variable with class probabilities p[(]i,[ℓ]1[)][,..,p]i,K[(][ℓ][)] [computed as softmax from][ ˜]p[(]i,[ℓ][)]
( )
trick (Gumbel, 1954)(Maddison et al., 2014)(Jang et al., 2017) provides a simple and efficient way
∶ [. The Gumbel-max]
to draw samples π[(][ℓ][)] _i_ :
( ) Π[(]i[ℓ][)] one-hot arg max _gi,k_ log p[(]i,k[ℓ][)][])][,]

_k_ 1,K

where gi,1,..,gi,K are i.i.d samples drawn from Gumbel= ( ∈[ ] [ +0, 1 . Given the clustering assignment
matrix Π[(][ℓ][)], the coarsened adjacency matrix (see Defs. 1 and 2) can be constructed as
Π[(][ℓ][)][T] Π[(][ℓ][)]. ( )
A[(][ℓ][−][1][)]
A[(][ℓ][)]


-----

It is desirable to have a balanced K-cluster partition in which clusters 1 _[,..,]_ [V]K[(][ℓ][)] [have similar sizes]
that are close to _K. The local encoders tend to generalize better for same-size subgraphs. We_
want the distribution of nodes into clusters to be close to the uniform distribution. We enforce V [(][ℓ][)]
the clustering procedure to produce a balanced cut by minimizing the following Kullback–Leibler ∣V [(][ℓ][)]∣/
divergence:

_K_
_KL_ _P_ _Q_ _P_ _k_ log _[P]_ [(][k][)] where _P_ 1 _K_ _Q_ [1] (2)
_k_ 1 _Q_ _k_ _K [,..,][ 1]K_

[∣] [∣]

D ( ∣∣ ) = ∑= ( ) = ( [∣V] [(][ℓ][)] = (

The whole construction of MGN is( equivariant) with respect to node permutations of∣V [(][ℓ][)]∣ _[,..,][ ∣V]∣V_ [(][(][ℓ][ℓ][)][)]∣ [)][,] . In the case[)][.]
of molecular property prediction, we want MGN to learn to predict a real value y R for each graph
while learning to find a balanced cut in each resolution to construct a hierarchical structure of G
latents and coarsen graphs. The total loss function is ∈
G _L_ _L_

2
MGN _,y_ _f_ _R_ _y_ 2 _λ[(][ℓ][)]_ KL _P_ _Q[(][ℓ][)]_ _,_ (3)
_ℓ_ 1 _ℓ_ 1

where f is a multilayer perceptron,L (G ) = ∣∣ ( ⊕= denotes the vector concatenation,(Z [(][ℓ][)])) − ∣∣ [+] ∑= D ( _R[(][ℓ] is a readout function that[)]∣∣_ )
produces a permutation invariant vector of size d given the latent at the ℓ-th resolution,
⊕
_λ[(][ℓ][)]_ R is a hyperparamter, and KL _P_ _Q[(][ℓ][)]_ is the balanced-cut loss as defined in Eq. 2.
Z [∣V] [(][ℓ][)][∣×][d]
∈ D ( [(][ℓ][)]∣∣ )

4 HIERARCHICAL GENERATIVE MODEL

In this section, we introduce our hierarchical generative model for multiresolution graph generation
based on variational principles.

4.1 BACKGROUND ON GRAPH VARIATIONAL AUTOENCODER

Suppose that we have input data consisting of m graphs (data points) 1,.., _m_ . The standard variational autoencoders (VAEs), introduced by Kingma & Welling (2014) have the following
generation process, in which each data graph _i for i_ 1, 2,..,m is generated independently: G = {G G }

1. Generate the latent variables 1,.., G _m_, where each ∈{ }i R[∣V][i][∣×][d][z] is drawn i.i.d. from a
prior distribution p0 (e.g., standard Normal distribution 0, 1 ).

2. Generate the data graph _i_ Z = {Zpθ _i_ _i_ from the model conditional distributionZ } Z ∈ _pθ._

N( )

We want to optimize θ to maximize the likelihood pθ _pθ_ _pθ_ _d_ . However, this re
G ∼ (G ∣Z )

quires computing the posterior distribution pθ _i_ 1 _[p][θ][(G][i][∣Z][i][)][, which is usually intractable.]_
Instead, VAEs apply the variational principle, proposed by Wainwright & Jordan (2005), to approxi-(G) = ∫ (Z) (G∣Z) Z
=
mate the posterior distribution as qφ (G∣Z) = ∏i 1 _[q][φ][(Z][i][∣G][m][i][)][ via amortized inference and maximize]_
the evidence lower bound (ELBO) that is a lower bound of the likelihood:
=
ELBO _φ,θ_ Eqφ (Z∣G) = ∏log pθ _[m]_ KL _qφ_ _p0_

_m_ (Z∣G) (4)

L ( ) = Eqφ [ _i_ _i_ log(G∣Z)] −D pθ _i_ _i_ ( (Z∣G)∣∣KL _qφ_ _i(Z))i_ _p0_ _i_ _._

_i_ 1
(Z ∣G )

The probabilistic encoder= q∑=φ [, the approximation to the posterior of the generative model[ (G ∣Z )] −D ( (Z ∣G )∣∣ (Z ))]
_pθ_ _,_, is modeled using equivariant graph neural networks (see Example 1) as follows. Assume the prior over the latent variables to be the centered isotropic multivariate Gaussian(Z∣G) _pθ_
(G ;0Z),I . We let qφ _i_ _i_ be a multivariate Gaussian with a diagonal covariance structure:
(Z) =
N(Z ) (Z ∣G ) log qφ _i_ _i_ log _i;_ _µi,σi[2][I][)][,]_ (5)

where µi,σi R[∣V][i][∣×][d][z] are the mean and standard deviation of the approximate posterior output by

(Z ∣G ) = N(Z

two equivariant graph encoders. We sample from the posterior qφ by using the reparameterization
trick: _i_ _µ ∈i_ _σi_ _ϵ, where ϵ_ 0,I and is the element-wise product.

On the another hand, the probabilistic decoder pθ _i_ _i_ defines a conditional distribution over the
entries of the adjacency matrix Z = + ⊙ ∼N(i: pθ )i _i_ ⊙ _u,v_ _i2_ _[p][θ][(A][iuv][ =][ 1][∣Z][iu][,]_ [Z][iv][)][. For example,]
Kipf & Welling (2016) suggests a simple dot-product decoder that is trivially equivariant:(G (∣Z ))∈V _pθ_ _iuv_
1 _iu,_ _iv_ _γ_ _iu[Z][iv][)][, where] A[ γ][ denotes the sigmoid function.](G_ ∣Z ) = ∏
(A =
∣Z Z ) = (Z _[T]_


-----

4.2 MULTIRESOLUTION VAES

Based on the construction of multiresolution graph network (see Sec. 3.1), the latent variables are partitioned into disjoint groups, _i_ _i_ _,_ _i_ _,..,_ _i_ where _i_ _i_ _k_

R[∣[V]i[(][ℓ][)] _k_ _dz_ _k is the set of latents at the ℓ-th resolution level in which the graph_ _i_ is partitioned
into a number of clusters] ∣× _i_ _k._ Z = {Z [(][1][)] Z [(][2][)] Z [(][L][)]} Z [(][ℓ][)] = {[Z [(][ℓ][)]] ∈
} G[(][ℓ][)]

In the area of normalzing flows (NFs), Wu et al. (2020) has shown that stochasticity (e.g., a chain

[G[(][ℓ][)]]

of stochastic sampling blocks) overcomes expressivity limitations of NFs. In general, our MGVAE
is a stochastic version of the deterministic MGN such that stochastic sampling is applied at each
resolution level in a bottom-up manner. The prior (Eq. 6) and the approximate posterior (Eq. 7) are
represented by

_L_ _L_
_p_ _i_ _p_ _i_ _p_ _i_ _k_ _,_ (6)
_ℓ_ 1 _ℓ_ 1 _k_

1

(Z ) = ∏= (Z [(][ℓ][)]) = ∏= ∏ ([Z [(][ℓ][)]] )

_qφ_ _i_ _i_ _qφ_ _i_ _i_ _qφ_ _i_ _i_ _,_ _i_ _,_ (7)
_ℓ_ _L_ 1

in which each conditional in the approximate posterior are in the form of factorial Normal distribu-(Z ∣G ) = (Z [(][L][)]∣G[(][L][)]) =∏− (Z [(][ℓ][)]∣Z [(][ℓ][+][1][)] G[(][ℓ][)])
tions, in particular

_qφ_ _i_ _i_ _,_ _i_ _qφ_ _i_ _k_ _i_ _,_ _i_ _k_ _,_
_k_

where each probabilistic encoder(Z [(][ℓ][)]∣Z [(][ℓ][+] q[1][)]φ G[(][ℓ]i[)]) = ∏k _i_ ([Z, [(][ℓ]i[)]] ∣Zk [(] operates on a subgraph[ℓ][+][1][)] [G[(][ℓ][)]] ) _i_ _k as_
follows:
([Z [(][ℓ][)]] ∣Z [(][ℓ][+][1][)] [G[(][ℓ][)]] ) [G[(][ℓ][)]]

-  The pooling network p[(][ℓ][+][1][)] shrinks the latent _i_ into the node features of _i_ as in
the construction of MGN (see Def. 4).
Z [(][ℓ][+][1][)] G[(][ℓ][)]

-  The local (deterministic) graph encoder d[(]local[ℓ][)] [encodes each subgraph][ [G]i[(][ℓ][)] _k into a mean_
vector and a diagonal covariance matrix (see Eq. 5). A second order encoder can produce a
positive semidefinite non-diagonal covariance matrix, that can be interpreted as a Gaussian]
Markov Random Fields (details in Sec. B). The new subgraph latent _i_ _k is sampled by_
the reparameterization trick.
[Z [(][ℓ][)]]

The prior can be either the isotropic Gaussian 0, 1 as in standard VAEs, or be implemented
as a parameterized Gaussian _µ,ˆ_ Σ[ˆ] where ˆµ and Σ[ˆ] are learnable equivariant functions (details in
Sec. C). The reparameterization trick for conventional N( ) 0, 1 prior is the same as in Sec. 4.1, while
the new one for the generalized and learnable prior N( ) _µ,ˆ_ Σ[ˆ] is given in Sec. B. On the another hand,
the probabilistic decoder pθ _i_ _,..,_ _i_ _i_ _,..,_ _i N( defines a conditional distribution over all)_
subgraph adjacencies at each resolution level: N( )
(G[(][1][)] G[(][L][)]∣Z [(][1][)] Z [(][L][)])

_pθ_ _i_ _,..,_ _i_ _i_ _,..,_ _i_ _pθ_ _i_ _i_ _pθ_ _i_ _i_ _k_ _._
_ℓ_ _ℓ_ _k_

Extending from Eq. 4, we write our multiresolution variational lower bound(G[(][1][)] G[(][L][)]∣Z [(][1][)] Z [(][L][)]) = ∏ (G[(][ℓ][)]∣Z [(][ℓ][)]) = ∏ ∏ ([A[(][ℓ][)][]][k][∣[Z]MGVAE[(][ℓ][)]] ) _φ,θ_ on
log p compactly as
L ( )

where the first term denotes the reconstruction loss (e.g.,LMGVAE(G) (φ,θ) = ∑i ∑ℓ [Eqφ(Zi(ℓ)∣Gi[(][ℓ][)])[[][log][ p][θ][(G]i[(][ℓ][)]∣Zi[(][ℓ][)])] −DiKL(qφ(Zi _i[(][ℓ][)]∣Gi[(][ℓ][)])∣∣p0i(Zi[(]is[ℓ][)]))]i itself,, (8)_
_i_ is the adjacency produced by MGN at level ℓ, and [ˆ]i are the reconstructed ones by the de
coders); and the second term is indeed KL _µ[(]i[ℓ][)][,]_ [Σ]i[(][ℓ][)] ∣∣A[(][ℓ]µ[)] [(]−[ℓ][)]A,[ˆ]Σ[ˆ][(][ℓ][(][)][ℓ][∣∣][)] [where] where[ A] µ[(][L][(]i[)][ℓ][)] GR[∣V]i[(][ℓ][)] _d_

A[(][ℓ][<][L][)] A[(][ℓ][)]

and Σ[(]i[ℓ][)] R[∣V]i[(][ℓ][)] _i_ _d are the mean and covariance tensors produced by the ℓ-th encoder for∣×_
graph _i, while ˆµ∣×∣V[(][ℓ][)]and[(][ℓ][)]∣×Σ[ˆ]_ [(][ℓ][)] are learnable ones in an equivariant manner as in Sec. C. In general, D (N( [)∣∣N(][ˆ] )) ∈
the overall optimization is given as follows:∈
G

min MGVAE _φ,θ;_ **_µˆ[(][ℓ][)],_** **Σ[ˆ]** [(][ℓ][)] _ℓ_ _λ[(][ℓ][)]_ KL _Pi_ _Q[(]i[ℓ][)]_ (9)
_φ,θ,_ **_µˆ[(][ℓ][)],Σ[ˆ]_** [(][ℓ][)] _ℓ_ _i,ℓ_

where φ denotes all learnable parameters of the encoders,{ } L ( { } ) + ∑ θ denotes all learnable parameters of theD ( [(][ℓ][)]∣∣ [)][,]
decoders, and KL _Pi_ _Q[(]i[ℓ][)]_
D ( [(][ℓ][)]∣∣ [)][ is the balanced-cut loss for graph][ G][i][ at level][ ℓ] [as defined in Sec. 3.3.]


-----

Figure 3: MGVAE generates molecules on QM9 (4 on the left) and ZINC (the rest) equivariantly.
There are many more examples of generated molecules in Sec. D.2. Both equivariant MGVAE and
autoregressive MGN generate high-quality molecules with complicated structures such as rings.

|Dataset|Method|Training size|Input features|Validity|Novelty|Uniqueness|
|---|---|---|---|---|---|---|
|QM9|GraphVAE|100K ∼|Graph|61.00%|85.00%|40.90%|
||CGVAE|||100%|94.35%|98.57%|
||MolGAN|||98.1%|94.2%|10.4%|
||Autoregressive MGN|10K||100%|95.01%|97.44%|
||All-at-once MGVAE|||100%|100%|95.16%|
|ZINC|GraphVAE|200K ∼|Graph|14.00%|100%|31.60%|
||CGVAE|||100%|100%|99.82%|
||JT-VAE|||100%|-|-|
||Autoregressive MGN|1K||100%|99.89%|99.69%|
||All-at-once MGVAE|10K|Chemical|99.92%|100%|99.34%|


**Dataset** **Method** **Training size Input features Validity Novelty Uniqueness**

GraphVAE 61.00% 85.00% 40.90%

CGVAE 100K 100% 94.35% 98.57%

QM9 MolGAN Graph 98.1% 94.2% 10.4%

**Autoregressive MGN** ∼ 100% 95.01% 97.44%
10K

**All-at-once MGVAE** 100% 100% 95.16%

GraphVAE 14.00% 100% 31.60%

CGVAE 200K 100% 100% 99.82%

ZINC Graph

JT-VAE 100% - 
**Autoregressive MGN** ∼ 1K 100% 99.89% 99.69%

**All-at-once MGVAE** 10K Chemical 99.92% 100% 99.34%


Table 1: Molecular graph generation results. GraphVAE results are taken from (Liu et al., 2018).

5 EXPERIMENTS

Many more experimental results and details are presented in the Sec. D of the Appendix.

5.1 MOLECULAR GRAPH GENERATION

We examine the generative power of MGN and MGVAE in the challenging task of molecule generation, in which the graphs are highly structured. We demonstrate that MGVAE is the first hierarchical
graph VAE model generating graphs in a permutation-equivariant manner that is competitive against
autoregressive results. We train on two datasets that are standard in the field:

1. QM9 (Ruddigkeit et al., 2012) (Ramakrishnan et al., 2014): contains 134K organic
molecules with up to nine atoms (C, H, O, N, and F) out of the GDB-17 universe of
molecules.

2. ZINC (Sterling & Irwin, 2015): contains 250K purchasable drug-like chemical compounds
with up to twenty-three heavy atoms.

We only use the graph features as the input, including the adjacency matrix, the one-hot vector
of atom types (e.g., carbon, hydrogen, etc.) and the bond types (single bond, double bond, etc.)
without any further domain knowledge from chemistry or physics. First, we train autoencoding
task of reconstructing the adjacency matrix and node features. We use a learnable equivariant prior
(see Sec. C) instead of the conventional 0, 1 . Then, we generate 5, 000 different samples from
the prior, and decode each sample into a generated graph (see Fig. 3). We implement our graph
construction (decoding) in two approaches: N( )

1. All-at-once: We reconstruct the whole adjacency matrix by running the probabilistic decoder (see Sec. 4). MGVAE enables us to generate a graph at any given resolution level ℓ.
In this particular case, we select the highest resolution ℓ _L. This approach of decoding_
preserves equivariance, but is harder to train. On ZINC, we extract several chemical/atomic
features from RDKit as the input for the encoders to reach a good convergence in training.=

2. Autoregressive: The graph is constructed iteratively by adding one edge in each iteration,
similarly to (Liu et al., 2018). But this approach does not respect permutation equivariance.

In our setting for small molecules, L 3 and K 2[ℓ][−][1] for the ℓ-th level. We compare our methods
with other graph-based generative models including GraphVAE (Simonovsky & Komodakis, 2018),
CGVAE (Liu et al., 2018), MolGAN (Cao & Kipf, 2018), and JT-VAE (Jin et al., 2018). We evaluate = =
the quality of generated molecules in three metrics: (i) validity, (ii) novelty and (iii) uniqueness as


-----

**COMMUNITY-SMALL** **EGO-SMALL**

**MODEL** DEGREE CLUSTER ORBIT DEGREE CLUSTER ORBIT

GRAPHVAE 0.35 0.98 0.54 0.13 0.17 0.05
DEEPGMG 0.22 0.95 0.4 0.04 0.10 0.02
GRAPHRNN 0.08 0.12 0.04 0.09 0.22 0.003
GNF 0.20 0.20 0.11 0.03 0.10 0.001
GRAPHAF 0.06 0.10 0.015 0.04 0.04 0.008

**MGVAE** **0.002** **0.01** **0.01** **1.74e-05** **0.0006** **6.53e-05**

Table 2: Graph generation results depicting MMD for various graph statistics between the test set
and generated graphs. MGVAE outperforms all competing methods.

the percentage of the generated molecules that are chemically valid, different from all molecules
in the training set, and not redundant, respectively. Because of high complexity, we only train on a
small random subset of examples while all other methods are trained on the full datasets. Our models
are equivalent with the state-of-the-art, even with a limited training set (see Table 1). Admittedly,
molecule generation is a somewhat subject task that can only be evaluated with objective numerical
measures up to a certain point. Qualitatively, however the molecules that MGVAE generates are as
good as the state of the art, in some cases better in terms of producing several high-quality drug-like
molecules with complicated functional groups and structures. Many further samples generated by
MGVAE and their analysis can be found in the Appendix.

5.2 GENERAL GRAPH GENERATION BY MGVAE

We further examine the expressive power of hierarchical latent structure of MGVAE in the task of
general graph generation. We choose two datasets from GraphRNN paper (You et al., 2018a):

1. Community-small: A synthetic dataset of 100 2-community graphs where 12 _V_ 20.
2. Ego-small: 200 3-hop ego networks extracted from the Citeseer network (Sen et al., 2008)
where 4 _V_ 18. ≤∣ ∣≤

The datasets are generated by the scripts from the GraphRNN codebase (You et al., 2018b). We keep

≤∣ ∣≤

80% of the data for training and the rest for testing. We evaluate our generated graphs by computing
Maximum Mean Discrepancy (MMD) distance between the distributions of graph statistics on the
test set and the generated set as proposed by (You et al., 2018a). The graph statistics are node
degrees, clustering coefficients, and orbit counts. As suggested by (Liu et al., 2019), we execute 15
runs with different random seeds, and we generate 1,024 graphs for each run, then we average the
results over 15 runs. We compare MGVAE against GraphVAE (Simonovsky & Komodakis, 2018),
DeepGMG (Li et al., 2018), GraphRNN (You et al., 2018a), GNF (Liu et al., 2019), and GraphAF
(Shi et al., 2020). The baselines are taken from GNF paper (Liu et al., 2019) and GraphAF paper (Shi
et al., 2020). In our setting of (all-at-once) MGVAE, we implement only L 2 levels of resolution
and K 2[ℓ] clusters for each level. Our encoders have 10 layers of message passing. Instead of
using a high order equivariant network as the global decoder for the bottom resolution, we only =
implement a simple fully connected network that maps the latent = R[∣V∣×][d][z] into an adjacency
matrix of size . For the ego dataset in particular, we implement the learnable equivariant
prior as in Sec. B and Sec.C. Table 2 includes our quantitative results in comparison with other Z [(][L][)] ∈
methods. MGVAE outperforms all competing methods. Figs. 10 11 show some generated examples ∣V∣× ∣V∣
and training examples on the 2-community and ego datasets.

6 CONCLUSION

We introduced MGVAE built upon MGN, the first generative model to learn and generate graphs in
a multiresolution and equivariant manner. The key idea of MGVAE is learning to construct a series
of coarsened graphs along with a hierarchy of latent distributions in the encoding process while
learning to decode each latent into the corresponding coarsened graph at every resolution level.
MGVAE achieves state-of-the-art results from link prediction to molecule and graph generation,
suggesting that accounting for the multiscale structure of graphs is a promising way to make graph
neural networks even more powerful.


-----

REFERENCES

Peter Battaglia, Razvan Pascanu, Matthew Lai, Danilo Jimenez Rezende, and Koray Kavukcuoglu.
Interaction networks for learning about objects, relations and physics. In Advances in Neural
_Information Processing Systems, volume 29. Curran Associates, Inc., 2016._

David M. Blei, Andrew Y. Ng, and Michael I. Jordan. Latent dirichlet allocation. J. Mach. Learn.
_Res., 3(null):9931022, March 2003. ISSN 1532-4435._

Nicola De Cao and Thomas Kipf. Molgan: An implicit generative model for small molecular graphs,
2018.

Xu Chen, Xiuyuan Cheng, and Stephane Mallat. Unsupervised deep haar scattering on graphs. In
_Advances in Neural Information Processing Systems, volume 27. Curran Associates, Inc., 2014._

Xiuyuan Cheng, Xu Chen, and St´ephane Mallat. Deep haar scattering networks. _CoRR,_
abs/1509.09187, 2015.

Kai-Yang Chiang, Joyce Jiyoung Whang, and Inderjit S. Dhillon. Scalable clustering of signed
networks using balance normalized cut. In Proceedings of the 21st ACM International Conference
_on Information and Knowledge Management, CIKM ’12, pp. 615624, New York, NY, USA, 2012._
Association for Computing Machinery. ISBN 9781450311564. doi: 10.1145/2396761.2396841.

Taco S. Cohen and Max Welling. Group equivariant convolutional networks. Proceedings of The
_33rd International Conference on Machine Learning, 48:2990–2999, 2016a._

Taco S. Cohen and Max Welling. Steerable cnns. ICLR’17, 2016b.

Ronald R. Coifman and Mauro Maggioni. Diffusion wavelets. Applied and Computational Har_monic Analysis, 21(1):53–94, 2006._ ISSN 1063-5203. doi: https://doi.org/10.1016/j.acha.
2006.04.004. [URL https://www.sciencedirect.com/science/article/pii/](https://www.sciencedirect.com/science/article/pii/S106352030600056X)
[S106352030600056X. Special Issue: Diffusion Maps and Wavelets.](https://www.sciencedirect.com/science/article/pii/S106352030600056X)

Hanjun Dai, Azade Nazi, Yujia Li, Bo Dai, and Dale Schuurmans. Scalable deep generative modeling for sparse graphs. In Proceedings of the 37th International Conference on Machine Learning,
volume 119 of Proceedings of Machine Learning Research, pp. 2302–2312. PMLR, 13–18 Jul
2020.

Inderjit Dhillon, Yuqiang Guan, and Brian Kulis. A fast kernel-based multilevel algorithm for graph
clustering. In Proceedings of the Eleventh ACM SIGKDD International Conference on Knowledge
_Discovery in Data Mining, KDD ’05, pp. 629634, New York, NY, USA, 2005. Association for_
Computing Machinery. ISBN 159593135X. doi: 10.1145/1081870.1081948.

Inderjit S. Dhillon, Yuqiang Guan, and Brian Kulis. Weighted graph cuts without eigenvectors a
multilevel approach. IEEE Transactions on Pattern Analysis and Machine Intelligence, 29(11):
1944–1957, 2007. doi: 10.1109/TPAMI.2007.1115.

Adji B. Dieng, Francisco J. R. Ruiz, David M. Blei, and Michalis K. Titsias. Prescribed generative
adversarial networks, 2019.

David K Duvenaud, Dougal Maclaurin, Jorge Iparraguirre, Rafael Bombarell, Timothy Hirzel, Alan
Aspuru-Guzik, and Ryan P Adams. Convolutional networks on graphs for learning molecular fingerprints. In Advances in Neural Information Processing Systems, volume 28. Curran Associates,
Inc., 2015.

Vijay Prakash Dwivedi, Chaitanya K. Joshi, Thomas Laurent, Yoshua Bengio, and Xavier Bresson.
Benchmarking graph neural networks. CoRR, abs/2003.00982, 2020.

Jack Edmonds and Richard M. Karp. Theoretical improvements in algorithmic efficiency for network flow problems. Journal of the ACM (JACM), 1972. doi: 10.1145/321694.321699.

Alex Fout, Jonathon Byrd, Basir Shariat, and Asa Ben-Hur. Protein interface prediction using graph
convolutional networks. In Proceedings of the 31st International Conference on Neural Informa_tion Processing Systems, NIPS’17, pp. 65336542, Red Hook, NY, USA, 2017. Curran Associates_
Inc. ISBN 9781510860964.


-----

J. Gilmer, S. S. Schoenholz, P. F. Riley, O. Vinyals, and G. E. Dahl. Neural message passing for
quantum chemistry. 70, 2017a.

Justin Gilmer, Samuel S. Schoenholz, Patrick F. Riley, Oriol Vinyals, and George E. Dahl. Neural
message passing for quantum chemistry. In Proceedings of the 34th International Conference on
_Machine Learning, volume 70 of Proceedings of Machine Learning Research, pp. 1263–1272._
PMLR, 06–11 Aug 2017b.

Xavier Glorot and Yoshua Bengio. Understanding the difficulty of training deep feedforward neural
networks. In Proceedings of the Thirteenth International Conference on Artificial Intelligence and
_Statistics, volume 9 of Proceedings of Machine Learning Research, pp. 249–256, Chia Laguna_
Resort, Sardinia, Italy, 13–15 May 2010. PMLR.

Ian Goodfellow, Jean Pouget-Abadie, Mehdi Mirza, Bing Xu, David Warde-Farley, Sherjil Ozair,
Aaron Courville, and Yoshua Bengio. Generative adversarial nets. In Advances in Neural Infor_mation Processing Systems, volume 27. Curran Associates, Inc., 2014._

E. J. Gumbel. Statistical theory of extreme values and some practical applications: a series of
lectures. US Govt. Print. Office, Number 33, 1954.

Rafael Gmez-Bombarelli, Jennifer N. Wei, David Duvenaud, Jos Miguel Hernndez-Lobato, Benjamn Snchez-Lengeling, Dennis Sheberla, Jorge Aguilera-Iparraguirre, Timothy D. Hirzel,
Ryan P. Adams, and Aln Aspuru-Guzik. Automatic chemical design using a data-driven continuous representation of molecules. ACS Central Science, 4(2):268–276, 2018. doi: 10.1021/
acscentsci.7b00572. PMID: 29532027.

David K. Hammond, Pierre Vandergheynst, and Rmi Gribonval. Wavelets on graphs via spectral graph theory. _Applied and Computational Harmonic Analysis, 30(2):129–150, 2011._
ISSN 1063-5203. doi: https://doi.org/10.1016/j.acha.2010.04.005. [URL https://www.](https://www.sciencedirect.com/science/article/pii/S1063520310000552)
[sciencedirect.com/science/article/pii/S1063520310000552.](https://www.sciencedirect.com/science/article/pii/S1063520310000552)

Martin Heusel, Hubert Ramsauer, Thomas Unterthiner, Bernhard Nessler, and Sepp Hochreiter.
Gans trained by a two time-scale update rule converge to a local nash equilibrium. In Advances
_in Neural Information Processing Systems, volume 30. Curran Associates, Inc., 2017._

P. Hohenberg and W. Kohn. Inhomogeneous electron gas. Phys. Rev., 136:864–871, 1964.

Truong Son Hy, Shubhendu Trivedi, Horace Pan, Brandon M. Anderson,, and Risi Kondor. Predicting molecular properties with covariant compositional networks. The Journal of Chemical
_Physics, 148, 2018._

John Ingraham and Debora Marks. Variational inference for sparse and undirected models. In Pro_ceedings of the 34th International Conference on Machine Learning, volume 70 of Proceedings_
_of Machine Learning Research, pp. 1607–1616. PMLR, 06–11 Aug 2017._

Eric Jang, Shixiang Gu, and Ben Poole. Categorical reparameterization with gumbel-softmax. In
_ICLR, 2017._

Wengong Jin, Regina Barzilay, and Tommi Jaakkola. Junction tree variational autoencoder for
molecular graph generation. In Proceedings of the 35th International Conference on Machine
_Learning, volume 80 of Proceedings of Machine Learning Research, pp. 2323–2332. PMLR,_
10–15 Jul 2018.

Steven Kearnes, Kevin McCloskey, Marc Berndl, Vijay Pande, and Patrick Riley. Molecular graph
convolutions: moving beyond fingerprints. Journal of Computer-Aided Molecular Design, 30(8):
595608, Aug 2016. ISSN 1573-4951. doi: 10.1007/s10822-016-9938-8.

Diederik P. Kingma and Jimmy Ba. Adam: A method for stochastic optimization. In 3rd Inter_national Conference on Learning Representations, ICLR 2015, San Diego, CA, USA, May 7-9,_
_2015, Conference Track Proceedings, 2015._

Diederik P. Kingma and Max Welling. Auto-Encoding Variational Bayes. In 2nd International
_Conference on Learning Representations, ICLR 2014, Banff, AB, Canada, April 14-16, 2014,_
_Conference Track Proceedings, 2014._


-----

Thomas N. Kipf and Max Welling. Variational graph auto-encoders, 2016.

Alexej Klushyn, Nutan Chen, Richard Kurle, Botond Cseke, and Patrick van der Smagt. Learning
hierarchical priors in vaes. In Advances in Neural Information Processing Systems, volume 32.
Curran Associates, Inc., 2019.

Daphne Koller and Nir Friedman. In Probabilistic Graphical Models: Principles and Techniques.
MIT Press, 2009.

Risi Kondor, Truong Son Hy, Horace Pan, Brandon M. Anderson, and Shubhendu Trivedi. Covariant
compositional networks for learning graphs, 2018.

Nils M. Kriege, Pierre-Louis Giscard, and Richard C. Wilson. On valid optimal assignment kernels
and applications to graph classification. In Proceedings of the 30th International Conference on
_Neural Information Processing Systems, NIPS’16, pp. 16231631, Red Hook, NY, USA, 2016._
Curran Associates Inc. ISBN 9781510838819.

Yann LeCun, Corinna Cortes, and Christopher J.C. Burges. The mnist database of handwritten
[digits. URL http://yann.lecun.com/exdb/mnist/.](http://yann.lecun.com/exdb/mnist/)

Yujia Li, Richard Zemel, Marc Brockschmidt, and Daniel Tarlow. Gated graph sequence neural
networks. In Proceedings of ICLR’16, April 2016.

Yujia Li, Oriol Vinyals, Chris Dyer, Razvan Pascanu, and Peter W. Battaglia. Learning deep generative models of graphs. ICML’18, abs/1803.03324, 2018.

Renjie Liao, Yujia Li, Yang Song, Shenlong Wang, Will Hamilton, David K Duvenaud, Raquel
Urtasun, and Richard Zemel. Efficient graph generation with graph recurrent attention networks.
In Advances in Neural Information Processing Systems, volume 32. Curran Associates, Inc., 2019.

Zinan Lin, Ashish Khetan, Giulia Fanti, and Sewoong Oh. Pacgan: The power of two samples in
generative adversarial networks. In Advances in Neural Information Processing Systems, volume 31. Curran Associates, Inc., 2018.

Jenny Liu, Aviral Kumar, Jimmy Ba, Jamie Kiros, and Kevin Swersky. Graph normalizing flows. In
_Advances in Neural Information Processing Systems, volume 32. Curran Associates, Inc., 2019._

Qi Liu, Miltiadis Allamanis, Marc Brockschmidt, and Alexander Gaunt. Constrained graph variational autoencoders for molecule design. In Advances in Neural Information Processing Systems,
volume 31. Curran Associates, Inc., 2018.

Jianxin Ma, Peng Cui, Kun Kuang, Xin Wang, and Wenwu Zhu. Disentangled graph convolutional
networks. In Proceedings of the 36th International Conference on Machine Learning, volume 97
of Proceedings of Machine Learning Research, pp. 4212–4221. PMLR, 09–15 Jun 2019.

Chris J Maddison, Daniel Tarlow, and Tom Minka. A* sampling. In Advances in Neural Information
_Processing Systems, volume 27. Curran Associates, Inc., 2014._

Haggai Maron, Heli Ben-Hamu, Hadar Serviansky, and Yaron Lipman. Provably powerful graph
networks. In Advances in Neural Information Processing Systems, volume 32. Curran Associates,
Inc., 2019a.

Haggai Maron, Heli Ben-Hamu, Nadav Shamir, and Yaron Lipman. Invariant and equivariant graph
networks. In International Conference on Learning Representations, 2019b.

Haggai Maron, Ethan Fetaya, Nimrod Segol, and Yaron Lipman. On the universality of invariant
networks. In Proceedings of the 36th International Conference on Machine Learning, volume 97
of Proceedings of Machine Learning Research, pp. 4363–4371. PMLR, 09–15 Jun 2019c.

Haggai Maron, Or Litany, Gal Chechik, and Ethan Fetaya. On learning sets of symmetric elements.
In Proceedings of the 37th International Conference on Machine Learning, volume 119 of Pro_ceedings of Machine Learning Research, pp. 6734–6744. PMLR, 13–18 Jul 2020._


-----

F. Monti, D. Boscaini, J. Masci, E. Rodola, J. Svoboda, and M. M. Bronstein. Geometric deep learning on graphs and manifolds using mixture model cnns. In 2017 IEEE Conference on Computer
_Vision and Pattern Recognition (CVPR), pp. 5425–5434, Los Alamitos, CA, USA, jul 2017. IEEE_
Computer Society. doi: 10.1109/CVPR.2017.576.

Kevin P. Murphy. Chapter 19: Undirected graphical models (markov random fields). In Machine
_Learning: A Probabilistic Perspective, pp. 663–707. MIT Press, 2012._

M. Niepert, M. Ahmed, and K. Kutzkov. Learning convolutional neural networks for graphs. In
_Proceedings of the International Conference on Machine Learning, 2016a._

Mathias Niepert, Mohamed Ahmed, and Konstantin Kutzkov. Learning convolutional neural networks for graphs. In Proceedings of The 33rd International Conference on Machine Learning,
volume 48 of Proceedings of Machine Learning Research, pp. 2014–2023, New York, New York,
USA, 20–22 Jun 2016b. PMLR.

B. Perozzi, R. Al-Rfou, and S. Skiena. Deepwalk: Online learning of social representations. Pro_ceedings of the 20th ACM SIGKDD International Conference on Knowledge Discovery and Data_
_Mining (KDD), pp. 701–710, 2014._

Alec Radford, Luke Metz, and Soumith Chintala. Unsupervised representation learning with deep
convolutional generative adversarial networks, 2016.

Raghunathan Ramakrishnan, Pavlo O Dral, Matthias Rupp, and O Anatole von Lilienfeld. Quantum
chemistry structures and properties of 134 kilo molecules. Scientific Data, 1, 2014.

Rajesh Ranganath, Dustin Tran, and David Blei. Hierarchical variational models. In Proceedings of
_The 33rd International Conference on Machine Learning, volume 48 of Proceedings of Machine_
_Learning Research, pp. 324–333, New York, New York, USA, 20–22 Jun 2016. PMLR._

Lars Ruddigkeit, Ruud van Deursen, Lorenz C. Blum, and Jean-Louis Reymond. Enumeration
of 166 billion organic small molecules in the chemical universe database gdb-17. Journal of
_Chemical Information and Modeling, 52(11):2864–2875, 2012. doi: 10.1021/ci300415d. PMID:_
23088335.

Havard Rue and Leonhard Held. Gaussian markov random fields: Theory and applications. In
_Monographs on Statistics and Applied Probability, volume 104, London, 2005. Chapman & Hall._

Raif Rustamov and Leonidas J Guibas. Wavelets on graphs via deep learning. In Advances in Neural
_Information Processing Systems, volume 26. Curran Associates, Inc., 2013._

Franco Scarselli, Marco Gori, Ah Chung Tsoi, Markus Hagenbuchner, and Gabriele Monfardini.
The graph neural network model. IEEE Transactions on Neural Networks, 20(1):61–80, 2009.
doi: 10.1109/TNN.2008.2005605.

[Maximilian Seitzer. pytorch-fid: FID Score for PyTorch. https://github.com/mseitzer/](https://github.com/mseitzer/pytorch-fid)
[pytorch-fid, August 2020. Version 0.1.1.](https://github.com/mseitzer/pytorch-fid)

P. Sen, G. M. Namata, M. Bilgic, L. Getoor, B. Gallagher,, and T. Eliassi-Rad. Collective classification in network data. AI Magazine, 29(3):93–106, 2008.

Hadar Serviansky, Nimrod Segol, Jonathan Shlomi, Kyle Cranmer, Eilam Gross, Haggai Maron,
and Yaron Lipman. Set2graph: Learning graphs from sets. In Advances in Neural Information
_Processing Systems, volume 33, pp. 22080–22091. Curran Associates, Inc., 2020._

Nino Shervashidze, Pascal Schweitzer, Erik Jan van Leeuwen, Kurt Mehlhorn, and Karsten M.
Borgwardt. Weisfeiler-lehman graph kernels. Journal of Machine Learning Research, 12(77):
2539–2561, 2011.

Chence Shi, Minkai Xu, Zhaocheng Zhu, Weinan Zhang, Ming Zhang, and Jian Tang. Graphaf:
a flow-based autoregressive model for molecular graph generation. In International Confer_[ence on Learning Representations, 2020. URL https://openreview.net/forum?id=](https://openreview.net/forum?id=S1esMkHYPr)_
[S1esMkHYPr.](https://openreview.net/forum?id=S1esMkHYPr)


-----

Donghyuk Shin, Si Si, and Inderjit S. Dhillon. Multi-scale link prediction. In Proceedings of
_the 21st ACM International Conference on Information and Knowledge Management, CIKM_
’12, pp. 215224, New York, NY, USA, 2012. Association for Computing Machinery. ISBN
9781450311564. doi: 10.1145/2396761.2396792.

Si Si, Donghyuk Shin, Inderjit S Dhillon, and Beresford N Parlett. Multi-scale spectral decomposition of massive graphs. In Advances in Neural Information Processing Systems, volume 27.
Curran Associates, Inc., 2014.

Martin Simonovsky and Nikos Komodakis. Graphvae: Towards generation of small graphs using
variational autoencoders. CoRR, abs/1802.03480, 2018.

Akash Srivastava, Lazar Valkov, Chris Russell, Michael U. Gutmann, and Charles Sutton. Veegan: Reducing mode collapse in gans using implicit variational learning. In Advances in Neural
_Information Processing Systems, volume 30. Curran Associates, Inc., 2017._

Teague Sterling and John J. Irwin. Zinc 15 ligand discovery for everyone. Journal of Chemical
_Information and Modeling, 55(11):2324–2337, 2015. doi: 10.1021/acs.jcim.5b00559. PMID:_
26479676.

L. Tang and H. Liu. Leveraging social media networks for classification. Data Mining and Knowl_edge Discovery, 23(3):447–478, 2011._

Erik Henning Thiede, Truong Son Hy, and Risi Kondor. The general theory of permutation equivarant neural networks and higher order graph variational encoders. CoRR, abs/2004.03990, 2020.

Arash Vahdat and Jan Kautz. Nvae: A deep hierarchical variational autoencoder. In Advances in
_Neural Information Processing Systems, volume 33, pp. 19667–19679. Curran Associates, Inc.,_
2020.

Petar Velikovi, Guillem Cucurull, Arantxa Casanova, Adriana Romero, Pietro Li, and Yoshua Bengio. Graph attention networks. In International Conference on Learning Representations, 2018.

Martin J. Wainwright and Michael I. Jordan. A variational principle for graphical models. New
_Directions in Statistical Signal Processing, 2005._

Hao Wu, Jonas K¨ohler, and Frank Noe. Stochastic normalizing flows. In Advances in Neural
_Information Processing Systems, volume 33, pp. 5933–5944. Curran Associates, Inc., 2020._

Bingbing Xu, Huawei Shen, Qi Cao, Yunqi Qiu, and Xueqi Cheng. Graph wavelet neural network.
In International Conference on Learning Representations, 2019.

Yiding Yang, Zunlei Feng, Mingli Song, and Xinchao Wang. Factorizable graph convolutional
networks. In Advances in Neural Information Processing Systems, volume 33, pp. 20286–20296.
Curran Associates, Inc., 2020.

Zhitao Ying, Jiaxuan You, Christopher Morris, Xiang Ren, Will Hamilton, and Jure Leskovec. Hierarchical graph representation learning with differentiable pooling. In Advances in Neural Infor_mation Processing Systems, volume 31. Curran Associates, Inc., 2018._

Jiaxuan You, Rex Ying, Xiang Ren, William Hamilton, and Jure Leskovec. GraphRNN: Generating realistic graphs with deep auto-regressive models. In Proceedings of the 35th International
_Conference on Machine Learning, volume 80 of Proceedings of Machine Learning Research, pp._
5708–5717. PMLR, 10–15 Jul 2018a.

Jiaxuan You, Rex Ying, Xiang Ren, William L. Hamilton, and Jure Leskovec. Code for graphrnn:
Generating realistic graphs with deep auto-regressive model. 2018b.

Manzil Zaheer, Satwik Kottur, Siamak Ravanbakhsh, Barnabas Poczos, Russ R Salakhutdinov, and
Alexander J Smola. Deep sets. In Advances in Neural Information Processing Systems, volume 30. Curran Associates, Inc., 2017.

Dawei Zhou, Lecheng Zheng, Jiejun Xu, and Jingrui He. Misc-gan: A multi-scale generative model
for graphs. Frontiers in Big Data, 2:3, 2019. ISSN 2624-909X. doi: 10.3389/fdata.2019.00003.
[URL https://www.frontiersin.org/article/10.3389/fdata.2019.00003.](https://www.frontiersin.org/article/10.3389/fdata.2019.00003)


-----

A BASIC TENSOR OPERATIONS

In order to build higher order equivariant networks, we revisit some basic tensor operations: tensor
product (see Def. 10) and tensor contraction (see Def. 11). It can be shown that these tensor operations respect permutation equivariance. Based on them, we build our second order message passing
networks.
**Definition 10. The tensor product of A** R[n][a] _with B_ R[n][b] _yields a tensor C_ _A_ _B_ R[n][a][+][b]

_where_
_Ci1,i2,..,ia_ _b ∈_ _Ai1,i2,..,ia_ _B ∈ia_ 1,ia 2,..,ia _b_ = ⊗ ∈
+ + + +

**Definition 11. The contraction of A** R[n][a] _along the pair of dimensions_ _x,y_ _(assuming x_ _y)_

=

_yields a_ _a_ 2 _-th order tensor_
∈ { } <

_Ci1,..,ix_ 1,j,ix 1,..,iy 1,j,iy 1,..,ia _Ai1,..,ia_

( − )

− + − + _ix,iy_

_where we assume that ix and iy have been removed from amongst the indices of = ∑_ _C. Using Einstein_
_notation, this can be written more compactly as_

_C_ _i1,i2,..,ia_ _ix,iy_ _Ai1,i2,..,ia_ _δ[i][x][,i][y]_

_where δ is the Kronecker delta. In general, the contraction of{_ }∖{ } _A along dimensions_ _x1,..,xp_ _yields_

=

_a tensor C_ _A_ _x1,..,xp_ R[n][a][−][p] _where_
{ }
↓
= ∈ _A_ _x1,..,xp_ _..._ _Ai1,i2,..,ia_

_ix1_ _ix2_ _ixp_
↓

_or compactly as_ = ∑ ∑ ∑
_A_ _x1,..,xp_ _Ai1,i2,..,ia_ _δ[i][x][1]_ _[,i][x][2]_ _[,..,i][xp]_ _._
↓

B MARKOV RANDOM FIELDS =

Undirected graphical models have been widely applied in the domains spatial or relational data,
such as image analysis and spatial statistics. In general, k-th order graph encoders encode an undirected graph _,_ into a k-th order latent z R[n][k][×][d][z], with learnable parameters θ, can be
represented as a parameterized Markov Random Field (MRF) or Markov network. Based on the
Hammersley-Clifford theorem (Murphy, 2012) (Koller & Friedman, 2009), a positive distribution G = (V E) ∈
_p_ **_z_** 0 satisfies the conditional independent properties of an undirected graph iff p can be
represented as a product of potential functions ψ, one per maximal clique, i.e.,
( ) > G

1
_p_ **_z_** **_θ_** _ψc_ _zc_ _θc_ (10)

_Z_ **_θ_** _c_

where is the set of all the (maximal) cliques of( ∣ ) = ( ), and[∏]∈C _Z(_ **_θ∣_** is the) _partition function to ensure the_
overall distribution sums to 1, and given by
C G ( )

_Z_ **_θ_** _ψc_ _zc_ _θc_
**_z_** _c_

Eq. 10 can be further written down as ( ) = ∑ ∏∈C ( ∣ )

_p_ **_z_** **_θ_** _ψv_ _zv_ **_θ_** _ψst_ _zst_ **_θ_** _ψc_ _zc_ **_θ_**
_v_ _s,t_ _c_ _i1,..,ik_ _k_

where ψv, ψst, and( _ψ∣_ _c) ∝∏ are the first order, second order and∈V_ ( ∣ ) ∏( )∈E ( ∣ ) ⋅⋅⋅ k=(-th order outputs of the encoder, cor-∏ )∈C ( ∣ )
responding to every vertex in, every edge in and every clique of size k in _k, respectively._
However, factorizing a graph into set of maximal cliques has an exponential time complexity, since
the problem of determining if there is a clique of size V E _k in a graph is known as an NP-complete C_
problem. Thus, the factorization based on Hammersley-Clifford theorem is intractable. The second
order encoder relaxes the restriction of maximal clique into edges, that is called as pairwise MRF:

_p_ **_z_** **_θ_** _ψst_ _zs,zt_
_s_ _t_
( ∣ ) ∝∏∼ ( )


-----

Our second order encoder inherits Gaussian MRF introduced by (Rue & Held, 2005) as pairwise
MRF of the following form

_p_ **_z_** **_θ_** _ψst_ _zs,zt_ _ψt_ _zt_
_s_ _t_ _t_

where ψst _zs,zt_ exp 2 _[z][s][Λ]([st]∣[z][t]) ∝∏[)][ is the edge potential, and]∼_ ( )∏ ( _[ ψ])_ _[t][(][z][t][) =][ exp]_ [( −] [1]2 [Λ][tt][z]t[2]

is the vertex potental. The joint distribution can be written in the information form of a multivariate
Gaussian in which( ) = ( − [1] [+][ η][t][z][t][)]
**Λ** **Σ[−][1]**

**_η_** **Λµ**

=

_p_ **_z_** **_θ_** exp =η[T] **_z_** (11)

2 **_[z][T][ Λ][z][)]_**

Sampling z from p **_z_** **_θ_** in Eq. 11 is the same as sampling from the multivariate Gaussian( ∣ ) ∝ ( − [1] **_µ,_** **Σ** .
To ensure end-to-end equivariance, we set the latent layer to be two tensors µ R[n][×][d][z] and Σ
R[n][×][n][×][d][z] that corresponds to( ∣ ) _dz multivariate Gaussians, whose first index, and second index are N(_ )
first order and second order equivariant with permutations. Computation of Σ is trickier than ∈ **_µ ∈,_**
simply because Σ must be invertible to be a covariance matrix. Thus, our second order encoder
produces tensor L as the second order activation, and set Σ **_LL[T]_** . The reparameterization trick
from Kingma & Welling (2014) is changed to
=

**_z_** **_µ_** **_Lϵ,_** **_ϵ_** 0, 1
= + ∼N( )

C EQUIVARIANT LEARNABLE PRIOR

The original VAE published by Kingma & Welling (2014) limits each covariance matrix Σ to be
diagonal and the prior to be 0, 1 . Our second order encoder removes the diagonal restriction on
the covariance matrix. Furthermore, we allow the prior **_µˆ_** _,_ **Σ[ˆ]** to be learnable in which ˆµ and
**Σˆ** are parameters optimized by back propagation in a data driven manner. Importantly, N( ) **Σ[ˆ]** cannot
be learned directly due to the invertibility restriction. Instead, similarly to the second order encoder, N( )
a matrix **_L[ˆ] is optimized, and the prior covariance matrix is constructed by setting_** **Σ[ˆ]** **_L[ˆ]L[ˆ]T . The_**
Kullback-Leibler divergence between the two distributions **_µ,_** **Σ** and **_µˆ,_** **Σ[ˆ]** is as follows:
=

1 1 **Σ**

_KL_ **_µ,_** **Σ** **_µˆ,_** **Σ[ˆ]** [1] **Σ** **Σ** **_µˆ_** **_µ_** N(Σ **_µˆ_** ) **_µ_** N(n ln ) (12)

2 [(][tr][(] **[ˆ]** − − det **Σ** [))]

Even thoughD (N( Σ is invertible, but gradient computation through the KL-divergence loss can be nu-)∣∣N( )) = ) + ( − )[T][ ˆ] ( − ) − + ( [det][ ˆ]
merical instable because of Cholesky decomposition procedure in matrix inversion. Thus, we add
neglectable noise ϵ 10[−][4] to the diagonal of both covariance matrices.

Importantly, during training, the KL-divergence loss breaks the permutation equivariance. Suppose

=

the set of vertices are permuted by a permutation matrix P σ for σ Sn. Since µ and Σ are the first
order and second order equivariant outputs of the encoder, they are changed to P σµ and P σΣP _[T]σ_
accordingly. But ∈

KL **_µ,_** **Σ** **_µˆ,_** **Σ[ˆ]** KL **_P σµ,_** **_P σΣP_** _[T]σ_ **_µ,_** **Σ[ˆ]**

To address the equivariance issue, we want to solve the following convex optimization problem thatD (N( )∣∣N( )) ≠D (N( [)∣∣N(][ˆ] ))
is our new equivariant loss function

min KL **_P σµ,_** **_P σΣP_** _[T]σ_ **_µ,_** **Σ[ˆ]** (13)
_σ_ Sn

However, solving the optimization based on Eq. 13 is computationally expensive. One solution∈ D (N( [)∣∣N(][ˆ] ))
is to solve the minimum-cost maximum-matching in a bipartite graph (Hungarian matching) with
the cost matrix C _ij_ _µi_ _µˆj_ by O _n[4]_ algorithm published by Edmonds & Karp (1972), that
can be still improved further into O _n[3]_ . The Hungarian matching preserves equiariance, but is
still computationally expensive. In practice, instead of finding a optimal permutation, we apply = ∣∣ − ∣∣ ( )
a free-matching scheme to find an assignment matrix( ) **Π such that: Πij** 1 if and only if j[∗]
arg minj _µi_ _µˆj_, for each i 1,n . The free-matching scheme preserves equivariance and can∗
be done efficiently in a simple O _n[2]_ algorithm that is also suitable for GPU computation.= =
∣∣ − ∣∣ ∈[ ]
( )


-----

D EXPERIMENTS

D.1 LINK PREDICTION ON CITATION GRAPHS

We demonstrate the ability of the MGVAE models to learn meaningful latent embeddings on a link
prediction task on popular citation network datasets Cora and Citeseer (Sen et al., 2008). At training
time, 15% of the citation links (edges) were removed while all node features are kept, the models are
trained on an incomplete graph Laplacian constructed from the remaining 85% of the edges. From
previously removed edges, we sample the same number of pairs of unconnected nodes (non-edges).
We form the validation and test sets that contain 5% and 10% of edges with an equal amount of
non-edges, respectively.

We compare our model MGVAE against popular methods in the field:

1. Spectral clustering (SC) (Tang & Liu, 2011)
2. Deep walks (DW) (Perozzi et al., 2014)
3. Variational graph autoencoder (VGAE) (Kipf & Welling, 2016)

on the ability to correctly classify edges and non-edges using two metrics: area under the ROC curve
(AUC) and average precision (AP). Numerical results of SC and DW are experimental settings are
taken from (Kipf & Welling, 2016). We reran the implementation of VGAE as in (Kipf & Welling,
2016).

For MGVAE, we initialize weights by Glorot initialization (Glorot & Bengio, 2010). We repeat the
experiments with 5 different random seeds and calculate the average AUC and AP along with their
standard deviations. The number of message passing layers ranges from 1 to 4. The size of latent
representation is 128. The number of coarsening levels is L 3, 7 . In the ℓ-th coarsening level, we
partition the graph into 2[ℓ] (for L 7) or 4[ℓ] (for L 3) clusters. We train for 2,048 epochs using
Adam optimization (Kingma & Ba, 2015) with a starting learning rate of ∈{ } 0.01. Hyperparameters
optimization (e.g. number of layers, dimension of the latent representation, etc.) is done on the G[(][ℓ][)] = =
validation set. MGVAE outperforms all other methods (see Table 3).

We propose our learning to cluster algorithm to achieve the balanced K-cut at every resolution level.
Besides, we also implement two fixed clustering algorithms:

1. Spectral: It is similar to the one implemented in (Rustamov & Guibas, 2013).

-  First, we embed each node i into R[n][max] as _ξ1_ _i_ _λ1_ _i_ _,..,ξnmax_ _i_ _λnmax_ _i_,
where _λn,ξn_ _n_ 0 are the eigen-pairs of the graph Laplacian where
_ii_ _j_ _ij. We assume that ∈V λ0_ _.._ _λnmax_ . In this case, ( ( )/ ( ) nmax 10(.)/ ( ))
=

-  At the { ℓ-th resolution level, we apply the K-Means clustering algorithm based on the}[n][max] L = D[−][1](D−A)
above node embedding to partition graphD = ∑ A ≤ ≤ . =

2. K-Means:

G[(][ℓ][)]

-  First, we apply PCA to compress the sparse word frequency vectors (of size 1,433 on
Cora and 3,703 on Citeseer) associating with each node into 10 dimensions.



-  We use the compressed node embedding for the K-Means clustering.

Tables 4 and 5 show that our learning to cluster algorithm returns a much more balanced cut on
the highest resolution level comparing to both Spectral and K-Means clusterings. For instance, we

Table 3: Citation graph link prediction results (AUC & AP)

± ± ± ±

|Dataset|Cora|Col3|Citeseer|Col5|
|---|---|---|---|---|
|Method|AUC (ROC)|AP|AUC (ROC)|AP|
|SC|84.6 0.01 ±|88.5 0.00 ±|80.5 0.01 ±|85.0 0.01 ±|
|DW|83.1 0.01 ±|85.0 0.00 ±|80.5 0.02 ±|83.6 0.01 ±|
|VGAE|90.97 0.77 ±|91.88 0.83 ±|89.63 1.04 ±|91.10 1.02 ±|
|MGVAE (Spectral)|91.19 0.76 ±|92.27 0.73 ±|90.55 1.17 ±|91.89 1.27 ±|
|MGVAE (K-Means)|93.07 5.61 ±|92.49 5.77 ±|90.81 1.19 ±|91.98 1.02 ±|
|MGVAE|95.67 3.11 ±|95.02 3.36 ±|93.93 5.87 ±|93.06 6.33 ±|


**Dataset** **Cora** **Citeseer**

**Method** AUC (ROC) AP AUC (ROC) AP

SC 84.6 0.01 88.5 0.00 80.5 0.01 85.0 0.01

DW 83.1 0.01 85.0 0.00 80.5 0.02 83.6 0.01

VGAE 90.97 ± 0.77 91.88 ± 0.83 89.63 ± 1.04 91.10 ± 1.02

**MGVAE (Spectral)** 91.19 ± 0.76 92.27 ± 0.73 90.55 ± 1.17 91.89 ± 1.27

**MGVAE (K-Means)** 93.07 ± 5.61 92.49 ± 5.77 90.81 ± 1.19 91.98 ± 1.02

**MGVAE** **95.67 ± 3.11** **95.02 ± 3.36** **93.93 ± 5.87** **93.06 ± 6.33**

± ± ± ±


-----

|Method|Min|Max|STD|KL divergence|
|---|---|---|---|---|
|Spectral|1|2020|177.52|3.14|
|K-Means|1|364|40.17|0.84|
|Learn to cluster|10|36|4.77|0.02|


**Method** Min Max STD KL divergence

Spectral 1 2020 177.52 3.14

K-Means 1 364 40.17 0.84

Learn to cluster 10 36 4.77 **0.02**

Table 4: Learning to cluster algorithm returns balanced cuts on Cora.

|Method|Min|Max|STD|KL divergence|
|---|---|---|---|---|
|Spectral|1|3320|292.21|4.51|
|K-Means|1|326|41.69|0.74|
|Learn to cluster|11|38|4.93|0.01|


**Method** Min Max STD KL divergence

Spectral 1 3320 292.21 4.51

K-Means 1 326 41.69 0.74

Learn to cluster 11 38 4.93 **0.01**


Table 5: Learning to cluster algorithm returns balanced cuts on Citeseer.

have L 7 resolution levels and we partition the ℓ-th resolution into K 2[ℓ] clusters. Thus, on the
bottom levels, we have 128 clusters. If we distribute nodes into clusters uniformly, the expected
number of nodes in a cluster is = 21.15 and 25.99 on Cora (2, 708 nodes) and Citeseer ( = 3, 327 nodes),
respectively. We measure the minimum, maximum, standard deviation of the numbers of nodes in
128 clusters. Furthermore, we measure the Kullback–Leibler divergence between the distribution of
nodes into clusters and the uniform distribution. Our learning to cluster algorithm achieves low KL
losses of 0.02 and 0.01 on Cora and Citeseer, respectively.

D.2 MOLECULAR GRAPH GENERATION

In this case, MGVAE and MGN are implemented with L 3 resolution levels, and the ℓ-th resolution
graph is partitioned into K 2[ℓ][−][1] clusters. On each resolution level, the local encoders and local
decoders are second-order Sn-equivariant networks with up to 4 equivariant layers. The number of =
channels for each node latent = dz is set to 256. We apply two approaches for graph decoding:

1. All-at-once: MGVAE reconstructs all resolution adjacencies by equivariant decoder networks. Furthermore, we apply learnable equivariant prior as in Sec. C. Our second order
encoders are interpreted as Markov Random Fields (see Sec. B). This approach preserves
permutation equivariance. In addition, we implement a correcting process: the decoder
network of the highest resolution level returns a probability for each edge, we sort these
probabilities in a descending order and gradually add the edges in that order to satisfy all
chemical constraints. Furthermore, we investigate the expressive power of the second order
Sn-equivariant decoder by replacing it by a multilayer perceptron (MLP) decoder with 2
hidden layers of size 512 and sigmoid nonlinearity. We find that the higher order decoder
outperforms the MLP decoder given the same encoding architecture. Table 6 shows the
comparison between the two decoding models.

2. Autoregressive: This decoding process is constructed in an autoregressive manner similarly to (Liu et al., 2018). First, we sample each vertex latent z independently. We randomly
select a starting node v0, then we apply Breath First Search (BFS) to determine a particular
node ordering from the node v0, however that breaks the permutation equivariance. Then
iteratively we add/sample new edge to the existing graph _t at the t-th iteration (given a_
randomly selected node v0 as the start graph 0) until completion. We apply second-order
MGN with gated recurrent architecture to produce the probability of edge G _u,v_ where one
vertex u is in the existing graph _t and the another one is outside; and also the probability G_
of its label. Intuitively, the decoding process is a sequential classification. ( )
G

We randomly select 10,000 training examples for QM9; and 1,000 (autoregressive) and 10,000 (allat-once) training examples for ZINC. It is important to note that our training sets are much smaller
comparing to other methods. For all of our generation experiments, we only use graph features as
the input for the encoder such as one-hot atomic types and bond types. Since ZINC molecules are
larger then QM9 ones, it is more difficult to train with the second order Sn-equivariant decoders
(e.g., the number of bond/non-bond predictions or the number of entries in the adjacency matrices
are proportional to squared number of nodes). Therefore, we input several chemical/atomic features


-----

|Dataset|Method|Validity|Novelty|Uniqueness|
|---|---|---|---|---|
|QM9|MLP decoder|100%|99.98%|77.62%|
||Sn decoder|100%|100%|95.16%|


**Dataset** **Method** **Validity** **Novelty** **Uniqueness**

MLP decoder 100% 99.98% 77.62%
QM9

Sn decoder 100% 100% 95.16%

Table 6: All-at-once MGVAE with MLP decoder vs. second order decoder.

|Feature|Type|Number|Description|
|---|---|---|---|
|GetAtomicNum|Integer|1|Atomic number|
|IsInRing|Boolean|1|Belongs to a ring?|
|IsInRingSize|Boolean|9|Belongs to a ring of size k 1, .., 9 ? ∈{ }|
|GetIsAromatic|Boolean|1|Aromaticity?|
|GetDegree|Integer|1|Vertex degree|
|GetExplicitValance|Integer|1|Explicit valance|
|GetFormalCharge|Integer|1|Formal charge|
|GetIsotope|Integer|1|Isotope|
|GetMass|Double|1|Atomic mass|
|GetNoImplicit|Boolean|1|Allowed to have implicit Hs?|
|GetNumExplicitHs|Integer|1|Number of explicit Hs|
|GetNumImplicitHs|Integer|1|Number of implicit Hs|
|GetNumRadicalElectrons|Integer|1|Number of radical electrons|
|GetTotalDegree|Integer|1|Total degree|
|GetTotalNumHs|Integer|1|Total number of Hs|
|GetTotalValence|Integer|1|Total valance|


**Feature** **Type** **Number** **Description**

GetAtomicNum Integer 1 Atomic number

IsInRing Boolean 1 Belongs to a ring?

IsInRingSize Boolean 9 Belongs to a ring of size k 1, .., 9 ?

GetIsAromatic Boolean 1 Aromaticity?

GetDegree Integer 1 Vertex degree ∈{ }

GetExplicitValance Integer 1 Explicit valance

GetFormalCharge Integer 1 Formal charge

GetIsotope Integer 1 Isotope

GetMass Double 1 Atomic mass

GetNoImplicit Boolean 1 Allowed to have implicit Hs?

GetNumExplicitHs Integer 1 Number of explicit Hs

GetNumImplicitHs Integer 1 Number of implicit Hs

GetNumRadicalElectrons Integer 1 Number of radical electrons

GetTotalDegree Integer 1 Total degree

GetTotalNumHs Integer 1 Total number of Hs

GetTotalValence Integer 1 Total valance


Table 7: The list of chemical/atomic features used for the all-at-once MGVAE on ZINC. We denote
each feature by its API in RDKit.

computed from RDKit for the all-at-once MGVAE on ZINC (see Table 7). We concatenate all these
features into a vector of size 24 for each atom.

We train our models with Adam optimization method (Kingma & Ba, 2015) with the initial learning
rate of 10[−][3]. Figs. 4 and 5 show some selected examples out of 5,000 generated molecules on QM9
by all-at-once MGVAE, while Fig. 6 shows the molecules generated by autoregressive MGN. Qualitatively, both the decoding approaches capture similar molecular substructures (bond structures).
Fig. 9 shows an example of interpolation on the latent space on ZINC with the all-at-once MGVAE.
Fig. 7 shows some generated molecules on ZINC by the all-at-once MGVAE. Fig. 8 and table 8
show some generated molecules by the autoregressive MGN on ZINC dataset with high Quantitative Estimate of Drug-Likeness (QED) computed by RDKit and their SMILES strings. On ZINC,
the average QED score of the generated molecules is 0.45 with standard deviation 0.21. On QM9,
the QED score is 0.44 0.07.

D.3 GENERAL GRAPH GENERATION BY ± MGVAE

Figs. 10 11 show some generated examples and training examples on the 2-community and ego
datasets, respectively.

D.4 UNSUPERVISED MOLECULAR PROPERTIES PREDICTION ON QM9

Density Function Theory (DFT) is the most successful and widely used approach of modern quantum chemistry to compute the electronic structure of matter, and to calculate many properties of
molecular systems with high accuracy (Hohenberg & Kohn, 1964). However, DFT is computationally expensive (Gilmer et al., 2017a), that leads to the use of machine learning to estimate the
properties of compounds from their chemical structure rather than computing them explicitly with
DFT (Hy et al., 2018). To demonstrate that MGVAE can learn a useful molecular representations
and capture important molecular structures in an unsupervised and variational autoencoding manner,
we extract the highest resolution latents (at ℓ _L) and use them as the molecular representations_
for the downstream tasks of predicting DFT’s molecular properties on QM9 including 13 learning
=


-----

Figure 4: Some generated examples on QM9 by the all-at-once MGVAE with second order Snequivariant decoders.

Figure 5: Some generated examples on QM9 by the all-at-once MGVAE with a MLP decoder
instead of the second order Sn-equivariant one. It generates more tree-like structures.


-----

Figure 6: Some generated examples on QM9 by the autoregressive MGN.

Figure 7: Some generated examples on ZINC by the all-at-once MGVAE with second order Snequivariant decoders. In addition of graph features such as one-hot atomic types, we include several
chemical features computed from RDKit (as in Table 7) as the input for the encoders. A generated
example can contain more than one connected components, each of them is a valid molecule.


-----

QED = 0.711 QED = 0.715 QED = 0.756 QED = 0.751

QED = 0.879 QED = 0.805 QED = 0.742 QED = 0.769

QED = 0.710 QED = 0.790 QED = 0.850 QED = 0.859

QED = 0.730 QED = 0.901 QED = 0.786 QED = 0.729

QED = 0.703 QED = 0.855 QED = 0.895 QED = 0.809

Figure 8: Some generated molecules on ZINC by the autoregressive MGN with high QED (druglikeness score).


-----

|Row|Column|SMILES|
|---|---|---|
|1|1|O=C1NC(CCCF)c2[nH]nnc21|
||2|OCC(OSBr)c1ccc(-c2cccc(Cl)c2Cl)[nH]1|
||3|C=CC1=CC=c2c(cc(=C3ONC(Cl)=C3Cl)[nH]c2=O)O1|
||4|COC(=CN1NC=CN1)C=C1C=CC(Cl)=CO1|
|2|1|[NH-]C(CNS1(=O)=NNc2c(F)cccc21)C1CC1|
||2|CS(=O)N1CC[SH](C)C(CNCc2ccccc2)C1|
||3|C=C(Cl)[SH](=O)(NC)C1c2ccc(Cl)c(n2)CC1O|
||4|CC(F)C(=C1C[NH2+]C([O-])N1)S(=O)Cc1ccccc1|
|3|1|CC1(NC2=CONN2c2ccccc2)C=C(C=O)N[N-]1|
||2|CC(=O)NN1N=C(C(O)c2cccc3ccoc23)C(=O)C1=O|
||3|C=CC(C)=C1C(F)=CC(C=C2ONN=NS2=O)=C1SCl|
||4|CCN1ON=C(C=C(Cl)c2ccco2)C(F)(F)C1=O|
|4|1|O=C(CCN(c1[nH+]cc(S)s1)c1ccc2cc1SC2)C1=NCC=C1Cl|
||2|CC=CNC1=C2Oc3ccccc3C(C)S2=S=N1|
||3|O=C(SC1=CC=NS1(=O)=O)c1ccc(Cl)cc1S1=NN=NN=N1|
||4|COCCNCc1cc2ccccc2[nH]1|
|5|1|ClC=C1CON=C(c2ncno2)N1CC(Cl)(Br)Br|
||2|CS(=O)(=O)c1ccc[nH+]c1SNCc1ccccc1Cl|
||3|O=S1(=O)CNS(=O)(N(c2ccccc2F)c2ccccc2Cl)=N1|
||4|O=C1NS(c2ccccc2Cl)=S2(=NSN=N2)O1|


**Row** **Column** **SMILES**

1 O=C1NC(CCCF)c2[nH]nnc21

2 OCC(OSBr)c1ccc(-c2cccc(Cl)c2Cl)[nH]1

1

3 C=CC1=CC=c2c(cc(=C3ONC(Cl)=C3Cl)[nH]c2=O)O1

4 COC(=CN1NC=CN1)C=C1C=CC(Cl)=CO1

1 [NH-]C(CNS1(=O)=NNc2c(F)cccc21)C1CC1

2 CS(=O)N1CC[SH](C)C(CNCc2ccccc2)C1

2

3 C=C(Cl)[SH](=O)(NC)C1c2ccc(Cl)c(n2)CC1O

4 CC(F)C(=C1C[NH2+]C([O-])N1)S(=O)Cc1ccccc1

1 CC1(NC2=CONN2c2ccccc2)C=C(C=O)N[N-]1

2 CC(=O)NN1N=C(C(O)c2cccc3ccoc23)C(=O)C1=O

3

3 C=CC(C)=C1C(F)=CC(C=C2ONN=NS2=O)=C1SCl

4 CCN1ON=C(C=C(Cl)c2ccco2)C(F)(F)C1=O

1 O=C(CCN(c1[nH+]cc(S)s1)c1ccc2cc1SC2)C1=NCC=C1Cl

2 CC=CNC1=C2Oc3ccccc3C(C)S2=S=N1

4

3 O=C(SC1=CC=NS1(=O)=O)c1ccc(Cl)cc1S1=NN=NN=N1

4 COCCNCc1cc2ccccc2[nH]1

1 ClC=C1CON=C(c2ncno2)N1CC(Cl)(Br)Br

2 CS(=O)(=O)c1ccc[nH+]c1SNCc1ccccc1Cl

5

3 O=S1(=O)CNS(=O)(N(c2ccccc2F)c2ccccc2Cl)=N1

4 O=C1NS(c2ccccc2Cl)=S2(=NSN=N2)O1

Table 8: [SMILES of the generated molecules included in Fig. 8. Online drawing tool: https:](https://pubchem.ncbi.nlm.nih.gov//edit3/index.html)
[//pubchem.ncbi.nlm.nih.gov//edit3/index.html](https://pubchem.ncbi.nlm.nih.gov//edit3/index.html)

Figure 9: Interpolation on the latent space: we randomly select two molecules from ZINC and we
reconstruct the corresponding molecular graphs on the interpolation line between the two latents.

Figure 10: The top row includes generated examples and the bottom row includes training examples
on the synthetic 2-community dataset.


-----

Generated examples

Training examples

Figure 11: EGO-SMALL.


-----

|Target|Unit|Mean|STD|Description|
|---|---|---|---|---|
|α|bohr3|75.2808|8.1729|Norm of the static polarizability|
|C v|cal/mol/K|31.6204|4.0674|Heat capacity at room temperature|
|G|eV|-70.8352|9.4975|Free energy of atomization|
|gap|eV|6.8583|1.2841|Difference between HOMO and LUMO|
|H|eV|-77.0167|10.4884|Enthalpy of atomization at room temperature|
|HOMO|eV|-6.5362|0.5977|Highest occupied molecular orbital|
|LUMO|eV|0.3220|1.2748|Lowest unoccupied molecular orbital|
|µ|D|2.6729|1.5034|Norm of the dipole moment|
|ω 1|cm−1|3504.1155|266.8982|Highest fundamental vibrational frequency|
|R2|bohr2|1189.4091|280.4725|Electronic spatial extent|
|U|eV|-76.5789|10.4143|Atomization energy at room temperature|
|U 0|eV|-76.1145|10.3229|Atomization energy at 0 K|
|ZPVE|eV|4.0568|0.9016|Zero point vibrational energy|


**Target** **Unit** **Mean** **STD** **Description**

_α_ bohr[3] 75.2808 8.1729 Norm of the static polarizability

Cv cal/mol/K 31.6204 4.0674 Heat capacity at room temperature

G eV -70.8352 9.4975 Free energy of atomization

gap eV 6.8583 1.2841 Difference between HOMO and LUMO

H eV -77.0167 10.4884 Enthalpy of atomization at room temperature

HOMO eV -6.5362 0.5977 Highest occupied molecular orbital

LUMO eV 0.3220 1.2748 Lowest unoccupied molecular orbital

_µ_ D 2.6729 1.5034 Norm of the dipole moment

_ω1_ cm[−][1] 3504.1155 266.8982 Highest fundamental vibrational frequency

R[2] bohr[2] 1189.4091 280.4725 Electronic spatial extent

U eV -76.5789 10.4143 Atomization energy at room temperature

U0 eV -76.1145 10.3229 Atomization energy at 0 K

ZPVE eV 4.0568 0.9016 Zero point vibrational energy

Table 9: Description and statistics of 13 learning targets on QM9.

|Col1|alpha|Cv|G|gap|H|HOMO|LUMO|mu|omega1|R2|U|U0|ZPVE|
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
|WL|3.75|2.39|4.84|0.92|5.45|0.38|0.89|1.03|192|154|5.41|5.36|0.51|
|NGF|3.51|1.91|4.36|0.86|4.92|0.34|0.82|0.94|168|137|4.89|4.85|0.45|
|PSCN|1.63|1.09|3.13|0.77|3.56|0.30|0.75|0.81|152|61|3.54|3.50|0.38|
|CCN 2D|1.30|0.93|2.75|0.69|3.14|0.23|0.67|0.72|120|53|3.02|2.99|0.35|
|MGVAE|2.83|0.91|1.78|0.66|1.87|0.34|0.58|0.95|195|90|1.89|1.90|0.14|


alpha Cv G gap H HOMO LUMO mu omega1 R2 U U0 ZPVE

WL 3.75 2.39 4.84 0.92 5.45 0.38 0.89 1.03 192 154 5.41 5.36 0.51

NGF 3.51 1.91 4.36 0.86 4.92 0.34 0.82 0.94 168 137 4.89 4.85 0.45

PSCN 1.63 1.09 3.13 0.77 3.56 0.30 0.75 0.81 152 61 3.54 3.50 0.38

CCN 2D **1.30** 0.93 2.75 0.69 3.14 **0.23** 0.67 **0.72** **120** **53** 3.02 2.99 0.35

**MGVAE** 2.83 **0.91 1.78 0.66 1.87** 0.34 **0.58** 0.95 195 90 **1.89 1.90** **0.14**


Table 10: Unsupervised molecular representation learning by MGVAE to predict molecular properties calculated by DFT on QM9 dataset.

targets. For the training, we normalize all learning targets to have mean 0 and standard deviation 1.
The name, physical unit, and statistics of these learning targets are detailed in Table 9.

The implementation of MGVAE is the same as detailed in Sec. D.2. MGVAE is trained to reconstruct
the highest resolution (input) adjacency, its coarsening adjacencies and the node atomic features. In
this case, we do not use any chemical features: the node atomic features are just one-hot atomic
types. After MGVAE is converged, to obtain the Sn-invariant molecular representation, we average
the node latents at the L-th level into a vector of size 256. Finally, we apply a simple Multilayer
Perceptron with 2 hidden layers of size 512, sigmoid nonlinearity and a linear layer on top to predict
the molecular properties based on the extracted molecular representation. We compare the results
in Mean Average Error (MAE) in the corresponding physical units with four methods on the same
split of training and testing from (Hy et al., 2018):

1. Support Vector Machine on optimal-assignment Weisfeiler-Lehman (WL) graph kernel
(Shervashidze et al., 2011) (Kriege et al., 2016)

2. Neural Graph Fingerprint (NGF) (Duvenaud et al., 2015)

3. PATCHY-SAN (PSCN) (Niepert et al., 2016a)

4. Second order Sn-equivariant Covariant Compositional Networks (CCN 2D) (Kondor et al.,
2018) (Hy et al., 2018)

Our unsupervised results show that MGVAE is able to learn a universal molecular representation in
an unsupervised manner and outperforms WL in 12, NGF in 10, PSCN in 8, and CCN 2D in 8 out
of 13 learning targets, respectively (see Table 10). There are other recent methods in the field that
use several chemical and geometric information but comparing to them would be unfair.

D.5 SUPERVISED MOLECULAR PROPERTIES PREDICTION ON ZINC

To further demonstrate the comprehensiveness of MGN, we apply our model in a supervised regression task to predict the solubility (LogP) on the ZINC dataset. We use the same split of 10K/1K/1K


-----

|Method|MLP|GCN|GAT|MoNet|DiscenGCN|FactorGCN|GatedGCN E|MGN|
|---|---|---|---|---|---|---|---|---|
|MAE|0.667|0.503|0.479|0.407|0.538|0.366|0.363|0.290|


**Method** MLP GCN GAT MoNet DiscenGCN FactorGCN GatedGCNE **MGN**

**MAE** 0.667 0.503 0.479 0.407 0.538 0.366 0.363 **0.290**

Table 11: Supervised MGN to predict solubility on ZINC dataset.

for training/validation/testing as in (Dwivedi et al., 2020). The implementation of MGN is almost
the same as detailed in Sec. D.2, except we include the latents of all resolution levels into the prediction. In particular, in each resolution level, we average all the node latents into a vector of size
256; then we concatentate all these vectors into a long vector of size 256 _L and apply a linear layer_
for the regression task. The baseline results are taken from (Yang et al., 2020) including:
×

1. Multilayer Perceptron (MLP),

2. Graph Convolution Networks (GCN),

3. Graph Attention Networks (GAT) (Velikovi et al., 2018),

4. MoNet (Monti et al., 2017),

5. Disentangled Graph Convolutional Networks (DisenGCN) (Ma et al., 2019),

6. Factorizable Graph Convolutional Networks (FactorGCN) (Yang et al., 2020),

7. GatedGCNE (Dwivedi et al., 2020) that uses additional edge information.

Our supervised result shows that MGN outperforms the state-of-the-art models in the field with a
margin of 20% (see Table 11).

D.6 GRAPH-BASED IMAGE GENERATION BY MGVAE

In this additional experiment, we apply MGVAE into the task of image generation. Instead of matrix
representation, an image I R[H][×][W] is represented by a grid graph of H _W nodes in which each_
node represents a pixel, each edge is between two neighboring pixels, and each node feature is the
corresponding pixel’s color (e.g., ∈ R[1] in gray scale, and R[3] in RGB scale). Fig. 12 demonstrates an ⋅
exmaple of graph representation for images. Since images have natural spatial clustering, instead of
learning to cluster, we implement a fixed clustering procedure as follows:

-  For the ℓ-th resolution level, we divide the grid graph of size H [(][ℓ][)] _W_ into clusters
of size h _w that results into a grid graph of size_ _[H]h[(][ℓ][)]_ _[W]w[ (][ℓ][)]_ [, supposingly][(][ℓ][)] _[ h][ and][ w][ are]_
×

divisible by H [(][ℓ][)] and W, respectively. Each resolution is associated with an image I [(][ℓ][)]

× ×

that is a zoomed out version of[(][ℓ][)] _I_ [(][ℓ][+][1][)].

-  The global encoder e[(][ℓ][)] is implemented with 10 layers of message passing that operates
on the whole H [(][ℓ][)] _W_ grid graph. We sum up all the node latents into a single latent
vectornetwork architecture of the generator of DCGAN model (Radford et al., 2016) to map Z [(][ℓ][)] R[d][z] . The global decoder× [(][ℓ][)] **_d[(][ℓ][)]_** is implemented by the convolutional neural Z [(][ℓ][)]

into an approximated image∈ _I[ˆ][(][ℓ][)]. The Sn-invariant pooler p[(][ℓ][)]_ is a network operating on
each small h _w grid graph to produce the corresponding node feature for the next level_
_ℓ_ 1. MGVAE is trained to reconstruct all resolution images. Fig. 13 shows an example of
reconstruction at each resolution on a test image of MNIST (after the network converged). ×
+

We evaluate our MGVAE architecture on the MNIST dataset (LeCun et al.) with 60,000 training
examples and 10,000 testing examples. The original image size is 28 28. We pad zero pixels to get
the image size of 2[5] 2[5] (e.g., H [(][5][)] _W_ 32). Each cluster is a small grid graph of size 2 2
(e.g., h _w_ 2). Accordingly, the image sizes for all resolutions are 32× 32, 16 16, 8 8, etc. In

[(][5][)]
this case, the whole network architecture is a 2-dimensional quadtree. The latent size× = = _dz is selected ×_
as 256. We train our model for 256 epochs by Adam optimizer (Kingma & Ba, 2015) with the initial = = × × ×
learning rate 10[−][3]. In the testing process, for the ℓ-th resolution, we sample a random vector of size
_dz from prior_ 0, 1 and use the decoder d[(][ℓ][)] to decode the corresponding image. We generate
10,000 examples for each resolution. We compute the Frechet Inception Distance (FID) proposed by
(Heusel et al., 2017) between the testing set and the generated set as the metric to evaluate the quality N( )


-----

|Method|FID (32 32) ↓ ×|FID (16 16) ↓ ×|FID (8 8) ↓ ×|
|---|---|---|---|
|DCGAN|113.129|N/A|N/A|
|VEEGAN|68.749|||
|PACGAN|58.535|||
|PresGAN|42.019|||
|MGVAE|39.474|64.289|39.038|


**Method** FID (32 32) FID (16 16) FID (8 8)

DCGAN 113.129↓ ↓ ↓

VEEGAN 68.749 × × ×
N/A N/A

PACGAN 58.535

PresGAN 42.019

**MGVAE** **39.474** 64.289 39.038

Table 12: Quantitative evaluation of the generated set by FID metric for each resolution level on
MNIST. It is important to note that the generation for each resolution is done separately: for the ℓ-th
resolution, we sample a random vector of size dz 256 from 0, 1, and use the global decoder
**_d[(][ℓ][)]_** to decode into the corresponding image size. The baselines are taken from (Dieng et al., 2019).
= N( )

Figure 12: An image of digit 8 from MNIST (left) and its grid graph representation at 16 16
resolution level (right).
×

of our generated examples. We use the FID implementation from (Seitzer, 2020). We compare our
MGVAE against variants of Generative Adversarial Networks (GANs) (Goodfellow et al., 2014)
including DCGAN (Radford et al., 2016), VEEGAN (Srivastava et al., 2017), PacGAN (Lin et al.,
2018), and PresGAN (Dieng et al., 2019). Table 12 shows our quantitative results in comparison
with other competing generative models. The baseline results are taken from Prescribed Generative
_Adversarial Networks paper (Dieng et al., 2019). MGVAE outperforms all the baselines for the_
highest resolution generation. Figs. 14 and 15 show some generated examples of the 32 32 and
16 16 resolution, respectively.
×
×


-----

32 32 16 16 8 8
× × ×

Target

Reconstruction

Figure 13: An example of reconstruction on each resolution level for a test image in MNIST.

Figure 14: Generated examples at the highest 32 32 resolution level.
×


-----

Figure 15: Generated examples at the 16 16 resolution level.
×


-----

