$\matrixelement{\mathrm{GHZ}_{\pm}}{e_{\bm{m}}}{\varphi_{\pm}(g)}$ are
exponentially small themselves. However, one could argue that they are
parametrically larger than the error term in Eq. (184). Let
$k=|\text{Supp}(e_{\bm{m}})|$, and assume $k/N=a$ for $0<a<\frac{1}{2}$. It is
straightforward to use (183) and Eq. (B) to show that, for both $P=\pm 1$
sectors, $c^{(k^{\prime})}_{\bm{m}}=0$ for $k^{\prime}<k$, but
$c^{(k)}_{\bm{m}}\neq 0$. This in turn suggests that
$\matrixelement{\mathrm{GHZ}_{\pm}}{e_{\bm{m}}}{\varphi_{\pm}(g)}\sim g^{k}$,
and hence,

$\displaystyle\matrixelement{\mathrm{GHZ}_{-}}{e_{\bm{m}}}{\varphi_{-}(g)}=\matrixelement{\mathrm{GHZ}_{+}}{e_{\bm{m}}}{\varphi_{+}(g)}(1+\mathcal{O}(g^{N/2-k})).$
(185)

Figure 8: The difference
$|\frac{\matrixelement{\mathrm{GHZ}_{+}}{e_{\bm{m}}}{\varphi_{+}(g)}}{\matrixelement{\mathrm{GHZ}_{-}}{e_{\bm{m}}}{\varphi_{-}(g)}}-1|$
as a function of $N$, when $e_{\bm{m}}$ is taken to be $\Pi_{j=1}^{k}X_{j}$,
for few different choices of $k$.

Fig.8 shows the absolute value of
$\frac{\matrixelement{\mathrm{GHZ}_{+}}{e_{\bm{m}}}{\varphi_{+}(g)}}{\matrixelement{\mathrm{GHZ}_{-}}{e_{\bm{m}}}{\varphi_{-}(g)}}-1$
as a function of $N$, when $e_{\bm{m}}=\Pi_{j=1}^{k}X_{j}$ for few different
choices of $k=O(N)$. As is clear from the figure, the difference goes to zero
exponentially fast in system size. A similar line of argument shows that

$\displaystyle\innerproduct{\varphi_{+}(g)}{\varphi_{+}(g^{\prime})}=\innerproduct{\varphi_{-}(g)}{\varphi_{-}(g^{\prime})}+\mathcal{O}(g_{\text{max}}^{N}),$
(186)

where $g_{\text{max}}=\max(g,g^{\prime})$. To see this, note that the first
term in the perturbative expansion of $\ket{\varphi_{+}(g)}$ that is has a
different coefficient from the corresponding term in $\ket{\varphi_{-}(g)}$,
is $c_{\bm{m}}^{(N/2)}$ with $k=|\text{Supp}(e_{\bm{m}})|=N/2$. Therefore the
first term in the perturbative expansion of
$\innerproduct{\varphi_{+}(g)}{\varphi_{+}(g^{\prime})}$ which has a different
value from the corresponding term in the expansion of
$\innerproduct{\varphi_{-}(g)}{\varphi_{-}(g^{\prime})}$ is of order
$(g^{\prime}g)^{N/2}$. At the next order, one can see that
$c_{\bm{m}}^{(N/2+1)}$ would be different for odd and even parity sectors,
even for $k=|\text{Supp}(e_{\bm{m}})|=N/2-1$. This then result in a terms of
order $g^{\prime N/2-1}g^{N/2+1}$ and $g^{\prime N/2+1}g^{N/2-1}$ in the
expansion of
$\innerproduct{\varphi_{+}(g)}{\varphi_{+}(g^{\prime})}-\innerproduct{\varphi_{-}(g)}{\varphi_{-}(g^{\prime})}$.
Continuing the same line of argument up to $N$th order shows that the terms
that appear in the perturbative expansion of
$\innerproduct{\varphi_{+}(g)}{\varphi_{+}(g^{\prime})}-\innerproduct{\varphi_{-}(g)}{\varphi_{-}(g^{\prime})}$,
are of the order $g^{\prime N/2-m}g^{N/2+m}$ for $m=0,\cdots,N/2$. For fixed
$g$,$g^{\prime}$ and large $N$, one can see that the largest contribution
comes from $g_{\text{max}}^{N}$. Fig.9 shows
$|\innerproduct{\varphi_{+}(g)}{\varphi_{+}(g^{\prime})}-\innerproduct{\varphi_{-}(g)}{\varphi_{-}(g^{\prime})}|$
computed numerically for different choices of $g$ and $g^{\prime}$, which
clearly shows the exponential drop with $N$.

Figure 9: The difference
$|\innerproduct{\varphi_{+}(g)}{\varphi_{+}(g^{\prime})}-\innerproduct{\varphi_{-}(g)}{\varphi_{-}(g^{\prime})}|$
as a function of $N$ for few different choices of $g$ and $g^{\prime}$.

##  Appendix C Adiabatic Drift and Recovery for the Non-Integrable TFIC

Here we consider the problem of recoverability in the phase diagram of the
following Hamiltonian,

$\displaystyle H(g)=$
$\displaystyle-\sum_{i}Z_{i}Z_{i+1}-\sum_{i}Z_{i}Z_{i+2}$
$\displaystyle-g\sum_{i}X_{i}-g\sum_{i}X_{i}X_{i+1}.$ (187)

This model no longer maps to free fermions under JW transformation.
Nonetheless, the $X\longleftrightarrow ZZ$ duality still pins the critical
point at $b=1$, separating the ferromagnetic phase ($g<1$) from the
paramagnetic phase ($g>1$). Moreover, the ground space subspace at the fixed
point $g=0$ still coincides with the the codespace of the repetition code. We
use this model to show that the perfect recoverability in Section II.1 does
not depend on the integrability of the underlying Hamiltonian.

We use matrix product state (MPS) representation to numerically find the odd
and even parity ground states $\ket{\varphi_{\pm}(g)}$ of Hamitonian (C) on an
open chain of $L$ qubits. To fix the relative phase of
$\ket{\varphi_{\pm}(g)}$ according to the adiabatic evolution, we fix
$\innerproduct{\mathrm{GHZ}_{\pm}}{\varphi_{\pm}(g)}$ to be real, which
according to Eq. (184) only introduces an error which is exponentially small
in $L$. We uses the same initial code state
$\ket{\psi}=1/\sqrt{2}(\ket{\mathrm{GHZ}_{+}}+\ket{\mathrm{GHZ}_{-}})$, and we
use the same recovery map after the adiabatic evolution as in Section II.1. As
can be seen from Fig. 10 the recovery succeeds perfectly in the thermodynamic
limit. The inset shows the data collapse with $g_{c}=1$ and $\nu=1$, which are
the same as for the integrable case. Nonetheless, note that the fidelity at
the critical point is lower from its critical value in Fig. 2.

Figure 10: Fidelity of the recovered state with respect to the original
encoded state as a function of $g$, for the nonintegrable model described by
the Hamiltonian (C). The inset shows the data collapse of the same data
points, using $g_{c}=1$ and $\nu=1$.

##  Appendix D Computing $\mathcal{Z}_{i}(\beta)$

In this Section we compute $\mathcal{Z}_{i}(\beta)$ appearing in Section II.2.
For the $\ket{\varphi_{i}(\beta)}$ states to be normalized we have,

$\displaystyle\mathcal{Z}_{i}(\beta)=\expectationvalue{\exp(\beta\sum_{j}Z_{j})}{\varphi_{i}(0)}.$
(188)

We first compute the expectation value for $\ket{\varphi_{0}(0)}$. First note
that one may write this state as follows,

$\displaystyle\ket{\varphi_{0}(0)}$
$\displaystyle=\frac{1}{2^{\frac{M}{2}}}\prod_{s}(1+A_{s})\ket{0}$ (189)
$\displaystyle=\frac{1}{2^{\frac{M}{2}}}\sum_{\sigma\in\\{-1,1\\}^{M}}\prod_{s}A_{s}^{(1-\sigma_{s})/2}\ket{0},$
(190)

where $\ket{0}=\ket{0}^{\otimes N}$ is the state where all qubits are in the
$\ket{0}$ of Pauli-$Z$, and $M$ is the number of stars in $\Lambda$. We can
think of $\sigma_{s}$ as a classical spin sitting on the star $s$. If the
qubit $j$ neighbors the stars $\mu$ and $\nu$ as in Fig. 4(a), then,

$\displaystyle
Z_{j}\prod_{s}A_{s}^{(1-\sigma_{s})/2}\ket{0}=\sigma_{\mu}\sigma_{\nu}\prod_{s}A_{s}^{(1-\sigma_{s})/2}\ket{0},$
(191)

Therefore, any configuration of classical spins $\\{\sigma_{s}\\}$ corresponds
to a specific computational basis of qubits, with the state of the qubit given
by the product of the neighboring classical spins. Fixing the boundary
classical spins to be in the $\sigma=1$ state ensures this to be true for
boundary qubits as well. Also note that the states represented by different
classical spins are orthogonal. Therefor we find,

$\displaystyle\mathcal{Z}_{0}(\beta)=\frac{1}{2^{M}}\sum_{\sigma\in\\{-1,1\\}^{M}}\exp(\beta\sum_{\langle\mu,\nu\rangle}\sigma_{\mu}\sigma_{\nu}),$
(192)

which is simply the partition function of the 2D-classical Ising model.

As for $\mathcal{Z}_{1}(\beta)$, one can easily verify that the same
calculation works for this case as well, if one replaces $\ket{0}$ by
$\overline{X}\ket{0}$, where all qubits are in the $\ket{0}$ state, except the
ones on the support of $\overline{X}$, which are flipped to $\ket{1}$. This
then implies that for the state
$\prod_{s}A_{s}^{(1-\sigma_{s})/2}\overline{X}\ket{0}$, we have
$Z_{j}=\sigma_{\mu}\sigma_{\nu}$ everywhere except the qubits on the support
of $\overline{X}$, for which we have $Z_{j}=-\sigma_{\mu}\sigma_{\nu}$. Thus
$\mathcal{Z}_{1}(\beta)$ would be the partition function of the 2D-classical
Ising model with a domain wall inserted along the support of $\overline{X}$
(see Fig. 4(b)), meaning that the coupling along the support of $\overline{X}$
is flipped to be anti-ferromagnetic.

##  Appendix E Explicit summations of Eq. (45) and Eq. (46)

In this section, we explain how one can arrive at Eq. (45) and Eq. (46) by
exlicitly performing the sums on the left hand sides.

###  E.1 relation between $P[e]$ and the partition function of the random-bond
Ising model(RBIM)

First we describe the relation between the function $P([e])$ and the partition
function of the random-bond Ising model, as was noted in Ref.[18]. define $J$
as

$\displaystyle e^{-2J}=\frac{q}{1-q},$ (193)

so we have,

$\displaystyle P[e]=(1-q)^{N}\sum_{e^{\prime}\in[e]}e^{-2J|e^{\prime}|}.$
(194)

Note that Eq. (193) is just the equality which defines the Nishimori line[40].
For completeness we repeat Eq. (40) here as well, which gives $q$ in terms of
$\beta$,

$\displaystyle q=\frac{1-e^{-\beta}}{2},$ (195)

If $e^{\prime}\in[e]$, it means that $C=e+e^{\prime}$ represents a
$Z$-stabilizer and hence $Z(C)$ can be written as the product of plaquette
operators. If we introduce a spin $\tau_{p}=\pm 1$ on each plaquette, we can
represent any product of plaquette stabilizers by a spin configuration
$\\{\tau_{p}\\}$, such that $Z(C)=\prod_{p}B_{p}^{(1-\tau_{p})/2}$. It is easy
to see that the corresponding $C$ has the following elements,

$\displaystyle C_{j}=(1-\tau_{\mu}\tau_{\nu})/2,$ (196)

where $\mu$ and $\nu$ are the plaquettes that are neighboring the edge $j$.
Since the product of all plaquettes is identity, spin configurations
$\\{\tau_{p}\\}$ and $\\{-\tau_{p}\\}$ correspond to the same loop $C$. Given
$C$, we can construct $e^{\prime}$ as $e^{\prime}_{j}=e_{j}+C_{j}$, however
the sum here is modulo $2$. To take care of that, we introduce variable
$\eta_{j}=(1-2e_{j})=\pm 1$. It is more convenient to write $\nu_{j}$ as
$\eta_{\mu\nu}$ where $\mu$ and $\nu$ are the plaquettes neighboring the edge
$j$. Then we have,

$\displaystyle e^{\prime}_{j}$
$\displaystyle=\big{[}1-(1-2e_{j})(1-2C_{j})\big{]}/2$
$\displaystyle=(1-\eta_{\mu\nu}\tau_{\mu}\tau_{\nu})/2$ (197)

where now all the summations are ordinary sums. Accordingly we find,

$\displaystyle P([e])$
$\displaystyle=\frac{(1-q)^{N}}{2}\sum_{\tau\in\\{1,-1\\}^{M^{\prime}}}e^{-J\sum_{\langle\mu,\nu\rangle}(1-\eta_{\mu\nu}\tau_{\mu}\tau_{\nu})}$
$\displaystyle=\frac{[q(1-q)]^{N/2}}{2}\sum_{\tau\in\\{1,-1\\}^{M^{\prime}}}e^{J\sum_{\langle\mu,\nu\rangle}\eta_{\mu\nu}\tau_{\mu}\tau_{\nu}}$
$\displaystyle=\frac{[q(1-q)]^{N/2}}{2}\mathcal{Z}_{\text{RBIM}}(e,J).$ (198)

where the factor of $\frac{1}{2}$ accounts for the double redundancy in
representation of $e^{\prime}$ in terms of $\tau$ spins. The $e$ error
determines the disorder $\eta$ in $\mathcal{Z}_{\text{RBIM}}$. Note that if
$e^{\prime}\in[e]$ then
$\mathcal{Z}_{\text{RBIM}}(e,J)=\mathcal{Z}_{\text{RBIM}}(e^{\prime},J)$.

###  E.2 Relation between $\overline{\mathcal{Z}_{\text{RBIM}}(e,J)}$ and
$\mathcal{Z}_{\text{Ising}}$

Here we related the average of partition function of random bond Ising model
the partition function of Ising model. It goes as this,

