# Photonic counterdiabatic quantum optimization algorithm

Pranav Chandarana<EMAIL_ADDRESS>Department of Physical
Chemistry, University of the Basque Country UPV/EHU, Apartado 644, 48080
Bilbao, Spain  EHU Quantum Center, University of the Basque Country UPV/EHU,
Barrio Sarriena, s/n, 48940 Leioa, Spain  Koushik Paul  Department of Physical
Chemistry, University of the Basque Country UPV/EHU, Apartado 644, 48080
Bilbao, Spain  EHU Quantum Center, University of the Basque Country UPV/EHU,
Barrio Sarriena, s/n, 48940 Leioa, Spain  Mikel Garcia-de-Andoin  Department
of Physical Chemistry, University of the Basque Country UPV/EHU, Apartado 644,
48080 Bilbao, Spain  EHU Quantum Center, University of the Basque Country
UPV/EHU, Barrio Sarriena, s/n, 48940 Leioa, Spain  TECNALIA, Basque Research
and Technology Alliance (BRTA), 48160 Derio, Spain  Yue Ban  TECNALIA, Basque
Research and Technology Alliance (BRTA), 48160 Derio, Spain  Mikel Sanz
Department of Physical Chemistry, University of the Basque Country UPV/EHU,
Apartado 644, 48080 Bilbao, Spain  EHU Quantum Center, University of the
Basque Country UPV/EHU, Barrio Sarriena, s/n, 48940 Leioa, Spain  IKERBASQUE,
Basque Foundation for Science, 48009 Bilbao, Spain  Basque Center for Applied
Mathematics BCAM, 48009 Bilbao, Spain  Xi Chen<EMAIL_ADDRESS>Department of Physical Chemistry, University of the Basque Country UPV/EHU,
Apartado 644, 48080 Bilbao, Spain  EHU Quantum Center, University of the
Basque Country UPV/EHU, Barrio Sarriena, s/n, 48940 Leioa, Spain

###### Abstract

We propose a hybrid quantum-classical approximate optimization algorithm for
photonic quantum computing, specifically tailored for addressing continuous-
variable optimization problems. Inspired by counterdiabatic protocols, our
algorithm significantly reduces the required quantum operations for
optimization as compared to adiabatic protocols. This reduction enables us to
tackle non-convex continuous optimization and countably infinite integer
programming within the near-term era of quantum computing. Through
comprehensive benchmarking, we demonstrate that our approach outperforms
existing state-of-the-art hybrid adiabatic quantum algorithms in terms of
convergence and implementability. Remarkably, our algorithm offers a practical
and accessible experimental realization, bypassing the need for high-order
operations and overcoming experimental constraints. We conduct proof-of-
principle experiments on an eight-mode nanophotonic quantum chip, successfully
showcasing the feasibility and potential impact of the algorithm.

##  I Introduction

Harnessing usefulness from current noisy intermediate-scale quantum (NISQ) [1]
computers has emerged as the main objective of the quantum computing
community. Variational quantum algorithms (VQAs) are the leading candidates to
achieve this goal, making use of the limited quantum resources that existing
NISQ computers offer [2]. These hybrid algorithms aim to solve computationally
demanding tasks with enhanced efficiency by synergistically combining the
computational power of classical systems. In the past few years, VQAs have
demonstrated significant potential in addressing numerous challenges in
contemporary science, such as problems involving many-body quantum
Hamiltonians [3, 4], quantum chemistry [5, 6, 7], combinatorial optimization
[8, 9], and others [10, 11].

The fundamental structure of any VQA typically consists of three essential
elements; encoding, processing, and decoding. The first step involves
transforming a logical optimization problem into a Hamiltonian, encoding the
solution in its ground state. The processing phase refers to designing a
parameterized quantum circuit or the circuit ansatz to manipulate quantum
states, whereas the decoding step includes a measurement scheme to evaluate
the cost function. This cost function is minimized using classical
optimization routines that eventually produce an optimal quantum state.
Currently, most research efforts are focused on encoding discrete optimization
problems using qubit-based approaches, which are well-suited for
implementation in superconducting circuits and trapped ions.

In contrast, photonic quantum computing (PQC) incorporates the continuous
variable (CV) formalism. This allows quantum information to be encoded in the
quadrature amplitudes of an electromagnetic field [12] known as qumodes. This
encoding provides benefits in representing continuous optimization problems
that are expensive to encode with qubits. Recently, there has been a growing
interest in developing CV quantum algorithms for solving various problems in
the PQC paradigm. For instance, a CV-based quantum approximate optimization
algorithm (CV-QAOA) has recently been proposed and benchmarked with the
minimization of the non-convex Styblinski-Tang function [13]. Likewise, a CV
adiabatic quantum algorithm was also proposed that investigated mixed-integer
programming problems using Fock encoding [14]. In addition, efforts have been
made to encode graph problems, imaginary time evolution for quantum field
theories, Grover search and instantaneous quantum polytime circuits over
continuous spaces [15, 16, 17, 18]. Apart from optimization, the CV regime is
also utilized as a tool for error-correcting codes [19, 20] and for quantum
state learning [21].

One of the challenges encountered while implementing adiabatic algorithms
using the CV approach is the inherent bottleneck arising from the quadratic
nature of Gaussian operations in phase space. Generally, complex optimization
problems require non-Gaussian operations and intricate gate decomposition for
effective time evolution. The experimental feasibility of these algorithms in
solving high-degree problems is limited, as recent experiments have only been
able to solve simple quadratic function of $x$ in the case of CV-QAOA [22].
Therefore, there is a pressing need for algorithms that operate within low
degrees while still enabling experimental exploration of complex high-degree
problems using near-term photonic devices.

In this article, we propose the photonic counterdiabatic quantum optimization
(PCQO) algorithm to address this challenge. PCQO is a hybrid quantum-classical
algorithm designed to solve problems suitable for currently available photonic
devices by utilizing a circuit ansatz and a classical optimization routine.
This circuit ansatz is designed from a pool of Gaussian and non-Gaussian
operations that are obtained by drawing inspiration from counterdiabatic (CD)
protocols [23]. These CD protocols, known as shortcuts-to-adiabaticity methods
[24, 25], accelerate the adiabatic process and circumvent the typically slow
evolution mandated by the adiabatic theorem [26]. Previous applications of
these methods have demonstrated substantial improvements in QAOA and digitized
adiabatic evolution [27, 28, 29, 30, 31].

We investigate the performance of this algorithm for (a) phase-space encoding,
which encodes information in the $\hat{x}$ quadrature, and (b) Hilbert-space
encoding, which represents information using Fock states $\ket{n}$. The former
includes classical non-convex continuous optimization problems [32, 33] and
the latter includes integer programming problems [34]. Additionally, we show
that PCQO outperforms state-of-the-art quantum algorithms such as CV-QAOA in
terms of performance and implementability. Lastly, we provide considerations
for implementing PCQO with NISQ devices and show a proof-of-principle
experiment that solves a simple problem in an eight-mode photonic quantum
computer.

This article is structured as follows. Section. II provides a discussion of
variational quantum algorithms in the context of qubits and qumodes. In
Section. III, we present a detailed description of the PCQO algorithm and
benchmarked the algorithm in Section. IV for non-convex continuous
optimization problems using phase-space encoding. In Section. V, we apply PCQO
to countably infinite discrete optimization like unbounded knapsack problem
and Maxclique problem using Hilbert-space encoding. Section. VII provides a
proof-of-principle experiment to further establish the feasibility of our
algorithm in the near-term era. Finally, we conclude in Section. VIII and
discuss possible future directions.

##  II Variational quantum algorithms

In the qubit-based regime, the encoding is generally done by considering a
physical system such as a molecule or a spin chain. A problem Hamiltonian
$H_{p}$ corresponding to this system is found such that its ground state
entails the information of the solution. The processing phase employs a
circuit ansatz comprising various parameterized gates. The design of this
ansatz is critical to the performance of the VQA, as it directly affects the
energy landscape and consequently, the convergence and success rate of the
optimization process.

Circuit ansatzes can be broadly divided into two categories: problem-inspired
and hardware efficient. In problem-inspired ansatzes, the parameterized
unitary is constructed by taking information from the $H_{p}$ which usually
corresponds to time evolutions. For example, in the quantum approximate
optimization algorithm (QAOA), the parameterized unitary is given by
$U(\alpha)=e^{-i\alpha H_{p}}$, where $\alpha$ is an optimizable parameter
[35]. On the other hand, the hardware-efficient ansatzes are specifically
designed to take into account the connectivity constraints of the underlying
quantum hardware [7]. Both problem-inspired and hardware-efficient ansatzes
have distinct roles in VQA design, providing varying trade-offs between
performance and feasibility in practical implementations. The choice of ansatz
depends on factors like the problem characteristics, available resources, and
desired accuracy.

Ideally, problem-inspired ansatzes should be prioritized over hardware-
efficient ansatzes because they narrow down the search space, enhancing
trainability. However, problem-inspired ansatzes often result in increased
circuit depths and unfavorable connections. Consequently, there is significant
interest in developing algorithms containing problem-inspired ansatzes that
can be implemented on near-term devices [36].

Lastly, the decoding step involves measurements, typically in the
computational basis, to evaluate a cost function. This cost function maps the
optimizable parameters to real numbers. In many cases, the cost function
corresponds to the expectation value of the $H_{p}$. However, alternative
metrics such as fidelity or conditional value at risk [37] can also be
considered.

In the qumode-based regime, the components required for designing a VQA differ
intrinsically from those in the qubit-based regime. In the encoding step, we
have the flexibility to choose between two formulations: the phase space
formulation and the Hilbert space formulation. These correspond to the wave-
like and particle-like nature of light, respectively. In the phase space
picture, the state of a single qumode is represented using
$(\hat{x},\hat{p})$, which are the position and momentum operators,
respectively. The problem Hamiltonian can be expressed as
$H_{p}=F(\hat{x},\hat{p})$. On the other hand, qumode states can also be
represented in a countably infinite Hilbert space spanned by Fock states
$\ket{n_{i=0,1,\dots}}$. Consequently, the problem Hamiltonian can be written
as $H_{p}=F(\hat{n})$.

