# Solving the max-3-cut problem with coherent networks

S. L. Harrison<EMAIL_ADDRESS>School of Physics and Astronomy,
University of Southampton, Southampton, SO17 1BJ, United Kingdom  H.
Sigurdsson  School of Physics and Astronomy, University of Southampton,
Southampton, SO17 1BJ, United Kingdom  Skolkovo Institute of Science and
Technology, Bolshoy Boulevard 30, bld. 1, Moscow, 121205, Russia  S. Alyatkin
Skolkovo Institute of Science and Technology, Bolshoy Boulevard 30, bld. 1,
Moscow, 121205, Russia  J. D. Töpfer  School of Physics and Astronomy,
University of Southampton, Southampton, SO17 1BJ, United Kingdom  Skolkovo
Institute of Science and Technology, Bolshoy Boulevard 30, bld. 1, Moscow,
121205, Russia  P. G. Lagoudakis<EMAIL_ADDRESS>Skolkovo Institute
of Science and Technology, Bolshoy Boulevard 30, bld. 1, Moscow, 121205,
Russia  School of Physics and Astronomy, University of Southampton,
Southampton, SO17 1BJ, United Kingdom

###### Abstract

Many computational problems are intractable through classical computing and,
as Moore’s law is drawing to a halt, demand for finding alternative methods in
tackling these problems is growing. Here, we realize a liquid light machine
for the NP-hard max-3-cut problem based on a network of synchronized exciton-
polariton condensates. We overcome the binary limitation of the decision
variables in Ising machines using the continuous-phase degrees of freedom of a
coherent network of polariton condensates. The condensate network dynamical
transients provide optically-fast annealing of the XY Hamiltonian. We apply
the Goemans and Williamson random hyperplane technique, discretizing the XY
ground state spin configuration to serve as ternary decision variables for an
approximate optimal solution to the max-3-cut problem. Applications of the
presented coherent network are investigated in image-segmentation tasks and in
circuit design.

##  I Introduction

Complexity in nature is as widespread as it is diverse, and with the turn of
the age there has been a rapid rise in non-Boolean strategies designed to
tackle complex computational problems too cumbersome for conventional Turing-
based computers which are limited by the Von-Neumann Bottleneck, as well as
CMOS architectures approaching their limits. The need for new computational
methods is underscored by many scientific fields such as those devoted to
climate change [1, 2], drug design [3], development of new materials and
batteries [4], and so on, that are dependent on solving complex problems. A
class of such problems is the computationally intractable NP-complete class,
where an estimated solution can be verified in polynomial time, though no
algorithm exists to calculate an exact solution. As these problems are so
widely encountered, they are often approached using approximation algorithms
or heuristic methods such as semi-definite programming [5, 6], genetic
algorithms [7], and nature inspired heuristic algorithms [8, 9]. However,
instead of building an approximate optimization algorithm concerned with
minimizing a cost function, an alternative is mapping the problem to a
physical system that relaxes to the ground state of its energy landscape
corresponding to the global minimum of the cost function. Coherent networks,
which are being regarded as the next possible generation of both quantum [10]
and classical computational devices [11], have generated much interest with
photonic-based classical annealers already realized for both Ising [12, 13,
14, 15, 16, 17, 18, 19, 20] and XY spin Hamiltonians [21, 22, 23, 24, 25],
waveguide networks for the subset sum problem [26], and digital degenerate
cavity laser for the phase retrieval problem [27].

Here, we test the concept of a liquid light machine based on planar networks
of exciton-polariton condensates in semiconductor microcavities and
demonstrate their ability to optimize a maximum-3-cut computational problem
compared to the brute force method with focus on two types of real world
applications, image-segmentation and constrained-via-minimization in circuit
design. Utilizing the recent developments connecting the dissipative nature of
polariton condensate dynamics to minimization of the XY model [22, 28, 29] we
apply a random hyperplane technique to bin the XY ground state spins obtained
from our polariton network into ternary decision variables [5, 30]. We proceed
to map the NP-hard max-3-cut (M3C) optimization problem [31] to the energy
minimization of a ternary phase-discretized XY model, which is then near-
optimally solved using the decision variables obtained from the standard XY
model. Variants on the max-cut problem have many applications, including
social network modelling [32], statistical physics [33], portfolio risk
analysis [34], circuit design [35], image-segmentation [36], and more [37, 38,
39].

Unlike Ising machines and quantum annealers, which can only be mapped to the
max-2-cut problem [33, 14, 40], the continuous degree of freedom of the
variables (spins) in XY systems makes their approximate partition into higher
dimensional decision variables possible, with a direct mapping to the M3C
problem. It also offers a perspective on whether higher order max-$k$-cuts can
be approximated using such continuous phase systems. Our method can be applied
to any system of interacting oscillators defined by their relative phases such
as electronic circuits, laser systems, and condensates. However, polariton
condensate systems have potential advantages over other light-based optimizers
[41, 42, 25, 43, 44, 45] due to the continuous-phase-locking capabilities of
the condensates, strong nonlinearities, the ability to arbitrarily control
polariton graph geometries and coupling strengths [46], their long coherence
length (reported up to 120 $\upmu$m [47]), their robust network
synchronization resulting from the ballistic nature of the condensate
interaction mechanism [48, 49], plus the robustness of the inorganic
semiconductor microcavity [50].

##  II Mapping the M3C Problem to a Ternary XY Model

The M3C problem is as follows: given an undirected graph
$\mathcal{G}=(\mathcal{V},\mathcal{E})$ consisting of vertices $\mathcal{V}$
and weighted edges $\mathcal{E}$, the M3C is the partition of $\mathcal{V}$
into three subsets, such that the sum of all edge weights that connect between
different subsets (labelled the “weight of cut”) is maximized. Following
similar arguments presented by Barahona for an Ising spin system [35] it can
be shown that a map from maximization of a M3C problem to minimizing the
energy of a ternary XY spin system,

$H_{T}=-\sum_{ij}J_{ij}\mathbf{s}_{i}\cdot\mathbf{s}_{j},$ (1)

where

$\displaystyle\mathbf{s}$ $\displaystyle=\begin{pmatrix}\cos(\theta)\\\
\sin(\theta)\\\
\end{pmatrix},\qquad\theta\in\left\\{0,\frac{2\pi}{3},\frac{4\pi}{3}\right\\}.$
(2)

exists.

Let us define an undirected graph $\mathcal{G}=(\mathcal{V},\mathcal{E})$ with
edges $\mathcal{E}$ connecting vertices $\mathcal{V}$, where $\mathcal{V}$
contain the ternary spins $\mathbf{s}$ from Eq. (2). Each edge connecting
vertex $\mathcal{V}_{i}$ and $\mathcal{V}_{j}$, with corresponding spins
$\mathbf{s}_{i}$ and $\mathbf{s}_{j}$, is assigned a weight $J_{ij}=J_{ji}$.
We define three sets of spins corresponding to the three orientations of
$\theta_{i}$ given by Eq. (2),

$\mathcal{V}_{n}=\left\\{i\in\mathcal{V}\ |\ \theta_{i}=\frac{2\pi
n}{3}\right\\},\qquad n=0,1,2.$ (3)

Let us define $\mathcal{E}_{n}$ as the set of edges connecting spins within
each set $\mathcal{V}_{n}$ and $\delta\mathcal{E}$ as the set of edges
connecting spins between different sets of vertices. We can then rewrite Eq.
(1) in the following manner,

$H_{T}=-\sum_{i,j\in\mathcal{E}_{n}}J_{ij}+\frac{1}{2}\sum_{i,j\in\delta\mathcal{E}}J_{ij}.$
(4)

By defining the sum total of all the edges $C=\sum_{ij\in\mathcal{E}}J_{ij}$,
we then have,

$H_{T}+C=\frac{3}{2}\sum_{i,j\in\delta\mathcal{E}}J_{ij}.$ (5)

Since $C$ is invariant on the spin configuration, it can be seen that
minimizing Eq. (1) is the same as maximizing a M3C problem by redefining
$J_{ij}=-\mathcal{W}_{ij}$, where $\mathcal{W}_{ij}$ is the weight of the edge
connecting $\mathcal{V}_{i}$ and $\mathcal{V}_{j}$,