$\displaystyle\overline{\mathcal{Z}_{\text{RBIM}}(e,J)}=\overline{\sum_{\tau\in\\{1,-1\\}^{M^{\prime}}}e^{J\sum_{\langle\mu,\nu\rangle}\eta_{\mu\nu}\tau_{\mu}\tau_{\nu}}}$
(199)
$\displaystyle=\sum_{\tau\in\\{1,-1\\}^{M^{\prime}}}\prod_{{\langle\mu,\nu\rangle}}\overline{\cosh(J)+\sinh(J)\eta_{\mu\nu}\tau_{\mu}\tau_{\nu}}$
(200)
$\displaystyle=\sum_{\tau\in\\{1,-1\\}^{M^{\prime}}}\prod_{{\langle\mu,\nu\rangle}}\cosh(J)+\sinh(J)(1-2q)\tau_{\mu}\tau_{\nu}$
(201)
$\displaystyle=\Bigg{[}\frac{\cosh(J)}{\cosh(\tilde{J})}\Bigg{]}^{N}\sum_{\tau\in\\{1,-1\\}^{M^{\prime}}}e^{\tilde{J}\sum_{\langle\mu,\nu\rangle}\tau_{\mu}\tau_{\nu}}$
(202)
$\displaystyle=\Bigg{[}\frac{\cosh(J)}{\cosh(\tilde{J})}\Bigg{]}^{N}\mathcal{Z}_{\text{Ising}}(\tilde{J}).$
(203)

with $\tilde{J}$ defined as,

$\displaystyle\tanh(\tilde{J})=(1-2q)\tanh(J).$ (204)

It is also straightforward to see that

$\displaystyle\overline{\mathcal{Z}_{\text{RBIM}}(e+\ell_{1},J)}=\Bigg{[}\frac{\cosh(J)}{\cosh(\tilde{J})}\Bigg{]}^{N}\mathcal{Z}_{\text{Ising},\ell_{1}}(\tilde{J}).$
(205)

where $\mathcal{Z}_{\text{Ising},\ell_{1}}$ is the Ising model with a domain
wall along $\ell_{1}$, meaning that the coupling along $\ell_{1}$ are anti-
ferromagnetic. This is simply because
$\mathcal{Z}_{\text{RBIM}}(e+\ell_{1},J)$ can be thought of as
$\mathcal{Z}_{\text{RBIM}}(e,J)$ but where interactions along the links on
$\ell_{1}$ are $-\eta_{\mu\nu}\tau_{\mu}\tau_{\nu}$ instead of
$+\eta_{\mu\nu}\tau_{\mu}\tau_{\nu}$.

###  E.3 The relation between $\mathcal{Z}_{\text{Ising}}$ and
$\mathcal{Z}_{0}$ and $\mathcal{Z}_{1}$

(a)

(b)

Figure 11: (a) a loop configuration $c\in\mathcal{L}_{\text{even}}$ and (b)
the corresponding loop configuration
$c^{\prime}=c+\ell_{1}\in\mathcal{L}_{\text{odd}}$ where $\ell_{1}$ is taken
to be the line that runs vertically in the middle.

Note that $\mathcal{Z}_{\text{Ising}}$ is the partition function for $\tau$
spins which has free open boundary conditions on top and bottom of cylinder.
But $\mathcal{Z}_{0}$ and $\mathcal{Z}_{1}$ are partition function for
$\sigma$ spins, with fixed boundary conditions on top and bottom. As we will
see these partition functions are related to each other via the high-
temperature low-temperature duality of the Ising model. We start by high-
temperature expansion of $\mathcal{Z}_{0}$,

$\displaystyle\mathcal{Z}_{0}(\beta)=\frac{1}{2^{M}}\sum_{\sigma\in\\{-1,1\\}^{M}}\exp(\beta\sum_{\langle\mu,\nu\rangle}\sigma_{\mu}\sigma_{\nu})$
(206)
$\displaystyle=\frac{1}{2^{M}}\sum_{\sigma\in\\{-1,1\\}^{M}}\prod_{\langle\mu,\nu\rangle}\Bigg{[}\cosh(\beta)+\sinh(\beta)\sigma_{\mu}\sigma_{\nu}\Bigg{]}$
(207)
$\displaystyle=\frac{1}{2^{M}}\cosh(\beta)^{N}2^{M}\sum_{c\in\mathcal{L}}\tanh(\beta)^{|c|}$
(208)

where $\mathcal{L}$ is the set of loop configurations with no open end inside
the bulk, but they can end on the boundaries (see Fig. 11). This is because we
are not summing over the spins at the boundary so those terms are also allowed
in addition to loop configurations. We can partition $\mathcal{L}$ into
$\mathcal{L}=\mathcal{L}_{\text{even}}\cup\mathcal{L}_{\text{odd}}$ where
$\mathcal{L}_{\text{even}}$ ($\mathcal{L}_{\text{odd}}$) is the set of allowed
loop configurations with even (odd) lines at each boundary (see Fig. 11(a) and
Fig. 11(b)). Therefore we have,

$\displaystyle\mathcal{Z}_{0}(\beta)=\cosh(\beta)^{N}\Bigg{[}\sum_{c\in\mathcal{L}_{\text{even}}}\tanh(\beta)^{|c|}+\sum_{c\in\mathcal{L}_{\text{odd}}}\tanh(\beta)^{|c|}\Bigg{]}$
(209)
$\displaystyle=\cosh(\beta)^{N}\Bigg{[}\sum_{c\in\mathcal{L}_{\text{even}}}\tanh(\beta)^{|c|}+\sum_{c\in\mathcal{L}_{\text{even}}}\tanh(\beta)^{|c+\ell_{1}|}\Bigg{]}$
(210)

On the other hand, if one consider the partition function
$\mathcal{Z}_{\text{Ising}}(\tilde{\beta})$ of the Ising model on $\tau$ spins
that live on the plaquettes, with free boundary conditions on the top and
bottom, the domain wall of each $\\{\tau\\}$ spin configuration correspond to
a loop configuration $c\in\mathcal{L}_{\text{even}}$, whose energy
$\tilde{\beta}\sum_{\langle\mu,\nu\rangle}\tau_{\mu}\tau_{\nu}=\tilde{\beta}\,N-2\tilde{\beta}|c|$.Therefore,

$\displaystyle\mathcal{Z}_{\text{Ising}}(\tilde{\beta})$
$\displaystyle=\sum_{\tau\in\\{1,-1\\}^{M^{\prime}}}e^{\tilde{J}\sum_{\langle\mu,\nu\rangle}\tau_{\mu}\tau_{\nu}}$
(211)
$\displaystyle=2e^{\tilde{\beta}\,N}\sum_{c\in\mathcal{L}_{\text{even}}}\Bigg{(}e^{-2\tilde{\beta}}\Bigg{)}^{|c|},$
(212)

where the factor of $2$ is there to account for the fact that each loop
configuration corresponds to two spin configurations. One can similarly
observe that $\mathcal{Z}_{\text{Ising},\ell_{1}}(\tilde{\beta})$ can be
written as,

$\displaystyle\mathcal{Z}_{\text{Ising},\ell_{1}}(\tilde{\beta})$
$\displaystyle=2e^{\tilde{\beta}\,N}\sum_{c\in\mathcal{L}_{\text{even}}}\Bigg{(}e^{-2\tilde{\beta}}\Bigg{)}^{|c+\ell_{1}|}.$
(213)

Therefore, we see that if we define $\tilde{\beta}$ such that,

$\displaystyle e^{-2\tilde{\beta}}=\tanh{\beta},$ (214)

then we have,

$\displaystyle\mathcal{Z}_{0}=\frac{\cosh(\beta)^{N}}{2e^{N\tilde{\beta}}}\Bigg{[}\mathcal{Z}_{\text{Ising}}(\tilde{\beta})+\mathcal{Z}_{\text{Ising},\ell_{1}}(\tilde{\beta})\Bigg{]}$
(215)

As for $\mathcal{Z}_{1}$, note that one can move the domain wall around by
flipping the $\sigma$ spins without changing $\mathcal{Z}_{1}$. Therefore, one
can get rid of the domain wall completely at the expense of flipping all the
fixed spins at the bottom boundary, such that they get fixed at $\sigma=-1$.
Then, one can do the same loop expansion as in Eq. (209), but with the
difference that the loops $c\in\mathcal{L}_{\text{odd}}$ appear with a minus
sign. therefore, we find that,

$\displaystyle\mathcal{Z}_{1}=\frac{\cosh(\beta)^{N}}{2e^{N\tilde{\beta}}}\Bigg{[}\mathcal{Z}_{\text{Ising}}(\tilde{\beta})-\mathcal{Z}_{\text{Ising},\ell_{1}}(\tilde{\beta})\Bigg{]}$
(216)

###  E.4 Deriving Eq. (45) and Eq. (46)

The sums in Eq. (45) and Eq. (46) have two parts,

$\displaystyle\sum_{\bm{m}}(P([e_{\bm{m}}])+(-1)^{j}P([e_{\bm{m}}+\ell_{1}]))^{2}$
$\displaystyle=\sum_{{\bm{m}},k}P([e_{\bm{m}}+\ell_{k}])^{2}+(-1)^{j}2\sum_{{\bm{m}}}P([e_{\bm{m}}])P([e_{\bm{m}}+\ell_{1}]).$
(217)

As for the first part we have,

$\displaystyle\sum_{{\bm{m}},k}P([e_{\bm{m}}+\ell_{k}])^{2}$
$\displaystyle=\sum_{{\bm{m}},k}\sum_{e^{\prime}\in[e_{\bm{m}}+\ell_{k}]}(1-q)^{N-|e^{\prime}|}q^{|e^{\prime}|}P([e_{\bm{m}}+\ell_{k}])$
$\displaystyle=\sum_{{\bm{m}},k}\sum_{e^{\prime}\in[e_{\bm{m}}+\ell_{k}]}(1-q)^{N-|e^{\prime}|}q^{|e^{\prime}|}P([e^{\prime}])$
(218)
$\displaystyle=\sum_{e^{\prime}}(1-q)^{N-|e^{\prime}|}q^{|e^{\prime}|}P([e^{\prime}])=\overline{P([e])}$
$\displaystyle=\frac{[q(1-q)]^{N/2}}{2}\overline{\mathcal{Z}_{\text{RBIM}}(e,J)}.$
(219)

As for the second sum, we find that,

$\displaystyle\sum_{{\bm{m}}}P([e_{\bm{m}}])P([e_{\bm{m}}+\ell_{1}])$
$\displaystyle=\frac{1}{2}\sum_{{\bm{m}},k}P([e_{\bm{m}}+\ell_{k}])P([e_{\bm{m}}+\ell_{k}+\ell_{1}])$
$\displaystyle=\frac{1}{2}\sum_{{\bm{m}},k}\sum_{e^{\prime}\in[e_{\bm{m}}+\ell_{k}]}(1-q)^{N-|e^{\prime}|}q^{|e^{\prime}|}P([e^{\prime}+\ell_{1}])$
(220)
$\displaystyle=\frac{1}{2}\sum_{e^{\prime}}(1-q)^{N-|e^{\prime}|}q^{|e^{\prime}|}P([e^{\prime}+\ell_{1}])=\frac{1}{2}\overline{P([e+\ell_{1}])}$
$\displaystyle=\frac{[q(1-q)]^{N/2}}{4}\overline{\mathcal{Z}_{\text{RBIM}}(e+\ell_{1},J)}.$
(221)

We note that $\mathcal{Z}_{\text{RBIM}}(e+\ell_{1},J)$ can be equivalently
thought of as the random bond Ising model partition function when there is a
domain wall inserted along the support of $\overline{Z}$. So we find that,

$\displaystyle\sum_{\bm{m}}(P([e_{\bm{m}}])+(-1)^{j}P([e_{\bm{m}}+\ell_{1}]))^{2}$
$\displaystyle=\frac{[q(1-q)]^{N/2}}{2}\Bigg{[}\overline{\mathcal{Z}_{\text{RBIM}}(e,J)}+(-1)^{j}\,\overline{\mathcal{Z}_{\text{RBIM}}(e+\ell_{1},J)}\Bigg{]}$
(222)

. We may use Eq. (203) to substitute
$\overline{\mathcal{Z}_{\text{RBIM}}(e,J)}$ with
$\mathcal{Z}_{\text{Ising}}(\tilde{J})$. However, first note that when on the
Nishimori line (see Eq. (193)), one has $\tanh(J)=(1-2q)$. Then, by using Eq.
(195) and Eq. (204), one can find $\tilde{J}$ in terms of $\beta$. After few
steps of algebra, we find that $\tilde{J}=\tilde{\beta}$ which is defined in
Eq. (214). Therefore, by substituting
$\overline{\mathcal{Z}_{\text{RBIM}}(e,J)}$ with
$\mathcal{Z}_{\text{Ising}}(\tilde{\beta})$, and then using Eq. (215) and Eq.
(216) to write the resulting expression in terms of $\mathcal{Z}_{j}$ one
finds,

$\displaystyle\sum_{\bm{m}}(P([e_{\bm{m}}])+(-1)^{j}P([e_{\bm{m}}+\ell_{1}]))^{2}=$
$\displaystyle=\frac{[q(1-q)]^{N/2}}{2}\times\Bigg{[}\frac{\cosh(J)}{\cosh(\tilde{\beta})}\Bigg{]}^{N}\times\frac{2e^{N\tilde{\beta}}}{\cosh(\beta)^{N}}\mathcal{Z}_{j}$
(223) $\displaystyle=e^{-N\beta}\mathcal{Z}_{j}=\mathcal{N}_{j},$ (224)

where from second line to third line we have used Eq. (195), Eq. (193) and Eq.
(214) to write $q$, $J$ and $\tilde{\beta}$ in terms of $\beta$.

##  Appendix F Syndrome correlations in the perturbed toric code

This section contains the derivation of Eq. (51). Let state $\ket{\psi}$
denote a state in the ground state subspace of $H(\beta)$. Since $m_{s}$ and
$m_{s^{\prime}}$ are measurement outcomes of $A_{s}$ and $A_{s^{\prime}}$
operators, the connected correlation function $C$ can be written as,

$\displaystyle
C(s,s^{\prime})=\expectationvalue{A_{s}A_{s^{\prime}}}-\expectationvalue{A_{s}}\expectationvalue{A_{s^{\prime}}},$
(225)

where $\expectationvalue{}$ is the expectation value with respect to
$\ket{\psi}$. As noted in Ref.[14], the ground states of the perturbed toric
code are annihilated by $Q_{s}=A_{s}-\exp(-\beta\sum_{j\in s}Z_{j})$ for any
star $s$ in lattice $\Lambda$, and thus
$A_{s}\ket{\varphi_{0}(\beta)}=\exp(-\beta\sum_{j\in
s}Z_{j})\ket{\varphi_{0}(\beta)}$. Therefore,

$\displaystyle C(s,s^{\prime})=$
$\displaystyle\expectationvalue{e^{-\beta\sum_{j\in s,s^{\prime}}Z_{j}}}$
$\displaystyle-\expectationvalue{e^{-\beta\sum_{j\in
s}Z_{j}}}\expectationvalue{e^{-\beta\sum_{j\in s^{\prime}}Z_{j}}}.$ (226)

