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