# Polynomial unconstrained binary optimization inspired by optical simulation

D.A. Chermoshentsev  e-mail<EMAIL_ADDRESS>Russian
Quantum Center, Skolkovo, Moscow 143025, Russia  Moscow Institute of Physics
and Technology, Moscow Region 141701, Russia  Skolkovo Institute of Science
and Technology, Skolkovo, Moscow, 143025, Russia  QBoard, Skolkovo, Moscow
143025, Russia  A.O. Malyshev  Clarendon Laboratory, University of Oxford,
Parks Road, Oxford OX1 3PU, UK  M. Esencan  Clarendon Laboratory, University
of Oxford, Parks Road, Oxford OX1 3PU, UK  E.S. Tiunov  Russian Quantum
Center, Skolkovo, Moscow 143025, Russia  Moscow Institute of Physics and
Technology, Moscow Region 141701, Russia  D. Mendoza  Department of Chemistry,
University of Toronto, 80 St. George Street, Toronto, ON M5S 3H6, Canada
Department of Chemistry and Chemical Biology, Harvard University, 12 Oxford
Street, Cambridge, MA 02138, USA  A. Aspuru-Guzik  Department of Chemistry,
University of Toronto, 80 St. George Street, Toronto, ON M5S 3H6, Canada
Canadian Institute for Advanced Research (CIFAR) Senior Fellow, Toronto, ON
M5S 1M1, Canada  Canada CIFAR AI Chair Lebovic Fellow, Vector Institute,
Toronto, ON M5S 1M1, Canada  Department of Computer Science, University of
Toronto, 40 St. George Street, Toronto, ON M5S 2E4, Canada  A.K. Fedorov
Russian Quantum Center, Skolkovo, Moscow 143025, Russia  QBoard, Skolkovo,
Moscow 143025, Russia  National University of Science and Technology “MISIS”,
Moscow 119049, Russia  A.I. Lvovsky  Russian Quantum Center, Skolkovo, Moscow
143025, Russia  Clarendon Laboratory, University of Oxford, Parks Road, Oxford
OX1 3PU, UK

###### Abstract

We propose an algorithm inspired by optical coherent Ising machines to solve
the problem of polynomial unconstrained binary optimization (PUBO). We
benchmark the proposed algorithm against existing PUBO algorithms on the
extended Sherrington-Kirkpatrick model and random third-degree polynomial
pseudo-Boolean functions, and observe its superior performance. We also
address instances of practically relevant computational problems such as
protein folding and electronic structure calculations with problem sizes not
accessible to existing quantum annealing devices. The application of our
algorithm to protein folding and quantum chemistry problems sheds light on the
shortcomings of approximating the electronic structure problem by a PUBO
problem, which, in turn, puts into question the applicability of the
unconstrained binary optimization formulation, such as that of quantum
annealers and coherent Ising machines, in this context.

##  I Introduction

A notable example of a combinatorial optimization problem is the problem of
polynomial unconstrained binary/spin optimization (PUBO/PUSO), which is
formulated as follows.

######  Problem 1.

Find the global minimum point of a pseudo-Boolean “energy” function
$H_{k}:\mathbf{s}=\\{-1,+1\\}^{N}\rightarrow\mathbb{R}{}$ such that:

$H_{k}(\mathbf{s})=\sum_{i_{1}}J_{i_{1}}^{(1)}s_{i_{1}}+\ldots+\sum_{i_{1}<\ldots<i_{k}}J_{i_{1}\ldots
i_{k}}^{(k)}s_{i_{1}}\ldots s_{i_{k}}.$ (1)

Here $N$ is the dimension of a problem instance (problem size), $k$ is the
function degree, the summations run from 1 to $N$, and the coupling tensors
$J^{(\cdot)}$ with real coefficients are given.

Historically the dominant subject of the PUBO-related research was its
particular case ($k=2$), also known as the _quadratic unconstrained binary
optimization (QUBO)_ problem. Among naturally arising QUBO problems of
practical value are finding the ground state of a spin-glass Wu  _et al._
(1991) and the problem of constrained via minimisation in very-large-scale
integration design Ben-Ameur  _et al._ (2014). The interest to the QUBO
problem is also motivated by the fact that any PUBO problem can be reduced to
a QUBO problem by increasing the problem dimension Biamonte (2008). In
addition, the NP-complete problem MAX-CUT can be reduced to a QUBO problem
instance Ben-Ameur  _et al._ (2014); Karp (1972). Therefore, if one has an
efficient procedure to solve the QUBO problem, one may efficiently solve any
other problem in the complexity class NP Karp (1972).

Traditionally, QUBO instances were addressed by means of ordinary
combinatorial optimization methods, which include the branch-and-cut method
Barahona  _et al._ (1989), semi-definite programming Poljak and Rendl (1995);
Goemans and Williamson (1994), polynomial time approximation schemes Fernandez
de la Vega (1996), simulated annealing Alkhamis  _et al._ (1998), etc. An
alternative perspective to solve QUBO problems is to build a special-purpose
computing device, and this approach has been widely investigated over the last
decades Boixo  _et al._ (2013); Yamaoka  _et al._ (2015); Puri  _et al._
(2017); Tsukamoto  _et al._ (2017); Pierangeli  _et al._ (2019, 2020); Roques-
Carmes  _et al._ (2020), see Vadlamani et al. Vadlamani  _et al._ (2020) for a
review. The idea is to relate a second degree pseudo-Boolean function to the
energy landscape of some physical system governed by the Ising Hamiltonian,
and to let the system evolve into its ground state.

A special family of special-purpose devices is the class of D-Wave quantum
annealers, which use superconducting qubits coupled by magnetic fields as the
underlying platform Boixo  _et al._ (2013). D-Wave annealers were used for a
number of combinatorial optimization tasks ranging from logistics to quantum
chemistry and material simulation Neukart  _et al._ (2017); Streif  _et al._
(2019a); Li  _et al._ (2018); Perdomo-Ortiz  _et al._ (2012); Boev  _et al._
(2020); Amin  _et al._ (2018); Biamonte  _et al._ (2017); Venturelli  _et al._
(2015); Adachi and Henderson (2015). A downside of these processors is limited
connectivity between qubits: a typical state-of-the-art D-Wave Advantage
processor has 5000 physical qubits, but each is only coupled to 15 others
D-Wave Systems Inc. (2020). The ability of D-Wave processors to demonstrate
quantum computation speed-up is a subject of ongoing research Rønnow  _et al._
(2014).

Another special-purpose devices are coherent Ising machines (CIM), which store
information about optimization variables in optical pulses and use an
optoelectronic feedback loop to implement the couplings between them Wang  _et
al._ (2013); McMahon  _et al._ (2016); Inagaki  _et al._ (2016). CIMs have no
restrictions on the connections between variables and are capable of dealing
with QUBO functions which act on up to 2048 variables and are “dense”, i.e.
have $~{}\Omega(N^{2})$ non-zero terms Inagaki  _et al._ (2016). Despite the
fact that quantum properties of these devices are also subject to debate, it
was recently shown that CIMs significantly outperform D-Wave processors in
dealing with dense QUBO functions Hamerly  _et al._ (2019), although this
claim was disputed by the D-Wave team McGeoch  _et al._ (2018).

An endeavour to overcome disadvantages of existing special-purpose devices
resulted in the development of so-called quantum-inspired optimization
algorithms Arrazola  _et al._ (2019); Goto  _et al._ (2019, 2021); Tatsumura
_et al._ (2021). For example, in Ref. Tiunov  _et al._ (2019), some of us
proposed a “SimCIM” algorithm to solve QUBO problems by classicaly simulating
the optimization procedure of CIMs. SimCIM can be implemented using GPU-
accelerated matrix-vector multiplications and achieves solution speed and
quality that is comparable to that of CIM and higher than many competing
software algorithms.

