Gaussian boson sampling and the Hafnian

“If you need to wait exponential time for [your single photon sources to emit simultaneously], then there would seem to be no advantage over classical computation. This is the reason why so far, boson sampling has only been demonstrated with 3-4 photons. When faced with these problems, until recently, all we could do was shrug our shoulders.” - Scott Aaronson

In this tutorial, we will walk through the application of the Gaussian boson sampling. It is also strongly recommended to read the boson sampling tutorial, which introduces a few concepts also used in Gaussian boson sampling.

Background theory

While boson sampling allows the experimental implementation of a sampling problem that is countably hard classically, one of the main issues it has in experimental setups is one of scalability, due to its dependence on an array of simultaneously emitting single photon sources. Currently, most physical implementations of boson sampling make use of a process known as spontaneous parametric down-conversion to generate the single-photon source inputs. However, this method is non-deterministic — as the number of modes in the apparatus increases, the average time required until every photon source emits a simultaneous photon increases exponentially.

In order to simulate a deterministic single-photon source array, several variations on boson sampling have been proposed; the most well-known being scattershot boson sampling [1]. However, a recent boson sampling variation by Hamilton et al. [2] mitigates the need for single photon Fock states altogether, by showing that incident Gaussian states — in this case, single mode squeezed states — can produce problems in the same computational complexity class as boson sampling. Even more significantly, this negates the scalability problem with single photon sources, as single mode squeezed states can be deterministically generated simultaneously.

The Gaussian boson sampling scheme remains, on initial observation, quite similar to that of boson sampling:

• \(N\) single mode squeezed states \(\ket{z}\), with squeezing parameter \(z=re^{i\phi}\), enter an \(N\) mode linear interferometer described by unitary \(U\) simultaneously.

