# Gaussian boson sampling tutorial¶

*Section author: Josh Izaac <josh@xanadu.ai>*

Note

This tutorial is also available in the form of an interactive Jupyter Notebook `GaussianBosonSampling.ipynb`

In this short tutorial, we will walk through the application of the Gaussian boson sampling Blackbird example provided in the Gaussian boson sampling section of the Quantum algorithms page. Be sure to read through that section to familiarise yourself with the concepts behind Gaussian boson sampling, as well as the introductory teleportation tutorial, before attempting this tutorial. It is also strongly recommended to read the boson sampling tutorial, which introduces a few concepts also used in Gaussian boson sampling.

## Circuit construction and simulation¶

A 4-mode Gaussian boson sampling circuit is given by

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, create a new text file with the name `gaussian_boson_sampling.py`

, and open it with an editor. Copy and paste the following Strawberry Fields code:

```
#!/usr/bin/env python3
import numpy as np
import strawberryfields as sf
from strawberryfields.ops import *
# initialise the engine and register
eng, q = sf.Engine(4)
# 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]
])
with eng:
# 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
# end circuit
# run the engine
state = eng.run('gaussian')
# 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 = state.fock_prob(i)
print("|{}>: {}".format("".join(str(j) for j in i), prob))
```

Note

This example program is included with Strawberry Fields as `examples/gaussian_boson_sampling.py`

.

A couple of things to note in this particular example:

- To prepare the input single mode squeezed vacuum state \(\ket{z} = 1\), we apply a squeezing operator to each modes (initially in the vacuum state). This is done via the
`Sgate`

operator.

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 tutorial. As a result, the decomposed beamsplitter and rotation gate parameters will also be identical.

- Unlike the boson sampling tutorial, the lack of Fock states means we can now use the Gaussian simulator backend, indicated with the argument
`'gaussian'`

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

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.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 Gaussian states. 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.

The simulation can be run by executing the file from the command line, resulting in the following output:

```
$ python3 gaussian_boson_sampling.py
|0000>: 0.176378447614135
|1100>: 0.0685595637122246
|0101>: 0.002056097258977398
|1111>: 0.00834294639986785
|2000>: 0.01031294525345511
```

## Equally squeezed inputs¶

Recall that, as shown in Ref. [13] and reviewed in the Gaussian boson sampling section of the quantum algorithms page, the formula for calculating the output Fock state probability,

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

## Analysis¶

To analyse the results, it is convenient to now move to a Python console or interactive environment, such as iPython or Jupyter Notebook. In the following, Python input will be specified with the prompt `>>>`

, and output will follow.

First, import some useful libraries, such as NumPy, as well as the `multi_dot`

and `block_diag`

functions from NumPy and SciPy respectively.:

```
>>> import numpy as np
>>> from numpy.linalg import multi_dot
>>> from scipy.linalg import block_diag
```

As we are using the interferometer as in the Boson sampling tutorial, we do not have to recalculate the unitary, we can use the same \(U\) computed there, or copy the `U`

defined in the file `gaussian_boson_sampling.py`

:

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 following custom function:

```
>>> from itertools import permutations
>>> from scipy.special import factorial
>>> def Haf(M):
... n=len(M)
... m=int(n/2)
... haf=0.0
... for i in permutations(range(n)):
... prod=1.0
... for j in range(m):
... prod*=M[i[2*j],i[2*j+1]]
... haf+=prod
... return haf/(factorial(m)*(2**m))
```

Note

This function is written based on the basic definition of the hafnian,

Notice that this function counts each term in the definition multiple times, and renormalizes to remove the multiple counts by dividing by a factor \(n!2^n\). **This function is extremely slow!**

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\),

```
>>> B[:,[0,1]][[0,1]]
[[-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:>>> 1/np.cosh(1)**4 0.17637844761413471

Compare this to the Strawberry Fields result

`0.176378447614135`

**Measuring**\(\ket{1,1,0,0}\)**at the output**>>> B = (np.dot(U,U.T) * np.tanh(1))[:, [0,1]][[0,1]] >>> np.abs(Haf(B))**2 / np.cosh(1)**4 0.068559563712223492

Compare this to the Strawberry Fields result

`0.0685595637122246`

**Measuring**\(\ket{0,1,0,1}\)**at the output**>>> B = (np.dot(U,U.T) * np.tanh(1))[:, [1,3]][[1,3]] >>> np.abs(Haf(B))**2 / np.cosh(1)**4 0.0020560972589773979

Compare this to the Strawberry Fields result

`0.002056097258977398`

**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)) >>> np.abs(Haf(B))**2 / np.cosh(1)**4 0.0083429463998674833

Compare this to the Strawberry Fields result

`0.00834294639986785`

**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]] >>> np.abs(Haf(B))**2 / (2*np.cosh(1)**4) 0.010312945253454881

Compare this to the Strawberry Fields result

`0.01031294525345511`

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

Exercises

Repeat this tutorial with

- A Fock backend such as Fock simulator backend, 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.