In Stroev and Berloff (2021) authors make first steps towards quantum-inspired
polynomial optimization — they propose an algorithm adopted from operation of
networks of nonequilibrium polariton condensates. In this paper we advance
further along this direction and introduce _PolySimCIM_ — a generalisation of
SimCIM Tiunov  _et al._ (2019) that enables direct solving of PUBO problems of
degrees $k>2$. This extends the application domain of our algorithm to
problems initially formulated as PUBO, such as finding the tertiary structure
of proteins, calculating the electronic structure of molecules, and solving
box constrained Diophantine equations Perdomo-Ortiz  _et al._ (2019); Babbush
_et al._ (2013); Xia  _et al._ (2017); He  _et al._ (2013). While any PUBO
problem can in principle be reduced to QUBO, this reduction is associated with
an overhead which in some cases can lead to an increase of the problem size
$N$ by several orders of magnitude, as we demonstrate in Sec. III.2.
PolySimCIM thus helps avoiding this overhead. We benchmark our algorithm on
synthetic data of random third-degree pseudo-Boolean functions and reveal its
advantage compared to other state-of-the-art optimization algorithms
(including the one presented in Stroev and Berloff (2021)). Additionally, we
test its performance on the protein folding and electronic structure problems
Babbush  _et al._ (2014); Xia  _et al._ (2017) and demonstrate that the
proposed approach allows dealing with problem sizes previously unreachable for
special-purpose QUBO devices (in particular, the D-Wave machine) or their
simulators.

##  II SimCIM for polynomial binary optimization

We begin by recapping the functionality of SimCIM Tiunov  _et al._ (2019). As
mentioned, this algorithm _sim_ ulates _CIM_ s, which are special-purpose
devices tailored to solve the QUBO problem. Physically each CIM consists of an
optical parametric oscillator formed by degenerate parametric amplifier
(single-mode squeezer) inside a fiber loop, with $N$ optical pulses travelling
in that loop. In each roundtrip, each pulse is subjected to (i) measurement of
the position quadrature and (ii) phase-space displacement along the position
axis. The displacement is computed from the measured position quadratures of
other pulses and the coefficients of the function $H_{2}(\mathbf{s})$ that is
being optimised. After multiple roundtrips, each pulse is amplified to a
coherent state of certain intensity. Thanks to the phase-sensitive nature of
the parametric amplification, the phase $\varphi_{i}$ of each pulse stabilizes
at either to $0$ or to $\pi$. The sequence of these phases is used to obtain
the sequence of $s_{i}=\pm 1$ corresponding to the solution produced by CIM
(1).

To fully simulate a CIM, one should consider evolution of both complex
position and momentum quadratures as well as linear and nonlinear losses
present in the system. However, SimCIM restricts analysis to only the real
part of the position quadrature and simplifies the effect of noise and
nonlinear losses. Each iteration consists of two steps as follows.

Step 1. Calculate the displacement and apply it to the vector of continuous
quadratures (set to zero before the start of simulation):

$\Delta
x_{i_{1}}(t)=\nu(t)x_{i_{1}}+\xi\left(\sum_{i_{2}}J_{{i_{1}}{i_{2}}}^{(2)}x_{i_{2}}+J_{i_{1}}^{(1)}\right)+f_{i_{1}}(t)$
(2)

Here $\nu(t)$ accounts for combination of (time-dependent) pump amplitude and
linear loss, $\xi$ is equivalent to learning rate in the conventional gradient
descent algorithm, and $f_{i}$ is Gaussian noise with variance $\sigma^{2}$.
To accelerate the convergence, we use the momentum method Qian (1999) with the
momentum parameter of $\alpha$ = 0.99.

Step 2. Apply an activation threshold function to each variable to account for
the saturation:

$x_{i_{1}}=\begin{cases}x_{i_{1}}&\mbox{for}\ |x_{i_{1}}|<x_{\rm sat},\\\
x_{\rm sat}&\mbox{otherwise}.\end{cases}$ (3)

In the above equations, $\nu(t)$, $\xi$, $\sigma$, and $x_{\rm sat}$ are
hyperparameters of the algorithm and require fine-tuning for each problem
instance. In the original paper, $\nu(t)$ grows with time according to the
hyperbolic tangent law. After the last iteration, the achieved “solution” is
evaluated as $s_{i}:=\operatorname{sign}{x_{i}}$.

The second term in Eq. (2) is basically a partial derivative with respect to
$x_{i_{1}}$ of a continuous function $H_{k}$($\mathbf{x}$) obtained from
$H_{k}(\mathbf{s})$ by extending the domain of the latter from
$\\{-1,1\\}^{N}$ to $\mathbb{R}^{N}$ (here $k=2$). Hence, Step 1 can be
understood as calculating the displacement vector as

$\Delta\mathbf{x}(t)=\nu(t)\cdot\mathbf{x}+\xi\cdot\nabla
H_{k}(\mathbf{x})+\mathbf{f}(t).$ (4)

with

$\nabla H_{k}(\mathbf{x})=\frac{\partial H_{k}(\mathbf{x})}{\partial
x_{i_{1}}}=\sum_{i_{2}}J_{{i_{1}}{i_{2}}}^{(2)}x_{i_{2}}+J_{i_{1}}^{(1)}$

for the SimCIM case.

Similarly to SimCIM, PolySimCIM follows the gradient of the energy function
(of an arbitrary degree) supported with some noise — so that the algorithm can
get out of local optima. More specifically, Step 1 of SimCIM is replaced in
PolySimCIM by

Step 1′. Calculate the displacement and apply it to the vector of quadratures
according to Eq. (4) with

$\nabla
H_{k}(\mathbf{x})=J_{i_{1}}^{(1)}+\ldots+\sum_{i_{2}<\ldots<i_{k}}J_{i_{1}\ldots
i_{k}}^{(k)}x_{i_{2}}\ldots x_{i_{k}}.$ (5)

The annealing parameter $\nu(t)$ used in calculations is a shifted hyperbolic
tangent which is determined by three hyperparameters $O,D,S$. The extended
form of $\nu(t)$ is as follows:

$\nu(t)=O\cdot\tanh\left(S\left(\frac{t}{N}-0.5\right)\right)-D,$ (6)

where $O$, $D$, $S$ are hyperparameters, $N$ is the number of steps in a
single run.

##  III Benchmarking

###  III.1 Comparison with PUBO algorithms

We begin by benchmarking PolySimCIM against existing optimization algorithms
on pseudo-Boolean functions of third degree:

$H_{3}(\mathbf{s})=\sum_{i_{1}<i_{2}<i_{3}}J_{i_{1}i_{2}i_{3}}s_{i_{1}}s_{i_{2}}s_{i_{3}}.$
(7)

We consider three classes of such functions based on the structure of tensor
$J_{i_{1}i_{2}i_{3}}$. Functions in Class I are those with each element
$J_{i_{1}i_{2}i_{3}}$ randomly set to $1$ or $-1$ with probability $0.5$. In a
sense, they are generalisation of the paradigmatic Sherrington-Kirkpatrick
model Fu and Anderson (1986) to the case of third degree pseudo-Boolean
functions. Functions in Class II also have dense tensors, but each element
$J_{i_{1}i_{2}i_{3}}$ is drawn uniformly from $[-1,1]$ interval. Class III
includes pseudo-Boolean functions with random _sparse_ tensors, which are in
practice obtained by taking a function from Class II and setting each tensor
element to zero (i.e. “dropping” it) with probability $0.9$. We considered
problem sizes $N$ ranging from 10 to 100 in steps of 10. For each problem
dimension, we generated ten pseudo-Boolean functions.

Figure 1: Performance comparison of polynomial optimization algorithms on the
generalised Sherrington-Kirckpatrick, dense random and sparse random tensors.
Top two rows: maximum (a–c) and median (d–f) achieved figure of merit $\eta$
over 2000 runs of each algorithm for 10 problem instances corresponding to a
specific problem class and size. Lines and shaded regions show, respectively,
the medians and the 25-75 percentiles of these maxima and medians over each
set of problem instances. Bottom row (g–i): histograms with kernel density
estimation (KDE)

of $\eta$ for each algorithm aggregated over all runs and problem instances of
size $N=100$.

