Present address: D-Wave Systems, Burnaby, British Columbia, Canada]

# Zeeman-driven parity transitions in an Andreev quantum dot

A. M. Whiticar  [  A. Fornieri  Center for Quantum Devices, Niels Bohr
Institute, University of Copenhagen and Microsoft Quantum Lab–Copenhagen,
Universitetsparken 5, 2100 Copenhagen, Denmark  A. Banerjee  Center for
Quantum Devices, Niels Bohr Institute, University of Copenhagen and Microsoft
Quantum Lab–Copenhagen, Universitetsparken 5, 2100 Copenhagen, Denmark  A. C.
C. Drachmann  Center for Quantum Devices, Niels Bohr Institute, University of
Copenhagen and Microsoft Quantum Lab–Copenhagen, Universitetsparken 5, 2100
Copenhagen, Denmark  S. Gronin  Department of Physics and Astronomy and
Microsoft Quantum Lab–Purdue, Purdue University, West Lafayette, Indiana 47907
USA  Birck Nanotechnology Center, Purdue University, West Lafayette, Indiana
47907 USA  G. C. Gardner  Department of Physics and Astronomy and Microsoft
Quantum Lab–Purdue, Purdue University, West Lafayette, Indiana 47907 USA
Birck Nanotechnology Center, Purdue University, West Lafayette, Indiana 47907
USA  T. Lindemann  Department of Physics and Astronomy and Microsoft Quantum
Lab–Purdue, Purdue University, West Lafayette, Indiana 47907 USA  Birck
Nanotechnology Center, Purdue University, West Lafayette, Indiana 47907 USA
M. J. Manfra  Department of Physics and Astronomy and Microsoft Quantum
Lab–Purdue, Purdue University, West Lafayette, Indiana 47907 USA  Birck
Nanotechnology Center, Purdue University, West Lafayette, Indiana 47907 USA
School of Materials Engineering, Purdue University, West Lafayette, Indiana
47907 USA  School of Electrical and Computer Engineering, Purdue University,
West Lafayette, Indiana 47907 USA  C. M. Marcus<EMAIL_ADDRESS>Center for
Quantum Devices, Niels Bohr Institute, University of Copenhagen and Microsoft
Quantum Lab–Copenhagen, Universitetsparken 5, 2100 Copenhagen, Denmark

###### Abstract

The Andreev spectrum of a quantum dot embedded in a hybrid semiconductor-
superconductor interferometer can be modulated by electrostatic gating,
magnetic flux through the interferometer, and Zeeman splitting from in-plane
magnetic field. We demonstrate parity transitions in the embedded quantum dot
system, and show that the Zeeman-driven transition is accompanied by a 0-$\pi$
transition in the superconducting phase across the dot. We further demonstrate
that flux through the interferometer modulates both dot parity and 0-$\pi$
transitions.

##  I Introduction

The interplay of confinement, spin, and superconductivity leads to a rich
variety of mesoscopic phenomena Eschrig (2011); De Franceschi  _et al._
(2010); Klapwijk (2004) that can be investigated in semiconductor-
superconductor hybrid materials coupled via the proximity effect Klapwijk
(2004); Ryazanov  _et al._ (2001). Recent advances in epitaxial growth of such
hybrids have demonstrated highly transparent heterointerfaces in several
material platforms Krogstrup  _et al._ (2015); Shabani  _et al._ (2016);
Kjaergaard  _et al._ (2017); Lutchyn  _et al._ (2018).

An important application is semiconducting Josephson junctions (JJs), where a
semiconducting normal (N) region is bounded by two superconductors (S), giving
rise to a spectrum of Andreev bound states (ABSs) in the N region at energies
below the gap, $\Delta$, of the superconductors Beenakker (1991). ABS energies
$E$ depend on the superconducting phase difference, $\varphi$, across the
junction and generate a supercurrent, $I_{\rm
s}(\varphi)=-(2e/h)\,dE/d\varphi$, where $e$ is the unit of charge and $h$ is
Planck’s constant Yokoyama  _et al._ (2014); Nichele  _et al._ (2020).
Semiconducting S-N-S junctions have been used as voltage-controlled transmons
(gatemons) de Lange  _et al._ (2015); Larsen  _et al._ (2015); Casparis  _et
al._ (2018); Kringhøj  _et al._ (2020); Bargerbos  _et al._ (2020) and Andreev
qubits Zazunov  _et al._ (2003); Janvier  _et al._ (2015); Hays  _et al._
(2018); Tosi  _et al._ (2019).

A quantum dot (QD) embedded in a Josephson junction (S-QD-S) can result in a
competition between superconductivity and spin in a confined system De
Franceschi  _et al._ (2010); Eichler  _et al._ (2007); Sand-Jespersen  _et
al._ (2007); Kiršanskas  _et al._ (2015); He  _et al._ (2020). The charging
energy of a weakly coupled QD typically stabilizes one of two spin states at
zero magnetic field, depending on dot occupancy: a spin-zero singlet $\ket{S}$
or a spin-$\frac{1}{2}$ doublet $\ket{D}$Kiršanskas  _et al._ (2015); Žitko
_et al._ (2015); Meng  _et al._ (2009). For even dot occupancy, the ground
state (GS) is typically a singlet for all coupling strengths; for odd
occupancy, the GS is either a $\ket{D}$ state for weak coupling, or a
delocalized singlet, where the spin of the dot is hybridized with spins in the
leads Žitko  _et al._ (2015). The hybridized odd-parity subgap spectrum
corresponds to Yu-Shiba-Rusinov states Jellinggaard  _et al._ (2016); Žitko
_et al._ (2015); Delagrange  _et al._ (2016), and the crossover between even
and odd parity Kiršanskas  _et al._ (2015); Meng  _et al._ (2009); Vecino  _et
al._ (2003); Rozhkov and Arovas (1999) is marked by a zero-energy crossing
Chang  _et al._ (2013); Lee  _et al._ (2014); Jellinggaard  _et al._ (2016);
van Gerven Oei  _et al._ (2017).

Bound-state (BS) excitations at energies $E_{\rm{BS}}$ correspond to the
differences between GS and first excited-state (ES) energies Lee  _et al._
(2014). When $E_{\rm{BS}}=0$, the GS and ES are degenerate and a fermionic
parity transition occurs Vecino  _et al._ (2003); Žitko  _et al._ (2015). The
even and odd GS parities can be distinguished by the phase dependence of
$E_{\rm{BS}}$. An odd GS leads to a superconducting phase difference of $\pi$
across the JJ, shifting the phase dependence of the subgap excitations by
$\pi$ that results in a negative supercurrent Rozhkov and Arovas (1999); Van
Dam  _et al._ (2006). The GS parity transition from even to odd is commonly
referred to as a $0{\text{-}}\pi$ phase transition and can be identified
spectroscopically by the observation of zero-bias crossings (ZBCs) and
$\pi$-shifted phase dependence of the subgap spectrum via tunneling of single
electrons into the junction from a weakly-coupled normal lead Kiršanskas  _et
al._ (2015); Vecino  _et al._ (2003); Pillet  _et al._ (2010); Chang  _et al._
(2013).

In this Article, we investigate an S-QD-S junction embedded in a
superconducting quantum interference device (SQUID). The device allowed
control of the superconducting phase across the QD by threading magnetic flux
through the SQUID loop, while BS energies were simultaneously measured via
tunneling spectroscopy into the QD with a third (normal) lead. The system was
fabricated from an epitaxial InAs-Al heterostructure patterned using
electrostatic gates. We investigate spectra of subgap excitations under the
influence of flux and Zeeman field across parity transitions identified by
ZBCs. We observe that the subgap spectrum acquires a $\pi$ shifted energy
dependence in the odd-parity GS, as expected for a $0{\text{-}}\pi$ phase
transition. We find that these transitions are controlled via magnetic field,
gate voltage, and superconducting phase difference, which demonstrates precise
control of Andreev states in semiconductor Josephson junctions.

Magnetic-field driven parity transitions in S-QD-S junctions have previously
been observed as zeros of a reentrant critical current Estrada Saldaña  _et
al._ (2019) or via direct spectroscopy of the QD Pillet  _et al._ (2010);
Chang  _et al._ (2013). Related measurements in N-QD-S devices showed gate
voltage and field-driven parity transition as ZBCs Lee  _et al._ (2014);
Jellinggaard  _et al._ (2016). Associated $0{\text{-}}\pi$ transitions were
detected as a supercurrent reversals when the S-QD-S junction was embedded in
a SQUID Van Dam  _et al._ (2006); Delagrange  _et al._ (2016); Li  _et al._
(2019); Razmadze  _et al._ (2020). In the present study, we spectroscopically
interrogate the ABS spectrum of a S-QD-S junction, which reveals the ability
for these parameters to work in concert with a superconducting phase
difference to cause parity transitions.

##  II Device

Figure 1: (a) False-color electron micrograph of an S-QD-S device, and (b)
device schematic. Device consists of a loop (blue) of epitaxial Al, and Ti/Au
electrostatic gates (yellow). The junction is formed by two Al leads defined
by gate voltage $V_{\rm S}$ and confined into a quantum dot (QD) by gate
voltage $V_{\rm{pg}}$. The tunnel barriers to the normal (N) lead and
superconducting leads are controlled by gate voltage $V_{\rm{t}}$. ac+dc bias
voltage $V_{\rm{sd}}$ is applied to the normal lead with the superconducting
loop grounded. Magnetic field directions $B_{\parallel}$ and $B_{\perp}$ are
shown, where $B_{\perp}$ is used to apply magnetic flux through the
superconducting loop.

Devices were fabricated from an InAs two-dimensional electron gas (2DEG)
heterostructure grown on InP with 8 nm of epitaxial Al deposited _in-situ_.
Details of the heterostructure stack and device fabrication are given in the
Appendix. Previous measurements on similar material revealed near-unity
transmission of the ABS in an S-N-S JJ Kjaergaard  _et al._ (2017); Nichele
_et al._ (2020).

We study two lithographically similar devices A and B. Figure 1 shows a
micrograph of the device A. The superconducting loop was selectively wet
etched from the Al film. A 15 nm HfO2 dielectric layer, grown by atomic layer
deposition, was then deposited over the entire device. Ti/Au top-gates,
patterned by electron-beam lithography, were then evaporated. The two S leads
were defined by a negative gate voltage $V_{\rm S}$, forming a ballistic JJ of
length 200 nm and connected by an Al loop to form a single-junction SQUID. The
QD between the two S leads was defined and controlled with a negative voltage
$V_{\rm{pg}}$. Typical QD charging energies varied between $U$ = 0.7 to 1 meV,
giving $U/\Delta\sim 4$ (see Supplementary Fig. S.2). Tunnel barriers to the N
and S leads were controlled by gate voltage $V_{\rm{t}}$ applied to both
barriers, providing tunneling spectroscopy of the S-QD-S junction. A voltage
bias $V_{\rm{sd}}$ consisting of ac and dc components was applied to the
normal semiconductor (N) lead and the resulting current $I$ and four-terminal
voltage $V_{\rm 4T}$ was measured using conventional lock-in techniques with
the S loop grounded. The in-plane magnetic fields $B_{\parallel}$ and
perpendicular field $B_{\perp}$ were applied using a three-axis vector magnet.
The superconducting phase difference $\varphi$ across the junction was
controlled by threading magnetic flux through the S loop, which has area of
1.8 $\mu\rm{m}^{2}$, so that $\sim 1.2$ mT corresponds to one flux quantum,
$\Phi_{0}=h/2e=2.07$ mT$\,\mu$m2.

##  III Model

Figure 2: Subgap excitation spectrum of an S-QD-S Josephson junction (JJ) as a
function of detuning, $\varepsilon$, centered on odd occupancy at
$\varepsilon=0$. (a) Bound-state energy $E_{\rm{BS}}(\varepsilon)$ normalized
by gap, $\Delta$ for strong ($g>1$, dashed) and weak ($g<1$, solid) coupling
to the superconductors. Grey shaded region indicates odd-parity ground state.
(b) $E_{\rm{BS}}$ as a function of superconducting phase difference $\varphi$
for an even (dashed) and odd (solid) parity ground state. (c) Dependence of
$E_{\rm{BS}}$ on Zeeman energy $E_{\rm{Z}}=|g^{*}|\mu_{\rm B}B_{\parallel}$
for an even-parity ground state and for $\varphi=0,\pi$ (solid and dashed
lines, respectively). A ground-state parity transition, from even to odd
parity, occurs when $E_{\rm{BS}}=0$ for a critical Zeeman energy $E_{\rm
Z,c}$. (d) $E_{\rm{BS}}(\varphi)$ for an intermediate Zeeman energy ($E_{\rm
Z,c}(\pi)<E_{\rm{Z}}<E_{\rm Z,c}(0)$) where two zero-energy crossings occur
near $\varphi=\pi$. Excitations are calculated based on a model introduced in
Ref. Kiršanskas  _et al._ (2015) with an asymmetric lead coupling of
$\theta\sim\pi/3$. Blue/red denote spin resolved subgap states. Black denotes
spin degenerate states.

Before presenting experimental results, we first discuss the expected
dependence of bound-state energies, $E_{\rm{BS}}$, of the QD on level
detuning, $\varepsilon$, normalized Zeeman energy, $E_{\rm Z}/\Delta$, and
phase difference, $\varphi$, across the dot, including parity and
$0{\text{-}}\pi$ transitions (see Model Details for further information).
Figure 2(a) shows that for an odd-occupied QD, $E_{\rm{BS}}$ is lowered as
$\varepsilon$ is tuned away from $\pm 1$, the QD charge degeneracy points. At
the particle-hole symmetry point, $\varepsilon=0$, the bound state energy
depends on the parameter $g$, which describes the exchange coupling strength
between the spin impurity and the S-lead as Kiršanskas  _et al._ (2015); Žitko
_et al._ (2015),

$E_{\rm{BS}}=\Delta\frac{1-g^{2}}{1+g^{2}}\quad.$ (1)

For strong coupling ($g>1$, dashed curves in Fig. 2a), $E_{\rm{BS}}$ does not
reach zero for any value $\varepsilon$, preserving an even-parity GS.
Physically, this is a consequence of strong screening of the unpaired spin of
the QD by a quasiparticle in the S leads, which together form a delocalized
singlet Žitko  _et al._ (2015). An increase in Coulomb interaction reduces $g$
that leads to a lowering of $E_{\rm{BS}}$. For $g<1$, $E_{\rm{BS}}$ crosses
zero energy, signalling a transition to an odd GS parity within the shaded
region of Fig. 2 Žitko  _et al._ (2015).

Phase dependence of even and odd-parity GSs are shown in Fig. 2b. Dashed
curves show the phase dependence for the even-parity GS (a zero-junction,
denoted 0-JJ), showing a $2\pi$ periodicity with $E_{\rm{BS}}$ minima at
$\varphi=\pi$. For the odd-parity GS, the phase dependence acquires a $\pi$
shift, with energy minima at $\varphi=0$ (solid curves), yielding a $\pi$
junction (denoted $\pi$-JJ). Notably, the $0{\text{-}}\pi$ transitions can be
identified from this characteristic phase dependence of $E_{\rm{BS}}$.

Zeeman coupling, for instance from an in-plane magnetic field, can induce
parity transitions by splitting an excited bound state doublet by the Zeeman
energy, $E_{\rm Z}=|g^{*}|\mu_{\rm B}B_{\parallel}$, giving
$E_{\rm{BS}}(B_{\parallel})=E_{\rm{BS}}(0)\pm E_{\rm Z}/2$ Meden (2019); Lee
_et al._ (2014); He  _et al._ (2020), where $g^{*}$ is the effective g-factor
and $\mu_{\rm B}$ is the Bohr magneton Žitko  _et al._ (2015); Yokoyama  _et
al._ (2014); Vecino  _et al._ (2003); Jellinggaard  _et al._ (2016). In this
scenario, at a critical Zeeman energy $E_{\rm Z,c}$, a zero-energy crossing
results in a transition to an odd-parity GS, as shown in Fig. 2c Lee  _et al._
(2014); Jellinggaard  _et al._ (2016). Increasing the magnetic field further
reopens a gap that stabilizes a magnetic doublet GS with a $\pi$-shifted phase
dispersion Wentzell  _et al._ (2016). Note that Zeeman-induced parity
transitions occur at reduced $E_{\rm Z}$ for nonzero $\varphi$ (dashed curves
in Fig. 2c).

Figure 2d shows the phase dependence of an even-parity GS for Zeeman coupling
in the range $E_{\rm Z,c}(\pi)<E_{\rm{Z}}<E_{\rm Z,c}(0)$. In this
intermediate range, $E_{\rm{BS}}$ is lowered such that a gap opens in the
vicinity of $\varphi=\pi$, marked by two zero-energy crossings, indicating an
odd-parity GS. The size of the resulting gap increases with increasing Zeeman
energy until $E_{\rm{Z}}$ reaches $E_{\rm Z,c}(0)$, resulting in a dispersion
with two minima, one at $\varphi=0$ and one at $\varphi=\pi$, denoted
$0^{\prime}$-JJ or $\pi^{\prime}$-JJ, depending on which minimum is deeper
Rozhkov and Arovas (1999); Kiršanskas  _et al._ (2015).

##  IV Experiment

###  IV.1 Gate voltage dependence

Tunneling spectroscopy of the S-QD-S junction was performed by creating a
tunnel barrier to the normal lead ($V_{\rm{t}}=-1.87$ V). In Fig. 3a, the
differential conductance $G=dI/dV_{\rm 4T}$ is shown as a function of bias
voltage $V_{\rm{sd}}$ and gate voltage $V_{\rm{pg}}$ used to tune the
occupancy of the QD. When $V_{\rm{pg}}$ was varied, a gap of
$V_{\rm{sd}}\sim\pm 140~{}\mu$V was observed, along with two subgap features.
The first feature occured at $V_{\rm{pg}}\sim-5.74$ V (green marker in Fig.
3a), where the gap is reduced to $\sim 100~{}\mu$V, indicating an odd QD
occupancy with an even GS parity due to a strong coupling $g$ (see dashed line
Fig 2a).

At the stronger subgap feature in Fig. 3a, around $V_{\rm{pg}}\sim$ -5.85 V,
the gap closed entirely, resulting in two ZBCs, indicated by purple markers.
Increasing $V_{\rm{t}}$ merged the ZBCs (black circle in Fig. 3b), then
removed the crossings entirely (Fig. 3c). The sequence demonstrates a gate-
voltage-induced GS parity transition of the type illustrated in Fig. 2a, where
$g$ is controlled by gate voltage $V_{\rm{t}}$. The position in $V_{\rm{pg}}$
of the two subgap features shift with $V_{\rm{t}}$ due to cross coupling.
Figure 3d shows the dependence of the ZBC on $B_{\parallel}$ at the merging
point, $V_{\rm{t}}=-1.82$ V. The splitting is roughly linear in
$B_{\parallel}$, yielding an effective g-factor $g^{*}\sim 5$.

Figure 3: Parity transition induced by voltage $V_{\rm{pg}}$ (device A).
Differential conductance $G$ as a function of bias $V_{\rm{sd}}$ and gate
voltage $V_{\rm{pg}}$ for tunnel-barrier gate voltage $V_{\rm{t}}=$ -1.87 V
(a), -1.82 V (b), and -1.8 V (c). (d) Zero-bias conductance, $G$, as a
function of in-plane magnetic field $B_{\parallel}$ and $V_{\rm{pg}}$ for
$V_{\rm{t}}=-1.82~{}V$. In panels a-c, $B_{\parallel}=\varphi=0$.

###  IV.2 Phase dependence

Figure 4:  Dependence of an even-parity ground state on superconducting phase
difference $\varphi$ (device A). (a) Differential conductance $G$ as a
function of $\varphi$ and bias voltage $V_{\rm{sd}}$ for in-plane magnetic
field (a) $B_{\parallel}=0$ T, (b) $B_{\parallel}=0.4$ T, (c)
$B_{\parallel}=0.6$ T, (d) $B_{\parallel}=0.8$ T. (e) Zero-bias $G$ as a
function of $\varphi$ and $B_{\parallel}$ with vertically offset line-cuts
shown in (f). The red marker indicates the $V_{\rm{pg}}$ position in Fig. 3
for reference.

We next examine the phase dependence of the even-parity GS at the location of
the red marker in Fig. 3b ($V_{\rm{pg}}=-5.69$ V, $V_{\rm{t}}=-1.82$ V) by
measuring the differential conductance $G$ as a function of $\varphi$ and
$V_{\rm{sd}}$, as shown in Fig. 4a. Tuning $\varphi$ from 0 to $2\pi$ lowers
$E_{\rm{BS}}$, eventually inducing a ZBC at $\varphi=\pi$. Applying an in-
plane field $B_{\parallel}$ caused a gap to open in the vicinity of
$\varphi=\pi$, as shown in Fig. 4b. This gap increased with $B_{\parallel}$,
while the gap at $\varphi=0$ decreased, as shown in Figs. 4b-d. Although the
induced superconducting gap is suppressed at $B_{\parallel}=0.8$ T (Fig. 4d),
it is evident that the minimum $E_{\rm{BS}}$ occurs at $\varphi=0$. We
interpret the gap opening at $\varphi=\pi$ as indicating $0^{\prime}$-JJ
behavior of the type illustrated in Fig. 2d.

The position of the ZBC in both $\varphi$ and $B_{\parallel}$ is captured by
measuring $G(V_{\rm{sd}}$=0), as shown in Fig. 4e. Increasing the field causes
the crossing at $\varphi=\pi$ to split and move towards $\varphi=0$, as shown
in Fig. 4f. The splitting, roughly linear at low fields, yields a g-factor
$g^{*}\sim 5$, consistent with the value found from Fig. 3d. Cuts in Fig. 4f
at $B_{\parallel}=0$ and 0.8 T show minima shifted by $\pi$, indicating a
$0{\text{-}}\pi$ transition driven by $B_{\parallel}$.

Away from the feature marked by the red dot in Fig. 3b, where the full gap is
observed ($V_{\rm{pg}}=-5.65$ V), ABSs do not cross zero-bias at $\varphi=\pi$
at zero magnetic field (see Supplemental Fig. S.1).

###  IV.3 Magnetic field dependence

Figure 5: Magnetic field dependence of an even-parity ground state (device B).
(a) Differential conductance $G$ as a function of gate voltage $V_{\rm{pg}}$
and bias voltage $V_{\rm{sd}}$ for $\varphi=B_{\parallel}=0$. (b) Conductance
$G$ at zero bias as a function of $\varphi$ and $B_{\parallel}$. $G$ as a
function of bias $V_{\rm{sd}}$ and in-plane field $B_{\parallel}$ for
superconducting phase difference (c) $\varphi=\pi$, (d) $\pi/2$, (e) $0$. The
plots in c-e are reconstructed from line cuts of Supplementary Fig. S.3 for
fixed $\varphi$ values. Dashed lines are guides to the eye.

We next investigate bound-state bias spectra as a continuous function of
$B_{\parallel}$, rather than for the discrete values of $B_{\parallel}$ shown
in Figs. 4a-d. Focusing now on device B, Fig. 5a shows a dip in $E_{\rm{BS}}$
as a function of gate voltage $V_{\rm{pg}}$ without ZBCs, measured at
$B_{\parallel}=0$. This indicates an even-parity GS and a doublet ES
throughout this range of $V_{\rm{pg}}$. Figure 5b shows that a ZBC is first
observed for $B_{\parallel}\sim 0.5$ T at a phase difference of $\pi$.
Increasing $B_{\parallel}$ further causes the ZBC to split and merge at
$\varphi=0$, similar to Fig. 4e. Comparing the phase dependence of the ZBC
between $B_{\parallel}=0.5$ T and $0.9$ T, it is clear that the position of
the ZBC in $\varphi$ is $\pi$-shifted, indicating a magnetic field induced
$0{\text{-}}\pi$ transition. We attribute the finite field needed to induce a
parity transition in device B in comparison to device A to reflect a different
coupling $g$ resulting from a different charging energy (see Fig. 8).

Figure 5c-e shows $E_{\rm{BS}}(B_{\parallel})$ for fixed $\varphi$. At
$\varphi=\pi$ (Fig. 5c), $E_{\rm{BS}}$ splits from its zero-field value of
$V_{\rm{sd}}=\pm 75~{}\mu{\rm V}$ moving linearly towards zero, crossing zero
bias at $B_{\parallel}=0.6$ T, indicating a field-driven GS parity transition.
At $\varphi=\pi/2$ (Fig. 5d) and $\varphi=0$ (Fig. 5e) the larger zero-field
splittings push the zero-bias crossing point to larger $B_{\parallel}$. A
g-factor $g^{*}\sim 4$, extracted from the slope of the lower ES, is
insensitive to $\varphi$. The dependence of the zero-bias crossing field on
$\varphi$ is consistent with expectations in Fig. 2c. Supplementary Fig. S.3
shows spectroscopy of $E_{\rm{BS}}(\varphi)$ for fixed $B_{\parallel}$,
displaying a continuous evolution of the $0{\text{-}}\pi$ transition.

###  IV.4 Odd-parity ground state

Figure 6: Phase dependence of the odd-parity ground state (device A).
Differential conductance $G$ as a function of superconducting phase difference
$\varphi$ and bias $V_{\rm{sd}}$ for (a) $B_{\parallel}=0$, (b) $0.2$ T, (c)
$B_{\parallel}=0.4$ T. Black marker indicates the $V_{\rm{pg}}$ position in
Fig. 3 for reference.

The odd-parity transition indicated by the black marker in Fig. 3b) is
investigated in Fig. 6. At zero field, $E_{\rm{BS}}(\varphi)$ is $\pi$-shifted
compared to the even-parity case (Fig. 4), with a minimum at $\varphi=0$,
indicating a $\pi$ junction. Near the odd-parity transition,
$E_{\rm{BS}}(\varphi)$ shows a reduced dependence on $\varphi$, as reported
previously Chang  _et al._ (2013); Pillet  _et al._ (2010). Increasing
$B_{\parallel}$ increases $E_{\rm{BS}}$, opening a gap for all $\varphi$ that
increases with field while retaining the $\pi$ phase shift. This behaviour is
consistent with theory (see Fig. 2) Wentzell  _et al._ (2016).