In the processing stage, the overall recipe remains the same, but PQC involves
different operations designed for qumodes. These operations can be categorized
into Gaussian operations and non-Gaussian operations. As the name suggests,
Gaussian operations map Gaussian states to themselves and are generated by
quadratic operations in $\hat{x}$ and $\hat{p}$. Single-mode Gaussian gates
include phase shifts, displacement, and squeezing. Two-mode gates include
beamsplitters. A combination of these gates can be used to implement gates
such as quadratic phase gates, controlled-phase gates, Mach-Zehender
interferometers, etc. On the other hand, non-Gaussian operations do not
preserve the Gaussian nature of the quadratic states. These include the cubic
phase gate, the Kerr gate, and the cross-Kerr gate. Expressions of all these
gates are given in Table 1. It has been shown that all these Gaussian
transformations combined with any single non-Gaussian transformation make a
universal gate set for PQC [38].

Table 1: Currently available qumode operations as gates. In each gate, the argument shows adjustable parameters. $i$ and $j$ denote two arbitrary qumodes. Gate |  Expression |  Type
---|---|---
Rotation or Phase-shift |  $R(\phi)=\exp[i\phi\hat{n}]$ |  Gaussian
Displacement |  $D(\alpha)=\exp[\alpha\hat{a}-\alpha^{*}\hat{a}^{\dagger}]$ |  Gaussian
Squeezing |  $S(r,\phi)=\exp\left[\frac{r}{2}(e^{-i\phi}\hat{a}^{2}-e^{i\phi}\hat{a}^{\dagger 2})\right]$ |  Gaussian
Beamsplitter |  $BS(\theta,\phi)=\exp\left[\theta\left(e^{i\phi}\hat{a}_{i}\hat{a}^{\dagger}_{j}-e^{-i\phi}\hat{a}^{\dagger}_{i}\hat{a}_{j}\right)\right]$ |  Gaussian
Quadratic Phase |  $P(s)=\exp\left[i\frac{s}{2\hbar}x^{2}\right]$ |  Gaussian (Decomposable)
Controlled-Phase |  $CZ(s)=\exp\left[is\frac{\hat{x}_{i}\hat{x}_{j}}{\hbar}\right]$ |  Gaussian (Decomposable)
Two-mode squeezing |  $S_{2}(z)=\exp\left[\ z\hat{a}^{\dagger}_{i}\hat{a}^{\dagger}_{j}-z^{*}\hat{a}_{i}\hat{a}_{j}\right]$ |  Gaussian (Decomposable)
Cubic Phase |  $V(\gamma)=\exp\left[i\frac{\gamma}{3\hbar}\hat{x}^{3}\right]$ |  Non-Gaussian
Kerr |  $K(\kappa)=\exp\left[i\kappa\hat{n}^{2}\right]$ |  Non-Gaussian
Cross-Kerr |  $CK(\kappa)=\exp\left[i\kappa\hat{n}_{i}\hat{n}_{j}\right]$ |  Non-Gaussian

In the decoding stage, the measurement scheme can be homodyne, heterodyne, or
photon number-resolving measurements. The choice of measurement scheme is
determined by the nature of the encoding and the specific variables of
interest, either position or number operators.

Based on the gates available in qumode-based architectures, one of the
challenges in implementing problem-inspired algorithms is to implement
problems represented by polynomials of degree $d\geq 3$ as all the Gaussian
operations are quadratic in terms of phase-space quadratures. Even for $d=3$,
implementing a cubic-phase gate is required, which is a non-Gaussian
operation. For $d>3$, we have to decompose the operations with higher degrees
in terms of available gates [39]. This challenge applies to both phase space
and Hilbert space encodings. Consequently, there is a need for an algorithm
that can efficiently handle functions of any degree using a low number of
gates, making it feasible for implementation on currently available photonic
devices. In the following sections, we explore how counterdiabatic protocols
can address this requirement and enable the efficient implementation of high-
degree functions in PQC. We accomplish this by providing examples illustrating
how both formulations can be used to encode different types of problems and
demonstrate the effectiveness of our algorithm in both cases.

##  III Photonic counterdiabatic quantum optimization

Figure 1: A schematic diagram illustrating the $p=1$ photonic counterdiabatic
quantum optimization (PCQO) algorithm. (a) Encoding phase: A logical problem,
denoted as $H_{p}$, is encoded into $F(\mathbf{n})$ or $F(\mathbf{x})$ based
on the problem type. The mixer $H_{m}$ is introduced satisfying the non-
commutativity relation, enabling the generation of the operator pool using the
nested commutator method. (b) Processing phase: The operators from
$\mathcal{A}=\\{A_{\lambda}^{(2)}\\}$ are exponentiated and employed as a
circuit ansatz $U_{cd}(\boldsymbol{\theta})$ with adjustable parameters
$\boldsymbol{\theta}=\\{\theta_{1},\theta_{2},\ldots,\theta_{Q}\\}$. $Q$ shows
the total number of parameters. The algorithm initiates with random parameter
values and iteratively updates them through classical optimization, aiming to
determine $F(\braket{\mathbf{n}})$ or $F(\braket{\mathbf{x}})$ until
convergence is achieved. (c) Decoding phase: Performing measurements and
extracting solutions from the minimum values, $F_{min}(\braket{\mathbf{n}})$
or $F_{min}(\braket{\mathbf{x}})$, enables the representation of solutions in
the form of the mean photon number $\braket{\mathbf{n}}$ or the mean
quadrature values $\braket{\mathbf{x}}$.

To design the algorithm, we start with an adiabatic quantum Hamiltonian
$H_{a}(t)$, given by

$H_{a}(t)=(1-\lambda(t))H_{m}+\lambda(t)H_{p}$ (1)

where $\lambda(t)$ is a scheduling function such that $\lambda(0)=0$ and
$\lambda(T)=1$ and $T$ is the total evolution time. $H_{p}$ is a problem
Hamiltonian whose ground state we need to find and $H_{m}$ is a mixer
Hamiltonian, such that $[H_{p},H_{m}]\neq 0$. According to the adiabatic
theorem, if the system is prepared in an eigenstate of $H_{m}$, it remains in
the instantaneous eigenstate during the evolution, provided that the evolution
is slow such that $|~{}\dot{\lambda}~{}|\ll 1$. A universal quantum computing
paradigm, called adiabatic quantum computing, was developed based on this
Hamiltonian [40]. This paradigm was extended to digital quantum computing with
experiments on a superconducting circuit using Trotterization [41]. One of the
VQAs that take inspiration from this is QAOA. In QAOA, two non-commuting
unitaries called the mixer term $U_{b}(\beta)$ and the Hamiltonian term
$U_{c}(\gamma)$ are applied iteratively for $p$ layers to an initial state
$\ket{\psi_{0}}$, which is the ground state of $H_{m}$. $\gamma$ and $\beta$
are optimizable parameters, and the parameterized unitary looks like

$U(\vec{\gamma},\vec{\beta})=\prod_{k=1}^{p}U_{b}(\beta_{k})U_{c}(\gamma_{k}).$
(2)

where, $U_{b}(\beta)=e^{-i\beta H_{m}}$, $U_{c}(\gamma)=e^{-i\gamma H_{p}}$.
The mixer Hamiltonian $H_{m}$ can take various forms depending on the specific
requirements of the problem being solved [42]. The difficulty that comes with
adiabatic algorithms is the requirement of deep circuits to satisfy the slow
evolution condition. If the evolution is not sufficiently slow, there will be
diabatic transitions that will reduce the probability of finding the system in
the ground state of $H_{p}$. In order to get fast evolutions, CD protocols are
implemented. Here, the task is to add velocity-dependent terms
$A_{\lambda}^{(l)}$ to minimize the non-adiabatic transitions. This results in
the Hamiltonian $H_{cd}$ given by

$H_{cd}(t)=H_{a}(t)+\dot{\lambda}(t)A_{\lambda}^{(l)}.$ (3)

The calculation of the exact CD term requires full spectral information [43,
44]. This information may not always be available hence approximate CD terms
can be used instead [45]. One of the ways to obtain these terms is to utilize
adiabatic gauge potentials that can be calculated using the nested commutator
method [46] given by

$A_{\lambda}^{(l)}=i\sum_{k=1}^{l}\alpha_{k}(t)\underbrace{[H_{a},[H_{a},\dots,[H_{a},}_{2k-1}\partial_{\lambda}H_{a}]]],$
(4)

where $l$ is the order of expansion and $\alpha_{k}(t)$ are the CD
coefficients that need to be optimized. This can be done analytically [46] or
by using variational circuits [47].

Now, let us assume that there exists a $\lambda$ such that
$|~{}\dot{\lambda}~{}|\gg 1$ at the beginning of the evolution. For this
condition, the $H_{a}$ term from Eq. (3) can be neglected because only
contributions that will dominate will be from the CD terms
$A_{\lambda}^{(l)}$. Thus, we can devise a parameterized unitary evolution
that will look like

$U(\boldsymbol{\theta})=\prod_{k=1}^{p}U_{cd}(\boldsymbol{\theta}_{k}),$ (5)

where $U_{cd}(\boldsymbol{\theta})=e^{-i\boldsymbol{\theta}\mathcal{A}}$.
$\mathcal{A}$ is a set of operators obtained by Eq. (4) with low order of
expansion and $\boldsymbol{\theta}_{k}$ is a set of tunable parameters. This
method has been shown to reduce the circuit depths drastically as compared to
other adiabatic algorithms in qubit-based technologies [48].

