Fast optimization of parametrized quantum optical circuits
Filippo M. Miatto
1,2
and Nicol
´
as Quesada
3
1
Institut Polytechnique de Paris
2
T
´
el
´
ecom Paris, LTCI, 19 Place Marguerite Perey 91120 Palaiseau
3
Xanadu, Toronto, ON, M5G 2C8, Canada
November 5, 2020
Parametrized quantum optical circuits are a
class of quantum circuits in which the carriers of
quantum information are photons and the gates
are optical transformations. Classically optimiz-
ing these circuits is challenging due to the infi-
nite dimensionality of the photon number vector
space that is associated to each optical mode.
Truncating the space dimension is unavoidable,
and it can lead to incorrect results if the gates
populate photon number states beyond the cut-
off. To tackle this issue, we present an algorithm
that is orders of magnitude faster than the cur-
rent state of the art, to recursively compute the
exact matrix elements of Gaussian operators and
their gradient with respect to a parametrization.
These operators, when augmented with a non-
Gaussian transformation such as the Kerr gate,
achieve universal quantum computation. Our
approach brings two advantages: first, by com-
puting the matrix elements of Gaussian opera-
tors directly, we don’t need to construct them
by combining several other operators; second,
we can use any variant of the gradient descent
algorithm by plugging our gradients into an au-
tomatic differentiation framework such as Ten-
sorFlow or PyTorch. Our results will find appli-
cations in quantum optical hardware research,
quantum machine learning, optical data process-
ing, device discovery and device design.
1 Introduction
In recent years the success of machine learning (ML)
in statistics, data science and artificial intelligence has
sparked much interest in the fields of quantum infor-
mation and quantum computing. Nowadays, quantum
machine learning is a term that encapsulates various
techniques and ideas — including both optimization
of quantum circuits via ML techniques [1, 2], as well
as running subroutines of ML algorithms on quantum
hardware [3–7]. As a parallel thread, the quantum
optics community has investigated how to compound
quantum optical primitives to generate complex quan-
tum states [8–12]. For this program to continue, it be-
comes necessary to develop new and efficient tools to
perform fast simulation of quantum optical circuits.
Parametrized quantum circuits are information pro-
cessing devices whose elements depend on one or more
parameters, and therefore one can think of them as
a particular class of neural networks [13–18]. Among
parametrized quantum circuits, quantum optical cir-
cuits are those in which the carriers of quantum infor-
mation are photons and the gates are optical transfor-
mations [9–11]. Working in the optical domain offers
several advantages such as high processing speed, no
demand for vacuum or cold temperatures and the po-
tential to build compact devices via integrated optics
[19, 19–21].
A typical architecture of quantum optical circuits al-
ternates between Gaussian transformations and non-
linear ones such as the Kerr gate, which is necessary
to achieve universality [22]. Although the evolution of
Gaussian states under Gaussian transformations can be
computed solely from their finite-dimensional symplec-
tic matrix [23, 24], Kerr gates transform Gaussian states
into non-Gaussian ones and so they make it necessary
to work in Fock space, which is infinite-dimensional and
requires one to impose a cutoff.
Optimizing a parametrized quantum optical circuit
is a challenging task, because one needs to simulate ex-
actly the evolution of an input state as it propagates
through the optical components, as well as the gradient
of the output with respect to the circuit parameters.
As each component can be represented in matrix form,
the transformation of the initial state is a simple ma-
trix multiplication and as such it is efficient. However,
the difficult part is to generate the transformation ma-
trix from the parameters of the gate: optical transfor-
mations (even the fundamental ones) are expressed as
matrix exponentials of a linear combination of infinite-
dimensional non-commuting operators [25]. Therefore
(unless the exponent is diagonal) one cannot simply
truncate the matrix at the exponent and subsequently
Accepted in Quantum 2020-11-05, click title to verify. Published under CC-BY 4.0. 1
arXiv:2004.11002v5 [quant-ph] 10 Nov 2020
compute the matrix exponential:
trunc(exp(A)) 6= exp(trunc(A)). (1)
For some special transformations such as the dis-
placement [26] and squeezing operator [27], exact for-
mulas for their matrix elements are known. For some
other transformations, such as two-mode squeezing and
beamsplitters, one can resort to disentangling theorems
which allow one to express the exponential in LDU form
(i.e., as the product of a lower-triangular, a diagonal and
an upper-triangular matrices), which can be safely trun-
cated before multiplication, or alternatively allows the
matrix elements to be expressed as polynomials of the
parameters of the transformations [25, 28, 29]. Bosonic
representations of the special unitary group SU(n) have
also been studied in the mathematical physics literature
[30]: in the quantum optical setting, these are the ma-
trix elements of a passive linear optical transformation
corresponding to an interferometer. In the context of
Franck-Condon transitions, the chemical-physics com-
munity [31–36] has studied the matrix elements of Gaus-
sian transformations where the displacement, squeezing
parameters and unitary matrices representing interfer-
ometers are real, corresponding to so-called point trans-
formations where canonical positions and momenta of
the harmonic oscillator do not mix [37–39]. Until this
work, no exact method was known for more general
Gaussian transformations and typically one expressed
them as a product of fundamental truncated gates, ac-
cepting that truncation errors could propagate.
Our method to perform the optimization of a
parametrized quantum optical circuit has several ad-
vantages over previous techniques:
1. It generates recursively the matrix elements of gen-
eral Gaussian transformations exactly up to the
desired cutoff dimension without the need for de-
compositions. These results are derived in Sec. 2
for general Gaussian gates and applied to specific
cases in Sec. 3.
2. It constructs gradients exactly, without the need
for computational graphs describing the functional
dependencies between the variables. This is ex-
plained in Sec. 4.
3. It is numerically stable and orders of magnitude
faster than previous methods. This is detailed in
Sec. 5 where we report on several numerical exper-
iments.
Even though we focus on Gaussian gates in the main
text, in Appendices C, D, E and F we derive the an-
alytical form of the matrix elements of the cubic- and
quartic- phase gates together with recurrence relations
for their calculation and the calculation of their gra-
dients. These operators are non-Gaussian gates that,
unlike the Kerr gate, are not diagonal in the Fock basis.
The algorithms we present in this manuscript are part of
the fock_gradients module of The Walrus [40] since
version 0.12 released in Dec. 2019.
Finally, note that in this work we construct recursively
the matrix representation of a Gaussian operator. If in-
stead one were interested in just a few specific matrix
elements, one could map this to the calculation of loop
hafnians [41–43].
2 The generating function method
2.1 Notation and definitions
We summarize here our notational conventions, which
are meant to simplify the presence of multiple indices
in the multimode formulas, and the definition of our
generating function.
We adopt the following short-hand notation:
β
n
≡
`
Y
i=1
β
n
i
i
, (2)
n! ≡
`
Y
i=1
n
i
!. (3)
We indicate tensor products of `-mode coherent and
Fock states as follows:
|βi = |β
1
i ⊗ . . . ⊗ |β
`
i, (4)
|ni = |n
1
i ⊗ . . . ⊗ |n
`
i, (5)
which we can use for instance to express a multimode
coherent state:
|βi = e
−
1
2
||β||
2
∞
X
n=0
β
n
√
n!
|ni, (6)
where we used the shorthand notation
P
∞
n=0
=
P
∞
n
1
=0
. . .
P
∞
n
`
=0
. Note that tensor indices are num-
bered starting from 1 (as in Eq. (2)) while the indices
themselves have values that begin from 0 (as in Eq. (6)).
We also use interchangeably the following two notations
for derivatives, depending on convenience and on the
available space:
∂
k
∂β
k
= ∂
k
β
≡
Y
i
∂
k
i
∂β
i
k
i
. (7)
The key insight of the generating function method is
that when we compute the inner product hα
∗
|G|βi of
any operator G (not just a Gaussian one) with a pair
of coherent states, we obtain a function of the coherent
Accepted in Quantum 2020-11-05, click title to verify. Published under CC-BY 4.0. 2
amplitudes that can be treated as a generating function
for the matrix elements of G:
Γ(α, β) = e
1
2
(||α||
2
+||β||
2
)
hα
∗
|G|βi (8)
=
∞
X
m,n=0
α
m
β
n
√
m!n!
hm|G|ni. (9)
The matrix element hm|G|ni is therefore contained in
the coefficient of the factor α
m
β
n
in the Taylor series
and we can isolate it by computing the derivatives of
the appropriate order, at zero:
G
mn
≡ hm|G|ni =
∂
m
α
∂
n
β
Γ(α, β)
√
m!n!
α=β=0
. (10)
We will derive recursion relations for G
mn
directly, in-
stead of calculating the numerator and denominator
separately in the ratio in the last equation. Thus our
methods are numerically more stable than previous ones
where divisions by factorials are common.
We now write the amplitudes and Fock indices as
direct sums of the bra and ket labels
ν =
α
β
∈ C
2`
, (11)
k =
m
n
∈ N
2`
0
, (12)
so that we can express all of our calculations in terms of
a single vector of amplitudes, including the generating
function, which we now indicate as
Γ(ν) =
∞
X
k=0
ν
k
√
k!
G
k
. (13)
where we take advantage of the definition of the vec-
tor k in Eq. (11) and write G
k
for G
mn
. As shown in
Appendix A, the generating function Γ(ν) for a Gaus-
sian operator G is in exponential form Γ(ν) = C e
Q(ν)
.
This implies that we can express high-order derivatives
recursively, using only the derivatives at zero of the ex-
ponent Q(ν), as shown in subsection 2.3. Furthermore,
the exponent Q(ν) is a polynomial of degree 2, which
we can write as:
Q(ν) = µ
T
ν −
1
2
ν
T
Σν. (14)
This means that all of its mixed derivatives of order
higher than 2 vanish:
∂
k
ν
Q(ν)|
ν=0
= 0 if
2`
X
i=1
k
i
> 2. (15)
Because of this, all the derivatives of the exponent that
we will ever need are those of degree 1 and 2:
∂
∂ν
i
Q(ν)|
ν=0
= µ
i
, (16)
∂
2
∂ν
i
∂ν
j
Q(ν)|
ν=0
= −Σ
ij
. (17)
Finally, note that although G
k
is a tensor of rank 2`,
for familiarity with the terminology we still refer to its
elements as “matrix elements”.
2.2 Generating function for multimode Gaussian
transformations
As a corollary of the Bloch-Messiah decomposition [23,
44], a general `−mode Gaussian transformation can be
parametrized as
G = G(γ, W , ζ, V ) = D(γ)U(W )S(ζ)U(V ), (18)
where we used the notation D(γ) =
N
`
i=1
D
i
(γ
i
),
S(ζ) =
N
`
i=1
S
i
(ζ
i
) for tensor products of single-mode
displacement and squeezing operators defined by
D
j
(γ
j
) = exp
γ
j
a
†
j
− H.c.
, γ
j
∈ C, (19)
S
j
(ζ
j
) = exp
ζ
∗
j
2
a
2
j
− H.c.
, ζ
j
= r
j
e
iδ
j
∈ C, (20)
where r
j
and δ
j
are the phase and modulus of the com-
plex squeezing parameter ζ
j
.
In the last equation we used the annihilation a
j
and
creation a
†
j
operators of the ` modes satisfying the
canonical commutation relations
[a
j
, a
†
l
] = δ
j,l
, [a
j
, a
l
] = [a
†
j
, a
†
l
] = 0. (21)
The Hilbert space operator U(V ) represents a gen-
eral multimode passive transformation (physically cor-
responding to an interferometer) and is parametrized
by a unitary matrix V that transforms the creation op-
erators as
U
†
(V )a
†
i
U(V ) =
`
X
l=1
V
∗
il
a
†
l
, U(V )a
†
l
U
†
(V ) =
`
X
i=1
V
il
a
†
i
.
(22)
Using these definitions, in Appendix A we show the
following:
e
1
2
[
||α||
2
+||β||
2
]
hα
∗
|D(γ)U(W )S(ζ)U(V )|βi =
C exp
µ
T
ν −
1
2
ν
T
Σν
, (23)
Accepted in Quantum 2020-11-05, click title to verify. Published under CC-BY 4.0. 3
where
C =
exp
−
1
2
||γ||
2
+ γ
†
W diag(e
iδ
tanh r)W
T
γ
∗
q
Q
`
i=1
cosh r
i
,
(24)
µ
T
= (25)
γ
†
W diag(e
iδ
tanh r)W
T
+ γ
T
, −γ
†
W diag(sech r)V
,
Σ = (26)
"
W diag(e
iδ
tanh r)W
T
−W diag(sech r)V
−V
T
diag(sech r)W
T
−V
T
diag(e
−iδ
tanh r)V
#
.
Note that µ is a vector of dimension 2`. The matrix Σ,
on the other hand is built out of four ` × ` blocks and
has therefore dimension 2` × 2`.
With these results we are now ready to write explicit
recurrence relations for the matrix elements of a general
multimode Gaussian gate.
2.3 From generating functions to matrix ele-
ments
We use the recent result in [45] to express the matrix
elements of the Gaussian operator recursively (see Ap-
pendix B):
G
0
= C, (27)
G
k+1
i
=
1
√
k
i
+ 1
k
X
j=0
s
k!
j!
1
(k − j)!
G
j
∂
k+1
i
−j
ν
Q(ν)|
ν=0
,
(28)
where 1
i
is a vector with all zeros and a single 1 at
position i, and therefore to obtain the vector k + 1
i
we
take the vector k and we increment k
i
by 1.
We stress that Eq. (28) works for any generating
function in exponential form. However, as mentioned
above, in our case the exponent Q is a polynomial of
degree 2, so the recurrence relation simplifies consider-
ably because the derivative of Q is zero unless j = k or
j = k − 1
l
for some index l:
G
k+1
i
=
1
√
k
i
+ 1
G
k
µ
i
−
2`
X
l=1
p
k
l
G
k−1
l
Σ
il
!
. (29)
This, together with Eq. (92) which shows how to ob-
tain gradients, is the primary result of our work. As
an example of Eq. (29), we derive the two recurrence
relations for the single mode case, i.e., ` = 1. For con-
venience and ease of vectorization when programming,
we can use the first relation to build the first column
of the matrix and the second relation to subsequently
build the rows:
G
m+1,n
=
1
√
m + 1
G
m,n
µ
1
−
√
mG
m−1,n
Σ
11
−
√
nG
m,n−1
Σ
12
, (30)
G
m,n+1
=
1
√
n + 1
G
m,n
µ
2
−
√
mG
m−1,n
Σ
21
−
√
nG
m,n−1
Σ
22
. (31)
We represent this example pictorially in Fig. 1.
For Gaussian operators on ` modes, the recurrence
relations allow us to combine 2` + 1 neighbouring ele-
ments in 2` different ways (by incrementing each of the
2` indices) and generate 2` new elements.
Note that these recurrence relations can be rewritten
in an alternative form by using the mode operators:
a
i
G = µ
i
G −
`
X
l=1
a
†
l
Σ
il
!
G − G
`
X
l=1
a
l
Σ
i,l+`
!
, (32)
Ga
†
i
= µ
i+`
G −
`
X
l=1
a
†
l
Σ
i+`,l
!
G − G
`
X
l=1
a
l
Σ
i+`,l+`
!
.
(33)
With these compact relations we can formulate the com-
mutation relations between a general multimode Gaus-
sian operator and the mode operators a
i
and a
†
i
, for
instance for a single-mode Gaussian operator, we have:
[G, a] = −µ
1
G + (Σ
12
+ 1)Ga + Σ
11
a
†
G, (34)
[G, a
†
] = µ
2
G − (Σ
21
+ 1)a
†
G − Σ
22
Ga. (35)
As an example, we can write the special cases of the
commutation relations between the mode operators and
the squeezing and displacement operators (using the
values from subsections 3.1.1 and 3.1.2):
[S(ζ), a] = (1 − sech(r))S(ζ)a + e
iδ
tanh(r)a
†
S(ζ),
(36)
[S(ζ), a
†
] = (sech(r) − 1)a
†
S(ζ) + e
−iδ
tanh(r)S(ζ)a,
(37)
[D(γ), a] = −γD(γ), (38)
[D(γ), a
†
] = −γ
∗
D(γ). (39)
3 Applications: single- and two-mode
gates and passive Gaussian transforma-
tions
All Gaussian single- and two-mode transformations will
be expressed as a combination of four elementary opera-
tions: single-mode phase shifts, single-mode squeezing,
Accepted in Quantum 2020-11-05, click title to verify. Published under CC-BY 4.0. 4
single-mode displacement and beamsplitters. These el-
ements are parametrized as follows: the phase rota-
tion gate R(φ) depends on an angle φ, the single-mode
squeezer S(ζ) depends on a single complex parameter
which we express in polar form ζ = re
iδ
, the displace-
ment gate D(γ) depends on a single complex parame-
ter γ which we leave as is, and finally the beamsplit-
ter B(θ, ϕ) depends on two angles (θ, ϕ). Each gate
is expressed as the exponential of mode operators. In
particular, for the rotation and beamsplitter gates we
write
R(φ) = exp
iφ a
†
a
, (40)
B(θ, ϕ) = exp
h
θ
e
iϕ
a
1
a
†
2
− e
−iϕ
a
†
1
a
2
i
. (41)
For example, the 50/50 beam splitter with real coeffi-
cients corresponds to (θ, ϕ) = (π/4, 0). Note that the
last two gates are special cases of a passive linear opti-
cal transformation, thus they can be written as U(W )
where
W = e
iφ
, (42a)
W =
"
cos θ −e
−iϕ
sin θ
e
iϕ
sin θ cos θ
#
, (42b)
respectively.
3.1 Single-mode Gaussian gate
A simple way to parametrize a general single-mode
Gaussian operation is as follows:
G
(1)
= G
(1)
(γ, φ, ζ) = D(γ)R(φ)S(ζ). (43)
By using Eq. (42a) and the general result from the last
section we easily find
C =
exp
−
1
2
|γ|
2
+ γ
∗2
e
i(δ+2φ)
tanh r
√
cosh r
, (44)
µ
T
= [γ
∗
e
i(δ+2φ)
tanh r + γ, −γ
∗
e
iφ
sech r], (45)
Σ =
"
e
i(δ+2φ)
tanh r −e
iφ
sech r
−e
iφ
sech r −e
−iδ
tanh r
#
. (46)
These values enter the recurrence relations (30) and
(31), which are depicted in Fig. 1.
3.1.1 Special case: single-mode squeezing
For a single-mode squeezing operator, we have φ = 0
and γ = 0, and thus
C =
√
sech r, (47)
µ
T
= [0, 0], (48)
Σ =
"
e
iδ
tanh r −sech r
−sech r −e
−iδ
tanh r
#
. (49)
G
m,n
<latexit sha1_base64="92F8+FGlFZdnU9CMsELx5IgD7eI=">AAAC1nicjVHLSgMxFD0dX7W+qi7dDBbBhZSpD3RZdKFLBWsLVUomjRqcFzMZRUrdiVt/wK1+kvgH+hfexBR8IJphZk7Oveck914/CWSmPO+l4AwNj4yOFcdLE5NT0zPl2bmjLM5TLho8DuK05bNMBDISDSVVIFpJKljoB6LpX+zoePNSpJmMo0N1nYiTkJ1F8lRypojqlGePQ6bOOQt6u/1OL1yJ+p1yxat6Zrk/Qc2CCuzaj8vPOEYXMThyhBCIoAgHYMjoaaMGDwlxJ+gRlxKSJi7QR4m0OWUJymDEXtD3jHZty0a0156ZUXM6JaA3JaWLJdLElJcS1qe5Jp4bZ83+5t0znvpu1/T3rVdIrMI5sX/pBpn/1elaFE6xZWqQVFNiGF0dty656Yq+ufupKkUOCXEadymeEuZGOeizazSZqV33lpn4q8nUrN5zm5vjTd+SBlz7Ps6f4Gi1WlurbhysV+rbdtRFLGARyzTPTdSxh300yPsKD3jEk9Nybpxb5+4j1SlYzTy+LOf+HayslrA=</latexit>
G
m1,n
<latexit sha1_base64="H0nnIDr8uUrGo2qCzSF9MrkZF6w=">AAAC2HicjVHLSsNAFD2Nr1pf1S7dBIvgQkviA10WXeiygm1FW8pkOmpoXiQToZSCO3HrD7jVLxL/QP/CO2MKPhCdkOTMufecmXuvE3luIi3rJWeMjU9MTuWnCzOzc/MLxcWlRhKmMRd1HnphfOqwRHhuIOrSlZ44jWLBfMcTTad3oOLNaxEnbhicyH4k2j67DNwLlzNJVKdYavlMXnHmDQ6HnYG/Ya8Hw06xbFUsvcyfwM5AGdmqhcVntNBFCI4UPgQCSMIeGBJ6zmHDQkRcGwPiYkKujgsMUSBtSlmCMhixPfpe0u48YwPaK89Eqzmd4tEbk9LEKmlCyosJq9NMHU+1s2J/8x5oT3W3Pv2dzMsnVuKK2L90o8z/6lQtEhfY0zW4VFOkGVUdz1xS3RV1c/NTVZIcIuIU7lI8Jsy1ctRnU2sSXbvqLdPxV52pWLXnWW6KN3VLGrD9fZw/QWOzYm9Vdo63y9X9bNR5LGMFazTPXVRxhBrq5N3HAx7xZJwZN8atcfeRauQyTQlflnH/DuS6lyI=</latexit>
G
m,n1
<latexit sha1_base64="JD7FF4c61h6TdgCkFdSmSBOUGB8=">AAAC2HicjVHLSsNAFD2Nr1pf1S7dBIvgQkviA10WXeiygm1FW8pkOmpoXiQToZSCO3HrD7jVLxL/QP/CO2MKPhCdkOTMufecmXuvE3luIi3rJWeMjU9MTuWnCzOzc/MLxcWlRhKmMRd1HnphfOqwRHhuIOrSlZ44jWLBfMcTTad3oOLNaxEnbhicyH4k2j67DNwLlzNJVKdYavlMXnHmDQ6HnYG/HmzYw06xbFUsvcyfwM5AGdmqhcVntNBFCI4UPgQCSMIeGBJ6zmHDQkRcGwPiYkKujgsMUSBtSlmCMhixPfpe0u48YwPaK89Eqzmd4tEbk9LEKmlCyosJq9NMHU+1s2J/8x5oT3W3Pv2dzMsnVuKK2L90o8z/6lQtEhfY0zW4VFOkGVUdz1xS3RV1c/NTVZIcIuIU7lI8Jsy1ctRnU2sSXbvqLdPxV52pWLXnWW6KN3VLGrD9fZw/QWOzYm9Vdo63y9X9bNR5LGMFazTPXVRxhBrq5N3HAx7xZJwZN8atcfeRauQyTQlflnH/DuUylyI=</latexit>
G
m,n+1
<latexit sha1_base64="6dRXhpSoZjxQER8Ie+5Rkk+YZCU=">AAAC2HicjVHLSsNAFD2Nr1pf1S7dBIsgKCXxgS6LLnRZwbaiLWUyHTU0L5KJUErBnbj1B9zqF4l/oH/hnTEFH4hOSHLm3HvOzL3XiTw3kZb1kjPGxicmp/LThZnZufmF4uJSIwnTmIs6D70wPnVYIjw3EHXpSk+cRrFgvuOJptM7UPHmtYgTNwxOZD8SbZ9dBu6Fy5kkqlMstXwmrzjzBofDzsDfCNbtYadYtiqWXuZPYGegjGzVwuIzWugiBEcKHwIBJGEPDAk957BhISKujQFxMSFXxwWGKJA2pSxBGYzYHn0vaXeesQHtlWei1ZxO8eiNSWlilTQh5cWE1WmmjqfaWbG/eQ+0p7pbn/5O5uUTK3FF7F+6UeZ/daoWiQvs6RpcqinSjKqOZy6p7oq6ufmpKkkOEXEKdykeE+ZaOeqzqTWJrl31lun4q85UrNrzLDfFm7olDdj+Ps6foLFZsbcqO8fb5ep+Nuo8lrGCNZrnLqo4Qg118u7jAY94Ms6MG+PWuPtINXKZpoQvy7h/B+BulyA=</latexit>
G
m+1,n
<latexit sha1_base64="duFh6l+/otGnlbKlTHGjOWaYRzE=">AAAC2HicjVHLSsNAFD2Nr1pf1S7dBIsgKCXxgS6LLnRZwbaiLWUyHTU0L5KJUErBnbj1B9zqF4l/oH/hnTEFH4hOSHLm3HvOzL3XiTw3kZb1kjPGxicmp/LThZnZufmF4uJSIwnTmIs6D70wPnVYIjw3EHXpSk+cRrFgvuOJptM7UPHmtYgTNwxOZD8SbZ9dBu6Fy5kkqlMstXwmrzjzBofDzsBftzeCYadYtiqWXuZPYGegjGzVwuIzWugiBEcKHwIBJGEPDAk957BhISKujQFxMSFXxwWGKJA2pSxBGYzYHn0vaXeesQHtlWei1ZxO8eiNSWlilTQh5cWE1WmmjqfaWbG/eQ+0p7pbn/5O5uUTK3FF7F+6UeZ/daoWiQvs6RpcqinSjKqOZy6p7oq6ufmpKkkOEXEKdykeE+ZaOeqzqTWJrl31lun4q85UrNrzLDfFm7olDdj+Ps6foLFZsbcqO8fb5ep+Nuo8lrGCNZrnLqo4Qg118u7jAY94Ms6MG+PWuPtINXKZpoQvy7h/B9/ylyA=</latexit>
Figure 1: The three neighbouring matrix elements shown in
light purple (G
m−1,n
, G
m,n
and G
m,n−1
) can be linearly com-
bined in two different ways to generate G
m+1,n
or G
m,n+1
.
These recurrence rules are also applicable at the edge of the
matrix by considering “outer elements” as zeros. For Gaussian
operators on ` modes, the recurrence rules allow us to combine
2` + 1 neighbouring elements in 2` different ways and generate
2` new elements.
The recurrence relations simplify accordingly, and we
obtain the matrix elements of the single-mode squeezer
S
m,n
= hm|S(ζ)|ni:
S
0,0
=
√
sech r, (50)
S
m+1,0
= −
r
m
m + 1
S
m−1,0
e
iδ
tanh r, (51)
S
m,n+1
=
1
√
n + 1
√
mS
m−1,n
sech r+
√
nS
m,n−1
e
−iδ
tanh r
, (52)
where we have set n = 0 in the first recurrence relation
so that it builds the first column, and we can use the
second relation to build the rows.
Thanks to µ being zero, we do not need to include
the middle matrix element at each step while filling out
the matrix, as depicted in Fig. 2. This gives rise to the
typical checkerboard pattern of the squeezer matrix.
3.1.2 Special case: displacement
For a single-mode displacement operator, we have φ = 0
and ζ = 0 and therefore tanh r = 0 and sech r = 1 and
we find
C = e
−
1
2
|γ|
2
, (53)
µ
T
= [γ, −γ
∗
], (54)
Σ =
"
0 −1
−1 0
#
. (55)
The recurrence relations simplify accordingly, and we
obtain the matrix elements of the single-mode displace-
Accepted in Quantum 2020-11-05, click title to verify. Published under CC-BY 4.0. 5
S
m,n−1
S
m−1,n
S
m,n+1
S
m+1,n
Figure 2: When generating the single-mode squeezer matrix,
two matrix elements are sufficient to generate another two at
each step. Notice how this produces the checkerboard pattern
of zeros that is typical of the squeezer matrix.
ment operator D
m,n
= hm|D(γ)|ni:
D
0,0
= e
−
|γ|
2
2
, (56)
D
m+1,0
=
γ
√
m + 1
D
m,0
, (57)
D
m,n+1
= −
γ
∗
√
n + 1
D
m,n
+
r
m
n + 1
D
m−1,n
. (58)
Thanks to the diagonal of Σ being zero, also in this case
each step is simplified, as depicted in Fig. 3.
D
m−1,n
D
m,n
D
m,n+1
D
m+1,n
D
m,n−1
Figure 3: When generating the single-mode displacement ma-
trix, we only need two matrix elements to generate each new
element along columns or rows.
3.2 Two-mode Gaussian gate
We can write the most general two-mode Gaussian
transformation as
G
(2)
= D(γ)R(φ)B(θ
0
, ϕ
0
)S(ζ)B(θ, ϕ), (59)
where φ = [φ
1
, φ
2
], γ = [γ
1
, γ
2
] and ζ = [ζ
1
, ζ
2
]
and where the single-mode gates factorize over the two
modes. This gives a total of 14 real parameters (count-
ing double for complex parameters), in accordance with
Ref. [46]. The equations above can be recast into
the form of Eq. (18) by writing U(V ) = B(θ, ϕ) and
U(W ) = R(φ)B(θ
0
, ϕ
0
), with
V =
"
cos θ −e
−iϕ
sin θ
e
iϕ
sin θ cos θ
#
, (60)
W =
"
e
iφ
1
0
0 e
iφ
2
#"
cos θ
0
−e
−iϕ
0
sin θ
0
e
iϕ
0
sin θ
0
cos θ
0
#
.
(61)
3.2.1 Special case: Two-mode squeezer
The two-mode squeezer is defined as S
(2)
(ζ) =
exp
ζa
†
1
a
†
2
− H.c.
and can be decomposed in terms of
beamsplitters and single-mode squeezing operations as
follows:
S
(2)
(ζ) = B(−π/4, 0) [S(ζ) ⊗ S(−ζ)] B(π/4, 0). (62)
From the decomposition above one easily finds C =
sech r, µ
T
= [0, 0, 0, 0] and
Σ = (−1)× (63)
0 e
iδ
tanh r sech r 0
e
iδ
tanh r 0 0 sech r
sech r 0 0 −e
−iδ
tanh r
0 sech r −e
−iδ
tanh r 0
.
Note that a two-mode squeezing operation creates pho-
tons in pairs, which implies that the difference of photon
number between the two modes is conserved. This ob-
servation is mathematically equivalent to the fact that
that the operation has an SU(1, 1) symmetry [47]. If we
write S
(2)
m,n,p,q
= hm, n|S
2
(ζ)|p, qi then the only nonzero
elements of this tensor satisfy m − n = p − q. Matrix
elements that do not satisfy this selection rule are zero
and we can construct the whole rank-4 tensor by looping
over at most 3 indices:
S
(2)
0,0,0,0
= sech r, (64)
S
(2)
n,n,0,0
= − S
(2)
n−1,n−1,0,0
Σ
2,1
, (65)
S
(2)
m,n,m−n,0
= −
r
m
m − n
S
(2)
m−1,n,m−n−1,0
Σ
3,1
,
(66)
S
(2)
m,n,p,p−(m−n)
=
−1
p
p − (m −n)
× (67)
√
nS
(2)
m,n−1,p,p−(m−n)−1
Σ
4,2
+
√
pS
(2)
m,n,p−1,p−(m−n)−1
Σ
4,3
.
Accepted in Quantum 2020-11-05, click title to verify. Published under CC-BY 4.0. 6
3.3 Passive Gaussian transformations: Interfer-
ometers
These `-mode transformations, which correspond to in-
terferometers parametrized by an ` × ` unitary matrix
V , can be obtained from our general results by setting
ζ = δ = 0 and selecting, without loss of generality,
W = I to obtain
C = 1, (68)
µ
T
= 0, (69)
Σ = −
"
0 V
V
T
0
#
. (70)
We can write hm|U(V )|ni = U
m,n
and find the follow-
ing recurrence relations
U
0,0
= 1, (71a)
U
m+1
i
,n
=
1
√
m
i
+ 1
`
X
j=1
V
i,j
√
n
j
U
m,n−1
j
, (71b)
U
m,n+1
i
=
1
√
n
i
+ 1
`
X
j=1
V
j,i
√
m
j
U
m−1
j
,n
. (71c)
These equations have a particular simple physical inter-
pretation: they tell us that to scatter off one additional
photon in the bra side m + 1
i
one photon has to be sent
into the interferometer in the ket side n − 1
j
. In par-
ticular, we find that the probability amplitude that a
single photon injected in port j ends up in port i is pre-
cisely the (i, j) entry of the unitary matrix describing
the interferometer:
U
1
i
,1
j
= V
i,j
. (72)
Note that U
m,n
is nonzero only if
P
`
i=1
m
i
=
P
`
i=1
n
i
since interferometers do not create or destroy particles
(mathematically, this is because they have U(n) symme-
try). Matrix elements that do not satisfy this selection
rule are zero and we can construct the whole rank-2`
tensor by looping over at most 2` − 1 indices. Let us
illustrate this for a beamsplitter, which is a special case
of a general passive transformation, and for which V
is of size 2 × 2 and given in Eq. (60). The recurrence
relations we find are:
B
0,0,0,0
=1, (73)
B
m,n,m+n,0
=
1
√
m + n
√
mV
1,1
B
m−1,n,m+n−1,0
+
√
nV
2,1
B
m,n−1,m+n−1,0
,
(74)
B
m,n,p,m+n−p
=
1
√
m + n −p
× (75)
√
mV
1,2
B
m−1,n,p,m+n−p−1
+
√
nV
2,2
B
m,n−1,p,m+n−p−1
,
where hm, n|B(θ, φ)|p, qi = B
m,n,p,q
.
Whenever a gate has a selection rule, not only its
computation is sped up, but also its application when
updating a tensor representing a state. For example, if
we write the state of two modes up to cutoff N as
|ψi =
N−1
X
k,l=0
c
k,l
|k, li, (76)
then the state after the application of a beamsplitter is
given by
|ψ
0
i =B(θ, φ)|ψi =
N−1
X
m,n=0
c
0
n,m
|n, mi, (77)
and because of the U(2) symmetry of the beamsplitter
one can obtain each entry of the updated tensor c
0
using
a single sum:
c
0
n,m
=
min(n+m,N−1)+1
X
k=max(1+n+m−N,0)
B
n,m,k,n+m−k
c
k,n+m−k
.
(78)
Similarly, one can use the SU(1, 1) symmetry of
the two-mode squeezer to show that if |ψ
00
i =
P
N−1
nm=0
c
00
n,m
|n, mi = S
(2)
(z)|ψi then the coefficients
can again be updated in a single sum as follows:
c
00
n,m
=
min(n−m,0)+N
X
k=max(n−m,0)
S
(2)
n,m,k,k+m−n
c
k,k+m−n
. (79)
These computational savings are doubled when consid-
ering density matrices, where updating a single element
of a two-mode state requires in general four sums, but
these sums can be cut in half for the beamsplitter or
the two-mode squeezer.
Accepted in Quantum 2020-11-05, click title to verify. Published under CC-BY 4.0. 7
4 Recursive gradients from generating
functions
In this section we present recurrence relations to con-
struct the gradients of Gaussian operators with respect
to their parameters. Gradients are essential in optimiza-
tion procedures, such as the optimization of a whole
circuit.
In order to carry out the optimization, we define a
function that evaluates the cost or “loss” associated
with the current parameter values. By design, the opti-
mal choice of parameters corresponds to the global min-
imum of the loss function. Our goal is then to minimize
the loss by slowly varying the parameters in a sequence
of steps via a gradient descent algorithm. In order to
apply the gradient descent algorithm, we need to back-
propagate the gradient of the loss function using the
chain rule until we reach each of the parameters. Doing
so gives us the rate of change of the loss with respect to
the parameters, which we can use to update their val-
ues and decrease the loss. The optimization procedure
starts by initializing the parameters randomly (usually
with very small values) and then slowly varies them
until we obtain a circuit that is similar enough to the
optimal one.
In the next subsections we explain how to correctly
handle the gradient of complex and real parameters,
we show that gradients can be easily computed from
the matrix representation of the operators and we give
some examples.
4.1 Gradients with respect to complex parame-
ters
For a complex parameter ξ, the gradient descent update
step should use the partial derivative of a real loss func-
tion L with respect to the conjugate of the parameter
[48, 49]:
ξ ← ξ − γ
∂L
∂ξ
∗
. (80)
Clearly, this update rule falls back to the regular rule in
the case of a real parameter. To compute the gradient
for the update, we need to treat complex variables and
their conjugate as independent variables, which allows
us to compute gradients of non-holomorphic functions
[48]. The chain rule then looks as follows:
∂L
∂ξ
∗
=
X
k
∂L
∂G
∗
k
∂G
∗
k
∂ξ
∗
+
∂L
∂G
k
∂G
k
∂ξ
∗
. (81)
In an automatic differentiation framework such as Ten-
sorFlow or PyTorch, if we wish to customize the com-
putation of the gradient of a new operation (e.g., when
the new operation makes use of compiled code which is
treated like a black box) we are supplied with the up-
stream gradient tensor ∂L/∂G
∗
k
. It is our task to com-
bine it with the local gradients ∂G
k
/∂ξ
∗
and ∂G
∗
k
/∂ξ
∗
as prescribed by the chain rule, to produce the down-
stream gradient ∂L/∂ξ
∗
.
In order to obtain the various parts of Eq. (81) we
proceed as follows: for the first part of (81), the up-
stream gradient tensor ∂L/∂G
∗
k
is given to us by the
software so we only have to compute the local gradient:
∂G
∗
k
∂ξ
∗
=
∂G
k
∂ξ
∗
. (82)
For the second part of (81) we can conjugate the up-
stream gradient tensor (as the loss function is real):
∂L
∂G
k
=
∂L
∂G
∗
k
∗
, (83)
but as the matrix elements of G
k
are complex and non-
holomorphic functions of ξ, their derivatives with re-
spect to ξ and ξ
∗
are independent and we need to com-
pute ∂G
k
/∂ξ
∗
separately. On the other hand, for a real
parameter x computing a single gradient tensor
∂G
k
∂x
would suffice:
∂L
∂x
=
X
k
∂L
∂G
∗
k
∂G
∗
k
∂x
+
∂L
∂G
k
∂G
k
∂x
(84)
= 2<
X
k
∂L
∂G
k
∂G
k
∂x
!
, (85)
as in this case
∂G
∗
k
∂x
and
∂G
k
∂x
are the conjugate of each
other. If one prefers to have only real parameters then
there is an additional step: one can write each complex
parameter in Cartesian or polar form, and then treat
each of these as an additional function of two real vari-
ables, e.g., for ξ = re
iφ
:
∂L
∂r
= 2<
∂L
∂ξ
∗
∂ξ
∗
∂r
= 2<
∂L
∂ξ
∗
e
−iφ
, (86)
∂L
∂φ
= 2<
∂L
∂ξ
∗
∂ξ
∗
∂φ
= −2<
∂L
∂ξ
∗
iξ
∗
. (87)
Furthermore, if the gate depends on a single parameter
(real or complex) more simplifications can be made as
shown in Appendix F.
4.2 Computing gradients with the generating
function method
An `-mode Gaussian gate has 2`
2
+ 3` real parameters
(counting double for complex parameters), therefore in
order to obtain all the parameter updates we need to
Accepted in Quantum 2020-11-05, click title to verify. Published under CC-BY 4.0. 8
Figure 4: Benchmarks for the computation of the squeezing S(ζ), displacement D(γ) and beamsplitter B(θ, ϕ) gates using
Strawberry Fields (SF) 0.12.1 and the methods from this work as implemented in The Walrus (TW). Note that for a cutoff
of N = 30 our implementations are two orders of magnitude faster. Finally, we also benchmarked the two-mode squeezing gate;
the times we find for our implementation are visually indistinguishable from the ones for our beamsplitter, reflecting the fact that
both gates can be implemented by looping over three indices despite being rank-4 tensors. The benchmarks in this figure were
performed on a single core of an Intel-i5 processor @ 2.50 GHz.
compute 2`
2
+3` derivatives, some of which will be com-
bined to form the parameter updates of the complex
parameters.
Another key insight of the generating function
method is that to compute the gradients of the ma-
trix elements we can simply differentiate the generating
function with respect to the parameters before comput-
ing the high-order derivatives (here with respect to a
generic parameter ξ):
∂
ξ
G
k
=
∂
k
ν
√
k!
∂
ξ
Γ(ν)|
ν=0
=
∂
k
ν
√
k!
∂
ξ
Ce
Q(ν)
|
ν=0
(88)
=
∂
k
ν
√
k!
Γ(ν)
∂
ξ
C
C
+ ∂
ξ
Q(ν)
|
ν=0
(89)
=
∂
k
ν
√
k!
Γ(ν)
∂
ξ
C
C
+
∂µ
T
∂ξ
ν −
1
2
ν
T
∂Σ
∂ξ
ν
ν=0
.
(90)
Expression (90) contains terms proportional to ν
n
Γ(ν),
whose derivative can be expressed using the elements of
G
k
:
∂
k
ν
√
k!
ν
n
Γ(ν)|
ν=0
=
p
(k)
n
G
k−n
, (91)
where (k)
n
= k(k − 1)(k − 2) . . . (k − n + 1) is the
Pochhammer symbol and (k)
n
=
Q
i
(k
i
)
n
i
(with the
convention that (k)
0
= 1). So once we compute the ten-
sor representation of a Gaussian operator, we can use it
to compute its gradient with respect to a parametriza-
tion. Continuing from Eq. (90) we can write:
∂G
k
∂ξ
=
∂
ξ
C
C
G
k
+
X
i
∂µ
i
∂ξ
p
k
i
G
k−1
i
−
X
i>j
∂Σ
ij
∂ξ
p
k
i
k
j
G
k−1
i
−1
j
−
1
2
X
i
∂Σ
ii
∂ξ
p
k
i
(k
i
− 1)G
k−2
i
. (92)
This, together with Eq. (29), is the primary result of
our work.
As an example, we show the derivatives of the single
mode Gaussian operator with respect to the complex
displacement parameter γ and γ
∗
(recall that we must
treat them as independent variables). We begin by dif-
Accepted in Quantum 2020-11-05, click title to verify. Published under CC-BY 4.0. 9
ferentiating C, µ and Σ:
∂
γ
C
C
= −
γ
∗
2
, (93)
∂
γ
∗
C
C
= −
γ
2
− γ
∗
e
i(δ+2φ)
tanh r, (94)
∂µ
∂γ
= [1, 0], (95)
∂µ
∂γ
∗
= [e
i(δ+2φ)
tanh r, −e
iφ
sech r], (96)
∂Σ
∂γ
= 0, (97)
∂Σ
∂γ
∗
= 0. (98)
Therefore, the gradients of the matrix elements are:
∂G
m,n
∂γ
= −
γ
∗
2
G
m,n
+
√
m G
m−1,n
, (99)
∂G
m,n
∂γ
∗
= −
γ
2
+
1
2
γ
∗
e
i(δ+2φ)
tanh r
G
m,n
(100)
+ e
i(δ+2φ)
√
mG
m−1,n
tanh r
− e
iφ
√
nG
m,n−1
sech r.
We then combine both of these expressions with the
upstream gradient ∂L/∂G
∗
mn
and its conjugate as pre-
scribed in Eq. (81) to obtain the single gradient for the
update, ∂L/∂γ
∗
. We proceed in an analogous way for
all the other parameters.
5 Numerical experiments
In this section we present a series of experiments that
show the effectiveness of our method. First, we bench-
mark our methods and implementations for the dis-
placement, squeezing and beamsplitter gate against
the implementations of Strawberry Fields (SF) ver-
sion 0.12.1. We note that as of version 0.13, SF uses
the methods presented here and implemented in The
Walrus. The results of the benchmarks are summa-
rized in Fig. 4; our implementations are two orders of
magnitude faster than the ones in SF 0.12.1, and are
also more stable and memory efficient. For example, if
one tries to calculate a displacement gate by using its
known expansion in terms of Laguerre polynomials [26]
it is necessary to perform divisions by the factorial of
the number of photons n, which results in numerical
overflows for n ∼ 100. Our explicit recurrence relations
avoid completely this issue for any number of modes.
Second, we compare our methods against the im-
plementations of the beamsplitter in the package
Bosonic[11]. This package has a hard-coded cutoff of
11. We benchmark up to this cutoff in Fig. 5. Even for
Figure 5: Benchmarks for the computation of the beamsplit-
ter BS(θ, ϕ) using Bosonic and the methods from this work
as implemented in The Walrus. Note that for a given cutoff
N, Bosonic calculates B
(N)
i,j
= hN − i, i|BS(θ, ϕ)|N − j, ji
and that in their implementation the constraint N ≤ 11
is hard-coded. For a fair comparison we benchmark their
code against a wrapper function that calculates the tensor
hm, n|BS(θ, ϕ)|k, li and then constructs the relevant matrix
B
(N)
i,j
from it. Note that our methods offer a very significant
improvement in speed even for the modest cutoff considered
here. The benchmarks in this figure were performed on a sin-
gle core of an Intel-i5 processor @ 2.50 GHz.
these modest sizes our implementation achieves almost
two orders of magnitude in improvement.
Third, we test a circuit optimization task. We start
by setting up a circuit built by alternating single-mode
Gaussian gates and single-mode Kerr gates. A Kerr
gate is a non-linear optical interaction which varies the
phase of a Fock state according to the square of the
occupation number:
K(κ) = exp
iκ(a
†
a)
2
, (101)
and is therefore diagonal in the photon number basis.
Such single-mode circuit corresponds to a unitary trans-
formation which we can write as a stack of M layers:
U(γ, φ, ζ, κ) =
M
Y
m=1
G
(1)
(γ
m
, φ
m
, ζ
m
)K(κ
m
). (102)
We then input the single-mode vacuum state |0i and op-
timize the parameters such that the unitary U(γ, φ, ζ)
produces an output as close as possible to the desired
target state. We choose three target states: a single
photon state |1i, an “ON” state, i.e. a superposition of
the vacuum and a Fock state, in our case (|0i+ |9i)/
√
2
and finally a Hex GKP state, which is a logical state
of the single-mode GKP error correcting code [50]. To
evaluate the quality of the state produced by the circuit
we define a loss function between the output U|0i and
Accepted in Quantum 2020-11-05, click title to verify. Published under CC-BY 4.0. 10
Hyperparams. Single photon ON Hex GKP
Cutoff N 100 (8) 100 (14) 100 (50)
Layers M 8 20 35 (25)
GD Steps 1500 (5000) 2500 (5000) 5000 (10000)
Results
Fidelity 99.998% 99.995 (99.93)% 99.83 (99.60)%
Runtime (s) 50 (65) 260 (436) 720 (6668)
Table 1: Optimization results for state generation. We indicate in parenthesis the values from Table I in Ref. [8] if they differ
from ours. Thanks to the improved speed of our recursive algorithm we can push the cutoff to N = 100 and we can increase the
number of layers for the generation of the Hex GKP state. Compared to the results in Table I of Ref. [8] we can either reach
the same fidelity with fewer steps or overall better results (higher fidelity, fewer steps, shorter runtime). All of these optimizations
have been run on a single core of an Intel Core i5 @ 3,1 GHz. For comparison the results in Ref. [8] are obtained in hardware of
comparable speed for the first two target states and on a 20-core Intel Xeon CPU @ 2.4GHz with 252GB RAM for the Hex GKP
state.
the target:
L(|ψ
out
i) = −|hψ
out
|ψ
target
i|. (103)
For the optimization we run the adaptive gradient de-
scent algorithm Adam. We note that the range over
which the parameters are optimized is unconstrained,
except for the fact that the parameters φ
m
and κ
m
are
real.
The results are presented in Table 1 and can be com-
pared with the ones in table I of Ref. [8] (reported in
parenthesis where they differ from ours). Our method
allows us to use much higher cutoff dimensions, leading
to equal or better fidelities in a significantly shorter time
or fewer gradient descent steps (GD Steps in the table).
The optimization of the Hex GKP state is particularly
striking, as it reached better results in a fraction of the
time while running on much slower hardware.
We point out that even if the final output has to
contain few photons, using a large cutoff is beneficial
because it allows intermediate layers to explore higher
excitation numbers, which may be necessary to reach
the desired output with high fidelity.
6 Conclusions
We have presented a unified procedure for obtaining
the Fock representation of a quantum optical Gaus-
sian transformation using generating functions. From
the generating functions we obtain recurrence rela-
tions which allow us to compute the matrix elements
of any given Gaussian transformation when written in
terms of its Bloch-Messiah decomposition. Further-
more, we have developed general techniques to also ob-
tain the gradients of a given transformation once it is
parametrized.
As particular examples we derived explicit recurrence
relations for the displacement, single- and two-mode
squeezing gates and the beamsplitter. These recurrence
relations are part of the fock_gradients module of
The Walrus written in pure Python using Numpy and
are sped up using the just-in-time compiling capabilities
of Numba. Our highly portable implementation makes
our algorithms ideal for high performance simulation of
quantum optical circuits using both CPUs and GPUs.
Furthermore, we expect them to accelerate the research
on quantum hardware, quantum machine learning, opti-
cal data processing, device discovery and device design.
Note added – While preparing this manuscript we be-
came aware of related work by J. Huh [51].
Acknowledgements
N.Q. thanks Safwan Hossain, Theodor Isacsson, Josh
Izaac, Nathan Killoran, John E. Sipe and Antal Sz´ava
for valuable discussions. F.M. thanks Electra Elefthe-
riadou for support throughout this work, Yuan Yao for
helping with the correctness of the final version of the
manuscript, Fr´ed´eric Grosshans, Steve Barnett, Ger-
ardo Adesso and Alessio Serafini for helpful inputs. The
authors also thank the open source scientific computing
community, in particular the developers of Numpy [52],
Scipy [53], Jupyter [54], Matplotlib [55] and Numba
[56], without whom this research would not have been
possible.
Accepted in Quantum 2020-11-05, click title to verify. Published under CC-BY 4.0. 11
A Multimode generating function
In Barnett and Radmore [25] (cf. Sec. 3.7 and Appendix 5) the following identities are proven:
exp
θa
†
a
= : exp
[exp(θ) − 1] a
†
a
:=
∞
X
n=0
[exp(θ) − 1]
n
a
†n
a
n
n!
, (104)
and also
S(ζ) = exp
1
2
ζ
∗
a
2
− H.c.
= exp
1
2
re
−iδ
a
2
− H.c.
= exp
−
tanh r
2
e
iδ
a
†2
exp
−
a
†
a +
1
2
log cosh r
exp
tanh r
2
e
iδ
a
2
. (105)
where we wrote ζ = re
iδ
. Putting these two identities together we can normal order a squeezing operator
exp
1
2
re
−iδ
a
2
− H.c.
= exp
−
tanh r
2
e
iδ
a
†2
×
: exp
[sech r − 1] a
†
a
:
√
cosh r
× exp
tanh r
2
e
−iδ
a
2
. (106)
Using this expression one can find the matrix element of the squeezing operation between two coherent states
hα
∗
|exp
1
2
re
−iδ
a
2
− H.c.
|βi =
1
√
cosh r
hα
∗
|βiexp
−
tanh r
2
e
iδ
α
2
+ [sech r − 1] αβ +
tanh r
2
e
−iδ
β
2
, (107)
and noting that
hα
∗
|βi = exp
−
1
2
|α|
2
+ |β|
2
− 2αβ
, (108)
one can write
hα
∗
|exp
1
2
re
iδ
a
2
− H.c.
|βi =
exp
−
1
2
|α|
2
+ |β|
2
√
cosh r
exp
−
1
2
[α β]
"
e
iδ
tanh r −sech r
−sech r −e
−iδ
tanh r
#"
α
β
#!
.
(109)
This proof generalizes in a straightforward way to
hα
∗
|S(ζ)|βi =
exp
−
1
2
||α||
2
+ ||β||
2
q
Q
`
i=1
cosh r
i
exp
−
1
2
[α β]
"
diag(e
iδ
tanh r) −diag(sech r)
−diag(sech r) −diag(e
−iδ
tanh r)
#"
α
β
#!
,
where ζ
j
= r
j
e
iδ
j
.
From the definitions in Eq. (22) it is direct to show that
U(V )|βi = |V βi and hα
∗
|U(W ) = hW
†
α
∗
|. (110)
Finally, let us note that
hα
∗
|D(γ) = hα
∗
− γ|exp
1
2
α
T
γ − γ
†
α
∗
. (111)
To show this last identity recall that hα
∗
| = h0|D(−α
∗
) and use the composition rule for displacement operators
(cf. Eq. 3.6.30 of Ref. [25]). With this setup we are finally ready to prove Eq. 23 as follows:
hα
∗
|D(γ)U(W )S(ζ)U(V )|βi (112a)
= exp
1
2
α
T
γ − γ
†
α
∗
hα
∗
− γ|U(W )S(ζ)|V βi (112b)
= exp
1
2
α
T
γ − γ
†
α
∗
hW
†
(α
∗
− γ)|S(ζ)|V βi (112c)
=
exp
1
2
α
T
γ − γ
†
α
∗
q
Q
`
i=1
cosh r
i
exp
−
1
2
||α
∗
− γ||
2
+ |β|
2
× (112d)
exp
−
1
2
(α
T
− γ
†
)W , β
T
V
T
"
diag(e
iδ
tanh r) −diag(sech r)
−diag(sech r) −diag(e
−iδ
tanh r)
#"
W
T
(α − γ
∗
)
V β
#!
.
Accepted in Quantum 2020-11-05, click title to verify. Published under CC-BY 4.0. 12
At this point we introduce the symmetric, complex-unitary matrix
Σ =
"
W diag(e
iδ
tanh r)W
T
−W diag(sech r)V
−V
T
diag(sech r)W
T
−V
T
diag(e
−iδ
tanh r)V
#
(113)
=
"
W 0
0 V
T
#"
diag(e
iδ
tanh r) −diag(sech r)
−diag(sech r) −diag(e
−iδ
tanh r)
#"
W 0
0 V
T
#
T
, (114)
and the complex vector
ν =
"
α
β
#
∈ C
2`
, (115)
to finally write
hα
∗
|D(γ)U(W )S(ζ)U(V )|βi = exp
−
1
2
||α||
2
+ ||β||
2
C exp
µ
T
ν −
1
2
ν
T
Σν
, (116)
where
C =
exp
−
1
2
||γ||
2
+ γ
†
W diag(e
iδ
tanh r)W
T
γ
∗
q
Q
`
i=1
cosh r
i
, (117)
µ
T
=
γ
†
W diag(e
iδ
tanh r)W
T
+ γ
T
, −γ
†
W diag(sech r)V
, (118)
which is the expression quoted in the main text. Note that a related expression for the matrix elements of a general
Gaussian transformations evaluated in the coordinate basis was worked out by Moshinsky and Quesne [37].
B Recurrence relation for the high-order derivatives of an exponential function
We rewrite Eq. (15) of Ref. [45] in our notation as follows:
Γ
k+1
n
=
k
X
j=0
k
j
Q
(k−j+1
n
)
Γ
j
, (119)
where we write ∂
k
ν
f(ν)|
ν=0
= f
(k)
= f
k
for the k
th
derivative of f at zero and define
P
k
j=0
k
j
=
P
k
1
j
1
···
P
k
n
j
n
k
1
j
1
. . .
k
n
j
n
. In our case Q is a quadratic in ν and we write
Q(ν) = µ
T
ν −
1
2
ν
T
Σν. (120)
We can write easily all its nonzero derivatives:
Q
(1
i
)
= µ
i
, (121)
Q
(1
i
+1
j
)
= −Σ
i,j
. (122)
We know that most of the terms in Eq. (119) will not survive the summation, indeed only the ones that satisfy
k − j = 1
i
or k − j = 1
i
+ 1
l
will contribute to the sum. We can use this to our advantage and write
Γ
k+1
n
= Q
(1
n
)
Γ
k
+
X
l
k
l
Q
(1
n
+1
l
)
Γ
k−1
l
= µ
n
Γ
k
−
X
l
k
l
Σ
l,n
Γ
k−1
l
. (123)
Now let us derive this in a different way. Start with the generating function of the multidimensional Hermite
polynomials [57–59]:
exp
µ
T
ν −
1
2
ν
T
Σν
=
∞
X
k=0
ν
k
k!
G
(Σ)
k
(µ). (124)
Accepted in Quantum 2020-11-05, click title to verify. Published under CC-BY 4.0. 13
They satisfy the recurrence relation
G
(Σ)
k+1
i
(µ) = µ
i
G
(Σ)
k
(µ) −
X
j=1
Σ
i,j
k
j
G
(Σ)
k−1
j
(µ),
from which we identify that for a quadratic function Q,
Γ
k
= G
(Σ)
k
(µ) =
√
k!G
k
. (125)
By replacing this last expression in the recurrence relation Eq. (123), and simplifying the factorials one arrives at
the recurrence relation Eq. (28) for the matrix elements G
k
.
C Recursive matrix elements from commutation relations
In this section we use an alternative approach for obtaining the matrix elements and gradients of a gate by using
the commutation relations with the destruction operator of a given mode. To illustrate this method we will focus
on the k
th
-order-phase gate defined as
V
(k)
(η) = exp
i
η
k~
q
k
, (126)
where the canonical position is defined as q =
q
~
2
(a+a
†
) and its canonical momentum is given by p = i
q
~
2
a
†
− a
.
In the Heisenberg picture the k
th
-order-phase gate transforms the operators according to
V
(k)†
(η)qV
(k)
(η) = q, (127)
V
(k)†
(η)pV
(k)
(η) = p + ηq
k−1
, (128)
V
(k)†
(η)aV
(k)
(η) = a + i
η
~
~
2
k/2
a + a
†
k−1
. (129)
One can write the matrix elements of the k
th
-order-phase gate in terms of integrals as follows
hm|V
(k)
(η)|ni = V
(k)
m,n
=
1
√
π2
n+m
n!m!
Z
dx exp
−x
2
+ i
η~
k/2
~k
x
k
H
n
(x)H
m
(x), (130)
where H
n
(x) is a Hermite polynomial of degree n [25]. The last equation makes explicit that the matrix elements
V
(k)
m,n
are symmetric V
(k)
m,n
= V
(k)
n,m
and moreover, given the parity of the Hermite polynomials H
n
(−x) = (−1)
n
H
n
(x),
that the k
th
-order-phase gate also has parity, i.e.,
V
(k)
m,n
(η) = (−1)
m+n
V
(k)
m,n
(−η) if k is odd, (131)
V
(k)
m,n
(η) = 0 if k is even and n + m is odd. (132)
Finally, note that the matrix elements of the conjugate k
th
-order-phase gate
˜
V
(k)
(η) = exp
i
η
~k
p
k
, (133)
are related to the ones for the k
th
-order-phase gate via
hm|
˜
V
(k)
(η)|ni =
˜
V
(k)
m,n
= (i)
m+n
V
(k)
m,n
, (134)
which can be easily confirmed by noticing that the rotation gate by angle θ = π/2 transforms q into p in the
Heisenberg picture and that when acted on a Fock state it gives R(π/2)|ni = (i)
n
|ni.
Having established these basic facts about the k
th
-order-phase gate we are now ready to write a recurrence
relation for its matrix elements. To this end, we write the following commutator
[a, V
(k)
(η)] = aV
(k)
(η) − V
(k)
(η)a = V
(k)
(η)V
(k)†
(η)aV
(k)
(η) − V
(k)
(η)a = V
(k)
(
i
η
~
~
2
k/2
(a + a
†
)
k−1
)
.
(135)
Accepted in Quantum 2020-11-05, click title to verify. Published under CC-BY 4.0. 14
If we multiply the last equation on the left by hm| and on the right |ni we obtain
√
m + 1V
(k)
m+1,n
−
√
nV
(k)
m,n−1
= hm|V
(k)
(η)
i
η
~
~
2
k/2
(a + a
†
)
k−1
!
|ni. (136)
Note that for any k one can expand the right hand side of the last equation and write (a + a
†
)
k−1
|ni =
P
n+k−1
k=n−(k−1)
c
k
|k + ni allowing us to obtain the desired recurrence relation. In the next subsection we special-
ize these results for the values k = 3 and k = 4. Note that the cases k = 1, 2 are Gaussian gates and thus can be
dealt with using the methods of the main text.
C.1 The cubic-phase gate
The cubic phase gate, introduced in Ref. [50], is the simplest non-Gaussian operation one can consider when
studying polynomial generators of continuous-variable quantum gates. This gate is obtained by setting k = 3
in Eq. (126) and can in principle be used to decompose any higher order gate generated by a polynomial in the
canonical position q and momentum p of a continuous-variable system by employing approximate commutator
expansions [22, 24, 60] or, for many important applications, finding exact decompositions [61, 62]. Significant effort
has been directed into realizing this gate either by using optical nonlinearities [63] or measurements [50, 64–67].
To obtain the recurrence relation for the cubic phase gate, we set k = 3 in Eq. (136) and rearrange to obtain
V
(3)
m,n+2
=
1
p
(n + 1)(n + 2)
i2
√
2
√
~η
n
√
nV
(3)
m,n−1
−
√
m + 1V
(3)
m+1,n
o
− (2n + 1)V
(3)
m,n
−
p
n(n − 1)V
(3)
m,n−2
. (137)
We can obtain the first two rows of V
(3)
m,n
by setting m = 0 or 1 in the last equation and using the following initial
seeds (cf. Appendix D)
V
(3)
0,0
= 2
√
πe
2y
6
3
|y|Ai
y
4
, (138)
V
(3)
1,1
= −2i
√
2y
3
V
(3)
0,1
= −2i
√
2y
3
V
(3)
1,0
= −8
√
πe
2y
6
3
|y|
5
Ai
0
y
4
+ y
2
Ai
y
4
, (139)
where y
3
= 1/(
√
~η) and Ai(x) and Ai
0
(x) are the Airy function and its first derivative [68, 69]. Since V
(k)
m,n
is
symmetric in m, n once the first two rows are obtained one also obtains the first two columns which serve as the
initial conditions to recursively obtain any other row and thus the whole matrix.
Even though the numerical recursive relation derived are exact, they are not numerically stable in the regime
where η < 1; notice that any small numerical error in the difference
√
nV
(3)
m,n−1
−
√
m + 1V
(3)
m+1,n
inside the curly
braces in Eq. 137 will be amplified by a factor of 1/η that will cause an uncotrollable loss of precision; because of
this, the recurrence relation can only be safely used in the regime where η > 1.
C.2 The quartic-phase gate
The quartic phase interaction arises in the context of self-interactions in scalar quantum field theories [70, 71] and
can also be use to achieve universality in GKP encoded qubits[50, 72]. We can obtain a recursive relation for the
matrix elements of this gate by setting k = 4 in Eq. (136) to obtain
V
(4)
m+3,n
=
1
p
(n + 1)(n + 2)(n + 3)
4i
η~
n
√
nV
(4)
m,n−1
−
√
m + 1V
(4)
m+1,n
o
−
p
n(n − 1)(n −2)V
(4)
m,n−3
− 3n
3
2
V
(4)
m,n−1
− 3(n + 1)
3
2
V
(4)
m,n+1
. (140)
Accepted in Quantum 2020-11-05, click title to verify. Published under CC-BY 4.0. 15
We can iterate this relation using the following initial values
V
(4)
0,0
= f(w, λ)|
λ=1
, (141)
V
(4)
1,1
= −2
∂
∂λ
f(w, λ)
λ=1
, (142)
V
(4)
0,2
=V
(4)
2,0
=
1
√
2
V
(4)
1,1
− V
(4)
0,0
, (143)
V
(4)
2,2
=
1
2
V
(4)
0,0
− V
(4)
1,1
+ 2
∂
2
∂λ
2
f(w, λ)
λ=1
, (144)
where f(w, λ) is defined in Eq. (153) and recalling that V
(4)
0,1
= V
(4)
1,0
= V
(4)
1,2
= V
(4)
2,1
= 0 by parity. With these nine
initial values we can fill the first three rows of V
(4)
m,n
using the recurrence relation. Then we can use the symmetry
to fill the first three columns and finally complete the whole matrix. Note however that, like its cubic counterpart,
the recurrence relation in Eq. (140) is numerically unstable for η < 1.
D Analytic matrix elements of the cubic-phase gate
In this appendix we obtain the matrix elements of the cubic-phase gate. We start from Eq. (130) with k = 3 and
y
3
= 1/(
√
~η) and doing the change of variables x → yk − iy
3
, and assuming without loss of generality that y > 0,
we obtain
V
(3)
m,n
=
ye
2
3
y
6
√
π2
n+m
n!m!
Z
∞
−∞
dk exp
iy
4
k + ik
3
/3
H
n
(yk − iy
3
)H
m
(yk − iy
3
). (145)
To make progress we recall the identity H
n
(x + y) =
P
n
k=0
n
k
H
k
(x)(2y)
(n−k)
, to obtain
V
(3)
m,n
=
ye
2
3
y
6
√
π2
n+m
n!m!
Z
∞
−∞
dk exp
iy
4
k + ik
3
/3
n
X
k=0
n
k
H
k
(−iy
3
)(2yk)
(n−k)
m
X
l=0
m
l
H
l
(−iy
3
)(2yk)
(m−l)
(146)
=
ye
2
3
y
6
√
π2
n+m
n!m!
n
X
k=0
m
X
l=0
(2y)
m+n−k−l
n
k
H
k
(−iy
3
)
m
l
H
l
(−iy
3
)
Z
∞
−∞
dk exp
iy
4
k + ik
3
/3
k
n+m−k−l
(147)
=
2
√
πye
2
3
y
6
√
2
n+m
n!m!
n
X
k=0
m
X
l=0
(−2iy)
m+n−k−l
n
k
H
k
(−iy
3
)
m
l
H
l
(−iy
3
)
d
n+m−l−k
dw
n+m−l−k
Ai(w)
w=y
4
. (148)
In the last equation, we used the well known Fourier expansion of the Airy function to write [69]
Ai(x) =
1
2π
Z
∞
−∞
dk exp
i
k
3
3
+ ikx
, (149)
d
`
dx
`
Ai(x) =
1
2π
Z
∞
−∞
dk(ik)
`
exp
i
k
3
3
+ ikx
. (150)
Note that any higher order derivative of the Airy function can be expressed in terms of polynomials of its arguments,
the Airy function itself and it first derivative as follows
d
`
dx
`
Ai(x) = p
`
(x) Ai(x) + q
`
(x)Ai
0
(x), (151)
where the polynomials p
`
and q
`
satisfy the recurrence relations
p
`+1
(x) = xq
`
(x) +
d
dx
p
`
(x), q
`+1
(x) = p
`
(x) +
d
dx
q
`
(x). (152)
By iterating the recurrence relations in Eq. (152) one can in principle obtain the polynomials p
`
(x) and q
`
(x) and
then any derivative of the Airy function and thus any matrix element V
(3)
n,m
.
Accepted in Quantum 2020-11-05, click title to verify. Published under CC-BY 4.0. 16
E Analytic matrix elements of the quartic-phase gate
We are now interested in calculating
V
(4)
m,n
= hm|V
(4)
(η)|ni =
1
√
π2
n+m
n!m!
Z
∞
−∞
dx exp
−x
2
+
i
4
wx
4
H
n
(x)H
m
(x),
where H
n
(x) is the n
th
order Hermite polynomial and we introduced w = ~η. First consider the related integral
(cf. Eq. 3.968 of Ref. [73])
f(w, λ) =
1
√
π
Z
∞
−∞
dx exp
−λx
2
+
i
4
wx
4
(153)
= −
i
4
√
πe
iλ
2
8w
+
iπ
8
r
λ
w
J
1
4
λ
2
8w
− iY
1
4
λ
2
8w
= −
i
4
√
πe
iλ
2
8w
+
iπ
8
r
λ
w
H
(2)
1
4
λ
2
8w
, (154)
where J
α
(x) is a Bessel function of the first kind, Y
α
(x) is a Bessel function of the second kind (sometimes also
called a Neumann function) and H
(2)
α
(x) is a Hankel function of the second kind (not to be confused with the
Hermite polynomials H
n
(x)). Note that f(w, 1) = V
(4)
0,0
(w). To obtain any other matrix element first notice that
these elements are nonzero if and only if n + m is even, which implies that H
n
(x)H
m
(x) contains only even powers
of x. Thus we only need to know how to evaluate
(−1)
k
√
π
Z
∞
−∞
dx x
2k
exp
−λx
2
+
i
4
wx
4
=
∂
k
∂λ
k
f(w, λ). (155)
In turn the derivatives of this function can be written as derivatives of the Hankel functions of the second kind
which satisfy the recurrence relation
∂
∂x
H
(2)
i
(x) =
1
2
H
(2)
i−1
(x) + H
(2)
i+1
(x)
. (156)
Thus, any matrix element V
(4)
m,n
can be obtained in terms of exponentials, polynomials and Hankel functions of the
second kind in the variable w.
F Gradients of single parameter gates
We outline here an alternative procedure to obtain gradients of single-parameter gates. Let us consider a gate
parametrized by a single complex number u = xe
i
,
G(u) = G(xe
i
) = exp
xe
i
h − H.c.
. (157)
The case of a gate parametrized by a single real number can be obtained by letting = π/2 and writing G(x) =
exp(ix(h + h
†
)) = exp(Hx) with H = i(h + h
†
) = −H
†
. For example, for a two-mode squeezing gate, h = a
†
1
a
†
2
and
for the cubic phase gate, H = iq
3
/(k~). Now it is direct to verify that
∂
x
hm|G|ni = hm|G
e
i
h − H.c.
|ni. (158)
To obtain the gradient with respect to the amplitude x, one only needs to express (e
i
h − H.c.)|ni in terms of a
linear combination of Fock states to reduce the gradient to a linear combination of the matrix elements of the gate
itself. Again using the two-mode squeezing gate and the cubic-phase gate as examples we easily find
∂
r
hm|S
(2)
|ni = hm
1
, m
2
|S
(2)
n
a
†
1
a
†
2
e
i
− H.c.
o
|n
1
, n
2
i (159)
=
p
(n
1
+ 1)(n
2
+ 1)e
i
S
(2)
m
1
,m
2
,n
1
+1,n
2
+1
−
√
n
1
n
2
e
−i
S
(2)
m
1
,m
2
,n
1
−1,n
2
−1
,
∂
η
hm|V
(3)
|ni = hm|iV
(3)
√
~
3 × 2
3/2
(a + a
†
)
3
|ni (160)
= i
√
~
3 × 2
3/2
p
n(n − 1)(n −2)V
(3)
m,n−3
+ 3n
3
2
V
(3)
m,n−1
+ 3(n + 1)
3
2
V
(3)
m,n+1
+
p
(n + 1)(n + 2)(n + 3)V
(3)
m,n+3
.
Accepted in Quantum 2020-11-05, click title to verify. Published under CC-BY 4.0. 17
Gate Parameter u = xe
i
Generator h Value of s
D(γ) γ a
†
1
1
S(ζ) ζ = re
iδ
−
1
2
a
†2
1
1
2
S
(2)
(ζ) ζ = re
iδ
a
†
1
a
†
2
1
B(θ, ϕ) θe
iϕ
a
1
a
†
2
−1
Table 2: Parametrizations of different common Gaussian gates and their s-parameter necessary to obtain the gradient with respect
to their phase.
Let us consider the gradient with respect to the phase: one would hope that a formula similar to Eq. (158) should
hold. Unfortunately this is not the case, the reason being that two gates with the same phase parameter commute
with each other, [G(xe
i
), G(x
0
e
i
)] = 0 but the same commutativity does not hold for different values of the phase,
[G(xe
i
), G(xe
i
0
)] 6= 0. To make progress we note that for any of the single complex-parameter gates one can obtain
the following simple decomposition:
G(xe
i
) = R
1
(s)G(x)R
†
1
(s), (161)
where R
1
(s) = e
isa
†
1
a
1
is a rotation gate. We summarize the values of the different parameters in this decomposition
in Table 2 for the single- and two-mode gates studied in the main text. With this decomposition we can now take
derivatives of matrix elements with respect to the phase argument of the complex parameter,
∂
hm|G|ni =hm|(∂
R
1
(s))G(x)R
†
1
() + R
1
(s)G(x)(∂
R
†
1
(s))|ni
=hm|(isa
†
1
a
1
)R
1
(s)G(x)R
†
1
(s) + R
1
(s)G(x)R
†
1
(s)(−isa
†
1
a
1
)|ni
=is(m
1
− n
1
)hm|G|ni. (162)
References
[1] Maria Schuld, Ville Bergholm, Christian Gogolin,
Josh Izaac, and Nathan Killoran. Evaluating an-
alytic gradients on quantum hardware. Phys.
Rev. A, 99(3):032331, 2019. DOI: 10.1103/Phys-
RevA.99.032331.
[2] Marcello Benedetti, Erika Lloyd, Stefan Sack, and
Mattia Fiorentini. Parameterized quantum cir-
cuits as machine learning models. Quantum Sci.
Technol., 4(4):043001, 2019. DOI: 10.1088/2058-
9565/ab4eb5.
[3] Vedran Dunjko, Jacob M Taylor, and Hans J
Briegel. Quantum-enhanced machine learning.
Phys. Rev. Lett., 117(13):130501, 2016. DOI:
10.1103/PhysRevLett.117.130501.
[4] Iordanis Kerenidis, Jonas Landman, Alessandro
Luongo, and Anupam Prakash. q-means: A quan-
tum algorithm for unsupervised machine learning.
In Advances in Neural Information Processing Sys-
tems, pages 4136–4146, 2019.
[5] Marco Cerezo, Alexander Poremba, Lukasz Cin-
cio, and Patrick J Coles. Variational quantum fi-
delity estimation. Quantum, 4:248, 2020. DOI:
10.22331/q-2020-03-26-248.
[6] Ryan LaRose, Arkin Tikku,
´
Etude O’Neel-Judy,
Lukasz Cincio, and Patrick J Coles. Variational
quantum state diagonalization. npj Quantum Inf.,
5(1):1–10, 2019. DOI: 10.1038/s41534-019-0167-6.
[7] Carlos Bravo-Prieto, Ryan LaRose, Marco Cerezo,
Yigit Subasi, Lukasz Cincio, and Patrick J Coles.
Variational quantum linear solver: A hybrid
algorithm for linear systems. arXiv preprint
arXiv:1909.05820, 2019.
[8] Juan Miguel Arrazola, Thomas R Bromley, Josh
Izaac, Casey R Myers, Kamil Br´adler, and Nathan
Killoran. Machine learning method for state prepa-
ration and gate synthesis on photonic quantum
computers. Quantum Sci. Technol., 4(2):024004,
2019. DOI: 10.1088/2058-9565/aaf59e.
[9] Nathan Killoran, Thomas R Bromley, Juan Miguel
Arrazola, Maria Schuld, Nicol´as Quesada, and Seth
Lloyd. Continuous-variable quantum neural net-
works. Phys. Rev. Research, 1(3):033063, 2019.
DOI: 10.1103/PhysRevResearch.1.033063.
[10] Nathan Killoran, Josh Izaac, Nicol´as Quesada,
Ville Bergholm, Matthew Amy, and Christian
Weedbrook. Strawberry Fields: A software plat-
form for photonic quantum computing. Quantum,
3:129, 2019. DOI: 10.22331/q-2019-03-11-129.
Accepted in Quantum 2020-11-05, click title to verify. Published under CC-BY 4.0. 18
[11] Gregory R Steinbrecher, Jonathan P Olson, Dirk
Englund, and Jacques Carolan. Quantum optical
neural networks. npj Quantum Inf., 5(1):1–9, 2019.
DOI: 10.1038/s41534-019-0174-7.
[12] N Quesada, LG Helt, J Izaac, JM Arrazola,
R Shahrokhshahi, CR Myers, and KK Sabapathy.
Simulating realistic non-Gaussian state prepara-
tion. Phys. Rev. A, 100(2):022341, 2019. DOI:
10.1103/PhysRevA.100.022341.
[13] Jonathan Romero, Jonathan P Olson, and Alan
Aspuru-Guzik. Quantum autoencoders for effi-
cient compression of quantum data. Quantum Sci.
Technol., 2(4):045001, 2017. DOI: 10.1088/2058-
9565/aa8072.
[14] Maria Schuld, Alex Bocharov, Krysta M Svore,
and Nathan Wiebe. Circuit-centric quantum clas-
sifiers. Phys. Rev. A, 101(3):032308, 2020. DOI:
10.1103/PhysRevA.101.032308.
[15] Maria Schuld and Nathan Killoran. Quantum ma-
chine learning in feature Hilbert spaces. Phys. Rev.
Lett., 122(4):040504, 2019. DOI: 10.1103/Phys-
RevLett.122.040504.
[16] Vojtˇech Havl´ıˇcek, Antonio D C´orcoles, Kristan
Temme, Aram W Harrow, Abhinav Kandala,
Jerry M Chow, and Jay M Gambetta. Su-
pervised learning with quantum-enhanced feature
spaces. Nature, 567(7747):209–212, 2019. DOI:
10.1038/s41586-019-0980-2.
[17] Peter JJ O’Malley, Ryan Babbush, Ian D
Kivlichan, Jonathan Romero, Jarrod R McClean,
Rami Barends, Julian Kelly, Pedram Roushan,
Andrew Tranter, Nan Ding, et al. Scalable
quantum simulation of molecular energies. Phys.
Rev. X, 6(3):031007, 2016. DOI: 10.1103/Phys-
RevX.6.031007.
[18] Alberto Peruzzo, Jarrod McClean, Peter Shadbolt,
Man-Hong Yung, Xiao-Qi Zhou, Peter J Love, Al´an
Aspuru-Guzik, and Jeremy L O’Brien. A vari-
ational eigenvalue solver on a photonic quantum
processor. Nature Commun., 5:4213, 2014. DOI:
10.1038/ncomms5213.
[19] Alberto Politi, Jonathan CF Matthews, Mark G
Thompson, and Jeremy L O’Brien. Integrated
quantum photonics. IEEE J. Sel. Top. Quan-
tum Electron., 15(6):1673–1684, 2009. DOI:
10.1109/JSTQE.2009.2026060.
[20] Y Zhang, M Menotti, K Tan, VD Vaidya,
DH Mahler, L Zatti, M Liscidini, B Morrison,
and Z Vernon. Single-mode quadrature squeez-
ing using dual-pump four-wave mixing in an in-
tegrated nanophotonic device. arXiv preprint
arXiv:2001.09474, 2020.
[21] VD Vaidya, B Morrison, LG Helt, R Shahrokshahi,
DH Mahler, MJ Collins, K Tan, J Lavoie, A Re-
pingon, M Menotti, et al. Broadband quadrature-
squeezed vacuum and nonclassical photon num-
ber correlations from a nanophotonic device. Sci.
Adv., 6(39):eaba9186, 2020. DOI: 10.1126/sci-
adv.aba9186.
[22] Seth Lloyd and Samuel L Braunstein. Quantum
computation over continuous variables. In Quan-
tum information with continuous variables, pages
9–17. Springer, 1999. DOI: 10.1007/978-94-015-
1258-9˙2.
[23] Alessio Serafini. Quantum continuous variables: a
primer of theoretical methods. CRC Press, 2017.
[24] Christian Weedbrook, Stefano Pirandola, Ra´ul
Garc´ıa-Patr´on, Nicolas J Cerf, Timothy C Ralph,
Jeffrey H Shapiro, and Seth Lloyd. Gaussian quan-
tum information. Rev. Mod. Phys., 84(2):621,
2012. DOI: 10.1103/RevModPhys.84.621.
[25] Stephen Barnett and Paul M Radmore. Methods
in theoretical quantum optics, volume 15. Oxford
University Press, 2002.
[26] Kevin E Cahill and Roy J Glauber. Ordered
expansions in boson amplitude operators. Phys.
Rev., 177(5):1857, 1969. DOI: 10.1103/Phys-
Rev.177.1857.
[27] P Kr´al. Displaced and squeezed Fock states.
J. Mod. Opt., 37(5):889–917, 1990. DOI:
10.1080/09500349014550941.
[28] Xin Ma and William Rhodes. Multimode squeeze
operators and squeezed states. Phys. Rev. A, 41
(9):4625, 1990. DOI: 10.1103/PhysRevA.41.4625.
[29] N. Quesada. Very Nonlinear Quantum Optics.
PhD thesis, University of Toronto, 2015.
[30] Ish Dhand, Barry C Sanders, and Hubert de Guise.
Algorithms for SU(n) boson realizations and D-
functions. J. Math. Phys., 56(11):111705, 2015.
DOI: 10.1063/1.4935433.
[31] EV Doktorov, IA Malkin, and VI Man’ko. Dy-
namical symmetry of vibronic transitions in poly-
atomic molecules and the Franck-Condon princi-
ple. J. Mol. Spectrosc, 64(2):302–326, 1977. DOI:
10.1016/0022-2852(75)90199-X.
[32] Daniel Gruner and Paul Brumer. Efficient evalua-
tion of harmonic polyatomic Franck-Condon fac-
tors. Chem. Phys. Lett., 138(4):310–314, 1987.
DOI: 10.1016/0009-2614(87)80389-5.
[33] R Berger, C Fischer, and M Klessinger. Calculation
of the vibronic fine structure in electronic spectra
at higher temperatures. 1. benzene and pyrazine.
J. Phys. Chem. A, 102(36):7157–7167, 1998. DOI:
10.1021/jp981597w.
[34] Vadim Mozhayskiy, Samer Gozem, and Anna I.
Krylov. ezspectrum v3.0. http://iopenshell.
usc.edu/downloads/, 2016.
Accepted in Quantum 2020-11-05, click title to verify. Published under CC-BY 4.0. 19
[35] Joonsuk Huh. Unified description of vibronic
transitions with coherent states. PhD thesis, Jo-
hann Wolfgang Goethe-Universit¨at in Frankfurt
am Main, 2011.
[36] Scott M Rabidoux, Victor Eijkhout, and John F
Stanton. A highly-efficient implementation of the
Doktorov recurrence equations for Franck–Condon
calculations. J. Chem. Theory Comput., 12(2):728–
739, 2016. DOI: 10.1021/acs.jctc.5b00560.
[37] Marcos Moshinsky and Christiane Quesne. Linear
canonical transformations and their unitary repre-
sentations. J. Math. Phys., 12(8):1772–1780, 1971.
DOI: 10.1063/1.1665805.
[38] Kurt Bernardo Wolf. Canonical transforms. i. com-
plex linear transforms. J. Math. Phys., 15(8):1295–
1301, 1974. DOI: 10.1063/1.1666811.
[39] P Kramer, Marcos Moshinsky, and TH Seligman.
Complex extensions of canonical transformations
and quantum mechanics. In Group theory and its
applications, pages 249–332. Elsevier, 1975. DOI:
10.1016/B978-0-12-455153-4.50011-3.
[40] Brajesh Gupt, Josh Izaac, and Nicol´as Quesada.
The Walrus: a library for the calculation of haf-
nians, Hermite polynomials and Gaussian boson
sampling. J. Open Source Softw., 4(44):1705, 2019.
DOI: 10.21105/joss.01705.
[41] Andreas Bj¨orklund, Brajesh Gupt, and Nicol´as
Quesada. A faster hafnian formula for complex ma-
trices and its benchmarking on a supercomputer.
ACM J. Exp. Algorithmics, 24(1):11, 2019. DOI:
10.1145/3325111.
[42] Nicol´as Quesada. Franck-Condon factors by
counting perfect matchings of graphs with loops.
J. Chem. Phys., 150(16):164113, 2019. DOI:
10.1063/1.5086387.
[43] Leonardo Banchi, Nicol´as Quesada, and
Juan Miguel Arrazola. Training Gaussian boson
sampling distributions. Phys. Rev. A, 102:012417,
2020. DOI: 10.1103/PhysRevA.102.012417.
[44] Claude Bloch and Albert Messiah. The canonical
form of an antisymmetric tensor and its application
to the theory of superconductivity. Nucl. Phys., 39:
95–106, 1962. DOI: 10.1016/0029-5582(62)90377-
2.
[45] Filippo M Miatto. Recursive multivariate deriva-
tives of e
f(x
1
,...,x
n
)
of arbitrary order. arXiv
preprint arXiv:1911.11722, 2019.
[46] Gianfranco Cariolaro and Gianfranco Pierobon.
Bloch-Messiah reduction of Gaussian unitaries by
Takagi factorization. Phys. Rev. A, 94(6):062109,
2016. DOI: 10.1103/PhysRevA.94.062109.
[47] Andrei B Klimov and Sergei M Chumakov.
A group-theoretical approach to quantum optics:
models of atom-field interactions. John Wiley &
Sons, 2009.
[48] Raphael Hunger. An introduction to complex dif-
ferentials and complex differentiability. Technical
report, Munich University of Technology, Inst. for
Circuit Theory and Signal Processing, 2007.
[49] Chu Guo and Dario Poletti. A scheme for au-
tomatic differentiation of complex loss functions.
arXiv preprint arXiv:2003.04295, 2020.
[50] Daniel Gottesman, Alexei Kitaev, and John
Preskill. Encoding a qubit in an oscillator. Phys.
Rev. A, 64(1):012310, 2001. DOI: 10.1103/Phys-
RevA.64.012310.
[51] Joonsuk Huh. Multimode Bogoliubov transfor-
mation and Husimi’s q-function. arXiv preprint
arXiv:2004.05766, 2020.
[52] St´efan van der Walt, S Chris Colbert, and Gael
Varoquaux. The NumPy array: a structure for ef-
ficient numerical computation. Comput. Sci. Eng.,
13(2):22–30, 2011. DOI: 10.1109/MCSE.2011.37.
[53] Pauli Virtanen, Ralf Gommers, Travis E Oliphant,
Matt Haberland, Tyler Reddy, David Courna-
peau, Evgeni Burovski, Pearu Peterson, Warren
Weckesser, Jonathan Bright, et al. SciPy 1.0:
fundamental algorithms for scientific computing in
Python. Nat. Methods, 17(3):261–272, 2020. DOI:
10.1038/s41592-019-0686-2.
[54] Thomas Kluyver, Benjamin Ragan-Kelley, Fer-
nando P´erez, Brian E Granger, Matthias Busson-
nier, Jonathan Frederic, Kyle Kelley, Jessica B
Hamrick, Jason Grout, Sylvain Corlay, et al.
Jupyter notebooks-a publishing format for repro-
ducible computational workflows. In ELPUB,
pages 87–90, 2016. DOI: 10.3233/978-1-61499-649-
1-87.
[55] John D Hunter. Matplotlib: A 2d graphics environ-
ment. Comput. Sci, Eng., 9(3):90–95, 2007. DOI:
10.1109/MCSE.2007.55.
[56] Siu Kwan Lam, Antoine Pitrou, and Stanley Seib-
ert. Numba: A LLVM-based Python JIT com-
piler. In Proceedings of the Second Workshop on
the LLVM Compiler Infrastructure in HPC, pages
1–6, 2015. DOI: 10.1145/2833157.2833162.
[57] S Berkowitz and FJ Garner. The calculation of
multidimensional Hermite polynomials and Gram-
Charlier coefficients. Math. Comput., 24(111):537–
545, 1970. DOI: 10.2307/2004829.
[58] Pieter Kok and Samuel L Braunstein. Multi-
dimensional Hermite polynomials in quantum op-
tics. J. Phys. A: Math. Gen., 34(31):6185, 2001.
DOI: 10.1088/0305-4470/34/31/312.
[59] Maurice M Mizrahi. Generalized Hermite polyno-
mials. J. Comput. Appl. Math., 1(3):137–140, 1975.
DOI: 10.1016/0771-050X(75)90031-5.
Accepted in Quantum 2020-11-05, click title to verify. Published under CC-BY 4.0. 20
[60] Samuel L Braunstein and Peter Van Loock.
Quantum information with continuous variables.
Rev. Mod. Phys., 77(2):513, 2005. DOI:
10.1103/RevModPhys.77.513.
[61] Timjan Kalajdzievski, Christian Weedbrook, and
Patrick Rebentrost. Continuous-variable gate de-
composition for the Bose-Hubbard model. Phys.
Rev. A, 97(6):062311, 2018. DOI: 10.1103/Phys-
RevA.97.062311.
[62] Timjan Kalajdzievski and Juan Miguel Arrazola.
Exact gate decompositions for photonic quantum
computing. Phys. Rev. A, 99(2):022341, 2019.
DOI: 10.1103/PhysRevA.99.022341.
[63] Ryotatsu Yanagimoto, Tatsuhiro Onodera, Edwin
Ng, Logan G Wright, Peter L McMahon, and
Hideo Mabuchi. Engineering a Kerr-based de-
terministic cubic phase gate via Gaussian opera-
tions. Phys. Rev. Lett., 124(24):240503, 2020. DOI:
10.1103/PhysRevLett.124.240503.
[64] Mile Gu, Christian Weedbrook, Nicolas C
Menicucci, Timothy C Ralph, and Peter van
Loock. Quantum computing with continuous-
variable clusters. Phys. Rev. A, 79(6):062318, 2009.
DOI: 10.1103/PhysRevA.79.062318.
[65] Kevin Marshall, Raphael Pooser, George Siopsis,
and Christian Weedbrook. Repeat-until-success
cubic phase gate for universal continuous-variable
quantum computation. Phys. Rev. A, 91(3):
032321, 2015. DOI: 10.1103/PhysRevA.91.032321.
[66] Krishna Kumar Sabapathy and Christian Weed-
brook. ON states as resource units for univer-
sal quantum computation with photonic architec-
tures. Phys. Rev. A, 97(6):062315, 2018. DOI:
10.1103/PhysRevA.97.062315.
[67] Krishna Kumar Sabapathy, Haoyu Qi, Josh Izaac,
and Christian Weedbrook. Production of pho-
tonic universal quantum gates enhanced by ma-
chine learning. Phys. Rev. A, 100(1):012326, 2019.
DOI: 10.1103/PhysRevA.100.012326.
[68] Milton Abramowitz and Irene A Stegun. Handbook
of mathematical functions with formulas, graphs,
and mathematical tables, volume 55. US Govern-
ment printing office, 1948.
[69] Vallee Olivier and Soares Manuel. Airy functions
and applications to physics. World Scientific Pub-
lishing Company, 2010.
[70] Kevin Marshall, Raphael Pooser, George Siop-
sis, and Christian Weedbrook. Quantum simu-
lation of quantum field theory using continuous
variables. Phys. Rev. A, 92:063825, 2015. DOI:
10.1103/PhysRevA.92.063825.
[71] Anthony Zee. Quantum field theory in a nutshell,
volume 7. Princeton university press, 2010.
[72] Laura Garc´ıa-
´
Alvarez, Cameron Calcluth,
Alessandro Ferraro, and Giulia Ferrini. Effi-
cient simulatability of continuous-variable circuits
with large Wigner negativity. arXiv preprint
arXiv:2005.12026, 2020.
[73] Izrail Solomonovich Gradshteyn and Iosif Moisee-
vich Ryzhik. Table of integrals, series, and prod-
ucts. Academic press, 1980.
Accepted in Quantum 2020-11-05, click title to verify. Published under CC-BY 4.0. 21