# Silicon spin chains at finite temperature: dynamics of Si(553)-Au

Steven C. Erwin<EMAIL_ADDRESS>Center for Computational Materials
Science, Naval Research Laboratory, Washington, DC 20375  P.C. Snijders
Materials Science and Technology Division, Oak Ridge National Laboratory, Oak
Ridge, TN 37830

###### Abstract

When gold is deposited on Si(553), the surface self-assembles to form a
periodic array of steps with nearly perfect structural order. In scanning
tunneling microscopy these steps resemble quasi-one-dimensional atomic chains.
At temperatures below $\sim$50 K the chains develop a tripled periodicity. We
recently predicted, on the basis of density-functional theory calculations at
$T=0$, that this tripled periodicity arises from the complete polarization of
the electron spin on every third silicon atom along the step; in the ground
state these linear chains of silicon spins are antiferromagnetically ordered.
Here we explore, using ab-initio molecular dynamics and kinetic Monte Carlo
simulations, the behavior of silicon spin chains on Si(553)-Au at finite
temperature. Thermodynamic phase transitions at $T>0$ in one-dimensional
systems are prohibited by the Mermin-Wagner theorem. Nevertheless we find that
a surprisingly sharp onset occurs upon cooling—at about 30 K for perfect
surfaces and at higher temperature for surfaces with defects—to a well-ordered
phase with tripled periodicity, in good agreement with experiment.

##  I Introduction

Linear atomic chains of metal atoms on semiconductor surfaces offer, in
principle, the physical realization of phenomena predicted theoretically for
one-dimensional model systems. In practice, however, unanticipated
interactions can often complicate the picture and lead to behavior not easily
explained by simple models. In this article we demonstrate theoretically how
the complex interactions among polarized electron spins in silicon surface
states determine the observed behavior of a well-studied atomic chain system,
Si(553)-Au, over a wide range of temperatures. The methods developed here and
the resulting predictions—which are qualitatively and quantitatively
consistent with experimental observations—are also likely to apply more
broadly to other vicinal Si/Au chain systems, such as Si(557)-Au.

The Si(553)-Au surface was first investigated in Ref. Crain  _et al._ , 2003,
an experimental study which established that the electronic band dispersion
and fermi surface were indeed those of a nearly one-dimensional metal. Since
then, numerous lines of research have emerged. Efforts to determine the basic
atomic structure of the surface have been based on data from diffraction
experimentsGhose  _et al._ (2005); Takayama  _et al._ (2009); Voegeli  _et
al._ (2010) and on the results of theoretical total-energy
calculations.Riikonen and Sanchez-Portal (2005, 2006, 2008); Krawiec (2010)
These were greatly aided by the first definitive determination of the coverage
of Au atoms on Si(553)-Au.Barke  _et al._ (2009) Other investigations have
explored the properties of finite-length chainsCrain and Pierce (2005); Crain
_et al._ (2006) as well as various native defectsOkino  _et al._ (2007a) and
foreign adsorbates Ryang  _et al._ (2007); Okino  _et al._ (2007b); Ahn  _et
al._ (2008); Kang  _et al._ (2009a, a, b); Nita  _et al._ (2011); Krawiec and
Jalochowski (2013) on the nominally clean Si(553)-Au surface.

One particularly interesting line of research has focused on the collective
behavior in Si(553)-Au that emerges at low temperature. Ideal one-dimensional
metals with partially filled bands exhibit a broken symmetry at low
temperature, namely a charge-density wave arising from the Peierls
instability. Indeed, broken symmetries in Si(553)-Au were observed using
scanning tunneling microscopy (STM) in Refs. Ahn  _et al._ , 2005 and Snijders
_et al._ , 2006. Images acquired at room temperature showed alternating bright
and dim rows with unit periodicity $a_{0}$ along the rows. Below $\sim$50 K
these rows separately developed higher-order periodicity: a tripled period
(3$a_{0}$) along the bright rows and a doubled period (2$a_{0}$) along the dim
rows. Subsequent review articles have discussed possible explanations for
these observations.Snijders and Weitering (2010); Hasegawa (2010)

Notwithstanding the fact that Peierls instabilities lead to higher-order
periodicity, a completely different theoretical explanation for the coexisting
triple and double periodicities in Si(553)-Au was proposed in Ref. Erwin and
Himpsel, 2010. The key idea, which was based on the results of density-
functional theory (DFT) calculations, was that the ground state of Si(553)-Au
is spin polarized. In the DFT ground state, the silicon atoms that comprise
the steps on this vicinal surface have dangling bonds, every third of which is
occupied by a single fully polarized electron while the other two are doubly
occupied. The bright rows seen in empty-state STM images arise from these
step-edge silicon atoms. At low temperature the $3a_{0}$ peaks that appear in
this row are from the spin-polarized atoms, which have precisely this
periodicity. The DFT ground state also reveals a period doubling within the
row of Au atoms. Both of these higher-order periodicities disappear if spin
polarization is suppressed in the calculation, providing compelling evidence
that spin polarization is the primary mechanism underlying the observed
symmetry breaking in Si(553)-Au.

Experiments were subsequently carried out to look for a spectroscopic
signature of this predicted spin-polarized ground state. The DFT calculations
showed that an unoccupied state should exist 0.5 eV above the fermi level and
be localized at the polarized silicon atoms.Erwin and Himpsel (2010) The
existence and spectral and spatial location of this state was indeed confirmed
by two-photon photoemission Biedermann  _et al._ (2012) and by scanning
tunneling spectroscopy.Snijders  _et al._ (2012)

The predictions of Ref. Erwin and Himpsel, 2010 only addressed the zero-
temperature ground state of Si(553)-Au. Left unanswered in that work was the
question of how the broken-symmetry ground state evolves to have normal
$a_{0}$ periodicity above $\sim$50 K. This article addresses that question
from a theoretical and computational perspective. Although it may seem obvious
that thermal fluctuations are important, the nature of these fluctuations
turns out to be unexpectedly subtle. Nevertheless, we derive here a number of
detailed qualitative as well as quantitative predictions that can easily be
tested experimentally. The results of these tests will furnish additional
evidence for evaluating the validity of the basic mechanism proposed in Ref.
Erwin and Himpsel, 2010.

##  II Ground state configuration

The physical and magnetic structure of Si(553)-Au in its ground state were
first proposed and discussed in Ref. Erwin and Himpsel, 2010 and for reference
are reproduced in Fig. 1. This is a stepped surface consisting of (111)
terraces and bilayer steps, and is stabilized by Au atoms that substitute for
Si atoms in the surface layer of the terrace. The steps themselves consist of
Si atoms organized into a thin graphitic strip of honeycomb hexagons (the
green atoms in Fig. 1).

Figure 1: (Color online) (a,b) Perspective and top views of Si(553)-Au in its
electronic ground state. Yellow atoms are Au, all others are Si. The Au atoms
are embedded in flat terraces, which are separated by steps consisting of Si
atoms arranged as a honeycomb chain (green). Every third Si atom (red, blue)
at the step has a spin magnetic moment of one Bohr magneton ($S=1/2$, arrows)
from the complete polarization of the electron occupying the dangling-bond
orbital. The sign of the polarization (red vs. blue) alternates along the
step. The six atoms in the outlined box are the focus of the ab-initio
molecular dynamics discussed in Sec. III.

The surface electronic structure of Si(553)-Au has two main contributions. The
first consists of two intense quasi-1D parabolic electron bands centered at
the boundary of the surface Brillouin zone. These “Au bands” arise from the
bonding and antibonding combinations of Au 6$s$ and subsurface Si orbitals
(purple atoms). The bonding Au band is approximately half-filled and the
antibonding band approximately one-fourth filled.

The second contribution arises from the very edge of the Si honeycomb chain,
which consists of threefold-coordinated Si atoms. The unpassivated $sp^{3}$
orbitals of these atoms can in principle be occupied by zero, one, or two
electrons. The Si atoms themselves supply, on average, one electron per
orbital. The step edge does not necessarily maintain this average occupancy,
because electronic charge can also be transferred to or from the Au bands.
Indeed, DFT calculations predict that the lowest energy configuration has one
electron in every third orbital (red and blue atoms in Fig. 1) and double
occupancy everywhere else (green atoms). The singly occupied orbitals are
completely spin-polarized and hence have local spin moments of 1 bohr magneton
each. Physically, these atoms relax slightly downward, by 0.3 Å, compared to
their nonpolarized neighbors. The sign of the spins alternates along the step
edge, with antiferromagnetic order favored by 15 meV per spin relative to
ferromagnetic order. Therefore the magnetic periodicity is 6$a_{0}$, where
$a_{0}$ is the Si surface lattice constant. This is also the smallest period
that allows for the coexistence of 3$a_{0}$ spacing of the spins and 2$a_{0}$
spacing (period doubling) within the Au chain. This coexistence was first
observed in STM experimentsAhn  _et al._ (2005); Snijders  _et al._ (2006) and
emerges naturally in DFT calculations—but only when the spin degree of freedom
is unconstrained.Erwin and Himpsel (2010)

##  III Finite temperature dynamics

The remainder of this article explores excitations of Si(553)-Au from its
ground state due to finite temperature. Two main theoretical tools were used:
ab-initio molecular dynamics (MD) and kinetic Monte Carlo (kMC) simulations.
The first was used to identify the most important low-energy activated
processes and to determine their activation barriers. Because of the
complexity of the system only small time scales (tens of ps) and a small
(1$\times$6) simulation cell could be addressed using ab-initio MD. To reach
much longer time scales (tens of ns) and larger system sizes (up to 128 spins)
we constructed a one-dimensional kMC model based on the processes and barriers
determined from ab-initio MD. In particular, the kMC model allowed us to
investigate finite-temperature behavior in the presence of pinning
defects—providing useful insight into temperature-dependent results from
scanning probe experiments, where defects often play a critical role.

Two simplifying assumptions were used throughout this work. (1) Electronic
excitations were not considered, and consequently the system stays on the
Born-Oppenheimer surface. This assumption is reasonable in view of the modest
temperatures—room temperature and lower—considered here. (2) Spin flips were
not allowed. Although, as we will see below, the spins can diffuse among the
Si step-edge atoms, their signs and ordering remained that of the original
antiferromagnetic ordering. Although a different initial spin ordering might
affect some details of the simulation, the overall qualitative findings would
be very similar.

###  III.1 Ab-initio molecular dynamics

The MD simulations were performed using the same basic geometry and
computational parameters described in Ref. Erwin and Himpsel, 2010. The
Si(553)-Au surface was represented by six layers of Si plus the reconstructed
top surface layer and a vacuum region of 10 Å. All atoms were free to move
during the simulation except the bottom Si layer and its passivating hydrogen
layer. Total energies and forces were calculated within the generalized-
gradient approximation of Perdew, Burke, and Ernzerhof to DFT using projector-
augmented wave potentials, as implemented in VASP.Kresse and Hafner (1993);
Kresse and Furthmüller (1996); Blochl (1994); Kresse and Joubert (1999) The
plane-wave cutoff was 200 eV and only the $\Gamma$ point was used. The
dynamics simulations were performed in the canonical ensemble using a Nosé
thermostat and a time step of 3 fs. Five temperatures, equally spaced in
$1/T$, were used (57, 67, 80, 110, 133 K). For each temperature a
thermalization run of 10 ps was first performed, followed by a dynamics run of
20 ps.