The adiabatic theorem is constructed under the assumption that both $H_{m}$
and $H_{p}$ have discrete spectra and there exists a finite minimum spectral
energy gap for all times $t$ [49]. However, in the CV regime, this is not
always guaranteed since the evolution cannot be bounded with a well-defined
spectral gap. This motivated proposals for equivalent adiabatic theorem for
continuous systems [50]. In any case, we can always discretize the time and
implement the evolution in terms of a gate-based model. The main motivation
for this work is to develop hybrid CD-inspired protocols for CV systems,
specifically PQC. The expected advantages of doing this are two-fold. Firstly,
these methods should result in a reduction in the required number of
operations. Since the operations in photonics are beam-splitters and
interferometers which are imperfect, this reduction should improve the
performance to a great extent. Secondly, since at a finite order $l$, we get a
pool of CD operators from which we can choose suitable operations that are
natively available in the device. This makes the algorithm more flexible for
NISQ devices. As we will see later, this freedom will also allow us to attempt
complicated problems without the need for decompositions. This is crucial
since decompositions are non-trivial and require lots of resources in the PQC
regime.

Now, we devise a CD-inspired algorithm for PQC. To do so, we start by
considering a problem Hamiltonian $H_{p}$. This can be a function of
$\mathbf{x}=(x_{1},\ldots,x_{N})\in\mathbb{R}^{N}$ or a function of
$\mathbf{n}=(n_{1},\ldots,n_{N})\in\mathbb{W}^{N}$ depending upon the
encoding. Then we define a mixer Hamiltonian $H_{m}$ whose ground-state is
easy to prepare and satisfies the condition $[H_{p},H_{m}]\neq 0$. With
$H_{m}$ and $H_{p}$, we now have $H_{a}$ given by Eq. (1). From this, we can
compute $A_{\lambda}^{(l)}$ using Eq. (4) with a specific order $l$ to get a
pool of operators. We heuristically select operators from this pool and
parameterize them. This selection is made based on the requirements of the
problem and the hardware. Starting from a specific initial state, these
parameters can be optimized by classical optimization routines to minimize
$F(\braket{\mathbf{x}})$ or $F(\braket{\mathbf{n}})$. A schematic diagram of
the algorithm is shown in Fig. 1.

The cost function used in this algorithm focuses on minimizing the mean values
of the operators $\braket{\hat{x}}$ or $\braket{\hat{n}}$ instead of
optimizing the expectation value $\braket{F(\mathbf{x})}$ or
$\braket{F(\mathbf{n})}$. This provides a useful approach for hybrid
algorithms seeking approximate solutions. By considering the mean values along
with the optimal state, the cost function takes into account the probability
of obtaining near-optimal states. Moreover, this cost function can exhibit
robustness against noise since it relies on the mean of local operators’
expectation values rather than the specific state itself. This robustness can
be advantageous in practical implementations where noise and imperfections are
inevitable. For example, if we want to solve $F(n)=(n-1)^{2}$, the circuit
ansatz will optimize the parameters such that $\braket{\hat{n}}=1$. Hence,
along with $\ket{n=1}$, $\ket{n=0}$ and $\ket{n=2}$ will also have finite
probabilities.

##  IV Results: Phase-space encoding

We examine the performance of the algorithm by applying it to two non-convex
classical optimization problems $F(\mathbf{x})$, both represented as
polynomials of degree $d$. For both problems, the initial state preparation
involves setting all the qumodes to vacuum states. The mixer Hamiltonian
$H_{m}=\sum_{i}(\hat{p_{i}}-p_{0})^{2}$ was selected, where $p_{0}$ is a
constant. This resembles the kinetic part of the harmonic oscillator and the
potential part would be $F(\mathbf{x})$. This mixer satisfies the non-
commutativity condition with $F(\mathbf{x})$. The next step involves finding
the operator pool $\mathcal{A}$ by evaluating $A_{\lambda}^{(l)}$ which
enables us to choose suitable parameterized gates from $\mathcal{A}$ as a
circuit ansatz. We utilize homodyne measurements to compute
$\braket{\mathbf{x}}$. Since this is a numerical analysis, we have to define a
cutoff dimension of the Fock space which was set to $D=15$ for the toy
function, $D=10$ for the Rosenbrock function due to computational limitations
and we also set $\hbar=2$. We considered two polynomials, one with $d=4$ and
$N=4$ variables and the other with $d=6$ and $N=3$ variables.

We start with the Rosenbrock function [32], defined as

$F(\mathbf{x})=\sum_{i=1}^{N-1}\left[100(x_{i+1}-x_{i}^{2})^{2}+(1-x_{i})^{2}\right]$
(6)

where $\mathbf{x}=(x_{1},\ldots,x_{N})\in\mathbb{R}^{N}$. Due to its non-
convex nature and difficulty in reaching global optima, this function is
widely used as a benchmark for optimization algorithms. We consider $N=4$ case
with global minimum $F_{min}(\mathbf{x}_{opt})=0$ at
$\mathbf{x}_{opt}=(1,1,1,1)$. We limit ourselves to the $A_{\lambda}^{(2)}$
and get
$\mathcal{A}=\\{\hat{p},~{}\hat{x}\hat{p}+\hat{p}\hat{x},~{}\hat{x}^{3},\ldots\\}$.
This will contain many higher-order terms but we heuristically select the
first three terms. These correspond to the $X$ gate, $S_{2}$ gate, and $V$
gate (See Table 1). This choice was made to include single-mode Gaussian, two-
mode Gaussian, and non-Gaussian gates. Since we will perform Gaussian
measurements, these non-Gaussian gates will make sure that the ansatz is not
effectively simulated classically [51]. Since $S_{2}$ is a two-mode gate, we
have to decide the connectivity. For simplicity, we keep this connectivity to
the nearest neighbors but this can be further fine-tuned based on hardware
constraints. Therefore, the circuit ansatz will have $X$ gates applied to all
qumodes, $S_{2}$ gates applied to nearest-neighbor qumodes, and $V$ gates
applied to all qumodes and thus the number of parameters required will be
$Q=2N+(N-1)=3N-1$.

The second benchmark function is a toy problem given by

$\displaystyle F(\mathbf{x})=$
$\displaystyle(x_{1}^{3}+x_{2}^{3}+x_{3}^{3}-1x_{1}+2x_{2}-3x_{3})^{2}+$ (7)
$\displaystyle(-x_{1}+x_{2}+x_{3})^{2}+0.01x_{1}+0.01x_{2}+0.01x_{3}.$

For this problem, the solution is $F_{min}(\mathbf{x}_{opt})=-0.028457$ at
$\mathbf{x}_{opt}=(-1.42212,-0.127017,-1.29723)$. This function is selected
because its degree is much larger than the degree of the operations
implementable in current technologies (See Table. 1). Furthermore, this
function forces correlations between the variables, so the problem cannot be
solved independently for each variable. The linear terms with low coefficients
ensure that the global minima are reached only at a single point. This
function will give the same lower-order operators in the $\mathcal{A}$ pool as
the Rosenbrock function hence we implement the same circuit ansatz as before.
In both cases, the mean values of $F(\braket{\mathbf{x}})$ (termed as energy)
over five random initializations across several iterations for $p=1$ layer of
the ansatz are computed and results are shown in Fig. 2 along with the
standard error. Energy variation corresponding to the best instance is also
plotted. For parameter optimization, we have used the Adam optimizer [52].

Figure 2: Energy (values of $F(\braket{\mathbf{x}})$) profiles for (a) the
Rosenbrock function and (b) the toy function (Eq. 7) with a $p=1$ ansatz.
Dashed lines represent the mean of five random initializations, solid lines
indicate the best outcome among the five instances, and shaded regions show
the standard error. The green line represents the exact energy. The inset
plots display the minimum energies achieved with the five instances.

Fig. 2(a) illustrates the convergence of the mean energy over the first 200
iteration steps out of 500, for the Rosenbrock function, with some iterations
skipped to disregard initial fluctuations. It can be observed that the mean
energy is slightly higher than the exact energy required to solve the problem.
This indicates that the algorithm’s performance is affected by the initial
parameters chosen. This observation is supported by the fact that the standard
error is high in later iterations. In the best-instance run, the algorithm
achieves the exact solution. The inset shows the minimum energy achieved
during the optimization with the five instances. We can see that, for most
instances, the algorithm results in solutions very close to the exact
solution. This demonstrates that if the objective is to obtain approximate
solutions, the algorithm exhibits good performance.

Fig. 2(b) displays the convergence of the mean energy over the first 200
iterations out of 1000, for the toy function in Eq. 7. It can be observed that
the mean energy converges close to the exact energy within the first 200
iterations. Initially, there is a high standard error, but in later
iterations, the algorithm successfully finds approximate solutions. At the
200th iteration, the best instance achieves $F(\mathbf{x})=-0.0037$ while the
mean energy is $F(\mathbf{x})=0.0029$. Both values are sufficiently close to
the exact solution. The inset plot shows the minimum energy obtained during
the 1000 iterations for all instances. The best solution we get is
$F(\mathbf{x})=-0.0280$.

For both cases, we noticed that optimization becomes slower after about 200
iterations. Due to this, a larger number of iterations to achieve more
accurate solutions is required. A possible explanation for this is that the
ansatz reaches a point where the energy landscape becomes flat which results
in low gradient values. Consequently, more iterations are needed to reach the
solution. This issue can potentially be addressed by adding more terms to the
ansatz, although this may introduce a more rugged solution landscape with
potential local minima. Nonetheless, for obtaining approximate solutions, the
algorithm performs well even with random initializations.

The results demonstrate that the proposed PCQO algorithm can obtain good
approximate solutions with just $p=1$ layer of the ansatz. This algorithm can
be extended to handle polynomials of any degree, and it would be intriguing to
explore its performance on higher-degree polynomials. Additionally,
investigating problems with a larger number of variables would provide further
insights into the algorithm’s capabilities and potential applications.

