# Demonstration of weighted graph optimization on a Rydberg atom array using
local light-shifts

A. G. de Oliveira<EMAIL_ADDRESS>Department of Physics and SUPA,
University of Strathclyde, Glasgow G4 0NG, UK  E. Diamond-Hitchcock
Department of Physics and SUPA, University of Strathclyde, Glasgow G4 0NG, UK
D. M. Walker  Department of Physics and SUPA, University of Strathclyde,
Glasgow G4 0NG, UK  M. T. Wells-Pestell  Department of Physics and SUPA,
University of Strathclyde, Glasgow G4 0NG, UK  G. Pelegrí  Department of
Physics and SUPA, University of Strathclyde, Glasgow G4 0NG, UK  C.J. Picken
M Squared Lasers Limited, 1 Kelvin Campus, West of Scotland Science Park,
Glasgow, G20 0SP, UK  G.P.A. Malcolm  M Squared Lasers Limited, 1 Kelvin
Campus, West of Scotland Science Park, Glasgow, G20 0SP, UK  A. J. Daley
Department of Physics and SUPA, University of Strathclyde, Glasgow G4 0NG, UK
J. Bass  Department of Physics and SUPA, University of Strathclyde, Glasgow G4
0NG, UK  J. D. Pritchard<EMAIL_ADDRESS>Department of
Physics and SUPA, University of Strathclyde, Glasgow G4 0NG, UK

###### Abstract

Neutral atom arrays have emerged as a versatile platform towards scalable
quantum computation and optimization. In this paper we present first
demonstrations of weighted graph optimization on a Rydberg atom array using
annealing with local light-shifts. We verify the ability to prepare weighted
graphs in 1D and 2D arrays, including embedding a five vertex non-unit disk
graph using nine physical qubits. We find common annealing ramps leading to
preparation of the target ground state robustly over a substantial range of
different graph weightings. This work provides a route to exploring large-
scale optimization of non-planar weighted graphs relevant for solving relevant
real-world problems.

Optimization of combinatorially hard problems has been identified as an early
potential application of quantum computing hardware [1], with significant
effort invested in developing protocols such as quantum annealing algorithms
(QAA) [2, 3, 4, 5] or variational based approaches such as Quantum Approximate
Optimization Algorithms [6, 7]. Despite these efforts, hardware capable of
demonstrating a practical quantum advantage in this domain remains elusive [8,
9, 10, 11].

Neutral atom quantum computers based on arrays of individual optical tweezers
[12, 13, 14, 15] provide a scalable and versatile platform for quantum
computing, with the ability to generate arrays of over 1000 qubits [16, 17,
18, 19] and perform high-fidelity single- [20] and two-qubit [21, 22, 23] gate
operations, enabling realization of small scale quantum algorithms [24]. This
can be extended to implementation of logical qubit operations [25] by
exploiting dynamical qubit reconfiguration [26]. Alongside digital operations,
neutral atom arrays provide access to programmable spin models for performing
quantum simulation [27] enabling exploration of phase diagrams and dynamics in
1D [28] and 2D [29, 30] and observation of novel topological phenomena [31,
32].

The ability to implement programmable spin models extends the utility of
neutral atom arrays to the regime of analogue quantum computing [33], where in
the context of graph optimization coupling atoms to the Rydberg states enables
native embedding of the Maximum Independent Set (MIS) problem for planar 2D
unit-disk graphs (UDG) [34, 35]. Here the Rydberg blockade mechanism [36]
prevents more than a single Rydberg excitation occuring within a radius
$r<r_{\mathrm{B}}$, enforcing the constraint that no two connected vertices
should be simultaneously part of the same MIS configuration. Experimental
demonstrations of solving MIS on 2D-UDG using QAA with hundreds of qubits
demonstrate the feasibility of this approach [37, 38], however recent analysis
shows for the native Kings-graph classical solvers are able to find solutions
to graphs with 1000 vertices in minutes [39]. One approach for extending
beyond planar 2D graph is to use ancillary qubits to create quantum wires [40]
which can be routed through neighboring trapping planes in 3D to implement
non-local graph connectivity [41]. This approach has been used to demonstrate
embedded integer factorization [42] and quadratic unconstrained binary
optimization (QUBO) problems [43] onto atom arrays using only global control.

Figure 1: Weighted graph optimization with local light-shifts. (a) The target
graph problem is embedded onto a neutral atom array as a Maximum Weighted
Independent Set (MWIS) on a Unit Disk Graph (UDG), using global Rydberg
excitation lasers ($\Omega_{1,2}$) and an additional set of tweezer beams for
implementing local light-shifts with a programmable shift on each qubit site.
The Rydberg blockade mechanism prevents two atoms within distance $r_{B}$
being simultaneously excited. (b) Atoms without light-shift experience an
effective two-photon Rabi frequency $\Omega$ and global detuning $\Delta$. The
light-shift tweezers induce AC shift
$\delta_{\mathrm{AC}}^{i}=w_{i}\delta_{\mathrm{AC}}$ on atom $i$ to implement
a local detuning $\Delta_{i}=\Delta+\delta_{\mathrm{AC}}^{i}$ with relative
weight $w_{i}$. (c) Annealing protocol using a two-stage process with the
global Rydberg laser ramped from an initial negative $\Delta$ to resonance,
then positive detunings defined by controlling power in light-shift tweezers
to scale $\delta_{\mathrm{AC}}$ resulting in local detunings with fixed
relative ratios defined by weightings $w_{i}$. (d) Ground-state solutions
obtained via closed-loop optimization of the annealing profile to minimize the
expectation value of the classical MWIS cost function.

Another approach to extend the range of problems that can be encoded using
neutral atom arrays is to introduce local control fields. This permits the
introduction of weighting to enable solution of Maximum Weighted Independent
Set (MWIS) problems on a 2D UDG, providing a route to map a variety of graph
and QUBO problems onto neutral atom arrays using this MWIS-UDG encoding [44]
with at worst a quadratic overhead in physical qubit number. Similar scaling
is obtained using a related parity-based encoding for graph problems [45],
whilst moving beyond the blockade regime to remove the UDG restriction offers
prospects for linear scalings in qubit overhead [46]. To date, local light-
shifts have been used on neutral atom arrays in the field of quantum
simulation for generating large-scale entangled states [47] and implementing
fast, arbitrary local control of interacting trimers [48].

In this paper we present first demonstrations of weighted graph optimization
on a Rydberg atom array using annealing with local light-shifts. Our approach
for weighting combines holographic techniques for generating programmable
local light-shifts compatible with use on an arbitrary array geometry with the
ability to implement graph annealing using global laser control fields,
simplifying the number of available parameters when performing optimization.
We verify the ability to prepare weighted graphs on both 1D and 2D arrays,
including embedding a 5-qubit non-UDG using 9 physical qubits. In each case we
find a common annealing ramp able to prepare the target ground state solutions
for arbitrary graph weightings. This work demonstrates the feasibility of
implementing weighted optimization using neutral atom arrays, and provides a
route to embedding non-planar graphs solving relevant real-world problems
[49].

_Weighted Graph Optimization_.— The Maximum Weighted Independent Set (MWIS)
problem on graph $G(V,E)$ comprized of $E$ edges and $V$ vertices, each of
weight $w_{i}>0$, is defined as finding the subset of vertices offering the
highest total summed weights whilst ensuring that no pair of adjacent vertices
(i.e. connected by an edge) are included. Introducing the binary variables
$n_{i}=\\{0,1\\}$ to indicate whether vertex $i$ belongs to the MWIS set, the
classical cost function that must be minimized to solve the MWIS problem is

$H_{\mathrm{MWIS}}=-\sum_{i\in V}w_{j}n_{j}+\sum_{(i,j)\in
E}U_{ij}n_{i}n_{j},$ (1)

where $U_{ij}>\mathrm{max}_{k}(w_{k})$. Here the first term rewards inclusion
of high-weight vertices, whilst the second term acts to penalize inclusion of
vertices connected by an edge. For the case that $G$ is a UDG, this cost
function can be natively mapped to the case of a neutral atom array which has
the Hamiltonian given by

$H_{\mathrm{Ryd}}=\sum_{i}(\frac{\Omega}{2}\sigma_{i}^{x}-\Delta_{i}\hat{n}_{i})+\sum_{j>i}V(|\bm{r}_{i}-\bm{r}_{j}|)\hat{n}_{i}\hat{n}_{j},$
(2)

where $\Omega$ is the Rabi frequency,
$\sigma_{i}^{x}=\ket{r}_{i}\bra{g}+\ket{g}_{i}\bra{r}$, $\Delta_{i}$ is a site
dependent detuning, $\hat{n}_{i}=\ket{r}_{i}\langle r|$ is the projector onto
the Rydberg state of atom $i$, $\bm{r}_{i}$ is a position vector and
$V(r)=C_{6}/r^{6}$ is the dipole-dipole interaction strength between two
atoms. In the limit $\Omega=0$, the ground state of $H_{\mathrm{Ryd}}$ encodes
the solution to the MWIS problem if $\Delta_{i}=w_{i}$ and
$V(|\bm{r}_{i}-\bm{r}_{j}|)>\mathrm{max}_{k}(w_{k})$ for all $i,j$ sharing an
edge, enabling solution of UDG-MWIS problems by annealing atoms to the ground
state configuration [44].

For arbitrary problems, such as non-UDG graphs or QUBO, the graph (and
associated constraints) must first be encoded into a weighted UDG-MWIS. Ref.
[44] provides one approach to implementing this mapping by introducing regular
sub-graph gadgets to realize arbitrary non-local connections between vertices,
however the focus of this present work is to verify the experimental
feasibility for preparing the MWIS on a neutral atom array.

Figure 2: 1D Weighted Graphs. (a) For a uniformly weighted, odd-length 1D
graph the ground state is the $Z_{2}$-ordered phase with Rydberg excitations
on odd sites which corresponds to the unweighted MIS. Introducing a weighting
with $w_{i}=2$ on even sites results in an MWIS ground state with Rydberg
excitations localized to the even sites which is no longer equivalent to the
MIS solution. (b) Optimal annealing ramp for preparing the weighted ground
state obtained via closed-loop optimization for $N=9$ atoms spaced by
$a=7~{}\mu$m. Output state probability (c) and time evolution (d) for the
unweighted graph showing the odd-ordered target ground state is prepared with
$19(1)\%$ probability. Output state probability (e) and time evolution (f) for
the weighted graph showing even-ordered ground state is also prepared with
$19(1)\%$ probability.