We compare PolySimCIM with four optimization algorithms, which are also based
on iterative updates of variables (Table 1). The first three of them use
continuous representation of optimization variables and update them based on
the energy function gradient. The first algorithm adapts higher-order
Hopfield-Tank neural networks and is described in Ref. Joya  _et al._ (2002).
The second approach (Leleu) generalises the method of Ref. Leleu  _et al._
(2019) to $k>2$: it uses an additional vector of time-dependent “error
variables” $\mathbf{e}$ that helps escaping local minima. We used the
realisation of the algorithm described in section “Benchmark results on the
G-SET” in Ref. Leleu  _et al._ (2019), straightforwardly adapting the
injection term to the polynomial case. The third algorithm (TGD+CC) Stroev and
Berloff (2021) permits the amplitude vector $\mathbf{x}$ to take on complex
values and makes the gain parameter $\nu(t)$ dependent on these amplitudes.
Additionally, whenever the algorithm appears to reach a steady state, it adds
complex parts to a small number of randomly chosen elements of
$J_{i_{1}i_{2}i_{3}}$, and removes them when the system leaves the local
minimum. The fourth algorithm is discrete greedy search: at each iteration it
moves to the neighbouring vertex of $\\{-1,1\\}^{N}$ hypercube which has the
smallest value of $H_{k}(\mathbf{s})$.

Algorithm |  Variable evolution law
---|---
PolySimCIM |  $\Delta\mathbf{x}(t)=\nu(t)\cdot\mathbf{x}+\xi\cdot\nabla H_{k}(\mathbf{x})+\mathbf{f}(t)$
Hopfield-Tank |  $\Delta\mathbf{x}(t)=-\mathbf{x}+\xi\cdot\nabla H_{k}(\tanh(\mathbf{x}/\beta))$
Leleu |  $\Delta\mathbf{x}(t)=\nu(t)\cdot\mathbf{x}-\mathbf{x}^{3}+\xi\mathbf{e}(t)\nabla H_{k}(\mathbf{x})$
|  $\Delta\mathbf{e}(t)=-\beta(t)(\mathbf{x}^{2}-a(t))\mathbf{e}$
TGD+CC |  $\Delta x_{i_{1}}(t)=\xi\big{(}x_{i_{1}}(\nu_{i_{1}}(t)-|x_{i_{1}}|^{2})+\sum\limits_{i_{2}<\ldots<i_{k}}^{N}J^{(k)}_{i_{1}i_{2}\ldots i_{k}}x_{i_{2}}x_{i_{3}}\ldots x^{*}_{i_{k}}\big{)}$ $\Delta\nu_{i_{1}}(t)=\xi\varepsilon(\rho_{\rm th}-|x_{i_{1}}|^{2})$
Greedy |  ${\mathbf{s}}(t+1)=\arg\displaystyle\min_{\mathbf{s}:|\mathbf{s}-\mathbf{s}(t)|=2}H_{k}(\mathbf{s})$
Table 1: Update rules for the considered polynomial optimization algorithms.
All arithmetic operations (i.e. addition, multiplication, exponentiation etc.)
are element-wise. All variables in equations except for evolving vectors and
pseudo-Boolean functions are hyperparameters of the corresponding algorithms.
After each iteration step of the PolySimCIM the activation threshold function
determined by the Eq. (3) is applied to the amplitudes.

For each problem instance, we performed 2000 simultaneous runs of each
algorithm. The hyperparameters of PolySimCIM, Leleu, Hopfield-Tank, and TGD+CC
algorithms were optimized with a machine-learning online optimization package
M-LOOP Wigley  _et al._ (2016) (see Appendix B for details). The realisation
of TGD+CC was provided by the authors of Ref. Stroev and Berloff (2021). The
number of steps was equal to 1000 for each method.

A set of hyperparameters was found for each class of energy functions on a
single problem instance of each dimension and then used for all other problems
of the same class and dimension. The performance of each algorithm could be
further optimized by picking a set of hyperparameters for each problem
instance individually. Even deeper optimization could be reached by choosing
the type of nonlinearity for a given algorithm and problem instance Böhm  _et
al._ (2021). However, such in-depth individual optimization would slow down
the process if multiple problem instances of similar types are to be solved,
so we choose not to use this strategy.

The starting configurations of the Hopfield-Tank and Leleu algorithms were
randomly initialised from a normal distribution $\mathcal{N}(0,10^{-4})$, and
for the Greedy algorithm the starting configurations were chosen randomly and
uniformly from $\\{-1,1\\}^{N}$. Different runs of the Hopfield-Tank, Leleu
and Greedy algorithms produced different results due to these random initial
conditions. For the PolySimCIM and TGD+CC algorithms, different results across
runs were obtained due to the random noise present in their update routines.
Each run of each algorithm produced an estimate $\hat{\mathbf{s}}_{\rm min}$
for the optimal spin configuration. Next, we picked the minimum among all
$5\cdot 2000=10000$ obtained values of $\hat{\mathbf{s}}_{\rm min}$ (i.e.
among the results produced by all algorithms), which we denote as
$\mathbf{s}_{\rm min}$. As a figure of merit for each run we use the following
ratio:

$\eta=\frac{H_{3}(\hat{\mathbf{s}}_{\rm min})}{H_{3}(\mathbf{s}_{\rm min})}.$

The energy functions can take both positive and negative values and so does
$\eta$. However, in all our experiments $H_{3}(\mathbf{s}_{\rm min})<0$, which
implies that $\eta\in(-\infty,1]$, and therefore higher solution quality
corresponds to $\eta$ being closer to unity.

To compare the performance, we aggregate the information about values of
$\eta$ achieved in different runs in three possible ways, as depicted in Fig.
1. First, for each problem instance, we pick the best value of $\eta$ achieved
by the given algorithm among all runs on a given problem instance and we take
the median of these values among all problem instances of a given problem size
[Fig. 1(a–c)]. Second, we find the median of $\eta$’s achieved by the given
algorithm among all runs on a given problem instance, and then plot the
medians of these medians over the 10 problem instances for each problem size
[Fig. 1(d–f)]. Finally, we plot the histograms of the $\eta$ values obtained
for all problems of size $N=100$ with each algorithm [Fig. 1(g–i)].

Figures 1(a–c) show that the Hopfield-Tank and Greedy algorithms perform worse
than others: they systematically fail to achieve the best known solutions for
large problem sizes, except for the Class I pseudo-Booleans, for which the
Greedy algorithm successfully reaches $\eta=1$ [Fig. 1(a–c)]. In contrast,
PolySimCIM manages to achieve $\eta=1$ for each problem instance of each
problem size and for all classes of pseudo-Boolean functions. Furthermore, as
evident from the histograms in Fig. 1(g–i), the quality of PolySimCIM
solutions is consistently high for different runs of the algorithm and
different problem instances. In this latter respect, PolySimCIM surpasses
Leleu, which is also able to find the optimal solution in most cases, but
shows higher variances among the runs and problem instances.

Algorithm |  Execution time
---|---
PolySimCIM |  2.7 s
Hopfield-Tank |  1.2 s
Leleu |  4.4 s
TGD+CC |  24.3 s
Greedy |  8522.8 s
Table 2: Average execution time of the algorithms for tensors with dimension
$N=100$ on the NVIDIA 2080Ti GPU. PolySimCIM, Hopfield-Tank, Leleu make 1000
iterations. TGD+CC makes 1000 iterations before complex coupling switching and
1000 iterations after it. Each algorithm made 2000 simultaneous runs.

###  III.2 Comparison with quadratic solvers

Since all the algorithms benchmarked above are still an active area of
research, it is of interest to compare results with commonly used blackbox
solvers. However, to our knowledge, there exist no such solvers specifically
optimized for PUBO with $k>2$. Existing solvers are limited to the quadratic
domain, so, in order to use them, we had to convert our problems to the QUBO
form. To this end, we used two methods: one based on the paper by Xia, Bian,
Kais Xia  _et al._ (2018) and the other from the open-source Python library
Qubovert Tiosue (2020). While the Xia-Bian-Kais method was able to produce
QUBOs with fewer variables, the transformations took considerably longer to
execute. However, for QUBOs with fewer qubits, quadratic solvers execute more
rapidly, compensating for the pre-processing time. Both methods showed similar
performance in terms of the minimum energies found, so here we report the
results obtained with Qubovert pre-processing (Fig. 2).