###  IV.5 Zero-bias crossings

The contribution of in-plane magnetic field, phase difference, and gate
voltage on GS parity and $0{\text{-}}\pi$ transitions is identified by
measuring $G(V_{\rm{pg}},\varphi)$ at zero-bias in device A (see Fig. 7). This
allows for the phase dependence of the ZBCs to be highlighted at specific
$V_{\rm{pg}}$ values. Figure 7a shows two distinct values of $V_{\rm{pg}}$
where ZBCs occur at $B_{\parallel}=0$ (see red and black markers).

At $V_{\rm{pg}}=-5.7$ V, ZBCs are observed at $\varphi=\pi$, marking the
position of the even-parity GS investigated in Fig. 4. We interpret the
limited range of this ZBC in $V_{\rm{pg}}$ to reflect the energy dependence of
$E_{\rm{BS}}(\varepsilon)$ illustrated in Fig 2a, with the red marker
signifying $\varepsilon=0$. For increasing magnetic field the ZBC at
$\varphi=\pi$ splits in both $\varphi$ and $V_{\rm{pg}}$, stabilizing an odd-
parity GS (see Fig. 7d).

Figure 7: Evolution of zero-bias crossings (device A). (a) Differential
conductance $G$ at zero bias as a function of phase difference $\varphi$ and
gate voltage $V_{\rm{pg}}$ for (a) $B_{\parallel}=0$. Red and black markers
show positions of line cuts in (b). Note opposite behavior along the two line
cuts, as discussed in the text. (b) Line-cuts of (a) for an even-parity (red)
and odd-parity (black) ground states. (c) Same as (a) for $B_{\parallel}=0.2$
T, (d) Same as (a) for $B_{\parallel}=0.6$ T.

