# Electron-phonon coupling in a magic-angle twisted-bilayer graphene device

Andreij C. Gadelha1,2,†† Viet-Hung Nguyen3,†† Eliel G. Neto4 Fabiano Santana1
Markus B. Raschke2 Michael Lamparski5 Vincent Meunier5 Jean-Christophe
Charlier3<EMAIL_ADDRESS>Ado Jorio1
<EMAIL_ADDRESS>1Physics Department, Universidade Federal de Minas
Gerais, Belo Horizonte, MG, 31270-901, Brazil.  2Department of Physics, and
JILA, University of Colorado at Boulder, Boulder, Colorado 80309, USA.
3Institute of Condensed Matter and Nanosciences, Université catholique de
Louvain (UCLouvain), Louvain-la-Neuve, Belgium.  4Physics Institute,
Universidade Federal da Bahia, Salvador - BA, 40170-115 Brazil.  5Department
of Physics, Applied Physics, and Astronomy, Jonsson Rowland Science Center,
Troy, NY, USA.  †† These two authors contributed equally to this work.

###### Abstract

The importance of phonons in the strong correlation phenomena observed in
twisted bilayer graphene (TBG) at the so-called magic-angle is under debate.
Here we apply gate-dependent micro-Raman spectroscopy to monitor the G band
linewidth in TBG devices of twist angles $\theta=0^{\circ}$, $\sim
1.1^{\circ}$ (magic-angle) and $\sim 7^{\circ}$ (large angle). The results
show a broad and p/n-asymmetric doping behavior at the magic-angle, in clear
contrast to the behavior observed in twist angles above and below. Atomistic
modeling reproduces the experimental observations, revealing how the unique
electronic structure of magic-angle TBGs influences the electron-phonon
coupling and, consequently, the G band linewidth. Our findings indicate a
connection between electron-phonon coupling and experimental observations of
strongly correlated phenomena in magic-angle TBG.

The unusual electronic structure of graphene, with massless Dirac fermions,
along with the relatively light C atoms bound together by strong covalent
$\sigma$ bonds is responsible for electron-phonon scattering that rules the
lifetime of phonons, in departure from the broadly adopted adiabatic
Born–Oppenheimer approximation Pisana et al. (2007). The breakdown fo the
adiabatic Born–Oppenheimer approximation generates a Kohn anomaly Lazzeri and
Mauri (2006), which has been studied in graphite Bonini et al. (2007), single-
layer graphene Pisana et al. (2007); Bonini et al. (2007); Lazzeri and Mauri
(2006); Hasdeo et al. (2016), double-layer graphene with the graphite-like AB
(Bernal) stacking Yan et al. (2008); Mafra et al. (2009); Das et al. (2009);
Mafra et al. (2012), and in metallic carbon nanotubes Lazzeri et al. (2006);
Tsang et al. (2007); Farhat et al. (2007) using Raman spectroscopy. Electrons
dress the phonon quasi-particles and, thus, modify phonon frequencies,
decrease their life-times and, ultimately, causes corresponding spectral
broadening Hasdeo et al. (2016). These effects in the Raman spectra were shown
to be dependent on the materials Fermi level ($E_{\mathrm{F}}$), because the
electron-phonon coupling can be changed by suppressing some Raman scattering
paths due to the Pauli exclusion principle, thereby, offering a unique way to
monitor accurately $E_{\mathrm{F}}$ in the family of $sp^{2}$ carbon
nanostructures Dresselhaus et al. (2010), including the influence of different
environments, such as different substrates Allard and Wirtz (2010); Forster et
al. (2013); Wang et al. (2017).

Twisted bilayer graphene (TBG) represents a new class of $sp^{2}$ carbon
nanostructures for the development of graphene-based twistronics, with novel
properties, such as (i) twist-angle $\theta$-dependent van Hove singularities
in the electronic density of states (DOS) for large-angle TBGs
($2^{\circ}<\theta<30^{\circ}$) Li et al. (2010); Jorio and Cançado (2013);
(ii) unusual superconductivity at the magic-angle ($\theta\sim 1.1^{\circ}$)
Cao et al. (2018); and (iii) lattice reconstruction with the formation of
strain soliton and topological points for low-angle TBGs
($\theta<1.1^{\circ}$) Yoo et al. (2019); Gargiulo and Yazyev (2017); Gadelha
et al. (2021a). Recent spectroscopic and imaging resolved with nanometer
spatial resolution (nano-Raman) measurements show that, in the low-angle
regime, the electronic and phononic structures feature significantly localized
properties Gadelha et al. (2021a). Consequently, the C–C stretching energy
uncertainty, as measured by the Raman G-band phonon linewidth
($\Gamma_{\rm{G}}$), varies spatially, following the moiré lattice structure
governed by atomic reconstruction. Additionally, it was shown, at the micro-
metric scale, that $\Gamma_{\rm{G}}$ reaches a maximum at the magic-angle,
decreasing towards the expected graphene linewidth far below and above
$1^{\circ}$, both in the large- and low-angle domains. These findings motivate
exploring the Kohn anomaly in TBGs using Raman spectroscopy gate-doping
experiments. These experiments enable the quantification of the importance of
electron-phonon in $\Gamma_{\rm{G}}$, as contrasted to other possible
structural effects Cocemasov et al. (2013); Lamparski et al. (2020); Gadelha
et al. (2021b), thus shedding light into the role of phonons in TBG strongly-
correlated phenomena Einenkel and Efetov (2011); Profeta et al. (2012); Wu et
al. (2018); Lian et al. (2019); Wu et al. (2019); Lewandowski et al. (2021a);
Choi and Choi (2021).

Figure 1 shows the schematics of these electrochemical TBG devices (a), and an
optical image of a representative device (b). Electrical contacts were pre-
patterned on a coverslip by optical lithography, followed by deposition of a 1
nm chromium adhesive layer and a 50 nm gold layer. TBG samples were prepared
by a dry “tear-and-stack” method, using a polydimethylsiloxane (PDMS) semi-
pyramidal stamp covered with a polycarbonate (PC) sheet, as described in
Gadelha et al. (2021c). The TBGs were then transferred onto the pre-patterned
electrical contacts for in situ electrical and optical studies. An
electrochemical gate composed of a phosphate-buffered saline (PBS) solution is
utilized (NaCl, 0.137 M; KCl, 0.0027 M; Na2HPO4, 0.01 M; KH2PO4, 0.0018 M). A
potential ($V_{\rm ion}$) is applied to a free electrode that is not connected
electrically to graphene, with $V_{\rm ion}$ kept low enough to avoid
electrical breakdown, while the electrical contacts connected to graphene are
grounded. The ionic liquid behavior is calibrated utilizing a Bernal bilayer
graphene deposited on 285 nm thick SiO2, which has a known capacitance
Taychatanapat and Jarillo-Herrero (2010).

Micro-Raman scattering experiments were performed on an inverted microscope,
with an oil-immersion objective ($N.A.=1.4$), excited with a HeNe laser
($\lambda=633$ nm). The back-scattered light is collected by an spectrometer
equipped with a charge-coupled device (CCD) and a 600 g/l grating (Andor
Solis). Figure 1(c) shows the spectral evolution of the magic-angle TBG Raman
G band spectra as a function of the carrier concentration ($\Delta n$).

Figure 1: (a) Schematic of the TBG electrochemical device with bottom
illumination and Raman detection. We exhibit the contacts (yellow stripes)
with the electrical wiring, and the ionic liquid drop (transparent blue)
covering the entire device. (b) Optical images of a magic-angle TBG device:
(main panel) G band Raman intensity; (inset) photo from the eyepiece of the
microscope. (c) magic-angle TBG G band spectra vs. $\Delta n$. Spectra
obtained with 3 accumulations of 20 s (for Bernal and large angle TBGs,
3$\times$10 s accumulation). We use laser powers lower than 1 mW in the sample
to avoid laser-induced sample heating.

Figs. 2(a-c) display the gate-dependent G-band full-width at half-maximum
($\Gamma_{\rm G}$), for (a) AB (Bernal) stacked bilayer graphene
($\theta=0^{\circ}$), (b) magic-angle TBG ($\theta\sim 1.1^{\circ}$) and (c)
large-angle TBG ($\theta\sim 7^{\circ}$). We fit the G bands with a single
Lorentzian to extract $\Gamma_{\rm G}$, and the band linewidths values
reported here are shifted by $-2$ cm-1 from the measured values to account of
the spectrometer resolution, measured with an Ar/Ne lamp. Although gate
experiments change the intensities, frequencies, and linewidths of different
Raman bands in graphene-related systems Hasdeo et al. (2016), here we focus on
the $\Gamma_{\rm{G}}$ behavior to collect information related to the electron-
phonon coupling.

Figure 2: (a-c) gate dependent G-band full-width at half-maximum, $\Gamma_{\rm
G}$, and (d-f) electronic bandstructures for Bernal (a,d –
$\theta=0^{\circ}$), magic-angle (b,e – 1.1∘) and large angle (c,f –
$7^{\circ}$) TBGs. In the a-c panels, dots are experimental data while red
lines and stars are theoretical predictions, up-shifted according to the
$\Delta$ value displayed in each panel.

The phonon frequency is given by
$\omega_{q}^{\nu}=\omega_{q}^{(0)\nu}+\omega_{q}^{(2)\nu}$, where
$\omega_{q}^{(0)\nu}$ and $\omega_{q}^{(2)\nu}$ are, respectively, the
unperturbed phonon energy and the frequency correction due to renormalization
by the electron-phonon interaction (Kohn anomaly) Park et al. (2008); Hasdeo
et al. (2016). The phonon linewidth is directly obtained from the imaginary
part of $\omega_{q}^{(2)\nu}$:

$\displaystyle\Gamma_{q}^{ph,\nu}=\Im(\omega_{q}^{(2)\nu})=4\pi\sum_{n,m}\int\frac{dk}{A_{BZ}}|M_{ep}^{\nu}(k,k+q)|^{2}\times$
$\displaystyle\left[f(E_{k}^{n})-f(E_{k+q}^{m})\right]\delta(\hbar\omega_{q}^{(0)\nu}-E^{m}_{k+q}+E^{n}_{k})\,$
(1)

where $E_{k}^{n}$ and $E_{k+q}^{m}$ are the energies for the electronic states
before and after the $q$-phonon mediated scattering, $f(E_{k}^{n})$ is the
Fermi distribution, $M_{ep}^{\nu}(k,k+q)$ is the electron-phonon matrix
element, and $A_{\rm BZ}$ is the area of the first Brillouin zone where the
$k$-integral is performed.

Generally, previous work on $sp^{2}$ carbon systems Bonini et al. (2007);
Pisana et al. (2007); Lazzeri and Mauri (2006); Yan et al. (2008); Das et al.
(2009); Mafra et al. (2012, 2009); Tsang et al. (2007); Farhat et al. (2007);
Dresselhaus et al. (2010); Lazzeri et al. (2006) considered
$M_{ep}^{\nu}(k,k+q)$ as weakly varying (almost constant over the experimental
range), and the phonon behavior (frequency and linewidth) can be reduced to
the calculation of the joint density of states (JDOS), while taking the Pauli
exclusion principle into account: $\Gamma_{\rm
G}^{ph}\propto\textrm{JDOS}(\hbar\omega_{G})\left[f(-\hbar\omega_{G}/2)-f(\hbar\omega_{G}/2)\right]$,
where $\omega_{G}$ is the frequency of G band phonons. This simplified
expression reproduces the behavior of single-layer, AB-bilayer, and large-
angle TBG, but it fails when applied to magic-angle TBG, as shown in this
work. In the former cases, the simple electronic structures, with linear or
parabolic energy dispersions, as illustrated in Fig. 2(d,f), are obtained in a
large range of energies around the $E_{\mathrm{F}}$ and, therefore, their low-
energy bands can be accurately modeled using a single-particle tight-binding
model, even when the system is doped. However, around and below the magic-
angle TBGs, strong electron localization occurs near the $E_{\mathrm{F}}$
Nguyen et al. (2021), therefore leading to strong Coulomb interactions
Kerelsky et al. (2019). These Coulomb interactions have been shown to smooth
out the spatial distribution of charge carriers and to lead to significant
changes in the electronic structure upon doping Cea et al. (2019); Rademaker
et al. (2019); Cea and Guinea (2020); Goodwin et al. (2020); Lewandowski et
al. (2021b); Cea and Guinea (2021). Due to these doping effects on the
electronic structure and the effects of atomic reconstruction on phonon
properties in small-angle TBG Lamparski et al. (2020), the doping dependence
of electron-phonon interactions in magic-angle TBG cannot be described simply
as in the cases of non-twisted systems and large angle TBG, as detailed below.

In order to take into account the doping effects, we first use an adjusted
electronic model, i.e., a single-particle tight-binding Hamiltonian modified
by adding a Hartree potential term Cea and Guinea (2020); Goodwin et al.
(2020), that depends on doping concentration SI . We then perform atomistic
calculations of the frequency $\omega_{q}^{(2)\nu}$ in Eq.(1), with the
electron-phonon matrix elements for the G band defined as Gunst et al. (2016);
Choi and Choi (2018)

$\displaystyle
M_{ep}^{G}(k,n,m)=\sqrt{\frac{\hbar}{2M_{C}\omega^{0}_{G}}}\sum_{i\alpha}{e^{G}_{i\alpha}\bra{m,k}\partial_{i\alpha}{\hat{H}}\ket{n,k}}$
(2)

where $M_{C}$ is the mass of a carbon atom, $\partial_{i\alpha}{\hat{H}}$ is
the derivative of the tight-binding Hamiltonian with respect to the position
of the i-th atom along the $\alpha$-axis, $\ket{n,k}$ are the eigenfunctions
of $\hat{H}$ with eigenvalues $E^{n}_{k}$, and the polarization vectors
$e^{G}_{\kappa\alpha}$ come from phonon calculations Lamparski et al. (2020).
The results obtained using Eq. 1 with the matrix elements from Eq. 2 are
included in Fig. 2 (red lines and stars in the (a-c) panels). We include few
stars in Fig. 2(b) because the model utilized here is limited to low doping.
The calculated results are offset from the experimental values by $\Delta$
(see values in the figure). For Bernal bilayer graphene ($\theta=0^{\circ}$),
$\Delta=+7$ cm-1, consistent with values reported in the literature and
attributed to inhomogeneous broadening and non-electron-phonon coupling
effects (anharmonicity) Das et al. (2009); Mafra et al. (2012); Hasdeo et al.
(2016). For the other TBGs, the values of $\Delta$ are higher, with
$\Delta=+10$ cm-1 for large angle TBG ($\theta\sim 7^{\circ}$) and
$\Delta=+14$ cm-1 for magic-angle TBG ($\theta\sim 1.1^{\circ}$).
Inhomogeneity, mistacking, reconstruction, strain and device quality are
factors that cannot be ruled out to explain the measured $\Delta$ for the
magic-angle TBG. The striking result is the considerably broader and
asymmetric $\Gamma_{\rm G}$-dependence on $\Delta n$ for magic-angle TBG, as
compared to AB bilayer and large angle TBG.

Figure 3: Doping effect on the electronic structure of magic-angle twisted
bilayer graphene: (left) n-doped, (center) undoped [zoomed in from Fig. 2(e)],
and (right) p-doped. The arrows indicate the G band energy (197 meV), enabling
the connection between the flat bands near the neutrality point (highlighted
by the transparent green stripe) and valence/conduction bands (transparent
orange stripes), depending on doping, as dictated by Pauli’s exclusion
principle.

Our results confirm that the electron-phonon coupling strength $\lambda$ is
larger for magic-angle TBG, when compared to the other structures in Fig. 2,
as predicted theoretically Choi and Choi (2018). The differences between the
TBGs electron-phonon coupling strengths are direct consequences of their
electronic bandstructures (see Figs. 2(d-f)). The broader dependence on
$\Delta n$ is due to the flat band near the neutrality point in magic-angle
TBGs, which requires a considerable amount of carriers to be filled or
emptied, associated with the presence of other electronic bands within the G
band phonon energy range (see Figs. 2(d-f) and Fig. 3). Although the results
here are limited to the Raman active in-plane optical phonon at $\Gamma$ (G
band), calculations indicate that the largest contribution to the total
coupling strength comes from the in-plane optical phonon modes near
$\omega_{\rm D}=167$ and $\omega_{\rm G}=197$ meV, with momenta corresponding
to $q=K,K^{\prime}$ and $\Gamma$, respectively. Together, they contribute to
about 70% of $\lambda$ ($\lambda_{\rm D}=0.914$ and $\lambda_{\rm G}=0.648$)
Choi and Choi (2021). These results also suggest that the electron-phonon
coupling is strong ($\lambda>1$) for TBG and twisted tri-layer graphene (TTG)
at their magic-angles, but it is an order of magnitude weaker for twisted
double bilayer graphene (TDBG) and twisted monolayer-bilayer graphene (TMBG),
since Bernal-stacked layers in TDBG and TMBG induce sublattice polarization in
at-band states, suppressing inter-sublattice electron-phonon matrix elements
Choi and Choi (2021).

The electron-phonon coupling depends more strongly on hole doping than on
electron doping Choi and Choi (2018), which is a direct manifestation of the
electron-hole asymmetry of low-energy bands in the relaxed magic-angle TBG.
However, our experiments find an opposite electron-hole asymmetry behavior,
i.e., the electron-phonon interaction is larger in the n-doped regime than in
the p-doped case, as also confirmed by our atomistic simulations. In contrast
to previous results Choi and Choi (2018), in this work we include the
smoothing of charge inhomogeneities due to strong Coulomb interactions in
magic-angle TBG. Besides the charge inhomogeneity smoothing in the moiré
super-cell, which leads to changes in the low-energy flat bands, significant
changes are also obtained on the high-energy bands when the system is doped.
First, as illustrated in Fig. 3, these changes in the electronic structure
lead to a reduction of available electronic states which mediate the G-phonon
scatterings, thus explaining the reduction of phonon linewidth when magic-
angle TBG is doped. Second, the phonon linewidth reduction is observed to be
stronger in the p-doped regime than in the n-doped regime (see Fig. 3), in
agreement with the experimentally observed electron-hole asymmetry [Fig.
2(b)].

We forsaken explain the sharp peak (kink) observed at zero doping for large-
angle TBG (see Fig. 2(c)). The $\Gamma_{\rm G}$ dependence on the
$E_{\mathrm{F}}$ in this case, as well as for Bernal stacking, can be
approximately described by the function
$f(-\hbar\omega_{G}/2)-f(\hbar\omega_{G}/2)$ that smoothly varies with the
$E_{\mathrm{F}}$ (including the region around $E_{\mathrm{F}}=0$). In the
Bernal bilayer system, $\Delta n$ scales as $\Delta n\propto E_{\mathrm{F}}$
because of its parabolic low energy bands. The large-angle TBGs, however,
exhibits a Dirac-like linear energy dispersion, similar to that observed in
monolayer graphene, and therefore $\Delta n\propto E_{\mathrm{F}}^{2}$ Kim et
al. (2012) and, as a consequence, a kink at zero doping is observed for large-
angle TBG. Note that the simple descriptions above are no longer valid in the
case of magic-angle TBG, because its flat electronic bands and doping effects
makes the model significantly more complicated (see the electronic band
structures illustrated in Fig. 2(d-f)).