_MWIS Annealing Protocol_.— Starting with a target MWIS graph problem, the
graph is first mapped to a UDG-MWIS which can be geometrically embedded onto a
neutral atom array using techniques from [44] with a spacing $a$ between
neighbouring atoms chosen to be below the blockade radius $r_{\mathrm{B}}$ to
ensure the blockade constraint is realized. Graph weighting is implemented
with an additional light-shifting potential that is holographically projected
onto the atoms using a spatial light modulator (SLM) to apply local light-
shifts with arbitrary geometry and programmable relative weighting. These
additional tweezer beams, shown schematically in Fig. 1, cause a differential
AC Stark shift $\delta_{\mathrm{AC}}^{i}$ on each atom proportional to the
relative power in each tweezer, resulting in a total shift
$\Delta_{i}=\Delta+\delta_{\mathrm{AC}}^{i}$, where $\Delta$ is the detuning
of the global Rydberg lasers used to realize a homogeneous two-photon coupling
between ground and Rydberg states. In the experiments below, the light-shift
laser is blue-detuned for both ground and Rydberg states, ensuring that
$\delta_{\mathrm{AC}}^{i}\geq 0$ as required. For a given graph weighting, the
relative light-shifts are calibrated to realize
$\delta_{\mathrm{AC}}^{i}=w_{i}\delta_{\mathrm{AC}}$, where
$\delta_{\mathrm{AC}}$ is the global unit scale ($w=1$) of the light-shift and
$w_{i}$ are the relative weights of the embedded graph problem. To avoid
violating the blockade condition, it is necessary to constrain the maximum
value of $\delta_{\mathrm{AC}}$ such that
$\mathrm{max}_{i}(\delta_{\mathrm{AC}}^{i})<V(a)$, where $a$ is the
interatomic spacing.

To adiabatically prepare atoms in the ground state encoding the solution of
the UDG-MWIS problem, a cubic detuning profile is used [50], as illustrated in
Fig. 1(c). Starting with the atoms initially prepared in $\ket{g}$, the global
Rabi frequency $\Omega$ is ramped up using a fixed negative detuning after
which the detuning is swept towards resonance. Once the global detuning
reaches resonance, $\Delta=0$, the global detuning is held fixed and the
remaining detuning profile is applied by dynamically varying the intensity of
the light-shifting potential to control $\delta_{\mathrm{AC}}(t)$. By deriving
the local light-shifts from a single laser, it is possible to use global
intensity control for implementing arbitrary detuning ramps whilst maintaining
fixed relative weightings on the graph. The Rabi frequency $\Omega$ is then
ramped off whilst maintaining the light-shift, after which the encoded
bitstrings $\\{n_{i}\\}$ are read out by fluorescence imaging. Atoms excited
to the Rydberg state are ejected from the trap, whilst atoms in the ground
state survive and appear in the images. Using a closed-loop optimization
process, the cubic ramp parameters are adjusted to minimize the classical cost
function $\langle H_{\mathrm{MWIS}}\rangle$ to maximize the probability of
preparing the system in the MWIS ground state, as shown in Fig. 1(d).
Following this optimization process, the solution of the target MWIS graph
problem can then be determined from the most probable output bitstrings.

_Experimental Setup_.— The experimental setup is based on arrays of
holographically trapped Cs atoms. Atoms are initially prepared in the
$|g\rangle=|4,0\rangle$ state, and are globally excited to the
$\ket{r}=\ket{80s_{1/2},m_{j}=+1/2}$ Rydberg state via two-photon excitation
using the $7p_{1/2}$ intermediate state using light at 459 nm and 1039 nm. We
operate with a detuning of $\Delta^{\prime}/2\pi=$500 MHz from the
intermediate state, and achieve a Rydberg Rabi frequency of
$\Omega/2\pi=2.70(2)$ MHz. Global detuning is controlled by adjusting the
frequency using a double-passed acousto-optic modulator (AOM) in the 459 nm
path. For $80s_{1/2}$ we calculate $C_{6}=-3376~{}\mathrm{GHz}\,\mu$m6 [51],
resulting in a resonant blockade radius
$r_{B}=\sqrt[6]{|C_{6}|/\Omega}=10.4~{}\mu$m.

Local light-shifts are realized using light at 800 nm derived from a Ti:Sapph
laser that is projected from a spatial light modulator (SLM) and overlapped
with the main optical trapping array. To suppress mechanical forces coming
from the repulsive optical potential, the light-shift spot arrays are re-
imaged onto the atoms with a $1/e^{2}$ radius of $5.1(2)~{}\mu$m. For a given
UDG-MWIS, the arrays of local light-shift beams are generated using a
Gerchberg-Saxon algorithm with the desired relative weighting [52]. To verify
the light-shifts experienced by the atoms, we perform spectroscopy on the
$\ket{g}\rightarrow\ket{50s_{1/2},m_{j}=+1/2}$ transition to suppress
interactions, and perform a closed-loop feedback on the SLM potentials to
achieve target weightings with $<2\%$ RMS errors in the relative weightings
after 5 iterations [50]. The resulting differential light-shift is calibrated
as $\delta_{\mathrm{AC}}=0.28(2)$ MHz$/$mW, which is dominated by the ground
state light-shift as the Rydberg AC shift calculated using the polarizability
of a free electron corresponds to only 40 kHz [53].

Figure 3: 2D Weighted Graphs (a) Target graph is a 5 vertex non-UDG MWIS,
which can be mapped to the 9 atom UDG-MWIS graph [44] by introducing four
ancilla qubits with weightings $w_{\alpha}$ and $w_{\beta}$ to implement the
edges from vertex 5 to vertices 1 and 3. (b) Optimized annealing profile for
the graph instance with $(w_{1},w_{2},w_{3},w_{4},w_{5})=(2,1,2,1,1)$ which
results in $\langle H_{\mathrm{MWIS}}\rangle=-4.0(1)$, compared to the
idealized cost of $-8$. (c) Time-evolution data for the weighted graph showing
adiabatic evolution during the ramp to reach the correct target state. (d)
Final output state probability distribution from 1000 experimental
realizations showing the correct MWIS solution is observed with $11(1)\%$
probability, with all other output strings observed with $\leq 2.5\%$
probability.

_1D Weighted Graphs_.— As an artificial first demonstration of the feasibility
to perform weighted optimization using the dual-stage protocol defined above,
we begin with a simple and well studied problem of a 1D weighted graph as
illustrated in Fig. 2(a). For a uniformly weighted set of $N$ atoms with
$w_{i}=1$ spaced such that $a<r_{\mathrm{B}}<2a$, for odd $N$ the ground state
corresponds to a $Z_{2}$ ordered phase with Rydberg excitations localized on
the odd sites as previously observed using global control only [28].
Introducing a relative weighting of $w_{i}=2$ on the even sites results in a
modification of the MWIS ground state to prepare the inverted output string
with Rydberg excitations now localized to the even sites.

We perform experiments on weighted 1D graphs using two 20 $\mu$m separated
parallel lines of $N=9$ atoms spaced by $a=7~{}\mu$m, with one row set to have
uniform light-shifts and the second row having $w_{i}=2$ on the even sites.
Starting with a linear profile, we apply the global light-shift on both graphs
simultaneously, and perform closed-loop feedback to minimize $\langle
H_{\mathrm{MWIS}}\rangle$ for the weighted graph using the classical cost
function defined in Eq. 1. Whilst encoding the same effective ground state
energy as $H_{\mathrm{Ryd}}$, using the classical cost function avoids issues
with the optimizer trivially driving the system to large negative detunings
$\Delta_{i}$ as a route to further minimising the energy of the observed
output configurations. The resulting optimized profile is shown in Fig. 2(b),
corresponding to a duration of 3 $\mu$s with light-shifts turned on from 1.5
$\mu$s when the global Rydberg laser detuning reaches zero.

The resulting distributions of output strings are shown for both uniform and
weighted graphs in Fig. 2(c,e) respectively. Using a common annealing profile,
we observe the correct ground states as the most probable state outputs in
both cases, with probabilities of $19(1)\%)$ in both cases. The additional 9%
peak in Fig. 2(e) with state index 93 corresponds to the target output state
with the first qubit also lost. This appears higher than other states due to
slightly increased single atom loss measured for this array site. Time
evolution data is shown for each graph in Fig. 2(d,f) respectively, which
clearly demonstrates the adiabatic nature of graph preparation at longer times
and the emergence of the ground state ordering as expected.

_2D Weighted Graphs_.— Moving beyond the trivial 1D case, we consider the
example of a non-UDG 2D MWIS shown in Fig. 3(a). This graph, featuring five
vertices, has MWIS solutions corresponding to either qubits 1 and 3 being
excited, or qubits 2, 4 and 5. To map this to a UDG-MWIS graph for
implementation on our neutral atom array, we follow the protocol of Ref. [44]
resulting in a 9-qubit graph requiring four ancilla qubits to implement the
edges between vertices 1 and 2 and vertex 5. By symmetry this introduces two
additional weightings, $w_{\alpha}$ and $w_{\beta}$ on the ancilla qubits. We
use the weightings $w_{\alpha}=(w_{1}+w_{3})/2$ and
$w_{\beta}=(w_{2}+w_{4}+w_{5})/2$ which guarantees that for any arbitrary
weighting with $w_{i}>0$ we realize the MWIS solution as the ground state of
the Rydberg interaction Hamiltonian.

Initially we consider the weighted graph instance indicated by the inset of
Fig. 3(d) with $(w_{1},w_{2},w_{3},w_{4},w_{5})=(2,1,2,1,1)$ using a $3\times
3$ arrangement of atoms with spacing $a=8~{}\mu$m. As before, we perform
closed-loop optimization of the annealing profile, resulting in a shorter
annealing time of only 1.65 $\mu$s with the light-shift turning on after 0.6
$\mu$s, as shown in Fig. 3(b). For this optimization process we minimize
$\langle H_{\mathrm{MWIS}}\rangle$ for the embedded 9 qubit graph (rather than
simply just considering the non-embedded 5-qubit graph), as whilst the values
of the ancilla qubits are not relevant to the problem solution, there are many
possible states encoded in the 9-qubit bitstrings that feature the correct
5-vertex solution state but only one of which corresponds to the energetic
ground state. This ensures the optimizer favours solutions that more closely
follow adiabatic preparation of the ground states, which can be verified by
observing the dynamical evolution in preparing the correct output bitstring.
Fig. 3(c) shows the preparation probability as a function of time, revealing a
smooth accumulation of population after 1.2 $\mu$s and no fast oscillations
observed when exploiting shorter annealing profile durations. Following
optimization of ramp parameters we obtain $\langle
H_{\mathrm{MWIS}}\rangle=-4.0(1)$, compared to an idealized cost of -8.