The solvers used for benchmarking were as follows.

  1. i

Simulated annealing, implemented with D-Wave’s open source code Inc. (2021);

  2. ii

Gurobi’s mixed integer programming solver Gurobi optimisation, LLC (2022).
Gurobi utilizes a plethora of methods such as branch-and-cut, deterministic
parallel, non-traditional search, heuristics, solution improvement, cutting
planes, and symmetry breaking algorithms.

  3. iii

IBM CPLEX IBM ILOG CPLEX optimisation Studio (2022) which utilizes methods
such as simplex, barrier interior point method, and second-order cone
programming.

To match execution times with the PUBO solvers, the time limit parameter for
(ii) and (iii) was set to two minutes. However the solvers did not conform to
this time limit in all cases: some executions took as long as ten hours.
Because of this slow performance, we collected only a single sample for each
problem instance with these algorithms. For (i), since the time limit
parameter setting was not available, other parameters were tuned to have the
execution time close to two minutes. All QUBO solvers were executed on a 2.3
GHz Quad-Core Intel Core i7 processor.

Figure 2: Comparison of the overhead of transformations with Qubovert and Xia-
Bian-Kais method in locality reduction of original PUBO problems: (a) problem
size (number of bits); (b) total number of terms in the energy function.

As seen in Fig. 3, the blackbox QUBO solvers performed significantly poorer
than most of the PUBO solvers discussed in Sec. III. This is likely due to the
significant overhead introduced by the locality reduction scheme [Fig. 2(a)].
A plausible inference from these results is that solving optimization problems
in their native form generally leads to higher quality solutions. Even though
the QUBO solvers we used are state-of-the-art in their own domain, they were
unable to perform competitively outside it. A perhaps more speculative
conclusion is that solving problems of higher than the second degree with
current quantum annealers is not a advisable approach, as the latter are
limited to the quadratic domain.

Figure 3: Performance comparison of PolySimCIM against QUBO solvers.
PolySimCIM and simulated annealing produced 2000 samples per problem instance,
of which the best $\eta$ value is chosen. Gurobi and CPLEX produced only one
sample and were run only once per problem instance due to their slow
convergence. Similar to Fig. 1, the median of these values over 10 problem
instances corresponding to a specific problem class and size is shown. Shaded
regions show the 25-75 percentiles of these values over each set of problem
instances.

##  IV Application to real-world problems

We now apply PolySimCIM to solve protein folding and electronic structure
problems. Both problems can be formulated in the PUBO framework. In previous
research, they were subsequently converted to QUBO and then solved using
quantum annealers, such as D-Wave 2000Q Perdomo-Ortiz  _et al._ (2012); Xia
_et al._ (2017). Applying PolySimCIM allows us to eliminate this second
conversion, and solve the problem directly as a PUBO, thereby achieving
superior performance. This being said, the “PUBO route” is not the only, and
oftentimes not the optimal, method to tackle them McArdle  _et al._ (2020a);
David Sherrill and Schaefer (1999); Thachuk  _et al._ (2007); Guo  _et al._
(2017), as we see below.

###  IV.1 Lattice protein folding

Being initially a mere chain of amino-acid residues, soon after the synthesis
each protein molecule folds into a unique spatial structure, corresponding to
the minimum of free energy. This structure determines most of the protein
behaviour, while proteins with damaged or incorrect folds are unable to
function properly. As proteins form arguably the most versatile and ubiquitous
class of biomolecules, the problem of finding the spatial configuration of a
protein from its primary amino-acid sequence — also known as _the protein
folding problem_ — is of paramount importance in computational biology Dill
_et al._ (2008); Šali  _et al._ (1994); Dill and MacCallum (2012); Englander
and Mayne (2014).

Over the last half-century, this problem has been approached with a variety of
methods. One family of approaches considers interaction field forces between
different amino-acid residues and solves equations of motion for the latter,
thus recovering both the folding process and its outcome Hansmann and Okamoto
(1999). The second large family are Monte Carlo methods, where random trial
moves of a protein chain are iteratively tested. These moves are then either
accepted or rejected based on the Boltzmann weight of the newly obtained
conformation Zhang and Chou (1992); Kolinski and Skolnick (1994). Finally,
techniques based on machine learning have recently come to the fore, with a
deep neural network predicting the distances between amino-acid pairs and the
angles between chemical bonds connecting those pairs Noé  _et al._ (2020);
Senior  _et al._ (2020). Note that this approach relies heavily on the
extensive database of known protein structures acquired over the last decades.

Figure 4: Lattice model of protein folding. Interactions between neighbouring
amino-acids are shown as dark cyan dotted lines. Encoding of the first four
turns is shown with last two red digits indicating the turn direction. The
figure is specific to the PHPPHPPHPPH sequence in the PH model.

However, despite the remarkable progress achieved so far, no “silver bullet”
algorithm to fully predict the three-dimensional structure of a protein is yet
developed. To gain insight into the mechanism of protein folding, simplified
models are sometimes considered, such as the 2D square lattice model (Fig. 4).
Babbush et al. showed that such a model can be reduced to PUBO, which paves
the way to applying quantum annealers and their simulators in the protein
folding problem Babbush  _et al._ (2014). Although lattice protein folding is
a toy-model, it has been shown to predict model protein tertiary structures to
remarkable accuracy Babej  _et al._ (2018). In what follows, we briefly
discuss the idea behind this reduction.

Suppose a protein consists of $M$ amino-acid residues and is embedded in a 2D
square lattice. To fully define a planar conformation of the protein, one
should specify directions of $M-1$ turns of the protein chain. There are only
four possible turn directions. Under the conventions of the “Turn Ancilla”
model Babbush  _et al._ (2014), each turn is encoded in a pair of binary
variables as follows: 01 — “right”, 11 — “up”, 10 — “left”, 00 — “down” (Fig.
4). Hence the fold of a protein can be described with $2(M-1)$ binary
variables. In fact, the first three bits can be set to arbitrary values
without loss of generality, and so the fold configuration in this model is
described by $2M-5$ bits. The configuration energy $E_{\rm pair}(\mathbf{s})$
is then made up of interaction energies between pairs of neighbouring amino-
acids. There exist a variety of models for evaluating these energies. In a
simple hydrophobic-polar (HP) model all amino-acids are classified as
hydrophobic or polar, and the interaction strength depends on which of these
two classes the interacting molecules belong to. In a more intricate Miyazawa-
Jernigan (MJ) model, the energy depends on the specific types of the
interacting amino-acids. Both models are detailed in Table 3.

Two further modifications are needed to express the interaction energy in the
PUBO form. First, unphysical fold configurations (i.e. those in which chain
fragments overlap or intersect) are penalised by an additional energy term
$E_{\rm penalty}(\mathbf{s})$. Second, the initial vector of $2M-5$ variables
is supplemented with ancillary variables, which are necessary to determine
which amino-acids are neighbours. The details of these transformations are
beyond the scope of our paper, but detailed in Ref. Babbush  _et al._ (2014).
The resulting energy function $E(\mathbf{s})=E_{\rm pair}(\mathbf{s})+E_{\rm
penalty}(\mathbf{s})$ is then of fourth degree. It can be either reduced to
the second-degree (QUBO) function at the cost of adding more binary variables
(Fig. 5) or treated directly using PolySimCIM.

Figure 5: QUBO and PUBO problem sizes in the protein folding problem under the
lattice model as a function of the amino-acid sequence length. The calculation
is outlined in Appendix A.

We applied PolySimCIM to four proteins: three of lengths 6, 7 and 8 in the
framework of the MJ model, and one of length 11 under the HP model. The choice
of the first chain (PSVKMA) has been motivated by an existing result obtained
with a D-Wave quantum annealer Perdomo-Ortiz  _et al._ (2012). In this work,
the ground state of the PSVKMA sequence was found with a success probability
of $0.13\%$ on the D-Wave quantum annealer using 81 qubits. The next two
chains (PSVKMAP and PSVKMAPS) are obtained from the first one by sequentially
appending arbitrary amino-acids. The largest chain (PHPPHPPHPPH) was chosen
from Ref. Rostoker and Lotun (2005).

