Jet Discrimination with Quantum Complete Graph Neural Network

Graphs are ubiquitous data structures that represent relationships and connections between entities. Analyzing and extracting valuable information from graph-structured data is a fundamental challenge in modern data science and machine learning. To address this challenge, Graph Neural Networks (GNNs) have emerged as a powerful and versatile framework for learning from graph-structured data [86, 87, 88].

A graph $G$ is described by its set of nodes and edges. Let $N$ denote the number of nodes, and $E_{ij}$ the edge from the $i$-th node to the $j$-th node. Throughout this paper, the graphs are assumed to be undirected ($E_{ij}=E_{ji}$) and unweighted (all edges have equal weights). Furthermore, a graph is considered complete [81] if all pairs of nodes are connected. A complete graph that is undirected and unweighted can also be viewed as a set. Each node is associated with a feature vector $\mathbf{x}\in\mathbb{R}^{d}$, where $d$ is the dimension of the features.

Graphs are permutation invariant, meaning that reordering the indices of the nodes does not alter the information contained in the graph. The general structure of permutation-invariant models has been studied in [89, 34]. In this work, we focus on a specific type of GNN, the Message Passing Graph Neural Network (MPGNN) [90], which provides a simple and intuitive approach to designing GNNs. The MPGNN is constructed by iterating the following formula¹¹1The formula is adapted from the PyTorch Geometric documentation https://pytorch-geometric.readthedocs.io/en/latest/tutorial/create_gnn.html. The $e_{j,i}$ term is omitted here as we only consider undirected and unweighted graphs.

$$\mathbf{x}^{(k)}_{i}=\gamma^{(k)}\left(\mathbf{x}^{(k-1)}_{i},\bigoplus_{\begin{subarray}{c}j\in\mathcal{N}(i)\end{subarray}}\Phi^{(k)}(\mathbf{x}^{(k-1)}_{i},\mathbf{x}^{(k-1)}_{j})\right),$$

(1)

where $\Phi^{(k)}$ extracts information from the neighboring nodes $\mathcal{N}(i)$ of the $i$-th node, and $\gamma^{(k)}$ updates the node features in the $k$-th iteration. As pointed out in [89], summation over a set is the sufficient and necessary condition for satisfying permutation invariance. This concept can be generalized by the aggregation function $\bigoplus$, which is typically chosen as SUM, MEAN, MAX, or MIN²²2Note that MAX and MIN can be expressed as summation over the $p$-norm, which can be absorbed into $\Phi$., among others. A brief discussion on the relation between state-of-the-art models and the MPGNN is provided in Appendix A.

In the noisy intermediate-scale quantum (NISQ) era [92], variational quantum algorithms [78, 79, 80] present an intuitive approach to implementing quantum neural networks. These networks utilize a variational quantum circuit (VQC) with tunable parameters updated via classical optimization routines. Typically, an $n$-qubit VQC can be expressed as:

$$f(\mathbf{x},\bm{\theta})=\bra{0}^{\otimes n}U^{\dagger}(\mathbf{x},\bm{\theta})PU(\mathbf{x},\bm{\theta})\ket{0}^{\otimes n}$$

for some Pauli string observable $P$ and unitary operator $U$. The unitary operator encodes the data $\mathbf{x}$ into the quantum state and includes several tunable parameters $\bm{\theta}$. These parameters are optimized via gradient descent [93, 94, 95, 96] using an appropriate loss function. Typically, $U$ can be divided into encoding and parametrized operators, denoted as $U_{\text{ENC}}$ and $U_{\text{PARAM}}$, respectively.

The QCGNN architecture comprises two qubit registers: an index register (IR) and a neural network register (NR), with $n_{I}$ and $n_{Q}$ qubits, respectively. Fig. 1 illustrates an example of a QCGNN circuit with $n_{I}=2$ and $n_{Q}=4$. To encode the information of an undirected and unweighted complete graph with $N$ nodes, we set $n_{I}=\lceil\log_{2}(N)\rceil$. For graphs with different number of nodes, $n_{I}$ could be different, e.g., a 4-particle jet needs $n_{I}=2$ qubits, while a 5-particle jet requires $n_{I}=3$. We will show that the output of QCGNN justifies that QCGNN can handle variable size of graphs. The quantum state in the IR is initialized to uniform basis states through a uniform state oracle (USO) and evolves as