To tackle these problems with qubit-based algorithms, the system size will
depend upon the bit resolution that needs to be achieved. This will lead to
large resource requirements and the cost functions will be many-body Ising
Hamiltonians [53]. Since the requirement to discretize the solution here is
elevated, any arbitrary precision can be reached using the PCQO algorithm with
a linear encoding with respect to the variables. These advantages make PCQO
far more suitable for tackling continuous optimization problems.

Despite the advantages, some challenges need to be addressed. For instance,
clever optimization techniques need to be developed to find the optimal
parameters faster. Strategies to efficiently choose the circuit ansatz from
$\mathcal{A}$ ensuring trainability and expressibility need to be developed
since there is no performance guarantee of the heuristic choice of ansatz.
Nevertheless, PCQO works extremely well for phase-space encoding. In the next
section, we will investigate the performance of the algorithm using Hilbert-
space encoding where the states are countably infinite and discrete.

##  V Results: Hilbert-space encoding

In Hilbert-space encoding, the variables are represented by Fock states
$\ket{n}$. This encoding scheme naturally lends itself to tackling integer
programming problems, which are a class of optimization problems where some or
all variables are constrained to be integers. Integer programming can be
further categorized into three subclasses: linear integer programming, non-
linear integer programming, and mixed-integer programming. The classification
depends on the linearity of the cost function and constraints. Here, we focus
on solving linear integer and non-linear programming problems, but the same
techniques can be readily extended to mixed-integer programming as well.

Similar to the phase-space encoding, the initial state of all the qumodes is
prepared in the vacuum state. The mixer we selected was
$H_{m}=\sum_{i}~{}(x_{i}-x_{0})^{2}+(p_{i}-p_{0})^{2}$ where $x_{0}$ and
$p_{0}$ are constants. The mixer resembles a shifted harmonic oscillator in
both quadratures. Instead of performing homodyne measurements, number-
resolving measurements are used to determine the values of
$\braket{\mathbf{n}}$. These measurements provide information about the number
of photons in each mode, which is essential in the Hilbert space encoding.

###  V.1 Unbounded knapsack problem

Figure 3: Energy evolution over iteration steps for (a) $N=3$ and (b) $N=4$
qumode knapsack problems using a layered ansatz with $p=1$. Dashed lines
represent the mean energy of 5 random initializations, while shaded regions
indicate the corresponding standard error. The solid lines depict the best
energy obtained among the 5 random instances. The green lines indicate the
exact energy minimization. Vertical black lines denote the iteration at which
$10^{-2}$ accuracy relative to the exact energy was reached during
optimization. Insets show the probability distribution of Fock states around
the solution up to the cutoff.

We commence with investigating a small instance of the unbounded knapsack
problem (UKP) [54]. In UKP, we consider a set of different types of items $i$,
each with a value $v_{i}$ and weight $w_{i}$. The objective is to maximize the
total value while ensuring that the total weight of the selected items does
not exceed the capacity $C$ of the knapsack. Unlike the bounded knapsack
problem, the UKP allows for an unlimited number of items of the same type to
be included in the knapsack. If the number of items of type $i$ that can be
included in the knapsack is given by $n_{i}$, then the optimization problem
looks like

$\displaystyle\min_{n_{i}}\ -\sum_{i=1}^{N}v_{i}n_{i}\quad\text{subject
to}\quad\sum_{i=1}^{N}w_{i}n_{i}\leq C,\ n_{i}\geq 0\ \forall i,$ (8)

where $N$ shows the total items. This problem is classified as an integer
linear programming problem, where both the cost function and constraints are
linear functions of $n_{i}$. This can be converted into a minimization problem
of

$F(\mathbf{n})=-\sum_{i=1}^{N}v_{i}n_{i}+\delta\sum_{i=1}^{N}(w_{i}n_{i}-C)^{2}.$
(9)

Here, $\delta$ represents the penalty term. Usually, an auxiliary variable is
added to account for the inequality. However, as this algorithm is aimed at
finding approximate solutions, we add the inequality as a ‘soft’ constraint in
the problem Hamiltonian. We utilized the $p=1$ layered PCQO algorithm to
address the UKP for $N=3$ and $N=4$ qumodes.

Table 2: The data considered for the UKP instances where $N$ are the total items, $v_{i}$ is the value of the item, $w_{i}$ is the weight of the item, $C$ is the total capacity of the knapsack, $F_{min}(\mathbf{n}_{opt})$ shows the optimal cost function at $\mathbf{n}_{opt}$ optimal values.  $N$ | $v$ | $w$ | $C$ | $F_{min}(\mathbf{n}_{opt})$ | $\mathbf{n}_{opt}$
---|---|---|---|---|---
$3$ | $[3,4,1]$ | $[9,5,8]$ | $10$ | $-8$ | $(0,2,0)$
$4$ | $[3,4,1,3]$ | $[2,7,6,6]$ | $10$ | $-15$ | $(5,0,0,0)$

We considered $l=2$ order nested commutator, which will result in
$\mathcal{A}=\\{\hat{x},~{}\hat{p},~{}\hat{x}\hat{p}^{2},~{}\hat{x}\hat{x},\ldots\\}$.
From this, we choose $\hat{x}$ and $\hat{x}\hat{x}$, whose exponentiation will
correspond to $X$ gates and $CZ$ gates respectively. We do not include the
non-Gaussian gates because the measurement is non-Gaussian, which will
restrict the effective classical simulation [51]. Thus the circuit ansatz
consists of $X$ gates applied to all qumodes and $CZ$ gates applied to nearest
neighbor connections. Hence, the number of parameters required will be
$Q=N+(N-1)=2N-1$. The cutoff dimension was chosen as $D=10$, the penalty term
was set to $\delta=4$, and we keep $\hbar=2$. Mean energy (values of
$F(\braket{\mathbf{n}})$) and the best instance as a function of the number of
iterations for both $N=3$ and $N=4$ qumode case are shown in Fig. 3. In both
cases, we skip some iterations to neglect the initial fluctuations. As before,
classical optimization was performed by using the Adam optimizer. The data
considered for the problems are given in Table. 2.

Fig. 3(a) depicts the energy over 200 iterations for the $N=3$ UKP. The mean
energy converges to a value higher than the exact solution. This is due to
initialization, which can cause solutions to fall into local minima. This
claim is supported by the high standard error even in later iterations.
Efficient initialization strategies are crucial for mitigating the dependence
on initial points in VQAs. However, the best instance achieves the exact
energy. The black vertical line represents the instance where $10^{-2}$
accuracy was first attained, around the 170th iteration. This indicates that
with optimal initial parameterization, the algorithm performs exceptionally
well and reaches the exact energy. Fig. 3(b) presents the mean energy
variation across 700 iterations for the $N=4$ UKP. Same as before, we observe
that the mean energy is higher than the exact energy, and the best case
reaches the exact energy. However, in contrast to the previous case, it takes
350 iterations to achieve a tolerance of $10^{-2}$, which is significantly
higher.

It is important to note a subtle caveat in our approach. As we compute
$\langle F(\mathbf{n})\rangle$ as the cost function, the solution space is
continuous because $\langle\hat{n}\rangle$ can take any non-negative values.
Consequently, post-processing techniques are necessary to extract the
solution. This can be accomplished in two ways. Approximating the solution to
the nearest integer or analyzing the probability distribution of the Fock
states to determine the likelihood of obtaining the desired state. The inset
plots in the figure display the probability distribution of states near the
optimal solution for the best instance. In the case of $N=3$, where the
optimal solution is represented by $\mathbf{n}_{opt}=(0,2,0)$, the
probabilities were plotted for states where the 2nd qumode varies from
$\ket{n_{1}}=0$ to $\ket{n_{1}}=D-1$. Similarly, for the $N=4$ case, with the
optimal solution of $\mathbf{n}_{opt}=(5,0,0,0)$. The darkest color in the
plots indicates states that achieve high probabilities, representing the
optimal solutions. However, we can also observe that approximate solutions
have relatively high probabilities also, depicted by lighter colors. This
occurrence arises because our cost function optimizes the mean, allowing sub-
optimal states to have finite probabilities. This characteristic is crucial
for approximate optimization algorithms since their objective is to find
approximate solutions rather than exact ones. For the $N=3$ case, the sum of
probabilities for states $\ket{n}=(0,1,0)$, $\ket{n}=(0,2,0)$, and
$\ket{n}=(0,3,0)$ amounts to approximately 72%. Similarly, in the $N=4$ case,
the sum of probabilities for states $\ket{n}=(4,0,0,0)$, $\ket{n}=(5,0,0,0)$,
and $\ket{n}=(6,0,0,0)$ totals around 45%. So, even when the mean photon
number is minimized, we can extract the solution successfully.

It is straightforward to notice that there is an improvement in terms of the
quantum resources required to encode a problem using qubits. For encoding all
possible solutions to this problem, we need the number of qubits
$N_{\text{qubits}}=\sum_{i=1}^{N}\bigg{\lceil}\log_{2}\left(\bigg{\lfloor}\frac{C}{w_{i}}\bigg{\rfloor}\right)\bigg{\rceil}\geq
N$. This scales worse than linearly with the number of items.
$N_{\text{qubits}}$ corresponds to the number of qubits necessary to encode
trivial solutions in which we only include a single item until we fill the
knapsack, giving an idea of the order of the resources required. On the other
hand, when employing our approach, we can trivially see that the number of
qumodes scales linearly with the number of items.

In conclusion, the PCQO algorithm works considerably well for the UKP problem
cases we considered. The energy profile motivates the investigation of
techniques that find suitable initial parameters to increase the performance
of the algorithm. In QAOA-like algorithms, we have to implement non-Gaussian
operations in the ansatz whereas PCQO can perform the optimization using only
Gaussian operations. The only non-Gaussian element that is introduced is the
number-resolving measurement. As the non-Gaussian operations are relatively
hard to realize experimentally, the PCQO algorithm becomes a preferable
candidate for near-term photonic devices solving integer programming problems.

