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 [_]. In this tutorial, we will walk through the application of the boson sampling. Background theory ----------------- Introduced by Aaronson and Arkhipov [_], 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 [_]), 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 [_]. 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 [_] .. 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 [_], 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) [_], 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 [_]. 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 Fock(1) | q Vac | q Fock(1) | q # rotation gates Rgate(0.5719) | q Rgate(-1.9782) | q Rgate(2.0603) | q Rgate(0.0644) | q # beamsplitter array BSgate(0.7804, 0.8578) | (q, q) BSgate(0.06406, 0.5165) | (q, q) BSgate(0.473, 0.1176) | (q, q) BSgate(0.563, 0.1517) | (q, q) BSgate(0.1323, 0.9946) | (q, q) BSgate(0.311, 0.3231) | (q, q) BSgate(0.4348, 0.0798) | (q, q) BSgate(0.4368, 0.6157) | (q, q) ###################################################################### # 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([], BSunitaries, []) UBS3 = block_diag(*BSunitaries[3:5]) UBS4 = block_diag([], BSunitaries, []) 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.op.p 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 Fock(1) | q Vac | q Fock(1) | q 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 [_]; # 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 # ---------- # # ..  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. # # ..  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. # # ..  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. # # ..  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. # # ..  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. # # ..  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. # # ..  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.