Figure 2 shows the resulting atomic trajectories during the entire run of 104
MD time steps for the lowest temperature studied, 57 K. The six curves are for
the six Si step-edge atoms in the outlined box of Fig. 1(b). The upper and
lower panels show the relative heights of the atoms and their local spin
moments, respectively. After thermalization was achieved the system settled
into its ground state configuration with two spin-polarized atoms (red and
blue) sitting $\sim$0.3 Å lower than their four non-polarized neighbors.

Figure 2: (Color online) (a) Ab-initio molecular dynamics trajectories of the
six Si step-edge atoms outlined in Fig. 1, at 57 K. Red and blue curves denote
the red and blue atoms, which are initially spin-polarized. Other colors
(magenta, cyan, dark green, light green) denote initially non-polarized atoms.
Thermalization is completed by about 10 ps. Upper panel: height of each atom,
relative to the average height of nonpolarized atoms. Lower panel: local spin
moment of each atom. (b) Expanded view of two spin hops occurring at 15.44 ps
(from the red atom to the cyan atom) and at 15.71 ps (from the blue atom to
the magenta atom).

The expanded view in Fig. 2(b) focuses on two events that occurred between 15
and 16 ps. At 15.44 ps the magnitude of the moment on the spin-up red atom
went rapidly to zero while, concurrently, a spin-up moment rapidly developed
on the neighboring cyan atom. At essentially the same time the height of the
red atom increased by 0.3 Å to that of a non-polarized atom, while the cyan
atom moved down by the same amount. In summary, the spin-up moment that was
localized on the red atom hopped to one of its neighbors.

Very soon after, a second hop occurred at 15.71 ps. This hop was made by the
other spin (with the opposite sign) which moved from the blue atom to the
magenta atom. It is not a coincidence that this hop occurred so soon after the
first. The first hop changed the minimum spacing between spins from 3$a_{0}$
to 2$a_{0}$, incurring an energy penalty (discussed in detail below). This
increase in energy in turn reduced the barrier for any hop that restores the
spacing to its optimal value. The cyan and magenta atoms are indeed separated
by 3$a_{0}$, and thus after two rapid spin hops the system was restored to an
equivalent ground state configuration, in which it remained for the rest of
the simulation.

The very small number of hops observed at 57 K makes it clear that ab-initio
MD simulations of Si(553)-Au at still lower temperatures, where many of the
relevant experiments are conducted, are not feasible. Instead we turn to
higher temperatures and ask how the frequency of hopping events depends on
temperature. This information will be useful in Sec. III.2 for calibrating and
validating our kMC model in a temperature range accessible to both methods.

Figure 3 shows the resulting time-averaged hopping rate for a single spin,
versus inverse temperature. The rates for low temperatures have large
statistical uncertainties (not shown) and hence it is reasonable to describe
these results by a simple linear Arrhenius fit, as shown. The attempt
frequency, $2.0\times 10^{13}$ s-1, is on the order of a surface vibrational
frequency, as expected. The activation energy, 12 meV, represents a
characteristic average of the individual barriers for spin hops weighted by
their relative probability of occurrence.

Figure 3:  Temperature dependence of the hopping rate for Si spins along the
Si(553)-Au step edge. The rates are time averages extracted from ab-initio
molecular dynamics (MD) trajectories, and are compared to rates from kinetic
Monte Carlo (kMC) simulations. The linear fits describe Arrhenius behavior.
The fit to ab-initio MD rates (thick line) gives a pre-exponential factor
$2.0\times 10^{13}$ s-1 and activation barrier 12 meV. For kMC rates (thin
line) the values are $2.2\times 10^{13}$ s-1 and 13 meV.

###  III.2 Kinetic Monte Carlo model

To construct the kMC model one first needs to enumerate all the relevant spin
hops and their individual rates. We used DFT results obtained from the full
Si(553)-Au system for this task. In the spirit of simplicity we constructed
the kMC model itself to be strictly one-dimensional, with an arbitrarily large
unit cell and periodic boundary conditions. Thus the kMC simulations inherit
much of the accuracy of the DFT calculations but make the additional
approximation that spin hops along different step edges are independent.

Figure 4 shows the DFT potential energy surface for the spin hop observed in
Fig. 2(b) at 15.44 ps into the MD simulation. Because the spins and the
heights of the atoms are tightly linked, the reaction coordinate $x$ is
approximately given by the relative heights $h$ of the red and cyan atoms,

$x\approx[1-(h_{\rm cyan}-h_{\rm red})/\Delta h]/2,$ (1)

where $\Delta h=0.3$ Å is the equilibrium height difference between spin-
polarized and non-polarized atoms. To definitively determine the detailed
reaction pathway and potential energy surface we used the nudged elastic-band
method.

Figure 4:  DFT potential energy surface for a single spin hopping from the red
atom to the neighboring cyan atom. The initial state (0) is the ground state
depicted in Fig. 1. The final state (1) is the metastable state, fully
relaxed, that exists between 15.44 and 15.71 ps in the MD simulation of Fig.
2. The activation barrier is 30 meV for the forward reaction and 5 meV for the
reverse reaction. Atom colors correspond to the trajectories in Fig. 2.

This potential energy surface confirms the assertion, made in Sec. III.1, that
the red-to-cyan (forward) spin hop incurs an energy penalty that leads to a
smaller barrier for the cyan-to-red (reverse) hop. Specifically, the
activation barrier for the forward hop is 30 meV, the resulting energy penalty
is 25 meV, and the barrier for the reverse hop is 5 meV. These two types of
hops, and their calculated barriers, are two of the three fundamental
processes included in our kMC model.

For convenience we define here a more compact notation for enumerating the
different types of spin hops. Careful examination of the ab-initio MD
trajectories reveals that all spin hops were to an adjacent site; there were
no double hops. Hence we can label every hop as either leftward ($\leftarrow$)
or rightward ($\rightarrow$). We assume that the barrier for a spin hop
depends only on the spin’s immediate environment, that is, on the distances
$ma_{0}$ and $na_{0}$ to the left and right neighboring spins, respectively,
measured before making the hop. Using this notation we can express the
barriers for the two hops shown in Fig. 4 as $b(3,3,\leftarrow)=30$ meV and
$b(2,4,\rightarrow)=5$ meV, where the two numerical arguments denote $m$ and
$n$, respectively.

The third important hop we considered occurs when $m+n=5$, rather than 6 as
depicted in Fig. 4. From DFT nudged elastic-band calculations we find
$b(3,2,\leftarrow)$ = $b(2,3,\rightarrow)=14$ meV (the barriers are equal by
symmetry). As expected from the distances to the neighboring spins, this
barrier is in between the previous two.

The MD trajectories also show that two spins never occupy adjacent sites.
Because of this, our enumeration of the possible spin hops is already complete
for all cases with $m+n\leq 6$. (It is worth noting that the configuration in
which a spin has both neighbors at $2a_{0}$ is allowed, but because its
adjacent sites cannot be occupied this spin cannot hop until one of its
neighbors does.) The cases with $m+n\geq 7$ are difficult to treat within DFT
but occur more rarely and thus are less important. For this reason we treated
the effect of neighbors beyond $3a_{0}$ as negligible, used the barrier of 30
meV for any hop that brings a spin within $2a_{0}$ of its neighbor, and
assigned a single (arbitrary) barrier of 10 meV to hops that maintain larger
separations than this. This completes our enumeration.

To finish the construction of the kMC model, we assumed that the rates for all
allowed spin hops are given by $r=a\exp(-b/kT)$, where $a$ is a common
prefactor and $b=b(m,n,\leftarrow)$ and $b(m,n,\rightarrow)$ are the DFT
barriers.

To determine the optimal value of $a$ and compare the predictions of the kMC
model to the ab-initio MD results, we applied the model to the system
discussed in Sec. III.1—two spin-polarized atoms in a six-atom unit cell with
periodic boundary conditions. The resulting kMC spin hopping rates obtained
using $a=6\times 10^{12}$ s-1 are plotted in Fig. 3 for direct comparison with
the rates from ab-initio MD. The kMC rates have negligible statistical errors
and it is clear that a simple Arrhenius fit describes them very well.
Moreover, the fitted attempt frequency, $2.2\times 10^{13}$ s-1, and
activation energy, 12 meV, are within 10% of the MD values. This confirms that
the kMC model accurately reproduces the ab-initio results within the
temperature range considered.

###  III.3 Spins at finite temperature near a defect

In real systems, the behavior of collective phenomena is often controlled by
defects that pin the phase of a low-temperature state having broken symmetry.
On the Si(553)-Au surface, a variety of defects—missing atoms, absorbates,
etc.—have been observed to act as pinning sites that locally stabilize the
1$\times$3 ground state.Snijders  _et al._ (2006); Kang  _et al._ (2009b);
Hasegawa (2010); Shin  _et al._ (2012) The important role played by such
pinning defects motivates our first application of the kMC model.

We prepared a system consisting of 128 independent spins, with periodic
boundary conditions, initially arranged in the antiferromagnetic ground state
with uniform 3$a_{0}$ spacing. One of the spins (at position 0) was pinned in
place throughout the simulation, thus representing a generic immobile defect.
The spins were allowed to hop stochastically among the 3$\times$128 lattice
sites according to probabilities defined by the hopping rates $r$.

Figure 5 displays the resulting trajectories at 300 K of all the spins over
the first 10 ns of the simulation. For clarity every tenth trajectory trace is
colored. As the system evolved, each spin explored a region of the lattice
around its initial position. These explorations were relatively small for
spins near the pinning defect and became progressively larger for spins
farther away.

Figure 5: (Color online) Kinetic Monte Carlo trajectories of 128 spins at 300
K in the presence of a pinning defect at the origin. Every tenth trace is
colored for clarity. Right panel: histogram of the positions occupied by each
spin, weighted by the time spent there, obtained over a simulation time of 1
$\mu$s. The heavy curve is the envelope function $d^{-2/3}$ describing the
decay of histogram heights with distance $d$ from the pinning defect.

The right panel in Fig. 5 examines this thermally induced wandering in greater
detail. For each of the 128 spins a histogram was made representing the
position of that spin at 300 K. At this temperature a simulation time of 1
$\mu$s was sufficient to obtain the steady-state distribution. It is readily
apparent from examining the colored histograms that each is well described by
a gaussian function centered on the spin’s initial position. Thus each of
these gaussians is entirely specified by its variance $\sigma^{2}$, whose
value depends on the distance $d$ to the pinning defect. To deduce this
dependence we first note that the area under each gaussian is by construction
the same. Hence the height of each gaussian is proportional to $1/\sigma$. We
find empirically that the dependence of these heights on distance is given
with excellent accuracy as $d^{-2/3}$. An envelope function with this
dependence is shown on the histogram plot as a heavy black curve. From this
dependence we thus deduce that the thermally induced widths $w$, defined here
as $2\sigma$, increase with distance from a pinning defect as $w\sim d^{2/3}$.

Now we move on to explore how temperature affects the thermal wandering of
spins near a pinning defect. We repeated the kMC simulation and analysis in
Fig. 5 for a series of temperatures between 10 and 300 K. We focus on the
variation of the thermal widths $w$ as a function of temperature $T$.

Figure 6: (Color online) Temperature dependence of the thermal widths $w(T)$
of every tenth spin, in the presence of a pinning defect. Labels indicate the
spin’s distance from the defect, in units of 3$a_{0}$. Colors correspond to
the trajectories in Fig. 5. The light gray shaded area is bounded by
logarithmic fits to the lower (red) and upper (orange) data points. The
characteristic temperatures $T_{0}$, $T_{1}$, and $T_{3}$ describe different
criteria by which thermal wandering of spins is expected to be either
eliminated or suppressed; see discussion in text.