The final output state probabilities for the weighted graph are shown in Fig.
3(d), with the correct 9-qubit output state corresponding to Rydberg
excitation of qubits mapping onto vertices 1 and 3 being observed with $\geq
10\%$ probability, significantly higher than any other output bitstrings and
showing demonstration of solving an embedded MWIS graph.

One further advantage of identifying an adiabatic annealing profile is that it
should be robust against changes in graph weightings. To verify this for the
graph shown in Fig. 3, we directly implement the same annealing profile but
with different relative weightings of $(1,2,1,2,2)$ and $(7/4,5/4,9/4,3/2,1)$.
The results are shown in Fig. 4 where for both cases the annealing process is
able to prepare the system in the correct ground state, even for Fig. 4(a)
which corresponds to the inverted solution state to that targeted previously
to optimize the original annealing profile.

Figure 4: Annealing alternative graph instances. We use the optimised ramp
found on the $(2,1,2,1,1)$ graph and apply it to (a) $(1,2,1,2,2)$ and (b)
$(7/4,5/4,9/4,3/2,1)$ graph instances. In both cases, the ramp provides the
correct solution to the MWIS graph even for the case of the inverted ground
state for (a), resulting in target state output probabilities of $17(2)\%$ and
$10(1)\%$ with $\langle H_{\mathrm{MWIS}}\rangle=-6.9(2)$ and $-3.5(1)$
respectively.

_Discussion and Outlook_ In this paper we have demonstrated a new approach to
implementing UDG-MWIS using arbitrary local light-shifts implemented using a
SLM combined with a dual stage annealing profile that combines the advantages
of local light-shifts with global frequency and intensity control for
implementing a quantum annealing algorithm. We demonstrate the ability to
prepare ground states of small-scale 1D and 2D graphs, including demonstrating
the embedding of a non-UDG MWIS problem. In both cases, we were able to
identify robust annealing protocols suitable for use with different relative
weightings providing the possibility to design universal ramps that are
suitable for solving a large number of graphs without the required overhead of
full closed-loop optimization of the annealing ramp each time.

This approach for local light-shifts is highly scalable, with SLMs able to
generate arrays of $>\\!\\!1000$ individually controllable spots and the
option to implement in-situ realtime feedback using a camera calibrated
against the measured atomic shift previously demonstrated for normalization of
arbitrary trapping potentials [25]. Currently at 800 nm we have up to 4 W
available at the experiment, providing sufficient power for $>\\!\\!100$
weighted sites. In future, this can be scaled further using coherent
combination or frequency doubling telecoms fiber lasers to increase the
available power and combined with techniques for beam-shaping to implement a
flat-top point spread function on the SLM [54] to reduce the size of the
light-shift spots and further reduce mechanical effects caused by using a
repulsive optical potential to enable longer annealing runs.

For larger graph problems, the quadratic overhead introduced by the UDG-MWIS
encoding [44] introduces interesting open questions around scalability of this
approach, whilst recent theoretical work exploring annealing of unweighted
graph shows even for relatively simple graph constructions it is possible to
have minimum energy gaps that decay superexponentially with system size [55]
requiring use of optimal control techniques or quenches to prepare atoms in
the target ground states. In future work we will explore this question
further, and provide a comparison against other encoding protocols requiring
local light-shifts as a route to identifying useful applications of weighted
graph optimization on neutral atom arrays [45, 46].

###### Acknowledgements.

We thank H. Pichler for helpful discussions. This work is supported by the
EPSRC Prosperity Partnership _SQuAre_ (Grant No. EP/T005386/1) with funding
from M Squared Lasers Ltd. The data presented in this work are available at
[56].