To compute the expectation values, we follow similar reasoning as in Appendix
D. First, assume $\ket{\psi}$ is the $\ket{\varphi_{0}(\beta)}$ state. Then,

$\displaystyle\expectationvalue{e^{-\beta\sum_{j\in s}Z_{j}}}_{0}$
$\displaystyle=\frac{1}{\mathcal{Z}_{0}}\expectationvalue{e^{-\beta\sum_{j\in
s}Z_{j}}e^{\beta\sum_{j}Z_{j}}}{\varphi_{0}(0)}$ (227)
$\displaystyle=\overline{\exp(-\beta\sigma_{s}\sum_{\langle\mu,s\rangle}\sigma_{\mu})},$
(228)

where the over line is the average with respect to the 2D Ising model, and
$\mu$ sums over nearest neighbors of $s$. One would find a similar expression
for the case when $\ket{\psi}=\ket{\varphi_{1}(\beta)}$, but the average is
taken with respect to the Ising model with a domain-wall along the support of
$\overline{X}$. More generally, when
$\ket{\psi}=\alpha_{0}\ket{\varphi_{0}(\beta)}+\alpha_{1}\ket{\varphi_{1}(\beta)}$,
one has
$\expectationvalue{A_{s}}=|\alpha_{0}|^{2}\expectationvalue{A_{s}}_{0}+|\alpha_{1}|^{2}\expectationvalue{A_{s}}_{1}$,
since $[A_{s},\overline{Z}]=0$, where the subscripts $1,2$ indicate the
expectation values is taken with respect to $\ket{\varphi_{1,2}(\beta)}$. When
$\beta<\beta_{c}$, the system has a finite correlation length with no long
range order. Moreover, note that we can always choose the support of
$\overline{X}$ to be far from $s$ and $s^{\prime}$, so the domain wall in the
Ising model can be chosen to be $O(L)$ distance apart from $s$ and
$s^{\prime}$. Thus we expect
$\expectationvalue{A_{s}}_{1}=\expectationvalue{A_{s}}_{2}$, up to errors
which are suppressed in the system size. Therefore, if we define
$\epsilon_{s}=-\sigma_{s}\sum_{\langle\mu,s\rangle}\sigma_{\mu}$ to denote the
local energy at site $s$, we find that for any general state in the ground
state subspace of $H(\beta)$ we have,

$\displaystyle
C(s,s^{\prime})=\overline{\exp(\beta\epsilon_{s})\exp(\beta\epsilon_{s^{\prime}})}-\overline{\exp(\beta\epsilon_{s})}\times\overline{\exp(\beta\epsilon_{s^{\prime}})},$
(229)

up to errors which vanish in the thermodynamic limit, where the average is
taken with respect to the 2D classical Ising model.

##  Appendix G Berry connection inside the gapped phase

In this section we argue that under physically motivated assumptions the Berry
connection vanishes in a neighborhood of the fixed point in the gapped phase.
Let $\ket{\varphi_{i}({\bm{\lambda}})}$ denote degenerate ground states at
${\bm{\lambda}}$. Importantly, we assume that
$\ket{\varphi_{i}({\bm{\lambda}})}$ is continuous functions of
${\bm{\lambda}}$ in a neighborhood $\mathcal{A}$ around ${\bm{\lambda}}=0$.
Our goal is to show that
$\matrixelement{\varphi_{j}({\bm{\lambda}})}{\partial_{\bm{\lambda}}}{\varphi_{i}({\bm{\lambda}})}$
vanishes in $\mathcal{A}$ in the thermodynamic limit.

Let
$\delta{\bm{\lambda}}=(\delta\lambda_{1},\delta\lambda_{2},,\cdots,\delta\lambda_{n})$
denote a small change in parameter ${\bm{\lambda}}$. We further assume that
$\ket{\varphi_{i}({\bm{\lambda}}+\delta{\bm{\lambda}})}$ is related to
$\ket{\varphi_{i}({\bm{\lambda}})}$ via a unitary
$u_{\delta{\bm{\lambda}}}=1+i\sum_{k}\delta\lambda_{k}D_{k}+O(\delta\lambda^{2})$,
where $D_{k}$ is a sum of $O(N)$ local Hermitian operators with bounded norm.
This assumption is justified by the fact that the states in a gapped phase are
related to each other by a constant time local unitary evolution. Therefore,

$\displaystyle\ket{\varphi_{i}({\bm{\lambda}}+\delta{\bm{\lambda}})}=$
$\displaystyle\ket{\varphi_{i}({\bm{\lambda}})}+i\sum_{k}\delta\lambda_{k}D_{k}\ket{\varphi_{i}({\bm{\lambda}})}+O(\delta\lambda^{2}).$
(230)

This shows that the components of the Berry connection are given as,

$\displaystyle A_{ij;k}({\bm{\lambda}})$
$\displaystyle=\matrixelement{\varphi_{j}({\bm{\lambda}})}{\partial_{\lambda_{k}}}{\varphi_{i}({\bm{\lambda}})}=i\matrixelement{\varphi_{j}({\bm{\lambda}})}{D_{k}}{\varphi_{i}({\bm{\lambda}})}$
$\displaystyle=d_{k}({\bm{\lambda}})\delta_{ij}+O(\exp(-L)).$ (231)

To arrive at the second line we used the local indistinguishability of the
ground states inside the gapped phase. This shows that the quantum holonomy
$\mathcal{U}_{\gamma}$ in Eq. (59) is just a global phase.

##  Appendix H Repetition code standard decoder and the phase diagram of the
rotated Ising Model

In this section we show that the standard decoder of the repetition code
(described in Section II.1.1) can be used to successfully recover the code
states which were encoded in the ground state subspace of the Ising model
$H_{0}=\sum_{i=1}^{N}Z_{i}Z_{i+1}$ but evolved adiabatically to the ground
state subspace of $H_{\theta}$, with

$\displaystyle H_{\theta}=U_{\theta}H_{0}U^{\dagger}_{\theta},\qquad
U_{\theta}=\prod_{i=1}^{N}e^{-i\frac{\theta}{2}X_{i}},$ (232)

for $-\pi/2<\theta<\pi/2$.

If the initial state has been

$\displaystyle\ket{\psi_{0}}=\alpha\ket{\mathrm{GHZ}_{+}}+\beta\ket{\mathrm{GHZ}_{-}},$
(233)

the adiabatic evolution would transform it into,

$\displaystyle\ket{\psi_{\theta}}=\alpha
U_{\theta}\ket{\mathrm{GHZ}_{+}}+\beta U_{\theta}\ket{\mathrm{GHZ}_{+}},$
(234)

up to a global phase. Imagine one measures all $Z_{i}Z_{i+1}$ stabilizers on
this state, finding the syndrome ${\bm{m}}$, and hence inferring error
$e_{\bm{m}}$ (see Section II.1.1). The probability of seeing syndrom
${\bm{m}}$ is given by,

$\displaystyle\expectationvalue{\Pi_{\bm{m}}}{\psi_{\theta}}=$
$\displaystyle|\alpha|^{2}\expectationvalue{U^{\dagger}_{\theta}\Pi_{\bm{m}}U_{\theta}}{\mathrm{GHZ}_{+}}$
$\displaystyle+|\beta|^{2}\expectationvalue{U^{\dagger}_{\theta}\Pi_{\bm{m}}U_{\theta}}{\mathrm{GHZ}_{-}}.$
(235)

To compute
$\expectationvalue{U^{\dagger}_{\theta}\Pi_{\bm{m}}U_{\theta}}{\mathrm{GHZ}_{\pm}}$,
we note that
$U_{\theta}=\Pi_{j}\big{[}\cos(\frac{\theta}{2})+i\sin(\frac{\theta}{2})X_{j}\big{]}$.
Assuming $|e_{\bm{m}}|=k$, we have,

$\displaystyle\expectationvalue{U^{\dagger}_{\theta}\Pi_{\bm{m}}U_{\theta}}{\mathrm{GHZ}_{\pm}}=$
$\displaystyle\Bigg{|}\cos(\theta/2)^{N-k}(i\sin(\theta/2))^{k}\pm\cos(\theta/2)^{k}(i\sin(\theta/2))^{N-k}\Bigg{|}^{2}$
$\displaystyle=$
$\displaystyle\cos(\theta/2)^{2(N-k)}\sin(\theta/2)^{2k}\times\Big{|}1\pm(i\tan(\theta/2))^{N-2k}\Big{|}^{2},$
(236)

Note that $k/N\leq 1/2$. Let us focus on $k/N<1/2-\delta$ for some positive
constant $\delta>0$, so we may drop the $\tan(\theta/2)$ term inside the
absolute value for large $N$. In that case, the probability that the inferred
error has support on $k$ qubits would simplify to,

$\displaystyle p_{k}={N\choose
k}(\sin(\theta/2)^{2})^{k}(\cos(\theta/2)^{2})^{N-k}$ (237)

which is a binomial distribution with mean,

$\displaystyle k/N=\sin(\theta/2)^{2}$ (238)

. For large $N$, the standard deviation of $k$ grows as $\sqrt{N}$, hence
$k/N$ would be peaked around its mean with a width that goes to zero as
$1/\sqrt{N}$. This justifies the assumption $k/N<1/2-\delta$ (Note that
$-\pi/2<\theta<\pi/2$). Moreover, the above analysis also shows that after
measuring the stabilizers and seeing syndrome ${\bm{m}}$, the state has been
collapsed into,

$\displaystyle\alpha e_{\bm{m}}\ket{\mathrm{GHZ}_{+}}+\beta
e_{\bm{m}}\ket{\mathrm{GHZ}_{-}}+O(\tan(\theta/2)^{N-2k})$ (239)

up to normalization. The error term is what results in the discrepancy between
the recovered state and the original encoded state. However, since the
probability distribution is peaked around $k/N=\sin(\theta/2)^{2}<1/2$, the
typical error term vanishes exponentially in system size. Hence the average
fidelity of the recovered state approaches $1$ in the thermodynamic limit.

##  Appendix I Decoding Adiabatic Noise in a
$\mathbb{Z}_{2}\times\mathbb{Z}_{2}$ SPT State

Consider the Hamiltonian for a one-dimensional cluster state in a transverse
field,

$\displaystyle H_{c}(g)=-\sum_{i=2}^{N-1}Z_{i-1}X_{i}Z_{i+1}-g\sum_{i}X_{i},$
(240)

on a chain of $N$ qubits, where $N$ is assumed to be even. For $b<1$, the
system is in the SPT phase protected by the
$\mathbb{Z}_{2}\times\mathbb{Z}_{2}$ symmetry generated by,

$\displaystyle G_{o}=\prod_{i}X_{2i+1},\quad G_{e}=\prod_{i}X_{2i}.$ (241)

On an open chain the ground state is four-fold degenerate, corresponding to
two logical qubits. At $g=0$, the codespace is stabilized by all terms in the
Hamiltonian, i.e. $Z_{i-1}X_{i}Z_{i+1}$ for $i=2,\cdots,N-1$. The logical
operators at $g=0$ thus could be taken to be,

$\displaystyle\overline{Z}_{1}=Z_{1},\quad\overline{X}_{1}=G_{o}=\prod_{i}X_{2i+1},$
(242)
$\displaystyle\overline{Z}_{2}=Z_{N},\quad\overline{X}_{2}=G_{e}=\prod_{i}X_{2i}.$
(243)

Note that $Z_{1}Z_{N-1}X_{N}$ is a non-trivial symmetric logical operator,
which is equal to $\overline{Z}_{1}\overline{X}_{2}$ up to multiplication by
stabilizers of the code. It means that even if one restricts to symmetric
errors, the code distance is of order $\mathcal{O}(1)$. However, such
symmetric errors are geometrically non-local and they are exponentially
suppressed in any locally symmetric noise model, leading to the possibility of
a finite error-threshold. The standard decoder for the cluster model is
basically two copies of the standard decoder for the repetition code, one for
even lattice sites and the other for the odd lattice sites.

Since $[H_{c}(g),G_{e}]=[H_{c}(g),G_{o}]=0$, $H_{c}(g)$ can be turned into a
block diagonal form, where each block correspond to specific charges under
$G_{e}$ and $G_{o}$. Let $\ket{\varphi_{\pm,\pm}(g)}$, denote the four ground
sates in the four charge sectors, where the first (second) sign denotes the
symmetry charge under $G_{o}$ ($G_{e}$) symmetry. The global phase of
$\ket{\varphi_{\pm,\pm}(g)}$ is fixed by demanding
$\expectationvalue{\partial_{g}}{\varphi_{\pm,\pm}(g)}$ to vanish. Similar to
the non-integrable model studied in Appendix C, we satisfy this condition by
fixing $\innerproduct{\varphi_{\pm,\pm}(0)}{\varphi_{\pm,\pm}(g)}$ to be real.
Fig.12 shows the average fidelity of the recovered state, and the originally
encoded state which we have taken to be
$\ket{\psi_{0}}=\frac{1}{2}\sum_{c_{1}=\pm,c_{2}=\pm}\ket{\varphi_{c_{1},c_{2}}(0)}$.
The inset shows the corresponding data collapse with $\nu=1$ and $g_{c}=1$. As
can be seen from the plot, perfect recovery is possible via the standard
decoder throughout the entire SPT phase.

Figure 12: Fidelity of the recovered state with respect to the original
encoded state (equal supper position of the four ground states of $H_{c}(0)$
with definite symmetries on an open chain) as a function of $g$. The inset
shows the data collapse of the same data points, using $g_{c}=1$ and $\nu=1$.

The success of the standard decoder can be understood b following the same
line of argument as for the repetition code. Note that the SPT phase is
characterized by non-vanishing of string order parameters
$\expectationvalue{s_{i,j}}$ when $|i-j|\gg 1$. Moreover, the string order
parameters $s_{i,j}$ are basically the product of the stabilizers
$Z_{i-1}X_{i}Z_{i+1}$ from $i$ to $j$. This is analogous to
$\expectationvalue{Z_{i}Z_{j}}$ order parameter in the ferromagnetic phase and
the fact that $Z_{i}Z_{j}$ is just the product of the stabilizers of the
repetition code. As such, non-zero expectation value of string order
parameters alongside the indistinguishably of the ground states throughout the
SPT phase ensures that the syndrome data would not reveal any information
about the logical sector and the decoding step would not implement any
erroneous logical operation. We expect the same result to hold for any 1D SPT
phase whenever the symmetry group is finite and Abelian and the SPT phase
corresponds to a maximally non-commutative representation of it, since each
step of the above argument can be generalized in that case.


# Diffusion-based Negative Sampling on Graphs for Link Prediction

