# Sample efficient graph classification using binary Gaussian boson sampling

Amanuel Anteneh<EMAIL_ADDRESS>Department of Computer Science,
University of Virginia, Charlottesville, Virginia 22903, USA  Olivier Pfister
<EMAIL_ADDRESS>Department of Physics, University of Virginia,
Charlottesville, Virginia 22903, USA

###### Abstract

We present a variation of a quantum algorithm for the machine learning task of
classification with graph-structured data. The algorithm implements a feature
extraction strategy that is based on Gaussian boson sampling (GBS) a near term
model of quantum computing. However, unlike the currently proposed algorithms
for this problem, our GBS setup only requires binary (light/no light)
detectors, as opposed to photon-number-resolving detectors. Binary detectors
are technologically simpler and can operate near room temperature, making our
algorithm much less complex and costly to implement physically. We also
investigate the connection between graph theory and the Torontonian matrix
function which characterizes the probabilities of binary GBS detection events.

††preprint: APS/123-QED

##  I Introduction

Graphs are one of the most versatile data structures used in computing, and
developing machine learning methods for working with graph-structured data has
been a growing sub-field of machine learning research. Graph classification,
in particular, has useful applications in fields such as bioinformatics,
network science and computer vision as many of the objects studied in these
fields can easily be represented as graphs. However, using graph-structured
data with machine learning models is not a straightforward task. This is
because one of the most common ways of representing a graph for computational
applications, i.e., as an adjacency matrix, cannot be easily used as an input
to machine learning classifiers which primarily take vector-valued data as
their inputs. Therefore, a common way of working with graph-structured data is
by defining a feature map $\phi$ that maps a graph $G$ to a vector in a
Hilbert space called a feature space. From there a function $\kappa$, called a
kernel, is defined that measures the similarity of two graphs in the feature
space. An example of a feature map from
$\mathbb{R}^{2}\rightarrow\mathbb{R}^{3}$ is shown in Fig. 1.

Kernel methods refer to machine learning algorithms that learn by comparing
pairs of data points using this similarity measure. In our context we have a
set of graphs $\mathbb{G}$ and we call a kernel $\kappa$ a graph kernel if it
is a function of the form
$\kappa:\mathbb{G}\times\mathbb{G}\rightarrow\mathbb{R}$ [1, 2]. The most
common example of a kernel function is the feature space’s inner product
$\kappa(x,x^{\prime})=\langle\phi(x),\phi(x^{\prime})\rangle$. The goal of
such methods is to construct mappings to feature vectors whose entries (the
features) relate to relevant information about the graphs. Using a Gaussian
boson sampling (GBS) device to construct graph kernels was an idea first
proposed by Schuld et al. in Ref. 3.

Figure 1: In the original input space $\mathbb{R}^{2}$ the data points, which
belong either to the class ‘red’ or ‘blue’, are not separable by a linear
function (the decision boundary) but after mapping the points to feature
vectors in a higher dimensional space $\mathbb{R}^{3}$ a linear function is
able to separate the two classes. This linear decision boundary can be
calculated by supervised machine learning models such as a support vector
machine. In our case the input space is the set of all undirected graphs which
we denote as $\mathbb{G}$.

Boson sampling was first proposed by Aaronson and Arkhipov [4] as a
task—generating the photon-counting outcomes of the “quantum Galton board”
constituted by an $M\times M$ optical interferometer fed with single photons
into some of its input ports—that is strongly believed to be intractable to
classical computers. The reason for this intractability is that calculating
the probability distribution for generating random outcomes using Monte Carlo
simulations requires calculating the permanent of an $M\times M$ matrix.
Calculating the permanent of a general matrix is known to be $\\#\P$-complete
[5] which is a class of problems comparable to the class of $\NP$-complete
problems in their difficulty. Gaussian boson sampling [6] is a variant of
boson sampling in which the single-photon inputs are replaced with single-mode
squeezed states, as produced, for example, by two-photon-emitting optical
parametric amplifiers [7]. The GBS probability distribution is governed by the
Hafnian of an $M\times M$ matrix. Calculating the Hafnian of a general square
matrix can be reduced to the task of calculating permanents therefore
calculating the Hafnian is also $\\#\P$-complete. In both cases, a quantum
machine implementing boson sampling or GBS can easily sample from these hard-
to-calculate probability distributions, just because they are “wired-in,” and
this constitutes the “quantum advantage” that was recently demonstrated in
optical experiments [8, 9]. Note also that the initial “quantum supremacy”
result obtained by Google on a superconducting qubit array [10] was a quantum
(circuit) sampling result as well.

Beyond these necessary initial steps of demonstrating that quantum hardware
can indeed reach regions inaccessible to classical hardware, a subsequent
question is that of the utility of a sampling task. Whereas the usefulness of
sampling in and of itself is far from established, we know that the histograms
produced by statistically significant sampling constitute empirical
probability distributions that tend toward the true, classically intractable
probability distributions for sample numbers linear in the number of possible
outcomes [11]. The problem is that this very number of possible outcomes grows
exponentially with $M$ in a $M$-qubit quantum circuit in general 111Note that
this is not related to the number of possible output quantum states, which
scales with the number of parameters governing the quantum evolution, e.g.
parameters of a simulated Hamiltonian. Obviously, no quantum advantage can be
obtained for $M$-qubit Hamiltonians that have $\mathcal{O}(2^{M})$ parameters
but all classically intractable $M$-qubit Hamiltonians of physical interest
are local and have parameter numbers polynomial in $M$ [40], which validates
Feynman’s proposed advantage for quantum simulation [41]. An $M\times M$
optical interferometer has $M^{2}$ parameters, for example. However, even
though any useful quantum computer will explore but a
$\mathcal{O}(M^{k})$-dimensional region of an $\mathcal{O}(2^{M})$-dimensional
Hilbert space, the number of measurement outcomes will still scale like
$\mathcal{O}(2^{M})$ a priori, simply because we do not know the adequate
measurement basis that best contains the $\mathcal{O}(M^{k})$ output states.
This is the well known exponential overhead of quantum state tomography., and
exponentially or super-exponentially with $M$ in an $M$-optical-mode boson or
Gaussian boson sampler, which dispels any notion of quantum advantage for
calculating the corresponding quantum probability distributions.

One direction that has been explored out of this conundrum is the binning of
GBS measurements results into outcome classes whose cardinality scales
favorably (e.g. polynomially) with the problem size (the GBS mode number). The
immediate downside of such an approach is loss of information it entails,
which impacts usefulness. However, graph classification using feature vectors
and coarse-graining might provide advantageous GBS applications. This was
first pointed out by Schuld et al. [3].

In this paper, we show that a technologically simpler version of GBS, which we
term binary GBS, can achieve comparable or better performance. The paper is
structured as follows. In Sec.II we give broad reminders about GBS and graph
theory (with details in Appendix A) and the current GBS graph kernel from Ref.
3. We then present our graph kernel in Sec.III along with results from
numerical experiments and analyses of its complexity, features and advantages.

Figure 2: Example of a 3-mode Gaussian boson sampler. Mode $i\in\\{1,2,3\\}$
starts in the vacuum state $\ket{0}$, is then squeezed by $\hat{S}(r_{i})$ and
passes through the network of two beamsplitters (the interferometer) before
the number of photons in each mode is measured by the detectors
$D_{i\in\\{1,2,3\\}}$.

##  II  Reminders about Gaussian Boson Sampling (GBS) and graph theory

###  II.1  Gaussian Boson Sampling

As mentioned above, an $M$-mode GBS devise comprises $M$ single-mode-squeezing
(SMS) inputs, an $M\times M$ optical interferometer, and $M$ photon-number-
resolving (PNR) detectors, see Fig.2 for an example. The latter have come of
age in superconducting devices such as transition edge sensors [13] and
superconducting nanowire single-photon detectors [14]. Both the former and the
latter have recently been used to make PNR measurements of as many as 100
photons [15, 16].

An $M$-mode Gaussian boson sampler prepares a Gaussian (Wigner function)
quantum state by the $M$ squeezers and the interferometer. The squeezers
output squeezed light into the interferometer and the photons are then passed
through the interferometer after which the $M$ detectors detect what modes the
photons end up in resulting in a detection event. We denote a detection event
as $\textbf{n}=(n_{1},...,n_{M})$, where $n_{i}$ is the photon count in the
$i$th mode and the total number of photons is $n=\sum_{i=1}^{M}n_{i}$.

We now consider binary detectors, such as single-photon avalanche photodiodes,
which are single-photon sensitive but aren’t PNR and give the same signal
however many photons were absorbed. In this case, we have $n_{i}\in\\{0,1\\}$
where $n_{i}=0$ indicated zero photons were detected in that mode and
$n_{i}=1$ indicates that at least one photon was detected. When using binary
detectors we no longer know the total photon number $n$ so we use $N$ to
denote the number of detectors that detect photons leading to
$\sum_{i=1}^{M}n_{i}=N\leq M$.

An $M$-mode Gaussian state is fully described by a covariance matrix
$\bm{\Sigma}\in\mathbb{R}^{2M\times 2M}$ and a displacement vector
$\textbf{d}\in\mathbb{R}^{2M}$ [17].

###  II.2  Graph theory

In this paper we define a graph $G=(V,E)$ as a set of vertices
$V=\\{v_{1},v_{2},...\\}$ and a set of edges
$E=\\{(v_{1},v_{1}),(v_{1},v_{2}),...(v_{i},v_{j}),...\\}$ that connect
vertices if the edge value is not zero. A graph can be unweighted, with all
nonzero edge weights equal to 1, or weighted, for example with real edge
weights in GBS. For undirected graphs, which is what we will exclusively work
with in this paper, $(v_{i},v_{j})=(v_{j},v_{i})$, $\forall i,j$. The size of
a graph is equal to the cardinality $|V|$ of its vertex set. The degree of a
vertex $v$ is the number of edges that are connected to it. The maximum degree
of a graph is the largest degree of a vertex in its vertex set.

Graphs can be represented in a number of ways such as a diagram, Fig.3(a), or
a more computationally useful way as an adjacency matrix, Fig.3(b). The
adjacency matrix of an undirected graph $G$ with $|V|$ vertices is a
$|V|\times|V|$ symmetric matrix $A$ with entries $a_{ij}$ where $a_{ij}$ is
the weight of the edge connecting vertices $i$ and $j$.

(a) Undirected weighted 4-vertex graph with 4 edges

(b) Adjacency matrix of graph

Figure 3: 4-vertex graph and its corresponding adjacency matrix

A substantial amount of work has been done on the connection between graph
theory and Gaussian boson sampling with PNR detectors [18, 19, 20]. In
Appendix A, we recall some details of this work, namely how a given graph
adjacency matrix can be encoded in a GBS experiment.

###  II.3  Sample complexity of GBS

####  II.3.1 The problem with using GBS beyond sampling

The sample complexity of a machine learning algorithm refers to the number of
samples or amount of data required to learn some target function. In the case
of GBS applications it refers to the number of samples we need to generate
from the GBS device to learn or approximate a probability distribution over
some set of the photon detection events. This complexity type is extremely
important to examine for any applications of GBS as it could potentially
render certain applications of GBS intractable for larger problem sizes.

For example it was shown that the GBS device utilizing PNR detectors can
encode the graph isomorphism problem [19]. This is done by encoding two graphs
into two GBS devices and sampling each $S$ times. The $S$ samples could then
be used, in principle, to reconstruct the probability distribution over all
possible detection events n for a given $M$ and $n$. However, this cannot be
done efficiently enough to provide a quantum advantage. Indeed, we know from
Refs. 11, 21 that reconstructing a probability distribution $D$ over a
discrete finite set $\Omega$ of cardinality $|\Omega|$ from an empirical
distribution $\hat{D}$ constructed from samples from $D$ we require

$S=\left\lceil\frac{2(\ln(2)|\Omega|+\ln(\frac{1}{\delta}))}{\epsilon^{2}}\right\rceil$
(1)

samples to guarantee that

$p(||D-\hat{D}||_{1}\geq\epsilon)\leq\delta$ (2)

where $||D-\hat{D}||_{1}$ denotes the $L_{1}$ distance between $D$ and
$\hat{D}$. In other words we require

$\mathcal{O}\left(\frac{|\Omega|+\ln(\frac{1}{\delta})}{\epsilon^{2}}\right)$
(3)