$\text{min}{\left[H_{T}\right]}\quad\leftrightarrow\quad\text{max}{\left[\sum_{i,j\in\delta\mathcal{E}}\mathcal{W}_{ij}\right]}.$
(6)

For comparison, the spins in an Ising system would be binary in orientation
(i.e., $\theta\in\\{0,\pi\\}$). For any two ternary spins we have
$\mathbf{s}_{i}\cdot\mathbf{s}_{j}=\cos(\theta_{i}-\theta_{j})\in$
$\left\\{1,-\frac{1}{2}\right\\}$. This model is also known as the $q=3$
vector Potts model which has been studied in statistical mechanics and quite
recently in the context of exciton-polariton condensates [51]. When
$\theta_{i}\in[-\pi,\pi)$ then Eq. (1) is just the standard XY model,

$H_{XY}=-\sum_{ij}J_{ij}\cos{(\theta_{i}-\theta_{j})},\qquad\theta_{i}\in[-\pi,\pi).$
(7)

In the polariton network, the interacting quantities in question are nonlinear
oscillators (the condensates) each characterized by a complex-valued number
$\psi_{i}=\rho_{i}e^{i\theta_{i}}$. In this sense, the phasors of the
condensates then take the role of interacting two-dimensional pseudospins
$\mathbf{s}$ which dynamically experience a gradient descent towards the
ground state of $H_{XY}$ [22].

As pointed out by Frieze and Jerrum [52], a mapping exists between a
max-$k$-cut problem and the ground state of a system of spins in the vertices
of an equilateral simplex in $\mathbb{R}^{k-1}$. For the max-2-cut problem the
corresponding simplex is a line with vertices $\theta=\\{0,\pi\\}$ which is
the reason an Ising system ground state maps directly to the max-2-cut [33].
For the M3C problem the simplex is an equilateral triangle inscribed by the
unit circle like in Eq. (1) and shown in Fig. 1(a). For the max-4-cut the
problem maps to the ground state of the Heisenberg spin system where the
simplex is an equilateral tetrahedron inscribed by the unit sphere, and so on.

###  II.1 The random hyperplane technique

To approach the ground state of the ternary Hamiltonian [Eq. (1)] the spins of
the XY Hamiltonian [Eq. (7)] ground state configuration
$\\{\mathbf{s}_{1},\mathbf{s}_{2},\dots\\}$, which can be approximately
obtained using a polariton network [22, 28], are projected (or binned) onto
their closest ternary counterpart corresponding to the vertices of the
inscribed triangle [Fig. 1(a)]. The method of projecting continuous spin
variables onto binary and ternary decision variables forms a semi-definite
relaxation program which was studied in computer science by Goemans and
Williamson [5, 30] and referred to as the random hyperplane technique. Here,
we focus on the ternary problem where a binning boundary is depicted with
dashed white lines in Fig. 1(a) like used in [30] which projects the
continuous spins of our polariton network (see example black arrows) into
their closest corner of the circumscribed triangle. An approximate solution of
the ternary spin Hamiltonian is obtained through sampling many random
orientations of the binning boundary with respect to the horizontal axis [30]
(hence the name “random” hyperplane technique). We point out that the
minimization of the ternary [Eq. (1)] and continuous [Eq. (7)] XY system share
the same relaxation procedure to a semi-definite program [5, 6] underlining
the common point of finding the ground states of the two systems.

Figure 1: (a) Schematic of a chosen binning boundary (dashed white lines)
projecting spins into the corners of the circumscribed triangle. (b) The XY
ground state of the house configuration with numbers representing the angles
(phases) $\theta_{i}$ of the vertices in radians; (c) the ternary mapping of
the phases in (b); (d, e) histograms of the M3C weights and XY energy
respectively from the corresponding house graph of (blue) 1000 random samples
of partitioning, (yellow) 2DGPE simulations and (purple) experiment. Dashed
black lines show the XY ground state on spectra, not discretized to the
histogram bins (white and gray background stripes). Weight of cuts in (d)
represent the optimum cut for each data set sampled across 100 uniformally
spaced binning boundaries (different rotations of the unit triangle) and the
average error, $\bar{S}_{W}$, of the simulated and experimental cuts are shown
in yellow and purple respectively. (f, i) Optical pump profile
$P(\mathbf{r})$; (g, j) condensate density $|\Psi(\mathbf{r})|^{2}$, and (h,
k) condensate phase map $\theta(\mathbf{r})$ from 2DGPE simulation and
experiment, respective to each row, of optically trapped polariton condensates
[29] in the AFM house configuration, with arrows in (h, k) showing the
discretized XY phase into subsets $\mathcal{V}_{0,1,2}$.

###  II.2 The house graph

An example procedure is shown visually in Fig. 1(b) and 1(c) using the XY
ground state of the anti-ferromagnetic (AFM) house graph, as previously
studied in [25] using degenerate laser cavities. Here, AFM refers to
$J_{ij}<0$ in Eqs. (1) and (7) with preferential antiparallel spin alignement
$\theta_{i}-\theta_{j}=\pi$ between two condensates/oscillators. The contrary,
ferromagnetic (FM) alignment refers to preferential in-phase
condensate/oscillator synchronization $\theta_{i}-\theta_{j}=0$, corresponding
to $J_{ij}>0$ interactions. The AFM house graph consists of vertices [colored
discs in Fig. 1(b)] arranged in the illustrated fashion with equally weighted
AFM edges $J_{ij}=J<0$. For brevity, we will use $J=-1$ in dimensionless
units. We first calculate the XY ground state using the conventional basin
hopping optimization method [53] with the ground state angles (phases)
$\theta_{i}$ given in radians inside the discs. An appropriately oriented
binning boundary bins the angles into their ternary counterparts shown inside
the discs in Fig. 1(c). We then apply this ternary XY outcome to the M3C
problem through Eq. (6) which dictates that the partitions (and their cuts)
follow the colored regions shown in Fig. 1(c). This is precisely the maximum
cut of the M3C problem for the house graph, with weight $=6$. It is worth
noting that the minimum weight of a cut here is $=2$.

To illustrate that the maximum weight has been obtained, we plot in Fig. 1(d)
the M3C weight (dashed black line) against the distribution of possible house
graph weights obtained from 1000 random samples of $\theta_{i}$ tested against
100 random binning boundaries (cyan bars). The results show that stochastic
sampling gives a wide spread in weights with a maximum around weight $=5$,
underlining that even when many different boundaries are tested it does not
guarantee good performance.

##  III Polariton dynamics

To investigate the applicability of our method in polaritonic systems, we
start by performing generalized two-dimensional Gross-Pitaevskii (2DGPE)
simulations [see Eqs. (12)-(13) in Appendix A] on five exciton-polariton
condensates in the AFM house configuration using the technique described in
[29]. For each steady-state realization, 100 random binning boundaries of the
condensate phases are applied to find the M3C. We plot the resulting M3C
weights and continuous XY energies (yellow bars) in Figs. 1(d,e), showing that
simulations give a correct M3C weight $=6$ while coming close to the XY ground
state as indicated by the dashed black line. We additionally plot the
stochastically sampled distribution of $H_{XY}$ energies in cyan bars.

Next, we experimentally inject five exciton-polariton condensates in the house
configuration in a semiconductor microcavity, similar to the methods detailed
in [29], and measure the output phase configuration of the interacting
condensates using the techniques described in [46] [see Appendix B]. We record
60 experimental realizations and calculate each time the XY energy of the
condensate network, obtaining the distribution given by the purple bars in
Fig. 1(e). Our observations confirm that interacting polariton condensates
favor phase configurations towards small $H_{XY}$ energies as pointed out in
[22], but might not necessarily reach the ground state and instead get stuck
in local minima of the energy landscape, or converge into a nonstationary
state where $\theta_{i}$ has no meaning. This explains why the theoretical and
experimental distributions are maximal around $H_{XY}\approx-6$ instead of the
ground state -8.7419. Nevertheless, the resulting M3C weights from the
experiment shown in Fig. 1(d) (purple bars) indicate very good performance,
implying that approximate solutions to the continuous variable problem of
minimizing $H_{XY}$ can indeed give good results to the M3C. This important
result underscores that coherent networks of dissipatively coupled
oscillators, like polariton condensates, can potentially perform as heuristic
solvers for the NP-hard M3C.