Figure 6 summarizes the resulting temperature dependence. The six datasets
show $w(T)$ for every tenth spin of the 128-spin simulation cell. For
reference, the six values at $T=300$ K correspond to the six gaussian widths
in the upper half of the histogram panel in Fig. 5. We find empirically that
the dependence of each dataset on temperature is close to logarithmic, as
shown by the light gray shaded area. This implies that the dependence of the
thermal widths on distance and temperature can be separated and written as

$w(d,T)=w_{0}(d)\ln(T/T_{0}),$ (2)

where $w_{0}(d)\sim d^{2/3}$ and the characteristic temperature $T_{0}$ has
the fitted value 21 K. All thermal wandering is, by definition, completely
eliminated at $T_{0}$. But two less restrictive criteria may be more relevant
for interpreting the experimentally observed transition to the period-tripled
ground state. At $T_{1}=27$ K the thermal widths for all 128 spins become
smaller than the width of a single lattice site. Hence, below this temperature
the spins will in effect be frozen into place on every third lattice site. At
still higher temperature, $T_{3}=42$ K, all thermal widths are less than or
equal to the average spacing (three lattice sites) between spins. Hence the
spins will first become distinguishable as the system is cooled below this
temperature.

To generalize this result and make predictions that can be tested by
experiment, we first assume that real systems can be characterized by a known
average concentration $c$ of pinning defects. Because the defects are
distributed in 1D, the characteristic distance $d$ from a spin to the nearest
defect scales as $1/c$. By inverting Eq. 2 we then immediately obtain a simple
result: the temperatures $T_{1,3}$ at which all spins in the system become
either frozen into place or distinguishable will scale as
$T_{1,3}\sim\exp(c^{2/3})$. Thus we predict this scaling to describe the
temperature at which the period-tripled ground state of Si(553)-Au is first
observed.

By inserting into this qualitative relationship the appropriate constants
obtained from the kMC simulations, we derive a quantitative prediction for the
freezing temperature,

$T_{1}=T_{0}\exp(k\,c^{2/3}),$ (3)

where $k=11.8$ is a dimensionless constant and $c$ is expressed in the
dimensionless units of defects per lattice site. The corresponding equation
for $T_{3}$ can be obtained from Eq. 3 by multiplying the argument of the
exponential by three. A useful guide to understanding the importance of
defects in Si(553)-Au is provided by linearizing Eq. 3 around a physically
plausible value (10-2) for the defect concentration. This leads to the result
that a change in the defect concentration will raise the freezing temperature
by $\Delta T_{1}=\gamma\Delta c$, with proportionality constant $\gamma=1330$
K. Thus, for samples with approximately one defect every 100 lattice sites, a
doubling of this concentration will increase the freezing temperature by 13 K.

###  III.4 Spins at finite temperature in the absence of defects

Although a system completely free of defects is obviously unrealistic, the
behavior of such an idealized system nevertheless offers complementary insight
into the thermal wandering of spins when the concentration of defects is very
low. As we show below, despite the absence of defects, the statistical
behavior of the spins still exhibits a sudden and qualitative change at about
30 K.

Figure 7: Trajectory of a single spin at 100 K in the absence of pinning
defects. Gray curve shows the theoretical average displacement versus time for
a isolated random walker in one dimension, $\langle d\rangle\propto\sqrt{t}$.

We constructed a periodic system of 64 spins similar to that described in Sec.
III.3, but now without a pinning defect. Thus each spin executed a random walk
in 1D. Figure 7 shows a typical trajectory at 100 K for one of the 64 spins.
Despite the stochastic nature of this single trajectory, it is already
plausible that the average displacement $\langle d\rangle$ depends on the time
$t$ according to $\langle d\rangle\propto\sqrt{t}$, which is the well-known
result for a single unbiased random walker in one dimension.

Figure 8:  Average displacement versus time of a single spin in the absence of
defects, at the indicated temperatures. Circles are statistical averages over
2000 kMC simulations. Straight lines are fits to $t^{H}$, where $H$ is the
Hurst exponent. The dotted line indicates $H$=1/2.

To analyze this behavior more systematically, we performed many independent
kMC simulations and computed the average displacements as a function of time.
Figure 8 shows these averages on a log-log scale for a range of temperatures.
At high temperatures we indeed obtain the behavior $\langle d\rangle\propto
t^{1/2}$, which is indicated by the dotted line. This behavior persists until
the temperature reaches the range 30–35 K, where it still exhibits power-law
behavior $\langle d\rangle\propto t^{H}$ but with a progressively larger
exponent $H>1/2$.

Figure 9:  Hurst exponent versus temperature for spins in the absence of a
defect. The solid curve is a fit to a generalized susceptibility with a
characteristic temperature $T_{c}=28$ K.

In general, a system for which the long-time displacements are characterized
by a Hurst exponent $H=1/2$ is said to be uncorrelated, while $H>1/2$
indicates that long-time correlations are present.Rangarajan and Ding (2000)
Figure 9 shows the Hurst exponent $H$, obtained by fitting the time-averaged
displacements, as a function of temperature. Above $\sim$40 K the system
displays uncorrelated behavior. Below this temperature we observe the rapid
onset of correlated behavior. To quantify the temperature of this transition
we fit the Hurst exponents to a generalized susceptibility of the form
$H=(1/2)/[1-(T_{c}/T)^{\nu}]$ and obtain a characteristic temperature
$T_{c}=28$ K. Hence as the clean system is cooled toward $T_{c}$ the random
walks described by individual spins rapidly lose their independent character.
This transition occurs with a characteristic temperature comparable to that
obtained, $T_{0}=21$ K, by extrapolating from the behavior in the presence of
defects. Thus these two complementary approaches lead to qualitatively as well
as quantitatively similar conclusions.

##  IV Discussion and conclusions

A guiding principle in one-dimensional physics is provided by the Mermin-
Wagner theorem, which states that phase transitions cannot occur above $T=0$
if the interactions are short-ranged.Mermin and Wagner (1966) Exploring the
different manifestations of this theorem in real systems can yield new and
unanticipated insights. For example, a previous publication by one of us
demonstrated theoretically that when the interactions are $not$ short-ranged,
as for $1/r$ Coulomb interactions, then a well-defined thermodynamic phase
transition can indeed occur—and likely does occur for a system of Ba
adsorbates on Si(111).Erwin and Hellberg (2005); Bruch (2005) In the present
article the new insights into 1D physics are of a different kind: we have
shown that even when the interactions are short-ranged, and the system purely
1D, a well-ordered phase with the clear signature of a broken symmetry can
form well above $T=0$. Moreover, an approximate transition temperature can be
readily identified using realistic simulations and straightforward statistical
analysis.

For Si(553)-Au the precise nature of the interacting entities is unexpected
and somewhat subtle: our MD simulations showed that they are neither simple
vibrations of atoms, nor spins on a simple fixed lattice, but rather a tightly
coupled combination of both. In this sense a description of Si(553)-Au based
on spin polarons is appropriate. Our kMC simulations showed that as the
temperature of the system is raised, the ground state $3a_{0}$ crystal formed
by these polarons melts at a temperature we estimate to be $\sim$30 K for
perfectly clean systems, and higher for systems with pinning defects. A direct
experimental test of this description is afforded by Eq. 3, which predicts how
the transition temperature varies with the concentration of defects.

It is important to acknowledge some limitations of this work. Because our
focus has been on the behavior of Si(553)-Au at low temperature, we have
assumed that the spin-polarized silicon states remain polarized at higher
temperatures. Our preliminary calculations indicate that this assumption is
completely justified at the temperatures of interest here. However, at much
higher temperatures the thermally induced vibrations of the step edge atoms
increasingly render the system nonpolarized for part of the time. For example,
at room temperature the average spin moment is reduced to roughly 2/3 of its
low-temperature value of 1 bohr magneton. Future theoretical investigations
into the behavior of Si(553)-Au near room temperature will have to account for
this thermal suppression of the spin polarization.

Finally, as mentioned in Sec. III, we have throughout assumed for simplicity
that the ordering of the spins remains antiferromagnetic, as in the ground
state. Preliminary calculations show that the barriers for spin hopping depend
quantitatively, although not qualitatively, on the signs of the neighboring
spins. A generalization of our kMC model that includes spin flips would be a
very interesting direction to pursue, but we anticipate that the qualitative
conclusions drawn here—as well as the overall consistency between our findings
and those of existing experiments—would be largely unchanged.

###### Acknowledgements.

Many discussions with F.J. Himpsel are gratefully acknowledged. This work was
supported by the Office of Naval Research (SCE) and the Department of Energy,
Basic Energy Sciences, Materials Sciences and Engineering Division (PCS).
Computations were performed at the DoD Major Shared Resource Centers at AFRL
and ERDC.