samples to ensure with probability at most $\delta$ that the sum of the
absolute values of the errors on the empirical probability distribution is
$\epsilon$ or greater.

This means the number of samples we need to approximate a probability
distribution scales linearly with the number of elements in it’s sample space
i.e. the number of outcomes. In the case where $D$ is the probability
distribution over the set of all possible PNR detection events the number of
such events for a given number of modes $M$ and maximum number of photons $n$
is

$|\Omega|={n+M-1\choose M-1}=\frac{(n+M-1)!}{n!(M-1)!}$ (4)

which in number theory is also known as the formula for the number of weak
compositions of an integer $n$ into $M$ parts. As shown in appendices B and C
under the assumption that the number of modes scales quadratically with the
number of photons, $M\in\mathcal{O}(n^{2})$, this quantity grows super
exponentially with $M$ and in general scales as $\mathcal{O}((n+M-1)^{M-1})$
meaning that as the size of the graphs increase, and therefore as the number
of modes of the GBS device increase, we require an exponential number of
samples to ensure the algorithm can give us the correct result within a
certain probability. Therefore while the algorithm may in principle be able to
decide graph isomorphism, it is sample inefficient to an exponential degree
making it intractable to implement even with a fault tolerant quantum
computer.

####  II.3.2 Coarse graining of sample distributions

However a method was suggested in [19] to coarse-grain the probability
distribution by combining outcomes into groups called orbits. Coarse-graining
in this sense means to construct a new probability distribution over the set
of these groups, the cardinality of which is less than the original set of all
possible detection events. An orbit $\textrm{O}_{\textbf{n}}$ consists of a
detection event n and all of its permutations. For example the orbit that
contains the detection event $\textbf{n}=(1,2,2)$ also contains the detection
events $(2,1,2)$ and $(2,2,1)$. The number of orbits for a 4-mode GBS device
is equal to the number of ways one can write $n_{1}+n_{2}+n_{3}+n_{4}=n$,
where the order of the summands does not matter. This is called the number of
partitions of the integer $n$ into $M$ parts and from the number theory
literature [22] it is known to behave asymptotically as

$|\Omega|\approx\frac{e^{\pi\sqrt{\frac{2(n-M)}{3}}}}{4\sqrt{3}(n-M)},M\leq
n\leq 2M.$ (5)

If we assume the number of photons grows linearly with the number of modes,
$n\in\Theta(M)\rightarrow n=2M$, we have the following asymptotic bound on the
number of orbits

$\frac{1}{4\sqrt{3}M}e^{\pi\sqrt{\frac{2M}{3}}}\in\mathcal{O}(\frac{e^{\pi\sqrt{\frac{2M}{3}}}}{M}).$
(6)

This means the number of orbits, which is now the number of outcomes
$|\Omega|$ from Eq. 1, grows like $M^{-1}e^{\pi\sqrt{2M/3}}$ meaning we would
have a sample complexity of

$\mathcal{O}\left(\frac{M^{-1}e^{\pi\sqrt{2M/3}}+\ln(\frac{1}{\delta})}{\epsilon^{2}}\right)$
(7)

which is subexponential but still intractable for large $M$.

####  II.3.3  Sample complexity of previously proposed GBS graph kernels

The first GBS based graph kernel proposed in [3] maps a graph $G$ to feature
vectors in a feature space
$\phi:G\rightarrow\textbf{f}=(f_{1},f_{2},...,f_{D})\in\mathbb{R}^{D}$. Where
$f_{i}=p(\textrm{O}_{\textbf{n}}^{i})$ is the probability of detecting a
detection event from the orbit $\textrm{O}_{\textbf{n}}^{i}$. This kernel was
shown to perform well against three of the four classical kernels we use as
benchmarks in this paper. However a shortcoming of this method is that the
sample complexity is $\mathcal{O}(M^{-1}e^{\pi\sqrt{2M/3}})$.

The second GBS kernel is of the form
$\phi:G\rightarrow\textbf{f}=(f_{1},f_{2},...,f_{D})\in\mathbb{R}^{D}$, with
$f_{i}=p(\mathcal{M}_{n,\Delta_{s}}^{i})$ where
$p(\mathcal{M}_{n,\Delta_{s}}^{i})$ is the probability of detecting a
detection event that belongs to the “meta-orbit”
$\mathcal{M}_{n,\Delta_{s}}^{i}$. A meta-orbit $\mathcal{M}_{n,\Delta_{s}}$ is
uniquely defined by a total photon number $n$ and $\Delta_{s}$ which is
defined as

$\Delta_{s}=\\{\textbf{n}\textrm{ : }\sum_{i}n_{i}=n\;\land\;\forall i\textrm{
: }n_{i}\leq s\;\\}.$ (8)

Therefore a meta-orbit consists of all detection events where total photon
number is equal to $n$, where no detector counts more than $s$ photons. It is
claimed that this strategy partitions the set of all PNR detection events into
a polynomial number of subsets in $n$ [23].

##  III The Algorithm

###  III.1  GBS with binary detectors and its relation to graph theory

While the relationship between GBS with PNR detectors and graph theory has
been thoroughly explored, there has been little exploration of how GBS with
binary detectors fits into the picture. In this section we shed some light on
the relationship between the two. As stated before when using binary detectors
the detection outcomes are of the form
$\textbf{n}_{\textrm{bin}}=(n_{1},...,n_{M})$ where $n_{i}\in\\{0,1\\}$
$\forall i$ and $n_{i}=1$ indicates the $i$th detector detected one or more
photons. The probability of detecting a detection outcome with binary
detectors is characterized by a matrix function called the Torontonian, to
which the same arguments for classical intractability as for the Hafnian can
be extended [24]. The probability of a given binary detection event
$\textbf{n}_{\textrm{bin}}$ is given by

$p(\textbf{n}_{\textrm{bin}})=\frac{\textrm{Tor}(O_{\textbf{n}_{\textrm{bin}}})}{\sqrt{\textrm{det}(Q)}}=\frac{\textrm{Tor}(X\tilde{A}_{\textbf{n}_{\textrm{bin}}})}{\sqrt{\textrm{det}(Q)}}$
(9)

where

$\displaystyle\tilde{A}$ $\displaystyle=(A\oplus A),$ (10) $\displaystyle X$
$\displaystyle=\begin{bmatrix}0&\mathbb{I}\\\ \mathbb{I}&0\end{bmatrix},$ (11)
$\displaystyle Q$ $\displaystyle=(\mathbb{I}_{2M}-X\tilde{A})^{-1},$ (12)
$\displaystyle O$ $\displaystyle=\mathbb{I}-Q^{-1}$ (13)

and Tor() is the Torontonian of a $2N\times 2N$ matrix $A$ defined as

$\textrm{Tor}(A)=\sum_{Z\in
P([N])}(-1)^{|Z|}\frac{1}{\sqrt{\textrm{det}(\mathbb{I}-A_{Z})}}.$ (14)

Where $P([N])$ is the power set, the set of all possible subsets, of the set
$[N]=\\{1,2,...,N\\}$. The probability of a PNR detection event n can be
written in terms of the matrix $O$ as

$p(\textbf{n})=\frac{1}{\sqrt{\textrm{det}(Q)}}\frac{\textrm{Haf}(\tilde{A}_{\textbf{n}})}{\textbf{n}!}=\frac{1}{\sqrt{\textrm{det}(Q)}}\frac{\textrm{Haf}(XO_{\textbf{n}})}{\textbf{n}!}.$
(15)

The probability of a binary GBS detection event is simply the sum of all
probabilities of the corresponding PNR detection events. A useful example to
illustrate this is a 4-mode Gaussian boson sampler programmed according to
some adjacency matrix $A$ of a graph $G$. Suppose we use binary detectors and
measure the detection event $\textbf{n}_{\textrm{bin}}=(1,0,1,0)$. The
corresponding detection events when using PNR detectors would be of the form
$\textbf{n}=(n_{1},0,n_{3},0)$ where $n_{1},n_{3}>0$ and $n_{1}+n_{3}$ is
even. We will define $\mathcal{N}$ to be the set of all possible $4$-mode PNR
detection events with 0’s in the 2nd and 4th index, i.e. only the 2nd and 4th
detectors detect no photons. From this we have

$\begin{split}p((1,0,1,0))=\frac{\textrm{Tor}(X\tilde{A}_{(1,0,1,0)})}{\sqrt{\textrm{det}(Q)}}\\\
=\sum_{\textbf{n}\in\mathcal{N}}p(\textbf{n})=\sum_{\textbf{n}\in\mathcal{N}}\frac{\textrm{Haf}^{2}(A_{\textbf{n}})}{\textbf{n}!\sqrt{\textrm{det}(Q)}}.\end{split}$
(16)

This means the Torontonian of $X\tilde{A}$ is proportional to an infinite sum
of Hafnians as there are an infinite number of integer lists of the form
$(n_{1},0,n_{3},0)$ where $n_{1},n_{3}>0$. In a real GBS experiment, however,
the energy is finite and therefore the measured probabilities of these events
would be equal to a finite version of this sum where all detection events with
total photon number greater than some cutoff photon number vanish from the
series.

In terms of graph theory this means the probability of detecting
$\textbf{n}_{\textrm{bin}}=(1,0,1,0)$ is proportional to the sum of the
squared Hafnians of all possible subgraphs of $G$ of unbounded and even size
with their 2nd and 4th vertices removed. But again in practice the maximum
size of the subgraphs will always be bounded by some maximum photon number for
a real GBS experiment. More generally we have

$p(\textbf{n}_{\textrm{bin}})=\frac{\textrm{Tor}(O_{\textbf{n}_{\textrm{bin}}})}{\sqrt{\textrm{det}(Q)}}=\frac{\textrm{Tor}(X\tilde{A}_{\textbf{n}_{\textrm{bin}}})}{\sqrt{\textrm{det}(Q)}}=\sum_{\textbf{n}\in\mathcal{N}}\frac{\textrm{Haf}^{2}(A_{\textbf{n}})}{\textbf{n}!\sqrt{\textrm{det}(Q)}}$
(17)

where $\mathcal{N}$ is the set of all PNR events that correspond to the binary
detection event $\textbf{n}_{\textrm{bin}}$ 222The Torontonian was also shown
to be a generating function for the Hafnian given by
$\textrm{Haf}(A)=\frac{1}{M!}\frac{d^{M}}{dz^{M}}\textrm{Tor}(zXA)\big{|}_{z=0}$
(18) where $A$ is a $2M\times 2M$ matrix. .

###  III.2 Constructing the feature vectors

Once the GBS device is programmed we generate $S$ samples from the device. For
our algorithm we use binary detectors so each sample is a list of length $M$
with entries either 0 or 1. Once we have these samples we use them to
construct the feature vector of which we have two definitions based on two
coarse-graining strategies.

The first is based on what we call the $\mu$ coarse-graining strategy where we
group together detection events that contain exactly $i$ detector ‘clicks’ or
ones. For example the detection events $(1,0,0)$ and $(0,0,1)$ would be
grouped together since they both contain exactly 1 detector click. These
groups can also be thought of as ‘binary orbits’ since they contain a
detection event and all its permutations. This strategy partitions the set of
all binary detection events into a linear number of disjoint subsets in $N$.
Using this strategy we can define the feature map as
$\phi:G\rightarrow\textbf{f}=(f_{0},f_{1},...,f_{N})\in\mathbb{R}^{N}$. Where
$N$ is the maximum number of detector clicks and $f_{i}=\frac{S_{i}}{S}$ with
$S_{i}$ being the number of samples which contain exactly $i$ ones.
Equivalently this is the probability of detecting an event where exactly $i$
detectors detect a photon.

The second feature map is based on what we call the $\nu$ coarse-graining
strategy. For a 5 mode boson sampler utilizing binary detectors with maximum
click number 5 there are $|\Omega|=32$ possible detection outcomes. This
coarse-graining strategy groups together detection events whose first 5 modes
are one of these 32 outcomes. For example the detection event
$\textbf{n}_{\textrm{bin}}=(0,1,0,0,1,0,1)$ belongs in the group associated
with the detection event $(0,1,0,0,1)$ since they are equal if one is only
concerned with the first 5 modes. This strategy partitions the set of all
detection events of $5$ or more modes into a constant number of subsets, i.e.
32. The feature map based on this strategy is defined as
$\phi:G\rightarrow\textbf{f}=(f_{[0,0,0,0,0]},f_{[1,0,0,0,0]},...,f_{\textbf{n}})\in\mathbb{R}^{32}$.
Where $f_{\textbf{n}}$ is the probability of detecting an event where the
first 5 modes correspond to one of the 32 possible detection outcomes. For
example $f_{[1,0,0,0,0]}$ is the probability that the first detector detects
photons and the following 4 detectors detect vacuum.