Trung-Kien Nguyen∗ University of Southern CaliforniaUnited States
<EMAIL_ADDRESS>and  Yuan Fang  Singapore Management UniversitySingapore
<EMAIL_ADDRESS>
(2024)

###### Abstract.

Link prediction is a fundamental task for graph analysis with important
applications on the Web, such as social network analysis and recommendation
systems, etc. Modern graph link prediction methods often employ a contrastive
approach to learn robust node representations, where negative sampling is
pivotal. Typical negative sampling methods aim to retrieve hard examples based
on either predefined heuristics or automatic adversarial approaches, which
might be inflexible or difficult to control. Furthermore, in the context of
link prediction, most previous methods sample negative nodes from existing
substructures of the graph, missing out on potentially more optimal samples in
the latent space. To address these issues, we investigate a novel strategy of
_multi-level_ negative sampling that enables negative node generation with
flexible and controllable “hardness” levels from the latent space. Our method,
called Conditional Diffusion-based Multi-level Negative Sampling (DMNS),
leverages the Markov chain property of diffusion models to generate negative
nodes in multiple levels of variable hardness and reconcile them for effective
graph link prediction. We further demonstrate that DMNS follows the sub-linear
positivity principle for robust negative sampling. Extensive experiments on
several benchmark datasets demonstrate the effectiveness of DMNS.

link prediction, negative sampling, diffusion models, graph neural networks

∗Work partly done at Singapore Management University.

††journalyear: 2024††copyright: rightsretained††conference: Proceedings of the
ACM Web Conference 2024; May 13–17, 2024; Singapore, Singapore††booktitle:
Proceedings of the ACM Web Conference 2024 (WWW ’24), May 13–17, 2024,
Singapore, Singapore††doi: 10.1145/3589334.3645650††isbn:
979-8-4007-0171-9/24/05††ccs: Computing methodologies Learning latent
representations††ccs: Information systems Data mining

##  1\. Introduction

Graph, which consists of nodes and links between them, is a ubiquitous data
structure for real-world networks and systems. Link prediction (Lü and Zhou,
2011) is a fundamental problem in graph analysis, aiming to model the
probability that two nodes relate to each other in a network or system.
Alternatively, given a query node, link prediction aims to rank other nodes
based on their probability of linking to the query node. Graph link prediction
enables a wide range of applications on the Web, such as friend suggestions in
social networks (Chen et al., 2020), products recommendation in e-commerce
platforms (Ying et al., 2018), and knowledge graph completion (Ahrabian et
al., 2020) for Web-scale knowledge bases.

A prominent approach that has been extensively studied for link prediction is
graph representation learning. It trains an encoder to produce low-dimensional
node embeddings that capture the original graph topology in a latent space.
Toward link prediction, many graph representation learning methods (Tang et
al., 2015; Yang et al., 2020; Hamilton et al., 2017) follow the noise
contrastive estimation approach (Gutmann and Hyvärinen, 2010), which resorts
to sampling a set of positive and negative nodes for any query node.
Specifically, the encoder is trained to capture graph topology by bringing
positive node pairs closer while separating the negative pairs in the
embedding space. While sampling the positive examples for link prediction is
relatively straightforward (e.g., one-hop neighbors of the query node),
sampling the negative examples involves a huge search space that is
potentially quadratic in the number of nodes and a significant fraction of
unlinked node pairs could be false negatives since not all links may be
observed. Therefore, studying negative sampling for contrastive link
prediction on graphs is a crucial research problem.

Many strategies have been proposed for negative sampling on graphs, yet it is
often challenging to flexibly model and control the quality of negative nodes.
While uniform negative sampling (Tang et al., 2015; Hamilton et al., 2017) is
simple, it ignores the quality of negative examples: Harder negative examples
can often contribute more to model training than easier ones. Many studies
explore predefined heuristics to select hard negative nodes, such as
popularity (Mikolov et al., 2013), dynamic selections based on current
predictions (Zhang et al., 2013), $k$-hop neighborhoods (Ahrabian et al.,
2020), Personalized PageRank (Ying et al., 2018), etc. However, heuristics not
only need elaborate designs, but also tend to be inflexible and may not extend
to different kinds of graph. For instance, homophilous and heterophilous
graphs exhibit different connectivity patterns, which means a good heuristic
for one would not work well for the other. Besides heuristics, automatic
methods leveraging generative adversarial networks (GANs) (Cai and Wang, 2017;
Wang et al., 2018) are also popular. They aim to learn the underlying
distribution of the nodes and retrieve harder ones automatically. However, it
is still difficult to flexibly control the “hardness” of the negative examples
for more optimal contrastive learning. Furthermore, it has been shown that the
hardest negative examples may impair the performance (Faghri et al., 2017;
Xuan et al., 2020), when they are nearly indistinguishable from the positive
ones particularly in the early phase of training. To overcome this, we propose
a strategy of _multi-level_ negative sampling, where we can flexibly control
the hardness level of the negative examples according to the need. For
instance, easier negative examples can be used to warm-up model training,
while harder ones are more critical to refining the decision boundary.
Overall, a well-controlled mixture of easy and hard examples are expected to
improve learning.

The idea of multi-level negative sampling immediately brings up the second
question: Where do we find enough negative examples of variable hardness? Most
negative sampling approaches (Tang et al., 2015; Hamilton et al., 2017;
Mikolov et al., 2013; Zhang et al., 2013; Ahrabian et al., 2020; Ying et al.,
2018) are limited to sampling nodes from the observed graph. However, observed
graphs in real-world scenarios are often noisy or incomplete, which are not
ideal sources for directly sampling a sufficient mixture of negative examples
at different levels of hardness. To tackle this issue, we propose to
synthesize more negative nodes to complement the real ones in the latent
space. The latent space can potentially provide infinite negative examples
with arbitrary hardness, to aid the generation of multi-level negative
examples in a flexible and controllable manner. Although a few studies (Hu et
al., 2019; Li et al., 2023) also utilize GANs to generate additional samples
in the latent space, they are not designed to synthesize multi-level examples.

To materialize our multi-level strategy, a natural choice is diffusion models.
Recently, diffusion models have achieved promising results in generation tasks
(Ho et al., 2020; Song and Ermon, 2019), particularly in visual applications
(Ho et al., 2022b; Dhariwal and Nichol, 2021). While GAN-based models are
successful in generating high-fidelity examples, they face many issues in
training such as vanishing gradients and mode collapse (Arjovsky and Bottou,
2017; Arjovsky et al., 2017). In contrast, diffusion models follow a
reconstruction mechanism that offers a stable training process (Luo, 2022).
More specific to our context, a desirable property of diffusion model is that
it utilizes a Markov chain with multiple time steps to denoise random input,
where the sample generated at each time step is conditioned on the sample in
the immediately preceding step. Hence, during the generation process, we can
naturally access the generated samples at different denoised time steps to
achieve our multi-level strategy. Hence, we propose a novel diffusion-based
framework to generate negative examples for graph link prediction, named
Conditional Diffusion-based Multi-level Negative Sampling (DMNS). On one hand,
we employ a diffusion model to learn the one-hop connectivity distribution of
any given query node, i.e., the positive distribution in link prediction. The
diffusion model allows us to flexibly control the hardness of each negative
example by looking up a specified time step $t$, ranging from a virtually
positive example (or indistinguishable from the positives) when $t=0$ and
harder examples as $t\to 0$, to easier ones amounting to random noises as
$t\to\infty$. On the other hand, we adopt a conditional diffusion model (Ho et
al., 2022b; Dhariwal and Nichol, 2021), which is designed to explicitly
condition the generation on side information. In graph link prediction, we
condition the positive distribution of a specific query node on the node
itself, and thus obtain query-specific diffusion models. Theoretically, we
show that the density function of our negative examples obeys the sub-linear
positivity principle (Yang et al., 2020) under some constraint, ensuring
robust negative sampling on graphs.

Our contributions in this work are summarized as follows. 1) We investigate
the strategy of multi-level negative sampling on graphs for contrastive link
prediction. 2) We propose a novel framework DMNS based on a conditional
diffusion model, which generates multi-level negative examples that can be
flexibly controlled to improve training. 3) We show that the distribution
function of our negative examples follows the sub-linear positivity principle
under a defined constraint. 4) We conduct extensive experiments on several
benchmark datasets, showing that our model outperforms state-of-the-art
baselines on graph link prediction.

##  2\. Related Work

In this section, we briefly review and discuss related work from three
perspectives, namely, link prediction, negative sampling, and diffusion
models.

#### Link prediction

Classical methods on link prediction hinge on feature engineering to measure
the connection probability between two nodes, such as based on the proximity
(Page et al., 1999), or the similarity of neighborhoods (Adamic and Adar,
2003; Zhou et al., 2009). The success of deep learning has motivated extensive
studies on graph representation learning (Kipf and Welling, 2017; Veličković
et al., 2018; Hamilton et al., 2017). On learned graph representations, a
simple strategy for link prediction is to employ a node-pair decoder (Grover
and Leskovec, 2016). More sophisticated approach exploits the representations
of the enclosing subgraphs. For example, SEAL (Zhang and Chen, 2018) proposes
the usage of local subgraphs based on the $\gamma$-decaying heuristic theory.
Meanwhile, ScaLED (Louis et al., 2022) utilizes random walks to efficiently
sample subgraphs for large-scale networks. On the other hand, BUDDY
(Chamberlain et al., 2022) proposes subgraph sketches to approximate essential
features of the subgraphs without constructing explicit ones.

#### Negative sampling

Graph representation learning approaches for link prediction commonly employ
contrastive strategies (Ying et al., 2018; You et al., 2020). The contrastive
methods require effective negative sampling to train the graph encoders.
Various negative sampling heuristics have been proposed, such as based on the
popularity of examples (Mikolov et al., 2013), current prediction scores
(Zhang et al., 2013) and selecting high-variance samples (Ding et al., 2020).
On graphs, SANS (Ahrabian et al., 2020) select negative examples from $k$-hop
neighborhoods. MixGCF (Huang et al., 2021) synthesizes hard negative examples
by hop and positive mixing. MCNS (Yang et al., 2020) develops a Metropolis-
Hasting sampling approach based on the proposed sub-linear positivity theory.
Another line of research utilizes generative adversarial networks (GANs) for
automatic negative sampling. Among them, GraphGAN (Wang et al., 2018) and
KBGAN (Cai and Wang, 2017) learns a distribution over negative candidates,
while others generate new examples not found in the original graph (Liu et
al., 2021a; Hu et al., 2019; Li et al., 2023). However, these methods cannot
easily control the hardness of the negative examples, which is the key
motivation of our multi-level strategy.

#### Diffusion models

Diffusion models (Ho et al., 2020; Sohl-Dickstein et al., 2015; Song and
Ermon, 2019; Song et al., 2020) have become state-of-the-art generative
models. They gradually inject noises into the data and then learn to reverse
this process for sampling. Denoising diffusion probabilistic models (DDPMs)
(Ho et al., 2020; Sohl-Dickstein et al., 2015) aims to predict noises from the
diffused output at arbitrary time steps. Score-based Generative Models (SGMs)
(Song and Ermon, 2019, 2020; Song et al., 2020) aims to predict the score
function $\nabla\log(p_{x_{t}})$. They have different approaches but are shown
to be equivalent to optimizing a diffusion model. Conditional diffusion models
(Ho et al., 2022b; Dhariwal and Nichol, 2021) permit explicit control of
generated samples via additional conditions, which enables a wide range of
applications such as visual generation (Saharia et al., 2022), NLP (Gong et
al., 2022), multi-modal generation (Nichol et al., 2021; Ho et al., 2022a),
etc. Recently, some studies have adopted diffusion models in graph generation
tasks. EDM (Hoogeboom et al., 2022) learns a diffusion process to work on 3d
molecule generation, while GDSS (Jo et al., 2022) learns to generate both node
features and adjacency matrices. However, they do not aim to generate multi-
level node samples for graph link prediction. We refer readers to a
comprehensive survey (Yang et al., 2022).

##  3\. Preliminaries

In this section, we briefly introduce the background of graph neural networks
and diffusion models, which are the foundation of our proposed DMNS.

#### Graph neural networks

Message-passing GNNs usually resort to multi-layer neighborhood aggregation,
in which each node recursively aggregates information from its neighbors.
Specifically, the representation of node $v$ in the $l$-th layer,
$\mathbf{h}^{l}_{v}\in\mathbb{R}^{d_{h}^{l}}$, can be constructed as

(1)
$\displaystyle\mathbf{h}^{l}_{v}=\sigma\left(\textsc{Aggr}(\mathbf{h}^{l-1}_{v},\\{\mathbf{h}^{l-1}_{i}:i\in\mathcal{N}_{v}\\};\omega^{l})\right),$

where $d_{h}^{l}$ is the dimension of node representations in the $l$-th
layer, $\textsc{Aggr}(\cdot)$ denotes an aggregation function such as mean-
pooling (Kipf and Welling, 2017), self-attention (Veličković et al., 2018) or
concatenation (Hamilton et al., 2017), $\sigma$ is an activation function,
$\mathcal{N}_{v}$ denotes the set of neighbors of $v$, and $\omega^{l}$
denotes the learnable parameters in layer $l$.

#### Diffusion models

Denoising diffusion model (DDPM) (Ho et al., 2020) boils down to learning the
Gaussian transitions of Markov chains. Specifically, DDPM consists of two
Markov chains: a forward diffusion process, and a backward denoising process
(also known as the “reverse” process). The forward chain transforms input data
to complete noise by gradually adding Gaussian noise at each step:

(2) $\displaystyle q(\mathbf{x}_{1:T}|\mathbf{x}_{0})$
$\displaystyle=\textstyle\prod^{T}_{t=1}q(\mathbf{x}_{t}|\mathbf{x}_{t-1}),$
(3) $\displaystyle q(\mathbf{x}_{t}|\mathbf{x}_{t-1})$
$\displaystyle=\mathcal{N}(\mathbf{x}_{t};\sqrt{1-\beta_{t}}\mathbf{x}_{t-1},\beta_{t}\mathbf{I}),$

where $T$ is the total number of time steps of the diffusion process, and
$\beta_{t}$ denotes the variances in the diffusion process, which can be
learnable or fixed constants via some scheduling strategy. By a
reparameterization trick, the closed form of output $\mathbf{x}_{t}$ at
arbitrary time step $t$ can be obtained as

(4)
$\displaystyle\mathbf{x}_{t}=\sqrt{\bar{\alpha}_{t}}\mathbf{x}_{0}+\sqrt{1-\bar{\alpha}_{t}}\epsilon_{t},$