## References

  * Crain  _et al._ (2003) J. N. Crain, A. Kirakosian, K. N. Altmann, C. Bromberger, S. C. Erwin, J. McChesney, J.-L. Lin, and F. J. Himpsel, Phys. Rev. Lett. 90, 176805 (2003).
  * Ghose  _et al._ (2005) S. K. Ghose, I. K. Robinson, P. A. Bennett, and F. J. Himpsel, Surf. Sci. 581, 199 (2005).
  * Takayama  _et al._ (2009) T. Takayama, W. Voegeli, T. Shirasawa, K. Kubo, M. Abe, T. Takahashi, K. Akimoto, and H. Sugiyama, e-Journal of Surface Science and Nanotechnology 7, 533 (2009).
  * Voegeli  _et al._ (2010) W. Voegeli, T. Takayama, T. Shirasawa, M. Abe, K. Kubo, T. Takahashi, K. Akimoto, and H. Sugiyama, Physical Review B 82, 075426 (2010).
  * Riikonen and Sanchez-Portal (2005) S. Riikonen and D. Sanchez-Portal, Nanotechnology 16, S218 (2005).
  * Riikonen and Sanchez-Portal (2006) S. Riikonen and D. Sanchez-Portal, Surf. Sci. 600, 1201 (2006).
  * Riikonen and Sanchez-Portal (2008) S. Riikonen and D. Sanchez-Portal, Physical Review B 77, 165418 (2008).
  * Krawiec (2010) M. Krawiec, Physical Review B 81, 115436 (2010).
  * Barke  _et al._ (2009) I. Barke, F. Zheng, S. Bockenhauer, K. Sell, V. v. Oeynhausen, K. H. Meiwes-Broer, S. C. Erwin, and F. J. Himpsel, Physical Review B 79, 155301 (2009).
  * Crain and Pierce (2005) J. N. Crain and D. T. Pierce, Science 307, 703 (2005).
  * Crain  _et al._ (2006) J. N. Crain, M. D. Stiles, J. A. Stroscio, and D. T. Pierce, Phys. Rev. Lett. 96, 156801 (2006).
  * Okino  _et al._ (2007a) H. Okino, I. Matsuda, S. Yamazaki, R. Hobara, and S. Hasegawa, Physical Review B 76, 035424 (2007a).
  * Ryang  _et al._ (2007) K.-D. Ryang, P. G. Kang, H. W. Yeom, and S. Jeong, Physical Review B 76, 205325 (2007).
  * Okino  _et al._ (2007b) H. Okino, I. Matsuda, R. Hobara, S. Hasegawa, Y. Kim, and G. Lee, Physical Review B 76, 195418 (2007b).
  * Ahn  _et al._ (2008) J. R. Ahn, P. G. Kang, J. H. Byun, and H. W. Yeom, Physical Review B 77, 035401 (2008).
  * Kang  _et al._ (2009a) P.-G. Kang, H. Jeong, and H. W. Yeom, Physical Review B 79, 113403 (2009a).
  * Kang  _et al._ (2009b) P.-G. Kang, J. S. Shin, and H. W. Yeom, Surf. Sci. 603, 2588 (2009b).
  * Nita  _et al._ (2011) P. Nita, M. Jalochowski, M. Krawiec, and A. Stepniak, Phys. Rev. Lett. 107, 026101 (2011).
  * Krawiec and Jalochowski (2013) M. Krawiec and M. Jalochowski, Physical Review B 87, 5445 (2013).
  * Ahn  _et al._ (2005) J. R. Ahn, P. G. Kang, K. D. Ryang, and H. W. Yeom, Phys. Rev. Lett. 95, 196402 (2005).
  * Snijders  _et al._ (2006) P. C. Snijders, S. Rogge, and H. H. Weitering, Phys. Rev. Lett. 96, 076801 (2006).
  * Snijders and Weitering (2010) P. C. Snijders and H. H. Weitering, Reviews of Modern Physics 82, 307 (2010).
  * Hasegawa (2010) S. Hasegawa, J. Phys.: Condens. Matter 22, 084026 (2010).
  * Erwin and Himpsel (2010) S. C. Erwin and F. J. Himpsel, Nature Communications 1, 58 (2010).
  * Biedermann  _et al._ (2012) K. Biedermann, S. Regensburger, T. Fauster, F. J. Himpsel, and S. C. Erwin, Physical Review B 85, 245413 (2012).
  * Snijders  _et al._ (2012) P. C. Snijders, P. S. Johnson, N. P. Guisinger, S. C. Erwin, and F. J. Himpsel, New Journal of Physics 14, 103004 (2012).
  * Kresse and Hafner (1993) G. Kresse and J. Hafner, Phys. Rev. B 47, 558 (1993).
  * Kresse and Furthmüller (1996) G. Kresse and J. Furthmüller, Phys. Rev. B 54, 11169 (1996).
  * Blochl (1994) P. E. Blochl, Phys. Rev. B 50, 17953 (1994).
  * Kresse and Joubert (1999) G. Kresse and D. Joubert, Physical Review B 59, 1758 (1999).
  * Shin  _et al._ (2012) J. S. Shin, K.-D. Ryang, and H. W. Yeom, Physical Review B 85, 073401 (2012).
  * Rangarajan and Ding (2000) G. Rangarajan and M. Z. Ding, Physical Review E 61, 4991 (2000).
  * Mermin and Wagner (1966) N. D. Mermin and H. Wagner, Physical Review Letters 17, 1133 (1966).
  * Erwin and Hellberg (2005) S. C. Erwin and C. S. Hellberg, Surface Science 585, L171 (2005).
  * Bruch (2005) L. W. Bruch, Surface Science 585, 135 (2005).


# Graph Classification Based on Skeleton and Component Features

Xue Liu  Wei Wei<EMAIL_ADDRESS>Xiangnan Feng  Xiaobo Cao  Dan Sun  Beijing
System Design Institute of Electro-Mechanic Engineering, Beijing, 100854,
China  School of Mathematical Sciences, Beihang University, Beijing, 100191,
China  Key Laboratory of Mathematics, Informatics and Behavioral Semantics,
Ministry of Education, 100191, China  Beijing Advanced Innovation Center for
Big Data and Brain Computing, Beihang University, Beijing, 100191, China  Peng
Cheng Laboratory, Shenzhen, Guangdong, 518066, China

###### keywords:

Graph representation, Graph classification, Feature learning

##  1 Introduction

Graph classification to distinguish the class labels of graphs in a dataset is
an important task with practical applications in a large spectrum of fields
(e.g., bioinformatics [1], social network analysis [2] and chemoinformatics
[3]). In these areas, data can be usually represented as graphs with labels.
For example, in bioinformatics, a protein molecule can be represented as a
graph whose nodes corresponds to atoms, and edges signify there exits chemical
bonds or not between atoms. The graphs are allocated with different labels
based on having specific function or not. To make classification in this task,
we usually make a common assumption that protein molecules with similar
structure have similar functional properties.

More recently, there has been a surge of approaches that seek to learn
representations or embeddings that encode features about the graphs and then
make classification. The idea behind these learning approaches focuses on
graph structure representation and learning a mapping that embeds nodes or
entire (sub)graphs, into a low-dimensional vector. Most of these methods can
be classified into two categories: (1) neural networks manners [4] that learn
the large-scale structures of target graph, (2) kernel methods [5] that learn
small-size structures of target graph. Different structures of graph imply
dissimilar features.

Graph neural networks (GNNs) [6] use a recurrent network framework to transmit
information from a calculated node to another new node until reaching a stop
situation. Analogous to image-based convolutional neural networks (CNNs) [7],
PATCHY-SAN (PSCN) [4] is motivated to operate on locally connected regions of
the input to learn graph embeddings. Graph convolutional networks (GCNs) [8]
operate directly on graph data using spectral filters to exploit local areas
and then extract local meaningful features shared with the entire graph to get
a large-size graph structure representation. The success of neural networks
relies on enormous amount of data, and usually uses iterative calculations to
spread information, so that the local information of the graph get coupled and
is integrated into the overall embedding.

Subgraph isomorphism has been proven to be NP-complete, however graph
isomorphism problem is in NP and has been neither proven NP-complete nor could
be solved by a polynomial-time algorithm [9]. Graph kernels [5] differentiate
two graphs by recursively decomposing them into substructures and defining a
function on graph to make classification based on graphs similarity measures
in an unsupervised way. They bridge the gap between graph data and a wide
range of machine learning methods such as Support Vector Machines (SVM),
regression, clustering and Principal Components Analysis (PCA), etc. Several
different graph kernels approaches are usually divided into two classes: walk-
based patterns and limited-size subgraph methods. In random walk pattern,
graph kernels count matched random walks pairs between two graphs [10]. The
shortest path kernels count pairs of shortest paths having the same beginning
node and sink labels and the same length in two graphs [11]. Graphlet kernels
discuss graph isomorphism by counting the occurrences of all types of fixed
size subgraphs [12]. Graphs are regarded similar if they share lots of common
subgraphs. [3] sets the source nodes of two graphs at a fixed distance from
each other, and then finds subgraphs containing nodes up to a certain distance
from the root and calculate the number of identical pairs of subgraphs.
Weisfeiler-Lehman graph kernels are highly efficient kernels to express graph
isomorphism based on comparing subtree-like patterns [5].

Figure 1:  Overview of topological graphs (as shown in (a)), molecular
chemical structures (as shown in (b)), their corresponding functional groups
and subgraphs (as shown in (c)) of phenol
($\textrm{C}_{6}\textrm{H}_{5}\textrm{O}\textrm{H}$) and methylbenzene
($\textrm{C}_{6}\textrm{H}_{5}\textrm{C}\textrm{H}_{3}$). In (a), we do random
walks with walk length $l=4$ from each node on every topological graph of
phenol and methylbenzene. The frequency $p$ of each edge covered in random
walks is $p=\frac{m_{e}}{M}$, where $m_{e}$ is the frequence of edge $e$ being
covered and $M$ is the total number of edges being covered. Edges with
different frequencies are represented in different colors. Phenol and
methylbenzene in (c) can be decomposed into three functional groups (hydroxy,
benzene and methyl), which are regarded as 2-order subgraph (2-hop path),
hexagon loop (6-order subgraph) and cross separately.

There are two critical limitations of graph kernels: (1) Many of them do not
provide explicit embeddings, so that kernels are unusable for many proposed
machine learning algorithms which operate on vector directly. (2) The
substructures (i.e., random walks, subgraphs, etc.) need to be manually set
priority, for instances, the length of random walk or the shape of subgraphs.
This makes some substructures not in specified shape be ignored easily.

Anonymous Walk Embeddings (AWE) proposes a novel graph embedding method that
relies on distribution of special random walks named anonymous random walks to
get graph representations [13]. In analogy to graph kernels and to avoid
sparse distribution, AWE method uses random walks in anonymous manner to catch
“skeletons” of the whole graph. [14] has shown that anonymous walks are Markov
processes from starting node. After adequate sampling, anonymous walks are
capable to reconstruct original graph. Two graphs with similar distributions
of anonymous walks are regarded topological similar [13].

However, AWE may ignore important subgraphs which actually determine graph
properties because of low distribution. As illustrated in Figure 1, we take
phenol ($\textrm{C}_{6}\textrm{H}_{5}\textrm{O}\textrm{H}$) and methylbenzene
($\textrm{C}_{6}\textrm{H}_{5}\textrm{C}\textrm{H}_{3}$) as simple examples.
Each of them can be decomposed into two main parts, the main structure
(benzene ring) as a skeleton, and functional groups (i.e., hydroxy
“$-\textrm{O}\textrm{H}$” and methyl “$-\textrm{C}\textrm{H}_{3}$”), and both
of them determine the molecular properties and functions. Taking these two
molecules as topological graphs and doing random walks on them, experiment
shows that the frequency of benzene ring edges being covered is twice more
than that of hydroxy or methyl edges being covered. In AWE, these two graphs
are structure similar since random walks catch same main structures as
skeletons. However, phenol and methylbenzene have quite different chemical
properties determined by different functional groups, which are in low
distributions. In addition, AWE lacks a hierarchical representation of the
entire graph structure. Structures AWE gets have identical size (walk length),
so that substructures whose size not equal to anonymous walks structures can
not be represented efficiently.

Our approach. We design a novel data-driven framework named GraphCSC that
represents entire graph in low-dimension vectors containing both main
structure features and component information in a global view. Our methodology
uses anonymous walks to represent skeleton information of a graph inspired by
the success of AWE; to learn graph component representation, finding all
different order subgraphs and checking subgraph isomorphism needs great
computational complexity [15]. Thus, our strategy is to embed graph using a
distribution of special subgraphs, ie., frequent subgraphs. Frequent subgraphs
are determined by a fixed threshold hyperparameter, which indicates what kinds
of subgraphs can be regarded frequent. In frequent-based subgraphs, it not
only contains the underlying semantics within an individual graph but also the
relationships among graphs. To learn graph embedding, motivated by a novel
supervised document method PV-DBOW (paragraph vector-distributed bag of words)
[16], anonymous walks and frequent subgraphs in our model are treated as words
and each graph is regarded as a document. Two graphs are similar in embedding
space if their skeletons and components information are similar.

Our contribution. To the best of our knowledge, GraphCSC is a new framework
that learns embeddings consisting both skeleton features and component
information compared to other existing embedding approaches. GraphCSC actually
studies the representation of graph structures from horizontal and vertical
perspectives: in horizontal perspective, we attempt to characterize relatively
high-order structures using anonymous random walks with same walk length,
which determines the skeletons of graphs; in vertical perspective, our model
focuses on what kinds of subgraphs with different sizes or shapes they have.
We use a NLP (Natural Language Processing) training framework with skeletons
and components together as inputs to learn graph representation with combined
information. Through empirical evaluation on multiple real-world datasets,
experiments show that our model is competitive than various established
baselines.

##  2 PROBLEM STATEMENT