In conclusion, gate-dependent Raman spectroscopy measurements and atomistic
simulations in TBGs show a broad and p/n-asymmetric $\Gamma_{\rm G}$ doping
behavior at the magic-angle, in clear contrast to the behavior observed for
both larger and smaller twist angles. These results support previous
theoretical studies that investigated the relevance of electron-phonon
coupling in magic-angle TBGs in order to accurately understand strongly-
correlated phenomena, such as superconductivity and its robustness in twisted
layered systems Choi and Choi (2018, 2021).

## References

  * Pisana et al. (2007) S. Pisana, M. Lazzeri, C. Casiraghi, K. S. Novoselov, A. K. Geim, A. C. Ferrari, and F. Mauri, Nature Materials 6, 198 (2007).
  * Lazzeri and Mauri (2006) M. Lazzeri and F. Mauri, Physical Review Letters 97, 266407 (2006).
  * Bonini et al. (2007) N. Bonini, M. Lazzeri, N. Marzari, and F. Mauri, Physical Review Letters 99, 176802 (2007).
  * Hasdeo et al. (2016) E. H. Hasdeo, A. R. Nugraha, M. S. Dresselhaus, and R. Saito, Physical Review B 94, 075104 (2016).
  * Yan et al. (2008) J. Yan, E. A. Henriksen, P. Kim, and A. Pinczuk, Physical Review Letters 101, 136804 (2008).
  * Mafra et al. (2009) D. Mafra, L. Malard, S. Doorn, H. Htoon, J. Nilsson, A. C. Neto, and M. Pimenta, Physical Review B 80, 241414 (2009).
  * Das et al. (2009) A. Das, B. Chakraborty, S. Piscanec, S. Pisana, A. Sood, and A. Ferrari, Physical Review B 79, 155417 (2009).
  * Mafra et al. (2012) D. Mafra, P. Gava, L. Malard, R. Borges, G. Silva, J. Leon, F. Plentz, F. Mauri, and M. Pimenta, Carbon 50, 3435 (2012).
  * Lazzeri et al. (2006) M. Lazzeri, S. Piscanec, F. Mauri, A. Ferrari, and J. Robertson, Physical Review B 73, 155426 (2006).
  * Tsang et al. (2007) J. Tsang, M. Freitag, V. Perebeinos, J. Liu, and P. Avouris, Nature Nanotechnology 2, 725 (2007).
  * Farhat et al. (2007) H. Farhat, H. Son, G. G. Samsonidze, S. Reich, M. Dresselhaus, and J. Kong, Physical Review Letters 99, 145506 (2007).
  * Dresselhaus et al. (2010) M. Dresselhaus, A. Jorio, and R. Saito, Annu. Rev. Condens. Matter Phys. 1, 89 (2010).
  * Allard and Wirtz (2010) A. Allard and L. Wirtz, Nano Letters 10, 4335 (2010).
  * Forster et al. (2013) F. Forster, A. Molina-Sanchez, S. Engels, A. Epping, K. Watanabe, T. Taniguchi, L. Wirtz, and C. Stampfer, Physical Review B 88, 085419 (2013).
  * Wang et al. (2017) X. Wang, J. W. Christopher, and A. K. Swan, Scientific reports 7, 1 (2017).
  * Li et al. (2010) G. Li, A. Luican, J. L. Dos Santos, A. C. Neto, A. Reina, J. Kong, and E. Andrei, Nature Physics 6, 109 (2010).
  * Jorio and Cançado (2013) A. Jorio and L. G. Cançado, Solid State Communications 175, 3 (2013).
  * Cao et al. (2018) Y. Cao, V. Fatemi, S. Fang, K. Watanabe, T. Taniguchi, E. Kaxiras, and P. Jarillo-Herrero, Nature 556, 43 (2018).
  * Yoo et al. (2019) H. Yoo, R. Engelke, S. Carr, S. Fang, K. Zhang, P. Cazeaux, S. H. Sung, R. Hovden, A. W. Tsen, T. Taniguchi, et al., Nature Materials 18, 448 (2019).
  * Gargiulo and Yazyev (2017) F. Gargiulo and O. V. Yazyev, 2D Materials 5, 015019 (2017).
  * Gadelha et al. (2021a) A. C. Gadelha, D. A. Ohlberg, C. Rabelo, E. G. Neto, T. L. Vasconcelos, J. L. Campos, J. S. Lemos, V. Ornelas, D. Miranda, R. Nadas, et al., Nature 590, 405 (2021a).
  * Cocemasov et al. (2013) A. I. Cocemasov, D. L. Nika, and A. A. Balandin, Physical Review B 88, 035428 (2013).
  * Lamparski et al. (2020) M. Lamparski, B. Van Troeye, and V. Meunier, 2D Materials 7, 025050 (2020).
  * Gadelha et al. (2021b) A. C. Gadelha, T. L. Vasconcelos, L. G. Can¸cado, and A. Jorio, The Journal of Physical Chemistry Letters 12, 7625 (2021b).
  * Einenkel and Efetov (2011) M. Einenkel and K. B. Efetov, Physical Review B 84, 214508 (2011).
  * Profeta et al. (2012) G. Profeta, M. Calandra, and F. Mauri, Nature Physics 8, 131 (2012).
  * Wu et al. (2018) F. Wu, A. MacDonald, and I. Martin, Physical Review Letters 121, 257001 (2018).
  * Lian et al. (2019) B. Lian, Z. Wang, and B. A. Bernevig, Physical Review Letters 122, 257002 (2019).
  * Wu et al. (2019) F. Wu, E. Hwang, and S. D. Sarma, Physical Review B 99, 165112 (2019).
  * Lewandowski et al. (2021a) C. Lewandowski, D. Chowdhury, and J. Ruhman, Physical Review B 103, 235401 (2021a).
  * Choi and Choi (2021) Y. W. Choi and H. J. Choi, Physical Review Letters 127, 167001 (2021).
  * Gadelha et al. (2021c) A. C. Gadelha, D. A. Ohlberg, F. C. Santana, G. S. Eliel, J. S. Lemos, V. Ornelas, D. Miranda, R. B. Nadas, K. Watanabe, T. Taniguchi, et al., ACS Applied Nano Materials 4, 1858 (2021c).
  * Taychatanapat and Jarillo-Herrero (2010) T. Taychatanapat and P. Jarillo-Herrero, Phys. Rev. Lett. 105, 166601 (2010).
  * Park et al. (2008) C.-H. Park, F. Giustino, M. L. Cohen, and S. G. Louie, Nano Letters 8, 4229 (2008).
  * Nguyen et al. (2021) V.-H. Nguyen, D. Paszko, M. Lamparski, B. Van Troeye, V. Meunier, and J.-C. Charlier, 2D Materials 8, 035046 (2021).
  * Kerelsky et al. (2019) A. Kerelsky, L. J. McGilly, D. M. Kennes, L. Xian, M. Yankowitz, S. Chen, K. Watanabe, T. Taniguchi, J. Hone, C. Dean, et al., Nature 572, 1476 (2019).
  * Cea et al. (2019) T. Cea, N. R. Walet, and F. Guinea, Phys. Rev. B 100, 205113 (2019).
  * Rademaker et al. (2019) L. Rademaker, D. A. Abanin, and P. Mellado, Phys. Rev. B 100, 205114 (2019).
  * Cea and Guinea (2020) T. Cea and F. Guinea, Phys. Rev. B 102, 045107 (2020).
  * Goodwin et al. (2020) Z. A. H. Goodwin, V. Vitale, X. Liang, A. A. Mostofi, and J. Lischner, Electronic Structure 2, 034001 (2020).
  * Lewandowski et al. (2021b) C. Lewandowski, S. Nadj-Perge, and D. Chowdhury, npj Quantum Materials 6, 82 (2021b).
  * Cea and Guinea (2021) T. Cea and F. Guinea, Proceedings of the National Academy of Sciences 118 (2021), ISSN 0027-8424, eprint https://www.pnas.org/content/118/32/e2107874118.full.pdf.
  * (43) See Supplemental Material at [URL will be inserted by publisher] for simulation details.
  * Gunst et al. (2016) T. Gunst, T. Markussen, K. Stokbro, and M. Brandbyge, Physical Review B 93, 035414 (2016).
  * Choi and Choi (2018) Y. W. Choi and H. J. Choi, Physical Review B 98, 241412 (2018).
  * Kim et al. (2012) S. Kim, I. Jo, D. C. Dillen, D. A. Ferrer, B. Fallahazad, Z. Yao, S. K. Banerjee, and E. Tutuc, Phys. Rev. Lett. 108, 116404 (2012).
  * Trambly de Laissardière et al. (2012) G. Trambly de Laissardière, D. Mayou, and L. Magaud, Phys. Rev. B 86, 125413 (2012).

## Acknowledgments

The authors thank Prof Leaonardo C. Campos for helpful discussions. This work
was supported by CNPq (302775/2018-8), CAPES (RELAII and 88881.198744/2018-01)
and FAPEMIG, Brazil. A.C.G. acknowledges partial support from DoE Award No.
DE-SC0008807. V-H.N. and J.-C.C. acknowledge financial support from the
Fédération Wallonie-Bruxelles through the ARC on 3D nano-architecturing of 2D
crystals (16/21-077), from the European Union’s Horizon 2020 Research Project
and Innovation Program — Graphene Flagship Core3 (881603), from the Flag-Era
JTC project “TATTOOS” (R.8010.19), and from the Belgium FNRS through the
research project (T.0051.18). V.M., M.L acknowledge support from NY State
Empire State Development’s Division of Science, Technology and Innovation
(NYSTAR). Computational resources have been provided by the CISM
supercomputing facilities of UCLouvain and the CECI consortium funded by
F.R.S.-FNRS of Belgium (2.5020.11).