The results are summarized in Table 4. In all four cases the fold
configuration with the lowest known energy has been found successfully, which
was confirmed by a brute force search. The evolution of variables in a
successful run of PolySimCIM for the PSVKMAP sequence is shown in Fig. 6(a).
One can see that, while most variables saturate relatively quickly, a few
others need a considerable number of iterations to “commit” to a binary value.
As a rule, the variables that take longer to stabilise are the ancillary ones,
rather than the $2M-5$ variables encoding the spatial structure of the
protein. This is evident in Fig. 6(b,c): at $t=750$ the protein molecule has
already folded correctly, but $E(\mathbf{s})=9$, which differs from the final
optimal value $E_{\min}=-6$. In other words, the algorithm may be stopped
before all variables are stabilised and still yield the optimal conformation.

Hydrophobic-polar | Miyazawa-Jernigan
---|---
HH |  HP |  PP | PS | PK | PA | SM | VA | VS
$-1$ |  0 |  0 | $-0.5$ | $-1$ | $-2$ | $-3$ | $-4$ | $-5$
Table 3: Values of interaction strengths between amino-acids for the
Hydrophobic-polar and Miyazawa-Jernigan models. In the framework of the latter
P, S, V, K, M, A stand for Proline, Serine, Valine, Lysine, Methionine and
Alanine, respectively. Miyazawa-Jernigan model
---
Primary sequence |  Energy |  $N_{\rm PUBO}$ |  $N_{\rm QUBO}$ |  $p_{\rm succ}$
PSVKMA |  $-6$ |  19 |  28 |  1.0
PSVKMAP |  $-6$ |  33 |  49 |  1.0
PSVKMAPS |  $-9$ |  48 |  73 |  $10^{-4}$
Hydrophobic-polar model
PHPPHPPHPPH |  $-4$ |  104 |  168 |  $2.5\cdot 10^{-5}$
Table 4: Results of PolySimCIM applied to four amino-acid sequences. The
columns are the energy of the optimal fold, number of bits in the QUBO and
PUBO settings and the success probability $p_{\rm succ}$, defined as the
number of runs where the algorithm achieved the ground state energy divided by
the total number of runs.  Figure 6: (a) Evolution of amplitudes in a
successful run of PolySimCIM for PSVKMAP amino-acid sequence; (b) Evolution of
$E(\mathbf{x})$ in the same run ($E_{\min}=-6$); (c) Protein conformations at
five moments in algorithm runtime [red circles in Fig. 6(b)].

###  IV.2 The electronic structure problem

Most of the molecule properties can be derived from its ground state wave
function (also known as _the electronic structure_), and thus an ability to
calculate the latter — or at least the ground state energy — is crucial for
quantum chemists. To find it, one usually writes down the Schrödinger equation
for a given molecule, transforms it to a second quantised form, and then
applies one of the quantum chemistry methods developed over the last half-
century — we refer the reader to a classic textbook on quantum chemistry Szabo
and Ostlund (2012) for more details. Alas, the Schrödinger equation in the
second quantised form is computationally tractable only for small molecules as
the size of system Hilbert space grows exponentially with the number of
electrons.

An approach to facilitate electronic structure calculations with quantum
annealers was proposed in Ref. Xia  _et al._ (2017). Its first step is to
rewrite the second quantized Hamiltonian in the qubit form using any of
relevant transformations developed so far (e.g. Jordan-Wigner, Bravyi-Kitaev
etc.) McArdle  _et al._ (2020b). The Hamiltonian then takes the form

$\hat{H}=\sum_{i,\alpha}h_{\alpha}^{i}\hat{\sigma}_{\alpha}^{i}+\sum_{\begin{subarray}{c}i,j\\\
\alpha,\beta\end{subarray}}h_{\alpha\beta}^{ij}\hat{\sigma}_{\alpha}^{i}\hat{\sigma}_{\beta}^{j}+\sum_{\begin{subarray}{c}i,j,k\\\
\alpha,\beta,\gamma\end{subarray}}h_{\alpha\beta\gamma}^{ijk}\hat{\sigma}_{\alpha}^{i}\hat{\sigma}_{\beta}^{j}\hat{\sigma}_{\gamma}^{k}+\ldots.$
(8)

Here Latin alphabet indices enumerate the qubits and range from $1$ to $M$,
while Greek letters point which of Pauli operators $\left\\{I,X,Y,Z\right\\}$
constitute the term.

The Hamiltonian (8) includes all three Pauli operators and is hence inherently
quantum. Ref. Xia  _et al._ (2017) makes it compatible with the classical PUBO
problem by mapping it onto a Hamiltonian

$H=\sum_{\mathcal{S}\subset\\{1,\ldots,rM\\}}h_{\mathcal{S}}\prod_{i\in\mathcal{S}}s^{i},$
(9)

which is a function of $rM$ bits that can take on values $s\in\\{-1,1\\}$.
Here $r$ is some integer number which allows fine-tuning the method accuracy.
The details of this mapping are outside the scope of this paper, but can be
found in Ref. Xia  _et al._ (2017). Note that the case of $r=1$ corresponds to
the Hartree-Fock solution.

To employ quantum annealers, Xia  _et al._ Xia  _et al._ (2017) propose to
subsequently reduce this Hamiltonian to a 2-local (QUBO) form. This reduction
however introduces a dramatic scaling overhead. The authors estimate the
number of terms in Hamiltonian (9) to scale as
$\mathcal{O}\left(2^{M}r^{2}M^{4}\right)$ and the number of bits in the final
QUBO Hamiltonian as $\mathcal{O}\left(2^{M}r^{2}M^{7}\right)$. Streif et al.
report that due to such prohibitive scaling only very small molecules can
actually be addressed with quantum annealers Streif  _et al._ (2019b). In
particular, the authors found the ground state of H2 ($M=2$) with the D-Wave
2000Q quantum annealer, requiring the scaling factor value as high as $r=16$
to obtain a good agreement with exact quantum calculation. However, the full
treatment of a somewhat more complex molecule LiH ($M=10$) proved to be
impossible on this machine as the final QUBO Hamiltonian required more qubits
than was available in D-Wave 2000Q. To embed this problem onto this machine,
the authors resorted to restricted active space methods and limited the number
of orbitals which can be occupied by electrons, but then they were not able to
obtain states with energies lower than those given by the Hartree-Fock method
Streif  _et al._ (2019b).

Figure 7: Comparison of QUBO and PUBO approaches to electronic structure
calculations: (a) problem size (number of bits); (b) number of terms in the
Hamiltonians obtained according to Xia  _et al._ (2017). The calculation for
LiH is limited to $r\leq 5$ because higher values of $r$ required
prohibitively long compute time. Figure 8: a) Potential energy surfaces
restored by PolySimCIM for (a) H2, $r\in\\{2,3,16\\}$; (b) LiH, $r=2$. Solid
black curves represent exact potential energy surfaces, solid red curves
correspond to the Hartree-Fock solution. FCI stands for Full Configuration
Interaction Szabo and Ostlund (2012), a quantum chemistry method, which
enables exact determination of the ground state energy.

These circumstances motivated us to apply PolySimCIM directly to the
Hamiltonian (9) in order to avoid exponential scaling of the number of bits.
We attempted PUBO optimization for H2 and LiH molecules. These molecules are
linear, and we optimised the ground state energy as a function of the
internuclear distance. We used the quantum chemistry package _Psi4_ Parrish
_et al._ (2017) to calculate two- and four-body integrals in the fermionic
molecular Hamiltonians. We also used a software package _OpenFermion_ McClean
_et al._ (2020) and applied symmetry-conserving Bravyi-Kitaev transformation
Bravyi  _et al._ (2017) to both molecules so that the resulting Hamiltonians
(8) had 2 and 10 qubits correspondingly.