$$\ket{\psi_{0}}=\frac{1}{\sqrt{N}}\sum^{N-1}_{i=0}\ket{i}\ket{0}^{\otimes n_{Q}},$$

(2)

where $\ket{0}^{\otimes n_{Q}}$ is the initial quantum state in the NR with all qubits in $\ket{0}$ state. The decimal representation, e.g., $\ket{0}$, $\ket{1}$, $\ket{2}$, $\ket{3}$, is used here, equivalent to the binary representation $\ket{00}$, $\ket{01}$, $\ket{10}$, $\ket{11}$. The context should clarify whether the decimal or binary representation is being used. If $N$ is a power of 2, i.e., $n_{I}=\log_{2}N$, the USO can be constructed using Hadamard gates. Otherwise, other methods described in [97, 98] can be employed.

The node features $\mathbf{x}_{i}$ of the $i$-th node are encoded through a series of unitary operators $U_{\text{ENC}}(\mathbf{x}_{i})$, controlled by the corresponding basis state in IR, where the control condition is the binary representation of $i$ with $n_{I}$ digits, denoted as $C^{n_{I}}_{\text{bin}(i)}\Bigl{(}U_{\text{ENC}}(\mathbf{x}_{i})\Bigr{)}$. This controlled operator acts on the quantum state as follows:

$$\begin{split}C^{n_{I}}_{\text{bin}(i)}&\Bigl{(}U_{\text{ENC}}(\mathbf{x}_{i})\Bigr{)}\Bigl{[}\ket{j}\ket{0}^{\otimes n_{Q}}\Bigr{]}=\\ &\delta_{ij}\ket{j}U_{\text{ENC}}(\mathbf{x}_{i})\ket{0}^{\otimes n_{Q}}+(1-\delta_{ij})\ket{j}\ket{0}^{\otimes n_{Q}},\end{split}$$

(3)

where $\delta_{ij}$ is the Kronecker delta. The controlled operator $C^{n_{I}}_{\text{bin}(i)}\Bigl{(}U_{\text{ENC}}(\mathbf{x}_{i})\Bigr{)}$ on the left-hand side acts on both the IR and NR, while $U_{\text{ENC}}(\mathbf{x}_{i})$ on the right-hand side acts only on the NR. After encoding the information of all $N$ nodes, the quantum state evolves as

$$\begin{split}\ket{\psi_{0}}&\rightarrow\ket{\psi_{1}}=\Biggl{[}\bigotimes^{N-1}_{i=0}C^{n_{I}}_{\text{bin}(i)}\Bigl{(}U_{\text{ENC}}(\mathbf{x}_{i})\Bigr{)}\Biggr{]}\ket{\psi_{0}}\\ &=\frac{1}{\sqrt{N}}\sum^{N-1}_{j=0}\Biggl{[}\bigotimes^{N-1}_{i=0}C^{n_{I}}_{\text{bin}(i)}\Bigl{(}U_{\text{ENC}}(\mathbf{x}_{i})\Bigr{)}\Biggr{]}\Bigl{(}\ket{j}\ket{0}^{\otimes n_{Q}}\Bigr{)}\\ &=\frac{1}{\sqrt{N}}\sum^{N-1}_{i=0}\ket{i}U_{\text{ENC}}(\mathbf{x}_{i})\ket{0}^{\otimes n_{Q}},\\ \end{split}$$

(4)

where it is used that $\ket{j}$ must be acted upon by one of the $C^{n_{I}}_{\text{bin}(i)}\Bigl{(}U_{\text{ENC}}(\mathbf{x}_{i})\Bigr{)}$ as $i$ ranges from $0$ to $N-1$. For example, in Fig. 1, where $N=3$, $n_{Q}=4$, and $n_{I}=\lceil\log_{2}N\rceil=2$, the quantum state $\ket{\psi_{1}}$ can be expressed as

$$\begin{split}\ket{\psi^{N=3}_{1}}&=\frac{1}{\sqrt{3}}\biggr{[}\ket{0}U_{\text{ENC}}(\mathbf{x}_{0})\ket{0000}\\ &+\ket{1}U_{\text{ENC}}(\mathbf{x}_{1})\ket{0000}+\ket{2}U_{\text{ENC}}(\mathbf{x}_{2})\ket{0000}\biggr{]},\end{split}$$

(5)

