A quantum hardware-induced graph kernel based on Gaussian Boson Sampling

TL;DR

This paper introduces a novel quantum hardware-based graph kernel using Gaussian Boson Sampling, demonstrating its effectiveness for graph similarity and providing theoretical insights into its performance, with promising applications in quantum computing.

Contribution

It presents a new quantum hardware-inspired graph kernel leveraging Gaussian Boson Sampling, linking quantum sampling distributions to graph matchings for improved similarity measures.

Findings

Kernel performs well on benchmark datasets

Links Gaussian Boson Sampling to graph matchings

Provides theoretical motivation for kernel effectiveness

Abstract

A device called a 'Gaussian Boson Sampler' has initially been proposed as a near-term demonstration of classically intractable quantum computation. As recently shown, it can also be used to decide whether two graphs are isomorphic. Based on these results we construct a feature map and graph similarity measure or 'graph kernel' using samples from the device. We show that the kernel performs well compared to standard graph kernels on typical benchmark datasets, and provide a theoretical motivation for this success, linking the distribution of a Gaussian Boson Sampler to the number of matchings in subgraphs. Our results contribute to a new way of thinking about kernels as a (quantum) hardware-efficient feature mapping, and lead to an interesting application for near-term quantum computing.

Figures

Tables

Table 1. Table 2: Mean test accuracy of the Support Vector Machine with different datasets and different graph kernels, with the standard deviation between 10 10 10 repetitions of the double cross-validation. GS, RW, and SM are three standard classical graph kernels described in the text. GBS refers to the postprocessing strategy of associating orbit probabilities with features, while GBS + summarises some orbits to meta-orbits (see text). ∗ Runtime > 20 absent 20 >20 days.

Dataset	GBS ($d_{0.0}$)	GBS ($d_{0.25}$)	GBS+ ($d_{0.0}$)	GBS+ ($d_{0.25}$)	GS	RW	SM
AIDS	$99.60\pm 0.05$	$99.62\pm 0.03$	$99.58\pm 00.06$	$99.61\pm 0.05$	$98.44\pm 0.09$	$56.95\pm 7.99$	$79.20\pm 0.68$
BZR_MD	$62.73\pm 0.71$	$62.13\pm 1.44$	$62.01\pm 1.43$	$63.16\pm 2.11$	$60.60\pm 1.77$	$49.88\pm 3.74$	$61.90\pm 1.21$
COX2_MD	$44.98\pm 1.80$	$50.11\pm 0.97$	$57.84\pm 4.04$	$57.89\pm 2.62$	$55.04\pm 3.33$	$57.72\pm 3.26$	$66.94\pm 1.22$
ENZYMES	$22.29\pm 1.60$	$28.01\pm 1.83$	$25.72\pm 2.60$	$40.42\pm 2.02$	$35.87\pm 2.19$	$21.13\pm 1.91$	$36.70\pm 2.83$
ER_MD	$70.36\pm 0.78$	$70.41\pm 0.47$	$71.01\pm 1.26$	$71.05\pm 0.83$	$65.65\pm 1.06$	$68.75\pm 0.53$	$68.21\pm 0.99$
FINGERPRINT	$65.42\pm 0.49$	$\text{65.85}\pm 0.36$	$66.19\pm 00.84$	$66.26\pm 4.29$	$64.10\pm 1.52$	$47.69\pm 0.21$	$47.14\pm 0.62$
IMDB-BIN	$64.09\pm 0.34$	$\text{68.71}\pm 0.59$	$68.14\pm 0.71$	$67.60\pm 0.75$	$68.37\pm 0.62$	$66.38\pm 0.21$	out of time∗
MUTAG	$86.41\pm 0.33$	$85.58\pm 0.59$	$85.64\pm 0.78$	$84.46\pm 0.44$	$81.08\pm 0.93$	$83.02\pm 1.08$	$83.14\pm 0.24$
NCI1	$63.61\pm 0.00$	$62.79\pm 0.00$	$63.59\pm 0.17$	$63.11\pm 0.93$	$49.96\pm 3.27$	$52.36\pm 2.63$	$51.36\pm 1.88$
PROTEINS	$66.88\pm 0.22$	$66.14\pm 0.48$	$65.73\pm 0.69$	$66.16\pm 0.76$	$65.91\pm 1.29$	$56.27\pm 1.23$	$63.03\pm 0.84$
PTC_FM	$53.84\pm 0.96$	$52.45\pm 1.78$	$59.14\pm 1.72$	$56.25\pm 2.04$	$59.48\pm 1.95$	$51.97\pm 2.68$	$54.92\pm 2.94$

Equations

$$\tilde{A}=c\left(\begin{array}{cc}A& 0\\ 0& A\end{array}\right)=c(A\oplus A),$$

$$\sigma =Q-\text{mathbbm}1/2,\text{ with }Q=(\text{mathbbm}1-X\tilde{A}{)}^{-1},X=\left(\begin{array}{cc}0& \text{mathbbm}1\\ \text{mathbbm}1& 0\end{array}\right).$$

$$p(\mathbf{n})=\frac{1}{\sqrt{\det (Q)}\mathbf{n}!}{\mathrm{haf}}^2(A_{\mathbf{n}}),$$

$$\mathrm{haf}(C)=\sum _{\pi \in P_N^{\{2\}}}\prod _{(u,v)\in \pi }{C}_{u,v}.$$