Once we construct the feature vector for each graph in the data set we input
them to a machine learning classifier such as a support vector machine.

Table 1:  Graph data set statistics after prepossessing. A more detailed description of these data sets can be found in appendix B of Ref. 3.  Data set | # of graphs | # of classes | avg. # of vertices | avg. # of edges
---|---|---|---|---
AIDS | $1723$ | $2$ | $11.11$ | $11.29$
BZR_MD | $257$ | $2$ | $20.10$ | $197.69$
COX2_MD | $118$ | $2$ | $23.90$ | $274.40$
ENZYMES | $204$ | $6$ | $18.56$ | $36.30$
ER_MD | $357$ | $2$ | $19.27$ | $185.15$
FINGERPRINT | $1080$ | $3$ | $10.58$ | $9.10$
IMDB-BINARY | $806$ | $2$ | $15.98$ | $63.32$
MUTAG | $179$ | $2$ | $17.48$ | $19.23$
NCI1 | $1853$ | $2$ | $19.77$ | $21.27$
PROTEINS | $515$ | $2$ | $15.77$ | $29.37$
PTC_FM | $284$ | $2$ | $13.64$ | $13.99$

###  III.3  Complexity analysis

In this section we discuss, in addition to the time and space complexity, the
sample complexity of our algorithm.

####  III.3.1 Sample Complexity

Since the $n_{i}$’s for binary detection events can be either $0$ or $1$ we
can think of the detection outcomes as binary strings of length $M$ with at
most $M$ ones. The number of binary strings of length $M$ with exactly $i$
ones is ${M\choose i}$. So the number of possible binary detection events, the
number of binary strings of length $M$ with at most $M$ ones, is given by

$|\Omega|=\sum_{i=0}^{M}{M\choose i}.$ (19)

We can show this function grows like $2^{M}$ using the binomial expansion

$2^{M}=(1+1)^{M}=\sum_{i=0}^{M}{M\choose
i}1^{M-i}1^{i}=\sum_{i=0}^{M}{M\choose i}.$ (20)

Therefore we could not simply use the probability of the individual detection
events as features without coarse-graining even when using binary detectors as
we would still need a prohibitively large number of samples to approximate
their probabilities to within a constant error. This was the reason for
introducing the $\nu$ and $\mu$ coarse-graining strategies.

Since the number of outcomes of the $\mu$ distribution scales linearly with
$N$ which is $\leq M$ the sample complexity of approximating the $\mu$ coarse-
grained probability distribution is

$\mathcal{O}\left(\frac{M+\ln(\frac{1}{\delta})}{\epsilon^{2}}\right)$ (21)

which reduces to $\mathcal{O}(M)$ for constant $\epsilon$ and $\delta$. The
sample complexity of approximating the $\nu$ coarse-grained probability
distribution is

$\mathcal{O}\left(\frac{32+\ln(\frac{1}{\delta})}{\epsilon^{2}}\right)$ (22)

which reduces to $\mathcal{O}(1)$ for constant $\epsilon$ and $\delta$.

####  III.3.2 Space Complexity

The size of the $\nu$ feature vectors is constant with respect to the graph
size so the space required is $\mathcal{O}(1)$ and for the $\mu$ feature
vectors the size grows linearly with $N$ which is $\leq M$ so the space
required is $\mathcal{O}(M)$. However storing the adjacency matrix of the
graphs requires $\mathcal{O}(M^{2})$ space complexity.

####  III.3.3 Time Complexity

The time complexity is determined by the most computationally time intensive
step of the algorithm which is encoding the adjacency matrix into the GBS
device. This is the case because the encoding process requires taking the
Takagi decomposition of the matrix $A$ which for a $M\times M$ matrix has time
complexity $\mathcal{O}(M^{3})$ as it is a special case of the singular value
decomposition [26]. However there do exist quantum algorithms for computing
the singular value decomposition of a matrix with complexity that is
polylogarithmic in the size of the matrix [27]. In particular the quantum
singular value estimation algorithm for a $m\times n$ matrix presented in [28]
has complexity $\mathcal{O}(\textrm{polylog}(mn)/\epsilon)$ where $\epsilon$
is an additive error.

##  IV  Numerical experiments

###  IV.1 Implementation details

We used The Walrus python library to classically sample from the GBS output
distribution when running our experiments and the GraKel python library to
fetch the data sets and simulate the classical graph kernels [29, 30].
Classically sampling from a GBS output distribution is very time intensive
even when using binary detectors so we choose to follow the choice made in [3]
and discard graphs with greater than 25 and less than 6 vertices for each data
set. Before sampling from the GBS device we have four parameters we can set:
the maximum number of detector clicks allowed $N$, the average photon number
$\bar{n}$, the displacement on each mode of the GBS device $d$ and lastly the
number of samples generated by the GBS device $S$. We set $N=6,\bar{n}=5$ and
$d=0$ for our results reported here leading to probability distribution of 32
outcomes using the $\nu$ coarse-graining strategy and $7$ outcomes using the
$\mu$ coarse-graining strategy. Using Eq. 1 with $\delta=0.01$ and
$\epsilon=0.06$ we require about $S=15000$ samples for the $\nu$ feature
vectors and about $S=6000$ samples for the $\mu$ feature vectors.

For the machine learning classifier we use a support vector machine with an
RBF kernel $\kappa_{\textrm{rbf}}$. We obtain the accuracies in Table 2 by
running a double 10-fold cross-validation 10 times. The inner fold performs a
grid search through the discrete set of values
$[10^{-4},10^{-3},...,10^{2},10^{3}]$ on the $C$ hyper-parameter of the SVM
which controls the penalty on misclassifications.

###  IV.2 Numerical results from GBS simulation and subsequent classification

We tested our graph kernel on the same data sets used in [3]. We also ignored
vertex labels, vertex attributes and edge attributes and converted all
adjacency matrices to be unweighted.

Four classical graph kernels were used as a benchmark for our algorithms
classification accuracy. The subgraph matching kernel (SM) with time
complexity $\mathcal{O}(kM^{k+1})$ where $M$ is the number of vertices and $k$
the size of the subgraphs being considered [31], the graphlet sampling kernel
(GS) with worst case time complexity $\mathcal{O}(M^{k})$ which can be
optimized to $\mathcal{O}(Md^{k-1})$ for graphs of bounded degree with the
restriction that $k\in\\{3,4,5\\}$, where $k$ is the graphlet size and $d$ is
the maximum degree of the graph [32], the random walk kernel (RW) with time
complexity $\mathcal{O}(M^{3})$ [33] and the shortest path kernel (SP) with
time complexity $\mathcal{O}(M^{4})$ [34]. For the graphlet sampling kernel we
set maximum graphlet size to $k=5$ and draw $5174$ samples, for the random
walk kernel we use fast computation and a geometric kernel type with the decay
factor set to $\lambda=10^{-3}$, for the subgraph matching kernel we set
maximum subgraph size to $k=5$ and for the shortest path kernel we used the
Floyd–Warshall algorithm to calculate shortest paths. The accuracies of all
four classical kernels and our kernel are shown in the top portion of Table 2.
The values in the bottom portion of Table 2 are the accuracies of the original
GBS graph kernels and are taken from [3] where the features vectors were
constructed with $n=6$. Some accuracies for the subgraph matching kernel are
not reported due to it’s $\mathcal{O}(M^{6})$ time complexity for $k=5$ which
required longer than 7 days of computation time for data sets with a high
average number of vertices.

We can see from the table that our kernel is very competitive with both the
classical and PNR based GBS graph kernels and in fact achieves the highest
accuracy on the ENZYMES dataset. The PNR based kernel only obtains a test
accuracy significantly higher than random guessing ($\approx 16\%$) for
ENZYMES only when displacement is applied. Our kernel can be seen as even more
feasible in this regard since we do not require the extra operation of
displacement to reach our level of accuracy. From Fig. 5 we see that for some
datasets such as AIDS and MUTAG there is a strong size imbalance amongst
graphs of different classes. For the AIDS dataset graphs belonging to class 1
are of much smaller size than those in class 0 while for MUTAG the converse is
true. This size imbalance also exists to a lesser extent for the ER_MD, BZR_MD
and PTC_FM datasets and for the FINGERPRINT dataset class 1 has a graph size
distribution significantly different from the other two. Both our GBS kernel
and the PNR based kernel perform well on these size imbalanced datasets which
indicates that graph size is a property that GBS based kernels are sensitive
to and that this sensitivity persists even when binary detectors are used.

Table 2:  Average test accuracies of the support vector machine with different
data sets and graph kernels. The values in bold in the upper/lower section of
the table are the best accuracy obtained for that section. The values in
parenthesis are the standard deviation across the 10 repeats of double cross
validation. GS, RW, SM and SP refer to the graphlet sampling, random walk,
subgraph matching and shortest path kernels respectively.
$\textrm{GBS}^{\textrm{bin}}_{\nu}$ and $\textrm{GBS}^{\textrm{bin}}_{\mu}$
denotes our GBS kernel with binary detectors that use the $\nu$ and $\mu$
coarse-graining strategies to construct the feature vectors respectively.
$\textrm{GBS}^{\textrm{bin+}}_{\nu}$ denotes that the feature associated with
detecting vacuum $[0,0,0,0,0]$ in the first 5 modes was dropped from all
feature vectors. $\textrm{GBS}^{\textrm{PNR}}$ and
$\textrm{GBS}^{\textrm{PNR+}}$ refer to the original GBS kernels with PNR
detectors that use orbit and meta-orbit probabilities as features respectively
with a displacement of $d$ on each mode. *Runtime $>$ 7 days

Data set | $\rm GBS^{\textrm{bin+}}_{\nu}$ | $\rm GBS^{\textrm{bin}}_{\nu}$ | $\rm GBS^{\textrm{bin}}_{\mu}$ | GS | RW | SM | SP
---|---|---|---|---|---|---|---
AIDS | $98.47(0.10)$ | $98.74(0.20)$ | $\textbf{99.53}(0.05)$ | $99.30(0.07)$ | $53.11(11.90)$ | $77.85(2.44)$ | $99.34(0.09)$
BZR_MD | $60.14(1.28)$ | $61.73(0.89)$ | $58.79(1.17)$ | $51.42(3.51)$ | $\textbf{64.54}(0.36)$ | time out* | $50.82(1.76)$
COX2_MD | $\textbf{51.62}(2.76)$ | $50.18(2.96)$ | $51.30(3.86)$ | $49.01(3.18)$ | $48.98(4.78)$ | time out* | $48.11(4.30)$
ENZYMES | $\textbf{48.10}(1.18)$ | $41.75(2.35)$ | $19.83(1.43)$ | $34.59(2.54)$ | $19.50(2.29)$ | $37.38(1.60)$ | $22.15(1.88)$
ER_MD | $67.74(0.94)$ | $69.19(0.33)$ | $68.84(0.50)$ | $48.88(4.53)$ | $\textbf{70.32}(0.02)$ | time out* | $45.23(4.35)$
FINGERPRINT | $64.45(0.78)$ | $\textbf{65.53}(0.86)$ | $63.56(0.67)$ | $65.25(1.30)$ | $33.63(3.57)$ | $46.89(0.56)$ | $46.22(1.02)$
IMDB-BINARY | $60.69(0.84)$ | $61.35(0.98)$ | $67.34(0.38)$ | $\textbf{68.49}(0.63)$ | $67.78(0.38)$ | time out* | $65.50(0.27)$
MUTAG | $84.63(0.91)$ | $\textbf{85.94}(0.98)$ | $81.37(0.90)$ | $80.80(0.91)$ | $83.22(0.04)$ | $83.24(1.27)$ | $82.74(1.65)$
NCI1 | $\textbf{63.45}(0.57)$ | $56.99(1.69)$ | $59.09(1.02)$ | $50.34(3.22)$ | $50.96(3.58)$ | time out* | $53.40(2.25)$
PROTEINS | $\textbf{65.95}(1.03)$ | $63.38(0.73)$ | $63.11(0.55)$ | $65.75(0.94)$ | $56.91(1.39)$ | $62.93(0.83)$ | $63.63(0.41)$
PTC_FM | $52.63(3.95)$ | $57.47(2.72)$ | $59.17(1.58)$ | $\textbf{60.74}(1.48)$ | $50.95(3.68)$ | $56.36(2.66)$ | $55.38(4.04)$