where $\ket{0000}$ is equivalent to $\ket{0}^{\otimes 4}$. At first glance, Eq. 4 and Eq. 5 appear to depend on the node ordering. However, as we will demonstrate, by selecting the appropriate observables, the final output of QCGNN is permutation invariant.

Next, a series of parametrized unitary operators $U_{\text{PARAM}}(\bm{\theta})$ with tunable parameters $\bm{\theta}$ are applied to the NR, evolving the quantum state to

$$\begin{split}\ket{\psi_{1}}\rightarrow\ket{\psi_{2}}&=U_{\text{PARAM}}(\bm{\theta})\ket{\psi_{1}}\\ &=\frac{1}{\sqrt{N}}\sum^{N-1}_{i=0}\ket{i}U_{\text{PARAM}}(\bm{\theta})U_{\text{ENC}}(\mathbf{x}_{i})\ket{0}^{\otimes n_{Q}}.\end{split}$$

(6)

To increase the expressive power of VQCs, the data re-uploading technique [91] can be employed. The idea of data re-uploading is to encode the data multiple times, allowing the quantum state to interact with the data in a more complex manner. The data re-uploading technique can be implemented by alternately applying the encoding and parametrized operators several times, but with different parameters. As the encoding operators in $\ket{\psi_{2}}$ correspond to the node indices $\ket{i}$ in the IR, the data re-uploading technique can be implemented straightforwardly without making each particle information entangled. After re-uploading $R$ times, the final quantum state $\ket{\psi}$ evolves as

$$\ket{\psi}=\frac{1}{\sqrt{N}}\sum^{N-1}_{i=0}\ket{i}\Biggl{[}\bigotimes^{R}_{r=1}U_{\text{PARAM}}(\bm{\theta}^{(r)})U_{\text{ENC}}(\mathbf{x}_{i})\Biggr{]}\ket{0}^{\otimes n_{Q}}\\ ,$$

(7)

where the parameters $\bm{\theta}$ may differ across each $U_{\text{PARAM}}$ operator, denoted by the superscript $r$. We denote the quantum state of the NR as

$$\ket{\mathbf{x}_{i},\bm{\theta}}=\Biggl{[}\bigotimes^{R}_{r=1}U_{\text{PARAM}}(\bm{\theta}^{(r)})U_{\text{ENC}}(\mathbf{x}_{i})\Biggr{]}\ket{0}^{\otimes n_{Q}}.$$

Given a Pauli string observable $P$ applied to the NR, the expectation value of the measurement is given by

$$\begin{split}\braket{\psi}{P}{\psi}&=\frac{1}{N}\sum^{N-1}_{i=0}\sum^{N-1}_{j=0}\braket{i}{j}\bra{\mathbf{x}_{i},\bm{\theta}}P\ket{\mathbf{x}_{j},\bm{\theta}}\\ &=\frac{1}{N}\sum^{N-1}_{i=0}\sum^{N-1}_{j=0}\delta_{ij}\bra{\mathbf{x}_{i},\bm{\theta}}P\ket{\mathbf{x}_{j},\bm{\theta}}\\ &=\frac{1}{N}\sum^{N-1}_{i=0}\bra{\mathbf{x}_{i},\bm{\theta}}P\ket{\mathbf{x}_{i},\bm{\theta}},\end{split}$$

(8)

which results in a sum over each node. Now, consider an additional $2^{n_{I}}\times 2^{n_{I}}$ Hermitian matrix $J$ filled with ones, used as an observable of the IR. Observing that

$$\begin{split}\braket{\psi}{J\otimes P}{\psi}&=\frac{1}{N}\sum^{N-1}_{i=0}\sum^{N-1}_{j=0}\bra{i}J\ket{j}\bra{\mathbf{x}_{i},\bm{\theta}}P\ket{\mathbf{x}_{j},\bm{\theta}}\\ &=\frac{1}{N}\sum^{N-1}_{i=0}\sum^{N-1}_{j=0}\bra{\mathbf{x}_{i},\bm{\theta}}P\ket{\mathbf{x}_{j},\bm{\theta}}\\ &=\frac{1}{N}\sum^{N-1}_{i=0}\sum^{N-1}_{j=0}h(\mathbf{x}_{i},\mathbf{x}_{j};P),\end{split}$$

(9)

where $\bra{i}J\ket{j}=J_{ij}=1$, and we define