Figure 2: (a,b,c) Performance of $|\mathcal{V}|=6,12,18$ vertex graphs,
respectively, showing an expected drop with graph size yet still maintaining
peak probability around zero error ($S_{W}=0$). The mean errors are
$\bar{S}_{W}=0.53\%$, $1.38\%$, and $1.86\%$. With more random binning
boundaries tested [different rotations of the triangle in Fig. 1(a)] the
probability of finding the correct solution increases. (d) Mean error
$\bar{S}_{W}$ (whole lines) plotted against the ground state energy gap,
$\text{min}{(H_{XY})}-\text{min}{(H_{T})}$, indicating weak dependence when
the energies between the ternary and continuous spin systems are different.
Shaded area denotes the standard deviation. Here, $\text{min}{(H_{XY})}$ is
estimated using the Stuart Landau network (as opposed to e.g. the basin
hopping method) for computational speediness, whereas $\text{min}{(H_{T})}$ is
found using brute force. Horizontal axis is given in units of
$|\mathcal{V}|\sigma$ where $\sigma$ is the standard deviation of
$\mathcal{W}_{ij}$.

In Figs. 1(f,i) we show a real-space map of the non-resonant laser intensity
used to excite the polariton condensates in simulation and experiment,
respectively. The laser intensity is shaped to form five rings with central
coordinates $\mathbf{r}_{i}=(x_{i},y_{i})$ which keep the condensates
localized at the vertices of the house graph. The couplings $J_{ij}=J<0$
between the condensates which enter into Eqs. (1) and (7) are determined by
choosing specific coordinates $\mathbf{r}_{i}$, as previously studied in Ref.
[29]. In Figs. 1(g,h) and 1(j,k) we show the density $|\Psi(\mathbf{r})|^{2}$
and the phase $\theta(\mathbf{r})=\arg{(\Psi(\mathbf{r}))}$ of the condensate
wavefunction $\Psi(\mathbf{r})$ from simulation and experiment respectively,
where Fig. 1(j) is averaged over many condensate realizations [see Fig. 7(c)
for the single shot real-space realization of the house graph PL]. The clear
interference pattern observed in both experiment [Fig. 1(j)] and simulation
[Fig. 1(g)] implies phase synchronization between condensates. With the
condensates synchronized, we can extract their phases at the location of the
vertices through interferometric techniques [46] such that
$\theta(\mathbf{r}_{i})=\theta_{i}$. The presented phase maps in Figs. 1(h)
and 1(k) both give a ternary phase configuration (overlaid) that matches the
M3C shown in Fig. 1(c).

###  III.1 Benchmarking phase-and-amplitude oscillators

The results shown in Fig. 1 underline the promise of applying dissipative
oscillatory systems to solve the M3C problem. But in order to gain a better
understanding on the quality of the method we look into its statistics by
testing many different graph configurations $\mathcal{G}$ of randomly chosen
weights $\mathcal{W}_{ij}$. In Fig. 2 we test the performance of our method to
solve the M3C by simulating a network of dissipatively coupled phase-and-
amplitude (Stuart-Landau) oscillators [see Eq. (15)] which form a very general
setting of coupled nonlinear oscillators. We compare the weight obtained from
the oscillator network $W_{\text{SL}}$ against the correct maximum/minimum
weight $W_{\text{BF}}^{\text{max(min)}}$ belonging to a 3-cut in the graph
$\mathcal{G}$, found using a brute force method. We define the normalized
error as,

$S_{W}=\frac{W_{\text{BF}}^{\text{max}}-W_{\text{SL}}}{W_{\text{BF}}^{\text{max}}-W_{\text{BF}}^{\text{min}}}.$
(8)

With this metric $S_{W}=0$ means that the system has found the best cut
whereas $S_{W}=1$ is the worst cut. We use 10000 dense random configurations
$\mathcal{W}_{ij}$ (normally distributed with zero mean) and numerically solve
the Stuart-Landau network dynamics from stochastic initial conditions for each
configuration to obtain the steady (synchronized) state, corresponding to
fixed points of $\psi_{n}^{*}\psi_{m}$. Increasing the number of random
binning boundaries tested results in increased performance which can be
intuitively understood from the fact that while Eq. (7) is independent of
global rotation of the spins $\theta_{i}\to\theta_{i}+\phi$ the procedure of
binning the XY spins to evaluate Eq. (1) is not. Therefore several different
orientations of the binning boundary should be attempted in order to obtain
the best value to the M3C. The results show very good performance with mean
error of $\bar{S}_{W}=0.53\%,1.38\%$, and $1.86\%$ for graphs of sizes
$|\mathcal{V}|=6,\ 12$, and $18$, well below today’s best error guarantee of
16.40% in [30] and 19.98% in [52].

We additionally investigate whether the performance of the system depends on
the energy gap between the ground states of the ternary spin system and
continuous $XY$ spin system, i.e. $\text{min}{(H_{XY})}-\text{min}{(H_{T})}$.
Results in Fig. 2(d) show the mean error for different energy gaps, indicating
good performance with no significant dependence on the ground state energy gap
between the two Hamiltonians.

Figure 3: (a) Each small square represents a pixel in a colored image with $i$
labelled and the surrounding patch colored according to RGB color layer.
Radius $r$ shows the maximum distance of pixel connectivity where only pixels
within the ring have an edge weight connecting $\mathcal{V}_{i}$ and
$\mathcal{V}_{j}$, otherwise $\mathcal{W}_{ij}=0$; (b) $m=12$ random pixels
are selected and transformed into a graph (c).

##  IV Image-Segmentation

The M3C problem can be used to segment an image into 3 objects or regions
[54]. When formulating a graph from an image, the vertices and edges represent
pixels and their relative similarity respectively. When selecting vertices,
all pixels or a smaller sample can be included, with connectivity to
neighbours within a given radius $r$, as shown in Fig. 3(a). To define the
edge weights $\mathcal{W}_{ij}$, both local and global properties between two
sampled pixels $i$ and $j$ such as color, brightness, texture, and spatial
proximity [54, 55, 56] can be used through variety of weight estimating
methods [see Eqs. (17-21) in Appendix C] [57, 58] to form a graph with
vertices representing the sampled pixels [Fig. 3(b,c)].

Figure 4: Images of (a) an apple and (g) a tree, with image-segmentation
results (b-f, h-l) using the Stuart-Landau network for methods Eqs. (17)-(21)
respectively, with $m=200$, $r=400$, $q=0.1$ for (b-f) and $q=0.2$ for (h-l)
[see Appendix C]. The cyan, magenta and yellow arrows represent the subsets
$\mathcal{V}_{0,1,2}$ respectively.

We show 3-way image-segmentation for pictures of an apple and a tree by
solving the M3C problem for each graph using the same procedure as described
in Figs. 1-2 with the Stuart-Landau network [see Eq. (15)]. In Fig. 4 we show
the partitioned spins overlaid on example images of an apple and a tree for
five different methods of estimating $\mathcal{W}_{ij}$ that all successfully
locate objects within each image. Methods which combine both local (near
neighbour) pixel and global properties [Figs. 4(c,f,i,l)] best locate a single
object and a background by predominantly cutting the graph into just two
subsets. This is a result of the similar-valued edge weights between pixels
within the red, black and green regions. However, the methods that only
consider local pixel values [Figs. 4(b,e,h,k)] segment the images into
multiple same-spin regimes. This is seen more clearly for the simpler apple
image, where the local methods are able to locate the white background, black
outline, red body and green leaf. As there are four main block colors making
up the apple image but only three subsets for the segmentation, a single
subset is representing multiple objects, such as $\mathbf{\mathcal{V}_{0}}$ in
Fig. 4(b) representing both the red and black regions of the apple. By
considering only global image properties [Figs. 4(d,j)], the images also
segment into object and background by locating the dominant colors of the
image, though some objects are located by multiple subsets. Again, this is an
artefact due to having only three subsets to segment a four object image with.
In Fig. 4(d) for example, the red of the apple is represented by both
$\mathcal{V}_{0}$ and $\mathcal{V}_{2}$ which indicates that the steady state
phases of these oscillators fell about the binning boundary between these two
segmentations. Through adjusting the segmentation parameters and choosing a
specific method, image segmentation using dissipative coupled oscillators can
be achieved to match a wide range of segmentation requirements.