Data set |  $\textrm{GBS}^{\textrm{PNR}}$ ($d=0$) |  $\textrm{GBS}^{\textrm{PNR}}$ ($d=0.25$) |  $\textrm{GBS}^{\textrm{PNR+}}$ ($d=0$) |  $\textrm{GBS}^{\textrm{PNR+}}$ ($d=0.25$)
---|---|---|---|---
AIDS | $99.60(0.05)$ | $\textbf{99.62}(0.03)$ | $99.58(0.06)$ | $99.61(0.05)$
BZR_MD | 62.73(0.71) | 62.13(1.44) | 62.01(1.43) |  63.16(2.11)
COX2_MD | 44.98(1.80) | 50.11(0.97) | 57.84(4.04) |  57.89(2.62)
ENZYMES | 22.29(1.60) | 28.01(1.83) | 25.72(2.60) |  40.42(2.02)
ER_MD | 70.36(0.78) | 70.41(0.47) | 71.01(1.26) |  71.05(0.83)
FINGERPRINT | 65.42(0.49) |  65.85(0.36) | 66.19(0.84) | 66.26(4.29)
IMDB-BINARY | 64.09(0.34) |  68.71(0.59) | 68.14(0.71) | 67.60(0.75)
MUTAG |  86.41(0.33) | 85.58(0.59) | 85.64(0.78) | 84.46(0.44)
NCI1 |  63.61(0.00) | 62.79(0.00) | 63.59(0.17) | 63.11(0.93)
PROTEINS |  66.88(0.22) | 66.14(0.48) | 65.73(0.69) | 66.16(0.76)
PTC_FM | 53.84(0.96) | 52.45(1.78) |  59.14(1.72) | 56.25(2.04)

###  IV.3  Feature analysis

Figure 4: Results of the principal component analysis (PCA) on the $\nu$
feature vector entries for the ENZYMES, MUTAG, IMDB_BINARY and FINGERPRINT
datasets. The heatmaps show the weight/coefficient associated with each
feature with regard to the first four principal components. Figure 5:
Distribution of graph sizes according to class for each dataset.

Fig. 4 shows the results of performing a principal component analysis on the
feature vectors generated using the $\nu$ coarse-graining strategy for various
datasets. The analysis shows that the feature associated with vacuum
$[0,0,0,0,0]$ contributes by far the most in the support of the first
principal component. The analysis also suggests that in some cases the first
10 or so features contribute the most to the support of all of the first four
principal components but in other cases, such as with FINGERPRINT, most
features contribute more or less equally.

###  IV.4  Comparison to classical kernels

Our graph kernel has a time complexity that is equivalent to the random walk
kernel and better than the shortest path kernel by a factor of $M$ while
outperforming both on most data sets. Furthermore the time complexity of our
kernel is not exponential in the size of the subgraphs we are probing like the
subgraph matching kernel. The graphlet sampling kernel does have a more
favorable complexity of $\mathcal{O}(Md^{k-1})$ for graphs with maximum degree
$d$. However it’s important to note that many real world graphs are what are
called ‘scale-free networks’ and from the network science literature [35] the
maximum degree of these graphs grows polynomially with the graph size.
Therefore it is possible that the maximum degree of these graphs grows
linearly with the graph size i.g. $d\in\mathcal{O}(M)$ which would lead to a
complexity of $\mathcal{O}(M^{k})$ for the graphlet sampling kernel. What is
also interesting is that GBS kernels seems to provide more distinguishing
power than some classical kernels for graphs with no vertex and edge labels
like those used in our simulations. Take for example the ENZYMES dataset for
which the binary GBS kernel achieves a classification accuracy of $\approx
48\%$ while the shortest path kernel reaches about 23%. If we instead choose
to not ignore vertex labels we found the shortest path kernel gives a
classification accuracy of about $50\%$. Since the GBS features are related to
Hafnians this suggests that features related to the number of perfect
matchings of a graph could be more useful for distinguishing graphs of
different classes when one has no information about the attributes of the
graph nodes.

##  V Conclusion

We proposed a variation of an algorithm for the machine learning task of
classification with graph-structured data that uses a Gaussian boson sampler
utilizing only binary detectors. We show that our algorithm outperforms four
classical graph kernels for the task of graph classification on many data
sets. This is most evident with regard to the ENZYMES data set where the $\nu$
feature map outperforms all methods. The feature corresponding to detecting
vacuum in the first 5 modes plays a particularly important role as shown by
the principal component analysis as it is related to the Hafnian of all
possible subgraphs of $G$ with their first 5 vertices removed. We also show
that the kernel is sample efficient, a major issue for applications of GBS,
and has a time complexity that is comparable with the classical strategies.

The fact that a GBS kernel using only binary detectors produces such
accuracies suggests that technologically more feasible—binary detectors such
as SPADs do not operate at cryogenic temperatures such as superconducting PNR
ones—GBS devices could have useful applications for machine learning with
graph-structured data. We believe that GBS with PNR detectors should also be
explored more for this application with particular attention given to coarse-
graining strategies that both reduce the sample complexity as well as provide
features that capture useful information about the graphs.

A number of questions remain open for investigation such as how vertex and
edge labels can be encoded into the GBS device. Also as stated earlier it is
known that the existence of a polynomial-time classical algorithm for exact
sampling from the output probability distribution of a boson sampling or
Gaussian boson sampling device would imply the collapse of the polynomial
hierarchy to the third level and thus the existence of such an algorithm is
believed to be very unlikely 333Although this has been proven rigorously for
the exact sampling case [4, 42] the proof pertaining to the approximate
sampling case rests on the assumption of two conjectures known as the
Permanent-of-Gaussians Conjecture and the Permanent Anti-Concentration
Conjecture which are as of now still unproven.. This result can also be
extended to GBS with binary detectors [24]. However it is not known, although
some work has been done in this area [23], if such arguments exist for
algorithms that sample from coarse-grained versions of these probability
distributions such as those defined in [3] or our work. It is important to
know if such arguments exist as they would imply these quantum kernels are
also likely hard to simulate classically.

###### Acknowledgements.

We thank Maria Schuld, Kamil Brádler, Scott Aaronson, Ignacio Cirac, Miller
Eaton, Nicolás Quesada, Andrew Blance, and Sefonias Maereg for useful advice
and discussions. We thank Research Computing at the University of Virginia for
providing access to, and support with, the Rivanna computing cluster. This
work was supported by NSF grant PHY-2112867.

##  Appendix A Reminders about standard GBS

###  A.1 GBS with PNR detectors

There has been substantial work done already on the connection between graph
theory and Gaussian boson sampling with PNR detectors [18, 19, 20]. Here we
present the important concepts. Any undirected graph $G$ with no self-loops
and $|V|=M$ vertices can be encoded into a $M$-mode GBS setup consisting of a
set of $M$ squeezers followed by an interferometer of beamsplitters according
to its adjacency matrix $A$. Once the graph is encoded into the GBS device the
probability of detecting a specific detection event
$\textbf{n}=(n_{1},...,n_{M})$ is equal to

$p(\textbf{n})=\frac{1}{\sqrt{\textrm{det}(Q)}}\frac{\textrm{Haf}(\tilde{A}_{\textbf{n}})}{\textbf{n}!}=\frac{1}{\sqrt{\textrm{det}(Q)}}\frac{\textrm{Haf}^{2}(A_{\textbf{n}})}{\textbf{n}!}$
(23)

with

$Q=(\mathbb{I}_{2M}-X\tilde{A})^{-1},\quad X=\begin{bmatrix}0&\mathbb{I}\\\
\mathbb{I}&0\end{bmatrix},$ (24)

$\textbf{n}!=n_{1}!\times...\times n_{M}!$, $\tilde{A}=(A\oplus A)$ and Haf()
denoting the Hafnian of a $2M\times 2M$ matrix. The Hafnian is a matrix
function defined mathematically as

$\textrm{Haf}(A)=\sum_{\pi\in S_{M}}\prod_{(u,v)\in\pi}A_{u,v},$ (25)

where $S_{M}$ is the partition of the set $\\{1,2,...,2M\\}$ into unordered
disjoint pairs. For example if $M=2$ then
$S_{M}=(\\{(1,2),(3,4)\\},\\{(1,4),(2,3)\\},\\{(1,3),(2,4)\\})$. If $A$ is the
adjacency matrix of an unweighted graph then the Hafnian is equal to the
number of perfect matchings of the vertices of the graph. A perfect matching
is a partition of the vertex set of a graph into pairs such that each vertex
is connected to exactly one edge from the edge set. All perfect matchings of a
complete 4-vertex graph are shown in Fig.6.

(a) Complete graph of 4 vertices

(b) The three perfect matchings of the complete 4-vertex graph

Figure 6: The complete graph of 4 vertices and its corresponding perfect
matching

$A_{\textbf{n}}$ is the $n\times n$ submatrix of $A$ induced according to the
photon detection event n. $A_{\textbf{n}}$ is obtained by repeating the $i$th
row and column according to the measurement pattern n. If $n_{i}=0$ then the
$i$th row and column are deleted from $A$ but if $n_{i}>0$ then the $i$th row
and column are repeated $n_{i}$ times. For example the probability of
detecting the event where each mode has exactly one photon
$\textbf{n}=(1,1,...,1)$ would be proportional to the Hafnian of the original
matrix $A$ since $A_{\textbf{n}}=A$. What this means in terms of the graph is
that vertex $i$ and all its edges are either deleted if $n_{i}=0$ or
duplicated $n_{i}$ times if $n_{i}>0$. Therefore the probability of a
detection event n is proportional to the squared Hafnian of the subgraph
$G_{\textbf{n}}$ corresponding to the induced adjacency matrix
$A_{\textbf{n}}$. Examples of different detection events and their
corresponding induced subgraphs are shown in Fig.7.

Figure 7:  Different photon detection events n and the corresponding subgraphs
$G_{\textbf{n}}$ they induce and the value of the squared Hafnians of those
subgraphs. The probability of the detection event where each detector detects
one photon corresponds to the Hafnian of the graph encoded into the GBS. We
can see in the third graph from the top when a detector detects 2 photons the
corresponding vertices and their respective edges are duplicated.

These induced subgraphs are of even size since the number of photons detected
is always even due to the fact that the inputs are squeezed states. However
when displacement is applied to the modes of the GBS the probability of
detecting an odd number of photons is in general not zero anymore and the
probability of individual detection events is characterized by the loop
Hafnian lHaf() as opposed to the Hafnian [37, 38]. We do not apply
displacement for the numerical experiments done in this paper.

###  A.2 Encoding a graph into a GBS device

To map a graph to a feature vector we must first program the GBS device, by
setting the squeezing parameters and beamsplitter angles of the device,
according to the adjacency matrix $A$ of the graph. Any adjacency matrix
$A\in\mathbb{R}^{M\times M}$ of an undirected graph of $M$ vertices can be
mapped to a symmetric, positive definite $2M\times 2M$ covariance matrix
$\bm{\Sigma}$ of a pure Gaussian state of $M$ modes via the following
procedure. First a doubled adjacency matrix $\tilde{A}$ is constructed,

$\tilde{A}=c\begin{bmatrix}A&0\\\ 0&A\end{bmatrix}=c(A\oplus A),$ (26)

where $c$ is a rescaling constant chosen such that
$0<c<1/\lambda_{\textrm{max}}$ where $\lambda_{\textrm{max}}$ is the maximum
singular value of $A$ [3]. We use $\tilde{A}$ as, unlike $A$, it is guaranteed
to map to a covariance matrix of a pure Gaussian state which is easier to
prepare than a mixed one [18]. This also has the advantage of allowing us to
utilize the identity $\textrm{Haf}(A\oplus A)=\textrm{Haf}^{2}(A)$ to relate
$\tilde{A}$ to $A$. To map $\tilde{A}$ to a covariance matrix $\bm{\Sigma}$ we
use the following matrix equations

$\bm{\Sigma}=Q-\mathbb{I}_{2M}/2,\textrm{ with
}Q=(\mathbb{I}_{2M}-X\tilde{A})^{-1},\quad X=\begin{bmatrix}0&\mathbb{I}\\\
\mathbb{I}&0\end{bmatrix}.$ (27)