# CommPOOL: An Interpretable Graph Pooling Framework for Hierarchical Graph
Representation Learning

Haoteng Tang,1 Guixiang Ma, 2 Lifang He 3 Heng Huang 1 Liang Zhan 1

###### Abstract

Recent years have witnessed the emergence and flourishing of hierarchical
graph pooling neural networks (HGPNNs) which are effective graph
representation learning approaches for graph level tasks such as graph
classification. However, current HGPNNs do not take full advantage of the
graph’s intrinsic structures (e.g., community structure). Moreover, the
pooling operations in existing HGPNNs are difficult to be interpreted. In this
paper, we propose a new interpretable graph pooling framework - CommPOOL, that
can capture and preserve the hierarchical community structure of graphs in the
graph representation learning process. Specifically, the proposed community
pooling mechanism in CommPOOL utilizes an unsupervised approach for capturing
the inherent community structure of graphs in an interpretable manner.
CommPOOL is a general and flexible framework for hierarchical graph
representation learning that can further facilitate various graph-level tasks.
Evaluations on five public benchmark datasets and one synthetic dataset
demonstrate the superior performance of CommPOOL in graph representation
learning for graph classification compared to the state-of-the-art baseline
methods, and its effectiveness in capturing and preserving the community
structure of graphs.

## Introduction

In recent years, Graph Neural Network (GNN) has emerged and been broadly used
as a generalized deep learning architecture for graph representation learning
in many fields, such as social network analysis (Chen, Ma, and Xiao 2018;
Huang et al. 2018), chemical molecule studies (Dai, Dai, and Song 2016;
Duvenaud et al. 2015; Gilmer et al. 2017) and brain network analysis (Ma et
al. 2019a; Liu et al. 2019). Generally, GNN models learn node embeddings by
passing, transforming and aggregating node features across the graph. The
generated node representations can then be forwarded to further layers for
specific learning tasks, such as node classification (Kipf and Welling 2016a;
Veličković et al. 2017) and link prediction (Kipf and Welling 2016b)..

Most of the existing GNN models (e.g., GCN (Kipf and Welling 2016a), GAT
(Veličković et al. 2017), GraphSage (Hamilton, Ying, and Leskovec 2017)) focus
on node-level representation learning and only propagate information across
edges of the graph in a flat way. When applying these GNNs for graph-level
tasks such as graph classifications, existing works usually apply simple
global pooling strategies (i.e., a summation over the learned node
representations) to obtain the graph-level embedding and use it for graph
label prediction (Li et al. 2015; Vinyals, Bengio, and Kudlur 2015; Zhang et
al. 2018) or graph similarity learning (Ma et al. 2019b). One main drawback in
these GNNs is that the hierarchical structure, often existing in graphs, is
ignored during the global pooling process, which makes the models less
effective for graph-level tasks. Hierarchical structure is a very important
structure for many graphs in various domains. For example, the hierarchical
community structure shown in Figure 1 is a typical pattern that often appears
in social networks (Girvan and Newman 2002; Long et al. 2019), chemical
molecule networks (Spirin and Mirny 2003) and brain networks (Kong and Yu
2014; Meunier et al. 2009). Therefore, preserving these community structures
is critical for better understanding and analyzing these graphs.

Figure 1: An example of hierarchical community structures in graphs

Some recent works proposed hierarchical graph pooling neural networks (HGPNNs)
to address the hierarchical structure representation issue by introducing the
hierarchical pooling operations (Ying et al. 2018; Lee et al. 2019; Gao and Ji
2019; Zhang et al. 2019). Generally, these HGPNNs consist of two components:
the GNN backbone which is used to embed the graph nodes and local structures,
and the pooling operation which represents graph structure in a hierarchical
way. These HGPNNs have demonstrated the necessity of adding hierarchical
pooling operations in GNNs to better preserve the graph hierarchical
structure.

However, a critical limitation of the existing hierarchical graph pooling
(HGP) strategies is that few of the pooling operations in the models are
interpretable. In many real applications, it is desirable to have an
interpretable model, where human can understand the cause of a decision made
by the model (Miller 2019; Molnar 2020). Moreover, an interpretable model is
more robust under adversarial attacks (Dai et al. 2018; Zügner, Akbarnejad,
and Günnemann 2018; Zügner and Günnemann 2020; Tang et al. 2020). A few of
recent works (Ying et al. 2019; Hou et al. 2019; Yuan et al. 2020) interpreted
the node feature embedding via GNN as a neighborhood aggregation scheme.
Particularly, they stated that the GNN embed the local feature of each node
$v$ within two steps: (1) neighbor node features aggregation and (2) node
feature transformation. However, the interpretability of pooling operations is
still not well solved. (Details are discussed in the Related Work). In order
to make the HGP operation interpretable, three questions should be considered:
Q1: How to capture the graph hierarchical structures in an interpretable way?
Q2: How to scale down the graph representation while preserving the structures
via an interpretable process?
Q3: What do we obtain after the pooling operation?

To address these challenges, we propose a Community-Based HGP framework,
COMMUNITY-POOL or CommPOOL. We aim to encode the hierarchical community
structure in graphs, which is a natural structure in many graphs, where nodes
within each community are more densely connected than the nodes across
different communities. Specifically, we propose a community-based hierarchical
pooling operation which aggregates and synthesizes the node features based on
the detected communities, such that the community structure of graphs can be
preserved during the pooling process. Moreover, we introduce a GNN-based
framework with the proposed community-based hierarchical pooling operation for
learning latent graph representations, where both local node features and the
hierarchical community-structure information are encoded and preserved. Our
contributions here can be summarized as:

  * •

We propose a community-based HGP framework (CommPOOL) for learning graph
representation in a hierarchical way that can preserve both the local node
features and the hierarchical community structure of graphs.

  * •

The proposed hierarchical community pooling strategy relies on the community
structure which is explicitly detected from the graphs, therefore the pooling
operation can capture the intrinsic community-level latent representation of
graphs and the pooling process is inherently interpretable.

  * •

We evaluate our CommPOOL framework for the whole graph classification task on
multiple public benchmark datasets. The results demonstrate the superior
performance of our model compared to several state-of-the-art graph pooling
neural networks.

  * •

Evaluations on synthetic graphs with community ground-truth labels show that
our proposed CommPOOL can capture and preserve the intrinsic community
structure of graphs during the learning process.

## Related Work

### Graph Pooling

Graph pooling operation is a strategy aiming to scale down the size of input
graphs. It can not only help to avoid model overfitting and reduce the
computational cost but also generate graph-level representations (Wu et al.
2020). In the early works (Henaff, Bruna, and LeCun 2015; Levie et al. 2018;
Dhillon, Guan, and Kulis 2007; Vinyals, Bengio, and Kudlur 2015), the graph
pooling methods simply compute the mean/max/sum of all graph node features as
the representation of the whole graph. Such a primitive pooling strategy is
named as global pooling. Later on, a few advanced techniques (e.g. attention
mechanisms (Li et al. 2015; Gilmer et al. 2017; Tran, Navarin, and Sperduti
2018), feature sorted (Zhang et al. 2018)) are proposed to improve the
performance and efficiency of the global pooling. However, the global pooling
methods do not learn the hierarchical representations, which are crucial for
capturing the structural information of graphs. Therefore, HGPNNs are
proposed.

### Interpretability of HGPNNs

Most HGP operations in the current HGPNNs (Zhang et al. 2019; Ying et al.
2018; Lee, Lee, and Kang 2019; Gao and Ji 2019; Kefato and Girdzijauskas 2020;
Bianchi, Grattarola, and Alippi 2020) show little interpretability and are
difficult to be understood by the users. Moreover, very few studies present
the interpretabiliy of their HGP operations in the paper, which may be
accounted for by the following two issues: (1). Hardly any clear definition or
analysis can be found to explain what is the captured graph structure.
Therefore, the model users are lack of heuristic knowledge to understand the
pooling operation. (2). Although some studies (Ying et al. 2018) present the
visualization of hierarchical clusters captured by the model, no quantitative
analysis is provided to examine whether the captured clusters of nodes are
aligned with the intrinsic clusters in the original graph. Apart from these,
most recent studies unfold the HGP operation as a neural network layer with
trainable parameters. The black-box nature of neural networks may also raise
extra difficulties to interpret the models in a way.

## Preliminaries

### Graph Notation

We consider the graph classification problem on attributed graphs with
different numbers of nodes. Let $G=(A,H)$ be any of the attributed graph with
$N$ nodes, where $A\in\\{0,1\\}^{N\times N}$ is the graph adjacency matrix and
$H\in\mathcal{R}^{N\times d}$ is the node feature matrix assuming that each
node has $d$ features. Also, $Z=[Z_{1},...,Z_{N}]^{T}$ is defined as the node
latent feature matrix where $Z_{i}$ is the latent feature vector for the node
$i$. Given a set of labeled data
$\mathcal{D}=\\{(G_{1},y_{1}),(G_{2},y_{2}),(G_{3},y_{3}),...\\}$ where
$y_{i}\in\mathcal{Y}$ is the classification label to the corresponding graph
$G_{i}\in\mathcal{G}$. The graph classification task can be formulated as
learning a mapping, $f$: $\mathcal{G}\rightarrow\mathcal{Y}$.

### Graph Neural Network

Graph Neural Network (GNN) is an effective message-passing architecture for
embedding the graph nodes and their local structures. Generally, GNN can be
formulated as:

$Z^{(k)}=F(A^{(k-1)},Z^{(k-1)};\theta^{(k)}),$ (1)

where $k$ denotes the layer $k$ of GNN. $A^{(k-1)}$ is the graph adjacency
matrix computed by layer $(k-1)$ of the GNN. $\theta^{(k)}$ is the trainable
parameters in the layer $k$ of the GNN. Particularly, $Z^{0}=H$.

$F(\cdot)$ is the forward function to combine and transform the messages
across the nodes. Many different versions of forward functions $F(\cdot)$ are
proposed in the previous studies (Gilmer et al. 2017; Hamilton, Ying, and
Leskovec 2017) such as Graph Convolutional Neural Network (GCN) (Kipf and
Welling 2016a) and Graph Attention Network (GAT) (Veličković et al. 2017). The
GCN linearly combines the neighborhoods as the node the representation. And
the GAT computes node representations in entire neighborhoods based on
attention mechanisms (Bruna et al. 2013).

Figure 2: Framework of the CommPOOL for graph classification. (A) is the
$1^{st}$ Embedding-Pooling (EP) module and (B) is the $2^{nd}$ EP module. In
each module, we embed the graph into the latent space by using VGAE. In the
latent space, we scale down the graph representation based on the detected
communities. (C) is the MLP for graph classification.

## The Proposed Framework

### Model Architecture

Our goal is to provide a general graph pooling framework that can capture and
preserve the hierarchical community structure of graphs in the representation
learning process of GNNs. The framework should be interpretable and it should
be able to facilitate further graph-level learning tasks, for example, graph
classification. To achieve this goal, we propose a community-based
hierarchical graph pooling (HGP) framework: CommPOOL, which is composed of $k$
cascaded Embedding-Pooling (EP) modules to learn the graph representation in a
hierarchical way. Each EP module consists of (1) an Embedding stage, where a
GNN model is employed to get the latent node representations (i.e., node
embeddings) of the input graph, and (2) a Pooling stage, where a newly
proposed community pooling mechanism is conducted on the node embeddings to
detect communities from the graph and obtain a scaled-down graph-level
representation that encodes both the local node features and the community
structure of the graph. The output of the last EP module will be the final
graph-level representation that preserves the overall hierarchical community
structure of the graph. Figure 2(A, B) shows an instance of the proposed
framework with two cascaded EP modules. In real applications of our framework,
the choice of value for $k$ is flexible and it can be decided based on the
practical needs or domain knowledge for the specific application (e.g., domain
evidence about how many community hierarchies exist in the graphs). In this
paper, we set $k=2$ and use the architecture given in Figure 2(A, B) for
illustrating our framework and we use the MLP shown in Figure 2(C) for
evaluating the CommPOOL in graph representation learning for facilitating
graph classification task.

In the following subsections, we introduce the two main parts in the proposed
EP module for CommPOOL: (1) the GNN-based Graph Node Embedding, and (2) the
Community Pooling Operation.

### GNN-based Node Embedding

We aim at a general GNN-based model to embed the graph nodes into the latent
feature space $\mathcal{Z}$ that well preserves the inherent graph structures.
On the one hand, the desired node latent features should well encode the node
information and the information between the node and its neighbors. On the
other hand, the latent features should preserve the intrinsic structures of
the graphs without task-specific influences or supervised information. On
account of the above considerations, we choose the Variational Graph Auto-
Encoders (VGAE) (Kipf and Welling 2016b) to embed nodes into the latent space
by reconstructing the graph itself.

#### Encoder

In the VGAE, we need to learn a Gaussian distribution
$q(Z|H,A)=\mathcal{N}(Z|\mu,\sigma^{2})$ which is used to approximate the
Gaussian prior $p(Z)=\mathcal{N}(Z|0,I)$. Particularly, we utilize two GNN
layers to compute the $\mu$ and $\sigma^{2}$ parameters of $q$. In the first
layer, $\mu$ and $\sigma^{2}$ share the same GNN encoder. And in the second
layer, two separate GNNs are used to generate $\mu$ and $\sigma^{2}$
respectively. The approximation can be achieved by maximizing the
Kullback–Leibler ($KL$) loss between $p$ and $q$:

$\displaystyle\mathcal{L}_{KL}=KL(q(Z|H,A)||p(Z))$ (2)

The latent features $Z$ can be obtained by resampling from the optimal
$q(Z|H,A)$.

#### Decoder

After we obtain the latent features $Z$, we reconstruct the original graph
adjacency matrix by:

$\displaystyle\hat{A}=sigmoid(ZZ^{T}).$ (3)

We define $+$ and $-$ as the edges and non-edges position index in $A$. So we
reconstruct the adjacency matrix by minimizing the $\mathcal{L}_{A}$:

$\displaystyle\mathcal{L}_{A}$ $\displaystyle=$
$\displaystyle\mathcal{L}^{+}+\mathcal{L}^{-}$ $\displaystyle=$
$\displaystyle-\frac{1}{E_{1}}\Sigma(log(\hat{A}^{+}))-\frac{1}{E_{2}}\Sigma(log(1-\hat{A}^{-}))$

where $E_{1},E_{2}$ is the number of edges and non-edges. The overall
objective function of VGAE is:

$\displaystyle\underset{Z\in\mathcal{Z}}{\text{minimize}}\quad\mathcal{L}_{A}-\mathcal{L}_{KL}$
(5)

In our CommPOOL, we use GCN (Kipf and Welling 2016a) to build up the basic
encoder layers and use GAT (Veličković et al. 2017) as the encoder variations.

### Community Pooling

#### Community Capturing

After embedding the graph nodes into the latent space $\mathcal{Z}$, we adopt
an unsupervised clustering method Partitioning Around Medoids (PAM) (Kaufmann
1987; Lee et al. 2019) on the node latent feature vectors to group the graph
nodes into $L$ different communities, where $L$ is a parameter denoting the
number of communities in the graph. Our community partition problem can be
defined as: given all the $N$ nodes in graph $G$ with their latent feature
vectors set $V=\\{Z_{1},...,Z_{N}\\}$, find $L$ different nodes with their
latent features $Z_{C}=\\{Z_{C_{1}},...,Z_{C_{L}}\\}\subset V$ from the $N$
nodes as the optimal community centers, and assign the other nodes into these
$L$ communities based on the distances between their latent feature vectors
($O=V\setminus Z_{C}$ ) and $Z_{C}$. PAM realize the community partition
problem via the following four steps.

  * •

Step 1. Initialization: Randomly select $L$ nodes with their features $Z_{C}$
as the community medoid nodes.

  * •

Step 2. Clustering: Compute the $L_{1}$ distances between the medoid nodes and
the rest nodes based on their feature vectors, and assign each non-medoid node
to its closest community; Calculate the value for the below cost function,
which computes the total distance between the non-medoid node feature vectors
$O_{j}\in O$ and their community medoid feature vectors by:

$\displaystyle Cost=\Sigma_{j=1}^{N-L}|O_{j}-Z_{C_{x}}|_{L_{1}},$ (6)

where $Z_{C_{x}}\in Z_{C}$ is the corresponding medoid of $O_{j}$.

  * •

Step 3. Adjusting: Swap each medoid node by all other non-medoids and
calculate the total cost for current configuration referring to Step 2.
Compare the cost of current and previous configuration and keep the
configuration with the smaller total cost.

  * •

Step 4. Optimization: Repeat Step 2 and 3 until the configuration does not
change.

Table 1: Average graph classification test accuracy $\pm$ standard deviation ($\%$). Models | BZR | Synthie | FRANKENSTEIN | PROTEINS | AIDS
---|---|---|---|---|---
Set2Set | $80.50\pm 1.03$ | $22.50\pm 0.86$ | $60.62\pm 0.27$ | $68.08\pm 0.56$ | $88.80\pm 0.45$
SortPool | $77.00\pm 1.24$ | $32.50\pm 1.24$ | $59.86\pm 1.22$ | $70.11\pm 0.04$ | $86.00\pm 2.42$
DIFFPOOL | $80.50\pm 1.48$ | $57.00\pm 2.62$ | $60.60\pm 1.62$ | $72.43\pm 0.26$ | $93.50\pm 1.00$
SAG-POOL | $82.00\pm 2.13$ | $45.00\pm 4.21$ | $61.73\pm 0.76$ | $71.86\pm 0.97$ | $93.50\pm 1.00$
HGP-SL | $83.00\pm 4.30$ | $54.00\pm 0.04$ | $59.51\pm 1.50$ | $\textbf{84.91}\pm\textbf{1.62}$ | $95.50\pm 1.00$
CommPOOL |  86.00 $\pm$ 1.23 |  66.50 $\pm$ 0.38 |  62.15 $\pm$ 0.37 |  $74.74$ $\pm$ $0.06$ |  98.50 $\pm$ 0.05

#### Community Pooling

In order to preserve the captured community structure during the pooling
process for the entire-graph representation learning, we propose a new pooling
mechanism called “community pooling”, which summarizes the learned node
representations based on the detected community structure. Suppose
$Z_{M_{i}}=\\{Z_{M_{i}}^{1},...,Z_{M_{i}}^{W}\\}$ is the set consisting of the
latent feature vectors of all $W$ community member nodes except for the
community center nodes $Z_{C_{i}}\in Z_{C}$ in the community$-i$. Our
community pooling problem can be defined as: given a community center feature
$Z_{C_{i}}\in Z_{C}$, and the corresponding $W$ community member features
$Z_{M_{i}}$, compute the community representation $Z_{Comm_{i}}$. The
community pooling operation computes the community$-i$’s representation by:

$\displaystyle Z_{Comm_{i}}=Z_{C_{i}}+\Sigma_{w=1}^{W}\
Sim(Z_{M_{i}}^{w},Z_{C_{i}})Z_{M_{i}}^{w},$ (7)

where $Sim(\cdot)$ is a function to measure the similarity between each member
$Z_{M_{i}}^{w}$ and the community center $Z_{C_{i}}$. In our model, we mainly
define $Sim(\cdot)$ based on $L_{1}$ distance:

$\displaystyle
Sim(Z_{M_{i}}^{w},Z_{C_{i}})=\frac{1}{\|Z_{M_{i}}^{w}-Z_{C_{i}}\|_{L1}}$ (8)

When each community representation $Z_{Comm_{i}}$ is computed, we replace the
center node feature $Z_{C{i}}$ by $Z_{Comm_{i}}$ and remove other community
member nodes. As for the graph topology structure, the preserved center nodes
are connected if and only if they are connected in the original graph. To sum
up, during the pooling, the community structure information and the node
features are preserved onto the community center nodes. And the graph
structures among the communities are presented as the topology structure of
down-scaled graph with $M<N$ nodes.

### CommPOOL for Graph Classification

When the community representations
$Z^{(K)}_{Comm}=[Z^{(K)}_{Comm_{1}},...,Z^{(K)}_{Comm_{L}}]^{T}$ are obtained
from the last Embedding-Pooling module ($k=K$), a global readout operation is
used to generate the whole graph representation $Z_{graph}$ by averaging
$Z^{(K)}_{Comm}$. Finally, an Multilayer Perceptron (MLP) utilizes $Z_{graph}$
to make predictions for graph classification. The training procedure of
CommPOOL for the graph classification task is summarized in Algorithm 1.

Input : graph: $G=(A,H)$, classification label: $y$, $K$

Output : prediction: $\hat{y}$

for  _$k=1,2,...,K$_ do

Step 1: Use $G$ to train the VGAE

Step 2: Obtain the latent feature using trained VGAE

Step 3: Community Pooling on latent features and generate down-scaled graph
$G^{(k)}=(A^{(k)},Z_{Comm}^{(k)})$. Set $G=G^{(k)}$.

end for

Step 4: $Z_{graph}=GlobalReadout(Z_{Comm}^{(K)})$

Step 5: Train $MLP$ to generate $\hat{y}=MLP(Z_{graph})$

Algorithm 1 Training Procedure

## Experiment

In this section, we evaluate our CommPOOL framework using graph classification
tasks. We present our experiment results in the following four subsections:
(1) We introduce the dataset used in the experiments. (2) We compare the graph
classification performance between CommPOOL and several competing HGPNN
models. (3) We provide some variations of the CommPOOL. (4) We test our model
on the simulation data to evaluate whether CommPOOL can accurately preserve
the community structures in the graph.

##### Dataset.

Five graph dataset are selected from the public benchmark graph data
collection (Kersting et al. 2016). Table 2 summarizes the statistics of all
dataset.

PROTEINS and Synthie (Borgwardt et al. 2005; Dobson and Doig 2003; Morris et
al. 2016) are two sets of graphs representing the protein structure. The nodes
are some amino acid features such as secondary structure content and amino
acid propensities. Nodes are linked by edges if the amino acid is an amino
acid sequence. FRANKENSTEIN (Orsini, Frasconi, and De Raedt 2015) is a set of
graphs representing the molecules with or without mutagenicity. The nodes
represent different chemical atoms and the edges are the chemical bonds type.
BZR (Sutherland, O’brien, and Weaver 2003) is a set of graphs representing the
ligands for the benzodiazepine receptor. And AIDS (Riesen and Bunke 2008) is
set of graphs representing molecular compounds with activity against HIV or
not. The molecules are converted into graphs by representing chemical atoms as
nodes and the bonds as edges.

Table 2: Dataset Statistics: $\mathcal{V}$ and the $\mathcal{E}$ represent the nodes and edges in graph $G$. $c$ represents graph classes. Dataset | $\\#|G|$ | Ave.$|\mathcal{V}|$ | Ave.$|\mathcal{E}|$ | $\\#|c|$
---|---|---|---|---
BZR | 405 | 35.75 | 38.36 | 2
Synthie | 400 | 95.00 | 172.93 | 4
FRANKENSTEIN | 4337 | 16.90 | 17.88 | 2
PROTEINS | 1113 | 39.06 | 72.82 | 2
AIDS | 2000 | 15.69 | 16.20 | 2

### Graph Classification

#### Baseline Methods

Our baseline methods include: two graph global pooling models (Set2Set
(Vinyals, Bengio, and Kudlur 2015) and SortPool (Zhang et al. 2018) ), and
three HGP models (DIFFPOOL (Ying et al. 2018), SAGPOOL (Lee et al. 2019) and
HGP-SL (Zhang et al. 2019)). For fair comparisons, we set two embedding-
pooling modules for all HGP models including three baseline HGPs and our
CommPOOL. For the baselines, we follow the original hyperparameter search
strategies provided in the related papers.

#### Experiments Setting

Following previous works (Ma et al. 2019c; Ying et al. 2018; Zhang et al.
2019), we randomly split the whole dataset into training ($80\%$) set,
validation ($10\%$) set and testing ($10\%$) set. We repeat this randomly
splitting process 10 times, and the average test performance with standard
derivation is reported in Table 1. We optimize the model via Pytorch Adam
optimizer. For the VGAE in the first module, the learning rate (lr) and the
weight decay (wd) are searched in $\\{0.0001,0.001,0.005,0.01,0.05,0.1\\}$.
The dimension of two latent GNN layers are $32$ and $16$. For the VGAE in the
second module, the lr and wd are searched in $\\{0.0001,0.001,0.005,0.01\\}$
and the dimension of two latent GNN layers are $64$ and $32$. In the community
pooling operation, the number of communities is searched in
$\\{40\%,50\%,60\%\\}$ of the number of graph nodes ($N$). The MLP consists of
two fully connected layers with $64$ and $32$ neurons and a softmax output
layer. The lr for training the MLP is searched in $\\{0.001,0.005,0.01\\}$. We
stop training if the validation loss does not decrease for 50 epochs. Our
experiments are deployed on NVIDIA Tesla P100 GPUs. We implement all baselines
and CommPOOL using PyTorch (Paszke et al. 2017) and the torch geometric
library (Fey and Lenssen 2019).

#### Summary of Results

Table 1 summarizes the classification performances of six models on five
public datasets. Our CommPOOL outperforms all baselines in the graph
classification task on almost all datasets, especially on the four-class data
Synthie. For example, our CommPOOL shows about $5.11\%$ improvement in the
classification accuracy comparing to all baselines on BZR data. This
superiority of CommPOOL may be credited to its advanced mechanism for
capturing and preserving the community structure in the pooling operation.
Also, these results indicate that the community is a crucial hierarchical
structure for learning the whole graph representation.

Moreover, Table 1 shows that hierarchical pooling methods generally perform
better than global pooling methods, which verifies that the hierarchical
pooling can better capture the graph global representations. Among all
baseline models, HGP-SL relatively performs better than others, which may be
attributed to the structure learning (SL) operations in the model. On
PROTEINS, HGP-SL performs the best among all baseline methods and even better
than ours, which implies that the structural learning strategy in HGP-SL might
be specifically suitable for PROTEINS data.

### Model Variations

To show the flexibility of CommPool, we compare several variations of CommPOOL
on PROTEINS and FRANKENSTEIN data. As noted in The Proposed Framework section,
GAT (Veličković et al. 2017) is used to replace GCN as a VGAE encoder
variation. Moreover, instead of using the reciprocal of $L_{1}$ distance, we
adopt the cosine-similarity as $Sim(\cdot)$ to measure the similarity between
community members $Z_{M_{i}}$ and the corresponding community center
$Z_{C_{i}}$ in community$-i$:

$Sim(Z^{w}_{M_{i}},Z_{C_{i}})=\frac{Z^{w}_{M_{i}}Z_{C_{i}}}{\|Z^{w}_{M_{i}}\|\|Z_{C_{i}}\|}$
(9)

The performance of CommPool with different encoders and similarity measures
are listed in Table 3, which indicates that GAT, compared to GCN, has a better
performance as the encoder in CommPOOL to embed the graph nodes. In addition,
Table 3 shows that $L_{1}$ distance is better than $cosine$ distance when
measuring the similarity between the latent features of community member nodes
and the community center nodes. A possible explanation is that $L_{1}$
distance is used in the PAM clustering. Therefore, it may be better to use the
same distance metric in the community partition process.

Table 3: Performance (%) of CommPOOL with different encoder settings and different similarity measures CommPOOL | PROTEINS | FRANKENSTEIN
---|---|---
GCN | $L_{1}$ | $74.74\pm 0.06$ | $62.15\pm 0.37$
$cosine$ | $73.84\pm 0.13$ | $60.18\pm 0.42$
GAT | $L_{1}$ | $78.84\pm 0.02$ | $63.48\pm 0.52$
$cosine$ | $76.01\pm 0.21$ | $62.32\pm 0.39$

### Community Evaluation on Simulation Data

Since no community ground-truth is provided in any publicly graph
classification datasets, we simulate a set of graphs with the known community
ground-truth and evaluate how CommPOOL preserves the intrinsic community
structures on these simulation graphs.

Table 4: Average graph classification accuracy $\pm$ standard deviation (%) on the simulation data. Models | Classification Accuracy
---|---
Set2Set | $46.54\pm 3.85$
SortPOOL | $51.29\pm 0.61$
DIFF-POOL | $67.14\pm 2.16$
SAG-POOL | N/A
HGP-SL | $72.70\pm 1.95$
CommPOOL | $80.14\pm 2.15$