## References

  * Abbas  _et al._ [2023] A. Abbas, A. Ambainis, B. Augustino, A. Bärtschi, H. Buhrman, C. Coffrin, G. Cortiana, V. Dunjko, D. J. Egger, B. G. Elmegreen, N. Franco, F. Fratini, B. Fuller, J. Gacon, C. Gonciulea, S. Gribling, S. Gupta, S. Hadfield, R. Heese, G. Kircher, T. Kleinert, T. Koch, G. Korpas, S. Lenk, J. Marecek, V. Markov, G. Mazzola, S. Mensa, N. Mohseni, G. Nannicini, C. O’Meara, E. P. Tapia, S. Pokutta, M. Proissl, P. Rebentrost, E. Sahin, B. C. B. Symons, S. Tornow, V. Valls, S. Woerner, M. L. Wolf-Bauwens, J. Yard, S. Yarkoni, D. Zechiel, S. Zhuk, and C. Zoufal, arXiv:2312.02279 [quant-ph] (2023).
  * Farhi  _et al._ [2000] E. Farhi, J. Goldstone, S. Gutmann, and M. Sipser, arXiv:quant-ph/0001106 [quant-ph] (2000).
  * Farhi  _et al._ [2001] E. Farhi, J. Goldstone, S. Gutmann, J. Lapan, A. Lundgren, and D. Preda, Science 292, 472 (2001).
  * Santoro and Tosatti [2006] G. E. Santoro and E. Tosatti, Journal of Physics A: Mathematical and General 39, R393 (2006).
  * Albash and Lidar [2018] T. Albash and D. A. Lidar, Rev. Mod. Phys. 90, 015002 (2018).
  * Farhi  _et al._ [2014] E. Farhi, J. Goldstone, and S. Gutmann, 1411.4028 [quant-ph] (2014).
  * Zhou  _et al._ [2020] L. Zhou, S.-T. Wang, S. Choi, H. Pichler, and M. D. Lukin, Phys. Rev. X 10, 021067 (2020).
  * Guerreschi and Matsuura [2019] G. G. Guerreschi and A. Y. Matsuura, Scientific Reports 9, 6903 (2019).
  * Lykov  _et al._ [2023] D. Lykov, J. Wurtz, C. Poole, M. Saffman, T. Noel, and Y. Alexeev, npj Quantum Information 9, 73 (2023).
  * Shaydulin  _et al._ [2023] R. Shaydulin, C. Li, S. Chakrabarti, M. DeCross, D. Herman, N. Kumar, J. Larson, D. Lykov, P. Minssen, Y. Sun, Y. Alexeev, J. M. Dreiling, J. P. Gaebler, T. M. Gatterman, J. A. Gerber, K. Gilmore, D. Gresh, N. Hewitt, C. V. Horst, S. Hu, J. Johansen, M. Matheny, T. Mengle, M. Mills, S. A. Moses, B. Neyenhuis, P. Siegfried, R. Yalovetzky, and M. Pistoia, arXiv:2308.02342 [quant-ph] (2023).
  * Bauza and Lidar [2024] H. M. Bauza and D. A. Lidar, arXiv:2401.07184 [quant-ph] (2024).
  * Saffman [2016] M. Saffman, J. Phys. B 49, 202001 (2016).
  * Henriet  _et al._ [2020] L. Henriet, L. Beguin, A. Signoles, T. Lahaye, A. Browaeys, G.-O. Reymond, and C. Jurczak, Quantum 4, 327 (2020).
  * Wu  _et al._ [2021] J.-L. Wu, Y. Wang, J.-X. Han, S.-L. Su, Y. Xia, Y. Jiang, and J. Song, Phys. Rev. A 103, 012601 (2021).
  * Morgado and Whitlock [2021] M. Morgado and S. Whitlock, AVS Quantum Science 3, 023501 (2021).
  * Huft  _et al._ [2022] P. Huft, Y. Song, T. M. Graham, K. Jooya, S. Deshpande, C. Fang, M. Kats, and M. Saffman, Phys. Rev. A 105, 063111 (2022).
  * Norcia  _et al._ [2024] M. A. Norcia, H. Kim, W. B. Cairncross, M. Stone, A. Ryou, M. Jaffe, M. O. Brown, K. Barnes, P. Battaglino, T. C. Bohdanowicz, A. Brown, K. Cassella, C. A. Chen, R. Coxe, D. Crow, J. Epstein, C. Griger, E. Halperin, F. Hummel, A. M. W. Jones, J. M. Kindem, J. King, K. Kotru, J. Lauigan, M. Li, M. Lu, E. Megidish, J. Marjanovic, M. McDonald, T. Mittiga, J. A. Muniz, S. Narayanaswami, C. Nishiguchi, T. Paule, K. A. Pawlak, L. S. Peng, K. L. Pudenz, D. R. Perez, A. Smull, D. Stack, M. Urbanek, R. J. M. van de Veerdonk, Z. Vendeiro, L. Wadleigh, T. Wilkason, T. Y. Wu, X. Xie, E. Zalys-Geller, X. Zhang, and B. J. Bloom, arXiv:2401.16177 [quant-ph] (2024).
  * Gyger  _et al._ [2024] F. Gyger, M. Ammenwerth, R. Tao, H. Timme, S. Snigirev, I. Bloch, and J. Zeiher, arXiv:2402.04994 [quant-ph] (2024).
  * Manetsch  _et al._ [2024] H. J. Manetsch, G. Nomura, E. Bataille, K. H. Leung, X. Lv, and M. Endres, arXiv:2403.12021 [quant-ph] (2024).
  * Nikolov  _et al._ [2023] B. Nikolov, E. Diamond-Hitchcock, J. Bass, N. L. R. Spong, and J. D. Pritchard, Phys. Rev. Lett. 131, 030602 (2023).
  * Levine  _et al._ [2019] H. Levine, A. Keesling, G. Semeghini, A. Omran, T. T. Wang, S. Ebadi, H. Bernien, M. Greiner, V. Vuletić, H. Pichler, and M. D. Lukin, Phys. Rev. Lett. 123, 170503 (2019).
  * Evered  _et al._ [2023] S. J. Evered, D. Bluvstein, M. Kalinowski, S. Ebadi, T. Manovitz, H. Zhou, S. H. Li, A. A. Geim, T. T. Wang, N. Maskara, H. Levine, G. Semeghini, M. Greiner, V. Vuletić, and M. D. Lukin, Nature 622, 268 (2023).
  * Ma  _et al._ [2023] S. Ma, G. Liu, P. Peng, B. Zhang, S. Jandura, J. Claes, A. P. Burgers, G. Pupillo, S. Puri, and J. D. Thompson, Nature 622, 279 (2023).
  * Graham  _et al._ [2022] T. M. Graham, Y. Song, J. Scott, C. Poole, L. Phuttitarn, K. Jooya, P. Eichler, X. Jiang, A. Marra, B. Grinkemeyer, , M. Kwon, , M. Ebert, J. Cherek, M. T. Lichtman, M. Gillette, J. Gilbert, D. Bowman, T. Ballance, C. Campbell, E. D. Dahl, O. Crawford, N. S. Blunt, B. Rogers, T. Noel, , and M. Saffman, Nature 604, 457 (2022).
  * Bluvstein  _et al._ [2024] D. Bluvstein, S. J. Evered, A. A. Geim, S. H. Li, H. Zhou, T. Manovitz, S. Ebadi, M. Cain, M. Kalinowski, D. Hangleiter, J. P. Bonilla Ataides, N. Maskara, I. Cong, X. Gao, P. Sales Rodriguez, T. Karolyshyn, G. Semeghini, M. J. Gullans, M. Greiner, V. Vuletić, and M. D. Lukin, Nature 626, 58 (2024).
  * Bluvstein  _et al._ [2022] D. Bluvstein, H. Levine, G. Semeghini, T. T. Wang, S. Ebadi, M. Kalinowski, A. Keesling, N. Maskara, H. Pichler, M. Greiner, V. Vuletic̀, and M. D. Lukin, Nature 604, 451 (2022).
  * Browaeys and Lahaye [2020] A. Browaeys and T. Lahaye, Nature Phys. 16, 132 (2020).
  * Bernien  _et al._ [2017] H. Bernien, S. Schwartz, A. Keesling, H. Levine, A. Omran, H. Pichler, S. Choi, A. S. Zibrov, M. Endres, M. Greiner, V. Vuletić, and M. D. Lukin, Nature 551, 579 (2017).
  * Ebadi  _et al._ [2021] S. Ebadi, T. T. Wang, H. Levine, A. Keesling, G. Semeghini, A. Omran, D. Bluvstein, R. Samajdar, H. Pichler, W. W. Ho, S. Choi, S. Sachdev, M. Greiner, V. Vuletić, , and M. D. Lukin, Nature 595, 227 (2021).
  * Scholl  _et al._ [2021] P. Scholl, M. Schuler, H. J. Williams, A. A. Eberharter, D. Barredo, K.-N. Schymik, V. Lienhard, L.-P. Henry, T. C. Lang, T. Lahaye, A. M. Läuchli, and A. Browaeys, Nature 595, 233 (2021).
  * de Léséleuc  _et al._ [2019] S. de Léséleuc, V. Lienhard, P. Scholl, D. Barredo, S. Weber, N. Lang, H. Buchler, T. Lahaye, and A. Browaeys, Science 365, 775 (2019).
  * Semeghini  _et al._ [2021] G. Semeghini, H. Levine, A. Keesling, S. Ebadi, T. T. Wang, D. Bluvstein, R. Verresen, H. Pichler, M. Kalinowski, R. Samajdar, A. Omran, S. Sachdev, A. Vishwanath, M. Greiner, V. Vuletic, and M. D. Lukin, Science 374, 1242 (2021).
  * Kim  _et al._ [2023a] M. Kim, J. Ahn, Y. Song, J. Moon, and H. Jeong, Journal of the Korean Physical Society 82, 827 (2023a).
  * Pichler  _et al._ [2018a] H. Pichler, S.-T. Wang, L. Zhou, S. Choi, and M. D. Lukin, 1808.10816 [quant-ph] (2018a).
  * Pichler  _et al._ [2018b] H. Pichler, S.-T. Wang, L. Zhou, S. Choi, and M. D. Lukin, 1809.04954 [quant-ph] (2018b).
  * Jaksch  _et al._ [2000] D. Jaksch, J. I. Cirac, P. Zoller, S. L. Rolston, R. Côté, and M. D. Lukin, Phys. Rev. Lett. 85, 2208 (2000).
  * Ebadi  _et al._ [2022] S. Ebadi, A. Keesling, M. Cain, T. T. Wang, H. Levine, D. Bluvstein, G. Semeghini, A. Omran, J.-G. Liu, R. Samajdar, X.-Z. Luo, B. Nash, X. Gao, B. Barak, E. Farhi, S. Sachdev, N. Gemelke, L. Zhou, S. Choi, H. Pichler, S.-T. Wang, M. Greiner, V. Vuletic̀, and M. D. Lukin, Science 376, 1209 (2022).
  * Kim  _et al._ [2023b] K. Kim, M. Kim, J. Park, A. Byun, and J. Ahn, arXiv:2311.13803 [quant-ph] (2023b).
  * Andrist  _et al._ [2023] R. S. Andrist, M. J. A. Schuetz, P. Minssen, R. Yalovetzky, S. Chakrabarti, D. Herman, N. Kumar, G. Salton, R. Shaydulin, Y. Sun, M. Pistoia, and H. G. Katzgraber, Phys. Rev. Res. 5, 043277 (2023).
  * Kim  _et al._ [2022] M. Kim, K. Kim, J. Hwang, E.-G. Moon, and J. Ahn, Nature Phys. 18, 755 (2022).
  * Dalyac  _et al._ [2023] C. Dalyac, L.-P. Henry, M. Kim, J. Ahn, and L. Henriet, arXiv:2306.13373 [quant-ph] (2023).
  * Park  _et al._ [2024] J. Park, S. Jeong, M. Kim, K. Kim, A. Byun, L. Vignoli, L.-P. Henry, L. Henriet, and J. Ahn, arXiv:2312.08703 [quant-ph] (2024).
  * Byun  _et al._ [2023] A. Byun, J. Jung, K. Kim, M. Kim, S. Jeong, H. Jeong, and J. Ahn, Rydberg-atom graphs for quadratic unconstrained binary optimization problems (2023), arXiv:2309.14847 [quant-ph] .
  * Nguyen  _et al._ [2023] M.-T. Nguyen, J.-G. Liu, J. Wurtz, M. D. Lukin, S.-T. Wang, and H. Pichler, PRX Quantum 4, 010316 (2023).
  * Lanthaler  _et al._ [2023] M. Lanthaler, C. Dlaska, K. Ender, and W. Lechner, Phys. Rev. Lett. 130, 220601 (2023).
  * Goswami  _et al._ [2023] K. Goswami, R. Mukherjee, H. Ott, and P. Schmelcher, arXiv:2308.07798 [quant-ph] (2023).
  * Omran  _et al._ [2019] A. Omran, H. Levine, A. Keesling, G. Semeghini, T. T. Wang, S. Ebadi, H. Bernien, A. S. Zibrov, H. Pichler, S. Choi, J. Cui, M. Rossignolo, P. Rembold, S. Montangero, T. Calarco, M. Endres, M. Greiner, V. Vuletić, and M. D. Lukin, Science 365, 570 (2019).
  * Bornet  _et al._ [2024] G. Bornet, G. Emperauger, C. Chen, F. Machado, S. Chern, L. Leclerc, B. Gély, D. Barredo, T. Lahaye, N. Y. Yao, and A. Browaeys, arXiv:2402.11056 [quant-ph] (2024).
  * Wurtz  _et al._ [2024] J. Wurtz, P. L. S. Lopes, C. Gorgulla, N. Gemelke, A. Keesling, and S. Wang, arXiv:2205.08500 [quant-ph] (2024).
  * [50] See Supplemental Material.
  * S̆ibalić  _et al._ [2017] N. S̆ibalić, J. D. Pritchard, C. S. Adams, and K. J. Weatherill, Comp. Phys. Comm. 220, 319 (2017).
  * Kim  _et al._ [2019] D. Kim, A. Keesling, A. Omran, H. Levine, H. Bernien, M. Greiner, M. D. Lukin, and D. R. Englund, Opt. Lett. 44, 3178 (2019).
  * Zhang  _et al._ [2011] S. Zhang, F. Robicheaux, and M. Saffman, Phys. Rev. A 84, 043408 (2011).
  * Schroff  _et al._ [2023] P. Schroff, A. La Rooij, E. Haller, and S. Kuhr, Scientific Reports 13, 3252 (2023).
  * Schiffer  _et al._ [2024] B. F. Schiffer, D. S. Wild, N. Maskara, M. Cain, M. D. Lukin, and R. Samajdar, Phys. Rev. Res. 6, 013271 (2024).
  * [56] DOI to be added in proof.
  * Sortais  _et al._ [2007] Y. R. P. Sortais, H. Marion, C. Tuchendler, A. M. Lance, M. Lamare, P. Fournet, C. Armellin, R. Mercier, G. Messin, A. Browaeys, and P. Grangier, Phys. Rev. A 75, 013406 (2007).
  * de Léséleuc  _et al._ [2018] S. de Léséleuc, D. Barredo, V. Lienhard, A. Browaeys, and T. Lahaye, Phys. Rev. A 97, 053803 (2018).
  * Legaie  _et al._ [2018] R. Legaie, C. J. Picken, and J. D. Pritchard, J. Opt. Soc. Am. B 35, 892 (2018).
  * Preuschoff  _et al._ [2020] T. Preuschoff, M. Schlosser, and G. Birkl, Rev. Sci. Inst. 91, 083001 (2020).
  * Wigley  _et al._ [2016] P. B. Wigley, P. J. Everitt, A. van den Hengel, J. W. Bastian, M. A. Sooriyabandara, G. D. McDonald, K. S. Hardman, C. D. Quinlivan, P. Manju, C. C. N. Kuhn, I. R. Petersen, A. N. Luiten, J. J. Hope, N. P. Robins, and M. R. Hush, Scientific Reports 6, 25890 (2016).

Supplemental Material for Demonstration of weighted graph optimization on a
Rydberg atom array using local light-shifts

Figure 5: Experimental Setup. (a) Atom arrays are generated using 1064 nm
trapping light projected into the chamber using an SLM and a high-numerical
aperture in-vacuum aspheric lens. Local light-shifts are implemented using a
second SLM to overlap spots at 800 nm with the 1064 nm tweezers. Atom imaging
is performed by collecting fluoresence onto an sCMOS camera. Rydberg
excitation is performed using counter-propagating beams at 459 nm and 1039 nm
with orthogonal circular polarisiations which globally illuminate the array.
(b) Two-photon Rydberg excitation scheme from the $6S_{1/2}$ clock state
$|g\rangle=|4,0\rangle$ to $|r\rangle=|ns_{1/2},m_{j}=+1/2\rangle$ via the
intermediate $7P_{1/2}$ level. Lasers are locked with a detuning of
$\Delta^{\prime}/2\pi=500$ MHz from the $7P_{1/2}$ transition.

##  I Experiment Setup

###  I.1 Neutral Atom Arrays

A schematic of the experimental setup is shown in Fig. 5(a), based on the
setup previously introduced in Ref. [20]. Atom arrays are created by loading
individual Cs atoms into 1064 nm optical tweezers generated using a spatial
light modulator (SLM). Traps are focused to a $1/e^{2}$ waist of 1.5 $\mu$m
using an in-vacuum aspheric lens with NA=0.45 [57] which has an ITO coating to
suppress charge build up. This same lens is used for qubit readout by
performing fluorescence imaging with an sCMOS (Teledyne Photometrics Prime
BSI) camera.