To program the GBS device to sample from the probability distribution
corresponding to the covariance matrix $\bm{\Sigma}$ of the pure Gaussian
state we need the unitary matrix $U$ that characterizes the interferometer of
the device as well as the squeezing parameters $r_{1},...,r_{M}$ of each of
the $M$ squeezers. We can obtain these values by taking the Takagi
decomposition of $A$ which is of the form

$A=U\textrm{diag}(\lambda_{1},...,\lambda_{M})U^{T}.$ (28)

The squeezing parameters are determined by the singular values
$\lambda_{1},...,\lambda_{M}$ and $c$ via the relationship
$r_{i}=\tanh^{-1}(c\lambda_{i})$. The singular values and $c$ also uniquely
determine the mean photon number $\bar{n}$ of the device according to

$\bar{n}=\sum_{i=1}^{M}\frac{(c\lambda_{i})^{2}}{1-(c\lambda_{i})^{2}}=\sum_{i=1}^{M}\sinh^{2}(r_{i}).$
(29)

The rescaling constant $c$ can be used to adjust $\bar{n}$ as multiplying $A$
by $c$ scales it’s singular values without changing the structure of the graph
other than scaling all edge weights by $c$. The matrix $U$ can be decomposed
to give the parameters of the beamsplitter gates of the interferometer [39].

The GBS device, if using PNR detectors, now samples from the probability
distribution

$p(\textbf{n})=\frac{1}{\sqrt{\textrm{det}(Q)}}\frac{\textrm{Haf}(\tilde{A}_{\textbf{n}})}{\textbf{n}!}=\frac{1}{\sqrt{\textrm{det}(Q)}}\frac{\textrm{Haf}^{2}(A_{\textbf{n}})}{\textbf{n}!}.$
(30)

##  Appendix B Super Exponential Growth of GBS Detection Events for
$M\in\mathcal{O}(n^{2})$

######  Lemma 1.

$\frac{(n+M-1)!}{n!(M-1)!}\in\omega(\lfloor\sqrt{M}\rfloor^{\lfloor{\sqrt{M}}\rfloor})$
for $n=\lfloor\sqrt{M}\rfloor$

###### Proof.

$\displaystyle\frac{(n+M-1)!}{n!(M-1)!}$
$\displaystyle\xrightarrow{n=\lfloor\sqrt{M}\rfloor}\frac{(\lfloor\sqrt{M}\rfloor+M-1)!}{(\lfloor\sqrt{M}\rfloor)!(M-1)!}$
$\displaystyle\frac{(\lfloor\sqrt{M}\rfloor+M-1)!}{(\lfloor\sqrt{M}\rfloor)!(M-1)!}$
$\displaystyle=\frac{[\prod_{i=1}^{\lfloor\sqrt{M}\rfloor}(M-1+i)](M-1)!}{[\prod_{i=1}^{\lfloor\sqrt{M}\rfloor}i](M-1)!}$
$\displaystyle=\frac{[\prod_{i=1}^{\lfloor\sqrt{M}\rfloor}(M-1+i)]}{[\prod_{i=1}^{\lfloor\sqrt{M}\rfloor}i]}$
$\displaystyle=\prod_{i=1}^{\lfloor\sqrt{M}\rfloor}[\frac{M-1}{i}+1]$
$\displaystyle>\prod_{i=1}^{\lfloor\sqrt{M}\rfloor}[\frac{M-1}{\lfloor\sqrt{M}\rfloor}+1]$
$\displaystyle=(\frac{M-1}{\lfloor\sqrt{M}\rfloor}+1)^{\lfloor\sqrt{M}\rfloor}$
$\displaystyle=(\frac{M}{\lfloor\sqrt{M}\rfloor}+1-\frac{1}{\lfloor\sqrt{M}\rfloor})^{\lfloor\sqrt{M}\rfloor}$
$\displaystyle\geq{\lfloor\sqrt{M}\rfloor}^{\lfloor\sqrt{M}\rfloor}$

Therefore
$\frac{(n+M-1)!}{n!(M-1)!}\in\omega(\lfloor\sqrt{M}\rfloor^{\lfloor{\sqrt{M}}\rfloor})$
for $n=\lfloor\sqrt{M}\rfloor$. ∎

##  Appendix C Induction Proof for ${n\choose k}\in\Theta(n^{k})$

######  Lemma 2.

${n\choose k}\in\Theta(n^{k})$

###### Proof.

Base Case: $k=2$

$\displaystyle{n\choose 2}$ $\displaystyle=\frac{n(n-1)}{2!}$
$\displaystyle\lim_{n\to\infty}\frac{\frac{n(n-1)}{2!}}{n^{2}}=\frac{1}{2!}$
$\displaystyle 0<\frac{1}{2!}<\infty$ $\displaystyle\therefore{n\choose
2}\in\Theta(n^{2})$

Assume result holds up to $k=\ell$

$\displaystyle{n\choose\ell}$
$\displaystyle=\frac{n(n-1)(n-2)\cdots(n-\ell+1)}{\ell!}\in\Theta(n^{\ell})$

Inductive Step: $k=\ell+1$

$\displaystyle{n\choose\ell+1}$
$\displaystyle=\frac{n(n-1)(n-2)\cdots(n-\ell)}{(\ell+1)!}$
$\displaystyle\lim_{n\to\infty}\frac{\frac{n(n-1)(n-2)\cdots(n-\ell)}{(\ell+1)!}}{n^{\ell+1}}$
$\displaystyle=$
$\displaystyle\lim_{n\to\infty}\frac{\frac{n(n-1)(n-2)\cdots(n-\ell+1)}{\ell!}}{n^{\ell}}\frac{\frac{(n-\ell)}{\ell+1}}{n}$
$\displaystyle=$
$\displaystyle\lim_{n\to\infty}\frac{\frac{n(n-1)(n-2)\cdots(n-\ell+1)}{\ell!}}{n^{\ell}}\lim_{n\to\infty}\frac{\frac{(n-\ell)}{\ell+1}}{n}$
$\displaystyle\equiv$ $\displaystyle\frac{1}{\ell!}\frac{1}{\ell+1}$
$\displaystyle=$ $\displaystyle\frac{1}{(\ell+1)!}$ $\displaystyle
0<\frac{1}{(\ell+1)!}<\infty$
$\displaystyle\therefore{n\choose\ell+1}\in\Theta(n^{\ell+1})$

∎

## References

  * Nikolentzos  _et al._ [2021] G. Nikolentzos, G. Siglidis, and M. Vazirgiannis, Graph kernels: A survey, Journal of Artificial Intelligence Research 72, 943 (2021).
  * Kriege  _et al._ [2020] N. M. Kriege, F. D. Johansson, and C. Morris, A survey on graph kernels, Applied Network Science 5, 1 (2020).
  * Schuld  _et al._ [2020] M. Schuld, K. Brádler, R. Israel, D. Su, and B. Gupt, Measuring the similarity of graphs with a gaussian boson sampler, Physical Review A 101, 032314 (2020).
  * Aaronson and Arkhipov [2010] S. Aaronson and A. Arkhipov, The computational complexity of linear optics, Electronic Colloquium on Computational Complexity Report No. 170, 1 (2010).
  * Valiant [1979] L. G. Valiant, The complexity of computing the permanent, Theoretical computer science 8, 189 (1979).
  * Hamilton  _et al._ [2017] C. S. Hamilton, R. Kruse, L. Sansoni, S. Barkhofen, C. Silberhorn, and I. Jex, Gaussian boson sampling, Phys. Rev. Lett. 119, 170501 (2017).
  * Bachor and Ralph [2019] H.-A. Bachor and T. C. Ralph,  _A Guide to Experiments in Quantum Optics_ , 3rd ed. (Wiley-VCH, 2019).
  * Zhong  _et al._ [2020] H.-S. Zhong, H. Wang, Y.-H. Deng, M.-C. Chen, L.-C. Peng, Y.-H. Luo, J. Qin, D. Wu, X. Ding, Y. Hu, P. Hu, X.-Y. Yang, W.-J. Zhang, H. Li, Y. Li, X. Jiang, L. Gan, G. Yang, L. You, Z. Wang, L. Li, N.-L. Liu, C.-Y. Lu, and J.-W. Pan, Quantum computational advantage using photons, Science 370, 1460 (2020), https://science.sciencemag.org/content/370/6523/1460.full.pdf .
  * Madsen  _et al._ [2022] L. S. Madsen, F. Laudenbach, M. F. Askarani, F. Rortais, T. Vincent, J. F. F. Bulmer, F. M. Miatto, L. Neuhaus, L. G. Helt, M. J. Collins, A. E. Lita, T. Gerrits, S. W. Nam, V. D. Vaidya, M. Menotti, I. Dhand, Z. Vernon, N. Quesada, and J. Lavoie, Quantum computational advantage with a programmable photonic processor, Nature 606, 75 (2022).
  * Arute  _et al._ [2019] F. Arute, K. Arya, R. Babbush, D. Bacon, J. C. Bardin, R. Barends, R. Biswas, S. Boixo, F. G. S. L. Brandao, D. A. Buell, B. Burkett, Y. Chen, Z. Chen, B. Chiaro, R. Collins, W. Courtney, A. Dunsworth, E. Farhi, B. Foxen, A. Fowler, C. Gidney, M. Giustina, R. Graff, K. Guerin, S. Habegger, M. P. Harrigan, M. J. Hartmann, A. Ho, M. Hoffmann, T. Huang, T. S. Humble, S. V. Isakov, E. Jeffrey, Z. Jiang, D. Kafri, K. Kechedzhi, J. Kelly, P. V. Klimov, S. Knysh, A. Korotkov, F. Kostritsa, D. Landhuis, M. Lindmark, E. Lucero, D. Lyakh, S. Mandrà, J. R. McClean, M. McEwen, A. Megrant, X. Mi, K. Michielsen, M. Mohseni, J. Mutus, O. Naaman, M. Neeley, C. Neill, M. Y. Niu, E. Ostby, A. Petukhov, J. C. Platt, C. Quintana, E. G. Rieffel, P. Roushan, N. C. Rubin, D. Sank, K. J. Satzinger, V. Smelyanskiy, K. J. Sung, M. D. Trevithick, A. Vainsencher, B. Villalonga, T. White, Z. J. Yao, P. Yeh, A. Zalcman, H. Neven, and J. M. Martinis, Quantum supremacy using a programmable superconducting processor, Nature 574, 505 (2019).
  * Weissman  _et al._ [2003] T. Weissman, E. Ordentlich, G. Seroussi, S. Verdu, and M. J. Weinberger, Inequalities for the $l_{1}$ deviation of the empirical distribution, HP Technical Reports HPL-2003-97, R1 (2003).
  * Note [1] Note that this is not related to the number of possible output quantum states, which scales with the number of parameters governing the quantum evolution, e.g. parameters of a simulated Hamiltonian. Obviously, no quantum advantage can be obtained for $M$-qubit Hamiltonians that have $\mathcal{O}(2^{M})$ parameters but all classically intractable $M$-qubit Hamiltonians of physical interest are local and have parameter numbers polynomial in $M$ [40], which validates Feynman’s proposed advantage for quantum simulation [41]. An $M\times M$ optical interferometer has $M^{2}$ parameters, for example. However, even though any useful quantum computer will explore but a $\mathcal{O}(M^{k})$-dimensional region of an $\mathcal{O}(2^{M})$-dimensional Hilbert space, the number of measurement outcomes will still scale like $\mathcal{O}(2^{M})$ a priori, simply because we do not know the adequate measurement basis that best contains the $\mathcal{O}(M^{k})$ output states. This is the well known exponential overhead of quantum state tomography.
  * Lita  _et al._ [2008] A. E. Lita, A. J. Miller, and S. W. Nam, Counting near-infrared single-photons with 95% efficiency, Opt. Expr. 16, 3032 (2008).
  * Cahall  _et al._ [2017] C. Cahall, K. L. Nicolich, N. T. Islam, G. P. Lafyatis, A. J. Miller, D. J. Gauthier, and J. Kim, Multi-photon detection using a conventional superconducting nanowire single-photon detector, Optica 4, 1534 (2017).
  * Eaton  _et al._ [2022] M. Eaton, A. Hossameldin, R. J. Birrittella, P. M. Alsing, C. C. Gerry, H. Dong, C. Cuevas, and O. Pfister, Resolution of 100 photons and quantum generation of unbiased random numbers, Nature Photonics 10.1038/s41566-022-01105-9 (2022).
  * Cheng  _et al._ [2022] R. Cheng, Y. Zhou, S. Wang, M. Shen, T. Taher, and H. X. Tang, A 100-pixel photon-number-resolving detector unveiling photon statistics, Nature Photonics 10.1038/s41566-022-01119-3 (2022).
  * Weedbrook  _et al._ [2012] C. Weedbrook, S. Pirandola, R. García-Patrón, N. J. Cerf, T. C. Ralph, J. H. Shapiro, and S. Lloyd, Gaussian quantum information, Reviews of Modern Physics 84, 621 (2012).
  * Brádler  _et al._ [2018] K. Brádler, P.-L. Dallaire-Demers, P. Rebentrost, D. Su, and C. Weedbrook, Gaussian boson sampling for perfect matchings of arbitrary graphs, Physical Review A 98, 032310 (2018).
  * Brádler  _et al._ [2021] K. Brádler, S. Friedland, J. Izaac, N. Killoran, and D. Su, Graph isomorphism and gaussian boson sampling, Special Matrices 9, 166 (2021).
  * Arrazola and Bromley [2018] J. M. Arrazola and T. R. Bromley, Using gaussian boson sampling to find dense subgraphs, Physical review letters 121, 030503 (2018).
  * Canonne [2020] C. L. Canonne, A short note on learning discrete distributions, arXiv preprint arXiv:2002.11457 (2020).
  * Oruç [2016] A. Y. Oruç, On number of partitions of an integer into a fixed number of positive integers, Journal of Number Theory 159, 355 (2016).
  * Bradler  _et al._ [2019] K. Bradler, R. Israel, M. Schuld, and D. Su, A duality at the heart of gaussian boson sampling, arXiv preprint arXiv:1910.04022 (2019).
  * Quesada  _et al._ [2018] N. Quesada, J. M. Arrazola, and N. Killoran, Gaussian boson sampling using threshold detectors, Physical Review A 98, 062322 (2018).
  * Note [2] The Torontonian was also shown to be a generating function for the Hafnian given by  $\textrm{Haf}(A)=\frac{1}{M!}\frac{d^{M}}{dz^{M}}\textrm{Tor}(zXA)\big{|}_{z=0}$ (31) where $A$ is a $2M\times 2M$ matrix.
  * Hahn [2006] T. Hahn, Routines for the diagonalization of complex matrices, arXiv preprint physics/0607103 (2006).
  * Gu  _et al._ [2019] L. Gu, X. Wang, and G. Zhang, Quantum higher order singular value decomposition (IEEE, 2019) pp. 1166–1171.
  * Kerenidis and Prakash [2016] I. Kerenidis and A. Prakash, Quantum recommendation systems, arXiv preprint arXiv:1603.08675 (2016).
  * Gupt  _et al._ [2019] B. Gupt, J. Izaac, and N. Quesada, The walrus: a library for the calculation of hafnians, hermite polynomials and gaussian boson sampling, Journal of Open Source Software 4, 1705 (2019).
  * Siglidis  _et al._ [2020] G. Siglidis, G. Nikolentzos, S. Limnios, C. Giatsidis, K. Skianis, and M. Vazirgiannis, Grakel: A graph kernel library in python., J. Mach. Learn. Res. 21, 1 (2020).
  * Kriege and Mutzel [2012] N. Kriege and P. Mutzel, Subgraph matching kernels for attributed graphs, arXiv preprint arXiv:1206.6483 (2012).
  * Shervashidze  _et al._ [2009] N. Shervashidze, S. Vishwanathan, T. Petri, K. Mehlhorn, and K. Borgwardt, Efficient graphlet kernels for large graph comparison, in  _Artificial intelligence and statistics_ (PMLR, 2009) pp. 488–495.
  * Vishwanathan  _et al._ [2010] S. V. N. Vishwanathan, N. N. Schraudolph, R. Kondor, and K. M. Borgwardt, Graph kernels, Journal of Machine Learning Research 11, 1201 (2010).
  * Borgwardt and Kriegel [2005] K. M. Borgwardt and H.-P. Kriegel, Shortest-path kernels on graphs, in  _Fifth IEEE international conference on data mining (ICDM’05)_ (IEEE, 2005) pp. 8–pp.
  * Barabási and Pósfai [2016] A.-L. Barabási and M. Pósfai,  _Network Science_ (Cambridge University Press, Cambridge, 2016).
  * Note [3] Although this has been proven rigorously for the exact sampling case [4, 42] the proof pertaining to the approximate sampling case rests on the assumption of two conjectures known as the Permanent-of-Gaussians Conjecture and the Permanent Anti-Concentration Conjecture which are as of now still unproven.
  * Quesada [2019] N. Quesada, Franck-condon factors by counting perfect matchings of graphs with loops, The Journal of chemical physics 150, 164113 (2019).
  * Bulmer  _et al._ [2022] J. F. Bulmer, S. Paesani, R. S. Chadwick, and N. Quesada, Threshold detection statistics of bosonic states, arXiv preprint arXiv:2202.04600 (2022).
  * Clements  _et al._ [2016] W. R. Clements, P. C. Humphreys, B. J. Metcalf, W. S. Kolthammer, and I. A. Walmsley, Optimal design for universal multiport interferometers, Optica 3, 1460 (2016).
  * Lloyd [1996] S. Lloyd, Universal quantum simulators, Science 273, 1073 (1996).
  * Feynman [1982] R. P. Feynman, Simulating physics with computers, Int. J. Theor. Phys. 21, 467 (1982).
  * Kruse  _et al._ [2019] R. Kruse, C. S. Hamilton, L. Sansoni, S. Barkhofen, C. Silberhorn, and I. Jex, Detailed study of gaussian boson sampling, Physical Review A 100, 032326 (2019).