###  V.2 Maxclique problem

Figure 4: Energy as a function of iteration steps for the $N=5$ and $N=6$
qumode Maxclique problems with a $p=1$ layer ansatz. Solid lines denote the
mean energy of 5 random initializations, while shaded regions represent the
standard error. (a) $N=5$ qumode Maxclique problem, (b) $N=6$ qumode Maxclique
problem. Insets show corresponding graphs and bars for occurrences of the
degenerate solutions. Green lines show the exact energy.

The final benchmark we consider is the Maxclique problem [55]. Consider a
graph $\mathcal{G}=(V,E)$ with $V$ vertices and $E$ edges that can be
represented by its adjacency matrix $Z$, where $Z_{ij}=1$ if there exists an
edge between vertices $i$ and $j$, and $Z_{ij}=0$ otherwise. The Maxclique
problem involves finding the largest subset of vertices $V_{s}\subseteq V$ in
graph $\mathcal{G}$, where all vertices in $V_{s}$ are mutually connected.

To encode the Maxclique problem in terms of Fock states, we introduce binary
variables $n_{i}\in\\{0,1\\}$, representing the selection of a vertex $i$ in
the maximum clique. This leads to

$\displaystyle\min_{n_{i}}\ -\sum_{i\in V}n_{i}\quad\text{subject
to}\quad\sum_{i,j\in V}(\mathbf{I}-Z_{ij})~{}n_{i}n_{j}=0,$ (10)

where $\mathbf{I}$ is an identity matrix. Here, the objective is to maximize
the sum of selected vertices $n_{i}$, indicating the size of the maximum
clique. The constraint ensures that if vertices $i$ and $j$ are not connected
by an edge, they cannot both be selected in the maximum clique. Hence the
problem Hamiltonian in terms of number operators will be

$F(\mathbf{n})=-\sum_{i\in V}n_{i}+\delta_{1}\sum_{i,j\in
V}(\mathbf{I}-Z_{ij})~{}n_{i}n_{j}+\delta_{2}\sum_{i\in V}n_{i}(n_{i}-1)$ (11)

where the first two terms correspond to the Eq. (10) and the third term
ensures that the search of the Fock space is restricted to one. Penalty terms
$(\delta_{1},\delta_{2})$ are applied to change the weights of the constraints
depending upon the problem instance. These problems come under the class of
nonlinear integer programming problems where the cost function and constraints
can both be nonlinear functions of integers.

The PCQO algorithm was applied to find the maximum clique in graphs with $N=5$
and $N=6$ nodes. The mean energy convergence for the $p=1$ layer ansatz is
shown in Fig. 4(a) for $N=5$ and Fig. 4(b) for $N=6$, with the respective
graphs shown as insets. The cutoff dimension was $D=5$ due to computational
limitations. The penalty terms were set to $\delta_{1}=10$ and $\delta_{2}=1$.
Adam optimizer was used as a classical optimizer. The maximum cliques obtained
were $V_{s}=\\{0,1,3\\}$ and $V_{s}=\\{0,2,3\\}$ for both $N=5$ and $N=6$
nodes graph. The exact energy was found in both cases regardless of the
initial parameters chosen. However, for the $N=5$ case, 200 iterations were
needed to reach the exact energy. The inset bar plots show that depending on
the initial parameters, the optimal converges to one of the two degenerate
solutions for both $N=5$ and $N=6$. It is worth noting that since the Fock
states are constrained to be 0 or 1, this algorithm can be implemented with
the same resources in qubit-based technologies as well. Therefore, a
performance comparison between qumode-based and qubit-based approaches would
be interesting.

##  VI Comparison with CV-QAOA

Figure 5: Energy as a function of 200 iteration steps comparing QAOA ($Q=2$,
$Q=10$) with PCQO ($Q=5$) for a $N=3$ UKP case. The best instance out of 5
randomly initialized instances is shown. Different markers depict the energy
convergence for different algorithms and the green solid line shows the exact
energy for the solution.

In this section, we compare the performance of the PCQO algorithm with two
variants of CV-QAOA for a $N=3$ qumode UKP case. In CV-QAOA, we define
$H_{m}=\sum_{i}p_{i}^{2}$ and prepare the initial state as a squeezed state
for all qumodes with squeezing parameter $r=1$. The corresponding mixer term
is given by $U_{b}(\beta)=\exp(-i\beta\sum_{i}p_{i}^{2})$, which can be
implemented using a custom gate $P_{z}$ defined as

$P_{z}(s)=R\left(-\frac{\pi}{2}\right)P(s)R\left(\frac{\pi}{2}\right),$ (12)

where the $R$ gates act as a Fourier transform, rotating the position
quadrature into the momentum quadrature. For the Hamiltonian term, we have
$U_{c}(\gamma)=\exp(-i\gamma F(\mathbf{n}))$, where $F(\mathbf{n})$ is
determined by Eq. (9). Implementing this term involves combining $R$ gates,
$K$ gates, and $CK$ gates (See Table 1). In PCQO, we initialize with a vacuum
state on all qumodes and use the same ansatz as described in the previous
section, consisting of $X$ gates followed by nearest-neighbor $CZ$ gates. In
both algorithms, we perform number-resolving measurements and optimize the
same cost function with cutoff dimension $D=10$.

Regarding parameterization, conventional QAOA has one parameter per unitary,
so for the $p=1$ layer, we have $Q=2$. However, in PCQO, we set one parameter
per gate, resulting in $Q=5$ parameters for the $p=1$ layer because of the
nearest neighbor two-mode gates. To ensure a fair comparison, we also consider
a variant of QAOA called MA-QAOA, where each gate has its free parameter [56].
For the $p=1$ layer, this leads to $Q=10$ parameters due to all-to-all
connected two-mode gates. For classical optimization, we have implemented the
Adam optimizer for all the cases.

In Fig. 5, we show the energy as a function of iteration steps for the $p=1$
layer in three different algorithms: PCQO and QAOA with $Q=2$ and $Q=10$. The
energy values shown correspond to the best outcome out of five randomly
initialized instances. We observe that PCQO outperforms both variants of QAOA,
achieving the exact energy within 200 iterations. This implies that the
operator pool calculated by the nested commutator method serves as a better
ansatz for low-layered algorithms compared to the QAOA ansatz. Additionally,
the performance of QAOA with $Q=10$ surpasses that of QAOA with $Q=2$ due to
the increased degree of freedom given by optimizable parameters.

As QAOA resembles adiabatic evolution, we might require a high-depth circuit
for optimal solutions. It is worth mentioning that the implementation of non-
Gaussian gates is approximate at finite cutoff dimensions, and increasing the
cutoff could potentially enhance the performance of QAOA in our simulations.
Nevertheless, the PCQO algorithm demonstrates superior performance and is
particularly suitable for current near-term devices, as it utilizes native
Gaussian operations, which are easier to implement experimentally. Regarding
phase-space encoding, it becomes apparent that for polynomial functions with
degrees higher than three, decomposing them into lower-order gates would be
necessary for optimization using the available gates. Such decomposition would
require substantial resources that are often unavailable, reinforcing the
preference for PCQO as it allows truncating the operator pool to match the
available gates. For instance, decomposing $e^{-ix^{4}}$ requires 29 quadratic
gates to be decomposed exactly [39].

Figure 6: (a) Schematic diagram for eight-qumode nanophotonic chip. The chip
is divided into a pair of identical qumodes (0,1,2,3) and (4,5,6,7) utilizing
$S_{2}(r)$ gate. This gate can be decomposed as a $S(r)$ gate ($S(r,\phi=0)$
in Table 1) with $r=0$ or $r=1$ and a $BS\left(\frac{\pi}{4},0\right)$ gate as
shown. Then, an arbitrary $U_{4}$ unitary is applied to the pair followed by
number-resolving measurements. (b) shows the PCQO ansatz considered for the
experiment. This includes $R(\phi)$ gates applied to all the qumodes and
$BS(\theta,0)$ gates applied to nearest-neighbor connections.
$\boldsymbol{\theta}=\\{\theta_{1},\theta_{2},\ldots\\}$ are optimizable
parameters. (c) shows the mean photon number for all the qumodes obtained with
the optimal circuit solving $F(\mathbf{n})=(n_{0}+n_{2}-0.75)^{2}$. The
results shown are of a numerical simulation with $D=3$ cutoff and the
experiment with 1000 shots. (d) shows the mean photon number obtained by
taking an average of the identical qumodes for both numerical simulation and
experiment.

##  VII Experimental considerations

One of the notable advantages of PQC is its operability at room temperatures
[57]. In the gate-based approach, photonic qumodes are prepared as Gaussian
states and manipulated using single or two-mode operations, which can be
Gaussian or non-Gaussian [58, 59]. These operations can be implemented using
optical waveguides or integrated photonic systems. The PCQO algorithm,
designed considering current hardware constraints, follows a specific ansatz.
It begins with a vacuum state, applies one and two-mode Gaussian and non-
Gaussian operations, and concludes with a homodyne or number-resolving
measurement. Implementing single-mode operations such as phase shifts and
beamsplitters is straightforward using passive linear optics. Displacement
gates can be achieved by utilizing an ancillary qumode [60]. However,
generating squeezing is challenging as it requires second-order nonlinearity,
which can be accomplished using an optical parametric oscillator [61]. Optical
parametric oscillators can also generate multimode squeezing and entanglement
[62]. Measurement-induced squeezing has also been proposed [63], and
experimental implementations of gates like the $CX$ gate have been achieved
[64]. Additionally, $CZ$ gates can be implemented experimentally using linear
optics and ancillary squeezed vacuum states [60].

Including a non-Gaussian operation becomes crucial when all other components
are Gaussian, as a purely Gaussian state, operations, and measurements can be
classically simulated [51]. However, implementation of non-Gaussian gates
poses experimental challenges in terms of performance [65, 66, 67].
Nevertheless, there are proposals for deterministic cubic-phase gates [68].
Recent proposals have also suggested experimentally feasible continuous-
variable quantum neural networks, where non-linearity is achieved through
repeat-until-success measurements on ancillary qumodes [69].

