Note

Click here to download the full example code

# Boson sampling and the permanent¶

“The extended Church-Turing thesis states that

‘Any algorithmic process can be simulated efficiently using a probabilistic Turing machine’.This ad-hoc modification of the Church-Turing thesis should leave you feeling rather queasy.” — Nielsen and Chuang [1].

In this tutorial, we will walk through the application of the boson sampling.

## Background theory¶

Introduced by Aaronson and Arkhipov [2], boson sampling presented a slight
deviation from the general approach in quantum computation. Rather than presenting a theoretical
model of universal quantum computation (i.e., a framework that enables quantum simulation of any
arbitrary Hamiltonian [1]), boson sampling-based devices are instead an example
of an **intermediate quantum computer**, designed to experimentally implement a computation that
is thought to be intractable classically [3].

Boson sampling proposes the following quantum linear optics scheme. An array of single-photon sources is set up, designed to simultaneously emit single photon states into a multimode linear interferometer; the results are then generated by sampling from the probability of single photon measurements from the output of the linear interferometer.

For example, consider \(N\) single photon Fock states, \(\ket{\psi}=\ket{m_1,m_2,\dots,m_N}\), composed of \(b=\sum_i m_i\) photons, incident on an \(N\)-mode linear interferometer, which performs the following linear transformation of the input mode creation and annihilations operators:

Here, the unitary \(U\) completely describes the interferometer. Thus, the probability of detecting \(n_j\) photons at the \(j\text{th}\) output mode is given by

where \(W\) represents the action of \(U\) on the Fock basis (\(W\) is simply a homomorphism of \(U\)). The remarkable nature of the boson sampling problem to challenge the extended Church-Turing thesis lies in the fact that [2]

i.e., the sampled single photon probability distribution is proportional to the **permanent** of
\(U_{st}\), a submatrix of the interferometer unitary, dependent upon the input and output Fock
states.

Note

The permanent of a matrix, defined by

where \(S_N\) is the set of all permutations of \(N\) elements, is calculated in a similar fashion to the determinant, but unlike the determinant, the signatures of the permutations are not taken into account - every permutation is taken as a positive quantity.

In graph theory, the permanent calculates the number of perfect matchings in a bipartite graph with adjacency matrix \(A\).