For $V_{\rm{pg}}=-5.8$ V, a bright vertical band is observed that indicates
the odd-parity GS examined in Fig. 6. In Fig. 7b the phase dependence of the
two GS locations is compared that reveals a $\pi$-shifted
$E_{\rm{BS}}(\varphi)$ dependence for the odd-parity GS. Increasing
$B_{\parallel}=0.2$ T causes the odd occupancy ZBC to split in $V_{\rm{pg}}$
while retaining a $\pi$-shifted phase dependence with $g^{*}\sim 5$. The
results of Fig. 7 reveals how the combination of gate-voltage, magnetic field,
and phase difference can control subgap excitations of the system and induce
GS parity transitions.

##  V Conclusion

To summarize, we have measured the subgap spectrum of an S-QD-S Josephson
junction under the influence of gate voltage, in-plane magnetic field, and
superconducting phase difference. We found that odd QD occupancies were not
always accompanied by parity transitions or $\pi$-shifted Andreev spectrum.
However, by controlling either the coupling, magnetic field, or phase
difference, subgap excitations could be lowered to zero bias, inducing a
parity transition. Furthermore, we showed that by applying a finite phase
difference across the junction, parity transitions can occur at lower magnetic
fields.

These results may have important implications for semiconductor based
superconducting qubits Casparis  _et al._ (2018), which recently showed that
an unintentional QD resonance resulted in a suppressed charge dispersion
Kringhøj  _et al._ (2020); Bargerbos  _et al._ (2020). We demonstrate that
highly tunable QDs can be intentionally placed in the weak link that may
enable an alternative mechanism for charge noise suppression while retaining
large qubit anharmonicity. Moreover, these results introduce novel means of
manipulating the spin of the ABS that could be used for controlling Andreev
qubits Hays  _et al._ (2018); Tosi  _et al._ (2019).

Our results demonstrate both the high material quality and device design
flexibility offered by the InAs-Al heterostructure material platform. The
S-QD-S device design studied here is a promising candidate for investigating
the hybridization of a QD with Majorana zero modes in the pursuit of parity
readout of a topological qubit Karzig  _et al._ (2017).

## References

  * Eschrig (2011) M. Eschrig, Phys. Today 64, 43 (2011).
  * De Franceschi  _et al._ (2010) S. De Franceschi, L. Kouwenhoven, C. Schönenberger, and W. Wernsdorfer, Nature Nanotechnology 5, 703 (2010).
  * Klapwijk (2004) T. Klapwijk, Journal of Superconductivity 17, 593 (2004).
  * Ryazanov  _et al._ (2001) V. V. Ryazanov, V. A. Oboznov, A. Y. Rusanov, A. V. Veretennikov, A. A. Golubov, and J. Aarts, Phys. Rev. Lett. 86, 2427 (2001).
  * Krogstrup  _et al._ (2015) P. Krogstrup, N. Ziino, W. Chang, S. Albrecht, M. Madsen, E. Johnson, J. Nygård, C. Marcus, and T. Jespersen, Nature Materials 14, 400 (2015).
  * Shabani  _et al._ (2016) J. Shabani, M. Kjaergaard, H. J. Suominen, Y. Kim, F. Nichele, K. Pakrouski, T. Stankevic, R. M. Lutchyn, P. Krogstrup, R. Feidenhans’l, S. Kraemer, C. Nayak, M. Troyer, C. M. Marcus, and C. J. Palmstrøm, Phys. Rev. B 93, 155402 (2016).
  * Kjaergaard  _et al._ (2017) M. Kjaergaard, H. J. Suominen, M. P. Nowak, A. R. Akhmerov, J. Shabani, C. J. Palmstrøm, F. Nichele, and C. M. Marcus, Phys. Rev. Applied 7, 034029 (2017).
  * Lutchyn  _et al._ (2018) R. M. Lutchyn, E. P. Bakkers, L. P. Kouwenhoven, P. Krogstrup, C. M. Marcus, and Y. Oreg, Nature Reviews Materials 3, 52 (2018).
  * Beenakker (1991) C. W. J. Beenakker, Phys. Rev. Lett. 67, 3836 (1991).
  * Yokoyama  _et al._ (2014) T. Yokoyama, M. Eto, and Y. V. Nazarov, Phys. Rev. B 89, 195407 (2014).
  * Nichele  _et al._ (2020) F. Nichele, E. Portolés, A. Fornieri, A. M. Whiticar, A. C. C. Drachmann, S. Gronin, T. Wang, G. C. Gardner, C. Thomas, A. T. Hatke, M. J. Manfra, and C. M. Marcus, Phys. Rev. Lett. 124, 226801 (2020).
  * de Lange  _et al._ (2015) G. de Lange, B. van Heck, A. Bruno, D. J. van Woerkom, A. Geresdi, S. R. Plissard, E. P. A. M. Bakkers, A. R. Akhmerov, and L. DiCarlo, Phys. Rev. Lett. 115, 127002 (2015).
  * Larsen  _et al._ (2015) T. W. Larsen, K. D. Petersson, F. Kuemmeth, T. S. Jespersen, P. Krogstrup, J. Nygård, and C. M. Marcus, Phys. Rev. Lett. 115, 127001 (2015).
  * Casparis  _et al._ (2018) L. Casparis, M. R. Connolly, M. Kjaergaard, N. J. Pearson, A. Kringhøj, T. W. Larsen, F. Kuemmeth, T. Wang, C. Thomas, S. Gronin, G. C. Gardner, M. J. Manfra, C. M. Marcus, and K. D. Petersson, Nature Nanotechnology 13, 915 (2018).
  * Kringhøj  _et al._ (2020) A. Kringhøj, B. van Heck, T. W. Larsen, O. Erlandsson, D. Sabonis, P. Krogstrup, L. Casparis, K. D. Petersson, and C. M. Marcus, Phys. Rev. Lett. 124, 246803 (2020).
  * Bargerbos  _et al._ (2020) A. Bargerbos, W. Uilhoorn, C.-K. Yang, P. Krogstrup, L. P. Kouwenhoven, G. de Lange, B. van Heck, and A. Kou, Phys. Rev. Lett. 124, 246802 (2020).
  * Zazunov  _et al._ (2003) A. Zazunov, V. S. Shumeiko, E. N. Bratus’, J. Lantz, and G. Wendin, Phys. Rev. Lett. 90, 087003 (2003).
  * Janvier  _et al._ (2015) C. Janvier, L. Tosi, L. Bretheau, Ç. Girit, M. Stern, P. Bertet, P. Joyez, D. Vion, D. Esteve, M. Goffman, _et al._ , Science 349, 1199 (2015).
  * Hays  _et al._ (2018) M. Hays, G. de Lange, K. Serniak, D. J. van Woerkom, D. Bouman, P. Krogstrup, J. Nygård, A. Geresdi, and M. H. Devoret, Phys. Rev. Lett. 121, 047001 (2018).
  * Tosi  _et al._ (2019) L. Tosi, C. Metzger, M. F. Goffman, C. Urbina, H. Pothier, S. Park, A. L. Yeyati, J. Nygård, and P. Krogstrup, Phys. Rev. X 9, 011010 (2019).
  * Eichler  _et al._ (2007) A. Eichler, M. Weiss, S. Oberholzer, C. Schönenberger, A. Levy Yeyati, J. C. Cuevas, and A. Martín-Rodero, Phys. Rev. Lett. 99, 126602 (2007).
  * Sand-Jespersen  _et al._ (2007) T. Sand-Jespersen, J. Paaske, B. M. Andersen, K. Grove-Rasmussen, H. I. Jørgensen, M. Aagesen, C. B. Sørensen, P. E. Lindelof, K. Flensberg, and J. Nygård, Phys. Rev. Lett. 99, 126603 (2007).
  * Kiršanskas  _et al._ (2015) G. Kiršanskas, M. Goldstein, K. Flensberg, L. I. Glazman, and J. Paaske, Phys. Rev. B 92, 235422 (2015).
  * He  _et al._ (2020) J. He, D. Pan, G. Yang, M. Liu, J. Ying, Z. Lyu, J. Fan, X. Jing, G. Liu, B. Lu, D. E. Liu, J. Zhao, L. Lu, and F. Qu, Phys. Rev. B 102, 075121 (2020).
  * Žitko  _et al._ (2015) R. Žitko, J. S. Lim, R. López, and R. Aguado, Phys. Rev. B 91, 045441 (2015).
  * Meng  _et al._ (2009) T. Meng, S. Florens, and P. Simon, Phys. Rev. B 79, 224521 (2009).
  * Jellinggaard  _et al._ (2016) A. Jellinggaard, K. Grove-Rasmussen, M. H. Madsen, and J. Nygård, Phys. Rev. B 94, 064520 (2016).
  * Delagrange  _et al._ (2016) R. Delagrange, R. Weil, A. Kasumov, M. Ferrier, H. Bouchiat, and R. Deblock, Phys. Rev. B 93, 195437 (2016).
  * Vecino  _et al._ (2003) E. Vecino, A. Martín-Rodero, and A. L. Yeyati, Phys. Rev. B 68, 035105 (2003).
  * Rozhkov and Arovas (1999) A. V. Rozhkov and D. P. Arovas, Phys. Rev. Lett. 82, 2788 (1999).
  * Chang  _et al._ (2013) W. Chang, V. E. Manucharyan, T. S. Jespersen, J. Nygård, and C. M. Marcus, Phys. Rev. Lett. 110, 217005 (2013).
  * Lee  _et al._ (2014) E. J. Lee, X. Jiang, M. Houzet, R. Aguado, C. M. Lieber, and S. De Franceschi, Nature Nanotechnology 9, 79 (2014).
  * van Gerven Oei  _et al._ (2017) W.-V. van Gerven Oei, D. Tanasković, and R. Žitko, Phys. Rev. B 95, 085115 (2017).
  * Van Dam  _et al._ (2006) J. A. Van Dam, Y. V. Nazarov, E. P. Bakkers, S. De Franceschi, and L. P. Kouwenhoven, Nature 442, 667 (2006).
  * Pillet  _et al._ (2010) J. Pillet, C. Quay, P. Morfin, C. Bena, A. L. Yeyati, and P. Joyez, Nature Physics 6, 965 (2010).
  * Estrada Saldaña  _et al._ (2019) J. C. Estrada Saldaña, R. Žitko, J. P. Cleuziou, E. J. H. Lee, V. Zannier, D. Ercolani, L. Sorba, R. Aguado, and S. De Franceschi, Science Advances 5 (2019), 10.1126/sciadv.aav1235.
  * Li  _et al._ (2019) C. Li, B. de Ronde, J. de Boer, J. Ridderbos, F. Zwanenburg, Y. Huang, A. Golubov, and A. Brinkman, Phys. Rev. Lett. 123, 026802 (2019).
  * Razmadze  _et al._ (2020) D. Razmadze, E. C. T. O’Farrell, P. Krogstrup, and C. M. Marcus, Phys. Rev. Lett. 125, 116803 (2020).
  * Meden (2019) V. Meden, Journal of Physics: Condensed Matter (2019).
  * Wentzell  _et al._ (2016) N. Wentzell, S. Florens, T. Meng, V. Meden, and S. Andergassen, Phys. Rev. B 94, 085151 (2016).
  * Karzig  _et al._ (2017) T. Karzig, C. Knapp, R. M. Lutchyn, P. Bonderson, M. B. Hastings, C. Nayak, J. Alicea, K. Flensberg, S. Plugge, Y. Oreg, C. M. Marcus, and M. H. Freedman, Phys. Rev. B 95, 235305 (2017).

