Enhancing Graph Neural Networks with Quantum Computed Encodings
Transformers are increasingly employed for graph data, demonstrating
competitive performance in diverse tasks. To incorporate graph information into
these models, it is essential to enhance node and edge features with positional
encodings. In this work, we propose novel families of positional encodings
tailored for graph transformers. These encodings leverage the long-range
correlations inherent in quantum systems, which arise from mapping the topology
of a graph onto interactions between qubits in a quantum computer. Our
inspiration stems from the recent advancements in quantum processing units,
which offer computational capabilities beyond the reach of classical hardware.
We prove that some of these quantum features are theoretically more expressive
for certain graphs than the commonly used relative random walk probabilities.
Empirically, we show that the performance of state-of-the-art models can be
improved on standard benchmarks and large-scale datasets by computing tractable
versions of quantum features. Our findings highlight the potential of
leveraging quantum computing capabilities to potentially enhance the
performance of transformers in handling graph data.
Slimane Thabet, Romain Fouilland, Mehdi Djellabi, Igor Sokolov, Sachin Kasture, Louis-Paul Henry, Loïc Henriet
# Enhancing Graph Neural Networks with

###### Abstract

Transformers are increasingly employed for graph data, demonstrating competitive performance in diverse tasks. To incorporate graph information into these models, it is essential to enhance node and edge features with positional encodings. In this work, we propose novel families of positional encodings tailored for graph transformers. These encodings leverage the long-range correlations inherent in quantum systems, which arise from mapping the topology of a graph onto interactions between qubits in a quantum computer. Our inspiration stems from the recent advancements in quantum processing units, which offer computational capabilities beyond the reach of classical hardware. We prove that some of these quantum features are theoretically more expressive for certain graphs than the commonly used relative random walk probabilities. Empirically, we show that the performance of state-of-the-art models can be improved on standard benchmarks and large-scale datasets by computing tractable versions of quantum features. Our findings highlight the potential of leveraging quantum computing capabilities to potentially enhance the performance of transformers in handling graph data.

## 1 Introduction

Graph machine learning (GML) is an expanding field of research with applications in chemistry (Gilmer et al., 2017), biology (Zitnik et al., 2018), drug design (Konaklieva, 2014), social networks (Scott, 2011), computer vision (Harchaoui and Bach, 2007) and science (Sanchez-Gonzalez et al., 2020; Xu et al., 2018). In the past few years, significant effort has been put into the design of Graph Neural Networks (GNNs) (Hamilton). The objective is to learn suitable representations that enable efficient solutions to the original problem.

To that end, a large number of models have been developed in the past few years (Kipf and Welling, 2016; Hamilton et al., 2018; Velickovic et al., 2018). While the prevalent approach for constructing GNNs relies on the Message Passing (MP) mechanism (Gilmer et al., 2017), this approach exhibits several recognized limitations, with the most significant being its theoretical expressivity. Indeed, two graphs that are indistinguishable via the Weisfeiler-Lehman (WL) test will lead to the same MP Neural Network (MPNN) output (Morris et al., 2019). Another limitation arises from the fact that MPNNs are more effective when dealing with homophilic data. This is based on the underlying assumption that nodes that are similar, either in structure or features, are more likely to be related. A study by (Zhu et al., 2020) demonstrates that MPNNs encounter difficulties when applied to heterophilic graphs. Finally, MPNNs are prone to over-smoothing (Chen et al., 2020) as well as over-squashing (Topping et al., 2021). These latter aspects constitute serious limitations for the datasets that exhibit long-range dependencies (Dwivedi et al., 2022).

The research community is actively exploring solutions to address these limitations. The key idea is to expand aggregation beyond neighbouring nodes by incorporating information related to the entire graph or a more extensive portion of it. Graph Transformers were created according to these requirements, with success on standard benchmarks (Ying et al., 2021; Rampasek et al., 2022). Among the myriad of proposed architectures, the Graph Inductive Bias Transformer (GRIT) (Ma
et al., 2023) stands out for its impressive generalization capacity. This stems from its independence from the MP mechanism and its utilization of multiple positional encodings in its architecture. While these features make it a good candidate for overcoming the aforementioned limitations, the authors relied on discrete \(k\)-step random walks to initialize the PE tensor. These random walks constitute to this day the most widely adopted choice (He et al., 2023; Dwivedi et al., 2022a), since the main alternative, the matrix of the Laplacian eigenvectors, is invariant by sign flip of each vector, resulting in \(2^{k}\) possible choices.

The goal of this work is to leverage new types of structural features emerging from quantum physics as positional encodings. The rapid development of quantum computers during the previous years provides the opportunity to compute features that would be otherwise intractable. These features contain complex topological characteristics of the graph, and their inclusion has the potential to enhance the model's quality, reduce training or inference time, and decrease energy consumption.

The paper is organized as follows: In Section 2, we provide a concise overview of the existing research on graph transformers, along with references to the latest developments in quantum graph machine learning. Section 3 delves into the core theoretical aspects of this work. It covers quantum mechanics basics for readers unfamiliar with the topic, details the way to construct a quantum state from a graph and explains why quantum states can provide relevant information that is hard to compute with a classical computer. Additionally, we introduce our central proposal, a framework for both static and dynamic positional encoding based on quantum correlations. Finally, Section 4 presents the outcomes of our numerical experiments and includes discussions of the results.