where $\alpha_{t}=1-\beta_{t}$, $\bar{\alpha}_{t}=\prod^{t}_{i=1}\alpha_{i}$,
and $\epsilon_{t}\sim\mathcal{N}(\mathbf{0},\mathbf{I})$. The backward
denoising process learns to reconstruct the input from noise:

(5) $\displaystyle p_{\theta}(\mathbf{x}_{0:T})$
$\displaystyle=p(\mathbf{x}_{T})\textstyle\prod^{T}_{t=1}p_{\theta}(\mathbf{x}_{t-1}|\mathbf{x}_{t})$
(6) $\displaystyle p_{\theta}(\mathbf{x}_{t-1}|\mathbf{x}_{t})$
$\displaystyle=\mathcal{N}(\mathbf{x}_{t-1};\mu_{\theta}(\mathbf{x}_{t},t),\Sigma_{\theta}(\mathbf{x}_{t},t))$
(7) $\displaystyle p(\mathbf{x}_{T})$
$\displaystyle=\mathcal{N}(\mathbf{x}_{T};\mathbf{0},\mathbf{I}),$

where $\theta$ parameterizes the diffusion model. To model the joint
distribution $p_{\theta}(\mathbf{x}_{t-1}|\mathbf{x}_{t})$, DDPM sets
$\Sigma_{\theta}(\mathbf{x}_{t},t)=\beta_{t}\mathbf{I}$ as scheduled
constants, and derive $\mu_{\theta}(\mathbf{x}_{t},t)$ by a reparameterization
trick as

(8)
$\displaystyle\mu_{\theta}(\mathbf{x}_{t},t)=\textstyle\frac{1}{\sqrt{\alpha_{t}}}\mathbf{x}_{t}-\frac{1-\alpha_{t}}{\sqrt{\alpha_{t}}\sqrt{1-\bar{\alpha}_{t}}}\epsilon_{\theta}(\mathbf{x}_{t},t),$

where $\epsilon_{\theta}(\mathbf{x}_{t},t)$ is a function to estimate the
source noise $\epsilon$ and can be implemented as a neural network. The
objective now learns to predict added noises instead of means, which has been
shown to achieve better performance.

##  4\. Proposed Model: DMNS

In this section, we introduce the proposed method DMNS to generate multi-level
negative examples for graph link prediction. Before we begin, we first present
the overall framework of DMNS in Fig. 1. We employ a standard GNN encoder to
obtain node embeddings, which captures content and structural neighborhood
information. Next, we train a diffusion model to learn the neighborhood
distribution conditioned on the query node. From the model we sample several
output embeddings at different time steps, to serve as negative examples at
multi-level hardness for contrastive learning.

Figure 1. Overall framework of DMNS.

###  4.1. Multi-level Negative Sample Generation

Given a query node, we aim to learn the distribution of its 1-hop neighbors,
i.e., its positive distribution for the link prediction task. Vanilla
diffusion model only generates generic examples without the ability to
personalize for the query node. To this end, we leverage the conditional
diffusion model (Ho et al., 2022b; Dhariwal and Nichol, 2021), taking query
node embeddings as additional information for sample generation. From that, we
can perform multi-level negative sampling by extracting generated node
embeddings from multiple time steps. The choice of time steps empowers us to
control the hardness of negative examples. For instance, the embedding output
from the final step of the denoise (reverse) process ($t=0$) is the hardest to
distinguish from the positive nodes in the latent space, while its
counterparts from earlier steps (larger $t$’s) can be taken as progressively
easier negative nodes. In general, the time step $t$ is a proxy to the
hardness of negative examples: a smaller $t$ gives harder examples. Thus, we
can automatically incorporate such multi-level negative examples into training
the link prediction task.

Specifically, for a query node $v$, we obtain its embedding
$\mathbf{h}_{v}\in\mathbb{R}^{d_{h}}$ from a GNN encoder, as well as its
neighbors’ embeddings $\mathbf{h}_{u},\ \forall u\in\mathcal{N}_{v}$ where
$\mathcal{N}_{v}$ denotes the neighbor set of $v$. The node embeddings will be
fed into the diffusion model to learn the neighbor (positive) distribution of
the query node. In the forward diffusion process, we gradually add bite-sized
noise to the neighbor node $u$’s embedding. Following the reparameterization
trick (Ho et al., 2020; Luo, 2022), we can obtain the closed form of the
output $\mathbf{h}_{u,t}$ at an arbitrary time step $t$ without relying on the
intermediate output:

(9)
$\displaystyle\mathbf{h}_{u,t}=\sqrt{\bar{\alpha}_{t}}\mathbf{h}_{u}+\sqrt{1-\bar{\alpha}_{t}}\epsilon_{t},\quad\forall
u\in\mathcal{N}_{v},$

In the denoise process, we aim to predict the added noise from the diffused
node $u$ at time step $t$ given the query node $v$, which is formulated as

(10)
$\displaystyle\epsilon_{t,\theta|v}=\tau(\mathbf{h}_{u,t},\mathbf{t},\mathbf{h}_{v};\theta),\quad\forall
u\in\mathcal{N}_{v},$

where $\mathbf{t}\in\mathbb{R}^{d_{h}}$ is a continuous embedding vector for
time $t$, and $\theta$ denotes learnable parameters. Following earlier work
(Ho et al., 2020), we employ the sinusoidal positional encoding for the time
steps, such that $[\mathbf{t}]_{2i}=\sin(t/10000^{\frac{2i}{d_{h}}})$ and
$[\mathbf{t}]_{2i+1}=\cos(t/10000^{\frac{2i}{d_{h}}})$, followed by a
multilayer perceptron (MLP). $\tau(\cdot;\theta)$ is a learnable
transformation function, which estimates the transformation of diffused
embeddings at time $t$ to predict the noise. In our model, we implement $\tau$
as a feature-wise linear modulation (FiLM) (Perez et al., 2018) layer, which
is conditioned on both the time step $t$ and query node embedding
$\mathbf{h}_{v}$. Specifically,

(11)
$\displaystyle\epsilon_{t,\theta|v}=(\gamma+\mathbf{1})\odot\mathbf{h}_{u,t}+\eta,$

where $\gamma$ and $\eta\in\mathbb{R}^{d_{h}}$ are scaling and shifting
vectors to transform the diffused node embedding, respectively; $\mathbf{1}$
denotes a vector of ones to center the scaling around one, and $\odot$ denotes
element-wise multiplication. The transformation vectors are learnable,
conditioned on the diffusion time embedding $\mathbf{t}$ and the query node
embedding $\mathbf{h}_{v}$. We implement $\gamma$ and $\beta$ as fully
connected layers (FCLs), dependent upon the conditional embeddings as follows.

(12)
$\displaystyle\gamma=\textsc{FCL}(\mathbf{t}+\mathbf{h}_{v};\theta_{\gamma}),\quad\eta=\textsc{FCL}(\mathbf{t}+\mathbf{h}_{v};\theta_{\eta}),$

where $\theta_{\gamma},\theta_{\eta}$ are parameters of the FCLs. We aggregate
the time step and query node information by summation, yet other methods such
as a neural network can also be considered. Also note that in practical
implementations, multiple FiLM layers can be stacked to enhance model
capacity.

###  4.2. Training Objective

The diffusion model and GNN encoder are trained simultaneously in an
alternating manner. In the outer loop, the GNN encoder is updated for link
prediction, taking into account the multi-level negative nodes generated by
the diffusion model conditioned on the query node. In the inner loop, the
diffusion model is updated to generate the positive neighbor distribution of
the query node, based on the current GNN encoder. More details of the training
process are outlined in Appendix B. In the following, we discuss the loss
functions for the diffusion and GNN-based link prediction.

#### Diffusion loss

We employ a mean squared error (MSE) between sampled noises from the forward
process and predicted noises from the reverse process (Ho et al., 2020). That
is, at time step $t$,

(13) $\mathcal{L}_{D}=\|\epsilon_{t}-\epsilon_{t,\theta|v}\|^{2}$

After the diffusion model is updated, we utilize it to synthesize node
embeddings for the main link prediction task. Starting from complete noises
($t=T$), we obtain generated samples by sequentially removing predicted noises
at each time step.

(14) $\displaystyle\mathbf{h}_{d,T|v}$
$\displaystyle\sim\mathcal{N}(\mathbf{0},\mathbf{I}),$ (15)
$\displaystyle\quad\mathbf{h}_{d,t-1|v}$
$\displaystyle=\textstyle\frac{1}{\sqrt{\alpha_{t}}}\left(\mathbf{h}_{d,t|v}-\frac{1-\alpha_{t}}{\sqrt{1-\bar{\alpha}_{t}}}\epsilon_{t,\theta|v}\right)+\sigma_{t}\mathbf{z},$

where the standard deviation $\sigma_{t}=\sqrt{\beta_{t}}$ and
$\mathbf{z}\sim\mathcal{N}(\mathbf{0},\mathbf{I})$.

The sequence of generated samples $\\{\mathbf{h}_{d,t|v}:0\leq t\leq T\\}$
represents multi-level negative nodes in the latent space, w.r.t. a query node
$v$. We can easily control the number and hardness requirement of the negative
sampling, by choosing certain time steps between $0$ and $T$. On the one hand,
sampling from too many smaller time steps consumes more memory but brings
little diversity. On the other hand, sampling from larger time steps may bring
in trivial noises which are easier to distinguish. In our implementation, we
balance the multi-level strategy by choosing the output from a range of well-
spaced steps to form our generated negative set for the query node $v$:
$D_{v}=\\{\mathbf{h}_{d,t|v}:t=\frac{T}{10},\frac{T}{8},\frac{T}{4},\frac{T}{2}\\}$,
which are both efficient and robust.

#### Link prediction loss

We adopt the log-sigmoid loss (Hamilton et al., 2017) for the link prediction
task. Consider a quadruplet $(v,u,u^{\prime},D_{v})$, where $v$ is a query
node, $u$ is a positive node linked to $v$, $u^{\prime}$ is an existing
negative node randomly sampled from the graph, and $D_{v}$ is a set of multi-
level negative nodes w.r.t. the query $v$, which are latent node embeddings
generated from the diffusion model at chosen time steps. Note that we still
employ existing nodes from the graph as negative nodes, to complement the
samples generated in the latent space. Then, the loss is formulated as

$\displaystyle\mathcal{L}=$
$\displaystyle-\log\sigma(\mathbf{h}^{\top}_{v}\mathbf{h}_{u})-\log\sigma(-\mathbf{h}^{\top}_{v}\mathbf{h}_{u^{\prime}})$
(16) $\displaystyle-\textstyle\sum_{d_{i}\in
D_{v}}w_{i}\log\sigma(-\mathbf{h}^{\top}_{v}\mathbf{h}_{d_{i}}))$

where $\sigma(\cdot)$ is the sigmoid activation, and $w_{i}$ is the weight
parameter for each negative example generated at different time steps. The
idea is, given different levels of hardness associated with different time
steps, the importance of the examples from different steps also varies. While
various strategies for $w$ can be applied, we use a simple linearly decayed
sequence of weights with the assumption that closer steps (harder examples)
are more important. Further investigation on the choice and weighting of
samples will be discussed in Sect. 5.

###  4.3. Theoretical Analysis

We present a theoretical analysis to justify the samples generated from our
conditional diffusion model. Specifically, a Sub-linear Positivity Principle
(Yang et al., 2020) has been established earlier for robust negative sampling
on graph data. The principle states that the negative distribution should be
_positively_ but _sub-linearly correlated_ to the positive distribution, which
has been shown to be able to balance the trade-off between the embedding
objective and expected risk. Here, we show that the negative examples
generated from the diffusion model are in fact drawn from a negative
distribution that follows the principle.

###### Theorem 1 (Sub-linear Positivity Diffusion).

Consider a query node $v$. Let
$\mathbf{x}_{n}\sim\mathcal{N}(\mu_{t,\theta},\Sigma_{t,\theta})$ and
$\mathbf{x}_{p}\sim\mathcal{N}(\mu_{0,\theta},\Sigma_{0,\theta})$ represent
samples drawn from the negative and positive distributions of node $v$,
respectively. Suppose the parameters of the two distributions are specified by
a diffusion model $\theta$ conditioned on the query node $v$ at time $t>0$ and
$0$, respectively. Then, the density function of the negative samples $f_{n}$
is sub-linearly correlated to that of the positive samples $f_{p}$:

(17) $\displaystyle f_{n}(\mathbf{x}_{n}|v)\propto
f_{p}(\mathbf{x}_{p}|v)^{\lambda},\quad\text{ for some }0<\lambda<1,$

as long as $\Psi\geq 0$, which is a random variable given by
$\Psi=2\Delta^{\top}\sqrt{\bar{\alpha}_{t}}(\mathbf{x}_{0}-\mu_{0})+\Delta^{\top}\Delta\geq
0$, where
$\Delta=\sqrt{\bar{\alpha}_{t}}\mu_{0}+\sqrt{1-\bar{\alpha}_{t}}\epsilon_{0}-\mu_{t}$,
$\mathbf{x}_{0}$ is generated by the model $\theta$ at time 0, and
$\epsilon_{0}\sim\mathcal{N}(\mathbf{0},\mathbf{I})$. $\Box$

###### Proof.

Note that conditional diffusion $\theta$ aims to learn neighbor (positive)
distribution of query node $v$, the samples generated from time step 0 can be
regarded as positive samples while their counterparts from larger $t>0$
treated as negative ones. Then the density functions of the generated positive
and negative sampling distributions are as follows:

(18) $\displaystyle
f_{p}(\mathbf{x}_{p}|v)=\mathcal{N}(\mathbf{x}_{0};\mu_{0},\Sigma^{2}_{1})$
$\displaystyle=\frac{1}{2\pi^{k/2}\det(\Sigma_{1})^{1/2}}\exp\biggl{\\{}-\frac{1}{2}(\mathbf{x}_{0}-\mu_{0})^{\top}\Sigma^{-1}_{1}(\mathbf{x}_{0}-\mu_{0})\biggl{\\}}$
$\displaystyle=\frac{1}{2\pi^{k/2}\beta_{1}^{k/2}}\exp\biggl{\\{}-\frac{1}{2\beta_{1}}(\mathbf{x}_{0}-\mu_{0})^{\top}(\mathbf{x}_{0}-\mu_{0})\biggl{\\}}$

(19) $\displaystyle
f_{n}(\mathbf{x}_{n}|v)=\mathcal{N}(\mathbf{x}_{t};\mu_{t},\Sigma^{2}_{t+1})$
$\displaystyle=\frac{1}{2\pi^{k/2}\det(\Sigma_{t+1})^{1/2}}\exp\biggl{\\{}-\frac{1}{2}(\mathbf{x}_{t}-\mu_{t})^{\top}\Sigma^{-1}_{t+1}(\mathbf{x}_{t}-\mu_{t})\biggl{\\}}$
$\displaystyle=\frac{1}{2\pi^{k/2}\beta_{t+1}^{k/2}}\exp\biggl{\\{}-\frac{1}{2\beta_{t+1}}(\mathbf{x}_{t}-\mu_{t})^{\top}(\mathbf{x}_{t}-\mu_{t})\biggl{\\}}$