$$h(\mathbf{x}_{i},\mathbf{x}_{j};P)=\frac{1}{2}\Bigl{[}\bra{\mathbf{x}_{i},\bm{\theta}}P\ket{\mathbf{x}_{j},\bm{\theta}}+\bra{\mathbf{x}_{j},\bm{\theta}}P\ket{\mathbf{x}_{i},\bm{\theta}}\Bigr{]},$$

which is symmetric, i.e., $h(\mathbf{x}_{i},\mathbf{x}_{j};P)=h(\mathbf{x}_{j},\mathbf{x}_{i};P)$. Notice that Eq. 9 computes the average value over all possible pairs, leading to the permutation invariance of the final output. In practice, the observable $J$ can be decomposed as

$$J=\bigotimes^{n_{I}}_{q=1}\left[\begin{matrix}1&1\\ 1&1\end{matrix}\right]=\bigotimes^{n_{I}}_{q=1}(I^{\text{IR}}_{q}+X^{\text{IR}}_{q}),$$

(10)

where $X^{\text{IR}}_{q}$ refers to the Pauli-X observable of the $q$-th qubit in IR, and $I^{\text{IR}}_{q}$ is the $2\times 2$ identity matrix. The expansion of Eq. 10 on the right-hand side yields $2^{n_{I}}\approx N$ different combinations of Pauli string observables. The value of Eq. 9 can be obtained by summing the expectation values of all combinations of Pauli strings in Eq. 10, corresponding to the SUM operation in the classical MPGNN’s aggregation function. In certain cases, one might wish to exclude contributions from nodes themselves, these values can be simply removed by considering

$$J\longrightarrow J-\bigotimes^{n_{I}}_{q=1}I^{\text{IR}}_{q},$$

which is equivalent to subtracting the value of Eq. 8 from Eq. 9.

In this subsection, we examine the computational complexity of both MPGNN and QCGNN in the context of complete graphs. We focus on the simplest case where both models map a set of $d$-dimensional features to a scalar, i.e., $\mathbb{R}^{d}\rightarrow\mathbb{R}$. For instance, the MPGNN might employ a feed-forward neural network terminating in a single neuron, while the QCGNN could utilize one qubit in NR, measured in the $Z$ basis ($n_{Q}=1$ and $P=Z$). Additionally, we assume that the computational cost of obtaining a scalar output in classical models is roughly equivalent to the cost of measuring a Pauli string observable in quantum models.

To compute all pairwise information, MPGNN requires $O(N^{2})$ computations, with each pair passing through the neural network once. In contrast, due to the quantum parallelism inherent in QCGNN, $U_{\text{PARAM}}$ can process all pairs of nodes simultaneously. To aggregate the final result, QCGNN requires only measuring on $O(2^{n_{I}})\approx O(N)$ Pauli string observables, as indicated by Eq. 10. This suggests that QCGNN could offer a polynomial speedup over MPGNN. However, in the case of QCGNN, additional costs associated with multi-controlled operators and the USO should be taken into account. Although various methods exist for decomposing multi-controlled operators, such as those discussed in [99, 100, 101], for simplicity, we adhere to a basic approach outlined in [102], which requires an additional $O(n_{I})\approx O(\log_{2}N)$ ancilla qubits and Toffoli gates. Based on the results in [97, 98], preparing a uniform quantum state necessitates $O(\log_{2}N)$ gates. If the parametrized operators sufficiently deep, these additional costs may become negligible, enabling QCGNN to achieve an $O(N)$ speedup over MPGNN.

Certain traditional VQC ansatz, such as quantum kernel methods [103, 104], also share similarities with QCGNN. Quantum kernel methods compute the kernel function $k(\mathbf{x}_{i},\mathbf{x}_{j})=|\braket{\mathbf{x}_{i},\bm{\theta}}{\mathbf{x}_{j},\bm{\theta}}|^{2}$ and are typically constructed using $U_{\text{ENC}}$ and $U_{\text{PARAM}}$. These methods have a computational complexity of $O(N^{2})$, as they still require computing each pair individually. Again, if $U_{\text{PARAM}}$ is sufficiently deep, the additional costs from data encoding may be negligible, allowing QCGNN to achieve an $O(N)$ speedup over quantum kernel methods. This advantage arises because QCGNN computes $O(N^{2})$ pairwise information simultaneously, with only $O(2^{n_{I}})\approx O(N)$ additional measurement costs, as given by Eq. 10.