Whilst the determinant can be calculated efficiently on classical computers, computing the permanent
belongs to the computational complexity class of #P-Hard problems [4], which are
strongly believed to be classically hard to calculate. (Surprisingly, even calculating the permanent
in an *approximate* manner is a member of #P and intractable classically).

This implies that simulating boson sampling cannot be done efficiently on a classical computer, providing a potential challenge to the extended Church-Turing thesis, and demonstrating the power of (non-universal) quantum computation.

## Circuit construction and simulation¶

In quantum linear optics, the multimode linear interferometer is commonly decomposed into two-mode
beamsplitters (`BSgate`

) and single-mode phase shifters
(`Rgate`

)
[5], 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 quantum photonics circuit is given by

In the above, the detectors perform Fock state measurements, and the parameters of the beamsplitters and the rotation gates determines the unitary \(U\). Note that, in order to allow for arbitrary linear unitaries for \(m\) imput modes, we must have a minimum of \(m+1\) columns in the beamsplitter array [6].

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 *
# initialize a 4 mode program
boson_sampling = sf.Program(4)
with boson_sampling.context as q:
# 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) | q[0]
Rgate(-1.9782) | q[1]
Rgate(2.0603) | q[2]
Rgate(0.0644) | q[3]
# 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])
```

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'}\).

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.If we wish to simulate Fock measurements, we can additionally include

MeasureFock() | q

after the beamsplitter array. After constructing the circuit and running the engine, the values of the Fock state measurements will be available within the

`samples`

attribute of the`Result`

object returned by the engine.

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.

As only Fock backends support the non-Gaussian initial states used in boson sampling,
we will use the NumPy `'fock'`

backend, along with 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 can now execute the program with the engine:

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,

```
# extract the joint Fock probabilities
probs = results.state.all_fock_probs()
# print the joint Fock state probabilities
print(probs[1, 1, 0, 1])
print(probs[2, 0, 0, 1])
```

Out:

```
0.17468916048563926
0.1064419272464234
```

## Calculating the unitary¶

Let’s explore these results further. To do so, however, we’ll need to compute the unitary that is being applied by our beamsplitters and rotation gates.

### By hand¶

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

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*1j)*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])
print(np.round(U, 4))
```

Out:

```
[[ 0.2195-0.2565j 0.6111+0.5242j -0.1027+0.4745j -0.0273+0.0373j]
[ 0.4513+0.6026j 0.457 +0.0123j 0.1316-0.4504j 0.0353-0.0532j]
[ 0.0387+0.4927j -0.0192-0.3218j -0.2408+0.5244j -0.4584+0.3296j]
[-0.1566+0.2246j 0.11 -0.1638j -0.4212+0.1836j 0.8188+0.068j ]]
```

We find that

### The Gaussian compiler¶

While we have done it by hand, Strawberry Fields also supports a Gaussian unitary compiler, that allows us to compile our program into a single Gaussian unitary.

Note that we must create a new program without the initial Fock states,
as the Gaussian unitary compiler only works with *Gaussian operations*.
Compiling the program:

```
prog_unitary = sf.Program(4)
prog_unitary.circuit = boson_sampling.circuit[4:]
prog_compiled = prog_unitary.compile(compiler="gaussian_unitary")
```

Printing `prog_compiled`

, we see it now consists of a single
`GaussianTransform`

operation, consisting of a single
symplectic matrix:

Out:

```
GaussianTransform([[ 0.2195 0.6111 -0.1027 -0.0273 0.2565 -0.5242 -0.4745 -0.0373]
[ 0.4513 0.457 0.1316 0.0353 -0.6026 -0.0123 0.4504 0.0532]
[ 0.0387 -0.0192 -0.2408 -0.4584 -0.4927 0.3218 -0.5244 -0.3296]
[-0.1566 0.11 -0.4212 0.8188 -0.2246 0.1638 -0.1836 -0.068 ]
[-0.2565 0.5242 0.4745 0.0373 0.2195 0.6111 -0.1027 -0.0273]
[ 0.6026 0.0123 -0.4504 -0.0532 0.4513 0.457 0.1316 0.0353]
[ 0.4927 -0.3218 0.5244 0.3296 0.0387 -0.0192 -0.2408 -0.4584]
[ 0.2246 -0.1638 0.1836 0.068 -0.1566 0.11 -0.4212 0.8188]]) | (q[0], q[1], q[2], q[3])
```

We can easily extract this symplectic matrix, and rewrite it as a unitary matrix:

```
S = prog_compiled.circuit[0].op.p[0]
U = S[:4, :4] + 1j*S[4:, :4]
print(U)
```

Out:

```
[[ 0.2195-0.2565j 0.6111+0.5242j -0.1027+0.4745j -0.0273+0.0373j]
[ 0.4513+0.6026j 0.457 +0.0123j 0.1316-0.4504j 0.0353-0.0532j]
[ 0.0387+0.4927j -0.0192-0.3218j -0.2408+0.5244j -0.4584+0.3296j]
[-0.1566+0.2246j 0.11 -0.1638j -0.4212+0.1836j 0.8188+0.068j ]]
```

which agrees with the result above.

### The Interferometer operation¶

Strawberry Fields supports the `Interferometer`

operation, which
allows numeric unitary matrices to be directly embedded in programs and decomposed into
beamsplitters and rotation gates:

```
boson_sampling = sf.Program(4)
with boson_sampling.context as q:
# prepare the input fock states
Fock(1) | q[0]
Fock(1) | q[1]
Vac | q[2]
Fock(1) | q[3]
Interferometer(U) | q
```

Compiling this for the Fock backend, and printing the result:

```
boson_sampling.compile(compiler="fock").print()
```

Out:

```
Fock(1) | (q[0])
Fock(1) | (q[1])
Vac | (q[2])
Fock(1) | (q[3])
Rgate(-3.124) | (q[0])
BSgate(0.9465, 0) | (q[0], q[1])
Rgate(2.724) | (q[2])
BSgate(0.09485, 0) | (q[2], q[3])
Rgate(-0.9705) | (q[1])
BSgate(0.7263, 0) | (q[1], q[2])
Rgate(-1.788) | (q[0])
BSgate(0.8246, 0) | (q[0], q[1])
Rgate(5.343) | (q[0])
Rgate(2.93) | (q[1])
Rgate(3.133) | (q[2])
Rgate(0.07904) | (q[3])
BSgate(-0.533, 0) | (q[2], q[3])
Rgate(2.45) | (q[2])
BSgate(-0.03962, 0) | (q[1], q[2])
Rgate(2.508) | (q[1])
```

## Comparing to the permanent¶

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.

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

```
print(probs[2,0,0,1])
```

Out:

```
0.1064419272464234
```

Before we can calculate the right-hand-side of the 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 Walrus library, installed
alongside Strawberry Fields:

```
from thewalrus import perm
```

Finally, how do we determine the submatrix \(U_{st}\)? This isn’t too hard [7];
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]]
```

Out:

```
array([[ 0.2195-0.2565j, 0.6111+0.5242j, -0.0273+0.0373j],
[ 0.4513+0.6026j, 0.457 +0.0123j, 0.0353-0.0532j],
[ 0.0387+0.4927j, -0.0192-0.3218j, -0.4584+0.3296j],
[-0.1566+0.2246j, 0.11 -0.1638j, 0.8188+0.068j ]])
```

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]]
```