We also consider a simpler colored image with 25 pixels in Fig. 5 to
demonstrate image-segmentation using a planar graph. We first find the image-
segmentation again using the Stuart-Landau network, with the all-to-all
coupling between 25 pixels [see Fig. 5(a)], and a random sample of 5 pixels
[see Fig. 2(b)], which correctly locates the different block colors in the
image. As the sampled pixels are sparsely connected in the latter approach,
the correct image-segmentation can be solved using a small planar graph of
polariton condensates like in Fig. 1 which, through optimization of the M3C,
converge to a phase map (Fig. 5(c)) representing the correct image-
segmentation.

Figure 5: Simple colored image (a,b) showing image-segmentation using the
Stuart-Landau network for the phase-discretized XY ground state of the graphs
following Eq. (20), with $r=4$, $q=0.01$ for (a) $m=25$ and (b) $m=5$. (c)
Shows the condensate phase map $\arg{(\Psi(\mathbf{r}))}$ obtained from 2DGPE
simulation solving the set of pixels given by (b) with transparency
proportional to the polariton density $|\Psi(\mathbf{r})|^{2}$. Projecting the
phase of each condensate into its ternary counterpart is given by the cyan,
magenta and yellow arrows, representing the subsets $\mathcal{V}_{0,1,2}$
respectively. Solid and dashed black lines represent $\mathcal{W}_{ij}=1.0$
and $\mathcal{W}_{ij}=0.8$ (rounded to 1 decimal place) respectively.

##  V Constrained-Via-Minimization

Here we discuss the application of the M3C for constrained-via-minimization
(CVM) in circuit design using coherent networks. In order to optimize space
used in commercial applications, complex circuits are often split over
multiple layers of a circuit board. This is achieved by drilling holes, known
as “vias”, that are lined with a conductive coating, allowing tracks to
connect between multiple layers. These additional vias increase production
time, complexity, and cost, making it desirable to minimize their number.

For CVM, all cells (gray areas in Fig. 6(a)) are pre-placed and vertical and
horizontal tracks are routed with the assumption that all pin connections are
bipartite (i.e., a single track connects two pins), but layer assignment is
not yet performed. Segments of track which overlap are labelled as critical
segments (solid lines) which cannot be on the same layer of circuit board.
Free segments (dashed lines) of track have no overlap with other tracks and
these are the regions in which vias can be placed. We demonstrate CVM through
reducing the task to a M3C problem [34], which we then solve through
simulation on a network of polariton condensates. We demonstrate CVM to a
maximum of three layers of circuit board through reducing the task to a M3C
problem. It is worth noting that this exact CVM problem is tackled in the work
of [35] with the max-2-cut problem, so we know that the minimum number of vias
for the circuit in Fig. 6(a) requires just two circuit board layers, such that
the solution to the max-3-cut problem will only incorporate two out of the
three possible spin subsets.

Figure 6: (a) A circuit with routing between pre-placed cells (gray areas),
with critical segments of tracks numbered and shown in solid black line and
free segments showed in gray dashes; (b) the layout graph of circuit (a) with
critical edges shown in solid black and continuation edges in gray dashes; (c)
the reduced layout graph of (b) with edges labelled
($\alpha_{ij},\beta_{ij}$); (d) the reduced layout graph of (b) with edge
weights $w_{ij}=\alpha_{ij}-\beta_{ij}$, and (e,f) ground state discretized XY
phase for the M3C of graph (d), solved using the Stuart-Landau network and
2DGPEs respectively, with black numbers representing the graph edge weights
$w_{ij}$. The cyan, magenta and yellow arrows represent the subsets
$\mathcal{V}_{0,1,2}$ respectively and in (f) the transparency is proportional
to the polariton density $|\Psi(\mathbf{r})|^{2}$.

For a circuit [see Fig. 6(a)], we define the layout graph
$\mathcal{G}=(\mathcal{V},\mathcal{E})$. Each critical segment is represented
by a vertex in set $\mathcal{V}$, where pairs of vertices are connected by
either a conflict edge $\mathcal{A}$ (when a pair of critical segments cross
paths) or by a continuation edge $\mathcal{B}$ (when a pair of critical
segments are connected by a free segment), such that
$\mathcal{E}=\mathcal{A}\cup\mathcal{B}$. The layout graph of the example
circuit is shown in Fig. 6(b). The reduced layout graph
$\mathcal{R}=(\mathcal{S},\mathcal{T})$ arbitrarily selects a vertex $v_{i}$
in $\mathcal{V}_{i}$ to represent critical region $i$, such that
$\mathcal{S}=\\{v_{i},\dots,v_{z}\\}$. $\mathcal{T}$ contains the edges
linking $v_{i}$ and $v_{j}$ for $i\neq j$, if and only if $\mathcal{G}$
contains a continuous edge connecting some vertex in $\mathcal{V}_{i}$ to some
vertex in $\mathcal{V}_{j}$. As continuation edges do not cross, $\mathcal{R}$
is a planar graph. The edges of $\mathcal{R}$ in Fig. 6(c) have weights
($\alpha_{ij}$, $\beta_{ij}$), as contained in $\mathcal{T}$, such that:

$\alpha_{ij}=$ |  Sum of free segments between $\mathcal{V}_{i}$ and $\mathcal{V}_{j}$ connecting two critical segments with different orientations ;
---|---
$\beta_{ij}=$ |  Sum of free segments between $\mathcal{V}_{i}$ and $\mathcal{V}_{j}$ connecting two critical segments with the same orientation.

The partition of this graph into $\mathcal{S}_{n}$, $n=0,1,2$, corresponds to
the assignment of each critical region to 3-layer $n$. For such a partition,
the number of vias required is

$\text{VIA}(\mathcal{S}_{n})=\sum_{v_{i}v_{j}\in\mathcal{T}_{n}}\alpha_{ij}+\sum_{v_{i}v_{j}\in\delta\mathcal{T}}\beta_{ij}$
(9)

where edges $\mathcal{T}_{n}$ connect critical regions assigned the same layer
and $\delta\mathcal{T}$ connect critical regions assigned different layers. By
defining $A=\sum_{v_{i}v_{j}\in\mathcal{T}}\alpha_{ij}$, then

$\text{VIA}(\mathcal{S}_{n})-A=\sum_{v_{i}v_{j}\in\delta\mathcal{T}}(\beta_{ij}-\alpha_{ij}).$
(10)

As $A$ is invariant on the layer assignment, and by redefining the edge
weights $\mathcal{W}_{ij}=\alpha_{ij}-\beta_{ij}$ as in Fig. 6(d), the problem
is reduced to a M3C problem,

$\max\sum_{v_{i}v_{j}\in\delta\mathcal{T}}\mathcal{W}_{ij}.$ (11)

The phase-discretized XY solution coming from the polariton network for our
example reduced layout graph [see Figs. 6(e,f)] partitions the vertices into
just two subsets, showing that the minimum via configurations only requires
two layers of circuit board. In this example, the minimized number of vias is
2, as $A-\sum_{v_{i}v_{j}\in\delta\mathcal{T}}\mathcal{W}_{ij}=2-0$ and thus a
correct solution has been found by the polariton network.

##  VI Conclusions

We have investigated both theoretically and experimentally the potential of
using nonlinear optical oscillatory networks, specifically exciton-polariton
condensates, in approximating the solutions to the NP-hard max-3-cut
optimization problem. Our study is motivated by recent works showing that
networks of exciton-polariton condensates undergo a gradient descent towards
equalibria of synchronized states with phasor configurations that correlate
with the ground state of the XY Hamiltonian [22, 28, 51, 29]. We exploit this
dynamical feature to approximate the max-3-cut through two methods: First, we
apply a semi-definite relaxation program by Goemans and Williamson known as a
random hyperplane method [5, 30] which projects (bins) the condensate phasors
into ternary decision variables. Second, we use the direct mapping between the
$q=3$ vector Potts model (which we refer to as the “ternary XY model”) and
max-3-cut, as originally presented by Frieze and Jerrum [52], which can then
be approximated using the previously obtained ternary decision variables. Our
study provides experimental evidence that polariton condensate networks can
potentially serve as optical annealers for the max-3-cut, and opens
perspectives on the role of coherent nonlinear optical networks as fast
approximate analogue solvers for complex combinatorial problems. While
interactions between polariton condensates are inherently planar work has
begin to address the possibility of all-to-all connectivity [59, 60, 61]. We
have also studied the more complex applications of image-segmentation and CVM
in circuit design which we heuristically solve using our proposed method and
simulations on both Stuart-Landau oscillator networks and polariton condensate
networks.