Figure 1: Summary of our method. **(a)** Our hybrid quantum-classical framework utilizes a classical computer for parameter optimization and employs a hybrid model using a Quantum Processing Unit (QPU) and a CPU and/or GPU, denoted as classical Processing Unit (cPU). In our quantum graph NN, we initialize QPU at a quantum state \(|\psi_{0}\rangle\), apply a mixing Hamiltonian \(\hat{\mathcal{H}}_{\mathrm{M}}\) evolution for a duration \(\theta\), and utilize a Hamiltonian \(\hat{\mathcal{H}}_{\mathcal{G}}\) evolution for the graph feature map with a duration \(t\). \(K\) layers are used to obtain a sufficiently expressive quantum model. Finally, the output is obtained by measuring correlators, e.g., \(\langle Z_{i}Z_{j}\rangle\). See Section 3.1 for details. **(b)** Static or trainable PE is constructed for a graph \(\mathcal{G}\) via **(c)** (quantum) random walk (static PE) or a quantum graph NN (static/trainable PE), which computes quantum correlations. Note that our PEs are not restricted to classical models (such as the transformer studied in this work) but are also applicable to all quantum models.
## 2 Related works

### Graph Transformers

Efforts have been made in the community to go beyond MPNNs due to several issues (Zhu et al., 2020; Chen et al., 2020; Topping et al., 2021). Inspired by the success of transformers in natural language processing (Vaswani et al., 2017; Alayrac et al., 2022), new architectures of GNNs have been proposed to allow an all-to-all aggregation between the nodes of the graphs, and called graph transformers (GT) (Dwivedi & Bresson, 2020; Dwivedi et al., 2021; Ramapask et al., 2022; Kreuzer et al., 2021; Zhang et al., 2023; Ma et al., 2023). However, due to the quadratic cost of computing the attention process, they are not applicable to large-scale graphs of millions of nodes and more. It has been shown that GTs that include graph inductive biases such as MP modules perform better than those that do not (Rampasek et al., 2022; Ma et al., 2023).

### Positional and Structural Encoding

Positional or structural embeddings are features computed from the graph that are concatenated to original node or edge features to enrich GNN architectures (either MPNN or GT). These two terms are used interchangeably in the literature and we denote them as "positional encodings" (PEs) in the rest of this work. PEs can include random walk probabilities (Rampasek et al., 2022; Ma et al., 2023), spectral information (Dwivedi et al., 2020; Rampasek et al., 2022; Kreuzer et al., 2021), shortest path distances (Li et al., 2018), or heat kernels (Mialon et al., 2021). They can also be learned (Dwivedi et al., 2021). We detail below the most common PEs used in the literature.

**Laplacian Eigenvectors.** The spectral information of the graph can be used as PE, more precisely the eigenvectors of the Laplacian matrix. If one takes a line graph, it almost corresponds to positional embeddings in the transformer architecture for sequences. The main issue of this encoding is to ensure that the model remains invariant by changing the sign of eigenvectors, which has been solved by (Lim et al., 2022).

**Relative Random Walk Probabilities (RRWP).** The authors of (Ma et al., 2023) introduced the RRWP with which they initialize their model. For a graph \(\mathcal{G}\), let \(A\) be the adjacency matrix and \(D\) the degree matrix. Let \(P\) be a 3 dimensional tensor such that \(P_{k,i,j}=(M^{k})_{ij}\) with \(M=D^{-1}A\). For each pair of node \((i,j)\), we associate the vector \(P_{i,i,j}\), i.e., the concatenation of the probabilities for all \(k\) to get from node \(i\) to node \(j\) in \(k\) steps in a random walk. \(P_{i,i}\) is the same as the Random Walk Structural Encodings (RWSE) defined in (Rampasek et al., 2022). The authors of (Ma et al., 2023) highlight the benefits of RRWP. They prove that the Generalized Distance WL (GD-WL) test introduced by (Zhang et al., 2023) with RRWP is strictly more powerful than GD-WL test with the shortest path distance, and they prove universal approximation results of multi-layer perceptrons (MLP) initialized with RRWP.

### Quantum Computing for Graph Machine Learning

Using quantum computing for Machine Learning on graphs has already been proposed in several works, as reviewed in (Tang et al., 2022). The authors of (Verdon et al., 2019) realized learning tasks by using a parameterized quantum circuit depending on a Hamiltonian whose interactions share the topology of an input graph. Comparable ideas were used to build graph kernels from the output of quantum procedures, for photonic (Schuld et al., 2020) as well as neutral atom quantum processors (Henry et al., 2021). The latter was successfully implemented on quantum hardware (Albrecht et al., 2023). The architectures proposed in these papers were entirely quantum and only relied on classical computing for the optimization of variational parameters. By contrast, in what we propose here quantum dynamics only plays a role in the aggregation phase of a larger entirely classical architecture. Such a hybrid model presents the advantage of gaining access to hard-to-access graph topological features through quantum dynamics while benefiting from the power of well-known existing classical architectures.
Methods and theory

In this section, we outline the process of mapping graphs to a quantum state of a QPU. To extract graph features, we introduce correlators and define the concept of the ground state for a quantum graph representation. Finally, we explore an alternative approach for extracting graph features using quantum random walks (QRW) and their advantages over classical analogues.

### Quantum Graph Machine Learning

