# Fast optimization of parametrized quantum optical circuits

Filippo M. Miatto  Institut Polytechnique de Paris  Télécom Paris, LTCI, 19
Place Marguerite Perey 91120 Palaiseau  Nicolás Quesada  Xanadu, Toronto, ON,
M5G 2C8, Canada

(November 5, 2020)

###### Abstract

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 optimizing these circuits is challenging due to
the infinite 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 cutoff. To tackle this issue, we present an algorithm
that is orders of magnitude faster than the current 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 computing
the matrix elements of Gaussian operators 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 automatic
differentiation framework such as TensorFlow or PyTorch. Our results will find
applications in quantum optical hardware research, quantum machine learning,
optical data processing, 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 information 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, 4, 5, 6, 7]. As a
parallel thread, the quantum optics community has investigated how to compound
quantum optical primitives to generate complex quantum states [8, 9, 10, 11,
12]. For this program to continue, it becomes necessary to develop new and
efficient tools to perform fast simulation of quantum optical circuits.

Parametrized quantum circuits are information processing devices whose
elements depend on one or more parameters, and therefore one can think of them
as a particular class of neural networks [13, 14, 15, 16, 17, 18]. Among
parametrized quantum circuits, quantum optical circuits are those in which the
carriers of quantum information are photons and the gates are optical
transformations [9, 10, 11]. Working in the optical domain offers several
advantages such as high processing speed, no demand for vacuum or cold
temperatures and the potential to build compact devices via integrated optics
[19, 19, 20, 21].

