# Vector Ising Spin Annealer for Minimizing Ising Hamiltonians

James S. Cummins and Natalia G. Berloff<EMAIL_ADDRESS>Department of Applied Mathematics and Theoretical Physics, University of
Cambridge, Wilberforce Road, Cambridge CB3 0WA, United Kingdom

(September 3, 2024)

###### Abstract

We introduce the Vector Ising Spin Annealer (VISA), a framework in gain-based
computing that harnesses light-matter interactions to solve complex
optimization problems encoded in spin Hamiltonians. Traditional driven-
dissipative systems often select excited states due to limitations in spin
movement. VISA transcends these constraints by enabling spins to operate in a
three-dimensional space, offering a robust solution to minimize Ising
Hamiltonians effectively. Our comparative analysis reveals VISA’s superior
performance over conventional single-dimension spin optimizers, demonstrating
its ability to bridge substantial energy barriers in complex landscapes.
Through detailed studies on cyclic and random graphs, we show VISA’s
proficiency in dynamically evolving the energy landscape with time-dependent
gain and penalty annealing, illustrating its potential to redefine
optimization in physical systems.

##  I Introduction

In pursuing advancements in digital computing, mainly aimed at addressing the
complexities inherent in AI and large-scale optimization problems, the
inherent limitations of the traditional von Neumann architecture come to the
forefront. These limitations are characterized by the escalating costs
required for incremental performance improvements, prompting a pivotal shift
toward more sustainable computational paradigms [1]. In this evolving
landscape, two primary paths have emerged: refining algorithms within the
existing computational frameworks and exploring novel hardware paradigms
grounded in physics principles, with a notable focus on exploiting light-
matter interactions.

Algorithmic enhancements, while valuable, offer incremental improvements. In
contrast, the exploration of alternative hardware paradigms holds the promise
of a substantial shift, leveraging the fundamental principles of physics—such
as minimizing entropy, energy, and dissipation [2]—and integrating quantum
phenomena like superposition and entanglement [3]. This innovative approach
aims to transcend the conventional computing model’s limitations, tapping into
the intrinsic computational potential of physical systems to address complex
optimization challenges currently beyond traditional methods’ reach.

Central to this innovative trajectory is the integration of analogue, physics-
based algorithms and hardware, which involve translating complex optimization
problems into universal spin Hamiltonians, like the classical Ising and XY
models [4, 5, 6]. This translation process involves embedding the structure of
the problem into the spin Hamiltonian’s coupling strengths and targeting the
ground state as the solution, with the physical system tasked with discovering
this state. The efficiency and accuracy of this mapping process are vital, as
they ensure the scalability of computational efforts, enabling these problems
to remain manageable despite increasing complexity [7].

Gain-based computing (GBC) using coupled light-matter emerges as a novel
disruptive computing platform distinct from gate-based and quantum or
classical annealing approaches. GBC combines optical light manipulations with
advancements in laser technology and spatial light modulators (SLMs),
facilitating parallel processing across multiple channels with significant
nonlinearities and high-energy efficiency. The operational principle of
GBC—increasing pumping power (annealing), followed by symmetry breaking and
gradient descent—relies on wave coherence and synchronization to lead the
system to a state that minimizes losses naturally. Recent developments in
physics-based hardware exploiting GBC principles have showcased diverse
technologies. These range from optical parametric oscillator (OPO) based
coherent Ising machines (CIMs) [8, 9, 10, 11], memristors [12], and laser
systems [13, 14, 15], to photonic simulators [16, 17], polaritons [5, 6],
photon condensates [18], and novel computational architectures like
Microsoft’s analogue iterative machine [19] and Toshiba’s simulated
bifurcation machine [20]. These platforms have been instrumental in minimizing
programmable spin Hamiltonians, demonstrating efficacy across a spectrum of
complex optimization challenges, including machine learning [21], financial
analytics [22], and biophysical modelling [23, 24].

The Ising Hamiltonian, originating from statistical physics, is a natural
model to be implemented in these platforms. It describes interactions between
spins in a lattice, providing a fundamental model for understanding complex
systems and can be written as

$H_{\rm I}=-\sum_{ij}J_{ij}s_{i}s_{j}-\sum_{i}h_{i}s_{i},$ (1)

where $s_{i}=\pm 1$ represents the spin, $J_{ij}$ the interaction strength
between $i$-th and $j$-th spins, and $h_{i}$ the external magnetic field. The
Ising model’s exploration of ground states in physical systems parallels the
search for minimum cost functions in optimization, while its energy landscapes
mirror the loss landscapes in machine learning, offering a unique perspective
on optimization problems.

At the core of the hardware representation of Ising Hamiltonians lies the
utilization of scalar soft-spins within soft-spin Ising models (SSIMs), which
effectively change energy barriers inherent in classical hard-spin Ising
Hamiltonians, $H_{\rm I}$. The transition from the hard-spin Ising
Hamiltonians $H_{\rm I}$ to the SSIMs consists of considering binary spins
$s_{i}$ as signs of continuous variables $x_{i}$ \- amplitudes - that
bifurcate from vacuum state ${\bf x}=0$ guided by the gain increase. While
facilitating enhanced problem-solving capabilities, this approach occasionally
encounters the obstacle of trajectory trapping within local minima due to the
escalating energy barriers as gain increases.

Therefore, the challenge of navigating energy landscapes with numerous local
minima, especially during amplitude bifurcation, necessitates an innovative
approach to avoid the system being trapped in local minima. In this paper, we
introduce the Vector Ising Spin Annealer (VISA), a model that integrates the
benefits of multidimensional spin systems and soft-spin gain-based evolution.
Unlike traditional single-dimension semi-classical spin models, VISA utilizes
three soft modes to represent vector components of an Ising spin, allowing for
movement in three dimensions and offering a robust framework for accurately
determining ground states in complex optimization problems. Our comparative
studies underscore the enhanced performance of the VISA compared to
established models with scalar spins such as Hopfield-Tank networks, coherent
Ising machines, and the spin-vector Langevin model (SVL). VISA showcases a
notable ability to overcome significant energy barriers and effectively
navigate through intricate energy landscapes. We illustrate VISA’s
effectiveness by tackling Ising Hamiltonian problems across a range of graph
structures. Our examples include analytically solvable 3-regular circulant
graphs, more complex circulant graphs, and random graphs where minimizing the
Ising Hamiltonian is known to be an NP-hard problem. In our discussions, when
we refer to “solving a graph A” or “minimizing a graph A”, we are using
shorthand for the more detailed process of “finding the global minimum of the
Ising Hamiltonian on graph A”. This terminology simplifies our reference to
the comprehensive task of determining the lowest possible energy state for the
Ising model applied to a specific graph structure.

##  II Vector Ising Spin Annealer

The VISA model is a semi-classical, three-dimensional soft-spin Ising model.
It employs annealing of the loss landscape, symmetry-breaking, bifurcation,
gradient descent, and mode selection to drive the system toward the global
minimum of the Ising Hamiltonian. Here, we represent Ising soft-spins as
continuous vectors in three-dimensional space ${\bf
x}_{i}=(x_{1i},x_{2i},x_{3i})$ that are free to move in that space. VISA may
be physically realized using amplitudes of network of coupled optical
oscillators. In these optical-based Ising machines, vectors of Ising spins are
represented by three soft-spin amplitudes. In total, $3N$ amplitudes are
required to minimize an $N$-spin Ising Hamiltonian. This approach is somewhat
complimentary to the recently proposed “dimensionality annealing” where soft
Ising amplitudes are considered as coordinates of the multidimensional spins
such as XY or Heisenberg spins [25, 26]. While the main idea is to exploit the
advantages of the higher dimensionality, VISA is aiming at the Ising
minimization rather than hyper-dimensional spin systems [25, 26].

