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:

\[\begin{split}&\hat{a}^\dagger_{out_k} = \sum_{j=0}^N U_{kj}\hat{a}^\dagger_{in_j}\\ &\hat{a}_{out_k} = \sum_{j=0}^N U_{jk}^\dagger\hat{a}_{in_j}\end{split}\]

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

\[\left|\braketT{n_1,n_2,\dots,n_N}{W}{\psi}\right|^2\]

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]

\[\left|\braketT{n_1,n_2,\dots,n_N}{W}{\psi}\right|^2 = \frac{\left|\text{Per}(U_{st})\right|^2}{m_1!m_2!\cdots m_N!n_1!n_2!\cdots n_N!}\]

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

\[\text{Per}(A) = \sum_{\sigma=S_N}\prod_{i=1}^N A_{i\sigma(i)}\]

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:

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

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

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

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

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

eng = sf.Engine(backend="fock", backend_options={"cutoff_dim": 7})

We can now execute the program with the engine:

results = eng.run(boson_sampling)

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

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

\[R(\phi)\hat{a} = \hat{a}e^{i\phi},\]

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

\[\begin{split}U_\phi = \left[\begin{matrix} e^{i\phi_1} & 0 & \cdots\\ 0 & e^{i\phi_2} & \cdots \\ \vdots & \vdots & \ddots \\ \end{matrix}\right]\end{split}\]

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:

\[\begin{split}BS(\theta, \phi) (\hat{a}_1,\hat{a}_2) = \left[\begin{matrix} t & -r^*\\ r & t\\ \end{matrix}\right] \left[\begin{matrix} \hat{a}_1\\ \hat{a}_2 \end{matrix}\right]\end{split}\]

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

\[\begin{split}U = \left[\begin{matrix} 0.2195-0.2565i & 0.6111+0.5242i & -0.1027+0.4745i & -0.0273+0.0373i\\ 0.4513+0.6026i & 0.4570+0.0123i & 0.1316-0.4504i & 0.0353-0.0532i\\ 0.0387+0.4927i & -0.0192-0.3218i & -0.2408+0.5244i & -0.4584+0.3296i\\ -0.1566+0.2246i & 0.1100-0.1638i & -0.4212+0.1836i & 0.8188+0.068i \end{matrix}\right]\end{split}\]

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:

prog_compiled.print()

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,

\[\left|\braketD{n_1,n_2,\dots,n_N}{\psi'}\right|^2 = \frac{\left|\text{Perm}(U_{st})\right|^2}{m_1!m_2!\cdots m_N!n_1!n_2!\cdots n_N!},\]

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

print(100*np.abs(BS-SF)/BS)

Out:

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

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

Out:

0.0009458483347132492
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.17468916048563932
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.

demos/run_boson_sampling

 

Download Python script

 

Download Notebook