Our work goes beyond previously studied Ising machines to solve the max-2-cut
problem [14] by exploiting the continuous phase degree-of-freedom in
oscillatory systems. Their natural tendency in minimizing the XY model can be
applied to solve the M3C problem through Goemans and Williamson inspired semi-
definite program, where we further explore the projection of phases to 2 and 4
bins to solve the max-2-cut and max-4-cut in Appendix D. This technique can be
applied to any dissipative oscillatory system such as laser networks and
photonic condensates, but in this study we have taken steps towards realizing
a liquid light machine of interacting polariton condensates as a practical
coherent network computational device by exploiting the ultrafast temporal
dynamics, parallel interactive nature, and continuous degree of freedom that
can now be readily accessed in state-of-the-art experiments in polariton
lattices [46, 47].

## Acknowledgements

S.L.H., H.S., J.D.T. and P.G.L. acknowledge the support of the UK’s
Engineering and Physical Sciences Research Council (grant EP/M025330/1 on
Hybrid Polaritonics), and S.A. acknowledges the funding of the Russian
Foundation for Basic Research (RFBR) within the joint RFBR and CNR project No.
20-52-7816. H.S. and P.G.L. also acknowledge the European Union’s Horizon 2020
program, through a FET Open research and innovation action under the grant
agreement No. 899141 (PoLLoC). S.L.H. acknowledges the use of the IRIDIS High
Performance Computing Facility, and associated support services at the
University of Southampton, in the completion of this work.

## Data Availability

The data presented in this manuscript can be accessed on the University of
Southampton data repository by following: doi.org/10.5258/SOTON/D2104.

##  Appendix A Polariton Theory

The ground state of the XY model is obtained numerically through previously
studied methods using the generalized GPE equation describing the dynamics of
interacting polariton condensates [22]. The laser driven microcavity system is
simulated close, but above, the condensation threshold where the interacting
condensates (corresponding to their respective laser spots) are found to
spontaneously self-organize into a phase configuration which maximizes the
particle number of the condensate which can be regarded as minimization of an
effective XY model [23]. The polariton condensate wavefunction,
$\Psi(\mathbf{r},t)$, is described by the generalized GPE coupled to a
reservoir of excitons $n(\mathbf{r},t)$ which experience bosonic stimulated
scattering into the condensates [62].

$i\frac{\partial\Psi}{\partial
t}=\bigg{[}-\frac{\hbar\nabla^{2}}{2m}+G\left(n+\frac{P(\mathbf{r})}{W}\right)\\\
+\alpha|\Psi|^{2}+\frac{i}{2}\left(Rn-\gamma\right)\bigg{]}\Psi,$ (12)

$\displaystyle\frac{\partial n}{\partial t}$
$\displaystyle=-\left(\Gamma+R|\Psi|^{2}\right)n+P(\mathbf{r}).$ (13)

Here, $m$ is the effective mass of a polariton in the lower dispersion branch,
$\alpha$ is the interaction strength of two polaritons in the condensate, $G$
is the polariton-reservoir interaction strength, $R$ is the rate of stimulated
scattering of polaritons into the condensate from the reservoir, $\gamma$ is
the polariton decay rate, $\Gamma$ is the decay rate of the the reservoir
excitons, $W$ quantifies the population ratio between low-momentum excitons
that scatter into the condensate and those which reside at higher momenta (so
called inactive excitons), and $P(\mathbf{r})$ is the non-resonant continuous
wave (CW) pump profile given by

$P(\mathbf{r})=P_{0}\sum_{i}p(\mathbf{r}-\mathbf{r}_{i}).$ (14)

Here, $P_{0}$ denotes the laser power density and the function $p(\mathbf{r})$
corresponds to the 2D annular shaped profile of a single laser incident onto
the microcavity plane and the coordinates $\mathbf{r}_{i}$ are the locations
of the vertices in the polariton graph. The annular shaped pump profiles
optically trap each condensate, yet allow for coherent transport of particles
between nearest neighbours, consistent with those described in [29].

In the 2DGPE simulations, the parameters are taken such that the polariton
mass and lifetime are based on the properties of a laboratory InGaAs
microcavity sample: $m=0.28$ meV ps2 $\upmu$m-2 and $\gamma=\frac{1}{5.5}$
ps-1. We choose values of interaction strengths typical of InGaAs based
systems: $\hbar\alpha=7$ $\upmu$eV $\upmu$m2, $G=10\alpha$. The reservoir
decay rate is taken comparable to the condensate decay rate $\Gamma=\gamma$
due to fast thermalization to the exciton background. The final two parameters
are then found by fitting to experimental results where we use the values
$\hbar R=98.9$ $\upmu$eV $\upmu$m2, and $W=0.035$ ps-1.

We also consider a more simple and general model which does not depend on
complicated spatial degrees of freedom. It describes dissipative coupling
between nonlinear oscillators $\psi_{n}(t)$ and is known as a Stuart-Landau
network,

$\frac{d\psi_{n}}{dt}=\left[P-|\psi_{n}|^{2}\right]\psi_{n}+\sum_{m}J_{nm}\psi_{m}.$
(15)

Here, $P$ denotes the “gain” of each oscillator and $J_{nm}$ is the coupling
strength. Such networks have shown good performance in minimizing the XY
Hamiltonian [28] and share similarities to adiabatic bifurcation networks
[63]. Our simulations always start from random initial conditions and once the
oscillators $\psi_{n}$ have converged to a steady state for a given $P$ we
extract their phases $\theta_{n}=\text{arg}{(\psi_{n})}$ to obtain the
approximate ground state energy of Eq. (7). We perform numerical integration
of Eqs. (12)-(13) and Eq. (15) in time using a linear multistep method. We
point out that, to the best of our knowledge, there is no analytical estimate
on the “performance guarantee” of the simulated Stuart-Landau network in
finding the XY ground state and therefore the performance guarantee of finding
the M3C cannot be ascertained at the current stage except through numerical
methods.

In the Stuart-Landau network simulations we increase $P$ linearly in time from
$P(t=0)=-\lambda_{\text{max}}$ to $P(t=T)=\lambda_{\text{max}}$ where $T$ is
the total integration time and $\lambda_{\text{max}}>0$ is the largest
eigenvalue of the coupling matrix
$\mathbf{J}=\\{J_{nm}\\}\in\mathbb{R}^{N\times N}$. This physically replicates
a slow ramp-up in laser power beyond the polariton condensation threshold
resulting in measurable photoluminescence from the system. We note that in the
terms of amplitude oscillator models, the condensation threshold is a
bifurcation point marking the departure of the condensate (the oscillator)
from the stable $|\psi_{n}|=0$ solution. In Fig. 8 we set $P=0$ to underline
that our results do not critically depend on slow pump rise times.

##  Appendix B Experiment

The GaAs-based microcavity is cooled down to 4 K in a closed-cycle helium
crysotation and pumped non-resonantly with a CW single-mode laser at a
negative exciton-photon detuning of -4.2 meV. In order to shape the excitation
profile we utilize a programmable reflective phase-only spatial light
modulator. The laser pump pattern is focused onto the sample with a 50x
microscope objective of NA=0.42. The separation distance between pumping rings
(20.5 $\upmu$m) and their diameter (8.1 $\upmu$m) is chosen such that any pair
of the nearest trapped polariton condensates in a square geometry demonstrates
AFM coupling [see Figs. 7(a,b)]. The excitation geometry defines the range of
the excitation intensity per single pumping ring, where the trapped
condensates maintain AFM coupling. Therefore, for the house configuration the
excitation power is proportionally increased to ensure AFM coupling between
nearest condensates. The results presented in Fig. 1 correspond to excitation
conditions supporting only single energy (i.e. stationary) states above
condensation threshold. Time-averaged measurements of real space polariton
photoluminescence are performed under CW excitation, acousto-optically
modulated in time with square pulses at a frequency of 5 kHz and duty cycle of
1$\%$. To implement the relative phase readout between the nodes in the house
configuration we utilize a homodyne interferometric technique [46]. Each
reconstructed phase map is extracted from a single-shot measurement in a Mach-
Zenhder interferometer. The excitation pulsewidth in all single-shot
measurements is 100 $\upmu$s.