In terms of measurements, developments have been made in performing homodyne
measurements [70, 71] and number-resolving measurements [72, 73]. In addition
to the circuit ansatz, the classical training part plays a crucial role. To
facilitate a gradient-based optimization routine, the gradients need to be
obtained from the circuit, which can be achieved using the parameter shift
rule [74]. Apart from this, recent experiments have demonstrated in situ
backpropagation for deep learning [75]. Moreover, a quadratic speedup in the
optimization of noisy quantum optical circuits has been investigated [76].
This indicates that the PCQO algorithm can be readily implemented
experimentally using currently available devices and can be further extended
when large-scale quantum computers with continuous-variable technologies
become available [77].

To validate the aforementioned claim, we conducted proof-of-principle
experiments utilizing a state-of-the-art eight qumode fully-programmable
nanophotonic chip [59]. This 10 mm $\times$ 4 mm chip incorporates two-mode
squeezed vacuum states as the initial states, effectively dividing the system
into a pair of four independent squeezed vacuum qumodes. The squeezing
parameter can be chosen as a binary option of $r=1$ or $r=0$. Thereafter, a
programmable SU(4) transformation is applied identically to each of the pair
of qumodes. This transformation is based on a network of six Mach-Zehnder
interferometers arranged in a rectangular configuration [78]. Each of the
interferometers has two tunable parameters which give a total of $Q=12$
parameters. Lastly, number-resolving measurements are simultaneously applied
to all qumodes. A schematic diagram of the circuit is depicted in Fig. 6(a).

For the experiment, we considered a two-mode toy problem in the Hilbert-space
encoding defined as

$F(\mathbf{n})=(n_{0}+n_{2}-0.75)^{2}.$ (13)

For the initial state, the two-mode squeezing parameter was kept to $r=1$ for
$(0,4)$ and $(2,6)$ qumodes, and the other two pairs were kept to zero. This
was done because the $F(\mathbf{n})$ is a function of $n_{0}$ and $n_{2}$, so
the squeezing will help in getting a larger mean-photon number. Similar to the
Section. V, the mixer was chosen to be
$H_{m}=\sum_{i}~{}(x_{i}-x_{0})^{2}+(p_{i}-p_{0})^{2}$. For simplicity, we
considered a four qumode circuit but it is important to note that the same
circuit structure was applied to the other four qumodes as well. To construct
the ansatz for our experiment, we considered the operator pool $\mathcal{A}$
obtained from $A_{\lambda}^{(l)}$ with $l=2$. Among many, this pool included
the $R$ gate and the $BS$ gate, which are native to the hardware. Thus, we
selected these two gates as the building blocks for our circuit ansatz. To
enhance the expressive power of the circuit, we incorporated two ancillary
qumodes in addition to the two qumodes required for the problem. This allowed
us to exploit the full SU(4) transformation. Therefore, the circuit ansatz
consisted of a series of $R$ gates applied to four qumodes, followed by $BS$
gates applied to nearest-neighbor qumodes. A graphical representation of the
ansatz is shown in Fig. 6(b).

Regarding the parameterization, each $R$ gate in the circuit had its own
independent free parameter. The $BS$ gate possesses a transmissivity angle and
a phase angle, as outlined in Table 1. However, in our specific experiments,
we set the phase angle to zero and treated the transmissivity angle as a free
parameter. Consequently, our ansatz was characterized by $Q=7$ free
parameters. Due to limited access to the physical hardware, we ran the
optimization loop as a numerical simulation with a cutoff $D=3$ and
implemented the circuit with the optimized parameters on the chip to obtain
measurement outcomes from 1000 shots. For the optimization, we implemented a
$p=1$ layer circuit with a gradient-free optimizer called COBYLA [79].

Fig. 6(c) illustrates the mean photon numbers obtained from both the simulator
and the experimental setup for all qumodes. Remarkably, even with a moderate
value of $D=3$, we successfully obtained the exact solution to the problem. In
the numerical simulation, the mean photon numbers for qumodes $(0,1,2,3)$
would be exactly the same as qumodes $(4,5,6,7)$ since the operations are
identical. However, in the experimental setting, the distribution deviates due
to inherent limitations such as the finite number of shots, noise, and losses
in the chip. Additionally, the utilization of a low cutoff in the numerical
simulation may have resulted in suboptimal parameters for the actual chip.
Despite these factors, a notable resemblance is observed between the
experimental results and the numerical simulations. To further analyze the
agreement, we computed the average of the mean photon numbers for identical
qumodes $(0,4)$, $(1,5)$. $(2,6)$, and $(3,7)$ and the results are shown in
Fig. 6(d). Encouragingly, this analysis demonstrates a high level of
concurrence between the experimental outcomes and the numerical simulations.
Therefore, the PCQO algorithm provides a promising circuit ansatz that can be
readily implemented using currently available hardware.

It is important to note that the implemented problem possesses relative
simplicity as the chip exclusively incorporates fixed squeezing and lacks
displacement operations. When encountering problems with large integer
solutions, achieving the desired mean photon number becomes challenging
without variable squeezing and displacement. Our current experimental results
focus on demonstrating the feasibility of the PCQO algorithm through simple
experiments. However, future advances enabling variable squeezing or
displacement operations will facilitate tackling more complex problem
instances.

##  VIII Discussions and future work

We proposed a hybrid quantum-classical optimization algorithm for photonic
quantum computing to tackle complex problems with the currently available
technologies. The circuit ansatz for this algorithm is a problem-inspired
ansatz computed by utilizing shortcuts-to-adiabaticity techniques,
specifically counterdiabatic protocols. We investigated the performance of the
algorithm by considering two non-convex continuous-variable optimization
problems up to four variables and with a degree of six. We also considered two
integer programming problems, specifically the unbounded knapsack problem for
up to four system sizes, and the max-clique problem for up to six node graphs.
We observed that the PCQO algorithm successfully finds good approximate
solutions to these problems using a few gates. To showcase the practical
feasibility of PCQO, we conducted experiments on an eight-mode nanophotonic
chip. These experiments substantiated that PCQO can be implemented on near-
term photonic chips, thereby providing a promising avenue for utilizing
photonic quantum computing to solve optimization problems.

We have considered relevant optimization problems for the industry, but this
algorithm can be extended to study physical problems as well. PCQO is a hybrid
algorithm but purely quantum counterdiabatic algorithms can be developed in
the future to study the performance from a point of view of shortcuts-to-
adiabaticity. The backend for this algorithm is a photonic system but this can
be extended to any bosonic systems as well and the performance analysis would
be interesting in this regard. Advanced machine learning techniques like
reinforcement learning [80] and adaptive techniques [81] are other aspects
that can be incorporated in PCQO to select the circuit ansatz in a better way.
Also, finding effective initialization strategies would be interesting for
future work. We believe that this work will serve as a benchmark for designing
more advanced hybrid qubit-bosonic optimization algorithms [82]. In summary,
this work introduces the PCQO algorithm as a compelling approach for
addressing hard optimization problems using photonic quantum computing. The
successful application of PCQO to various problem domains, combined with its
potential for further advancements and extensions, positions photonic quantum
computing as a competitive candidate alongside qubit-based technologies for
tackling challenging optimization tasks.

###### Acknowledgements.

We acknowledge the use of Strawberryfields Library [83] for performing the
simulations and the experiment. The authors acknowledge Tasio Gonzalez-Raya,
Narendra Hegade, and Martin Larocca for useful discussions. This work is
supported by EU FET Open Grant EPIQUS (899368), and the Basque Government
through Grant No. IT1470-22, the project grant PID2021-126273NB-I00 funded by
MCIN/AEI/10.13039/501100011033 and by “ERDF A way of making Europe” and “ERDF
Invest in your Future”, the Spanish CDTI through Plan complementario
Comunicación cuántica (EXP. 2022/01341)(A/20220551), and project
OpenSuperQ+100 (101113946) of the EU Flagship on Quantum Technologies, and the
IKUR Strategy under the collaboration agreement between Ikerbasque Foundation
and BCAM on behalf of the Department of Education of the Basque Government.
M.S. acknowledges support from Spanish Ramón y Cajal Grant RYC-2020-030503-I.
X.C. acknowledges ayudas para contratos Ramón y Cajal–2015-2020
(RYC-2017-22482). MGdA acknowledges support from the UPV/EHU and TECNALIA 2021
PIF contract call.