**The graph as a quantum state.** We explain in this subsection how to create a quantum state that contain relevant information about the graph. More details about quantum information processing can be found in (Nielsen and Chuang, 2002). The _quantum state_\(\ket{\psi}\) of a system of \(N\) qubits can be represented as a vector of unit norm in \(\mathbb{C}^{2^{N}}\). Quantum states are modified through the action of _operators_, that can be represented as hermitian matrices of size \(2^{N}\times 2^{N}\). Its dynamics obeys the Schrodinger equation \(-i\frac{d\ket{\psi}}{dt}=\hat{\mathcal{H}}\ket{\psi}\), where the operator \(\hat{\mathcal{H}}\) is the _Hamiltonian_ of the system, with solution \(\ket{\psi(t)}=\mathcal{T}\exp\left[-i\int_{0}^{t}\hat{\mathcal{H}}(\tau)d\tau \right]\ket{\psi(0)}\), with \(\mathcal{T}\) the time-ordering operator. An operator \(\hat{\mathcal{O}}\) that can be measured is called an _observable_, and its eigenvalues correspond to possible outcome of its measurement. Its expectation value on the quantum state \(\ket{\psi}\) is the scalar \(\bra{\hat{\mathcal{O}}}=\bra{\psi}\hat{\mathcal{O}}\ket{\psi}\), where \(\bra{\psi}\) is the conjugate transpose of \(\ket{\psi}\). The _Pauli matrices_ are defined as follows: \(I=\left(\begin{smallmatrix}1&0\\ 0&1\end{smallmatrix}\right)\), \(X=\left(\begin{smallmatrix}0&1\\ 1&0\end{smallmatrix}\right)\), \(Y=\left(\begin{smallmatrix}0&-i\\ i&0\end{smallmatrix}\right)\), \(Z=\left(\begin{smallmatrix}1&0&-1\\ 0&-1\end{smallmatrix}\right)\). They form a basis of hermitian matrices of size \(2\times 2\). A _Pauli string_ of size \(N\) is an operator that can be written as the Kronecker product of \(N\) Pauli matrices. We will note Pauli strings by their non-trivial Pauli operations. For instance, in a system of 5 qubits, \(X_{0}Y_{3}=X\otimes I\otimes I\otimes Y\otimes I\). We associate a graph \(\mathcal{G}(\mathcal{V},\mathcal{E})\), to a quantum state \(\ket{\psi_{\mathcal{G}}}\) of \(\ket{\mathcal{V}}\) qubits containing information about \(\mathcal{G}\) via a hamiltonian \(\hat{\mathcal{H}}_{\mathcal{G}}\) of the form

\[\hat{\mathcal{H}}_{\mathcal{G}}=\sum_{(i,j)\in\mathcal{E}}\hat{\mathcal{H}}_{ij} \tag{1}\]