Acknowledgments This work was supported by Microsoft Corporation, the Danish
National Research Foundation, and the Villum Foundation. We thank Karsten
Flensberg, Jens Paaske, and Jens Schulenborg for useful discussions.

## Appendix

###  V.1 Wafer structure

The wafers used for fabricating the devices were grown by molecular beam
epitaxy. The material stack consists of an InP substrate with a 100-nm-thick
$\rm{In_{0.52}Al_{0.48}As}$ lattice matched buffer, a 1-$\mathrm{\mu m}$-thick
step-graded buffer realized with alloy steps from $\rm{In_{0.52}Al_{0.48}As}$
to $\rm{In_{0.89}Al_{0.11}As}$ (20 steps, 50 nm/step), a $58~{}\rm{nm}$
$\rm{In_{0.82}Al_{0.18}As}$ layer, a $4~{}\rm{nm}$ $\rm{In_{0.75}Ga_{0.25}As}$
bottom barrier, a $7~{}\rm{nm}$ InAs quantum well, a $10~{}\rm{nm}$
$\rm{In_{0.75}Ga_{0.25}As}$ top barrier, two monolayers of GaAs and a
$7~{}\rm{nm}$ film of epitaxial Al deposited in-situ without breaking the MBE
chamber vacuum.

Hall bar device geometries (where the Al was removed) were used to
characterize the two-dimensional electron gas and revealed an electron
mobility peak $\mu=43,000~{}\rm{cm^{2}V^{-1}s^{-1}}$ for an electron density
$n=8\times 10^{11}~{}\rm{cm^{-2}}$, corresponding to an electron mean free
path of $l_{\rm e}\sim 600$ nm.

###  V.2 Fabrication Details

Devices were fabricated using standard electron beam lithography and wet
etching techniques. The devices were electrically isolated using a two-step
mesa etch by first removing the top Al film with Al etchant Transene D, and
then a deep $\sim 300$ nm III-V chemical wet etch
$\rm{H_{2}O:C_{6}H_{8}O_{7}:H_{3}PO_{4}:H_{2}O_{2}}$ (220:55:3:3). In a
following lithography step, the Al film on the mesa was selectively etched
into a SQUID with Al etchant Transene D at $50^{\circ}\mathrm{C}$. A
$15~{}\rm{nm}$ thick layer of insulating $\rm{HfO_{2}}$ was grown over the
entire sample by atomic layer deposition at a temperature of
$90^{\circ}\mathrm{C}$. Finally, top gates of Ti/Au (5/25$~{}\rm{nm}$) were
deposited by electron beam evaporation and connected to bonding pads with
leads of Ti/Au (5/300$~{}\rm{nm}$).

###  V.3 Measurement Details

Electrical measurements were performed in a dilution refrigerator at a base
temperature of $20~{}\rm{mK}$. Using conventional lock-in techniques at 166
Hz, an ac excitation voltage of 3 $\mu\rm{V}$ and a variable dc bias voltage
$V_{\rm{sd}}$ was applied to the normal lead ohmic as shown in Fig. 1. The
resulting current across the device was recorded by grounding the
superconducting loop ohmic via a low-impedance current-to-voltage converter,
and the four terminal voltage was measured by an ac voltage amplifier with an
input impedance of $500~{}\rm{M\Omega}$.

###  V.4 Model Details

The energy of the subgap excitations in an S-QD-S system was theoretically
examined under the influence of QD level detuning $\varepsilon$, coupling to
the superconductor $g$, magnetic field $B_{\parallel}$, and superconducting
phase difference $\varphi$ with the model proposed by Kiršanskas, G. et al.
Kiršanskas  _et al._ (2015). This is an Anderson-type model describing a
single Coulomb blockaded QD level that is coupled to two S leads with
superconducting gaps $\Delta\exp(\pm i\varphi/2)$ in the limit of
$U\gg\Delta$. The excitations in this model are Yu-Shiba-Rosinov states
resulting from spinful odd QD occupancies.

In Fig. 2 we examine the energy of bound-state excitations $E_{\rm{BS}}$ in an
S-QD-S JJ. The bound-state energies are calculated from Kiršanskas  _et al._
(2015),

$\begin{split}&E_{\pm,\sigma}=\frac{1}{2}\sigma E_{\rm{Z}}-\frac{\sigma
c_{\pm}\Delta}{\sqrt{(1+u)^{2}+4g^{2}}}\Big{[}(1+u)(1+\chi u)\\\ &+2g^{2}\pm
2g\sqrt{g^{2}+u(1-\chi)(1+\chi u)}\Big{]}^{1/2}\end{split}$ (2)

where the following shorthand notation is used,

$\begin{split}&\chi=1-\sin^{2}(2\theta)\sin^{2}(\phi/2),\quad
u=w^{2}-g^{2},\quad c_{+}=1\\\ &c_{-}=\rm{sign}(1+\chi
u),\quad\tan(\theta)=t_{R}/t_{L}\quad.\end{split}$ (3)

The exchange scattering amplitude $g$ and the potential scattering amplitude
$w$ both depend on the position of the QD level detuning $\varepsilon$. The
spin of the bound-states $\sigma=\pm 1$ is either aligned or anti-aligned with
respect to the spin of the QD. An angle $\theta$ is introduced to account for
an asymmetry between the left and right tunnel barriers ($t_{R/L}$) to the
superconducting leads, where $\theta=\pi/4$ represents a symmetrical coupling.
In Ref. Kiršanskas  _et al._ (2015) these bound-state energies are used to
calculate the conductance with a weakly coupled normal lead (similar setup as
in Fig. 1b), where a good agreement between the simulated conductance and the
subgap spectra shown in Fig. 2 is found.

In Fig. 8(a-d) the dependence of a varying charging energy $U$ is shown. By
decreasing the charging energy, a gap at $\varphi=\pi$ opens up due to an
increased $g$. This shifts the critical Zeeman energy to higher fields (see
Fig. 8 b, d). Figure 8(e,f) shows the effect of coupling asymmetry on the
phase dispersion. Asymmetric left/right coupling can open a gap at
$\varphi=\pi$, which can be closed by symmetrizing the coupling or applying a
magnetic field.

Figure 8: Dependence of bound-state spectra on superconducting phase
difference $\varphi$ and Zeeman energy $E_{\rm{Z}}$ for (a,b) charging energy
$U$ = 2 and (c,d) $U$ = 4 for asymmetric coupling to the superconducting leads
$\theta=\pi/3$. (e,f) Dependence of asymmetric coupling $\varphi$ on phase
dispersion and Zeeman energy of the bound-state excitations.

In the model of Kiršanskas et al. Kiršanskas  _et al._ (2015), a polarized
spin approximation on the QD is employed to derive Eq. 2. Therefore, the
Zeeman energy does not influence the QD but induces spin splitting in the
superconducting leads. Experimentally we interpret the observed magnetic field
dependence to reflect Zeeman splitting of the doublet ground states as
discussed theoretically in Refs. Jellinggaard  _et al._ (2016); Žitko  _et
al._ (2015) and experimentally in Ref. Lee  _et al._ (2014). Experimentally it
is challenging to differentiate between the two Zeeman splitting mechanisms
since they contribute different g-factor values as discussed in Ref. van
Gerven Oei  _et al._ (2017). We therefore assume an effective g-factor $g^{*}$
that accounts for a contribution from both mechanisms.

##  VI Supplementary Information

Figure S.1: Dependence of a even-parity ground state on a superconducting
phase difference $\varphi$ for intermediate values of $V_{\rm{pg}}$ in device
A. (a-d) Differential conductance $G$ as a function of $\varphi$ and bias
voltage $V_{\rm{sd}}$ for magnetic field $B_{\parallel}=0$ T (a), 0.2 T (b),
0.4 T (c), 0.6 T (d) for $V_{\rm{pg}}=-5.73$ V. (e) $G$ as a function of
$V_{\rm{sd}}$ and $V_{\rm{pg}}$ for tunnel barrier gate voltage $V_{\rm{t}}=$
-1.82 V. (f-i) $G$ as a function of $\varphi$ and $V_{\rm{sd}}$ for
$B_{\parallel}=0$ T (f), 0.4 T (g), 0.6 T (h), 0.8 T (i) for
$V_{\rm{pg}}=-5.65$ V.  Figure S.2: Coulomb blockade in device A. (a-d)
Differential conductance $G$ as a function of bias voltage $V_{\rm{sd}}$ and
$V_{\rm{pg}}$ in the normal state (a,c) and in the superconducting state
(b,d). Panels (a) and (b) are measured in a similar gate configuration shown
in Fig. 3b. Panels (c,d) are measured with more negative $V_{\rm S}$ gate
voltages to allow for clearer Coulomb blockade features.  Figure S.3:
Evolution of phase dispersion in magnetic field for device B. (a-l)
Differential conductance $G$ as a function of $\varphi$ and bias voltage
$V_{\rm{sd}}$ for increasing magnetic field $B_{\parallel}$.


# Graph Classification Gaussian Processes via Spectral Features

Felix L. Opolka Department of Computer Science and Technology
University of Cambridge
Cambridge, UK  Yin-Cong Zhi  Department of Engineering Science
University of Oxford
Oxford, UK  Pietro Liò  Department of Computer Science and Technology
University of Cambridge
Cambridge, UK  Xiaowen Dong  Department of Engineering Science
University of Oxford
Oxford, UK

###### Abstract

Graph classification aims to categorise graphs based on their structure and
node attributes. In this work, we propose to tackle this task using tools from
graph signal processing by deriving spectral features, which we then use to
design two variants of Gaussian process models for graph classification. The
first variant uses spectral features based on the distribution of energy of a
node feature signal over the spectrum of the graph. We show that even such a
simple approach, having no learned parameters, can yield competitive
performance compared to strong neural network and graph kernel baselines. A
second, more sophisticated variant is designed to capture multi-scale and
localised patterns in the graph by learning spectral graph wavelet filters,
obtaining improved performance on synthetic and real-world data sets. Finally,
we show that both models produce well calibrated uncertainty estimates,
enabling reliable decision making based on the model predictions.

##  1 Introduction

Data that are collected in a network environment, hence supported by a graph
structure, have become pervasive in modern data analysis and processing tasks.
This poses the new task of graph classification, which, similar to image
classification, aims at classifying graph-structured data into different
classes. For example, representing protein structures as graphs, one may wish
to classify whether they are toxic or not; modelling information propagation
cascades on social media platforms as graphs, one may wish to detect whether
the originating post of each cascade corresponds to fake news or not;
considering urban transportation networks as graphs, one may wish to identify
whether a particular area is likely to lead to traffic congestion.

Such graph-level classification problems are non-trivial generalisation of
classical classification problems: graphs are irregular structures and
traditional techniques defined in the Euclidean domains, such as the Fourier
transform, cannot be applied directly; data collected on the nodes (or edges)
of the graph are often continuous measurements, therefore traditional tools in
graph analysis need to be adapted; finally, real-world graphs often come with
extremely large size, making scalability of the algorithms a challenge.

Recent efforts to tackle the problem of graph classification mainly fall into
two categories. First, graph kernels [Nikolentzos et al., 2021, Kriege et al.,
2020, Borgwardt et al., 2020], as traditional ways of comparing graphs or
compute distance between them, have been adapted for graph classification.
However, it is still a challenge for these methods to handle multi-dimensional
and continuous node features as well as graphs of different sizes. Second,
graph neural networks [Bruna et al., 2014, Defferrard et al., 2016, Kipf and
Welling, 2017], generalisations of neural networks to deal with graph-
structured data, can also be utilised for graph classification. In these
frameworks, a read-out function is often deployed after the neural network
layers to summarise node representations into a single graph representation,
for example using summation, averaging, or pooling [Dai et al., 2016, Duvenaud
et al., 2015, Gilmer et al., 2017, Ying et al., 2018]. This addressed the
issue of comparing graphs of different sizes; however, these architectures are
often trained with large amount of data, and the predictions are not easily
interpretable.