# On a conjecture of Iizuka

Azizul Hoque<EMAIL_ADDRESS>
###### Abstract.

For a given odd positive integer $n$ and an odd prime $p$, we construct an
infinite family of quadruples of imaginary quadratic fields
$\mathbb{Q}(\sqrt{d})$, $\mathbb{Q}(\sqrt{d+1})$, $\mathbb{Q}(\sqrt{d+4})$ and
$\mathbb{Q}(\sqrt{d+4p^{2}})$ with $d\in\mathbb{Z}$ such that the class number
of each of them is divisible by $n$. Subsequently, we show that there is an
infinite family of quintuples of imaginary quadratic fields
$\mathbb{Q}(\sqrt{d})$, $\mathbb{Q}(\sqrt{d+1})$, $\mathbb{Q}(\sqrt{d+4})$,
$\mathbb{Q}(\sqrt{d+36})$ and $\mathbb{Q}(\sqrt{d+100})$ with $d\in\mathbb{Z}$
whose class numbers are all divisible by $n$. Our results provide a complete
proof of Iizuka’s conjecture (in fact a generalization of it) for the case
$m=1$. Our results also affirmatively answer a weaker version of (a
generalization of) Iizuka’s conjecture for $m\geq 4$.

###### Key words and phrases:

Imaginary quadratic field, Class number, Iizuka’s conjecture, Exponent

###### 2020 Mathematics Subject Classification:

11R11; 11R29

##  1\. Introduction

It has been proved that there are infinitely many real (resp. imaginary)
quadratic fields with class numbers divisible by a given positive integer (see
[1, 3, 21, 24]). An analogous problem for tuples of quadratic fields arises
from Scholz’s Spiegelungssatz [19]. In [13], Komatsu studied this problem for
a pair of quadratic fields and proved that there are infinitely many pairs of
quadratic fields $\mathbb{Q}(\sqrt{d})$ and $\mathbb{Q}(\sqrt{md})$ with
$m,d\in\mathbb{Z}$ whose class numbers are divisible by $3$. Later, he
generalized this result in [14] to $n$-divisibility of the class numbers of
pairs of imaginary quadratic fields. On the other hand, Iizuka [11] studied a
slight variant of this problem and construct an infinite family of pairs of
imaginary quadratic fields $\mathbb{Q}(\sqrt{d})$ and $\mathbb{Q}(\sqrt{d+1})$
with $d\in\mathbb{Z}$ whose class numbers are divisible by $3$. Further, he
posed the following conjecture in the same paper.

######  Conjecture 1.1.

For any prime number $p$ and any positive integer $m$, there is an infinite
family of $m+1$ successive real (or imaginary) quadratic fields

$\mathbb{Q}(\sqrt{d}),\mathbb{Q}(\sqrt{d+1}),\cdots,\mathbb{Q}(\sqrt{d+m})$

with $d\in\mathbb{Z}$ whose class numbers are divisible by $p$.

In [4], Chattopadhyay and Muthukrishnan extended Iizuka’s result from pairs to
certain triples of imaginary quadratic fields following the methods used in
[11]. In other words, they gave an affirmative answer of a weaker version of
Conjecture 1.1 for $p=3$. It follows from a recent result of Iizuka, Konomi
and Nakano [12] that Conjecture 1.1 is true for $m=1$ when $p=3,5,7$.

Very recently, Krishnamoorthy and Pasupulati [15] cleverly used [16, Theorem
1] and an extended version of [9, Theorem 3.2] to settled Conjecture 1.1 for
$m=1$.

In this paper, for a given odd prime $p$ we construct an infinite family of
quadruples of imaginary quadratic fields $\mathbb{Q}(\sqrt{d})$,
$\mathbb{Q}(\sqrt{d+1})$, $\mathbb{Q}(\sqrt{d+4})$ and
$\mathbb{Q}(\sqrt{d+4p^{2}})$ with $d\in\mathbb{Z}$ whose class numbers are
all divisible by a given odd integer $n\geq 3$. This extends the results of
[4, 12, 15] in both directions; from pairs to quadruples/quituples of fields
and from primes to odd integers. It also gives a proof of a weaker version of
Conjecture 1.1 for any prime $p\geq 3$ (in fact for any odd integer $n\geq
3$). The precise statement of our first result is the following:

######  Theorem 1.1.

For any odd positive integer $n$ and any odd prime $p$, there are infinitely
many quadruples of imaginary quadratic fields $\mathbb{Q}(\sqrt{d})$,
$\mathbb{Q}(\sqrt{d+1})$, $\mathbb{Q}(\sqrt{d+4})$ and
$\mathbb{Q}(\sqrt{d+4p^{2}})$ whose class numbers are all divisible by $n$.

Note that we can construct an infinite family of quintuples or higher tuples
of imaginary quadratic fields whose class numbers are all divisible by $n$ by
choosing different values for the prime $p$. For instance, utilizing Corollary
3.1 and following the proof of Theorem 1.1, we get the following:

######  Theorem 1.2.

For a given odd integer $n\geq 1$, there are infinitely many quintuples of
imaginary quadratic fields

$\left(\mathbb{Q}(\sqrt{d}),\mathbb{Q}(\sqrt{d+1}),\mathbb{Q}(\sqrt{d+4}),\mathbb{Q}(\sqrt{d+36}),\mathbb{Q}(\sqrt{d+100})\right)$

with $d\in\mathbb{Z}$ whose class numbers are all divisible by $n$.

Our method relies on the prominent result of Bilu, Hanrot and Voutier [2]
concerning the primitive divisors of Lehmer numbers as well as on the
solvability of certain Lebesgue-Ramanujan-Nagell type equations considered in
[8]. This method does not allow us to include an imaginary quadratic field of
the form $\mathbb{Q}(\sqrt{d+m})$ in the tuple, when $m$ is non-square
integer.

##  2\. $n$-divisibility of the class-numbers of
$\mathbb{Q}\left(\sqrt{1-4U^{n}}\right)$ and
$\mathbb{Q}\left(\sqrt{1-V^{n}}\right)$

Here, we recall some results concerning the $n$-divisibility of class numbers
of the imaginary quadratic fields $\mathbb{Q}\left(\sqrt{1-4U^{n}}\right)$ and
$\mathbb{Q}\left(\sqrt{1-V^{n}}\right)$. These results will be used in the
proof of Theorem 1.1.

######  Theorem A.

If $n\geq 3$ is an odd integer, then for any integer $U\geq 2$ the ideal class
groups of the imaginary quadratic fields $\mathbb{Q}(\sqrt{1-4U^{n}})$ contain
an element of order $n$.