$$P_4^{\{2\}}=\{(1,2),(3,4)\},\{(1,3),(2,4)\},\{(1,4),(2,3)\},$$

$$P_4^{\{1,2\}}=\{$$

$$\{(2,3),(1),(4)\},\{(2,4),(1),(3)\},\{(3,4),(1),(2)\},$$

$$\{(1,2),(3,4)\},\{(1,3),(2,4)\},\{(1,4),(2,3)\},$$

$$\{(1),(2),(3),(4)\}$$

$$\mathrm{haf}({\tilde{A}}_{\mathbf{n}-\{i_1,\dots ,i_{2n}\}})=\mathrm{haf}(A_{\mathbf{n}-\{j_1,\dots ,j_s\}})\mathrm{haf}(A_{\mathbf{n}-\{k_1,\dots ,k_{s^{\prime }}\}}).$$

$$p(\mathbf{n})=\alpha {\left[\sum _{n=0}^M\sum _{\{i_1..i_n\}\subseteq {\mathcal{I}}_M}{b}_{i_1}\dots b_{i_n}\mathrm{haf}(A_{\mathbf{n}-\{i_1..i_n\}})\right]}^2$$

$$\alpha =\frac{e^{-\frac{1}{2}{\mathbf{d}}^{†}{Q}^{-1}\mathbf{d}}}{\sqrt{\det (Q)}\mathbf{n}!},$$

$$p(O_{\mathbf{n}}):=\sum _{\mathbf{n}\in O_{\mathbf{n}}}p(\mathbf{n}).$$

$$\Delta_{s}=\big{\{}\mathbf{n}:\sum_{i}n_{i}=|\mathbf{n}|,(\forall i)(n_{i}\leq s),(\forall\mathbf{n}\exists n_{i}\in\mathbf{n})(n_{i}=s)\big{\}}.$$

$$p({\mathcal{M}}_{|\mathbf{n}|,{\Delta }_s}):=\sum _{\mathbf{n}\in {\Delta }_s}p(O_{\mathbf{n}}).$$

$$S=\lceil \frac{2(\log (2)D+\log (\frac{1}{\delta }))}{{ϵ}^2}\rceil$$

$$f_i=p(O_{\mathbf{n}}^i),\text{ or }{f}_i=p({\mathcal{M}}_{|\mathbf{n}|,{\Delta }_s}^i),$$

$${\kappa }_{\mathrm{lin}}(G,G^{\prime })$$

$${\kappa }_{\mathrm{rbf}}(G,G^{\prime })$$

$$\mu (G)=\sum _{r=0}^{\lceil M/2\rceil }(-1{)}^{r}m(G,r)x^{M-2r}.$$

$$g(G,r)=\sum _{\mathbf{n}\in O_{[1,\dots ,1,0,\dots ]}}{\mathrm{haf}}^2(A_{\mathbf{n}}).$$

$$f_g\propto \sum _{\mathbf{n}\in O_{[1,..,1,0,\dots ]}}\text{mathbbm}{1}_{g\cong G_{\mathbf{n}}},$$

$$f_i=f_{{\mathbf{n}}_i^{\ast }}\propto \sum _{\mathbf{n}\in O_{{\mathbf{n}}^{\ast }}}{\left(\sum _{g\in {\mathcal{P}}^{|\mathbf{n}|}}\text{mathbbm}{1}_{g\cong G_{\mathbf{n}}}\right)}^2,$$

$$\sigma (\nu )=(1-\nu )\sigma +\frac{\nu }{2}\text{mathbbm}{1}_M,$$

$$\tilde{C}=X_{2M}\big{(}\mathbbm{1}_{2M}-W^{-1}\mathrm{diag}\left[{\Lambda_{1}-1\over\Lambda_{1}\nu-1},\dots,{\Lambda_{2M}-1\over\Lambda_{2M}\nu-1}\right]W\big{)},$$

$$\mathrm{haf}(\text{mathaccentV}tilde07E{A}_{\mathbf{n}-\{i_1,\dots ,i_{2n}\}})=\mathrm{haf}(A_{\mathbf{n}-\{j_1,\dots ,j_s\}})\mathrm{haf}(A_{\mathbf{n}-\{k_1,\dots ,k_{s^{\prime }}\}}).$$

Full text

I Introduction

Measuring the similarity of two graphs for practical applications is notoriously difficult. Firstly, there are many different notions of similarity, and practical tasks crucially depend on what property of the graph is exploited in the comparison. Secondly, even the task of determining whether two graphs are exactly the same can be computationally extremely costly. This is due to the fact that a representation of a graph is not unique: Different ways of enumerating its nodes and edges can give rise to the same structure. The complexity of deciding whether two graphs are isomorphic is unknown; neither a polynomial-time algorithm nor NP-completeness proof has been discovered yet Kobler et al. (2012). Existing algorithms for graph isomorphism McKay et al. (1981) and graph similarity Ghosh et al. (2018) are efficient in practice, but are still costly for large graphs and may require exponential time for some problem instances.

In this paper we suggest the use of quantum hardware to map a graph $G$ to a feature vector which represents $G$ in Euclidean space. Standard distance measures, such as taking the inner product of two feature vectors, then result in a distance measure between graphs mitigated by the feature embedding. The quantum device we investigate is a Gaussian Boson Sampling (GBS) setup Hamilton et al. (2017); Lund et al. (2014); Kruse et al. (2018). GBS is a generalization of Boson Sampling Tillmann et al. (2013); Broome et al. (2013), which has originally been proposed as a classically intractable problem to demonstrate the power of near-term quantum hardware Aaronson and Arkhipov (2011). An optical GBS device prepares a quantum state of $M$ optical modes and counts the photons in each mode. Some of the authors have previously shown how a graph can be encoded into the quantum state of light Brádler et al. (2018), so that the photon measurement statistics give rise to a complete set of graph isomorphism invariants Brádler et al. (2018).

Here we extend this result and study the graph similarity measure derived from a GBS device for a practical application, namely for classification for machine learning. Graph-structured data plays an increasingly important role in this field, for example to predict properties of a social media network given a dataset of networks for which the properties are known. In machine learning, a similarity measure between data is called a kernel, and lots of methods for pattern recognition – such as support vector machines and Gaussian processes – are built around this concept. Mapping graphs to feature vectors or graph embeddings Zhang et al. (2018); Goyal and Ferrara (2018); Grover and Leskovec (2016) is a well known strategy, and graph kernels from explicit feature vectors Kriege et al. (2014) have been studied in detail.

The connection between kernel methods for machine learning and quantum computing has recently been made in Refs. Schuld and Killoran (2019); Havlíček et al. (2019). Any positive-definite kernel can be formally understood as the inner product of two feature vectors that represent the data points in a Hilbert space Scholkopf and Smola (2001). Hence, the Hilbert space of a quantum system can be interpreted as a feature space, in which a subroutine can compute inner products “coherently”. By using measurement samples from the quantum hardware to construct low-dimensional feature vectors that can be stored and further processed on a classical computer, we follow a different, even more minimalistic route to define a “quantum feature map”, and ultimately a quantum kernel. The advantage in using quantum hardware this way is that device performs a combinatorial computation that is very resource-intense – possibly even intractable – for classical computers. In fact, we show that the GBS feature map is related to a class of classical graph kernels which count subgraphs Shervashidze et al. (2009), but instead of only considering subgraphs of constant size, the sampling statistics reveal information on all possible subgraphs, as well as subgraphs constructed from copying nodes and their edges. The resulting features contain information about the number of $r$-matchings of the original graph. Numerical experiments reveal that graph kernels from a GBS-induced feature map can outperform classical graph kernels in classification task for standard benchmark datasets, results that can be further improved by using displaced light modes.

II Turning GBS samples into features

An optical Gaussian Boson Sampler is a device where a special quantum state (a so-called Gaussian state) is prepared by the optical squeezing of $M$ displaced light modes, followed by an interferometer of beamsplitters. Such a Gaussian state is fully described by a covariance matrix $\sigma\in\mathbb{R}^{2M\times 2M}$ as well as a displacement vector $\mathbf{d}\in\mathbb{R}^{2M}$ Weedbrook et al. (2012). Photon number resolving detectors count the photons in each mode.

In this section we describe the mathematical details of the quantum hardware-induced feature map (see also Figure 1), summarizing what has been described in Brádler et al. Brádler et al. (2018), and adding the effect of displacement as well as a further step of turning samples to feature vectors through what we will call “meta-orbits”. The scheme works for simple graphs, i.e. undirected graphs without self-loops or multiple edges. While edge weights can be treated on the same footing as unweighted edges, we leave the inclusion of categorical edge labels or node labels for future studies. Mindful of readers from fields other than quantum optics we will only highlight some important aspects of Gaussian Boson Sampling and refer to Refs. Hamilton et al. (2017); Kruse et al. (2018); Brádler et al. (2018) for more detail.

II.1 Encoding graphs into the GBS device

As outlined in Brádler et al. (2018), a quantum state prepared by a GBS device can encode a graph $G=(V,E)$ with an adjacency matrix $A$ of entries $A_{ij}$ that are one if the edge $(i,j)$ exists in $G$ and zero else. The entries of $A$ can also represent continuous “edge weights” that denote the strength of a connection. In the latter case we will speak of a “weighted adjacency matrix”.

In order to associate $A$ with the symmetric, positive definite $2M$-dimensional covariance matrix of a Gaussian state of $M$ modes, we have to construct a “doubled adjacency matrix”

[TABLE]

where the rescaling constant $c$ is chosen so that $0<c<1/s_{\mathrm{max}}$, and $s_{\mathrm{max}}$ is the maximum singular value of $A$ Brádler et al. (2018, 2018).111As long as it fulfills the above inequality, $c$ can be treated as a hyperparameter of the feature map, which may also be influenced by hardware constraints since it relates ultimately to the amount of squeezing required. For simplicity we will always rescale all adjacency matrices with a factor $1/(s^{\{G\}}_{\mathrm{max}}+10^{-8})$ where $s^{\{G\}}_{\mathrm{max}}$ is the largest singular value among all graphs in the data set under consideration. As a result we will assume that $c=1$ and $\tilde{A}=A\oplus A$ can be encoded into a GBS device. We call this the “doubled encoding strategy”.

The matrix $\tilde{A}$ can now be associated with a quantum state’s covariance matrix $\sigma$ by setting the squeezing as well as the beamsplitter angles of the interferometer so that

[TABLE]

II.2 Sampling photon counting events

After embedding $A$ via $\tilde{A}$ into the quantum state of the GBS, each measurement of the photon number resolving detectors returns a photon event $\mathbf{n}=[n_{1},\dots,n_{M}]$, with $n_{i}\in\mathbb{N}$ indicating the number of photons measured in the $i$-th mode. Assuming for now that the displacement $\mathbf{d}$ is zero, the probability of measuring a given photon counting event is

[TABLE]

where $\mathbf{n}!=n_{1}!n_{2}!\cdots n_{M}!$.

Let us go through this nontrivial equation bit by bit. The Hafnian $\mathrm{haf}()$ is a matrix operation similar to the determinant or permanent. For a general symmetric matrix $C\in\mathbb{R}^{N}\times\mathbb{R}^{N}$ it reads

[TABLE]

Here, $P_{N}^{\{2\}}$ is the set of all $N!/((N/2)!2^{N/2})$ ways to partition the index set $\{1,2,\dots,N\}$ into $N/2$ unordered pairs of size $2$, such that each index only appears in one pair. The Hafnian is zero for odd $N$. As an example, for the index set $\{1,2,3,4\}$ we have $P_{4}^{\{2\}}=\{(1,2),(3,4)\},\{(1,3),(2,4)\},\{(1,4),(2,3)\}$.

If $C$ is interpreted as an adjacency matrix containing the edges of a graph, the set $P_{N}^{\{2\}}$ contains edge-sets of all possible perfect matchings on $G$. A perfect matching is a subset of edges such that every node is covered by exactly one of the edges. The Hafnian therefore sums the products of the edge weights in all perfect matchings. If all edge weights are constant, it simply counts the number of perfect matchings in $G$ (see also Figure 2). Note that in Eq. (3) we used the fact that for real and symmetric $A$, $\text{haf}(\tilde{A})=\text{haf}(A\oplus A)=\mathrm{haf^{2}}(A)$. In other words, the doubled encoding strategy leads to a square factor which will play a profound role in the quantum feature map we are aiming to construct.

Eq. (3) does not depend on the Hafnian of the adjacency matrix $A$, but on a matrix $A_{\mathbf{n}}$. $A_{\mathbf{n}}$ contains $n_{j}$ duplicates of the $j$th row and column in $A$. If $n_{j}=0$, the $j$th row/column in $A$ does not appear in $A_{\mathbf{n}}$. Effectively, this constructs a new graph $G_{\mathbf{n}}$ from $A$ according to the following rules (see also Table 1):

If all $n_{j},j=1,\dots,M$ are one (i.e., each detector counted exactly one photon), $A_{\mathbf{n}}=A$. 2. 2.

If some $n_{j}$ are zero and others one (i.e., these detectors report no photons), $A_{\mathbf{n}}$ describes an induced subgraph $G_{\mathbf{n}}$ of $G$, in which nodes that correspond to detectors with zero count were deleted together with any edge that connected them to other nodes. 3. 3.

If some $n_{j}$ are larger than one (i.e., these detectors count more than one photon), $A_{\bm{N}}$ describes what we call an extended induced subgraph in which the corresponding nodes and all their connections are duplicated $n_{j}$ times.

In short, the probability of a photon event to be measured by the GBS device is proportional to the square of the (weighted) number of perfect matchings in a -possibly extended - induced subgraph of the graph encoded into the interferometer.

Computing the Hafnian of a general matrix is in complexity class #P, and formally reduces to the task of computing permanents Valiant (1979). If no entry in the matrix is negative, efficient approximation heuristics are known, although their success is only guaranteed under specific circumstances Barvinok (2016); Rudelson et al. (2016).

II.3 The effect of displacement

The Gaussian Boson Sampling setup underlying Eq. (3) consists of squeezing and interferometers. But a Gaussian quantum state can also be manipulated by a third operation: displacement. Displacement changes the mean of the $M$-mode Gaussian state while leaving the covariance matrix (and therefore the encoding strategy) as before. A non-zero mean changes Eq. (3) in an interesting, but non-trivial manner.

Without going into the details Kruse et al. (2018), if considering nonzero displacement, instead of summing over $P_{N}^{\{2\}}$ in Eq. (4), we have to sum over $P_{N}^{\{1,2\}}$, or the set of partitions of the index set $\{1,\dots,N\}$ into subsets of size up to $2$. For the index set $\{1,2,3,4\}$, we had

[TABLE]

which now becomes

[TABLE]

Instead of the Hafnian in Eq. (3), we therefore get a mixture of Hafnians of $A_{\mathbf{n}}$’s submatrices (stemming from the pairs) and other factors (stemming from the size-1 sets).

Assume that displacement is applied to both the $\hat{x}$ and $\hat{p}$ quadratures of each mode, described by a $2M$-dimensional displacement vector $\mathbf{d}=(d_{1},\dots,d_{M},d_{1}^{*},\dots,d_{M}^{*})^{T}$. The effect on Eq. (3) is as follows. Let $Q$ the $2M\times 2M$ matrix from Eq. (2), and $\mathbf{b}=\textbf{d}^{\dagger}Q^{-1}$. We get222To derive Eq. (5) from the analysis in Kruse et al. (2018), one uses the fact that for $\tilde{A}$ being a direct sum $A\oplus A$, the index set $i_{1},\dots,i_{2n}\in\mathcal{I}_{2M}$ can be divided into two index sets: $j_{1},\dots,j_{s}$ which contains all $s$ indices from the ‘first subspace’ (i.e., the first $M$ dimensions) of $\tilde{A}$, and $k_{1},\dots,k_{s^{\prime}}$ containing the $s^{\prime}$ indices from the ‘second subspace’, with $s+s^{\prime}=2n$. The fact that $\mathrm{haf}(A\oplus B)=\mathrm{haf}(A)\mathrm{haf}(B)$, allows us to express the Hafnian of reduced versions of $\tilde{A}_{\mathbf{n}}$ as a product of reduced versions of matrix $\tilde{A}_{\mathbf{n}}$,

$\mathrm{haf}(\tilde{A}_{\mathbf{n}-\{i_{1},\dots,i_{2n}\}})=\mathrm{haf}(A_{\mathbf{n}-\{j_{1},\dots,j_{s}\}})\mathrm{haf}(A_{\mathbf{n}-\{k_{1},\dots,k_{s^{\prime}}\}}).$

[TABLE]

with

[TABLE]

where $\mathcal{I}_{2M}$ is the index set $\{1,\dots 2M\}$. In this notation we assume $\{i_{1},\dots,i_{0}\}=\{\}$ and $b_{i_{1}}\dots b_{i_{0}}=1$. The “reduced” Hafnians of the form $A_{\mathbf{n}-\{i,j\}},A_{\mathbf{n}-\{i,j,k,l\}}\dots$ are constructed by “deleting” rows and columns $\{i,j\},\{i,j,k,l\},...$ in $\tilde{A}_{\mathbf{n}}$. The expression in the brackets of Eq. (5) is also known as a “loop Hafnian” of a matrix $\tilde{A}_{\mathbf{n}}$ that carries $b_{1},...,b_{2M}$ on its diagonal Quesada (2019).

One can see that displacement explores substructures of extended subgraphs, adding another layer of “resolution” to the photon number distribution. An important effect of displacement is that $p(\mathbf{n})$ for odd total photon numbers $|\mathbf{n}|$ is not necessarily zero any more, since the sum in Eq. (5) contains Hafnians of even-sized subgraphs.

II.4 Turning samples into features

The basic idea of how to turn samples of photon counting events into feature vectors is to associate the probability of a certain measurement result with a feature. To estimate the probability of measurement outcomes, one simply divides the number of times a result has been measured by the total number of measurements. However, if we simply used the probabilities of photon events $p(\mathbf{n})$ as features, we would face a very fast – more precisely, a doubly factorial – explosion of the number of features with the total number of photons, while almost all events become vanishingly unlikely for realistic amounts of squeezing. In practice we will truncate the total number of photons at a fixed value $k$ and discard all measurement results with $|\mathbf{n}|>k$ in the construction of the feature vector, but even then the sampling task quickly becomes unfeasible.

We therefore define the probability of certain types of photon events as features, thereby “coarse-graining” the probability distribution. As a compromise between experimental feasibility and expressive power, we consider two different coarse-graining strategies here. The first one follows Brádler et al.’s Brádler et al. (2018) suggestion to coarse-grain the distribution of photon counting events by summarizing them to sets called orbits (see Table 1). An orbit $O_{\mathbf{n}}=\{\text{perm}(\mathbf{n})\}$ contains permutations of the detection event $\mathbf{n}$. For example, $[2,1,1,0]$ is in the same orbit as $[0,1,2,1]$, but not $[2,2,0,0]$. The photon counting event $\mathbf{n}$ in the index is therefore an arbitrary “representative” of the photon counting events in an orbit. The probability of detecting a photon counting event of orbit $O_{\mathbf{n}}$ is given by the sum of the individual probabilities,

[TABLE]

The number of orbits $O_{\mathbf{n}}$ containing events of up to $k$ photons in total is equal to the number of ways that the integers of $1,...,k$ can be partitioned into a sum of at most $M$ terms. In practice we usually have $k\ll M$, in which case there are $2,4,7,12,19,30,45,67$ orbits for $k=1,\dots,8$, respectively333See also A000070 in the Online Encyclopedia of Integer Sequences, https://oeis.org/A000070.. In a real GBS setup, the energy is finite and high photon counts therefore become very unlikely. 444The energy of a Gaussian quantum state, and hence the average photon number, is determined by the squeezing and displacement operations.

The second post-processing strategy builds on top of the first, and summarizes orbits to ‘meta-orbits’ $\mathcal{M}_{|\mathbf{n}|,\Delta_{s}}$, where

[TABLE]

In words, a meta-orbit contains all orbits of $|\mathbf{n}|$ photons, which have at least one detector counting $s$ photons, but no detector counts more than $s$ photons (see also Table 1). The probability of detecting an event from a meta-orbit is given by

[TABLE]

From here on, when using meta-orbit features, we refer to the GBS as “GBS+”.

It is interesting to estimate how many samples are needed to estimate a feature vector. In Ref Shervashidze et al. (2009) we find that we can approximate a probability distribution of $D$ possible outcomes, with probability at most $\delta$ that the sum of absolute values of the errors in the empirical probabilities of the outcomes is $\epsilon$ or more, using

[TABLE]

samples. For orbits up to $k=8$ photons, there are $D=67$ features. Setting $\epsilon=0.05$ and $\delta=0.05$ and assuming a perfect GBS device, we need $39,550$ samples. Since current-day photon number resolving detectors can accumulate about $10^{5}$ samples of photon counting events per second Vaidya et al. (2019), it takes in principle only a fraction of a second for the orbit probabilities to be estimated by the physical hardware. The number of samples does not grow with the graph size, but of course the GBS device itself grows linearly in the number of nodes.

While hardware implementations of Gaussian Boson Samplers are rapidly advancing, in this paper we still resort to simulations. Sampling from photon event distributions is still a topic of active research, and to ensure that the results are not influenced by approximation errors we will use exact calculations here. This limits the scope of the experiments to graphs of the order of $25$ nodes.

II.5 Constructing a similarity measure

Summarizing the above, the feature map implemented by a GBS device maps a graph to a feature vector, $G\rightarrow\mathbf{f}\in\mathbb{R}^{D}$, where the entries $f_{i},i=1,\ldots,D$ of $\mathbf{f}$ are the probabilities of detecting certain types of photon events that we called orbits and meta-orbits,

[TABLE]

and the probability of the $i$’th (meta-)orbit is fully defined by Eqs. (6) and (7) (while ordering in the feature vector does not matter).

Assuming that the maximum number $k$ of photons we consider is smaller or equal to the number of detectors, or $k\leq M$ for all graphs, the size $D$ of the feature vector is solely determined by $k$, which is a hyperparameter of the feature map. Another hyperparameter is the displacement that can be applied to the light modes. We will assume here that the displacement applied to all modes is a constant value $d$.

Once constructed, the feature vectors can be used for various applications. In the context of machine learning, they can be directly fed into neural network classifiers. Here we are interested in constructing a similarity measure or kernel that computes the similarity between two graphs $G$ and $G^{\prime}$. A standard choice is to use the feature vectors in a ‘linear’ and ‘rbf’ kernel (with a hyperparameter $\delta$)

[TABLE]

both of which are well known to be positive semi-definite so that the results of kernel theory apply to the “GBS kernel” constructed here.

III The GBS graph features

In this section we will analyze the features of the first post-processing strategy in more detail; we discuss their intimate relation to the coefficients of a graph property called a “matching polynomial”, the relation of photon event probabilities to higher-order moments of multivariate normal distributions, the connection between the GBS and graphlet sampling kernel, and we finally discuss the devastating effect of photon loss on the features.

III.1 Single-photon features and $r$-matchings

It turns out that the probabilities of ‘single-photon’ orbits (i.e., each detector counts either zero or one photon) are related to a graph property called the “matching polynomial” of $G$ Farrell (1979); Godsil and Gutman (1981); Heilmann and Lieb (1972),

[TABLE]

The coefficients $m(G,r)$ of the matching polynomial count the number of $r$-matchings or “independent edge sets” in $G$ – sets of $r$ edges that have no vertex in common (see Figure 3). In the language of Hafnians, the $r$ matching can be written as $m(G,r)=\sum_{\mathbf{n}\in O_{[1,\dots,1,0,\dots]}}\mathrm{haf}(A_{\mathbf{n}})$ (where $[1,\dots,1,0,\dots]$ contains $2r$ single photon detections). Hence, if it were not for the square of the Hafnian in Eq. (3), the probability $p(O_{\mathbf{n}})$ of a single-photon orbit would be proportional to a $|\mathbf{n}|/2$-matching $m(G,|\mathbf{n}|/2)$ of $G$. The square gives rise to a new object

[TABLE]

Replacing $m$ with $g$ in Eq. (9) leads to a new type of polynomial $\gamma(G)$ which we call a GBS polynomial.

This definition opens up a range of interesting questions, for example whether the GBS polynomial has advantages over a standard polynomial, or how multi-photon events and displacement fits into this interpretation. We will investigate these questions in separate works.

An interesting observation for the context of machine learning occurs for the feature corresponding to orbit $O_{[1,1,0,\dots]}$ (see for example Table 1). Since there are only two options – the two nodes are connected and have therefore exactly one perfect matching, or they are not and have none – the square does not have any effect, and the probability of the orbit is proportional to the number of $1$-matchings of this graph, which is in turn equal to its number of edges. Hence, we have that $p(O_{[1,1,0,\dots]})\propto|E|$, and the hardware natively returns an “edge counting” feature.

III.2 Higher-order moments

The probability of measuring a given photon counting event $\mathbf{n}=[n_{1},\dots,n_{M}]$ can also be interpreted from a slightly different, more physically motivated viewpoint. The $M$ nodes of a graph can be associated with $M$ random variables drawn from a multivariate normal distribution $N(\xi,\Sigma)$, where the covariance matrix $\Sigma$ corresponds to the doubled adjacency matrix $\tilde{A}$, and $\xi$ is the mean vector related to displacement via $\xi=Q^{-1}\mathbf{d}^{\dagger}$. The higher-order moments $E[X_{1}^{(1)}\dots X_{1}^{(n_{1})}\dots X_{M}^{(1)}\dots X_{M}^{(n_{M})}]$ of this distribution are proportional to $\mathrm{haf}(A_{\mathbf{n}})$, which in turn is related to the probability of a photon event via Eq. (3). This result follows from Isserlis’ theorem Isserlis (1918), which decomposes the higher order moments into sums of products of covariances $E[X_{a}X_{b}]$. In short, the GBS device turns a graph into a multivariate normal distribution and samples from its moments.

Using this picture, the first-order moments of the ‘graph-induced distribution’ correspond to photon events of the form $[1,0,\dots]$ and their probability is indeed proportional to the mode means as apparent from Eq. (5). The second-order moments correspond to photon events of the form $[1,1,0,\dots]$ and their probability is proportional to the entries of the adjacency matrix – the edge weights. Consistent with this observation, we stated before that orbits with $2$ non-zero detectors “measure” the edge count of a graph.

While the doubled encoding strategy as well as the presence of multi-photon events somewhat obscure interpretations of features in terms of $r$-matchings and higher-order moments, we found in numerical experiments not reported in this paper that they can be a blessing in disguise, making very similar graphs distinguishable by smaller maximum photon numbers $k$.

III.3 Comparison to Graphlet Sampling kernel

Counting subgraphs in a larger graph is a concept used in various classical graph kernels. Graphlet Sampling kernels Shervashidze et al. (2009) bear the most striking similarity to GBS feature maps, since the features count how often graphlets of size $|V|=3,4,5,\dots$ appear in a graph $G$. In the language developed here we can express the feature $f_{g}$ which counts graphlet $g$ via

[TABLE]

using an indicator function $\mathbbm{1}_{g\cong G_{\mathbf{n}}}$ that is one if graphlet $g$ is isomorphic to the subgraph $G_{\mathbf{n}}$ and zero else, as well as the orbit represented by $[1,..,1,0,\dots]$ counting $|V|$ single photons. In comparison, rewriting Eq. (8) in a similar way, the GBS features are

[TABLE]

where $\mathcal{P}^{|\mathbf{n}|}$ is the set of all perfect matchings of size $|\mathbf{n}|$. As a result, instead of counting graphlets, the GBS feature map sums squares of perfect matching counts in graphlets. Also, GBS feature map does not restrict the size of the graphlet probed.

III.4 Errors due to photon loss

One of the main sources of errors in a realistic GBS device is a photon loss in the linear interferometer, and we demonstrate here that loss is a serious problem for applications of a GBS for graph similarity as proposed in this paper. Methods of dealing with this kind of errors will be discussed in upcoming work. Here we show the effect of the loss on the coarse-grained probabilities with a numerical example.

The effect of loss is described by the action of the lossy bosonic channel on a pure covariance matrix $\sigma$ resulting in

[TABLE]

where $\nu=1-\eta$ and $\eta$ is the overall transmissivity. One way of viewing this is that the matrix $\tilde{A}$ from Eq. ((1)) does not have the block-diagonal structure $c(A\oplus A)$ any more, but is of the form

[TABLE]

where $WX_{2M}\tilde{A}W^{-1}=\bm{\Lambda}$ is the eigendecomposition of $X_{2M}\tilde{A}$. Figure 4 shows the effect of this loss model on the probability distribution $p(O_{\mathbf{n}})$ over orbits for a random unweigthed graph $G$ on ten vertices. It is apparent that loss introduces errors in the distribution, populating orbits which have a zero probability in the zero-displacement case, and distorting the remaining probabilities significantly. In the remainder of the paper we will consider only a lossless GBS device, but remark herewith that loss mitigation strategies are crucial for practical applications of GBS feature maps.

IV Experiments

Finally, we provide some numerical results to investigate the GBS graph kernel in practice. Benchmarks suggest that it is well competitive to standard “classical” graph kernels, at least in the hypothetical case of a perfect device. We furthermore show that displacement may improve classification accuracy by shifting weight into the higher-order orbits, and that orbits with photon numbers smaller or equal to $2$ contribute most to the result.

IV.1 Benchmarking

To benchmark the GBS feature map, we use a setup that has become a standard in testing graph kernels: A C-Support Vector Machine (SVM) with a precomputed kernel. The test accuracies in Table 2 are obtained by running $10$ repeats of a double $10$-fold cross-validation. The inner fold extracts the best model by adjusting the $C$-parameter of the SVM – which controls the penalty on misclassifications – via grid search between values $[10^{-4},10^{3}]$, and the best model is then used to get the accuracy of the test set in the outer cross-validation loop. The GBS feature vectors were used in conjunction with a ‘rbf’ kernel $\kappa_{\mathrm{rbf}}$.

For the GBS graph kernel, we chose a gentle displacement of $d=0.25$ on every mode and $k=6$, leading to $30$-dimensional feature vectors. We used exact simulations based on the hafnian library Björklund et al. (2018). These are computationally very expensive, which is why we only consider small datasets. Three classical graph kernels are benchmarked for comparison: The Graphlet Sampling kernel Shervashidze et al. (2009) (GS) with maximum graphlet size of $k=5$ and $5174$ samples drawn, the Random Walk kernel Gärtner et al. (2003) (RW) with fast computation and a geometric kernel type, and the Subgraph Matching kernel (SM) Kriege and Mutzel (2012). The three classical kernels were simulated using Python’s grakel library Siglidis et al. (2018).555Experiments were run on IBM’s cloud platform using four 2.8GHz Intel Xeon-IvyBridge Ex (E7-4890-V2-PentadecaCore) processors with 15 CPU cores each, as well as on Oak Ridge’s Titan supercomputer.

The datasets are taken from the repository of the Technical University of Dortmund Kersting et al. (2016) (see Figure 5). Preprocessing of the benchmarking datasets includes these three steps:

Graph selection: Graphs which have less than $6$ or more than $25$ nodes are excluded to keep the feature vectors constant and to limit the time of simulations. The share of excluded graphs is displayed in Figure (3) in the main paper, and ranges from $5\%$ to $55\%$. 2. 2.

Labels and attributes: Potential node labels, node attributes and edge attributes are ignored. The edge labels in BZR_MD, COX2_MD, ER_MD, MUTAG and PTC_FM were translated to the following weights: [math] - no chemical bond, $1$ - single bond/double bond/triple bond/aromatic bond. The edge labels in AIDS where translated into the weights: [math] - no edge $1$ - valence of zero, one or two. In FINGERPRINT, only graphs of the three dominant classes $0,4,5$ were considered, since the other classes did not contain a sufficient number of samples after graph selection. 3. 3.

Rescaling: The final (weighed or unweighed) adjacency matrix is divided by a normalization constant $c=1/(\lambda^{\{G\}}_{\mathrm{max}}+10^{-8})$ that is slightly larger than the largest singular value $s^{\{G\}}_{\mathrm{max}}$ of any adjacency matrix in the dataset, as explained in Section 2.1 of the main paper.

All datasets were chosen before the first experiments were run, to avoid a post-selection bias in favour of the GBS kernel.

As Table 2 shows, the GBS kernel performs well and outperforms the other methods visibly for MUTAG and NCI1, while still leading for AIDS, BZR_MD, ER_MD, FINGERPRINT and PROTEINS. Displacement increases the performance of the GBS kernel significantly for COX2_MD, ENZYMES and IMDB-BIN, but not for other data sets. The GBS kernel does well on datasets where the distribution of node and edge numbers differs strongly between classes. However, we confirmed that excluding the ‘edge counting features’ $[1,1,0..],[2,2,0..],\dots$ does not influence classification performance. While the graph size is considered by the GBS kernel, it seems to be only one of many properties that enters the notion of similarity.

IV.2 Displacement and feature importance

The hyperparameters of the GBS and GBS+ graph kernels are the constant displacement $d$ which adiministered to each node, as well as the maximum photon number $k$. Since simulations restrict the value of $k$ at this stage, we focus on the effect of displacement, using the orbit-features (i.e., the GBS kernel). Displacement can change the similarity measure significantly. For example, comparing graphs of size $|V|=3$, one finds that the fully disconnected graph is closer to the fully connected graph than a graph with two edges for $d=1$, but vice versa for $d=0$.

Figure 6 uses the example of IMDB-BIN and MUTAG to investigate the GBSI or “orbit” features for $d=0,d=0.25$ and $d=1$. The feature averages show that the general distribution of the feature vector is similar for both classes, but still visually distinguishable.666Standarization of the feature vectors to emphasize their mutual differences improved classification accuracy in some cases, but deteriorated it in others. Consistent with the theory, increasing displacement shifts the features towards higher-order orbits, and populates features that are zero when $d=0$. Features associated with orbits $[1,1,0,...],[1,1,1,1,0...]$ and $[1,1,1,1,1,1]$, as well as $[2,1,0,...]$ and $[2,1,1,1,0,...]$ seem to be particularly important in the support of principal components, and get high weights when training a perceptron on the GBS features. Where displacement renders them nonzero, uneven orbits such as $[1,1,1,0,...],[1,1,1,1,1,0...]$ follow suit. During our investigations we confirmed that dropping features with high single-detector photon numbers did not have a huge influence on classification. Consistent with the results from Table 2, MUTAG has ‘richer’ features for $d=0$ than IMDB-BIN for classification with a perceptron, an advantage that IMDB-BIN equalizes with growing displacement.

The feature analysis suggests that features related to subgraphs of all sizes (here $1$ to $6$) are important for the classification results, and that duplication of a single node in the subgraphs may be beneficial – a feature that Graphlet Sampling kernels do not explore. The effect of displacement varies with the dataset, and $d$ should therefore be kept as a hyperparameter for model selection.

V Conclusion

We proposed a new type of feature extraction strategy for graph-structured data based on the quantum technique of Gaussian Boson Sampling. We suggested that the success of the method is related to the fact that such a system samples from distributions that are related to useful graph properties. For classical machine learning, this method presents a potentially powerful extension to the gallery of graph kernels, each of which has strengths on certain data sets. For quantum machine learning, this proposes the first application of a “quantum kernel”.

A lot of questions are still open for further investigation, for example regarding the role and interpretation of displacement, how GBS performs with weighted adjacency matrices, how node and edge labels can be considered, as well as whether the feature vectors are useful in combination with other methods such as neural networks. We expect that the rapid current development of numeric GBS samplers as well as quantum hardware will help answering these questions in the near future.

Acknowledgements

We thank Christopher Morris and Nicolas Quesada for valuable advice, as well as the authors of Python’s GraKel library. This research used resources of the Oak Ridge Leadership Computing Facility, which is a DOE Office of Science User Facility supported under Contract DE-AC05-00OR22725.