r"""
.. role:: html(raw)
:format: html
.. _boson_tutorial:
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 :math:`N` single photon Fock states,
:math:`\ket{\psi}=\ket{m_1,m_2,\dots,m_N}`, composed of :math:`b=\sum_i m_i` photons, incident on an
:math:`N`-mode linear interferometer, which performs the following linear transformation of the
input mode creation and annihilations operators:
.. math::
&\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}
Here, the unitary :math:`U` completely describes the interferometer. Thus, the probability of
detecting :math:`n_j` photons at the :math:`j\text{th}` output mode is given by
.. math:: \left|\braketT{n_1,n_2,\dots,n_N}{W}{\psi}\right|^2
where :math:`W` represents the action of :math:`U` on the Fock basis (:math:`W` is simply a
homomorphism of :math:`U`). The remarkable nature of the boson sampling problem to challenge the
extended Church-Turing thesis lies in the fact that [[2]_]
.. math::
\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
:math:`U_{st}`, a submatrix of the interferometer unitary, dependent upon the input and output Fock
states.
.. note::
The permanent of a matrix, defined by
.. math:: \text{Per}(A) = \sum_{\sigma=S_N}\prod_{i=1}^N A_{i\sigma(i)}
where :math:`S_N` is the set of all permutations of :math:`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 :math:`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 (:class:`~strawberryfields.ops.BSgate`) and single-mode phase shifters
(:class:`~strawberryfields.ops.Rgate`)
[[5]_], allowing for a straightforward translation into a CV quantum circuit.
For example, in the case of a 4 mode interferometer, with arbitrary :math:`4\times 4` unitary
:math:`U`, the quantum photonics circuit is given by
.. image:: /tutorials/images/boson_sampling_ex.svg
:align: center
:width: 70%
:target: javascript:void(0);
:html:`
`
In the above, the detectors perform Fock state measurements, and the parameters of the beamsplitters
and the rotation gates determines the unitary :math:`U`. Note that, in order to allow for arbitrary
linear unitaries for :math:`m` imput modes, we must have a minimum of :math:`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 :math:`\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 :class:`~strawberryfields.ops.Fock`
# operator. Mode 2 is initialized as a vacuum state using the ``Vac`` operator (a shortcut to
# :class:`~strawberryfields.ops.Vacuum`). This is *optional* - modes are automatically created in
# the vacuum state upon engine initialization.
#
# ..
#
# 2. Next we apply the rotation gates, :class:`~strawberryfields.ops.Rgate`, to each mode. The
# resulting rotation in the phase space occurs in the anticlockwise direction, with angle
# :math:`\phi`.
#
# ..
#
# 3. Finally, we apply the array of beamsplitters, using the :class:`~strawberryfields.ops.BSgate`
# operator, with arguments ``(theta,phi)``.
#
# * The transmission amplitude is then given by :math:`t=\cos\theta`
# * The reflection amplitude is given by :math:`r=e^{i\phi}\sin{\theta}`
#
# ..
#
# 4. The rotation gate and beamsplitter parameters have been chosen at random, leading to an
# arbitrary unitary :math:`U` acting on the input modes annihilation and creation operators.
# After leaving the beamsplitter array, we will denote the output state by :math:`\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
#
# .. code-block:: python
#
# 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
# :class:`~strawberryfields.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
# :math:`\ket{n}`, :math:`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 :meth:`~strawberryfields.backends.BaseFockState.all_fock_probs` returns a
# :math:`\overbrace{D\times D\times\cdots\times D}^{\text{num modes}}` array, where :math:`D` is the
# Fock basis cutoff truncation, containing the marginal Fock state probabilities - in this example,
# it would have dimension :math:`7\times 7\times 7\times 7`. We can then access the probability of
# measuring a particular output state by indexing. For example,
#
# .. math:: \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])
######################################################################
# 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,
#
# .. math:: R(\phi)\hat{a} = \hat{a}e^{i\phi},
#
# and thus the column of rotation gates has the following block diagonal form:
#
# .. math::
# U_\phi = \left[\begin{matrix}
# e^{i\phi_1} & 0 & \cdots\\
# 0 & e^{i\phi_2} & \cdots \\
# \vdots & \vdots & \ddots \\
# \end{matrix}\right]
#
# 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 :math:`(\hat{a}_1,\hat{a}_2)`, instead acts as
# follows:
#
# .. math::
#
# 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]
#
# where :math:`t=\cos(\theta)` and :math:`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 :math:`t` and :math:`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
# :math:`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 :math:`4\times 4` unitary via matrix
# multiplication; :math:`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))
######################################################################
# We find that
#
# .. math::
# 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]
#
# 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
# :class:`~strawberryfields.ops.GaussianTransform` operation, consisting of a single
# symplectic matrix:
prog_compiled.print()
######################################################################
# 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)
######################################################################
# which agrees with the result above.
#
# The Interferometer operation
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#
# Strawberry Fields supports the :class:`~strawberryfields.ops.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()
######################################################################
# Comparing to the permanent
# --------------------------
#
# Now that we have the interferometer unitary transformation :math:`U`, as well as the 'experimental'
# results, let's compare the two, and see if the boson sampling result,
#
# .. math::
#
# \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 :math:`\ket{2,0,0,1}`. Extracting
# :math:`\left|\braketT{2,0,0,1}{\psi'}\right|^2` from the output data, we see that
print(probs[2,0,0,1])
######################################################################
# 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 :math:`U_{st}`? This isn't too hard [[7]_];
# first, we calculate the submatrix :math:`U_s` by taking :math:`m_k` copies of the :math:`k\text{th}`
# **columns** of :math:`U`, where :math:`m_k` is the photon number of the :math:`k\text{th}` input
# state.
#
# Since our input state is :math:`\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 :math:`n_k` copies of the :math:`k\text{th}` **row** of :math:`U_s`, where
# :math:`n_k` is the photon number of the :math:`k\text{th}` **output** state that is measured.
# Here, our measurement is :math:`\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 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
# :math:`0!=1!=1`, we only need to take into account the case :math:`2!=2`:
print(np.abs(perm(U[:, [0,1,3]][[0,0,3]]))**2 / 2)
######################################################################
# 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)
######################################################################
# 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:
#
# :math:`\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)
######################################################################
# :math:`\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)
######################################################################
# .. 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] 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] 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.