In 1978, Gross and Rohrlich gave the outline of a proof of this theorem (see
[6, Theorem 5.3 and Remark 2]). Their method of proof was based upon the
affine points on the Fermat curve $x^{p}+y^{p}=1$ over the imaginary quadratic
field $\mathbb{Q}(\sqrt{1-4U^{n}})$. Later, Louboutin [16] gave a complete
proof of this theorem using number theoretic technique. It follows from
Siegel’s theorem (see [20]) that for each integer $d>1$ there are at most
finitely many positive integers $U$ such that $1-4U^{n}=-dX^{2}$. This ensures
the infinitude of the above family of imaginary quadratic fields.

The $n$-divisibility of the class numbers of the family of imaginary quadratic
fields $\mathbb{Q}(\sqrt{1-V^{n}})$ was studied by Nagell [18, Theorem 25] for
any odd integers $V\geq 3$ and $n\geq 3$. Later, Murty [17, Theorems 1 and 2]
proved that the class group of the imaginary quadratic field
$\mathbb{Q}(\sqrt{1-V^{n}})$ has an element of order $n$ when either $V^{n}-1$
is square-free or its square part is $<V^{n/2}/8$ for any odd integers $V\geq
5$ and $n\geq 3$. However, it follows from the fact $DX^{2}+1=V^{n}$ has no
integer solution when both $V$ and $n$ are odd, except for
$(V,n)=(5,3),(7,3),(13,3)$ (see [7]), that the above conditions are no longer
required. This fact also confirms that there are infinitely many such
imaginary quadratic fields. Finally, this was elucidated by Cohn in [5,
Corollary 1] as follows:

######  Theorem B.

Assume that $n\geq 3$ and $V\geq 3$ are odd integers. Then the class number of
the imaginary quadratic field $\mathbb{Q}(\sqrt{1-V^{n}})$ is divisible by
$n$, except for $(n,V)=(5,3)$.

##  3\. The divisibility of the class number of
$\mathbb{Q}(\sqrt{p^{2}-\ell^{n}})$

Many special cases of the divisibility of the class number of the imaginary
quadratic field $\mathbb{Q}(\sqrt{x^{2}-y^{n}})$ have been studied with some
restrictions on $x,y$ and $n$. One of such restrictions is that $y$ is an odd
prime (see [3] and the references therein), and hence none of the known
results can be used to complete the proof of Theorem (1.1). Thus, we consider
a family of imaginary quadratic fields of the above form where $y$ is not a
prime and it will be a useful ingredient in the proof of Theorem 1.1. Here, we
mainly prove:

######  Theorem 3.1.

Let $\ell>1$ and $n>1$ be odd integers, and $p$ an odd prime such that
$\ell\equiv 3\pmod{4},~{}\gcd(\ell,p)=1$ and $p^{2}<\ell^{n}$. Assume that
$-d$ is the square-free part of $p^{2}-\ell^{n}$. If $p\not\equiv\pm
1\pmod{d}$, then the class number of $\mathbb{Q}(\sqrt{p^{2}-\ell^{n}})$ is
divisible by $n$.

Theorem 3.1 extends [3, Theorem 1.1], where the authors assumed that $\ell$ is
an odd prime. This primality condition on $\ell$ restricts us to apply [3,
Theorem 1.1] in the proof of Theorem 1.1. Further, following the proof of [3,
Theorem 1.2], we make the following remark.

######  Remark 3.1.

The family of imaginary quadratic fields discussed in Theorem 3.1 has
infinitely many members.

Now for $p=3,5$, the condition ‘$p\not\equiv\pm 1\pmod{d}$’ can be removed by
applying 3.7 except for $(\ell,n)=(3,3)$, and thus we have the following
straightforward corollary.

######  Corollary 3.1.

Let $\ell>1$ and $n>1$ be as in Theorem 3.1. For $x=3,5$ with
$\gcd(\ell,x)=1$, the class number of $\mathbb{Q}(\sqrt{x^{2}-\ell^{n}})$ is
divisible by $n$ except the case $(\ell,n)=(3,3)$.

The proof of Theorem 3.1 relies on the prominent result of Bilu, Hanrot and
Voutier [2] on existence of primitive divisors of Lehmer numbers.

###  3.1. Lehmer numbers and their primitive divisors

A pair $(\alpha,\beta)$ of algebraic integers is said to be a Lehmer pair if
$(\alpha+\beta)^{2}$ and $\alpha\beta$ are two non-zero coprime rational
integers, and $\alpha/\beta$ is not a root of unity. For a given positive
integer $n$, the Lehmer numbers correspond to the pair $(\alpha,\beta)$ are
defined as

$\mathcal{L}_{n}(\alpha,\beta)=\begin{cases}\dfrac{\alpha^{n}-\beta^{n}}{\alpha-\beta},&\text{
if $n$ is odd,}\vspace{1mm}\\\
\dfrac{\alpha^{n}-\beta^{n}}{\alpha^{2}-\beta^{2}},&\text{ if $n$ is
even}.\end{cases}$

It is known that all Lehmer numbers are non-zero rational integers. Two Lehmer
pairs $(\alpha_{1},\beta_{1})$ and $(\alpha_{2},\beta_{2})$ are said to be
equivalent if $\alpha_{1}/\alpha_{2}=\beta_{1}/\beta_{2}\in\\{\pm
1,\pm\sqrt{-1}\\}$. A prime divisor $p$ of $\mathcal{L}_{n}(\alpha,\beta)$ is
said to be primitive if
$p\nmid(\alpha^{2}-\beta^{2})^{2}\mathcal{L}_{1}(\alpha,\beta)\mathcal{L}_{2}(\alpha,\beta)\cdots\mathcal{L}_{n-1}(\alpha,\beta)$.
The following classical result was proved in [2, Theorem 1.4].

######  Theorem C.

The Lehmer number $\mathcal{L}_{n}(\alpha,\beta)$ has primitive divisors for
any integer $n>30$.

Given a Lehmer pair $(\alpha,\beta)$, let $a=(\alpha+\beta)^{2}$ and
$b=(\alpha-\beta)^{2}$. Then $\alpha=(\sqrt{a}\pm\sqrt{b})/2$ and
$\beta=(\sqrt{a}\mp\sqrt{b})/2$. The pair $(a,b)$ is called the parameters
corresponding to the Lehmer pair $(\alpha,\beta)$. The following lemma is
extracted from [22, Theorem 1].

######  Lemma 3.1.

Let $t$ be an odd integer such that $7\leq t\leq 29$. If the Lehmer numbers
$\mathcal{L}_{t}(\alpha,\beta)$ have no primitive divisor, then up to
equivalence, the parameters $(a,b)$ of the corresponding pair $(\alpha,\beta)$
are as follows:

  * (i)

$(a,b)=(1,-7),(1,-19),(3,-5),(5,-7),(13,-3),(14,-22)$, when $t=7$;

  * (ii)

$(a,b)=(5,-3),(7,-1),(7,-5)$, when $t=9$;

  * (iii)

$(a,b)=(1,-7)$, when $t=13$;

  * (iv)

$(a,b)=(7,-1),(10,-2)$, when $t=15$.

Let $F_{k}$ (resp. $L_{k}$) denote the $k$-th term in the Fibonacci (resp.
Lucas) sequence defined by $F_{0}=0,F_{1}=1$, and $F_{k+2}=F_{k}+F_{k+1}$
(resp. $L_{0}=2,L_{1}=1$, and $L_{k+2}=L_{k}+L_{k+1}$), where $k\geq 0$ is an
integer. The following lemma is a part of [2, Theorem 1.3].

######  Lemma 3.2.

For $p=3,5$, let the Lehmer numbers $\mathcal{L}_{p}(\alpha,\beta)$ have no
primitive divisor. Then up to equivalence, the parameters $(a,b)$ of the
corresponding pair $(\alpha,\beta)$ are:

  * (i)

For $p=3,(a,b)=\begin{cases}(1+t,1-3t)\text{ with }t\neq 1,\\\
(3^{k}+t,3^{k}-3t)\text{ with }t\not\equiv 0\pmod{3},(k,t)\neq(1,1);\\\
\end{cases}$

  * (ii)

For
$p=5,(a,b)=\begin{cases}(F_{k-2\varepsilon},F_{k-2\varepsilon}-4F_{k})\text{
with }k\geq 3,\\\ (L_{k-2\varepsilon},L_{k-2\varepsilon}-4L_{k})\text{ with
}k\neq 1;\end{cases}$

where $t\neq 0$ and $k\geq 0$ are any integers and $\varepsilon=\pm 1$.

###  3.2. Two important lemmas