(a) Graph signal $\boldsymbol{\mathrm{x}}$ consisting of three concatenated
sine-waves of different frequency.

(b) Signal under view fully localised in space (original signal
$\boldsymbol{\mathrm{x}}$)

(c) Signal under view partially localised both in space and frequency (Wavelet
transform of $\boldsymbol{\mathrm{x}}$)

(d) Signal under view fully localised in frequency (Fourier transform of
$\boldsymbol{\mathrm{x}}$)

Figure 1: Visualisation of the different space-frequency resolutions of the
spatial view of a signal (b), its wavelet transform (c), and its Fourier
transform (d) for the example of concatenated sine-waves (a) on a path graph.
A path graph forms a line of nodes where each node, except for the two end
nodes, has exactly two neighbours.

We address the above limitations in this work, and our main contributions are
as follows. First, inspired by the image segmentation literature [Porter and
Canagarajah, 1996] as well as recent development in the field of graph signal
processing [Shuman et al., 2013, Ortega et al., 2018, Dong et al., 2020,
Ortega, 2022], we propose to consider multi-dimensional and continuous node
features as graph signals, and compute spectral features using the graph
Fourier transform, i.e., energy distribution of the Fourier coefficients in
different frequency bands. We then utilise these features in a Gaussian
process framework, which has the advantage of not requiring a validation set
and being more interpretable. We show that this simple method already achieves
surprisingly competitive performance compared to baselines relying on graph
kernels and graph neural networks. Second, we further derive a second method
based on the spectral graph wavelets [Hammond et al., 2011], which possesses
the additional benefits of capturing multi-scale and localisation information
and leads to further improvement in classification performance. Finally, the
proposed Gaussian process based methods allow us to quantify uncertainty
information in the graph classification results, which to our knowledge has
not been considered so far in the literature.

##  2 Preliminaries

The models proposed here will exploit spectral features of graphs within the
framework of learning with Gaussian processes. In the following, we will
therefore give a brief introduction to the techniques from graph signal
processing and Gaussian process inference that we will use in our methodology.

###  2.1 Graph Signal Processing

Graph signal processing (see Ortega [2022] for a general introduction) offers
a range of tools for analysing the spectral properties of graph signals. One
of its key tools is the generalisation of the Fourier transform to the graph
domain, as described by Shuman et al. [2013]. It breaks down a signal into
components of different frequency, giving rise to a complementary view of the
signal on the frequency domain as opposed to the spatial domain. Moreover, the
Fourier transform is useful for analysing learned transforms of signals, which
form the basis for machine learning. Spectral graph theory [Chung, 1997] has
generalised the Fourier transform for signals on the Euclidean domain to
signals on the more general graph domain. In general, the Fourier transform of
a signal on any domain is defined as the decomposition of that signal into the
basis functions of the Laplace operator. On the graph domain, the Laplace
operator is given by the symmetric and positive-semidefinite Laplacian matrix
$\boldsymbol{\mathrm{L}}=\boldsymbol{\mathrm{D}}-\boldsymbol{\mathrm{A}}$,
where $\boldsymbol{\mathrm{A}}\in\mathbb{R}^{N\times N}$ is the adjacency
matrix of the graph with $N$ nodes and
$\boldsymbol{\mathrm{D}}\in\mathbb{R}^{N\times N}$ is the diagonal degree
matrix of the graph such that $D_{ii}=\sum_{j=1}^{N}A_{ij}$. Using the
eigendecomposition of the graph Laplacian
$\boldsymbol{\mathrm{L}}=\boldsymbol{\mathrm{U}}\boldsymbol{\mathrm{\Lambda}}\boldsymbol{\mathrm{U}}^{\top}$,
we can define the Fourier transform of a graph signal
$\boldsymbol{\mathrm{x}}\in\mathbb{R}^{N}$ as
$\hat{\boldsymbol{\mathrm{x}}}=\boldsymbol{\mathrm{U}}^{\top}\boldsymbol{\mathrm{x}}$.

The original signal $\boldsymbol{\mathrm{x}}$ and its Fourier transform
$\hat{\boldsymbol{\mathrm{x}}}$ form the two extremes of a trade-off
relationship between resolution in space and resolution in frequency:
$\boldsymbol{\mathrm{x}}$ is fully localised in space but not localised in
frequency, whereas $\hat{\boldsymbol{\mathrm{x}}}$ is fully localised in
frequency but not localised in space. The wavelet transform is a convenient
technique for analysing signals with resolution in both space and frequency.
It decomposes a signal into a set of basis functions obtained by scaling and
shifting a so-called mother wavelet. Hammond et al. [2011] derive a graph
signal wavelet transform for signal $\boldsymbol{\mathrm{x}}$ using mother
wavelet $b(\lambda)$ as

$\boldsymbol{\mathrm{w}}=\boldsymbol{\mathrm{U}}b(\beta\boldsymbol{\mathrm{\Lambda}})\boldsymbol{\mathrm{U}}^{\top}\boldsymbol{\mathrm{x}},$
(1)

where $\beta$ is the scale parameter. Adjusting this scale parameter allows
examining the signal at varying frequency ranges, which correspond to
different localisation behaviour in the spatial node domain.

The relationship between the original spatial signal, Fourier transform, and
wavelet transform is best visualised for a signal on a path graph, i.e. a
graph that simply forms a line of nodes, as it resembles an interval on the
real line. We plot a signal on such a graph consisting of three sine-waves of
different frequencies in Figure 1(a) along with the spectrograms of the signal
in the spatial domain (1(b)), in the frequency domain (1(d)), and for the
wavelet transform (1(c)), which demonstrate how the wavelet transform trades
off between spatial and spectral resolution.

To analyse a signal at multiple scales simultaneously, a low-pass filter
$h(\alpha\lambda)$ with parameter $\alpha$ and multiple scaled versions of the
mother wavelet can be combined into a more complex wavelet filter function
$g_{\theta}(\lambda)=h(\alpha\lambda)+\sum_{l=1}^{L}b(\beta_{l}\lambda)$ with
the set of scale parameters $\theta=\\{\alpha,\beta_{1},\beta_{2},\ldots\\}$.

###  2.2 Gaussian Processes

Gaussian processes (see Rasmussen and Williams [2005] for a general
introduction) are stochastic processes that can be considered multivariate
normal distributions extended to infinitely many dimensions. Their properties
make them a convenient choice for prior distributions in Bayesian machine
learning models. A GP prior on a latent function
$f:\mathbb{R}^{D}\rightarrow\mathbb{R}$ is given by

$f(\boldsymbol{\mathrm{x}})\sim\mathcal{GP}\left(m(\boldsymbol{\mathrm{x}}),k_{\theta}(\boldsymbol{\mathrm{x}},\boldsymbol{\mathrm{x}}^{\prime})\right),$
(2)

where $m(\boldsymbol{\mathrm{x}})$ is the mean function of the process and
often set to 0 and
$k_{\theta}(\boldsymbol{\mathrm{x}},\boldsymbol{\mathrm{x}}^{\prime})$ is its
kernel function with a set of kernel hyperparameters $\theta$. Given input
data
$\boldsymbol{\mathrm{X}}=[\boldsymbol{\mathrm{x}}_{1},\boldsymbol{\mathrm{x}}_{2},\ldots,\boldsymbol{\mathrm{x}}_{N}]^{\top}\in\mathbb{R}^{N\times
D}$ and labels
$\boldsymbol{\mathrm{y}}=[y_{1},y_{2},\ldots,y_{N}]^{\top}\in\mathbb{R}$, the
GP model can be used for performing probabilistic inference. If the model
specifies a Gaussian likelihood, the posterior distribution
$p(\boldsymbol{\mathrm{f}}|\boldsymbol{\mathrm{y}})$ is analytically tractable
and also follows a Gaussian process distribution. Furthermore, the marginal
likelihood $p(\boldsymbol{\mathrm{y}})$ is tractable and can be maximised to
optimise the kernel hyperparameters $\theta$.

Inference in a GP model becomes prohibitive for data sets of large size $N$,
as computing the posterior distribution requires inverting an $N\times N$
covariance matrix. Sparse Gaussian processes alleviate this computational
burden by constructing a smaller pseudo-data set of $M$ so-called inducing
points, where $M\ll N$. The inducing points are chosen such that the GP
posterior they induce is similar to the actual posterior using all $N$ data
points. Titsias [2009] proposes a way of learning the inducing inputs as
variational parameters by optimising a lower bound to the marginal likelihood.

In case the model specifies a non-Gaussian likelihood such as a Categorical
likelihood—as would be the case for classification tasks—the posterior
distribution is no longer analytically tractable and needs to be approximated.
This can be achieved using a variational approximation to the posterior
distribution where the variational family is chosen to be a Gaussian process.
Hensman et al. [2015] show how a variational approximation for non-Gaussian
likelihoods can be designed for sparse GPs. The resulting model has the
benefit that the lower bound to the marginal likelihood, referred to as the
Evidence Lower Bound (ELBO), has a data term that factorises over data points
and can therefore be maximised using stochastic gradient descent [Hensman et
al., 2013].

##  3 Spectral Feature Learning for Graph Classification

In the following, we present two GP models for graph classification. The first
model is focused on simplicity, while the second model trades off simplicity
against higher expressive power. The former is based on the graph Fourier
transform of the node feature signal and is therefore referred to as FT-GP,
whereas the latter employs the spectral graph wavelet transform and is
therefore referred to as WT-GP. The key idea in common for both approaches is
that alternatives to the spatial view of the node feature graph signals, as
provided by the Fourier and wavelet transform, are better starting points for
designing graph kernels. Consequently, we hypothesise that even simple
transformations and representations of the signals under these views, mainly
focusing on how the energy of these signals are distributed in different parts
of the spectrum and/or space, are sufficiently expressive for distinguishing
attributed graphs.

Both models are used for typical graph classification tasks where we assume to
be given a set of $G$ training set graphs
$\mathcal{T}=\\{\mathcal{G}_{1},\mathcal{G}_{2},\ldots,\mathcal{G}_{G}\\}$ of
varying sizes $N_{1},N_{2},\ldots,N_{G}$, alongside their corresponding class
labels $y_{1},y_{2},\ldots y_{G}$. For each graph $i$, we are given an
adjacency matrix $\boldsymbol{\mathrm{A}}^{(i)}\in\mathbb{R}^{N_{i}\times
N_{i}}$ and $D$-dimensional node features
$\boldsymbol{\mathrm{X}}^{(i)}\in\mathbb{R}^{N_{i}\times D}$. For the
following exposition, we initially assume one-dimensional node features, i.e.
$D=1$, and describe how to address the case of general $D$ later on. The goal
of graph classification is to learn a mapping from the adjacency matrix and
node features of a graph to its class label such that it generalises to graphs
outside the training set.

###  3.1 Spectral features for graph-level prediction

(a) Signal plotted for graph A

(b) Cumulative energy function for graph A

(c) Signal plotted for graph B

(d) Cumulative energy function for graph B

Figure 2: Comparison of two graphs from the ENZYME data set Borgwardt et al.
[2005] in terms of cumulative energy function. The same node feature is
plotted for both graphs in (a) and (c) where the node colour indicates the
node feature value. Figures (b) and (d) show the corresponding cumulative
energy functions, revealing the distinct energy profiles of both graphs: in
particular, graph A has a larger low-frequency component compared to graph B.

For each graph $\mathcal{G}^{(i)}$ with adjacency matrix
$\boldsymbol{\mathrm{A}}^{(i)}\in\mathbb{R}^{N_{i}\times N_{i}}$ we can
compute the symmetrically normalised graph Laplacian matrix
$\boldsymbol{\mathrm{L}}^{(i)}=\boldsymbol{\mathrm{I}}-(\boldsymbol{\mathrm{D}}^{(i)})^{-1/2}\boldsymbol{\mathrm{A}}^{(i)}(\boldsymbol{\mathrm{D}}^{(i)})^{-1/2}$.
Its eigendecomposition is
$\boldsymbol{\mathrm{L}}^{(i)}=\boldsymbol{\mathrm{U}}^{(i)}\boldsymbol{\mathrm{\Lambda}}^{(i)}\boldsymbol{\mathrm{U}}^{(i)\top}$,
where $\boldsymbol{\mathrm{U}}^{(i)}\in\mathbb{R}^{N_{i}\times N_{i}}$ denotes
the eigenvector matrix and
$\boldsymbol{\mathrm{\Lambda}}^{(i)}=\text{diag}(\lambda_{1}^{(i)},\ldots,\lambda_{N_{i}}^{(i)})\in\mathbb{R}^{N_{i}\times
N_{i}}$ denotes the diagonal eigenvalue matrix. The spectrum of the
symmetrically normalised graph Laplacian is in $[0,2]$, i.e.
$0\leq\lambda_{j}^{(i)}\leq 2\,\forall j=1,\ldots,N_{i}$ [Shuman et al.,
2013]. Notably, this holds regardless of the size $N_{i}$ of the graph.