##### Simulation Graphs.

We create $3$ classes of simulation graphs using different graph generating
methods, including the Random Partition Graphs, the Relax Caveman Graphs, and
the Gaussian Random Community Graphs (Brandes, Gaertler, and Wagner 2003;
Fortunato 2010). Each class contains $300$ graphs and each graph has $4$
communities with the average size of $6$ nodes. A community label is assigned
to each graph node. Meanwhile, we randomly sample from the normal distribution
$\mathcal{N}(0,I)$ as node features. We evaluate CommPOOL on the simulation
graphs to predict their class labels. Table 4 compares the graph
classification performance of CommPool with the baseline models. The results
show that, on the simulation data, the CommPOOL can also outperform all the
baseline models. N/A in Table 4 indicates the SAG-POOL cannot achieve an
optimal point in reachable epochs.

More importantly, in order to evaluate if CommPool can capture the community
structures, we compare the node community label assigned by PAM clustering in
the $1^{st}$ EP module to the community ground-truth labels. Specifically, we
compute the Normalized Mutual Information (NMI) (Strehl and Ghosh 2002)
between distribution of community labels predicted by the model and given by
the ground-truth for each graph. Figure 3a is a histogram presenting the
distribution of NMI scores for all $900$ simulation graphs. Statistically,
$79.44\%$ graphs have an NMI score larger than $0.9$ and the mean NMI score is
$0.952\pm 0.098$.

Figure 3: (a) is the NMI histogram between the distributions of community
labels of ground-truth and CommPOOL prediction. (b) is a positive example of
community structure captured by CommPOOL leading to a correct classification.
(c) is a negative example which leads to a graph misclassification. Different
colors represent different communities.

## Evaluation and Discussion

In this section, we firstly discuss the interpretibility of our proposed
community pooling operation. And then we analyze the importance of community
structure to the graph classification task.

### Interpretability of Community Pooling

CommPOOL is a hierarchical graph pooling framework with an interpretable
pooling operation. The user can transparently understand the pooling results
by monitoring the pooling operation. An interpretable pooling operation should
be capable of clearly answering three questions mentioned in the Introduction
section. Our CommPOOL provides the heuristic and knowledgeable answers for the
questions in the following way:

  * •

Q1: How to capture the graph hierarchical structures in an interpretable way?
The CommPOOL considers the communities as the basic graph hierarchical
structure. In the community pooling operation, we adopt PAM to group graph
nodes into different communities based on the distances among their features.
Such a clustering-based community capture method is transparent enough to be
interpreted. The community capture ability of our pooling operation has been
shown in the previous Community Evaluation on Simulation Data section.

  * •

Q2: How to scale down the graph while preserving the structures via an
interpretable process?
To scale down the graph, we choose the community medoid node as the
representation of the whole community, which can be understood like the
centroid can be used to represent the whole mass. Meanwhile, without loss of
necessary graph structure information, an interpretable structure preservation
process is introduced during downscaling the graph. The community pooling
achieves the preservation via gathering the nodal and structure information of
the community member nodes as the features of the community medoid node.

  * •

Q3: What do we obtain after the pooling operation?
From the graph topology view, the community pooling generates community-based
sub-graphs of the original graph since the pooling operation does not generate
any new graph nodes and edges. Each node in the sub-graph contains the
corresponding community information.

### Community Effect on Graph Classification

We design a further experiment named semi-random pooling to show that a solid
community preservation is important to the graph classification. Instead of
randomly partitioning the graph into multiple communities, we only randomly
select the community center nodes. After determining the community center
nodes, we assign each other node to the closest community based on the
similarity of node features. Such a semi-random partition method can generate
a few node cliques in graphs. These cliques, though are not the optimal
communities, can still maintain the hierarchical information to some degree.
We replace the PAM clustering by the semi-random partition in the pooling
operation. Table 5 indicates that the community pooling has significant
improvements in the graph classification tasks comparing to the semi-random
pooling, which demonstrates that the success of community capture and
preservation is crucial to the graph classification. To visualize, we select
two simulation graphs to show (1) a positive example of community structure
captured by the CommPOOL (Figure 3b); and (2) a negative example of community
structure captured by the CommPOOL, which eventually leads to the graph’s
misclassification (Figure 3c). In addition, the performance of semi-random
pooling does not decrease a lot comparing with the community pooling, which is
beyond our expectations in a way. A reasonable explanation is that although
unable to preserve the optimal community structure, the semi-random pooling
method can still capture some degree of graph hierarchical structure, which
again justifies that the significance of the community structure in the graph.

Table 5: Graph Classification Accuracy of Semi-random Pooling vs. Community-based Pooling Dataset | semi-rand. | community-based
---|---|---
PROTEINS | $64.90\pm 2.45$ | $74.74\pm 0.06$
BZR | $81.50\pm 2.82$ | $86.00\pm 1.23$
Synthie | $59.00\pm 5.89$ | $66.50\pm 0.38$
Simulation | $70.34\pm 1.26$ | $80.14\pm 2.15$

## Conclusion

In this paper, we propose CommPOOL, a new interpretable hierarchical graph
pooling framework. CommPOOL is designed for being able to capture and preserve
the inherent hierarchical community structures in graphs during the graph
representation learning and scaling-down process. Moreover, CommPOOL is a
general graph representation learning framework that can facilitate various
graph-level tasks. Experiments on both real-world graph datasets from
different domains and synthetic graph data have shown that CommPOOL
outperforms the state-of-the-art methods in graph representation learning for
the graph classification task. In future work, we will explore leveraging
CommPOOL for other graph-level tasks, such as graph regression.