Given an integer $D\equiv 0,1\pmod{4}$, assume that $h^{*}(D)$ is the class
number of binary quadratic primitive forms with discriminant $D$. Also for a
square-free integer $d$, let $h(d)$ denote the class number of
$\mathbb{Q}(\sqrt{d})$. Then we have (cf. [10, §16.13; p. 444] the following:

######  Lemma 3.3.

Let $d\equiv 2\pmod{4}$ be a square-free positive integer. Then
$h(-d)=h^{*}(-4d).$

The following lemma is a special case of [8, Theorem 6.2] when
$(D_{1},D_{2})=(1,-d)$.

######  Lemma 3.4.

Let $d>3$ and $\ell>1$ be integers such that $\gcd(\ell,2d)=1$. If the
equation

$x^{2}+dy^{2}=\ell^{z},~{}~{}x,y,z\in\mathbb{N},~{}\gcd(x,y)=1$

has a solution, then all the solutions $(x,y,z)$ of this equation can be
expressed as

$x+y\sqrt{-d}=\varepsilon(a+\mu b\sqrt{-d})^{t},~{}~{}z=st,$

where $a,b,s,t$ are positive integers satisfying

$a^{2}+db^{2}=\ell^{s},~{}~{}\gcd(a,b)=1\text{ and }s\mid h^{*}(-4d)$

and $\varepsilon,\mu\in\\{-1,1\\}$.

###  3.3. Proof of Theorem 3.1

Let $d$ be the square-free part of $\ell^{n}-p^{2}$. Then
$p^{2}-\ell^{n}=-dr^{2}$ for some $r\in\mathbb{N}$, and thus $(x,y,z)=(p,r,n)$
is a positive integer solution of the equation

$x^{2}+dy^{2}=\ell^{z},~{}\gcd(x,y)=1.$

Thus by Lemma 3.4, we get

$p+r\sqrt{-d}=\varepsilon(a+\mu
b\sqrt{-d})^{t},~{}~{}\varepsilon,\mu\in\\{-1,1\\}$ (3.1)

with

$n=st,~{}~{}s,t\in\mathbb{N},$ (3.2)

where $a$ and $b$ are positive integers satisfying

$a^{2}+db^{2}=\ell^{s},~{}~{}\gcd(a,b)=1$ (3.3)

and

$s\mid h^{*}(-4d).$ (3.4)

Since $\ell\equiv 3\pmod{4}$ and $n$ is odd, so that $p^{2}-\ell^{n}=-dr^{2}$
gives $d\equiv 2\pmod{4}$ and $r$ is odd. Also both $s$ and $t$ are odd as $n$
is odd. Further reading (3.3) modulo $4$, we get $a^{2}+2b^{2}\equiv
3\pmod{4}$ as $\ell\equiv 3\pmod{4}$ and $s$ is odd, which ensures that both
$a$ and $b$ are odd.

We now equate the real parts from both sides in (3.1) to get

$p=\varepsilon
a\sum^{\frac{t-1}{2}}_{j=0}\binom{t}{2j}a^{t-2j-1}(-db^{2})^{j}.$ (3.5)

This implies $a\mid p$ and thus $a=1,p$. If $a=1$, then it becomes

$\sum^{\frac{t-1}{2}}_{j=0}\binom{t}{2j}(-db^{2})^{j}=p\varepsilon=\pm p.$
(3.6)

Reading (3.6) modulo $d$, we get $p\equiv\pm 1\pmod{d}$, which contradicts to
the assumption. Therefore $a=p$ and thus (3.5) becomes

$\sum^{\frac{t-1}{2}}_{j=0}\binom{t}{2j}p^{t-2j-1}(-db^{2})^{j}=\varepsilon=\pm
1.$ (3.7)

As $a=p$, so that (3.1) reduces to

$p+r\sqrt{-d}=\varepsilon(p+\mu
b\sqrt{-d})^{t},~{}~{}\varepsilon,\mu\in\\{-1,1\\}.$ (3.8)

We now assume that $\alpha=\mu b\sqrt{-d}+p$ and $\beta=\mu b\sqrt{-d}-p$.
Then both $\alpha$ and $\beta$ are algebraic integers. Clearly,
$(\alpha+\beta)^{2}=-4db^{2}$ and $\alpha\beta=-p^{2}-db^{2}=-\ell^{s}$ (by
(3.3)) are coprime rational integers. Furthermore, it follows from the
following identity

$\frac{4db^{2}}{\ell^{s}}=\frac{(\alpha+\beta)^{2}}{\alpha\beta}=\frac{\alpha}{\beta}+\frac{\beta}{\alpha}+2$

that

$\ell^{s}\left(\frac{\alpha}{\beta}\right)^{2}+2(\ell^{s}-2db^{2})\frac{\alpha}{\beta}+\ell^{s}=0.$

Since $\ell>1$ and
$\gcd(\ell^{s},2(\ell^{s}-2db^{2}))=\gcd(\ell^{s},4db^{2})=\gcd(p^{2}+db^{2},4db^{2})=1$,
so that the last equation shows that $\dfrac{\alpha}{\beta}$ is not a root of
unity. Therefore $(\alpha,\beta)$ is a Lehmer pair with parameters
$(-4db^{2},4p^{2})$ and thus the corresponding Lehmer number for $t$ is

$\mathcal{L}_{t}(\alpha,\beta)=\frac{\alpha^{t}-\beta^{t}}{\alpha-\beta}$

as $t$ is odd. Utilizing (3.8), we get

$|\mathcal{L}_{t}(\alpha,\beta)|=1.$

This confirms that the Lehmer number $\mathcal{L}_{t}(\alpha,\beta)$ has no
primitive divisor, and hence Theorem C and Lemma 3.1 (utilizing the fact that
$(-4db^{2},4p^{2})$ is the parameters) ensure that $t\in\\{1,3,5\\}$.

In case of $t=5$, we get by Lemma 3.2 that $-4db^{2}=F_{k-2\varepsilon}$ or
$-4db^{2}=L_{k-2\varepsilon}$. Clearly, none of these is possible.

Finally for $t=3$, (3.7) implies that $p^{2}-3db^{2}=\pm 1$. Reading it modulo
$4$, we see that ‘$+$’ sign is not possible, and thus $p^{2}-3db^{2}=-1$. This
is not possible by reading it modulo $3$.

Therefore $t=1$, and thus (3.2) and (3.4) together imply that $n\mid
h^{*}(-4d)$. Thus, we complete the proof by Lemma 3.3.

##  4\. Proof of Theorem 1.1

We first fix an odd integer $n\geq 3$. We now define the set

$\mathcal{N}_{n}=\left\\{k\in\mathbb{N}:n\mid h(1-4k^{n})\right\\}.$

Then by Theorem A the set $\mathcal{N}_{n}$ is infinite.

Now for any $k\in\mathcal{N}_{n}$, we set $d=4(1-4k^{n})^{n}$. Then
$\mathbb{Q}(\sqrt{d})=\mathbb{Q}(\sqrt{1-4k^{n}})$ as $n$ is odd. Thus, by
Theorem A there are infinitely many such $d$ satisfying $n\mid h(d)$. In other
words,
$\mathcal{F}(\mathcal{N}_{n})=\\{\mathbb{Q}(\sqrt{1-4k^{n}}):k\in\mathcal{N}_{n}\\}$
is an infinite set.

Now we assume that $U=4k^{n}-1$ with $k\in\mathcal{N}_{n}$. Then
$1-4U^{n}=1-4(4k^{n}-1)^{n}=4(1-4k^{n})^{n}+1=d+1$, and thus by Theorem A, we
have $n\mid h(d+1)$.

Again for $k\in\mathcal{N}_{n}$, let us assume that $V=4k^{n}-1$. Then $V\geq
3$ and is odd, and thus by Theorem B, we get $n\mid h(1-V^{n})$. Since
$4(1-V^{n})=4-4(4k^{n}-1)^{n}=4+4(1-4k^{n})^{n}=d+4$ and
$\mathbb{Q}(\sqrt{4(1-V^{n})})=\mathbb{Q}(\sqrt{(1-V^{n})})$, so that $n\mid
h(d+4)$.

Finally for any $k\in\mathcal{N}_{n}$, let $\ell=4k^{n}-1$. Then $\ell\equiv
3\pmod{4}$ and hence by utilizing Theorem 3.1, we have $n\mid
h(p^{2}-\ell^{n})$ for any odd prime $p$ satisfying $p\not\equiv\pm
1\pmod{d}$. Here, $d$ is the square-free part of $\ell^{n}-p^{2}$. Now
$4(p^{2}-\ell^{n})=4p^{2}-4\ell^{n}=4p^{2}-4(4k^{n}-1)^{n}=d+4p^{2}$, which
implies that
$\mathbb{Q}(\sqrt{d+4p^{2}})=\mathbb{Q}(\sqrt{4(p^{2}-\ell^{n})})$, and thus
$n\mid h(d+4p^{2})$. This completes the proof of Theorem 1.1.

##  5\. Concluding remarks

In [23], Xie and Chao studied Conjecture 1.1 and proved the following result
using Yamamoto’s [24] construction.

######  Theorem D.

For any odd positive integer $n$ and any positive integer $m$, there are
infinitely many pairs of imaginary fields $\mathbb{Q}(\sqrt{d})$ and
$\mathbb{Q}(\sqrt{d+m})$ whose class groups have an element of order $n$
respectively.

Theorem D can be viewed as a weaker variant of a generalization of Conjecture
1.1. For $m=1$, it provides a generalization of the main result of [15] though
[23] appeared before [15]. In other words, Theorem D gives a complete proof of
the following generalization of Conjecture 1.1 for $m=1$ and a proof of a
weaker version of the same for $m\geq 2$.

######  Conjecture 5.1.

For any odd integer $n\geq 3$ and any integer $m\geq 1$, there is an infinite
family of $m+1$ successive imaginary (or real) quadratic fields

$\mathbb{Q}(\sqrt{d}),\mathbb{Q}(\sqrt{d+1}),\cdots,\mathbb{Q}(\sqrt{d+m})$

with $d\in\mathbb{Z}$ whose class numbers are all divisible by $n$.

Theorem 1.1 offers a constructive proof of Conjecture 5.1 for $m=1$. This
theorem also offers a proof of a weaker version of Conjecture 5.1 for $m=4$,
which has missed the families of imaginary quadratic fields
$\mathbb{Q}(\sqrt{d+2})$ and $\mathbb{Q}(\sqrt{d+3})$ from the complete proof.
When $m=4p^{2}$ with $p$ an odd prime, Theorem 1.1 presents a proof of a
weaker version of Conjecture 5.1. We complete this paper by the following
remark.

######  Remark 5.1.

For a given positive integer $m$, let $p_{m}$ denote the largest prime less
than or equal to $m$ and $\pi(m)$ the prime-counting function. Then for a
given positive odd integer $n$, our construction gives an infinite family of
at least $(\pi(m)+2)$-tuples of imaginary quadratic fields,

$\left(\mathbb{Q}(\sqrt{d}),\mathbb{Q}(\sqrt{d+1}),\mathbb{Q}(\sqrt{d+4}),\mathbb{Q}(\sqrt{d+36}),\cdots\mathbb{Q}(\sqrt{d+4p_{m}^{2}})\right)$

with $d\in\mathbb{Z}$ whose class numbers are all divisible by $n$.

## Acknowledgements

The author is grateful to Professor K. Chakraborty and Professor Yasuhiro
Kishi for their valuable comments on the paper. The author is thankful to
Professor Y. Iizuka for providing a copy of [12]. The author gratefully
acknowledges the anonymous referee for his/her valuable remarks that immensely
improved the results as well as the presentation of the paper. This work was
supported by SERB-NPDF (PDF/2017/001958), Govt. of India.

## References

  * [1] N. C. Ankeny and S. Chowla, On the divisibility of the class number of quadratic fields, Pacific J. Math. 5 (1955), 321–324.
  * [2] Y. Bilu, G. Hanrot and P. M. Voutier, Existence of primitive divisors of Lucas and Lehmer numbers (with an appendix by M. Mignotte), J. Reine Angew. Math. 539 (2001), 75–122.
  * [3] K. Chakraborty, A. Hoque, Y. Kishi and P. P. Pandey, Divisibility of the class numbers of imaginary quadratic fields, J. Number Theory 185 (2018), 339–348.
  * [4] J. Chattopadhyay and S. Muthukrishnan, On the simultaneous $3$-divisibility of class numbers of triples of imaginary quadratic fields, Acta Arith. 197 (2021), no. 1, 105–110.
  * [5] J. H. E. Cohn, On the Diophantine equation $x^{n}=Dy^{2}+1$, Acta Arith. 106 (2003), no. 1, 73–83.
  * [6] B. H. Gross and D. E. Rohrlich, Some results on the Mordell-Weil group of the Jacobian of the Fermat curve, Invent. Math. 44 (1978), 201–224.
  * [7] E. Herrmann, I. Járási and A. Pethő, Note on J. H. E. Cohn’s paper “The Diophantine equation $x^{n}=Dy^{2}+1$”, Acta Arith. 113 (2004), no. 1, 69–76.
  * [8] C. Heuberger and M. Le, On the generalized Ramanujan-Nagell equation $x^{2}+D=p^{Z}$, J. Number Theory 78 (1999), no. 2, 312–331.
  * [9] A. Hoque and K. Chakraborty, Divisibility of class numbers of certain families of quadratic fields, J. Ramanujan Math. Soc. 34 (2019), no. 3, 281–289.
  * [10] L. K. Hua, Introduction to number theory, Springer-Verlag, New York, 1982.
  * [11] Y. Iizuka, On the class number divisibility of pairs of imaginary quadratic fields, J. Number Theory 184 (2018), 122–127.
  * [12] Y. Iizuka, Y. Konomi and S. Nakano, An application of the arithmetic of elliptic curves to the class number problem for quadratic fields, Tokyo J. Math. (2021). doi: 10.3836/tjm/1502179314
  * [13] T. Komatsu, An infinite family of pairs of quadratic fields $\mathbb{Q}(\sqrt{D})$ and $\mathbb{Q}(\sqrt{mD})$ whose class numbers are both divisible by $3$, Acta Arith. 104 (2002), 129–136.
  * [14] T. Komatsu, An infinite family of pairs of imaginary quadratic fields with ideal classes of a given order, Int. J. Number Theory 13 (2017), no. 2, 253–260.
  * [15] S. Krishnamoorthy and S. Pasupulati, Note on the $p$-divisibility of class numbers of an infinite family of imaginary quadratic fields, Glasgow Math. J. (2021). doi: 10.1017/S001708952100015X
  * [16] S. R. Louboutin, On the divisibility of the class number of imaginary quadratic number fields, Proc. Amer. Math. Soc. 137 (2009), no. 12, 4025–4028.
  * [17] M. R. Murty, The ABC conjecture and exponents of class groups of quadratic fields, Contemp. Math. 210 (1998), 85–95.
  * [18] T. Nagell, Über die Klassenzahl imaginär quadratischer, Zählkörper, Abh. Math. Sem. Univ. Hambg. 1 (1922), 140–150.
  * [19] A. Scholz, Über die Beziehung der Klassenzahlen quadratischer Körper zueinander, J. Reine Angew. Math. 166 (1932), 201–203.
  * [20] C. L. Siegel, Uber einige Anwendungen Diophantischer Approximationen, Abh. Preuss. Akad. Wiss. Phys. Math. Kl. 1 (1929), 1-70; Ges. Abh., Band 1, 209–266.
  * [21] K. Soundararajan, Divisibility of class numbers of imaginary quadratic fields, J. London Math. Soc. 61 (2000), 681–690.
  * [22] P. M. Voutier, Primitive divisors of Lucas and Lehmer sequences, Math. Comp. 64 (1995), 869–888.
  * [23] C. -F. Xie and C. F. Chao, On the divisibility of class numbers of imaginary quadratic fields $\left(\mathbb{Q}(\sqrt{D}),\mathbb{Q}(\sqrt{D+m})\right)$, Ramanujan J. 53 (2020), 517–528.
  * [24] Y. Yamamoto, On unramified Galois extensions of quadratic number fields, Osaka J. Math. 7 (1970), 57–76.