Figure 7(a) shows the characteristic numbers of bits in the PUBO and QUBO
Hamiltonians for the two molecules. For H2, these numbers are different by
only $1.5$, but for LiH this gap reaches several orders of magnitude. While
this comparison appears favourable for PUBO, the problem complexity depends
not only on the number of bits, but also on the number of terms in the
Hamiltonian. These numbers are comparable for the two cases and scale
exponentially with $r$ [Fig. 7(b)]. Even for $r=2$, the PUBO and QUBO
Hamiltonians have, respectively, $\sim 2\cdot 10^{5}$ and $\sim 10^{6}$ terms,
compared to 631 in the initial molecular Hamiltonian (8). This scaling arises
because, unlike the problems studied in previous sections, the degree of terms
in the PUBO Hamiltonian (9) is limited only by the number of bits. As a
result, while the problem still appears simpler in the PUBO formulation, it
remains computationally expensive to solve.

The results are presented in Fig. 8. The calculations for H2 have been
performed for $r\in\\{2,3,16\\}$, while for LiH we considered only $r=2$ as
higher values of $r$ would result in prohibitive computation overheads. For
each setting, PolySimCIM made 500 independent runs of 2000 update steps. The
same hyperparameters were used for all internuclear distances. In all cases
studied, PolySimCIM produced solutions equal to the “brute force” ones
obtained by exact optimization of the corresponding PUBO Hamiltonians (9); it
was possible to calculate the latter as the Hamiltonians had at most 32 bits.

However, these solutions were not always consistent with the ground states of
the quantum Hamiltonian (8). For the H2 molecule, increasing the value of $r$
leads to consistent improvement of the approximation. For $r=16$, the computed
energy lies within the chemical accuracy from the surface obtained by direct
diagonalization of the Hamiltonian (8). For LiH, the calculation with $r=2$
did not provide better solution than Hartree-Fock. Thus, even though
PolySimCIM allowed us to perform calculations for system sizes that were
previously unachievable using quantum annealers, the approximating the quantum
molecular Hamiltonian by a PUBO Hamiltonian does not appear to be a viable
approach for electronic structure calculations. The reasons are (i)
exponential scaling of the number of terms in the PUBO Hamiltonian and (ii)
high values of $r$ required to obtain the desired accuracy.

##  V Discussion

We have extended the quantum-inspired SimCIM algorithm to cover polynomial
unconstrained binary optimization with degrees $k>2$. We benchmarked
PolySimCIM on a variety of PUBO graphs in comparison with (1) other QUBO
algorithms modified to solve PUBO and (2) blackbox QUBO algorithms, for which
the PUBO problems were converted to the QUBO format. In both cases, PolySimCIM
proved superior.

We have also applied PolySimCIM to the protein folding and electronic
structure problems, surpassing the performance of quantum annealers and their
simulators. However, despite this success, large-scale electronic structure
calculations are unlikely to succeed with this approach. The fault appears to
be not in the algorithm itself, but rather in mapping this problem onto PUBO:
the resulting energy function scales exponentially with the molecule size and
quickly becomes intractable. This reinforces our hunch from Sec. III.2 that
optimization problems are best solved in their original format.

The main conclusions of this work are that (1) PolySimCIM is a state-of-the-
art algorithm for PUBO problems, and (2) for not “native PUBO” combinatorial
or quantum optimization problems, PolySimCIM provides a good testbed to find
out whether conversion to PUBO is the optimal route towards solution.

We suppose that the performance of PolySimCIM can be improved by adding time-
dependent correction to amplitude inhomogeneity akin to Ref. Goto  _et al._
(2021). This modification can be subject of future investigations.

The extension from QUBO to PUBO is also possible in hardware. For example, in
the optoelectronic implementation of a fiber optic coherent Ising machine
Honjo  _et al._ (2021); McMahon  _et al._ (2016), the gradient of pseudo-
Boolean energy function is computed on a FPGA by multiplying the matrix
$J_{i_{1}}i_{2}^{(2)}$ by a vector of measured quadratures. This operation can
be extended to PUBO by multiplying the PUBO tensor $J_{i_{1}\ldots
i_{k}}^{(k)}x_{i_{2}}$ by multiple copies of the measured quadrature vector.
Moreover, PUBO can be implemented by pure analogue means, without any digital
computing. Aside from the aforementioned Ref. Stroev and Berloff (2021), in
which an implementation via polariton condensates is proposed, a high-order
coherent Ising machine can be based on the frequency conversion in a nonlinear
crystal Kumar  _et al._ (2020). Such an experimental setup can solve
polynomial pseudo-Boolean functions with all-to-all connections.

###### Acknowledgements.

This work was supported by Russian Science Foundation (19-71-10092). The part
on quantum chemistry is also partially supported by Nissan Research (QUBO
analysis). D.A. acknowledges Nikita Stroev for fruitful discussions and help
with TGD+CC algorithm realisation. A.F. is supported by the Priority 2030
program at the National University of Science and Technology “MISIS” under the
project K1- 2022-027 (analysis of the method).

##  Appendix A Number of bits in the protein folding problem.

We briefly recap the method for calculating the number of bits in a protein
molecule consisting of $M$ interacting amino acids Babbush  _et al._ (2014).
Three terms $n_{\rm phys}$, $n_{\rm penalty}$, $n_{\rm pair}$ and $n_{\rm
reduction}$ contribute to the total number of bits in the PUBO model. The
first term corresponds to the number of bits $n_{\rm phys}$ encoding the
spatial configuration of the protein. The second term $n_{\rm penalty}$
emerges from the penalty part of the pseudo-Boolean function designed to avoid
physically impossible conformations. The third term $n_{\rm pair}$ is the
number of ancilla variables, which is equal to a number of potential
interactions between amino acids in the sequence. Let us note that the
interactions between amino acids are possible only if the difference between
their positions in the primary sequence is greater than 3. It means that first
and third aminoacids cannot interact, but the first and fourth can. When the
PUBO pseudo-boolean function is reduced to the QUBO one, an additional term
$n_{\rm reduction}$ is added to the number of bits:

$\begin{gathered}N_{\rm PUBO}=n_{\rm phys}+n_{\rm penalty}+n_{\rm pair};\\\
N_{\rm QUBO}=n_{\rm PUBO}+n_{\rm reduction}.\end{gathered}$ (10)

This reduction can be realised in a variety of ways, and choosing the most
efficient procedure is an NP-hard problem Boros and Hammer (2002). In this
manuscript, we use the method of Ref. Babbush  _et al._ (2014) due to its high
efficiency in the case of the “Turn Ancilla” model. The explicit formulae for
calculating the number of bits are the following:

$\begin{gathered}n_{\rm phys}=2M-5;\\\ n_{\rm
penalty}=\sum_{i=4}^{M-4}\sum_{j=i+4}^{M}\lceil\log_{2}(i-j)^{2}\rceil[(1+i-j)\mod
2];\\\ n_{\rm pair}=\sum_{i=1}^{M-3}\sum_{j=i+3}^{M}[(i-j)\mod 2];\\\ n_{\rm
reduction}=\sum_{i=1}^{2M-7}\sum_{j=i+2}^{2M-5}[(i-j+1)\mod 2].\end{gathered}$
(11)

##  Appendix B Hyperparameter Search

To find optimal hyperparameters for the algorithms, we use the M-LOOP machine
learning package Wigley  _et al._ (2016). Within M-LOOP, we use the
“differential evolution” method for sampling the data set to train the neural
network and the “neural net” controller for the optimization itself. The cost
function has the following form of weighted sum:

$\mathrm{cost}=0.25\min{\left[\frac{\langle
H(\mathbf{s})\rangle}{\left|\langle
H_{0}(\mathbf{s})\rangle\right|},1\right]}+\min\left[\frac{\min_{N}H(\mathbf{s})}{\left|\min_{N}H_{0}(\mathbf{s})\right|},1\right]$

Here $\\{H(\mathbf{s})\\}$ is the set of values of the energy function
obtained from all runs of the algorithm. The first term corresponds to the
ratio between the mean value of the energy function over all runs of the
algorithm in each M-LOOP step and the modulus of the mean value in the first
step. The second term corresponds to the ratio between the minimum (best)
obtained value of the energy function over all runs and the modulus of the
minimum value in the first step of M-LOOP. The cost function is bounded by the
value 1.25 from above. This normalization allows us to improve the training of
the M-LOOP neural networks and to speed up the process of finding the optimal
hyperparameters. 1500 steps of M-LOOP were realized for each algorithm. The
parameters are searched on the logarithmic grid in the interval
$[10^{-8},10^{3}]$ with the exception for $\sigma$, which was searched in the
interval $[10^{-8},10^{-0.1}]$.

