# Boson sampling tutorial¶

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

In this short tutorial, we will walk through the application of the boson sampling Blackbird example provided in the Boson sampling section of the Quantum algorithms page. Be sure to read through that section to familiarize yourself with the concepts behind boson sampling, as well as the introductory teleportation tutorial, before attempting this tutorial.

## Circuit construction and simulation¶

A 4-mode 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 `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)
with eng:
# prepare the input fock states
Fock(1) | q[0]
Fock(1) | q[1]
Vac | q[2]
Fock(1) | q[3]
# rotation gates
Rgate(0.5719)
Rgate(-1.9782)
Rgate(2.0603)
Rgate(0.0644)
# beamsplitter array
BSgate(0.7804, 0.8578) | (q[0], q[1])
BSgate(0.06406, 0.5165) | (q[2], q[3])
BSgate(0.473, 0.1176) | (q[1], q[2])
BSgate(0.563, 0.1517) | (q[0], q[1])
BSgate(0.1323, 0.9946) | (q[2], q[3])
BSgate(0.311, 0.3231) | (q[1], q[2])
BSgate(0.4348, 0.0798) | (q[0], q[1])
BSgate(0.4368, 0.6157) | (q[2], q[3])
# end circuit
# run the engine
state = eng.run('fock', cutoff_dim=7)
# extract the joint Fock probabilities
probs = state.all_fock_probs()
# print the joint Fock state probabilities
print(probs[1,1,0,1])
print(probs[2,0,0,1])
```

Note

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

.

A couple of things to note in this particular example:

- We initialize the input state \(\ket{\psi} = \ket{1,1,0,1}\), by creating a single photon Fock state in modes 0, 1 and 3. This is done via the
`Fock`

operator. Mode 2 is initialized as a vacuum state using the`Vac`

operator (a shortcut to`Vacuum`

). This is*optional*- modes are automatically created in the vacuum state upon engine initialization.

- Next we apply the rotation gates,
`Rgate`

, to each mode. The resulting rotation in the phase space occurs in the anticlockwise direction, with angle \(\phi\).

- Finally, we apply the array of beamsplitters, using the
`BSgate`

operator, with arguments`(theta,phi)`

.- The transmission amplitude is then given by \(t=\cos\theta\)
- The reflection amplitude is given by \(r=e^{i\phi}\sin{\theta}\)