Following stochastic loading of the array, we rearrange the atoms into the
target spatial configuration with a mobile tweezer beam at 1039 nm created
using a dual-axis acousto-optic deflector (AOD) driven by a fast micro-
controller. Atom sorting is performed using a trap depth of 1 mK and a tweezer
of 9.8 mK, using a customised algorithm to calculate moves in real time that
aims to minimise the number of atoms movements required to fill the target
array sites. We achieve single move efficiencies $\geq 97\%$, and typical
success rates $\geq 70\%$ for the graphs shown above using a single round of
sorting. For larger graphs, multiple sorting iterations can be performed.

After rearrangement, atoms are optically pumped into the $\ket{4,0}$ hyperfine
state using $\pi$-polarized light on the Cs $D_{1}$ transition resonant with
the $F=4\rightarrow F^{\prime}=4$, achieving an optical pumping efficiency of
$0.971(9)$ measured using microwave benchmarking [20].

For the experiments in this paper, prior to Rydberg operations the optical
dipole traps are ramped down to 15 $\mu$K, cooling the atoms to temperatures
of 2 $\mu$K. Atoms are then released for 5 $\mu$s to minimise the mechanical
force of the tweezer beams on the Rydberg states, and atoms are then
recaptured at 1 mK to facilitate efficient Rydberg state detection [58].

Figure 6: Feedback on local light-shift potentials. (a) Spectroscopy of
initial hologram generation for the $3x3$ atom arrangement with
$(w_{1},w_{2},w_{3},w_{4},w_{5})=(1,2,1,2,2)$ requiring $w_{\alpha}=2$ and
$w_{\beta}=3$ shows significant variation in the shift experienced by each
atom. (b) Performing closed-loop feedback on the SLM hologram generation using
the measured light shifts enables rapid convergence onto the target light
shifts with RMS below $<2\%$. (c) Spectroscopy measurement performed after
iteration 5 showing close convergence of the relative shifts experienced by
each atom. The dashed line shows the targeted resonance for each atom.

###  I.2 Rydberg Excitation

Rydberg excitation is performed using two-photon excitation via the
intermediate excited state as shown in Fig. 5(b) using light at 459 nm and
1039 nm to drive transitions from $6P_{1/2}\rightarrow 7P_{1/2}$ and
$7P_{1/2}\rightarrow nS_{1/2}$ respectively. Both beams are derived from
Ti:Sapph lasers (M Squared SolsTiS) at 918 nm and 1039 nm which are frequency
stabilised to a commercial ULE cavity (M Squared) offering finesse
$\mathcal{F}>40$ k at both wavelengths. The lasers are offset-locked using a
Pound-Drever-Hall with fiber EOMs [59] to provide controllable detuning with
respect to the 3 GHz spaced cavity modes, with few kHz linewidths measured via
a beatnote between the two lasers.

The 459 nm light is frequency doubled (M Squared ECDF) and the detuning
controlled using a double pass accousto-optic modulator (AOM,AA
MQ180-A0.25-VIS) in a cat-eye geometry to enable fast frequency control. The
AOM is driven using a high bandwidth arbitrary waveform generator (AWG,
Spectrum M4i6631-x8). This is delivered to the experiment by fiber, and
overlapped onto the array with a 125 $\mu$m $1/e^{2}$ waist. The 1039 nm laser
seeds a fiber amplifier (Azur ALS-1039-20-CP) which is intensity controlled
using a single-pass AOM and coupled into a high-power fiber. The beam is then
imaged onto the atoms using cylindrical optics to create a beam with a
$1/e^{2}$ waist size of $60\times 20~{}\mu$m. For the experiments presented
here, Rydberg excitation is performed using an intermediate detuning of
$\Delta^{\prime}/2\pi=500$ MHz with respect to the $7P_{1/2}$ centre of mass
energy, with powers set to provide a Rabi frequency of $\Omega/2\pi=2.70(0)$
MHz over the central rows of the array.

To stabilise the intensity, prior to each experimental cycle a 1 ms sample
pulse is applied allowing closed-loop feedback to an AOM based noise eater on
each laser [60], with the intensity held constant for the fast Rydberg pulses.
We achieve $1\%$ shot-to-shot stability stable over hours using this
technique.

###  I.3 Local Light-shifts

To apply local light-shifts for implementing graph weighting, an 800 nm laser
(M Squared SolsTiS) is used with a second spatial light modulator to create a
secondary tweezer array that is overlapped with the underlying 1064 nm traps
as shown in Fig. 5(a). Due to the chromatic shifts of the in-vacuum aspheric,
it is not possible to image both 1064 nm and 800 nm arrays overlapped after
the chamber focused in a common plane. Instead, to overlap the arrays we
exploit the fact that ground-state atoms are anti-trapped by this blue detuned
light. Using this effect, we initially generate arrays with a 2 $\mu$m waist
and apply the 800 nm potential onto the trapped atoms for 1-10 ms, resulting
in ejection of atoms from sites that are overlapped between the two beams.
Using $10\times 10$ arrays with 8 $\mu$m spacing, we perform multiple
iterations to adjust relative spacing, rotation and defocus of the 800 spots
with respect to the 1064 nm array using the SLM to maximise the relative
alignment of the two arrays by minimising atom survival. This also ensures the
800 nm light is focused on the same plane as the 1064 nm traps.

For performing local-light shifts, we suppress the mechanical effects of the
repulsive light shift on the ground-state atoms by changing the input beam
size on the SLM to create spots with an effective waist size of 5.1 $\mu$m,
resulting in $<2\%$ additional atom loss when applying local shifts of up to
10 MHz for periods of up to 3 $\mu$s.

Light-shift potentials are generated using the Gerchberg-Saxton algorithm
[52], with the modification that the relative weighting of each spot is
adjusted to reflect the target weighting $w_{i}$ of the specific graph
problem. Following the initial hologram generation, we perform Rydberg
spectroscopy at $n=50$ to minimise shifts due to interactions to measure the
resulting $\delta_{\mathrm{AC}}^{i}$ for each site. An example spectroscopy
measurement is shown in Fig. 6(a) for the case of the
$(w_{1},w_{2},w_{3},w_{4},w_{5})=(1,2,1,2,2)$ graph requiring $w_{\alpha}=2$
and $w_{\beta}=3$. This data shows that whilst sites initially defined with
higher weightings have larger shifts, there is a large RMS error between the
initial hologram generation and the target light shifts. To rectify this, we
implement closed-loop feedback on the relative weightings using the same
approach used for normalising relative intensity in trap arrays [52, 29].
Following each spectroscopy measurement, an additional hologram generation
step is taken using updated weights to drive convergence towards the target
light-shifts. Fig. 6(b) shows this can be used to obtain RMS errors $<2\%$
after 5 iterations, with the measured light-shifts after the final iteration
shown in Fig. 6(c) highlighting the ability to obtain accurate light shift
potentials.

###  I.4 Weighted-graph Annealing

Optimization of weighted graphs is performed using a dual-stage annealing
profile as illustrate in Fig. 1(c). To simplify optimization of the annealing
profile, we parameterize the detuning using a cubic ramp of the form [28]

$\delta(t)=at^{3}+bt+c,$ (3)

with ramp parameters $a=8sc/\tau^{3}$, $b=2c/\tau-a\tau^{2}/4$ and
$c=(\Delta_{\mathrm{max}}-\Delta_{\mathrm{min}})/2$. Here $0\leq s\leq 1$ acts
as a shape parameter, with $s=0$ equivalent to a linear ramp, and $\tau$ is
the total ramp duration.

This profile is then used to build piece-wise profiles, with the global
Rydberg laser detuning $\Delta(t)$ and 800 nm light-shift
$\delta_{\mathrm{AC}}$ given by