To formulate the Hamiltonian that is capable of representing the dynamics of
the GBC, we incorporate the term that is convex and dominates when the
effective gain (loss) is large and negative, the Ising term that is minimized
when gain is large and positive, and a term that aligns the spins at the end
of the process so that their direction can be associated with the binary spin
$s_{i}=\pm 1$. We write, therefore, the VISA Hamiltonian as the sum of three
terms $H_{\rm VISA}=H_{1}+H_{2}+H_{3}$, where

$\displaystyle H_{1}$
$\displaystyle=\frac{1}{4}\sum_{i=1}^{N}(\gamma_{i}(t)-|{\bf
x}_{i}|^{2})^{2},$ (2) $\displaystyle H_{2}$
$\displaystyle=-\frac{1}{2}\sum_{i,j=1}^{N}J_{ij}{\bf x}_{i}\cdot{\bf x}_{j},$
(3) $\displaystyle H_{3}$ $\displaystyle=\frac{P(t)}{2}\sum_{i,j=1}^{N}|{\bf
x}_{i}\times{\bf x}_{j}|^{2}.$ (4)

As the effective gain $\gamma_{i}(t)$ increases with time $t$ from negative
(effective losses) to positive (effective gain) values, $H_{1}$ anneals
between a convex function with minimum at ${\bf x}_{i}={\bf 0}$ to nonzero
amplitudes. $H_{2}$ coincides with the Ising Hamiltonian when all ${\bf
x}_{i}$ have unit magnitude when projected along the same direction. Finally,
$H_{3}$ is a penalty term with time-dependent magnitude $P(t)$ to enforce
collinearity between ${\bf x}_{i}$ and ${\bf x}_{j}$ at the end of the gain-
induced landscape change. At which time, the condition on the amplitudes
$|{\bf x}_{i}|=1$ is imposed by the feedback on the gain realised by
$\dot{\gamma}_{i}=\varepsilon(1-|{\bf x}_{i}|^{2})$ [27]. Analogously to CIM
operation, as $\gamma_{i}(t)$ grows from negative to positive values, $H_{\rm
VISA}$ anneals from the dominant convex term $H_{1}$ that is minimized at
${\bf x}_{i}=(0,0,0)$ for all $i$, to the minimum of $H_{2}+H_{3}$ with $|{\bf
x}_{i}|=1$ via symmetry-breaking bifurcation. Concurrently, as $P(t)$
increases from $P(0)=0$ to sufficiently large $P(T)>0$, the soft-spins vectors
become collinear. The Ising spins are calculated by projecting the three-
dimensional vectors along the axis ${\bf k}$ they have centred around and
taking signs of the resultant scalar $s_{i}=\mathrm{sgn}({\bf x}_{i}\cdot{\bf
k})$. At $t=T$, the target hard-spin Ising Hamiltonian $H_{\rm I}$ is
minimized.

The Hamiltonian $H_{\rm VISA}$ is $4$-local due to the $H_{3}$ term. Some
optical hardware can directly encode such interactions [28]. At the same time,
it is possible to reduce the $4$-local Hamiltonian to quadratic Hamiltonian by
substituting the product of two variables with a new one while adding the
appropriate penalty terms. For instance, the quartic interaction term
$x_{i}x_{j}x_{k}x_{m}$ can be reduced to a triplet by introducing a new
variable $y$ instead of $x_{k}x_{m}$ and the term
$q_{km}(y)=Q(x_{k}x_{m}-2x_{k}y-2x_{m}y+3y)$ (known as Rosenberg polynomial
[29]) in the Hamiltonian where $Q>0$ is the penalty constant. Note that
$q_{km}=0$ if and only if $x_{k}x_{m}=y$ and is strictly positive otherwise.
The reduction of the triplet $x_{i}x_{j}y$ is done similarly by replacing
$x_{i}x_{j}$ with a new spin variable $z$ and introducing the penalty term
$q_{ij}(z)$ in the Hamiltonian [30].

The governing equations $\dot{{\bf x}}_{i}=-\nabla_{i}H_{\rm VISA}$ represent
the gradient descent combined with the temporal change of annealing parameters
$\gamma_{i}(t)$ and $P(t)$ as

$\begin{split}{\bf\dot{x}}_{i}&={\bf x}_{i}(\gamma_{i}(t)-|{\bf
x}_{i}|^{2})+\sum_{j=1}^{N}J_{ij}{\bf x}_{j}\\\ &-P(t)\Big{(}{\bf
x}_{i}\sum_{j=1}^{N}|{\bf x}_{j}|^{2}-\sum_{j=1}^{N}{\bf x}_{j}({\bf
x}_{i}\cdot{\bf x}_{j})\Big{)},\end{split}$ (5)

where in the last two terms, the vector triple product has been invoked. The
operation of VISA, therefore, relies on the gradient descent of a gain-driven
losses landscape.

We illustrate a toy model of VISA for two dimensions in Fig. (1),
demonstrating how the energy landscape evolves as the gain $\gamma_{i}(t)$ and
penalty $P(t)$ are evolved, and how the minimizers ${\bf x}^{*}$ and minima of
$H_{\rm VISA}$ change during this process.

Figure 1: Energy landscape of the two-dimensional VISA Hamiltonian for $N=2$
with spin-spin coupling terms $J_{11}=J_{22}=0$ and $J_{12}=J_{21}=-1$. The
normalized minimizer is ${\bf x}_{1}=(1,-1)/\sqrt{2}$ and ${\bf
x}_{2}=(-1,1)/\sqrt{2}$, and therefore in this illustration we choose state
vectors ${\bf x}_{1}=(x,y)$ while setting ${\bf x}_{2}=(-1,1)/\sqrt{2}$. The
self-interacting term $H_{1}$ is shown in (a)-(c) for $\gamma_{i}=\gamma=-0.5$
(a), $0.5$ (b), and $1$ (c), illustrating the symmetry-breaking as the gain
increases. $H_{1}+H_{2}$ is shown in (d)-(f) which seeks to minimise
interaction energy between soft spins; in (f) the minimum is achieved as the
the angle between vectors ${\bf x}_{1}$ and ${\bf x}_{2}$ is $\pi$; in a more
general coupling matrix the vector spins will tend to adjust the angles as to
minimise the pairwise interactions. $H_{1}+H_{2}+H_{3}$ is shown in (g)-(i)
for the same gain parameter as in (a)-(c) and $P(t)=$ $1.0$ (g), $1.5$ (h),
and $2.0$ (i). The vector spins are forced to be collinear with the direction
that can be identified with the hard Ising spins in (i).

##  III Numerical Evolution of Vector Ising Spin Annealers

We consider two popular models to introduce the concept of a scalar-based GBC
that VISA will be benchmarked against. (A) Hopfield-Tank (HT) neural networks
[31, 32] have energy landscape given by the Lyapunov function

$E_{H}=-p(t)\sum_{i=1}^{N}\int_{0}^{v_{i}}g^{-1}(x)\mathrm{d}x-\frac{1}{2}\sum_{i,j=1}^{N}J_{ij}v_{i}v_{j},$
(6)

with nonlinear activation function $v_{i}=g(x_{i})$ and real soft-spin
variables $x_{i}(t)$ describing the network state. At any time $t$ the Ising
state can be obtained from $x_{i}$ by associating spins $s_{i}$ with the sign
of the soft-spin variables $s_{i}=\mathrm{sgn}(x_{i})$. The governing equation
for HT neural networks is $\dot{x}_{i}=p(t)x_{i}+\sum_{j}J_{ij}v_{j}$ with
annealing parameter $p(t)$. This first-order equation can be momentum-enhanced
and replaced with the second order equation leading to Microsoft analogue
iterative machine [19] or Toshiba bifurcation machine [33]. We will use these
enhancements below. (B) CIM using the degenerate optical parametric
oscillators (DOPO) has an energy function

$E_{\rm
CIM}=\frac{1}{4}\sum_{i=1}^{N}(p(t)-x_{i}^{2})^{2}-\frac{1}{2}\sum_{i,j=1}^{N}J_{ij}x_{i}x_{j},$
(7)