where
$\mu_{0}=\mu_{\theta}(\mathbf{x}_{1},1,\mathbf{v}),\mu_{t}=\mu_{\theta}(\mathbf{x}_{t+1},t+1,\mathbf{v})$,
$\Sigma_{i}=\beta_{i}\mathbf{I}$, $k$ is the vector dimension.

By applying reparameterization trick (Ho et al., 2020; Luo, 2022) we derive:

(20) $\displaystyle\mathbf{x}_{t}$
$\displaystyle=\sqrt{\bar{\alpha}_{t}}\mathbf{x}_{0}+\sqrt{1-\bar{\alpha}_{t}}\epsilon_{0}$

where $\epsilon_{0}\sim\mathcal{N}(\mathbf{0},\mathbf{I})$. Replacing $x_{t}$
by $x_{0}$ in $f_{n}(\mathbf{x}_{n}|v)$ we obtain:

(21) $\displaystyle
f_{n}(\mathbf{x}_{n}|v)=\frac{1}{2\pi^{k/2}\beta_{t+1}^{k/2}}\exp\biggl{\\{}-\frac{1}{2\beta_{t+1}}(\sqrt{\bar{\alpha}_{t}}\mathbf{x}_{0}+\sqrt{1-\bar{\alpha}_{t}}\epsilon_{0}-\mu_{t})^{\top}$
$\displaystyle(\sqrt{\bar{\alpha}_{t}}\mathbf{x}_{0}+\sqrt{1-\bar{\alpha}_{t}}\epsilon_{0}-\mu_{t})\biggl{\\}}$
(22)
$\displaystyle=\frac{1}{2\pi^{k/2}\beta_{t+1}^{k/2}}\exp\biggl{\\{}-\frac{1}{2\beta_{t+1}}\Big{[}\sqrt{\bar{\alpha}_{t}}(\mathbf{x}_{0}-\mu_{0})+\Delta\Big{]}^{\top}$
$\displaystyle\Big{[}\sqrt{\bar{\alpha}_{t}}(\mathbf{x}_{0}-\mu_{0})+\Delta\Big{]}\biggl{\\}}$
(23)
$\displaystyle=\frac{1}{2\pi^{k/2}\beta_{t+1}^{k/2}}\exp\biggl{\\{}-\frac{1}{2\beta_{t+1}}\Big{[}\bar{\alpha}_{t}(\mathbf{x}_{0}-\mu_{0})^{\top}(\mathbf{x}_{0}-\mu_{0})$
$\displaystyle+2\Delta^{\top}\sqrt{\bar{\alpha}_{t}}(\mathbf{x}_{0}-\mu_{0})+\Delta^{\top}\Delta\Big{]}\biggl{\\}}$

with
$\Delta=\sqrt{\bar{\alpha}_{t}}\mu_{0}+\sqrt{1-\bar{\alpha}_{t}}\epsilon_{0}-\mu_{t}$.
We denote
$\Psi=2\Delta^{\top}\sqrt{\bar{\alpha}_{t}}(\mathbf{x}_{0}-\mu_{0})+\Delta^{\top}\Delta$.
If $\Psi\geq 0$, then:

(24)
$\bar{\alpha}_{t}(\mathbf{x}_{0}-\mu_{0})^{\top}(\mathbf{x}_{0}-\mu_{0})+\Psi\geq\bar{\alpha}_{t}(\mathbf{x}_{0}-\mu_{0})^{\top}(\mathbf{x}_{0}-\mu_{0})$
(25)
$\equiv-\frac{1}{2\beta_{t+1}}\Big{[}\bar{\alpha}_{t}(\mathbf{x}_{0}-\mu_{0})^{\top}(\mathbf{x}_{0}-\mu_{0})+\Psi\Big{]}$
(26)
$\leq-\frac{1}{2\beta_{t+1}}\bar{\alpha}_{t}(\mathbf{x}_{0}-\mu_{0})^{\top}(\mathbf{x}_{0}-\mu_{0})$
(27) $\displaystyle\equiv
f_{n}(\mathbf{x}_{n}|v)\leq\frac{1}{2\pi^{k/2}\beta_{t+1}^{k/2}}\exp\biggl{\\{}-\frac{1}{2\beta_{t+1}}\bar{\alpha}_{t}(\mathbf{x}_{0}-\mu_{0})^{\top}(\mathbf{x}_{0}-\mu_{0})\biggl{\\}}$
(28)
$\displaystyle\leq\frac{\beta_{1}^{k/2}}{\beta_{t+1}^{k/2}}\frac{1}{2\pi^{k/2}\beta_{1}^{k/2}}\exp\biggl{\\{}-\frac{\beta_{1}}{2\beta_{1}\beta_{t+1}}\bar{\alpha}_{t}(\mathbf{x}_{0}-\mu_{0})^{\top}(\mathbf{x}_{0}-\mu_{0})\biggl{\\}}$
(29)
$\displaystyle\leq\frac{\beta_{1}^{k/2}}{\beta_{t+1}^{k/2}}\frac{1}{2\pi^{k/2}\beta_{1}^{k/2}}\exp\biggl{\\{}-\frac{1}{2\beta_{1}}(\mathbf{x}_{0}-\mu_{0})^{\top}(\mathbf{x}_{0}-\mu_{0})\biggl{\\}}^{\lambda}$
(30) $\displaystyle\propto f_{p}(\mathbf{x}_{p}|v)^{\lambda}$

where
$0<\lambda=\big{(}\frac{\beta_{1}}{\beta_{t+1}}\big{)}\bar{\alpha}_{t}<1$
$\big{(}\beta_{1}<\beta_{t+1}$ through variances scheduling and
$0<\bar{\alpha}_{t}<1\big{)}$. Therefore, the density function of negative
samples is sub-linearly correlated to that of positive samples under the
constraint
$\Psi=2\Delta^{\top}\sqrt{\bar{\alpha}_{t}}(x_{0}-\mu_{0})+\Delta^{\top}\Delta\geq
0$. ∎

Specifically, we run an experiment to verify the distribution of $\Psi$. We
use our diffusion model to generate a large number of examples at time $0$ and
a given $t$, to compute their mean $\mu_{0}$ and $\mu_{t}$, respectively. We
present the empirical distributions (histograms) of $\Psi$ on the four
datasets, across time steps
${\frac{T}{10},\frac{T}{8},\frac{T}{4},\frac{T}{2}}$ in Fig. 2. The
probabilities that $\Psi\geq 0$ averaged over samples from 4 time steps on
Cora, Citeseer, Coauthor-CS and Actor are 80.14%, 81.62%, 81.65% and 84.31%,
respectively. The results indicate that the majority of generated examples
from DMNS adhere to the sub-linear positivity theorem in practice.

Figure 2. Empirical distributions (histograms) of $\Psi$ on (a1–a4) Cora,
(b1–b4) Citeseer, (c1–c4) Coauthor-CS, (d1–d4) Actor, across different time
steps.

###  4.4. Algorithm and Complexity

Algorithm 1 Model Training for DMNS

1:Graph $G=(V,E)$, training triplets $T=(v,u,u^{\prime})$

2:GNN model parameters $\omega$, Diffusion model parameters $\theta$.

3:initialize parameters $\omega$, $\theta$

4:while not converged do

5: sample a batch of triplets $T_{\text{b}}\subset T$;

6: for each node $v$ in the batch $T_{\text{b}}$ do

7:
$\mathbf{h}_{v}=\sigma\left(\textsc{Aggr}(\mathbf{h}_{v},\\{\mathbf{h}_{i}:i\in\mathcal{N}_{v}\\};\omega)\right)$;

8: //Train diffusion

9: while not converged do

10: for each query node $v$ in the batch $T_{\text{b}}$ do

11: $u\sim N_{v}$, $\epsilon\sim\mathcal{N}(\mathbf{0},\mathbf{I})$;

12: $t\sim Uniform({1,..,T})$, $\mathbf{t}=MLP(t)$;

13:
$\epsilon_{t,\theta|v}=\tau(\mathbf{h}_{u,t},\mathbf{t},\mathbf{h}_{v};\theta)$;

14: $\mathcal{L}_{D}=\|\epsilon_{t}-\epsilon_{t,\theta|v}\|^{2}$;

15: update $\theta$ by minimizing $\mathcal{L}_{D}$ with $\omega$ fixed;

16: //Sampling multi-level generated negative nodes

17: for each query node $v$ in the batch $T_{\text{b}}$ do

18: $\mathbf{h}_{d,T|v}\sim\mathcal{N}(\mathbf{0},\mathbf{I})$,

19: for $t$ = T-1,.., 0 do

20: $\mathbf{z}\sim\mathcal{N}(\mathbf{0},\mathbf{I})$ if $t>1$ else
$\mathbf{z}=0$ ;

21:
$\mathbf{h}_{d,t-1|v}=\textstyle\frac{1}{\sqrt{\alpha_{t}}}\left(\mathbf{h}_{d,t|v}-\frac{1-\alpha_{t}}{\sqrt{1-\bar{\alpha}_{t}}}\epsilon_{t,\theta|v}\right)+\sigma_{t}\mathbf{z}$;

22:
$D_{v}=\\{\mathbf{h}_{d,t|v}:t=\frac{T}{10},\frac{T}{8},\frac{T}{4},\frac{T}{2}\\}$

23: Calculate $\mathcal{L}$ as Eqn. (15);

24: update $\omega$ by minimizing $\mathcal{L}$ with $\theta$ fixed;

25:return $\omega$, $\theta$.

We outline the model training for DMNS in Algorithm 1. In line 1, we
initialize the model parameters. In line 3, we sample a batch of triplets from
training data. In lines 4–6, we obtain embeddings for all nodes in training
batch by the GNN encoder. In lines 8–16, we train conditional diffusion model
to learn the neighborhood distribution of given query node $v$. Specifically,
we calculate the predicted noise at arbitrary time step $t$ conditioned on
query node in line 12. We compute the diffusion loss and update diffusion
parameters in lines 14–15. In lines 18–25, we sample a set multi-level
negative nodes for query node by diffusion model. In lines 26-27, we compute
the main link prediction loss and update the parameters of GNN encoder.

We analyze the complexity of DMNS for one node. Taking GCN as base encoder,
the neighborhood aggregation for one node in the $l$-th layer has complexity
$O(d_{l}d_{l-1}\bar{d})$, where $d_{l}$ is the dimension of the $l$-th layer
and $\bar{d}$ is the average node degree. The computation for diffusion model
includes time embeddings module and noise estimation module. The time
embeddings module employs a MLP of $L_{1}$ layers, where a $l$-layer has
complexity $O(d_{l}d_{l-1})$, with $d_{l}$ is the layer dimension. The noise
estimation module employs a neural network of $L_{2}$ FiLM layers, where
$l$-th layer has complexity of $O(2d_{l}d_{l-1})$ with $d_{l}$ is the layer
dimension. Thus, the diffusion training takes $N$ iterations has complexity
$O(NL_{1}L_{2}d_{l}d_{l-1})$. After that, the sampling runs the reverse
process of $T$ steps (lines 20–23) has complexity $O(2Td_{l}d_{l-1})$.
Compared to the base GCN, the overhead of diffusion part has complexity
$O\big{(}(2T+NL_{1}L_{2})d_{l}d_{l-1}\big{)}$.

##  5\. Experiments

In this section, we conduct extensive experiments to evaluate the
effectiveness of DMNS 111Code & data at https://github.com/ntkien1904/DMNS.git
on several benchmark datasets, and analyze several key aspects of the model.

###  5.1. Experimental Setup

Table 1. Summary of datasets. Datasets | Nodes | Edges | Features | Property
---|---|---|---|---
Cora | 2708 | 5429 | 1433 | homophilous
Citeseer | 3327 | 4732 | 3703 | homophilous
Coauthor-CS | 18333 | 163788 | 6805 | homophilous
Actor | 7600 | 30019 | 932 | heterophilous
Table 2. Evaluation of link prediction against baselines using GCN as the base encoder. Methods | Cora | Citeseer | Coauthor-CS | Actor
---|---|---|---|---
MAP | NDCG | MAP | NDCG | MAP | NDCG | MAP | NDCG
GCN | .742 ± .003 | .805 ± .003 | .735 ± .011 | .799 ± .008 | .823 ± .004 | .867 ± .003 | .521 ± .004 | .634 ± .003
GVAE |  .783 ± .003 |  .835 ± .002 | .743 ± .004 | .805 ± .003 | .843 ± .011 | .882 ± .008 |  .587 ± .004 |  .684 ± .003
PNS | .730 ± .008 | .795 ± .006 | .748 ± .006 | .809 ± .005 | .817 ± .004 | .863 ± .003 | .517 ± .006 | .631 ± .006
DNS | .735 ± .007 | .799 ± .005 |  .777 ± .005 |  .831 ± .004 | .845 ± .003 | .883 ± .002 | .558 ± .006 | .663 ± .005
MCNS | .756 ± .004 | .815 ± .003 | .750 ± .006 | .810 ± .004 | .824 ± .004 | .868 ± .004 | .555 ± .005 | .659 ± .004
GraphGAN | .739 ± .003 | .802 ± .002 | .740 ± .011 | .803 ± .008 | .818 ± .007 | .863 ± .005 | .534 ± .007 | .644 ± .005
ARVGA | .732 ± .011 | .797 ± .009 | .689 ± .005 | .763 ± .004 | .811 ± .003 | .858 ± .002 | .526 ± .012 | .638 ± .009
KBGAN | .615 ± .004 | .705 ± .003 | .568 ± .006 | .668 ± .005 |  .852 ± .002 |  .888 ± .002 | .472 ± .003 | .596 ± .002
SEAL | .751 ± .007 | .812 ± .005 | .718 ± .002 | .784 ± .002 | .850 ± .001 | .886 ± .001 | .536 ± .001 | .641 ± .001
ScaLed | .676 ± .004 | .752 ± .003 | .630 ± .004 | .712 ± .003 | .828 ± .001 | .869 ± .001 | .459 ± .001 | .558 ± .001
DMNS |  .793 ± .003 |  .844 ± .002 |  .790 ± .004 |  .841 ± .003 |  .871 ± .002 |  .903 ± .001 |  .600 ± .002 |  .696 ± .002
∗Best is bolded and runner-up underlined.
Table 3. Evaluation of link prediction on DMNS with various base encoders. Methods | Cora | Citeseer | Coauthor-CS | Actor
---|---|---|---|---
MAP | NDCG | MAP | NDCG | MAP | NDCG | MAP | NDCG
GAT | .766 ± .006 | .824 ± .004 | .767 ± .007 | .763 ± .062 | .833 ± .003 | .874 ± .002 | .479 ± .004 | .603 ± .003
DMNS-GAT |  .813 ± .004 |  .859 ± .003 |  .788 ± .007 |  .840 ± .006 |  .851 ± .002 |  .889 ± .002 |  .573 ± .007 |  .675 ± .005
SAGE | .598 ± .014 | .668 ± .013 | .622 ± .012 | .713 ± .009 | .768 ± .005 | .826 ± .004 | .486 ± .004 | .604 ± .003
DMNS-SAGE |  .700 ± .007 |  .773 ± .005 |  .669 ± .013 |  .749 ± .010 |  .843 ± .004 |  .883 ± .003 |  .582 ± .017 |  .682 ± .013