## References

  * Wu  _et al._ (1991) W. Wu, B. Ellman, T. Rosenbaum, G. Aeppli, and D. Reich, Physical Review Letters 67, 2076 (1991).
  * Ben-Ameur  _et al._ (2014) W. Ben-Ameur, A. R. Mahjoub, and J. Neto, in  _Paradigms of Combinatorial Optimization: Problems and New Approaches_ , edited by V. T. Paschos (John Wiley & Sons, Incorporated, 2014) Chap. 6, pp. 131–172, 2nd ed.
  * Biamonte (2008) J. D. Biamonte, Physical Review A 77, 052331 (2008).
  * Karp (1972) R. M. Karp, in  _Complexity of computer computations_ (Springer, 1972) pp. 85–103.
  * Barahona  _et al._ (1989) F. Barahona, M. Jünger, and G. Reinelt, Mathematical Programming 44, 127 (1989).
  * Poljak and Rendl (1995) S. Poljak and F. Rendl, SIAM Journal on Optimization 5, 467 (1995).
  * Goemans and Williamson (1994) M. X. Goemans and D. P. Williamson, in  _Proceedings of the twenty-sixth annual ACM symposium on Theory of computing_ (1994) pp. 422–431.
  * Fernandez de la Vega (1996) W. Fernandez de la Vega, Random Structures & Algorithms 8, 187 (1996).
  * Alkhamis  _et al._ (1998) T. M. Alkhamis, M. Hasan, and M. A. Ahmed, European Journal of Operational Research 108, 641 (1998).
  * Boixo  _et al._ (2013) S. Boixo, T. Albash, F. M. Spedalieri, N. Chancellor, and D. A. Lidar, Nature Communications 4, 10.1038/ncomms3067 (2013), arXiv:1212.1739 .
  * Yamaoka  _et al._ (2015) M. Yamaoka, C. Yoshimura, M. Hayashi, T. Okuyama, H. Aoki, and H. Mizuno, IEEE Journal of Solid-State Circuits 51, 303 (2015).
  * Puri  _et al._ (2017) S. Puri, C. K. Andersen, A. L. Grimsmo, and A. Blais, Nature communications 8, 1 (2017).
  * Tsukamoto  _et al._ (2017) S. Tsukamoto, M. Takatsu, S. Matsubara, and H. Tamura, Fujitsu Sci. Tech. J 53, 8 (2017).
  * Pierangeli  _et al._ (2019) D. Pierangeli, G. Marcucci, and C. Conti, Phys. Rev. Lett. 122, 213902 (2019).
  * Pierangeli  _et al._ (2020) D. Pierangeli, G. Marcucci, and C. Conti, Optica 7, 1535 (2020).
  * Roques-Carmes  _et al._ (2020) C. Roques-Carmes, Y. Shen, C. Zanoci, M. Prabhu, F. Atieh, L. Jing, T. Dubček, C. Mao, M. R. Johnson, V. Čeperić, J. D. Joannopoulos, D. Englund, and M. Soljačić, Nat. Commun. 11, 1 (2020).
  * Vadlamani  _et al._ (2020) S. K. Vadlamani, T. P. Xiao, and E. Yablonovitch, Proceedings of the National Academy of Sciences 117, 26639 (2020), https://www.pnas.org/doi/pdf/10.1073/pnas.2015192117 .
  * Neukart  _et al._ (2017) F. Neukart, G. Compostella, C. Seidel, D. von Dollen, S. Yarkoni, and B. Parney, Frontiers in ICT 4, 1 (2017), arXiv:1708.01625 .
  * Streif  _et al._ (2019a) M. Streif, F. Neukart, and M. Leib, Lecture Notes in Computer Science (including subseries Lecture Notes in Artificial Intelligence and Lecture Notes in Bioinformatics) 11413 LNCS, 111 (2019a), arXiv:1811.05256 .
  * Li  _et al._ (2018) R. Y. Li, R. Di Felice, R. Rohs, and D. A. Lidar, npj Quantum Information 4, 14 (2018).
  * Perdomo-Ortiz  _et al._ (2012) A. Perdomo-Ortiz, N. Dickson, M. Drew-Brook, G. Rose, and A. Aspuru-Guzik, Scientific Reports 2, 1 (2012), arXiv:1204.5485 .
  * Boev  _et al._ (2020) A. S. Boev, A. S. Rakitko, S. R. Usmanov, A. N. Kobzeva, I. V. Popov, V. V. Ilinsky, E. O. Kiktenko, and A. K. Fedorov, Genome assembly using quantum and quantum-inspired annealing (2020), arXiv:2004.06719 .
  * Amin  _et al._ (2018) M. H. Amin, E. Andriyash, J. Rolfe, B. Kulchytskyy, and R. Melko, Physical Review X 8, 021050 (2018).
  * Biamonte  _et al._ (2017) J. Biamonte, P. Wittek, N. Pancotti, P. Rebentrost, N. Wiebe, and S. Lloyd, Nature 549, 195 (2017).
  * Venturelli  _et al._ (2015) D. Venturelli, D. J. J. Marchand, and G. Rojo, Quantum annealing implementation of job-shop scheduling (2015), arXiv:1506.08479 .
  * Adachi and Henderson (2015) S. H. Adachi and M. P. Henderson, Application of quantum annealing to training of deep neural networks (2015), arXiv:1510.06356 .
  * D-Wave Systems Inc. (2020) D-Wave Systems Inc., Advantage datasheet (2020).
  * Rønnow  _et al._ (2014) T. F. Rønnow, Z. Wang, J. Job, S. Boixo, S. V. Isakov, D. Wecker, J. M. Martinis, D. A. Lidar, and M. Troyer, Science 345, 420 (2014).
  * Wang  _et al._ (2013) Z. Wang, A. Marandi, K. Wen, R. L. Byer, and Y. Yamamoto, Physical Review A - Atomic, Molecular, and Optical Physics 88, 10.1103/PhysRevA.88.063853 (2013), arXiv:arXiv:1311.2696v1 .
  * McMahon  _et al._ (2016) P. L. McMahon, P. L. McMahon, A. Marandi, Y. Haribara, R. Hamerly, C. Langrock, S. Tamate, T. Inagaki, H. Takesue, S. Utsunomiya, K. Aihara, R. L. Byer, M. M. Fejer, H. Mabuchi, and Y. Yamamoto, Science 5178 (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, O. Tadanaga, H. Takenouchi, K. Aihara, K. I. Kawarabayashi, K. Inoue, S. Utsunomiya, and H. Takesue, Science 354, 603 (2016).
  * 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, K. Enbutsu, T. Umeki, R. Kasahara, S. Utsunomiya, S. Kako, K. I. Kawarabayashi, R. L. Byer, M. M. Fejer, H. Mabuchi, D. Englund, E. Rieffel, H. Takesue, and Y. Yamamoto, Science Advances 5, 1 (2019), arXiv:1805.05217 .
  * McGeoch  _et al._ (2018) C. C. McGeoch, W. Bernoudy, and J. King, arXiv preprint arXiv:1807.00089 (2018).
  * Arrazola  _et al._ (2019) J. M. Arrazola, A. Delgado, B. R. Bardhan, and S. Lloyd, arXiv preprint arXiv:1905.10415 , 1 (2019).
  * Goto  _et al._ (2019) H. Goto, K. Tatsumura, and A. R. Dixon, Science Advances 5, 1 (2019).
  * 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, 1 (2021).
  * Tatsumura  _et al._ (2021) K. Tatsumura, M. Yamasaki, and H. Goto, Nature Electronics 4, 208 (2021).
  * Tiunov  _et al._ (2019) E. S. Tiunov, A. E. Ulanov, and A. I. Lvovsky, Optics Express 27, 10288 (2019), arXiv:1901.08927 .
  * Stroev and Berloff (2021) N. Stroev and N. G. Berloff, Physical Review Letters 126, 050504 (2021).
  * Perdomo-Ortiz  _et al._ (2019) A. Perdomo-Ortiz, A. Feldman, A. Ozaeta, S. V. Isakov, Z. Zhu, B. O’Gorman, H. G. Katzgraber, A. Diedrich, H. Neven, J. de Kleer, B. Lackey, and R. Biswas, Physical Review Applied 12, 014004 (2019).
  * Babbush  _et al._ (2013) R. Babbush, B. O’Gorman, and A. Aspuru‐Guzik, Annalen der Physik 525, 877 (2013).
  * Xia  _et al._ (2017) R. Xia, T. Bian, and S. Kais, The Journal of Physical Chemistry B 122, 3384 (2017).
  * He  _et al._ (2013) S. He, Z. Li, and S. Zhang, Journal of the Operations Research Society of China 1, 3 (2013).
  * Babbush  _et al._ (2014) R. Babbush, A. Perdomo-Ortiz, B. O’Gorman, W. Macready, and A. Aspuru-Guzik, Advances in Chemical Physics 155, 201 (2014), arXiv:arXiv:1211.3422v2 .
  * Qian (1999) N. Qian, Neural Networks 12, 145 (1999).
  * Fu and Anderson (1986) Y. Fu and P. W. Anderson, Journal of Physics A: Mathematical and General 19, 1605 (1986).
  * Joya  _et al._ (2002) G. Joya, M. A. Atencia, and F. Sandoval, Neurocomputing 43, 219 (2002).
  * Leleu  _et al._ (2019) T. Leleu, Y. Yamamoto, P. L. McMahon, and K. Aihara, Physical Review Letters 122, 10.1103/PhysRevLett.122.040607 (2019), arXiv:1810.12565 .
  * 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).
  * Böhm  _et al._ (2021) F. Böhm, T. V. Vaerenbergh, G. Verschaffelt, and G. Van der Sande, Communications Physics 4, 149 (2021).
  * Xia  _et al._ (2018) R. Xia, T. Bian, and S. Kais, The Journal of Physical Chemistry B 122, 3384 (2018), pMID: 29099600, https://doi.org/10.1021/acs.jpcb.7b10371 .
  * Tiosue (2020) J. Tiosue, qubovert, https://github.com/jtiosue/qubovert (2020).
  * Inc. (2021) D.-W. S. Inc., dwave-neal, https://github.com/dwavesystems/dwave-neal (2021).
  * Gurobi optimisation, LLC (2022) Gurobi optimisation, LLC, Gurobi Optimizer Reference Manual (2022).
  * IBM ILOG CPLEX optimisation Studio (2022) IBM ILOG CPLEX optimisation Studio , CPLEX User’s Manual (2022).
  * McArdle  _et al._ (2020a) S. McArdle, S. Endo, A. Aspuru-Guzik, S. C. Benjamin, and X. Yuan, Rev. Mod. Phys. 92, 015003 (2020a).
  * David Sherrill and Schaefer (1999) C. David Sherrill and H. F. Schaefer (Academic Press, 1999) pp. 143–269.
  * Thachuk  _et al._ (2007) C. Thachuk, A. Shmygelska, and H. H. Hoos, BMC Bioinformatics 8, 342 (2007).
  * Guo  _et al._ (2017) Y. Guo, F. Tao, Z. Wu, and Y. Wang, BMC Systems Biology 11, 93 (2017).
  * Dill  _et al._ (2008) K. A. Dill, S. B. Ozkan, M. S. Shell, and T. R. Weikl, Annual review of biophysics 37, 289 (2008).
  * Šali  _et al._ (1994) A. Šali, E. Shakhnovich, and M. Karplus, Nature 369, 248 (1994).
  * Dill and MacCallum (2012) K. A. Dill and J. L. MacCallum, Science 338, 1042 LP (2012).
  * Englander and Mayne (2014) S. W. Englander and L. Mayne, Proceedings of the National Academy of Sciences of the United States of America 111, 15873 (2014).
  * Hansmann and Okamoto (1999) U. H. Hansmann and Y. Okamoto, Current Opinion in Structural Biology 9, 177 (1999).
  * Zhang and Chou (1992) C. Zhang and K. Chou, Biophysical Journal 63, 1523 (1992).
  * Kolinski and Skolnick (1994) A. Kolinski and J. Skolnick, Proteins: Structure, Function, and Genetics 18, 338 (1994).
  * Noé  _et al._ (2020) F. Noé, G. De Fabritiis, and C. Clementi, Current Opinion in Structural Biology 60, 77 (2020).
  * Senior  _et al._ (2020) A. W. Senior, R. Evans, J. Jumper, J. Kirkpatrick, L. Sifre, T. Green, C. Qin, A. Žídek, A. W. Nelson, A. Bridgland, H. Penedones, S. Petersen, K. Simonyan, S. Crossan, P. Kohli, D. T. Jones, D. Silver, K. Kavukcuoglu, and D. Hassabis, Nature 577, 706 (2020).
  * Babej  _et al._ (2018) T. Babej, C. Ing, and M. Fingerhuth, arXiv preprint arXiv:1811.00713 (2018).
  * Rostoker and Lotun (2005) C. Rostoker and R. Lotun, Iterative Monte Carlo Protein Design, https://www.researchgate.net/publication/228701618_Iterative_Monte_Carlo_Protein_Design (2005).
  * Szabo and Ostlund (2012) A. Szabo and N. S. Ostlund,  _Modern quantum chemistry: introduction to advanced electronic structure theory_ (Courier Corporation, 2012).
  * McArdle  _et al._ (2020b) S. McArdle, S. Endo, A. Aspuru-Guzik, S. C. Benjamin, and X. Yuan, Rev. Mod. Phys. 92, 015003 (2020b).
  * Streif  _et al._ (2019b) M. Streif, F. Neukart, and M. Leib, in  _Quantum Technology and Optimization Problems_ , edited by S. Feld and C. Linnhoff-Popien (Springer International Publishing, Cham, 2019) pp. 111–122.
  * Parrish  _et al._ (2017) R. M. Parrish, L. A. Burns, D. G. A. Smith, A. C. Simmonett, A. E. DePrince, E. G. Hohenstein, U. Bozkaya, A. Y. Sokolov, R. Di Remigio, R. M. Richard, J. F. Gonthier, A. M. James, H. R. McAlexander, A. Kumar, M. Saitow, X. Wang, B. P. Pritchard, P. Verma, H. F. Schaefer, K. Patkowski, R. A. King, E. F. Valeev, F. A. Evangelista, J. M. Turney, T. D. Crawford, and C. D. Sherrill, Journal of Chemical Theory and Computation 13, 3185 (2017), pMID: 28489372, https://doi.org/10.1021/acs.jctc.7b00174 .
  * McClean  _et al._ (2020) J. R. McClean, N. C. Rubin, K. J. Sung, I. D. Kivlichan, X. Bonet-Monroig, Y. Cao, C. Dai, E. S. Fried, C. Gidney, B. Gimby, P. Gokhale, T. Häner, T. Hardikar, V. Havlíček, O. Higgott, C. Huang, J. Izaac, Z. Jiang, X. Liu, S. McArdle, M. Neeley, T. O’Brien, B. O’Gorman, I. Ozfidan, M. D. Radin, J. Romero, N. P. D. Sawaya, B. Senjean, K. Setia, S. Sim, D. S. Steiger, M. Steudtner, Q. Sun, W. Sun, D. Wang, F. Zhang, and R. Babbush, Quantum Science and Technology 5, 034014 (2020).
  * Bravyi  _et al._ (2017) S. Bravyi, J. M. Gambetta, A. Mezzacapo, and K. Temme, arXiv preprint arXiv:1701.08213 (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, K. I. Kawarabayashi, and H. Takesue, Science Advances 7, 10.1126/sciadv.abh0952 (2021).
  * Kumar  _et al._ (2020) S. Kumar, H. Zhang, and Y. P. Huang, Communications Physics 3, 1 (2020), arXiv:2001.05680 .
  * Boros and Hammer (2002) E. Boros and P. L. Hammer, Discrete Applied Mathematics 123, 155 (2002).