QCGNN can also be extended to weighted graphs, but the additional cost might render it impractical. Consider a simple case, where an undirected, weighted graph has an adjacency matrix $A$ that can be expressed as the outer product of a vector $\ket{w}=\sum_{i}w_{i}\ket{i}$, with edge weight $A_{ij}=w_{i}w_{j}$. Instead of initializing the IR uniformly, we initialize the quantum state as

$$\sum^{N-1}_{i=0}\frac{1}{\sqrt{N}}\ket{i}\ket{0}^{\otimes n_{Q}}\longrightarrow\sum^{N-1}_{i=0}\frac{w_{i}}{\sqrt{\braket{w}{w}}}\ket{i}\ket{0}^{\otimes n_{Q}},$$

so that the terms in Eq. 9 are modified as

$$\frac{1}{N}h(\mathbf{x}_{i},\mathbf{x}_{j};P)\rightarrow\frac{w_{i}w_{j}}{\braket{w}{w}}h(\mathbf{x}_{i},\mathbf{x}_{j};P)=\frac{A_{ij}}{\operatorname{Tr}(A)}h(\mathbf{x}_{i},\mathbf{x}_{j};P).$$

Instead of initializing with a uniform state, the quantum state $\ket{w}$ can be initialized using AMPLITUDE EMBEDDING, where the information of the weights $w_{i}$ is embedded in the amplitude of $\ket{i}$.

To generalize to directed and weighted graphs, note that any matrix can be decomposed into symmetric and skew-symmetric matrices. Since both types are normal matrices, they are diagonalizable according to the spectral theorem. One can apply the method described above for each eigenbasis individually and multiply by a factor proportional to the corresponding eigenvalue. However, the additional computational cost associated with diagonalizing matrices and AMPLITUDE EMBEDDING might be substantial, potentially negating the advantages of QCGNN.

For practical applications, we primarily consider the use of QCGNN for undirected, unweighted, and complete graphs. The added complexity of handling weighted, directed, and incomplete graphs may diminish the computational benefits of QCGNN, making it less feasible for real-world applications without further optimizations.

We demonstrate the feasibility of QCGNN using two publicly available Monte Carlo simulated datasets³³3In the first published version, we used the dataset generated by ourselves using [105, 106, 107, 108, 109]. For the revised version, we switched to other existing public datasets since the number of data is much more sufficient. for jet discrimination: the Top dataset [82] and the JetNet dataset [84]. The jets in both datasets are clustered using the anti-$k_{T}$ algorithm [110, 111] with a distance parameter $R=0.8$.

The Top dataset [82] is used for binary classification, distinguishing signal jets from top quarks (Top) and background jets from mixed quark-gluon interactions (QCD). The transverse momentum of the jets is in the range $[550,650]$ GeV. The dataset is divided into 1.2 million training samples, 400 thousand validation samples, and 400 thousand testing samples. Further details of the Top dataset can be found in [83].

The JetNet dataset [84] is used for multi-class classification, with jets originating from gluons (g), top quarks (t), light quarks (q), $W$ bosons (w), and $Z$ bosons (z). Each class of jet has a transverse momentum of approximately 1 TeV, with around 170 thousand samples. For each jet event, only the top 30 particles with the highest transverse momentum are retained if the number of particles exceeds 30. Further details of the JetNet dataset can be found in [85]

In our approach, each jet is represented as a complete graph. Each node corresponds to a particle in the jet, with node features related to particle flow information. For the $i$-th particle, the input features $\mathbf{x}^{(0)}_{i}$ include the transverse momentum fraction $z_{i}={p_{T}}_{i}/{p_{T}}_{jet}$, the relative pseudorapidity $\Delta\eta_{i}=\eta_{i}-\eta_{jet}$, and the relative azimuthal angle $\Delta\phi_{i}=\phi_{i}-\phi_{jet}$. For QCGNN, these input features are further preprocessed as follows:

$$\begin{split}z_{i}&\longrightarrow\tan^{-1}(z_{i}),\\ \Delta\eta_{i}&\longrightarrow\frac{\pi}{2}\frac{\Delta\eta_{i}}{R},\\ \Delta\phi_{i}&\longrightarrow\frac{\pi}{2}\frac{\Delta\phi_{i}}{R},\end{split}$$

(11)

since rotation gates are used for data encoding (see Sec. III.2). Note that the indices of the particles are arbitrary due to the use of permutation-invariant models for graphs.

