\n```\nSimulating this circuit using PennyLane is easy; we can simply read off\nthe gates from left to right, and convert it into a QNode.\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "import numpy as np\n\n# set the random seed\nnp.random.seed(42)\n\n# import PennyLane\nimport pennylane as qml" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We must define the unitary matrix we would like to embed in the circuit.\nWe will use SciPy to generate a Haar-random unitary:\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "from scipy.stats import unitary_group\n\n# define the linear interferometer\nU = unitary_group.rvs(4)\nprint(U)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can now use this to construct the circuit, choosing a compatible\ndevice. For the simulation, we can use the Strawberry Fields Gaussian\nbackend. This backend is perfectly suited for simulation of GBS, as the\ninitial states are Gaussian, and all gates transform Gaussian states to\nother Gaussian states.\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "n_wires = 4\ncutoff = 10\n\ndev = qml.device(\"strawberryfields.gaussian\", wires=n_wires, cutoff_dim=cutoff)\n\n\n@qml.qnode(dev)\ndef gbs_circuit():\n # prepare the input squeezed states\n for i in range(n_wires):\n qml.Squeezing(1.0, 0.0, wires=i)\n\n # linear interferometer\n qml.InterferometerUnitary(U, wires=range(n_wires))\n return qml.probs(wires=range(n_wires))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "A couple of things to note in this particular example:\n\n1. To prepare the input single mode squeezed vacuum state\n $\\ket{re^{i\\phi}}$, where $r = 1$ and $\\phi=0$, we apply a squeezing\n gate (`~pennylane.Squeezing`{.interpreted-text role=\"class\"}) to\n each of the wires (initially in the vacuum state).\n2. Next we apply the linear interferometer to all four wires using\n `~pennylane.Interferometer`{.interpreted-text role=\"class\"} and the\n unitary matrix `U`. This operator decomposes the unitary matrix\n representing the linear interferometer into single-mode rotation\n gates (`~pennylane.PhaseShift`{.interpreted-text role=\"class\"}) and\n two-mode beamsplitters (`~pennylane.Beamsplitter`{.interpreted-text\n role=\"class\"}). After applying the interferometer, we will denote\n the output state by $\\ket{\\psi'}$.\n3. GBS takes place physically in an infinite-dimensional Hilbert space,\n which is not practical for simulation. We need to set an upper limit\n on the maximum number of photons we can detect. This is the `cutoff`\n value we defined above; we will only be considering detection events\n containing 0 to 9 photons per mode.\n\nWe can now execute the QNode, and extract the resulting probability\ndistribution:\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "probs = gbs_circuit().reshape([cutoff] * n_wires)\nprint(probs.shape)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For example, element `[1,2,0,1]` represents the probability of detecting\n1 photon on wire `0` and wire `3`, and 2 photons at wire `1`, i.e., the\nvalue\n\n$$\\text{prob}(1,2,0,1) = \\left|\\braketD{1,2,0,1}{\\psi'}\\right|^2.$$\n\nLet\\'s extract and view the probabilities of measuring various Fock\nstates.\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Fock states to measure at output\nmeasure_states = [(0, 0, 0, 0), (1, 1, 0, 0), (0, 1, 0, 1), (1, 1, 1, 1), (2, 0, 0, 0)]\n\n# extract the probabilities of calculating several\n# different Fock states at the output, and print them out\nfor i in measure_states:\n print(f\"|{''.join(str(j) for j in i)}>: {probs[i]}\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The GBS Distribution\n====================\n\nHamilton et al. showed that the probability of measuring a final state\ncontaining only 0 or 1 photons per mode is given by\n\n$$\\left|\\left\\langle{n_1,n_2,\\dots,n_N}\\middle|{\\psi'}\\right\\rangle\\right|^2 =\n\\frac{\\left|\\text{Haf}[(U(\\bigoplus_i\\tanh(r_i))U^T)]_{st}\\right|^2}{\\prod_{i=1}^N \\cosh(r_i)}$$\n\ni.e., the sampled single-photon probability distribution is proportional\nto the **hafnian** of a submatrix of $U(\\bigoplus_i\\tanh(r_i))U^T$.\n\n::: {.note}\n::: {.title}\nNote\n:::\n\nThe hafnian of a matrix is defined by\n\n$$\\text{Haf}(A) = \\sum_{\\sigma \\in \\text{PMP}_{2N}}\\prod_{i=1}^N A_{\\sigma(2i-1)\\sigma(2i)},$$\n\nwhere $\\text{PMP}_{2N}$ is the set of all perfect matching permutations\nof $2N$ elements. In graph theory, the hafnian calculates the number of\nperfect\n[matchings](https://en.wikipedia.org/wiki/Matching_(graph_theory)) in a\ngraph with adjacency matrix $A$.\n\nCompare this to the permanent, which calculates the number of perfect\nmatchings on a *bipartite* graph. Notably, the permanent appears in\nvanilla Boson Sampling in a similar way that the hafnian appears in GBS.\nThe hafnian turns out to be a generalization of the permanent, with the\nrelationship\n\n$$\\begin{aligned}\n\\text{Per(A)} = \\text{Haf}\\left(\\left[\\begin{matrix}\n 0&A\\\\ A^T&0\n\\end{matrix}\\right]\\right).\n\\end{aligned}$$\n:::\n\nAs any algorithm that could calculate (or even approximate) the hafnian\ncould also calculate the permanent\\-\\--a [\\#P-hard\nproblem](https://en.wikipedia.org/wiki/%E2%99%AFP)\\-\\--it follows that\ncalculating or approximating the hafnian must also be a classically hard\nproblem. This lies behind the classical hardness of GBS.\n\nIn this demo, we will use the same squeezing parameter, $z=r$, for all\ninput states; this allows us to simplify this equation. To start with,\nthe hafnian expression simply becomes $\\text{Haf}[(UU^T\\tanh(r))]_{st}$,\nremoving the need for the direct sum.\n\nThus, we have\n\n$$\\left|\\left\\langle{n_1,n_2,\\dots,n_N}\\middle|{\\psi'}\\right\\rangle\\right|^2 =\n\\frac{\\left|\\text{Haf}[(UU^T\\tanh(r))]_{st}\\right|^2}{n_1!n_2!\\cdots n_N!\\cosh^N(r)}.$$\n\nNow that we have the theoretical formulas, as well as the probabilities\nfrom our simulated GBS QNode, we can compare the two and see whether\nthey agree.\n\nIn order to calculate the probability of different GBS events\nclassically, we need a method for calculating the hafnian. For this, we\nwill use [The Walrus](https://the-walrus.readthedocs.io) library (which\nis installed as a dependency of the PennyLane-SF plugin):\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "from thewalrus import hafnian as haf" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now, for the right-hand side numerator, we first calculate the submatrix\n$A = [(UU^T\\tanh(r))]_{st}$:\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "A = np.dot(U, U.T) * np.tanh(1)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In GBS, we determine the submatrix by taking the rows and columns\ncorresponding to the measured Fock state. For example, to calculate the\nsubmatrix in the case of the output measurement\n$\\left|{1,1,0,0}\\right\\rangle$, we have\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "print(A[:, [0, 1]][[0, 1]])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "i.e., we consider only the rows and columns where a photon was detected,\nwhich gives us the submatrix corresponding to indices $0$ and $1$.\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Comparing to simulation\n=======================\n\nNow that we have a method for calculating the hafnian, let\\'s compare\nthe output to that provided by the PennyLane QNode.\n\n**Measuring** $\\ket{0,0,0,0}$ **at the output**\n\nThis corresponds to the hafnian of an *empty* matrix, which is simply 1:\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "print(1 / np.cosh(1) ** 4)\nprint(probs[0, 0, 0, 0])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Measuring** $\\ket{1,1,0,0}$ **at the output**\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "A = (np.dot(U, U.T) * np.tanh(1))[:, [0, 1]][[0, 1]]\nprint(np.abs(haf(A)) ** 2 / np.cosh(1) ** 4)\nprint(probs[1, 1, 0, 0])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Measuring** $\\ket{0,1,0,1}$ **at the output**\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "A = (np.dot(U, U.T) * np.tanh(1))[:, [1, 3]][[1, 3]]\nprint(np.abs(haf(A)) ** 2 / np.cosh(1) ** 4)\nprint(probs[0, 1, 0, 1])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Measuring** $\\ket{1,1,1,1}$ **at the output**\n\nThis corresponds to the hafnian of the full matrix $A=UU^T\\tanh(r)$:\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "A = np.dot(U, U.T) * np.tanh(1)\nprint(np.abs(haf(A)) ** 2 / np.cosh(1) ** 4)\nprint(probs[1, 1, 1, 1])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Measuring** $\\ket{2,0,0,0}$ **at the output**\n\nSince we have two photons in mode `q[0]`, we take two copies of the\nfirst row and first column, making sure to divide by $2!$:\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "A = (np.dot(U, U.T) * np.tanh(1))[:, [0, 0]][[0, 0]]\nprint(np.abs(haf(A)) ** 2 / (2 * np.cosh(1) ** 4))\nprint(probs[2, 0, 0, 0])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The PennyLane simulation results agree (with almost negligible numerical\nerror) to the expected result from the Gaussian boson sampling equation!\n\nThis demo provides an entry-level walkthrough to the ideas behind GBS,\nproviding you with the basic code needed for exploring the ideas behind\nthe photonic quantum advantage paper. Try changing the number of modes,\nthe number of injected squeezed states, or the cutoff dimension, and see\nhow each of these affect the classical computation time. If you\\'re\ninterested in learning more about GBS, or about photonic quantum\ncomputing in general, the [Strawberry Fields\nwebsite](https://strawberryfields.ai/) is a great resource.\n\nReferences\n==========\n\nAbout the authors\n=================\n" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.9.15" } }, "nbformat": 4, "nbformat_minor": 0 }