## References

  * Preskill [2018] J. Preskill, Quantum Computing in the NISQ era and beyond, Quantum 2, 79 (2018).
  * Cerezo  _et al._ [2021] M. Cerezo, A. Arrasmith, R. Babbush, S. C. Benjamin, S. Endo, K. Fujii, J. R. McClean, K. Mitarai, X. Yuan, L. Cincio, and P. J. Coles, Variational quantum algorithms, Nature Reviews Physics 3, 625 (2021).
  * Cade  _et al._ [2020] C. Cade, L. Mineh, A. Montanaro, and S. Stanisic, Strategies for solving the fermi-hubbard model on near-term quantum computers, Phys. Rev. B 102, 235122 (2020).
  * Bravo-Prieto  _et al._ [2020] C. Bravo-Prieto, J. Lumbreras-Zarapico, L. Tagliacozzo, and J. I. Latorre, Scaling of variational quantum circuit depth for condensed matter systems, Quantum 4, 272 (2020).
  * Colless  _et al._ [2018] J. I. Colless, V. V. Ramasesh, D. Dahlen, M. S. Blok, M. E. Kimchi-Schwartz, J. R. McClean, J. Carter, W. A. de Jong, and I. Siddiqi, Computation of molecular spectra on a quantum processor with an error-resilient algorithm, Phys. Rev. X 8, 011021 (2018).
  * Grimsley  _et al._ [2019a] H. R. Grimsley, S. E. Economou, E. Barnes, and N. J. Mayhall, An adaptive variational algorithm for exact molecular simulations on a quantum computer, Nature Communications 10, 3007 (2019a).
  * Kandala  _et al._ [2017] A. Kandala, A. Mezzacapo, K. Temme, M. Takita, M. Brink, J. M. Chow, and J. M. Gambetta, Hardware-efficient variational quantum eigensolver for small molecules and quantum magnets, Nature 549, 242 (2017).
  * Anschuetz  _et al._ [2019] E. Anschuetz, J. Olson, A. Aspuru-Guzik, and Y. Cao, Variational quantum factoring, Quantum Technology and Optimization Problems , 74 (2019).
  * Karamlou  _et al._ [2021] A. H. Karamlou, W. A. Simon, A. Katabarwa, T. L. Scholten, B. Peropadre, and Y. Cao, Analyzing the performance of variational quantum factoring on a superconducting quantum processor, npj Quantum Information 7, 156 (2021).
  * Robert  _et al._ [2021] A. Robert, P. K. Barkoutsos, S. Woerner, and I. Tavernelli, Resource-efficient quantum algorithm for protein folding, npj Quantum Information 7, 38 (2021).
  * de Andoin  _et al._ [2022] M. G. de Andoin, E. Osaba, I. Oregi, E. Villar-Rodriguez, and M. Sanz, Hybrid quantum-classical heuristic for the bin packing problem, Proceedings of the Genetic and Evolutionary Computation Conference Companion GECCO ’22, 2214–2222 (2022).
  * Killoran  _et al._ [2019a] N. Killoran, T. R. Bromley, J. M. Arrazola, M. Schuld, N. Quesada, and S. Lloyd, Continuous-variable quantum neural networks, Phys. Rev. Res. 1, 033063 (2019a).
  * Verdon  _et al._ [2019] G. Verdon, J. M. Arrazola, K. Brádler, and N. Killoran, A quantum approximate optimization algorithm for continuous problems, arXiv preprint arXiv:1902.00409 (2019).
  * Khosravi  _et al._ [2021] F. Khosravi, A. Scherer, and P. Ronagh, Mixed-integer programming using a bosonic quantum computer, arXiv preprint arXiv:2112.13917 (2021).
  * Mezher  _et al._ [2023] R. Mezher, A. F. Carvalho, and S. Mansfield, Solving graph problems with single-photons and linear optics, arXiv preprint arXiv:2301.09594 (2023).
  * Yeter-Aydeniz  _et al._ [2022] K. Yeter-Aydeniz, E. Moschandreou, and G. Siopsis, Quantum imaginary-time evolution algorithm for quantum field theories with continuous variables, Phys. Rev. A 105, 012412 (2022).
  * Pati  _et al._ [2000] A. K. Pati, S. L. Braunstein, and S. Lloyd, Quantum searching with continuous variables, arXiv preprint quant-ph/0002082 (2000).
  * Douce  _et al._ [2017] T. Douce, D. Markham, E. Kashefi, E. Diamanti, T. Coudreau, P. Milman, P. van Loock, and G. Ferrini, Continuous-variable instantaneous quantum computing is hard to sample, Phys. Rev. Lett. 118, 070503 (2017).
  * Gottesman  _et al._ [2001] D. Gottesman, A. Kitaev, and J. Preskill, Encoding a qubit in an oscillator, Phys. Rev. A 64, 012310 (2001).
  * Terhal  _et al._ [2020] B. M. Terhal, J. Conrad, and C. Vuillot, Towards scalable bosonic quantum error correction, Quantum Science and Technology 5, 043001 (2020).
  * Arrazola  _et al._ [2019] J. M. Arrazola, T. R. Bromley, J. Izaac, C. R. Myers, K. Brádler, and N. Killoran, Machine learning method for state preparation and gate synthesis on photonic quantum computers, Quantum Science and Technology 4, 024004 (2019).
  * Enomoto  _et al._ [2022] Y. Enomoto, K. Anai, K. Udagawa, and S. Takeda, Continuous-variable quantum approximate optimization on a programmable photonic quantum processor, arXiv preprint arXiv:2206.07214 (2022).
  * del Campo [2013] A. del Campo, Shortcuts to adiabaticity by counterdiabatic driving, Phys. Rev. Lett. 111, 100502 (2013).
  * Chen  _et al._ [2011a] X. Chen, E. Torrontegui, and J. G. Muga, Lewis-riesenfeld invariants and transitionless quantum driving, Phys. Rev. A 83, 062116 (2011a).
  * Torrontegui  _et al._ [2013] E. Torrontegui, S. Ibáñez, S. Martínez-Garaot, M. Modugno, A. del Campo, D. Guéry-Odelin, A. Ruschhaupt, X. Chen, and J. G. Muga, Chapter 2 - shortcuts to adiabaticity, Advances in Atomic, Molecular, and Optical Physics Advances In Atomic, Molecular, and Optical Physics, 62, 117 (2013).
  * Chen  _et al._ [2011b] X. Chen, E. Torrontegui, and J. G. Muga, Lewis-riesenfeld invariants and transitionless quantum driving, Phys. Rev. A 83, 062116 (2011b).
  * Chandarana  _et al._ [2022a] P. Chandarana, N. N. Hegade, K. Paul, F. Albarrán-Arriagada, E. Solano, A. del Campo, and X. Chen, Digitized-counterdiabatic quantum approximate optimization algorithm, Phys. Rev. Res. 4, 013141 (2022a).
  * Hegade  _et al._ [2021a] N. N. Hegade, K. Paul, Y. Ding, M. Sanz, F. Albarrán-Arriagada, E. Solano, and X. Chen, Shortcuts to adiabaticity in digitized adiabatic quantum computing, Phys. Rev. Applied 15, 024038 (2021a).
  * Hegade  _et al._ [2021b] N. N. Hegade, K. Paul, F. Albarrán-Arriagada, X. Chen, and E. Solano, Digitized adiabatic quantum factorization, Phys. Rev. A 104, L050403 (2021b).
  * Hegade  _et al._ [2022a] N. N. Hegade, P. Chandarana, K. Paul, X. Chen, F. Albarrán-Arriagada, and E. Solano, Portfolio optimization with digitized counterdiabatic quantum algorithms, Phys. Rev. Res. 4, 043204 (2022a).
  * Hegade  _et al._ [2022b] N. N. Hegade, X. Chen, and E. Solano, Digitized counterdiabatic quantum optimization, Phys. Rev. Res. 4, L042030 (2022b).
  * Rosenbrock [1960] H. H. Rosenbrock, An Automatic Method for Finding the Greatest or Least Value of a Function, The Computer Journal 3, 175 (1960).
  * Styblinski and Tang [1990] M. Styblinski and T.-S. Tang, Experiments in nonconvex optimization: Stochastic approximation with function smoothing and simulated annealing, Neural Networks 3, 467 (1990).
  * Kolman and Beck [1995] B. Kolman and R. E. Beck, 4 - integer programming, Elementary Linear Programming with Applications (Second Edition) , 249 (1995).
  * Farhi  _et al._ [2014] E. Farhi, J. Goldstone, and S. Gutmann, A quantum approximate optimization algorithm, arXiv preprint arXiv:1411.4028 (2014).
  * Blekos  _et al._ [2023] K. Blekos, D. Brand, A. Ceschini, C.-H. Chou, R.-H. Li, K. Pandya, and A. Summer, A review on quantum approximate optimization algorithm and its variants, arXiv preprint arXiv:2306.09198 (2023).
  * Barkoutsos  _et al._ [2020] P. K. Barkoutsos, G. Nannicini, A. Robert, I. Tavernelli, and S. Woerner, Improving Variational Quantum Optimization using CVaR, Quantum 4, 256 (2020).
  * Lloyd and Braunstein [1999] S. Lloyd and S. L. Braunstein, Quantum computation over continuous variables, Phys. Rev. Lett. 82, 1784 (1999).
  * Kalajdzievski and Arrazola [2019] T. Kalajdzievski and J. M. Arrazola, Exact gate decompositions for photonic quantum computing, Phys. Rev. A 99, 022341 (2019).
  * Albash and Lidar [2018] T. Albash and D. A. Lidar, Adiabatic quantum computation, Rev. Mod. Phys. 90, 015002 (2018).
  * Barends  _et al._ [2016] R. Barends, A. Shabani, L. Lamata, J. Kelly, A. Mezzacapo, U. L. Heras, R. Babbush, A. G. Fowler, B. Campbell, Y. Chen, Z. Chen, B. Chiaro, A. Dunsworth, E. Jeffrey, E. Lucero, A. Megrant, J. Y. Mutus, M. Neeley, C. Neill, P. J. J. O’Malley, C. Quintana, P. Roushan, D. Sank, A. Vainsencher, J. Wenner, T. C. White, E. Solano, H. Neven, and J. M. Martinis, Digitized adiabatic quantum computing with a superconducting circuit, Nature 534, 222 (2016).
  * Hadfield  _et al._ [2019] S. Hadfield, Z. Wang, B. O’Gorman, E. G. Rieffel, D. Venturelli, and R. Biswas, From the quantum approximate optimization algorithm to a quantum alternating operator ansatz, Algorithms 12 (2019).
  * Demirplak and Rice [2003] M. Demirplak and S. A. Rice, Adiabatic population transfer with control fields, The Journal of Physical Chemistry A 107, 9937 (2003).
  * Berry [2009] M. V. Berry, Transitionless quantum driving, Journal of Physics A: Mathematical and Theoretical 42, 365303 (2009).
  * Sels and Polkovnikov [2017] D. Sels and A. Polkovnikov, Minimizing irreversible losses in quantum systems by local counterdiabatic driving, Proceedings of the National Academy of Sciences 114, E3909 (2017).
  * Claeys  _et al._ [2019] P. W. Claeys, M. Pandey, D. Sels, and A. Polkovnikov, Floquet-engineering counterdiabatic protocols in quantum many-body systems, Phys. Rev. Lett. 123, 090602 (2019).
  * Sun  _et al._ [2022] D. Sun, P. Chandarana, Z.-H. Xin, and X. Chen, Optimizing counterdiabaticity by variational quantum circuits, Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 380, 20210282 (2022).
  * Chandarana  _et al._ [2022b] P. Chandarana, N. N. Hegade, I. Montalban, E. Solano, and X. Chen, Digitized-counterdiabatic quantum algorithm for protein folding, arXiv preprint arXiv:2212.13511 (2022b).
  * Farhi  _et al._ [2000] E. Farhi, J. Goldstone, S. Gutmann, and M. Sipser, Quantum computation by adiabatic evolution, arXiv preprint quant-ph/0001106 (2000).
  * Maamache and Saadi [2008] M. Maamache and Y. Saadi, Adiabatic theorem and generalized geometrical phase in the case of continuous spectra, Phys. Rev. Lett. 101, 150407 (2008).
  * Bartlett  _et al._ [2002] S. D. Bartlett, B. C. Sanders, S. L. Braunstein, and K. Nemoto, Efficient classical simulation of continuous variable quantum information processes, Phys. Rev. Lett. 88, 097904 (2002).
  * Kingma and Ba [2014] D. P. Kingma and J. Ba, Adam: A method for stochastic optimization, arXiv preprint arXiv:1412.6980 (2014).
  * Stein  _et al._ [2023] J. Stein, F. Chamanian, M. Zorn, J. Nüßlein, S. Zielinski, M. Kölle, and C. Linnhoff-Popien, Evidence that pubo outperforms qubo when solving continuous optimization problems with the qaoa, arXiv preprint arXiv:2305.03390 (2023).
  * Lueker [1975] G. Lueker, Two np-complete problems in nonnegative integer programming, Princeton University. Department of Electrical Engineering (1975).
  * Bomze  _et al._ [1999] I. M. Bomze, M. Budinich, P. M. Pardalos, and M. Pelillo, The maximum clique problem, Handbook of Combinatorial Optimization: Supplement Volume A , 1 (1999).
  * Herrman  _et al._ [2022] R. Herrman, P. C. Lotshaw, J. Ostrowski, T. S. Humble, and G. Siopsis, Multi-angle quantum approximate optimization algorithm, Scientific Reports 12, 6781 (2022).
  * Rudolph [2017] T. Rudolph, Why I am optimistic about the silicon-photonic route to quantum computing, APL Photonics 2 (2017).
  * O’Brien  _et al._ [2009] J. L. O’Brien, A. Furusawa, and J. Vučković, Photonic quantum technologies, Nature Photonics 3, 687 (2009).
  * Arrazola  _et al._ [2021] J. M. Arrazola, V. Bergholm, K. Brádler, T. R. Bromley, M. J. Collins, I. Dhand, A. Fumagalli, T. Gerrits, A. Goussev, L. G. Helt, J. Hundal, T. Isacsson, R. B. Israel, J. Izaac, S. Jahangiri, R. Janik, N. Killoran, S. P. Kumar, J. Lavoie, A. E. Lita, D. H. Mahler, M. Menotti, B. Morrison, S. W. Nam, L. Neuhaus, H. Y. Qi, N. Quesada, A. Repingon, K. K. Sabapathy, M. Schuld, D. Su, J. Swinarton, A. Száva, K. Tan, P. Tan, V. D. Vaidya, Z. Vernon, Z. Zabaneh, and Y. Zhang, Quantum circuits with many photons on a programmable nanophotonic chip, Nature 591, 54 (2021).
  * Filip  _et al._ [2005] R. Filip, P. Marek, and U. L. Andersen, Measurement-induced continuous-variable quantum interactions, Phys. Rev. A 71, 042308 (2005).
  * Wu  _et al._ [1986] L.-A. Wu, H. J. Kimble, J. L. Hall, and H. Wu, Generation of squeezed states by parametric down conversion, Phys. Rev. Lett. 57, 2520 (1986).
  * Ou  _et al._ [1992] Z. Y. Ou, S. F. Pereira, H. J. Kimble, and K. C. Peng, Realization of the einstein-podolsky-rosen paradox for continuous variables, Phys. Rev. Lett. 68, 3663 (1992).
  * Yoshikawa  _et al._ [2007] J.-i. Yoshikawa, T. Hayashi, T. Akiyama, N. Takei, A. Huck, U. L. Andersen, and A. Furusawa, Demonstration of deterministic and high fidelity squeezing of quantum information, Phys. Rev. A 76, 060301 (2007).
  * Yoshikawa  _et al._ [2008] J.-i. Yoshikawa, Y. Miwa, A. Huck, U. L. Andersen, P. van Loock, and A. Furusawa, Demonstration of a quantum nondemolition sum gate, Phys. Rev. Lett. 101, 250501 (2008).
  * Shapiro [2006] J. H. Shapiro, Single-photon kerr nonlinearities do not help quantum computation, Phys. Rev. A 73, 062305 (2006).
  * Gea-Banacloche [2010] J. Gea-Banacloche, Impossibility of large phase shifts via the giant kerr effect with single-photon wave packets, Phys. Rev. A 81, 043823 (2010).
  * He and Scherer [2012] B. He and A. Scherer, Continuous-mode effects and photon-photon phase gate performance, Phys. Rev. A 85, 033814 (2012).
  * Yanagimoto  _et al._ [2020] R. Yanagimoto, T. Onodera, E. Ng, L. G. Wright, P. L. McMahon, and H. Mabuchi, Engineering a kerr-based deterministic cubic phase gate via gaussian operations, Phys. Rev. Lett. 124, 240503 (2020).
  * Bangar  _et al._ [2022] S. Bangar, G. Siopsis, and K. Yeter-Aydeniz, Experimentally realizable continuous-variable quantum neural networks, Quantum 2.0 Conference and Exhibition , QTu2A.4 (2022).
  * Yuen and Chan [1983] H. P. Yuen and V. W. Chan, Noise in homodyne and heterodyne detection., Opt Lett 8, 177 (1983).
  * Raffaelli  _et al._ [2018] F. Raffaelli, G. Ferranti, D. H. Mahler, P. Sibson, J. E. Kennard, A. Santamato, G. Sinclair, D. Bonneau, M. G. Thompson, and J. C. F. Matthews, A homodyne detector integrated onto a photonic chip for measuring quantum states and generating random numbers, Quantum Science and Technology 3, 025003 (2018).
  * Calkins  _et al._ [2013] B. Calkins, P. L. Mennea, A. E. Lita, B. J. Metcalf, W. S. Kolthammer, A. Lamas-Linares, J. B. Spring, P. C. Humphreys, R. P. Mirin, J. C. Gates, P. G. R. Smith, I. A. Walmsley, T. Gerrits, and S. W. Nam, High quantum-efficiency photon-number-resolving detector for photonic on-chip information processing, Opt. Express 21, 22657 (2013).
  * Höpker  _et al._ [2019] J. P. Höpker, T. Gerrits, A. Lita, S. Krapick, H. Herrmann, R. Ricken, V. Quiring, R. Mirin, S. W. Nam, C. Silberhorn, and T. J. Bartley, Integrated transition edge sensors on titanium in-diffused lithium niobate waveguides, APL Photonics 4 (2019).
  * Schuld  _et al._ [2019] M. Schuld, V. Bergholm, C. Gogolin, J. Izaac, and N. Killoran, Evaluating analytic gradients on quantum hardware, Phys. Rev. A 99, 032331 (2019).
  * Pai  _et al._ [2023] S. Pai, Z. Sun, T. W. Hughes, T. Park, B. Bartlett, I. A. D. Williamson, M. Minkov, M. Milanizadeh, N. Abebe, F. Morichetti, A. Melloni, S. Fan, O. Solgaard, and D. A. B. Miller, Experimentally realized in situ backpropagation for deep learning in photonic neural networks, Science 380, 398 (2023).
  * De Prins  _et al._ [2023] R. De Prins, Y. Yao, A. Apte, and F. M. Miatto, A quadratic speedup in the optimization of noisy quantum optical circuits, arXiv preprint arXiv:2303.08879 (2023).
  * Fukui and Takeda [2022] K. Fukui and S. Takeda, Building a large-scale quantum computer with continuous-variable optical technologies, Journal of Physics B: Atomic, Molecular and Optical Physics 55, 012001 (2022).
  * Clements  _et al._ [2016] W. R. Clements, P. C. Humphreys, B. J. Metcalf, W. S. Kolthammer, and I. A. Walmsley, Optimal design for universal multiport interferometers, Optica 3, 1460 (2016).
  * Powell [1994] M. J. D. Powell, A direct search optimization method that models the objective and constraint functions by linear interpolation, Advances in Optimization and Numerical Analysis , 51 (1994).
  * Yao  _et al._ [2021] J. Yao, L. Lin, and M. Bukov, Reinforcement learning for many-body ground-state preparation inspired by counterdiabatic driving, Phys. Rev. X 11, 031070 (2021).
  * Grimsley  _et al._ [2019b] H. R. Grimsley, S. E. Economou, E. Barnes, and N. J. Mayhall, An adaptive variational algorithm for exact molecular simulations on a quantum computer, Nature Communications 10, 3007 (2019b).
  * Stavenger  _et al._ [2022] T. J. Stavenger, E. Crane, K. C. Smith, C. T. Kang, S. M. Girvin, and N. Wiebe, C2qa - bosonic qiskit, 2022 IEEE High Performance Extreme Computing Conference (HPEC) , 1 (2022).
  * Killoran  _et al._ [2019b] N. Killoran, J. Izaac, N. Quesada, V. Bergholm, M. Amy, and C. Weedbrook, Strawberry Fields: A software platform for photonic quantum computing, Quantum 3, 129 (2019b).