where $x_{i}$, $p(t)$, and $J_{ij}$ represent DOPO quadrature, effective laser
pumping power, and coupling strength, respectively. The system evolves as

$\dot{x}_{i}=p(t)x_{i}-x_{i}^{3}+\sum_{j}J_{ij}x_{j}.$ (8)

In the scalar-based GBC models, as the gain $p(t)$ increases from negative
values (representing effective losses) to large positive values (large
effective gain), the amplitudes undergo Hopf bifurcation and reach
$x_{i}=\pm\sqrt{p(T)}$ as $t\rightarrow T,$ $p(t)=p(T)={\rm const}$ for $t>T$.
At the fixed point, the second term in Eq. (6) dominates, which is the target
Ising Hamiltonian scaled by $p(T)$.

The energy landscapes during amplitude bifurcation have many local minima from
the increased degrees of freedom of the soft-mode systems. Moreover, the
ground state of the soft-mode system may correspond to the excited state of
the hard Ising model. At the bifurcation, the system trajectories may get
trapped at these minima and could not transition to the correct global minimum
at higher gains especially when a global cluster spin flips are required [34].
By enhancing the dimensionality of the energy landscape, VISA may allow the
system to evolve from the global minimum of SSIM to the global minimum of HISM
as we now illustrate.

###  III.1 $J$-Möbius Ladder Graph

First, we examine the Ising Hamiltonian minimization on simple cyclic graphs.
These graphs offer analytically solvable benchmarks with distinct and
identifiable obstacles in finding ground states. The process of finding Ising
ground states on such graphs is significantly influenced by the eigenvalues of
the coupling matrix $\mathbf{J}$, especially the relationship between the
eigenvectors’ component signs and the Ising Hamiltonian’s global minimum [34].
For instance, HT networks adjust spin amplitudes to favour the principle
eigenvector component signs [35].

We start by using a Möbius ladder-type graph, a cyclic structure with an even
number of vertices arranged in a ring, featuring variable couplings between
adjacent and opposite vertices. To explore non-trivial ground states, we
introduce equal antiferromagnetic couplings ($J_{ij}=-1$) among nearest
neighbors and variable cross-ring antiferromagnetic couplings ($J_{ij}=-J$).
We refer to the Möbius graphs with these couplings as $J$-Möbius graph.
Following Ref. [34], we define $S_{0}$ as the alternating spin state around
the ring and $S_{1}$ as the state with spins alternating except at two
opposite ring points with frustrated spins, as depicted in Fig. (5)(a) and
(b). When $N/2$ is even, the energies and principal eigenvalues of $S_{0}$ and
$S_{1}$ intersect at $J_{\rm crit}$ and $J_{\rm e}$, respectively. In the
range $J_{\rm e}<J<J_{\rm crit}$, $S_{0}$’s largest eigenvalues are smaller
than those of $S_{1}$, with an eigenvalue gap $\Delta=2\cos(2\pi/N)+2J-2$,
even though $S_{0}$ represents the ground state with lower energy. The
dynamics of semi-classical soft-spin models, including these considerations,
are juxtaposed with quantum annealing approaches in the minimization of the
Ising Hamiltonian on $J$-Möbius ladder graphs, as detailed in Ref. [34].

The presence of amplitude heterogeneities enables the soft-spin model to
acquire and maintain its state achieved at the bifurcation even when it
diverges from the HSIH ground state, potentially complicating the optimization
process when $J_{\rm e}<J<J_{\rm crit}$. To mitigate this issue, the manifold
reduction CIM (MR-CIM) technique was developed, incorporating an additional
feedback mechanism to regulate soft spin amplitudes, ensuring they remain
close to their mean value [34]. The amplitude adjustment after each update at
time step $t$ is governed by

$x_{i}\rightarrow(1-\delta)x_{i}+\delta\frac{Rx_{i}}{|x_{i}|},$ (9)

where $0<\delta<1$. This adjustment draws the spins nearer to the mean, with
the average defined by the squared radius of the quadrature
$R(\mathbf{x})=\sum_{i}x_{i}^{2}/N$, thus aligning the local and global minima
of the soft and hard spin models more closely. Comparative analyses between
various CIM modes and quantum annealing have been conducted, demonstrating
that, despite quantum annealing’s ability to utilize quantum entanglement and
inter-spin correlations to identify ground states, it exhibits heightened
sensitivity to diminishing energy gaps near $J_{\rm crit}$. Consequently,
quantum annealing demands longer annealing schedules to accurately determine
ground state solutions as $J$ nears $J_{\rm crit}$ [34].

Figure 2: Evolution of the VISA model in an $N=8$ $J$-Möbius ladder network
with coupling strength $J=0.4$, as governed by Eq. (5). This figure
illustrates the dynamics of VISA’s soft-spin components: $x_{1}=x_{1i}$ (a),
$x_{2}=x_{2i}$ (b), and $x_{3}=x_{3i}$ (c), with vertical dashed lines
indicating the snapshots at times $t_{1}$ and $t_{2}$. Panels (d) and (e)
depict the corresponding Bloch spheres, highlighting the orientation of spin
vectors at $t_{1}$–influenced by $H_{1}+H_{2}$–and at $t_{2}$, where spins
align in a unified direction due to a high $P(t_{2})$, revealing the HSIM’s
global minimum. The final orientation of the vectors at $t_{2}$ is
spontaneous; thus, the Bloch sphere is adjusted such that the $z$-axis aligns
with this direction. Panel (f) shows the trajectory of the VISA Hamiltonian,
$H_{\rm VISA}$, over time as it converges to the ground state. The numerical
integration of Eq. (5) utilized the following parameters $\varepsilon=0.03$,
$\gamma_{i}(0)=-0.5$, $P(t)=t/200$, and initial conditions $x_{ij}(0)$
uniformly distributed within $[-0.01,0.01]$. The fourth-order Runge-Kutta
method with fixed time step $\Delta t=0.1$ was employed for solving Eq. (5).

Figure (2) illustrates the application of VISA to a $J$-Möbius ladder graph
instance characterized by cross-circle couplings where $J_{\rm e}<J<J_{\rm
crit}$. During the process, the continuous spin components experience an
Aharonov-Hopf bifurcation, exploiting the minimal energy barriers inherent in
soft-spin models. We particularly focus on the intermediary time $t_{1}$
following the bifurcation, noting that at this time, the spin amplitudes have
yet to achieve unit magnitude, and the spins are not collinear. By a later
time $t_{2}>t_{1}$, the spin vectors stabilize as the minimization of the VISA
Hamiltonian progresses. These soft-spin vectors then exhibit Ising spin
characteristics, including unit magnitude and (anti-)parallel alignment, while
the coupling term $H_{2}$ guarantees the identification of the target Ising
Hamiltonian’s ground state.

Figure 3: Bloch spheres illustrate the evolution of VISA spins for an $N=8$
$J$-Möbius graph with coupling constant $J=0.4$. The spins reach final states
adhering to the Ising model’s collinear and unit magnitude constraints. The
spontaneous orientation of vectors at time $t_{2}$ leads to an adjustment of
the Bloch sphere, aligning the $z$-axis with the resultant direction.
Parameters used are $\varepsilon=0.04$, $\gamma_{i}(0)=-0.5$,
$P(t)=\frac{t}{25},$ and initial conditions $x_{pi}(0)$ uniformly distributed
within $[-0.01,0.01]$.

In Fig. (3), we depict typical trajectories of VISA vectors on Bloch spheres,
drawing a parallel to similar quantum representations. To facilitate
visualization, we orient the spontaneously aligning spins along the $z$-axis.
With a fixed interaction strength where $J_{\rm e}<J<J_{\rm crit}$ and
considering $N=8$, we demonstrate the bifurcation dynamics originating from
the center. This visualization shows how the three-dimensional nature of VISA
allows spins to traverse paths that connect various minima, ultimately
converging on the ground state $S_{0}$. The final vector states are indicated
by black arrows, while red points mark the terminal positions of spin vectors
at sequential time steps $t$. These trajectories elucidate the process by
which spins increasingly satisfy collinearity and unit magnitude constraints
through external annealing of $P(t)$ and the feedback mechanism on
$\gamma_{i}$, steering the spin magnitudes to the threshold value of $|{\bf
x}_{i}|=1$.

