Note

Click here to download the full example code

# 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

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

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

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:

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).

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.

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.

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

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.17637844761413474
|1100>: 0.06855956371224507
|0101>: 0.002056097258972289
|1111>: 0.008342946399881935
|2000>: 0.010312945253440309
```

## Equally squeezed inputs¶

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

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

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.17637844761413474
```

**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.06855956371214537
0.06855956371224507
```

**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.002056097258972814
0.002056097258972289
```

**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.008342946399869364
0.008342946399881935
```

**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.010312945253440309
```

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

A Fock backend such as

`'fock'`

instead of the Gaussian backend.Different beamsplitter and rotation parameters.

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.

**Total running time of the script:** ( 0 minutes 0.017 seconds)

## Contents

## Downloads

## Related tutorials