Definition 1 (Graph Classification). Given a set of graphs with labels,
$\mathbb{G}$ = $\\{\emph{G}_{1},\ldots,\emph{G}_{N}\\}$, the goal is to learn
a function $\varphi:\mathbb{G}\rightarrow\mathbb{L}$, where $\mathbb{G}$ is
the input space of graphs and $\mathbb{L}$ is the set of graph labels. Each
weighted graph _G_ i with label $l_{i}$, $i=1,\ldots,N$, is a tuple $G_{i}$ =
$(V_{i},E_{i},\Omega_{i})$, where $V_{i}$ is the set of $n_{i}$ vertices from
$G_{i}$, $E_{i}\subseteq V_{i}\times V_{i}$ is the set of edges from $G_{i}$,
and $\Omega_{i}$ is the set of edge weights from $G_{i}$.

Definition 2 (Graph Embedding). Given a set of graphs $\mathbb{G}$ =
$\\{G_{1},\ldots,G_{N}\\}$, the goal is to learn a graph embedding matrix
$\boldsymbol{X}_{N\times d}$, where each $i$-th row is a $d$-dimensions vector
of graph $G_{i}$ which is learned by mapping
$\psi:\mathbb{G}\rightarrow\mathbb{R}^{1\times d}$. Graph embeddings capture
the graph similarity between $G_{i}$ and $G_{j}$ in the sense that vector
$\psi(G_{i})$ and $\psi(G_{j})$ are close in embedding space.

##  3 BACKGROUND

We will leverage two techniques Skipgram and PV-DBOW which have achieved
success in NLP to learn graph representation. Before we propose our approach,
we review these powerful models and state them as background of our model.

###  3.1 word2vec and Skipgram

How to get continuous-valued word vector representations is a core task in NLP
applications. word2vec [17] uses Skipgram to learn low-dimension embeddings of
words that capture rich semantic relationships between words. Skipgram maps
words contained in similar sentences to “near” positions in embedding space,
i.e., their representation vectors are similar.

Given the target word $w_{t}$ from vocabulary set $V$ and a sequence of words
$w_{1},\ldots,w_{t},\ldots,w_{T}$, a context
$w_{t-c},\ldots,w_{t},\ldots,w_{t+c}$ is defined as a fixed number of words
surrounding $w_{t}$ within a window $c$. Skipgram maximizes the co-occurrence
probability among words that appear in context:

$\sum\limits_{t=1}^{T}\log\mathbb{Pr}(w_{t-c},\ldots,w_{t+c}|w_{t}).$ (1)

The conditional probability $\mathbb{Pr}(w_{t-c},\ldots,w_{t+c}|w_{t})$ is
approximated under the following independence assumption:

$\mathbb{Pr}(w_{t-c},\ldots,w_{t+c}|w_{t})=\prod\limits_{j=t-c,j\neq{t}}^{t+c}\mathbb{Pr}(w_{j}|w_{t}).$
(2)

###  3.2 Softmax and negative sampling

To learn such a posterior distribution $\mathbb{Pr}(w_{t+j}|w_{t})$,
conventional classifier such as logistic regression requires vast
computational resources since the number of labels equals to vocabulary size
$|V|$. To avoid heavy calculation, conditional probability distribution
$\mathbb{Pr}(w_{j}|w_{t})$ is defined by a Hierarchical Softmax [18]:

$\mathbb{Pr}(w_{j}|w_{t})=\frac{\exp(\boldsymbol{w}_{t}\cdot\boldsymbol{w}_{j})}{\sum\limits_{i=1}^{|V|}\exp(\boldsymbol{w}_{t}\cdot\boldsymbol{w}_{i})},$
(3)

where $\boldsymbol{w}_{t}$ and $\boldsymbol{w}_{t+j}$ are embedding vectors of
word $w_{t}$ and $w_{t+j}$. Thus the co-occurrence probability (1) is written
as:

$\sum\limits_{t=1}^{T}\log\prod\limits_{j=t-c,j\neq{t}}^{t+c}\frac{\exp(\boldsymbol{w}_{t}\cdot\boldsymbol{w}_{j})}{\sum\limits_{i=1}^{|V|}\exp(\boldsymbol{w}_{t}\cdot\boldsymbol{w}_{i})}.$
(4)

To speed up the training process of Skipgram, negative sampling [19] method
random samples a small set of words as negative samples which are not involved
in context. Then only target word and negative samples are updated instead of
the whole words from vocabulary set $V$ in the process of iteration training.
This strategy would be efficient especially for situations when tasks face
huge computational pressure.

###  3.3 doc2vec and PV-DBOW

In analogy to word2vec, doc2vec [16] uses PV-DBOW to learn representations of
arbitrary size document in a document set. More specifically, given a document
set $\mathbb{D}=\\{D_{1},\ldots,D_{N}\\}$ with a set of words
$V=\\{w_{1},\ldots,w_{|V|}\\}$ and the target document $D_{t}\in\mathbb{D}$
contains a sequence of words $\\{w_{1},\ldots,w_{l}\\}$, the goal is to learn
a low-dimension vector $\textbf{D}_{t}$ of document $D_{t}$ by maximize the
following log probability of words $w_{1},\ldots,w_{l}$ contained in $D_{t}$:

$\sum\limits_{i=1}^{l}\log\mathbb{Pr}(w_{j}|D_{t}).$ (5)

The conditional probability $\mathbb{Pr}(w_{j}|D_{t})$ above is defined as:

$\mathbb{Pr}(w_{j}|D_{t})=\frac{\exp(\boldsymbol{D}_{t}\cdot\boldsymbol{w}_{j})}{\sum\limits_{i=1}^{|V|}\exp(\boldsymbol{D}_{t}\cdot\boldsymbol{w}_{i})},$
(6)

where $\boldsymbol{D}_{t}$ and $\boldsymbol{w}_{j}$ are corresponding
representation of $D_{t}$ and $w_{j}$, $|V|$ is the number of all words across
all documents in $\mathbb{D}$. The log probability (5) could be approximated
efficiently using negative sampling.

##  4 PROPOSED MODEL

Our proposed model GraphCSC has two main modules, skeleton module and
component module as summarized in Figure 5. Skeleton module (in Algorithm 1)
and component module (in Algorithm 2) mine corresponding substructures
separately but synchronously. Finally, GraphCSC integrates these two modules
and optimizes overall loss function by gradient descent.

###  4.1 Skeleton Module

####  4.1.1 Anonymous Random Walks

Definition 3 (Random Walk). In graph $G$, a random walk $w$ is defined as a
finite sequence $(v_{0},v_{1},\dots,v_{l})$ with length $l$, where $v_{0}$ is
the root node and node $v_{i+1}$ is sampled independently among the neighbors
of node $v_{i}$.

Random walks are regarded as Markov processes, recently anonymous random walks
have been proven to be capable to learn graphs structural properties and
reconstruct graphs with full descriptions of every node’s state appeared in
the random walk process in its own label space instead of global label space
[14]. This makes random walks in anonymous experiment more flexible and
compact.

Definition 4 (Position Function). Let $w=(v_{0},v_{1},\ldots,v_{l})$ be a
random walk with length _l_ on graph. The position function _f_ is defined as
$f(v_{i})=|(v_{0},\ldots,v_{i^{\prime}})|$, where $i’$ is the smallest integer
such that $v_{i^{\prime}}=v_{i}$.

Definition 5 (Anonymous Random Walk). Let $w$ = $(v_{0},v_{1},\ldots,v_{l})$
be a random walk with length _l_ on graph _G_ i. The corresponding anonymous
random walk is defined as a sequence of integers
$a=g(w)=(f(v_{0}),f(v_{1}),\ldots,f(v_{l}))$, where _f_ is the position
function.

Figure 2:  We sample four random walks with length $3$ (marked in different
colors) in $G_{1}$, $G_{2}$ and illustrate their relation with corresponding
anonymous random walks.

Example. Random walks $(A,B,C,D)$, $(B,C,D,A)$, $(A,A,B,C)$ and $(A,A,D,C)$
sampled from $G_{1}$ and $G_{2}$ in Figure 2 are completely different.
However, in anonymous random walks view, the new label for each node is
redefined as the position of the first occurrence of node with same label in
the random walk sequence. This makes the initial four different random walks
be changed into two anonymous random walks, $(1,2,3,4)$ and $(1,1,2,3)$.

As shown in Figure 3, all different random walks with length $3$ in $G_{1}$
and $G_{2}$ from Figure 2 will be converted into only $8$ anonymous random
walks. Random walks record the labels information of nodes traveled, under
sufficient samples, they can accurately capture structure traits and
reconstructure the original graph. But over precise information will not be
efficient to capture graph structure features because of sparse distribution
over all random walks. Anonymous random walks translate random walks into a
sequence of integers recording first appearing positions.

Figure 3:  All kinds of different random walks with length $3$ in $G_{1}$ and
$G_{2}$ from Figure 2 and their corresponding anonymous random walks.

####  4.1.2 Skeleton

We draw independently a set of $\xi$ random walks
$W_{G_{i}}=\\{w_{1},\ldots,w_{\xi}\\}$ with length $l$ on $G_{i}$ in graphs
set $\mathbb{G}=\\{G_{1},\ldots,G_{i},\ldots,G_{N}\\}$, and calculate its
$\mu$ anonymous random walks $A_{G_{i}}=\\{a_{1},\ldots a_{\mu}\\}$.

For a large graph $G_{i}$, to count all possible anonymous random walks needs
vast computational resources. However, sampling $\zeta$ random walks with
length $l$ to approximate actual distribution of anonymous random walks, the
overall computing running time will be $\emph{O}(\zeta l)$ [16]. The relation
of estimation of samples number $\zeta$ and the number of anonymous random
walks $\lambda$ is determined by hyperparameters $\varepsilon$ and $\delta$
[20]:

$\zeta=[\frac{2}{\varepsilon^{2}}(log(2^{\lambda}-2)-log(\delta))],\varepsilon>0,0\leq\delta\leq
1.$ (7)

Next we tend to leverage vectors to represent whether some specific anonymous
random walks contained in a graph or not.

Definition 6 (Skeleton). Given graphs set
$\mathbb{G}=\\{G_{1},\ldots,G_{i},\ldots,G_{N}\\}$, we sample random walks
with length $l$ on each graph in $\mathbb{G}$ and integrate them as set
$W=\\{w_{1},\ldots,w_{m},\ldots,w_{\xi}\\}$. The correspond anonymous random
walks set $A=\\{a_{1},\ldots,a_{s},\ldots,a_{\mu}\\}$ has $\mu$ unequal
elements. The skeleton of graph $G_{i}$ is defined as an $1\times\mu$ shape
vector $\boldsymbol{s}_{G_{i}}=[\xi_{s}]_{1\times\mu}$, where