Out:

```
array([[ 0.2195-0.2565j, 0.6111+0.5242j, -0.0273+0.0373j],
[ 0.2195-0.2565j, 0.6111+0.5242j, -0.0273+0.0373j],
[-0.1566+0.2246j, 0.11 -0.1638j, 0.8188+0.068j ]])
```

Now, using the permanent function we imported 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\):

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

Out:

```
0.10644192724642326
```

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]
print(100*np.abs(BS-SF)/BS)
```

Out:

```
1.3037896031031125e-13
```

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}\):

```
print(probs[3,0,0,0])
print(np.abs(perm(U[:, [0,1,3]][[0,0,0]]))**2 / 6)
```

Out:

```
0.000945848334713249
0.0009458483347132489
```

\(\ket{1,1,0,1}\bra{1,1,0,1}\):

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

Out:

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

## References¶

- 1(1,2)
M.A. Nielsen and I.L. Chuang. Quantum Computation and Quantum Information. Cambridge University Press, 2010. ISBN 9780511992773. URL: https://books.google.ca/books?id=JRz3jgEACAAJ.

- 2(1,2)
Scott Aaronson and Alex Arkhipov. The computational complexity of linear optics. Theory of Computing, 9(1):143–252, 2013. doi:10.4086/toc.2013.v009a004.

- 3
Max Tillmann, Borivoje Dakić, René Heilmann, Stefan Nolte, Alexander Szameit, and Philip Walther. Experimental boson sampling. Nature Photonics, 7(7):540–544, May 2013. doi:10.1038/nphoton.2013.102.

- 4
L.G. Valiant. The complexity of computing the permanent. Theoretical Computer Science, 8(2):189–201, 1979. doi:10.1016/0304-3975(79)90044-6.

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

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

- 7
Max Tillmann, Borivoje Dakić, René Heilmann, Stefan Nolte, Alexander Szameit, and Philip Walther. Experimental boson sampling. Nature Photonics, 7(7):540–544, May 2013. doi:10.1038/nphoton.2013.102.

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

## Contents

## Downloads

## Related tutorials