Gaussian boson sampling tutorial


This tutorial is also available in the form of an interactive Jupyter Notebook GaussianBosonSampling.ipynb

In this short tutorial, we will walk through the application of the Gaussian boson sampling Blackbird example provided in the Gaussian boson sampling section of the Quantum algorithms page. Be sure to read through that section to familiarise yourself with the concepts behind Gaussian boson sampling, as well as the introductory teleportation tutorial, before attempting this tutorial. It is also strongly recommended to read the boson sampling tutorial, which introduces a few concepts also used in Gaussian boson sampling.

Circuit construction and simulation

A 4-mode Gaussian boson sampling circuit is given by


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, create a new text file with the name, and open it with an editor. Copy and paste the following Strawberry Fields code:

#!/usr/bin/env python3
import numpy as np
import strawberryfields as sf
from strawberryfields.ops import *

# initialize engine and program objects
eng = sf.Engine(backend="gaussian")
gbs = sf.Program(4)

# define the linear interferometer
U = np.array([
 [ 0.219546940711-0.256534554457j, 0.611076853957+0.524178937791j,
 [ 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.156619083736+0.224568570065j, 0.109992223305-0.163750223027j,
    -0.421179844245+0.183644837982j, 0.818769184612+0.068015658737j]

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
    # end circuit

# run the engine
results =

# 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))


This example program is included with Strawberry Fields as examples/

A couple of things to note in this particular example:

  1. To prepare the input single mode squeezed vacuum state \(\ket{z}\) where \(z = 1\), we apply a squeezing gate Sgate to each of the modes (initially in the vacuum state).
  1. Next we apply the linear interferometer to all four modes, using the decomposition operator Interferometer, and the unitary matrix U. This operator decomposes the unitary matrix representing the linear interferometer into single mode rotation gates Rgate, and two-mode beamsplitters BSgate. After applying the interferometer, we will denote the output state by \(\ket{\psi'}\).


    You can view the decomposed beamsplitters and rotation gates which correspond to the linear interferometer U by calling eng.print_applied() after running the engine.


    The interferometer applied here is identical to that from the boson sampling tutorial. As a result, the decomposed beamsplitter and rotation gate parameters will also be identical.

  1. Unlike the 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.
  1. 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.

    The state method 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

    \[\text{prob}(1,2,0,1) = \left|\braketD{1,2,0,1}{\psi'}\right|^2\]

    The Fock state method 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.

The simulation can be run by executing the file from the command line, resulting in the following output:

$ python3
|0000>: 0.176378447614135
|1100>: 0.0685595637122246
|0101>: 0.002056097258977398
|1111>: 0.00834294639986785
|2000>: 0.01031294525345511

Equally squeezed inputs

Recall that, as shown in Ref. [13] and reviewed in the Gaussian boson sampling section of the quantum algorithms page, the formula for calculating the output Fock state probability,

\[\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 \(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, \(z=r\), for all input states - this allows us to simplify this equation. To start with, the hafnian expression simply becomes \(\text{Haf}[(UU^T\tanh(r))]_{st}\), removing the need for the tensor sum.

Thus, we have

\[\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)}.\]


To analyse the results, it is convenient to now move to a Python console or interactive environment, such as iPython or Jupyter Notebook. In the following, Python input will be specified with the prompt >>>, and output will follow.

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

As we are using the interferometer as in the Boson sampling tutorial, we do not have to recalculate the unitary, we can use the same \(U\) computed there, or copy the U defined in the file

\[\begin{split}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]\end{split}\]

Now that we have the interferometer unitary transformation \(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 haf

Now, for the right hand side numerator, we first calculate the submatrix \([(UU^T\tanh(r))]_{st}\):

>>> B = (, 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 \(\left|{1,1,0,0}\right\rangle\),

>>> B[:,[0, 1]][[0, 1]]
[[-0.10219728 + 0.32633851j,  0.55418347 + 0.28563583j],
[ 0.55418347 + 0.28563583j, -0.10505237 + 0.32960794j]]

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 \(\ket{0,0,0,0}\) at the output

    This corresponds to the hafnian of an empty matrix, which is simply 1:

    >>> 1 / np.cosh(1) ** 4

    Compare this to the Strawberry Fields result 0.176378447614135

  • Measuring \(\ket{1,1,0,0}\) at the output

    >>> B = (, U.T) * np.tanh(1))[:, [0, 1]][[0, 1]]
    >>> np.abs(haf(B)) ** 2 / np.cosh(1) ** 4

    Compare this to the Strawberry Fields result 0.0685595637122246

  • Measuring \(\ket{0,1,0,1}\) at the output

    >>> B = (, U.T) * np.tanh(1))[:, [1, 3]][[1, 3]]
    >>> np.abs(haf(B)) ** 2 / np.cosh(1) ** 4

    Compare this to the Strawberry Fields result 0.002056097258977398

  • Measuring \(\ket{1,1,1,1}\) at the output

    This corresponds to the hafnian of the full matrix \(B=UU^T\tanh(r)\):

    >>> B = (, U.T) * np.tanh(1))
    >>> np.abs(haf(B)) ** 2 / np.cosh(1) ** 4

    Compare this to the Strawberry Fields result 0.00834294639986785

  • Measuring \(\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 \(2!\):

    >>> B = (, U.T) * np.tanh(1))[:, [0, 0]][[0, 0]]
    >>> np.abs(haf(B)) ** 2 / (2 * np.cosh(1) ** 4)

    Compare this to the Strawberry Fields result 0.01031294525345511

They Strawberry Field simulation results agree (with almost negligible numerical error) to the expected result from the Gaussian boson sampling equation!


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 \(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.