## References

  * Bianchi, Grattarola, and Alippi (2020) Bianchi, F. M.; Grattarola, D.; and Alippi, C. 2020.  Spectral Clustering with Graph Neural Networks for Graph Pooling.  _arXiv: Learning_ .
  * Borgwardt et al. (2005) Borgwardt, K. M.; Ong, C. S.; Schönauer, S.; Vishwanathan, S.; Smola, A. J.; and Kriegel, H.-P. 2005.  Protein function prediction via graph kernels.  _Bioinformatics_ 21(suppl_1): i47–i56.
  * Brandes, Gaertler, and Wagner (2003) Brandes, U.; Gaertler, M.; and Wagner, D. 2003.  Experiments on graph clustering algorithms.  In _European Symposium on Algorithms_ , 568–579. Springer.
  * Bruna et al. (2013) Bruna, J.; Zaremba, W.; Szlam, A.; and LeCun, Y. 2013.  Spectral networks and locally connected networks on graphs.  _arXiv preprint arXiv:1312.6203_ .
  * Chen, Ma, and Xiao (2018) Chen, J.; Ma, T.; and Xiao, C. 2018.  Fastgcn: fast learning with graph convolutional networks via importance sampling.  _arXiv preprint arXiv:1801.10247_ .
  * Dai, Dai, and Song (2016) Dai, H.; Dai, B.; and Song, L. 2016.  Discriminative embeddings of latent variable models for structured data.  In _International conference on machine learning_ , 2702–2711.
  * Dai et al. (2018) Dai, H.; Li, H.; Tian, T.; Huang, X.; Wang, L.; Zhu, J.; and Song, L. 2018.  Adversarial attack on graph structured data.  _arXiv preprint arXiv:1806.02371_ .
  * Dhillon, Guan, and Kulis (2007) Dhillon, I. S.; Guan, Y.; and Kulis, B. 2007.  Weighted graph cuts without eigenvectors a multilevel approach.  _IEEE transactions on pattern analysis and machine intelligence_ 29(11): 1944–1957.
  * Dobson and Doig (2003) Dobson, P. D.; and Doig, A. J. 2003.  Distinguishing enzyme structures from non-enzymes without alignments.  _Journal of molecular biology_ 330(4): 771–783.
  * Duvenaud et al. (2015) Duvenaud, D. K.; Maclaurin, D.; Iparraguirre, J.; Bombarell, R.; Hirzel, T.; Aspuru-Guzik, A.; and Adams, R. P. 2015.  Convolutional networks on graphs for learning molecular fingerprints.  In _Advances in neural information processing systems_ , 2224–2232.
  * Fey and Lenssen (2019) Fey, M.; and Lenssen, J. E. 2019.  Fast graph representation learning with PyTorch Geometric.  _arXiv preprint arXiv:1903.02428_ .
  * Fortunato (2010) Fortunato, S. 2010.  Community detection in graphs.  _Physics reports_ 486(3-5): 75–174.
  * Gao and Ji (2019) Gao, H.; and Ji, S. 2019.  Graph u-nets.  _arXiv preprint arXiv:1905.05178_ .
  * Gilmer et al. (2017) Gilmer, J.; Schoenholz, S. S.; Riley, P. F.; Vinyals, O.; and Dahl, G. E. 2017.  Neural message passing for quantum chemistry.  _arXiv preprint arXiv:1704.01212_ .
  * Girvan and Newman (2002) Girvan, M.; and Newman, M. E. 2002.  Community structure in social and biological networks.  _Proceedings of the national academy of sciences_ 99(12): 7821–7826.
  * Hamilton, Ying, and Leskovec (2017) Hamilton, W.; Ying, Z.; and Leskovec, J. 2017.  Inductive representation learning on large graphs.  In _Advances in neural information processing systems_ , 1024–1034.
  * Henaff, Bruna, and LeCun (2015) Henaff, M.; Bruna, J.; and LeCun, Y. 2015.  Deep convolutional networks on graph-structured data. arXiv 2015.  _arXiv preprint arXiv:1506.05163_ .
  * Hou et al. (2019) Hou, Y.; Zhang, J.; Cheng, J.; Ma, K.; Ma, R. T.; Chen, H.; and Yang, M.-C. 2019\.  Measuring and improving the use of graph information in graph neural networks.  In _International Conference on Learning Representations_.
  * Huang et al. (2018) Huang, W.; Zhang, T.; Rong, Y.; and Huang, J. 2018.  Adaptive sampling towards fast graph representation learning.  In _Advances in neural information processing systems_ , 4558–4567.
  * Kaufmann (1987) Kaufmann, L. 1987.  Clustering by means of medoids.  In _Proc. Statistical Data Analysis Based on the L1 Norm Conference, Neuchatel, 1987_ , 405–416.
  * Kefato and Girdzijauskas (2020) Kefato, Z. T.; and Girdzijauskas, S. 2020.  Graph neighborhood attentive pooling.  _arXiv preprint arXiv:2001.10394_ .
  * Kersting et al. (2016) Kersting, K.; Kriege, N. M.; Morris, C.; Mutzel, P.; and Neumann, M. 2016.  Benchmark Data Sets for Graph Kernels.  URL http://graphkernels.cs.tu-dortmund.de.
  * Kipf and Welling (2016a) Kipf, T. N.; and Welling, M. 2016a.  Semi-supervised classification with graph convolutional networks.  _arXiv preprint arXiv:1609.02907_ .
  * Kipf and Welling (2016b) Kipf, T. N.; and Welling, M. 2016b.  Variational graph auto-encoders.  _arXiv preprint arXiv:1611.07308_ .
  * Kong and Yu (2014) Kong, X.; and Yu, P. S. 2014.  Brain network analysis: a data mining perspective.  _ACM SIGKDD Explorations Newsletter_ 15(2): 30–38.
  * Lee et al. (2019) Lee, C.-F.; Shen, J.-J.; Hou, K.-L.; and Hsu, F.-W. 2019.  A High-performance Computing Method for Photographic Mosaics upon the Hadoop Framework.  _Journal of Internet Technology_ 20(5): 1343–1358.
  * Lee, Lee, and Kang (2019) Lee, J.; Lee, I.; and Kang, J. 2019.  Self-attention graph pooling.  _arXiv preprint arXiv:1904.08082_ .
  * Levie et al. (2018) Levie, R.; Monti, F.; Bresson, X.; and Bronstein, M. M. 2018.  Cayleynets: Graph convolutional neural networks with complex rational spectral filters.  _IEEE Transactions on Signal Processing_ 67(1): 97–109.
  * Li et al. (2015) Li, Y.; Tarlow, D.; Brockschmidt, M.; and Zemel, R. 2015.  Gated graph sequence neural networks.  _arXiv preprint arXiv:1511.05493_ .
  * Liu et al. (2019) Liu, J.; Ma, G.; Jiang, F.; Lu, C.-T.; Philip, S. Y.; and Ragin, A. B. 2019.  Community-preserving graph convolutions for structural and functional joint embedding of brain networks.  In _2019 IEEE International Conference on Big Data (Big Data)_ , 1163–1168. IEEE.
  * Long et al. (2019) Long, Q.; Wang, Y.; Du, L.; Song, G.; Jin, Y.; and Lin, W. 2019.  Hierarchical community structure preserving network embedding: A subspace approach.  In _Proceedings of the 28th ACM International Conference on Information and Knowledge Management_ , 409–418.
  * Ma et al. (2019a) Ma, G.; Ahmed, N. K.; Willke, T. L.; Sengupta, D.; Cole, M. W.; Turk-Browne, N. B.; and Yu, P. S. 2019a.  Deep graph similarity learning for brain data analysis.  In _Proceedings of the 28th ACM International Conference on Information and Knowledge Management_ , 2743–2751.
  * Ma et al. (2019b) Ma, G.; Ahmed, N. K.; Willke, T. L.; and Yu, P. S. 2019b.  Deep Graph Similarity Learning: A Survey.  _arXiv preprint arXiv:1912.11615_ .
  * Ma et al. (2019c) Ma, Y.; Wang, S.; Aggarwal, C. C.; and Tang, J. 2019c.  Graph convolutional networks with eigenpooling.  In _Proceedings of the 25th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining_, 723–731.
  * Meunier et al. (2009) Meunier, D.; Lambiotte, R.; Fornito, A.; Ersche, K.; and Bullmore, E. T. 2009.  Hierarchical modularity in human brain functional networks.  _Frontiers in neuroinformatics_ 3: 37.
  * Miller (2019) Miller, T. 2019.  Explanation in artificial intelligence: Insights from the social sciences.  _Artificial Intelligence_ 267: 1–38.
  * Molnar (2020) Molnar, C. 2020.  _Interpretable Machine Learning_.  Lulu. com.
  * Morris et al. (2016) Morris, C.; Kriege, N. M.; Kersting, K.; and Mutzel, P. 2016.  Faster kernels for graphs with continuous attributes via hashing.  In _2016 IEEE 16th International Conference on Data Mining (ICDM)_ , 1095–1100. IEEE.
  * Orsini, Frasconi, and De Raedt (2015) Orsini, F.; Frasconi, P.; and De Raedt, L. 2015.  Graph invariant kernels.  In _Proceedings of the twenty-fourth international joint conference on artificial intelligence_ , volume 2015, 3756–3762. IJCAI-INT JOINT CONF ARTIF INTELL.
  * Paszke et al. (2017) Paszke, A.; Gross, S.; Chintala, S.; Chanan, G.; Yang, E.; DeVito, Z.; Lin, Z.; Desmaison, A.; Antiga, L.; and Lerer, A. 2017.  Automatic differentiation in pytorch .
  * Riesen and Bunke (2008) Riesen, K.; and Bunke, H. 2008.  IAM graph database repository for graph based pattern recognition and machine learning.  In _Joint IAPR International Workshops on Statistical Techniques in Pattern Recognition (SPR) and Structural and Syntactic Pattern Recognition (SSPR)_ , 287–297. Springer.
  * Spirin and Mirny (2003) Spirin, V.; and Mirny, L. A. 2003.  Protein complexes and functional modules in molecular networks.  _Proceedings of the national Academy of sciences_ 100(21): 12123–12128.
  * Strehl and Ghosh (2002) Strehl, A.; and Ghosh, J. 2002.  Cluster ensembles—a knowledge reuse framework for combining multiple partitions.  _Journal of machine learning research_ 3(Dec): 583–617.
  * Sutherland, O’brien, and Weaver (2003) Sutherland, J. J.; O’brien, L. A.; and Weaver, D. F. 2003.  Spline-fitting with a genetic algorithm: A method for developing classification structure- activity relationships.  _Journal of chemical information and computer sciences_ 43(6): 1906–1915.
  * Tang et al. (2020) Tang, H.; Ma, G.; Chen, Y.; Guo, L.; Wang, W.; Zeng, B.; and Zhan, L. 2020.  Adversarial Attack on Hierarchical Graph Pooling Neural Networks.  _arXiv preprint arXiv:2005.11560_ .
  * Tran, Navarin, and Sperduti (2018) Tran, D. V.; Navarin, N.; and Sperduti, A. 2018.  On filter size in graph convolutional networks.  In _2018 IEEE Symposium Series on Computational Intelligence (SSCI)_ , 1534–1541. IEEE.
  * Veličković et al. (2017) Veličković, P.; Cucurull, G.; Casanova, A.; Romero, A.; Lio, P.; and Bengio, Y. 2017.  Graph attention networks.  _arXiv preprint arXiv:1710.10903_ .
  * Vinyals, Bengio, and Kudlur (2015) Vinyals, O.; Bengio, S.; and Kudlur, M. 2015.  Order matters: Sequence to sequence for sets.  _arXiv preprint arXiv:1511.06391_ .
  * Wu et al. (2020) Wu, Z.; Pan, S.; Chen, F.; Long, G.; Zhang, C.; and Yu, S. P. 2020.  A comprehensive survey on graph neural networks.  _IEEE Transactions on Neural Networks and Learning Systems_ .
  * Ying et al. (2019) Ying, Z.; Bourgeois, D.; You, J.; Zitnik, M.; and Leskovec, J. 2019.  Gnnexplainer: Generating explanations for graph neural networks.  In _Advances in neural information processing systems_ , 9244–9255.
  * Ying et al. (2018) Ying, Z.; You, J.; Morris, C.; Ren, X.; Hamilton, W.; and Leskovec, J. 2018.  Hierarchical graph representation learning with differentiable pooling.  In _Advances in neural information processing systems_ , 4800–4810.
  * Yuan et al. (2020) Yuan, H.; Tang, J.; Hu, X.; and Ji, S. 2020.  XGNN: Towards Model-Level Explanations of Graph Neural Networks.  _arXiv preprint arXiv:2006.02587_ .
  * Zhang et al. (2018) Zhang, M.; Cui, Z.; Neumann, M.; and Chen, Y. 2018.  An end-to-end deep learning architecture for graph classification.  In _Thirty-Second AAAI Conference on Artificial Intelligence_.
  * Zhang et al. (2019) Zhang, Z.; Bu, J.; Ester, M.; Zhang, J.; Yao, C.; Yu, Z.; and Wang, C. 2019.  Hierarchical graph pooling with structure learning.  _arXiv preprint arXiv:1911.05954_ .
  * Zügner, Akbarnejad, and Günnemann (2018) Zügner, D.; Akbarnejad, A.; and Günnemann, S. 2018.  Adversarial attacks on neural networks for graph data.  In _Proceedings of the 24th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining_, 2847–2856.
  * Zügner and Günnemann (2020) Zügner, D.; and Günnemann, S. 2020.  Certifiable Robustness of Graph Convolutional Networks under Structure Perturbations.  In _Proceedings of the 26th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining_, 1656–1665.