The results of the single shot measurements are shown in Figs. 1(d,e), on the
purple histograms of the M3C weights and XY energy respectively. Here, in Fig.
7(c) we show the single shot polariton PL in real space for the house graph
that gives the best result in minimizing the XY energy of the system. The
corresponding phase map for this realization is shown in Fig. 1(k).

Figure 7: Time integrated (a) real-space and (b) reciprocal-space polariton PL
for the square cell of trapped condensates sustaining AFM coupling. (c) Single
shot realization of the polariton PL in real space for the AFM house graph.

##  Appendix C Image-Segmentation

For a colored image, we randomly sample $m$ pixels. For each color layer, we
consider each sampled pixel $i$ and define a local patch encompassing it and
its (up to) eight-way nearest neighbours, as depicted in Fig. 3(a). We define
the sum of all $N$ pixels in the patch as,

$p_{i}=\sum_{i=1}^{N}i,$ (16)

where $p_{i}$ is calculated for all pixels in the image and are normalized
over each color layer [58]. In addition to these local weights, we assign a
global color weight, $c_{i}=M_{A}/M$, to each pixel according to its color
frequency in the image, where $A$ is the color of pixel $i$, $M_{A}$ is the
number of pixels with color $A$, and $M$ is the total number of pixels. We
define the edge weight connecting pixels $i$ and $j$ as $\mathcal{W}_{ij}$,
which is equal to 0 if the pixel separation is greater than $r$. Otherwise, we
define five methods for enumerating $\mathcal{W}_{ij}$ between two sampled
pixels based on a variety of techniques described in literature [54, 58, 64,
57] that utilize the difference between their pixel values, patch values and
global weighting of their colors, as well as the difference between the
product of their global color weighting with the patch or pixel values,

$\hskip 23.84592pt\text{Method 1:}\hskip
4.33382pt\mathcal{W}_{ij}=\exp\bigg{(}\frac{|p_{i}-p_{j}|^{q}}{\sigma}\bigg{)};$
(17)

$\hskip 34.69038pt\text{Method 2:}\hskip
4.33382pt\mathcal{W}_{ij}=\exp\bigg{(}\frac{|c_{i}p_{i}-c_{j}p_{j}|^{q}}{\sigma}\bigg{)};$
(18)

$\hskip 21.68231pt\text{Method 3:}\hskip
4.33382pt\mathcal{W}_{ij}=\exp\bigg{(}\frac{|c_{i}-c_{j}|^{q}}{\sigma}\bigg{)};$
(19)

$\text{Method 4:}\hskip
4.33382pt\mathcal{W}_{ij}=\exp\bigg{(}\frac{|i-j|^{q}}{\sigma}\bigg{)};$ (20)

$\hskip 28.18636pt\text{Method 5:}\hskip
4.33382pt\mathcal{W}_{ij}=\exp\bigg{(}\frac{|ic_{i}-jc_{j}|^{q}}{\sigma}\bigg{)}.$
(21)

Here, $\sigma$ is the standard deviation in brightness across each patch, $q$
is a free parameter and $\mathcal{W}_{ij}$ is averaged over each color layer.
To find the results of the image-segmentation, we find the M3C of the graph
$\mathcal{G}=(\mathcal{V},\mathcal{W})$, with $|\mathcal{V}|=m$ vertices
representing the sampled pixels with weights $\mathcal{W}_{ij}$ connecting
vertices $i$ and $j$. We will show that the partition of the continuous phases
in a polariton condensate graph (just like demonstrated in Fig. 1 for the
house graph) into a ternary phase configuration, acting as decision variables
for the M3C, can segment different objects within an image.

##  Appendix D Other maximum cuts

It is worth noting that the continuous-phase spins of the XY model can also be
easily projected onto binary angles 0 and $\pi$ (i.e. to 1-dimensional spins)
which leads to an Ising energy function whose ground state can be mapped to
the max-2-cut (M2C) problem [35]. Such a binary projection strategy has been
previously explored with coherent photonic Ising machines [12, 13, 14, 15, 16,
17, 20]. However, to the best of our knowledge, continuous phase photonic
machines have not been used to explore the mapping between the ground state of
the ternary spin Hamiltonian and the M3C problem (i.e., projection onto [0,
2$\pi$/3, 4$\pi$/3] spins).

Figure 8: Obtained weights of (a) max-2-cut, (b) max-3-cut and (c) max-4-cut
from 160 random graphs using the Stuart-Landau model ($W_{SL}$) with
$|\mathcal{V}|=10$ and 100 unique binning boundaries and the brute force
method ($W_{BF}$). Black line indicates $W_{SL}=W_{BF}$. Comparison of the
mean error $\bar{S}_{W}$ obtained by averaging over 160 random graphs of
different max-$k$-cut tasks as a function of (d) network size $|\mathcal{V}|$
(tested for 100 unique binning boundaries) and (e) number of binning
boundaries with $|\mathcal{V}|=10$.

We show the weight of cuts solved using the Stuart-Landau model compared to
the brute force method for the M2C and M3C in Fig. 8(a,b). In the former case,
the error between $W_{SL}$ and $W_{BF}$ is negligible. As there are only two
possible phase binnings for the M2C problem, every oscillator has a window of
$180^{\circ}$ in which to be correctly binarized. In the latter case, this
segment is reduced to $120^{\circ}$, leading to less room for error in the
binning process and thus an increase in error between $W_{SL}$ and $W_{BF}$.

Furthermore, as we pointed out in Sec. II, mapping the ground state
configurations of spin Hamiltonians to the max-$k$-cut problem requires the
spins to belong to the corners of the unit $\mathbb{R}^{k-1}$ simplex. This
means that the mutual angle between any two spins is either 0 or a constant
value (e.g., for the triangle it is $120^{\circ}$ and for the tetrahedron it
is $\approx 70.53^{\circ}$). For this reason, projecting the XY ground state
spins to quaternary decision variables corresponding to the angles
$\\{0,\pi/2,\pi,3\pi/2\\}$ on the unit circle, with the possible set of values
$\cos{(\theta_{i}-\theta_{j})}\in$ $\left\\{1,0,-1\right\\}$, is expected to
perform worse in solving the max-4-cut (M4C). Indeed, Fig. 8(c) shows that the
weights found by the Stuart-Landau network have a greater spread away from the
optimal weight. The difference between the $k=2,3,4$ max-$k$-cuts becomes even
more apparent when we plot the normalized mean error in Fig. 8(d) as a
function of graph size. Interestingly, the error for the M2C stays practically
negligible indicating that Stuart-Landau systems can compete with photonic
Ising machines. The larger error for the M3C is expected due to the higher
partition complexity (i.e., 3 binning options) but amazingly stays practically
invariant with graph size in contrast to the growing error of the M4C. We
finally point out that even though the error of the M4C is increasing with
system size it remained below $\bar{S}_{W}<0.05$ for $|\mathcal{V}|=14$ vertex
graphs which possess $4^{13}=67108864$ different partitions. This opens an
exciting perspective in using the condensate network phase dynamics to go
beyond optimizing the M3C even if a direct mapping no longer exists. We also
plot the mean error of the different max-$k$-cut problems while scanning the
number of binning boundaries tested [Fig. 8(e)], which shows that the error
converges to some minimum asymptotic value as the number of binning
boundaries. Notably, the asymptotic error increases with the number of cuts.