#### Datasets

We employ four public graph datasets, summarized in Table 1. Cora and Citeseer
(Yang et al., 2016) are two citation networks, where each node is a document
and the edges represent citation links. Coauthor-CS (Shchur et al., 2018) is a
co-authorship network, where each node is an author and an edge exists if they
co-authored a paper. Actor (Pei et al., 2020) is an actor co-occurrence
network, where each node denotes an actor and each edge connects two actors
both occurring on the same Wikipedia page.

#### Baselines

We employ a comprehensive suite of baselines for the link prediction task. (1)
_Classical GNNs_ : GCN (Kipf and Welling, 2017) and VGAE (Kipf and Welling,
2016), which are classical GNNs models for link prediction. (2) _Heuristic
negative sampling_ : PNS (Mikolov et al., 2013), DNS (Zhang et al., 2013) and
MCNS (Yang et al., 2020), which employ various heuristics to retrieve hard
negative examples. (3) _Generative adversarial methods_ : GraphGAN (Wang et
al., 2018), ARVGA (Pan et al., 2018) and KBGAN (Cai and Wang, 2017), which
leverages GANs to learn the negative distribution and select hard examples.
(4) _Subgraph-based GNNs_ : SEAL (Zhang and Chen, 2018) and ScaLed (Louis et
al., 2022), which utilizes local subgraphs surrounding the candidate nodes.
See Appendix A for a more detailed description of the baselines.

#### Task setup and evaluation

On each graph, we randomly split its links for training, validation and
testing following the proportions 90%:5%:5%. Note that the graphs used in
training are reconstructed from only the training links. We adopt a ranking-
based link prediction during testing. Given a query node $v$, we sample a
positive candidate $u$ such that $(v,u)$ is a link in the test set, and
further sample 9 nodes that are not linked to $v$ as negative candidates. For
evaluation, we rank the 10 candidate nodes based on their dot product with the
query node $v$. Based on the ranked list, we report two ranking-based metrics,
namely, NDCG and MAP (Liu et al., 2021b), averaged over five runs.

#### Parameters and settings

For our model DMNS, by default, we employ GCN (Kipf and Welling, 2017) as the
base encoder. To further evaluate the effectiveness of DMNS on different
encoders, we also conduct experiments using GAT (Veličković et al., 2018) and
GraphSAGE (SAGE) (Hamilton et al., 2017). For GCN, we employ two layers with
dimension 32. For the diffusion model, we use two FiLM layers with output
dimension 32, set the total time step as $T=50$, and assign the variances
$\beta$ as constants increasing linearly from $10^{-4}$ to 0.02. We set the
weights of negative examples to $\\{1,0.9,0.8,0.7\\}$ in Eq. (16),
corresponding to time steps
$\\{\frac{T}{10},\frac{T}{8},\frac{T}{4},\frac{T}{2}\\}$. For all baselines,
we adopt the same GNN architecture and settings as in DMNS for a fair
comparison. They also employ a log-sigmoid objective (Hamilton et al., 2017)
consistent with our link prediction loss. For those baselines which do not
propose a negative sampling method, we perform the uniform negative sampling
if required for training. Additional settings of our method and the baselines
can be found in Appendix B.

###  5.2. Evaluation on Link Prediction

We evaluate the performance of DMNS on link prediction against various
baselines and with different base encoders.

#### Comparison with baselines.

We report the results of DMNS and the baselines in Table 2. Overall, DMNS
significantly outperforms competing baselines across all datasets and metrics.
The results indicate that our strategy of multi-level negative sampling is
effective. We make three further observations. First, heuristic negative
sampling methods generally improve the performance upon the base encoder with
uniform sampling (i.e., GCN), showing the utility of harder examples. Among
them, DNS and MCNS are more robust as they leverage more sophisticated
heuristics to select better negative examples. Second, GAN-based methods are
often worse than classical GNNs, which is potentially due to the difficulty in
training GANs. For instance, GVAE achieves a relatively competitive
performance, while ARVGA—its variant with an adversarial regularizer—suffers
from a considerable drop. Third, subgraph-based GNNs achieve mixed results,
suggesting that the local structures can be effective to some extent yet high-
quality negative examples are still needed for training.

#### Evaluation with other encoders

We further utilize GAT and SAGE as the base encoders, and the corresponding
model DMNS-GAT and DMNS-SAGE, respectively. From the results reported in Table
3, we observe that DMNS achieves significant improvements compared to its
corresponding base encoders in all cases, which demonstrates the flexibility
of DMNS to effectively work with various encoders.

###  5.3. Model Analyses

Finally, we investigate various aspects of DMNS through ablation studies,
parameter sensitivity analysis and visualization on the four datasets.

#### Ablation studies.

The first study investigates the ablation on model design, to demonstrate the
effectiveness of each module in DMNS. We compare with two variants: (1)
_unconditional_ diffusion, by removing the condition on the query node (and
its associated neighboring nodes) such that the diffusion model now generates
arbitrary node embeddings from noises; (2) _unweighted_ negative examples, by
setting all $w_{i}=1$ in Eq. (16). We report the results in Fig. 3a, and make
the following observations. First, unconditional diffusion performs
significantly worse than its conditional counterpart, since the generated
negative examples are not optimized for a specific query node. Second, the
unweighted variant also exhibits a drop in performance, implying that negative
examples from different levels of hardness have different importance.

The second study involves the ablation on the sampling choice. Instead of
sampling a mixture of examples from well-spaced time steps, we now only use
output from a single time step for each query node. We vary the time step
between $\frac{T}{10}$ and $\frac{T}{2}$, as shown in Fig. 3b. The performance
of each single time step varies but all are worse than combining them together
(i.e., DMNS), demonstrating the effectiveness of multi-level sampling. It is
also observed that smaller time steps ($\frac{T}{10},\frac{T}{8}$) often
outperform larger time steps ($\frac{T}{4},\frac{T}{2}$), indicating that
harder examples could be more useful.

Figure 3. Ablation studies. Figure 4. Parameter sensitivity. Figure 5.
Histogram of embedding distances from query.

#### Parameter studies

We showcase the impact of two important hyperparameters, including the total
time step $T$ for the diffusion process and the scheduling for diffusion
variances $\beta$. In Fig. 4a, we observe that increasing total time steps may
not boost the performance while requiring more computation. In Fig. 4b,
results indicate that linear scheduling achieves the best performance among
several standard scheduling policies, including cosine and sigmoid scheduling
(Nichol and Dhariwal, 2021).

#### Visualization

We further investigate the quality of the negative examples generated by DMNS.
For each query node, we construct a positive set consisting of its neighbors,
a negative set consisting of negative examples from DMNS from time
$\frac{T}{10}$ to $\frac{T}{2}$, and a second negative set consisting of
uniformly sampled nodes from the graph. We calculate the Euclidean distance
between the query node embedding and the embeddings in each set, and plot the
empirical distribution (histogram) of the distances in Fig. 5. For the
positive set, we naturally expect smaller distances to the query node. For the
negative sets, the distance can be regarded as a proxy to hardness: Smaller
distances from the query node imply harder examples. As we can see, examples
of DMNS from $\frac{T}{10}$ to $\frac{T}{4}$ are harder than uniform sampling,
but are not too hard (i.e., not closer to the query node than the positives)
to impair the performance (Faghri et al., 2017; Xuan et al., 2020), while
using larger steps from $\frac{T}{2}$ may not bring more performance gain than
uniform sampling. Overall, utilizing multi-level samples allows to capture a
wide range of hardness levels for negative sampling.

##  6\. Conclusion

In this paper, we investigated a novel strategy of multi-level negative
sampling for graph link prediction. Existing methods aim to retrieve hard
examples heuristically or adversarially, but they are often inflexible or
difficult to control the “hardness”. In response, we proposed a novel sampling
method named DMNS based on a conditional diffusion model, which empowers the
sampling of negative examples at different levels of hardness. In particular,
the hardness can be easily controlled by sampling from different time steps of
the denoise process within the diffusion model. Moreover, we showed that DMNS
largely obeys the sub-linear positivity principle for robust negative
sampling. Finally, we conducted extensive experiments to demonstrate the
effectiveness of DMNS. A limitation of our work is the focus on the
effectiveness and robustness of negative sampling for link prediction. Hence,
future direction can aim to optimize the sampling process for efficiency on
large-scale graphs. Alternatively, we can investigate the potential of
diffusion model in other graph learning tasks such as node classification.

###### Acknowledgements.

This research / project is supported by the Ministry of Education, Singapore,
under its Academic Research Fund Tier 2 (Proposal ID: T2EP20122-0041). Any
opinions, findings and conclusions or recommendations expressed in this
material are those of the author(s) and do not reflect the views of the
Ministry of Education, Singapore.