Using the definition of the graph Fourier transform as presented in Section
2.1, we can compute the Fourier coefficients of the node feature signal for
graph $\mathcal{G}^{(i)}$ as
$\hat{\boldsymbol{\mathrm{x}}}^{(i)}=\boldsymbol{\mathrm{U}}^{(i)\top}\boldsymbol{\mathrm{x}}^{(i)}$.
Based on the eigenvalues of the graph and the Fourier coefficients of the node
feature signal, we can define a cumulative energy function for graph $i$ as

$e^{(i)}(z)=\sum\limits_{j=1}^{N^{(i)}}\hat{x}^{(i)2}_{j}\mathds{1}_{\\{\lambda_{j}^{(i)}\leq
z\\}}.$ (3)

This function $e^{(i)}(z)$ represents the energy of the node feature signal
that is contained in the spectrum up to a frequency $z$ and it thereby encodes
both graph structure and node feature information. We hypothesise that the
energy profile formulated in this way is expressive enough to distinguish
attributed graphs of different classes. Figure 2 visualises the cumulative
energy function for two sample graphs from the ENZYMES data set Borgwardt et
al. [2005]. We find that the two cumulative energy functions differ clearly,
both as a result of different eigenvalue locations and distinct Fourier
coefficients.

We can derive a feature vector
$\boldsymbol{\mathrm{e}}^{(i)}\in\mathbb{R}^{M}$ for graph $i$ from
$e^{(i)}(z)$ by evaluating $e^{(i)}(z)$ at a sequence of $M$ evaluation points
$[h_{1},\ldots,h_{M}]$, therefore

$e^{(i)}_{m}=e^{(i)}(h_{m})=\sum\limits_{j=1}^{N^{(i)}}\hat{x}^{(i)2}_{j}\mathds{1}_{\\{\lambda_{j}^{(i)}\leq
h_{m}\\}}.$ (4)

In practice, we will often choose $[h_{1},\ldots,h_{M}]$ to be linearly spaced
on the $[0,2]$ interval of the spectrum. Crucially, this guarantees that all
graph representations
$\boldsymbol{\mathrm{e}}^{(1)},\ldots,\boldsymbol{\mathrm{e}}^{(G)}$ are of
size $M$ regardless of individual graph sizes.

The above derivation assumes one-dimensional node features ($D=1$). The graph
representations can be generalised to multi-dimensional node features ($D>1$)
by defining a cumulative energy function for each feature dimension and
concatenating their discretisations, resulting in graph representations
$\boldsymbol{\mathrm{e}}^{(i)}$ of size $M\times D$.

The resulting graph representations can now serve as input into a base kernel
of choice and we use the radial basis function (RBF) kernel
$k(\boldsymbol{\mathrm{x}},\boldsymbol{\mathrm{x}}^{\prime})=\exp(-\frac{l}{2}\lVert\boldsymbol{\mathrm{x}}-\boldsymbol{\mathrm{x}}^{\prime}\rVert_{2})$
with lengthscale $l$.

###  3.2 Spectral wavelet features for graph-level prediction

A key limitation of the FT-GP approach stems from the Fourier transform
providing full resolution in frequency but no resolution in space. As a
result, we expect the model to under-perform in correctly classifying graphs
whose class is determined by localised patterns (cf. Figure 1). This is
evaluated in more detail in experiments with synthetic data in Section 5.1. We
can partially alleviate this limitation by employing the wavelet transform to
derive kernels, allowing us to better trade off between localisation in space
and frequency in a flexible manner.

We design wavelet filters consisting of a single low-pass filter and multiple
band-pass filters to obtain filtered signals. A particular wavelet filter
offers a distinct view of the attributed graph at hand, hence using $K$ of
those wavelet filters allows us to obtain a more diverse view of the graph.
Moreover, we can seize the capability of GPs to optimise hyper-parameters to
find wavelet scales that better distinguish the classes of graphs of the
particular data set at hand. Figure 3 plots the wavelet transformed signal for
a particular graph and two filters of different scales, to showcase the
distinct view of the attributed graph each filter offers.

(a) Wavelet filter A

(b) Wavelet filter B

(c) Wavelet atom localised at red node (filter A)

(d) Wavelet atom localised at red node (filter B)

(e) Signal filtered with filter A

(f) Signal filtered with filter B

Figure 3: Comparison of two different filters for the same graph and node
feature (as a graph signal) from the ENZYME data set. The band-pass filter of
large wavelet scale (a) highlights low frequencies, captures a larger node
neighbourhood (c), and produces a smooth signal (e). The band-pass filter of
small wavelet scale (b) highlights high frequencies, captures a smaller node
neighbourhood (d), and produces a less smooth signal (f).

We define each of the $K$ wavelet filter functions as a sum of atomic filters
including a single low-pass filter and $L$ band-pass filters

$g_{\theta}(\lambda)=h(\alpha\lambda)+\sum\limits_{l=1}^{L}b(\beta_{l}\lambda),$
(5)

where $h(\alpha\lambda)$ refers to the low-pass filter function with scale
$\alpha$, $b(\beta_{l}\lambda)$ refers to the $l$-th band-pass filter with
scale $\beta_{l}$, and $\theta=\\{\alpha,\beta_{1},\beta_{2},\ldots\\}$ is the
set of scale parameters. Reminiscent of the signal energy computed for the FT-
GP model proposed earlier, we compute the magnitude of signal filtered with
filter $k$

$w^{(i)}_{k}=\lVert\boldsymbol{\mathrm{U}}^{(i)}g_{\theta_{k}}(\boldsymbol{\mathrm{\Lambda}}^{(i)})\boldsymbol{\mathrm{U}}^{(i)\top}\boldsymbol{\mathrm{x}}^{(i)}\rVert_{2}$
(6)

to express how much energy of the signal is captured by the $k$th filter. The
final feature vector $\boldsymbol{\mathrm{w}}^{(i)}\in\mathbb{R}^{K\cdot D}$
for graph $i$ stacks the filter coefficients for each of the $K$ wavelet
filters and each of the $D$ node feature signals. As before, the resulting
signal representations are fed into a radial basis function kernel to obtain
the covariance between graphs. The kernel hyper-parameters are optimised
alongside the wavelet scale parameters $\theta$ in a data-driven way employing
type II maximum likelihood estimation.

###  3.3 Relationship to graph neural networks

The GP models presented here operate in two stages, one focused on
transformations of the node feature signal and the second on deriving a graph-
level summary. In this regard, they share similarities with graph neural
networks for graph classification, making it worthwhile to compare the two
approaches more thoroughly.

Both GP methods, in the first step, translate the node feature signal from the
spatial domain to a domain that is either fully localised in frequency or
partially localised in both space and frequency. The resulting coefficients
thereby combine information from all nodes across the graph, overcoming the
spatial locality of each node. In case of WT-GP, the wavelet transform can be
considered an aggregation operation of node features in a neighbourhood around
each node whose shape is determined by the particular wavelet filter in use
[Opolka et al., 2022]. This aggregation of information across graphs is
reminiscent of the aggregation operation inherent to graph convolutional
networks (GCN) [Kipf and Welling, 2017]. In fact, Wu et al. [2019] describe
how the GCN aggregation operation acts as a low-pass filter on the node
feature graph signal.

In a second step, the GP methods aggregate coefficients across the whole graph
via binning in case of FT-GP and summation in case of WT-GP. Neural network
methods apply similar so-called “read out” operations to summarise node
representations into a single graph representation via simple summation or
averaging [Dai et al., 2016, Duvenaud et al., 2015, Gilmer et al., 2017], as
well as more complicated pooling operations [Ying et al., 2018]. In contrast
to GCN-based approaches however, the aggregation approach presented here has a
straightforward spectral interpretation and allows aggregating information in
a way that goes beyond averaging low-pass filtered values, making it possible
to capture multi-scale information in the graphs. Furthermore, the GP
framework provides the flexibility of a non-parametric approach and outputs
uncertainty estimates for predictions, which are examined in more detail in
Section 5.4.

###  3.4 Scalability

In terms of scalability of the GP models, the quantities of concern are the
number of graphs in a data set and the size of those graphs. Within the
computation of the kernel of either proposed approach, the computational
complexity is dominated by the necessary eigendecomposition of the graph
Laplacian, which is in $\mathcal{O}(N_{i}^{3})$, where $N_{i}$ is the size of
a graph $i$. Graphs in graph classification tasks are typically small (cf.
Table 2), therefore the number of nodes is not usually an obstacle for
scalability. For larger graphs, one can reduce the computational complexity by
resorting to the approximate Fourier transform of a graph signal [Le Magoarou
et al., 2018].

The total number of graphs in the data set may pose another limit to the
scalability of the method when computing the GP posterior. Sparse GPs as
discussed in Section 2.2 alleviate this issue by enabling stochastic
optimisation in mini-batches when the data set size requires it.

##  4 Related Work

Applications of Spectral and Wavelet Energy: Our work is related to building
classifiers in the graph spectral domain, and wavelets are an extension to
such approaches with the benefits of multi-scale properties and better
localisation. Solving problems on graphs by making use of the energy
distribution in the spectral domain has been demonstrated effectively for the
task of image segmentation in classical signal processing [Porter and
Canagarajah, 1996], and for geographical mobility prediction in graph signal
processing [Dong et al., 2013]. Both suggest that the energy distribution
clearly contains a significant level of information which we aim to utilise in
this work for graph classification.

Graph Kernels: One of the first studies into kernel functions acting on the
graph level is the graph kernels summarised in the surveys of Nikolentzos et
al. [2021], Kriege et al. [2020], Borgwardt et al. [2020]. Graph kernels
generally have multiple definitions, of which a number of choices have in
common the existence of a double sum of a base kernel. To compute the kernel
between two graphs, a standard base kernel is chosen and takes as inputs a
node feature from each graph, and the double summation then covers all
possible pairs of nodes between the two graphs. Similar definition can be
applied to graphs that contain edge features. Such designs are limited as any
edge connections are ignored in the double sum, and the kernel therefore boils
down to a comparison of cross-products of the node features. Another common
design is the computation of a product graph, where a new graph is constructed
based on the two input graphs. When we simulate a random walk on the product
graph, we compute the number of matching walks on the two individual graphs.
Random walk kernels do take into account the graph structure, but is limited
by scalability as the computation of the kernel over a product graph leads to
a $\mathcal{O}(N_{i}^{3}N_{j}^{3})$ complexity.

Graph kernels has also been developed from the Weisfeiler-Lehman (WL) test for
graph isomorphism [Huang and Villar, 2021], which also acts as a graph
similarity. The WL test produces a series of node “colouring” from a
neighbourhood gathering and relabelling procedure. In the WL kernel
[Shervashidze et al., 2011], the “colourings” are passed through a hashing
function (or histogram mapping) to form inputs to a base kernel. As an
alternative, the Wasserstein distance between two sets of labels can be used
instead of the hashing function as shown in Togninalli et al. [2019]. Compared
to our kernel design, WL-based kernels generally work with the neighbourhoods
in the spatial domain, and so they ignore any information in the spectral
domain. Additionally, if the graph contains high dimensional node features,
the WL does not have an efficient way to encode them into the embeddings for
the WL algorithm.

The use of Fourier and wavelets basis are Laplacian-based models. Wavelets in
particular, allow for aggregation over a continuous neighbourhood, giving the
model the ability to operate on a multi-scale level. Though multi-scale
Laplacian kernels exist such as Kondor and Pan [2016], they only operate in
the spatial domain by computing kernels between sub-graphs of various
neighbourhood sizes. On the other hand, the work of Pineau [2019] does make
use of the Laplacian spectral information, but like the previous kernel, it
ignores the node features that may come with the graphs, while we focus on the
spectral information of the node features.

Graph Neural Network Models: Lastly, graph neural networks (GNNs) have also
been applied as a test for graph isomorphism, and as a result can be applied
to graph classification. The message passing step of a GNN layer (examples
including Duvenaud et al. [2015], Li et al. [2012], Murphy et al. [2019]) is
comparable to the neighbourhood gathering step in the WL algorithm. Analysis
of certain GNN models showed that they are at best as powerful as the WL test,
and the graph isomophism network (GIN) proposed in Xu et al. [2019] achieves
the theoretical guarantee of the WL. Meanwhile, simpler models such as GCN
[Kipf and Welling, 2017] and GraphSAGE [Hamilton et al., 2017] have been shown
to be unable to distinguish certain types of graphs that GIN can handle.

##  5 Experiments

Our work aims to investigate two core empirical questions with regards to the
models presented. Firstly, whether the energy profile as captured by either of
the presented kernels is sufficiently expressive to classify real-world
graphs. Secondly, whether the wavelet transform approach of WT-GP improves
upon the FT-GP model based on the Fourier transform. We begin with the latter
question by comparing the two proposed models on synthetic data sets.