where \(\hat{\mathcal{H}}_{ij}\) is an Pauli string acting non-trivially on \(i\) and \(j\) only. We will be focusing on the Ising hamiltonian \(\hat{\mathcal{H}}^{I}=\sum_{(i,j)\in\mathcal{E}}Z_{i}Z_{j}\) and the XY hamiltonian \(\hat{\mathcal{H}}^{XY}=\sum_{(i,j)\in\mathcal{E}}X_{i}X_{j}+Y_{i}Y_{j}\). We will note \(\ket{0}\) and \(\ket{1}\) the two eigenstates (or eigenvectors) of \(Z\) with respective eigenvalues 1 and -1, and we will use \(\left\{\ket{\mathbf{b}}=\bigotimes_{i=1}^{N}\ket{b_{i}}\right\}_{\mathbf{b} \in\{0,1\}^{N}}\) as a basis of the \(2^{N}\)-dimensional space of quantum states. We consider the quantum state obtained by alternated action of \(p\) layers of \(\hat{\mathcal{H}}_{\mathcal{G}}\) and a _mixing_ hamiltonian \(\hat{\mathcal{H}}_{M}\) (that doesn't commute with \(\hat{\mathcal{H}}_{\mathcal{G}}\), for instance \(\hat{\mathcal{H}}_{M}\propto\sum_{i}Y_{i}\))

\[\ket{\psi_{\mathcal{G}}(\mathbf{\theta})}=\prod_{k=1}^{p}\left(\mathrm{e}^{-i \hat{\mathcal{H}}_{M}\theta_{k}}\mathrm{e}^{-i\hat{\mathcal{H}}_{G}t_{k}} \right)\mathrm{e}^{-i\hat{\mathcal{H}}_{M}\theta_{0}}\ket{\psi_{0}}, \tag{2}\]

where \(\mathbf{\theta}=(\theta_{0},t_{0},\theta_{1},t_{1},\ldots\theta_{p})\) is a real vector of parameters. The choice of these states is motivated by their similarity with the _Trotterized_ dynamics of a lot of quantum systems(Suzuki, 1976).

**Correlation.** The correlations (or _correlators_) \(C_{ij}\) of local operators \(\hat{\mathcal{O}}_{i}\) and \(\hat{\mathcal{O}}_{j}\) acting respectively on qubits \(i\) and \(j\) can be defined either as the expectation value of their product \(\langle\hat{\mathcal{O}}_{i}\hat{\mathcal{O}}_{j}\rangle\), or their covariance \(\langle\hat{\mathcal{O}}_{i}\hat{\mathcal{O}}_{j}\rangle-\langle\hat{\mathcal{ O}}_{i}\rangle\langle\hat{\mathcal{O}}_{j}\rangle\) (note that the orders matters if \(\hat{\mathcal{O}}_{i}\) and \(\hat{\mathcal{O}}_{j}\) don't commute). In the rest of the paper, we will indifferently call correlation the two former expressions, and give precisions when necessary. We will be focusing on the case where \(\hat{\mathcal{O}}_{i}\) is a Pauli string of length 1 (i.e. \(X_{i}\), \(Y_{i}\) or \(Z_{i}\)).

**Ground state.** The ground state of a system is defined as the lowest-energy eigenstate of its hamiltonian (when it is degenerate, one considers the _ground state manifold_\(\mathbb{H}_{GS}\)). Ground state properties are widely studied in many-body physics and their properties depend on the topology of the graph. Preparing this state is the purpose of quantum annealing (Das and Chakrabarti, 2008). When using neutral atom quantum processors (Henriet et al., 2020), one can natively address hamiltonians of the form \(\hat{\mathcal{H}}_{\mathcal{G}}=\sum_{(i,j)\in\mathcal{E}}J_{ij}(Z_{i}-\alpha_{i }I)(Z_{j}-\alpha_{j}I)\), with \(\alpha_{i}\) real coefficients. Its eigenstates are
the basis states \(|\mathbf{b}\rangle\) described above. In the case where \(\alpha_{i}=1-\delta/(2z_{i})\) with \(z_{i}=\sum_{j|(i,j)\in\mathcal{E}}J_{ij}\) and \(J_{ij}=1/4\), the eigenenergies (or eigenvalues) are \(E(\mathbf{b})=\sum_{i,j\in\mathcal{E}}b_{i}b_{j}-\delta\sum_{i=1}^{N}b_{i}\). When \(0<\delta<1\), this is the cost function associated with the maximum independent set problem, a NP-hard problem (Garey & Johnson, 1979). In the absence of degeneracy-lifting or symmetry-breaking effects, a quantum annealing scheme would prepare a symmetric, equal-weight superposition of all maximum independent sets. With that in mind, we will call _ground state of the graph_ the state \(|\psi_{GS}\rangle=\frac{1}{\sqrt{|\mathbb{H}_{GS}|}}\sum_{\mathbf{b}\in\mathbb{ H}_{GS}}|\mathbf{b}\rangle\).

**Classical and Quantum Walks.** Quantum walks, as introduced by (Aharonov et al., 1993), differ fundamentally from classical random walks by evolving through unitary processes, allowing for interference between different trajectories. These walks manifest in two primary types: _continuous-time quantum walks_ (CQRW) (Farhi & Gutmann, 1998; Rossi et al., 2017) and _discrete-time quantum walks_ (DQRW) (Lourett et al., 2010). Discrete classical random walks on \(\mathcal{G}(\mathcal{V},\mathcal{E})\) use the probability matrix \(M=D^{-1}A\) for node transitions over a walk of length \(K\), resulting in the probability distribution \(P_{K}=M^{K}P_{0}\)(Aharonov et al., 2001). Note that this approach is utilized in RRWP encodings (Ma et al., 2023). In the continuous case, CQRW can be viewed as a natural extension of continuous-time classical random walks (CRW). In CRW, the probability of a walker being at vertex \(i\) and time \(t\) is represented as \(p_{i}(t)\), which follows the differential equation \(\frac{\mathrm{d}}{\mathrm{d}t}p_{i}(t)=-\sum_{j}G_{ij}p_{j}(t)\). Here, the infinitesimal generator \(G_{ij}=-\gamma\) if an edge exists between nodes \(i\) and \(j\), and \(0\) otherwise, with diagonal elements \(G_{ii}=k_{i}\gamma\) determined by the node degree \(k_{i}\). Considering now a quantum evolution with a graph Hamiltonian \(\hat{\mathcal{H}}_{\mathcal{G}}\), given a \(2^{N}\)-dimensional Hilbert space of \(N\) qubits, the Schrodinger equation which governs the evolution of a quantum state \(|\psi_{\mathcal{G}}\rangle\) when projected onto a state \(|i\rangle\) is given as

\[\underbrace{i\frac{\mathrm{d}}{\mathrm{d}t}\langle i|\psi_{\mathcal{G}}(t) \rangle=\sum_{j}\langle i|\hat{\mathcal{H}}_{\mathcal{G}}|j\rangle\langle j| \psi_{\mathcal{G}}(t)\rangle}_{\mathrm{Quantum}}\longleftrightarrow\underbrace{ \frac{\mathrm{d}}{\mathrm{d}t}p_{i}(t)=-\sum_{j}G_{ij}p_{j}(t)}_{\mathrm{ Classical}}. \tag{3}\]

Note the similarity between the differential equations of CQRW and CRW. A quantum analogue of CRW can be obtained by taking \(\langle i|\hat{\mathcal{H}}_{\mathcal{G}}|j\rangle=G_{ij}\). The probabilities are preserved as the sum of amplitude squared, \(\sum_{i}|\langle i|\psi_{\mathcal{G}}(t)\rangle|^{2}=1\), in the quantum case, instead of \(\sum_{i}p_{i}(t)=1\) in the classical case. This difference between evolution of probabilities(which are real) and evolution of amplitudes(which are complex) leads to interesting differences between dynamics of classical and quantum walks. Using this formalism, any quantum evolution can be thought of as a CQRW (Childs et al., 2002). Notably, quantum walks have demonstrated exponential hitting time advantage for graphs like hypercubes (Kempe, 2002) and glued binary trees (Childs et al., 2003). These results have been recently extended for more general hierarchical graphs (Balasubramanian et al., 2023). For an overview, refer to (Kempe, 2003).

### Positional encodings with quantum features

In this section, we detail our proposals to incorporate quantum features in GNN models, and we discuss the potential benefits and drawbacks. Our methods can be roughly divided in two categories: quantum features that are used as static positional encodings and are precomputed at the begining of the procedure and quantum features that can be dynamically trained.

#### 3.2.1 Static position encoding

**Eigenvectors of the correlation on the ground state.** We propose to use the correlation matrix \(C_{ij}=\langle Z_{i}Z_{j}\rangle\) on the ground state of the graph defined in Sec. 3.1. Since this matrix is symmetric with non-negative eigenvalues, it can formally be used in the same place as the Laplacian matrix in graph learning models. Hence, we use the eigenvectors of this correlation matrix in the same way Laplacian eigenvectors (LE) are used in other architectures of graph transformers. Instead of taking the eigenvectors with the lowest eigenvalues as for the Laplacian eigenmaps, we take the ones with highest eigenvalues, since they are the ones in which most of the information about the correlation matrix is contained. We expect to face the same challenges about the sign ambiguity (Dwivedi et al., 2021; Kreuzer et al., 2021), and to implement the same techniques to alleviate them (Lim et al., 2022).
\(k\)**-particles quantum random walks (\(k\)-QRW).** In this work, we introduce the \(k\)-particles (or walkers) random walk positional encoding, that can be can be obtained using \(\hat{\mathcal{H}}^{XY}\). We note \(\hat{\mathcal{H}}^{XY}_{k}\) the XY hamiltonian restricted to the \(k\) particles subspace \(\mathbb{H}_{k}\) (_i.e._ the span of states \(\left|\mathbf{b}\right\rangle\) of hamming weight \(k\), noted \(\left|i_{1}\ldots i_{k}\right\rangle\), parameterized by \(k\) integers \(i_{1}\ldots i_{k}\in\{0,1\}^{k}\)). For a 1-particle QRW, we calculate the probability \([X^{(1)}(t)]_{ij}=|\left\langle j\right|e^{-i\hat{\mathcal{H}}^{XY}_{k}t}|i \rangle|^{2}\) to find particle at node \(j\) coming from node \(i\) after time \(t\). Similarly for a 2-particle QRW, we calculate \([X^{(2)}(t)]_{ij}=|\left\langle ij\right|e^{-i\hat{\mathcal{H}}^{XY}_{k}t}| \psi_{i}\rangle|^{2}\), where \(|i,j\rangle\in\hat{\mathcal{H}}^{XY}_{k}\) is the state with walkers at nodes \(i\) and \(j\) and \(|\psi_{i}\rangle\in\hat{\mathcal{H}}^{XY}_{k}\) the initial state. As choices for the initial distribution we propose to use some localised state \(\left|\psi_{\text{init}}\right\rangle\propto|ij\rangle\), or the uniform distribution over all pairs of nodes \(\left|\psi_{\text{init}}\right\rangle\propto\sum_{(i,j)\in\mathcal{V}^{2}|ij \rangle}j)\), or the uniform distribution over the edges of the original graph \(\left|\psi_{\text{init}}\right\rangle\propto\sum_{(i,j)\in\mathcal{E}}|ij\rangle\). From these we obtain the positional encodings using \(\mathbf{P}_{ij}=[I,X^{n_{w}}(t_{1}),X^{n_{w}}(t_{2})...X^{n_{w}}(t_{K})]_{ij}\), where \(n_{w}=1,2\) is the number of walkers. From the symmetries of \(\hat{\mathcal{H}}^{XY}\), this 2-QRW can be viewed as a 1-QRW on the set of 2-particle states (see A.3.2). From there, we consider a _discrete_ 2-particle quantum-inspired RW (2-QiRW) encoding that reads \(\mathbf{P}_{ij}=\left[\langle ij|\,(\hat{\mathcal{H}}^{XY}_{2})^{k}\,|\,\psi_{ \text{init}}\rangle\,|k\in[0,K]\right]_{ij}\).

#### 3.2.2 Learnable positional encodings

Here we consider a specific case of equation 2, with \(p=1\), \(\theta_{0}=-\theta_{1}=\theta\), \(\hat{\mathcal{H}}_{M}\propto\sum_{i\in\mathcal{V}}Y_{i}\) and \(\hat{\mathcal{H}}_{\mathcal{G}}=\sum_{(i,j)\in\mathcal{E}}Z_{i}Z_{j}-\delta \sum_{i\in\mathcal{V}}Z_{i}\). A similar setting was implemented on neutral atom QPU in (Albrecht et al., 2023), where \(\hat{\mathcal{H}}_{M}\propto\sum_{i\in\mathcal{V}}X_{i}+\varepsilon\hat{ \mathcal{H}}_{\mathcal{G}}\) (with \(\varepsilon\lesssim 1\)) due to hardware constraints. We then use the covariance matrix of the number of occupation observable \(\hat{\mathcal{O}}_{i}=\frac{1}{2}(I-Z_{i})\), which equals one when the \(i-\)th atom is in its excited state, and 0 otherwise. This is one of the few particular cases where we can recover a closed formula for the correlation matrix used as PE, see Appendix A.2.1 for its full expression.

The goal in this section is to learn the positional encoding by training a GNN, or any permutation equivariant NN, to find an optimised value of the parameters (\(\theta\), \(t\) and \(\delta\)) involved in our PE. The training of this module is carried out jointly with that of the transformer layers, and the input PE is updated after each backward pass. This allows a custom value of the parameters for each graph in the dataset. To this end, we train a GNN \(P_{\mathbf{W}}(\mathbf{A},\mathbf{X})\), with \(\mathbf{A}\) being the adjacency matrix, \(\mathbf{X}\) the feature matrix of the nodes in the graph, and \(\mathbf{W}\) the parameter set of the NN. We obtain the PE parameters as \((\theta||t||\delta)=P_{\mathbf{W}}(\mathbf{A},\mathbf{X})\in\mathbb{R}^{3\times k}\) (we learn \(k\) triples, encoding as many correlations matrices, which we concatenate in one tensor as it is done in the original GRIT paper). Here we chose \(\mathbf{X}\) as the pairwise graph distance matrix between the nodes, such that we only consider its structural features. Since the positional encoding discussed here is obtained from a special case of QAOA, a key aspect of this approach concerns the initial values taken by \(\theta,t\) and \(\delta\)(Egger et al., 2021). Among the many possible initialization protocols, the authors in (Jain et al., 2022) used GNN to find a warm start in the non-convex energy landscape, prior to an optimization through quantum annealing, for the Max-Cut problem. In our case, the training is carried out classically (only because we recover a closed formula), and its extension as a quantum NN module in a transformer-like architecture is yet to be investigated in future works.

### GNN models with quantum aggregation

In the previous subsection, we explained how to integrate positional encodings coming from a quantum processing unit into a graph transformer of graph neural network model. In this subsection, we propose to directly include the quantum correlations as a trainable part of the model. Given \(H^{l}\) the node features matrix at the layer \(l\), the node features matrix at the next layer is computed with the following formula :

\[H^{l+1}=\sigma((A(\theta)H^{l}||H^{l})W) \tag{4}\]

where \(A(\theta)\) is the quantum attention matrix. Given a quantum graph state, we compute for every pair of nodes \((i,j)\) the vector of 2-bodies observables \(C_{ij}=[\langle Z_{i}Z_{j}\rangle,\langle X_{i}X_{j}\rangle,\langle Y_{i}Y_{j} \rangle,\langle X_{i}Z_{j}\rangle,\langle X_{i}Y_{j}\rangle,\langle Y_{i}Z_{j} \rangle,\langle X_{j}Z_{i}\rangle,\langle X_{j}Y_{i}\rangle,\langle Y_{j}Z_{i} \rangle]^{T}\). The quantum atten
tion matrix is computed by taking a linear combination of the previous correlation vector and optionally a softmax over the lines. More details are given in the Appendix B.

### Theoretical arguments for QRWs

**Two-interacting-particle QRW are more expressive than RRWP.** WL tests (Weisfeiler and Leman, 1968) and their extensions (Morris et al., 2019; Grohe and Otto, 2015; Grohe, 2017) are widely used algorithms for distinguishing graphs (Graph Isomorphism). A standard WL test generates node colorings as a function of the node neighbourhood. An extension of this called the generalized distance WL (GD-WL) test had been introduced by (Zhang et al., 2023) as a generalization of the WL isomorphism test. Let \(\mathcal{G}\) be a graph on a set of vertices \(\mathcal{V}\) and let \(d(u,v)\) be the distance between nodes \(u,v\). Just like the WL test, the GD-WL test starts with an initial color \(h_{0}(u)\) for a node u. Then at iteration \(l\), the color \(h_{l}(u)\) is computed by \(h_{l}(u)=\texttt{HASH}(\{(h_{-1}(v),d(u,v)),\,v\in\mathcal{V}\})\) where HASH is a hash function. The algorithm stops when the colors don't change after an iteration. The authors of Zhang et al., 2023 show that if equipped with some distances such as the shortest path distance (SPD), the GD-WL test is strictly more powerful than the WL test. The authors of (Ma et al., 2023) show that the GD-WL test equipped with RRWP is strictly more powerful than the GD-WL test with the SPD distance. We can show however (see Appendix A.3.1 for a proof), that GD-WL test with RRWP embeddings cannot distinguish non-isomorphic strongly regular graphs (see (Gamble et al., 2010) for the definition of strongly regular graphs) :

**Theorem 1**: _GD-WL test with RRWP embedding fails to distinguish non isomorphic strongly regular graphs. 2 particle-QRW can distinguish some of them._

## 4 Experiments

We performed several experiments to assess the capabilities of quantum features to improve existing GNN models. We report in the following subsections the most significant results. Less conclusive experiments are detailed in Appendices C.4, C.3.

### Experiments on RW models

In this subsection, we test concatenating the QRW encodings to the RRWP in the GRIT model (Ma et al., 2023). We compute the (continuous) 1-CQRW for \(K\) random times and the discrete 2-QRW for \(K\) steps. Those encodings are computed numerically since they are still tractable for graphs below 200 nodes compared to the higher order \(k\)-QRW ones. We benchmark our method on 7 datasets from (Dwivedi et al., 2020), following the experimental setup of (Rampasek et al., 2022) and (Ma et al., 2023). Our method is compared to many other architectures and the results directly taken from (Ma et al., 2023). We do not perform an extensive hyperparameter search for each architecture and only run ourselves the GRIT model by taking the same hyperparameters as the authors. The experiments are done by building on the codebase of (Ma et al., 2023) which is itself built on (Rampasek et al., 2022). More details about the protocol and hyperparameters can be found in C.1, and more details about the datasets can be found in Appendix D. The results are included in Table 1. Our methods performs better on ZINC, MNIST and CIFAR10 than all others, and comes second for PATTERN and CLUSTER. We also benchmark our methods on large-scale datasets, ZINC-full (a bigger version of ZINC (Irwin et al., 2012)) and PCQM4MV2 (Hu et al., 2021). As before we only compute GRIT and report other results from (Ma et al., 2023). Our methods performs the best among all others.

### Experiments on learning the positional encoding

Here we discuss results related to the learned parameters of the quantum PE, which were mainly run on the ZINC dataset. We observe from the comparison of these results, displayed in the right part of table 2, and those of random parameters (table 1) that the performances are significantly lower in this case. The main reason for this is the incapacity of the proposed model, to efficiently explore the space of the encoded parameters \(\theta,t\) and \(\delta\). More details about this are provided in the Appendix A.1 Given such limitations, the question arises as to whether the proposed positional encoding, via learned or randomly fixed parameters, plays any role in the training of the transformer model, or whether the attention scheme proposed in (Ma et al., 2023) along with the graph external features is enough to
obtain the same performances. For that we compare our results with three different GRIT+RRWP models on randomized Zinc datasets. In the first one, we remove the structural information while keeping the degree sequence in each graph (configuration model of the graph (Newman, 2010)), in the second we remove any trace of structural information by randomly replacing each graph in the dataset by a random graph with the same number of nodes and edges, and in the final case we keep the structural information but randomly permute the feature vectors of each graph, such that we isolate the contribution of the structural features from that of the external features. In each case, we benchmark these results using a 2 layered graph neural network (with GAT convolutional layers), to encode the parameters of the PE tensor. The results of these comparisons are showcased in the right part of table 2, and further details about the randomization process in Appendix A.2

\begin{table}
\begin{tabular}{l c c c c c} \hline \hline
**Model** & **ZINC** & **MNIST** & **CIFAR10** & **PATTERN** & **CLUSTER** \\ \cline{2-6}  & **MAE\(\downarrow\)** & **Accuracy\(\uparrow\)** & **Accuracy\(\uparrow\)** & **Accuracy\(\uparrow\)** & **Accuracy\(\uparrow\)** \\ \hline GCN & \(0.367\pm 0.011\) & \(90.705\pm 0.218\) & \(55.710\pm 0.381\) & \(71.892\pm 0.334\) & \(68.498\pm 0.976\) \\ GIN & \(0.526\pm 0.051\) & \(96.485\pm 0.252\) & \(55.255\pm 1.527\) & \(58.387\pm 0.136\) & \(64.716\pm 1.553\) \\ GAT & \(0.384\pm 0.007\) & \(55.35\pm 0.205\) & \(64.223\pm 0.455\) & \(78.271\pm 0.186\) & \(70.587\pm 0.447\) \\ GateGCN & \(0.282\pm 0.015\) & \(97.340\pm 0.143\) & \(67.312\pm 0.311\) & \(85.568\pm 0.088\) & \(73.840\pm 0.326\) \\ DCN & \(0.168\pm 0.003\) & \(-\) & \(72.838\pm 0.417\) & \(86.680\pm 0.034\) & \(-\) \\ GIN-AK+ & \(0.080\pm 0.001\) & \(-\) & \(72.19\pm 0.13\) & \(86.850\pm 0.057\) & \(-\) \\ SAN & \(0.139\pm 0.006\) & \(-\) & \(-\) & \(86.581\pm 0.037\) & \(76.961\pm 0.65\) \\ K-Sabgraph SAT & \(0.094\pm 0.008\) & \(98.715\pm 0.087\) & \(68.702\pm 0.409\) & \(86.821\pm 0.020\) & \(79.232\pm 0.348\) \\ EGT & \(0.108\pm 0.009\) & \(98.051\pm 0.126\) & \(72.298\pm 0.336\) & \(86.685\pm 0.059\) & \(78.016\pm 0.180\) \\ GPS & \(0.059\pm 0.002\) & \(98.108\pm 0.111\) & \(76.468\pm 0.881\) & \(87.196\pm 0.076\) & \(80.026\pm 0.277\) \\ GRIT (our run) & \(0.060\pm 0.002\) & \(98.164\pm 0.054\) & \(76.198\pm 0.744\) & \(90.405\pm 0.232\) & \(79.856\pm 0.156\) \\ \hline
**GRIT 1-CQRW** & \(0.058\pm 0.002\) & \(98.108\pm 0.111\) & \(76.347\pm 0.704\) & \(87.205\pm 0.040\) & \(78.895\pm 0.1145\) \\
**GRIT 2-QQRW** & \(0.059\pm 0.004\) & \(98.204\pm 0.048\) & \(76.442\pm 1.07\) & \(90.165\pm 0.446\) & \(79.777\pm 0.171\) \\ \hline \hline \end{tabular}
\end{table}
Table 1: Test performance in five benchmarks from (Dwivedi et al., 2020). We show the mean \(\pm\) s.d. of 4 runs with different random seeds as in (Ma et al., 2023). Highlighted are the top first, second, and third results. Models are restricted to \(\sim 500K\) parameters for ZINC, PATTERN, CLUSTER \(\sim 100K\) for MNIST and CIFAR10. We compare our model to our run of GRIT and indicate the results obtained by the authors for information. Figures other than the last 3 lines are taken from (Ma et al., 2023). Models in **bold** are our models.

\begin{table}
\begin{tabular}{c l c|c c} \hline \hline
**Method** & **Model** & **ZINC-full** & **PCQMAMv2** & **Model** & **Zine** \\  & & (MAE \(\downarrow\)) & (MAE \(\downarrow\)) & & (MAE \(\downarrow\)) \\ \hline \multirow{5}{*}{MPNNs} & GIN & \(0.088\pm 0.002\) & \(0.1195\) & \(-\) & \(-\) \\  & GraphSAGE & \(0.126\pm 0.003\) & \(-\) & \(-\) & \(-\) \\  & GAT & \(0.111\pm 0.002\) & \(-\) & \(-\) & \(-\) \\  & GCN & \(0.113\pm 0.002\) & \(0.1195\) & \(-\) & \(-\) \\ \hline PE-GNN & SignNet & \(0.024\pm 0.003\) & \(-\) & \(-\) & \(-\) \\ \hline \multirow{3}{*}{
\begin{tabular}{c} Graph \\ Transformers \\ \end{tabular} } & Graphormer-URPE & \(0.052\pm 0.005\) & \(0.0864\) & **Rand. fests +RRWP** & \(0.393\pm 0.012\) \\  & Graphormer-URPE & \(0.028\pm 0.002\) & \(-\) & **Rand. strvet +RRWP** & \(0.245\pm 0.009\) \\  & Graphormer-GD & \(0.025\pm 0.004\) & & **Config. model+RRWP** & \(0.156\pm 0.004\) \\  & GPS-medium & \(-\) & \(0.0858\) & **2 layers GAT+QCorr** & \(0.134\pm 0.017\) \\  & GRIT (MA \(\downarrow\), 2023) & \(0.023\pm 0.001\) & \(0.0859\) & **2 layers SAGE+QCorr** & \(0.130\pm 0.003\) \\  & GRIT (our run) & \(0.025\pm 0.002\) & \(0.0842\) & **2 layers Transf. +QCorr** & \(0.127\pm 0.014\) \\  & **GRIT 1-CQRW (ours)** & \(0.025\pm 0.003\) & \(0.0947^{*}\) & **2 layers GCN +QCorr** & \(0.111\pm 0.003\) \\  & **GRIT 2-QQRW (ours)** & \(0.023\pm 0.002\) & \(0.0838\) & \(-\) & \(-\) \\ \hline \hline \end{tabular}
\end{table}
Table 2: Left part : test performance on ZINC-full (Irwin et al., 2012) and PCQM4MV2. For ZINC-full, we show the mean and s.d of 4 runs with different random seeds and we limit the model to \(\sim 500K\) parameters. For PCQM4MV2 we show the output of a single run due to computation time. We compare our model to our run of GRIT and indicate the results obtained by the authors for information. Figures other than GRIT are taken from (Ma et al., 2023). Highlighted are the top first, second, and third results. Right part: results related to the isolation of the effects of structural and external graph features in the GRIT model, on the small zinc dataset. For the result with \({}^{*}\) only 55 epochs were used.
### Synthetic experiments

In this section, we provide two examples of datasets with a binary graph classification task for which the use of the correlation matrix on the ground state as defined in 3.2.1 is more powerful than other commonly used features like the eigenvectors of the laplacian matrix or the diagonal of the random walk matrix. We name our datasets special pattern (S-PATTERN) and cross ladder (C-LADDER), and we provide a detail of their construction in appendices D.1.1 and D.1.2. The idea is to construct graphs that will exhibit very different Ising ground states but similar spectral properties or random walk transition probabilities. We illustrate the differences between the encodings in Appendix C.2. We train classical models on these datasets with LEs and random walk embeddings (RWSE) as node features, and we compare it to the same models with eigenvectors of the correlation on the ground state. We experiment with GCN and GPS with many hyperparameters combinations. More details on the protocols can be found in C.2 We also benchmark the GRIT model with RRWP. The results are shown in table 3. The quantum encoding models achieve 100% accuracy in both cases whereas the models with LE or RWSE achieve only 57% maximum. The GRIT model achieves 78% for the S-PATTERN, and 100% for the C-LADDER.

### Discussion

We performed several experiments comparing the quantum encodings to the classical ones. Including the quantum walk features into state of the art models improves their performances on most of the datasets tested. It is not surprising that the method works well for datasets for which random walks are known to be relevant features like ZINC (Rampasek et al., 2022). We only limited ourselves to versions of quantum features that are efficiently computable, and we were able to show a small gain in performances compared to state of the art models. It is then believable that using quantum features that cannot be classically accessible could lead to a great improvement of models, if quantum hardware can be made widely available. We were able to engineer artificial datasets with for which classical approaches have difficulties to perform the associated binary classification tasks. GPS model fails both of the tasks whereas GRIT model is successful on the S-PATTERN task and mildly successful on the C-LADDER task. We think GRIT is successful on the S-PATTERN dataset because the degree distribution is different for the two classes, and GRIT specifically uses a degree scaler in its architecture. We have addressed, albeit at a preliminary stage, the issue of learning parameters related to the positional encoding. We have shown that even with non-optimal parameters (noted QCor on the right-hand side of the table 2), the model still performs better than when undergoing structural randomization, which underlines the discriminative capacity of the QCor positional encoder, and that the performance of the GRIT-QCorr model is not exclusively due to the graph external features.

## 5 Conclusion

In this paper we have investigated how quantum computing architectures can be used to construct new families of graph neural networks. This study involved measuring observables like correlations and probabilities for a quantum system whose hamiltonian has the same topology as the graph of interest. We then integrated these observables as positional encodings and used them in different classical graph neural network architectures. We also used them as attention layers in graph transformers. We proved that some positional encodings that use quantum features are theoretically more expressive than ones based on simple random walk, on certain classes of graphs. Our experiments show that state of the art models can already be enhanced with restricted quantum features that are classically efficient to compute. This study provides strong indications that the full leverage of quantum hardware can

\begin{table}
\begin{tabular}{c c c c c} \hline \hline
**Dataset** & **(GCN/GPS)-LE** & **(GCN/GPS)-RWSE** & **GRIT-RRWP** & **(GCN/GPS)-Q** \\ \hline C-LADDER & \(57\) & \(56\) & \(78\) & \(100\) \\ S-PATTERN & \(56\) & \(55\) & \(100\) & \(100\) \\ \hline \hline \end{tabular}
\end{table}
Table 3: Results on synthetic data. We show the accuracy on the test set for each datasets. For each positional encoding, we show the score of the best model among all the combinations of hyperparameters tested. LE : Laplacian Eignevectors, RWSE: Random Walk Structural Encodings, Q: Quantum, eigenvectors of the correlation on the ground state.
lead to development of high-performance architectures for certain tasks. In particular, Neutral Atom quantum hardware (Henriet et al., 2020) would be particularly suited to the type of time-dependent Hamiltonian we described here. Furthermore, we can create artificial classification tasks that are easily solvable with quantum enhanced models while classical models struggle. While the exact capabilities of our approach have to be explored, the results we obtain show that quantum enhanced GNNs are a promising family of models that could be fully exploited with near term quantum hardware.

## Acknowledgements

We would like to thank Jacob Bamberger, Constantin Dalyac and Vincent Elfving for useful discussions and comments. We especially thank Jacob Bamberger for introducing us a method to generate non-isomorphic graphs indistinguishable by the WL test.

## Reproducibility statement

The code to reproduce the experiments discussed in the main text is included in the supplementary materials with indications to run it. We also detail in the appendices C.1 C.2 C.4 the details of the protocols.