- The rotation gate and beamsplitter parameters have been chosen at random, leading to an arbitrary unitary \(U\) acting on the input modes annihilation and creation operators. After leaving the beamsplitter array, we will denote the output state by \(\ket{\psi'}\).

- As only Fock backends support boson sampling, we are using the Fock simulator backend, indicated with the argument
`'fock'`

, along with a 4 mode register and a Fock state truncation/cutoff dimension of 7 (i.e., all information of quantum amplitudes of Fock states \(\ket{n}\), \(n\geq 7\), is discarded).

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

`all_fock_probs()`

returns a \(\overbrace{D\times D\times\cdots\times D}^{\text{num modes}}\) array, where \(D\) is the Fock basis cutoff truncation, containing the marginal Fock state probabilities - in this example, it would have dimension \(7\times 7\times 7\times 7\). We can then access the probability of measuring a particular output state by indexing. For example,\[\texttt{probs[1, 1, 1, 0]} = \left|\braketD{1,1,1,0}{\psi'}\right|^2\]

After saving the file, you can then run the simulation by executing on the command line:

```
$ python3 boson_sampling.py
```

## Calculating the unitary¶

To explore further, add the following line to the bottom of the file `boson_sampling.py`

and re-run it (alternatively, you can run the full program in an interactive shell/notebook):

```
np.save('boson_fock_output', probs)
```

You should now have the file `boson_fock_output.npy`

in the same directory containing the joint Fock state probabilities of the output modes; to analyze these 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, and load the Fock state probabilities:

```
>>> import numpy as np
>>> from numpy.linalg import multi_dot
>>> from scipy.linalg import block_diag
>>> probs = np.load('boson_fock_output.npy')
```

To start with, let’s generate the matrix representing the unitary transformation of the input mode annihilation and creation operators. The rotation gates simply act as follows,

and thus the column of rotation gates has the following block diagonal form:

Generating this in NumPy:

```
>>> Uphase = np.diag([np.exp(0.5719*1j),np.exp(-1.9782*1j),np.exp(2.0603*1j),np.exp(0.0644*1j)])
```

A single beamsplitter, acting on two input modes \((\hat{a}_1,\hat{a}_2)\), instead acts as follows:

where \(t=\cos(\theta)\) and \(r=e^{i\phi}\sin(\theta)\). Again, like the rotation gate, they combine as block diagonal matrices.

First of all, we need to convert the `BSgate`

arguments, `theta`

and `phi`

(reproduced below for convenience),

```
>>> BSargs = [
... (0.7804, 0.8578),
... (0.06406, 0.5165),
... (0.473, 0.1176),
... (0.563, 0.1517),
... (0.1323, 0.9946),
... (0.311, 0.3231),
... (0.4348, 0.0798),
... (0.4368, 0.6157)
... ]
```

into transmission and reflection amplitudes \(t\) and \(r\):

```
>>> t_r_amplitudes = [(np.cos(q), np.exp(p)*np.sin(q)) for q,p in BSargs]
```

Next, we generate the individual beamsplitter unitary transformations,

```
>>> BSunitaries = [np.array([[t, -np.conj(r)], [r, t]]) for t,r in t_r_amplitudes]
```

before using the SciPy function `scipy.linalg.block_diag()`

to calculate the resulting \(U_{BS_i}\), i.e., the unitary corresponding to each column of ‘independent’ beamsplitters in the above circuit diagram:

```
>>> UBS1 = block_diag(*BSunitaries[0:2])
>>> UBS2 = block_diag([[1]], BSunitaries[2], [[1]])
>>> UBS3 = block_diag(*BSunitaries[3:5])
>>> UBS4 = block_diag([[1]], BSunitaries[5], [[1]])
>>> UBS5 = block_diag(*BSunitaries[6:8])
```

Finally, we combine the unitaries to form a single \(4\times 4\) unitary via matrix multiplication; \(U = U_{BS_5}U_{BS_4}U_{BS_3}U_{BS_2}U_{BS_1}U_{\phi}\). Since `numpy.dot()`

only supports matrix multiplication of two arrays, we instead use `numpy.linalg.multi_dot()`

:

```
>>> U = multi_dot([UBS5, UBS4, UBS3, UBS2, UBS1, Uphase])
```

We find that

## Analysis¶

Now that we have the interferometer unitary transformation \(U\), as well as the ‘experimental’ results, let’s compare the two, and see if the boson sampling result,

holds. First, import the Strawberry Fields generated joint Fock state probability data:

```
>>> probs = np.loadtxt('boson_fock_output.txt')
```

For this example, we’ll consider the output state \(\ket{2,0,0,1}\). Extracting \(\left|\braketT{2,0,0,1}{\psi'}\right|^2\) from the output data, we see that

```
>>> probs[2,0,0,1]
0.10644192724642336
```

Before we can calculate the right hand side of equation, we need a method of calculating the permanent. Since the permanent is classically hard to compute, it is not provided in either NumPy *or* SciPy, so we will use the definition provided in SciPy Issue #7151 on GitHub:

```
>>> def perm(M):
... n = M.shape[0]
... d = np.ones(n)
... j = 0
... s = 1
... f = np.arange(n)
... v = M.sum(axis=0)
... p = np.prod(v)
... while (j < n-1):
... v -= 2*d[j]*M[j]
... d[j] = -d[j]
... s = -s
... prod = np.prod(v)
... p += s*prod
... f[0] = 0
... f[j] = f[j+1]
... f[j+1] = j+1
... j = f[0]
... return p/2**(n-1)
```

Note

This function makes use of the Balasubramanian-Bax-Franklin-Glynn formula, which scales as \(\mathcal{O}(2^{n-1}n^2)\) for an \(n\times n\) matrix.

Finally, how do we determine the submatrix \(U_{st}\)? This isn’t too hard [3]; first, we calculate the submatrix \(U_s\) by taking \(m_k\) copies of the \(k\text{th}\) **columns** of \(U\), where \(m_k\) is the photon number of the \(k\text{th}\) input state.

Since our input state is \(\ket{\psi}=\ket{1,1,0,1}\), we simply take single copies of the first, second, and fourth columns:

```
>>> U[:,[0,1,3]]
```

Next, we take \(n_k\) copies of the \(k\text{th}\) **row** of \(U_s\), where \(n_k\) is the photon number of the \(k\text{th}\) **output** state that is measured. Here, our measurement is \(\ket{2,0,0,1}\bra{2,0,0,1}\) so we take two copies of the first row, and one copy of the last row:

```
>>> U[:,[0,1,3]][[0,0,3]]
```

Now, using the permanent function we defined above, we can take the absolute value squared of the permanent. Finally, we divide by the product of the input and output state number factorials. Since \(0!=1!=1\), we only need to take into account the case \(2!=2\):

```
>>> np.abs(perm(U[:,[0,1,3]][[0,0,3]]))**2 / 2
0.10644192724642332
```

Comparing this to the result from Strawberry Fields, we can see that they only differ at the **17th decimal point**. Calculating the exact percentage difference,

```
>>> BS = np.abs(perm(U[:,[0,1,3]][[0,0,3]]))**2 / 2
>>> SF = probs[2,0,0,1]
>>> 100*np.abs(BS-SF)/BS
3.9113688093093357e-14
```

They agree with almost negligable error! This is due to the high accuracy of the Fock backend, despite the Fock state truncation/cutoff.

This close result stands for any other output Fock state measurement that preserves the photon number, for example:

- \(\ket{3,0,0,0}\bra{3,0,0,0}\):
>>> probs[3,0,0,0] 0.00094584833471324822 >>> np.abs(perm(U[:,[0,1,3]][[0,0,0]]))**2 / 6 0.00094584833471324887

- \(\ket{1,1,0,1}\bra{1,1,0,1}\):
>>> probs[1,1,0,1] 0.17468916048563926 >>> np.abs(perm(U[:,[0,1,3]][[0,1,3]]))**2 / 1 0.17468916048563934

Note

Although returning only an **approximation** of the Fock state joint probability, and thus only approximating the various submatrix permanents, the Fock backends are still computing a classically intractable problem.

This is because approximating the matrix permanent remains a countably hard classical problem.