A typical architecture of quantum optical circuits alternates 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
symplectic 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 exactly 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 matrix multiplication and as such it is efficient. However, the
difficult part is to generate the transformation matrix from the parameters of
the gate: optical transformations (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 compute the
matrix exponential:

$\displaystyle\mathrm{trunc}(\exp(A))\neq\exp(\mathrm{trunc}(A)).$ (1)

For some special transformations such as the displacement [26] and squeezing
operator [27], exact formulas 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 truncated before multiplication, or
alternatively allows the matrix elements to be expressed as polynomials of the
parameters of the transformations [28, 25, 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 matrix
elements of a passive linear optical transformation corresponding to an
interferometer. In the context of Franck-Condon transitions, the chemical-
physics community [31, 32, 33, 34, 35, 36] has studied the matrix elements of
Gaussian transformations where the displacement, squeezing parameters and
unitary matrices representing interferometers are real, corresponding to so-
called point transformations where canonical positions and momenta of the
harmonic oscillator do not mix [37, 38, 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, accepting that
truncation errors could propagate.

Our method to perform the optimization of a parametrized quantum optical
circuit has several advantages over previous techniques:

  1. 1.

It generates recursively the matrix elements of general Gaussian
transformations exactly up to the desired cutoff dimension without the need
for decompositions. These results are derived in Sec. 2 for general Gaussian
gates and applied to specific cases in Sec. 3.

  2. 2.

It constructs gradients exactly, without the need for computational graphs
describing the functional dependencies between the variables. This is
explained in Sec. 4.

  3. 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 experiments.

Even though we focus on Gaussian gates in the main text, in Appendices C, D, E
and F we derive the analytical form of the matrix elements of the cubic- and
quartic- phase gates together with recurrence relations for their calculation
and the calculation of their gradients. 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 instead one were interested in just
a few specific matrix elements, one could map this to the calculation of loop
hafnians [41, 42, 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:

$\displaystyle\bm{\beta}^{\bm{n}}$
$\displaystyle\equiv\prod_{i=1}^{\ell}\beta_{i}^{n_{i}},$ (2)
$\displaystyle\bm{n}!$ $\displaystyle\equiv\prod_{i=1}^{\ell}n_{i}!.$ (3)

We indicate tensor products of $\ell$-mode coherent and Fock states as
follows:

$\displaystyle|\bm{\beta}\rangle$
$\displaystyle=|\beta_{1}\rangle\otimes\ldots\otimes|\beta_{\ell}\rangle,$ (4)
$\displaystyle|\bm{n}\rangle$
$\displaystyle=|n_{1}\rangle\otimes\ldots\otimes|n_{\ell}\rangle,$ (5)

which we can use for instance to express a multimode coherent state:

$\displaystyle|\bm{\beta}\rangle$
$\displaystyle=e^{-\frac{1}{2}||\bm{\beta}||^{2}}\sum_{\bm{n}=\bm{0}}^{\bm{\infty}}\frac{\bm{\beta}^{\bm{n}}}{\sqrt{\bm{n}!}}|\bm{n}\rangle,$
(6)

where we used the shorthand notation
$\sum_{\bm{n}=\bm{0}}^{\bm{\infty}}=\sum_{n_{1}=0}^{\infty}\ldots\sum_{n_{\ell}=0}^{\infty}$.
Note that tensor indices are numbered 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:

$\displaystyle\frac{\partial^{\bm{k}}}{\partial{\bm{\beta}}^{\bm{k}}}=\partial_{\bm{\beta}}^{\bm{k}}\equiv\prod_{i}\frac{\partial^{k_{i}}}{\partial{\beta_{i}}^{k_{i}}}.$
(7)

The key insight of the generating function method is that when we compute the
inner product $\langle\bm{\alpha}^{*}|\mathcal{G}|\bm{\beta}\rangle$ of any
operator $\mathcal{G}$ (not just a Gaussian one) with a pair of coherent
states, we obtain a function of the coherent amplitudes that can be treated as
a generating function for the matrix elements of $\mathcal{G}$:

$\displaystyle\Gamma(\bm{\alpha},\bm{\beta})$
$\displaystyle=e^{\frac{1}{2}(||\bm{\alpha}||^{2}+||\bm{\beta}||^{2})}\langle\bm{\alpha}^{*}|\mathcal{G}|\bm{\beta}\rangle$
(8)
$\displaystyle=\sum_{\bm{m},\bm{n}=\bm{0}}^{\bm{\infty}}\frac{\bm{\alpha}^{\bm{m}}\bm{\beta}^{\bm{n}}}{\sqrt{\bm{m}!\bm{n}!}}\langle\bm{m}|\mathcal{G}|\bm{n}\rangle.$
(9)

The matrix element $\langle\bm{m}|\mathcal{G}|\bm{n}\rangle$ is therefore
contained in the coefficient of the factor
$\bm{\alpha}^{\bm{m}}\bm{\beta}^{\bm{n}}$ in the Taylor series and we can
isolate it by computing the derivatives of the appropriate order, at zero:

$\displaystyle\mathcal{G}_{\bm{m}\bm{n}}\equiv\langle\bm{m}|\mathcal{G}|\bm{n}\rangle=\frac{\partial^{\bm{m}}_{\bm{\alpha}}\partial^{\bm{n}}_{\bm{\beta}}\Gamma(\bm{\alpha},\bm{\beta})}{\sqrt{\bm{m}!\bm{n}!}}\bigg{|}_{\bm{\alpha}=\bm{\beta}=0}.$
(10)

We will derive recursion relations for $\mathcal{G}_{\bm{mn}}$ directly,
instead 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

$\displaystyle\bm{\nu}$ $\displaystyle=\left[\begin{array}[]{c}\bm{\alpha}\\\
\bm{\beta}\end{array}\right]\in\mathbb{C}^{2\ell},$ (13) $\displaystyle\bm{k}$
$\displaystyle=\left[\begin{array}[]{c}\bm{m}\\\
\bm{n}\end{array}\right]\in\mathbb{N}_{0}^{2\ell},$ (16)

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

$\displaystyle\Gamma(\bm{\nu})=\sum_{\bm{k}=\bm{0}}^{\bm{\infty}}\frac{\bm{\nu}^{\bm{k}}}{\sqrt{\bm{k}!}}\mathcal{G}_{\bm{k}}.$
(17)

where we take advantage of the definition of the vector $\bm{k}$ in Eq. (13)
and write $\mathcal{G}_{\bm{k}}$ for $\mathcal{G}_{\bm{m}\bm{n}}$. As shown in
Appendix A, the generating function $\Gamma(\bm{\nu})$ for a Gaussian operator
$\mathcal{G}$ is in exponential form $\Gamma(\bm{\nu})=C\,e^{Q(\bm{\nu})}$.
This implies that we can express high-order derivatives recursively, using
only the derivatives at zero of the exponent $Q(\bm{\nu})$, as shown in
subsection 2.3. Furthermore, the exponent $Q(\bm{\nu})$ is a polynomial of
degree 2, which we can write as:

$\displaystyle
Q(\bm{\nu})=\bm{\mu}^{T}\bm{\nu}-\frac{1}{2}\bm{\nu}^{T}\bm{\Sigma}\bm{\nu}.$
(18)

This means that all of its mixed derivatives of order higher than 2 vanish:

$\displaystyle\partial^{\bm{k}}_{\bm{\nu}}Q(\bm{\nu})|_{\bm{\nu}=\bm{0}}=0\text{
if }\sum_{i=1}^{2\ell}k_{i}>2.$ (19)

Because of this, all the derivatives of the exponent that we will ever need
are those of degree 1 and 2:

$\displaystyle\frac{\partial}{\partial\nu_{i}}Q(\bm{\nu})|_{\bm{\nu}=\bm{0}}$
$\displaystyle=\mu_{i},$ (20)
$\displaystyle\frac{\partial^{2}}{\partial\nu_{i}\partial\nu_{j}}Q(\bm{\nu})|_{\bm{\nu}=\bm{0}}$
$\displaystyle=-\Sigma_{ij}.$ (21)

Finally, note that although $\mathcal{G}_{\bm{k}}$ is a tensor of rank
$2\ell$, 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 [44, 23], a general
$\ell-$mode Gaussian transformation can be parametrized as

$\displaystyle\mathcal{G}=\mathcal{G}(\bm{\gamma},\bm{W},\bm{\zeta},\bm{V})=\mathcal{D}(\bm{\gamma})\mathcal{U}(\bm{W})\mathcal{S}(\bm{\zeta})\mathcal{U}(\bm{V}),$
(22)

where we used the notation
$\mathcal{D}(\bm{\gamma})=\bigotimes_{i=1}^{\ell}D_{i}(\gamma_{i})$,
$\mathcal{S}(\bm{\zeta})=\bigotimes_{i=1}^{\ell}S_{i}(\zeta_{i})$ for tensor
products of single-mode displacement and squeezing operators defined by

$\displaystyle D_{j}(\gamma_{j})$
$\displaystyle=\exp\left(\gamma_{j}a^{\dagger}_{j}-\text{H.c.}\right),\quad\gamma_{j}\in\mathbb{C},$
(23) $\displaystyle S_{j}(\zeta_{j})$
$\displaystyle=\exp\left(\tfrac{\zeta_{j}^{*}}{2}a_{j}^{2}-\text{H.c.}\right),\quad\zeta_{j}=r_{j}e^{i\delta_{j}}\in\mathbb{C},$
(24)

where $r_{j}$ and $\delta_{j}$ are the phase and modulus of the complex
squeezing parameter $\zeta_{j}$.

In the last equation we used the annihilation $a_{j}$ and creation
$a_{j}^{\dagger}$ operators of the $\ell$ modes satisfying the canonical
commutation relations

$\displaystyle[a_{j},a_{l}^{\dagger}]=\delta_{j,l},\quad[a_{j},a_{l}]=[a^{\dagger}_{j},a^{\dagger}_{l}]=0.$
(25)

The Hilbert space operator $\mathcal{U}(\bm{V})$ represents a general
multimode passive transformation (physically corresponding to an
interferometer) and is parametrized by a unitary matrix $\bm{V}$ that
transforms the creation operators as

$\displaystyle\mathcal{U}^{\dagger}(\bm{V})a^{\dagger}_{i}\mathcal{U}(\bm{V})=\sum_{l=1}^{\ell}V^{*}_{il}a_{l}^{\dagger},\quad\mathcal{U}(\bm{V})a^{\dagger}_{l}\mathcal{U}^{\dagger}(\bm{V})=\sum_{i=1}^{\ell}V_{il}a_{i}^{\dagger}.$
(26)

Using these definitions, in Appendix A we show the following:

$\displaystyle
e^{\frac{1}{2}\left[||\bm{\alpha}||^{2}+||\bm{\beta}||^{2}\right]}\langle\bm{\alpha}^{*}|\mathcal{D}(\bm{\gamma})\mathcal{U}(\bm{W})\mathcal{S}(\bm{\zeta})\mathcal{U}(\bm{V})|\bm{\beta}\rangle=$
$\displaystyle\quad
C\exp\left(\bm{\mu}^{T}\bm{\nu}-\frac{1}{2}\bm{\nu}^{T}\bm{\Sigma}\bm{\nu}\right),$
(27)

where

$\displaystyle
C=\frac{\exp\left(-\tfrac{1}{2}\left[||\bm{\gamma}||^{2}+\bm{\gamma}^{\dagger}\bm{W}\text{diag}(e^{i\bm{\delta}}\tanh\bm{r})\bm{W}^{T}\bm{\gamma}^{*}\right]\right)}{\sqrt{\prod_{i=1}^{\ell}\cosh
r_{i}}},$ (28) $\displaystyle\bm{\mu}^{T}=$ (29)
$\displaystyle\left[\bm{\gamma}^{\dagger}\bm{W}\text{diag}(e^{i\bm{\delta}}\tanh\bm{r})\bm{W}^{T}+\bm{\gamma}^{T},-\bm{\gamma}^{\dagger}\bm{W}\text{diag}(\operatorname{sech}\bm{r})\bm{V}\right],$
$\displaystyle\bm{\Sigma}=$ (30)
$\displaystyle\left[\begin{array}[]{c|c}\bm{W}\text{diag}(e^{i\bm{\delta}}\tanh\bm{r})\bm{W}^{T}&-\bm{W}\text{diag}(\operatorname{sech}\bm{r})\bm{V}\\\
\hline\cr-\bm{V}^{T}\text{diag}(\operatorname{sech}\bm{r})\bm{W^{T}}&-\bm{V}^{T}\text{diag}(e^{-i\bm{\delta}}\tanh\bm{r})\bm{V}\end{array}\right].$
(33)

Note that $\bm{\mu}$ is a vector of dimension $2\ell$. The matrix $\Sigma$, on
the other hand is built out of four $\ell\times\ell$ blocks and has therefore
dimension $2\ell\times 2\ell$.

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 elements

We use the recent result in [45] to express the matrix elements of the
Gaussian operator recursively (see Appendix B):

$\displaystyle\mathcal{G}_{\bm{0}}$ $\displaystyle=C,$ (34)
$\displaystyle\mathcal{G}_{\bm{k}+1_{i}}$
$\displaystyle=\frac{1}{\sqrt{k_{i}+1}}\sum_{\bm{j}=\bm{0}}^{\bm{k}}\sqrt{\frac{\bm{k}!}{\bm{j}!}}\frac{1}{(\bm{k}-\bm{j})!}\mathcal{G}_{\bm{j}}\partial_{\bm{\nu}}^{\bm{k}+1_{i}-\bm{j}}Q(\bm{\nu})|_{\bm{\nu}=\bm{0}},$
(35)

where $1_{i}$ is a vector with all zeros and a single 1 at position $i$, and
therefore to obtain the vector $\bm{k}+1_{i}$ we take the vector $\bm{k}$ and
we increment $k_{i}$ by 1.

We stress that Eq. (35) 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 considerably
because the derivative of $Q$ is zero unless $\bm{j}=\bm{k}$ or
$\bm{j}=\bm{k}-1_{l}$ for some index $l$:

$\displaystyle\mathcal{G}_{\bm{k}+1_{i}}=\frac{1}{\sqrt{k_{i}+1}}\left(\mathcal{G}_{\bm{k}}\mu_{i}-\sum_{l=1}^{2\ell}\sqrt{k_{l}}\mathcal{G}_{\bm{k}-1_{l}}\Sigma_{il}\right).$
(36)

This, together with Eq. (118) which shows how to obtain gradients, is the
primary result of our work. As an example of Eq. (36), we derive the two
recurrence relations for the single mode case, i.e., $\ell=1$. For convenience
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:

$\displaystyle\mathcal{G}_{m+1,n}=$
$\displaystyle\frac{1}{\sqrt{m+1}}\biggl{(}\mathcal{G}_{m,n}\mu_{1}-\sqrt{m}\mathcal{G}_{m-1,n}\Sigma_{11}$
$\displaystyle-\sqrt{n}\mathcal{G}_{m,n-1}\Sigma_{12}\biggr{)},$ (37)
$\displaystyle\mathcal{G}_{m,n+1}=$
$\displaystyle\frac{1}{\sqrt{n+1}}\biggl{(}\mathcal{G}_{m,n}\mu_{2}-\sqrt{m}\mathcal{G}_{m-1,n}\Sigma_{21}$
$\displaystyle-\sqrt{n}\mathcal{G}_{m,n-1}\Sigma_{22}\biggr{)}.$ (38)

We represent this example pictorially in Fig. 1.

For Gaussian operators on $\ell$ modes, the recurrence relations allow us to
combine $2\ell+1$ neighbouring elements in $2\ell$ different ways (by
incrementing each of the $2\ell$ indices) and generate $2\ell$ new elements.

Note that these recurrence relations can be rewritten in an alternative form
by using the mode operators:

$\displaystyle a_{i}\mathcal{G}$
$\displaystyle=\mu_{i}\mathcal{G}-\left(\sum_{l=1}^{\ell}a_{l}^{\dagger}\Sigma_{il}\right)\mathcal{G}-\mathcal{G}\left(\sum_{l=1}^{\ell}a_{l}\Sigma_{i,l+\ell}\right),$
(39) $\displaystyle\mathcal{G}a_{i}^{\dagger}$
$\displaystyle=\mu_{i+\ell}\mathcal{G}-\left(\sum_{l=1}^{\ell}a_{l}^{\dagger}\Sigma_{i+\ell,l}\right)\mathcal{G}-\mathcal{G}\left(\sum_{l=1}^{\ell}a_{l}\Sigma_{i+\ell,l+\ell}\right).$
(40)

With these compact relations we can formulate the commutation relations
between a general multimode Gaussian operator and the mode operators $a_{i}$
and $a_{i}^{\dagger}$, for instance for a single-mode Gaussian operator, we
have:

$\displaystyle[\mathcal{G},a]$
$\displaystyle=-\mu_{1}\mathcal{G}+(\Sigma_{12}+1)\mathcal{G}a+\Sigma_{11}a^{\dagger}\mathcal{G},$
(41) $\displaystyle[\mathcal{G},a^{\dagger}]$
$\displaystyle=\mu_{2}\mathcal{G}-(\Sigma_{21}+1)a^{\dagger}\mathcal{G}-\Sigma_{22}\mathcal{G}a.$
(42)

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

$\displaystyle[S(\zeta),a]$
$\displaystyle=(1-\mathrm{sech}(r))S(\zeta)a+e^{i\delta}\tanh(r)a^{\dagger}S(\zeta),$
(43) $\displaystyle[S(\zeta),a^{\dagger}]$
$\displaystyle=(\mathrm{sech}(r)-1)a^{\dagger}S(\zeta)+e^{-i\delta}\tanh(r)S(\zeta)a,$
(44) $\displaystyle[D(\gamma),a]$ $\displaystyle=-\gamma D(\gamma),$ (45)
$\displaystyle[D(\gamma),a^{\dagger}]$ $\displaystyle=-\gamma^{*}D(\gamma).$
(46)

##  3 Applications: single- and two-mode gates and passive Gaussian
transformations

All Gaussian single- and two-mode transformations will be expressed as a
combination of four elementary operations: single-mode phase shifts, single-
mode squeezing, single-mode displacement and beamsplitters. These elements are
parametrized as follows: the phase rotation gate $R(\phi)$ depends on an angle
$\phi$, the single-mode squeezer $S(\zeta)$ depends on a single complex
parameter which we express in polar form $\zeta=re^{i\delta}$, the
displacement gate $D(\gamma)$ depends on a single complex parameter $\gamma$
which we leave as is, and finally the beamsplitter $B(\theta,\varphi)$ depends
on two angles $(\theta,\varphi)$. Each gate is expressed as the exponential of
mode operators. In particular, for the rotation and beamsplitter gates we
write

$\displaystyle R(\phi)$ $\displaystyle=\exp\left(i\phi\,a^{\dagger}a\right),$
(47) $\displaystyle B(\theta,\varphi)$
$\displaystyle=\exp\left[\theta\left(e^{i\varphi}a_{1}a_{2}^{\dagger}-e^{-i\varphi}a_{1}^{\dagger}a_{2}\right)\right].$
(48)

For example, the 50/50 beam splitter with real coefficients corresponds to
$(\theta,\varphi)=(\pi/4,0)$. Note that the last two gates are special cases
of a passive linear optical transformation, thus they can be written as
$\mathcal{U}(\bm{W})$ where

$\displaystyle\bm{W}$ $\displaystyle=e^{i\phi},$ (49a) $\displaystyle\bm{W}$
$\displaystyle=\begin{bmatrix}\cos\theta&-e^{-i\varphi}\sin\theta\\\
e^{i\varphi}\sin\theta&\cos\theta\end{bmatrix},$ (49b)

respectively.

###  3.1 Single-mode Gaussian gate

A simple way to parametrize a general single-mode Gaussian operation is as
follows:

$\displaystyle\mathcal{G}^{(1)}=\mathcal{G}^{(1)}(\gamma,\phi,\zeta)=D(\gamma)R(\phi)S(\zeta).$
(50)

By using Eq. (49a) and the general result from the last section we easily find

$\displaystyle C$
$\displaystyle=\frac{\exp\left(-\tfrac{1}{2}\left[|\gamma|^{2}+\gamma^{*2}e^{i(\delta+2\phi)}\tanh
r\right]\right)}{\sqrt{\cosh r}},$ (51) $\displaystyle\bm{\mu}^{T}$
$\displaystyle=[\gamma^{*}e^{i(\delta+2\phi)}\tanh
r+\gamma,-\gamma^{*}e^{i\phi}\operatorname{sech}r],$ (52)
$\displaystyle\bm{\Sigma}$
$\displaystyle=\left[\begin{array}[]{c|c}e^{i(\delta+2\phi)}\tanh
r&-e^{i\phi}\operatorname{sech}r\\\
\hline\cr-e^{i\phi}\operatorname{sech}r&-e^{-i\delta}\tanh r\\\
\end{array}\right].$ (55)

These values enter the recurrence relations (37) and (38), which are depicted
in Fig. 1.

Figure 1: The three neighbouring matrix elements shown in light purple
($\mathcal{G}_{m-1,n}$, $\mathcal{G}_{m,n}$ and $\mathcal{G}_{m,n-1}$) can be
linearly combined in two different ways to generate $\mathcal{G}_{m+1,n}$ or
$\mathcal{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 $\ell$ modes, the recurrence rules allow us to combine $2\ell+1$
neighbouring elements in $2\ell$ different ways and generate $2\ell$ new
elements.

####  3.1.1 Special case: single-mode squeezing

For a single-mode squeezing operator, we have $\phi=0$ and $\gamma=0$, and
thus

$\displaystyle C$ $\displaystyle=\sqrt{\operatorname{sech}r},$ (56)
$\displaystyle\bm{\mu}^{T}$ $\displaystyle=[0,0],$ (57)
$\displaystyle\bm{\Sigma}$
$\displaystyle=\left[\begin{array}[]{c|c}e^{i\delta}\tanh
r&-\operatorname{sech}r\\\ \hline\cr-\operatorname{sech}r&-e^{-i\delta}\tanh
r\\\ \end{array}\right].$ (60)

The recurrence relations simplify accordingly, and we obtain the matrix
elements of the single-mode squeezer $S_{m,n}=\langle m|S(\zeta)|n\rangle$:

$\displaystyle S_{0,0}=$ $\displaystyle\sqrt{\operatorname{sech}r},$ (61)
$\displaystyle S_{m+1,0}=$
$\displaystyle-\sqrt{\frac{m}{m+1}}S_{m-1,0}e^{i\delta}\tanh r,$ (62)
$\displaystyle S_{m,n+1}=$
$\displaystyle\frac{1}{\sqrt{n+1}}\bigl{(}\sqrt{m}S_{m-1,n}\operatorname{sech}r+$
$\displaystyle\sqrt{n}S_{m,n-1}e^{-i\delta}\tanh r\bigr{)},$ (63)

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 $\bm{\mu}$ 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.

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.

####  3.1.2 Special case: displacement

For a single-mode displacement operator, we have $\phi=0$ and $\zeta=0$ and
therefore $\tanh r=0$ and $\operatorname{sech}r=1$ and we find

$\displaystyle C$ $\displaystyle=e^{-\frac{1}{2}|\gamma|^{2}},$ (64)
$\displaystyle\bm{\mu}^{T}$ $\displaystyle=[\gamma,-\gamma^{*}],$ (65)
$\displaystyle\bm{\Sigma}$ $\displaystyle=\left[\begin{array}[]{c|c}0&-1\\\
\hline\cr-1&0\\\ \end{array}\right].$ (68)

The recurrence relations simplify accordingly, and we obtain the matrix
elements of the single-mode displacement operator $D_{m,n}=\langle
m|D(\gamma)|n\rangle$:

$\displaystyle D_{0,0}$ $\displaystyle=e^{-\frac{|\gamma|^{2}}{2}},$ (69)
$\displaystyle D_{m+1,0}$ $\displaystyle=\frac{\gamma}{\sqrt{m+1}}D_{m,0},$
(70) $\displaystyle D_{m,n+1}$
$\displaystyle=-\frac{\gamma^{*}}{\sqrt{n+1}}D_{m,n}+\sqrt{\frac{m}{n+1}}D_{m-1,n}.$
(71)

Thanks to the diagonal of $\bm{\Sigma}$ being zero, also in this case each
step is simplified, as depicted in Fig. 3.

Figure 3: When generating the single-mode displacement matrix, 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

$\displaystyle\mathcal{G}^{(2)}$
$\displaystyle=\mathcal{D}(\bm{\gamma})\mathcal{R}(\bm{\phi})B(\theta^{\prime},\varphi^{\prime})\mathcal{S}(\bm{\zeta})B(\theta,\varphi),$
(72)

where $\bm{\phi}=[\phi_{1},\phi_{2}]$, $\bm{\gamma}=[\gamma_{1},\gamma_{2}]$
and $\bm{\zeta}=[\zeta_{1},\zeta_{2}]$ and where the single-mode gates
factorize over the two modes. This gives a total of 14 real parameters
(counting double for complex parameters), in accordance with Ref. [46]. The
equations above can be recast into the form of Eq. (22) by writing
$\mathcal{U}(\bm{V})=B(\theta,\varphi)$ and
$\mathcal{U}(\bm{W})=\mathcal{R}(\bm{\phi})B(\theta^{\prime},\varphi^{\prime})$,
with

$\displaystyle\bm{V}$
$\displaystyle=\left[\begin{array}[]{cc}\cos\theta&-e^{-i\varphi}\sin\theta\\\
e^{i\varphi}\sin\theta&\cos\theta\\\ \end{array}\right],$ (75)
$\displaystyle\bm{W}$
$\displaystyle=\left[\begin{array}[]{cc}e^{i\phi_{1}}&0\\\ 0&e^{i\phi_{2}}\\\
\end{array}\right]\left[\begin{array}[]{cc}\cos\theta^{\prime}&-e^{-i\varphi^{\prime}}\sin\theta^{\prime}\\\
e^{i\varphi^{\prime}}\sin\theta^{\prime}&\cos\theta^{\prime}\\\
\end{array}\right].$ (80)

####  3.2.1 Special case: Two-mode squeezer

The two-mode squeezer is defined as $S^{(2)}(\zeta)=\exp\left(\zeta
a^{\dagger}_{1}a^{\dagger}_{2}-\text{H.c.}\right)$ and can be decomposed in
terms of beamsplitters and single-mode squeezing operations as follows:

$\displaystyle S^{(2)}(\zeta)=B(-\pi/4,0)\left[S(\zeta)\otimes
S(-\zeta)\right]B(\pi/4,0).$ (81)

From the decomposition above one easily finds $C=\operatorname{sech}r$,
$\bm{\mu}^{T}=[0,0,0,0]$ and

$\displaystyle\bm{\Sigma}=(-1)\times$ (82)
$\displaystyle\left[\begin{array}[]{cccc}0&e^{i\delta}\tanh
r&\operatorname{sech}r&0\\\ e^{i\delta}\tanh r&0&0&\operatorname{sech}r\\\
\operatorname{sech}r&0&0&-e^{-i\delta}\tanh r\\\
0&\operatorname{sech}r&-e^{-i\delta}\tanh r&0\\\ \end{array}\right].$ (87)

Note that a two-mode squeezing operation creates photons in pairs, which
implies that the difference of photon number between the two modes is
conserved. This observation 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}=\langle m,n|S_{2}(\zeta)|p,q\rangle$ 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:

$\displaystyle S^{(2)}_{0,0,0,0}=$ $\displaystyle\operatorname{sech}r,$ (88)
$\displaystyle S^{(2)}_{n,n,0,0}=$
$\displaystyle-S^{(2)}_{n-1,n-1,0,0}\Sigma_{2,1},$ (89) $\displaystyle
S^{(2)}_{m,n,m-n,0}=$
$\displaystyle-\sqrt{\frac{m}{m-n}}S^{(2)}_{m-1,n,m-n-1,0}\Sigma_{3,1},$ (90)
$\displaystyle S^{(2)}_{m,n,p,p-(m-n)}=$
$\displaystyle\frac{-1}{\sqrt{p-(m-n)}}\times$ (91) $\displaystyle\
\biggl{(}\sqrt{n}S^{(2)}_{m,n-1,p,p-(m-n)-1}\Sigma_{4,2}$
$\displaystyle+\sqrt{p}S^{(2)}_{m,n,p-1,p-(m-n)-1}\Sigma_{4,3}\biggr{)}.$

###  3.3 Passive Gaussian transformations: Interferometers

These $\ell$-mode transformations, which correspond to interferometers
parametrized by an $\ell\times\ell$ unitary matrix $\bm{V}$, can be obtained
from our general results by setting $\bm{\zeta}=\bm{\delta}=\bm{0}$ and
selecting, without loss of generality, $\bm{W}=\mathbb{I}$ to obtain

$\displaystyle C$ $\displaystyle=1,$ (92) $\displaystyle\bm{\mu}^{T}$
$\displaystyle=\bm{0},$ (93) $\displaystyle\bm{\Sigma}$
$\displaystyle=-\left[\begin{array}[]{c|c}\bm{0}&\bm{V}\\\
\hline\cr\bm{V}^{T}&\bm{0}\end{array}\right].$ (96)

We can write
$\langle\bm{m}|\mathcal{U}(\bm{V})|\bm{\bm{n}}\rangle=U_{\bm{m},\bm{n}}$ and
find the following recurrence relations

$\displaystyle{U}_{\bm{0},\bm{0}}$ $\displaystyle=1,$ (97a)
$\displaystyle{U}_{\bm{m}+1_{i},\bm{n}}$
$\displaystyle=\frac{1}{\sqrt{m_{i}+1}}\sum_{j=1}^{\ell}V_{i,j}\sqrt{n_{j}}\
{U}_{\bm{m},\bm{n}-1_{j}},$ (97b) $\displaystyle{U}_{\bm{m},\bm{n}+1_{i}}$
$\displaystyle=\frac{1}{\sqrt{n_{i}+1}}\sum_{j=1}^{\ell}V_{j,i}\sqrt{m_{j}}\
{U}_{\bm{m}-1_{j},\bm{n}}.$ (97c)

These equations have a particular simple physical interpretation: they tell us
that to scatter off one additional photon in the bra side $\bm{m}+1_{i}$ one
photon has to be sent into the interferometer in the ket side $\bm{n}-1_{j}$.
In particular, we find that the probability amplitude that a single photon
injected in port $j$ ends up in port $i$ is precisely the $(i,j)$ entry of the
unitary matrix describing the interferometer:

$\displaystyle{U}_{1_{i},1_{j}}=V_{i,j}.$ (98)

Note that ${U}_{\bm{m},\bm{n}}$ is nonzero only if
$\sum_{i=1}^{\ell}m_{i}=\sum_{i=1}^{\ell}n_{i}$ since interferometers do not
create or destroy particles (mathematically, this is because they have $U(n)$
symmetry). Matrix elements that do not satisfy this selection rule are zero
and we can construct the whole rank-$2\ell$ tensor by looping over at most
$2\ell-1$ indices. Let us illustrate this for a beamsplitter, which is a
special case of a general passive transformation, and for which $\bm{V}$ is of
size $2\times 2$ and given in Eq. (75). The recurrence relations we find are:

$\displaystyle B_{0,0,0,0}=$ $\displaystyle 1,$ (99) $\displaystyle
B_{m,n,m+n,0}=$
$\displaystyle\frac{1}{\sqrt{m+n}}\biggl{(}\sqrt{m}V_{1,1}B_{m-1,n,m+n-1,0}$
$\displaystyle\quad\quad+\sqrt{n}V_{2,1}B_{m,n-1,m+n-1,0}\biggr{)},$ (100)
$\displaystyle B_{m,n,p,m+n-p}=$ $\displaystyle\frac{1}{\sqrt{m+n-p}}\times$
(101) $\displaystyle\biggl{(}\sqrt{m}V_{1,2}B_{m-1,n,p,m+n-p-1}$
$\displaystyle+\sqrt{n}V_{2,2}B_{m,n-1,p,m+n-p-1}\biggr{)},$

where $\langle m,n|B(\theta,\phi)|p,q\rangle=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

$\displaystyle|\psi\rangle=\sum_{k,l=0}^{N-1}c_{k,l}|k,l\rangle,$ (102)

then the state after the application of a beamsplitter is given by

$\displaystyle|\psi^{\prime}\rangle=$ $\displaystyle
B(\theta,\phi)|\psi\rangle=\sum_{m,n=0}^{N-1}c^{\prime}_{n,m}|n,m\rangle,$
(103)

and because of the $U(2)$ symmetry of the beamsplitter one can obtain each
entry of the updated tensor $c^{\prime}$ using a single sum:

$\displaystyle
c^{\prime}_{n,m}=\sum\limits_{k=\text{max}(1+n+m-N,0)}^{\text{min}(n+m,N-1)+1}B_{n,m,k,n+m-k}\
c_{k,n+m-k}.$ (104)

Similarly, one can use the $SU(1,1)$ symmetry of the two-mode squeezer to show
that if
$|\psi^{\prime\prime}\rangle=\sum_{nm=0}^{N-1}c^{\prime\prime}_{n,m}|n,m\rangle=S^{(2)}(z)|\psi\rangle$
then the coefficients can again be updated in a single sum as follows:

$\displaystyle
c^{\prime\prime}_{n,m}=\sum\limits_{k=\text{max}(n-m,0)}^{\text{min}(n-m,0)+N}S^{(2)}_{n,m,k,k+m-n}\
c_{k,k+m-n}.$ (105)

These computational savings are doubled when considering 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.

##  4 Recursive gradients from generating functions

In this section we present recurrence relations to construct the gradients of
Gaussian operators with respect to their parameters. Gradients are essential
in optimization 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 optimal choice of parameters corresponds to the global minimum 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 backpropagate 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 values 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 parameters

For a complex parameter $\xi$, the gradient descent update step should use the
partial derivative of a real loss function $L$ with respect to the _conjugate_
of the parameter [48, 49]:

$\displaystyle\xi\leftarrow\xi-\gamma\frac{\partial L}{\partial\xi^{*}}.$
(106)

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:

$\displaystyle\frac{\partial L}{\partial\xi^{*}}$
$\displaystyle=\sum_{\bm{k}}\frac{\partial
L}{\partial\mathcal{G}_{\bm{k}}^{*}}\frac{\partial\mathcal{G}_{\bm{k}}^{*}}{\partial\xi^{*}}+\frac{\partial
L}{\partial\mathcal{G}_{\bm{k}}}\frac{\partial\mathcal{G}_{\bm{k}}}{\partial\xi^{*}}.$
(107)

In an automatic differentiation framework such as TensorFlow or PyTorch, if we
wish to customize the computation 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 upstream gradient tensor $\partial
L/\partial\mathcal{G}_{\bm{k}}^{*}$. It is our task to combine it with the
local gradients $\partial\mathcal{G}_{\bm{k}}/\partial\xi^{*}$ and
$\partial\mathcal{G}_{\bm{k}}^{*}/\partial\xi^{*}$ as prescribed by the chain
rule, to produce the downstream gradient $\partial L/\partial\xi^{*}$.

In order to obtain the various parts of Eq. (107) we proceed as follows: for
the first part of (107), the upstream gradient tensor $\partial
L/\partial\mathcal{G}_{\bm{k}}^{*}$ is given to us by the software so we only
have to compute the local gradient:

$\displaystyle\frac{\partial\mathcal{G}_{\bm{k}}^{*}}{\partial\xi^{*}}=\left(\frac{\partial\mathcal{G}_{\bm{k}}}{\partial\xi}\right)^{*}.$
(108)

For the second part of (107) we can conjugate the upstream gradient tensor (as
the loss function is real):

$\displaystyle\frac{\partial
L}{\partial\mathcal{G}_{\bm{k}}}=\left(\frac{\partial
L}{\partial\mathcal{G}_{\bm{k}}^{*}}\right)^{*},$ (109)

but as the matrix elements of $\mathcal{G}_{\bm{k}}$ are complex and non-
holomorphic functions of $\xi$, their derivatives with respect to $\xi$ and
$\xi^{*}$ are independent and we need to compute
$\partial\mathcal{G}_{\bm{k}}/\partial\xi^{*}$ separately. On the other hand,
for a real parameter $x$ computing a single gradient tensor
$\frac{\partial\mathcal{G}_{\bm{k}}}{\partial x}$ would suffice:

$\displaystyle\frac{\partial L}{\partial x}$
$\displaystyle=\sum_{\bm{k}}\frac{\partial
L}{\partial\mathcal{G}_{\bm{k}}^{*}}\frac{\partial\mathcal{G}_{\bm{k}}^{*}}{\partial
x}+\frac{\partial
L}{\partial\mathcal{G}_{\bm{k}}}\frac{\partial\mathcal{G}_{\bm{k}}}{\partial
x}$ (110) $\displaystyle=2\Re\left(\sum_{\bm{k}}\frac{\partial
L}{\partial\mathcal{G}_{\bm{k}}}\frac{\partial\mathcal{G}_{\bm{k}}}{\partial
x}\right),$ (111)

as in this case $\frac{\partial\mathcal{G}_{\bm{k}}^{*}}{\partial x}$ and
$\frac{\partial\mathcal{G}_{\bm{k}}}{\partial 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 variables,
e.g., for $\xi=re^{i\phi}$:

$\displaystyle\frac{\partial L}{\partial r}$
$\displaystyle=2\Re\left(\frac{\partial
L}{\partial\xi^{*}}\frac{\partial\xi^{*}}{\partial
r}\right)=2\Re\left(\frac{\partial L}{\partial\xi^{*}}e^{-i\phi}\right),$
(112) $\displaystyle\frac{\partial L}{\partial\phi}$
$\displaystyle=2\Re\left(\frac{\partial
L}{\partial\xi^{*}}\frac{\partial\xi^{*}}{\partial\phi}\right)=-2\Re\left(\frac{\partial
L}{\partial\xi^{*}}i\xi^{*}\right).$ (113)

Furthermore, if the gate depends on a single parameter (real or complex) more
simplifications can be made as shown in Appendix F.

Figure 4: Benchmarks for the computation of the squeezing $S(\zeta)$,
displacement $D(\gamma)$ and beamsplitter $B(\theta,\varphi)$ 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.

###  4.2 Computing gradients with the generating function method

An $\ell$-mode Gaussian gate has $2\ell^{2}+3\ell$ real parameters (counting
double for complex parameters), therefore in order to obtain all the parameter
updates we need to compute $2\ell^{2}+3\ell$ derivatives, some of which will
be combined 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 matrix elements we can simply differentiate the generating
function with respect to the parameters _before_ computing the high-order
derivatives (here with respect to a generic parameter $\xi$):

$\displaystyle\partial_{\xi}\mathcal{G}_{\bm{k}}$
$\displaystyle=\frac{\partial_{\bm{\nu}}^{\bm{k}}}{\sqrt{\bm{k}!}}\partial_{\xi}\Gamma(\bm{\nu})|_{\bm{\nu}=\bm{0}}=\frac{\partial_{\bm{\nu}}^{\bm{k}}}{\sqrt{\bm{k}!}}\partial_{\xi}Ce^{Q(\bm{\nu})}|_{\bm{\nu}=\bm{0}}$
(114)
$\displaystyle=\frac{\partial_{\bm{\nu}}^{\bm{k}}}{\sqrt{\bm{k}!}}\Gamma(\bm{\nu})\left(\frac{\partial_{\xi}C}{C}+\partial_{\xi}Q(\bm{\nu})\right)|_{\bm{\nu}=\bm{0}}$
(115)
$\displaystyle=\left.\frac{\partial_{\bm{\nu}}^{\bm{k}}}{\sqrt{\bm{k}!}}\Gamma(\bm{\nu})\left(\frac{\partial_{\xi}C}{C}+\frac{\partial\bm{\mu}^{T}}{\partial\xi}{\bm{\nu}}-\frac{1}{2}{\bm{\nu}}^{T}\frac{\partial\bm{\Sigma}}{\partial\xi}{\bm{\nu}}\right)\right|_{\bm{\nu}=\bm{0}}.$
(116)

Expression (116) contains terms proportional to
$\bm{\nu}^{\bm{n}}\Gamma(\bm{\nu})$, whose derivative can be expressed using
the elements of $\mathcal{G}_{\bm{k}}$:

$\displaystyle\frac{\partial^{\bm{k}}_{\bm{\nu}}}{\sqrt{\bm{k}!}}{\bm{\nu}}^{\bm{n}}\Gamma(\bm{\nu})|_{\bm{\nu}=\bm{0}}=\sqrt{(\bm{k})_{\bm{n}}}\mathcal{G}_{\bm{k}-\bm{n}},$
(117)

where $(k)_{n}=k(k-1)(k-2)\dots(k-n+1)$ is the Pochhammer symbol and
$(\bm{k})_{\bm{n}}=\prod_{i}(k_{i})_{n_{i}}$ (with the convention that
$(k)_{0}=1$). So once we compute the tensor representation of a Gaussian
operator, we can use it to compute its gradient with respect to a
parametrization. Continuing from Eq. (116) we can write:

$\displaystyle\frac{\partial\mathcal{G}_{\bm{k}}}{\partial\xi}=$
$\displaystyle\frac{\partial_{\xi}C}{C}\mathcal{G}_{\bm{k}}+\sum_{i}\frac{\partial\mu_{i}}{\partial\xi}\sqrt{k_{i}}\mathcal{G}_{\bm{k}-1_{i}}$
$\displaystyle-\sum_{i>j}\frac{\partial\Sigma_{ij}}{\partial\xi}\sqrt{k_{i}k_{j}}\mathcal{G}_{\bm{k}-1_{i}-1_{j}}$
$\displaystyle-\frac{1}{2}\sum_{i}\frac{\partial\Sigma_{ii}}{\partial\xi}\sqrt{k_{i}(k_{i}-1)}\mathcal{G}_{\bm{k}-2_{i}}.$
(118)

This, together with Eq. (36), 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 $\gamma$ and $\gamma^{*}$
(recall that we must treat them as independent variables). We begin by
differentiating $C$, $\bm{\mu}$ and $\bm{\Sigma}$:

$\displaystyle\frac{\partial_{\gamma}C}{C}$
$\displaystyle=-\frac{\gamma^{*}}{2},$ (119)
$\displaystyle\frac{\partial_{\gamma^{*}}C}{C}$
$\displaystyle=-\frac{\gamma}{2}-\gamma^{*}e^{i(\delta+2\phi)}\tanh r,$ (120)
$\displaystyle\frac{\partial\bm{\mu}}{\partial\gamma}$ $\displaystyle=[1,0],$
(121) $\displaystyle\frac{\partial\bm{\mu}}{\partial\gamma^{*}}$
$\displaystyle=[e^{i(\delta+2\phi)}\tanh r,-e^{i\phi}\operatorname{sech}r],$
(122) $\displaystyle\frac{\partial\bm{\Sigma}}{\partial\gamma}$
$\displaystyle=\bm{0},$ (123)
$\displaystyle\frac{\partial\bm{\Sigma}}{\partial\gamma^{*}}$
$\displaystyle=\bm{0}.$ (124)

Therefore, the gradients of the matrix elements are:

$\displaystyle\frac{\partial\mathcal{G}_{m,n}}{\partial\gamma}=$
$\displaystyle-\frac{\gamma^{*}}{2}\,\mathcal{G}_{m,n}+\sqrt{m}\,\mathcal{G}_{m-1,n},$
(125) $\displaystyle\frac{\partial\mathcal{G}_{m,n}}{\partial\gamma^{*}}=$
$\displaystyle-\left(\frac{\gamma}{2}+\frac{1}{2}\gamma^{*}e^{i(\delta+2\phi)}\tanh
r\right)\mathcal{G}_{m,n}$ (126)
$\displaystyle+e^{i(\delta+2\phi)}\sqrt{m}\mathcal{G}_{m-1,n}\tanh r$
$\displaystyle-e^{i\phi}\sqrt{n}\mathcal{G}_{m,n-1}\operatorname{sech}r.$

We then combine both of these expressions with the upstream gradient $\partial
L/\partial\mathcal{G}_{mn}^{*}$ and its conjugate as prescribed in Eq. (107)
to obtain the single gradient for the update, $\partial L/\partial\gamma^{*}$.
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 benchmark our methods and implementations for the
displacement, squeezing and beamsplitter gate against the implementations of
Strawberry Fields (SF) version 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 summarized 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\sim 100$. Our explicit recurrence
relations avoid completely this issue for any number of modes.

Figure 5: Benchmarks for the computation of the beamsplitter
$BS(\theta,\varphi)$ using Bosonic and the methods from this work as
implemented in The Walrus. Note that for a given cutoff $N$, Bosonic
calculates $\mathcal{B}^{(N)}_{i,j}=\langle
N-i,i|BS(\theta,\varphi)|N-j,j\rangle$ and that in their implementation the
constraint $N\leq 11$ is hard-coded. For a fair comparison we benchmark their
code against a wrapper function that calculates the tensor $\langle
m,n|BS(\theta,\varphi)|k,l\rangle$ and then constructs the relevant matrix
$\mathcal{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 single core of an Intel-i5
processor @ 2.50 GHz.

Second, we compare our methods against the implementations 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 these modest sizes our
implementation achieves almost two orders of magnitude in improvement.

_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.

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:

$\displaystyle K(\kappa)=\exp\left[i\kappa(a^{\dagger}a)^{2}\right],$ (127)

and is therefore diagonal in the photon number basis. Such single-mode circuit
corresponds to a unitary transformation which we can write as a stack of $M$
layers:

$\displaystyle
U(\bm{\gamma},\bm{\phi},\bm{\zeta},\bm{\kappa})=\prod_{m=1}^{M}\mathcal{G}^{(1)}(\gamma_{m},\phi_{m},\zeta_{m})K(\kappa_{m}).$
(128)

We then input the single-mode vacuum state $|0\rangle$ and optimize the
parameters such that the unitary $U(\bm{\gamma},\bm{\phi},\bm{\zeta})$
produces an output as close as possible to the desired target state. We choose
three target states: a single photon state $|1\rangle$, an “ON” state, i.e. a
superposition of the vacuum and a Fock state, in our case
$(|0\rangle+|9\rangle)/\sqrt{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|0\rangle$ and the target:

$\displaystyle
L(|\psi_{\mathrm{out}}\rangle)=-|\langle\psi_{\mathrm{out}}|\psi_{\mathrm{target}}\rangle|.$
(129)

For the optimization we run the adaptive gradient descent algorithm `Adam`. We
note that the range over which the parameters are optimized is unconstrained,
except for the fact that the parameters $\phi_{m}$ and $\kappa_{m}$ are real.

The results are presented in Table 1 and can be compared 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 Gaussian transformation using generating functions. From the
generating functions we obtain recurrence relations which allow us to compute
the matrix elements of any given Gaussian transformation when written in terms
of its Bloch-Messiah decomposition. Furthermore, we have developed general
techniques to also obtain 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, optical data
processing, device discovery and device design.

_Note added –_ While preparing this manuscript we became 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áva for valuable discussions. F.M. thanks Electra
Eleftheriadou for support throughout this work, Yuan Yao for helping with the
correctness of the final version of the manuscript, Frédéric Grosshans, Steve
Barnett, Gerardo 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.

##  Appendix A Multimode generating function

In Barnett and Radmore [25] (cf. Sec. 3.7 and Appendix 5) the following
identities are proven:

$\displaystyle\exp\left(\theta a^{\dagger}a\right)$ $\displaystyle=\
:\exp\left(\left[\exp(\theta)-1\right]a^{\dagger}a\right):=\sum_{n=0}^{\infty}\frac{\left[\exp(\theta)-1\right]^{n}a^{\dagger
n}a^{n}}{n!},$ (130)

and also

$\displaystyle{S}(\zeta)=$
$\displaystyle\exp\left(\frac{1}{2}\zeta^{*}a^{2}-\text{H.c.}\right)=\exp\left(\frac{1}{2}re^{-i\delta}a^{2}-\text{H.c.}\right)$
$\displaystyle=$ $\displaystyle\exp\left(-\tfrac{\tanh
r}{2}e^{i\delta}a^{\dagger
2}\right)\exp\left(-\left[a^{\dagger}a+\tfrac{1}{2}\right]\log\cosh
r\right)\exp\left(\tfrac{\tanh r}{2}e^{i\delta}a^{2}\right).$ (131)

where we wrote $\zeta=re^{i\delta}$. Putting these two identities together we
can normal order a squeezing operator

$\displaystyle\exp\left(\frac{1}{2}re^{-i\delta}a^{2}-\text{H.c.}\right)=$
$\displaystyle\exp\left(-\tfrac{\tanh r}{2}e^{i\delta}a^{\dagger
2}\right)\times\frac{:\exp\left(\left[\operatorname{sech}r-1\right]\
a^{\dagger}a\right):}{\sqrt{\cosh r}}\times\exp\left(\tfrac{\tanh
r}{2}e^{-i\delta}a^{2}\right).$ (132)

Using this expression one can find the matrix element of the squeezing
operation between two coherent states

$\displaystyle\langle\alpha^{*}|\exp\left(\frac{1}{2}re^{-i\delta}a^{2}-\text{H.c.}\right)|\beta\rangle=\frac{1}{\sqrt{\cosh
r}}\langle\alpha^{*}|\beta\rangle\exp\left(-\tfrac{\tanh
r}{2}e^{i\delta}\alpha^{2}+\left[\operatorname{sech}r-1\right]\
\alpha\beta+\tfrac{\tanh r}{2}e^{-i\delta}\beta^{2}\right),$ (133)

and noting that

$\displaystyle\langle\alpha^{*}|\beta\rangle=\exp\left(-\tfrac{1}{2}\left[|\alpha|^{2}+|\beta|^{2}-2\alpha\beta\right]\right),$
(134)

one can write

$\displaystyle\langle\alpha^{*}|$
$\displaystyle\exp\left(\frac{1}{2}re^{i\delta}a^{2}-\text{H.c.}\right)|\beta\rangle=\frac{\exp\left(-\tfrac{1}{2}\left[|\alpha|^{2}+|\beta|^{2}\right]\right)}{\sqrt{\cosh
r}}\exp\left(-\frac{1}{2}\left[\alpha\
\beta\right]\left[\begin{array}[]{cc}e^{i\delta}\tanh
r&-\operatorname{sech}r\\\ -\operatorname{sech}r&-e^{-i\delta}\tanh r\\\
\end{array}\right]\left[\begin{array}[]{c}\alpha\\\ \beta\\\
\end{array}\right]\right).$ (139)

This proof generalizes in a straightforward way to

$\displaystyle\langle\bm{\alpha}^{*}|\mathcal{S}(\bm{\zeta})|\bm{\beta}\rangle=\frac{\exp\left(-\tfrac{1}{2}\left[||\bm{\alpha}||^{2}+||\bm{\beta}||^{2}\right]\right)}{\sqrt{\prod_{i=1}^{\ell}\cosh
r_{i}}}\exp\left(-\frac{1}{2}\left[\bm{\alpha}\
\bm{\beta}\right]\left[\begin{array}[]{c|c}\text{diag}(e^{i\bm{\delta}}\tanh\bm{r})&-\text{diag}(\operatorname{sech}\bm{r})\\\
\hline\cr-\text{diag}(\operatorname{sech}\bm{r})&-\text{diag}(e^{-i\bm{\delta}}\tanh\bm{r})\\\
\end{array}\right]\left[\begin{array}[]{c}\bm{\alpha}\\\ \bm{\beta}\\\
\end{array}\right]\right),$ (144)

where $\zeta_{j}=r_{j}e^{i\delta_{j}}$.

From the definitions in Eq. (26) it is direct to show that

$\displaystyle\mathcal{U}(\bm{V})|\bm{\beta}\rangle=|\bm{V}\bm{\beta}\rangle\text{
and
}\langle\bm{\alpha}^{*}|\mathcal{U}(\bm{W})=\langle\bm{W}^{\dagger}\bm{\alpha}^{*}|.$
(145)

Finally, let us note that

$\displaystyle\langle\bm{\alpha}^{*}|\mathcal{D}(\bm{\gamma})=\langle\bm{\alpha}^{*}-\bm{\gamma}|\exp\left(\tfrac{1}{2}\left[\bm{\alpha}^{T}\bm{\gamma}-\bm{\gamma}^{\dagger}\bm{\alpha}^{*}\right]\right).$
(146)

To show this last identity recall that
$\langle\bm{\alpha}^{*}|=\langle\bm{0}|\mathcal{D}(-\bm{\alpha}^{*})$ 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. 2.2 as follows:

$\displaystyle\langle\bm{\alpha}^{*}$
$\displaystyle|\mathcal{D}(\bm{\gamma})\mathcal{U}(\bm{W})\mathcal{S}(\bm{\zeta})\mathcal{U}(\bm{V})|\bm{\beta}\rangle$
(147a)
$\displaystyle=\exp\left(\tfrac{1}{2}\left[\bm{\alpha}^{T}\bm{\gamma}-\bm{\gamma}^{\dagger}\bm{\alpha}^{*}\right]\right)\langle\bm{\alpha}^{*}-\bm{\gamma}|\mathcal{U}(\bm{W})\mathcal{S}(\bm{\zeta})|\bm{V}\bm{\beta}\rangle$
(147b)
$\displaystyle=\exp\left(\tfrac{1}{2}\left[\bm{\alpha}^{T}\bm{\gamma}-\bm{\gamma}^{\dagger}\bm{\alpha}^{*}\right]\right)\langle\bm{W}^{\dagger}(\bm{\alpha}^{*}-\bm{\gamma})|\mathcal{S}(\bm{\zeta})|\bm{V}\bm{\beta}\rangle$
(147c)
$\displaystyle=\frac{\exp\left(\tfrac{1}{2}\left[\bm{\alpha}^{T}\bm{\gamma}-\bm{\gamma}^{\dagger}\bm{\alpha}^{*}\right]\right)}{\sqrt{\prod_{i=1}^{\ell}\cosh
r_{i}}}\exp\left(-\tfrac{1}{2}\left[||\bm{\alpha}^{*}-\bm{\gamma}||^{2}+|\bm{\beta}|^{2}\right]\right)\times$
(147d)
$\displaystyle\quad\quad\exp\left(-\tfrac{1}{2}\left[(\bm{\alpha}^{T}-\bm{\gamma}^{\dagger})\bm{W},\
\
\bm{\beta}^{T}\bm{V}^{T}\right]\left[\begin{array}[]{c|c}\text{diag}(e^{i\bm{\delta}}\tanh\bm{r})&-\text{diag}(\operatorname{sech}\bm{r})\\\
\hline\cr-\text{diag}(\operatorname{sech}\bm{r})&-\text{diag}(e^{-i\bm{\delta}}\tanh\bm{r})\\\
\end{array}\right]\left[\begin{array}[]{c}\bm{W}^{T}(\bm{\alpha}-\bm{\gamma}^{*})\\\
\bm{V}\bm{\beta}\\\ \end{array}\right]\right).$ (147i)

At this point we introduce the symmetric, complex-unitary matrix

$\displaystyle\bm{\Sigma}$
$\displaystyle=\left[\begin{array}[]{c|c}\bm{W}\text{diag}(e^{i\bm{\delta}}\tanh\bm{r})\bm{W}^{T}&-\bm{W}\text{diag}(\operatorname{sech}\bm{r})\bm{V}\\\
\hline\cr-\bm{V}^{T}\text{diag}(\operatorname{sech}\bm{r})\bm{W}^{T}&-\bm{V}^{T}\text{diag}(e^{-i\bm{\delta}}\tanh\bm{r})\bm{V}\end{array}\right]$
(150) $\displaystyle=\left[\begin{array}[]{c|c}\bm{W}&\bm{0}\\\
\hline\cr\bm{0}&\bm{V}^{T}\end{array}\right]\left[\begin{array}[]{c|c}\text{diag}(e^{i\bm{\delta}}\tanh\bm{r})&-\text{diag}(\operatorname{sech}\bm{r})\\\
\hline\cr-\text{diag}(\operatorname{sech}\bm{r})&-\text{diag}(e^{-i\bm{\delta}}\tanh\bm{r})\end{array}\right]\left[\begin{array}[]{c|c}\bm{W}&\bm{0}\\\
\hline\cr\bm{0}&\bm{V}^{T}\end{array}\right]^{T},$ (157)

and the complex vector

$\displaystyle\bm{\nu}$ $\displaystyle=\left[\begin{array}[]{c}\bm{\alpha}\\\
\bm{\beta}\end{array}\right]\in\mathbb{C}^{2\ell},$ (160)

to finally write

$\displaystyle\langle\bm{\alpha}^{*}|\mathcal{D}(\bm{\gamma})\mathcal{U}(\bm{W})\mathcal{S}(\bm{\zeta})\mathcal{U}(\bm{V})|\bm{\beta}\rangle=\exp\left(-\tfrac{1}{2}\left[||\bm{\alpha}||^{2}+||\bm{\beta}||^{2}\right]\right)C\exp\left(\bm{\mu}^{T}\bm{\nu}-\frac{1}{2}\bm{\nu}^{T}\bm{\Sigma}\bm{\nu}\right),$
(161)

where

$\displaystyle C$
$\displaystyle=\frac{\exp\left(-\tfrac{1}{2}\left[||\bm{\gamma}||^{2}+\bm{\gamma}^{\dagger}\bm{W}\text{diag}(e^{i\bm{\delta}}\tanh\bm{r})\bm{W}^{T}\bm{\gamma}^{*}\right]\right)}{\sqrt{\prod_{i=1}^{\ell}\cosh
r_{i}}},$ (162) $\displaystyle\bm{\mu}^{T}$
$\displaystyle=\left[\bm{\gamma}^{\dagger}\bm{W}\text{diag}(e^{i\bm{\delta}}\tanh\bm{r})\bm{W}^{T}+\bm{\gamma}^{T},-\bm{\gamma}^{\dagger}\bm{W}\text{diag}(\operatorname{sech}\bm{r})\bm{V}\right],$
(163)

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].

##  Appendix B Recurrence relation for the high-order derivatives of an
exponential function

We rewrite Eq. (15) of Ref. [45] in our notation as follows:

$\displaystyle\Gamma_{\bm{k}+1_{n}}=\sum_{\bm{j}=\bm{0}}^{\bm{k}}\binom{\bm{k}}{\bm{j}}Q^{(\bm{k}-\bm{j}+1_{n})}\Gamma_{\bm{j}},$
(164)

where we write
$\partial_{\bm{\nu}}^{\bm{k}}f(\bm{\nu})|_{\bm{\nu}=\bm{0}}=f^{(\bm{k})}=f_{\bm{k}}$
for the $\bm{k}^{\text{th}}$ derivative of $f$ at zero and define
$\sum_{\bm{j}=\bm{0}}^{\bm{k}}\binom{\bm{k}}{\bm{j}}=\sum_{j_{1}}^{k_{1}}\dots\sum_{j_{n}}^{k_{n}}\binom{k_{1}}{j_{1}}\dots\binom{k_{n}}{j_{n}}$.
In our case $Q$ is a quadratic in $\bm{\nu}$ and we write

$\displaystyle
Q(\bm{\nu})=\bm{\mu}^{T}\bm{\nu}-\frac{1}{2}\bm{\nu}^{T}\bm{\Sigma}\bm{\nu}.$
(165)

We can write easily all its nonzero derivatives:

$\displaystyle Q^{(1_{i})}$ $\displaystyle=\mu_{i},$ (166) $\displaystyle
Q^{(1_{i}+1_{j})}$ $\displaystyle=-\Sigma_{i,j}.$ (167)

We know that most of the terms in Eq. (164) will not survive the summation,
indeed only the ones that satisfy $\bm{k}-\bm{j}=1_{i}$ or
$\bm{k}-\bm{j}=1_{i}+1_{l}$ will contribute to the sum. We can use this to our
advantage and write

$\displaystyle\Gamma_{\bm{k}+1_{n}}$
$\displaystyle=Q^{(1_{n})}\Gamma_{\bm{k}}+\sum_{l}k_{l}Q^{(1_{n}+1_{l})}\Gamma_{\bm{k}-1_{l}}=\mu_{n}\Gamma_{\bm{k}}-\sum_{l}k_{l}\Sigma_{l,n}\Gamma_{\bm{k}-1_{l}}.$
(168)

Now let us derive this in a different way. Start with the generating function
of the multidimensional Hermite polynomials [57, 58, 59]:

$\displaystyle\exp\left(\bm{\mu}^{T}\bm{\nu}-\tfrac{1}{2}\bm{\nu}^{T}\bm{\Sigma}\bm{\nu}\right)=\sum_{\bm{k}=\bm{0}}^{\bm{\infty}}\frac{\bm{\nu}^{\bm{k}}}{\bm{k}!}G_{\bm{k}}^{(\bm{\Sigma})}(\bm{\mu}).$
(169)

They satisfy the recurrence relation

$\displaystyle
G_{\bm{k}+1_{i}}^{(\bm{\bm{\Sigma}})}(\bm{\mu})=\mu_{i}G_{\bm{k}}^{(\bm{\bm{\Sigma}})}(\bm{\mu})-\sum_{j=1}\Sigma_{i,j}k_{j}G_{\bm{k}-1_{j}}^{(\bm{\bm{\Sigma}})}(\bm{\mu}),$

from which we identify that for a quadratic function $Q$,

$\displaystyle\Gamma_{\bm{k}}=G^{(\bm{\Sigma})}_{\bm{k}}(\bm{\mu})=\sqrt{\bm{k}!}\mathcal{G}_{\bm{k}}.$
(170)

By replacing this last expression in the recurrence relation Eq. (168), and
simplifying the factorials one arrives at the recurrence relation Eq. (35) for
the matrix elements $\mathcal{G}_{\bm{k}}$.

##  Appendix 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^{\text{th}}$-order-phase gate defined as

$\displaystyle V^{(k)}(\eta)=\exp\left(i\frac{\eta}{k\hbar}q^{k}\right),$
(171)

where the canonical position is defined as
$q=\sqrt{\frac{\hbar}{2}}(a+a^{\dagger})$ and its canonical momentum is given
by $p=i\sqrt{\tfrac{\hbar}{2}}\left(a^{\dagger}-a\right)$. In the Heisenberg
picture the $k^{\text{th}}$-order-phase gate transforms the operators
according to

$\displaystyle V^{(k)\dagger}(\eta)qV^{(k)}(\eta)$ $\displaystyle=q,$ (172)
$\displaystyle V^{(k)\dagger}(\eta)pV^{(k)}(\eta)$ $\displaystyle=p+\eta
q^{k-1},$ (173) $\displaystyle V^{(k)\dagger}(\eta)aV^{(k)}(\eta)$
$\displaystyle=a+i\frac{\eta}{\hbar}\left(\frac{\hbar}{2}\right)^{k/2}\left(a+a^{\dagger}\right)^{k-1}.$
(174)

One can write the matrix elements of the $k^{\text{th}}$-order-phase gate in
terms of integrals as follows

$\displaystyle\langle m|$ $\displaystyle
V^{(k)}(\eta)|n\rangle=V_{m,n}^{(k)}=\frac{1}{\sqrt{\pi 2^{n+m}n!m!}}\int
dx\exp\left(-x^{2}+i\frac{\eta\hbar^{k/2}}{\hbar
k}x^{k}\right)~{}H_{n}(x)H_{m}(x),$ (175)

where $H_{n}(x)$ is a Hermite polynomial of degree $n$ [25]. The last equation
makes explicit that the matrix elements $V_{m,n}^{(k)}$ are symmetric
$V_{m,n}^{(k)}=V_{n,m}^{(k)}$ and moreover, given the parity of the Hermite
polynomials $H_{n}(-x)=(-1)^{n}H_{n}(x)$, that the $k^{\text{th}}$-order-phase
gate also has parity, i.e.,

$\displaystyle V_{m,n}^{(k)}(\eta)$
$\displaystyle=(-1)^{m+n}V_{m,n}^{(k)}(-\eta)\text{ if $k$ is odd,}$ (176)
$\displaystyle V_{m,n}^{(k)}(\eta)$ $\displaystyle=0\text{ if $k$ is even and
$n+m$ is odd.}$ (177)

Finally, note that the matrix elements of the conjugate $k^{\text{th}}$-order-
phase gate

$\displaystyle\tilde{V}^{(k)}(\eta)=\exp\left(i\frac{\eta}{\hbar
k}p^{k}\right),$ (178)

are related to the ones for the $k^{\text{th}}$-order-phase gate via

$\displaystyle\langle
m|\tilde{V}^{(k)}(\eta)|n\rangle=\tilde{V}^{(k)}_{m,n}=(i)^{m+n}V_{m,n}^{(k)},$
(179)

which can be easily confirmed by noticing that the rotation gate by angle
$\theta=\pi/2$ transforms $q$ into $p$ in the Heisenberg picture and that when
acted on a Fock state it gives $R(\pi/2)|n\rangle=(i)^{n}|n\rangle$.

Having established these basic facts about the $k^{\text{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

$\displaystyle[a,V^{(k)}(\eta)]=aV^{(k)}(\eta)-V^{(k)}(\eta)a=V^{(k)}(\eta)V^{(k)\dagger}(\eta)aV^{(k)}(\eta)-V^{(k)}(\eta)a=V^{(k)}\left\\{i\frac{\eta}{\hbar}\left(\frac{\hbar}{2}\right)^{k/2}(a+a^{\dagger})^{k-1}\right\\}.$
(180)

If we multiply the last equation on the left by $\langle m|$ and on the right
$|n\rangle$ we obtain

$\displaystyle\sqrt{m+1}V^{(k)}_{m+1,n}-\sqrt{n}V^{(k)}_{m,n-1}=\langle
m|V^{(k)}(\eta)\left(i\frac{\eta}{\hbar}\left(\frac{\hbar}{2}\right)^{k/2}(a+a^{\dagger})^{k-1}\right)|n\rangle.$
(181)

Note that for any $k$ one can expand the right hand side of the last equation
and write
$(a+a^{\dagger})^{k-1}|n\rangle=\sum_{k=n-(k-1)}^{n+k-1}c_{k}|k+n\rangle$
allowing us to obtain the desired recurrence relation. In the next subsection
we specialize 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. (171)
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, 60, 24]
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, 65, 66, 67].

To obtain the recurrence relation for the cubic phase gate, we set $k=3$ in
Eq. (181) and rearrange to obtain

$\displaystyle
V^{(3)}_{m,n+2}=\frac{1}{\sqrt{(n+1)(n+2)}}\biggl{(}\frac{i2\sqrt{2}}{\sqrt{\hbar}\eta}\left\\{\sqrt{n}V^{(3)}_{m,n-1}-\sqrt{m+1}V^{(3)}_{m+1,n}\right\\}-(2n+1)V^{(3)}_{m,n}-\sqrt{n(n-1)}V^{(3)}_{m,n-2}\biggr{)}.$
(182)

We can obtain the first two rows of $V_{m,n}^{(3)}$ by setting $m=0$ or $1$ in
the last equation and using the following initial seeds (cf. Appendix D)

$\displaystyle V_{0,0}^{(3)}$
$\displaystyle=2\sqrt{\pi}e^{\frac{2y^{6}}{3}}|y|\text{Ai}\left(y^{4}\right),$
(183) $\displaystyle V_{1,1}^{(3)}$
$\displaystyle=-2i\sqrt{2}y^{3}V_{0,1}^{(3)}=-2i\sqrt{2}y^{3}V_{1,0}^{(3)}=-8\sqrt{\pi}e^{\frac{2y^{6}}{3}}|y|^{5}\left(\text{Ai}^{\prime}\left(y^{4}\right)+y^{2}\text{Ai}\left(y^{4}\right)\right),$
(184)

where $y^{3}=1/(\sqrt{\hbar}\eta)$ and $\text{Ai}(x)$ and
$\text{Ai}^{\prime}(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 $\eta<1$; notice that any small
numerical error in the difference
$\sqrt{n}V^{(3)}_{m,n-1}-\sqrt{m+1}V^{(3)}_{m+1,n}$ inside the curly braces in
Eq. 182 will be _amplified_ by a factor of $1/\eta$ that will cause an
uncotrollable loss of precision; because of this, the recurrence relation can
only be safely used in the regime where $\eta>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[72, 50]. We can obtain a recursive relation
for the matrix elements of this gate by setting $k=4$ in Eq. (181) to obtain

$\displaystyle V^{(4)}_{m+3,n}=\frac{1}{\sqrt{(n+1)(n+2)(n+3)}}\biggl{(}$
$\displaystyle\frac{4i}{\eta\hbar}\left\\{\sqrt{n}V^{(4)}_{m,n-1}-\sqrt{m+1}V^{(4)}_{m+1,n}\right\\}-\sqrt{n(n-1)(n-2)}V^{(4)}_{m,n-3}$
$\displaystyle-3n^{\tfrac{3}{2}}V^{(4)}_{m,n-1}-3(n+1)^{\tfrac{3}{2}}V^{(4)}_{m,n+1}\biggr{)}.$
(185)

We can iterate this relation using the following initial values

$\displaystyle V^{(4)}_{0,0}=$
$\displaystyle\left.f(w,\lambda)\right|_{\lambda=1},$ (186) $\displaystyle
V^{(4)}_{1,1}=$
$\displaystyle\left.-2\frac{\partial}{\partial\lambda}f(w,\lambda)\right|_{\lambda=1},$
(187) $\displaystyle V_{0,2}^{(4)}=$ $\displaystyle
V_{2,0}^{(4)}=\frac{1}{\sqrt{2}}\left(V_{1,1}^{(4)}-V_{0,0}^{(4)}\right),$
(188) $\displaystyle V_{2,2}^{(4)}=$
$\displaystyle\frac{1}{2}V_{0,0}^{(4)}-V_{1,1}^{(4)}+2\left.\frac{\partial^{2}}{\partial\lambda^{2}}f(w,\lambda)\right|_{\lambda=1},$
(189)

where $f(w,\lambda)$ is defined in Eq. (198) and recalling that
$V_{0,1}^{(4)}=V_{1,0}^{(4)}=V_{1,2}^{(4)}=V_{2,1}^{(4)}=0$ by parity. With
these nine initial values we can fill the first three rows of $V_{m,n}^{(4)}$
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. (C.2) is numerically
unstable for $\eta<1$.

##  Appendix 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. (175) with $k=3$ and $y^{3}=1/(\sqrt{\hbar}\eta)$ and doing the
change of variables $x\to yk-iy^{3}$, and assuming without loss of generality
that $y>0$, we obtain

$\displaystyle V^{(3)}_{m,n}=\frac{ye^{\frac{2}{3}y^{6}}}{\sqrt{\pi
2^{n+m}n!m!}}\int_{-\infty}^{\infty}dk\exp\left(iy^{4}k+ik^{3}/3\right)~{}H_{n}(yk-
iy^{3})H_{m}(yk-iy^{3}).$ (190)

To make progress we recall the identity
$H_{n}(x+y)=\sum_{k=0}^{n}\binom{n}{k}H_{k}(x)(2y)^{(n-k)}$, to obtain

$\displaystyle V^{(3)}_{m,n}$
$\displaystyle=\frac{ye^{\frac{2}{3}y^{6}}}{\sqrt{\pi
2^{n+m}n!m!}}\int_{-\infty}^{\infty}dk\exp\left(iy^{4}k+ik^{3}/3\right)~{}\sum_{k=0}^{n}\binom{n}{k}H_{k}(-iy^{3})(2yk)^{(n-k)}\sum_{l=0}^{m}\binom{m}{l}H_{l}(-iy^{3})(2yk)^{(m-l)}$
(191) $\displaystyle=\frac{ye^{\frac{2}{3}y^{6}}}{\sqrt{\pi
2^{n+m}n!m!}}\sum_{k=0}^{n}\sum_{l=0}^{m}(2y)^{m+n-k-l}\binom{n}{k}H_{k}(-iy^{3})\binom{m}{l}H_{l}(-iy^{3})\int_{-\infty}^{\infty}dk\exp\left(iy^{4}k+ik^{3}/3\right)k^{n+m-k-l}$
(192)
$\displaystyle=\frac{2\sqrt{\pi}ye^{\frac{2}{3}y^{6}}}{\sqrt{2^{n+m}n!m!}}\sum_{k=0}^{n}\sum_{l=0}^{m}(-2iy)^{m+n-k-l}\binom{n}{k}H_{k}(-iy^{3})\binom{m}{l}H_{l}(-iy^{3})\Bigg{.}\left[\frac{d^{n+m-l-k}}{dw^{n+m-l-k}}\operatorname{Ai}(w)\right]\Bigg{|}_{w=y^{4}}.$
(193)

In the last equation, we used the well known Fourier expansion of the Airy
function to write [69]

$\displaystyle\operatorname{Ai}(x)$
$\displaystyle=\frac{1}{2\pi}\int_{-\infty}^{\infty}dk\exp\left(i\frac{k^{3}}{3}+ikx\right),$
(194) $\displaystyle\frac{d^{\ell}}{dx^{\ell}}\operatorname{Ai}(x)$
$\displaystyle=\frac{1}{2\pi}\int_{-\infty}^{\infty}dk(ik)^{\ell}\exp\left(i\frac{k^{3}}{3}+ikx\right).$
(195)

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

$\displaystyle\frac{d^{\ell}}{dx^{\ell}}\operatorname{Ai}(x)=p_{\ell}(x)\operatorname{Ai}(x)+q_{\ell}(x)\mathrm{Ai^{\prime}}(x),$
(196)

where the polynomials $p_{\ell}$ and $q_{\ell}$ satisfy the recurrence
relations

$\displaystyle p_{\ell+1}(x)=xq_{\ell}(x)+\frac{d}{dx}p_{\ell}(x),\quad
q_{\ell+1}(x)=p_{\ell}(x)+\frac{d}{dx}q_{\ell}(x).$ (197)

By iterating the recurrence relations in Eq. (197) one can in principle obtain
the polynomials $p_{\ell}(x)$ and $q_{\ell}(x)$ and then any derivative of the
Airy function and thus any matrix element $V^{(3)}_{n,m}$.

##  Appendix E Analytic matrix elements of the quartic-phase gate

We are now interested in calculating

$\displaystyle V^{(4)}_{m,n}=\langle
m|V^{(4)}(\eta)|n\rangle=\frac{1}{\sqrt{\pi
2^{n+m}n!m!}}\int_{-\infty}^{\infty}dx\exp\left(-x^{2}+\tfrac{i}{4}wx^{4}\right)~{}H_{n}(x)H_{m}(x),$

where $H_{n}(x)$ is the $n^{\text{th}}$ order Hermite polynomial and we
introduced $w=\hbar\eta$. First consider the related integral (cf. Eq. 3.968
of Ref. [73])

$\displaystyle f(w,\lambda)$
$\displaystyle=\frac{1}{\sqrt{\pi}}\int_{-\infty}^{\infty}dx\exp\left(-\lambda
x^{2}+\tfrac{i}{4}wx^{4}\right)$ (198)
$\displaystyle=-\frac{i}{4}\sqrt{\pi}e^{\frac{i\lambda^{2}}{8w}+\frac{i\pi}{8}}\sqrt{\frac{\lambda}{w}}\left(J_{\frac{1}{4}}\left(\frac{\lambda^{2}}{8w}\right)-iY_{\frac{1}{4}}\left(\frac{\lambda^{2}}{8w}\right)\right)=-\frac{i}{4}\sqrt{\pi}e^{\frac{i\lambda^{2}}{8w}+\frac{i\pi}{8}}\sqrt{\frac{\lambda}{w}}H^{(2)}_{\frac{1}{4}}\left(\frac{\lambda^{2}}{8w}\right),$
(199)

where $J_{\alpha}(x)$ is a Bessel function of the first kind, $Y_{\alpha}(x)$
is a Bessel function of the second kind (sometimes also called a Neumann
function) and $H_{\alpha}^{(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_{0,0}^{(4)}(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

$\displaystyle\frac{(-1)^{k}}{\sqrt{\pi}}\int_{-\infty}^{\infty}dx\
x^{2k}\exp\left(-\lambda
x^{2}+\tfrac{i}{4}wx^{4}\right)=\frac{\partial^{k}}{\partial\lambda^{k}}f(w,\lambda).$
(200)

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

$\displaystyle\frac{\partial}{\partial
x}H_{i}^{(2)}(x)=\frac{1}{2}\left(H_{i-1}^{(2)}(x)+H_{i+1}^{(2)}(x)\right).$
(201)

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$.

##  Appendix 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\epsilon}$,

$\displaystyle
G(u)=G(xe^{i\epsilon})=\exp\left(xe^{i\epsilon}h-\text{H.c.}\right).$ (202)

The case of a gate parametrized by a single _real_ number can be obtained by
letting $\epsilon=\pi/2$ and writing $G(x)=\exp(ix(h+h^{\dagger}))=\exp(Hx)$
with $H=i(h+h^{\dagger})=-H^{\dagger}$. For example, for a two-mode squeezing
gate, $h=a_{1}^{\dagger}a^{\dagger}_{2}$ and for the cubic phase gate,
$H=iq^{3}/(k\hbar)$. Now it is direct to verify that

$\displaystyle\partial_{x}\langle\bm{m}|G|\bm{n}\rangle=\langle\bm{m}|G\left\\{e^{i\epsilon}h-\text{H.c.}\right\\}|\bm{n}\rangle.$
(203)

To obtain the gradient with respect to the amplitude $x$, one only needs to
express $(e^{i\epsilon}h-\text{H.c.})|\bm{n}\rangle$ 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

$\displaystyle\partial_{r}\langle\bm{m}|S^{(2)}|\bm{n}\rangle=\langle
m_{1},m_{2}|S^{(2)}\
\left\\{a_{1}^{\dagger}a_{2}^{\dagger}e^{i\epsilon}-\text{H.c.}\right\\}|n_{1},n_{2}\rangle$
(204)
$\displaystyle\quad=\sqrt{(n_{1}+1)(n_{2}+1)}e^{i\epsilon}S^{(2)}_{m_{1},m_{2},n_{1}+1,n_{2}+1}-\sqrt{n_{1}n_{2}}e^{-i\epsilon}S^{(2)}_{m_{1},m_{2},n_{1}-1,n_{2}-1},$
$\displaystyle\partial_{\eta}\langle m|V^{(3)}|n\rangle=\langle
m|iV^{(3)}\frac{\sqrt{\hbar}}{3\times 2^{3/2}}(a+a^{\dagger})^{3}|n\rangle$
(205) $\displaystyle\quad=i\frac{\sqrt{\hbar}}{3\times
2^{3/2}}\left(\sqrt{n(n-1)(n-2)}V^{(3)}_{m,n-3}+3n^{\tfrac{3}{2}}V^{(3)}_{m,n-1}+3(n+1)^{\tfrac{3}{2}}V^{(3)}_{m,n+1}+\sqrt{(n+1)(n+2)(n+3)}V^{(3)}_{m,n+3}\right).$

Gate | Parameter $u=xe^{i\epsilon}$ | Generator $h$ | Value of $s$
---|---|---|---
$D(\gamma)$ | $\gamma$ | $a^{\dagger}_{1}$ | 1
$S(\zeta)$ | $\zeta=re^{i\delta}$ | $-\tfrac{1}{2}a_{1}^{\dagger 2}$ | $\tfrac{1}{2}$
$S^{(2)}(\zeta)$ | $\zeta=re^{i\delta}$ | $a_{1}^{\dagger}a_{2}^{\dagger}$ | 1
$B(\theta,\varphi)$ | $\theta e^{i\varphi}$ | $a_{1}a^{\dagger}_{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. (203) 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\epsilon}),G(x^{\prime}e^{i\epsilon})]=0$ but the same
commutativity does not hold for different values of the phase,
$[G(xe^{i\epsilon}),G(xe^{i\epsilon^{\prime}})]\neq 0$. To make progress we
note that for any of the single complex-parameter gates one can obtain the
following simple decomposition:

$\displaystyle
G(xe^{i\epsilon})=R_{1}(s\epsilon)G(x)R_{1}^{\dagger}(s\epsilon),$ (206)

where $R_{1}(s\epsilon)=e^{is\epsilon a_{1}^{\dagger}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,

$\displaystyle\partial_{\epsilon}\langle\bm{m}|G|\bm{n}\rangle=$
$\displaystyle\langle\bm{m}|(\partial_{\epsilon}R_{1}(s\epsilon))G(x)R_{1}^{\dagger}(\epsilon)+R_{1}(s\epsilon)G(x)(\partial_{\epsilon}R_{1}^{\dagger}(s\epsilon))|\bm{n}\rangle$
$\displaystyle=$
$\displaystyle\langle\bm{m}|(isa_{1}^{\dagger}a_{1})R_{1}(s\epsilon)G(x)R_{1}^{\dagger}(s\epsilon)+R_{1}(s\epsilon)G(x)R_{1}^{\dagger}(s\epsilon)(-isa_{1}^{\dagger}a_{1})|\bm{n}\rangle$
$\displaystyle=$ $\displaystyle is(m_{1}-n_{1})\langle\bm{m}|G|\bm{n}\rangle.$
(207)

## References

  * Schuld et al. [2019] Maria Schuld, Ville Bergholm, Christian Gogolin, Josh Izaac, and Nathan Killoran.  Evaluating analytic gradients on quantum hardware.  _Phys. Rev. A_ , 99(3):032331, 2019.  doi: 10.1103/PhysRevA.99.032331.
  * Benedetti et al. [2019] Marcello Benedetti, Erika Lloyd, Stefan Sack, and Mattia Fiorentini.  Parameterized quantum circuits as machine learning models.  _Quantum Sci. Technol._ , 4(4):043001, 2019.  doi: 10.1088/2058-9565/ab4eb5.
  * Dunjko et al. [2016] 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.
  * Kerenidis et al. [2019] Iordanis Kerenidis, Jonas Landman, Alessandro Luongo, and Anupam Prakash.  q-means: A quantum algorithm for unsupervised machine learning.  In _Advances in Neural Information Processing Systems_ , pages 4136–4146, 2019.
  * Cerezo et al. [2020] Marco Cerezo, Alexander Poremba, Lukasz Cincio, and Patrick J Coles.  Variational quantum fidelity estimation.  _Quantum_ , 4:248, 2020.  doi: 10.22331/q-2020-03-26-248.
  * LaRose et al. [2019] Ryan LaRose, Arkin Tikku, Étude 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.
  * Bravo-Prieto et al. [2019] 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.
  * Arrazola et al. [2019] Juan Miguel Arrazola, Thomas R Bromley, Josh Izaac, Casey R Myers, Kamil Brádler, and Nathan Killoran.  Machine learning method for state preparation and gate synthesis on photonic quantum computers.  _Quantum Sci. Technol._ , 4(2):024004, 2019.  doi: 10.1088/2058-9565/aaf59e.
  * Killoran et al. [2019a] Nathan Killoran, Thomas R Bromley, Juan Miguel Arrazola, Maria Schuld, Nicolás Quesada, and Seth Lloyd.  Continuous-variable quantum neural networks.  _Phys. Rev. Research_ , 1(3):033063, 2019a.  doi: 10.1103/PhysRevResearch.1.033063.
  * Killoran et al. [2019b] Nathan Killoran, Josh Izaac, Nicolás Quesada, Ville Bergholm, Matthew Amy, and Christian Weedbrook.  Strawberry Fields: A software platform for photonic quantum computing.  _Quantum_ , 3:129, 2019b.  doi: 10.22331/q-2019-03-11-129.
  * Steinbrecher et al. [2019] 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.
  * Quesada et al. [2019] N Quesada, LG Helt, J Izaac, JM Arrazola, R Shahrokhshahi, CR Myers, and KK Sabapathy.  Simulating realistic non-Gaussian state preparation.  _Phys. Rev. A_ , 100(2):022341, 2019.  doi: 10.1103/PhysRevA.100.022341.
  * Romero et al. [2017] Jonathan Romero, Jonathan P Olson, and Alan Aspuru-Guzik.  Quantum autoencoders for efficient compression of quantum data.  _Quantum Sci. Technol._ , 2(4):045001, 2017.  doi: 10.1088/2058-9565/aa8072.
  * Schuld et al. [2020] Maria Schuld, Alex Bocharov, Krysta M Svore, and Nathan Wiebe.  Circuit-centric quantum classifiers.  _Phys. Rev. A_ , 101(3):032308, 2020.  doi: 10.1103/PhysRevA.101.032308.
  * Schuld and Killoran [2019] Maria Schuld and Nathan Killoran.  Quantum machine learning in feature Hilbert spaces.  _Phys. Rev. Lett._ , 122(4):040504, 2019.  doi: 10.1103/PhysRevLett.122.040504.
  * Havlíček et al. [2019] Vojtěch Havlíček, Antonio D Córcoles, Kristan Temme, Aram W Harrow, Abhinav Kandala, Jerry M Chow, and Jay M Gambetta.  Supervised learning with quantum-enhanced feature spaces.  _Nature_ , 567(7747):209–212, 2019.  doi: 10.1038/s41586-019-0980-2.
  * O’Malley et al. [2016] 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/PhysRevX.6.031007.
  * Peruzzo et al. [2014] Alberto Peruzzo, Jarrod McClean, Peter Shadbolt, Man-Hong Yung, Xiao-Qi Zhou, Peter J Love, Alán Aspuru-Guzik, and Jeremy L O’Brien.  A variational eigenvalue solver on a photonic quantum processor.  _Nature Commun._ , 5:4213, 2014.  doi: 10.1038/ncomms5213.
  * Politi et al. [2009] Alberto Politi, Jonathan CF Matthews, Mark G Thompson, and Jeremy L O’Brien.  Integrated quantum photonics.  _IEEE J. Sel. Top. Quantum Electron._ , 15(6):1673–1684, 2009.  doi: 10.1109/JSTQE.2009.2026060.
  * Zhang et al. [2020] Y Zhang, M Menotti, K Tan, VD Vaidya, DH Mahler, L Zatti, M Liscidini, B Morrison, and Z Vernon.  Single-mode quadrature squeezing using dual-pump four-wave mixing in an integrated nanophotonic device.  _arXiv preprint arXiv:2001.09474_ , 2020.
  * Vaidya et al. [2020] VD Vaidya, B Morrison, LG Helt, R Shahrokshahi, DH Mahler, MJ Collins, K Tan, J Lavoie, A Repingon, M Menotti, et al.  Broadband quadrature-squeezed vacuum and nonclassical photon number correlations from a nanophotonic device.  _Sci. Adv._ , 6(39):eaba9186, 2020.  doi: 10.1126/sciadv.aba9186.
  * Lloyd and Braunstein [1999] Seth Lloyd and Samuel L Braunstein.  Quantum computation over continuous variables.  In _Quantum information with continuous variables_ , pages 9–17. Springer, 1999.  doi: 10.1007/978-94-015-1258-9˙2.
  * Serafini [2017] Alessio Serafini.  _Quantum continuous variables: a primer of theoretical methods_.  CRC Press, 2017.
  * Weedbrook et al. [2012] Christian Weedbrook, Stefano Pirandola, Raúl García-Patrón, Nicolas J Cerf, Timothy C Ralph, Jeffrey H Shapiro, and Seth Lloyd.  Gaussian quantum information.  _Rev. Mod. Phys._ , 84(2):621, 2012.  doi: 10.1103/RevModPhys.84.621.
  * Barnett and Radmore [2002] Stephen Barnett and Paul M Radmore.  _Methods in theoretical quantum optics_ , volume 15.  Oxford University Press, 2002.
  * Cahill and Glauber [1969] Kevin E Cahill and Roy J Glauber.  Ordered expansions in boson amplitude operators.  _Phys. Rev._ , 177(5):1857, 1969.  doi: 10.1103/PhysRev.177.1857.
  * Král [1990] P Král.  Displaced and squeezed Fock states.  _J. Mod. Opt._ , 37(5):889–917, 1990.  doi: 10.1080/09500349014550941.
  * Ma and Rhodes [1990] Xin Ma and William Rhodes.  Multimode squeeze operators and squeezed states.  _Phys. Rev. A_ , 41(9):4625, 1990.  doi: 10.1103/PhysRevA.41.4625.
  * Quesada [2015] N. Quesada.  _Very Nonlinear Quantum Optics_.  PhD thesis, University of Toronto, 2015.
  * Dhand et al. [2015] 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.
  * Doktorov et al. [1977] EV Doktorov, IA Malkin, and VI Man’ko.  Dynamical symmetry of vibronic transitions in polyatomic molecules and the Franck-Condon principle.  _J. Mol. Spectrosc_ , 64(2):302–326, 1977.  doi: 10.1016/0022-2852(75)90199-X.
  * Gruner and Brumer [1987] Daniel Gruner and Paul Brumer.  Efficient evaluation of harmonic polyatomic Franck-Condon factors.  _Chem. Phys. Lett._ , 138(4):310–314, 1987.  doi: 10.1016/0009-2614(87)80389-5.
  * Berger et al. [1998] 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.
  * Mozhayskiy et al. [2016] Vadim Mozhayskiy, Samer Gozem, and Anna I. Krylov.  ezspectrum v3.0.  http://iopenshell.usc.edu/downloads/, 2016.
  * Huh [2011] Joonsuk Huh.  _Unified description of vibronic transitions with coherent states_.  PhD thesis, Johann Wolfgang Goethe-Universität in Frankfurt am Main, 2011.
  * Rabidoux et al. [2016] 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.
  * Moshinsky and Quesne [1971] Marcos Moshinsky and Christiane Quesne.  Linear canonical transformations and their unitary representations.  _J. Math. Phys._ , 12(8):1772–1780, 1971.  doi: 10.1063/1.1665805.
  * Wolf [1974] Kurt Bernardo Wolf.  Canonical transforms. i. complex linear transforms.  _J. Math. Phys._ , 15(8):1295–1301, 1974.  doi: 10.1063/1.1666811.
  * Kramer et al. [1975] 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.
  * Gupt et al. [2019] Brajesh Gupt, Josh Izaac, and Nicolás Quesada.  The Walrus: a library for the calculation of hafnians, Hermite polynomials and Gaussian boson sampling.  _J. Open Source Softw._ , 4(44):1705, 2019.  doi: 10.21105/joss.01705.
  * Björklund et al. [2019] Andreas Björklund, Brajesh Gupt, and Nicolás Quesada.  A faster hafnian formula for complex matrices and its benchmarking on a supercomputer.  _ACM J. Exp. Algorithmics_ , 24(1):11, 2019.  doi: 10.1145/3325111.
  * Quesada [2019] Nicolás Quesada.  Franck-Condon factors by counting perfect matchings of graphs with loops.  _J. Chem. Phys._ , 150(16):164113, 2019.  doi: 10.1063/1.5086387.
  * Banchi et al. [2020] Leonardo Banchi, Nicolás Quesada, and Juan Miguel Arrazola.  Training Gaussian boson sampling distributions.  _Phys. Rev. A_ , 102:012417, 2020.  doi: 10.1103/PhysRevA.102.012417.
  * Bloch and Messiah [1962] 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.
  * Miatto [2019] Filippo M Miatto.  Recursive multivariate derivatives of $e^{f(x_{1},\dots,x_{n})}$ of arbitrary order.  _arXiv preprint arXiv:1911.11722_ , 2019.
  * Cariolaro and Pierobon [2016] 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.
  * Klimov and Chumakov [2009] Andrei B Klimov and Sergei M Chumakov.  _A group-theoretical approach to quantum optics: models of atom-field interactions_.  John Wiley & Sons, 2009.
  * Hunger [2007] Raphael Hunger.  An introduction to complex differentials and complex differentiability.  Technical report, Munich University of Technology, Inst. for Circuit Theory and Signal Processing, 2007.
  * Guo and Poletti [2020] Chu Guo and Dario Poletti.  A scheme for automatic differentiation of complex loss functions.  _arXiv preprint arXiv:2003.04295_ , 2020.
  * Gottesman et al. [2001] Daniel Gottesman, Alexei Kitaev, and John Preskill.  Encoding a qubit in an oscillator.  _Phys. Rev. A_ , 64(1):012310, 2001.  doi: 10.1103/PhysRevA.64.012310.
  * Huh [2020] Joonsuk Huh.  Multimode Bogoliubov transformation and Husimi’s q-function.  _arXiv preprint arXiv:2004.05766_ , 2020.
  * Walt et al. [2011] Stéfan van der Walt, S Chris Colbert, and Gael Varoquaux.  The NumPy array: a structure for efficient numerical computation.  _Comput. Sci. Eng._ , 13(2):22–30, 2011.  doi: 10.1109/MCSE.2011.37.
  * Virtanen et al. [2020] Pauli Virtanen, Ralf Gommers, Travis E Oliphant, Matt Haberland, Tyler Reddy, David Cournapeau, 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.
  * Kluyver et al. [2016] Thomas Kluyver, Benjamin Ragan-Kelley, Fernando Pérez, Brian E Granger, Matthias Bussonnier, Jonathan Frederic, Kyle Kelley, Jessica B Hamrick, Jason Grout, Sylvain Corlay, et al.  Jupyter notebooks-a publishing format for reproducible computational workflows.  In _ELPUB_ , pages 87–90, 2016.  doi: 10.3233/978-1-61499-649-1-87.
  * Hunter [2007] John D Hunter.  Matplotlib: A 2d graphics environment.  _Comput. Sci, Eng._ , 9(3):90–95, 2007.  doi: 10.1109/MCSE.2007.55.
  * Lam et al. [2015] Siu Kwan Lam, Antoine Pitrou, and Stanley Seibert.  Numba: A LLVM-based Python JIT compiler.  In _Proceedings of the Second Workshop on the LLVM Compiler Infrastructure in HPC_ , pages 1–6, 2015.  doi: 10.1145/2833157.2833162.
  * Berkowitz and Garner [1970] 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.
  * Kok and Braunstein [2001] Pieter Kok and Samuel L Braunstein.  Multi-dimensional Hermite polynomials in quantum optics.  _J. Phys. A: Math. Gen._ , 34(31):6185, 2001.  doi: 10.1088/0305-4470/34/31/312.
  * Mizrahi [1975] Maurice M Mizrahi.  Generalized Hermite polynomials.  _J. Comput. Appl. Math._ , 1(3):137–140, 1975\.  doi: 10.1016/0771-050X(75)90031-5.
  * Braunstein and Van Loock [2005] 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.
  * Kalajdzievski et al. [2018] Timjan Kalajdzievski, Christian Weedbrook, and Patrick Rebentrost.  Continuous-variable gate decomposition for the Bose-Hubbard model.  _Phys. Rev. A_ , 97(6):062311, 2018.  doi: 10.1103/PhysRevA.97.062311.
  * Kalajdzievski and Arrazola [2019] 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.
  * Yanagimoto et al. [2020] Ryotatsu Yanagimoto, Tatsuhiro Onodera, Edwin Ng, Logan G Wright, Peter L McMahon, and Hideo Mabuchi.  Engineering a Kerr-based deterministic cubic phase gate via Gaussian operations.  _Phys. Rev. Lett._ , 124(24):240503, 2020.  doi: 10.1103/PhysRevLett.124.240503.
  * Gu et al. [2009] 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.
  * Marshall et al. [2015a] 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, 2015a.  doi: 10.1103/PhysRevA.91.032321.
  * Sabapathy and Weedbrook [2018] Krishna Kumar Sabapathy and Christian Weedbrook.  ON states as resource units for universal quantum computation with photonic architectures.  _Phys. Rev. A_ , 97(6):062315, 2018.  doi: 10.1103/PhysRevA.97.062315.
  * Sabapathy et al. [2019] Krishna Kumar Sabapathy, Haoyu Qi, Josh Izaac, and Christian Weedbrook.  Production of photonic universal quantum gates enhanced by machine learning.  _Phys. Rev. A_ , 100(1):012326, 2019.  doi: 10.1103/PhysRevA.100.012326.
  * Abramowitz and Stegun [1948] Milton Abramowitz and Irene A Stegun.  _Handbook of mathematical functions with formulas, graphs, and mathematical tables_ , volume 55.  US Government printing office, 1948.
  * Olivier and Manuel [2010] Vallee Olivier and Soares Manuel.  _Airy functions and applications to physics_.  World Scientific Publishing Company, 2010.
  * Marshall et al. [2015b] Kevin Marshall, Raphael Pooser, George Siopsis, and Christian Weedbrook.  Quantum simulation of quantum field theory using continuous variables.  _Phys. Rev. A_ , 92:063825, 2015b.  doi: 10.1103/PhysRevA.92.063825.
  * Zee [2010] Anthony Zee.  _Quantum field theory in a nutshell_ , volume 7.  Princeton university press, 2010.
  * García-Álvarez et al. [2020] Laura García-Álvarez, Cameron Calcluth, Alessandro Ferraro, and Giulia Ferrini.  Efficient simulatability of continuous-variable circuits with large Wigner negativity.  _arXiv preprint arXiv:2005.12026_ , 2020.
  * Gradshteyn and Ryzhik [1980] Izrail Solomonovich Gradshteyn and Iosif Moiseevich Ryzhik.  _Table of integrals, series, and products_.  Academic press, 1980.