Figure 4: Panel (a) compares the trajectory evolution of VISA (blue) and CIM
(red) on an $N=8$ $J$-Möbius ladder graph with a coupling constant $J=0.4$.
Here, VISA successfully reaches the ground state $S_{0}$, while CIM remains in
an excited state $S_{1}$. Both systems start from the same initial conditions,
with CIM beginning at $x_{i}(0)=a$ and VISA at ${\bf x}_{i}(0)=(a,b,b)$, where
$a$ and $b$ are uniformly chosen from the ranges $[-0.01,0.01]$ and
$[-0.0001,0.0001]$ respectively. Panels (b)-(d) for CIM and (e)-(g) for VISA
show the spin states at three different times $t_{1}$, $t_{2}$, and $t_{3}$,
marked by vertical dashed lines in panel (a), illustrating CIM’s inability to
overcome the energy barrier to reach the ground state, in contrast to VISA’s
effective navigation in three-dimensional space. Panel (h) compares the VISA
Hamiltonian $H_{\rm VISA}$ against the CIM energy $E_{\rm CIM}$ as calculated
from Eq. (7), while panel (i) focuses on the Ising energy, with the ground
state $S_{0}$ indicated by a black dashed line. Lastly, panel (j) outlines the
regions in the $\gamma-P$ space that correspond to different global minima of
$H_{\rm VISA}$ for $N=8$, showing $S_{1}$ in blue and $S_{0}$ in varying
colors (yellow, orange, red) for $J$ values of $0.35$, $0.40$, and $0.45$,
respectively.

In Fig. (4), we compare VISA with CIM under equivalent starting conditions,
highlighting that while CIM cannot find the ground state within the range
$J_{\rm e}<J<J_{\rm crit}$, VISA excels by leveraging the dynamics of spins in
three dimensions to identify the global minimum. Soft-spins contribute to this
success in both models by facilitating the escape from local minima through
reduced energy barriers, a feature not present in hard-spin models. Unlike
CIM, constrained by the principal eigenvalue of ${\bf J}$ to an excited state,
VISA uses its multidimensional advantage to bridge minima unreachable by CIM’s
one-dimensional approach. Additionally, VISA’s mechanism allows for spin
flipping at energy costs lower than that required by CIM, thanks to its
ability to navigate through three-dimensional space to find the most energy-
efficient path to the global minimum. Figures (4)(a)-(i) depict how VISA
effectively minimizes Ising energy, contrasting with CIM’s stagnation in an
excited state. For $J_{\rm e}<J<J_{\rm crit}$, and as the gain $\gamma_{i}$
and penalty $P(t)$ increase, VISA goes towards the ground state $S_{0}$,
aligning with the lowest energy state indicated by the HSIH, as shown in Fig.
(4)(j). Nonetheless, it’s important to note that increasing spin amplitudes
too fast may also inadvertently heighten energy barriers, potentially impeding
state transitions. Further insight into the structure of the VISA energy
landscape and its comparison with the energy landscape of the scalar models
can be gained by analysing the critical points (see Supplementary
Information). VISA energy landscape shows much diminished energy barriers.

We compare VISA to two second-order scalar networks methods, namely (i)
momentum-enhanced Hopfield-Tank (ME-HT), and (ii) spin-vector Langevin. With
rigorous exploration phases of hyper-parameters, ME-HT outperforms parallel
tempering, simulated annealing, and commercial solver Gurobi at various QUBO
benchmarks [19]. The ME-HT governing equation is

$m\ddot{x}_{i}+\gamma\dot{x}_{i}-\beta(t)x_{i}-\alpha\frac{\partial
H(g(\mathbf{x}))}{\partial x_{i}}=0,$ (10)

with effective mass $m$, momentum term $\gamma\in[0,1)$, and hyperparameters
$\alpha$ and $\beta(t)$, the latter of which undergoes an annealing protocol
given by $\beta(t)=\beta_{0}(1-t/T)$. Equation (10) combines gradient descent
with annealed non-conservative dissipation, and the addition of momentum
distinguishes it from regular first-order HT networks with energy landscape
Eq. (6). Momentum, or the heavy-ball method, aims to overcome the pitfalls of
pathological curvature in deep learning and accelerates standard gradient
descender optimizers [36, 37].

VISA can be further compared and contrasted with the spin-vector Langevin
(SVL) model that was proposed as a classical analogue of quantum annealing
description using stochastic Langevin time evolution governed by the
fluctuation-dissipation theorem [38]. SVL is based on the time-dependent
Hamiltonian used in quantum annealing $H(t)=A(t)H_{0}+B(t)H_{P}$, where
initial Hamiltonian $H_{0}=\sum_{i}\sigma_{i}^{x}$, and problem Hamiltonian
$H_{P}=-\sum_{i,j}J_{ij}\sigma_{i}^{z}\sigma_{j}^{z}$, with Pauli operator
${\bf\sigma}_{i}$ acting on the $i$-th variable. Real annealing functions
satisfy boundary conditions $A(0)=B(T)=1$ and $A(T)=B(0)=0$, where $T$ is the
temporal length of the annealing schedule. If the rate of change of the
functions is slow enough, the system stays in the ground state of the
instantaneous Hamiltonian so that at $t=T$ the Ising Hamiltonian is minimized.
Quantum annealing has shown competitive results in quadratic unconstrained
binary optimization (QUBO) problems such as subset sum, vertex cover, graph
coloring, and travelling salesperson [39]. The SVL model replaces Pauli
operators with real functions of continuous angle
$\sigma_{i}^{z}\rightarrow\sin\theta_{i}$,
$\sigma_{i}^{x}\rightarrow\cos\theta_{i}$, and is therefore a classical
annealing Hamiltonian using continuous-valued spins $\sin\theta_{i}$. SVL
dynamics is described by a system of coupled stochastic equations

$m\ddot{\theta}_{i}+\gamma\dot{\theta}_{i}+\frac{\partial
H(\mathbf{\theta})}{\partial\theta_{i}}+\xi_{i}(t)=0,$ (11)

where $m$ is the effective mass, $\gamma$ is the damping constant, and
$\xi_{i}(t)$ is an iid Gaussian noise. For long annealing times, the minima of
HSIHs are obtained through the transformation
$s_{i}=\mathrm{sgn}(\sin\theta_{i})$. The gradient term in Eq. (11) is

$\frac{\partial
H(\mathbf{\theta})}{\partial\theta_{i}}=-B(t)\sum_{j=1}^{N}J_{ij}\cos\theta_{i}\sin\theta_{j}+A(t)\sin\theta_{i},$
(12)

which in conjunction with fluctuation-dissipation relations
$\langle\xi_{i}(t)\rangle=0$ and
$\langle\xi_{i}(t)\xi_{j}(t^{\prime})\rangle=\delta_{ij}\delta(t-t^{\prime})$,
give $2N$ stochastic differential equations:
$\mathrm{d}\theta_{i}=(p_{i}/m)\mathrm{d}t$ and

$\mathrm{d}p_{i}=\Big{(}\frac{\partial
H(\mathbf{\theta})}{\partial\theta_{i}}+\frac{\gamma}{m}p_{i}\Big{)}\mathrm{d}t+\mathrm{d}W_{i},$
(13)

where $\mathrm{d}W_{i}$ represents a real-valued continuous-time stochastic
Wiener process [38]. The characteristic amplitude bifurcation of scalar spins
according to ME-HT and SVL are given in Supplementary Information.