The classical model for benchmarking is based on MPGNN from Eq. 1, with the aggregation function chosen to be SUM. The function $\Phi$ is implemented as a feed-forward neural network consisting of linear layers and the ReLU activation functions [113], while $\gamma$ is simply the summation of $\Phi$ only, i.e., $\mathbf{x}^{(1)}_{i}=\sum_{j\neq i}\Phi(\mathbf{x}^{(0)}_{i},\mathbf{x}^{(0)}_{j})$, where $j$ ranges from 0 to $N-1$ except $i$. The input to $\Phi$ is simply the concatenation of $\mathbf{x}^{(0)}_{i}$ and $\mathbf{x}^{(0)}_{j}$, requiring only 6 neurons in the input layer of $\Phi$. Consequently, the graph feature, denoted as $\mathbf{x}^{C}$, is computed through

$$\begin{split}\mathbf{x}^{C}&=\sum^{N-1}_{i=0}\mathbf{x}^{(1)}_{i}\\ &=\sum^{N-1}_{i=0}\sum_{j\neq i}\Phi(\mathbf{x}^{(0)}_{i},\mathbf{x}^{(0)}_{j})\\ &=\sum^{N-1}_{i=0}\sum^{N-1}_{j=0}\Phi(\mathbf{x}^{(0)}_{i},\mathbf{x}^{(0)}_{j})-\sum^{N-1}_{i=0}\Phi(\mathbf{x}^{(0)}_{i},\mathbf{x}^{(0)}_{i})\end{split}$$

(12)

The structure of MPGNN is similar to the Particle Flow Network (PFN) in [34], with the distinction that PFN calculates $\Phi(\mathbf{x}^{(0)}_{i})$ for each particle, whereas MPGNN calculates the pairwise information $\Phi(\mathbf{x}^{(0)}_{i},\mathbf{x}^{(0)}_{j})$ between particles.

The quantum model is based on QCGNN, which consists of encoding operators and parametrized operators. The data re-uploading technique [91] is employed 2 times before the final measurements (indicated by the dashed box in Fig. 1 with $R=2$). For simplicity, we use single-angle rotation gates, defined as

$$\begin{split}R_{x}(\theta)&=\begin{bmatrix}\cos(\theta/2)&-i\sin(\theta/2)\\ -i\sin(\theta/2)&\cos(\theta/2)\end{bmatrix},\\ R_{y}(\theta)&=\begin{bmatrix}\cos(\theta/2)&-\sin(\theta/2)\\ \sin(\theta/2)&\cos(\theta/2)\end{bmatrix},\\ R_{z}(\theta)&=\begin{bmatrix}e^{-i\theta/2}&0\\ 0&e^{i\theta/2}\end{bmatrix},\\ \end{split}$$

(13)

and triple-angle rotation gates, defined as

$$\begin{split}R(\alpha,\beta,\gamma)&=R_{z}(\gamma)R_{y}(\beta)R_{z}(\alpha)\\ &=\begin{bmatrix}e^{-i\frac{\alpha+\gamma}{2}}\cos(\frac{\beta}{2})&-e^{-i\frac{\alpha-\gamma}{2}}\sin(\frac{\beta}{2})\\ e^{-i\frac{\alpha-\gamma}{2}}\sin(\frac{\beta}{2})&e^{-i\frac{\alpha+\gamma}{2}}\cos(\frac{\beta}{2})\end{bmatrix}\end{split}$$

(14)

to encode the particle flow information. The parametrized operators are constructed with strongly entangling layers [112] using rotation gates and CNOT gates. The ansatz for encoding $i$-th particle features and the strongly entangling layers are shown in Fig. 2 and Fig. 3 respectively. The $q$-th component of the graph feature $\mathbf{x}^{Q}\in\mathbb{R}^{n_{Q}}$ is computed by

$$\begin{split}\mathbf{x}^{Q}_{q}&=N\Bigl{[}\braket{\psi}{J\otimes Z^{\text{NR}}_{q}}{\psi}-\braket{\psi}{Z^{\text{NR}}_{q}}{\psi}\Bigr{]}\\ &=\sum^{N-1}_{i=0}\sum^{N-1}_{j=0}h(\mathbf{x}^{(0)}_{i},\mathbf{x}^{(0)}_{j};Z^{\text{NR}}_{q})-\sum^{N-1}_{i=0}h(\mathbf{x}^{(0)}_{i},\mathbf{x}^{(0)}_{i};Z^{\text{NR}}_{q})\\ \end{split}$$