In all our experiments, we use the same experimental setup, unless explicitly
stated otherwise. The FT-GP discretises the cumulative energy function by
evaluating it at $30$ evaluation points. The WT-GP uses $10$ filters
consisting of a single low-pass and three band-pass filters. Each low-pass
scale of each filter is uniformly randomly initialised to a value between
$4.0$ and $6.0$. The band-pass scales of each filter are uniformly randomly
initialised to a value between $0.1$ and $5.0$. A non-sparse variational
Gaussian process as described in Section 2.2 is trained using the L-BFGS-B
optimiser [Liu and Nocedal, 1989, Byrd et al., 1995] until the ELBO
convergences. All results are cross-validated using 10-folds under a
stratified split where in each fold $80\%$ of the data is used for training
and $10\%$ for testing. The remaining $10\%$ are set aside for validation to
make the evaluation comparable to results in related work, although neither of
the two GP models make use of the validation set. This validation procedure is
suggested by Errica et al. [2020] to overcome weaknesses in the evaluation of
graph classification methods in previous work.

###  5.1 Synthetic experiments

(a) Class $0$

(b) Class $1$

Figure 4: Visualisation of example graphs of the two different classes in the
synthetic ring-vs-clique data set. All graphs consist of a graph component
sampled from an Erdős-Rényi (ER) model, which is connected to either a ring
graph component in case of class $0$ or a fully connected graph component in
case of class $1$ (highlighted in blue). Both components of each graph can
vary in the number of nodes.

We design two synthetic data sets to highlight the differences between the two
GP models presented here. The first data set, referred to as ring-vs-clique,
contains two classes. The classes differ in that graphs in class $0$ are
guaranteed to contain a subgraph that forms a ring of at least size $5$, while
graphs in class $1$ are guaranteed to contain a subgraph that forms a complete
graph (or clique) of at least size $5$. We generate the class-balanced data
set by sampling $200$ graphs using the following procedure. In a first step, a
graph of uniformly random size between $10$ and $30$ nodes is sampled from an
Erdős-Rényi model [Erdös and Rényi, 1959]. In the second step, depending on
the class label, a ring or clique of uniformly random size between $5$ and
$10$ nodes is constructed and connected to a randomly chosen node in the graph
sampled in the first step through a single edge. The two different classes are
sketched in Figure 4. The data set is designed to test how well a model can
distinguish graphs when class labels are determined only by a relatively
small, localised subgraph.

The second data set, referred to as sbm, contains graphs drawn from a
stochastic block model [Holland et al., 1983], where nodes of the same block
are connected with a probability of $80\%$ and nodes of different blocks are
connected with $10\%$ probability. Graphs of class $0$ are drawn from a model
with $2$ blocks and graphs of class $1$ are drawn from a model with $3$
blocks. All graphs have a uniform random size between $10$ and $30$ nodes.

The results for the binary classification tasks on these data sets for the two
GP models are shown in Table 1. We find that on both tasks, the wavelet-based
GP model strongly outperforms the baseline FT-GP model. In fact, FT-GP
performs only slightly better than random guessing on both data sets whereas
WT-GP achieves a near perfect accuracy on ring-vs-clique and a high accuracy
on the sbm data set. The results confirm that designing graph kernels based on
wavelet filtered node feature signals leads to more expressive kernels
compared to a kernel based on Fourier transformed node features.

| ring-vs-clique | sbm
---|---|---
FT-GP |  $62.5$ $\pm 7.5$ |  $58.5$ $\pm 15.2$
WT-GP |  $\boldsymbol{\mathrm{99.5}}$ $\pm 1.5$ |  $\boldsymbol{\mathrm{91.0}}$ $\pm 5.4$
Table 1: Classification accuracy of the FT-GP model compared to the WT-GP
model on the two synthetic binary classification data sets.

Data set | # Graphs | # Classes | Avg # Nodes | Avg # Edges | # Node Attr
---|---|---|---|---|---
ENZYMES | $600$ | $6$ | $32.63$ | $62.14$ | $21$
MUTAG | $188$ | $2$ | $17.93$ | $19.79$ | $7$
NCI1 | $4{,}110$ | $2$ | $29.87$ | $32.30$ | $37$
IMDB-BIN | $1{,}000$ | $2$ | $19.77$ | $96.53$ | $-$
IMDB-MUL | $1{,}500$ | $3$ | $13.00$ | $65.94$ | $-$

Table 2: Statistics of the data sets used in the empirical evaluation of the
proposed models. We use three data sets with and two data sets without node
features.

###  5.2 Real-world experiments

To examine whether the proposed GP models based on the energy profiles of
attributed graphs are sufficiently expressive for classifying real-world
graphs, we conduct a number of experiments on benchmark data sets and compare
the model performance to popular baseline methods for graph classification.

The data sets in our empirical evaluation are shown in Table 2 along with an
overview of their statistics. The ENZYMES data set [Borgwardt et al., 2005]
contains protein graphs of enzymes that require classification into the Enzyme
Commission top level enzyme classes. The MUTAG data set [Debnath et al., 1991]
consists of molecular graphs and the task is to detect an effect on the
Salmonella typhimurium bacterium. Similarly, the NCI1 data set [Wale and
Karypis, 2006] holds molecular graphs that need to be classified based on
whether they are active against certain types of cancer. Finally, the IMDB-
BINARY and IMDB-MULTI data sets [Yanardag and Vishwanathan, 2015] are derived
from the Internet Movie Database (imdb.com) and consist of graphs of actors
and actresses who have co-starred together in a film. A classifier has to
predict the genre of each graph belonging to a particular actor or actress.
The graphs in the first three data sets come with node attributes while the
IMDB data sets are unattributed. We add one-hot encoded node degrees as
features to ensure all graphs have node attributes.

We compare the proposed GP models to a number of neural network and graph
kernel baselines. The neural network baselines include DGCNN [Zhang et al.,
2018], GraphSAGE [Hamilton et al., 2017], DiffPool [Ying et al., 2018], and
GIN [Xu et al., 2019]. Among graph kernels we compare to the Shortest Path
(SP) kernel [Borgwardt and Kriegel, 2005], the Weisfeiler-Lehman (WL) kernel
[Shervashidze et al., 2011], and the Multiscale Laplacian (ML) kernel [Kondor
and Pan, 2016].

| ENZYMES | MUTAG | NCI1 | IMDB-BINARY | IMDB-MULTI
---|---|---|---|---|---
DGCNN |  $38.9$ $\pm 5.7$ |  $84.0$ $\pm 7.1$ |  $76.4$ $\pm 1.7$ |  $69.2$ $\pm 3.0$ |  $45.6$ $\pm 3.4$
GraphSAGE |  $58.2$ $\pm 6.0$ |  $83.6$ $\pm 9.6$ |  $76.0$ $\pm 1.8$ |  $68.8$ $\pm 4.5$ |  $47.6$ $\pm 3.5$
DiffPool |  $59.5$ $\pm 5.6$ |  $79.8$ $\pm 6.7$ |  $76.9$ $\pm 1.9$ |  $68.4$ $\pm 3.3$ |  $45.6$ $\pm 3.4$
GIN |  $\boldsymbol{\mathrm{59.6}}$ $\pm 4.5$ |  $84.7$ $\pm 6.7$ |  $\boldsymbol{\mathrm{80.0}}$ $\pm 1.4$ |  $\boldsymbol{\mathrm{71.2}}$ $\pm 3.9$ |  $\boldsymbol{\mathrm{48.5}}$ $\pm 3.3$
SP | timeout |  $82.4$ $\pm 5.5$ |  $72.5$ $\pm 2.0$ |  $58.2$ $\pm 4.7$ |  $39.2$ $\pm 2.3$
WL |  $50.7$ $\pm 7.3$ |  $\boldsymbol{\mathrm{86.7}}$ $\pm 7.3$ |  $\boldsymbol{\mathrm{85.2}}$ $\pm 2.2$ |  $70.7$ $\pm 6.8$ |  $\boldsymbol{\mathrm{51.3}}$ $\pm 4.4$
ML |  $33.2$ $\pm 5.8$ |  $\boldsymbol{\mathrm{87.2}}$ $\pm 7.5$ |  $\boldsymbol{\mathrm{79.7}}$ $\pm 1.8$ |  $69.9$ $\pm 4.8$ |  $47.7$ $\pm 3.2$
FT-GP |  $\boldsymbol{\mathrm{60.7}}$ $\pm 4.3$ |  $85.7$ $\pm 6.2$ |  $77.7$ $\pm 1.6$ |  $\boldsymbol{\mathrm{72.7}}$ $\pm 3.9$ |  $\boldsymbol{\mathrm{48.8}}$ $\pm 2.8$
WT-GP |  $\boldsymbol{\mathrm{63.8}}$ $\pm 5.3$ |  $\boldsymbol{\mathrm{87.3}}$ $\pm 4.8$ |  $78.1$ $\pm 2.1$ |  $\boldsymbol{\mathrm{74.6}}$ $\pm 4.1$ |  $48.4$ $\pm 2.9$

Table 3: Comparison of the proposed GP models with common neural network and
graph kernel baselines in terms of classification accuracy. Colours indicate
the best, second-, and third-best result.

The classification accuracy of each model on all data sets is shown in Table 3
where the results for baseline methods are obtained from the empirical
evaluation by Nikolentzos et al. [2021]. Mirroring the results on the
synthetic data set, we find that WT-GP performs at least as well as FT-GP on
all data sets except on IMDB-MULTI where both perform roughly similarly.
Moreover, both models yield competitive performance compared to the baselines,
with WT-GP achieving the highest accuracy on three of the five data sets and
FT-GP being among the best two models also on three out of five data sets. We
highlight the surprising effectiveness of the Fourier features of FT-GP in
comparison to the baselines, especially since no learning is involved in the
feature construction.

While overall outperforming the neural network baselines, the trend across
data sets is comparable between the GP models and the neural networks. In
fact, where the kernel methods outperform the GNNs, they also tend to be more
comparable to the GP models. This appears to provide evidence for the
similarity between the methods proposed here and GNNs, as outlined in Section
3.3.

###  5.3 Robustness analysis

While the proposed GP models have fewer hyper-parameters than comparable
neural network models, a small number of values needs to be selected prior to
model fitting. For the FT-GP this is primarily the number of evaluation points
$M$ and for WT-GP this is mainly the number of filters $K$.

We evaluate FT-GP for a range of different numbers of evaluation points $M$ on
two data sets, one with and one without node features. The results are shown
in Table 4. We find that FT-GP is overall robust to varying the number of
evaluation points. We note that larger $M$ linearly increase time and memory
complexity for training and inference. Notably, however, larger $M$ do not
increase the number of learned parameters of FT-GP (which does not have
learned parameters) but merely increase the “resolution” of the approximation
of the cumulative energy function and we therefore do not expect FT-GP to
begin to over-fit for larger $M$.

In a similar vein, the performance of WT-GP for varying numbers of filters $K$
is presented in Table 5. The results indicate that WT-GP is robust to
different numbers of filters. Similar to $M$, increasing $K$ linearly
increases the time and memory complexity for training and inference. Unlike
$M$, however, increasing $K$ means a larger number of scale parameters need to
be estimated using MLE and therefore we do expect the model to over-fit for
very large values of $K$. For the values that are computationally feasible for
the given data sets, we have not yet observed over-fitting.

$M$ | ENZYMES | IMDB-BINARY
---|---|---
$\boldsymbol{\mathrm{20}}$ |  $61.2$ $\pm 5.1$ |  $73.1$ $\pm 4.0$
$\boldsymbol{\mathrm{25}}$ |  $60.5$ $\pm 4.5$ |  $73.5$ $\pm 3.7$
$\boldsymbol{\mathrm{30}}$ |  $60.7$ $\pm 4.3$ |  $72.7$ $\pm 3.9$
$\boldsymbol{\mathrm{35}}$ |  $60.3$ $\pm 4.9$ |  $73.6$ $\pm 3.5$
$\boldsymbol{\mathrm{40}}$ |  $61.0$ $\pm 4.9$ |  $73.0$ $\pm 3.6$
$\boldsymbol{\mathrm{45}}$ |  $61.0$ $\pm 4.9$ |  $73.0$ $\pm 4.0$
$\boldsymbol{\mathrm{50}}$ |  $61.0$ $\pm 4.7$ |  $72.9$ $\pm 4.2$
Table 4: Performance of FT-GP on a data set with node features (ENZYMES) and a data set without node features (IMDB-BINARY) for varying number of evaluation points $M$. $K$ | ENZYMES | IMDB-BINARY
---|---|---
$\boldsymbol{\mathrm{5}}$ |  $64.3$ $\pm 5.0$ |  $74.6$ $\pm 3.8$
$\boldsymbol{\mathrm{10}}$ |  $63.8$ $\pm 5.3$ |  $74.6$ $\pm 4.1$
$\boldsymbol{\mathrm{15}}$ |  $64.8$ $\pm 5.4$ |  $74.3$ $\pm 3.0$
$\boldsymbol{\mathrm{20}}$ |  $65.5$ $\pm 5.5$ |  $74.2$ $\pm 4.1$
Table 5: Performance of WT-GP on a data set with node features (ENZYMES) and a
data set without node features (IMDB-BINARY) for varying number of filters
$K$.