• Each output mode of the interferometer (denoted by state \(\ket{\psi'}\)) is then measured in the Fock basis, \(\bigotimes_i n_i\ket{n_i}\bra{n_i}\).

Without loss of generality, we can absorb the squeezing phase parameter \(\phi\) into the interferometer, and set \(\phi=0\) for convenience.

Using phase space methods, Hamilton et al. [2] showed that the probability of measuring a Fock state containing only 0 or 1 photons per mode is given by

\[\left|\left\langle{n_1,n_2,\dots,n_N}\middle|{\psi'}\right\rangle\right|^2 = \frac{\left|\text{Haf}[(U(\bigoplus_i\tanh(r_i))U^T)]_{st}\right|^2}{\prod_{i=1}^N \cosh(r_i)}\]

i.e., the sampled single photon probability distribution is proportional to the hafnian of a submatrix of \(U(\bigoplus_i\tanh(r_i))U^T\), dependent upon the output covariance matrix.

Note

The hafnian of a matrix is defined by

\[\text{Haf}(A) = \frac{1}{n!2^n}\sum_{\sigma \in S_{2N}}\prod_{i=1}^N A_{\sigma(2i-1)\sigma(2i)}\]

where \(S_{2N}\) is the set of all permutations of \(2N\) elements. In graph theory, the hafnian calculates the number of perfect matchings in an arbitrary graph with adjacency matrix \(A\).

Compare this to the permanent, which calculates the number of perfect matchings on a bipartite graph - the hafnian turns out to be a generalization of the permanent, with the relationship

\[\begin{split}\text{Per(A)} = \text{Haf}\left(\left[\begin{matrix} 0&A\\ A^T&0 \end{matrix}\right]\right)\end{split}\]

As any algorithm that could calculate (or even approximate) the hafnian could also calculate the permanent - a #P-hard problem - it follows that calculating or approximating the hafnian must also be a classically hard problem.

Circuit construction and simulation

As with the boson sampling problem, the multimode linear interferometer can be decomposed into two-mode beamsplitters (BSgate) and single-mode phase shifters (Rgate) [3], allowing for a straightforward translation into a CV quantum circuit.

For example, in the case of a 4 mode interferometer, with arbitrary \(4\times 4\) unitary \(U\), the CV quantum circuit for Gaussian boson sampling is given by

In the above, the single mode squeeze states all apply identical squeezing \(z=r\), the parameters of the beamsplitters and the rotation gates determine the unitary \(U\), and finally the detectors perform Fock state measurements on the output modes. As with boson sampling, for \(N\) input modes, we must have a minimum of \(N+1\) columns in the beamsplitter array [4].

Simulating this circuit using Strawberry Fields is easy; we can simply read off the gates from left to right, and convert it into the Blackbird circuit language.

To begin, we create the boson sampling quantum program using Strawberry Fields:

import numpy as np

# set the random seed
np.random.seed(42)

# import Strawberry Fields
import strawberryfields as sf
from strawberryfields.ops import *

Unlike the Boson sampling and the permanent, we will directly apply a Gaussian unitary to the circuit using the Interferometer operation. First, we must define the unitary matrix we would like to embed in the circuit:

# define the linear interferometer
U = np.array([
 [ 0.219546940711-0.256534554457j, 0.611076853957+0.524178937791j,
    -0.102700187435+0.474478834685j,-0.027250232925+0.03729094623j],
 [ 0.451281863394+0.602582912475j, 0.456952590016+0.01230749109j,
    0.131625867435-0.450417744715j, 0.035283194078-0.053244267184j],
 [ 0.038710094355+0.492715562066j,-0.019212744068-0.321842852355j,
    -0.240776471286+0.524432833034j,-0.458388143039+0.329633367819j],
 [-0.156619083736+0.224568570065j, 0.109992223305-0.163750223027j,
    -0.421179844245+0.183644837982j, 0.818769184612+0.068015658737j]
])

We can use this to now construct the circuit:

# create the 4 mode Strawberry Fields program
gbs = sf.Program(4)

with gbs.context as q:
    # prepare the input squeezed states
    S = Sgate(1)
    S | q[0]
    S | q[1]
    S | q[2]
    S | q[3]

    # linear interferometer
    Interferometer(U) | q

A couple of things to note in this particular example:

1. To prepare the input single mode squeezed vacuum state \(\ket{z}\) where \(z = 1\), we apply a squeezing gate Sgate to each of the modes (initially in the vacuum state).

2. Next we apply the linear interferometer to all four modes, using the decomposition operator Interferometer, and the unitary matrix U. This operator decomposes the unitary matrix representing the linear interferometer into single mode rotation gates Rgate, and two-mode beamsplitters BSgate. After applying the interferometer, we will denote the output state by \(\ket{\psi'}\).

   Note

   You can view the decomposed beamsplitters and rotation gates which correspond to the linear interferometer U by calling eng.print_applied() after running the engine.

   Note

   The interferometer applied here is identical to that from the Boson sampling and the permanent. As a result, the decomposed beamsplitter and rotation gate parameters will also be identical.

3. We are not performing Fock measurements at the output; this is to ensure the state is preserved, so we can extract the joint Fock state probabilities after the beamsplitter array.

Unlike the boson sampling tutorial, the lack of Fock states means we can now use the Numpy-based 'gaussian' backend, along with a 4-mode register. The Gaussian backend is perfectly suited for simulation of Gaussian boson sampling, as all initial states are Gaussian, and all the required operators transform Gaussian states to other Gaussian states.

Now that the program is created, we can initialize an engine, and execute the program on one of the built-in Strawberry Fields state simulators.

eng = sf.Engine(backend="gaussian")
results = eng.run(gbs)

The state method fock_prob() accepts a list or a tuple containing the Fock state to be measured and returns the probability of that measurement. For example, [1,2,0,1] represents the measurement resulting in the detection of 1 photon at mode q[0] and mode q[3], and 2 photons at mode q[1], and would return the value

\[\text{prob}(1,2,0,1) = \left|\braketD{1,2,0,1}{\psi'}\right|^2\]

The Fock state method all_fock_probs(), used previously to return all Fock state probabilities as an array, is not supported by the Gaussian backend. This is because computing the Fock probabilities of states in the Gaussian representation has exponential scaling - while this is fine for computing particular Fock basis probabilities, it becomes computationally demanding to return all Fock state probabilities using the Gaussian backend.

Let’s use fock_prob() in a for-loop to extract the probabilities of measuring various Fock states:

# Fock states to measure at output
measure_states = [[0,0,0,0], [1,1,0,0], [0,1,0,1], [1,1,1,1], [2,0,0,0]]

# extract the probabilities of calculating several
# different Fock states at the output, and print them to the terminal
for i in measure_states:
    prob = results.state.fock_prob(i)
    print("|{}>: {}".format("".join(str(j) for j in i), prob))

Out:

|0000>: 0.17637844761413482
|1100>: 0.0685595637122452
|0101>: 0.00205609725897229
|1111>: 0.008342946399881956
|2000>: 0.010312945253440318

Equally squeezed inputs

As shown earlier, the formula for calculating the output Fock state probability in Gaussian boson sampling is given by

\[\left|\left\langle{n_1,n_2,\dots,n_N}\middle|{\psi'}\right\rangle\right|^2 = \frac{\left|\text{Haf}[(U\bigoplus_i\tanh(r_i)U^T)]_{st}\right|^2}{n_1!n_2!\cdots n_N! \cosh(r_i)}\]

where \(U\) is the rotation/beamsplitter unitary transformation on the input and output mode annihilation and creation operators.

However, in this particular example, we are using the same squeezing parameter, \(z=r\), for all input states - this allows us to simplify this equation. To start with, the hafnian expression simply becomes \(\text{Haf}[(UU^T\tanh(r))]_{st}\), removing the need for the tensor sum.

Thus, we have

\[\left|\left\langle{n_1,n_2,\dots,n_N}\middle|{\psi'}\right\rangle\right|^2 = \frac{\left|\text{Haf}[(UU^T\tanh(r))]_{st}\right|^2}{n_1!n_2!\cdots n_N!\cosh^N(r)}.\]

Now that we have the interferometer unitary transformation \(U\), as well as the ‘experimental’ results, let’s compare the two, and see if the Gaussian boson sampling result in the case of equally squeezed input modes, agrees with the Strawberry Fields simulation probabilities.

Calculating the hafnian

Before we can calculate the right hand side of the Gaussian boson sampling equation, we need a method of calculating the hafnian. Since the hafnian is classically hard to compute, it is not provided in either NumPy or SciPy, so we will use The Walrus library, installed alongside Strawberry Fields:

from thewalrus import hafnian as haf

Now, for the right hand side numerator, we first calculate the submatrix \([(UU^T\tanh(r))]_{st}\):

B = (np.dot(U, U.T) * np.tanh(1))

Unlike the boson sampling case, in Gaussian boson sampling, we determine the submatrix by taking the rows and columns corresponding to the measured Fock state. For example, to calculate the submatrix in the case of the output measurement \(\left|{1,1,0,0}\right\rangle\),

print(B[:, [0, 1]][[0, 1]])

Out:

[[-0.10219728+0.32633851j  0.55418347+0.28563583j]
 [ 0.55418347+0.28563583j -0.10505237+0.32960794j]]

Comparing to Strawberry Fields

Now that we have a method for calculating the hafnian, let’s compare the output to that provided by Strawberry Fields.

Measuring \(\ket{0,0,0,0}\) at the output

This corresponds to the hafnian of an empty matrix, which is simply 1:

print(1 / np.cosh(1) ** 4)

print(results.state.fock_prob([0, 0, 0, 0]))

Out:

0.1763784476141347
0.17637844761413482

Measuring \(\ket{1,1,0,0}\) at the output

B = (np.dot(U, U.T) * np.tanh(1))[:, [0, 1]][[0, 1]]
print(np.abs(haf(B)) ** 2 / np.cosh(1) ** 4)

print(results.state.fock_prob([1, 1, 0, 0]))

Out:

0.0685595637121454
0.0685595637122452

Measuring \(\ket{0,1,0,1}\) at the output

B = (np.dot(U, U.T) * np.tanh(1))[:, [1, 3]][[1, 3]]
print(np.abs(haf(B)) ** 2 / np.cosh(1) ** 4)

print(results.state.fock_prob([0, 1, 0, 1]))

Out:

0.0020560972589728135
0.00205609725897229

Measuring \(\ket{1,1,1,1}\) at the output

This corresponds to the hafnian of the full matrix \(B=UU^T\tanh(r)\):

B = (np.dot(U, U.T) * np.tanh(1))
print(np.abs(haf(B)) ** 2 / np.cosh(1) ** 4)

print(results.state.fock_prob([1, 1, 1, 1]))

Out:

0.008342946399869365
0.008342946399881956

Measuring \(\ket{2,0,0,0}\) at the output

Since we have two photons in mode q[0], we take two copies of the first row and first column, making sure to divide by \(2!\):

B = (np.dot(U, U.T) * np.tanh(1))[:, [0, 0]][[0, 0]]
print(np.abs(haf(B)) ** 2 / (2 * np.cosh(1) ** 4))

print(results.state.fock_prob([2, 0, 0, 0]))

Out:

0.010312945253479155
0.010312945253440318

The Strawberry Field simulation results agree (with almost negligible numerical error) to the expected result from the Gaussian boson sampling equation!

Warning

Repeat this tutorial with

1. A Fock backend such as 'fock' instead of the Gaussian backend.

2. Different beamsplitter and rotation parameters.

3. Input states with differing squeezed values \(r_i\). You will need to modify the code to take into account the fact that the output covariance matrix determinant must now be calculated.

References

1

A. P. Lund, A. Laing, S. Rahimi-Keshari, T. Rudolph, J. L. O’Brien, and T. C. Ralph. Boson sampling from a gaussian state. Physical Review Letters, 113:100502, Sep 2014. doi:10.1103/PhysRevLett.113.100502.

2(1,2)

Craig S. Hamilton, Regina Kruse, Linda Sansoni, Sonja Barkhofen, Christine Silberhorn, and Igor Jex. Gaussian boson sampling. Physical Review Letters, 119:170501, Oct 2017. arXiv:1612.01199, doi:10.1103/PhysRevLett.119.170501.

3

Michael Reck, Anton Zeilinger, Herbert J. Bernstein, and Philip Bertani. Experimental realization of any discrete unitary operator. Physical Review Letters, 73(1):58–61, Jul 1994. doi:10.1103/physrevlett.73.58.

4

William R Clements, Peter C Humphreys, Benjamin J Metcalf, W Steven Kolthammer, and Ian A Walsmley. Optimal design for universal multiport interferometers. Optica, 3(12):1460–1465, 2016. doi:10.1364/OPTICA.3.001460.

demos/run_gaussian_boson_sampling

 

Download Python script

 

Download Notebook