(15)

where the observable $Z^{\text{NR}}_{q}$ refers to the Pauli-Z measurement of the $q$-th qubit in NR only, and the summations are computed through Eq. 8 and Eq. 9. To be clarified, the subscript $i$ of $\mathbf{x}^{(0)}_{i}$ corresponds to the $i$-th node, while the subscript $q$ of $\mathbf{x}^{Q}_{q}$ refers to the $q$-th component of QCGNN output. Note that all qubits in NR can be measured simultaneously, but the measurement output from other qubits is ignored when calculating the expectation value over $Z^{\text{NR}}_{q}$. This setup can be thought of as a classical feed-forward neural network with $n_{Q}$ neurons in the output layer. Notice how Eq. 15 resembles Eq. 12, indicating the permutation invariance of the final output.

Eventually, both $\mathbf{x}^{C}$ and $\mathbf{x}^{Q}$ are followed by another feed-forward neural network consisting of linear layers and ReLU activation functions. The full setup of the classical and quantum models is depicted in Fig. 4. For binary classification using the Top dataset, the output layer has a single neuron followed by a Sigmoid function and is trained with binary cross-entropy loss. For multi-class classification using the JetNet dataset, the output layer has five neurons followed by a Softmax function and is trained with multi-class cross-entropy loss.

The complete training process was conducted with 5 different random seeds, each for 30 epochs. The Top and the JetNet datasets comprise 2 and 5 classes, respectively. For each class, we selected 25,000 training samples, 2,500 validation samples, and 2,500 testing samples. This limited data selection is due to the extensive training time required for the QCGNN, as discussed in Sec. III.4. For each random seed, the data were randomly sampled from the original dataset. To balance between demonstrating the training performance and computational demands, particles with transverse momentum less than 2.5% of ${p_{T}}_{jet}$, i.e., $z_{i}<0.025$, were discarded, so that the majority of the distribution of the number of particles per jet is between 4 and 16. The histograms of the number of particles per jet for the original and preprocessed Top and JetNet datasets are shown in Fig. 5. To mitigate the extensive training time associated with simulating QML, events with fewer than 4 or more than 16 particles were discarded. These choices strike a balance between performance and the amount of training data, as discussed in Appendix B.

Due to limited computational resources, the number of qubits $n_{Q}$ in NR was tested with $n_{Q}=3$ and $n_{Q}=6$, with the number of strongly entangling layers in each parametrized operator ($U_{\text{PARAM}}$) set to $n_{Q}/3$. Given that the maximum number of particles in a jet is 16, we used $n_{I}=\lceil\log_{2}16\rceil=4$ qubits for IR. The number of hidden neurons $n_{M}$ in the MPGNN was set to 3 or 6 for comparison with the QCGNN, ensuring that both models have a comparable number of parameters, as discussed in Appendix C.

We also evaluated the performance of classical state-of-the-art models, including the Particle Flow Network (PFN) [34], Particle Net (PNet) [35], and Particle Transformer (ParT) [36]. The structure and hyperparameters of PFN, PNet, and ParT were configured according to their respective original publications. Notably, we excluded mass information from the interaction matrix of ParT, as only particle flow information $(p_{T},\Delta\eta,\Delta\phi)$ was used. The $k$-nearest neighbor method used in the original PNet was configured with $k=3$, given that the minimum number of particles per jet was 4.

The classical models were implemented using PyTorch [114] and PyTorch Geometric [115], while the quantum circuit of QCGNN was simulated using PennyLane [116]. The cross-entropy loss was optimized using the Adam optimizer [117] with a learning rate of $10^{-3}$ for all models. The batch size was set to 64, the maximum allowable due to memory constraints, as simulating quantum circuits requires substantial memory resources.

Unlike classical models, parameter gradients for real quantum computers cannot be computed using traditional methods such as finite difference methods. Instead, the parameter-shift rule (PSR) [93, 94, 95, 96] can be employed to calculate gradients. However, applying PSR on quantum computers necessitates extensive requests and long queue times with actual quantum devices. Furthermore, the noise in current quantum computers is insufficiently low to enable stable training of quantum neural networks, often resulting in training failures.

