r"""
.. role:: html(raw)
:format: html
.. _gaussian_boson_tutorial:
Gaussian boson sampling and the Hafnian
=======================================
"If you need to wait exponential time for [your single photon sources to emit simultaneously],
then there would seem to be no advantage over classical computation. This is the reason why so
far, boson sampling has only been demonstrated with 3-4 photons. When faced with these problems,
until recently, all we could do was shrug our shoulders." - `Scott Aaronson
`_
In this tutorial, we will walk through the application of the Gaussian boson sampling.
It is also strongly recommended to read the :doc:`boson sampling tutorial
`, which introduces a few concepts also used in Gaussian boson sampling.
Background theory
-----------------
While boson sampling allows the experimental implementation of a sampling problem that is countably
hard classically, one of the main issues it has in experimental setups is one of **scalability**,
due to its dependence on an array of simultaneously emitting single photon sources. Currently, most
physical implementations of boson sampling make use of a process known as `spontaneous parametric
down-conversion `_ to generate
the single-photon source inputs. However, this method is non-deterministic --- as the number of
modes in the apparatus increases, the average time required until every photon source emits a
simultaneous photon increases exponentially.
In order to simulate a deterministic single-photon source array, several variations on boson
sampling have been proposed; the most well-known being scattershot boson sampling [[1]_].
However, a recent boson sampling variation by Hamilton et al. [[2]_] mitigates the
need for single photon Fock states altogether, by showing that incident Gaussian states
--- in this case, single mode squeezed states --- can produce problems in the same computational
complexity class as boson sampling. Even more significantly, this negates the scalability
problem with single photon sources, as single mode squeezed states can be deterministically
generated simultaneously.
The Gaussian boson sampling scheme remains, on initial observation, quite similar to that of boson
sampling:
* :math:`N` single mode squeezed states :math:`\ket{z}`, with squeezing parameter
:math:`z=re^{i\phi}`, enter an :math:`N` mode linear interferometer described by unitary :math:`U`
simultaneously.
* Each output mode of the interferometer (denoted by state :math:`\ket{\psi'}`) is then measured in
the Fock basis, :math:`\bigotimes_i n_i\ket{n_i}\bra{n_i}`.
Without loss of generality, we can absorb the squeezing phase parameter :math:`\phi` into the
interferometer, and set :math:`\phi=0` for convenience.
Using phase space methods, Hamilton et al. [[2]_] showed that the probability of
measuring a Fock state containing only 0 or 1 photons per mode is given by
.. math::
\left|\left\langle{n_1,n_2,\dots,n_N}\middle|{\psi'}\right\rangle\right|^2 =
\frac{\left|\text{Haf}[(U(\bigoplus_i\tanh(r_i))U^T)]_{st}\right|^2}{\prod_{i=1}^N \cosh(r_i)}
i.e., the sampled single photon probability distribution is proportional to the **hafnian** of a
submatrix of :math:`U(\bigoplus_i\tanh(r_i))U^T`, dependent upon the output covariance matrix.
.. note::
The hafnian of a matrix is defined by
.. math:: \text{Haf}(A) = \frac{1}{n!2^n}\sum_{\sigma \in S_{2N}}\prod_{i=1}^N A_{\sigma(2i-1)\sigma(2i)}
where :math:`S_{2N}` is the set of all permutations of :math:`2N` elements. In graph theory, the
hafnian calculates the number of perfect `matchings
`_ in an **arbitrary graph** with
adjacency matrix :math:`A`.
Compare this to the permanent, which calculates the number of perfect matchings on a *bipartite*
graph - the hafnian turns out to be a generalization of the permanent, with the relationship
.. math::
\text{Per(A)} = \text{Haf}\left(\left[\begin{matrix}
0&A\\ A^T&0
\end{matrix}\right]\right)
As any algorithm that could calculate (or even approximate) the hafnian could also calculate the
permanent - a #P-hard problem - it follows that calculating or approximating the hafnian must also
be a classically hard problem.
Circuit construction and simulation
-----------------------------------
As with the boson sampling problem, the multimode linear interferometer can be decomposed into
two-mode beamsplitters (:class:`~strawberryfields.ops.BSgate`) and single-mode phase shifters
(:class:`~strawberryfields.ops.Rgate`)
[[3]_], 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 CV quantum circuit for Gaussian boson sampling is given by
.. image:: /tutorials/images/gaussian_boson_sampling.svg
:align: center
:width: 70%
:target: javascript:void(0);
.. raw:: html
In the above, the single mode squeeze states all apply identical squeezing :math:`z=r`, the
parameters of the beamsplitters and the rotation gates determine the unitary :math:`U`, and finally
the detectors perform Fock state measurements on the output modes. As with boson sampling, for
:math:`N` input modes, we must have a minimum of :math:`N+1` columns in the beamsplitter array
[[4]_].
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 *
######################################################################
# Unlike the :doc:`run_boson_sampling`, we will directly apply a Gaussian unitary
# to the circuit using the :class:`~strawberryfields.ops.Interferometer` operation.
# First, we must define the unitary matrix we would like to embed in the circuit:
# 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]
])
######################################################################
# We can use this to now construct the circuit:
# create the 4 mode Strawberry Fields program
gbs = sf.Program(4)
with gbs.context as q:
# 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
######################################################################
# A couple of things to note in this particular example:
#
# ..
#
# 1. To prepare the input single mode squeezed vacuum state :math:`\ket{z}` where :math:`z = 1`, we
# apply a squeezing gate :class:`~strawberryfields.ops.Sgate` to each of the modes (initially in
# the vacuum state).
#
# ..
#
# 2. Next we apply the linear interferometer to all four modes, using the decomposition operator
# :class:`~strawberryfields.ops.Interferometer`, and the unitary matrix ``U``. This operator
# decomposes the unitary matrix representing the linear interferometer into single mode
# rotation gates :class:`~strawberryfields.ops.Rgate`, and two-mode beamsplitters
# :class:`~strawberryfields.ops.BSgate`. After applying the interferometer, we will denote the
# output state by :math:`\ket{\psi'}`.
#
# .. note::
#
# You can view the decomposed beamsplitters and rotation gates which correspond to the linear
# interferometer ``U`` by calling :meth:`eng.print_applied()
# ` after running the engine.
#
# .. note::
#
# The interferometer applied here is **identical** to that from the
# :doc:`run_boson_sampling`. As a result, the decomposed beamsplitter and rotation gate
# parameters will also be identical.
#
# ..
#
# 3. 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.
#
# Unlike the :doc:`boson sampling tutorial `, the lack of Fock states means we
# can now use the Numpy-based ``'gaussian'`` backend, 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.
#
# 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.
eng = sf.Engine(backend="gaussian")
results = eng.run(gbs)
######################################################################
# The state method :meth:`~strawberryfields.backends.BaseState.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
#
# .. math:: \text{prob}(1,2,0,1) = \left|\braketD{1,2,0,1}{\psi'}\right|^2
#
# The Fock state method :meth:`~strawberryfields.backends.BaseFockState.all_fock_probs`, used
# previously to return *all* Fock state probabilities as an array, is **not supported** by the
# Gaussian backend. 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.
#
# Let's use :meth:`~strawberryfields.backends.BaseState.fock_prob` in a for-loop to extract
# the probabilities of measuring various Fock states:
# 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 = results.state.fock_prob(i)
print("|{}>: {}".format("".join(str(j) for j in i), prob))
######################################################################
# Equally squeezed inputs
# -----------------------
#
# As shown earlier, the formula for calculating the output
# Fock state probability in Gaussian boson sampling is given by
#
# .. math::
#
# \left|\left\langle{n_1,n_2,\dots,n_N}\middle|{\psi'}\right\rangle\right|^2 =
# \frac{\left|\text{Haf}[(U\bigoplus_i\tanh(r_i)U^T)]_{st}\right|^2}{n_1!n_2!\cdots n_N!
# \cosh(r_i)}
#
# where :math:`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, :math:`z=r`, for
# all input states - this allows us to simplify this equation. To start with, the hafnian expression
# simply becomes :math:`\text{Haf}[(UU^T\tanh(r))]_{st}`, removing the need for the tensor sum.
#
# Thus, we have
#
# .. math::
#
# \left|\left\langle{n_1,n_2,\dots,n_N}\middle|{\psi'}\right\rangle\right|^2 =
# \frac{\left|\text{Haf}[(UU^T\tanh(r))]_{st}\right|^2}{n_1!n_2!\cdots n_N!\cosh^N(r)}.
#
# 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 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 Walrus
# `_ library, installed alongside Strawberry Fields:
from thewalrus import hafnian as haf
######################################################################
# Now, for the right hand side numerator, we first calculate the submatrix
# :math:`[(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 :math:`\left|{1,1,0,0}\right\rangle`,
print(B[:, [0, 1]][[0, 1]])
######################################################################
# 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** :math:`\ket{0,0,0,0}` **at the output**
#
# This corresponds to the hafnian of an *empty* matrix, which is simply 1:
print(1 / np.cosh(1) ** 4)
print(results.state.fock_prob([0, 0, 0, 0]))
######################################################################
# **Measuring** :math:`\ket{1,1,0,0}` **at the output**
B = (np.dot(U, U.T) * np.tanh(1))[:, [0, 1]][[0, 1]]
print(np.abs(haf(B)) ** 2 / np.cosh(1) ** 4)
print(results.state.fock_prob([1, 1, 0, 0]))
######################################################################
# **Measuring** :math:`\ket{0,1,0,1}` **at the output**
B = (np.dot(U, U.T) * np.tanh(1))[:, [1, 3]][[1, 3]]
print(np.abs(haf(B)) ** 2 / np.cosh(1) ** 4)
print(results.state.fock_prob([0, 1, 0, 1]))
######################################################################
# **Measuring** :math:`\ket{1,1,1,1}` **at the output**
#
# This corresponds to the hafnian of the full matrix :math:`B=UU^T\tanh(r)`:
B = (np.dot(U, U.T) * np.tanh(1))
print(np.abs(haf(B)) ** 2 / np.cosh(1) ** 4)
print(results.state.fock_prob([1, 1, 1, 1]))
######################################################################
# **Measuring** :math:`\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 :math:`2!`:
B = (np.dot(U, U.T) * np.tanh(1))[:, [0, 0]][[0, 0]]
print(np.abs(haf(B)) ** 2 / (2 * np.cosh(1) ** 4))
print(results.state.fock_prob([2, 0, 0, 0]))
######################################################################
# The Strawberry Field simulation results agree (with almost negligible numerical error) to the
# expected result from the Gaussian boson sampling equation!
#
# .. warning::
#
# Repeat this tutorial with
#
# 1. A Fock backend such as ``'fock'`` instead of the Gaussian backend.
#
# 2. Different beamsplitter and rotation parameters.
#
# 3. Input states with *differing* squeezed values
# :math:`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.
#
# References
# ----------
#
# .. [1] A. P. Lund, A. Laing, S. Rahimi-Keshari, T. Rudolph, J. L. O’Brien, and T. C. Ralph.
# Boson sampling from a gaussian state. Physical Review Letters, 113:100502, Sep 2014.
# doi:10.1103/PhysRevLett.113.100502.
#
# .. [2] Craig S. Hamilton, Regina Kruse, Linda Sansoni, Sonja Barkhofen, Christine Silberhorn,
# and Igor Jex. Gaussian boson sampling. Physical Review Letters, 119:170501, Oct 2017.
# arXiv:1612.01199, doi:10.1103/PhysRevLett.119.170501.
#
# .. [3] 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.
#
# .. [4] 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.