The key feature distinguishing VISA from the discussed models lies in its
novel gain-based and dimensionality annealing strategy, applied across
multiple dimensions. In the analysis of $J$-Möbius ladder graphs, as shown in
Fig. (5), we compute the ground state probability for VISA alongside SVL, ME-
HT, CIM, and the Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm [40].
Within the range $J_{\rm e}<J<J_{\rm crit}$, VISA consistently identifies the
ground state $S_{0}$ with a higher probability $p_{\rm GS}$ compared to the
other models.

Figure 5: $N=8$ $J$-Möbius ladder graphs with circle (blue) and cross-circle
(red) couplings for the (a) $S_{0}$ and (b) $S_{1}$ configurations. (c) Ground
state probability for VISA, SVL, ME-HT, CIM, and BFGS for an $N=8$ $J$-Möbius
ladder graph with varying cross-circle couplings $J$. One thousand runs are
used to calculate the probability of finding the ground state $p_{\rm GS}$ for
each value of $J$. Dashed lines corresponding to $J_{\rm e}$ (black) and
$J_{\rm crit}$ (red) show the values of $J$ for which the leading eigenvalue
changes and the ground state configuration changes, respectively. For VISA,
$\varepsilon=0.02$, $\gamma_{i}(0)=-0.5$, and $P(t)=\alpha t$. For SVL,
$m=1.0$, $\xi\sim N(\mu=0,\sigma=0.1)$, and $\gamma=0.99$. For momentum-
enhanced HT, $m=1.0$, and $\gamma=0.99$. For CIM we used
$p(t)=(1-p_{0})\tanh(\varepsilon t)+p_{0}$ with $p_{0}=J-2$, and
$\varepsilon=0.003$. For each value of $J$, optimal values of $\alpha$ and
$\beta(0)$ are chosen based on sets of preliminary runs in which they are
varied. Inset (d) illustrates the structure of Möbius ladder coupling matrices
${\bf J}$ with $J_{ij}=0$ (white), $-1$ (blue), $-J$ (red).

###  III.2 $J-G$ Cyclic Graphs

Next we consider the $J-G$ cyclic graphs that are a variant of $J$-Möbius
ladder graphs that include additional couplings, connecting each vertex $i$ to
vertices $i\pm k$ with a weight of $-G$, where $1<k<N/2$. This modification
aims to explore more complex interaction patterns and their impact on the
system’s ground state behavior and eigenvalue distributions, particularly
focusing on how these factors evolve in the context of optimization and the
search for ground states in varied graph structures. The interaction strength
range is expanded to accommodate both ferromagnetic and antiferromagnetic
interactions, with $-1\leq J,G\leq 1$. This yields a 5-regular circulant graph
characterized by weights $J_{ij}\in\\{-1,-J,-G\\}$, depicted in Fig. (6).
Cyclic graphs maintain their local and global topological properties under
rotational transformations, encapsulating all connectivity information within
any row of ${\bf J}$. By selecting the first row $J_{1,j}$, we compute
eigenvalues as $\lambda_{n}=\sum_{j=1}^{N}J_{1,j}\cos\left[2\pi
n(j-1)/N\right]$, leading to $\lambda_{n}=-2\cos(2\pi
n/N)-J(-1)^{n}-2G\cos(2\pi kn/N)$ [41]. We then deduce boundaries by observing
eigenvalues in the $J-G$ plane, identifying where leading eigenvalues and
corresponding eigenvectors change their correspondence with the ground and
excited states. The ground state boundaries are identified within the $J-G$
space in the Supplementary Information. Figure (6) depicts regions where the
global minimum diverges from the hypercube corner of the projected eigenvector
with the highest eigenvalue.

Figure 6: $J-G$ cyclic graphs for $N=8$ with couplings (a) $k=2$, and (b)
$k=3$. Panel (c) depicts the matrix structure ${\bf J}$ for $k=3$, where the
color coding represents $J_{ij}=0$ (white), $-1$ (blue), $-J$ (red), and $-G$
(green). The ground states and leading eigenvalue states for $k=2$ are
displayed in panels (d) and (e) respectively, within the $J-G$ space. Panel
(f) delineates the energy and eigenvalue boundaries, with the pink region
indicating discrepancies between ground and maximal eigenvalue states. The
sequence from panels (g) to (i) extends this analysis to $k=3$, highlighting
the ground state probabilities across easy-hard transitions in the $J-G$
space, indicated by solid red lines in panels (f) and (i).

We investigate ground state probabilities $p_{\rm GS}$ over transition regions
in these graphs with both $J$ and $G$ cross-circle couplings. We choose values
in $J-G$ parameter space that demonstrate transitions between regions in which
the eigenvector corresponding to the leading eigenvalue does not match the
ground state solution (exact calculations of these regions are presented in
Supp. Inf.). Ground state probabilities for VISA, SVL, ME-HT, CIM, and BFGS
are illustrated in Fig. (7). Four cases are analyzed, split into two sets:
$k=2$ and $k=3$. We consider perpendicular lines in the two-dimensional $J-G$
space for each set. Specifically, we vary (fix) $J$ and fix (vary) $G$. For
$k=2$ and $G=0.5$, an easy-hard-easy transition emerges as $J$ increases, akin
to $J$-Möbius ladder previously studied. Indeed, for $J_{\rm e}<J<J_{\rm
crit}$, $p_{GS}$ decreases as the eigenvalue gap $\Delta$ increases. If,
instead, we fix $J=-0.5$, a transition occurs, centred at the change between
ground states given by $G_{\rm crit}=0.5$. Here, the hard region is bounded by
eigenvalue crossing points $G_{\rm e_{1}}=1-\sqrt{2}/2$ and $G_{\rm
e_{2}}=1/\sqrt{2}$. For $k=3$, $p_{\rm GS}$ is less sensitive to the magnitude
of $\Delta$ for regions where $S_{3}$ is the ground state. This is due to the
proximity between the leading eigenvalue state $S_{1}$ and ground solution
$S_{3}$ in topological spin space. More precisely, the transformation from
$S_{1}$ to $S_{3}$ requires only a single spin flip, representing a nominal
energy barrier for soft-spin models. Therefore, hardness in $J-G$ cyclic
graphs derives not just from the eigenvalue gap magnitude but additionally
from a distance between hypercube corners of the ground and leading eigenvalue
states.

Figure 7: Ground state probability for VISA, SVL, ME-HT, CIM, and BFGS for
$N=8$ $J-G$ cyclic graphs with varying cross-circle couplings $J$ and $G$. (a)
$k=2$, $J\in[-0.5,0.5]$, and $G=0.5$. (b) $k=2$, $J=0.5$, and $G\in[0,1]$. (c)
$k=3$, $J\in[-1,0.5]$, and $G=-0.8$. (d) $k=3$, $J=-0.5$, and $G\in[-1,0]$.
One thousand runs are used to calculate the probability of finding the ground
state $p_{\rm GS}$ for each value of $J$. Dashed red and black lines show
where the ground state energies and leading eigenvalues change, respectively.
For VISA, $\varepsilon=0.04$, $\gamma_{i}(0)=-0.5$, and $P(t)=\alpha t$. For
SVL, $m=1.0$, $\xi\sim N(\mu=0,\sigma=0.1)$, and $\gamma=0.99$. For momentum-
enhanced HT, $m=1.0$, and $\gamma=0.99$. For CIM, the general annealing
protocol is used with $p_{0}=J-2$, and $\varepsilon=0.003$. For each value of
$J$, optimal values of $\alpha$ and $\beta(0)$ are chosen based on sets of
preliminary runs in which they are varied.

###  III.3 Random Graphs