$\quad\xi_{s}=\left\\{\begin{aligned} 1&,\ if\ a_{s}\in G_{i}\\\ 0&,\
else\end{aligned}\right.,s=1,\ldots,\mu.$ (8)

Algorithm 1 Skeleton Module

0: $\mathbb{G}$ = $\\{G_{1},\ldots,G_{N}\\}$: graphs set; $l$: random walk
length; $T$: random walk times.

0: $A$: anonymous random walks set; $\boldsymbol{s}_{G_{i}}$, $i=1,\ldots,N$:
skeleton for each graph $G_{i}$ in $\mathbb{G}$.

1: % _Mining Anonymous Walks_

2: $A=\emptyset$

3: for each $G_{i}$ in $\mathbb{G}$ do

4: $A_{G_{i}}=\emptyset$

5: for each node $v_{j}$ in $V(G_{i})$ do

6: for $k=1$; $k<T$; $k++$ do

7: $w=RandomWalk(G_{i},v_{j},l)$

8: $a=g(w)$

9: if $a$ not in $A_{G_{i}}$ then

10: $A_{G_{i}}=A_{G_{i}}\cup\\{a\\}$

11: end if

12: end for

13: end for

14: end for

15: $A=\bigcup\limits_{i=1}^{N}A_{G_{i}}$

16: % _Getting_ skeleton _s_

17: for each $G_{i}$ in $\mathbb{G}$ do

18: initialize $\boldsymbol{s}_{G_{i}}=[s_{k}]_{1\times|A|}$, $s_{k}=0$

19: for each anonymous random walk $a_{k}$ in $A$ do

20: if $a_{k}$ in $A_{i}$ then

21: $s_{k}=1$

22: end if

23: end for

24: end for

25: return $A$; $\boldsymbol{s}_{G_{i}}$, $i=1,\ldots,N$

###  4.2 Component Module

Component module uses a more flexible approach to capture similarity between
two graphs by frequent patterns, for example, itemsets, subsequences, or
subgraphs, which appear in a data set with frequency no less than a user-
specified minimum support ($min\\_sup$) threshold. Component module basically
includes three steps: (1) mining frequent subgraphs, (2) selecting features,
(3) learning component features.

####  4.2.1 Frequent Subgraph

Definition 7 (Subgraph). A graph $sg$ is called a subgraph of a graph $G$ if
nodes satisfy V($sg$) $\subseteq$ V($G$) and edges satisfy E($sg$) $\subseteq$
E($G$).

Figure 4: Different kinds of subgraphs mined from $G_{1}$ and $G_{2}$ in
Figure 2.

Two graphs are similar if their subgraphs are similar since subgraphs have
been recognized as fundamental units and are building blocks for complex
networks [21]. But usually only a fraction of the large amount number of
subgraphs are actually relevant to data mining problems. Our methodology focus
on mining subgraphs of different proportion in graph set, and then embed these
distribution features into graph component representation.

Definition 8 (Frequent Subgraph). Let $sg$ be a subgraph in subgraphs set $SG$
mined from $\mathbb{G}=\\{G_{1},\ldots,G_{i},\ldots,G_{N}\\}$, and will be
called frequent subgraph $fsg$ if $\frac{|\\{G_{i}\in\mathbb{G}|sg\in
G_{i}\\}|}{|\mathbb{G}|}$ $\geq\theta$, where $\theta$ is the $min\\_sup$
threshold, $0\leq\theta\leq 1$, and $\frac{|\\{G_{i}\in\mathbb{G}|sg\in
G_{i}\\}|}{|\mathbb{G}|}$ is called the $support$ of $G_{i}$.

Example. We still take $G_{1}$, $G_{2}$ in Figure 2 and subgraphs $sg_{1}$,
$sg_{2}$, $sg_{3}$ in Figure 4 as examples. If we set $min\\_sup$ threshold
$\theta=1$, only subgraphs $sg_{1}$ can be regarded as frequent subgraph.
However, if we set $min\\_sup$ threshold $\theta=0.5$, all of $sg_{1}$,
$sg_{2}$ and $sg_{3}$ will be regarded frequent.

Frequent pattern, as a form of non-linear feature combinations over the set of
different subgraphs, has higher discriminative power than that of single kind
of subgraph because they capture more underlying semantics of the data. The
key point is to specify the hyperparameter $min\\_sup$ threshold $\theta$ used
in model in frequent pattern, and we will study the influence of $\theta$ for
machine learning tasks in Section 5. If an infrequent feature is used, the
model cannot generalize well to the test data since it is built based on
statistically minor observations, hence the discriminative power of low-
support features will be limited [22].

Figure 5:  Overview of GraphCSC which can be decomposed into two main parts,
the skeleton module and component module. Skeleton module outputs anonymous
random walks matrix $\mathbb{M}_{A}$ in $\mu\times d$ size and skeletons for
graphs from $\mathbb{G}$. Meanwhile, component module gets frequent subgraphs
matrix $\mathbb{M}_{FSG}$ in $\nu\times d$ size and componnets for graphs from
$\mathbb{G}$. Finally, GraphCSC integrates these two modules and yields global
optimization object (14).

####  4.2.2 Component

To mine all frequent subgraphs will face two challenges: (1) The time
complexity $O(2^{|V|})$ which is dominated by the overall nodes number of each
graph, makes it computational infeasible to find all subgraphs. (2) In
analogue of graph isomorphism checking in which two graphs have different
sizes, subgraph isomorphism checking has been proven to be NP-complete [23]. A
slightly less restrictive measure of similarity can be deﬁned based on the
size of the largest common subgraph in two graphs, but unfortunately the
problem of ﬁnding the largest common subgraph of two graphs is NP-complete as
well [23].

To tackle with such two challenges above, gSpan [24] algorithm builds a new
lexicographic order among graphs, and maps each graph to a unique minimum DFS-
code as its canonical label. Based on this lexicographic order among DFS-code,
gSpan adopts the depth-first search strategy to discover frequent subgraphs
efficiently. Repeating this procedure until either when the $support$ of a
graph is less than $min\\_sup$ threshold $\theta$, or its code is not a
minimum code, which means this graph and all its descendants have been
generated and discovered before. Now we give the definition of component to
indicate what frequent subgraphs mined by gSpan graph $G_{i}$ has.

Definition 9 (Component). Given a set of graphs
$\mathbb{G}=\\{G_{1},\ldots,G_{i},\ldots,G_{N}\\}$ and $min\\_sup$ threshold
$\theta$, all different frequent subgraphs obtained by gSpan are contained in
set $FSG$ = $\\{fsg_{1},\ldots,fsg_{t},\ldots,fsg_{\nu}\\}$. The component of
$G_{i}$ is defined as an $1\times\nu$ shape vector
$\boldsymbol{c}_{G_{i}}=[\eta_{t}]_{1\times\nu}$, where

$\quad\\\ \eta_{t}=\left\\{\begin{aligned} 1&,\ if\ fsg_{t}\in G_{i}\\\ 0&,\
else\end{aligned}\right.,t=1,\ldots,\nu.$ (9)

Algorithm 2 Component Module

0: $\mathbb{G}$ = $\\{G_{1},\ldots,G_{N}\\}$: graphs set; $\theta$:
$min\\_sup$ threshold.

0: $FSG$: frequent subgraphs set; $\boldsymbol{c}_{G_{i}}$, $i=1,\ldots,N$:
component for each graph $G_{i}$ in $\mathbb{G}$

1: % _Mining Frequent Subgraphs_

2: $FSG=\textbf{gSpan}(\mathbb{G},\theta)$

3: for each $G_{i}$ in $\mathbb{G}$ do

4: $FSG_{G_{i}}=\emptyset$

5: for each frequent subgraph $fsg_{j}$ in $FSG$ do

6: if $fsg_{j}$ in $G_{i}$ then

7: $FSG_{G_{i}}=FSG_{G_{i}}\cup\\{fsg_{j}\\}$

8: end if

9: end for

10: end for

11: % _Getting_ component _s_

12: for each $G_{i}$ in $\mathbb{G}$ do

13: initialize $\boldsymbol{c}_{G_{i}}=[c_{k}]_{1\times|FSG|}$, $c_{k}=0$

14: for each frequent subgraph $fsg_{k}$ in $FSG$ do

15: if $fsg_{k}$ in $FSG_{i}$ then

16: $s_{k}=1$

17: end if

18: end for

19: end for

20: return $FSG$; $\boldsymbol{c}_{G_{i}}$, $i=1,\ldots,N$

###  4.3 TRAINING

For graph data set $\mathbb{G}=\\{G_{1},\ldots,G_{i},\ldots,G_{N}\\}$,
supposed that all anonymous random walks with length $l$ mined are in
$A=\\{a_{1},\ldots,a_{\mu}\\}$ and all frequent subgraphs mined under
$min\\_sup$ threshold $\theta$ are in $FSG=\\{fsg_{1},\ldots,fsg_{\nu}\\}$ ,
we leverage PV-DBOW to learn graphs embedding matrix $\mathbb{X}$ whose each
row is an embedding for each graph, anonymous random walks matrix
$\mathbb{M}_{A}$ whose each row is an embedding for each anonymous random walk
and frequent subgraphs matrix $\mathbb{M}_{FSG}$ whose each row is an
embedding for each frequent subgraph. Vectors $\boldsymbol{s}_{G_{i}}$ and
$\boldsymbol{c}_{G_{i}}$ are corresponding skeleton and component for graph
$G_{i}$, $i=1,\ldots,N$.

Graphs embedding matrix $\mathbb{X}$ has $N\times d$ size, matrix
$\mathbb{M}_{A}$ is in $\mu\times d$ size and matrix $\mathbb{M}_{FSG}$ is in
$\nu\times d$ size, where $d$ is embedding dimension. PV-DBOW treats graph
data set $\mathbb{G}$ as a documents set, each graph $G_{i}$ in $\mathbb{G}$
as a document and each substructure mined in $G_{i}$ as a word contained in a
document.

To begin with, we focus on anonymous random walks and seek to optimize the
following objective function, which maximizes the log-probability of
predicting anonymous random walks that appear in graph $G_{i}$:

$\max\sum\limits_{s=1}^{\mu}\boldsymbol{s}_{G_{i}}\cdot
e^{{}^{\prime}}_{s}\log\mathbb{Pr}(a_{s}|\emph{G}_{i}),$ (10)

where $e_{s}^{{}^{\prime}}$ is an $1\times\mu$ binary vector whose $s$-th
column element is $1$ and $j$-th column element equals to $0$ if $j\neq s$.

The probability $\mathbb{Pr}(a_{s}|\emph{G}_{i})$ is defined as a softmax unit
parametrized by a dot product of $\boldsymbol{a}_{s}$ and $\boldsymbol{G}_{i}$
which are embedding vectors of $a_{s}$ and $G_{i}$:

$\mathbb{Pr}(a_{s}|\emph{G}_{i})=\frac{\exp(\boldsymbol{a}_{s}\cdot\boldsymbol{G}_{i})}{\sum\limits_{p=1}^{\mu}\exp(\boldsymbol{a}_{p}\cdot\boldsymbol{G}_{i})}.$
(11)

Next, we use same method to maximize the log-probability of predicting the
frequent subgraphs that appear in graph $G_{i}$:

$\max\sum\limits_{t=1}^{\nu}\boldsymbol{c}_{G_{i}}\cdot
e^{{}^{\prime\prime}}_{t}\log\mathbb{Pr}(fsg_{t}|\emph{G}_{i}),$ (12)

where $e_{t}^{{}^{\prime\prime}}$ is an $1\times\nu$ binary vector whose
$t$-th column element is $1$ and $j$-th column element equals to $0$ if $j\neq
t$.

The probability $\mathbb{Pr}(fsg_{t}|\emph{G}_{i})$ is defined as a softmax
function parametrized by a dot product of $\boldsymbol{fsg}_{t}$ and
$\boldsymbol{G}_{i}$ which are embedding vectors of subgraph $fsg_{t}$ and
$G_{i}$:

$\mathbb{Pr}(fsg_{t}|\emph{G}_{i})=\frac{\exp(\boldsymbol{fsg}_{t}\cdot\boldsymbol{G}_{i})}{\sum\limits_{q=1}^{\nu}\exp(\boldsymbol{fsg}_{q}\cdot\boldsymbol{G}_{i})}.$
(13)

To get graph representation satisfies (10) and (12), the global optimization
object is:

$\max\sum\limits_{s=1}^{\mu}\boldsymbol{s}_{G_{i}}\cdot
e^{{}^{\prime}}_{s}\log\mathbb{Pr}(a_{s}|\emph{G}_{i})+\max\sum\limits_{t=1}^{\nu}\boldsymbol{c}_{G_{i}}\cdot
e^{{}^{\prime\prime}}_{t}\log\mathbb{Pr}(fsg_{t}|\emph{G}_{i}).$ (14)

We yield its loss function as following:

$\begin{split}\mathbb{L}=&\sum\limits_{s=1}^{\mu}\boldsymbol{s}_{G_{i}}\cdot
e^{{}^{\prime}}_{s}\log\sigma(\boldsymbol{a}_{s}\cdot\boldsymbol{G}_{i})+\sum\limits_{t=1}^{\nu}\boldsymbol{c}_{G_{i}}\cdot
e^{{}^{\prime\prime}}_{t}\log\sigma(\boldsymbol{fsg}_{t}\cdot\boldsymbol{G}_{i})\\\
&\sum\limits_{p=1}^{\mu}\log(\sigma(-\boldsymbol{a}_{p}\cdot\boldsymbol{G}_{i}))+\sum\limits_{q=1}^{\nu}\log(\sigma(-\boldsymbol{fsg}_{q}\cdot\boldsymbol{G}_{i})),\end{split}$
(15)

where $\sigma(x)=\frac{1}{1+exp(-x)}$ is sigmoid function.

The last two items in equation (15) which sum over all anonymous random walks
and subgraphs directly, are too expensive since $\mu$ and $\nu$ usually tend
to be very large. Hence we proceed with an approximation by negative sampling
to make the optimization problem tractable. The normalization terms from the
softmax are replaced by $K_{1}$ anonymous random walks negative samples
$\\{{a^{\prime}_{1},\ldots,a^{\prime}_{K_{1}}}\\}$ from $A$ but not contained
in $G_{i}$ and $K_{2}$ frequent subgraphs negative samples
$\\{{fsg^{\prime}_{1},\ldots,fsg^{\prime}_{K_{2}}}\\}$ from $F$ but not
contained in $G_{i}$. Thus equation (15) can be rewritten as:

$\begin{split}\mathbb{L}=&\sum\limits_{s=1}^{\mu}\boldsymbol{s}_{G_{i}}\cdot
e^{{}^{\prime}}_{s}\log\sigma(\boldsymbol{a}_{s}\cdot\boldsymbol{G}_{i})+\sum\limits_{t=1}^{\nu}\boldsymbol{c}_{G_{i}}\cdot
e^{{}^{\prime\prime}}_{t}\log\sigma(\boldsymbol{fsg}_{t}\cdot\boldsymbol{G}_{i})\\\
&\sum\limits_{p=0}^{K_{1}}log(\sigma(-\boldsymbol{a^{\prime}}_{p}\cdot\boldsymbol{G}_{i}))+\sum\limits_{q=0}^{K_{2}}log(\sigma(-\boldsymbol{fsg^{\prime}}_{q}\cdot\boldsymbol{G}_{i})),\end{split}$
(16)

where $\boldsymbol{a^{\prime}}_{p}$ and $\boldsymbol{fsg^{\prime}}_{q}$ are
the embedding vectors of samples $a^{\prime}_{p}$ and $fsg^{\prime}_{q}$
respectively, and $a^{\prime}_{p}$ belongs to $G_{i}$ when $p=0$,
$fsg^{\prime}_{q}$ is in $G_{i}$ when $q=0$.

Finally, we optimize loss function (16) with stochastic gradient descent and
update $\boldsymbol{G}_{i}$. After the learning process finishes, two graphs
are near in embedding space if they have similar skeleton and component, and
we summarize this training process in Algorithm 3.

Algorithm 3 Training

0: $\mathbb{G}$ = $\\{G_{1},\ldots,G_{N}\\}$: graphs set; $A$: anonymous
random walks set; $\boldsymbol{s}_{G_{i}}$, $i=1,\ldots,N$: skeleton for each
graph $G_{i}$ in $\mathbb{G}$; $FSG$: frequent subgraphs set;
$\boldsymbol{c}_{G_{i}}$, $i=1,\ldots,N$: component for each graph $G_{i}$ in
$\mathbb{G}$; $e^{{}^{\prime}}_{1},\dots,e^{{}^{\prime}}_{|A|}$: $|A|$ binary
vectors for skeletons;
$e^{{}^{\prime\prime}}_{1},\dots,e^{{}^{\prime\prime}}_{|FSG|}$: $|FSG|$
binary vectors for components; $d$: embedding dimension; $\alpha$: learning
rate.

0: $\mathbb{X}^{N\times d}$: graphs representation matrix

1: initialize $\mathbb{X}=[x_{ij}]_{N\times d}$,
$x_{ij}\sim\emph{{N}}(0,0.001)$

2: for each anonymous walk $a_{s}$ in $A$ do

3: for each frequent subgraph $fsg_{t}$ in $FSG$ do

4: $\mathbb{L}(\mathbb{X})=-\boldsymbol{s}_{G_{i}}\cdot
e^{{}^{\prime}}_{s}\log\mathbb{Pr}(a_{s}|G_{i})-\boldsymbol{c}_{G_{i}}\cdot
e^{{}^{\prime\prime}}_{t}\log\mathbb{Pr}(fsg_{t}|G_{i})$

5:
$\mathbb{X}=\mathbb{X}-\alpha\frac{\partial\mathbb{L}(\mathbb{X})}{\mathbb{X}}$

6: end for

7: end for

8: return $\mathbb{X}$

##  5 EXPERIMENTS

In this section, to quantitatively evaluate classification capability of our
model, we conduct extensive experiments on a variety of widely-used datasets
to compare with several state-of-the-art baselines.

###  5.1 DATASETS

We evaluate our proposed method on binary classification task using seven
real-world graph datasets whose statictics are summarized in Table 1. MUTAG
[25] is a dataset of aromatic and heteroaromatic nitro compounds labeled
according to whether or not they have a mutagenic effect on bacteria. PROTEINS
[26] is a set of proteins graphs where nodes represent secondary structure
elements and edges indicate neighborhood in the amino-acid sequence or in
3-dimension space. ENZYMES [26] consists of protein tertiary structures
obtained from the BRENDA enzyme database. DD [27] is a dataset of protein
structures where nodes represent amino acids and edges indicate spatial
closeness, which are classiﬁed into enzymes or non-enzymes. PTC-MR [28]
consists of graph representations of chemical molecules labeled according to
carcinogenicity on rodents. NCI-1, NCI-109 [1] are datasets of chemical
compounds divided by the anti-cancer property (active or negative). These
datasets have been made publicly available by the National Cancer Institute
(NCI).

Table 1: Statistics of the benchmark graph datasets. The columns are: Name of dataset, Number of graphs, Number of classes (maximum number of graphs in a class), Average number of nodes/edges. Datasets | Graph $\sharp$ | Class $\sharp$ | Average Node $\sharp$ | Average Edge $\sharp$
---|---|---|---|---
MUTAG | 188 | 2 | 17.93 | 19.79
DD | 1178 | 2 | 284.32 | 715.66
PTC-MR | 344 | 2 | 14.29 | 14.69
NCI-1 | 4110 | 2 | 29.87 | 32.30
NCI-109 | 4127 | 2 | 29.68 | 32.13
PROTEINS | 1113 | 2 | 39.06 | 72.82
ENZYMES | 600 | 6 | 32.63 | 64.14

###  5.2 BASELINES

In order to demonstrate the effectiveness of our proposed approach, we compare
it with several baseline methods, all of which utilize the entire graph for
feature extraction. These competitors can be categorized into four main
groups:

  * 1.

Graph kernels based methods: The shortest path (SP) [11] kernel measures the
similarity of a pair of graphs by comparing the distance of the shortest paths
between nodes in the graphs. Graphlet kernel (GK) [20] measures graph
similarity by counting the number of different graphlets and Deep GK [29] is
deep graphlet kernel. Weisfeiler-Lehman kernel (WL) [5] uses subtree pattern
to mine structure information and Deep WL [29] is deep Weisfeiler-Lehman
kernel.

  * 2.

Unsupervised graph embedding methods: node2vec [4] is an unsupervised task
agnostic method that learns entire graph embedding. It proposes every graph
into a fixed size vector containing distributed representation of graph
structures.

  * 3.

Supervised graph embedding methods: PSCN [8] is a convolutional neural network
algorithm which has achieved strong classiﬁcation accuracy in many datasets.

  * 4.

Unsupervised graph embedding methods: AWE [13] uses an anonymous random walks
approach to embed entire graphs in an unsupervised manner.

Table 2: Classiﬁcation accuracy (standard deviation) of our method GraphCSC
and state-of-the-art baselines on benchmark datasets. The last two rows denote
the $\theta$ values and the $dim$ values used by our method for each dataset;
these values are determined in Fig.5 and Fig.6. “-” means the classiﬁcation
accuracy (standard deviation) is not available in the original papers.

Algorithm | MUTAG | DD | PTC-MR | NCI-1 | NCI-109 | PROTEINS | ENZYMES
---|---|---|---|---|---|---|---
WL | 80.63 (3.07) | 77.95 (0.70) | 56.97 (2.01) | 80.13 (0.50) | 80.22 (0.34) | 72.92 (0.56) | 53.15 (1.14)
Deep WL | 82.95 (1.96) | - | 59.17 (1.56) | 80.31 (0.46) | 80.32 (0.33) | 73.30 (0.82) | 53.43 (0.91)
GK | 81.66 (2.11) | 78.45 (0.26) | 57.26 (1.41) | 62.28 (0.29) | 62.60 (0.19) | 71.67 (0.55) | 26.61 (0.99)
Deep GK | 82.66 (1.45) | - | 57.32 (1.13) | 62.48 (0.25) | 62.69 (0.23) | 71.68 (0.50) | 27.08 (0.79)
graph2vec | 83.15 (9.25) | 58.64 (0.01) | 60.17 (6.86) | 73.22 (1.81) | 74.26 (1.47) | 73.30 (2.05) | 44.33 (0.09)
PSCN | 92.63 (4.21) | 77.12 (2.41) | 60.00 (4.82) | 78.59 (1.89) | - | 75.89 (2.76) | -
AWE | 87.87 (9.76) | 71.51 (4.02) | 59.14 (1.83) | 62.72 (1.67) | 63.21 (1.42) | 70.01 (2.52) | 35.77 (5.93)
GraphCSC | 88.42 (6.47) | 89.38 (2.73) | 64.04 (3.61) | 85.70 (4.29) | 85.46 (3.66) | 76.71 (3.06) | 57.21 (5.70)
$\theta$ | 0.15 | 0.25 | 0.75 | 0.20 | 0.20 | 0.60 | 0.55
dim | 128 | 8 | 128 | 16 | 16 | 32 | 128

Figure 6: Here we show parameter sensitivity of $min\\_sup$ threshold $\theta$
in graph classification on MUTAG, DD, PTC-MR, NCI-1, PROTEINS and ENZYMES with
fixed embedding dimension $dim=128$. The best expressions marked by red points
are used in our experiments and the $(\pm)$standard deviation of each result
is indicated with gray error bar. Figure 7: Parameter sensitivity of embedding
dimension $dim$ in graph classification on MUTAG, DD, PTC-MR, NCI-1, PROTEINS
and ENZYMES with best $min\\_sup$ threshold $\theta$ selected in Fig.5. The
best expressions marked by red points are used in our experiments. The
$(\pm)$standard deviation of each result is indicated with gray error bar.

###  5.3 Evaluation Metrics

To evaluate the performance of GraphCSC, we randomly split the data into $10$
roughly equal-size parts and perform a $10$-fold cross validation on each
dataset, in which $9$ folds for training and $1$ fold as validation for
testing. This process is repeated $10$ times and an average accuracy is
reported as prediction. Since we focus on graph embedding not on a classifier,
we feed the embedding vectors to Support Vector Machine (SVM) with RBF kernel
function with parameter $\sigma$ varing from the range
$[10^{-4},10^{-3},10^{-2},10^{-1},1,10]$.

###  5.4 Parameter sensitivity

Since we tend to prove that the performance of our model GraphCSC may have a
gain on experiments results than method only considers skeletons, thus to
explore how the hyperparameters in skeleton module affect tasks performace
will not be our priority. For the need of brevity, the length $l$ is set as
$10$ to generate a corpus of co-occurred anonymous random walks for all given
datasets. To approximate actual distribution of anonymous random walks, we
follow [13] to set the sampling hyperparameters $\varepsilon=1$ and
$\delta=0.05$. In order to evaluate how the parameters sensitivity of
$min\\_sup$ threshold $\theta$ and embedding dimensions $dim$ affect the
classification performance of GraphCSC on datasets, all parameters are assumed
to be default except $\theta$ and $dim$.

We first conduct experiments with $dim=128$, then assess the classification
accuracy as a function of $min\\_sup$ threshold $\theta$ for different
datasets. Best performance is indicated with red mark in Figure 6. Then in
Figure 7, experiments examine the influence of varying the $dim$ from
$[8,16,32,64,128,256,512,1024]$ with the best $\theta$s obtained from Fig.5.
The best accuracy is also signed in red as shown in each figure.

###  5.5 Results and Discussion

The average classification accuracy (standard deviation) over 10-fold cross-
validation of GraphCSC (with $\theta$ and $dim$ that lead to optimal result
from Figure 6 and Figure 7, and baselines on seven real-world datasets are
summarized in Table 2. From the results, it is evidently shown that GraphCSC
is always the best in terms of performance on 6 datasets with exception on
MUTAG, where GraphCSC gets second best result. More specifically, the proposed
model achieves $3.41\%-30.60\%$ improvement over graph kernals based methods
(WL, Deep WL, GK and Deep GK), and is competitive against graph embedding
method (graph2vec) with $3.41\%-30.74\%$ gain in accuracy. However, supervised
graph embedding method PSCN is more superior with 4.21% higher in accuracy in
MUTAG dataset classification. It is also obvious that GraphCSC outperforms AWE
approach in every single dataset. This significantly demonstrates the
effectiveness of the proposed model on classification tasks due to having
access to skeleton and component features, which enable GraphCSC to get more
complex structures information, while AWE processes only fixed size skeletons
instead.

Figure 6 presents the classification accuracy of increasing $min\\_sup$
threshold $\theta$. In MUTAG and NCI-1, experiments get best performance when
$\theta=0.15$ and $0.2$ separately. Then results decrease slowly and maintain
stable, after $\theta$ reaching $0.7$ (for MUTAG) and $0.6$ (for NCI-1), lines
drop rapidly until getting to another flat states. A common phenomenon is that
for DD, PROTEINS and ENZYMES, accuracy lines increase from the beginnings and
decrease slowly after meeting the tops. Another observation is that in PTC-MR,
performance has 2 closed top best results when $\theta=0.25$ and $0.75$, while
the result at $\theta=0.75$ has only $1\%$ gain over that at $\theta=0.25$. It
is very interesting that $min\\_sup$ threshold $\theta$ centers around a low
value for most datasets. A possible explanation is that low $\theta$ will lead
to more sufficient frequent subgraphs which will provide adequate
complementary structure information for pattern just uses skeletons. As a
consequence, this enforces model’s discriminate power.

Figure 7 examines the effects of embedding dimension hyperparameter $dim$. As
we can see from the figure, in MUTAG and PROTEINS, there is a relatively
gentle ﬁrst-increasing and then-decreasing accuracy line when $dim$ increases.
We note that for DD, the performance line goes down slowly and maintains
stable at around $83\%$ and then again falls quickly. Accuracy in ENZYMES
changes roughly, indicating that classification result significantly suffers
from embedding dimension $dim$. While it seems that $dim$ has no clearly
effects on classification results for PTC-MR and NCI-1. The beat $dim$s for
MUTAG and ENZYMES center at 128, but the $dim$ values tend to be small
(usually no more than 32) in DD, PTC-MR, NCI-1 and PROTEINS, we can infer that
small $dim$s would reflect more obvious features for these 4 datasets.

##  6 Conclusion

In this paper, we focus to study the problem of graph classification and
propose GraphCSC, a new methodology that learns graph representation from
horizontal and vertical perspectives to mine graph structures, i.e., skeletons
and components. And our model is demonstrated to be superior to approach which
utilizes only skeletons. Our method treats a graph as a document and different
substructures or subgraphs in it as words by NLP framework PV-DBOW, thus the
graph embedding GraphCSC learns integrate various size structures information.
Several real-world graph classification tasks show that our model can achieve
promising performances than a list of state-of-art baselines.

For future work, an interesting problem needs to be study further. We tend to
fucus on more complex graph structures, for example, clusters. Since a cluster
usually contains dozens or even hundreds of nodes and edges, there may show
more valuable features which may reflect more implicit graph properties. We
would also like to investigate graph embedding with mixed structures from low-
order subgraphs to big clusters, and to verify whether this would be helpful
for better results than GraphCSC.

### Acknowledgements

This work is supported by the Fundamental Research Funds for the Central
Universities, the National Natural Science Foundation of China (No.11201019),
the International Cooperation Project No.2010DFR00700, Fundamental Research of
Civil Aircraft No. MJ-F-2012-04, the Beijing Natural Science Foundation
(1192012, Z180005) and National Natural Science Foundation of China
(No.62050132).

## References

  * [1] J. B. Lee, R. Rossi, X. Kong, Graph classification using structural attention, in: Proceedings of the 24th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining, 2018, pp. 1666–1674.
  * [2] L. Backstrom, J. Leskovec, Supervised random walks: predicting and recommending links in social networks, in: the Fourth International Conference on Web Search and Web Data Mining, 2011, p. 635–644.
  * [3] F. Costa, K. D. Grave, Fast neighborhood subgraph pairwise distance kernel, in: International Conference on International Conference on Machine Learning, 2010\.
  * [4] M. Niepert, M. Ahmed, K. Kutzkov, Learning convolutional neural networks for graphs, in: International conference on machine learning, 2016, pp. 2014–2023.
  * [5] N. Shervashidze, P. Schweitzer, E. Jan, V. Leeuwen, K. M. Borgwardt, Weisfeiler-lehman graph kernels, Journal of Machine Learning Research 1 (3) (2010) 1–48.
  * [6] Z. Wu, S. Pan, F. Chen, G. Long, P. S. Yu, A comprehensive survey on graph neural networks, IEEE Transactions on Neural Networks and Learning Systems PP (99) (2020) 1–21.
  * [7] A. Krizhevsky, I. Sutskever, G. E. Hinton, Imagenet classification with deep convolutional neural networks, in: Advances in neural information processing systems, 2012, pp. 1097–1105.
  * [8] T. N. Kipf, M. Welling, Semi-supervised classification with graph convolutional networks, arXiv preprint arXiv:1609.02907.
  * [9] J. Hartmanis, Computers and intractability: a guide to the theory of np-completeness (michael r. garey and david s. johnson), Siam Review 24 (1) (1982) 90.
  * [10] H. Kashima, K. Tsuda, A. Inokuchi, Marginalized kernels between labeled graphs, in: Proceedings of the 20th international conference on machine learning (ICML-03), 2003, pp. 321–328.
  * [11] K. M. Borgwardt, H.-P. Kriegel, Shortest-path kernels on graphs, in: Fifth IEEE international conference on data mining (ICDM’05), IEEE, 2005, pp. 8–pp.
  * [12] R. Kondor, N. Shervashidze, K. M. Borgwardt, The graphlet spectrum, in: Proceedings of the 26th Annual International Conference on Machine Learning, 2009, pp. 529–536.
  * [13] S. Ivanov, E. Burnaev, Anonymous walk embeddings, arXiv preprint arXiv:1805.11921.
  * [14] S. Micali, Z. A. Zhu, Reconstructing markov processes from independent and anonymous experiments, Discrete Applied Mathematics 200 (2016) 108–122.
  * [15] D. Nguyen, W. Luo, T. D. Nguyen, S. Venkatesh, D. Phung, Learning graph representation via frequent subgraphs, in: Proceedings of the 2018 SIAM International Conference on Data Mining, SIAM, 2018, pp. 306–314.
  * [16] Q. Le, T. Mikolov, Distributed representations of sentences and documents, in: International conference on machine learning, 2014, pp. 1188–1196.
  * [17] T. Mikolov, Distributed representations of words and phrases and their compositionality, Advances in Neural Information Processing Systems 26 (2013) 3111–3119.
  * [18] A. Mnih, G. E. Hinton, A scalable hierarchical distributed language model, in: Advances in neural information processing systems, 2009, pp. 1081–1088.
  * [19] F. Rousseau, E. Kiagias, M. Vazirgiannis, Text categorization as a graph classification problem, in: Proceedings of the 53rd Annual Meeting of the Association for Computational Linguistics and the 7th International Joint Conference on Natural Language Processing (Volume 1: Long Papers), 2015, pp. 1702–1712.
  * [20] N. Shervashidze, S. Vishwanathan, T. Petri, K. Mehlhorn, K. Borgwardt, Efficient graphlet kernels for large graph comparison, in: Artificial Intelligence and Statistics, 2009, pp. 488–495.
  * [21] R. Milo, S. Shen-Orr, S. Itzkovitz, N. Kashtan, D. Chklovskii, U. Alon, Network motifs: simple building blocks of complex networks, Science 298 (5594) (2002) 824–827.
  * [22] C. Hong, Towards accurate and efficient classification: A discriminative and frequent pattern-based approach, Ph.D. thesis, University of Illinois at Urbana-Champaign. (2008).
  * [23] M. J. Zaki, W. Meira, Data mining and analysis: fundamental concepts and algorithms, Cambridge University Press, 2014.
  * [24] X. Yan, J. Han, gspan: Graph-based substructure pattern mining, in: 2002 IEEE International Conference on Data Mining, 2002. Proceedings., IEEE, 2002, pp. 721–724.
  * [25] A. K. Debnath, R. L. L. D. Compadre, G. Debnath, A. J. Shusterman, C. 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) 786–797.
  * [26] K. M. Borgwardt, O. C. Soon, S. Stefan, S. V. N. Vishwanathan, A. J. Smola, K. Hans-Peter, Protein function prediction via graph kernels, Bioinformatics 21 (suppl_1) (2005) i47–i56.
  * [27] P. D. Dobson, A. J. Doig, Distinguishing enzyme structures from non-enzymes without alignments, Journal of Molecular Biology 330 (4) (2003) 771–783.
  * [28] C. Helma, R. D. King, S. Kramer, A. Srinivasan, The predictive toxicology challenge 2000-2001, Bioinformatics 17 (1) (2001) 107–108.
  * [29] P. Yanardag, S. Vishwanathan, Deep graph kernels, in: Proceedings of the 21th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, 2015, pp. 1365–1374.