Strawberry Fields:
A Software Platform for Photonic Quantum Computing
Nathan Killoran, Josh Izaac, Nicolás Quesada, Ville Bergholm, Matthew Amy, and
Christian Weedbrook
Xanadu, 372 Richmond St W, Toronto, M5V 1X6, Canada
We introduce Strawberry Fields, an open-source quantum programming architecture for light-based quantum computers,
and detail its key features. Built in Python, Strawberry Fields is a full-stack library for design, simulation, optimization, and
quantum machine learning of continuous-variable circuits. The platform consists of three main components: (i) an API for
quantum programming based on an easy-to-use language named Blackbird; (ii) a suite of three virtual quantum computer
backends, built in NumPy and TensorFlow, each targeting specialized uses; and (iii) an engine which can compile Black-
bird programs on various backends, including the three built-in simulators, and – in the near future – photonic quantum
information processors. The library also contains examples of several paradigmatic algorithms, including teleportation,
(Gaussian) boson sampling, instantaneous quantum polynomial, Hamiltonian simulation, and variational quantum circuit
optimization.
Introduction
The decades-long worldwide quest to build practical
quantum computers is currently undergoing a critical pe-
riod. During the next few years, a number of different
quantum devices will become available to the public. While
fault-tolerant quantum computers will one day provide sig-
nificant computational speedups for problems like factor-
ing [1], search [2], or linear algebra [3], near-term quan-
tum devices will be noisy, approximate, and sub-universal
[4]. Nevertheless, these emerging quantum processors are
expected to be strong enough to show a computational ad-
vantage over classical computers for certain problems, an
achievement known as quantum computational supremacy.
As we approach this milestone, work is already under-
way exploring how such quantum processors might best
be leveraged. Popular techniques include variational quan-
tum eigensolvers [5, 6], quantum approximate optimiza-
tion algorithms [7, 8], sampling from computationally hard
probability distributions [9–17], and quantum annealing
[18, 19]. Notably, these methodologies can be applied to
tackle important practical problems in chemistry [20–23],
finance [24], optimization [25, 26], and machine learning
[27–34]. These known applications are very promising, yet
it is perhaps the unknown future applications of quantum
computers that are most tantalizing. We may not know the
best applications of quantum computers until these devices
become available more widely to researchers, students, en-
trepreneurs, and programmers worldwide.
To this end, a nascent quantum software ecosystem has
recently begun to develop [35–47]. However, a prevail-
ing theme for these software efforts is to target qubit-
based quantum devices. In reality, there are several com-
peting models of quantum computing which are equiva-
lent in computational terms, but which are conceptually
quite distinct. One prominent approach is the continuous
variable (CV) model of quantum computing [48–50]. In
the CV model, the basic information-processing unit is an
infinite-dimensional bosonic mode, making it particularly
well-suited for implementations and applications based on
light. The CV model retains the computational power of
the qubit model (cf. Chap. 4 of Ref. [51]), while offering a
number of unique features. For instance, the CV model is a
natural fit for simulating bosonic systems (electromagnetic
fields, trapped atoms, harmonic oscillators, Bose-Einstein
condensates, phonons, or optomechanical resonators) and
for settings where continuous quantum operators – such as
position & momentum – are present. Even in classical com-
puting, recent advances from deep learning have demon-
strated the power and flexibility of a continuous represen-
tation of computation [52, 53] in comparison to the discrete
computational model which has historically dominated.
Here we introduce Strawberry Fields
1
, an open-source
software architecture for photonic quantum computing.
Strawberry Fields is a full-stack quantum software plat-
form, implemented in Python, specifically targeted to the
CV model. Its main element is a new quantum program-
ming language named Blackbird. To lay the groundwork
for future photonic quantum computing hardware, Straw-
1
This document refers to Strawberry Fields version 0.9. Full documenta-
tion is available online at strawberryfields.readthedocs.io.
arXiv:1804.03159v2 [quant-ph] 4 Mar 2019
2
berry Fields also includes a suite of three CV quantum simu-
lator backends implemented using NumPy [54] and Tensor-
Flow [55]. Strawberry Fields comes with a built-in engine
to convert Blackbird programs to run on any of the simula-
tor backends or, when they are available, on photonic quan-
tum computers. To accompany the library, an online service
for interactive exploration and simulation of CV circuits is
available at strawberryfields.ai.
Aside from being the first quantum software frame-
work to support photonic quantum computation with
continuous-variables, Strawberry Fields provides additional
computational features not presently available in the quan-
tum software ecosystem:
1. We provide two numeric simulators; a Gaussian back-
end, and a Fock-basis backend. These two formula-
tions are unique to the CV model of quantum compu-
tation due to the use of an infinite Hilbert space, and
came with their own technical challenges.
(a) The Gaussian backend provides state-of-the-art
methods and functions for calculating the fi-
delity and Fock state probabilities, involving cal-
culations of the classically intractable hafnian
[56].
(b) The Fock backend allows operations such as
squeezing and beamsplitters to be performed
in the Fock-basis, a computationally intensive
calculation that has been highly vectorized and
benchmarked for performance.
2. We provide a suite of important circuit decomposi-
tions appearing in quantum photonics – such as the
Williamson, Bloch-Messiah, and Clements decompo-
sitions.
3. The Fock-basis backend written using the TensorFlow
machine learning library allows for symbolic calcu-
lations, automatic differentiation, and backpropaga-
tion through CV quantum simulations. As far as we
are aware, this is the first quantum simulation library
written using a high-level machine learning library,
with support for dataflow programming and auto-
matic differentiation.
The remainder of this white paper is structured as fol-
lows. Before presenting Strawberry Fields, we first provide
a brief overview of the key ingredients for CV quantum com-
putation, specifically the most important states, gates, and
measurements. We then introduce the Strawberry Fields ar-
chitecture in full, presenting the Blackbird quantum assem-
bly language, outlining how to use the library for numerical
simulation, optimization, and quantum machine learning.
Finally, we discuss the three built-in simulators and the in-
ternal representations that they employ. In the Appendices,
we give further mathematical and software details and pro-
vide full example code for a number of important CV quan-
tum computing tasks.
Quantum Computation with Continuous
Variables
Many physical systems in nature are intrinsically contin-
uous, with light being the prototypical example. Such sys-
tems reside in an infinite-dimensional Hilbert space, offer-
ing a paradigm for quantum computation which is distinct
from the discrete qubit model. This continuous-variable
model takes its name from the fact that the quantum op-
erators underlying the model have continuous spectra. It
is possible to embed qubit-based computations into the CV
picture [57], so the CV model is as powerful as its qubit
counterparts.
From Qubits to Qumodes
A high-level comparison of CV quantum computation
with the qubit model is depicted in Table I. In the remainder
of this section, we will provide a basic presentation of the
key elements of the CV model. A more detailed technical
overview can be found in Appendix A. Readers experienced
with CV quantum computing can safely skip to the next sec-
tion.
CV Qubit
Basic element Qumodes Qubits
Relevant
operators
Quadratures
ˆ
x,
ˆ
p
Mode operators
ˆ
a,
ˆ
a
†
Pauli operators
ˆ
σ
x
,
ˆ
σ
y
,
ˆ
σ
z
Common states
Coherent states |α〉
Squeezed states |z〉
Number states |n〉
Pauli eigenstates
|0/1〉, |±〉, |±i〉
Common gates
Rotation,
Displacement,
Squeezing,
Beamsplitter, Cubic
Phase
Phase shift,
Hadamard, CNOT,
T-Gate
Common
measurements
Homodyne |x
φ
〉〈x
φ
|,
Heterodyne
1
π
|α〉〈α|,
Photon-counting |n〉〈n|
Pauli eigenstates
|0/1〉〈0/1|, |±〉〈±|,
|±i〉〈±i|
Table I: Basic comparison of the CV and qubit settings.
The most elementary CV system is the bosonic harmonic
oscillator, defined via the canonical mode operators
ˆ
a and
ˆ
a
†
. These satisfy the well-known commutation relation
[
ˆ
a,
ˆ
a
†
] = I. It is also common to work with the quadra-
ture operators (also called the position & momentum oper-
3
ators)
2
,
ˆ
x :=
v
t
ħh
2
(
ˆ
a +
ˆ
a
†
), (1)
ˆ
p := −i
v
t
ħh
2
(
ˆ
a −
ˆ
a
†
), (2)
where [
ˆ
x,
ˆ
p] = iħh
ˆ
I. We can picture a fixed harmonic oscil-
lator mode (say, within an optical fibre or a waveguide on a
photonic chip) as a single ‘wire’ in a quantum circuit. These
qumodes are the fundamental information-carrying units of
CV quantum computers. By combining multiple qumodes –
each with corresponding operators
ˆ
a
i
and
ˆ
a
†
i
– and interact-
ing them via sequences of suitable quantum gates, we can
implement a general CV quantum computation.
CV States
The dichotomy between qubit and CV systems is perhaps
most evident in the basis expansions of quantum states:
Qubit |φ〉 = φ
0
|0〉+ φ
1
|1〉, (3)
Qumode |ψ〉 =
Z
d x ψ(x) |x〉. (4)
For qubits, we use a discrete set of coefficients; for CV sys-
tems, we can have a continuum. The states |x〉are the eigen-
states of the
ˆ
x quadrature,
ˆ
x |x〉 = x |x〉, with x ∈ R. These
quadrature states are special cases of a more general family
of CV states, the Gaussian states, which we now introduce.
Gaussian states
Our starting point is the vacuum state |0〉. Other states
can be created by evolving the vacuum state according to
|ψ〉 = exp(−itH) |0〉, (5)
where H is a bosonic Hamiltonian (i.e., a function of the
operators
ˆ
a
i
and
ˆ
a
†
i
) and t is the evolution time. States
where the Hamiltonian H is at most quadratic in the oper-
ators
ˆ
a
i
and
ˆ
a
†
i
(equivalently, in
ˆ
x
i
and
ˆ
p
i
) are called Gaus-
sian. For a single qumode, Gaussian states are parameter-
ized by two continuous complex variables: a displacement
parameter α ∈ C and a squeezing parameter z ∈ C (of-
ten expressed as z = r exp(iφ), with r ≥ 0). Gaussian
states are so-named because we can identify each Gaussian
state with a corresponding Gaussian distribution. For single
qumodes, the identification proceeds through its displace-
ment and squeezing parameters. The displacement gives
2
It is common to picture ħh as a (dimensionless) scaling parameter for
the
ˆ
x and
ˆ
p operators rather than a physical constant. However, there
are several conventions for the scaling value in common use [58].
These self-adjoint operators are proportional to the Hermitian and anti-
Hermitian parts of the operator
ˆ
a. Strawberry Fields allows the user to
specify this value, with the default ħh = 2.
FIG. 1: Schematic representation of a Gaussian state for a
single mode. The shape and orientation are parameterized
by the displacement α and squeezing z = r exp(iφ).
the centre of the distribution, while the squeezing deter-
mines the variance and rotation of the distribution (see Fig.
1). Multimode Gaussian states, on the other hand, are pa-
rameterized by a vector of displacements
¯
r and a covariance
matrix V. Many important pure states in the CV model are
special cases of the pure Gaussian states; see Table II for a
summary.
State family Displacement Squeezing
Vacuum state |0〉 α = 0 z = 0
Coherent states |α〉 α ∈ C z = 0
Squeezed states |z〉 α = 0 z ∈ C
Displaced squeezed
states |α, z〉
α ∈ C z ∈ C
ˆ
x eigenstates |x〉
α ∈ C,
x = 2
q
ħh
2
Re(α)
φ = 0, r → ∞
ˆ
p eigenstates |p〉
α ∈ C,
p = 2
q
ħh
2
Im(α)
φ = π, r → ∞
Fock states |n〉 N.A. N.A.
Table II: Common single-mode pure states and their relation
to the displacement and squeezing parameters. All listed
families are Gaussian, except for the Fock states. The n = 0
Fock state is also the vacuum state.
Fock states
Complementary to the continuous Gaussian states are the
discrete Fock states (or number states) |n〉, where n are non-
negative integers. These are the eigenstates of the number
operator
ˆ
n =
ˆ
a
†
ˆ
a. The Fock states form a discrete (count-
able) basis for qumode systems. Thus, each of the Gaussian
states considered in the previous section can be expanded
in the Fock-state basis. For example, coherent states have
4
the form
|α〉 = exp
−
|α|
2
2
∞
X
n=0
α
n
p
n!
|n〉, (6)
while (undisplaced) squeezed states only have even number
states in their expansion:
|z〉 =
1
p
cosh r
∞
X
n=0
p
(2n)!
2
n
n!
[−e
iφ
tanh(r)]
n
|2n〉. (7)
Mixed states
Mixed Gaussian states are also important in the CV pic-
ture, for instance, the thermal state
ρ(n) :=
∞
X
n=0
n
n
(1 + n)
n+1
|n〉〈n|, (8)
which is parameterized via the mean photon number n :=
Tr(ρ(n)
ˆ
n). Starting from this state, we can consider a
mixed-state-creation process similar to Eq. (5), namely
ρ = exp(−itH)ρ(n) exp( i tH). (9)
Analogously to pure states, by applying Hamiltonians of
second-order (or lower) to thermal states, we generate the
family of Gaussian mixed states.
CV Gates
Unitary operations can be associated with a generating
Hamiltonian H via the recipe (cf. Eqs. (5) & (9))
U := exp (−i tH). (10)
For convenience, we classify unitaries by the degree of their
generator. A CV quantum computer is said to be universal
if it can implement, to arbitrary precision and with a finite
number of steps, any unitary which is polynomial in the
mode operators [48]. We can build a multimode unitary by
applying a sequence of gates from a universal gate set, each
of which acts only on one or two modes. We focus on a
universal set made from the following two subsets:
Gaussian gates: Single and two-mode gates which are at
most quadratic in the mode operators, e.g., Displace-
ment, Rotation, Squeezing, and Beamsplitter gates.
Non-Gaussian gate: A single-mode gate which is degree 3
or higher, e.g., the Cubic phase gate.
A number of fundamental CV gates are presented in Ta-
ble III. Many of the Gaussian states from the previous sec-
tion are connected to a corresponding Gaussian gate. Any
multimode Gaussian gate can be implemented through a
suitable combination of Displacement, Rotation, Squeezing,
and Beamsplitter Gates [50], making these gates sufficient
for constructing all quadratic unitaries. The cubic phase
gate is presented as an exemplary non-Gaussian gate, but
any other non-Gaussian gate could also be used to achieve
universality. A number of other useful CV gates are listed in
Appendix B.
Gate Unitary Symbol
Displacement D
i
(α) = exp (α
ˆ
a
†
i
−α
∗
ˆ
a
i
)
D
Rotation R
i
(φ) = exp (iφ
ˆ
n
i
)
R
Squeezing S
i
(z) = exp (
1
2
(z
∗
ˆ
a
2
i
−z
ˆ
a
†2
i
))
S
Beamsplitter
BS
i j
(θ, φ) =
exp (θ(e
iφ
ˆ
a
i
ˆ
a
†
j
− e
−iφ
ˆ
a
†
i
ˆ
a
j
))
BS
Cubic phase
V
i
(γ) = exp
i
γ
3ħh
ˆ
x
3
i
V
Table III: Some important CV model gates. All listed gates
except the cubic phase gate are Gaussian.
CV Measurements
As with CV states and gates, we can distinguish between
Gaussian and non-Gaussian measurements. The Gaus-
sian class consists of two (continuous) types: homodyne
and heterodyne measurements, while the key non-Gaussian
measurement is photon counting. These are summarized in
Table IV.
Homodyne measurements
Ideal homodyne detection is a projective measurement
onto the eigenstates of the quadrature operator
ˆ
x. These
states form a continuum, so homodyne measurements are
inherently continuous, returning values x ∈ R. More gen-
erally, we can consider projective measurement onto the
eigenstates
x
φ
of the Hermitian operator
ˆ
x
φ
:= cos φ
ˆ
x + sin φ
ˆ
p. (11)
This is equivalent to rotating the state clockwise by φ and
performing an
ˆ
x-homodyne measurement. If we have a
multimode Gaussian state and we perform homodyne mea-
surement on one of the modes, the conditional state of the
unmeasured modes remains Gaussian.
Heterodyne measurements
Whereas homodyne detection is a measurement of
ˆ
x,
heterodyne detection can be seen as a simultaneous mea-
surement of both
ˆ
x and
ˆ
p. Because these operators do not
commute, they cannot be simultaneously measured without
some degree of uncertainty. Equivalently, we can picture
heterodyne measurement as projection onto the coherent
states, with measurement operators
1
π
|α〉〈α|. Because the
coherent states are not orthogonal, there is a corresponding
lack of sharpness in the measurements. If we perform het-
erodyne measurement on one mode of a multimode state,
the conditional state on the remaining modes stays Gaus-
sian.
5
Measurement
Measurement
Operators
Measurement
values
Homodyne
|x
φ
〉〈x
φ
|
x ∈ R
Heterodyne
1
π
|α〉〈α|
α ∈ C
Photon counting |n〉〈n| n ∈ N
Table IV: Key measurement types for the CV model. The
‘-dyne’ measurements are Gaussian, while photon-counting is
non-Gaussian.
Photon Counting
Photon counting (also known as as photon-number re-
solving measurement), is a complementary measurement
method to the ‘-dyne’ measurements, revealing the particle-
like, rather than the wave-like, nature of qumodes. This
measurement projects onto the number eigenstates |n〉, re-
turning non-negative integer values n ∈ N. Except for the
outcome n = 0, a photon-counting measurement on a sin-
gle mode of a multimode Gaussian state will cause the re-
maining modes to become non-Gaussian. Thus, photon-
counting can be used as an ingredient for implementing
non-Gaussian gates. A related process is photodetection,
where a detector only resolves the vacuum state from non-
vacuum states. This process has only two measurement op-
erators, namely |0〉〈0| and I −|0〉〈0|.
The Strawberry Fields Software Platform
The Strawberry Fields library has been designed with sev-
eral key goals in mind. Foremost, it is a standard-bearer for
the CV model, laying the groundwork for future photonic
quantum computers. As well, Strawberry Fields is designed
to be simple to use, giving entry points for as many users
as possible. Finally, since the potential applications of near-
term quantum computers are still being worked out, it is
important that Strawberry Fields provides powerful tools
to easily explore many different use-cases and applications.
Strawberry Fields has been implemented in Python, a
modern language with a gentle learning curve which is al-
ready familiar to many programmers and scientific practi-
tioners. The accompanying quantum simulator backends
are built upon the widely used Python packages NumPy
and TensorFlow. All Strawberry Fields code is open source.
Strawberry Fields can be accessed programmatically as a
Python package, or via a browser-based interface for de-
signing quantum circuits.
A pictorial outline of Strawberry Fields’ key elements and
their interdependencies is presented in Fig. 2. Conceptu-
ally, the software stack is separated into two main pieces: a
user-facing frontend layer and a lower-level backends com-
ponent. The frontend encompasses the Strawberry Fields
Python API and the Blackbird quantum assembly language.
These elements provide access points for users to design
quantum circuits. These circuits are then linked to a back-
end via a quantum compiler engine. For a backend, the
engine can currently target one of three included quantum
computer simulators. When CV quantum processors be-
come available in the near future, the engine will also build
and run circuits on those devices. Further, high-level quan-
tum computing applications can be built by leverging the
Strawberry Fields frontend API. Existing examples include
the Strawberry Fields Interactive website, the Quantum Ma-
chine Learning Toolbox (for streamlining the training of
variational quantum circuits), and SFOpenBoson (an inter-
face between the electronic structure library OpenFermion
[45] and Strawberry Fields).
In the remainder of this section, the key elements of
Strawberry Fields will be presented in more detail. Proceed-
ing through a series of examples, we show how CV quantum
computations can be defined using the Blackbird language,
then compiled and run on a quantum computer backend.
We also outline how to use Strawberry Fields for optimiza-
tion and machine learning on quantum circuits. Finally, we
discuss the suite of quantum computer simulators included
within Strawberry Fields.
Blackbird: A Quantum Programming Language
As classical computers have become progressively more
powerful, the languages used to program them have also
undergone considerable paradigmatic changes. Machine
code gave way to human-readable assembly languages, fol-
lowed by higher-level procedural and object-oriented lan-
guages. With each generation, the trend has been towards
higher levels of abstraction, separating the programmer
more and more from details of the actual computer hard-
ware. Quantum computers are still at an early stage of
development, so while we can imagine what higher-level
quantum programming might look like, in the near term we
first need to build languages which are conceptually closer
to the quantum hardware.
Blackbird is a standalone domain specific language (DSL)
for continuous-variable quantum computation. With a well-
defined grammar in extended Backus-Naur form, and both
Python and C++ parsers available, Blackbird provides op-
erations that match the basic CV states, gates, and mea-
surements, and maps directly to low-level hardware instruc-
tions. The abstract syntax keeps a close connection between
code and the quantum operations that they implement; this
syntax is modeled after that of ProjectQ [41], but special-
ized to the CV setting. Blackbird can be used as part of
the Strawberry Fields stack, but also directly with photonic
quantum computing hardware systems.
Within the Strawberry Fields framework, we have built
an implementation of Blackbird using Python 3 as the em-
bedding language — an ‘embedded’ DSL. This ‘Python-
enhanced’ Blackbird language provides the same core oper-
ations and follows the same grammar and syntactical rules
as the standalone DSL, but, by nature, may also contain
valid Python constructs. Furthermore, Strawberry Fields’
‘Python-enhanced’ Blackbird provides users with additional
6
COMPILER
ENGINE
BACK
-
ENDS
FRONT
-
ENDS
Strawberry Fields
Interactive
Serve
r
Applications
Strawberry
Fields API
Blackbird
Quantum
Programming
Language
Quantum Processor
Simulators
Fock Representation
Gaussian Representation
NumPy
Tensorflow
NumPy
FIG. 2: Outline of the Strawberry Fields software stack. The Strawberry Fields Interactive server is available online at
strawberryfields.ai.
quantum operations that are decomposed into lower-level
Blackbird assembly commands. We will introduce the ele-
ments of Blackbird through a series of basic examples, dis-
cussing more technical aspects as they arise.
Operations
Quantum computations consist of four main ingredients:
state preparations, application of gates, performing mea-
surements, and adding/removing subsystems. In Blackbird,
these are all considered as Operations, and share the same
basic syntax. In the following code examples, we use the
variable q for a set of qumodes (more specifically, a quan-
tum register), the details of which are deferred until the next
section.
Our first considered Operation is state preparation. By
default, qumodes are initialized in the vacuum state. Var-
ious other important CV states can be created with simple
Blackbird commands.
1 # Create the vacuum state in qumode 0
2 Vac | q[0]
3
4 # Create a coherent state in qumode 1
5 alpha = 2.0 + 1j
6 Coherent(alpha) | q[1]
7
8 # Create squeezed states in qumodes 0 & 1
9 S = Squeezed(2.0)
10 S | q[0]
11 S | q[1]
12
13 # Create a Fock state in qumode 1
14 Fock(4) | q[1]
Codeblock 1: Blackbird code for creating various CV
quantum states.
Blackbird state preparations such as those used in Code-
block 1 implicitly reset the existing state of the qumodes.
Conceptually, the vertical bar symbol ‘|’ separates Opera-
tions – like state preparation – from the registers that they
act upon. Notice that we can use Operations inline, or con-
struct them separately and reuse them several times.
After creating states, we will want to transform these us-
ing quantum gates.
1 # Apply the Displacement gate to qumode 0
2 alpha = 2.0 + 1j
3 Dgate(alpha) | q[0]
4
5 # Apply the Rotation gate
6 phi = 1.157
7 Rgate(phi) | q[0]
8
9 # Apply the Squeezing gate
10 Sgate(2.0, 0.17) | q[0]
11
12 # Apply the Beamsplitter gate to qumodes 0 & 1
13 BSgate(0.314, 0.223) | (q[0], q[1])
14
15 # Apply Cubic phase gate (VGate) to qumode 0
16 gamma = 0.1
17 Vgate(gamma) | q[0]
18
19 # Apply Hermitian conjugate of a gate
20 V = Vgate(gamma)
21 V.H | q[0]
Codeblock 2: Blackbird code for applying various CV gates.
Blackbird supports all of the gates listed in the previous sec-
tion as well as a number of composite gates, each of which
can be decomposed using the universal gates. The sup-
ported composite gates are: controlled X (CXgate), con-
trolled Z (CZgate), quadratic phase (Pgate), and two-
mode squeezing (S2gate). A full list of gates currently
supported in Blackbird can be found in Appendix B.
Finally, we can specify measurement Operations using
Blackbird.
7
1 # Homodyne measurement at angle phi
2 phi = 0.785
3 MeasureHomodyne(phi) | q[0]
4
5 # Special homodyne measurements
6 MeasureX | q[0]
7 MeasureP | q[1]
8
9 # Heterodyne measurement
10 MeasureHeterodyne() | q[0]
11 MeasureHD | q[1] # shorthand
12
13 # Number state measurements of various qumodes
14 MeasureFock() | q[0]
15 MeasureFock() | (q[1], q[2]) # multiple modes
16 Measure | q[3] # shorthand
Codeblock 3: Blackbird code for carrying out CV
measurements.
Measurements have several effects. For one, the numeri-
cal result of the measurement is placed in a classical regis-
ter. As well, the state of all remaining qumodes is projected
to the (normalized) conditional state for that measurement
value. Finally, the state of the measured qumode is reset
to the vacuum state. This is the typical behaviour of pho-
tonic hardware, where measurements absorb all the energy
of the measured qumode.
Running Blackbird Programs in Python
Within Python, Blackbird programs are managed by an
Engine. The function Engine in the Strawberry Fields API
will instantiate an Engine, returning both the Engine and
its corresponding quantum register. The Engine is used as
a Python context manager, providing a convenient way to
encapsulate Blackbird programs.
1 # Create Engine and quantum register
2 import strawberryfields as sf
3 eng, q = sf.Engine(num_subsystems=2)
4
5 # The register is also available via
eng.register,→
6 assert q == eng.register
7
8 # Put Blackbird Operations in namespace
9 from strawberryfields.ops import *
10
11 # Declare a Blackbird program
12 from math import pi
13 z = 4.
14 S = Sgate(z)
15 B = BSgate(pi / 4, 0)
16 with eng:
17 S | q[0]
18 S.H | q[1]
19 B | q
20 MeasureP | q[0]
21 MeasureP | q[1]
22
23 # Execute Blackbird program and extract values
24 eng.run(backend="gaussian")
25 vals = [reg.val for reg in q]
Codeblock 4: Code for declaring and running Blackbird
programs using the Strawberry Fields library.
The above code example is runnable and carries out a com-
plete quantum computation, namely the preparation and
measurement of an EPR entangled state. Note that Oper-
ations can be declared outside of the Engine, but their ac-
tion on the quantum registers must come within the En-
gine context. Also notice that our register has a length
of 2, so any single-mode Operations must act on specific
elements, i.e., q[i], while two-mode Operations can act
on q directly. Finally, the user must specify a backend –
as well as any backend-dependent settings – when calling
eng.run(). We will discuss the Strawberry Fields back-
ends further in a later section.
Quantum and Classical Registers
When a Strawberry Fields Engine is constructed, the
user must specify the number of qumode subsystems to be-
gin with. This number is required for the initialization of
the Engine, but may change within a computation (e.g.,
when temporary ancilla modes are used). Qumodes can
be added/deleted by using the New and Del Operations.
1 # A Blackbird circuit where gates
2 # are added and deleted
3 alice = q[0]
4 with eng:
5 Sgate(1) | alice
6 bob, charlie = New(2)
7 BSgate(0.5) | (alice, bob)
8 CXgate(1) | (alice, charlie)
9 Del | alice
10 S2gate(0.4) | (charlie, bob)
11
12 # Attempting to act on registers which have
13 # been removed will raise an IndexError
14 try:
15 with eng:
16 Dgate(0.1) | alice
17 except Exception as e:
18 assert isinstance(e, IndexError)
Codeblock 5: Adding and deleting qumode subsystems.
An Engine maintains a unique numeric indexing for the
quantum registers based on the order they were added.
When a subsystem is deleted from the circuit, no further
gates can act on that register.
As stated earlier, measurement Operations produce clas-
sical information (the measurement result) and manipu-
late the corresponding register. Compared to previously
released quantum programming frameworks, such as Pro-
jectQ, PyQuil, and Qiskit, Strawberry Fields has been de-
signed so that the notation q[i], while principally denoting
the quantum register, may also encapsulate classical infor-
mation or a classical register. This behaviour is contextual,
8
and unique to Strawberry Fields. For example, prior to mea-
surement, q[i] simply references a quantum register. Once
a measurement operation is performed, q[i] continues to
represent a quantum register — now reset to the vacuum
state — as well as storing the numerical value of the mea-
surement, accessible via the attribute q[i].val. Note that
this numerical value is only available if a computation has
been run up to the point of measurement. We may also use
a classical measurement result symbolically as a parameter
in later gates without first running the computation. To do
this, we simply pass the measured register (e.g., q[i]) ex-
plicitly as an argument to the required gate. As before, the
Strawberry Fields quantum register object is contextual —
when passed as a gate argument, Strawberry Fields implic-
itly accesses the encapsulated classical register.
1 # Numerical evaluation of a measurement
2 # result using eng.run()
3 with eng:
4 MeasureX | q[0]
5 eng.run("gaussian")
6 val = q[0].val
7
8 # Use a measured register symbolically
9 # in another gate
10 with eng:
11 MeasureX | q[0]
12 Dgate(q[0]) | q[1]
13 eng.run("gaussian")
Codeblock 6: Evaluating measurement results numerically
and using them symbolically.
In quantum algorithms, it is common to process a measure-
ment result classically and use the post-processed value as
a parameter for further operations in a circuit. Strawberry
Fields provides the convert decorator to transform a user-
specified numerical function into one which acts on regis-
ters.
1 @sf.convert
2 def neg(x):
3 return -x
4
5 # A Blackbird computation using classical
6 # data processing
7 with eng:
8 MeasureX | q[0]
9 Dgate(neg(q[0])) | q[1]
10 eng.run("gaussian")
Codeblock 7: Symbolically processing a measured value
before using it.
Post-selection
The measurement Operations in Strawberry Fields are
stochastic in nature, with outcomes determined by some
underlying quantum probability distribution. Often it is
convenient to select specific values for these measurements
rather than sampling them. For instance, we might want to
explore the conditional state created by a specific value, de-
termine the measurement-dependent corrections we need
to make in a teleportation circuit, or even design an algo-
rithm which inherently contains post-selection. This func-
tionality is supported in Strawberry Fields through the op-
tional keyword argument select, which can be supplied
for any measurement Operation. The measurement out-
come will then return exactly this value, while the remain-
ing modes will be projected into the conditional state cor-
responding to this value
3
.
1 with eng:
2 Fock(3) | q[0]
3 Fock(2) | q[1]
4 BSgate() | (q[0], q[1])
5 MeasureFock(select=4) | q[0]
6 MeasureFock() | q[1]
7 eng.run("fock", cutoff_dim=6)
8 assert q[0].val == 4
Codeblock 8: Selecting a specific desired measurement
outcome.
Decompositions
In addition to the core CV operations discussed above,
Strawberry Fields also provides support for some impor-
tant decompositions frequently used in quantum optics.
These include the (a) Williamson decomposition [59],
for decomposing arbitrary Gaussian states to a symplectic
transformation acting on a thermals state, (b) the Bloch-
Messiah decomposition [60–62], for decomposing the ac-
tion of symplectic transformations to interferometers and
single-mode squeezing, and (c) the Clements decomposi-
tion [63], for decomposing multi-mode linear interferom-
eters into arrays of beamsplitters and rotations of fixed
depth. In all cases, the resulting decomposition into the
universal CV gate set may be viewed via the engine method
eng.print_applied(). Strawberry Fields thus provides
a natural environment for embedding graphs and matrices
in quantum optical circuits, and viewing the resulting phys-
ical components.
1 U = np.array([[1-1j, np.sqrt(2)],
2 [-np.sqrt(2), 1+1j]])/2
3
4 eng, q = sf.Engine(2)
5
6 with eng:
7 Squeezed(0.43) | q[0]
8 Interferometer(U) | (q[0], q[1])
9
10 eng.run("gaussian")
11 eng.print_applied()
3
Users should be careful to avoid post-selection on measurement values
which have no probability of occuring given the current circuit state. In
this case, the expected behaviour of a backend is not defined.
9
12 # >> Squeezed(0.43, 0) | (q[0])
13 # >> Rgate(2.356) | (q[0])
14 # >> BSgate(0.7854, 0) | (q[0], q[1])
15 # >> Rgate(-3.142) | (q[0])
16 # >> Rgate(0.7854) | (q[1])
Codeblock 9: Using the in-built Clements decomposition to
decompose a 2 ×2 Interferometer into beamsplitters
and phase rotations.
Optimization and Quantum Machine Learning
Strawberry Fields can perform quantum circuit simula-
tions using both numerical and symbolic representations.
Numerical computation is the default operating mode and
is supported by all three supplied backends. Symbolic com-
putation is enabled only for the TensorFlow backend. In
this section, we outline the main TensorFlow functionali-
ties accessible through the Strawberry Fields frontend in-
terface. More details about the corresponding TensorFlow
backend are discussed in the next section. TensorFlow [55]
models calculations abstractly using a computational graph,
where individual operations are represented as nodes and
their dependencies by directed edges. This viewpoint sep-
arates the symbolic representation of a computation from
its numerical evaluation, and makes optimization and ma-
chine learning more amenable. On top of this, TensorFlow
provides a number of advanced functionalities, including
automatic gradient computation, GPU utilization, built-in
optimization algorithms, and various other machine learn-
ing tools. Note that the term ‘quantum machine learning’
will be used here in a hybrid sense, i.e., applying conven-
tional machine learning methods to quantum systems.
To build a TensorFlow computational graph using Straw-
berry Fields, we instantiate an Engine, declare a circuit in
Blackbird code, then execute eng.run on the TensorFlow
("tf") backend. To keep the underlying simulation fully
symbolic, the extra argument eval=False must be given.
1 # Create Engine and run symbolic computation
2 import strawberryfields as sf
3 from strawberryfields.ops import *
4 import tensorflow as tf
5
6 eng, q = sf.Engine(num_subsystems=1)
7 with eng:
8 Dgate(0.5) | q[0]
9 MeasureX | q[0]
10 eng.run("tf",
11 cutoff_dim=5,
12 eval=False)
13
14 # Registers contain symbolic measurements
15 print(q[0].val)
16 # >> Tensor("Measure_homodyne/Meas_result:0",
shape=(), dtype=float64),→
Codeblock 10: Creating a TensorFlow computational graph
for a quantum circuit.
When we do this, any registers measured in the circuit will
be populated with unevaluated Tensor objects rather than
numerical values (without eval=False, the TensorFlow
backend returns purely numerical results, similar to the
other simulators). These Tensors can still be evaluated nu-
merically by running them in a TensorFlow Session. In
this case, measurement results will be determined stochas-
tically on each evaluation.
18 # Evaluate Tensors numerically
19 sess = tf.Session()
20 sess.run(tf.global_variables_initializer())
21 print(sess.run(q[0].val))
22 # >> a numerical measurement value
23 print(sess.run(q[0].val))
24 # >> a (different) numerical measurement value
Codeblock 11: Numerically evaluating Tensors.
When specifying a circuit in Blackbird, we can make
use of various special symbolic TensorFlow classes, such as
Variable, Tensor, placeholder, or constant.
26 # Declare circuits using Tensorflow objects
27 alpha = tf.Variable(0.1)
28 phi = tf.constant(0.5)
29 theta_bs = tf.placeholder(tf.float64)
30 phi_bs = tf.nn.relu(phi) # a Tensor object
31
32 eng, q = sf.Engine(num_subsystems=2)
33 with eng:
34 Coherent(alpha) | q[0]
35 Measure | q[0]
36 BSgate(theta_bs, phi_bs) | (q[0], q[1])
37 MeasureHomodyne(phi) | q[1]
38 eng.run("tf",
39 cutoff_dim=5,
40 eval=False)
41
42 sess = tf.Session()
43 sess.run(tf.global_variables_initializer())
44 feed_dict = {theta_bs: 0.5, q[0].val: 2}
45 print(sess.run(q[1].val, feed_dict=feed_dict))
46 # >> a numerical measurement value
Codeblock 12: Using abstract TensorFlow classes as circuit
parameters.
In the above example, we supplied an additional
feed_dict argument when evaluating. This is a Python
dictionary which specifies the numerical values (typically,
coming from a dataset) for every placeholder that
appears in a circuit. However, as can be seen from the
example, it is also possible to substitute desired values for
other nodes in the computation, including the values stored
in quantum registers. This allows us to easily post-select
measurement values and explore the resulting conditional
states.
We can also perform post-processing of measurement re-
sults when working with TensorFlow objects. In this case,
the functions decorated by convert should be written us-
ing TensorFlow Tensors, Variables, and operations.
10
48 # Use a simple neural network as
49 # a processing function
50 @sf.convert
51 def NN(x):
52 weight = tf.Variable(0.5, dtype=tf.float64)
53 bias = tf.Variable(0.1, dtype=tf.float64)
54 return tf.sigmoid(weight * x + bias)
55
56 eng, q = sf.Engine(num_subsystems=2)
57 with eng:
58 MeasureP | q[0]
59 Dgate(NN(q[0])) | q[1]
60 MeasureX | q[1]
61 eng.run("tf",
62 cutoff_dim=5,
63 eval=False)
64 print([r.val for r in q])
65 # >> [<tf.Tensor
'Measure_homodyne_2/Meas_result:0'
shape=() dtype=float64>, <tf.Tensor
'Measure_homodyne_3/Meas_result:0'
shape=() dtype=float64>]
,→
,→
,→
,→
66
Codeblock 13: Processing a measurement result using a
neural network. For compactness, the example uses a width
1 perceptron, but any continuous processing function
supported by TensorFlow can be used.
The TensorFlow backend additionally supports the use
of batched processing, allowing for many evaluations of
a quantum circuit to potentially be computed in paral-
lel. Scalars are automatically broadcast to the specified
batch size. Finally, we can easily run circuit simulations on
special-purpose hardware like GPUs or TPUs.
68 batch_size = 3
69 eng, q = sf.Engine(num_subsystems=1)
70
71 alpha = tf.Variable([0.1] * batch_size)
72 # scalars are automatically cast to batch size
73 beta = tf.Variable(0.5)
74
75 sess = tf.Session()
76 sess.run(tf.global_variables_initializer())
77
78 with tf.device("/gpu:0"):
79 with eng:
80 Dgate(alpha) | q[0]
81 Dgate(beta) | q[0]
82 MeasureX | q[0]
83 eng.run("tf",
84 cutoff_dim=10,
85 batch_size=batch_size,
86 eval=False)
87 print(q[0].val)
88 # >>
Tensor("Measure_homodyne_4/Meas_result:0",
shape=(3,), dtype=float64,
device=/device:GPU:0)
,→
,→
,→
Codeblock 14: Running a batched computation and
explicitly placing the computation on a GPU.
By taking advantage of these additional functionalities of
the TensorFlow backend, we can straightforwardly perform
optimization and machine learning on quantum circuits in
Strawberry Fields [64, 65]. A complete code example for
optimization of a quantum circuit is located in Appendix C.
Strawberry Fields’ Quantum Simulators
The ultimate goal is for Blackbird programs to be carried
out on photonic quantum computers. To lay the ground-
work for these emerging devices, Strawberry Fields comes
with a suite of three CV quantum computer simulators spe-
cially designed for the CV model. These simulators tar-
get different use-cases and support different functionality.
For example, many important algorithms in the CV formal-
ism involve only Gaussian states, operations, and measure-
ments. We can take advantage of this structure to more
efficiently simulate such computations. Other circuits have
inherently non-Gaussian elements to them; for these, the
Fock basis provides the standard description. These repre-
sentations are available in Strawberry Fields in the Gaussian
backend and the Fock backend, respectively. The third built-
in backend is the TensorFlow backend. Also using the Fock
representation, this backend is geared primarily towards
optimization and machine learning applications.
Most Blackbird operations are supported across all three
backends. A small subset, however, are not supported uni-
formly due to mathematical incompatibility. For example,
the Gaussian backend does not generally support non-
Gaussian gates or the preparation/measurement of Fock
states. Sometimes it is also useful to work with a backend
directly. To allow this, Strawberry Fields provides a back-
end API, giving access to additional methods and properties
which are not part of the streamlined frontend API. A stan-
dalone backend can be created in Strawberry Fields using
strawberryfields.backend.load_backend(name)
where name is one of "gaussian", "fock", or "tf" (for
comparison, the backend associated to an Engine eng is
available via eng.backend).
Three important backend methods, common to all simu-
lators, are begin_circuit, reset, and state. The first
command instantiates the circuit simulation, the second re-
sets the simulation back to an initial vacuum state, clearing
all previous operations, and the third returns a class which
encapsulates the current quantum state of the simulator.
In addition to containing the numerical (or symbolic) state
data, state classes also contain a number of useful meth-
ods and attributes for further exploring the quantum state,
such as fidelity, mean_photon, or wigner. As a conve-
nience for the user, all simulations carried out via eng.run
will return a state class representing the final circuit state
(see Appendix C for examples).
Gaussian Backend
This backend, written in NumPy, uses the symplectic for-
malism to represent CV systems. At a high level, this repre-
11
sentation tracks the quantum state of an N -mode quantum
system using two Gaussian components: a 2N-dimensional
displacement vector
¯
r and a 2N × 2N-dimensional covari-
ance matrix V (a deeper technical overview is located
in Appendix A). After we have created a Gaussian state
(either via state = eng.run(backend="gaussian")
or by directly calling the state method of a Gaussian
backend), we can access
¯
r and V via state.means and
state.cov, respectively. Other useful Gaussian state
methods are displacement and squeezing, which re-
turn the Gaussian parameters associated to the underlying
state.
The scaling of the symplectic representation with the
number of modes is O(N
2
). On one hand, this is quite pow-
erful. It allows us to to efficiently simulate any computa-
tions which are fully Gaussian. On the other, the formalism
is not expressive enough to simulate more general quan-
tum computations. Only a small number of non-Gaussian
methods are available for this backend. These are auxiliary
methods where we extract some non-Gaussian information
from a Gaussian state, but do not update the state of the
circuit. One such method is fock_prob, which is imple-
mented using an optimized – yet still exponentially scaling –
algorithm. This method enables simulation of the Gaussian
boson sampling algorithm [10] using the Gaussian backend;
see Appendix C for a complete code example.
Fock Backend
This backend, also written in NumPy, uses a fundamen-
tally different description for qumodes than the Gaussian
representation. As discussed in the introductory sections,
the Fock representation encodes quantum computation in
a countably infinite-dimensional Hilbert space. This repre-
sentation is faithful, allowing a precise description of CV
systems, in particular non-Gaussian circuits. Yet simulating
infinite-dimensional systems leads to some computational
tradeoffs which are not present for qubit simulators. Most
importantly, we impose a cutoff dimension D for simula-
tions (chosen by the user), so the Fock backend only covers
a restricted set of number states {|0〉, ..., |D −1〉} for each
mode. The size of simulated quantum systems thus depends
on both the number of subsystems N and the cutoff, being
O(D
N
). We contrast this with qubit systems, where the base
is fixed, i.e., O(2
N
). While these scalings are both expo-
nential, in practice simulating qumode systems for D > 2 is
more computationally demanding than qubits. This is be-
cause the (truncated) qumode subsystems have a higher di-
mension and thus encode more information than their qubit
counterparts.
Physically, imposing a cutoff is a reasonable strategy since
higher photon-number states must have higher energy and,
in practice, quantum-optical systems will have bounded en-
ergy (e.g., limited by the power of a laser). On the other
hand, there are certainly states which can be easily pre-
pared in the lab, yet would not fit accurately on the sim-
ulator. Thus, some care must be taken to trade off between
the numerical cutoff value, the number of modes, and the
energy scale of the circuit. If the energy scale is sufficiently
low that all states fit within the specified cutoff, then sim-
ulations with the Fock and Gaussian backends will be in
numerical agreement.
Like the Gaussian backend, the Fock backend has a
state method which encapsulates the numerical state,
while also providing a number of methods and attributes
specific to the Fock representation (such as ket, trace,
and all_fock_probs.). Unlike the Gaussian representa-
tion, mixed state simulations take up more resources than
pure states. Pure states are represented in the Fock back-
end by an D
N
-dimensional complex vector and mixed states
by a D
2N
-dimensional density matrix. Because of this extra
overhead, by default the Fock backend will internally rep-
resent a quantum circuit as long as possible as a pure state,
switching to the mixed state representation only when it
becomes necessary. Most importantly, for N > 2 qumodes,
all state preparation Operations (Vacuum, Squeezed, Fock,
etc.) cause the representation to become mixed. This is be-
cause the mode where the state is prepared could be entan-
gled with other modes. To keep physically consistent, the
Fock backend will first trace out the relevant mode, neces-
sitating a mixed state representation. When possible, it is
recommended to apply gates to the (default) vacuum state in
order to efficiently prepare pure states. If desired, a mixed
state simulation can be enforced by passing the argument
pure=False when calling begin_circuit.
TensorFlow Backend
The other built-in backend for Strawberry Fields is coded
using TensorFlow. As a simulator, it uses the same internal
representation as the Fock backend (Fock basis, finite cut-
off, pure vs. mixed state representations, etc.) and has the
same methods. It can operate as a numerical simulator sim-
ilar to the other backends. Its main purpose, however, is to
leverage the many powerful tools provided by TensorFlow
to enable optimization and machine learning on quantum
circuits. Much of this functionality was presented in the
previous section, so we will not repeat it here.
Like the other simulators, users can query the TensorFlow
backend’s state method to access the internal representa-
tion of a circuit’s quantum state. This functions similarly
to the Fock backend’s state method, except that the state
returned can be an unevaluated Tensor object when the
keyword argument eval is set to False. This state Ten-
sor can be combined with any supported TensorFlow op-
erations (norm, self_adjoint_eig, inv, etc.) to enable
optimization and machine learning on various properties of
quantum circuits and quantum states.
Conclusions
We have introduced Strawberry Fields, a multi-faceted
software platform for continuous-variable quantum com-
puting. The main components of this library – a custom
quantum programming language (Blackbird), a compiler
engine, and a suite of quantum simulators targeting distinct
12
applications – have been presented in detail. Further infor-
mation is available in both the Appendices and the Straw-
berry Fields online documentation.
The stage is now set for the broader community to use
Strawberry Fields for exploration, research, and develop-
ment of new quantum algorithms, specialized circuits, and
machine learning models. We anticipate the creation of
further software applications and backend modules for the
Strawberry Fields platform (developed both internally and
externally), providing advanced functionality and applica-
tions for quantum computing and quantum machine learn-
ing.
Acknowledgements
We thank our colleagues at Xanadu for testing Strawberry
Fields, reviewing this white paper, and providing helpful
feedback. In particular, we thank Patrick Rebentrost for
valuable discussions and suggestions.
13
Appendix A: The CV model
In this Appendix, we provide more technical and mathe-
matical details about the CV model, the quantum comput-
ing paradigm underlying Strawberry Fields.
Universal Gates for CV Quantum Computing
In discrete qubit systems, the notion of a universal gate set
has the following meaning: given a set of universal gates,
we can approximate an arbitrary unitary by composition of
said gates. In the CV setting, we have a similar situation: we
can approximate a broad set of generators – i.e., the Hamil-
tonians appearing in Eq. (10) – by combining elements of
a CV universal gate set. However, unlike the qubit case, we
do not try to approximate all conceivable unitaries. Rather,
we seek to create all generators that are a polynomial func-
tion of the quadrature (or mode) operators of the system
[48, 49]. Remember that generators of second-degree or
lower belong to the class of Gaussian operations, while all
higher degrees are non-Gaussian.
We can create a higher-order generator out of lower-
order generators
ˆ
A and
ˆ
B by using the following two con-
catenation identities [48]:
e
−i
ˆ
Aδt
e
−i
ˆ
Bδt
e
i
ˆ
Aδt
e
i
ˆ
Bδt
= e
[
ˆ
A,
ˆ
B]δt
2
+ O(δt
3
), (A1a)
e
i
ˆ
Aδt/2
e
i
ˆ
Bδt/2
e
i
ˆ
Bδt/2
e
i
ˆ
Aδt/2
= e
i(
ˆ
A+
ˆ
B)δt
+ O(δt
2
). (A1b)
If we have two second-degree generators, such as
ˆ
u =
ˆ
x
2
+
ˆ
p
2
(the generator for the rotation gate) and
ˆ
s =
ˆ
x
ˆ
p +
ˆ
p
ˆ
x (the generator for the squeezing gate), and a
third-degree (or higher) generator, such as
ˆ
x
3
(the gen-
erator for the cubic phase gate), we can easily construct
generators of all higher-degree polynomials, e.g.,
ˆ
x
4
=
−[
ˆ
x
3
, [
ˆ
x
3
,
ˆ
u]]/(18ħh
2
). This reasoning can be extended by
induction to any finite-degree polynomial in
ˆ
x and
ˆ
p (equiv-
alently, in
ˆ
a and
ˆ
a
†
) [48, 49].
In the above argument, it is important that at least one
of the generators is third-degree or higher. Indeed, com-
mutators of second-degree polynomials of
ˆ
x and
ˆ
p are also
second-degree polynomials and thus their composition us-
ing Eq. (A1) cannot generate higher-order generators. The
claim can be easily extended to N-mode systems and mul-
tivariate polynomials of the operators
ˆ
r = (
ˆ
x
1
, . . . ,
ˆ
x
N
,
ˆ
p
1
, . . . ,
ˆ
p
N
). (A2)
Combining single-mode universal gates (including at least
one of third-degree or higher) with some multimode in-
teraction, e.g., the beamsplitter interaction generated by
ˆ
b
i,j
=
ˆ
p
i
ˆ
x
j
−
ˆ
x
i
ˆ
p
j
, we can construct arbitrary-degree poly-
nomials of the quadrature operators of an N-mode system.
With the above discussion in mind, we can combine
the set of single-qumode gates generated by {
ˆ
x
i
,
ˆ
x
2
i
,
ˆ
x
3
i
,
ˆ
u
i
}
and the two-mode gate generated by
ˆ
b
i,j
for all pairs of
modes into a universal gate set. The first, third and fourth
single-mode generators correspond to the displacement, cu-
bic phase and rotation gates in Table III, while the two-
mode generator corresponds to the beamsplitter gate. Fi-
nally, the
ˆ
x
2
generator corresponds to a quadratic phase
gate,
ˆ
P(s) = exp(is
ˆ
x
2
/(2ħh)). This gate can be written in
terms of the single-mode squeezing gate and the rotation
gate as follows:
ˆ
P(s) =
ˆ
R(θ)
ˆ
S(re
iφ
), where cosh(r) =
p
1 + (s/2)
2
, tan(θ) = s/2, φ = −θ − sign(s)π/2. Equiv-
alently, we could have included the squeezing generator
ˆ
x
ˆ
p +
ˆ
p
ˆ
x in place of the quadratic phase and still had a uni-
versal set.
We have just outlined an efficient method to construct
any gate of the form exp (−iH t), where the generator H is
a polynomial of the quadrature (or mode) operators. How
can this be used for quantum computations? As shown in
Eq. (4), the eigenstates of the
ˆ
x quadrature form an (or-
thogonal) basis for representing qumode states. Thus, these
states are often taken as the computational basis of the CV
model (though other choices are also available). By apply-
ing gates as constructed above and performing measure-
ments in this computational basis (i.e., homodyne measure-
ments), we can carry out a CV quantum computation. One
primitive example [49] is to compute the product of two
numbers – whose values are stored in two qumode regis-
ters – and store the result in the state of a third qumode.
Consider the generator
ˆ
x
1
ˆ
x
2
ˆ
p
3
/ħh, which will cause the “po-
sition” operators to evolve according to
ˆ
x
1
→
ˆ
x
1
,
ˆ
x
2
→
ˆ
x
2
,
ˆ
x
3
→
ˆ
x
3
+
ˆ
x
1
ˆ
x
2
t. (A3)
A measurement in the computational basis of mode 3 will
reveal the value of the product x
1
x
2
. Note that these encod-
ings of classical continuous degress of freedom into quan-
tum registers allows for the generalization of neural net-
works into the quantum regime [66]. In the following sec-
tions we show how more complicated algorithms are con-
structed.
Multiport Gate Decompositions
One important set of gates for which it is critical to de-
rive a decomposition in terms of universal gates is the set of
multiport interferometers. A multiport interferometer, rep-
resented by a unitary operator
ˆ
U acting on N modes, will
map (in the Heisenberg picture) the annihilation operator
of mode i into a linear combination of all other modes
ˆ
a
i
→
ˆ
U
†
ˆ
a
i
ˆ
U =
X
j
(U
i j
ˆ
a
j
). (A4)
In order to preserve the commutation relations of different
modes, the matrix U must also be unitary U U
†
= U
†
U = I
N
.
Note that every annihilation operator is mapped to a linear
combination of annihilation operators and thus annihilation
and creation operators are not mixed. Because of this, mul-
tiport interferometers generate all passive (in the sense that
no photons are created or destroyed) linear optics transfor-
mations. In a pioneering work by Reck et al. [67], it was
shown that any multiport interferometer can be realized us-
ing N(N −1)/2 beamsplitter gates distributed over 2N −3
14
layers. Recently this result was improved upon by Clements
et al. [63], who showed that an equivalent decomposition
can be achieved with the same number of beamsplitters but
using only N + 1 layers.
Gaussian Operations
As mentioned in the previous subsections, generators
that are at most quadratic remain closed under composi-
tion of their associated unitaries. In the Heisenberg picture
these quadratic generators will produce all possible linear
transformations between the quadrature (or mode) opera-
tors,
ˆ
a
i
→
X
j
U
i j
ˆ
a
j
+ V
i j
ˆ
a
†
j
. (A5)
These operations are known as Gaussian operations. All
gates in Table III are Gaussian operations except for the cu-
bic phase gate. Pure Gaussian states are the set of states
that can be obtained from the (multimode) vacuum state
by Gaussian operations [50, 68]. Mixed Gaussian states
are obtained by applying Gaussian operations to thermal
states or marginalizing pure Gaussian states. Analogous to
Gaussian states, we can also define Gaussian measurements
as the set of measurements whose positive-operator valued
measure (POVM) elements can be obtained from vacuum
via Gaussian transformations. Homodyne and heterodyne
measurements are prominent examples of Gaussian mea-
surements, whereas photon counting and photodetection
are prominent examples of non-Gaussian measurements.
An important result for the CV formalism is that Gaussian
quantum computation, i.e., computation that occurs with
Gaussian states, operations and measurements, can be effi-
ciently simulated on a classical computer (this is the foun-
dation for the Gaussian backend in Strawberry Fields). This
result is the CV version [69] of the Gottesman-Knill theorem
of discrete-variable quantum information [49]. Hence we
need non-Gaussian operations in order to achieve quantum
supremacy in the CV model. Interestingly, even in the re-
stricted case where all states and gates are Gaussian, with
only the final measurements being non-Gaussian, there is
strong evidence that such a circuit cannot be efficiently sim-
ulated classically [10, 70]. More discussion and example
code for this situation (known as Gaussian boson sampling)
is provided in Appendix C.
Symplectic Formalism
In this section we review the symplectic formalism which
lies at the heart of the Gaussian backend of Strawberry
Fields. The symplectic formalism is an elegant and compact
description of Gaussian states in terms of covariance matri-
ces and mean vectors [50, 68]. To begin, the commutation
relations of the 2N position and momentum operators of
Eq. (A2) can be easily summarized as [
ˆ
r
i
,
ˆ
r
j
] = ħhiΩ
i j
where
Ω =
0 I
N
−I
N
0
is the symplectic matrix. Using the symplectic
matrix, we can define the Weyl operator D(ξ) (a multimode
displacement operator) and the characteristic function χ(ξ)
of a quantum N mode state ρ:
ˆ
D(ξ) = exp (i
ˆ
rΩξ) , χ(ξ) = 〈
ˆ
D(ξ)〉
ρ
, (A6)
where ξ ∈ R
2N
. We can now consider the Fourier transform
of the characteristic function to obtain the Wigner function
of the state
ˆ
ρ
W (r) =
Z
R
2N
d
2N
ξ
(2π)
2N
exp(−irΩξ)χ(ξ). (A7)
The 2N real arguments r of the Wigner function are the
eigenvalues of the quadrature operators from
ˆ
r.
The above recipe maps an N-mode quantum state liv-
ing in a Hilbert space to the real symplectic space K :=
(R
2N
, Ω), which is called phase space. The Wigner function
is an example of a quasiprobability distribution. Like a prob-
ability distribution over this phase space, the Wigner func-
tion is normalized to one; however, unlike a probability dis-
tribution, it may take negative values. Gaussian states have
the special property that their characteristic function (and
hence their Wigner function) is a Gaussian function of the
variables r. In this case, the Wigner function takes the form
W (r) =
exp
−
1
2
(r −
¯
r)V
−1
(r −
¯
r)
(2π)
N
p
detV
(A8)
where
¯
r = 〈
ˆ
r〉
ρ
= Tr (
ˆ
r
ˆ
ρ) is the displacement or mean vector
and V
i j
=
1
2
〈∆r
i
∆r
j
+ ∆r
i
∆r
j
〉
ρ
with ∆
ˆ
r =
ˆ
r −
¯
r. Note that
the only pure states that have non-negative Wigner func-
tions are the pure Gaussian states [58].
Each type of Gaussian state has a specific form of co-
variance matrix V and mean vector
¯
r. For the single-mode
vacuum state, we have V =
ħh
2
I
2
and
¯
r = (0, 0)
T
. A ther-
mal state (Eq. (8)) has the same (zero) displacement but
a covariance matrix V = (2
¯
n + 1)
ħh
2
I
2
, where
¯
n is the mean
photon number. A coherent state (Eq. (6)), obtained by
displacing vacuum, has the same V as the vacuum state
but a nonzero displacement vector
¯
r = 2
q
ħh
2
(Re(α), Im(α)).
Lastly, a squeezed state (Eq. (7)) has zero displacement
and covariance matrix V =
ħh
2
diag(e
−2r
, e
2r
). In the limit
r → ∞, the squeezed state’s variance in the
ˆ
x quadrature
becomes zero and the state becomes proportional to the
ˆ
x-
eigenstate |x〉 with eigenvalue 0. Consistent with the un-
certainty principle, the squeezed state’s variance in
ˆ
p blows
up. We can also consider the case r → −∞, where we find
a state proportional to the eigenstate |p〉of the
ˆ
p quadrature
with eigenvalue 0. In the laboratory and in numerical simu-
lation we must approximate every quadrature eigenstate us-
ing a finitely squeezed state (being careful that the variance
of the relevant quadrature is much smaller than any other
uncertainty relevant to the system). Any other quadrature
eigenstate can be obtained from the x = 0 eigenstate by
applying suitable displacement and rotation operators. Fi-
nally, note that Gaussian operations will transform the vec-
tor of means via affine transformations and the covariance
matrix via congruence transformations; for a detailed dis-
15
cussion of these transformations, see Sec. 2 of [50].
Given a 2N × 2N real symmetric matrix, how can we
check that it is a valid covariance matrix? If it is valid,
which operations (displacement, squeezing, multiport in-
terferometers) should be performed to prepare the corre-
sponding Gaussian state? To answer the first question: a
2N × 2N real symmetric matrix
˜
V corresponds to a Gaus-
sian quantum state if and only if
˜
V + i
ħh
2
Ω ≥ 0 (the matrix
inequality is understood in the sense that the eigenvalues
of the quantity
˜
V + i
ħh
2
Ω are nonnegative). The answer to
the second question is provided by the Bloch-Messiah reduc-
tion [60–62]. This reduction shows that any N -mode Gaus-
sian state (equivalently any covariance matrix and vector
of means) can be constructed by starting with a product of
N thermal states
N
i
ρ
i
(
¯
n
i
) (with potentially different mean
photon numbers), then applying a multiport interferometer
ˆ
U , followed by single-mode squeezing operations
N
i
S
i
(z
i
),
followed by another multiport
ˆ
V, followed by single-mode
displacement operations
N
i
D
i
(α
i
). Explicitly,
ρ
Gaussian
=
ˆ
W
O
i
ρ
i
(
¯
n
i
)
ˆ
W
†
, (A9)
ˆ
W =
O
i
D
i
(α
i
)
ˆ
V
O
i
S
i
(z
i
)
ˆ
U . (A10)
Note that if the Gaussian state is pure (which happens
if and only if det(V) = (ħh/2)
2N
), the occupation number
of the thermal states in the Bloch-Messiah decomposition
are all zero and the first interferometer will turn the vac-
uum to vacuum again. Thus for pure Gaussian states we
need only generate N single-mode squeezed states and send
them through a single multiport interferometer
ˆ
V before
displacing. For a recent discussion of this decomposition
see Ref. [ 71, 72]. More generally, the occupation numbers
of the different thermal states in Eq. (A9) n
i
= (ν
i
− 1)/2
can be obtained by calculating the symplectic eigenvalues
ν
i
of the covariance matrix V . The symplectic eigenvalues
come in pairs and are just the standard eigenvalues of the
matrix |i(2/ħh)ΩV| where the modulus is understood in the
operator sense (see Sec. II.C.1. of Ref. [50]).
16
Appendix B: Strawberry Fields Operations
In this Appendix, we present a complete list of the CV states, gates, and measurements available in Strawberry Fields.
Operation Name Definition
Vacuum() Vacuum state
The vacuum state |0〉, representing zero photons
Coherent(a) Coherent state
A displaced vacuum state, |α〉 = D(α) |0〉, α ∈ C
Squeezed(r,phi) Squeezed state
A squeezed vacuum state, |z〉 = S(z) |0〉,
where z = re
iφ
, r, φ ∈ R, r ≥ 0, φ ∈ [0, 2π)
DisplacedSqueezed(a,r,phi) Displaced squeezed state
A squeezed then displaced vacuum state,
|α, z〉 = D(α)S(z) |0〉
Thermal(n) Thermal state
ρ(
¯
n) =
P
∞
n=0
[n
n
/(1 + n)
n+1
]|n〉〈n|,
where
¯
n ∈ R
+
is the mean photon number
Fock(n)* Fock state or number state
|n〉, where n ∈ N
0
represents the photon number
Catstate(a,p)* Cat state
1
p
N
(|α〉 + exp(iπp) |−α〉), where p = 0, 1 gives an
even/odd cat state and N is the normalization
Ket(x)* Arbitrary Fock-basis ket
Prepare an arbitrary multi-mode pure state, repre-
sented by array x, in the Fock basis.
DensityMatrix(x)* Arbitrary Fock basis state
Prepare an arbitrary multi-mode mixed state, repre-
sented by a density matrix array x in the Fock basis.
Table V: State preparations available in Strawberry Fields. Those indicated with an asterisk (*) are non-Gaussian.
Operation Name Definition
Dgate(a) Displacement gate
D(α) = exp(α
ˆ
a
†
−α
∗
ˆ
a)
= exp
−i(Re(α)
ˆ
p −Im(α)
ˆ
x)
Æ
2/ħh
, α ∈ C
Xgate(x) Position displacement gate
X (x) = D(x/
p
2ħh) = exp (−ix
ˆ
p/ħh), x ∈ R
Zgate(p) Momentum displacement gate
Z(p) = D(i p/
p
2ħh) = exp (ip
ˆ
x/ħh), p ∈ R
Sgate(r,phi) Squeezing gate
S(z) = exp
(z
∗
ˆ
a
2
−z
ˆ
a
†
2
)/2
, where
z = r exp (iφ), r, φ ∈ R, r ≥ 0, φ ∈ [0, 2π)
Rgate(theta) Rotation gate
R(θ) = exp
iθ
ˆ
a
†
ˆ
a
, where θ ∈ [0, 2π)
Fouriergate() Fourier gate
F = R(π/2) = exp
i(π/2)
ˆ
a
†
ˆ
a
Pgate(s) Quadratic phase gate
P(s) = exp
is
ˆ
x
2
/(2ħh)
, where s ∈ R
Vgate(g)* Cubic phase gate
V (γ) = exp
iγ
ˆ
x
3
/(3ħh)
, where γ ∈ R
Kgate(k)* Kerr interaction gate
K(κ) = exp
iκ
ˆ
a
†
ˆ
a
2
, where κ ∈ R
Table VI: Single mode gate operations available in Strawberry Fields. Those indicated with an asterisk (*) are non-Gaussian
and can only be used with a backend that uses the Fock representation.
17
Operation Name Definition
BSgate(theta,phi) Beamsplitter
B(θ, φ) = exp
θ(e
iφ
ˆ
a
1
ˆ
a
†
2
− e
−iφ
ˆ
a
†
1
ˆ
a
2
)
,
where the transmissivity and reflectivity
amplitudes are t = cos θ, r = e
iφ
sin θ
S2gate(r,p) Two-mode squeezing gate
S
2
(z) = exp
z
∗
ˆ
a
1
ˆ
a
2
−z
ˆ
a
†
1
ˆ
a
†
2
,
where z = r exp (iφ)
CXgate(s) Controlled-X or addition gate C X (s) = exp (−is
ˆ
x
1
ˆ
p
2
/ħh), s ∈ R
CZgate(s) Controlled phase shift gate C Z(s) = exp (is
ˆ
x
1
ˆ
x
2
/ħh), s ∈ R
CKgate(k)* Controlled Kerr interaction gate
CK(κ) = exp
iκ
ˆ
a
†
1
ˆ
a
1
ˆ
a
†
2
ˆ
a
2
, κ ∈ R
Table VII: Two-mode gate operations available in Strawberry Fields.
Operation Name Definition
Gaussian(cov,mu) Gaussian state preparation
Prepares an arbitrary N-mode Gaussian state de-
fined by covariance matrix V ∈ R
2N ×2N
and means
vector µ ∈ R
2N
using the Williamson decomposition
GaussianTransform(S) Gaussian transformation
Applies an N-mode Gaussian transformation de-
fined by symplectic matrix S ∈ R
2N ×2N
using the
Bloch-Messiah decomposition
Interferometer(U) Multi-mode linear interferometer
Applies an N-mode interferometer defined by
unitary matrix U ∈ C
N×N
using the Clements
decomposition
Table VIII: Multi-mode Gaussian decompositions available in Strawberry Fields.
Operation Name Definition
MeasureHomodyne(phi) Homodyne measurement
Projects the state onto |x
φ
〉〈x
φ
|
where
ˆ
x
φ
= cos φ
ˆ
x + sin φ
ˆ
p
MeasureHeterodyne() Heterodyne measurement
Projects the state onto the coherent states, sam-
pling from the joint Husimi distribution
1
π
〈α|ρ|α〉
MeasureFock()* Photon counting Projects the state onto |n〉〈n|
Table IX: Measurement operations available in Strawberry Fields. Those indicated with an asterisk (*) are non-Gaussian and
can only be used with a backend that uses the Fock representation.
Gate Decompositions
In addition, the Strawberry Fields frontend can be used to
provide decompositions of certain compound gates. The
following gate decompositions are currently supported.
Quadratic phase gate
The quadratic phase shift gate Pgate(s) is decomposed
into a squeezing and a rotation,
P(s) = R(θ )S(re
iφ
),
where cosh(r) =
p
1 + (s/2)
2
, tan(θ ) = s/2, and φ =
−sign(s)
π
2
−θ .
Two-mode squeeze gate
The two-mode squeeze gate S2gate(z) is decomposed
into a combination of beamsplitters and single-mode
squeezers
S
2
(z) = B
†
(π/4, 0) [S(z) ⊗S(−z)] B(π/4, 0)
18
Controlled addition gate
The controlled addition or controlled-X gate CXgate(s) is
decomposed into a combination of beamsplitters and single-
mode squeezers
CX (s) = B(φ, 0) [S(r) ⊗S(−r)] B(π/2 + φ, 0),
where sin(2φ) = −1/cosh(r), cos(2φ) = −tanh(r), and
sinh(r) = −s/2.
Controlled phase gate
The controlled phase shift gate CZgate(s) is decomposed
into a controlled addition gate, with two rotation gates act-
ing on the target mode,
C Z(s) = [I ⊗R(π/2)] C X (s)
I ⊗R(π/2)
†
19
Appendix C: Quantum Algorithms
In this Appendix, we present full example code for sev-
eral important algorithms, subroutines, and use-cases for
Strawberry Fields. These examples are presented in more
detail in the online documentation located at strawberry-
fields.readthedocs.io.
Quantum Teleportation
Quantum teleportation — sometimes referred to as state
teleportation to avoid confusion with gate teleportation —
is the reliable transfer of an unknown quantum state across
spatially separated qubits or qumodes, through the use of a
classical transmission channel and quantum entanglement
[73]. Considered a fundamental quantum information pro-
tocol, it has applications ranging from quantum commu-
nication to enabling distributed information processing in
quantum computation [74].
In general, all quantum teleportation circuits work on
the same basic principle. Two distant operators, Alice and
Bob, share a maximally entangled quantum state (in dis-
crete variables, any one of the four Bell states, and in CVs,
a maximally entangled state for a fixed energy), and have
access to a classical communication channel. Alice, in pos-
session of an unknown state which she wishes to transfer
to Bob, makes a joint measurement of the unknown state
and her half of the entangled state, by projecting onto the
Bell or quadrature basis. By transmitting the results of her
measurement to Bob, Bob is then able to transform his half
of the entangled state to an accurate replica of the origi-
nal unknown state, by performing a conditional phase flip
(for qubits) or displacement (for qumodes) [75]. The CV
teleportation circuit is shown in Fig. 3.
q
0
: |ψ〉
BS
•
|x〉〈x| = m
1
q
1
: |0〉
P
BS
•
|p〉〈p| = m
2
q
2
: |0〉
P
X Z
|ψ〉
FIG. 3: State teleportation of state |ψ〉 from mode q
0
to
mode q
2
.
1 import strawberryfields as sf
2 from strawberryfields.ops import *
3 from strawberryfields.utils import scale
4 from numpy import pi, sqrt
5
6 eng,q = sf.Engine(3)
7
8 with eng:
9 psi, alice, bob = q[0], q[1], q[2]
10
11 # state to be teleported:
12 Coherent(1+0.5j) | psi
13
14 # 50-50 beamsplitter
15 BS = BSgate(pi/4, 0)
16
17 # maximally entangled states
18 Squeezed(-2) | alice
19 Squeezed(2) | bob
20 BS | (alice, bob)
21
22 # Alice performs the joint measurement
23 # in the maximally entangled basis
24 BS | (psi, alice)
25 MeasureX | psi
26 MeasureP | alice
27
28 # Bob conditionally displaces his mode
29 # based on Alice's measurement result
30 Xgate(scale(psi, sqrt(2))) | bob
31 Zgate(scale(alice, sqrt(2))) | bob
32 # end circuit
33
34 state = eng.run("gaussian")
35 # view Bob's output state and fidelity
36 print(q[0].val,q[1].val)
37 print(state.displacement([2]))
38 print(state.fidelity_coherent([0,0,1+0.5j]))
Codeblock 15: Quantum teleportation of a coherent state in
Strawberry Fields. Note: while the CV quantum teleportation
algorithm relies on infinitely squeezed resource states, in
practice a squeezing magnitude of -2 (∼ 18 dB) is sufficient.
Gate Teleportation
In the quantum state teleportation algorithm discussed in
the previous section, the quantum state is transferred from
the sender to the receiver. However, quantum teleportation
can be used in a much more powerful manner, by simulta-
neously transforming the teleported state — this is known
as gate teleportation.
In gate teleportation, rather than applying a quantum
unitary directly to the state prior to teleportation, the
unitary is applied indirectly, via the projective measure-
ment of the first subsystem onto a particular basis. This
measurement-based approach provides significant advan-
tages over applying unitary gates directly, for example by
reducing resources, and in the application of experimen-
tally hard-to-implement gates [74]. In fact, gate teleporta-
tion forms a universal quantum computing primitive, and is
a precursor to cluster state models of quantum computation
[76].
First described by Gottesman and Chuang [77] in the case
of qubits, gate teleportation was generalised for the CV case
by Bartlett and Munro [78] (see Fig. 4). In an analogous
process to the discrete-variable case, we begin with the al-
gorithm for local state teleportation. Unlike the spatially-
separated quantum state teleportation we considered in the
previous section, local teleportation can transport the state
using only two qumodes; the state we are teleporting is en-
tangled directly with the squeezed vacuum state in the mo-
mentum space through the use of a controlled phase gate.
20
FIG. 4: Gate teleportation of a quadratic phase gate P onto
input state |ψ〉.
To recover the teleported state exactly, we must perform
Weyl-Heisenberg corrections to the second mode; here, that
would be F
†
X (m)
†
, where m is the correction based on the
Homodyne measurement. However, for convenience and
simplicity, it is common to simply write the circuit without
the corrections applied explicitly.
1 import strawberryfields as sf
2 from strawberryfields.ops import *
3
4 eng, q = sf.Engine(4)
5
6 with eng:
7 # create initial states
8 Squeezed(0.1) | q[0]
9 Squeezed(-2) | q[1]
10 Squeezed(-2) | q[2]
11
12 # apply the gate to be teleported
13 Pgate(0.5) | q[1]
14
15 # conditional phase entanglement
16 CZgate(1) | (q[0], q[1])
17 CZgate(1) | (q[1], q[2])
18
19 # projective measurement onto
20 # the position quadrature
21 Fourier.H | q[0]
22 MeasureX | q[0]
23 Fourier.H | q[1]
24 MeasureX | q[1]
25 # compare against the expected output
26 # X(q1).F.P(0.5).X(q0).F.|z>
27 # not including the corrections
28 Squeezed(0.1) | q[3]
29 Fourier | q[3]
30 Xgate(q[0]) | q[3]
31 Pgate(0.5) | q[3]
32 Fourier | q[3]
33 Xgate(q[1]) | q[3]
34 # end circuit
35
36 state = eng.run("gaussian")
37 print(state.reduced_gaussian([2]))
38 print(state.reduced_gaussian([3]))
Codeblock 16: Gate teleportation of a quadratic phase gate
in Strawberry Fields.
Note that additional gates can be added to the gate tele-
portation scheme described above simply by introducing
additional qumodes with the appropriate projective mea-
q
0
: |n
1
〉
R
BS BS BS
ˆ
n
q
1
: |n
2
〉
R
BS BS
ˆ
n
q
2
: |n
3
〉
R
BS BS BS
ˆ
n
q
3
:
n
4
R
ˆ
n
FIG. 5: 4-mode boson sampling.
surements, all ‘stacked vertically’ (i.e., coupled to each con-
secutive qumode via a conditional phase gate). It is from
this primitive that the model of cluster state quantum com-
putation can be derived [76].
Boson Sampling
Introduced by Aaronson and Arkhipov [9], boson sam-
pling presented a slight deviation from the general ap-
proach in quantum computation. Rather than presenting a
theoretical model of universal quantum computation (i.e.,
a framework that enables quantum simulation of any arbi-
trary Hamiltonian [51]), boson sampling-based devices are
instead an example of an intermediate quantum computer,
designed to experimentally implement a computation that
is intractable classically [79–82].
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 gen-
erated by sampling from the probability of single photon
measurements from the output of the linear interferometer.
For example, consider N single photon Fock states, |ψ〉 =
|m
1
, m
2
, . . . , m
N
〉, composed of b =
P
i
m
i
photons, incident
on an N-mode linear interferometer described by the uni-
tary U acting on the input mode creation and annihilation
operators. It was shown that the probability of detecting n
j
photons at the jth output mode is given by [9]
Pr(n
1
, . . . , n
N
) =
|Per(U
st
)|
2
m
1
! ···m
N
!n
1
! ···n
N
!
; (C1)
i.e., the sampled single photon probability distribution is
proportional to the permanent of U
st
, a submatrix of the
interferometer unitary, dependent upon the input and out-
put Fock states. Whilst the determinant can be calculated
efficiently on classical computers, calculation of the per-
manent belongs to the computational complexity class #P-
Hard problems [83], which are strongly believed to be clas-
sically hard to calculate. This implies that simulating bo-
son sampling is an intractable task for classical comput-
ers, providing an avenue for the demonstration of quantum
supremacy.
Continuous-variable quantum computation is ideally
suited to the simulation and demonstration of the boson
sampling scheme, due to its grounding in quantum optics.
In quantum linear optics, the multimode linear interferom-
eter is commonly decomposed into two-mode beamsplitters
(BSgate) and single-mode phase shifters (Rgate) [67], al-
lowing for a straightforward translation into a CV quantum
21
q
0
: |0〉
S
R
BS BS BS
ˆ
n
q
1
: |0〉
S
R
BS BS
ˆ
n
q
2
: |0〉
S
R
BS BS BS
ˆ
n
q
3
: |0〉
S
R
ˆ
n
FIG. 6: 4-mode Gaussian boson sampling.
circuit. In order to allow for arbitrary linear unitaries on m
imput modes, we must have a minimum of m + 1 columns
in the beamsplitter array [63].
1 import numpy as np
2 import strawberryfields as sf
3 from strawberryfields.ops import *
4
5 # initialise the engine and register
6 eng, q = sf.Engine(4)
7
8 with eng:
9 # prepare the input fock states
10 Fock(1) | q[0]
11 Fock(1) | q[1]
12 Vac | q[2]
13 Fock(1) | q[3]
14
15 # rotation gates
16 Rgate(0.5719)
17 Rgate(-1.9782)
18 Rgate(2.0603)
19 Rgate(0.0644)
20
21 # beamsplitter array
22 BSgate(0.7804, 0.8578) | (q[0], q[1])
23 BSgate(0.06406, 0.5165) | (q[2], q[3])
24 BSgate(0.473, 0.1176) | (q[1], q[2])
25 BSgate(0.563, 0.1517) | (q[0], q[1])
26 BSgate(0.1323, 0.9946) | (q[2], q[3])
27 BSgate(0.311, 0.3231) | (q[1], q[2])
28 BSgate(0.4348, 0.0798) | (q[0], q[1])
29 BSgate(0.4368, 0.6157) | (q[2], q[3])
30 # end circuit
31
32 # run the engine
33 state = eng.run("fock", cutoff_dim=7)
34
35 # extract the joint Fock probabilities
36 probs = state.all_fock_probs()
37
38 # print the joint Fock state probabilities
39 print(probs[1,1,0,1])
40 print(probs[2,0,0,1])
Codeblock 17: 4-mode boson sampling example in
Strawberry Fields. Parameters are chosen arbitrarily.
Gaussian Boson Sampling
While boson sampling allows the experimental imple-
mentation of a sampling problem that is hard to simulate
classically, one of the setbacks in experimental setups is
scalability, due to its dependence on an array of simultane-
ously emitting single photon sources. Currently, most phys-
ical implementations of boson sampling make use of a pro-
cess known as Spontaneous Parametric Down-Conversion
(SPDC) to generate the single-photon source inputs. How-
ever, 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 in-
creases exponentially.
In order to simulate a deterministic single photon source
array, several variations on boson sampling have been pro-
posed; the most well-known being scattershot boson sam-
pling [70]. However, a recent boson sampling variation by
[10] negates the need for single photon Fock states alto-
gether, by showing that incident Gaussian states – in this
case, single mode squeezed states – can produce problems
in the same computational complexity class as boson sam-
pling. Even more significantly, this mitigates the scalabil-
ity problem with single photon sources, as single mode
squeezed states can be simultaneously generated experi-
mentally.
With an input ensemble of N single-mode squeezed states
with squeezing parameter z = r ∈ R, incident on a linear-
interferometer described by unitary U, it can be shown
that the probability of detecting an output photon pattern
(n
1
, . . . , n
N
), where n
k
∈ {0,1}, is given by [10]
Pr(n
1
, . . . , n
N
) =
Haf[(U
L
i
tanh(r
i
)U
T
)
S
]
2
Q
i
cosh(r
i
)
, (C2)
where S denotes the subset of modes where a photon was
detected and Haf[·] is the Hafnian [84, 85]. That is, the
sampled single photon probability distribution is propor-
tional to the hafnian of a submatrix of U
L
i
tanh(r
i
)U
T
.
The hafnian is known to be a generalization of the per-
manent, and can be used to count the number of perfect
matchings of an arbitrary graph [86]. The formula above
can be generalized to pure Gaussian states with finite dis-
placements by using the loop hafnian function which counts
the number of perfect matchings of graphs with loops[56].
Since any algorithm that could calculate the hafnian could
also calculate the permanent, it follows that calculating the
hafnian remains a classically hard problem; indeed, the best
known classical algorithm for the calculation of a hafnian of
an arbitrary symmetric complex matrix of size N scales like
O(N
3
2
N/2
) [87]. The hardness of approximate GBS, under
imperfections such as loss, is a subject of current research
[88, 89].
1 import numpy as np
2 import strawberryfields as sf
3 from strawberryfields.ops import *
4
5 # initialise the engine and register
6 eng, q = sf.Engine(4)
7
8 with eng:
9 # prepare the input squeezed states
22
10 S = Sgate(1)
11 S | q[0]
12 S | q[1]
13 S | q[2]
14 S | q[3]
15
16 # rotation gates
17 Rgate(0.5719) | q[0]
18 Rgate(-1.9782) | q[1]
19 Rgate(2.0603) | q[2]
20 Rgate(0.0644) | q[3]
21
22 # beamsplitter array
23 BSgate(0.7804, 0.8578) | (q[0], q[1])
24 BSgate(0.06406, 0.5165) | (q[2], q[3])
25 BSgate(0.473, 0.1176) | (q[1], q[2])
26 BSgate(0.563, 0.1517) | (q[0], q[1])
27 BSgate(0.1323, 0.9946) | (q[2], q[3])
28 BSgate(0.311, 0.3231) | (q[1], q[2])
29 BSgate(0.4348, 0.0798) | (q[0], q[1])
30 BSgate(0.4368, 0.6157) | (q[2], q[3])
31 # end circuit
32
33 # run the engine
34 state = eng.run("gaussian")
35
36 # Fock states to measure at output
37 measure_states = [[0,0,0,0], [1,1,0,0],
[0,1,0,1], [1,1,1,1], [2,0,0,0]],→
38
39 # extract the probabilities of calculating
40 # several different Fock states at the output
41 # and print them to the terminal
42 for i in measure_states:
43 prob = state.fock_prob(i)
44 print("|{}>: {}".format("".join(str(j) for
j in i), prob)),→
Codeblock 18: 4-mode Gaussian boson sampling example
in Strawberry Fields. Parameters are chosen arbitrarily.
Instantaneous Quantum Polynomial (IQP)
Like boson sampling and Gaussian boson sampling before
it, the instantaneous quantum polynomial (IQP) protocol is
a restricted, non-universal model of quantum computation,
designed to implement a computation that is classically in-
tractable and verifiable by a remote adjudicator. First intro-
duced by Shepherd and Bremner [90], IQP circuits are de-
fined by the following conditions: (i) there are N inputs in
state |0〉acted on by Hadamard gates, (ii) the resulting com-
putation is diagonal in the computational basis by randomly
selecting from the set U = {R(π/4),
p
C Z} (hence the term
‘instantaneous’; since all gates commute, they can be ap-
plied in any temporal order), and (iii) the output qubits are
measured in the Pauli-X basis. Efficient classical sampling
of the resulting probability distribution H
⊗N
UH
⊗N
|0〉
⊗N
—
even approximately [13] or in the presence of noise [91] —
has been shown to be #P-hard, and would result in the col-
lapse of the polynomial hierarchy to the third level [12, 92].
Unlike boson sampling and Gaussian boson sampling,
|0〉
S
•
Z
•
V
ˆ
p
|0〉
S
• •
V
ˆ
p
|0〉
S
V
Z
• •
Z
ˆ
p
|0〉
S
Z Z Z
ˆ
p
FIG. 7: 4-mode instantaneous quantum polynomial (IQP)
CV circuit example.
however, the IQP protocol was not constructed with quan-
tum linear optics in mind. Nevertheless, the IQP model
was recently extended to the CV formulation of quantum
computation by Douce et al. [17], taking advantage of
the ability to efficiently prepare large resource states, and
the higher efficiencies afforded by homodyne detection.
Furthermore, the computational complexity results of the
discrete-variable case apply equally to the so-called CV-IQP
model, assuming a specific input squeezing parameter de-
pendent on the circuit size. The CV-IQP model is defined as
follows:
1. inputs are squeezed momentum states |z〉, where z =
r ∈ R and r < 0;
2. the unitary transformation is diagonal in the
ˆ
x
quadrature basis, by randomly selecting from the set
of gates U = {Z(p), C Z(s), V (γ)};
3. and finally, homodyne measurements are performed
on all modes in the
ˆ
p quadrature.
1 import strawberryfields as sf
2 from strawberryfields.ops import *
3
4 # initialise the engine and register
5 eng, q = sf.Engine(4)
6
7 with eng:
8 # prepare the input squeezed states
9 S = Sgate(-1)
10 S | q[0]
11 S | q[1]
12 S | q[2]
13 S | q[3]
14
15 # gates diagonal in the x quadrature
16 CZgate(0.5719) | (q[0], q[1])
17 Vgate(-1.9782) | q[2]
18 Zgate(2.0603) | q[3]
19 Zgate(-0.7804) | q[0]
20 Zgate(0.8578) | q[2]
21 CZgate(1.321) | (q[1], q[2])
22 Zgate(0.473) | q[3]
23 CZgate(0.9946) | (q[0], q[2])
24 Zgate(0.1223) | q[3]
25 Vgate(0.6157) | q[0]
26 Vgate(0.3110) | q[1]
27 Zgate(-1.243) | q[2]
23
28 # end circuit
29
30 # run the engine
31 eng.run("fock", cutoff_dim=5)
Codeblock 19: 4-mode instantaneous quantum polynomial
(IQP) example in Strawberry Fields.
Moreover, the resulting probability distributions have
been shown to be given by integrals of oscillating functions
in large dimensions, which are believed to be intractable
to compute by classical computers. This leads to the result
that even in the case of finite squeezing and reduced mea-
surement precision, approximate sampling from the output
CV-IQP model remains classically intractable [17, 93].
Hamiltonian Simulation
The simulation of atoms, molecules and other biochemi-
cal systems is another application uniquely suited to quan-
tum computation. For example, the ground state energy
of large systems, the dynamical behaviour of an ensemble
of molecules, or complex molecular behaviour such as pro-
tein folding, are often computationally hard or downright
impossible to determine via classical computation or exper-
imentation [94, 95].
In the discrete-variable qubit model, efficient methods of
Hamiltonian simulation have been discussed at length, pro-
viding several implementations depending on properties of
the Hamiltonian, and resulting in a linear simulation time
[96, 97]. Efficient implementations of Hamiltonian simula-
tion also exist in the CV formulation [98], with specific ap-
plication to Bose-Hubbard Hamiltonians (describing a sys-
tem of interacting bosonic particles on a lattice of orthog-
onal position states [99]). As such, this method is ideally
suited to photonic quantum computation.
q
0
: |2〉
BS
K R
BS
K R
q
1
: |0〉
K R K R
FIG. 8: Bose-Hubbard Hamiltonian simulation for a 2-node
lattice, with k = 2. The error is of the order O(t
2
/k).
Considering a lattice composed of two adjacent nodes
characterised by Adjacency matrix A, the Bose-Hubbard
Hamiltonian is given by
H = J
X
i
X
j
A
i j
ˆ
a
†
i
ˆ
a
j
+
1
2
U
X
i
ˆ
n
i
(
ˆ
n
i
−1)
= J (
ˆ
a
†
1
ˆ
a
2
+
ˆ
a
†
2
ˆ
a
1
) +
1
2
U(
ˆ
n
2
1
−
ˆ
n
1
+
ˆ
n
2
2
−
ˆ
n
2
), (C3)
where J represents the transition of the boson between
nodes, and U is the on-site interaction potential. Applying
the Lie product formula, we find that
e
iH t
=
exp
−i
J t
k
(
ˆ
a
†
1
ˆ
a
2
+
ˆ
a
†
2
ˆ
a
1
)
exp
−i
U t
2k
ˆ
n
2
1
exp
−i
U t
2k
ˆ
n
2
2
exp
i
U t
2k
ˆ
n
1
exp
i
U t
2k
ˆ
n
2
k
+ O
t
2
/k
, (C4)
where O
t
2
/k
is the order of the error term, derived from
the Lie product formula. Comparing this to the form of var-
ious gates in the CV circuit model, we can write this as the
product of beamsplitters, Kerr gates, and rotation gates:
e
iH t
= {BS (θ , φ) [K(r)R(−r) ⊗ K(r)R(−r)]}
k
+ O
t
2
/k
(C5)
where θ = −J t/k, φ = π/2, and r = −U t/2k. Using J = 1,
U = 1.5, k = 20, and a timestep of t = 1.086, this can be
easily implemented in Strawberry Fields using only 20×3 =
60 gates.
1 import strawberryfields as sf
2 from strawberryfields.ops import *
3 from numpy import pi
4
5 # initialise the engine and register
6 eng, q = sf.Engine(2)
7
8 # set the Hamiltonian parameters
9 J = 1 # hopping transition
10 U = 1.5 # on-site interaction
11 k = 20 # Lie product decomposition
terms,→
12 t = 1.086 # timestep
13 theta = -J*t/k
14 r = -U*t/(2*k)
15
16 with eng:
17 # prepare the initial state
18 Fock(2) | q[0]
19
20 # Two node tight-binding
21 # Hamiltonian simulation
22
23 for i in range(k):
24 BSgate(theta, pi/2) | (q[0], q[1])
25 Kgate(r) | q[0]
26 Rgate(-r) | q[0]
27 Kgate(r) | q[1]
28 Rgate(-r) | q[1]
29 # end circuit
30
31 # run the engine
32 state = eng.run("fock", cutoff_dim=5)
33 # the output state probabilities
34 print(state.fock_prob([0,2]))
35 print(state.fock_prob([1,1]))
36 print(state.fock_prob([2,0]))
Codeblock 20: Tight-binding Hamiltonian simulation for a
2-node lattice in Strawberry Fields.
24
For more complex Hamiltonian CV decompositions, in-
cluding those with nearest-neighbour, see Kalajdzievski
et al. [98]. This decomposition is also implemented in
the SFOpenBoson plugin [45], which provides an interface
between OpenFermion, the quantum electronic structure
package, and Strawberry Fields. This allows arbitrary Bose-
Hubbard Hamiltonians, generated in OpenFermion, to be
simulated using Strawberry Fields.
Optimization of Quantum Circuits
One of the unique features of Strawberry Fields is that it
has been designed from the start to support modern com-
putational methods like automatic gradients, optimization,
and machine learning by leveraging a TensorFlow backend.
We have already given an overview in the main text of how
these features are accessed in Strawberry Fields. We present
here a complete example for the optimization of a quantum
circuit. Our goal in this circuit is to find the Dgate param-
eter which leads to the highest probability of a n = 1 Fock
state output. This simple baseline example can be straight-
forwardly extended to much more complex circuits. As op-
timization is a key ingredient of machine learning, this ex-
ample can also serve as a springboard for more advanced
data-driven modelling tasks.
We note that optimization here is in a variational sense,
i.e., we choose an ansatz for the optimization by fixing the
discrete structure of our circuits (the selection and arrange-
ment of specific states, gates, and measurements). The op-
timization then proceeds over the parameters of these op-
erations. Finally, we emphasize that for certain circuits and
objective functions, the optimization might be non-convex
in general. Thus we should not assume (without proof) that
the solution obtained via a Strawberry Fields optimization
is always the global optimum. However, it is often the case
in machine learning that local optima can still provide ef-
fective solutions, so this may not be an issue, depending on
the application.
1 import strawberryfields as sf
2 from strawberryfields.ops import *
3 import tensorflow as tf
4
5 eng, q = sf.Engine(1)
6
7 alpha = tf.Variable(0.1)
8 with eng:
9 Dgate(alpha) | q[0]
10 state = eng.run("tf", cutoff_dim=7,
eval=False),→
11
12 # loss is probability for the Fock state n=1
13 prob = state.fock_prob([1])
14 loss = -prob # negative sign to maximize prob
15
16 # Set up optimization
17 optimizer = tf.train.GradientDescentOptimizer(
learning_rate=0.1),→
18 minimize_op = optimizer.minimize(loss)
19
20 # Create TF Session and initialize variables
21 sess = tf.Session()
22 sess.run(tf.global_variables_initializer())
23
24 # Carry out optimization
25 for step in range(50):
26 prob_val, _ = sess.run([prob, minimize_op])
27 print("Value at step {}: {}".format(step,
prob_val)),→
Codeblock 21: Example of optimizing quantum circuit
parameters in Strawberry Fields. This can be extended to a
machine learning setting by incorporating data and a more
sophisticated loss function.
[1] P. W. Shor. Polynomial-time algorithms for prime fac-
torization and discrete logarithms on a quantum com-
puter. SIAM review, 41(2):303–332, 1999. doi:
10.1137/S0036144598347011.
[2] L. K. Grover. A fast quantum mechanical algorithm for
database search. In Proceedings of the twenty-eighth annual
ACM Symposium on Theory of Computing, pages 212–219.
ACM, 1996.
[3] A. W. Harrow, A. Hassidim, and S. Lloyd. Quantum algorithm
for linear systems of equations. Physical Review Letters, 103
(15):150502, 2009. doi:10.1103/PhysRevLett.103.150502.
[4] J. Preskill. Quantum Computing in the NISQ era and beyond.
Quantum, 2:79, 2018. doi:10.22331/q-2018-08-06-79.
[5] A. Peruzzo, J. McClean, P. Shadbolt, M.-H. Yung, X.-Q. Zhou,
P. J. Love, A. Aspuru-Guzik, and J. L. O’Brien. A variational
eigenvalue solver on a photonic quantum processor. Nature
Communications, 5, 2014. doi:10.1038/ncomms5213.
[6] J. R. McClean, J. Romero, R. Babbush, and A. Aspuru-
Guzik. The theory of variational hybrid quantum-classical
algorithms. New Journal of Physics, 18(2):023023, 2016.
doi:10.1088/1367-2630/18/2/023023.
[7] E. Farhi, J. Goldstone, and S. Gutmann. A quan-
tum approximate optimization algorithm. arXiv preprint
arXiv:1411.4028, 2014.
[8] E. Farhi and A. W. Harrow. Quantum supremacy through
the quantum approximate optimization algorithm. arXiv
preprint arXiv:1602.07674, 2016.
[9] S. Aaronson and A. Arkhipov. The computational complexity
of linear optics. In Proceedings of the forty-third annual ACM
symposium on Theory of computing, pages 333–342. ACM,
2011. doi:10.1145/1993636.1993682.
[10] C. S. Hamilton, R. Kruse, L. Sansoni, S. Barkhofen,
C. Silberhorn, and I. Jex. Gaussian boson sam-
pling. Physical Review Letters, 119:170501, 2017. doi:
10.1103/PhysRevLett.119.170501.
[11] L. Chakhmakhchyan, R. Garcia-Patron, and N. J. Cerf. Boson
sampling with Gaussian measurements. Physical Review A,
96:032326, 2017. doi:10.1103/PhysRevA.96.032326.
25
[12] M. J. Bremner, R. Jozsa, and D. J. Shepherd. Classical simu-
lation of commuting quantum computations implies collapse
of the polynomial hierarchy. In Proceedings of the Royal Soci-
ety of London A: Mathematical, Physical and Engineering Sci-
ences, page rspa20100301. The Royal Society, 2010. doi:
10.1098/rspa.2010.0301.
[13] M. J. Bremner, A. Montanaro, and D. J. Shepherd. Average-
case complexity versus approximate simulation of commut-
ing quantum computations. Physical Review Letters, 117(8):
080501, 2016. doi:10.1103/PhysRevLett.117.080501.
[14] S. Boixo, S. V. Isakov, V. N. Smelyanskiy, R. Babbush, N. Ding,
Z. Jiang, M. J. Bremner, J. M. Martinis, and H. Neven. Char-
acterizing quantum supremacy in near-term devices. Nature
Physics, 14(6):595, 2018. doi:10.1038/s41567-018-0124-x.
[15] S. Aaronson and L. Chen. Complexity-Theoretic Foundations
of Quantum Supremacy Experiments. In R. O’Donnell, edi-
tor, 32nd Computational Complexity Conference (CCC 2017),
volume 79 of Leibniz International Proceedings in Informat-
ics (LIPIcs), pages 22:1–22:67. Schloss Dagstuhl–Leibniz-
Zentrum fuer Informatik, 2017. ISBN 978-3-95977-040-8.
doi:10.4230/LIPIcs.CCC.2017.22.
[16] C. Neill, P. Roushan, K. Kechedzhi, S. Boixo, S. V. Isakov,
V. Smelyanskiy, A. Megrant, B. Chiaro, A. Dunsworth,
K. Arya, et al. A blueprint for demonstrating quan-
tum supremacy with superconducting qubits. Science, 360
(6385):195–199, 2018. doi:10.1126/science.aao4309.
[17] T. Douce, D. Markham, E. Kashefi, E. Diamanti, T. Coudreau,
P. Milman, P. van Loock, and G. Ferrini. Continuous-
variable instantaneous quantum computing is hard to
sample. Physical Review Letters, 118(7), 2017. doi:
10.1103/PhysRevLett.118.070503.
[18] A. Finnila, M. Gomez, C. Sebenik, C. Stenson, and J. Doll.
Quantum annealing: a new method for minimizing multi-
dimensional functions. Chemical Physics Letters, 219(5-6):
343–348, 1994. doi:10.1016/0009-2614(94)00117-0.
[19] M. W. Johnson, M. H. Amin, S. Gildert, T. Lanting, F. Hamze,
N. Dickson, R. Harris, A. J. Berkley, J. Johansson, P. Bunyk,
et al. Quantum annealing with manufactured spins. Nature,
473(7346):194–198, 2011. doi:10.1038/nature10012.
[20] M.-H. Yung, J. Casanova, A. Mezzacapo, J. McClean,
L. Lamata, A. Aspuru-Guzik, and E. Solano. From transistor
to trapped-ion computers for quantum chemistry. Scientific
Reports, 4:3589, 2014. doi:10.1038/srep03589.
[21] P. O’Malley, R. Babbush, I. Kivlichan, J. Romero, J. Mc-
Clean, R. Barends, J. Kelly, P. Roushan, A. Tranter, N. Ding,
et al. Scalable quantum simulation of molecular en-
ergies. Physical Review X, 6(3):031007, 2016. doi:
10.1103/PhysRevX.6.031007.
[22] Y. Shen, X. Zhang, S. Zhang, J.-N. Zhang, M.-H. Yung,
and K. Kim. Quantum implementation of the unitary
coupled cluster for simulating molecular electronic struc-
ture. Physical Review A, 95(2):020501, 2017. doi:
10.1103/PhysRevA.95.020501.
[23] A. Kandala, A. Mezzacapo, K. Temme, M. Takita, M. Brink,
J. M. Chow, and J. M. Gambetta. Hardware-efficient varia-
tional quantum eigensolver for small molecules and quan-
tum magnets. Nature, 549(7671):242–246, 2017. doi:
10.1038/nature23879.
[24] G. Rosenberg, P. Haghnegahdar, P. Goddard, P. Carr, K. Wu,
and M. L. de Prado. Solving the optimal trading trajectory
problem using a quantum annealer. IEEE Journal of Selected
Topics in Signal Processing, 10(6):1053–1060, 2016. doi:
10.1109/JSTSP.2016.2574703.
[25] A. Lucas. Ising formulations of many NP problems. Frontiers
in Physics, 2:5, 2014. doi:10.3389/fphy.2014.00005.
[26] F. Neukart, D. Von Dollen, G. Compostella, C. Seidel,
S. Yarkoni, and B. Parney. Traffic flow optimization using
a quantum annealer. Frontiers in ICT, 4:29, 2017. doi:
10.3389/fict.2017.00029.
[27] H. Neven, V. S. Denchev, G. Rose, and W. G. Macready. Train-
ing a large scale classifier with the quantum adiabatic algo-
rithm. arXiv preprint arXiv:0912.0779, 2009.
[28] K. L. Pudenz and D. A. Lidar. Quantum adiabatic machine
learning. Quantum Information Processing, 12(5):2027–
2070, 2013. doi:10.1007/s11128-012-0506-4.
[29] D. Crawford, A. Levit, N. Ghadermarzy, J. S. Oberoi, and
P. Ronagh. Reinforcement learning using quantum Boltz-
mann machines. Quantum Information & Computation, 18
(1-2):0051–0074, 2018. doi:10.26421/QIC18.1-2.
[30] M. H. Amin, E. Andriyash, J. Rolfe, B. Kulchytskyy, and
R. Melko. Quantum boltzmann machine. Phys. Rev. X, 8:
021050, May 2018. doi:10.1103/PhysRevX.8.021050.
[31] D. Ristè, M. P. Da Silva, C. A. Ryan, A. W. Cross, A. D. Cór-
coles, J. A. Smolin, J. M. Gambetta, J. M. Chow, and B. R.
Johnson. Demonstration of quantum advantage in machine
learning. npj Quantum Information, 3(1):16, 2017. doi:
10.1038/s41534-017-0017-3.
[32] G. Verdon, M. Broughton, and J. Biamonte. A quantum al-
gorithm to train neural networks using low-depth circuits.
arXiv preprint arXiv:1712.05304, 2017.
[33] J. Otterbach, R. Manenti, N. Alidoust, A. Bestwick, M. Block,
B. Bloom, S. Caldwell, N. Didier, E. S. Fried, S. Hong, et al.
Unsupervised machine learning on a hybrid quantum com-
puter. arXiv preprint arXiv:1712.05771, 2017.
[34] M. Schuld and N. Killoran. Quantum machine learning in
feature hilbert spaces. Physical Review Letters, 122:040504,
Feb 2019. doi:10.1103/PhysRevLett.122.040504.
[35] A. S. Green, P. L. Lumsdaine, N. J. Ross, P. Selinger, and B. Val-
iron. Quipper: a scalable quantum programming language.
In ACM SIGPLAN Notices, volume 48, pages 333–342. ACM,
2013. doi:10.1145/2491956.2462177.
[36] D. Wecker and K. M. Svore. Liqui|>: A software design ar-
chitecture and domain-specific language for quantum com-
puting. arXiv preprint arXiv:1402.4467, 2014.
[37] A. JavadiAbhari, S. Patil, D. Kudrow, J. Heckey, A. Lvov, F. T.
Chong, and M. Martonosi. ScaffCC: Scalable compilation
and analysis of quantum programs. Parallel Computing, 45:
2–17, 2015. doi:10.1016/j.parco.2014.12.001.
[38] M. Smelyanskiy, N. P. Sawaya, and A. Aspuru-Guzik. qHiP-
STER: the quantum high performance software testing envi-
ronment. arXiv preprint arXiv:1601.07195, 2016.
[39] S. Pakin. A quantum macro assembler. In High Performance
Extreme Computing Conference (HPEC), 2016 IEEE, pages 1–
8. IEEE, 2016. doi:10.1109/HPEC.2016.7761637.
[40] R. S. Smith, M. J. Curtis, and W. J. Zeng. A practi-
cal quantum instruction set architecture. arXiv preprint
arXiv:1608.03355, 2016.
[41] D. S. Steiger, T. Häner, and M. Troyer. ProjectQ: an open
source software framework for quantum computing. Quan-
tum, 2:49, 2018. doi:10.22331/q-2018-01-31-49.
[42] A. W. Cross, L. S. Bishop, J. A. Smolin, and J. M. Gam-
betta. Open quantum assembly language. arXiv preprint
arXiv:1707.03429, 2017.
[43] E. S. Fried, N. P. Sawaya, Y. Cao, I. D. Kivlichan, J. Romero,
and A. Aspuru-Guzik. qtorch: The quantum tensor con-
traction handler. PloS one, 13(12):e0208510, 2018. doi:
10.1371/journal.pone.0208510.
[44] A. J. McCaskey, E. F. Dumitrescu, D. Liakh, M. Chen, W.-c.
26
Feng, and T. S. Humble. Extreme-scale programming model
for quantum acceleration within high performance comput-
ing. arXiv preprint arXiv:1710.01794, 2017.
[45] J. R. McClean, I. D. Kivlichan, D. S. Steiger, Y. Cao, E. S.
Fried, C. Gidney, T. Häner, V. Havlí
ˇ
cek, Z. Jiang, M. Neeley,
et al. OpenFermion: The electronic structure package for
quantum computers. arXiv preprint arXiv:1710.07629, 2017.
[46] S. Liu, X. Wang, L. Zhou, J. Guan, Y. Li, Y. He, R. Duan,
and M. Ying. Q|SI〉: A quantum programming environment.
In Symposium on Real-Time and Hybrid Systems, pages 133–
164. Springer, 2018. doi:10.1007/978-3-030-01461-2_8.
[47] K. Svore, A. Geller, M. Troyer, J. Azariah, C. Granade,
B. Heim, V. Kliuchnikov, M. Mykhailova, A. Paz, and M. Roet-
teler. Q#: Enabling scalable quantum computing and de-
velopment with a high-level DSL. In Proceedings of the Real
World Domain Specific Languages Workshop 2018, page 7.
ACM, 2018. doi:10.1145/3183895.3183901.
[48] S. Lloyd and S. L. Braunstein. Quantum computation over
continuous variables. Physical Review Letters, 82(8):1784,
1999. doi:10.1103/PhysRevLett.82.1784.
[49] S. L. Braunstein and P. van Loock. Quantum information with
continuous variables. Reviews of Modern Physics, 77:513–
577, 2005. doi:10.1103/RevModPhys.77.513.
[50] C. Weedbrook, S. Pirandola, R. García-Patrón, N. J. Cerf, T. C.
Ralph, J. H. Shapiro, and S. Lloyd. Gaussian quantum infor-
mation. Reviews of Modern Physics, 84(2):621, 2012. doi:
10.1103/RevModPhys.84.621.
[51] M. A. Nielsen and I. Chuang. Quantum computation and
quantum information. Cambridge University Press, 2002.
[52] A. Graves, G. Wayne, and I. Danihelka. Neural Turing ma-
chines. arXiv preprint arXiv:1410.5401, 2014.
[53] A. Graves, G. Wayne, M. Reynolds, T. Harley, I. Danihelka,
A. Grabska-Barwi
´
nska, S. G. Colmenarejo, E. Grefenstette,
T. Ramalho, J. Agapiou, et al. Hybrid computing using a
neural network with dynamic external memory. Nature, 538
(7626):471–476, 2016. doi:10.1038/nature20101.
[54] S. van der Walt, S. C. Colbert, and G. Varoquaux. The NumPy
array: a structure for efficient numerical computation. Com-
puting in Science & Engineering, 13(2):22–30, 2011. doi:
10.1109/MCSE.2011.37.
[55] M. Abadi, A. Agarwal, P. Barham, E. Brevdo, Z. Chen,
C. Citro, G. S. Corrado, A. Davis, J. Dean, M. Devin, et al.
Tensorflow: Large-scale machine learning on heterogeneous
distributed systems. arXiv preprint arXiv:1603.04467, 2016.
[56] N. Quesada. A faster calculation of franck-condon factors
and fock matrix elements of gaussian unitaries using loop
hafnians. arXiv preprint arXiv:1811.09597, 2018.
[57] D. Gottesman, A. Kitaev, and J. Preskill. Encoding a qubit in
an oscillator. Physical Review A, 64(1):012310, 2001. doi:
10.1103/PhysRevA.64.012310.
[58] A. Ferraro, S. Olivares, and M. G. Paris. Gaussian states
in continuous variable quantum information. arXiv preprint
quant-ph/0503237, 2005.
[59] J. Williamson. On the algebraic problem concerning the nor-
mal forms of linear dynamical systems. American Journal of
Mathematics, 58(1):141, 1936. doi:10.2307/2371062.
[60] C. Bloch and A. Messiah. The canonical form of an anti-
symmetric tensor and its application to the theory of su-
perconductivity. Nuclear Physics, 39:95–106, 1962. doi:
10.1016/0029-5582(62)90377-2.
[61] S. L. Braunstein. Squeezing as an irreducible re-
source. Physical Review A, 71(5):055801, 2005. doi:
10.1103/PhysRevA.71.055801.
[62] R. Simon, N. Mukunda, and B. Dutta. Quantum-noise ma-
trix for multimode systems: U(n) invariance, squeezing, and
normal forms. Physical Review A, 49(3):1567, 1994. doi:
10.1103/PhysRevA.49.1567.
[63] W. R. Clements, P. C. Humphreys, B. J. Metcalf, W. S.
Kolthammer, and I. A. Walsmley. Optimal design for uni-
versal multiport interferometers. Optica, 3(12):1460–1465,
2016. doi:10.1364/OPTICA.3.001460.
[64] J. M. Arrazola, T. R. Bromley, J. Izaac, C. R. Myers, K. Bradler,
and N. Killoran. Machine learning method for state prepa-
ration and gate synthesis on photonic quantum computers.
Quantum Science and Technology, 2018. doi:10.1088/2058-
9565/aaf59e.
[65] N. Quesada and A. M. Bra
´
nczyk. Gaussian functions
are optimal for waveguided nonlinear-quantum-optical pro-
cesses. Physical Review A, 98:043813, 2018. doi:
10.1103/PhysRevA.98.043813.
[66] N. Killoran, T. R. Bromley, J. M. Arrazola, M. Schuld, N. Que-
sada, and S. Lloyd. Continuous-variable quantum neural net-
works. arXiv preprint arXiv:1806.06871, 2018.
[67] M. Reck, A. Zeilinger, H. J. Bernstein, and P. Bertani.
Experimental realization of any discrete unitary opera-
tor. Physical Review Letters, 73:58–61, 1994. doi:
10.1103/PhysRevLett.73.58.
[68] A. Serafini. Quantum Continuous Variables: A Primer of The-
oretical Methods. CRC Press, 2017.
[69] S. D. Bartlett, B. C. Sanders, S. L. Braunstein, and K. Nemoto.
Efficient classical simulation of continuous variable quantum
information processes. Physical Review Letters, 88:097904,
2002. doi:10.1103/PhysRevLett.88.097904.
[70] A. Lund, A. Laing, S. Rahimi-Keshari, T. Rudolph, J. L.
O’Brien, and T. Ralph. Boson sampling from a Gaussian
state. Physical Review Letters, 113(10):100502, 2014. doi:
10.1103/PhysRevLett.113.100502.
[71] G. Cariolaro and G. Pierobon. Bloch-Messiah reduction of
Gaussian unitaries by Takagi factorization. Physical Review
A, 94(6):062109, 2016. doi:10.1103/PhysRevA.94.062109.
[72] G. Cariolaro and G. Pierobon. Reexamination of Bloch-
Messiah reduction. Physical Review A, 93(6):062115, 2016.
doi:10.1103/PhysRevA.93.062115.
[73] C. H. Bennett, G. Brassard, C. Crépeau, R. Jozsa, A. Peres,
and W. K. Wootters. Teleporting an unknown quantum
state via dual classical and Einstein-Podolsky-Rosen chan-
nels. Physical Review Letters, 70:1895–1899, 1993. doi:
10.1103/PhysRevLett.70.1895.
[74] A. Furusawa and P. van Loock. Quantum Teleportation and
Entanglement: A Hybrid Approach to Optical Quantum Infor-
mation Processing. Wiley, 2011.
[75] W. Steeb and Y. Hardy. Problems and Solutions in Quan-
tum Computing and Quantum Information. World Scientific,
2006.
[76] M. Gu, C. Weedbrook, N. C. Menicucci, T. C. Ralph, and
P. van Loock. Quantum computing with continuous-variable
clusters. Physical Review A, 79:062318, 2009. doi:
10.1103/PhysRevA.79.062318.
[77] D. Gottesman and I. L. Chuang. Demonstrating the viabil-
ity of universal quantum computation using teleportation
and single-qubit operations. Nature (London), 402:390–393,
1999. doi:10.1038/46503.
[78] S. D. Bartlett and W. J. Munro. Quantum teleportation of
optical quantum gates. Physical Review Letters, 90:117901,
2003. doi:10.1103/PhysRevLett.90.117901.
[79] M. Tillmann, B. Daki
´
c, R. Heilmann, S. Nolte, A. Szameit,
and P. Walther. Experimental boson sampling. Nature Photon-
ics, 7(7):540–544, 2013. doi:10.1038/nphoton.2013.102.
27
[80] N. Spagnolo, C. Vitelli, M. Bentivegna, D. J. Brod, A. Crespi,
F. Flamini, S. Giacomini, G. Milani, R. Ramponi, P. Mataloni,
R. Osellame, E. F. Galvão, and F. Sciarrino. Experimental
validation of photonic boson sampling. Nature Photonics, 8
(8):615–620, 2014. doi:10.1038/nphoton.2014.135.
[81] A. Crespi, R. Osellame, R. Ramponi, D. J. Brod, E. F. Galvão,
N. Spagnolo, C. Vitelli, E. Maiorino, P. Mataloni, and F. Scia-
rrino. Integrated multimode interferometers with arbitrary
designs for photonic boson sampling. Nature Photonics, 7(7):
545–549, 2013. doi:10.1038/nphoton.2013.112.
[82] J. B. Spring, B. J. Metcalf, P. C. Humphreys, W. S. Koltham-
mer, X.-M. Jin, M. Barbieri, A. Datta, N. Thomas-Peter, N. K.
Langford, D. Kundys, J. C. Gates, B. J. Smith, P. G. R.
Smith, and I. A. Walmsley. Boson sampling on a pho-
tonic chip. Science, 339(6121):798–801, 2012. doi:
10.1126/science.1231692.
[83] L. Valiant. The complexity of computing the permanent.
Theoretical Computer Science, 8(2):189–201, 1979. doi:
10.1016/0304-3975(79)90044-6.
[84] E. R. Caianiello. Combinatorics and renormalization in quan-
tum field theory, volume 38. Benjamin, 1973.
[85] A. Barvinok. Combinatorics and complexity of partition func-
tions, volume 274. Springer, 2016.
[86] K. Brádler, P.-L. Dallaire-Demers, P. Rebentrost, D. Su, and
C. Weedbrook. Gaussian boson sampling for perfect match-
ings of arbitrary graphs. Physical Review A, 98:032310, Sep
2018. doi:10.1103/PhysRevA.98.032310.
[87] A. Björklund, B. Gupt, and N. Quesada. A faster hafnian
formula for complex matrices and its benchmarking on the
titan supercomputer. arXiv preprint arXiv:1805.12498, 2018.
[88] N. Quesada, J. M. Arrazola, and N. Killoran. Gaussian boson
sampling using threshold detectors. Physical Review A, 98:
062322, 2018. doi:10.1103/PhysRevA.98.062322.
[89] B. Gupt, J. M. Arrazola, N. Quesada, and T. R. Bromley. Clas-
sical benchmarking of gaussian boson sampling on the titan
supercomputer. arXiv preprint arXiv:1810.00900, 2018.
[90] D. Shepherd and M. J. Bremner. Temporally unstruc-
tured quantum computation. Proceedings of the Royal
Society of London A: Mathematical, Physical and En-
gineering Sciences, 465(2105):1413–1439, 2009. doi:
10.1098/rspa.2008.0443.
[91] M. J. Bremner, A. Montanaro, and D. J. Shepherd. Achieving
quantum supremacy with sparse and noisy commuting quan-
tum computations. Quantum, 1:8, 2017. doi:10.22331/q-
2017-04-25-8.
[92] A. P. Lund, M. J. Bremner, and T. C. Ralph. Quantum sam-
pling problems, BosonSampling and quantum supremacy.
npj Quantum Information, 3(1), 2017. doi:10.1038/s41534-
017-0018-2.
[93] J. M. Arrazola, P. Rebentrost, and C. Weedbrook. Quantum
supremacy and high-dimensional integration. arXiv preprint
arXiv:1712.07288, 2017.
[94] A. Aspuru-Guzik, A. D. Dutoi, P. J. Love, and M. Head-
Gordon. Simulated quantum computation of molecular
energies. Science, 309(5741):1704–1707, 2005. doi:
10.1126/science.1113479.
[95] J. D. Whitfield, J. Biamonte, and A. Aspuru-Guzik. Simu-
lation of electronic structure Hamiltonians using quantum
computers. Molecular Physics, 109(5):735–750, 2011. doi:
10.1080/00268976.2011.552441.
[96] A. M. Childs and N. Wiebe. Hamiltonian simulation using
linear combinations of unitary operations. Quantum Infor-
mation and Computation, 12(11-12):901–924, 2012. doi:
10.26421/QIC12.11-12.
[97] D. W. Berry, G. Ahokas, R. Cleve, and B. C. Sanders. Effi-
cient quantum algorithms for simulating sparse Hamiltoni-
ans. Communications in Mathematical Physics, 270(2):359–
371, 2006. doi:10.1007/s00220-006-0150-x.
[98] T. Kalajdzievski, C. Weedbrook, and P. Rebentrost.
Continuous-variable gate decomposition for the bose-
hubbard model. Physical Review A, 97:062311, Jun 2018.
doi:10.1103/PhysRevA.97.062311.
[99] T. Sowi
´
nski, O. Dutta, P. Hauke, L. Tagliacozzo, and
M. Lewenstein. Dipolar molecules in optical lat-
tices. Physical Review Letters, 108:115301, 2012. doi:
10.1103/PhysRevLett.108.115301.