## References

  * Ibsen-Jensen  _et al._ [2015] R. Ibsen-Jensen, K. Chatterjee, and M. A. Nowak, Computational complexity of ecological and evolutionary spatial dynamics, Proceedings of the National Academy of Sciences 112, 15636 (2015).
  * Zhang and Dincer [2016] X. Zhang and I. Dincer,  _Energy Solutions to Combat Global Warming_ (Springer, 2016).
  * Jayaraj  _et al._ [2016] P. B. Jayaraj, K. Rahamathulla, and G. Gopakumar, A gpu based maximum common subgraph algorithm for drug discovery applications, in  _2016 IEEE International Parallel and Distributed Processing Symposium Workshops (IPDPSW)_ (2016) pp. 580–588.
  * Pierce and Winfree [2002] N. A. Pierce and E. Winfree, Protein Design is NP-hard, Protein Engineering, Design and Selection 15, 779 (2002).
  * Goemans and Williamson [1995] M. X. Goemans and D. P. Williamson, Improved approximation algorithms for maximum cut and satisfiability problems using semidefinite programming, J. ACM 42, 1115–1145 (1995).
  * Zhang and Huang [2006] S. Zhang and Y. Huang, Complex quadratic optimization and semidefinite programming, SIAM Journal on Optimization 16, 871 (2006).
  * Such  _et al._ [2017] F. P. Such, V. Madhavan, E. Conti, J. Lehman, K. O. Stanley, and J. Clune, Deep neuroevolution: Genetic algorithms are a competitive alternative for training deep neural networks for reinforcement learning, CoRR abs/1712.06567 (2017), arXiv:1712.06567 .
  * Yang [2014] X.-S. Yang,  _Nature-Inspired Optimization Algorithms_ (Elsevier, Amsterdam, Netherlands, 2014).
  * Tovey [2018] C. A. Tovey, Nature-Inspired Heuristics: Overview and Critique, in  _Recent Advances in Optimization and Modeling of Contemporary Problems_, INFORMS TutORials in Operations Research (INFORMS, 2018) pp. 158–192.
  * Aspuru-Guzik and Walther [2012] A. Aspuru-Guzik and P. Walther, Photonic quantum simulators, Nature Physics 8, 285 (2012).
  * Sun  _et al._ [2015] C. Sun, M. T. Wade, Y. Lee, J. S. Orcutt, L. Alloatti, M. S. Georgas, A. S. Waterman, J. M. Shainline, R. R. Avizienis, S. Lin, et al., Single-chip microprocessor that communicates directly using light, Nature 528, 534 (2015).
  * Marandi  _et al._ [2014] A. Marandi, Z. Wang, K. Takata, R. L. Byer, and Y. Yamamoto, Network of time-multiplexed optical parametric oscillators as a coherent ising machine, Nature Photonics 8, 937 (2014).
  * McMahon  _et al._ [2016] P. L. McMahon, A. Marandi, Y. Haribara, R. Hamerly, C. Langrock, S. Tamate, T. Inagaki, H. Takesue, S. Utsunomiya, K. Aihara, et al., A fully programmable 100-spin coherent Ising machine with all-to-all connections, Science 354, 614 (2016).
  * Inagaki  _et al._ [2016a] T. Inagaki, Y. Haribara, K. Igarashi, T. Sonobe, S. Tamate, T. Honjo, A. Marandi, P. L. McMahon, T. Umeki, K. Enbutsu, et al., A coherent Ising machine for 2000-node optimization problems, Science 354, 603 (2016a).
  * Inagaki  _et al._ [2016b] T. Inagaki, K. Inaba, R. Hamerly, K. Inoue, Y. Yamamoto, and H. Takesue, Large-scale ising spin network based on degenerate optical parametric oscillators, Nature Photonics 10, 415 (2016b).
  * Kyriienko  _et al._ [2019] O. Kyriienko, H. Sigurdsson, and T. C. H. Liew, Probabilistic solving of $np$-hard problems with bistable nonlinear optical networks, Phys. Rev. B 99, 195301 (2019).
  * Pierangeli  _et al._ [2019] D. Pierangeli, G. Marcucci, and C. Conti, Large-scale photonic ising machine by spatial light modulation, Phys. Rev. Lett. 122, 213902 (2019).
  * Böhm  _et al._ [2019] F. Böhm, G. Verschaffelt, and G. Van der Sande, A poor man’s coherent ising machine based on opto-electronic feedback systems for solving optimization problems, Nature Communications 10, 3538 (2019).
  * Luo  _et al._ [2020] S. Luo, L. Liao, Z. Zhang, J. Wang, X. Shen, and Z. Chen, Classical Spin Chains Mimicked by Room-Temperature Polariton Condensates, Phys. Rev. Applied 13, 044052 (2020).
  * Roques-Carmes  _et al._ [2020] C. Roques-Carmes, Y. Shen, C. Zanoci, M. Prabhu, F. Atieh, L. Jing, T. Dubček, C. Mao, M. R. Johnson, V. Čeperić, et al., Heuristic recurrent algorithms for photonic ising machines, Nature Communications 11, 249 (2020).
  * Nixon  _et al._ [2013] M. Nixon, E. Ronen, A. A. Friesem, and N. Davidson, Observing Geometric Frustration with Thousands of Coupled Lasers, Phys. Rev. Lett. 110, 184102 (2013).
  * Berloff  _et al._ [2017] N. G. Berloff, M. Silva, K. Kalinin, A. Askitopoulos, J. D. Töpfer, P. Cilibrizzi, W. Langbein, and P. G. Lagoudakis, Realizing the classical XY Hamiltonian in polariton simulators, Nature Materials 16, 1120 (2017).
  * Lagoudakis and Berloff [2017] P. G. Lagoudakis and N. G. Berloff, A polariton graph simulator, New J. Phys. 19, 125008 (2017).
  * Takeda  _et al._ [2017] Y. Takeda, Y. Takeda, S. Tamate, Y. Yamamoto, Y. Yamamoto, H. Takesue, T. Inagaki, S. Utsunomiya, and S. Utsunomiya, Boltzmann sampling for an XY model using a non-degenerate optical parametric oscillator network, in  _Frontiers in Optics 2017 (2017), paper FM4E.3_ (Optical Society of America, 2017) p. FM4E.3.
  * Gershenzon  _et al._ [2020] I. Gershenzon, G. Arwas, S. Gadasi, C. Tradonsky, A. Friesem, O. Raz, and N. Davidson, Exact mapping between a laser network loss rate and the classical XY Hamiltonian by laser loss control, Nanophotonics , 20200137 (2020).
  * Xu  _et al._ [2020] X.-Y. Xu, X.-L. Huang, Z.-M. Li, J. Gao, Z.-Q. Jiao, Y. Wang, R.-J. Ren, H. P. Zhang, and X.-M. Jin, A scalable photonic computer solving the subset sum problem, Science Advances 6, 10.1126/sciadv.aay5853 (2020).
  * Tradonsky  _et al._ [2019] C. Tradonsky, I. Gershenzon, V. Pal, R. Chriki, A. A. Friesem, O. Raz, and N. Davidson, Rapid laser solver for the phase retrieval problem, Science Advances 5, 10.1126/sciadv.aax4530 (2019).
  * Kalinin and Berloff [2018a] K. P. Kalinin and N. G. Berloff, Global optimization of spin hamiltonians with gain-dissipative systems, Scientific Reports 8, 17791 (2018a).
  * Harrison  _et al._ [2020] S. L. Harrison, H. Sigurdsson, and P. G. Lagoudakis, Synchronization in optically trapped polariton Stuart-Landau networks, Phys. Rev. B 101, 155402 (2020).
  * Goemans and Williamson [2004] M. X. Goemans and D. P. Williamson, Approximation algorithms for max-3-cut and other problems via complex semidefinite programming, Journal of Computer and System Sciences 68, 442 (2004).
  * Garey and Johnson [1979] M. R. Garey and D. S. Johnson,  _Computers and Intractability: A Guide to the Theory of NP-Completeness_ (W. H. Freeman & Co., USA, 1979).
  * Harary [1959] F. Harary, On the measurement of structural balance, Behavioral Science 4, 316 (1959).
  * Barahona [1982] F. Barahona, On the computational complexity of Ising spin glass models, J. Phys. A: Math. Gen. 15, 3241 (1982).
  * Harary  _et al._ [2002] F. Harary, M.-H. Lim, and D. C. Wunsch, Signed graphs for portfolio analysis in risk management, IMA J Management Math 13, 201 (2002).
  * Barahona  _et al._ [1988] F. Barahona, M. Grötschel, M. Jünger, and G. Reinelt, An Application of Combinatorial Optimization to Statistical Physics and Circuit Layout Design, Operations Research 36, 493 (1988).
  * de Sousa  _et al._ [2013] S. de Sousa, Y. Haxhimusa, and W. G. Kropatsch, Estimation of Distribution Algorithm for the Max-Cut Problem, in  _Graph-Based Representations in Pattern Recognition_, Lecture Notes in Computer Science, edited by W. G. Kropatsch, N. M. Artner, Y. Haxhimusa, and X. Jiang (Springer, Berlin, Heidelberg, 2013) pp. 244–253.
  * Barahona [1996] F. Barahona, Network Design Using Cut Inequalities, SIAM J. Optim. 6, 823 (1996).
  * White [1973] A. T. White, Chapter 8 map-coloring problems, in  _Graphs, Groups and Surfaces_, North-Holland Mathematics Studies, Vol. 8, edited by A. T. White (North-Holland, 1973) pp. 101 – 124\.
  * Deb  _et al._ [2001] S. K. Deb, B. Bhattacharyya, and S. K. Sorkhel, Development of intelligent mathematical modeling for facilities layout design, Proceedings of the National Conference on Mathematical and Computational Models , 235 (2001).
  * Zhou  _et al._ [2020] L. Zhou, S.-T. Wang, S. Choi, H. Pichler, and M. D. Lukin, Quantum Approximate Optimization Algorithm: Performance, Mechanism, and Implementation on Near-Term Devices, Phys. Rev. X 10, 021067 (2020).
  * Yamamoto  _et al._ [2017] Y. Yamamoto, K. Aihara, T. Leleu, K.-i. Kawarabayashi, S. Kako, M. Fejer, K. Inoue, and H. Takesue, Coherent Ising machines—optical neural networks operating at the quantum limit, npj Quantum Information 3, 49 (2017).
  * Hamerly  _et al._ [2019] R. Hamerly, T. Inagaki, P. L. McMahon, D. Venturelli, A. Marandi, T. Onodera, E. Ng, C. Langrock, K. Inaba, T. Honjo, et al., Experimental investigation of performance differences between coherent Ising machines and a quantum annealer, Science Advances 5, 10.1126/sciadv.aau0823 (2019).
  * Vretenar  _et al._ [2021] M. Vretenar, B. Kassenberg, S. Bissesar, C. Toebes, and J. Klaers, Controllable Josephson junction for photon Bose-Einstein condensates, Phys. Rev. Research 3, 023167 (2021), publisher: American Physical Society.
  * Reifenstein  _et al._ [2021] S. Reifenstein, S. Kako, F. Khoyratee, T. Leleu, and Y. Yamamoto, Coherent Ising Machines with Optical Error Correction Circuits, Advanced Quantum Technologies 4, 2100077 (2021).
  * Reddy  _et al._ [2022] A. N. K. Reddy, A. N. K. Reddy, S. Mahler, S. Mahler, A. Goldring, V. Pal, A. A. Friesem, and N. Davidson, Phase locking of lasers with Gaussian coupling, Opt. Express, OE 30, 1114 (2022), .
  * Alyatkin  _et al._ [2020] S. Alyatkin, J. D. Töpfer, A. Askitopoulos, H. Sigurdsson, and P. G. Lagoudakis, Optical control of couplings in polariton condensate lattices, Phys. Rev. Lett. 124, 207402 (2020).
  * Töpfer  _et al._ [2021] J. D. Töpfer, I. Chatzopoulos, H. Sigurdsson, T. Cookson, Y. G. Rubo, and P. G. Lagoudakis, Engineering spatial coherence in lattices of polariton condensates, Optica 8, 106 (2021).
  * Ohadi  _et al._ [2016] H. Ohadi, R. Gregory, T. Freegarde, Y. Rubo, A. Kavokin, N. Berloff, and P. Lagoudakis, Nontrivial Phase Coupling in Polariton Multiplets, Phys. Rev. X 6, 031032 (2016), publisher: American Physical Society.
  * Töpfer  _et al._ [2020] J. D. Töpfer, H. Sigurdsson, L. Pickup, and P. G. Lagoudakis, Time-delay polaritonics, Communications Physics 3, 2 (2020).
  * Cilibrizzi  _et al._ [2014] P. Cilibrizzi, A. Askitopoulos, M. Silva, F. Bastiman, E. Clarke, J. M. Zajac, W. Langbein, and P. G. Lagoudakis, Polariton condensation in a strain-compensated planar microcavity with InGaAs quantum wells, Appl. Phys. Lett. 105, 191118 (2014), publisher: American Institute of Physics.
  * Kalinin and Berloff [2018b] K. P. Kalinin and N. G. Berloff, Simulating ising and $n$-state planar potts models and external fields with nonequilibrium condensates, Phys. Rev. Lett. 121, 235302 (2018b).
  * Frieze and Jerrum [1997] A. Frieze and M. Jerrum, Improved approximation algorithms for maxk-cut and max bisection, Algorithmica 18, 67 (1997).
  * Wales and Doye [1997] D. J. Wales and J. P. K. Doye, Global optimization by basin-hopping and the lowest energy structures of lennard-jones clusters containing up to 110 atoms, The Journal of Physical Chemistry A 101, 5111 (1997).
  * Felzenszwalb and Huttenlocher [2004] P. F. Felzenszwalb and D. P. Huttenlocher, Efficient Graph-Based Image Segmentation, International Journal of Computer Vision 59, 167 (2004).
  * Rouco  _et al._ [2016] J. Rouco, E. Azevedo, and A. Campilho, Automatic Lumen Detection on Longitudinal Ultrasound B-Mode Images of the Carotid Using Phase Symmetry, Sensors 16, 350 (2016).
  * Martin  _et al._ [2001] D. Martin, C. Fowlkes, D. Tal, and J. Malik, A database of human segmented natural images and its application to evaluating segmentation algorithms and measuring ecological statistics, Proceedings Eighth IEEE International Conference on Computer Vision. ICCV 2001 2, 416 (2001).
  * Boykov and Funka-Lea [2006] Y. Boykov and G. Funka-Lea, Graph Cuts and Efficient N-D Image Segmentation, Int J Comput Vision 70, 109 (2006).
  * Wang  _et al._ [2013] X. Wang, C. Zhu, C.-E. Bichot, and S. Masnou, Graph-based Image Segmentation Using Weighted Color Patch, in  _IEEE International Conference on Image Processing (ICIP)_ (Melbourne, Australia, 2013) pp. 4064–4068.
  * Sigurdsson  _et al._ [2019] H. Sigurdsson, O. Kyriienko, K. Dini, and T. C. H. Liew, All-to-all intramodal condensate coupling by multifrequency excitation of polaritons, ACS Photonics 6, 123 (2019).
  * Kalinin  _et al._ [2020] K. P. Kalinin, A. Amo, J. Bloch, and N. G. Berloff, Polaritonic xy-ising machine, Nanophotonics 9, 4127 (2020).
  * Harrison  _et al._ [2021] S. L. Harrison, H. Sigurdsson, and P. G. Lagoudakis, Minor embedding with Stuart-Landau oscillator networks, arXiv e-prints , arXiv:2109.10142 (2021), arXiv:2109.10142 [quant-ph] .
  * Wouters and Carusotto [2007] M. Wouters and I. Carusotto, Excitations in a nonequilibrium bose-einstein condensate of exciton polaritons, Phys. Rev. Lett. 99, 140402 (2007).
  * Goto  _et al._ [2019] H. Goto, K. Tatsumura, and A. R. Dixon, Combinatorial optimization by simulating adiabatic bifurcations in nonlinear hamiltonian systems, Science Advances 5, 10.1126/sciadv.aav2372 (2019).
  * Peng  _et al._ [2013] B. Peng, L. Zhang, and D. Zhang, A survey of graph theoretical approaches to image segmentation, Pattern Recognition 46, 1020 (2013).