To circumvent these issues during the NISQ era [92], we trained the QCGNNs on classical computers using PennyLane [116] quantum circuit simulators with zero noise. Nonetheless, simulating quantum circuits is highly time-consuming, even for a few qubits. Although PennyLane supports QML on GPUs, speed improvements over CPUs are significant only with many qubits, typically more than 20 qubits⁴⁴4The benchmark of quantum simulation on GPU can be found in PennyLane’s blog: ”Lightning-fast simulations with PennyLane and the NVIDIA cuQuantum SDK”. For this study, we used CPUs to train the QCGNNs. Training a 10-qubit ($n_{I}=4$ and $n_{Q}=6$) QCGNN with 10,000 samples and a batch size of 64 takes approximately 1,000 seconds per epoch. Consequently, the training over 5 random seeds for 30 epochs required nearly a month.

The performance of the classical and quantum models on the Top and the JetNet datasets is summarized in Table. 1. The inference scores of the MPGNN and QCGNN are comparable when the number of parameters is in roughly the same order of amount. We anticipate that QCGNN has the potential to achieve performance on par with MPGNN as the number of qubits increases. When training with a smaller number of parameters on the multi-class classification JetNet dataset, we observe that QCGNN is more stable than MPGNN, with the latter exhibiting a larger standard deviation.

Among the state-of-the-art models, MPGNN with $n_{M}=64$ has higher inference scores compared to others. This is likely due to the information of jets being lost during the data preprocessing, where only 4 to 16 particles per jet are utilized. When training with the full information from the original jet dataset, i.e., without discarding information from soft particles, other state-of-the-art models can compete with MPGNN, even better. Details of the state-of-the-art models trained on the original jet dataset are provided in Appendix B.

Although implementing full training on quantum computers is impractical in the NISQ era, we can still evaluate the performance of the pre-trained QCGNN on IBMQ real devices [118]. To minimize the noise effects caused by real quantum gates, we select events with only four particles from the Top dataset, i.e., using $n_{I}=2$ qubits in IR, thereby reducing the number of gates required for initial state preparation and data encoding. In this setup, the USO can be efficiently implemented using Hadamard gates for each qubit in IR. On IBMQ real devices, only 1-qubit and 2-qubit gates are available, and the multi-controlled gates used in data encoding are decomposed using methods described in [102].

We selected ibm_brussels with 1024 shots to test the performance of QCGNN on an IBMQ real device. However, the quantum computers in the NISQ era are currently too noisy to yield usable results. For binary classification, the inference of QCGNN on ibm_brussels results in approximately $0.5$ AUC and $0.5$ accuracy, which equates to random guessing. To assess how noise affects the performance of QCGNN, we perform an extrapolation over noise using PennyLane simulators, with the results shown in Fig. 6. We simulate quantum noise, including depolarizing error and amplitude damping, occurring after each quantum operation with a certain probability. As indicated in Fig. 6, the noise probability must be reduced to below $10^{-3}$ to achieve reliable results⁵⁵5The noise probability here corresponds to the simulated noise in the quantum circuit simulation..

To validate the time complexity analysis discussed in Sec. II.3, we initialized untrained QCGNNs and executed them on various IBMQ backends, including ibm_nazca and ibm_strasbourg, with $n_{Q}=100$ and 1024 shots. We set the number of nodes to 2, 4, and 8, such that only Hadamard gates are required for the initial state preparation. To determine the quantum gate runtime for encoding and parametrized operators, we first ran QCGNN without any operators to measure the runtime $T_{0}$ for quantum state initialization and measurement. We then applied encoding operators with 10 times of re-uploading to obtain the runtime $T_{1}$. Finally, we applied parametrized operators, constructed with 10 strongly entangling layers and 10 times of re-uploading (resulting in 100 strongly entangling layers in total), to measure the runtime $T_{2}$. Each runtime measurement was averaged over 10 executions. The runtime of encoding operators was computed as follows:

$$T_{\text{ENC}}=\frac{T_{1}-T_{0}}{10},$$

and the runtime of each strongly entangling layer in parametrized operators is calculated by

$$T_{\text{PARAM}}=\frac{T_{2}-T_{1}}{100}.$$

The results presented in Table. 2 indicate that the runtime of encoding operators scales approximately linearly with the number of particles per jet, while the runtime of parametrized operators remains approximately constant, as expected. As discussed in Sec. II.3, when the parametrized operators are sufficiently deep, the runtime will be dominated by these operators, making the additional computational cost associated with data encoding negligible.