$\Delta(t)=\left\\{\begin{array}[]{ll}\delta(t)&\text{if }\delta(t)<0,\\\
0&\text{otherwise,}\end{array}\right.,\delta_{\mathrm{AC}}(t)=\left\\{\begin{array}[]{ll}0&\text{if
}\delta(t)<0,\\\ \delta(t)&\text{otherwise.}\end{array}\right.$ (4)

To synchronise the pulses, the parameters are sent to the same dual channel
AWG with one channel controlling amplitude and frequency of the 459 nm laser,
and the second channel connected to the 800 nm AOM situated before the SLM
that provides intensity control.

For the graphs shown above, annealing profile parameters $s$, $\tau$,
$\Delta_{\mathrm{min}}$ and $\Delta_{\mathrm{max}}$ are optimized using MLOOP
[61], an open source machine learning package interfaced through Python. For
each set of parameters, we perform 200-300 repeated measurements to evaluate
the averaged classical cost function $\langle H_{\mathrm{MWIS}}\rangle$ from
the distribution of experimentally observed bitstrings, and use this cost to
control the evolution of the optimization process. Following optimization, we
perform 1000 shot repeats to verify the final output probabilities presented
in the paper.


# A degree reduction method for an efficient QUBO formulation for the graph
coloring problem

Namho Hong1 ,  Hyunwoo Jung1 ,  Hyosang Kang2,∗ ,  Hyunjin Lim1 ,  Chaehwan
Seol3 and  Seokhyun Um1 1Taegu Science High School (TSHS), Daegu 42110, South
Korea.  2Daegu Gyeongbuk Institute of Science and Technology (DGIST), Daegu
42988, South Korea.  3Gwangju Science Academy (GSA), Gwangju 61005, South
Korea.  ∗ Corresponding author<EMAIL_ADDRESS>
###### Abstract.

We introduce a degree reduction method for symmetric polynomials on binary
variables. We also design an degree reduction algorithm for general
polynomials on binary variables, simulated on the graph coloring problem for
random graphs, and compared the results with the conventional method. The data
shows that our method produces quadratic polynomial of less variables than the
conventional method. The algorithm for our new degree reduction method is
robust, and applies to any QUBO formulation for quantum annealing systems.

This work was supported by TSHS R&E grant at Daegu Gyeongbuk Institute of
Science and Technology (DGIST). (grant number 2022030163)

Key words. Degree reduction, Graph coloring, QUBO, Quantum annealing

##  1\. Introduction

A graph is a 1-dimensional object that consists of vertices and edges. Two
vertices are called adjacent if they are connected by an edge, and a graph is
called simple if there is no multiple edges that joins adjacent vertices. A
vertex coloring is an assignment of “colors” to vertices such that no adjacent
vertices are assigned as the same color. The graph coloring problem is finding
the minimum number of colors that properly colors the graph. The graph
coloring problem has a wide range of applications such as scheduling [9],
register allocation in compilers [3], and frequency assignment in wireless
communications [2].

The graph coloring problem is NP-hard, meaning that there is no known
polynomial-time algorithm that solves the problem [5, 6]. One of the heuristic
methods for solving the graph coloring problem is the quantum annealing, which
uses quantum tunneling effect [1, 8]. To use quantum annealing system such as
D-Wave’s for solving graph coloring problems, we should first formulate the
problem as QUBO (quadratic unconstrained binary optimization) problem [10,
11]. This means the utility functions should be given as a polynomial on
binary variables of degree two. Since the utility functions are of higher
degree in general, we need to apply degree reduction methods to obtain a
quadratic polynomial [12].

The utility function for the graph coloring problem is formulated as follows.
Let $V$ and $E$ be the sets of vertices and edges in the graph $G$. Suppose
that we color each vertex $v\in V$ by a $m$-digit binary number
$\mathbf{x}_{v}=[x_{v,0}\cdots x_{v,m-1}]_{(2)}$. Given two vertices $v$ and
$w$, define the polynomial $P(v,w)$ as follows:

(1) $P_{v,w}=\prod_{k=0}^{m-1}(x_{v,k}x_{w,k}+(1-x_{v,k})(1-x_{w,k})).$

This polynomial is $0$ if and only if $x_{v,k}\neq x_{w,k}$ for any
$k=0,\cdots,m-1$. This means that it is zero only when $v$ and $w$ are colored
differently. Define the utility polynomial $Q(G)$ as follows:

(2) $Q_{G}=\sum_{(v,w)\in E}P_{v,w}.$

The minimum value of $Q_{G}$ is $0$ if and only if the graph $G$ is properly
colored. The degree of $Q_{G}$ is a polynomial is always greater than $2$
except for the trivial case, when the graph can be colored by a single color.
So we need to apply degree reduction methods on $Q_{G}$.

To formulate QUBO problem for the graph coloring problem on the graph $G$, we
need to apply degree reduction methods on the utility function $Q_{G}$. A
degree reduction on a polynomial $p$ of binary variables $\mathbf{x}$ is to
find a polynomial $q$ of lower degree that attains the same values of $p$ at
the same values for $\mathbf{x}$. More precisely, when
$\mathbf{x}=(x_{1},\cdots,x_{n})$ and $\mathbf{w}=(w_{1},\cdots,w_{d})$ are
binary variables, we say the polynomial $q(\mathbf{x},\mathbf{w})$ is a
reduced polynomial from $p(\mathbf{x})$ if $\deg q<\deg p$ and

(3) $p(\mathbf{x})=\min_{\begin{subarray}{c}w_{i}=0,1\\\
i=1,\ldots,d\end{subarray}}q(\mathbf{x},\mathbf{w}).$

We will call the variables in $\mathbf{w}$ as auxiliary variables. A method of
obtaining such $q$ is called a degree reduction method. Since we are
interested in the QUBO formulation, any degree reduction refers to the degree
reduction to quadratic polynomials. Freedman [4] and Ishikawa [7] introduced
degree reduction methods that reduces any monomials to degree $2$ polynomials.
Freedman’s method applies to negative monomials and Ishikawa’s method applies
to positive monomials. Ishikawa’s method splits into two cases when the
degrees are even and odd.

(4) $\displaystyle-\prod_{j=1}^{d}x_{j}$
$\displaystyle=\min_{w=0,1}w\left(d-\sum_{j=0}^{d}x_{j}\right),\ $ (5)
$\displaystyle\prod_{j=1}^{2d}x_{j}$ $\displaystyle=\sum_{1\leq i<j\leq
2d}x_{i}x_{j}+\min_{\begin{subarray}{c}w_{i}=0,1\\\
i=1,\ldots,d-1\end{subarray}}\left[\sum_{i=1}^{d-1}(4i-1)w_{i}-2\left(\sum_{i=1}^{d-1}w_{i}\right)\left(\sum_{j=1}^{2d}x_{j}\right)\right],$
(6) $\displaystyle\prod_{j=1}^{2d+1}x_{j}$ $\displaystyle=\sum_{1\leq i<j\leq
2d+1}x_{i}x_{j}+$ (7)
$\displaystyle\qquad\min_{\begin{subarray}{c}w_{i}=0,1\\\
i=0,\ldots,d\end{subarray}}\left[\sum_{i=1}^{d}(4i-1)w_{i}-2dw_{d-1}+\left(w_{d-1}-2\sum_{i=1}^{d}w_{i}\right)\left(\sum_{j=0}^{2d+1}x_{j}\right)\right].$

Given any polynomial, we can apply their methods to each monomial to obtain a
reduced polynomial of degree $2$. We will call this degree reduction method as
monomial reduction (MonoRed for short). Each time a degree reduction method is
applied to a monomial, new auxiliary variables are produced. The more
auxiliary variables we have, heavier loads are put on a quantum annealing
system, and the number of auxiliary variables easily exceeds the number of
original variables. For example, the utility polynomial $Q_{G}$ for the
complete graph $G$ of vertex size $8$ has $24$ variables and $1429$ monomials.
After the monomial degree reduction, there there are $1156$ new auxiliary
variables (c.f. Table 3).

In this paper, we propose a new method called the symmetric reduction (SymmRed
for short) that produces less number of auxiliary variables than the monomial
reduction. The main idea is to reduce symmetric polynomials. In $\S$2, we
state two main results (Theorem 1, 2) on reducing the symmetric polynomials of
positive and negative coefficients. We prove them in $\S$3, 4. In $\S$5, we
describe the algorithm of the symmetric reduction (c.f. Table 1, 2), and show
the efficiency on lowering the number of auxiliary variables in the symmetric
reduction. We tested the monomial degree reduction method and the symmetric
reduction method on random graphs and complete graphs of various vertex sizes.

##  2\. The main results

Let us redefine a binomial coefficient symbol $n\choose m$ for integer $n,m$
as followings:

(8) ${n\choose
m}=\left\\{\begin{array}[]{cl}\displaystyle\frac{n!}{m!(n-m)!}&\textrm{if
}m\leq n\\\ 0&\textrm{otherwise}\end{array}.\right.$

Let $S_{n}^{(m)}$ be the collection of all subsets of $\\{1,\cdots,n\\}$
consists of $m$ distinct integers. Let $P_{n}^{(m)}(\mathbf{x})$ be the
homogeneous polynomial of degree $m$ on $n$ variables
$\mathbf{x}=(x_{1},\ldots,x_{n})$:

$P_{n}^{(m)}(\mathbf{x})=\sum_{\\{i_{1},\cdots,i_{m}\\}\in
S_{n}^{(m)}}x_{i_{1}}\cdots x_{i_{m}}.$

We will call $P_{n}^{(m)}$ the symmetric $m$-polynomial on the variable
$\mathbf{x}$. For example, the symmetric $3$-polynomial on $5$ variables
$x_{1},\cdots,x_{5}$ is the following.

$\displaystyle P_{5}^{(3)}(\mathbf{x})$
$\displaystyle=x_{1}x_{2}x_{3}+x_{1}x_{2}x_{4}+x_{1}x_{2}x_{5}+x_{1}x_{3}x_{4}+x_{1}x_{3}x_{5}$
(9)
$\displaystyle\qquad+x_{1}x_{4}x_{5}+x_{2}x_{3}x_{4}+x_{2}x_{3}x_{5}+x_{2}x_{4}x_{5}+x_{3}x_{4}x_{5}.$

The following two theorems are the keys to to the symmetric reduction.

######  Theorem 1.

Let $n>m>2$ and

(10) $L(n,m)=\displaystyle{n-2\choose m-2},\quad
d=\displaystyle\left\lfloor\frac{n-1}{2}\right\rfloor.$

For $1\leq i\leq d$, define

(11) $\displaystyle a_{i}(n,m)=$ $\displaystyle(4i-1){n-2\choose
m-2}-2i{2i\choose m-1}+2i{2i-2\choose m-1}+{2i-2\choose m-2},$ (12)
$\displaystyle b_{i}(n,m)=$ $\displaystyle-2{n-2\choose m-2}+{2i\choose
m-1}-{2i-2\choose m-1}.$

Then

(13) $P_{n}^{(m)}(\mathbf{x})=L(n,m)\sum_{1\leq i<j\leq
n}x_{i}x_{j}+\min_{w_{i}\in\\{0,1\\}}\left[\sum_{i=1}^{d}w_{i}\left(a_{i}(n,m)+b_{i}(n,m)\sum_{j=1}^{n}x_{j}\right)\right].$

######  Theorem 2.

Let $n>m>2$ and $d=\displaystyle\left\lfloor\frac{n-m+2}{2}\right\rfloor$. For
$0\leq i\leq d$, define

(14) $\displaystyle a_{i}(m)$ $\displaystyle=(m-1)\left({m+2i-1\choose
m}-{m+2i-3\choose m}\right),$ (15) $\displaystyle b_{i}(m)$
$\displaystyle={m+2i-4\choose m-1}-{m+2i-2\choose m-1}.$

Then

(16)
$-P_{n}^{(m)}(\mathbf{x})=\min_{w_{i}\in\\{0,1\\}}\sum_{i=1}^{d}w_{i}\left(a_{i}(m)+b_{i}(m)\sum_{j=1}^{n}x_{j}\right).$

$l(\mathbf{x})$ | $P_{5}^{(3)}(\mathbf{x})$ | $\displaystyle\sum_{i<j}x_{i}x_{j}$ | $A_{l}$ | $\displaystyle\min_{w_{1}=0,1}w_{1}(7-5l(\mathbf{x}))$ | $\displaystyle\min_{w_{2}=0,1}w_{2}(3-l(\mathbf{x}))$
---|---|---|---|---|---
$0$ | $0$ | $0$ | $0$ | $0$ | $0$
$1$ | $0$ | $0$ | $0$ | $0$ | $0$
$2$ | $0$ | $1$ | $-3$ | $-3$ | $0$
$3$ | $1$ | $3$ | $-8$ | $-8$ | $0$
$4$ | $4$ | $6$ | $-14$ | $-13$ | $-1$
$5$ | $10$ | $10$ | $-20$ | $-18$ | $-2$
Table 1. The sequence $A_{l}$ in Equation (18) when $c=3$ and the relative
arithmetic progressions.

Let us explain the underlying ideas in Theorem 1 with the example of
$P_{5}^{(3)}(\mathbf{x})$. The goal is to find a quadratic polynomial that
attains the same value of $P_{5}^{(3)}$ with some auxiliary variables
$w_{1},\ldots,w_{d}$. The only possible way to formulate such quadratic
polynomial is the following:

(17) $P_{5}^{(3)}(\mathbf{x})=\sum_{1\leq i<j\leq
5}c_{ij}x_{i}x_{j}+\sum_{i}d_{i}x_{i}+\min_{w_{j}=0,1}\sum_{j=1}^{d}w_{j}\left(a_{j}+\sum_{i=1}^{5}b_{i,j}x_{i}\right).$

In fact, this is the general form of any degree reduction. We omitted the
constant term in Equation (17) since the minimum $P_{n}^{(m)}$ is always $0$.
Let us assume that $b_{i,j}=b_{j}$ only depends on $i$, and $c_{i,j}=c$,
$d_{i}=0$ are constants. Let $l=l(\mathbf{x})=x_{1}+\ldots x_{5}$ be the
arithmetic sum (not the binary sum) of all values in the variable
$\mathbf{x}$. Then the sequence $A_{l}$ defined by

(18) $A_{l}=P_{5}^{(3)}(\mathbf{x})-c\sum_{i,j}x_{i}x_{j}$

is the sum of arithmetic progressions on $l$:

(19)
$A_{l}=\min_{w_{j}=0,1}\left[\sum_{j=1}^{d}w_{j}\left(a_{j}+b_{j}l\right)\right].$

We observe that

(20) $P_{5}^{(3)}(\mathbf{x})={l\choose 3},\quad\sum_{i<j}x_{i}x_{j}={l\choose
2}.$

When $c=3$, the sequence $A_{0},A_{1},\ldots,A_{5}$ becomes a decreasing
sequence of non-positive integers as in Table 1. Such sequence can always be
expressed as the sum of multiple arithmetic progressions of decreasing non-
positive integers. With $\min$ and the auxiliary variables $w_{j}$, the value
of arithmetic progression $a_{j}+b_{j}l$ never exceed $0$.

Table 1 shows the values of $A_{l}$ and the arithmetic progressions that
exhaust $A_{l}$ when added together. The first two nonzeros $-3,-8$ of $A_{l}$
follows the arithmetic progression $7-5l$. By subtracting the sequence $\min
w_{1}(7-5l)$ from $A_{l}$ index wise, we get a sequence that has two negatives
$-1,-2$ at $l=4,5$. We can eliminate this by the sequence $\min w_{2}(3-l)$.
Thus we get the following equation.

(21)
$P_{5}^{(3)}(\mathbf{x})=3\sum_{i<j}x_{i}x_{j}+\min_{w_{1},w_{2}=0,1}\left[w_{1}(7-5\sum_{i=1}^{5}x_{i})+w_{2}(3-\sum_{i=1}^{5}x_{i})\right].$

$l(\mathbf{x})$ | $-P_{5}^{(3)}(\mathbf{x})$ | $\displaystyle\min_{w_{1}=0,1}w_{1}(8-3l(\mathbf{x}))$ | $\displaystyle\min_{w_{2}=0,1}w_{2}(12-3l(\mathbf{x}))$
---|---|---|---
$0$ | $0$ | $0$ | $0$
$1$ | $0$ | $0$ | $0$
$2$ | $0$ | $0$ | $0$
$3$ | $-1$ | $-1$ | $0$
$4$ | $-4$ | $-4$ | $0$
$5$ | $-10$ | $-7$ | $-3$
Table 2. The values of $-P_{5}^{(3)}(\mathbf{x})$ and the relative arithmetic
progressions.

Now let us explain the idea of Theorem 2 with the example of
$-P_{5}^{(3)}(\mathbf{x})$. Let $l=l(\mathbf{x})$ be the arithmetic sum of
$\mathbf{x}$ as before, and define $B_{l}=-P_{5}^{(3)}(\mathbf{x})$. The
sequence $B_{l}$ is already a decreasing sequence of negative integers as
shown in Table 2. The first two negatives $-1,-4$ follows the arithmetic
progression $8-3l$. Subtracting the sequence $\min w_{1}(8-3l)$ from $B_{l}$
leaves two values $0,-3$ at $l=4,5$. This can be exhausted by the sequence
$\min w_{2}(12-3l)$. Thus we get the following equation.

(22)
$-P_{5}^{(3)}(\mathbf{x})=\min_{w_{1},w_{2}=0,1}\left[w_{1}(8-3\sum_{i=1}^{5}x_{i})+w_{2}(12-3\sum_{i=1}^{5}x_{i}).\right]$

In the next section, we will prove Theorem 1, 2 for the general symmetric
polynomial $P_{n}^{(m)}$. It is worth to point two facts.

  1. (1)

There are $n-1$ negatives in the sequence $A_{l}$. Thus we need
$\displaystyle\lfloor(n-1)/2\rfloor$ auxiliary variables for $A_{l}$. Also,
the coefficients $a_{j}$ and $b_{j}$ are dependent to both $n$ and $m$.

  2. (2)

There are $n-m+2$ negatives in $B_{l}$, so we need
$\displaystyle\lfloor(n-m+2)/2\rfloor$ auxiliary variables. We do not need the
constant $n$ to determine $a_{j}$ and $b_{j}$, because the sequence $B_{l}$ is
determined by $m$ only.

##  3\. Proof of Theorem 1

Let $\mathbf{x}=(x_{1},\cdots,x_{n})$ and

(23) $l=l(\mathbf{x})=\sum_{i=1}^{n}x_{i}$

be the sum of all values in the variables $x_{1},\ldots,x_{n}$. Our first task
is to find a value $c=c(n,m)$ so that

(24) $A_{l}=P_{n}^{(m)}(\mathbf{x})-c\displaystyle\sum_{i<j}x_{i}x_{j}$

becomes a decreasing sequence on $l$. For each $i\geq 0$ satisfying $n-i-2\geq
m$, we observe that

(25) $\displaystyle A_{n-i}$ $\displaystyle={n-i\choose m}-c{n-i\choose 2},$
(26) $\displaystyle A_{n-i-1}$ $\displaystyle={n-i-1\choose m}-c{n-i-1\choose
2},$ (27) $\displaystyle A_{n-i-2}$ $\displaystyle={n-i-2\choose
m}-c{n-i-2\choose 2}.$

We want these three values become a decreasing sequence for any such $i$. This
is possible only when

(28) $\displaystyle(A_{n-i}-A_{n-i-1})-(A_{n-i-1}-A_{n-i-2})$ (29)
$\displaystyle=\left({n-i-1\choose m-1}-{n-i-2\choose
m-1}\right)-c\left({n-i-1\choose 1}-{n-i-2\choose 1}\right)$ (30)
$\displaystyle={n-i-2\choose m-2}-c\leq 0.$

Thus we get

(31) $c=\max_{i}{n-i-2\choose m-2}={n-2\choose m-2}.$

Next, let us find $a_{j},b_{j}$ satisfying Equation (19). For $i\geq 1$,

(32) $\displaystyle{2i\choose m}-{n-2\choose m-2}{2i\choose 2}$
$\displaystyle=\sum_{k=1}^{i}a_{k}+\sum_{k=1}^{i}2kb_{1},$ (33)
$\displaystyle{2i+1\choose m}-{n-2\choose m-2}{2i+1\choose 2}$
$\displaystyle=\sum_{k=1}^{i}a_{k}+\sum_{k=1}^{i}(2k+1)b_{1}.$

Thus we get

(34) $\sum_{k=1}^{i}b_{k}={2i\choose m-1}-2i{n-2\choose m-2}.$

This gives the formula (12).

Next, let us take (32) + (33) to get

(35) ${2i\choose m}+{2i+1\choose m}-{n-2\choose m-2}\left({2i\choose
2}+{2i+1\choose 2}\right)=2\sum_{k=1}^{i}a_{k}+(4i+1)\sum_{k=1}^{i}b_{k}$

Applying Equation (34), we get

(36) $\sum_{k=1}^{i}a_{k}=2{n-2\choose m-2}{2i+1\choose m}-(m-1){2i+1\choose
m}.$

This derives Equation (11).

##  4\. Proof of Theorem 2

With the same notation for $l$ in Equation (23), let us define
$B_{l}=-P_{n}^{(m)}(\mathbf{x})$. For each $j\geq 0$, We want to find
$a_{j},b_{j}$ that satisfies

(37)
$B_{l}=\min_{w_{j}=0,1}\left[\sum_{j=1}^{d}w_{j}\left(a_{j}+b_{j}l\right)\right].$

For each $i\geq 1$,

(38) $\displaystyle-{m+2i-2\choose m}=$
$\displaystyle\sum_{k=1}^{i}a_{k}+(m+2j-2)\sum_{k=1}^{i}b_{k},$ (39)
$\displaystyle-{m+2i-1\choose m}=$
$\displaystyle\sum_{k=1}^{i}a_{k}+(m+2j-1)\sum_{k=1}^{i}b_{k}.$

By taking (38)-(39), we get

(40) $-{m+2j-2\choose m-1}=\sum_{k=1}^{j}b_{k}.$

Thus we get Equation (15). By taking (38)+(39), we get

(41) $2\sum_{k=1}^{i}a_{k}+(2m+4i-3)\sum_{k=1}^{i}b_{k}=-{m+2i-1\choose
m}-{m+2i-2\choose m}.$

Applying Equation (40), we get

(42) $\sum_{k=1}^{i}a_{k}=(m-1){m+2i-1\choose m}.$

Therefore, we get Equation (14).

##  5\. Simulated results

In this section, we introduce two algorithms, MaxSymm and SymmRed. MaxSymm is
an algorithm that finds a maximal symmetric polynomial in a given polynomial.
SymmRed is a degree reduction algorithm that reply on MaxSymm and MonoRed
(c.f. $\S\ref{sec:intro}$).

Algorithm 1 MaxSymm algorithm

1:procedure MaxSymm($p$, $m$)$\triangleright$ $p$ is a polynomial, $m\geq 3$
is a degree

2: Define $A:=B:=C:=$ the list of all $m$-monomials in $p$

3: while $C$ is not empty do

4: $B\leftarrow C$

5: $C\leftarrow$ the empty list (of polynomials)

6: for $q$ in $B$ do

7: for $r$ in $A$ do

8: $S\leftarrow$ the set of all variables in $q$ and $r$

9: $s\longleftarrow$ the symmetric $m$-polynomial on $S$

10: if all monomials in $s$ lie in $p$ then

11: $C$.append($s$)

12: end if

13: end for

14: end for

15: end while

16: return $B$

17:end procedure

1:procedure SymmRed($p$)$\triangleright$ $p$ is a polynomial

2: $q\leftarrow 0$

3: for $m\leftarrow 3$ to $\deg p$ do

4: $s\leftarrow$ MaxSymm($p$, $m$)

5: while $s$ is not $0$ do

6: $a\leftarrow$ the most common coefficient of monomials in $p$ and $s$

7: $p\leftarrow p-a\cdot s$

8: if $a>0$ then

9: $r\leftarrow$ apply Theorem 1 to $s$

10: else

11: $r\leftarrow$ apply Theorem 2 to $-s$

12: end if

13: $q\leftarrow q+|a|\cdot r$

14: $s\leftarrow$ MaxSymm($p$, $m$)

15: end while

16: end for

17: $q\leftarrow q+$MonoRed(p) $\triangleright$ do the monomial reduction

18:end procedure

Algorithm 2 The Symmetric Reduction algorithm

Algorithm 1 is the pseudocode for MaxSymm algorithm. We will say a polynomial
$q$ lies in a polynomial $p$ if all monomials in $q$ are monomials in $p$. We
will say a symmetric $m$-polynomial $q$ is maximal in $p$ if $q$ lies in $p$
and there is no larger symmetric $m$-polynomial that lies in $p$. In this
section, we introduce algorithm to find a maximal symmetric polynomials in
$p$.

Algorithm 2 is the pseudocode for SymmRed algorithm. For each symmetric
polynomial found by this algorithm, we apply Theorem 1 or 2. (c.f. Line 9, 11
in Algorithm 2) We apply the monomial degree reduction (by Freedman and
Ishikawa) on the rest of the parts.

Line 6 in Algorithm 2) defines a coefficient $a$ of a symmetric polynomial. We
choose $a$ so that we can eliminate as many monomials as possible in Line 7.
However, we may have more one candidate for $a$. For example, the polynomial
below give two choices for $a$ of different signs.

$\displaystyle-
x_{0}x_{1}x_{2}-x_{0}x_{1}x_{3}+x_{0}x_{2}x_{3}+x_{1}x_{2}x_{3}$
$\displaystyle=\underbrace{-x_{0}x_{1}x_{2}-x_{0}x_{1}x_{3}-x_{0}x_{2}x_{3}-x_{1}x_{2}x_{3}}_{\textrm{Apply
Theorem
}\ref{thm:main-1}}+\underbrace{2x_{0}x_{2}x_{3}+2x_{1}x_{2}x_{3}}_{\textrm{monomial-
wise reduction}}$ (43)
$\displaystyle=\underbrace{x_{0}x_{1}x_{2}+x_{0}x_{1}x_{3}+x_{0}x_{2}x_{3}+x_{1}x_{2}x_{3}}_{\textrm{Apply
Theorem
}\ref{thm:main-2}}+(\underbrace{-2x_{0}x_{1}x_{2}-2x_{0}x_{1}x_{3}}_{\textrm{monomial-
wise reduction}})$

We tested the Symmetric Reduction algorithm on random $p$-graphs. We varied
the number of vertices $V=3,4,5,6,7,8$ and the probabilities $p=0.75$, $0.80$,
$0.85$, $0.90$, $0.95$, and $1.00$. We obtained the average number of
variables $r_{i}$ and monomials $N_{i}$ of the reduced polynomials for both
reduction methods where $i=1$ for symmetric reduction and $i=2$ for monomial
reduction. Note that when $p=1.00$ we get the complete graph, so the numbers
are fixed. In the simulation, we set the number of colors to be
$d=\displaystyle\lceil\log_{2}V\rceil$.

The results are shown in Table 3 and Figure 1. On average, the numbers of
variables produced by the symmetric reduction method are approximately
$38.90\%\sim 35.95\%$ of that by the monomial reduction method. The reduced
polynomial obtained by the symmetric reduction method approximately contains
$41.12\%\sim 43.35\%$ of monomials that produced by the monomial reduction
method.

(a) $p=0.75$

(b) $p=0.80$

(c) $p=0.85$

(d) $p=0.90$

(e) $p=0.95$

(f) $p=1.00$

Figure 1. Graphs of the average numbers of variables (vertical) on the numbers of vertices (horizontal) of random $p$-graphs. The squared-dots are for the monomial reduction where as the rounded-dots are for the the symmetric reduction. Table 3. The average numbers of variables $r_{i}$ and monomials $N_{i}$ from the symmetric ($i=1$) and monomial ($i=2$) reduction methods on random $p$-graphs with $V$ vertices and the ratios. $p$ | $V$ | Symmetric reduction | Monomial reduction | $r_{1}/r_{2}$ | $N_{1}/N_{2}$
---|---|---|---|---|---
$r_{1}$ | $N_{1}$ | $r_{2}$ | $N_{2}$
$0.75$ | 3 | 10.01 | 40.41 | 16.67 | 64.81 | 60.07% | 62.34%
4 | 16.92 | 76.14 | 30.49 | 125.90 | 55.49% | 60.47%
5 | 132.63 | 706.06 | 329.15 | 1604.47 | 40.30% | 44.01%
6 | 191.24 | 1033.05 | 483.05 | 2367.48 | 39.59% | 43.63%
7 | 265.02 | 1446.57 | 675.62 | 3324.94 | 39.23% | 43.51%
8 | 344.87 | 1896.55 | 886.97 | 4376.71 | 38.90% | 43.35%
$0.80$ | 3 | 10.59 | 43.35 | 17.81 | 69.83 | 59.45% | 62.08%
4 | 17.54 | 80.14 | 31.95 | 132.98 | 54.90% | 60.27%
5 | 140.79 | 757.57 | 350.27 | 1710.45 | 40.19% | 44.00%
6 | 204.09 | 1107.13 | 517.53 | 2540.70 | 39.44% | 43.58%
7 | 277.52 | 1518.67 | 709.22 | 3493.69 | 39.13% | 43.47%
8 | 368.10 | 2030.41 | 948.97 | 4688.19 | 38.79% | 43.31%
$0.85$ | 3 | 10.95 | 45.38 | 18.57 | 73.34 | 58.94% | 61.87%
4 | 18.20 | 84.57 | 33.55 | 140.83 | 54.26% | 60.05%
5 | 147.75 | 792.82 | 368.94 | 1804.23 | 40.05% | 43.94%
6 | 215.48 | 1173.98 | 550.00 | 2703.81 | 39.18% | 43.42%
7 | 296.99 | 1629.53 | 759.03 | 3743.98 | 39.13% | 43.52%
8 | 387.75 | 2144.74 | 1003.50 | 4962.17 | 38.64% | 43.22%
0.90 | 3 | 11.33 | 47.66 | 19.42 | 77.32 | 58.35% | 61.64%
4 | 18.80 | 88.62 | 35.00 | 148.05 | 53.70% | 59.86%
5 | 155.46 | 837.38 | 389.48 | 1907.39 | 39.92% | 43.90%
6 | 225.05 | 1231.30 | 579.07 | 2849.89 | 38.86% | 43.21%
7 | 314.98 | 1732.93 | 806.38 | 3981.91 | 39.06% | 43.52%
8 | 409.07 | 2269.40 | 1063.36 | 5262.93 | 38.47% | 43.12%
0.95 | 3 | 11.68 | 49.84 | 20.22 | 81.15 | 57.76% | 61.41%
4 | 19.39 | 92.74 | 36.48 | 155.38 | 53.16% | 59.68%
5 | 163.19 | 881.79 | 409.50 | 2008.00 | 39.85% | 43.91%
6 | 231.97 | 1278.27 | 608.69 | 2998.73 | 38.11% | 42.63%
7 | 328.62 | 1813.48 | 845.34 | 4177.61 | 38.87% | 43.41%
8 | 425.13 | 2374.63 | 1126.03 | 5577.85 | 37.90% | 42.69%
1.00 | 3 | 12 | 52 | 21 | 85 | 57.14% | 61.18%
4 | 20 | 97 | 38 | 163 | 52.63% | 59.51%
5 | 171 | 927 | 430 | 2111 | 39.77% | 43.91%
6 | 234 | 1306 | 639 | 3151 | 36.62% | 41.45%
7 | 345 | 1909 | 889 | 4397 | 38.81% | 43.42%
8 | 424 | 2405 | 1180 | 5849 | 35.93% | 41.12%

##  6\. Conclusion

One of the main limitation that the current quantum annealing systems have in
solving real-world problems is the number of qubits available on the system.
Thus the optimizing the number of variables in a QUBO formulation is
essential. We showed that there is alternative degree reduction method that
produces less number of variables during the QUBO formulation. The algorithm
we implemented for degree reduction is robust. The algorithm works better on
polynomials that has higher symmetry. We expect that we can apply this method
to any QUBO formulations of various problems.

## References

  * [1] B. Apolloni, C. Carvalho, and D. de Falco.  Quantum stochastic optimization.  Stochastic Processes and their Applications, 33(2):233–244, 1989\.
  * [2] Balabhaskar Balasundaram and Sergiy Butenko.  Graph Domination, Coloring and Cliques in Telecommunications, pages 865–890.  Springer US, Boston, MA, 2006.
  * [3] G. J. Chaitin.  Register allocation & spilling via graph coloring.  SIGPLAN Not., 17(6):98–101, jun 1982.
  * [4] D. Freedman and P. Drineas.  Energy minimization via graph cuts: settling what is possible.  In 2005 IEEE Computer Society Conference on Computer Vision and Pattern Recognition (CVPR’05), volume 2, pages 939–946 vol. 2, 2005.
  * [5] Michael R. Garey and David S. Johnson.  Computers and Intractability; A Guide to the Theory of NP-Completeness.  W. H. Freeman & Co., USA, 1990.
  * [6] M.R. Garey, D.S. Johnson, and L. Stockmeyer.  Some simplified np-complete graph problems.  Theoretical Computer Science, 1(3):237–267, 1976.
  * [7] Hiroshi Ishikawa.  Transformation of general binary mrf minimization to the first-order case.  IEEE Transactions on Pattern Analysis and Machine Intelligence, 33(6):1234–1249, 2011.
  * [8] Tadashi Kadowaki and Hidetoshi Nishimori.  Quantum annealing in the transverse ising model.  Phys. Rev. E, 58:5355–5363, Nov 1998.
  * [9] Frank Thomson Leighton.  A graph coloring algorithm for large scheduling problems.  Journal of research of the National Bureau of Standards, 84 6:489–506, 1979.
  * [10] Carla Silva, Ana Aguiar, Priscila M. V. Lima, and Inês Dutra.  Mapping graph coloring to quantum annealing.  Quantum Machine Intelligence, 2(2):16, 2020.
  * [11] Z. Tabi, K. H. El-Safty, Z. Kallus, P. Haga, T. Kozsik, A. Glos, and Z. Zimboras.  Quantum optimization for the graph coloring problem with space-efficient embedding.  In 2020 IEEE International Conference on Quantum Computing and Engineering (QCE), pages 56–62, Los Alamitos, CA, USA, oct 2020. IEEE Computer Society.
  * [12] Salvador E. Venegas-Andraca, William Cruz-Santos, Catherine McGeoch, and Marco Lanzagorta.  A cross-disciplinary introduction to quantum annealing-based algorithms.  Contemporary Physics, 59(2):174–197, 2018.