## References

  * (1)
  * Adamic and Adar (2003) Lada A Adamic and Eytan Adar. 2003.  Friends and neighbors on the web.  _Social networks_ 25, 3 (2003), 211–230.
  * Ahrabian et al. (2020) Kian Ahrabian, Aarash Feizi, Yasmin Salehi, William L Hamilton, and Avishek Joey Bose. 2020.  Structure aware negative sampling in knowledge graphs.  _arXiv preprint arXiv:2009.11355_ (2020).
  * Arjovsky and Bottou (2017) Martin Arjovsky and Léon Bottou. 2017.  Towards principled methods for training generative adversarial networks.  _arXiv preprint arXiv:1701.04862_ (2017).
  * Arjovsky et al. (2017) Martin Arjovsky, Soumith Chintala, and Léon Bottou. 2017.  Wasserstein generative adversarial networks. In _International conference on machine learning_. PMLR, 214–223.
  * Cai and Wang (2017) Liwei Cai and William Yang Wang. 2017.  Kbgan: Adversarial learning for knowledge graph embeddings.  _arXiv preprint arXiv:1711.04071_ (2017).
  * Chamberlain et al. (2022) Benjamin Paul Chamberlain, Sergey Shirobokov, Emanuele Rossi, Fabrizio Frasca, Thomas Markovich, Nils Hammerla, Michael M Bronstein, and Max Hansmire. 2022.  Graph Neural Networks for Link Prediction with Subgraph Sketching.  _arXiv preprint arXiv:2209.15486_ (2022).
  * Chen et al. (2020) Liang Chen, Yuanzhen Xie, Zibin Zheng, Huayou Zheng, and Jingdun Xie. 2020.  Friend recommendation based on multi-social graph convolutional network.  _IEEE Access_ 8 (2020), 43618–43629.
  * Dhariwal and Nichol (2021) Prafulla Dhariwal and Alexander Nichol. 2021.  Diffusion models beat gans on image synthesis.  _Advances in Neural Information Processing Systems_ 34 (2021), 8780–8794.
  * Ding et al. (2020) Jingtao Ding, Yuhan Quan, Quanming Yao, Yong Li, and Depeng Jin. 2020.  Simplify and robustify negative sampling for implicit collaborative filtering.  _Advances in Neural Information Processing Systems_ 33 (2020), 1094–1105.
  * Faghri et al. (2017) Fartash Faghri, David J Fleet, Jamie Ryan Kiros, and Sanja Fidler. 2017.  Vse++: Improving visual-semantic embeddings with hard negatives.  _arXiv preprint arXiv:1707.05612_ (2017).
  * Gong et al. (2022) Shansan Gong, Mukai Li, Jiangtao Feng, Zhiyong Wu, and LingPeng Kong. 2022.  Diffuseq: Sequence to sequence text generation with diffusion models.  _arXiv preprint arXiv:2210.08933_ (2022).
  * Grover and Leskovec (2016) Aditya Grover and Jure Leskovec. 2016.  node2vec: Scalable feature learning for networks. In _Proceedings of the 22nd ACM SIGKDD international conference on Knowledge discovery and data mining_. 855–864.
  * Gutmann and Hyvärinen (2010) Michael Gutmann and Aapo Hyvärinen. 2010.  Noise-contrastive estimation: A new estimation principle for unnormalized statistical models. In _Proceedings of the thirteenth international conference on artificial intelligence and statistics_. JMLR Workshop and Conference Proceedings, 297–304.
  * Hamilton et al. (2017) Will Hamilton, Zhitao Ying, and Jure Leskovec. 2017.  Inductive representation learning on large graphs. In _NeurIPS_. 1024–1034.
  * Ho et al. (2022a) Jonathan Ho, William Chan, Chitwan Saharia, Jay Whang, Ruiqi Gao, Alexey Gritsenko, Diederik P Kingma, Ben Poole, Mohammad Norouzi, David J Fleet, et al. 2022a.  Imagen video: High definition video generation with diffusion models.  _arXiv preprint arXiv:2210.02303_ (2022).
  * Ho et al. (2020) Jonathan Ho, Ajay Jain, and Pieter Abbeel. 2020.  Denoising diffusion probabilistic models.  _Advances in Neural Information Processing Systems_ 33 (2020), 6840–6851.
  * Ho et al. (2022b) Jonathan Ho, Chitwan Saharia, William Chan, David J Fleet, Mohammad Norouzi, and Tim Salimans. 2022b.  Cascaded Diffusion Models for High Fidelity Image Generation.  _J. Mach. Learn. Res._ 23, 47 (2022), 1–33.
  * Hoogeboom et al. (2022) Emiel Hoogeboom, Vıctor Garcia Satorras, Clément Vignac, and Max Welling. 2022.  Equivariant diffusion for molecule generation in 3d. In _International Conference on Machine Learning_. PMLR, 8867–8887.
  * Hu et al. (2019) Binbin Hu, Yuan Fang, and Chuan Shi. 2019.  Adversarial learning on heterogeneous information networks. In _Proceedings of the 25th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining_. 120–129.
  * Huang et al. (2021) Tinglin Huang, Yuxiao Dong, Ming Ding, Zhen Yang, Wenzheng Feng, Xinyu Wang, and Jie Tang. 2021.  Mixgcf: An improved training method for graph neural network-based recommender systems. In _Proceedings of the 27th ACM SIGKDD Conference on Knowledge Discovery & Data Mining_. 665–674.
  * Jo et al. (2022) Jaehyeong Jo, Seul Lee, and Sung Ju Hwang. 2022.  Score-based generative modeling of graphs via the system of stochastic differential equations. In _International Conference on Machine Learning_. PMLR, 10362–10383.
  * Kipf and Welling (2016) Thomas N Kipf and Max Welling. 2016.  Variational graph auto-encoders.  _arXiv preprint arXiv:1611.07308_ (2016).
  * Kipf and Welling (2017) Thomas N Kipf and Max Welling. 2017.  Semi-supervised classification with graph convolutional networks. In _ICLR_.
  * Li et al. (2023) Jianxin Li, Xingcheng Fu, Shijie Zhu, Hao Peng, Senzhang Wang, Qingyun Sun, S Yu Philip, and Lifang He. 2023.  A robust and generalized framework for adversarial graph embedding.  _IEEE Transactions on Knowledge and Data Engineering_ (2023).
  * Liu et al. (2021a) Zemin Liu, Yuan Fang, Yong Liu, and Vincent W Zheng. 2021a.  Neighbor-anchoring adversarial graph neural networks.  _IEEE Transactions on Knowledge and Data Engineering_ 35, 1 (2021), 784–795.
  * Liu et al. (2021b) Zemin Liu, Trung-Kien Nguyen, and Yuan Fang. 2021b.  Tail-gnn: Tail-node graph neural networks. In _Proceedings of the 27th ACM SIGKDD Conference on Knowledge Discovery & Data Mining_. 1109–1119.
  * Louis et al. (2022) Paul Louis, Shweta Ann Jacob, and Amirali Salehi-Abari. 2022.  Sampling Enclosing Subgraphs for Link Prediction. In _Proceedings of the 31st ACM International Conference on Information & Knowledge Management_. 4269–4273.
  * Lü and Zhou (2011) Linyuan Lü and Tao Zhou. 2011.  Link prediction in complex networks: A survey.  _Physica A: statistical mechanics and its applications_ 390, 6 (2011), 1150–1170.
  * Luo (2022) Calvin Luo. 2022.  Understanding diffusion models: A unified perspective.  _arXiv preprint arXiv:2208.11970_ (2022).
  * Mikolov et al. (2013) Tomas Mikolov, Ilya Sutskever, Kai Chen, Greg S Corrado, and Jeff Dean. 2013.  Distributed representations of words and phrases and their compositionality.  _Advances in neural information processing systems_ 26 (2013).
  * Nguyen et al. (2023) Trung-Kien Nguyen, Zemin Liu, and Yuan Fang. 2023.  Link Prediction on Latent Heterogeneous Graphs. In _Proceedings of the ACM Web Conference 2023_. 263–273.
  * Nichol et al. (2021) Alex Nichol, Prafulla Dhariwal, Aditya Ramesh, Pranav Shyam, Pamela Mishkin, Bob McGrew, Ilya Sutskever, and Mark Chen. 2021.  Glide: Towards photorealistic image generation and editing with text-guided diffusion models.  _arXiv preprint arXiv:2112.10741_ (2021).
  * Nichol and Dhariwal (2021) Alexander Quinn Nichol and Prafulla Dhariwal. 2021.  Improved denoising diffusion probabilistic models. In _International Conference on Machine Learning_. PMLR, 8162–8171.
  * Page et al. (1999) Lawrence Page, Sergey Brin, Rajeev Motwani, and Terry Winograd. 1999.  _The PageRank citation ranking: Bringing order to the web._ Technical Report. Stanford infolab.
  * Pan et al. (2018) Shirui Pan, Ruiqi Hu, Guodong Long, Jing Jiang, Lina Yao, and Chengqi Zhang. 2018.  Adversarially regularized graph autoencoder for graph embedding.  _arXiv preprint arXiv:1802.04407_ (2018).
  * Pei et al. (2020) Hongbin Pei, Bingzhe Wei, Kevin Chen-Chuan Chang, Yu Lei, and Bo Yang. 2020.  Geom-gcn: Geometric graph convolutional networks.  _arXiv preprint arXiv:2002.05287_ (2020).
  * Perez et al. (2018) Ethan Perez, Florian Strub, Harm De Vries, Vincent Dumoulin, and Aaron Courville. 2018.  Film: Visual reasoning with a general conditioning layer. In _Proceedings of the AAAI Conference on Artificial Intelligence_ , Vol. 32.
  * Saharia et al. (2022) Chitwan Saharia, William Chan, Huiwen Chang, Chris Lee, Jonathan Ho, Tim Salimans, David Fleet, and Mohammad Norouzi. 2022.  Palette: Image-to-image diffusion models. In _ACM SIGGRAPH 2022 Conference Proceedings_. 1–10.
  * Shchur et al. (2018) Oleksandr Shchur, Maximilian Mumme, Aleksandar Bojchevski, and Stephan Günnemann. 2018.  Pitfalls of graph neural network evaluation.  _arXiv preprint arXiv:1811.05868_ (2018).
  * Sohl-Dickstein et al. (2015) Jascha Sohl-Dickstein, Eric Weiss, Niru Maheswaranathan, and Surya Ganguli. 2015.  Deep unsupervised learning using nonequilibrium thermodynamics. In _International Conference on Machine Learning_. PMLR, 2256–2265.
  * Song and Ermon (2019) Yang Song and Stefano Ermon. 2019.  Generative modeling by estimating gradients of the data distribution.  _Advances in neural information processing systems_ 32 (2019).
  * Song and Ermon (2020) Yang Song and Stefano Ermon. 2020.  Improved techniques for training score-based generative models.  _Advances in neural information processing systems_ 33 (2020), 12438–12448.
  * Song et al. (2020) Yang Song, Jascha Sohl-Dickstein, Diederik P Kingma, Abhishek Kumar, Stefano Ermon, and Ben Poole. 2020.  Score-based generative modeling through stochastic differential equations.  _arXiv preprint arXiv:2011.13456_ (2020).
  * Tang et al. (2015) Jian Tang, Meng Qu, Mingzhe Wang, Ming Zhang, Jun Yan, and Qiaozhu Mei. 2015.  Line: Large-scale information network embedding. In _Proceedings of the 24th international conference on world wide web_. 1067–1077.
  * Veličković et al. (2018) Petar Veličković, Guillem Cucurull, Arantxa Casanova, Adriana Romero, Pietro Lio, and Yoshua Bengio. 2018.  Graph attention networks. In _ICLR_.
  * Wang et al. (2018) Hongwei Wang, Jia Wang, Jialin Wang, Miao Zhao, Weinan Zhang, Fuzheng Zhang, Xing Xie, and Minyi Guo. 2018.  Graphgan: Graph representation learning with generative adversarial nets. In _Proceedings of the AAAI conference on artificial intelligence_ , Vol. 32.
  * Xuan et al. (2020) Hong Xuan, Abby Stylianou, Xiaotong Liu, and Robert Pless. 2020.  Hard negative examples are hard, but useful. In _Computer Vision–ECCV 2020: 16th European Conference, Glasgow, UK, August 23–28, 2020, Proceedings, Part XIV 16_. Springer, 126–142.
  * Yang et al. (2022) Ling Yang, Zhilong Zhang, Yang Song, Shenda Hong, Runsheng Xu, Yue Zhao, Yingxia Shao, Wentao Zhang, Bin Cui, and Ming-Hsuan Yang. 2022.  Diffusion models: A comprehensive survey of methods and applications.  _arXiv preprint arXiv:2209.00796_ (2022).
  * Yang et al. (2016) Zhilin Yang, William Cohen, and Ruslan Salakhudinov. 2016.  Revisiting semi-supervised learning with graph embeddings. In _International conference on machine learning_. PMLR, 40–48.
  * Yang et al. (2020) Zhen Yang, Ming Ding, Chang Zhou, Hongxia Yang, Jingren Zhou, and Jie Tang. 2020.  Understanding negative sampling in graph representation learning. In _Proceedings of the 26th ACM SIGKDD international conference on knowledge discovery & data mining_. 1666–1676.
  * Ying et al. (2018) Rex Ying, Ruining He, Kaifeng Chen, Pong Eksombatchai, William L Hamilton, and Jure Leskovec. 2018.  Graph convolutional neural networks for web-scale recommender systems. In _Proceedings of the 24th ACM SIGKDD international conference on knowledge discovery & data mining_. 974–983.
  * You et al. (2020) Yuning You, Tianlong Chen, Yongduo Sui, Ting Chen, Zhangyang Wang, and Yang Shen. 2020.  Graph contrastive learning with augmentations.  _Advances in neural information processing systems_ 33 (2020), 5812–5823.
  * Zhang and Chen (2018) Muhan Zhang and Yixin Chen. 2018.  Link prediction based on graph neural networks.  _Advances in neural information processing systems_ 31 (2018).
  * Zhang et al. (2013) Weinan Zhang, Tianqi Chen, Jun Wang, and Yong Yu. 2013.  Optimizing top-n collaborative filtering via dynamic negative item sampling. In _Proceedings of the 36th international ACM SIGIR conference on Research and development in information retrieval_. 785–788.
  * Zhou et al. (2009) Tao Zhou, Linyuan Lü, and Yi-Cheng Zhang. 2009.  Predicting missing links via local information.  _The European Physical Journal B_ 71 (2009), 623–630.

## Appendices

###  A. Details of Baselines

In this section, we describe each baseline in more details.

  * •

_Classical GNNs_ : GCN (Kipf and Welling, 2017) applies mean-pooling
aggregates information for the target node by over its neighbors. GVAE (Kipf
and Welling, 2016) proposes an unsupervised training for graph representation
learning by reconstructing graph structure.

  * •

_Heuristic negative sampling_ : PNS is adopted from word2vec (Mikolov et al.,
2013), where the distribution for negative samples is calculated by the
normalized node degree. DNS (Zhang et al., 2013) dynamically select hard
negative examples from candidates ranked by the current prediction scores.
MCNS (Yang et al., 2020) proposes a Metropolis-Hasting method for negative
sampling based on the sub-linear positivity principle.

  * •

_Generative adversarial methods_ : GraphGAN (Wang et al., 2018) samples hard
negative examples from the learned connectivity distribution through
adversarial training. KBGAN (Cai and Wang, 2017) samples high quality negative
examples for knowledge graph embeddings by the generator that learns to
produce distribution for negative candidates conditioned on the positives.
ARVGA (Pan et al., 2018) integrates an adversarial regularizer to GVAE to
enforce the model on generating realistic samples.

  * •

_Subgraph-based GNNs_ : SEAL (Zhang and Chen, 2018) proposes to sample k-hops
local subgraphs surrounding two candidates for link prediction. ScaLED (Louis
et al., 2022) improves SEAL by utilizing random walks for efficient subgraphs
sampling.

For other base GNN encoders, GAT (Veličković et al., 2018) employs self-
attention mechanism to produce learnable weights to each neighbor of target
node during aggregation. SAGE (Hamilton et al., 2017) first aggregates
information from target node neighbors, then concatenates with the node itself
to obtain the node embedding.

###  B. Model Settings

For all the approaches, we use a two layers GCN as encoder with output
dimension as 32 for fair comparison to conduct link prediction. For our model
DMNS, we set learning rate as 0.01, dropout ratio as 0.1. For each query nodes
$v$, we sample a neighbor set $N_{v}=20$ for diffusion training.

For GCN, we set the hidden dimension as 32. For GAT, we use three attention
heads with hidden dimension for each head as 32. For SAGE, we use mean
aggregator for concatenation and set its hidden dimension as 32. In PNS, we
utilize node degree to calculate the popularity distribution and normalize it
to 3/4${}^{\text{th}}$ power. For DNS, we set the number of negative
candidates as 10 and select 1 negative sample as the highest ranked node. For
MCNS, we set the DFS length as 5 and the proposal distribution $q(y|x)$ is
mixed betwwen uniform and k-nearest nodes sampling wih $p=0.5$. In these
methods we set learning rate as 0.01.

For GraphGAN, we use pretrained GCN to obtain initialized node embeddings and
set the learning rates for discriminator $D$ and generator $G$ as 0.0001. In
each iteration, we set the number of gradient updates for both $D$ and $G$ as
1 for best results. For KBGAN, we adapt it to the homogeneous graph setting by
neglecting all relation types and use a single embedding to represent a
universal relation. We further set the number of negative examples as 20 and
number of gradient updates for both $D$ and $G$ as 1. For ARVGA, learning
rates for encoder $E$ and discriminator $D$ are 0.005 and 0.001 respectively,
hidden dimension for $D$ is 64, the number of gradient updates $D$ for each
iteration is 5.

For Subgraph-based GNNs approaches, we replace their original objective Binary
Cross Entropy with our Log Sigmoid and do not utilize additional node
structural embeddings for fair evaluation. For SEAL, we set learning rate as
0.0001, the number of hop $k$ equals 1. On ScaLed, we set random walk length
as 3 and number of walks as 20.

###  C. Environment

All experiments are conducted on a workstation with a 12-core CPU, 128GB RAM,
and 2 RTX-A5000 GPUs. We implemented DMNS using Python 3.8 and Pytorch 1.13 on
Ubuntu-20.04.

###  D. Impact of Generated Samples

To evaluate effectiveness of generated samples, we conduct additional ablation
studies on the two separate sets with heuristic and generated ones. Results
from Fig. 6 show that generated samples can achieve better results than
heuristic ones, and combining both sets yields the best performance.

Figure 6. Ablation on generated samples.

###  E. Scalability

To analyze time complexity, we sample five subgraphs from the large-scale OGB
graph (Nguyen et al., 2023) with sizes ranging from 20k to 100k nodes. We
report the training time (in second) of one epoch of DMNS in the Table 5. We
observe that the training time of DMNS increases linearly with the graph size
in terms of the number of nodes, demonstrating its scalablity to larger
graphs.

Nodes | Edges | Time per epoch
---|---|---
20k | 370k | 1.61s
40k | 810k | 3.18s
60k | 1.2M | 4.82s
80k | 1.6M | 6.27s
100k | 1.8M | 7.51s
Table 4. Train time

Methods | #Parameters
---|---
MCNS | 5k
GVAE | 48k
SEAL | 78k
ARVGA | 56k
DMNS | 13k
Table 5. Model Size

Regarding space complexity, we provide the model size of DMNS (i.e., number of
parameters) compared to several representative baselines, all using a base
encoder of 2-layer GCN with hidden dimension as 32. Table 5 shows that our
model is comparatively efficient in model size compared to other advanced
baselines like GVAE, SEAL or ARVGA.