###  5.4 Quality of uncertainty estimates

FT-GP

(a) MUTAG

(b) IMDB-BINARY (c) MUTAG (d) IMDB-BINARY

WT-GP

Figure 5: Accuracy of predictions when rejecting predictions that have a
variance above a certain threshold. As the threshold becomes lower, i.e. more
and more strict (from right to left on x-axis), the prediction accuracy
increases.

One of the core advantages of a Bayesian treatment of graph classification is
the availability of uncertainty estimates that can be a crucial part of the
real-world application of a model. Uncertainty estimates allow downstream
users of the model to weigh the reliability of its predictions and thus more
confidently make decisions based on those predictions. The GP models proposed
here make uncertainty predictions in the form of the variance of the
variational posterior predictive distribution. We can assess the quality of
those uncertainty estimates by simulating an experiment where downstream users
are allowed to “reject” predictions if their variance is above a certain
threshold. We can then plot the accuracy of the predictions that are not
rejected. If the uncertainty estimates are well calibrated, we expect the
accuracy of the remaining, low-variance (i.e. high certainty) predictions to
increase. We plot the results for both our models and two of the data sets in
Figure 5. We find that, as expected, the accuracy increases for stricter
variance thresholds. As the variance threshold increases, the accuracy rises
from the results reported in Table 3 for the complete set of predictions, to
$1.0$ for only the most confident predictions.

##  6 Conclusions

We have proposed two GP models using spectral features to classify graphs. The
first approach constructs the energy profile of an attributed graph and
compares across features. We find that even though the model requires no
learning to obtain these spectral features, it performs competitively to graph
neural network and kernel baselines. The second approach learns more complex
wavelet filters and compares graphs based on the corresponding filtered node
features. It outperforms the first approach both on real-world and synthetic
data sets. Our work indicates that spectral features constitute a powerful
basis for graph classification both on their own and even more so when
combined with the learning capabilities of GPs.

###### Acknowledgements.

FLO acknowledges funding from the Huawei Hisilicon Studentship at the
Department of Computer Science and Technology of the University of Cambridge.
X.D. acknowledges support from the Oxford-Man Institute of Quantitative
Finance and the EPSRC (EP/T023333/1).

## References

  * Borgwardt et al. [2020] Karsten Borgwardt, Elisabetta Ghisu, Felipe Llinares-López, Leslie O’Bray, Bastian Rieck, et al.  Graph kernels: State-of-the-art and future challenges.  _Foundations and Trends® in Machine Learning_ , 13(5-6), 2020.
  * Borgwardt and Kriegel [2005] Karsten M. Borgwardt and Hans-Peter Kriegel.  Shortest-path kernels on graphs.  In _Fifth IEEE International Conference on Data Mining_ , 2005.
  * Borgwardt et al. [2005] Karsten M. Borgwardt, Cheng Soon Ong, Stefan Schönauer, S. V. N. Vishwanathan, Alex J. Smola, and Hans-Peter Kriegel.  Protein function prediction via graph kernels.  _Bioinformatics_ , 21, 2005.
  * Bruna et al. [2014] Joan Bruna, Wojciech Zaremba, Arthur Szlam, and Yann LeCun.  Spectral networks and locally connected networks on graphs.  In _2nd International Conference on Learning Representations_ , 2014\.
  * Byrd et al. [1995] Richard H. Byrd, Peihuang Lu, Jorge Nocedal, and Ciyou Zhu.  A limited memory algorithm for bound constrained optimization.  _SIAM Journal on Scientific Computing_ , 16(5), 1995.
  * Chung [1997] Fan R. K. Chung.  _Spectral Graph Theory_.  American Mathematical Society, 1997.
  * Dai et al. [2016] Hanjun Dai, Bo Dai, and Le Song.  Discriminative embeddings of latent variable models for structured data.  In _Proceedings of the 33rd International Conference on Machine Learning_ , 2016.
  * Debnath et al. [1991] Asim Kumar Debnath, Rosa L. Lopez de Compadre, Gargi Debnath, Alan J. Shusterman, and Corwin Hansch.  Structure-activity relationship of mutagenic aromatic and heteroaromatic nitro compounds. correlation with molecular orbital energies and hydrophobicity.  _Journal of Medicinal Chemistry_ , 34(2), 1991.
  * Defferrard et al. [2016] Michaël Defferrard, Xavier Bresson, and Pierre Vandergheynst.  Convolutional neural networks on graphs with fast localized spectral filtering.  In _Advances in Neural Information Processing Systems_ , 2016.
  * Dong et al. [2013] Xiaowen Dong, Antonio Ortega, Pascal Frossard, and Pierre Vandergheynst.  Inference of mobility patterns via spectral graph wavelets.  In _2013 IEEE International Conference on Acoustics, Speech and Signal Processing_ , 2013.
  * Dong et al. [2020] Xiaowen Dong, Dorina Thanou, Laura Toni, Michael M. Bronstein, and Pascal Frossard.  Graph signal processing for machine learning: A review and new perspectives.  _IEEE Signal Processing Magazine_ , 37(6), 2020.
  * Duvenaud et al. [2015] David K Duvenaud, Dougal Maclaurin, Jorge Iparraguirre, Rafael Bombarell, Timothy Hirzel, Alán Aspuru-Guzik, and Ryan P Adams.  Convolutional networks on graphs for learning molecular fingerprints.  In _Advances in neural information processing systems_ , 2015.
  * Erdös and Rényi [1959] Paul Erdös and Alfréd Rényi.  On random graphs I.  _Publicationes Mathematicae Debrecen_ , 1959.
  * Errica et al. [2020] Federico Errica, Marco Podda, Davide Bacciu, and Alessio Micheli.  A fair comparison of graph neural networks for graph classification.  In _8th International Conference on Learning Representations_ , 2020\.
  * Gilmer et al. [2017] 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_ , 2017.
  * Hamilton et al. [2017] Will Hamilton, Zhitao Ying, and Jure Leskovec.  Inductive representation learning on large graphs.  In _Advances in neural information processing systems_ , 2017.
  * Hammond et al. [2011] David K Hammond, Pierre Vandergheynst, and Rémi Gribonval.  Wavelets on graphs via spectral graph theory.  _Applied and Computational Harmonic Analysis_ , 30(2), 2011\.
  * Hensman et al. [2013] James Hensman, Nicolò Fusi, and Neil D. Lawrence.  Gaussian processes for big data.  In _Proceedings of the 29th Conference on Uncertainty in Artificial Intelligence_ , 2013.
  * Hensman et al. [2015] James Hensman, Alexander G. de G. Matthews, and Zoubin Ghahramani.  Scalable variational gaussian process classification.  In _Proceedings of the 18th International Conference on Artificial Intelligence and Statistics_ , 2015.
  * Holland et al. [1983] Paul W. Holland, Kathryn Blackmond Laskey, and Samuel Leinhardt.  Stochastic blockmodels: First steps.  _Social Networks_ , 5(2), 1983.
  * Huang and Villar [2021] Ningyuan Teresa Huang and Soledad Villar.  A short tutorial on the weisfeiler-lehman test and its variants.  In _IEEE International Conference on Acoustics, Speech and Signal Processing_ , 2021.
  * Kipf and Welling [2017] Thomas N. Kipf and Max Welling.  Semi-supervised classification with graph convolutional networks.  In _5th International Conference on Learning Representations_ , 2017\.
  * Kondor and Pan [2016] Risi Kondor and Horace Pan.  The multiscale laplacian graph kernel.  In _Advances in neural information processing systems_ , 2016.
  * Kriege et al. [2020] Nils M Kriege, Fredrik D Johansson, and Christopher Morris.  A survey on graph kernels.  _Applied Network Science_ , 5(1), 2020.
  * Le Magoarou et al. [2018] Luc Le Magoarou, Rémi Gribonval, and Nicolas Tremblay.  Approximate fast graph fourier transforms via multilayer sparse approximations.  _IEEE Transactions on Signal and Information Processing over Networks_ , 4(2), 2018.
  * Li et al. [2012] Bin Li, Xingquan Zhu, Lianhua Chi, and Chengqi Zhang.  Nested subtree hash kernels for large-scale graph classification over streams.  In _2012 IEEE 12th International Conference on Data Mining_ , 2012\.
  * Liu and Nocedal [1989] Dong C. Liu and Jorge Nocedal.  On the limited memory BFGS method for large scale optimization.  _Mathematical Programming_ , 45(3, (Ser. B)), 1989.
  * Murphy et al. [2019] Ryan Murphy, Balasubramaniam Srinivasan, Vinayak Rao, and Bruno Ribeiro.  Relational pooling for graph representations.  In _Proceedings of the 36th International Conference on Machine Learning_ , 2019.
  * Nikolentzos et al. [2021] Giannis Nikolentzos, Giannis Siglidis, and Michalis Vazirgiannis.  Graph kernels: A survey.  _Journal of Artificial Intelligence Research_ , 72, 2021.
  * Opolka et al. [2022] Felix Opolka, Yin-Cong Zhi, Pietro Liò, and Xiaowen Dong.  Adaptive gaussian processes on graphs via spectral graph wavelets.  In _Proceedings of The 25th International Conference on Artificial Intelligence and Statistics_ , 2022.
  * Ortega [2022] Antonio Ortega.  _Introduction to Graph Signal Processing_.  Cambridge University Press, 2022.
  * Ortega et al. [2018] Antonio Ortega, Pascal Frossard, Jelena Kovačević, José MF Moura, and Pierre Vandergheynst.  Graph signal processing: Overview, challenges, and applications.  _Proceedings of the IEEE_ , 106(5), 2018.
  * Pineau [2019] Edouard Pineau.  Using laplacian spectrum as graph feature representation.  _arXiv preprint arXiv:1912.00735_ , 2019.
  * Porter and Canagarajah [1996] Robert Porter and Nishan Canagarajah.  A robust automatic clustering scheme for image segmentation using wavelets.  _IEEE transactions on image processing_ , 5(4), 1996.
  * Rasmussen and Williams [2005] Carl Edward Rasmussen and Christopher K. I. Williams.  _Gaussian Processes for Machine Learning (Adaptive Computation and Machine Learning)_.  The MIT Press, 2005.
  * Shervashidze et al. [2011] Nino Shervashidze, Pascal Schweitzer, Erik Jan van Leeuwen, Kurt Mehlhorn, and Karsten M. Borgwardt.  Weisfeiler-lehman graph kernels.  _Journal of Machine Learning Research_ , 12, 2011.
  * Shuman et al. [2013] David I Shuman, Sunil K Narang, Pascal Frossard, Antonio Ortega, and Pierre Vandergheynst.  The emerging field of signal processing on graphs: Extending high-dimensional data analysis to networks and other irregular domains.  _IEEE signal processing magazine_ , 30(3), 2013.
  * Titsias [2009] Michalis Titsias.  Variational learning of inducing variables in sparse gaussian processes.  In _Proceedings of the 12th International Conference on Artificial Intelligence and Statistics_ , 2009.
  * Togninalli et al. [2019] Matteo Togninalli, Elisabetta Ghisu, Felipe Llinares-López, Bastian Rieck, and Karsten Borgwardt.  Wasserstein weisfeiler-lehman graph kernels.  _Advances in Neural Information Processing Systems_ , 32, 2019.
  * Wale and Karypis [2006] Nikil Wale and George Karypis.  Comparison of descriptor spaces for chemical compound retrieval and classification.  In _Sixth International Conference on Data Mining_ , 2006.
  * Wu et al. [2019] Felix Wu, Amauri Souza, Tianyi Zhang, Christopher Fifty, Tao Yu, and Kilian Weinberger.  Simplifying graph convolutional networks.  In _Proceedings of the 36th International Conference on Machine Learning_ , 2019.
  * Xu et al. [2019] Keyulu Xu, Weihua Hu, Jure Leskovec, and Stefanie Jegelka.  How powerful are graph neural networks?  In _7th International Conference on Learning Representations_ , 2019\.
  * Yanardag and Vishwanathan [2015] Pinar Yanardag and S.V.N. Vishwanathan.  Deep graph kernels.  In _Proceedings of the 21th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining_ , 2015.
  * Ying et al. [2018] Rex Ying, Jiaxuan You, Christopher Morris, Xiang Ren, William L. Hamilton, and Jure Leskovec.  Hierarchical graph representation learning with differentiable pooling.  In _Advances in Neural Information Processing Systems_ , 2018.
  * Zhang et al. [2018] Muhan Zhang, Zhicheng Cui, Marion Neumann, and Yixin Chen.  An end-to-end deep learning architecture for graph classification.  _Proceedings of the AAAI Conference on Artificial Intelligence_ , 32(1), 2018.