We extend VISA’s evaluation to QUBO instances renowned for their computational
intensity as they scale: dense fully connected graphs and sparse three-regular
graphs with normally distributed random couplings $J_{ij}$, representing
Sherrington-Kirkpatrick (SK) and weighted three-regular Max-Cut problems,
respectively. These models are pivotal for benchmarking physical simulators
[42, 43, 44, 45] and are categorized under the NP-hard complexity class [46,
4]. Unlike cyclic graphs, these instances do not have analytically known
ground states, necessitating the use of the Gurobi optimization suite for
ground state estimations. Gurobi employs advanced pre-processing and heuristic
enhancements to expedite branch-and-bound algorithms [47]. Figure (8)(a)
compares the ground state approximations achieved by Gurobi with those by VISA
for $N=100$ SK and weighted three-regular Max-Cut graphs. Alongside, we
include comparisons with SVL, MR-CIM, and CIM, where MR-CIM and CIM adapt
their laser intensities following a general pumping scheme
$p(t)=(1-p_{0})\tanh(\varepsilon t)+p_{0}$. MR-CIM further applies an
additional feedback mechanism as per Eq. (9) on top of Eq. (8), controlling
soft-spin amplitudes and the dimensionality landscape. We define the quality
improvement of VISA over another method $X$ in terms of objective values $O$
as $(O_{\rm VISA}-O_{X})/O_{\rm VISA}$, showcasing these metrics for SK and
three-regular problems in Fig. (8)(b), where $X$ represents the best-
performing competing method for each instance.

Figure 8: (a) The proximity gap, defined as the ratio of the found objective
to the Gurobi objective, for VISA, SVL, manifold reduction CIM, and CIM
methods. Gurobi was given a 30-second time limit for each of the 100 problem
instances with size $N=100$. The instances are divided equally between two
graph topologies: dense, fully connected, and sparse three-regular graphs. In
both cases, the matrix weight elements are drawn from the Gaussian
distribution with zero mean and unit variance, resulting in instances
belonging to the Sherrington-Kirkpatrick and weighted 3-regular Max-Cut
problems. (b) Violin plots demonstrate the distribution of VISA’s quality
improvement performance compared to the best solution found by competing
methods {SVL, MR-CIM, CIM} across the SK and 3-Regular QUBO benchmarks.

##  IV Conclusions

This paper introduces the Vector Ising Spin Annealer (VISA), a model that
capitalizes on the advantages of multidimensional spin systems, gain-based
operation and soft-spin annealing techniques to optimize Ising Hamiltonians on
a variety of graph structures. VISA distinguishes itself by enabling more
efficient navigation through the solution space, enhancing spin mobility in
higher-dimensional spaces, and providing a robust framework for connecting
local minima and reducing energy barriers.

A key focus of VISA is its ability to recover ground states effectively, even
in scenarios where these states do not correspond to the principal eigenvector
of the coupling matrix. The model’s performance was numerically evaluated
against other methods, demonstrating its superior ability to find ground
states across different graph types and complex QUBO instances, thereby
highlighting its potential to address NP-hard problems.

The future research could explore the role of defects in spin models, such as
topological defects, domain walls, and vortex rings and their role in
achieving the ground state during gain-based operation. In higher-dimensional
systems, such as those utilized by VISA, vortices may exhibit more efficient
annihilation properties. This is potentially due to the additional spatial
degrees of freedom, which could facilitate the merging or cancellation of
vortices and anti-vortices, a phenomenon less constrained than in two-
dimensional spaces. This enhanced annihilation could lead to a smoother energy
landscape, aiding the system in avoiding local minima and more effectively
converging to the ground state.

By leveraging multidimensional spins, VISA opens new avenues for the
development of analogue optimization machines, potentially incorporating
quantum effects to enhance computational capabilities. The prospects of
applying dimensionality annealing techniques in optical-based Ising machines
suggest an exciting future for speed-of-light computation and accurate ground
state recovery, marking a significant advancement in the field of optimization
technologies.

##  V Acknowledgements

J.S.C. acknowledges the PhD support from the EPSRC, N.G.B. acknowledges the
support from the Julian Schwinger Foundation Grant No. JSF-19-02-0005, the
HORIZON EIC-2022-PATHFINDERCHALLENGES-01 HEISINGBERG project 101114978, and
Weizmann-UK Make Connection grant 142568.

##  VI Supplementary Information

###  VI.1 Critical Points of VISA Hamiltonian

Figure 9: Critical points of $H_{\rm VISA}$ for $N=8$ $J$-Möbius ladder graph
with $J=0.4$. Each point is characterized by its energy $E$ and Euclidean
distance from the origin $d$. The number of soft spins normalizes distance
$d$. Blue and red circles denote minima for states $S_{1}$ and $S_{0}$
respectively, while black circles indicate other minima and saddle points of
$H_{\rm VISA}$. The closest saddle point to excited state $S_{1}$ is
highlighted in green. Progressing from figures (a) to (f), we systematically
increase the gain $\gamma_{i}=\gamma$ for all $i$ and the penalty $P$,
following a typical annealing schedule. Notably, at (c), with $\gamma=-0.087$
and $P=0.32$, we reach a critical juncture that delineates the transition
between global minima $S_{0}$ and $S_{1}$, as illustrated in Fig. (4)(j),
where $S_{0}$ and $S_{1}$ attain equal energy levels $E$. Adjacent bars to the
main plots list distance $\Delta d$ for journey $S_{1}\rightarrow{\rm
SP}\rightarrow S_{0}$ for VISA (blue) and CIM (red), where for CIM we take
$p=\gamma$.

To determine the critical points of $H_{\rm VISA}$, we set $\frac{\partial
E}{\partial x_{i}}=0$ for all $i$, where $i$ ranges from $1$ to $N$. Figure
(9) depicts these critical points, with the minima corresponding to the states
$S_{0}$ and $S_{1}$ emphasized. As the gain $\gamma_{i}(t)$ and the
collinearity penalty $P(t)$ are systematically increased, a transition in the
ground state from $S_{1}$ to $S_{0}$ is observed. Notably, there is an
infinite continuum of $(\gamma_{i},P)$ pairs where $S_{0}$ and $S_{1}$ share
the same energy level. This continuum forms a demarcation line in the
$\gamma-P$ space, illustrated in Fig. (4)(j). At a specific point on this
boundary, as depicted in Fig. (9)(c), both $S_{0}$ and $S_{1}$ emerge as
ground states. Furthermore, Fig. (9) quantifies the average Euclidean distance
$\Delta d$ that each soft-spin $x_{pi}$ needs to traverse to transition from
the local minimum $S_{1}$ to the global minimum $S_{0}$, passing through the
nearest saddle point. In every scenario analyzed, the VISA model demonstrates
a shorter requisite distance compared to the analogous CIM model, underscoring
its efficiency in navigating the energy landscape.

###  VI.2 Basins of Attraction

Figure 10: Basins of attraction of the VISA model landscape on $J$-Möbius
ladder graph. We take $N=8$ and $J=0.4$ with various gains
$\gamma_{i}=\gamma\forall i$ and penalties $P$ as they anneal in the VISA
method. 4000 randomly distributed starting points $x_{pi}(0)$ in $[-1,1]$ are
plotted with colors corresponding to minima reached via gradient descent using
Newton’s method. $(\gamma,P)$ (a) $(-0.5,0)$, (b) $(-0.25,0.15)$, (c)
$(-0.087,0.32)$, (d) $(0.25,0.5)$, (e) $(0.5,0.75)$, and (f) $(0.75,1.0)$. To
characterize points, the magnetization magnitude $|{\bf
m}|=\sqrt{m_{x}^{2}+m_{y}^{2}+m_{z}^{2}}$ and correlation magnitude $|{\bf
X_{\rm corr}}|=\sqrt{X_{x}^{2}+X_{y}^{2}+X_{z}^{2}}$ are used, where
$m_{p}=\sum_{i}x_{pi}/N$ and
$X_{p}=\sum_{i}(x_{pi}-m_{p})(x_{p(i+1)}-m_{p})/\sum_{i}(x_{pi}-m_{p})^{2}$
are the average magnetization and correlations along the $p$-axis,
respectively. For small $\gamma$ and $P$, the basin of attraction is dominated
by the $S_{1}$ state. As $\gamma$ and $P$ grow, the volume of the basin of
attraction for state $S_{0}$ increases. For $\gamma$ and $P$ which lie on the
critical boundary, shown here in (c), both $S_{0}$ and $S_{1}$ states have the
same energy. As $\gamma$ and $P$ increase further, $S_{0}$ becomes the ground
state. This is indicated by the switch between excited (red) and ground (blue)
states.

Figure (10) illustrates the basins of attraction for a range of gain
$\gamma_{i}=\gamma\forall i$ and penalty $P$ values on $J$-Möbius ladder
graphs with $J=0.4$. Panels (a) to (f) in the figure show the evolution of
these basins as $\gamma$ and $P$ are incrementally increased, following a
typical gain-based and penalty annealing schedule in the VISA model. The
basins of attraction are defined by the initial states $x_{pi}(0)$, uniformly
distributed across the interval $[-1,1]$, which converge to distinct minima
through gradient descent. Initially, at $(\gamma,P)=(-0.5,0)$, the global
minimum is at $S_{1}$, which has the largest basin of attraction, capturing
all initial states towards $S_{1}$. However, as $\gamma$ and $P$ are
progressively increased, the basin of attraction associated with the excited
state $S_{0}$ expands. Beyond the critical threshold in the $\gamma-P$
parameter space, demonstrated in Fig. (4)(j), $S_{0}$ transitions to become
the new ground state. Further annealing of $\gamma$ and $P$ beyond this
threshold reveals additional higher-energy states through the process of
gradient descent.

###  VI.3 Second-Order Scalar Networks Minimizing the Ising Hamiltonian on
$J$-Möbius Ladder Graphs

Figure 11: Evolution of momentum enhanced Hopfield-Tank method for an $N=8$
$J$-Möbius ladder network with $J=1$. (a) The amplitudes connected by the
frustrated edges are lower than in the rest of the system and are shown in
red. (b) The corresponding Ising Hamiltonian decreases in time, and the ground
state is found, with energy given by the dashed red line. (c) Schematic
representation of state $S_{1}$ realised by the soft-spin momentum enhanced HT
model. Here, $m=1.0$, $\alpha=2.5/\lambda_{\rm max}$, $\beta=1.5(1-t/T)$,
$\gamma=0.7$, and $x_{i}(0)$ is chosen uniformly at random from the range
$[-0.01,0.01]$. Fourth order Runge-Kutta is used with fixed time step $\Delta
t=0.1$ to solve Eq. (10).

Figure (11) illustrates the Aharonov-Hopf bifurcation in soft-spins within the
ME-HT model. For each time step $t$, we apply a clipping function $[{\bf
x}_{i}]_{t}=\text{sgn}([{\bf x}_{i}]_{t})$ when $|[{\bf x}_{i}]_{t}|>1$,
alongside a nonlinear activation function $g({\bf x}_{i})=\tanh({\bf x}_{i})$.
The parameter $\alpha$ is adjusted according to the largest positive
eigenvalue $\lambda_{\rm max}$ of the adjacency matrix ${\bf J}$.
Concurrently, Fig. (12) displays the temporal evolution of scalar amplitudes
as dictated by the SVL model, applied to the same Möbius graph with $J=1$.

Figure 12: Evolution of the spin-vector Langevin (SVL) model. Application of
the SVL model on an $N=8$ $J$-Möbius ladder network with coupling strength
$J=1$. (a) depicts the amplitudes connected by frustrated edges, shown in red,
which are lower than those in the rest of the system. (b) illustrates the real
annealing functions $A(t)$ in blue and $B(t)$ in red, satisfying the boundary
conditions $A(0)=B(T)=1$ and $A(T)=B(0)=0$, over an annealing period $T=100$.
(c) shows the total Hamiltonian $H(t)=A(t)H_{0}+B(t)H_{P}$ dynamics: it
decreases initially as the ground state of $H_{0}$ is reached, then increases
as $A(t)$ and $B(t)$ are annealed. Post a critical juncture, the SVL model
aims to minimize the problem Hamiltonian $H_{P}$, identifying the ground state
of the $J$-Möbius ladder, indicated by the dashed red line. Parameters used
are $m=1.0$, $\alpha=1/\lambda_{\rm max}$, $\gamma=1.0$, and noise $\xi=0$.
Annealing schedules are defined as $A(t)=1-t/T$ and $B(t)=t/T$. Initial phases
$\theta_{i}(0)$ are randomly selected from $[-\pi/2,\pi/2]$. Equation (11) is
solved using a fourth-order Runge-Kutta method with a fixed time step of
$\Delta t=0.1$.

###  VI.4 Eigenvalues, Eigenvectors, and Ground States of $J-G$ Cyclic Graphs

Here we refer to $J-G$ cyclic graph defined in the main text, and illustrated
in Figs. (6)(a)-(b) for $N=8$. For general $N$ with couplings
$J_{ij}\in\\{-1,-J,-G\\}$, the $N$ eigenvalues are given by

$\lambda_{n}=-2\cos(2\pi n/N)-J(-1)^{n}-2G\cos(2\pi kn/N),$ (14)

where $n=1,2,\ldots,N$. Equation (14) follows from substituting $J_{1,j}$ into
the general form of matrix eigenvalues for cyclic graphs given in the main
text. As $J$ and $G$ vary within $[-1,1]$, the leading eigenvalue changes.
Constrained to this domain, for $N=8$ and $k=2$, the leading eigenvalues are
given by $n=4,5,6$. Substituting these values into Eq. (14) and equating each
pair gives three boundaries defined by $J=1-\sqrt{2}/2-G$, $J=G-1/\sqrt{2}$,
and $G=1/2$. Similarly for $k=3$, we obtain two unique leading eigenvalues
($n=4,5$) in the $J-G$ domain, with boundary $J=1-\sqrt{2}/2+(1+1/\sqrt{2})G$.
Eigenvectors corresponding to leading eigenvalues are used to deduce analogous
Ising states by projecting the eigenvector onto the nearest hypercube corner
$[-1,1]^{N}$. These Ising states are given in Figs. (6)(e) and (h). Gurobi is
used to obtain ground states in $J-G$ space. Four unique configurations are
found with $N=8$: namely $S_{0}$, $S_{1}$, and $S_{2}$ for $k=2$, and $S_{0}$,
$S_{1}$, and $S_{3}$ for $k=3$. $S_{2}$ and $S_{3}$ are defined as the
configurations of Ising spins given by $S_{2}=(1,1,-1,-1,1,1,-1,-1)$ and
$S_{3}=(1,1,-1,1,1,-1,1,-1)$. For even $N/2$, these four states have energies

$\displaystyle E_{0}(J,G)$ $\displaystyle=(J-2-2G(-1)^{k+1})N/2,$ (15)
$\displaystyle E_{1}(J,G)$ $\displaystyle=4-(J+2)N/2+(-1)^{k}(N-4k)G,$ (16)
$\displaystyle E_{2}(J,G)$
$\displaystyle=((-1)^{N/4}J+(1+(-1)^{k})(-1)^{d}G)N/2,$ (17) $\displaystyle
E_{3}(J,G)$ $\displaystyle=4-N+(2-N/2)J+(-1)^{k}(N-4k)G,$ (18)

where $d=k/2$ if $k$ is even, and $0$ otherwise. Equating the relevant
energies for $k=2$ gives three boundaries defined by separators $J=1/2-G$,
$J=G-1/2$, and $G=1/2$. Similarly for $k=3$, the boundary edges are given by
$J=1/2+3/2G$, $J=2/3+2G$, and $J=0$. Eigenvalues and ground states for $J-G$
cyclic graphs are shown in Fig. (13) for values of $J$ and $G$ considered in
the main text.

Figure 13: Eigenvalues of $N=8$ $J-G$ cyclic graphs as functions of $J$ and
$G$. Each plot corresponds to one of the red lines in Figs. (6)(f) and (i).
Ground state probabilities are calculated on these lines in $J-G$ space for
various methods in the main text. (a) $k=2$ and $G=0.5$. (b) $k=2$ and $J=0$.
(c) $k=3$ and $G=-0.8$. (d) $k=3$ and $J=-0.5$. Vertical dashed red and black
lines show where ground state energies and leading eigenvalues change,
respectively. Background colors denote ground states in specific regions, with
states and corresponding colors listed in the key.

## References

  * Thompson and Spanuth [2021] N. C. Thompson and S. Spanuth, Communications of the ACM 64, 64 (2021).
  * Vadlamani  _et al._ [2020] S. K. Vadlamani, T. P. Xiao, and E. Yablonovitch, Proceedings of the National Academy of Sciences 117, 26639 (2020).
  * Farhi  _et al._ [2000] E. Farhi, J. Goldstone, S. Gutmann, and M. Sipser, arXiv preprint quant-ph/0001106 (2000).
  * Lucas [2014] A. Lucas, Frontiers in physics 2, 5 (2014).
  * 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, Nature materials (2017).
  * Kalinin  _et al._ [2020] K. P. Kalinin, A. Amo, J. Bloch, and N. G. Berloff, Nanophotonics 9, 4127 (2020).
  * De las Cuevas and Cubitt [2016] G. De las Cuevas and T. S. Cubitt, Science 351, 1180 (2016).
  * 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._ , Science 354, 614 (2016).
  * Inagaki  _et al._ [2016] T. Inagaki, Y. Haribara, K. Igarashi, T. Sonobe, S. Tamate, T. Honjo, A. Marandi, P. L. McMahon, T. Umeki, K. Enbutsu, _et al._ , Science 354, 603 (2016).
  * Yamamoto  _et al._ [2017] Y. Yamamoto, K. Aihara, T. Leleu, K.-i. Kawarabayashi, S. Kako, M. Fejer, K. Inoue, and H. Takesue, npj Quantum Information 3, 1 (2017).
  * Honjo  _et al._ [2021] T. Honjo, T. Sonobe, K. Inaba, T. Inagaki, T. Ikuta, Y. Yamada, T. Kazama, K. Enbutsu, T. Umeki, R. Kasahara, _et al._ , Science advances 7, eabh0952 (2021).
  * Cai  _et al._ [2020] F. Cai, S. Kumar, T. Van Vaerenbergh, X. Sheng, R. Liu, C. Li, Z. Liu, M. Foltin, S. Yu, Q. Xia, _et al._ , Nature Electronics 3, 409 (2020).
  * Babaeian  _et al._ [2019] M. Babaeian, D. T. Nguyen, V. Demir, M. Akbulut, P.-A. Blanche, Y. Kaneda, S. Guha, M. A. Neifeld, and N. Peyghambarian, Nature communications 10, 1 (2019).
  * Pal  _et al._ [2020] V. Pal, S. Mahler, C. Tradonsky, A. A. Friesem, and N. Davidson, Physical Review Research 2, 033008 (2020).
  * Parto  _et al._ [2020] M. Parto, W. Hayenga, A. Marandi, D. N. Christodoulides, and M. Khajavikhan, Nature materials 19, 725 (2020).
  * Pierangeli  _et al._ [2019] D. Pierangeli, G. Marcucci, and C. Conti, Physical review letters 122, 213902 (2019).
  * 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._ , Nature communications 11, 249 (2020).
  * Vretenar  _et al._ [2021] M. Vretenar, B. Kassenberg, S. Bissesar, C. Toebes, and J. Klaers, Physical Review Research 3, 023167 (2021).
  * Kalinin  _et al._ [2023] K. Kalinin, G. Mourgias-Alexandris, H. Ballani, N. G. Berloff, J. H. Clegg, D. Cletheroe, C. Gkantsidis, I. Haller, V. Lyutsarev, F. Parmigiani, L. Pickup, _et al._ , arXiv preprint arXiv:2304.12594 (2023).
  * Goto  _et al._ [2021] H. Goto, K. Endo, M. Suzuki, Y. Sakai, T. Kanao, Y. Hamakawa, R. Hidaka, M. Yamasaki, and K. Tatsumura, Science Advances 7, eabe7953 (2021).
  * Date  _et al._ [2021] P. Date, D. Arthur, and L. Pusey-Nazzaro, Scientific reports 11, 10029 (2021).
  * Gilli  _et al._ [2019] M. Gilli, D. Maringer, and E. Schumann,  _Numerical methods and optimization in finance_ (Academic Press, 2019).
  * Pierce and Winfree [2002] N. A. Pierce and E. Winfree, Protein engineering 15, 779 (2002).
  * Dill  _et al._ [2008] K. A. Dill, S. B. Ozkan, M. S. Shell, and T. R. Weikl, Annu. Rev. Biophys. 37, 289 (2008).
  * Calvanese Strinati and Conti [2022] M. Calvanese Strinati and C. Conti, Nature Communications 13, 7248 (2022).
  * Strinati and Conti [2024] M. C. Strinati and C. Conti, Physical Review Letters 132, 017301 (2024).
  * Kalinin and Berloff [2018] K. P. Kalinin and N. G. Berloff, New Journal of Physics 20, 113023 (2018).
  * Stroev and Berloff [2021] N. Stroev and N. G. Berloff, Physical Review Letters 126, 050504 (2021).
  * Rosenberg [1975] I. G. Rosenberg, (1975).
  * Boros and Hammer [2002] E. Boros and P. L. Hammer, Discrete applied mathematics 123, 155 (2002).
  * Hopfield [1982] J. J. Hopfield, Proceedings of the national academy of sciences 79, 2554 (1982).
  * Hopfield and Tank [1985] J. J. Hopfield and D. W. Tank, Biological cybernetics 52, 141 (1985).
  * Goto [2016] H. Goto, Scientific reports 6, 1 (2016).
  * Cummins  _et al._ [2023] J. S. Cummins, H. Salman, and N. G. Berloff, arXiv preprint arXiv:2311.17359 (2023).
  * Kalinin and Berloff [2022] K. P. Kalinin and N. G. Berloff, Communications Physics 5, 20 (2022).
  * Orvieto and Lucchi [2019] A. Orvieto and A. Lucchi, Advances in Neural Information Processing Systems 32 (2019).
  * Saab Jr  _et al._ [2022] S. Saab Jr, S. Phoha, M. Zhu, and A. Ray, Machine Learning 111, 3245 (2022).
  * Subires  _et al._ [2022] D. Subires, F. J. Gómez-Ruiz, A. Ruiz-García, D. Alonso, and A. Del Campo, Physical Review Research 4, 023104 (2022).
  * Jiang and Chu [2022] J.-R. Jiang and C.-W. Chu, in  _2022 IEEE 4th ECICE_ (IEEE, 2022) pp. 406–411.
  * [40] The BFGS algorithm is implemented with minimizer FindMinimum in Wolfram Mathematica. VISA Hamiltonian $H_{\rm VISA}$ is subject to BFGS with fixed $\gamma_{i}=\gamma=1.0\forall i$ and $P=0.8$.
  * Gancio and Rubido [2022] J. Gancio and N. Rubido, Chaos, Solitons & Fractals 158, 112001 (2022).
  * Haribara  _et al._ [2017] Y. Haribara, H. Ishikawa, S. Utsunomiya, K. Aihara, and Y. Yamamoto, Quantum Science and Technology 2, 044002 (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._ , Science advances 5, eaau0823 (2019).
  * Harrigan  _et al._ [2021] M. P. Harrigan, K. J. Sung, M. Neeley, K. J. Satzinger, F. Arute, K. Arya, J. Atalaya, J. C. Bardin, R. Barends, S. Boixo, _et al._ , Nature Physics 17, 332 (2021).
  * Böhm  _et al._ [2019] F. Böhm, G. Verschaffelt, and G. Van der Sande, Nature communications 10, 1 (2019).
  * Arora  _et al._ [2005] S. Arora, E. Berger, H. Elad, G. Kindler, and M. Safra, in  _46th Annual IEEE FOCS Symposium_ (IEEE, 2005) pp. 206–215.
  * Gurobi Optimization, LLC [2023] Gurobi Optimization, LLC, “Gurobi Optimizer Reference Manual,”  (2023).